# A New High-Order Stable Numerical Method for Matrix Inversion

Authors: F. Khaksar Haghani; F. Soleymani
Journal: The Scientific World Journal
 (2014)

---

## Abstract

A stable numerical method is proposed for matrix inversion. The new method is accompanied by theoretical proof to illustrate twelfth-order convergence. A discussion of how to achieve the convergence using an appropriate initial value is presented. The application of the new scheme for finding Moore-Penrose inverse will also be pointed out analytically. The efficiency of the contributed iterative method is clarified on solving some numerical examples.


---

## Body

## 1. Introduction and Preliminary Notes

It is well known that the inverse of a square matrixA
m
×
m, which is also known as a reciprocal matrix, is a matrix A
-
1 such that A
A
-
1
=
I, where I is the identity matrix. A regular nonsingular matrix A
m
×
m can be inverted using methods such as the Gauss-Jordan elimination or Gaussian elimination method. Such schemes fall in the category of direct methods for this purpose.The direct methods cannot properly handle sparse matrices possessing sparse inverses arising in the numerical solution of integral equations [1]. On the other hand, methods such as conjugate gradient for symmetric positive definite matrices and GMRES are effective for large sparse linear systems. However, there is a problem when the coefficient matrix (when solving a linear system of equations) is ill-conditioned. To remedy this, one may apply a preconditioner to the system in which its construction is not an easy task [2].An iterative method for preconditioning is the SPAI (sparse approximate inverse preconditioner) algorithm [3]. Given a sparse matrix A the SPAI algorithm computes a sparse approximate inverse M by minimizing ∥
A
M
-
I
∥ in the Frobenius norm. Then, the approximate inverse is computed explicitly and can be applied as a preconditioner to an iterative method.There are other types of schemes, which can be considered as iteration methods while they have different structures; see, for example, [4, 5]. In such iterative methods, at each iteration an approximate inverse of a matrix (if it is rectangular, one can find Moore-Penrose inverse) may easily be attained. And consequently, the users have the ability to solve the linear systems (with multiple right-hand side vectors) iteratively or use the approximate inverses in sensitivity analysis and the preconditioning of a linear system. This type of methods is in focus here.Several known methods were proposed for approximating matrix inverse, such as those based on the so-called minimum residual iterations and Hotelling-Bodewig algorithm [6]. The Hotelling-Bodewig algorithm is defined by
(1)
V
n
+
1
=
V
n
(
2
I
-
A
V
n
)
,
n
=
0,1
,
2
,
…
,
where I is the identity matrix. Schulz in [7] found that the eigenvalues of I
-
A
V
0 must have magnitudes less than 1 to ensure the convergence, which is a key element in designing higher efficient Schulz-type iterative methods.In 2011, Li et al. in [8] theoretically investigated
(2)
V
n
+
1
=
V
n
(
3
I
-
A
V
n
(
3
I
-
A
V
n
)
)
,
n
=
0,1
,
2
,
…
,
and also proposed another third-order iterative method for finding A
-
1 as follows:
(3)
V
n
+
1
=
[
I
+
1
4
(
I
-
V
n
A
)
(
3
I
-
V
n
A
)
2
]
V
n
,
n
=
0,1
,
2
,
…
.The iterative method (2) can also be found in Chapter 5 of the textbook [9]. As an another example from this primary source, the authors provided the following twelfth-order method:(4)
V
n
+
1
=
V
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
(
I
+
Y
n
)
)
)
)
)
)
)
)
)
)
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
n
=
0,1
,
2
,
…
,in which Y
n
=
I
-
A
V
n. For further reading, one may refer to [10–12].In this paper, we will propose an efficient iterative method for findingA
-
1 numerically. The theoretical convergence of the method will also be studied. We also discuss the application of the new scheme in finding Moore-Penrose inverse (also known as pseudoinverse) for rectangular or singular matrices. It is also proven analytically that the new method has asymptotical stability in general. Some large-scale sparse matrices will be taken into account as some examples to put on show a clear reduction of the execution time when the new algorithm is applied.The rest of the paper is organized as follows. The main contribution of this paper is given in Sections2-3. Section 2 is devoted to the analysis of convergence which shows that the method can be considered for the pseudoinverse as well. Section 3 thoroughly and for the first time studies the stability of this high-order Schulz-type iterative method for finding generalized inverses. Section 4 covers the matter of initial guess/value in order to preserve the convergence order. Subsequently, the method is examined in Section 5 numerically. And finally, concluding remarks are presented in Section 6.


## 2. A New Method and Its Convergence Study

By applying the following nonlinear equation solver (to see the new developments on root-finding methods, refer to [13])
(5)
y
n
=
x
n
-
2
-
1
f
′
(
x
n
)
-
1
f
(
x
n
)
,
z
n
=
x
n
-
f
(
y
n
)
-
1
f
(
x
n
)
,
u
n
=
z
n
-
[
[
z
n
-
x
n
]
-
1
(
f
(
z
n
)
-
f
(
x
n
)
)
]
-
1
f
(
z
n
)
,
g
n
=
u
n
-
f
(
u
n
)
-
1
f
(
u
n
)
,
x
n
+
1
=
g
n
-
[
[
g
n
-
u
n
]
-
1
(
f
(
g
n
)
-
f
(
u
n
)
)
]
-
1
f
(
g
n
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
n
=
0,1
,
2
,
…
,
on the nonlinear matrix equation A
V
=
I, we obtain a fixed-point iteration for matrix inversion using ψ
n
=
A
V
n as follows:
(6)
V
n
+
1
=
1
64
V
n
ϕ
(
V
n
)
=
1
64
V
n
(
816
I
-
4812
ψ
n
+
17393
ψ
n
2
-
43044
ψ
n
3
h
h
h
h
h
h
+
77154
ψ
n
4
-
103312
ψ
n
5
+
105039
ψ
n
6
h
h
h
h
h
h
-
81576
ψ
n
7
+
48268
ψ
n
8
-
21516
ψ
n
9
h
h
h
h
h
h
+
7071
ψ
n
10
-
1652
ψ
n
11
+
258
ψ
n
12
h
h
h
h
h
h
-
24
ψ
n
13
+
ψ
n
14
)
.
Simplifying (6) by proper factorizing yields
(7)
ζ
n
=
17
I
+
ψ
n
(
-
28
I
+
ψ
n
(
22
I
+
ψ
n
(
-
8
I
+
ψ
n
)
)
)
,
κ
n
=
ψ
n
ζ
n
,
V
n
+
1
=
1
64
V
n
ζ
n
(
48
I
+
κ
n
(
-
12
I
+
κ
n
)
)
,
n
=
0,1
,
2
,
…
,
wherein the sequence of iterates {
V
n
}
n
=
0
n
=
∞ converges to A
-
1 for a good initial value. Such a guess, V
0, will be discussed in Section 4.Theorem 1.
LetA
=
[
a
i
j
]
m
×
m be a nonsingular real or complex matrix. If the initial approximation V
0 satisfies
(8)
∥
I
-
A
V
0
∥
<
1
,
then the iterative method (7) converges with at least twelfth order to A
-
1.Proof.
Let∥
I
-
A
V
0
∥
<
1, and for the sake of simplicity assume that E
0
=
I
-
A
V
0 and E
n
=
I
-
A
V
n
=
I
-
ψ
n stand for the symmetric residual matrix. It is straightforward to have
(9)
E
n
+
1
=
I
-
A
V
n
+
1
=
I
-
A
(
1
64
V
n
(
816
I
-
4812
ψ
n
+
17393
ψ
n
2
-
43044
ψ
n
3
h
h
h
h
h
h
h
h
h
h
h
h
+
77154
ψ
n
4
-
103312
ψ
n
5
+
105039
ψ
n
6
h
h
h
h
h
h
h
h
h
h
h
h
-
81576
ψ
n
7
+
48268
ψ
n
8
-
21516
ψ
n
9
h
h
h
h
h
h
h
h
h
h
h
h
+
7071
ψ
n
10
-
1652
ψ
n
11
h
h
h
h
h
h
h
h
h
h
h
h
+
258
ψ
n
12
-
24
ψ
n
13
+
ψ
n
14
17393
ψ
n
2
)
1
64
)
=
-
1
64
(
-
4
I
+
ψ
n
)
3
(
-
I
+
ψ
n
)
12
=
1
64
(
3
I
+
I
-
ψ
n
)
3
(
I
-
ψ
n
)
12
=
1
64
(
3
I
+
E
n
)
3
E
n
12
=
1
64
(
27
E
n
12
+
27
E
n
13
+
9
E
n
14
+
E
n
15
)
.
Hence, we attain
(10)
∥
E
n
+
1
∥
≤
1
64
(
27
∥
E
n
∥
12
+
27
∥
E
n
∥
13
+
9
∥
E
n
∥
14
+
∥
E
n
∥
15
)
.
In addition, since ∥
E
0
∥
<
1, by relation (10) and using mathematical induction, we obtain that ∥
E
1
∥
≤
(
1
/
64
)
(
27
∥
E
0
∥
12
+
27
∥
E
0
∥
13
+
9
∥
E
0
∥
14
+
∥
E
0
∥
15
)
≤
∥
E
0
∥
12
<
1. If we consider ∥
E
n
∥
<
1, therefore
(11)
∥
E
n
+
1
∥
≤
1
64
(
27
∥
E
n
∥
12
+
27
∥
E
n
∥
13
+
9
∥
E
n
∥
14
+
∥
E
n
∥
15
)
≤
∥
E
n
∥
12
.
Furthermore, we get that
(12)
∥
E
n
+
1
∥
≤
∥
E
n
∥
12
≤
⋯
≤
∥
E
0
∥
1
2
n
+
1
<
1
.
That is, I
-
A
V
n
→
0, when n
→
∞, and thus V
n
→
A
-
1, as n
→
∞.
Now, we must show the twelfth order using the sequence{
V
n
}
n
=
0
n
=
∞. To do this, we denote ε
n
=
V
n
-
A
-
1 as the error matrix in the iterative procedure (7). We have
(13)
I
-
A
V
n
+
1
=
1
64
[
27
(
I
-
A
V
n
)
12
+
27
(
I
-
A
V
n
)
13
+
9
(
I
-
A
V
n
)
14
+
(
I
-
A
V
n
)
15
]
.
Hence, we could get that
(14)
A
(
A
-
1
-
V
n
+
1
)
=
1
64
[
27
A
12
(
A
-
1
-
V
n
)
12
+
27
A
13
(
A
-
1
-
V
n
)
13
+
9
A
14
(
A
-
1
-
V
n
)
14
+
A
15
(
A
-
1
-
V
n
)
15
]
.
By multiplying A
-
1 by the left side, we have(15)
A
-
1
-
V
n
+
1
=
1
64
[
27
A
11
(
A
-
1
-
V
n
)
12
+
27
A
12
(
A
-
1
-
V
n
)
13
+
9
A
13
(
A
-
1
-
V
n
)
14
+
A
14
(
A
-
1
-
V
n
)
15
]
,
which now by taking an arbitrary matrix norm results in
(16)
∥
ε
n
+
1
∥
≤
1
64
[
27
∥
A
∥
11
∥
ε
n
∥
12
+
27
∥
A
∥
12
∥
ε
n
∥
13
+
9
∥
A
∥
13
∥
ε
n
∥
14
+
∥
A
∥
14
∥
ε
n
∥
15
]
.
And thus
(17)
∥
ε
n
+
1
∥
=
(
1
64
[
27
∥
A
∥
11
+
27
∥
A
∥
12
∥
ε
n
∥
1
+
9
∥
A
∥
13
∥
ε
n
∥
2
+
∥
A
∥
14
∥
ε
n
∥
3
]
1
64
[
27
∥
A
∥
11
+
27
∥
A
∥
12
∥
ε
n
∥
1
)
∥
ε
n
∥
12
.
That is, the iteration (7) converges with at least twelfth order to A
-
1. This concludes the proof.At this time, we discuss an application of (7) for finding the generalized inverses. The Moore-Penrose inverse of a complex matrix A
∈
ℂ
m
×
k (also called pseudoinverse), denoted by A
†
∈
ℂ
k
×
m, is a unique matrix V
∈
ℂ
k
×
m satisfying the following four Penrose equations:
(18)
A
V
A
=
A
,
V
A
V
=
V
,
(
A
V
)
*
=
A
V
,
(
V
A
)
*
=
V
A
,
wherein A
* is the conjugate transpose of A. Ben-Israel and his colleagues in [14, 15] used the method (1) with the starting value
(19)
V
0
:
=
α
A
*
,
where 0
<
α
<
2
/
ρ
(
A
*
A
) and ρ
(
·
) denotes the spectral radius. The authors in [15] further investigated that the sequence generated by
(20)
V
n
=
α
∑
i
=
0
n
‍
A
*
(
I
-
α
A
A
*
)
i
,
n
=
0,1
,
2
,
…
,
converges to the pseudoinverse.In the following theorem, we show analytically that in case of having singular or rectangular matrices, scheme (7) by considering the initial approximation (19) converges to the Moore-Penrose generalized inverse.Theorem 2.
For the sequence{
V
n
}
n
=
0
n
=
∞ generated by the iterative Schulz-type method (7), and any n
≥
0, it holds that
(21)
(
A
V
n
)
*
=
A
V
n
,
(
V
n
A
)
*
=
V
n
A
,
V
n
A
A
†
=
V
n
,
A
†
A
V
n
=
V
n
.Proof.
We will prove the conclusion by induction onn. For n
=
0, and considering (19), the first two equations of (21) can be demonstrated simply. And thus, we only give a verification to the last two equations as follows:
(22)
V
0
A
A
†
=
α
A
*
A
A
†
=
α
A
*
(
A
A
†
)
*
=
α
A
*
(
A
†
)
*
A
*
=
α
(
A
A
†
A
)
*
=
α
A
*
=
V
0
,
A
†
A
V
0
=
(
A
†
A
)
α
A
*
=
α
(
A
†
A
)
*
A
*
=
α
A
*
(
A
†
)
*
A
*
=
α
(
A
A
†
A
)
*
=
α
A
*
=
V
0
.
Assume now that the conclusion holds for some n
>
0. We now show that it continues to hold for n
+
1. Using the iterative method (7), one has
(23)
(
A
V
n
+
1
)
*
=
[
A
(
1
64
V
n
ϕ
(
V
n
)
)
]
*
=
1
64
(
816
ψ
n
*
-
4812
(
ψ
n
*
)
2
+
17393
(
ψ
n
*
)
3
h
h
h
h
h
-
43044
(
ψ
n
*
)
4
+
77154
(
ψ
n
*
)
5
-
103312
(
ψ
n
*
)
6
h
h
h
h
h
+
105039
(
ψ
n
*
)
7
-
81576
(
ψ
n
*
)
8
+
48268
(
ψ
n
*
)
9
h
h
h
h
h
-
21516
(
ψ
n
*
)
10
+
7071
(
ψ
n
*
)
11
-
1652
(
ψ
n
*
)
12
h
h
h
h
h
+
258
(
ψ
n
*
)
13
-
24
(
ψ
n
*
)
14
+
(
ψ
n
*
)
15
)
=
1
64
[
816
ψ
n
-
4812
ψ
n
2
+
17393
ψ
n
3
-
43044
ψ
n
4
h
h
h
h
h
+
77154
ψ
n
5
-
103312
ψ
n
6
+
105039
ψ
n
7
h
h
h
h
h
-
81576
ψ
n
8
+
48268
ψ
n
9
-
21516
ψ
n
10
h
h
h
h
h
+
7071
ψ
n
11
-
1652
ψ
n
12
+
258
ψ
n
13
-
24
ψ
n
14
+
ψ
n
15
]
=
A
(
1
64
V
n
ϕ
(
V
n
)
)
=
A
V
n
+
1
,
where the following fact (
A
V
n
)
*
=
A
V
n has been used. Thus, the first equality in (21) holds for n
+
1, and the second equality can be proved in a similar way. For the third equality in (21), using the assumption that V
n
A
A
†
=
V
n and the iterative method (7), we could write down
(24)
V
n
+
1
A
A
†
=
[
1
64
V
n
ϕ
(
V
n
)
]
A
A
†
=
1
64
(
816
V
n
A
A
†
-
4812
V
n
ψ
n
A
A
†
+
17393
V
n
ψ
n
2
A
A
†
h
h
h
h
h
-
43044
V
n
ψ
n
3
A
A
†
+
77154
V
n
ψ
n
4
A
A
†
h
h
h
h
h
-
103312
V
n
ψ
n
5
A
A
†
+
105039
V
n
ψ
n
6
A
A
†
h
h
h
h
h
-
81576
V
n
ψ
n
7
A
A
†
+
48268
V
n
ψ
n
8
A
A
†
h
h
h
h
h
-
21516
V
n
ψ
n
9
A
A
†
+
7071
V
n
ψ
n
10
A
A
†
h
h
h
h
h
-
1652
V
n
ψ
n
11
A
A
†
+
258
V
n
ψ
n
12
A
A
†
h
h
h
h
h
-
24
V
n
ψ
n
13
A
A
†
+
V
n
ψ
n
14
A
A
†
)
=
1
64
(
816
V
n
-
4812
V
n
ψ
n
+
17393
V
n
ψ
n
2
h
h
h
h
h
-
43044
V
n
ψ
n
3
+
77154
V
n
ψ
n
4
-
103312
V
n
ψ
n
5
h
h
h
h
h
+
105039
V
n
ψ
n
6
-
81576
V
n
ψ
n
7
+
48268
V
n
ψ
n
8
h
h
h
h
h
-
21516
V
n
ψ
n
9
+
7071
V
n
ψ
n
10
-
1652
V
n
ψ
n
11
h
h
h
h
h
+
258
V
n
ψ
n
12
-
24
V
n
ψ
n
13
+
V
n
ψ
n
14
)
=
1
64
V
n
ϕ
(
V
n
)
=
V
n
+
1
.
Hence, the third equality in (21) holds for n
+
1. The fourth equality can similarly be proved, and the desired result follows.The iterative method (7) is a matrix multiplication rich scheme. So, in order to reduce the computational load of matrix multiplications, it is enough to use SparseArray[mat] to avoid unnecessary multiplications to the nonzero elements when dealing with large sparse matrices.Remark 3.
The (inverse finder) informational efficiency index is defined asp
/
η where p and η stand for the local order and the number of matrix-matrix products per cycle. The proposed method of this paper requires 8 matrix-matrix multiplications to achieve the convergence order 12. This implies that 12
/
8
=
1.5 as its informational index, which is much better than 2
/
2
=
1, 3
/
3
=
1, 3
/
4
=
0.75, and 12
/
12
=
1 of schemes (1), (2), (3), and (4), respectively.Before stating the main theorem for finding Moore-Penrose inverse, it is required to recall that forA
∈
ℂ
m
×
k with the singular values σ
1
>
σ
2
>
⋯
σ
r
>
0 and the initial approximation V
0
=
α
A
* with 0
<
α
<
2
/
σ
1
2, it holds that
(25)
∥
A
(
V
0
-
A
†
)
∥
<
1
.
We are about to use this fact in the following theorem so as to find the theoretical order of the proposed method (7) for finding the Moore-Penrose inverse (see [16] for more details).Theorem 4.
ForA
∈
ℂ
m
×
k, with the singular values σ
1
>
σ
2
>
⋯
σ
r
>
0, the sequence {
V
n
}
n
=
0
n
=
∞ generated by (7) and using the initial approximation V
0
=
α
A
* converges to the Moore-Penrose inverse A
† in twelfth order provided that 0
<
α
<
2
/
σ
1
2.Proof.
Set𝔼
n
=
V
n
-
A
† and E
n
=
I
-
A
V
n; we have
(26)
A
𝔼
n
+
1
=
A
V
n
+
1
-
A
A
†
=
A
V
n
+
1
-
I
+
I
-
A
A
†
=
-
E
n
+
1
+
I
-
A
A
†
=
-
1
16
[
27
E
n
12
+
27
E
n
13
+
9
E
n
14
+
E
n
15
]
+
I
-
A
A
†
.
On the other hand, from the definitions of the Moore-Penrose inverse A
†, we have
(27)
(
I
-
A
A
†
)
k
=
I
-
A
A
†
,
k
=
1,2
,
…
;
(
I
-
A
A
†
)
A
𝔼
n
=
0
.
The use of these relations implies that
(28)
A
𝔼
n
+
1
=
1
16
[
-
27
(
A
𝔼
n
)
12
+
27
(
A
𝔼
n
)
13
h
h
h
h
-
9
(
A
𝔼
n
)
14
+
(
A
𝔼
n
)
15
]
.
So, for any matrix norm ∥
·
∥, we obtain
(29)
∥
A
𝔼
n
+
1
∥
≤
1
16
[
27
∥
A
𝔼
n
∥
12
+
27
∥
A
𝔼
n
∥
13
h
h
h
h
+
9
∥
A
𝔼
n
∥
14
+
∥
A
𝔼
n
∥
15
]
.
Applying (25), which implies that ∥
A
𝔼
0
∥
<
1, and a similar reasoning as in (10)–(12), one can obtain
(30)
∥
A
𝔼
n
+
1
∥
≤
1
16
[
27
∥
A
𝔼
n
∥
12
+
27
∥
A
𝔼
n
∥
13
+
9
∥
A
𝔼
n
∥
14
+
∥
A
𝔼
n
∥
15
]
≤
∥
A
𝔼
n
∥
12
≤
∥
A
∥
12
∥
𝔼
n
∥
12
.
Finally, using the properties of the Moore-Penrose inverseA
† and Theorem 2, it would be now easy to find the error inequality of the new scheme (7) as follows:
(31)
∥
V
n
+
1
-
A
†
∥
=
∥
A
†
A
V
n
+
1
-
A
†
A
A
†
∥
≤
∥
A
†
∥
∥
A
V
n
+
1
-
A
A
†
∥
=
∥
A
†
∥
∥
A
𝔼
n
+
1
∥
≤
∥
A
†
∥
∥
A
∥
12
∥
𝔼
n
∥
12
.
Thus, ∥
V
n
-
A
†
∥
→
0; that is, the sequence of (7) converges to the Moore-Penrose inverse in twelfth order as n
→
+
∞. This ends the proof.


## 3. Stability

We investigate the stability of (7) for finding A
† (or the simplified case A
-
1) in a neighborhood of the solution of equation A
V
=
I. Note that if the iteration is not self-correcting, that is, if errors made at one stage are not subsequently damped, then the inevitable rounding errors introduced into the iteration may accumulate to the point where they overwhelm the answer. Thus, we should either show that the proposed method is self-correcting or must furnish an analysis showing that rounding errors remain under control. This will be done in what follows. In fact, we analyze how a small perturbation at the nth iterate is amplified or damped along the iterates. Note that this procedure has recently been applied on a general family of methods for matrix inversion in [17].Theorem 5.
The sequence{
V
n
}
n
=
0
n
=
∞ generated by (7) with the initial approximation (19) is asymptotically stable for finding the Moore-Penrose generalized inverse.Proof.
LetΔ
V
n be the numerical perturbation introduced at the nth iterate of (7). Next, one has
(32)
V
~
n
=
V
n
+
Δ
V
n
.
Here, we perform a first-order error analysis; that is, we formally neglect quadratic or higher terms such as (
Δ
V
n
)
2. This formal manipulation is meaningful if Δ
V
n is sufficiently small and further yields to V
n
·
Δ
V
n
≈
Δ
V
n
·
V
n. We have
(33)
V
~
n
+
1
=
1
64
V
~
n
ϕ
(
V
~
n
)
=
1
64
V
~
n
(
816
I
-
4812
A
V
~
n
+
17393
(
A
V
~
n
)
2
-
43044
(
A
V
~
n
)
3
h
h
h
h
h
h
+
77154
(
A
V
~
n
)
4
-
103312
(
A
V
~
n
)
5
+
105039
(
A
V
~
n
)
6
h
h
h
h
h
h
-
81576
(
A
V
~
n
)
7
+
48268
(
A
V
~
n
)
8
-
21516
(
A
V
~
n
)
9
h
h
h
h
h
h
+
7071
(
A
V
~
n
)
10
-
1652
(
A
V
~
n
)
11
+
258
(
A
V
~
n
)
12
h
h
h
h
h
h
-
24
(
A
V
~
n
)
13
+
(
A
V
~
n
)
14
)
=
1
64
(
V
n
+
Δ
V
n
)
×
(
816
I
-
4812
A
(
V
n
+
Δ
V
n
)
+
17393
(
A
(
V
n
+
Δ
V
n
)
)
2
h
h
h
-
43044
(
A
(
V
n
+
Δ
V
n
)
)
3
+
77154
(
A
(
V
n
+
Δ
V
n
)
)
4
h
h
h
-
103312
(
A
(
V
n
+
Δ
V
n
)
)
5
+
105039
(
A
(
V
n
+
Δ
V
n
)
)
6
h
h
h
-
81576
(
A
(
V
n
+
Δ
V
n
)
)
7
+
48268
(
A
(
V
n
+
Δ
V
n
)
)
8
h
h
h
-
21516
(
A
(
V
n
+
Δ
V
n
)
)
9
+
7071
(
A
(
V
n
+
Δ
V
n
)
)
10
h
h
h
-
1652
(
A
(
V
n
+
Δ
V
n
)
)
11
+
258
(
A
(
V
n
+
Δ
V
n
)
)
12
h
h
h
-
24
(
A
(
V
n
+
Δ
V
n
)
)
13
+
(
A
(
V
n
+
Δ
V
n
)
)
14
)
≈
1
64
(
V
n
+
Δ
V
n
)
×
(
(
(
A
V
n
)
14
+
14
A
(
A
V
n
)
13
Δ
V
n
)
816
I
-
4812
(
A
V
n
+
A
Δ
V
n
)
h
h
h
+
17393
(
(
A
V
n
)
2
+
2
A
(
A
V
n
)
1
Δ
V
n
)
h
h
h
-
43044
(
(
A
V
n
)
3
+
3
A
(
A
V
n
)
2
Δ
V
n
)
h
h
h
+
77154
(
(
A
V
n
)
4
+
4
A
(
A
V
n
)
3
Δ
V
n
)
h
h
h
-
103312
(
(
A
V
n
)
5
+
5
A
(
A
V
n
)
4
Δ
V
n
)
h
h
h
+
105039
(
(
A
V
n
)
6
+
6
A
(
A
V
n
)
5
Δ
V
n
)
h
h
h
-
81576
(
(
A
V
n
)
7
+
7
A
(
A
V
n
)
6
Δ
V
n
)
h
h
h
+
48268
(
(
A
V
n
)
8
+
8
A
(
A
V
n
)
7
Δ
V
n
)
h
h
h
-
21516
(
(
A
V
n
)
9
+
9
A
(
A
V
n
)
8
Δ
V
n
)
h
h
h
+
7071
(
(
A
V
n
)
10
+
10
A
(
A
V
n
)
9
Δ
V
n
)
h
h
h
-
1652
(
(
A
V
n
)
11
+
11
A
(
A
V
n
)
10
Δ
V
n
)
h
h
h
+
258
(
(
A
V
n
)
12
+
12
A
(
A
V
n
)
11
Δ
V
n
)
h
h
h
-
24
(
(
A
V
n
)
13
+
13
A
(
A
V
n
)
12
Δ
V
n
)
h
h
h
+
(
(
A
V
n
)
14
+
14
A
(
A
V
n
)
13
Δ
V
n
)
)
≈
V
n
+
1
+
51
4
V
n
-
1203
16
(
A
V
n
)
1
V
n
+
17393
64
(
A
V
n
)
2
V
n
-
10761
16
(
A
V
n
)
3
V
n
+
38577
32
(
A
V
n
)
4
V
n
-
6457
4
(
A
V
n
)
5
V
n
+
105039
64
(
A
V
n
)
6
V
n
-
10197
8
(
A
V
n
)
7
V
n
+
12067
16
(
A
V
n
)
8
V
n
-
5379
16
(
A
V
n
)
9
V
n
+
7071
64
(
A
V
n
)
10
V
n
-
413
16
(
A
V
n
)
11
V
n
+
129
32
(
A
V
n
)
12
V
n
-
3
8
(
A
V
n
)
13
V
n
+
1
64
(
A
V
n
)
14
V
n
+
51
4
Δ
V
n
-
1203
8
A
V
n
Δ
V
n
+
52179
64
(
A
V
n
)
2
Δ
V
n
-
10761
4
(
A
V
n
)
3
Δ
V
n
+
192885
32
(
A
V
n
)
4
Δ
V
n
-
19371
2
(
A
V
n
)
5
Δ
V
n
+
735273
64
(
A
V
n
)
6
Δ
V
n
-
10197
(
A
V
n
)
7
Δ
V
n
+
108603
16
(
A
V
n
)
8
Δ
V
n
-
26895
8
(
A
V
n
)
9
Δ
V
n
+
77781
64
(
A
V
n
)
10
Δ
V
n
-
1239
4
(
A
V
n
)
11
Δ
V
n
+
1677
32
(
A
V
n
)
12
Δ
V
n
-
21
4
(
A
V
n
)
13
Δ
V
n
+
15
64
(
A
V
n
)
14
Δ
V
n
,
where (
Δ
V
n
)
i
≈
0, i
≥
2 has been used, since they are very close to the zero (matrix). After some algebraic manipulation and using Δ
V
n
+
1
≈
V
~
n
+
1
-
V
n
+
1, we have
(34)
Δ
V
n
+
1
≈
3
64
(
4
I
-
A
V
n
)
2
(
I
-
A
V
n
)
11
(
17
I
-
5
A
V
n
)
Δ
V
n
≈
3
64
[
(
3
I
+
I
-
A
V
n
)
2
(
I
-
A
V
n
)
11
(
12
I
+
5
(
I
-
A
V
n
)
)
]
Δ
V
n
,
≈
3
64
[
108
(
I
-
A
A
†
)
11
+
117
(
I
-
A
A
†
)
12
+
42
(
I
-
A
A
†
)
13
+
5
(
I
-
A
A
†
)
14
]
Δ
V
n
,
by using the fact that for enough large n, we have V
n
≈
A
†. We attain
(35)
∥
Δ
V
n
+
1
∥
≤
3
64
[
108
∥
I
-
A
A
†
∥
11
+
117
∥
I
-
A
A
†
∥
12
+
42
∥
I
-
A
A
†
∥
13
+
5
∥
I
-
A
A
†
∥
14
]
∥
Δ
V
n
∥
.
From (35), we can conclude that the perturbation at the iterate n
+
1, is bounded. Therefore, the sequence {
V
n
}
n
=
0
n
=
∞ generated by (7) is asymptotically stable. This ends the proof.Corollary 6.
By using the matrix identity in (27), it would be possible to further simplify bound (35) as follows:
(36)
∥
Δ
V
n
+
1
∥
≤
(
51
4
∥
I
-
A
A
†
∥
)
∥
Δ
V
n
∥
.
Using (36) recursively, one may attain the following very simple bound:
(37)
∥
Δ
V
n
+
1
∥
≤
[
(
51
4
)
n
+
1
∥
I
-
A
A
†
∥
]
∥
Δ
V
0
∥
.Remark 7.
In case of finding the regular inverse of nonsingular matrices, that is, whenA
†
=
A
-
1, according to (37), we have Δ
V
n
+
1
≈
0, and so the matrix method is strongly numerically stable. Consequently, in case of finding the A
†, the matrix method (7), is asymptotically stable. Of course, since the iteration is not self-correcting in the general case, proceeding beyond convergence may cause a serious increase in error.


## 4. Initial Value

The iterative methods that have been discussed up to now are sensitive upon choosing the initial guess/value to start the process. As a matter of fact, the high accuracy and efficiency of such types of iterative algorithms are guaranteed only if the initial value satisfies the appropriate condition given in Theorem1. Thus, in order to preserve the convergence order, we present some ways from the literature to remedy this flaw, although an efficient way for square or rectangular matrices is the way (19).For a symmetric positive definite (SPD) matrixA, one can easily use the Householder-John Theorem to attain
(38)
V
0
=
P
-
1
as the initial value, wherein the matrix P is any of the matrices such that P
+
P
T
-
A is SPD [8].If the square matrixA is diagonally dominant, one may apply the approach given in [18] and use
(39)
V
0
=
diag
⁡
(
1
a
11
,
1
a
22
,
…
,
1
a
n
n
)
,
wherein a
i
i is the diagonal entry of A. Note that this choice is so much fruitful in solving PDEs resulting from discretizations. Some further generalizations of such an initial matrix are given in [19].Although the two abovementioned ways are efficient, they cannot be applied for finding an initial guess/value to the inverse of general input matrices. For instance, for large-scale matrices which do not satisfy the above structures, they may fail to provide the convergence. Hence, we here take into account the suboptimal way of producingV
0 as given by Pan and Schreiber in [20] as follows:
(40)
V
0
=
A
T
∥
A
∥
1
∥
A
∥
∞
.We should note that choosing the initial value as given above can satisfy the necessary condition of arriving to the convergence phase. Some ways for updating the initial matrix for sparse matrices are brought forward by [21].In what follows, we provide an algorithm for improving an initial matrix for square matrices rapidly. In fact, the derivation of LU factorization for almost all kinds of square nonsingular matrices could be done rapidly in the linear algebra programming packages. In the Mathematica, the one argument commandLinearSolve[ ] provides an LU factorization of A too much fast. Then, by applying the LU decomposition on the columns of a identity matrix recursively, one could update the columns of a derived initial matrix produced by other strategies such as (40). Such a procedure is illustrated in Algorithm 1. We summarized this idea as in Algorithm 2.Algorithm 1:An algorithm for constructing a rapid and robust initial approximation for A
-
1.
(1)Given A
∈
ℂ
m
×
m
(2) construct the matrixV
0
=
(
1
/
(
∥
A
∥
F
2
)
)
A
*
(3) obtain the LU decomposition of the matrixA, for example L
U = LinearSolve[A]
(4) forl
=
1
,
…
,
j
(
j
≪
n
), and when e
l is the lth column of the identity matrix
(5) update the columns ofV
0 by finding the columns of the real A
-
1 as follows: v
l
=
LU
[
e
l
]
(
6
)
end for.Algorithm 2:
Two-argument function written in the Mathematica environment.
i
nitial1
[
A_, num_
]
∶
=
Quiet@Module
[
{
n = Dimensions
[
A
]
[
[
1
]
]
, i = 1, LU = LinearSolve
[
A
]
}
,
Id = SparseArray
[
{
{
k_, k_}
-> 1.},{n, n
}
]
;
mat = (1/Norm
[
A,“Frobenius”]
∧
2
)*ConjugateTranspose
[
A
]
;
While
[
i
<
= num,{ith = LU
[
Id
[
[
All, n + 1 - i
]
]
]
;
mat
[
[
All, n + 1 - i
]
]
= ith; i++;
}
]
; mat
]
;
initial2
[
A_, num_
]
∶
=
Quiet@Module
[
{
n = Dimensions
[
A
]
[
[
1
]
]
, i = 1},
Id = SparseArray
[
{
{
k_, k
_}
-
> 1
.},{n, n}, 0
]
; LU = LinearSolve
[
A
]
;
mat = Quiet[initial1
[
A, num
]
]
;
While
[
i <= num,{ith = LU
[
Id
[
[
All, i
]
]
]
;
mat
[
[
All, i
]
]
= ith; i++;
}
]
; mat
]
;Theinitial1[A_,num_] takes the nonsingular matrix A and the number of columns that users wish to update from the real inverse into the approximate inverse from the left, while the function initial2[A_,num_] works doubly. That is, if, for example, num=10, it updates the first and last 10 columns of the approximate inverse as rapidly as possible.Next, we conduct some numerical tests to support the theoretical results given in Section2 using the initial value discussed herein.


## 5. Computational Aspects

In this section, some experiments are presented to demonstrate the capability of the proposed method. The programming packagemathematica 8 [22] has been used in the demonstrations. We work on the numerical aspects of the methods in machine precision.It is clear that large sparse matrices cannot be handled easily and their storage needs to be done in sparse form as in the input form to be accessible and economic in real applications. Methods like (7) are powerful in finding an approximate inverse or a robust approximate inverse preconditioner in low number of steps and computational time, in which the output form of the approximate inverse is also sparse.In this paper, as the programs were running, we found the running time using the commandAbsoluteTiming[ ] to report the elapsed CPU time (in seconds) for this experiment. In this paper, the computer specifications are Microsoft Windows XP Intel(R), Pentium(R) 4, and CPU 3.20 GHz, with 4 GB of RAM.Experiment 1.
This test is devoted to the application of the Schulz-type iterative methods in finding the pseudoinverse of 30 large random complex matrices defined as shown in Algorithm3 (I
=
-
1).Algorithm 3
m
= 1500; k = 1800; number = 30; SeedRandom
[
1234
]
;
Table
[
A
[
l
]
= SparseArray
[
{
Band
[
{
400, 10},{m, k
}
]
-> Random
[
]
- I,
Band
[
{
10, 200},{m, k
}
]
->{1.1, -Random
[
]
}
,
Band
[
{
-60, 100
}
]
-> -3.02, Band
[
{
-100, 500
}
]
-> 3.1 I},{m, k}, 0.
]
;,{l, number
}
]
;
Id = SparseArray
[
{
{
i
_
,
i
_}
-> 1.},{m, m}, 0.
]
;The results of comparisons for these random matrices of the sizem
×
k
=
1500
×
1800 are reported in Figures 1, 2, and 3 in terms of the number of iterations and the computational time. The compared methods are (1) denoted by “Schulz,” (2) denoted by “Chebyshev,” (4) denoted by “KMS,” and the new iterative scheme (7) denoted by “PM." A saving in the elapsed time by considering the stopping criteria as ∥
V
n
+
1
-
V
n
∥
∞
≤
1
0
-
6 and ∥
V
n
+
1
-
V
n
∥
F
≤
1
0
-
6 can be observed for the studied method (7). In this test, the initial matrix has been computed for each random matrix by V0[j]=ConjugateTranspose[A[j]] *(1./((SingularValueList[A[j],1][[1]])2)), while the maximum number of iterations is set to 100.Figure 1
The results of comparisons for Experiment1 in terms of the number of iterations.Figure 2
The results of comparisons for Experiment1 in terms of the elapsed time using ∥
V
n
+
1
-
V
n
∥
∞
≤
1
0
-
6.Figure 3
The results of comparisons for Experiment1 in terms of the elapsed time using ∥
V
n
+
1
-
V
n
∥
F
≤
1
0
-
6.


## 6. Concluding Remarks

It is well known that matrix inverse and generalized inverse are important in applied fields of nature science, such as the solution to various systems of linear and nonlinear equations, eigenvalue problems, and the linear least square problems. Iterative methods are often effective especially for large-scale systems with sparsity and Schulz-type methods are great tools for preconditioning such systems or in finding pseudoinverses.Hotelling-Bodewig algorithm is simple to describe and analyze and is numerically stable. This was the idea of developing an iterative method of this type in this paper.In this paper, we have shown that the suggested method (7) reaches twelfth order of convergence. The stability of the new method was also studied in detail and established that the new scheme is asymptotically stable. The efficacy of the new scheme was illustrated numerically in Section 5. Finally, and based on the numerical results obtained, one can conclude that the presented method is useful. Further extensions of the new scheme for other generalized inverses (such as the ones in [23, 24]) can be done for future works.



---

*Source: 830564-2014-02-06.xml*