An improved FFT-based version of Ramaswami’s formula

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

Abstract

Ramaswami’s formula allows a numerically stable computation of the probability invariant vector of a stochastic matrix of M/G/1 type. We observe that this formula consists in solving, by means of forward substitution, a block lower triangular block Toeplitz infinite system. By using the natural isomorphism between matrix polynomials and block Toeplitz matrices, we derive a new fast algorithm, based on the use of block FFT’s, for the computation of Ramaswami’s formula. The proposed algorithm is particularly effective when many block components of the probability invariant vector must be computed.

Key words. M/G/1 type matrices, Ramaswami’s formula, Toeplitz matrices, Fast Fourier Transform.

1 INTRODUCTION

In this paper, in order to agree with the customary notation of polynomial computation theory, we adopt transposed notations with respect to those standard in Markov chain theory.

A stochastic matrix of M/G/1 type [14] is an infinite row stochastic matrix P^{T} of the form

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

where B_{i}, A_{i}, i ≥ 0, are k × k matrices. Observe that, except for the first block column, P is a block Toeplitz matrix, i.e., a matrix with block entries that depend only on the difference of their indices. Suppose that P is irreducible and positive recurrent [7] and denote by π its invariant probability vector, that is, the solution to the equation

P π = π e^{T} π = 1,  (2)

where here and hereafter e denotes the vector, of suitable dimension, having all the components equal to 1. According to the structure of P, we partition the vector π into the vectors π_{0}, π_{1}, π_{2}, . . .. The block components π_{1}, π_{2}, π_{3}, … can be computed by directly solving the above system (for example by means of the methods proposed in [2], [17]), or by applying a numerically stable recursive formula, once the solution of a suitable nonlinear matrix equation has been calculated. More precisely, let G be the column stochastic matrix given by the minimal nonnegative solution of the matrix equation

G = ∑_{i=0}^{+∞} G^{i}A_{i},  (3)

and define the matrices

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

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

B̄_{0}π_{0} = π_{0}, k^{T} π_{0} = 1,

where

k^{T} = e^{T} + e^{T} (I − ∑_{i=1}^{+∞} Ā_{i})^{−1} ∑_{i=1}^{+∞} B̄_{i},

the vectors π_{i}, for i ≥ 1, are given by Ramaswami’s formula [15]

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

As there exist many efficient numerical methods to solve the matrix equation (3) (see [3], [8], [10], [13]) and since the matrix G is useful in calculating performance measures [14], Ramaswami’s formula is an effective tool for the computation of the probability invariant vector associated with M/G/1 type Markov chains.

Formula (5) was first proved by Ramaswami in [15], by using results from Markov renewal theory (Ramaswami’s proof may also be found in [14]). In [16], formula (5) was proved by using a different technique, based on the computation of the Schur complement of a leading principal submatrix of I − P. In both the approaches, the proof consists in building, by suitably truncating the matrix P of (1), the sequence {P^{(n)}}_{n} of column stochastic (n+1)×(n+1) block matrices, whose probability invariant vector is given, up to a multiplicative constant, by the corresponding components of the infinite vector π. Specifically, in [15] and [16], by using different tools, the authors prove the following result, which underlies Ramaswami’s formula.

Theorem 1.1 Let n ∈ N and consider the matrix

P^{(n)} = ( \begin{matrix}B_{0} & A_{0} \\ B_{1} & A_{1} & A_{0} & & ◯ \\ ⋮ & ⋮ & ... & ... \\ B_{n−1} & A_{n−1} & … & A_{1} & A_{0} \\ B̄_{n} & Ā_{n} & … & Ā_{2} & Ā_{1} \end{matrix} ) ,

where Ā_{i} and B̄_{i} are defined in (4). Then,

– the matrix P^{(n)} is an irreducible and positive recurrent column stochastic matrix;

– the probability invariant vector y^{(n)} associated with P^{(n)}, made up by the block components y_{0}^{(n)}, y_{1}^{(n)}, …, y_{n}^{(n)}, satisfies y_{i}^{(n)} = γ_{n}π_{i}, 0 ≤ i ≤ n, for a suitable constant γ_{n} > 0.

By using generating functions associated with triangular Toeplitz matrices, we give a novel interpretation of Ramaswami’s formula and a new fast algorithm, based on the use of FFT techniques, for its computation. More precisely, by means of the natural isomorphism between formal power series and triangular Toeplitz matrices, we find a block UL factorization (U is a block upper triangular matrix and L is a block lower triangular matrix) of the block Hessenberg matrix

H = ( \begin{matrix}I − A_{1} & −A_{0} & & ◯ \\ −A_{2} & I − A_{1} & −A_{0} \\ −A_{3} & −A_{2} & I − A_{1} & −A_{0} \\ ⋮ & ⋮ & ... & ... & ... \end{matrix} ) .  (6)

Then we derive the formula of Ramaswami by solving the system (I − P)π = 0 by means of that block UL factorization. The computation of the vector π is so reduced to solving a block lower triangular block Toeplitz infinite system. If that system is solved by means of forward substitution we obtain formula (5). By using this technique, the computational cost to calculate p block components of the probability invariant vector π grows in the variable p as 2pm_{0}k^{2} arithmetic operations, where m_{0} is an integer such that for i > m_{0} the matrices A_{i} are negligible. If the system is solved by FFT-based computations, we obtain a new algorithm, with a substantial reduction of the asymptotic computational cost. More precisely, the computational cost to calculate p block components of the probability invariant vector π grows in p at most as p(7k^{2} + 40 log m_{0}) arithmetic operations. We observe that the advantages are particularly strong in the case where the blocks A_{i} are negligible only for large values of i and in the case where many block components of the vector π need to be computed.

With regards to numerical stability, the customary Ramaswami formula is strongly numerically stable, since all the computations require only additions of nonnegative numbers, multiplications and divisions. The FFT-based version is weakly numerically stable [11], i.e., we are able to give upper bound only to the norm of the relative errors generated in the computations. However, in many numerical experiments performed, we detected no appreciable differences between the results computed with the two algorithms.

The paper is organized as follows. In Section 2 we derive the Ramaswami formula by means of the block UL factorization of the matrix H. In Section 3 we introduce the basic tools for FFT’s computations and we propose the new fast algorithm for the computation of the vector π. In Section 4 we present some numerical results.

2 RAMASWAMI’S FORMULA AND BLOCK UL FACTORIZATION

Let A be an infinite block lower triangular block Toeplitz matrix, defined by the k × k matrices α_{0}, α_{1}, α_{2}, …, i.e.

A = ( \begin{matrix}α_{0} \\ α_{1} & α_{0} & & ◯ \\ α_{2} & α_{1} & α_{0} \\ α_{3} & α_{2} & α_{1} & α_{0} \\ ⋮ & ... & ... & ... & ... \end{matrix} ) .

We associate the matrix A with the formal k × k matrix power series φ_{A}(z) = ∑_{i=0}^{+∞} α_{i}z^{i}, called the generating function for the matrix A. Generating functions allow us to represent row by column multiplication and inversion of block lower triangular block Toeplitz matrices. In fact, the following results can be easily shown

Proposition 2.1 Let A and B be two infinite block lower triangular block Toeplitz matrices, with blocks of the same dimension. Let φ_{A}(z) and φ_{B}(z), respectively, be their matrix-generating functions. Then the matrix C = AB is a block lower triangular block Toeplitz matrix whose matrix-generating function φ_{C}(z) satisfies φ_{C}(z) = φ_{A}(z)φ_{B}(z), the matrix product of φ_{A}(z) and φ_{B}(z).

Corollary 2.2 Let A be an infinite block lower triangular block Toeplitz matrix with matrix-generating function φ_{A}(z). If A is nonsingular, i.e., det α_{0} ≠ 0, then the inverse A^{−1} is a block lower triangular block Toeplitz matrix and φ_{A^{−1}}(z) satisfies φ_{A}(z)φ_{A^{−1}}(z) = I.

Consider the block Toeplitz block lower triangular matrix, defined by the second block column of the matrix I −P, where P is given by (1), and associate with it the matrix generating function φ(z) = zI − ∑_{i=0}^{+∞} A_{i}z^{i}. The following result holds:

Proposition 2.3 Let G be the minimal nonnegative solution of the matrix equation (3) and let Ā_{i}, for i ≥ 1, be defined in (4). Then

φ(z) = (zI − G)(I − ∑_{i=0}^{+∞} Ā_{i+1}z^{i}).  (7)

Proof. Immediate, by comparing the coefficients of z^{i} in both sides of (7), since A_{i} = Ā_{i} − GĀ_{i+1}, for i ≥ 1. □

For Proposition 2.1, the functional decomposition given in Proposition 2.3 can be rephrased in terms of block lower triangular block Toeplitz matrices factorization.

Corollary 2.4 The matrix H of (6) can be factorized as the product H = UL, where

U = ( \begin{matrix}I & −G & & & ◯ \\ & I & −G \\ & & I & −G \\ ◯ & & & ... & ... \end{matrix} ) L = ( \begin{matrix}I − Ā_{1} & & & ◯ \\ −Ā_{2} & I − Ā_{1} \\ −Ā_{3} & −Ā_{2} & I − Ā_{1} \\ ⋮ & ... & ... & ... \end{matrix} ) .

Proof. The proof readily follows by taking the matrix equation corresponding to decomposition (7) and by discarding the first block row. □

Ramaswami’s formula can now be derived by solving the system

UL ( π_{1} ; π_{2} ; π_{3} ; ⋮ ) = ( B_{1} ; B_{2} ; B_{3} ; ⋮ ) π_{0}.

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

L ( π_{1} ; π_{2} ; π_{3} ; ⋮ ) = ( B̄_{1} ; B̄_{2} ; B̄_{3} ; ⋮ ) π_{0}.  (8)

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

3 RAMASWAMI’S FORMULA COMPUTATION BY MEANS OF FFT-BASED TECHNIQUES

The block lower triangular block Toeplitz system (8) can be solved by using FFT-based techniques, thus obtaining a fast and weakly numerically stable method for the computation of the probability invariant vector of (2). Before explaining the method, we introduce some basic tools for Toeplitz matrices computations by means of FFT’s.

3.1 TOEPLITZ MATRICES COMPUTATION

Matrix polynomials computation and block Toeplitz matrices computation are closely related. Specifically,

• the product of an (m/2)×(m/2) block Toeplitz matrix and a block vector can be performed by means of a matrix polynomial multiplication modulo z^{m} − 1;

• the inversion of an m × m block lower triangular block Toeplitz matrix can be performed by inverting a matrix polynomial modulo z^{m}.

We will therefore present fast algorithms for matrix polynomial computations, by extending to the block case the FFT-based algorithms for scalar polynomials (see [4], [3]).

Throughout this section we will denote y = (y_{0}, …, y_{d−1}) = DFT(x) the Discrete Fourier Transform of the vector x = (x_{0}, …, x_{d−1}) defined by y_{i} = ∑_{j=0}^{d−1} ω^{ij}x_{j}, i = 0, …, d − 1, where ω = cos\frac{2π}{d} + i sin\frac{2π}{d}, and i is the imaginary unit such that i^{2} = −1. Analogously x = IDFT(y) denotes the Inverse Discrete Fourier Transform such that x_{i} = \frac{1}{d} ∑_{j=0}^{d−1} ω^{−ij}y_{j}, i = 0, …, d − 1. It is well known that if d = 2^{M}, where M is a positive integer, the computation of DFT and IDFT can be performed with O(d log d) real arithmetic operations (hereafter denoted by ops) by means of the base-2 FFT algorithms [4]. More precisely, in the case of real input the cost of DFT and IDFT is \frac{5}{2}d log d + O(1) ops and \frac{5}{2}d log d + d + O(1) ops, respectively, if we do not count the cost of computing the d-th roots of 1.

Similarly, we extend the above notation to block vectors. Given the block vector X = (X_{0}, …, X_{d−1}), where X_{s} = (x_{i,j}^{(s)}) are k×h matrices, we define Y = (Y_{0}, …, Y_{d−1}) = DFT(X) the block vector such that Y_{s} = (y_{i,j}^{(s)}) are k × h matrices and (y_{i,j}^{(0)}, …, y_{i,j}^{(d−1)}) = DFT(x_{i,j}^{(0)}, …, x_{i,j}^{(d−1)}), i = 1, …, k, j = 1, …, h . Analogously we set X = IDFT(Y) if (x_{i,j}^{(0)}, …, x_{i,j}^{(d−1)}) = IDFT(y_{i,j}^{(0)}, …, y_{i,j}^{(d−1)}). In this way the computation of DFT(X) and IDFT(Y) is reduced to the computation of kh DFT’s and IDFT’s, respectively, for the cost of O(khd log d) ops.

Consider the problem of the matrix polynomial multiplication modulo z^{m}−1. Let P_{i}, i = 0, …, m − 1, be k × l matrices, and let Q_{i}, i = 1, …, m − 1, be l × h matrices. Let P (z) = ∑_{i=0}^{m−1} P_{i}z^{i}, Q(z) = ∑_{i=0}^{m−1} Q_{i}z^{i}, R(z) = ∑_{i=0}^{m−1} R_{i}z^{i} be matrix polynomials such that R(z) = P (z)Q(z) mod (z^{m} − 1). Then the k × h matrix coefficients of R(z) can be computed, given the matrix coefficients of P (z) and Q(z), by means of the following algorithm.

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

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

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

Computation.

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

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

2. Compute the m matrix products W_{s} = U_{s}V_{s}, s = 0, …, m − 1, such

that W_{s} = P (ω^{s})Q(ω^{s}) = R(ω^{s}).

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

The cost of Algorithm 3.1 is O(klhm + (kl + lh + kh)m log m) ops. More precisely, for real input matrices the cost is less than 3(kl + lh + kh)m log m + m(3klh + 2kh). Indeed, due to the real input we have W_{s} = W̄_{m−s}, s = 1, …, m/2 − 1, where W̄_{m−s} is the complex conjugate of W_{m−s}, moreover W_{0}, W_{m/2} are real matrices. Thus the computation at stage 2 is reduced to m/2 − 1 complex matrix products and 2 real matrix products. Complex matrix products are computed by means of the algorithm which uses 3 multiplications and 5 additions [6].

Observe that Algorithm 3.1 can be also used to compute efficiently the product of a block Toeplitz matrix and a vector. In fact, if P (z) has degree m−1 and Q(z) has degree m/2−1, then by comparing the terms of degree m/2, …, m−1 in the equation R(z) = P (z)Q(z) mod (z^{m} − 1), we get the following equation involving a block Toeplitz matrix

( R_{m/2} ; ⋮ ; R_{m−1} ) = ( \begin{matrix}P_{m/2} & … & P_{1} \\ ⋮ & ... & ⋮ \\ P_{m−1} & … & P_{m/2} \end{matrix} ) ( Q_{0} ; ⋮ ; Q_{m/2−1} ) .

Thus, that product can be computed by means of Algorithm 3.1, i.e., by evaluating a product of two matrix polynomials. The polynomials P (z) and Q(z) have coefficients P_{0}, P_{1}, …, P_{m−1}, Q_{0}, Q_{1}, …, Q_{m−1}, where P_{0} = Q_{m/2} = … = Q_{m−1} = O.

Consider now the inversion of a polynomial modulo z^{m}. Suppose that P_{i} and Q_{i}, i = 0, …, m − 1, are square matrices of order k. Comparing the coefficients of z^{i} in both sides of the equation P (z)Q(z) = I mod z^{m} yields the equation involving a lower block triangular block Toeplitz matrix

( \begin{matrix}P_{0} & & & ◯ \\ P_{1} & P_{0} \\ ⋮ & ... & ... \\ P_{m−1} & … & P_{1} & P_{0} \end{matrix} ) ( Q_{0} ; Q_{1} ; ⋮ ; Q_{m−1} ) = ( I ; O ; ⋮ ; O ) .  (9)

Observe that, since the inverse of a lower block triangular block Toeplitz matrix is also lower block triangular block Toeplitz, the block vector (Q_{0}^{T} , Q_{1}^{T} , …, Q_{m−1}^{T})^{T} determines the inverse of the matrix in (9).

Denote by T_{m} the matrix in (9), assume for simplicity that m = 2^{M}, M positive integer, and partition T_{m} as

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

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

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

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

Algorithm 3.2 Solution of the congruence P (z)Q(z) = I mod z^{m}, or, equivalently, of the system (9).

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

Output. The k × k matrices Q_{0}, Q_{1}, …, Q_{m−1} satisfying (9), or equivalently, such that the polynomial Q(z) = ∑_{i=0}^{m−1} Q_{i}z^{i} solves the congruence

P (z)Q(z) = I mod z^{m}.

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 (10) with m = 2^{i+1} in the following way: by applying Algorithm 3.1 compute the products W = H_{2^{i}}U and V = T_{2^{i}}^{−1}W.

Observe that, at the i-th stage of Algorithm 3.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 k × k matrices. The overall cost is about 6mk^{3} + \frac{25}{2}k^{2}m log m ops.

3.2 COMPUTATION OF RAMASWAMI’S FORMULA

Observe that, since the matrices ∑_{i=0}^{+∞} A_{i} and ∑_{i=0}^{+∞} B_{i} are stochastic, the matrices A_{i} and B_{i} are negligible for large values of i. Therefore, the blocks A_{i} and B_{i} can be set equal to zero for indices i > m_{0}, where m_{0} is such that ∑_{i=0}^{m_{0}} A_{i} and ∑_{i=0}^{m_{0}} B_{i} are numerically stochastic (a nonnegative matrix A is said numerically stochastic if e^{T} (I − A) < μe^{T}, where μ is the machine precision). The matrices Ā_{i} and B̄_{i} are then also equal to zero for i > m_{0}, and (8) reduces to a block banded system with bandwidth m_{0}. Solving the banded system (8) by forward substitution leads to the Ramaswami algorithm:

Algorithm 3.3 CRF – Customary Ramaswami’s Formula

Input. Integers p and m_{0}, the matrix G solution of (3) and the matrices A_{i}, B_{i}, i = 0, …, m_{0}, such that ∑_{i=0}^{m_{0}} A_{i}, ∑_{i=0}^{m_{0}} B_{i} are numerically stochastic.

Output. Components π_{0}, π_{1}, …, π_{p} of the probability invariant vector π solving (2) and the value σ(p) = ∑_{i=0}^{p} e^{T} π_{i}.

Computation.

1. Computation of matrices Ā_{i}, i = 0, …, m_{0}, B̄_{i}, i = 1, …, m_{0}:

(a) Set Ā_{m_{0}} = A_{m_{0}}, B̄_{m_{0}} = B_{m_{0}}.

(b) For i = m_{0} − 1, …, 1 compute Ā_{i} = A_{i} + GĀ_{i−1} and B̄_{i} = B_{i} + GB̄_{i−1}.

(c) Compute B̄_{0} = B_{0} + GB̄_{1}.

2. Computation of π_{0}:

(a) Solve the system B̄_{0}π_{0} = π_{0}, k^{T} π_{0} = 1, where k^{T} = e^{T} +e^{T} (I − ∑_{i=1}^{m_{0}} Ā_{i})^{−1}(∑_{i=1}^{m_{0}} B̄_{i}).

(b) Set σ(p) = e^{T} π_{0}.

3. Computation of π_{1}, …, π_{p}:

(a) Compute π_{1} = (I − Ā_{1})^{−1}B̄_{1}π_{0};

(b) For i = 2, …, p, compute

\begin{matrix}π_{i} = (I − Ā_{1})^{−1} [B̄_{i}π_{0} + ∑_{j=j_{0}}^{i−1} Ā_{i+1−j}π_{j}] \\ σ(p) = σ(p) + e^{T} π_{i} \end{matrix}

where j_{0} = min(1, i + 1 − m_{0}) and B̄_{i} = O for i > m_{0}.

The computational cost of step 1 and 2 is about 4m_{0}k^{3} ops; the computational cost of step 3, if p > m_{0}, is about k^{3} + k^{2}m_{0}^{2} + 2k^{2}m_{0}(p − m_{0}) ops. Therefore, the asymptotic computational cost to calculate p components of π, where p > m_{0}, by means of Algorithm 3.3 is 2pk^{2}m_{0}. That computational cost can be substantially reduced by solving system (8) by the FFT-based algorithms for Toeplitz computation described in Section 3.1. Suppose, without loss of generality, that m_{0} = 2^{M} (if m_{0} is not a power of 2, we truncate the sequences A_{i}, B_{i} at the least power of 2 greater than m_{0}). Let us partition the matrix L of (8) into blocks of dimension km_{0}, where each block is itself partitioned into m_{0} sub-blocks of dimension k. So, the matrix L can be rewritten as

L = ( \begin{matrix}I − R & & & ◯ \\ −S & I − R \\ & −S & I − R \\ ◯ & & ... & ... \end{matrix} ) ,  (11)

where

R = ( \begin{matrix}Ā_{1} & & & & ◯ \\ Ā_{2} & Ā_{1} \\ Ā_{3} & Ā_{2} & Ā_{1} \\ ⋮ & ... & ... & ... \\ Ā_{m_{0}} & Ā_{m_{0}−1} & … & Ā_{2} & Ā_{1} \end{matrix} ) , S = ( \begin{matrix}O & Ā_{m_{0}} & … & Ā_{3} & Ā_{2} \\ & O & Ā_{m_{0}} & … & Ā_{3} \\ & & ... & ... & ⋮ \\ & & & O & Ā_{m_{0}} \\ ◯ & & & & O \end{matrix} )

are block triangular block Toeplitz matrices. Let us partition the vector π according to the structure (11) of L, i.e., π^{T} = (π̃_{0}^{T} , π̃_{1}^{T} , π̃_{2}^{T} , …)^{T}, where π̃_{0} = π_{0}, π̃_{i}^{T} = (π_{m_{0}(i−1)+1}^{T}, π_{m_{0}(i−1)+2}^{T}, …, π_{m_{0}(i−1)+m_{0}}^{T})^{T}. With these notations the system (8) can be solved by means of the equations

\begin{matrix}π̃_{1} = (I − R)^{−1}b \\ π̃_{i} = (I − R)^{−1}Sπ̃_{i−1}, i ≥ 2, \end{matrix}  (12)

where

b = ( B̄_{1} ; B̄_{2} ; ⋮ ; B̄_{m_{0}} ) π_{0}.  (13)

The matrices involved in (12) are block triangular block Toeplitz matrices. The computations can be performed by the FFT-based methods of Section 3.1, leading to the following fast algorithm for π:

Algorithm 3.4 FRF – Fast Ramaswami’s Formula

Input. Integers n and m_{0}, the matrix G solution of (3) and the matrices A_{i}, B_{i}, i = 0, …, m_{0}, such that ∑_{i=0}^{m_{0}} A_{i}, ∑_{i=0}^{m_{0}} B_{i} are numerically stochastic.

Output. Components π_{0}, π_{1}, …, π_{p} of the probability invariant vector π, solving (2), where p = nm_{0}, and the value σ(p) = ∑_{i=0}^{p} e^{T} π_{i}.

Computation.

1. Computation of matrices Ā_{i}, i = 0, …, m_{0}, B̄_{i}, i = 1, …, m_{0}:

(a) Set Ā_{m_{0}} = A_{m_{0}}, B̄_{m_{0}} = B_{m_{0}}.

(b) For i = m_{0} − 1, …, 1 compute Ā_{i} = A_{i} + GĀ_{i−1}, B̄_{i} = B_{i} + GB̄_{i−1}.

(c) Compute B̄_{0} = B_{0} + GB̄_{1}.

2. Computation of π_{0}:

(a) Solve the system B̄_{0}π_{0} = π_{0}, k^{T} π_{0} = 1, where k^{T} = e^{T} +e^{T} (I − ∑_{i=1}^{m_{0}} Ā_{i})^{−1}(∑_{i=1}^{m_{0}} B̄_{i}).

(b) Set σ(p) = e^{T} π_{0}.

3. Computation of π_{1}, …, π_{p}:

(a) Compute (I − R)^{−1} by Algorithm 3.2.

(b) Compute, by Algorithm 3.1, π̃_{1} = (I − R)^{−1}b, where b is given in (13), and compute σ(p) = σ(p) + e^{T} π̃_{1}.

(c) For i = 2, …, n, compute, by Algorithm 3.1,

s = Sπ̃_{i−1}, π̃_{i} = (I − R)^{−1}s,

and σ(p) = σ(p) + e^{T} π̃_{i}.

Steps 1 and 2 in Algorithm 3.4 are the same as Steps 1 and 2 in Algorithm 3.3, and have the computational cost of 4m_{0}k^{3} ops. At Step 3, the computational cost of calculating stage (a) and (b) is about 6m_{0}k^{3} +35/2k^{2}m_{0} log m_{0}; at stage (c), once the DFT of the matrix S is computed (with less than 10k^{2}m_{0} log m_{0} ops.), the computational cost of the generic block π̃_{i} is at most 7k^{2}m_{0} + 20km_{0} log m_{0} + 20km_{0} ≤ 7k^{2}m_{0} + 40km_{0} log m_{0} ops. Whence the overall cost to compute p = nm_{0} components is at most 10m_{0}k^{3} + 39m_{0} log m_{0} + n(7k^{2}m_{0} + 40km_{0} log m_{0}), i.e., the computational cost grows linearly in p at most as p(7k^{2} + 40k log m_{0}). Therefore Algorithm 3.4 is asymptotically (as p tends to infinity) 2k^{2}m_{0}/(7k^{2} + 40 log m_{0}) times faster than Algorithm 3.3. Moreover Algorithm 3.4 is particularly faster in cases where m_{0} is large, that is, in cases where the blocks A_{i} and B_{i} are negligible only for large values of i.

Concerning numerical stability, all computations in Algorithm 3.3 consist in performing additions of nonnegative numbers, multiplications and divisions, so that Algorithm 3.3 is strongly numerically stable, i.e., the relative errors generated by the floating point arithmetic can be bounded entry-wise. Algorithm 3.4, due to the use of FFT’s, is weakly numerically stable, i.e., it is only possible to give very good relative error bounds in terms of the norm (see [11], [4]). However, in spite of the weak stability of FFT, in numerous experiments we have performed, we detected no difference between the results of the algorithms.

4 NUMERICAL RESULTS AND COMPARISONS

We implemented Algorithms 3.3 and 3.4 in Fortran 77; the programs were run on an Alpha Workstation using standard double precision IEEE arithmetic. We tested the algorithms on a BMAP/G/1 queue (Example 1) and on a problem arising from modeling a Metaring MAC Protocol [1] (Example 2).

Example 1 We consider an 8-state BMAP/G/1 queue with constant service times. The blocks A_{i} and B_{i} have dimension 8 and have the structure described in [12], that is,

A_{i}^{T} = ∑_{j=0}^{+∞} e^{−θ} \frac{θ^{j}}{j!} K_{i}^{(j)}, B_{i}^{T} = −D_{0}^{−1} ∑_{j=0}^{i} D_{j+1}A_{i−j}^{T},

where θ = max_{h}{−(D_{0})_{hh}} and the {K_{i}^{(j)}} are recursively defined by K_{0}^{(0)} = I, K_{i}^{(0)} = O, i ≥ 1, and

\begin{matrix}K_{0}^{(j+1)} = K_{0}^{(j)}(I + θ^{−1}D_{0}), \\ K_{i}^{(j+1)} = θ^{−1} ∑_{h=0}^{i−1} K_{h}^{(j)}D_{i−h} + K_{i}^{(j)}(I + θ^{−1}D_{0}), i ≥ 1. \end{matrix}

In our particular example,

(D_{i})_{jh} = e^{−λ_{jh}} \frac{λ_{jh}^{i}}{i!} , i ≥ 1,

(D_{0})_{jh} = e^{−λ_{jh}}, j ≠ h,

(D_{0})_{hh} = − ∑_{j≠h}(D_{0})_{jh} − ∑_{i=1}^{+∞} ∑_{j=1}^{8}(D_{i})_{jh},

for suitable nonnegative parameters λ_{jh}. The integer m_{0} for which the matrices ∑_{i=0}^{m_{0}} A_{i}, ∑_{i=0}^{m_{0}} B_{i} are numerically stochastic is 64. Table 1 reports, for different numbers p of computed block components, the corresponding value σ(p) = ∑_{i=0}^{p} e^{T} π_{i}, the CPU times (in seconds) needed by Algorithm 3.3 and Algorithm 3.4, respectively, and their ratios. Figure 1 shows the CPU time

\begin{matrix}p & σ(p) & CRF (sec.) & F RF (sec.) & CRF/F RF \\ 2^{6} & 0.077379 & 0.01 & 0.07 & 0.1 \\ 2^{7} & 0.148238 & 0.03 & 0.08 & 0.4 \\ 2^{8} & 0.274049 & 0.05 & 0.09 & 0.6 \\ 2^{9} & 0.472667 & 0.11 & 0.12 & 0.9 \\ 2^{10} & 0.721746 & 0.22 & 0.17 & 1.3 \\ 2^{11} & 0.922526 & 0.41 & 0.27 & 1.5 \\ 2^{12} & 0.993994 & 0.84 & 0.59 & 1.4 \\ 2^{13} & 0.999963 & 1.77 & 1.03 & 1.7 \\ 2^{14} & 0.999999 & 3.73 & 1.74 & 2.1 \end{matrix}

Table 1: BMAP/G/1 Queue

needed by Algorithm 3.4 and Algorithm 3.3 to compute the same number of block components of the vector π.

Example 2 We consider a problem arising from modeling a metropolitan queueing problem (Metaring MAC Protocol [1]). The dimension k of the blocks A_{i}

Figure 1: BMAP/G/1 Queue

and B_{i} is 80, and the integer m_{0} such that the matrices ∑_{i=0}^{m_{0}} A_{i}, ∑_{i=0}^{m_{0}} B_{i} are numerically stochastic is 256. For a description of the model and the structure of the blocks we refer to [1]. This particular example is characterized by large block bandwidth of the matrix L in system (8), by a large dimension of the blocks A_{i} and B_{i}, and by a long queue, i.e., the blocks π_{i} tend to zero slowly. Table 2 reports, for different numbers p of block components, the corresponding value σ(p) = ∑_{i=0}^{p} e^{T} π_{i}, the CPU times (in seconds) needed by Algorithm 3.3 and Algorithm 3.4, respectively, and their ratios. Figure 2 shows the CPU time needed by Algorithm 3.4 and Algorithm 3.3 to compute the same number of block components of the vector π.

Observe that, in both the examples, the customary Ramaswami formula is faster than Fast Ramaswami formula in the case where only few block components of the vector π need to be computed. Fast Ramaswami formula really shows its effectiveness when many block components must be computed. When the block dimension k is large and the number m_{0} of nonnegligible blocks A_{i} and B_{i} is large, the advantages of the FRF method are clear.

References

[1] Anastasi G., Lenzini L., Meini B., “Performance evaluation of a worst

\begin{matrix}p & σ(p) & CRF (sec.) & F RF (sec.) & CRF/F RF \\ 2^{9} & 0.285677 & 181 & 374 & 0.5 \\ 2^{10} & 0.457895 & 385 & 399 & 1.0 \\ 2^{11} & 0.687780 & 797 & 424 & 1.9 \\ 2^{12} & 0.896434 & 1681 & 454 & 3.7 \\ 2^{13} & 0.988604 & 3236 & 487 & 6.6 \\ 2^{14} & 0.999862 & 6632 & 652 & 10.2 \\ 2^{15} & 0.999999 & 13322 & 829 & 16.1 \end{matrix}

Table 2: Metaring MAC Protocol

case model of the Metaring MAC Protocol with global fairness”, Performance Evaluation, forthcoming, 1997.

[2] Bini D., Meini B. “On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems”, Computations with Markov Chains, W. J. Stewart Ed., pp. 21–38, Kluwer Academic Publishers, 1996.

[3] Bini D., Meini B., “On the solution of a nonlinear matrix equation arising in queueing problems”, SIAM J. Matrix Anal. Appl., 17 (1996), to appear.

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

[5] Bini D., Pan V., “Polynomial division and its computational complexity”, J. Complexity, 6 (1986), pp. 179–203.

[6] Borodin A., Munro I., The Computational Complexity of Algebraic and Numeric Problems, American Elsevier, New York, 1975.

[7] Çinlar E., Introduction to Stochastic Processes, Prentice-Hall, Englewood Cliffs, N.J., 1975.

[8] Latouche G., “Algorithms for evaluating the matrix G in Markov chains of PH/G/1 type”, Bellcore Technical Report, 1992.

[9] Latouche G., Ramaswami V., “A logarithmic reduction algorithm for quasi-birth-death processes”, J. Appl. Prob., 30 (1993), pp. 650–674.

[10] Latouche G., Stewart G. W., “Numerical methods for M/G/1 type queues”, Computations with Markov Chains, W. J. Stewart Ed., pp. 571– 581, Kluwer Academic Publishers, 1996.

[11] Linzer E., “On the stability of transform-based circular deconvolution”, SIAM J. Numer. Anal., 29 (1992), pp. 1482–1492.

Figure 2: Metaring MAC Protocol

[12] Lucantoni D. M., “New results on the single server queue with a batch Markovian arrival process”, Commun. Statist. Stochastic Models, 7 (1991), pp. 1–46.

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

[14] Neuts M. F., Structured Stochastic Matrices of M/G/1 Type and their Applications, Marcel Dekker Inc., New York, 1989.

[15] Ramaswami V., “A stable recursion for the steady state vector in Markov chains of M/G/1 type”, Commun. Statist. Stochastic Models, 4 (1988), pp. 183-188.

[16] Schellhaas H., “On Ramaswami’s algorithm for the computation of the steady state vector in Markov chains of M/G/1 type”, Commun. Statist. Stochastic Models, 6 (1990), pp. 541–550.

[17] Stewart G. W., “On the solution of block Hessenberg systems”, Numer. Lin. Algebra with Appl., 2 (1995), pp. 287–296.