Vectorization of a double for loop
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Geoffrey Dillon
el 9 de Nov. de 2017
Comentada: Geoffrey Dillon
el 13 de Nov. de 2017
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.
3 comentarios
Respuesta aceptada
Roger Stafford
el 10 de Nov. de 2017
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
Más respuestas (1)
Roger Stafford
el 9 de Nov. de 2017
Editada: Roger Stafford
el 10 de Nov. de 2017
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 comentarios
dpb
el 10 de Nov. de 2017
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
el 10 de Nov. de 2017
@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.
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!