Daniel B Szyld
Classical iterative methods for the solution of Generalized Lyapunov Equations

Department of Mathematics
Temple University
1804 N Broad St
Philadelphia
PA 19122
szyld@temple.edu

There has been a flurry of activity in recent years in the area of solution of matrix equations. In particular, a good understanding has been reached on how to approach the solution of large scale Lyapunov equations; see, e.g., [1]-[3]. An effective way to solve Lyapunov equations of the form $A^{T}X + XA +C^{T}C = 0$ , where $A$ and $X$ are $n \times n$ , is to use Galerkin projection with appropriate extended or rational Krylov subspaces. These methods work in part because the solution is known to be symmetric positive definite with rapidly decreasing singular values, and therefore it can be approximated by a low rank matrix $X_k= Z_k Z_k^T$ . Thus the computations are performed usually with storage which is lower rank, i.e., much lower than order of $n^2$ .

Generalized Lyapunov equations have additional terms. In this paper, we concentrate on equations of the following form

$$A^{T}X + XA + \sum_{j=1}^m N_j X N_j^T + C^{T}C = 0,$$ (1)

Such equations arise for example in stochastic control [4]. It has been proposed, e.g., in [5], to use approximations to conjugate gradients or to BiCGStab, appropriately preconditioned, where the basis vectors (matrices) and iterates are ``truncated" throughout the algorithm to keep all these elements represented by low-rank matrices.

In the present work, we propose a return to classical iterative methods, and consider instead stationary iterations, where the iterate $X_{k+1}= Z_{k+1} Z_{k+1}^T$ is the solution of

$$P(X):= A^{T}X + XA = Q(X_k),    k=0,1, \ldots$$ (2)

where $Q(X) = - \sum_{j=1}^m N_j X N_j^T - C^{T}C$ ; cf. [6]. The classical theory of splittings applies here, while imposing certain conditions on $P(X)$ and $Q(X)$ .

One of the advantages of this classical approach is that only the data and the low-rank factors of the old and new iterates $X_k$ and $X_{k+1}$ need to be kept in storage. Furthermore, the solutions of the Lyapunov equations () can be performed with the Galerkin projection methods mentioned above, where the growth of rank can usually be well contained.

In our work, we developed a general srategy for augmented Krylov projection methods for sequences of Lyapunov equations with converging right-hand sides, and we apply it to the sequence (), whose solutions converge to the solution of ().

Numerical experiments show the competitiveness of the proposed approach.

This is joint work with Stephen D. Shank and Valeria Simoncini.