matrix inverse results different between r2023b and 2021b

2024 6월 11
1 답변
답변 채택됨
업데이트 시간: 2024 6월 14
조회 수: 16 (30일)

1 개 추천

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'

이 질문에 답변하려면 로그인하십시오.

 채택된 답변

David Goodmanson
David Goodmanson 2024년 6월 11일
편집: David Goodmanson 2024년 6월 13일

0 개 추천

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.

댓글 수: 3
이전 댓글 1개 표시 이전 댓글 1개 숨기기

Qiang Li
Qiang Li 2024년 6월 12일
편집: Qiang Li 2024년 6월 12일
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
David Goodmanson 2024년 6월 13일
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
Qiang Li 2024년 6월 14일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

도움말 센터 및 File Exchange에서 Gaussian Process Regression에 대해 자세히 알아보기

제품

릴리스

R2023b

태그

질문:

Qiang Li
Qiang Li
2024년 6월 11일

댓글:

Qiang Li
Qiang Li
2024년 6월 14일

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by