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

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

[Jan]

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.

[Xiaohan Du]

but Matlab warns me if I use

C = P' * M' * M * P

parenthesize the multiplication of M and its transpose to ensure the result is Hermitian

So is it just a different rounding?

[Jan]

I do not get such a hint in R2016b. Maybe a newer Matlab version recognizes the symmetry of the formula and can apply a specific BLAS/LAPACK method.

See the modified code above, which displays the difference of the results and run the tests with your Matlab version also. What do you get as output?

[Xiaohan Du]

Tested in R2017b using your code, speed test results are:

Elapsed time is 0.160442 seconds.
Elapsed time is 0.159922 seconds.
Elapsed time is 0.009256 seconds.
Elapsed time is 0.007466 seconds.

tested accuracy using:

norm((C - C1) ./ C, 'fro')
norm((C - C2) ./ C, 'fro')
norm((C - C3) ./ C, 'fro')
norm((C - C4) ./ C, 'fro')

results are:

ans =
     0
ans =
     0
ans =
   3.2822e-16 
ans =
   7.4978e-16

It seems that test 3 and 4 results are slightly different (rounding effect? not a wrong result?), while in test 3, Matlab warns me:

parenthesize the multiplication of M and its transpose to ensure the result is Hermitian

It's weird you do not have this warning, cause I found this post from 2012: https://stackoverflow.com/questions/10724549/warning-parethesize-the-multiplication-of-d-and-its-transpose-to-ensure-the

[John BG]

Hi Xiaohan Du

this is John BG <mailto:[e-mail] [e-mail]>

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.

I think Jan Simon deserves this answer accepted, would you consider marking his answer as the accepted answer to your question?

John BG

<mailto:[e-mail] [e-mail]>

[Xiaohan Du]

Thank you, I tested your data, got roughly the same elapsed time with you.

[Jan]

@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.