mlsl/mlsl.c

   1 /* Copyright (c) 2007-2010 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 /* Multi-Level Single-Linkage (MLSL) algorithm for
  24    global optimization by random local optimizations (a multistart algorithm
  25    with "clustering" to avoid repeated detection of the same local minimum),
  26    modified to optionally use a Sobol' low-discrepancy sequence (LDS) instead
  27    of pseudorandom numbers.  See:
  28
  29    A. H. G. Rinnooy Kan and G. T. Timmer, "Stochastic global optimization
  30    methods," Mathematical Programming, vol. 39, p. 27-78 (1987).
  31        [ actually 2 papers -- part I: clustering methods (p. 27), then
  32                               part II: multilevel methods (p. 57) ]
  33
  34    and also:
  35
  36    Sergei Kucherenko and Yury Sytsko, "Application of deterministic
  37    low-discrepancy sequences in global optimization," Computational
  38    Optimization and Applications, vol. 30, p. 297-318 (2005).
  39
  40    Compared to the above papers, I made a couple other modifications
  41    to the algorithm besides the use of a LDS.
  42
  43    1) first, the original algorithm suggests discarding points
  44       within a *fixed* distance of the boundaries or of previous
  45       local minima; I changed this to a distance that decreases with,
  46       iteration, proportional to the same threshold radius used
  47       for clustering.  (In the case of the boundaries, I make
  48       the proportionality constant rather small as well.)
  49
  50    2) Kan and Timmer suggest using fancy "spiral-search" algorithms
  51       to quickly locate nearest-neighbor points, reducing the
  52       overall time for N sample points from O(N^2) to O(N log N)
  53       However, recent papers seem to show that such schemes (kd trees,
  54       etcetera) become inefficient for high-dimensional spaces (d > 20),
  55       and finding a better approach seems to be an open question.  Therefore,
  56       since I am mainly interested in MLSL for high-dimensional problems
  57       (in low dimensions we have DIRECT etc.), I punted on this question
  58       for now.  In practice, O(N^2) (which does not count the #points
  59       evaluated in local searches) does not seem too bad if the objective
  60       function is expensive.
  61
  62 */
  63
  64 #include <stdlib.h>
  65 #include <math.h>
  66 #include <string.h>
  67
  68 #include "redblack.h"
  69 #include "mlsl.h"
  70
  71 /* data structure for each random/quasirandom point in the population */
  72 typedef struct {
  73      double f; /* function value at x */
  74      int minimized; /* if we have already minimized starting from x */
  75      double closest_pt_d; /* distance^2 to closest pt with smaller f */
  76      double closest_lm_d; /* distance^2 to closest lm with smaller f*/
  77      double x[1]; /* array of length n (K&R struct hack) */
  78 } pt;
  79
  80 /* all of the data used by the various mlsl routines...it's
  81    not clear in hindsight that we need to put all of this in a data
  82    structure since most of the work occurs in a single routine,
  83    but it doesn't hurt us */
  84 typedef struct {
  85      int n; /* # dimensions */
  86      const double *lb, *ub;
  87      nlopt_stopping *stop; /* stopping criteria */
  88      nlopt_func f; void *f_data;
  89
  90      rb_tree pts; /* tree of points (k == pt), sorted by f */
  91      rb_tree lms; /* tree of local minimizers, sorted by function value
  92                      (k = array of length d+1, [0] = f, [1..d] = x) */
  93
  94      nlopt_sobol s; /* sobol data for LDS point generation, or NULL
  95                        to use pseudo-random numbers */
  96
  97      double R_prefactor, dlm, dbound, gamma; /* parameters of MLSL */
  98      int N; /* number of pts to add per iteration */
  99 } mlsl_data;
 100
 101 /* comparison routines to sort the red-black trees by function value */
 102 static int pt_compare(rb_key p1_, rb_key p2_)
 103 {
 104      pt *p1 = (pt *) p1_;
 105      pt *p2 = (pt *) p2_;
 106      if (p1->f < p2->f) return -1;
 107      if (p1->f > p2->f) return +1;
 108      return 0;
 109 }
 110 static int lm_compare(double *k1, double *k2)
 111 {
 112      if (*k1 < *k2) return -1;
 113      if (*k1 > *k2) return +1;
 114      return 0;
 115 }
 116
 117 /* Euclidean distance |x1 - x2|^2 */
 118 static double distance2(int n, const double *x1, const double *x2)
 119 {
 120      int i;
 121      double d = 0.;
 122      for (i = 0; i < n; ++i) {
 123           double dx = x1[i] - x2[i];
 124           d += dx * dx;
 125      }
 126      return d;
 127 }
 128
 129 /* find the closest pt to p with a smaller function value;
 130    this function is called when p is first added to our tree */
 131 static void find_closest_pt(int n, rb_tree *pts, pt *p)
 132 {
 133      rb_node *node = rb_tree_find_lt(pts, (rb_key) p);
 134      double closest_d = HUGE_VAL;
 135      while (node) {
 136           double d = distance2(n, p->x, ((pt *) node->k)->x);
 137           if (d < closest_d) closest_d = d;
 138           node = rb_tree_pred(node);
 139      }
 140      p->closest_pt_d = closest_d;
 141 }
 142
 143 /* find the closest local minimizer (lm) to p with a smaller function value;
 144    this function is called when p is first added to our tree */
 145 static void find_closest_lm(int n, rb_tree *lms, pt *p)
 146 {
 147      rb_node *node = rb_tree_find_lt(lms, &p->f);
 148      double closest_d = HUGE_VAL;
 149      while (node) {
 150           double d = distance2(n, p->x, node->k+1);
 151           if (d < closest_d) closest_d = d;
 152           node = rb_tree_pred(node);
 153      }
 154      p->closest_lm_d = closest_d;
 155 }
 156
 157 /* given newpt, which presumably has just been added to our
 158    tree, update all pts with a greater function value in case
 159    newpt is closer to them than their previous closest_pt ...
 160    we can ignore already-minimized points since we never do
 161    local minimization from the same point twice */
 162 static void pts_update_newpt(int n, rb_tree *pts, pt *newpt)
 163 {
 164      rb_node *node = rb_tree_find_gt(pts, (rb_key) newpt);
 165      while (node) {
 166           pt *p = (pt *) node->k;
 167           if (!p->minimized) {
 168                double d = distance2(n, newpt->x, p->x);
 169                if (d < p->closest_pt_d) p->closest_pt_d = d;
 170           }
 171           node = rb_tree_succ(node);
 172      }
 173 }
 174
 175 /* given newlm, which presumably has just been added to our
 176    tree, update all pts with a greater function value in case
 177    newlm is closer to them than their previous closest_lm ...
 178    we can ignore already-minimized points since we never do
 179    local minimization from the same point twice */
 180 static void pts_update_newlm(int n, rb_tree *pts, double *newlm)
 181 {
 182      pt tmp_pt;
 183      rb_node *node;
 184      tmp_pt.f = newlm[0];
 185      node = rb_tree_find_gt(pts, (rb_key) &tmp_pt);
 186      while (node) {
 187           pt *p = (pt *) node->k;
 188           if (!p->minimized) {
 189                double d = distance2(n, newlm+1, p->x);
 190                if (d < p->closest_lm_d) p->closest_lm_d = d;
 191           }
 192           node = rb_tree_succ(node);
 193      }
 194 }
 195
 196 static int is_potential_minimizer(mlsl_data *mlsl, pt *p,
 197                                   double dpt_min,
 198                                   double dlm_min,
 199                                   double dbound_min)
 200 {
 201      int i, n = mlsl->n;
 202      const double *lb = mlsl->lb;
 203      const double *ub = mlsl->ub;
 204      const double *x = p->x;
 205
 206      if (p->minimized)
 207           return 0;
 208
 209      if (p->closest_pt_d <= dpt_min*dpt_min)
 210           return 0;
 211
 212      if (p->closest_lm_d <= dlm_min*dlm_min)
 213           return 0;
 214
 215      for (i = 0; i < n; ++i)
 216           if (x[i] - lb[i] <= dbound_min || ub[i] - x[i] <= dbound_min)
 217                return 0;
 218
 219      return 1;
 220 }
 221
 222 #define K2PI (6.2831853071795864769252867665590057683943388)
 223
 224 /* compute Gamma(1 + n/2)^{1/n} ... this is a constant factor
 225    used in MLSL as part of the minimum-distance cutoff*/
 226 static double gam(int n)
 227 {
 228      /* use Stirling's approximation:
 229         Gamma(1 + z) ~ sqrt(2*pi*z) * z^z / exp(z)
 230         ... this is more than accurate enough for our purposes
 231             (error in gam < 15% for d=1, < 4% for d=2, < 2% for d=3, ...)
 232             and avoids overflow problems because we can combine it with
 233             the ^{1/n} exponent */
 234      double z = n/2;
 235      return sqrt(pow(K2PI * z, 1.0/n) * z) * exp(-0.5);
 236 }
 237
 238 static pt *alloc_pt(int n)
 239 {
 240      pt *p = (pt *) malloc(sizeof(pt) + (n-1) * sizeof(double));
 241      if (p) {
 242           p->minimized = 0;
 243           p->closest_pt_d = HUGE_VAL;
 244           p->closest_lm_d = HUGE_VAL;
 245      }
 246      return p;
 247 }
 248
 249 /* wrapper around objective function that increments evaluation count */
 250 static double fcount(unsigned n, const double *x, double *grad, void *p_)
 251 {
 252      mlsl_data *p = (mlsl_data *) p_;
 253      p->stop->nevals++;
 254      return p->f(n, x, grad, p->f_data);
 255 }
 256
 257 static void get_minf(mlsl_data *d, double *minf, double *x)
 258 {
 259      rb_node *node = rb_tree_min(&d->pts);
 260      if (node) {
 261           *minf = ((pt *) node->k)->f;
 262           memcpy(x, ((pt *) node->k)->x, sizeof(double) * d->n);
 263      }
 264      node = rb_tree_min(&d->lms);
 265      if (node && node->k[0] < *minf) {
 266           *minf = node->k[0];
 267           memcpy(x, node->k + 1, sizeof(double) * d->n);
 268      }
 269 }
 270
 271 #define MIN(a,b) ((a) < (b) ? (a) : (b))
 272
 273 #define MLSL_SIGMA 2. /* MLSL sigma parameter, using value from the papers */
 274 #define MLSL_GAMMA 0.3 /* MLSL gamma parameter (FIXME: what should it be?) */
 275
 276 nlopt_result mlsl_minimize(int n, nlopt_func f, void *f_data,
 277                            const double *lb, const double *ub, /* bounds */
 278                            double *x, /* in: initial guess, out: minimizer */
 279                            double *minf,
 280                            nlopt_stopping *stop,
 281                            nlopt_opt local_opt,
 282                            int Nsamples, /* #samples/iteration (0=default) */
 283                            int lds) /* random or low-discrepancy seq. (lds) */
 284 {
 285      nlopt_result ret = NLOPT_SUCCESS;
 286      mlsl_data d;
 287      int i;
 288      pt *p;
 289
 290      if (!Nsamples)
 291           d.N = 4; /* FIXME: what is good number of samples per iteration? */
 292      else
 293           d.N = Nsamples;
 294      if (d.N < 1) return NLOPT_INVALID_ARGS;
 295
 296      d.n = n;
 297      d.lb = lb; d.ub = ub;
 298      d.stop = stop;
 299      d.f = f; d.f_data = f_data;
 300      rb_tree_init(&d.pts, pt_compare);
 301      rb_tree_init(&d.lms, lm_compare);
 302      d.s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
 303
 304      nlopt_set_min_objective(local_opt, fcount, &d);
 305      nlopt_set_lower_bounds(local_opt, lb);
 306      nlopt_set_upper_bounds(local_opt, ub);
 307      nlopt_set_stopval(local_opt, stop->minf_max);
 308
 309      d.gamma = MLSL_GAMMA;
 310
 311      d.R_prefactor = sqrt(2./K2PI) * pow(gam(n) * MLSL_SIGMA, 1.0/n);
 312      for (i = 0; i < n; ++i)
 313           d.R_prefactor *= pow(ub[i] - lb[i], 1.0/n);
 314
 315      /* MLSL also suggests setting some minimum distance from points
 316         to previous local minimiza and to the boundaries; I don't know
 317         how to choose this as a fixed number, so I set it proportional
 318         to R; see also the comments at top.  dlm and dbound are the
 319         proportionality constants. */
 320      d.dlm = 1.0; /* min distance/R to local minima (FIXME: good value?) */
 321      d.dbound = 1e-6; /* min distance/R to ub/lb boundaries (good value?) */
 322
 323
 324      p = alloc_pt(n);
 325      if (!p) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 326
 327      /* FIXME: how many sobol points to skip, if any? */
 328      nlopt_sobol_skip(d.s, (unsigned) (10*n+d.N), p->x);
 329
 330      memcpy(p->x, x, n * sizeof(double));
 331      p->f = f(n, x, NULL, f_data);
 332      stop->nevals++;
 333      if (!rb_tree_insert(&d.pts, (rb_key) p)) {
 334           free(p); ret = NLOPT_OUT_OF_MEMORY;
 335      }
 336      if (nlopt_stop_forced(stop)) ret = NLOPT_FORCE_STOP;
 337      else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 338      else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 339      else if (p->f < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 340
 341      while (ret == NLOPT_SUCCESS) {
 342           rb_node *node;
 343           double R;
 344
 345           get_minf(&d, minf, x);
 346
 347           /* sampling phase: add random/quasi-random points */
 348           for (i = 0; i < d.N && ret == NLOPT_SUCCESS; ++i) {
 349                p = alloc_pt(n);
 350                if (!p) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 351                if (d.s) nlopt_sobol_next(d.s, p->x, lb, ub);
 352                else { /* use random points instead of LDS */
 353                     int j;
 354                     for (j = 0; j < n; ++j) p->x[j] = nlopt_urand(lb[j],ub[j]);
 355                }
 356                p->f = f(n, p->x, NULL, f_data);
 357                stop->nevals++;
 358                if (!rb_tree_insert(&d.pts, (rb_key) p)) {
 359                     free(p); ret = NLOPT_OUT_OF_MEMORY;
 360                }
 361                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCE_STOP;
 362                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 363                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 364                else if (p->f < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 365                else {
 366                     find_closest_pt(n, &d.pts, p);
 367                     find_closest_lm(n, &d.lms, p);
 368                     pts_update_newpt(n, &d.pts, p);
 369                }
 370           }
 371
 372           /* distance threshhold parameter R in MLSL */
 373           R = d.R_prefactor
 374                * pow(log((double) d.pts.N) / d.pts.N, 1.0 / n);
 375
 376           /* local search phase: do local opt. for promising points */
 377           node = rb_tree_min(&d.pts);
 378           for (i = (int) (ceil(d.gamma * d.pts.N) + 0.5);
 379                node && i > 0 && ret == NLOPT_SUCCESS; --i) {
 380                p = (pt *) node->k;
 381                if (is_potential_minimizer(&d, p,
 382                                           R, d.dlm*R, d.dbound*R)) {
 383                     nlopt_result lret;
 384                     double *lm;
 385                     double t = nlopt_seconds();
 386
 387                     if (nlopt_stop_forced(stop)) {
 388                          ret = NLOPT_FORCE_STOP; break;
 389                     }
 390                     if (nlopt_stop_evals(stop)) {
 391                          ret = NLOPT_MAXEVAL_REACHED; break;
 392                     }
 393                     if (stop->maxtime > 0 &&
 394                         t - stop->start >= stop->maxtime) {
 395                          ret = NLOPT_MAXTIME_REACHED; break;
 396                     }
 397                     lm = (double *) malloc(sizeof(double) * (n+1));
 398                     if (!lm) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 399                     memcpy(lm+1, p->x, sizeof(double) * n);
 400                     lret = nlopt_optimize_limited(local_opt, lm+1, lm,
 401                                                   stop->maxeval - stop->nevals,
 402                                                   stop->maxtime -
 403                                                   (t - stop->start));
 404                     p->minimized = 1;
 405                     if (lret < 0) { free(lm); ret = lret; goto done; }
 406                     if (!rb_tree_insert(&d.lms, lm)) {
 407                          free(lm); ret = NLOPT_OUT_OF_MEMORY;
 408                     }
 409                     else if (nlopt_stop_forced(stop)) ret = NLOPT_FORCE_STOP;
 410                     else if (*lm < stop->minf_max)
 411                          ret = NLOPT_MINF_MAX_REACHED;
 412                     else if (nlopt_stop_evals(stop))
 413                          ret = NLOPT_MAXEVAL_REACHED;
 414                     else if (nlopt_stop_time(stop))
 415                          ret = NLOPT_MAXTIME_REACHED;
 416                     else
 417                          pts_update_newlm(n, &d.pts, lm);
 418                }
 419
 420                /* TODO: additional stopping criteria based
 421                   e.g. on improvement in function values, etc? */
 422
 423                node = rb_tree_succ(node);
 424           }
 425      }
 426
 427      get_minf(&d, minf, x);
 428
 429  done:
 430      nlopt_sobol_destroy(d.s);
 431      rb_tree_destroy_with_keys(&d.lms);
 432      rb_tree_destroy_with_keys(&d.pts);
 433      return ret;
 434 }