[Help-gsl] BFGS optmizer, how to use ?

[Alex Brussee]

Hi All,

I am trying to use the bfgs minimizer to optimize A3, A2, A1, A0 in the next 
formula:
f(x) = A3*x*x*x + A2*x*x + A1*x + A0. In my program, the objective is not to 
optimize the function directly, but by comparing results from the current 
vector of function parameters (A3, A2, A1, A0) for a number of X-values to a 
set of precalculated values made with a target vector of function parameters 
(10,20,15,25 in my case). The BFGS alg. seems to be converging, but after a 
number steps it gives NaN's. I do not know what i am doing wrong, does 
anybody see what is missing from my code ?

N.B. I used http://sources.redhat.com/ml/gsl-discuss/2003-q2/msg00214.html 
as example for my code. Solving the formula directly did work.

My code:
#include <gsl/gsl_multimin.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>

#define NTARGETS 10 // there are six targets fo be fitted
#define NCOEFF 4 // there are four (a3, a2, a1, a0) to be fitted
#define STARTGRAD 0.0001 // initial search gradient

using namespace std;

typedef struct
{
        double XTargets[NTARGETS];
        double YTargets[NTARGETS];
} Data;

void PrintErf(gsl_vector *, Data *);
void PrintState(int, gsl_vector *, gsl_vector *, gsl_vector *);

double CubicFn(const gsl_vector *Params, double x)
{
        double a3 = gsl_vector_get(Params, 0); // first parameter = A3 !!!
        double a2 = gsl_vector_get(Params, 1);
        double a1 = gsl_vector_get(Params, 2);
        double a0 = gsl_vector_get(Params, 3);

        return(a3*x*x*x + a2*x*x + a1*x + a0);
}

double Erf(const gsl_vector *Params, Data *Input)
{
        double fx = 0;
        double result = 0;

        for(int i = 0;i < NTARGETS;i++)
        {
                fx = CubicFn(Params, Input->XTargets[i]) - Input->YTargets[i];
                fx = fx * fx;
                result += fx;
        }

        return result;
}

double Sum(const gsl_vector *Params)
{
        double result = 0;

        for(int i = 0;i < NCOEFF;i++)
        {
                result += fabs(gsl_vector_get(Params, i));
        }

        return result;
}

double func(const gsl_vector *Params, void *Targets)
{
   Data *Input = (Data *)Targets;

        return Erf(Params, Input);
}

void func_df (const gsl_vector *Params, void *Targets, gsl_vector 
*ParamGradients)
{
   Data *Input = (Data *)Targets;
        double ParamSum = Sum(ParamGradients);
        double Fx = 0;

        for(int i = 0;i < NCOEFF;i++)
        {
		Fx = (CubicFn(Params, Input->XTargets[i]) - Input->YTargets[i]) * 
(CubicFn(Params, Input->XTargets[i]) - Input->YTargets[i]);
                gsl_vector_set(ParamGradients, i, Fx / Erf(Params, Input));
        }
}

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

void InitData(Data *Input, const double *XVals, const gsl_vector *Params)
{
        for(int i = 0;i < NTARGETS;i++)
        {
                Input->XTargets[i] = XVals[i];
                Input->YTargets[i] = CubicFn(Params, Input->XTargets[i]);
        }
}



int main(int argv, char **argc)
{
   size_t iter = 0;
   int status = 0;

   const gsl_multimin_fdfminimizer_type *T;
   gsl_multimin_fdfminimizer *s;
   gsl_multimin_function_fdf mf;
   gsl_vector *StartParams;
   gsl_vector *Solution;

        double TargetX[NTARGETS];

        for(int i = 0;i < NTARGETS;i++) TargetX[i] = i;

        Data *Input = new Data;

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

   StartParams = gsl_vector_alloc(NCOEFF);
   Solution = gsl_vector_alloc(NCOEFF);

        gsl_vector_set(Solution, 0, 10);
        gsl_vector_set(Solution, 1, 20);
        gsl_vector_set(Solution, 2, 15);
        gsl_vector_set(Solution, 3, 25);

        InitData(Input, TargetX, Solution);

        cout << "#usage: bfgs A3 A2 A1 A0 (doubles)\n";
        cout << "#searching for solution (A3, A2, A1, A0):\n";
        cout << "#minizing f(x) = " << gsl_vector_get(Solution, 0) << "*x^3 +"
                                                                << gsl_vector_get(Solution, 1) 
<< "*x^2 +"
                                                                << gsl_vector_get(Solution, 2) 
<< "*x +"
                                                                << gsl_vector_get(Solution, 3) 
<< "\n";
        if(argv > 4)
        {
		for(int i = 0;i < NCOEFF;i++) gsl_vector_set(StartParams, i, atof(argc[i + 
1]));
        }       else for(int i = 0;i < NCOEFF;i++) gsl_vector_set(StartParams, 
i, 1);

        cout << "#Starting search using next formula:\n";
	cout << "#starting with: f(x) = " << gsl_vector_get(StartParams, 0) << 
"*x^3 +"
                                                                << gsl_vector_get(StartParams, 
1) << "*x^2 +"
                                                                << gsl_vector_get(StartParams, 
2) << "*x +"
                                                                << gsl_vector_get(StartParams, 
3) << "\n\n";

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

   gsl_multimin_fdfminimizer_set(s, &mf, StartParams, STARTGRAD, 1e-4);

        PrintErf(s->x, Input);

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

                PrintState(iter, s->x, s->dx, s->gradient);

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

   gsl_multimin_fdfminimizer_free(s);
   gsl_vector_free(StartParams);
   gsl_vector_free(Solution);

   return 0;
}


void PrintState(int iter, gsl_vector *Params, gsl_vector *dx, gsl_vector 
*Gradient)
{
//      cout    << "Iter: " << iter <<  endl;
        cout    << gsl_vector_get(Params, 0) << "\t"
                        << gsl_vector_get(Params, 1) << "\t"
                        << gsl_vector_get(Params, 2) << "\t"
                        << gsl_vector_get(Params, 3) << endl;
/*      cout    << gsl_vector_get(dx, 0) << " "
                        << gsl_vector_get(dx, 1) << " "
                        << gsl_vector_get(dx, 2) << " "
                        << gsl_vector_get(dx, 3) << endl;
        cout    << gsl_vector_get(Gradient, 0) << "\t"
                        << gsl_vector_get(Gradient, 1) << "\t"
                        << gsl_vector_get(Gradient, 2) << "\t"
                        << gsl_vector_get(Gradient, 3) << endl;*/

}

void PrintErf(gsl_vector *Params, Data *Input)
{
        double fx = 0;
        double result = 0;

        for(int i = 0;i < NTARGETS;i++)
        {
                fx = CubicFn(Params, Input->XTargets[i]) - Input->YTargets[i];
                fx = fx * fx;
		cout << Input->XTargets[i] << ": " << CubicFn(Params, Input->XTargets[i]) 
<< "\t - " << Input->YTargets[i] << "\t => " << fx << endl;
                result += fx;
        }
        cout << " total: " << result << endl << endl;
}

_________________________________________________________________
MSN Search, for accurate results! http://search.msn.nl