help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Fwd: [Help-gsl] Non-linear Multi-parameter Least-Squares Example


From: John Gehman
Subject: Fwd: [Help-gsl] Non-linear Multi-parameter Least-Squares Example
Date: Thu, 19 Jul 2007 15:07:02 +1000


Hello Gordan,

That block of code is verbatim from the example, and it has worked for me, so I reckon there's something wrong with either the component functions/variables of your gsl_multifit_function_fdf (i.e. f, df, fdf, n, p, params), or your gsl_multifit_fdfsolver_...

Cheers,
john

---------------------------------------------------------

Dr John Gehman (address@hidden)
Research Fellow; Separovic Lab
School of Chemistry
University of Melbourne (Australia)


On 18/07/2007, at 8:11 PM, Gordan Bobic wrote:

Thanks for that, really appreciate your help.

The problem I seem to have now is that the main fitting iteration loop
seems to just bail out on me immediately:

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

        printf ("status = %s\n", gsl_strerror (status));

        print_state (iter, s);

        if (status) break;

        status = gsl_multifit_test_delta (s->dx, s->x,
                                        1e-4, 1e-4);
}
while (status == GSL_CONTINUE && iter < 500);

"if (status)" always evaluates to true (status == -2, I checked) but the
error string is actually "the iteration has not converged yet". If I
change this to:

if (status != -2), then it just keeps going until the maximum number of iterations are up. So, status == GSL_CONTINUE, but the worrying thing is that at each pass the numbers output by print_state() are the same. It is
as if there is no gradient descent happening.

Am I missing something obvious? It doesn't seem to make any difference
what I put in x_init[], be it the output of the guesstimator, 1s 0s or
some other random number. It never seems to descend.

I removed all the gsl_rng* stuff because I have a fixed data block I use
for testing so I can make meaningful comparisons between different
algorithms and libraries (i.e. content of y[] is pre-generated). This
couldn't be the cause of the problem, could it?

My Jacobian setting is:

for (i = 0; i < n; i++)
{
        // Yi = a * sin(X/b + c) + d

        gsl_matrix_set (J, i, 0, sinf(b / i + c) / sigma[i]);
        gsl_matrix_set (J, i, 1, a * cosf(b / i + c) / sigma[i] / i);
        gsl_matrix_set (J, i, 2, a * cosf(b / i + c) / sigma[i]);
        gsl_matrix_set (J, i, 3, 1 / sigma[i]);
}

But even if this is wrong, surely it should still descend _somewhere_,
given random starting points, should it not?

Gordan


On Wed, 18 Jul 2007, John Gehman wrote:

G'day Gordan,

The Jacobian matrix J is the differential of your objective equation
with respect to each parameter (corresponding  to the second index)
at each data point (corresponding to the first index).

The relevant equation is  [a * sin(b / t + c) + d - yi ] / sigma[i],
i.e. you're trying to minimize the difference between the model and
the data (yi), where some points may be more reliable than others
(given sigma[i]),  by optimizing the amplitude, period, phase, and
offset of your sine wave. The Jacobian provides the analytic
equations to inform the system of the sensitivity of the residuals in
the evolving fit to changes in each of the parameters.

Taking the partial derivative of your equation with respect to each
of the floating parameters, and setting them appropriately into the
matrix J, assuming your vector of floating parameters elsewhere is
ordered (a,b,c,d), and retaining your s = sigma[i], I get:


gsl_matrix_set (J, i, 0, sin(b / t + c) / s);           # derivative
with respect to a
gsl_matrix_set (J, i, 1, cos(b / t + c) * a/(t*s) ); # derivative
with respect to b
gsl_matrix_set (J, i, 2, cos(b / t + c) * a/s);       # derivative
with respect to c
gsl_matrix_set (J, i, 3, 1/s);                                #
derivative with respect to d

The analytic derivatives are usually my problem, however, so please
confirm them!




_______________________________________________
Help-gsl mailing list
address@hidden
http://lists.gnu.org/mailman/listinfo/help-gsl




reply via email to

[Prev in Thread] Current Thread [Next in Thread]