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      double infeasibility;
 167
 168      sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
 169      if (!sigma) return NLOPT_OUT_OF_MEMORY;
 170      dfdx = sigma + n;
 171      dfdx_cur = dfdx + n;
 172      xcur = dfdx_cur + n;
 173      xprev = xcur + n;
 174      xprevprev = xprev + n;
 175      fcval = xprevprev + n;
 176      fcval_cur = fcval + m;
 177      rhoc = fcval_cur + m;
 178      gcval = rhoc + m;
 179      dual_lb = gcval + m;
 180      dual_ub = dual_lb + m;
 181      y = dual_ub + m;
 182      dfcdx = y + m;
 183      dfcdx_cur = dfcdx + m*n;
 184
 185      dd.n = n;
 186      dd.x = x;
 187      dd.lb = lb;
 188      dd.ub = ub;
 189      dd.sigma = sigma;
 190      dd.dfdx = dfdx;
 191      dd.dfcdx = dfcdx;
 192      dd.fcval = fcval;
 193      dd.rhoc = rhoc;
 194      dd.xcur = xcur;
 195      dd.gcval = gcval;
 196
 197      for (j = 0; j < n; ++j) {
 198           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 199                sigma[j] = 1.0; /* arbitrary default */
 200           else
 201                sigma[j] = 0.5 * (ub[j] - lb[j]);
 202      }
 203      rho = 1.0;
 204      for (i = 0; i < m; ++i) {
 205           rhoc[i] = 1.0;
 206           dual_lb[i] = y[i] = 0.0;
 207           dual_ub[i] = HUGE_VAL;
 208      }
 209
 210      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 211      stop->nevals++;
 212      memcpy(xcur, x, sizeof(double) * n);
 213
 214      feasible = 1; infeasibility = 0;
 215      for (i = 0; i < m; ++i) {
 216           fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
 217           feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i]));
 218           if (fcval[i] > infeasibility) infeasibility = fcval[i];
 219      }
 220      /* For non-feasible initial points, set a finite (large)
 221         upper-bound on the dual variables.  What this means is that,
 222         if no feasible solution is found from the dual problem, it
 223         will minimize the dual objective with the unfeasible
 224         constraint weighted by 1e40 -- basically, minimizing the
 225         unfeasible constraint until it becomes feasible or until we at
 226         least obtain a step towards a feasible point.
 227
 228         Svanberg suggested a different approach in his 1987 paper, basically
 229         introducing additional penalty variables for unfeasible constraints,
 230         but this is easier to implement and at least as efficient. */
 231      if (!feasible)
 232           for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
 233
 234      while (1) { /* outer iterations */
 235           double fprev = fcur;
 236           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 237           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 238           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 239           if (ret != NLOPT_SUCCESS) goto done;
 240           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 241           memcpy(xprev, xcur, sizeof(double) * n);
 242
 243           while (1) { /* inner iterations */
 244                double min_dual, infeasibility_cur;
 245                int feasible_cur, inner_done, save_verbose;
 246                int new_infeasible_constraint;
 247                nlopt_result reti;
 248
 249                /* solve dual problem */
 250                dd.rho = rho; dd.count = 0;
 251           dual_solution:
 252                save_verbose = mma_verbose;
 253                mma_verbose = 0;
 254                reti = nlopt_minimize(
 255                     dual_alg, m, dual_func, &dd,
 256                     dual_lb, dual_ub, y, &min_dual,
 257                     -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
 258                     stop->maxtime - (nlopt_seconds() - stop->start));
 259                mma_verbose = save_verbose;
 260                if (reti == NLOPT_FAILURE && dual_alg != NLOPT_LD_MMA) {
 261                     /* LBFGS etc. converge quickly but are sometimes
 262                        very finicky if there are any rounding errors in
 263                        the gradient, etcetera; if it fails, try again
 264                        with MMA called recursively for the dual */
 265                     dual_alg = NLOPT_LD_MMA;
 266                     if (mma_verbose)
 267                          printf("MMA: switching to recursive MMA for dual\n");
 268                     goto dual_solution;
 269                }
 270                if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 271                     ret = reti;
 272                     goto done;
 273                }
 274
 275                dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 276                if (mma_verbose) {
 277                     printf("MMA dual converged in %d iterations to g=%g:\n",
 278                            dd.count, dd.gval);
 279                     for (i = 0; i < MIN(mma_verbose, m); ++i)
 280                          printf("    MMA y[%d]=%g, gc[%d]=%g\n",
 281                                 i, y[i], i, dd.gcval[i]);
 282                }
 283
 284                fcur = f(n, xcur, dfdx_cur, f_data);
 285                stop->nevals++;
 286                feasible_cur = 1; infeasibility_cur = 0;
 287                new_infeasible_constraint = 0;
 288                inner_done = dd.gval >= fcur;
 289                for (i = 0; i < m; ++i) {
 290                     fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
 291                                       fc_data + fc_datum_size * i);
 292                     if (!isnan(fcval_cur[i])) {
 293                          feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
 294                          if (!isnan(fcval[i]))
 295                               inner_done = inner_done &&
 296                                    (dd.gcval[i] >= fcval_cur[i]);
 297                          else if (fcval_cur[i] > 0)
 298                               new_infeasible_constraint = 1;
 299                          if (fcval_cur[i] > infeasibility_cur)
 300                               infeasibility_cur = fcval_cur[i];
 301                     }
 302                }
 303
 304                if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
 305                     || (!feasible && infeasibility_cur < infeasibility)) {
 306                     if (mma_verbose && !feasible_cur)
 307                          printf("MMA - using infeasible point?\n");
 308                     dd.fval = *minf = fcur;
 309                     infeasibility = infeasibility_cur;
 310                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 311                     memcpy(x, xcur, sizeof(double)*n);
 312                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 313                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 314
 315                     /* once we have reached a feasible solution, the
 316                        algorithm should never make the solution infeasible
 317                        again (if inner_done), although the constraints may
 318                        be violated slightly by rounding errors etc. so we
 319                        must be a little careful about checking feasibility */
 320                     if (feasible_cur) {
 321                          if (!feasible) /* reset upper bounds to infinity */
 322                               for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
 323                          feasible = 1;
 324                     }
 325                     else if (new_infeasible_constraint) feasible = 0;
 326
 327                }
 328                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 329                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 330                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 331                if (ret != NLOPT_SUCCESS) goto done;
 332
 333                if (inner_done) break;
 334
 335                if (fcur > dd.gval)
 336                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 337                for (i = 0; i < m; ++i)
 338                     if (!isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i])
 339                          rhoc[i] =
 340                               MIN(10*rhoc[i],
 341                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 342                                          / dd.wval));
 343
 344                if (mma_verbose)
 345                     printf("MMA inner iteration: rho -> %g\n", rho);
 346                for (i = 0; i < MIN(mma_verbose, m); ++i)
 347                     printf("                 MMA rhoc[%d] -> %g\n", i,rhoc[i]);
 348           }
 349
 350           if (nlopt_stop_ftol(stop, fcur, fprev))
 351                ret = NLOPT_FTOL_REACHED;
 352           if (nlopt_stop_x(stop, xcur, xprev))
 353                ret = NLOPT_XTOL_REACHED;
 354           if (ret != NLOPT_SUCCESS) goto done;
 355
 356           /* update rho and sigma for iteration k+1 */
 357           rho = MAX(0.1 * rho, MMA_RHOMIN);
 358           if (mma_verbose)
 359                printf("MMA outer iteration: rho -> %g\n", rho);
 360           for (i = 0; i < m; ++i)
 361                rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
 362           for (i = 0; i < MIN(mma_verbose, m); ++i)
 363                printf("                 MMA rhoc[%d] -> %g\n", i, rhoc[i]);
 364           if (k > 1) {
 365                for (j = 0; j < n; ++j) {
 366                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 367                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 368                     sigma[j] *= gam;
 369                     if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
 370                          sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 371                          sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 372                     }
 373                }
 374                for (j = 0; j < MIN(mma_verbose, n); ++j)
 375                     printf("                 MMA sigma[%d] -> %g\n",
 376                            j, sigma[j]);
 377           }
 378      }
 379
 380  done:
 381      free(sigma);
 382      return ret;
 383 }