Re: gsl_multimin question

[Rajarshi Guha <rajarshi at presidency dot com>]

On Tue, 2003-06-03 at 13:13, Fabrice Rossi wrote:
> Hi.
> 
> Could you show the code? Are you sure the derivatives are correct? 

I dont have the actual code that was giving the problem, but I modified
the example program given in the manual to use a function similar to the
one in my original code. The same error occurs.

I've attached the C program - essentially I'm trying to minimize a sum
of squares error function. In the attached code the form is

E = sqrt( \sum x_{i} - X_{i} )

where X_{i} is my target value (observed value).

I've selected the initial argument vector randomly, but the values are
close to the target vector. I expected that the values printed in each
iteration should converge to my target vector - but they dont and after
a few iterations all I see is NaN's.

As far as I can see the derivatives are properly defined. Since you
mentioned thgat your collegue gets it to work I must be doing something
wrong somewhere but I cant see it :(

Any pointers would be appreciated.
-------------------------------------------------------------------
Rajarshi Guha <rajarshi@presidency.com> <http://jijo.cjb.net>
GPG Fingerprint: 0CCA 8EE2 2EEB 25E2 AB04 06F7 1BB9 E634 9B87 56EE
-------------------------------------------------------------------
A meeting is an event at which the minutes are kept and the hours are
lost.

#include <gsl/gsl_multimin.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double func(const gsl_vector *v, void *param)
{
    int i;
    double sum = 0;
    double x,y;
    double *dp = (double *)param;

    for (i = 0; i < 4;i++)
        sum += (gsl_vector_get(v,i) - dp[i]) * (gsl_vector_get(v,i) - dp[i]);
    sum = sqrt(sum);
    return (sum);
}

void 
func_df (const gsl_vector *v, void *params, 
        gsl_vector *df)
{
    int i;
    double sum = 0;
    double x, y;
    double *dp = (double *)params;
    
    for (i = 0; i < 4;i++)
        sum += (gsl_vector_get(v,i) - dp[i]) * (gsl_vector_get(v,i) - dp[i]);
    sum = sqrt(sum);

    for (i = 0; i < 4; i++)
        gsl_vector_set(df,i, (gsl_vector_get(v,i) - dp[i]) /  sum); 
}

void 
func_fdf (const gsl_vector *x, void *params, 
        double *f, gsl_vector *df) 
{
  *f =func(x, params); 
  func_df(x, params, df);
}

int main(void)
{
    size_t iter = 0;
    int status = 0;
    
    const gsl_multimin_fdfminimizer_type *T;
    gsl_multimin_fdfminimizer *s;
    gsl_multimin_function_fdf mf;
    gsl_vector *x;

    double p[4] = {1.0,2.0,3.0,4.0};

    mf.f = &func;
    mf.df = &func_df;
    mf.fdf = &func_fdf;
    mf.n = 4;
    mf.params = (void*)p;


    x = gsl_vector_alloc (4);
    gsl_vector_set (x, 0, 5.4);
    gsl_vector_set (x, 1, 4.3);
    gsl_vector_set (x, 2, 2.5);
    gsl_vector_set (x, 3, -3.0);

    T = gsl_multimin_fdfminimizer_vector_bfgs;
    s = gsl_multimin_fdfminimizer_alloc (T, 4);

    gsl_multimin_fdfminimizer_set (s, &mf, x, 0.01, 1e-4);

    do
      {
        iter++;
        status = gsl_multimin_fdfminimizer_iterate (s);

        if (status)
            break;

        status = gsl_multimin_test_gradient (s->gradient, 1e-3);

        if (status == GSL_SUCCESS)
            printf ("Minimum found at:\n");

        printf ("%5d %.5f %.5f  %.5f  %.5f %10.5f\n", iter,
                gsl_vector_get (s->x, 0), 
                gsl_vector_get (s->x, 1), 
                gsl_vector_get (s->x, 2), 
                gsl_vector_get (s->x, 3), 
                s->f);

      }
    while (status == GSL_CONTINUE && iter < 1000);

    gsl_multimin_fdfminimizer_free (s);
    gsl_vector_free (x);

    return 0;
}