src/algs/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_p);                                                    \
  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])) {
 157               nlopt_stop_msg(stop, "starting step size led to simplex that was too small in dimension %d: %g is too close to x[%d]=%g",
 158                              i, pt[1+i], i, x[i]);
 159               ret=NLOPT_FAILURE;
 160               goto done;
 161           }
 162           pt[0] = f(n, pt+1, NULL, f_data);
 163           CHECK_EVAL(pt+1, pt[0]);
 164      }
 165
 166  restart:
 167      for (i = 0; i < n + 1; ++i)
 168           if (!rb_tree_insert(&t, pts + i*(n+1))) {
 169                ret = NLOPT_OUT_OF_MEMORY;
 170                goto done;
 171           }
 172
 173      while (1) {
 174           rb_node *low = rb_tree_min(&t);
 175           rb_node *high = rb_tree_max(&t);
 176           double fl = low->k[0], *xl = low->k + 1;
 177           double fh = high->k[0], *xh = high->k + 1;
 178           double fr;
 179
 180           *fdiff = fh - fl;
 181
 182           if (init_diam == 0) /* initialize diam. for psi convergence test */
 183                for (i = 0; i < n; ++i) init_diam += fabs(xl[i] - xh[i]);
 184
 185           if (psi <= 0 && nlopt_stop_ftol(stop, fl, fh)) {
 186                ret = NLOPT_FTOL_REACHED;
 187                goto done;
 188           }
 189
 190           /* compute centroid ... if we cared about the perfomance of this,
 191              we could do it iteratively by updating the centroid on
 192              each step, but then we would have to be more careful about
 193              accumulation of rounding errors... anyway n is unlikely to
 194              be very large for Nelder-Mead in practical cases */
 195           memset(c, 0, sizeof(double)*n);
 196           for (i = 0; i < n + 1; ++i) {
 197                double *xi = pts + i*(n+1) + 1;
 198                if (xi != xh)
 199                     for (j = 0; j < n; ++j)
 200                          c[j] += xi[j];
 201           }
 202           for (i = 0; i < n; ++i) c[i] *= ninv;
 203
 204           /* x convergence check: find xcur = max radius from centroid */
 205           memset(xcur, 0, sizeof(double)*n);
 206           for (i = 0; i < n + 1; ++i) {
 207                double *xi = pts + i*(n+1) + 1;
 208                for (j = 0; j < n; ++j) {
 209                     double dx = fabs(xi[j] - c[j]);
 210                     if (dx > xcur[j]) xcur[j] = dx;
 211                }
 212           }
 213           for (i = 0; i < n; ++i) xcur[i] += c[i];
 214           if (psi > 0) {
 215                double diam = 0;
 216                for (i = 0; i < n; ++i) diam += fabs(xl[i] - xh[i]);
 217                if (diam < psi * init_diam) {
 218                     ret = NLOPT_XTOL_REACHED;
 219                     goto done;
 220                }
 221           }
 222           else if (nlopt_stop_x(stop, c, xcur)) {
 223                ret = NLOPT_XTOL_REACHED;
 224                goto done;
 225           }
 226
 227           /* reflection */
 228           if (!reflectpt(n, xcur, c, alpha, xh, lb, ub)) {
 229                ret=NLOPT_XTOL_REACHED; goto done;
 230           }
 231           fr = f(n, xcur, NULL, f_data);
 232           CHECK_EVAL(xcur, fr);
 233
 234           if (fr < fl) { /* new best point, expand simplex */
 235                if (!reflectpt(n, xh, c, gamm, xh, lb, ub)) {
 236                     ret=NLOPT_XTOL_REACHED; goto done;
 237                }
 238                fh = f(n, xh, NULL, f_data);
 239                CHECK_EVAL(xh, fh);
 240                if (fh >= fr) { /* expanding didn't improve */
 241                     fh = fr;
 242                     memcpy(xh, xcur, sizeof(double)*n);
 243                }
 244           }
 245           else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */
 246                memcpy(xh, xcur, sizeof(double)*n);
 247                fh = fr;
 248           }
 249           else { /* new worst point, contract */
 250                double fc;
 251                if (!reflectpt(n,xcur,c, fh <= fr ? -beta : beta, xh, lb,ub)) {
 252                     ret=NLOPT_XTOL_REACHED; goto done;
 253                }
 254                fc = f(n, xcur, NULL, f_data);
 255                CHECK_EVAL(xcur, fc);
 256                if (fc < fr && fc < fh) { /* successful contraction */
 257                     memcpy(xh, xcur, sizeof(double)*n);
 258                     fh = fc;
 259                }
 260                else { /* failed contraction, shrink simplex */
 261                     rb_tree_destroy(&t);
 262                     rb_tree_init(&t, simplex_compare);
 263                     for (i = 0; i < n+1; ++i) {
 264                          double *pt = pts + i * (n+1);
 265                          if (pt+1 != xl) {
 266                               if (!reflectpt(n,pt+1, xl,-delta,pt+1, lb,ub)) {
 267                                    ret = NLOPT_XTOL_REACHED;
 268                                    goto done;
 269                               }
 270                               pt[0] = f(n, pt+1, NULL, f_data);
 271                               CHECK_EVAL(pt+1, pt[0]);
 272                          }
 273                     }
 274                     goto restart;
 275                }
 276           }
 277
 278           high->k[0] = fh;
 279           rb_tree_resort(&t, high);
 280      }
 281
 282 done:
 283      rb_tree_destroy(&t);
 284      return ret;
 285 }
 286
 287 nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data,
 288                              const double *lb, const double *ub, /* bounds */
 289                              double *x, /* in: initial guess, out: minimizer */
 290                              double *minf,
 291                              const double *xstep, /* initial step sizes */
 292                              nlopt_stopping *stop)
 293 {
 294      nlopt_result ret;
 295      double *scratch, fdiff;
 296
 297      *minf = f(n, x, NULL, f_data);
 298      ++ *(stop->nevals_p);
 299      if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
 300      if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 301      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 302      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 303
 304      scratch = (double*) malloc(sizeof(double) * ((n+1)*(n+1) + 2*n));
 305      if (!scratch) return NLOPT_OUT_OF_MEMORY;
 306
 307      ret = nldrmd_minimize_(n, f, f_data, lb, ub, x, minf, xstep, stop,
 308                             0.0, scratch, &fdiff);
 309      free(scratch);
 310      return ret;
 311 }