Matrix multiply result different from loops

I know I'm having numeric precision issues with a matrix mutiplication operation, but I can't reproduce the same result using for loops. Using the two matricies found in the attached mat file, run the following script. No matter which method I used to calculate the resulting matrix, I can't reproduce the same result as the matrix multiply. How is the Matlab matrix multiply computed?

load matricies
C = A*B; % normal matix multiply A is 54 x 3 and B is 3 x 54, so C is 54 x 54
Ci1 = zeros(54);
Ci2 = zeros(54);
Ci3 = zeros(54);
for i=1:54
  for j=1:54
     for k=1:3
        Ci1(i,j) = Ci1(i,j)+A(i,k)*B(k,j);
     end 
     for k=3:-1:1
        Ci2(i,j) = Ci2(i,j)+A(i,k)*B(k,j);
     end
     Ci3(i,j) = A(i,:)*B(:,j);
   end 
end
diff = C - Ci1;
fprintf('max error using k=1:3: %d\n', max(diff(:)))
diff = C - Ci2;
fprintf('max error with k=3:-1:1 %d\n', max(diff(:)))
diff = C - Ci3;
fprintf('max error with A(i,:)*B(:,j) %d\n', max(diff(:)))
eps_max_orig = eps(max(abs(C(:))));
eps_min_orig = eps(min(abs(C(:))));
fprintf('range of eps is [%d, %d]\n', eps_min_orig, eps_max_orig)

Answer 1

[James Tursa]

MATLAB calls 3rd party BLAS library code to do matrix multiply. This is a highly optimized multi-threaded library. The ordering of the operations is not published, but likely depends on size of the matrix, number of cores used, cache sizes for your CPU, etc. You might get lucky with a guess at operation order for your particular matrix size and your particular machine, but this wouldn't necessarily tell you what would happen with other matrix sizes or on a different machine.

[Davis Pan]

Thanks for your quick answer. I tried every permutation of the ordering for k (script below) and still could not get the same results as the matrix multiply, so I assume the operation ordering can change for calculation of each matrix element?

load matricies

C = A*B; % normal matix multiply A is 54 x 3 and B is 3 x 54, so C is 54 x 54

Ci1 = zeros(54);

Ci2 = zeros(54);

Ci3 = zeros(54);

Ci4 = zeros(54);

Ci5 = zeros(54);

Ci6 = zeros(54);

Ci7 = zeros(54);

for i=1:54

for j=1:54

for k=1:3 % order 1,2,3

Ci1(i,j) = Ci1(i,j)+A(i,k)*B(k,j);

end

for k=3:-1:1 % order 3,2,1

Ci2(i,j) = Ci2(i,j)+A(i,k)*B(k,j);

end

Ci3(i,j) = A(i,:)*B(:,j);

% order 1,3,2

Ci4(i,j) = Ci4(i,j)+A(i,1)*B(1,j);

Ci4(i,j) = Ci4(i,j)+A(i,3)*B(3,j);

Ci4(i,j) = Ci4(i,j)+A(i,2)*B(2,j);

% order 2,1,3

Ci5(i,j) = Ci5(i,j)+A(i,2)*B(2,j);

Ci5(i,j) = Ci5(i,j)+A(i,1)*B(1,j);

Ci5(i,j) = Ci5(i,j)+A(i,3)*B(3,j);

% order 2,3,1

Ci6(i,j) = Ci6(i,j)+A(i,2)*B(2,j);

Ci6(i,j) = Ci6(i,j)+A(i,3)*B(3,j);

Ci6(i,j) = Ci6(i,j)+A(i,1)*B(1,j);

% order 3,1,2

Ci7(i,j) = Ci7(i,j)+A(i,3)*B(3,j);

Ci7(i,j) = Ci7(i,j)+A(i,1)*B(1,j);

Ci7(i,j) = Ci7(i,j)+A(i,2)*B(2,j);

end

end

diff = C - Ci1;

fprintf('max error using k=1,2,3: %d\n', max(diff(:)))

diff = C - Ci2;

fprintf('max error with k=3,2,1 %d\n', max(diff(:)))

diff = C - Ci3;

fprintf('max error with A(i,:)*B(:,j) %d\n', max(diff(:)))

diff = C - Ci4;

fprintf('max error using k=1,3,2: %d\n', max(diff(:)))

diff = C - Ci5;

fprintf('max error using k=2,1,3: %d\n', max(diff(:)))

diff = C - Ci6;

fprintf('max error using k=2,3,1: %d\n', max(diff(:)))

diff = C - Ci7;

fprintf('max error using k=3,1,2: %d\n', max(diff(:)))

eps_max_orig = eps(max(abs(C(:))));

eps_min_orig = eps(min(abs(C(:))));

fprintf('range of eps is [%d, %d]\n', eps_min_orig, eps_max_orig)

[James Tursa]

So, the BLAS library is likely using combined operations in the background using higher precision intermediate results than you are using in your looping. E.g.,

>> format long
>> load matricies
>> C = A*B;
>> A1 = vpa(A(1,:));
>> B1 = vpa(B(:,1));
>> C(1,1)
ans =
      3.307031567995865e-05 + 8.712100175967504e-04i
>> double(double(A1(1)*B1(1) + A1(2)*B1(2)) + A1(3)*B1(3))
ans =
      3.307031567995865e-05 + 8.712100175967504e-04i    % matches
>> A(1,:)*B(:,1)
ans =
      3.307031568011709e-05 + 8.712100175967546e-04i    % does not match

So, by doing intemediate operations in higher precision and selectively converting intermediate results into double, I can match the MATLAB BLAS result for the C(1,1) spot.

Bottom line: Unless you can match all of the higher precision combined operations in the same way that the BLAS is doing the calculations, you should not expect to get the exact results as the BLAS library.