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      if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */
 219
 220      while (1) { /* outer iterations */
 221           double fprev = fcur;
 222           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 223           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 224           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 225           if (ret != NLOPT_SUCCESS) goto done;
 226           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 227           memcpy(xprev, xcur, sizeof(double) * n);
 228
 229           while (1) { /* inner iterations */
 230                double min_dual;
 231                int feasible_cur, inner_done, save_verbose;
 232                int new_infeasible_constraint;
 233                nlopt_result reti;
 234
 235                /* solve dual problem */
 236                dd.rho = rho; dd.count = 0;
 237           dual_solution:
 238                save_verbose = mma_verbose;
 239                mma_verbose = 0;
 240                reti = nlopt_minimize(
 241                     dual_alg, m, dual_func, &dd,
 242                     dual_lb, dual_ub, y, &min_dual,
 243                     -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
 244                     stop->maxtime - (nlopt_seconds() - stop->start));
 245                mma_verbose = save_verbose;
 246                if (reti == NLOPT_FAILURE && dual_alg != NLOPT_LD_MMA) {
 247                     /* LBFGS etc. converge quickly but are sometimes
 248                        very finicky if there are any rounding errors in
 249                        the gradient, etcetera; if it fails, try again
 250                        with MMA called recursively for the dual */
 251                     dual_alg = NLOPT_LD_MMA;
 252                     if (mma_verbose)
 253                          printf("MMA: switching to recursive MMA for dual\n");
 254                     goto dual_solution;
 255                }
 256                if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 257                     ret = reti;
 258                     goto done;
 259                }
 260
 261                dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 262                if (mma_verbose) {
 263                     printf("MMA dual converged in %d iterations to g=%g:\n",
 264                            dd.count, dd.gval);
 265                     for (i = 0; i < MIN(mma_verbose, m); ++i)
 266                          printf("    MMA y[%d]=%g, gc[%d]=%g\n",
 267                                 i, y[i], i, dd.gcval[i]);
 268                }
 269
 270                fcur = f(n, xcur, dfdx_cur, f_data);
 271                stop->nevals++;
 272                feasible_cur = 1;
 273                new_infeasible_constraint = 0;
 274                inner_done = dd.gval >= fcur;
 275                for (i = 0; i < m; ++i) {
 276                     fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
 277                                       fc_data + fc_datum_size * i);
 278                     if (!isnan(fcval_cur[i])) {
 279                          feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
 280                          if (!isnan(fcval[i]))
 281                               inner_done = inner_done &&
 282                                    (dd.gcval[i] >= fcval_cur[i]);
 283                          else if (fcval_cur[i] > 0)
 284                               new_infeasible_constraint = 1;
 285                     }
 286                }
 287
 288                if (fcur < *minf && (inner_done || feasible_cur || !feasible)) {
 289                     if (mma_verbose && !feasible_cur)
 290                          printf("MMA - using infeasible point?\n");
 291                     dd.fval = *minf = fcur;
 292                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 293                     memcpy(x, xcur, sizeof(double)*n);
 294                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 295                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 296
 297                     /* once we have reached a feasible solution, the
 298                        algorithm should never make the solution infeasible
 299                        again (if inner_done), although the constraints may
 300                        be violated slightly by rounding errors etc. so we
 301                        must be a little careful about checking feasibility */
 302                     if (feasible_cur) feasible = 1;
 303                     else if (new_infeasible_constraint) feasible = 0;
 304
 305                }
 306                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 307                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 308                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 309                if (ret != NLOPT_SUCCESS) goto done;
 310
 311                if (inner_done) break;
 312
 313                if (fcur > dd.gval)
 314                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 315                for (i = 0; i < m; ++i)
 316                     if (!isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i])
 317                          rhoc[i] =
 318                               MIN(10*rhoc[i],
 319                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 320                                          / dd.wval));
 321
 322                if (mma_verbose)
 323                     printf("MMA inner iteration: rho -> %g\n", rho);
 324                for (i = 0; i < MIN(mma_verbose, m); ++i)
 325                     printf("                 MMA rhoc[%d] -> %g\n", i,rhoc[i]);
 326           }
 327
 328           if (nlopt_stop_ftol(stop, fcur, fprev))
 329                ret = NLOPT_FTOL_REACHED;
 330           if (nlopt_stop_x(stop, xcur, xprev))
 331                ret = NLOPT_XTOL_REACHED;
 332           if (ret != NLOPT_SUCCESS) goto done;
 333
 334           /* update rho and sigma for iteration k+1 */
 335           rho = MAX(0.1 * rho, MMA_RHOMIN);
 336           if (mma_verbose)
 337                printf("MMA outer iteration: rho -> %g\n", rho);
 338           for (i = 0; i < m; ++i)
 339                rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
 340           for (i = 0; i < MIN(mma_verbose, m); ++i)
 341                printf("                 MMA rhoc[%d] -> %g\n", i, rhoc[i]);
 342           if (k > 1) {
 343                for (j = 0; j < n; ++j) {
 344                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 345                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 346                     sigma[j] *= gam;
 347                     if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
 348                          sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 349                          sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 350                     }
 351                }
 352                for (j = 0; j < MIN(mma_verbose, n); ++j)
 353                     printf("                 MMA sigma[%d] -> %g\n",
 354                            j, sigma[j]);
 355           }
 356      }
 357
 358  done:
 359      free(sigma);
 360      return ret;
 361 }