Hi, I'm trying to replicate the Levenberg-Marquardt optimization on an exponential curve from here: www .csse.uwa.edu.au/programming/gsl-1.0/gsl-ref_35.html in alglib with F#.
When I run what I have, it calculates the Jacobian and then terminates saying it met the gradient condition(4).
I figure I'm either calculating the Jacobian incorrectly or I'm using some other part of the library incorrectly. I've been stuck on this for a day now, and I figure it should be something obvious about how to setup the end conditions or a fundamental misunderstanding of the optimizer.
Here's the alglib code in F# I'm using now:
Code:
let expbF (expData : float[]) (x : float[]) (fi : float[]) myobj =
for i in 0..fi.Length-1 do
fi.[i] <- ((x.[0] * System.Math.Exp(-x.[1] * float i) + x.[2]) - expData.[i]) / 0.1
let expbDf (x : float[]) (fi : double[]) (jac : double[,]) myobj =
let a = x.[0]
let lambda = x.[1]
for i in 0..fi.Length-1 do
let t = float i
let e = System.Math.Exp(-lambda * t)
jac.[i, 0] <- e / 0.1
jac.[i, 1] <- -t * a * e / 0.1
jac.[i, 2] <- 1.0 / 0.1
let testDV rawExpData =
let mutable fitExp : float[] = [|1.0; 0.0; 0.0; |]
let mutable state : alglib.minlmstate = (alglib.minlmstate())
let mutable rep : alglib.minlmreport = alglib.minlmreport()
alglib.minlmcreatevj(Array.length rawExpData, fitExp, &state)
alglib.minlmsetcond(state, 0.0, 0.0, 0.0, 0)
alglib.minlmoptimize(state, alglib.ndimensional_fvec(expbF rawExpData), expbDf, null, null )
alglib.minlmresults(state, &fitExp, &rep)
rep.terminationtype, alglib.ap.format(fitExp, 2)
I'm calling testDV with an array with the same values for the data that I'm generating in C++, but it doesn't matter because it never even calls expbF. I tried to follow the source, but I got really lost in the many different cases/gotos of minlmiteration.
And I'll in line my slightly modified GSL C++ code that I'm running to compare as well (VS2015):
Code:
namespace GSL {
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_fit.h>
}
struct data {
size_t n;
const double * y;
};
int expb_f(const GSL::gsl_vector * x, void *data, GSL::gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
const double *y = ((struct data *)data)->y;
double A = gsl_vector_get(x, 0);
double lambda = gsl_vector_get(x, 1);
double b = gsl_vector_get(x, 2);
size_t i;
for (i = 0; i < n; i++)
{
/* Model Yi = A * exp(-lambda * i) + b */
double t = i;
double Yi = A * exp(-lambda * t) + b;
gsl_vector_set(f, i, (Yi - y[i]) / 0.1);
}
return GSL::GSL_SUCCESS;
}
int expb_df(const GSL::gsl_vector * x, void *data, GSL::gsl_matrix * J)
{
size_t n = ((struct data *)data)->n;
double A = GSL::gsl_vector_get(x, 0);
double lambda = GSL::gsl_vector_get(x, 1);
size_t i;
for (i = 0; i < n; i++)
{
/* Jacobian matrix J(i,j) = dfi / dxj, */
/* where fi = (Yi - yi)/sigma[i], */
/* Yi = A * exp(-lambda * i) + b */
/* and the xj are the parameters (A,lambda,b) */
double t = i;
double s = 0.1;
double e = exp(-lambda * t);
gsl_matrix_set(J, i, 0, e / s);
gsl_matrix_set(J, i, 1, -t * A * e / s);
gsl_matrix_set(J, i, 2, 1 / s);
}
return GSL::GSL_SUCCESS;
}
int expb_fdf(const GSL::gsl_vector * x, void *data, GSL::gsl_vector * f, GSL::gsl_matrix * J)
{
expb_f(x, data, f);
expb_df(x, data, J);
return GSL::GSL_SUCCESS;
}
double* getTestData() {
const size_t n = 40;
auto y = new double[n];
GSL::gsl_rng_env_setup();
auto type = GSL::gsl_rng_default;
GSL::gsl_rng * r;
r = GSL::gsl_rng_alloc(type);
for (size_t i = 0; i < n; i++)
{
double t = i;
y[i] = 1.0 + 5 * exp(-0.1 * t) + gsl_ran_gaussian(r, 0.1);
};
return y;
}
int main(void)
{
double* const y = ExpFit::getTestData();
const GSL::gsl_multifit_fdfsolver_type * T = GSL::gsl_multifit_fdfsolver_lmsder;
GSL::gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(T, 40, 3);
struct data d = { 40, y };
GSL::gsl_multifit_function_fdf f;
f.f = &expb_f;
f.df = &expb_df;
f.fdf = &expb_fdf;
f.n = 40;
f.p = 3;
f.params = &d;
double x_init[3] = { 1.0, 0.0, 0.0 };
GSL::gsl_vector_view x = GSL::gsl_vector_view_array(x_init, 3);
gsl_multifit_fdfsolver_set(s, &f, &x.vector);
int status;
unsigned int iter = 0;
do
{
iter++;
status = gsl_multifit_fdfsolver_iterate(s);
if (status)
{
break;
}
status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
} while (status == GSL::GSL_CONTINUE && iter < 500);
gsl_multifit_fdfsolver_free(s);
delete[](y)
return 0;
}
Thoughts?
At the very least this code should be a good example of this for others when it's done. I assumed this sort of thing would be more common.
Login
Register
FAQ
Search
