mma/mma.c

   1 /* Copyright (c) 2007-2008 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 #include <stdlib.h>
  24 #include <math.h>
  25 #include <string.h>
  26 #include <stdio.h>
  27
  28 #include "mma.h"
  29 #include "config.h"
  30
  31 int mma_verbose = 0; /* > 0 for verbose output */
  32
  33 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  34 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  35
  36 #ifndef HAVE_ISNAN
  37 static int my_isnan(double x) { return x != x; }
  38 #  define isnan my_isnan
  39 #endif
  40
  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
  47
  48 /***********************************************************************/
  49 /* function for MMA's dual solution of the approximate problem */
  50
  51 typedef struct {
  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) */
  60 } dual_data;
  61
  62 static double sqr(double x) { return x * x; }
  63
  64 static double dual_func(int m, const double *y, double *grad, void *d_)
  65 {
  66      dual_data *d = (dual_data *) d_;
  67      int n = d->n;
  68      const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
  69           *dfdx = d->dfdx;
  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;
  75      int i, j;
  76      double val;
  77
  78      d->count++;
  79
  80      val = d->gval = fval;
  81      d->wval = 0;
  82      for (i = 0; i < m; ++i)
  83           val += y[i] * (gcval[i] = isnan(fcval[i]) ? 0 : fcval[i]);
  84
  85      for (j = 0; j < n; ++j) {
  86           double u, v, dx, denominv, c, sigma2, dx2;
  87
  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
  95              is given by:
  96                      (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2))
  97              (which goes to zero as u -> 0). */
  98
  99           u = dfdx[j];
 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];
 104           }
 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);
 109           }
 110           else
 111                dx = (v/u)*sigma2 * (-1 + sqrt(fabs(1 - sqr(u/(v*sigma[j])))));
 112           xcur[j] = x[j] + dx;
 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];
 117           dx = xcur[j] - x[j];
 118
 119           /* function value: */
 120           dx2 = dx * dx;
 121           denominv = 1.0 / (sigma2 - dx2);
 122           val += (u * dx + v * dx2) * denominv;
 123
 124           /* update gval, wval, gcval (approximant functions) */
 125           c = sigma2 * dx;
 126           d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2)
 127                * denominv;
 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)
 132                     * denominv;
 133      }
 134
 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];
 139      return -val;
 140 }
 141
 142 /***********************************************************************/
 143
 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. */
 147
 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 */
 153                           double *minf,
 154                           nlopt_stopping *stop,
 155                           nlopt_algorithm dual_alg,
 156                           double dual_tolrel, int dual_maxeval)
 157 {
 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;
 162      int i, j, k = 0;
 163      char *fc_data = (char *) fc_data_;
 164      dual_data dd;
 165      int feasible;
 166
 167      sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
 168      if (!sigma) return NLOPT_OUT_OF_MEMORY;
 169      dfdx = sigma + n;
 170      dfdx_cur = dfdx + n;
 171      xcur = dfdx_cur + n;
 172      xprev = xcur + n;
 173      xprevprev = xprev + n;
 174      fcval = xprevprev + n;
 175      fcval_cur = fcval + m;
 176      rhoc = fcval_cur + m;
 177      gcval = rhoc + m;
 178      dual_lb = gcval + m;
 179      dual_ub = dual_lb + m;
 180      y = dual_ub + m;
 181      dfcdx = y + m;
 182      dfcdx_cur = dfcdx + m*n;
 183
 184      dd.n = n;
 185      dd.x = x;
 186      dd.lb = lb;
 187      dd.ub = ub;
 188      dd.sigma = sigma;
 189      dd.dfdx = dfdx;
 190      dd.dfcdx = dfcdx;
 191      dd.fcval = fcval;
 192      dd.rhoc = rhoc;
 193      dd.xcur = xcur;
 194      dd.gcval = gcval;
 195
 196      for (j = 0; j < n; ++j) {
 197           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 198                sigma[j] = 1.0; /* arbitrary default */
 199           else
 200                sigma[j] = 0.5 * (ub[j] - lb[j]);
 201      }
 202      rho = 1.0;
 203      for (i = 0; i < m; ++i) {
 204           rhoc[i] = 1.0;
 205           dual_lb[i] = y[i] = 0.0;
 206           dual_ub[i] = HUGE_VAL;
 207      }
 208
 209      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 210      stop->nevals++;
 211      memcpy(xcur, x, sizeof(double) * n);
 212
 213      feasible = 1;
 214      for (i = 0; i < m; ++i) {
 215           fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
 216           feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i]));
 217      }
 218      /* For non-feasible initial points, set a finite (large)
 219         upper-bound on the dual variables.  What this means is that,
 220         if no feasible solution is found from the dual problem, it
 221         will minimize the dual objective with the unfeasible
 222         constraint weighted by 1e40 -- basically, minimizing the
 223         unfeasible constraint until it becomes feasible or until we at
 224         least obtain a step towards a feasible point.
 225
 226         Svanberg suggested a different approach in his 1987 paper, basically
 227         introducing additional penalty variables for unfeasible constraints,
 228         but this is easier to implement and at least as efficient. */
 229      if (!feasible)
 230           for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
 231
 232      while (1) { /* outer iterations */
 233           double fprev = fcur;
 234           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 235           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 236           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 237           if (ret != NLOPT_SUCCESS) goto done;
 238           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 239           memcpy(xprev, xcur, sizeof(double) * n);
 240
 241           while (1) { /* inner iterations */
 242                double min_dual;
 243                int feasible_cur, inner_done, save_verbose;
 244                int new_infeasible_constraint;
 245                nlopt_result reti;
 246
 247                /* solve dual problem */
 248                dd.rho = rho; dd.count = 0;
 249           dual_solution:
 250                save_verbose = mma_verbose;
 251                mma_verbose = 0;
 252                reti = nlopt_minimize(
 253                     dual_alg, m, dual_func, &dd,
 254                     dual_lb, dual_ub, y, &min_dual,
 255                     -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
 256                     stop->maxtime - (nlopt_seconds() - stop->start));
 257                mma_verbose = save_verbose;
 258                if (reti == NLOPT_FAILURE && dual_alg != NLOPT_LD_MMA) {
 259                     /* LBFGS etc. converge quickly but are sometimes
 260                        very finicky if there are any rounding errors in
 261                        the gradient, etcetera; if it fails, try again
 262                        with MMA called recursively for the dual */
 263                     dual_alg = NLOPT_LD_MMA;
 264                     if (mma_verbose)
 265                          printf("MMA: switching to recursive MMA for dual\n");
 266                     goto dual_solution;
 267                }
 268                if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 269                     ret = reti;
 270                     goto done;
 271                }
 272
 273                dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 274                if (mma_verbose) {
 275                     printf("MMA dual converged in %d iterations to g=%g:\n",
 276                            dd.count, dd.gval);
 277                     for (i = 0; i < MIN(mma_verbose, m); ++i)
 278                          printf("    MMA y[%d]=%g, gc[%d]=%g\n",
 279                                 i, y[i], i, dd.gcval[i]);
 280                }
 281
 282                fcur = f(n, xcur, dfdx_cur, f_data);
 283                stop->nevals++;
 284                feasible_cur = 1;
 285                new_infeasible_constraint = 0;
 286                inner_done = dd.gval >= fcur;
 287                for (i = 0; i < m; ++i) {
 288                     fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
 289                                       fc_data + fc_datum_size * i);
 290                     if (!isnan(fcval_cur[i])) {
 291                          feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
 292                          if (!isnan(fcval[i]))
 293                               inner_done = inner_done &&
 294                                    (dd.gcval[i] >= fcval_cur[i]);
 295                          else if (fcval_cur[i] > 0)
 296                               new_infeasible_constraint = 1;
 297                     }
 298                }
 299
 300                if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
 301                     || (!feasible && feasible_cur)) {
 302                     if (mma_verbose && !feasible_cur)
 303                          printf("MMA - using infeasible point?\n");
 304                     dd.fval = *minf = fcur;
 305                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 306                     memcpy(x, xcur, sizeof(double)*n);
 307                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 308                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 309
 310                     /* once we have reached a feasible solution, the
 311                        algorithm should never make the solution infeasible
 312                        again (if inner_done), although the constraints may
 313                        be violated slightly by rounding errors etc. so we
 314                        must be a little careful about checking feasibility */
 315                     if (feasible_cur) {
 316                          if (!feasible) /* reset upper bounds to infinity */
 317                               for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
 318                          feasible = 1;
 319                     }
 320                     else if (new_infeasible_constraint) feasible = 0;
 321
 322                }
 323                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 324                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 325                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 326                if (ret != NLOPT_SUCCESS) goto done;
 327
 328                if (inner_done) break;
 329
 330                if (fcur > dd.gval)
 331                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 332                for (i = 0; i < m; ++i)
 333                     if (!isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i])
 334                          rhoc[i] =
 335                               MIN(10*rhoc[i],
 336                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 337                                          / dd.wval));
 338
 339                if (mma_verbose)
 340                     printf("MMA inner iteration: rho -> %g\n", rho);
 341                for (i = 0; i < MIN(mma_verbose, m); ++i)
 342                     printf("                 MMA rhoc[%d] -> %g\n", i,rhoc[i]);
 343           }
 344
 345           if (nlopt_stop_ftol(stop, fcur, fprev))
 346                ret = NLOPT_FTOL_REACHED;
 347           if (nlopt_stop_x(stop, xcur, xprev))
 348                ret = NLOPT_XTOL_REACHED;
 349           if (ret != NLOPT_SUCCESS) goto done;
 350
 351           /* update rho and sigma for iteration k+1 */
 352           rho = MAX(0.1 * rho, MMA_RHOMIN);
 353           if (mma_verbose)
 354                printf("MMA outer iteration: rho -> %g\n", rho);
 355           for (i = 0; i < m; ++i)
 356                rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
 357           for (i = 0; i < MIN(mma_verbose, m); ++i)
 358                printf("                 MMA rhoc[%d] -> %g\n", i, rhoc[i]);
 359           if (k > 1) {
 360                for (j = 0; j < n; ++j) {
 361                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 362                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 363                     sigma[j] *= gam;
 364                     if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
 365                          sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 366                          sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 367                     }
 368                }
 369                for (j = 0; j < MIN(mma_verbose, n); ++j)
 370                     printf("                 MMA sigma[%d] -> %g\n",
 371                            j, sigma[j]);
 372           }
 373      }
 374
 375  done:
 376      free(sigma);
 377      return ret;
 378 }