Non-linear least squares fitting: Levemberg-Marquardts

Given a group of M pairs of experimental data (x₁, y₁), (x₂, y₂), ... , (x_M, y_M), having to fulfill certain function, and which depend on the variable x as well as on a group of N parameters a_j

$y (x) = y (x, a_{1}, a_{2}, . . ., a_{N})$

If the dependency of the a_j parameters is non-linear, the y(x) function is expanded by means of a Taylor series of first order around the point a_j⁰ (the initial values of the coefficients)

$y (x) = y_{0} (x) + \sum_{j = 1}^{N} [\frac{\partial y_{0} (x)}{\partial a_{j}} δ a_{j}]$

where

$δ a_{j} = a_{j} - a_{j}^{0}$

and y_j represents the value of the function for a_j.

The least squares method requires to minimize the expression:

$χ^{2} = \sum_{i = 1}^{M} [\frac{1}{σ_{i}^{2}} {\{y_{i} - y_{0} (x_{i}) - \sum_{j = 1}^{N} \frac{\partial y_{0} (x_{i})}{\partial a_{j}} δ a_{j}\}}^{2}]$

with respect to each one of the coefficients, i.e.:

$\frac{\partial χ^{2}}{\partial δ a_{k}} = 0 k = 1, . . ., N$

which leads to the following system of equations:

$\sum_{i = 1}^{M} \sum_{j = 1}^{N} \frac{1}{σ_{i}^{2}} \frac{\partial y_{o} (x_{i})}{\partial a_{j}} \frac{\partial y_{o} (x_{i})}{\partial a_{k}} δ a_{j} = \sum_{i = 1}^{M} \frac{1}{σ_{i}^{2}} [y_{i} - y_{0} (x_{i})] \frac{\partial y_{o} (x_{i})}{\partial a_{k}} k = 1, . . ., N$

which can be written in matrix notation:

$α δa = β$

where:

$α_{j k} = \sum_{i = 1}^{M} [\frac{1}{σ_{i}^{2}} \frac{\partial y_{0} (x_{i})}{\partial a_{j}} \frac{\partial y_{0} (x_{i})}{\partial a_{k}}]$

$β_{k} = \sum_{i = 1}^{M} [\frac{1}{σ_{i}^{2}} [y_{i} - y_{0} (x_{i})] \frac{\partial y_{0} (x_{i})}{\partial a_{k}}]$

The method of Marquardt exchanges matrix α in the previous equation by matrix α':

$α' δa = β$

$α'_{j j} = α_{j j} (1 + λ)$

$α'_{i j} = α_{i j} (i \neq j)$

Library functions (Numerical Recipes):

mrqmin():Go to function call

it takes twelve arguments; vectors x[], y[] (which contain the experimental points), sig[] (which contains the experimental points' uncertainties), the number of data points ndata, vectors a[] (which contains the coefficients of the function) and ia[] (which contains the values "true" or "false" to fit the corresponding parameter or keep its value unchanged), the variable covar[][] (which is the number of arguments), the matrices covar[][] (used in the calculation process, where the covariance matrix is returned) and alpha[][] (which contains the curvature matrix α'), the variable chisq[] (where the value of χ² is returned), the variable funcs() (which contains the memory address of the user function), and variable alamda (which contains the λ parameter of the Marquardt function).

mrqcof():Go to function call

It's a function used by mrqmin(): and takes as arguments the vectors x[], y[], sig[], a[] and ia[], the alpha[] matrix, the variable func and the vector beta[] (which contains the values of β_k).

func():Go to function call

The user must implement the fitting function with the following arguments: the x variable (which contains the value of the abcissa x_i where the function will be evaluated), the vector a[], the variable y (where the value of the function for x_i will be returned) and the vector dyda[] (which will return the values of the derivatives with respect to each of the evaluated parameters):

$\frac{\partial σ}{\partial σ_{0}} = \frac{1}{(E - E_{r})^{2} + \frac{γ^{2}}{4}}$

$\frac{\partial σ}{\partial E_{r}} = \frac{2 σ_{0} (E - E_{r})}{{[(E - E_{r})^{2} + \frac{γ^{2}}{4}]}^{2}}$

$\frac{\partial σ}{\partial γ} = - \frac{σ_{0} γ}{2 {[(E - E_{r})^{2} + \frac{γ^{2}}{4}]}^{2}}$

The files for this calculation are in folder breitwigner/3_Marquardt. It uses the Numerical Recipes libraries mentioned above, plus covsrt.c and gausj.c.

The program marquardt.c reads the experimental data from data.txt and writes the results in an output.txt file. The format of this file is given in three columns, where the first represents the energy, the second, the cross-section measured at that energy, and the third, the uncertainties. The values are taken from the experimental data of neutron scattering cross-sections provided in the main page.

In this case we have 3 parameters, σ, E_r and γ. We have to provide initial values for them to start the iteration, and they will be adjusted by the mrqmin() function. The program fits all the data points by default, with ia[i] = 1.

	/***************************************************
	 * Initial values for the parameter vectors
	 * For example:
	 *    a[1] = 70000.0;
	 *    a[2] = 75.0;
	 *    a[3] = 60.0;
	 ***************************************************/
	for (i=1; i<=N; i++) {
		printf("Enter initial value of parameter %d: ", i);
		scanf("%lf", &a[i]);
	}

	// Fit all
	for (i=1; i<=N; i++)
		ia[i] = 1;

Once we have all the inputs, we can call the mrqmin() function (you may want to check the definition again).

	/***************************************************
	 * ITERATION LOOP
	 ***************************************************/
	char answer;
	do {
		printf("\nHow many iterations?: ");
		scanf("%d", &k);

		for(i=1; i<=k; i++) {
			iter++;
			ochisq = chisq;

			mrqmin(x, y, sig, M, a, ia, N, covar, alpha, &chisq, funcs, &alamda);

			if (chisq > ochisq)
				n=0;
			else if (fabs(ochisq - chisq) < 0.1)
				n++;

			if(n >= 4)
				break;
		}

		if(n < 4) {
			printf("\nIterate again? (Y/N)?: ");
			scanf("%c", &answer);
		}
	} while (answer == 's' || answer == 'S');

Finally, we can use the user defined function func() to recalculate the experimental points (see func() definition ).

	printf(     "\nData points recalculated from the fit:\n");
	fprintf(fp, "\nData points recalculated from the fit:\n");
	for (i=1; i<=M; i++) {

		(*funcs)(x[i], a, &ymod[i], dyda, N);
	
		sum = 0.0;
		for (j=1; j<=N; j++)
			sum += pow(dyda[j]*sqrt(covar[j][j]), 2);

		printf(     "\nSigma[%d] = %3.3lf\t\tdSigma[%d] = %3.3lf", i, ymod[i], i, sqrt(sum));
		fprintf(fp, "\nSigma[%d] = %3.3lf\t\tdSigma[%d] = %3.3lf", i, ymod[i], i, sqrt(sum));
	}

Results

The following values with their errors are obtained for the parameters of interest:

σ₀ = 66900.3
E_r = 77.5 MeV
γ = 56.2 MeV

Breit-Wigner curve obtained with the method of least-squares