Result: Shift techniques and canonical factorizations in the solution of M/G/1-type Markov chains

□ By using properties of canonical factorizations, we prove that under very mild assumptions, the shifted cyclic reduction method (SCR) can be applied for solving QBD problems with no breakdown and that it always converges. For general M/G/1 type Markov chains we prove that SCR always converges if no breakdown is encountered. Numerical experiments showing the acceleration provided by SCR versus cyclic reduction are presented.

AN0017539968;i3f01may.05;2019Feb27.14:25;v2.2.500

Shift Techniques and Canonical Factorizations in the Solution of M/G/1-Type Markov Chains. 

By using properties of canonical factorizations, we prove that under very mild assumptions, the shifted cyclic reduction method (SCR) can be applied for solving QBD problems with no breakdown and that it always converges. For general M/G/1 type Markov chains we prove that SCR always converges if no breakdown is encountered. Numerical experiments showing the acceleration provided by SCR versus cyclic reduction are presented.

Keywords: Canonical factorization; Cyclic reduction; Markov chains; Matrix equations; Shift technique; Primary 15A24, 60J22; Secondary 47B35, 60K25

1. INTRODUCTION AND PRELIMINARIES

Matrix equations of the kind

Graph

where A<subs>i</subs>, for i = − 1, 0, 1,..., are nonnegative m×m matrices such that A<subs>i</subs> is row stochastic, are encountered in the solution of Markov chains of the M/G/1 and G/M/1 type which model large parts of queuing problems (Neuts[12].

For this kind of equation the main interest is to compute the minimal nonnegative solutionG of (1), that is, such that 0 ≤ G ≤ X for any other non-negative solution X. Here, matrix inequality is meant component-wise. A nice feature of this kind of problems is that G is also the spectral minimal solution of (1); that is ρ(G) ≤ ρ(X) for any other solution X.

The design of efficient algorithms for solving (1) is a crucial step in the solution of queuing models arising in applications. Functional iterations provide a reliable tool for many problems, however, in many circumstances, their slow (linear or even sublinear) convergence makes this class of algorithms almost useless because of the huge number of iterations needed to provide reasonable approximations. In critical, or close to critical cases, quadratically convergent algorithms like Logarithmic Reduction and Cyclic Reduction (CR) are very effective means to overcome these difficulties (Latouche[11]; Bini[5]). In particular, for positive recurrent or transient Markov chains, they provide a matrix sequence {X<sups>(n)</sups>}<subs>n≥0</subs>, which quadratically converges to G, while in the null recurrent case the convergence turns to linear (Guo[8][9]).

To be more specific, consider the M/G/1 case modelled by the stochastic matrix

Graph

where B<subs>i</subs> ≥ 0 and B<subs>i</subs> is row stochastic, and observe that the matrix valued function A(z) = z<sups>i+1</sups>A<subs>i</subs> associated with this Markov chain is analytic for |z| < 1 and continuous for |z| ≤ 1.

Here, we make slightly stronger assumptions, which are usually satisfied in the applications. They are listed below together with the non-negativity properties of the matrices A<subs>i</subs> and the stochasticity of A<subs>i</subs>.

Assumption 1

The matrix function

Graph

where A<subs>i</subs> are non-negative m×m matrices such that A<subs>i</subs> is row stochastic, is analytic for |z| < r for some r > 1, and the associated Markov chains are positive recurrent. Moreover, there exists a zero of modulus greater than 1 of det(zI − A(z)).

Assumption 2

The Markov chain on ℤ×{1, 2,..., m} with transition matrix A = (A<subs>j−i</subs>)<subs>i, j∈ℤ</subs>, where we assume A<subs>i</subs> = 0 if i < − 1, has only one final class ℤ×&#55349;&#56494;, where &#55349;&#56494; ⊆ {1, 2,..., m}. Every other state is on a path to the final class.

Under the above assumptions, the zeros ξ<subs>i</subs>, i = 1, 2,..., of det(zI − A(z)) with modulus less than r are such that

Graph

moreover, CR can be applied with no breakdown and for any matrix norm &verbar; ⋅ &verbar;, the approximation error &verbar;X<sups>(n)</sups> − G&verbar; is such that

Graph

for a positive constant γ. We refer the reader to Bini[4] and to Gail[6] for more details in this regard.

A suitable technique for accelerating the convergence of CR has been introduced in He[10] and generalized in Bini[2]. It is based on replacing the original equation (1) with a slightly different one

Graph

where the new coefficients Aˆ<subs>i</subs> for i = − 1, 0,..., are generated according to the formulas

Graph

With this choice, the matrix power series is still analytic for |z| < r and the zeros η<subs>i</subs>, i = 1, 2,..., of with modulus less than r, are such that η<subs>1</subs> = 0, η<subs>i</subs> = ξ<subs>i−1</subs>, i = 2,..., m, η<subs>i</subs> = ξ<subs>i</subs> for i > m and therefore satisfy the relations

Graph

Moreover, the minimal non-negative solution G of the original equation (1) is obtained from the spectral minimal solution Gˆ of the new equation (4) by means of the simple expression

Graph

where a spectral minimal solution is a solution with minimal spectral radius.

In principle, CR can be applied to (4) to approximate Gˆ so that G = Gˆ + E but neither convergence nor applicability is guaranteed. In He[10] it has been proved that for QBD processes, where A<subs>i</subs> = 0 for i > 1, CR is convergent provided that no breakdown occurs. More precisely, the sequence {Xˆ<sups>(n)</sups>}<subs>n≥0</subs> is generated such that

Graph

for a positive constant hatγ. The acceleration that is obtained with this technique becomes particularly effective when |ξ<subs>m+1</subs>| is close to 1, this case occurs for problems close to being null recurrent where ξ<subs>m</subs> = ξ<subs>m+1</subs> = 1. Moreover, in many practical situations the acceleration is dramatic.

Observe that the zeros of det(zI − Aˆ(z)) are the same zeros of det(zI − A(z)) except for ξ<subs>m</subs> = 1, which is shifted to η<subs>1</subs> = 0. For this reason we refer to CR applied to solve (4) and complemented by (6) as Shifted Cyclic Reduction or SCR for short and to the new equation (4) as the shifted equation.

In this paper we analyze the applicability and the convergence of SCR under Assumptions 1 and 2. In particular we show that SCR always converges if no breakdown occurs and that SCR is asymptotically applicable in a sense that will be specified later on. Moreover, we prove the quadratic convergence and provide an explicit expression of the asymptotic behavior of the approximation error in terms of the zeros of det(zI − A(z)). Finally, we consider the case of QBD processes and show that breakdown never occurs. The latter property is particularly important since, as it is shown in Ramaswami[13] and Bini[4], any M/G/1 or G/M/1 problem can be simply transformed into a QBD.

The proof of these results relies on the functional interpretation of SCR and on a theorem on canonical factorizations.

The paper is organized in the following way: in section 2 we recall some theoretical results concerning matrix equations and canonical factorizations, and prove a result which, under suitable assumptions, relates the existence of left and right canonical factorization of a matrix function F(z) with the nonsingularity of the constant term H<subs>0</subs> of the Laurent series H(z) = z<sups>i</sups>H<subs>i</subs> = F(z)<sups>−1</sups>. This result is fundamental to prove convergence and applicability of SCR under very mild conditions. In Section 3 we recall the CR and SCR iterations, report their main properties and provide a functional description of both iterations. In section 4, we analyze convergence properties of SCR. We prove that the function F(z) = I − z<sups>−1</sups>A(z) has a right and a left canonical factorization so that the block H<subs>0</subs> in the Laurent expansion of H(z) = F(z)<sups>−1</sups> is nonsingular. The nonsingularity of H<subs>0</subs> is the key property that allows us to prove the convergence of SCR for M/G/1-type Markov chains and to evaluate its convergence rate. Asymptotic applicability of SCR, which is a direct consequence of the nonsigularity of H<subs>0</subs>, is used in section 5 to prove that for QBD problems breakdown of SCR never occurs. Section 6 reports the results of some numerical experiments that show the acceleration and the effectiveness of SCR with respect to CR.

2. THEORETICAL RESULTS

Throughout this section, F<subs>i</subs> for i = − 1, 0, 1,..., are real m×m matrices even though part of the results reported here still hold if F<subs>i</subs> are bounded operators. Moreover, given a matrix A = (a<subs>i, j</subs>) we denote |A| = (|a<subs>i, j</subs>|).

Definition 1

The Wiener algebra is the set &#55349;&#56498; of complex matrix valued functions F(z) = z<sups>i</sups>F<subs>i</subs> such that |F<subs>i</subs>| is finite.

Observe that under our assumptions, the function A(z) of (2) is in the Wiener algebra.

Definition 2

Let F(z) = z<sups>i</sups>F<subs>i</subs> ∈ &#55349;&#56498;. A factorization of the kind

Graph

is called left canonical factorization of F(z) if in the representation above U(z) = z<sups>i</sups>U<subs>i</subs> ∈ &#55349;&#56498;, L(z) = z<sups>−i</sups>L<subs>−i</subs> ∈ &#55349;&#56498;, and if U(z), L(z<sups>−1</sups>) are invertible for |z| ≤ 1, where U<subs>i</subs> and L<subs>−i</subs> are matrices of the same size as F<subs>i</subs>.

A factorization of the kind

Graph

is called right canonical factorization of F(z) if U(z) and L(z) satisfy the same properties as above.

Observe that the invertibility of U(z) and L(z<sups>−1</sups>) for |z| ≤ 1 implies that U<subs>0</subs> and L<subs>0</subs> are nonsingular.

We say that F(z) admits a weak left/right canonical factorization if in the representation above U(z) and L(z<sups>−1</sups>) are invertible for |z| < 1 and either U(z) or L(z) is singular for some z of modulus 1.

Canonical factorizations are closely related to matrix equations of the kind (1). In this regard we mention the following results (Bini[4]).

Theorem 2.1

Let F(z) = z<sups>i</sups>F<subs>i</subs> ∈ &#55349;&#56498;. If there exists a left canonical factorization

Graph

then L(z) = L<subs>0</subs> + z<sups>−1</sups>L<subs>−1</subs>, and G = − L<subs>−1</subs> is the unique solution of the matrix equation

Graph

such that ρ(G) < 1. In particular, G is the spectral minimal solution. Conversely, if there exists a solution G of the matrix equation (7), such that ρ(G) < 1, and if hboxdet(zF(z)) has exactly m roots in the open unit disk and is nonsingular for |z| = 1, then F(z) has a left canonical factorization

Graph

where U<subs>i</subs> = F<subs>j</subs>G<sups>j−i</sups>, i ≥ 0.

An analogous result relates the right canonical factorization and the minimal solution of the equation X<sups>i+1</sups>F<subs>i</subs> = 0 (see Bini[4]).

Weak canonical factorizations apply to the cases where F(z) is not invertible for some z of modulus 1 like in the Markov chains problems for the function F(z) = I − z<sups>−1</sups>A(z), which is singular for z = 1.

Theorem 2.2

Let F(z) = z<sups>i</sups>F<subs>i</subs> ∈ &#55349;&#56498;. If there exists a weak left canonical factorization

Graph

such that G = − L<subs>−1</subs> is power bounded, i.e., G<sups>i</sups> is uniformly bounded from above for any i > 0, then G is a spectral minimal solution of the equation (7) such that ρ(G) ≤ 1. Conversely, if F'(z) ∈ &#55349;&#56498;, if there exists a power bounded solution G of the matrix equation (7) such that ρ(G) = 1, and if all the zeros of hboxdet(zF(z)) in the open unit disk are eigenvalues of G then there exists a weak left canonical factorization of the form (8).

Now we are ready to prove the fundamental result that enables us to demonstrate convergence and applicability properties of SCR.

Theorem 2.3

Let F(z) = z<sups>i</sups>F<subs>i</subs> ∈ &#55349;&#56498; be invertible for |z| = 1 and denote H(z) = F(z)<sups>−1</sups>, H(z) = z<sups>i</sups>H<subs>i</subs>. If there exist a left and a right canonical factorization of F(z) then H<subs>0</subs> is invertible.

Proof

Left canonical factorization of F(z), if exists, has the form

Graph

with the matrix G having spectral radius strictly less than 1. Taking the inverses in both sides,

Graph

From here it follows, in particular, that

Graph

and

Graph

Suppose now that the coefficient H<subs>0</subs> is singular. Then there exists a non-zero vector <bold>x</bold> for which H<subs>0</subs><bold>x</bold> = 0. From (9) it follows that also H<subs>−i</subs><bold>x</bold> = 0 for all i = 1, 2,..., so that the vector function defined by

Graph

is analytic in the unit disk. On the other hand, the vector function

Graph

is analytic outside the unit disk, vanishes at infinity, and together with ϕ<sups>+</sups>(z) satisfies the boundary condition

Graph

Thus, the homogeneous Riemann-Hilbert problem (10) with the matrix coefficient H(z) has non-trivial solutions. This means (see, e.g.,[7]) that H(z) does not admit a left canonical factorization. Equivalently, its inverse F(z) does not admit a right canonical factorization. The contradiction obtained shows that H<subs>0</subs> is invertible.

Observe for completeness that condition on the Laurent series expansion of F(z) in Theorem 2.3, i.e., that the coefficients of z<sups>i</sups>, i < − 1 are zero, is essential. Namely, there even exist scalar functions F(z) in &#55349;&#56498; admitting a canonical factorization – both left and right, for which F(z)<sups>−1</sups> has a vanishing constant term. As a concrete example, let H(z) = − 2z<sups>−4</sups> + 2z<sups>−1</sups> − 3z. Since | − 2| + |2| + | − 3| = 7 < + ∞, this function belongs to &#55349;&#56498;, moreover, since H(z) is non-vanishing on the unit circle, then H(z)<sups>−1</sups> ∈ &#55349;&#56498;. Now, let us set F(z) = H(z)<sups>−1</sups> = z<sups>i</sups>F<subs>i</subs> and observe that, since the winding number of H(z) is equal to zero ([7]), the function H(z) = F(z)<sups>−1</sups> (and therefore its inverse F(z)) admits a canonical factorization. Nevertheless, the constant term H<subs>0</subs> of F(z)<sups>−1</sups> = H(z) is indeed zero.

3. SHIFTED CYCLIC REDUCTION

Define the matrix power series

Graph

where A(z) = z<sups>i+1</sups>A<subs>i</subs> satisfies Assumptions 1 and 2, and set F(z) = z<sups>−1</sups>ϕ(z) = z<sups>i</sups>F<subs>i</subs>. Then there exists the minimal solution G of (1) such that G is stochastic and therefore, G is power bounded and ρ(G) = 1 (Neuts[12]). Since F(z) is analytic for 0 < |z| < r then also F'(z) is analytic for 0 < |z| < r, moreover, since F<subs>i</subs> ≤ 0 for i > 0, then both F(z) and F'(z) belong to &#55349;&#56498;. Therefore, we may apply Theorem 2.2 so that there exists the weak left canonical factorization

Graph

Since the zeros in the closed unit disk of detϕ(z) are the eigenvalues of G, then U(z) is nonsingular for |z| ≤ 1. Moreover,

Graph

Define the matrix power series hatϕ(z) = zI − z<sups>i+1</sups>Aˆ<subs>i</subs> where , defined in (5), are associated with the "shifted" matrix equation (4) and set Fˆ(z) = z<sups>−1</sups>hatϕ(z). It is a simple matter to show that hatϕ(z) = ϕ(z)(I − z<sups>−1</sups>E)<sups>−1</sups> and

Graph

Moreover, since (I − z<sups>−1</sups>G)(I − z<sups>−1</sups>E)<sups>−1</sups> = I − z<sups>−1</sups>Gˆ, where Gˆ = G − E, and ρ(Gˆ) < 1 (Bini et al.[2]), we deduce that

Graph

with Uˆ(z) = U(z), is a left canonical factorization of Fˆ(z). Moreover, it holds

Graph

In Bini et al.[4] it is shown that hatϕ(z) is analytic for |z| < r and hatϕ(z) ∈ &#55349;&#56498; so that the hypotheses of Theorem 2.1 hold for Fˆ(z). Therefore, Gˆ is the spectral minimal solution of (4).

Throughout the rest of the paper, given a function f(z), we define its even and odd parts as

Graph

Cyclic reduction applied to the matrix power series ϕ(z) consists in generating the sequences {ϕ<sups>(n)</sups>(z)}<subs>n≥0</subs> and {θ<sups>(n)</sups>(z)}<subs>n≥0</subs> of matrix power series and the sequence {X<sups>(n)</sups>}<subs>n≥0</subs> of matrices by means of

Graph

for n = 0, 1,..., starting with ϕ<sups>(0)</sups>(z) = ϕ(z) and θ<sups>(0)</sups>(z) = z<sups>i</sups>A<subs>i</subs>, where we set θ<sups>(n)</sups>(z) = z<sups>i</sups> (see Bini and Meini[5], Bini et al.[4]). Moreover we have the following fundamental relation

Graph

which, after post-multiplying both sides by G, yields

Graph

The latter equation provides the following expression for the error G − X<sups>(n)</sups>:

Graph

The sequence X<sups>(n)</sups> can be computed with a low computational cost, we refer the reader to Bini and Meini[5] and to Bini et al.[4] concerning the computational details of CR. Observe that the matrix power series ϕ<sups>(n)</sups>(z) and θ<sups>(n)</sups>(z) generated by CR in (13) are well defined if [ϕ<sups>(n)</sups>(z)]<subs>odd</subs> is invertible for |z| ≤ 1 so that its inverse is still a matrix power series.

In Bini et al.[4] it is shown that under Assumptions 1 and 2, the matrix power series ϕ<sups>(n)</sups>(z), θ<sups>(n)</sups>(z) and [ϕ<sups>(n)</sups>(z)]<subs>odd</subs> belong to &#55349;&#56498; and [ϕ<sups>(n)</sups>(z)]<subs>odd</subs> is invertible for |z| ≤ 1 so that CR can be applied with no breakdown. Moreover, I −  is invertible and the sequence X<sups>(n)</sups> quadratically converges to G. More precisely, for any norm &verbar; ⋅ &verbar; and for any σ such that || < σ < 1, there exists a constant γ > 0 such that &verbar;G − X<sups>(n)</sups>&verbar; ≤ γσ<sups>2n</sups>.

A better functional interpretation of CR can be given by introducing the matrix valued functions ψ<sups>(n)</sups>(z) = (z<sups>−1</sups>ϕ<sups>(n)</sups>(z))<sups>−1</sups>, for n ≥ 0. In Bini et al.[4] it is shown that

Graph

so that, by setting ψ<sups>(0)</sups>(z) = H(z) = z<sups>i</sups>H<subs>i</subs>, one has

Graph

Moreover, by using (13), we may show that (see Bini et al.[4]) θ<sups>(n+1)</sups>(z<sups>2</sups>)ψ<sups>(n+1)</sups>(z<sups>2</sups>) = [θ<sups>(n)</sups>(z)ψ<sups>(n)</sups>(z)]<subs>even</subs> so that we obtain

Graph

where

Graph

Similarly, the shifted cyclic reduction (SCR) algorithm, that is, CR applied to hatϕ(z), consists in generating the sequences {hatϕ<sups>(n)</sups>(z)}<subs>n≥0</subs> and {hatθ<sups>(n)</sups>(z)}<subs>n≥0</subs> of matrix power series and the sequence {Xˆ<sups>(n)</sups>}<subs>n≥0</subs> of matrices by means of

Graph

for n = 0, 1,..., starting with hatϕ<sups>(0)</sups>(z) = hatϕ(z) and . Observe that the expression of Xˆ<sups>(n)</sups> is obtained from that of X<sups>(n)</sups>, by replacing with hat.

Moreover we have the following fundamental relation

Graph

which, after multiplying both sides on the right by Gˆ, yields

Graph

Similarly, for the matrix functions hatψ<sups>(n)</sups>(z) = (z<sups>−1</sups>hatϕ<sups>(n)</sups>(z))<sups>−1</sups> it holds

Graph

so that, by setting one has

Graph

Moreover, it holds

Graph

where

Graph

Remark

A nice consequence of (17) is that the sequence {hatψ<sups>(n)</sups>(z)} can be defined in terms of hatψ<sups>(0)</sups>(z) for any n ≥ 0 even though for some n the matrix [hatϕ<sups>(n)</sups>(z)]<subs>odd</subs> is singular at some complex value z of modulus less than or equal to 1. That is, hatψ<sups>(n)</sups>(z) exist for any n even in the case where the SCR step cannot be applied and hatϕ<sups>(n)</sups>(z) cannot be defined. Moreover, hatϕ<sups>(n)</sups>(z) = zhatψ<sups>(n)</sups>(z)<sups>−1</sups> is defined for any z for which hatψ<sups>(n)</sups>(z) is invertible. This means that, even though breakdown occurs in SCR, the functions hatϕ<sups>(n)</sups>(z) still exist.

4. CONVERGENCE

In this section we provide convergence properties for the iteration SCR of (15). The key tool for proving convergence is given by Theorem 2.3. Therefore, the first step of our analysis is to prove that the assumptions of Theorem 2.3 are satisfied by the function Fˆ(z) = z<sups>−1</sups>hatϕ(z). We have the followings

Theorem 4.1

Assume that Assumptions 1, 2 are satisfied and let <bold>u</bold> be a positive vector such that <bold>u</bold><sups>T</sups><bold>e</bold> = 1. Let Fˆ(z) = I − z<sups>−1</sups>Aˆ(z), where is the matrix power series whose coefficients Aˆ<subs>i</subs>, i ≥ − 1 are defined by (5). Then there exists the spectral minimal solution Rˆ of the equation

Graph

and the function Fˆ(z) has the right canonical factorization

Graph

where Lˆ(z) = z<sups>i</sups>Lˆ<subs>i</subs>.

Proof

Since the M/G/1 Markov chain is positive recurrent, then there exists the minimal non-negative solution R of the equation X = X<sups>i+1</sups>A<subs>i</subs>, and the matrix R has spectral radius ρ(R) = 1 (Neuts[12]). Moreover ρ(R) is a simple eigenvalue and is the unique eigenvalue of modulus 1 of R. In light of Theorem 2.2 there exists a weak right canonical factorization

Graph

where S(z) = z<sups>i</sups>S<subs>i</subs>, S<subs>0</subs> is nonsingular, and S(z) is analytic for |z| < r. Observe that since 0 = F(1)<bold>e</bold> = (I − R)S(1)<bold>e</bold>, then <bold>v</bold> = S(1)<bold>e</bold> is a right eigenvector of R corresponding to the eigenvalue ρ(R) = 1. Moreover, <bold>v</bold> is non-negative since R is non-negative. Let us denote by <bold>w</bold> a non-negative vector such that <bold>w</bold><sups>T</sups><bold>v</bold> = 1. Rewrite (22) as

Graph

A simple calculation shows that (I − z<sups>−1</sups>R)(I − z<sups>−1</sups><bold>v</bold><bold>w</bold><sups>T</sups>)<sups>−1</sups> = I − z<sups>−1</sups>Rˆ for Rˆ = R − <bold>v</bold><bold>w</bold><sups>T</sups>. Therefore, since Fˆ(z) = F(z)(I − z<sups>−1</sups>E)<sups>−1</sups>, for E = <bold>e</bold><bold>u</bold><sups>T</sups>, we have

Graph

Now, taking determinants of zI − Rˆ = (zI − R)(I − z<sups>−1</sups><bold>v</bold><bold>w</bold><sups>T</sups>)<sups>−1</sups> yields $ . This shows that 0 is an eigenvalue of Rˆ together with all the eigenvalues of R different from 1 and that, since 1 is a simple eigenvalue of R, it cannot be eigenvalue of Rˆ. This implies that ρ(Rˆ) < 1 and that the matrix function I − z<sups>−1</sups>Rˆ is invertible for |z| ≥ 1. Now we prove that the rightmost factor Lˆ(z) = (I − z<sups>−1</sups><bold>v</bold><bold>w</bold><sups>T</sups>)S(z)(I − z<sups>−1</sups>E)<sups>−1</sups> can be extended for z = 1 by continuity, is a matrix power series in the Wiener algebra &#55349;&#56498; and is invertible for |z| ≤ 1. We first analyze the product B(z) = S(z)(I − z<sups>−1</sups>E)<sups>−1</sups> = z<sups>i</sups>B<subs>i</subs>. It holds

Graph

in fact, if i < 0 then B<subs>i</subs> = S<subs>j</subs>E = S(1)E = S(1)<bold>e</bold><bold>u</bold><sups>T</sups> = <bold>v</bold><bold>u</bold><sups>T</sups>. Now we are ready to show that for a suitable <bold>w</bold> the product

Graph

is a matrix power series in the Wiener algebra &#55349;&#56498;. Indeed, since Lˆ(z) = z<sups>i</sups>B<subs>i</subs> − z<sups>i</sups>(<bold>v</bold><bold>w</bold><sups>T</sups>B<subs>i+1</subs>), we obtain

Graph

and from (23) we deduce that Lˆ<subs>−1</subs> = B<subs>−1</subs> − <bold>v</bold><bold>w</bold><sups>T</sups>B<subs>0</subs> = <bold>v</bold>(<bold>u</bold><sups>T</sups> − <bold>w</bold><sups>T</sups>B<subs>0</subs>). Therefore, since  ≥ 0 (Neuts[12]), then <bold>u</bold><sups>T</sups><bold>v</bold> > 0, and choosing <bold>w</bold><sups>T</sups> = (1/<bold>u</bold><sups>T</sups><bold>v</bold>)<bold>u</bold><sups>T</sups> yields <bold>w</bold><sups>T</sups><bold>v</bold> = 1 and

Graph

where we used the property S<subs>j</subs><bold>e</bold> = <bold>v</bold>. This implies that hatL<subs>−1</subs> = 0. From (23) and (24) it follows that Lˆ<subs>i</subs> = 0 for i < − 1, that is Lˆ(z) is a matrix power series. Moreover, |Lˆ<subs>i</subs>| ≤ |B<subs>i</subs>| + |<bold>v</bold><bold>w</bold><sups>T</sups>B<subs>i+1</subs>|. Therefore, in order to prove that Lˆ(z) ∈ &#55349;&#56498; it is sufficient to show that |B<subs>i</subs>| is finite. To prove this property it is sufficient to observe that B(z) = S(z)(I − z<sups>−1</sups>E)<sups>−1</sups> is analytic for 1 < |z| < r since F(z) is analytic for 1 < |z| < r. Therefore by the Cauchy integral theorem, the coefficients B<subs>i</subs> for i > 0 have entries bounded in modulus from above by γσ<sups>−i</sups> where 1 < σ < r and γ is a constant so that |B<subs>i</subs>| is convergent. Finally, in order to complete the proof we have to show that Lˆ(z) is invertible for |z| ≤ 1. This property holds since det(zFˆ(z)) has exactly m zeros in the open unit disk, and no zeros of modulus 1. Therefore, the only values of z of modulus less than 1 which make Fˆ(z) noninvertible are the zeros of det(I − z<sups>−1</sups>Rˆ) that is the eigenvalues of Rˆ, and the factorization implies the nonsingularity of Lˆ(z) for |z| ≤ 1.

Observe that Fˆ(z) has the left canonical factorization (12) and the right canonical factorization (21) so that the hypotheses of Theorem 3 are satisfied. Therefore, is such that detHˆ<subs>0</subs> ≠ 0.

Observe also that since Fˆ(z) is analytic and invertible in the annulus

Graph

then also Hˆ(z) is analytic in &#55349;&#56476;.

Theorem 4.2

Let A(z) be the generating function associated with an M/G/1-type Markov chain which satisfies Assumptions 1 and 2. Let ξ<subs>i</subs>, i ≥ 1 be the roots of ϕ(z) = zI − A(z) such that (3) holds. Let <bold>u</bold> > 0 be such that <bold>u</bold><sups>T</sups><bold>e</bold> = 1, set E = <bold>e</bold><bold>u</bold><sups>T</sups>, hatϕ(z) = ϕ(z)(I − z<sups>−1</sups>E)<sups>−1</sups>. Then for any ϵ > 0 such that |ξ<subs>m−1</subs>| + ϵ < 1 < |ξ<subs>m+1</subs>| − ϵ, there exist a real 0 < β < 1 and a positive integer n<subs>0</subs> such that for any n ≥ n<subs>0</subs> the function hatψ<sups>(n)</sups>(z) defined in (17) is analytic and invertible in the annulus calA<subs>n</subs>(β, ϵ) = {z ∈ ℂ:β<sups>−1</sups>(|ξ<subs>m−1</subs>| + ϵ)<sups>2n</sups> < |z| < β(|ξ<subs>m+1</subs>| − ϵ)<sups>2n</sups>}, and has a uniformly bounded inverse for z ∈ calA<subs>n</subs>(β, ϵ). Moreover, for any n ≥ n<subs>0</subs>, the function

Graph

is analytic and invertible in &#55349;&#56476;<subs>n</subs>(β, ϵ), and for any operator norm &verbar; ⋅ &verbar; there exist positive constants c<subs>i</subs>, i ∈ ℤ, such that

Graph

for any n ≥ n<subs>0</subs>.

Proof

Let ϵ > 0 be such that |ξ<subs>m−1</subs>| + ϵ < 1 < |ξ<subs>m+1</subs>| − ϵ. By using an inductive argument on n, from (17) it follows that hatψ<sups>(n)</sups>(z) is analytic in the annulus {z ∈ ℂ:(|ξ<subs>m−1</subs>| + ϵ)<sups>2n</sups> < |z| < (|ξ<subs>m+1</subs>| − ϵ)<sups>2n</sups>} which contains &#55349;&#56476;<subs>n</subs>(β, ϵ) for any 0 < β < 1.

Now we prove that there exist n<subs>0</subs> > 0 and 0 < β < 1 such that, for any n ≥ n<subs>0</subs>, the matrix hatψ<sups>(n)</sups>(z) is nonsingular for z ∈ &#55349;&#56476;<subs>n</subs>(β, ϵ) and has a uniformly bounded inverse.

Recall that z<sups>−1</sups>hatϕ(z) has a left canonical factorization and, by Theorem 4.1, a right canonical factorization. Therefore we may apply Theorem 2.3 and deduce that the matrix $ is nonsingular. Thus, we may write where

Graph

In order to prove the nonsingularity of hatψ<sups>(n)</sups>(z) and the boundedness of its inverse for z ∈ &#55349;&#56476;<subs>n</subs>(β, ϵ) it is sufficient to prove that there exists an operator norm such that &verbar;W<sups>(n)</sups>(z)&verbar; ≤ 1/2 for z ∈ &#55349;&#56476;<subs>n</subs>(β, ϵ). In fact, the latter bound implies that ρ(W<sups>(n)</sups>(z)) ≤ 1/2 for z ∈ calA<subs>n</subs>(β, ϵ) so that det(I + W<sups>(n)</sups>(z)) ≠ 0. Moreover, since and since &verbar;(I + W<sups>(n)</sups>(z))<sups>−1</sups>&verbar; ≤ 1/(1 − &verbar;W<sups>(n)</sups>(z)&verbar;), then the condition &verbar;W<sups>(n)</sups>(z)&verbar; ≤ 1/2 for z ∈ calA<subs>n</subs>(β, ϵ), implies &verbar;(I + W<sups>(n)</sups>(z))<sups>−1</sups>&verbar; ≤ 2 for z ∈ &#55349;&#56476;<subs>n</subs>(β, ϵ), that is, the uniform boundedness of hatψ<sups>(n)</sups>(z)<sups>−1</sups>.

Let us prove that there exists an operator norm &verbar; ⋅ &verbar; and there exists 0 < β < 1 and n<subs>0</subs> > 0 such that for any n ≥ n<subs>0</subs>, &verbar;W<sups>(n)</sups>(z)&verbar; < 1/2 for z ∈ &#55349;&#56476;<subs>n</subs>(β, ϵ). Since hatψ(z) is analytic in the annulus &#55349;&#56476; of (25), then by the Cauchy integral theorem for any operator norm &verbar; ⋅ &verbar; there exists a constant c such that

Graph

Let n<subs>0</subs> be such that for any n ≥ n<subs>0</subs>. Then from (28) and (27) we obtain that

Graph

where the latter inequality is valid for any n ≥ n<subs>0</subs> provided that , with γ = 4c&verbar;&verbar;. Therefore, it is sufficient to prove that there exists a constant 0 < β < 1 such that

Graph

and

Graph

Since (|ξ<subs>m−1</subs>| + ϵ)/(|ξ<subs>m+1</subs>| − ϵ) < 1, equation (30) is satisfied for . Moreover, (29) is satisfied for any β ≤ 1/(4γ). Choosing β = 1/(4γ) and n<subs>0</subs> such that then (30) and (29) are satisfied for any n ≥ n<subs>0</subs>.

Concerning the last part of the theorem, recall that hatϕ<sups>(n)</sups>(z) = zhatψ<sups>(n)</sups>(z)<sups>−1</sups> and hatψ<sups>(n)</sups>(z)<sups>−1</sups> is analytic and bounded in calA<subs>n</subs>(β, ϵ) for any n ≥ n<subs>0</subs>. Therefore, from the Cauchy integral theorem applied to hatψ<sups>(n)</sups>(z)<sups>−1</sups> one deduces that there exist positive constants γ and δ < 1 such that &verbar;hat&verbar; ≤ γδ<sups>−i+1</sups>(|ξ<subs>m+1</subs>| − ϵ)<sups>−(i−1)2n</sups> for i ≥ 1, &verbar;hat&verbar; ≤ γδ<sups>i−1</sups>(|ξ<subs>m−1</subs>| + ϵ)<sups>−(i−1)2n</sups> for i ≤ 0. Therefore, choosing c<subs>i</subs> = γδ<sups>−i+1</sups> for i ≥ 1 and c<subs>i</subs> = γδ<sups>i−1</sups> for i ≤ 0 yields (26).

The above theorem does not exclude that break-down may occur in SCR. It rather shows that the functions hatψ<sups>(n)</sups>(z), which are defined by means of (17) for any n ≥ 0, are invertible for n ≥ n<subs>0</subs> sufficiently large so that also hatϕ<sups>(n)</sups>(z) = zhatψ<sups>(n)</sups>(z)<sups>−1</sups> are defined and invertible for n ≥ n<subs>0</subs>.

This shows the existence of hatϕ<sups>(n)</sups>(z) for n ≥ n<subs>0</subs> independently of the possible break-down of SCR due to the singularity of some [hatϕ<sups>(n)</sups>(z)]<subs>odd</subs> for certain n < n<subs>0</subs>. The assumption that SCR can be applied with no break-down is needed in order to provide an algorithmic construction of the sequence hatϕ<sups>(n)</sups>(z). In the following we assume this condition.

The asymptotic invertibility of the functions [hatϕ<sups>(n)</sups>(z)]<subs>odd</subs> is summarized in the following:

Theorem 4.3

SCR is asymptotically applicable, that is, the function [hatϕ<sups>(n)</sups>(z)]<subs>odd</subs> can be singular, and therefore hatϕ<sups>(n+1)</sups>(z) is undefined, only for a finite number of values of n.

From Theorem 4.2 we deduce the following result concerning the convergence of the matrix power series hatθ<sups>(n)</sups>(z) = z<sups>i</sups>hat generated by SCR (15).

Theorem 4.4

Assume that the hypotheses of Theorem 4.2 are satisfied and that SCR can be carried out with no break-down. If det(I − hat) ≠ 0 and if &verbar;(I − hat)<sups>−1</sups>&verbar; is uniformly bounded from above, then łim<subs>n</subs>Xˆ<sups>(n)</sups> = G, where Xˆ<sups>(n)</sups> = (I − hat)<sups>−1</sups>A<subs>−1</subs>. Moreover,

Graph

and for any matrix norm and for any ϵ > 0 such that |ξ<subs>m−1</subs>| + ϵ < 1 < |ξ<subs>m+1</subs>| − ϵ, there exist γ > 0 and σ<subs>i</subs> > 0, i ≥ 1, such that

Graph

and

Graph

Proof

Let us prove (32). From (19), denoting V(z) = hatθ(z)hatψ(z), one has hatθ<sups>(n)</sups>(z)hatψ<sups>(n)</sups>(z) = V<sups>(n)</sups>(z), where V<sups>(n)</sups>(z) = z<sups>i</sups>V<subs>i2n</subs>. Therefore,

Graph

Comparing the coefficients in the above expression yields

Graph

Taking norms of both sides yields

Graph

Because V(z) is analytic for |ξ<subs>m−1</subs>| < |z| < |ξ<subs>m+1</subs>|, then for any ϵ > 0 such that |ξ<subs>m−1</subs>| + ϵ < 1 < |ξ<subs>m+1</subs>| − ϵ, there exists β > 0 such that &verbar;V<subs>k⋅2n</subs>&verbar; ≤ β(|ξ<subs>m+1</subs>| − ϵ)<sups>−k⋅2n</sups> for any n, k > 0 and &verbar;V<subs>k⋅2n</subs>&verbar; ≤ β(|ξ<subs>m−1</subs>| + ϵ)<sups>−k⋅2n</sups> for any k < 0, n > 0. Moreover, from (26) we have &verbar;hat&verbar; ≤ c<subs>i</subs>(|ξ<subs>m+1</subs>| − ϵ)<sups>−(i−1)2n</sups>, i ≥ 1. Therefore we find that there exists σ<subs>k</subs> > 0 such that

Graph

for any n, k > 0. Equations (31)s can be proved by induction. Concerning the convergence of Xˆ<sups>(n)</sups> to G, since

Graph

from equation (31) one has (I − hat)G = A<subs>−1</subs> + hatGˆ<sups>i⋅2n</sups>. Therefore we obtain

Graph

Since ρ(Gˆ) = |ξ<subs>m−1</subs>|, for any ϵ > 0 there exists a matrix norm &verbar; ⋅ &verbar; such that . Therefore, for the equivalence of matrix norms, for any matrix norm there exists a constant β > 0 such that &verbar;Gˆ<sups>i2n</sups>&verbar; ≤ β(|ξ<subs>m−1</subs>| + ϵ)<sups>i2n</sups>, so that

Graph

Because &verbar;(I − hat)<sups>−1</sups>&verbar; is bounded from above by a constant, from (32) we have for a suitable γ.

5. APPLICABILITY

If F(z) = − A<subs>−1</subs>z<sups>−1</sups> + I − A<subs>0</subs> − A<subs>1</subs>z and , then we may show that CR applied to hatϕ(z) can be carried out with no break-down provided that Assumption 1 is satisfied. We first recall that in principle CR can encounter a break-down if the matrix [ϕ<sups>(n)</sups>(z)]<subs>odd</subs> = (I − ) is singular. It is shown in Bini et al.[3] that det(I − ) is equal to zero if and only if the block tridiagonal matrix T<subs>2n+1−1</subs> is singular, where we denote by T<subs>k</subs> the mk×mk leading principal submatrix of the semi-infinite block tridiagonal block Toeplitz matrix

Graph

Moreover, in Bini et al.[4] it is shown that under Assumption 1, CR can be carried out with no breakdown so that detT<subs>2n+1−1</subs> ≠ 0 for any n ≥ 0.

In order to show that SCR can be carried out with no break-down it is sufficient to prove that detTˆ<subs>2n+1−1</subs> ≠ 0 for any n ≥ 0 where Tˆ<subs>k</subs> is the mk×mk leading principal submatrix of the semi-infinite block tridiagonal block Toeplitz matrix

Graph

In order to prove this property we need some preliminary results.

First of all we prove that, if CR can be applied to ϕ(z), then detT<subs>k</subs> ≠ 0 for any k > 1. Assume by absurd that there exists some k<subs>0</subs> such that detT<subs>k0</subs> = 0 then ρ(I − T<subs>k0</subs>) = 1. Let k be any integer such that k > k<subs>0</subs>. Denoting with U<subs>k</subs> the mk×mk matrix having I − T<subs>k0</subs> as leading principal submatrix and null entries elsewhere, it holds 0 ≤ U<subs>k</subs> ≤ I − T<subs>k</subs>. From the Perron–Frobenius theorem we get 1 = ρ(U<subs>k</subs>) ≤ ρ(I − T<subs>k</subs>) ≤ 1. Thus, we obtain that ρ(I − T<subs>k</subs>) = 1 and applying once again the Perron–Frobenius theorem we find that there exists an eigenvalue of I − T<subs>k</subs> equal to 1, therefore the matrix T<subs>k</subs> is singular for any k > k<subs>0</subs>. This would imply that CR cannot be carried out in contradiction with our assumptions.

Theorem 5.1

Let T<subs>k</subs> be nonsingular, set H<subs>k</subs> =  and denote with P<subs>k</subs> the m×m block of H<subs>k</subs> in the lower rightmost corner. Then the matrix Tˆ<subs>k</subs> is singular if and only if <bold>u</bold><sups>T</sups>P<subs>k</subs>A<subs>1</subs><bold>e</bold> = 1.

Proof

The equation hatϕ(z) = ϕ(z)(I − z<sups>−1</sups>E)<sups>−1</sups> can be equivalently rewritten in the matrix form as

Graph

whence we deduce that

Graph

This implies that Tˆ<subs>k</subs> is singular if and only if the matrix in the right-hand side of the above expression is singular, that is, if and only if

Graph

The latter determinant is given by 1 − <bold>u</bold><sups>T</sups>P<subs>k</subs>A<subs>1</subs><bold>e</bold>.

Theorem 5.2

Let <bold>u</bold> > 0 and detT<subs>k</subs> ≠ 0, then <bold>u</bold><sups>T</sups>P<subs>k</subs>A<subs>1</subs><bold>e</bold> = 1 if and only if P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold>.

Proof

It is clear that if P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold> then <bold>u</bold><sups>T</sups>P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>u</bold><sups>T</sups><bold>e</bold> = 1. In order to prove the reverse implication, observe that, since (A<subs>−1</subs> + A<subs>0</subs> + A<subs>1</subs>)<bold>e</bold> = <bold>e</bold>, then

Graph

Whence we obtain

Graph

Therefore from the last m rows in the above expression we get

Graph

where Q<subs>k</subs> is the m×m submatrix in the lower leftmost corner of H<subs>k</subs>. Since H<subs>k</subs> is the inverse of a nonsingular M-matrix, it is nonnegative (Berman and Plemmons[1], Varga[14]) so that Q<subs>k</subs>A<subs>−1</subs><bold>e</bold> ≥ 0 and 0 ≤ P<subs>k</subs>A<subs>1</subs><bold>e</bold> ≤ <bold>e</bold>. This implies that if P<subs>k</subs>A<subs>1</subs><bold>e</bold> ≠ <bold>e</bold> then <bold>u</bold><sups>T</sups>P<subs>k</subs>A<subs>1</subs><bold>e</bold> − <bold>e</bold> ≠ 0, which completes the proof.

Now we are ready to prove the main applicability result:

Theorem 5.3

Let A<subs>−1</subs>, A<subs>0</subs> and A<subs>1</subs> be m×m non-negative matrices. Let us assume that A<subs>−1</subs> + A<subs>0</subs> + A<subs>1</subs> is row stochastic and that CR can be applied to ϕ<subs>k</subs>(z). If <bold>u</bold> > 0 and if I − A<subs>0</subs> − A<subs>1</subs>E is nonsingular, then SCR can be applied with no break-down.

Proof

We prove that detTˆ<subs>k</subs> ≠ 0 for k ≥ 1. If k = 1 then Tˆ<subs>1</subs> = I − Aˆ<subs>0</subs> = I − A<subs>0</subs> − A<subs>1</subs>E. Therefore the nonsingularity of Tˆ<subs>1</subs> follows by assumption.

Let us consider the case k > 1. Since CR can be applied, then detT<subs>k</subs> ≠ 0 for any k. Apply Theorems 5.1 and 5.2 and deduce that detTˆ<subs>k</subs> = 0 if and only if P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold>. Assume by absurd that detTˆ<subs>k0</subs> = 0 for some k<subs>0</subs>, so that P<subs>k0</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold>. Then we prove that P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold> for any k ≥ k<subs>0</subs> so that, in view of Theorems 5.1 and 5.2, the matrix Tˆ<subs>k</subs> would be singular for any k ≥ k<subs>0</subs>. This property contradicts the asymptotic applicability guaranteed by Theorem 4.3. In order to prove that P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold>, consider the Schur complement S<subs>k</subs> of T<subs>k−1</subs> in T<subs>k</subs> and recall that since T<subs>k</subs> and T<subs>k−1</subs> are nonsingular then S<subs>k</subs> exists, is nonsingular and  = P<subs>k</subs>. Moreover, from the definition of the Schur complement and from the block tridiagonal structure of T<subs>k</subs> and T<subs>k+1</subs> it holds that S<subs>k+1</subs> = (I − A<subs>0</subs>) − A<subs>−1</subs>A<subs>1</subs> = (I − A<subs>0</subs>) − A<subs>−1</subs>P<subs>k</subs>A<subs>1</subs>. Therefore, if P<subs>k</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold> then S<subs>k+1</subs><bold>e</bold> = (I − A<subs>0</subs>)<bold>e</bold> − A<subs>−1</subs>A<subs>1</subs><bold>e</bold> = <bold>e</bold> − A<subs>0</subs><bold>e</bold> − A<subs>−1</subs><bold>e</bold> = A<subs>1</subs><bold>e</bold>, where the latter equation holds since A<subs>−1</subs> + A<subs>0</subs> + A<subs>1</subs> is stochastic. Thus we get <bold>e</bold> = P<subs>k+1</subs>A<subs>1</subs><bold>e</bold> which completes the proof.

Theorem 5.3 guarantees that no break-down occurs in SCR, but it does not guarantee from near-breakdown situations where the nonsingular matrix T<subs>k</subs> is close to be singular (ill-conditioned). The analysis of near-breakdown situations is problem of great interest which is under investigation.

6. NUMERICAL TESTS

In order to test the effectiveness of the shifted cyclic reduction technique, we implemented the formulas which provide the matrices Aˆ<subs>i</subs> given the matrices A<subs>i</subs>, for i = − 1, 0,..., N, as a Fortran 95 subroutine. Then we applied the program of Bini and Meini[5] to the original data A<subs>i</subs> and to the "shifted" data Aˆ<subs>i</subs>; this program implements the point-wise cyclic reduction and is available at www.dm.unipi.it/meini. In the shifted case we have modified the stop condition since the one used in the original program does not apply to this case. In fact the iteration is halted if &verbar;hat − hat&verbar;<subs>∞</subs> ≤ ϵ, for ϵ = 10<sups>−12</sups>. We used the fortran 95 compiler Lahey-Fujitsu L6.20a and we ran it on a laptop with a Pentium 4 cpu and with the linux operating system. In our experiments we reported the number of iterations it needed to arrive at convergence, the residual error err, the maximum number of interpolation points inter needed by the point-wise technique, and the cpu time cpu needed for the computation. We denoted with PWCR the algorithm of Bini and Meini[5] based on point-wise cyclic reduction and with SPWCR the algorithm obtained by applying point-wise cyclic reduction after shift.

We have considered the following problems.

Problem 1

The matrices defining this problem are A<subs>−1</subs> = (1 − 2α)C, A<subs>1</subs> = (α/m)E, A<subs>2</subs> = (α/m)E, and A<subs>i</subs> = 0 elsewhere, where 0 < α ≤ 2/5, E is the m×m matrix with all the entries equal to 1, and C is the m×m circulant matrix

Graph

The roots of a(z) = det(A<subs>−1</subs> − zI + z<sups>2</sups>A<subs>1</subs> + z<sups>3</sups>A<subs>2</subs>) are the roots of the poly\-nomials (1 − 2α) − z + αz<sups>2</sups> + αz<sups>3</sups> and (1 − 2α)ω<sups>i</sups> − z, for i = 1,..., m − 1, where , and <bold>i</bold><sups>2</sups> = − 1. If α < 2/5 the problem is positive recurrent, that is, 1 is a simple zero of a(z), whereas if α = 2/5 the associated Markov chain is null recurrent and 1 is a double zero of a(z). For this problem the data are such that the shifted cyclic reduction algorithm converges in just one step.

We have chosen α = 0.4, α = 0.2 and α = 0.1. In Table 1 we report the results concerning the experiments for this problem.

TABLE 1 Problem 1: number it of iterations, residual error err, number of interpolation points inter and cpu time in seconds for the algorithms PWCR and SPWCR

<table><thead valign="bottom"><tr><td /><td>PWCR</td><td>SPWCR</td></tr><tr><td>&#945;</td><td>it</td><td>err</td><td>inter</td><td>cpu</td><td>it</td><td>err</td><td>inter</td><td>cpu</td></tr></thead><tbody><tr><td>0.1</td><td>3</td><td>1.4e-08</td><td>8</td><td>0.11</td><td>1</td><td>1.8e-15</td><td>8</td><td>0.05</td></tr><tr><td>0.2</td><td>5</td><td>4.7e-11</td><td>16</td><td>0.20</td><td>1</td><td>3.2e-15</td><td>16</td><td>0.13</td></tr><tr><td>0.4</td><td>18</td><td>1.2e-11</td><td>16</td><td>0.67</td><td>1</td><td>2.4e-15</td><td>32</td><td>0.24</td></tr></tbody></table>

Problem 2

The matrices defining this problem are A<subs>−1</subs> = (1 − 2α)C, A<subs>1</subs> = αC, A<subs>2</subs> = (α/m)E, and A<subs>i</subs> = 0 elsewhere, where 0 < α ≤ 2/5. The roots of a(z) = det(A<subs>−1</subs> − zI + z<sups>2</sups>A<subs>1</subs> + z<sups>3</sups>A<subs>2</subs>) are the roots of the polynomials (1 − 2α) − z + αz<sups>2</sups> + αz<sups>3</sups> and (1 − 2α)ω<sups>i</sups> − z + αω<subs>i</subs>z<sups>2</sups>, for i = 1,..., m − 1.

If α < 2/5 the problem is positive recurrent, that is, 1 is a simple zero of a(z), whereas if α = 2/5 the associated Markov chain is null recurrent and 1 is a double zero of a(z).

In Table 2 we report the results concerning the experiments for this problem.

TABLE 2 Problem 2: number it of iterations, residual error err, number of interpolation points inter and cpu time in seconds for the algorithms PWCR and SPWCR

<table><thead valign="bottom"><tr><td /><td>PWCR</td><td>SPWCR</td></tr><tr><td>&#945;</td><td>it</td><td>err</td><td>inter</td><td>cpu</td><td>it</td><td>err</td><td>inter</td><td>cpu</td></tr></thead><tbody><tr><td>0.1</td><td>4</td><td>1.4e-08</td><td>8</td><td>0.12</td><td>3</td><td>1.0e-15</td><td>8</td><td>0.10</td></tr><tr><td>0.2</td><td>5</td><td>4.7e-11</td><td>16</td><td>0.22</td><td>5</td><td>1.8e-13</td><td>16</td><td>0.29</td></tr><tr><td>0.4</td><td>18</td><td>1.2e-11</td><td>16</td><td>0.69</td><td>5</td><td>8.4e-16</td><td>32</td><td>0.85</td></tr></tbody></table>

From the numerical experiments it results that the shift techniques generally provide a more precise approximation of the solution. The reduction of the number of iterations is more evident for problems which are (close to) null recurrent. In certain cases, where the reduction of the number of iterations is stronger, the shifting techniques allow a reduction of the cpu time. In other cases, the reduction of the number of iterations is set off by an increase of the number of interpolation points; this leaves the cpu time almost unchanged. The latter behaviour is an issue to be better investigated.

For null-recurrent problems the convergence of SCR seems to be quadratic. This property has been proved by Guo[9] for QBD probelms. For the M/G/1 case the shifted function Aˆ(z) is such that the zeros η<subs>i</subs>, i = 1, 2,..., of satisfy |η<subs>m</subs>| < 1 and η<subs>m+1</subs> = 1, therefore |η<subs>m</subs>| < |η<subs>m+1</subs>|. This is a necessary condition in order to have quadratic convergence. However, even though we cannot apply the arguments used in section 4 to prove convergence since η<subs>m+1</subs> = 1 is in the unit circle, we believe that this property is valid in general. This is an issue of future investigation.

A different approach in order to deal with the null-recurrent case is to apply a double step of the shift technique as described in section 3.6 of Bini et al.[4]. More precisely, after that ξ<subs>m</subs> = 1 has been shifted to 0, also ξ<subs>m+1</subs> = 1 is shifted away from the unit circle, in this case to infinity. This technique can be applied if A(z) is a matrix polynomial. The double shift technique leads to a nicer splitting of the zeros η<subs>i</subs> of the polynomial det(zI − Aˆ(z)) where |η<subs>n</subs>| = |ξ<subs>m−1</subs>| < 1 < |ξ<subs>m+1</subs>| = |η<subs>m</subs>|, so that the function zI − Aˆ(z) is nonsingular for |z| = 1. This approach is the subject of future research.

ACKNOWLEDGMENTS

The authors wish to thank two anonymous referees who have provided valuable suggestions for improving the presentation of this paper.

By DarioA. Bini; Beatrice Meini and IlyaM. Spitkovsky

Reported by Author; Author; Author