Levenberg-Marquandt code implementation

Hi, I'm using the Levenberg-Marquandt algorithm to solve a non-linear equations system. Using fsolve() in debug mode, i follow step by step the matlab implementation, but I didn't find what I have expected. The implemented algorithm calculate the levenberg step with the following formula

AugJacobian=[Jacobian;diag(sqrt(lambda))]; AugResidual=[Residual;zeros];
step=AugJac \ AugRes;

In theory the Levenberg step is calculated with

step = (Hessian+diag(lambda)) \ (Jacobian'*Residual)

The results of the two formulas are different.

Why the implemented algorithm is not the one described in theory http://www.mathworks.it/it/help/optim/ug/least-squares-model-fitting-algorithms.html ?

or I'm missing something?

Answer 1

[Matt J]

Because of missing parentheses, I imagine

 step = (Hessian+diag(lambda)) \ (Jacobian'*Residual)

[RiccardoTisseur]

I'm sorry, probably I haven't explained well my problem. I like to know why in fsolve the lenveberg step is calculated with

AugJacobian=[Jacobian;diag(sqrt(lambda))]; AugResidual=[Residual;zeros];
step=AugJac \ AugRes;

and not with strictly the theoretical formula

step = (Hessian+diag(lambda)) \ (Jacobian'*Residual)

[Matt J]

Theoretically, the two are equivalent (I assume you understand why), but numerically, the latter is not as well conditioned.

As an example, consider the following A*x=b system, the solution to which is x=[1;1]

>> A=[rand(2); eye(2)]; A(1)=1e6; b=sum(A,2);

Using the first solution approach you've shown above, we obtain a certain error

   >> A\b  - [1;1]

ans =

     1.0e-15 *

      0.2220
           0

whereas the second approach is equivalent to the following, which yields much higher error

   >> (A.'*A)\(A.'*b)  - [1;1]

ans =

     1.0e-10 *

     -0.0000
      0.8179

[RiccardoTisseur]

I think you have my solution :-) The problem is that I don't understand why the two are equivalent. The second approach that you have written is the Newton step evaluation

(A'*A)\(A'*b)

what i have expected for the levenberg is more like

((A'*A)+lambda*diag)\(A'*b)

the difference is that in the first the Hessian is approximated with first order terms so equal to J'J, and in the second the Hessian is approximated with an augmented diagonal term -> J'J + lambda*diag .

Thanks in advance for the help

[Matt J]

When A has the form [J; c*D] where D is diagonal, and b has the form [r;zeros] then

 (A.'*A) = J.'*J + c^2*D^2

and

   A.'*b= J.'*r

Thus, (A'*A)\(A'*b) reduces exactly to the expression you're looking for with c^2=lambda and D^2=diag.

[RiccardoTisseur]

thank you very much for the help