MATLAB Answers

Vectorization problem - using vectors as inputs for another matrix?

조회 수: 4(최근 30일)
Tae Lim
Tae Lim 24 Mar 2017
댓글: Tae Lim 27 Mar 2017
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

  댓글 수: 7

표시 이전 댓글 수: 4
Tae Lim
Tae Lim 27 Mar 2017
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
John BG 27 Mar 2017
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
dpb
dpb 27 Mar 2017
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.

Sign in to comment.

답변(2개)

dpb
dpb 24 Mar 2017
편집: dpb 24 Mar 2017
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.

  댓글 수: 1

Tae Lim
Tae Lim 27 Mar 2017
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!

Sign in to comment.


Jan
Jan 24 Mar 2017
편집: Jan 24 Mar 2017
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.

  댓글 수: 2

dpb
dpb 25 Mar 2017
"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
Tae Lim 27 Mar 2017
Thank you Jan this is extremely helpful!

Sign in to comment.

태그

제품


Translated by