1 /* Copyright (c) 2007-2008 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.
31 int mma_verbose = 0; /* > 0 for verbose output */
33 #define MIN(a,b) ((a) < (b) ? (a) : (b))
34 #define MAX(a,b) ((a) > (b) ? (a) : (b))
37 static int my_isnan(double x) { return x != x; }
38 # define isnan my_isnan
41 /* magic minimum value for rho in MMA ... the 2002 paper says it should
42 be a "fixed, strictly positive `small' number, e.g. 1e-5"
43 ... grrr, I hate these magic numbers, which seem like they
44 should depend on the objective function in some way ... in particular,
45 note that rho is dimensionful (= dimensions of objective function) */
46 #define MMA_RHOMIN 1e-5
48 /***********************************************************************/
49 /* function for MMA's dual solution of the approximate problem */
52 int count; /* evaluation count, incremented each call */
53 int n; /* must be set on input to dimension of x */
54 const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
55 const double *dfcdx; /* m-by-n array of fc gradients */
56 double fval, rho; /* must be set on input */
57 const double *fcval, *rhoc; /* arrays of length m */
58 double *xcur; /* array of length n, output each time */
59 double gval, wval, *gcval; /* output each time (array length m) */
62 static double sqr(double x) { return x * x; }
64 static double dual_func(int m, const double *y, double *grad, void *d_)
66 dual_data *d = (dual_data *) d_;
68 const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
70 const double *dfcdx = d->dfcdx;
71 double rho = d->rho, fval = d->fval;
72 const double *rhoc = d->rhoc, *fcval = d->fcval;
73 double *xcur = d->xcur;
74 double *gcval = d->gcval;
82 for (i = 0; i < m; ++i)
83 val += y[i] * (gcval[i] = isnan(fcval[i]) ? 0 : fcval[i]);
85 for (j = 0; j < n; ++j) {
86 double u, v, dx, denominv, c, sigma2, dx2;
88 /* first, compute xcur[j] for y. Because this objective is
89 separable, we can minimize over x analytically, and the minimum
90 dx is given by the solution of a quadratic equation:
91 u dx^2 + 2 v sigma^2 dx + u sigma^2 = 0
92 where u and v are defined by the sums below. Because of
93 the definitions, it is guaranteed that |u/v| <= sigma,
94 and it follows that the only dx solution with |dx| <= sigma
96 (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2))
97 (which goes to zero as u -> 0). */
100 v = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
101 for (i = 0; i < m; ++i) if (!isnan(fcval[i])) {
102 u += dfcdx[i*n + j] * y[i];
103 v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
105 u *= (sigma2 = sqr(sigma[j]));
106 if (fabs(u) < 1e-3 * (v*sigma[j])) { /* Taylor exp. for small u */
107 double a = u / (v*sigma[j]);
108 dx = -sigma[j] * (0.5 * a + 0.125 * a*a*a);
111 dx = (v/u)*sigma2 * (-1 + sqrt(fabs(1 - sqr(u/(v*sigma[j])))));
113 if (xcur[j] > ub[j]) xcur[j] = ub[j];
114 else if (xcur[j] < lb[j]) xcur[j] = lb[j];
115 if (xcur[j] > x[j]+0.9*sigma[j]) xcur[j] = x[j]+0.9*sigma[j];
116 else if (xcur[j] < x[j]-0.9*sigma[j]) xcur[j] = x[j]-0.9*sigma[j];
119 /* function value: */
121 denominv = 1.0 / (sigma2 - dx2);
122 val += (u * dx + v * dx2) * denominv;
124 /* update gval, wval, gcval (approximant functions) */
126 d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2)
128 d->wval += 0.5 * dx2 * denominv;
129 for (i = 0; i < m; ++i) if (!isnan(fcval[i]))
130 gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j]
131 + 0.5*rhoc[i]) * dx2)
135 /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
136 we only need the partial derivative with respect to y, and
137 we negate because we are maximizing: */
138 if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
142 /***********************************************************************/
144 /* note that we implement a hidden feature not in the standard
145 nlopt_minimize_constrained interface: whenever the constraint
146 function returns NaN, that constraint becomes inactive. */
148 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
149 int m, nlopt_func fc,
150 void *fc_data_, ptrdiff_t fc_datum_size,
151 const double *lb, const double *ub, /* bounds */
152 double *x, /* in: initial guess, out: minimizer */
154 nlopt_stopping *stop,
155 nlopt_algorithm dual_alg,
156 double dual_tolrel, int dual_maxeval)
158 nlopt_result ret = NLOPT_SUCCESS;
159 double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
160 double *dfcdx, *dfcdx_cur;
161 double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
163 char *fc_data = (char *) fc_data_;
167 sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
168 if (!sigma) return NLOPT_OUT_OF_MEMORY;
173 xprevprev = xprev + n;
174 fcval = xprevprev + n;
175 fcval_cur = fcval + m;
176 rhoc = fcval_cur + m;
179 dual_ub = dual_lb + m;
182 dfcdx_cur = dfcdx + m*n;
196 for (j = 0; j < n; ++j) {
197 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
198 sigma[j] = 1.0; /* arbitrary default */
200 sigma[j] = 0.5 * (ub[j] - lb[j]);
203 for (i = 0; i < m; ++i) {
205 dual_lb[i] = y[i] = 0.0;
206 dual_ub[i] = HUGE_VAL;
209 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
211 memcpy(xcur, x, sizeof(double) * n);
214 for (i = 0; i < m; ++i) {
215 fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
216 feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i]));
218 /* For non-feasible initial points, set a finite (large)
219 upper-bound on the dual variables. What this means is that,
220 if no feasible solution is found from the dual problem, it
221 will minimize the dual objective with the unfeasible
222 constraint weighted by 1e40 -- basically, minimizing the
223 unfeasible constraint until it becomes feasible or until we at
224 least obtain a step towards a feasible point.
226 Svanberg suggested a different approach in his 1987 paper, basically
227 introducing additional penalty variables for unfeasible constraints,
228 but this is easier to implement and at least as efficient. */
230 for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
232 while (1) { /* outer iterations */
234 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
235 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
236 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
237 if (ret != NLOPT_SUCCESS) goto done;
238 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
239 memcpy(xprev, xcur, sizeof(double) * n);
241 while (1) { /* inner iterations */
243 int feasible_cur, inner_done, save_verbose;
244 int new_infeasible_constraint;
247 /* solve dual problem */
248 dd.rho = rho; dd.count = 0;
250 save_verbose = mma_verbose;
252 reti = nlopt_minimize(
253 dual_alg, m, dual_func, &dd,
254 dual_lb, dual_ub, y, &min_dual,
255 -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
256 stop->maxtime - (nlopt_seconds() - stop->start));
257 mma_verbose = save_verbose;
258 if (reti == NLOPT_FAILURE && dual_alg != NLOPT_LD_MMA) {
259 /* LBFGS etc. converge quickly but are sometimes
260 very finicky if there are any rounding errors in
261 the gradient, etcetera; if it fails, try again
262 with MMA called recursively for the dual */
263 dual_alg = NLOPT_LD_MMA;
265 printf("MMA: switching to recursive MMA for dual\n");
268 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
273 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
275 printf("MMA dual converged in %d iterations to g=%g:\n",
277 for (i = 0; i < MIN(mma_verbose, m); ++i)
278 printf(" MMA y[%d]=%g, gc[%d]=%g\n",
279 i, y[i], i, dd.gcval[i]);
282 fcur = f(n, xcur, dfdx_cur, f_data);
285 new_infeasible_constraint = 0;
286 inner_done = dd.gval >= fcur;
287 for (i = 0; i < m; ++i) {
288 fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
289 fc_data + fc_datum_size * i);
290 if (!isnan(fcval_cur[i])) {
291 feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
292 if (!isnan(fcval[i]))
293 inner_done = inner_done &&
294 (dd.gcval[i] >= fcval_cur[i]);
295 else if (fcval_cur[i] > 0)
296 new_infeasible_constraint = 1;
300 if (fcur < *minf && (inner_done || feasible_cur || !feasible)) {
301 if (mma_verbose && !feasible_cur)
302 printf("MMA - using infeasible point?\n");
303 dd.fval = *minf = fcur;
304 memcpy(fcval, fcval_cur, sizeof(double)*m);
305 memcpy(x, xcur, sizeof(double)*n);
306 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
307 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
309 /* once we have reached a feasible solution, the
310 algorithm should never make the solution infeasible
311 again (if inner_done), although the constraints may
312 be violated slightly by rounding errors etc. so we
313 must be a little careful about checking feasibility */
315 if (!feasible) /* reset upper bounds to infinity */
316 for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
319 else if (new_infeasible_constraint) feasible = 0;
322 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
323 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
324 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
325 if (ret != NLOPT_SUCCESS) goto done;
327 if (inner_done) break;
330 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
331 for (i = 0; i < m; ++i)
332 if (!isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i])
335 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
339 printf("MMA inner iteration: rho -> %g\n", rho);
340 for (i = 0; i < MIN(mma_verbose, m); ++i)
341 printf(" MMA rhoc[%d] -> %g\n", i,rhoc[i]);
344 if (nlopt_stop_ftol(stop, fcur, fprev))
345 ret = NLOPT_FTOL_REACHED;
346 if (nlopt_stop_x(stop, xcur, xprev))
347 ret = NLOPT_XTOL_REACHED;
348 if (ret != NLOPT_SUCCESS) goto done;
350 /* update rho and sigma for iteration k+1 */
351 rho = MAX(0.1 * rho, MMA_RHOMIN);
353 printf("MMA outer iteration: rho -> %g\n", rho);
354 for (i = 0; i < m; ++i)
355 rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
356 for (i = 0; i < MIN(mma_verbose, m); ++i)
357 printf(" MMA rhoc[%d] -> %g\n", i, rhoc[i]);
359 for (j = 0; j < n; ++j) {
360 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
361 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
363 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
364 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
365 sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
368 for (j = 0; j < MIN(mma_verbose, n); ++j)
369 printf(" MMA sigma[%d] -> %g\n",