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