Non-linear least squares

Detailed Description

Non-linear least squares fitting is a mathematical optimization technique that attempts to find a "best fit" to a set of data by attempting to minimize the sum of the squares of the differences (called residuals) between the fitted function and the data. The Levenberg-Marquardt algorithm is an iterative procedure that provides a numerical solution to non-linear least squares fitting. Please refer to the Wikipedia article on the Levenberg-Marquardt algorithm for more information. (http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm)

These robust routines are implemented using the GNU Scientific Library. By using the fit_method struct, the user can choose which of the two Levenberg-Marquardt algorithms available in the GSL to use.

Functions

int lstsq_main (Parratt *d)
The main function for Non-linear least squares.
int lstsq_fdf (const gsl_vector *params, void *data, gsl_vector *f, gsl_matrix *J)
Required by the GSL library.
int lstsq_f (const gsl_vector *params, void *data, gsl_vector *f)
Calculates the residuals.
int lstsq_df (const gsl_vector *params, void *data, gsl_matrix *J)
Calculates the Jacobian matrix.
double lstsq_derivatives1 (double p, void *data)
Numerically calculate the residual gradients.
double lstsq_derivatives2 (double p, void *data)
Numerically calculate the residual gradients.
void lstsq_status (int i, gsl_multifit_fdfsolver *s, gsl_matrix *c, Parratt *d)
Updates the fit to stdout.

Function Documentation

 int lstsq_main ( Parratt * d )

The method minimizes the sum of residuals , with respect to a set of parameters , which define the layered structure. The function to be minimized is given by,

where the sum is over the data points .

The iteration can be stopped by an error, or by successful completion. There are two stopping criterion, and the iteration stops on the success of either one. The first is when the gradient falls below a tolerance level. The test can be summarized as,

The second stop test compares the residuals, , of the last step with the residuals of the current step according to,

Perhaps a later version will allow user configuration of the tolerances.

As in the other fitting methods, the GSL library takes a (void *) struct as the fitting parameters, so this is set to be the Parratt struct. In this way the data is totally accessible. The parameters are overwritten with this technique, however, since the Parratt struct is used as workspace.

The code requires one external function, calcReflectivity(), which takes the form

double calcReflectivity(double q, double *a, Parratt *p)

Parameters:
 q The value to evaluate . a The set of all adjustable parameters in one array. p The Parratt struct containing all the other information necessary for calculation of .
Returns:
The reflectivity .
00023 {
00024     unsigned int  iter;
00025     int           status, status2, ndat;
00026     gsl_vector                         *guess;
00027     gsl_vector                         *gradt;
00028     gsl_matrix                         *covar;
00029     gsl_multifit_function_fdf           my_func;
00030     gsl_multifit_fdfsolver             *solv;
00031
00032     /* Initial guess */
00033     guess = gsl_vector_alloc(d->m.num);
00034     for (iter=0; iter<d->m.num; iter++)
00035     {
00036         gsl_vector_set(guess, iter, d->m.par[iter]);
00037     }
00038
00039     /* Need for stop test #2 */
00040     gradt = gsl_vector_alloc(d->m.num);
00041
00042     /* Co-variant matrix      */
00043     covar = gsl_matrix_alloc(d->m.num,d->m.num);
00044
00045     my_func.f      = &lstsq_f;
00046     my_func.df     = &lstsq_df;
00047     my_func.fdf    = &lstsq_fdf;
00048     /* Number of data points */
00049     if (d->s.pntp == POLARIZED) ndat = d->d1.num + d->d2.num;
00050     else                        ndat = d->d1.num;
00051     my_func.n      = ndat;
00052     /* Number of parameters */
00053     my_func.p      = d->m.num;
00054     my_func.params = (void *)d;
00055
00056     solv = gsl_multifit_fdfsolver_alloc(fit_method.lstsq_type, \
00057                                         ndat, d->m.num);
00058
00059     /*==================================================*/
00060     iter = 0;
00061     gsl_multifit_fdfsolver_set(solv, &my_func, guess);
00062     gsl_multifit_covar(solv->J, 0.0, covar);
00063     lstsq_status(iter, solv, covar, d);
00064
00065     do
00066     {
00067         iter++;
00068         status = gsl_multifit_fdfsolver_iterate(solv);
00069         gsl_multifit_covar(solv->J, 0.0, covar);
00070         lstsq_status(iter, solv, covar, d);
00071
00072         if (status)
00073         {
00074             fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
00075             break;
00076         }
00077
00078         /* Stop test #1 */
00079         status = gsl_multifit_test_delta(solv->dx, solv->x, 1e-7, 1e-7);
00080
00081         /* Stop test #2 */
00082         gsl_multifit_gradient(solv->J, solv->f, gradt);
00083         status2 = gsl_multifit_test_gradient(gradt, 1e-7);
00084
00085         if ((status==GSL_SUCCESS) || (status2==GSL_SUCCESS))
00086         {
00087             fprintf(stderr,"STATUS1: %s\n", gsl_strerror(status));
00088             fprintf(stderr,"STATUS2: %s\n", gsl_strerror(status2));
00089         }
00090
00091     } while ((status == GSL_CONTINUE)  && \
00092              (status2 == GSL_CONTINUE) && \
00093              (iter < MAX_ITER));
00094
00095     /*==================================================*/
00096     gsl_multifit_covar(solv->J, 0.0, covar);
00097     lstsq_status(-1, solv, covar, d);
00098
00099     /*==================================================*/
00100
00101     gsl_multifit_fdfsolver_free(solv);
00102     gsl_matrix_free(covar);
00103     gsl_vector_free(guess);
00104     gsl_vector_free(gradt);
00105
00106     return GSL_SUCCESS;
00107 }


 int lstsq_fdf ( const gsl_vector * params, void * data, gsl_vector * f, gsl_matrix * J )

This function is required by the GSL, and merely calls the functions that fills the vector f with the residuals and the matrix J with the derivatives.

00110 {
00111     lstsq_f(params, data, f);
00112     lstsq_df(params, data, J);
00113     return GSL_SUCCESS;
00114 }


 int lstsq_f ( const gsl_vector * params, void * data, gsl_vector * f )

The array f is the residuals .

00117 {
00118     int     i, j;
00119     double  tmp;
00120     Parratt    *d;
00121
00122     d = (Parratt *)data;
00123
00124     for (i=0; i<d->m.num; i++)
00125     {
00126         d->m.par[i] = gsl_vector_get(params,i);
00127     }
00128
00129     d->g.plz = POLARIZED_UP;
00130     for(i=0; i<d->d1.num; i++)
00131     {
00132         tmp = (calcReflectivity(d->d1.xpt[i], d->m.par, (Parratt *)data) \
00133                               - d->d1.ypt[i]) \
00134                               / d->d1.sig[i];
00135         gsl_vector_set(f, i, tmp);
00136     }
00137
00138     if (d->s.pntp == POLARIZED)
00139     {
00140         d->g.plz = POLARIZED_DOWN;
00141         for(i=0, j=d->d1.num; i<d->d2.num; i++, j++)
00142         {
00143             tmp = (calcReflectivity(d->d2.xpt[i], d->m.par, (Parratt *)data) \
00144                                   - d->d2.ypt[i]) \
00145                                   / d->d2.sig[i];
00146             gsl_vector_set(f, j, tmp);
00147         }
00148     }
00149
00150     return GSL_SUCCESS;
00151 }


 int lstsq_df ( const gsl_vector * params, void * data, gsl_matrix * J )

The gsl_matrix J is the Jacobian matrix, whose elements are

This is calculated using numerical derivatives. If the parameter is to be held fixed, then the element is set to zero.

00154 {
00155     int i, j, k;
00156     gsl_function F;
00157     double res, err;
00158     Parratt   *d;
00159
00160     d = (Parratt *)data;
00161     F.params   = data;
00162
00163     for (i=0; i<d->m.num; i++)
00164     {
00165         d->m.par[i] = gsl_vector_get(params,i);
00166     }
00167
00168
00169     F.function = &lstsq_derivatives1;
00170     d->g.plz = POLARIZED_UP;
00171     for(i=0; i<d->d1.num; i++)
00172     {
00173         d->d1.idx = i;
00174
00175         for (j=0; j<d->m.num; j++)
00176         {
00177             d->m.idx = j;
00178             if (d->m.fix[j]) res=0.0;
00179             else gsl_deriv_central(&F, d->m.par[j], 1e-4, &res, &err);
00180             gsl_matrix_set(J, i, j, res);
00181         }
00182     }
00183
00184     if (d->s.pntp == POLARIZED)
00185     {
00186         F.function = &lstsq_derivatives2;
00187         d->g.plz = POLARIZED_DOWN;
00188         for(i=0, k=d->d1.num; i<d->d2.num; i++, k++)
00189         {
00190             d->d2.idx = i;
00191
00192             for (j=0; j<d->m.num; j++)
00193             {
00194                 d->m.idx = j;
00195                 if (d->m.fix[j]) res=0.0;
00196                 else gsl_deriv_central(&F, d->m.par[j], 1e-4, &res, &err);
00197                 gsl_matrix_set(J, k, j, res);
00198             }
00199         }
00200     }
00201
00202     return GSL_SUCCESS;
00203 }


 double lstsq_derivatives1 ( double p, void * data )

This function is needed by gsl_deriv_central to return the value the residual at slightly changing values of the parameter .

In order for this trick to work, the parameter that the derivative is with respect to is flagged by d->m.idx and the data point the derivative is taken at is d->d.idx . The value of the other parameters are unchanged in this derivative, since they were put into d->m.par , back in the function mdmin_df() , they are available to fill the vector a.

00206 {
00207     int     i;
00208     double  tmp, a[MAX_PARAM];
00209     Parratt *d;
00210
00211     d = (Parratt *)data;
00212
00213     for (i=0; i<d->m.num; i++)
00214     {
00215         if (d->m.idx==i) a[i] = p;
00216         else a[i] = d->m.par[i];
00217     }
00218
00219     tmp = (calcReflectivity(d->d1.xpt[d->d1.idx], a, (Parratt *)data) \
00220                           - d->d1.ypt[d->d1.idx]) \
00221                           / d->d1.sig[d->d1.idx];
00222     return tmp;
00223 }


 void lstsq_status ( int i, gsl_multifit_fdfsolver * s, gsl_matrix * c, Parratt * d )

Parameters:
 i Current iteration s Current state of the solution
This is a report of the current state of the solution.
00246 {
00247     int      j;
00248     double   chisq;
00249
00250     if (i >= 0)
00251     {
00252         fprintf(stderr,"ITER:%3u ", i);
00253         fprintf(stderr,"===================\n");
00254     }
00255     else
00256     {
00257         fprintf(stderr,"RESULT ");
00258         fprintf(stderr,"===================\n");
00259     }
00260
00261     /*
00262      * Chisq
00263      */
00264     if (d->s.pntp == POLARIZED)
00265     {
00266         chisq = gsl_blas_dnrm2(s->f) / \
00267                (double)((d->d1.num)+(d->d2.num)-(d->m.num));
00268     }
00269     else
00270     {
00271         chisq = gsl_blas_dnrm2(s->f) / \
00272                (double)((d->d1.num)-(d->m.num));
00273     }
00274
00275     /*
00276      *  Report results
00277      */
00278     fprintf(stderr,"ChiSq: %12.6e\n", chisq);
00279     for (j=0; j<d->m.num; j++)
00280     {
00281         if (d->m.fix[j] == 0)
00282         {
00283             fprintf(stderr,"%3d: %12.6f +/- %8.6f\n",
00284                     j,
00285                     gsl_vector_get(s->x, j),
00286                     sqrt(gsl_matrix_get(c, j, j)));
00287             d->m.par[j] = gsl_vector_get(s->x, j);
00288         }
00289     }
00290     fprintf(stderr,"\n");
00291
00292
00293     return;
00294 }


Generated on Wed Mar 14 13:24:56 2007 for Parratt by  1.4.7