added my own implementation of controlled random search (CRS) algorithm

darcs-hash:20070907015248-c8de0-5f583ad3dac1a4d688b992a894c8ce75745bbbc7.gz

Makefile.am

@@ -3,13 +3,13 @@ lib_LTLIBRARIES = libnlopt.la
 
 ACLOCAL_AMFLAGS=-I ./m4
 
-SUBDIRS= util subplex direct cdirect stogo praxis luksan api octave . test
+SUBDIRS= util subplex direct cdirect stogo praxis luksan crs api octave . test
 EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4
 
 libnlopt_la_SOURCES = 
 libnlopt_la_LIBADD = subplex/libsubplex.la direct/libdirect.la	\
 cdirect/libcdirect.la stogo/libstogo.la praxis/libpraxis.la	\
-luksan/libluksan.la api/libapi.la util/libutil.la
+luksan/libluksan.la crs/libcrs.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)/luksan -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)/luksan -I$(top_srcdir)/crs -I$(top_srcdir)/util
 
 include_HEADERS = nlopt.h
 noinst_LTLIBRARIES = libapi.la

api/nlopt.c

@@ -41,7 +41,7 @@ void nlopt_version(int *major, int *minor, int *bugfix)
 
 /*************************************************************************/
 
-static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = {
+static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = {
      "DIRECT (global, no-derivative)",
      "DIRECT-L (global, no-derivative)",
      "Randomized DIRECT-L (global, no-derivative)",
@@ -61,6 +61,9 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = {
      "Truncated Newton with restarting (local, derivative-based)",
      "Preconditioned truncated Newton (local, derivative-based)",
      "Preconditioned truncated Newton with restarting (local, derivative-based)",
+     "Controlled random search (CRS2) with local mutation (global, no-derivative)",
+     "Controlled quasi-random search (CRS2) with local mutation (global, no-derivative), using Sobol' LDS",
+     
 };
 
 const char *nlopt_algorithm_name(nlopt_algorithm a)
@@ -125,6 +128,8 @@ static double f_direct(int n, const double *x, int *undefined, void *data_)
 
 #include "luksan.h"
 
+#include "crs.h"
+
 /*************************************************************************/
 
 /* for "hybrid" algorithms that combine global search with some
@@ -323,6 +328,9 @@ static nlopt_result nlopt_minimize_(
 				 1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
 				 1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
 
+	 case NLOPT_GN_CRS2_LM:
+	      return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, 0);
+
 	 default:
 	      return NLOPT_INVALID_ARGS;
      }

api/nlopt.h

@@ -50,6 +50,8 @@ typedef enum {
      NLOPT_LD_TNEWTON_PRECOND,
      NLOPT_LD_TNEWTON_PRECOND_RESTART,
 
+     NLOPT_GN_CRS2_LM,
+
      NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */
 } nlopt_algorithm;
 

configure.ac

@@ -181,6 +181,7 @@ AC_CONFIG_FILES([
    subplex/Makefile
    praxis/Makefile
    luksan/Makefile
+   crs/Makefile
    test/Makefile
 ])
 

crs/Makefile.am

+AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api
+
+noinst_LTLIBRARIES = libcrs.la
+libcrs_la_SOURCES = crs.c crs.h
+
+EXTRA_DIST = README

crs/README

+This is my implementation of the "controlled random search" (CRS2) algorithm
+with the "local mutation" modification, as defined by:
+
+       P. Kaelo and M. M. Ali, "Some variants of the controlled random
+       search algorithm for global optimization," J. Optim. Theory Appl.
+       130 (2), 253-264 (2006).
+
+It is under the same MIT license as the rest of my code in NLopt (see
+../COPYRIGHT).
+
+Steven G. Johnson
+September 2007

crs/crs.c

+#include <stdlib.h>
+#include <string.h>
+
+#include "crs.h"
+#include "redblack.h"
+
+/* Controlled Random Search 2 (CRS2) with "local mutation", as defined
+   by:
+       P. Kaelo and M. M. Ali, "Some variants of the controlled random
+       search algorithm for global optimization," J. Optim. Theory Appl.
+       130 (2), 253-264 (2006).
+*/
+
+typedef struct {
+     int n; /* # dimensions */
+     const double *lb, *ub;
+     nlopt_stopping *stop; /* stopping criteria */
+     nlopt_func f; void *f_data;
+
+     int N; /* # points in population */
+     double *ps; /* population array N x (n+1) of tuples [f(x), x] */
+     double *p; /* single point array (length n+1), for temp use */
+     rb_tree t; /* red-black tree of population, sorted by f(x) */
+     nlopt_sobol s; /* sobol data for LDS point generation, or NULL
+		       to use pseudo-random numbers */
+} crs_data;
+
+/* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
+static int crs_compare(double *k1, double *k2)
+{
+     if (*k1 < *k2) return -1;
+     if (*k1 > *k2) return +1;
+     return k1 - k2; /* tie-breaker */
+}
+
+/* set x to a random trial value, as defined by CRS:
+     x = 2G - x_n, 
+   where x_0 ... x_n are distinct points in the population
+   with x_0 the current best point and the other points are random,
+   and G is the centroid of x_0...x_{n-1} */
+static void random_trial(crs_data *d, double *x, rb_node *best)
+{
+     int n = d->n, n1 = n+1, i, k, i0, jn;
+     double *ps = d->ps, *xi;
+
+     /* initialize x to x_0 = best point */
+     memcpy(x, best->k + 1, sizeof(double) * n);
+     i0 = (best->k - ps) / n1;
+
+     jn = nlopt_iurand(n); /* which of remaining n points is "x_n",
+			      i.e. which to reflect through ...
+			      this is necessary since we generate
+			      the remaining points in order, so
+			      just picking the last point would not
+			      be very random */
+
+     /* use "method A" from
+	
+           Jeffrey Scott Vitter, "An efficient algorithm for
+	   sequential random sampling," ACM Trans. Math. Soft. 13 (1),
+	   58--67 (1987).  
+
+        to randomly pick n distinct points out of the remaining N-1 (not 
+        including i0!).  (The same as "method S" in Knuth vol. 2.)
+        This method requires O(N) time, which is fine in our case
+        (there are better methods if n << N). */
+     {
+	  int Nleft = d->N - 1, nleft = n;
+	  int Nfree = Nleft - nleft;
+	  i = 0; i += i == i0;
+	  while (nleft > 1) {
+	       double q = ((double) Nfree) / Nleft;
+	       double v = nlopt_urand(0., 1.);
+	       while (q > v) {
+		    ++i; i += i == i0;
+		    --Nfree; --Nleft;
+		    q = (q * Nfree) / Nleft;
+	       }
+	       xi = ps + n1 * i + 1;
+	       if (jn-- == 0) /* point to reflect through */
+		    for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
+	       else /* point to include in centroid */
+		    for (k = 0; k < n; ++k) x[k] += xi[k];
+	       ++i; i += i == i0;
+	       --Nleft; --nleft;
+	  }
+	  i += nlopt_iurand(Nleft); i += i == i0;
+	  xi = ps + n1 * i + 1;
+	  if (jn-- == 0) /* point to reflect through */
+	       for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
+	  else /* point to include in centroid */
+	       for (k = 0; k < n; ++k) x[k] += xi[k];
+     }
+     for (k = 0; k < n; ++k) {
+	  x[k] *= 2.0 / n; /* renormalize */
+	  if (x[k] > d->ub[k]) x[k] = d->ub[k];
+	  else if (x[k] < d->lb[k]) x[k] = d->lb[k];
+     }
+}
+
+#define NUM_MUTATION 1 /* # "local mutation" steps to try if trial fails */
+
+static nlopt_result crs_trial(crs_data *d)
+{
+     rb_node *best = rb_tree_min(&d->t);
+     rb_node *worst = rb_tree_max(&d->t);
+     int mutation = NUM_MUTATION;
+     int i, n = d->n;
+     random_trial(d, d->p + 1, best);
+     do {
+     	  d->p[0] = d->f(n, d->p + 1, NULL, d->f_data);
+	  d->stop->nevals++;
+	  if (d->p[0] < worst->k[0]) break;
+	  if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED;
+	  if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED;
+	  if (mutation) {
+	       for (i = 0; i < n; ++i) {
+		    double w = nlopt_urand(0.,1.);
+		    d->p[1+i] = best->k[1+i] * (1 + w) - w * d->p[1+i];
+		    if (d->p[1+i] > d->ub[i]) d->p[1+i] = d->ub[i];
+		    else if (d->p[1+i] < d->lb[i]) d->p[1+i] = d->lb[i];
+	       }
+	       mutation--;
+	  }
+	  else {
+	       random_trial(d, d->p + 1, best);
+	       mutation = NUM_MUTATION;
+	  }
+     } while (1);
+     memcpy(worst->k, d->p, sizeof(double) * (n+1));
+     rb_tree_resort(&d->t, worst);
+     return NLOPT_SUCCESS;
+}
+
+static void crs_destroy(crs_data *d)
+{
+     nlopt_sobol_destroy(d->s);
+     rb_tree_destroy(&d->t);
+     free(d->ps);
+}
+
+static nlopt_result crs_init(crs_data *d, int n, const double *x,
+			     const double *lb, const double *ub,
+			     nlopt_stopping *stop, nlopt_func f, void *f_data,
+			     int lds)
+{
+     int i;
+
+     /* TODO: how should we set the initial population size? 
+	the Kaelo and Ali paper suggests 10*(n+1), but should
+	we add more random points if maxeval is large, or... ? */
+     d->N = 10 * (n + 1); /* heuristic initial population size */
+
+     d->n = n;
+     d->stop = stop;
+     d->f = f; d->f_data = f_data;
+     d->ub = ub; d->lb = lb;
+     d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1));
+     if (!d->ps) return NLOPT_OUT_OF_MEMORY;
+     d->p = d->ps + d->N * (n+1);
+     rb_tree_init(&d->t, crs_compare);
+
+     /* we can either use pseudorandom points, as in the original CRS
+	algorithm, or use a low-discrepancy Sobol' sequence ... I tried
+	the latter, however, and it doesn't seem to help, probably
+	because we are only generating a small number of random points
+	to start with */
+     d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
+     nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
+
+     /* generate initial points randomly, plus starting guess x */
+     memcpy(d->ps + 1, x, sizeof(double) * n);
+     d->ps[0] = f(n, x, NULL, f_data);
+     stop->nevals++;
+     rb_tree_insert(&d->t, d->ps);
+     if (d->ps[0] < 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;
+     for (i = 1; i < d->N; ++i) {
+	  double *k = d->ps + i*(n+1);
+	  if (d->s) 
+	       nlopt_sobol_next(d->s, k + 1, lb, ub);
+	  else {
+	       int j;
+	       for (j = 0; j < n; ++j) 
+		    k[1 + j] = nlopt_urand(lb[j], ub[j]);
+	  }
+	  k[0] = f(n, k + 1, NULL, f_data);
+	  stop->nevals++;
+	  rb_tree_insert(&d->t, k);
+	  if (k[0] < 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;	  
+     }
+
+     return NLOPT_SUCCESS;;
+}
+
+nlopt_result crs_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,
+			  int lds) /* random or low-discrepancy seq. (lds) */
+{
+     nlopt_result ret;
+     crs_data d;
+     rb_node *best;
+
+     ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, lds);
+     if (ret < 0) return ret;
+     
+     best = rb_tree_min(&d.t);
+     *minf = best->k[0];
+     memcpy(x, best->k + 1, sizeof(double) * n);
+
+     while (ret == NLOPT_SUCCESS) {
+	  if (NLOPT_SUCCESS == (ret = crs_trial(&d))) {
+	       best = rb_tree_min(&d.t);
+	       if (best->k[0] < *minf) {
+		    if (best->k[0] < stop->minf_max)
+			 ret = NLOPT_MINF_MAX_REACHED;
+		    else if (nlopt_stop_f(stop, best->k[0], *minf))
+			 ret = NLOPT_FTOL_REACHED;
+		    else if (nlopt_stop_x(stop, best->k + 1, x))
+			 ret = NLOPT_XTOL_REACHED;
+		    *minf = best->k[0];
+		    memcpy(x, best->k + 1, sizeof(double) * n);
+	       }
+	       if (ret != NLOPT_SUCCESS) {
+		    if (nlopt_stop_evals(stop)) 
+			 ret = NLOPT_MAXEVAL_REACHED;
+		    else if (nlopt_stop_time(stop)) 
+			 ret = NLOPT_MAXTIME_REACHED;
+	       }
+	  }
+     }
+     crs_destroy(&d);
+     return ret;
+}

util/sobolseq.c

@@ -207,7 +207,7 @@ extern void nlopt_sobol_destroy(nlopt_sobol s)
 /* next vector x[sdim] in Sobol sequence, with each x[i] in (0,1) */
 void nlopt_sobol_next01(nlopt_sobol s, double *x)
 {
-     if (!s || !sobol_gen(s, x)) {
+     if (!sobol_gen(s, x)) {
 	  /* fall back on pseudo random numbers in the unlikely event
 	     that we exceed 2^32-1 points */
 	  unsigned i;
@@ -232,9 +232,11 @@ void nlopt_sobol_next(nlopt_sobol s, double *x,
    points equal to the largest power of 2 smaller than n */
 void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
 {
-     unsigned k = 1;
-     while (k*2 < n) k *= 2;
-     while (k-- > 0) sobol_gen(s, x);
+     if (s) {
+	  unsigned k = 1;
+	  while (k*2 < n) k *= 2;
+	  while (k-- > 0) sobol_gen(s, x);
+     }
 }
 
 /************************************************************************/