transform from "for loop" into "matrix" form

Hi all,

I have a complicated problem about "for loop" and "matrix", 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, and Q be a N1xN2 matrix. I like compute the matrix Q as follows:

for i1=1:N1

for i2=1:N2
        Q(i1,i2)=Q(i1,i2)+trapz(X2(p2(1):p2(3)),trapz(X1(p1(1):p1(3)),Q(p1(1):p1(3),p2(1):p2(3))...
                                                        .*repmat(F1(i1,p1(1):p1(3))',1,p2(3)-p2(1)+1)).*F2(i2,p2(1):p2(3)));

end

end

where p1 and p2 are index vectors.

This program with "for loop" gives me a good result but takes time. I would like to transform it into "matrix" form to accelerate the computation.

The problem is that the output Q is taken as input in the body of "for loop". If it is not taken as the input, I can transform the above "for loop" as follows:

I=squeeze(trapz(X1(p1(1):p1(3)),repmat(Q(p1(1):p1(3),p2(1):p2(3)),[1 1 N1])...
                                                 .*rot90_3D(rot90_3D(repmat(F1(:,p1(1):p1(3)),[1 1 p2(3)-p2(1)+1]),3,3),1,3)));

Q=Q+squeeze(trapz(X2(p2(1):p2(3)),repmat(I,[1 1 N2]).*rot90_3D(rot90_3D(repmat(F2(:,p2(1):p2(3)),[1 1 N1]),3,3),1,3)));

It works but it is not exactly what I want.

So, do you have any idea to deal with the problem? Thank you in advance for you help.

Tuan

Answer 1

[Jan]

Start with simplifying the loops (Q is pre-allocated, isn't it?!):

p11 = p1(1);
p13 = p1(3);
p21 = p2(1);
p23 = p2(3);
c1 = X2(p21:p23);
c2 = p23 - p21 + 1;
c3 = F2(i2, p21:p23);
c4 = X1(p11:p13);
c5 = repmat(F1(i1, p11:p13)', 1, c2);

for i1=1:N1
  for i2=1:N2
    v = Q(p11:p13, p21:p23);
    Q(i1,i2) = Q(i1,i2) + trapz(c1, trapz(c4, v .* c5) .* c3);
  end
end

Is this faster already?

I guess that the most time is spent in trapz() such that optimizing or vectorizing the loop will not help very much.

[Tuan]

Dear Jan Simon,

Thank you so much for your reply. I have simplified the code according to your suggestion. Since c3 and c5 depend on the index i2 and i1, we have to take them into the loops. And the code that I have tried is as follows

 p2321 = p23 - p21 + 1;
 c1 = X2(p21:p23);
 c4 = X1(p11:p13);
 for i1=1:N1
     c5 = repmat(F1(i1, p11:p13)', 1, p2321);
     for i2=1:N2
         c3 = F2(i2, p21:p23);
         v = Q(p11:p13, p21:p23);
         Q(i1,i2) = Q(i1,i2) + trapz(c1, trapz(c4, v .* c5) .* c3);
     end
 end

It works and is faster than the original one around 30%.

As in your comment, the most time is spent in trapz. This is why I would like to transform into "matrix" form. Because with "for loop" we have to do 2*N1*N2 trapz operation, while with "matrix" form, we just do it 2 times. But the problem is always that the output Q (partial) is taken as input in the body of "for loop". And I do not know how to transform it into "matrix" form.

Thank you once again for your help and hope to receive your comments.

Tuan

[Jan]

Reducing the number of trapz calls would be an advantage. Creating a large temporary matrix as input for trapz would be a disadvantage.

But updating Q(i1,i2) can change the value of v=Q(p11:p13, p21:p23) for the next iteration in theory. Therefore a matrix version is not possible even for the inner loop.

[Tuan]

Thank you!