Matrix multiply result different from loops
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Davis Pan
el 18 de Ag. de 2020
Comentada: James Tursa
el 18 de Ag. de 2020
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)
0 comentarios
Respuesta aceptada
James Tursa
el 18 de Ag. de 2020
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 comentarios
James Tursa
el 18 de Ag. de 2020
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.
Más respuestas (0)
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!