Polyfit polynomial badly conditioned on Linux whereas on Mac no warnings

I have a code which has to interpolate some simple data grid using polyfit.

However, when I run the code on my computer (MacBook, Matlab 2021a) there are no warnings or errors. Whereas, when I run the same identical code on the server I work with (which runs on Linux, Matlab 2019b) Matlab gives me warnings about the polynomial that is "badly conditioned".

I have also checked one-by-one the results before the polyfit and they are identical between the two systems.

Do you have any idea?

[Jan]

Without seeing the data it would be pure guessing.

If the polynomial is badly conditioned, why don't to improve the condition using the 3rd output of polyfit?

[EDITED] I wrote "2nd output", but for normalizing it is the 3rd output.

[Torsten]

Matlab 2019b and Matlab 2021a are different versions - so there should be some differences.

[Luigi De Maria]

@Jan the data are coordinates with respect to time instant (3 positions, 3 velocities).

The data to be interpolated are:

x = [-0.4297; -0.2478; -0.0659; 0.1160] * 1.0e-08

Interpolation nodes:

tspan = [0.0537; 0.0786; 0.1035; 0.1284] *1.0e-03

The interpolation is called as:

p1 = polyfit(tspan, x, 3);

For the second output we have:

df = 0

normr = 1.2451e-24

[Walter Roberson]

The data is probably not what you have posted. You have posted the "format short" version of the data, which is only 4 significant figures. You should

fprintf('%.16g ', x); fprintf('\n');
fprintf('%.16g ', tspan); fprintf('\n');

to show us the exact data.

[Luigi De Maria]

@Walter Roberson I am so sorry, I am in a bit of a hurry due to a deadline and I didn't realize I posted short format.

Here are data in the format you required:

x = [-4.297441637883469e-09; -2.478438427644686e-09; -6.594352312665492e-10; 1.159567962615356e-09]

tspan = [-3.73762709838521e-05; -1.245875699462662e-05; 1.245875699462662e-05; 3.73762709838521e-05]

[Walter Roberson]

What shows up for

version('-blas')

on the two systems?

Chances are high that the newer version of MATLAB uses a newer MKL

[Luigi De Maria]

@Walter Roberson the command gives:

Matlab 2019b: Intel(R) Math Kernel Library Version 2018.0.3 Product Build 20180406 for Intel(R) 64 architecture applications, CNR branch auto

Matlab 2021a: Intel(R) Math Kernel Library Version 2019.0.3 Product Build 20190125 for Intel(R) 64 architecture applications, CRN branch AVX2

[Luigi De Maria]

As a conclusion, if you experience the same problem, consider to call all three outputs of polyfit as adviced in the documentation. This will automatically introduce a proper scaling and recentering of the nodes.

Answer 1

[Jan]

The warning is serious. Although the results look find in your case, scaling the input is the way to go for a scientifically stable solution. It is cheap and easy: If the 3rd output MU is requested from polyfit, x is scaled by:

MU = [mean(x), std(x)]
x  = (x - MU(1)) / MU(2);

The same is applied in polyval also before running the horner scheme.

Scaling is important in numerical maths to control the effects of the limited precision. This concerns leats square fits, integrations and optimizations. The argument "the results are okay so I do not need scaling" or "I suppress the warnings" is a bad decision and should not appear in serious scientific work.

In numerical maths the final value can be called "result" only, if the accuracy is controlled. Otherwise the program is a weak random number generator with a certian probability < 1.0 to reply a final value near to the expected one. I've seen to many PhD theses, in which method of numerical maths have been applied without knowing, how they work and which limitations have to be considered. Typical discussions contain statements like:

"I've used ODE45 because everybody uses it. Matlab is slow so I did not care that the integration runs for 2 days. What is stiffness and scaling? What is sensitivity to variation of initial values and parameters?"

A polynomial fit is instable without scaling and the residuals must be controlled and considered in the following computations. Warnings are very important. For modern polyfit implementations they appear a little bit later, but the problem is the same.

[Luigi De Maria]

I am totally aware of what you are talking about. I had to simplify what I wrote in the question and comments because the code behind it is much bigger and I wanted to focus the attention on polyfit.

However, I implemented the scaling months ago because at the beginning the warnings were also on my Matlab 2021a and I spent days trying different solutions. The current scaling works fine.

How can I say it works fine?

Because I am reproducing a test case to validate the code (which will be applied on a new case). Moreover, this code is inside a Monte-Carlo simulation. I don't have just the results of this MC but I also have a Subset Simulation, a Gaussian Mixtures Simulation and another method (called Chan).

So I decided to just turn off warning not because "the results are okay so I do not need scaling", but because:

• I already checked the code (and hence polynomial fitting part) using other similar method (which do not involve polyfit)
• I am analyzing a test case which is certified and the results have always been identical to it
• I had already introduced a proper scaling months ago
• Residuals have always been very low (about 1e-12)

I thank you for all your advices and care, but don't worry if I am not 99.999% sure the code is mathematically correct and computationally fine, I never trust it.

[Jan]

Fine. The correct application of numerical maths is not trivial.

Answer 2

[Bruno Luong]

In R2022a the Vandermonde matrix V is built the condition number is aroun 1e12, and the resolution is carrired out by

[p, rankV, QRfactor, perm] = matlab.internal.math.leastSquaresFit(V,y1);

I don't have the version handly but if someone have, can you check what is the resolution engine used by R2019a by editing polyfit.m, (line #72 in R2022a)

[Bruno Luong]

In the release R2018b the resolution is carried out by QR

[Q,R] = qr(V,0);
...
    p = R\(Q'*y);               % Same as p = V\y
...

I guess it's the same code in R2019b. So here we go the difference. matlab.internal.math.leastSquaresFit(V,y1) won't throw a waring and R\(Q'*y) does.

[Luigi De Maria]

So the issue is related to the computation of the least square problem?

In the 2019b the method used is the one you wrote in the comment here above, whereas on the 2021a is the one wrote in your answer.

Indeed, I know it is not a good practice but I had a very long computation and very few time remained, I turned off the warnings and actually the result is correct. I tried with different test cases and everything was fine even if Matlab gave me warnings about polynomial conditioning.

[Bruno Luong]

The high condition number is due to the scale of your tspan of the order of 1e-4.

As it stands it's OK, but if you have smaller tspan, increase the polynom order, ... I would rather use normalization polyfit as Jan's suggestion.

[Bruno Luong]

The function matlab.internal.math.leastSquaresFit uses QR decomposition WITH permutation, so somewhat R is slighly better conditionned. It is also solved with specified rank so the warning kicks in a little later.

To me it is a slightly better algorithm, but honestly the non-permutation QR used by polyfit in older version should work OK.

[Luigi De Maria]

Thank you very much to all of you for your kind replies. It has been a very formative discover this difference.