On inexact alternating direction implicit iteration for continuous Sylvester equations

In this paper, we study the alternating direction implicit (ADI) iteration for solving the continuous Sylvester equation AX + XB = C, where the coefficient matrices A and B are assumed to be positive semi‐definite matrices (not necessarily Hermitian), and at least one of them to be positive definite. We first analyze the convergence of the ADI iteration for solving such a class of Sylvester equations, then derive an upper bound for the contraction factor of this ADI iteration. To reduce its computational complexity, we further propose an inexact variant of the ADI iteration, which employs some Krylov subspace methods as its inner iteration processes at each step of the outer ADI iteration. The convergence is also analyzed in detail. The numerical experiments are given to illustrate the effectiveness of both ADI and inexact ADI iterations.


INTRODUCTION
Consider the iterative solution to the continuous Sylvester equation by the alternating direction implicit (ADI)-like iterations, where A ∈ C m×m , B ∈ C n×n , C ∈ C m×n are large sparse matrices. For definiteness, throughout this paper, both A and B in (1) are assumed to be positive semi-definite * and at least one of them to be positive definite.
It is known that the Sylvester equation ( Sylvester equations play important roles in numerous applications such as matrix eigen-decompositions, control theory, model reduction, numerical solution of matrix differential Riccati equations, image processing, and many more, see for example, References 3-8 and a large literature therein. The Sylvester equation (1) is mathematically equivalent to the larger linear system of the form where = I m ⊗ A + B T ⊗ I n with ⊗ denoting the standard Kronecker product symbol, x, c are two column-stacking vectors of the matrices X and C, respectively. It is useful to treat Equation (2) as a general linear system in theoretical analysis, but impractical in numerical solution to the continuous Sylvester equation (1), because Equation (2) is costly to solve and can be ill-conditioned.
The common approaches to solving (1) are the Bartels-Stewart 9 and the Hessenberg-Schur 7,10 methods, each of which needs to transform A and B into triangular or Hessenberg form by an orthogonal similarity transformation and then solving the resulting system directly by Gaussian elimination with partial pivoting. Those methods are usually referred to as the direct methods. The direct methods are mainly applicable for small and medium model problems but often too expensive to be practical for large sparse problems.
When A and B are large sparse matrices, an alternative for solving (1) is iterative methods. The Smith method 11 and the ADI iterative method 5,[12][13][14][15][16] are popular ones for solving large sparse Sylvester equations. Much more attention was paid to the case of the coefficient matrices A and B being either Hermitian positive definite matrices or M-matrices. However, little attention was focused on the case of the coefficient matrices A and B being non-Hermitian positive definite matrices. Recently, some authors applied the state-of-the-art iterative methods, such as the Hermitian and skew-Hermitian splitting (HSS) iteration, 17 the positive-definite and skew-Hermitian splitting (PSS) iteration 18 and the normal and skew-Hermitian splitting (NSS) iteration 19 for solving non-Hermitian positive definite linear systems † , to develop HSS, PSS, and NSS iteration solvers for (1). Of those methods, one also needs to solve two Sylvester equations with certain structural coefficient matrices at each inner iteration, those special structures allow the uses of the numerically stable Cholesky/Bunch-Parlett factorizations (as for the direct methods) or the short-term recurrence CG/GMRES (CGNE, CGNR) methods (as for the iteration methods). The resulting methods converge unconditionally and are efficient and robust numerically, see for example References 4,21,22. Nevertheless, the coefficient matrices may be dense (for instance, when the matrix A is an upper Hessenberg matrix, H and S in the HSS splittings are still very dense), see Reference 18. This motivates us to further study the ADI-like iterations.
In this paper, we revisit the ADI iteration for solving (1). We first analyze its convergence and then derive an upper bound for its contraction factor. To reduce the computational complexity, we further establish an inexact variant of the ADI (IADI for short) iteration. The convergence of the IADI method is also analyzed in detail. Numerical experiments show that those methods are efficient and robust solvers and have a better performance than the HSS iteration solver for (1), and the IADI iteration is usually superior to the ADI iteration in terms of computation efficiency. Moreover, the IADI iteration as well as the ADI iteration outperform the IHSS iteration in terms of both the number of iterations and the computation efficiency.
The organization of this paper is as follows. After analyze the convergence rate of ADI iterative method for solving (1) in next section, we then establish the IADI iteration to improve the computing efficiency of the ADI iteration in section 3. Numerical experiments are illustrated in section 4 to show the effectiveness and robustness of our methods.

THE ADI ITERATION
Let us begin with some basic notations. For convenience, throughout this paper, we denote by k ∈ C n 2 the column-stacking vector of the matrix K ∈ C n×n , and we denote by (K), (K), ||K|| 2 and ||K|| F the spectrum, spectral radius, the 2-norm, and the Frobenius norm of the matrix K ∈ C n×n , respectively. A matrix sequence {Y (k) } ∞ k=0 ⊆ C n×n is said to be convergent to a matrix Y ∈ C n×n if the corresponding column-stacking vector sequence of {y (k) } ∞ k=0 ⊆ C n 2 of {Y (k) } ∞ k=0 is convergent to the corresponding column-stacking vector y ∈ C n 2 of Y .
The classical ADI iterative method for solving (1) is as follows. † In fact, such kind of methods are analogues to the classical ADI iteration introduced by Peaceman and Rachford for solving partial differential equations, also see Reference 20.
In this case, ADI iteration (3) reduces to Smith iteration, [ℋ ( )] is bounded bŷ with ′ min and ′ max denoting the lower and the upper bounds of the real part of the eigenvalues of the matrices A and B, and ′′ min and ′′ max denoting the lower and the upper bounds of the absolute values of the imaginary part of the eigenvalues of the matrices A and B. The parameter * is chosen such that the above estimate can be minimized. This fact is precisely stated as the following theorem.

Lemma 3. The minimizer of̂( ) over all positive is attained at
and the corresponding minimum value is equal tô The proof is a verbatim of one of corollary 2.3 in Reference 17 and therefore omitted.
Remark 1. We remark that the condition − ′ min < Δ < ′ min is sufficient but not necessary. That is to say that even if − ′ min < Δ < ′ min does not hold, we cannot guarantee that the upper bound is less than 1, but it does not mean the spectral radius [ℋ ( , )] cannot be less than 1.
In fact, if the imaginary parts of the eigenvalues are zero, for example, when A and B are both Hermitian positive definite matrices, then condition − ′ min < Δ < ′ min always holds. When A and B are both non-Hermitian positive definite matrices, we cannot guarantee that − ′ min < Δ < ′ min . From (10), however, we observe that if one of the following cases: (1) ′′ max and ′′ max are small enough, holds, then we can find a Δ such that − ′ min < Δ < ′ min holds. This phenomenon appears in our numerical tests. Therefore, we have the following corollary. Corollary 1. Under the assumption of Theorem 1, the iterative sequence {X (k) } determined by (3) converges to the exact solution X * of (1), if ′′ max and ′′ max are small enough, || ′′ max | − | ′′ max || is sufficient small, or ′ max (and ′ max ) is much larger than ′′ max (and ′′ max ). Corollary 2. If = , i.e., Δ = 0, then the iteration (3) converges to the exact solution unconditionally for all > 0. (1), then (B) = (A). From (10), we therefore have Δ = 0. In this case we obtain the following result.

THE IADI ITERATION
The two half-steps at each step of the ADI iteration for solving continuous Sylvester equation (1) need to solve two matrix equations like where Y 1 and Y 2 are known. This may be very costly and impractical in actual implementations. To further improve the computational efficiency of the ADI iteration, we can solve the two subproblems in (13) inexactly by employing some state of the art iterative methods such as GMRES, which results in the basic framework of the IADI iterative method for solving (1).
Given an initial guess X (0) ∈ C n×n , for k = 1,2, … , until {X (k) } ∞ k=0 ⊆ C n×n satisfies the stopping criterion, solve X by employing an inner iteration (e.g., the GMRES) with X (k) as the initial guess; then solve X (k + 1) approximately from by employing an inner iteration (e.g., GMRES) with X as the initial guess, where , are given positive constants. To simplify numerical implementation and convergence analysis, based on the residual-updating form 28 we may rewrite the above IADI iteration as the following equivalent scheme.
The IADI Iteration. Given an initial guess X (0) ∈ C n×n , for k = 1,2, … , until {X (k) } ∞ k=0 ⊆ C n×n converges: 1. approximate the solution of is such that the residual and then compute where { k } and { k } are prescribed tolerances used to control the accuracies of the inner iterations.
In order to analyze the convergence of the above IADI iteration, we need to recall the following convergence theorem of iterative solution to a general linear system Ax = b by an inexact two-step splitting iterative method. Let ||v|| (for all v ∈ C n ) be a general vector norm and M be a nonsingular matrix. We define a new vector norm by |||v||| M = ||Mv|| (for all v ∈ C n ), then its induced matrix norm is |||K||| M = ||MKM −1 || (for all K ∈ C n×n ).

Lemma 4 (17, theorem 3.1). Let A ∈ C n×n and A
, and , then {x (k) } is of the form ).
Moreover, if x * ∈ C n is the exact solution of the linear system Ax = b, then we have In particular, if then the iterative sequence {x (k) } converges to x * ∈ C n , where max = max k k and max = max k k .
Now we can demonstrate the following convergence result concerning the above IADI iterative method.

Theorem 2.
Let {X (k) } ∞ k=0 ⊆ C n×n be an iteration sequence generated by the IADI iterative method and let X * ∈ C n×n be the exact solution of (1). Under the assumption of Theorem 1, we have that In particular, if then the iteration sequence {X (k) } ∞ k=0 converges to X * , where max = max k k and max = max k k .
Proof. 1,2) are two splittings of . Again, by making use of the kronecker product, we can rewrite the above-described IADI iteration in the following matrix-vector form: with r (k) = c − x (k) , where z (k) is the approximate solution such that the residual satisfies ||p (k) || 2 ≤ k ||r (k) || 2 , and is the approximate solution such that the residual By Lemma 4, we can easily obtain Now, taking |||y||| ℳ 2 = ||ℳ 2 y|| 2 , then we have Hence, we can equivalently rewrite the inequality (19) as Thus we complete the proof of this theorem. ▪ According to Theorem 2, we want to choose tolerances { k } and { k } so that the computational work of the IADI iterative method is minimized. In fact, the tolerances { k } and { k } are not required to approach zero as k increases in order to ensure the convergence of the IADI iteration but are required to approach zero in order to asymptotically recover the original convergence rate (cf. Theorem 1) of the ADI iteration. How to arrive at a tradeoff between the computational complexity and the convergence rate is a difficult optimal problem, it deserves further in-depth study.
The following theorem presents one possible way of choosing the tolerances { k } and { k } such that the original convergence rate of the ADI iteration can be asymptotically recovered. Theorem 3. Let the assumptions in Theorem 2 be satisfied. Suppose that both { 1 (k)} and { 2 (k)} are nondecreasing and positive sequences satisfying 1 (k) ≥ 1, 2 (k) ≥ 1, and lim k→∞ sup 1 (k) = lim k→∞ sup 2 (k) = +∞, and that both 1 and 2 are real constants on the interval (0,1) satisfying where t 1 and t 2 are nonnegative constants. Then we have In particular, we have i.e. the convergence rate of the IADI iterative method is asymptotically the same as that of the ADI iterative method.
Proof. The proof is a verbatim of theorems 3.3 and 3.4 in Reference 17 and thus omitted. ▪

NUMERICAL EXAMPLES
In this section, we use some examples to illustrate the effectiveness of ADI and IADI iterations for solving the Sylvester equation (1).
In actual computations, all iterations are started from the zero matrix, performed in MATLAB with machine precision 10 −16 , and stopped if the norm of the current residual matrix satisfies ||R (k) || F / || R (0) || F ≤ 10 −6 , where, following the definition, For convenience, throughout our numerical experiments, we denote by IT the number of the iteration steps, by CPU the computing time (in seconds), and by exp and exp the experimentally found local optimal values § of the iteration parameters of the ADI iteration, respectively.  Such a class of problems arise frequently in the preconditioned Krylov subspace iterative methods, see Reference 4 and reference therein. For comparison, we also test the HSS solver. The numerical results obtained by the ADI method and the HSS solver are listed in Table 1 with r = 1, Table 2 with r = 0.1, Table 3 with r = 0.01, where let exp = exp , and exp is the numerically found optimal value of the shift parameter of the HSS iterations in Reference 4. We can see that both the number of iteration steps and the runtime by ADI iteration are much less than those by HSS iteration in all cases.
In Figure 1, we depict the convergence behavior of the ADI and HSS iterations for Example 1, which are denoted by +++ and • • • curves, respectively. From Figure 1, we can see that the ADI iteration has a better convergence behavior than the HSS iteration. We can also observe that the HSS is not sensitive to the nonsymmetric part, as mentioned in, 4 where the author pointed out that the convergence depends only on the spectrums of the Hermitian parts (symmetric parts in our example).

Example 2. We consider the Sylvester equation (1) and the matrices
with L the strictly lower triangular matrix having ones in the lower triangle part and t is a problem parameter to be specified in actual computations.
The Sylvester equation in Example 2 is solved by the ADI, IADI, and IHSS iterative methods, respectively. Parameters 1 exp and 1 exp are the numerically found optimal values for IHSS iteration. The corresponding results are listed in Table 4, where we set t = n, r = 1 n , and we use the GMRES as the inner iteration scheme and set k = k = 0.01, for k = 0,1,2, … . From Table 4, we can observe that the IADI iteration is usually superior to the ADI iteration in terms of computation efficiency. Moreover, the IADI iteration as well as the ADI iteration outperform the IHSS iteration in terms of both the number of iterations and the computation efficiency, especially when the order of the coefficient matrices is large enough.

CONCLUSIONS
In this paper, we have revisited ADI iteration for solving the continuous Sylvester equation AX + XB = C, where the coefficient matrices A and B are assumed to be positive semi-definite matrices (not necessarily Hermitian), and at least one of them to be positive definite. For such a class of Sylvester equations, we have analyzed the convergence of the ADI iteration and derived an upper bound of its contraction factor. To reduce the computational complexity, we have also proposed the IADI iteration whose convergence has been analyzed in detail. We have presented some numerical experiments to illustrate the effectiveness of both ADI and IADI iterations.
In order to get the practical choices of the iteration parameters and , we first get the maximal and minimal eigenvalues of A and B via power method and inverse iteration, then get and Δ according to formulas (2.6) to (2.10), and finally get the approximations to and . The scalars exp and exp in text are obtained by searching the optimal values of the iteration parameters for the ADI iteration in two intervals ( − 1, + 1) and ( − 1, + 1) with stepsize 0.1, respectively. We notice that those iteration parameters exp and exp used in our numerical experiments for the ADI iteration are the experimentally locally optimal values, not the globally optimal values. Therefore, we can expect the ADI and IADI iterations with exact optimal iteration parameters to have better convergence behavior and computational efficiency. However, it is an important and hard task to find the optimal and which strongly depend on the specific structures and properties of the coefficient matrices A and B and need further in-depth study from the viewpoint of both theory and computations.
We remark here that our convergence theory regarding ADI iterations with two parameters for solving Sylvester equations can be easily extended to solve the general positive definite linear system, Kx = b, if K = U + V with U and V being positive definite. In this sense, we can say we have generalized the convergence theory in Reference 12 concerning ADI iterations with two parameters for solving Hermitian positive definite linear system Kx = b to the case of K being non-Hermitian positive definite.