Relaxed functional iteration techniques for the numerical solution of M/G/1 type Markov chains

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

Abstract We introduce a new iterative method for the computation of the minimal nonnegative solution G of the matrix equation X = ∑_{i=0}^{+∞} X^{i}A_{i}, arising in the numerical solution of M/G/1 type Markov chains. The idea consists in applying a relaxation technique to customarily used functional iteration formulas. The proposed method is easy to implement and outperforms, in terms of number of iterations and execution time, the standard functional iteration techniques.

AMS subject classification. 65U05, 65F15, 15A51, 15A24, 60J10, 60K25. Key Words. Functional iterations, relaxation, Markov chains, M/G/1 type matrices.

1 Introduction.

Let P be the infinite column stochastic matrix

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

defined by the k × k blocks B_{i+1}, A_{i}, i ≥ 0. A nonnegative matrix M (denoted with M ≥ O), possibly infinite, is called column stochastic if e^{T} M = e^{T}, where e is the vector having all the entries equal to 1. Matrices of the structure (1) are known in literature as stochastic matrices of M/G/1 type [18] and arise in a wide variety of queueing problems modeled by Markov chains M, where P^{T} is the transition matrix associated with M. One of the major problems related to Markov chains is the computation of the nonnegative vector π such that

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

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

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

†Dip. Matematica, via Buonarroti 2, 56127 Pisa, Italy. E-mail: meini@dm.unipi.it

If system (2) has a unique solution, P is called positive recurrent and π is called the probability invariant vector associated with P. In the case where P has the structure (1) the computation of π can be reduced (compare [18]) to the computation of the minimal nonnegative solution G of the nonlinear matrix equation

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

where X is a k × k matrix. In the case where the matrix P is irreducible and positive recurrent [8], equation (3) has a unique nonnegative solution, which is column stochastic and has eigenvalues λ_{1}, …, λ_{k}, with |λ_{i}| < 1, i = 1, …, k −1, and λ_{k} = 1 [18].

Once the matrix G is known, an arbitrary number of components of the vector π can be recovered by means of a recursive numerically stable formula, called Ramaswami’s formula [20], which involves the block entries A_{i}, B_{i} of the matrix P in (1), and G.

The problem of the computation of the nonnegative solution G of (3) has been investigated by several authors. New numerical methods have been recently proposed in [5], [6], [13], [15], [16], [2]: in [5], [6], [15] a method based on cyclic reduction is presented; in [13] a modification of Newton’s method is devised; in [16] the evaluation of the matrix G is reduced to solving, by means of a doubling technique, a block Hessenberg system; in [2] a suitable invariant subspace is computed in order to approximate G. All these methods are quadratically convergent, however some of them are not easy to implement ([5], [6], [16]), or apply only to specific problems ([2], [15]), or have a high computational cost for iteration ([13], [2], [16]). The classical methods [14], [19], customarily used in the applications and easy to implement, are functional iterations based on the recursion

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

where X_{0} is a nonnegative matrix. The matrix function F (·) belongs to the class F of matrix functions such that:

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

where H(X) = ∑_{i=0}^{+∞} X^{i}H_{i}, K(X) = ∑_{i=0}^{+∞} X^{i}K_{i} verify H_{i} ≥ O, K_{i} ≥ O, i ≥ 0, and the formal power series A(z) = ∑_{i=0}^{+∞} z^{i}A_{i} is such that A(z) = z ∑_{i=0}^{+∞} z^{i}H_{i} + ∑_{i=0}^{+∞} z^{i}K_{i}.

The sequence {X_{n}} of (4), with F (·) belonging to F, converges linearly to the matrix G. A detailed analysis of the speed of convergence, in relation with the particular function F (·) and with the starting matrix X_{0}, has been provided in [10], [17]. In particular, it is shown that:

– setting X_{0} = Q, with Q any stochastic matrix, leads to a sequence (4) that converges more rapidly than the sequence obtained by starting with X_{0} = 0; – the functional iteration method generated by choosing

F (X) = A_{0} (I − ∑_{i=1}^{+∞} X^{i−1}A_{i})^{−1}  (6)

is the fastest one in the class of method defined by F ∈ F, in the case where X_{0} = 0.

In this paper we consider functional iteration methods based on the recursion

X_{n+1} = (1 − ω)X_{n} + ωF (X_{n}), n ≥ 0  (7)

where X_{0} is a stochastic matrix, F (·) is given by (5) and ω is a relaxation parameter. In this way we propose a new method for the computation of the matrix G, based on a relaxation technique applied to customary functional iterations, which is still easy to implement and has a higher convergence rate. More precisely, by performing a convergence analysis of the sequence (7), we prove that the asymptotic rate of convergence is bounded above by the spectral radius of a suitable matrix W_{ω}. We then characterize, in terms of the matrix G, the parameter ω^{∗} which minimizes the spectral radius of the matrix W_{ω}. This theoretical result, even if interesting in itself, does not allow us to compute the asymptotic optimal relaxation parameter ω^{∗} without knowing the matrix G. We overcome this problem by dynamically computing, at the general n–th step, an optimal relaxation parameter ω_{n} and generating the sequence of matrices

X_{n+1} = (1 − ω_{n})X_{n} + ω_{n}F (X_{n}), n ≥ 0.  (8)

At each step, the parameter ω_{n} minimizes the spectral radius of a suitable matrix W_{n,ω} that relates the error of the approximation X_{n+1} to the error of the approximation X_{n}. The parameter ω_{n} can be evaluated by only knowing the approximation X_{n} and the eigenvalues of the matrix G, which can be approximated by numerically computing the zeros of the function det(zI − ∑_{i=0}^{+∞} A_{i}z^{i}) lying in the closed unit circle of the complex plane.

Our algorithm has been implemented in Fortran 77 and tested with several input data. From the many numerical experiments that we performed, we observed that the number of iterations required by the relaxed functional iteration formulas is about half of the number of iterations required by the standard algorithm, moreover, the computational cost needed to estimate the relaxation parameters ω_{n} is negligible with respect to the overall cost.

The paper is organized as follows. In Section 2 we analyze the convergence properties of the relaxed functional iteration formulas (7) and we characterize the asymptotic optimal relaxation parameter ω^{∗}. In Section 3 we propose the dynamical relaxation technique of (8) and we give some methods to approximate the eigenvalues of the matrix G. In Section 4 we present the results of the performed numerical experiments.

Here and hereafter, Re x and Im x denote the real and imaginary part of the complex number x, ρ(A) denotes the spectral radius of the matrix A, vec(A) denotes the function that associates with the matrix A the vector obtained by column-wise arranging its entries, vec^{−1}(·) denotes the inverse function of vec(·) and A ⊗ B denotes the Kronecker product of the matrices A and B, i.e., A ⊗ B is the block matrix defined by the blocks a_{i,j}B, where A = {a_{i,j}}_{i,j}.

We also introduce the vector norm ∥v∥_{1} = ∑_{i} |v_{i}|, for vectors v = {v_{i}}_{i} and the corresponding induced matrix norm ∥A∥_{1} = max_{j} ∑_{i} |a_{i,j}|, for matrices A = {a_{i,j}}_{i,j}.

Throughout, we suppose that the matrix P is irreducible and positive recurrent, and that the matrix ∑_{i=0}^{+∞} A_{i} is irreducible.

2 Convergence properties of relaxed functional iteration methods.

In this section we analyze the convergence properties of relaxed functional iteration method defined by (7), where X_{0} is a column stochastic matrix and F (X) is given by (5). In particular, we characterize the optimal relaxation parameter ω^{∗}, such that the sequence (7) has the lowest asymptotic convergence bound. The convergence proof is based on the following general result:

Theorem 1 Let Ω ⊂ R^{m} an open set, g : Ω → R^{m} a C^{1}(Ω) function and α ∈ Ω a fixed point of g(x). Denote by J(x) the m × m Jacobian matrix associated with g(x). Let S_{1} be an m_{1}–dimensional linear space and S_{2} the linear space S_{2} = S_{1}^{⊥}, generated by the orthonormal columns of the matrices S_{1} ∈ R^{m×m_{1}} and S_{2} ∈ R^{m×m_{2}}, respectively, where m_{2} = m − m_{1}. Consider the sequence

x_{n+1} = g(x_{n}), n ≥ 0  (9)

and suppose that:

i) (x_{0} − α) ∈ S_{2} ⇒ (x_{n} − α) ∈ S_{2} for any n

ii) w^{T} J(α) ∈ S_{1} for any w ∈ S_{1}.

If ρ(S_{2}^{T} J(α)S_{2}) < 1 then there exists a neighborhood V ⊂ Ω of α such that the sequences x_{n} of (9) converge to α for any x_{0} ∈ V such that x_{0} − α ∈ S_{2}.

Proof. Let ϵ be such that 0 < ϵ < (1 − ρ(S_{2}^{T} J(α)S_{2}))/2 and let T an m_{2} × m_{2} nonsingular matrix such that T S_{2}^{T} J(α)S_{2}T^{−1} is the Jordan canonical form of S_{2}^{T} J(α)S_{2}. Then, it is well known [21] that the vector norm on R^{m_{2}}

||x||_{ϵ,m_{2}} = ||D_{ϵ}^{−1}T x||_{∞},

where D_{ϵ} = Diag(1, ϵ, ϵ^{2}, …, ϵ^{m_{2}−1}), is such that the corresponding induced matrix norm satisfies

∥S_{2}^{T} J(α)S_{2}∥_{ϵ,m_{2}} ≤ ρ(S_{2}^{T} J(α)S_{2}) + ϵ.

Denote by M the m × m matrix

M = ( \begin{matrix}I & O \\ O & D_{ϵ}^{−1}T \end{matrix} )

and define the vector norms on R^{m}:

\begin{matrix}||x||_{ϵ,m} = ||Mx||_{∞} \\ ||x||_{S} = ||S^{T} x||_{ϵ,m}, \end{matrix}

where S = (S_{1}|S_{2}). Since g ∈ C^{1}(Ω) then there exists a neighborhood V ⊂ Ω of α such that

∥g(x) − g(α) − J(α)(x − α)∥_{S} ≤ ϵ∥x − α∥_{S}, ∀x ∈ V.

We prove that, if x_{0} ∈ V and (x_{0} − α) ∈ S_{2}, then, for each n, x_{n} ∈ V and ||x_{n+1} − α||_{S} ≤ σ||x_{n} − α||_{S}, for a suitable 0 < σ < 1. This inequality can be proved by observing that:

∥x_{n+1} − α∥_{S} = ∥S^{T} (x_{n+1} − α)∥_{ϵ,m} ≤
∥S^{T} ((x_{n+1} − α) − J(α)(x_{n} − α))∥_{ϵ,m} + ∥S^{T} J(α)SS^{T} (x_{n} − α)∥_{ϵ,m};

on the other hand, from hypothesis i) and ii), it can be easily verified that

S^{T} J(α)SS^{T} (x_{n} − α) = ( \begin{matrix}O \\ S_{2}^{T} J(α)S_{2}S_{2}^{T} (x_{n} − α) \end{matrix} ) ,

whence

∥x_{n+1} − α∥_{S} ≤ ϵ∥x_{n} − α∥_{S} + ||S_{2}^{T} J(α)S_{2}S_{2}^{T} (x_{n} − α)||_{ϵ,m_{2}} ≤

ϵ∥x_{n}−α∥_{S}+∥S_{2}^{T} J(α)S_{2}∥_{ϵ,m_{2}}∥S_{2}^{T} (x_{n}−α)∥_{ϵ,m_{2}} ≤ (2ϵ+ρ(S_{2}^{T} J(α)S_{2}))∥x_{n}−α∥_{S}.

For the choice of ϵ, it is σ = 2ϵ + ρ(S_{2}^{T} J(α)S_{2}) < 1. □

Let F_{ω}(X) = (1 − ω)X + ωF (X) and let g_{F,ω} : R^{k^{2}} → R^{k^{2}} such that

g_{F,ω}(x) = vec(F_{ω}(vec^{−1}(x)))  (10)

In order to apply Theorem 1 to our problem the following lemma is needed:

Lemma 2 Denote by J_{F,ω}(x) the k^{2} × k^{2} Jacobian matrix associated with the function g_{F,ω} defined in (10) and let g = vec(G). Then J_{F,ω}(g) = (1 − ω)I + ωR̂_{F} where

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

and

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

Proof. The proof directly follows from the properties of Kronecker product, see [12]. □

Define, for each integer n ≥ 0, the matrices E_{n} = G − X_{n}, which represent the error at step n of the sequence {X_{n}} of (7). Let ê_{n} = vec(E_{n}) be the associated k^{2}-dimensional error vector and let S_{1} be the linear space generated by the orthogonal vectors

v_{1} = \frac{1}{\sqrt{k}} (e_{1} ⊗ e), v_{2} = \frac{1}{\sqrt{k}} (e_{2} ⊗ e), …, v_{k} = \frac{1}{\sqrt{k}} (e_{k} ⊗ e),  (12)

where e is the k-dimensional vector having all the entries equal to 1 and e_{i} is the k-dimensional vector having the i-th entry equal to 1 and the remaining ones equal to zero. It can be easily observed that, since the column sum of the matrices X_{n} and G is equal to 1, the error vector ê_{n} belongs to the linear space S_{1}^{⊥}, for any n. Moreover, by taking into account equation (11), it is easy to verify that w^{T} J_{F,ω}(g) ∈ S_{1} for any w ∈ S_{1} . So the following convergence result is a straightforward consequence of Theorem 1.

Theorem 3 Let F_{ω}(X) and J_{F,ω}(x) be defined as in Lemma 2. Let S_{1} the linear space spanned by vectors of (12) and let S_{2} the matrix whose columns are an orthonormal basis of the linear space S_{2} = S_{1}^{⊥}. If ω is such that ρ(S_{2}^{T} J_{F,ω}(g)S_{2}) < 1 then there exists a neighborhood V ⊂ R^{k×k} of G such that the sequences X_{n+1} = F_{ω}(X_{n}) converge to G for any column stochastic matrix X_{0} ∈ V.

With the same arguments of [17], we can prove the following general recursive relations for the error vectors ê_{n}:

Theorem 4 It holds that

ê_{n+1} = ((1 − ω)I + ωR̂_{n})ê_{n}, n ≥ 0  (13)

where

\begin{matrix}R̂_{n} = ∑_{i=1}^{+∞} (X_{n}^{i−1}K_{i}(I − H(X_{n}))^{−1})^{T} ⊗ I+ \\ ∑_{j=2}^{+∞} ∑_{i=j}^{+∞} (X_{n}^{i−j}A_{i}(I − H(X_{n}))^{−1})^{T} ⊗ G^{j−1}. \end{matrix}

Observe that the matrix R̂_{n}, for n ≥ 0, can be expressed in the form

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

where

\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}  (15)

Moreover, under the hypothesis that sequence (7) is convergent, the following theorem gives an upper bound to the rate of convergence:

Theorem 5 Let

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

be the mean asymptotic rate of convergence of the sequence (7) where F (·) is defined by (5) and ∥ · ∥ is any vector norm. Then

r ≤ ρ(S_{2}^{T} J_{F,ω}(g)S_{2}) = ρ((1 − ω)I + ωW_{F})

where W_{F} = S_{2}^{T} R̂_{F} S_{2}.

Proof. The proof is analogous to that of Theorem 15 in [17]

□

Hereafter we use the notations

W_{ω} = (1 − ω)I + ωW_{F}

and

f_{μ}(x) = \frac{2(Re x − Re μ)}{2(Re x − Re μ) + |μ|^{2} − |x|^{2}} .

In the following theorem we characterize the value ω^{∗} for which the upper bound to the asymptotic rate of convergence is minimum :

Theorem 6 Let X_{n} = (1 − ω)X_{n−1} + ωF (X_{n−1}), where X_{0} is a stochastic matrix and F (·) is defined by (5). Let Λ be the set of the eigenvalues of the matrix W_{F} and Γ = {α ∈ Λ : Re α = max_{λ∈Λ} Re λ}. Then the parameter ω^{∗} such that ρ(W_{ω^{∗}}) = min_{ω>0} ρ(W_{ω}) is given by

ω^{∗} = min \{f_{μ}(ν), \frac{1 − Re μ}{1 − 2Re μ + |μ|^{2}} \}  (16)

where:

- μ ∈ Γ is such that Im μ = max_{α∈Γ} Im α - ν ∈ Λ is such that f_{μ}(ν) = min_{λ∈Λ−Γ} f_{μ}(λ).

Proof. Given the matrix W_{F} the value

ρ^{2}(W_{ω}) = max_{λ∈Λ} p_{λ}(ω) = max_{λ∈Λ} [ω^{2}(1 + |λ|^{2} − 2Re λ) − 2ω(1 − Re λ) + 1]

as a function of ω is the maximum of a sheaf of k^{2} − k parabolas with center in (0, 1) and vertices in

(x_{λ}, y_{λ}) = (\frac{1 − Re λ}{1 + |λ|^{2} − 2Re λ}, \frac{Im ^{2}λ}{1 + |λ|^{2} − 2Re λ}) , λ ∈ Λ.

Since ρ(W_{F}) < 1 (compare [17]), for any λ ∈ Λ we have Re λ < 1, whence x_{λ} > 0 and 0 < y_{λ} < 1. Moreover we have p_{λ}(0) = 1 and p′_{λ}(0) = 2(Re λ − 1). Hence, for ω in a right neighborhood of zero, the value max_{λ∈Λ} p_{λ}(ω) is reached for μ ∈ Γ such that Im μ = max_{α∈Γ} Im α. It can be observed that

min_{ω>0} (max_{λ∈Λ} p_{λ}(ω)) = p_{μ}(ω^{∗})

where

ω^{∗} = min (x_{μ}, ω̂),

with

ω̂ = min {ω > 0 : p_{λ}(ω) = p_{μ}(ω), λ ≠ μ}

that is

ω̂ = min_{λ∈Λ−Γ} f_{μ}(λ).

□

3 A dynamical relaxation technique.

The asymptotic optimal relaxation parameter ω^{∗} characterized by Theorem 6 cannot be evaluated without knowing the matrix G. We overcome this problem by dynamically computing an optimal relaxation parameter by minimizing, at each step, the spectral radius of a suitable matrix which relates the error of two successive approximations. More precisely, since the error vector ê_{n} belongs to the linear space S_{2} for any n, the recurrences (13) can be rewritten (compare [17]), by premultiplying both sides by the matrix S^{T}, in the form

S_{2}^{T} ê_{n+1} = ((1 − ω)I + ωW_{n}))(S_{2}^{T} ê_{n}),  (17)

where W_{n} = S_{2}^{T} R̂_{n}S_{2}. Based on these considerations, we generate the sequence of matrices

X_{n+1} = (1 − ω_{n})X_{n} + ω_{n}F (X_{n}), n ≥ 0,  (18)

where the value of ω_{n} is such that

ρ(W_{n,ω_{n}}) = min_{ω>0} ρ(W_{n,ω}),

and W_{n,ω} = (1 − ω)I + ωW_{n}. Indeed, the following results hold:

Theorem 7 Let X_{n} be the approximation of the matrix G obtained at step n. Let Λ_{n} be the set of the eigenvalues of the matrix W_{n} and Γ_{n} = {α ∈ Λ_{n} : Re α = max_{λ∈Λ_{n}} Re λ}. If ρ(W_{n}) < 1, then the parameter ω_{n} such that

ρ(W_{n,ω_{n}}) = min_{ω>0} ρ(W_{n,ω})

is given by

ω_{n} = min \{f_{μ}(ν), \frac{1 − Re μ}{1 − 2Re μ + |μ|^{2}} \}

where:

- μ ∈ Γ_{n} is such that Im μ = max_{α∈Γ_{n}} Im α

- ν ∈ Λ_{n} is such that f_{μ}(ν) = min_{λ∈Λ_{n}−Γ_{n}} f_{μ}(λ).

Proof. The proof can be carried out analogously to the proof of Theorem 6. □

In the case where ρ(W_{n}) ≥ 1 we set ω_{n} = 1, i.e. we apply the standard functional iteration formulas; however, since the non relaxed sequence converges to the matrix G, after a finite number of steps the condition ρ(W_{n}) < 1 is satisfied.

Observe that the optimal parameter ω_{n} depends on the eigenvalues of the matrix W_{n}, that can be computed by only knowing the approximation X_{n} and the eigenvalues of the matrix G, as stated by the following proposition:

Proposition 8 Let λ_{1}, λ_{2}, …, λ_{k} be the eigenvalues of the matrix G, with λ_{k} = 1. For i = 1, …, k, set

Q_{i,n} = ∑_{j=0}^{+∞} λ_{i}^{j}Y_{j,n},

where the matrices Y_{j,n} are defined in (15). Then the set Λ_{n} of the eigenvalues of the matrix W_{n} is given by

Λ_{n} = ⋃_{i=1,...,k−1}{μ : μ is eigenvalue of Q_{i,n}}

Proof. Let us first consider the eigenvalues of the matrix R̂_{n} of (14). Observe that the matrix R̂_{n} is similar to the matrix R_{n} = ∑_{j=0}^{+∞} G^{j} ⊗ Y_{j,n}. Let T = U^{H}GU be the Schur canonical form of the matrix G, where U is a unitary matrix and T is an upper triangular matrix whose diagonal entries are the eigenvalues λ_{1}, …, λ_{k} of G. Then, we have

(U ⊗ I)^{−1}R_{n}(U ⊗ I) = ∑_{j=0}^{+∞} T^{j} ⊗ Y_{j,n},

that is, the matrix R_{n} is similar to a block upper triangular matrix, whose diagonal blocks are the matrices Q_{i,n} = ∑_{j=0}^{+∞} λ_{i}^{j}Y_{j,n}. Whence the eigenvalues of the matrix R̂_{n} are the eigenvalues of the matrices Q_{i,n}, for i = 1, …, k. On the other hand the matrix R̂_{n} is also similar to the matrix

S^{T} R̂_{n}S = ( \begin{matrix}V_{n} & O \\ T_{n} & W_{n} \end{matrix} )

where the matrix V_{n} coincides with the matrix Q_{k,n}, corresponding to λ_{k} = 1. Therefore the eigenvalues of the matrix W_{n} are given by the set Λ_{n}. □

Propositions 7 and 8 give a complete characterization of the parameter ω_{n}, by means of the eigenvalues of the matrix G only. Moreover, such eigenvalues can be simply computed without calculating the matrix G. Indeed, it is well known [11] that the eigenvalues of the matrix G are given by the zeros of the function

φ(z) = det(zI − ∑_{i=0}^{+∞} A_{i}z^{i}) = det Φ(z)  (19)

lying in the closed unit circle of the complex plane.

In the case where A_{i} = 0 for i > 2 (QBD, Quasi–Birth–Death, case), the eigenvalues of the matrix G can be simply computed by approximating the eigenvalues of a suitable 2k × 2k matrix, by means of a standard routine, with the cost of O(k^{3}) arithmetic operations. Indeed, as shown in [9], the eigenvalues of G are given by the reciprocal of the eigenvalues of the matrix

( \begin{matrix}−(A_{1} − I)A_{0}^{−1} & −A_{2}A_{0}^{−1} \\ I & 0 \end{matrix} )  (20)

having modulus greater than or equal to 1. In the general case, where the blocks A_{i} are negligible only for large values of i, in order to compute the zeros of (19), we can apply a generalization of the Aberth method [1], which has been originally introduced for the simultaneous computation of the zeros of a polynomial. Aberth’s method, applied to compute the zeros of the polynomial P (z) of degree N, consists in generating the sequence of the N approximations of the zeros

z_{i}^{(n+1)} = z_{i}^{(n)} − \frac{1}{\frac{P′(z_{i}^{(n)})}{P (z_{i}^{(n)})} − ∑_{{j=1 \atop j≠i}}^{N} \frac{1}{z_{i}^{(n)} − z_{j}^{(n)}}} , i = 1, …, N,  (21)

where z_{1}^{(0)}, …, z_{N}^{(0)} are the initial approximations. If the zeros are pairwise distinct, the method has a local cubic convergence. If, instead of a polynomial, we take an analytic function with distinct zeros, the order of convergence is still cubic. The effectiveness of Aberth’s method for the computation of the zeros of an analytic function in the closed unit circle has been already tested in [7].

In our case we have to find the zeros of the analytic function φ(z) of (19) in the closed unit circle. In the iterative formulae (21) the most complicated computation consists in evaluating \frac{φ′(z)}{φ(z)}; such computation can be simplified by means of the following:

Lemma 9 Let S(z) = ∑_{i=0}^{+∞} V_{i}z^{i} be a series having matrix coefficients V_{i}. Then

\frac{(det S(z))′}{det S(z)} = tr(S(z)^{−1}S′(z))

where tr(M) denotes the trace of the matrix M and S′(z) denotes the derivative of the series S(z).

Proof. Denote by s_{i,j}(z) the entries of the matrices S(z). Observe that

(det S(z))′ = ∑_{i=1}^{k} ∑_{j=1}^{k} \frac{∂ det S}{∂s_{i,j}} \frac{ds_{i,j}(z)}{dz},

where k is the dimension of the matrices V_{i}. On the other hand, for any h ∈ {1, …, k}, we have

det S(z) = ∑_{l=1}^{k}(−1)^{h+l}s_{h,l}(z)[S(z)]_{h,l},

where [S(z)]_{h,l} denotes the determinant of the (k − 1) × (k − 1) matrix obtained by deleting the h-th row and the l-th column from the matrix S(z). Therefore

\frac{∂ det S}{∂s_{i,j}} = (−1)^{i+j}[S(z)]_{i,j},

for any i, j. Whence

(det S(z))′ = ∑_{i=1}^{k} ∑_{j=1}^{k}(−1)^{i+j}[S(z)]_{i,j} \frac{ds_{i,j}(z)}{dz} =

det S(z) ∑_{i=1}^{k} ∑_{j=1}^{k}(S(z)^{−1})_{j,i} \frac{ds_{i,j}(z)}{dz} = det S(z)tr(S(z)^{−1}S′(z)).

□

In practice, the series ∑_{i=0}^{+∞} A_{i}z^{i} is truncated at the index m such that the matrix ∑_{i=0}^{m} A_{i} is numerically stochastic. A nonnegative matrix A is said numerically stochastic if e^{T} − e^{T} A < τe^{T}, where τ is a lower bound to the machine precision. Based on (21) and Lemma 9, the following algorithm for the computation of the eigenvalues of the matrix G can be carried out:

Algorithm 3.1 Computation of the eigenvalues of the matrix G by means of Aberth’s method.

Input. Nonnegative k × k matrices A_{i}, i = 0, 1, …, m, such that ∑_{i=0}^{m} A_{i} is numerically stochastic, starting approximations λ_{i}^{(0)}, i = 1, …, k − 1 of the eigenvalues, and a real number δ > 0.

Output. An approximation of the eigenvalues λ_{1}, …, λ_{k} of the matrix G.

Computation.

1. Set n = 0 and λ_{k}^{(j)} = 1, for any nonnegative integer j.

2. Repeat

(a) Compute, for i = 1, …, k − 1,

λ_{i}^{(n+1)} = λ_{i}^{(n)} − \frac{1}{tr(Φ(λ_{i}^{(n)})^{−1}Φ′(λ_{i}^{(n)})) − ∑_{{j=1 \atop j≠i}}^{k} \frac{1}{λ_{i}^{(n)} − λ_{j}^{(n)}}} .

(b) Set n = n + 1.

Until

The conditioning of the matrix Φ(λ_{i}^{(n)}) is greater than δ, for i = 1, …, k − 1.

The most expensive computations of each step n in the above algorithm are the evaluation of the matrix polynomials Φ(λ_{i}^{(n)}), Φ′(λ_{i}^{(n)}) and the inversion of the matrices Φ(λ_{i}^{(n)}), for i = 1, …, k − 1. Whence, it can be easily verified that the computational cost of each step is O(k^{4} + mk^{3}) arithmetic operations.

Once the eigenvalues of the matrix G are available, based on Propositions 7 and 8, the following algorithm for the computation of the matrix G by means of a relaxation technique can be carried out:

Algorithm 3.2 Computation of the matrix G.

Input. Nonnegative k × k matrices A_{i}, i = 0, 1, …, m, such that ∑_{i=0}^{m} A_{i} is numerically stochastic and an error bound ϵ > 0.

Output. The matrix G approximating the solution of the equation (3).

Computation.

1. Compute the eigenvalues of the matrix G.

2. Set X_{0} = I and n = 0.

3. Repeat (a) Compute the eigenvalues of the matrix W_{n} and:

i. if ρ(W_{n}) ≥ 1 then set ω_{n} = 1; ii. if ρ(W_{n}) < 1 then compute the optimal relaxation parameter

ω_{n} of Proposition 7.

(b) Compute X_{n+1} = (1 − ω_{n})X_{n} + ω_{n}F (X_{n}).

(c) Set n = n + 1.

Until

||X_{n} − X_{n−1}||_{1} < ϵ

The computational cost to evaluate the optimal relaxation parameter at each step n, i.e., to perform stage (a) of the above algorithm, is O(k^{4} + mk^{3}) arithmetic operations. Moreover, the cost to compute the new approximation of the solution G, i.e., to perform stage (b), is O(mk^{3}) arithmetic operations. Therefore the overall computational cost of each step of the above algorithm is O(k^{4} + mk^{3}) arithmetic operations. Observe that the most expensive stage consists in evaluating, at each step n, the optimal relaxation parameter ω_{n}. On the other hand, from the many numerical experiments that we performed, we observed that the sequence ω_{n} rapidly converges to the asymptotic optimal relaxation parameter ω^{∗} of Theorem 6. This suggested us to stop the computation of the relaxation parameter ω_{n} and set ω_{n} = ωn̄, for n > n̄, where

|ωn̄ − ω_{n̄−1}| < σ_{1}, ||X_{n̄+1} − Xn̄||_{1} < σ_{2}

for suitable positive σ_{1}, σ_{2}. In this way we arrive at the following modification of Algorithm 3.2:

Algorithm 3.3 Computation of the matrix G.

Input. Nonnegative k × k matrices A_{i}, i = 0, 1, …, m, such that ∑_{i=0}^{m} A_{i} is numerically stochastic and error bounds ϵ, σ_{1}, σ_{2} > 0.

Output. The matrix G approximating the solution of the equation (3).

Computation.

1. Compute the eigenvalues of the matrix G.

2. Set X_{0} = I and n = 0.

3. Repeat

(a) If n ≥ 2, |ω_{n−1} − ω_{n−2}| < σ_{1} and ||X_{n} − X_{n−1}||_{1} < σ_{2} then set ω_{n} = ω_{n−1}; otherwise, compute the eigenvalues of the matrix W_{n} and i. if ρ(W_{n}) ≥ 1 then set ω_{n} = 1; ii. if ρ(W_{n}) < 1 then compute the optimal relaxation parameter

ω_{n} of Proposition 7.

(b) Compute X_{n+1} = (1 − ω_{n})X_{n} + ω_{n}F (X_{n}).

(c) Set n = n + 1.

Until

||X_{n} − X_{n−1}||_{1} < ϵ

In the above algorithm, in order to estimate the optimal relaxation parameter ω_{n} at step n, we perform only few times the computation of the eigenvalues of the matrix W_{n}. In this way we have the same asymptotic computational cost per step as the standard functional iterations, i.e. O(mk^{3}) operations; at the same time the rate of convergence is substantially improved.

We refer the reader to next section for a discussion of the numerical results.

4 Numerical results.

We have performed numerical experimentations of Algorithm 3.3 by choosing as function F (·) the one defined by formula (6) and as starting approximation X_{0} = I. In this case, in order to compute at each step the matrix Z_{1} = ∑_{i=1}^{m} X_{n}^{i−1}A_{i}, we applied the Horner scheme, thus computing the matrices:

\begin{matrix}Z_{m} = A_{m} \\ Z_{j} = X_{n}Z_{j+1} + A_{j}, for j = m − 1, …, 1. \end{matrix}  (22)

In the case where we need the evaluation of the eigenvalues of the matrix W_{n}, i.e., of the matrices Q_{i,n} = ∑_{j=0}^{+∞} λ_{i}^{j}Y_{j,n}, for i = 1, …, k − 1, (compare Proposition 8), the Horner scheme (22) allows us to compute the matrices Y_{j,n} by means of the relations:

Y_{0,n} = 0
Y_{j,n}^{T} = Z_{j+1}(I − Z_{1})^{−1}, j = 1, …, m − 1
Y_{j,n} = 0, j ≥ m.

We have implemented the algorithms in Fortran 77 and the programs have been run on a Sparc Workstation using the standard double precision IEEE arithmetic.

We have tested the dynamical relaxation technique on two families of M/G/1 type problems performing comparisons with the standard functional iteration techniques. For Algorithm 3.3 we have chosen ϵ = 10^{−12}, σ_{1} = σ_{2} = 10^{−4}; for standard functional iteration technique, we have used the same stopping criterium of Algorithm 3.3 based on the difference between two successive approximations and the same error bound ϵ = 10^{−12}.

The first family of problems is taken from [9], where the authors model a teletraffic system as a continuous-time QBD Markov chain, defined by the 24×24 blocks A′_{0}, A′_{1}, A′_{2}. The blocks A′_{0} and A′_{2} are diagonal, with A′_{0} = 192ρ_{d}I, (A′_{2})_{j,j} = 192(1 − j/24), j = 0, …, 23. The matrix A′_{1} is tridiagonal with (A′_{1})_{j,j+1} = ar(M − j)/M, j = 0, …, 22, and (A′_{1})_{j,j−1} = jr, j = 1, …, 23,

Relaxed FIF  Standard FIF  Standard/Relaxed FIF
M  Steps Time (s.) ω_{∗}  Steps Time (s.) Steps  Time
64 31792  1107  2.086 63776  1832  2.01  1.65
128 46081  1571  2.086 92209  2620  2.00  1.67
256 55402  1850  2.086 110728  3165  2.00  1.71
512 60650  2024  2.086 121097  3440  2.00  1.70
1024 63404  2097  2.086 126593  3603  2.00  1.72
2048 64819  2139  2.086 129342  3690  2.00  1.72

Figure 1: Teletraffic system

where the quantities ρ_{d}, a, r and M are parameters of the system. It is well known that the matrix G solves the matrix equation X = A_{0} + XA_{1} + X^{2}A_{2}, where A_{0}^{T} = −A′_{1}^{−1}A′_{2}, A_{1} = 0 and A_{2}^{T} = −A′_{1}^{−1}A′_{0}. By following the ideas of [9], we have computed the eigenvalues of the matrix G by approximating the eigenvalues of the auxiliary matrix (20) by means of the routine cg.f of Eispack Library. We have applied Algorithm 3.3 and the standard functional iteration algorithm, to the above problem by setting r = 1/300, a = 18.244, ρ_{d} = 0.28 and varying M.

Figure 1 reports the execution time (in seconds) and the number of iterations needed by the two algorithms to reach the required accuracy, the ratios between the execution times and the number of iterations of the two algorithms and the asymptotic optimal relaxation parameter ω_{∗}. The time needed to compute the eigenvalues of the matrix G is not reported in the figure, since, for each value of the parameter M, it is smaller than one second. Figure 2 reports, for the value M = 64, the logarithm (to the base 10) of the error ∥X_{n} − ∑_{i=0}^{+∞} X_{n}^{i}A_{i}∥_{1} of the approximations obtained at step n by the two algorithms.

The second family of problems arises from the modeling of metropolitan queueing network (Metaring MAC Protocol [3]). For a description of the model and the structure of the blocks we refer the reader to [3]. The dimension k of the blocks A_{i} is equal to 16. We have tested Algorithm 3.3 and standard algorithm for different values of the parameter ρ, 0 < ρ < 1, which represents the stability condition of the associated Markov chain [18]: as ρ tends to one, the problem of computing the solution G of (3) by using customary techniques becomes more difficult; in fact, for ρ = 1 the Markov chain is not positive recurrent [18]. We have approximated the eigenvalues of the matrix G by means of Algorithm 3.1, based on Aberth’s method, by setting δ = 10^{12} and by randomly choosing λ_{i}^{(0)}, i = 1, …, 15, on the circle of radius 1/2. The integer m such that the matrices ∑_{i=0}^{m} A_{i} are numerically stochastic is equal to 240, for ρ = 0.8, and to 264, for the other tested values of ρ. Figure 3 reports the number of iterations and the time (in seconds) needed to approximate the eigenvalues of G by means of

Figure 2: Teletraffic system

Aberth’s method, the number of iterations and the time (in seconds) needed by the relaxed functional iteration formulas to reach the required accuracy, the asymptotic optimal relaxation parameter, and the total time. Figure 4 reports the number of iterations and the time (in seconds) needed by the standard functional iteration formulas to reach the required accuracy, and the ratios between the execution times and the number of iterations of standard and relaxed algorithms. Figure 5 reports, for the value ρ = 0.97, the logarithm (to the base 10) of the error ∥X_{n} − ∑_{i=0}^{+∞} X_{n}^{i}A_{i}∥_{1} of the approximations obtained at step n by the two algorithms.

From the above examples, we observe that the number of iterations required by the relaxed functional iteration formulas is about half of the number of iterations required by the standard algorithm, and the amount of computations needed to estimate the relaxation parameter does not significantly increase the ratio between the execution times.

The relaxation technique, applied to functional iteration formulas, thus leads to a method for the computation of the matrix G that is easy to implement and that converges faster, even if still linearly, than the standard one.

Aberth’s method  Relaxed FIF
ρ  Steps Time (s.) Steps Time(s.) ω_{∗}  Total time (s.)
0.8  13  11  88  67  1.716  78
0.9  14  13  213  153  1.801  166
0.91  16  13  251  177  1.810  190
0.92  16  13  305  210  1.819  223
0.93  15  12  389  251  1.829  263
0.94  16  13  538  348  1.838  361
0.95  16  13  861  525  1.847  538
0.96  16  13  1763  1053  1.856  1066
0.97  16  13  1848  1109  1.860  1122

Figure 3: Metaring MAC–Protocol

Standard FIF  Standard/Relaxed FIF
ρ  Steps Time (s.) Steps  Time
0.8  148  80  1.68  1.02
0.9  373  224  1.75  1.46
0.91 440  269  1.75  1.41
0.92 537  337  1.76  1.52
0.93 689  417  1.77  1.66
0.94 956  605  1.78  1.68
0.95 1534  928  1.78  1.77
0.96 3157  1911  1.79  1.79
0.97 3336  2031  1.80  1.81

Figure 4: Metaring MAC–Protocol

Figure 5: Metaring MAC–Protocol, ρ = 0.97

References

[1] O. Aberth, “Iteration methods for finding all zeros of a polynomial simultaneously”, Math. Comp, \underline{27} (1973), pp. 339–344.

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

[3] G. Anastasi, L. Lenzini, B. Meini, “Performance evaluation of a worst case model of the Metaring MAC Protocol with global fairness”, Performance Evaluation, \underline{29} (1997), pp. 127–151.

[4] A. Berman, R.J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, New York, 1979.

[5] D. Bini, B. Meini, “On the solution of a nonlinear matrix equation arising in queueing problems”, SIAM J. Matrix Anal. Appl., \underline{4} (1996), pp. 906–926.

[6] D. Bini, B. Meini, “Improved cyclic reduction for solving queueing problems”, Numer. Alg., \underline{15} (1997), pp. 57–74.

[7] D. Bini, B. Meini, “Solving certain queueing problems modeled by Toeplitz matrices”, Calcolo, \underline{30} (1993), pp. 395–420.

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

[9] J.N. Daigle, D.M. Lucantoni, “Queueing systems having phasedependent arrival and service rates”, in Numerical Solution of Markov Chains, ed. W.J. Stewart, pp. 161-202, Marcel Dekker, New York, 1991.

[10] P. Favati, B. Meini, “On functional iteration methods for solving nonlinear matrix equations arising in queueing problems”, submitted for publication, 1997.

[11] H.R. Gail, S.L. Hantler, B.A. Taylor, “Solutions of the basic matrix equation for M/G/1 and G/M/1 type Markov chains”, Commun. Statist. Stochastic Models, \underline{10} (1994), pp. 1–43.

[12] A. Graham, Kronecker Products and Matrix Calculus with Applications, Ellis Horwood Limited, Chichester, 1981.

[13] G. Latouche, “Newton’s iteration for non-linear equations in Markov chains”, IMA J. of Numerical Analysis, \underline{14} (1994), pp. 583–598.

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

[15] G. Latouche, V. Ramaswami, “A logarithmic reduction algorithm for Quasi-Birth-Death processes”, J. Appl. Prob., \underline{30} (1993), pp. 650–674.

[16] G. Latouche, G.W. Stewart, “Numerical methods for M/G/1 type queues”, in Proceedings of the Second International Workshop on Numerical Solution of Markov Chains, Raleigh, North Carolina, 1995, W.J. Stewart ed., pp. 571–581, 1995.

[17] B. Meini, “New convergence results on functional iteration techniques for the numerical solution of M/G/1 type Markov chains”, Numer. Math., \underline{78} (1997), pp. 39–58.

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

[19] V. Ramaswami, “Nonlinear matrix equations in applied probability - Solution techniques and open problems”, SIAM Review, \underline{30} (1988), pp. 256–263.

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

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