added MMA implementation

darcs-hash:20080729061836-c8de0-7119b5077ac025515b357e848c140f2f081ff831.gz

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@
 

api/Makefile.am

-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

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;
      }

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;
 

configure.ac

@@ -200,6 +200,7 @@ AC_CONFIG_FILES([
    luksan/Makefile
    crs/Makefile
    mlsl/Makefile
+   mma/Makefile
    test/Makefile
 ])
 

mma/Makefile.am

+AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api
+
+noinst_LTLIBRARIES = libmma.la
+libmma_la_SOURCES = mma.c mma.h
+
+EXTRA_DIST = README

mma/README

+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

mma/mma.c

+#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;
+}

mma/mma.h

+#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
+