From 721f1b8b16d6f7a16b8d7a9a67b883a4add3256f Mon Sep 17 00:00:00 2001 From: stevenj Date: Tue, 29 Jul 2008 02:18:36 -0400 Subject: [PATCH] added MMA implementation darcs-hash:20080729061836-c8de0-7119b5077ac025515b357e848c140f2f081ff831.gz --- Makefile.am | 10 ++--- api/Makefile.am | 2 +- api/nlopt.c | 13 ++++-- api/nlopt.h | 2 + configure.ac | 1 + mma/Makefile.am | 6 +++ mma/README | 20 +++++++++ mma/mma.c | 107 ++++++++++++++++++++++++++++++++++++++++++++++++ mma/mma.h | 23 +++++++++++ 9 files changed, 175 insertions(+), 9 deletions(-) create mode 100644 mma/Makefile.am create mode 100644 mma/README create mode 100644 mma/mma.c create mode 100644 mma/mma.h diff --git a/Makefile.am b/Makefile.am index 1a3d1b4..3570f88 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 lbfgs api . octave test +SUBDIRS= util subplex direct cdirect $(CXX_DIRS) praxis luksan crs mlsl mma lbfgs api . octave test EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4 if WITH_NOCEDAL @@ -16,10 +16,10 @@ NOCEDAL_LBFGS=lbfgs/liblbfgs.la endif libnlopt_la_SOURCES = -libnlopt_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 api/libapi.la \ -util/libutil.la +libnlopt_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 \ +api/libapi.la util/libutil.la libnlopt_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@ diff --git a/api/Makefile.am b/api/Makefile.am index 0823aaa..3a616b1 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)/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)/util include_HEADERS = nlopt.h noinst_LTLIBRARIES = libapi.la diff --git a/api/nlopt.c b/api/nlopt.c index a39010c..6ad36f3 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -59,9 +59,9 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "StoGO randomized (NOT COMPILED)", #endif #ifdef WITH_NOCEDAL_LBFGS - "original NON-FREE L-BFGS implementation by Nocedal et al. (local, deriv.-based)" + "original NON-FREE L-BFGS code by Nocedal et al. (local, deriv.-based)", #else - "original NON-FREE L-BFGS implementation by Nocedal et al. (NOT COMPILED)" + "original NON-FREE L-BFGS code by Nocedal et al. (NOT COMPILED)", #endif "Low-storage BFGS (LBFGS) (local, derivative-based)", "Principal-axis, praxis (local, no-derivative)", @@ -75,7 +75,8 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "Multi-level single-linkage (MLSL), random (global, no-derivative)", "Multi-level single-linkage (MLSL), random (global, derivative)", "Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)", - "Multi-level single-linkage (MLSL), quasi-random (global, derivative)" + "Multi-level single-linkage (MLSL), quasi-random (global, derivative)", + "Method of Moving Asymptotes (MMA) (local, derivative)" }; const char *nlopt_algorithm_name(nlopt_algorithm a) @@ -148,6 +149,9 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) #include "crs.h" +#include "mlsl.h" +#include "mma.h" + /*************************************************************************/ /* for "hybrid" algorithms that combine global search with some @@ -398,6 +402,9 @@ static nlopt_result nlopt_minimize_( local_search_maxeval, algorithm >= NLOPT_GN_MLSL_LDS); + case NLOPT_LD_MMA: + return mma_minimize(n, f, f_data, lb, ub, x, minf, &stop); + default: return NLOPT_INVALID_ARGS; } diff --git a/api/nlopt.h b/api/nlopt.h index b6e19a5..edd04c4 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -59,6 +59,8 @@ typedef enum { NLOPT_GN_MLSL_LDS, NLOPT_GD_MLSL_LDS, + NLOPT_LD_MMA, + NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; diff --git a/configure.ac b/configure.ac index 8a3b8a1..4f98d81 100644 --- a/configure.ac +++ b/configure.ac @@ -200,6 +200,7 @@ AC_CONFIG_FILES([ luksan/Makefile crs/Makefile mlsl/Makefile + mma/Makefile test/Makefile ]) diff --git a/mma/Makefile.am b/mma/Makefile.am new file mode 100644 index 0000000..9d4811e --- /dev/null +++ b/mma/Makefile.am @@ -0,0 +1,6 @@ +AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api + +noinst_LTLIBRARIES = libmma.la +libmma_la_SOURCES = mma.c mma.h + +EXTRA_DIST = README diff --git a/mma/README b/mma/README new file mode 100644 index 0000000..07fa109 --- /dev/null +++ b/mma/README @@ -0,0 +1,20 @@ +Implementation of the globally-convergent method-of-moving-asymptotes (MMA) +algorithm for gradient-based local optimization, as described in: + + Krister Svanberg, "A class of globally convergent optimization + methods based on conservative convex separable approximations," + SIAM J. Optim. 12 (2), p. 555-573 (2002). + +However, in the case of NLopt, because we have no constraints other +than the bound constraints on the input variables, MMA reduces to an +extremely simple algorithm. This doesn't mean it's a bad algorithm -- +it may be especially suited to topology optimization problems where +the optimal design variables are almost always slammed against their +constraints -- just that this implementation doesn't really do justice +to the power of MMA for nonlinearly constrained problems. + +It is under the same MIT license as the rest of my code in NLopt (see +../COPYRIGHT). + +Steven G. Johnson +July 2008 diff --git a/mma/mma.c b/mma/mma.c new file mode 100644 index 0000000..23f9ee4 --- /dev/null +++ b/mma/mma.c @@ -0,0 +1,107 @@ +#include +#include +#include + +#include "mma.h" + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* magic minimum value for rho in MMA ... the 2002 paper says it should + be a "fixed, strictly positive `small' number, e.g. 1e-5" + ... grrr, I hate these magic numbers, which seem like they + should depend on the objective function in some way */ +#define MMA_RHOMIN 1e-5 + +nlopt_result mma_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, + nlopt_stopping *stop) +{ + nlopt_result ret = NLOPT_SUCCESS; + double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur; + int j, k = 0; + + sigma = (double *) malloc(sizeof(double) * 6*n); + if (!sigma) return NLOPT_OUT_OF_MEMORY; + dfdx = sigma + n; + dfdx_cur = dfdx + n; + xcur = dfdx_cur + n; + xprev = xcur + n; + xprevprev = xprev + n; + for (j = 0; j < n; ++j) + sigma[j] = 0.5 * (ub[j] - lb[j]); + rho = 1.0; + + fcur = *minf = f(n, x, dfdx, f_data); + memcpy(xcur, x, sizeof(double) * n); + stop->nevals++; + while (1) { /* outer iterations */ + double fprev = fcur; + if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; + else if (nlopt_stop_time(stop)) ret = NLOPT_MAXEVAL_REACHED; + else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; + if (ret != NLOPT_SUCCESS) goto done; + if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n); + memcpy(xprev, xcur, sizeof(double) * n); + + while (1) { /* inner iterations */ + double gcur = *minf, w = 0; + for (j = 0; j < n; ++j) { + /* because we have no constraint functions, minimizing + the MMA approximate function can be done analytically */ + double dx = -sigma[j]*sigma[j]*dfdx[j] + / (2*sigma[j]*fabs(dfdx[j]) + rho); + xcur[j] = x[j] + dx; + if (xcur[j] > x[j] + 0.9*sigma[j]) + xcur[j] = x[j] + 0.9*sigma[j]; + else if (xcur[j] < x[j] - 0.9*sigma[j]) + xcur[j] = x[j] - 0.9*sigma[j]; + if (xcur[j] > ub[j]) xcur[j] = ub[j]; + else if (xcur[j] < lb[j]) xcur[j] = lb[j]; + dx = xcur[j] - x[j]; + gcur += (sigma[j]*sigma[j]*dfdx[j]*dx + + sigma[j]*fabs(dfdx[j])*dx*dx + + 0.5*rho*dx*dx) / (sigma[j]*sigma[j]-dx*dx); + w += 0.5*rho*dx*dx / (sigma[j]*sigma[j]-dx*dx); + } + + fcur = f(n, xcur, dfdx_cur, f_data); + stop->nevals++; + if (fcur < *minf) { + *minf = fcur; + memcpy(x, xcur, sizeof(double)*n); + memcpy(dfdx, dfdx_cur, sizeof(double)*n); + } + if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; + else if (nlopt_stop_time(stop)) ret = NLOPT_MAXEVAL_REACHED; + else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; + if (ret != NLOPT_SUCCESS) goto done; + + if (gcur >= fcur) break; + rho = MIN(10*rho, 1.1 * (rho + (fcur - gcur) / w)); + } + + if (nlopt_stop_ftol(stop, fcur, fprev)) + ret = NLOPT_FTOL_REACHED; + if (nlopt_stop_x(stop, xcur, xprev)) + ret = NLOPT_XTOL_REACHED; + if (ret != NLOPT_SUCCESS) goto done; + + /* update rho and sigma for iteration k+1 */ + rho = MAX(0.1 * rho, MMA_RHOMIN); + if (k > 1) + for (j = 0; j < n; ++j) { + double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]); + double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1); + sigma[j] *= gam; + sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j])); + sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j])); + } + } + + done: + free(sigma); + return ret; +} diff --git a/mma/mma.h b/mma/mma.h new file mode 100644 index 0000000..d492de8 --- /dev/null +++ b/mma/mma.h @@ -0,0 +1,23 @@ +#ifndef MMA_H +#define MMA_H + +#include "nlopt.h" +#include "nlopt-util.h" + +#ifdef __cplusplus +extern "C" +{ +#endif /* __cplusplus */ + +nlopt_result mma_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, + nlopt_stopping *stop); + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif + -- 2.30.2