[Top][All Lists]

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

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

From: John Gehman
Subject: Re: [Help-gsl] Non-linear Multi-parameter Least-Squares Example
Date: Fri, 20 Jul 2007 13:47:00 +1000

Hi Gordan --

This will make you spew --

In using integers from 0 to 40 for your xdata, and computing data with (b/c) as part of your function, you're introducing an undefined number into the whole process. Hence the residual will be undefined, and the sum of the squares of the residuals will be undefined, and the gradients at that point, etc etc. At a minimum try float f = i +1.0 where the data is built, and in the _f and _df functions.

That's surely the biggest problem. But you also might want to be careful about mixing your float data types with the doubles assumed by gsl structures. I didn't look closely, but was left with the impression that they were mixed a bit in your code. You might be getting away with it anyway, but just a word of caution.

Also, with regard to the ultimate intent of fitting the data, I presume this was just an incremental step in your development of the fitting procedure, but insofar as this goes, once you clean up the above, you'll still get a shite fit, as the harmonic function, is not described well by the data sampling. Have a look by plotting the randomized-generated data vs the theoretical function it was based on.


John Gehman
Research Fellow
School of Chemistry
University of Melbourne

On 19/07/2007, at 6:50 PM, Gordan Bobic wrote:

I know, that's why I can't figure out what's going wrong. Attached is my
modified sinfit.c, adapted from expfit.c and merged with nlfit.c

The expfit example works, the sinfit doesn't. :-(

I'm sure I'm missing something important, but from my limited
understanding of GSL, I can't see what.



On Thu, 19 Jul 2007, John Gehman wrote:

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_...



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:

        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?


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

reply via email to

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