From: stevenj Date: Sat, 8 Nov 2008 21:56:21 +0000 (-0500) Subject: added nelder-mead simplex algorithm (in preparation for re-implementing subplex,... X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=3176a5eb11e4fa8bae6898172cd95df080d2a8ee;p=nlopt.git added nelder-mead simplex algorithm (in preparation for re-implementing subplex, and possibly recent provably convergent variants of nelder-meade) darcs-hash:20081108215621-c8de0-91e9b421d37162aa2b8f4013f1c88a4003d05215.gz --- diff --git a/Makefile.am b/Makefile.am index 7045113..52c438a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -8,7 +8,7 @@ CXX_DIRS = stogo CXX_LIBS = stogo/libstogo.la endif -SUBDIRS= util subplex direct cdirect $(CXX_DIRS) praxis luksan crs mlsl mma cobyla newuoa lbfgs api . octave test +SUBDIRS= util subplex direct cdirect $(CXX_DIRS) praxis luksan crs mlsl mma cobyla newuoa lbfgs neldermead api . octave test EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4 if WITH_NOCEDAL @@ -19,7 +19,7 @@ libnlopt@NLOPT_SUFFIX@_la_SOURCES = libnlopt@NLOPT_SUFFIX@_la_LIBADD = subplex/libsubplex.la \ direct/libdirect.la cdirect/libcdirect.la $(CXX_LIBS) \ praxis/libpraxis.la $(NOCEDAL_LBFGS) luksan/libluksan.la crs/libcrs.la \ -mlsl/libmlsl.la mma/libmma.la cobyla/libcobyla.la newuoa/libnewuoa.la api/libapi.la util/libutil.la +mlsl/libmlsl.la mma/libmma.la cobyla/libcobyla.la newuoa/libnewuoa.la neldermead/libneldermead.la api/libapi.la util/libutil.la libnlopt@NLOPT_SUFFIX@_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@ diff --git a/api/Makefile.am b/api/Makefile.am index 1ec9af2..8754302 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -1,4 +1,4 @@ -AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/lbfgs -I$(top_srcdir)/luksan -I$(top_srcdir)/crs -I$(top_srcdir)/mlsl -I$(top_srcdir)/mma -I$(top_srcdir)/cobyla -I$(top_srcdir)/newuoa -I$(top_srcdir)/util +AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/lbfgs -I$(top_srcdir)/luksan -I$(top_srcdir)/crs -I$(top_srcdir)/mlsl -I$(top_srcdir)/mma -I$(top_srcdir)/cobyla -I$(top_srcdir)/newuoa -I$(top_srcdir)/neldermead -I$(top_srcdir)/util include_HEADERS = nlopt.h nlopt.f noinst_LTLIBRARIES = libapi.la diff --git a/api/nlopt.c b/api/nlopt.c index 9109a71..ea2a438 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -101,7 +101,8 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "Method of Moving Asymptotes (MMA) (local, derivative)", "COBYLA (Constrained Optimization BY Linear Approximations) (local, no-derivative)", "NEWUOA unconstrained optimization via quadratic models (local, no-derivative)", - "Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)" + "Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)", + "Nelder-Mead simplex algorithm" }; const char *nlopt_algorithm_name(nlopt_algorithm a) @@ -223,6 +224,7 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) #include "mma.h" #include "cobyla.h" #include "newuoa.h" +#include "neldermead.h" /*************************************************************************/ @@ -486,6 +488,16 @@ static nlopt_result nlopt_minimize_( return newuoa(n, 2*n+1, x, lb, ub, initial_step(n, lb, ub, x), &stop, minf, f_noderiv, &d); + case NLOPT_LN_NELDERMEAD: { + nlopt_result ret; + double *scale = (double *) malloc(sizeof(double) * n); + if (!scale) return NLOPT_OUT_OF_MEMORY; + for (i = 0; i < n; ++i) + scale[i] = initial_step(1, lb+i, ub+i, x+i); + ret = nldrmd_minimize(n, f,f_data, lb,ub, x, minf, scale, &stop); + free(scale); + return ret; + } default: return NLOPT_INVALID_ARGS; diff --git a/api/nlopt.h b/api/nlopt.h index 7692bb0..327f67e 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -90,6 +90,8 @@ typedef enum { NLOPT_LN_NEWUOA, NLOPT_LN_NEWUOA_BOUND, + NLOPT_LN_NELDERMEAD, + NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; diff --git a/configure.ac b/configure.ac index b461fdc..0f50a2d 100644 --- a/configure.ac +++ b/configure.ac @@ -214,6 +214,7 @@ AC_CONFIG_FILES([ mma/Makefile cobyla/Makefile newuoa/Makefile + neldermead/Makefile test/Makefile ]) diff --git a/neldermead/Makefile.am b/neldermead/Makefile.am new file mode 100644 index 0000000..527cf7f --- /dev/null +++ b/neldermead/Makefile.am @@ -0,0 +1,6 @@ +AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api + +noinst_LTLIBRARIES = libneldermead.la +libneldermead_la_SOURCES = nldrmd.c neldermead.h + +EXTRA_DIST = README diff --git a/neldermead/README b/neldermead/README new file mode 100644 index 0000000..6d71e9a --- /dev/null +++ b/neldermead/README @@ -0,0 +1,49 @@ +Nelder-Mead and variations thereof. Possibly the algorithms: + +----------------- + +First, the original Nelder-Mead algorithm. + +----------------- + +Second, the provably convergent variant of Nelder-Mead described in: + + C. J. Price, I. D. Coope, and D. Byatt, "A convergent variant + of the Nelder-Mead algorithm," J. Optim. Theory Appl. 113 (1), + p. 5-19 (2002). + +And/or possibly the (claimed superior) one in: + + A. Burmen, J. Puhan, and T. Tuma, "Grid restrained Nelder-Mead + algorithm," Computational Optim. Appl. 34(3), 359-375 (2006). + +----------------- + +My own independent implemention of Tom Rowan's Subplex algorithm (a +more-efficient variant of Nelder-Mead simplex), which was described at: + + http://www.netlib.org/opt/subplex.tgz + + T. Rowan, "Functional Stability Analysis of Numerical Algorithms", + Ph.D. thesis, Department of Computer Sciences, University of Texas + at Austin, 1990. + +I would have liked to use Rowan's original implementation, but its +legal status is unfortunately unclear. Rowan didn't include any +license statement at all with the original code, which makes it +technically illegal to redistribute. I contacted Rowan about getting +a clear open-source/free-software license for it, and he was very +agreeable, but he said he had to think about the specific license +choice and would get back to me. Unfortunately, a year later I still +haven't heard from him, and his old email address no longer seems to +work, so I don't know how to contact him for permission. + +Although I now have other derivative-free optimization routines in +NLopt, the subplex algorithm is nice to have because it is somewhat +tolerant of discontinuous and/or noisy objectives, which may make it a +good choice for some problems. + +Tom Rowan expressed a preference that modified versions of his code +use a different name from "subplex". Since this is a complete +from-scratch re-implementation, I figured that he would want a +different name too, so I am calling it "sbplx". diff --git a/neldermead/neldermead.h b/neldermead/neldermead.h new file mode 100644 index 0000000..0b19d8d --- /dev/null +++ b/neldermead/neldermead.h @@ -0,0 +1,46 @@ +/* Copyright (c) 2007-2008 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#ifndef NELDERMEAD_H +#define NELDERMEAD_H + +#include "nlopt.h" +#include "nlopt-util.h" + +#ifdef __cplusplus +extern "C" +{ +#endif /* __cplusplus */ + +nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, + const double *xstep, /* initial step sizes */ + nlopt_stopping *stop); + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif + diff --git a/neldermead/nldrmd.c b/neldermead/nldrmd.c new file mode 100644 index 0000000..2a2b5b9 --- /dev/null +++ b/neldermead/nldrmd.c @@ -0,0 +1,209 @@ +/* Copyright (c) 2007-2008 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include + +#include "neldermead.h" +#include "redblack.h" + +/* Nelder-Mead simplex algorithm, used as a subroutine for the Rowan's + subplex algorithm. Modified to handle bound constraints ala + Richardson and Kuester (1973), as mentioned below. */ + +/* heuristic "strategy" constants: */ +static const double alpha = 1, beta = 0.5, gamm = 2, delta = 0.5; + +/* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */ +static int simplex_compare(double *k1, double *k2) +{ + if (*k1 < *k2) return -1; + if (*k1 > *k2) return +1; + return k1 - k2; /* tie-breaker */ +} + +/* pin x to lie within the given lower and upper bounds; this is + probably a suboptimal way to handle bound constraints, but I don't + know a better way and it is the method suggested by J. A. Richardson + and J. L. Kuester, "The complex method for constrained optimization," + Commun. ACM 16(8), 487-489 (1973). */ +static void pin(int n, double *x, const double *lb, const double *ub) { + int i; + for (i = 0; i < n; ++i) { + if (x[i] < lb[i]) x[i] = lb[i]; + else if (x[i] > ub[i]) x[i] = ub[i]; + } +} + + +#define CHECK_EVAL(xc,fc) \ + if ((fc) <= *minf) { \ + *minf = (fc); memcpy(x, (xc), n * sizeof(double)); \ + if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; } \ + } \ + stop->nevals++; \ + if (nlopt_stop_evals(stop)) { ret=NLOPT_MAXEVAL_REACHED; goto done; } \ + if (nlopt_stop_time(stop)) { ret=NLOPT_MAXTIME_REACHED; goto done; } + +nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, + const double *xstep, /* initial step sizes */ + nlopt_stopping *stop) +{ + double *pts; /* (n+1) x (n+1) array of n+1 points plus function val [0] */ + double *c; /* centroid * n */ + double *xcur; /* current point */ + rb_tree t; /* red-black tree of simplex, sorted by f(x) */ + int i, j; + double ninv = 1.0 / n; + nlopt_result ret = NLOPT_SUCCESS; + + pts = (double*) malloc(sizeof(double) * (n+1) * (n+1)); + if (!pts) return NLOPT_OUT_OF_MEMORY; + c = (double*) malloc(sizeof(double) * n * 2); + if (!c) { free(pts); return NLOPT_OUT_OF_MEMORY; } + xcur = c + n; + + /* initialize the simplex based on the starting xstep */ + memcpy(pts+1, x, sizeof(double)*n); + *minf = pts[0] = f(n, pts+1, NULL, f_data); + CHECK_EVAL(pts+1, pts[0]); + for (i = 0; i < n; ++i) { + double *pt = pts + (i+1)*(n+1); + memcpy(pt+1, x, sizeof(double)*n); + if (x[i] + xstep[i] <= ub[i]) + pt[1+i] += xstep[i]; + else if (x[i] - xstep[i] >= lb[i]) + pt[1+i] -= xstep[i]; + else if (ub[i] - x[i] > x[i] - lb[i]) + pt[1+i] += 0.5 * (ub[i] - x[i]); + else + pt[1+i] -= 0.5 * (x[i] - lb[i]); + pt[0] = f(n, pt+1, NULL, f_data); + CHECK_EVAL(pt+1, pt[0]); + } + + rb_tree_init(&t, simplex_compare); + restart: + for (i = 0; i < n + 1; ++i) + if (!rb_tree_insert(&t, pts + i*(n+1))) { + ret = NLOPT_OUT_OF_MEMORY; + goto done; + } + + while (1) { + rb_node *low = rb_tree_min(&t); + rb_node *high = rb_tree_max(&t); + double fl = low->k[0], *xl = low->k + 1; + double fh = high->k[0], *xh = high->k + 1; + double fr; + + if (nlopt_stop_f(stop, fl, fh)) ret = NLOPT_FTOL_REACHED; + + /* compute centroid ... if we cared about the perfomance of this, + we could do it iteratively by updating the centroid on + each step, but then we would have to be more careful about + accumulation of rounding errors... anyway n is unlikely to + be very large for Nelder-Mead in practical cases */ + memset(c, 0, sizeof(double)*n); + for (i = 0; i < n + 1; ++i) { + double *xi = pts + i*(n+1) + 1; + if (xi != xh) + for (j = 0; j < n; ++j) + c[j] += xi[j]; + } + for (i = 0; i < n; ++i) c[i] *= ninv; + + /* x convergence check: find xcur = max radius from centroid */ + memset(xcur, 0, sizeof(double)*n); + for (i = 0; i < n + 1; ++i) { + double *xi = pts + i*(n+1) + 1; + for (j = 0; j < n; ++j) { + double dx = fabs(xi[j] - c[j]); + if (dx > xcur[j]) xcur[j] = dx; + } + } + for (i = 0; i < n; ++i) xcur[i] += c[i]; + if (nlopt_stop_x(stop, c, xcur)) ret = NLOPT_XTOL_REACHED; + + /* reflection */ + for (i = 0; i < n; ++i) xcur[i] = c[i] + alpha * (c[i] - xh[i]); + pin(n, xcur, lb, ub); + fr = f(n, xcur, NULL, f_data); + CHECK_EVAL(xcur, fr); + + if (fr < fl) { /* new best point, expand simplex */ + for (i = 0; i < n; ++i) xh[i] = c[i] + gamm * (c[i] - xh[i]); + pin(n, xh, lb, ub); + fh = f(n, xh, NULL, f_data); + CHECK_EVAL(xh, fh); + if (fh >= fr) { /* expanding didn't improve */ + fh = fr; + memcpy(xh, xcur, sizeof(double)*n); + } + } + else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */ + memcpy(xh, xcur, sizeof(double)*n); + fh = fr; + } + else { /* new worst point, contract */ + double fc; + if (fh <= fr) + for (i = 0; i < n; ++i) xcur[i] = c[i] - beta*(c[i]-xh[i]); + else + for (i = 0; i < n; ++i) xcur[i] = c[i] + beta*(c[i]-xh[i]); + pin(n, xcur, lb, ub); + fc = f(n, xcur, NULL, f_data); + CHECK_EVAL(xcur, fc); + if (fc < fr && fc < fh) { /* successful contraction */ + memcpy(xh, xcur, sizeof(double)*n); + fh = fc; + } + else { /* failed contraction, shrink simplex */ + rb_tree_destroy(&t); + rb_tree_init(&t, simplex_compare); + for (i = 0; i < n+1; ++i) { + double *pt = pts + i * (n+1); + if (pt+1 != xl) { + for (j = 0; j < n; ++j) + pt[1+j] = xl[j] + delta * (pt[1+j] - xl[j]); + pt[0] = f(n, pt+1, NULL, f_data); + CHECK_EVAL(pt+1, pt[0]); + } + } + goto restart; + } + } + + high->k[0] = fh; + rb_tree_resort(&t, high); + } + +done: + rb_tree_destroy(&t); + free(c); + free(pts); + return ret; +}