mma/mma.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include <string.h>
   4 #include <stdio.h>
   5
   6 #include "mma.h"
   7
   8 int mma_verbose = 0; /* > 0 for verbose output */
   9
  10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  12
  13 /* magic minimum value for rho in MMA ... the 2002 paper says it should
  14    be a "fixed, strictly positive `small' number, e.g. 1e-5"
  15    ... grrr, I hate these magic numbers, which seem like they
  16    should depend on the objective function in some way ... in particular,
  17    note that rho is dimensionful (= dimensions of objective function) */
  18 #define MMA_RHOMIN 1e-5
  19
  20 /***********************************************************************/
  21 /* function for MMA's dual solution of the approximate problem */
  22
  23 typedef struct {
  24      int n; /* must be set on input to dimension of x */
  25      const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
  26      const double *dfcdx; /* m-by-n array of fc gradients */
  27      double fval, rho; /* must be set on input */
  28      const double *fcval, *rhoc; /* arrays of length m */
  29      double *xcur; /* array of length n, output each time */
  30      double gval, wval, *gcval; /* output each time (array length m) */
  31 } dual_data;
  32
  33 static double sqr(double x) { return x * x; }
  34
  35 static double dual_func(int m, const double *y, double *grad, void *d_)
  36 {
  37      dual_data *d = (dual_data *) d_;
  38      int n = d->n;
  39      const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
  40           *dfdx = d->dfdx;
  41      const double *dfcdx = d->dfcdx;
  42      double rho = d->rho, fval = d->fval;
  43      const double *rhoc = d->rhoc, *fcval = d->fcval;
  44      double *xcur = d->xcur;
  45      double *gcval = d->gcval;
  46      int i, j;
  47      double val;
  48
  49      val = d->gval = fval;
  50      d->wval = 0;
  51      for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]);
  52
  53      for (j = 0; j < n; ++j) {
  54           double u, v, dx, denominv, c, sigma2, dx2;
  55
  56           /* first, compute xcur[j] for y.  Because this objective is
  57              separable, we can minimize over x analytically, and the minimum
  58              dx is given by the solution of a quadratic equation:
  59                      u dx^2 + 2 v sigma^2 dx + u sigma^2 = 0
  60              where u and v are defined by the sums below.  Because of
  61              the definitions, it is guaranteed that |u/v| <= sigma,
  62              and it follows that the only dx solution with |dx| <= sigma
  63              is given by:
  64                      (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2))
  65              (which goes to zero as u -> 0). */
  66
  67           u = dfdx[j];
  68           v = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
  69           for (i = 0; i < m; ++i) {
  70                u += dfcdx[i*n + j] * y[i];
  71                v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
  72           }
  73           u *= (sigma2 = sqr(sigma[j]));
  74           dx = u==0 ? 0 : (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j]))));
  75           xcur[j] = x[j] + dx;
  76           if (xcur[j] > ub[j]) xcur[j] = ub[j];
  77           else if (xcur[j] < lb[j]) xcur[j] = lb[j];
  78           if (xcur[j] > x[j]+0.9*sigma[j]) xcur[j] = x[j]+0.9*sigma[j];
  79           else if (xcur[j] < x[j]-0.9*sigma[j]) xcur[j] = x[j]-0.9*sigma[j];
  80           dx = xcur[j] - x[j];
  81
  82           /* function value: */
  83           dx2 = dx * dx;
  84           denominv = 1.0 / (sigma2 - dx2);
  85           val += (u * dx + v * dx2) * denominv;
  86
  87           /* update gval, wval, gcval (approximant functions) */
  88           c = sigma2 * dx;
  89           d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2)
  90                * denominv;
  91           d->wval += 0.5 * dx2 * denominv;
  92           for (i = 0; i < m; ++i)
  93                gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j]
  94                                                 + 0.5*rhoc[j]) * dx2)
  95                     * denominv;
  96      }
  97
  98      /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
  99         we only need the partial derivative with respect to y, and
 100         we negate because we are maximizing: */
 101      if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
 102      return -val;
 103 }
 104
 105 /***********************************************************************/
 106
 107 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
 108                           int m, nlopt_func fc,
 109                           void *fc_data_, ptrdiff_t fc_datum_size,
 110                           const double *lb, const double *ub, /* bounds */
 111                           double *x, /* in: initial guess, out: minimizer */
 112                           double *minf,
 113                           nlopt_stopping *stop,
 114                           nlopt_algorithm dual_alg,
 115                           double dual_tolrel, int dual_maxeval)
 116 {
 117      nlopt_result ret = NLOPT_SUCCESS;
 118      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
 119      double *dfcdx, *dfcdx_cur;
 120      double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
 121      int i, j, k = 0;
 122      char *fc_data = (char *) fc_data_;
 123      dual_data dd;
 124      int feasible;
 125
 126      sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
 127      if (!sigma) return NLOPT_OUT_OF_MEMORY;
 128      dfdx = sigma + n;
 129      dfdx_cur = dfdx + n;
 130      xcur = dfdx_cur + n;
 131      xprev = xcur + n;
 132      xprevprev = xprev + n;
 133      fcval = xprevprev + n;
 134      fcval_cur = fcval + m;
 135      rhoc = fcval_cur + m;
 136      gcval = rhoc + m;
 137      dual_lb = gcval + m;
 138      dual_ub = dual_lb + m;
 139      y = dual_ub + m;
 140      dfcdx = y + m;
 141      dfcdx_cur = dfcdx + m*n;
 142
 143      dd.n = n;
 144      dd.x = x;
 145      dd.lb = lb;
 146      dd.ub = ub;
 147      dd.sigma = sigma;
 148      dd.dfdx = dfdx;
 149      dd.dfcdx = dfcdx;
 150      dd.fcval = fcval;
 151      dd.rhoc = rhoc;
 152      dd.xcur = xcur;
 153      dd.gcval = gcval;
 154
 155      for (j = 0; j < n; ++j) {
 156           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 157                sigma[j] = 1.0; /* arbitrary default */
 158           else
 159                sigma[j] = 0.5 * (ub[j] - lb[j]);
 160      }
 161      rho = 1.0;
 162      for (i = 0; i < m; ++i) {
 163           rhoc[i] = 1.0;
 164           dual_lb[i] = y[i] = 0.0;
 165           dual_ub[i] = HUGE_VAL;
 166      }
 167
 168      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 169      stop->nevals++;
 170      memcpy(xcur, x, sizeof(double) * n);
 171
 172      feasible = 1;
 173      for (i = 0; i < m; ++i) {
 174           fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
 175           feasible = feasible && (fcval[i] <= 0);
 176      }
 177      if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */
 178
 179      while (1) { /* outer iterations */
 180           double fprev = fcur;
 181           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 182           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 183           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 184           if (ret != NLOPT_SUCCESS) goto done;
 185           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 186           memcpy(xprev, xcur, sizeof(double) * n);
 187
 188           while (1) { /* inner iterations */
 189                double min_dual;
 190                int feasible_cur, inner_done;
 191                nlopt_result reti;
 192
 193                /* solve dual problem */
 194                dd.rho = rho;
 195                reti = nlopt_minimize(
 196                     dual_alg, m, dual_func, &dd,
 197                     dual_lb, dual_ub, y, &min_dual,
 198                     -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
 199                     stop->maxtime - (nlopt_seconds() - stop->start));
 200                if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 201                     ret = reti;
 202                     goto done;
 203                }
 204
 205                dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 206
 207                fcur = f(n, xcur, dfdx_cur, f_data);
 208                stop->nevals++;
 209                inner_done = dd.gval >= fcur;
 210                feasible_cur = 1;
 211                for (i = 0; i < m; ++i) {
 212                     fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
 213                                       fc_data + fc_datum_size * i);
 214                     feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
 215                     inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
 216                }
 217
 218                if (fcur < *minf && (feasible_cur || !feasible)) {
 219                     feasible = feasible_cur;
 220                     dd.fval = *minf = fcur;
 221                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 222                     memcpy(x, xcur, sizeof(double)*n);
 223                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 224                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 225                }
 226                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 227                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 228                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 229                if (ret != NLOPT_SUCCESS) goto done;
 230
 231                if (inner_done) break;
 232
 233                if (fcur > dd.gval)
 234                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 235                for (i = 0; i < m; ++i)
 236                     if (fcval_cur[i] > dd.gcval[i])
 237                          rhoc[i] =
 238                               MIN(10*rhoc[i],
 239                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 240                                          / dd.wval));
 241
 242                if (mma_verbose)
 243                     printf("MMA inner iteration: rho -> %g\n", rho);
 244                for (i = 0; i < mma_verbose; ++i)
 245                     printf("                     rhoc[%d] -> %g\n", i,rhoc[i]);
 246           }
 247
 248           if (nlopt_stop_ftol(stop, fcur, fprev))
 249                ret = NLOPT_FTOL_REACHED;
 250           if (nlopt_stop_x(stop, xcur, xprev))
 251                ret = NLOPT_XTOL_REACHED;
 252           if (ret != NLOPT_SUCCESS) goto done;
 253
 254           /* update rho and sigma for iteration k+1 */
 255           rho = MAX(0.1 * rho, MMA_RHOMIN);
 256           if (mma_verbose)
 257                printf("MMA outer iteration: rho -> %g\n", rho);
 258           for (i = 0; i < m; ++i)
 259                rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
 260           for (i = 0; i < mma_verbose; ++i)
 261                printf("                     rhoc[%d] -> %g\n", i, rhoc[i]);
 262           if (k > 1) {
 263                for (j = 0; j < n; ++j) {
 264                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 265                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 266                     sigma[j] *= gam;
 267                     if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
 268                          sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 269                          sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 270                     }
 271                }
 272                for (j = 0; j < mma_verbose; ++j)
 273                     printf("                     sigma[%d] -> %g\n",
 274                            j, sigma[j]);
 275           }
 276      }
 277
 278  done:
 279      free(sigma);
 280      return ret;
 281 }