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:
7 A. H. G. Rinnooy Kan and G. T. Timmer, "Stochastic global optimization
8 methods," Mathematical Programming, vol. 39, p. 27-78 (1987).
9 [ actually 2 papers -- part I: clustering methods (p. 27), then
10 part II: multilevel methods (p. 57) ]
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).
18 Compared to the above papers, I made a couple other modifications
19 to the algorithm besides the use of a LDS.
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.)
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.
49 /* data structure for each random/quasirandom point in the population */
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) */
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 */
63 int n; /* # dimensions */
64 const double *lb, *ub;
65 nlopt_stopping *stop; /* stopping criteria */
66 nlopt_func f; void *f_data;
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) */
72 nlopt_sobol s; /* sobol data for LDS point generation, or NULL
73 to use pseudo-random numbers */
75 nlopt_algorithm local_alg; /* local search algorithm */
76 int local_maxeval; /* max # local iterations (0 if unlimited) */
78 double R_prefactor, dlm, dbound, gamma; /* parameters of MLSL */
79 int N; /* number of pts to add per iteration */
82 /* comparison routines to sort the red-black trees by function value */
83 static int pt_compare(rb_key p1_, rb_key p2_)
87 if (p1->f < p2->f) return -1;
88 if (p1->f > p2->f) return +1;
91 static int lm_compare(double *k1, double *k2)
93 if (*k1 < *k2) return -1;
94 if (*k1 > *k2) return +1;
98 /* Euclidean distance |x1 - x2|^2 */
99 static double distance2(int n, const double *x1, const double *x2)
103 for (i = 0; i < n; ++i) {
104 double dx = x1[i] - x2[i];
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)
114 rb_node *node = rb_tree_find_lt(pts, (rb_key) p);
115 double closest_d = HUGE_VAL;
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);
121 p->closest_pt_d = closest_d;
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)
128 rb_node *node = rb_tree_find_lt(lms, &p->f);
129 double closest_d = HUGE_VAL;
131 double d = distance2(n, p->x, node->k+1);
132 if (d < closest_d) closest_d = d;
133 node = rb_tree_pred(node);
135 p->closest_lm_d = closest_d;
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)
145 rb_node *node = rb_tree_find_gt(pts, (rb_key) newpt);
147 pt *p = (pt *) node->k;
149 double d = distance2(n, newpt->x, p->x);
150 if (d < p->closest_pt_d) p->closest_pt_d = d;
152 node = rb_tree_succ(node);
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)
166 node = rb_tree_find_gt(pts, (rb_key) &tmp_pt);
168 pt *p = (pt *) node->k;
170 double d = distance2(n, newlm+1, p->x);
171 if (d < p->closest_lm_d) p->closest_lm_d = d;
173 node = rb_tree_succ(node);
177 static int is_potential_minimizer(mlsl_data *mlsl, pt *p,
183 const double *lb = mlsl->lb;
184 const double *ub = mlsl->ub;
185 const double *x = p->x;
190 if (p->closest_pt_d <= dpt_min*dpt_min)
193 if (p->closest_lm_d <= dlm_min*dlm_min)
196 for (i = 0; i < n; ++i)
197 if (x[i] - lb[i] <= dbound_min || ub[i] - x[i] <= dbound_min)
203 #define K2PI (6.2831853071795864769252867665590057683943388)
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)
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 */
216 return sqrt(pow(K2PI * z, 1.0/n) * z) * exp(-0.5);
219 static pt *alloc_pt(int n)
221 pt *p = (pt *) malloc(sizeof(pt) + (n-1) * sizeof(double));
224 p->closest_pt_d = HUGE_VAL;
225 p->closest_lm_d = HUGE_VAL;
230 /* wrapper around objective function that increments evaluation count */
231 static double fcount(int n, const double *x, double *grad, void *p_)
233 mlsl_data *p = (mlsl_data *) p_;
235 return p->f(n, x, grad, p->f_data);
238 static void get_minf(mlsl_data *d, double *minf, double *x)
240 rb_node *node = rb_tree_min(&d->pts);
242 *minf = ((pt *) node->k)->f;
243 memcpy(x, ((pt *) node->k)->x, sizeof(double) * d->n);
245 node = rb_tree_min(&d->lms);
246 if (node && node->k[0] < *minf) {
248 memcpy(x, node->k + 1, sizeof(double) * d->n);
252 #define MIN(a,b) ((a) < (b) ? (a) : (b))
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?) */
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 */
261 nlopt_stopping *stop,
262 nlopt_algorithm local_alg,
264 int lds) /* random or low-discrepancy seq. (lds) */
266 nlopt_result ret = NLOPT_SUCCESS;
272 d.lb = lb; d.ub = ub;
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;
280 d.gamma = MLSL_GAMMA;
281 d.N = 4; /* FIXME: what is good number of samples per iteration? */
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);
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?) */
297 if (!p) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
299 /* FIXME: how many sobol points to skip, if any? */
300 nlopt_sobol_skip(d.s, (unsigned) (10*n+1), p->x);
302 memcpy(p->x, x, n * sizeof(double));
303 p->f = f(n, x, NULL, f_data);
305 if (!rb_tree_insert(&d.pts, (rb_key) p)) {
306 free(p); ret = NLOPT_OUT_OF_MEMORY;
308 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
309 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
310 else if (p->f < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
312 while (ret == NLOPT_SUCCESS) {
316 get_minf(&d, minf, x);
318 /* sampling phase: add random/quasi-random points */
319 for (i = 0; i < d.N && ret == NLOPT_SUCCESS; ++i) {
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 */
325 for (j = 0; j < n; ++j) p->x[j] = nlopt_urand(lb[j],ub[j]);
327 p->f = f(n, p->x, NULL, f_data);
329 if (!rb_tree_insert(&d.pts, (rb_key) p)) {
330 free(p); ret = NLOPT_OUT_OF_MEMORY;
332 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
333 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
334 else if (p->f < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
336 find_closest_pt(n, &d.pts, p);
337 find_closest_lm(n, &d.lms, p);
338 pts_update_newpt(n, &d.pts, p);
342 /* distance threshhold parameter R in MLSL */
344 * pow(log((double) d.pts.N) / d.pts.N, 1.0 / n);
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) {
351 if (is_potential_minimizer(&d, p,
352 R, d.dlm*R, d.dbound*R)) {
355 double t = nlopt_seconds();
357 if (nlopt_stop_evals(stop)) {
358 ret = NLOPT_MAXEVAL_REACHED; break;
360 if (stop->maxtime > 0 &&
361 t - stop->start >= stop->maxtime) {
362 ret = NLOPT_MAXTIME_REACHED; break;
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,
370 stop->ftol_rel, stop->ftol_abs,
371 stop->xtol_rel, stop->xtol_abs,
374 stop->maxeval - stop->nevals)
375 : stop->maxeval - stop->nevals,
376 stop->maxtime - (t - stop->start));
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;
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;
389 pts_update_newlm(n, &d.pts, lm);
392 /* TODO: additional stopping criteria based
393 e.g. on improvement in function values, etc? */
395 node = rb_tree_succ(node);
399 get_minf(&d, minf, x);
402 nlopt_sobol_destroy(d.s);
403 rb_tree_destroy_with_keys(&d.lms);
404 rb_tree_destroy_with_keys(&d.pts);