A Preconditioned Iteration Method for Solving Sylvester Equations

Abstract

A preconditioned gradient-based iterative method is derived by judicious selection of two auxil- iary matrices. The strategy is based on the Newton’s iteration method and can be regarded as a generalization of the splitting iterative method for system of linear equations. We analyze the convergence of the method and illustrate that the approach is able to considerably accelerate the convergence of the gradient-based iterative method.

1. Introduction

In this paper, we consider preconditioned iterative methods for solving Sylvester equations of form where , , and , with generally. A Lyapunov equation is a special case of (1.1) with , , and . Such kind of problems frequently arise from many areas of applications in control and system theory [1], stability of linear systems [2], analysis of bilinear systems [3], power systems [4], signal and image processing [5], and so forth.

Throughout the paper, we assume that Sylvester equation (1.1) possess a unique solution, that is, where and denote the spectra of and , respectively [6]. In theory, the exact solution of (1.1) can be computed by “linearization,” that is, by solving an equivalent system of linear equations of form where , with , and is the Kronecker product. However, the above direct method requires considerable computational efforts, due to the high dimension of the problem.

For small-to-medium-scale problems of the form (1.1), direct approaches such as Bartels-Stewart method [7, 8] and Hessenberg-Schur method [6, 9] have been the methods of choice. The main idea of these approaches is to transform the original linear system into a structured system that can be solved efficiently by forward or backward substitutions.

In the numerical linear community, iterative methods are becoming more and more popular. Several iterative schemes for Sylvester equations have been proposed; see, for example, [10–15]. Recently, Some gradient based iterative methods [3–5, 16–26] have been investigated for solving general coupled matrix equations and general matrix equations. For Sylvester equations of form (1.1), the gradient based iterative methods use a hierarchical identification principle to compute the approximate solution. The convergence condition of these methods is investigated in [16–18]. It is proved that the gradient based iterative methods are convergent under certain conditions. However, we observe that the convergence speed of the gradient based iterative methods is generally very slow, which is similar to the behavior of classical iterative methods applied to systems of linear equations. In this paper, we consider preconditioning schemes for solving Sylvester equations of form (1.1). We illustrated that the preconditioned gradient based iterative methods can be derived by selecting two auxiliary matrices. The selection of preconditioners is natural from the view point of splitting iteration methods for systems of linear equations. The convergent property of the preconditioned method is proved and the optimal relaxation parameter is derived. The performance of the method is compared with the original method in several examples. Numerical results show that preconditioning is able to considerably speed up the convergence of the gradient based iterative method.

The paper is organized as follows. In Section 2, a gradient based iterative method is recalled, and the preconditioned gradient based method is introduced and analyzed. In Section 3, performance of the preconditioned gradient based method is compared with the unpreconditioned one, and the influence of an iterative parameter is experimentally studied. Finally, we conclude the paper in Section 4.

2. A Brief Review of the Gradient Based Iterative Method

We firstly recall an iterative method proposed by Ding and Chen [18] for solving (1.1). The basic idea is regarding (1.1) as two linear matrix equations: Then define two recursive sequences where is the iterative step size. The above procedures can be regarded as two separate iterative procedures for solving two matrix equations in (2.1).

With and at hand, then the th approximate solution can be defined by taking the average of two approximate solutions, that is, By selecting an appropriate initial approximate solution , and using to substitute in (2.2) and in (2.3), then the above (2.2)-(2.3) constitute the gradient based iterative method proposed in [18]. It is shown [18] that the gradient based iterative algorithm converges as long as where is the largest eigenvalue of .

3. Preconditioned Gradient Based Iterative Method

We start with a general algebraic equation [27]

Suppose is an approximate solution for (3.1), for example, (). Then it follows that the accuracy of can be improved by the following scheme: with being a Lagrangian multiplier. It is well known that the optimal value of is determined by From the above condition, the Newton's iteration follows:

Now, let us consider a general system of linear equations of form where is a general square matrix. Let , and be an initial approximate solution. Then by one-step Newton's iteration (3.4), we can theoretically find the exact solution However, the previous procedure has the same difficulty as solving (3.5) by direct inverse.

By utilizing the idea of preconditioning, the iteration formula (3.6) can be rebuilt by replacing by an approximate matrix . Taking as a starting vector and , then it follows that which is the preconditioned iterative method. Formula (3.7) can also be derived from the splitting iteration [28] with More generally, a relaxation parameter can be introduced as follows:

Recall that Sylvester equations (1.1) can be rewritten as the following two fictitious matrix equations: where and .

By (3.9),we define preconditioned iterations as follows: where and are well-chosen preconditioners for and , respectively. The th approximate solution can be defined by taking the average of two approximate solutions [16, 18], for example, By selecting two initial approximate solutions, the above procedures (3.11)–(3.13) constitute the framework of the Newton's iteration method.

The above process can be accomplished by the following algorithm.

Algorithm 3.1. Preconditioned gradient based iterative algorithm (PGBI) can be given as follows. (1)Give two initial approximate solutions and . (2)for , until converges (3)(4)(5)(6)end.

Remark 3.2. The above algorithm follows the framework of the gradient based iterative method proposed in [16–18]. However, in [16, 18] the matrix is chosen as (), and in [17] the matrix is chosen as (). The previous selection of the matrices is obviously not wise. In this paper, we have illustrated that and should be the reasonable preconditioner for and , respectively.
Replacing by , then the above algorithm can be used to solve generalized Sylvester equation of form We will present an example of form (3.14) in the last section of the paper.

Lemma 3.3 (see [29]). Leting , then is convergent if and only if .

Based on Theorem 3.4 in [18], we have the following theorem.

Theorem 3.4. Suppose Sylvester equation (1.1) has a unique solution , and then the iterative sequence generated by Algorithm 3.1 converges to , that is, converges to zero for any initial value .

Proof. Since , it follows that
Letting and submitting into the above formula, we have
Using to subtract both sides of the above equation, we have
Let be defined as the vector formed by stacking the columns of on the top of one another, that is,
Then (3.18) is equivalent to where . Therefore, the iteration is convergent if and only if , that is, The proof is complete.

Remark 3.5. The choice of parameter is an important issue. We will experimentally study its influence on the convergence. However, the parameter is problem dependent; so seeking a parameter that is suitable for a broad range of problems is a difficult task.

4. Numerical Examples

In the following tests, the parameter is set to be in the gradient based iterative method (GBI), and is chosen experimentally in the preconditioned gradient based iterative method (PGBI). In all the test, the stopping criterion is set to be . The exact solution is set as such that the right hand side , and the initial approximate solution is set to be . The relative error norms are recorded and plotted by Matlab command semilogy [30]. We always choose ILU(0) [31] as the preconditioner. However, other suitable preconditioners can also be adapted.

Example 4.1. The coefficient matrices used in this example are taken from [18]. We outline the matrix for completeness: In the matrices, there is a parameter which can be used to change the weight of the diagonal entries. The case with . For , the behavior of the error norms is recorded in Figure 1(a). From this figure, we can see that the convergence of GBI method tends to slow down after some iteration steps. The PGBI method converges linearly, and much faster than GBI method. In Figure 1(b), the convergence curves with different parameter are recorded. For this problem, we can see that the optimal value of is very close to 0.5. By comparing the convergence of GBI method, we can see that PGBI is able to converge within 300 iteration steps, whereas GBI needs more than 1000 iteration steps. Therefore, even not with the optimal , the convergence of PGBI method is also much faster than that of GBI method.

(a)

(b)

(a)
(b)

Figure 1 

Comparison of convergence curves using GBI and PGBI (a); the influence of parameter on PGBI method (b).

Example 4.2. In this example, we test a problem with , where coefficient matrix is generated from the discretization of the following two-dimensional Poisson problem: posed on the unit square with Dirichlet boundary conditions Discretizing this problem by the standard second-order finite difference (FE) scheme on a uniform grid with step size in each direction, we obtain a coefficient matrix with and , . By choosing , we compared the performance of GBI method and the preconditioned GBI method. The convergence curves are recorded in Figure 2(a). Figure 2(b), the influence of parameter is investigated. From this figure we can see that the optimal is close to 0.4. By comparing the convergence of GBI method, it is easy to see that for a wide range of the preconditioned GBI method is much better than GBI method.

(a)

(b)

(a)
(b)

Figure 2 

Comparison of convergence curves using GBI and PGBI (a); the influence of parameter on PGBI method (b).

Example 4.3. In this example, we consider the convection diffusion equation with Dirichlet boundary conditions [32] The equation is discretized by using central finite differences on , with mesh size in the X-direction, and in the Y-direction, produces two tridiagonal matrices and given by By taking , , and , we have tested the performance of GBI method and the preconditioned GBI method. The convergence curves are recorded in Figure 3(a). From this figure, we can see that the GBI method converges very slowly and nearly stagnate. The preconditioned GBI method has much better performance. We also investigate the influence of parameter on this problem. The convergence behavior with different is recorded in Figure 3(b). From this figure, we can see that the parameter becomes very sensitive when it is larger than the optimal value. Therefore, a relative small parameter is more reliable.

(a)

(b)

(a)
(b)

Figure 3 

Comparison of convergence curves using GBI and PGBI (a); the influence of parameter on PGBI method (b).

Example 4.4. In this example, we intend to test the algorithm for solving generalized Sylvester equation (3.14). All the initializations are the same except setting . The coefficient matrix has the following structure: and . We set in this example. The tested results are shown in Figure 4. The faster convergence behavior is observable from Figure 4(a).

(a)

(b)

(a)
(b)

Figure 4 

Comparison of convergence curves using GBI and PGBI (a); the influence of parameter on PGBI method (b).

5. Conclusions

In this paper, a preconditioned gradient based iteration (PGBI) method is proposed. The convergence of PGBI is analyzed. The choice of parameter is an important issue, and its influence is experimentally studied. The principle idea of this paper can be extended to the more general setting like coupled Sylvester matrix equations.

Acknowledgment

This work is supported by NSF no. 11101204, and XJTLU RDF.

References

1. R. Bhatia and P. Rosenthal, “How and why to solve the operator equation $AX-XB=Y$,” The Bulletin of the London Mathematical Society, vol. 29, no. 1, pp. 1–21, 1997. View at: Publisher Site | AX-XB=Y&author=R. Bhatia &author=P. Rosenthal&publication_year=1997" target="_blank">Google Scholar | Zentralblatt MATH
2. B. N. Datta, Numerical Methods for Linear Control Systems, Elsevier Academic Press, 2003.
3. L. Xie, Y. J. Liu, and H.-Z. Yang, “Gradient based and least squares based iterative algorithms for matrix equations $AXB+C{X}^{T}D=F$,” Applied Mathematics and Computation, vol. 217, no. 5, pp. 2191–2199, 2010. View at: Publisher Site | AXB+CXTD=F&author=L. Xie&author=Y. J. Liu&author=&author=H.-Z. Yang&publication_year=2010" target="_blank">Google Scholar | Zentralblatt MATH
4. L. Xie, J. Ding, and F. Ding, “Gradient based iterative solutions for general linear matrix equations,” Computers & Mathematics with Applications, vol. 58, no. 7, pp. 1441–1448, 2009. View at: Publisher Site | Google Scholar | Zentralblatt MATH
5. J. Ding, Y. J. Liu, and F. Ding, “Iterative solutions to matrix equations of the form $A_{i}X{B}_i=F_i$,” Computers & Mathematics with Applications, vol. 59, no. 11, pp. 3500–3507, 2010. View at: Publisher Site | AiXBi=Fi&author=J. Ding&author=Y. J. Liu&author=&author=F. Ding&publication_year=2010" target="_blank">Google Scholar | Zentralblatt MATH
6. G. H. Golub, S. Nash, and C. Van Loan, “A Hessenberg-Schur method for the problem $AX-XB=C$,” IEEE Transactions on Automatic Control, vol. 24, no. 6, pp. 909–913, 1979. View at: Publisher Site | AX-XB=C&author=G. H. Golub&author=S. Nash&author=&author=C. Van Loan&publication_year=1979" target="_blank">Google Scholar | Zentralblatt MATH
7. R. H. Bartels and G. W. Stewart, “Algorithm 432: solution of the matrix equation $AX-XB=C$,” Communications of the ACM, vol. 15, pp. 820–826, 1972. View at: AX-XB=C&author=R. H. Bartels &author=G. W. Stewart&publication_year=1972" target="_blank">Google Scholar
8. D. C. Sorensen and Y. K. Zhou, “Direct methods for matrix Sylvester and Lyapunov equations,” Journal of Applied Mathematics, vol. 2003, no. 6, pp. 277–303, 2003. View at: Publisher Site | Google Scholar | Zentralblatt MATH
9. W. H. Enright, “Improving the efficiency of matrix operations in the numerical solution of stiff ordinary differential equations,” ACM Transactions on Mathematical Software, vol. 4, no. 2, pp. 127–136, 1978. View at: Publisher Site | Google Scholar | Zentralblatt MATH
10. D. Calvetti and L. Reichel, “Application of ADI iterative methods to the restoration of noisy images,” SIAM Journal on Matrix Analysis and Applications, vol. 17, no. 1, pp. 165–186, 1996. View at: Publisher Site | Google Scholar | Zentralblatt MATH
11. D. Y. Hu and L. Reichel, “Krylov-subspace methods for the Sylvester equation,” Linear Algebra and Its Applications, vol. 172, no. 15, pp. 283–313, 1992. View at: Publisher Site | Google Scholar | Zentralblatt MATH
12. P. Kirrinnis, “Fast algorithms for the Sylvester equation $AX-X{B}^T=C$,” Theoretical Computer Science, vol. 259, no. 1-2, pp. 623–638, 2001. View at: Publisher Site | AX-XBT=C&author=P. Kirrinnis&publication_year=2001" target="_blank">Google Scholar | Zentralblatt MATH
13. G. Starke and W. Niethammer, “SOR for $AX-XB=C$,” Linear Algebra and Its Applications, vol. 154/156, pp. 355–375, 1991. View at: Publisher Site | AX-XB=C&author=G. Starke &author=W. Niethammer&publication_year=1991" target="_blank">Google Scholar | Zentralblatt MATH
14. Q. Niu, X. Wang, and L.-Z. Lu, “A relaxed gradient based algorithm for solving Sylvester equations,” Asian Journal of Control, vol. 13, no. 3, pp. 461–464, 2011. View at: Publisher Site | Google Scholar
15. Z.-Z. Bai, “On Hermitian and Skew-Hermitian splitting iteration methods for continuous Sylvester equations,” Journal of Computational Mathematics, vol. 29, no. 2, pp. 185–198, 2011. View at: Google Scholar
16. F. Ding and T. Chen, “On iterative solutions of general coupled matrix equations,” SIAM Journal on Control and Optimization, vol. 44, no. 6, pp. 2269–2284, 2006. View at: Publisher Site | Google Scholar | Zentralblatt MATH
17. F. Ding and T. Chen, “Iterative least-squares solutions of coupled Sylvester matrix equations,” Systems & Control Letters, vol. 54, no. 2, pp. 95–107, 2005. View at: Publisher Site | Google Scholar | Zentralblatt MATH
18. F. Ding and T. Chen, “Gradient based iterative algorithms for solving a class of matrix equations,” IEEE Transactions on Automatic Control, vol. 50, no. 8, pp. 1216–1221, 2005. View at: Publisher Site | Google Scholar
19. F. Ding, P. X. Liu, and J. Ding, “Iterative solutions of the generalized Sylvester matrix equations by using the hierarchical identification principle,” Applied Mathematics and Computation, vol. 197, no. 1, pp. 41–50, 2008. View at: Publisher Site | Google Scholar | Zentralblatt MATH
20. F. Ding and T. Chen, “Gradient-based identification methods for Hammerstein nonlinear ARMAX models,” Nonlinear Dynamics, vol. 45, no. 1-2, pp. 31–43, 2006. View at: Publisher Site | Google Scholar | Zentralblatt MATH
21. F. Ding and T. Chen, “Gradient and least-squares based identification methods for OE and OEMA systems,” Digital Signal Processing, vol. 20, no. 3, pp. 664–677, 2010. View at: Google Scholar
22. Y. J. Liu, J. Sheng, and R. F. Ding, “Convergence of stochastic gradient estimation algorithm for multivariable ARX-like systems,” Computers & Mathematics with Applications, vol. 59, no. 8, pp. 2615–2627, 2010. View at: Publisher Site | Google Scholar | Zentralblatt MATH
23. Y. J. Liu, Y. S. Xiao, and X. L. Zhao, “Multi-innovation stochastic gradient algorithm for multiple-input single-output systems using the auxiliary model,” Applied Mathematics and Computation, vol. 215, no. 4, pp. 1477–1483, 2009. View at: Publisher Site | Google Scholar | Zentralblatt MATH
24. F. Ding, G. J. Liu, and X. P. Liu, “Parameter estimation with scarce measurements,” Automatica, vol. 47, no. 8, pp. 1646–1655, 2011. View at: Google Scholar
25. F. Ding and T. W. Chen, “Performance analysis of multi-innovation gradient type identification methods,” Automatica, vol. 43, no. 1, pp. 1–14, 2007. View at: Publisher Site | Google Scholar | Zentralblatt MATH
26. F. Ding, Y. J. Liu, and B. Bao, “Gradient based and least squares based iterative estimation algorithms for multi-input multi-output systems,” Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 226, no. 1, pp. 43–55, 2012. View at: Google Scholar
27. M. Y. Waziri, W. J. Leong, and M. Mamat, “A two-step matrix-free secant method for solving large-scale systems of nonlinear equations,” Journal of Applied Mathematics, vol. 2012, Article ID 348654, 9 pages, 2012. View at: Publisher Site | Google Scholar
28. L. A. Hageman and D. M. Young, Applied Iterative Methods, Academic Press, New York, NY, USA, 1981.
29. R. S. Varga, Matrix Iterative Analysis, Springer, Berlin, Germany, 2nd edition, 2000. View at: Publisher Site
30. The MathWorks, Inc. MATLAB 7, September 2004.
31. Y. Saad, Iterative Methods for Sparse Linear Systems, PWS, Boston, Mass, USA, 1996.
32. A. Bouhamidi and K. Jbilou, “A note on the numerical approximate solutions for generalized Sylvester matrix equations with applications,” Applied Mathematics and Computation, vol. 206, no. 2, pp. 687–694, 2008. View at: Publisher Site | Google Scholar | Zentralblatt MATH

Copyright

Copyright © 2012 Jituan Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.