Vectorization problem - using vectors as inputs for another matrix?

Hello,

I have been struggling to vectorize the following code but need some help. The goal is to make this code run faster. I was able to vectorize part of the code, for example nested for loop, but wasn't able to vectorize some key sections such as

B(row)=B(row)-Ini(i,j-k+1)*temp;

where j and k are vectors of different size (see below). I would greatly appreciate some help. Thank you in advance!

N = 10;
Y = 100;
T= 60;
P = rand(N,T+1);
Ini = 10*rand(N,T+1);
A = sparse(Y*N*T, N*Y*(T+1));
B = zeros(Y*N*T,1);
for k=1:Y
  for j=1:T
    for i=1:N
      row= i+N*(j-1)+N*T*(k-1);
      if j>=k
        temp = 1;
        for m=j-k+1:j
          temp = temp*P(i,m);
        end
        B(row)=B(row)-Ini(i,j-k+1)*temp;
        start = N*Y;
        for m=1:k-1
          temp = 1;
          for n=m+j-k+1:j
            temp = temp*P(i,n);
          end
          col = start + i+N*(m+j-k-1)+N*T*(m-1);
          A(row,col) = A(row, col) - temp;
        end
        A(row,start +i+N*(j-1)+N*T*(k-1)) = A(row, start +i+N*(j-1)+N*T*(k-1)) -1;
      else %j<k
        % something  similar as in j>=k case
      end
    end
  end
end

[Matt J]

What is the code doing?

[Jan]

"for m=j-k+1;j" ??? Do you mean "for m=j-k+1:j" with a colon instead of a semicolon? This occurres twice.

[Tae Lim]

Matt - This code generates A and B matrix, which will be used for further analysis. This is a beginning portion of a much longer code.

Jan - You are correct. I've corrected the error.

Thanks!

[dpb]

Tae -- I think Matt is asking for a description of what A, B are, conceptually or general expression as it's apparent that the code is doing a lot of index manipulation to store elements into the two arrays. But, it's a lot of work for an outsider to try to decipher such code with no comments or other bread crumbs other than to find just syntax issues or the like. There might be clever algebraic manipulations that could be used or other specialized matrix constructs but to find those just by looking is tough...

[Tae Lim]

dpb - This is the beginning section of a stock-and-flow type of optimization model. A and B contains both current and past information of technologies of interest. In the original code, A and B sets the initial stock of technologies, which is then used to predict future stocks based on linear optimization and constraints. The else case (j<k) contains very similar code as in if case (j>=k) and can be vectorized in the same manner (and thus omission).

I hope this helps! I greatly appreciate your time and comments

[John BG]

Tae

this piece of code is faulty.

It's not about having it run faster, but first it shouldn't crash.

I say this because when ctrl+R the if clause inside the 3rd for loop, you can start getting some timing, but there the conditions in the if clause and the 4th for loop conditions keep it running endlessly.

Do you have a version of this code that stops?

John BG

<mailto:[e-mail] [e-mail]>

[dpb]

Well, w/o knowing which row/col are populated on the other case, not much (like anything) anybody can do with it at all and it pretty well stymies serious efforts overall since can't investigate overall patterns.

The "what" is interesting; what would be beneficial is the logic behind which are populated and the rules behind that instead of just trying to read code. As said, that's a hard nut to crack when that's all there is to go on.

Answer 1

[dpb]

Well, outside the above comment, starting from inside out,

temp = 1;
for n=m+j-k+1:j
  temp = temp*P(i,n);
end

can be replaced with

temp = prod(P(i,m+j-k+1:j);

You can begin working your way outwards from there to see what else can be done mechanicalistically, but the previous remark of needing to see the defining rules underling the implementation is probably the only way anybody can help too much by revising the algorithm. Now, that may or may not be possible, but think it's probably only real shot to make significant changes.

[Tae Lim]

This is very helpful, thank you! I have a further question on your comment. The revision you made works great with nested for loops but how can I can further vectorize this and avoid nested for loops? I want to replace nested for loops with:

k=1:Y
j=1:T
i=1:N
 % followed by vectorized codes....

This also relates to my question in the original post:

B(row)=B(row)-Ini(i,j-k+1)*temp;

Since both j,k arrays are used as inputs to another matrix, this requires different vectorization technique to avoid nested for loops. Thank you in advance for your time!

Answer 2

[Jan]

Most of the run time is used by the line

A(row,col) = A(row, col) - temp;

There is an MLint warning in the editor, that this kind of sparse indexing operations, which change the number of non-zero elements, is slow. Using full matrices instead would create a 27GB matrix, which is beyond my available RAM, such that I cannot compare the speed.

Using the suggested prod will improve the readability of the code (as an auto-indentation also), but has only tiny effects to the speed.

You cann accelerate the code by replacing

A = sparse(Y*N*T, N*Y*(T+1));

by

A = spalloc(Y*N*T, N*Y*(T+1), 1e7);

With Y=30 this reduces the runtime from 34 to 16 seconds.

[dpb]

"the suggested prod ... has only tiny effects to the speed."

Yeah, that's just the tip of the iceberg, so to speak. The spalloc to avoid reallocation as you point out is big.

BUT, fixing up the innermost loops is just the beginning; the end result may not be sparse at all depending on what is in the unprovided else block; the row index calculation of

row=i+N*(j-1)+N*T*(k-1);

is the same as

row=0;
for k
  for j
    for i
      row=row+1;
      ...

and only the j>=k clause prevents storing something every row for B, anyways. So, it may well be depending upon what B is for j<k that the two cases should be separate entirely and the j loop should just be for j=k:T and eliminate the branch.

Then, working thru the rest, one could likely compute the indexing beginning for each section being stored and write a vector expression for B; certainly the

B(row)=B(row)-

is superfluous as B is zero so a simple assignment is all that's needed for RHS.

Some pencil-pushing could figure out a lot here and using a small demo array to study patterns would help immeasurably; OP can do that better himself as he (at least should) know innately what the sparsity looks like (altho looks like probably will end up pretty dense in the end, irregardless).

It's a lot to expect of volunteers with no more assistance than provided as comments note....

[Tae Lim]

Thank you Jan this is extremely helpful!