From 2a873952cee721218e16cd65cf73818cf7d68e20 Mon Sep 17 00:00:00 2001 From: stevenj Date: Mon, 10 Nov 2008 17:03:44 -0500 Subject: [PATCH] added sbplx darcs-hash:20081110220344-c8de0-a64c0038ba3f2220e4c130c577d772286bbf3a53.gz --- api/nlopt.c | 12 ++- api/nlopt.h | 1 + neldermead/Makefile.am | 2 +- neldermead/neldermead.h | 7 ++ neldermead/sbplx.c | 222 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 240 insertions(+), 4 deletions(-) create mode 100644 neldermead/sbplx.c diff --git a/api/nlopt.c b/api/nlopt.c index 757b53f..bc191e5 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -102,7 +102,8 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "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)", - "Nelder-Mead simplex algorithm" + "Nelder-Mead simplex algorithm (local, no-derivative)", + "Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)" }; const char *nlopt_algorithm_name(nlopt_algorithm a) @@ -488,13 +489,18 @@ 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: { + case NLOPT_LN_NELDERMEAD: + case NLOPT_LN_SBPLX: + { 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); + if (algorithm == NLOPT_LN_NELDERMEAD) + ret = nldrmd_minimize(n,f,f_data,lb,ub,x,minf,scale,&stop); + else + ret = sbplx_minimize(n,f,f_data,lb,ub,x,minf,scale,&stop); free(scale); return ret; } diff --git a/api/nlopt.h b/api/nlopt.h index 327f67e..88d4578 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -91,6 +91,7 @@ typedef enum { NLOPT_LN_NEWUOA_BOUND, NLOPT_LN_NELDERMEAD, + NLOPT_LN_SBPLX, NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; diff --git a/neldermead/Makefile.am b/neldermead/Makefile.am index 527cf7f..1c9ed8a 100644 --- a/neldermead/Makefile.am +++ b/neldermead/Makefile.am @@ -1,6 +1,6 @@ AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api noinst_LTLIBRARIES = libneldermead.la -libneldermead_la_SOURCES = nldrmd.c neldermead.h +libneldermead_la_SOURCES = nldrmd.c neldermead.h sbplx.c EXTRA_DIST = README diff --git a/neldermead/neldermead.h b/neldermead/neldermead.h index 8ec8bfc..f1f35f5 100644 --- a/neldermead/neldermead.h +++ b/neldermead/neldermead.h @@ -46,6 +46,13 @@ nlopt_result nldrmd_minimize_(int n, nlopt_func f, void *f_data, nlopt_stopping *stop, double psi, double *scratch); +nlopt_result sbplx_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 *xstep0, /* initial step sizes */ + nlopt_stopping *stop); + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ diff --git a/neldermead/sbplx.c b/neldermead/sbplx.c new file mode 100644 index 0000000..cc39941 --- /dev/null +++ b/neldermead/sbplx.c @@ -0,0 +1,222 @@ +/* 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 + +#include "neldermead.h" + +/* subplex strategy constants: */ +static const double psi = 0.25, omega = 0.1; +static const int nsmin = 2, nsmax = 5; + +int sbplx_verbose = 0; /* for debugging */ + +/* qsort_r comparison function for sorting indices into delta-x array */ +static int p_compare(void *dx_, const void *i_, const void *j_) +{ + const double *dx = (const double *) dx_; + int i = *((const int *) i_), j = *((const int *) j_); + double dxi = fabs(dx[i]), dxj = fabs(dx[j]); + return (dxi > dxj ? -1 : (dxi < dxj ? +1 : 0)); +} + +typedef struct { + const int *p; /* subspace index permutation */ + int is; /* starting index for this subspace */ + int n; /* dimension of underlying space */ + double *x; /* current x vector */ + nlopt_func f; void *f_data; /* the "actual" underlying function */ +} subspace_data; + +/* wrapper around objective function for subspace optimization */ +static double subspace_func(int ns, const double *xs, double *grad, void *data) +{ + subspace_data *d = (subspace_data *) data; + int i, is = d->is; + const int *p = d->p; + double *x = d->x; + + (void) grad; /* should always be NULL here */ + for (i = is; i < is + ns; ++i) x[p[i]] = xs[i-is]; + return d->f(d->n, x, NULL, d->f_data); +} + +nlopt_result sbplx_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 *xstep0, /* initial step sizes */ + nlopt_stopping *stop) +{ + nlopt_result ret = NLOPT_SUCCESS; + double *xstep, *xprev, *dx, *xs, *lbs, *ubs, *xsstep, *scratch; + int *p; /* permuted indices of x sorted by decreasing magnitude |dx| */ + int i; + subspace_data sd; + double fprev; + + *minf = f(n, x, NULL, f_data); + stop->nevals++; + if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED; + if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED; + if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED; + + xstep = (double*)malloc(sizeof(double) * (n*3 + nsmax*4 + + (nsmax+1)*(nsmax+1)+2*nsmax)); + if (!xstep) return NLOPT_OUT_OF_MEMORY; + xprev = xstep + n; dx = xprev + n; + xs = dx + n; xsstep = xs + nsmax; + lbs = xsstep + nsmax; ubs = lbs + nsmax; + scratch = ubs + nsmax; + p = (int *) malloc(sizeof(int) * n); + if (!p) { free(xstep); return NLOPT_OUT_OF_MEMORY; } + + memcpy(xstep, xstep0, n * sizeof(double)); + memset(dx, 0, n * sizeof(double)); + + sd.p = p; + sd.n = n; + sd.x = x; + sd.f = f; + sd.f_data = f_data; + + while (1) { + double normi = 0; + double normdx = 0; + int ns, nsubs = 0; + int nevals = stop->nevals; + + memcpy(xprev, x, n * sizeof(double)); + fprev = *minf; + + /* sort indices into the progress vector dx by decreasing + order of magnitude |dx| */ + for (i = 0; i < n; ++i) p[i] = i; + nlopt_qsort_r(p, (size_t) n, sizeof(int), dx, p_compare); + + /* find the subspaces, and perform nelder-mead on each one */ + for (i = 0; i < n; ++i) normdx += fabs(dx[i]); /* L1 norm */ + for (i = 0; i + nsmin < n; i += ns) { + /* find subspace starting at index i */ + int k, nk; + double ns_goodness = -HUGE_VAL, norm = normi; + nk = i+nsmax > n ? n : i+nsmax; /* max k for this subspace */ + for (k = i; k < i+nsmin-1; ++k) norm += fabs(dx[p[k]]); + ns = nsmin; + for (k = i+nsmin-1; k < nk; ++k) { + double goodness; + norm += fabs(dx[p[k]]); + /* remaining subspaces must be big enough to partition */ + if (n-(k+1) < nsmin) continue; + /* maximize figure of merit defined by Rowan thesis: + look for sudden drops in average |dx| */ + if (k+1 < n) + goodness = norm/(k+1) - (normdx-norm)/(n-(k+1)); + else + goodness = normdx/n; + if (goodness > ns_goodness) { + ns_goodness = goodness; + ns = (k+1)-i; + } + } + for (k = i; k < i+ns; ++k) normi += fabs(dx[p[k]]); + /* do nelder-mead on subspace of dimension ns starting w/i */ + sd.is = i; + for (k = i; k < i+ns; ++k) { + xs[k-i] = x[p[k]]; + xsstep[k-i] = xstep[p[k]]; + lbs[k-i] = lb[p[k]]; + ubs[k-i] = ub[p[k]]; + } + ++nsubs; + nevals = stop->nevals; + ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf, + xsstep, stop, psi, scratch); + if (sbplx_verbose) + printf("%d NM iterations for (%d,%d) subspace\n", + stop->nevals - nevals, sd.is, ns); + for (k = i; k < i+ns; ++k) x[p[k]] = xs[k-i]; + if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; } + if (ret != NLOPT_XTOL_REACHED) goto done; + } + /* nelder-mead on last subspace */ + ns = n - i; + sd.is = i; + for (; i < n; ++i) { + xs[i-sd.is] = x[p[i]]; + xsstep[i-sd.is] = xstep[p[i]]; + lbs[i-sd.is] = lb[p[i]]; + ubs[i-sd.is] = ub[p[i]]; + } + ++nsubs; + nevals = stop->nevals; + ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf, + xsstep, stop, psi, scratch); + if (sbplx_verbose) + printf("%d NM iterations for (%d,%d) subspace\n", + stop->nevals - nevals, sd.is, ns); + for (i = sd.is; i < n; ++i) x[p[i]] = xs[i-sd.is]; + if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; } + if (ret != NLOPT_XTOL_REACHED) goto done; + + /* termination tests: */ + if (nlopt_stop_f(stop, *minf, fprev)) { + ret = NLOPT_FTOL_REACHED; + goto done; + } + if (nlopt_stop_x(stop, x, xprev)) { + ret = NLOPT_XTOL_REACHED; + goto done; + } + + /* compute change in optimal point */ + for (i = 0; i < n; ++i) dx[i] = x[i] - xprev[i]; + + /* setting stepsizes */ + { + double scale; + if (nsubs == 1) + scale = psi; + else { + double stepnorm = 0, dxnorm = 0; + for (i = 0; i < n; ++i) { + stepnorm += fabs(xstep[i]); + dxnorm += fabs(dx[i]); + } + scale = dxnorm / stepnorm; + if (scale < omega) scale = omega; + if (scale > 1/omega) scale = 1/omega; + } + for (i = 0; i < n; ++i) + xstep[i] = (dx[i] == 0) ? -(xstep[i] * scale) + : copysign(xstep[i] * scale, dx[i]); + } + } + + done: + free(p); + free(xstep); + return ret; +} -- 2.30.2