Vectorizing nested loops ---- indexing problem

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

[per isakson]

• 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

[Jan]

Duplicate threads joined.

Answer 1

[Teja Muppirala]

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);

[Josh Parks]

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]

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

Answer 2

[Matt J]

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

[Josh Parks]

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!

Answer 3

[Jan]

[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
  ...