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?

  • $\begingroup$ What is your understanding of both algorithms? Why do you expect them to give you similar results? $\endgroup$ – MBaz May 20 '18 at 19:27
  • $\begingroup$ @MBaz In my view least squares method should give the same result regardless of the optimization algorithm that you use to minimize squared residuals. Am I wrong? $\endgroup$ – Georgiy Manuilov May 20 '18 at 19:38
  • $\begingroup$ Well, they could be converging to different local minima. Are the two solutions "good"? $\endgroup$ – MBaz May 20 '18 at 19:49
  • $\begingroup$ @MBaz Well, Matlab solution is definitely "good" because it corresponds to the physics of my problem, but I can't say that the C++ solution is wrong from the mathematical point of view. However I need to obtain similar to Matlab solution in my C++ code. How can I do that? Using Alglib is not necessary if there is more appropriate solution. $\endgroup$ – Georgiy Manuilov May 20 '18 at 20:00

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:

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.
  • $\begingroup$ Wow, thank you for such a comprehensive answer and great advices. I think the Ceres Solver is the thing I've looked for. I guess I'll try to use it. $\endgroup$ – Georgiy Manuilov Aug 16 '18 at 9:40

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.

  • $\begingroup$ Hi: I can't tell from the code because I am not familar with matlab but, it's generally a good idea to supply the gradients. ( and even the hessian if the alg uses one ). I'm also not familiar with levenberg-marquardt but do have experience with the variable metric algs such as BFGS. If you have access to that one, supply the gradient and see if that helps. $\endgroup$ – mark leeds May 21 '18 at 5:48
  • $\begingroup$ P.S: I also agree with Stanley that, in an ideal world. it's best to roll your own if possible but the convergence criteria piece can get very tricky even rolling your own so that may not help. Using libraries that are well tested is probably a more practical approach. If you have access to R, I highly recommend the Rvvmin package.It's essentially an improved version of BFGS and well tested and written by an expert in the field. $\endgroup$ – mark leeds May 21 '18 at 5:52
  • $\begingroup$ I share your point of you, but I'm not very good at optimization theory to write my own code for algorithms. I think I just try to use Matlab C++ API, which I didn't want to use due to performance issues. $\endgroup$ – Georgiy Manuilov May 21 '18 at 6:16
  • $\begingroup$ @markleeds I've already tried to pass both gradient and Hessian to lsfitfit with my function, but it had no effect. Thanks for advice about BFGS, I guess I'll try it. The problem is that there are not so many well tested libraries for C++ for optimization with understandable API. I've already implemented the other part of my algorithm which uses Fourier and Hilbert transofrms, and some complex numbers mathematics in C++ rewriting it from Matlab and it was piece a cake. So I expected the same with rewriting the next part, but it seems that this is not an easy problem. $\endgroup$ – Georgiy Manuilov May 21 '18 at 6:26
  • $\begingroup$ Hi: If you can rewrite it in R and try Rvmmin and that doesn't work, then, if you send the code to the R-list, I bet someone over there ( I'll read it but don't know if I can help ), might say something useful. $\endgroup$ – mark leeds May 21 '18 at 15:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.