Multidimensional minimization

Detailed Description

Multi-dimensional minimization is a collection of iterative algorithms for finding minima of arbitrary functions. In this case, we are minimizing the Chi-squared measure of the fit for a given set of parameters. All algorithms are either some form of the conjugate gradient or steepest descent methods, and numerically calculate the gradient of with respect to the parameters to choose the direction of minimization.

The routine is implemented using the GNU Scientific library. The user can choose which of the four complex conjugate algorithms available in the GSL to use. None of the routines, however, are quite as robust as the least squares method.

No fitting errors are returned for this method.

Functions

int mdmin_main (void *void_gsl_fit_struct)
The main function for performing multidimensional minimization.
void mdmin_fdf (const gsl_vector *parameters, void *void_gsl_fit_struct, double *f, gsl_vector *df)
A routine which calculates the value and the gradient of together.
double mdmin_f (const gsl_vector *parameters, void *void_gsl_fit_struct)
A routine which calculates the .
void mdmin_df (const gsl_vector *parameters, void *void_gsl_fit_struct, gsl_vector *df)
A routine which calculates the gradient of .
double mdmin_derivatives (double p, void *void_gsl_fit_struct)
A routine which calculates with a new parameter.
double mdmin_chisq (gsl_fit *gft)
A routine which calculates .
void mdmin_status (gsl_fit *gft, int iteration, gsl_multimin_fdfminimizer *s)
A fitting status update printed to stderr.

Function Documentation

 int mdmin_main ( void * void_gsl_fit_struct )

Parameters:
 void_gsl_fit_struct A void pointer, expected to be a gsl_fit.
Returns:
GSL error code.
The method minimizes , with respect to a set of parameters . The function is given by,

where the sum is over the data points , is the number of data points, and the total number of parameters. This is not exactly correct, since should be the number of parameters not being held fixed.

The iteration can be stopped by an error, or by successful completion. The stopping criterion is the gradient falls below a tolerance level. The test can be summarized as,

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

One difficulty is the initial step size, since there is only one that must be shared among all the variables.

Definition at line 10 of file gsl_mdmin.c.

00011 {
00012   int                         count;
00013   int                         status;
00014   gsl_vector                 *guess;
00015   gsl_multimin_function_fdf   mdmin_func;
00016   gsl_multimin_fdfminimizer  *solvr;
00017   gsl_fit                    *gft;
00018
00019   gft = (gsl_fit *)void_gsl_fit_struct;
00020
00021   /* Inital guess */
00022   guess = gsl_vector_alloc(gft->para.num);
00023   for (count=0; count<gft->para.num; count++)
00024   {
00025     gsl_vector_set(guess, count, gft->para.par[count]);
00026   }
00027
00028   mdmin_func.f      = &mdmin_f;
00029   mdmin_func.df     = &mdmin_df;
00030   mdmin_func.fdf    = &mdmin_fdf;
00031   mdmin_func.n      = gft->para.num;
00032   mdmin_func.params = (void *)gft;
00033
00034   solvr = gsl_multimin_fdfminimizer_alloc(gft->meth.mdmin_type, \
00035                                           gft->para.num);
00036
00037   /*==================================================*/
00038   count = 0;
00039   gsl_multimin_fdfminimizer_set(solvr, &mdmin_func, guess, 0.01, 1e-6);
00040   mdmin_status(gft, count, solvr);
00041
00042   do
00043   {
00044     count++;
00045
00046     status = gsl_multimin_fdfminimizer_iterate(solvr);
00047     if (status != GSL_SUCCESS)
00048     {
00049       fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
00050       break;
00051     }
00052
00053     mdmin_status(gft, count, solvr);
00054
00056
00057     if (status == GSL_SUCCESS)
00058     {
00059       fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
00060     }
00061
00062   } while ((status == GSL_CONTINUE) && (count < gft->maxi));
00063
00064   mdmin_status(gft, -1, solvr);
00065
00066   gsl_multimin_fdfminimizer_free(solvr);
00067   gsl_vector_free(guess);
00068
00069   return GSL_SUCCESS;
00070 }


 void mdmin_fdf ( const gsl_vector * parameters, void * void_gsl_fit_struct, double * f, gsl_vector * df )

Parameters:
 parameters A gsl_vector containing the current state of the fit parameters, . void_gsl_fit_struct A void pointer, expected to be a gsl_fit. f A double evaluated to be . df A gsl_vector, where the derivative of the model function with respect to each parameter: df[j].
Returns:
Void.
This routine calculates both the value, and the gradient of with respect to the parameters. This function does nothing but call the functions mdmin_f and mdmin_df to do the work.

Definition at line 72 of file gsl_mdmin.c.

00075 {
00076   *f = mdmin_f(parameters, void_gsl_fit_struct);
00077   mdmin_df(parameters, void_gsl_fit_struct, df);
00078   return;
00079 }


 double mdmin_f ( const gsl_vector * parameters, void * void_gsl_fit_struct )

Parameters:
 parameters A gsl_vector containing the current state of the fit parameters, . void_gsl_fit_struct A void pointer, expected to be a gsl_fit.
Returns:
as a double.
This routine calculates . This function does nothing but call the function mdmin_chisq to do the work of calculating .

In this function, the parameters are loaded into the gsl_fit before mdmin_chisq is called. This may not be necessary. It means that the parameters in the gsl_fit are overwritten at every iteration.

Definition at line 81 of file gsl_mdmin.c.

00083 {
00084   int     i;
00085   double  tmp;
00086   gsl_fit   *gft;
00087
00088   gft = (gsl_fit *)void_gsl_fit_struct;
00089
00090   /*
00091   **  Load the params back into the tmp parameter array
00092   */
00093   for(i=0; i<gft->para.num; i++)
00094   {
00095     gft->para.tmp[i] = gsl_vector_get(parameters, i);
00096   }
00097
00098   tmp = mdmin_chisq(gft);
00099
00100   return tmp;
00101 }


 void mdmin_df ( const gsl_vector * parameters, void * void_gsl_fit_struct, gsl_vector * df )

Parameters:
 parameters A gsl_vector containing the current state of the fit parameters, . void_gsl_fit_struct A void pointer, expected to be a gsl_fit. df A gsl_vector, where the derivative of the model function with respect to each parameter: df[j].
Returns:
Void.
This routine calculates the gradient of with respect to the parameters.

It uses gsl_deriv_central as a numerical method to calculate the derivative. In order for this to work, we need a sub-fucntion that calculates with the new parameter. This is what mdmin_derivatives does. In order for mdmin_derivatives to know which parameter is currently being tested, the index gfc->para.idx is set here.

If a parameter is to be held fixed, then .

In this function, the parameters are loaded into the gsl_fit before mdmin_chisq is called. This may not be necessary. It means that the parameters in the gsl_fit are overwritten at every iteration.

Definition at line 103 of file gsl_mdmin.c.

00106 {
00107   int                i;
00108   double             res, err;
00109   gsl_function       F;
00110   gsl_fit *gft;
00111
00112   gft = (gsl_fit *)void_gsl_fit_struct;
00113
00114   /*
00115   **  Load the params back into the parameter array
00116   */
00117   for(i=0; i<gft->para.num; i++)
00118   {
00119     gft->para.par[i] = gsl_vector_get(parameters, i);
00120   }
00121
00122   /*
00123   ** Pass the gsl_fit as a void pointer to the
00124   ** function calculating the derivative.
00125   */
00126   F.params   = void_gsl_fit_struct;
00127   F.function = &mdmin_derivatives;
00128
00129   /*
00130   **  d(chisq)/d(p_i)
00131   **  => i is the parameter index
00132   ** If the parameter is to be held fixed, then
00133   ** set the derivative to zero.
00134   */
00135   for (i=0; i<gft->para.num; i++)
00136   {
00137     gft->para.idx = i;
00138     if (gft->para.fix[i]) res=0.0;
00139     else gsl_deriv_central(&F, gft->para.par[i], 1e-4, &res, &err);
00140     gsl_vector_set(df, i, res);
00141   }
00142
00143   return;
00144 }


 double mdmin_derivatives ( double p, void * void_gsl_fit_struct )

Parameters:
 p The new parameter to evaluate with. void_gsl_fit_struct A void pointer, expected to be a gsl_fit.
Returns:
as a double.
This routine calculates with a slightly changing parameter.

This routine is called by gsl_deriv_central as a numerical method to calculate the derivative. The parameter is currently being tested is found by the index _para.idx. A complete set of parameters is loded into _para:tmp, and passed to mdmin_chisq.

Definition at line 146 of file gsl_mdmin.c.

00147 {
00148   int                i;
00149   double             tmp;
00150   gsl_fit *gft;
00151
00152   gft = (gsl_fit *)void_gsl_fit_struct;
00153
00154   /*
00155   **  We need a complete set of paramters,
00156   **  with the derivative parameter 'p'
00157   **  replacing the parameter at index 'd->para.idx'.
00158   */
00159   for(i=0; i<gft->para.num; i++) gft->para.tmp[i] = gft->para.par[i];
00160   gft->para.tmp[gft->para.idx] = p;
00161
00162   tmp  = mdmin_chisq(gft);
00163
00164   return tmp;
00165 }


 double mdmin_chisq ( gsl_fit * gft )

Parameters:
 gft A gsl_fit.
Returns:
as a double.
This routine calculates based on the temporary parameters in gfc->para.tmp .

Definition at line 167 of file gsl_mdmin.c.

00168 {
00169   int    i, j, k;
00170   double tmp;
00171   double residual = 0.0;
00172   int    numdata;
00173
00174   /*
00175   ** Loop over the data sets with index k.
00176   */
00177   numdata = 0;
00178   for (k=0; k<gft->numd; k++)
00179   {
00180     /*
00181     ** Loop over the data points with index i
00182     */
00183     for (i=0; i<gft->data[k].num; i++)
00184     {
00185       /*
00186       **   Loop over the functions with index j
00187       */
00188       tmp = 0.0;
00189       for (j=0; j<gft->numf; j++)
00190       {
00191         if (gft->type == MULTIPLE)
00192         {
00193           tmp += gft->func[j].fun_m(gft->data[k].xpt[i], \
00194                                     gft->para.tmp + gft->func[j].idx, \
00195                                     k, j);
00196         }
00197         else
00198         {
00199           tmp += gft->func[j].fun_s(gft->data[k].xpt[i], \
00200                                     gft->para.tmp + gft->func[j].idx);
00201         }
00202       }
00203       tmp = (tmp - gft->data[k].ypt[i]) / gft->data[k].sig[i];
00204       residual += tmp * tmp;
00205       numdata++;
00206     }
00207   }
00208
00209   tmp =  residual / (numdata - gft->para.num);
00210   return tmp;
00211 }


 void mdmin_status ( gsl_fit * gft, int iteration, gsl_multimin_fdfminimizer * s )

Parameters:
 gft A gsl_fit. iteration The current iteration of the fitting method s A gsl_multimin_fdfminimizer containing the current state of the fit.
Returns:
Void
Prints the status of the fit to stderr

Definition at line 213 of file gsl_mdmin.c.

00216 {
00217   int      j;
00218   double   chisq;
00219
00220   if (iteration >= 0)
00221   {
00222     fprintf(stderr,"ITERATION:%3u ", iteration);
00223     fprintf(stderr,"===================\n");
00224   }
00225   else
00226   {
00227     fprintf(stderr,"COMPLETE ");
00228     fprintf(stderr,"===================\n");
00229   }
00230
00231   /*===========*/
00232   for (j=0; j<gft->para.num; j++)
00233   {
00234     gft->para.par[j] = gsl_vector_get(s->x, j);
00235     gft->para.tmp[j] = gsl_vector_get(s->x, j);
00236   }
00237   chisq = mdmin_chisq(gft);
00238   fprintf(stderr,"ChiSq: %12.6e\n", chisq);
00239   for (j=0; j<gft->para.num; j++)
00240   {
00241     fprintf(stderr,"%3d: %12.6f\n", j, gsl_vector_get(s->x, j));
00242   }
00243   fprintf(stderr,"\n");
00244
00245   /*===========*/
00246   if (iteration < 0)
00247   {
00248     for (j=0; j<gft->para.num; j++)
00249     {
00250       gft->para.par[j] = gsl_vector_get(s->x, j);
00251       gft->para.err[j] = 0.0;
00252     }
00253     gft->chis = chisq;
00254   }
00255
00256   return;
00257 }


Generated on Fri Jan 19 14:54:27 2007 for gsl_fit by  1.4.7