Matlab produces different matrix product results depending on computation form

Question

In Matlab, I compute a rotation of a collection of 2D points in two ways: one by a regular matrix-matrix product, and the other by iterative vector-matrix products, as follows.

>> points = % read data from some file, Nx2 matrix
>> R = [cosd(-18), sind(-18); -sind(-18), cosd(-18)];  % My rotation matrix
>> prod1 = points * R;
>> numpt = size(points, 1);
>> for k=1:numpt, prod2(k,:) = points(k,:) * R; end;

I am using a "regular" (intel-based) PC with Windows 10 OS.

It turns out that on some computers, prod1 ~= prod2, and other computers, prod1 == prod2. This can be checked by

>> max(max(abs(prod2 - prod1)))

ans =

   1.1102e-16

This difference is equal to 0 on the "weaker" computers and nonzero on my "powerful" computer.

I suppose that the reason for this happening on some computers but not others is that where it happens, there is some H/W acceleration of matrix multiplication (maybe involving madd ternary operations, notorious for this kind of difference).

Is this some known issue, like a "bug"? Is there a workaround, for example to disable or suspend this sort H/W acceleration?

I am seeking to obtain identical outcomes of the computation on different computers, as part of unit test. I can settle for "near equality". But I should not if I can get true equality.

EDIT 1

I stress that the core issue is that the exact same syntactical expression produces different results on different computers, and the apparent cause is different computational optimizations done on different computers. Bit identity is a requirement that cannot be waved off. I would like both platforms, which are 64-bit intel-based Windows 10, to compute exactly the same outcome for exactly the same input and expression.

EDIT 2

I have tried to extract specific elements of the input on which the outcome differs (i.e. prod2(k,:) ~= prod1(k,:). I collected them into a matrix to repeat the computation just on them. And wow! This time the outcome is not the same as in prod1. It's more similar to the explicitly iterative case.

So I'm sorry, I can't quite reproduce the problem in the scope of this discussion. But we can discuss some more.

I'll check on Zhang's thread issue. And I will even try to repeat the computation several times, see if the differences appear at the same location. If it's about multithreading, I expect the differences to show up in non-deterministic indexes.

Answer 1

Each variable in each step need to be verified to answer all the questions. In general, the result should be the same for well-defined math functions across platforms. Meanwhile, you can leverage on the symbolic toolbox to increase precision of numeric calculations.

Answer 2

The official response from MathWorks tech support points to these links:

1. An answer to a similar question explaining that the details of computation order are in the LAPACK library, and, therefore, not tunable from Matlab.
2. An example from the Symbolic Toolbox that shows how to improve overall precision, with the price being the Symbolic Toolbox itself, a new code that uses it, and computation time increase.

Given this, my preference is to round the computation outcome. As explained in the conversations, my goal is repeatability, not precision. And it looks like the most cost effective answer here is rounding.

Answer 3

Floating-point operations incur a rounding error. The order of operations, which changes depending on platform, compiler, compiler options, etc., affects the rounding error. So, internal MATLAB operations will produce slightly different results from one version of MATLAB to another, and from one platform to another, and between different ways of expressing the same value. You simply cannot expect two floating-point arrays to be identical.

So, instead of comparing two arrays using the == operation, which requires them to be identical, you want to compare them with a tolerance. One way to do so is:

tol = 10
all( abs(prod1 - prod2) < tol * eps(max(abs(prod1),[],'all')) ,'all')
%    ^^^^^^^^^^^^^^^^^^   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%    absolute difference              tolerance

Here we are ensuring that the absolute difference between the two arrays is less than some tolerance, for all array elements.

The tolerance is constructed by taking the eps for the maximum absolute value of the array, and scaling that with a parameter tol that you set yourself. eps(x) is the smallest value you can add to x such that it is no longer equal to x. That is, it is the precision of the value x. A rounding error is always within eps, but rounding errors can accumulate, so you might want to start with tol = 10 as in the example above.

Since you are talking about unit tests, you might be using MATLAB's Testing Framework, in which case you'd use verifyEqual.