Difference between these 2 operations: D = M' * M then P' * D * P and P' * (M' * M) * P.
7 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Xiaohan Du
el 25 de Feb. de 2018
Comentada: Jan
el 26 de Feb. de 2018
Hi all,
Imagine I have
M = rand(N, N);
P = rand(N, 3);
where N is a VERY large scalar, then for these 2 functions: 1.
function C = test1(M, P)
D = M' * M;
C = P' * D * P;
and 2.
function C = test2(M, P)
C = P' * (M' * M) * P;
so the output C are both 3-by-3 matrices. Is the 2nd function more efficient and cost less storage than the 1st one? I think so because in the 2nd fucntion I do not physically compute M' * M, is this correct?
0 comentarios
Respuesta aceptada
Jan
el 25 de Feb. de 2018
Editada: Jan
el 25 de Feb. de 2018
No, the assumption is not correct. The matrix in the parentheses is calculated and stored as a temporary array. There is no way to avoid the explicit calculation and in my speed measurements both takes exactly the same time.
By the way, this is 8 times faster:
function C = test3(M, P)
C = P' * M' * M * P;
And even faster:
function C = test4(M, P)
D = M * P;
C = D' * D;
The result differ by rounding effects only.
Some code for testing:
L = 10; % Number of loops
N = 1000;
M = rand(N, N);
P = rand(N, 3);
% Warm up! Ignore!
for k = 1:L
C = P' * (M' * M) * P;
end
tic;
for k = 1:L
D = M' * M;
C1 = P' * D * P;
end
toc
tic;
for k = 1:L
C2 = P' * (M' * M) * P;
end
toc
tic;
for k = 1:L
C3 = P' * M' * M * P;
end
toc
tic;
for k = 1:L
D = M * P;
C4 = D' * D;
end
toc
Results:
Elapsed time is 1.596403 seconds.
Elapsed time is 1.606413 seconds.
Elapsed time is 0.159568 seconds.
Elapsed time is 0.071261 seconds.
Check the results:
(C - C1) ./ C
(C - C2) ./ C
(C - C3) ./ C
(C - C4) ./ C
All relative differences are in the magnitude of EPS - for some random input data. This might be different for your case, so check this.
6 comentarios
Jan
el 26 de Feb. de 2018
@John BG: Your C5 is exactly the same as my C4: D=M*P;C4=D'*D;.
@Xiaohan Du: A relative difference of 3.2822e-16 is tiny and cannot be avoided and caused by rounding, but the results are correct. A mathematical proof is easy: The matrix multiplication is associative: (A*B)*C == A*(B*C). And (A'*B')'==B*A. Numerically there are some differences due to rounding. A test with your real values is useful to estimate the true rounding effects. If C4 differs from C1 by much more than EPS(max(C(:)), this does not mean that either C1 or C4 is wrong, but only, that the calculations are highly sensitive. It is a valuable strategy to program both methods and compare the results to estimate roughly, how susceptible the calculations are for rounding effects.
The warning message (which is magically missing on my R2016b - I'm going to ask the technical support) means, that A'*(B'*B)*A is computed such, that the result is Hermitian. If you really need this property, you can assure it cheaper by using C4 and: C=0.5*(C + C') - but I'd expect C4 to be (guaranteed to be) Hermitean already, while C3 is not.
Más respuestas (0)
Ver también
Categorías
Más información sobre Test Execution 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!