From: Lisa on
In my simulations these loops need to be run through almost 200000 times. With my clumsy programming I estimate that this will take almost 20 days! I would really appreciate some help in vectorising and hopefully optimising this chunck of code (below)

Thanks!

%Y, A and O are matrices
%c and r are vectors
%d is scalar

for i = 1:final
d = SMALL_VALUE;
for j = 1:final
Y(i,j)=A(i,j)*r(j);
d = d + Y(i,j);
end
for j = 1:final
O(i,j) = (c(i)/d) * Y(i,j);
end
end
for j = 1:final
d = SMALL_VALUE;
for i =1:final
d = d + (A(i,j)*r(j) - O(i,j))* A(i,j);
end
r(j)= r(j) - d/p(j);
end
From: Roger Stafford on
"Lisa " <lwillats(a)unimelb.edu.au> wrote in message <hvuged$omr$1(a)fred.mathworks.com>...
> In my simulations these loops need to be run through almost 200000 times. With my clumsy programming I estimate that this will take almost 20 days! I would really appreciate some help in vectorising and hopefully optimising this chunck of code (below)
>
> Thanks!
>
> %Y, A and O are matrices
> %c and r are vectors
> %d is scalar
>
> for i = 1:final
> d = SMALL_VALUE;
> for j = 1:final
> Y(i,j)=A(i,j)*r(j);
> d = d + Y(i,j);
> end
> for j = 1:final
> O(i,j) = (c(i)/d) * Y(i,j);
> end
> end
> for j = 1:final
> d = SMALL_VALUE;
> for i =1:final
> d = d + (A(i,j)*r(j) - O(i,j))* A(i,j);
> end
> r(j)= r(j) - d/p(j);
> end
- - - - - - - - - -
Y = A.*repmat(r,final,1); % r a row vector
O = Y.*repmat(c./(sum(Y,2)+SMALL_VALUE),1,final); % c a col. vector
r = r - (sum((Y-O).*A,1)+SMALL_VALUE)./p; % p a row vector

This assumes that r and p are row vectors and c is a column vector.

Roger Stafford
From: Matt Fig on
Did you pre-allocate your arrays? How large are your input arrays? If you did pre-allocate, I would be very surprised it this took 20 days to run. Your first FOR loop can be vectorized as follows (I assume c and r are row vectors):



Y2 = bsxfun(@times,A,r);
D = SMALL_VALUE + cumsum(Y2,2);
D = D(:,final);
O2 = bsxfun(@times,Y2,c.'./D);


I don't have time to do the second loop, but I still wonder about your speeds. If you had pre-allocated Y and O, your code should still be pretty fast - depending on the sizes involved.
From: Matt Fig on
Alright, so I do have time. This is similar to Roger's solution except for the use of bsxfun. If you don't need the final values of d, this could be simplified even further.


Y2 = bsxfun(@times,A,r2);
D = SMALL_VALUE + sum(Y2,2);
O2 = bsxfun(@times,Y2,c'./D);
D = sum((bsxfun(@times,A,r2) - O).*A) + SMALL_VALUE;
r2 = r2 - D./p;
From: Lisa on
Thanks for your time! This is really helpful. I forgot to mention in the orginal post that 'p' is a matrix, so p(j) refers to the jth column of the matrix. I assume this will affect the last line of your code (and Roger's)?

Cheers.

"Matt Fig" <spamanon(a)yahoo.com> wrote in message <hvul0j$i6l$1(a)fred.mathworks.com>...
> Alright, so I do have time. This is similar to Roger's solution except for the use of bsxfun. If you don't need the final values of d, this could be simplified even further.
>
>
> Y2 = bsxfun(@times,A,r2);
> D = SMALL_VALUE + sum(Y2,2);
> O2 = bsxfun(@times,Y2,c'./D);
> D = sum((bsxfun(@times,A,r2) - O).*A) + SMALL_VALUE;
> r2 = r2 - D./p;