Error in Levenberg-Marquardt implementation

[Dirk@LV]

I believe that the implementation of weights in the LM algorithm is erroneous: In particular, in NI_Gmath.lvlib:CNR abx reweighted.vi, the gradient descent steps and the Hessian are calculated with a weighting matrix W build from the data weights w, which is erroneously applied twice for the Hessian. By this, the parameter step sizes from one iteration to the next are wrong as soon as one uses weights, especially if the w's are far away from 1. As a consequence, convergence may become extremly poor even for simple models. And the convergence scales with the scaling of the problem: if you would re-scale your data from meters to millimeters, so that measurement errors (sigma, thus weight=1/sigma^2) change from 10e-3m to 1mm, you could drastically increase your convergence by a factor of 10-100.

In the original implementation, the Hessian is calculated from the Jacobian J as (WJ)^T(WJ) instead of J^TWJ.

For security, I am only posting a screenshot here of a corrected version, because the original VI is in NI-gmath.lib (also I did not connect the Jacobian out wire, which is not used in the Bounded LM VIs but may be used elsewhere). The modified version converges much faster and is independend of the relative scale of the problem.

[Dirk@LV]

I forgot to mention: this is for LV2014 SP1 and older. Did not check if this had been modified in LV2015.

D.

[annav]

Hi,

can you post an example where you can see the convergence difference?

Best,

anna

[Dirk@LV]

Hi Anna,

meanwhile I found that the error seems to be only in the "Constraint Nonlinear Curve Fit", not in the unconstraint. Both versions use a different VI for the gradient step and Hession calculation. The CNR abx reweighted.vi of the constraint fit includes a suspicious "scale factors" input which seems to be left over from the implementation phase and might indicate that this VI has not been properly finished. As can be seen in some other places of the LM implementation, too. See for example this forum discussion:

http://forums.ni.com/t5/LabVIEW/Wrong-math-in-quot-LM-get-new-coefficients-vi-quot-of-Levenberg/m-p/...

Attached you can find a very simple model function, which is just fitting a constant y=f(x)=a to some noisy data points (generated with the same model function VI, but just some extra noise). If you generate data for a parameter a=100, with a standard deviation SDV=1, and initial fit parameter of 120, it takes 30 function calls to find the result.

Now, rescale this to a=1, SDV=0.01, initial parameter=1.2: more than 2000 function calls are required!

I included a corrected version of the VI, too, (with new_ prefix). If you want to use it (on your own risk), you could open the original in the gmath lib and copy/paste the improved section into it.

The unconstraint version seems to work correctly, though.

Dirk

[Henrik_Volkers]

If confirmed, CAR?

[louis57]

HI,

Please can you post a version for LV2013.

Thank you

[Dirk@LV]

Here come the files down-saved for LV13, hope it works.

Dirk

[louis57]

Thank you very much

[DSPmchen]

Hi Dirk,

Thanks for reporting this issue. 

I check this VI. It does contain a bug. However, I think the bug is in the calculation of gradient, not Hessian.

The weight input passes through a Square Root primitive.  So the Hessian calculation in the original implementation is J^TWJ. But the gradient is J^T*sqrt(W)*R, which should be J^TWR.

To fix this issue, a simple solution is like follows - change one wire in calculating the gradient. By doing so, I can pass the test in your attached zip file.

Please let me know if you have any questions or comments.

Best Regards,

Michael

[Henrik_Volkers]

Thanks Michael for confirming this!

Is there a CAR to track the official fix to a new version?