How to vectorize for loop of partial arrays?

I am tryig to accelerate this part of my code, as it consumes the most of the runtime. I posted a simplified snippet which has the core functionality. Is there any way to use arrayfun or bsxfun, as my attempts have failed so far. See code below. Thanks a bunch. Ravi

 X=zeros(length(k),length(k));

 for i=1:length(f)
     X=X+(a(:,i)*a(:,i)')/b(i);
 end

[Geoff Hayes]

Ravi - what are the dimensions of a? Is it just a two-dimensional array whereby you multiply each column by itself and then divide each column by ith element of b? But if that were true, then you wouldn't need to initialize X as a two dimensional array (it would just be a single column).

Please provide some details concerning your a, b and why X is initialized using the length of k but you iterate over the length of f.

[Ravi Goyal]

Geoff, X becomes a square matrix, for instance 4x4. However, as the function gets called many times, the matrix grows and can be up 80x80. I multiply each column of a with its transposed and dived that by a single element of vector b, which yields me the quadratic matrix X. Each time the matrix X is calculated for all frequencies f. a has the size(k,f). k can also change as the function is recursive and grows up to f, mostly to terminate earlier though as an optimum is found later on in the function evaluation solving X*k=b for a vector k.

Answer 1

[Lessmann]

Hi,

this is a vectorized version of your loop.

k = ones(1000,1);
b = 10*rand(25,1);
f = ones(25,1);
a = 10*rand(1000,25);

c = 1./repmat(b',length(k),1);
X = (a.*c)*a';