清空筛选条件
清空筛选条件

Vectorization of a double for loop

2 次查看（过去 30 天）
Geoffrey Dillon
Geoffrey Dillon 2017-11-9
评论： Geoffrey Dillon 2017-11-13
Hello all, I am trying to speed up the assembly of a matrix with a complicated structure. Right now, I have several sets of nested loops such as
for i=M+2:M+N-2
for j=i+1:M+N-1
BL1(j,i) = gamm/(gbp1*h(i)^(1-beta))*...
((j-i+3/2)^beta-2*(j-i+1/2)^beta+(j-i-1/2)^beta);
end
end
Here gamm, beta, are constants, h is a vector of length M+N and BL1 is a (dense) matrix of size (M+N-1) x (M+N-1). Is it possible to vectorize this? I've tried several things, but none seem to work.
Thanks for taking a look.

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

采纳的回答

Roger Stafford
Roger Stafford 2017-11-10
The following code should correct the error in my first answer. This version also avoids the repetition involved in your code for the expression
((j-i+3/2)^beta-2*(j-i+1/2)^beta+(j-i-1/2)^beta)
and hopefully it should therefore take less execution time.
BL1 = zeros(M+N-1);
T = (1:N-3)';
X = (T+3/2).^beta - 2*(T+1/2).^beta + (T-1/2).^beta;
for i = M+2:M+N-2
BL1(T(1:M+N-1-i)+i,i) = gamm/(gbp1*h(i)^(1-beta))*X(1:M+N-1-i);
end

更多回答（1 个）

Roger Stafford
Roger Stafford 2017-11-9
编辑：Roger Stafford 2017-11-10
You can save unneeded repetition by computing one of the factors outside the loop.
BL1 = zeros(M+N-1);
T = (1:M+N-1)';
X = (T+3/2).^beta - 2*(T+1/2).^beta + (T-1/2).^beta;
for i=M+2:M+N-2
BL1(T+i,i) = gamm/(gbp1*h(i)^(1-beta)) * X;
end
This answer is incorrect. Read my correction in the revised answer below. RAS
  5 个评论
显示 3更早的评论
dpb
dpb 2017-11-10
As written for the above M,N --> BL1 will be M+N-1 x M+N-2 --> 4127x4126 as those are the upper limits for the i,j loops respectively and the indices for BL1 are those indices identically. Every row below M+2 (=4098) will be identically zero and all columns diagonally to the left will also be zero beyond that row.
Is that the intent? But if it is to be square and 4127x4127, then the i upper index is wrong or the column subscript expression isn't correct.
Did you preallocate BL1 before timing the loop solution? Awaiting the answer to the above observation I've not looked at whether there's an efficient vectorization to be done or not; often the looping solution can be pretty effective if do the obvious.
Roger Stafford
Roger Stafford 2017-11-10
@Geoffrey: My apologies. I mistakenly read the line
for j=i+1:M+N-1
as though it had the parentheses:
for j=i+(1:M+N-1)
I have written a second “answer” which should correct this error.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by