4.4

4.4 | 17 ratings Rate this file 126 Downloads (last 30 days) File Size: 858.31 KB File ID: #17534
image thumbnail

LMFnlsq - Solution of nonlinear least squares

by Miroslav Balda

 

15 Nov 2007 (Updated 20 Dec 2011)

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

| Watch this File

File Information
Description

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

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
LMFsolve.m: Levenberg-Marquardt-Fletcher algorithm for nonlinear least squares problems
This submission has inspired the following:
LMFsolve.m: Levenberg-Marquardt-Fletcher algorithm for nonlinear least squares problems

MATLAB release MATLAB 7.3 (R2006b)
Other requirements The testing script LMFnlsqtest exploids functions inp, fig and separator from FEX. Inform me on bugs in LMFnlsq or possible improvements by e-mail.
Tags for This File  
Everyone's Tags
Tags I've Applied
Add New Tags Please login to tag files.
Comments and Ratings (25)
07 Jul 2008 Ender Wiggin

The solution does not work. I get an error whenever I copy and paste the code than run it on MATLAB

08 Jul 2008 Miroslav Balda

I tested the function LMFnlsq once again after reading your comments and an evaluation. I am sorry with you that you are not able to copy and paste ascii lines without errors. You did not try the script LMFtest, which is included in the zip file together with LMFnlsq, where you could get results without any copying. Was it sou difficult to look at it and find how to use the function LMFnlsq. I am afraid of the reaction of other people, who could use even better function than that they may find in the Optimization Toolbox, but seeing your evaluation they will pass it. It is not fair. I thing that you evaluated you yourself.

15 Jul 2008 Jon Jackson

Works well, thanks.

For others who may be looking for the LMFtest dependencies, they are file ID numbers 9033, 9035 and 11725 (as documented in the pdf file appendix).

24 Sep 2008 Levente Hunyadi  
25 Sep 2008 C Schwalm  
28 Oct 2008 Steven White

This is a great bit of code that does exactly what it says. I used it to replace lsqcurvefit as I no longer have the optimisation toolbox.

I would suggest that the author also considers using additional examples, as the ones given are very complicated to understand if you come from a very different field. However, after a little search I found this by the author, which I found much more useful:
http://mathforum.org/kb/message.jspa?messageID=6283753&tstart=0

It tells you how to use the function to fit a curve to some data. This should be easily modified for fitting other functions/data etc.

27 Jan 2009 Polarman  
24 May 2009 Kevin

Professor Blada:
I truly apprecate you helping me out to solve nonlinear equation.
Your code worked beautifully for me and thanks again for your prompt feedback in weekend time.
Thanks again, Dr. Balda.

03 Jun 2009 jeyasenthil

i try to run the Example Rosenbrock's function ( SQP)with "xy" in the output arguments ..but it not work ...it shows xy =[ ] empty matrix...can any body clear this to me?...

function ros = rosen(x)
 d = sqrt(x(1)^2+x(2)^2)-1.5;
 r=1.5;
  ros = [ 10*(x(2)-x(1)^2) % first part, f1(x)
1-x(1) % second part, f2(x)
(d>0)*d*1000 ];

 
x0=[-1.9,2];
>> options=LMFnlsq;
>> [x,ssq,cnt,nfJ,XY] = LMFnlsq('rosen',x0,options)

03 Jun 2009 Hossein SOLEIMANI

Hi

how can I give an interval for the results, i don't need negative values for answers or I need some answers higher
than some values,

you also helped me before Mira

thanks.

04 Jun 2009 Miroslav Balda

This part of the submission is intended for Comments and Ratings. More over, I pleased anybody to write me by e-mail on errors and improvements, see the section Other requirements. This is the reason why I'll not communicate here anymore.

11 Jun 2009 thomas  
15 Aug 2009 Jon Doan

Very nice

01 Oct 2009 Shin-ichi Uehara

Hi,

I want to do curve fitting of some laboratory results with a model.
The model which I want to use can not explained by a single function,
but is explained by a series of differential equations.
In this case, I do not know how to apply LMFnlsq.
It is very helpful if you give any advice.

Thanks.

(I am applying "fminsearch" to this problem so far, but
the description of "fminsearch" mentions that this code is
not suitable to this kind of problems.)

12 Oct 2009 Miroslav Balda

Read my answer dated 04 jun 2009. I'll answer you later by e-mail.

09 Nov 2010 Neil Dalchau

This is an excellent improvement to LMFsolve, which I only became aware of today. This upgrade has considerably improved some of my existing analyses. Convergence is now 100% for my problem.

Many thanks Miroslav!

16 Feb 2011 Xylar Asay-Davis

This function was extremely helpful for the problem I'm working on: least-squares fit to multidimensional rational polynomials.

Thanks very much!

09 Mar 2011 Ryan Muir

This program is excellent. I found that it was able to do my nonlinear regression very quickly, with multiple quick and easy options for both analysis and displaying the data correctly. Further, once you figure out how to use the function, performing a nonlinear regression was quite easy.

However, trying to initially figure out how to use the software was quite hard =P. I initially attempted to follow the instructions posted here, which are not complete or correct in its explanation or equations. Once I discovered the included test file, which is correct and will run, I was able to copy and paste the relevant section of code (example 3) into matlab and get instant results. I was then able to pick apart and transform the example for my own purposes. Figuring all that out took quite a while. I took off one star for that upset, but it's really an excellent software package!

23 Jul 2011 Ryan Muir

Works fast and well! However, it has a HUGE learning curve if you don't write code at a very high level...the included examples I found to be less than completely understandable to me, and don't work unless you know how to use them. Instead, here is a VERY SIMPLE program that uses this algorithm to fit a function to a simulated dataset. I took it and modified it from a forum post elsewhere on the internets. Recall that chi^2 = sum(((x - mu)/sigma).^2), and the point of this function is to determine what parameters in mu minimize chi^2:

% Ender 2008-07-08
a = 10, b = -8, c = -.2 %Model parameters
t = (0:.01:1)'; %should be a column vector

data = a + b*(1-exp(t/c)) + .1*randn(size(t)); %Function plus random error; simulates data

res = @(x) data - (x(1) + x(2)*(1-exp(t/x(3))) ); %Numerator of chi^2, with quantity squared implied for fitting in LMFnlsq; also known as residuals

x0 = [5,3,-.1]; %rather poor initial guess of parameters

[x,ssq,cnt] = LMFnlsq(res,x0) %Returns the parameter values which minimized the value of the numerator

%Plot the data
hold on
plot(t,data)
plot(t,(x(1) + x(2)*(1-exp(t/x(3)))),'r'), grid

02 Dec 2011 David Zhang

This doesn't work for all the least square tests established by the National Institute of Standard and Technology (NIST).

http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml

So this can't be called the industry standard.

02 Dec 2011 David Zhang

By the way, the failures occurs in

MGH17
BoxBOD
MGH10
Eckerle4
Bennett5

14 Dec 2011 Miroslav Balda

@ David:
I thank you for thorough testing of my function LMFnlsq. You have found that 5 functions of the highest difficulty out of total 27 failed. It is true only partially. They failed only for the most bad initial estimates of the solution. The good solution has been obtained with the more close estimate for all above examples but Bennett5. In this case only one digit of sum of squares was valid, however, when sum of squares dropped by 8 orders! All those results have been obtained with embedded function for finite difference Jacobian matrix estimates, what strongly diminishes precission of results. Only analytical Jacobian matrix ensures top result. I am affraid that you also used finite difference approximates.

The NIST collection of testing problems is only a recommendation for software authors. In contrary to your opinion, I am satisfied by the results I have obtained. Especially, after testing a less difficult problem Hahn1, which converged for both estimates of the solution, while the function lsqcurvefit from the Optimization Toolbox (OT) did not find any solution. Would you evaluate it also by one star? Not me, because OT is much better. However, my function is smaler and takes fewer iterations to solve the problem. This is the reason, why I regard your evaluation as a bad joke.
One star means: "withdraw your function from FEX!". I will not obey your evaluation, because there are hundreds and may be thousands of people who downloaded and successfully applied my function for solving their problems. Some of them gave the function the top evaluation, even that I never did not declare it as an "industry standard".

14 Dec 2011 David Zhang

@Miroslav:

I was using analytical Jacobian. I also will not obey your evaluation that a failed lsqcurvefit does not deserve a 1 star. The NIST solutions are tested using different methods (hence the standard deviation is also available), so failing even one of them means there are codes out there that performs better, even using the poor starting condition.

Furthermore, 3 of the 5 failure cases are observed data. They are not some cooked up example designed to cause the optimization to fail. Therefore, if used in real world examples, there is a possibility that it will not work. I am simply letting everybody know a caveat cantor that the program solution should be treated with a grain of salt, even if the sum of squares dropped by 8 orders. Although as the website says, passing all tests is not a pre-requisite to be a good optimization software, but don't you think not passing the tests from an agency that establishes industry standards would not be considered as "industry standard"?

I admit giving it 1 star is too aggressive. I do this to get your attention to respond. Indeed it worked. Unfortunately I can't go back and change it. Because this program works for all easy to intermediate difficulty problems, it deserves a 4 star.

17 Dec 2011 Shayan

Man thank you so much, it was amazing!

02 Jan 2012 Sadik

I thank you so much for providing this code. Implementation has its own challenges but your code worked like a charm to show how beautifully theory and practice meet. With good initial estimates, the algorithm converged in only a few iterations.

FYI: I used it for T1 estimation in MR where the equations in our case were much more complicated than single exponential recovery equations.

Thanks.

Please login to add a comment or rating.
Updates
15 Nov 2007

Formal changes of the description for better outlook.

06 Dec 2007

Improved setting of Options, added an example in LMFtest.m and description and help.

09 Jan 2009

Improved functionality, removed bug in assembling analytical Jacobian matrix, complemented example for curve fitting.

27 Jan 2009

Implemented new subfunction for printing intermediate results, modifiedcation both the function and a testing script code. The description has been complemented by a description of the LMF method.

20 Dec 2011

Improved description (help), improved printing of results, complemented results of NIST test by BoxBOD problem.

20 Dec 2011

Forgotten zip-file of functions appended

Tag Activity for this File
Tag Applied By Date/Time
optimization Miroslav Balda 22 Oct 2008 09:35:12
levenberg Miroslav Balda 22 Oct 2008 09:35:12
marquardt Miroslav Balda 22 Oct 2008 09:35:12
least squares Miroslav Balda 22 Oct 2008 09:35:12
fletcher Miroslav Balda 22 Oct 2008 09:35:12
optimization Cristina McIntire 28 Jan 2009 10:21:29
fig and separator Jack cao 31 Mar 2011 04:30:18

Contact us at files@mathworks.com