forum.alglib.net

ALGLIB forum
It is currently Sun Nov 18, 2018 6:41 am

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  [ 2 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: 833
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  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 2 posts ] 

All times are UTC


Who is online

Users browsing this forum: Bing [Bot] and 7 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