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

조회 수: 4(최근 30일)
Xiaohan Du 25 Feb 2018
Commented: Jan 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?

로그인 to comment.

### 채택된 답변

Jan 25 Feb 2018
Jan 님이 편집함. 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.

#### 댓글 수: 6

표시 이전 댓글 수: 3
John BG 25 Feb 2018
Hi Xiaohan Du
I have repeated Jan Simon's simulation
1.- with 10 times more data
2.- with just 2 matrix products
L = 10; % Number of loops
N = 10000;
M = rand(N);
P = rand(N, 3);
..
tic
for k=1:L
D=M*P;
C5=D'*D;
end
toc
and it's slightly faster than all the above so far
Elapsed time is 102.358362 seconds.
Elapsed time is 100.317063 seconds.
Elapsed time is 1.529945 seconds.
Elapsed time is 0.809471 seconds.
Elapsed time is 0.773689 seconds.
At the end, what consumes CPU time is the amount of operations, if same function can be performed with significant less operations, the time saving should show up.
John BG
Xiaohan Du 26 Feb 2018
Thank you, I tested your data, got roughly the same elapsed time with you.
Jan 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.

로그인 to comment.