chiark / gitweb /
added MMA implementation
authorstevenj <stevenj@alum.mit.edu>
Tue, 29 Jul 2008 06:18:36 +0000 (02:18 -0400)
committerstevenj <stevenj@alum.mit.edu>
Tue, 29 Jul 2008 06:18:36 +0000 (02:18 -0400)
darcs-hash:20080729061836-c8de0-7119b5077ac025515b357e848c140f2f081ff831.gz

Makefile.am
api/Makefile.am
api/nlopt.c
api/nlopt.h
configure.ac
mma/Makefile.am [new file with mode: 0644]
mma/README [new file with mode: 0644]
mma/mma.c [new file with mode: 0644]
mma/mma.h [new file with mode: 0644]

index 1a3d1b457a3db2f9da7acda9109414971bf2be14..3570f88ad527ee9ce28c9461ac933aa5d62cb48b 100644 (file)
@@ -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@
 
index 0823aaa40ad767c6bc2604a0c058666c60d9a392..3a616b1fdf042eb202765bfc35fc6fb56f4c48a2 100644 (file)
@@ -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
index a39010c3fdd52856b546e3dc7ce2b3987e9eb062..6ad36f38cb9ba8811ec60c9c436ccb69e14fbc70 100644 (file)
@@ -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;
      }
index b6e19a5bd8e4ddf3dc9ae3ef0cb7b58a490c3070..edd04c41233f88f9d0bc92a45122765d0eac4a09 100644 (file)
@@ -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;
 
index 8a3b8a196a18c2f0f69b7410657d704800ef7d38..4f98d81e0054656e229e782de3adfa54323f22ba 100644 (file)
@@ -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 (file)
index 0000000..9d4811e
--- /dev/null
@@ -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 (file)
index 0000000..07fa109
--- /dev/null
@@ -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 (file)
index 0000000..23f9ee4
--- /dev/null
+++ b/mma/mma.c
@@ -0,0 +1,107 @@
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#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 (file)
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
+