Near Inversion of Generalized Sylvester Matrices - University of Waterloo
Computation of Numerical Pad´ e-Hermite and Simultaneous Pad´ e Systems I: Near Inversion of Generalized Sylvester Matrices Stan Cabay∗, Anthony R. Jones†, and George Labahn‡ Abstract. We present new formulae for the “near” inverses of striped Sylvester and mosaic Sylvester matrices. The formulae assume computation over floating-point rather than exact arithmetic domains. The near inverses are expressed in terms of numerical Pad´ e-Hermite systems and simultaneous Pad´ e systems. These systems are approximants for the power series determined from the coefficients of the Sylvester matrices. The inverse formulae provide good estimates for the condition numbers of these matrices, and serve as primary tools in a companion paper for the development of a fast, weakly stable algorithm for the computation of Pad´ e-Hermite and simultaneous Pad´ e systems and, thereby, also for the numerical inversion of striped and mosaic Sylvester matrices. Key Words. striped Sylvester inverses, mosaic Sylvester inverses, Pad´ e-Hermite approximants, simultaneous Pad´ e approximants, numerical stability AMS(MOS) Subject Classification. 41A21, 65F05, 65G05
1. Introduction. Let n = [n0 , . . . , nk ] be a vector of integers with nβ ≥ 0, 0 ≤ β ≤ k. A striped Sylvester matrix of order knk is given by (0) (0) a0 ak .. .. . . . . , (0) (0) . . Mn = (1) . a0 ··· . ak .. .. . . (knk−1) (knk−n0 ) (knk−1) (knk−nk ) a0 · · · a0 ak · · · ak (`)
where knk = n0 + . . . + nk and where the coefficients aβ
∈ R, the field of real
(0) a0
numbers. Assume that 6= 0. In this paper, we present a formula for the inverse of Mn expressed in terms of Pad´e-Hermite and simultaneous Pad´e systems. Pad´e-Hermite and simultaneous Pad´e systems [7, 9] are approximants for the vector At (z) = [a0 (z), . . . , ak (z)] of power series associated with the coefficients of Mn , namely, knk−1
aβ (z) =
X
(`)
aβ z ` , with 0 ≤ β ≤ k.
`=0
They provide necessary and sufficient conditions for Mn to be nonsingular, and generalize the notions of Pad´e-Hermite and simultaneous Pad´e approximants. Pad´eHermite and simultaneous Pad´e approximants were introduced by Hermite in the last ∗ Department of Computing Science, University of Alberta, Edmonton, Alberta, Canada, T6G 2H1. The research of this author as partially supported by Natural Sciences and Engineering Research Council of Canada grant A8035. † Bell Northern Research, P.O. Box 3511, Station C, Ottawa, Ontario, Canada, K1Y 4H7 ‡ Department of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, N2L3G1. The research of this author was partially supported by Natural Sciences and Engineering Research Council of Canada grant FS1525C.
1
century and has been widely studied by several authors (for a bibliography, see, for example [1, 2, 3, 4, 15]). The inverse formula given here is intended to be applied in a numerical setting; it accommodates errors that may have been made in the computation of Pad´e-Hermite and simultaneous Pad´e systems. That is, the formula gives the “near” inverse for Mn since it expresses the inverse in terms of Pad´e-Hermite and simultaneous Pad´e systems which are computed using floating-point arithmetic. There are other closed formulae (c.f. [12, 16, 18, 19, 20]) for M−1 n . The formula given here is different in that it expresses the inverse directly in terms of numerical Pad´e-Hermite and simultaneous Pad´e systems. The near inverse formula depends on the computation of both Pad´e systems. It is possible to determine a simultaneous Pad´e system from its “dual” Pad´e-Hermite system via the adjoint operation [6, 15]. In a numerical setting, however, this does not provide enough control over the resulting floating-point errors [14]. Instead, simultaneous Pad´e systems can be computed independently. Whereas a Pad´e-Hermite system can be obtained by solving a system of linear equations with M n as the coefficient matrix, a simultaneous Pad´e system can be similarly and independently obtained with a coefficient matrix that is a specialized mosaic Sylvester matrix. This specialized mosaic Sylvester matrix of order kknk is defined to be ∗ ∗ S0,1 · · · S0,k .. , M∗n = ... (2) . ∗ ∗ Sk,1 · · · Sk,k ∗ where Sα,β are matrices of size (knk − nα ) × knk, with (0) (knk−1) aβ ··· aβ .. ∗ .. S0,β = − . . (0) (n0 ) aβ ··· aβ
∗ Sβ,β =
(0)
a0
(knk−1)
··· ..
a0
.. .
. (0)
a0
···
(nβ )
a0
,
,
∗ = 0. The matrix M∗n is closely related to for 1 ≤ β ≤ k, and with the remaining Sα,β Mn . Indeed, in the special case when k = 1 the matrix Mn and the transpose of M∗n coincide, up to a sign and a permutation of the rows and columns. Similarly, when k ≥ 1 and a0 (z) = 1, the matrix Mn and the transpose of M∗n are both obtained by a suitable block extension of the same matrix. In this paper, we also provide a “near” inverse formula for the matrix M∗n , again in terms of numerical Pad´e-Hermite and simultaneous Pad´e systems. The inverse formulae provide “good” estimates for the condition numbers of M n and M∗n and these enable the stable numerical computation of Pad´e-Hermite and simultaneous Pad´e systems described in the companion paper [6]. In return, the accurate computation of Pad´e-Hermite and simultaneous Pad´e systems by the algorithm in [6] enables the effective inversion of generalized Sylvester matrices by the formulae given in this paper, as well as the solution of linear systems of equations with these as the coefficient matrices.
2
This paper is organized as follows. Preliminary definitions and basic facts about Pad´e-Hermite and simultaneous Pad´e systems are given in the next two sections. §3 gives a near commutativity relationship between these two systems in floating-point arithmetic, while §4 and §5 give the approximate inversion formulae for striped and mosaic Sylvester matrices. The final section gives a summary and some conclusions. We conclude this section by defining some norms which are used in §4 and §5. Let a(z) =
∞ X
a(`) z ` ∈ R [[z]] ,
`=0
where R[[z]] is the domain of power series with coefficients from R, and define the bounded power series ¯ ) ( ∞ ¯ X ¯ (`) B |a | < ∞ . R [[z]] = a(z) ¯ a(z) ∈ R [[z]] , ¯ `=0
A norm for a(z) ∈ RB [[z]] is
ka(z)k =
∞ X
|a(`) |.
`=0
RB [[z]] includes the domain of polynomials R[z]. So, for s(z) =
∂ X
s(`) z ` ∈ R [z] ,
`=0
we use the norm ks(z)k =
∂ X
| s(`) | .
`=0
For vectors and matrices over RB [[z]], we use the 1-norm unless otherwise specified. So, for example, the norm for At (z) is kAt (z)k = max {kaβ (z)k} 0≤β≤k
and the norm for S(z) ∈ R(k+1)×(k+1) [z] is kS(z)k = max
0≤β≤k
(
k X
α=0
)
kSα,β (z)k .
It is easy to verify that various compatibility conditions are satisfied. For example, kAt (z) · S(z)k ≤ kAt (z)k · kS(z)k and ka(z) · b(z)k ≤ ka(z)k · kb(z)k, 3
where b(z) is also a bounded power series. In addition, for S ∗ (z) ∈ R(k+1)×(k+1) [z] and A∗ (z) ∈ RB (k+1)×k [[z]], kS ∗ (z) · A∗ (z)k ≤ kS ∗ (z)k · kA∗ (z)k, kS(z) · S ∗ (z)k ≤ kS(z)k · kS ∗ (z)k. In the subsequent development, we also make use of the inequality ka(z) (mod z knk+1 )k ≤ ka(z)k, where a(z) (mod z knk+1 ) =
knk X
a(`) z ` +
`=0
∞ X
0 · z ` ∈ RB [[z]]
`=knk+1
2. Multi-dimensional Pad´ e Systems. In this section, we introduce the notion of Pad´e-Hermite and simultaneous Pad´e systems. Let n = [n0 , . . . , nk ] and define knk = n0 + · · · + nk . Also let At (z) = [a0 (z), . . . , ak (z)] , where aβ (z) =
∞ X
(`)
aβ z ` ,
β = 0, . . . , k,
`=0
(0)
(`)
with aβ ∈ R, the field of real numbers. Assume that a0 6= 0, which means that t knk+1 )k = 1, a−1 0 (z) exists. Assume also that A (z) is scaled so that kaβ (z) (mod z 0 ≤ β ≤ k. The (k + 1) × (k + 1) matrix of polynomials zp(z) u1 (z) · · · uk (z) · ¸ zq1 (z) v1,1 (z) · · · v1,k (z) zp(z) U t (z) S(z) = (3) = .. .. .. zQ(z) V (z) . . . zqk (z) vk,1 (z) · · · vk,k (z) is a numerical Pad´e-Hermite system (NPHS) [9] of type n for A(z) if the following conditions are satisfied. I. (Degree conditions): For 1 ≤ α, β ≤ k, (4)
p(z) =
nX 0 −1
p(`) z ` ,
uβ (z) =
`=0
qα (z) =
nX α −1
n0 X `=0
qα(`) z ` ,
vα,β (z) =
`=0
nα X `=0
II. (Order condition): (5)
(`)
uβ z ` ,
At (z)S(z) = z knk T t (z) + δT t (z), 4
(`)
vα,β z ` .
where T t (z) = [r(z), zW t (z)] with W t (z) = [w1 (z), . . . , wk (z)] is the residual, and where δT t (z) = [z δr(z), δW t (z)] is the residual error, with δW t (z) = [δw1 (z), . . . , δwk (z)] and with knk−2
δr(z) =
X
δr
(`)
`
z , δwβ (z) =
`=0
knk X
(`)
δwβ z ` .
`=0
If δT t (z) = 0, then S(z) is an exact, rather than a numerical, Pad´e-Hermite system. III. (Nonsingularity condition): The constant term of V (z) is a diagonal matrix, (6)
V (0) = diag [γ1 , . . . , γk ] ,
and (7)
(0)
γ ≡ (a0 )−1
k Y
γα 6= 0,
α=0
where γ0 = r(0). Remark 1: The nonsingularity condition III is equivalent to the condition that r(0) 6= 0 and that V (0) be a nonsingular diagonal matrix. Remark 2: The NPHS is said to be normalized [9] if the nonsingularity condition III is replaced by r(0) = 1 and V (0) = Ik . This can be achieved by multiplying S(z) on the right by Γ−1 , where (8)
Γ = diag [γ0 , . . . , γk ] .
The NPHS is said to be scaled [14] if each column of S(z) has norm equal to 1 and if, in addition, γβ > 0, 0 ≤ β ≤ k. Here, also, scaling a NPHS is accomplished by multiplying it on the right by an appropriate diagonal matrix. Remark 3: The nonsingularity condition III, namely γ 6= 0, refers to the nonsingularity of the associated striped Sylvester matrix Mn defined in (1); in [9] it is shown that an exact NPHS exists iff Mn is nonsingular. Remark 4: In [5, 6, 9], rather than S(z), the Pad´e-Hermite system is defined to be S(z) · diag[z, 1, . . . , 1]. The notation used here is consistent with that of [16], and simplifies the presentation of some of the results. A Pad´e-Hermite system gives an approximation to a vector of formal power series using matrix multiplication on the right. A simultaneous Pad´e system is a similar approximation using matrix multiplication on the left and with degree constraints that can be thought of as being “dual” to the degree constraints of a Pad´e-Hermite system.
5
Let1
−a1 (z) a0 (z) A∗ (z) = 0
(9)
−ak (z) 0 a0 (z)
··· ..
.
be a (k + 1) × k matrix of power series. The (k + 1) × (k + 1) matrix of polynomials ∗ ··· u∗k (z) u∗1 (z) v (z) · ∗ ¸ zq ∗ (z) zp∗ (z) · · · zp∗ (z) 1,1 1,k v (z) U ∗t (z) 1 S ∗ (z) = (10) = . . .. ∗ ∗ . . zQ (z) zP (z) . . . zqk∗ (z) zp∗k,1 (z) · · · zp∗k,k (z)
is a numerical simultaneous Pad´e system (NSPS) [7, 9] of type n for A∗ (z) if the following conditions are satisfied. I. (Degree conditions): For 1 ≤ α, β ≤ k, knk−nβ
knk−n0
(11)
v ∗ (z) =
v ∗ (`) z ` ,
X
uβ z ` ,
`=0
`=0
knk−nβ −1
knk−n0 −1
qα∗ (z) =
∗(`)
X
u∗β (z) =
X
qα∗(`) z ` ,
X
p∗α,β (z) =
∗(`)
pα,β z ` .
`=0
`=0
II. (Order condition): S ∗ (z)A∗ (z) = z knk T ∗ (z) + δT ∗ (z),
(12)
where T ∗t (z) = [z W ∗ (z)|R∗t (z)] with R∗ (z) a k × k is the residual, and where δT ∗t (z) = [δW ∗ (z)|z δR∗t (z)] is the residual error, with δR∗ (z) a k × k matrix and δwβ∗ (z) =
knk X
knk−2 ∗(`)
δwβ
∗ z ` , δrα,β (z) =
X
∗(`)
δrα,β z ` .
`=0
`=0
If δT ∗ (z) = 0, then S ∗ (z) is an exact NSPS. III. (Nonsingularity condition): The constant term of R ∗ (z) is a diagonal matrix R∗ (0) = diag [γ1∗ , . . . , γk∗ ] ,
(13) and
(0)
γ ∗ ≡ (a0 )−1
(14)
k Y
γα∗ 6= 0,
α=0 1
More generally,
A∗ (z) =
a∗0,1 (z) a∗1,1 (z) . . . ∗ ak,1 (z)
··· ···
···
a∗0,k (z) a∗1,k (z) . . . ∗ ak,k (z)
h =
B ∗t (z) C ∗ (z)
i
with C ∗ (0) nonsingular, but the restriction to (9), which algebraically is made without loss of generality, permits us to establish in §3 a duality relationship between Pad´ e-Hermite systems and simultaneous Pad´ e systems. 6
where γ0∗ = v ∗ (0). Remark 5: The NSPS is said to be normalized [7] if the nonsingularity condition III is replaced by v ∗ (0) = 1 and R∗ (0) = Ik . This can be achieved by multiplying S ∗ (z) on the left by Γ∗−1 , where Γ∗ = diag [γ0∗ , . . . , γk∗ ] .
(15)
The NSPS is said to be scaled when each row of S ∗ (z) has norm equal to 1 and if, in addition, γα∗ > 0, 0 ≤ α ≤ k. Here, also, scaling a NSPS is accomplished by multiplying it on the left by an appropriate diagonal matrix. Remark 6: The nonsingularity condition III, namely γ ∗ 6= 0, refers to the nonsingularity of the associated mosaic Sylvester matrix M∗n defined in (2); in [7] it is shown that an exact NSPS exists iff M∗n is nonsingular. Remark 7: In [5, 6, 9], rather than S ∗ (z), the simultaneous Pad´e system is defined to be diag[1, z, . . . , z] · S ∗ (z). This is for notational convenience only. 3. Duality. Theorem 1 below gives a relationship between Pad´e-Hermite and simultaneous Pad´e systems which is crucial to the results of the subsequent sections. It generalizes earlier results of Mahler and their extensions to block matrices [10, 15, 17, 21]. Theorem 1. If S(z) is a NPHS of type n for A(z) and S ∗ (z) is a NSPS of type n for A∗ (z), then (0)
S ∗ (z) · S(z) = z knk (a0 )−1 Γ∗ Γ + θI (z),
(16) where θI (z)
=
a−1 0 (z)
½·
v ∗ (z) zQ∗ (z)
¸
∗
t
δT (z) + δT (z)
£
zQ(z)
V (z)
¤
¾
(mod z knk+1 ).
Proof. The theorem (in the case that δT (z) = 0 and δT ∗ (z) = 0) follows from [15]. The arguments used in the following proof, however, are considerably simpler. Let B t (z) = [a1 (z), · · · , ak (z)]. Then, using (5) and (12), (17)
a0 (z) S ∗ (z) · S(z) ½· ∗ ¸ · ∗t ¸ ¾ £ ¤ ¤ v (z) U (z) £ t zp(z) U (z) zQ(z) V (z) = a0 (z) + zQ∗ (z) zP ∗ (z) · ∗ ¸ ¤ £ ¤ª © £ v (z) = a0 (z) zp(z) U t (z) + B t (z) zQ(z) V (z) ∗ zQ (z) ½ · ∗t ¸ · ∗ ¸ ¾ £ ¤ U (z) v (z) t zQ(z) V (z) + a0 (z) − B (z) ∗ ∗ zP (z) zQ (z) 7
= =
·
v ∗ (z) zQ∗ (z) ½·
¸
£ ¤ At (z)S(z) + S ∗ (z)A∗ (z) zQ(z) V (z) ¸ ¸ · £ ¤ v ∗ (z) zW ∗t (z) £ t r(z) W (z) zQ(z) z knk + zQ∗ (z) R∗ (z) ¸ · ∗ £ ¤ v (z) δT t (z) + δT ∗ (z) zQ(z) V (z) . + ∗ zQ (z)
V (z)
¤
¾
But, from (4) and (11), the degrees of S ∗ (z)S(z) are bounded component-wise by knk. It then follows from (17) that ¸ · ∗ 0 v (0)r(0) ∗ knk (0) −1 + θI (z) S (z)S(z) = z (a0 ) 0 R∗ (0)V (0) =
(0)
z knk (a0 )−1 Γ∗ Γ + θI (z),
which is (16). Corollary 2. If S(z) is a normalized NPHS of type n for A(z) and S ∗ (z) is a normalized NSPS of type n for A∗ (z), then for sufficiently small2 δT (z) and δT ∗ (z) (0)
S(z) · S ∗ (z) = z knk (a0 )−1 Ik+1 + θII (z),
(18) where
(0)
θII (z) = S(z) θI (z) [z knk (a0 )−1 Ik+1 + θI (z)]−1 S ∗ (z). Proof. If θI (z) is so small, that is, if δT (z) and δT ∗ (z) are so small, that (0) z knk (a0 )−1 Ik+1 + θI (z) is nonsingular, then from (16) (0)
S ∗−1 (z) = S(z) · [z knk (a0 )−1 Ik+1 + θI (z)]−1 . So, (0)
S(z)S ∗ (z) − z knk (a0 )−1 Ik+1 n o (0) = S(z) − z knk (a0 )−1 S ∗−1 (z) S ∗ (z) n o (0) (0) = S(z) Ik+1 − z knk (a0 )−1 [z knk (a0 )−1 Ik+1 + θI (z)]−1 S ∗ (z) =
(0)
S(z) θI (z) [z knk (a0 )−1 Ik+1 + θI (z)]−1 S ∗ (z).
Corollary 3. The residuals T (z) for a normalized NPHS of type n for A(z) and T ∗ (z) for a normalized NSPS of type n for A∗ (z) satisfy (0)
t T t (z) S ∗ (z) = (a0 )−1 At (z) + θIII (z),
(19) where
2
© ª t θIII (z) = At (z)θII (z) − δT t (z)S ∗ (z)) /z knk .
It is adequate, for example, that condition (34) of Corollary 6 be satisfied. 8
Proof. From (5) and (18), it follows that n o z knk T t (z) + δT t (z) S ∗ (z) = At (z) S(z) S ∗ (z) n o (0) = At (z) z knk (a0 )−1 Ik+1 + θII (z)
and so (19) is true.
4. The Inverse of a Striped Sylvester Matrix. In this section, a formula is given for the inverse of Mn expressed in terms of both S(z) and S ∗ (z). This enables estimating the condition number of Mn without explicitly computing M−1 n . Associated with the NPHS S(z), define the order knk matrices (0) (n −1) (0) (n −1) q1 · · · q1 1 qk · · · qk k p(0) · · · p(n0 −1) .. .. .. . . . .. .. .. 0 0 . . 0 . (n1 −1) (nk −1) .. .. p(n0 −1) . . . . . q · · · q 1 k P= . . . . . . 0 . 0 . 0 . .. .. .. . . . 0
0 ··· 0 (20) and, for β = 1, 2, . . . , k, (1) (n ) u · · · uβ 0 β. . . .. 0 . (n0 ) . . . u (21) Uβ = β .. 0 . .. . 0 ··· 0
(1)
v1,β .. . (n )
v1,β1
···
0
··· . .. . ..
v1,β1
0 P∞
(n )
···
(1)
vk,β .. .
0 ··· .. .
0 .. .
0
0
(n )
vk,βk
···
··· . .. . ..
0 .. .
0 .. . 0
(n )
vk,βk
···
0
0
.
(`) ` Also, for any power£ series a(z) ¤ = `=0 a z , and any integer function f (i, j), (f (i,j)) i, j = 1, 2, . . . , let a ] denote a matrix of order knk whose element in position (i, j) is a(f (i,j)) . The main result of this section is Theorem 4 below which gives the inverse of M n in terms of the NPHS S(z) and the NSPS S ∗ (x) of types n for A(z).
Theorem 4. In terms of the normalized NPHS S(z) and the normalized NSPS S ∗ (x) of types n for A(z), the inverse of Mn satisfies k h o nh i i X h i (0) (i−j) ∗(knk−i−j) P t v ∗(knk−i−j+1) + a0 + θIV = a0 (22) M−1 Uβt qβ , n β=1
where
θIV
=
(0) a0
( h
(i−j)
(θIII )0
i
−
k h i ih X (i−j) (knk+i−j) (θII )α,0 aα
α=0
9
h
+ δr(i+j−2)
ih
v ∗(knk−i−j+1)
i
k h ih i X (i+j−1) ∗(knk−i−j) . δwβ + qβ β=1
Proof. The coefficient of z i+j−2 , for i, j = 1, 2, . . . , knk, in the first component of (5), namely, a0 (z) p(z) +
k X
aα (z) qα (z) = z knk−1 r(z) + δr(z),
α=1
is n0 X
(i+j−`−2) (`)
a0
p
+
k nX α −1 X
α=1 `=0
`=0
This is the (i, j)th component of £ (−knk+i+j−1) ¤ £ (i+j−2) ¤ (23)
(i+j−`−2) (`) qα = r(−knk+i+j−1) + δr (i+j−2) . aα
+ δr
r
h
=
(knk+i−j)
a0
k X £
+
i£
p(−knk+i+j−2)
a(knk+i−j) α
α=1
¤£
¤
qα(−knk+i+j−2) + Mn P t .
¤
Similarly, the coefficient of z i+j−1 , for i, j = 1, 2, . . . , knk, in the (β + 1)st component, β = 1, . . . , k, of (5), namely, a0 (z) uβ (z) +
k X
aα (z) vα,β (z) = z knk+1 wβ (z) + δwβ (z),
α=1
is n0 X
(i+j−`−1) (`) uβ
a0
+
nα k X X
(i+j−1)
+ δwβ
.
α=1 `=0
`=0
This is the (i, j)th component of i i h h (24)
(−knk+i+j−2)
(`)
(i+j−`−1) aα vα,β = wβ
(−knk+i+j−2)
wβ
(i+j−1)
+ δwβ
= +
h
(knk+i−j)
a0
k X £
ih
a(knk+i−j) α
α=1
(−knk+i+j−1)
uβ
¤h
i
(−knk+i+j−1)
vα,β
i
+ Mn Uβt .
Next, the coefficient of z i−j−1 for i, j = 1, . . . , knk, in the first row and first column of (18) for a normalized NPHS and a normalized NSPS, namely, p(z) v ∗ (z) +
k X
(0)
uβ (z)qβ∗ (z) = z knk−1 (a0 )−1 + z −1 (θII )0,0 (z),
β=1
is nX 0 −1 `=0
v
∗(i−j−`−1)
p
(l)
+
n0 k X X
β=1 `=0
10
∗(i−j−`−1) (`) uβ
qβ
(i−j)
= (θII )0,0 .
This is the (i, j)th component of (25)
i ih h p(−knk+i+j−2) v ∗(knk−i−j+1)
k h X
+
ih
(−knk+i+j−1)
uβ
β=1
h
=
(i−j)
(θII )0,0
i
∗(knk−i−j)
qβ
i
.
The coefficient of z i−j−1 in the first column and the (α + 1)st row, α = 1, . . . , k, of (18), namely, qα (z) v ∗ (z) +
k X
vα,β (z)qβ∗ (z) = z −1 (θII )α,0 (z)
β=1
is nα X
v ∗(i−j−`−1) qα(l) +
nα k X X
∗(i−j−`−1) (`) vα,β
qβ
(i−j)
= (θII )α,0 .
β=1 `=0
`=0
This is the (i, j)th component of (26)
h
qα(−knk+i+j−2)
ih
v ∗(knk−i−j+1)
i
+
k h X
(−knk+i+j−1)
vα,β
β=1
=
h
(i−j)
(θII )α,0
i
ih
∗(knk−i−j)
qβ
i
.
Also, the coefficient of z i−j , for i, j = 1, . . . , knk in the first component of (19) for a normalized NPHS and NSPS, namely, ∗
r(z)v (z) + z
2
k X
(0)
wβ (z)qβ∗ (z) = (a0 )−1 a0 (z) + (θIII )0 (z).
β=1
is the (i, j)th component of i i h h (i−j) (0) (i−j) (27) = + (θIII )0 (a0 )−1 a0 +
h
r(−knk+i+j−1)
k h X
ih
v ∗(knk−i−j+1)
(−knk+i+j−2)
wβ
β=1
ih
i
∗(knk−i−j)
qβ
i
.
We are finally ready to prove the theorem. ¿From (23), (24), (25), (26) and (27),
h
i
k X
h
∗(knk−i−j)
i
P t v ∗(knk−i−j+1) + Uβt qβ β=1 i i h ih i h nh (knk+i−j) p(−knk+i+j−2) = r(−knk+i+j−1) + δr(i+j−2) − a0
Mn
−
k h X
α=1
(knk+i−j) aα
ih
qα(−knk+i+j−2) 11
io h
v ∗(knk−i−j+1)
i
+
k nh X
(−knk+i+j−2)
wβ
β=1
−
k h X
(knk+i−j) aα
α=1
=
ih
i
h i h ih i (i+j−1) (knk+i−j) (−knk+i+j−1) + δwβ − a0 uβ
(−knk+i+j−1)
vα,β
io h
∗(knk−i−j)
qβ
i
k h h ih i X ih i (−knk+i+j−2) ∗(knk−i−j) r(−knk+i+j−1) v ∗(knk−i−j+1) + wβ qβ β=1
+
k h h ih i X ih i (i+j−1) ∗(knk−i−j) δwβ δr(i+j−2) v ∗(knk−i−j+1) + qβ β=1
−
k X
α=0
=
h
(knk+i−j) aα
ih
(i−j)
(θII )α,0
h i (0) (i−j) (a0 )−1 a0 + θIV .
i
The result (22) now follows. Corollary 5 below drops the requirement in Theorem 4 that S(z) and S ∗ (z) be normalized. In particular, the corollary is valid in the case that S(z) and S ∗ (z) are scaled. Corollary 5. In terms of the (unnormalized) NPHS S(z) of type n for A(z) and the (unnormalized) NSPS S ∗ (z) of type n for A∗ (z), the inverse of Mn is given by o nh i (i−j) (28) M−1 a0 + θ¨IV n k h i X i h ∗(knk−i−j) (0) (γβ γβ∗ )−1 Uβt qβ , (γ0 γ0∗ )−1 P t v ∗(knk−i−j+1) + = a0 β=1
where
(29)
θ¨IV
=
(0) a0
+ +
(
(γ0 γ0∗ )−1 k X
β=1
(30)
t θ¨III (z)
=
(31)
θ¨II (z)
=
(32)
θ¨I (z)
=
n
k h h i X ih i (i−j) (i−j+1) (knk+i−j) ¨ (θIII )0 − aα (θ¨II )α,0 α=0
h
δr
(i+j−2)
h
ih
v ∗(knk−i−j+1)
(i+j−1)
(γβ γβ∗ )−1 δwβ
ih
i
∗(knk−i−j)
qβ
i
,
o At (z)θ¨II (z) − δT t (z)(Γ∗ Γ)−1 S ∗ (z)) /z knk ,
S(z) (Γ∗ Γ)−1 θ¨I (z) (Γ∗ Γ)−1 (0) ·[z knk (a0 )−1 Ik+1 + θ¨I (z)(Γ∗ Γ)−1 ]−1 S ∗ (z), ½· ∗ ¸ v (z) a−1 (z) δT t (z) 0 zQ∗ (z) £ ¤ª + δT ∗ (z) zQ(z) V (z) (mod z knk+1 ) 12
Proof. The normalized NPHS is obtained from an unnormalized one by multiplying it on the right by the diagonal matrix diag[γ0−1 , . . . , γk−1 ]. Similarly, the normalized NSPS is obtained from an unnormalized one by multiplying it on the left by the diagonal matrix diag[γ0∗−1 , . . . , γk∗−1 ]. The result now follows directly from (22). Let (33)
κ=
k X
(γβ γβ∗ )−1 .
β=0
In the corollary below, we give a bound for M−1 n in terms of κ. Corollary 6. If the residual errors δT t (z) and δT ∗ (z) associated with the scaled S(z) and the scaled S ∗ (z) are not too large, so that h
(34)
then
kM−1 n k1
(35)
(0)
knk+1 (κ + 1)(k + 2)|a0 |(ka−1 )k + 1) 0 (z) (mod z £ ¤ · (k + 2)kδT t (z)k + kδT ∗ (z)k ≤ 1/8,
≤
i2
(0)
knk+1 2κ · |a0 | · ka−1 )k. 0 (z) (mod z
Proof. We use Corollary 5 with S(z) and S ∗ (z) scaled. We begin by finding a bound for θ¨IV appearing in the inverse formula (28) for Mn . A bound for θ¨IV depends on bounds for θ¨I (z), θ¨II (z) and θ¨III (z). ¿From (16), © ª knk+1 kθ¨I (z)k ≤ ka−1 )k · (k + 1)kδT t (z)k + kδT ∗ (z)k , 0 (z) (mod z
since kS(z)k = 1 and kS ∗ (z)k ≤ k + 1. ¿From (16) and (32), note that θ¨I (z) is a matrix polynomial of at most degree knk and so, using (34), (0) (0) (0) ka0 z knk θ¨I (z −1 )(Γ∗ Γ)−1 k = ka0 θ¨I (z)(Γ∗ Γ)−1 k ≤ κ · |a0 | · kθ¨I (z)k ≤ 1/2, (0)
since k(Γ∗ Γ)−1 k ≤ κ. So as in Stewart [22, page 187], the inverse of (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 exists and o−1 n (0) k (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 k
(0)
≤
|a0 | (0) knk ¨ 1 − ka z θI (z −1 )(Γ∗ Γ)−1 k
≤
2|a0 |.
0
(0)
To determine a bound for θ¨II (z) in (31), let N = max0≤β≤k {nβ } and observe from 18 that θ¨II (z) is also a matrix polynomial, now of degree at most knk + N . Consequently, (36)
kθ¨II (z)k
kz knk+N θ¨II (z −1 )k n o © ª = k z N S(z −1 ) (Γ∗ Γ)−1 z knk θ¨I (z −1 ) (Γ∗ Γ)−1
=
13
≤
≤ ≤
o−1 n (0) · (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 [z knk S ∗ (z −1 )]k
κ2 kz N S(z −1 )k · kz knk θ¨I (z −1 )k · kz knk S ∗ (z −1 )k n o−1 (0) ·k (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 k (0) 2κ2 |a0 | · kS(z)k · kθ¨I (z)k · kS ∗ (z)k (0) 2κ2 (k + 1) · |a | · kθ¨I (z)k. 0
In addition, it now follows that a bound for θ¨III (z) in (30) is given by (0) t (z)k ≤ 2κ2 (k + 1) · |a0 | · kθ¨I (z)k + κ(k + 1) · kδT t (z)k. kθ¨III
We are now ready to give a bound for θ¨IV appearing in the inverse formula (28). But, first observe that h i (i−j) k (θ¨III ) k1 ≤ kθ¨t (z)k 0
III
and that k
k h X
(knk+i−j) aα
α=0
ih
(i+j)
(θ¨II )α,0
i
k1 ≤
k X
k(θ¨II )α,0 (z)k ≤ k(θ¨II )(z)k.
α=0
Thus, kθ¨IV k1
=
(0) ka0
+
(
(γ0 γ0∗ )−1
≤ ≤ ≤
≤
α=0
h
δr
(i+j−2)
ih
v ∗(knk−i−j+1)
i
h ih i (i+j−1) ∗(knk−i−j) + (γβ γβ∗ )−1 δwβ qβ k1 β=1 k X (0) t (γβ γβ∗ )−1 kδT t (z)k (z)k + kθ¨II (z)k + (γ0 γ0∗ )−1 kδT t (z)k + |a0 | · kθ¨III β=1 n o (0) t |a0 | kθ¨III (z)k + kθ¨II (z)k + κkδT t (z)k n o (0) (0) |a0 | κ(k + 1)kδT t (z)k + 4κ2 (k + 1)|a0 | · kθ¨I (z)k + κkδT t (z)k (0) © κ(k + 2) · |a0 | kδT t (z)k £ ¤o (0) knk+1 t ∗ +4κ|a0 | · ka−1 (z) (mod z )k (k + 1)kδT (z)k + kδT (z)k 0 h i2 £ ¤ (0) knk+1 4ka−1 (z) (mod z )k (κ + 1)(k + 2)|a | (k + 2)kδT t (z)k + kδT ∗ (z)k . 0 0 k X
≤
k h i ih i X h (i−j) (i−j) (knk+i−j) aα (θ¨II )α,0 − (θ¨III )0
It then follows from (34) that
h i−1 (i−j) knk+1 k a0 θ¨IV k1 ≤ ka−1 )k · kθ¨IV k1 ≤ 1/2, 0 (z) (mod z 14
i−1 h (i−j) and so Iknk + a0 θ¨IV is invertible. In addition, k
nh
(i−j)
a0
i
+ θ¨IV
o−1
k1
¾−1 h ½ i−1 i−1 h (i−j) (i−j) a0 k1 k Iknk + a0 θ¨IV
≤
h i−1 (i−j) k a0 k1 i−1 h (i−j) θ¨IV k1 1 − k a0 i−1 h (i−j) 2k a0 k1
≤
≤
knk+1 2ka−1 )k. 0 (z) (mod z
≤
Therefore, a bound for M−1 n in (28) is given kM−1 n k1
o−1
h i (γ0 γ0∗ )−1 Pnt v ∗(knk−i−j+1) k h i X ∗(knk−i−j) t qβ (γβ γβ∗ )−1 Un,β + k 1
≤
k
nh
(i−j)
a0
i
+ θ¨IV
(0)
k · ka0
n
β=1
(0)
knk+1 2κ|a0 | · ka−1 )k. 0 (z) (mod z
≤
¿From (35), it follows that a bound for the 1-norm condition number of Mn is (0)
−1 knk+1 kMn k1 · kM−1 )k, n k1 ≤ 2κ|a0 | · ka0 (z) (mod z
since it is assumed that each aβ (z) is scaled. 5. The Inverse of a Mosaic Sylvester Matrix. In this section, a formula is given for the inverse of M∗n expressed in terms of both S(z) and S ∗ (z). This enables estimating the condition number of M∗n without explicitly computing M∗−1 n . Associated with the NPHS S(z) and the NSPS S ∗ (z), for β = 1, 2, . . . , k, define the knk × kknk matrices (knk−1) (0) (knk−1) (0) v1,β · · · v1,β vk,β · · · vk,β .. .. . . , Vβ = .. .. ··· . . (0) (0) vk,β v1,β
Q=
(knk−2)
q1
.. . (0)
q1 0
··· . .. . ..
(0)
q1
(knk−2)
0
qk ···
.. . 0)
qk 0
15
··· . .. . ..
(0)
qk
0
,
v ∗(1) .. . v ∗(η0 ) V∗ = 0 .. . 0
and
∗(0)
qβ .. .
∗(η −1) q 0 ∗ β Qβ = 0 .. . 0
v ∗(η0 )
··· . .. ..
0
.
··· . .. ..
∗(η1 )
u1
.. . ···
··· . .. . ..
∗(1)
u1 .. .
∗(η0 −1)
0
···
∗(0)
pβ,1 .. .
∗(η −1)
pβ,11 .. .
···
..
∗(ηk )
∗(ηk )
uk
0
. .. .
0 .. .
···
0
··· . .. . ..
pβ,11
···
··· . ..
uk
0
∗(η −1)
···
∗(0)
pβ,k .. .
0 ··· .. .
0 .. . 0
0
0 .. .
0
qβ
∗(1)
uk .. .
.
0 .. .
0
∗(η1 )
u1
∗(η −1)
pβ,kk
··· . .. . ..
∗(η −1)
pβ,kk 0 .. .
0 .. . 0
0
0
···
0
,
where ηβ = knk − nβ . For β = 1, 2, . . . , k, also define the knk × kknk residual error matrices ¯ β , 0n , . . . , 0 n ] δWβ = [δ W 1 k and ¯ 0n , . . . , 0n ], δR = [δ R, 1 k where
¯ δ Wβ =
(knk−1)
δwβ
···
.. .
(0) δwβ
.. (0)
δwβ
(n0 )
δwβ .. .
.
¯= , δR
δr(knk−2)
δr (n0 −1) .. . (0) δr , 0 .. . 0
···
.. . .. δr(0) 0
.
. .. ···
and 0nβ is a knk × knk − nβ matrix of zeroes. Also, let θ0,0 · · · θ0,k .. , θ = ... . θk,0 · · · θk,k
where each θα,β is an (knk − nα ) × (knk − nβ ) matrix given by (2knk−nβ ) (knk+1) (θII )α,β ··· (θII )α,β .. .. θα,β = . . (knk+nα −nβ +1) (nα +2) (θII )α,β · · · (θII )α,β 16
i h (i−j) denote an order knk, with θII (z) the error appearing in (18). Finally, let a0 lower triangular, matrix as in §4. The main result of this section is Theorem 7 below which gives the inverse of M ∗n in terms of the NPHS S(z) and the NSPS S ∗ (z) of types n for A(z). Theorem 7. In terms of the normalized NPHS S(z) and the normalized NSPS S ∗ (x) of types n for A(z), the inverse of M∗n satisfies k o h h n i−1 i−1 X (i−j) (i−j) (0) −1 t ∗ t ∗ = Q a a (37) M∗−1 (a ) I + θ V Q∗β , V + kknk IV β n 0 0 0 β=1
where ∗ θIV
(38)
=
k i−1 h i−1 h X (i−j) (i−j) Q∗β δWβt a0 V∗ − θ − δRt a0 β=1
Proof. Let
¯ Q=
p(knk−2)
···
.. .
p(n0 −1) .. .
(knk−2)
q1
.. .
p(0) ..
p(0) 0
.
.
···
..
(0)
q1 0
.
.
(knk−2)
qk
.. . (0)
..
0
(n1 −1)
q1
q1
···
.. .
...
.. . (0)
..
0
..
(nk −1)
qk
(0)
.
qk 0
.
qk
0
..
Then, the order condition (5) for an NPHS implies that h i ¯ t · a(i−j) − δRt . (39) M∗n · Qt = Q 0
To see this, note the (i, j)th component, 1 ≤ i ≤ knk − n0 , 1 ≤ j ≤ knk, of (39) is the coefficient of z knk−i−j in a0 (z) p(z) +
k X
aα (z) qα (z) = z knk−1 r(z) + δr(z).
α=1
The remaining components of (39) are obvious identities. Similarly, for 1 ≤ β ≤ k, let (knk−1) (knk−1) (n ) (n ) · · · v1,β1 uβ · · · uβ 0 v1,β .. .. . . .. .. ¯ (0) (0) Vβ = . uβ . v1,β . . . . . .. .. (0)
uβ
(0)
v1,β
17
(knk−1)
vk,β
···
.. . .. (0)
vk,β
.
(n ) vk,βk .. . (0) . vk,β
.
Then, the coefficient of z knk−i−j+1 , 1 ≤ i ≤ knk − n0 , 1 ≤ j ≤ knk, in the order condition (5) for an NPHS, namely, a0 (z) uβ (z) +
k X
aα (z) vα,β (z) = z knk+1 wβ (z) + δwβ (z),
α=1
gives the (i, j)th component of i h (i−j) − δW t . M∗n · Vβt = V¯βt · a0
(40)
The remaining components of (40) are easy to verify. Next, observe that Theorem 1 and Corollary 2 imply that ¯t · V ∗ + Q
(41)
k X
(0) V¯βt · Q∗β = (a0 )−1 Ikknk + θ.
β=1
Combining (39), (40) and (41), we obtain the result (37). Corollary 8 below drops the requirement in Theorem 7 that S(z) and S ∗ (z) be normalized. In particular, the results of the corollary apply when S(z) and S ∗ (z) are scaled. Corollary 8. In terms of the NPHS S(z) (unnormalized) of type n for A(z) and the NSPS S ∗ (z) (unnormalized) of type n for A∗ (z), the inverse of M∗n is given by o n (0) ∗ (42) (a0 )−1 Ikknk + θ¨IV M∗−1 n =
(γ0 γ0∗ )−1 Qt
k i−1 h i−1 h X (i−j) (i−j) ∗ Q∗β , (γβ γβ∗ )−1 Vβt a0 a0 V + β=1
where ∗ θ¨IV
=
k h i−1 h i−1 X (i−j) (i−j) V∗ − θ¨ − (γ0 γ0∗ )−1 δRt a0 (γβ γβ∗ )−1 δWβt a0 Q∗β β=1
and θ¨0,0 θ¨ = ... θ¨k,0
with θ¨α,β
(knk+1)
(θ¨II )α,β .. = . (nα +2) ¨ (θII )α,β
··· ···
··· ···
θ¨0,k .. . ¨ θk,k (2knk−nβ ) (θ¨II )α,β .. . (knk+n −n +1) α β ¨ (θII )α,β
Proof. The normalized NPHS is obtained from an unnormalized one by multiplying it on the right by the diagonal matrix diag[γ0−1 , . . . , γk−1 ]. Similarly, the 18
normalized NSPS is obtained from an unnormalized one by multiplying it on the left by the diagonal matrix diag[γ0∗−1 , . . . , γk∗−1 ]. The result now follows directly from (37). Corollary 9. If the conditions of Corollary 6 are satisfied, then3 kM∗−1 n k∞
(43)
≤
(0)
knk+1 2κ · |a0 | · ka−1 )k. 0 (z) (mod z
Proof. ¿From (36), ¨ ∞ kθk
≤
(k + 1)kθ¨II (z)k
≤
© ª (0) knk+1 2κ2 (k + 1)2 · |a0 | · ka−1 )k · (k + 1)kδT t (z)k + kδT ∗ (z)k . 0 (z) (mod z
Thus, ∗ k∞ kθ¨IV
=
k i−1 i−1 h h X (i−j) (i−j) V∗ − Q∗β k∞ kθ¨ − (γ0 γ0∗ )−1 δRt a0 (γβ γβ∗ )−1 δWβt a0 β=1
≤
¨ ∞+ kθk +
k X
(γ0 γ0∗ )−1 (k
knk+1 + 1) · kδT (z)k · ka−1 )k · kS ∗ (z)k 0 (z) (mod z t
knk+1 (γβ γβ∗ )−1 (k + 1) · kδT t (z)k · ka−1 )k · kS ∗ (z)k 0 (z) (mod z
β=1 knk+1 κ(k + 1)2 · ka−1 )k 0 (z) (mod z n £ ¤o (0) · kδT t (z)k + 2κ|a0 | · (k + 1)kδT t (z)k + kδT ∗ (z)k h i2 (0) knk+1 ≤ 4|a0 | · (κ + 1)(k + 2)(ka−1 (z) (mod z )k + 1) 0 £ ¤ t ∗ · (k + 2)kδT (z)k + kδT (z)k .
≤
Therefore, using the assumption (34),
n o−1 (0) (0) ∗ k (a0 )−1 Ikknk + θ¨IV k∞ ≤ 2|a0 |
and so kM∗−1 n k∞
≤
h n i−1 o−1 (0) (i−j) ∗ k (a0 )−1 Ikknk + θ¨IV k∞ · k(γ0 γ0∗ )−1 Qt a0 V∗ +
k X
β=1
≤
i−1 h (i−j) (γβ γβ∗ )−1 Vβt a0 Q∗β k∞
(0)
2κ|a0 | · ka0−1 (z) (mod z knk+1 )k.
6. Conclusions. In this paper we have presented new formulae for the “near” inverses of striped and mosaic Sylvester matrices. The formulae are given in terms of numerical Pad´e-Hermite and simultaneous Pad´e systems. They are important for 3
The ∞-norm, rather than the 1-norm, is used here because it is more suitable for purposes in
[6]. 19
numerical computation since they incorporate errors caused by floating-point arithmetic. In particular, the formulae can be used to determine good estimates for the condition numbers of these matrices. Our primary motivation for obtaining these formulae is the numerically stable computation of Pad´e-Hermite and simultaneous Pad´e approximants, the subject of the companion paper [6]. As such we have restricted our attention to a striped and a specific mosaic Sylvester matrices. We conjecture that a similar approach can also be used for determining near inverse formulae of other structured matrices, for example, of mosaic Hankel, Toeplitz or Sylvester matrices [13, 16]. Some preliminary work on this topic has already been done in [8]. Together with the results of [6], we believe that the formulae given in this paper can be used to stably invert striped and mosaic Sylvester matrices and to stably solve systems of linear equations with these as coefficient matrices. This matter requires formal verification, such as that reported in [11] for the case k = 1 and a0 (z) = 1. Acknowledgement. We are very greatful to a referee who contributed much in terms of the correctness of results and the clarity of presentation. REFERENCES [1] M. V. Barel and A. Bultheel, The computation of non-perfect Pad´ e-Hermite approximants, Numerical Algorithms, 1 (1991), pp. 285–304. [2] B. Beckermann, Zur Interpolation mit polynomialen Linearkombinationen beliebiger Funktionen, PhD thesis, Institut f¨ ur Angewandte Mathematik, Universit¨ at Hannover, 1988. [3] B. Beckermann, A reliable method for computing M-Pad´ e approximants on arbitrary staircases, Journal of Computational and Applied Mathematics, 40 (1992), pp. 19–42. [4] B. Beckermann and G. Labahn, A uniform approach for the fast computation of matrix-type Pad´ e approximants, SIAM Journal on Matrix Analysis and Applications, (1994), pp. 804– 823. [5] S. Cabay, A. Jones, and G. Labahn, A stable algorithm for multi-dimensional Pad´ e systems and the inversion of generalized Sylvester matrices, Tech. Report TR 94-07, Dept. Comp. Sci., Univ. Alberta, 1994. , Computation of numerical Pad´ e-Hermite and simultaneous Pad´ e systems II: A weakly[6] stable algorithm, (in preparation). [7] S. Cabay and G. Labahn, A superfast algorithm for multi-dimensional Pad´ e systems, Numerical Algorithms, 2 (1992), pp. 201–224. [8] , Fast,stable inversion of mosaic Hankel matrices, Systems and Networks: Mathematical Theory and Applications, Proceedings of MTNS 93, (1994), pp. 625–630. [9] S. Cabay, G. Labahn, and B. Beckermann, On the theory and computation of non-perfect Pad´ e-Hermite approximants, Journal of Computational and Applied Mathematics, 39 (1992), pp. 295–313. [10] S. Cabay and R. Meleshko, A weakly stable algorithm for Pad´ e approximants and inversion of Hankel matrices, SIAM Journal on Matrix Analysis and Applications, 14 (1993), pp. 735– 765. [11] M. H. Gutknecht and M. Hochbruck, The stability of inversion formulas for Toeplitz matrices, Tech. Report IPS 93–13, IPS-Z¨ urich, 1993. [12] G. Heinig and K. Rost, Algebraic Methods for Toeplitz-like Matrices and Operators, Birkhauser Verlag, Basel, 1984. [13] G. Heinig and A. Tewodros, On the inverses of Hankel and Toeplitz mosaic matrices, Seminar Analysis Operator Equation and Numerical Analysis 1987/1988, (1988), pp. 53–65. [14] T. Jones, The numerical computation of Pad´ e-Hermite systems, master’s thesis, Dept. Comp. Sci., Univ. Alberta, 1992. [15] G. Labahn, Inversion components for block Hankel-like matrices, Linear Algebra and Its Applications, 177 (1992), pp. 7–48. [16] G. Labahn, B. Beckermann, and S. Cabay, Inversion of mosaic Hankel matrices via matrix polynomial systems, to appear in Linear Algebra and its Applications, (1995). [17] G. Labahn, D. K. Choi, and S. Cabay, The inverses of block Hankel and block Toeplitz matrices, SIAM J. Comput., 19 (1990), pp. 98–123. 20
[18] G. Labahn and T. Shalom, Inversion of Toeplitz structured matrices using only standard equations, Linear Algebra and Its Applications, (1994), pp. 49–70. [19] L. Lerer and M. Tismenetsky, Generalized Bezoutians and the inversion problem for block matrices, Integral Equations and Operator Theory, 9 (1986), pp. 790–819. [20] , Toeplitz classification of matrices and inversion formulas, ii. block-Toeplitz and perturbed block-Toeplitz matrices, Tech. Report 88.197, IBM-Israel Scientific Center, Haifa, 1986. [21] K. Mahler, Perfect systems, Compositio Math., 19 (1968), pp. 95–166. [22] G. W. Stewart, Introduction to Matrix Computations, Academic Press, 1973.
21
1. Introduction. Let n = [n0 , . . . , nk ] be a vector of integers with nβ ≥ 0, 0 ≤ β ≤ k. A striped Sylvester matrix of order knk is given by (0) (0) a0 ak .. .. . . . . , (0) (0) . . Mn = (1) . a0 ··· . ak .. .. . . (knk−1) (knk−n0 ) (knk−1) (knk−nk ) a0 · · · a0 ak · · · ak (`)
where knk = n0 + . . . + nk and where the coefficients aβ
∈ R, the field of real
(0) a0
numbers. Assume that 6= 0. In this paper, we present a formula for the inverse of Mn expressed in terms of Pad´e-Hermite and simultaneous Pad´e systems. Pad´e-Hermite and simultaneous Pad´e systems [7, 9] are approximants for the vector At (z) = [a0 (z), . . . , ak (z)] of power series associated with the coefficients of Mn , namely, knk−1
aβ (z) =
X
(`)
aβ z ` , with 0 ≤ β ≤ k.
`=0
They provide necessary and sufficient conditions for Mn to be nonsingular, and generalize the notions of Pad´e-Hermite and simultaneous Pad´e approximants. Pad´eHermite and simultaneous Pad´e approximants were introduced by Hermite in the last ∗ Department of Computing Science, University of Alberta, Edmonton, Alberta, Canada, T6G 2H1. The research of this author as partially supported by Natural Sciences and Engineering Research Council of Canada grant A8035. † Bell Northern Research, P.O. Box 3511, Station C, Ottawa, Ontario, Canada, K1Y 4H7 ‡ Department of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, N2L3G1. The research of this author was partially supported by Natural Sciences and Engineering Research Council of Canada grant FS1525C.
1
century and has been widely studied by several authors (for a bibliography, see, for example [1, 2, 3, 4, 15]). The inverse formula given here is intended to be applied in a numerical setting; it accommodates errors that may have been made in the computation of Pad´e-Hermite and simultaneous Pad´e systems. That is, the formula gives the “near” inverse for Mn since it expresses the inverse in terms of Pad´e-Hermite and simultaneous Pad´e systems which are computed using floating-point arithmetic. There are other closed formulae (c.f. [12, 16, 18, 19, 20]) for M−1 n . The formula given here is different in that it expresses the inverse directly in terms of numerical Pad´e-Hermite and simultaneous Pad´e systems. The near inverse formula depends on the computation of both Pad´e systems. It is possible to determine a simultaneous Pad´e system from its “dual” Pad´e-Hermite system via the adjoint operation [6, 15]. In a numerical setting, however, this does not provide enough control over the resulting floating-point errors [14]. Instead, simultaneous Pad´e systems can be computed independently. Whereas a Pad´e-Hermite system can be obtained by solving a system of linear equations with M n as the coefficient matrix, a simultaneous Pad´e system can be similarly and independently obtained with a coefficient matrix that is a specialized mosaic Sylvester matrix. This specialized mosaic Sylvester matrix of order kknk is defined to be ∗ ∗ S0,1 · · · S0,k .. , M∗n = ... (2) . ∗ ∗ Sk,1 · · · Sk,k ∗ where Sα,β are matrices of size (knk − nα ) × knk, with (0) (knk−1) aβ ··· aβ .. ∗ .. S0,β = − . . (0) (n0 ) aβ ··· aβ
∗ Sβ,β =
(0)
a0
(knk−1)
··· ..
a0
.. .
. (0)
a0
···
(nβ )
a0
,
,
∗ = 0. The matrix M∗n is closely related to for 1 ≤ β ≤ k, and with the remaining Sα,β Mn . Indeed, in the special case when k = 1 the matrix Mn and the transpose of M∗n coincide, up to a sign and a permutation of the rows and columns. Similarly, when k ≥ 1 and a0 (z) = 1, the matrix Mn and the transpose of M∗n are both obtained by a suitable block extension of the same matrix. In this paper, we also provide a “near” inverse formula for the matrix M∗n , again in terms of numerical Pad´e-Hermite and simultaneous Pad´e systems. The inverse formulae provide “good” estimates for the condition numbers of M n and M∗n and these enable the stable numerical computation of Pad´e-Hermite and simultaneous Pad´e systems described in the companion paper [6]. In return, the accurate computation of Pad´e-Hermite and simultaneous Pad´e systems by the algorithm in [6] enables the effective inversion of generalized Sylvester matrices by the formulae given in this paper, as well as the solution of linear systems of equations with these as the coefficient matrices.
2
This paper is organized as follows. Preliminary definitions and basic facts about Pad´e-Hermite and simultaneous Pad´e systems are given in the next two sections. §3 gives a near commutativity relationship between these two systems in floating-point arithmetic, while §4 and §5 give the approximate inversion formulae for striped and mosaic Sylvester matrices. The final section gives a summary and some conclusions. We conclude this section by defining some norms which are used in §4 and §5. Let a(z) =
∞ X
a(`) z ` ∈ R [[z]] ,
`=0
where R[[z]] is the domain of power series with coefficients from R, and define the bounded power series ¯ ) ( ∞ ¯ X ¯ (`) B |a | < ∞ . R [[z]] = a(z) ¯ a(z) ∈ R [[z]] , ¯ `=0
A norm for a(z) ∈ RB [[z]] is
ka(z)k =
∞ X
|a(`) |.
`=0
RB [[z]] includes the domain of polynomials R[z]. So, for s(z) =
∂ X
s(`) z ` ∈ R [z] ,
`=0
we use the norm ks(z)k =
∂ X
| s(`) | .
`=0
For vectors and matrices over RB [[z]], we use the 1-norm unless otherwise specified. So, for example, the norm for At (z) is kAt (z)k = max {kaβ (z)k} 0≤β≤k
and the norm for S(z) ∈ R(k+1)×(k+1) [z] is kS(z)k = max
0≤β≤k
(
k X
α=0
)
kSα,β (z)k .
It is easy to verify that various compatibility conditions are satisfied. For example, kAt (z) · S(z)k ≤ kAt (z)k · kS(z)k and ka(z) · b(z)k ≤ ka(z)k · kb(z)k, 3
where b(z) is also a bounded power series. In addition, for S ∗ (z) ∈ R(k+1)×(k+1) [z] and A∗ (z) ∈ RB (k+1)×k [[z]], kS ∗ (z) · A∗ (z)k ≤ kS ∗ (z)k · kA∗ (z)k, kS(z) · S ∗ (z)k ≤ kS(z)k · kS ∗ (z)k. In the subsequent development, we also make use of the inequality ka(z) (mod z knk+1 )k ≤ ka(z)k, where a(z) (mod z knk+1 ) =
knk X
a(`) z ` +
`=0
∞ X
0 · z ` ∈ RB [[z]]
`=knk+1
2. Multi-dimensional Pad´ e Systems. In this section, we introduce the notion of Pad´e-Hermite and simultaneous Pad´e systems. Let n = [n0 , . . . , nk ] and define knk = n0 + · · · + nk . Also let At (z) = [a0 (z), . . . , ak (z)] , where aβ (z) =
∞ X
(`)
aβ z ` ,
β = 0, . . . , k,
`=0
(0)
(`)
with aβ ∈ R, the field of real numbers. Assume that a0 6= 0, which means that t knk+1 )k = 1, a−1 0 (z) exists. Assume also that A (z) is scaled so that kaβ (z) (mod z 0 ≤ β ≤ k. The (k + 1) × (k + 1) matrix of polynomials zp(z) u1 (z) · · · uk (z) · ¸ zq1 (z) v1,1 (z) · · · v1,k (z) zp(z) U t (z) S(z) = (3) = .. .. .. zQ(z) V (z) . . . zqk (z) vk,1 (z) · · · vk,k (z) is a numerical Pad´e-Hermite system (NPHS) [9] of type n for A(z) if the following conditions are satisfied. I. (Degree conditions): For 1 ≤ α, β ≤ k, (4)
p(z) =
nX 0 −1
p(`) z ` ,
uβ (z) =
`=0
qα (z) =
nX α −1
n0 X `=0
qα(`) z ` ,
vα,β (z) =
`=0
nα X `=0
II. (Order condition): (5)
(`)
uβ z ` ,
At (z)S(z) = z knk T t (z) + δT t (z), 4
(`)
vα,β z ` .
where T t (z) = [r(z), zW t (z)] with W t (z) = [w1 (z), . . . , wk (z)] is the residual, and where δT t (z) = [z δr(z), δW t (z)] is the residual error, with δW t (z) = [δw1 (z), . . . , δwk (z)] and with knk−2
δr(z) =
X
δr
(`)
`
z , δwβ (z) =
`=0
knk X
(`)
δwβ z ` .
`=0
If δT t (z) = 0, then S(z) is an exact, rather than a numerical, Pad´e-Hermite system. III. (Nonsingularity condition): The constant term of V (z) is a diagonal matrix, (6)
V (0) = diag [γ1 , . . . , γk ] ,
and (7)
(0)
γ ≡ (a0 )−1
k Y
γα 6= 0,
α=0
where γ0 = r(0). Remark 1: The nonsingularity condition III is equivalent to the condition that r(0) 6= 0 and that V (0) be a nonsingular diagonal matrix. Remark 2: The NPHS is said to be normalized [9] if the nonsingularity condition III is replaced by r(0) = 1 and V (0) = Ik . This can be achieved by multiplying S(z) on the right by Γ−1 , where (8)
Γ = diag [γ0 , . . . , γk ] .
The NPHS is said to be scaled [14] if each column of S(z) has norm equal to 1 and if, in addition, γβ > 0, 0 ≤ β ≤ k. Here, also, scaling a NPHS is accomplished by multiplying it on the right by an appropriate diagonal matrix. Remark 3: The nonsingularity condition III, namely γ 6= 0, refers to the nonsingularity of the associated striped Sylvester matrix Mn defined in (1); in [9] it is shown that an exact NPHS exists iff Mn is nonsingular. Remark 4: In [5, 6, 9], rather than S(z), the Pad´e-Hermite system is defined to be S(z) · diag[z, 1, . . . , 1]. The notation used here is consistent with that of [16], and simplifies the presentation of some of the results. A Pad´e-Hermite system gives an approximation to a vector of formal power series using matrix multiplication on the right. A simultaneous Pad´e system is a similar approximation using matrix multiplication on the left and with degree constraints that can be thought of as being “dual” to the degree constraints of a Pad´e-Hermite system.
5
Let1
−a1 (z) a0 (z) A∗ (z) = 0
(9)
−ak (z) 0 a0 (z)
··· ..
.
be a (k + 1) × k matrix of power series. The (k + 1) × (k + 1) matrix of polynomials ∗ ··· u∗k (z) u∗1 (z) v (z) · ∗ ¸ zq ∗ (z) zp∗ (z) · · · zp∗ (z) 1,1 1,k v (z) U ∗t (z) 1 S ∗ (z) = (10) = . . .. ∗ ∗ . . zQ (z) zP (z) . . . zqk∗ (z) zp∗k,1 (z) · · · zp∗k,k (z)
is a numerical simultaneous Pad´e system (NSPS) [7, 9] of type n for A∗ (z) if the following conditions are satisfied. I. (Degree conditions): For 1 ≤ α, β ≤ k, knk−nβ
knk−n0
(11)
v ∗ (z) =
v ∗ (`) z ` ,
X
uβ z ` ,
`=0
`=0
knk−nβ −1
knk−n0 −1
qα∗ (z) =
∗(`)
X
u∗β (z) =
X
qα∗(`) z ` ,
X
p∗α,β (z) =
∗(`)
pα,β z ` .
`=0
`=0
II. (Order condition): S ∗ (z)A∗ (z) = z knk T ∗ (z) + δT ∗ (z),
(12)
where T ∗t (z) = [z W ∗ (z)|R∗t (z)] with R∗ (z) a k × k is the residual, and where δT ∗t (z) = [δW ∗ (z)|z δR∗t (z)] is the residual error, with δR∗ (z) a k × k matrix and δwβ∗ (z) =
knk X
knk−2 ∗(`)
δwβ
∗ z ` , δrα,β (z) =
X
∗(`)
δrα,β z ` .
`=0
`=0
If δT ∗ (z) = 0, then S ∗ (z) is an exact NSPS. III. (Nonsingularity condition): The constant term of R ∗ (z) is a diagonal matrix R∗ (0) = diag [γ1∗ , . . . , γk∗ ] ,
(13) and
(0)
γ ∗ ≡ (a0 )−1
(14)
k Y
γα∗ 6= 0,
α=0 1
More generally,
A∗ (z) =
a∗0,1 (z) a∗1,1 (z) . . . ∗ ak,1 (z)
··· ···
···
a∗0,k (z) a∗1,k (z) . . . ∗ ak,k (z)
h =
B ∗t (z) C ∗ (z)
i
with C ∗ (0) nonsingular, but the restriction to (9), which algebraically is made without loss of generality, permits us to establish in §3 a duality relationship between Pad´ e-Hermite systems and simultaneous Pad´ e systems. 6
where γ0∗ = v ∗ (0). Remark 5: The NSPS is said to be normalized [7] if the nonsingularity condition III is replaced by v ∗ (0) = 1 and R∗ (0) = Ik . This can be achieved by multiplying S ∗ (z) on the left by Γ∗−1 , where Γ∗ = diag [γ0∗ , . . . , γk∗ ] .
(15)
The NSPS is said to be scaled when each row of S ∗ (z) has norm equal to 1 and if, in addition, γα∗ > 0, 0 ≤ α ≤ k. Here, also, scaling a NSPS is accomplished by multiplying it on the left by an appropriate diagonal matrix. Remark 6: The nonsingularity condition III, namely γ ∗ 6= 0, refers to the nonsingularity of the associated mosaic Sylvester matrix M∗n defined in (2); in [7] it is shown that an exact NSPS exists iff M∗n is nonsingular. Remark 7: In [5, 6, 9], rather than S ∗ (z), the simultaneous Pad´e system is defined to be diag[1, z, . . . , z] · S ∗ (z). This is for notational convenience only. 3. Duality. Theorem 1 below gives a relationship between Pad´e-Hermite and simultaneous Pad´e systems which is crucial to the results of the subsequent sections. It generalizes earlier results of Mahler and their extensions to block matrices [10, 15, 17, 21]. Theorem 1. If S(z) is a NPHS of type n for A(z) and S ∗ (z) is a NSPS of type n for A∗ (z), then (0)
S ∗ (z) · S(z) = z knk (a0 )−1 Γ∗ Γ + θI (z),
(16) where θI (z)
=
a−1 0 (z)
½·
v ∗ (z) zQ∗ (z)
¸
∗
t
δT (z) + δT (z)
£
zQ(z)
V (z)
¤
¾
(mod z knk+1 ).
Proof. The theorem (in the case that δT (z) = 0 and δT ∗ (z) = 0) follows from [15]. The arguments used in the following proof, however, are considerably simpler. Let B t (z) = [a1 (z), · · · , ak (z)]. Then, using (5) and (12), (17)
a0 (z) S ∗ (z) · S(z) ½· ∗ ¸ · ∗t ¸ ¾ £ ¤ ¤ v (z) U (z) £ t zp(z) U (z) zQ(z) V (z) = a0 (z) + zQ∗ (z) zP ∗ (z) · ∗ ¸ ¤ £ ¤ª © £ v (z) = a0 (z) zp(z) U t (z) + B t (z) zQ(z) V (z) ∗ zQ (z) ½ · ∗t ¸ · ∗ ¸ ¾ £ ¤ U (z) v (z) t zQ(z) V (z) + a0 (z) − B (z) ∗ ∗ zP (z) zQ (z) 7
= =
·
v ∗ (z) zQ∗ (z) ½·
¸
£ ¤ At (z)S(z) + S ∗ (z)A∗ (z) zQ(z) V (z) ¸ ¸ · £ ¤ v ∗ (z) zW ∗t (z) £ t r(z) W (z) zQ(z) z knk + zQ∗ (z) R∗ (z) ¸ · ∗ £ ¤ v (z) δT t (z) + δT ∗ (z) zQ(z) V (z) . + ∗ zQ (z)
V (z)
¤
¾
But, from (4) and (11), the degrees of S ∗ (z)S(z) are bounded component-wise by knk. It then follows from (17) that ¸ · ∗ 0 v (0)r(0) ∗ knk (0) −1 + θI (z) S (z)S(z) = z (a0 ) 0 R∗ (0)V (0) =
(0)
z knk (a0 )−1 Γ∗ Γ + θI (z),
which is (16). Corollary 2. If S(z) is a normalized NPHS of type n for A(z) and S ∗ (z) is a normalized NSPS of type n for A∗ (z), then for sufficiently small2 δT (z) and δT ∗ (z) (0)
S(z) · S ∗ (z) = z knk (a0 )−1 Ik+1 + θII (z),
(18) where
(0)
θII (z) = S(z) θI (z) [z knk (a0 )−1 Ik+1 + θI (z)]−1 S ∗ (z). Proof. If θI (z) is so small, that is, if δT (z) and δT ∗ (z) are so small, that (0) z knk (a0 )−1 Ik+1 + θI (z) is nonsingular, then from (16) (0)
S ∗−1 (z) = S(z) · [z knk (a0 )−1 Ik+1 + θI (z)]−1 . So, (0)
S(z)S ∗ (z) − z knk (a0 )−1 Ik+1 n o (0) = S(z) − z knk (a0 )−1 S ∗−1 (z) S ∗ (z) n o (0) (0) = S(z) Ik+1 − z knk (a0 )−1 [z knk (a0 )−1 Ik+1 + θI (z)]−1 S ∗ (z) =
(0)
S(z) θI (z) [z knk (a0 )−1 Ik+1 + θI (z)]−1 S ∗ (z).
Corollary 3. The residuals T (z) for a normalized NPHS of type n for A(z) and T ∗ (z) for a normalized NSPS of type n for A∗ (z) satisfy (0)
t T t (z) S ∗ (z) = (a0 )−1 At (z) + θIII (z),
(19) where
2
© ª t θIII (z) = At (z)θII (z) − δT t (z)S ∗ (z)) /z knk .
It is adequate, for example, that condition (34) of Corollary 6 be satisfied. 8
Proof. From (5) and (18), it follows that n o z knk T t (z) + δT t (z) S ∗ (z) = At (z) S(z) S ∗ (z) n o (0) = At (z) z knk (a0 )−1 Ik+1 + θII (z)
and so (19) is true.
4. The Inverse of a Striped Sylvester Matrix. In this section, a formula is given for the inverse of Mn expressed in terms of both S(z) and S ∗ (z). This enables estimating the condition number of Mn without explicitly computing M−1 n . Associated with the NPHS S(z), define the order knk matrices (0) (n −1) (0) (n −1) q1 · · · q1 1 qk · · · qk k p(0) · · · p(n0 −1) .. .. .. . . . .. .. .. 0 0 . . 0 . (n1 −1) (nk −1) .. .. p(n0 −1) . . . . . q · · · q 1 k P= . . . . . . 0 . 0 . 0 . .. .. .. . . . 0
0 ··· 0 (20) and, for β = 1, 2, . . . , k, (1) (n ) u · · · uβ 0 β. . . .. 0 . (n0 ) . . . u (21) Uβ = β .. 0 . .. . 0 ··· 0
(1)
v1,β .. . (n )
v1,β1
···
0
··· . .. . ..
v1,β1
0 P∞
(n )
···
(1)
vk,β .. .
0 ··· .. .
0 .. .
0
0
(n )
vk,βk
···
··· . .. . ..
0 .. .
0 .. . 0
(n )
vk,βk
···
0
0
.
(`) ` Also, for any power£ series a(z) ¤ = `=0 a z , and any integer function f (i, j), (f (i,j)) i, j = 1, 2, . . . , let a ] denote a matrix of order knk whose element in position (i, j) is a(f (i,j)) . The main result of this section is Theorem 4 below which gives the inverse of M n in terms of the NPHS S(z) and the NSPS S ∗ (x) of types n for A(z).
Theorem 4. In terms of the normalized NPHS S(z) and the normalized NSPS S ∗ (x) of types n for A(z), the inverse of Mn satisfies k h o nh i i X h i (0) (i−j) ∗(knk−i−j) P t v ∗(knk−i−j+1) + a0 + θIV = a0 (22) M−1 Uβt qβ , n β=1
where
θIV
=
(0) a0
( h
(i−j)
(θIII )0
i
−
k h i ih X (i−j) (knk+i−j) (θII )α,0 aα
α=0
9
h
+ δr(i+j−2)
ih
v ∗(knk−i−j+1)
i
k h ih i X (i+j−1) ∗(knk−i−j) . δwβ + qβ β=1
Proof. The coefficient of z i+j−2 , for i, j = 1, 2, . . . , knk, in the first component of (5), namely, a0 (z) p(z) +
k X
aα (z) qα (z) = z knk−1 r(z) + δr(z),
α=1
is n0 X
(i+j−`−2) (`)
a0
p
+
k nX α −1 X
α=1 `=0
`=0
This is the (i, j)th component of £ (−knk+i+j−1) ¤ £ (i+j−2) ¤ (23)
(i+j−`−2) (`) qα = r(−knk+i+j−1) + δr (i+j−2) . aα
+ δr
r
h
=
(knk+i−j)
a0
k X £
+
i£
p(−knk+i+j−2)
a(knk+i−j) α
α=1
¤£
¤
qα(−knk+i+j−2) + Mn P t .
¤
Similarly, the coefficient of z i+j−1 , for i, j = 1, 2, . . . , knk, in the (β + 1)st component, β = 1, . . . , k, of (5), namely, a0 (z) uβ (z) +
k X
aα (z) vα,β (z) = z knk+1 wβ (z) + δwβ (z),
α=1
is n0 X
(i+j−`−1) (`) uβ
a0
+
nα k X X
(i+j−1)
+ δwβ
.
α=1 `=0
`=0
This is the (i, j)th component of i i h h (24)
(−knk+i+j−2)
(`)
(i+j−`−1) aα vα,β = wβ
(−knk+i+j−2)
wβ
(i+j−1)
+ δwβ
= +
h
(knk+i−j)
a0
k X £
ih
a(knk+i−j) α
α=1
(−knk+i+j−1)
uβ
¤h
i
(−knk+i+j−1)
vα,β
i
+ Mn Uβt .
Next, the coefficient of z i−j−1 for i, j = 1, . . . , knk, in the first row and first column of (18) for a normalized NPHS and a normalized NSPS, namely, p(z) v ∗ (z) +
k X
(0)
uβ (z)qβ∗ (z) = z knk−1 (a0 )−1 + z −1 (θII )0,0 (z),
β=1
is nX 0 −1 `=0
v
∗(i−j−`−1)
p
(l)
+
n0 k X X
β=1 `=0
10
∗(i−j−`−1) (`) uβ
qβ
(i−j)
= (θII )0,0 .
This is the (i, j)th component of (25)
i ih h p(−knk+i+j−2) v ∗(knk−i−j+1)
k h X
+
ih
(−knk+i+j−1)
uβ
β=1
h
=
(i−j)
(θII )0,0
i
∗(knk−i−j)
qβ
i
.
The coefficient of z i−j−1 in the first column and the (α + 1)st row, α = 1, . . . , k, of (18), namely, qα (z) v ∗ (z) +
k X
vα,β (z)qβ∗ (z) = z −1 (θII )α,0 (z)
β=1
is nα X
v ∗(i−j−`−1) qα(l) +
nα k X X
∗(i−j−`−1) (`) vα,β
qβ
(i−j)
= (θII )α,0 .
β=1 `=0
`=0
This is the (i, j)th component of (26)
h
qα(−knk+i+j−2)
ih
v ∗(knk−i−j+1)
i
+
k h X
(−knk+i+j−1)
vα,β
β=1
=
h
(i−j)
(θII )α,0
i
ih
∗(knk−i−j)
qβ
i
.
Also, the coefficient of z i−j , for i, j = 1, . . . , knk in the first component of (19) for a normalized NPHS and NSPS, namely, ∗
r(z)v (z) + z
2
k X
(0)
wβ (z)qβ∗ (z) = (a0 )−1 a0 (z) + (θIII )0 (z).
β=1
is the (i, j)th component of i i h h (i−j) (0) (i−j) (27) = + (θIII )0 (a0 )−1 a0 +
h
r(−knk+i+j−1)
k h X
ih
v ∗(knk−i−j+1)
(−knk+i+j−2)
wβ
β=1
ih
i
∗(knk−i−j)
qβ
i
.
We are finally ready to prove the theorem. ¿From (23), (24), (25), (26) and (27),
h
i
k X
h
∗(knk−i−j)
i
P t v ∗(knk−i−j+1) + Uβt qβ β=1 i i h ih i h nh (knk+i−j) p(−knk+i+j−2) = r(−knk+i+j−1) + δr(i+j−2) − a0
Mn
−
k h X
α=1
(knk+i−j) aα
ih
qα(−knk+i+j−2) 11
io h
v ∗(knk−i−j+1)
i
+
k nh X
(−knk+i+j−2)
wβ
β=1
−
k h X
(knk+i−j) aα
α=1
=
ih
i
h i h ih i (i+j−1) (knk+i−j) (−knk+i+j−1) + δwβ − a0 uβ
(−knk+i+j−1)
vα,β
io h
∗(knk−i−j)
qβ
i
k h h ih i X ih i (−knk+i+j−2) ∗(knk−i−j) r(−knk+i+j−1) v ∗(knk−i−j+1) + wβ qβ β=1
+
k h h ih i X ih i (i+j−1) ∗(knk−i−j) δwβ δr(i+j−2) v ∗(knk−i−j+1) + qβ β=1
−
k X
α=0
=
h
(knk+i−j) aα
ih
(i−j)
(θII )α,0
h i (0) (i−j) (a0 )−1 a0 + θIV .
i
The result (22) now follows. Corollary 5 below drops the requirement in Theorem 4 that S(z) and S ∗ (z) be normalized. In particular, the corollary is valid in the case that S(z) and S ∗ (z) are scaled. Corollary 5. In terms of the (unnormalized) NPHS S(z) of type n for A(z) and the (unnormalized) NSPS S ∗ (z) of type n for A∗ (z), the inverse of Mn is given by o nh i (i−j) (28) M−1 a0 + θ¨IV n k h i X i h ∗(knk−i−j) (0) (γβ γβ∗ )−1 Uβt qβ , (γ0 γ0∗ )−1 P t v ∗(knk−i−j+1) + = a0 β=1
where
(29)
θ¨IV
=
(0) a0
+ +
(
(γ0 γ0∗ )−1 k X
β=1
(30)
t θ¨III (z)
=
(31)
θ¨II (z)
=
(32)
θ¨I (z)
=
n
k h h i X ih i (i−j) (i−j+1) (knk+i−j) ¨ (θIII )0 − aα (θ¨II )α,0 α=0
h
δr
(i+j−2)
h
ih
v ∗(knk−i−j+1)
(i+j−1)
(γβ γβ∗ )−1 δwβ
ih
i
∗(knk−i−j)
qβ
i
,
o At (z)θ¨II (z) − δT t (z)(Γ∗ Γ)−1 S ∗ (z)) /z knk ,
S(z) (Γ∗ Γ)−1 θ¨I (z) (Γ∗ Γ)−1 (0) ·[z knk (a0 )−1 Ik+1 + θ¨I (z)(Γ∗ Γ)−1 ]−1 S ∗ (z), ½· ∗ ¸ v (z) a−1 (z) δT t (z) 0 zQ∗ (z) £ ¤ª + δT ∗ (z) zQ(z) V (z) (mod z knk+1 ) 12
Proof. The normalized NPHS is obtained from an unnormalized one by multiplying it on the right by the diagonal matrix diag[γ0−1 , . . . , γk−1 ]. Similarly, the normalized NSPS is obtained from an unnormalized one by multiplying it on the left by the diagonal matrix diag[γ0∗−1 , . . . , γk∗−1 ]. The result now follows directly from (22). Let (33)
κ=
k X
(γβ γβ∗ )−1 .
β=0
In the corollary below, we give a bound for M−1 n in terms of κ. Corollary 6. If the residual errors δT t (z) and δT ∗ (z) associated with the scaled S(z) and the scaled S ∗ (z) are not too large, so that h
(34)
then
kM−1 n k1
(35)
(0)
knk+1 (κ + 1)(k + 2)|a0 |(ka−1 )k + 1) 0 (z) (mod z £ ¤ · (k + 2)kδT t (z)k + kδT ∗ (z)k ≤ 1/8,
≤
i2
(0)
knk+1 2κ · |a0 | · ka−1 )k. 0 (z) (mod z
Proof. We use Corollary 5 with S(z) and S ∗ (z) scaled. We begin by finding a bound for θ¨IV appearing in the inverse formula (28) for Mn . A bound for θ¨IV depends on bounds for θ¨I (z), θ¨II (z) and θ¨III (z). ¿From (16), © ª knk+1 kθ¨I (z)k ≤ ka−1 )k · (k + 1)kδT t (z)k + kδT ∗ (z)k , 0 (z) (mod z
since kS(z)k = 1 and kS ∗ (z)k ≤ k + 1. ¿From (16) and (32), note that θ¨I (z) is a matrix polynomial of at most degree knk and so, using (34), (0) (0) (0) ka0 z knk θ¨I (z −1 )(Γ∗ Γ)−1 k = ka0 θ¨I (z)(Γ∗ Γ)−1 k ≤ κ · |a0 | · kθ¨I (z)k ≤ 1/2, (0)
since k(Γ∗ Γ)−1 k ≤ κ. So as in Stewart [22, page 187], the inverse of (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 exists and o−1 n (0) k (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 k
(0)
≤
|a0 | (0) knk ¨ 1 − ka z θI (z −1 )(Γ∗ Γ)−1 k
≤
2|a0 |.
0
(0)
To determine a bound for θ¨II (z) in (31), let N = max0≤β≤k {nβ } and observe from 18 that θ¨II (z) is also a matrix polynomial, now of degree at most knk + N . Consequently, (36)
kθ¨II (z)k
kz knk+N θ¨II (z −1 )k n o © ª = k z N S(z −1 ) (Γ∗ Γ)−1 z knk θ¨I (z −1 ) (Γ∗ Γ)−1
=
13
≤
≤ ≤
o−1 n (0) · (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 [z knk S ∗ (z −1 )]k
κ2 kz N S(z −1 )k · kz knk θ¨I (z −1 )k · kz knk S ∗ (z −1 )k n o−1 (0) ·k (a0 )−1 Ik+1 + z knk θ¨I (z −1 )(Γ∗ Γ)−1 k (0) 2κ2 |a0 | · kS(z)k · kθ¨I (z)k · kS ∗ (z)k (0) 2κ2 (k + 1) · |a | · kθ¨I (z)k. 0
In addition, it now follows that a bound for θ¨III (z) in (30) is given by (0) t (z)k ≤ 2κ2 (k + 1) · |a0 | · kθ¨I (z)k + κ(k + 1) · kδT t (z)k. kθ¨III
We are now ready to give a bound for θ¨IV appearing in the inverse formula (28). But, first observe that h i (i−j) k (θ¨III ) k1 ≤ kθ¨t (z)k 0
III
and that k
k h X
(knk+i−j) aα
α=0
ih
(i+j)
(θ¨II )α,0
i
k1 ≤
k X
k(θ¨II )α,0 (z)k ≤ k(θ¨II )(z)k.
α=0
Thus, kθ¨IV k1
=
(0) ka0
+
(
(γ0 γ0∗ )−1
≤ ≤ ≤
≤
α=0
h
δr
(i+j−2)
ih
v ∗(knk−i−j+1)
i
h ih i (i+j−1) ∗(knk−i−j) + (γβ γβ∗ )−1 δwβ qβ k1 β=1 k X (0) t (γβ γβ∗ )−1 kδT t (z)k (z)k + kθ¨II (z)k + (γ0 γ0∗ )−1 kδT t (z)k + |a0 | · kθ¨III β=1 n o (0) t |a0 | kθ¨III (z)k + kθ¨II (z)k + κkδT t (z)k n o (0) (0) |a0 | κ(k + 1)kδT t (z)k + 4κ2 (k + 1)|a0 | · kθ¨I (z)k + κkδT t (z)k (0) © κ(k + 2) · |a0 | kδT t (z)k £ ¤o (0) knk+1 t ∗ +4κ|a0 | · ka−1 (z) (mod z )k (k + 1)kδT (z)k + kδT (z)k 0 h i2 £ ¤ (0) knk+1 4ka−1 (z) (mod z )k (κ + 1)(k + 2)|a | (k + 2)kδT t (z)k + kδT ∗ (z)k . 0 0 k X
≤
k h i ih i X h (i−j) (i−j) (knk+i−j) aα (θ¨II )α,0 − (θ¨III )0
It then follows from (34) that
h i−1 (i−j) knk+1 k a0 θ¨IV k1 ≤ ka−1 )k · kθ¨IV k1 ≤ 1/2, 0 (z) (mod z 14
i−1 h (i−j) and so Iknk + a0 θ¨IV is invertible. In addition, k
nh
(i−j)
a0
i
+ θ¨IV
o−1
k1
¾−1 h ½ i−1 i−1 h (i−j) (i−j) a0 k1 k Iknk + a0 θ¨IV
≤
h i−1 (i−j) k a0 k1 i−1 h (i−j) θ¨IV k1 1 − k a0 i−1 h (i−j) 2k a0 k1
≤
≤
knk+1 2ka−1 )k. 0 (z) (mod z
≤
Therefore, a bound for M−1 n in (28) is given kM−1 n k1
o−1
h i (γ0 γ0∗ )−1 Pnt v ∗(knk−i−j+1) k h i X ∗(knk−i−j) t qβ (γβ γβ∗ )−1 Un,β + k 1
≤
k
nh
(i−j)
a0
i
+ θ¨IV
(0)
k · ka0
n
β=1
(0)
knk+1 2κ|a0 | · ka−1 )k. 0 (z) (mod z
≤
¿From (35), it follows that a bound for the 1-norm condition number of Mn is (0)
−1 knk+1 kMn k1 · kM−1 )k, n k1 ≤ 2κ|a0 | · ka0 (z) (mod z
since it is assumed that each aβ (z) is scaled. 5. The Inverse of a Mosaic Sylvester Matrix. In this section, a formula is given for the inverse of M∗n expressed in terms of both S(z) and S ∗ (z). This enables estimating the condition number of M∗n without explicitly computing M∗−1 n . Associated with the NPHS S(z) and the NSPS S ∗ (z), for β = 1, 2, . . . , k, define the knk × kknk matrices (knk−1) (0) (knk−1) (0) v1,β · · · v1,β vk,β · · · vk,β .. .. . . , Vβ = .. .. ··· . . (0) (0) vk,β v1,β
Q=
(knk−2)
q1
.. . (0)
q1 0
··· . .. . ..
(0)
q1
(knk−2)
0
qk ···
.. . 0)
qk 0
15
··· . .. . ..
(0)
qk
0
,
v ∗(1) .. . v ∗(η0 ) V∗ = 0 .. . 0
and
∗(0)
qβ .. .
∗(η −1) q 0 ∗ β Qβ = 0 .. . 0
v ∗(η0 )
··· . .. ..
0
.
··· . .. ..
∗(η1 )
u1
.. . ···
··· . .. . ..
∗(1)
u1 .. .
∗(η0 −1)
0
···
∗(0)
pβ,1 .. .
∗(η −1)
pβ,11 .. .
···
..
∗(ηk )
∗(ηk )
uk
0
. .. .
0 .. .
···
0
··· . .. . ..
pβ,11
···
··· . ..
uk
0
∗(η −1)
···
∗(0)
pβ,k .. .
0 ··· .. .
0 .. . 0
0
0 .. .
0
qβ
∗(1)
uk .. .
.
0 .. .
0
∗(η1 )
u1
∗(η −1)
pβ,kk
··· . .. . ..
∗(η −1)
pβ,kk 0 .. .
0 .. . 0
0
0
···
0
,
where ηβ = knk − nβ . For β = 1, 2, . . . , k, also define the knk × kknk residual error matrices ¯ β , 0n , . . . , 0 n ] δWβ = [δ W 1 k and ¯ 0n , . . . , 0n ], δR = [δ R, 1 k where
¯ δ Wβ =
(knk−1)
δwβ
···
.. .
(0) δwβ
.. (0)
δwβ
(n0 )
δwβ .. .
.
¯= , δR
δr(knk−2)
δr (n0 −1) .. . (0) δr , 0 .. . 0
···
.. . .. δr(0) 0
.
. .. ···
and 0nβ is a knk × knk − nβ matrix of zeroes. Also, let θ0,0 · · · θ0,k .. , θ = ... . θk,0 · · · θk,k
where each θα,β is an (knk − nα ) × (knk − nβ ) matrix given by (2knk−nβ ) (knk+1) (θII )α,β ··· (θII )α,β .. .. θα,β = . . (knk+nα −nβ +1) (nα +2) (θII )α,β · · · (θII )α,β 16
i h (i−j) denote an order knk, with θII (z) the error appearing in (18). Finally, let a0 lower triangular, matrix as in §4. The main result of this section is Theorem 7 below which gives the inverse of M ∗n in terms of the NPHS S(z) and the NSPS S ∗ (z) of types n for A(z). Theorem 7. In terms of the normalized NPHS S(z) and the normalized NSPS S ∗ (x) of types n for A(z), the inverse of M∗n satisfies k o h h n i−1 i−1 X (i−j) (i−j) (0) −1 t ∗ t ∗ = Q a a (37) M∗−1 (a ) I + θ V Q∗β , V + kknk IV β n 0 0 0 β=1
where ∗ θIV
(38)
=
k i−1 h i−1 h X (i−j) (i−j) Q∗β δWβt a0 V∗ − θ − δRt a0 β=1
Proof. Let
¯ Q=
p(knk−2)
···
.. .
p(n0 −1) .. .
(knk−2)
q1
.. .
p(0) ..
p(0) 0
.
.
···
..
(0)
q1 0
.
.
(knk−2)
qk
.. . (0)
..
0
(n1 −1)
q1
q1
···
.. .
...
.. . (0)
..
0
..
(nk −1)
qk
(0)
.
qk 0
.
qk
0
..
Then, the order condition (5) for an NPHS implies that h i ¯ t · a(i−j) − δRt . (39) M∗n · Qt = Q 0
To see this, note the (i, j)th component, 1 ≤ i ≤ knk − n0 , 1 ≤ j ≤ knk, of (39) is the coefficient of z knk−i−j in a0 (z) p(z) +
k X
aα (z) qα (z) = z knk−1 r(z) + δr(z).
α=1
The remaining components of (39) are obvious identities. Similarly, for 1 ≤ β ≤ k, let (knk−1) (knk−1) (n ) (n ) · · · v1,β1 uβ · · · uβ 0 v1,β .. .. . . .. .. ¯ (0) (0) Vβ = . uβ . v1,β . . . . . .. .. (0)
uβ
(0)
v1,β
17
(knk−1)
vk,β
···
.. . .. (0)
vk,β
.
(n ) vk,βk .. . (0) . vk,β
.
Then, the coefficient of z knk−i−j+1 , 1 ≤ i ≤ knk − n0 , 1 ≤ j ≤ knk, in the order condition (5) for an NPHS, namely, a0 (z) uβ (z) +
k X
aα (z) vα,β (z) = z knk+1 wβ (z) + δwβ (z),
α=1
gives the (i, j)th component of i h (i−j) − δW t . M∗n · Vβt = V¯βt · a0
(40)
The remaining components of (40) are easy to verify. Next, observe that Theorem 1 and Corollary 2 imply that ¯t · V ∗ + Q
(41)
k X
(0) V¯βt · Q∗β = (a0 )−1 Ikknk + θ.
β=1
Combining (39), (40) and (41), we obtain the result (37). Corollary 8 below drops the requirement in Theorem 7 that S(z) and S ∗ (z) be normalized. In particular, the results of the corollary apply when S(z) and S ∗ (z) are scaled. Corollary 8. In terms of the NPHS S(z) (unnormalized) of type n for A(z) and the NSPS S ∗ (z) (unnormalized) of type n for A∗ (z), the inverse of M∗n is given by o n (0) ∗ (42) (a0 )−1 Ikknk + θ¨IV M∗−1 n =
(γ0 γ0∗ )−1 Qt
k i−1 h i−1 h X (i−j) (i−j) ∗ Q∗β , (γβ γβ∗ )−1 Vβt a0 a0 V + β=1
where ∗ θ¨IV
=
k h i−1 h i−1 X (i−j) (i−j) V∗ − θ¨ − (γ0 γ0∗ )−1 δRt a0 (γβ γβ∗ )−1 δWβt a0 Q∗β β=1
and θ¨0,0 θ¨ = ... θ¨k,0
with θ¨α,β
(knk+1)
(θ¨II )α,β .. = . (nα +2) ¨ (θII )α,β
··· ···
··· ···
θ¨0,k .. . ¨ θk,k (2knk−nβ ) (θ¨II )α,β .. . (knk+n −n +1) α β ¨ (θII )α,β
Proof. The normalized NPHS is obtained from an unnormalized one by multiplying it on the right by the diagonal matrix diag[γ0−1 , . . . , γk−1 ]. Similarly, the 18
normalized NSPS is obtained from an unnormalized one by multiplying it on the left by the diagonal matrix diag[γ0∗−1 , . . . , γk∗−1 ]. The result now follows directly from (37). Corollary 9. If the conditions of Corollary 6 are satisfied, then3 kM∗−1 n k∞
(43)
≤
(0)
knk+1 2κ · |a0 | · ka−1 )k. 0 (z) (mod z
Proof. ¿From (36), ¨ ∞ kθk
≤
(k + 1)kθ¨II (z)k
≤
© ª (0) knk+1 2κ2 (k + 1)2 · |a0 | · ka−1 )k · (k + 1)kδT t (z)k + kδT ∗ (z)k . 0 (z) (mod z
Thus, ∗ k∞ kθ¨IV
=
k i−1 i−1 h h X (i−j) (i−j) V∗ − Q∗β k∞ kθ¨ − (γ0 γ0∗ )−1 δRt a0 (γβ γβ∗ )−1 δWβt a0 β=1
≤
¨ ∞+ kθk +
k X
(γ0 γ0∗ )−1 (k
knk+1 + 1) · kδT (z)k · ka−1 )k · kS ∗ (z)k 0 (z) (mod z t
knk+1 (γβ γβ∗ )−1 (k + 1) · kδT t (z)k · ka−1 )k · kS ∗ (z)k 0 (z) (mod z
β=1 knk+1 κ(k + 1)2 · ka−1 )k 0 (z) (mod z n £ ¤o (0) · kδT t (z)k + 2κ|a0 | · (k + 1)kδT t (z)k + kδT ∗ (z)k h i2 (0) knk+1 ≤ 4|a0 | · (κ + 1)(k + 2)(ka−1 (z) (mod z )k + 1) 0 £ ¤ t ∗ · (k + 2)kδT (z)k + kδT (z)k .
≤
Therefore, using the assumption (34),
n o−1 (0) (0) ∗ k (a0 )−1 Ikknk + θ¨IV k∞ ≤ 2|a0 |
and so kM∗−1 n k∞
≤
h n i−1 o−1 (0) (i−j) ∗ k (a0 )−1 Ikknk + θ¨IV k∞ · k(γ0 γ0∗ )−1 Qt a0 V∗ +
k X
β=1
≤
i−1 h (i−j) (γβ γβ∗ )−1 Vβt a0 Q∗β k∞
(0)
2κ|a0 | · ka0−1 (z) (mod z knk+1 )k.
6. Conclusions. In this paper we have presented new formulae for the “near” inverses of striped and mosaic Sylvester matrices. The formulae are given in terms of numerical Pad´e-Hermite and simultaneous Pad´e systems. They are important for 3
The ∞-norm, rather than the 1-norm, is used here because it is more suitable for purposes in
[6]. 19
numerical computation since they incorporate errors caused by floating-point arithmetic. In particular, the formulae can be used to determine good estimates for the condition numbers of these matrices. Our primary motivation for obtaining these formulae is the numerically stable computation of Pad´e-Hermite and simultaneous Pad´e approximants, the subject of the companion paper [6]. As such we have restricted our attention to a striped and a specific mosaic Sylvester matrices. We conjecture that a similar approach can also be used for determining near inverse formulae of other structured matrices, for example, of mosaic Hankel, Toeplitz or Sylvester matrices [13, 16]. Some preliminary work on this topic has already been done in [8]. Together with the results of [6], we believe that the formulae given in this paper can be used to stably invert striped and mosaic Sylvester matrices and to stably solve systems of linear equations with these as coefficient matrices. This matter requires formal verification, such as that reported in [11] for the case k = 1 and a0 (z) = 1. Acknowledgement. We are very greatful to a referee who contributed much in terms of the correctness of results and the clarity of presentation. REFERENCES [1] M. V. Barel and A. Bultheel, The computation of non-perfect Pad´ e-Hermite approximants, Numerical Algorithms, 1 (1991), pp. 285–304. [2] B. Beckermann, Zur Interpolation mit polynomialen Linearkombinationen beliebiger Funktionen, PhD thesis, Institut f¨ ur Angewandte Mathematik, Universit¨ at Hannover, 1988. [3] B. Beckermann, A reliable method for computing M-Pad´ e approximants on arbitrary staircases, Journal of Computational and Applied Mathematics, 40 (1992), pp. 19–42. [4] B. Beckermann and G. Labahn, A uniform approach for the fast computation of matrix-type Pad´ e approximants, SIAM Journal on Matrix Analysis and Applications, (1994), pp. 804– 823. [5] S. Cabay, A. Jones, and G. Labahn, A stable algorithm for multi-dimensional Pad´ e systems and the inversion of generalized Sylvester matrices, Tech. Report TR 94-07, Dept. Comp. Sci., Univ. Alberta, 1994. , Computation of numerical Pad´ e-Hermite and simultaneous Pad´ e systems II: A weakly[6] stable algorithm, (in preparation). [7] S. Cabay and G. Labahn, A superfast algorithm for multi-dimensional Pad´ e systems, Numerical Algorithms, 2 (1992), pp. 201–224. [8] , Fast,stable inversion of mosaic Hankel matrices, Systems and Networks: Mathematical Theory and Applications, Proceedings of MTNS 93, (1994), pp. 625–630. [9] S. Cabay, G. Labahn, and B. Beckermann, On the theory and computation of non-perfect Pad´ e-Hermite approximants, Journal of Computational and Applied Mathematics, 39 (1992), pp. 295–313. [10] S. Cabay and R. Meleshko, A weakly stable algorithm for Pad´ e approximants and inversion of Hankel matrices, SIAM Journal on Matrix Analysis and Applications, 14 (1993), pp. 735– 765. [11] M. H. Gutknecht and M. Hochbruck, The stability of inversion formulas for Toeplitz matrices, Tech. Report IPS 93–13, IPS-Z¨ urich, 1993. [12] G. Heinig and K. Rost, Algebraic Methods for Toeplitz-like Matrices and Operators, Birkhauser Verlag, Basel, 1984. [13] G. Heinig and A. Tewodros, On the inverses of Hankel and Toeplitz mosaic matrices, Seminar Analysis Operator Equation and Numerical Analysis 1987/1988, (1988), pp. 53–65. [14] T. Jones, The numerical computation of Pad´ e-Hermite systems, master’s thesis, Dept. Comp. Sci., Univ. Alberta, 1992. [15] G. Labahn, Inversion components for block Hankel-like matrices, Linear Algebra and Its Applications, 177 (1992), pp. 7–48. [16] G. Labahn, B. Beckermann, and S. Cabay, Inversion of mosaic Hankel matrices via matrix polynomial systems, to appear in Linear Algebra and its Applications, (1995). [17] G. Labahn, D. K. Choi, and S. Cabay, The inverses of block Hankel and block Toeplitz matrices, SIAM J. Comput., 19 (1990), pp. 98–123. 20
[18] G. Labahn and T. Shalom, Inversion of Toeplitz structured matrices using only standard equations, Linear Algebra and Its Applications, (1994), pp. 49–70. [19] L. Lerer and M. Tismenetsky, Generalized Bezoutians and the inversion problem for block matrices, Integral Equations and Operator Theory, 9 (1986), pp. 790–819. [20] , Toeplitz classification of matrices and inversion formulas, ii. block-Toeplitz and perturbed block-Toeplitz matrices, Tech. Report 88.197, IBM-Israel Scientific Center, Haifa, 1986. [21] K. Mahler, Perfect systems, Compositio Math., 19 (1968), pp. 95–166. [22] G. W. Stewart, Introduction to Matrix Computations, Academic Press, 1973.
21
