- Email: [email protected]
Successive projection iterative method for solving matrix equation AX=B
Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Successive projection iterative method for solving matrix equation AX = B Fan-Liang Li a,∗ , Li-sha Gong b , Xi-Yan Hu c , Lei Zhang c a
Institute of Mathematics and Physics, School of Sciences, Central South University of Forestry and Technology, Changsha, 410004, PR China
b
School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu, Sichuan, 610054, PR China
c
College of Mathematics and Econometrics, Hunan University, Changsha, 410082, PR China
article
info
Article history: Received 31 August 2009 Received in revised form 31 December 2009 MSC: 15A24 11Y40
abstract In this paper, we present a new iterative method (successive projection iterative method) to solve matrix equation AX = B, where A is a symmetric positive definite (SPD) matrix. Based on this method an algorithm is proposed and proved to be convergent. In addition, analysis of the algorithm and numerical experiments are also given to show the efficiency of the method. © 2010 Elsevier B.V. All rights reserved.
Keywords: Matrix equation Symmetric positive definite matrix A-orthogonal matrix group Successive projection iterative method
1. Introduction The problem of solving matrix equation is a hot topic in the field of numerical algebra in recent years. Many authors studied this problem by using the special structure of matrix and the decomposition of matrix pairs, and a series of meaningful results were achieved [1–6]. For example, Golub [1] and Zhou [2] discussed the solution of matrix equation with the singular value decomposition (SVD) of matrix; Loan [3] and Chu [4] put forward and discussed the generalized singular value decomposition (GSVD) of matrix pairs to solve matrix equation; Golub [5] and Higham [6] considered matrix equation with the canonical correlation decomposition of matrix pairs (CCD) and the quotient singular value decomposition (QSVD) of matrix pairs, respectively. The method in these papers is called direct method. With this method, the authors can present the solvability conditions and the expression of the solution for matrix equations directly. However, in these cases, the solvability conditions and the expression of the solution for matrix equations are usually very complex which may cause huge error during calculating. Furthermore, some matrix equations cannot be solved with direct method at all. For example, it is difficult to obtain the centrosymmetric solution, the reflexive solution and the double symmetric solution of matrix equation AXB = C with direct method. In this paper, based on the Gauss–Seidel iterative method and successive projection method, which are used to solve linear equations Ax = b [7,8], we propose an improved method to solve matrix equation AX = B. Let Rn×m be the set of all n × m real matrices, In be the identity matrix of order n, SRn+×n denote the set of all n × n SPD matrices, M T , M + , r (M ) and tr(M ) represent the transpose, the Moore–Penrose generalized inverse, rank and trace of M,
∗
Corresponding author. Tel.: +86 15973179263; fax: +86 073185010892. E-mail address: [email protected] (F.-L. Li).
0377-0427/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2010.03.008
2406
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
respectively. For M = (mij ), N = (nij ) ∈ Rn×m , hM , N i = tr(N T M ) denote the inner product of matrix M and N. The induced 1
1
matrix norm is called Frobenius norm, i.e. kM k = hM , M i 2 = (tr(M T M )) 2 , then Rn×m is a Hilbert inner product space. Definition 1. Let M , N ∈ Rn×m , if hM , N i = tr(N T M ) = 0, then M , N is called orthogonal; for a given SPD matrix A, if hM , N iA = tr(N T AM ) = 0, then M , N is called A-orthogonal. Definition 2. Let A ∈ SRn+×n , if M1 , M2 , . . . , Mp ∈ Rn×m , and Mi 6= 0, i = 1, 2, . . . , p, for any Mi , Mj , hMi , Mj i = tr(MjT Mi ) =
0(hMi , Mj iA = tr(MjT AMi ) = 0), i 6= j, i, j = 1, 2, . . . , p, then M1 , M2 , . . . , Mp are called orthogonal (A-orthogonal) matrix group. It is obviously that orthogonal (A-orthogonal) matrix group M1 , M2 , . . . , Mp are linearly independent. That is, if c1 M1 + c2 M2 + · · · + cp Mp = 0 holds, then c1 = c2 = · · · = cp = 0. Definition 3. Let A ∈ SRn+×n , if M1 , M2 , . . . , Mp ∈ Rn×m are orthogonal (A-orthogonal) matrix group, and kMi k = 1(kMi kA = tr(MiT AMi )
12
= 1), i = 1, 2, . . . , p, then M1 , M2 , . . . , Mp are called orthonormal (A-orthonormal) matrix group.
In this paper, we consider the following problem. Problem. Giving A ∈ SRn+×n , B ∈ Rn×m , find X ∈ Rn×m such that AX = B.
(1.1)
This paper is organized as follows. In Section 2, successive projection iterative method to solve (1.1) is presented. In Section 3, an algorithm is derived by giving an A-orthonormal matrix group and proved to be convergent. In the end, analysis of the algorithm and numerical experiments are given. 2. Successive projection iterative method Let A ∈ SRn+×n , if M1 , M2 , . . . , Mp ∈ Rn×m is an A-orthonormal matrix group, and let L = K = span{M1 , M2 , . . . , Mp }, p ≤ m. Then the method of successive projection iterative can be described as follows: let Xk , Xk+1 be the k, k + 1 step iterative approximate solution of (1.1), respectively, then Xk+1 satisfies the following conditions. Xk+1 ∈ Xk + K ,
Rk+1 = B − AXk+1 ⊥ L.
From this, we can suppose that (k)
(k)
Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp(k) Mp , (k)
(k)
(k)
where α1 , α2 , . . . , αp are constants. Hence, the k + 1 step residual matrix (k)
(k)
Rk+1 = B − AXk+1 = B − A(Xk + α1 M1 + α2 M2 + · · · + αp(k) Mp )
= Rk − α1(k) AM1 − α2(k) AM2 − · · · − αp(k) AMp , where Rk = B − AXk . It is easy to see hRk+1 , Mi i = 0, i = 1, 2, . . . , p from Rk+1 = B − AXk+1 ⊥ L, and L = K = span{M1 , M2 , . . . , Mp }. Which implies
(k) α1 hAM1 , M1 i + α2(k) hAM2 , M1 i + · · · + αp(k) hAMp , M1 i = hRk , M1 i, α (k) hAM1 , M2 i + α (k) hAM2 , M2 i + · · · + α (k) hAMp , M2 i = hRk , M2 i, p 1 2 .. . (k) α1 hAM1 , Mp i + α2(k) hAM2 , Mp i + · · · + αp(k) hAMp , Mp i = hRk , Mp i. Noticing that M1 , M2 , . . . , Mp ∈ Rn×m is an A-orthonormal matrix group, we have
α1(k) = tr(M1T Rk ), α2(k) = tr(M2T Rk ), . . . , αp(k) = tr(MpT Rk ). Up to now, we can design the following method to solve (1.1). Successive projection iterative method. 1. Choose an initial matrix X1 ∈ Rn×m , compute R1 = B − AX1 , k := 1. 2. If kRk k ≤ error bound, then X = Xk , stop iterating; otherwise go to 3.
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
(k)
2407
(k)
(k)
3. Give an A-orthonormal matrix group M1 , M2 , . . . , Mp ∈ Rn×m , compute Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp Mp , (k)
where αi = tr( ), i = 1, 2, . . . , p, p ≤ m. 4. Let k := k + 1, and go to 2. MiT Rk
Theorem 1. In the above method, if Xk , Xk+1 are the k, k + 1 step iterative approximate solution of (1.1) respectively, and X is the exact solution of (1.1). Then
kDk kA ≥ kDk+1 kA , where Dk = X − Xk , Dk+1 = X − Xk+1 . (In fact, Dk is the difference between the exact solution and k step iterative approximate solution, and we call it the k step error.). (k)
(k)
(k)
Proof. According to Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp Mp , we have (k)
(k)
Dk+1 = X − Xk+1 = X − (Xk + α1 M1 + α2 M2 + · · · + αp(k) Mp )
= Dk − (α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp ). This implies
hADk+1 , Dk+1 i = hADk − (α1(k) AM1 + α2(k) AM2 + · · · + αp(k) AMp ), Dk − (α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp )i = hADk , Dk i − 2hADk , α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp i +
p X (αi(k) )2 . i=1
Combining ADk = AX − AXk = B − AXk = Rk , it is easy to see
hADk+1 , Dk+1 i = hADk , Dk i − 2hRk , α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp i +
p X (αi(k) )2 . i =1
(k)
Noticing that αi
= hRk , Mi i, i = 1, 2, . . . , p, we have
hADk+1 , Dk+1 i = hADk , Dk i −
p X (αi(k) )2 , i =1
i.e.
kDk k2A − kDk+1 k2A =
p X (αi(k) )2 ≥ 0,
(2.1)
i=1
and this gives the conclusion.
(αi(k) )2 = 0, i.e. αi(k) = 0, i = 1, . . . , p, then Xk+1 = Xk , and we do not have any reduction Pp (k) 2 2 2 of k + 1 step error comparing to k step error in A-norm meaning. But if i=1 (αi ) > 0, i.e. kDk kA > kDk+1 kA , then the Theorem 1 shows that if
Pp
i=1
error is reduced. Theorem 1 also shows that the sequence {kDk k2A } is convergent. Hence, we can suppose as follows.
kDk k2A −→ a (k −→ ∞),
(2.2)
where a ≥ 0. Combining (2.1) and (2.2) gives the following results.
αi(k) −→ 0 (k −→ ∞).
(2.3)
From the theory of projection, it is easy to derive the following theorem. Theorem 2. In the above method, if Xk , Xk+1 are the k, k + 1 step iterative approximate solution of (1.1), respectively, then for any Xk˜+1 ∈ Xk + K , Xk+1 satisfies
kB − AXk+1 k =
min ∀Xk˜+1 ∈Xk +K
kB − AXk˜+1 k.
Theorem 3. In the above method, if Xk , Xk+1 are the k, k + 1 step iterative approximate solution of (1.1), respectively, then for any Xk˜+1 ∈ Xk + K , Xk+1 satisfies
kX − Xk+1 kA =
min ∀Xk˜+1 ∈Xk +K
kX − Xk˜+1 kA ,
where X is the exact solution of (1.1).
2408
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
Proof. Definition 3 shows that kX − Xk˜+1 kA = tr[(X − Xk˜+1 )T A(X − Xk˜+1 )], and
kB − AXk˜+1 k = tr[(B − AXk˜+1 )T (B − AXk˜+1 )] = tr[(X − Xk˜+1 )T AT A(X − Xk˜+1 )] = tr[(X − Xk˜+1 )T A2 (X − Xk˜+1 )]. A ∈ SRn+×n implies that there exists an n × n orthogonal matrix P such that
λ1
0 A=P .. .
λ2 .. .
··· ··· .. .
0
0
···
0
0 0
T .. P , . λn
where λi > 0, i = 1, 2, . . . , n, are eigenvalues of matrix A. Thus A2 can be written as
2 λ1 0 A2 = P . ..
λ22 .. .
··· ··· .. .
0
0
···
0
0 0
T .. P . . λ2n
This implies
λ1
0 ˜ T kX − Xk˜+1 kA = tr (X − Xk+1 ) P .. .
λ2 .. .
··· ··· .. .
0
0
···
λ 0 kB − AXk˜+1 k = tr (X − Xk˜+1 )T P . ..
0
λ22 .. .
0
0
··· ··· .. . ···
2 1
0
0 0
T ˜ .. P (X − Xk+1 ) , . λn 0 0
T .. P (X − Xk˜+1 ) . . 2 λn
Let F (Xk˜+1 ) = kX − Xk˜+1 kA , G(Xk˜+1 ) = kB − AXk˜+1 k, then it can be verified without difficulty that function F (Xk˜+1 ) and G(Xk˜+1 ) have the same point of minimum. Combining Theorem 2, we can obtain the conclusion of Theorem 3. Theorem 3 shows that Xk+1 is the result of projection of Xk onto K (orthogonal to L) iff it minimizes the A-norm of the k + 1 step error over Xk + K . Theorem 3 also shows that if K1 = L1 ⊆ K2 = L2 , then the reduction of A-norm of the k step error obtained by subspaces K2 = L2 is more or equal to that of subspaces K1 = L1 . Hence by increasing the value of p, the convergence rate may increase. 3. An algorithm and numerical results Analysis of the above method Suppose A ∈ SRn+×n , the key problem of the above method is how to obtain an A-orthonormal matrix group M1 , M2 , . . . , Mp ∈ Rn×m in step 3. In the following section, we give a method to obtain a special A-orthonormal matrix (k)
(k)
(k)
(k)
(k)
group M1 , M2 , . . . , Mp ∈ Rn×m . Let Rk = (rij )n×m , and suppose |ri1 j1 |, |ri2 j2 |, . . . , |rip jp | be the p largest values in all |rij | and (k)
(k)
(k)
satisfy |ri1 j1 | ≥ |ri2 j2 | ≥ · · · ≥ |rip jp |, jl 6= jh , jl , jh = j1 , . . . , jp . Let M1 = √
1
1 1 Ei1 j1 , M2 = √ Ei2 j2 , . . . , Mp = √ Eip jp , ai1 i1 ai 2 i 2 aip ip
(3.1)
where Eij is an n × m matrix with the (i, j) element equal to 1, and others equal to zero, ai1 i1 , ai2 i2 , . . . , aip ip are the i1 th, i2 th, . . . , ip th element of A in the main diagonal line. It is can be obtained immediately that M1 , M2 , . . . , Mp is an A-orthonormal matrix group with p ≤ m, and (k)
ri j
αl(k) = tr(MlT Rk ) = √ l l , ail il
l = 1 , 2 , . . . , p.
Combining (2.3) gives Rk −→ 0 (k −→ ∞), which implies Xk −→ X
(k −→ ∞),
(3.2)
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
2409
where, X is the exact solution of (1.1). The above analysis shows that if the A-orthonormal matrix group M1 , M2 , . . . , Mp is given by (3.1), then we can provide the following algorithm which ensures that iterative sequence {Xk } converge to X . Algorithm. 1. Input an initial matrix X1 ∈ Rn×m , compute R1 = B − AX1 , k := 1. 2. If kRk k ≤ error bound, then X = Xk , stop iterating; otherwise go to 3. 3. Compute M1 , M2 , . . . , Mp according to (3.1). (k)
(k)
(k)
(k)
4. Compute Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp Mp , where αi 5. Let k := k + 1, and go to 2.
= tr(MiT Rk ), i = 1, 2, . . . , p, p ≤ m.
Analysis of the algorithm In the above algorithm, step 1 requires n2 m arithmetic operations (additions and multiplications), step 2 requires n2 m arithmetic operations (additions and multiplications), step 3 requires n arithmetic operations (multiplications), and step 4 requires nmp arithmetic operations (additions and multiplications). Numerical examples (n = 10, m = 6) Let
40 10 0.5 0.5 0.5 A= 0.5 0.5 0.5 0.5 0.5
0.5 10 40 10 0.5 0.5 0.5 0.5 0.5 0.5
10 40 10 0.5 0.5 0.5 0.5 0.5 0.5 0.5
54 63.5 63.5 63.5 63.5 B= 63.5 63.5 63.5 63.5
54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
54
0.5 0.5 10 40 10 0.5 0.5 0.5 0.5 0.5
54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
0.5 0.5 0.5 10 40 10 0.5 0.5 0.5 0.5 54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
0.5 0.5 0.5 0.5 10 40 10 0.5 0.5 0.5
0.5 0.5 0.5 0.5 0.5 10 40 10 0.5 0.5
54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
0.5 0.5 0.5 0.5 0.5 0.5 10 40 10 0.5
0.5 0.5 0.5 0.5 0.5 0.5 0.5 10 40 10
0.5 0.5 0.5 0.5 0.5 , 0.5 0.5 0.5 10 40
54 63.5 63.5 63.5 63.5 , 63.5 63.5 63.5 63.5 54
it is clear that the exact solution X of (1.1) is
1 1 1 X = 1 T
1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 . 1 1 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 , 0 0 0
Input an initial guess X1
0 X1T
0 0 = 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
and use kRk = B − AXk k < 0.01 as the stopping criterion. By running some MATLAB codes, we obtain the following table with p = 4, 5, 6. Table 1 confirms that the above algorithm is convergent and by increasing the value of p, the convergence rate may increase. Conclusion In this paper, based on the Gauss–Seidel iterative method and successive projection method, which are used to solve linear equation Ax = b in [7,8], we propose a successive projection iterative method to solve matrix equation AX = B. Compared to [7,8], this paper has two important achievements. One is that the successive projection method is used to solve matrix equation AX = B in this paper instead of linear equation Ax = b in [7,8]. The other is we provide the method
2410
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
Table 1 Numerical results for example. p
Step
kRk k
Cputime
4 5 6
104 83 70
kR104 k < 0.01 kR83 k < 0.01 kR70 k < 0.01
1.6250 0.6560 0.2820
to compute an A-orthonormal matrix group M1 , M2 , . . . , Mp as the basis of iteration according to (3.1), and the algorithm yielded by it is proved to be convergent. Certainly, it is possible to define different A-orthonormal matrix group and design another algorithm to solve matrix equation (1.1). Acknowledgements The authors are very grateful to the referee for their valuable comments. They also thank Professor Michael Ng for his helpful suggestions. This research was supported by the Excellent Youth Foundation of Educational Committee of Hunan Province (09B113), by National natural Science Foundation of China (10571047), by Specialized Research Fund for the Doctoral Program of Higher Education (20060532014). References [1] [2] [3] [4] [5] [6] [7] [8]
G.H. Golub, C.F. Van Loan, Matrix Computations, Johns Hopkins UP, Baltimore, 1996, pp. 53–644. S.S. Zhou, H. Dai, Inverse Problems of Algebraic Eigenvalue, Henan Science and Technology Press, Zhengzhou, 1991, pp. 6–93. C.F. Van Loan, Generalizing the singular value decomposition, SIAM J. Numer. Anal. 13 (1976) 76–83. M.T. Chu, R.E. Funderlic, G.H. Golub, On a variational formulation of generalized singular value decomposition, SIAM J. Matrix Anal. Appl. 18 (1) (1997) 1082–1092. G.H. Golub, H.Y. Zha, Perturbation analysis of the canonical correlations of matrix pairs, Linear Algebra Appl. 210 (1994) 3–28. N.J. Higham, Computing a nearest symmetric positive semidefinite matrix, Linear Algebra Appl. 103 (1988) 103–118. N. Ujevic, A new iterative method for solving linear systems, Appl. Math. Comput. 179 (2006) 725–730. Y.F. Jing, T.Z. Huang, On a new iterative method for solving linear systems and comparison results, J. Comput. Appl. Math. 220 (2008) 74–84.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Successive projection iterative method for solving matrix equation AX = B Fan-Liang Li a,∗ , Li-sha Gong b , Xi-Yan Hu c , Lei Zhang c a
Institute of Mathematics and Physics, School of Sciences, Central South University of Forestry and Technology, Changsha, 410004, PR China
b
School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu, Sichuan, 610054, PR China
c
College of Mathematics and Econometrics, Hunan University, Changsha, 410082, PR China
article
info
Article history: Received 31 August 2009 Received in revised form 31 December 2009 MSC: 15A24 11Y40
abstract In this paper, we present a new iterative method (successive projection iterative method) to solve matrix equation AX = B, where A is a symmetric positive definite (SPD) matrix. Based on this method an algorithm is proposed and proved to be convergent. In addition, analysis of the algorithm and numerical experiments are also given to show the efficiency of the method. © 2010 Elsevier B.V. All rights reserved.
Keywords: Matrix equation Symmetric positive definite matrix A-orthogonal matrix group Successive projection iterative method
1. Introduction The problem of solving matrix equation is a hot topic in the field of numerical algebra in recent years. Many authors studied this problem by using the special structure of matrix and the decomposition of matrix pairs, and a series of meaningful results were achieved [1–6]. For example, Golub [1] and Zhou [2] discussed the solution of matrix equation with the singular value decomposition (SVD) of matrix; Loan [3] and Chu [4] put forward and discussed the generalized singular value decomposition (GSVD) of matrix pairs to solve matrix equation; Golub [5] and Higham [6] considered matrix equation with the canonical correlation decomposition of matrix pairs (CCD) and the quotient singular value decomposition (QSVD) of matrix pairs, respectively. The method in these papers is called direct method. With this method, the authors can present the solvability conditions and the expression of the solution for matrix equations directly. However, in these cases, the solvability conditions and the expression of the solution for matrix equations are usually very complex which may cause huge error during calculating. Furthermore, some matrix equations cannot be solved with direct method at all. For example, it is difficult to obtain the centrosymmetric solution, the reflexive solution and the double symmetric solution of matrix equation AXB = C with direct method. In this paper, based on the Gauss–Seidel iterative method and successive projection method, which are used to solve linear equations Ax = b [7,8], we propose an improved method to solve matrix equation AX = B. Let Rn×m be the set of all n × m real matrices, In be the identity matrix of order n, SRn+×n denote the set of all n × n SPD matrices, M T , M + , r (M ) and tr(M ) represent the transpose, the Moore–Penrose generalized inverse, rank and trace of M,
∗
Corresponding author. Tel.: +86 15973179263; fax: +86 073185010892. E-mail address: [email protected] (F.-L. Li).
0377-0427/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2010.03.008
2406
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
respectively. For M = (mij ), N = (nij ) ∈ Rn×m , hM , N i = tr(N T M ) denote the inner product of matrix M and N. The induced 1
1
matrix norm is called Frobenius norm, i.e. kM k = hM , M i 2 = (tr(M T M )) 2 , then Rn×m is a Hilbert inner product space. Definition 1. Let M , N ∈ Rn×m , if hM , N i = tr(N T M ) = 0, then M , N is called orthogonal; for a given SPD matrix A, if hM , N iA = tr(N T AM ) = 0, then M , N is called A-orthogonal. Definition 2. Let A ∈ SRn+×n , if M1 , M2 , . . . , Mp ∈ Rn×m , and Mi 6= 0, i = 1, 2, . . . , p, for any Mi , Mj , hMi , Mj i = tr(MjT Mi ) =
0(hMi , Mj iA = tr(MjT AMi ) = 0), i 6= j, i, j = 1, 2, . . . , p, then M1 , M2 , . . . , Mp are called orthogonal (A-orthogonal) matrix group. It is obviously that orthogonal (A-orthogonal) matrix group M1 , M2 , . . . , Mp are linearly independent. That is, if c1 M1 + c2 M2 + · · · + cp Mp = 0 holds, then c1 = c2 = · · · = cp = 0. Definition 3. Let A ∈ SRn+×n , if M1 , M2 , . . . , Mp ∈ Rn×m are orthogonal (A-orthogonal) matrix group, and kMi k = 1(kMi kA = tr(MiT AMi )
12
= 1), i = 1, 2, . . . , p, then M1 , M2 , . . . , Mp are called orthonormal (A-orthonormal) matrix group.
In this paper, we consider the following problem. Problem. Giving A ∈ SRn+×n , B ∈ Rn×m , find X ∈ Rn×m such that AX = B.
(1.1)
This paper is organized as follows. In Section 2, successive projection iterative method to solve (1.1) is presented. In Section 3, an algorithm is derived by giving an A-orthonormal matrix group and proved to be convergent. In the end, analysis of the algorithm and numerical experiments are given. 2. Successive projection iterative method Let A ∈ SRn+×n , if M1 , M2 , . . . , Mp ∈ Rn×m is an A-orthonormal matrix group, and let L = K = span{M1 , M2 , . . . , Mp }, p ≤ m. Then the method of successive projection iterative can be described as follows: let Xk , Xk+1 be the k, k + 1 step iterative approximate solution of (1.1), respectively, then Xk+1 satisfies the following conditions. Xk+1 ∈ Xk + K ,
Rk+1 = B − AXk+1 ⊥ L.
From this, we can suppose that (k)
(k)
Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp(k) Mp , (k)
(k)
(k)
where α1 , α2 , . . . , αp are constants. Hence, the k + 1 step residual matrix (k)
(k)
Rk+1 = B − AXk+1 = B − A(Xk + α1 M1 + α2 M2 + · · · + αp(k) Mp )
= Rk − α1(k) AM1 − α2(k) AM2 − · · · − αp(k) AMp , where Rk = B − AXk . It is easy to see hRk+1 , Mi i = 0, i = 1, 2, . . . , p from Rk+1 = B − AXk+1 ⊥ L, and L = K = span{M1 , M2 , . . . , Mp }. Which implies
(k) α1 hAM1 , M1 i + α2(k) hAM2 , M1 i + · · · + αp(k) hAMp , M1 i = hRk , M1 i, α (k) hAM1 , M2 i + α (k) hAM2 , M2 i + · · · + α (k) hAMp , M2 i = hRk , M2 i, p 1 2 .. . (k) α1 hAM1 , Mp i + α2(k) hAM2 , Mp i + · · · + αp(k) hAMp , Mp i = hRk , Mp i. Noticing that M1 , M2 , . . . , Mp ∈ Rn×m is an A-orthonormal matrix group, we have
α1(k) = tr(M1T Rk ), α2(k) = tr(M2T Rk ), . . . , αp(k) = tr(MpT Rk ). Up to now, we can design the following method to solve (1.1). Successive projection iterative method. 1. Choose an initial matrix X1 ∈ Rn×m , compute R1 = B − AX1 , k := 1. 2. If kRk k ≤ error bound, then X = Xk , stop iterating; otherwise go to 3.
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
(k)
2407
(k)
(k)
3. Give an A-orthonormal matrix group M1 , M2 , . . . , Mp ∈ Rn×m , compute Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp Mp , (k)
where αi = tr( ), i = 1, 2, . . . , p, p ≤ m. 4. Let k := k + 1, and go to 2. MiT Rk
Theorem 1. In the above method, if Xk , Xk+1 are the k, k + 1 step iterative approximate solution of (1.1) respectively, and X is the exact solution of (1.1). Then
kDk kA ≥ kDk+1 kA , where Dk = X − Xk , Dk+1 = X − Xk+1 . (In fact, Dk is the difference between the exact solution and k step iterative approximate solution, and we call it the k step error.). (k)
(k)
(k)
Proof. According to Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp Mp , we have (k)
(k)
Dk+1 = X − Xk+1 = X − (Xk + α1 M1 + α2 M2 + · · · + αp(k) Mp )
= Dk − (α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp ). This implies
hADk+1 , Dk+1 i = hADk − (α1(k) AM1 + α2(k) AM2 + · · · + αp(k) AMp ), Dk − (α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp )i = hADk , Dk i − 2hADk , α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp i +
p X (αi(k) )2 . i=1
Combining ADk = AX − AXk = B − AXk = Rk , it is easy to see
hADk+1 , Dk+1 i = hADk , Dk i − 2hRk , α1(k) M1 + α2(k) M2 + · · · + αp(k) Mp i +
p X (αi(k) )2 . i =1
(k)
Noticing that αi
= hRk , Mi i, i = 1, 2, . . . , p, we have
hADk+1 , Dk+1 i = hADk , Dk i −
p X (αi(k) )2 , i =1
i.e.
kDk k2A − kDk+1 k2A =
p X (αi(k) )2 ≥ 0,
(2.1)
i=1
and this gives the conclusion.
(αi(k) )2 = 0, i.e. αi(k) = 0, i = 1, . . . , p, then Xk+1 = Xk , and we do not have any reduction Pp (k) 2 2 2 of k + 1 step error comparing to k step error in A-norm meaning. But if i=1 (αi ) > 0, i.e. kDk kA > kDk+1 kA , then the Theorem 1 shows that if
Pp
i=1
error is reduced. Theorem 1 also shows that the sequence {kDk k2A } is convergent. Hence, we can suppose as follows.
kDk k2A −→ a (k −→ ∞),
(2.2)
where a ≥ 0. Combining (2.1) and (2.2) gives the following results.
αi(k) −→ 0 (k −→ ∞).
(2.3)
From the theory of projection, it is easy to derive the following theorem. Theorem 2. In the above method, if Xk , Xk+1 are the k, k + 1 step iterative approximate solution of (1.1), respectively, then for any Xk˜+1 ∈ Xk + K , Xk+1 satisfies
kB − AXk+1 k =
min ∀Xk˜+1 ∈Xk +K
kB − AXk˜+1 k.
Theorem 3. In the above method, if Xk , Xk+1 are the k, k + 1 step iterative approximate solution of (1.1), respectively, then for any Xk˜+1 ∈ Xk + K , Xk+1 satisfies
kX − Xk+1 kA =
min ∀Xk˜+1 ∈Xk +K
kX − Xk˜+1 kA ,
where X is the exact solution of (1.1).
2408
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
Proof. Definition 3 shows that kX − Xk˜+1 kA = tr[(X − Xk˜+1 )T A(X − Xk˜+1 )], and
kB − AXk˜+1 k = tr[(B − AXk˜+1 )T (B − AXk˜+1 )] = tr[(X − Xk˜+1 )T AT A(X − Xk˜+1 )] = tr[(X − Xk˜+1 )T A2 (X − Xk˜+1 )]. A ∈ SRn+×n implies that there exists an n × n orthogonal matrix P such that
λ1
0 A=P .. .
λ2 .. .
··· ··· .. .
0
0
···
0
0 0
T .. P , . λn
where λi > 0, i = 1, 2, . . . , n, are eigenvalues of matrix A. Thus A2 can be written as
2 λ1 0 A2 = P . ..
λ22 .. .
··· ··· .. .
0
0
···
0
0 0
T .. P . . λ2n
This implies
λ1
0 ˜ T kX − Xk˜+1 kA = tr (X − Xk+1 ) P .. .
λ2 .. .
··· ··· .. .
0
0
···
λ 0 kB − AXk˜+1 k = tr (X − Xk˜+1 )T P . ..
0
λ22 .. .
0
0
··· ··· .. . ···
2 1
0
0 0
T ˜ .. P (X − Xk+1 ) , . λn 0 0
T .. P (X − Xk˜+1 ) . . 2 λn
Let F (Xk˜+1 ) = kX − Xk˜+1 kA , G(Xk˜+1 ) = kB − AXk˜+1 k, then it can be verified without difficulty that function F (Xk˜+1 ) and G(Xk˜+1 ) have the same point of minimum. Combining Theorem 2, we can obtain the conclusion of Theorem 3. Theorem 3 shows that Xk+1 is the result of projection of Xk onto K (orthogonal to L) iff it minimizes the A-norm of the k + 1 step error over Xk + K . Theorem 3 also shows that if K1 = L1 ⊆ K2 = L2 , then the reduction of A-norm of the k step error obtained by subspaces K2 = L2 is more or equal to that of subspaces K1 = L1 . Hence by increasing the value of p, the convergence rate may increase. 3. An algorithm and numerical results Analysis of the above method Suppose A ∈ SRn+×n , the key problem of the above method is how to obtain an A-orthonormal matrix group M1 , M2 , . . . , Mp ∈ Rn×m in step 3. In the following section, we give a method to obtain a special A-orthonormal matrix (k)
(k)
(k)
(k)
(k)
group M1 , M2 , . . . , Mp ∈ Rn×m . Let Rk = (rij )n×m , and suppose |ri1 j1 |, |ri2 j2 |, . . . , |rip jp | be the p largest values in all |rij | and (k)
(k)
(k)
satisfy |ri1 j1 | ≥ |ri2 j2 | ≥ · · · ≥ |rip jp |, jl 6= jh , jl , jh = j1 , . . . , jp . Let M1 = √
1
1 1 Ei1 j1 , M2 = √ Ei2 j2 , . . . , Mp = √ Eip jp , ai1 i1 ai 2 i 2 aip ip
(3.1)
where Eij is an n × m matrix with the (i, j) element equal to 1, and others equal to zero, ai1 i1 , ai2 i2 , . . . , aip ip are the i1 th, i2 th, . . . , ip th element of A in the main diagonal line. It is can be obtained immediately that M1 , M2 , . . . , Mp is an A-orthonormal matrix group with p ≤ m, and (k)
ri j
αl(k) = tr(MlT Rk ) = √ l l , ail il
l = 1 , 2 , . . . , p.
Combining (2.3) gives Rk −→ 0 (k −→ ∞), which implies Xk −→ X
(k −→ ∞),
(3.2)
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
2409
where, X is the exact solution of (1.1). The above analysis shows that if the A-orthonormal matrix group M1 , M2 , . . . , Mp is given by (3.1), then we can provide the following algorithm which ensures that iterative sequence {Xk } converge to X . Algorithm. 1. Input an initial matrix X1 ∈ Rn×m , compute R1 = B − AX1 , k := 1. 2. If kRk k ≤ error bound, then X = Xk , stop iterating; otherwise go to 3. 3. Compute M1 , M2 , . . . , Mp according to (3.1). (k)
(k)
(k)
(k)
4. Compute Xk+1 = Xk + α1 M1 + α2 M2 + · · · + αp Mp , where αi 5. Let k := k + 1, and go to 2.
= tr(MiT Rk ), i = 1, 2, . . . , p, p ≤ m.
Analysis of the algorithm In the above algorithm, step 1 requires n2 m arithmetic operations (additions and multiplications), step 2 requires n2 m arithmetic operations (additions and multiplications), step 3 requires n arithmetic operations (multiplications), and step 4 requires nmp arithmetic operations (additions and multiplications). Numerical examples (n = 10, m = 6) Let
40 10 0.5 0.5 0.5 A= 0.5 0.5 0.5 0.5 0.5
0.5 10 40 10 0.5 0.5 0.5 0.5 0.5 0.5
10 40 10 0.5 0.5 0.5 0.5 0.5 0.5 0.5
54 63.5 63.5 63.5 63.5 B= 63.5 63.5 63.5 63.5
54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
54
0.5 0.5 10 40 10 0.5 0.5 0.5 0.5 0.5
54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
0.5 0.5 0.5 10 40 10 0.5 0.5 0.5 0.5 54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
0.5 0.5 0.5 0.5 10 40 10 0.5 0.5 0.5
0.5 0.5 0.5 0.5 0.5 10 40 10 0.5 0.5
54 63.5 63.5 63.5 63.5 63.5 63.5 63.5 63.5 54
0.5 0.5 0.5 0.5 0.5 0.5 10 40 10 0.5
0.5 0.5 0.5 0.5 0.5 0.5 0.5 10 40 10
0.5 0.5 0.5 0.5 0.5 , 0.5 0.5 0.5 10 40
54 63.5 63.5 63.5 63.5 , 63.5 63.5 63.5 63.5 54
it is clear that the exact solution X of (1.1) is
1 1 1 X = 1 T
1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 . 1 1 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 , 0 0 0
Input an initial guess X1
0 X1T
0 0 = 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
and use kRk = B − AXk k < 0.01 as the stopping criterion. By running some MATLAB codes, we obtain the following table with p = 4, 5, 6. Table 1 confirms that the above algorithm is convergent and by increasing the value of p, the convergence rate may increase. Conclusion In this paper, based on the Gauss–Seidel iterative method and successive projection method, which are used to solve linear equation Ax = b in [7,8], we propose a successive projection iterative method to solve matrix equation AX = B. Compared to [7,8], this paper has two important achievements. One is that the successive projection method is used to solve matrix equation AX = B in this paper instead of linear equation Ax = b in [7,8]. The other is we provide the method
2410
F.-L. Li et al. / Journal of Computational and Applied Mathematics 234 (2010) 2405–2410
Table 1 Numerical results for example. p
Step
kRk k
Cputime
4 5 6
104 83 70
kR104 k < 0.01 kR83 k < 0.01 kR70 k < 0.01
1.6250 0.6560 0.2820
to compute an A-orthonormal matrix group M1 , M2 , . . . , Mp as the basis of iteration according to (3.1), and the algorithm yielded by it is proved to be convergent. Certainly, it is possible to define different A-orthonormal matrix group and design another algorithm to solve matrix equation (1.1). Acknowledgements The authors are very grateful to the referee for their valuable comments. They also thank Professor Michael Ng for his helpful suggestions. This research was supported by the Excellent Youth Foundation of Educational Committee of Hunan Province (09B113), by National natural Science Foundation of China (10571047), by Specialized Research Fund for the Doctoral Program of Higher Education (20060532014). References [1] [2] [3] [4] [5] [6] [7] [8]
G.H. Golub, C.F. Van Loan, Matrix Computations, Johns Hopkins UP, Baltimore, 1996, pp. 53–644. S.S. Zhou, H. Dai, Inverse Problems of Algebraic Eigenvalue, Henan Science and Technology Press, Zhengzhou, 1991, pp. 6–93. C.F. Van Loan, Generalizing the singular value decomposition, SIAM J. Numer. Anal. 13 (1976) 76–83. M.T. Chu, R.E. Funderlic, G.H. Golub, On a variational formulation of generalized singular value decomposition, SIAM J. Matrix Anal. Appl. 18 (1) (1997) 1082–1092. G.H. Golub, H.Y. Zha, Perturbation analysis of the canonical correlations of matrix pairs, Linear Algebra Appl. 210 (1994) 3–28. N.J. Higham, Computing a nearest symmetric positive semidefinite matrix, Linear Algebra Appl. 103 (1988) 103–118. N. Ujevic, A new iterative method for solving linear systems, Appl. Math. Comput. 179 (2006) 725–730. Y.F. Jing, T.Z. Huang, On a new iterative method for solving linear systems and comparison results, J. Comput. Appl. Math. 220 (2008) 74–84.