mlsl/mlsl.c

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