Fast matrix computation of a 'riccati like' equation in a 'for' loop

4 次查看（过去 30 天）
Djabba
Djabba 2014-8-20
评论： Matz Johansson Bergström 2014-8-21
hello, i'd like to optimize the following computation :
for k=1:N
P = A*(P + Q)*A' + W
end
with N very large (>1e6) ; A,P,Q and W are [50x50] real matrices.
I found that the following code was slighty better (few %):
for k=1:N
Ptmp = P+Q;
Ptmp = A*tmp;
Ptmp = Ptm*A';
P = Ptmp + W;
end
I tried the 'mtimesx' function, but it was slower than the previous code (I used all the optiosn : matlab,speed, speedomp).
Is there a better way to do this computation? Thank you

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2014-8-20
编辑：Matt J 2014-8-20
Assuming nothing in the loop is k-dependent (see My Comment), then you can pre-compute some things
At=A.';
Z=A*Q*At+W;
for k=1:N
P = A*P*At + Z;
end
  1 个评论
Djabba
Djabba 2014-8-20
thanks for your comment & answer. Unfortunately, the A matrix is k dependent : I shall re-write my code as follow :
for k=1:N
Ptmp = P+Q;
Ptmp = A(k)*tmp;
Ptmp = Ptm*A(k)';
P = Ptmp + W;
end
Actually, A(k) is extracted from a 50x50xN matrix B with the following code :
A(k) = B(:,:,k)
So if I understand correctly your answer, the point is to pre-compute the Z 3D-matrix [50x50xN] in a way that would be equivalent to :
Z = B*Q*B'+W;
with B 3D-matrix & Q,W 2D-matrix. I don't know if it's feasible : I'll try with mtimesx. thanks again!

请先登录，再进行评论。

更多回答（2 个）

Matz Johansson Bergström
I would like to preface my answer by saying that your code is essentially as fast as it gets. Matlab is very fast when it comes to Matrix-Matrix mutliplication. It was designed to. Without knowing the structure of the matrices I can only assume that we can apply ordinary matrix laws to break them up.
In the first expression we have two multiplications and two additions to do in a loop. On my computer, your first code takes 37.492 s. The second code takes 37.412 s. If we avoid to recompute multiplications and transposes, we can get it down to 35 s by just breaking up the expression and moving the computation around. This is only a 7 percent speedup and I think that this is as good as it gets.
We have
P = A*(P + Q)*A' + W
P = (A*P + A*Q)*A' + W
P = A*P*A' + A*Q*A' + W
This is four multiplications and two additions, but since we know A*Q*A'+W is never changing, we store it before the iterations. The resulting expression requires two multiplications and one addition. So, we only saved one addition, but the time spent is accumulated for each iteration.
I also store A', which provides some slight speedup. I also avoid using Ptmp, because we don't need to copy the data for each iteration, we can use P directly.
The code
%Preprocessing:
At = A';
con = A*Q*At + W;
for k = 1:N
P = A*P;
P = P*At;
P = P + con;
end
I believe it is correct. I run it on 100 iterations becase by running it for 1e6 iterations P consists of Inf elements.
So, I only managed to remove one addition, but I believe it is a good as it gets. Hopefully someone else can find a better way to do it.
  2 个评论
Djabba
Djabba 2014-8-20
thanks for your answer. As I wrote in thy comment to Matt's answer, A is k-dependent. So I've to find a way to do the pre-computation with 3D-matrix (see my previous comment), and then I should be able to use your code. thanks again!

请先登录，再进行评论。

Matt J
Matt J 2014-8-20
编辑：Matt J 2014-8-20
If A is symmetric and everything is k-independent, you can do as follows. I find it to be 15 times faster than the for loop for N=1e6 and 50x50 input matrices.
Z=A*Q*At+W;
An=A^N;
[R,E]=eig(A);
e=bsxfun(@power,diag(E),0:N-1);
Rk=kron(R,R); Ek=e*e.';
Z=bsxfun(@times,Rk, Ek(:).')*(Z(:).'*Rk).';
Z=reshape(Z,size(A));
P=An*P0*An.' + Z;
  2 个评论
Djabba
Djabba 2014-8-20
A is not symmetric (that's why I need to do the transpose). I must admit I don't really understand your code, but it looks efficient!
Matt J
Matt J 2014-8-20
Well, never mind. It also assumes A is k-independent. I don't think there's much that can be done to accelerate your code unless there is simplifying structure in A, Q and W and their k-dependence. If there is, it would be advisable to update your post with that information. Are Q and W also k-dependent?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by