A Higher Order Iterative Method for Computing the Drazin Inverse

Abstract

A method with high convergence rate for finding approximate inverses of nonsingular matrices is suggested and established analytically. An extension of the introduced computational scheme to general square matrices is defined. The extended method could be used for finding the Drazin inverse. The application of the scheme on large sparse test matrices alongside the use in preconditioning of linear system of equations will be presented to clarify the contribution of the paper.

1. Introduction

Computing the matrix inverse of nonsingular matrices of higher sizes is difficult and is a time consuming task. Application of higher order algorithms to solve this problem is very desirable. Generally speaking, in wide variety of topics, one must compute the inverse or particularly the generalized inverses to comprehend and realize significant features of the involved problems [1]. An example could be in phased-array radar whereas the target tracking is a recursive prediction correction process, when Kalman filtering is extensively consumed; see [2, 3]. Target equations are modeled explicitly such that the position and velocity and potentially higher derivatives of each measurement are approximated by the track filter as a state vector. The approximated error with the state vector is modeled by taking into account a covariance matrix, which is then used in subsequent computations. To be more precise, this matrix gets updated in each iteration of the track filter. Finding the inverse in the next iteration could make use of the inverse in the present iteration. In this circumstance, fast and efficient iterative algorithms are required.

There are some techniques to tackle this problem, which are basically divided into two main parts: the direct solvers such as Gaussian Elimination with Partial Pivoting (GEPP), which requires a massive load of computations and memory for large scale problems, and the iterative methods of the class of Schulz-type iterations, in which an approximation of the matrix inverse (by using a threshold) can be found per step up to the desired accuracy.

Almost all of the direct methods for matrix inversion require high accuracy in the computations to attain proper solutions as they are not tolerant to errors in the computed matrices. In contrast, iterative method compensates for individual and accumulation of round-off errors as it is a process of successive refinement.

In this paper, we focus on presenting and demonstrating a new iterative method to find approximate inverse matrices as fast as possible with a close attention in reducing the computational time. Toward this goal, a theoretical discussion will also be given to show the behavior of the proposed scheme. An interesting point in the contribution is that it could easily be applied to complex matrices as well as for finding the Drazin inverse.

To clarify the procedure, we now remind of some of the well-known methods in what follows. Perhaps, the most common technique to compute the inverse of a nonsingular complex matrix 𝐴, is the Schulz method given in [4] as follows: 𝑉_𝑛+1=𝑉_𝑛2𝐼−𝐴𝑉_𝑛,𝑛=0,1,2,…,(1) where 𝐼 is the identity matrix with the same dimension that of the matrix 𝐴. The scheme (1) has become popular in the 1980s due to the emerging of parallel machines.

Such iterative methods are sensitive for the initial guess/value (𝑉₀) to start the process and converge to 𝐴⁻¹. In practice, the Schulz-type methods are efficient (specially for structured matrices) but a difficulty lies in the initial approximation of the 𝐴⁻¹. This need was fulfilled by providing some appropriate initial approximations in the literature. For example, Rajagopalan in [5] gave some initial approximations for the inverse matrix by considering different norms as: 𝑉₀=𝐴‖𝐴‖²_∞.(2) In fact, by choosing (2), we attain 𝑉_0∞𝐴⁻¹_∞=‖𝐴‖_∞‖𝐴‖²_∞𝐴⁻¹_∞.(3) Based on the fact that 𝜅(𝐴)=‖𝐴‖‖𝐴⁻¹‖≥‖𝐼‖≥1, we could have 𝑉_0∞𝐴⁻¹_∞=1𝜅_∞(𝐴)≤1,(4) which suggests that by choosing (2), the iterative scheme (1) will be almost always convergent for regular matrices.

Some other ways to choose the initial approximation 𝑉₀ have been listed in Table 1, where 𝐴^𝑇 and 𝐴^∗ are the transpose and the conjugate transpose of the complex matrix 𝐴, respectively, and 𝑁 stands for the size of the square matrix.

Table 1 

Some of the general ways to choose 𝑉₀.

A vast discussion on choosing the initial approximation 𝑉₀ is given in [6]. For instance, Pan and his coworkers discussed the possible ways to reduce the computational load of Schulz-type iteration methods for structured matrices such as Toeplitz or Hankel matrices; see, for example, [7]. To illustrate further, for a symmetric positive definite (SPD) matrix 𝐴, one can choose 𝑉₀ as follows: 𝑉₀=𝐼‖𝐴‖_𝐹.(5)

Another interesting choice is based on [8] for diagonally dominant matrices, which is fruitful when dealing with large scale sparse systems arising from the discretization of differential equations as follows: 𝑉₀1=diag𝑎₁₁,1𝑎₂₂1,…,𝑎_𝑁𝑁,(6) with 𝑎_𝑖𝑖 as the diagonal elements of 𝐴.

Let us now review some of the high order iterative methods for matrix inversion. The perception of the need for higher order methods is the fact that (1) is too slow at the beginning of the process before arriving at the convergence phase for general matrices, and this would increase the computational load of the whole algorithm used for matrix inversion [9].

Li et al. in [10] investigated an iteration of the form 𝑉_𝑛+1=𝑉_𝑛3𝐼−𝐴𝑉_𝑛3𝐼−𝐴𝑉_𝑛,𝑛=0,1,2,…(7) of third-order and also proposed another iterative method for finding 𝐴⁻¹ of the same order as it is given in 𝑉_𝑛+1=𝑉_𝑛1𝐼+2𝐼−𝐴𝑉_𝑛𝐼+2𝐼−𝐴𝑉_𝑛²,𝑛=0,1,2,….(8)

Note that a general procedure for constructing such methods was given in ([11], Chapter 5). Krishnamurthy and Sen provided the following fourth-order method: 𝑉_𝑛+1=𝑉_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛,𝑛=0,1,2,…,(9) in which 𝑌_𝑛=𝐼−𝐴𝑉_𝑛. As another example, a ninth-order method could be presented as 𝑉_𝑛+1=𝑉_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛𝐼+𝑌_𝑛,𝑛=0,1,2,….(10)

In the sequel, we present a new iteration for matrix inversion. Actually, the following section covers our main contribution as a new ninth-order efficient inverse-finding iterative method. We also prove the main theorem therein. Next in Section 3, we discuss the complexity of the iterative methods to theoretically find the most efficient method. In Section 4, we analytically discuss the application of the new algorithm in the computation of the Drazin inverse, which is of interest in numerical analysis. Section 5 applies the suggested iteration in finding robust approximate inverses for large sparse matrices with real or complex entries in details, while the application of the new scheme in preconditioning of the practical problems will be also given. A clear reduction in the computational time to attain the desired accuracy will be observed therein. Finally, in Section 6 our concluding remarks will be furnished.

2. A New Method for Matrix Inversion

A new scheme must be designed by applying an efficient nonlinear (scalar) equation solver to the matrix equation 𝐹(𝑉)=𝑉⁻¹−𝐴. Then, one may obtain an iterative process using proper factorization. In this way, by applying the iterative scheme 𝑌_𝑛=𝑉_𝑛−𝐹𝑉_𝑛𝐹^′𝑉_𝑛,𝑍_𝑛=𝑉_𝑛−𝐹𝑉_𝑛21𝐹^′𝑉_𝑛+1𝐹^′𝑌_𝑛,𝑉_𝑛+1=𝑍_𝑛−𝐹𝑍_𝑛2𝐹^′𝑍_𝑛𝐹2+^′′𝑍_𝑛𝐹𝑍_𝑛𝐹^′𝑍_𝑛²,𝑛=0,1,2,…(11) on the matrix equation 𝑉⁻¹−𝐴=0, we are able to find the fixed-point iteration𝑉_𝑛+11=−8𝑉_𝑛−7𝐼+9𝐴𝑉_𝑛−5𝐴𝑉_𝑛²+𝐴𝑉_𝑛³×12𝐼−42𝐴𝑉_𝑛+103𝐴𝑉_𝑛²−156𝐴𝑉_𝑛³+157𝐴𝑉_𝑛⁴−104𝐴𝑉_𝑛⁵+43𝐴𝑉_𝑛⁶−10𝐴𝑉_𝑛⁷+𝐴𝑉_𝑛⁸.(12) Now, by simplification, we suggest the following matrix iteration: 𝜒_𝑛=−7𝐼+𝜓_𝑛9𝐼+𝜓_𝑛−5𝐼+𝜓_𝑛,𝜗_𝑛=𝜓_𝑛𝜒_𝑛,𝑉_𝑛+11=−8𝑉_𝑛𝜒_𝑛12𝐼+𝜗_𝑛6𝐼+𝜗_𝑛,𝑛=0,1,2,…,(13) wherein 𝜓_𝑛=𝐴𝑉_𝑛,  𝐼 is the identity matrix, and the sequence of iterates {𝑉_𝑛}^𝑛=∞_𝑛=0 converges to 𝐴⁻¹ under some condition.

In numerical mathematics, it is very useful and essential to know the behavior of an approximate method. Therefore, we are about to prove its order of convergence in Theorem 1.

Theorem 1. Let 𝐴=[𝑎_𝑖𝑗]∈ℂ^𝑁×𝑁 be a nonsingular complex matrix. If the initial approximation 𝑉₀ satisfies 𝐼−𝐴𝑉₀<1,(14) then the iterative method (13) converges with ninth order to 𝐴⁻¹.

Proof. We use notations that 𝐸₀=𝐼−𝐴𝑉₀ and subsequently 𝐸_𝑛=𝐼−𝐴𝑉_𝑛. Then, 𝐸_𝑛+1=𝐼−𝐴𝑉_𝑛+1−1=𝐼−𝐴8𝑉_𝑛−7𝐼+𝐴𝑉_𝑛9𝐼+𝐴𝑉_𝑛−5𝐼+𝐴𝑉_𝑛×12𝐼+𝐴𝑉_𝑛−7𝐼+𝐴𝑉_𝑛9𝐼+𝐴𝑉_𝑛−5𝐼+𝐴𝑉_𝑛×6𝐼+𝐴𝑉_𝑛−7𝐼+𝐴𝑉_𝑛9𝐼+𝐴𝑉_𝑛−5𝐼+𝐴𝑉_𝑛−1=𝐼−𝐴8𝑉_𝑛−7𝐼+9𝐴𝑉_𝑛−5𝐴𝑉_𝑛²+𝐴𝑉_𝑛³×12𝐼−42𝐴𝑉_𝑛¹+103𝐴𝑉_𝑛²−156𝐴𝑉_𝑛³+157𝐴𝑉_𝑛⁴−104𝐴𝑉_𝑛⁵+43𝐴𝑉_𝑛⁶−10𝐴𝑉_𝑛⁷+𝐴𝑉_𝑛⁸=18−2𝐼+𝐴𝑉_𝑛³−𝐼+𝐴𝑉_𝑛⁹=18𝐼+𝐼−𝐴𝑉_𝑛³𝐼−𝐴𝑉_𝑛⁹=18𝐼+𝐸_𝑛³𝐸⁹_𝑛=18𝐸⁹_𝑛+3𝐸_𝑛¹⁰+3𝐸_𝑛¹¹+𝐸_𝑛¹².(15)Hence, by taking arbitrary matrix norm on both sides of (15), we attain 𝐸_𝑛+1≤18𝐸_𝑛⁹+3𝐸_𝑛¹⁰+3𝐸_𝑛¹¹+𝐸_𝑛¹².(16) In addition, since ‖𝐸₀‖<1, by relation (16), we obtain that 𝐸₁≤18𝐸₀⁹+3𝐸₀¹⁰+3𝐸₀¹¹+𝐸₀¹²≤𝐸₀⁹<1.(17) Now if we consider ‖𝐸_𝑛‖<1, therefore 𝐸_𝑛+1≤18𝐸_𝑛⁹+3𝐸_𝑛¹⁰+3𝐸_𝑛¹¹+𝐸_𝑛¹²≤𝐸_𝑛⁹.(18) Using mathematical induction, we obtain 𝐸_𝑛+1≤𝐸_𝑛⁹,𝑛≥0.(19) Furthermore, we get that 𝐸_𝑛+1≤𝐸_𝑛⁹≤⋯≤𝐸₀^9𝑛+1.(20) That is, 𝐼−𝐴𝑉_𝑛→0, when 𝑛→∞ and 𝑉_𝑛⟶𝐴⁻¹,as𝑛⟶∞.(21) Thus, the new method (13) 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 nine. Let 𝜀_𝑛 denotes the error matrix 𝜀_𝑛=𝐴⁻¹−𝑉_𝑛; afterwards 𝐴𝜀_𝑛=𝐼−𝐴𝑉_𝑛=𝐸_𝑛.(22) The identity (22) in conjunction with (15) implies that 𝐴𝜀_𝑛+1=18𝐴𝜀_𝑛⁹+3𝐴𝜀_𝑛¹⁰+3𝐴𝜀_𝑛¹¹+𝐴𝜀_𝑛¹².(23) Therefore, using invertibility of 𝐴, it follows immediately that 𝜀_𝑛+1=18𝜀_𝑛𝐴𝜀_𝑛⁸+3𝜀_𝑛𝐴𝜀_𝑛⁹+3𝜀_𝑛𝐴𝜀_𝑛¹⁰+𝜀_𝑛𝐴𝜀_𝑛¹¹.(24) By taking any subordinate norm of (24), we obtain 𝜀_𝑛+1≤18‖𝐴‖⁸+3‖𝐴‖⁹𝜀_𝑛+3‖𝐴‖¹⁰𝜀_𝑛²+‖𝐴‖¹¹𝜀_𝑛³𝜀_𝑛⁹.(25) Consequently, it is proved that the iterative formula (13) converges to 𝐴⁻¹, and the order of this method is at least nine.

The Schulz-type iterations are strongly numerically stable, that is, they have the self-correcting characteristic and are essentially based upon matrix multiplication per an iterative step. Multiplication is effectively parallelizable for structured matrices represented in compressed form.

The iterative scheme (13) could efficiently be combined with sparse techniques in order to reduce the computational load of matrix-by-matrix multiplications per step. We should also point out that even if the matrix 𝐴 is singular, the Schulz-type methods, including the scheme (13), converge to the Moore-Penrose inverse using a proper initial matrix. A full discussion on this feature of this type of iterative methods has been given in [12].

Note that {𝑉_𝑛}^𝑛=∞_𝑛=0 produced from (13), under a certain condition (when 𝐴𝑉₀=𝑉₀𝐴), may be applied not only to the left preconditioned linear system 𝑉_𝑛𝐴𝑥=𝑉_𝑛𝑏 but also to the right preconditioned linear system 𝐴𝑉_𝑛𝑦=𝑏, where 𝑦=𝑉_𝑛𝑥. In fact, an important application of the new method (13) is in preconditioning of the linear system of equations. Practically, experimental results in Section 5 will show that the preconditioner obtained from (13) may lead to nicely clustered eigenvalue distributions of the preconditioned matrices and, hence, results in fast convergence of the preconditioned Krylov subspace iteration methods, such as GMRES and BiCGSTAB for solving some classes of large sparse system of linear equations.

3. Complexity of the Methods

Let us consider the computational complexity of the existing iterative processes (1), (7), (8), (9), (10), and (13), since they are all convergent to 𝐴⁻¹ in the same condition. From a theoretical analysis, and by assuming a uniform cost for the arithmetic operations, typical of the floating point computations, we consider the inverse-finder informational efficiency index. It uses 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 𝜌IIEI=𝜅.(26)

Hence, a favorable method in theoretical point of view must reach an order 𝜌 with fewer matrix multiplications 𝜅, (i.e., 𝜅≤𝜌).

In Table 2, we furnish a comparison on the order along the number of matrix multiplications, rate of convergence, and the index (26) for different methods. The results show that the new established method in Section 2 is better than the others. In fact, by comparing these results, one can see that the iterative process (13) 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 2 

Comparison of the computational complexity for different methods.

4. Application in Finding the Drazin Inverse

In 1958, Drazin in [13] introduced a different kind of generalized inverse in associative rings and semigroups that does not have the reflexivity property but commutes with the element. The importance of this kind of inverse and its computation was later expressed away fully by Wilkinson in [14]. This was the motivation of many authors to develop direct or iterative methods for this important problem; see, for example, [1].

Definition 2. The smallest nonnegative integer 𝑘, such that rank(𝐴^𝑘+1)=rank(𝐴^𝑘), is called the index of 𝐴 and denoted by ind(𝐴).

Definition 3. Let 𝐴 be an 𝑁×𝑁 complex matrix; the Drazin inverse of 𝐴, denoted by 𝐴^𝐷, is the unique matrix 𝑉 satisfying 1^𝑘𝐴^𝑘𝑉𝐴=𝐴^𝑘,(2)𝑉𝐴𝑉=𝑉,(5)𝐴𝑉=𝑉𝐴,(27) where 𝑘=ind(𝐴) is the index of 𝐴.

Note that if ind(𝐴)=1, the matrix 𝑉 is called the group inverse of 𝐴. Also, if 𝐴 is nonsingular, then it is easily seen that ind(𝐴)=0 and 𝐴^𝐷=𝐴⁻¹. Note that the idempotent matrix 𝐴𝐴^𝐷 is the projector on ℛ(𝐴^𝑘) along 𝒩(𝐴^𝑘), where ℛ(𝐴^𝑘) denotes the range of 𝐴^𝑘 and 𝒩(𝐴^𝑘) is the null space of 𝐴^𝑘.

Wei in [15] proved that the general solution of the square singular linear system 𝐴𝑥=𝑏 can be obtained using the Drazin inverse as 𝑥=𝐴^𝐷𝑏+(𝐼−𝐴𝐴^𝐷)𝑧, where 𝑧∈ℛ(𝐴^𝑘−1)+𝒩(𝐴).

In 2004, Li and Wei in [16] proved that the matrix method of Schulz (1) can be used for finding the Drazin inverse of square matrices both possessing real or complex spectra. They proposed the following initial matrix: 𝑉₀=𝑊₀=𝛼𝐴^𝑙,𝑙≥ind(𝐴)=𝑘,(28) where the parameter 𝛼 must be chosen so that the condition ‖𝐼−𝐴𝑉₀‖<1 is satisfied. Using the initial matrix of the form (28) yields to a matrix method for finding the famous Drazin inverse with quadratical convergence.

As a consequence, we could present the iterative method of the form (13) with ninth order of convergence for finding the Drazin inverse, where the initial approximation is chosen as 𝑉₀=𝑊₀=2Tr𝐴^𝑘+1𝐴^𝑘,(29) wherein Tr(⋅) stands for the trace of an arbitrary square matrix.

In what follows, we use the following auxiliary results.

Proposition 4 (see [17]). Let 𝑀∈ℂ^𝑁×𝑁 and 𝜀>0 be given. There is at least one matrix norm ‖⋅‖ such that 𝜌(𝑀)≤‖𝑀‖≤𝜌(𝑀)+𝜖,(30) where 𝜌(𝑀) denotes the spectral radius of 𝑀.

Proposition 5 (see [18]). If 𝑃_𝐿,𝑀 denotes the projector on a space 𝐿 along a space 𝑀, then(i)𝑃_𝐿,𝑀𝑄=𝑄 if and only if ℛ(𝑄)⊆𝐿;(ii)𝑄𝑃_𝐿,𝑀=𝑄 if and only if 𝒩(𝑄)⊇𝑀.

Theorem 6. Let 𝐴∈ℂ^𝑁×𝑁 be singular square matrix. Also, suppose that the initial approximation 𝑊₀ is chosen by means of (28). Then the sequence {𝑊_𝑛}^𝑛=∞_𝑛=0 defined by the iterative method (13) satisfies the following error estimate when finding the Drazin inverse: 𝐴^𝐷−𝑊_𝑛𝐴^𝐷≤𝐼−𝐴𝑊₀^9𝑛.(31)

Proof. Consider the notation 𝐹₀=𝐼−𝐴𝑊₀ and subsequently the residual matrix as 𝐹_𝑛=𝐼−𝐴𝑊_𝑛. Then similarly as in (15), we get𝐹_𝑛+1=𝐼−𝐴𝑊_𝑛+1−1=𝐼−𝐴8𝑊_𝑛−7𝐼+𝐴𝑊_𝑛9𝐼+𝐴𝑊_𝑛−5𝐼+𝐴𝑊_𝑛×12𝐼+𝐴𝑊_𝑛−7𝐼+𝐴𝑊_𝑛9𝐼+𝐴𝑊_𝑛−5𝐼+𝐴𝑊_𝑛×6𝐼+𝐴𝑊_𝑛−7𝐼+𝐴𝑊_𝑛9𝐼+𝐴𝑊_𝑛−5𝐼+𝐴𝑊_𝑛=18𝐼+𝐼−𝐴𝑊_𝑛³𝐼−𝐴𝑊_𝑛⁹=18𝐹⁹_𝑛+3𝐹_𝑛¹⁰+3𝐹_𝑛¹¹+𝐹_𝑛¹².(32)By taking an arbitrary matrix norm on both sides of (32), we attain 𝐹_𝑛+1≤18𝐹_𝑛⁹+3𝐹_𝑛¹⁰+3𝐹_𝑛¹¹+𝐹_𝑛¹².(33) In addition, since ‖𝐹₀‖<1, by relation (33), we obtain that ‖𝐹₁‖≤(1/8)(‖𝐹₀‖⁹+3‖𝐹₀‖¹⁰+3‖𝐹₀‖¹¹+‖𝐹₀‖¹²)≤‖𝐹₀‖⁹. Similarly, ‖𝐹_𝑛+1‖≤(1/8)(‖𝐹_𝑛‖⁹+3‖𝐹_𝑛‖¹⁰+3‖𝐹_𝑛‖¹¹+‖𝐹_𝑛‖¹²)≤‖𝐹_𝑛‖⁹. Using mathematical induction, we obtain ‖𝐹_𝑛+1‖≤‖𝐹_𝑛‖⁹ for each 𝑛≥0. This implies that 𝐹_𝑛⁹≤𝐹₀^9𝑛,𝑛≥0.(34)
Since 𝑊₀ is chosen as in (28), ℛ(𝑊₀)⊆ℛ(𝐴^𝑘) immediately follows. Further, using this fact in conjunction with (12), which implies that ℛ(𝑊_𝑛)⊆ℛ(𝑊_𝑛−1), we conclude that ℛ𝑊_𝑛⊆ℛ𝐴^𝑘,𝑛≥0.(35) Similarly, if we rewrite (13) in the form 𝑊_𝑛+11=−8−7𝐼+9𝑊_𝑛𝐴−5𝑊_𝑛𝐴²+𝑊_𝑛𝐴³×12𝐼−42𝐴𝑊_𝑛+103𝑊_𝑛𝐴²−156𝑊_𝑛𝐴³+157𝑊_𝑛𝐴⁴−104𝑊_𝑛𝐴⁵+43𝑊_𝑛𝐴⁶−10𝑊_𝑛𝐴⁷+𝑊_𝑛𝐴⁸𝑊_𝑛.(36) it is not difficult to verify that 𝒩𝑊_𝑛⊇𝒩𝐴^𝑘,𝑛≥0.(37) Now, an application of the well-known results from [12] 𝐴𝐴^𝐷=𝐴^𝐷𝐴=𝑃_ℛ(𝐴^𝑘_),𝒩(𝐴^𝑘₎(38) in conjunction with Proposition 5 and (35), (37), immediately follows 𝑊_𝑛𝐴𝐴^𝐷=𝑊_𝑛=𝐴^𝐷𝐴𝑊_𝑛,𝑛≥0.(39) Therefore, the error matrix 𝛿_𝑛=𝐴^𝐷−𝑊_𝑛 satisfies 𝛿_𝑛=𝐴^𝐷−𝑊_𝑛=𝐴^𝐷−𝐴^𝐷𝐴𝑊_𝑛=𝐴^𝐷𝐼−𝐴𝑊_𝑛=𝐴^𝐷𝐹_𝑛.(40) From the last identity and (34) we have 𝛿_𝑛=𝐴^𝐷𝐹_𝑛≤𝐴^𝐷𝐹₀^9𝑛,(41) which is a confirmation of (31).

The following result is a consequence of Theorem 6.

Corollary 7. If the conditions of Theorem 6 are satisfied and the initial iteration 𝑊₀ is chosen such that 𝐹₀=𝐼−𝐴𝑊₀<1,(42) the iterative method (13) converges to 𝐴^𝐷.

Therefore, our goal is to find initial approximations 𝑊₀ satisfying (42). In accordance with Proposition 4, 𝑊₀ must satisfy the following inequality to ensure the convergence in the Drazin inverse case: max_{1≤𝑖≤𝑡}1−𝜆_𝑖𝐴𝑊₀<1,(43) where rank(𝐴𝑊₀)=𝑡 and 𝜆_𝑖(𝐴𝑊₀),𝑖=1,…,𝑡 are eigenvalues of 𝐴𝑊₀.

Theorem 8. Let 𝐴∈ℂ^𝑁×𝑁 be a singular square matrix with ind(𝐴)=𝑘, and the sequences 𝐹_𝑛, 𝛿_𝑛 are defined as in Theorem 6. The sequence {𝑊_𝑛}^𝑛=∞_𝑛=0 defined by the iterative method (13) converges with ninth order to 𝐴^𝐷 if the initial approximation 𝑊₀ is in accordance with (42).

Proof. Now, by considering 𝛿_𝑛=𝐴^𝐷−𝑊_𝑛 as the error matrix for finding the Drazin inverse, we have 𝐴𝛿_𝑛+1=𝐴𝐴^𝐷−𝐴𝑊_𝑛+1=𝐴𝐴^𝐷−𝐼+𝐼−𝐴𝑊_𝑛+1=𝐴𝐴^𝐷−𝐼+𝐹_𝑛+1.(44) Taking into account (32) and using elementary algebraic transformations, we further derive 𝐴𝛿_𝑛+1=𝐴𝐴^𝐷1−𝐼+8𝐹⁹_𝑛+3𝐹_𝑛¹⁰+3𝐹_𝑛¹¹+𝐹_𝑛¹²=18𝐹⁹_𝑛+3𝐹_𝑛¹⁰+3𝐹_𝑛¹¹+𝐹_𝑛¹²+8𝐴𝐴^𝐷=1−𝐼8𝐹⁹_𝑛+𝐴𝐴^𝐷−𝐼+3𝐹_𝑛¹⁰+𝐴𝐴^𝐷𝐹−𝐼+3_𝑛¹¹+𝐴𝐴^𝐷−𝐼+𝐹_𝑛¹²+𝐴𝐴^𝐷−𝐼.(45) Now, using the idempotent property (𝐼−𝐴𝐴^𝐷)^𝑡=(𝐼−𝐴𝐴^𝐷), 𝑡≥1, and the following consequence of (39) 𝐼−𝐴𝐴^𝐷𝐴𝛿_𝑛=𝐼−𝐴𝐴^𝐷𝐴𝐴^𝐷−𝑊_𝑛=𝑊_𝑛−𝐴𝐴^𝐷𝑊_𝑛=0,(46) we obtain for each 𝑡≥1 the following: 𝐹_𝑛^𝑡+𝐴𝐴^𝐷−𝐼=𝐼−𝐴𝑊_𝑛^𝑡+𝐴𝐴^𝐷=−𝐼𝐼−𝐴𝐴^𝐷+𝐴𝐴^𝐷−𝐴𝑊_𝑛^𝑡+𝐴𝐴^𝐷=−𝐼𝐼−𝐴𝐴^𝐷+𝐴𝛿_𝑛^𝑡+𝐴𝐴^𝐷−𝐼=𝐼−𝐴𝐴^𝐷+𝐴𝛿_𝑛^𝑡+𝐴𝐴^𝐷=−𝐼𝐴𝛿_𝑛^𝑡.(47) From (47) and (45), 𝐴𝛿_𝑛+1=18𝐴𝛿_𝑛⁹+3𝐴𝛿_𝑛¹⁰+3𝐴𝛿_𝑛¹¹+𝐴𝛿_𝑛¹²(48) Therefore, 𝐴𝛿_𝑛+1≤18𝐴𝛿_𝑛⁹+3𝐴𝛿_𝑛¹⁰+3𝐴𝛿_𝑛¹¹+𝐴𝛿_𝑛¹²≤‖𝐴𝛿_𝑛‖⁹.(49)
Finally, applying (39), it is now easy to find the error inequality of the new scheme (13) using (49) and the second condition of (27), when finding the Drazin inverse, as follows: 𝛿_𝑛+1=𝑊_𝑛+1−𝐴^𝐷=𝐴^𝐷𝐴𝑊_𝑛+1−𝐴^𝐷𝐴𝐴^𝐷=𝐴^𝐷𝐴𝑊_𝑛+1−𝐴𝐴^𝐷≤𝐴^𝐷𝐴𝛿_𝑛+1≤𝐴^𝐷‖𝐴‖⁹𝛿_𝑛⁹.(50)
Therefore, since (42) is satisfied, from (31) follows 𝛿_𝑛→0. Furthermore, the inequalities in (50) immediately lead to the conclusion that 𝑊_𝑛→𝐴^𝐷 as 𝑛→+∞ with the ninth order of convergence.

5. Numerical Aspects

Using the programming package Mathematica 8 [19] in this section, we apply our iterative method on some practical numerical tests and compare it with the existing methods in order to manifest the applicability and the consistency of numerical results with the theoretical aspects illustrated in Sections 2–4.

For numerical comparisons in this section, we have used the methods (1), (7), (8), (9), (10), and (13) denoted by “Schulz”, “Li et al. I”, “Li et al. II”, “KMS4”, “KMS9,” and the “Proposed method”, respectively. We have carried out the numerical tests with machine precision on a computer with Pentium 4. In fact, the computer characteristics are Microsoft Windows XP Intel(R), Pentium(R) 4 CPU, 3.20 GHz with 4 GB of RAM. In all computations, the running time in seconds using AbsoluteTiming[] was attained.

In sparse-matrix algebra, the iteration methods such as (1), (7), (8), (9), (10) and (13) should be coded in sparse form using some well-known commands such as SparseArray[] to reduce the computational burden and preserve the sparsity feature of the approximate inverse per computing step. This is done herein along with a threshold by applying Chop[] on each approximate inverse 𝑉_𝑖.

Now, we apply the above inverse finders for finding the approximate inverses of the following large sparse matrices.

Test Problem 1. In this test, 10 large sparse random complex matrices of the dimension 5000 are considered as shown in Algorithm 1.

n  = 5000; number = 10; SeedRandom [123] Table[A[j]  = SparseArray[{Band[{−100, 1100}]  -> RandomReal[20],      Band[{1, 1}]  -> 2., Band[{1000, −50}, {n −  20, n −  25}]  -> {2.8, RandomReal[]  + I},      Band[{600, 150}, {n −  100, n −  400]  -> {-RandomReal[3], 3. + 3 I}}, {n, n},           0.], {j, number}];

Algorithm 1 

Note that 𝐼=√−1. In this test, the stopping criterion is ‖𝑉_𝑛+1−𝑉_𝑛‖₁≤10⁻⁶ and the maximum number of iterations allowed is set to 75. Note that in this test the initial choice has been constructed using different available ways as discussed in Table 1.

The results of comparisons for this test have been presented in Figures 1, 2, 3, and 4. In Figure 1, we show the number of iterations required by different methods to attain the desired accuracy when 𝑉₀=𝐴/‖𝐴‖²₁. As it is obvious that higher order methods require lower number of iterations to converge, we then put our focus on the computational time needed to satisfy the desired tolerance using three different initial approximations 𝑉₀=𝐴/‖𝐴‖²₁, 𝑉₀=𝐴/‖𝐴‖²_∞, and 𝑉₀=𝐴/‖𝐴‖²_𝐹. As could be observed for all the three forms of the initial guesses and in all the 10 test problems (given in Figures 2–4), our iterative method (13) beats the other existing schemes, which immediately follows the theoretical results of Table 2.

Figure 1 

Comparison of the number of evaluations for solving the Test Problem 1 using 𝑉0=𝐴/‖𝐴‖21.

Figure 2 

Comparison of the elapsed time for solving the Test Problem 1 using 𝑉0=𝐴/‖𝐴‖21.

Figure 3 

Comparison of the elapsed time for solving the Test Problem 1 using 𝑉0=𝐴/‖𝐴‖2∞.

Figure 4 

Comparison of the elapsed time for solving the Test Problem 1 using 𝑉0=𝐴/‖𝐴‖2𝐹.

The attained results have reverified the robustness of the proposed iterative method (13) by a clear reduction in the number of iterations and the elapsed time.

The practical application of the new scheme (13) falls within many problems as discussed in Section 1. For instance, in solving second-kind integral equations by Wavelet-like approach, the problem will be reduced to find that the inverse of a large sparse matrix possesses a sparse inverse as well [20]. In the rest of this section, we apply our new iterative method as a robust technique to produce accurate preconditioners for accelerating modern iterative solvers such as GMRES or BiCGSTAB for solving large scale sparse linear systems; see, for example, [21].

Test Problem 2. Consider solving the following Boundary Value Problem (BVP) using discretization. In such problems, and in order to capture all the shocks and the behavior of the solution, a very fine grid of points is required in 𝔽𝔻 approach: 𝑦^′′=3𝑦−2𝑦^′[],𝑦(𝑎),𝑡∈𝑎,𝑏=𝑒³,𝑦(𝑏)=𝑒⁻³,(51) where the domain of the solution is 𝑎=0 and 𝑏=2. To solve (51) using 3-point 𝔽𝔻 discretization, we assume that in the grid points 𝑡_𝑖, the exact solution is denoted by 𝑦(𝑡_𝑖) and the approximate solution is defined by 𝑤_𝑖, for any 𝑖=0,1,…,𝑛, 𝑛+1.

Thus, (51) can be written in discretized form in MATHEMATICA language (see Algorithm 2), wherein 𝑛 is the number of grid points which result in an 𝑛×𝑛 sparse linear system of equations. In this test, we consider the tolerance as 10⁻⁸ for the final solution of the linear systems. The results of solving this system “without preconditioning” and “with left preconditioning” are given in Tables 3 and 4. Note that in Tables 3 and 4, for example, PBiCGSTAB-(8)-II stands for the left preconditioned linear system (𝑉𝐴𝑥=𝑉𝑏) applying the preconditioner (approximative inverse) obtained by the scheme (8) after 2 iterations, while it is solved by the iterative solver BiCGSTAB.

Table 3 

Comparison of the computational time in solving the linear system resulting of discretization of (51) when 𝑛=1500.

Table 4 

Comparison of the computational time in solving the linear system resulting of discretization of (51) when 𝑛=2000.

h  = (2. −  0)/(n + 1); Subscript[t, 0]  = 0; Subscript[w, 0]  =          E∧3; Subscript[w, n + 1]  = Exp[−3]; Table[N[(Subscript[w, i −  1]  −  2 Subscript[w, i]  + Subscript[w, i + 1])/       h∧2 −  3 Subscript[w, i]  +      2 (Subscript[w, i + 1]  −  Subscript[w, i −  1])/(2 h) == 0], {i, 1, n}];

Algorithm 2 

The numerical results clearly support the efficiency of the method (13). A clear reduction in the elapsed time is observable. Even for the case of GMRES solver which had failed (not convergent after 1500 iterations to the considered tolerance), a simple preconditioner obtained by the new method (13) significantly improved the problem. Note that in this test, we have constructed 𝑉₀ for the compared methods using (6).

Test Problem 3. The aim of this example is to apply the discussions of Section 4, for finding the Drazin inverse of the following square matrix (taken from [16]): 𝐴=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,(52) with 𝑘=ind(𝐴)=3. To simplify the process, we write a general code in the programming package Mathematica for the iterative process (13), to find the (pseudo-)inverse or the Drazin inverse of arbitrary matrices (see Algorithm 3).

PM[X_]  := With[{Id = SparseArray[{{i_, i_}  -> 1.}, {n, n}]],          X1 = A.X; X2 = −7 Id + X1.(9 Id + X1.(−5 Id + X1)); X3 = X1.X2;          (−1/8) X.X2.(12 Id + X3.(6 Id + X3))           ]; InitialMatrix[A_]  := 1/SingularValueList[A, 1] [[1]]∧2  ConjugateTranspose[A]; InitialDrazin[A_]  := 2/Tr[MatrixPower[A, k + 1]]  MatrixPower[A, k]; DrazinInverse[A_, tolerance_]  := If[k == 0,       Module[{X0 = InitialMatrix[A]}, FixedPoint[(PM[#]  &), X0,         SameTest -> (Norm[#1 - #2, Infinity]  <= tolerance  &)]],       Module[{X0 = InitialDrazin[A], FixedPoint[(PM[#]  &), X0,         SameTest -> (Norm[#1 −  #2, Infinity]  <= tolerance  &)]]       ];

Algorithm 3 

The two-argument function DrazinInverse[A_,tolerance_] takes the arbitrary matrix 𝐴 and the tolerance from the user to obtain its Drazin inverse by knowing the index 𝑘. In this case, by choosing tolerance=10⁻⁸, we obtain𝐴^𝐷=^{0.25−0.250.0.0.0.0.0.0.0.0.0.1.251.250.0.0.0.0.0.0.0.0.0.−1.66406−0.9921870.25−0.250.0.0.0.−0.0625−0.06250.0.15625−1.19531−0.679687−0.250.250.0.0.0.−0.06250.18750.68751.34375−2.76367−1.04492−1.875−1.25−1.251.251.251.251.484382.578133.320316.64063−2.76367−1.04492−1.875−1.25−1.251.251.251.251.484382.578134.570318.51563}_{14.10946.300786.6253.3755.−3.−5.−5.−4.1875−8.5−10.5078−22.4609−19.3242−8.50781−9.75−5.25−7.54.57.57.56.37512.562515.976633.7891−0.625−0.31250.0.0.0.0.0.0.25−0.25−0.875−1.625−1.25−0.93750.0.0.0.0.0.−0.250.25−0.875−1.6250.0.0.0.0.0.0.0.0.0.1.251.250.0.0.0.0.0.0.0.0.0.−0.250.25}.(53)

Checking the conditions of Definition 3 yields to 𝐴^𝑘+1𝐴^𝐷−𝐴^𝑘_∞=1.48415×10⁻¹²,𝐴^𝐷𝐴𝐴^𝐷−𝐴^𝐷_∞=1.20264×10⁻¹⁰,𝐴𝐴^𝐷−𝐴^𝐷𝐴_∞=8.93836×10⁻¹¹,(54) which supports the theoretical discussions.

6. Conclusion

In many engineering applications, extracting the diagonal or the whole entries of the inverse of a given matrix (basically large and sparse) is an important part of the computation, for example, in electronic structure calculation and especially for models based on effective one-electron Hamiltonians such as the tight-binding models, or in modern statistics, and thus, developing high-order efficient Schulz-type methods is of practical interest.

In this paper, we have studied a high-order iteration method (13) for matrix inversion. Convergence analysis of our iterative algorithm has been studied and a discussion on the choice of the initial value in order to start the process and preserve the convergence order has been given. We also discussed that under what conditions a new method could be applied for finding the Drazin inverse of square matrices having real or complex spectra.

As a result, the total time consuming of the suggested iteration (13) is remarkably low in contrast with the existing methods of this type in the case of constructing approximate inverses and in preconditioning.

Working on the extension of the proposed method (13) for interval matrix inversion [22] can be considered as future works in this field of study.

Conflict of Interests

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

Acknowledgment

The second author gratefully acknowledges support from the Research Project 174013 of the Serbian Ministry of Science.