A shifted cyclic reduction algorithm for QBDs

C. He^{∗}  B. Meini^{†}  N.H. Rhee^{‡}

May 9, 2000

Abstract

The problem of the computation of the rate matrix G associated with discrete-time QBD Markov chains is analyzed. We present a shifted cyclic reduction algorithm and show that the speed of convergence of the latter modified algorithm is always faster than that of the original cyclic reduction.

1 Introduction

The block tridiagonal and block Toeplitz structure of the probability transition matrix P associated with QBD problems, i.e.,

P = [ \begin{matrix}B_{0} & B_{1} & & 0 \\ A_{0} & A_{1} & A_{2} \\ & A_{0} & A_{1} & ... \\ 0 & & ... & ... \end{matrix} ] ,

allowed the design of fast and reliable methods for computation of the rate matrix G [16, 11, 2, 3, 4, 12, 5]. Here, A_{0}, A_{1}, A_{2}, B_{0}, B_{1} are k × k nonnegative matrices such that A_{0} + A_{1} + A_{2} is irreducible, and P is irreducible, stochastic and positive recurrent. These fast methods rely on the property that the matrix G, which solves the nonlinear matrix equation

G = A_{0} + A_{1}G + A_{2}G^{2},  (1)

∗Sprint Corporation, Network Planning and Design, 7171 W. 95th Street Overland Park, KS 66212, charlie.he@mail.sprint.com

†Dipartimento di Matematica, Università di Pisa, Italy, meini@dm.unipi.it

‡Department of Mathematics and Statistics, University of Missouri at Kansas City, Kansas City, MO 64110, rhee@cctr.umkc.edu

can be computed by solving the infinite block tridiagonal, block Toeplitz system

[ \begin{matrix}I − A_{1} & −A_{2} & & 0 \\ −A_{0} & I − A_{1} & −A_{2} \\ & −A_{0} & I − A_{1} & ... \\ 0 & & ... & ... \end{matrix} ] [ G ; G^{2} ; G^{3} ; ⋮ ] = [ A_{0} ; 0 ; 0 ; ⋮ ] .  (2)

Recently, in [2, 3, 4], Bini and Meini have devised a new quadratically convergent and numerically stable algorithm for the computation of G, based on a functional representation of cyclic reduction, which applies to general M/G/1 type Markov chains [16] and which extends the method of Latouche and Ramaswami [11].

The aim of this paper is to introduce a shifted cyclic reduction algorithm for QBDs and to show that, under the nonrestrictive assumption that G has only one eigenvalue of modulus one, the speed of convergence of a shifted cyclic reduction algorithm is always faster than that of the original cyclic reduction algorithm.

More precisely, since G has a known eigenvalue (equal to one) and a known corresponding right eigenvector, we may apply a shifting technique, that consists in moving the eigenvalue 1 to 0, and in maintaining the remaining ones. This trick leads to a new quadratic matrix equation, having as solution a singular matrix H, such that G = H + eu^{T}, where e is the vector having all the entries equal to 1, and u is a known arbitrary vector, with positive entries, such that u^{T} e = 1. Thus the problem of the computation of G is reduced to the problem of the computation of H, that shares with G all the eigenvalues, except for the eigenvalue 1, that is moved to 0.

In order to compute H we apply the cyclic reduction algorithm, and we show that the shifting leads to a faster convergence. More precisely, in standard cyclic reduction, the error of the approximation to G, for j → ∞, goes to zero as O ((1/σ)^{2^{j}} ), where σ = min{|λ| : |λ| > 1, ϕ(λ) = 0}, and ϕ(λ) = det(−A_{0} + (I − A_{1})λ − A_{2}λ^{2}). With the shifting technique, the error goes to zero as O ((θ/σ)^{2^{j}} ), where θ is any real number such that max{|λ| : |λ| < 1, ϕ(λ) = 0} < θ < 1. Thus the convergence speed can be much increased; in particular the improvement is more appreciated if σ ≈ 1 and θ << 1.

Finally, we introduce a means to measure the conditioning of the quadratic matrix equations, and we prove that the shifted equation is better conditioned than the original one. Indeed, the shifting technique leads to a better rate of convergence, but destroys the nonnegativity and the M-matrix properties of the blocks generated at each step of standard cyclic reduction, and in principle this fact could lead to a loss of accuracy of the results.

We have performed several numerical experiments, which show that the shifted cyclic reduction algorithm is fast and numerically accurate.

The paper is organized as follows. In Section 2 we recall the cyclic reduction algorithm for QBDs. In Section 3 we apply the shifting technique, and we show how the solutions of the shifted matrix equation are related to the solution of the original one. In Section 4 we analyze the convergence properties of the cyclic reduction algorithm applied to the shifted matrix equation. In Section 5 we study the conditioning of the two matrix equations. In Section 6 we present some numerical results.

2 The cyclic reduction algorithm

In this section we recall the cyclic reduction algorithm for QBDs, described in [3].

Let us consider the system (2). By recursively applying block cyclic reduction, i.e., an odd-even permutation of block rows and block columns, followed by one step of Gaussian elimination, the following sequence of infinite block tridiagonal systems is generated:

[ \begin{matrix}I − Â_{1}^{(j)} & −A_{2}^{(j)} & & 0 \\ −A_{0}^{(j)} & I − A_{1}^{(j)} & −A_{2}^{(j)} \\ & −A_{0}^{(j)} & I − A_{1}^{(j)} & ... \\ 0 & & ... & ... \end{matrix} ] [ G ; G^{2^{j}+1} ; G^{2·2^{j}+1} ; ⋮ ] = [ A_{0} ; 0 ; 0 ; ⋮ ] , j ≥ 0,

(3)

where A_{0}^{(0)} = A_{0}, Â_{1}^{(0)} = A_{1}^{(0)} = A_{1}, A_{2}^{(0)} = A_{2}, and the blocks Â_{1}^{(j)}, A_{i}^{(j)}, i = 0, 1, 2, j ≥ 1, are defined by the recurrences:

\begin{matrix}A_{0}^{(j+1)} = A_{0}^{(j)} (I − A_{1}^{(j)})^{−1} A_{0}^{(j)} \\ A_{1}^{(j+1)} = A_{1}^{(j)} + A_{0}^{(j)} (I − A_{1}^{(j)})^{−1} A_{2}^{(j)} + A_{2}^{(j)} (I − A_{1}^{(j)})^{−1} A_{0}^{(j)} \\ A_{2}^{(j+1)} = A_{2}^{(j)} (I − A_{1}^{(j)})^{−1} A_{2}^{(j)} \\ Â_{1}^{(j+1)} = Â_{1}^{(j)} + A_{2}^{(j)} (I − A_{1}^{(j)})^{−1} A_{0}^{(j)}. \end{matrix}  (4)

The sequences of matrices generated by the above relations allow the fast computation of the matrix G. Indeed, from (3) it follows that

(I − Â_{1}^{(j)})G − A_{2}^{(j)}G^{2^{j}+1} = A_{0}.  (5)

On the other hand, if we denote by R the minimal nonnegative solution of the matrix equation R = A_{2} + RA_{1} + R^{2}A_{0}, then the following equation is verified (see [3, 11])

−A_{2}^{(j)} + R^{2^{j}}(I − A_{1}^{(j)}) − R^{2·2^{j}}A_{0}^{(j)} = 0.

Moreover, the maximum modulus eigenvalue of R is real, simple and unique, and it is equal to 1/σ (see [8]), where σ = min{|λ| : |λ| > 1, ϕ(λ) = 0}, and ϕ(λ) = det(−A_{0} + (I − A_{1})λ − A_{2}λ^{2}). Thus, there exists an operator norm || · || such that ||R|| = ρ(R), where ρ(R) denotes the spectral radius of R; hence, we obtain

||A_{2}^{(j)}|| ≤ ρ(R)^{2^{j}}(||I − A_{1}^{(j)}|| + ρ(R)^{2^{j}}||A_{0}^{(j)}||).

Since the matrices A_{1}^{(j)} and A_{0}^{(j)} are bounded in norm [3], it follows that ||A_{2}^{(j)}|| = O (ρ(R)^{2^{j}} ) and from (5), since G^{2^{j}+1} is bounded, that

(I − Â_{1}^{(j)})G − A_{0} = O ((\frac{1}{σ} )^{2^{j}} ) .  (6)

Hence an approximation of the matrix G is given by (I − Â_{1}^{(j)})^{−1} A_{0}, for a sufficiently large value of j.

Due to the double exponential convergence to zero of the sequence A_{2}^{(j)}, just a small number of steps can be sufficient to reach a good approximation of G.

The rate of convergence is given by 1/σ, and it is therefore related to the closeness to the unit circle of the smallest zero of ϕ(λ) of modulus larger than one.

In the next section we present a trick to improve the rate of convergence; indeed, we show that by applying a deflating technique, that consists in removing the zero λ = 1 of ϕ(λ), the rate of convergence can be reduced to θ/σ, where θ is any real number such that max{|λ| : |λ| < 1, ϕ(λ) = 0} < θ < 1.

3 A shifted matrix equation

In the sequel we assume, without loss of generality, that λ = 1 is the only zero of the function det(−A_{0} + (I − A_{1})λ − A_{2}λ^{2}) = 0 on the unit circle [7].

Since G has a known eigenvalue 1 and a known eigenvector e where e is the k-dimensional vector having all the entries equal to 1, we can modify G so that the eigenvalue 1 is shifted to 0. This shifting technique improves the convergence speed of the cyclic reduction algorithm.

Set H = G−eu^{T}, where u is any vector whose elements are positive and such that u^{T} e = 1. Then the eigenvalues of H are those of G except that in H the eigenvalue 1 of G is replaced by 0. Moreover, e is an eigenvector of H corresponding to the eigenvalue 0, and hence He = 0. So we have

G = H + eu^{T} and G^{2} = (H + eu^{T})^{2} = H^{2} + eu^{T} H + eu^{T} .

By replacing G by H in (1) we obtain that H solves the following shifted equation

B_{0} + B_{1}H + B_{2}H^{2} = H,  (7)

where

\begin{matrix}B_{0} = A_{0} + (A_{1} + A_{2} − I)eu^{T} = A_{0}(I − eu^{T}) \\ B_{1} = A_{1} + A_{2}eu^{T} \\ B_{2} = A_{2}. \end{matrix}  (8)

So the computation of the matrix G can be reduced to the computation of the matrix H = G − eu^{T}, that solves the nonlinear matrix equation (7).

3.1 Spectral properties of the shifted matrix polynomial

The rank one corrections of the matrices A_{0} and A_{1} have the effect to keep unchanged all the zeros of the polynomial ϕ(λ) = det A(λ), A(λ) = −A_{0} + (I −A_{1})λ−A_{2}λ^{2}, except for the zero λ = 1, that is moved to λ = 0. Indeed, the following result holds:

Theorem 1 The zeros of the polynomial ψ(λ) = det B(λ), B(λ) = −B_{0} + (I − B_{1})λ − B_{2}λ^{2}, are

{λ : det A(λ) = 0, λ ≠ 1} ∪ {0}.

To prove the above theorem we need to introduce some notations and some results on pairs of matrices (see, for example, [9]).

Definition 1 Let A and B be n × n complex matrices. We define Λ(A, B) by

Λ(A, B) = {λ ∈ C : det(A − λB) = 0}.

If λ ∈ Λ(A, B), λ is called an eigenvalue of the pair (A, B). If λ is an eigenvalue of the pair (A, B), then there exists a nonzero vector x such that Ax = λBx, and such an x is called an eigenvector of the pair (A, B) corresponding to the eigenvalue λ.

Theorem 2 If A and B are n×n complex matrices, then there exist unitary matrices U and V such that U^{H}AV = R and U^{H}BV = R̃ where R and R̃ are upper triangular. If for some i, r_{ii} and r̃_{ii} are both zero, then Λ(A, B) = C. Otherwise

Λ(A, B) = {r_{ii}/r̃_{ii} : r̃_{ii} ≠ 0, i = 1, …, n}.

Note that v_{1}, the first column of V, is an eigenvector of the pair (A, B) corresponding to the eigenvalue r_{11}/r̃_{11}.

Definition 2 Two pairs of matrices (A, B) and (C, D) are said to be equivalent if there exist nonsingular matrices L and M such that C = LAM and D = LBM. It is easy to verify that if (A, B) and (C, D) are equivalent pairs, then Λ(A, B) = Λ(C, D).

Let us now prove Theorem 1: Proof of Theorem 1. Let

Λ(A) = {λ ∈ C : det A(λ) = 0}, Λ(B) = {λ ∈ C : det B(λ) = 0}.

If we set

S = [ \begin{matrix}I − A_{1} & −A_{0} \\ I & 0 \end{matrix} ]  and  T = [ \begin{matrix}A_{2} & 0 \\ 0 & I \end{matrix} ]

it is easy to verify that Λ(A) = Λ(S, T). Similarly, if we set

S̃ = [ \begin{matrix}I − B_{1} & −B_{0} \\ I & 0 \end{matrix} ]

then Λ(B) = Λ(S̃, T).

Now if we let

L = [ \begin{matrix}I & −A_{2}eu^{T} \\ 0 & I \end{matrix} ]  and  M = [ \begin{matrix}I & eu^{T} \\ 0 & I \end{matrix} ]

then by direct computation one can verify that

S̃ = L(S − T [ \begin{matrix}e \\ e \end{matrix} ] [0, u^{T}])M  and  T = LT M.

So

Λ(B) = Λ(S̃, T) = Λ(S − T [ \begin{matrix}e \\ e \end{matrix} ] [0, u^{T}], T).

Note that the vector [e^{T} , e^{T}]^{T} is an eigenvector of the pair (S, T) corresponding to the eigenvalue 1. So there exist unitary matrices U and V such that

U^{H}SV = R  and  U^{H}T V = R̃

with

v_{1} = \frac{1}{\sqrt{2k}} [ \begin{matrix}e \\ e \end{matrix} ] .

Since v_{1} is an eigenvector of the pair (S, T) corresponding to the eigenvalue 1, we have r_{11}/r̃_{11} = 1. By Theorem 2 we know that

\begin{aligned}Λ(A) & = & Λ(S, T) \\ & = & \{\frac{r_{11}}{r̃_{11}} , \frac{r_{22}}{r̃_{22}} , …, \frac{r_{l,l}}{r̃_{l,l}} \} \\ & = & \{1, \frac{r_{22}}{r̃_{22}} , …, \frac{r_{l,l}}{r̃_{l,l}} \}. \end{aligned}  (9)

Here we assumed that r̃_{i,i} = 0 for i = l + 1, …, 2k.

Since U^{H}T [e^{T} , e^{T}]^{T} = R̃V^{H}(\sqrt{2k}v_{1}) = \sqrt{2k}r̃_{11}e_{1}, where e_{1} is the first standard unit vector of length 2k, and the first component of the row vector [0, u^{T}]V is 1/\sqrt{2k}, it follows that the matrix

U^{H} \{S − T [ \begin{matrix}e \\ e \end{matrix} ] [0, u^{T}]\} V = R − U^{H}T [ \begin{matrix}e \\ e \end{matrix} ] [0, u^{T}]V

is upper triangular and has r_{11} −r̃_{11}, r_{22}, …, r_{2k,2k} as its diagonal elements. Hence by Theorem 2 we have

\begin{aligned}Λ(B) & = & Λ(S̃, T) \\ & = & Λ(S − T [ \begin{matrix}e \\ e \end{matrix} ] [0, u^{T}], T) \\ & = & \{\frac{r_{11}−r̃_{11}}{r̃_{11}} , \frac{r_{22}}{r̃_{22}} , …, \frac{r_{l,l}}{r̃_{l,l}} \} \\ & = & \{0, \frac{r_{22}}{r̃_{22}} , …, \frac{r_{l,l}}{r̃_{l,l}} \}. \end{aligned}  (10)

From (9) and (10) we see that Λ(B) is the same as Λ(A) except that in Λ(B) 1 ∈ Λ(A) is replaced by 0.

3.2 Solutions of the shifted matrix equation

The blocks B_{0}, B_{1}, B_{2} are such that H = G − eu^{T} solves the nonlinear matrix equation (7). The following theorem shows how the minimal nonnegative solutions of the nonlinear matrix equations (11) are transformed by the shifting:

Theorem 3 Let u be positive, and let G, R, S and F be the minimal nonnegative solutions of the matrix equations

\begin{matrix}G = A_{0} + A_{1}G + A_{2}G^{2} \\ S = A_{0}S^{2} + A_{1}S + A_{2} \\ R = R^{2}A_{0} + RA_{1} + A_{2} \\ F = A_{0} + F A_{1} + F^{2}A_{2}, \end{matrix}  (11)

respectively. Then H = G − eu^{T}, K = (I − eu^{T} S)S(I − eu^{T} S)^{−1}, T = R and V = F − wu^{T} (I − A_{1} − A_{0}S)^{−1}, where w = (1 − u^{T} Se)^{−1}A_{0}(I − S)e, solve the matrix equations:

\begin{matrix}H = B_{0} + B_{1}H + B_{2}H^{2} \\ K = B_{0}K^{2} + B_{1}K + B_{2} \\ T = T^{2}B_{0} + T B_{1} + B_{2} \\ V = B_{0} + V B_{1} + V^{2}B_{2}. \end{matrix}  (12)

Proof. First observe that if u is positive then u^{T} Se < 1, since Se ≤ e and S is substochastic [10, 14], hence I − eu^{T} S is nonsingular.

The matrix H solves the first equation of (12) by construction. For the matrix K we have:

B_{0}K^{2} + B_{1}K + B_{2} =

B_{0}(I − eu^{T} S)S^{2}(I − eu^{T} S)^{−1} + B_{1}(I − eu^{T} S)S(I − eu^{T} S)^{−1} + B_{2} =

(A_{0}(I − eu^{T})S^{2} + (A_{1} + (A_{0}e − e)u^{T} S + A_{2}eu^{T} ) S+
+A_{2}(I − eu^{T} S)) (I − eu^{T} S)^{−1} =
(S − eu^{T} S^{2})(I − eu^{T} S)^{−1} = K.

For the matrix T = R, since [10] R = A_{2}(I − A_{1} − A_{2}G)^{−1} and

A_{2}e = R(I − A_{1} − RA_{0})e = R (I − A_{1} − A_{2}(I − A_{1} − A_{2}G)^{−1}A_{0}) e =
R(I − A_{1} − A_{2})e = RA_{0}e,

we obtain

B_{2} + RB_{1} + R^{2}B_{0} = A_{2} + R(A_{1} + A_{2}eu^{T}) + R^{2}A_{0}(I − eu^{T}) =
A_{2} + RA_{1} + R^{2}A_{0} + RA_{2}eu^{T} − R^{2}A_{0}eu^{T} =
R + R(A_{2}e − RA_{0}e)u^{T} = R.

Concerning the matrix V, observe that I − B_{1} − B_{0}K is nonsingular, since I − A_{1} − A_{0}S − B_{2}G is nonsingular [10, 14], and I − A_{1} − A_{0}S − B_{2}G = (I − B_{1} − B_{0}K − B_{2}H)(I − eu^{T} S) = (I − B_{1} − B_{0}K)(I − KH)(I − eu^{T} S). It can be shown by direct substitution that V = B_{0}(I − B_{1} − B_{0}K)^{−1} solves the last matrix equation in (12). By rewriting the matrix V in terms of the blocks A_{0}, A_{1}, A_{2} we obtain that:

V = A_{0}(I − eu^{T})(I − A_{1} − A_{0}S − A_{2}eu^{T})^{−1}.

On the other hand it holds:

(I − A_{1} − A_{0}S − A_{2}eu^{T})^{−1} =
(I − A_{1} − A_{0}S)^{−1} +

(1 − u^{T} (I − A_{1} − A_{0}S)^{−1}A_{2}e)^{−1}(I − A_{1} − A_{0}S)^{−1}A_{2}eu^{T} (I − A_{1} − A_{0}S)^{−1} =

(I − A_{1} − A_{0}S)^{−1} + (1 − u^{T} Se)^{−1}Seu^{T} (I − A_{1} − A_{0}S)^{−1} =
(I − (1 − u^{T} Se)^{−1}Seu^{T} ) (I − A_{1} − A_{0}S)^{−1}.

From the latter equation and from the relation F = A_{0}(I − A_{1} − A_{0}S)^{−1} [10], we obtain:

V = F + (1 − u^{T} Se)^{−1}A_{0}Seu^{T} (I − A_{1} − A_{0}S)^{−1} −
− (A_{0}e + (u^{T} Se)(1 − u^{T} Se)^{−1}A_{0}e) u^{T} (I − A_{1} − A_{0}S)^{−1} =
F − wu^{T} (I − A_{1} − A_{0}S)^{−1},

where w = (1 − u^{T} Se)^{−1}A_{0}(I − S)e.

The spectral properties of the solutions of the shifted matrix equations allow us to prove that the rate of convergence is improved, with respect to standard cyclic reduction. Let us define

\begin{matrix}θ̄ = max\{|λ| : |λ| < 1, ϕ(λ) = 0\} \\ σ = min\{|λ| : |λ| > 1, ϕ(λ) = 0\}. \end{matrix}  (13)

The convergence properties that will be proved in the next section rely on the fact that θ̄ < 1, σ = 1/ρ(R) > 1, S and R have the same eigenvalues (see [14]), and thus also K and R, and ρ(K) = 1/σ < 1.

4 A shifted cyclic reduction algorithm

In this section we apply the cyclic reduction algorithm to solve the shifted matrix equation (7), and we show that the rate of convergence is improved, with respect to the same algorithm applied to the original matrix equation (1).

By following the approach described in section 2, we generate by means of cyclic reduction the sequence of infinite block tridiagonal systems

[ \begin{matrix}I − B̂_{1}^{(j)} & −B_{2}^{(j)} & & 0 \\ −B_{0}^{(j)} & I − B_{1}^{(j)} & −B_{2}^{(j)} \\ & −B_{0}^{(j)} & I − B_{1}^{(j)} & ... \\ 0 & & ... & ... \end{matrix} ] [ H ; H^{2^{j}+1} ; H^{2·2^{j}+1} ; ⋮ ] = [ B_{0} ; 0 ; 0 ; ⋮ ] , j ≥ 0,

(14)

where B_{0}^{(0)} = B_{0}, B̂_{1}^{(0)} = B_{1}^{(0)} = B_{1}, B_{2}^{(0)} = B_{2}, and the blocks B̂_{1}^{(j)}, B_{i}^{(j)}, i = 0, 1, 2, j ≥ 1, are defined by the recurrences:

\begin{matrix}B_{0}^{(j+1)} = B_{0}^{(j)} (I − B_{1}^{(j)})^{−1} B_{0}^{(j)} \\ B_{1}^{(j+1)} = B_{1}^{(j)} + B_{0}^{(j)} (I − B_{1}^{(j)})^{−1} B_{2}^{(j)} + B_{2}^{(j)} (I − B_{1}^{(j)})^{−1} B_{0}^{(j)} \\ B_{2}^{(j+1)} = B_{2}^{(j)} (I − B_{1}^{(j)})^{−1} B_{2}^{(j)} \\ B̂_{1}^{(j+1)} = B̂_{1}^{(j)} + B_{2}^{(j)} (I − B_{1}^{(j)})^{−1} B_{0}^{(j)}. \end{matrix} (15)

In the shifted case, the problem of the nonsingularity of the blocks I − B_{1}^{(j)} must be considered. In fact, in the original case, the matrices I − A_{1}^{(j)} are nonsingular M-matrices for any j, since I − A_{1}^{(j)} can be viewed as a Schur complement of the block (2^{j+1} − 1) × (2^{j+1} − 1) matrix obtained by truncating the infinite matrix (2) at the block size 2^{j+1} − 1 (see [1, 6]), and this finite matrix is a nonsingular M-matrix. Analogously, the matrix I − B_{1}^{(j)} can be viewed as a Schur complement of the matrix Q_{2^{j+1}−2}, where Q_{n} is the (n + 1) × (n + 1) block matrix

Q_{n} = [ \begin{matrix}I − B_{1} & −B_{2} & & 0 \\ −B_{0} & ... & ... \\ & ... & ... & −B_{2} \\ 0 & & −B_{0} & I − B_{1} \end{matrix} ] .  (16)

Thus, the (j +1)-st step of cyclic reduction can be performed, i.e., I −B_{1}^{(j)} is nonsingular, if and only if Q_{2^{j+1}−2} is nonsingular [1]. Based on this property, we prove the following result:

Theorem 4 Let u be any positive vector such that u^{T} e = 1 and let H = G − ue^{T}, K = (I − eu^{T} S)S(I − eu^{T} S)^{−1}, where G and S are the minimal nonnegative solutions of the matrix equations G = A_{0} + A_{1}G + A_{2}G^{2}, S = A_{0}S^{2} + A_{1}S + A_{2}. Let B_{0}^{(j)}, B_{1}^{(j)}, B_{2}^{(j)}, B̂_{1}^{(j)} be the blocks generated at the j-th step of cyclic reduction, j ≥ 1. Then, the following conditions are equivalent:

1. I − (HK^{2^{j+1}−1})(KH^{2^{j+1}−1}) is nonsingular;

2. I − (KH^{2^{j+1}−1})(HK^{2^{j+1}−1}) is nonsingular;

3. I − B_{1}^{(j)} is nonsingular.

Proof. The matrix I − B_{1}^{(j)} is nonsingular if and only if Q_{2^{j+1}−2} is singular, where Q_{n} is defined in (16). For simplicity of notations, set n = 2^{j+1}−2. Observe that the matrix W = I − B_{1} − B_{0}K − B_{2}H is nonsingular, since W = (I − A_{1} − A_{0}S − A_{2}G)(I − eu^{T} S) and I − A_{1} − A_{0}S − B_{2}G is nonsingular (see [10, 14]). Let x = (x_{i})_{i=0,...,n} such that Q_{n}x = 0. Therefore, since I − B_{1} − B_{0}K − B_{2}H is nonsingular, then

x_{i} = H^{i}r + K^{n−i}s, i = 0, …, n

where r and s are suitable vectors, such that the boundary conditions

\begin{matrix}(I − B_{1})x_{0} − B_{2}x_{1} = 0 \\ −B_{0}x_{n−1} + (I − B_{1})x_{n} = 0 \end{matrix}  (17)

are satisfied (see [15]). By imposing the above equalities, we obtain that r and s must solve the following homogeneous linear system:

[ \begin{matrix}I − B_{1} − B_{2}H & ((I − B_{1})K − B_{2}) K^{n−1} \\ ((I − B_{1})H − B_{0}) H^{n−1} & I − B_{1} − B_{0}K \end{matrix} ] [ \begin{matrix}r \\ s \end{matrix} ] = 0,

that can be written as

[ \begin{matrix}I − B_{1} − B_{2}H & B_{0}K^{n+1} \\ B_{2}H^{n+1} & I − B_{1} − B_{0}K \end{matrix} ] [ \begin{matrix}r \\ s \end{matrix} ] = 0,  (18)

since H and K solve the matrix equations (1), (7). Hence, x = 0 if and only if r = s = 0. Now, the block diagonal entries in the matrix of (18) are nonsingular since I − B_{1} − B_{0}K − B_{2}H = (I − B_{1} − B_{2}H)(I − HK) = (I − B_{1} − B_{0}K)(I − KH), and I − B_{1} − B_{0}K − B_{2}H is nonsingular. By computing the Schur complement of I − B_{1} − B_{2}H in I − B_{1} − B_{0}K we obtain the matrix

S_{1} = I − B_{1} − B_{0}K − B_{2}H^{n+1}(I − B_{1} − B_{2}H)^{−1}B_{0}K^{n+1} =
I − B_{1} − B_{0}K − B_{2}H^{n+2}K^{n+1}.

Thus, S_{1} (and hence Q_{n}) is singular if and only if (I − B_{1} − B_{0}K)^{−1}S_{1} = I − KH^{n+2}K^{n+1} is nonsingular. By taking the Schur complement of I − B_{1} − B_{0}K in I − B_{1} − B_{0}K we complete the proof of the theorem.

The above theorem gives necessary and sufficient conditions for the applicability of cyclic reduction. Since ρ(H) < 1 and ρ(K) < 1, the matrix Z_{j} = I − (HK^{2^{j+1}−1})(KH^{2^{j+1}−1}) is nonsingular for sufficiently large values of j. From the numerical experiments that we have performed, and from the fact that I − (GS^{2^{j+1}−1})(SG^{2^{j+1}−1}) is nonsingular for any value of j, we conjecture that also Z_{j} is nonsingular for any value of j. However, if it were not, the cyclic reduction algorithm can be still applied, by performing a different permutation of block rows and columns, that allows one to skip the steps that cannot be performed for the singularity of the blocks I − B_{1}^{(j)} (we refer to the paper [1] for details of this subject).

Here and hereafter we assume that the matrices I − B_{1}^{(j)} are nonsingular for any j.

The nice convergence properties stated by the next theorem enable us to efficiently compute the matrix H:

Theorem 5 Let B_{0}^{(j)}, B_{1}^{(j)}, B_{2}^{(j)} be the blocks generated at the j-th step of cyclic reduction. Then B_{1}^{(j)}, (I − B_{1}^{(j)})^{−1} are bounded, and for any operator norm it holds ||B_{2}^{(j)}|| = O ((1/σ)^{2^{j}} ) and ||B_{0}^{(j)}|| = O (θ^{2^{j}} ), for any θ̄ < θ < 1, where θ̄ and σ are defined in (13).

Proof. From (14), at each step j the following equations hold:

H^{2^{j}} = B_{0}^{(j)} + B_{1}^{(j)}H^{2^{j}} + B_{2}^{(j)}H^{2·2^{j}},  (19)
K^{2^{j}} = B_{0}^{(j)}K^{2·2^{j}} + B_{1}^{(j)}K^{2^{j}} + B_{2}^{(j)},  (20)

and therefore

(I − B_{1}^{(j)})^{−1}B_{0}^{(j)} + H^{2^{j}} + (I − B_{1}^{(j)})^{−1}B_{2}^{(j)}H^{2·2^{j}} = 0,  (21)
(I − B_{1}^{(j)})^{−1}B_{2}^{(j)} + K^{2^{j}} + (I − B_{1}^{(j)})^{−1}B_{0}^{(j)}K^{2·2^{j}} = 0.  (22)

We first prove that the matrices C_{0}^{(j)} = (I − B_{1}^{(j)})^{−1}B_{0}^{(j)} and C_{2}^{(j)} = (I − B_{1}^{(j)})^{−1}B_{2}^{(j)} are bounded in norm.

Since the maximum modulus eigenvalue of R, and thus of K, is real, simple and unique [8] then there exists an operator norm || · ||_{K} such that ||K||_{K} = ρ(K). Moreover, let ϵ > 0 such that ρ(H) + ϵ < 1, and let || · ||_{H,ϵ} be an operator norm such that ||H||_{H,ϵ} ≤ ρ(H) + ϵ (see [9]).

Let α_{j} = ||C_{0}^{(j)}||_{H,ϵ}. If α_{j} is not bounded, then there exists a subsequence α_{j_{h}} such that α_{j_{h}} diverges to infinity. From (21) we have

α_{j_{h}} ≤ ||C_{2}^{(j_{h})}||_{H,ϵ}||H^{2·2^{j_{h}}} ||_{H,ϵ} + ||H^{2^{j_{h}}} ||_{H,ϵ} ≤
(ρ(H) + ϵ)^{2^{j_{h}}} ((ρ(H) + ϵ)^{2^{j_{h}}} ||C_{2}^{(j_{h})}||_{H,ϵ} + 1) .

Thus,

||C_{2}^{(j_{h})}||_{H,ϵ} ≥ (\frac{α_{j_{h}}}{(ρ(H) + ϵ)^{2^{j_{h}}}} − 1) \frac{1}{(ρ(H) + ϵ)^{2^{j_{h}}}}

and, since α_{j_{h}} diverges to infinity, there exists a constant c > 0 such that

||C_{2}^{(j_{h})}||_{H,ϵ} ≥ c\frac{α_{j_{h}}}{(ρ(H) + ϵ)^{2^{j_{h}}}} .

For the equivalence of the operator norms, there exist constants c′ > 0 and c′′ > 0 such that

||C_{2}^{(j_{h})}||_{K} ≥ c′ \frac{α_{j_{h}}}{(ρ(H) + ϵ)^{2^{j_{h}}}} ,  (23)
||C_{0}^{(j)}||_{K} ≤ c′′||C_{0}^{(j)}||_{H,ϵ} = c′′α_{j}.  (24)

On the other hand, from (22) we have

||C_{2}^{(j_{h})}||_{K} ≤ ||C_{0}^{(j_{h})}||_{K}||K^{2·2^{j_{h}}} ||_{K} + ||K^{2^{j_{h}}} ||_{K}.

Thus, from the latter inequality and (24) we have

||C_{2}^{(j_{h})}||_{K} ≤ c′′α_{j_{h}}ρ(K)^{2·2^{j_{h}}} + ρ(K)^{2^{j_{h}}} .

Hence, from (23), we obtain

c′ \frac{α_{j_{h}}}{(ρ(H) + ϵ)^{2^{j_{h}}}} ≤ c′′α_{j_{h}}ρ(K)^{2·2^{j_{h}}} + ρ(K)^{2^{j_{h}}} ,

that contradicts the assumption that α_{j_{h}} goes to infinity. By using a similar argument we can prove that also C_{2}^{(j)} is bounded in norm. From (21), we obtain

||C_{0}^{(j)}||_{H,ϵ} ≤ (ρ(H) + ϵ)^{2^{j}} + ||C_{2}^{(j)}||_{H,ϵ}(ρ(H) + ϵ)^{2·2^{j}}.

Hence, since ||C_{2}^{(j)}||_{H,ϵ} is bounded, then there exists a constant γ > 0 such that

||C_{0}^{(j)}||_{H,ϵ} ≤ γ(ρ(H) + ϵ)^{2^{j}}.

Thus, for the equivalence of the operator norms, for any operator norm || · || there exists a constant γ such that

||C_{0}^{(j)}|| ≤ γ(ρ(H) + ϵ)^{2^{j}}.  (25)

Similarly, for any operator norm || · || there exists a constant γ′ such that

||C_{2}^{(j)}|| ≤ γ′ρ(K)^{2^{j}}.  (26)

From (25) and (26), since

I − B_{1}^{(j+1)} = (I − B_{1}^{(j)})(I − C_{0}^{(j)}C_{2}^{(j)} + C_{2}^{(j)}C_{0}^{(j)}),

we have

||I − B_{1}^{(j+1)}|| ≤ ||I − B_{1}^{(j)}||(1 + σ_{j})

where σ_{j} = O ((ρ(H) + ϵ)^{2^{j}}ρ(K)^{2^{j}} ). Thus the matrices I − B_{1}^{(j)}, and hence

the matrices B_{1}^{(j)}, are bounded in norm. Similarly, it holds

||(I − B_{1}^{(j+1)})^{−1}|| ≤ \frac{1}{1 − σ_{j}} ||(I − B_{1}^{(j)})^{−1}||,

thus ||(I − B_{1}^{(j)})^{−1}|| is bounded.

Now, from the boundness of ||I − B_{1}^{(j)}||, and from the relations

−B_{0}^{(j)} + (I − B_{1}^{(j)})H^{2^{j}} − B_{2}^{(j)}H^{2·2^{j}} = 0,
−B_{2}^{(j)} + (I − B_{1}^{(j)})K^{2^{j}} − B_{0}^{(j)}K^{2·2^{j}} = 0,

derived from (19) and (20), by applying the same argument used to derive (25) and (26), we can show that ||B_{2}^{(j)}|| and ||B_{0}^{(j)}|| are bounded, and thus that ||B_{0}^{(j)}|| = O ((ρ(H) + ϵ)^{2^{j}} ), ||B_{2}^{(j)}|| = O (ρ(K)^{2^{j}} ).

The matrix G can be directly recovered, without computing the matrix H, according to the following result:

Theorem 6 For any operator norm || · || and for any θ̄ < θ < 1 it holds

A_{0} − (I − B̂_{1}^{(j)})G = O ((\frac{θ}{σ} )^{2^{j}} ) .  (27)

Proof. From the relation

(I − B̂_{1}^{(j)})H − B_{2}^{(j)}H^{2^{j}+1} − B_{0} = 0

and from Theorem 5 it follows that for any operator norm || · || and for any θ̄ < θ < 1 it holds

B_{0} − (I − B̂_{1}^{(j)})H = O ((\frac{θ}{σ} )^{2^{j}} ) .

On the other hand, it can be easily proved by induction that (I − B̂_{1}^{(j)})e = A_{0}e for any j ≥ 0. Thus, by replacing H with G − eu^{T} and B_{0} with A_{0} − A_{0}eu^{T}, we arrive at (27).

From the above result, compared with (6), it follows that the shifted cyclic reduction can be much faster than the original one. Indeed, if the second largest modulus eigenvalue θ̄ of G is far from the unit circle, then the rate of convergence is much improved.

Thus, the deflating technique leads to a better rate of convergence, but destroys the nonnegativity and M-matrix properties of the blocks generated at each step. Indeed, in general neither I − B_{1}^{(j)} is an M-matrix, nor B_{0}^{(j)}, B_{2}^{(j)} are nonnegative matrices (from the numerical experiments it seems that I − B̂_{1}^{(j)} is an M-matrix for j = 0, 1, …). In principle this fact could lead to a loss of accuracy of the results obtained with the shifting technique. In practice, we have not observed any differences, in terms of accuracy, between the results obtained with the two algorithms. Furthermore, as we will prove in the next section, the shifted equation is better conditioned than the original one.

5 Conditioning of the shifted matrix equation

In this section we introduce a measure of the conditioning of the matrix equation, and we show that the shifted equation is better conditioned than the original one.

Consider the matrix equation (1) and a perturbed equation

G+△G = (A_{0}+△A_{0})+(A_{1}+△A_{1})(G+△G)+(A_{2}+△A_{2})(G+△G)^{2}. (28)

Using (1) this perturbed equation simplifies up to the first oder in △G

(I − A_{1} − A_{2}G)△G − A_{2}(△G)G = △A,  (29)

where

△A = △A_{0} + (△A_{1})G + (△A_{2})G^{2}.

Note that (29) can be written

W vec(△G) = vec(△A),  (30)

where

W = I ⊗ (I − A_{1} − A_{2}G) + G^{T} ⊗ (−A_{2})

and vec(A) is the k^{2}-dimensional vector obtained by column-wise arranging the entries of the matrix A.

Let Λ(M) denote the set of all eigenvalues of any square matrix M.

Theorem 7 The k^{2} × k^{2} matrix W is nonsingular. Furthermore

min{|λ| : λ ∈ Λ(W)} = 1 − ρ(A_{1} + A_{2}G + A_{2}).

Proof. Let S be the Schur canonical form of G^{T}. Then the matrix W is similar to

W̃ = I ⊗ (I − A_{1} − A_{2}G) + S ⊗ (−A_{2}).

Due to the upper triangular structure of S and the fact that Λ(S) = Λ(G) it is easy to see that

Λ(W) = Λ(W̃) = ⋃_{λ∈Λ(G)} Λ(I − A_{1} − A_{2}G − λA_{2}).

Note that whenever λ ∈ Λ(G)

ρ(A_{1} + A_{2}G + λA_{2}) ≤ ρ(|A_{1} + A_{2}G + λA_{2}|)
≤ ρ(A_{1} + A_{2}G + |λ|A_{2})
≤ ρ(A_{1} + A_{2}G + A_{2})
< 1.

In the last step we used the fact that since A_{0} + A_{1} + A_{2} is irreducible and positive recurrent the matrix A_{1} + A_{2}G + A_{2} has spectral radius less than 1. It follows that 0 does not belong to Λ(W), and hence W is nonsingular. Clearly

min{|λ| : λ ∈ Λ(W)} = 1 − ρ(A_{1} + A_{2}G + A_{2})

Since W is nonsingular

vec(△G) = W^{−1}vec(△A).

Hence

∥△G∥_{F} = ∥vec(△G)∥_{2}
≤ ∥W^{−1}∥_{2}∥vec(△A)∥_{2}
= \frac{1}{σ_{min}(W)}∥△A_{0} + (△A_{1})G + (△A_{2})G^{2}∥_{F}
≤ \frac{\sqrt{k}}{σ_{min}(W)}(∥△A_{0}∥_{F} + ∥△A_{1}∥_{F} + ∥△A_{2}∥_{F}).

Here σ_{min}(W) is the minimum singular value of the matrix W. So we may view 1/σ_{min}(W) as a condition number of the equation (1). Even though 1/ min{|λ| : λ ∈ Λ(W)} ≤ 1/σ_{min}(W), we may consider

\frac{1}{min\{|λ| : λ ∈ Λ(W)\}} = \frac{1}{1 − ρ(A_{1} + A_{2}G + A_{2})}

as a number which reflects the conditioning of the equation (1).

Consider now the shifted matrix equation (7) and the perturbed equation

H + △H = (B_{0} + △B_{0}) + (B_{1} + △B_{1})(H + △H) + (B_{2} + △B_{2})(H + △H)^{2}.

As before this perturbed equation simplifies up to the first oder in △H

(I − B_{1} − B_{2}H)△H − B_{2}(△H)H = △B,  (31)

where

△B = △B_{0} + (△B_{1})H + (△B_{2})H^{2}.

Note that (31) can be written

Qvec(△H) = vec(△B),

where

Q = I ⊗ (I − B_{1} − B_{2}H) + H^{T} ⊗ (−B_{2}).

Note that

B_{1} + B_{2}H = A_{1} + A_{2}eu^{T} + A_{2}(G − eu^{T}) = A_{1} + A_{2}G.

Hence

Q = I ⊗ (I − A_{1} − A_{2}G) + H^{T} ⊗ (−A_{2}).

Note that if we denote Λ(G) by

Λ(G) = {λ_{1}, λ_{2}, …, λ_{k}}

with

1 = λ_{1} > |λ_{2}| ≥ … ≥ |λ_{k}|

then

Λ(H) = {λ_{2}, λ_{3}, …, λ_{k}, 0}

and |λ_{2}| = θ̄. Theorem 8 The k^{2} × k^{2} matrix Q is nonsingular. Furthermore

min{|λ| : λ ∈ Λ(Q)} ≥ 1 − ρ(A_{1} + A_{2}G + θ̄A_{2})

and the equality holds if λ_{2} is a positive real number.

Proof. As in the proof of Theorem 7 it is easy to see that

Λ(Q) = Λ(Q̃) = ⋃_{λ∈Λ(H)} Λ(I − A_{1} − A_{2}G − λA_{2}).

Note that since ρ(H) = θ̄ = |λ_{2}| whenever λ ∈ Λ(H)

ρ(A_{1} + A_{2}G + λA_{2}) ≤ ρ(|A_{1} + A_{2}G + λA_{2}|)
≤ ρ(A_{1} + A_{2}G + |λ|A_{2})
≤ ρ(A_{1} + A_{2}G + |λ_{2}|A_{2})
< 1.

In the last step we used the fact that

ρ(A_{1} + A_{2}G + |λ_{2}|A_{2}) < ρ(A_{1} + A_{2}G + A_{2}) < 1.

It follows that 0 does not belong to Λ(Q), and hence Q is nonsingular. Clearly

min{|λ| : λ ∈ Λ(Q)} ≥ 1 − ρ(A_{1} + A_{2}G + |λ_{2}|A_{2})

and the equality holds if λ_{2} is a positive real number.

Since Q is nonsingular

vec(△H) = Q^{−1}vec(△B).

Hence

∥△H∥_{F} = ∥vec(△H)∥_{2}
≤ ∥Q^{−1}∥_{2}∥vec(△B)∥_{2}
= \frac{1}{σ_{min}(Q)}∥△B_{0} + (△B_{1})H + (△B_{2})H^{2}∥_{F}
≤ \frac{2\sqrt{k}}{σ_{min}(Q)}(∥△B_{0}∥_{F} + ∥△B_{1}∥_{F} + ∥△B_{2}∥_{F}).

Here σ_{min}(Q) is the minimum singular value of the matrix Q. So we may view 1/σ_{min}(Q) as a condition number of the shifted equation (7). Even though 1/ min{|λ| : λ ∈ Λ(Q)} ≤ 1/σ_{min}(Q), we may consider 1/ min{|λ| : λ ∈ Λ(Q)} as a number which reflects the conditioning of the shifted equation (7). Since 1/ min{|λ| : λ ∈ Λ(Q)} ≤ 1/(1 − ρ(A_{1} + A_{2}G + θ̄A_{2})),

\frac{1}{1 − ρ(A_{1} + A_{2}G + θ̄A_{2})}

is a number which reflects the conditioning of the shifted equation (7). The inequality

\frac{1}{1 − ρ(A_{1} + A_{2}G + θ̄A_{2})} < \frac{1}{1 − ρ(A_{1} + A_{2}G + A_{2})},

suggests that the shifted equation (7) has a better conditioning than the original equation (1).

6 Numerical results

We have tested the cyclic reduction algorithm and the shifted cyclic reduction algorithm on the examples in [13] using Matlab. In the case of the cyclic reduction algorithm we stopped when

∥Â_{1}^{(j)} − Â_{1}^{(j−1)}∥_{∞} ≤ 10^{−12}

and we accepted

G̃ = (I − Â_{1}^{(j)})^{−1} A_{0}

as an approximation of G. Similarly, in the case of a shifted cyclic reduction we stopped when

∥B̂_{1}^{(j)} − B̂_{1}^{(j−1)}∥_{∞} ≤ 10^{−12}

Cyclic Reduction  Shifted Cyclic Reduction

\begin{matrix}p \\ 10^{−1} \\ 10^{−2} \\ 10^{−3} \\ 10^{−4} \\ 10^{−5} \\ 10^{−6} \\ 10^{−7} \\ 10^{−8} \\ 10^{−9} \\ 10^{−10} \end{matrix} Iter. Res.  Stoc.

9 0.0  6.7 · 10^{−16}
13 0.0  7.9 · 10^{−15}
16 0.0  6.0 · 10^{−15}
20 1.1 · 10^{−16} 3.4 · 10^{−13}
23 2.2 · 10^{−16} 1.3 · 10^{−11}
26 2.2 · 10^{−16} 4.9 · 10^{−11}
30 2.2 · 10^{−16} 3.1 · 10^{−10}
33 0.0  7.9 · 10^{−9}
36 1.1 · 10^{−16} 3.0 · 10^{−8}
39 0.0  2.8 · 10^{−7}

Iter. Res.  Stoc.
2 0.0  2.2 · 10^{−16}
2 0.0  1.1 · 10^{−16}
2 0.0  0.0
2 0.0  0.0
2 2.2 · 10^{−16} 0.0
2 2.2 · 10^{−16} 1.1 · 10^{−16}
2 5.6 · 10^{−17} 2.2 · 10^{−16}
2 0.0  0.0
2 5.6 · 10^{−17} 1.1 · 10^{−16}
2 0.0  0.0

Table 1: Example 1

and we accepted

G̃ = (I − B̂_{1}^{(j)})^{−1} A_{0}

as an approximation of G. For each example we have reported tables with the number of iterations, the residual error, Res., which is defined by ∥G̃ − A_{0} − A_{1}G̃ − A_{2}G̃^{2}∥_{∞}, and the closeness to stochasticity, Stoc., which is defined by ∥G̃e − e∥_{∞}, for standard and shifted cyclic reduction. For particular examples we have also reported tables with the values of θ̄, of σ, 1/σ and θ̄/σ, to show the reduction of the rate of convergence, and tables with 1 − ρ(A_{1} + A_{2}G + A_{2}), 1 − ρ(A_{1} + A_{2}G + θ̄A_{2}), σ_{min}(W), σ_{min}(Q), that show the conditioning of the original and shifted matrix equations.

Example 1 Here the blocks, having dimension 2, depend on a parameter p > 0, and are defined as

A_{0} = [ \begin{matrix}1 − p & 0 \\ 0 & 0 \end{matrix} ] , A_{1} = [ \begin{matrix}0 & p \\ 2p & 0 \end{matrix} ] , A_{2} = [ \begin{matrix}0 & 0 \\ 0 & 1 − 2p \end{matrix} ]

We have tested with ten different p values and obtained the following results summarized in Tables 1, 2, and 3.

Example 2 In this example, A_{0}, A_{1}, and A_{2} are defined as follows. First we define 24 × 24 matrices A′_{0}, A′_{1}, and A′_{2} as follows:

(A′_{0}) = \{ \begin{matrix}192(1 − i/24), & i = j \\ 0, & i ≠ j \end{matrix} (A′_{2}) = \{ \begin{matrix}192ρ_{d}, & i = j \\ 0, & i ≠ j \end{matrix}

\begin{matrix}p & θ̄ & σ & 1/σ & θ̄/σ \\ 10^{−1} & 0 & 1.125 & 0.8889 & 0 \\ 10^{−2} & 0 & 1.0102 & 0.9899 & 0 \\ 10^{−3} & 0 & 1.001 & 0.9990 & 0 \\ 10^{−4} & 0 & 1.0001 & 0.9999 & 0 \\ 10^{−5} & 0 & 1.00001 & 0.99999 & 0 \end{matrix}

Table 2: Example 1

\begin{matrix}p & 1 − ρ(A_{1} + A_{2}G + A_{2}) & 1 − ρ(A_{1} + A_{2}G + θ̄A_{2}) & σ_{min}(W) & σ_{min}(Q) \\ 10^{−1} & 0.09 & 0.68 & 0.05 & 0.35 \\ 10^{−6} & 1.0 · 10^{−6} & 0.999 & 5.0 · 10^{−7} & 0.30 \\ 10^{−10} & 1.0 · 10^{−10} & 1.0 & 5.0 · 10^{−11} & 0.44 \end{matrix}

Table 3: Example 1

and

(A′_{1}) = \{ \begin{matrix}ar(M − i)/M, & i − j = −1 \\ ir, & i − j = 1 \\ 0, & elsewhere \end{matrix}

Here i and j are integers between 0 and 23, and a, r, M, and ρ_{d} are parameters. Now A_{0}, A_{1}, and A_{2} are defined such that A_{0} = sA′_{0}, A_{1} = sA′_{1}, and A_{2} = sA′_{2} where s is a scalar which makes A_{0} + A_{1} + A_{2} stochastic. Tables 4, 5, and 6 report the results obtained by choosing different values of M while we fix r = 1/300, a = 18.244, and ρ_{d} = 0.280. Tables 7, 8, and 9 report the results obtained by choosing different values of ρ_{d} while we fix r = 1/100, a = 18.244, and M = 512.

Example 3 In this example we construct a QBD problem defined by the k × k matrices A_{0} = R + δI, A_{1} = A_{2} = R, where R is a matrix having null diagonal entries and constant off-diagonal entries, and 0 < δ < 1. As was observed in [13], the rate ρ = p^{T} (A_{1} + 2A_{2})e, where p^{T} (A_{0} + A_{1} + A_{2}) = p^{T}, p^{T} e = 1, is exactly 1−δ. We have tested with eight different δ values. Tables 10, 11, and 12 report the results obtained with k = 16 and Tables 13 and 14 report the results obtained with sizes k = 32 and k = 64, respectively.

Cyclic Reduction  Shifted Cyclic Reduction

\begin{matrix}M \\ 64 \\ 128 \\ 256 \\ 512 \\ 1024 \\ 2048 \\ 4096 \\ 8192 \\ 16384 \\ 32768 \\ 65536 \end{matrix} Iter. Res.  Stoc.

19 1.6 · 10^{−16} 3.2 · 10^{−12}
20 2.2 · 10^{−16} 8.4 · 10^{−13}
21 3.0 · 10^{−16} 5.9 · 10^{−12}
22 2.2 · 10^{−16} 1.4 · 10^{−12}
23 2.4 · 10^{−16} 2.1 · 10^{−11}
24 3.0 · 10^{−16} 6.1 · 10^{−11}
25 2.6 · 10^{−16} 5.4 · 10^{−11}
26 3.1 · 10^{−16} 3.9 · 10^{−10}
27 2.2 · 10^{−16} 2.0 · 10^{−11}
29 1.3 · 10^{−16} 1.1 · 10^{−9}
34 2.2 · 10^{−16} 5.5 · 10^{−8}

Iter. Res.  Stoc.
18 5.1 · 10^{−16} 4.4 · 10^{−16}
19 6.9 · 10^{−16} 4.4 · 10^{−16}
19 4.6 · 10^{−16} 4.4 · 10^{−16}
19 4.7 · 10^{−16} 4.4 · 10^{−16}
19 5.3 · 10^{−16} 5.6 · 10^{−16}
19 5.0 · 10^{−16} 4.4 · 10^{−16}
19 5.0 · 10^{−17} 6.7 · 10^{−16}
19 4.8 · 10^{−16} 6.7 · 10^{−16}
19 5.9 · 10^{−17} 6.7 · 10^{−16}
19 3.5 · 10^{−16} 4.4 · 10^{−16}
19 5.0 · 10^{−16} 6.7 · 10^{−16}

Table 4: Example 2: r = 1/300, a = 18.244, ρ_{d} = 0.280

\begin{matrix}M & θ̄ & σ & 1/σ & θ̄/σ \\ 64 & 0.9999 & 1.0002 & 0.9998 & 0.9997 \\ 128 & 0.9999 & 1.0001 & 0.9999 & 0.9998 \\ 256 & 0.9999 & 1.000039 & 0.999961 & 0.9998 \\ 512 & 0.9999 & 1.000019 & 0.999981 & 0.9998 \\ 1024 & 0.9999 & 1.000009 & 0.999991 & 0.9998 \\ 2048 & 0.9999 & 1.000004 & 0.999996 & 0.9998 \\ 4096 & 0.9999 & 1.000002 & 0.999998 & 0.9999 \\ 8192 & 0.9999 & 1.000001 & 0.999999 & 0.9999 \\ 16384 & 0.9999 & 1.0000004 & 0.9999996 & 0.9999 \\ 32768 & 0.9999 & 1.0000001 & 0.9999999 & 0.9999 \\ 65536 & 0.9999 & 1.000000002 & 0.999999998 & 0.9999 \end{matrix}

Table 5: Example 2: r = 1/300, a = 18.244, ρ_{d} = 0.280

\begin{matrix}M & 1 − ρ(A_{1} + A_{2}G + A_{2}) & 1 − ρ(A_{1} + A_{2}G + θ̄A_{2}) & σ_{min}(W) & σ_{min}(Q) \\ 64 & 1.4 · 10^{−4} & 2.2 · 10^{−4} & 2.7 · 10^{−5} & 3.8 · 10^{−5} \\ 8192 & 7.1 · 10^{−7} & 1.1 · 10^{−4} & 1.5 · 10^{−7} & 1.5 · 10^{−5} \\ 65536 & 1.6 · 10^{−9} & 1.0 · 10^{−4} & 3.5 · 10^{−10} & 1.5 · 10^{−5} \end{matrix}

Table 6: Example 2: r = 1/300, a = 18.244, ρ_{d} = 0.280

Cyclic Reduction  Shifted Cyclic Reduction

\begin{matrix}ρ_{d} \\ 0.01 \\ 0.025 \\ 0.05 \\ 0.075 \\ 0.1 \\ 0.12 \\ 0.14 \\ 0.16 \\ 0.18 \\ 0.2 \\ 0.22 \\ 0.24 \\ 0.26 \\ 0.28 \\ 0.29 \\ .29568 \end{matrix} Iter. Res.  Stoc.

6 2.2 · 10^{−16} 4.4 · 10^{−16}
7 3.3 · 10^{−16} 6.6 · 10^{−16}
10 1.5 · 10^{−16} 8.9 · 10^{−16}
12 2.2 · 10^{−16} 9.0 · 10^{−15}
14 2.2 · 10^{−16} 2.2 · 10^{−15}
14 2.2 · 10^{−16} 2.8 · 10^{−14}
15 2.2 · 10^{−16} 2.1 · 10^{−13}
16 2.2 · 10^{−16} 1.5 · 10^{−14}
16 2.2 · 10^{−16} 5.9 · 10^{−14}
17 2.2 · 10^{−16} 1.0 · 10^{−14}
18 2.2 · 10^{−16} 4.3 · 10^{−13}
18 2.3 · 10^{−16} 7.8 · 10^{−13}
19 1.9 · 10^{−16} 8.6 · 10^{−13}
20 2.3 · 10^{−16} 3.5 · 10^{−12}
22 2.6 · 10^{−16} 1.8 · 10^{−11}
31 2.6 · 10^{−16} 1.9 · 10^{−8}

Iter. Res.  Stoc.
6 3.9 · 10^{−16} 2.2 · 10^{−16}
7 5.7 · 10^{−16} 4.4 · 10^{−16}
10 3.0 · 10^{−16} 4.4 · 10^{−16}
12 6.1 · 10^{−16} 4.4 · 10^{−16}
14 4.2 · 10^{−16} 4.4 · 10^{−16}
14 3.9 · 10^{−16} 4.4 · 10^{−16}
15 5.1 · 10^{−17} 4.4 · 10^{−16}
16 5.3 · 10^{−16} 4.4 · 10^{−16}
16 4.7 · 10^{−17} 3.3 · 10^{−16}
16 4.5 · 10^{−16} 6.7 · 10^{−16}
17 6.0 · 10^{−16} 3.3 · 10^{−16}
17 3.7 · 10^{−16} 6.7 · 10^{−16}
17 5.0 · 10^{−16} 4.4 · 10^{−16}
17 5.1 · 10^{−16} 3.3 · 10^{−16}
17 3.7 · 10^{−16} 2.2 · 10^{−16}
17 6.0 · 10^{−16} 4.4 · 10^{−16}

Table 7: Example 2: r = 1/100, a = 18.244, M = 512

\begin{matrix}ρ_{d} & θ̄ & σ & 1/σ & θ̄/σ \\ 0.01 & 0.9998 & 4.3182 & 0.2316 & 0.2315 \\ 0.025 & 0.9998 & 1.7714 & 0.5645 & 0.5644 \\ 0.05 & 0.9998 & 1.0814 & 0.9248 & 0.9246 \\ 0.075 & 0.9998 & 1.0170 & 0.9833 & 0.9831 \\ 0.1 & 0.9998 & 1.0063 & 0.9938 & 0.9936 \\ 0.12 & 0.9998 & 1.0034 & 0.9966 & 0.9964 \\ 0.14 & 0.9997 & 1.0021 & 0.9979 & 0.9977 \\ 0.16 & 0.9997 & 1.0013 & 0.9987 & 0.9984 \\ 0.18 & 0.9997 & 1.0009 & 0.9991 & 0.9988 \\ 0.2 & 0.9997 & 1.0006 & 0.9994 & 0.9991 \\ 0.22 & 0.9997 & 1.0004 & 0.9996 & 0.9993 \\ 0.24 & 0.9996 & 1.0002 & 0.9998 & 0.9994 \\ 0.26 & 0.9996 & 1.0001 & 0.9999 & 0.9995 \\ 0.28 & 0.9996 & 1.0001 & 0.9999 & 0.9995 \\ 0.29 & 0.9996 & 1.000019 & 0.99998 & 0.9995 \\ 0.29568 & 0.9995 & 1.00000003 & 0.99999997 & 0.9995 \end{matrix}

Table 8: Example 2: r = 1/100, a = 18.244, M = 512

\begin{matrix}p & 1 − ρ(A_{1} + A_{2}G + A_{2}) & 1 − ρ(A_{1} + A_{2}G + θ̄A_{2}) & σ_{min}(W) & σ_{min}(Q) \\ 0.01 & 0.6265 & 0.6265 & 0.6259 & 0.6175 \\ 0.28 & 4.0 · 10^{−5} & 3.4 · 10^{−4} & 8.7 · 10^{−6} & 5.3 · 10^{−5} \\ 0.29568 & 1.9 · 10^{−8} & 3.2 · 10^{−4} & 3.6 · 10^{−9} & 4.4 · 10^{−5} \end{matrix}

Table 9: Example 2: r = 1/100, a = 18.244, M = 512

Cyclic Reduction  Shifted Cyclic Reduction

\begin{matrix}δ \\ 10^{−1} \\ 10^{−2} \\ 10^{−3} \\ 10^{−4} \\ 10^{−5} \\ 10^{−6} \\ 10^{−7} \\ 10^{−8} \end{matrix} Iter. Res.  Stoc.

8 5.7 · 10^{−16} 1.3 · 10^{−15}
11 6.5 · 10^{−16} 7.2 · 10^{−15}
14 6.7 · 10^{−16} 9.7 · 10^{−14}
17 6.4 · 10^{−16} 8.8 · 10^{−13}
21 1.1 · 10^{−15} 1.2 · 10^{−11}
24 6.5 · 10^{−16} 5.8 · 10^{−11}
27 6.3 · 10^{−16} 1.4 · 10^{−9}
29 7.0 · 10^{−16} 3.5 · 10^{−9}

Iter. Res.  Stoc.
5 2.8 · 10^{−16} 3.3 · 10^{−16}
4 3.7 · 10^{−16} 0.0
4 1.9 · 10^{−16} 2.2 · 10^{−16}
4 4.5 · 10^{−16} 2.2 · 10^{−16}
4 3.5 · 10^{−16} 4.4 · 10^{−16}
4 2.5 · 10^{−16} 8.9 · 10^{−16}
4 2.5 · 10^{−16} 6.7 · 10^{−16}
4 2.5 · 10^{−16} 4.4 · 10^{−16}

Table 10: Example 3: k = 16

\begin{matrix}δ & θ̄ & σ & 1/σ & θ̄/σ \\ 10^{−1} & 0.0783 & 1.3333 & 0.75 & 0.0587 \\ 10^{−2} & 0.0174 & 1.0303 & 0.9706 & 0.0140 \\ 10^{−3} & 0.0207 & 1.0030 & 0.9970 & 0.0207 \\ 10^{−4} & 0.0216 & 1.0003 & 0.9997 & 0.0216 \\ 10^{−5} & 0.0217 & 1.00003 & 0.99997 & 0.0217 \end{matrix}

Table 11: Example 3: k = 16

\begin{matrix}δ & 1 − ρ(A_{1} + A_{2}G + A_{2}) & 1 − ρ(A_{1} + A_{2}G + θ̄A_{2}) & σ_{min}(W) & σ_{min}(Q) \\ 10^{−1} & 0.1 & 0.3765 & 0.1 & 0.3234 \\ 10^{−5} & 1.0 · 10^{−5} & 0.3261 & 1.0 · 10^{−5} & 0.2310 \\ 10^{−8} & 1.0 · 10^{−8} & 0.3261 & 1.0 · 10^{−8} & 0.2721 \end{matrix}

Table 12: Example 3: k = 16

Cyclic Reduction  Shifted Cyclic Reduction

\begin{matrix}δ \\ 10^{−1} \\ 10^{−2} \\ 10^{−3} \\ 10^{−4} \\ 10^{−5} \\ 10^{−6} \\ 10^{−7} \\ 10^{−8} \end{matrix} Iter. Res.  Stoc.

8 1.3 · 10^{−15} 1.1 · 10^{−15}
11 9.5 · 10^{−16} 1.1 · 10^{−14}
14 1.4 · 10^{−15} 8.0 · 10^{−14}
17 1.6 · 10^{−15} 3.7 · 10^{−12}
21 1.1 · 10^{−15} 1.0 · 10^{−11}
24 1.2 · 10^{−15} 2.1 · 10^{−10}
27 1.1 · 10^{−15} 7.9 · 10^{−10}
29 8.2 · 10^{−16} 1.4 · 10^{−8}

Iter. Res.  Stoc.
5 1.2 · 10^{−15} 8.9 · 10^{−16}
4 4.1 · 10^{−16} 7.8 · 10^{−16}
4 5.4 · 10^{−16} 1.1 · 10^{−15}
4 6.8 · 10^{−16} 1.8 · 10^{−15}
4 5.1 · 10^{−16} 6.6 · 10^{−16}
4 5.5 · 10^{−16} 8.9 · 10^{−16}
4 6.5 · 10^{−16} 1.0 · 10^{−15}
4 7.1 · 10^{−16} 7.8 · 10^{−16}

Table 13: Example 3: k = 32

Cyclic Reduction  Shifted Cyclic Reduction

\begin{matrix}δ \\ 10^{−1} \\ 10^{−2} \\ 10^{−3} \\ 10^{−4} \\ 10^{−5} \\ 10^{−6} \\ 10^{−7} \\ 10^{−8} \end{matrix} Iter. Res.  Stoc.

8 2.8 · 10^{−15} 4.9 · 10^{−15}
11 3.1 · 10^{−15} 1.1 · 10^{−14}
14 3.1 · 10^{−15} 1.8 · 10^{−13}
17 2.4 · 10^{−15} 4.6 · 10^{−12}
21 2.3 · 10^{−15} 2.9 · 10^{−11}
24 3.3 · 10^{−15} 2.8 · 10^{−10}
27 2.4 · 10^{−15} 1.9 · 10^{−10}
29 2.6 · 10^{−15} 1.7 · 10^{−8}

Iter. Res.  Stoc.
5 2.0 · 10^{−15} 1.8 · 10^{−15}
4 1.4 · 10^{−15} 2.7 · 10^{−15}
4 1.3 · 10^{−15} 1.1 · 10^{−15}
4 1.2 · 10^{−15} 2.9 · 10^{−15}
4 1.6 · 10^{−15} 4.4 · 10^{−16}
4 1.4 · 10^{−15} 2.4 · 10^{−15}
4 1.3 · 10^{−15} 3.8 · 10^{−15}
4 1.2 · 10^{−15} 3.3 · 10^{−15}

Table 14: Example 3: k = 64

References

[1] D. A. Bini, L. Gemignani, and B. Meini. Factorization of analytic functions by means of Koenig’s theorem and Toeplitz computations. Numer. Math., 2000. To appear.

[2] D. A. Bini and B. Meini. On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems. In W. J. Stewart, editor, Computations with Markov Chains, pages 21–38. Kluwer Academic Publisher, 1995.

[3] 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.

[4] D. A. Bini and B. Meini. Improved cyclic reduction for solving queueing problems. Numerical Algorithms, 15:57–74, 1997.

[5] D. A. Bini and B. Meini. Inverting block Toeplitz matrices in block Hessenberg form by means of displacement operators: application to queueing problems. Linear Algebra Appl., 272:1–16, 1998.

[6] D. A. Bini and B. Meini. Effective methods for solving banded Toeplitz systems. SIAM J. Matrix Anal. Appl., 20:700–719, 1999.

[7] H. R. Gail, S. L. Hantler, and B. A. Taylor. Spectral analysis of M/G/1 and G/M/1 type Markov chains. Adv. in Appl. Probab., 28:114–165, 1996.

[8] 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.

[9] G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, 1989.

[10] G. Latouche. A note on two matrices occurring in the solution of quasibirth-and-death processes. Commun. Statist. Stochastic Models, 3:251– 257, 1987.

[11] G. Latouche and V. Ramaswami. A logarithmic reduction algorithm for Quasi-Birth-Death processes. J. Appl. Probability, 30:650–674, 1993.

[12] G. Latouche and G.W. Stewart. Numerical methods for M/G/1 type queues. In W. J. Stewart, editor, Computations with Markov Chains, pages 571–581. Kluwer Academic Publishers, 1995.

[13] B. Meini. Solving QBD problems: the cyclic reduction algorithm versus the invariant subspace method. Adv. Perf. Anal., 1:215–225, 1998.

[14] V. Naoumov. Matrix-multiplicative approach to quasi-birth-and-death processes analysis. In A. S. Alfa and S. Chakravarthy, editors, Matrix-analytic methods in stochastic models (Flint, MI), pages 85– 106. Dekker, New York, 1997.

[15] V. Naoumov, U. R. Krieger, and D. Wagner. Analysis of a multiserver delay-loss system with a general Markovian arrival process. In A. S. Alfa and S. Chakravarthy, editors, Matrix-analytic methods in stochastic models (Flint, MI), pages 43–66. Dekker, New York, 1997.

[16] M. F. Neuts. Structured Stochastic Matrices of M/G/1 Type and Their Applications. Marcel Dekk., New York, 1989.