Infinite loop in GSL Levenberg-Marquardt fitting

ROOT
1

The question will be hard to answer, for I cannot prepare a small piece of code with the problem. I decided to give GSL fitting methods a try, however they are rather not documented in root docs, so I am not sure if I do everything right. So in my code, which works fine with minuit, I make:

ROOT::Math::MinimizerOptions::SetDefaultMinimizer("GSLMultiFit", ""); proj->Fit(psf_fin, "RV");

everything seems to be initialized properly. In the first step, I just look for proper centre, scale and background for my 2D function. Works fine. Than I release one parameter responsible for shape change, and very quickly I get an infinite loop with fitter trying just to values of the parameter, not changing other parameters. This is my output:

Change in parameter 6. Old value: -8.41e-05, new value: -8.40963430686444e-05 other: 0 39.3975452251773 16.0088352035183 0.0653320018192862
Recalculating profile...
Warning in <TFile::Append>: Replacing existing TH1: psf_histo (Potential memory leak).
-8.40963430686444e-05
Change in parameter 6. Old value: -8.40963430686444e-05, new value: -8.41e-05 other: 0 39.3975452251773 16.0155785493089 0.0653320018192862
Recalculating profile...

Two values of parameter on which it loops change with the limits I set for the parameter. I understand that you cannot examine the problem without a small code snippet, but maybe you can answer the following questions:

  1. Am I switching to GSL and using it in a proper way?
  2. Is this algorithm is confirmed to be working in root?
  3. Did you ever encounter similar problem and can give me a clue?

Thanks for any answer.

2

Hi,

you are using the algorithm in the proper way and it is supposed to work.
The Levenberg-Marquardt algorithm is a good and efficient algorithm for non linear least square fitting. It works well and it is very fast if you are not too far away from the solution.
Your problem might be caused by a numerical error in the derivative calculation, but I would need to understand it better.
Normally it cannot go in infinite loop because there is a maximum number of iteration set to something like ~ 100. I would eventually decrease this value with ROOT::Math::MinimizerOPtions::SetDefaultMaxIterations(n).

Could you send me also the output you get from the GSL fitter when running with option “V” ?

Best Regards

Lorenzo

3

[quote=“moneta”]Hi,

you are using the algorithm in the proper way and it is supposed to work.
The Levenberg-Marquardt algorithm is a good and efficient algorithm for non linear least square fitting. It works well and it is very fast if you are not too far away from the solution.
[/quote]

Hehe, according to wikipedia, this algorithm: “The LMA is more robust than the GNA, which means that in many cases it finds a solution even if it starts very far off the final minimum.”

However, it’s not the first time when wiki is wrong. Anyway, I just wanted to try and perhaps if there is a bug, help to find it.

[quote]
Your problem might be caused by a numerical error in the derivative calculation, but I would need to understand it better.
Normally it cannot go in infinite loop because there is a maximum number of iteration set to something like ~ 100. I would eventually decrease this value with ROOT::Math::MinimizerOPtions::SetDefaultMaxIterations(n).

Could you send me also the output you get from the GSL fitter when running with option “V” ? [/quote]

Indeed, there may be a problem in derivative or a kind of a bug in function itself, a bug which does not (at least visibly) affect other fitting methods. The function is quite complicated (numerical integration, than numerical convolution, everything with a finite number of points) and there may be some special, unstable cases. However, I have doubts that any fitting algorithm should behave like that, which I describe below.

The output when fitting just 2D centre position, scaling parameter and constant background is like that:

<TCanvas::MakeDefCanvas>: created default TCanvas with name c1
Fitting!
Minimize using GSLNLSMinimizer lmsder
----------> Iteration 0 / 100 status success
            FVAL = 435.719635354222021
            X Values :  x = 0 y = 600 dx = 0 dy = 40.8317670854359 z = 0.085 ast = 0 coma = -8.41e-05 spher = 0 spher1 = -0 scale = 37.9333080069766 bg = 0.0495324967614449
          after Gradient and Delta tests:  the iteration has not converged yet

than, after a number of iterations it does not converge (which is strange but normal for my function for any fitting algorithm), but finishes with a solution. However, than I start a next fit, releasing one of the shape parameters, and it does not even come to “Minimize using GSLNLSMinimizer lmsder” message. It just, depending on the limits, starts with just 2 parameter values which it changes one to another in a loop, or quickly comes to this stage (I know that from the output from the fitted function). There is no output from the fitter (fcourse with “V” option). I guess the loop is infinite - I had it running for a few hours without change:

<TCanvas::MakeDefCanvas>: created default TCanvas with name c1
Fitting!
Change in parameter 6. Old value: -8.41e-05, new value: -8.40963452435962e-05 other: 0 41 60 0.0101667202426
Change in parameter 6. Old value: -8.40963452435962e-05, new value: -8.41e-05 other: 0 41 60.01200044998 0.0101667202426
Change in parameter 6. Old value: -8.41e-05, new value: -8.40963452435962e-05 other: 0 41 60 0.0101667202426
Change in parameter 6. Old value: -8.40963452435962e-05, new value: -8.41e-05 other: 0 41 60.01200044998 0.0101667202426

If it helps, I can send the output from the position+scale+background fit which comes to an end. I can also send my code, which is rather long and dirty, but I would not want to waste your time debugging my function… :slight_smile:

4

Hi,

first of all, concerning the remark from the wikipedia, "The LMA is more robust than the GNA, which means that in many cases it finds a solution even if it starts very far off the final minimum."
This is correct if applied to standard Gauss-Newton algorithms for least square problem.

I was referring in particular to another class of algorithms, the variable metric ones (en.wikipedia.org/wiki/Quasi-Newton_method). Minuit, for example, belongs to that category.

Concerning your problem, if I have understood well, you go in infinite loop before getting the message, "Minimize using GSLNLSMinimizer lmsder " ?

This is very strange, no calls to the function should have been made before. If this is the case, it would help a back trace obtained with gdb to understand the call stack of this infinite loop

Best Regards

Lorenzo

5

Yes, exactly that thing happens. To be sure that these are no other calls to function, my code is now:

cerr << "starting fitting" << endl;
		proj->Fit(psf_fin, "RV");

The “starting fitting” message is the last message before the output from the fitted function itself.

I’ll try to obtain proper output from gdb.

6

Hi,

after having looked more in detailed at the code, I see that the function is called before fitting to calculate the initial derivatives at each fit point.
Each time the parameter value is changed by 10^-4, the stepsize which it is at the moment hardcoded inside (I will change this for the next release).

Could it be that for some parameters value and some point the evaluation of your function goes in an infinite loop ?

Lorenzo

7

That is indeed a lot for some of my parameters, but according to what the function prints out, the step of parameter varries…

I am not sure I understand :slight_smile: The printout from my function showing what change of parameter was called, is printed one per function call with this parameter change (when only other point of the function is requested, there is no printout). Therefore the next printout means another function call and that the previous has returned - so there is no infinite loop in the function itself. Is that what you ment?

8

Let me understand better, are you seeing that the function is called an infinite number of times, i.e. you have an infinite printout
or you observe that the time between function calls goes to infinite ?

Lorenzo

9

[quote=“moneta”]Let me understand better, are you seeing that the function is called an infinite number of times, i.e. you have an infinite printout
or you observe that the time between function calls goes to infinite ?
[/quote]

The function is called an infinite number of times, with the shape parameter change just between two values.

10

How many fit points are you having ? The function should be called for every fit point scanning each parameter for two values to calculate the derivatives with respect each parameter at each point.
It should be called npoints * ( npar + 1) times.
I noticed however that in the current version some extra and unneeded call are done for each point. So the total number of function calls done before starting the minimization (before the message is printed) is now
npoints * (npar + 3)

Can you please check if in your case you are calling many more times than this number ?

I will anyway try to remove these unneeded calls for the next version

Best Regards

Lorenzo

11

Your ask for number of calls helped - I think I found where the problem is. It seems, that the fitter chooses one x,y, then performs those two-value changes of all free parameters. Than it takes another x,y, and perform again the same set of changes. I didn’t print the x,y requested before, now I did. So I guess the loop is not really infinite - just very long. Single calculation of this gradient in this order takes number_of_points*all_pair_of_parameter_changes number of function calculations. My function is not fast, takes about 1-2 seconds for each recalculations, so the gradient calculation takes ages.

It seems in other fitting methods the order of changes is different - one set of parameters, request for all x,y, than change of parameters and again request for all x,y function values. This way I calculate the function itself only once for a set of parameter values, which is ~10000 less calculations. Would it be possible to change the order in L-M method too? :slight_smile:

12

In the current interface it is not really possible to change the order. The value of each element of the chi2 (the residual) and its gradient is needed for each point (x,y,…) and for a set of parameter values.
To use L-M (and also for Fumili and Fumili2) the DataElement method of interface FitMethodFunction needs to be implemented
(see project-mathlibs.web.cern.ch/pro … ction.html )

You could probably re-implement this method caching inside the function evaluation. However, I will think about how to improve this interface and
provide a way to do this more transparently.

Thank you for your valuable feedback

Best Regards

Lorenzo

13

Yes, a long time ago I had some kind of caching. Howver I guess that problem is a problem for all cases of some kind of convolution, where to get a point of function you need to have a few points of other, long to calculate function. Without caching one calculates the same points of the other function many times.

Thanks for your answer - I hope that at least those small bugs found in the code during the search for a big unexisitng bug were worth the time spent on it :slight_smile: