Levenberg-Marquandt code implementation

2014 1 28
1 个回答
回答已采纳
更新时间：2014 1 29
1 次查看（30 天）

0 个投票

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?

请先登录，再回答此问题。

 采纳的回答

Matt J
Matt J 2014-1-28
编辑：Matt J 2014-1-28

2 个投票

Because of missing parentheses, I imagine
step = (Hessian+diag(lambda)) \ (Jacobian'*Residual)

5 个评论
显示 3更早的评论 隐藏 3更早的评论

RiccardoTisseur
RiccardoTisseur 2014-1-28
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
Matt J 2014-1-28
编辑：Matt J 2014-1-28
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
RiccardoTisseur 2014-1-29
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
Matt J 2014-1-29
编辑：Matt J 2014-1-29
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
RiccardoTisseur 2014-1-29
thank you very much for the help

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by