forum.alglib.net

ALGLIB forum
It is currently Thu Mar 28, 2024 1:55 pm

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 5 posts ] 
Author Message
 Post subject: Running Levenberg-Marquardt Non linear algorithm
PostPosted: Mon Sep 10, 2018 3:49 pm 
Offline

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]))

Image


Top
 Profile  
 
 Post subject: Re: Running Levenberg-Marquardt Non linear algorithm
PostPosted: Tue Sep 11, 2018 9:48 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
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).


Top
 Profile  
 
 Post subject: Hidden Bug. Re: Running Levenberg-Marquardt Non linear algor
PostPosted: Sun Jan 27, 2019 8:42 pm 
Offline

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.


Top
 Profile  
 
 Post subject: Re: Running Levenberg-Marquardt Non linear algorithm
PostPosted: Mon Feb 18, 2019 9:03 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 903
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.


Top
 Profile  
 
 Post subject: Re: Running Levenberg-Marquardt Non linear algorithm
PostPosted: Tue Mar 26, 2019 2:22 pm 
Offline

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.


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 5 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 50 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group