This paper deals with the numerical computation of Lagrangian invariant subspaces of large-scale Hamiltonian matrices in the context of the solution of Lyapunov equations.
Given matrices \(F,W=W^T\in\mathbb{R}^{n\times n}\), a Lyapunov equation is a linear matrix equation of the form \[ F^TX+XF+W=0, \] to be solved for a symmetric matrix \(X\). This equation is closely connected to the algebraic Riccati equations \[ F^TX+XF+W-XGX=0, \] where \(G=G^T\in\mathbb{R}^{n\times n}\), to be solved for \(X=X^T\) such that \(F-GX\) has all its eigenvalues in the open left half plane.
For small and medium-size problems (\(n\approx 1000\)), the most usual methods in practice for solving the Riccati equation use a reformulation as an invariant subspace problem: Find \(S\in\mathbb{R}^{2n\times n}\) with full column rank and \(T\in\mathbb{R}^{n\times n}\) with eigenvalues in the left half plane such that \(HS=ST\) with \[ H=\begin{bmatrix} F&G\\ -W&-F^T\end{bmatrix}. \] The solution of this problem is related to the Riccati equation: Setting \[ S=\begin{bmatrix} S_1\\ S_2\end{bmatrix},\;\;\text{with}\;\;S_1,S_2\in\mathbb{R}^{n\times n} \] the solution of the Riccati equation is \(X=S_2S_1^{-1}\). The difference in accuracy between the approaches based on the invariant subspace problem \(HS=ST\) and those based on solving Riccati equation directly is mostly evident in the case in which matrix \(S_1\) is ill-conditioned or almost singular, corresponding to a solution \(X\) with very large entries. This situation is described as close to uncontrollability.
On the other hand, a property that holds in many cases is that the eigenvalues of \(X\) decay very rapidly, so that \(X\) can be approximated accurately by a low-rank matrix \(ZZ^T\), \(Z\in\mathbb{R}^{n\times k}\), \(k\ll n\). So, \[ X\approx ZZ^T,\;\mathcal{S} = \operatorname{im}S,\;\text{with}\;S\approx\begin{bmatrix} I\\ ZZ^T\end{bmatrix}, \] is a data-sparse format to represent and compute the solution.
In this paper, the authors consider the low-rank alternating direction iteration (LR-ADI) method for Lyapunov equations and propose a way to combine it with ideas from inverse-free techniques to represent a basis of invariant subspace \(\mathcal{S}=\operatorname{Im}\begin{bmatrix} I\\ X\end{bmatrix}\) that is well-conditioned even if \(\|X\|>1\) and the problem is close to uncontrollability. The inverse-free representation uses a pair \((A,E)\) such that \(Z=AE^{-1}\) in the low-rank factor \(Z\) of \(X\), and they show how to adapt the main LR-ADI loop to work with it.
At termination of the LR-ADI, the authors show how \((A,E)\) can be used to construct directly a data-sparse representation of \(\mathcal{S}\) \[ \mathcal{S} = \operatorname{im}U,\;\text{with}\;U=\begin{bmatrix} I_n-U_1U_2^T\\ U_3U_4^T\end{bmatrix}, \] where \(U_1,U_2,U_3,U_4\in\mathbb{R}^{n\times k}\) with \(k\ll n\), for which the condition number \(\kappa(U)\) can be bounded explicitly. This achieves both goals of data sparsity and stability in the case of close to uncontrollable problems.
Finally, a partial error analysis is presented and some numerical examples are provided to confirm the theoretical results and to compare the proposed methods with other known ones.