The multiplication result by sum() and matrix multiplication are not the same

2 次查看（过去 30 天）
Wei Zhao
Wei Zhao 2023-7-30
编辑： Bruno Luong 2023-7-31
Ran in:
Assume 2 random complex matrices and we want to multiply them.
A = rand(4, 1) + 1j * rand(4, 1)
A =
0.9628 + 0.0224i 0.2966 + 0.3287i 0.1870 + 0.6470i 0.8037 + 0.7556i
B = rand(4, 8) + 1j * rand(4, 8)
B =
0.1255 + 0.3132i 0.8355 + 0.9817i 0.5459 + 0.3694i 0.1726 + 0.9117i 0.3396 + 0.6657i 0.8552 + 0.6826i 0.0278 + 0.3613i 0.0696 + 0.4578i 0.0828 + 0.4257i 0.6748 + 0.4067i 0.7103 + 0.6342i 0.5629 + 0.9616i 0.8074 + 0.4321i 0.5599 + 0.0596i 0.2606 + 0.6650i 0.9858 + 0.5410i 0.5360 + 0.6049i 0.7389 + 0.9686i 0.3796 + 0.8984i 0.0507 + 0.5859i 0.1426 + 0.9283i 0.4000 + 0.1699i 0.2336 + 0.1639i 0.5021 + 0.1212i 0.8022 + 0.6056i 0.6300 + 0.1437i 0.4667 + 0.6784i 0.4866 + 0.4009i 0.7399 + 0.3622i 0.1861 + 0.0170i 0.5266 + 0.7786i 0.8595 + 0.6314i
C1 = sum(A.*B,1)
C1 =
-0.1056 + 2.0105i 0.7581 + 2.5571i -0.1283 + 2.1011i -0.2848 + 2.1842i 0.1565 + 2.1580i 1.0562 + 1.3229i -0.3501 + 1.8368i 0.4004 + 2.4312i
C2 = A.'*B
C2 =
-0.1056 + 2.0105i 0.7581 + 2.5571i -0.1283 + 2.1011i -0.2848 + 2.1842i 0.1565 + 2.1580i 1.0562 + 1.3229i -0.3501 + 1.8368i 0.4004 + 2.4312i
sum(C1 == C2, "all")
ans = 3
C1 - C2
ans =
1.0e-15 * 0.0139 - 0.4441i 0.0000 - 0.4441i 0.0555 - 0.4441i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.2220 + 0.0000i 0.0000 + 0.0000i 0.0555 + 0.0000i
The weird thing happens. The result is not the same.
May I ask the reason?
Which one is right?
How to prevent this from happening?
Thanks a lot!
  2 个评论
Stephen23
Stephen23 2023-7-30
编辑：Stephen23 2023-7-30
"May I ask the reason?"
Different numeric operations performed in a different order. No surprises there.
"Which one is right?"
Neither one is more "right" than the other.
"How to prevent this from happening?"
You cannot prevent binary floating point numbers from behaving like binary floating point numbers.
You can learn that binary floating point numbers have a limited precision and that operations on them accumulate floating point error. This means that all computations that use binary floating point numbers must take that limited precision and floating point error into account, e.g. not assuming exact equivalence between two different algorithms:
https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero
https://www.mathworks.com/matlabcentral/answers/316889-why-two-equal-numbers-are-not-equal
http://www.mathworks.com/matlabcentral/answers/102419-how-do-i-determine-if-the-error-in-my-answer-is-the-result-of-round-off-error-or-a-bug
http://www.mathworks.com/matlabcentral/answers/140656-about-numerical-difficulty-of-equality-comparison-operation-between-two-double-real-numbers-in-matlab
This is worth reading as well:
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
Bruno Luong
Bruno Luong 2023-7-30
编辑：Bruno Luong 2023-7-30
Ran in:
Stephen23 reply should be in answers
To illustrate the order of operation mater, a simple example of suming an array of 100 elements in two different orders
for ntest = 1:10
a = rand(1,100);
s1 = sum(a);
s2 = sum(a(randperm(end)));
fprintf('s1 = %1.16f, s2 = %1.16f\n', s1, s2);
fprintf('\ts1-s2 = %1.16f\n', s1-s2);
end
s1 = 50.5746657174914631, s2 = 50.5746657174914560
s1-s2 = 0.0000000000000071
s1 = 50.0609477280682142, s2 = 50.0609477280682214
s1-s2 = -0.0000000000000071
s1 = 49.7007348624369030, s2 = 49.7007348624369030
s1-s2 = 0.0000000000000000
s1 = 50.3659969116610000, s2 = 50.3659969116610142
s1-s2 = -0.0000000000000142
s1 = 44.3158313068657037, s2 = 44.3158313068657037
s1-s2 = 0.0000000000000000
s1 = 50.0573409775813971, s2 = 50.0573409775813971
s1-s2 = 0.0000000000000000
s1 = 46.3840663641301845, s2 = 46.3840663641301774
s1-s2 = 0.0000000000000071
s1 = 48.4878190738828110, s2 = 48.4878190738828181
s1-s2 = -0.0000000000000071
s1 = 44.9918199491275743, s2 = 44.9918199491275672
s1-s2 = 0.0000000000000071
s1 = 50.6535878109839928, s2 = 50.6535878109839928
s1-s2 = 0.0000000000000000
You just have to accept this when working with finite precision floating point numbers and make the code robust to round-off error.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Walter Roberson
Walter Roberson 2023-7-30
Ran in:
format long g
A = rand(4, 1) + 1j * rand(4, 1)
A =
0.436623964048917 + 0.822164039093503i 0.488439227712175 + 0.633998025984701i 0.888980048931903 + 0.711911618725595i 0.360295102265367 + 0.502544697727889i
B = rand(4, 8) + 1j * rand(4, 8)
B =
Columns 1 through 3 0.0566523672499591 + 0.127739189477948i 0.440193960894806 + 0.358676840152069i 0.831007582073232 + 0.154686156693833i 0.130751893609349 + 0.134919277784539i 0.139794683225943 + 0.702720459382693i 0.621594600648509 + 0.889235380682137i 0.221647466555571 + 0.601467418925383i 0.101566918535911 + 0.807592614267102i 0.717370309162311 + 0.789958018186872i 0.734238258245452 + 0.944387635313803i 0.110040345713131 + 0.811944894656452i 0.239720351892684 + 0.689231189863182i Columns 4 through 6 0.914975858281857 + 0.0158215359399975i 0.986931904120533 + 0.263658695854385i 0.529912834899685 + 0.578143306898728i 0.229160081857441 + 0.501641924562142i 0.585123867371859 + 0.433747384303353i 0.475105105028833 + 0.170953714435258i 0.0764326576749834 + 0.79996325316813i 0.741536067276155 + 0.642763015085329i 0.348002500998525 + 0.215646024698359i 0.235830872979799 + 0.242155001648183i 0.765612881066247 + 0.0880723247953943i 0.944129213078834 + 0.913087607163485i Columns 7 through 8 0.00421695171173275 + 0.859831645625063i 0.245487025367385 + 0.364608214917568i 0.855464396612252 + 0.129877132349552i 0.170715658769038 + 0.131633454263892i 0.758988020151318 + 0.898411718210346i 0.0702563535057346 + 0.149101170604057i 0.442072389008343 + 0.669075060662297i 0.846749956236221 + 0.802113071649759i
C1 = sum(A.*B,1)
C1 =
Columns 1 through 3 -0.543167007195224 + 1.65287956604633i -1.33296934955559 + 2.08846469258086i -0.20915367022258 + 3.16094941498314i Columns 4 through 6 -0.357897758492674 + 2.12080442766428i 0.658156488650943 + 3.02516459958456i -0.0831361750164952 + 2.31572285832962i Columns 7 through 8 -0.511409053660225 + 2.78691803830161i -0.334361538355948 + 1.43064701412489i
C2 = A.'*B
C2 =
Columns 1 through 3 -0.543167007195224 + 1.65287956604633i -1.33296934955559 + 2.08846469258086i -0.20915367022258 + 3.16094941498314i Columns 4 through 6 -0.357897758492674 + 2.12080442766428i 0.658156488650943 + 3.02516459958456i -0.0831361750164952 + 2.31572285832962i Columns 7 through 8 -0.511409053660225 + 2.78691803830161i -0.334361538355948 + 1.43064701412489i
sum(C1(:) == C2(:))
ans =
2
C1 - C2
ans =
Columns 1 through 3 0 - 2.22044604925031e-16i 0 + 0i 0 - 4.44089209850063e-16i Columns 4 through 6 -5.55111512312578e-17 + 0i 1.11022302462516e-16 + 4.44089209850063e-16i 0 + 0i Columns 7 through 8 0 + 4.44089209850063e-16i 5.55111512312578e-17 + 0i
As = sym(A, 'f')
As = 
Bs = sym(B, 'f')
Bs = 
C1s = sum(As.*Bs,1)
C1s = 
C2s = As.'*Bs
C2s = 
sum(C1s(:) == C2s(:))
ans = 
C1s - C2s
ans = 
C1s - C1
ans = 
double(ans)
ans =
Columns 1 through 3 1.06798582793158e-16 + 8.68260406164711e-17i 9.55856578796629e-17 + 1.07481835083934e-16i 1.77518227272422e-17 + 2.57767199975482e-16i Columns 4 through 6 -4.11627663808118e-17 + 1.11900343514023e-17i -1.09064614719157e-16 - 1.2559753058777e-16i -4.57013833099403e-18 + 1.27226375380194e-16i Columns 7 through 8 -1.51210486886328e-17 - 2.46976885627192e-16i -2.15644517837071e-17 + 5.3277460502552e-17i
C2s - C2
ans = 
double(ans)
ans =
Columns 1 through 3 1.06798582793158e-16 - 1.3521856430856e-16i 9.55856578796629e-17 + 1.07481835083934e-16i 1.77518227272422e-17 - 1.8632200987458e-16i Columns 4 through 6 -9.66739176120697e-17 + 1.11900343514023e-17i 1.9576877433582e-18 + 3.18491679262293e-16i -4.57013833099403e-18 + 1.27226375380194e-16i Columns 7 through 8 -1.51210486886328e-17 + 1.97112324222871e-16i 3.39466994475507e-17 + 5.3277460502552e-17i
This tells you that if you take the numbers and convert them to the exact equivalent symbolic fractions, and do the calculations both ways with the symbolic fractions, that the two ways come out the same. And then from the fact that neither C1 nor C2 is exactly equal to the symbolic calculation, it follows that neither of the C1 or C2 calculation is "right"
double(C1s) - C1
ans =
Columns 1 through 3 1.11022302462516e-16 + 0i 0 + 0i 2.77555756156289e-17 + 4.44089209850063e-16i Columns 4 through 6 -5.55111512312578e-17 + 0i -1.11022302462516e-16 + 0i 0 + 0i Columns 7 through 8 0 - 4.44089209850063e-16i 0 + 0i
norm(ans)
ans =
6.50333926369705e-16
double(C2s) - C2
ans =
Columns 1 through 3 1.11022302462516e-16 - 2.22044604925031e-16i 0 + 0i 2.77555756156289e-17 + 0i Columns 4 through 6 -1.11022302462516e-16 + 0i 0 + 4.44089209850063e-16i 0 + 0i Columns 7 through 8 0 + 0i 5.55111512312578e-17 + 0i
norm(ans)
ans =
5.24426158823621e-16
In my test runs, the norms of the differences are not consistently better for one version compared to the other, so you cannot really come out with a "best" solution for the two.
Why does it happen? It happens because MATLAB calls high-performance library routines, and there happens to be a specific library routine for A' * B that does the calculations internally without ever explicitly calculating A' -- but the internal order of the calculations is not necessarily going to be exactly the same compared to the routine that calculates .* -- potentially it might even involve using different hardware instructions.
  4 个评论
显示 2更早的评论
Bruno Luong
Bruno Luong 2023-7-30
编辑：Bruno Luong 2023-7-30
How double(sym) works? For examle in this case:
x = sym(2);
ys = exp(x)
y =double(ys)
I imagine symbolic represents somehow infinite series
ys = 1 + 2/2! + ... 2^n/n! + ...
When I need y =double(ys)
how MATLAB do the conversion, obviously it cannot compute inifinite sum, so it must truncate. Does the calculation is initite precision on truncation, then find the closest FP?
Does it compute the closest FP on each terms then sum the truncated series?
Do something else?
EDIT: according to the doc https://www.mathworks.com/help/symbolic/double.html double command do internal calculation by evaluate symbolic expression in 32 digits by default. So about twice as numerical double. But it can suffer from cancellation and round-off error as well (see examole High-Precision Conversion)

请先登录，再进行评论。

Askic V
Askic V 2023-7-30
In addition what Walter and Stephen said, if you have some specific application and you need to compare the floating point numbers to check whether they're equal,, I would suggest you to look for "isequal with tolerance":
https://www.mathworks.com/matlabcentral/fileexchange/113375-isequal-function-with-tolerance
  1 个评论
Walter Roberson
Walter Roberson 2023-7-30
ismembertol can also be used for "isequal with tolerance", if you can restrict one of the parameters to being a scalar.
It can handle (scalar, array) [is this number approximately equal to any of these numbers?] or (array, scalar) [are any of these numbers approximately equal to this number), but by itself is not suitable for "is each array element approximately equal to the corresponding array element", which is where the File Exchange contributions are useful.

请先登录，再进行评论。

Bruno Luong
Bruno Luong 2023-7-31
编辑：Bruno Luong 2023-7-31
Ran in:
You cannot prevent the roundoff error to occur and operation order dependency.
Even for the same computer the result can change by setting different number of computational threads as showed here on MATLAB online server
format long
a=randn(1,1e6);
s1=sum(a);
default_nthread = maxNumCompThreads(1)
default_nthread =
2
s2=sum(a);
maxNumCompThreads(default_nthread)
ans =
1
d = s1-s2;
fprintf('s1 = %1.16f, s2 = %1.16f\n', s1, s2);
s1 = 790.9046894595330741, s2 = 790.9046894595336425
fprintf('\ts1-s2 = %1.16f\n', s1-s2);
s1-s2 = -0.0000000000005684
and sum order
s1 = sum(a);
s3 = sum(a(randperm(end)));
fprintf('s1 = %1.16f, s3 = %1.16f\n', s1, s3);
s1 = 790.9046894595330741, s3 = 790.9046894595322783
fprintf('\ts1-s3 = %1.16f\n', s1-s3);
s1-s3 = 0.0000000000007958
Multiplication of matrices have infinity of methods, implementations with different ways to accumulate the result https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm so it is expected that two different functions will not lead to the same numerical result.

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by