help-gsl
[Top][All Lists]

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

 From: Alex Brussee Subject: [Help-gsl] BFGS optmizer, how to use ? Date: Sat, 22 May 2004 12:29:04 +0000

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

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

PrintErf(s->x, Input);

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

if (status)
break;

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

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

}

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

```