MATLAB Answers

Vectorizing nested loops ---- indexing problem

2 views (last 30 days)
Hi all,
so I've been fiddling with some nested for loops a coworker wrote and I've got the inner loops vectorized now, but I'm really at a loss for how to vectorize the outer loop, here is where my code stands:
data = [1 4 2 ....... 5 2]; %n length series of random points
length = length(data);
start = 1;
endd = length;
kmax = endd-1;
size = endd-start;
corr = zeros(1,kmax);
sigma = zeros(1,kmax);
for k = 1 : kmax
Fi = 0;
Fj = 0;
Fiik = 0;
Fi2ik2 = 0;
Fiik = data(start:(endd-k))^2.*data((start+k):endd)^2;
Fi2ik2 = sum(Fiik.^2);
Fiik = sum(Fiik);
Fj = sum(data(start:(endd-k)));
Fi = sum(data((start+k):endd));
corr(k) = (size-k)*Fiik/(Fi*Fj);
sigma(k) = sqrt(Fi2ik2/(size-k)-Fiik^2/(size-k)^2/sqrt(size-k)*Fi*Fj/(size-k)^2);
I found this vectorizing Loops but it seems that the index matrices pA and pB are created uniformly, which is this case the indexing is a little more complicated and I need to use actual values (i realize I can use a logical matrix for this as long as I can get the indexing correct). Any help is much appreciated, as I've spent the last 3 hours trying to find a sneaky way to get this working. As of now, the computation with actual data sets can be rather lengthy.
P.S. I would like to avoid MEX if possible, that is why I post this question
Jan on 12 Oct 2012
Duplicate threads joined.

Sign in to comment.

Accepted Answer

Teja Muppirala
Teja Muppirala on 12 Oct 2012
The sorts of operations you are doing, where you sweep one vector through another, multiply and then add them, can be done very efficiently by using convolution (CONV):
data = randn(1,10000); %n length series of random points
len = length(data);
start = 1;
endd = len;
kmax = endd-1;
sze = endd-start;
% Make the necessary vectors
data2 = data.^2;
data2c = conv(data2(end:-1:1),data2);
data22c = conv(data2(end:-1:1).^2,data2.^2);
% Make some more necessary vectors
FiikVec = data2c(kmax:-1:1);
Fi2ik2Vec = data22c(kmax:-1:1);
FjVec = cumsum(data(1:end-1));
FjVec = FjVec(end:-1:1);
FiVec = cumsum(data(end:-1:2));
FiVec = FiVec(end:-1:1);
kVec = 1:kmax;
% The vectorized calculation
C = (sze-kVec).*FiikVec./(FiVec.*FjVec);
sigma = sqrt(Fi2ik2Vec./(sze-kVec)-FiikVec.^2./...
Josh Parks
Josh Parks on 2 Nov 2012
for reference, I used this solution and it is considerably faster than building your own vectorized code for large data sets (as it the latter options takes a huge amount of RAM and then peters out).
Thanks again Teja

Sign in to comment.

More Answers (2)

Matt J
Matt J on 11 Oct 2012
Edited: Matt J on 11 Oct 2012
I doubt it's worth vectorizing this loop, but here are some tips,
(1) We can help you better if you send us code that runs error-free. The code you posted does not, with obvious errors in lines like
length = length(data); %length used both as a function and a variable name
data(start:(endd-k))^2.*data((start+k):endd)^2; %matrix operation ^2 instead of element-wise .^2
(2) Avoid using variable names like 'size' which are also function names. This deprives you of the use of that function.
Aside from the above, you can improve the speed of the loop by avoiding repeat computations. Here is one suggestion:
for k = 1 : kmax
range1=data( start:(endd-k) );
range2= data( (start+k):endd );
Fiik = (range1.*range2).^2;
Fi2ik2 = sum(Fiik.^2);
Fiik = sum(Fiik);
Fj = sum(range1);
Fi = sum(range2);
corr(k) = szk*Fiik/(Fi*Fj);
sigma(k) = sqrt( Fi2ik2/szk - Fiik^2/Fi*Fj/szk^(5/2) );
  1 Comment
Josh Parks
Josh Parks on 11 Oct 2012
yes, apologies for the non-error-free code. I will resist the temptation of coding a question at midnight from now on. Thanks for the helpful reply in despite my blunder!

Sign in to comment.

Jan on 12 Oct 2012
Edited: Jan on 12 Oct 2012
[Copied from the duplicate thread]
No a vectorization, but some speedup:
data2 = data .* data; % Avoid squaring
for k = 1 : kmax
% Remove the useless: Fi = 0; Fj = 0; Fiik = 0; Fi2ik2 = 0;
Fiik = data2(start:(endd-k)) .* data2((start+k):endd);
Fiik2 = Fiik .^ 2;
Fi2ik2 = Fiik * Fiik'; % Implicite sum by BLAS dot product

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!