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]

Re: new methode to initialize nm_simplex


Am Freitag, den 17.03.2006, 16:13 +0000 schrieb Brian Gough:
> Ivo Alxneit-Kamber writes:
>  > i fear you missunderstood (maybe my comments were not really sufficient)
>  > or i do not really get your point here.
>  > 
> 
> Just that you don't need to modify the API if your algorithm only
> needs an initial scalar step instead of a vector -- use the magnitude
> of the vector (and ignore the direction if necessary, or put it to
> some use).  That fits the existing API.
> 
ok, done so (plus your previous remarks)
- supplied own (vax) random number generator.
- API kept (old initializer still supported).
- updated test case. both initializer are tested now.
- updated documentation (brian, could you please have a look at the
  language).

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

--- /tmp/gsl-1.7/multimin/simplex.c	2005-06-26 15:25:35.000000000 +0200
+++ gsl-1.7/multimin/simplex.c	2006-03-20 08:51:58.000000000 +0100
@@ -32,7 +32,7 @@
 */
 
 #include <config.h>
-#include <stdlib.h>
+//#include <stdlib.h>
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_multimin.h>
 
@@ -45,6 +45,287 @@
 }
 nmsimplex_state_t;
 
+
+/* The following is a copy of gsl_rng_type vax. */
+
+typedef struct
+  {
+    unsigned long int x;
+  }
+my_rng_state_t;
+
+static void
+my_rng_set(my_rng_state_t *s, const unsigned long int seed)
+{
+ s->x = seed;
+}
+
+static double
+my_rng_get(my_rng_state_t *s)
+{
+  s->x = (69069 * s->x + 1) & 0xffffffffUL;
+  return (s->x / 4294967296.0);
+}
+
+
+static int
+old_init(void *vstate, gsl_multimin_function * f,
+         const gsl_vector * x, const gsl_vector * step_size)
+{
+  int status;
+  size_t i;
+  double val;
+  gsl_vector *xtemp;
+
+  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
+
+  xtemp = state->ws1;
+
+  /* first point is the original x0 */
+
+  val = GSL_MULTIMIN_FN_EVAL (f, x);
+  gsl_matrix_set_row (state->x1, 0, x);
+  gsl_vector_set (state->y1, 0, val);
+
+  /* following points are initialized to x0 + step_size */
+
+  for (i = 0; i < x->size; i++)
+    {
+      status = gsl_vector_memcpy (xtemp, x);
+
+      if (status != 0)
+	{
+	  GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
+	}
+
+      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);
+    }
+  return GSL_SUCCESS;
+}
+
+
+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 e_i: 0<i<=n in n dimensions
+
+       e_0(a0, 0, 0, ...)
+       e_1(a0, b1, 0, ..)
+       e_2(a0, a1, b2, .)
+    e_n-1(a0, ......, bn)
+      e_n(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 (e_i, e_i) = 1
+       all (e_i, e_j) = 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);
+    }
+
+  /* Initialize e_0. */
+
+    gsl_matrix_set(e, 0, 0, 1.0);
+
+  /* Initialize 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);
+    }
+
+  /* Initialize 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;
+  my_rng_state_t *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);
+    }
+
+  r = (my_rng_state_t *)malloc(sizeof(my_rng_state_t));
+  if (r == NULL)
+    {
+      gsl_matrix_free(basis);
+      gsl_matrix_free(a);
+      return (NULL);
+    }
+  
+  val = getenv("GSL_NM_SIMPLEX_SEED");
+
+  if (val)
+      seed = strtoul(val, 0, 0);
+  else
+      seed = 0;
+
+  my_rng_set(r, seed); 
+
+
+  /* Get n random vectors of length n. */
+  
+  for (i = 0; i < n; i++)
+    for (j = 0; j < n; j++)
+      {
+        double value = my_rng_get(r);
+        gsl_matrix_set(a, i, j, 2.0 * value - 1.0);
+      }
+
+  free(r);
+  
+  /* Orthogonalize matrix a, result in matrix basis.
+	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 int
+new_init(void *vstate, gsl_multimin_function * f,
+         const gsl_vector * x, const double step_size)
+{
+  size_t i;
+  gsl_matrix *n_pod, *basis;
+
+  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
+
+  size_t n = x->size;
+
+  n_pod = init_npod(n);
+  if (n_pod == NULL)
+    {
+      GSL_ERROR("failed to allocate space for n_pod", GSL_ENOMEM);
+    }
+
+  /* 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);
+    }
+
+  /* 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);
+
+  gsl_matrix_free(basis);
+  gsl_matrix_free(n_pod);
+
+  /* Now state->x1 contains n+1 vectors forming a n+1 dimensional, randomly
+     oriented simplex of size gsl_vector_get(step_size, 0) centered at the
+     origin. we now offset it by the initial vector x and initialize the
+     corresponding y values in state->y1. */
+     
+  for (i = 0; i <= n; i++)
+    {
+      double val;
+      gsl_vector_view v = gsl_matrix_row(state->x1, i);
+
+      gsl_vector_add(&v.vector, x);
+
+      val = GSL_MULTIMIN_FN_EVAL(f, &v.vector);
+      gsl_vector_set(state->y1, i, val);
+    }
+
+  return GSL_SUCCESS;
+
+}
+
+
 static double
 nmsimplex_move_corner (const double coeff, const nmsimplex_state_t * state,
                        size_t corner, gsl_vector * xc,
@@ -221,35 +502,53 @@
                double *size, const gsl_vector * step_size)
 {
   int status;
-  size_t i;
-  double val;
+  size_t i, n_non_zero;
 
   nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
 
-  gsl_vector *xtemp = state->ws1;
-
-  /* first point is the original x0 */
-
-  val = GSL_MULTIMIN_FN_EVAL (f, x);
-  gsl_matrix_set_row (state->x1, 0, x);
-  gsl_vector_set (state->y1, 0, val);
-
-  /* following points are initialized to x0 + step_size */
+  size_t n = x->size;
+  
+  /* Find out which method to use for setting the initial simplex.
+     If only one element of step_size is non-zero we set the initial
+     simplex to a randomly oriented n-pod. */
+     
+  n_non_zero = 0;
+  for(i=0; i<n; i++)
+    {
+      if(gsl_vector_get(step_size, i) != 0.0)
+        {
+          n_non_zero++;
+	}
+    }
 
-  for (i = 0; i < x->size; i++)
+  if(n_non_zero == n)
     {
-      status = gsl_vector_memcpy (xtemp, x);
+      /* Use old methode. */
+      
+      status = old_init(vstate, f, x, step_size);
+      if (status != GSL_SUCCESS)
+        {
+	  return status;
+	}
 
-      if (status != 0)
+    }
+  else if(n_non_zero == 1)
+    {
+      /* New initialization method.
+         Prepare regular, normalized n_pod:
+             equilateral triangle 2D
+             equilateral tetrahedron 3D. */
+	     
+      status = new_init(vstate, f, x, gsl_vector_get(step_size, 0));
+      if (status != GSL_SUCCESS)
         {
-          GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
-        }
+	  return status;
+	}
 
-      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);
+    }
+  else
+    {
+      GSL_ERROR("invalid step size", GSL_EINVAL);
     }
 
   /* Initialize simplex size */
--- /tmp/gsl-1.7/multimin/test.c	2005-06-26 15:25:35.000000000 +0200
+++ gsl-1.7/multimin/test.c	2006-03-18 19:20:24.000000000 +0100
@@ -33,7 +33,8 @@
          initpt_function initpt, const gsl_multimin_fdfminimizer_type *T);
 
 int
-test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt);
+test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt,
+       const int method);
 
 int
 main (void)
@@ -59,9 +60,15 @@
       T++;
     }
 
-  test_f("Roth", &roth_fmin, roth_initpt);
-  test_f("Wood", &wood_fmin, wood_initpt);
-  test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt);
+/* test old initialization method */
+  test_f("Roth", &roth_fmin, roth_initpt, 0);
+  test_f("Wood", &wood_fmin, wood_initpt, 0);
+  test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt, 0);
+
+/* test new initialization method */
+  test_f("Roth", &roth_fmin, roth_initpt, 1);
+  test_f("Wood", &wood_fmin, wood_initpt, 1);
+  test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt, 1);
 
   T = fdfminimizers;
 
@@ -133,7 +140,8 @@
 }
 
 int
-test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt)
+test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt,
+       const int method)
 {
   /* currently this function tests only nmsimplex */
 
@@ -143,15 +151,23 @@
 
   gsl_vector *x = gsl_vector_alloc (f->n);
 
-  gsl_vector *step_size = gsl_vector_alloc (f->n);
+  gsl_vector *step_size = gsl_vector_calloc (f->n);
 
   gsl_multimin_fminimizer *s;
 
   (*initpt) (x);
 
-  for (i = 0; i < f->n; i++) 
-    gsl_vector_set (step_size, i, 1);
-
+  if (method == 0)
+    {
+      for (i = 0; i < f->n; i++) 
+        gsl_vector_set (step_size, i, 1);
+    }
+    
+  if (method == 1)
+    {
+      gsl_vector_set(step_size, 0, 1);
+    }
+  
   s = gsl_multimin_fminimizer_alloc(T, f->n);
 
   gsl_multimin_fminimizer_set (s, f, x, step_size);
--- /tmp/gsl-1.7/doc/multimin.texi	2005-05-21 15:28:05.000000000 +0200
+++ gsl-1.7/doc/multimin.texi	2006-03-20 09:19:27.000000000 +0100
@@ -457,7 +457,18 @@
 @end ifinfo
 @noindent
 These vectors form the @math{n+1} vertices of a simplex in @math{n}
-dimensions.  On each iteration the algorithm tries to improve
+dimensions.  Alternatively, if only the first component of @var{step_size}
+is non zero, a randomly oriented regular @var{n+1}-pod of the size given
+in @var{step_size}, centered at @var{x} is used as simplex to start
+with.  The orientation of this initial simplex can be controlled by
+setting the environment variable @code{GSL_NM_SIMPLEX_SEED} to some
+integer value.  The initial simplex then corresponds to
+an equilateral triangle (two dimensional problem) or to
+a regular tetrahedron (three dimensional problem).  Note, that
+the parameters to be determined should be are of similar magnitude
+if this initialize is applied.
+ 
+On each iteration the algorithm tries to improve
 the parameter vector @math{p_i} corresponding to the highest
 function value by simple geometrical transformations.  These
 are reflection, reflection followed by expansion, contraction and multiple

Attachment: signature.asc
Description: Dies ist ein digital signierter Nachrichtenteil


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