This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

new methode to initialize nm_simplex


hi all,

currently the nm_simplex minimizer is initialized by giving it a
starting point (vector) and a vector of initial stepsizes. from this
nm_simplex_set constructs the initial simplex.
i had sometimes difficulties with this initializer because it results in
a "very regular" initial simplex. the first trial steps often only
change one parameter. (just as an example. i tried to optimize the
contour of a mirror defined by a spline interpolation of several points.
the current minimizer would start by trying to shift only individual
points. this, of course were not successful steps. it then contracted
for a long time until it finally took off correctly sometimes).

so my proposal (see patch) is to change the call to
gsl_multimin_fminimizer_set() to use a starting vector and a single
value for the initial step size. nm_simplex_set() then contructs a
regular n-pod of size initial_step_size with the starting point at its
center. you can fine-tune the orientation of this initial simplex by
using different values for the environment variable GSL_NM_SIMPLEX_SEED.

notes:
-make test in multimin fails to build because it needs to link
with ../rng/libgslrng.la (you have to edit the makefile to make it work.
i just did not find where to configure this thing)
-because this is a change in the api it may make sense to change it even
more drastic i.e. to only use one generic gsl_multimin_minimizer_set()
(the one currently used for the fdf minimizer) and ignore the parameters
that do not apply for the given minimizer. this would be similar to how
some ode solvers make use of the jacobian and some do not.

cheers
-- 
Dr. Ivo Alxneit
Laboratory for Solar Technology   phone: +41 56 310 4092
Paul Scherrer Institute             fax: +41 56 310 2688
CH-5232 Villigen                   http://solar.web.psi.ch
Switzerland                   gnupg key: 0x515E30C7

diff -Naur /usr/local/src/gsl-1.7/multimin/fminimizer.c multimin/fminimizer.c
--- /usr/local/src/gsl-1.7/multimin/fminimizer.c	2005-06-26 15:25:35.000000000 +0200
+++ multimin/fminimizer.c	2006-02-23 16:16:42.000000000 +0100
@@ -74,14 +74,14 @@
 gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * s,
                              gsl_multimin_function * f,
                              const gsl_vector * x,
-                             const gsl_vector * step_size)
+                             double step_size)
 {
   if (s->x->size != f->n)
     {
       GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN);
     }
   
-  if (x->size != f->n || step_size->size != f->n) 
+  if (x->size != f->n) 
     {
       GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN);
     }  
@@ -89,8 +89,10 @@
   s->f = f;
 
   gsl_vector_memcpy (s->x,x);
-
-  return (s->type->set) (s->state, s->f, s->x, &(s->size), step_size);
+  
+  s->size = step_size;
+  
+  return (s->type->set) (s->state, s->f, s->x, step_size);
 }
 
 void
diff -Naur /usr/local/src/gsl-1.7/multimin/gsl_multimin.h multimin/gsl_multimin.h
--- /usr/local/src/gsl-1.7/multimin/gsl_multimin.h	2005-06-26 15:25:35.000000000 +0200
+++ multimin/gsl_multimin.h	2006-02-23 16:16:33.000000000 +0100
@@ -84,8 +84,7 @@
   int (*alloc) (void *state, size_t n);
   int (*set) (void *state, gsl_multimin_function * f,
               const gsl_vector * x, 
-              double * size,
-              const gsl_vector * step_size);
+              double step_size);
   int (*iterate) (void *state, gsl_multimin_function * f, 
                   gsl_vector * x, 
                   double * size,
@@ -117,7 +116,7 @@
 gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * s,
                              gsl_multimin_function * f,
                              const gsl_vector * x,
-                             const gsl_vector * step_size);
+                             double step_size);
 
 void
 gsl_multimin_fminimizer_free(gsl_multimin_fminimizer *s);
diff -Naur /usr/local/src/gsl-1.7/multimin/simplex.c multimin/simplex.c
--- /usr/local/src/gsl-1.7/multimin/simplex.c	2005-06-26 15:25:35.000000000 +0200
+++ multimin/simplex.c	2006-02-23 16:16:21.000000000 +0100
@@ -33,7 +33,9 @@
 
 #include <config.h>
 #include <stdlib.h>
+#include <math.h>
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_rng.h>
 #include <gsl/gsl_multimin.h>
 
 typedef struct
@@ -45,6 +47,169 @@
 }
 nmsimplex_state_t;
 
+static double a_m(const int m, const double c)
+{
+  return (c * sqrt(1 - c) / sqrt((1 + (m - 1) * c) * (1 + m * c)));
+}
+
+static double b_m(const size_t m, const double c)
+{
+  return (sqrt((1 - c) * (1 + m * c) / (1 + (m - 1) * c)));
+}
+
+static gsl_matrix *init_npod(const size_t n)
+{
+/*
+ * given the following n+1 vectors ei: 0<i<=n in n dimensions
+ * 
+ *    e0(a0, 0, 0, ...)
+ *    e1(a0, b1, 0, ..)
+ *    e2(a0, a1, b2, .)
+ * en-1(a0, ......, bn)
+ *   en(a0, ....., -bn)
+ * 
+ * that form a regular normalized n+1_pod in n dimensions, i.e.
+ * a regular triangle in 2D or a regulart tetrahedron in 3D
+ * 
+ * requesting that 
+ * all (ei, ei) = 1
+ * all (ei, ej) = c
+ */
+
+  size_t i, j;
+  const double c = -1.0 / n;
+  gsl_vector_view v;
+  gsl_matrix *e;
+
+
+  e = gsl_matrix_calloc(n + 1, n);
+  if (e == NULL)
+    {
+      return (NULL);
+    }
+
+/*
+ * init e_0    
+ */
+    gsl_matrix_set(e, 0, 0, 1.0);
+
+/*
+ * init e_1 to e_n-1
+ */
+  for (i = 1; i < n; i++)
+    {
+      for (j = 0; j < i; j++)
+        gsl_matrix_set(e, i, j, a_m((int) j, c));
+
+      gsl_matrix_set(e, i, j, b_m(j, c));
+
+      for (j = i + 1; j < n; j++)
+        gsl_matrix_set(e, i, j, 0.0);
+    }
+
+/*
+ * init e_n
+ */
+  v = gsl_matrix_row(e, n - 1);
+  gsl_matrix_set_row(e, n, &v.vector);
+
+  gsl_matrix_set(e, n, n - 1, -gsl_matrix_get(e, n - 1, n - 1));
+
+  return (e);
+
+}
+
+static gsl_matrix *init_basis(const size_t n)
+{
+/*
+ * orthonolmal basis, randomly oriented
+ */
+  size_t i, j;
+  char *val;
+  unsigned int seed;
+  gsl_rng *r;
+  gsl_matrix *a, *basis;
+
+  a = gsl_matrix_alloc(n, n);
+  if (a == NULL)
+    {
+      return (NULL);
+    }
+
+  basis = gsl_matrix_alloc(n, n);
+  if (basis == NULL)
+    {
+      gsl_matrix_free(a);
+      return (NULL);
+    }
+
+  val = getenv("GSL_NM_SIMPLEX_SEED");
+  if (val)
+      seed = strtoul(val, 0, 0);
+  else
+      seed = 0;
+
+  r = gsl_rng_alloc(gsl_rng_taus);
+  if (r == NULL)
+    {
+      gsl_matrix_free(a);
+      gsl_matrix_free(basis);
+      return (NULL);
+    }
+  gsl_rng_set(r, seed);
+
+
+/*
+ * get n random vectors of length n
+ */
+  for (i = 0; i < n; i++)
+    for (j = 0; j < n; j++)
+      gsl_matrix_set(a, i, j, 2.0 * gsl_rng_uniform(r) - 1.0);
+
+  gsl_rng_free(r);
+
+/*
+ * orthogonalize
+ * 
+ * 	FIXME: we do not check for colinear vectors
+ */
+  gsl_matrix_memcpy(basis, a);
+
+  for (i = 0; i < n; i++)
+    {
+      gsl_vector_view v1 = gsl_matrix_row(a, i);
+      gsl_vector_view v2 = gsl_matrix_row(basis, i);
+
+      for (j = 0; j < i; j++)
+        {
+          double t1, t2;
+          gsl_vector_view v3 = gsl_matrix_row(basis, j);
+
+          gsl_blas_ddot(&v1.vector, &v3.vector, &t1);
+          gsl_blas_ddot(&v3.vector, &v3.vector, &t2);
+
+          gsl_blas_daxpy(-t1 / t2, &v3.vector, &v2.vector);
+
+        }
+    }
+
+  gsl_matrix_free(a);
+
+/*
+ * normalize basis vectors
+ */
+  for (i = 0; i < n; i++)
+    {
+      gsl_vector_view v1 = gsl_matrix_row(basis, i);
+      double t1 = gsl_blas_dnrm2(&v1.vector);
+
+      gsl_vector_scale(&v1.vector, 1 / t1);
+
+    }
+
+    return (basis);
+
+}
 static double
 nmsimplex_move_corner (const double coeff, const nmsimplex_state_t * state,
                        size_t corner, gsl_vector * xc,
@@ -216,47 +381,65 @@
 }
 
 static int
-nmsimplex_set (void *vstate, gsl_multimin_function * f,
-               const gsl_vector * x,
-               double *size, const gsl_vector * step_size)
+nmsimplex_set(void *vstate, gsl_multimin_function * f,
+	      const gsl_vector * x, double step_size)
 {
   int status;
   size_t i;
-  double val;
+  gsl_matrix *n_pod, *basis;
 
   nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
 
-  gsl_vector *xtemp = state->ws1;
+  size_t n = x->size;
+
+/*
+ * prepare regular, normalized n_pod
+ * equilateral triangle 2D
+ * equilateral tetrahedron 3D
+ */
+  n_pod = init_npod(n);
+  if (n_pod == NULL)
+    {
+      GSL_ERROR("failed to allocate space for n_pod", GSL_ENOMEM);
+    }
 
-  /* first point is the original x0 */
+/*
+ * construct a randomly oriented, orthonormal basis in n dimensions
+ */
+  basis = init_basis(n);
+  if (basis == NULL)
+    {
+      gsl_matrix_free(n_pod);
+      GSL_ERROR("failed to allocate space for basis", GSL_ENOMEM);
+    }
 
-  val = GSL_MULTIMIN_FN_EVAL (f, x);
-  gsl_matrix_set_row (state->x1, 0, x);
-  gsl_vector_set (state->y1, 0, val);
+/*
+ * scale n_pod by step_size and express it in basis. result in state->x1
+ */
+  gsl_blas_dgemm(CblasNoTrans, CblasTrans, step_size, n_pod, basis, 0.0,
+		   state->x1);
 
-  /* following points are initialized to x0 + step_size */
+  gsl_matrix_free(basis);
+  gsl_matrix_free(n_pod);
 
-  for (i = 0; i < x->size; i++)
+/*
+ * now state->x1 contains n+1 vectors forming a n+1 dimensional, randomly
+ * oriented simplex of size *size centered at the origin.
+ * we now offset if by the initial vector x and initialize the
+ * corresponding y values in state->y1.
+ */
+  for (i = 0; i <= n; i++)
     {
-      status = gsl_vector_memcpy (xtemp, x);
+      double val;
+      gsl_vector_view v = gsl_matrix_row(state->x1, i);
 
-      if (status != 0)
-        {
-          GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
-        }
+      gsl_vector_add(&v.vector, x);
 
-      val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i);
-      gsl_vector_set (xtemp, i, val);
-      val = GSL_MULTIMIN_FN_EVAL (f, xtemp);
-      gsl_matrix_set_row (state->x1, i + 1, xtemp);
-      gsl_vector_set (state->y1, i + 1, val);
+      val = GSL_MULTIMIN_FN_EVAL(f, &v.vector);
+      gsl_vector_set(state->y1, i, val);
     }
 
-  /* Initialize simplex size */
-
-  *size = nmsimplex_size (state);
-
-  return GSL_SUCCESS;
+    return GSL_SUCCESS;
 }
 
 static void
diff -Naur /usr/local/src/gsl-1.7/multimin/test.c multimin/test.c
--- /usr/local/src/gsl-1.7/multimin/test.c	2006-02-22 16:27:14.000000000 +0100
+++ multimin/test.c	2006-02-24 08:43:43.000000000 +0100
@@ -58,7 +58,7 @@
       test_fdf("Rosenbrock", &rosenbrock, rosenbrock_initpt,*T);
       T++;
     }
-printf("simplex\n");
+
   test_f("Roth", &roth_fmin, roth_initpt);
   test_f("Wood", &wood_fmin, wood_initpt);
   test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt);
@@ -143,15 +143,12 @@
 
   gsl_vector *x = gsl_vector_alloc (f->n);
 
-  gsl_vector *step_size = gsl_vector_alloc (f->n);
+  double step_size = 1.0;
 
   gsl_multimin_fminimizer *s;
 
   (*initpt) (x);
 
-  for (i = 0; i < f->n; i++) 
-    gsl_vector_set (step_size, i, 1);
-
   s = gsl_multimin_fminimizer_alloc(T, f->n);
 
   gsl_multimin_fminimizer_set (s, f, x, step_size);
@@ -185,7 +182,6 @@
 
   gsl_multimin_fminimizer_free(s);
   gsl_vector_free(x);
-  gsl_vector_free(step_size);
 
   return status;
 }

Attachment: signature.asc
Description: This is a digitally signed message part


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]