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

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!

[Stephen23]

"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]

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.

Answer 1

[Walter Roberson]

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.

[Bruno Luong]

@Walter Roberson When compare symbolic with double array, there are two ways, the first one err2_1 is what you use in your post. But there is a second way that give different results (err2_2). Why do you use the first one, is there a legitimate explanation for correctness of the first one?

format long
A = rand(4, 1) + 1j * rand(4, 1);
B = rand(4, 8) + 1j * rand(4, 8);
C2 = A.'*B;

As = sym(A, 'f');
Bs = sym(B, 'f');
C2s = As.'*Bs;

err2_1 = double(C2s - C2)

err2_1 = 
   1.0e-15 *

  Columns 1 through 4

  0.063541347713746 + 0.110278369769590i  0.080906056511433 + 0.059930338286829i  0.164673801607527 + 0.025396397557960i -0.004536241505032 - 0.107223471321709i

  Columns 5 through 8

 -0.019339478112225 + 0.207422136056787i -0.070507706160179 - 0.084393411859533i  0.074925847251870 + 0.151440939527197i  0.145998754128206 + 0.052051741422148i

err2_2 = double(C2s) - C2

err2_2 = 1×8

1.0e-15 * 0.055511151231258 0 0.222044604925031 0 0 -0.111022302462516 0.055511151231258 0.222044604925031

% difference of two results
all(err2_1 == err2_2)

ans = logical
   0

err2_1 - err2_2

ans = 
   1.0e-15 *

  Columns 1 through 4

  0.008030196482488 + 0.110278369769590i  0.080906056511433 + 0.059930338286829i -0.057370803317504 + 0.025396397557960i -0.004536241505032 - 0.107223471321709i

  Columns 5 through 8

 -0.019339478112225 + 0.207422136056787i  0.040514596302337 - 0.084393411859533i  0.019414696020612 + 0.151440939527197i -0.076045850796826 + 0.052051741422148i

[Walter Roberson]

I used double(symbolic) in order to find the closest floating point number to the indefinite-precision symbolic numbers. double(symbolic) is the "ideal" floating point number -- as if the calculation were somehow being done internally with some hypothetical extended floating point precision that was then truncated to standard floating point precision, with the hypothetical extended precision being "big enough" to be sure of giving bitwise-exact answers for these operations.

[Bruno Luong]

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)

[Walter Roberson]

See also the discussion of guard digits at digits

Answer 2

[Askic V]

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

[Walter Roberson]

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.

Answer 3

[Bruno Luong]

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.