Did mldivide Implementation Change in R2022b?

When I run this code here on 2022b (actually, 2023a now, but I saw the same effect on 2022b) I get:

load('mldivdivideex')
format long e
clsys.A\clsys.B

ans = 5×2

-1.000000000000000e+00 2.047143402298295e-17 -7.194341130875050e-12 7.868810611894586e-13 -1.000000000000001e+00 1.999999999991291e-02 -9.738334588696467e-13 8.679999999981491e+01 -2.428571428833528e-03 3.004408164584185e-04

But when I run on my local installation using 2022a I get

>> clsys.A\clsys.B

ans =

-1.000000000000000e+00 2.168317634010083e-17

-6.812536955189428e-12 1.287988200135117e-13

-1.000000000000000e+00 1.999999999991291e-02

-9.670661834891370e-13 8.679999999981482e+01

-2.428571428833527e-03 3.004408164584182e-04

The second row is obviously different, as are the (4,1) and (1,2) elements. Even the "large" elements can be different in the last digit.

The 2022b release notes do mention "mldivide and pagemldivide Functions: Improved performance with small matrices." But the release notes don't say anything about expecting different answers.

[Rik]

This doesn't actually answer your question, but I mistrust most arrays that span 17 orders of magnitude.

I suppose even performance increase worth mentioning will require a change in algortithm, which may or may not be numerically stable. It would have been nice to see a mention of that in the release notes, but I can understand why they left it where they did.

[Star Strider]

‘This doesn't actually answer your question ...’

@Rik — It actually does, albeit indirectly. See the Answer @Matt J posted.

[Bruno Luong]

I get warining when I read the mat file

>> load('C:\Users\b.luong\Downloads\mldivdivideex.mat')
Warning: Variable 'clsys' originally saved as a ss cannot be instantiated as an object and will be read in as a uint32. 
>> clsys.A
Dot indexing is not supported for variables of this type.

Does it requires a toolbox?

[Matt J]

@Bruno Luong Indeed a toolbox is needed to support the 'ss' object type.

load('mldivdivideex')

class(clsys)

ans = 'ss'

@Paul You might consider posting your example using more fundamental Matlab variable types.

[Bruno Luong]

Yes it requires Control System Toolbox apparently

[Paul]

@Matt J @Bruno Luong

The A and B matrices are in the .mat file attached to this comment.

[Bruno Luong]

It seems that the huge condition number is due to the discrepency of scaling between rows 1/5 compared to rows 2/4

format long
load('mldivideexample.mat')
rn = vecnorm(A,2,2)

rn = 5×1

1.0e+08 * 0.000000010000000 9.800028893975272 0.000061376869422 9.800113494010919 0.000000010000000

[a,b]=bounds(rn); b/a

ans = 
     9.800113494010919e+08

[Paul]

Maybe this an example that motivates your desire to rehabilitate inv.

Assuming the sym/vpa solution is as good as can be done ...

load(websave('mldivideexample.mat','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1328525/mldivideexample.mat'));
format long e

C1 = double(inv(sym(A))*sym(B))

C1 = 5×2

-1.000000000000000e+00 0 0 0 -1.000000000000000e+00 1.999999999991288e-02 -4.041939973822555e-13 8.679999999981486e+01 -2.428571428833525e-03 3.004408164584182e-04

ans = logical
   1

The exact -1's and 0's are appealing, and I believe theoretically expected for this problem.

C2 = inv(A)*B

C2 = 5×2

-1.000000000000000e+00 0 0 0 -1.000000000000000e+00 1.999999999991291e-02 -9.738334588696467e-13 8.679999999981491e+01 -2.428571428833528e-03 3.004408164584187e-04

The exact -1's and 0's are appealing

C3 = A\B

C3 = 5×2

-1.000000000000000e+00 2.047143402298295e-17 -7.194341130875050e-12 7.868810611894586e-13 -1.000000000000001e+00 1.999999999991291e-02 -9.738334588696467e-13 8.679999999981491e+01 -2.428571428833528e-03 3.004408164584185e-04

norm(C1 - C2)

ans = 
     5.724687657581802e-13

norm(C1 - C3)

ans = 
     7.259849696053097e-12

[Bruno Luong]

You have good memory @Paul

Answer 1

[Matt J]

The condition number of clsys.A is very poor, so you should not expect any consistency in the result among different implementations of the mldivide operation nor among different CPUs.

load('mldivdivideex')

size(clsys.A)

ans = 1×2

5 5

cond(clsys.A)

ans = 1.9600e+09

[Paul]

Is 2e9 really that big in the context of this particular observation? I truly have no feel for such things.

I guess I can understand how running on different CPUs might change results, but is whether or not such differences are observed (not the amount of difference, just is there any difference at all) be dependent on condition number?

The heart of my concern is whether or not the change in 2022B to make mldvide more efficient with small matrices is also a fundamental change in the algorithm, i.e., different results would obtain on 2022a and 2022b on the same hardware. Unfortunately, I'm not set up to do that test.

[Matt J]

Is 2e9 really that big in the context of this particular observation?

2e9 is really big in any context. It means that if you perturb the right hand side B by 1%, you can expect a change in the solution of up to 2e9%. Basically, your matrix is no different from a singular matrix in the world of finite precision math.

The heart of my concern is whether or not the change in 2022B to make mldvide more efficient with small matrices is also a fundamental change in the algorithm, i.e., different results would obtain on 2022a and 2022b on the same hardware.

Yes, any change in the computational steps of an algorithm is likely to result in some change on the same hardware, due to differences in floating point error accumulation. The question is, how large do we expect the change to be? With a matrix of condition number of 1, you would expect very little change, but with 2e9, very big changes make quite a bit of sense.

[Bruno Luong]

"Is 2e9 really that big in the context of this particular observation? I truly have no feel for such things."

A reasonable condition number one can hope reliable inversion is about below 1e5 to my explerience for DOUBLE matrix. Beyond that sometime it works sometime it is not, it's lottery game.

Not only it's depends on the A, but it is also affected by b (whareas b is projected on the smallest singular vector of A)

The heart of my concern is whether or not the change in 2022B to make mldvide more efficient with small matrices is also a fundamental change in the algorithm, i.e., different results would obtain on 2022a and 2022b on the same hardware.

Fix your condition number if you can (sometime it is just a question of scaling observation and/or decision variables) or add regularization penalty. Don't blame any algorithm wuth such Cond number, even for the same hardware but different multi-thread setup can result different result (for instant if someone run parallel MATLAB with some other task on computer, the developer can decide how many CPU cores MATLAB can use as upperbound). If you want to be strick about the stability of the result for such fragile numerical code, you should not change MATLAB version, OS, Hardware, then you always get the same result. But still you have no idea about the correctness of the result.

[Steven Lord]

Let's take a sample matrix with a condition number of 1e9.

format longg
A = [1 0; 0 1e9]

A = 2×2

1.0e+00 * 1 0 0 1000000000

cond(A)

ans = 
                1000000000

Now let's multiply A by two vectors that are very close to one another.

x1 = [1; 1]

x1 = 2×1

1 1

y1 = A*x1

y1 = 2×1

1.0e+00 * 1 1000000000

x2 = [1; 1.1]

x2 = 2×1

1 1.1

y2 = A*x2

y2 = 2×1

1.0e+00 * 1 1100000000

How different are y1 and y2? The first elements are the same. The second elements:

difference = y2(2)-y1(2)

difference = 
   100000000

The differences in the second element of the x vectors was only 0.1. The differences in the second elements of the y vectors was much, much larger. Granted this was a contrived example, but hopefully you can see that badly conditioned matrices can magnify differences or errors quite substantially.

From the Wikipedia entry on condition number: "As a rule of thumb, if the condition number , then you may lose up to k digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods."

[Bruno Luong]

"As a rule of thumb, if the condition number is 10^k, then you may lose up to k digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods."

Hmmm I admit I don't understand this sentence from Wikipedia. Or at least I'm not sure how to interpret it correctly.