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

Question

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?

Answer 1

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.