A High-Order Iterate Method for Computing 𝐴⁽²⁾_𝑇,𝑆

Abstract

We investigate a new higher order iterative method for computing the generalized inverse 𝐴⁽²⁾_𝑇,𝑆 for a given matrix 𝐴. We also discuss how the new method could be applied for finding approximate inverses of nonsingular square matrices. Analysis of convergence is included to show that the proposed scheme has at least fifteenth-order convergence. Some tests are also presented to show the superiority of the new method.

1. Introduction

The traditional generalized inverses, the Moore-Penrose inverse, the weighted Moore-Penrose inverse, the Drazin inverse, the group inverse, the Bott-Duffin inverse, and so forth are of special interest in matrix theory. They are extensively used in statistics, control theory, power systems, nonlinear equations, optimization and numerical analysis, and so on. Most of these generalized inverses are outer inverses of the form 𝐴⁽²⁾_𝑇,𝑆 having the prescribed range 𝑇 and null space 𝑆. For a given complex matrix 𝐴∈ℂ^𝑚×𝑛, the unique matrix 𝑋∈ℂ^𝑛×𝑚 such that 𝑋𝐴𝑋=𝑋, 𝑅(𝑋)=𝑇, and 𝑁(𝑋)=𝑆 is known as the outer inverse or 2-inverse denoted by 𝐴⁽²⁾_𝑇,𝑆 of 𝐴 with the prescribed range 𝑇 and null space 𝑆. This also plays an important role in singular differential and difference equations, the stable approximations of ill-posed problems, and in linear and nonlinear problems [1, 2]. As a result, it has been extensively studied by many researchers and many methods [3–16] are proposed in the literature for computing them. There are some techniques to tackle this problem, which are basically divided into two parts: the direct solvers such as Gaussian elimination with partial pivoting (GEPP). However, direct methods usually need much cost in both time and space in order to achieve desired results; sometimes they may not be able to work at all. In contrast, iterative method compensates for individual and accumulation of round-off errors as it is a process of successive refinement.

The well-known second order Schulz method in [17] firstly presents Newton’s method: 𝑋_𝑘+1=𝑋_𝑘2𝐼−𝐴𝑋_𝑘, (1) where 𝐼 is the identity matrix. In 2011, Li et al. in [18] gave the following locally cubically convergent scheme: 𝑋_𝑘+1=𝑋_𝑘3𝐼−𝐴𝑋_𝑘3𝐼−𝐴𝑋_𝑘 (2) and also proposed another iterative method: 𝑋_𝑘+1=𝑋_𝑘1𝐼+2𝐼−𝐴𝑋_𝑘𝐼+2𝐼−𝐴𝑋_𝑘. (3) Krishnamurthy and Sen in [19] provided the following fourth-order method: 𝑋_𝑘+1=𝑋_𝑘𝐼+𝑌_𝑘𝐼+𝑌_𝑘𝐼+𝑌_𝑘, (4) in which 𝑌_𝑘=𝐼−𝐴𝑋_𝑘. As another example from this primary source, the authors provided the following twelfth-order method: 𝑋_𝑘+1=𝑋_𝑘𝐼+𝑌_𝑘𝐼+𝑌_𝑘×𝐼+𝑌_𝑘𝐼+𝑌_𝑘×𝐼+𝑌_𝑘𝐼+𝑌_𝑘×𝐼+𝑌_𝑘𝐼+𝑌_𝑘. (5)

It is known that the Schulz iteration is numerically stable. Unfortunately, the method is too slow at the beginning of the process for finding 𝐴⁽²⁾_𝑇,𝑆. To remedy this drawback and lead to the better equilibrium between the high speed and the operational cost, we define a new higher order iterative method.

The aim of the present paper is to introduce and investigate new higher order iterative methods for computing the generalized inverse 𝐴⁽²⁾_𝑇,𝑆 for a given matrix 𝐴. We also discussed how the new method could be applied for finding approximate inverses of nonsingular square matrices. Analysis of convergence is included to show that the proposed scheme has at least fifteenth-order convergence. Some tests are also presented to show the superiority of the new method.

The paper is organized as follows. Section 1 is the introduction. In Section 2, some preliminaries involving notations are given. We analytically discuss the application of the new algorithm in the computation of 𝐴⁽²⁾_𝑇,𝑆 in Section 3. The proposed iterative method for computing matrix inversion for its convergence analysis is also given. In Section 4, we discuss the complexity of the iterative methods to theoretically find the most efficient method. Some numerical examples matrices are worked out in Section 5; the results obtained are compared with existing methods.

2. Preliminaries

In this section, we shall describe some concepts used in this paper. Let the symbols ℂ^𝑚×𝑛, 𝑇 and 𝑆 denote the set of 𝑚×𝑛 complex matrices, the subspace of 𝐶^𝑛, and the subspace of 𝐶^𝑚, respectively. In addition, the symbols 𝑅(𝑥), 𝑁(𝑥), 𝐴^∗, 𝜌(𝑋), and ‖𝑋‖ denote the range, the null space, the conjugate transpose, spectral radius, and Frobenius norm of 𝑋. It is well known that 𝐴 converges; that is, lim_𝑛→∞𝐴^𝑛=0 if and only if 𝜌(𝐴)<1.

The following lemmas are needed in what follows.

Lemma 1. If 𝑃_𝐿,𝑀 denotes the projector on a space 𝐿 along a space 𝑀, then (i)𝑃_𝐿,𝑀𝑄=𝑄 if and only if 𝑅(𝑄)⊆𝐿;(ii)𝑄𝑃_𝐿,𝑀=𝑄 if and only if 𝑁(𝑄)⊇𝑀.

Lemma 2. Let 𝐴∈ℂ^𝑚×𝑛; let 𝑇 and 𝑆 be the subspace of 𝐶^𝑛 and 𝐶^𝑚, respectively, with dim𝑇=dim𝑆^⟂=𝑡≤𝑟. Then, 𝐴 has a unique outer inverse 𝑋 denoting 𝐴⁽²⁾_𝑇,𝑆 such that 𝑅(𝑋)=𝑇 and 𝑁(𝑋)=𝑆 if and only if 𝐴𝑇⊕𝑆=𝐶^𝑚. (6)

It is well known that, for 𝐴∈ℂ^𝑚×𝑛, the Moore-Penrose inverse 𝐴^†, the Drazin inverse 𝐴^𝐷, the weighted Moore-Penrose inverse 𝐴^†_𝑀,𝑁, and the weighted Drazin inverse 𝐴^𝐷,𝑊 can be represented by(i)𝐴^†=𝐴⁽²⁾_{ℛ(𝐴∗),𝒩(𝐴∗)};(ii)𝐴^𝐷=𝐴⁽²⁾_{ℛ(𝐴𝑘),𝒩(𝐴𝑘)}, where 𝑘=Ind(𝐴);(iii)𝐴^†_𝑀,𝑁=𝐴⁽²⁾_{ℛ(𝑁−1𝐴∗𝑀),𝒩(𝑁−1𝐴∗𝑀)};(iv)𝐴^𝐷,𝑊=(𝑊𝐴𝑊)⁽²⁾_{ℛ((𝐴𝑊)𝑘𝐴),𝒩((𝐴𝑊)𝑘𝐴)}, where 𝑘=Ind(𝐴𝑊).

3. A New Method for Finding the Generalized Inverse 𝐴⁽²⁾_𝑇,𝑆

Inspired and motivated by the classical method (1) for finding the inverse of a matrix 𝐴, it should be pointed out that Newton’s method for the matrix equation [17] is as follows: 𝐹(𝑉)∶=𝑉⁻¹−𝐴=0. (7) To construct a new matrix iteration, we must find a nonlinear equation solver, which is more efficient than Newton’s method when it is applied to (7). Toward this goal, we apply the following rational iteration function: 𝑌_𝑘=𝑉_𝑘−𝐹^′𝑉_𝑘⁻¹𝐹𝑉_𝑘,𝑍_𝑘=𝑉_𝑘−2⁻¹𝐹𝑉_𝑘𝐹^′𝑉_𝑘⁻¹+𝐹^′𝑌_𝑘⁻¹,𝑊_𝑘=𝑍_𝑘−𝑍_𝑘−𝑌_𝑘⁻¹𝐹(𝑍_𝑘)−𝐹(𝑌_𝑘)⁻¹𝐹𝑍_𝑘,𝑉_𝑘+1=𝑊_𝑘−2⁻¹𝐹𝑊_𝑘𝐹^′𝑊_𝑘⁻¹×2+𝐹^′′𝑊_𝑘𝐹𝑊_𝑘𝐹^′𝑊_𝑘⁻² (8) on the matrix equation (7), to obtain the following efficient high-order method after simplifying: 𝑉_𝑘+11=−8𝑉_𝑘−11𝐼+25𝐴𝑉_𝑘−30𝐴𝑉_𝑘²+20𝐴𝑉_𝑘³−7𝐴𝑉_𝑘⁴+𝐴𝑉_𝑘⁵×12𝐼−66𝐴𝑉_𝑘+271𝐴𝑉_𝑘²−730𝐴𝑉_𝑘³+2060𝐴𝑉_𝑘⁴−1982𝐴𝑉_𝑘⁵+1260𝐴𝑉_𝑘⁶−1572𝐴𝑉_𝑘⁷+870𝐴𝑉_𝑘⁸−340𝐴𝑉_𝑘⁹+89𝐴𝑉_𝑘¹⁰−14𝐴𝑉_𝑘¹¹+𝐴𝑉_𝑘¹². (9) Now, by simplification, we suggest the following matrix iteration: 𝜓_𝑘=𝐴𝑉_𝑘,𝜒_𝑘=−11𝐼+𝜓_𝑘25𝐼+𝜓_𝑘−30𝐼+𝜓_𝑘20𝐼+𝜓_𝑘−7𝐼+𝜓_𝑘,𝜃_𝑘=𝜓_𝑘𝜒_𝑘,𝑉_𝑘+11=−8𝑉_𝑘𝜒_𝑘12𝐼+𝜃_𝑘6𝐼+𝜃_𝑘𝑛=0,1,2,…, (10) where 𝐼 is the identity matrix.

Note that it is also known that for 𝑌∈ℂ_𝑟^𝑛×𝑚 with 𝑅(𝑌)=𝑇 and 𝑁(𝑌)=𝑆, and the initial approximation 𝑉₀=𝛼𝑌, with 0<𝛼<2/(𝜌(𝐴𝑌)), it holds that 𝐼−𝐴𝑉₀<1. (11) We are about to use this fact in the following theorem so as to find the theoretical order of the proposed method (10) for finding the generalized inverse 𝐴⁽²⁾_𝑇,𝑆.

Theorem 3. Let 𝐴∈ℂ_𝑟^𝑚×𝑛, let 𝑇 be a subspace of 𝐶^𝑛 of dimension 𝑠≤𝑟, and let 𝑆 be a subspace of 𝐶^𝑚 of dimension 𝑚−𝑠. Suppose 𝑌∈ℂ_𝑟^𝑛×𝑚 with 𝑅(𝑌)=𝑇 and 𝑁(𝑌)=𝑆 and the initial approximation 𝑉₀=𝛼𝑌; the sequence {𝑉_𝑘}^𝑘=∞_𝑘=0 generalized by (9) satisfies the following error estimate when finding 𝐴⁽²⁾_𝑇,𝑆: 𝐴⁽²⁾_𝑇,𝑆−𝑉_𝑘𝐴⁽²⁾_𝑇,𝑆≤𝐼−𝐴𝑉₀^15𝑘,𝑘=0,1,2,…. (12)

Proof. We use notations that 𝐸₀=𝐼−𝐴𝑉₀ and subsequently 𝐸_𝑘=𝐼−𝐴𝑉_𝑘. Then, 𝐸_𝑘+1=𝐼−𝐴𝑉_𝑘+11=𝐼−8𝐴𝑉_𝑘−11𝐼+25𝐴𝑉_𝑘−30𝐴𝑉_𝑘²+20𝐴𝑉_𝑘³−7𝐴𝑉_𝑘⁴+𝐴𝑉_𝑘⁵×12𝐼−66𝐴𝑉_𝑘+271𝐴𝑉_𝑘²−730𝐴𝑉_𝑘³+1405𝐴𝑉_𝑘⁴−1982𝐴𝑉_𝑘⁵+2060𝐴𝑉_𝑘⁶−1472𝐴𝑉_𝑘⁷+870𝐴𝑉_𝑘⁸−340𝐴𝑉_𝑘⁹+89𝐴𝑉_𝑘¹⁰−14𝐴𝑉_𝑘¹¹+𝐴𝑉_𝑘¹²=182𝐼−𝐴𝑉_𝑘³𝐼−𝐴𝑉_𝑘¹⁵=18𝐼+𝐸_𝑘³𝐸_𝑘¹⁵=18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸. (13) By taking an arbitrary matrix norm on both sides of (13), we attain 𝐸_𝑘+1≤18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸. (14) In addition, since ‖𝐸₀‖≤1, by relation (14), we obtain that 𝐸₁≤18𝐸₀¹⁵+3𝐸₀¹⁶+3𝐸₀¹⁷+𝐸₀¹⁸≤𝐸₀¹⁵. (15) Now if we consider ‖𝐸_𝑘‖<1, therefore 𝐸_𝑘+1≤18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸≤𝐸_𝑘¹⁵. (16) Using mathematical induction, we obtain 𝐸_𝑘+1≤𝐸₀^15𝑘,𝑘≥0. (17) Using this fact 𝑉₀=𝛼𝑌 with 𝑅(𝑌)=𝑇 and 𝑁(𝑌)=𝑆 in conjunction with (9), we conclude that 𝑅𝑉_𝑘⊆𝑇,𝑁𝑉_𝑘⊇𝑆,𝑘≥0. (18) It is not difficult to verify that 𝐴⁽²⁾_𝑇,𝑆𝐴𝑉_𝑘=𝑉_𝑘,𝑉_𝑘𝐴𝐴⁽²⁾_𝑇,𝑆=𝑉_𝑘,𝑘≥0. (19) Therefore, the error matrix 𝐹_𝑘=𝐴⁽²⁾_𝑇,𝑆−𝑉_𝑘 satisfies 𝐹_𝑘=𝐴⁽²⁾_𝑇,𝑆−𝑉_𝑘=𝐴⁽²⁾_𝑇,𝑆−𝐴⁽²⁾_𝑇,𝑆𝐴𝑉_𝑘=𝐴⁽²⁾_𝑇,𝑆𝐼−𝐴𝑉_𝑘=𝐴⁽²⁾_𝑇,𝑆𝐸_𝑘. (20) From the last identity and (17) we have 𝐹_𝑘≤𝐴⁽²⁾_𝑇,𝑆𝐸_𝑘≤𝐴⁽²⁾_𝑇,𝑆𝐸₀^15𝑘, (21) which is a confirmation of (12).

Theorem 4. Let 𝐴∈ℂ_𝑟^𝑚×𝑛, let 𝑇 be a subspace of 𝐶^𝑛 of dimension 𝑠≤𝑟, and let 𝑠 be a subspace of 𝐶^𝑚 of dimension 𝑚−𝑠. Suppose 𝑌∈ℂ_𝑟^𝑛×𝑚 with 𝑅(𝑌)=𝑇 and 𝑁(𝑌)=𝑆. If the initial approximation 𝑉₀=𝛼𝑌 is in accordance with (11). The sequence {𝑉_𝑘}^𝑘=∞_𝑘=0 generalized by (9) converges with fifteenth order to 𝐴⁽²⁾_𝑇,𝑆.

Proof. From (13), we know 𝐸_𝑘+1=18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸. (22) Set 𝐹_𝑘=𝑉_𝑘−𝐴⁽²⁾_𝑇,𝑆. We have 𝐴𝐹_𝑘+1=𝐴𝑉_𝑘+1−𝐴𝐴⁽²⁾_𝑇,𝑆=𝐴𝑉_𝑘+1−𝐼+𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆=−𝐸_𝑘+1+𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆1=−8𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸+𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆. (23) On the other hand, from conditions (18) and (19), we can obtain the equation 𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆^𝑡=𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆,𝑡=1,2,3,…,𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆𝐴𝐹_𝑘=𝐼−𝐴𝐴⁽²⁾_𝑇,𝑆𝐴𝑉_𝑘−𝐴𝐴⁽²⁾_𝑇,𝑆=𝐴𝑉_𝑘−𝐴𝐴⁽²⁾_𝑇,𝑆−𝐴𝐴⁽²⁾_𝑇,𝑆𝐴𝑉_𝑘+𝐴𝐴⁽²⁾_𝑇,𝑆𝐴𝐴⁽²⁾_𝑇,𝑆=𝐴𝑉_𝑘−𝐴𝐴⁽²⁾_𝑇,𝑆−𝐴𝑉_𝑘+𝐴𝐴⁽²⁾_𝑇,𝑆=0. (24)(25) According to the relations (22)–(25), we have 𝐴𝐹_𝑘+1=18𝐹_𝑘¹⁵+3𝐹_𝑘¹⁶+3𝐹_𝑘¹⁷+𝐹_𝑘¹⁸. (26) Therefore 𝐴𝐹_𝑘+1≤18𝐴𝐹_𝑘¹⁵+3𝐴𝐹_𝑘¹⁶+3𝐴𝐹_𝑘¹⁷+𝐴𝐹_𝑘¹⁸≤𝐴𝐹_𝑘¹⁵. (27) It would be now easy to find the error inequality of the high-order iterative as follows: 𝐹_𝑘+1=𝑉_𝑘+1−𝐴⁽²⁾_𝑇,𝑆=𝐴⁽²⁾_𝑇,𝑆𝐴𝑉_𝑘+1−𝐴⁽²⁾_𝑇,𝑆𝐴𝐴⁽²⁾_𝑇,𝑆≤𝐴⁽²⁾_𝑇,𝑆𝐴𝑉_𝑘+1−𝐴𝐴⁽²⁾_𝑇,𝑆=𝐴⁽²⁾_𝑇,𝑆𝐴𝐹_𝑘+1≤𝐴⁽²⁾_𝑇,𝑆‖𝐴‖¹⁵𝐹_𝑘¹⁵. (28)
Thus, ‖𝑉_𝑘+1−𝐴⁽²⁾_𝑇,𝑆‖→0; that is, the sequence {𝑉_𝑘}^𝑘=∞_𝑘=0 converges to the generalized inverse 𝐴⁽²⁾_𝑇,𝑆 in at least fifteenth-order as 𝑛→0. This ends the proof.

Since the Moore-Penrose inverse, the weighted Moore-Penrose inverse, the Drazin inverse, and the weighted Drazin inverse of a matrix over ℂ always exist, the below corollary follows immediately by Theorems 4.

Corollary 5. Suppose the condition of Theorem 3 holds. Let 𝑉_𝑘 be defined by Theorem 3. Denote 𝑉_∞=lim_𝑘→∞𝑉_𝑘, the sequence {𝑉_𝑘}^𝑘=∞_𝑘=0 generalized by (9) and using the initial approximation 𝑉₀=𝛼𝑌 and such that 0<𝛼<2/(𝜌(𝐴𝑌)), if and only if one of the following holds:(i)𝑉_∞=𝐴^† with 𝑌=𝐴^∗;(ii)𝑉_∞=𝐴^𝐷 with 𝑌=𝐴^𝑙, 𝑙≥Ind(𝐴);(iii)𝑉_∞=𝐴^†_𝑀,𝑁 with 𝑌=𝑁⁻¹𝐴^∗𝑀, where 𝑀 and 𝑁 are Hermite positive definite matrices of the orders 𝑚 and 𝑛, respectively;(iv)𝑉_∞=𝐴^𝐷,𝑊 with 𝑌=(𝐴𝑊)^𝑙𝐴, 𝑙≥Ind(𝐴𝑊), where 𝑊∈ℂ^𝑛×𝑚.

By taking special matrix, we can get some desired results. We discuss the applications of (9) in computing matrix inversion 𝐴⁻¹ as follows.

Corollary 6. Let 𝐴=[𝑎_𝑖𝑗]∈ℂ^𝑛×𝑛 be a nonsingular real or complex matrix. If the initial approximation 𝑉₀ satisfies 𝐼−𝐴𝑉₀<1, (29) then the iterative method (9) converges with at least fifteenth-order to 𝐴⁻¹.

Proof. From (13), we know 𝐸_𝑘+1=𝐼−𝐴𝑉_𝑘+1=18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸. (30) Hence, by taking arbitrary matrix norm on both sides of (30), we attain 𝐸_𝑘+1≤18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸. (31) In addition, since ‖𝐸₀‖<1 by relation (13), we obtain that 𝐸₁≤18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸≤𝐸_𝑜¹⁵<1. (32) Now if we consider ‖𝐸_𝑘‖<1, therefore 𝐸_𝑘+1≤18𝐸_𝑘¹⁵+3𝐸_𝑘¹⁶+3𝐸_𝑘¹⁷+𝐸_𝑘¹⁸≤𝐸_𝑘¹⁵<1. (33) Using mathematical induction, we obtain 𝐸_𝑘+1≤𝐸_𝑘¹⁵,𝑘≥0. (34) Furthermore, we get that 𝐸_𝑘+1≤𝐸_𝑘¹⁵≤⋯≤𝐸₀^15𝑘+1. (35) That is, 𝐼−𝐴𝑉_𝑘→0, when 𝑘→∞ and 𝑉_𝑘⟶𝐴⁻¹,as𝑘⟶∞. (36)
Thus, the new method (9) converges to the inverse of the matrix 𝐴 in the case 𝜌(𝐴𝑉₀)<1, where 𝜌 is the spectral radius. Now, we prove that the order of convergence for the sequence {𝑉_𝑘}^𝑘=∞_𝑘=0 is at least fifteen. Let 𝜎_𝑘 denote the error matrix 𝜎_𝑛=𝐴⁻¹−𝑉_𝑘; afterwards 𝐸_𝑘=𝐴𝜎_𝑘=𝐼−𝐴𝑉_𝑘. (37) The identity (13) in conjunction with (37) implies that 𝐴𝜎_𝑘+1=18𝐴𝜎_𝑘¹⁵+3𝐴𝜎_𝑘¹⁶+3𝐴𝜎_𝑘¹⁷+𝐴𝜎_𝑘¹⁸. (38) Therefore, using invertibility of 𝐴, it follows immediately that 𝜎_𝑘+1=18𝐴𝜎_𝑘¹⁴+3𝐴𝜎_𝑘¹⁵+3𝐴𝜎_𝑘¹⁶+𝐴𝜎_𝑘¹⁷𝜎_𝑘. (39) By taking any subordinate norm of (39), we obtain 𝜎_𝑘+1=18‖𝐴‖¹⁴+3‖𝐴‖¹⁵𝜎_𝑘+3‖𝐴‖¹⁶𝜎_𝑘²+‖𝐴‖¹⁷𝜎_𝑘³𝜎_𝑘¹⁵. (40) Consequently, it is proved that the iterative formula (9) converges to 𝐴⁻¹, and the order of this method is at least fifteen.

4. Complexity of the Methods

Let us consider two parameters 𝜌 and 𝜅 which stand for the rate of convergence and the number of matrix-by-matrix multiplications in floating point arithmetics, respectively. Then the comparative index could be expressed by 𝜌𝐸=𝜅. (41) According to Table 1, we can see that the iterative process (9) reduces the computational complexity by using less number of basic operations and leads to the better equilibrium between the high speed and the operational cost.

Table 1 

Comparison of the computational complexity for different methods.

5. Numerical Examples

In this section, we will illustrate a numerical experiment on our method and give some theoretical and numerical comparison with the existing solution methods. All these numerical experiments are executed in the Matlab 7.13.0.564 (R2011b) (Build 7601: Service Pack 1) of the software used on an intel(R) core (TM)2 Duo CPU T6500 2.10 GHz with 2 GB RAM memory running on the windows XP Professional Version 2002 Service Pack 3.

Example 1. Consider the singular M-matrix 𝐴 of order 6×6 in [20] 𝐴=.1−10000−110000−1−11−100−1−1−1100−1−1−102−1−1−10−1−12 (42) Clearly, 𝑘=ind(𝐴)=2. Now, taking 𝑌=𝐴^𝑘 given by 𝑌=,2−20000−220000002−20000−220000−325−4002−3−45 (43) such that 𝑅(𝑌)=𝑇, 𝑁(𝑌)=𝑆. We must take 𝛼 which satisfies 0<𝛼<0.074074 to satisfy the convergence condition 0<𝛼<2/(𝜌(𝐴𝑌)). The Drazin inverse 𝐴^𝑑 is obtained as 𝐴^𝐷=14−14−1000041410000004−1410000−41450000−−712212313700−−512112323. (44)

Table 2 compares the number of iterations (𝑁) and error bounds of the existing methods [17, 18] and the fastest method of [20] with the stopping criteria ‖𝑋_𝑘+1−𝑋_𝑘‖≤10⁻¹⁰. It is found that our method takes fewer numbers of iterations and gives accurate estimations of 𝐴^𝑑 in comparison to the method compared.

Table 2 

Comparison of 𝑁 and error bounds for different values of 𝛼.

Example 2. Consider the following 12×12 matrix investigated in ([21] Example 3.2): 𝐴=.20.40000000000−20.40000000000−1−11−10000−1000−1−1−1100000000000011−1−100−10000011−1−10000000−1−20.4000000000020.40000000−10000001−1−1−100000000−11−1−100000000000.4−200000000000.42 (45) Then 𝑘=ind(𝐴)=3. Now, 𝑌=𝐴^𝑘 with 𝛼=2/tr(𝐴^𝑘+1).

According to Table 3, our method algorithm only needs 5 iterations to converge with ‖𝑋_𝑘−𝑋_𝑘−1‖≤10⁻⁸. We obtain 𝐴^𝐷=[𝐴₁₁𝐴₁₂]. Consider 𝐴₁₁=,𝐴0.25−0.2500001.251.250000−1.6641−0.99220.25−0.2500−1.1953−0.6797−0.250.2500−2.7637−1.0449−1.875−1.25−1.251.25−2.7637−1.0449−1.875−1.25−1.251.2514.10946.30086.6253.3755−3−19.3242−8.5078−9.75−5.25−7.54.5−0.625−0.31250000−1.25−0.93750000000000000000₁₂=.00000000000000−0.0625−0.062500.156300−0.06250.18750.68751.34381.251.251.4843842.57813.32036.64061.251.251.48442.57814.57038.5156−5−5−4.1875−8.5−10.5078−22.46097.57.56.37512.562515.976633.7891000.25−0.25−0.875−1.62500−0.250.25−0.875−1.62500001.251.250000−0.250.25 (46)

Table 3 

Comparisons.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Xiaoji Liu and Zemeng Zuo were supported by the Guangxi Natural Science Foundation (2013GXNSFAA019008), the Key Project of Education Department of Guangxi (201202ZD031), Project supported by the National Science Foundation of China (11361009), and Science Research Project 2013 of the China-ASEAN Study Center (Guangxi Science Experiment Center) of Guangxi University for Nationalities.