Preconditioned Galerkin and minimal residual methods for solving Sylvester equations

Applied Mathematics and Computation 181 (2006) 1208–1214 www.elsevier.com/locate/amc

Preconditioned Galerkin and minimal residual methods for solving Sylvester equations A. Kaabi *, F. Toutounian, A. Kerayechian Department of Mathematics, Ferdowsi University of Mashhad, Mashhad 54646446, Iran

Abstract This paper presents preconditioned Galerkin and minimal residual algorithms for the solution of Sylvester equations AX  XB = C. Given two good preconditioner matrices M and N for matrices A and B, respectively, we solve the Sylvester equations MAXN  MXBN = MCN. The algorithms use the Arnoldi process to generate orthonormal bases of certain Krylov subspaces and simultaneously reduce the order of Sylvester equations. Numerical experiments show that the solution of Sylvester equations can be obtained with high accuracy by using the preconditioned versions of Galerkin and minimal residual algorithms and this versions are more robust and more efficient than those without preconditioning.  2006 Elsevier Inc. All rights reserved. Keywords: Sylvester matrix equations; Preconditioning; Galerkin method; Minimal residual method; Krylov subspace

1. Introduction Matrix Sylvester equations are very important in control theory and many other branches of engineering [5,9,11,12,15,20,21]. In this paper, we focus on numerical solution of the Sylvester equations AX  XB ¼ C; nn

ð1Þ mm

nm

nm

where A 2 R , B 2 R , C 2 R , and X 2 R is the solution matrix sought. Without loss of the generality, throughout this paper, we suppose that m = n. The necessary and sufficient condition for (1) to have a unique solution is that SðAÞ \ SðBÞ ¼ ;; where S(A) and S(B) are the spectrums of A and B, respectively [8,10,14–17,21]. The Bartels and Stewart methods provided the first numerically stable way to systematically solve the Sylvester equation (1). When the matrices A and B are dense, the Bartels–Stewart [1] and Golub–Nash–Van Loan algorithms [10] are attractive. When the matrices A and B are large and sparse, iterative solution of (1) by the alternating-direction-implicit

*

Corresponding author. E-mail address: [email protected] (A. Kaabi).

0096-3003/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.02.021

A. Kaabi et al. / Applied Mathematics and Computation 181 (2006) 1208–1214

1209

(ADI) method may be attractive; see [3,4,6,7,18–20,22,23]. The Krylov-subspaces Galerkin and minimal residual algorithms for the solution of Sylvester equation (1), are presented by Hu and Reichel [13]. As we know, the rate of convergence of the iterative methods for solving a linear system is strongly depended on the spectral properties of coefficient matrix of the system. Hence, iterative methods usually involve a second matrix that transforms the coefficient matrix into one with a more favorable spectrum. The transformation matrix is called a preconditioner. Let M and N be two approximate inverse preconditioners of A and B, respectively, such that MA  I and NB  I, we can then convert Eq. (1) to the equation MAXN  MXBN ¼ MCN : This equation can be written as a linear system of equations ðN T  MA  N T BT  MÞx ¼ ðN T  MÞc;

ð2Þ n2

n2

where  denote the kronecker product, and x 2 R and c 2 R are vectors whose entries are the elements of the matrices X and C, respectively. More precisely, eTj Xek ¼ eTjþnðk1Þ x;

1 6 j; k 6 n;

and the same formula also relates the elements of C to those of c. Throughout this paper ej denotes the jth axis vector of suitable dimension. The main propose of this paper is to present the preconditioned versions of Galerkin and minimal residual algorithms for solving the Sylvester equation (1). Given an initial approximate solution x0, we replace the Krylov subspace Km(NT  MA  NTBT  M, r0), where r0 = (NT  M)c  (NT  MA  NTBT  M)x0, with a subspace of the form Km(NTBT, g)  Km(MA, f) for certain vectors f ; g 2 Rn . This reduces the order of the preconditioned Sylvester equations, the computer storage and arthumetic work required. The smaller preconditioned Sylvester equations obtained can be solved with a direct method such that the Bartels–Stwart algorithm [1]. This paper is organized as follows. In Sections 2 we describe a preconditioned Galerkin method for solution (2), or equivalenty (1). A preconditioned minimal residual algorithm is presented in Section 3. These algorithms will be derived from FOM and GMRES algorithms in the same way Galerkin and minimal residual algorithms are derived in [13]. In Section 4 some numerical examples are tested. Finally, Section 5 summarizes the main conclusion of this paper. 2. A preconditioned Galerkin algorithm for Sylvester equations In this section we derived a preconditioned version of the Galerkin method described in [13] for solving Sylvester equation (1). 2 Let M and N be good preconditioner matrices for matrices A and B, respectively, and x0 2 Rn be an initial approximate solution of linear system (2), then the residual vector is given by r0 ¼ ðN T  MÞc  ðN T  MA  N T BT  MÞx0 : Let, we define matrix PM,N and vector cM,N as follows: P M;N ¼ N T  MA  N T BT  M;

cM;N ¼ ðN T  MÞc:

The full orthogonalization method seeks an approximate solution xm from the affine subspace x0 + Km(PM,N, r0) of dimension m by imposing the Galerkin condition cM;N  P M;N xm ? K m ðP M;N ; r0 Þ: This method uses the Arnoldi process to compute an orthonormal basis of Krylov-subspace Km(PM,N, r0). The matrix PM,N in application can be very large, so, it can be difficult to store a basis of the Krylov-subspace Km(PM,N, r0) in fast computer memory for dimension m that yields an acceptable rate of convergence. For remedy this problem we replace the Krylov-subspace Km(PM,N, r0) with a subspace of the form Km(NTBT, g)  Km(MA, f) for certain vectors f ; g 2 Rn . Strategies for selecting f and g have been discussed in [12]. Let for the moment f and g be arbitrary but fixed in Rn , and use modified Gram Schmidt algorithm mþ1 mþ1 to generate orthonormal bases fvj gj¼1 and fwj gj¼1 of Km+1(MA, f) and Km+1(NTBT, g), respectively. By

1210

A. Kaabi et al. / Applied Mathematics and Computation 181 (2006) 1208–1214

introducing the matrices Vj = [v1, v2, . . . , vj] and Wj = [w1, w2, . . . , wj], for j 2 {m, m + 1}, modified Gram Schmidt algorithm yields the entries of the upper Hessenberg matrices HMA and HNB and Hessenberg-like e MA and H e NB that satisfy, matrices H H MA ¼ V Tm MAV m ;

H B ¼ W Tm N T BT W m ;

and e MA ; MAV m ¼ V mþ1 H

e NB : N T BT W m ¼ W mþ1 H

ð3Þ

The columns of the matrix Wm  Vm form an orthonormal basis of Km(NTBT, g)  Km(MA, f). We wish to determine a correction z0 2 Km(NTBT, g)  Km(MA, f) of x0 and obtain a new approximate solution x1 = x0 + z0 of (2), such that the residual vector r1 ¼ cM;N  P M;N x1 T

ð4Þ T

is orthogonal to Km(N B , g)  Km(MA, f). The correction z0 can be written as z0 ¼ ðW m  V m Þy 0

ð5Þ

m2

T

for y 0 2 R . Then from the orthogonality relation, (Wm  Vm) r1 = 0, we have  T T  W m N W m  H MA  H NB  V Tm MV m y 0 ¼ ~r0 ;

ð6Þ

where ~r0 ¼ ðW m  V m ÞT r0 : The system (6) is a small Sylvester equations, which can be solve with a direct method such the Bartels–Stewart algorithm [1]. We outline how this algorithm proceeds when applied to (6). First determine the Schur factorizations H MA ¼ QMA RMA QMA ;

H NB ¼ QNB RNB QNB ;

where QMA, QNB 2 Cm·m are unitary and RMA, RNB 2 Cm·m are upper triangular. Here * denotes transposition and complex conjugation. Then (6) is equivalent with    ð7Þ QNB W Tm N T W m QNB  RMA  RNB  QMA V Tm MV m QMA y 0 ¼ r0 ; where y 0 ¼ ðQNB  QMA Þ y 0

ð8Þ

and r0 ¼ ðQNB  QMA Þ~r0 : Having computed y0, we can determine z0 by using (5). We now can compute x1 = x0 + z0 and the corresponding residual vector r1. The computation of r1 from its definition can be expensive unless the matrices MA and NB are very sparse. It is often cheaper to evaluate the expression e MA Þy 0 þ ðW mþ1 H e NB  MV m Þy 0 : r1 ¼ r0  ðN T W m  V mþ1 H

ð9Þ

If kr1k2 is not sufficiently small, then we let x0 = x1 and r0 = r1, and apply the preconditioned Galerkin method again with new vectors f ; g 2 Rn . This is the restarted preconditioned Galerkin algorithm. In the restarted algorithm the residual vector r1 should be evaluated periodically by (4) because of the possible loss of accuracy when (9) is used. Algorithm 1. Restarted preconditioned Galerkin method (PGal) for the solution of (2) 1. Choose an initial approximate solution x0, tolerance  > 0 and integer m. 2. Compute r0 = cM,N  PM,Nx0. 3. If kr0k2 <  then exit.

A. Kaabi et al. / Applied Mathematics and Computation 181 (2006) 1208–1214

1211

4. Choose f ; g 2 Rn , as in [13]. Apply Arnoldi process to compute orthonormal bases of Km+1(MA, f) and e MA , H e NB , Vm, Vm+1, Wm, and Wm+1. Km+1(NTBT, g). We obtain matrices H T 5. Compute ~r0 ¼ ðW m  V m Þ r0 . e MA , R e NB , R e MA from H e MA and QNB, Q e NB from H e NB . 6. Determine QMA, Q  e e 7. Compute r0 ¼ ð Q NB  Q MA Þ ~r0 . 8. Solve the small Sylvester equation (7). 9. Compute y 0 ¼ ðQNB  QMA Þy 0 . 10. Compute z0 = (Wm  Vm)y0. 11. Compute x0 = x0 + z0. e MA Þy 0 þ ðW mþ1 H e NB  MV m Þy 0 . 12. Compute r0 ¼ r0  ðN T W m  V mþ1 H 13. Go to 3. 3. A preconditioned minimal-residual algorithm for Sylvester equations In this section we derived a preconditioned version of the minimal residual method described in [13] for solving Sylvester equation (1). Let x0 be a given initial approximate solution of (2), and the residual vector is r0 = cM,N  PM,Nx0. We would like to determine a new approximate solution xm = x0 + z0 with z0 2 Km(NTBT, g)  Km(MA, f) chosen so that the norm of the residual vector (4) associate with x1 is minimal. An algorithm for the case when r0 2 Km+1(NTBT, g)  Km+1(MA, f) can be derived as follows. Let ~r0 ¼ ðW mþ1  V mþ1 ÞT r0 : Then r0 ¼ ðW mþ1  V mþ1 Þ~r0 :

ð10Þ

The correction z0 can be expressed by (5). Substitution (3), (5), and (10) into r1 ¼ r0  ðN T  MA  N T BT  MÞz0 yields

h i e MA  H e NB  V T MV m Þy 0 : r1 ¼ ðW mþ1  V mþ1 Þ ~r0  ðW Tmþ1 N T W m  H mþ1

Therefore, minimizing kr1k2 over all z0 2 Km(NTBT, g)  Km(MA, f) is equivalent with the task   e MA  H e NB  V T MV m Þy 0  min ~r0  ðW Tmþ1 N T W m  H mþ1 2

ð11Þ

m2

over all y 0 2 R . Recalling the Schur factorizations of HMA and HNB, we define unitary (m + 1) · (m + 1) matrices with QMA and QNB as leading principal submatrices     QMA 0 QNB 0 e e Q MA  ; Q NB  ; 0 1 0 1 and introduce the matrices e H e MA QMA ; e MA ¼ Q R MA

e H e e NB ¼ Q R NB NB QNB :

e MA is RMA, and (m + 1)st row of R e MA is given by Then the leading m · m principal submatrix of R T e MA ej ¼ eT H e eTmþ1 R mþ1 MA em em QMA ej ;

ð12Þ

1 6 j 6 m:

e NB are given by a e NB is RNB, and the elements of row m + 1 of R The leading m · m principal submatrix of R e NB  Q e MA Þ~r0 and y 0 be defined by (8). Then (11) is equivalent with formula similar to (12). Let r0 ¼ ð Q    e  W T N T W m QNB  R e  V T MV m QMA Þy 0  e MA  R e NB  Q ð13Þ min r0  ð Q  NB mþ1 MA mþ1 2

1212

A. Kaabi et al. / Applied Mathematics and Computation 181 (2006) 1208–1214 2

over all y 0 2 Rm . We can determine y 0 by solving this least squares problem and then computing y0, z0, x1, and r1. The computation of r1 from its definition can be expensive unless the matrices MA and NB are very sparse. It is often cheaper to evaluate the expression e MA Þy 0 þ ðW mþ1 H e NB  MV m Þy 0 : r1 ¼ r0  ðN T W m  V mþ1 H If kr1k2 is not sufficiently small, then we let x0 = x1 and r0 = r1, and apply the preconditioned minimal residual method again with new vectors f ; g 2 Rn . This is the restarted preconditioned minimal residual algorithm. Algorithm 2. Restarted preconditioned minimal residual method (PMR) for the solution of (2) 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.

Choose an initial approximate solution x0, tolerance  > 0 and integer m. Compute r0 = cM,N  PM,Nx0. If kr0k2 <  then exit. Choose f ; g 2 Rn , as in [13]. Apply Arnoldi process to compute orthonormal bases of Km+1(MA, f) and e MA , H e NB , Vm, Vm+1, Wm, and Wm+1. Km+1(NTBT, g). We obtain matrices H T Compute ~r0 ¼ ðW mþ1  V mþ1 Þ r0 . e MA , R e NB , R e MA from H e MA and QNB, Q e NB from H e NB . Determine QMA, Q  e NB  Q e MA Þ ~r0 . Compute r0 ¼ ð Q Solve the least-squares problem (13). Denote the solution by y 0 . Compute y 0 ¼ ðQNB  QMA Þy 0 . Compute z0 = (Wm  Vm)y0. Compute x0 = x0 + z0. e MA Þy 0 þ ðW mþ1 H e NB  MV m Þy 0 . Compute r0 ¼ r0  ðN T W m  V mþ1 H Go to 3.

4. Numerical experiments In this section, we compared the performance of the Galerkin (Gal), the preconditioned Galerkin (PGal), the minimal residual (MR) and the preconditioned minimal residual (PMR) algorithms. For all the examples, we used the stopping criterion krk k2 6 106 ; kr0 k2 and the maximum number of iterations allowed set to 1000. The number of iterations are presented in Tables 1–4. For all the experiments, the initial guess was X0 = 0. The right-hand-side matrix C is chosen so that Xi, j = 1, 1 6 i 6 n, 1 6 j 6 n solves Eq. (1). For producing a sparse inverse preconditioner, we used the biconjugation algorithm [2]. For preserving sparsity, we used the following dropping strategy [10]. • Dropping entries of mj, columns of M, (nj, columns of N) (i.e., replaced by zero) which are less than the relative tolerance sj(cj) obtained by multiplying s by the original 2-norm of the jth column of A(B) corresponding to mj(nj). Table 1 m=0 n·n s

10 · 10 107

15 · 15 106

20 · 20 107

25 · 25 107

30 · 30 1010

MR PMR

75 5

190 8

340 10

460 7

800 10

Gal PGal

70 5

170 8

290 10

450 12

880 15

A. Kaabi et al. / Applied Mathematics and Computation 181 (2006) 1208–1214

1213

Table 2 m = 10 n·n s

10 · 10 103

15 · 15 103

20 · 20 103

25 · 25 103

30 · 30 103

MR PMR

90 40

150 57

220 65

350 105

500 140

Gal PGal

85 35

170 63

230 70

310 92

430 162

Table 3 m=4 n·n s

10 · 10 103

15 · 15 103

20 · 20 103

25 · 25 104

30 · 30 103

MRM PMRM

10 5

40 15

60 20

95 35

130 30

GM PGM

12 3

35 4

60 21

80 28

105 30

Table 4 m=2 A B n·n s

fidap005 hilbert 10 · 10 1011

pores1 hilbert 15 · 15 1013

fidapm05 hilbert 20 · 20 1015

bcsstk01 bcsstm05 25 · 25 1020

bcsstk01 hilbert 30 · 30 1013

MRM PMRM

–a 6

270 3

–a 3

600 2

280 9

GM PGM

–a 2

220 3

–a 2

700 3

250 9

a

No solution has been obtained after 1000 iterations.

As a first test problem, we consider the ellipitic operator LðuÞ ¼ Du þ 2mux þ 2muy with Dirichlet boundary conditions. The operator was discretized using central finite differences on [0, 1] · [0, 1]. The mesh sizes in both the x- and y-directions is h = 1/(n + 1). This yields two tridiagonal matrices A and B given by A ¼ tridiagð1  mh; 2; 1 þ mhÞ and

B ¼ tridiagð1  mk; 2; 1 þ mkÞ:

Then we have to solve a Sylvester matrix equation AX  XB = C. The four algorithms were run in a restarted mode m = 2. The results obtained are reported in Tables 1 and 2 with different values of n, m, and s. These tables show that the number of iterations PMR and PGal algorithms are very smaller than those of MR and Gal algorithms. As a second test problem, we consider the coefficient matrix G = A  I + I  BT of the linear system Gx = b arising in solving the Poisson equation on the rectangular domain by using nine point disctretization formula. In this case the n · n matrices in the Sylvester equation (1) are, respectively, A ¼ pentadiagð1; 16; 30; 16; 1Þ and

B ¼ pentadiagð1; 16; 30; 16; 1Þ:

The results obtained are presented in Table 3. These results show that PMR and PGal algorithms performed better than MR and Gal algorithms, respectively.

1214

A. Kaabi et al. / Applied Mathematics and Computation 181 (2006) 1208–1214

In the third set of experiments, we used the matrices A and B of the Sylvester equation (1) from Harwell– Boing collection. The results obtained are presented in Table 4. We mention that the results of column 5 of Table 4 have been obtained with N = I, because the matrix bcsstm05 is singular and we can not obtain an approximate inverse preconditioner. As shown, all Sylvester equations were solved in a relatively small number of iterations by PMR and PGal algorithms. 5. Conclusion We have presented the preconditioned versions of Galerkin and minimal residual algorithms for solving the Sylvester equations. The algorithms use the Arnoldi process to generate orthonormal bases of certain Krylov subspaces and reduce the order of the Sylvester equations, the computer storage and arithmetic work required. The experiments presented in this paper showed that the solution of Sylvester equations can be obtained with high accuracy by using the preconditioned versions of Galerkin and minimal residual algorithms and this versions are more robust and more efficient than those without preconditioning. References [1] R.H. Bartels, G.W. Stewart, Algorithm 432: Solution of matrix equation AX + XB = C, Circ. Syst. Signal Process. 13 (1994) 820–828. [2] M. Benzi, M. Tuma, A sparse approximate inverse preconditioner for nonsymmetric linear systems, SIAM J., Sci. Comput. 19 (3) (1998) 968–994. [3] R. Bhatia, C. Davis, A. Mcintosh, Pertubation of spectoral subspaces and solution of linear operator equations, Linear Algebra Appl. 52/53 (1983) 45–67. [4] G. Birkhoff, R.S. Varga, D. Young, Alternating direction implicit methods, Advances in Computing, vol. 3, Academic, New York, 1962, pp. 189–273. [5] E. Chow, Y. Saad, Approximate inverse preconditioner via sparse-sparse iterations, SIAM J. Sci. Comput. 19 (1998) 995–1023. [6] B.N. Datta, K. Datta, Theoretical and computational aspects of some linear algebra problems in control theory, in: C.I. Byrnes, A. Lindquist (Eds.), Computational and Combinatorial Methods in Systems Theory, Elsevier, Amsterdam, 1986, pp. 201–212. [7] N.S. Ellner, E.L. Wachspress, Alternating direction implicit iteration for systems with complex spectra, SIAM J. Numer. Anal. 28 (1991) 859–870. [8] D.J. Evans, E. Galligani, A parallel additive preconditioner for conjugate gradient method for AX + XB = C, Parallel Comput. 20 (1994) 1055–1064. [9] E. Gallopoulos, Y. Saad, On the parallel solution of parabolic equations, in: R. De Groot (Ed.), Proceedings of the International Conference on Supercomputing 1989, Heraklion, Crete, June 5–9, 1989, ACM Press, New York, 1989. [10] G.H. Goulb, S. Nash, G.F. Van Loan, A Hessenberg–Schur method for the problem AX  XB = C, IEEE Trans. Automat. Control AC-24 (1979) 909–913. [11] G.H. Goulb, G.F. Van Loan, Matrix Computations, second ed., Johns Hoplins U.P., Baltimore, 1989. [12] M. Hochbruck, C. Lubich, On the Krylov subspace approximation to the matrix exponential operator, SIAM. J. Numer. Anal. 34 (1997) 1911–1925. [13] D.Y. Hu, L. Reichel, Krylov-subspace methods for the Sylvester equation, Linear Algebra Appl. 172 (1992) 283–313. [14] K. Jbilou, On the Krylov subspace methods for solving the matrix equation AX  XB = C, Technical report LMPA (71), Universite du Littoral, Calais, France, 1998. [15] A. Kaebi, A. Kerayehchian, F. Toutounian, Approximate inverse preconditioner by computing approximate solution of Sylvester equations, Appl. Math. Comput. 170 (2005) 1067–1076. [16] M. Rosenblum, On the operator equation BX  XA = Q, Duke Math. J. 23 (1956) 263–269. [17] Y. Saad, Numerical solution of large Lyapanov equations, in: M.A. Kaashoek, J.H. van Schuppen, A.C.M. ran (Eds.), Signal Processing, Scattering, Operator Theory and Numerical Methods, Birkhauser, Boston, 1990, pp. 503–511. [18] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Press, New York, 1995. [19] Y. Saad, M.H. Schultz, GMRES, A generalized minimal residual algorithm for nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 7 (1986) 856–869. [20] G. Starke, R.S. Varga, A hybrid Arnoldi-Faber iterative method for nonsymmetric systems of linear equations, Numer. Math. 64 (1993) 213–240. [21] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1962. [22] E.L. Wachspress, Iterative solution of the Lyapanov matrix equation, Appl. Math. Lett. 1 (1988) 87–90. [23] E.L. Wachspress, in: D.R. Kincaid, L.J. Hayes (Eds.), The ADI Minimax Problem for Complex Spectra in Iterative Methods for Large Linear Systems, Academic, San Diego, 1990, pp. 251–271.