# Execute Elementary Row and Column Operations on the Partitioned Matrix to Compute M-P InverseA†

Authors: Xingping Sheng
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We first study the complexity of the algorithm presented in Guo and Huang (2010). After that, a new explicit formula for computational of the Moore-Penrose inverseA† of a singular or rectangular matrix A. This new approach is based on a modified Gauss-Jordan elimination process. The complexity of the new method is analyzed and presented and is found to be less computationally demanding than the one presented in Guo and Huang (2010). In the end, an illustrative example is demonstrated to explain the corresponding improvements of the algorithm.


---

## Body

## 1. Introduction

Throughout this paper we use the following notation. LetCn and Crm×n be the n dimensional complex space and the set of m×n complex matrices with rank r. For a matrix A∈Cm×n, R(A) and N(A) are the range and null space of A; r(A) and A* denote the rank and the the conjugate transpose of A, while A† and ∥A∥ denote the M-P inverse and Frobenius norm, respectively.In 1920, Moore [1] defined a new inverse of a matrix by projection matrices. Moore’s definition of the generalized inverse of an m×n matrix A is equivalent to the existence of an n×m matrix G such that
(1)AG=PR(A),GA=PR(G),
where PR(A) is the orthogonal projector on R(A). Unaware of Moore’s work, In 1955 Penrose [2] showed that there exists a unique matrix X satisfying the four conditions
(2)AXA=A,XAX=X,(AX)*=AX,(XA)*=XA,
where * denotes conjugate transpose. These conditions are equivalent to Moore’s conditions. The unique matrix X that satisfied these conditions was known as the Moore-Penrose inverse (abbreviated M-P) and is denoted by A†. For a subset {i,j,…,k} of the set {1,2,3,4}, the set of n×m matrices satisfied the equations (i),(j),…,(k), from among (2) is denoted by A{i,j,…,k}. These concepts can be found in Ben-Israel and Greville’s famous book [3] or Campell and Meyer’s book [4]. In their famous books [3, 4], the next statement is valid for a rectangular matrix.Lemma 1 (see [3, 4]).
LetA∈Cm×n, X∈Cn×m be the M-P inverse of A if and only if X is a {2}-inverse of A with range R(A*) and null space N(A*).In the latest fifty years, there have been many famous specialists and scholars, who investigated the generalized inverse and published many articles and books. Its perturbation theories were introduced in [5–7], and the algebraical perturbations were in [8, 9]. Some minors, Cramer rulers, and sums of A† can be seen in [10–13]. There are a large number of papers [9, 14–18] using various methods, iterative or not, for computing A†.One handy method of computing the inverse of a nonsingular matrixA is the Gauss-Jordan elimination procedure by executing elementary row operations on the pair (A,I) to transform it into (I,A-1). Moreover Gauss-Jordan elimination can be used to determine whether a matrix is nonsingular or not. However, one cannot directly use this method on a generalized inverse of a rectangular matrix or a square singular matrix A.Recently, Author [19] proposed a Gauss-Jordan elimination algorithm to compute A†, which required 3n3 multiplications and divisions. More recently, Ji [20] improved author’s algorithm [19] and pointed that only 2n3 multiplications and divisions are required. Following these lines, Stanimirovic´ and Petkovic´ in [21] extended the method of [20] to AT,S(2). But these three algorithms need also switching. Guo and Huang [22] executed row and column elementary transformations for computing M-P inverse A† by applying the rank equalities of matrix A. They did not analyze the complexity of their algorithm. In this paper, we first study the total number arithmetic operation of GH-algorithm, then improve it, and present an alternative explicit formula for the M-P inverse of a matrix; the improvements save the total number arithmetic operation. We must point that GH-algorithm and our algorithm do not need to switch blocks of certain matrix in the process of computation.The paper is organized as follows. The computational complexity of GH-Algorithm3 for computing M-P inverse A† is surveyed in the next section. In Section 3, we derive a novel explicit expression for A†, propose a new Gauss-Jordan elimination procedure for A† based on the formula, and study the computational complexity of the new approach and Algorithm 7. In Section 4, an illustrative example is presented to explain the corresponding improvements of the algorithm.


## 2. The Computational Complexity of GH-Algorithm

In [22], Guo and Huang gave a method of elementary transformation for computing M-P inverse A† by applying the rank equalities of matrix A.Lemma 2.
Suppose thatA∈Cm×n, X∈Ck×l, 1≤k≤n, 1≤l≤m. If
(3)r(A*AA*A*(Il0)(Ik0)A*X)=r(A).
Then X=(Ik0)A†(Il0). In particualr when k=n, l=m and X=A†.In [22], the authors also considered an algorithm based on Lemma 2, which was stated as follows.Algorithm 3.
M-P inverse GH-Algorithm is as follows.(1)
Compute partitioned matrixN1=(A*AA*A*A*0).(2)
Make the block matrixA*AA* of N1(1,1) become Ir(A) by applying elementary transformations. Meanwhile, the block matrices A* of N1(1,2) and N1(2,1) are accordingly transformed. A new partitioned matrix N2=(Ir(A)0B1000C100) is obtained.(3)
Make the block matrices ofN2(1,2) and N2(2,1) be zero matrices by applying matrix Ir(A), which yields
(4)N3=(Ir(A)0000000-C1B1)Then A†=C1B1Nevertheless, Guo and Huang [22] did not analyze the complexity of the numerical algorithm. In the following theorem, we will study the total number of arithmetic operations.Theorem 4.
LetA∈Crm×n; the total number of multiplications and divisions required in Algorithm 3 to compute M-P inverse A† is
(5)N(m,n,r)=2mn2+4m-r-12nr+(m-r)nr+mnr.
Moreover, N(m,n,r) is bounded above by 6mn2-(3/2)n3-(1/2)n2.Proof.
It needs2mn2 multiplications to compute A*AA*. r row pivoting steps and column pivoting steps are needed to transform the partitioned matrix N1 into N2 following the r(A*AA*)=r. First row pivoting step involves 2m nonzero columns in (A*AA*A*). Thus, it needs 2m-1 divisions and (2m-1)(n-1) multiplications with a total of (2m-1)n multiplications and divisions. On the second row pivoting step, there is one less column in the first part of the pair. There are 2m-1 nonzero columns to deal with. These pivoting steps require (2m-2)n operations. Following the same idea, the ith (1≤i≤r) pivoting steps require (2m-i)n operations. So it requires (2m-1)n+(2m-2)n+⋯+(2m-r)n=((4m-r-1)/2)nr.
For simplicity, assume thatB1=(IrB12). Following the same line, this requires (m-r)nr multiplications and divisions on the r column pivoting steps.
Then resume elementary row and columns operations on the matrixN2 to transform it into N3, which requires mnr multiplications, which is the count of C1B1.
Therefore, the total number of operations needed for computation ofA† is
(6)N(m,n,r)=2mn2+4m-r-12nr+(m-r)nr+mnr.
Define a functionf(r)=2mn2+((4m-r-1)/2)nr+(m-r)nr+mnr for fixed m and n. Since 0≤r≤min⁡{m,n}=n, we have
(7)f′(r)=4m-r-12n-12nr+(m-r)n-nr+mn=(4m-12-3r)n≥(m-12)n>0
which means that f(r) is monotonically increasing over [0,n] when n≤m. Therefore f(r) takes its maximum value when r=n, which implies
(8)f(r)≤f(n)=6mn2-32n3-12n2.


## 3. Main Results

The Gauss-Jordan row and column elimination procedure for the M-P inverseA† of a matrix by Guo and Huang is based on the partitioned matrix N1=(A*AA*A*A*0). In this section, we will first propose a modified Gauss-Jordan elimination process to compute A†and then summarize an algorithm of this method. Finally, the complexity of the algorithm is analyzed.Theorem 5.
LetA∈Crm×n; there exist an elementary row operation matrix E1=(E11E12)nrn-r and an elementary column operation matrix F1=(F11F12)rm-rm, such that
(9)B=E1A*=(E11E12)A*=(B10),C=A*F1=A*(F11F12)=(C10),
where (B10) and (C10) are the row and column reduced echelon form of A*, respectively.
Further, there exists an elementary row operation matrixE2∈Cr×r such that
(10)E2B1AC1=Ir.
Then
(11)A†=C1E2B1=C1(B1AC1)-1B1=A*F11E2E11A*.Proof.
Forr(A)=r(A*)=r, there exist two elementary row and column operation matrices E1=(E11E12)nrn-r and F1=(F11F12)rm-rm, such that
(12)B=E1A*=(E11A*E12A*)=(B10),C=A*F1=(A*F11A*F12)=(C10).
It is easy to check thatB1∈Crr×n and C1∈Crm×r; then the matrix B1AC1 is nonsingular, which implies that there exists another elementary row operation matrix E2∈Cr×r such that
(13)E2B1AC1=Ir.
The above formula also shows that E2=(B1AC1)-1.
Denote thatX=C1E2B1=A*F11E2E11A*; it is obvious that R(X)⊃R(A*) and N(X)⊂N(A*). If we can prove r(X)=r(A*)=r, then R(X)=R(A*) and N(X)=N(A*).
In fact,C1 is a full column rank matrix and E2 is an invertible matrix, which implies that r(X)=r(B1)=r(A*)=r.
By deducing, we obtain that(14)XAX=C1E2B1AC1E2B1=C1E2B1=X.
This means thatX is a 2-inverse of A with R(X)=R(A*) and N(X)=N(A*). From Lemma 1, we know that X=A†.Remark 6.
The representation ofA† in Theorem 5 is consistent with the one in [3], although we use Gauss-Jordan elimination procedure.
According to the representation introduced in Theorem5, we summarize the following algorithm for computing M-P inverse A†.Algorithm 7.
M-P inverse-Sheng algorithm is as follows.(1)
Input:A∈Crm×n.(2)
Execute elementary row operations on firstn rows of the partitioned matrix M1=(0A*A*0) into M2=(0BA*0), where B=(B10) is a reduced row-echelon matrix.(3)
Perform elementary column operations on firstm columns of the partitioned matrix M2 into M3=(0BC0), where matrix C=(C10) has a reduced column-echelon form.(4)
ComputeD=BAC=(B1AC1000) and form the block matrix
(15)M4=(DBC0)=(B1AC10B1000C100).(5)
Execute the elementary row operations on firstr rows of the partitioned matrix M4 into M5=(Ir0B1~000C100).(6)
Make the block matrices ofM5(1,3) and M5(3,1) be zero matrices by applying elementary row and column transformations, respectively, through matrix Ir, which yields
(16)M6=(Ir0000000-C1B1~).Then A†=C1B1~.
According to Algorithm7, the next theorem will analyze the computational complexity of it.Theorem 8.
The total number of multiplications and divisions required for Algorithm7 to compute M-P inverse A† of a matrix A∈Crm×n is
(17)T(m,n,r)=(4mn-m+n2)r+m-n2r2-r3.
Further, the upper bound ofT(m,n,r) is less than (7/2)mn2-(1/2)n3-((m+n)/2)n when n≤m.Proof.
For a matrixA with rank r, r pivoting steps are needed to make the partitioned matrix M1 into M2. First pivoting step involves m nonzero columns in A*. Thus, it needs m-1 divisions and (m-1)(n-1) multiplications with a total number of (m-1)n multiplications and divisions. On the second pivoting step, there is one less column in the first part of the pair. There are m-1 nonzero columns to deal with. This pivoting step requires (m-2)n operations. Following the same idea, the ith(1≤i≤r) pivoting step requires (m-i)n operations. So these r pivoting steps require
(18)(m-1)n+(m-2)n+⋯+(m-r)n=2m-r-12nr
multiplications and divisions to reach the matrix M2.
Similarly, it needs(n-1)m+(n-2)m+⋯+(n-r)m=((2n-r-1)/2)mr multiplications and divisions to change the matrix M2 into M3.
For simplicity, assume thatB1=(IrB11) and C1=(IrC11), which follows from that B and C are row-echelon and column-echelon reduced matrix, respectively.
(
m
-
r
)
r
n
+
(
n
-
r
)
r
2
=
(
m
n
-
r
2
)
r multiplications are required to form M4=(B1AC10B1000C100) under the above assumption. Since every row of the partitioned matrix (B1AC10B1) has m+1 nonzero elements, each pivoting step needs mr multiplications and divisions. Thus, it requires mr2 multiplications and divisions to obtain the matrix M5.
Then resume elementary row and columns operations on the matrixM5 to transform it into M6. The complexity of this process is mnr multiplications, which is the count to compute C1B1~.
Hence, the total number of complexity of Algorithm7 is
(19)T(m,n,r)=2m-r-12nr+2n-r-12mr+(mn-r2)r+mr2+mnr=(4mn-m+n2)r+m-n2r2-r3.
With fixedm and n, define a function g(r)=(4mn-(m+n)/2)r+((m-n)/2)r2-r3 for 0≤r≤min⁡{m,n}=n. Then we have
(20)g′(r)=(4mn-m+n2)+(m-n)r-3r2=3(mn-r2)+(m-n)r+(mn-m+n2)≥(m-n)+(mn-m+n2)>0
which implies that g(r) is also monotonically increasing over [0,n] when n≤m.
Therefore, whenr=n, g(r) obtains its maximum value, which yields
(21)g(r)≤g(n)=72mn2-12n3-m+n2n.
Furthermore, we give two remarks: one is explaining the computation speed and the other is how to improve the accuracy of Algorithm7.Remark 9.
The algorithm in this paper does not need to switch block of certain matrix in the process computation, unlike the existing algorithm in [19–21]. The higher computational complexity is about 3n3 multiplications and divisions, that is, less than GH-algorithm [22], which requires 4.5n3 multiplications and divisions, when they are applied to the case of r≈n=m for A†.Remark 10.
In order to improve the accuracy of the algorithm, we must select nonzero entries in pivot row and column in each step of the Gauss-Jordan elimination. This improvement is based on the fact that Gauss-Jordan elimination is applied on matrices containing nonnegligible number zero elements.


## 4. Numerical Examples

In this section, we will use a numerical example to demonstrate our results. A handy method is used to computeA† on a low order matrix.Example 1.
Use Algorithm7 to compute the M-P inverse A† of the matrix in [21], where
(22)A=(101112002211).
Solution. Execute elementary row operations on the first four rows of the partitioned matrix M1=(0A*A*0); we have
(23)M1=(0A*A*0)=(000112000022000101000101112000022000101000101000)⟶M2=(000101000011000000000000112000022000101000101000).
Then perform elementary column operations on the first three columns of matrix M2, which yields
(24)M2=(000101000011000000000000112000022000101000101000)⟶M3=(0001010000110000000000001000000200001-100001-10000).
Denote
(25)B=(101011000000),C=(1000101-1201-120).
By computing, we have
(26)D=BAC=(101011000000)(101112002211)(1000101-1201-120)=(700530000000).
We execute elementary row operations on the first two rows of the partitioned matrixM4=(DBC0) again; we have
(27)M4=(7001015300110000000000001000000100001-1200001-120000)⟶M5=(10017017010-521132210000000000001000000100001-1200001-120000).
One then resumes elementary row and column operations onM5, which results in
(28)M5=(10017017010-521132210000000000001000000100001-1200001-120000)⟶M6=(100000010000000000000000000-170-17000521-13-221000-114216-221000-114216-221).
Then we can obtain
(29)A†=(17017-521132211142-162211142-16221).



---

*Source: 596049-2014-03-10.xml*