Comparison between MKL Nonlinear Least Squares and Matlab

[Schimi]

Hi,

I work on a project that needs a nonlinear least squares solver. So far I worked in Matlab, but the final code needs to run in C/C++ (but could be also Fortran if necessary). Therefore I would like to replace the Matlab solver with the MKL one. 

Generally both algorithms should be similar, as they are based on the trust region algorithm, but implementation wise they are different. I created an working solution with the provided Intel example that at a first glance worked very well. After some testing though, I found out that MKL gets more unstable when the initial guesses are further away from the best solution. Matlab instead still gets to the best solution with wide varying start points and boundaries.

I have attached an example that already shows my problem (see figure at the bottom) and is based on the census example in Matlab. I created two projects, with the first creating an executable that spits outs the parameters to solve the equation y=a*(x-b)^2 or a .mex file project can be used that can be called from Matlab to have a direct comparison.

As I'm new to the MKL libraries I might miss some important considerations or tips that I need to take into account. Moreover the implementation might be not 100 % correct and I hope you can have a look and guide me in the right direction. 

Thanks,

Stefan

[mecej4]

I posted a reply yesterday (11/21/2023); when that did not appear in the thread after several minutes, I posted again (slight differences in wording). This morning, I did not see either reply at 8:00 AM CST (USA). An hour and a half later, both messages appeared.

I have replaced my first reply with this explanatory note. For the OP, the next reply should suffice.

[mecej4]

The expression that you are trying to fit has a defect: it gives y = 0 when x = b, but the data does not match that, as you can see in the portion of your plot around year 1800, where the data points deviate from the fitted curve.

A second objection that I should raise is that there is no need to approach this fitting problem as one in the class of nonlinear least squares problems.

Both objections can be answered by using a polynomial, y = a.x^2 + b.x + c, and using a linear least squares method. Linear regression problems can be solved without needing to provide initial guesses for the regression coefficients.

Here is Fortran code to solve the linear regression problem using the Lapack routine DGELS, which is included in MKL.

program qfit
   implicit none
   integer ndat,i,ndeg,lwork,info
   double precision, allocatable :: yr(:), pop(:), popf(:), Amat(:,:),c(:),work(:)
   
   ndat = 21
   ndeg = 2
   lwork = (ndeg+1)*ndat*2
   allocate(yr(ndat),pop(ndat),popf(ndat),AMAT(ndat,ndeg+1),c(ndeg+1),work(lwork))
   
   yr = [(1790.0d0+i*10.0d0,i=0,ndat-1)]
   pop = [3.9,5.3,7.2,9.6,12.9,17.1,23.1,31.4,38.6,50.2,62.9,76.0,92.0,105.7,122.8,131.7,150.7,179.0,205.0,226.5,248.7]
   
   AMAT(:,1) = 1.0d0
   AMAT(:,2) = yr
   AMAT(:,3) = yr*yr
   
   call dgels('N',ndat,ndeg+1,1,AMAT,ndat,pop,ndat,work,lwork,info)
   
   print '(A,3es12.5)','Coefficients: ',pop(1:3)
end program

The results agree quite well with the results from the following Matlab script.

yr=1790:10:1990;
pop = [3.9,5.3,7.2,9.6,12.9,17.1,23.1,31.4,38.6,50.2,62.9,76.0,92.0,105.7,122.8,131.7,150.7,179.0,205.0,226.5,248.7];
c=polyfit(yr,pop,2);
popf = polyval(c,yr); 
c
plot(yr,pop,'o',yr,popf,'r')

[Schimi]

Hi @mecej4 ,

at first thanks for your reply!

This example can be solved with a polynomial equation true, but it should be only something that shows my problem with the nonlinear solver. Actually I try to fit a generalized logistic equation through data points and the MKL results are much more unstable than the Matlab results. Please have a look on the picture below. The dashed lines represent the MKL results and the solid lines the fits from Matlab. In the optimum case they should follow at the left top chart each other and are close to 0. The fits are okay I would say from the R² value, but my determination of the inversion parameters aren't that robust. That's why I showed the simple example above, as I have already seen it there and are not allowed to post the data that I use for my project.

I tried around with the parameters that I can set, but I haven't found a solution that suffice what I want to achieve. So it is more a question about to setup the MKL solver correctly or to get some tips, tricks, considerations to achieve similar results as Matlab is offering.

Thanks,

Stefan

[Schimi]

Hi @mecej4 ,

thanks for your reply, my one wasn't posted as well yesterday, so I retry now.

You're absolutely right that this problem could be solved with a polynomial approach, but this example was more to show the problem that I actually have with MKL in comparison to Matlab. I haven't posted my real problem yet, but to foster the conversation I will do it now. 

I try to fit a curve that is based on a generalized logistic equation through data points. This means I need to fit a curve that is described by 4 or 5 parameters, which needs to be constrained by an upper and lower limit. The figure below shows the comparison between the Matlab (solid line) and MKL (dashed line) results. Ideally both curves should be on top of each other and close to 0 on the top left chart (I know this isn't the case for the left hand samples, but this is another topic and can be neglected). You can see that the MKL solver is more unstable and does not find every time the best fit, which can be seen from the R² values. 

Since I'm a bit tied up from my company to share this data with you, I have chosen the simple Matlab example to illustrate my problem that setting the same initial values and boundaries can lead to different results between the two solvers. My questions is therefore more if there are some tips and tricks, comments, considerations, etc. that I need to take into account to achieve similar results with MKL? 

Hope this makes my situation a bit clearer .

Thanks,

Stefan

[mecej4]

If you cannot provide a test case (one for which you provide code and data with no restrictions on exchange and discussion), there is little to write about. Your zip file contains about 64 MB of files, and that is far beyond what is reasonable for a test case, especially for a rather simple NLS problem.

All that we need to see is (a) the expression(s) to fit (b) the x-y data to fit (c) bounds on the coefficients, if any and (d) starting values of the coefficients.

A nonlinear least squares problem may have multiple solutions. Which of those solutions is found can depend on the algorithm as well as the initial guesses that are provided.

I have used the MKL trust-region solver in the past. When applied to the NIST NLS test problems, the (unconstrained) solver worked very well. For the same test problem set, there is Matlab code on the Matlab File Exchange site.

[Schimi]

Sorry, I forgot to tidy up the folders by removing the unnecessary compilation files. Please find the condensed version attached.

Indeed, I don't expect a unique solution, but when I have a look at the difference in the output of my first post between Matlab and MKL for the same initial guess and limits, there must be something in the solver that I need to trigger to avoid landing in a local minima that doesn't lead to a fit at all.

Please have a look CensusTest.m file in the MatlabCompare folder that does the comparison between the Matlab fit and the MKL one. I have created three fits, a second order polynomial,  one with y=a*(x-b)^2 and another one with y=a*(x-b)^3. For all fits I see that if I move away from initial values that are close to the fit that Matlab provides, I get problems with MKL to have an acceptable fit at all. 

Maybe I do something wrong and therefore I have added my Mex project to call the MKL solver in Matlab and an .exe project to test it directly in C++. Again, I know it is a bit unlucky that I cannot share my real problem with the community, but as this simple Matlab problem already shows differences that are major, It would be nice to understand if and how I can setup the MKL solver to solve or mitigate these differences.

I hope this makes it clearer, if not I can try to explain it even further.

Thanks,

Stefan

[mecej4]

Your file  NLSTestCensus.c is almost "there". The only defect is that your guesses for X[0] and X[1] are too far from the solution. You have chosen 1.0 for both. By looking at your graph above, we can make a better guess for X[1] = 1800 (French Revolution just completed). We can then take the coordinates of the rightmost point on the graph and use the expression being fitted to solve for X[0] = 0.0065. Now put these estimates in your code and run.

The output after making these changes (and changing the bounds as described below):

S:\LANG\MKL>NLSTestCensus.exe
Iterations : 3
Final residual : 1.568930e+01
Stop-criteria : 6
Parameter a : 6.092417e-03
Parameter b : 1.788581e+03

This agrees with the  results that one obtains from the Matlab code.

Your lower and upper bounds are so wide as to be immaterial. If you wish, you can reduce the bounds on X[0] to (0.005, 0.01) and the bounds on X[1] to (1500, 2000).

It may be worthwhile to do a preliminary NLS fit without bounds if there are no hard bounds in the specification of the problem. In a subsequent run, you can add in bounds guided by the results from the unconstrained run.

As far as hooking up C, MKL and Matlab is concerned, I think that it would be quite counterproductive and difficult to debug until the user has developed a feel for what are reasonable estimates of the fit coefficients. Get the pieces to work separately, and then, if you wish, hook them together.

[ShanmukhS_Intel]

Hi Stefan,

Thanks for posting in Intel Communities.

Has the information provided earlier helped? Could you please let us know if you have any updates on your issue?

Best Regards,

Shanmukh.SS

[Schimi · Accepted solution]

Sorry I was caught by a cold and not able to work for some days.

My colleague and I had a look in the referenced book from the MKL library describing the trust region algorithm and found a solution that at least show closer result for our problem.

We did a normalization on the initial data and the upper an lower limits to have values that are close to 1. Moreover we tried to have lower and upper bounds which are separated by 1 to have more control with the RS value. 

Anyway, thanks for your reply's @mecej4, it just showed that you need to be certain what you do and that you need to understand the sensitivity of the algorithm.

Thanks,

Stefan

[ShanmukhS_Intel]

Hi,

It’s great to know that the issue has been resolved, in case you run into any other issues please feel free to create a new thread.

Best Regards,

Shanmukh.SS