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

[Gordan Bobic]

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!