Riccati solvers with S parameter

[pmli]

Feature

Low-rank solvers and the dense solver from pymess should be extended to support the S parameter.

Notes

The Riccati equation $$A X E^T + E X A^T - \left(E X C^T + S^T\right) R^{-1} \left(C X E^T + S\right) + B B^T = 0$$ can be written as $$\left(A - S^T R^{-1} C\right) X E^T + E X \left(A - S^T R^{-1} C\right)^T - E X C^T R^{-1} C X E^T + \left(B B^T - S^T R^{-1} S\right) = 0.$$ Therefore, the substitution

$$ \begin{align*} A & \to A - S^T R^{-1} C, \\ B B^T & \to B B^T - S^T R^{-1} S, \\ S & \to 0, \end{align*} $$

leads to an equivalent Riccati equation. The matrix $A - S^T R^{-1} C$ can be represented using a LowRankUpdatedOperator. The matrix $B B^T - S^T R^{-1} S$ can be written as

$$ \begin{bmatrix} B & S^T \end{bmatrix} \begin{bmatrix} I & 0 \\ 0 & -R^{-1} \end{bmatrix} \begin{bmatrix} B & S^T \end{bmatrix}^T, $$

so computing the SVD of $[B \ S^T] = U \Sigma V^T$ and the Cholesky decomposition of

$$ \Sigma V^T \begin{bmatrix} I & 0 \\ 0 & -R^{-1} \end{bmatrix} V \Sigma $$

could lead to an updated (low-rank) factorization.

[drittelhacker]

The last step could be tricky, as I assume R is positive definite and thus the central matrix is indefinite. However, all algorithms in pymess can be formulated for $LDL^T$ respecting this indefiniteness. So maybe this is worth a feature request upstream, if $LDL^T$ is not yet supported there.