1 /* Copyright (c) 2007-2014 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) {
238 nlopt_stop_msg(stop, "dual optimizer has wrong dimension %d != %d",
239 nlopt_get_dimension(dual_opt), m);
240 return NLOPT_INVALID_ARGS;
242 sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
243 if (!sigma) return NLOPT_OUT_OF_MEMORY;
248 xprevprev = xprev + n;
249 fcval = xprevprev + n;
250 fcval_cur = fcval + m;
251 rhoc = fcval_cur + m;
254 dual_ub = dual_lb + m;
257 dfcdx_cur = dfcdx + m*n;
270 dd.pre = pre; dd.pre_data = f_data;
271 dd.prec = NULL; dd.prec_data = NULL;
275 dd.prec = (nlopt_precond *) malloc(sizeof(nlopt_precond) * m);
276 dd.prec_data = (void **) malloc(sizeof(void *) * m);
277 if (!dd.prec || !dd.prec_data) {
278 ret = NLOPT_OUT_OF_MEMORY;
281 for (i = ifc = 0; ifc < mfc; ++ifc) {
282 unsigned inext = i + fc[ifc].m;
283 for (; i < inext; ++i) {
284 dd.prec[i] = fc[ifc].pre;
285 dd.prec_data[i] = fc[ifc].f_data;
290 no_precond = pre == NULL;
292 for (i = 0; i < m; ++i)
293 no_precond = no_precond && dd.prec[i] == NULL;
296 dd.scratch = (double*) malloc(sizeof(double) * (4*n));
299 return NLOPT_OUT_OF_MEMORY;
301 pre_lb = dd.scratch + 2*n;
304 pre_opt = nlopt_create(nlopt_get_algorithm(dual_opt), n);
306 nlopt_stop_msg(stop, "failure creating precond. optimizer");
310 ret = nlopt_set_min_objective(pre_opt, g0, &dd);
311 if (ret < 0) goto done;
312 ret = nlopt_add_inequality_mconstraint(pre_opt, m, gi, &dd, NULL);
313 if (ret < 0) goto done;
314 ret = nlopt_set_ftol_rel(pre_opt, nlopt_get_ftol_rel(dual_opt));
315 if (ret < 0) goto done;
316 ret = nlopt_set_ftol_abs(pre_opt, nlopt_get_ftol_abs(dual_opt));
317 if (ret < 0) goto done;
318 ret = nlopt_set_maxeval(pre_opt, nlopt_get_maxeval(dual_opt));
319 if (ret < 0) goto done;
322 for (j = 0; j < n; ++j) {
323 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
324 sigma[j] = 1.0; /* arbitrary default */
326 sigma[j] = 0.5 * (ub[j] - lb[j]);
329 for (i = 0; i < m; ++i) {
331 dual_lb[i] = y[i] = 0.0;
332 dual_ub[i] = HUGE_VAL;
335 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
337 memcpy(xcur, x, sizeof(double) * n);
338 if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
340 feasible = 1; infeasibility = 0;
341 for (i = ifc = 0; ifc < mfc; ++ifc) {
342 nlopt_eval_constraint(fcval + i, dfcdx + i*n,
345 if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
347 for (i = 0; i < m; ++i) {
348 feasible = feasible && fcval[i] <= 0;
349 if (fcval[i] > infeasibility) infeasibility = fcval[i];
351 /* For non-feasible initial points, set a finite (large)
352 upper-bound on the dual variables. What this means is that,
353 if no feasible solution is found from the dual problem, it
354 will minimize the dual objective with the unfeasible
355 constraint weighted by 1e40 -- basically, minimizing the
356 unfeasible constraint until it becomes feasible or until we at
357 least obtain a step towards a feasible point.
359 Svanberg suggested a different approach in his 1987 paper, basically
360 introducing additional penalty variables for unfeasible constraints,
361 but this is easier to implement and at least as efficient. */
363 for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
365 nlopt_set_min_objective(dual_opt, dual_func, &dd);
366 nlopt_set_lower_bounds(dual_opt, dual_lb);
367 nlopt_set_upper_bounds(dual_opt, dual_ub);
368 nlopt_set_stopval(dual_opt, -HUGE_VAL);
369 nlopt_remove_inequality_constraints(dual_opt);
370 nlopt_remove_equality_constraints(dual_opt);
372 while (1) { /* outer iterations */
374 if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
375 else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
376 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
377 else if (feasible && *minf < stop->minf_max)
378 ret = NLOPT_MINF_MAX_REACHED;
379 if (ret != NLOPT_SUCCESS) goto done;
380 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
381 memcpy(xprev, xcur, sizeof(double) * n);
383 while (1) { /* inner iterations */
384 double min_dual, infeasibility_cur;
385 int feasible_cur, inner_done;
386 unsigned save_verbose;
390 /* solve dual problem */
391 dd.rho = rho; dd.count = 0;
392 save_verbose = ccsa_verbose;
393 ccsa_verbose = 0; /* no recursive verbosity */
394 reti = nlopt_optimize_limited(dual_opt, y, &min_dual,
399 ccsa_verbose = save_verbose;
400 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
405 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
409 for (j = 0; j < n; ++j) {
410 pre_lb[j] = MAX(lb[j], x[j] - sigma[j]);
411 pre_ub[j] = MIN(ub[j], x[j] + sigma[j]);
414 nlopt_set_lower_bounds(pre_opt, pre_lb);
415 nlopt_set_upper_bounds(pre_opt, pre_ub);
417 dd.rho = rho; dd.count = 0;
418 save_verbose = ccsa_verbose;
419 ccsa_verbose = 0; /* no recursive verbosity */
420 reti = nlopt_optimize_limited(pre_opt, xcur, &pre_min,
424 ccsa_verbose = save_verbose;
425 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
430 /* evaluate final xcur etc */
431 dd.gval = g0(n, xcur, NULL, &dd);
432 gi(m, dd.gcval, n, xcur, NULL, &dd);
436 printf("CCSA dual converged in %d iters to g=%g:\n",
438 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
439 printf(" CCSA y[%u]=%g, gc[%u]=%g\n",
440 i, y[i], i, dd.gcval[i]);
443 fcur = f(n, xcur, dfdx_cur, f_data);
445 if (nlopt_stop_forced(stop)) {
446 ret = NLOPT_FORCED_STOP; goto done; }
447 feasible_cur = 1; infeasibility_cur = 0;
448 inner_done = dd.gval >= fcur;
449 for (i = ifc = 0; ifc < mfc; ++ifc) {
450 nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
453 if (nlopt_stop_forced(stop)) {
454 ret = NLOPT_FORCED_STOP; goto done; }
456 for (i = ifc = 0; ifc < mfc; ++ifc) {
457 unsigned i0 = i, inext = i + fc[ifc].m;
458 for (; i < inext; ++i) {
459 feasible_cur = feasible_cur
460 && fcval_cur[i] <= fc[ifc].tol[i-i0];
461 inner_done = inner_done &&
462 (dd.gcval[i] >= fcval_cur[i]);
463 if (fcval_cur[i] > infeasibility_cur)
464 infeasibility_cur = fcval_cur[i];
468 if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
469 || (!feasible && infeasibility_cur < infeasibility)) {
470 if (ccsa_verbose && !feasible_cur)
471 printf("CCSA - using infeasible point?\n");
472 dd.fval = *minf = fcur;
473 infeasibility = infeasibility_cur;
474 memcpy(fcval, fcval_cur, sizeof(double)*m);
475 memcpy(x, xcur, sizeof(double)*n);
476 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
477 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
479 /* once we have reached a feasible solution, the
480 algorithm should never make the solution infeasible
481 again (if inner_done), although the constraints may
482 be violated slightly by rounding errors etc. so we
483 must be a little careful about checking feasibility */
484 if (infeasibility_cur == 0) {
485 if (!feasible) { /* reset upper bounds to infin. */
486 for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
487 nlopt_set_upper_bounds(dual_opt, dual_ub);
493 if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
494 else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
495 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
496 else if (feasible && *minf < stop->minf_max)
497 ret = NLOPT_MINF_MAX_REACHED;
498 if (ret != NLOPT_SUCCESS) goto done;
500 if (inner_done) break;
503 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
504 for (i = 0; i < m; ++i)
505 if (fcval_cur[i] > dd.gcval[i])
508 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
512 printf("CCSA inner iteration: rho -> %g\n", rho);
513 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
514 printf(" CCSA rhoc[%u] -> %g\n", i,rhoc[i]);
517 if (nlopt_stop_ftol(stop, fcur, fprev))
518 ret = NLOPT_FTOL_REACHED;
519 if (nlopt_stop_x(stop, xcur, xprev))
520 ret = NLOPT_XTOL_REACHED;
521 if (ret != NLOPT_SUCCESS) goto done;
523 /* update rho and sigma for iteration k+1 */
524 rho = MAX(0.1 * rho, CCSA_RHOMIN);
526 printf("CCSA outer iteration: rho -> %g\n", rho);
527 for (i = 0; i < m; ++i)
528 rhoc[i] = MAX(0.1 * rhoc[i], CCSA_RHOMIN);
529 for (i = 0; i < MIN(ccsa_verbose, m); ++i)
530 printf(" CCSA rhoc[%u] -> %g\n", i, rhoc[i]);
532 for (j = 0; j < n; ++j) {
533 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
534 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
536 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
537 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
538 /* use a smaller lower bound than Svanberg's
539 0.01*(ub-lb), which seems unnecessarily large */
540 sigma[j] = MAX(sigma[j], 1e-8*(ub[j]-lb[j]));
543 for (j = 0; j < MIN(ccsa_verbose, n); ++j)
544 printf(" CCSA sigma[%u] -> %g\n",
550 nlopt_destroy(pre_opt);
551 if (dd.scratch) free(dd.scratch);