Successive matrix squaring algorithm for ... - Semantic Scholar
Successive matrix squaring algorithm for computing outer inverses Predrag S. Stanimirovi´ c∗ and Dragana S. Cvetkovi´ c-Ili´ c University of Niˇs, Department of Mathematics, Faculty of Science, † P.O. Box 224, Viˇsegradska 33, 18000 Niˇs, Serbia E-mail: [email protected], [email protected]
Abstract In this paper we derive a successive matrix squaring (SMS) algorithm to approximate an outer generalized inverse with prescribed range and null space of a given matrix A ∈ Cm×n . We generalize the results from r the papers [3], [16], [18], and obtain an algorithm for computing various classes of outer generalized inverses of A. Instead of particular matrices used in these articles, we use an appropriate matrix R ∈ Cn×m , s ≤ r. s Numerical examples are presented. AMS Subj. Class.: 15A09. Key words: Generalized inverse; outer inverse; SMS algorithm; full rank factorization; matrix rank.
1
Intoduction and preliminaries
Let Cm×n and Cm×n denote the set of all complex m × n matrices and all r complex m × n matrices of rank r, respectively. In denotes the unit matrix of order n. By A∗ , R(A), rank(A) and N (A) we denote the conjugate transpose, the range, the rank and the null space of A ∈ Cm×n . By Re z and Im z we denote a real and imaginary part of a complex number z, respectively. For A ∈ Cm×n , the set of inner and outer generalized inverses are defined by the following, respectively: A{1} = {X ∈ Cn×m | AXA = A},
A{2} = {X ∈ Cn×m | XAX = X}.
The set of all outer inverses with prescribed rank s is denoted by A{2}s , 0 ≤ s ≤ r = rank(A). The symbols A− or A(1) stand for an arbitrary generalized ∗ Corresponding
author by Grant No. 144003 of the Ministry of Science, Technology and Development, Republic of Serbia. † Supported
1
2
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
inner inverse of A and by A(2) we denote an arbitrary generalized outer inverse of A. Also, the matrix X which satisfies AXA = A
and
XAX = X
is called the reflexive g-inverse of A and it is denoted by A(1,2) . The set of all reflexive g-inverses is denoted by A{1, 2}. Subsequently, the sets of {1, 2, 3} and {1, 2, 4} inverses of A are defined by A{1, 2, 3} = A{1, 2} ∩ {X| (AX)∗ = AX}, A{1, 2, 4} = A{1, 2} ∩ {X| (XA)∗ = XA}. By A† we denote the Moore-Penrose inverse of A, i.e. the unique matrix A† satisfying AA† A = A, A† AA† = A† , (AA† )∗ = AA† , (A† A)∗ = A† A. For A ∈ Cn×n the smallest nonnegative integer k such that rank(Ak+1 ) = rank(Ak ) is called the index of A and denoted by ind(A). If A ∈ Cn×n is a square matrix with ind(A) = k, then the matrix X ∈ Cn×n which satisfies the following conditions Ak XA = Ak , XAX = X, AX = XA, is called the Drazin inverse of A and it is denoted by AD . When ind(A) = 1, Drazin inverse AD is called the group inverse and it is denoted by A# . Suppose that M and N are Hermite positive definite matrices of the order m and n, respectively. Then there exists the unique matrix X ∈ Cn×m such that AXA = A, XAX = X, (M AX)∗ = M AX, (N XA)∗ = N XA. The matrix X is called the weighted Moore-Penrose inverse of A, and denoted by X = A†M,N . In particular, if M = Im and N = In , then A†M,N = A† . If A ∈ Cn×m and W ∈ Cm×n , then the unique solution X ∈ Cn×m of the equations (AW )k+1 XW = (AW )k ,
XW AW X = X,
AW X = XW A,
(1.1)
where k = ind(AW ), is called the W-weighted Drazin inverse of A and it is denoted by AD,W . If A ∈ Cm×n , T is a subspace of Cn of dimension t ≤ r and S is a subspace r m of C of dimension m − t, then A has a {2} inverse X such that R(X) = V and N (X) = U if and only if AV ⊕ U = Cm , (2)
in which case X is unique and we denote it by AV,U .
Successive matrix squaring algorithm for computing outer inverses
3
It is well-known that for A ∈ Cm×n , the Moore-Penrose A† , the weighted Moore-Penrose inverse A†M,N and the weighted Drazin inverse AD,W can be represented by: (2)
(i) A† = AR(A∗ ),N (A∗ ) , (2)
(ii) A†M,N = AR(A] ),N (A] ) , where A] = N −1 A∗ M, (2)
(iii) AD,W = (W AW )R(A(W A)k ),N (A(W A)k ) , where W ∈ Cn×n , k = ind(W A). Also, for A ∈ Cn×n , the Drazin inverse AD can be represented by: (2)
AD = AR(Ak ),N (Ak ) , where ind(A) = k. The following representations of {2, 3}, {2, 4}-inverses with prescribed rank s are restated from [11]: Proposition 1.1 Let A ∈ Cm×n and 0 < s < r be a chosen integer. Then the r following is valid: © ª (a) A{2, 4}s = (ZA)† Z | Z ∈ Cs×m , ZA ∈ Cs×n ; s © ª † n×s m×s (b) A{2, 3}s = Y (AY ) | Y ∈ C , AY ∈ Cs . General representations for various classes of generalized inverses can be found in [4, 8, 10, 12]. Some of these representations are restated here for the sake of completeness. Proposition 1.2 Let A ∈ Crm×n be an arbitrary matrix and A = P Q is a full-rank factorization of A. There are the following general representations for some classes of generalized inverses: A{2}s = {F (GAF )−1 G | F ∈ Cn×s, G ∈ Cs×m, rank(GAF ) = s}; r [ A{2} = A{2}s ; s=0
A{1, 2} = {F (GAF )−1 G | F ∈ Cn×r, G ∈ Cr×m, rank(GAF ) = r} = A{2}r ; A{1, 2, 3} = {F (P ∗ AF )−1 P ∗ | F ∈ Cn×r , rank(P ∗ AF ) = r}; A{1, 2, 4} = {Q∗ (GAQ∗ )−1 G | G ∈ Cr×m, rank(GAQ∗ ) = r}; A† = Q∗ (P ∗ AQ∗ )−1 P ∗ ; AD = PAl (QAl A PAl )−1 QAl , Al = PAl QAl , l ≥ ind(A). For other important properties of generalized inverses see [1, 2, 6, 13]. We will use the following well-known result: Lemma 1.1 [7] Let M ∈ Cn×n and ε > 0 be given. There is at least one matrix norm k · k such that ρ(M ) ≤ kM k ≤ ρ(M ) + ², (1.2) where ρ(M ) denotes the spectral radius of M .
4
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
The basic motivation for this paper were the paper of L. Chen et al. [3] and the papers of Y. Wei [16] and Y. Wei et al. [18]. In the paper of L. Chen et al. [3] a successive matrix squaring (SMS) algorithm which approximates the Moore-Penrose inverse of a given matrix A were considering, while in the papers of Y. Wei [16] and Y. Wei et al. [18] the authors considered two variants of SMS algorithm which approximate the Drazin inverse and the weighted MoorePenrose inverse of A, respectively. In this paper we derive a SMS algorithm to approximate an outer generalized inverse with prescribed range and null space of a given matrix A ∈ Cm×n , r which is based on successive squaring of an appropriate block matrix T . We generalized the results from the papers [3], [16] and [18]. Instead of matrices A∗ , Ak , k ≥ ind(A) and A] , used to compose the matrix T in these articles, , 0 ≤ s ≤ r. Furthermore, we use an universal appropriate matrix R ∈ Cn×m s we give an explicit form of the outer generalized inverse of A corresponding to chosen R. An explicit approximation of rank(X) is also given in the same iterative scheme. We also estimate the number of iterative steps necessary to obtain a prescribed precision. Numerical examples in the last section show that numerical approximation of rank(X) is important and that the number of SMS iterations is less than in the classical Euler-Knopp iterations.
2
Results
Let A ∈ Cm×n and R ∈ Cn×m , 0 ≤ s ≤ r be given. The general iterative r s scheme used in this paper is given by X1 = Q, Xk+1 = P Xk + Q, k ∈ N,
(2.1)
where P = I − βRA, Q = βR and β is a relaxation parameter. The iterative process (2.1) can be accelerated by means of (m + n) × (m + n) composite matrix T , partitioned in the block form · ¸ P Q T = . (2.2) 0 I The improvement of (2.1) can be done by computing the matrix power · Pk−1 i ¸ Pk i=0 P Q . Tk = 0 I Pk−1 It is not difficult to see that the iterative scheme (2.1) gives Xk = i=0 P i Q. k Hence, the matrix Xk is equal to the right upper block in T k . In turn, T 2 can be computed by k repeated squaring, that is T0 Ti+1
= =
T Ti2 , i = 0, 1, . . . , k − 1.
(2.3)
Successive matrix squaring algorithm for computing outer inverses It is clear that
" Tk = T
2k
P2 0
=
k
P2k −1 i=0
5
# P iQ
I
.
(2.4)
Now, we will state an auxiliary result: Lemma 2.1 If P ∈ Cn×n and S ∈ Cn×n are such that P = P 2 and P S = SP then ρ(P S) ≤ ρ(S). Proof. Suppose that λ is an eigenvalue of P S and x is a corresponding vector such that P Sx = λx. Then, λP x = λx which implies that λ = 0 or P x = x. If P x = x, then Sx = SP x = P Sx = λx, i.e. λ is an eigenvalue of S. Hence, the set of eigenvalues of P S is the subset of the set of eigenvalues of S union {0}, which implies that ρ(P S) ≤ ρ(S). The following theorem represents the analogue of Theorem 3.1 [15] which deals with V -norm, where V-norm is defined as kBkV = kBV k2 , for any B ∈ Cn×m and V is invertible such that V −1 AGV is the Jordan form of AG. The method used in our proof is different then the one used in Theorem 3.1 [15]. and R ∈ Csn×m , 0 ≤ s ≤ r be given such that Theorem 2.1 Let A ∈ Cm×n r AR(R) ⊕ N (R) = Cm .
(2.5)
The sequence of approximations X2k =
k 2X −1
(I − βRA)i βR
(2.6)
i=0
determined by the SMS algorithm (2.3) converges in the matrix norm k · k to the (2) outer inverse X = AR(R),N (R) of A if β is a fixed real number such that max |1 − βλi | < 1,
1≤i≤t
(2.7)
where rank(RA) = t, λi , i = 1, . . . , t are eigenvalues of RA and k · k satisfies condition (1.2) from Lemma 1.1 for M = I − βAR. In the case of convergence we have the following error estimate k kX − X2k k ≤ max |1 − βλi |2 + O(ε), k ≥ 0. 1≤i≤t kXk
(2.8)
Proof. From the condition (2.5), it follows the existence of the outer inverse (2) X = AR(R),N (R) of A. It is evident that XAR = R and R = RAX. Since R(X2k ) ⊆ R(R), we have that XAX2k = X2k .
6
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
So, we obtain kX − X2k k = kXAX − XAX2k k = kX(AX − AX2k )k ≤ kXk · kAX − AX2k k.
(2.9)
Now, let us prove that AX − AXk = (AX − AX1 )k .
(2.10)
We will use the mathematical induction. For k = 1, the equality (2.10) is true. Suppose that (2.10) holds for k = p−1 and prove that it holds for k = p. Using P = I − QA, we have AX − AXp
= = =
AX − A(P Xp−1 + Q) AX − A(I − QA)Xp−1 − AQ AX − AXp−1 + AQAXp−1 − AQ.
Since R = RAX implies Q = QAX, we get AX − AXp
= AXAX − AXAXp−1 + AQAXp−1 − AQAX = (AX − AQ)(AX − AXp−1 ) = (AX − AQ)p = (AX − AX1 )p .
Now, (2.9) and (2.10) imply k
kX − X2k k ≤ kXk · kAX − AX2k k ≤ kXk · kAX − AX1 k2 . If (2.7) is satisfied, we can choose ² such that max1≤i≤t |1 − βλi | + ε < 1. For such ² and M = I − βAR, there exists a matrix norm k · k satisfies (1.2). Now, µ ¶2k k kX − X2k k ≤ kAX − AX1 k2 ≤ max |1 − βλi | + ε . 1≤i≤t kXk The last inequality follows if we apply Lema 2.1 for P = AX and S = I −AX1 = I − βAR and by Lemma 1.1. Remark that a similar problem as in Theorem 3.1 [15] was discussed in [17] but in a more general case: If {Sn (x)} is a family of continuous real valued functions on open set Ω such that σ((RA)|T ) ⊂ Ω ⊂ (0, ∞) and limn→∞ Sn (x) = 1/x uniformly on σ((RA)|T ), then (2)
AT,S = lim Sn ((RA)|T )R n→∞
and (2)
kSn ((RA)|T )R − AT,S kV (2)
kAT,S kV
≤
max
x∈σ((RA)|T )
|Sn (x)x − 1| + O(ε), ² > 0.
The generalization to the bounded linear operator on Banach space of the mentioned result proved in [17] is given in [5].
Successive matrix squaring algorithm for computing outer inverses
7
Corollary 2.1 Under the assumption and the conditions of Theorem 2.1, (i) In the case m = n, R = Al , l ≥ ind(A), we have AD = lim X2k , k→∞
(ii) In the case R = A∗ , A† = lim X2k , k→∞
(iii) In the case R = A] = N −1 A∗ M , A†M,N = lim X2k , k→∞
(iv) In the case R = A(W A)k , AD,W = lim X2k . k→∞
Corollary 2.2 In the case rank(R) = rank(A), under the conditions of Theorem 2.1, we have X ∈ A{1, 2}. Proof. Using that rank(X) = rank(A), by Corollary 1 (pg. 19, [1]) it follows that X ∈ A{1, 2}. Theorem 2.2 Let A ∈ Rrm×n and R ∈ Rsn×m , 0 ≤ s ≤ r be real matrices and assume that β satisfies condition (2.7). Then the sequence of real numbers defined by k
k
r2k = n − Tr((I − βRA)2 ) = n − Tr(P 2 ), k ≥ 0
(2.11)
converges to rank(RA). Proof. It is sufficient to use k
Tr((I − βRA)2 ) = n − t +
t X k (1 − βλi )2 ,
(2.12)
i=1
where t = rank(RA) ≤ rank(R). Therefore, using (2.12), we get r 2k = t −
t X
k
(1 − βλi )2 .
(2.13)
i=1
Now, the proof can be completed using (2.7). Remark 2.1 Taking R = A∗ in Theorem 2.1 and Theorem 2.2 we get the results of Theorem 2 from [3]. Also, taking R = Al , l ≥ ind(A) and R = A] in Theorem 2.1, respectively, we get the results from [16] and [18]. The next theorem presents an explicit representation of the outer inverse of A which is a limit value of the sequence of approximations (2.1), for given A and R.
8
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
Theorem 2.3 Assume that A∈Cm×n . Consider an arbitrary matrix R∈Cn×m , r s s×m 0 ≤ s ≤ r, and its full rank factorization R = F G, F ∈ Cn×s , G ∈ C , such s s that GAF is invertible. The sequence of approximations X2k =
k 2X −1
(I − βRA)i βR
(2.14)
i=0
determined by the SMS algorithm (2.3) converges in the matrix norm to the outer inverse X = F (GAF )−1 G of A, if β is a fixed real number satisfying max |1−βλi | < 1,
(2.15)
1≤i≤s
where λi , i = 1, . . . , s are eigenvalues for F GA. In the case of convergence we have the error estimate k kX − X2k k ≤ max |1 − βλi |2 + O(ε), k ≥ 0, 1≤i≤s kXk
(2.16)
where the matrix norm k · k satisfies condition (1.2), for M = I − βAR. Proof. Since the matrices F and G satisfy rank(GAF ) = rank(F ) = rank(G), according to Theorem 1.3.8 (pg. 33, [13]) and [14], we conclude that X = F (GAF )−1 G is a {2} inverse of A having the range R(X) = R(F ) = V and null space N (X) = N (G) = U , i.e. (2)
(2)
X = AV,U = AR(F ),
N (G) .
Also, conditions conditions AV ∩ U = {0} and AV ⊕ U = Cm for the existence of X are ensured (Theorem 1.3.8, pg. 33, [13]). On the other hand, since rank(R) = s = rank(F G) = rank(F ) = rank(G) (2)
we obtain R(F ) = R(R), N (G) = N (R) which implies X = AR(R),N (R) . In view of (2.15) immediately follows that the conditions of Theorem 2.1 are satisfied. Therefore, we may conclude that the sequence of approximations (2.14) determined by the SMS algorithm (2.3) converges to X = F (GAF )−1 G, in the matrix norm satisfying (1.2) for M = I − βAR. Corollary 2.3 Assume that A ∈ Rm×n and R ∈ Rn×m , 0 ≤ s ≤ r, are two r s arbitrary chosen real matrices. Let R = F G, F ∈ Cn×s , G ∈ Cs×m be a s s full rank factorization of R, such that GAF is invertible. The sequence of real numbers defined by k
k
r2k = n − Tr((I − βRA)2 ) = n − Tr(P 2 ), k ≥ 0 converges to rank(X).
(2.17)
Successive matrix squaring algorithm for computing outer inverses
9
Proof. According to Theorem 2.2, r2k → rank(RA). Since rank(RA) = rank(GAF ) = s = rank(R) = rank(X), we complete the proof. Corollary 2.4 Consider A ∈ Cm×n and arbitrary chosen matrix R ∈ Cn×m , r s n×s s ≤ r. Assume that R = F G, F ∈ Cs , G ∈ Css×m is a full-rank factorization of R. Also, suppose that condition (2.15) of Theorem 2.3 are valid and that A = P Q is a full rank factorization of A. If the sequences of approximations X2k is defined by (2.14), the following hold: ½ ¾ s×m a) lim X2k | F ∈ Cn×s , G ∈ C , rank(GAF ) = s, 0 ≤ s ≤ r ⊂ A{2}; s s k→∞
½ b)
¾ n×r
lim X2k | F ∈ C
k→∞
,G ∈ C
r×m
, rank(GAF ) = r
½ c)
¾ n×r
lim X2k | F ∈ C
k→∞
∗
, G = P , rank(GAF ) = r
½ d)
⊂ A{1, 2};
⊂ A{1, 2, 3}; ¾
∗
r×m
lim X2k | F = Q , G ∈ C
k→∞
, rank(GAF ) = r
⊂ A{1, 2, 4};
satisfy rank(GAF ) = s, Proof. a) If the matrices F ∈ Csn×s and G ∈ Cs×m s applying Theorem 2.3, we conclude that X = F (GAF )−1 G is the limit value of the sequence of approximations (2.14). Using general representation of outer inverses form Proposition 1.2, we immediately obtain X ∈ A{2}. Parts b)-d) can be proved using Theorem 2.3 and general representations of pseudoinverses from Proposition 1.2. and integer s such that 0 ≤ s ≤ r. Corollary 2.5 Consider matrix A ∈ Cm×n r Then the following is valid: ( Ã k−1 ! ) 2P i ∗ ∗ n×s a) Y · lim (I − β(AY ) AY ) (AY ) |Y ∈C , rank(AY ) = s k→∞
i=0
⊂ A{2, 3}s ; ( Ã k−1 ! ) 2P i b) lim (I − β(ZA)∗ ZA) (ZA)∗ · Z | Z ∈ Cs×n , rank(ZA) = s k→∞
i=0
⊂ A{2, 4}s . Proof. Follows from Corollary 2.1, part (ii) and the general representations of the sets A{2, 3}s and A{2, 4}s given in Proposition 1.1. In the following theorem we investigate the limits of the relaxation parameter β by which the convergence of the sequence X2k is ensured. Eigenvalues of RA are complex numbers, in general.
10
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
Theorem 2.4 Let A ∈ Cm×n and R ∈ Cn×m , 0 ≤ s ≤ r be given. Assume that r s rank(RA) = t and λi , i = 1, . . . , t are eigenvalues of RA. Then the condition max |1 − βλi | < 1
(2.18)
1≤i≤t
is equivalent with the following: a) For each i ∈ {1, . . . , t}, Re λi > 0, and the relaxation parameter β satisfies
½ 0 < β < 2 min
¾ Re λi | i = 1, . . . , t, |λi | 6= 0 , |λi |2
(2.19)
or b) For each i ∈ {1, . . . , t}, Re λi < 0, and the relaxation parameter β satisfies
½ 2 max
Re λi | i = 1, . . . , t, |λi | 6= 0 |λi |2
¾ < β < 0.
(2.20)
Proof. According to (2.18), for each i ∈ {1, . . . , t} we get |1 − βλi | =
p (1 − β Re λi )2 + (β Im λi )2 < 1.
This implies β(β|λi |2 − 2 Re λi ) < 0.
(2.21)
There are two possible cases: a) If β > 0, then β|λi |2 − 2 Re λi < 0, so Re λi > 0, for every i ∈ {1, . . . , t} λi and β < 2 Re |λi |2 , for every λi such that |λi | 6= 0. b) If β < 0, then β|λi |2 − 2 Re λi > 0, so Re λi < 0, for every i ∈ {1, . . . , t} λi and β > 2 Re |λi |2 , for every λi such that |λi | 6= 0. Now we are looking for practical values of the parameter β. Let us denote by mRe = min{ Re λ1 , . . . , Re λt }, M Re = max{ Re λ1 , . . . , Re λt } and (M Im)2 = max{( Im λ1 )2 , . . . , ( Im λt )2 }. In accordance with conditions a) of Theorem 2.4 and condition (2.19), in the case Re λi > 0, for every i ∈ {1, . . . , t}, we can use the following value for β: mRe C βopt = . (2.22) 2 (M Re) + (M Im)2 Similarly, according to conditions b) of Theorem 2.4 and condition (2.20), in the case Re λi < 0, for every i ∈ {1, . . . , t}, we can use: C βopt =
M Re . (mRe)2 + (M Im)2
(2.23)
Successive matrix squaring algorithm for computing outer inverses
11
Theorem 2.5 Let A ∈ Crm×n and R ∈ Cn×m , 0 ≤ s ≤ r be given. Also, s suppose that rank(RA) = t and λi , i = 1, . . . , t are eigenvalues of RA. a) If Re λi > 0, for every i ∈ {1, . . . , t}, the following statements are valid: C (i) For βopt defined in (2.22) the iterative sequence (2.1) converges.
(ii) The following error estimate is valid: ¡ ¢2k−1 kX − X2k k C ≤ 1 − βopt · mRe + O(ε). kXk (iii) The number k of required iterations needful to attain the accuracy kX − X2k k ≤ δ, kXk
0 0, i = 1, . . . , t and 1 − βopt · mRe = 0 or Re λi < 0, i = 1, . . . , t and C 1 − βopt · M Re = 0 the sufficient number of iterative steps is k = 1. (2) Because of the accumulation of round off errors in the iterative steps of SMS method, we recommend to use the minimal number of steps from (2.25) or C · mRe 6= 0 or (2.27). Therefore, in the case Re λi > 0, i = 1, . . . , t and 1 − βopt C in the case Re λi < 0, i = 1, . . . , t and 1 − βopt · M Re 6= 0 we have ½ C · mRe)2 ) + F (δ)e, Re λi > 0, i = 1, . . . , t d2 − F ((1 − βopt (2.28) k= C d2 − F ((1 − βopt · M Re)2 ) + F (δ)e, Re λi < 0, i = 1, . . . , t.
3
Numerical examples
The implementation of the SMS method described in Theorem 2.1 is written in the programming package MATHEMATICA. About the package see, for example C [20]. In the next examples we use β = βopt defined in (2.22) or (2.23). Example 3.1 In this example we get exact outer inverse X of A in the first iteration as we note in the part (1) of Remark 2.2. Consider the following 6 × 5 matrix of rank 4, investigated in [9]: A=
1 1 2 3 4 6
2 3 3 4 5 6
3 4 4 5 6 7
4 6 5 6 7 7
1 2 3 4 6 8
.
In accordance with Theorem 2.3, we can choose appropriate matrices F and G in order to generate {2}-inverses F (GAF )−1 G of A. The matrices F and
Successive matrix squaring algorithm for computing outer inverses
13
G must satisfy rank(F G) = s ≤ 4, and dimensions of the matrix F G must be 6 × 5. For example, one can choose F and G using the following pattern: Fs 0
F =
,
Is
G=
0
,
where Fs As = Is and As is the principal minor which contains first s rows and first s columns of A. Let us choose s = 2. In this case the block F2 is equal to F2 =
3 −1
−2 1
.
C Applying (2.22), we get βopt = 1. Both nonzero eigenvalues of RA are equal to 1, so in this case is mRe = M Re = 1 and Im λ1 = Im λ2 = 0. As it is claimed in Remark 2.2, part (1), the outer inverse X1 is obtained by SMS algorithm (2.14) immediately, in the first iteration. Corresponding {2}-inverse X1 = F (GAF )−1 G = F G of rank 2 is equal to
X1 = X21 =
3 −1 0 0 0
−2 1 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
.
Example 3.2 Consider the matrix A from Example 3.1, and choose the following matrices F ∈ C 5×2 and G ∈ C 2×6 from [9]:
0 2 F = 3 5 1
0 1 2 3 0
,
G=
0 1
1 0
0 1
1 0
0 1
1 0
.
(3.1)
Since rank(F G) = 2 corresponding {2}-inverse X1 = F (GAF )−1 G for A of rank 2 is equal to A(2)
1 = X1 = 174
0 −21 60 39 −102
0 19 −46 −27 84
0 −21 60 39 −102
0 19 −46 −27 84
0 −21 60 39 −102
0 19 −46 −27 84
.
Let us choose the precision δ = 10−6 . According to heuristic (2.22) the C parameter β is equal to βopt = 9.208888 × 10−6 . In accordance with (2.28), the minimal number of iterative steps is equal to k = 23. An application of iterations (2.14) produces the following approximation of A(2) = X1 after 23 iterations:
X223 =
0. −0.12069 0.344828 0.224138 −0.586207
0. 0.109195 −0.264368 −0.155172 0.482759
0. −0.12069 0.344828 0.224138 −0.586207
0. 0.109195 −0.264368 −0.155172 0.482759
0. −0.12069 0.344828 0.224138 −0.586207
0. 0.109195 −0.264368 −0.155172 0.482759
.
The relative error kA(2) −X223 k/kA(2) k = kX1−X223 k/kX1k, where the matrix norm k · k represents the maximal singular value of the argument, is equal to 3.40765 × 10−11 .
14
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c Also, obtained r223 = 2., which is the exact approximation of rank(R) = 2.
Now, we generate a {1, 2, 3}-inverse of A. Applying MATHEMATICA function QRDecomposition, we obtain the following full-rank factorization A = P Q for A: P
=
Q
=
0.122169 0.122169 0.244339 0.366508 0.488678 0.733017
8.18535 0. 0. 0.
0.339377 0.7528 0.265331 0.191285 0.117239 −0.444275 9.65139 2.41883 0. 0.
0.753171 −0.393966 0.382379 0.0115872 −0.359205 0.046349
11.7283 3.64059 0.440315 0.
−0.229416 −0.114708 0.688247 −0.573539 0.344124 −0.114708
13.1943 6.05942 0.440315 0.
,
11.3618 0.555344 . −0.625711 0.458831
Using G = P T and the following, randomly chosen n × r = 5 × 4 matrix F =
1 2 −1 0 1
2 3 −1 1 2
3 4 2 −2 1
4 1 3 3 3
,
we obtain R = FG =
2.14277 4.04574 0.356548 −1.85521 0.865847
−0.0129617 0.812165 −2.00703 1.19661 0.889679
4.67513 3.50243 2.31983 1.56531 3.22212
−1.51032 0.779682 −2.25524 −1.55251 −0.959952
1.02204 0.236379 −0.291955 1.86802 1.39632
−0.475318 0.203896 −0.540167 −0.881097 −0.453308
.
According to (2.28), the optimal number of iterative steps which ensures the precision δ = 10−6 is equal to k = 20. In terms of (2.22) we compute C βopt = 0.000190063. Corresponding {1, 2, 3}-inverse of A, obtained after 20 iterations, is X220 =
0.152778 −0.652778 1.59722 −0.597222 −0.5
−0.541667 2.29167 −1.20833 −0.0416667 −0.25
1.06944 −6.56944 1.18056 1.81944 1.5
−0.722222 4.47222 −0.277778 −1.47222 −1.25
0.194444 −1.44444 −0.694444 0.944444 0.75
0.0138889 0.736111 0.236111 −0.486111 −0.25
,
and the approximation of its rank is equal to r220 = 4. If we apply standard MATHEMATICA function MatrixRank[] to compute the rank of X, we obtain MatrixRank[X]=4, which is approximated by r220 . Relative error estimate is kA(1,2,3) − X220 k/kA(1,2,3) k = kX1 − X220 k/kX1k = 1.39537 × 10−11 . In the sequel we verify the recommendation that the number of iterations should be as small as possible. After k = 45 > 20 iterations we obtain the following {1, 2, 3}-inverse of A
X245
0.151176 −0.651176 1.59882 = −0.598824 −0.5
−0.540797 2.2908 −1.2092 −0.040797 −0.25
1.06935 −6.56935 1.18065 1.81935 1.5
−0.722722 4.47272 −0.277278 −1.47272 −1.25
0.195448 −1.44545 −0.695448 0.945448 0.75
0.0135242 0.736476 0.236476 −0.486476 −0.25
.
Successive matrix squaring algorithm for computing outer inverses
15
Also, the approximation of rank(X245 ) is equal to r245 = 4.00037. MATHEMATICA function MatrixRank[] produces miscount result: MatrixRank[X245 ]=5. Since MatrixRank[R]=4, from Corollary 2.3 we conclude that rank of X is equal to 4, which is approximated by r245 . Therefore, we also show that the result proved in Corollary 2.3 is important. The Moore-Penrose inverse of A is equal to
4 −8 1 † P seudoInverse[A] = A = 10 8 −2 −4
−1 15 −13 3 −2
−8 −36 26 −2 12
7 23 −15 1 −10
−5 −5 1 1 6
3 3 −1 −1 −2
,
where Pseudoinverse[A] is standard MATHEMATICA function for computing A† (see, for example [20]) On the other hand, applying iterations (2.14) in the case R = F G = A∗ , we obtain an approximation of A† . In accordance with (2.28), the minimal number of iterative steps, corresponding to δ = 10−6 , is equal to k = 36, for C = 4.318278 × 10−8 . The approximation of the Moore-Penrose inverse of A βopt is X236 =
0.500002 −1. 1.25 −0.249997 −0.5
−0.125004 1.875 −1.625 0.374996 −0.25
−0.999991 −4.50001 3.24999 −0.249991 1.5
0.874994 2.87501 −1.87499 0.124994 −1.25
−0.624999 −0.625001 0.124999 0.125001 0.75
0.374999 0.375001 −0.124999 −0.125001 −0.25
.
In this example we again show that information r2k → rank(X) is important. An application of the standard MATHEMATICA function MatrixRank[X236 ] gives the result 5. However, since MatrixRank[R]=4, we conclude that MATHEMATICA is unable to compute the exact rank of X236 . On the other hand, we obtain r236 = 4., which is the desired value in accordance with Corollary 2.3. As a confirmation, we observe that MatrixRank[Pseudoinverse[A]]=4. Relative error estimate is kA† − X236 k/kA† k = 3.12183 × 10−6 . Example 3.3 Consider the following matrix: A=
1 −1 −1 −1 −1 −1
−1 1 −1 −1 −1 −1
0 0 1 −1 −1 0
0 0 −1 1 0 −1
0 0 0 0 2 −1
0 0 0 0 −1 2
.
Using R = A2 and β = 0.00137174, derived from (2.22), SMS method produces well known approximation with the precision δ = 10−6 of the Drazin inverse from [19] after k = 15 iterations: X215 =
0.25 −0.25 0. 0. 0. 0.
−0.25 0.25 0. 0. 0. 0.
0. 0. 0.25 −0.25 −0.416667 −0.583333
0. 0. −0.25 0.25 −0.583333 −0.416667
0. 0. 0. 0. 0.666667 0.333333
0. 0. 0. 0. 0.333333 0.666667
.
16
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
The exact Drazin inverse is equal to
A
D
=
1 4 − 14
0 0 0 0
− 14
0 0
0 0
0 0 0 0
1 4 − 14 5 − 12 7 − 12
− 14
1 4
1 4 7 − 12 5 − 12
0 0 0 0
0 0 0 0
2 3 1 3
1 3 2 3
.
The relative error kAD − X215 k/kAD k, where the matrix norm k · k represents the maximal singular value of the argument, is equal to 2.67633 × 10−14 . Let us mention that in [19] Euler-Knopp representation, corresponding to iterations (2.6) in the case R = Al , l ≥ ind(A), gives an analogous representation in 29697 iterations. Rank of the matrix R can be computed using r215 = 4. Example 3.4 In this example we verify SMS algorithm (2.14) in the case Re λi < 0, for each i = 1, . . . , t. For this purpose, we accomplished the SMS algorithm for the matrix −A, where A is defined as in Example 3.1. We used C βopt = −9.208888 × 10−6 defined by (2.23) and the number of iterative steps k = 23 conditioned by the second case of (2.28). Corresponding outer inverse, generated by the SMS algorithm is equal to −X223 , where X223 is the outer inverse of A generated in Example 3.1.
4
Conclusion
SMS algorithm is investigated in the papers [3], [16], [18]. In these articles SMS algorithm is based on successive squaring of the appropriate composite matrix T in the form (2.2). In the paper [3] using that P = I−βA∗ A, Q = βA∗ the MoorePenrose inverse is derived. The weighted Moore-Penrose inverse is derived in [18], using P = I − βA] A, Q = βA] , and finally the Drazin inverse is derived in [16] applying the SMS iterative process with the matrices P = I − βAl A, Q = βAl . In the present article we introduced and investigated SMS iterative scheme based on the matrices P = I − βRA and Q = βR. In this way, we derive an universal SMS algorithm containing all previous results as partial cases. Moreover, we prove that our method can be used to approximate various outer . Properties of the generated outer inverse inverses of a given matrix A ∈ Cm×n r X corresponding to chosen matrix R are examined. Moreover, in the case when a full-rank factorization R = F G is known, we give an explicit form of an outer generalized inverse of a A corresponding to the chosen R. An explicit approximation of rank(R) and rank(X) is also given. We implemented our algorithm in the programming package MATHEMATICA, and many numerical examples are presented in the last section. In these examples we show that sometimes MATHEMATICA is unable to compute the exact value of rank(X), so its iterative approximation is important.
Successive matrix squaring algorithm for computing outer inverses
17
References [1] A. Ben-Israel and T.N.E. Greville, Generalized inverses: theory and applications, Second Ed., Springer, 2003. [2] S.R.Caradus Generalized Inverses and Operator Theory, Queen’s Paper in Pure and Applied Mathematics, Queen’s University, Kingston, Ontario, 1978. [3] L. Chen, E.V. Krishnamurthy, I. Macleod, Generalized matrix inversion and rank compuation by successive matrix powering, Parallel Computing 20 (1994) 297–311. [4] D.S. Djordjevi´c and P.S. Stanimirovi´c, General representations of pseudoinverses, Matematiˇcki Vesnik 51 (1999), 69-76. [5] D.S. Djordjevi´c, P.S. Stanimirovic and Y. Wei, The representation and approximation of outer generalized inverses, Acta Math. Hungar, 104 (2004), 1–26. [6] R.E. Harte, Invertibility and singularity for bounded linear operator, Dekker 1988. [7] R.A. Horn and C.R. Johnson, Matrix Analysis, CambrodgrUniversity Press, Cambridge, New York, New Rochelle, Melbourne, Sydney, 1986. [8] C.R. Rao and S.K. Mitra, Generalized Inverse of Matrices and its Applications, John Wiley Sons, Inc, New York, London, Sydney, Toronto, 1971. ´ c, Adjoint mappings and inverses [9] P.S. Stanimirovi´c, S. Bogdanovi´c, M. Ciri´ of matrices, Algebra Colloquium 13(3) (2006) 421–432. [10] P.S. Stanimirovi´c, Block representations of {2}, {1,2} inverses and the Drazin inverse, Indian Journal Pure Appl. Math., 29 (1998) 1159–1176. [11] P.S. Stanimirovi´c, Applications of hyper-power method for computing matrix products, Publikacije Elektrotehniˇckog fakulteta, 15 (2004) 13–25. [12] P.S. Stanimirovi´c and D.S. Djordjevi´c, Full-rank and determinantal representation of the Drazin inverse, Linear Algebra Appl., 311 (2000) 31–51. [13] G. Wang, Y. Wei, S. Qiao, Generalized inverses: theory and computations, Science Press, 2003. (1,2)
[14] G. Wang, The representations of the generalized inverses (A ⊗ B)T,S and (2)
(A⊗ B)T,S and some applications, J. Shanghau Univ. (Natural Sciences) 24 (1995) 1–6. [15] Y. Wei, A characterization and representation for the generalized inverse (2) AT,S and its applications, Linear Algebra Appl., 280 (1998), 87–96.
18
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
[16] Y. Wei, Successive matrix squaring algorithm for computing Drazin inverse, Appl. Math. Comp. 108 (2000), 67–75. [17] Y. Wei and H. Wu, The representation and approximation for the general(2) ized inverse AT,S , Appl. Math. Comput., 135 (2003), 263–276. [18] Y. Wei, H. Wu, J. Wei Successive matrix squaring algorithm for parallel computing the weighted generalized inverse A†M N , Appl. Math. Comp. 116 (2000) 289–296. [19] Y. Wei and H. Wu, The representation and approximation for Drazin inverse, J. Comput. Appl. Math. 126 (2000), 417–423. [20] S. Wolfram, The Mathematica Book, 4th ed., Wolfram Media/Cambridge University Press, 1999.
Abstract In this paper we derive a successive matrix squaring (SMS) algorithm to approximate an outer generalized inverse with prescribed range and null space of a given matrix A ∈ Cm×n . We generalize the results from r the papers [3], [16], [18], and obtain an algorithm for computing various classes of outer generalized inverses of A. Instead of particular matrices used in these articles, we use an appropriate matrix R ∈ Cn×m , s ≤ r. s Numerical examples are presented. AMS Subj. Class.: 15A09. Key words: Generalized inverse; outer inverse; SMS algorithm; full rank factorization; matrix rank.
1
Intoduction and preliminaries
Let Cm×n and Cm×n denote the set of all complex m × n matrices and all r complex m × n matrices of rank r, respectively. In denotes the unit matrix of order n. By A∗ , R(A), rank(A) and N (A) we denote the conjugate transpose, the range, the rank and the null space of A ∈ Cm×n . By Re z and Im z we denote a real and imaginary part of a complex number z, respectively. For A ∈ Cm×n , the set of inner and outer generalized inverses are defined by the following, respectively: A{1} = {X ∈ Cn×m | AXA = A},
A{2} = {X ∈ Cn×m | XAX = X}.
The set of all outer inverses with prescribed rank s is denoted by A{2}s , 0 ≤ s ≤ r = rank(A). The symbols A− or A(1) stand for an arbitrary generalized ∗ Corresponding
author by Grant No. 144003 of the Ministry of Science, Technology and Development, Republic of Serbia. † Supported
1
2
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
inner inverse of A and by A(2) we denote an arbitrary generalized outer inverse of A. Also, the matrix X which satisfies AXA = A
and
XAX = X
is called the reflexive g-inverse of A and it is denoted by A(1,2) . The set of all reflexive g-inverses is denoted by A{1, 2}. Subsequently, the sets of {1, 2, 3} and {1, 2, 4} inverses of A are defined by A{1, 2, 3} = A{1, 2} ∩ {X| (AX)∗ = AX}, A{1, 2, 4} = A{1, 2} ∩ {X| (XA)∗ = XA}. By A† we denote the Moore-Penrose inverse of A, i.e. the unique matrix A† satisfying AA† A = A, A† AA† = A† , (AA† )∗ = AA† , (A† A)∗ = A† A. For A ∈ Cn×n the smallest nonnegative integer k such that rank(Ak+1 ) = rank(Ak ) is called the index of A and denoted by ind(A). If A ∈ Cn×n is a square matrix with ind(A) = k, then the matrix X ∈ Cn×n which satisfies the following conditions Ak XA = Ak , XAX = X, AX = XA, is called the Drazin inverse of A and it is denoted by AD . When ind(A) = 1, Drazin inverse AD is called the group inverse and it is denoted by A# . Suppose that M and N are Hermite positive definite matrices of the order m and n, respectively. Then there exists the unique matrix X ∈ Cn×m such that AXA = A, XAX = X, (M AX)∗ = M AX, (N XA)∗ = N XA. The matrix X is called the weighted Moore-Penrose inverse of A, and denoted by X = A†M,N . In particular, if M = Im and N = In , then A†M,N = A† . If A ∈ Cn×m and W ∈ Cm×n , then the unique solution X ∈ Cn×m of the equations (AW )k+1 XW = (AW )k ,
XW AW X = X,
AW X = XW A,
(1.1)
where k = ind(AW ), is called the W-weighted Drazin inverse of A and it is denoted by AD,W . If A ∈ Cm×n , T is a subspace of Cn of dimension t ≤ r and S is a subspace r m of C of dimension m − t, then A has a {2} inverse X such that R(X) = V and N (X) = U if and only if AV ⊕ U = Cm , (2)
in which case X is unique and we denote it by AV,U .
Successive matrix squaring algorithm for computing outer inverses
3
It is well-known that for A ∈ Cm×n , the Moore-Penrose A† , the weighted Moore-Penrose inverse A†M,N and the weighted Drazin inverse AD,W can be represented by: (2)
(i) A† = AR(A∗ ),N (A∗ ) , (2)
(ii) A†M,N = AR(A] ),N (A] ) , where A] = N −1 A∗ M, (2)
(iii) AD,W = (W AW )R(A(W A)k ),N (A(W A)k ) , where W ∈ Cn×n , k = ind(W A). Also, for A ∈ Cn×n , the Drazin inverse AD can be represented by: (2)
AD = AR(Ak ),N (Ak ) , where ind(A) = k. The following representations of {2, 3}, {2, 4}-inverses with prescribed rank s are restated from [11]: Proposition 1.1 Let A ∈ Cm×n and 0 < s < r be a chosen integer. Then the r following is valid: © ª (a) A{2, 4}s = (ZA)† Z | Z ∈ Cs×m , ZA ∈ Cs×n ; s © ª † n×s m×s (b) A{2, 3}s = Y (AY ) | Y ∈ C , AY ∈ Cs . General representations for various classes of generalized inverses can be found in [4, 8, 10, 12]. Some of these representations are restated here for the sake of completeness. Proposition 1.2 Let A ∈ Crm×n be an arbitrary matrix and A = P Q is a full-rank factorization of A. There are the following general representations for some classes of generalized inverses: A{2}s = {F (GAF )−1 G | F ∈ Cn×s, G ∈ Cs×m, rank(GAF ) = s}; r [ A{2} = A{2}s ; s=0
A{1, 2} = {F (GAF )−1 G | F ∈ Cn×r, G ∈ Cr×m, rank(GAF ) = r} = A{2}r ; A{1, 2, 3} = {F (P ∗ AF )−1 P ∗ | F ∈ Cn×r , rank(P ∗ AF ) = r}; A{1, 2, 4} = {Q∗ (GAQ∗ )−1 G | G ∈ Cr×m, rank(GAQ∗ ) = r}; A† = Q∗ (P ∗ AQ∗ )−1 P ∗ ; AD = PAl (QAl A PAl )−1 QAl , Al = PAl QAl , l ≥ ind(A). For other important properties of generalized inverses see [1, 2, 6, 13]. We will use the following well-known result: Lemma 1.1 [7] Let M ∈ Cn×n and ε > 0 be given. There is at least one matrix norm k · k such that ρ(M ) ≤ kM k ≤ ρ(M ) + ², (1.2) where ρ(M ) denotes the spectral radius of M .
4
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
The basic motivation for this paper were the paper of L. Chen et al. [3] and the papers of Y. Wei [16] and Y. Wei et al. [18]. In the paper of L. Chen et al. [3] a successive matrix squaring (SMS) algorithm which approximates the Moore-Penrose inverse of a given matrix A were considering, while in the papers of Y. Wei [16] and Y. Wei et al. [18] the authors considered two variants of SMS algorithm which approximate the Drazin inverse and the weighted MoorePenrose inverse of A, respectively. In this paper we derive a SMS algorithm to approximate an outer generalized inverse with prescribed range and null space of a given matrix A ∈ Cm×n , r which is based on successive squaring of an appropriate block matrix T . We generalized the results from the papers [3], [16] and [18]. Instead of matrices A∗ , Ak , k ≥ ind(A) and A] , used to compose the matrix T in these articles, , 0 ≤ s ≤ r. Furthermore, we use an universal appropriate matrix R ∈ Cn×m s we give an explicit form of the outer generalized inverse of A corresponding to chosen R. An explicit approximation of rank(X) is also given in the same iterative scheme. We also estimate the number of iterative steps necessary to obtain a prescribed precision. Numerical examples in the last section show that numerical approximation of rank(X) is important and that the number of SMS iterations is less than in the classical Euler-Knopp iterations.
2
Results
Let A ∈ Cm×n and R ∈ Cn×m , 0 ≤ s ≤ r be given. The general iterative r s scheme used in this paper is given by X1 = Q, Xk+1 = P Xk + Q, k ∈ N,
(2.1)
where P = I − βRA, Q = βR and β is a relaxation parameter. The iterative process (2.1) can be accelerated by means of (m + n) × (m + n) composite matrix T , partitioned in the block form · ¸ P Q T = . (2.2) 0 I The improvement of (2.1) can be done by computing the matrix power · Pk−1 i ¸ Pk i=0 P Q . Tk = 0 I Pk−1 It is not difficult to see that the iterative scheme (2.1) gives Xk = i=0 P i Q. k Hence, the matrix Xk is equal to the right upper block in T k . In turn, T 2 can be computed by k repeated squaring, that is T0 Ti+1
= =
T Ti2 , i = 0, 1, . . . , k − 1.
(2.3)
Successive matrix squaring algorithm for computing outer inverses It is clear that
" Tk = T
2k
P2 0
=
k
P2k −1 i=0
5
# P iQ
I
.
(2.4)
Now, we will state an auxiliary result: Lemma 2.1 If P ∈ Cn×n and S ∈ Cn×n are such that P = P 2 and P S = SP then ρ(P S) ≤ ρ(S). Proof. Suppose that λ is an eigenvalue of P S and x is a corresponding vector such that P Sx = λx. Then, λP x = λx which implies that λ = 0 or P x = x. If P x = x, then Sx = SP x = P Sx = λx, i.e. λ is an eigenvalue of S. Hence, the set of eigenvalues of P S is the subset of the set of eigenvalues of S union {0}, which implies that ρ(P S) ≤ ρ(S). The following theorem represents the analogue of Theorem 3.1 [15] which deals with V -norm, where V-norm is defined as kBkV = kBV k2 , for any B ∈ Cn×m and V is invertible such that V −1 AGV is the Jordan form of AG. The method used in our proof is different then the one used in Theorem 3.1 [15]. and R ∈ Csn×m , 0 ≤ s ≤ r be given such that Theorem 2.1 Let A ∈ Cm×n r AR(R) ⊕ N (R) = Cm .
(2.5)
The sequence of approximations X2k =
k 2X −1
(I − βRA)i βR
(2.6)
i=0
determined by the SMS algorithm (2.3) converges in the matrix norm k · k to the (2) outer inverse X = AR(R),N (R) of A if β is a fixed real number such that max |1 − βλi | < 1,
1≤i≤t
(2.7)
where rank(RA) = t, λi , i = 1, . . . , t are eigenvalues of RA and k · k satisfies condition (1.2) from Lemma 1.1 for M = I − βAR. In the case of convergence we have the following error estimate k kX − X2k k ≤ max |1 − βλi |2 + O(ε), k ≥ 0. 1≤i≤t kXk
(2.8)
Proof. From the condition (2.5), it follows the existence of the outer inverse (2) X = AR(R),N (R) of A. It is evident that XAR = R and R = RAX. Since R(X2k ) ⊆ R(R), we have that XAX2k = X2k .
6
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
So, we obtain kX − X2k k = kXAX − XAX2k k = kX(AX − AX2k )k ≤ kXk · kAX − AX2k k.
(2.9)
Now, let us prove that AX − AXk = (AX − AX1 )k .
(2.10)
We will use the mathematical induction. For k = 1, the equality (2.10) is true. Suppose that (2.10) holds for k = p−1 and prove that it holds for k = p. Using P = I − QA, we have AX − AXp
= = =
AX − A(P Xp−1 + Q) AX − A(I − QA)Xp−1 − AQ AX − AXp−1 + AQAXp−1 − AQ.
Since R = RAX implies Q = QAX, we get AX − AXp
= AXAX − AXAXp−1 + AQAXp−1 − AQAX = (AX − AQ)(AX − AXp−1 ) = (AX − AQ)p = (AX − AX1 )p .
Now, (2.9) and (2.10) imply k
kX − X2k k ≤ kXk · kAX − AX2k k ≤ kXk · kAX − AX1 k2 . If (2.7) is satisfied, we can choose ² such that max1≤i≤t |1 − βλi | + ε < 1. For such ² and M = I − βAR, there exists a matrix norm k · k satisfies (1.2). Now, µ ¶2k k kX − X2k k ≤ kAX − AX1 k2 ≤ max |1 − βλi | + ε . 1≤i≤t kXk The last inequality follows if we apply Lema 2.1 for P = AX and S = I −AX1 = I − βAR and by Lemma 1.1. Remark that a similar problem as in Theorem 3.1 [15] was discussed in [17] but in a more general case: If {Sn (x)} is a family of continuous real valued functions on open set Ω such that σ((RA)|T ) ⊂ Ω ⊂ (0, ∞) and limn→∞ Sn (x) = 1/x uniformly on σ((RA)|T ), then (2)
AT,S = lim Sn ((RA)|T )R n→∞
and (2)
kSn ((RA)|T )R − AT,S kV (2)
kAT,S kV
≤
max
x∈σ((RA)|T )
|Sn (x)x − 1| + O(ε), ² > 0.
The generalization to the bounded linear operator on Banach space of the mentioned result proved in [17] is given in [5].
Successive matrix squaring algorithm for computing outer inverses
7
Corollary 2.1 Under the assumption and the conditions of Theorem 2.1, (i) In the case m = n, R = Al , l ≥ ind(A), we have AD = lim X2k , k→∞
(ii) In the case R = A∗ , A† = lim X2k , k→∞
(iii) In the case R = A] = N −1 A∗ M , A†M,N = lim X2k , k→∞
(iv) In the case R = A(W A)k , AD,W = lim X2k . k→∞
Corollary 2.2 In the case rank(R) = rank(A), under the conditions of Theorem 2.1, we have X ∈ A{1, 2}. Proof. Using that rank(X) = rank(A), by Corollary 1 (pg. 19, [1]) it follows that X ∈ A{1, 2}. Theorem 2.2 Let A ∈ Rrm×n and R ∈ Rsn×m , 0 ≤ s ≤ r be real matrices and assume that β satisfies condition (2.7). Then the sequence of real numbers defined by k
k
r2k = n − Tr((I − βRA)2 ) = n − Tr(P 2 ), k ≥ 0
(2.11)
converges to rank(RA). Proof. It is sufficient to use k
Tr((I − βRA)2 ) = n − t +
t X k (1 − βλi )2 ,
(2.12)
i=1
where t = rank(RA) ≤ rank(R). Therefore, using (2.12), we get r 2k = t −
t X
k
(1 − βλi )2 .
(2.13)
i=1
Now, the proof can be completed using (2.7). Remark 2.1 Taking R = A∗ in Theorem 2.1 and Theorem 2.2 we get the results of Theorem 2 from [3]. Also, taking R = Al , l ≥ ind(A) and R = A] in Theorem 2.1, respectively, we get the results from [16] and [18]. The next theorem presents an explicit representation of the outer inverse of A which is a limit value of the sequence of approximations (2.1), for given A and R.
8
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
Theorem 2.3 Assume that A∈Cm×n . Consider an arbitrary matrix R∈Cn×m , r s s×m 0 ≤ s ≤ r, and its full rank factorization R = F G, F ∈ Cn×s , G ∈ C , such s s that GAF is invertible. The sequence of approximations X2k =
k 2X −1
(I − βRA)i βR
(2.14)
i=0
determined by the SMS algorithm (2.3) converges in the matrix norm to the outer inverse X = F (GAF )−1 G of A, if β is a fixed real number satisfying max |1−βλi | < 1,
(2.15)
1≤i≤s
where λi , i = 1, . . . , s are eigenvalues for F GA. In the case of convergence we have the error estimate k kX − X2k k ≤ max |1 − βλi |2 + O(ε), k ≥ 0, 1≤i≤s kXk
(2.16)
where the matrix norm k · k satisfies condition (1.2), for M = I − βAR. Proof. Since the matrices F and G satisfy rank(GAF ) = rank(F ) = rank(G), according to Theorem 1.3.8 (pg. 33, [13]) and [14], we conclude that X = F (GAF )−1 G is a {2} inverse of A having the range R(X) = R(F ) = V and null space N (X) = N (G) = U , i.e. (2)
(2)
X = AV,U = AR(F ),
N (G) .
Also, conditions conditions AV ∩ U = {0} and AV ⊕ U = Cm for the existence of X are ensured (Theorem 1.3.8, pg. 33, [13]). On the other hand, since rank(R) = s = rank(F G) = rank(F ) = rank(G) (2)
we obtain R(F ) = R(R), N (G) = N (R) which implies X = AR(R),N (R) . In view of (2.15) immediately follows that the conditions of Theorem 2.1 are satisfied. Therefore, we may conclude that the sequence of approximations (2.14) determined by the SMS algorithm (2.3) converges to X = F (GAF )−1 G, in the matrix norm satisfying (1.2) for M = I − βAR. Corollary 2.3 Assume that A ∈ Rm×n and R ∈ Rn×m , 0 ≤ s ≤ r, are two r s arbitrary chosen real matrices. Let R = F G, F ∈ Cn×s , G ∈ Cs×m be a s s full rank factorization of R, such that GAF is invertible. The sequence of real numbers defined by k
k
r2k = n − Tr((I − βRA)2 ) = n − Tr(P 2 ), k ≥ 0 converges to rank(X).
(2.17)
Successive matrix squaring algorithm for computing outer inverses
9
Proof. According to Theorem 2.2, r2k → rank(RA). Since rank(RA) = rank(GAF ) = s = rank(R) = rank(X), we complete the proof. Corollary 2.4 Consider A ∈ Cm×n and arbitrary chosen matrix R ∈ Cn×m , r s n×s s ≤ r. Assume that R = F G, F ∈ Cs , G ∈ Css×m is a full-rank factorization of R. Also, suppose that condition (2.15) of Theorem 2.3 are valid and that A = P Q is a full rank factorization of A. If the sequences of approximations X2k is defined by (2.14), the following hold: ½ ¾ s×m a) lim X2k | F ∈ Cn×s , G ∈ C , rank(GAF ) = s, 0 ≤ s ≤ r ⊂ A{2}; s s k→∞
½ b)
¾ n×r
lim X2k | F ∈ C
k→∞
,G ∈ C
r×m
, rank(GAF ) = r
½ c)
¾ n×r
lim X2k | F ∈ C
k→∞
∗
, G = P , rank(GAF ) = r
½ d)
⊂ A{1, 2};
⊂ A{1, 2, 3}; ¾
∗
r×m
lim X2k | F = Q , G ∈ C
k→∞
, rank(GAF ) = r
⊂ A{1, 2, 4};
satisfy rank(GAF ) = s, Proof. a) If the matrices F ∈ Csn×s and G ∈ Cs×m s applying Theorem 2.3, we conclude that X = F (GAF )−1 G is the limit value of the sequence of approximations (2.14). Using general representation of outer inverses form Proposition 1.2, we immediately obtain X ∈ A{2}. Parts b)-d) can be proved using Theorem 2.3 and general representations of pseudoinverses from Proposition 1.2. and integer s such that 0 ≤ s ≤ r. Corollary 2.5 Consider matrix A ∈ Cm×n r Then the following is valid: ( Ã k−1 ! ) 2P i ∗ ∗ n×s a) Y · lim (I − β(AY ) AY ) (AY ) |Y ∈C , rank(AY ) = s k→∞
i=0
⊂ A{2, 3}s ; ( Ã k−1 ! ) 2P i b) lim (I − β(ZA)∗ ZA) (ZA)∗ · Z | Z ∈ Cs×n , rank(ZA) = s k→∞
i=0
⊂ A{2, 4}s . Proof. Follows from Corollary 2.1, part (ii) and the general representations of the sets A{2, 3}s and A{2, 4}s given in Proposition 1.1. In the following theorem we investigate the limits of the relaxation parameter β by which the convergence of the sequence X2k is ensured. Eigenvalues of RA are complex numbers, in general.
10
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
Theorem 2.4 Let A ∈ Cm×n and R ∈ Cn×m , 0 ≤ s ≤ r be given. Assume that r s rank(RA) = t and λi , i = 1, . . . , t are eigenvalues of RA. Then the condition max |1 − βλi | < 1
(2.18)
1≤i≤t
is equivalent with the following: a) For each i ∈ {1, . . . , t}, Re λi > 0, and the relaxation parameter β satisfies
½ 0 < β < 2 min
¾ Re λi | i = 1, . . . , t, |λi | 6= 0 , |λi |2
(2.19)
or b) For each i ∈ {1, . . . , t}, Re λi < 0, and the relaxation parameter β satisfies
½ 2 max
Re λi | i = 1, . . . , t, |λi | 6= 0 |λi |2
¾ < β < 0.
(2.20)
Proof. According to (2.18), for each i ∈ {1, . . . , t} we get |1 − βλi | =
p (1 − β Re λi )2 + (β Im λi )2 < 1.
This implies β(β|λi |2 − 2 Re λi ) < 0.
(2.21)
There are two possible cases: a) If β > 0, then β|λi |2 − 2 Re λi < 0, so Re λi > 0, for every i ∈ {1, . . . , t} λi and β < 2 Re |λi |2 , for every λi such that |λi | 6= 0. b) If β < 0, then β|λi |2 − 2 Re λi > 0, so Re λi < 0, for every i ∈ {1, . . . , t} λi and β > 2 Re |λi |2 , for every λi such that |λi | 6= 0. Now we are looking for practical values of the parameter β. Let us denote by mRe = min{ Re λ1 , . . . , Re λt }, M Re = max{ Re λ1 , . . . , Re λt } and (M Im)2 = max{( Im λ1 )2 , . . . , ( Im λt )2 }. In accordance with conditions a) of Theorem 2.4 and condition (2.19), in the case Re λi > 0, for every i ∈ {1, . . . , t}, we can use the following value for β: mRe C βopt = . (2.22) 2 (M Re) + (M Im)2 Similarly, according to conditions b) of Theorem 2.4 and condition (2.20), in the case Re λi < 0, for every i ∈ {1, . . . , t}, we can use: C βopt =
M Re . (mRe)2 + (M Im)2
(2.23)
Successive matrix squaring algorithm for computing outer inverses
11
Theorem 2.5 Let A ∈ Crm×n and R ∈ Cn×m , 0 ≤ s ≤ r be given. Also, s suppose that rank(RA) = t and λi , i = 1, . . . , t are eigenvalues of RA. a) If Re λi > 0, for every i ∈ {1, . . . , t}, the following statements are valid: C (i) For βopt defined in (2.22) the iterative sequence (2.1) converges.
(ii) The following error estimate is valid: ¡ ¢2k−1 kX − X2k k C ≤ 1 − βopt · mRe + O(ε). kXk (iii) The number k of required iterations needful to attain the accuracy kX − X2k k ≤ δ, kXk
0 0, i = 1, . . . , t and 1 − βopt · mRe = 0 or Re λi < 0, i = 1, . . . , t and C 1 − βopt · M Re = 0 the sufficient number of iterative steps is k = 1. (2) Because of the accumulation of round off errors in the iterative steps of SMS method, we recommend to use the minimal number of steps from (2.25) or C · mRe 6= 0 or (2.27). Therefore, in the case Re λi > 0, i = 1, . . . , t and 1 − βopt C in the case Re λi < 0, i = 1, . . . , t and 1 − βopt · M Re 6= 0 we have ½ C · mRe)2 ) + F (δ)e, Re λi > 0, i = 1, . . . , t d2 − F ((1 − βopt (2.28) k= C d2 − F ((1 − βopt · M Re)2 ) + F (δ)e, Re λi < 0, i = 1, . . . , t.
3
Numerical examples
The implementation of the SMS method described in Theorem 2.1 is written in the programming package MATHEMATICA. About the package see, for example C [20]. In the next examples we use β = βopt defined in (2.22) or (2.23). Example 3.1 In this example we get exact outer inverse X of A in the first iteration as we note in the part (1) of Remark 2.2. Consider the following 6 × 5 matrix of rank 4, investigated in [9]: A=
1 1 2 3 4 6
2 3 3 4 5 6
3 4 4 5 6 7
4 6 5 6 7 7
1 2 3 4 6 8
.
In accordance with Theorem 2.3, we can choose appropriate matrices F and G in order to generate {2}-inverses F (GAF )−1 G of A. The matrices F and
Successive matrix squaring algorithm for computing outer inverses
13
G must satisfy rank(F G) = s ≤ 4, and dimensions of the matrix F G must be 6 × 5. For example, one can choose F and G using the following pattern: Fs 0
F =
,
Is
G=
0
,
where Fs As = Is and As is the principal minor which contains first s rows and first s columns of A. Let us choose s = 2. In this case the block F2 is equal to F2 =
3 −1
−2 1
.
C Applying (2.22), we get βopt = 1. Both nonzero eigenvalues of RA are equal to 1, so in this case is mRe = M Re = 1 and Im λ1 = Im λ2 = 0. As it is claimed in Remark 2.2, part (1), the outer inverse X1 is obtained by SMS algorithm (2.14) immediately, in the first iteration. Corresponding {2}-inverse X1 = F (GAF )−1 G = F G of rank 2 is equal to
X1 = X21 =
3 −1 0 0 0
−2 1 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
.
Example 3.2 Consider the matrix A from Example 3.1, and choose the following matrices F ∈ C 5×2 and G ∈ C 2×6 from [9]:
0 2 F = 3 5 1
0 1 2 3 0
,
G=
0 1
1 0
0 1
1 0
0 1
1 0
.
(3.1)
Since rank(F G) = 2 corresponding {2}-inverse X1 = F (GAF )−1 G for A of rank 2 is equal to A(2)
1 = X1 = 174
0 −21 60 39 −102
0 19 −46 −27 84
0 −21 60 39 −102
0 19 −46 −27 84
0 −21 60 39 −102
0 19 −46 −27 84
.
Let us choose the precision δ = 10−6 . According to heuristic (2.22) the C parameter β is equal to βopt = 9.208888 × 10−6 . In accordance with (2.28), the minimal number of iterative steps is equal to k = 23. An application of iterations (2.14) produces the following approximation of A(2) = X1 after 23 iterations:
X223 =
0. −0.12069 0.344828 0.224138 −0.586207
0. 0.109195 −0.264368 −0.155172 0.482759
0. −0.12069 0.344828 0.224138 −0.586207
0. 0.109195 −0.264368 −0.155172 0.482759
0. −0.12069 0.344828 0.224138 −0.586207
0. 0.109195 −0.264368 −0.155172 0.482759
.
The relative error kA(2) −X223 k/kA(2) k = kX1−X223 k/kX1k, where the matrix norm k · k represents the maximal singular value of the argument, is equal to 3.40765 × 10−11 .
14
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c Also, obtained r223 = 2., which is the exact approximation of rank(R) = 2.
Now, we generate a {1, 2, 3}-inverse of A. Applying MATHEMATICA function QRDecomposition, we obtain the following full-rank factorization A = P Q for A: P
=
Q
=
0.122169 0.122169 0.244339 0.366508 0.488678 0.733017
8.18535 0. 0. 0.
0.339377 0.7528 0.265331 0.191285 0.117239 −0.444275 9.65139 2.41883 0. 0.
0.753171 −0.393966 0.382379 0.0115872 −0.359205 0.046349
11.7283 3.64059 0.440315 0.
−0.229416 −0.114708 0.688247 −0.573539 0.344124 −0.114708
13.1943 6.05942 0.440315 0.
,
11.3618 0.555344 . −0.625711 0.458831
Using G = P T and the following, randomly chosen n × r = 5 × 4 matrix F =
1 2 −1 0 1
2 3 −1 1 2
3 4 2 −2 1
4 1 3 3 3
,
we obtain R = FG =
2.14277 4.04574 0.356548 −1.85521 0.865847
−0.0129617 0.812165 −2.00703 1.19661 0.889679
4.67513 3.50243 2.31983 1.56531 3.22212
−1.51032 0.779682 −2.25524 −1.55251 −0.959952
1.02204 0.236379 −0.291955 1.86802 1.39632
−0.475318 0.203896 −0.540167 −0.881097 −0.453308
.
According to (2.28), the optimal number of iterative steps which ensures the precision δ = 10−6 is equal to k = 20. In terms of (2.22) we compute C βopt = 0.000190063. Corresponding {1, 2, 3}-inverse of A, obtained after 20 iterations, is X220 =
0.152778 −0.652778 1.59722 −0.597222 −0.5
−0.541667 2.29167 −1.20833 −0.0416667 −0.25
1.06944 −6.56944 1.18056 1.81944 1.5
−0.722222 4.47222 −0.277778 −1.47222 −1.25
0.194444 −1.44444 −0.694444 0.944444 0.75
0.0138889 0.736111 0.236111 −0.486111 −0.25
,
and the approximation of its rank is equal to r220 = 4. If we apply standard MATHEMATICA function MatrixRank[] to compute the rank of X, we obtain MatrixRank[X]=4, which is approximated by r220 . Relative error estimate is kA(1,2,3) − X220 k/kA(1,2,3) k = kX1 − X220 k/kX1k = 1.39537 × 10−11 . In the sequel we verify the recommendation that the number of iterations should be as small as possible. After k = 45 > 20 iterations we obtain the following {1, 2, 3}-inverse of A
X245
0.151176 −0.651176 1.59882 = −0.598824 −0.5
−0.540797 2.2908 −1.2092 −0.040797 −0.25
1.06935 −6.56935 1.18065 1.81935 1.5
−0.722722 4.47272 −0.277278 −1.47272 −1.25
0.195448 −1.44545 −0.695448 0.945448 0.75
0.0135242 0.736476 0.236476 −0.486476 −0.25
.
Successive matrix squaring algorithm for computing outer inverses
15
Also, the approximation of rank(X245 ) is equal to r245 = 4.00037. MATHEMATICA function MatrixRank[] produces miscount result: MatrixRank[X245 ]=5. Since MatrixRank[R]=4, from Corollary 2.3 we conclude that rank of X is equal to 4, which is approximated by r245 . Therefore, we also show that the result proved in Corollary 2.3 is important. The Moore-Penrose inverse of A is equal to
4 −8 1 † P seudoInverse[A] = A = 10 8 −2 −4
−1 15 −13 3 −2
−8 −36 26 −2 12
7 23 −15 1 −10
−5 −5 1 1 6
3 3 −1 −1 −2
,
where Pseudoinverse[A] is standard MATHEMATICA function for computing A† (see, for example [20]) On the other hand, applying iterations (2.14) in the case R = F G = A∗ , we obtain an approximation of A† . In accordance with (2.28), the minimal number of iterative steps, corresponding to δ = 10−6 , is equal to k = 36, for C = 4.318278 × 10−8 . The approximation of the Moore-Penrose inverse of A βopt is X236 =
0.500002 −1. 1.25 −0.249997 −0.5
−0.125004 1.875 −1.625 0.374996 −0.25
−0.999991 −4.50001 3.24999 −0.249991 1.5
0.874994 2.87501 −1.87499 0.124994 −1.25
−0.624999 −0.625001 0.124999 0.125001 0.75
0.374999 0.375001 −0.124999 −0.125001 −0.25
.
In this example we again show that information r2k → rank(X) is important. An application of the standard MATHEMATICA function MatrixRank[X236 ] gives the result 5. However, since MatrixRank[R]=4, we conclude that MATHEMATICA is unable to compute the exact rank of X236 . On the other hand, we obtain r236 = 4., which is the desired value in accordance with Corollary 2.3. As a confirmation, we observe that MatrixRank[Pseudoinverse[A]]=4. Relative error estimate is kA† − X236 k/kA† k = 3.12183 × 10−6 . Example 3.3 Consider the following matrix: A=
1 −1 −1 −1 −1 −1
−1 1 −1 −1 −1 −1
0 0 1 −1 −1 0
0 0 −1 1 0 −1
0 0 0 0 2 −1
0 0 0 0 −1 2
.
Using R = A2 and β = 0.00137174, derived from (2.22), SMS method produces well known approximation with the precision δ = 10−6 of the Drazin inverse from [19] after k = 15 iterations: X215 =
0.25 −0.25 0. 0. 0. 0.
−0.25 0.25 0. 0. 0. 0.
0. 0. 0.25 −0.25 −0.416667 −0.583333
0. 0. −0.25 0.25 −0.583333 −0.416667
0. 0. 0. 0. 0.666667 0.333333
0. 0. 0. 0. 0.333333 0.666667
.
16
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
The exact Drazin inverse is equal to
A
D
=
1 4 − 14
0 0 0 0
− 14
0 0
0 0
0 0 0 0
1 4 − 14 5 − 12 7 − 12
− 14
1 4
1 4 7 − 12 5 − 12
0 0 0 0
0 0 0 0
2 3 1 3
1 3 2 3
.
The relative error kAD − X215 k/kAD k, where the matrix norm k · k represents the maximal singular value of the argument, is equal to 2.67633 × 10−14 . Let us mention that in [19] Euler-Knopp representation, corresponding to iterations (2.6) in the case R = Al , l ≥ ind(A), gives an analogous representation in 29697 iterations. Rank of the matrix R can be computed using r215 = 4. Example 3.4 In this example we verify SMS algorithm (2.14) in the case Re λi < 0, for each i = 1, . . . , t. For this purpose, we accomplished the SMS algorithm for the matrix −A, where A is defined as in Example 3.1. We used C βopt = −9.208888 × 10−6 defined by (2.23) and the number of iterative steps k = 23 conditioned by the second case of (2.28). Corresponding outer inverse, generated by the SMS algorithm is equal to −X223 , where X223 is the outer inverse of A generated in Example 3.1.
4
Conclusion
SMS algorithm is investigated in the papers [3], [16], [18]. In these articles SMS algorithm is based on successive squaring of the appropriate composite matrix T in the form (2.2). In the paper [3] using that P = I−βA∗ A, Q = βA∗ the MoorePenrose inverse is derived. The weighted Moore-Penrose inverse is derived in [18], using P = I − βA] A, Q = βA] , and finally the Drazin inverse is derived in [16] applying the SMS iterative process with the matrices P = I − βAl A, Q = βAl . In the present article we introduced and investigated SMS iterative scheme based on the matrices P = I − βRA and Q = βR. In this way, we derive an universal SMS algorithm containing all previous results as partial cases. Moreover, we prove that our method can be used to approximate various outer . Properties of the generated outer inverse inverses of a given matrix A ∈ Cm×n r X corresponding to chosen matrix R are examined. Moreover, in the case when a full-rank factorization R = F G is known, we give an explicit form of an outer generalized inverse of a A corresponding to the chosen R. An explicit approximation of rank(R) and rank(X) is also given. We implemented our algorithm in the programming package MATHEMATICA, and many numerical examples are presented in the last section. In these examples we show that sometimes MATHEMATICA is unable to compute the exact value of rank(X), so its iterative approximation is important.
Successive matrix squaring algorithm for computing outer inverses
17
References [1] A. Ben-Israel and T.N.E. Greville, Generalized inverses: theory and applications, Second Ed., Springer, 2003. [2] S.R.Caradus Generalized Inverses and Operator Theory, Queen’s Paper in Pure and Applied Mathematics, Queen’s University, Kingston, Ontario, 1978. [3] L. Chen, E.V. Krishnamurthy, I. Macleod, Generalized matrix inversion and rank compuation by successive matrix powering, Parallel Computing 20 (1994) 297–311. [4] D.S. Djordjevi´c and P.S. Stanimirovi´c, General representations of pseudoinverses, Matematiˇcki Vesnik 51 (1999), 69-76. [5] D.S. Djordjevi´c, P.S. Stanimirovic and Y. Wei, The representation and approximation of outer generalized inverses, Acta Math. Hungar, 104 (2004), 1–26. [6] R.E. Harte, Invertibility and singularity for bounded linear operator, Dekker 1988. [7] R.A. Horn and C.R. Johnson, Matrix Analysis, CambrodgrUniversity Press, Cambridge, New York, New Rochelle, Melbourne, Sydney, 1986. [8] C.R. Rao and S.K. Mitra, Generalized Inverse of Matrices and its Applications, John Wiley Sons, Inc, New York, London, Sydney, Toronto, 1971. ´ c, Adjoint mappings and inverses [9] P.S. Stanimirovi´c, S. Bogdanovi´c, M. Ciri´ of matrices, Algebra Colloquium 13(3) (2006) 421–432. [10] P.S. Stanimirovi´c, Block representations of {2}, {1,2} inverses and the Drazin inverse, Indian Journal Pure Appl. Math., 29 (1998) 1159–1176. [11] P.S. Stanimirovi´c, Applications of hyper-power method for computing matrix products, Publikacije Elektrotehniˇckog fakulteta, 15 (2004) 13–25. [12] P.S. Stanimirovi´c and D.S. Djordjevi´c, Full-rank and determinantal representation of the Drazin inverse, Linear Algebra Appl., 311 (2000) 31–51. [13] G. Wang, Y. Wei, S. Qiao, Generalized inverses: theory and computations, Science Press, 2003. (1,2)
[14] G. Wang, The representations of the generalized inverses (A ⊗ B)T,S and (2)
(A⊗ B)T,S and some applications, J. Shanghau Univ. (Natural Sciences) 24 (1995) 1–6. [15] Y. Wei, A characterization and representation for the generalized inverse (2) AT,S and its applications, Linear Algebra Appl., 280 (1998), 87–96.
18
P.S. Stanimirovi´c and D.S. Cvetkovi´c-Ili´c
[16] Y. Wei, Successive matrix squaring algorithm for computing Drazin inverse, Appl. Math. Comp. 108 (2000), 67–75. [17] Y. Wei and H. Wu, The representation and approximation for the general(2) ized inverse AT,S , Appl. Math. Comput., 135 (2003), 263–276. [18] Y. Wei, H. Wu, J. Wei Successive matrix squaring algorithm for parallel computing the weighted generalized inverse A†M N , Appl. Math. Comp. 116 (2000) 289–296. [19] Y. Wei and H. Wu, The representation and approximation for Drazin inverse, J. Comput. Appl. Math. 126 (2000), 417–423. [20] S. Wolfram, The Mathematica Book, 4th ed., Wolfram Media/Cambridge University Press, 1999.
