Analyzing M/G/1 paradigms through QBDs: the role of the block structure

in computing the matrix G

Dario A. Bini^{∗} Beatrice Meini^{∗}  V. Ramaswami^{†}

Abstract

Recently, it has been proved by V. Ramaswami that an M/G/1 type Markov chain can be rewritten in terms of a suitable QBD process, with infinite and strongly structured block entries [11]. This reformulation allows to derive new algorithms for the computation of the matrix G. We show that by applying the cyclic reduction algorithm to the QBD process we generate a sequence of QBDs, having constant displacement rank. More precisely, we show that these blocks belong to the algebra generated by a suitable block Frobenius matrix. This fact allows us to devise a new quadratically convergent algorithm, that makes use of Fast Fourier Transforms and of Toeplitz matrix computations, for approximating the matrix G. A complexity analysis is performed.

1 Introduction

In a recent paper [11] it has been proved by V.Ramaswami that an M/G/1 type Markov chain can be embedded into a suitable QBD process with an infinite number of phases. This fact allows one to apply, for the solution of an M/G/1 type Markov chain, any algorithm designed for solving QBD problems. Indeed, this fact provides new algorithmic possibilities and new approaches for solving general M/G/1 Markov chains and to reduce the solution of the matrix equation

G = ∑_{i=0}^{+∞} D_{i}G^{i}

∗Dipartimento di Matematica, Università di Pisa, 56127 Pisa, Italy.

†AT&T Labs, Holmdel, NJ, USA

associated to an M/G/1 type Markov chain, to the solution of

H = C + AH + BH^{2}  (1)

associated with the embedding QBD.

An important feature of this reduction is that the special QBD problem derived from the general one, has a very strong structure; therefore, despite the new infinite size of the blocks, one might expect a reduction of the complexity if the underlying structure is suitably exploited when applying the QBD solver.

In this paper we provide a structure analysis of the QBD problem obtained by this reduction and show how to exploit this structure in order to speed up the logarithmic reduction algorithm (also known as cyclic reduction and hereafter referred to as CR) introduced in [8] and modified in [4, 3].

In particular, we prove that the sequence of QBD’s generated by cyclic reduction are defined by blocks that are related to the matrix algebra defined by a suitable block Frobenius matrix. This fact allows us to prove that the displacement rank of each block is finite, constant and, consequently, that a single step of cyclic reduction can be performed with a low computational costs by means of a small number of FFT’s.

More precisely, in order to solve the matrix equation (1), the arithmetic cost of one step of CR is O(m^{3}N + m^{2}N log N) arithmetic operations (ops), where m is the size of D_{i} (number of phases), N is the numerical degree of the matrix power series ∑_{i=0}^{+∞} D_{i}z^{i}, that is, N is such that ∑_{i=0}^{N} D_{i}e < (1−ϵ)e where e = (1, …, 1)^{T} and ϵ is the machine precision used in the computation.

The algorithm obtained in this way has been implemented in fortran 90 and the numerical experiments, performed so far, show very good convergence properties even for problems having a certain complexity [1]. Situations where the new approach can be more favorable with respect to directly applying CR to M/G/1 type Markov chains are investigated.

Similarly to the recent applications of CR to factoring polynomials and power series [2], the new approach seems to provide new and more efficient tools for computing factors of small degree of a given polynomial or power series.

The paper is organized as follows. In section 2 we recall the reduction of [11] from general QBD to M/G/1 problems. In section 3 we synthesize the algorithm of cyclic reduction with its convergence and basic structural properties. In section 4 we prove the main displacement structures of the blocks generated by cyclic reduction that are used in section 5 for the explicit description of an algorithm for the computation of the matrix G. In section 6 we briefly comment on some numerical experiments made with a prototype implementation of our algorithm and on possible applications of the newly presented approach.

2 M/G/1 paradigms through QBDs: structure analysis

Let

P = [ \begin{matrix}E_{0} & E_{1} & E_{2} & … \\ D_{0} & D_{1} & D_{2} & … \\ 0 & D_{0} & D_{1} & ... \\ ⋮ & ... & ... & ... \end{matrix} ]  (2)

be the transition matrix of an M/G/1 type Markov chain, where E_{i}, D_{i}, i ≥ 0, are m × m matrices, and let G be the minimal nonnegative solution of G = ∑_{i=0}^{∞} D_{i}G^{i}. Consider the QBD chain obtained by embedding the M/G/1 type Markov chain, as in [11]:

Q = [ \begin{matrix}Ã & B \\ C & A & B \\ & C & A & ... \\ & & ... & ... \end{matrix} ]  (3)

where

\begin{matrix}Ã = [ \begin{matrix}E_{0} & E_{1} & E_{2} & … \\ 0 & 0 & 0 & … \\ ⋮ & ⋮ & ⋮ \end{matrix} ] , A = [ \begin{matrix}D_{1} & D_{2} & D_{3} & … \\ 0 & 0 & 0 & … \\ ⋮ & ⋮ & ⋮ \end{matrix} ] , \\ C = [ \begin{matrix}D_{0} & 0 & 0 & … \\ 0 & 0 & 0 & … \\ ⋮ & ⋮ & ⋮ \end{matrix} ] , B = [ \begin{matrix}0 & 0 & 0 & … \\ I & 0 & 0 & … \\ 0 & I & 0 & … \\ ⋮ & ... & ... & ... \end{matrix} ] . \end{matrix}  (4)

We recall the following result of [11] that allows us to reduce a general M/G/1-type Markov chain to a QBD process with infinite phases.

Theorem 1 Let G be the minimal nonnegative solution of G = ∑_{i=0}^{+∞} D_{i}G^{i}, and denote H the minimal nonnegative solution of the matrix equation H = C + AH + BH^{2}, where A, B, C are defined in (4). Then it holds

H = [ \begin{matrix}G & 0 & 0 & … \\ G^{2} & 0 & 0 & … \\ G^{3} & 0 & 0 & … \\ ⋮ & ⋮ & ⋮ \end{matrix} ] .

Let us denote A_{N}, B_{N}, C_{N} and H_{N} the block N × N matrices obtained by truncating A, B, C and H at the block size N. Then it is easy to verify that, if D_{i} = 0, for i > N then H_{N} = C_{N} + A_{N}H_{N} + B_{N}H_{N}^{2}.

Besides its theoretical interest, the above results allow us to apply algorithms designed for the QBD case in order to solve general M/G/1-type Markov chains. In particular we may consider quadratically convergent methods like the logarithmic reduction technique of [8] or the slightly modified version of [3] based on the cyclic reduction algorithm.

The main issue in the application of these methods is the exploitation of the structures of the blocks (4) and of the solution H that must be taken into account if we want to arrive at an efficient design of the solution algorithm.

Since ∑_{i=0}^{+∞} D_{i} is stochastic, it holds that D_{i} is negligible for i large enough. In practice, if the entries of D_{i} are less than the machine precision of the floating point arithmetic used in the computation for i > N, then we may assume that D_{i} = 0 for i > N. In this way without loss of generality, for practical applications, we may assume that the blocks A, B, C of (4) have finite size N. Here and hereafter we make this assumption.

3 Properties of cyclic reduction

In this section we briefly recall the cyclic reduction technique [3] for the solution of the matrix equation H = B +AH +CH^{2} and analyze its behavior in the case (4).

The cyclic reduction method generates the sequences of matrices A_{n}, B_{n}, C_{n}, Â_{n}, n ≥ 0, such that

\begin{matrix}Â_{0} = A_{0} = A, B_{0} = B, C_{0} = C, \\ A_{n+1} = A_{n} + B_{n}(I − A_{n})^{−1}C_{n} + C_{n}(I − A_{n})^{−1}B_{n}, \\ B_{n+1} = B_{n}(I − A_{n})^{−1}B_{n}, \\ C_{n+1} = C_{n}(I − A_{n})^{−1}C_{n}, \\ Â_{n+1} = Â_{n} + C_{n}(I − A_{n})^{−1}B_{n}, n ≥ 0. \end{matrix}  (5)

It holds

H = C + Â_{n}H + B_{n}H^{2^{n}+1}.  (6)

Moreover, we have the following convergence result [9, 5, 4, 3].

Theorem 2 If the QBD process (3) is positive recurrent then the polynomial f(z) = det(C +zA+z^{2}B −zI) has m zeros in the closed unit disk. Moreover, if the polynomial f(z) has m zeros in the closed unit disk then ||B_{n}|| = O(λ^{2^{n}}), where || · || is any induced matrix norm, λ = 1/μ < 1 and μ = min{|z| : f(z) = 0, |z| > 1}.

From this property and from equation (6) we find that (I − Â_{n})^{−1}C provides accurate approximations of H for n large enough.

3.1 Structure of the blocks

The aim of this section is to describe the structural properties of the blocks A_{n}, B_{n}, C_{n} and Â_{n} generated by the cyclic reduction (5). We recall that we are assuming that A, B, C are defined by (4) and have finite size N.

Let us denote by e_{1} the Nm × m matrix made up by the first m columns of the identity matrix of size Nm and by Z the block down-shift matrix, i.e., the N × N finite section of the matrix B in (4). Here and hereafter we call block vector any Nm × m matrix.

Theorem 3 The matrices A_{n}, B_{n}, C_{n} and Â_{n} satisfy the following relations

\begin{matrix}Â_{n} = e_{1}v_{0}^{T} + Zu_{n}e_{1}^{T} , \\ A_{n} = e_{1}v_{n}^{T} + Zu_{n}e_{1}^{T} , \\ B_{n+1} = B_{n}(I + û_{n}S_{n}v_{n}^{T})B_{n}, \\ C_{n} = e_{1}D_{0}^{(n)}e_{1}^{T} , \\ (I − A_{n})^{−1} = I + [ e_{1} Zu_{n} ] [ \begin{matrix}S_{n} & T_{n} \\ S_{n} & I + T_{n} \end{matrix} ] [ \begin{matrix}v & T n \\ e & T 1 \end{matrix} ] , n ≥ 0, \end{matrix}  (7)

where

\begin{matrix}B_{0} = Z, D_{0}^{(0)} = D_{0}, û_{0} = e_{1}, v_{0}^{T} = [D_{1}, …, D_{N}] , \\ D_{0}^{(n+1)} = D_{0}^{(n)}(S_{n}v_{n,1} + T_{n} + I)D_{0}^{(n)}, \\ û_{n+1} = û_{n} + B_{n}û_{n}(S_{n}v_{n,1} + T_{n} + I)D_{0}^{(n)}, \\ v_{n+1}^{T} = v_{n}^{T} + D_{0}^{(n)}S_{n}v_{n}^{T}B_{n}, n ≥ 0, \end{matrix}  (8)

and

S_{n} = (I − v_{n}^{T} û_{n})^{−1}, T_{n} = S_{n}v_{n}^{T}Zu_{n}, u_{n} = Z^{T} û_{n},

and v_{n,1} is the first block entry of the block vector v_{n}^{T}.

Proof. Proceed by induction on n. For n = 0 the theorem trivially holds. In the general case, by applying (5) we immediately obtain the representation of B_{n}, C_{n} and A_{n}. The representation of (I−A_{n})^{−1} follows from the Sherman-Morrison-Woodbury formula [6] applied to (I − A_{n}).

From the above theorem it follows that the blocks Â_{n}, A_{n} and C_{n} still keep the sparsity properties of Â_{0}, A_{0} and C_{0}, respectively. More precisely, the matrices Â_{n} and A_{n} have null block entries except for the blocks in the first block column and in the first block row; and the matrices C_{n} have null block entries except for the one in position (1,1). At a first glance, the matrices B_{n} have no sparsity nor structural properties; indeed, it is easy to verify that if n < log_{2} N, then B_{n} has the following structure

B_{n} = [ \begin{matrix}0 & 0 & … & 0 \\ ∗ & ∗ & … & ∗ \\ ⋮ & ⋮ & & ⋮ \\ ∗ & ∗ & … & ∗ \\ I \\ . & . & . \\ & & I \end{matrix} ]

where ∗ denotes a possibly nonzero block entry and where the number of nonzero block rows is 2^{n} − 1; if n ≥ log_{2} N, then B_{n} seems to be a general N ×N matrix without any structural properties. The lack of structure of B_{n} would lead to an inefficient implementation of the cyclic reduction algorithm, in fact each step of CR would cost O(N^{3}m^{3}) ops.

Fortunately, it is possible to prove that the matrix B_{n} has a very strong structure (the displacement structure) that allows us to implement one step of CR with O(m^{2}N log N +m^{3}N) ops. This property is not trivial and needs some sophisticated technique to be proved. We will do this in section 4.

3.2 Convergence results

Before giving the convergence properties of the cyclic reduction algorithm, we recall that the N × N block Frobenius matrix, associated with the monic matrix polynomial P (z) = ∑_{i=0}^{N−1} P_{i}z^{i} + z^{N}I, is the matrix

F = [ \begin{matrix}−P_{N−1} & −P_{N−2} & … & −P_{0} \\ I & 0 & … & 0 \\ & ... & ... & ⋮ \\ ◯ & & I & 0 \end{matrix} ];

moreover, it holds det P (z) = det(zI − F) (compare [7]).

Concerning the convergence results we have the following

Theorem 4 Let f(z) = det(C +zA+z^{2}B −zI), g(z) = det(∑_{i=0}^{N} D_{i}z^{i} −zI). Then we have f(z) = z^{m(N−1)}g(z).

Proof. Assume that D_{0} is nonsingular. It holds

f(z) = det ([ \begin{matrix}D_{0} & & & & & ◯ \\ & −zI \\ & & . & . & . \\ ◯ & & & & & −zI \end{matrix} ] (I − zF))

where F is the block Frobenius matrix defined by the blocks −D_{0}^{−1}(D_{1} − I), −D_{0}^{−1}D_{i}, i = 2, …, N. Therefore, f(z) = z^{m(N−1)} det D_{0} det(z(z^{−1}I − F)). In the light of the properties of the determinants of block Frobenius matrices we have f(z) = z^{m(N−1)} det(zI − ∑_{i=0}^{N} D_{i}z^{i}) = z^{m(N−1)}g(z). If the matrix D_{0} is singular we may apply a continuity argument in order to complete the proof.

The result expressed in the above theorem is fundamental for proving the convergence of cyclic reduction. In fact, if the M/G/1 Markov chain (2) is positive recurrent then the polynomial g(z) has m roots inside the unit disk [10] (compare also theorem 2), therefore the polynomial f(z) has mN roots of modulus at most 1 of which m(N − 1) are zero. In the light of theorem 2 the latter property implies the following convergence result:

Theorem 5 If the M/G/1 Markov chain (2) is positive recurrent then ||B_{n}|| = O(λ^{2^{n}}), where || · || is any induced matrix norm, λ = 1/μ and μ = min{|z| : g(z) = 0, |z| > 1}. Moreover, the sequences {A_{n}}, {C_{n}}, {Â_{n}} are convergent.

From the above theorem we deduce the following property

Theorem 6 If the matrix G is nonsingular then u_{∗,i} = G^{i}, i = 1, …, N − 1, where u_{∗,i} = lim_{n} u_{n,i} and u_{n,i} is the i-th block component of the vector u_{n}.

Proof. From (6) and from theorem 5 we obtain (I − Â_{∗})H = B where Â_{∗} = lim_{n} Â_{n}. Since Â_{∗} has the same structure of Â_{n} we deduce that u_{∗,i}G = G^{i+1}, i = 1, …, N − 1, whence, from the nonsingularity of G we obtain that u_{∗,i} = G^{i}.

Observe that if the matrix G is nonsingular, then it can be approximated by the matrix u_{n,1}, for sufficiently large n; on the other hand, if G is singular, then it can be approximated by the matrix D_{0}(I −D_{1}−∑_{i=2}^{N} D_{i}u_{n,i−1})^{−1} that can be easily derived from the first block column of (6).

4 Displacement structure

We recall the following properties concerning the displacement rank of a matrix. Let Δ(A) = AZ − ZA and consider the following decomposition

Δ(A) = ∑_{i=1}^{k} r_{i}s_{i}^{T}  (9)

where r_{i} and s_{i} are block vectors. We define the block displacement rank of A as the minimum integer k such that (9) holds. If (9) holds then (compare [5])

A = L(Ae_{1}) − ∑_{i=1}^{k} L(r_{i})L(Zs_{i})^{T}

where L(w) denotes the block lower triangular block Toeplitz matrix defined by the block column w having block components w_{i}, i = 1, …, N,

L(w) = [ \begin{matrix}w_{1} & & & O \\ w_{2} & ... \\ ⋮ & ... & ... \\ w_{N} & … & w_{2} & w_{1} \end{matrix} ] .

We recall that a block Toeplitz matrix A is a matrix having block entries (A)_{i,j} = A_{i−j}. It is immediate to check that a block Toeplitz matrix has block displacement rank at most 2 and a block lower triangular block Toeplitz matrix has block displacement rank 0.

We recall also that the product L(w)v, for a pair of block vectors w and v can be computed with O(m^{2}N log_{2} N +m^{3}N) ops by using the Fast Fourier Transforms (FFT) [5].

4.1 Frobenius matrices and displacement rank

Let F be an N × N block Frobenius matrix. Then the powers of F can be expressed in terms of block Toeplitz matrices by means of the following result:

Theorem 7 For any integer k ≥ 0 it holds

Δ(F^{k}) = e_{1}u^{T} F^{k} − F^{k}e_{1}u^{T}  (10)

where u^{T} = e_{1}^{T} F is the first block row of F.

Proof. Since Z = F − e_{1}u^{T} it holds Δ(F^{k}) = F^{k}Z − ZF^{k} = F^{k}(F − e_{1}u^{T}) − (F − e_{1}u^{T})F^{k} = e_{1}u^{T} F^{k} − F^{k}e_{1}u^{T}.

From the above theorem it follows that any power of a block Frobenius matrix has block displacement rank at most 2.

4.2 Displacement properties of the blocks

Now, for the sake of simplicity we assume that det D_{0} ≠ 0. We recall a suitable relation that expresses the blocks A_{n}, B_{n}, C_{n} in functional form. By using this relation we may prove the displacement structure of the blocks B_{n} at each step of the cyclic reduction.

Theorem 8 Let ϕ_{n}(z) = zI − B_{n} − zA_{n} − z^{2}C_{n}. If D_{0} is nonsingular, then at each step n of cyclic reduction the matrix e_{1}^{T} F^{2^{n}−1}e_{1} is nonsingular and it holds:

ϕ_{n}(z) = (I − w_{n}(z)H_{n}e_{1}^{T})(zI − F^{2^{n}})  (11)

where H_{n} = (e_{1}^{T} F^{2^{n}−1}e_{1})^{−1}, w_{n}(z) = F^{2^{n}−1}e_{1} + ze_{1}D_{0} and

F = [ \begin{matrix}−D_{0}^{−1}(D_{1} − I) & −D_{0}^{−1}D_{2} & … & −D_{0}^{−1}D_{N} \\ I & 0 & … & 0 \\ & ... & ... & ⋮ \\ ◯ & & I & 0 \end{matrix} ] .

In particular, for the block B_{n} we have

B_{n} = ϕ_{n}(0) = − (F^{2^{n}−1} − F^{2^{n}−1}e_{1}H_{n}e_{1}^{T} F^{2^{n}−1}) Z.  (12)

Proof. We prove, by induction on n, that ψ_{n}(z) = ϕ_{n}(z)^{−1} verifies

(zI − F^{2^{n}})ψ_{n}(z) = (I − v_{n}(z)e_{1}^{T} )  (13)

for all the values z for which ψ_{n}(z) exists, where v_{n}(z) = e_{1}+z^{−1}F^{2^{n}−1}e_{1}D_{0}^{−1}. For n = 0 the equation trivially holds. Assuming the equation true for n, evaluate (13) for the values z and −z, compute the sums and the differences of both members of the equations obtained this way. Recover ψ_{n}(z) + ψ_{n}(−z) from one equation and substitute this value in the other one, and recover (z^{2}I − F^{2^{n+1}})(ψ_{n}(z) − ψ_{n}(−z))/(2z). The induction step is completed since ψ_{n+1}(z^{2}) = (ψ_{n}(z) − ψ_{n}(−z))/(2z) (compare [5]). Applying the Sherman-Woodbury-Morrison formula to the matrix (I −v_{n}(z)e_{1}^{T}) we obtain that this matrix is nonsingular if and only if z^{−1}e_{1}^{T} F^{2^{n}−1}e_{1} is nonsingular. Since ϕ_{n}(z) is singular for a finite number of values of z, the matrix ψ_{n}(z) = ϕ_{n}(z)^{−1} exists and is nonsingular for an infinite number of values of z. Therefore from (13) it follows that there exists at least a value z_{0} ≠ 0 of z such that the matrix I − v_{n}(z_{0})e_{1}^{T} is nonsingular. Whence we obtain that e_{1}^{T} F^{2^{n}−1}e_{1} is nonsingular and it holds that (I − v_{n}(z)e_{1}^{T})^{−1} = I − w_{n}(z)H_{n}e_{1}^{T} for any z ≠ 0. Therefore from (13) we obtain ϕ_{n}(z) = (I − w_{n}(z)H_{n}e_{1}^{T} ) (zI − F^{2^{n}}).

By applying a continuity argument we deduce that the above relation holds for any z. In particular for z = 0 we obtain (12).

An immediate consequence of the above theorem is that B_{n} has displacement rank at most 2 whenever D_{0} is nonsingular.

Theorem 9 If D_{0} is nonsingular, the matrices B_{n}, n ≥ 0, generated by (5), are such that

Δ(B_{n}) = ((e_{1} − F^{2^{n}−1}e_{1}H_{n})u^{T} + ZF^{2^{n}−1}e_{1}H_{n}e_{1}^{T} ) F^{2^{n}−1}Z (14)

where H_{n} = (e_{1}^{T} F^{2^{n}−1}e_{1})^{−1}, thus they have block displacement rank at most 2.

Proof. For the sake of notational simplicity set k = 2^{n} − 1. We have

\begin{matrix}Δ(B_{n}) = Δ(F^{k} − F^{k}e_{1}H_{n}e_{1}^{T} F^{k})Z = \\ Δ(F^{k})Z − (F^{k}e_{1}H_{n}e_{1}^{T} F^{k}Z − ZF^{k}e_{1}H_{n}e_{1}^{T} F^{k})Z. \end{matrix}  (15)

Since Z = F − e_{1}u^{T} we obtain F^{k}e_{1}H_{n}e_{1}^{T} F^{k}Z = F^{k}e_{1}H_{n}u^{T} F^{k} − F^{k}e_{1}u^{T}. Plugging the latter equation into (15) we arrive at (14).

If det D_{0} = 0, since B_{n} is well defined through (5), by using a continuity argument we arrive at the following:

Corollary 9.1 The matrices B_{n} have block displacement rank at most 2 for any n ≥ 0.

Observe that the matrix B_{n} can be represented by means of five block vectors: the first block column b_{n} of B_{n} and the 4 block vectors, c_{1}^{(n)}, c_{2}^{(n)}, r_{1}^{(n)}, r_{2}^{(n)}, such that

Δ(B_{n}) = c_{1}^{(n)}r_{1}^{(n)T} + c_{2}^{(n)}r_{2}^{(n)T} ,

defining the block displacement rank of B_{n}. In the next section we describe an algorithm for their computation.

5 Computing the matrix G

Now we are ready to describe the main algorithm for performing one step of cyclic reduction.

Algorithm: Computing one step of cyclic reduction.

Input:

– The block vectors u, v, such that A = e_{1}v^{T} + Zue_{1}^{T}, Â = e_{1}v_{0}^{T} + Zue_{1}^{T};

– the block entry D such that C = e_{1}De_{1}^{T};

– the block vectors b = Be_{1}, and c_{1}, c_{2}, r_{1}, r_{2}, such that Δ(B) = c_{1}r_{1}^{T} + c_{2}r_{2}^{T}.

Output:

– The block vectors u′, v′, such that A′ = e_{1}v^{′T} + Zu′e_{1}^{T}, Â′ = e_{1}v_{0}^{T} + Zu′e_{1}^{T};

– the block entry D′ such that C′ = e_{1}D′e_{1}^{T};

– the block vectors b′ = B′e_{1}, and c′_{1}, c′_{2}, r′_{1}, r′_{2}, such that Δ(B′) = c′_{1}r_{1}^{′T} + c′_{2}r_{2}^{′T},

where A′, B′, C′, Â′ are the matrices obtained after one step of cyclic reduction applied to A, B, C, Â.

Computation:

1. Compute S = (I − v^{T} û)^{−1}, T = Sv^{T} Zu, that define (I − A)^{−1} by means of formula (7).

2. Compute the block entry D′ by means of (8).

3. Compute the vectors û′, v′, by means of (8), where the matrix B is written in the form

B = L(b) − L(r_{1})L(Zc_{1})^{T} − L(r_{2})L(Zc_{2})^{T} .  (16)

4. Compute the vector b′ = B′e_{1} = −B(I −A)^{−1}b, where the matrix B is written in the form (16).

5. Compute the (Nm) × (7m) matrices W and Y such that

W Y^{T} = −Δ(B)(I −A)^{−1}B −BΔ((I −A)^{−1})B −B(I −A)^{−1}Δ(B).

6. Compute a QR factorization of W and Y, i.e., W = Q_{1}R_{1}, Y = Q_{2}R_{2}, where Q_{1}, Q_{2} are (Nm) × (7m) matrices, and R_{1}, R_{2} are (7m) × (7m) matrices; in this way W Y^{T} = Q_{1}(R_{1}R_{2}^{T})Q_{2}^{T}.

7. Compute the singular value decomposition of R_{1}R_{2}^{T}, i.e., R_{1}R_{2}^{T} = UΣV^{T}, where Σ = diag(σ_{1}, …, σ_{7m}), σ_{1} ≥ ⋯ ≥ σ_{7m}.

8. Set

[ c′_{1} c′_{2} ] = Q_{1}U_{2m}diag(σ_{1}, …, σ_{2m}), [ \begin{matrix}r & {′T \atop 1} \\ r & {′T \atop 2} \end{matrix} ] = V_{2m}^{T}Q_{2}^{T} ,

where U_{2m}, V_{2m} are the (Nm) × (2m) matrices made up by the first 2m columns of U and V, respectively.

Observe that the matrices W and Y^{T} have size (Nm) × (7m) since Δ(B) has rank at most 2m and Δ((I − A)^{−1}) has rank at most 3m.

The computational cost of stage 1 is O(m^{3}N) ops and the cost of stage 2 is O(m^{3}) ops; at stage 3 we have to perform products between matrices with block displacement rank at most 2, and block vectors, that can be computed in O(m^{3}N + m^{2}N log N) ops; similar computations are involved in stages 4 and 5, with the same cost; at stage 6 the computation of the QR factorizations costs O(m^{3}N) ops, and the computation of the singular values decomposition at stage 7 costs O(m^{3}) ops. Thus, the overall computational cost of one step of cyclic reduction is O(m^{3}N + m^{2}N log N).

6 Numerical experiments.

The algorithm has been implemented, in a prototypal version, in fortran 90. The program has been run on a pentium II server with the Linux system for different input data. The experiments have been performed for testing the convergence properties of the technique and its numerical stability. In this preliminary version, the Toeplitz structure of the blocks B_{n} has not been exploited for the fast computation of the matrix-vector products. Therefore in the table below we report only the number of steps and the residual error of the computation. The problems that we have considered are the ones of [1] with different values of the parameter ρ. For ρ = 1 the chain is not positive recurrent. The size m of the blocks is 16, the value N depends on the parameter ρ. The results have been compared with the experiments reported in [4] (PWCR algorithm).

We observe that the number of steps and the residual errors of the method described in this paper are comparable with the correspondent values of the algorithm of [4].

By looking at the complexity bounds we may determine the cases where the new approach is more convenient than the point-wise cyclic reduction. In fact, the cost of the point-wise cyclic reduction of [4] depends on the numerical degree of the power series generated at the intermediate steps. At the first step this number is N and in the limit is 1, but in the intermediate

New Algorithm  PWCR
ρ N Steps Residual Steps Residual
0.1 256 10 6.8 · 10^{−15} 9 1.8 · 10^{−13}
0.8 256 13 1.2 · 10^{−14} 13 5.4 · 10^{−14}
0.9 512 14 3.6 · 10^{−14} 14 1.6 · 10^{−14}
0.97 512 20 3.2 · 10^{−14} 20 4.2 · 10^{−14}

Table 1: Number of steps and residual errors

steps this number could be larger than N. On the other hand for the new method the cost is a function of only the numerical degree of the input matrix power series. Therefore the new approach should be more convenient in the cases where the numerical degrees of the intermediate power series generated by the point-wise cyclic reduction exceed N.

Another interesting application of the new approach concerns the problem of approximating a factor of a polynomial or of a power series. In [2] CR has been successfully applied to split a polynomial of degree n into two factors of degrees n − k and k. The most expensive stage consists in inverting a Toeplitz like matrix of size max{n − k, k}. With the new approach we have to invert a matrix of size min{n − k, k}. This might be a great advantage if k << n/2.

References

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

[2] D. Bini, L. Gemignani and B. Meini. Factorization of analytic functions by means of Koenig’s theorem and Toeplitz computations. Technical Report, Dipartimento di Informatica, Università di Pisa, 1998.

[3] D. 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. Bini and B. Meini. Improved cyclic reduction for solving queueing problems. Numerical Algorithms, 15:57–74, 1997.

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

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

[7] P. Lancaster and M. Tismenetski. Matrix Polynomials. Academic Press Inc., New York, 1982.

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

[9] M. F. Neuts. Matrix-Geometric Solutions in Stochastic Models. Johns Hopkins University Press, Baltimore, 1981.

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

[11] V. Ramaswami. The generality of Quasi Birth-and-Death processes. In A. Alfa and S. Chakravarthy, editors, Advances in Matrix Analytic Methods for Stochastic Models - Proceedings of the 2nd international conference on matrix analytic methods, pages 93–113. Notable Publications Inc, NJ, 1998.