transform from "for loop" into "matrix" form

Question

Tuan 2012-11-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/53313-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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Jan 2012-11-9

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/53313-transform-from-for-loop-into-matrix-form#answer_64946

编辑：Jan 2012-11-9

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Tuan 2012-11-10

编辑：Tuan 2012-11-10

在 MATLAB Online 中打开

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 2012-11-10

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 2012-11-10

Thank you!

请先登录，再进行评论。

transform from "for loop" into "matrix" form

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

transform from "for loop" into "matrix" form

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论