Vectorization of a double for loop

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.

[dpb]

gbpi also a constant?

[Geoffrey Dillon]

Whoops, yes: gbpi is also a constant. Sorry!

[Jan]

"gbpi"? It is gbp1.

Answer 1

[Roger Stafford]

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

[Geoffrey Dillon]

Wow, that did it! Thank you so much for all of your help!

Answer 2

[Roger Stafford]

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

[Jan]

EDITED: The quote was not ', but a similar character. But not ´ . Strange.

[Roger Stafford]

@Jan: Thanks for the editing. I thought I typed in an ordinary, garden-variety single quote. Perhaps my computer is playing tricks on me.

[Geoffrey Dillon]

If M=2^5 and N=2^12, then BL1 should be 4127x4127. I placed a breakpoint right after the this new code and BL1 is now 8253x4127, which will eventually cause an error.

In addition, this new block of code takes significantly longer (10x) to execute than the double for loop!

[dpb]

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]

@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.