1 /* Copyright (c) 2007-2010 Massachusetts Institute of Technology
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:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
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.
28 #include "neldermead.h"
30 /* subplex strategy constants: */
31 static const double psi = 0.25, omega = 0.1;
32 static const int nsmin = 2, nsmax = 5;
34 int sbplx_verbose = 0; /* for debugging */
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_)
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));
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 */
53 /* wrapper around objective function for subspace optimization */
54 static double subspace_func(int ns, const double *xs, double *grad, void *data)
56 subspace_data *d = (subspace_data *) data;
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);
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 */
70 const double *xstep0, /* initial step sizes */
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| */
80 *minf = f(n, x, NULL, f_data);
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;
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; }
96 memcpy(xstep, xstep0, n * sizeof(double));
97 memset(dx, 0, n * sizeof(double));
109 int nevals = stop->nevals;
110 double fdiff, fdiff_max = 0;
112 memcpy(xprev, x, n * sizeof(double));
115 /* sort indices into the progress vector dx by decreasing
116 order of magnitude |dx| */
117 for (i = 0; i < n; ++i) p[i] = i;
118 nlopt_qsort_r(p, (size_t) n, sizeof(int), dx, p_compare);
120 /* find the subspaces, and perform nelder-mead on each one */
121 for (i = 0; i < n; ++i) normdx += fabs(dx[i]); /* L1 norm */
122 for (i = 0; i + nsmin < n; i += ns) {
123 /* find subspace starting at index i */
125 double ns_goodness = -HUGE_VAL, norm = normi;
126 nk = i+nsmax > n ? n : i+nsmax; /* max k for this subspace */
127 for (k = i; k < i+nsmin-1; ++k) norm += fabs(dx[p[k]]);
129 for (k = i+nsmin-1; k < nk; ++k) {
131 norm += fabs(dx[p[k]]);
132 /* remaining subspaces must be big enough to partition */
133 if (n-(k+1) < nsmin) continue;
134 /* maximize figure of merit defined by Rowan thesis:
135 look for sudden drops in average |dx| */
137 goodness = norm/(k+1) - (normdx-norm)/(n-(k+1));
140 if (goodness > ns_goodness) {
141 ns_goodness = goodness;
145 for (k = i; k < i+ns; ++k) normi += fabs(dx[p[k]]);
146 /* do nelder-mead on subspace of dimension ns starting w/i */
148 for (k = i; k < i+ns; ++k) {
150 xsstep[k-i] = xstep[p[k]];
155 nevals = stop->nevals;
156 ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
157 xsstep, stop, psi, scratch, &fdiff);
158 if (fdiff > fdiff_max) fdiff_max = fdiff;
160 printf("%d NM iterations for (%d,%d) subspace\n",
161 stop->nevals - nevals, sd.is, ns);
162 for (k = i; k < i+ns; ++k) x[p[k]] = xs[k-i];
163 if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
164 if (ret != NLOPT_XTOL_REACHED) goto done;
166 /* nelder-mead on last subspace */
170 xs[i-sd.is] = x[p[i]];
171 xsstep[i-sd.is] = xstep[p[i]];
172 lbs[i-sd.is] = lb[p[i]];
173 ubs[i-sd.is] = ub[p[i]];
176 nevals = stop->nevals;
177 ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
178 xsstep, stop, psi, scratch, &fdiff);
179 if (fdiff > fdiff_max) fdiff_max = fdiff;
181 printf("sbplx: %d NM iterations for (%d,%d) subspace\n",
182 stop->nevals - nevals, sd.is, ns);
183 for (i = sd.is; i < n; ++i) x[p[i]] = xs[i-sd.is];
184 if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
185 if (ret != NLOPT_XTOL_REACHED) goto done;
187 /* termination tests: */
188 if (nlopt_stop_ftol(stop, *minf, *minf + fdiff_max)) {
189 ret = NLOPT_FTOL_REACHED;
192 if (nlopt_stop_x(stop, x, xprev)) {
194 /* as explained in Rowan's thesis, it is important
195 to check |xstep| as well as |x-xprev|, since if
196 the step size is too large (in early iterations),
197 the inner Nelder-Mead may not make much progress */
198 for (j = 0; j < n; ++j)
199 if (fabs(xstep[j]) * psi > stop->xtol_abs[j]
200 && fabs(xstep[j]) * psi > stop->xtol_rel * fabs(x[j]))
203 ret = NLOPT_XTOL_REACHED;
208 /* compute change in optimal point */
209 for (i = 0; i < n; ++i) dx[i] = x[i] - xprev[i];
211 /* setting stepsizes */
217 double stepnorm = 0, dxnorm = 0;
218 for (i = 0; i < n; ++i) {
219 stepnorm += fabs(xstep[i]);
220 dxnorm += fabs(dx[i]);
222 scale = dxnorm / stepnorm;
223 if (scale < omega) scale = omega;
224 if (scale > 1/omega) scale = 1/omega;
227 printf("sbplx: stepsize scale factor = %g\n", scale);
228 for (i = 0; i < n; ++i)
229 xstep[i] = (dx[i] == 0) ? -(xstep[i] * scale)
230 : copysign(xstep[i] * scale, dx[i]);