[Math] A System of Matrix Equations (2 Riccati, 1 Lyapunov)

matricesmatrix equationsoptimal controloptimization

Setup:

Let $\gamma \in(0,1)$, ${\bf F},{\bf Q} \in \mathbb R^{n\times n}$, ${\bf H}\in \mathbb R^{n\times r}$, and ${\bf R}\in \mathbb R^{r\times r}$ be given and suppose that ${\bf P}$,${\bf W}$,${\bf X}\in \mathbb R^{n\times n}$, and ${\bf K}$,${\bf L}\in \mathbb R^{n\times r}$ satisfy

\begin{align}
{\bf P} &={\bf F}({\bf I}_{n}-{\bf K} {\bf H}^\top){\bf P}{\bf F}^\top+{\bf Q}, \;\;\;\;\;\;\;\;\text{where}\;\;\;\;{\bf K}\equiv {\bf P} {\bf H}\left({\bf H}^\top{\bf P} {\bf H}+ \frac{1}{\gamma}{\bf R} \right)^{-1} \tag1\\[4ex]
{\bf W} &={\bf F}({\bf I}_{n}-{\bf L} {\bf H}^\top){\bf W}{\bf F}^\top+{\bf Q}, \;\;\;\;\;\;\;\;\text{where}\;\;\;\;{\bf L}\equiv {\bf W} {\bf H}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})^{-1} \tag2\\[4ex]
{\bf X} &={\bf K}{\bf H}^\top {\bf W}+({\bf I}_n-{\bf K}{\bf H}^\top){\bf F}{\bf X}({\bf I}_n-{\bf H} {\bf L}^\top){\bf F}^\top \tag3\\
{\color{white}X}
\end{align}

Moreover, assume that ${\bf P}$, ${\bf W}$, ${\bf R}$, and ${\bf Q}$ are symmetric, and ${\bf P}$, ${\bf W}$ and ${\bf X}$ are invertible.

Want to prove:

$${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}){\bf L} \tag4$$

Some ideas and comments:

  • In the scalar case (with $n=r=1$) what works is to subtract $(2)$ from $(1)$ which eliminates ${\bf Q}$, then use that to solve for ${\bf F}^2$, i.e. (suppressing the $^\top$ notation)
    $$ {\bf F}^2 = \frac{(\gamma{\bf H}^2{\bf P}+{\bf R})({\bf H}^2{\bf W}+{\bf R})({\bf P}-{\bf W})}{{\bf R}^2({\bf P}-{\bf W})+(1-\gamma){\bf H}^2{\bf R}{\bf P}{\bf W})},$$
    substitute this into $(3)$ and solve for ${\bf X}$ which yields
    $${\bf X} = \frac{(\gamma {\bf R}({\bf P}-{\bf W})+\gamma(1-\gamma){\bf H}^2{\bf P}{\bf W}}{(1-\gamma)(\gamma {\bf H}^2{\bf P}+{\bf R}))}. $$
    Finally, solving for ${\bf X}$ using $(4)$ yields the same thing.

  • Notice that the equations $(1)$ and $(2)$ are Riccati equations on ${\bf P}$ and ${\bf W}$ respectively, which means that there are known methods (described here) to solve for ${\bf P}$ and ${\bf W}$. Though I was not able to make use of those solutions.

  • Equation $(3)$ is similar to a Lyapunov equation on ${\bf X}$, which means that a vectorization method (described here) allows one to obtain an equation for ${\bf X}$.

  • I think the proof will come from a procedure similar to the one that works for the scalar case. That is, subtracting $(2)$ from $(1)$, using this equation to eliminate ${\bf F}$ from $(3)$, then showing that the resulting equation implies $(4)$.

  • If you can prove it under additional assumptions that could also be helpful.

Here is a simple Matlab code that allows you to test the result:

    gamma = 0.5;

    F = [2, 2, 3; 4, 5, 6; 7, 8, 9];  % n x n
    H = [1, 2; 3, 1; 2, 1];           % n x r
    R = [3 , 1; 1, 3];                % r x r, symmetric
    Q = [2, 0, 0; 0, 5, 0; 0, 0, 9];  % n x n, symmetric

    n = length(F); 
    r = length(R); 

    I = eye(n);

    P = I;
    W = I;
    X = I;
    for i=1:1000
        K = (P*H)/(H'*P*H+(R/gamma));
        P = F*(P-K*H'*P)*F'+Q;
        L = (W*H)/(H'*W*H+R);    
        W = F*(W-L*H'*W)*F'+Q;
        X = K*H'*W+(I-K*H')*F*X*(I-H*L')*F';
    end

    disp(['K - (gamma*I+(1-gamma)*X*W^(-1))*L = ',...
          num2str(sum(sum(abs(K-(gamma*I+(1-gamma)*X*(W^(-1)))*L))))])

This question is a boiled down version of another question: Forecast equivalence between two steady-state Kalman filters (with ${\bf L} \equiv {\bf L}_1$, ${\bf W} \equiv {\bf W}_{11}$, and ${\bf X} \equiv {\bf W}_{12}$).

Best Answer

I posting it as another answer because I don't want to delete anything in previous yet and the other option is to make it too long.

4) The equation for $X$ is Sylvester equation, but you right this is closely related to Lyapunov. Lyapunov equation is a special case of Sylvester equation.

5) Remove linear algebra tag, the problem is non linear.

6) Add optimization and oprimal control tags.

7) Since we figured out that the problem with your code is well known numerical instability of matrix inverse (I should have noticed it instantly). I tried to write a code which bypasses the problem, i.e. do not invert matrices explicitly. I've ended with the following code which also give more hints on where to look for additional conditions about the matrices. Hopefully these conditions will help to find to prove the conjecture.

Briefly about the algorithm,

  • use matlab's algorithm dare to get the ricatti matrices $P$ and $W$
  • To obtain the values of $K$ and $L$ use zero finding algorithm (fsolve) using their formula in ($1$) and ($2$) respectively.
  • Use matlab algorithm dlyap to obtain $X$
  • To verify the result solve perturb the required equation so that there no matrix inverse, i.e. $${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}){\bf L}$$ $${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}) {\bf W} {\bf H}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})^{-1} $$ $${\bf K}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}) {\bf W} {\bf H} $$ $${\bf K}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})=(\gamma{\bf W}+(1-\gamma) {\bf X} ) {\bf H} $$

possible hints for the conjecture The algorithm is more stable with respect to matrix inverse. The finding of $K,L$ is somewhat a problem thought because the initial guess is not guided. In general, it may be essential difficulty in arriving to correct values. However when I run the algorithm $1000$ times I get $>600$ fails of dare and only few fails of fsolve. I think dare may shed some light on the original problem.

Thus, I've found that dare is very capricious algorithm. More precisely it doesn't has a solution for any input, but for narrowed family of matrices. The exact conditions as described in matlab aren't clear to me, so I even asked another question with hope that somebody will give an answer. I've found several unanswered questions about algebraic Ricatti equations, so I would suggest to go to literature to discover under what conditions to $F,Q,H,R$ the condition that the solution to Ricatti exists. The same about Silverster\Lyapunov. I would also investigate the link between algebraic and differential equations, perhaps the reduction to DE make the prove easier, so worth a try.

The code

n=9;
r=5;

F=round(10*rand(n),0);
H=round(10*rand(n,r),0);
R= round(10*rand(r),0);
R=R+R';
Q=round(10*rand(n),0);
Q=Q+Q';

P = dare(F',H,Q,R/gamma);
W = dare(F',H,Q,R);

options=optimset('Algorithm',{'levenbergmarquardt',0.0005});%,'Display','iter');
[K,fvalK,exitflagK] = fsolve(@(x) F*(P-x*H'*P)*F'+Q-P,zeros(n,r),options);
[L,fvalL,exitflagL] = fsolve(@(x) F*(W-x*H'*W)*F'+Q-W,zeros(n,r),options);

if exitflagK<=0 || exitflagL<=0 
        error('try again') 
end

X = dlyap((I-K*H')*F,(I-H*L')*F',K*H'*W);

tmp = K*(H'*W*H+R)-(gamma*W+(1-gamma)*X)*H;   
err=norm(tmp(:),inf);
disp(['K(H''W*H+R) - (gamma*W+(1-gamma)*X)*H) = ', num2str(err)])
Related Question