neldermead/nldrmd.c

   1 /* Copyright (c) 2007-2010 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 ((fc) <= *minf) {                                                     \
  82    *minf = (fc); memcpy(x, (xc), n * sizeof(double));                     \
  83    if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; } \
  84  }                                                                        \
  85  if (nlopt_stop_evals(stop)) { ret=NLOPT_MAXEVAL_REACHED; goto done; }    \
  86  if (nlopt_stop_time(stop)) { ret=NLOPT_MAXTIME_REACHED; goto done; }
  87
  88 /* Internal version of nldrmd_minimize, intended to be used as
  89    a subroutine for the subplex method.  Three differences compared
  90    to nldrmd_minimize:
  91
  92    *minf should contain the value of f(x)  (so that we don't have to
  93    re-evaluate f at the starting x).
  94
  95    if psi > 0, then it *replaces* xtol and ftol in stop with the condition
  96    that the simplex diameter |xl - xh| must be reduced by a factor of psi
  97    ... this is for when nldrmd is used within the subplex method; for
  98    ordinary termination tests, set psi = 0.
  99
 100    scratch should contain an array of length >= (n+1)*(n+1) + 2*n,
 101    used as scratch workspace.
 102
 103    On output, *fdiff will contain the difference between the high
 104    and low function values of the last simplex. */
 105 nlopt_result nldrmd_minimize_(int n, nlopt_func f, void *f_data,
 106                              const double *lb, const double *ub, /* bounds */
 107                              double *x, /* in: initial guess, out: minimizer */
 108                              double *minf,
 109                              const double *xstep, /* initial step sizes */
 110                              nlopt_stopping *stop,
 111                              double psi, double *scratch,
 112                              double *fdiff)
 113 {
 114      double *pts; /* (n+1) x (n+1) array of n+1 points plus function val [0] */
 115      double *c; /* centroid * n */
 116      double *xcur; /* current point */
 117      rb_tree t; /* red-black tree of simplex, sorted by f(x) */
 118      int i, j;
 119      double ninv = 1.0 / n;
 120      nlopt_result ret = NLOPT_SUCCESS;
 121      double init_diam = 0;
 122
 123      pts = scratch;
 124      c = scratch + (n+1)*(n+1);
 125      xcur = c + n;
 126
 127      rb_tree_init(&t, simplex_compare);
 128
 129      *fdiff = HUGE_VAL;
 130
 131      /* initialize the simplex based on the starting xstep */
 132      memcpy(pts+1, x, sizeof(double)*n);
 133      pts[0] = *minf;
 134      if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; }
 135      for (i = 0; i < n; ++i) {
 136           double *pt = pts + (i+1)*(n+1);
 137           memcpy(pt+1, x, sizeof(double)*n);
 138           pt[1+i] += xstep[i];
 139           if (pt[1+i] > ub[i]) {
 140                if (ub[i] - x[i] > fabs(xstep[i]) * 0.1)
 141                     pt[1+i] = ub[i];
 142                else /* ub is too close to pt, go in other direction */
 143                     pt[1+i] = x[i] - fabs(xstep[i]);
 144           }
 145           if (pt[1+i] < lb[i]) {
 146                if (x[i] - lb[i] > fabs(xstep[i]) * 0.1)
 147                     pt[1+i] = lb[i];
 148                else {/* lb is too close to pt, go in other direction */
 149                     pt[1+i] = x[i] + fabs(xstep[i]);
 150                     if (pt[1+i] > ub[i]) /* go towards further of lb, ub */
 151                          pt[1+i] = 0.5 * ((ub[i] - x[i] > x[i] - lb[i] ?
 152                                            ub[i] : lb[i]) + x[i]);
 153                }
 154           }
 155           if (close(pt[1+i], x[i])) { ret=NLOPT_FAILURE; goto done; }
 156           pt[0] = f(n, pt+1, NULL, f_data);
 157           CHECK_EVAL(pt+1, pt[0]);
 158      }
 159
 160  restart:
 161      for (i = 0; i < n + 1; ++i)
 162           if (!rb_tree_insert(&t, pts + i*(n+1))) {
 163                ret = NLOPT_OUT_OF_MEMORY;
 164                goto done;
 165           }
 166
 167      while (1) {
 168           rb_node *low = rb_tree_min(&t);
 169           rb_node *high = rb_tree_max(&t);
 170           double fl = low->k[0], *xl = low->k + 1;
 171           double fh = high->k[0], *xh = high->k + 1;
 172           double fr;
 173
 174           *fdiff = fh - fl;
 175
 176           if (init_diam == 0) /* initialize diam. for psi convergence test */
 177                for (i = 0; i < n; ++i) init_diam = fabs(xl[i] - xh[i]);
 178
 179           if (psi <= 0 && nlopt_stop_ftol(stop, fl, fh)) {
 180                ret = NLOPT_FTOL_REACHED;
 181                goto done;
 182           }
 183
 184           /* compute centroid ... if we cared about the perfomance of this,
 185              we could do it iteratively by updating the centroid on
 186              each step, but then we would have to be more careful about
 187              accumulation of rounding errors... anyway n is unlikely to
 188              be very large for Nelder-Mead in practical cases */
 189           memset(c, 0, sizeof(double)*n);
 190           for (i = 0; i < n + 1; ++i) {
 191                double *xi = pts + i*(n+1) + 1;
 192                if (xi != xh)
 193                     for (j = 0; j < n; ++j)
 194                          c[j] += xi[j];
 195           }
 196           for (i = 0; i < n; ++i) c[i] *= ninv;
 197
 198           /* x convergence check: find xcur = max radius from centroid */
 199           memset(xcur, 0, sizeof(double)*n);
 200           for (i = 0; i < n + 1; ++i) {
 201                double *xi = pts + i*(n+1) + 1;
 202                for (j = 0; j < n; ++j) {
 203                     double dx = fabs(xi[j] - c[j]);
 204                     if (dx > xcur[j]) xcur[j] = dx;
 205                }
 206           }
 207           for (i = 0; i < n; ++i) xcur[i] += c[i];
 208           if (psi > 0) {
 209                double diam = 0;
 210                for (i = 0; i < n; ++i) diam += fabs(xl[i] - xh[i]);
 211                if (diam < psi * init_diam) {
 212                     ret = NLOPT_XTOL_REACHED;
 213                     goto done;
 214                }
 215           }
 216           else if (nlopt_stop_x(stop, c, xcur)) {
 217                ret = NLOPT_XTOL_REACHED;
 218                goto done;
 219           }
 220
 221           /* reflection */
 222           if (!reflectpt(n, xcur, c, alpha, xh, lb, ub)) {
 223                ret=NLOPT_XTOL_REACHED; goto done;
 224           }
 225           fr = f(n, xcur, NULL, f_data);
 226           CHECK_EVAL(xcur, fr);
 227
 228           if (fr < fl) { /* new best point, expand simplex */
 229                if (!reflectpt(n, xh, c, gamm, xh, lb, ub)) {
 230                     ret=NLOPT_XTOL_REACHED; goto done;
 231                }
 232                fh = f(n, xh, NULL, f_data);
 233                CHECK_EVAL(xh, fh);
 234                if (fh >= fr) { /* expanding didn't improve */
 235                     fh = fr;
 236                     memcpy(xh, xcur, sizeof(double)*n);
 237                }
 238           }
 239           else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */
 240                memcpy(xh, xcur, sizeof(double)*n);
 241                fh = fr;
 242           }
 243           else { /* new worst point, contract */
 244                double fc;
 245                if (!reflectpt(n,xcur,c, fh <= fr ? -beta : beta, xh, lb,ub)) {
 246                     ret=NLOPT_XTOL_REACHED; goto done;
 247                }
 248                fc = f(n, xcur, NULL, f_data);
 249                CHECK_EVAL(xcur, fc);
 250                if (fc < fr && fc < fh) { /* successful contraction */
 251                     memcpy(xh, xcur, sizeof(double)*n);
 252                     fh = fc;
 253                }
 254                else { /* failed contraction, shrink simplex */
 255                     rb_tree_destroy(&t);
 256                     rb_tree_init(&t, simplex_compare);
 257                     for (i = 0; i < n+1; ++i) {
 258                          double *pt = pts + i * (n+1);
 259                          if (pt+1 != xl) {
 260                               if (!reflectpt(n,pt+1, xl,-delta,pt+1, lb,ub)) {
 261                                    ret = NLOPT_XTOL_REACHED;
 262                                    goto done;
 263                               }
 264                               pt[0] = f(n, pt+1, NULL, f_data);
 265                               CHECK_EVAL(pt+1, pt[0]);
 266                          }
 267                     }
 268                     goto restart;
 269                }
 270           }
 271
 272           high->k[0] = fh;
 273           rb_tree_resort(&t, high);
 274      }
 275
 276 done:
 277      rb_tree_destroy(&t);
 278      return ret;
 279 }
 280
 281 nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data,
 282                              const double *lb, const double *ub, /* bounds */
 283                              double *x, /* in: initial guess, out: minimizer */
 284                              double *minf,
 285                              const double *xstep, /* initial step sizes */
 286                              nlopt_stopping *stop)
 287 {
 288      nlopt_result ret;
 289      double *scratch, fdiff;
 290
 291      *minf = f(n, x, NULL, f_data);
 292      stop->nevals++;
 293      if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 294      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 295      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 296
 297      scratch = (double*) malloc(sizeof(double) * ((n+1)*(n+1) + 2*n));
 298      if (!scratch) return NLOPT_OUT_OF_MEMORY;
 299
 300      ret = nldrmd_minimize_(n, f, f_data, lb, ub, x, minf, xstep, stop,
 301                             0.0, scratch, &fdiff);
 302      free(scratch);
 303      return ret;
 304 }