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

[Gordan Bobic]

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.

Thanks.

Gordan

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

sinfit.c
Description: Text document