Get the diagonal without calculating the explicit matrix

Dear all:

I am trying to calculate a diagonal of a matrix (denoted A), which is formed by multiplying two large-dimensional matrices (denoted as B*C).

A naive way to do it is: first, calculating explicitly A = B*C, then get diagonal out from A. However, the first step takes forever to run due to the high-dimension of B and C. But the only thing I need is the diagonal of A.

Another straightforward way in my mind is: I could create a loop by calculating each element of the diagonal of A one by one. It will surely save a lot of time, but I am not sure if this is the most efficient way.

I am wondering if anyone knows a faster/smarter way to calculate it.

Thank you very much in advance!

Best,

Long

[Matt J]

The best approach will depend on the dimensions of the matrices, and whether they are of sparse-type or not.

Answer 1

[Matt J]

Assuming B*C results in a square matrix,

diagonal=sum(B.' .* C, 1);

[Matt J]

If B*C is not square, the same concept applies, just a bit messier:

N=min(size(B,1), size(C,2));

diagonal=sum(B(1:N,:).' .* C(:,1:N), 1);

[Long Hong]

Thank you Matt - The algorithm is simple! However, it does not seem faster than doing the loop. Maybe it is because making a transpose of a high-dimensional matrix takes a lot of time. What do you think?

[Matt J]

Or perhaps the outer dimensions of your matrices are much smaller than the inner dimensions.

[Matt J]

Note that if your matrices are sparse, the looping method will be much slower.

N=7000;
B=sprand(N,10*N,10/N); C=B.';


tic; sum(B.'.*C,1); toc  %Elapsed time is 0.021923 seconds.

tic;
for i=1:N
   B(i,:)*C(:,i); 
end
toc; %Elapsed time is 11.467097 seconds.

[the cyclist]

If the transpose in Matt's algorithm is really the time-consuming step, would it be possible to go upstream in your code and do all prior calculations in a way that the B you end up with is the transpose of the current one?

[Long Hong]

Thank you Matt! I have tested it in a relatively large subset of my original data, your algorithm is indeed faster. Thank you for the valuable advice!

[Long Hong]

Thank you the cyclist! Do you have any insight in doing the transpose quicker? I am a bit confused here.

[Matt J]

the cyclist means you might avoid the transpose by loading data column-wise instead of row-wise when you first build B.

Answer 2

[the cyclist]

Here is one way:

% Make up some inputs
N = 4;
B = rand(N);
C = rand(N);

% Calculate the diagonal
A_diag = 0;
for nr = 1:N
    A_diag = A_diag + B(:,nr).*C(nr,:)';
end

[Long Hong]

Thanks the cyclist! This is a method I have applied currently.