On functional iteration methods for solving nonlinear matrix equations arising in queueing

problems

Paola Favati^{∗}  Beatrice Meini^{†}

Abstract

The problem of the computation of the minimal nonnegative solution G of the nonlinear matrix equation X = ∑_{i=0}^{+∞} X^{i}A_{i} is considered. This problem arises in the numerical solution of M/G/1 type Markov chains, where A_{i}, i ≥ 0, are nonnegative k × k matrices such that ∑_{i=0}^{+∞} A_{i} is column stochastic. We analyze classical functional iteration methods, by estimating the rate of convergence, in relation with spectral properties of the starting approximation matrix X_{0}. Based on these new convergence results, we propose an effective method to choose a matrix X_{0}, which drastically reduces the number of iterations; the additional cost needed to compute X_{0} is much less than the overall savings achieved by reducing the number of iterations.

1 Introduction

Many queueing problems can be modeled by Markov chains of M/G/1 type, that is, Markov chains whose probability transition matrix has the structure

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

where B_{i+1}, A_{i}, i ≥ 0, are k × k nonnegative matrices such that ∑_{i=1}^{+∞} B_{i} and ∑_{i=0}^{+∞} A_{i} are column stochastic. One of the major problems in Markov chains is the computation of the probability invariant vector π, i.e., the solution of the infinite system π = P^{t}π, ∥π∥_{1} = 1. For the matrices of structure (1), the computation of π can be reduced (compare [16]) to the computation of the minimal nonnegative solution G of the nonlinear matrix equation

∗Istituto di Matematica Computazionale del C.N.R., via S. Maria 46, 56127 Pisa, Italy.

E-mail: favati@imc.pi.cnr.it

†Dipartimento di Matematica, via Buonarroti 2, 56127 Pisa, Italy. E-mail:

meini@dm.unipi.it

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

where X is a k × k matrix. Indeed, once G is known, the vector π can be recovered by means of the Ramaswami formula [17]. Moreover, the knowledge of the matrix G is very important, due to its role in calculating performance measures [16].

Recently, the analysis of numerical methods for computing the solution of the matrix equation (2) has been largely developed. New quadratically convergent methods have been proposed in [12], [13], [11], [3], [4], [5]; classical linearly convergent methods based on functional iterations have been improved and studied in details in [8], [14]. The latter methods, unlike most of the quadratically convergent methods, are very easy to implement. They are based on the recursion

X_{n+1} = F (X_{n}), n ≥ 0,  (3)

where

F (X) = K(X)(I − H(X))^{−1},  (4)

and H(X) = ∑_{i=0}^{+∞} X^{i}H_{i}, K(X) = ∑_{i=0}^{+∞} X^{i}K_{i} are such that H_{i} ≥ O, K_{i} ≥ O, i ≥ 0, and ∑_{i=0}^{+∞} z^{i}A_{i} = z ∑_{i=0}^{+∞} z^{i}H_{i} + ∑_{i=0}^{+∞} z^{i}K_{i} [8], [14]. It can be proved that, by suitable choosing the starting approximation X_{0}, the sequence (3) linearly converges to the solution G of (2); moreover, the rate of convergence is strongly related to the choice of X_{0}. More precisely, we associate with any function F (·) of (4) a suitable nonnegative k^{2} × k^{2} matrix R̂_{F}. In the case where X_{0} = 0 in [14] it is shown that the asymptotic rate of convergence of {X_{n}} coincides with the spectral radius of R̂_{F}. Moreover, in the case where X_{0} is any column stochastic matrix, it is shown that at each step n, the vector ê_{n}, obtained by the column-wise arrangement of the entries of the matrix G − X_{n}, is orthogonal to the nonnegative left eigenvector v^{T} associated with the spectral radius of the matrix R̂_{F}. From this property, it follows that the rate of convergence of X_{n} is less than or equal to the spectral radius of the matrix M_{1} obtained by projecting R̂_{F} onto a suitable subspace V_{1}.

From the above results it follows that, if the matrix R̂_{F} is irreducible, then the mean asymptotic rate of convergence of the sequence X_{n} is:

– equal to the spectral radius of R̂_{F}, if X_{0} = O;

– less than or equal to the second largest modulus eigenvalue of R̂_{F}, if X_{0} is a column stochastic matrix.

The better rate of convergence, when X_{0} is a stochastic matrix, is mainly due to the fact that both G and X_{0} have the eigenvalue 1 and associated left eigenvector e^{t} = (1, …, 1).

In this paper, in order to further improve the rate of convergence, we consider starting approximations X_{0} that share with G more eigenvalues, and the corresponding left eigenvectors, in addition to the eigenvalue 1 and the eigenvector e^{t}.

More precisely, let λ_{1} = 1, λ_{2}, …, λ_{h}, h ≤ k, be eigenvalues of the matrix G such that the corresponding left eigenvectors w_{1}^{t} = e^{t}, w_{2}^{t}, …, w_{h}^{t} are linearly independent. Suppose that the matrix X_{0} satisfies

w_{i}^{t}X_{0} = λ_{i}w_{i}^{t}, i = 1, …, h,  (5)

that is, it shares with G the eigenvalues λ_{1}, λ_{2}, …, λ_{h} and the corresponding left eigenvectors. Then we prove that the rate of convergence of the sequence X_{n} is bounded from above by ρ(M_{h}), where ρ(A) denotes the spectral radius of the matrix A and M_{h} is the matrix obtained by projecting R̂_{F} onto a suitable subspace V_{h}. Moreover, we show that ρ(M_{h}) ≤ ρ(M_{1}), for any h ∈ {1, …, k}. In this way the rate of convergence can be substantially reduced.

Once the eigenvalues λ_{1}, λ_{2}, …, λ_{h}, and the corresponding eigenvectors w_{1}^{t}, …, w_{h}^{t} are known, a matrix X_{0} satisfying (5) can simply be obtained by means of the equation X_{0} = V Diag(λ_{1}, …, λ_{h})W^{t}, where W is the k×h matrix whose columns are the vectors w_{1}, …, w_{h}, and V is the generalized inverse of W^{t}. When h = k, i.e., all the eigenvalues, and the corresponding left eigenvectors, are known, the matrix V coincides with the inverse of W^{t}. In this case some authors [9], [15] have proposed to approximate G by simply computing W^{−t}Diag(λ_{1}, …, λ_{k})W^{t}. In practice, due to the roundoff errors generated in the computation, the matrix G̃ obtained in this way can be a poor approximation of G; we propose to improve this approximation by choosing X_{0} = G̃ and by applying few steps of the functional iteration method (3).

We have performed several numerical experiments. We have first computed eigenvalues and corresponding left eigenvectors, for different values of h. Then a matrix X_{0} satisfying (5) has been computed and functional iteration has been applied. From the many results, we have observed a drastic reduction of the number of iterations needed to approximate G when h > 1, with respect to the case where h = 1, or X_{0} is a generic column stochastic matrix. On the other hand, the time needed to approximate the eigenvalues, and the corresponding left eigenvectors, when they are not known, is negligible with respect to the time needed to perform functional iterations.

The paper is organized as follows: in Section 2 we prove the main convergence results, by estimating the rate of convergence, in relation with the starting approximation X_{0}; in Section 3 we propose a method for computing a matrix X_{0} satisfying (5); in Section 4 we present numerical results and comparisons.

2 Convergence rate and starting approximation

In this section we extend the convergence results proved in [14] for functional iteration methods based on the recursion (3), for solving the nonlinear matrix equation (2); these new convergence results suggest how to choose the starting approximation matrix X_{0} in order to drastically reduce the number of iterations needed to approximate the solution G.

Let us denote by ê_{n} the vector obtained by column-wise arranging the entries of the matrix G − X_{n}. It can be shown [14] that the error ê_{n+1} is related to the error ê_{n} by means of the recursion

ê_{n+1} = R̂_{n}ê_{n}  (6)

where

R̂_{n} = ∑_{j=0}^{+∞} Y_{j,n} ⊗ G^{j}  (7)

and

\begin{matrix}Y_{0,n} = ∑_{i=1}^{+∞} (X_{n}^{i−1}K_{i}(I − H(X_{n}))^{−1})^{t} \\ Y_{j,n} = ∑_{i=j}^{+∞} (X_{n}^{i−j}A_{i+1}(I − H(X_{n}))^{−1})^{t} , j ≥ 1. \end{matrix}

Denote by

R̂_{F} = ∑_{j=0}^{+∞} Y_{j} ⊗ G^{j}

where

\begin{matrix}Y_{0} = ∑_{i=1}^{+∞} (G^{i−1}K_{i}(I − H(G))^{−1})^{t} \\ Y_{j} = ∑_{i=j}^{+∞} (G^{i−j}A_{i+1}(I − H(G))^{−1})^{t} , j ≥ 1. \end{matrix}

Observe that, if lim_{n} X_{n} = G, then R̂_{F} = lim_{n} R̂_{n}. In the following we denote by e_{i}, i = 1, …, k, the k-dimensional vector having the i-th entry equal to 1, and the remaining ones equal to 0.

Relation (6) allows us to prove the following convergence result [14]:

Theorem 1 Let

r = lim_{n} ^{n} \sqrt{∥ê_{n}∥}

be the mean asymptotic rate of convergence of the sequence X_{n+1} = F (X_{n}), where F (·) is defined by (4) and ∥ · ∥ is any vector norm. The following results hold:

– if X_{0} = 0 then r = ρ(R̂_{F});

– if X_{0} is a stochastic matrix then r ≤ ρ(T_{1}^{t}R̂_{F} T_{1}), where T_{1} is the k^{2} × (k^{2} − k) matrix whose columns complete the set of vectors

\frac{1}{\sqrt{k}} (e_{i} ⊗ e), i = 1, …, k,  (8)

to an orthonormal basis of R^{k^{2}}.

The reduction of the rate of convergence, in the case where X_{0} is a column stochastic matrix, is mainly due to the fact that X_{0} has the eigenvalue 1 and associated left eigenvector e^{t}, as well as the matrix G. In fact, from this property it follows that the vector ê_{n} is orthogonal to the nonnegative left eigenvector associated with the spectral radius of the matrix R̂_{F}.

This fact suggests us to choose a starting approximation X_{0} that shares with G some other eigenvalues and the corresponding left eigenvectors, in addition to the eigenvalue equal to 1 and the eigenvector e^{t}.

Denote by Λ = {λ_{1}, …, λ_{k}} the set of the eigenvalues of the matrix G, and by w_{1}^{t}, …, w_{k}^{t} the corresponding left eigenvectors, where λ_{1} = 1 and w_{1}^{t} = e^{t}. Let Λ_{h} be a subset of h eigenvalues of Λ, 1 ≤ h ≤ k, such that:

1. λ_{1} ∈ Λ_{h};

2. if λ ∈ Λ_{h}, then λ̄ ∈ Λ_{h}, where λ̄ denotes the complex conjugate of λ;

3. the left eigenvectors corresponding to the eigenvalues in Λ_{h} are linearly independent.

Throughout, we suppose, without loss of generality, that Λ_{h} = {λ_{1}, …, λ_{h}}. Denote by U_{h} the subspace generated by the k^{2}–dimensional vectors e_{i} ⊗w_{j},

i = 1, …, k, j = 1, …, h, and let V_{h} = U_{h}^{⊥}. Let S_{h} = I ⊗ E_{h} a k^{2} × hk matrix whose columns are an orthonormal basis of U_{h}, and let T_{h} = I ⊗ F_{h} a k^{2} × (k^{2} − hk) matrix whose columns are an orthonormal basis of V_{h}. Then the following result holds:

Theorem 2 Let X_{n}, n ≥ 0, be a convergent sequence generated by X_{n+1} = F (X_{n}), where F (·) is defined by (4). Denote by

r = lim_{n} ^{n} \sqrt{∥ê_{n}∥}

the mean asymptotic rate of convergence, where ∥ · ∥ is any vector norm. Let 1 ≤ h ≤ k and Λ_{h} ⊂ Λ satisfying properties 1, 2, 3. Suppose that X_{0} verifies

w^{t}X_{0} = λw^{t},  (9)

for any λ ∈ Λ_{h}, where w^{t} is the corresponding left eigenvector of G. Then it holds r ≤ ρ(T_{h}^{t}R̂_{F} T_{h}).

Proof. First observe that ê_{n} ∈ V_{h}, for n ≥ 0. In fact, for n = 0, this property follows from the relations w^{t}(G − X_{0}) = λw^{t} − λw^{t} = 0, for any λ ∈ Λ_{h}; for n ≥ 1, the property follows from (6) and from the fact that, if u ∈ U_{h}, then u^{t}R̂_{n} ∈ U_{h}, for any n ≥ 0.

Let P = (S_{h}|T_{h}). From (6) it follows that

P^{t}ê_{n+1} = ( \begin{matrix}0 \\ f_{n+1} \end{matrix} ) = (P^{t}R̂_{n}P)(P^{t}ê_{n}) = (P^{t}R̂_{n}P) ( \begin{matrix}0 \\ f_{n} \end{matrix} ) ,  (10)

where f_{n} = T_{h}^{t}e_{n}. Since S_{h}^{t}R̂_{n}T_{h} = O, it follows that

P^{t}R̂_{n}P = ( \begin{matrix}S_{h}^{t} \\ − \\ T_{h}^{t} \end{matrix} ) R̂_{n} (S_{h}|T_{h}) = ( \begin{matrix}V_{n} & O \\ T_{n} & Z_{n} \end{matrix} ) ,  (11)

where Z_{n} = T_{h}^{t}R̂_{n}T_{h}. From (10) and (11), the recurrence (6) can be rewritten in the form

f_{n+1} = Z_{n}f_{n}  (12)

and

r = lim_{n} ^{n} \sqrt{∥ê_{n}∥} = lim_{n} ^{n} \sqrt{∥f_{n}∥}.

By following the same lines of the proof of Theorem 15 in [14], we easily arrive at r ≤ ρ(T_{h}^{t}R̂_{F} T_{h}). □

From the above theorem it follows that, by starting with a matrix X_{0} sharing with G more than one eigenvalue, the rate of convergence can be reduced. In fact, the following result holds:

Theorem 3 Under the hypotheses of Theorem 2, the eigenvalues of the matrix T_{h}^{t}R̂_{F} T_{h} are given by the set

⋃_{i=h+1}^{k}{μ : μ is eigenvalue of Q_{i}}

where Q_{i} = ∑_{j=0}^{∞} λ_{i}^{j}Y_{j}, i = 1, …, k.

Proof. Observe that the matrix Z_{F} = T_{h}^{t}R̂_{F} T_{h} can be written as

Z_{F} = (I ⊗ F_{h}^{t})R̂_{F} (I ⊗ F_{h}) = ∑_{j=0}^{∞} Y_{j} ⊗ F_{h}^{t}G^{j}F_{h}.

Whence, the matrix Z_{F} is similar to the matrix ∑_{j=0}^{∞} F_{h}^{t}G^{j}F_{h} ⊗ Y_{j}, which is similar to ∑_{j=0}^{∞} U^{j} ⊗ Y_{j}, where U is the Schur canonical form of F_{h}^{t}GF_{h}. Since the eigenvalues of F_{h}^{t}GF_{h} are λ_{h+1}, …, λ_{k}, U is an upper triangular matrix whose diagonal entries are λ_{h+1}, …, λ_{k}. In this way we find that Z_{F} is similar to a block upper triangular matrix whose diagonal blocks are the matrices Q_{i}, i = h + 1, …, k. □

From the above theorem it follows that:

Corollary 4 Under the hypotheses of Theorem 2 it holds that

ρ(T_{h}^{t}R̂_{F} T_{h}) ≤ ρ(T_{1}^{t}R̂_{F} T_{1})

for any 1 ≤ h ≤ k.

In the above theorems we have presented conditions on the starting matrix X_{0} to increase the rate of convergence of the sequence X_{n+1} = F (X_{n}), under the hypothesis that X_{0} is such that the sequence X_{n} is convergent. In the case where the matrix X_{0} is stochastic, the sequence is obviously convergent. In the more general case where X_{0} is not necessarily nonnegative, the local convergence is guaranteed by the classical results on functional iteration methods [18], since ρ(R̂_{F}) < 1.

3 Choosing the starting approximation

Based on the results of Section 2, in order to find a good starting matrix X_{0}, we need first to know eigenvalues and left eigenvectors of the matrix G. The eigenvalues λ_{1}, …, λ_{k} of G are given by the zeros of the function

φ(z) = det(zI − ∑_{j=0}^{+∞} A_{j}z^{j})  (13)

lying in the closed unit circle of the complex plane (compare [9]), and the left eigenvectors are the solutions of the linear systems

w_{i}^{t}(λ_{i}I − ∑_{j=0}^{+∞} A_{j}λ_{i}^{j}) = 0, i = 1, …, k.  (14)

The zeros of the function φ(z) can be computed by means of any algorithm for approximating the zeros of a polynomial, for example by means of Aberth’s method [7], [1]; the left eigenvectors can be computed by solving the above homogeneous linear systems.

Concerning the choice of the number h of shared eigenvalues/eigenvectors, we are not able to establish in advance the optimal number, in order to obtain the best rate of convergence. In fact, the optimal number is strongly related to the spectral properties of G, and thus to the specific problems. For this reason, we suggest to compute a matrix X_{0}, that shares with G, in addition to 1, only the eigenvalues that can be computed in a few numbers of iterations. In this way the additional time needed to approximate the shared eigenvalues is kept low; moreover, a bad approximation of an eigenvalue, and thus of the corresponding eigenvector, may not to lead to an improvement of the rate of convergence.

Suppose we know a subset Λ_{h} = {λ_{1}, …, λ_{h}} of Λ satisfying properties 1, 2, 3 of Section 2; denote by w_{1}^{t}, …, w_{h}^{t} the corresponding left eigenvectors. Let W_{h} be the k × h matrix whose columns are the vectors w_{1}, …, w_{h} and let V_{h} be a k × h matrix such that

V_{h}^{t}W_{h} = I_{h},  (15)

where I_{h} is the h dimensional identity matrix. Then, the matrix

X_{0} = V_{h}Diag(λ_{1}, …, λ_{h})W_{h}^{t}  (16)

verifies conditions (9). A simple way to obtain a matrix V_{h} satisfying (15) is to compute the generalized inverse of W_{h}^{t}, say by means of QR or SVD method [10], or by means of the relation V_{h} = (W_{h}W_{h}^{t})^{−1}W_{h}, that holds since W_{h} is a full rank matrix. Since λ ∈ Λ_{h} implies λ̄ ∈ Λ_{h}, the matrix X_{0} obtained in this way is real. If, in addition, V_{h} is such that

V_{h}Diag(λ_{1}, …, λ_{h})W_{h}^{t} ≥ 0,  (17)

then X_{0} is stochastic, since X_{0} ≥ 0 and w_{1}^{t}X_{0} = w_{1}^{t}, where w_{1}^{t} = e^{t}.

A matrix V_{h} satisfying both conditions (15) and (17) can be viewed as a solution of a linear programming problem. However the computation of a solution, if it exists, can require a long time. On the other hand, in many numerical experiments, even if the starting matrix X_{0}, obtained by means of the generalized inverse, has negative entries after few steps the matrix X_{n} is nonnegative. Whence, the use of generalized inverse, leads to an advantageous choice of X_{0}.

Suppose we know all the eigenvalues λ_{1}, …, λ_{k} of G and the corresponding left eigenvectors w_{1}^{t}, …, w_{k}^{t}. If the above eigenvectors are linearly independent, then it is well known [9], [15] that an approximation of the matrix G can be obtained by means of the relation

G = W_{k}^{−t}Diag(λ_{1}, …, λ_{k})W_{k}^{t}.  (18)

However, due to errors generated in the computation, the matrix obtained by calculating W_{k}^{−t}Diag(λ_{1}, …, λ_{k})W_{k}^{t} can be a not sufficiently accurate approximation of G. In this case, in order to improve the approximation of G, we apply the functional iteration method, starting with X_{0} = W_{k}^{−t}Diag(λ_{1}, …, λ_{k})W_{k}^{t}. In most applications, the matrix X_{0} is nonnegative and generates a sequence X_{n} which rapidly converges to G.

4 Numerical results

We have performed numerical experiments in order to test the behavior of the rate of convergence with different starting approximations X_{0}. Since the matrix X_{0} of (16) is not necessarily nonnegative, the matrices X_{n} could have negative entries and the matrices I − H(X_{n}) of (4) could be singular, for small values of n. Thus, in order to avoid possibly singular matrix inversions, we use the standard functional iteration method X_{n+1} = ∑_{i=0}^{+∞} X_{n}^{i}A_{i}, n ≥ 0.

We have chosen matrices X_{0} sharing with G one, two, or more eigenvalues and corresponding left eigenvectors. We have stopped iterations when the residual error ||X_{n} − ∑_{i=0}^{+∞} X_{n}^{i}A_{i}||_{1} was less than ϵ = 10^{−11}.

We have implemented our method in Fortran 77; the program has been run on an Ultra Sparc Workstation.

We have considered two problems; Problem 1 arises from the mathematical modeling of a metropolitan queueing network (Metaring MAC Protocol) [2],

h Iterations Time (s.)
1  2763  197.3
2  34  4.9
4  26  4.4
5  25  4.4
6  20  4.0
8  20  4.0
9  20  4.0
10  20  4.1
12  20  4.1
14  20  4.1
15  19  4.0
16  19  4.0

Table 1: Metaring MAC Protocol

where the blocks A_{i} have dimension 16; Problems 2 models a teletraffic system as a continuous time QBD Markov chain, where the blocks A_{i} have dimension 24.

For Problem 1 we have approximated the eigenvalues of the matrix G by means of Aberth’s method, applied to the function (13), and the left eigenvectors by solving the linear systems (14). The time needed to approximate all the eigenvalues, and the corresponding left eigenvectors, was 2.5 seconds, and 0.02 seconds, respectively. We have sorted the eigenvalues in such a way that |λ_{1}| ≥ |λ_{2}| ≥ … ≥ |λ_{k}|, and we have chosen different subsets Λ_{h} = {λ_{1}, …, λ_{h}}, where h is such that properties 1, 2, and 3 of Section 2 are satisfied. Table 1 reports: the number h of eigenvalues shared by X_{0} with G; the number of iterations needed by the functional iteration method to reach the required accuracy; the total time (in seconds) needed to compute the matrix X_{0} and to perform the functional iterations. Figure 1 reports the logarithm (to the base 10) of the residual error for h = 1, 2, 16.

For this problem, the time needed to reach the required accuracy with the faster functional iteration X_{n+1} = A_{0}(I − ∑_{i=1}^{+∞} X_{n}^{i−1}A_{i})^{−1} [14] and X_{0} = I is 207.3 seconds, with 2943 iterations. It is worth pointing out the substantial reduction of the number of iterations (and of the total time) when h passes from 1 to 2, and the strong advantage with respect to the faster functional iteration method starting with the identity matrix.

For Problem 2 we have computed the eigenvalues of the matrix G, by approximating the eigenvalues of an auxiliary 2k × 2k matrix, as suggested in [6], and the eigenvectors, by solving the linear systems (14). The time needed was 0.01 seconds. The eigenvalues, which are real and distinct for this problem, are

Figure 1: Metaring MAC Protocol

sorted in such a way that |λ_{1}| ≥ |λ_{2}| ≥ … ≥ |λ_{k}|. We have chosen X_{0} with eigenvalues λ_{1}, …, λ_{h}, and left eigenvectors w_{1}^{t}, …, w_{h}^{t}, for different values of h in {1, …, k}. Table 2 reports: the number h of eigenvalues shared by X_{0} with G; the number of iterations needed by the functional iteration method to reach the required accuracy; the total time (in seconds) needed to compute the matrix X_{0} and to perform the functional iterations.

From Table 2 and Figure 2, we can observe a substantial reduction of the number of iterations for increasing values of h. When h = 24, i.e., when all the eigenvalues and the corresponding left eigenvectors are known, the matrix X_{0} is a poor approximation of G, that is refined in few steps of functional iteration.

References

[1] O. Aberth. Iteration methods for finding all zeros of a polynomial simultaneously. Math. Comp., 27, 1973.

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

h Iterations Time (s.)
1  83361  167.1
2  67271  134.1
5  39732  79.8
10  16718  33.5
12  11198  24.0
15  4736  9.9
18  340  0.8
24  340  0.8

Table 2: PH/PH/D Queue Figure 2: PH/PH/D Queue [3] 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, 1995.

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

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

[6] J. N. Daigle and D. M. Lucantoni. Queueing systems having phasedependent arrival and service rates. In W. J. Stewart, editor, Numerical solution of Markov Chains, pages 161–202. Marcel Dekker, New York, 1991.

[7] L. W. Ehrilch. A modified Newton method for polynomials. Comm. A.C.M., 10:107–108, 1967.

[8] P. Favati and B. Meini. Relaxed functional iteration techniques for the numerical solution of M/G/1 type Markov chains. BIT, 38, 1998. In printing.

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

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

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

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

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

[14] B. Meini. New convergence results on functional iteration techniques for the numerical solution of M/G/1 type Markov chains. Numer. Math., 78:39–58, 1997.

[15] G. R. Murthy, M. Kim, and E. J. Coyle. Equilibrium analysis of skip free Markov chains: nonlinear matrix equations. Commun. Statist. Stochastic Models, 7:547–571, 1991.

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

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

[18] R. S. Varga. Matrix Iterative Analysis. Prentice Hall, New Jersey, 1962.