- Email: [email protected]

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A family of iterative methods for computing the approximate inverse of a square matrix and inner inverse of a non-square matrix Weiguo Li *, Zhi Li School of Mathematics and Computational Sciences, China University of Petroleum, Dongying 257061, Shandong Province, PR China

a r t i c l e

i n f o

Keywords: Inverse of matrix Inner inverse of matrix Iterative method Convergence Linear systems

a b s t r a c t Based on a quadratical convergence method, a family of iterative methods to compute the approximate inverse of square matrix are presented. The theoretical proofs and numerical experiments show that these iterative methods are very effective. And, more importantly, these methods can be used to compute the inner inverse and their convergence proofs are given by fundamental matrix tools. Ó 2009 Elsevier Inc. All rights reserved.

1. Introduction It is well-known that matrix inverse and generalized inverse are important in applied ﬁelds of nature science, such as solution to various systems of linear equation, eigenvalue problems, the linear least square problems. Of course, many important methods to ﬁnd inverse and generalized inverse of matrix have been developed . In these methods, direct methods usually need much cost in both time and space in order to achieve the desirable results, sometimes it may not be able to work at all. So iterative methods are often effective especially for large scale systems with sparse matrix. In this paper, we consider a family of methods for computing approximate inverse of nonsingular matrix and generalized inverse of any m n matrix. In [1], the authors mentioned an iterative formula from [2], that is

V qþ1 ¼ V q ð2I AV q Þ

ð1:1Þ

and showed

kEqþ1 k 6 kAkkEq k2 ; where Eq ¼ I AV q , I is an identity matrix with the same dimension as matrix A, and kI AV 0 k < 1. In [3], the author gave some comments on convergence analysis in this iterative algorithm for the inverse of an nonsingular square matrix considered in paper [1]. Motivated by Saberi Najaﬁ et al. [1], Wu [3] and SMS algorithm for computing Drazin inverse in [4] and outer inverse in [5], we can get a high order method to ﬁnd the inverse of nonsingular square matrix and the inner inverse of any m n matrix. Indeed, there are many current and related references, some of which are interesting and beneﬁcial to our paper: computing Moore-Penrose inverses of matrices by Newton’s iteration [6] or by interval iterative methods [7], computing the group inverses of singular Toeplitz matrices by Modiﬁed Newton’s algorithm [8], computing generalized matrix inversion and rank by successive matrix powering method [9], approximating outer generalized inverse [10], Drazin inð2Þ verse [11] and the generalized inverse AT;S of a matrix or an operator A [12,13] by constructive methods. * Corresponding author. E-mail address: [email protected] (W. Li). 0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.10.038

3434

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

2. A family of iterative methods for computing approximate inverse Let A be nonsingular and V 0 as an initial approximation to A1 . With (1.1), we get a sequence of approximation fV q g to A . In [3], the author proved that V q generated from the iterative formula (1.1) is quadratically convergent to A1 . In this section, we ﬁrst introduce a third-order convergence iterative method for ﬁnding the inverse of a square matrix as follows: 1

V qþ1 ¼ V q ð3I 3AV q þ ðAV q Þ2 Þ:

ð2:1Þ

Now we have Theorem 2.1. If A is nonsingular and initial approximation V 0 satisﬁes

kE0 k ¼ kI AV 0 k < 1; then the iterative formula (2.1) is convergent and V q converges cubically to A1 . Proof. From (2.1) and let

Eq ¼ I AV q ; then

Eqþ1 ¼ I AV qþ1 ¼ I AV q ð3I 3AV q þ ðAV q Þ2 Þ ¼ ðI AV q Þ3 ¼ E3q :

ð2:2Þ

Since

kE0 k ¼ kI AV 0 k < 1: So q

kEq k 6 kEq1 k2 6 kE0 k3 ! 0; when q ! 1: Namely,

I AV q ! 0; when q ! 1; that is,

lim V q ¼ A1 :

q!1

Let

eq ¼ A1 V q ; then

Aeq ¼ I AV q ¼ Eq : From (2.2), we have

ðAeq ÞðAeq Þ2 ¼ E3q ¼ Eqþ1 : By Aeqþ1 ¼ Eqþ1 , we obtain

eqþ1 ¼ eq ðAeq Þ2 : Therefore it follows immediately that

keqþ1 k 6 keq ðAeq Þ2 k 6 kAk2 keq k3 : Consequently, it is proved that the iterative formula (2.1) is convergent and V q at least converges cubically to A1 . h Numerical experiments show the results of formula (1.1) are evidently improved by (2.1). See Examples 4.1 and 4.2. In fact, under the same conditions of Theorem 2.1, we can present a family of formulae as follows:

kðk 1Þ V qþ1 ¼ V q kI AV q þ þ ð1Þk1 ðAV q Þk1 ; k ¼ 2; 3; . . . : 2 We can obtain the following conclusion that is similar to Theorem 2.1. Theorem 2.2. If A is nonsingular and initial approximation V 0 satisﬁes

kE0 k ¼ kI AV 0 k < 1;

ð2:3Þ

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

3435

then V q generated from iterative formula (2.3) is convergent to A1 and the order of convergence of this formula is k. Proof. From (2.1) and let

Eq ¼ I AV q ; then

kðk 1Þ Eqþ1 ¼ I AV qþ1 ¼ I AV q kI AV q þ þ ð1Þk1 ðAV q Þk1 ¼ ðI AV q Þk ¼ Ekq : 2

ð2:4Þ

The rest of this proof is similar to that of Theorem 2.1. h Especially, when k ¼ 2, the formula (2.3) is just the formula (1.1) from [2], and when k ¼ 3, the formula (2.3) is just the formula (2.1). In fact, the iterative methods to ﬁnd the approximate inverse of nonsingular matrix are very powerful in practice but a difﬁcult problem is to ﬁnd a good initial approximate inverse besides the convergence of iterative sequence. In the following theorem, we will give an easy method to ﬁnd an initial approximate inverse of any nonsingular matrix. Theorem 2.3. For any nonsingular matrix A 2 C nn and any 0 < a < r22ðAÞ, where r1 ðAÞ ¼ kAk2 , and setting V 0 ¼ aA , where A is 1 the conjugate transpose matrix of A, then V q generated from the iterative formula (2.3) is convergent to A1 with k-order convergence. Proof. Since A is nonsingular, there exists an unitary matrix Q such that

Q A AQ ¼ diagðr21 ; . . . ; r2n Þ; where

r21 P P r2n > 0. Based on the condition 0 < a < r22ðAÞ, it is easy to see that 1

1 < 1 ar21 6 kj ðI aA AÞ 6 1 ar2n < 1: Let

d ¼ max j1 ar21 j; 1 ar2n ; then qðI aA AÞ ¼ d < 1, where qðI aA AÞ is the spectral radius of matrix I aA A. By d ¼ qðI V 0 AÞ < 1 and Lemma 5.6.10 in [14], there exists a constant > 0 and a matrix norm k k such that

kI V 0 Ak 6 qðI V 0 AÞ þ ¼ d þ < 1: Based on Theorem 2.2, we have gained our conclusion. h

3. Computing approximate inner (generalized) inverse of non-square matrix Much more numerical experiments show the family of formulae (2.3) can be used to compute approximate inner (generalized) inverse of a matrix. The following is a detailed deduction of this statement. Lemma 3.1 [15]. Let G 2 C nn . Then lim Gr exists if and only if the following three conditions hold at the same time: r!1

(1) qðGÞ 6 1, where qðGÞ is the spectral radius of matrix G; (2) if qðGÞ ¼ 1, then these eigenvalues whose modulus are 1 must be 1; (3) if qðGÞ ¼ 1, then these eigenvalues whose modulus are 1 must be semi-simple (i.e., corresponding Jordan block is 1 1). q

Lemma 3.2. Let A 2 C mn ; fV q g # C nm . Then lim ðI V q AÞ ¼ lim ðI V 0 AÞk exists if and only if the following two conditions hold. q!1

q!1

(1) j1 lj ðV 0 AÞj 6 1; j ¼ 1; . . . ; n; (2) if j1 lj ðV 0 AÞj ¼ 1, then there must be lj ðV 0 AÞ ¼ 0, and these eigenvalues are semi-simple, where I is identity matrix, V q comes from (2.3) and lj ðV 0 AÞ; j ¼ 1; . . . ; n, are eigenvalues of V 0 A. q

Proof. Let G ¼ I V 0 A and r ¼ k . From Lemma 3.1 and kj ðI V 0 AÞ ¼ 1 lj ðV 0 AÞ, j ¼ 1; . . . ; n, the conclusion is clear. h Note 3.1. The conditions of Lemma 3.2 are equivalent to that the eigenvalues of ðV 0 AÞ are in the circle on the complex plane, whose center is at ð1; 0Þ and whose radius is 1. Or the eigenvalues of ðV 0 AÞ are in the origin (these eigenvalues are semi-simple). See Fig. 3.1. In this situation,

3436

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

1.5

Range of Eigenvalues of VA 1

imaginary axis

0.5

0

0

1

2

−0.5

−1

−1.5 −0.5

0

0.5

1

1.5

2

2.5

real number axis Fig. 3.1. The eigenvalues of V 0 A are in the circle or at origin on the complex plane (Note 3.1). q

lim ðI V q AÞ ¼ lim ðI V 0 AÞk ¼ PBP 1 ;

q!1

q!1

ð3:1Þ

where P is a similar matrix of Jordan canonical form of matrix V 0 A; r ¼ rankðV 0 AÞ,

P1 ðV 0 AÞP ¼

0nr

0

0

Jr

and

B¼

Inr

0

0

0

:

Deﬁnition 3.1 [16]. Let A 2 C mn . If there exists a matrix V 2 C nm satisfying AVA ¼ A, then V is called the inner (generalized) inverse of A or {1} inverse of A. Note 3.2. Inner inverse is not unique. In general textbooks, the set of the inner inverse of the matrix A is denoted as A . Let A 2 C mn and rankðAÞ ¼ r. If there exist two invertible matrices P 2 C nn and R 2 C mm such that

RAP ¼

Ir

0

0

0

;

then V 2 A if and only if

V ¼P where U 2 C

Ir

U

W

T

rðmrÞ

R;

; W 2 C ðnrÞr ; T 2 C ðnrÞðmrÞ are arbitrary.

Lemma 3.3 [16]. If V 2 A , then x ¼ Vb is a special solution of compatible systems of linear equations Ax ¼ b. Consider the sequence fV q g generated from (2.3). We can choose V 0 ¼ aA; 0 < a < r22ðAÞ, where r1 ðAÞ is the maximum singular 1 value of matrix A. In the following, we will prove that this sequence satisﬁes lim AV q A ¼ A for any A 2 C mn and any k. That is, V k is q!1 an approximate inner inverse of matrix A. Lemma 3.4. For any A 2 C mn and any X 2 C nl ; AX ¼ 0 if and only if A AX ¼ 0. Proof. If AX ¼ 0, then A AX ¼ 0. Inversely, if A AX ¼ 0, then A AX j ¼ 0; j ¼ 1; . . . ; l, where X ¼ ½X 1 ; . . . ; X l . Therefore, X j A AX j ¼ 0; j ¼ 1; . . . ; l, that is, ðAX j Þ ðAX j Þ ¼ 0, or AX j ¼ 0; j ¼ 1; . . . ; l. So we get AX ¼ 0. h Now we prove the main theorem of this section.

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

3437

Theorem 3.1. For any A 2 C mn and any 0 < a < r22ðAÞ, let V 0 ¼ aA , then the sequence (2.3) satisﬁes 1

lim AV q A ¼ A

ð3:2Þ

q!1

and the order of convergence is k. Proof. Since 0 < a < r22ðAÞ, the eigenvalues of matrix aA A satisﬁes 1

0 6 lj ðaA AÞ < 2; and

j ¼ 1; . . . ; n;

lj ðaA AÞ are all semi-simple. Obviously, there exists an unitary matrix Q such that Q A AQ ¼

0nr

0

0

Kr

;

ð3:3Þ

where Kr ¼ diagðr21 ; . . . ; r2r Þ, and r2j > 0; j ¼ 1; . . . ; r are singular values of matrix A. Then the sequence (2.3) satisﬁes the conditions of Lemma 3.2. So we know that q

lim ðI V q AÞ ¼ lim ðI aA AÞk ¼ QBQ ¼ Q

q!1

Inr

0

0

0;

q!1

exists. Therefore, we have

lim ðAV q A AÞ ¼ lim AðI V q AÞ ¼ AQ

q!1

q!1

Inr 0

Q

0 Q : 0:

ð3:4Þ

By (3.3),

A AQ

Inr

0

0

0

¼Q

0nr

0

0

Kr

Inr

0

0

0

¼ 0:

ð3:5Þ

Again from Lemma 3.4,

AQ

Inr

0

0

0:

¼ 0:

ð3:6Þ

Therefore, from (3.4) and (3.6), we obtain lim ðAV q AÞ ¼ A: q!1

proof, we only need to prove that ðI V q AÞ In order to prove that AV q A has k-order convergence, from above converge kq ðI a K Þ converges to zero matrix with k-order. By converges to Q BQ with k-order, or to prove that r r q 0 < ar2j < 2; 1 < dj ¼ 1 ar2j < 1; j ¼ 1; . . . ; r, and ðdj Þk ! 0 when0 q ! 1, it is self-evident. h Note 3.3. From (3.2), we can say that V k is the approximation inner inverse of A. But we cannot prove that V k converges (in norm) to a ﬁxed inner inverse of A. (In fact, V k is weakly convergent to an inner inverse of A.) A lot of numerical experiments can show this fact.

4. Numerical experiments We ﬁrst present two examples to illustrate the efﬁciency of the new iterative method with the new initial approximation method to compute the approximate inverse of a nonsingular square matrix. We compare the results of formulae (2.3) when k ¼ 2 and k ¼ 3. In the following examples, we all set that termination parameter e ¼ 108 and maximum iterative number of step n ¼ 1000. Example 4.1. Let A ¼ randð100; 100Þ; a ¼

1 kAk22

and V 0 ¼ aA . The stop criterion is kI VAk2 < e. We have tested 50 times with

MATLAB 7 and carried out with double precision arithmetic on the computer pentium 4. The time consuming (s) and the iterative times are compared in Figs. 4.1 and 4.2, respectively. x-axis represents the times of choice random matrices, y-axis in Fig. 4.1 represents the computing time and y-axis in Fig. 4.2 represents the iterative times. The dashed line corresponds to the iteration by (1.1) and solid line to the iteration (2.1). We see that (2.1) has a more considerable improvement than (1.1) in both the computing time and the iterative times. 1 Example 4.2. Let A ¼ randnð300; 300Þ; a ¼ kAk 2 and V 0 ¼ aA . The stop criterion is kI VAk2 < e. We have tested 20 times 2

with MATLAB 7 and carried out with double precision arithmetic on the computer pentium 4. The time consuming (s) and the iterative times are compared in Figs. 4.3 and 4.4, respectively. From Figs. 4.3 and 4.4, we can get the same conclusion as that in Example 4.1. Example 4.3. We will compare iterative methods (1.1) and (2.1) for computing the inner inverse of non-square matrix A. We also consider the following compatible systems of linear algebraic equations

3438

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

0.5 2−order convergence 3−order convergence

iterative time (sencond)

0.45

0.4

0.35

0.3

0.25

0.2 0

5

10

15

20

25

30

35

40

45

50

random times Fig. 4.1. Comparison of consuming time (Example 4.1).

40 2−order convergence 3−order convergence

iterative times

35

30

25

20

15

0

5

10

15

20

25

30

35

40

45

50

random times Fig. 4.2. Comparison of iterative times (Example 4.1).

Ax ¼ b; A ¼ ðaij Þ;

i ¼ 1; . . . ; m;

j ¼ 1; . . . ; n;

ð4:1Þ

where

bi ¼

m X

ail ;

i ¼ 1; 2; . . . ; n:

ð4:2Þ

l¼1

It is easy to see that the accuracy solution of the system (4.1) is given by

x ¼ x ¼ ð1; 1; . . . ; 1ÞT :

randnð200; 210Þ 1 ; a ¼ kAk 2 and V 0 ¼ aA . The stop criterion is kA AVAk2 < e. We have tested 50 times with onesð40; 210Þ 2 MATLAB 7 and carried out with double precision arithmetic on the computer pentium 4. The time consuming (s), the iterative times and the absolute error between the approximate solution V q B and the accuracy solution x of the compatible systems are compared in Figs. 4.5, 4.6, and 4.7, respectively. x-axis represents the times of choice random matrices, y-axis in Let A ¼

3439

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

7 2−oredr convergence 3−order convergence

6.5

iterative time (sencond)

6

5.5

5

4.5

4

3.5

3

2.5 0

2

4

6

8

10

12

14

16

18

20

random times Fig. 4.3. Comparison of consuming time (Example 4.2).

35 2−oredr convergence 3−order convergence

iterative times

30

25

20

15

10 0

2

4

6

8

10

12

14

16

18

20

random times Fig. 4.4. Comparison of iterative times (Example 4.2).

Fig. 4.5 represents the computing time, y-axis in Fig. 4.6 represents the iterative times and y-axis in Fig. 4.7 represents the absolute error between the approximate solution V q B and the accuracy solution x of the compatible systems. The dashed line corresponds to the iteration by (1.1) and the solid line to the iteration (2.1). We see that the iteration method (2.1) has a more considerable improvement than (1.1) in all the computing time and the iterative times.

3440

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

1.9 2−order convergence 3−order convergence

1.8

iterative time (sencond)

1.7

1.6

1.5

1.4

1.3

1.2

1.1

1

0

5

10

15

20

25

30

35

40

45

50

random times Fig. 4.5. Comparison of consuming time (Example 4.3).

24 2−order convergence 3−order convergence

22

iterative times

20

18

16

14

12

10

8 0

5

10

15

20

25

30

35

40

45

50

random times Fig. 4.6. Comparison of iterative times (Example 4.3).

Example 4.4. As an important application, we consider the typical ill-conditioned systems of linear algebraic equations with 1 ; i; j ¼ 1; 2; . . . ; n. coefﬁcient matrix (Hilbert) A ¼ ðhij Þnn , where hij ¼ iþj1

Ax ¼ b; where

ð4:3Þ

3441

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442 −11

10

error between approximate and accuracy solution

2−order convergence 3−order convergence

−12

10

−13

10

−14

10

0

5

10

15

20

25

30

35

40

45

50

random times Fig. 4.7. Absolute error (Example 4.3).

Table 4.1 Range of Eigenvalues of VA. n/cond(A)

Consuming time (s)

Computing times

Absolute error of solution

By (1.1)

By (2.1)

By (1.1)

By (2.1)

By (1.1)

10=1:6 1013

0.015

0.011

49

31

7:14 105

7:20 105

18

0.112

0.079

50

32

100=1:5 1020

0.623

0.441

53

34

8:42 10 0.0015

4

7:35 104 0.0018

200=1:1 1020

3.29

2.48

52

33

0.0034

0.0033

300=8:5 1019

11.88

8.29

53

34

0.0047

0.0035

500=1:9 1020

45.38

31.49

53

34

0.0081

0.0073

50=9:4 10

bi ¼

n X

0:01 l hil ;

i ¼ 1; 2; . . . ; n:

By (2.1)

ð4:4Þ

l¼1

It is easy to see that the accuracy solution of the system (4.3) is given by

x ¼ x ¼ ð0:01; 0:02; . . . ; 0:01 nÞT : Letting n ¼ 10; 50; 100; 200; 300, and 500, the corresponding spectral condition numbers from MATLAB function condðAÞ and the absolute errors between the approximate solution V q b and the accuracy solution x of the compatible systems are shown in Table 4.1. 1 We compute the inverse of Hilbert matrices by the iterative methods with a ¼ kAk 2 and V 0 ¼ aA . The stop criteria is 2 kA AVAk2 < e. We have tested with MATLAB 7 and carried out with double precision arithmetic on the computer pentium 4. The time consuming (s), the iterative times and the absolute error between the approximate solution and the accuracy solution are compared in Table 4.1.

5. Conclusions In this paper, we have developed a family of k-order iterative convergence formulae (2.3), which has generalized the quadratically convergent method in [1–3]. Especially, a three-order method is tested. The numerical experiments show the efﬁciency and improvement to second-order methods. Furthermore, these formulae can be used to compute the inner

3442

W. Li, Z. Li / Applied Mathematics and Computation 215 (2010) 3433–3442

(generalized) inverse of any matrix. At the same time, we also have given an easy method to furnish an initial approximate inverse which is very effective. Acknowledgement This research is supported by National Natural Science Foundation of China (60971132). References [1] H. Saberi Najaﬁ, M. Shams Solary, Computational algorithms for computing the inverse of a square matrix, quasi-inverse of a nonsquare matrix and block matrices, Applied Mathematics and Computation 183 (2006) 539–550. [2] G.M. Phillips, P.J. Taylor, Theory and Applications of Numerical Analysis, Academic Press, 1980. [3] Xinyuan Wu, A note on computational algorithm for the inverse of a square matrix, Applied Mathematics and Computation 187 (2007) 962–964. [4] Y. Wei, Successive matrix squaring algorithm for computing Drazin inverse, Applied Mathematics and Computation 108 (2000) 67–75. [5] Predrag S. Stanimirovic, Dragana S. Cvetkovic-Ilic, Successive matrix squaring algorithm for computing outer inverses, Applied Mathematics and Computation 203 (2008) 19–29. [6] Yimin Wei, Jianfeng Cai, Michael K. Ng, Computing Moore-Penrose inverses of Toeplitz matrices by Newton’s iteration, Mathematical and Computer Modelling 40 (1–2) (2004) 181–191. [7] Xian Zhang, Jianfeng Cai, Yimin Wei, Interval iterative methods for computing Moore-Penrose inverse, Applied Mathematics and Computation 183 (1) (2006) 522–532. [8] Jian-feng Cai, Michael K. Ng, Yi-min Wei, Modiﬁed Newton’s algorithm for computing the group inverses of singular Toeplitz matrices, Journal of Computational Mathematics 24 (5) (2006) 647–656. [9] L. Chen, E.V. Krishnamurthy, I. Macleod, Generalized matrix inversion and rank computation by successive matrix powering, Parallel Computing 20 (1994) 297–311. [10] D.S. Djordjevic, P.S. Stanimirovic, Y. Wei, The representation and approximation of outer generalized inverses, Acta Mathematica Hungar 104 (2004) 1– 26. [11] Y. Wei, H. Wu, The representation and approximation for Drazin inverse, Journal of Computational and Applied Mathematics 126 (2000) 417–423. [12] Y. Wei, A characterization and representation for the generalized inverse A2T;S and its applications, Linear Algebra and its Applications 280 (1998) 87– 96. [13] Yaoming Yu, Yimin Wei, The representation and computational procedures for the generalized inverse Að2ÞT; S of an operator A in Hilbert spaces, Numerical Functional Analysis and Optimization 30 (1–2) (2009) 168–182. [14] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, New York, New Rochelle, Melbourne, Sydney, 1986. [15] James M. Ortega, Numerical Analysis: A Second Course, Academic Press, New York, 1973. [16] G. Wang, Y. Wei, S. Qiao, Generalized Inverses: Theory and Computations, Science Press, 2004.