Fast product of large matrices

Hello,

I need to calculate a matrix product like the following a lot of times:

A=rand(2^18,500)
AT=transpose(A)
x=rand(2^18,1)
B=AT*spdiags(x)*A

This should give a 500x500 matrix B (above matrices are just a random example). I need to do this quite often for changing x but constant A and it would take too long if I do it as above, is there any way to do it faster? I think the expression comes from the sum

B(k,l)=sum(A(:,k).*A(:,l).*x)

written in matrix form.

Thanks for your help!

[KALYAN ACHARJYA]

Apart from your question, is the dimensionality OK for matrix multiplication?

[Matt J]

So A is not sparse? What do you plan to do with B downstream in your code?

Answer 1

[Matt J]

If see about a 25% speed-up when done as B=AT*(x.*A),

N=2^18;
A=rand(N,500);
AT=A.';
x=rand(N,1);

tic;
B=AT*spdiags(x,0,N,N)*A;
toc%Elapsed time is 1.594161 seconds.

tic;
B=AT*(x.*A);
toc%Elapsed time is 1.195049 seconds.

I don't think you can hope for better than 50% speed-up. In the very best case, when x=ones(N,1), the compute time you are stuck with is,

tic;
B=AT*A;
toc;%Elapsed time is 0.899292 seconds.

[Matt J]

If you have the Parallel Computing Toolbox and a decent grpahics card, you could get some speed-up by doing this using gpuArrays.