# Successive Matrix Squaring Algorithm for Computing the Generalized InverseAT,S(2)

Authors: Xiaoji Liu; Yonghui Qin
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We investigate successive matrix squaring (SMS) algorithms for
computing the generalized inverseAT,S(2) of a given matrix A∈Cm×n.


---

## Body

## 1. Introduction

Throughout this paper, the symbolCm×n denotes a set of all m×n complex matrices. Let A∈Cm×n, and the symbols R(A),N(A),ρ(A), and ∥·∥ stand for the range, the null space, the spectrum of matrix A, and the matrix norm, respectively.A matrixB is called a {2}-inverse of matrix A if BAB=B holds. The symbols A†, Ind(A), and AD denote, respectively, the Moore-Penrose inverse, the index, and the Drazin inverse of A, and, obviously, rank(A†)=rank(A) (see [1] for details). Let A∈Crm×n,T⊂Cn, S⊂Cm, and dim(T)=t≤r and dim(S)=m-t, and there exists and unique matrix B∈Cn×m such that
(1.1)BAB=B,R(B)=T,N(B)=S
then B∈Cn×m is called {2}-inverse of A with the prescribed range T and null space S of A, denoted by AT,S(2).In [1], it is well known that the generalized inverse AT,S(2) of a given matrix A∈Cm×n with the prescribed range T and null space S is very important in applications of many mathematics branches such as stable approximations of ill-posed problems, linear and nonlinear problems involving rank-deficient generalized, and the applications to statistics [2]. In particular, the generalized inverse AT,S(2) plays an important role for the iterative methods for solving nonlinear equations [1, 2].In recent years, successive matrix squaring algorithms are investigated for computing the generalized inverse of a given matrixA∈Cm×n in [3–7]. In [3], the authors exhibit a deterministic iterative algorithm for linear system solution and matrix inversion based on a repeated matrix squaring scheme. Wei derives a successive matrix squaring (SMS) algorithm to approximate the Drazin inverse in [4]. Wei et al. in [5] derive a successive matrix squaring (SMS) algorithm to approximate the weighted generalized inverse AM,N†, which can be expressed in the form of successive squaring of a composite matrix T. Stanimirović and Cvetković-Ilić derive a successive matrix squaring (SMS) algorithm to approximate an outer generalized inverse with prescribed range and null space of a given matrix A∈Crm×n in [6]. In [7], authors introduce a new algorithm based on the successive matrix squaring (SMS) method and this algorithm uses the strategy of ϵ-displacement rank in order to find various outer inverses with prescribed ranges and null spaces of a square Toeplitz matrix.In this paper, based on [3–5], we investigate successive matrix squaring algorithms for computing the generalized inverse AT,S(2) of a matrix A in Section 2 and also give a numerical example for illustrating our results in Section 3.The following given lemma suggests that the generalized inverseAT,S(2) is unique.Lemma 1.1 (see [1, Theorem  2.14]).
LetA∈Cm×n with rank r, let T be a subspace of ℂn of dimension s≤r, and let S be a subspace of ℂm of dimension m-s. Then, A has a {2}-inverse X such that ℛ(X)=T and 𝒩(X)=S if and only if
(1.2)AT⊕S=ℂm
in which case X is unique.The following nations are stated in Banach space but they are true in the finite dimension space. Throughout this paper, letH,K denote the Banach space and let ℬ(H,K) stand for the set of all bounded linear operators from H to K, in particular ℬ(H,H)=ℬ(H).In the following, we state two lemmas which are given for Banach space but it can be used also for the finite dimension space.Lemma 1.2 (see [8, Section  4]).
LetA∈ℬ(H,K) and T and S, respectively, closed subspaces of H and K. Then the following statements are equivalent:(i)
Ahas a {2}-inverse B∈K,H such that R(B)=T and N(B)=S,(ii)
T is a complemented subspace of H, A|T:T→A(T) is invertible and A(T)⊕S=K.Lemma 1.3 (see [9, Section  3]).
Suppose that the conditions of Lemma1.2 are satisfied. If we take T1=N(AT,S(2)A), then H=T⊕T1 holds and A has the following matrix form:
(1.3)A=[A100A2]:[TT1]→[A(T)S],
where A1 is invertible. Moreover, AT,S(2) has the matrix following form:
(1.4)AT,S(2)=[A1-1000]:[A(T)S]→[TT1].From (1.5), we obtain the following projections (see [9]):
(1.5)PA(T),S=AAT,S(2)=[I000]:[A(T)S]→[A(T)S],PT,T1=AT,S(2)A=[I000]:[TT1]→[TT1].


## 2. Main Result

In this section, we consider successive matrix squaring (SMS) algorithms for computing the generalized inverseAT,S(2).LetA∈Cm×n and the sequence {Xn} in Cn×m, and we can define the iterative form as follows ([10, Theorem  2.2] for computing the generalized inverse AT,S(2) in the infinite space case):
(2.1)Rk=PA(T),S-PA(T),SAXk,Xk+1=X0Rk+Xk,k=0,1,2,….
From [10], the authors have proved that the iteration (2.1) converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1, where T⊂Cn and PA(T),S=AAT,S(2) (for the proof see [11] and [10, Theorem  2.1] when p=2).In the following, we give the algorithm for computing the generalized inverseAT,S(2) of a matrix A∈Cm×n.LetP=R0=PA(T),S-PA(T),SAX0 and Q=X0. It is not difficult to see that the above fact can be written as follows:
(2.2)M=[R00X0I]=[P0QI].
From (2.2) and letting Xk=Q∑i=0kPi, we have
(2.3)Mk=[Pk0Q∑i=0k-1‍PiI]=[Pk0Xk-1I].By (2.3), we prove that the iterative (2.1) Xk is equal to the right upper block in the matrix Mk. Note that we defined the new iterative form {Mk} as follows:
(2.4)M0=M,Mk+1=Mk2,k=0,1,2,….From the new iterative form (2.4), we arrive at
(2.5)Mk=M2k=[P2k0Q∑i=02k-1‍PiI]=[P2k0X2k-1I].Assume thatX2k-1=X^k, and by (2.5), we have
(2.6)Mk=[P2k0X^kI]=[P2k0X2k-1I]=[P2k0Q∑i=02k-1‍PiI].By (2.4)–(2.6), we have Algorithm 1.Algorithm 1:SMS algorithm for computing the generalized inverse AT,S(2).
Input: Input the initial value matricesA,X0,PA(T),S and the accurate value ϵ;
Output: The algorithm export the matrix, that isX≈AT,S(2);
Begin: Assignment the matrixQ by the initial value matrix X0, that is Q⇐X0;
Assigned the matrixP by PA(T),S-PA(T),SAQ, that is P⇐PA(T),S-PA(T),SAQ;
Computed matrixX1, that is X1⇐I+P;
Computed the error betweenX1 and X0, that is e⇐∥X1-X0∥;
Judged that whethere is lower than ϵ or not,
that is whilee<ϵ, do P⇐P·P;
Defined the loop function:Xk+1⇐Xk+P;
Computed the error betweenXk and Xk+1, that is e⇐∥Xk+1-Xk∥;
Finished the loop function
Thek+1 matrix Xk+1 multiplied by Q and assigned to X, that is X⇐QXk+1;
End the algorithm.From (2.4)–(2.6) and Algorithm 1, we obtain the following result.Theorem 2.1.
LetA∈Cm×n, and the sequence {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. In this case
(2.7)∥AT,S(2)-X^k∥≤q2k+1(1-q)-1∥X0∥,
where q=∥R0∥ and
(2.8)T⊂Cn,PA(T),S=AAT,S(2).Proof.
From the proof in [11] and [10, Theorem  2.1] when p=2 and according to (2.4), (2.5) and (2.6), we easily finish the proof of the former of the theorem. In the following, we only prove the last section, that is, prove that the inequality (2.7) holds.
By applying (2.5) and (2.6), we obtain
(2.9)X^k=X2k-1=∑i=02k-1‍PiQ.
By the iteration (2.4) and (2.9), we arrive at
(2.10)∥AT,S(2)-X^k∥=∥X0(I-R0)-1-X0∑i=02k-1‍R0i∥=∥X0∑i=0∞‍R0i-X0∑i=02k-1‍R0i∥=∥X0∑i=2k∞‍R0i∥=∥X0R02k∑i=0∞‍R0i∥≤∥R02k∥∑i=0∞‍∥R0i∥∥X0∥≤q2k+1(1-q)-1∥X0∥.The following corollary given the result is the same as theorem in [6, Theorem  2.3]. It also presents an explicit representation of the the generalized inverse AT,S(2) and the sequence (2.4) converges to a {2}-inverse of a given matrix A by its full-rank decomposition.Corollary 2.2.
LetA∈Crm×n,A=FG be full rank decomposition, and the sequence {X^k} converges to the {2}-inverse X=F(GAF)-1G if and only if ρ(R0)<1. In this case
(2.11)∥X-X^k∥≤q2k+1(1-q)-1∥X0∥,
where q=∥R0∥ and
(2.12)F∈Csm×s,G∈Css×n,PR(AX),N(AX)=AX.Proof.
From Theorem2.5 and by [6, Theorem  2.3], we have the result.In the following, we consider the improvement of the iterative form (2.1) (see [11] for computing the Moore-Penrose inverse and the Drazin inverse of the matrix case and [10, Theorem  2.2] for computing the generalized inverse AT,S(2) in the infinite space case):
(2.13)Rk=PA(T),S-PA(T),SAXk,Xk+1=Xk(I+Rk+⋯+Rkp-1),p≥2,k=0,1,2,….
Let M be a m×m block matrix and
(2.14)M=[Pm-10⋯0Pm-20⋯0*0P0⋯0QQ⋯I],
then
(2.15)M2=[P2m-10⋯0P2m-20⋯0*0Pm0⋯0Q∑i=0m-1‍PiQ⋯I].By induction ifMk-1 has the following form:
(2.16)Mk-1=[P(k-1)m-10⋯0P(k-1)m-20⋯0*0P(k-2)m0⋯0Q∑i=0(k-2)m-1‍PiQ⋯I],
then
(2.17)Mk=[Pkm-10⋯0Pkm-20⋯0*0P(k-1)m0⋯0Q∑i=0(k-1)m-1‍PiQ⋯I].Similarly to the iterative form (2.4), we also define the new iterative scheme {Mk}(2.18)M0=M,Mk+1=Mkp,k=0,1,2,….
Note that from (2.18)
(2.19)Mk=Mpk=[Ppkm-10⋯0Ppkm-20⋯0*0P(pk-1)m0⋯0Q∑i=0(pk-1)m-1‍PiQ⋯I]=[Ppkm-10⋯0Ppkm-20⋯0*0P(pk-1)m0⋯0X(pk-1)m-1Q⋯I].LetX(pk-1)m-1=X^k, and by (2.18), and (2.19), we arrive at
(2.20)Mk=[*0X(pk-1)m-1*]=[*0X^k*].From (2.14) to (2.20), we find that if one wants to compute the generalized inverse AT,S(2) then we only compute the element (m,1) of the matrix M2k. Similarly to Algorithm 1, we also obtain Algorithm 2.Algorithm 2:SMS algorithm for computing the generalized inverse AT,S(2).
Input: Input the matricesA,X0,PA(T),S and the accurate value ϵ;
Output: The algorithm export the matrix:X≈AT,S(2);
Begin: Assignment the matrixQ by the initial value matrix X0, that is Q⇐X0;
Assigned the matrixP1 by PA(T),S-PA(T),SAQ, that is P1⇐PA(T),S-PA(T),SAQ;
Computed the product ofP1 and P1, and assigned its value to P2. that is P2⇐P1·P1;
Similarly, we repeatedly do the computation for the productPi and P1 as well as above
the computation, wherei=2,…,m-2.
Computed the product of the matrixPm-1 and P1, and assigned its value to Pm as well
as above computations, that isPm⇐Pm-1·P1;
Assigned the matrixX1 by the sum of the matrices Pi, where i=0,1,2,…,m and
P
0
=
I. that is X1⇐I+P1+···+Pm-1+Pm;
Take the norm of∥X1-X0∥ and assigned its value to e. that is e⇐∥X1-X0∥;
whilee<ϵ do;
We need the iteration not to exceed 500 times. that isn=500; (In fact (pk-1)m-1=500)
Do 500 step repeatedly computations in the following.
that is Fori= 1 : n
Computed the product of the given matrixPm=Pm and the iteration matrix Pi, and
assigned its value to the new matrixPi. that is Pi⇐Pm·Pi;
From the iterationPi⇐Pm·Pi, we obtain the new matrix Pi and add its value to Xi,
and assigned the sum ofPi and Xi to the matrix Xi+1..that is Xi+1⇐Xi+Pi; After these,
return the the stepPi=Pm·Pi.
Finished theFor loop function that is end
Computed the error betweenXk and Xk+1, that is e⇐∥Xk+1-Xk∥;
Finished theWhile loop function. that is end
Thek+1 matrix Xk+1 multiplied by Q and assigned to X, that is X⇐QXk+1;
End the algorithm.Analogous to Theorem2.5 by Algorithm 2 and sequence (2.18), we also have the following theorem.Theorem 2.3.
LetA∈Cm×n, and the sequence {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. In this case
(2.21)∥AT,S(2)-X^k∥≤q(pk-1)m+1(1-q)-1∥X0∥,
where q=∥X0∥ and
(2.22)T⊂Cn,PA(T),S=AAT,S(2).Proof.
Similarly the proof in [10, Theorem  2.1], we can prove the former of this theorem. Analogous to the proof of Theorem 2.5, we finish the proof of the theorem.In the following, we extend the sequence (2.4) to
(2.23)M0=M,Mk+1=Mkt,k=0,1,2,…,foranyt≥2.
By (2.26) and by induction, we have
(2.24)Mk=Mtk=[Ptk0Q∑i=0tk-1‍PiI].Assume thatXtk-1=X^k, we easily have
(2.25)Mk=[Ptk0X^kI]=[Ptk0Xtk-1I]=[Ptk0Q∑i=0tk-1‍PiI].Similarly, from (2.23) and (2.25), we obtain the following result.Theorem 2.4.
LetA∈Cm×n, and the sequence {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. In this case
(2.26)∥AT,S(2)-X^k∥≤qtk+1(1-q)-1∥X0∥,
where q=∥R0∥ and
(2.27)T⊂Cn,PA(T),S=AAT,S(2).Proof.
From (2.25) and only using t instead of 2 in Theorem 2.1, we easily have that {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. Similarly to the formula (2.29), we obtain that
(2.28)∥AT,S(2)-X^k∥≤qtk+1(1-q)-1∥X0∥,
where q, T, and PA(T),S are the same as Theorem 2.5.In the following, we consider the dually iterative form.LetA∈Cm×n and the sequence {Xn} in Cn×m, and we can define the iterative form as follows (see [11] and [10, Theorem  2.3]):
(2.29)Rk=PT,T1-AXkPT,T1,Xk+1=RkX0+Xk,k=0,1,2,….LetP=R0=PT,T1-X0APT,T1 and Q=X0. It is not difficult to see that the above fact can be written as follows:
(2.30)M=[R00X0I]=[P0QI].From iterative forms (2.26) and (2.29), we have the following theorem.Theorem 2.5.
LetA∈Cm×n, and the sequence {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. In this case
(2.31)∥AT,S(2)-X^k∥≤q2k+1(1-q)-1∥X0∥,
where q=∥R0∥ and
(2.32)T⊂Cn,PA(T),S=AAT,S(2).Similarly to Corollary2.2, we have the result as follows.Corollary 2.6.
LetA∈Crm×n, A=FG full rank decomposition, and the sequence {X^k} converges to the {2}-inverse X=F(GAF)-1G if and only if ρ(R0)<1. In this case
(2.33)∥X-X^k∥≤q2k+1(1-q)-1∥X0∥,
where q=∥R0∥ and
(2.34)F∈Csm×s,G∈Css×n,PR(XA),N(XA)=XA.In the following, we consider the improvement of the iterative form (2.29) (see [11] for computing the Moore-Penrose inverse and the Drazin inverse of the matrix case and [10, Theorem  2.3] for computing the generalized inverse AT,S(2) in the infinite space case):
(2.35)Rk=PT,T1-AXkPT,T1,Xk+1=(I+Rk+⋯+Rkp-1)Xk,p≥2,k=0,1,2,….It is similar to (2.14), and we have
(2.36)M=[PmPm-1⋯PQ00⋯0Q*000⋯0Q00⋯0I].Analogous to Theorem2.5 by Algorithm 2 and from (2.36), we obtain the theorem in the following.Theorem 2.7.
LetA∈Cm×n, and the sequence {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. In this case
(2.37)∥AT,S(2)-X^k∥≤q(pk-1)m+1(1-q)-1∥X0∥,
where q=∥X0∥ and
(2.38)T⊂Cn,PT,T1=AT,S(2)A.Dually, we give the SMS algorithm for computing the generalized inverseAT,S(2) which are analogous to the iterative form (2.23) as follows and omit their proofs:
(2.39)Mk=[PtkX^k0I]=[PtkXtk-10I]=[PtkQ∑i=0tk-1‍Pi0I].Similarly Theorem2.4, from (2.35) and (2.39), we obtain the following result.Theorem 2.8.
LetA∈Cm×n, and the sequence {X^k} converges to the generalized inverse AT,S(2) if and only if R(X0)⊂T,ρ(R0)<1. In this case
(2.40)∥AT,S(2)-X^k∥≤qtk+1(1-q)-1∥X0∥,
where q=∥R0∥ and
(2.41)T⊂Cn,PT,T1=AT,S(2)A.


## 3. Example

Here is an example to verify the effectiveness of the SMS method.Example 3.1.
Let(3.1)A=[210200].
Let T∈C2; e=(0;0;1)T∈C3,S=span{e}.
Take(3.2)X0=[0.40000.40].
By (2.2), we have
(3.3)R0=[0.2-0.4000.20000].
From [10, 12], we easily have the generalized inverse AT,S(2) in
(3.4)AT,S(2)=[0.5-0.25000.50].Then, from Algorithm1, we obtain
(3.5)X1=[0.4800-0.1600000.48000],X2=[0.5600-0.3200000.56000].But by the iteration (2.1), we get
(3.6)X1=[0.4800-0.1600000.48000],X2=[0.4960-0.2240000.49600],X3=[0.4992-0.2432000.49920],X4=[0.4998-0.2483000.49980],X5=[0.5000-0.2496000.50000],X6=[0.5000-0.2499000.50000].From the data in (3.5) and (3.6), we obtain Table 1.Table 1
Method
Iteration (2.1)
Algorithm1
Steps
5
2From the above in (3.5), (3.6), and Table 1, we know that we only need two steps by Algorithm 1, but five steps by using iterative form (2.1).



---

*Source: 262034-2012-12-27.xml*