matrix inverse results different between r2023b and 2021b

I was solving a high-dimensional nonlinear optimization problem with gradient descent method and found out that R2023b and R2021b calculates matrix inverse differently, resulting in completely different gradients and convergence character of the algorithm.

2020b and 2021b gave the following gradient and converged after 18 iterations.

direction = [4.33783e-06 2.33596e-06 -5.33992e-06 -8.36479e-06 -9.19797e-05 1 ]

2023b gave the following gradient and didn't converge after 44 iterations.

direction = [4.17191e-06 2.24661e-06 -5.13617e-06 -7.53861e-06 -5.19154e-05 1 ]

The code is very slow, so I can't let it run forever and find out when it will converge.

I narrowed it down to the inverse of a 10by10 matrix

A_2stepblock =

   1.0e+06 *

  Columns 1 through 7

                   0                   0                   0                   0                   0                   0                   0
                   0                   0                   0                   0                   0                   0                   0
                   0                   0                   0                   0                   0                   0                   0
                   0                   0                   0                   0                   0                   0   0.000006000000000
                   0                   0                   0                   0                   0   0.000024000000000                   0
   1.953050781823359   0.390611805629631   0.078122690980311   0.015624604167218   0.003124934027730   0.000624989444415   0.000124998416659
   3.515506250666680   0.624981527842487   0.109372229170523   0.018749604166382   0.003124947222073   0.000499993666635   0.000074999366662
   5.624833750582380   0.874977833364188   0.131247229164674   0.018749683332437   0.002499968333177   0.000299997466649   0.000029999873332
   7.874800500277688   1.049977833317391   0.131247783327059   0.014999809999061   0.001499987333244   0.000119999493329   0.000006000000000
   9.449800499856519   1.049982266616475   0.104998669993428   0.008999923999464   0.000599997466643   0.000024000000000                   0

  Columns 8 through 10

                   0                   0   0.000001000000000
                   0   0.000001000000000                   0
   0.000002000000000                   0                   0
                   0                   0                   0
                   0                   0                   0
   0.000024999788887   0.000004999978889   0.000001000000000
   0.000009999957777   0.000001000000000                   0
   0.000002000000000                   0                   0
                   0                   0                   0
                   0                   0                   0

r2021b gives

inv(A_2stepblock)

ans =

  Columns 1 through 7

  -0.000035841361961  -0.000089603026574  -0.000096002837407  -0.000053334684477  -0.000013333614821   0.000035841361961  -0.000089603026574
   0.000806427239167   0.002048060531355   0.002240056748025   0.001280027022812   0.000333338963074  -0.000806427239167   0.001984058639750
  -0.006912204293321  -0.017920453984200  -0.020160425609289  -0.012000202670662  -0.003333375556299   0.006912204293321  -0.016640421556757
   0.026880680976300   0.071681513277471   0.084001418694634   0.053334008900781   0.016666807409588  -0.026880680976300   0.062721324117785
  -0.040320851218577  -0.112001891592844  -0.140001773364547  -0.100000844457525  -0.041666842594947   0.040320851218577  -0.089601513274271
                   0                   0                   0                   0   0.041666666666667                   0                   0
                   0                   0                   0   0.166666666666667                   0                   0                   0
                   0                   0   0.500000000000000                   0                   0                   0                   0
                   0   1.000000000000000                   0                   0                   0                   0                   0
   1.000000000000000                   0                   0                   0                   0                   0                   0

  Columns 8 through 10

   0.000096002837407  -0.000053334684477   0.000013333614821
  -0.002080052694595   0.001120023644960  -0.000266671170459
   0.016960358052258  -0.008800148625151   0.002000025333779
  -0.061601040376060   0.030667055117945  -0.006666722963834
   0.084001064018720  -0.040000337783004   0.008333368518987
                   0                   0                   0
                   0                   0                   0
                   0                   0                   0
                   0                   0                   0
                   0                   0                   0

and r2023b gives

inv(A_2stepblock)

ans =

  Columns 1 through 6

  -0.000035841361961  -0.000089603026574  -0.000096002837407  -0.000053334684477  -0.000013333614821   0.000035841361961
   0.000806427239167   0.002048060531355   0.002240056748026   0.001280027022812   0.000333338963074  -0.000806427239167
  -0.006912204293323  -0.017920453984205  -0.020160425609294  -0.012000202670665  -0.003333375556300   0.006912204293323
   0.026880680976306   0.071681513277487   0.084001418694652   0.053334008900791   0.016666807409591  -0.026880680976306
  -0.040320851218585  -0.112001891592865  -0.140001773364570  -0.100000844457540  -0.041666842594951   0.040320851218585
                   0                   0                   0                   0   0.041666666666667                   0
                   0                   0                   0   0.166666666666667                   0                   0
                   0                   0   0.500000000000000                   0                   0                   0
                   0   1.000000000000000                   0                   0                   0                   0
   1.000000000000000                   0                   0                   0                   0                   0

  Columns 7 through 10

  -0.000089603026574   0.000096002837407  -0.000053334684477   0.000013333614821
   0.001984058639750  -0.002080052694595   0.001120023644961  -0.000266671170459
  -0.016640421556761   0.016960358052262  -0.008800148625154   0.002000025333779
   0.062721324117800  -0.061601040376075   0.030667055117952  -0.006666722963836
  -0.089601513274291   0.084001064018740  -0.040000337783014   0.008333368518990
                   0                   0                   0                   0
                   0                   0                   0                   0
                   0                   0                   0                   0
                   0                   0                   0                   0
                   0                   0                   0                   0

I used

num2str(A,100)

to make sure it's the same matrix to the last digit in both versions.

From the release notes, they did changed quite a few matrix operations from 2022a.

So, which version is more accurate?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Edit after David Goodmanson's response, the second "r2021b" should be "r2023b", I corrected that.

Running DG's answer directly on my computer to show the difference between the versions, on the r2021b

K>> A = A_2stepblock;
K>> F = inv(A);
K>> max(abs(A*F-eye(10,10)),[],'all')
max(abs(F*A-eye(10,10)),[],'all')

ans =

     1.818989403545856e-12


ans =

     3.887682601841180e-11

K>> norm(A*F-eye(10,10))
norm(F*A-eye(10,10))

ans =

     5.914730386430228e-12


ans =

     4.379288551440263e-11

K>> num2hex(A)

ans =

  100×16 char array

    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '413dcd1ac825935f'
    '414ad2392015d885'
    '4155750070098aae'
    '415e0a3c20048cb5'
    '416206290ffed319'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4117d74f38f6f95a'
    '412312ab0e415ed4'
    '412ab3c3aaaeb5f6'
    '41300579d55449dd'
    '4130057e4440fa32'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '40f312ab0e415ed4'
    '40fab3c3aaaeb5f6'
    '41000579d55449dd'
    '4100057e4440fa32'
    '40f9a26ab84b0751'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '40ce844d5559f487'
    '40d24f66aaa978fc'
    '40d24f6bbbb7f95e'
    '40cd4be7ae0c9aa6'
    '40c193f6459d4bb4'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '40a869de38e1f6a6'
    '40a869e4fa4aa1d2'
    '40a387efc95dbc6f'
    '40976ff3077c64f0'
    '4082bffacfcaa3cf'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4038000000000000'
    '408387ea61d54e42'
    '407f3fe60efc60b1'
    '4072bff59f96b726'
    '405dfff7b2ddd2e6'
    '4038000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4018000000000000'
    '0000000000000000'
    '405f3fe60efc60b1'
    '4052bff59f96b726'
    '403dfff7b2ddd2e6'
    '4018000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4038fff22a1e4988'
    '4023fffa773e8c99'
    '4000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4013fffa773e8c99'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'

On the r2023b

K>> F = inv(A_2stepblock);
K>> A = A_2stepblock;
K>> max(abs(A*F-eye(10,10)),[],'all')
max(abs(F*A-eye(10,10)),[],'all')

ans =

     2.728484105318785e-12


ans =

     1.472269235603454e-11

K>> norm(A*F-eye(10,10))
norm(F*A-eye(10,10))

ans =

     7.261754408131635e-12


ans =

     2.136608381476120e-11
     
num2hex(A_2stepblock)

ans =

  100×16 char array

    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '413dcd1ac825935f'
    '414ad2392015d885'
    '4155750070098aae'
    '415e0a3c20048cb5'
    '416206290ffed319'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4117d74f38f6f95a'
    '412312ab0e415ed4'
    '412ab3c3aaaeb5f6'
    '41300579d55449dd'
    '4130057e4440fa32'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '40f312ab0e415ed4'
    '40fab3c3aaaeb5f6'
    '41000579d55449dd'
    '4100057e4440fa32'
    '40f9a26ab84b0751'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '40ce844d5559f487'
    '40d24f66aaa978fc'
    '40d24f6bbbb7f95e'
    '40cd4be7ae0c9aa6'
    '40c193f6459d4bb4'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '40a869de38e1f6a6'
    '40a869e4fa4aa1d2'
    '40a387efc95dbc6f'
    '40976ff3077c64f0'
    '4082bffacfcaa3cf'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4038000000000000'
    '408387ea61d54e42'
    '407f3fe60efc60b1'
    '4072bff59f96b726'
    '405dfff7b2ddd2e6'
    '4038000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4018000000000000'
    '0000000000000000'
    '405f3fe60efc60b1'
    '4052bff59f96b726'
    '403dfff7b2ddd2e6'
    '4018000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4038fff22a1e4988'
    '4023fffa773e8c99'
    '4000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '4013fffa773e8c99'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '3ff0000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'
    '0000000000000000'

Answer 1

[David Goodmanson]

Hi QL.

you called both the inverse matrices '2021b' so I called the first one F and the second one S (and the orginal matrix is A). So why not find out. Sounds like you were careful to get all the digits in place but there will still be a bit of uncertainty in copying numbers off of the display as opposed to having them in memory. Num2hex is the real solution because it shows exactly what is stored in memory**. But with that caveat and on the basis of how large [A*inverse(A) - I] (or the other way round) is, F is slightly better in this particular case. In both measures below, S*A is the outlier.

cond(A)
ans =
   1.4390e+07    % pretty close to singular

format compact

% largest nonzero element of the matrix that should be zero
max(abs(A*F-eye(10,10)),[],'all')
max(abs(F*A-eye(10,10)),[],'all')

max(abs(A*S-eye(10,10)),[],'all')
max(abs(S*A-eye(10,10)),[],'all')

ans =
2.7399e-09
ans =
5.1108e-09
ans =
2.9945e-09
ans =
9.9112e-09

% norm of the matrix that should be zero
norm(A*F-eye(10,10))
norm(F*A-eye(10,10))

norm(A*S-eye(10,10))
norm(S*A-eye(10,10))

ans =
9.2492e-09
ans =
6.7615e-09
ans =
8.6579e-09
ans =
1.3678e-08

** for example

num2hex(1/100)    ans =  '3f847ae147ae147b'
num2hex(1/128)    ans =  '3f80000000000000'

illustrates the difference between storing inverse powers of 10 and inverse powers of 2.

[Qiang Li]

Hi @David Goodmanson,

I've updated the post according to your answer. Both measures seem to show that both versions have smaller errors when multiplying the inverse on the left. Does that indicate anything? Thx.

[David Goodmanson]

HI Ql, It doesn't seem to. If you run the following example a few times there are cases where b*a is bigger in both measures, smaller in both measures, or one of each.

a = 2*rand(10,10) -1;
b = inv(a);
maxx(abs(a*b-eye(size(a))))
maxx(abs(b*a-eye(size(a))))
norm(a*b-eye(size(a)))
norm(b*a-eye(size(a)))

function b = maxx(a)
% maximum value of the elements of an n-dimensional array
%
% b = maxx(a)
b = max(a(:));
end

[Qiang Li]

I've improved my code over the past couple of days. The difference is still there, but the code is more robust now and I don't feel like that it matters that much any more. Thanks for all the help.