Solving M/G/1 type Markov chains: recent advances and applications

Beatrice Meini Dipartimento di Matematica Università di Pisa Via Buonarroti 2 56127 Pisa, Italy E-mail: meini@dm.unipi.it

Abstract

A survey on numerical methods, based on Toeplitz matrix computation, for the solution of M/G/1 type Markov chains is presented together with some new advances. We show how the block Toeplitz structure, shared by the transition matrices associated with M/G/1 type Markov chains, can be exploited as the basis for devising fast and numerically stable algorithms, for the computation of the probability invariant vector and for the solution of the nonlinear matrix equation X = ∑_{i=0}^{+∞} A_{i}X^{i}.

Key words. M/G/1 type matrices, Toeplitz matrices, Fast Fourier Transform.

1 INTRODUCTION

Discrete time Markov chains of M/G/1 type are defined by a transition probability matrix of the form

P = [ \begin{matrix}B_{1} & B_{2} & B_{3} & … \\ A_{0} & A_{1} & A_{2} & … \\ & A_{0} & A_{1} & ... \\ ◯ & & ... & ... \end{matrix} ]  (1)

where A_{i}, B_{i+1}, i ≥ 0, are m × m nonnegative matrices such that ∑_{i=0}^{+∞} A_{i} and ∑_{i=1}^{+∞} B_{i} are row stochastic.

These types of Markov chains, introduced by M.F. Neuts [23] model a very large variety of queueing problems. Due to the strong structure of P, it is possible to give an elegant representation of the probability invariant vector and of other performance measures of interest [23], in spite of the infinite size of the matrix.

A specific feature of M/G/1 type Markov chains is the block Hessenberg form of the matrix P and the block Toeplitz structure of the submatrix obtained by removing the first block row of P. Here we recall that a block matrix has a block Toeplitz structure if its blocks are constant along each block diagonal.

More general M/G/1 type Markov chains (non-skip-free) are defined by transition matrices of the form

P = [ \begin{matrix}B_{1,1} & B_{1,2} & B_{1,3} & … \\ ⋮ & ⋮ & ⋮ \\ B_{k,1} & B_{k,2} & B_{k,3} & … \\ A_{0} & A_{1} & A_{2} & … \\ & A_{0} & A_{1} & ... \\ ◯ & & ... & ... \end{matrix} ]  (2)

where A_{i}, B_{j,i+1}, i ≥ 0, j = 1, …, k, are m × m nonnegative matrices such that ∑_{i=0}^{+∞} A_{i} and ∑_{i=1}^{+∞} B_{j,i}, j = 1, …, k, are row stochastic.

The latter class of M/G/1 type Markov chains can be reduced to the block Hessenberg form (1) by partitioning the matrix P into blocks A_{i}, B_{i} of dimension p = mk, where

A_{0} = [ \begin{matrix}A_{0} & A_{1} & … & A_{k−1} \\ 0 & A_{0} & ... & ⋮ \\ ⋮ & ... & ... & A_{1} \\ 0 & … & 0 & A_{0} \end{matrix} ] , A_{i} = [ \begin{matrix}A_{ik} & A_{ik+1} & … & A_{ik+k−1} \\ A_{ik−1} & A_{ik} & ... & ⋮ \\ ⋮ & ... & ... & A_{ik+1} \\ A_{ik−k+1} & … & A_{ik−1} & A_{ik} \end{matrix} ] , i ≥ 1, (3)

B_{i} = [ \begin{matrix}B_{1,(i−1)k+1} & B_{1,(i−1)k+2} & … & B_{1,ik} \\ B_{2,(i−1)k+1} & B_{2,(i−1)k+2} & … & B_{2,ik} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ B_{k,(i−1)k+1} & B_{k,(i−1)k+2} & … & B_{k,ik} \end{matrix} ] , i ≥ 1.

A fundamental step in solving Markov chains of M/G/1 type is computing the stochastic matrix G that solves the nonlinear matrix equation

X = ∑_{i=0}^{+∞} A_{i}X^{i}  (4)

where X is an m × m matrix. The knowledge of G allows us to compute the probability invariant vector π^{T} = [π_{0}^{T} , π_{1}^{T} , …], partitioned according to the block structure of (1), by means of the recursive formula (Ramaswami’s formula [25]):

π_{i} = (I − Ā_{1}^{T})^{−1} [B̄_{i}^{T} π_{0} + ∑_{j=1}^{i−1} Ā_{i+1−j}^{T}π_{j}] , i ≥ 1  (5)

where

Ā_{i} = ∑_{j=i}^{+∞} A_{j}G^{j−i}, B̄_{i} = ∑_{j=i}^{+∞} B_{j}G^{j−i}, for i ≥ 1;  (6)

the vector π_{0} is the solution of the system

\{ \begin{matrix}B̄_{1}^{T} π_{0} = π_{0} \\ (e^{T} + e^{T} (I − ∑_{i=1}^{+∞} Ā_{i}^{T})^{−1} ∑_{i=2}^{+∞} B̄_{i}^{T} ) π_{0} = 1. \end{matrix}

In many applications, the dimension m of the blocks A_{i} can be very large. For this reason it is important to design efficient algorithms for solving the nonlinear matrix equation (4), for computing Ramaswami’s formula (5), and for solving the more general case of non-skip-free M/G/1 type Markov chains.

Recently, the analysis of numerical methods for solving the nonlinear matrix equation (4) has been developed by several authors. Classical linearly convergent methods based on functional iterations, proposed in [24], [16], have been studied in detail in [22], where an accurate analysis of the speed of convergence is provided. An adaptation of Newton’s method is proposed in [17]. An efficient quadratically convergent algorithm for solving QBD problems is presented in [18]. A doubling algorithm that solves a sequence of nested block Toeplitz block Hessenberg systems is analyzed in [19]. A method based on the use of invariant subspaces is proposed in [1]. In the case of non-skip-free M/G/1 type Markov chains, an elegant analysis of the structure of the matrix G associated with the reblocked transition matrix (2), is provided in [13]; the authors observe that the matrix G is uniquely defined by a k × k block Frobenius matrix with blocks of size m, that is uniquely defined by the first block row g of G; in order to evaluate g an iterative technique, providing a sequence g_{n} of block vectors that monotonically converges to g, is also proposed.

In this paper we show how the block Toeplitz structure of M/G/1 type Markov chains can be exploited as the basis for devising fast and numerically stable algorithms, for the computation of the probability invariant vector and for the solution of the nonlinear matrix equation (4), also in the general non-skip-free case. We show how the strong connection between polynomials and Toeplitz matrices allows us to devise a fast version of Ramaswami’s formula, based on the computation of the inverse of an infinite block triangular block Toeplitz matrix, by means of Fast Fourier Transforms (FFT’s). We present an improvement of the doubling algorithm of Latouche and Stewart [19] for solving (4), based on the concept of displacement rank: due to the Toeplitz structure of the problem, we show that only few block vectors need to be computed in order to solve the block Hessenberg block Toeplitz systems at each step of the algorithm, thus substantially reducing the computational cost. The classical cyclic reduction method is rephrased in functional form; this formulation allows us to devise a quadratically convergent and numerically stable algorithm, based on the use of FFT, for solving (4). An adaptation of this algorithm to QBD problems, leading to a method similar to that of Latouche and Ramaswami [18], is presented. Finally we consider the more general case of a non-skip-free M/G/1 type Markov chain; we present a new quadratically convergent algorithm based on cyclic reduction and on the concept of displacement rank; the properties of displacement rank allow us to fully exploit the specific Toeplitz structure of the blocks A_{i} of (3), thus dramatically reducing the computational cost of the cyclic reduction, while preserving the nice features of the quadratic convergence. Numerical results are presented.

The paper is organized as follows. In Section 2 we recall some basic tools for Toeplitz matrix computations. In Section 3 we present the fast version of Ramaswami’s formula. In Section 4 we analyze the quadratically convergent algorithm for solving (4) based on the doubling technique of [19] and on cyclic reduction. In Section 5 we consider the case of non-skip-free M/G/1 type Markov chains. In Section 6 we present numerical results. In Appendix A, details on fast polynomial arithmetic are provided.

2 TOEPLITZ MATRIX TOOLS

In this section we recall the definition of Toeplitz matrices and the basic properties used for devising fast algorithms for computing the vector π and for computing the matrix G.

A block (possibly infinite) matrix A is a block Toeplitz matrix if its block entries are constant along each block diagonal.

2.1 TOEPLITZ MATRICES AND POLYNOMIALS

Matrix polynomials, matrix power series and block Toeplitz matrices are intimately related. The formal relation among the coefficients of the matrix power series P (z) = ∑_{i=0}^{+∞} P_{i}z^{i}, Q(z) = ∑_{i=0}^{+∞} Q_{i}z^{i}, and R(z) = ∑_{i=0}^{+∞} R_{i}z^{i}, where R(z) = P (z)Q(z) and P_{i}, Q_{i}, R_{i}, i ≥ 0, are m × m matrices, can be easily rewritten in matrix form as

[ \begin{matrix}R_{0} & & & O \\ R_{1} & R_{0} \\ R_{2} & R_{1} & R_{0} \\ ⋮ & ... & ... & ... \end{matrix} ] = [ \begin{matrix}P_{0} & & & O \\ P_{1} & P_{0} \\ P_{2} & P_{1} & P_{0} \\ ⋮ & ... & ... & ... \end{matrix} ] [ \begin{matrix}Q_{0} & & & O \\ Q_{1} & Q_{0} \\ Q_{2} & Q_{1} & Q_{0} \\ ⋮ & ... & ... & ... \end{matrix} ] .

In this way computations involving block triangular block Toeplitz matrices can be reduced to manipulating formal matrix power series (for infinite matrices) or polynomials (for finite matrices), and viceversa. Moreover, if P (z) is a matrix polynomial of degree n − 1 and Q(z) is a matrix polynomial of degree n/2 − 1, then, by comparing the terms of degree n/2, …, n − 1 in the equation R(z) = P (z)Q(z) mod (z^{n} − 1), we arrive at the following equation involving a block Toeplitz matrix

[ R_{n/2} ; ⋮ ; R_{n−1} ] = [ \begin{matrix}P_{n/2} & … & P_{1} \\ ⋮ & ... & ⋮ \\ P_{n−1} & … & P_{n/2} \end{matrix} ] [ Q_{0} ; ⋮ ; Q_{n/2−1} ] .

Thus, the product between a finite block Toeplitz matrix and a block vector can be reduced to multiplying two matrix polynomials. Similarly, the inverse of a block triangular Toeplitz matrix is a block triangular Toeplitz matrix, uniquely defined by its first block column, that can be computed by inverting a matrix polynomial. In fact, comparing the coefficients of the term of degree i in both sides of the equation P (z)Q(z) = I mod z^{n} yields the matrix equation involving a block lower triangular block Toeplitz matrix

[ \begin{matrix}P_{0} & & & ◯ \\ P_{1} & P_{0} \\ ⋮ & ... & ... \\ P_{n−1} & … & P_{1} & P_{0} \end{matrix} ] [ Q_{0} ; Q_{1} ; ⋮ ; Q_{n−1} ] = [ I ; O ; ⋮ ; O ] .  (7)

In Appendix A we describe fast algorithms based on FFT for computing products and inverses of matrix polynomials. There, it is shown that due to this connection with matrix polynomial arithmetic, we may perform fast block Toeplitz matrix computation by using the technique of evaluation/interpolation at the roots of 1, based on FFT. More precisely, the computational cost needed to compute the product of an n×n block Toeplitz matrix and a block vector, or the product between two n × n block lower triangular block Toeplitz matrices, reduces to O(m^{3}n + m^{2}n log n) arithmetic operations; similarly, the same bound holds for the computational cost needed to compute the inverse of an n × n block Toeplitz matrix. Due to the use of FFT’s, block Toeplitz matrix computation is weakly numerically stable [20]; we say that an algorithm is weakly numerically stable if the relative rounding errors in the output can be only bounded in norm.

2.2 DISPLACEMENT RANK

The concept of displacement rank, introduced by Kailath et al. in [15], is fundamental in devising and analyzing algorithms related to Toeplitz matrices [9], [8], [3].

Define the n × n block down-shift matrix

Z_{n} = [ \begin{matrix}O \\ I & O \\ & ... & ... \\ & & I & O \end{matrix} ] ,

where the blocks have dimension m, and consider the following block displacement operator

Δ(H) = Z_{n}H − HZ_{n}

defined for any n × n block matrix H.

Let us also introduce the following notations: the block vectors E_{n}^{(1)} = [I, O, …, O]^{T}, E_{n}^{(2)} = [O, …, O, I]^{T} denote mn×m matrices made up by the first and the last m columns of the identity matrix of order mn. Throughout we call by the name n-dimensional block column vector any mn × m matrix; similarly we call by the name n-dimensional block row vector any m × mn matrix. Moreover we say that the block matrix H has block displacement rank r if r is the minimum integer such that there exist block column vectors U^{(i)}, and block row vectors V^{(i)}, i = 1, …, r, satisfying the equation Δ(H) = ∑_{i=1}^{r} U^{(i)}V^{(i)}.

It can be easily observed that, if H is a block Toeplitz matrix, than H has block displacement rank 2. We call by the name block Toeplitz-like matrix an n × n block matrix having block displacement rank r which is independent of the block dimension n.

A nice property of Δ(H) is that the block displacement of the inverse matrix can be explicitly related to the block displacement of the matrix itself, in fact, if H is nonsingular, then

Δ(H^{−1}) = −H^{−1}Δ(H)H^{−1}.  (8)

From this property it follows in particular that the inverse of a block Toeplitz, or Toeplitzlike, matrix is a block Toeplitz-like matrix.

The following result can be easily proved by extending to block matrices the same proof given in [9].

Theorem 1 Let K_{n} be an n × n block matrix such that Δ(K_{n}) = ∑_{i=1}^{r} U^{(i)}V^{(i)}, where U^{(i)} and V^{(i)} are n-dimensional block column and block row vectors, respectively. Then we have

K_{n} = L_{n}(K_{n}E_{n}^{(1)}) − ∑_{i=1}^{r} L_{n}(U^{(i)})L_{n}^{T}(Z_{n}V^{(i)T})

where L_{n}(W) denotes the n × n block lower triangular block Toeplitz matrix defined by its first block column W, and L_{n}^{T}(W) is its transposed matrix.

The above result allows one to represent any matrix K_{n} as a sum of products of block lower and upper triangular block Toeplitz matrices defined by the first block column of K_{n} and by the block vectors U^{(i)}, V^{(i)} associated with the block displacement of K_{n}.

If the matrix K_{n} is nonsingular then the above representation theorem can be applied to K_{n}^{−1} in the light of (8) yielding the following result

K_{n}^{−1} = L_{n}(K_{n}^{−1}E_{n}^{(1)}) + ∑_{i=1}^{r} L_{n}(K_{n}^{−1}U^{(i)})L_{n}^{T}(Z_{n}K_{n}^{−T} V^{(i)T}).  (9)

From the above results it follows that a block Toeplitz-like matrix, and its inverse, can be expressed by means of a sum of few products between block lower triangular and block upper triangular block Toeplitz matrices. This fact allows us to perform block Toeplitzlike matrix computations in terms of a few multiplications between matrix polynomials, in the light of the results of Section 2.1, with a substantial reduction of the computational cost.

3 FAST RAMASWAMI FORMULA

By means of the natural isomorphism between formal power series and triangular Toeplitz matrices outlined in Section 2.1, it is possible to give a novel interpretation of Ramaswami’s formula and to devise a new fast algorithm, based on the use of FFT techniques, for its computation. We refer the reader to [21] for details.

We show that the block Hessenberg matrix

H = [ \begin{matrix}I − A_{1}^{T} & −A_{0}^{T} & & ◯ \\ −A_{2}^{T} & I − A_{1}^{T} & −A_{0}^{T} \\ −A_{3}^{T} & −A_{2}^{T} & I − A_{1}^{T} & −A_{0}^{T} \\ ⋮ & ⋮ & ... & ... & ... \end{matrix} ]  (10)

can be factorized as the product H = UL, where

U = [ \begin{matrix}I & −G^{T} & & & ◯ \\ & I & −G^{T} \\ & & I & −G^{T} \\ ◯ & & & ... & ... \end{matrix} ] L = [ \begin{matrix}I − Ā_{1}^{T} & & & ◯ \\ −Ā_{2}^{T} & I − Ā_{1}^{T} \\ −Ā_{3}^{T} & −Ā_{2}^{T} & I − Ā_{1}^{T} \\ ⋮ & ... & ... & ... \end{matrix} ] ,

G is the solution of (4) and Ā_{i} are defined in (6). Ramaswami’s formula (5) can now be derived by solving the system

UL [ π_{1} ; π_{2} ; π_{3} ; ⋮ ] = [ B_{1}^{T} ; B_{2}^{T} ; B_{3}^{T} ; ⋮ ] π_{0}.

In fact, by multiplying the above equation by the left inverse of U, we obtain that

L [ π_{1} ; π_{2} ; π_{3} ; ⋮ ] = [ B̄_{1}^{T} ; B̄_{2}^{T} ; B̄_{3}^{T} ; ⋮ ] π_{0}.  (11)

Hence, by solving the block lower triangular system (11) by forward substitution, we arrive at (5).

The computation of the vector π is thus reduced to solving a block lower triangular block Toeplitz infinite system. In [21] a fast Ramaswami formula, based on the use of FFT for solving the above infinite system, has been devised with a substantial reduction of the asymptotic computational cost. More precisely, the computational cost to calculate q block components of the probability invariant vector π is roughly q(7m^{2} + 40 log n_{0}) arithmetic operations, where n_{0} is such that the blocks A_{i} are negligible for i > n_{0}. The computational cost of the customary Ramaswami formula is 2qn_{0}m^{2} arithmetic operations. Due to the use of FFT’s, the algorithm is weakly numerically stable.

4 COMPUTATION OF THE MATRIX G

In this section we present quadratically convergent methods for computing the matrix G. These methods consists of solving the infinite block Hessenberg block Toeplitz system

[X_{1}, X_{2}, X_{3}, …]H = [A_{0}^{T} , O, O, …]  (12)

where the matrix H is given in (10) and X_{i} are m × m matrices. In fact, it can be shown [12] that the solution of the above system is given by X_{1} = G^{T}, X_{i} = X_{1}^{i}, i ≥ 2.

4.1 A DOUBLING ALGORITHM

Consider the n × n block Toeplitz matrix in block Hessenberg form

H_{n} = [ \begin{matrix}I − A_{1}^{T} & −A_{0}^{T} & & & ◯ \\ −A_{2}^{T} & I − A_{1}^{T} & −A_{0}^{T} \\ ⋮ & ... & ... & ... \\ −A_{n−1}^{T} & ... & ... & I − A_{1}^{T} & −A_{0}^{T} \\ −A_{n}^{T} & −A_{n−1}^{T} & … & −A_{2}^{T} & I − A_{1}^{T} \end{matrix} ]  (13)

obtained by truncating the infinite matrix (10) and consider the finite system

[X_{1}^{(n)}, X_{2}^{(n)}, …, X_{n}^{(n)}]H_{n} = [A_{0}^{T} , O, …, O].  (14)

In [19] the authors observe that the sequence X_{1}^{(n)} monotonically converges, as n tends to infinity, to G^{T}. Hence, they propose a doubling algorithm for solving the finite system (14).

Let us assume for simplicity that n = 2^{q}, for a positive integer q. Partition the matrix H_{n} in the following way

H_{n} = [ \begin{matrix}H\frac{n}{2} & U\frac{n}{2} V T_{\frac{n}{2}} \\ T\frac{n}{2} & H\frac{n}{2} \end{matrix} ]

where

U\frac{n}{2} = [ O ; ⋮ ; O ; I ] , V\frac{n}{2} = [ −A_{0} ; O ; ⋮ ; O ] , T\frac{n}{2} = [ \begin{matrix}−A_{\frac{n}{2} +1}^{T} & −AT_{\frac{n}{2}} & … & −A_{2}^{T} \\ −A_{\frac{n}{2} +2}^{T} & −A_{\frac{n}{2} +1}^{T} & ... & −A_{3}^{T} \\ ⋮ & ... & ... & ⋮ \\ −A_{n}^{T} & −A_{n−1}^{T} & … & −A_{\frac{n}{2} +1}^{T} \end{matrix} ]

and I denotes the m × m identity matrix. By applying the Sherman-Woodbury-Morrison formula [14] to the decomposition

\begin{matrix}H_{n} = [ \begin{matrix}H\frac{n}{2} & O \\ T\frac{n}{2} & H\frac{n}{2} \end{matrix} ] + [ \begin{matrix}O & U\frac{n}{2} V T_{\frac{n}{2}} \\ O & O \end{matrix} ] = S_{n} + M_{n}N_{n}^{T} , \\ M_{n} = [ \begin{matrix}U\frac{n}{2} \\ O \end{matrix} ] , N_{n} = [ \begin{matrix}O \\ V\frac{n}{2} \end{matrix} ] \end{matrix}

we immediately find the following expression for the matrix inverse of H_{n}:

H_{n}^{−1} = S_{n}^{−1} − S_{n}^{−1}M_{n}(I + N_{n}^{T} S_{n}^{−1}M_{n})^{−1}N_{n}^{T} S_{n}^{−1},

where

S_{n}^{−1} = [ \begin{matrix}H−1_{\frac{n}{2}} & O \\ −H−1_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}} & H−1_{\frac{n}{2}} \end{matrix} ] .

By means of the above relations, in [19], [26], a doubling algorithm for the solution of the linear system (12), truncated at a finite system of block dimension n, has been devised, with a computational cost of O(m^{3}n log^{2} n) arithmetic operations. This algorithm can be much simplified by applying the tool of displacement rank of Section 2.2.

In fact, observe that for the matrix H_{n} of (13) it holds

Δ(H_{n}) = − [ I ; O ; ⋮ ; O ] [A_{0}^{T} , O, …, O] + [ O ; ⋮ ; O ; I ] [O, …, O, A_{0}^{T}].  (15)

Using this, the equation (9) can be much simplified by noting that due to the special structure (15) of Δ(H_{n}) we may easily prove that

H_{n}^{−1} = L_{n}(C_{n}^{(1)})L_{n}^{T}(E_{n}^{(1)} − Z_{n}R_{n}^{(1)T} A_{0}) + L_{n}(C_{n}^{(2)})L_{n}^{T}(Z_{n}R_{n}^{(2)T} A_{0})

where C_{n}^{(1)}, C_{n}^{(2)} denote the first and the last block column, respectively, of H_{n}^{−1}, while R_{n}^{(1)}, R_{n}^{(2)} denote the first and the last block row, respectively, of the matrix H_{n}^{−1}.

In this way the inverse of H_{n} is explicitly determined as a sum of products of block triangular block Toeplitz matrices defined by the block A_{0}, by the block rows R_{n}^{(1)}, R_{n}^{(2)}, and the block columns C_{n}^{(1)}, C_{n}^{(2)}.

The Sherman-Morrison-Woodbury formula allows us to relate the block vectors C_{n}^{(1)},

C_{n}^{(2)}, R_{n}^{(1)}, R_{n}^{(2)} defining the matrix H_{n}^{−1} and the block vectors C(1)_{\frac{n}{2}}, C(2)_{\frac{n}{2}}, R(1)_{\frac{n}{2}}, R(2)_{\frac{n}{2}}

defining the matrix H−1_{\frac{n}{2}}. Moreover, since such relations only involve operations between block Toeplitz matrices and block vectors, we may devise an efficient scheme for their implementation based on the use of FFT. In this way we arrive at a weakly numerically stable algorithm for the computation of the inverses of H_{2^{i}}, i = 0, 1, …, q, requiring O(m^{3}n + m^{2}n log n) arithmetic operations, for n = 2^{q}, and a low storage cost; in fact, only four auxiliary block vectors need to be allocated in order to carry out our algorithm. Another interesting feature is that in the PRAM model of parallel computation the algorithm can be carried out in O(log^{2} n) parallel steps using O(m^{3}n) processors.

The following result provides the desired relations among block rows and columns of

H_{n}^{−1} and H−1_{\frac{n}{2}}.

Theorem 2 The first and the last block columns C_{n}^{(1)}, C_{n}^{(2)}, and the first and the last block rows R_{n}^{(1)}, R_{n}^{(2)} of H_{n}^{−1} satisfy the following relations

C_{n}^{(1)} = [ \begin{matrix}C(1)_{\frac{n}{2}} \\ −H−1_{\frac{n}{2}}T\frac{n}{2} C(1)_{\frac{n}{2}} \end{matrix} ] + [ \begin{matrix}C(2)_{\frac{n}{2}} \\ −H−1_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}} \end{matrix} ] P R(1)_{\frac{n}{2}}T\frac{n}{2} C(1)_{\frac{n}{2}}

C_{n}^{(2)} = [ \begin{matrix}O \\ C(2)_{\frac{n}{2}} \end{matrix} ] − [ \begin{matrix}C(2)_{\frac{n}{2}} \\ −H−1_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}} \end{matrix} ] P R(1)_{\frac{n}{2}}E(2)_{\frac{n}{2}}

R_{n}^{(1)} = −A_{0}^{T} [R(1)_{\frac{n}{2}}, O] − A_{0}^{T} R(1)_{\frac{n}{2}}E(2)_{\frac{n}{2}}P [R(1)_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}}, −R(1)_{\frac{n}{2}}]

R_{n}^{(2)} = A_{0}^{T} [R(2)_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}}, −R(2)_{\frac{n}{2}}] + A_{0}^{T} R(2)_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}}P [R(1)_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}}, −R(1)_{\frac{n}{2}}],

where P = (I − A_{0}^{T} R(1)_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}})^{−1}A_{0}^{T}.

4.2 CYCLIC REDUCTION BASED ALGORITHM

In this section we propose a method for the numerical solution of (12), introduced in [5], [7], [6], which is quadratically convergent, weakly numerically stable and, unlike the method of Section 4.1, avoids the need of cutting the infinite system (12) to a finite size. This method, based on the properties of the Schur complement of block Toeplitz matrices that we describe in functional form, requires a small number of arithmetic operations. In fact, the most expensive computations at each stage of our algorithm are reduced to performing operations between matrix power series, and for these computations we may use fast algorithms based on the use of FFT.

Our method relies on the technique of successive state reduction (cyclic reduction [10]), applied to the infinite system (12).

By performing an odd-even permutation of the block rows and block columns in (12) we find that

[G^{T} , G^{T}^{3}, … | I, G^{T}^{2}, …] [ \begin{matrix}T_{1}^{(0)} & W^{(0)} \\ Z^{(0)} & T_{2}^{(0)} \end{matrix} ] = [O, O, … | A_{0}^{T} , O, …]  (16)

where T_{1}^{(0)}, T_{2}^{(0)} and Z^{(0)} are infinite block triangular block Toeplitz matrices, W^{(0)} is an infinite block Hessenberg block Toeplitz matrix. Applying one step of block-Gaussian elimination to the 2 × 2 block system (16) yields

[I, G^{T}^{2}, G^{T}^{4}, …]H^{(1)} = [A_{0}^{T} , O, O, …]  (17)

where H^{(1)} = T_{2}^{(0)} − Z^{(0)}T_{1}^{(0)^{−1}}W^{(0)} is the Schur complement of T_{2}^{(0)}. It is interesting to observe that Q^{(1)} = I − H^{(1)} is a nonnegative block Hessenberg matrix which, except for the first block column, has the block Toeplitz structure of the matrix of (12).

In other words, H^{(1)} is uniquely determined by the blocks Â_{i}^{(1)}, i ≥ 1, A_{i}^{(1)}, i ≥ 0, defining the first two block columns of Q^{(1)}.

We may recursively apply the same reduction (cyclic reduction) to the matrix equation (17) thus obtaining a sequence of matrices {H^{(j)}} such that Q^{(j)} = I−H^{(j)} is a nonnegative matrix having, except for the first block column, the same block Hessenberg and Toeplitz structure of the matrix in (12). Each matrix Q^{(j)} is determined by the blocks Â_{i}^{(j)}, i ≥ 1, A_{i}^{(j)}, i ≥ 0, defining its first two block columns.

The block entries {A_{i}^{(j+1)}}_{i≥0}, {Â_{i}^{(j+1)}}_{i≥1}, obtained at step j + 1, are related to the block entries {A_{i}^{(j)}}_{i≥0}, {Â_{i}^{(j)}}_{i≥1}, obtained at step j, by means of functional relations, that express in compact form the block odd-even permutation and the Gaussian elimination of one step of cyclic reduction. More precisely, for any j ≥ 0, let us associate with the matrix sequences {A_{i}^{(j)}}_{i≥0}, {Â_{i}^{(j)}}_{i≥1}, the formal matrix power series A^{(j)}(z) = ∑_{i=0}^{+∞} A_{i}^{(j)}z^{i}, Â^{(j)}(z) = ∑_{i=0}^{+∞} Â_{i+1}^{(j)}z^{i}, respectively, where A_{i}^{(0)} = Â_{i}^{(0)} = A_{i}^{T}, i ≥ 1, A_{0}^{(0)} = A_{0}^{T}. Then the following functional relations hold:

\{ \begin{matrix}A^{(j+1)}(z) = zA_{odd}^{(j)}(z) + A_{even}^{(j)}(z)(I − A_{odd}^{(j)}(z))^{−1}A_{even}^{(j)}(z) \\ Â^{(j+1)}(z) = Â_{odd}^{(j)}(z) + A_{even}^{(j)}(z)(I − A_{odd}^{(j)}(z))^{−1}Â_{even}^{(j)}(z) \end{matrix}  (18)

where

\begin{matrix}A_{even}^{(j)}(z) = ∑_{i=0}^{+∞} A_{2i}^{(j)}z^{i}, & A_{odd}^{(j)}(z) = ∑_{i=0}^{+∞} A_{2i+1}^{(j)}z^{i} \\ Â_{even}^{(j)}(z) = ∑_{i=0}^{+∞} Â_{2(i+1)}^{(j)}z^{i}, & Â_{odd}^{(j)}(z) = ∑_{i=0}^{+∞} Â_{2i+1}^{(j)}z^{i}. \end{matrix}

The relations (18), besides expressing in compact form the recursions which the cyclic reduction method is based on, provide a basic tool for the efficient computation of the matrices {A_{i}^{(j)}}_{i≥0}, {Â_{i}^{(j)}}_{i≥1}, by using FFT’s.

The following results are fundamental for devising an efficient algorithm based on cyclic reduction:

Theorem 3 The blocks {A_{i}^{(j)}}_{i≥0}, {Â_{i}^{(j)}}_{i≥1} are nonnegative matrices such that ∑_{i=0}^{+∞} A_{i}^{(j)}, A_{0}^{T} + ∑_{i=1}^{+∞} Â_{i}^{(j)} are column stochastic. Moreover, if the matrix I − ∑_{i=0}^{+∞} G^{T}^{i·2^{j}} Â_{i+1}^{(j)} is

nonsingular, then G^{T} = A_{0}^{T} (I − ∑_{i=0}^{+∞} G^{T}^{i·2^{j}} Â_{i+1}^{(j)})^{−1} .

Under mild conditions, usually satisfied in the applications, the block entries A_{i}^{(j)} and Â_{i}^{(j)} quadratically converge to zero, for i ≥ 2, for j → ∞, as stated by the following theorem [6], [7], [5]:

Theorem 4 Let G′ = lim_{j→∞} G^{T}^{j}. Then the following convergence properties hold:

1. lim_{j→∞} Â_{i}^{(j)} = 0, for i ≥ 2;

2. if the entries of the matrix (I − ∑_{i=1}^{+∞} A_{i}^{(j)})^{−1} are bounded above by a constant, then the sequence of matrices R^{(j)} = A_{0}^{(j)} (I − ∑_{i=1}^{+∞} A_{i}^{(j)})^{−1} converges quadratically to the matrix G′;

3. if the solution G of (4) is irreducible then lim_{j→∞} A_{0}^{(j)}(I − A_{1}^{(j)})^{−1} = G′ and lim_{j→∞} A_{i}^{(j)} = 0, i ≥ 2.

The use of FFT techniques allows us to reduce the computational cost of the j-th step of the cyclic reduction algorithm to O(m^{3}n_{j} + m^{2}n_{j} log n_{j}) arithmetic operations, where n_{j} is such that the matrices A_{i}^{(j)} are negligible for i > n_{j}.

In the case of QBD problems, where A_{i} = O for i > 2, functional relations (18) are simplified, leading to an algorithm similar to that of Latouche and Ramaswami [18], but with a lower computational cost:

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_{0}^{(j)}(I − A_{1}^{(j)})^{−1}Â_{2}^{(j)}
Â_{2}^{(j+1)} = A_{2}^{(j+1)}.

For details on the implementation of the algorithm we refer the reader to papers [7] and [5], where two different FFT based techniques are proposed for the computation of the coefficients of the matrix power series A^{(j)}(z) and Â^{(j)}(z). A Fortran 90 code implementing the version described in [5] can be obtained by anonymous ftp at morse.dm.unipi.it in the directory /pub/MARKOV.

5 SOLVING NON-SKIP-FREE M/G/1 TYPE MARKOV CHAINS

The properties of displacement rank and the nice features of the cyclic reduction of Section 4.2 allow us to devise a new efficient quadratically convergent algorithm for computing the matrix G associated with reblocked non-skip-free M/G/1 type Markov chains.

Let us introduce the following notations: the block vectors E_{1} = [I, O, …, O]^{T}, E_{k} = [O, …, O, I]^{T} denote mk × m matrices made up by the first and the last m columns of the identity matrix of order mk.

Consider the matrices A_{i} of (3) obtained by reblocking the matrix (2) and consider the sequence of matrix power series A^{(j)}(z), Â^{(j)}(z) obtained by means of functional relations (18), where each A is replaced with A. We may easily observe that the formal matrix power series H^{(0)}(z) = ∑_{i=0}^{+∞} H_{i}^{(0)}z^{i} = zI −A^{(0)}(z) and Ĥ^{(0)}(z) = ∑_{i=0}^{+∞} Ĥ_{i}^{(0)}z^{i} = I −Â^{(0)}(z) are block Toeplitz matrices. At a first glance, the Toeplitz structure looks lost by the matrix power series H^{(j)}(z) = ∑_{i=0}^{+∞} H_{i}^{(j)}z^{i} = zI − A^{(j)}(z) and Ĥ^{(j)}(z) = ∑_{i=0}^{+∞} Ĥ_{i}^{(j)}z^{i}I − Â^{(j)}(z), for j ≥ 1. However, the block Toeplitz-like structure is kept at each step j of cyclic reduction, as stated by the following:

Theorem 5 For the matrix power series H^{(j)}(z) = zI − A^{(j)}(z), Ĥ^{(j)}(z) = I − Â^{(j)}(z) generated at j-th step of cyclic reduction we have

\begin{matrix}Δ(H^{(j)}(z)) = (H^{(j)}(z)E_{1})U^{(j)}(z) − V^{(j)}(z)(E_{k}^{T} H^{(j)}(z)) \\ Δ(Ĥ^{(j)}(z)) = −E_{1}(E_{k}^{T} A_{0}^{(0)}) − (H^{(j)}(z)E_{1})Û^{(j)}(z) − V^{(j)}(z)(E_{k}^{T} Ĥ^{(j)}(z)) \end{matrix}  (19)

where U^{(j)}(z) = ∑_{i=0}^{+∞} U_{i}^{(j)}z^{i}, Û^{(j)}(z) = ∑_{i=0}^{+∞} Û_{i}^{(j)}z^{i} are block row vector power series, V^{(j)}(z) = ∑_{i=0}^{+∞} V_{i}^{(j)}z^{i} is a block column vector power series, defined by

U^{(j)}(z^{2}) = −U^{(j−1)}(z)(H^{(j−1)}(z) − H^{(j−1)}(−z))^{−1}H^{(j−1)}(−z) +
U^{(j−1)}(−z)(H^{(j−1)}(z) − H^{(j−1)}(−z))^{−1}H^{(j−1)}(z)
Û^{(j)}(z^{2}) = (Û^{(j−1)}(z) − Û^{(j−1)}(−z))/(2z) + (U^{(j−1)}(z) − U^{(j−1)}(−z))
(H^{(j−1)}(z) − H^{(j−1)}(−z))^{−1}(Ĥ^{(j−1)}(z) − Ĥ^{(j−1)}(−z))/(2z),
V^{(j)}(z^{2}) = −H^{(j−1)}(−z)(H^{(j−1)}(z) − H^{(j−1)}(−z))^{−1}V^{(j−1)}(z) +
H^{(j−1)}(z)(H^{(j−1)}(z) − H^{(j−1)}(−z))^{−1}V^{(j−1)}(−z)

for j ≥ 1, and U^{(0)}(z) = zE_{k}^{T}, Û^{(0)}(z) = E_{k}^{T}, V^{(0)}(z) = zE_{1}.

Proof. The proof follows by induction on j. □

From the above theorem, the following representations hold for the matrix power series H^{(j)}(z) and Ĥ^{(j)}(z):

Theorem 6 At each step j of cyclic reduction it holds that

H^{(j)}(z) = L(H^{(j)}(z)E_{1})(I − L^{T} (Z_{k}U^{(j)}(z)^{T})) + L(V^{(j)}(z))L^{T} (Z_{k}H^{(j)}(z)^{T} E_{k})
Ĥ^{(j)}(z) = L(Ĥ^{(j)}(z)E_{1}) + L^{T} (Z_{k}A_{0}^{(0)T} E_{k}) +
L(H^{(j)}(z)E_{1})L^{T} (Z_{k}Û^{(j)}(z)^{T}) + L(V^{(j)}(z))L^{T} (Z_{k}Ĥ^{(j)}(z)^{T} E_{k}).

Proof. The proof simply follows from Theorem 5 and from the results of Section 2.2. □

From the above results it follows that, at each step j, only the seven block vectors power series U^{(j)}(z), Û^{(j)}(z), V^{(j)}(z), C^{(j)}(z) = H^{(j)}(z)E_{1}, R^{(j)}(z) = E_{k}^{T} H^{(j)}(z), Ĉ^{(j)}(z) = Ĥ^{(j)}(z)E_{1}, R̂^{(j)}(z) = E_{k}^{T} Ĥ^{(j)}(z), need to be computed in order to represent the matrix power series A^{(j)}(z) and Â^{(j)}(z).

Explicit relations between the block vector power series C^{(j)}(z) = ∑_{i=0}^{+∞} C_{i}^{(j)}z^{i}, R^{(j)}(z) = ∑_{i=0}^{+∞} R_{i}^{(j)}z^{i}, Ĉ^{(j)}(z) = ∑_{i=0}^{+∞} Ĉ_{i}^{(j)}z^{i}, R̂^{(j)}(z) = ∑_{i=0}^{+∞} R̂_{i}^{(j)}z^{i} at two subsequent steps of cyclic reduction can be obtained by multiplying functional relations (18) on the left by E_{k}^{T}, and on the right by E_{1}.

Moreover the block vectors defining A^{(j)}(z) and Â^{(j)}(z) are related to the block vectors obtained at step j − 1, which define A^{(j−1)}(z) and Â^{(j−1)}(z), by means of functional relations that only involve sums of products between block triangular block Toeplitz matrices. Hence, the computational cost needed to perform j-th step of cyclic reduction by using FFT’s reduces to O(n_{j}m^{3}k log^{2} k + m^{2}kn_{j} log n_{j}), where n_{j} is such that A_{i}^{(j)} and Â_{i}^{(j)} are negligible for i > n_{j}. We refer the reader to [4] for details.

In the case where A_{i} = O for i > 2 the above results are further simplified. We can give explicit expressions of the blocks H_{0}^{(j)} = −A_{0}^{(j)}, H_{1}^{(j)} = I − A_{1}^{(j)}, H_{2}^{(j)} = −A_{2}^{(j)}, Ĥ_{1}^{(j)} = I − Â_{1}^{(j)} by means of the following theorems:

Theorem 7 If A_{i} = O for i > 2, then for the matrices H_{0}^{(j)}, H_{1}^{(j)}, H_{2}^{(j)}, Ĥ_{1}^{(j)}, j ≥ 0, generated by the cyclic reduction (18) we have

Δ(H_{0}^{(j)}) = C_{0}^{(j)}U_{0}^{(j)} − V_{0}^{(j)}R_{0}^{(j)}
Δ(H_{1}^{(j)}) = C_{0}^{(j)}U_{1}^{(j)} + C_{1}^{(j)}U_{0}^{(j)} − V_{0}^{(j)}R_{1}^{(j)} − V_{1}^{(j)}R_{0}^{(j)}
Δ(H_{2}^{(j)}) = S^{(j)}U_{1}^{(j)} − V_{1}^{(j)}R̂_{1}^{(j)}
Δ(Ĥ_{1}^{(j)}) = −E_{1}R_{0}^{(0)} − C_{0}^{(j)}U_{1}^{(j)} − V_{0}^{(j)}R̂_{1}^{(j)},

where S^{(0)} = C_{1}^{(0)}, S^{(j)} = S^{(j−1)} − H_{2}^{(j−1)}H_{1}^{(j−1)}^{−1}C_{0}^{(j−1)}, for j ≥ 1.

Proof. The proof follows by comparing the terms of the same degree in formula (19). □

Theorem 8 At each step j of cyclic reduction, the matrices H_{0}^{(j)}, H_{1}^{(j)}, H_{2}^{(j)}, H_{1}^{(j)^{−1}}, Ĥ_{1}^{(j)} can be rewritten as

H_{0}^{(j)} = L(C_{0}^{(j)})(I − L^{T} (Z_{k}U_{0}^{(j)T})) + L(V_{0}^{(j)})L^{T} (Z_{k}R_{0}^{(j)T})
H_{1}^{(j)} = L(C_{1}^{(j)})(I − L^{T} (Z_{k}U_{0}^{(j)T})) − L(C_{0}^{(j)})L^{T} (Z_{k}U_{1}^{(j)T}) + L(V_{0}^{(j)})L^{T} (Z_{k}R_{1}^{(j)T}) +
L(V_{1}^{(j)})L^{T} (Z_{k}R_{0}^{(j)T})
H_{2}^{(j)} = L(C_{2}^{(j)}) − L(S^{(j)})L^{T} (Z_{k}U_{1}^{(j)T}) + L(V_{1}^{(j)})L^{T} (Z_{k}R̂_{1}^{(j)T})
H_{1}^{(j)^{−1}} = L(H_{1}^{(j)^{−1}}E_{1}) + L^{T} (Z_{k}H_{1}^{(j)^{−T}} U_{0}^{(j)T}) + L(H_{1}^{(j)^{−1}}C_{0}^{(j)})L^{T} (Z_{k}H_{1}^{(j)^{−T}} U_{1}^{(j)T}) −
L(H_{1}^{(j)^{−1}}V_{1}^{(j)})L^{T} (Z_{k}H_{1}^{(j)^{−T}} R_{0}^{(j)T})
Ĥ_{1}^{(j)} = L(Ĉ_{1}^{(j)}) + L^{T} (Z_{k}R_{0}^{(0)T}) − L(C_{0}^{(j)})L^{T} (Z_{k}U_{1}^{(j)T}) + L(V_{0}^{(j)})L^{T} (Z_{k}R̂_{1}^{(j)T}).

Proof. The proof follows from Theorem 7 and by the results of Section 2.2. □

Explicit relations for the vectors U_{i}^{(j)}, V_{i}^{(j)} C_{i}^{(j)}, R_{i}^{(j)}, Ĉ_{i}^{(j)}, R̂_{i}^{(j)}, can be provided:

U_{0}^{(j)} = U_{0}^{(j−1)} − U_{1}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{0}^{(j−1)}
U_{1}^{(j)} = −U_{1}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{2}^{(j−1)}
V_{0}^{(j)} = V_{0}^{(j−1)} − H_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}V_{1}^{(j−1)}
V_{1}^{(j)} = −H_{2}^{(j−1)}H_{1}^{(j−1)^{−1}}V_{1}^{(j−1)}
C_{0}^{(j)} = −H_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}C_{0}^{(j−1)}
C_{1}^{(j)} = C_{1}^{(j−1)} − H_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}C_{2}^{(j−1)} − H_{2}^{(j−1)}H_{1}^{(j−1)^{−1}}C_{0}^{(j−1)}
C_{2}^{(j)} = −H_{2}^{(j−1)}H_{1}^{(j−1)^{−1}}C_{2}^{(j−1)}
R_{0}^{(j)} = −R_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{0}^{(j−1)}
R_{1}^{(j)} = R_{1}^{(j−1)} − R_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{2}^{(j−1)} − R_{2}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{0}^{(j−1)}
R_{2}^{(j)} = −R_{2}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{2}^{(j−1)}
R̂_{1}^{(j)} = R̂_{1}^{(j−1)} − H_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}R_{2}^{(j−1)}
Ĉ_{1}^{(j)} = Ĉ_{1}^{(j−1)} − C_{0}^{(j−1)}H_{1}^{(j−1)^{−1}}H_{2}^{(j−1)}

for j ≥ 1, where U_{0}^{(0)} = V_{0}^{(0)} = 0, U_{1}^{(0)} = E_{k}^{T}, V_{0}^{(0)} = E_{1}. For this particular case, the computational cost per step is reduced to O(m^{3}k log^{2} k) operations.

For the blocks A_{i}^{(j)}, Â_{i}^{(j)} the quadratically convergence properties of the blocks A_{i}^{(j)}, Â_{i}^{(j)} of Section 4.2 still hold. Hence, due to the displacement structure of the blocks A_{i} of (3), we have devised a quadratically convergent algorithm, with low computational cost, for non-skip-free M/G/1 type Markov chains. This new algorithm particularly shows its effectiveness, with respect to the linearly convergent method of [13], when the block dimension k is large; in fact, the computational cost of one step of the latter method is O(mk^{2}) and many iterations may need to be computed in order to reach a good approximation of G.

Numerical experiments which show the efficiency of our approach are presented in next section.

6 NUMERICAL RESULTS

We have tested our algorithms for computing the matrix G on a problem arising from the mathematical modeling of a Metaring MAC Protocol [2]. For this particular example the blocks A_{i} have dimension m = 16. We have tested the algorithm for different values of the parameter ρ, 0 < ρ < 1, which represents the stability condition of the associated Markov chain [23].

The integer n_{0} such that the blocks A_{i} are negligible for i > n_{0} is equal to 168 for ρ = 0.1, is equal to 240 for ρ = 0.8, and is equal to 264 for the remaining tested values of ρ ≥ 0.9. We have compared cyclic reduction algorithm with the functional iteration method based on the recursion

\begin{matrix}X_{0} = I, \\ X_{n+1} = A_{0} (I − ∑_{i=2}^{+∞} A_{i}X_{n}^{i−1})^{−1} , n ≥ 0, \end{matrix}  (20)

which is the fastest among the classical linearly convergent functional iterations [22].

The numerical experiments have been performed on an alpha workstation with a base 2 arithmetic endowed with 53 bits.

Table 1 reports the CPU time (in seconds), the number of iterations needed by the Cyclic Reduction (CR) algorithm to compute an approximation of the matrix G; the residual of the computed approximation is also reported. We compared cyclic reduction

ρ Time (s.) Iterations Residual
0.1  0.9  9  1.8 · 10^{−13}
0.8  1.5  13  5.4 · 10^{−14}
0.9  2.3  14  1.5 · 10^{−14}
0.95  2.4  16  1.8 · 10^{−14}
0.96  2.4  17  2.3 · 10^{−14}
0.97  2.5  20  4.2 · 10^{−14}

Table 1: Cyclic Reduction

with classical techniques based on functional iterations: Table 2 reports the CPU time (in seconds) and the number of iterations needed by the Functional Iteration Formula (FIF) defined in (20), the residual of the computed approximation, and the ratio between the times needed by this algorithm and by CR. Finally we have considered the problem corresponding to ρ = 0.8 and we have compared the CPU time needed by the doubling algorithm based on displacement rank of Section 4.1 and the Latouche-Stewart algorithm [19]. Table 3 reports the number of CPU seconds required by our algorithm (algorithm DR: Displacement Rank), together with the residual error, for the values of the block dimension n of (14) ranging from 256 to 4096. It reports also the time required by the Latouche–Stewart algorithm [19] (denoted by LS) to reach the same residual error of our algorithm, and the ratio between the times. The symbol “*” denotes that the algorithm is not able to reach the required accuracy for lack of memory.

ρ Time (s.) Iterations Residual FIF/CR
0.1  0.3  22  2.2 · 10^{−14}  0.3
0.8  3.9  148  1.1 · 10^{−13}  2.6
0.9  10.8  373  1.2 · 10^{−13}  4.7
0.95  44.4  1534 1.3 · 10^{−13} 18.5
0.96  91.4  3158 1.3 · 10^{−13} 38.1
0.97  96.8  3336 1.3 · 10^{−13} 38.7

Table 2: Functional Iteration Method

DR  LS
n time residual  time residual ratio
256 8.7 2.5 · 10^{−3}  38.3 2.5 · 10^{−3} 4.4
512 18.1 4.0 · 10^{−4} 103.5 4.0 · 10^{−4} 5.7
1024 39.2 1.7 · 10^{−5} 264.9 1.7 · 10^{−5} 6.8
2048 89.3 3.4 · 10^{−8}  *  *  *
4096 193.4 1.3 · 10^{−13}  *  *  *

Table 3: Fast and customary doubling algorithm

We have also tested the cyclic reduction based algorithm for computing the matrix G associated with reblocked non-skip-free M/G/1 type Markov chains. We have randomly generated blocks A_{i} of dimension m = 4, such that A_{i} = 0 for i > 2k, for different values of k. That is, we have considered the simple case where A_{i} = 0 for i > 2. We have compared the above algorithm (NSF-CR Algorithm) with cyclic reduction algorithm for QBD problem of Section 4.2 (CR Algorithm), and with the algorithm [13] of Gail, Hantler, Taylor (GHT Algorithm). Table 4 reports the CPU time (in seconds) and the number of iterations needed by NSF-CR, CR and GHT algorithms to compute an approximation of the matrix G; the residual of the computed approximation is also reported. The symbol “*” denotes that the algorithm is not able to run for lack of memory.

NSF-CR  CR  GHT
k Steps Time Residual Steps Time Residual Steps Time Residual
16  10  14.4 5.9 · 10^{−16} 10  1.1 9.9 · 10^{−16} 1534  6.1 9.9 · 10^{−11}
32  9  19.2 6.3 · 10^{−15}  9  9.8 1.4 · 10^{−15} 848  13.3 9.9 · 10^{−11}
64  9  37.7 5.4 · 10^{−13}  9  136.3 2.4 · 10^{−14} 980  60.4 9.9 · 10^{−11}
128 10  89.5 1.2 · 10^{−11}  *  *  *  2229 555.5 9.9 · 10^{−11}
256  9  146.4 1.9 · 10^{−15}  *  *  *  800 816.4 9.9 · 10^{−11}
512 12 448.1 3.1 · 10^{−11}  *  *  *  1380 5830.3 1.7 · 10^{−4}

Table 4: Non-skip-free Markov chains

A MATRIX POLYNOMIAL ARITHMETIC

Throughout this section i is the imaginary unit such that i^{2} = −1, ω = cos \frac{2π}{n} + i sin \frac{2π}{n} denotes a principal n-th root of 1, x̄ denotes the complex conjugate of x ∈ C. The vector y = (y_{j}) ∈ C^{n} such that y_{j} = ∑_{i=0}^{n−1} x_{i}ω^{ij} is called the discrete Fourier transform of the vector x = (x_{j}) ∈ C^{n}, and is denoted by y = DFT_{n}(x); the vector x = IDFT_{n}(y), x = (x_{j}), x_{j} = \frac{1}{n} ∑_{i=0}^{n−1} y_{i}ω̄^{ij} denotes the inverse discrete Fourier transform of y.

Given the block vector X = (X_{0}, …, X_{n−1}), where X_{s} = (x_{i,j}^{(s)}) are m × m matrices, we define Y = (Y_{0}, …, Y_{n−1}) = DFT(X) the block vector such that Y_{s} = (y_{i,j}^{(s)}) are m × m matrices and (y_{i,j}^{(0)}, …, y_{i,j}^{(n−1)}) = DFT(x_{i,j}^{(0)}, …, x_{i,j}^{(n−1)}), i, j = 1, …, m. Analogously we set X = IDFT(Y) if (x_{i,j}^{(0)}, …, x_{i,j}^{(n−1)}) = IDFT(y_{i,j}^{(0)}, …, y_{i,j}^{(n−1)}). In this way the computation of DFT(X) and IDFT(Y) is reduced to the computation of m^{2} DFT’s and IDFT’s, respectively, for the cost of O(m^{2}n log n) ops [11].

Let P (z) = ∑_{i=0}^{n−1} P_{i}z^{i}, Q(z) = ∑_{i=0}^{n−1} Q_{i}z^{i}, R(z) = ∑_{i=0}^{n−1} R_{i}z^{i} be matrix polynomials such that R(z) = P (z)Q(z) mod (z^{n} − 1). Then the matrix coefficients of R(z) can be computed, given the matrix coefficients of P (z) and Q(z), by means of the following algorithm based on evaluation/interpolation at the roots of 1.

Algorithm 1. Computation of matrix polynomial product modulo z^{n} − 1.

Input. Positive integers M, m and the m × m matrices P_{0}, …, P_{n−1}, Q_{0}, …, Q_{n−1}, coefficients of the matrix polynomials P (z) = ∑_{i=0}^{n−1} P_{i}z^{i}, Q(z) = ∑_{i=0}^{n−1} Q_{i}z^{i}, where

n = 2^{M}.

Output. The m×m matrices R_{0}, …, R_{n−1}, coefficients of the matrix polynomial R(z) = ∑_{i=0}^{n−1} R_{i}z^{i} such that R(z) = P (z)Q(z) mod (z^{n} − 1).

Computation. 1. (Evaluation) Compute the entries of the 2n matrices U_{s} = P (ω^{s}), V_{s} = Q(ω^{s}), s = 0, …, n − 1, ω = cos\frac{2π}{n} + i sin\frac{2π}{n}, in the following way

\begin{matrix}(U_{0}, …, U_{n−1}) = DFT(P_{0}, …, P_{n−1}), \\ (V_{0}, …, V_{n−1}) = DFT(Q_{0}, …, Q_{n−1}). \end{matrix}

2. Compute the n matrix products W_{s} = U_{s}V_{s}, s = 0, …, n − 1, such that W_{s} = P (ω^{s})Q(ω^{s}) = R(ω^{s}).

3. (Interpolation) Compute the entries of the matrices R_{0}, …, R_{n−1} by means of the equation (R_{0}, …, R_{n−1}) = IDFT(W_{0}, …, W_{n−1}).

The cost of Algorithm 1 is O(m^{3}n + m^{2}n log n) ops.

Consider now the problem of inverting a polynomial, i.e., of solving the congruence P (z)Q(z) = I mod z^{n}.

Let us denote by T_{n} the matrix in (7), assume for simplicity m = n = 2^{M}, M positive integer, and partition T_{n} as follows

T_{n} = ( \begin{matrix}T_{n/2} & O \\ H_{n/2} & T_{n/2} \end{matrix} )

where all the four blocks are (n/2) × (n/2) block Toeplitz matrices. Then the matrix T_{n}^{−1} can be easily written as

T_{n}^{−1} = ( \begin{matrix}T_{n/2}^{−1} & O \\ −T_{n/2}^{−1}H_{n/2}T_{n/2}^{−1} & T_{n/2}^{−1} \end{matrix} )  (21)

and the first block column (Q_{0}^{T} , …, Q_{n−1}^{T})^{T} of T_{n}^{−1}, which solves (7), can be computed from the first block column (Q_{0}^{T} , …, Q_{n/2−1}^{T})^{T} of T_{n/2}^{−1} by means of two multiplications between an (n/2) × (n/2) block Toeplitz matrix and a block vector. In this way we arrive at the following algorithm which extends to the block case the algorithm of Lafon and Sieveking-Kung (compare [9]).

Algorithm 2. Solution of the congruence P (z)Q(z) = I mod z^{n}, or, equivalently, of the system (7).

Input. A positive integer M, the m × m matrices P_{0}, P_{1}, …, P_{n−1}, n = 2^{M}, defining the matrix polynomial P (z) = ∑_{i=0}^{n−1} P_{i}z^{i}, det P_{0} ≠ 0.

Output. The matrices Q_{0}, Q_{1}, …, Q_{n−1} satisfying (7), or equivalently, such that the polynomial Q(z) = ∑_{i=0}^{n−1} Q_{i}z^{i} solves the congruence P (z)Q(z) = I mod z^{n}.

Computation. 1. Compute Q_{0} = P_{0}^{−1}.

2. For i = 0, …, M − 1, given the first column U = (Q_{0}^{T} , …, Q_{2^{i}−1}^{T})^{T} of T_{2^{i}}^{−1}, compute the block vector V = (Q_{2^{i}}^{T}, …, Q_{2^{i+1}−1}^{T})^{T}, which defines the remaining blocks of the first column of T_{2^{i+1}}^{−1}, by applying equation (21) with n = 2^{i+1} in the following way: by applying Algorithm 1 the products W = H_{2^{i}}U and V = T_{2^{i}}^{−1}W.

Observe that at the i-th stage of Algorithm 2 the computation of 3 DFT’s and 2 IDFT’s of order 2^{i+1} must be performed, together with the multiplication of about 2^{i+1} complex m × m matrices. The overall cost is O(m^{3}n + m^{2}n log n) ops.

References

[1] N. Akar and K. Sohraby. An invariant subspace approach in M/G/1 and G/M/1 type Markov chains. Commun. Statist. Stochastic Models, 1997. To appear.

[2] G. Anastasi, L. Lenzini, and B. Meini. Performance evaluation of a worst case model of the Metaring MAC Protocol with global fairness. Performance Evaluation, 29:127– 151, 1997.

[3] D. Bini. On a class of matrices related to Toeplitz matrices. Technical Report 83-5, SUNY University, Albany, New York, 1983.

[4] D. Bini and B. Meini. Effective methods for solving banded Toeplitz systems. SIAM J. Matrix Anal. Appl., 1997. Submitted for publication.

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

[6] D. 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, 1996.

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

[8] D. Bini and V. Pan. Improved parallel computations with Toeplitz-like and Hankel matrices. Linear Algebra Appl., 188,189:3–29, 1993.

[9] D. Bini and V. Pan. Matrix and Polynomial Computations, Vol. 1: Fundamental Algorithms. Birkhäuser, Boston, 1994.

[10] B. L. Buzbee, G. H. Golub, and C. W. Nielson. On direct methods for solving Poisson’s equation. SIAM J. Num. Anal., 7:627–656, 1970.

[11] D.F. Elliott and K.R. Rao. Fast Transform Algorithms, Analyses, Applications. Academic Press, New York, 1982.

[12] H. R. Gail, S. L. Hantler, and B.A. Taylor. Solutions of the basic matrix equation for M/G/1 and G/M/1 type Markov Chains. Commun. Statist. Stochastic Models, 10:1–43, 1994.

[13] H. R. Gail, S. L. Hantler, and B.A. Taylor. Non-skip-free M/G/1 and G/M/1 type Markov Chains. Adv. Appl. Prob., 1997. To appear.

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

[15] T. Kailath, A. Viera, and M. Morf. Inverses of Toeplitz operators, innovations and orthogonal polynomials. SIAM Review, 20:106–119, 1978.

[16] G. Latouche. Algorithms for infinite markov chains with repeating columns. In C.D. Meyer and R.J. Plemmons, editors, Linear Algebra, Queueing Models and Markov Chains, pages 231–265. Springer-Verlag, New York,, 1993.

[17] G. Latouche. Newton’s iteration for non-linear equations in Markov chains. IMA J. of Numerical Analysis, 14:583–598, 1994.

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

[19] 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, 1996.

[20] E. Linzer. On the stability of transform-based circular deconvolution. SIAM J. Num. Anal., 29:1482–1492, 1992.

[21] B. Meini. An improved FFT-based version of Ramaswami’s formula. Commun. Statist. Stochastic Models, 13:223-238, 1997.

[22] B. Meini. New convergence results on functional iteration techniques for the numerical solution of M/G/1 type Markov chains. Numer. Math., 1997. To appear.

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

[24] V. Ramaswami. Nonlinear matrix equations in applied probability - Solution techniques and open problems. SIAM Review, 30:256–263, 1988.

[25] V. Ramaswami. A stable recursion for the steady state vector in Markov chains of M/G/1 type. Commun. Statist. Stochastic Models, 4:183–188, 1988.

[26] G. W. Stewart. On the solution of block Hessenberg systems. Numer. Lin. Algebra with Appl., 2:287–296, 1995.