neldermead/nldrmd.c

   1 /* Copyright (c) 2007-2014 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <math.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26
  27 #include "neldermead.h"
  28 #include "redblack.h"
  29
  30 /* Nelder-Mead simplex algorithm, used as a subroutine for the Rowan's
  31    subplex algorithm.  Modified to handle bound constraints ala
  32    Richardson and Kuester (1973), as mentioned below. */
  33
  34 /* heuristic "strategy" constants: */
  35 static const double alpha = 1, beta = 0.5, gamm = 2, delta = 0.5;
  36
  37 /* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
  38 static int simplex_compare(double *k1, double *k2)
  39 {
  40      if (*k1 < *k2) return -1;
  41      if (*k1 > *k2) return +1;
  42      return k1 - k2; /* tie-breaker */
  43 }
  44
  45 /* return 1 if a and b are approximately equal relative to floating-point
  46    precision, 0 otherwise */
  47 static int close(double a, double b)
  48 {
  49      return (fabs(a - b) <= 1e-13 * (fabs(a) + fabs(b)));
  50 }
  51
  52 /* Perform the reflection xnew = c + scale * (c - xold),
  53    returning 0 if xnew == c or xnew == xold (coincident points), 1 otherwise.
  54
  55    The reflected point xnew is "pinned" to the lower and upper bounds
  56    (lb and ub), as suggested by J. A. Richardson and J. L. Kuester,
  57    "The complex method for constrained optimization," Commun. ACM
  58    16(8), 487-489 (1973).  This is probably a suboptimal way to handle
  59    bound constraints, but I don't know a better way.  The main danger
  60    with this is that the simplex might collapse into a
  61    lower-dimensional hyperplane; this danger can be ameliorated by
  62    restarting (as in subplex), however. */
  63 static int reflectpt(int n, double *xnew,
  64                      const double *c, double scale, const double *xold,
  65                      const double *lb, const double *ub)
  66 {
  67      int equalc = 1, equalold = 1, i;
  68      for (i = 0; i < n; ++i) {
  69           double newx = c[i] + scale * (c[i] - xold[i]);
  70           if (newx < lb[i]) newx = lb[i];
  71           if (newx > ub[i]) newx = ub[i];
  72           equalc = equalc && close(newx, c[i]);
  73           equalold = equalold && close(newx, xold[i]);
  74           xnew[i] = newx;
  75      }
  76      return !(equalc || equalold);
  77 }
  78
  79 #define CHECK_EVAL(xc,fc)                                                 \
  80  stop->nevals++;                                                          \
  81  if (nlopt_stop_forced(stop)) { ret=NLOPT_FORCED_STOP; goto done; }        \
  82  if ((fc) <= *minf) {                                                     \
  83    *minf = (fc); memcpy(x, (xc), n * sizeof(double));                     \
  84    if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; } \
  85  }                                                                        \
  86  if (nlopt_stop_evals(stop)) { ret=NLOPT_MAXEVAL_REACHED; goto done; }    \
  87  if (nlopt_stop_time(stop)) { ret=NLOPT_MAXTIME_REACHED; goto done; }
  88
  89 /* Internal version of nldrmd_minimize, intended to be used as
  90    a subroutine for the subplex method.  Three differences compared
  91    to nldrmd_minimize:
  92
  93    *minf should contain the value of f(x)  (so that we don't have to
  94    re-evaluate f at the starting x).
  95
  96    if psi > 0, then it *replaces* xtol and ftol in stop with the condition
  97    that the simplex diameter |xl - xh| must be reduced by a factor of psi
  98    ... this is for when nldrmd is used within the subplex method; for
  99    ordinary termination tests, set psi = 0.
 100
 101    scratch should contain an array of length >= (n+1)*(n+1) + 2*n,
 102    used as scratch workspace.
 103
 104    On output, *fdiff will contain the difference between the high
 105    and low function values of the last simplex. */
 106 nlopt_result nldrmd_minimize_(int n, nlopt_func f, void *f_data,
 107                              const double *lb, const double *ub, /* bounds */
 108                              double *x, /* in: initial guess, out: minimizer */
 109                              double *minf,
 110                              const double *xstep, /* initial step sizes */
 111                              nlopt_stopping *stop,
 112                              double psi, double *scratch,
 113                              double *fdiff)
 114 {
 115      double *pts; /* (n+1) x (n+1) array of n+1 points plus function val [0] */
 116      double *c; /* centroid * n */
 117      double *xcur; /* current point */
 118      rb_tree t; /* red-black tree of simplex, sorted by f(x) */
 119      int i, j;
 120      double ninv = 1.0 / n;
 121      nlopt_result ret = NLOPT_SUCCESS;
 122      double init_diam = 0;
 123
 124      pts = scratch;
 125      c = scratch + (n+1)*(n+1);
 126      xcur = c + n;
 127
 128      rb_tree_init(&t, simplex_compare);
 129
 130      *fdiff = HUGE_VAL;
 131
 132      /* initialize the simplex based on the starting xstep */
 133      memcpy(pts+1, x, sizeof(double)*n);
 134      pts[0] = *minf;
 135      if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; }
 136      for (i = 0; i < n; ++i) {
 137           double *pt = pts + (i+1)*(n+1);
 138           memcpy(pt+1, x, sizeof(double)*n);
 139           pt[1+i] += xstep[i];
 140           if (pt[1+i] > ub[i]) {
 141                if (ub[i] - x[i] > fabs(xstep[i]) * 0.1)
 142                     pt[1+i] = ub[i];
 143                else /* ub is too close to pt, go in other direction */
 144                     pt[1+i] = x[i] - fabs(xstep[i]);
 145           }
 146           if (pt[1+i] < lb[i]) {
 147                if (x[i] - lb[i] > fabs(xstep[i]) * 0.1)
 148                     pt[1+i] = lb[i];
 149                else {/* lb is too close to pt, go in other direction */
 150                     pt[1+i] = x[i] + fabs(xstep[i]);
 151                     if (pt[1+i] > ub[i]) /* go towards further of lb, ub */
 152                          pt[1+i] = 0.5 * ((ub[i] - x[i] > x[i] - lb[i] ?
 153                                            ub[i] : lb[i]) + x[i]);
 154                }
 155           }
 156           if (close(pt[1+i], x[i])) { ret=NLOPT_FAILURE; goto done; }
 157           pt[0] = f(n, pt+1, NULL, f_data);
 158           CHECK_EVAL(pt+1, pt[0]);
 159      }
 160
 161  restart:
 162      for (i = 0; i < n + 1; ++i)
 163           if (!rb_tree_insert(&t, pts + i*(n+1))) {
 164                ret = NLOPT_OUT_OF_MEMORY;
 165                goto done;
 166           }
 167
 168      while (1) {
 169           rb_node *low = rb_tree_min(&t);
 170           rb_node *high = rb_tree_max(&t);
 171           double fl = low->k[0], *xl = low->k + 1;
 172           double fh = high->k[0], *xh = high->k + 1;
 173           double fr;
 174
 175           *fdiff = fh - fl;
 176
 177           if (init_diam == 0) /* initialize diam. for psi convergence test */
 178                for (i = 0; i < n; ++i) init_diam += fabs(xl[i] - xh[i]);
 179
 180           if (psi <= 0 && nlopt_stop_ftol(stop, fl, fh)) {
 181                ret = NLOPT_FTOL_REACHED;
 182                goto done;
 183           }
 184
 185           /* compute centroid ... if we cared about the perfomance of this,
 186              we could do it iteratively by updating the centroid on
 187              each step, but then we would have to be more careful about
 188              accumulation of rounding errors... anyway n is unlikely to
 189              be very large for Nelder-Mead in practical cases */
 190           memset(c, 0, sizeof(double)*n);
 191           for (i = 0; i < n + 1; ++i) {
 192                double *xi = pts + i*(n+1) + 1;
 193                if (xi != xh)
 194                     for (j = 0; j < n; ++j)
 195                          c[j] += xi[j];
 196           }
 197           for (i = 0; i < n; ++i) c[i] *= ninv;
 198
 199           /* x convergence check: find xcur = max radius from centroid */
 200           memset(xcur, 0, sizeof(double)*n);
 201           for (i = 0; i < n + 1; ++i) {
 202                double *xi = pts + i*(n+1) + 1;
 203                for (j = 0; j < n; ++j) {
 204                     double dx = fabs(xi[j] - c[j]);
 205                     if (dx > xcur[j]) xcur[j] = dx;
 206                }
 207           }
 208           for (i = 0; i < n; ++i) xcur[i] += c[i];
 209           if (psi > 0) {
 210                double diam = 0;
 211                for (i = 0; i < n; ++i) diam += fabs(xl[i] - xh[i]);
 212                if (diam < psi * init_diam) {
 213                     ret = NLOPT_XTOL_REACHED;
 214                     goto done;
 215                }
 216           }
 217           else if (nlopt_stop_x(stop, c, xcur)) {
 218                ret = NLOPT_XTOL_REACHED;
 219                goto done;
 220           }
 221
 222           /* reflection */
 223           if (!reflectpt(n, xcur, c, alpha, xh, lb, ub)) {
 224                ret=NLOPT_XTOL_REACHED; goto done;
 225           }
 226           fr = f(n, xcur, NULL, f_data);
 227           CHECK_EVAL(xcur, fr);
 228
 229           if (fr < fl) { /* new best point, expand simplex */
 230                if (!reflectpt(n, xh, c, gamm, xh, lb, ub)) {
 231                     ret=NLOPT_XTOL_REACHED; goto done;
 232                }
 233                fh = f(n, xh, NULL, f_data);
 234                CHECK_EVAL(xh, fh);
 235                if (fh >= fr) { /* expanding didn't improve */
 236                     fh = fr;
 237                     memcpy(xh, xcur, sizeof(double)*n);
 238                }
 239           }
 240           else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */
 241                memcpy(xh, xcur, sizeof(double)*n);
 242                fh = fr;
 243           }
 244           else { /* new worst point, contract */
 245                double fc;
 246                if (!reflectpt(n,xcur,c, fh <= fr ? -beta : beta, xh, lb,ub)) {
 247                     ret=NLOPT_XTOL_REACHED; goto done;
 248                }
 249                fc = f(n, xcur, NULL, f_data);
 250                CHECK_EVAL(xcur, fc);
 251                if (fc < fr && fc < fh) { /* successful contraction */
 252                     memcpy(xh, xcur, sizeof(double)*n);
 253                     fh = fc;
 254                }
 255                else { /* failed contraction, shrink simplex */
 256                     rb_tree_destroy(&t);
 257                     rb_tree_init(&t, simplex_compare);
 258                     for (i = 0; i < n+1; ++i) {
 259                          double *pt = pts + i * (n+1);
 260                          if (pt+1 != xl) {
 261                               if (!reflectpt(n,pt+1, xl,-delta,pt+1, lb,ub)) {
 262                                    ret = NLOPT_XTOL_REACHED;
 263                                    goto done;
 264                               }
 265                               pt[0] = f(n, pt+1, NULL, f_data);
 266                               CHECK_EVAL(pt+1, pt[0]);
 267                          }
 268                     }
 269                     goto restart;
 270                }
 271           }
 272
 273           high->k[0] = fh;
 274           rb_tree_resort(&t, high);
 275      }
 276
 277 done:
 278      rb_tree_destroy(&t);
 279      return ret;
 280 }
 281
 282 nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data,
 283                              const double *lb, const double *ub, /* bounds */
 284                              double *x, /* in: initial guess, out: minimizer */
 285                              double *minf,
 286                              const double *xstep, /* initial step sizes */
 287                              nlopt_stopping *stop)
 288 {
 289      nlopt_result ret;
 290      double *scratch, fdiff;
 291
 292      *minf = f(n, x, NULL, f_data);
 293      stop->nevals++;
 294      if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
 295      if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 296      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 297      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 298
 299      scratch = (double*) malloc(sizeof(double) * ((n+1)*(n+1) + 2*n));
 300      if (!scratch) return NLOPT_OUT_OF_MEMORY;
 301
 302      ret = nldrmd_minimize_(n, f, f_data, lb, ub, x, minf, xstep, stop,
 303                             0.0, scratch, &fdiff);
 304      free(scratch);
 305      return ret;
 306 }