1 /* Copyright (c) 2007-2012 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.
23 /* In this file we implement Svanberg's CCSA algorithm with the
24 simple linear approximation + quadratic penalty function.
26 We also allow the user to specify an optional "preconditioner": an
27 approximate Hessian (which must be symmetric & positive
28 semidefinite) that can be added into the approximation. [X. Liang
29 and I went through the convergence proof in Svanberg's paper
30 and it does not seem to be modified by such preconditioning, as
31 long as the preconditioner eigenvalues are bounded above for all x.]
33 For the non-preconditioned case the trust-region subproblem is
34 separable and can be solved by a dual method. For the preconditioned
35 case the subproblem is still convex but in general is non-separable
36 so we solve by calling the same algorithm recursively, under the
37 assumption that the subproblem objective is cheap to evaluate.
46 #include "nlopt-util.h"
48 unsigned ccsa_verbose = 0; /* > 0 for verbose output */
50 #define MIN(a,b) ((a) < (b) ? (a) : (b))
51 #define MAX(a,b) ((a) > (b) ? (a) : (b))
53 /* magic minimum value for rho in CCSA ... the 2002 paper says it should
54 be a "fixed, strictly positive `small' number, e.g. 1e-5"
55 ... grrr, I hate these magic numbers, which seem like they
56 should depend on the objective function in some way ... in particular,
57 note that rho is dimensionful (= dimensions of objective function) */
58 #define CCSA_RHOMIN 1e-5
60 /***********************************************************************/
61 /* function for CCSA's dual solution of the approximate problem */
64 int count; /* evaluation count, incremented each call */
65 unsigned n; /* must be set on input to dimension of x */
66 const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
67 const double *dfcdx; /* m-by-n array of fc gradients */
68 double fval, rho; /* must be set on input */
69 const double *fcval, *rhoc; /* arrays of length m */
70 double *xcur; /* array of length n, output each time */
71 double gval, wval, *gcval; /* output each time (array length m) */
72 nlopt_precond pre; void *pre_data;
73 nlopt_precond *prec; void **prec_data; /* length = # constraints */
74 double *scratch; /* length = 2*n */
77 static double sqr(double x) { return x * x; }
79 static double dual_func(unsigned m, const double *y, double *grad, void *d_)
81 dual_data *d = (dual_data *) d_;
83 const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
85 const double *dfcdx = d->dfcdx;
86 double rho = d->rho, fval = d->fval;
87 const double *rhoc = d->rhoc, *fcval = d->fcval;
88 double *xcur = d->xcur;
89 double *gcval = d->gcval;
97 for (i = 0; i < m; ++i)
98 val += y[i] * (gcval[i] = fcval[i]);
100 for (j = 0; j < n; ++j) {
101 double u, v, dx, sigma2, dx2, dx2sig;
103 /* first, compute xcur[j] = x+dx for y. Because this objective is
104 separable, we can minimize over x analytically, and the minimum
105 dx is given by the solution of a linear equation
106 u dx + v sigma^2 = 0, i.e. dx = -sigma^2 v/u
107 where u and v are defined by the sums below. However,
108 we also have to check that |dx| <= sigma and that
111 if (sigma[j] == 0) { /* special case for lb[i] == ub[i] dims, dx=0 */
118 for (i = 0; i < m; ++i) {
120 v += dfcdx[i*n + j] * y[i];
122 dx = -(sigma2 = sqr(sigma[j])) * v/u;
124 /* if dx is out of bounds, we are guaranteed by convexity
125 that the minimum is at the bound on the side of dx */
126 if (fabs(dx) > sigma[j]) dx = copysign(sigma[j], dx);
128 if (xcur[j] > ub[j]) xcur[j] = ub[j];
129 else if (xcur[j] < lb[j]) xcur[j] = lb[j];
132 /* function value: */
134 val += v * dx + 0.5 * u * dx2 / sigma2;
136 /* update gval, wval, gcval (approximant functions) */
137 d->gval += dfdx[j] * dx + rho * (dx2sig = 0.5*dx2/sigma2);
139 for (i = 0; i < m; ++i)
140 gcval[i] += dfcdx[i*n+j] * dx + rhoc[i] * dx2sig;
143 /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
144 we only need the partial derivative with respect to y, and
145 we negate because we are maximizing: */
146 if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
150 /***********************************************************************/
152 /* compute g(x - x0) and its gradient */
153 static double gfunc(unsigned n, double f, const double *dfdx,
154 double rho, const double *sigma,
156 nlopt_precond pre, void *pre_data, double *scratch,
157 const double *x, double *grad)
159 double *dx = scratch, *Hdx = scratch + n;
163 for (j = 0; j < n; ++j) {
164 double sigma2inv = 1.0 / sqr(sigma[j]);
165 dx[j] = x[j] - x0[j];
166 val += dfdx[j] * dx[j] + (0.5*rho) * sqr(dx[j]) * sigma2inv;
167 if (grad) grad[j] = dfdx[j] + rho * dx[j] * sigma2inv;
171 pre(n, x0, dx, Hdx, pre_data);
172 for (j = 0; j < n; ++j)
173 val += 0.5 * dx[j] * Hdx[j];
175 for (j = 0; j < n; ++j)
182 static double g0(unsigned n, const double *x, double *grad, void *d_)
184 dual_data *d = (dual_data *) d_;
186 return gfunc(n, d->fval, d->dfdx, d->rho, d->sigma,
188 d->pre, d->pre_data, d->scratch,
193 static void gi(unsigned m, double *result,
194 unsigned n, const double *x, double *grad, void *d_)
196 dual_data *d = (dual_data *) d_;
198 for (i = 0; i < m; ++i)
199 result[i] = gfunc(n, d->fcval[i], d->dfcdx + i*n, d->rhoc[i],
202 d->prec ? d->prec[i] : NULL,
203 d->prec_data ? d->prec_data[i] : NULL,
209 /***********************************************************************/
211 nlopt_result ccsa_quadratic_minimize(
212 unsigned n, nlopt_func f, void *f_data,
213 unsigned m, nlopt_constraint *fc,
217 const double *lb, const double *ub, /* bounds */
218 double *x, /* in: initial guess, out: minimizer */
220 nlopt_stopping *stop,
223 nlopt_result ret = NLOPT_SUCCESS;
224 double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
225 double *dfcdx, *dfcdx_cur;
226 double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
227 double *pre_lb, *pre_ub;
228 unsigned i, ifc, j, k = 0;
231 double infeasibility;
234 nlopt_opt pre_opt = NULL;
236 m = nlopt_count_constraints(mfc = m, fc);
237 if (nlopt_get_dimension(dual_opt) != m) return NLOPT_INVALID_ARGS;
238 sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
239 if (!sigma) return NLOPT_OUT_OF_MEMORY;
244 xprevprev = xprev + n;
245 fcval = xprevprev + n;
246 fcval_cur = fcval + m;
247 rhoc = fcval_cur + m;
250 dual_ub = dual_lb + m;
253 dfcdx_cur = dfcdx + m*n;
266 dd.pre = pre; dd.pre_data = f_data;
267 dd.prec = NULL; dd.prec_data = NULL;
271 dd.prec = (nlopt_precond *) malloc(sizeof(nlopt_precond) * m);
272 dd.prec_data = (void **) malloc(sizeof(void *) * m);
273 if (!dd.prec || !dd.prec_data) {
274 ret = NLOPT_OUT_OF_MEMORY;
277 for (i = ifc = 0; ifc < mfc; ++ifc) {
278 unsigned inext = i + fc[ifc].m;
279 for (; i < inext; ++i) {
280 dd.prec[i] = fc[ifc].pre;
281 dd.prec_data[i] = fc[ifc].f_data;
286 no_precond = pre == NULL;
288 for (i = 0; i < m; ++i)
289 no_precond = no_precond && dd.prec[i] == NULL;
292 dd.scratch = (double*) malloc(sizeof(double) * (4*n));
295 return NLOPT_OUT_OF_MEMORY;
297 pre_lb = dd.scratch + 2*n;
300 pre_opt = nlopt_create(nlopt_get_algorithm(dual_opt), n);
301 if (!pre_opt) { ret = NLOPT_FAILURE; goto done; }
302 ret = nlopt_set_min_objective(pre_opt, g0, &dd);
303 if (ret < 0) goto done;
304 ret = nlopt_add_inequality_mconstraint(pre_opt, m, gi, &dd, NULL);
305 if (ret < 0) goto done;
306 ret = nlopt_set_ftol_rel(pre_opt, nlopt_get_ftol_rel(dual_opt));
307 if (ret < 0) goto done;
308 ret = nlopt_set_ftol_abs(pre_opt, nlopt_get_ftol_abs(dual_opt));
309 if (ret < 0) goto done;
310 ret = nlopt_set_maxeval(pre_opt, nlopt_get_maxeval(dual_opt));
311 if (ret < 0) goto done;
314 for (j = 0; j < n; ++j) {
315 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
316 sigma[j] = 1.0; /* arbitrary default */
318 sigma[j] = 0.5 * (ub[j] - lb[j]);
321 for (i = 0; i < m; ++i) {
323 dual_lb[i] = y[i] = 0.0;
324 dual_ub[i] = HUGE_VAL;
327 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
329 memcpy(xcur, x, sizeof(double) * n);
330 if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
332 feasible = 1; infeasibility = 0;
333 for (i = ifc = 0; ifc < mfc; ++ifc) {
334 nlopt_eval_constraint(fcval + i, dfcdx + i*n,
337 if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
339 for (i = 0; i < m; ++i) {
340 feasible = feasible && fcval[i] <= 0;
341 if (fcval[i] > infeasibility) infeasibility = fcval[i];
343 /* For non-feasible initial points, set a finite (large)
344 upper-bound on the dual variables. What this means is that,
345 if no feasible solution is found from the dual problem, it
346 will minimize the dual objective with the unfeasible
347 constraint weighted by 1e40 -- basically, minimizing the
348 unfeasible constraint until it becomes feasible or until we at
349 least obtain a step towards a feasible point.
351 Svanberg suggested a different approach in his 1987 paper, basically
352 introducing additional penalty variables for unfeasible constraints,
353 but this is easier to implement and at least as efficient. */
355 for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
357 nlopt_set_min_objective(dual_opt, dual_func, &dd);
358 nlopt_set_lower_bounds(dual_opt, dual_lb);
359 nlopt_set_upper_bounds(dual_opt, dual_ub);
360 nlopt_set_stopval(dual_opt, -HUGE_VAL);
361 nlopt_remove_inequality_constraints(dual_opt);
362 nlopt_remove_equality_constraints(dual_opt);
364 while (1) { /* outer iterations */
366 if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
367 else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
368 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
369 else if (feasible && *minf < stop->minf_max)
370 ret = NLOPT_MINF_MAX_REACHED;
371 if (ret != NLOPT_SUCCESS) goto done;
372 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
373 memcpy(xprev, xcur, sizeof(double) * n);
375 while (1) { /* inner iterations */
376 double min_dual, infeasibility_cur;
377 int feasible_cur, inner_done;
378 unsigned save_verbose;
382 /* solve dual problem */
383 dd.rho = rho; dd.count = 0;
384 save_verbose = ccsa_verbose;
385 ccsa_verbose = 0; /* no recursive verbosity */
386 reti = nlopt_optimize_limited(dual_opt, y, &min_dual,
391 ccsa_verbose = save_verbose;
392 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
397 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
401 for (j = 0; j < n; ++j) {
402 pre_lb[j] = MAX(lb[j], x[j] - sigma[j]);
403 pre_ub[j] = MIN(ub[j], x[j] + sigma[j]);
406 nlopt_set_lower_bounds(pre_opt, pre_lb);
407 nlopt_set_upper_bounds(pre_opt, pre_ub);
409 dd.rho = rho; dd.count = 0;
410 save_verbose = ccsa_verbose;
411 ccsa_verbose = 0; /* no recursive verbosity */
412 reti = nlopt_optimize_limited(pre_opt, xcur, &pre_min,
416 ccsa_verbose = save_verbose;
417 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
422 /* evaluate final xcur etc */
423 dd.gval = g0(n, xcur, NULL, &dd);
424 gi(m, dd.gcval, n, xcur, NULL, &dd);
428 printf("CCSA dual converged in %d iters to g=%g:\n",
430 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
431 printf(" CCSA y[%d]=%g, gc[%d]=%g\n",
432 i, y[i], i, dd.gcval[i]);
435 fcur = f(n, xcur, dfdx_cur, f_data);
437 if (nlopt_stop_forced(stop)) {
438 ret = NLOPT_FORCED_STOP; goto done; }
439 feasible_cur = 1; infeasibility_cur = 0;
440 inner_done = dd.gval >= fcur;
441 for (i = ifc = 0; ifc < mfc; ++ifc) {
442 nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
445 if (nlopt_stop_forced(stop)) {
446 ret = NLOPT_FORCED_STOP; goto done; }
448 for (i = ifc = 0; ifc < mfc; ++ifc) {
449 unsigned i0 = i, inext = i + fc[ifc].m;
450 for (; i < inext; ++i) {
451 feasible_cur = feasible_cur
452 && fcval_cur[i] <= fc[ifc].tol[i-i0];
453 inner_done = inner_done &&
454 (dd.gcval[i] >= fcval_cur[i]);
455 if (fcval_cur[i] > infeasibility_cur)
456 infeasibility_cur = fcval_cur[i];
460 if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
461 || (!feasible && infeasibility_cur < infeasibility)) {
462 if (ccsa_verbose && !feasible_cur)
463 printf("CCSA - using infeasible point?\n");
464 dd.fval = *minf = fcur;
465 infeasibility = infeasibility_cur;
466 memcpy(fcval, fcval_cur, sizeof(double)*m);
467 memcpy(x, xcur, sizeof(double)*n);
468 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
469 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
471 /* once we have reached a feasible solution, the
472 algorithm should never make the solution infeasible
473 again (if inner_done), although the constraints may
474 be violated slightly by rounding errors etc. so we
475 must be a little careful about checking feasibility */
476 if (infeasibility_cur == 0) {
477 if (!feasible) { /* reset upper bounds to infin. */
478 for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
479 nlopt_set_upper_bounds(dual_opt, dual_ub);
485 if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
486 else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
487 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
488 else if (feasible && *minf < stop->minf_max)
489 ret = NLOPT_MINF_MAX_REACHED;
490 if (ret != NLOPT_SUCCESS) goto done;
492 if (inner_done) break;
495 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
496 for (i = 0; i < m; ++i)
497 if (fcval_cur[i] > dd.gcval[i])
500 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
504 printf("CCSA inner iteration: rho -> %g\n", rho);
505 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
506 printf(" CCSA rhoc[%d] -> %g\n", i,rhoc[i]);
509 if (nlopt_stop_ftol(stop, fcur, fprev))
510 ret = NLOPT_FTOL_REACHED;
511 if (nlopt_stop_x(stop, xcur, xprev))
512 ret = NLOPT_XTOL_REACHED;
513 if (ret != NLOPT_SUCCESS) goto done;
515 /* update rho and sigma for iteration k+1 */
516 rho = MAX(0.1 * rho, CCSA_RHOMIN);
518 printf("CCSA outer iteration: rho -> %g\n", rho);
519 for (i = 0; i < m; ++i)
520 rhoc[i] = MAX(0.1 * rhoc[i], CCSA_RHOMIN);
521 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
522 printf(" CCSA rhoc[%d] -> %g\n", i, rhoc[i]);
524 for (j = 0; j < n; ++j) {
525 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
526 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
528 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
529 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
530 /* use a smaller lower bound than Svanberg's
531 0.01*(ub-lb), which seems unnecessarily large */
532 sigma[j] = MAX(sigma[j], 1e-8*(ub[j]-lb[j]));
535 for (j = 0; j < MIN(ccsa_verbose, n); ++j)
536 printf(" CCSA sigma[%d] -> %g\n",
542 nlopt_destroy(pre_opt);
543 if (dd.scratch) free(dd.scratch);