neldermead/sbplx.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 <stdio.h>
  26 #include <string.h>
  27
  28 #include "neldermead.h"
  29
  30 /* subplex strategy constants: */
  31 static const double psi = 0.25, omega = 0.1;
  32 static const int nsmin = 2, nsmax = 5;
  33
  34 int sbplx_verbose = 0; /* for debugging */
  35
  36 /* qsort_r comparison function for sorting indices into delta-x array */
  37 static int p_compare(void *dx_, const void *i_, const void *j_)
  38 {
  39      const double *dx = (const double *) dx_;
  40      int i = *((const int *) i_), j = *((const int *) j_);
  41      double dxi = fabs(dx[i]), dxj = fabs(dx[j]);
  42      return (dxi > dxj ? -1 : (dxi < dxj ? +1 : 0));
  43 }
  44
  45 typedef struct {
  46      const int *p; /* subspace index permutation */
  47      int is; /* starting index for this subspace */
  48      int n; /* dimension of underlying space */
  49      double *x; /* current x vector */
  50      nlopt_func f; void *f_data; /* the "actual" underlying function */
  51 } subspace_data;
  52
  53 /* wrapper around objective function for subspace optimization */
  54 static double subspace_func(int ns, const double *xs, double *grad, void *data)
  55 {
  56      subspace_data *d = (subspace_data *) data;
  57      int i, is = d->is;
  58      const int *p = d->p;
  59      double *x = d->x;
  60
  61      (void) grad; /* should always be NULL here */
  62      for (i = is; i < is + ns; ++i) x[p[i]] = xs[i-is];
  63      return d->f(d->n, x, NULL, d->f_data);
  64 }
  65
  66 nlopt_result sbplx_minimize(int n, nlopt_func f, void *f_data,
  67                             const double *lb, const double *ub, /* bounds */
  68                             double *x, /* in: initial guess, out: minimizer */
  69                             double *minf,
  70                             const double *xstep0, /* initial step sizes */
  71                             nlopt_stopping *stop)
  72 {
  73      nlopt_result ret = NLOPT_SUCCESS;
  74      double *xstep, *xprev, *dx, *xs, *lbs, *ubs, *xsstep, *scratch;
  75      int *p; /* permuted indices of x sorted by decreasing magnitude |dx| */
  76      int i;
  77      subspace_data sd;
  78      double fprev;
  79
  80      *minf = f(n, x, NULL, f_data);
  81      stop->nevals++;
  82      if (nlopt_stop_forced(stop)) return NLOPT_FORCE_STOP;
  83      if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
  84      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
  85      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
  86
  87      xstep = (double*)malloc(sizeof(double) * (n*3 + nsmax*4
  88                                                + (nsmax+1)*(nsmax+1)+2*nsmax));
  89      if (!xstep) return NLOPT_OUT_OF_MEMORY;
  90      xprev = xstep + n; dx = xprev + n;
  91      xs = dx + n; xsstep = xs + nsmax;
  92      lbs = xsstep + nsmax; ubs = lbs + nsmax;
  93      scratch = ubs + nsmax;
  94      p = (int *) malloc(sizeof(int) * n);
  95      if (!p) { free(xstep); return NLOPT_OUT_OF_MEMORY; }
  96
  97      memcpy(xstep, xstep0, n * sizeof(double));
  98      memset(dx, 0, n * sizeof(double));
  99
 100      sd.p = p;
 101      sd.n = n;
 102      sd.x = x;
 103      sd.f = f;
 104      sd.f_data = f_data;
 105
 106      while (1) {
 107           double normi = 0;
 108           double normdx = 0;
 109           int ns, nsubs = 0;
 110           int nevals = stop->nevals;
 111           double fdiff, fdiff_max = 0;
 112
 113           memcpy(xprev, x, n * sizeof(double));
 114           fprev = *minf;
 115
 116           /* sort indices into the progress vector dx by decreasing
 117              order of magnitude |dx| */
 118           for (i = 0; i < n; ++i) p[i] = i;
 119           nlopt_qsort_r(p, (size_t) n, sizeof(int), dx, p_compare);
 120
 121           /* find the subspaces, and perform nelder-mead on each one */
 122           for (i = 0; i < n; ++i) normdx += fabs(dx[i]); /* L1 norm */
 123           for (i = 0; i + nsmin < n; i += ns) {
 124                /* find subspace starting at index i */
 125                int k, nk;
 126                double ns_goodness = -HUGE_VAL, norm = normi;
 127                nk = i+nsmax > n ? n : i+nsmax; /* max k for this subspace */
 128                for (k = i; k < i+nsmin-1; ++k) norm += fabs(dx[p[k]]);
 129                ns = nsmin;
 130                for (k = i+nsmin-1; k < nk; ++k) {
 131                     double goodness;
 132                     norm += fabs(dx[p[k]]);
 133                     /* remaining subspaces must be big enough to partition */
 134                     if (n-(k+1) < nsmin) continue;
 135                     /* maximize figure of merit defined by Rowan thesis:
 136                        look for sudden drops in average |dx| */
 137                     if (k+1 < n)
 138                          goodness = norm/(k+1) - (normdx-norm)/(n-(k+1));
 139                     else
 140                          goodness = normdx/n;
 141                     if (goodness > ns_goodness) {
 142                          ns_goodness = goodness;
 143                          ns = (k+1)-i;
 144                     }
 145                }
 146                for (k = i; k < i+ns; ++k) normi += fabs(dx[p[k]]);
 147                /* do nelder-mead on subspace of dimension ns starting w/i */
 148                sd.is = i;
 149                for (k = i; k < i+ns; ++k) {
 150                     xs[k-i] = x[p[k]];
 151                     xsstep[k-i] = xstep[p[k]];
 152                     lbs[k-i] = lb[p[k]];
 153                     ubs[k-i] = ub[p[k]];
 154                }
 155                ++nsubs;
 156                nevals = stop->nevals;
 157                ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
 158                                       xsstep, stop, psi, scratch, &fdiff);
 159                if (fdiff > fdiff_max) fdiff_max = fdiff;
 160                if (sbplx_verbose)
 161                     printf("%d NM iterations for (%d,%d) subspace\n",
 162                            stop->nevals - nevals, sd.is, ns);
 163                for (k = i; k < i+ns; ++k) x[p[k]] = xs[k-i];
 164                if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
 165                if (ret != NLOPT_XTOL_REACHED) goto done;
 166           }
 167           /* nelder-mead on last subspace */
 168           ns = n - i;
 169           sd.is = i;
 170           for (; i < n; ++i) {
 171                xs[i-sd.is] = x[p[i]];
 172                xsstep[i-sd.is] = xstep[p[i]];
 173                lbs[i-sd.is] = lb[p[i]];
 174                ubs[i-sd.is] = ub[p[i]];
 175           }
 176           ++nsubs;
 177           nevals = stop->nevals;
 178           ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
 179                                  xsstep, stop, psi, scratch, &fdiff);
 180           if (fdiff > fdiff_max) fdiff_max = fdiff;
 181           if (sbplx_verbose)
 182                printf("sbplx: %d NM iterations for (%d,%d) subspace\n",
 183                       stop->nevals - nevals, sd.is, ns);
 184           for (i = sd.is; i < n; ++i) x[p[i]] = xs[i-sd.is];
 185           if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
 186           if (ret != NLOPT_XTOL_REACHED) goto done;
 187
 188           /* termination tests: */
 189           if (nlopt_stop_ftol(stop, *minf, *minf + fdiff_max)) {
 190                ret = NLOPT_FTOL_REACHED;
 191                goto done;
 192           }
 193           if (nlopt_stop_x(stop, x, xprev)) {
 194                int j;
 195                /* as explained in Rowan's thesis, it is important
 196                   to check |xstep| as well as |x-xprev|, since if
 197                   the step size is too large (in early iterations),
 198                   the inner Nelder-Mead may not make much progress */
 199                for (j = 0; j < n; ++j)
 200                     if (fabs(xstep[j]) * psi > stop->xtol_abs[j]
 201                         && fabs(xstep[j]) * psi > stop->xtol_rel * fabs(x[j]))
 202                          break;
 203                if (j == n) {
 204                     ret = NLOPT_XTOL_REACHED;
 205                     goto done;
 206                }
 207           }
 208
 209           /* compute change in optimal point */
 210           for (i = 0; i < n; ++i) dx[i] = x[i] - xprev[i];
 211
 212           /* setting stepsizes */
 213           {
 214                double scale;
 215                if (nsubs == 1)
 216                     scale = psi;
 217                else {
 218                     double stepnorm = 0, dxnorm = 0;
 219                     for (i = 0; i < n; ++i) {
 220                          stepnorm += fabs(xstep[i]);
 221                          dxnorm += fabs(dx[i]);
 222                     }
 223                     scale = dxnorm / stepnorm;
 224                     if (scale < omega) scale = omega;
 225                     if (scale > 1/omega) scale = 1/omega;
 226                }
 227                if (sbplx_verbose)
 228                     printf("sbplx: stepsize scale factor = %g\n", scale);
 229                for (i = 0; i < n; ++i)
 230                     xstep[i] = (dx[i] == 0) ? -(xstep[i] * scale)
 231                          : copysign(xstep[i] * scale, dx[i]);
 232           }
 233      }
 234
 235  done:
 236      free(p);
 237      free(xstep);
 238      return ret;
 239 }