Vectorizing nested loops ---- indexing problem

1 次查看（过去 30 天）
Josh Parks
Josh Parks 2012-10-11
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);
end
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.
Thanks,
Josh
P.S. I would like to avoid MEX if possible, that is why I post this question
  2 个评论
per isakson
per isakson 2012-10-11
  • length, size, corr, and sigma are names of functions. Avoid to use function names as names of variables.
  • The command , Fiik = data(...., causes an error. I cannot run the code
  • What should be the output of the vectorized code? corr and sigma?
  • There might not be a vectorized version

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Teja Muppirala
Teja Muppirala 2012-10-12
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./...
(sze-kVec).^2./sqrt(sze-kVec).*FiVec.*FjVec./(sze-kVec).^2);
  2 个评论
Josh Parks
Josh Parks 2012-10-13
编辑：Josh Parks 2012-10-13
electrical engineer teja? I had heard of convolution integrals but had never seen them implemented. Turns out this is exactly what my application needed...
Josh Parks
Josh Parks 2012-11-2
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

请先登录，再进行评论。

更多回答（2 个）

Matt J
Matt J 2012-10-11
编辑：Matt J 2012-10-11
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 );
szk=size-k;
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) );
end
  1 个评论
Josh Parks
Josh Parks 2012-10-11
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!

请先登录，再进行评论。

Jan
Jan 2012-10-12
编辑：Jan 2012-10-12
[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
...

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by