new methode to initialize nm_simplex

new methode to initialize nm_simplex

[Ivo Alxneit-Kamber]

[-- Attachment #1.1: Type: text/plain, Size: 1909 bytes --]

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


[-- Attachment #1.2: patch --]
[-- Type: text/x-patch, Size: 9650 bytes --]

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 #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 189 bytes --]

Re: new methode to initialize nm_simplex

[Brian Gough]

Ivo Alxneit-Kamber writes:
 > 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.

If you take the existing step_size vector, you could create a n-pod
oriented in that direction with size |step_size| -- that would avoid
any changes to the API, and allow it to be a new method (reusing most
of the existing functions).

To randomize the other orthogonal directions I'd suggest putting a
simple RNG like 'vax' inside the method itself, so there aren't any
dependencies on libgslrng.

-- 
Brian Gough

Re: new methode to initialize nm_simplex

[Ivo Alxneit-Kamber]

[-- Attachment #1: Type: text/plain, Size: 3951 bytes --]

On Thu, 2006-03-02 at 19:43 +0100, Brian Gough wrote:
> Ivo Alxneit-Kamber writes:
>  > 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.
> 
> If you take the existing step_size vector, you could create a n-pod
> oriented in that direction with size |step_size| -- that would avoid
> any changes to the API, and allow it to be a new method (reusing most
> of the existing functions).
i fear you missunderstood (maybe my comments were not really sufficient)
or i do not really get your point here.

i think i better explain how the initialization is supposed to work with
an example for a two parameter case (2D case, the simplex is a
triangle). the current linitialization gives you a triangle with:

s0=(x0, y0) 	-> starting point, i.e. vector x (x0,y0)
s1=(x0 + a, y0)	-> from start vector x
s2=(x0, y0 + b)	   and vector step_size (a,b)

so step_size defines

|s1-s0|=a and |s2-s0|=b i.e. the initial step size in the direction of
the two parameters. this allows that x_min and y_min are of largely
different magnitude and  that the initial steps the algorithm takes
honor this by e.g. taking large steps in one direction and small ones in
the other.

now the new initialization i propose will give you an equilateral
triangle its center at (x0,y0) (none of the vertices is at the center)
that is randomly oriented as initial simplex. here i only request that

|s0|=|s1|=s2|=1 and (s0,s1)=(s0,s2)=(s1,s2)=constant.

because the triangle is randomly oriented you do not know how how to
scale the individaul vertices to force the algorithm to take large steps
in one and small steps in an other direction. so i think the
initialization only makes sense with an equilateral triangle.

\begin{parentheses}
i think it is also generally accepted good practice to scale the
function to be minimized in a way to ensure that all the parameters are
of the same magnitude. so defining individual stepsizes is not really
necessary.

if you think breaking the existing api is a-bad-thing-to-do i propose
the following:
- if all elements of the vector step_size are different from zero
  the old (current) methode is used.
- if only the first element is different from zero the new (proposed)
  method is used.
this way old code does not break.
\end{parentheses}


> To randomize the other orthogonal directions I'd suggest putting a
> simple RNG like 'vax' inside the method itself, so there aren't any
> dependencies on libgslrng.
> 
would be using the standart c random number generator in stdlib.h
srandom(seed) to initialize and random() to pull a random number be also
ok?
thats fine with me. i will make the changes.

p.s. i will be off for skiing the next week, so do not expect
answers/comments/code right now.
-- 
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


[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 189 bytes --]

Re: new methode to initialize nm_simplex

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

-- 
Brian Gough

Re: new methode to initialize nm_simplex

[Alxneit-Kamber Ivo]

[-- Attachment #1.1: Type: text/plain, Size: 1016 bytes --]

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


[-- Attachment #1.2: simplex.patch --]
[-- Type: text/x-patch, Size: 8834 bytes --]

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

[-- Attachment #1.3: test.patch --]
[-- Type: text/x-patch, Size: 1986 bytes --]

--- /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);

[-- Attachment #1.4: doc.patch --]
[-- Type: text/x-patch, Size: 1224 bytes --]

--- /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 #2: Dies ist ein digital signierter Nachrichtenteil --]
[-- Type: application/pgp-signature, Size: 189 bytes --]