CODEWITHSUNDEEP
performancematlabfor-loopvectorization
huynh khactuan
huynh khactuan

Reputation: 35

Vectorize a for loop in which the input depends on the output

I have a complicated problem about vectorization of a dependence for loop, I would like to have some helps from you.

Let's X1 be a vector with length n1, X2 be a vector with length n2, F1 be a N1xn1 matrix, F2 be a N2xn2 matrix, Q be a N1xN2 matrix, and the notation p... are index vectors. ntrapz is a function for trapezoidal numerical integration. I would like to compute the matrix Q as follows:

for i1=1:N1
    F1_13tmp=F1(i1, p1_13)';                       % ' 
    F1_13=F1_13tmp(:,ones(n2,1));
    for i2=1:N2
        F2_13 = F2(i2, p2_13);
        Q_13_13 = Q(p1_13, p2_13);
        Q(i1,i2) = Q(i1,i2) + 
              ntrapz(X2(p2_13), ntrapz(X1(p1_13)', Q_13_13.*F1_13).*F2_13);
    end
end

The problem is that updating Q(i1,i2) changes the value of Q_13_13 = Q(p1_13, p2_13) for the next iteration. I would like to know if we can vectorize a such for loop. If not, any idea to speed up the code?

Thank you in advance for your help.

Upvotes: 2

Views: 255

Answers (1)

Dennis Jaheruddin
Dennis Jaheruddin

Reputation: 21563

Vectorization is typically the way to go for parallel operations. However, in your case it appears that you have sequential operations and those are generally not suitable for vectorization.

I am not familiar with the function you use, but the only way I can think of to possibly speed up your calculation is writing it in a non-dependand form. Note that this is more of a math exercise than a programming exercise, here is a simple example:

Suppose the formula is equal to:

Y(1)=0.5;
Y(t)=0.3*Y(t-1)+epsilon(t);
% Let as assume epsilon is also known
epsilon = rand(9,1)

The simple way to calculate Y(10) would be to do 10 calculation steps.

for t=2:10
    Y(t)=0.3*Y(t-1)+epsilon(t);
end

The fastest way to calculate Y(10) would be to do something like this (not sure if it is correct, but should be close)

Y(10)=Y(1)*0.3^9 + sum(epsilon.^(1:9))

So to conclude: The only way to "vectorize" the loop, is if you can help the computer significantly and basically let it do a different calculation. If this is not possible you can of course always try to squeeze out a bit of extra speed by making it a mex file.

Upvotes: 0

Related Questions