Why Is Non Linear Least Squares Method from MATLAB and Alglib Gives Different Results on the Same Data?

Question

i'm trying to rewrite my Matalab prototype for some DSP to C++ and encountering a displeasing problem. I'm trying to fit data to a function $y = a * (\pi / 2 + arctg(b * x))$. In Matlab it works well with the following code:

fo = fitoptions('Method', 'NonlinearLeastSquares',...
                'Lower', [0, 0],...
                'Upper', [Inf, max(x)],...
                'StartPoint', [1, 1]);
[curve2, gof2] = fit(x,y,ft);

In curve2 I get a and b coefficients, which are pretty well. I try to reproduce this in C++ with Alglib's lsfit using Levenberg–Marquardt algorithm:

alglib::real_1d_array y, c;
alglib::real_2d_array x;
... // filling x and y arrays with data
c = "[1.0, 1.0]";
alglib::ae_int_t maxits = 0;
alglib::ae_int_t info;
alglib::lsfitstate state;
alglib::lsfitreport rep;
alglib::lsfitcreatef(x, y, c, 1e-4, state);
alglib::lsfitsetcond(state, 1e-5, maxits);
lsfitsetbc(state, "[0.0, 0.0]", "[+inf, 5.0]");
alglib::lsfitfit(state, func);
lsfitresults(state, info, c, rep);

At the output (c) I get completely different coefficeints which are doesn't have any correlation with the same from Matlab.

x and y are the same in Matlab and C++ program. The only difference I see is the optimization algorithm: in Matlab I use trust-region method when in C++ I use Levenberg–Marquardt algorithm.

Could you explain me that strange behaviour?

Answer 1

When you solve Non Linear Least Squares problem of a non convex cost function the end solution (Which is guaranteed to be a Local Minimum) will depend on:

1. Method of Minimization.
2. Method Parameters.
3. Starting Point.

In the case above you set the starting point to be the same for both.
Yet in AlgLib you use the method of Levenberg Marquardt (Classic for Non Linear Least Squares).
MATLAB used to use Levenberg Marquardt as its default in the past.
Yet in recent versions it uses more modern method called Trust Region.
The trust region based methods limit their step size to be more conservative.
Basically like the LM they approximate the problem using Quadratic Function yet in each step the find the minimum with constrain about the domain (The Trusted Region) hence usually ends in the closest local minima while LM might skip it.

If you share the values of a and b in a single simulation we'll be able to look into the Cost Function (As it depends on 2 parameters).

If you are after more modern methods to solve the problem in C (Which I recommend), you should look for Trust Region solvers:

• Gnu Scientific Library (GSL) - Non Linear Least Squares Fitting.
• Ceres Solver - Solving Non Linear Least Squares.
• OptimLib - Conjugate Gradient (It's not optimized for Non Linear Least Squares, yet can be used).

If you look in GitHub for Non Linear Least Squares and Trust Region you'll find more.

Few tips:

• The problem isn't convex hence a good strategy is perturbate the starting point. This is a poor man Global Optimization strategy.
• If you're sensitive to the result (Must have the global minimizer) you should start with a rough grid search to get better starting points.
• At top level if needed use full Global Optimization framework like you have in MATLAB.

Answer 2

Many model fitting problems don't have unique solutions. One reason is that the underlying minimization isn't convex or they use different regularizations. MATLAB also usually has a number of options for each call. This is a common problem in feature extraction.

The MATALB doc interface often has a few reference papers near the bottom of the page for each function. You should look there.

If you really need MATALB and C++ to agree, roll your own code. Don't use libraries. You can also call MATLAB fro C++ or the other way around. This doesn't really address uniqueness.

You should perturb your problem a little bit and see how much change you get in your solutions. If you don't get good agreement, within MATLAB or C++, the problem is probably ill conditioned.

Another diagnostic would be the final solution from one as a starting point for the other. If either converged quickly to the initial starting point, multiple convex minima would explain your results.