# A Higher Order Iterative Method for Computing the Drazin Inverse

Authors: F. Soleymani; Predrag S. Stanimirović
Journal: The Scientific World Journal
 (2013)

---

## 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.


---

## Body

## 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 matrixA, is the Schulz method given in [4] as follows:
(1)Vn+1=Vn(2I-AVn),n=0,1,2,…,
where I is the identity matrix with the same dimension that of the matrix A. 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 (V0) to start the process and converge to A-1. In practice, the Schulz-type methods are efficient (specially for structured matrices) but a difficulty lies in the initial approximation of the A-1. 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)V0=A∥A∥∞2.
In fact, by choosing (2), we attain
(3)∥V0∥∞∥A-1∥∞=∥A∥∞∥A∥∞2∥A-1∥∞.
Based on the fact that κ(A)=∥A∥∥A-1∥≥∥I∥≥1, we could have
(4)∥V0∥∞∥A-1∥∞=1κ∞(A)≤1,
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 approximationV0 have been listed in Table 1, where AT and A* are the transpose and the conjugate transpose of the complex matrix A, respectively, and N stands for the size of the square matrix.Table 1
Some of the general ways to chooseV0.
Forms
Form 1
Form 2
Form 3
Form 4
Form 5
Initial formulation
A
/
∥
A
∥
1
2
A
/
∥
A
∥
∞
2
A
/
∥
A
∥
F
2
A
^
T
/
(
N
∥
A
∥
1
∥
A
∥
∞
)
A
^
*
/
∥
A
∥
2
2A vast discussion on choosing the initial approximationV0 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 A, one can choose V0 as follows:
(5)V0=I∥A∥F.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:
(6)V0=diag(1a11,1a22,…,1aNN),
with aii as the diagonal elements of A.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
(7)Vn+1=Vn(3I-AVn(3I-AVn)),n=0,1,2,…
of third-order and also proposed another iterative method for finding A-1 of the same order as it is given in
(8)Vn+1=Vn[I+12(I-AVn)(I+(2I-AVn)2)],n=0,1,2,….Note that a general procedure for constructing such methods was given in ([11], Chapter 5). Krishnamurthy and Sen provided the following fourth-order method:
(9)Vn+1=Vn(I+Yn(I+Yn(I+Yn))),n=0,1,2,…,
in which Yn=I-AVn. As another example, a ninth-order method could be presented as (10)Vn+1=Vn(I+Yn(I+Yn(I+Yn(I+Yn(I+Yn(I+Yn(I+Yn(I+Yn)))))))),n=0,1,2,….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 Section3, 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 equationF(V)=V-1-A. Then, one may obtain an iterative process using proper factorization. In this way, by applying the iterative scheme
(11)Yn=Vn-F(Vn)F′(Vn),Zn=Vn-F(Vn)2(1F′(Vn)+1F′(Yn)),Vn+1=Zn-F(Zn)2F′(Zn)(2+F′′(Zn)F(Zn)F′(Zn)2),n=0,1,2,…
on the matrix equation V-1-A=0, we are able to find the fixed-point iteration(12)Vn+1=-18Vn(-7I+9AVn-5(AVn)2+(AVn)3)×(12I-42AVn+103(AVn)2-156(AVn)3+157(AVn)4-104(AVn)5+43(AVn)6-10(AVn)7+(AVn)8).
Now, by simplification, we suggest the following matrix iteration:
(13)χn=-7I+ψn(9I+ψn(-5I+ψn)),ϑn=ψnχn,Vn+1=-18Vnχn(12I+ϑn(6I+ϑn)),n=0,1,2,…,
wherein ψn=AVn,  I is the identity matrix, and the sequence of iterates {Vn}n=0n=∞ converges to A-1 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 Theorem1.Theorem 1.
LetA=[aij]∈ℂN×N be a nonsingular complex matrix. If the initial approximation V0 satisfies
(14)∥I-AV0∥<1,
then the iterative method (13) converges with ninth order to A-1.Proof.
We use notations thatE0=I-AV0 and subsequently En=I-AVn. Then, (15)En+1=I-AVn+1=I-A(-18Vn(-7I+AVn(9I+AVn(-5I+AVn)))×(12I+AVn(-7I+AVn(9I+AVn(-5I+AVn)))×(6I+AVn(-7I+AVn(9I+AVn(-5I+AVn)))))18)=I-A(-18Vn(-7I+9AVn-5(AVn)2+(AVn)3)×(12I-42(AVn)1+103(AVn)2-156(AVn)3+157(AVn)4-104(AVn)5+43(AVn)6-10(AVn)7+(AVn)8)18)=18(-2I+AVn)3(-I+AVn)9=18(I+I-AVn)3(I-AVn)9=18(I+En)3En9=18(En9+3En10+3En11+En12).Hence, by taking arbitrary matrix norm on both sides of (15), we attain
(16)∥En+1∥≤18(∥En∥9+3∥En∥10+3∥En∥11+∥En∥12).
In addition, since ∥E0∥<1, by relation (16), we obtain that
(17)∥E1∥≤18(∥E0∥9+3∥E0∥10+3∥E0∥11+∥E0∥12)≤∥E0∥9<1.
Now if we consider ∥En∥<1, therefore
(18)∥En+1∥≤18(∥En∥9+3∥En∥10+3∥En∥11+∥En∥12)≤∥En∥9.
Using mathematical induction, we obtain
(19)∥En+1∥≤∥En∥9,n≥0.
Furthermore, we get that
(20)∥En+1∥≤∥En∥9≤⋯≤∥E0∥9n+1.
That is, I-AVn→0, when n→∞ and
(21)Vn⟶A-1,asn⟶∞.
Thus, the new method (13) converges to the inverse of the matrix A in the case ρ(AV0)<1, where ρ is the spectral radius. Now, we prove that the order of convergence for the sequence {Vn}n=0n=∞ is at least nine. Let εn denotes the error matrix εn=A-1-Vn; afterwards
(22)Aεn=I-AVn=En.
The identity (22) in conjunction with (15) implies that
(23)Aεn+1=18((Aεn)9+3(Aεn)10+3(Aεn)11+(Aεn)12).
Therefore, using invertibility of A, it follows immediately that
(24)εn+1=18(εn(Aεn)8+3εn(Aεn)9+3εn(Aεn)10+εn(Aεn)11).
By taking any subordinate norm of (24), we obtain
(25)∥εn+1∥≤(18(∥A∥8+3∥A∥9∥εn∥+3∥A∥10∥εn∥2+∥A∥11∥εn∥3)18)∥εn∥9.
Consequently, it is proved that the iterative formula (13) converges to A-1, 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 A 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{Vn}n=0n=∞ produced from (13), under a certain condition (when AV0=V0A), may be applied not only to the left preconditioned linear system VnAx=Vnb but also to the right preconditioned linear system AVny=b, where y=Vnx. 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 A-1 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
(26)IIEI=ρκ.Hence, a favorable method in theoretical point of view must reach an orderρ with fewer matrix multiplications κ, (i.e., κ≤ρ).In Table2, 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.
Methods
(1)
(7)
(8)
(9)
(10)
(13)
Rate of convergence
2
3
3
4
9
9
Number of matrix multiplications
2
3
4
4
9
7
IIEI
2
/
2
=
1
3
/
3
=
1
3
/
4
=
0.75
4
/
4
=
1
9
/
9
=
1
9
/
7
≈
1.285


## 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 integerk, such that rank(Ak+1)=rank(Ak), is called the index of A and denoted by ind(A).Definition 3.
LetA be an N×N complex matrix; the Drazin inverse of A, denoted by AD, is the unique matrix V satisfying
(27)(1k)AkVA=Ak,(2)VAV=V,(5)AV=VA,
where k=ind(A) is the index of A.Note that ifind(A)=1, the matrix V is called the group inverse of A. Also, if A is nonsingular, then it is easily seen that ind(A)=0 and AD=A-1. Note that the idempotent matrix AAD is the projector on ℛ(Ak) along 𝒩(Ak), where ℛ(Ak) denotes the range of Ak and 𝒩(Ak) is the null space of Ak.Wei in [15] proved that the general solution of the square singular linear system Ax=b can be obtained using the Drazin inverse as x=ADb+(I-AAD)z, where z∈ℛ(Ak-1)+𝒩(A).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:
(28)V0=W0=αAl,l≥ind(A)=k,
where the parameter α must be chosen so that the condition ∥I-AV0∥<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
(29)V0=W0=2Tr(Ak+1)Ak,
wherein Tr(·) stands for the trace of an arbitrary square matrix.In what follows, we use the following auxiliary results.Proposition 4 (see [17]).
LetM∈ℂN×N and ε>0 be given. There is at least one matrix norm ∥·∥ such that
(30)ρ(M)≤∥M∥≤ρ(M)+ϵ,
where ρ(M) denotes the spectral radius of M.Proposition 5 (see [18]).
IfPL,M denotes the projector on a space L along a space M, then(i)
PL,MQ=Q if and only if ℛ(Q)⊆L;(ii)
QPL,M=Q if and only if 𝒩(Q)⊇M.Theorem 6.
LetA∈ℂN×N be singular square matrix. Also, suppose that the initial approximation W0 is chosen by means of (28). Then the sequence {Wn}n=0n=∞ defined by the iterative method (13) satisfies the following error estimate when finding the Drazin inverse:
(31)∥AD-Wn∥∥AD∥≤∥I-AW0∥9n.Proof.
Consider the notationF0=I-AW0 and subsequently the residual matrix as Fn=I-AWn. Then similarly as in (15), we get(32)Fn+1=I-AWn+1=I-A(-18Wn(-7I+AWn(9I+AWn(-5I+AWn)))×(12I+AWn(-7I+AWn(9I+AWn(-5I+AWn)))×(6I+AWn(-7I+AWn(9I+AWn(-5I+AWn)))))18)=18(I+I-AWn)3(I-AWn)9=18(Fn9+3Fn10+3Fn11+Fn12).By taking an arbitrary matrix norm on both sides of (32), we attain
(33)∥Fn+1∥≤18(∥Fn∥9+3∥Fn∥10+3∥Fn∥11+∥Fn∥12).
In addition, since ∥F0∥<1, by relation (33), we obtain that ∥F1∥≤(1/8)(∥F0∥9+3∥F0∥10+3∥F0∥11+∥F0∥12)≤∥F0∥9. Similarly, ∥Fn+1∥≤(1/8)(∥Fn∥9+3∥Fn∥10+3∥Fn∥11+∥Fn∥12)≤∥Fn∥9. Using mathematical induction, we obtain ∥Fn+1∥≤∥Fn∥9 for each n≥0. This implies that
(34)∥Fn∥9≤∥F0∥9n,n≥0.
SinceW0 is chosen as in (28), ℛ(W0)⊆ℛ(Ak) immediately follows. Further, using this fact in conjunction with (12), which implies that ℛ(Wn)⊆ℛ(Wn-1), we conclude that
(35)ℛ(Wn)⊆ℛ(Ak),n≥0.
Similarly, if we rewrite (13) in the form
(36)Wn+1=-18(-7I+9WnA-5(WnA)2+(WnA)3)×(12I-42AWn+103(WnA)2-156(WnA)3+157(WnA)4-104(WnA)5+43(WnA)6-10(WnA)7+(WnA)8)Wn.
it is not difficult to verify that
(37)𝒩(Wn)⊇𝒩(Ak),n≥0.
Now, an application of the well-known results from [12]
(38)AAD=ADA=Pℛ(Ak),𝒩(Ak)
in conjunction with Proposition 5 and (35), (37), immediately follows
(39)WnAAD=Wn=ADAWn,n≥0.
Therefore, the error matrix δn=AD-Wn satisfies
(40)δn=AD-Wn=AD-ADAWn=AD(I-AWn)=ADFn.
From the last identity and (34) we have
(41)∥δn∥=∥AD∥∥Fn∥≤∥AD∥∥F0∥9n,
which is a confirmation of (31).The following result is a consequence of Theorem6.Corollary 7.
If the conditions of Theorem6 are satisfied and the initial iteration W0 is chosen such that
(42)∥F0∥=∥I-AW0∥<1,
the iterative method (13) converges to AD.Therefore, our goal is to find initial approximationsW0 satisfying (42). In accordance with Proposition 4, W0 must satisfy the following inequality to ensure the convergence in the Drazin inverse case:
(43)max1≤i≤t|1-λi(AW0)|<1,
where rank(AW0)=t and λi(AW0),i=1,…,t are eigenvalues of AW0.Theorem 8.
LetA∈ℂN×N be a singular square matrix with ind(A)=k, and the sequences Fn, δn are defined as in Theorem 6. The sequence {Wn}n=0n=∞ defined by the iterative method (13) converges with ninth order to AD if the initial approximation W0 is in accordance with (42).Proof.
Now, by consideringδn=AD-Wn as the error matrix for finding the Drazin inverse, we have
(44)Aδn+1=AAD-AWn+1=AAD-I+I-AWn+1=AAD-I+Fn+1.
Taking into account (32) and using elementary algebraic transformations, we further derive
(45)Aδn+1=AAD-I+18(Fn9+3Fn10+3Fn11+Fn12)=18((Fn9+3Fn10+3Fn11+Fn12)+8(AAD-I))=18((Fn9+AAD-I)+3(Fn10+AAD-I)+3(Fn11+AAD-I)+(Fn12+AAD-I)).
Now, using the idempotent property (I-AAD)t=(I-AAD), t≥1, and the following consequence of (39)
(46)(I-AAD)Aδn=(I-AAD)A(AD-Wn)=Wn-AADWn=0,
we obtain for each t≥1 the following:
(47)(Fn)t+AAD-I=(I-AWn)t+AAD-I=(I-AAD+AAD-AWn)t+AAD-I=((I-AAD)+Aδn)t+AAD-I=I-AAD+(Aδn)t+AAD-I=(Aδn)t.
From (47) and (45),
(48)Aδn+1=18((Aδn)9+3(Aδn)10+3(Aδn)11+(Aδn)12)
Therefore,
(49)∥Aδn+1∥≤18(∥Aδn∥9+3∥Aδn∥10+3∥Aδn∥11+∥Aδn∥12)≤∥Aδn∥9.
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:
(50)∥δn+1∥=∥Wn+1-AD∥=∥ADAWn+1-ADAAD∥=∥AD(AWn+1-AAD)∥≤∥AD∥∥Aδn+1∥≤∥AD∥∥A∥9∥δn∥9.
Therefore, since (42) is satisfied, from (31) follows δn→0. Furthermore, the inequalities in (50) immediately lead to the conclusion that Wn→AD as n→+∞ 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 Vi.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 Algorithm1.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}];Note thatI=-1. In this test, the stopping criterion is ∥Vn+1-Vn∥1≤10-6 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 Figures1, 2, 3, and 4. In Figure 1, we show the number of iterations required by different methods to attain the desired accuracy when V0=A/∥A∥12. 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 V0=A/∥A∥12, V0=A/∥A∥∞2, and V0=A/∥A∥F2. 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 usingV0=A/∥A∥12.Figure 2
Comparison of the elapsed time for solving the Test Problem 1 usingV0=A/∥A∥12.Figure 3
Comparison of the elapsed time for solving the Test Problem 1 usingV0=A/∥A∥∞2.Figure 4
Comparison of the elapsed time for solving the Test Problem 1 usingV0=A/∥A∥F2.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:
(51)y′′=3y-2y′,t∈[a,b],y(a)=e3,y(b)=e-3,
where the domain of the solution is a=0 and b=2. To solve (51) using 3-point 𝔽𝔻 discretization, we assume that in the grid points ti, the exact solution is denoted by y(ti) and the approximate solution is defined by wi, for any i=0,1,…,n, n+1.Thus, (51) can be written in discretized form in Mathematica language (see Algorithm 2), wherein n is the number of grid points which result in an n×n sparse linear system of equations. In this test, we consider the tolerance as 10-8 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 (VAx=Vb) 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 n=1500.
Methods
GMRES
PGMRES-(1)-VI
PGMRES-(7)-IV
PGMRES-(8)-III
PGMRES-(9)-III
PGMRES-(13)-II
Total time
Fail
3.50
1.78
1.53
2.62
1.31Table 4
Comparison of the computational time in solving the linear system resulting of discretization of (51) when n=2000.
Methods
BiCGSTAB
PBiCGSTAB-(1)-III
PBiCGSTAB-(7)-III
PBiCGSTAB-(8)-II
PBiCGSTAB-(9)-II
PBiCGSTAB-(13)-I
Total time
1.98
0.53
0.60
0.53
0.54
0.46Algorithm 2
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}];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 V0 for the compared methods using (6).Test Problem 3.
The aim of this example is to apply the discussions of Section4, for finding the Drazin inverse of the following square matrix (taken from [16]):
(52)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. 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).Algorithm 3
P
M
[
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&)]]
];The two-argument functionDrazinInverse[A_,tolerance_] takes the arbitrary matrix A and the tolerance from the user to obtain its Drazin inverse by knowing the index k. In this case, by choosing tolerance=10-8, we obtain(53)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].Checking the conditions of Definition3 yields to
(54)∥Ak+1AD-Ak∥∞=1.48415×10-12,∥ADAAD-AD∥∞=1.20264×10-10,∥AAD-ADA∥∞=8.93836×10-11,
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.



---

*Source: 708647-2013-09-30.xml*