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

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]

What within the for-loop depends on 'k' (the loop counter)?

Answer 1

[Matt J]

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

[Djabba]

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!

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

[Djabba]

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!

[Matz Johansson Bergström]

No problem and I was too slow to answer too ;-)

Answer 3

[Matt J]

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;

[Djabba]

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]

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?