mma/ccsa_quadratic.c

   1 /* Copyright (c) 2007-2011 Massachusetts Institute of Technology
   2  *
   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:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  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.
  21  */
  22
  23 /* In this file we implement Svanberg's CCSA algorithm with the
  24    simple linear approximation + quadratic penalty function.
  25
  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.]
  32
  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.
  38 */
  39
  40 #include <stdlib.h>
  41 #include <math.h>
  42 #include <string.h>
  43 #include <stdio.h>
  44
  45 #include "mma.h"
  46 #include "nlopt-util.h"
  47
  48 unsigned ccsa_verbose = 0; /* > 0 for verbose output */
  49
  50 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  51 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  52
  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
  59
  60 /***********************************************************************/
  61 /* function for CCSA's dual solution of the approximate problem */
  62
  63 typedef struct {
  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 */
  75 } dual_data;
  76
  77 static double sqr(double x) { return x * x; }
  78
  79 static double dual_func(unsigned m, const double *y, double *grad, void *d_)
  80 {
  81      dual_data *d = (dual_data *) d_;
  82      unsigned n = d->n;
  83      const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
  84           *dfdx = d->dfdx;
  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;
  90      unsigned i, j;
  91      double val;
  92
  93      d->count++;
  94
  95      val = d->gval = fval;
  96      d->wval = 0;
  97      for (i = 0; i < m; ++i)
  98           val += y[i] * (gcval[i] = fcval[i]);
  99
 100      for (j = 0; j < n; ++j) {
 101           double u, v, dx, sigma2, dx2, dx2sig;
 102
 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
 109              lb <= x+dx <= ub. */
 110
 111           if (sigma[j] == 0) { /* special case for lb[i] == ub[i] dims, dx=0 */
 112                xcur[j] = x[j];
 113                continue;
 114           }
 115
 116           u = rho;
 117           v = dfdx[j];
 118           for (i = 0; i < m; ++i) {
 119                u += rhoc[i] * y[i];
 120                v += dfcdx[i*n + j] * y[i];
 121           }
 122           dx = -(sigma2 = sqr(sigma[j])) * v/u;
 123
 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);
 127           xcur[j] = x[j] + dx;
 128           if (xcur[j] > ub[j]) xcur[j] = ub[j];
 129           else if (xcur[j] < lb[j]) xcur[j] = lb[j];
 130           dx = xcur[j] - x[j];
 131
 132           /* function value: */
 133           dx2 = dx * dx;
 134           val += v * dx + 0.5 * u * dx2 / sigma2;
 135
 136           /* update gval, wval, gcval (approximant functions) */
 137           d->gval += dfdx[j] * dx + rho * (dx2sig = 0.5*dx2/sigma2);
 138           d->wval += dx2sig;
 139           for (i = 0; i < m; ++i)
 140                gcval[i] += dfcdx[i*n+j] * dx + rhoc[i] * dx2sig;
 141      }
 142
 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];
 147      return -val;
 148 }
 149
 150 /***********************************************************************/
 151
 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,
 155                     const double *x0,
 156                     nlopt_precond pre, void *pre_data, double *scratch,
 157                     const double *x, double *grad)
 158 {
 159      double *dx = scratch, *Hdx = scratch + n;
 160      double val = f;
 161      unsigned j;
 162
 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;
 168      }
 169
 170      if (pre) {
 171           pre(n, x0, dx, Hdx, pre_data);
 172           for (j = 0; j < n; ++j)
 173                val += 0.5 * dx[j] * Hdx[j];
 174           if (grad)
 175                for (j = 0; j < n; ++j)
 176                     grad[j] = Hdx[j];
 177      }
 178
 179      return val;
 180 }
 181
 182 static double g0(unsigned n, const double *x, double *grad, void *d_)
 183 {
 184      dual_data *d = (dual_data *) d_;
 185      d->count++;
 186      return gfunc(n, d->fval, d->dfdx, d->rho, d->sigma,
 187                   d->x,
 188                   d->pre, d->pre_data, d->scratch,
 189                   x, grad);
 190 }
 191
 192
 193 static void gi(unsigned m, double *result,
 194                unsigned n, const double *x, double *grad, void *d_)
 195 {
 196      dual_data *d = (dual_data *) d_;
 197      unsigned i;
 198      for (i = 0; i < m; ++i)
 199           result[i] = gfunc(n, d->fcval[i], d->dfcdx + i*n, d->rhoc[i],
 200                             d->sigma,
 201                             d->x,
 202                             d->prec ? d->prec[i] : NULL,
 203                             d->prec_data ? d->prec_data[i] : NULL,
 204                             d->scratch,
 205                             x, grad);
 206 }
 207
 208
 209 /***********************************************************************/
 210
 211 nlopt_result ccsa_quadratic_minimize(
 212      unsigned n, nlopt_func f, void *f_data,
 213      unsigned m, nlopt_constraint *fc,
 214
 215      nlopt_precond pre,
 216
 217      const double *lb, const double *ub, /* bounds */
 218      double *x, /* in: initial guess, out: minimizer */
 219      double *minf,
 220      nlopt_stopping *stop,
 221      nlopt_opt dual_opt)
 222 {
 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;
 229      dual_data dd;
 230      int feasible;
 231      double infeasibility;
 232      unsigned mfc;
 233      unsigned no_precond;
 234      nlopt_opt pre_opt = NULL;
 235
 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;
 240      dfdx = sigma + n;
 241      dfdx_cur = dfdx + n;
 242      xcur = dfdx_cur + n;
 243      xprev = xcur + n;
 244      xprevprev = xprev + n;
 245      fcval = xprevprev + n;
 246      fcval_cur = fcval + m;
 247      rhoc = fcval_cur + m;
 248      gcval = rhoc + m;
 249      dual_lb = gcval + m;
 250      dual_ub = dual_lb + m;
 251      y = dual_ub + m;
 252      dfcdx = y + m;
 253      dfcdx_cur = dfcdx + m*n;
 254
 255      dd.n = n;
 256      dd.x = x;
 257      dd.lb = lb;
 258      dd.ub = ub;
 259      dd.sigma = sigma;
 260      dd.dfdx = dfdx;
 261      dd.dfcdx = dfcdx;
 262      dd.fcval = fcval;
 263      dd.rhoc = rhoc;
 264      dd.xcur = xcur;
 265      dd.gcval = gcval;
 266      dd.pre = pre; dd.pre_data = f_data;
 267      dd.prec = NULL; dd.prec_data = NULL;
 268      dd.scratch = NULL;
 269
 270      if (m) {
 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;
 275                goto done;
 276           }
 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;
 282                }
 283           }
 284      }
 285
 286      no_precond = pre == NULL;
 287      if (dd.prec)
 288           for (i = 0; i < m; ++i)
 289                no_precond = no_precond && dd.prec[i] == NULL;
 290
 291      if (!no_precond) {
 292           dd.scratch = (double*) malloc(sizeof(double) * (2*n));
 293           if (!dd.scratch) {
 294                free(sigma);
 295                return NLOPT_OUT_OF_MEMORY;
 296           }
 297           pre_lb = dual_lb;
 298           pre_ub = dual_ub;
 299
 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;
 310      }
 311
 312      for (j = 0; j < n; ++j) {
 313           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 314                sigma[j] = 1.0; /* arbitrary default */
 315           else
 316                sigma[j] = 0.5 * (ub[j] - lb[j]);
 317      }
 318      rho = 1.0;
 319      for (i = 0; i < m; ++i) {
 320           rhoc[i] = 1.0;
 321           dual_lb[i] = y[i] = 0.0;
 322           dual_ub[i] = HUGE_VAL;
 323      }
 324
 325      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 326      stop->nevals++;
 327      memcpy(xcur, x, sizeof(double) * n);
 328      if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 329
 330      feasible = 1; infeasibility = 0;
 331      for (i = ifc = 0; ifc < mfc; ++ifc) {
 332           nlopt_eval_constraint(fcval + i, dfcdx + i*n,
 333                                 fc + ifc, n, x);
 334           i += fc[ifc].m;
 335           if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 336      }
 337      for (i = 0; i < m; ++i) {
 338           feasible = feasible && fcval[i] <= 0;
 339           if (fcval[i] > infeasibility) infeasibility = fcval[i];
 340      }
 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.
 348
 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. */
 352      if (!feasible)
 353           for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
 354
 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);
 361
 362      while (1) { /* outer iterations */
 363           double fprev = fcur;
 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);
 372
 373           while (1) { /* inner iterations */
 374                double min_dual, infeasibility_cur;
 375                int feasible_cur, inner_done;
 376                unsigned save_verbose;
 377                nlopt_result reti;
 378
 379                if (no_precond) {
 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,
 385                                                   0,
 386                                                   stop->maxtime
 387                                                   - (nlopt_seconds()
 388                                                      - stop->start));
 389                     ccsa_verbose = save_verbose;
 390                     if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 391                          ret = reti;
 392                          goto done;
 393                     }
 394
 395                     dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 396                }
 397                else {
 398                     double pre_min;
 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]);
 402                          xcur[j] = x[j];
 403                     }
 404                     nlopt_set_lower_bounds(pre_opt, pre_lb);
 405                     nlopt_set_upper_bounds(pre_opt, pre_ub);
 406
 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,
 411                                                   0, stop->maxtime
 412                                                   - (nlopt_seconds()
 413                                                      - stop->start));
 414                     ccsa_verbose = save_verbose;
 415                     if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 416                          ret = reti;
 417                          goto done;
 418                     }
 419
 420                     /* evaluate final xcur etc */
 421                     dd.gval = g0(n, xcur, NULL, &dd);
 422                     gi(m, dd.gcval, n, xcur, NULL, &dd);
 423                }
 424
 425                if (ccsa_verbose) {
 426                     printf("CCSA dual converged in %d iters to g=%g:\n",
 427                            dd.count, dd.gval);
 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]);
 431                }
 432
 433                fcur = f(n, xcur, dfdx_cur, f_data);
 434                stop->nevals++;
 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,
 441                                           fc + ifc, n, xcur);
 442                     i += fc[ifc].m;
 443                     if (nlopt_stop_forced(stop)) {
 444                          ret = NLOPT_FORCED_STOP; goto done; }
 445                }
 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];
 455                     }
 456                }
 457
 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);
 468
 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);
 478                          }
 479                          feasible = 1;
 480                     }
 481
 482                }
 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;
 489
 490                if (inner_done) break;
 491
 492                if (fcur > dd.gval)
 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])
 496                          rhoc[i] =
 497                               MIN(10*rhoc[i],
 498                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 499                                          / dd.wval));
 500
 501                if (ccsa_verbose)
 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]);
 505           }
 506
 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;
 512
 513           /* update rho and sigma for iteration k+1 */
 514           rho = MAX(0.1 * rho, CCSA_RHOMIN);
 515           if (ccsa_verbose)
 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]);
 521           if (k > 1) {
 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);
 525                     sigma[j] *= gam;
 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]));
 531                     }
 532                }
 533                for (j = 0; j < MIN(ccsa_verbose, n); ++j)
 534                     printf("                 CCSA sigma[%d] -> %g\n",
 535                            j, sigma[j]);
 536           }
 537      }
 538
 539  done:
 540      nlopt_destroy(pre_opt);
 541      if (dd.scratch) free(dd.scratch);
 542      if (m) {
 543           free(dd.prec_data);
 544           free(dd.prec);
 545      }
 546      free(sigma);
 547      return ret;
 548 }