mma/mma.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include <string.h>
   4
   5 #include "mma.h"
   6
   7 #define MIN(a,b) ((a) < (b) ? (a) : (b))
   8 #define MAX(a,b) ((a) > (b) ? (a) : (b))
   9
  10 /* magic minimum value for rho in MMA ... the 2002 paper says it should
  11    be a "fixed, strictly positive `small' number, e.g. 1e-5"
  12    ... grrr, I hate these magic numbers, which seem like they
  13    should depend on the objective function in some way */
  14 #define MMA_RHOMIN 1e-5
  15
  16 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
  17                           const double *lb, const double *ub, /* bounds */
  18                           double *x, /* in: initial guess, out: minimizer */
  19                           double *minf,
  20                           nlopt_stopping *stop)
  21 {
  22      nlopt_result ret = NLOPT_SUCCESS;
  23      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
  24      int j, k = 0;
  25
  26      sigma = (double *) malloc(sizeof(double) * 6*n);
  27      if (!sigma) return NLOPT_OUT_OF_MEMORY;
  28      dfdx = sigma + n;
  29      dfdx_cur = dfdx + n;
  30      xcur = dfdx_cur + n;
  31      xprev = xcur + n;
  32      xprevprev = xprev + n;
  33      for (j = 0; j < n; ++j)
  34           sigma[j] = 0.5 * (ub[j] - lb[j]);
  35      rho = 1.0;
  36
  37      fcur = *minf = f(n, x, dfdx, f_data);
  38      memcpy(xcur, x, sizeof(double) * n);
  39      stop->nevals++;
  40      while (1) { /* outer iterations */
  41           double fprev = fcur;
  42           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
  43           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
  44           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
  45           if (ret != NLOPT_SUCCESS) goto done;
  46           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
  47           memcpy(xprev, xcur, sizeof(double) * n);
  48
  49           while (1) { /* inner iterations */
  50                double gcur = *minf, w = 0;
  51                for (j = 0; j < n; ++j) {
  52                     /* because we have no constraint functions, minimizing
  53                        the MMA approximate function can be done analytically */
  54                     double dx = -sigma[j]*sigma[j]*dfdx[j]
  55                          / (2*sigma[j]*fabs(dfdx[j]) + rho);
  56                     xcur[j] = x[j] + dx;
  57                     if (xcur[j] > x[j] + 0.9*sigma[j])
  58                          xcur[j] = x[j] + 0.9*sigma[j];
  59                     else if (xcur[j] < x[j] - 0.9*sigma[j])
  60                          xcur[j] = x[j] - 0.9*sigma[j];
  61                     if (xcur[j] > ub[j]) xcur[j] = ub[j];
  62                     else if (xcur[j] < lb[j]) xcur[j] = lb[j];
  63                     dx = xcur[j] - x[j];
  64                     gcur += (sigma[j]*sigma[j]*dfdx[j]*dx
  65                              + sigma[j]*fabs(dfdx[j])*dx*dx
  66                              + 0.5*rho*dx*dx) / (sigma[j]*sigma[j]-dx*dx);
  67                     w += 0.5*rho*dx*dx / (sigma[j]*sigma[j]-dx*dx);
  68                }
  69
  70                fcur = f(n, xcur, dfdx_cur, f_data);
  71                stop->nevals++;
  72                if (fcur < *minf) {
  73                     *minf = fcur;
  74                     memcpy(x, xcur, sizeof(double)*n);
  75                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
  76                }
  77                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
  78                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
  79                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
  80                if (ret != NLOPT_SUCCESS) goto done;
  81
  82                if (gcur >= fcur) break;
  83                rho = MIN(10*rho, 1.1 * (rho + (fcur - gcur) / w));
  84           }
  85
  86           if (nlopt_stop_ftol(stop, fcur, fprev))
  87                ret = NLOPT_FTOL_REACHED;
  88           if (nlopt_stop_x(stop, xcur, xprev))
  89                ret = NLOPT_XTOL_REACHED;
  90           if (ret != NLOPT_SUCCESS) goto done;
  91
  92           /* update rho and sigma for iteration k+1 */
  93           rho = MAX(0.1 * rho, MMA_RHOMIN);
  94           if (k > 1)
  95                for (j = 0; j < n; ++j) {
  96                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
  97                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
  98                     sigma[j] *= gam;
  99                     sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 100                     sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 101                }
 102      }
 103
 104  done:
 105      free(sigma);
 106      return ret;
 107 }