On Functional Iteration Methods for Solving M/G/1 Type Markov Chains

Paola Favati^{1} and Beatrice Meini^{2} 1Istituto di Matematica Computazionale del CNR, via S. Maria 46, 56127 Pisa, Italy; favati@imc.pi.cnr.it 2Dipartimento di Matematica, via Buonarroti 2, 56127 Pisa, Italy; meini@dm.unipi.it

Abstract: We consider functional iteration methods, based on the recursion X_{n+1} = F (X_{n}), n ≥ 0, for solving the nonlinear matrix equation X = ∑_{i=0}^{+∞} X^{i}A_{i} which arises in the numerical solution of M/G/1 type Markov chains. We propose two strategies for improving the rate of convergence of such iterative methods, based on the spectral properties of the solution G. The first strategy consists in choosing an initial approximation X_{0} which shares with G some eigenvalues and the corresponding left eigenvectors; the second one relies on a relaxation technique which modifies the spectral properties of the Jacobian matrix associated with the iteration function F. Numerical results show the effectiveness of these strategies.

1 INTRODUCTION

Many queueing problems can be modeled by Markov chains of M/G/1 type [21], 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} ) ,

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. It is well known that, in the case where P^{t} is positive recurrent, the problem of the computation of the probability invariant vector π, i.e., the solution of the infinite system π = P^{t}π, ∥π∥_{1} = 1, can be reduced (compare [21]) to the computation of the column stochastic matrix G, given by the minimal nonnegative solution of the nonlinear matrix equation

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

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

Recently, the analysis of numerical methods for computing the solution of the matrix equation (1) has been largely developed. New quadratically convergent methods have been proposed in [2], [4], [5], [6], [7], [16], [17], [18]; classical linearly convergent methods based on functional iterations have been improved and studied in details in [10], [11], [19]. 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,  (2)

where

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

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} (see [11], [19]).

It can be proved that, by suitably choosing the starting approximation X_{0}, the sequence generated by (2) linearly converges to the solution G of (1) and the rate of convergence is strongly related to the choice of X_{0} and to the function F. In the case where X_{0} = 0 in [19] it is shown that the mean asymptotic rate of convergence of {X_{n}} coincides with the spectral radius of the k^{2} × k^{2} nonnegative matrix R̂_{F}, given by the Jacobian matrix of the iteration function F computed in G. Moreover, in the case where X_{0} is any column stochastic matrix, 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 of R̂_{F} associated with the spectral radius. From these properties it follows that the mean asymptotic rate of convergence of the sequence X_{n} is less than or equal to the second largest modulus eigenvalue of R̂_{F}, if X_{0} is a column stochastic matrix and R̂_{F} is irreducible. In this way in [19] the rate of convergence is improved by suitably choosing the starting approximation X_{0}.

In this paper we provide a more precise analysis of the convergence speed of (2) and we propose two strategies to further improve the rate of convergence, based on the knowledge of the eigenvalues of the matrix G. These eigenvalues 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.

The first strategy is based on an extension of the results of [19] and consists in considering initial 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} = (1, …, 1) (shared by any column stochastic matrix).

The second strategy relies on a relaxation technique, applied in order to modify the spectral properties of the Jacobian matrix of the iteration function F.

Let λ_{1} = 1, λ_{2}, …, λ_{h} be 1 ≤ h < k 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,  (4)

that is, it shares with G the eigenvalues λ_{1}, λ_{2}, …, λ_{h} and the corresponding left eigenvectors. Then we prove that the vector ê_{n} belongs to the subspace S_{h}^{⊥}, where S_{h} is the kh-dimensional subspace generated by the vectors e_{i} ⊗ w_{j}, i = 1, …, k, j = 1, … h, and e_{i} is the k-dimensional vector having the i-th entry equal to 1, and the remaining ones equal to 0. Whence, we show that the rate of convergence of the sequence X_{n} is bounded from above by the spectral radius of the matrix obtained by projecting R̂_{F} onto the subspace T_{h} = S_{h}^{⊥}. In this way the rate of convergence can be substantially reduced. Once the eigenvalues λ_{1} = 1, λ_{2}, …, λ_{h}, and the corresponding eigenvectors w_{1}^{t}, …, w_{h}^{t} are known, a matrix X_{0} satisfying (4) can simply be obtained by

X_{0} = V Diag(λ_{1}, …, λ_{h})W^{t},  (5)

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

The idea to use the matrix X_{0} of (5) as starting point of functional iterations, is almost similar to the idea proposed in [20], [12] for the direct approximation of the matrix G. In fact, in the case where h = k, the matrix V coincides with the inverse of W^{t} and W^{−t}Diag(λ_{1}, …, λ_{k})W^{t} coincides with G. In practice, in the case where the matrix W is nearly singular, due to the roundoff errors generated in the computation, the matrix G̃ obtained in this way can be a poor approximation of G. Our approach provides the possibility of selecting those eigenvalues and eigenvectors which bring most of the information and which are not strongly affected by rounding errors, say, the eigenvectors which are not nearly linearly dependent.

In our second approach we consider functional iterations defined by the recursion

X_{n+1} = F_{ω}(X_{n}), n ≥ 0,  (6)

where F_{ω}(X) = (1 − ω)X + ωF (X) and ω is a relaxation parameter. Based on the relation between the convergence rate of the sequence (6) and the spectral properties of R̂_{F_{ω}}, we may characterize the parameter ω^{∗} which minimizes the spectral radius of the matrix obtained by projecting R̂_{F_{ω}} onto the subspace T_{h} = S_{h}^{⊥}. We prove that the parameter ω^{∗} depends on the eigenvalues of suitable matrices Q_{i}, i = 1, …, k, defined by means of the matrix G. This dependency in principle does not allow us to explicitly determine ω^{∗}. However, this characterization allows us to devise a strategy for the dynamic choice, at each step n, of a nearly optimal value ω_{n} of ω. This is obtained by replacing at each step n the matrix G with the current approximation X_{n}, thus generating the sequence of matrices

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

We present numerical results which confirm the theoretical properties proved in the paper and which show the effectiveness of the proposed ideas to improve the rate of convergence.

In section 2 we analyze the rate of convergence of functional iterations (2). In section 3 we consider the relaxed sequence (6).

2 ANALYSIS OF THE RATE OF CONVERGENCE

The rate of convergence of functional iterations (2) is strongly related to the spectral properties of the starting approximation X_{0}. Indeed, it was observed in [15], [14] that, in many numerical experiments, the number of iterations needed to reach the required accuracy by starting with X_{0} = I was much less than the number of iterations needed by starting with X_{0} = 0. This behavior has been fully explained in [19], where it is proved that the better rate of convergence is due to the fact that X_{0} shares with G the eigenvalue 1 and the corresponding left eigenvector e^{t} = (1, …, 1). We generalize this result, by considering starting approximations X_{0} which share with G more eigenvalues, and corresponding left eigenvectors.

Denote by Λ = {λ_{1}, …, λ_{k}} the set of the eigenvalues of the matrix G, and with 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 S_{h} the subspace generated by the k^{2}–dimensional vectors e_{i} ⊗ w_{j}, i = 1, …, k, j = 1, …, h, and let T_{h} = S_{h}^{⊥}. Let S_{h} = I ⊗ E_{h} be a k^{2} × hk matrix whose columns are an orthonormal basis of S_{h}, and let T_{h} = I ⊗F_{h} be a k^{2} ×(k^{2} −hk) matrix whose columns are an orthonormal basis of T_{h}.

Let us denote by ê_{n} the vector obtained by column-wise arranging the entries of the matrix G − X_{n}. The error ê_{n+1} can be explicitly related to the error ê_{n}, as shown in [19], by means of the recursion

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

where

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

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}  (8)

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

The mean asymptotic rate of convergence of the sequence ê_{n}, i.e., the value

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

where || · || is any vector norm, can be estimated, as shown in the following theorem:

Theorem 1. Let X_{n}, n ≥ 0, the sequence generated by X_{n+1} = F (X_{n}), where F (·) is defined by (3). Let 1 ≤ h < k and Λ_{h} ⊂ Λ satisfying properties 1, 2, 3. Suppose that X_{0} satisfies

w^{t}X_{0} = λw^{t},

for any λ ∈ Λ_{h}, where w^{t} is the corresponding left eigenvector of G. If the sequence X_{n} is convergent then it holds

r = lim_{n} ^{n} \sqrt{∥ê_{n}∥} ≤ ρ(T_{h}^{t}R̂_{F} T_{h}).

Proof. First observe that ê_{n} ∈ T_{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 (7) and from the fact that, if u ∈ S_{h}, then u^{t}R̂_{n} ∈ S_{h}, for any n ≥ 0.

Let Π = (S_{h}|T_{h}). From (7) it follows that

Π^{t}ê_{n+1} = ( \begin{matrix}& 0 \\ f & n+1 \end{matrix} ) = (Π^{t}R̂_{n}Π)(Π^{t}ê_{n}) = (Π^{t}R̂_{n}Π) ( \begin{matrix}0 \\ f & n \end{matrix} ) ,  (9)

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

Π^{t}R̂_{n}Π = ( \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} ) ,  (10)

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

f_{n+1} = Z_{n}f_{n}.

Let ϵ > 0 be fixed and let || · ||_{ϵ} be a matrix norm such that

∥T_{h}^{t}R̂_{F} T_{h}∥_{ϵ} ≤ ρ(T_{h}^{t}R̂_{F} T_{h}) + ϵ.

Since the sequence {Z_{n}} converges to T_{h}^{t}R̂_{F} T_{h}, it can be easily verified that there exists a positive constant c such that

∥f_{n}∥_{ϵ} ≤ c(ρ(T_{h}^{t}R̂_{F} T_{h}) + 2ϵ)^{n}.

Whence, we have

r = lim_{n} ^{n} \sqrt{∥ê_{n}∥} = lim_{n} ^{n} \sqrt{∥f_{n}∥} ≤ ρ(T_{h}^{t}R̂_{F} T_{h}).

The above result generalizes the convergence properties stated in [19], which can be obtained in the case where h = 1 and X_{0} is nonnegative, that is, X_{0} is a column stochastic matrix.

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 improved. In fact, the following result holds:

Theorem 2. Under the hypotheses of Theorem 1, the set Δ of the eigenvalues of the matrix T_{h}^{t}R̂_{F} T_{h} is given by

Δ = ⋃_{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.

In particular, from the above results, it follows that

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

for any 1 < h < k, that is, the rate of convergence can be improved, if the starting approximation X_{0} has h > 1 eigenvalues, and corresponding left eigenvectors, in common with G.

In the above theorems we have given conditions on the starting matrix X_{0} which allow to increase the rate of convergence of the sequence X_{n+1} = F (X_{n}), under the hypothesis that the sequence X_{n} is convergent. In the case where the matrix X_{0} is column stochastic, the sequence is obviously convergent. In the more general case where X_{0} is not nonnegative, the local convergence is guaranteed by the classical results on functional iteration methods [23], as stated by the next theorem. Let us first introduce some notations: for any k × k matrix A denote by a = vec(A) the k^{2}-dimensional vector obtained by column-wise arranging the entries of the matrix A; analogously, we denote A = vec^{−1}(a).

Theorem 3. Let g_{F} : R^{k^{2}} → R^{k^{2}} the function defined by

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

where F is given in (3), and denote by J_{F} (x) the k^{2} × k^{2} Jacobian matrix associated with the function g_{F}. Then, g_{F} (g) = g and R̂_{F} = J_{F} (g), where g = vec(G).

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

Since ρ(R̂_{F}) < 1 (compare [19]), the local convergence is guaranteed.

Once the eigenvalues λ_{1} = 1, λ_{2}, …, λ_{h}, and the corresponding eigenvectors w_{1}^{t}, …, w_{h}^{t} are known, a matrix X_{0} satisfying (4) can simply be obtained by

X_{0} = V Diag(λ_{1}, …, λ_{h})W^{t},  (11)

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

The eigenvalues λ_{1}, …, λ_{k} of G are given by the zeros of the function φ(z) = det(zI −∑_{j=0}^{+∞} A_{j}z^{j}) lying in the closed unit circle of the complex plane (compare [12]), 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. 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 method [9], [1]; the left eigenvectors can be computed by solving the corresponding homogeneous linear systems.

In the light of theorems 1, 2, 3, choosing an initial approximation satisfying (11) can lead to a strong reduction of the number of iterations needed to reach the required accuracy.

We have verified the effectiveness of the proposed strategy on a problem arising from the modeling of a Metaring MAC protocol [3], having block size k = 16, for different values of the system utilization ρ. For these problems, as ρ approaches 1, functional iteration methods converge very slowly, even if the starting approximation is a column stochastic matrix. We have tested our algorithm in the cases where ρ = 0.5, 0.8, 0.97. 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 this section are satisfied. Table 1 reports, for different values of ρ and h, the number of iterations and the CPU time (in seconds) needed by functional iterations in order to obtain an approximation X with a residual error ∥X − ∑_{i=0}^{+∞} X^{i}A_{i}∥_{1} less than 10^{−11}. We observe that our method shows its effectiveness especially when ρ is close to 1. In particular, in the cases ρ = 0.8 and ρ = 0.97, choosing a starting approximation that shares with G only 2 eigenvalues (and corresponding left eigenvectors), leads to a substantial reduction of the number of iterations.

Figure 1 reports the logarithm (to the base 10) of the residual error for h = 1, 2, 16, in the case where ρ = 0.97. It is worth pointing out the dramatic reduction of the number of iterations when h passes from 1 to 2. In the case where h = 16, i.e., all the eigenvalues and corresponding eigenvectors are used, the matrix X_{0} is a poor approximation of G, which can be refined by applying few steps of functional

ρ = 0.5  ρ = 0.8  ρ = 0.97
h  1 2 8 16 1  2 8 16  1  2 8 16
Iter. 76 32 24 10 204 33 20 8 2763 34 20 19
Time 5.0 2.1 1.6 0.7 15.0 2.4 1.5 0.6 203.3 2.5 1.5 1.4

Table 1: Metaring MAC Protocol

Figure 1: Metaring MAC Protocol, ρ = 0.97

iterations. However, from figure 1, it appears that the asymptotic rate of convergence is almost the same as in the case h = 2.

3 A RELAXATION TECHNIQUE

In the previous section we have observed that the rate of convergence of functional iteration methods is related to the spectral properties of the starting approximation X_{0} and to the Jacobian matrix R̂_{F} associated with the iteration function F. In this section we propose to improve the rate of convergence by applying a relaxation technique to customary functional iterations. More precisely, we consider functional iterations defined by the recursion

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

where F (·) is given by (3) and ω is a relaxation parameter.

It can be readily observed that, in this case, the error ê_{n+1} is related to the error ê_{n} by means of the recursion

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

From the above equation, it follows that the rate of convergence of the relaxed sequence is related to the parameter ω, as stated by the following theorem:

Theorem 4. Let ω > 0 a relaxation parameter. Let 1 ≤ h < k and Λ_{h} ⊂ Λ satisfying properties 1, 2, 3. Suppose that X_{0} satisfies

w^{t}X_{0} = λw^{t},

for any λ ∈ Λ_{h}, where w^{t} is the corresponding left eigenvector of G. If the sequence X_{n}, n ≥ 0 of (12) is convergent, then it holds

r = lim_{n} ^{n} \sqrt{∥ê_{n}∥} ≤ ρ((1 − ω)I + ωT_{h}^{t}R̂_{F} T_{h}).

Proof. The proof can be carried out by following the same lines of the proof of Theorem 2.

The above result suggests us to determine a value ω^{∗} of ω for which the upper bound ρ((1 − ω)I + ωT_{h}^{t}R̂_{F} T_{h}) to the rate of convergence is minimum.

Observe that, if ω = 1 the local convergence of the sequence X_{n} is guaranteed, as pointed out in the previous section, since ρ(T_{h}^{t}R̂_{F} T_{h}) ≤ ρ(R̂_{F}) < 1. By choosing ω = ω^{∗} the local convergence is still guaranteed since ρ((1 − ω^{∗})I + ω^{∗}T_{h}^{t}R̂_{F} T_{h}) ≤ ρ(T_{h}^{t}R̂_{F} T_{h}).

Hereafter we use the notations

W_{ω} = (1 − ω)I + ωT_{h}^{t}R̂_{F} T_{h}

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 5. Let X_{n} = (1 − ω)X_{n−1} + ωF (X_{n−1}), for n ≥ 1. Let 1 ≤ h < k and Λ_{h} ⊂ Λ satisfying properties 1, 2, 3 and suppose that X_{0} satisfies the hypothesis of theorem 4. Let Δ be the set of the eigenvalues of the matrix T_{h}^{t}R̂_{F} T_{h} and Γ = {α ∈ Δ : Re α = max_{λ∈Δ} Re λ}. Then the parameter ω^{∗} such that ρ(W_{ω^{∗}}) = min_{ω>0} ρ(W_{ω}) is given by

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

where:

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

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

Proof. Let us estimate the spectral radius of ρ(W_{ω}). It can be easily verified that the eigenvalues of the matrix W_{ω} are given by the set

{1 − ω + ωλ, λ ∈ Δ}.

Whence, we have

ρ(W_{ω})^{2} = max_{λ∈Δ} p_{λ}(ω),

where

p_{λ}(ω) = |1 − ω + ωλ|^{2} = ω^{2}(1 + |λ|^{2} − 2Re λ) − 2ω(1 − Re λ) + 1.

The function p_{λ}(ω) belongs to a sheaf of k^{2} − hk parabolas with center in (0, 1) and vertices in

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

Since ρ(T_{h}^{t}R̂_{F} T_{h}) < 1, 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_{μ}(λ).

From theorem 2 it follows that the entries of the set Δ, and consequently, the asymptotic optimal relaxation parameter ω^{∗}, 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. Observe that, if X_{0} satisfies the hypothesis of Theorem 4, then the error vector ê_{n} belongs to the linear space T_{h} for any n. Whence, by premultiplying both sides by the matrix Π^{t} as in the proof of Theorem 1, the recurrences (13) can be rewritten in the form

T_{h}^{t}ê_{n+1} = ((1 − ω)I + ωT_{h}^{t}R̂_{n}T_{h}))(T_{h}^{t}ê_{n}).

Based on these considerations, we generate the sequence of matrices

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

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

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

and W_{n,ω} = (1 − ω)I + ωT_{h}^{t}R̂_{n}T_{h}. Indeed, the following results hold:

Theorem 6. Let X_{n} = (1 − ω)X_{n−1} + ωF (X_{n−1}). Let 1 ≤ h < k and Λ_{h} ⊂ Λ satisfying properties 1, 2, 3 and suppose that X_{0} satisfies the hypothesis of theorem 4. Let Δ_{n} be the set of the eigenvalues of the matrix T_{h}^{t}R̂_{n}T_{h} and Γ_{n} = {α ∈ Δ_{n} : Re α = max_{λ∈Δ_{n}} Re λ}. If ρ(T_{h}^{t}R̂_{n}T_{h}) < 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} − Γ_{n} is such that f_{μ}(ν) = min_{λ∈Δ_{n}−Γ_{n}} f_{μ}(λ).

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

In the case where ρ(T_{h}^{t}R̂_{n}T_{h}) ≥ 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 T_{h}^{t}R̂_{n}T_{h}, that can be computed by only knowing the approximation X_{n} and the eigenvalues of the matrix G, as stated by the following result, which is analogous to Theorem 2:

Figure 2: Teletraffic system

Theorem 7. Let

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

where the matrices Y_{j,n} are defined in (8) and λ_{i}, i = 1, …, k, are the eigenvalues of the matrix G. Under the hypothesis of Theorem 1 the set Δ_{n} of the eigenvalues of the matrix T_{h}^{t}R̂_{n}T_{h} is given by

Δ_{n} = ⋃_{i=h+1,...,k}{μ : μ is eigenvalue of Q_{i,n}}

From the many numerical experiments that we performed, we observed that the sequence ω_{n} rapidly converges to the asymptotic optimal relaxation parameter ω^{∗} of Theorem 5. This suggested us to stop the computation of the dynamical relaxation parameter ω_{n} when the difference between the parameters obtained at two subsequent steps is smaller than a positive constant σ, depending on the used machine precision.

We have tested the relaxation technique on the family of teletraffic systems, modeled as a continuous-time QBD Markov chain and defined by blocks of dimension 24, described in [8].

Figure 2 reports, for one of these problems, the logarithm (to the base 10) of the residual error ∥X_{n} − ∑_{i=0}^{+∞} X_{n}^{i}A_{i}∥_{1} of the approximations obtained at step n, with the relaxed and standard functional iteration methods, starting with X_{0} = I.

For this problem, the asymptotic optimal relaxation parameter is ω^{∗} = 2.086. The time needed for the computation of the relaxation parameter does not significantly affect the overall computational time of the relaxed method. The ratio between the number of iterations (and the computational time) of the standard and relaxed method, needed to reach a residual error less than 10^{−12}, is about 2.

4 CONCLUDING REMARKS

In this paper we have proposed two strategies which allow to improve the rate of convergence of functional iterations. They are based on a different use of spectral properties of the matrix G and of related matrices (e.g., R̂_{F}, W_{ω}). For this reason, in general we cannot establish when a strategy is more effective than the other one. In the applications they can also be applied simultaneously, as proposed in Section 3.

REFERENCES

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

[2] N. Akar and K. Sohraby. An invariant subspace approach in M/G/1 and G/M/1 type Markov chains. Commun. Statist. Stochastic Models, 13:381–416, 1997.

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

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

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

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

[7] D. Bini and B. Meini. Inverting block Toeplitz matrices in block Hessenberg form by means of displacement operators: application to queueing problems. Linear Algebra Appl., 272:1–16, 1998.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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