Efficiently calculating the trace of a matrix product

I have two NxN square matrices, A and B, and I would like to calculate the trace of AB. Since the trace of AB only depends on its diagonal elements, it should hypothetically not be necessary to compute all of AB, thereby reducing the amount of operations from N^3 to N^2. My question is twofold:

• Does calling tr(AB) in MATLAB automatically exploit this fact?
• If not, is there an efficient way of doing this that doesn't involve calling for loops?

Thanks!

Answer 1

[Matt J]

Bt=B.';
traceProduct = A(:).'*Bt(:);

[Jan]

Addition: We do not know, if Matlab's JIT acceleration exploits this fact, because it is not documented. You can check this with measuring some timings:

len = 100:200:3000;
t1 = zeros(1, numel(len));  t2 = t1;  t3 = t1;  % Pre-allocate
for k = 1:numel(len)
  A = rand(len(k));
  B = rand(len(k));
  t1(k) = timeit(@() A*B);
  t2(k) = timeit(@() trace(A*B));
  t3(k) = timeit(@() A(:).' * reshape(B.', [], 1));
end
loglog(len, t1, 'r', len, t2, 'b', len, t3, 'c');
legend({'A*B', 'trace(A*B)', 'A(:).''*Bt(:)'})
p1 = polyfit(log(len), log(t1), 1)
p2 = polyfit(log(len), log(t2), 1)
p3 = polyfit(log(len), log(t3), 1)

Which slopes do you see?

Result: trace(A*B) needs almost the same time as A*B and much longer than A(:).' * Bt(:) - at least in Matlab online. Please run this on a local machine also. Interpretation: The JIT does not abbreviate the computations in: trace(A*B).

[James Tursa]

I think the only way the JIT could improve upon this would be to not form the transpose explicitly in memory, but do the dot product with fancy indexing on either A or B. It would be a pretty simple mex routine.

[Matt J]

Another way to test if the trace product is JIT optimized is to compare both implementations on a large matrix,

A=rand(3000); B=A;

tic;
version1 =trace(A*B);
toc;
%Elapsed time is 0.956904 seconds.

tic;
 version2 = A(:).'*reshape(B.',[],1);
toc;
%Elapsed time is 0.068032 seconds.

It is pretty clear that the direct implementation is not optimized.

[James Tursa]

I get the same results as Matt on various versions. And even if there is some version (maybe future) of MATLAB that does this optimization, it is still only the physical transpose part that could beat the direct implementation above ... i.e., compiled code that avoids physically forming the transpose.