Matrix multiply result different from loops

1 次查看（过去 30 天）
Davis Pan
Davis Pan 2020-8-18
评论： James Tursa 2020-8-18
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)

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

采纳的回答

James Tursa
James Tursa 2020-8-18
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.
  2 个评论
Davis Pan
Davis Pan 2020-8-18
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
James Tursa 2020-8-18
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by