neldermead/nldrmd.c

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