Matrix multiplication accuracy Matlab vs Python

Question

I am translating some Matlab code into Python and I having some problems regarding matrix multiplication accuracy. Assuming we have following data:

A: 6x6 matrix
B: 5x5 matrix
C: 2x2 matrix
D: 5x5 matrix
E: 6x5 matrix

In Matlab, my operation looks as follows:

R1 = A * (-( B*C(1,1) + D*C(2,1) ) * E.').'

Previous set of operations produces a 6x5 matrix (R1 matrix).

In my Python code I have the same matrices and the operation looks as follows:

R2 = np.matmul(A, np.matmul(-(np.multiply(B, C[0,0]) + np.multiply(D, C[1,0])), E.transpose()).transpose())

However, the same operation on same data produces different results in terms of norm-2, namely:

norm-2(R1-R2) = 7.4506e-09

I am not able to understand why the results are different. Does anyone know what the reason could be?

For the sake of clarity, I attach the data and the scripts here. The instructions are:

1. Run the Run_python.py script. This will perform the operations and will generate a python_results.mat file.
2. Run the Run_matlab script. This will perform the operations and will compare both results in terms of its norm-2.

Finally, my Python version is 3.6.3, with numpy (1.14.3) and scipy (1.1.0).

Answer 1

First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :)

Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very BLAS-specific reasons why this could happen:

• Matlab and Python are probably linked to different versions of the BLAS libraries. They are not guaranteed to yield the exact same result up to the last bit. Different optimization could rearrange sums in different ways (using associativity), and give slightly different results. See https://bebop.cs.berkeley.edu/reproblas/ for more details.
• Matlab includes special tricks to detect multiplications of the form A*B', A*B.', A*A' and map them to a single BLAS call. I'm not sure Python does that, too, and if it does it might map them to a slightly different call when something like (A*B.').' is used.
• Both languages may include optimizations to map A*B where A and/or B are symmetric to a specialized BLAS call (DSYMM). Again, this might be implemented in slightly different ways in the two languages, and possibly also depend on the result of previous computations.

Answer 2

Here is R1, as computed in MATLAB:

   1.0e+07 *
  -7.382605957465515  -9.599867106092937  -2.830412177259742  -0.000000000002830  -0.000000000002830
  -1.230434326244253  -1.599977851015490  -0.471735362876624  -0.000000000000472  -0.000000000000472
   3.691302978732758   4.799933553046468   1.415206088629871   0.000000000001415   0.000000000001415
  -5.056592907133381  -6.575268976533256  -1.938643647277773  -0.000000000001939  -0.000000000001939
  -2.842303293624223  -3.695948835845781  -1.089708688245001  -0.000000000001090  -0.000000000001090
  -2.522882542531216  -3.280594585722160  -0.967246188042229  -0.000000000000967  -0.000000000000967

norm(R1) = 1.76e8 (largest singular value)

Remaining singular values of R1, per MATLAB double precision, range from 1.48e-8 down to 1.57e-37. So R1 is extremely ill-conditioned.

Given norm(R1) = 1.76e8, you can't expect norm(R1) to be more accurate than your observed MATLAB to PYTHON norm(R1-R2) = 7e-9 (which already could be double the error in norm(R1) due to triangle inequality). All told, you've done about an order of magnitude better than worst case for double precision, which has machine precision 2.2e-16.