LMFnlsq - Solution of nonlinear least squares

Efficient and stable Levenberg-Marquard-Fletcher method for solving of nonlinear equations

The function The LMFnlsq.m serves for finding optimal solution of an overdetermined system of nonlinear equations in the least-squares sense. The standard Levenberg- Marquardt algorithm was modified by Fletcher and coded in FORTRAN many years ago (see the Reference). This version of LMFnlsq is its complete MATLAB implementation complemented by setting parameters of iterations as options. This part of the code has been strongly influenced by Duane Hanselman's function mmfsolve.m.

Calling of the function is rather simple and is one of the following:
LMFnlsq % for help output
Options = LMFnlsq('default');
Options = LMFnlsq(Name1,Value1,Name2,Value2,...);
x = LMFnlsq(Eqns,X0);
x = LMFnlsq(Eqns,X0,'Name',Value,...);
x = LMFnlsq(Eqns,X0,Options);
[x,ssq] = LMFnlsq(Eqns,...);
[x,ssq,cnt] = LMFnlsq(Eqns,...);
[x,ssq,cnt,nfJ] = LMFnlsq(Eqns,...);
[x,ssq,cnt,nfJ,XY] = LMFnlsq(Eqns,...);

In all cases, the applied variables have the following meaning:
% Eqns is a function name or a handle defining a set of equations,
% X0 is a vector of initial estimates of solutions,
% x is the least-squares solution,
% ssq is sum of squares of equation residuals,
% cnt is a number of iterations,
% nfJ is a sum of calls of Eqns and function for Jacobian matrix,
% xy is a matrix of iteration results for 2D problem [x(1), x(2)].
% Options is a list of Name-Value pairs, which may be set by the calls
Options = LMFnlsq; % for default values,
Options = LMFnlsq('Name',Value,...); % for users' chosen parameters,
Options = LMFnlsq(Options,'Name',Value,...); % for updating Options.
If no Options is defined, default values of options are used.

Field names 'Name' of the structure Options are:
% 'Display' for control of iteration results,
% 'MaxIter' for setting maximum number of iterations,
% 'ScaleD' for defining diagonal matrix of scales,
% 'FunTol' for tolerance of final function values,
% 'XTol' scalar or vector for tolerance of final solution increments,
% 'Trace' for control of iteration saving,
% 'Lambda' for setting of initial value of the parameter lambda.
% 'Jacobian' for a handle of function, which evaluates Jacobian matrix.
If no handle is declared, internal function for finite difference approximation of the matrix is used.

Example 1:
The general Rosenbrock's function has the form
f(x) = 100(x(2)-x(1)^2)^2 + (1-x(1))^2
Optimum solution gives f(x)=0 for x(1)=x(2)=1. Function f(x) can be expressed in the form
f(x) = f1(x)^2 + f2(x)^2, where f1(x) = 10(x(2)-x(1)^2), and f2(x) = 1-x(1).
Values of the functions f1(x) and f2(x) can be used as residuals. The parameter Eqns has a form of either named function:
% function r = rosen(x)
%% ROSEN Residuals of the Rosenbrock valey:
% r = [ 10*(x(2)-x(1)^2) % first part, f1(x)
% 1-x(1) % second part, f2(x)
% ];
or a handle of the anonymous function:
% rosen = @(x) [10*(x(2)-x(1)^2); 1-x(1)];
The calls are different:
[x,...] = LMFnlsq('rosen',...); % in case of named function, or
[x,...] = LMFnlsq(rosen,...); % in case of the function handle.
LMFnlsq finds the exact solution of this problem in 17 iterations.

Example 2:
Regression of experimental data.
Let us have experimental data simulated by the decaying function
y = c(1) + c(2)*exp(c(3)*x) + w*randn(size(x));
for column vector x. An initial guess of unknown parameters is obtained from approximation y(x) for x=0 and x->Inf as
c1 = y(end);
c2 = y(1)-c(1); and
c3=x(2:end-1)\log((y(2:end-1)-c1)/c2)
The anonymous function is defined for predefined column vector x as
res = @(c) c(1)+c(2)*exp(c(3)*x) - y;
and the call, say
[x,ssq,cnt] = LMFnlsq(res,[c1;c2;c3],'Display',1);
Provided w=0 (without errors), x=(0:.1:2)' with c=[1,2,-1], the initial estimates of unknown coefficients are
c0 = [1.2707, 1.7293, -1.6226].
The call
[c,ssq,cnt] = LMFnlsq(res,[c1,c2,c3])
gives the exact solution c=[1, 2, -1] in 9 iterations.

Notes:
* Users having old MATLAB versions earlier than 7, which has no anonymous functions implemented, have to call LMFnlsq with named functions for evaluation of residuals.
* See LMFnlsqtest for a short explanation and solved examples.
* If values of unknowns are different in orders, it is recommended to use variable scaling (see help),
* The script BoxBOD.m from NIST testing problems is complemented with the solution.

Reference:
Fletcher, R., (1971): A Modified Marquardt Subroutine for Nonlinear Least Squares. Rpt. AERE-R 6799, Harwell