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

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)])