Running Levenberg-Marquardt Non linear algorithm

lookinganswers · **Joined:** Mon Sep 10, 2018 3:34 pm **Posts:** 1

We have MATLAB LM (Levenberg-Marquardt) code which needs to be done in C++. I am interested in buying the commercial edition, but before that I would like to see if it is works with our input parameters or not.

To start with, I tried to run simple input values for LM and It failed to run. Error was -8 and the traget function was throwing infinity at 60th iteration - Refer 1 screenshot.
I have attached the problem with sample input values, which I tried to run it.

Please do let me know how to move forward.

Code:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "interpolation.h"

using namespace alglib;
void function_cx_1_func(const real_1d_array &c, const real_1d_array &x, double &func, void *ptr)
{
   // Exponential function of the form: x(1)*exp(x(2)*xdata);
   func = c[0] * exp(c[1] * x[0]);

}

int main(int argc, char **argv)
{
   //
   // In this example we demonstrate exponential fitting
   // by f(x) = exp(-c*x^2)
   // using function value only.
   //
   // Gradient is estimated using combination of numerical differences
   // and secant updates. diffstep variable stores differentiation step 
   // (we have to tell algorithm what step to use).
   //
   real_2d_array x = "[[0.9],[1.5],[13.8],[19.8],[24.1],[28.2],[35.2],[60.3],[74.6],[81.3]]";
   real_1d_array y = "[455.2, 428.6, 124.1 ,67.3, 43.2 ,28.1, 13.1, -0.4, -1.3, -1.5]";
   real_1d_array c = "[100,-1]";
  

   double epsx = 0.000001;
   ae_int_t maxits = 0;
   ae_int_t info;
   lsfitstate state;
   lsfitreport rep;
   double diffstep = 0.0001;

   //
   // Fitting without weights
   //
   lsfitcreatef(x, y, c, diffstep, state);
   lsfitsetcond(state, epsx, maxits);
    alglib::lsfitfit(state, function_cx_1_func);
   lsfitresults(state, info, c, rep);
   printf("%d\n", int(info)); // EXPECTED: 2
   printf("%s\n", c.tostring(1).c_str()); // EXPECTED: [498.8309   -0.1013]
     return 0;
}

Snapshot 1.
Format:
0.000000 -> (100.000000) (-1.000000) (81.300003)
func -> (c[0]) (c[1]) (x[0]))

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

The problem is that exponentials are really unstable - wild derivatives which enforce long steps which in turn result in numerical overflow. I debugged the solver, and it turned out that solver tried to evaluate c[1]=10 which overflowed.

There are exist several solutions to this problem which can be used together:

1. You can solve this problem by putting a limit on step length with lsfitsetstpmax(state, step_length_limit).

In my experiments step limit equal to 1.0 was enough to prevent the algorithm from the destabilization and [498.8309,-0.1013] was returned as result.

However, algorithm had to perform roughly 498 steps to reach this point due to limit on the step length. You can improve situation with performance by specifying per-variable scales, say you can set scaling to [100,1.0]. Thus algorithm will know that it can make long steps for variable #0, but #1 have to be probed carefully. lsfitsetscale(state, "[100,1]") works good enough.

2. Another option is to control what your function returns to the solver and to prevent it from returning too large values. Say, if the value being returned is outside of [-1e20,+1e20] range, it should be clipped. You can do so by adding if( fabs(func)>1.0e20 ) { func=fmax(func,-1.0e20); func=fmin(func,+1.0e20); } to the end of the function.

More tricky target functions may return NANs instead of infinities for "incorrect" inputs, so adding if( !isfinite(func) ) { func=1.0e50; return; } may be necessary.

You can use (1) and (2) together to get better results. Also, after studying your case I feel that it is good idea to modify the solver so it will be able to handle situations like this by automatically decreasing step length. It will be introduced in the next release (scheduled for October/November).

RockBrentwood · **Joined:** Mon Nov 20, 2017 11:14 pm **Posts:** 21

I found a gap in the code which is a hidden bug and might also be related to the problem described above.
The lsfitcreate* routines do not explicitly initialize the optstate field, before minlmcreate* is called on optstate and runs asserts on its fields (in particular on optstate.epsx). The problem is partially masked by your use of memsets for local structure variables. In my copy, I took out all the memsets, so hidden initialization bugs will show themselves more readily.

What I tried at first was to inherit optstate.epsx from epsx in lsfitsetcond(); state->optstate.epsx = state->epsx = epsx. But a simpler and more direct repair is to just reorder the routines in the minlmcreate* constructors from
minlm_lmprepare(), minlmsetacctype(), minlmsetcond()
to
minlmsetacctype(), minlmsetcond(), minlm_lmprepare()

There may be other explicit initializations that need to be done, even when the memsets are used; and similar problems may arise elsewhere, wherever you have nested state structures.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Hi! Sorry for delay with replying :(

I do not understand what is wrong with optstate initialization. This structure is initialized at the end of the every lsfitcreate* call, I just checked - nothing is missing. It is used only in lsfititeration, after correct initialization in the constructor code. So I simply see no place for the error like one you described.

pkuzmic · **Joined:** Tue Mar 26, 2019 2:09 pm **Posts:** 1

I can confirm that the answer posted by Sergey Bochkanov is of course correct. Here is what I got with the software package DynaFit (http://www.biokin.com/dynafit/):

Input script:

Code:

; http://forum.alglib.net/viewtopic.php?f=2&t=3864
; y = c[0] * exp(c[1] * x);
; c = "[100,-1]"; ... initial values
;______________________________________________________________________________
[task]
   task = fit
   data = generic
   algorithm = LM ; means ... Levenberg-Marquardt
[parameters]
   x, c0, c1
[model]
   c0 = 100 ?
   c1 = -1 ?
   y = c0 * exp(c1 * x)
[data]
   variable x
   set data
[output]
   directory ./test/_temp/output/fit-001
[settings]
{Constraints}
   AllParametersConstrained = n
[set:data]
x   y
0.9   455.2
1.5   428.6
13.8   124.1
19.8   67.3
24.1   43.2
28.2   28.1
35.2   13.1
60.3   -0.4
74.6   -1.3
81.3   -1.5
[end]

Best-fit parameters:

c[0] = 498.831 +/- 0.965729
c[1] = -0.101257 +/- 0.000462374

The "+/-" values are formal standard deviations.

DynaFit [1,2] is free to all academic users (in binary form).

References

[1] Kuzmic, P., DynaFit - A software package for enzymology, Meth. Enzymol. 467 (2009) 247-280.
[2] Kuzmic, P., Program DYNAFIT for the analysis of enzyme kinetic data: Application to HIV proteinase, Anal. Biochem. 237 (1996) 260-273.

forum.alglib.net

Forum rules

Running Levenberg-Marquardt Non linear algorithm

Who is online