Solving certain matrix equations by means of Toeplitz computations: algorithms and applications

Dario A. Bini, Luca Gemignani and Beatrice Meini, Dipartimento di Matematica Università di Pisa

 bini@dm.unipi.it, gemignan@dm.unipi.it, meini@dm.unipi.it

December 20, 2001

Abstract

We show that several classes of matrix equations can be reduced to solving a quadratic matrix equation of the kind AX^{2} + BX + C = 0 where A, B, C, X are m × m matrices or semi-infinite matrices. The problem of computing the minimal solution, if it exists, of the latter equation is reduced to computing the matrix coefficients H_{0} and H_{1} of the Laurent matrix series ∑_{i=−∞}^{+∞} z^{i}H_{i} such that H(z)(zA + B + z^{−1}C) = I. Known algorithms for this computation are revisited in terms of operations among block Toeplitz matrices and new algorithms are introduced. An application to the solution of Non-Skip-Free Markov chains is shown.

Keywords: Matrix equations, Toeplitz matrices, Markov chains, FFT, Graeffe iteration.

1 Introduction

Matrix equations are encountered in the mathematical modeling of several problems of the real world. In our analysis the main motivation comes from problems in queueing theory where the numerical solution of the Markov chains which model the problem is ultimately reduced to solving equations of the kind

∑_{i=0}^{n} X^{i}A_{i} = 0  (1)

where n is a positive integer, in certain cases n may be infinite, A_{i}, i = 0, 1, …, and X are m × m matrices. An important feature of this kind of problems is that for a suitable positive integer k, the matrix I + A_{k} as well as the remaining matrices A_{i} have nonnegative entries and I + ∑_{i} A_{i} is a column stochastic matrix, i.e., e^{T} ∑ A_{i} = e^{T}, for e^{T} = (1, 1, …). The main interest in this framework is in the minimal nonnegative solution G of (1), such that G ≥ 0 and G ≤ X for any other nonnegative solution X of (1). Here and hereafter, any inequality involving matrices means an entrywise inequality. A minimal solution has also the following property: its eigenvalues are the zeros of minimum modulus of the equation det(∑_{i=0}^{n} z^{i}A_{i}) = 0. We refer the reader to the books [16], [17], [13] and to the papers [10], [11] for more details on this subject and for conditions under which G exists.

Other applications where matrix equations play an important role are dumped vibrations, ladder networks, dynamic programming, stochastic filtering, statistics, and control theory. We refer the reader to [1], [21], [12], [9] for more information and references on this matter. For instance, in [9] the non linear equation

X + BX^{−1}A − C = 0,  (2)

is analyzed under the assumptions that B = A^{T} and both A and C are positive definite, where the main goal is to compute the extremal positive definite solutions X^{+}, X^{−} such that X − X^{−} and X^{+} − X are positive definite for any other positive definite solution X. A detailed analysis of this problem together with some effective solution algorithms can be found in [14], [21].

A generalization of (2), encountered in certain queueing problems called tree-processes [13], is the following

X + ∑_{i=1}^{k} B_{i}X^{−1}A_{i} − C = 0.  (3)

Algorithms for this class of equations are analyzed in [5].

In this paper we show that equations (2) and (1), with n ≥ 2 (included the case n = ∞), can be reduced to solving the quadratic equation

X^{2}A_{2} + X A_{1} + A_{0} = 0  (4)

where the matrices A_{2}, A_{1}, A_{0} are suitably related to the blocks A_{i} of (1) and A, B, C of (2). Then we show that the minimal solution of (4), that is the one having as eigenvalues the zeros of minimum modulus of det(A_{2}λ^{2} +A_{1}λ+A_{0}) is ultimately reduced to inverting the matrix Laurent polynomial A_{0}z^{−1}+A_{1}+A_{2}z and we introduce some algorithms for this computation. These results provide a sort of unifying framework where recent algorithms designed in the literature can be better understood [6], [7] and new ones can be easily devised. Finally we propose a new algorithm based on the evaluation/interpolation technique at the Fourier points for solving (4). The results are obtained by using the strict analogy between matrix Laurent polynomials and block Toeplitz matrices. We recall that a block n×n matrix T = (T_{i,j}) with m×m blocks T_{i,j}, i, j = 1, …, n is block Toeplitz if there are m × m matrices S_{k}, k = −n + 1, …, n − 1, such that T_{i,j} = S_{i−j}, i, j = 1, …, n.

The paper is organized as follows. In section 2 we describe how the different problems can be reduced to equation (4). In section 3 we prove suitable relations between the minimal solution G of the equation (4) and the matrix coefficients H_{i}, i ∈ ℤ of the Laurent series H(z) = ∑_{i∈ℤ} z^{i}H_{i} such that H(z)A(z) = I, under the condition that A(z) = z^{−1}A_{0} + A_{1} + zA_{2} is analytic in the annulus R = {z ∈ C : r_{1} < |z| < r_{2}} for r_{1} < 1 < r_{2}. These results are used in section 4 for devising two different algorithms for computing the minimal solution of (4). Section 5 is devoted to an application of these tools to the numerical solution of Markov chains, where the case of Non-Skip-Free problems [10] is solved in a more computational effective way. The paper ends with section 6 devoted to the conclusions.

2 Reductions

In this section we show that any solution of (1), for any value of n ≥ 2 (included n = +∞) and of (2) can be obtained from a corresponding solution of (4) for suitable values of the block coefficients. Moreover, the minimality features are kept unchanged in the reduction from (1) to (4). We say that a solution G of (1) is minimal if its m eigenvalues, counted with their multiplicities, coincide with the m zeros of smallest modulus of the polynomial det A(λ), A(λ) = ∑_{i=0}^{n} A_{i}λ^{i}. In fact, it is easy to verify that any eigenvalue of G is zero of the latter polynomial since the condition x^{T} G = λx^{T} implies that 0 = x^{T} (∑_{i=0}^{n} G^{i}A_{i}) = x^{T} (∑_{i=0}^{n} A_{i}λ^{i}), i.e., det A(λ) = 0. If n = +∞ the latter property still holds in the sense that if G is a solution of (2) then the power series det A(λ) has at least m zeros which coincide with the eigenvalues of G.

The reduction of (2) to (4) is almost immediate. Multiply (2) to the left by Y = BX^{−1} and obtain that Y^{2}A − Y C + B = 0. We refer the reader to [14] for more details concerning the minimality of the solutions.

The reduction of (1) to (4) for 2 < n < +∞ can be obtained in two different ways. A first technique can be derived by a trick introduced in [10], and relies on the following argument. If G is a solution of (1) then, by a direct inspection, we may see that

X = [ \begin{matrix}0 & … & … & 0 \\ ⋮ & … & … & ⋮ \\ 0 & … & … & 0 \\ G & G^{2} & … & G^{n−1} \end{matrix} ]  (5)

is a solution of (4) where

\begin{matrix}A_{0} = [ \begin{matrix}0 & … & … & 0 \\ ⋮ & … & … & ⋮ \\ 0 & … & … & 0 \\ A_{0} & 0 & … & 0 \end{matrix} ] , & A_{1} = [ \begin{matrix}A_{1} & A_{0} & & O \\ A_{2} & ... & ... \\ ⋮ & ... & ... & A_{0} \\ A_{n−1} & … & A_{2} & A_{1} \end{matrix} ] , \\ A_{2} = [ \begin{matrix}A_{n} & A_{n−1} & … & A_{2} \\ & A_{n} & ... & ⋮ \\ & & ... & A_{n−1} \\ O & & & A_{n} \end{matrix} ] , \end{matrix}

are matrices of size m(n−1). Moreover, observe that the nonzero eigenvalues of X are eigenvalues of G^{n−1} and any eigenvalue of G powered to n−1 is eigenvalue of X. Therefore if G is a minimal solution of (1) then X in (5) is a minimal solution of (4) and viceversa.

A different reduction can be obtained, by following [18], [8]. In fact, a direct inspection shows that if G is a solution of (2) then the matrix of size mn

X = [ \begin{matrix}G & G^{2} & … & G^{n} \\ 0 & 0 & … & 0 \\ ⋮ & ... & ... & ⋮ \\ 0 & … & … & 0 \end{matrix} ]  (6)

solves (4) where

\begin{matrix}A_{0} = [ \begin{matrix}A_{0} & 0 & … & 0 \\ 0 & 0 & ... & 0 \\ ⋮ & ... & ... & ⋮ \\ 0 & … & … & 0 \end{matrix} ] , & A_{1} = [ \begin{matrix}A_{1} & 0 & … & … & 0 \\ A_{2} & I & 0 & ... & 0 \\ A_{3} & 0 & I & ... & 0 \\ ⋮ & ... & ... & ... & 0 \\ A_{n} & 0 & … & 0 & I \end{matrix} ] , \\ A_{2} = [ \begin{matrix}0 & −I & 0 & … & 0 \\ & 0 & −I & ... & ⋮ \\ & & ... & ... & 0 \\ & & & 0 & −I \\ & ◯ & & & 0 \end{matrix} ] , \end{matrix}

are matrices of size mn. Moreover, if G is minimal then also X of (6) is minimal and viceversa since the eigenvalues of X coincide with the eigenvalues of G filled up with zeros.

Unlike the former reduction, this technique can be used for reducing equation (1), for n = +∞, to solving (4) where the matrices A_{2}, A_{1} and A_{0} become semiinfinite, that is, their entries have subscripts ranging in the set ℕ of natural numbers. For n = +∞ we assume that ∑_{i=0}^{+∞} |A_{i}| has finite entries, where |A_{i}| is the matrix whose entries are the moduli of the entries of A_{i}. In fact, a direct inspection shows that if G is solution of (1) then the semi-infinite matrix

X = [ \begin{matrix}G & G^{2} & G^{3} & … \\ 0 & 0 & 0 & … \\ ⋮ & ⋮ & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋮ \end{matrix} ]  (7)

is solution of A_{0} + X A_{1} + X^{2}A_{2} = 0 where

\begin{matrix}A_{0} = [ \begin{matrix}A_{0} & 0 & … \\ 0 & 0 & ... \\ ⋮ & ... & ... \end{matrix} ] , & A_{1} = [ \begin{matrix}A_{1} & & & O \\ A_{2} & I \\ A_{3} & 0 & I \\ ⋮ & ⋮ & ... & ... \end{matrix} ] , \\ A_{2} = [ \begin{matrix}0 & −I & & O \\ & 0 & −I \\ & & ... & ... \\ O & & & ... \end{matrix} ] , \end{matrix}  (8)

are semi-infinite matrices. We observe also that if the spectral radius ρ(G) of G is less than 1, then X is a bounded operator in ℓ^{1}(ℕ), i.e. in the linear space of sequences {a_{i}}_{i∈ℕ} such that ∑_{i=0}^{+∞} |a_{i}| is finite. Moreover, since ∑_{i=0}^{+∞} |A_{i}| is convergent, then also A_{0}, A_{1} and A_{2} of (8) are bounded operators. Observe also that that the nonnull eigenvalues of X coincide with the nonnull eigenvalues of G. This implies that G is minimal if and only if X is a solution with eigenvalues of minimum modulus.

No reduction is known of (3) to (4). In [5] it is shown how (3) can be reduced to a system of quadratic equations.

3 Reduction to inverting a matrix Laurent polynomial

Let us consider the equation (1) for n = 2 and for a finite value of m, define A(z) = (z^{−1}A_{0} + A_{1} + zA_{2}). Assume also that a(z) = det A(z) has zeros ξ_{i}, i = 1, …, 2m such that

|ξ_{1}| ≤ … ≤ |ξ_{m}| < r_{1} < 1 < r_{2} < |ξ_{m+1}| ≤ … ≤ |ξ_{2m}|  (9)

for suitable r_{1}, r_{2}, so that a(z) is analytic in the annulus R = {z ∈ C : r_{1} ≤ |z| ≤ r_{2}}. If G is a solution such that ρ(G) = |ξ_{m}| then ξ_{1}, …, ξ_{m} are the eigenvalues of G and we have the block Wiener-Hopf factorization

\begin{matrix}A(z) = (I − z^{−1}G)(B_{0} − zB_{1}) = U(z^{−1})L(z) \\ B_{1} = −A_{2}, B_{0} = GA_{2} + A_{1}, \end{matrix}  (10)

where det U(z^{−1}) is analytic for |z| > r_{1}, and det L(z) is analytic for |z| < r_{2}. Equation (10) can be rewritten in matrix form as

[ \begin{matrix}A_{1} & A_{0} & & O \\ A_{2} & A_{1} & A_{0} \\ O & ... & ... & ... \end{matrix} ] = [ \begin{matrix}I & −G & & O \\ & I & −G \\ O & & ... & ... \end{matrix} ] [ \begin{matrix}B_{0} & & O \\ −B_{1} & B_{0} \\ O & ... & ... \end{matrix} ] = UL

(11)

where U and L are block upper triangular and lower triangular semi-infinite block Toeplitz matrices, respectively.

Proposition 1 Assume that the Laurent polynomial a(z) = det A(z) has zeros ξ_{i}, i = 1, …, 2m such that (9) holds so that it is analytic in the annulus R = {z ∈ C : r_{1} ≤ |z| ≤ r_{2}}. Assume that there exists a matrix G such that ρ(G) < 1 and

G^{2}A_{2} + GA_{1} + A_{0} = 0,

so that the eigenvalues of G coincide with the m zeros of a(z) of modulus less than r_{1}. Then the following properties hold

i) det B_{0} ≠ 0, for B_{0} = GA_{2} + A_{1}.

ii) There exist, bounded, the operators U^{−1}, L^{−1} of (11)

iii) There exist the Laurent matrix polynomial H(z) = ∑_{i=−∞}^{+∞} z^{i}H_{i} such that A(z)H(z) = I

iv) For S = B_{0}^{−1}B_{1} = −(GA_{2} + A_{1})^{−1}A_{2}, it holds ρ(S) < 1/r_{2} < 1,

A_{2} + A_{1}S + A_{0}S^{2} = 0,

the sum ∑_{i=0}^{+∞} G^{i}B_{0}^{−1}S^{i} is finite and

H_{i} = \{ \begin{matrix}H_{0}G^{−i}, & i < 0, \\ ∑_{i=0}^{+∞} G^{i}B_{0}^{−1}S^{i} & i = 0, \\ S^{i}H_{0} & i > 0. \end{matrix}

v) If the matrix equation

Y^{2}A_{0} + Y A_{1} + A_{2} = 0

has a solution R such that ρ(R) < 1/r_{2} < 1, then det H_{0} ≠ 0, det(Y A_{0} + A_{1}) ≠ 0 and W = −(Y A_{0} + A_{1})^{−1}A_{0} solves the equation

A_{2}W^{2} + A_{1}W + A_{0} = 0,

moreover we have the following representation of H_{i}

H_{i} = \{ \begin{matrix}H_{0}G^{−i} & = & W^{−i}H_{0}, & i ≤ 0, \\ H_{0}R^{i} & = & S^{i}H_{0}, & i ≥ 0 \end{matrix}

therefore, G = H_{0}^{−1}H_{−1}, R = H_{0}^{−1}H_{1}, S = H_{1}H_{0}^{−1}, W = H_{−1}H_{0}^{−1}.

Proof. From (10) we deduce that det A(z) = det(I − z^{−1}G) det(B_{0} − zB_{1}). Therefore, since G has eigenvalues of modulus less than 1 which coincide with ξ_{1}, …, ξ_{m}, then L(z) = det(B_{0} + zB_{1}) is singular only for z = ξ_{i}, i = m + 1, …, 2m, whence we deduce that B_{0} = L(0) is nonsingular and ρ(S) < 1/r_{2} < 1 for S = B_{0}^{−1}B_{1}. Moreover, the upper block triangular semi-infinite block Toeplitz matrix Z defined by its first block row (I, G, G^{2}, …) is such that ZU = UZ = I and, since ρ(G) < r_{1} < 1, Z is a bounded operator and Z = U^{−1}. A similar argument applies to the lower block triangular block Toeplitz semiinfinite matrix V defined by its first column (B_{0}^{−1}, SB_{0}^{−1}, S^{2}B_{0}^{−1}, …), for S = B_{0}^{−1}B_{1}, since it holds that V L = LV = I and ρ(S) < 1/r_{2} < 1. Part iii) follows with H(z) = (I − zS)^{−1}B_{0}^{−1}(I − z^{−1}G)^{−1} which is analytic in the annulus R. From the latter equation we deduce that

H_{k} = \{ \begin{matrix}S^{k}(∑_{i=0}^{+∞} S^{i}B_{0}^{−1}G^{i}) = S^{k}H_{0} & for k > 0, \\ ∑_{i=0}^{+∞} S^{i}B_{0}^{−1}G^{i} & for k = 0, \\ (∑_{i=0}^{+∞} S^{i}B_{0}^{−1}G^{i})G^{−k} = H_{0}G^{−k} & for k < 0. \end{matrix}  (12)

The sum H_{0} = ∑_{i=0}^{+∞} S^{i}B_{0}^{−1}G^{i} has entries with finite values since, for any operator norm ∥ · ∥ it holds ∥H_{0}∥ ≤ ∑_{i=0}^{+∞} ∥S^{i}∥ · ∥B_{0}^{−1}∥ · ∥G^{i}∥; moreover, from the relations between spectral radius and operator norms we know that for any ϵ > 0 there exists a norm such that ∥G∥ ≤ ρ(G) + ϵ and from the continuity of norms, since lim S^{i} = 0, we deduce that ∥S^{i}∥ is bounded from above by a constant γ. Whence we obtain that ∥H_{0}∥ ≤ γ∥B_{0}^{−1}∥ ∑_{i=0}^{+∞}(ρ(G) + ϵ)^{i} which is bounded from above by a constant if ϵ is such that ρ(G) + ϵ < 1. Concerning part v), if there exists a solution R of the equation R^{2}A_{0} + RA_{1} + A_{2} = 0 such that ρ(R) < 1 then, by following the same arguments as before, we deduce that Y A_{0} + A_{1} is nonsingular and the matrix W = −(Y A_{0} + A_{1})^{−1}A_{0} solves the equation A_{2}W^{2} + A_{1}W + A_{0}. In particular we obtain that H_{i} = H_{0}R^{i} for i ≥ 0 and H_{i} = W^{−i}H_{0} for i < 0, so that H(z) = H_{0}(I + ∑_{i=1}^{+∞}(z^{−i}G^{i} + z^{i}R^{i})) = (I + ∑_{i=1}^{+∞}(z^{−i}W^{i} + z^{i}S^{i}))H_{0}. The singularity of H_{0} would imply that det H(z) = 0 for any z ∈ R which is absurd since for z ∈ R it holds det H(z) = 1/ det A(z) ≠ 0. ⊓ ⊔

By using the equivalence of norms over a finite dimensional space and using the relation between operator norms and spectral radius, we may deduce that for any operator norm ∥ · ∥ there exists a constant γ > 0 such that

∥H_{i}∥ ≤ \{ \begin{matrix}∥H_{0}∥ · ∥G^{−i}∥ ≤ γr_{1}^{−i} & i ≤ 0 \\ ∥H_{0}∥ · ∥R^{i}∥ ≤ γr_{2}^{−i} & i ≥ 0 \end{matrix}  (13)

that is, the coefficients of H_{i} decay exponentially to zero. This nice property allows us to look at H(z) as to a Laurent polynomial plus some truncation error

H(z) = ∑_{i=−N/2}^{N/2−1} z^{i}H_{i} + E(z),  (14)

where the norm of the block coefficients of E(z) is O(ρ^{N}) for ρ = max{r_{1}, r_{2}^{−1}}. If N is “large enough” the coefficients H_{−N/2}, …, H_{N/2−1} can be approxiamted by means of the coefficients of the Laurent matrix polynomial

H̃(z) = ∑_{i=−N/2}^{N/2−1} z^{i}H̃_{i}

which interpolates H(z) at the N = 2^{d} roots of 1, i.e., such that

\{ \begin{matrix}H̃(ω_{N}^{j}) = H(ω_{N}^{j}), \\ & j = 0, 1, …, N − 1 \\ H(ω_{N}^{j}) = [ω_{N}^{−j}A(ω_{N}^{j})]^{−1} \end{matrix}  (15)

where ω_{N} = cos(2π/N) + i sin(π/N) and i^{2} = −1.

In fact, we have the following estimate for the infinity norm ∥H_{i} − H̃_{i}∥_{∞}.

Proposition 2 For the coefficients H̃_{i} of the Laurent polynomial H̃(z) satisfying (15) it holds ∥H_{i} − H̃_{i}∥_{∞} ≤ max_{j} ∥R(ω_{N}^{j})∥_{∞}. Therefore ∥H_{i} − H̃_{i}∥_{∞} ≤ αρ^{N/2} for a positive constant α and for ρ = max{r_{1}, r_{2}^{−1}}.

Proof. From (14) it holds ∑_{j=0}^{N−1}(H̃_{j} − H_{j})ω_{N}^{ij} = E(ω_{N}^{i}), whence (H̃_{j} − H_{j}) = \frac{1}{N} ∑_{i=0}^{N−1} ω_{N}^{−ij}E(ω_{N}^{−i}). Taking the infinity norm of both sides of the latter equation yields |∥H̃_{j} − H_{j}∥_{∞} ≤ max_{i} ∥E(ω_{N}^{−1})∥_{∞}. Since E(z) = ∑_{i≥N/2} z^{i}H_{i} + ∑_{i≤−N/2−1} z^{i}H_{i}, from (13) we conclude that there exists a constant α such that ∥E(ω_{N}^{i})∥_{∞} ≤ αρ^{N/2} for ρ = max{r_{1}, r_{2}^{−1}}. ⊓ ⊔

3.1 The case of semi-infinite blocks

The case where the blocks A_{0}, A_{1}, A_{2} are semi-infinite is more tricky as the following example shows. Let Z = (z_{i,j}), z_{i+1,i} = 1, z_{i,j} = 0 for i + 1 ≠ j be a semi-infinite matrix and let A_{0} = −2I, A_{1} = Z, A_{2} = Z^{2}. Since Z^{T} Z = I, then X = Z^{T} solves the equation A_{0} + XA_{1} + X^{2}A_{2} = 0. On the other hand the spectrum of X is S_{X} = {λ ∈ C : |λ| ≤ 1}, whereas the spectrum of A_{0} + λA_{1} + λ^{2}A_{2} is S = {λ ∈ C : |λ| ≥ 1}. Therefore λ ∈ S_{X} does not necessarily imply that λ ∈ S.

However, under suitable assumptions, we may arrive at a characterization of the minimal solution of the quadratic equation A_{0} + XA_{1} + X^{2}A_{2} = 0 in terms of the blocks H_{1} and H_{−1} of the Laurent series H(z) = ∑_{i∈ℤ} z^{i}H_{i} such that H(z)A(z) = A(z)H(z) = I.

Let R = {z ∈ C : r_{1} ≤ |z| ≤ r_{2}} for given r_{1} < 1 and r_{2} > 1, denote ρ = max{r_{1}, r_{2}^{−1}}. Assume that A_{0}, A_{1}, A_{2} are semi-infinite bounded operators in ℓ^{2}(ℕ) such that the operator Laurent polynomial A(z) = z^{−1}A_{0} + A_{1} + zA_{2} is analytic and invertible with bounded inverse for z ∈ R. In this hypothesis [19] there exists an operator Laurent series H(z) = ∑_{i∈ℤ} z^{i}H_{i} which is analytic and bounded for z ∈ R such that H(z)A(z) = A(z)H(z) = I; moreover it holds

H_{k} = \frac{1}{2πı} ∫_{|z|=r} z^{−1−k}A(z)^{−1}dz

where r_{1} ≤ r ≤ r_{2}. So that, for the 2-norm ∥ · ∥, there exists a positive constant α such that

∥H_{k}∥ ≤ αρ^{|k|}, for k ∈ ℤ.  (16)

Assume also that the operator represented by the semi-infinite matrix

K = [ \begin{matrix}A_{1} & A_{0} & & O \\ A_{2} & A_{1} & A_{0} \\ O & ... & ... & ... \end{matrix} ]

is invertible and there exists a semi-infinite matrix G which solves the equation G^{2}A_{2} + GA_{1} + A_{0} = 0 and such that ∥G^{k}∥ ≤ βg^{k} for a positive constant β and a positive g < 1. Then the following factorization holds

K = [ \begin{matrix}I & −G & & O \\ & I & −G \\ O & & ... & ... \end{matrix} ] [ \begin{matrix}B_{0} & & O \\ −B_{1} & B_{0} \\ O & ... & ... \end{matrix} ] ,  (17)

where B_{0} = GA_{2} + A_{1}, B_{1} = −A_{2}. We have the following

Lemma 1 The inverse of K can be represented as

K^{−1} = [ \begin{matrix}H_{0} & H_{−1} & H_{−2} & … \\ H_{1} & H_{0} & H_{−1} & ... \\ ⋮ & ... & ... & ... \end{matrix} ] − [ H_{1} ; H_{2} ; ⋮ ] [ G G^{2} …] .

Proof. Multiply the above expression to the right by K, use the fact that [G, G^{2}, …]K = [−A_{0}, 0, …] and that H_{i+1}A_{0} + H_{i}A_{1} + H_{i−1}A_{2} is zero for i ≠ 0 and is the identity matrix for i = 0, and obtain that the result is the identity matrix. ⊓ ⊔

From the factorization (17) we find that B_{0} is invertible and, setting S = B_{0}^{−1}B_{1} we have

K^{−1} [ \begin{matrix}I & −G & & O \\ & I & −G \\ O & & ... & ... \end{matrix} ] [ \begin{matrix}B_{0}^{−1} & & O \\ & B_{0}^{−1} \\ O & & ... \end{matrix} ] = [ \begin{matrix}I & & & O \\ S & I \\ S^{2} & S & I \\ ⋮ & ... & ... & ... \end{matrix} ]

that is, in the light of Lemma 1, we get S^{k} = H_{k}B_{0}^{−1} − H_{k+1}GB_{0}^{−1}. Whence, from (16) it follows that ∥S^{k}∥ ≤ δρ^{k} for a positive δ. This implies the invertibility of I − zS and from the relation H(z) = (I − zS)^{−1}B_{0}^{−1}(I − z^{−1}G)^{−1} we deduce that (12) holds.

The invertibility of H_{0} follows from the invertibility of Φ(z) = I+∑_{i=1}^{+∞}(z^{−i}G^{i}+ z^{i}R^{i}) in R. In order to prove the invertibility of Φ(z) we observe that multiplying Φ(z) on the left by I − Rz we get

(I − Rz)Φ(z) = (I − RG) ∑_{i=0}^{+∞} G^{i}z^{−i} = (I − RG)(I − Gz)^{−1}.

Therefore, Φ(z) = (I − Rz)^{−1}(I − RG)(I − Gz^{−1})^{−1} is invertible as product of invertible operators.

The invertibility of H_{0} allows us to obtain the following representations of the semi-infinite matrices G and R

\begin{matrix}G = H_{0}^{−1}H_{−1} \\ R = H_{0}^{−1}H_{1}. \end{matrix}

Similar arguments can be applied for the solutions S and W of the matrix equations A_{2} + A_{1}S + A_{0}S^{2} = 0 and A_{2}W^{2} + A_{1}W + A_{0} = 0.

4 Algorithms

From Propositions 1 and 2 of section 3 we obtain the following algorithm for computing a minimal solution of the equation (1) for n = 2.

Algorithm 1 Solution of a quadratic matrix equation.

Input: The m × m matrices A_{0}, A_{1}, A_{2} such that the Laurent polynomial A(z) = z^{−1}A_{0} + A_{1} + zA_{2} has m zeros of modulus less than 1 and m zeros of modulus greater than 1.

Output: An approximation G̃ of a minimal solution G of the matrix equation A_{0} + XA_{1} + X^{2}A_{2} = 0 if it exists.

Computation:

Stage 1. Let N a sufficiently large integer and denote ω_{N} a primitive N-th root of 1. Compute W_{i} = ω_{N}^{−i}A(ω_{N}^{i}), i = 0, …, N − 1.

Stage 2. Compute V_{i} = W_{i}^{−1} for i = 0, …, N − 1.

Stage 3. Compute H̃_{i} = \frac{1}{N} ∑_{j=0}^{N−1} ω_{N}^{−ij}V_{j}, i = −1, 0.

Stage 4. Output the approximation X̃ = H_{0}^{−1}H_{−1}.

It is immediate to verify that the cost of the above algorithm is γ_{1}m^{3}N + γ_{2}m^{2}N for two suitable constants γ_{1} and γ_{2}. Moreover, due to Proposition 2 it is sufficient to chose N = O(log ϵ^{−1}/ log ρ) for ρ = max{r_{1}, 1/r_{2}} in order to have an approximation error at most ϵ in the solution.

Algorithm 1 can be implemented with a heuristic strategy for the automatic evaluation of the value of the integer N. In fact, it is possible to apply the above algorithms with different values of N, say, N = N_{0}, 2N_{0}, 4N_{0}, …, until the difference of two consecutive approximations of the solution has norm less than a given error bound ϵ. The quantities computed by the algorithm applied with N = 2^{k}N_{0} can be used in the computation with N = 2^{k+1}N_{0}.

Remark 1 On the case of structured blocks. Observe that if the quadratic equation is obtained through the reduction of the general equation (1) to (4), then the blocks A_{i} are structured. Therefore, the costs of stage 1 of the above algorithm is substantially reduced. In fact, these blocks belong to a linear space whose dimension is much lower than the square of the size of the blocks and the linear combinations involved in stage 1 can be simply restricted to a lower dimensional basis. We can get computational advantages even in stage 2 and 3 if the structure is somehow kept under inversion. A detailed analysis of this case is performed in section 5.1 concerning an application to the numerical solution of certain Markov chains.

Remark 2 On the case of semi-infinite blocks. In principle the above algorithm can be applied even in the case where m = +∞ and the matrices have the structure (8). In fact, in this case in order to compute G we have to compute only the leading principal submatrix of size m × m of X of (7). For this purpose we have to compute all the entries of H_{0} and the first m columns of H_{−1}. A crucial issue is the inversion of the semi-infinite matrices W_{i} = A_{0}ω_{N}^{−i} + A_{1} + A_{2}ω_{N}^{i}, i = 0, 1, …, N, which is needed at stage 2 as well as the computation of H_{0}^{−1}. Indeed, these computations must be done in approximate way and the analysis of the conditions on A(z) in order that the inverse of a finite section of W_{i} (or of H_{0}) provides the seeked approximation are under investigation.

4.1 Using Graeffe’s iteration

Assume that m is finite and consider the following sequence of Laurent matrix polynomials {A^{(i)}(z)}_{i∈ℕ}

\begin{matrix}A^{(0)}(z) = A(z) = z^{−1}A_{0} + A_{1} + z^{2}A_{2} \\ A^{(i+1)}(z^{2}) = A^{(i)}(z)[A_{1}^{(i)}]^{−1}A^{(i)}(−z), i = 0, 1, … \end{matrix}  (18)

Under the assumptions of Proposition 1, if the matrices A_{1}^{(i)}, i = 1, 2, …, are nonsingular then the following properties hold [4]

G^{2^{i}} solves A_{0}^{(i)} + XA_{1}^{(i)} + X^{2}A_{2}^{(i)} = 0
R^{2^{i}} solves X^{2}A_{0}^{(i)} + XA_{1}^{(i)} + A_{2}^{(i)} = 0
∥A_{0}^{(i)}∥ = O(r_{1}^{2^{i}}), ρ(G) < r_{1} < 1
∥A_{2}^{(i)}∥ = O(r_{2}^{−2^{i}}), ρ(R) < r_{2}^{−1} < 1
∥A_{1}^{(i)}∥, ∥[A_{1}^{(i)}]^{−1}∥ ≤ γ, γ > 0
lim_{i} A^{(i)}(z)[A_{1}^{(i)}]^{−1} = I

Moreover we have the following identity which can be proved by induction on i

\begin{matrix}A^{(i+1)}(z^{2^{i+1}})[A_{1}^{(i+1)}]^{−1} = A(z)H^{(i)} \\ H^{(i)} = \{[A_{1}]^{−1}A(−z)[A_{1}^{(1)}]^{−1}A^{(1)}(−z^{2})[A_{1}^{(2)}]^{−1}A^{(2)}(−z^{4}) ⋯ \\ A^{(i)}(−z^{2^{i}})\} [A_{1}^{(i+1)}]^{−1} \end{matrix}  (19)

whence we find that H^{(i)} converges to H(z) double exponentially.

The following algorithm approximates the minimal solution G of the equation G^{2}A_{2} + GA_{1} + A_{0} = 0 by means of the relation G = H_{0}^{−1}H_{−1} where H_{0} and H_{−1} are computed by means of Graeffe’s iteration.

Algorithm 2 Input: The m × m matrices A_{0}, A_{1}, A_{2} such that the assumptions of Proposition 1 are satisfied; a positive real ϵ.

Output: An approximation of the minimal solution G of the equation G^{2}A_{2} + GA_{1} + A_{0} = 0.

Computation:

Stage 1. (Graeffe iteration) Compute the Graeffe sequence A^{(i)}(z) i = 1, …, k of (18) until ∥A_{0}^{(k)}∥ ≤ ϵ, ∥A_{2}^{(k)}∥ ≤ ϵ.

Stage 2. (Back substitution)

\begin{matrix}B^{(k)}(z) = A^{(k)}(z) \\ B^{(i−1)}(z) = A^{(i−1)}(z)[A_{1}^{(i−1)}]^{−1}B^{(i)}(z^{2}), i = k, …, 1 \end{matrix}

Stage 3. Output [B_{0}^{(0)}]^{−1}B_{−1}^{(0)}.

Observe that in the computation of B^{(i)}(z) of the above algorithm only the blocks B_{−1}^{(i)}, B_{0}^{(i)}, B_{1}^{(i)} must be computed. The cost per step is just one matrix inversion, 10 matrix multiplications 4 matrix additions.

This algorithm is strictly related with the one introduced in [6] which is based on the cyclic reduction technique. Both the algorithms compute the coefficients of the Graeffe sequence. Algorithm 1 performs the back substitution stage separately whereas the algorithm of [6] performs the back substitution stage simultaneously with Graeffe’s iteration.

Remark 3 On the case of semi-infinte blocks. In principle, the above algorithm can be applied even in the case m = +∞. A crucial issue is the approximation of the inverses of the semi-infinite matrices A_{1}^{(i)} in the Graeffe sequence, with A_{1}^{(0)} = A_{1}, as well as the inversion of B_{0}^{(0)}. The analysis of conditions of A(z) under which this approximation can be performed by means of finite sections is under investigation.

5 The case of Markov chains

For matrix equations arising from the solution of Markov chains there is an additional difficulty. Namely, the Laurent polynomial det A(z) is zero for z = 1 so that the assumption of Proposition 1 are not satisfied anymore. More precisely, the zeros of det A(z) are

|ξ_{1}| ≤ ⋯ ≤ ξ_{m} = 1 < |ξ_{m+1}| ≤ ⋯ ≤ |ξ_{2m}|

so that ρ(G) = 1, moreover, e^{T} G = e^{T}, e^{T} = (1, …, 1). This drawback can be easily overcome by applying a simple idea introduced in [15] for accelerating the convergence rate.

Let u ∈ ℝ^{m}, e^{T} u = 1 and set V = G − ue^{T}. The eigenvalues of V are ξ_{1}, …, ξ_{m−1}, 0, in fact, e^{T} V = e^{T} G − e^{T} = 0, and if w ≠ 0 is a right eigenvector of G, i.e., such that Gw = λw, λ ≠ 1, then w is orthogonal to the left eigenvector e, i.e., w^{T} e = 0. Therefore V w = Gw − ue^{T} w = Gw = λw, that is, λ is eigenvalue of V. Moreover it holds G^{i} = V^{i} + ∑_{j=0}^{i−1} V^{j}ue^{T}. Thus the equation A_{0} + GA_{1} + G^{2}A_{2} = 0 is reduced to

Â_{0} + V Â_{1} + V^{2}Â_{2} = 0,  (20)

where

Â_{0} = A_{0} + ue^{T} (A_{1} + A_{2}),
Â_{1} = A_{1} + ue^{T} A_{2},
Â_{2} = A_{2}.

Furthermore det(Â_{0} + zÂ_{1} + z^{2}Â_{2}) has m zeros of modulus < 1 and m zeros of modulus > 1 so that the hypotheses of Proposition 1 are satisfied by the equation (20) and we may apply the machinery of section 4 in order to approximate V and then recover G as G = V + ue^{T}.

5.1 Non-Skip-Free Markov chains

A special case encountered in queueing problems is the class of Non-Skip-Free Markov chains [10] where we are given integers n, m, k and m × m matrices A_{0}, …, A_{n} such that the following properties hold

A_{k}+I ≥ 0, A_{i} ≥ 0, i ≠ k, and ∑_{i=0}^{n} A_{i}+I is stochastic, that is e^{T} ∑_{i=0}^{n} A_{i} = 0.

For A(z) = z^{−k} ∑_{i=0}^{n} z^{i}A_{i}, the Laurent polynomial det A(z) has

mk − 1 zeros of modulus < 1 one zero equal to 1 (n − k)m zeros of modulus > 1.

Moreover, there exists the Wiener-Hopf factorization

A(z) = (∑_{i=0}^{k} z^{−i}U_{i})(∑_{i=0}^{n−k} z^{i}L_{i})  (21)

where U_{i} and L_{i} are m × m matrices such that U_{0} = I, det(∑_{i=0}^{k} z^{−i}U_{i}) has mk zeros of modulus ≤ 1, det(∑_{i=0}^{n−k} z^{i}L_{i}) has m(n − k) zeros of modulus > 1.

Here the problem is to compute the matrices U_{0}, …, U_{k−1} in (21).

Without loss of generality we may assume that the Laurent polynomial det A(z) is nonzero for z = 1. In fact, similarly to the case of a quadratic equation, we may apply a deflation technique for removing the unit zero of det A(z). We show the generalization of this technique below.

Let u be any vector such that u^{T} e = 1 and consider the matrix E = ue^{T} such that E^{i} = E, i ≠ 0. Define

E(z) = ∑_{i=0}^{+∞} z^{i}E^{i} = I + E ∑_{i=1}^{+∞} z^{i}

and observe that E(z) = (I − zE)^{−1}. Since e^{T} ∑_{i=0}^{n} A_{i} = 0, it holds that E ∑_{i=0}^{n} A_{i} = 0. Therefore the function Â(z) = E(z^{−1})A(z) is a matrix Laurent polynomial, i.e.,

\begin{matrix}Â(z) = z^{−k} ∑_{i=0}^{n} Â_{i}z^{i}, \\ Â_{i} = A_{i} + E(∑_{j=i+1}^{n} A_{j}), i = 0, …, n. \end{matrix}  (22)

We may also write

A(z) = (I − z^{−1}E)Â(z).

Now, for the sake of simplicity assume that A_{0} and A_{n} are nonsingular (the case where det A_{0} = 0 or det A_{n} = 0 can be treated by means of a continuity argument). Since Â_{n} = A_{n} the Laurent polynomials det Â(z) and det A(z) have the same degree, and, therefore the same number of roots.

Since the matrix I − z^{−1}E is defined for z ≠ 0 and is singular only for z = 1, we deduce that if λ is zero of det Â(z) and λ ≠ 0, 1, then it is zero of det A(z) and viceversa. Moreover, det Â(0) = 0 since e^{T} Â(0) = eÂ_{0} = e(∑_{i=0}^{n} A_{0}) = 0. Thus, we obtain that det Â(z) has the same zeros of det A(z) except for z = 1 which is replaced in det Â(z) by z = 0.

Now observe that multiplying equation A(z) = U(z^{−1})L(z) to the left by E(z) yields

Â(z) = Û(z^{−1})L(z)
Û(z) = E(z)U(z) = ∑_{i=0}^{k} z^{i}Û_{i}
Û_{i} = U_{i} + E ∑_{j=1}^{i−1} U_{i}, i = 0, …, k.

Therefore, in order to compute U_{0}, …, U_{k} of the factorization (21), we may first evaluate the blocks Â_{i}, i = 0, …, n by means of (22), then compute the Wiener-Hopf factorization Â(z) = Û(z^{−1})L(z) of Â(z) and recover the blocks U_{i}, i = 0, …, k from the blocks Û_{i}, i = 0, …, k by means of the equation U(z) = E(z)^{−1}Û(z), i.e.,

\begin{matrix}U_{0} = Û_{0} \\ U_{i} = Û_{i} − EU_{i−1}, i = 1, …, k. \end{matrix}

Observe that the Laurent polynomial det Â(z) is nonzero for z = 1. Therefore it is no loss of generality to assume that the zeros ξ_{i} of det A(z) are such that

|ξ_{1}| ≤ … ≤ |ξ_{mk}| < 1 ≤ |ξ_{mk+1}| ≤ … ≤ |ξ_{m}|.  (23)

Now let us rewrite A(z) = U(z^{−1})L(z) in matrix form as

\begin{matrix}[ \begin{matrix}A_{k} & A_{k−1} & … & A_{0} \\ A_{k+1} & A_{k} & A_{k−1} & ... & A_{0} \\ ⋮ & ... & ... & ... & ... & ... \\ A_{n} & ... & ... & ... & ... & ... & ... \\ & ... & ... & ... & ... & ... & ... \end{matrix} ] = \\ [ \begin{matrix}U_{0} & U_{1} & … & U_{k} \\ & U_{0} & U_{1} & ... & U_{k} \\ & & ... & ... & ... & ... \end{matrix} ] & [ \begin{matrix}L_{0} \\ L_{1} & L_{0} \\ ⋮ & L_{1} & L_{0} \\ L_{n−k} & ... & ... & ... \\ & ... & ... & ... \end{matrix} ] \end{matrix}  (24)

If we partition the matrices in (24) into qm×qm blocks where q = max{k, n− k} we obtain the equivalent equation

[ \begin{matrix}A_{1} & A_{0} & & O \\ A_{2} & A_{1} & A_{0} \\ & ... & ... & ... \\ O & & ... & ... \end{matrix} ] = [ \begin{matrix}I & −G & & O \\ & I & −G \\ O & & ... & ... \end{matrix} ] [ \begin{matrix}B_{0} & & O \\ −B_{1} & B_{0} \\ O & ... & ... \end{matrix} ]

where

\begin{matrix}A_{0} = [ \begin{matrix}A_{k−q} & & O \\ ⋮ & ... \\ A_{k−1} & … & A_{k−q} \end{matrix} ] , A_{1} = [ \begin{matrix}A_{k} & … & A_{k−q+1} \\ ⋮ & ... & ⋮ \\ A_{k+q−1} & … & A_{k} \end{matrix} ] , \\ A_{2} = [ \begin{matrix}A_{k+q} & … & A_{k+1} \\ ⋮ & ... & ⋮ \\ A_{k+q−1} & … & A_{k+q} \end{matrix} ] \end{matrix}

and

G = − [ \begin{matrix}U_{q} & & O \\ ⋮ & ... \\ U_{1} & … & U_{q} \end{matrix} ] [ \begin{matrix}U_{0} & … & U_{q−1} \\ & ... & ⋮ \\ O & & U_{0} \end{matrix} ]^{−1} ,

B_{0} = [ \begin{matrix}U_{0} & … & U_{q−1} \\ & ... & ⋮ \\ O & & U_{0} \end{matrix} ]^{−1} [ \begin{matrix}L_{0} & & O \\ ⋮ & ... \\ L_{q−1} & … & L_{0} \end{matrix} ] ,

B_{1} = − [ \begin{matrix}U_{0} & … & U_{q−1} \\ & ... & ⋮ \\ O & & U_{0} \end{matrix} ]^{−1} [ \begin{matrix}L_{q} & … & L_{1} \\ ⋮ & ... & ⋮ \\ L_{2q−1} & … & L_{q} \end{matrix} ] .

and we assume that A_{i} = 0, U_{i} = 0, L_{i} = 0 if i is out of range.

A simple way for solving this computational problem is to rewrite it in terms of matrix equation as

A_{0} + GA_{1} + G^{2}A_{2} = 0

and to apply the algorithms of section 4 for its solution under the assumption (23).

It is important to point out that it is enough to compute the first block column of the matrix G which provides the seeked solution U_{1}, …, U_{k}.

For the algorithm 1 based on the evaluation/interpolation technique, we have the following simplifications.

Due to the structure of A_{0}, A_{1}, A_{2}, the evaluation of W_{i} = ω_{N}^{−i}A_{0}+A_{1}+ω_{N}^{i}A_{2} for i = 0, …, N − 1 at stage 1 can be restricted to the first m columns and to the first m rows so that the cost is O(Nm^{2}q) ops.

Since the matrices W_{i}, i = 0, …, N are ω_{N}^{i}-block circulant matrices [7], then their inversion at stage 2 can be performed by means of FFTs and costs O(Nm^{2}q log q + Nm^{3}) ops. We recall that the set of α-block circulant matrices is the algebra formed by matrices of the kind ∑_{i=0}^{q−1} C(α)^{i} ⊗ W_{i}, where

C(α) = [ \begin{matrix}0 & 1 & 0 & … & 0 \\ 0 & ... & ... & ... & ⋮ \\ ⋮ & ... & ... & ... & 0 \\ 0 & ... & ... & 0 & 1 \\ α & 0 & … & … & 0 \end{matrix} ] ,

⊗ denotes the Kronecker product and W_{i} are m × m matrices.

Since ω_{N}^{i} block circulant matrices are an algebra made up by block Toeplitz matrices, then the matrices V_{i} = W_{i}^{−1} are block Toeplitz and also the matrices H̃_{i} at stage 4 are block Toeplitz. Hence, their computation costs O(Nm^{2}q) ops.

Finally, at stage 4 we have to compute only the first m columns of X̃ and this computation can be accomplished by solving m block Toeplitz systems. The cost of this part is independent of N.

The overall cost of Algorithm 1 amounts to O(N(pm^{3} + m^{2}q log q)), where we ignore the part independent of N, which for k = O(n) reduces to O(km^{3} + km^{2} log k) This approach compares favourably with the customary techniques based on functional iterations [10] where the computational cost is O(m^{3}k^{2}) ops per step. The acceleration that we obtain this way is by a factor of k.

Similar accelerations can be obtained if the Graeffe iteration is applied by using the Toeplitzness of the blocks H_{−1}, H_{0}, H_{1}.

Another possibility of dealing with Non-Skip-Free Markov chains is to apply the reduction based on (8). The structure analysis of the matrices A_{0}, A_{1}, and A_{2} obtained in this way has not yet been done and is object of future research.

The last remark of this section concern the case where det A(z) has more than one zero of modulus 1. If ξ ≠ 1 is such that |ξ| = 1 and det A(ξ) = 0, and if we know a vector w such that w^{T} A(ξ) = 0, then we find that ŵÂ(ξ) = 0, for ŵ = w^{T} (I − ξ^{−1}E), therefore we may apply the same deflation technique in order to deflate the root ξ from the matrix Laurent polynomial Â(z). That is, it is enough to chose a vector v such that v^{T} ŵ = 1, define the matrix E_{1} = vw̃^{T} and the function E_{1}(z) = (I − zE_{1})^{−1} and finally consider the deflated matrix Laurent polynomial E_{1}(z^{−1})Â(z).

6 Conclusion

We have seen that different classes of matrix equations can be reduced to solving a quadratic matrix polynomial equation. The computation of the minimal solution of a quadratic matrix polynomial equation is reduced to inverting a matrix Laurent polynomial. Two algorithms for this computation have been introduced. One is based on the evaluation/interpolation technique at the roots of 1, the second one is based on applying the Graeffe iteration to matrix Laurent polynomials. An application to solving problems related to the numerical solution of Markov chains described by a matrix Laurent polynomial A(z) has been shown. A technique for deflating the unwanted unit zero of det A(z) has been presented and extended for the computation of the Wiener-Hopf factorization in the Non-Skip-Free model.

References

[1] W. N. Anderson Jr., T. D. Morley, and G. E. Trapp. Positive solutions to X = A − BX^{−1}B^{∗}. Linear Algebra Appl., 134:53–62, 1990.

[2] S. Barnett. Polynomials and Linear Control Systems. Number 77 in Textbooks in Pure and Applied Mathematics. Marcel Dekker, Inc., New York, 1983.

[3] D. A. Bini and A. Böttcher. Polynomial factorization through Toeplitz matrix computations. To appear on Linear Algebra Appl.

[4] D. A. Bini, L. Gemignani, and B. Meini. Computations with infinite Toeplitz matrices and polynomials. To appear on Linear Algebra Appl.

[5] D. A. Bini, G. Latouche, B. Meini, Solving nonlinear matrix equations arising in Tree-Like stochastic processes, submitted to Linear Algebra Appl.

[6] D. A. Bini and B. Meini. On the solution of a nonlinear matrix equation arising in queueing problems. SIAM J. Matrix Anal. Appl., 17:906–926, 1996.

[7] D. A. Bini, B. Meini, Solving block banded block Toeplitz systems with structured blocks: algorithms and applications, in Structured Matrices: recent Developments in Theory and Computation, D.A. Bini, E. Tyrtyshnikov and P. Yalamov, Editors, Nova Science Publisher Inc., New York 2001.

[8] D. A. Bini, B. Meini, and V. Ramaswami. Analyzing M/G/1 paradigms through QBDs: the role of the block structure in computing the matrix G. In G. Latouche and P. Taylor, editors, Advances in Algorithmic Methods for Stochastic Models, pages 73–86. Notable Publications, New Jersey, USA, 2000. Proceedings of the Third Conference on Matrix Analytic Methods.

[9] A. Ferrante and B. C. Levy. Hermitian solutions of the equation X = Q + NX^{−1}N^{∗}. Linear Algebra Appl., 247:359–373, 1996.

[10] H. R. Gail, S. L. Hantler, and B. A. Taylor. Non-skip-free M/G/1 and G/M/1 type Markov chains. Adv. Appl. Prob., 29:733–758, 1997.

[11] H. R. Gail, S. L. Hantler, and B. A. Taylor. Use of characteristics roots for solving infinite state Markov chains. In W. K. Grassmann, editor, Computational Probability, pages 205–255. Kluwer Academic Publishers, 2000.

[12] N. J. Higham and H.-M. Kim. Numerical analysis of a quadratic matrix equation. IMA J. of Numerical Analysis, 20:499–519, 2000.

[13] G. Latouche and V. Ramaswami. Introduction to Matrix Analytic Methods in Stochastic Modeling. ASA-SIAM Series on Statistics and Applied Probability 5. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1999.

[14] B. Meini. Efficient computation of the extreme solutions of X +A^{∗}X^{−1}A = Q, X − A^{∗}X^{−1}A = Q. Math. Comp., 2001. To appear.

[15] C. He, B. Meini, N. H. Rhe, A shifted cyclic reduction for quasi birth-death problems, SIAM J. Matrix Anal. Appl., to appear.

[16] M. F. Neuts. Structured Stochastic Matrices of M/G/1 Type and Their Applications. Dekker Inc., 1989.

[17] M. F. Neuts. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach. The Johns Hopkins University Press, Baltimore, MD, 1981.

[18] V. Ramaswami, The generality of Quasi Birth-and-Death processes, in Advances in Matrix Analytic Methods for Stochastic Models, A. Alfa and S. Chakravarthy Editors, pp. 93–113, Notable Publications Inc., NJ 1998.

[19] A. E. Taylor and D. C. Lay, Introduction to Functional Analysis, John Wiley& Sons, New York 1980.

[20] X. Zhan and J. Xie. On the matrix equation X + A^{T} X^{−1}A = I. Linear Algebra Appl., 247:337–345, 1996.

[21] X. Zhan. Computing the extremal positive definite solutions of a matrix equation. SIAM J. Sci. Comput., 17:1167–1174, 1996.