Block Diagonal matrix multiplication with a sparse matrix

Hello,

I have a matrix A of type [W, 0; 0 W], with W a 1600x1600 complex full matrix and 0 a 1600x1600 matrix of zeros, and another matrix B which is a sparse 3200x3200 matrix. I want to perform A' * B * A , but the time to compute is around 19 ~ 20 seconds, which for my purpose is too slow. Would there be any way in which I could compute this operation faster? I noticed that when computing A' * B, it takes ~ 1 second, but one the next computation has to be done, it takes around 18 sec.

Thank you.

Answer 1

[Teja Muppirala]

Although (A'*B) can be done quickly, the answer is no longer sparse. It is just some complex matrix. So then you are basically multiplying two full 3200x3200 complex matrices which will necessarily take a long time. There is no way around this. You can use the fact that A = [W 0; 0 W] to speed things up by roughly a factor of two, by doing the multiplication blockwise.

%%Make some random A and B
X = randn(1600) + 1i*randn(1600);
A = blkdiag(X,X);
B = sprand(3200,3200,0.01);

%%Do the multiplication blockwise
tic;
W = A(1:1600,1:1600);
B11 = B(1:1600,1:1600);
B12 = B(1:1600,1601:3200);
B21 = B(1601:3200,1:1600);
B22 = B(1601:3200,1601:3200);
F1 = [W'*B11*W W'*B12*W; W'*B21*W W'*B22*W];
toc;

tic; F2 = A'*B*A; toc;

%%Confirm that both yield the same answer
E = F1-F2;
max(abs(E(:))) %<-- Very small

If B has some special characteristics, you maybe be able to exploit them to improve this further.