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_;
166 double infeasibility;
168 sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
169 if (!sigma) return NLOPT_OUT_OF_MEMORY;
174 xprevprev = xprev + n;
175 fcval = xprevprev + n;
176 fcval_cur = fcval + m;
177 rhoc = fcval_cur + m;
180 dual_ub = dual_lb + m;
183 dfcdx_cur = dfcdx + m*n;
197 for (j = 0; j < n; ++j) {
198 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
199 sigma[j] = 1.0; /* arbitrary default */
201 sigma[j] = 0.5 * (ub[j] - lb[j]);
204 for (i = 0; i < m; ++i) {
206 dual_lb[i] = y[i] = 0.0;
207 dual_ub[i] = HUGE_VAL;
210 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
212 memcpy(xcur, x, sizeof(double) * n);
214 feasible = 1; infeasibility = 0;
215 for (i = 0; i < m; ++i) {
216 fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
217 feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i]));
218 if (fcval[i] > infeasibility) infeasibility = fcval[i];
220 /* For non-feasible initial points, set a finite (large)
221 upper-bound on the dual variables. What this means is that,
222 if no feasible solution is found from the dual problem, it
223 will minimize the dual objective with the unfeasible
224 constraint weighted by 1e40 -- basically, minimizing the
225 unfeasible constraint until it becomes feasible or until we at
226 least obtain a step towards a feasible point.
228 Svanberg suggested a different approach in his 1987 paper, basically
229 introducing additional penalty variables for unfeasible constraints,
230 but this is easier to implement and at least as efficient. */
232 for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
234 while (1) { /* outer iterations */
236 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
237 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
238 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
239 if (ret != NLOPT_SUCCESS) goto done;
240 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
241 memcpy(xprev, xcur, sizeof(double) * n);
243 while (1) { /* inner iterations */
244 double min_dual, infeasibility_cur;
245 int feasible_cur, inner_done, save_verbose;
246 int new_infeasible_constraint;
249 /* solve dual problem */
250 dd.rho = rho; dd.count = 0;
252 save_verbose = mma_verbose;
254 reti = nlopt_minimize(
255 dual_alg, m, dual_func, &dd,
256 dual_lb, dual_ub, y, &min_dual,
257 -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
258 stop->maxtime - (nlopt_seconds() - stop->start));
259 mma_verbose = save_verbose;
260 if (reti == NLOPT_FAILURE && dual_alg != NLOPT_LD_MMA) {
261 /* LBFGS etc. converge quickly but are sometimes
262 very finicky if there are any rounding errors in
263 the gradient, etcetera; if it fails, try again
264 with MMA called recursively for the dual */
265 dual_alg = NLOPT_LD_MMA;
267 printf("MMA: switching to recursive MMA for dual\n");
270 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
275 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
277 printf("MMA dual converged in %d iterations to g=%g:\n",
279 for (i = 0; i < MIN(mma_verbose, m); ++i)
280 printf(" MMA y[%d]=%g, gc[%d]=%g\n",
281 i, y[i], i, dd.gcval[i]);
284 fcur = f(n, xcur, dfdx_cur, f_data);
286 feasible_cur = 1; infeasibility_cur = 0;
287 new_infeasible_constraint = 0;
288 inner_done = dd.gval >= fcur;
289 for (i = 0; i < m; ++i) {
290 fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
291 fc_data + fc_datum_size * i);
292 if (!isnan(fcval_cur[i])) {
293 feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
294 if (!isnan(fcval[i]))
295 inner_done = inner_done &&
296 (dd.gcval[i] >= fcval_cur[i]);
297 else if (fcval_cur[i] > 0)
298 new_infeasible_constraint = 1;
299 if (fcval_cur[i] > infeasibility_cur)
300 infeasibility_cur = fcval_cur[i];
304 if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
305 || (!feasible && infeasibility_cur < infeasibility)) {
306 if (mma_verbose && !feasible_cur)
307 printf("MMA - using infeasible point?\n");
308 dd.fval = *minf = fcur;
309 infeasibility = infeasibility_cur;
310 memcpy(fcval, fcval_cur, sizeof(double)*m);
311 memcpy(x, xcur, sizeof(double)*n);
312 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
313 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
315 /* once we have reached a feasible solution, the
316 algorithm should never make the solution infeasible
317 again (if inner_done), although the constraints may
318 be violated slightly by rounding errors etc. so we
319 must be a little careful about checking feasibility */
321 if (!feasible) /* reset upper bounds to infinity */
322 for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
325 else if (new_infeasible_constraint) feasible = 0;
328 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
329 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
330 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
331 if (ret != NLOPT_SUCCESS) goto done;
333 if (inner_done) break;
336 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
337 for (i = 0; i < m; ++i)
338 if (!isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i])
341 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
345 printf("MMA inner iteration: rho -> %g\n", rho);
346 for (i = 0; i < MIN(mma_verbose, m); ++i)
347 printf(" MMA rhoc[%d] -> %g\n", i,rhoc[i]);
350 if (nlopt_stop_ftol(stop, fcur, fprev))
351 ret = NLOPT_FTOL_REACHED;
352 if (nlopt_stop_x(stop, xcur, xprev))
353 ret = NLOPT_XTOL_REACHED;
354 if (ret != NLOPT_SUCCESS) goto done;
356 /* update rho and sigma for iteration k+1 */
357 rho = MAX(0.1 * rho, MMA_RHOMIN);
359 printf("MMA outer iteration: rho -> %g\n", rho);
360 for (i = 0; i < m; ++i)
361 rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
362 for (i = 0; i < MIN(mma_verbose, m); ++i)
363 printf(" MMA rhoc[%d] -> %g\n", i, rhoc[i]);
365 for (j = 0; j < n; ++j) {
366 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
367 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
369 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
370 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
371 sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
374 for (j = 0; j < MIN(mma_verbose, n); ++j)
375 printf(" MMA sigma[%d] -> %g\n",