# Difference between these 2 operations: D = M' * M then P' * D * P and P' * (M' * M) * P.

5 views (last 30 days)
Xiaohan Du on 25 Feb 2018
Commented: Jan on 26 Feb 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?

Jan on 25 Feb 2018
Edited: Jan on 25 Feb 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.
Jan on 26 Feb 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.