neldermead/sbplx.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 <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 (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
  83      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
  84      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
  85
  86      xstep = (double*)malloc(sizeof(double) * (n*3 + nsmax*4
  87                                                + (nsmax+1)*(nsmax+1)+2*nsmax));
  88      if (!xstep) return NLOPT_OUT_OF_MEMORY;
  89      xprev = xstep + n; dx = xprev + n;
  90      xs = dx + n; xsstep = xs + nsmax;
  91      lbs = xsstep + nsmax; ubs = lbs + nsmax;
  92      scratch = ubs + nsmax;
  93      p = (int *) malloc(sizeof(int) * n);
  94      if (!p) { free(xstep); return NLOPT_OUT_OF_MEMORY; }
  95
  96      memcpy(xstep, xstep0, n * sizeof(double));
  97      memset(dx, 0, n * sizeof(double));
  98
  99      sd.p = p;
 100      sd.n = n;
 101      sd.x = x;
 102      sd.f = f;
 103      sd.f_data = f_data;
 104
 105      while (1) {
 106           double normi = 0;
 107           double normdx = 0;
 108           int ns, nsubs = 0;
 109           int nevals = stop->nevals;
 110
 111           memcpy(xprev, x, n * sizeof(double));
 112           fprev = *minf;
 113
 114           /* sort indices into the progress vector dx by decreasing
 115              order of magnitude |dx| */
 116           for (i = 0; i < n; ++i) p[i] = i;
 117           nlopt_qsort_r(p, (size_t) n, sizeof(int), dx, p_compare);
 118
 119           /* find the subspaces, and perform nelder-mead on each one */
 120           for (i = 0; i < n; ++i) normdx += fabs(dx[i]); /* L1 norm */
 121           for (i = 0; i + nsmin < n; i += ns) {
 122                /* find subspace starting at index i */
 123                int k, nk;
 124                double ns_goodness = -HUGE_VAL, norm = normi;
 125                nk = i+nsmax > n ? n : i+nsmax; /* max k for this subspace */
 126                for (k = i; k < i+nsmin-1; ++k) norm += fabs(dx[p[k]]);
 127                ns = nsmin;
 128                for (k = i+nsmin-1; k < nk; ++k) {
 129                     double goodness;
 130                     norm += fabs(dx[p[k]]);
 131                     /* remaining subspaces must be big enough to partition */
 132                     if (n-(k+1) < nsmin) continue;
 133                     /* maximize figure of merit defined by Rowan thesis:
 134                        look for sudden drops in average |dx| */
 135                     if (k+1 < n)
 136                          goodness = norm/(k+1) - (normdx-norm)/(n-(k+1));
 137                     else
 138                          goodness = normdx/n;
 139                     if (goodness > ns_goodness) {
 140                          ns_goodness = goodness;
 141                          ns = (k+1)-i;
 142                     }
 143                }
 144                for (k = i; k < i+ns; ++k) normi += fabs(dx[p[k]]);
 145                /* do nelder-mead on subspace of dimension ns starting w/i */
 146                sd.is = i;
 147                for (k = i; k < i+ns; ++k) {
 148                     xs[k-i] = x[p[k]];
 149                     xsstep[k-i] = xstep[p[k]];
 150                     lbs[k-i] = lb[p[k]];
 151                     ubs[k-i] = ub[p[k]];
 152                }
 153                ++nsubs;
 154                nevals = stop->nevals;
 155                ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
 156                                       xsstep, stop, psi, scratch);
 157                if (sbplx_verbose)
 158                     printf("%d NM iterations for (%d,%d) subspace\n",
 159                            stop->nevals - nevals, sd.is, ns);
 160                for (k = i; k < i+ns; ++k) x[p[k]] = xs[k-i];
 161                if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
 162                if (ret != NLOPT_XTOL_REACHED) goto done;
 163           }
 164           /* nelder-mead on last subspace */
 165           ns = n - i;
 166           sd.is = i;
 167           for (; i < n; ++i) {
 168                xs[i-sd.is] = x[p[i]];
 169                xsstep[i-sd.is] = xstep[p[i]];
 170                lbs[i-sd.is] = lb[p[i]];
 171                ubs[i-sd.is] = ub[p[i]];
 172           }
 173           ++nsubs;
 174           nevals = stop->nevals;
 175           ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
 176                                  xsstep, stop, psi, scratch);
 177           if (sbplx_verbose)
 178                printf("%d NM iterations for (%d,%d) subspace\n",
 179                       stop->nevals - nevals, sd.is, ns);
 180           for (i = sd.is; i < n; ++i) x[p[i]] = xs[i-sd.is];
 181           if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
 182           if (ret != NLOPT_XTOL_REACHED) goto done;
 183
 184           /* termination tests: */
 185           if (nlopt_stop_ftol(stop, *minf, fprev)) {
 186                ret = NLOPT_FTOL_REACHED;
 187                goto done;
 188           }
 189           if (nlopt_stop_x(stop, x, xprev)) {
 190                int j;
 191                /* as explained in Rowan's thesis, it is important
 192                   to check |xstep| as well as |x-xprev|, since if
 193                   the step size is too large (in early iterations),
 194                   the inner Nelder-Mead may not make much progress */
 195                for (j = 0; j < n; ++j)
 196                     if (fabs(xstep[j]) * psi > stop->xtol_abs[j]
 197                         && fabs(xstep[j]) * psi > stop->xtol_rel * fabs(x[j]))
 198                          break;
 199                if (j == n) {
 200                     ret = NLOPT_XTOL_REACHED;
 201                     goto done;
 202                }
 203           }
 204
 205           /* compute change in optimal point */
 206           for (i = 0; i < n; ++i) dx[i] = x[i] - xprev[i];
 207
 208           /* setting stepsizes */
 209           {
 210                double scale;
 211                if (nsubs == 1)
 212                     scale = psi;
 213                else {
 214                     double stepnorm = 0, dxnorm = 0;
 215                     for (i = 0; i < n; ++i) {
 216                          stepnorm += fabs(xstep[i]);
 217                          dxnorm += fabs(dx[i]);
 218                     }
 219                     scale = dxnorm / stepnorm;
 220                     if (scale < omega) scale = omega;
 221                     if (scale > 1/omega) scale = 1/omega;
 222                }
 223                for (i = 0; i < n; ++i)
 224                     xstep[i] = (dx[i] == 0) ? -(xstep[i] * scale)
 225                          : copysign(xstep[i] * scale, dx[i]);
 226           }
 227      }
 228
 229  done:
 230      free(p);
 231      free(xstep);
 232      return ret;
 233 }