Multidimensional and Multi-Parameter Fortran-Based Curve Fitting Tools

The Levenberg-Marquardt algorithm has become a popular method in nonlinear curve fitting works. In this paper, following the steps of Levenberg-Marquardt algorithm, we extend the framework of the algorithm to two and three dimensional real and complex functions. This work briefly describes the mathematics behind the algorithm, and also elaborates how to implement it using FORTRAN 95 programming language. The advantage of this algorithm, when it is extended to surfaces and complex functions, is that it makes researchers to have a better trust during fitting. It also improves the generalization and predictive performance of 2D and 3D real and complex functions.

Let us construct the chi-square function: is called residue function. The goal of the least square method is to determine the parameters P of the regression function ( ) P X f , so as to minimize the squared deviations between i f and ( ) P X f i , for all data points: If we assume that all measured values of i f are normally distributed with standard deviations given by , i σ then 'statistically-the-best' match would correspond to the minimal value of 2 χ . Thus, the suitable model is essentially the one which gives the minimum value of the chi-square with respect to the parameters. That is why the method itself is called the 'least-square' technique. Of course, the error bars are determined not only by a statistical noise, but also by systematic inaccuracies, which are very difficult to estimate and are not normally distributed. However, to move on, we assume that they are (2) Where J is called Jacobian matrix of the residue ) ( c i P r which is defined in Eqn. 1. The one-half coefficient is put to simplify the formulas. To improve the fit, we can shift the parameters , k kc kc The steepest descent strategy is justified, when one is far from the minimum, but suffers from slow convergence in the plateau close to the minimum, especially in the multi- The chi function (which is quadratic) to be minimized has almost parabolic shape. The Hessian matrix, which is proportional to the curvature of 2 χ , is given by (the one-half here is also added for the sake of simplicity). The components kl α of the Hessian matrix in Eqn. (7) depends both on the first derivative, The Second derivative can be ignored when it is zero, or small enough to be negligible when compared to the term involving the first derivative. In practice, this is quite often small enough to neglect. If one looks at Eqn. (7) carefully, the second derivative is . For the successful model, this term should just be the random measurement error of each point. This error can have either sign, and should in general be uncorrelated with the model. Therefore, the second derivative terms tend to cancel out when summed over time i . Inclusion of second derivative term can in fact be destabilizing if the model fits badly or is contaminated by outlier points that are unlikely to be offset by compensating points of opposite sign. So, instead of Eqn. (7) we shall define the α-matrix simply as: After computing, numerically or analytically, the gradient and Hessian matrices for the current set of parameters, one can immediately move to the minimum by shifting the parameters , where the displacement vector k p δ is determined from the linear system derived in Eqn. (5), i.e., One of the problems associated with Newton's method (Levenberg, 1944;Kelley, 1999;Madsen, et al., 2004;and Lawson & R.J. Hanson, 1974)

The Levenberg-Marquardt Algorithm
In order for the chi-square function to converge to a minimum rapidly, one needs a large step in the direction along with the low curvature (near the minimum) and a small step in the direction with the high curvature (i.e. a steep incline). The gradient descent and Gauss-Newton iterations provide additional advantages. The LM algorithm is based on the self-adjustable balance between the two minimizing strategies: the Vanilla Gradient Descent and the Inverse Hessian methods.
Coming back to the steepest descent technique 2 χ is dimensionless but k β has the same dimension as The update rule is used as follows. If the error goes down following an update, it implies that our quadratic assumption on 2 χ is working and reduce λ (usually by a factor of 10) to reduce the influence of gradient descent. On the other hand, if the error goes up, we would like to follow the gradient more and so λ is increased by the same factor. If the initial guess is good but 2 χ does not fall down to the required minimum value, we have to change the initial value of λ slightly.

IMPLEMENTATION OF THE LM ALGORITHM
In this paper Gauss's elimination and Gauss's Jordan matrix inversion methods are used to determine the shift parameters. Among the several tests made on real and complex non linear functions, only three examples are illustrated to see how much this method is effective and faster than the other methods.

Test on real three dimensional wave function
The first test is applied to two dimensional data coordinate ( )   (yellow) after iteration From the above results (Table.1), one can easily see that the data (surface) follows the wave function having the form We have then made a fitting, using the LM approach, in order to find the values of the parameters ( ) Fig. 3 (a) and (b)).
In this case the dimension is As one can see from the above results, the LM model is highly useful when it is implemented to complicated-shaped surfaces. What is also important here is here that selecting an appropriate type of function (such as sine, power, decay, etc functions) and lambda. The shift parameters are not that much changed by normalized random errors only minimum of chi-function increases. Hence, based on the above two figures (Figs. 3 (a) and (b)), one can conclude that new equations/relations and modifications to the already existing formulas can be obtained from experimental data having disturbed/complicated surfaces.

Test made on complex two dimensional function
In ellipsometery the complex ratio

Test on complex two dimensional power function
The third test was made on complex three dimensional power functions (their derivatives are logarithmic functions). Consider the following experimental data: Table 3. Experimental data on 2D power functions.    Based on the above results we can conclude that the LM algorithm is popular method and has the following advantages

(i)
The parameters converge rapidly around the minimum in multi dimensional surfaces with complicated landscapes.
(ii) Even though the initial guess is poor, LM fits partly/most of the parameters to make fresh start.

(iii)
The convergence speed needed to reach the minimum, is not significantly influenced by the number of parameters.
(iv) The shift parameters are not that much changed by normalized random errors.
Only the minimum of the chi-function increases.
(v) Normalized random errors do not bring much change on the convergence speed, etc. Like any other non-linear optimization techniques, the LM algorithm method