DOI QR코드

DOI QR Code

HIGHER ORDER ITERATIONS FOR MOORE-PENROSE INVERSES

  • Srivastava, Shwetabh (Department of Mathematics, IIT Kharagpur) ;
  • Gupta, D.K. (Department of Mathematics, IIT Kharagpur)
  • Received : 2013.03.10
  • Accepted : 2013.04.26
  • Published : 2014.01.30

Abstract

A higher order iterative method to compute the Moore-Penrose inverses of arbitrary matrices using only the Penrose equation (ii) is developed by extending the iterative method described in [1]. Convergence properties as well as the error estimates of the method are studied. The efficacy of the method is demonstrated by working out four numerical examples, two involving a full rank matrix and an ill-conditioned Hilbert matrix, whereas, the other two involving randomly generated full rank and rank deficient matrices. The performance measures are the number of iterations and CPU time in seconds used by the method. It is observed that the number of iterations always decreases as expected and the CPU time first decreases gradually and then increases with the increase of the order of the method for all examples considered.

Keywords

1. Introduction

The theory of generalized inverses has seen a substantial growth over the past few decades. Many applications of statistics, prediction theory, control analysis and numerical analysis often require computation of generalized inverses. The Moore-Penrose inverse is one of the most important generalized inverses of arbitrary singular square or rectangular (real or complex) matrix. It has been extensively studied by many researchers [4,2,6,9,5,3] and many methods are proposed in the literature. Accordingly, it is important both practically and theoretically to find good higher order algorithms for computing a Moore-Penrose inverse of a given arbitrary matrix. Let ℂm×n and denote the set of all complex (m × n) matrices and all complex (m × n) matrices with rank r, respectively. For A ∈ ℂm×n, it is denoted by A† and defined as

where, R(A) and PR(A)denote the range space of A and orthogonal projection on R(A) respectively. The unique matrix A† satisfies the following four equations

Both direct and iterative methods (cf. [11,10,2,3,9]) can be used to compute A†. One of the most commonly used direct methods is the Singular Value Decomposition (SVD) method. For A ∈ ℂm×n, the SVD method is a factorization of the form

where, U is an m×m complex unitary matrix, Σ is an m×n rectangular diagonal matrix with nonnegative real numbers on the diagonal and V is an n×n complex unitary matrix. The A† can then be written as

where, Σ† is the generalized inverse of Σ obtained by replacing every nonzero diagonal entry in Σ by its reciprocal and then transposing the resulting matrix. This method is very accurate but time-intensive since it requires a large amount of computational resources, especially in the case of large matrices. The most frequently used iterative method for approximating A−1 is the method studied by Householder [14] by analyzing successive improvements of a matrix X to solve AX = M, A nonsingular, using the equation

One version of the above has the matrices Ck project on the ith row, where i cycles through all the rows of A. By taking M = I and Ck = Vk, the famous Newton’s method originated in [21] is

As usual, I denotes the identity matrix of an appropriate order. If L is the desired limit matrix and Vk is the k-th estimate of L, then the convergence properties of examined iterative method can be studied with the aid of error matrix Ek = Vk − L. If an iterative method is expressible as a simple matrix formula, Ek+1 is sum of several terms

All suitable algorithms have a zero-order term equal to 0. Hence, the firstorder term determine the terminal convergence properties. The eigenvalues of are given by

It is further established that the eigenvalues of I − AV0 must have magnitude less than 1 to ensure the convergence of (3). Since the residuals Rk = I −AVk in (3) satisfies ∥Rk+1∥ ≤ ∥A∥∥Rk∥2, Newton’s method is a second order iterative method [7]. Ben-Israel [11] used (3) and the starting value X0 = αA∗, where α satisfies

The iterative scheme (3) is generalized by the iteration Uk+1 = Uk(2PR(A) − AUk), which converges to A†. Ben-Israel and Charnes [13] proved that the sequence

converges to A† under the assumption (5). Li et al.[6] established a family of iterative methods to compute the approximate inverse of a square matrix and inner inverse of a non-square matrix given by

Chen and Wang [4] showed that this can be extended to compute A† with higher orders. Katskis et.al.[8] developed a much faster computational method to calculate the A† of singular square matrices and of rectangular matrices A. This method has significantly better accuracy than the already proposed methods and works for full and sparse matrices. Weiguo et al.[5] improves on the generalized properties of a family of iterative methods to compute the approximate inverses of square matrices and fixed inner inverses of rectangular matrices proposed in [6]. They have also established that this fixed inverse is the Moore-Penrose inverse of the considered matrix. Vasilios et.al.[3] also presented a very fast and reliable method to compute Moore-Penrose inverse. By using a general framework where analytic functions of scalers are first developed and then matrices substituted for the scalers, Katsaggelos and Efstratiadis [12] produced a convergence faster than quadratic, for restricted initial estimates.

In this paper, a higher order iterative method to compute the Moore-Penrose inverses of arbitrary matrices using only the Penrose equation (ii) is developed by extending the iterative method described in [1]. Convergence properties as well as the error estimates of the method are studied. The efficacy of the method is demonstrated by working out four numerical examples, two involving a full rank simple and an ill-conditioned Hilbert matrix, whereas, the other two involving full rank and rank deficient randomly generated matrices. The performance measures are the number of iterations and CPU time in seconds used by the method. It is observed that the number of iterations always decreases as expected and the CPU time first decreases gradually and then increases with the increase of the order of the method for all examples considered.

The paper is organized as follows. Section 1 is the introduction. In Section 2, the higher order iterative method for computing the Moore-Penrose generalized inverse of an arbitrary rectangular matrix is described. Some Lemmas and the convergence analysis of the method are also established. In Section 3, the efficacy of the method is demonstrated by working out four numerical examples, two involving a full rank simple and an ill-conditioned Hilbert matrix, whereas, the other two involving full rank and rank deficient randomly generated matrices. Finally, conclusions are included in Section 4.

 

2. Higher order iteration with convergence Analysis

In this section, higher order iterative methods and their convergence analysis to compute the Moore-Penrose inverse A† of an arbitrary complex matrix A ∈ ℂm×n is described. Let At, A∗, μ(A), ν(A) and rank(A), represent the transpose, the conjugate transpose, the range space, the null space and the rank of the matrix A ∈ ℂm×n, respectively.

Definition 2.1 ([19]). Let A ∈ ℂm×n, with column vectors v1, v2, . . . , vn, the set

are called the Range (the Image, or the Column Space) of A. The set

is called the Null Space (or the Kernel) of A and for any A ∈ ℂm×n,B ∈ ℂs×t, we call

the Range of (A,B) and the Null Space of (A,B) respectively.

We shall now describe the higher order iterative method extending the method described in [1] for computing the Moore-Penrose generalized inverses. Petković et al. [1] proposed a quadratically convergent iterative method for computing A† based on Penrose equations (ii) and (iv) as follows. Let

Hence, for arbitrary β ∈ ℝ, we get

or equivalently,

This leads to the following iterative method

where, X0 = βA∗ for an appropriate real number β. For β < 1, the method (8) has a linear convergence while for β = 1 its convergence is quadratic. The first-order and the second-order error terms of (8) are

where, Ek = Xk−A† is the error matrix. Now, using only the Penrose equations (ii) given by X = XAX and for arbitrary β ∈ ℝ, we get

or, equivalently

This leads to the following third order extension of method (8)

for X0 = βA∗. Continuing in a similar manner, this can further be extended to the pth order for p ≥ 2, given by

for X0 = βA∗, where β is an appropriate real number. The following Lemmas will be used for establishing the convergence of these iterative methods.

Lemma 2.2. For all k ≥ 0, the sequence {Xk} generated by (10) and (11) satisfies

Proof. This Lemma can be proved by induction. Clearly, (1) holds for k = 0 as X0A = βA∗A = (X0A)∗. Assume that it holds for some k. To show that it also holds for k + 1, we consider

Hence it holds for all k ≥ 0. Likewise for pth order method (11), we get

Clearly, (2) trivially holds for k = 0. Let it also holds for some k. To show that it also holds for k + 1, we get

Likewise for pth order method (11), we get

Proceeding in a similar manner, (3) can easily be proved for (10) and (11). Hence, the Lemma 2.1 is proved for all k ≥ 0.

Theorem 2.3. Let 0 ≠ A ∈ ℂm×n, X = A†, the initial approximation X0 = βA∗, β ∈ (0, 1] and its residual R0 = (X0−X)A satisfy ∥R0∥ < 1. The sequence {Xk} generated by (10) starting with X0 = βA∗ converges to the Moore-Penrose inverse A†. It has linear convergence for β ≠ 1 and third order convergence for β = 1. Its first, second and the third order error terms are given by

where, Ek = Xk − A† denotes the error matrix.

Proof. Using Lemma 2.2 and substituting for Xk+1, we get

Using Lemma 2.2 and (10), we get

Thus, the sequence of residual matrices defined by Rk = XkA−XA satisfies the following recurrence relation

Let sk = ∥Rk∥. Now, for the convergence of the sequence {Xk}, we require that sk → 0 as k → ∞. This can be shown by mathematical induction. It trivially holds for k = 0, since, s0 = ∥R0∥ = ∥X0A − XA∥ < 1. Let it holds for some k, i.e., sk < 1. To show that it also holds for k + 1, we take norm on (13) and get

Thus, sk is a monotonically decreasing bounded sequence converging to s as k → ∞ and 0 ≤ s < 1. From (14), we get

This gives either s = 0 or s ≥ 1. Thus, s = 0. This completes the proof that sk → 0 when k → ∞. Thus, Xk → X as k → ∞. Substituting Xk = X + Ek in (10) and rearranging the terms, we get Ek+1 given by

This leads to

and

One can easily see that the method is linear convergent for β ≠ 1 and of third order convergent for β = 1 as error1 and error2 are equal to 0.

Theorem 2.4. Let 0 ≠ A ∈ ℂm×n, X = A†, the initial approximation X0 = βA∗, β ∈ (0, 1] and its residual R0 = (X0−X)A satisfy ∥R0∥ < 1. The sequence {Xk} generated by (11) starting with X0 = βA∗ converges to the Moore-Penrose inverse A†. It has linear convergence for β ≠ 1 and pth order convergence for β = 1, where, p ≥ 2 is a positive integer. All error terms of the method vanishes except the first and the pth order error terms given by

Proof. Using Lemma 2.2 and substituting for Xk+1, we get

Again using Lemma 2.2 and (11), this gives

Thus, the sequence of residual matrices defined by Rk = XkA − XA satisfies the following recurrence relation

Proceeding similar to theorem 2.3, it can be proved that Xk → X as k → ∞. Substituting Xk = X+Ek in (11) and rearranging the terms, we get Ek+1 given by

On simplification, we get all error terms equal to zero except error1 and errorp given by

One can easily see that the method is linearly convergent for β ≠ 1 and of pth order convergent for β = 1 as all other error terms vanishes. Hence, the theorem is proved.

We must note here that the convergence of the above methods (10) and (11) require the condition ∥(βA∗ − X)A∥ < 1. To verify this, we need an equivalent condition which does not contain the Moore-Penrose inverse X or A†. We follow exactly the same way as done in [1] by using the following Lemma.

Lemma 2.5 ([1]). Let the eigenvalues of matrix A∗A satisfies (4). Condition ρ((βA∗ − X)A) < 1 is satisfied under the assumption

To establish this Lemma, we need the following Lemmas.

Lemma 2.6 ([15]). Let M ∈ ℂn×n and ϵ > 0 be given. There is at least one matrix norm ∥.∥ such that

where ρ(M) = max{|λ1(M)|, . . . , |λn(M)|} denotes the spectral radius of M.

Lemma 2.7 ([18]). If P ∈ ℂn×n and S ∈ ℂn×n are such that P = P2 and PS = SP then

Remark 2.1. For any A ∈ ℂm×n, the sequence generated by our higher order iterative methods starting with X0 = βA∗ are convergent for any β satisfying 0 < β < 2/σ2(A), where σ(A) = ∥A∥2.

 

3. Numerical examples

In this section, four numerical examples are worked out to demonstrate the efficacy of our pth order method for various values of p. The mean CPU time in second and the number of iterations are measured as performance parameters. For full rank matrices the number of iterations used are compared by plotting figures with x-axis representing the values of p and y-axis representing the number of iterations. For randomly generated matrices, we have tested 50 times matrices and the mean CPU time in second are measured and tabulated for values of 2 ≤ p ≤ 7 and the number of iterations are compared by plotting figures with x-axis representing the times of choice of random matrices and y-axis representing the number of iterations only for p = 2 and p = 3. The value of ϵ = 10−7 and the maximum number of iterations equal to 100 are taken as the termination criteria. All the examples are run on an Intel core 2 Duo processor running at 2.80GHz and using MATLAB 7.4 (R2009b).

Example 3.1. Consider the matrix A of order (5 × 4) given by

The choice β = 0.60 satisfies the convergence criteria given by

since the eigenvalues of the matrix A∗A are

The iterative method (11) generates a sequence of iterates {Xk} converging to the Moore-Penrose generalized inverse A† given by

The comparison of number of iterations are plotted in Figure 1. It can be observed from Figure 1 that the iterative method (11) converges to the Moore-Penrose generalized inverse A† in 36 iterations for p = 2 and as expected as the order p increases, it reduces to 25 for p = 8.

FIGURE 1.Number of iterations versus the value of p for example 3.1

Example 3.2. Consider a (30 × 30) matrix A whose elements are generated randomly from [−0.2, 0.2]. Taking the termination criteria as

where, ∥.∥F stands for the Frobenius-norm of a matrix, Table 1 shows the comparison of mean CPU time in seconds for for 2 ≤ p ≤ 7. One can easily see that with the increase of p, the CPU time decreases initially and then starts increases after p > 4. The comparison of number of iterations are plotted in Figure 2 for p = 2 and p = 3. Similar results are also observed for higher order random matrices, for example, (40×40) matrix A whose elements are also randomly generated from [−0.2, 0.2].

TABLE 1.Comparison of CPU time

FIGURE 2.Comparison of number of iterations for example 3.2

Example 3.3. Consider a (100 × 50) rank deficient matrix A whose elements are generated randomly from [−0.2, 0.2]. Taking the termination criteria as

where, ∥.∥F stands for the Frobenius-norm of a matrix, Table 2 shows the comparison of mean CPU time in seconds for for 2 ≤ p ≤ 7. One can easily see that with the increase of p, the CPU time decreases initially and then starts increases after p ≥ 4. The comparison of number of iterations are plotted in Figure 3 for p = 2 and p = 3.

FIGURE 3.Comparison of number of iterations for example 3.3

TABLE 2.Comparison of CPU time

Example 3.4. Consider the ill-conditioned Hilbert matrix A of order (5 × 5) given by

The choice β = 0.80 satisfies the convergence criteria given by

since the eigenvalues of the matrix A∗A are

The iterative method (11) generates a sequence of iterates {Xk} converging to the Moore-Penrose generalized inverse A† given by

The comparison of number of iterations are plotted in Figure 4. It can be observed from Figure 4 that the iterative method (11) converges to the Moore-Penrose generalized inverse A† in 51 iterations for p = 2 and as expected as the order p increases, it reduces to 19 for p = 10.

FIGURE 4.Number of iterations versus the value of p for example 3.4

 

4. Conclusions

A quadratically convergent iterative method proposed by [1] is extended to a family of higher order iterative methods to compute the Moore-Penrose generalized inverses of arbitrary singular or rectangular (real or complex) matrices. This extension is carried out by using only the Penrose equation (ii) in place of (ii) and (iv) as used in [1]. Convergence properties as well as the error estimations are studied. The efficacy of the method is demonstrated by working out four numerical examples involving full rank and rank deficient randomly generated matrices and ill-conditioned Hilbert matrix. The performance in terms of computational time and the number of iterations are evaluated with respect to the order of the iterative methods. It is observed that the number of iterations always decreases as expected and the CPU time first decreases gradually and then increases with respect to the order of the iterative methods for all examples.

References

  1. Marko D. Petkovic and Predrag S.Stanimirovic, Iterative method for computing Moore-Penrose inverse based on Penrose equations, Journal of Computational and Applied Mathematics 235 (2011), 1604-1613. https://doi.org/10.1016/j.cam.2010.08.042
  2. W. Guo and T. Huang, Method of elementary transformation to compute Moore-Penrose inverse, Appl. Math. Comput. 216 (2010), 1614-1617. https://doi.org/10.1016/j.amc.2010.03.016
  3. Vasilios N. Katsikis, Dimitrios Pappas and Athanassios Petralias, An improved method for the computation of the Moore-Penrose inverse matrix, Appl. Math. Comput. 217 (2011), 9828-9834. https://doi.org/10.1016/j.amc.2011.04.080
  4. Haibin Chen and Yiju Wang, A Family of higher-order convergent iterative methods for computing the Moore-Penrose inverse, Appl. Math. Comput. 218 (2011), 4012-4016. https://doi.org/10.1016/j.amc.2011.05.066
  5. Li Weiguo, Li Juan and Qiao Tiantian, A family of iterative methods for computing Moore-Penrose inverse of a matrix, Linear Algebra and its Applications 438 (2013), 47-56. https://doi.org/10.1016/j.laa.2012.08.004
  6. Weiguo Li and Zhi Li, A family of iterative methods for computing the approximate inverse of a square matrix and inner inverse of a non-square matrix, Appl. Math. Comput. 215 (2010), 3433-3442 https://doi.org/10.1016/j.amc.2009.10.038
  7. V. Pan and R. Schreiber, An improved Newton iteration for the generalized inverse of a matrix with applications, SIAM J. Sci. Statist. Comput. 12 (1991), 1109-1130. https://doi.org/10.1137/0912058
  8. V.N. Katsikis, D. Pappas and T. Huang, Fast computing of the Moore-Penrose inverse matrix, Electronic Journal of Linear Algebra 17 (2008), 637-650.
  9. F. Toutounian and A. Ataei, A new method for computing Moore-Penrose inverse matrices, Journal of Computational and applied Mathematics 228 (2009), 412-417. https://doi.org/10.1016/j.cam.2008.10.008
  10. P. Courrieu, Fast Computation of Moore-Penrose Inverse matrices, Neural Information Processing-Letters and Reviews 8 (2005), 25-29.
  11. A. Ben-Israel and T. N. E. Greville, Generalized Inverses: Theory and Applications, Springer-Verlag, NewYork, 2003.
  12. A.K.Katsaggelos and S.N. Efstratiadis, A class of iterative signal restoration algorithms, IEE Trans. Acoust. Speech Signal Process, 38 (1990), 778-793. https://doi.org/10.1109/29.56022
  13. A. Ben-Israel and A.Charnes, Contribution to the theory of generalized inverses, J.Soc.Indust.Appl.Math. 11 (1963), 667-669. https://doi.org/10.1137/0111051
  14. A.S.Householder, The Theory of Matrices in Numerical Analysis, Dover Publications, Newyork, 1964.
  15. R.A. Horn and C.R. Johnson, Matrix Analysis Cambridge University Press, Cambridge, UK, 1986.
  16. F.J. Hall and C.D. Meryer Jr., Generalized inverses of the fundamental bordered matrix used in linear estimation, Sankhya Ser. A 37 (1975), 428-438.
  17. B. Nacevska, Iterative methods for computing generalized inverses and splittings of operators, Appl. Math. Comput., 208 (2009), 186-188. https://doi.org/10.1016/j.amc.2008.11.030
  18. P.S.Stanimirovic and D.S. Cvetkovicillic, Succesive matrix squaring algorithm for com-puting outer inverses, Appl. Math. Comput., 203 (2008), 19-29. https://doi.org/10.1016/j.amc.2008.04.037
  19. W. Sun and Y. Yuan, Optimization Theory and Methods- Nonlinear Programming, Springer, New York, 2006.
  20. Stephen L. Campbell, and Carl D. Meyer, Generalized Inverses of Linear Transformations, SIAM, 2009.
  21. G.Schultz, iterative Berechmmg der reziproken matrix, Angew. Math. Mech., 13 (1933), 57-59. https://doi.org/10.1002/zamm.19330130111