simplex minimization

simplex minimization

[Tim F]

I've been using gsl (1.3) to minimize a seven parameter function against 
some sample data I have, and it seems to give me good results based on fdf 
minimization (any of the fr, pr or bfgs algorithms yield similar 
numbers):

iteration: 6
params:
 -0.00131 gradient: -5610116.24709
  0.02920 gradient: -8297325.73283
 -0.05581 gradient: -184748.86988
 -0.00000 gradient:   0.00000
  0.00000 gradient:   0.00000
  0.00000 gradient:   0.00000
  0.20521 gradient: 48147814.97896
f(x): 1491335.000000 tot. gradient: 34055624.129164

these numbers are around what I expect (based on comparison with other
programs that carry out a similar form of minimization) even though the
gradients may seem large (the 4th-6th parameters are forced to zero in the
example I'm testing).  Out of curiousity, I tried to compare these values
with those determined by the simplex algorithm.  However, the minimization
always stops and claims to have succeeded before performing any iterations
(with a size stopping point of 0.01):

iteration: 1
params:
  0.00000
  0.00000
  0.00000
  0.00000
  0.00000
  0.00000
  0.20521
 f(x): 1566284.000000 tot. size: 0.000000000000

and no matter the starting values, this always seems the case - i.e. the
simplex size is always 0.00.  Might it have something to do with forcing
some of the parameters in the minimization to zero?  If code would be
helpful in nailing this problem down, let me know.

	Regards,
	Tim F

Re: simplex minimization

[Brian Gough]

Tim F writes:
 > and no matter the starting values, this always seems the case - i.e. the
 > simplex size is always 0.00.  Might it have something to do with forcing
 > some of the parameters in the minimization to zero?  If code would be
 > helpful in nailing this problem down, let me know.

Yes please send an example.  Thanks.

Re: simplex minimization

[Tim F]

On Wed, 1 Jan 2003, Brian Gough wrote:

>
> Yes please send an example.  Thanks.
> 

I hope this example isn't too terse... let me know if other information 
will help.

	Tim F

#include <stdio.h>
#include <stdlib.h>

#include <glib.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_vector.h>

// num iterations for solver (max)
#define NUM_ITER 100

struct func_data{
  gint low, high;
  gdouble *gstar;
  reflection_p data;
};

/*
 * perform the following fit of k and [M]:
 * Fo = k*Fc*e^-([T]^T[M][T])
 * (T = 2pi[hkl])
 */
void calc_scale_factor(void){
  gint i, j, status;
  // number of parameters
  const size_t p = 7;
  gdouble param_init[7];
  gsl_vector_view param;
  struct func_data d;

  // f only minimizer
  gdouble ssval;
  const gsl_multimin_fminimizer_type *T;
  gsl_multimin_fminimizer *s;
  gsl_multimin_function f;
  gsl_vector *ss;
 
  // our initial M is just 0
  // order: u11 u22 u33 u12 u13 u23
  for (i=0; i<6; i++)
    param_init[i] = 0.0;

  // initial k (<FoFc> / <Fc^2>)
  param_init[6] = initial_scale_factor(curdata->sfo_hkls);

  param = gsl_vector_view_array(param_init, p);

  // set up data for solver
  /*
   * gstar is a 3x3 symmetric tensor
   * sfo_hkls is the Fo and Fc data
   */
  d.low = curdata->low;
  d.high = curdata->high;
  d.gstar = curdata->gstar;
  d.data = curdata->sfo_hkls;

  // set up function
  f.f = &scale_f;
  f.n = p;
  f.params = (void *)&d;

  // set up solver
  T = gsl_multimin_fminimizer_nmsimplex;
  s = gsl_multimin_fminimizer_alloc(T, p);
  ss = gsl_vector_alloc(p);
  gsl_multimin_fminimizer_set(s, &f, &param.vector, ss);

  for (i=1; i<NUM_ITER; i++){
    status = gsl_multimin_fminimizer_iterate(s);
    if (status){
      g_print("GSL minimizer: %s\n", gsl_strerror(status));
      break;
    }

    status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), 1e-2);
    ssval = gsl_multimin_fminimizer_size(s);
    if (status == GSL_SUCCESS)
      g_print("minimum found at:\n");

    g_print("iteration: %d\nparams:\n", i);
    for (j=0; j<p; j++)
      g_print("%9.5f\n", gsl_vector_get(s->x, j));
    g_print(" f(x): %f tot. size: %.12f\n\n", s->fval, ssval);

    if (status != GSL_CONTINUE)
      break;
  }

  // set values based on minimization
  ...
  // free up minimizer
  gsl_multimin_fminimizer_free(s);
  gsl_vector_free(ss);
}

// function
gdouble scale_f(const gsl_vector *x, void *params){
  gint i;
  gint low = ((struct func_data *)params)->low;
  gint high = ((struct func_data *)params)->high;
  gdouble *gstar = ((struct func_data *)params)->gstar;
  reflection_p data = ((struct func_data *)params)->data;
  gdouble m[] = {gsl_vector_get(x, 0), gsl_vector_get(x, 1), 
                 gsl_vector_get(x, 2), gsl_vector_get(x, 3), 
                 gsl_vector_get(x, 4), gsl_vector_get(x, 5)};
  gdouble scale_k = gsl_vector_get(x, 6);

  gdouble u;
  guint64 sum;

  for (i=low, sum=0.0; i<=high; i++){
    if (!data[i].mate)
      continue;

    /*
     * in some cases, the off diagonal elements of gstar (elements 3-5)
     * may be zero, forcing the m[3] to m[5] parameters to zero
     * during the minimization
     */
    // calculate [T]^T[M][T]
    u = gsl_pow_2(data[i].hkl[0]) * gstar[0] * m[0] +
      gsl_pow_2(data[i].hkl[1]) * gstar[1] * m[1] +
      gsl_pow_2(data[i].hkl[2]) * gstar[2] * m[2] +
      2.0 * data[i].hkl[0] * data[i].hkl[1] * gstar[3] * m[3] +
      2.0 * data[i].hkl[0] * data[i].hkl[2] * gstar[4] * m[4] +
      2.0 * data[i].hkl[1] * data[i].hkl[2] * gstar[5] * m[5];

    // calculate sum: (Fo-k*Fc*e^(-[T]^T[M][T]))^2 / (sigFo)^2
    sum += (guint64) ((gsl_pow_2(fabs(data[i].data.fo[0]) - 
                       fabs(gsl_sf_exp_mult(-2.0*M_PI*M_PI*u, 
                       scale_k * gsl_complex_abs(data[i].mate->data.fcphi))))) / 
                       gsl_pow_2(data[i].data.fo[1]));
  }

  return((gdouble) sum);
}

Re: simplex minimization

[ivo.alxneit]

[-- Attachment #1: Type: TEXT/plain, Size: 999 bytes --]

On  1 Jan, Tim F wrote:
> On Wed, 1 Jan 2003, Brian Gough wrote:
> 
>>
>> Yes please send an example.  Thanks.
>> 
> 
> I hope this example isn't too terse... let me know if other information 
> will help.
> 
> 	Tim F
......deletion
> 
>   // set up solver
>   T = gsl_multimin_fminimizer_nmsimplex;
>   s = gsl_multimin_fminimizer_alloc(T, p);
>   ss = gsl_vector_alloc(p);
>   gsl_multimin_fminimizer_set(s, &f, &param.vector, ss);
> 
you did not initialize ss! ss should be initialized by a reasonable step
size for each of the parameters. this is needed to calculate n more
corners (n parameters) for the initial simplex. the elements of ss seem
to be all zero thus the initial size of the simplex is zero as well and
the minimizer stop right away
-- 
Dr. Ivo Alxneit
Laboratory for Solar Technology   phone: +41 56 310 4092
CH-5232 Villigen                    fax: +41 56 310 2624
Paul Scherrer Institute          http://solar.web.psi.ch
Switzerland                        gnupg key: 0x515E30C7

[-- Attachment #2: Type: APPLICATION/pgp-signature, Size: 189 bytes --]

Re: simplex minimization

[Tim F]

On Fri, 3 Jan 2003 ivo.alxneit@psi.ch wrote:

> >   // set up solver
> >   T = gsl_multimin_fminimizer_nmsimplex;
> >   s = gsl_multimin_fminimizer_alloc(T, p);
> >   ss = gsl_vector_alloc(p);
> >   gsl_multimin_fminimizer_set(s, &f, &param.vector, ss);
> > 
> you did not initialize ss! ss should be initialized by a reasonable step
> size for each of the parameters. this is needed to calculate n more
> corners (n parameters) for the initial simplex. the elements of ss seem
> to be all zero thus the initial size of the simplex is zero as well and
> the minimizer stop right away
> 

doh!  I knew it would be something dumb I was doing wrong...  Thanks -
that fixed it.

	Regards,
	Tim F