Vectorization problem - using vectors as inputs for another matrix?
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
John BG
2017-3-27
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
2017-3-27
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.
回答(2 个)
dpb
2017-3-24
编辑:dpb
2017-3-24
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.
Jan
2017-3-24
编辑:Jan
2017-3-24
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
2017-3-25
"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....
另请参阅
类别
在 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!