Speeding up a code involving nested for loops

[EDIT: Wed Jun 1 16:10:09 UTC 2011 - Reformat - MKF]

The following is a simplified version of a code that I need to run many times. The 'rand' function is replacing calculations that are of the same order of complexity. Any intelligent way of converting the nested loops (and the multiply-sum operation)into matrix (or faster) operations would be greatly appreciated.

M=1000;
N=1000;

PSI_conv = zeros (M,N);
PSI_ab = rand(M,N);

for s = 1 : M  % M = 1000 
    for t = 1 : N  % N = 1000
        PHI_sph = rand(M,N);
        PSI_conv(s,t) = sum(PSI_sph(:).* conj(PSI_ab(:)));
    end
end

Answer 1

[Sean de Wolski]

The idea is to minimize the number of computations inside the FOR-loop.

In this case, conj(PSI_ab(:)) doesn't change and thus only needs to be computed once. Why bother generating PSI_sph as an MxN matrix when you could just generate it as a vector and then not need the (:) operation?

M=1000; 
N=1000;

PSI_conv = zeros (M,N); 
PSI_ab = rand(M,N);
PSI_ab = conj(reshape(PSI_ab,numel(PSI_ab),1));
MN = M*N;

for s = 1 : M % M = 1000 
  for t = 1 : N % N = 1000 
      PHI_sph = rand(1,MN); %Edit per Jan's comment and my time test.
      PSI_conv(s,t) = PSI_sph*PSI_ab; 
  end
end

[Jan]

"rand(1,MN) * PSI_ab" might be faster.

[Mehdi]

Thanks Sean and Jan. The idea of removing the repetetive conj operation is definitely a time saver, thanks. I will also try converting matrices to vectors to see how much gain I can get. I was hoping that there is a smart way of leveraging matrix multiplications (like Jan's comment) to avoid one (or both) of the for loops.

[Sean de Wolski]

You would be correct. I'm seeing it take around 57% of the time.

t1 = 0;

t2 = 0;

mcol = rand(500*500,1);

mrow = rand(1,500*500);

m2 = rand(500*500,1);

for ii = 1:100;

tic

R1 = mcol.*m2;

t1 = t1+toc;

tic

R2 = mrow*m2;

t2 = t2+toc;

end

t2/t1

%{

ans =

0.55224

0.57016

0.57103

%}

[Jan]

"rand(1,MN)*PSI_ab" is a scalar already, such that SUM can be omitted.

[Sean de Wolski]

Dur-da-dur. Time for another coffee :)

[Matt Fig]

Another case where the (or a?) one liner vectorization is slower!

PSI_conv = reshape(sum(bsxfun(@times,rand(M*N,N,M),conj(rand(M*N,1)))),M,N);

[Sean de Wolski]

The 1000^4 element matrix must take some serious time to construct. I wonder if a single for-loop and 1000^3 matrix on each iteration would be faster.

[Matt Fig]

I posted a version with BSXFUN in a single loop and it was slower, so I deleted it...