# Improved Qrginv Algorithm for Computing Moore-Penrose Inverse Matrices

Authors: Alireza Ataei
Journal: ISRN Applied Mathematics
 (2014)

---

## Abstract

Katsikis et al. presented a computational method in order to calculate the Moore-Penrose inverse of an arbitrary matrix (including singular and rectangular) (2011). In this paper, an improved version of this method is presented for computing the pseudo inverse of anm
×
n real matrix A with rank r
>
0. Numerical experiments show that the resulting pseudoinverse matrix is reasonably accurate and its computation time is significantly less than that obtained by Katsikis et al.


---

## Body

## 1. Introduction

Letℝ
m
×
n denote the set of all m
×
n matrices over the field of real numbers, ℝ. The symbols A
T,rank(
A
T
) will stand for the transpose and rank of A
∈
ℝ
m
×
n, respectively. For a matrix A
∈
ℝ
m
×
n, the Moore-Penrose inverse of A, denoted by A
†, is the unique matrix X
∈
ℝ
n
×
m satisfying the following equations:(i)
A
X
A
=
A,
(ii)
X
A
X
=
X,
(iii)
(
A
X
)
T
=
A
X,
(iv)
(
X
A
)
T
=
X
A.A lot of works concerning generalized inverses have been carried out, in finite and infinite dimensions (e.g., [1–3]). There are several methods for computing the Moore-Penrose inverse matrix [2, 4–8]. In a recent article [9], an improved method for the computation of the Moore-Penrose inverse matrix provided.In this paper, we aim to improve their method so that it can be used for any kind of matrices, square or rectangular, full rank or not. The numerical examples show that our method is competitive in terms of accuracy and is much faster than the commonly used methods and can also be used for large sparse matrices.This paper is organized as follows. In Section2 the improved version of this method is presented for computing the pseudoinverse of an m
×
n real matrix A with rank r
>
0. In Section 3, the numerical results of some test matrices are given. Section 4 is devoted to the concluding remarks.


## 2. Improved Qrginv Method (IMqrginv)

In [9], a method (qrginv) for computing the Moore-Penrose inverse of an arbitrary matrix was presented. They made use of the QR-factorization, as well as an algorithm based on a known reverse order law for generalized inverse matrices, and also they apply a method (ginv), presented in [4], based on a full rank Cholesky factorization of possibly singular symmetric positive matrices.In the current paper, we improvedqrginv algorithm using the QR-factorization by Gram-Schmidt orthonormalization (GSO) and Theorem 1 for faster computing Moore-Penrose inverse of arbitrary matrices (including singular and rectangular). We should note that we invoke ginv algorithm. In order to support and state our achievement we need to remind Gram-Schmidt orthonormalization (GSO) and the QR-factorization.


### 2.1. The Gram-Schmidt Procedure

Let us remind a generalization of the Gram-Schmidt orthonormalization process (shortly GSO) which is applied for singular matrices. Let𝒜
=
{
a
1
,
a
2
,
…
,
a
n
}
⊆
ℝ
m be a set of vectors spanning a subspace V, this process generates a set of mutually orthonormal vectors such as 𝒬
=
{
q
1
,
q
2
,
…
,
q
r
}
⊆
ℝ
m having the property that 𝒬 is an orthonormalization basis for V. 𝒬 is obtained using the Gram-Schmidt orthonormalization process (shortly GSO) as follows:
(1)
q
1
=
a
c
1
∥
a
c
1
∥
,
if
a
c
1
≠
0
=
a
j
for
1
≤
j
<
c
1
,
q
^
j
=
a
j
-
∑
i
=
1
k
-
1
‍
(
a
j
,
q
i
)
q
i
,
j
=
c
k
-
1
+
1
,
c
k
-
2
+
2
,
…
,
c
k
,
q
k
=
q
^
j
∥
q
^
j
∥
,
if
q
^
c
k
≠
0
=
q
j
,
for
c
k
-
1
+
1
≤
j
<
c
k
,
k
=
2
,
…
,
r
.The integerr found by the GSO process is the dimension of the subspace V. The integers {
c
1
,
…
,
c
r
} are the indices of a maximal linearly independent subset {
a
c
1
,
…
,
a
c
r
} of 𝒜.


### 2.2. The QR-Factorization

Let us remind the QR-factorization for arbitrary matrices (including singular and rectangular).Let the orthonormal set{
q
1
,
q
2
,
…
,
q
r
} be obtained from the set of vectors {
a
1
,
a
2
,
…
,
a
n
} by the GSO process described in Section 2.1, and let
(2)
Q
~
=
[
q
1
,
q
2
,
…
,
q
r
]
∈
ℝ
m
×
r
,
A
=
[
a
1
,
a
2
,
…
,
a
n
]
∈
ℝ
m
×
n
,
with rank
⁡
(
A
)
=
r
>
0, be the corresponding matrices. Then there exist matrices A
~, Q, and R such that
(3)
A
~
=
A
P
=
Q
R
,
where (i)
P is a permutation matrix (therefore orthonormal);
(ii)
Q
=
[
Q
~
Z
]
∈
ℝ
m
×
m, where Q
~ and Z denote matrices whose columns are an orthonormal basis of R
(
A
) and N
(
A
T
), respectively;
(iii)
R
=
[
R
~
0
]
∈
ℝ
m
×
n, where R
~
∈
ℝ
r
×
n is upper triangular matrix with rank
⁡
(
R
~
)
=
r. One obtains from (3)
(4)
A
~
=
[
Q
~
Z
]
[
R
~
0
]
=
Q
~
R
~
.
It follows that A
~ has a Q
~
R
~-factorization. A nonzero matrix can be expressed as the product of a matrix of full column rank and a matrix of full row rank. In fact, for given A
∈
ℝ
m
×
n (rank
⁡
(
A
)
=
r
>
0) there exist matrices F
∈
ℝ
m
×
r and G
∈
ℝ
r
×
n such that A
=
F
G [2]. Such factorization, which is the so-called a full rank factorization, turns out to be a powerful tool in the study of generalized inverses.The following theorem is due to MacDuffe [10] who apparently was the first one to point out that a full rank factorization of a matrix A leads to an explicit formula for its Moore-Penrose inverse, A
†.Theorem 1.
IfA
∈
ℝ
m
×
n matrix, with rank
⁡
(
A
)
=
r
>
0, has a full rank factorization A
=
F
G, then
(5)
A
†
=
G
T
(
F
T
A
G
T
)
-
1
F
T
.As a direct consequence of Theorem1, we have the following.Corollary 2.
LetA
∈
ℝ
m
×
n, rank
⁡
(
A
)
=
r
>
0, and A
~
=
Q
~
R
~ be the Q
~
R
~-factorization of A
~. Then
(6)
A
†
=
P
T
(
R
~
T
)
(
R
~
R
~
T
)
-
1
Q
~
T
.Proof.
WithA
=
A
~, F
=
Q
~, and G
=
R
~ in Theorem 1 one obtains
(7)
(
A
~
)
†
=
(
A
P
)
†
=
R
~
T
(
R
~
R
~
T
)
-
1
Q
~
T
,
and consequently
(8)
A
†
=
P
R
~
T
(
R
~
R
~
T
)
-
1
Q
~
T
.The functionIMqrginv 4 provided all implementation details of the above solution, in MATLAB code. To calculate the rank of R
~, one needs only the number of its columns having at least one value above a tolerance level in absolute terms. This tolerance is set to be equal to 1
0
-
5, which is also used by Katsikis et al. [9], and turns out to provide accurate results.


## 2.1. The Gram-Schmidt Procedure

Let us remind a generalization of the Gram-Schmidt orthonormalization process (shortly GSO) which is applied for singular matrices. Let𝒜
=
{
a
1
,
a
2
,
…
,
a
n
}
⊆
ℝ
m be a set of vectors spanning a subspace V, this process generates a set of mutually orthonormal vectors such as 𝒬
=
{
q
1
,
q
2
,
…
,
q
r
}
⊆
ℝ
m having the property that 𝒬 is an orthonormalization basis for V. 𝒬 is obtained using the Gram-Schmidt orthonormalization process (shortly GSO) as follows:
(1)
q
1
=
a
c
1
∥
a
c
1
∥
,
if
a
c
1
≠
0
=
a
j
for
1
≤
j
<
c
1
,
q
^
j
=
a
j
-
∑
i
=
1
k
-
1
‍
(
a
j
,
q
i
)
q
i
,
j
=
c
k
-
1
+
1
,
c
k
-
2
+
2
,
…
,
c
k
,
q
k
=
q
^
j
∥
q
^
j
∥
,
if
q
^
c
k
≠
0
=
q
j
,
for
c
k
-
1
+
1
≤
j
<
c
k
,
k
=
2
,
…
,
r
.The integerr found by the GSO process is the dimension of the subspace V. The integers {
c
1
,
…
,
c
r
} are the indices of a maximal linearly independent subset {
a
c
1
,
…
,
a
c
r
} of 𝒜.


## 2.2. The QR-Factorization

Let us remind the QR-factorization for arbitrary matrices (including singular and rectangular).Let the orthonormal set{
q
1
,
q
2
,
…
,
q
r
} be obtained from the set of vectors {
a
1
,
a
2
,
…
,
a
n
} by the GSO process described in Section 2.1, and let
(2)
Q
~
=
[
q
1
,
q
2
,
…
,
q
r
]
∈
ℝ
m
×
r
,
A
=
[
a
1
,
a
2
,
…
,
a
n
]
∈
ℝ
m
×
n
,
with rank
⁡
(
A
)
=
r
>
0, be the corresponding matrices. Then there exist matrices A
~, Q, and R such that
(3)
A
~
=
A
P
=
Q
R
,
where (i)
P is a permutation matrix (therefore orthonormal);
(ii)
Q
=
[
Q
~
Z
]
∈
ℝ
m
×
m, where Q
~ and Z denote matrices whose columns are an orthonormal basis of R
(
A
) and N
(
A
T
), respectively;
(iii)
R
=
[
R
~
0
]
∈
ℝ
m
×
n, where R
~
∈
ℝ
r
×
n is upper triangular matrix with rank
⁡
(
R
~
)
=
r. One obtains from (3)
(4)
A
~
=
[
Q
~
Z
]
[
R
~
0
]
=
Q
~
R
~
.
It follows that A
~ has a Q
~
R
~-factorization. A nonzero matrix can be expressed as the product of a matrix of full column rank and a matrix of full row rank. In fact, for given A
∈
ℝ
m
×
n (rank
⁡
(
A
)
=
r
>
0) there exist matrices F
∈
ℝ
m
×
r and G
∈
ℝ
r
×
n such that A
=
F
G [2]. Such factorization, which is the so-called a full rank factorization, turns out to be a powerful tool in the study of generalized inverses.The following theorem is due to MacDuffe [10] who apparently was the first one to point out that a full rank factorization of a matrix A leads to an explicit formula for its Moore-Penrose inverse, A
†.Theorem 1.
IfA
∈
ℝ
m
×
n matrix, with rank
⁡
(
A
)
=
r
>
0, has a full rank factorization A
=
F
G, then
(5)
A
†
=
G
T
(
F
T
A
G
T
)
-
1
F
T
.As a direct consequence of Theorem1, we have the following.Corollary 2.
LetA
∈
ℝ
m
×
n, rank
⁡
(
A
)
=
r
>
0, and A
~
=
Q
~
R
~ be the Q
~
R
~-factorization of A
~. Then
(6)
A
†
=
P
T
(
R
~
T
)
(
R
~
R
~
T
)
-
1
Q
~
T
.Proof.
WithA
=
A
~, F
=
Q
~, and G
=
R
~ in Theorem 1 one obtains
(7)
(
A
~
)
†
=
(
A
P
)
†
=
R
~
T
(
R
~
R
~
T
)
-
1
Q
~
T
,
and consequently
(8)
A
†
=
P
R
~
T
(
R
~
R
~
T
)
-
1
Q
~
T
.The functionIMqrginv 4 provided all implementation details of the above solution, in MATLAB code. To calculate the rank of R
~, one needs only the number of its columns having at least one value above a tolerance level in absolute terms. This tolerance is set to be equal to 1
0
-
5, which is also used by Katsikis et al. [9], and turns out to provide accurate results.


## 3. Numerical Examples

In this section, we compare the performance of the proposed method (IMqrginv function) to that of Katsikis et al. [9] for the computation of Moore-Penrose inverse matrices. Testingqrginv andIMqrginv was performed separately for random singular and singular matrices with “large” condition number from the Matrix Computation Toolbox (see [11]). We also tested the proposed method for some sparse matrices and obtained very fast and accurate results. Specifically, MATLAb 7.11 (R2010b) Service Pack 3 version of the software was used on an Intel Core 2 (Duo) 8400 Processor running 2.26 GHz with 3 GB of RAM memory using Windows Vista professional 32-bit operating system.


### 3.1. Random Singular Matrices

We are computing the performance of the proposed methodIMqrginv to that of [9] (qrginv function). In the same way of [4] we tested on a series of random singular matrices of size m
×
n, with n
=
2
k
,
k
=
7
,
…
,
11, and m
=
2
n, which are rank deficient, with rank r
=
7
n
/
8. In addition, the accuracy of the results is examined with the matrix 2-norm in error matrices corresponding to the four properties characterizing the Moore-Penrose inverses shown in Table 1. The computation error is less than 1
0
-
12 per coefficient in the error matrices, in all cases. The computation time (in seconds) is reported in Table 1. We observe that the computation time ofIMqrginv method is substantially less than that of theqrginv method.Table 1
Error and computational time results for random singular matrices.
n
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
2
7
qrginv
0.0317
3.4169
e
-
014
0.2638
e
-
016
5.7671
e
-
015
3.2645
e
-
015
IMqrginv
0.0137
1.6175
e
-
014
3.7297
e
-
016
3.0344
e
-
015
2.5228
e
-
015
2
8
qrginv
0.1176
1.0052
e
-
013
9.9243
e
-
016
8.0286
e
-
015
4.0476
e
-
015
IMqrginv
0.0786
2.9804
e
-
014
3.4041
e
-
016
3.8500
e
-
015
3.1148
e
-
015
2
9
qrginv
1.0584
1.3649
e
-
013
1.0186
e
-
015
1.1006
e
-
014
4.9494
e
-
015
IMqrginv
0.8236
5.3456
e
-
014
5.0122
e
-
015
5.0122
e
-
015
3.8394
e
-
015
2
10
qrginv
10.8336
2.6266
e
-
013
1.0616
e
-
015
1.4943
e
-
014
5.9069
e
-
015
IMqrginv
8.2976
8.8502
e
-
014
2.9117
e
-
016
6.3837
e
-
015
4.8854
e
-
015
2
11
qrginv
81.1429
4.9931
e
-
013
1.1850
e
-
015
2.1858
e
-
014
7.2568
e
-
015
IMqrginv
55.0746
1.6676
e
-
013
2.5936
e
-
016
8.3743
e
-
015
6.5585
e
-
015


### 3.2. Singular Matrices

In this section we use a set of singular test matrices that includes singular matrices of size200
×
200, obtained from the function matrix in the Matrix Computation Toolbox [11] (which inculdes test matrices from Matlab itself). The condition number of test matrices ranges from order 1
0
16 to 1
0
135. Since the matrices are of relatively small size and so as to measure the time needed for each algorithm to compute the Moore-Penrose inverse accurately, each algorithm runs 100 distinct times. The reported time is the mean time over these 100 replications. For each case, the time responses together with the error results are presented in Tables 2, 3, 4, 5, 6, 7, 8, and 9. We observe that the computation time ofIMqrginv method is shorter than theqrginv method for all matrices and is proved to be more efficient.Table 2
Error and computational time results forchow (rank = 199, cond = 5.9407
e
+
135).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0376
3.6711
e
-
013
1.7331
e
-
013
4.4765
e
-
013
2.4702
e
-
013
IMqrginv
0.0187
4.0038
e
-
013
1.7331
e
-
013
2.4448
e
-
013
2.4741
e
-
013Table 3
Error and computational time results forcycol (rank = 50, cond = 3.7045
e
+
048).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.2816
1.5114
e
-
013
2.3446
e
-
017
1.1512
e
-
015
8.8889
e
-
016
IMqrginv
0.2334
8.1308
e
-
014
1.4034
e
-
017
1.0225
e
-
015
8.2643
e
-
016Table 4
Error and computational time results forgearmat (rank = 199, cond = 1.8074
e
+
016).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0377
4.6584
e
-
015
3.0532
e
-
013
1.1996
e
-
013
2.1253
e
-
014
IMqrginv
0.0165
2.8959
e
-
015
3.3357
e
-
013
7.7888
e
-
014
2.1380
e
-
014Table 5
Error and computational time results forkahan (rank = 199, cond = 1.9055
e
+
024).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0279
1.9877
e
-
005
4.1177
e
-
009
8.8330
e
-
001
5.4162
e
-
014
IMqrginv
0.0099
1.9877
e
-
005
3.8389
e
-
009
8.8330
e
-
001
1.0398
e
-
014Table 6
Error and computational time results forlotkin (rank = 19, cond = 3.8210
e
+
021).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0244
8.2512
e
-
006
1.2717
e
-
011
4.4898
e
-
002
3.9689
e
-
009
IMqrginv
0.0102
8.2512
e
-
006
3.2435
e
-
009
4.4898
e
-
002
1.2636
e
-
011Table 7
Error and computational time results forprolate (rank = 117, cond = 4.7489
e
+
017).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0276
1.3837
e
-
006
1.2998
e
-
007
4.7715
e
-
002
4.7317
e
-
011
IMqrginv
0.0138
1.3837
e
-
006
1.1842
e
-
007
4.7715
e
-
002
4.7401
e
-
011Table 8
Error and computational time results forhilb (rank = 20, cond = 4.4158
e
+
021).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0247
7.7880
e
-
006
1.1444
e
-
008
1.0053
e
-
001
5.5974
e
-
012
IMqrginv
0.0117
7.7880
e
-
006
1.1184
e
-
008
1.0053
e
-
001
5.5636
e
-
012Table 9
Error and computational time results formagic (rank = 3, cond = 5.3481
e
+
019).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0154
1.5837
e
-
009
4.4922
e
-
009
9.4641
e
-
014
6.3953
e
-
015
IMqrginv
0.0035
1.4929
e
-
009
4.8349
e
-
009
4.7537
e
-
014
6.0546
e
-
015


### 3.3. Matrix-Market Sparse Matrices

For sparse matrices, we have chosen some matrices from Matrix-Market collection [11]. We follow the same method as in [8], and we have the rank deficient matrices as
(9)
A
-
Z
≡
[
A
Z
]
,
where A is one of the chosen matrices and Z is a zero matrix of order m
×
100. The chosen matrices with their properties are shown in Table 10. The results of the methods are presented in Table 11. We observe that the Moore-Penrose inverses obtained byIMqrginv are reasonably accurate in all cases; the computation time required by theIMqrginv method is significantly less than the time required by theIMqrginv methods. On the other hand, we can see that the accuracy computation of theIMqrginv method is less than theqrginv method; however, in some cases, the accuracy of the results of both methods is low. We can conclude thatIMqrginv method is a robust and efficient tool for obtaining the Moore-Penrose inverse of large sparse and rank deficient matrices.Table 10
Test problem information.
Matrix∖property
m
n
nnz
Cond
WELL1033
1033
320
4732
Not available
WELL1850
1850
712
8758
Not available
ILCC1033
1033
320
4732
Not available
ILCC1850
10850
712
8758
Not available
WATT1
1856
1856
11360
5.4
e
+
009
GR-30-30
900
900
4322
3.8
e
+
002
ADD20
2395
2395
17319
1.76
e
+
004
NOSE3
960
960
8402
7.3
e
+
004
SHERMAN1
1000
1000
3750
2.3
e
+
004Table 11
Error and computational time results formagic (rank = 3, cond = 5.3481
e
+
019).
Method
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
WELL1033_Z
qrginv
0.6277
2.1186
e
-
013
3.6375
e
-
011
1.4406
e
-
011
1.8819
e
-
013
IMqrginv
0.3300
3.1283
e
-
014
2.5635
e
-
011
2.4875
e
-
012
9.1501
e
-
014
WELL1850_Z
qrginv
3.8271
1.0218
e
-
012
5.2163
e
-
011
5.6680
e
-
011
3.5487
e
-
013
IMqrginv
1.7559
4.0066
e
-
014
6.3726
e
-
012
1.9053
e
-
012
7.7633
e
-
014
ILCC1033_Z
qrginv
0.6119
1.4956
e
-
010
1.1845
e
-
005
1.1674
e
-
006
6.9918
e
-
011
IMqrginv
0.3786
2.3305
e
-
011
8.1774
e
-
006
1.5766
e
-
008
5.6012
e
-
010
ILCC1850_Z
qrginv
3.9197
1.9075
e
-
011
1.9254
e
-
008
9.4935
e
-
009
1.1185
e
-
011
IMqrginv
1.6791
2.2511
e
-
013
9.5637
e
-
009
1.2945
e
-
010
6.6275
e
-
012
WATT1_Z
qrginv
10.5110
7.1357
e
-
007
1.8762
e
-
005
6.1068
7.4875
e
-
009
IMqrginv
3.1018
7.1357
e
-
007
6.4114
e
-
006
6.1068
7.4875
e
-
009
GR-30-30_Z
qrginv
3.1095
2.3077
e
-
011
5.3574
e
-
011
3.5492
e
-
010
3.7356
e
-
013
IMqrginv
1.6834
3.2708
e
-
013
1.2685
e
-
011
4.9654
e
-
012
2.3249
e
-
013
ADD20_Z
qrginv
60.9176
2.1032
e
-
008
4.9866
e
-
004
2.7022
e
-
004
1.4014
e
-
008
IMqrginv
59.0933
9.3835
e
-
011
2.0676
e
-
005
8.6123
e
-
007
1.0353
e
-
009
NOS3_Z
qrginv
4.3836
1.6638
e
-
009
2.6309
e
-
009
5.4663
e
-
008
2.1929
e
-
011
IMqrginv
3.0214
8.6964
e
-
011
8.4269
e
-
011
9.9464
e
-
010
2.1778
e
-
011
SHERMAN1_Z
qrginv
4.3150
1.0363
e
-
009
4.5703
e
-
007
7.4492
e
-
007
3.0278
e
-
010
IMqrginv
1.3069
1.6398
e
-
012
3.2617
e
-
007
1.4722
e
-
009
5.6796
e
-
012


## 3.1. Random Singular Matrices

We are computing the performance of the proposed methodIMqrginv to that of [9] (qrginv function). In the same way of [4] we tested on a series of random singular matrices of size m
×
n, with n
=
2
k
,
k
=
7
,
…
,
11, and m
=
2
n, which are rank deficient, with rank r
=
7
n
/
8. In addition, the accuracy of the results is examined with the matrix 2-norm in error matrices corresponding to the four properties characterizing the Moore-Penrose inverses shown in Table 1. The computation error is less than 1
0
-
12 per coefficient in the error matrices, in all cases. The computation time (in seconds) is reported in Table 1. We observe that the computation time ofIMqrginv method is substantially less than that of theqrginv method.Table 1
Error and computational time results for random singular matrices.
n
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
2
7
qrginv
0.0317
3.4169
e
-
014
0.2638
e
-
016
5.7671
e
-
015
3.2645
e
-
015
IMqrginv
0.0137
1.6175
e
-
014
3.7297
e
-
016
3.0344
e
-
015
2.5228
e
-
015
2
8
qrginv
0.1176
1.0052
e
-
013
9.9243
e
-
016
8.0286
e
-
015
4.0476
e
-
015
IMqrginv
0.0786
2.9804
e
-
014
3.4041
e
-
016
3.8500
e
-
015
3.1148
e
-
015
2
9
qrginv
1.0584
1.3649
e
-
013
1.0186
e
-
015
1.1006
e
-
014
4.9494
e
-
015
IMqrginv
0.8236
5.3456
e
-
014
5.0122
e
-
015
5.0122
e
-
015
3.8394
e
-
015
2
10
qrginv
10.8336
2.6266
e
-
013
1.0616
e
-
015
1.4943
e
-
014
5.9069
e
-
015
IMqrginv
8.2976
8.8502
e
-
014
2.9117
e
-
016
6.3837
e
-
015
4.8854
e
-
015
2
11
qrginv
81.1429
4.9931
e
-
013
1.1850
e
-
015
2.1858
e
-
014
7.2568
e
-
015
IMqrginv
55.0746
1.6676
e
-
013
2.5936
e
-
016
8.3743
e
-
015
6.5585
e
-
015


## 3.2. Singular Matrices

In this section we use a set of singular test matrices that includes singular matrices of size200
×
200, obtained from the function matrix in the Matrix Computation Toolbox [11] (which inculdes test matrices from Matlab itself). The condition number of test matrices ranges from order 1
0
16 to 1
0
135. Since the matrices are of relatively small size and so as to measure the time needed for each algorithm to compute the Moore-Penrose inverse accurately, each algorithm runs 100 distinct times. The reported time is the mean time over these 100 replications. For each case, the time responses together with the error results are presented in Tables 2, 3, 4, 5, 6, 7, 8, and 9. We observe that the computation time ofIMqrginv method is shorter than theqrginv method for all matrices and is proved to be more efficient.Table 2
Error and computational time results forchow (rank = 199, cond = 5.9407
e
+
135).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0376
3.6711
e
-
013
1.7331
e
-
013
4.4765
e
-
013
2.4702
e
-
013
IMqrginv
0.0187
4.0038
e
-
013
1.7331
e
-
013
2.4448
e
-
013
2.4741
e
-
013Table 3
Error and computational time results forcycol (rank = 50, cond = 3.7045
e
+
048).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.2816
1.5114
e
-
013
2.3446
e
-
017
1.1512
e
-
015
8.8889
e
-
016
IMqrginv
0.2334
8.1308
e
-
014
1.4034
e
-
017
1.0225
e
-
015
8.2643
e
-
016Table 4
Error and computational time results forgearmat (rank = 199, cond = 1.8074
e
+
016).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0377
4.6584
e
-
015
3.0532
e
-
013
1.1996
e
-
013
2.1253
e
-
014
IMqrginv
0.0165
2.8959
e
-
015
3.3357
e
-
013
7.7888
e
-
014
2.1380
e
-
014Table 5
Error and computational time results forkahan (rank = 199, cond = 1.9055
e
+
024).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0279
1.9877
e
-
005
4.1177
e
-
009
8.8330
e
-
001
5.4162
e
-
014
IMqrginv
0.0099
1.9877
e
-
005
3.8389
e
-
009
8.8330
e
-
001
1.0398
e
-
014Table 6
Error and computational time results forlotkin (rank = 19, cond = 3.8210
e
+
021).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0244
8.2512
e
-
006
1.2717
e
-
011
4.4898
e
-
002
3.9689
e
-
009
IMqrginv
0.0102
8.2512
e
-
006
3.2435
e
-
009
4.4898
e
-
002
1.2636
e
-
011Table 7
Error and computational time results forprolate (rank = 117, cond = 4.7489
e
+
017).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0276
1.3837
e
-
006
1.2998
e
-
007
4.7715
e
-
002
4.7317
e
-
011
IMqrginv
0.0138
1.3837
e
-
006
1.1842
e
-
007
4.7715
e
-
002
4.7401
e
-
011Table 8
Error and computational time results forhilb (rank = 20, cond = 4.4158
e
+
021).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0247
7.7880
e
-
006
1.1444
e
-
008
1.0053
e
-
001
5.5974
e
-
012
IMqrginv
0.0117
7.7880
e
-
006
1.1184
e
-
008
1.0053
e
-
001
5.5636
e
-
012Table 9
Error and computational time results formagic (rank = 3, cond = 5.3481
e
+
019).
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
qrginv
0.0154
1.5837
e
-
009
4.4922
e
-
009
9.4641
e
-
014
6.3953
e
-
015
IMqrginv
0.0035
1.4929
e
-
009
4.8349
e
-
009
4.7537
e
-
014
6.0546
e
-
015


## 3.3. Matrix-Market Sparse Matrices

For sparse matrices, we have chosen some matrices from Matrix-Market collection [11]. We follow the same method as in [8], and we have the rank deficient matrices as
(9)
A
-
Z
≡
[
A
Z
]
,
where A is one of the chosen matrices and Z is a zero matrix of order m
×
100. The chosen matrices with their properties are shown in Table 10. The results of the methods are presented in Table 11. We observe that the Moore-Penrose inverses obtained byIMqrginv are reasonably accurate in all cases; the computation time required by theIMqrginv method is significantly less than the time required by theIMqrginv methods. On the other hand, we can see that the accuracy computation of theIMqrginv method is less than theqrginv method; however, in some cases, the accuracy of the results of both methods is low. We can conclude thatIMqrginv method is a robust and efficient tool for obtaining the Moore-Penrose inverse of large sparse and rank deficient matrices.Table 10
Test problem information.
Matrix∖property
m
n
nnz
Cond
WELL1033
1033
320
4732
Not available
WELL1850
1850
712
8758
Not available
ILCC1033
1033
320
4732
Not available
ILCC1850
10850
712
8758
Not available
WATT1
1856
1856
11360
5.4
e
+
009
GR-30-30
900
900
4322
3.8
e
+
002
ADD20
2395
2395
17319
1.76
e
+
004
NOSE3
960
960
8402
7.3
e
+
004
SHERMAN1
1000
1000
3750
2.3
e
+
004Table 11
Error and computational time results formagic (rank = 3, cond = 5.3481
e
+
019).
Method
Method
Time
∥
A
A
†
A
-
A
∥
2
∥
A
†
A
A
†
-
A
†
∥
2
∥
(
A
A
†
)
T
-
(
A
A
†
)
∥
2
∥
(
A
†
A
)
T
-
(
A
†
A
)
∥
2
WELL1033_Z
qrginv
0.6277
2.1186
e
-
013
3.6375
e
-
011
1.4406
e
-
011
1.8819
e
-
013
IMqrginv
0.3300
3.1283
e
-
014
2.5635
e
-
011
2.4875
e
-
012
9.1501
e
-
014
WELL1850_Z
qrginv
3.8271
1.0218
e
-
012
5.2163
e
-
011
5.6680
e
-
011
3.5487
e
-
013
IMqrginv
1.7559
4.0066
e
-
014
6.3726
e
-
012
1.9053
e
-
012
7.7633
e
-
014
ILCC1033_Z
qrginv
0.6119
1.4956
e
-
010
1.1845
e
-
005
1.1674
e
-
006
6.9918
e
-
011
IMqrginv
0.3786
2.3305
e
-
011
8.1774
e
-
006
1.5766
e
-
008
5.6012
e
-
010
ILCC1850_Z
qrginv
3.9197
1.9075
e
-
011
1.9254
e
-
008
9.4935
e
-
009
1.1185
e
-
011
IMqrginv
1.6791
2.2511
e
-
013
9.5637
e
-
009
1.2945
e
-
010
6.6275
e
-
012
WATT1_Z
qrginv
10.5110
7.1357
e
-
007
1.8762
e
-
005
6.1068
7.4875
e
-
009
IMqrginv
3.1018
7.1357
e
-
007
6.4114
e
-
006
6.1068
7.4875
e
-
009
GR-30-30_Z
qrginv
3.1095
2.3077
e
-
011
5.3574
e
-
011
3.5492
e
-
010
3.7356
e
-
013
IMqrginv
1.6834
3.2708
e
-
013
1.2685
e
-
011
4.9654
e
-
012
2.3249
e
-
013
ADD20_Z
qrginv
60.9176
2.1032
e
-
008
4.9866
e
-
004
2.7022
e
-
004
1.4014
e
-
008
IMqrginv
59.0933
9.3835
e
-
011
2.0676
e
-
005
8.6123
e
-
007
1.0353
e
-
009
NOS3_Z
qrginv
4.3836
1.6638
e
-
009
2.6309
e
-
009
5.4663
e
-
008
2.1929
e
-
011
IMqrginv
3.0214
8.6964
e
-
011
8.4269
e
-
011
9.9464
e
-
010
2.1778
e
-
011
SHERMAN1_Z
qrginv
4.3150
1.0363
e
-
009
4.5703
e
-
007
7.4492
e
-
007
3.0278
e
-
010
IMqrginv
1.3069
1.6398
e
-
012
3.2617
e
-
007
1.4722
e
-
009
5.6796
e
-
012


## 4. Conclusion

In this paper we have presented a new method, calledIMqrginv, for the fast computation of Moore-Penrose inverse of singular square, rectangular, full, or sparse matrices. This method is based on the GSO method and Theorem 1. Invoking the ginv function we improvedqrginv methods presented by Katsikis et al. [9]. We have compared the performance of the proposed methodIMqrginv to theqrginv method. Numerical examples (see Tables 1–9 and 11) show that the proposed method is reasonably accurate, and its computation time is less than that of pseudoinverses obtained by theqrginv. Hence we conclude that theIMqrginv algorithm is a robust and efficient tool for computing the Moore-Penrose inverse of arbitrary matrices (including singular and rectangular).



---

*Source: 641706-2014-03-12.xml*