Inexact methods for the low rank solution to large scale Lyapunov equations

Abstract

The rational Krylov subspace method (RKSM) and the low-rank alternating directions implicit (LR-ADI) iteration are established numerical tools for computing low-rank solution factors of large-scale Lyapunov equations. In order to generate the basis vectors for the RKSM, or extend the low-rank factors within the LR-ADI method, the repeated solution to a shifted linear system of equations is necessary. For very large systems this solve is usually implemented using iterative methods, leading to inexact solves within this inner iteration (and therefore to “inexact methods”). We will show that one can terminate this inner iteration before full precision has been reached and still obtain very good accuracy in the final solution to the Lyapunov equation. In particular, for both the RKSM and the LR-ADI method we derive theory for a relaxation strategy (e.g. increasing the solve tolerance of the inner iteration, as the outer iteration proceeds) within the iterative methods for solving the large linear systems. These theoretical choices involve unknown quantities, therefore practical criteria for relaxing the solution tolerance within the inner linear system are then provided. The theory is supported by several numerical examples, which show that the total amount of work for solving Lyapunov equations can be reduced significantly.

Appendices

Appendix A: Proof of Lemma 2.2

Proof

From $\omega \underline{H}_j=0$ we have for $k\le j$

which, for $k=j$, immediately leads to $f_j^{(j)}=-\omega _{1:j}/(h_{j+1,j}\omega _{j+1})$, the first equality in (2.21a). Similarly, it is easy to show that for $k<j$, a left null space vector $\hat{\omega }\in \mathbb {C}^{1\times k+1}$ of $\underline{H}_{k}$ is given by the first $k+1$ entries of the null vector $\omega$ of $\underline{H}_j$. Hence, $f^{(k)}_k=-\omega _{1:k}/(h_{k+1,k}\omega _{k+1})$ holds for all $k\le j$.

For computing $f_j^{(k)}= e_k^*H_j^{-1}$, $k<j$, we use the following partition of $H_j$ and consider the splitting $f^{(k)}_j = [u,y]$, $u\in \mathbb {C}^{1\times k}$, $y\in \mathbb {C}^{1\times k-j}$:

$$\begin{aligned}&e_k^*=\left[ \begin{array}{ccc} 0_{1,k-1}\big |1\big |&0_{1,j-k} \end{array}\right] =f^{(k)}_jH_j=[u,y] \left[ \begin{array}{c} \begin{array}{c} H_{k-1}\\ e_{k-1}^*h_{k,k-1}\\ 0_{j-k,k-1} \end{array}\Bigg | H_{:,k:j} \end{array}\right] \\&=\left[ \begin{array}{cc} u\left[ \begin{array}{c} H_{k-1}\\ e_{k-1}^*h_{k,k-1} \end{array}\right]&\Bigg |[u,y]H_{:,k:j}\end{array}\right] . \end{aligned}$$

This structure enforces conditions on [u, y], which we now explore. First, u has to be a multiple of $\omega _{1:k}$. Here we exploited that due to the Hessenberg structure, $\omega _{1:k}\underline{H}_{k-1}=0$. In particular,

$$\begin{aligned} u_{1:k-1}H_{k-1}+u_ke_{k-1}^*h_{k,k-1}=0, \end{aligned}$$

such that $u_kh_{k,k-1}e_{k-1}^*H_{k-1}^{-1}=-u_{1:k-1}$. Since $f^{(k-1)}_{k-1}=e_{k-1}^*H_{k-1}^{-1}$ we can infer $u_{1:k-1}=-u_kh_{k,k-1}f^{(k-1)}_{k-1}$ and, consequently, (2.21c) for $k>1$. Similarly, $[\omega _{1:k},y]$ has to satisfy

$$\begin{aligned} 0_{1,j-k}&=[\omega _{1:k},y]H_{:,k+1:j}=\omega _{1:k}H_{1:k,k+1:j}+yH_{k+1:j,k+1:j}=\omega H_{1:j+1,k+1:j}\\&=\omega _{1:k}H_{1:k,k+1:j}+\omega _{k+1:j+1}H_{k+1:j+1,k+1:j}, \end{aligned}$$

leading to

$$\begin{aligned} y=\omega _{k+1:j+1}H_{k+1:j+1,k+1:j}H_{k+1:j,k+1:j}^{-1}=\omega _{k+1:j} + \omega _{j+1} h_{j+1,j}e_{j-k}^*H_{k+1:j,k+1:j}^{-1}, \end{aligned}$$

where $e_{j-k}$ is a canonical vector of length $j-k$. Hence,

$$\begin{aligned} v_j^{(k)}=[\omega _{1:k},y]=\omega _{1:j}+[0_{1,k},[0_{1,j-k-1},h_{j+1,j}\omega _{j+1}]H_{k+1:j,k+1:j}^{-1}], \end{aligned}$$

leading to (2.21b). Finally, the normalization constant $\phi _{j}^{(k)}$ is obtained by the requirement $[u,y]H_{1:j,k}=1$. $\square$

Appendix B: Setup of example msd

The example msd represents a variant of one realization of the series of examples in [52]. Set $K:={\text {diag}}\!\left( I_3\otimes K_1,k_0+k_1+k_2+k_3\right) +k_+e_{3n}^T-e_{3n}k_+\in {\mathbb {R}}^{n_2\times n_2},~K_1:=\mathrm {tridiag}(-1,2,-1)\in {\mathbb {R}}^{n_1\times n_1}$ with $k_0=0.1$, $k_1=1$, $k_2=2$, $k_3=4$, $n_2=3n_1+1$, and $k_+:=[(k_1,k_2,k_3)\otimes e_{n_1}^T,0]^T$. Further, $M_1:={\text {diag}}\!\left( m_1I_{n_1},m_2I_{n_1},m_3I_{n_1}m_0\right) \in {\mathbb {R}}^{n_2\times n_2}$ with $m_0=1000$, $m_1=10$, $m_2=k_2$, $m_3=k_3$ and $D:=\alpha M+\beta (K+K(M^{-1}K)+K(M^{-1}K)^2+K(M^{-1}K)^3)+\nu [e_1,e_{n_1},e_{2n_1+1}][e_1,e_{n_1},e_{2n_1+1}]^T$, $\alpha =0.8$, $\beta =0.1$, $\nu =16$. In [52] a different matrix D was used involving a term $M^{\frac{1}{2}}\sqrt{M^{-\frac{1}{2}}KM^{-\frac{1}{2}}}M^{\frac{1}{2}}$ which was infeasible to set up in a large scale setting (the middle matrix square is a dense matrix). The version of D used here was similarly used in [10]. Construct the $n\times n$, $n:=2n_2$ block matrices $\hat{A}:=\left[ \begin{array}{ll} 0&{}\quad I_{n_2}\\ -K&{}\quad -D\end{array}\right]$, $\hat{M}:={\text {diag}}\!\left( I_{n_2},M_1\right)$ representing a linearization of the quadratic matrix pencil $\lambda ^2M_1+\lambda D+K$. The right hand side factor is set up as $\hat{B}=\left[ \begin{matrix} I_{2m}&{}0_{2m,n_2-m}&{}\left[ \begin{matrix}0_{m,n_2-m}&{}I_m\\ I_m&{}0_{m,n_2-m}\end{matrix}\right] \end{matrix}\right] ^T\in {\mathbb {R}}^{n\times 2m}$. In Sect. 4 we use $n_1=4000$, $m=2$ leading to $n=24{,}002$, $r=4$. Finally, $A:=P^T\hat{A}P$, $M:=P^T\hat{M}P$, $B:=P^T\hat{B}$ as in [29, Section 5.6], where P is a perfect shuffle permutation: This leads to banded matrices A, M, except for some sole off-band entries from the low-rank update of D, resulting in noticable computational savings when applying sparse direct solvers or computing sparse (incomplete) factorizations. An alternative to exploit the structure in $\hat{A},\hat{M}$ within RKSM, LR-ADI is described in, e.g., [10, 30].