# forum.alglib.net

ALGLIB forum
 It is currently Sat Aug 24, 2019 10:07 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.

 Page 1 of 1 [ 6 posts ]
 Print view Previous topic | Next topic
Author Message
 Post subject: Running Levenberg-Marquardt Non linear algorithmPosted: Mon Sep 10, 2018 3:49 pm

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

Top

 Post subject: Re: Running Levenberg-Marquardt Non linear algorithmPosted: Tue Sep 11, 2018 9:48 am

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

 Post subject: Hidden Bug. Re: Running Levenberg-Marquardt Non linear algorPosted: Sun Jan 27, 2019 8:42 pm

Joined: Mon Nov 20, 2017 11:14 pm
Posts: 16
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

 Post subject: Re: Running Levenberg-Marquardt Non linear algorithmPosted: Mon Feb 18, 2019 9:03 am

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

 Post subject: Re: Running Levenberg-Marquardt Non linear algorithmPosted: Tue Mar 26, 2019 2:22 pm

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
;______________________________________________________________________________
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

 Post subject: Running Levenberg Marquardt Non linear algorithmPosted: Mon Aug 19, 2019 3:05 pm

Joined: Mon Aug 19, 2019 2:09 am
Posts: 4
Location: Bolivia
Matlab had a solver for linear least squares minimisation with linear constraints, which sounds like what you looking for. If you dont have Matlab, then scipy has linear least squares with variable bounds as constraints, which might be too restrictive for you.

CVXpy says that it has an LS solver that recognises linear least squares ptoblems with linear constraints and solve them efficiently, but I have never used that.

Otherwise, just go with a quadratic solver with linear constraintslike cvxopt, those should still be much faster than a general nonlinear solver.

_________________
Games - https://torrent-katochka.com

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 6 posts ]

 All times are UTC

#### Who is online

Users browsing this forum: Brentonunfaivy, Google [Bot] and 1 guest

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

Search for: