Computationally enhanced projection methods for symmetric Sylvester and Lyapunov matrix equationsVersion of February 1, 2017. This research is supported in part by the FARB12SIMO grant of the Università di Bologna, and in part by INdAM-GNCS under the 2015 Project Metodi di regolarizzazione per problemi di ottimizzazione e applicazioni.

# Computationally enhanced projection methods for symmetric Sylvester and Lyapunov matrix equations††thanks: Version of February 1, 2017. This research is supported in part by the FARB12SIMO grant of the Università di Bologna, and in part by INdAM-GNCS under the 2015 Project Metodi di regolarizzazione per problemi di ottimizzazione e applicazioni.

Davide Palitta Dipartimento di Matematica, Università di Bologna, Piazza di Porta S. Donato, 5, I-40127 Bologna, Italy (davide.palitta3@unibo.it).    Valeria Simoncini Dipartimento di Matematica, Università di Bologna, Piazza di Porta S. Donato, 5, I-40127 Bologna, Italy (valeria.simoncini@unibo.it), and IMATI-CNR, Pavia, Italy.
###### Abstract

In the numerical treatment of large-scale Sylvester and Lyapunov equations, projection methods require solving a reduced problem to check convergence. As the approximation space expands, this solution takes an increasing portion of the overall computational effort. When data are symmetric, we show that the Frobenius norm of the residual matrix can be computed at significantly lower cost than with available methods, without explicitly solving the reduced problem. For certain classes of problems, the new residual norm expression combined with a memory-reducing device make classical Krylov strategies competitive with respect to more recent projection methods. Numerical experiments illustrate the effectiveness of the new implementation for standard and extended Krylov subspace methods.

S

ylvester equation, Lyapunov equation, projection methods, Krylov subspaces

{AMS}

47J20, 65F30, 49M99, 49N35, 93B52

## 1 Introduction

Consider the Sylvester matrix equation

 AX+XB+C1CT2=0,A∈Rn1×n1,B∈Rn2×n2,C1∈Rn1×s,C2∈Rn2×s (1)

where are very large and sparse, symmetric negative definite matrices, while are tall, that is . Under these hypotheses, there exists a unique solution matrix . This kind of matrix equation arises in many applications, from the analysis of continuous-time linear dynamical systems to eigenvalue problems and the discretization of self-adjoint elliptic PDEs; see, e.g., [1], and [20] for a recent survey. Although and are sparse, the solution is in general dense so that storing it may be unfeasible for large-scale problems. On the other hand, under certain hypotheses on the spectral distribution of and , the singular values of present a fast decay, see, e.g., [16], thus justifying the search for a low-rank approximation to so that only these two tall matrices are actually computed and stored. To simplify the presentation of what follows, from now on we will focus on the case of the Lyapunov matrix equation, that is () and , so that will be square, symmetric and positive semidefinite [21]. In later sections we will describe how to naturally treat the general case with and distinct and not necessarily with the same dimensions, and different .

For the Lyapunov equation, projection methods compute the numerical solution in a sequence of nested subspaces, , . The approximation, usually denoted by , is written as the product of low-rank matrices where and the columns of are far fewer than . The quality and effectiveness of the approximation process depend on how much spectral information is captured by , without the space dimension being too large. The matrix is determined by solving a related (reduced) problem, whose dimension depends on the approximation space dimension. To check convergence, the residual matrix norm is monitored at each iteration by using but without explicitly computing the large and dense residual matrix [20]. The solution of the reduced problem is meant to account for a low percentage of the overall computation cost. Unfortunately, this cost grows nonlinearly with the space dimension, therefore solving the reduced problem may become very expensive if a large approximation space is needed.

A classical choice for is the (standard) block Krylov subspace [8], whose basis can be generated iteratively by means of the block Lanczos procedure. Numerical experiments show that may need to be quite large before a satisfactory approximate solution is obtained [15],[19]. This large number of iterations causes high computational and memory demands. More recent alternatives include projection onto extended or rational Krylov subspaces [19],[6], or the use of explicit iterations for the approximate solution [15]; see the thorough presentation in [20]. Extended and more generally rational Krylov subspaces contain richer spectral information, that allow for a significantly lower subspace dimension, at the cost of more expensive computations per iteration, since system solves with the coefficients matrix are required at each iteration.

We devise a strategy that significantly reduces the computational cost of evaluating the residual norm for both and the extended Krylov subspace . In case of a “two-pass” strategy is implemented to avoid storing the whole basis ; see [12] for earlier use of this device in the same setting, and, e.g., [7] in the matrix function context.

Throughout the paper, Greek bold letters () will denote matrices, while roman capital letters () larger ones. In particular will denote the th block of columns of the identity matrix . Scalar quantities will be denoted by Greek letters ().

Here is a synopsis of the paper. In Section 2 the basic tools of projection methods for solving (1) are recalled. In Section 3.1 we present a cheap residual norm computation whose implementation is discussed in Section 3.2. The two-pass strategy for is examined in Section 3.3. In Section 4 we extend the residual computation to . Section 5 discusses the generalization of this procedure to the case of the Sylvester equation in (1). In particular, Section 5.1 analyzes the case when both coefficient matrices are large, while Section 5.2 discusses problems where one of them has small dimensions. Numerical examples illustrating the effectiveness of our strategy are reported in Section 6, while our conclusions are given in Section 7.

## 2 Galerkin projection methods

Consider a subspace spanned by the orthonormal columns of the matrix and seek an approximate solution to (1) of the form with symmetric and positive semidefinite, and residual matrix . With the matrix inner product

 ⟨Q,P⟩F:=trace(PTQ),Q,P∈Rn1×n2,

the matrix can be determined by imposing an orthogonality (Galerkin) condition on the residual with respect to this inner product,

 Rm⊥Km⇔VTmRmVm=0. (2)

Substituting into (2), we obtain , that is

 (VTmAVm)YmVTmVm+VTmVmYm(VTmAVm)+VTmCCTVm=0. (3)

We assume , that is for some nonsingular . Since has orthonormal columns, and equation (3) can be written as

 TmYm+YmTm+E1γγTET1=0, (4)

where is symmetric and negative definite. The orthogonalization procedure employed in building determines the sparsity pattern of . In particular, for , the block Lanczos process produces a block tridiagonal matrix with blocks of size ,

 Tm=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝τ11τ12τ21τ22τ23⋱⋱⋱⋱⋱τm−1,mτm,m−1τm,m⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

As long as is of moderate size, methods based on the Schur decomposition of the coefficient matrix can be employed to solve equation (4), see, e.g., [2], [9].

The last columns (or rows) of the solution matrix are employed to compute the residual norm. In particular, letting , it was shown in [11] that the norm of the residual in (2) satisfies

 ∥Rm∥F=√2∥YmT––TmEm+1∥F=√2∥YmEmτTm+1,m∥F. (5)

The matrix is determined by solving (4), and it is again symmetric and positive semidefinite. At convergence, the backward transformation is never explicitly computed or stored. Instead, we factorize as

 Ym=ˆYˆYT,ˆY∈Rsm×sm, (6)

from which a low-rank factor of is obtained as , . The matrix may be numerically rank deficient, and this can be exploited to further decrease the rank of . We write the eigendecomposition of , (with eigenvalues ordered non-increasingly) and discard only the eigenvalues below a certain tolerance, that is , with (in all our experiments we used ). Therefore, we define again , with , ; in this way, . Hence, we set . We notice that a significant rank reduction in is an indication that all relevant information for generating is actually contained in a subspace that is much smaller than . In other words, the generated Krylov subspace is not efficient in capturing the solution information and a much smaller space could have been generated to obtain an approximate solution of comparable accuracy.

Algorithm 1 describes the generic Galerkin procedure to determine and as grows, see, e.g., [20]. Methods thus differ for the choice of the approximation space. If the block Krylov space is chosen, the block Lanczos method can be employed in line of Algorithm 1. In exact arithmetic,

 Vmτm+1,m=AVm−1−Vm−1τm,m−Vm−2τm−1,m. (7)

Algorithm 2 describes this process at iteration , with , where the orthogonalization coefficients ’s are computed by the modified block Gram-Schmidt procedure (MGS), see, e.g., [17]; to ensure local orthogonality in finite precision arithmetic, MGS is performed twice (beside each command is the leading computational cost of the operation). To simplify the presentation, we assume throughout that the generated basis is full rank. Deflation could be implemented as it is customary in block methods whenever rank deficiency is detected.

We emphasize that only the last terms of the basis must be stored, and the computational cost of Algorithm 2 is fixed with respect to . In particular, at each iteration , Algorithm 2 costs flops.

As the approximation space expands, the principal costs of Algorithm 1 are steps 4 and 6.1. In particular, the computation of the whole matrix requires full matrix-matrix operations and a Schur decomposition of the coefficient matrix , whose costs are flops. Clearly, step 6.1 becomes comparable with step 4 in cost for , for instance if convergence is slow, so that .

Step 9 of Algorithm 1 shows that at convergence, the whole basis must be saved to return the factor . This represents a major shortcoming when convergence is slow, since may require large memory allocations.

## 3 Standard Krylov subspace

For the block space , we devise a new residual norm expression and discuss the two-pass strategy.

### 3.1 Computing the residual norm without the whole solution Ym

The solution of the projected problem (4) requires the Schur decomposition of . For real symmetric matrices, the Schur decomposition amounts to the eigendecomposition , , and the symmetric block tridiagonal structure of can be exploited so as to use only flops; see section 3.2 for further details. Equation (4) can thus be written as

 Λm˜Y+˜YΛm+QTmE1γγTET1Qm=0,where˜Y:=QTmYmQm. (8)

Since is diagonal, the entries of can be computed by substitution [20, Section 4], so that

 Ym=Qm˜YQTm=−Qm(eTiQTmE1γγTET1Qmejλi+λj)ijQTm, (9)

where denotes the th vector of the canonical basis of . It turns out that only the quantities within parentheses in (9) are needed for the residual norm computation, thus avoiding the cost of recovering .

{proposition}

Let denote the eigendecomposition of . Then

 ∥Rm∥2F = 2(∥eT1SmD−11Wm∥22+…+∥eTsmSmD−1smWm∥22), (10)

where , and for all .

{proof}

Exploiting (5) and the representation formula (9) we have

 ∥Rm∥2F=2∥YmEmτTm+1,m∥2F=2∥∥ ∥∥(eTiQTmE1γγTET1Qmejλi+λj)ijQTmEmτTm+1,m∥∥ ∥∥2F=2s∑k=1∥∥ ∥∥(eTiSmejλi+λj)ijWmek∥∥ ∥∥22. (11)

For all , we can write

 ∥∥ ∥∥(eTiSmejλi+λj)ijWmek∥∥ ∥∥22=(sm∑j=1eT1Smejλ1+λjeTjWmek)2+…+(sm∑j=1eTsmSmejλsm+λjeTjWmek)2=(eT1SmD−11Wmek)2+…+(eTsmSmD−1smWmek)2. (12)

Plugging (12) into (11) we have

 ∥Rm∥2F=2s∑k=1sm∑i=1(eTiSmD−1iWmek)2=2sm∑i=1s∑k=1(eTiSmD−1iWmek)2=2sm∑i=1∥∥eTiSmD−1iWm∥∥22.

### 3.2 The algorithm for the residual norm computation

Algorithm 3 summarizes the procedure that takes advantage of Proposition 3.1. Computing the residual norm by (11) has a leading cost of flops for standard Krylov (with ). This should be compared with the original procedure in steps and of Algorithm 1, whose cost is flops, with a large constant. Proposition 3.1 also shows that only the first and last components of the eigenvectors of are necessary in the residual norm evaluation and the computation of the complete eigendecomposition may be avoided. To this end, the matrix can be tridiagonalized, , explicitly computing only the first and last rows of the transformation matrix , namely and . The eigendecomposition is computed exploiting the tridiagonal structure of . The matrices and needed in (10) are then computed as , .

Once the stopping criterion in step 6.3 of Algorithm 1 is satisfied, the factor can be finally computed. Once again, this can be performed without explicitly computing , which requires the expensive computation . Indeed, the truncation strategy discussed around (6) can be applied to by computing the matrix , so that . This factorization further reduces the overall computational cost, since only flops are required to compute , with no loss of information at the prescribed accuracy. The solution factor is then computed as .

To make fair comparisons with state-of-the-art algorithms that employ LAPACK and SLICOT subroutines (see Section 6 for more details), we used a C-compiled mex-code cTri to implement Algorithm 3, making use of LAPACK and BLAS subroutines. In particular, the eigendecomposition is performed as follows. The block tridiagonal matrix is tridiagonalized, , by the LAPACK subroutine dsbtrd that exploits its banded structure. The transformation matrix is represented as a product of elementary reflectors and only its first and last rows, , , are actually computed. The LAPACK subroutine dstevr is employed to compute the eigendecomposition of the tridiagonal matrix . This routine applies Dhillon’s MRRR method [5] whose main advantage is the computation of numerically orthogonal eigenvectors without an explicit orthogonalization procedure. This feature limits to flops the computation of ; see [5] for more details. Since the residual norm computation (10) requires the first and last rows of the eigenvectors matrix , we compute only those components, that is and , avoiding the expensive matrix-matrix product .

### 3.3 A “two-pass” strategy

While the block Lanczos method requires the storage of only basis vectors, the whole is needed to compute the low-rank factor at convergence (step 9 of Algorithm 1). Since

 (13)

we suggest not to store during the iterative process but to perform, at convergence, a second Lanczos pass computing and adding the rank- term in (13) at the th step, in an incremental fashion. We point out that the orthonormalization coefficients are already available in the matrix , therefore is simply computed by repeating the three-term recurrence (7), which costs flops plus the multiplication by , making the second Lanczos pass cheaper than the first one.

## 4 Extended Krylov subspace

Rational Krylov subspaces have shown to provide dramatic performance improvements over classical polynomial Krylov subspaces, because they build spectral information earlier, thus generating a much smaller space dimension to reach the desired accuracy. The price to pay is that each iteration is more computationally involved, as it requires solves with the coefficient matrices. The overall CPU time performance thus depends on the data sparsity of the given problem; we refer the reader to [20] for a thorough discussion.

In this section we show that the enhanced procedure for the residual norm computation can be applied to a particular rational Krylov based strategy, the Extended Krylov subspace method, since also this algorithm relies on a block tridiagonal reduced matrix when data is symmetric. Different strategies for building the basis of the extended Krylov subspace can be found in the literature, see, e.g., [10],[13],[19]. An intuitive key fact is that the subspace expands in the directions of and at the same time. In the block case, a natural implementation thus generates two new blocks of vectors at the time, one in each of the two directions. Starting with , the next iterations generate the blocks by multiplication by and solve with , respectively, and then setting . As a consequence, the block Lanczos procedure described in Algorithm 2 can be employed with (with columns). The orthogonalization process determines the coefficients of the symmetric block tridiagonal matrix with blocks of size ,

 Hm=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ϑ11ϑ12ϑ21ϑ22ϑ23⋱⋱⋱⋱⋱ϑm−1,mϑm,m−1ϑm,m⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠∈R2sm×2sm,

such that . The coefficients ’s correspond to the ’s in Algorithm 2, however as opposed to the standard Lanczos procedure, . Nonetheless, a recurrence can be derived to compute the columns of from those of during the iterations; see [19, Proposition 3.2]. The computed is block tridiagonal, with blocks of size , and this structure allows us to use the same approach followed for the block standard Krylov method as relation (5) still holds. Algorithm 3 can thus be adopted to compute the residual norm also in the extended Krylov approach with . Moreover, it is shown in [19] that the off-diagonal blocks of have a zero lower block, that is

 τi,i−1=[¯¯¯τi,i−10],¯¯¯τi,i−1∈Rs×2s,i=1,…,m.

This observation can be exploited in the computation of the residual norm as

 ∥Rm∥F=√2∥YmEmτTm+1,m∥F=√2∥YmEm¯¯¯τTm+1,m∥F,

and can be passed as an input argument to cTri instead of the whole .

The extended Krylov subspace dimension grows faster than the standard one as it is augmented by vectors per iteration. In general, this does not create severe storage difficulties as the extended Krylov approach exhibits faster convergence than standard Krylov in terms of number of iterations. However, for hard problems the space may still become too large to be stored, especially for large . In this case, a two-pass-like strategy may be appealing. To avoid the occurrence of new system solves with , however, it may be wise to still store the second blocks, , , and only save half memory allocations, those corresponding to the matrices , .

Finally, we remark that if we were to use more general rational Krylov subspaces, which use rational functions other than and to generate the space [20], the projected matrix would lose the convenient block tridiagonal structure, so that the new strategy would not be applicable.

## 5 The case of the Sylvester equation

The strategy presented for the symmetric Lyapunov equation (1) can be extended to the Sylvester equation

 AX+XB+C1CT2=0,A∈Rn1×n1,B∈Rn2×n2,C1∈Rn1×s,C2∈Rn2×s, (14)

where the coefficient matrices are both symmetric and negative definite while are tall, that is .

### 5.1 Large A and large B

We consider the case when both and are large and sparse matrices. If their eigenvalue distributions satisfy certain hypotheses, the singular values of the nonsymmetric solution to (14) exhibit a fast decay, and a low-rank approximation to can be sought; see, e.g., [18, Th. 2.1.1], [20, Section 4.4].

Projection methods seek an approximate solution to (14) of the form where the orthonormal columns of and span suitable subspaces and respectively111The space dimensions of and are not necessarily equal, we limit our discussion to the same dimension for simplicity of exposition.. The construction of two approximation spaces is thus requested and, for the sake of simplicity, we limit our discussion to the standard Krylov method, that is and , with obvious generalization to the extended Krylov subspace. As in the Lyapunov case, is computed by imposing a Galerkin condition on the residual matrix , that is

 VTmRmUm=0. (15)

We assume , for some nonsingular , and a similar discussion to the one presented in Section 2 shows that condition (15) is equivalent to solving the reduced Sylvester problem

 TmYm+YmJm+E1γ1γT2ET1=0, (16)

where , . Similarly to the Lyapunov case, computing the eigendecompositions , , and , , the solution to (16) can be written as

 Ym=Qm˜YPTm=−Qm(eTiQTmE1γ1γT2ET1Pmejυi+λj)ijPTm. (17)

The last rows and columns of are employed in the residual norm calculation. Indeed, letting and , it can be shown that

 ∥Rm∥2F = ∥T––mYm∥2F+∥YmJ––Tm∥2F=∥τm+1,mETmYm∥2F+∥YmEmιTm+1,m∥2F, (18)

where and , see, e.g., [20],[3].

The same arguments of Section 3.1 can be applied to the factors in (18) leading to Algorithm 4 for the computation of the residual norm without explicitly assembling the matrix . The eigendecompositions in step 1 are not fully computed. In particular, only the spectrum and the first and last components of the eigenvectors of and are explicitly computed following the strategy presented in Section 3.2.

At convergence, the matrix can be computed by (17). Also in the Sylvester problem the matrix may be numerically singular. In this case, factors , , such that can be computed via the truncated singular value decomposition of the nonsymmetric matrix . The low-rank factors of , , are then computed as and .

If equation (14) is solved by the standard Krylov method, the two-pass strategy presented in Section 3.3 can be easily adapted to the Sylvester case. Indeed, denoting by and , the low-rank factors and can be written as

 Z1=Vm(QmˆY1)=m∑i=1ViETi(QmˆY1),Z2=Um(PmˆY2)=m∑i=1UiETi(PmˆY2).

As in the Lyapunov case, the factors , can be computed in a second Lanczos pass since the terms and do not require the whole basis to be available. Therefore, for the Sylvester problem (14), the “two-pass” strategy allows us to store only basis vectors, vectors for each of the two bases.

### 5.2 Large A and small B

In some applications, such as the solution of eigenvalues problems [23] or boundary value problems with separable coefficients [22], the matrices and in (14) could have very different dimensions. In particular, one of them, for instance, , could be of moderate size, that is . In this case, the projection method presented in Section 5.1 can be simplified. Indeed, a reduction of the matrix becomes unnecessary, so that a numerical solution to (14) of the form is sought, where the columns of span , as before. The Galerkin condition on the residual matrix thus becomes