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 /* Perform the reflection xnew = c + scale * (c - xold),
  46    returning 0 if xnew == c or xnew == xold (coincident points), 1 otherwise.
  47
  48    The reflected point xnew is "pinned" to the lower and upper bounds
  49    (lb and ub), as suggested by J. A. Richardson and J. L. Kuester,
  50    "The complex method for constrained optimization," Commun. ACM
  51    16(8), 487-489 (1973).  This is probably a suboptimal way to handle
  52    bound constraints, but I don't know a better way.  The main danger
  53    with this is that the simplex might collapse into a
  54    lower-dimensional hyperplane; this danger can be ameliorated by
  55    restarting (as in subplex), however. */
  56 static int reflectpt(int n, double *xnew,
  57                      const double *c, double scale, const double *xold,
  58                      const double *lb, const double *ub)
  59 {
  60      int equalc = 1, equalold = 1, i;
  61      for (i = 0; i < n; ++i) {
  62           double newx = c[i] + scale * (c[i] - xold[i]);
  63           if (newx < lb[i]) newx = lb[i];
  64           if (newx > ub[i]) newx = ub[i];
  65           equalc = equalc && (newx == c[i]);
  66           equalold = equalold && (newx == xold[i]);
  67           xnew[i] = newx;
  68      }
  69      return !(equalc || equalold);
  70 }
  71
  72 #define CHECK_EVAL(xc,fc)                                                 \
  73  stop->nevals++;                                                          \
  74  if ((fc) <= *minf) {                                                     \
  75    *minf = (fc); memcpy(x, (xc), n * sizeof(double));                     \
  76    if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; } \
  77  }                                                                        \
  78  if (nlopt_stop_evals(stop)) { ret=NLOPT_MAXEVAL_REACHED; goto done; }    \
  79  if (nlopt_stop_time(stop)) { ret=NLOPT_MAXTIME_REACHED; goto done; }
  80
  81 /* Internal version of nldrmd_minimize, intended to be used as
  82    a subroutine for the subplex method.  Three differences compared
  83    to nldrmd_minimize:
  84
  85    *minf should contain the value of f(x)  (so that we don't have to
  86    re-evaluate f at the starting x).
  87
  88    if psi > 0, then it *replaces* xtol and ftol in stop with the condition
  89    that the simplex diameter |xl - xh| must be reduced by a factor of psi
  90    ... this is for when nldrmd is used within the subplex method; for
  91    ordinary termination tests, set psi = 0.
  92
  93    scratch should contain an array of length >= (n+1)*(n+1) + 2*n,
  94    used as scratch workspace. */
  95 nlopt_result nldrmd_minimize_(int n, nlopt_func f, void *f_data,
  96                              const double *lb, const double *ub, /* bounds */
  97                              double *x, /* in: initial guess, out: minimizer */
  98                              double *minf,
  99                              const double *xstep, /* initial step sizes */
 100                              nlopt_stopping *stop,
 101                              double psi, double *scratch)
 102 {
 103      double *pts; /* (n+1) x (n+1) array of n+1 points plus function val [0] */
 104      double *c; /* centroid * n */
 105      double *xcur; /* current point */
 106      rb_tree t; /* red-black tree of simplex, sorted by f(x) */
 107      int i, j;
 108      double ninv = 1.0 / n;
 109      nlopt_result ret = NLOPT_SUCCESS;
 110      double init_diam = 0;
 111
 112      pts = scratch;
 113      c = scratch + (n+1)*(n+1);
 114      xcur = c + n;
 115
 116      rb_tree_init(&t, simplex_compare);
 117
 118      /* initialize the simplex based on the starting xstep */
 119      memcpy(pts+1, x, sizeof(double)*n);
 120      pts[0] = *minf;
 121      if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; }
 122      for (i = 0; i < n; ++i) {
 123           double *pt = pts + (i+1)*(n+1);
 124           memcpy(pt+1, x, sizeof(double)*n);
 125           if (x[i] + xstep[i] <= ub[i])
 126                pt[1+i] += xstep[i];
 127           else if (x[i] - xstep[i] >= lb[i])
 128                pt[1+i] -= xstep[i];
 129           else if (ub[i] - x[i] > x[i] - lb[i])
 130                pt[1+i] += 0.5 * (ub[i] - x[i]);
 131           else
 132                pt[1+i] -= 0.5 * (x[i] - lb[i]);
 133           if (pt[1+i] == x[i]) { ret=NLOPT_FAILURE; goto done; }
 134           pt[0] = f(n, pt+1, NULL, f_data);
 135           CHECK_EVAL(pt+1, pt[0]);
 136      }
 137
 138  restart:
 139      for (i = 0; i < n + 1; ++i)
 140           if (!rb_tree_insert(&t, pts + i*(n+1))) {
 141                ret = NLOPT_OUT_OF_MEMORY;
 142                goto done;
 143           }
 144
 145      while (1) {
 146           rb_node *low = rb_tree_min(&t);
 147           rb_node *high = rb_tree_max(&t);
 148           double fl = low->k[0], *xl = low->k + 1;
 149           double fh = high->k[0], *xh = high->k + 1;
 150           double fr;
 151
 152           if (init_diam == 0) /* initialize diam. for psi convergence test */
 153                for (i = 0; i < n; ++i) init_diam = fabs(xl[i] - xh[i]);
 154
 155           if (psi <= 0 && nlopt_stop_f(stop, fl, fh)) {
 156                ret = NLOPT_FTOL_REACHED;
 157                goto done;
 158           }
 159
 160           /* compute centroid ... if we cared about the perfomance of this,
 161              we could do it iteratively by updating the centroid on
 162              each step, but then we would have to be more careful about
 163              accumulation of rounding errors... anyway n is unlikely to
 164              be very large for Nelder-Mead in practical cases */
 165           memset(c, 0, sizeof(double)*n);
 166           for (i = 0; i < n + 1; ++i) {
 167                double *xi = pts + i*(n+1) + 1;
 168                if (xi != xh)
 169                     for (j = 0; j < n; ++j)
 170                          c[j] += xi[j];
 171           }
 172           for (i = 0; i < n; ++i) c[i] *= ninv;
 173
 174           /* x convergence check: find xcur = max radius from centroid */
 175           memset(xcur, 0, sizeof(double)*n);
 176           for (i = 0; i < n + 1; ++i) {
 177                double *xi = pts + i*(n+1) + 1;
 178                for (j = 0; j < n; ++j) {
 179                     double dx = fabs(xi[j] - c[j]);
 180                     if (dx > xcur[j]) xcur[j] = dx;
 181                }
 182           }
 183           for (i = 0; i < n; ++i) xcur[i] += c[i];
 184           if (psi > 0) {
 185                double diam = 0;
 186                for (i = 0; i < n; ++i) diam += fabs(xl[i] - xh[i]);
 187                if (diam < psi * init_diam) {
 188                     ret = NLOPT_XTOL_REACHED;
 189                     goto done;
 190                }
 191           }
 192           else if (nlopt_stop_x(stop, c, xcur)) {
 193                ret = NLOPT_XTOL_REACHED;
 194                goto done;
 195           }
 196
 197           /* reflection */
 198           if (!reflectpt(n, xcur, c, alpha, xh, lb, ub)) {
 199                ret=NLOPT_XTOL_REACHED; goto done;
 200           }
 201           fr = f(n, xcur, NULL, f_data);
 202           CHECK_EVAL(xcur, fr);
 203
 204           if (fr < fl) { /* new best point, expand simplex */
 205                if (!reflectpt(n, xh, c, gamm, xh, lb, ub)) {
 206                     ret=NLOPT_XTOL_REACHED; goto done;
 207                }
 208                fh = f(n, xh, NULL, f_data);
 209                CHECK_EVAL(xh, fh);
 210                if (fh >= fr) { /* expanding didn't improve */
 211                     fh = fr;
 212                     memcpy(xh, xcur, sizeof(double)*n);
 213                }
 214           }
 215           else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */
 216                memcpy(xh, xcur, sizeof(double)*n);
 217                fh = fr;
 218           }
 219           else { /* new worst point, contract */
 220                double fc;
 221                if (!reflectpt(n,xcur,c, fh <= fr ? -beta : beta, xh, lb,ub)) {
 222                     ret=NLOPT_XTOL_REACHED; goto done;
 223                }
 224                fc = f(n, xcur, NULL, f_data);
 225                CHECK_EVAL(xcur, fc);
 226                if (fc < fr && fc < fh) { /* successful contraction */
 227                     memcpy(xh, xcur, sizeof(double)*n);
 228                     fh = fc;
 229                }
 230                else { /* failed contraction, shrink simplex */
 231                     rb_tree_destroy(&t);
 232                     rb_tree_init(&t, simplex_compare);
 233                     for (i = 0; i < n+1; ++i) {
 234                          double *pt = pts + i * (n+1);
 235                          if (pt+1 != xl) {
 236                               if (!reflectpt(n,pt+1, xl,-delta,pt+1, lb,ub)) {
 237                                    ret = NLOPT_XTOL_REACHED;
 238                                    goto done;
 239                               }
 240                               pt[0] = f(n, pt+1, NULL, f_data);
 241                               CHECK_EVAL(pt+1, pt[0]);
 242                          }
 243                     }
 244                     goto restart;
 245                }
 246           }
 247
 248           high->k[0] = fh;
 249           rb_tree_resort(&t, high);
 250      }
 251
 252 done:
 253      rb_tree_destroy(&t);
 254      return ret;
 255 }
 256
 257 nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data,
 258                              const double *lb, const double *ub, /* bounds */
 259                              double *x, /* in: initial guess, out: minimizer */
 260                              double *minf,
 261                              const double *xstep, /* initial step sizes */
 262                              nlopt_stopping *stop)
 263 {
 264      nlopt_result ret;
 265      double *scratch;
 266
 267      *minf = f(n, x, NULL, f_data);
 268      stop->nevals++;
 269      if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 270      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 271      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 272
 273      scratch = (double*) malloc(sizeof(double) * ((n+1)*(n+1) + 2*n));
 274      if (!scratch) return NLOPT_OUT_OF_MEMORY;
 275
 276      ret = nldrmd_minimize_(n, f, f_data, lb, ub, x, minf, xstep, stop,
 277                             0.0, scratch);
 278      free(scratch);
 279      return ret;
 280 }