gsl_fit: Multidimensional 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.


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 $\chi^2$ value and the gradient of $\chi^2$ together.
double	mdmin_f (const gsl_vector parameters, void void_gsl_fit_struct)
	A routine which calculates the $\chi^2$ .
void	mdmin_df (const gsl_vector parameters, void void_gsl_fit_struct, gsl_vector *df)
	A routine which calculates the gradient of $\chi^2$ .
double	mdmin_derivatives (double p, void *void_gsl_fit_struct)
	A routine which calculates $\chi^2$ with a new parameter.
double	mdmin_chisq (gsl_fit *gft)
	A routine which calculates $\chi^2$ .
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 $\chi^2$ , with respect to a set of parameters $\mathbf{a}=[a_1,a_2,\ldots]$ . The function $\chi^2$ is given by,

$\chi^2 = \frac{1}{n-p} \sum_{i} \left(\frac{f(x_i;\mathbf{a}) - y_i}{\sigma_i}\right)^2$

where the sum is over the data points $(x_i,y_i,\sigma_i)$ , $n$ is the number of data points, and $p$ the total number of parameters. This is not exactly correct, since $p$ 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,

$\left|\frac{\partial\chi^2}{\partial a_i}\right| \leq \epsilon$

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     
00055     status = gsl_multimin_test_gradient(solvr->gradient, 1e-6);
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, $\mathbf{a}$ .
	void_gsl_fit_struct	A void pointer, expected to be a `gsl_fit`.
	f	A double evaluated to be $\chi^2$ .
	df	A gsl_vector, where the derivative of the model function with respect to each parameter: `df[j]` $=\partial \chi^2/\partial a_j$ .

Returns:: Void.

This routine calculates both the $\chi^2$ value, and the gradient of $\chi^2$ 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, $\mathbf{a}$ .
	void_gsl_fit_struct	A void pointer, expected to be a `gsl_fit`.

Returns:: $\chi^2$ as a double.

This routine calculates $\chi^2$ . This function does nothing but call the function mdmin_chisq to do the work of calculating $\chi^2$ .

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, $\mathbf{a}$ .
	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]` $=\partial \chi^2/\partial a_j$ .

Returns:: Void.

This routine calculates the gradient of $\chi^2$ 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 $\chi^2$ 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 $\partial \chi^2/\partial a_j=0$ .

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 $\chi^2$ with.
	void_gsl_fit_struct	A void pointer, expected to be a `gsl_fit`.

Returns:: $\chi^2$ as a double.

This routine calculates $\chi^2$ 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:: $\chi^2$ as a double.

This routine calculates $\chi^2$ 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 }

Multidimensional minimization

Detailed Description

Functions

Function Documentation