# An Improved Computationally Efficient Method for Finding the Drazin Inverse

Authors: Haifa Bin Jebreen; Yurilev Chalco-Cano
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

Drazin inverse is one of the most significant inverses in the matrix theory, where its computation is an intensive and useful task. The objective of this work is to propose a computationally effective iterative scheme for finding the Drazin inverse. The convergence is investigated analytically by applying a suitable initial matrix. The theoretical discussions are upheld by several experiments showing the stability and convergence of the proposed method.


---

## Body

## 1. Introductory Notes

Methods for calculating the generalized inverses are a topic of current investigation in computational mathematics (see, e.g., [1, 2]). Enormous amount of work in the topic of generalized inverses has been done during the past 63 years. E.H. Moore (1920) was a pioneer in the area of a generalized algebraic matrix inverse. But intense activities started since 1955 [3, 4]. Note that R. Penrose (1955, 1956), C.R. Rao and S.K. Mitra (1971), and G.H. Golub and W. Kahan (1965) are just a couple of names who have contributed significantly in the area of generalized matrix inverse and its applications leaving very little scope for others to contribute in this area.The terms “the minimum-norm least squares inverse,” “the Moore-Penrose inverse,” and “the pseudo-inverse” of a matrix are much more used than the term “Drazin inverse.” This unique inverse is globally used for solving linear systems and linear optimization problems.The fundamental partitioning method of Greville for calculating generalized inverses was discussed in [3]. This method demands many operations and clearly includes more round-off errors. A lot of computational schemes for computing the Drazin inverse lack numerical stability or slow convergence and accordingly it is requisite to study and propose new iterations for this aim.It is known that cumulative rounding errors must be canceled completely, which is only doable via “symbolic computations.” For such a case, variables are handled in the exact form, yielding no loss of precision in the process of computation. However, once the dimension of the input matrix is big, the calculation of its Drazin inverse via symbolic methods would be so time-consuming. This encouraged researchers to recommend computational stable and fast iterative methods for such a purpose; see, e.g., [5–7].The author in [8] investigated that, in case of having the following square singular linear system(1)Ax=b,the general solution can be represented as follows:(2)x=ADb+I-AADz,using the Drazin inverse, where z∈R(Ak-1)+N(A).LetCN×N and CrN×N indicate the collection of N×N complex matrices and ones having rank r, respectively. Furthermore, with A∗, R(A), rank(A), and N(A), we indicate the conjugate transpose, the range, the rank, and the null space of the matrix A∈CN×N, respectively.In the study of associative rings and semigroups, Drazin in the fundamental work [9] showed a different type of generalized inverse that does not possess the reflexivity feature while commuting with the entry/element. The significance of this type of inverse and its calculation was then completely provided by Wilkinson in [10].Before continuing, it is required to recall several useful definitions in what follows.Definition 1.
The smallest nonnegative integerk, so that(3)rankAk+1=rankAk,is named as the index of A and it is shown by ind(A).Definition 2.
Assume thatA is an n×n complex matrix; then the Drazin inverse AD of the matrix A is a unique matrix V that satisfies the following conditions(4)1kAkVA=Ak,2VAV=V,5AV=VA,wherein k=ind(A) is the matrix index of A.Recal that onceind(A)=1, then V is named as the group inverse of A. In addition, if A is nonsingular, thence it is straightforward to observe that ind(A)=0, and AD=A-1.Note that the idempotent matrixAAD is the projector on R(Ak) alongside N(Ak), whereas R(Ak) denotes the range of Ak and N(Ak) is the null space of Ak. Also, if A is nilpotent, then AD=0; see [11–14] for more.Matrix iterations are sensitive upon the choice of the first approximation (V0) to initiate the scheme and observe the convergence to AD. Practically speaking, the iterative methods (like the Schulz-type iterations) are effective (particularly for sparse matrices having sparse inverses) but an issue lies in the initial value. However, this requirement was lifted by furnishing several suitable initial values in the literature. A vast discussion on choosing the initial choice V0 is given in the work [15].In this work, we concentrate on proposing and investigating a novel matrix iteration for calculating the Drazin inverse quickly and efficiently with a clear concentration on decreasing the elapsed CPU times in contrast to several well-known competitors in the literature. To do this, an analytical discussion will be furnished to manifest the behavior of the presented method. It is proven that the novel approach has a high convergence order using fewer number of matrix-matrix multiplications; i.e., it has a better computational efficiency index.After a short introduction in Section1 and a brief review on the most common iterative schemes for calculating the Drazin inverse in Section 2, we propose a novel iterative method for computing the Drazin inverse. Section 3 unfolds our contributed method as a novel high order efficient method. Section 4 furnishes a discussion regarding its convergence rate. Next in Section 5, we study the complexity of the iterative methods to analytically select the best iterative expression. In Section 6, we numerically investigate the application and usefulness of the novel scheme in the calculation of the Drazin inverse. A clear decrease of the elapsed CPU time will be seen therein. Ultimately in Section 7, summary and comments are brought forward.


## 2. The Literature

One of the fundamental techniques to calculate the inverse of a nonsingular complex matrixA is the Schulz method given in [16] as follows:(5)Vn+1=Vn2I-AVn,n=0,1,2,…,wherein I is the unit matrix of the same dimension as A.In 2004, Li and Wei in [17] proved that the matrix method of Schulz (5) could be applied for calculating the Drazin inverse of square matrices having real or complex eigenvalues. They proposed the following initial matrix(6)V0=V0=αAl,l≥indA=k,wherein the parameter α should be selected such that the following criterion holds(7)I-AV0<1.Using the initial matrix of the form (6) along with (5) results in a quadratically convergent iterative scheme for calculating the well-known Drazin inverse.Let us now briefly review some of the higher order iteration schemes for such a purpose. The notion of the need for efficient schemes is the fact that (5) is slow at its initial stage of iterates, and this would increase the computational burdensome of the scheme applied for matrix inversion.Li et al. in [18] studied and discussed an iterative method in the following formulation(8)Vn+1=Vn3I-AVn3I-AVn,n=0,1,2,…,with cubic convergence rate and also presented another scheme for calculating A-1 (and similarly for AD using a suitable initial value) of the same order as it is provided in the coming formula:(9)Vn+1=VnI+12I-AVnI+2I-AVn2,n=0,1,2,….Recall that a general way for deriving similar iterations was furnished in [19, Chapter 5]. In fact, Krishnamurthy and Sen proposed the following quartically-convergent scheme:(10)Vn+1=VnI+RnI+RnI+Rn,n=0,1,2,…,in which Rn=I-AVn. As another instance, a ninth-order matrix iteration can be written as follows:(11)Vn+1=VnI+RnI+RnI+RnI+RnI+RnI+RnI+RnI+Rn,n=0,1,2,….Generally speaking, applying Schröder’s general iteration (often named as Schröder-Traub’s sequence [20]) to the nonlinear matrix equation(12)AV=I,one can obtain the following scheme [21]:(13)Vn+1=VnI+Rn+Rn2+⋯+Rnm-1=VnI+RnI+Rn⋯+Rn⋯,n=0,1,2,…,of order m, needing m Horner’s matrix products, whereas Rn=I-AVn.


## 3. Proposing a New Method

Now, we propose a high order scheme which is computationally efficient in terms of the number of matrix-matrix multiplications. Let us first consider a tenth-order method using (13) as follows:(14)Vn+1=VnI+Rn+Rn2+⋯+Rn9.To improve the performance of this scheme, we factorize (14) as much as possible so as to decrease the number of matrix-matrix products. First we can obtain(15)Vn+1=VnI+RnI+Rn+Rn2I+Rn3+Rn6,which includes seven matrix-matrix multiplications.By further factorizations and simplifications, we can propose the following iterative method:(16)Vn+1=VnI+RnI+αRn2+Rn4I+βRn2+Rn4,with only 6 matrix-matrix products per cycle where(17)α=121-5,β=121+5.Obtainingα and β so as to reduce the matrix matrix products is new and not easy. In fact, we should solve a system of equations in symbolic environment so as to do such a task.The Schulz-type iterations such as (16) are numerically stable; i.e., they have the self-correcting characteristic and are exactly based upon matrix multiplication per cycle. Multiplication is efficiently parallelizable for special matrices. The method (16) could be mixed with sparse-saving techniques so as to decrease the burdensome of matrix-matrix products per cycle.We can now apply (16) with tenth convergence rate for calculating the Drazin inverse when the first value is selected as follows:(18)V0=2TrAk+1Ak,or(19)V0=12A2k+1Ak,wherein Tr(·) stands for the trace of an arbitrary square matrix with index k.Before providing the main theorems regarding the convergence analysis and the rate of convergence of the proposed method, we recall the following lemmas.Proposition 3 (see [22]).
LetM∈Cn×n and ε>0 be given. There is at least one matrix norm · so that(20)ρM≤M≤ρM+ϵ,where ρ(M) indicates the set of all eigenvalues of M (in the maximum of absolute value sense).Proposition 4 (see [22]).
IfPL,M indicates the projector on a space L along a space M, then(i)
PL,MQ=Q iff R(Q)⊆L;(ii)
QPL,M=Q iff N(Q)⊇M.


## 4. Convergence Analysis

Theorem 5.
Assume thatA∈CN×N is a singular square matrix. In addition, let the initial value V0 be selected via (6) or (18). Thence, the matrices {Vn}n=0n=∞ generated via (16) satisfy the following error estimate for calculating the Drazin inverse:(21)AD-VnAD≤I-AV010n.Also the convergence speed is ten.Proof.
Let us write(22)Rn+1=I-AVn+1=I-AVnI+RnI+αRn2+Rn4I+βRn2+Rn4=I-AVnI+RnI+Rn+Rn2I+Rn3+Rn6=I-AVnI+RnI+Rn+Rn2I+Rn3+Rn6=I-AVn10=Rn10,wherein α and β are defined in (17). Applying an arbitrary matrix norm on relation (22), it is possible to write down(23)Rn+1≤Rn10.Since V0 is chosen as in (6) or (18), it follows that(24)RV0⊆RAk.This could now state that(25)RVn⊆RVn-1.Thus, we can conclude that(26)RVn⊆RAk,n≥0.Similarly if the novel scheme for the Drazin inverse is defined by left multiplying of Vn, we can state that(27)NVn⊇NAk,n≥0.Now, an application of definition of the Drazin inverse yields(28)AAD=ADA=PRAk,NAk.Proposition 4 along with (26), (27), and (28) could lead to(29)VnAAD=Vn=ADAVn,n≥0.To complete the proof, we proceed in what follows. The error matrix δn=AD-Vn satisfies (30)δn=AD-Vn=AD-ADAVn=ADI-AVn=ADRn.Using (23), we obtain(31)δn=ADRn≤ADR010n,which is a confirmation of (21). As a direct result of (31) and Proposition 4, we can obtain(32)Aδn+1=AAD-AVn+1=AAD-I+I-AVn+1=AAD-I+Rn+1.Considering (22) and applying several simplifications, we obtain that(33)Aδn+1=AAD-I+Rn10.Applying the idempotent property (I-AAD)t=(I-AAD),t≥1 being a positive integer from now on, and the following fact of (29):(34)I-AADAδn=I-AADAAD-Vn=Vn-AADVn=0,we obtain for each t≥1 that (here we use (34) in simplifications)(35)Rnt+AAD-I=I-AVnt+AAD-I=I-AAD+AAD-AVnt+AAD-I=I-AAD+Aδnt+AAD-I=I-AAD+Aδnt+AAD-I=Aδnt.From (35) and (33), we have(36)Aδn+1=Aδn10.Finally, we get that(37)Aδn+1≤Aδn10.It is now straightforward to calculate the error inequality of the proposed iteration (16) considering (37) and the second criterion of (4), when computing the Drazin inverse, as comes next(38)δn+1=Vn+1-AD=ADAVn+1-ADAAD=ADAVn+1-AAD≤ADAδn+1≤ADδn10.The inequalities in (38) lead to the fact that Vn→AD as n→+∞ with the tenth convergence speed. The proof is ended now.Remark 6.
An important issue is to find initial approximationsV0. In accordance with Proposition 3, V0 must read as the following relation to guarantee the convergence in the Drazin inverse case:(39)max1≤i≤t⁡1-λiAV0<1,where(40)rankAV0=h,and λi(AV0),i=1,2,…,h, are eigenvalues of AV0.


## 5. Efficiency Comparison

The computational complexity of the matrix inverse isO(mn-squared), where the order of the matrix is m×n (a rectangular matrix). Broadly speaking, if the complexity O(n-cubed) is brought down to O(nk), where k<3, say 2.9 for a general matrix, then it is definitely an achievement. Reducing the complexity by a linear factor is not attractive since what we have currently in 2018 is exa-flops computing speed and desktop and laptop computers are available in 100s of millions worldwide unlike the period during mid-20th centuries.Thus, we calculate and compare theoretically the computational efficiencies of various schemes (5), (8), (9), (10), (11), and (16), since they all could converge to AD upon the choice of a suitable initial value (6).By considering a unit cost for the arithmetic operations, typical of the floating point calculations, we can take into account theinverse-finder informational efficiency index. This index applies two values ρ and κ, which are for the convergence speed and the number of matrix-matrix products, respectively. Hence, this index is given by [20](41)IIEI=ρκ.Apart from this index, another useful one which is mainly in agreement with the CPU elapsed time of the iterative methods in numerical linear algebra is thecomputational efficiency index defined by [23](42)CEI=ρ1/κ.Therefore, a fruitful scheme in analytical standpoint should achieve the speed ρ with fewer matrix products κ in contrast to the competitors (viz, κ≤ρ).In Table1, we provide a comparison of the number of matrix products and the rate of convergence, accompanied by (41) and (42) for various schemes. The results indicate that the proposed scheme in Section 3 is more efficient than its competitors.Table 1
A comparison on the calculational burdensome of various methods applying different indices.
Methods (5) (8) (9) (10) (11) (16) ρ 2 3 3 4 9 6 κ 2 3 4 4 9 10 IIEI 2 2 = 1 3 3 = 1 3 4 = 0.75 4 4 = 1 9 9 = 1 10 6 ≈ 1.666 CEI 2 1 / 2 ≃ 1.414 3 1 / 3 ≃ 1.442 3 1 / 4 ≃ 1.316 4 1 / 4 ≃ 1.414 9 1 / 9 ≃ 1.276 1 0 1 / 6 ≃ 1.467As a matter of fact, one may observe that the proposed iteration (16) decreases the calculation burdensome via applying fewer operations and yields a better balance between the higher order and the computational load. This fact will be numerically investigated later in Section 6.


## 6. Experiments

This section studies issues associated with the computational accuracy and times of finding the Drazin inverse. The new iteration (16) is free from matrix power in its implementation and this allows one to apply it for finding generalized inverses easily. For computational comparisons, we have used the methods (5), (8), (9), (10), (11), and (16), denoted by “SM2," “CM3,” “LM3,” “KMS4,” “HM9,” and “PM10,” respectively.A couple of remarks are in order:(i)
The simulations are done in Mathematica 11.0 [24].(ii)
The time needed to obtain the approximate inverses is reported in seconds.(iii)
Whenever the elapsed times are reported, the compared methods are programmed in the same environment.Experiment 1.
The aim of this experiment is to calculate the Drazin inverse while the input matrix is [17](43)A=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,with k=ind(A)=3. The exact Drazin inverse is given by (44)AD=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.5156314.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.The stopping termination isVk+1-Vk∞≤ϵ with ϵ=10-8. Checking the conditions of Definition 2 for our proposed iterative method PM10 yields(45)Ak+1AD-Ak∞=3.69638×10-12,ADAAD-AD∞=8.43992×10-10,AAD-ADA∞=3.75205×10-10,which supports the theoretical discussions.The practical application of the new scheme (16) lies in several problems; see, for example, [25, 26]. For instance, in the process of resolving second-kind integral equations via Wavelet-like technique, the whole discretized problem will be yielded to calculate the inverse of a large sparse matrix [27].Additionally, we can apply/use (16), as an approach to provide good preconditioners for speeding-up modern Krylov methods, such as GMRES or BiCGSTAB, for solving large scale sparse linear systems; see, e.g., [28]. For such a task, we need to define a new initial matrix defined in Table 2 [29].Table 2
Several choices forV0 (for calculating A-1).
Different ways Form 1 Form 2 Form 3 Form 4 Form 5 Initial value A A 1 2 A A ∞ 2 A A F 2 A T N A 1 A ∞ A ∗ A 2 2In modern numerical linear algebra, schemes like (16) 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.It is necessary to test the behavior of the new method in a fair environment with a clear comparison taking into account various competitors. On the other hand, since the generation of random square matrices having Drazin inverses is difficult, in what follows, we compare various competitors in terms of the elapsed computational time so as to attain regular approximate inverses for large sparse matrices.Experiment 2.
The aim of this experiment is to compare the elapsed CPU times of different methods for the following 25 large sparse random complex matrices:SeedRandom[123]n = 5000; number = 25;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];
HereI=-1. In this test, the stopping criterion is Vn+1-Vn1≤10-6 and the maximum number of iterations allowed is set to 75. Moreover, the initial value is constructed via Form 3 of Table 2.The results of comparisons for this test are presented in Figure1. As is obvious that higher order methods require lower number of iterations to converge, thus we then put our focus on the computational time needed to satisfy the desired tolerance. As could be observed, our scheme (16) beats its competitors and also supports the analytical reports in Table 1.Figure 1
Comparison of the elapsed times for various methods in Experiment2.Experiment 3.
Here we rerun Experiment2 with the stopping termination Vn+1-Vn∞≤10-4 and initial value chosen by Form 2 of Table 2. The results are summarized in Figure 2. Furthermore, to check the usefulness of PM10, we rerun Experiment 2 with larger sizes matrices, i.e., when n=10000, the stopping termination Vn+1-Vn∞≤10-4 and initial value chosen by Form 1 of Table 2. The results of comparisons for this case are given in Figure 3.Figure 2
Comparison of the elapsed times for various methods in Experiment3.Figure 3
Comparison of the elapsed times for various methods in Experiment3 for larger sizes A10000×10000 and a different initial value.The attained results have reverified the robustness of the proposed iterative method (16) by a clear reduction in the elapsed CPU times.Note that the 10th-order method is better than a fourth-order method in CPU time. However, this is more tangible if a sharp initial approximation is chosen meaning that the 10th-order method arrives at the convergence phase quickly.We also emphasise that the construction of any higher order method is meaningful only if we observe an improvement in CEI as discussed in Section5. Accordingly, a member of the family (13) should be chosen such that we could arrive at the unknown parameters (e.g., like (17)) in order to improve CEI. An attempt to find an optimal iteration in this way is still a research topic in the field.Remark 7.
With the standard precision of 15 digits (most widely used globally, e.g., in Matlab or Mathematica), convergence of order 2 or 3 has been found to be computationally (amount of numerical computation) optimal [19].


## 7. Summary and Remarks

The Drazin inverse is investigated in the matrix theory (particularly in the topic of generalized inverses) and also in the ring theory; see, e.g., [30].In this work, we have investigated a higher-order matrix scheme (16) for calculating the Drazin inverse. Convergence analysis of our scheme has been discussed and an investigation on the selection of the initial approximation so as to initiate the iterates and keep the convergence speed was furnished. We also discussed under what conditions the novel scheme can be taken into account for calculating the Drazin inverse of square matrices.Furthermore, the total elapsed timings consuming of the proposed scheme (16) is low in contrast to the competitors of this type in the case of constructing approximate inverses.Tackling on the generalization of the new scheme (16) for interval matrix inversion or the application of such matrix methods in option pricing in order to act as a preconditioner to reduce the ill-conditioning of the large sized matrices occuring in the process of pricing (see, e.g., [31]) could be taken into account as future investigations in this field of study.



---

*Source: 6758302-2018-10-17.xml*