1 /* Copyright (c) 2007-2011 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_LD_CCSAQ, 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, 1e-12);
307 if (ret < 0) goto done;
308 ret = nlopt_set_maxeval(pre_opt, 100000);
309 if (ret < 0) goto done;
312 for (j = 0; j < n; ++j) {
313 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
314 sigma[j] = 1.0; /* arbitrary default */
316 sigma[j] = 0.5 * (ub[j] - lb[j]);
319 for (i = 0; i < m; ++i) {
321 dual_lb[i] = y[i] = 0.0;
322 dual_ub[i] = HUGE_VAL;
325 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
327 memcpy(xcur, x, sizeof(double) * n);
328 if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
330 feasible = 1; infeasibility = 0;
331 for (i = ifc = 0; ifc < mfc; ++ifc) {
332 nlopt_eval_constraint(fcval + i, dfcdx + i*n,
335 if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
337 for (i = 0; i < m; ++i) {
338 feasible = feasible && fcval[i] <= 0;
339 if (fcval[i] > infeasibility) infeasibility = fcval[i];
341 /* For non-feasible initial points, set a finite (large)
342 upper-bound on the dual variables. What this means is that,
343 if no feasible solution is found from the dual problem, it
344 will minimize the dual objective with the unfeasible
345 constraint weighted by 1e40 -- basically, minimizing the
346 unfeasible constraint until it becomes feasible or until we at
347 least obtain a step towards a feasible point.
349 Svanberg suggested a different approach in his 1987 paper, basically
350 introducing additional penalty variables for unfeasible constraints,
351 but this is easier to implement and at least as efficient. */
353 for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
355 nlopt_set_min_objective(dual_opt, dual_func, &dd);
356 nlopt_set_lower_bounds(dual_opt, dual_lb);
357 nlopt_set_upper_bounds(dual_opt, dual_ub);
358 nlopt_set_stopval(dual_opt, -HUGE_VAL);
359 nlopt_remove_inequality_constraints(dual_opt);
360 nlopt_remove_equality_constraints(dual_opt);
362 while (1) { /* outer iterations */
364 if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
365 else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
366 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
367 else if (feasible && *minf < stop->minf_max)
368 ret = NLOPT_MINF_MAX_REACHED;
369 if (ret != NLOPT_SUCCESS) goto done;
370 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
371 memcpy(xprev, xcur, sizeof(double) * n);
373 while (1) { /* inner iterations */
374 double min_dual, infeasibility_cur;
375 int feasible_cur, inner_done;
376 unsigned save_verbose;
380 /* solve dual problem */
381 dd.rho = rho; dd.count = 0;
382 save_verbose = ccsa_verbose;
383 ccsa_verbose = 0; /* no recursive verbosity */
384 reti = nlopt_optimize_limited(dual_opt, y, &min_dual,
389 ccsa_verbose = save_verbose;
390 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
395 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
399 for (j = 0; j < n; ++j) {
400 pre_lb[j] = MAX(lb[j], x[j] - sigma[j]);
401 pre_ub[j] = MIN(ub[j], x[j] + sigma[j]);
404 nlopt_set_lower_bounds(pre_opt, pre_lb);
405 nlopt_set_upper_bounds(pre_opt, pre_ub);
407 dd.rho = rho; dd.count = 0;
408 save_verbose = ccsa_verbose;
409 ccsa_verbose = 0; /* no recursive verbosity */
410 reti = nlopt_optimize_limited(pre_opt, xcur, &pre_min,
414 ccsa_verbose = save_verbose;
415 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
420 /* evaluate final xcur etc */
421 dd.gval = g0(n, xcur, NULL, &dd);
422 gi(m, dd.gcval, n, xcur, NULL, &dd);
426 printf("CCSA dual converged in %d iters to g=%g:\n",
428 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
429 printf(" CCSA y[%d]=%g, gc[%d]=%g\n",
430 i, y[i], i, dd.gcval[i]);
433 fcur = f(n, xcur, dfdx_cur, f_data);
435 if (nlopt_stop_forced(stop)) {
436 ret = NLOPT_FORCED_STOP; goto done; }
437 feasible_cur = 1; infeasibility_cur = 0;
438 inner_done = dd.gval >= fcur;
439 for (i = ifc = 0; ifc < mfc; ++ifc) {
440 nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
443 if (nlopt_stop_forced(stop)) {
444 ret = NLOPT_FORCED_STOP; goto done; }
446 for (i = ifc = 0; ifc < mfc; ++ifc) {
447 unsigned i0 = i, inext = i + fc[ifc].m;
448 for (; i < inext; ++i) {
449 feasible_cur = feasible_cur
450 && fcval_cur[i] <= fc[ifc].tol[i-i0];
451 inner_done = inner_done &&
452 (dd.gcval[i] >= fcval_cur[i]);
453 if (fcval_cur[i] > infeasibility_cur)
454 infeasibility_cur = fcval_cur[i];
458 if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
459 || (!feasible && infeasibility_cur < infeasibility)) {
460 if (ccsa_verbose && !feasible_cur)
461 printf("CCSA - using infeasible point?\n");
462 dd.fval = *minf = fcur;
463 infeasibility = infeasibility_cur;
464 memcpy(fcval, fcval_cur, sizeof(double)*m);
465 memcpy(x, xcur, sizeof(double)*n);
466 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
467 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
469 /* once we have reached a feasible solution, the
470 algorithm should never make the solution infeasible
471 again (if inner_done), although the constraints may
472 be violated slightly by rounding errors etc. so we
473 must be a little careful about checking feasibility */
474 if (infeasibility_cur == 0) {
475 if (!feasible) { /* reset upper bounds to infin. */
476 for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
477 nlopt_set_upper_bounds(dual_opt, dual_ub);
483 if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
484 else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
485 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
486 else if (feasible && *minf < stop->minf_max)
487 ret = NLOPT_MINF_MAX_REACHED;
488 if (ret != NLOPT_SUCCESS) goto done;
490 if (inner_done) break;
493 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
494 for (i = 0; i < m; ++i)
495 if (fcval_cur[i] > dd.gcval[i])
498 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
502 printf("CCSA inner iteration: rho -> %g\n", rho);
503 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
504 printf(" CCSA rhoc[%d] -> %g\n", i,rhoc[i]);
507 if (nlopt_stop_ftol(stop, fcur, fprev))
508 ret = NLOPT_FTOL_REACHED;
509 if (nlopt_stop_x(stop, xcur, xprev))
510 ret = NLOPT_XTOL_REACHED;
511 if (ret != NLOPT_SUCCESS) goto done;
513 /* update rho and sigma for iteration k+1 */
514 rho = MAX(0.1 * rho, CCSA_RHOMIN);
516 printf("CCSA outer iteration: rho -> %g\n", rho);
517 for (i = 0; i < m; ++i)
518 rhoc[i] = MAX(0.1 * rhoc[i], CCSA_RHOMIN);
519 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
520 printf(" CCSA rhoc[%d] -> %g\n", i, rhoc[i]);
522 for (j = 0; j < n; ++j) {
523 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
524 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
526 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
527 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
528 /* use a smaller lower bound than Svanberg's
529 0.01*(ub-lb), which seems unnecessarily large */
530 sigma[j] = MAX(sigma[j], 1e-8*(ub[j]-lb[j]));
533 for (j = 0; j < MIN(ccsa_verbose, n); ++j)
534 printf(" CCSA sigma[%d] -> %g\n",
540 nlopt_destroy(pre_opt);
541 if (dd.scratch) free(dd.scratch);