New convergence results on functional iteration techniques for the numerical solution of M/G/1

type Markov chains

Beatrice Meini^{∗} Università di Pisa

Abstract

By performing an accurate analysis of the convergence, we give a complete theoretical explanation of the experimental behaviour of functional iteration techniques 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 (here the A_{i}’s are nonnegative k×k matrices such that the matrix ∑_{i=0}^{+∞} A_{i} is column stochastic). Moreover, we introduce a general class of functional iteration methods, which includes the standard methods, and we give an optimality convergence result in this class.

AMS subject classification. 15A51, 15A24, 60J10, 60K25, 65U05.

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)

where the k × k blocks B_{i+1}, A_{i}, i ≥ 0, are nonnegative and ∑_{i=1}^{+∞} B_{i} and ∑_{i=0}^{+∞} A_{i} are column stochastic matrices. Matrices of the structure (1) are known in literature as stochastic matrices of M/G/1 type [9] and arise in a wide variety of queueing problems modelled by Markov chains M, where P^{T} is the transition matrix associated with M. The computation of the probability invariant vector associated with a stochastic matrix of M/G/1 type (1), i.e., the nonnegative infinite vector π such that π = P π, ∥π∥_{1} = 1, can be reduced (compare [9]) to the computation of the minimal nonnegative solution G of the nonlinear matrix equation

∗Dip. Matematica, v. 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. In the case where the matrix P is irreducible and positive recurrent [4], equation (2) has an unique nonnegative solution, which is column stochastic [9].

The problem of the computation of the nonnegative solution G of (2) has been investigated by several authors. New numerical methods have been recently proposed in [3] and in [7]: in [3] the authors propose a method based on cyclic reduction, which is numerically stable, quadratically convergent and with low computational cost; in [7] the evaluation of the matrix G is reduced to solving, by means of a doubling technique, a block Hessenberg system. However, the most commonly used algorithms for the computation of the matrix G are functional iteration methods based on the recursion

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

where X_{0} is a nonnegative matrix and F (·) is given by

F (X) = ∑_{i=0}^{+∞} X^{i}A_{i}  (4)

or

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

or

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

The sequence {X_{n}} of (3), with F (·) defined by (4), (5) or (6), converges linearly to the matrix G. However a detailed analysis of the speed of convergence, in relation with the particular function F (·) or with the starting matrix X_{0}, has not yet been provided. In [5] Latouche shows that, in the case where X_{0} = O, the sequence (3), where F (·) is (4), or (5), or (6), converges monotonically to the matrix G; moreover, he proves that the method based on (6) is faster than that based on (5), and the method based on (5) is faster than that based on (4), but the rate of convergence of the three sequences is not estimated. In [5] the author also shows that, by choosing as starting matrix X_{0} a general column stochastic matrix, the sequence (3), where F (·) is (4), or (5), or (6), is a sequence of stochastic matrices that still converges to the matrix G. Moreover, Latouche observes [6] that, in many numerical experimentations, the sequence obtained in this way converges much faster than the sequence obtained by starting with X_{0} = O.

In this paper we introduce a class of functional iterations (3) containing the ones described above; then we associate, with each function F defining a functional iteration of this class, a suitable nonnegative k × k matrix R_{F}. In the case where the starting matrix of (3) is X_{0} = O, we first show that the sequence (3) converges monotonically to the matrix G. Then we prove that, for the monotonic convergence of the matrices {X_{n}} to the solution G, the spectral radius ρ_{F} of the matrix R_{F} coincides with the mean asymptotic rate of convergence of the sequence {X_{n}}_{n}. On the other hand, for each function F (·) belonging to the above class, the matrix R_{F} can be obtained as a regular splitting (compare [10]) of the same M-matrix. This fact allows us to prove that the mean asymptotic rates of convergence r_{N}, r_{T}, r_{U} of the methods based on (4), (5), (6), respectively, are related by the following inequality:

r_{U} ≤ r_{T} ≤ r_{N}.

Moreover, we prove that the method based on (6) is the fastest method in the class of functional iterations that we introduced.

In the case where the starting matrix X_{0} of (3) is any column stochastic matrix, we associate with any function F (·) of the above class a nonnegative k^{2}× k^{2} matrix R̂_{F}. The increasing of the speed of convergence of the sequence {X_{n}}, with respect to the case X_{0} = O observed in [6], can be formally explained by means of the matrices R̂_{F}. In fact, we first observe that the spectral radius of the matrix R̂_{F} coincides with the spectral radius of the matrix R_{F}. Then we show that, at each step n, the vector ê_{n} obtained by the column-wise arrangement of the entries of the matrix G − X_{n}, which represents the error at step n, is orthogonal to the nonnegative left eigenvector w^{T} associated with the spectral radius of the matrix R̂_{F}. This fact allows us to prove that the asymptotic rate of convergence r of the sequence {X_{n}} coincides with the spectral radius of the matrix W_{F} obtained by projecting the matrix R̂_{F} onto a suitable linear space E. The linear space E is such that w is orthogonal to any vector belonging to E. Whence, in the case where the matrix R̂_{F} is irreducible the spectral radius of the matrix W_{F} is less than or equal to the second largest modulus eigenvalue of the matrix R̂_{F}; in the case where the matrix R̂_{F} is reducible, the spectral radius of the matrix W_{F} is less than or equal to the spectral radius of R̂_{F}.

The paper is organized as follows. In section 2 we analyze functional iteration methods (3) defined by (4), (5) and (6) with X_{0} = O. In section 3 we introduce the general class of functional iteration methods and we estimate the rate of convergence of the sequences generated in this class, in the case where X_{0} = O. In section 4 we analyze functional iteration methods having, as starting point X_{0}, a column stochastic matrix.

In what follows, e^{T} denotes the k-dimensional vector e^{T} = (1, 1, …, 1), ρ(M) denotes the spectral radius of the matrix M and Ā_{i}, for i ≥ 1, denotes the k × k matrix Ā_{i} = ∑_{j=i}^{+∞} G^{j−i}A_{i}. 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}. We suppose that the matrix P is irreducible and positive recurrent, so that the matrix equation (2) has a unique nonnegative solution. Moreover we suppose that the sequence {X_{n}} of (3) does not reach the solution G in a finite number of steps.

2 Convergence properties of standard functional iteration methods.

In this section we analyze functional iteration methods defined by (3), with X_{0} = O and F (X) given by (4), (5), (6). A convergence analysis of such iterative methods have already been performed by G. Latouche in [5]. More precisely, by following the notations introduced in [5], let us denote by {X_{n}^{(N)}} the sequence defined by (4) (Natural algorithm), by {X_{n}^{(T)}} the sequence defined by (5) (Traditional algorithm) and by {X_{n}^{(U)}} the sequence defined by (6) (algorithm based on the matrix U = ∑_{i=1}^{+∞} G^{i−1}A_{i}), with X_{0}^{(N)} = X_{0}^{(T)} = X_{0}^{(U)} = O. Then the following result holds [5]:

Theorem 1 The sequences {X_{n}^{(N)}}_{n}, {X_{n}^{(T)}}_{n}, {X_{n}^{(U)}}_{n} converge monotonically to the matrix G. Moreover, for any n ≥ 0, it holds

X_{n}^{(N)} ≤ X_{n}^{(T)} ≤ X_{n}^{(U)}.

Theorem 1 allows us to deduce [5] that the number of iterations I_{N}, I_{T}, I_{U} sufficient to obtain the same approximation of the matrix G by means of the sequences {X_{n}^{(N)}}_{n}, {X_{n}^{(T)}}_{n}, {X_{n}^{(U)}}_{n}, respectively, are such that

I_{U} ≤ I_{T} ≤ I_{N}.

Therefore the method based on (6) is faster than that based on (5), and the method based on (5) is faster than that based on (4). We give a more precise result, by estimating also the rate of convergence of the three sequences. Let us first prove the following

Lemma 2 Let X, Y be k × k matrices. Then, for any integer i ≥ 1, it holds that

X^{i} − Y^{i} = ∑_{j=1}^{i} X^{j−1}(X − Y)Y^{i−j}.

Proof. The equality can be simply proved by induction on i. □

Define, for each integer n, the matrices E_{n}^{(N)} = G − X_{n}^{(N)}, E_{n}^{(T)} = G − X_{n}^{(T)} and E_{n}^{(U)} = G − X_{n}^{(U)}, which represent the error at step n for the sequences (3) defined by (4), (5), (6), respectively, with X_{0} = O. Observe that, for the monotonicity stated by theorem 1, the errors at any step n are such that E_{n}^{(N)} ≥ O, E_{n}^{(T)} ≥ O, E_{n}^{(U)} ≥ O. Moreover it is possible to express them recursively, as shown by the following

Proposition 3 For every integer n ≥ 0 it holds that

e^{T} E_{n+1}^{(N)} = e^{T} E_{n}^{(N)}R_{n}^{(N)}, e^{T} E_{n+1}^{(T)} = e^{T} E_{n}^{(T)}R_{n}^{(T)}, e^{T} E_{n+1}^{(U)} = e^{T} E_{n}^{(U)}R_{n}^{(U)}

where, for n ≥ 0,

\begin{matrix}R_{n}^{(N)} = ∑_{i=1}^{+∞} ∑_{j=1}^{i} X_{n}^{(N)}^{i−j}A_{i} \\ R_{n}^{(T)} = (∑_{i=1}^{+∞} ∑_{j=1}^{i} X_{n}^{(T)}^{i−j}A_{i} − A_{1}) (I − A_{1})^{−1} \\ R_{n}^{(U)} = (∑_{i=2}^{+∞} ∑_{j=2}^{i} X_{n}^{(U)}^{i−j}A_{i}) (I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i})^{−1} . \end{matrix}  (7)

Proof. Let us first analyze the functional iteration method based on (4). From lemma 2 and from (2) and (3) it follows that

E_{n+1}^{(N)} = G − X_{n+1}^{(N)} = ∑_{i=1}^{+∞} (G^{i} − X_{n}^{(N)}^{i}) A_{i} = ∑_{i=1}^{+∞} ∑_{j=1}^{i} G^{j−1}E_{n}^{(N)}X_{n}^{(N)}^{i−j}A_{i}.

Whence, since G is column stochastic, we have

e^{T} E_{n+1}^{(N)} = e^{T} ∑_{i=1}^{+∞} ∑_{j=1}^{i} E_{n}^{(N)}X_{n}^{(N)}^{i−j}A_{i} = e^{T} E_{n}^{(N)}R_{n}^{(N)}.

Analogously, for the functional iteration method based on (5), we obtain

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

Hence, from the stochasticity of G, we obtain

e^{T} E_{n+1}^{(T)} = e^{T} E_{n}^{(T)} (∑_{i=2}^{+∞} ∑_{j=1}^{i} X_{n}^{(T)}^{i−j}A_{i}) (I − A_{1})^{−1} = e^{T} E_{n}^{(T)}R_{n}^{(T)}.

Let us now analyze the functional iteration method based on (6). From (3) and (2) we have

E_{n+1}^{(U)} =

G − X_{n+1}^{(U)} = A_{0} [(I − ∑_{i=1}^{+∞} G^{i−1}A_{i})^{−1} − (I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i})^{−1}] =

G [I − (I − ∑_{i=1}^{+∞} G^{i−1}A_{i})(I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i})^{−1}] .

Whence, from lemma 2, it follows that

E_{n+1}^{(U)} =

G [I − (I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i} − ∑_{i=2}^{+∞} ∑_{j=1}^{i−1} G^{j−1}E_{n}^{(U)}X_{n}^{(U)}^{i−1−j}A_{i})·

·(I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i})^{−1}] =

G(∑_{i=2}^{+∞} ∑_{j=1}^{i−1} G^{j−1}E_{n}^{(U)}X_{n}^{(U)}^{i−1−j}A_{i})(I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i})^{−1}.

Therefore, from the stochasticity of G, we have

\begin{matrix}e^{T} E_{n+1}^{(U)} = \\ e^{T} E_{n}^{(U)}(∑_{i=2}^{+∞} ∑_{j=1}^{i−1} X_{n}^{(U)}^{i−1−j}A_{i})(I − ∑_{i=1}^{+∞} X_{n}^{(U)}^{i−1}A_{i})^{−1} = e^{T} E_{n}^{(U)}R_{n}^{(U)} \end{matrix}

□

Proposition 3 allows us to express the norm of the error at step n, by means of the matrices R_{n}^{(N)}, R_{n}^{(T)}, R_{n}^{(U)} of (7):

Corollary 4 For every integer n ≥ 1 it holds that

∥E_{n}^{(N)}∥_{1} = ∥∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1}, ∥E_{n}^{(T)}∥_{1} = ∥∏_{i=0}^{n−1} R_{i}^{(T)}∥_{1}, ∥E_{n}^{(U)}∥_{1} = ∥∏_{i=0}^{n−1} R_{i}^{(U)}∥_{1},

where the matrices R_{i}^{(N)}, R_{i}^{(T)}, R_{i}^{(U)} are defined in (7) and ∏_{i=0}^{n−1} H_{i} = H_{0}H_{1} … H_{n−1}, for k × k matrices H_{i}.

Proof. For every integer n, since E_{n}^{(N)} ≥ O, we have ∥E_{n}^{(N)}∥_{1} = ∥e^{T} E_{n}^{(N)}∥_{1}. Whence, for proposition 3, ∥E_{n}^{(N)}∥_{1} = ∥e^{T} E_{0}^{(N)} ∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1}. On the other hand, since X_{0}^{(N)} = O, we have e^{T} E_{0}^{(N)} = e^{T} G = e^{T}. Therefore, ∥E_{n}^{(N)}∥_{1} = ∥e^{T} ∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1} = ∥ ∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1}. The analogous relations for E_{n}^{(T)} and E_{n}^{(U)} readily follow by using the same arguments. □

It is also possible to estimate the asymptotic rate of convergence of the sequences {X_{n}^{(N)}}, {X_{n}^{(T)}}, {X_{n}^{(U)}}. Let us first define the matrices

\begin{matrix}R^{(N)} = ∑_{i=1}^{+∞} Ā_{i} \\ R^{(T)} = (∑_{i=1}^{+∞} Ā_{i} − A_{1}) (I − A_{1})^{−1} \\ R^{(U)} = (∑_{i=2}^{+∞} Ā_{i}) (I − Ā_{1})^{−1} \end{matrix}  (8)

such that

R^{(N)} = lim_{n} R_{n}^{(N)}, R^{(T)} = lim_{n} R_{n}^{(T)}, R^{(U)} = lim_{n} R_{n}^{(U)}.  (9)

Then the following theorem holds:

Theorem 5 Let

r_{N} = lim_{n} ^{n} \sqrt{∥E_{n}^{(N)}∥}, r_{T} = lim_{n} ^{n} \sqrt{∥E_{n}^{(T)}∥}, r_{U} = lim_{n} ^{n} \sqrt{∥E_{n}^{(U)}∥},

be the mean asymptotic rates of convergence for the sequences {X_{n}^{(N)}}, {X_{n}^{(T)}}, {X_{n}^{(U)}}, respectively, where ∥ · ∥ is any matrix norm. Then it holds that

r_{N} = ρ(R^{(N)}), r_{T} = ρ(R^{(T)}), r_{U} = ρ(R^{(U)}),

where the matrices R^{(N)}, R^{(T)}, R^{(U)} are defined in (8).

Proof. We prove the theorem only for the sequence {X_{n}^{(N)}} and we leave the reader to complete the proof for the sequences {X_{n}^{(T)}} and {X_{n}^{(U)}}. For corollary 4 it holds that

r_{N} = lim_{n} ^{n} \sqrt{∥E_{n}^{(N)}∥_{1}} = lim_{n} ^{n} \sqrt{∥∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1}}.  (10)

For the monotonicity of the sequence {X_{n}^{(N)}} stated by theorem 1 and for the properties of the spectral radius (compare [10]), we find

lim_{n} ^{n} \sqrt{∥∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1}} ≤ lim_{n} ^{n} \sqrt{∥R^{(N)}^{n}∥_{1}} = ρ(R^{(N)}),

whence

r_{N} ≤ ρ(R^{(N)}).  (11)

Let us now prove the opposite inequality. Suppose first that the nonnegative vector e^{T} R_{0}^{(N)} has no null components. For any integer m < n, for the monotonicity of the sequence {X_{n}^{(N)}} and for the monotonicity of ∥ · ∥_{1}, it holds that

∥∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1} ≥ ∥ ∏_{i=0}^{m−1} R_{i}^{(N)}R_{m}^{(N)}^{n−m}∥_{1} ≥ ∥R_{0}^{(N)}^{m}R_{m}^{(N)}^{n−m}∥_{1}.  (12)

Hence, for any integer m, since the nonnegative vector e^{T} R_{0}^{(N)}^{m} also has no null components, it holds that

∥R_{0}^{(N)}^{m}R_{m}^{(N)}^{n−m}∥_{1} ≥ c_{m}∥R_{m}^{(N)}^{n−m}∥_{1}  (13)

where c_{m} > 0 is a suitable constant. From (9) it follows that, for any fixed ϵ > 0, there exists an integer m_{0} such that ρ(R_{m_{0}}^{(N)}) ≥ ρ(R^{(N)}) − ϵ. Whence, from (10), (12) and (13), we obtain

r_{N} ≥ lim_{n} ^{n} \sqrt{∥R_{0}^{(N)}^{m_{0}}R_{m_{0}}^{(N)}^{n−m_{0}} ∥_{1}} ≥ lim_{n} ^{n} \sqrt{c_{m_{0}} ∥R_{m_{0}}^{(N)}^{n−m_{0}} ∥_{1}} ≥
lim_{n} ^{n} \sqrt{c_{m_{0}}} ^{n} \sqrt{∥R_{m_{0}}^{(N)}^{n}∥_{1}} ^{n} \sqrt{∥R_{m_{0}}^{(N)}^{m_{0}} ∥_{1}^{−1}} = lim_{n} ^{n} \sqrt{∥R_{m_{0}}^{(N)}^{n}∥_{1}} =
ρ(R_{m_{0}}^{(N)}) ≥ ρ(R^{(N)}) − ϵ.

Therefore, for the arbitrarity of ϵ, the above inequality, together with (11), leads to the thesis. Suppose now that the first i components of the nonnegative vector e^{T} R_{0}^{(N)} = e^{T} E_{1}^{(N)} are zero and that the remaining components are nonnull. Whence, for the monotonic convergence of the sequence {X_{n}^{(N)}}, it follows that also the first i components of the nonnegative vector e^{T} E_{n}^{(N)} are zero for any n ≥ 1. Moreover, it can easily be proved (we leave the proof to the reader) that, for any n ≥ 0, the first i components of the vector e^{T} R_{n}^{(N)} are zero. Therefore, the matrices R_{n}^{(N)} have the structure

R_{n}^{(N)} = ( \begin{matrix}O & T_{n}^{(N)} \\ O & U_{n}^{(N)} \end{matrix} )

where U_{n}^{(N)} is a (k − i) × (k − i) matrix and T_{n}^{(N)} is an i × (k − i) matrix. Moreover, the matrix R^{(N)} has the structure

R^{(N)} = ( \begin{matrix}O & T^{(N)} \\ O & U^{(N)} \end{matrix} )

where U^{(N)} = lim_{n} U_{n}^{(N)}, T^{(N)} = lim_{n} T_{n}^{(N)}. If, except for the first i components, no other entry of the vector e^{T} E_{n}^{(N)} vanishes in a finite number of steps, then the matrices U_{n}^{(N)} have no null columns for any n ≥ 1. Let v_{1}^{T} be the k − i dimensional vector defined by the nonnull components of the vector e^{T} E_{1}^{(N)}. From the monotonicity of the sequence {X_{n}^{(N)}} and from the monotonicity of ∥ · ∥_{1}, it follows that, for any integer m < n

\begin{matrix}∥∏_{i=0}^{n−1} R_{i}^{(N)}∥_{1} = ∥v_{1}^{T} ∏_{i=1}^{n−1} U_{i}^{(N)}∥_{1} ≥ c∥∏_{i=1}^{n−1} U_{i}^{(N)}∥_{1} ≥ \\ c∥∏_{i=1}^{m−1} U_{i}^{(N)}U_{m}^{(N)}^{n−m}∥_{1} ≥ c∥U_{1}^{(N)}^{m−1}U_{m}^{(N)}^{n−m}∥_{1} ≥ c_{m}∥U_{m}^{(N)}^{n−m}∥_{1}, \end{matrix}

where c, c_{m} > 0 are suitable constants. Whence, for every ϵ > 0, there exists m_{0} such that

r_{N} ≥ lim_{n} ^{n} \sqrt{c_{m_{0}} ∥U_{m_{0}}^{(N)}^{n−m_{0}} ∥_{1}} = ρ(U_{m_{0}}^{(N)}) ≥ ρ(U^{(N)}) − ϵ = ρ(R^{(N)}) − ϵ,

which leads to the thesis, for the arbitrarity of ϵ. If the components i+1, …, i+i_{0} of the vector e^{T} E_{n}^{(N)} vanish in a finite number of steps and no other component vanishes, it follows that, for n ≥ 1, the matrices U_{n}^{(N)} have the structure

U_{n}^{(N)} = ( \begin{matrix}V̂_{n}^{(N)} & T̂_{n}^{(N)} \\ O & Û_{n}^{(N)} \end{matrix} )

where Û_{n}^{(N)} is a (k − i − i_{0}) × (k − i − i_{0}) matrix having no null columns and V̂_{n}^{(N)} is an upper triangular matrix with null diagonal entries. Therefore we apply the same arguments to the matrix Û_{n}^{(N)}, concluding the proof. □

The following proposition, together with theorem 5, allows us to compare the rate of convergence of the three sequences. In particular, it shows that the sequence {X_{n}^{(U)}} based on (6) is asymptotically faster than the sequences {X_{n}^{(N)}} and {X_{n}^{(T)}} based on (4) and (5), respectively:

Proposition 6 The spectral radii of the matrices R^{(N)}, R^{(T)} and R^{(U)} are related by the following inequality

ρ(R^{(U)}) ≤ ρ(R^{(T)}) ≤ ρ(R^{(N)}).

Proof. Observe that

R^{(N)} = C^{(N)}D^{(N)}^{−1}, R^{(T)} = C^{(T)}D^{(T)}^{−1}, R^{(U)} = C^{(U)}D^{(U)}^{−1},

where

\begin{matrix}C^{(N)} = ∑_{i=1}^{+∞} Ā_{i}, & D^{(N)} = I \\ C^{(T)} = ∑_{i=1}^{+∞} Ā_{i} − A_{1}, & D^{(T)} = I − A_{1} \\ C^{(U)} = ∑_{i=2}^{+∞} Ā_{i}, & D^{(U)} = I − Ā_{1} \end{matrix}

are such that

D^{(N)} − C^{(N)} = D^{(T)} − C^{(T)} = D^{(U)} − C^{(U)} = I − ∑_{i=1}^{+∞} Ā_{i},

that is, the matrices R^{(N)}, R^{(T)} and R^{(U)} are obtained by means of regular splittings (compare [10]) of the M-matrix I −∑_{i=1}^{+∞} Ā_{i}. Whence, for the Perron-Frobenius theorem (compare [2]), since

C^{(U)} ≤ C^{(T)} ≤ C^{(N)},

it follows that

ρ(R^{(U)}) ≤ ρ(R^{(T)}) ≤ ρ(R^{(N)}).
□

3 A general class of functional iteration methods.

Functional iteration methods (3) defined by (4), (5), (6) belong to a more general class where the function F (·) is such that the matrix G satisfies the equation X = F (X). More precisely, consider the formal power series A(z) = ∑_{i=0}^{+∞} z^{i}A_{i} (observe that A(z) converges in the closed unit circle) and suppose that A(z) = zH(z) + K(z), where H(z) = ∑_{i=0}^{+∞} z^{i}H_{i} and K(z) = ∑_{i=0}^{+∞} z^{i}K_{i} are such that H_{i} ≥ O, K_{i} ≥ O, i ≥ 0 (observe that A_{0} = K_{0} and A_{i} = H_{i−1} + K_{i}, for i ≥ 1). Then the following property holds

Proposition 7 The matrix G solves equation (2) if and only if it solves the matrix equation

X = K(X)(I − H(X))^{−1}  (14)

where X is a k × k matrix and K(X) = ∑_{i=0}^{+∞} X^{i}K_{i}, H(X) = ∑_{i=0}^{+∞} X^{i}H_{i}.

Proof. The matrix G solves equation (2) if and only if it solves equation

X(I − H(X)) = K(X).

Since H_{i} ≤ A_{i+1} for i ≥ 0 and since the matrix I − ∑_{i=0}^{+∞} G^{i}A_{i+1} is nonsingular (compare [9]), it follows that I − H(G) is nonsingular. Hence G solves (14). □

Equation (14) allows us to define the sequence of matrices

X_{n} = F (X_{n−1}), n ≥ 1

where

F (X) = K(X)(I − H(X))^{−1}.  (15)

By choosing H(z) = 0, H(z) = A_{1}, H(z) = ∑_{i=0}^{+∞} z^{i}A_{i+1}, we obtain the functional iteration methods defined by (4), (5), (6), respectively. The results stated in section 2 can easily be extended to iterations (3) defined by (15), with X_{0} = O. Let us first prove the following result, which extends theorem 1:

Theorem 8 The sequence X_{n} = F (X_{n−1}), X_{0} = O, where F (·) is given by (15), converges monotonically to the matrix G.

Proof. Let us prove the monotonicity of the sequence {X_{n}} by induction on n. For n = 1 the thesis is true, since X_{1} = F (X_{0}) ≥ O = X_{0}. Suppose that X_{n} ≥ X_{n−1} and let us prove that X_{n+1} ≥ X_{n}: from (15) and from the inductive hypothesis, we obtain

X_{n+1} = F (X_{n}) = K(X_{n})(I − H(X_{n}))^{−1} ≥ K(X_{n−1})(I − H(X_{n−1}))^{−1} = X_{n}.

□

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}}. Then the following general recursive relations for the matrices {E_{n}} hold:

Proposition 9 Let X_{n} = F (X_{n−1}), X_{0} ≥ O, where F (·) is defined by (15). Then for every integer n ≥ 0 it holds that

E_{n+1} = (∑_{i=1}^{+∞} ∑_{j=1}^{i}(G^{j−1}E_{n}X_{n}^{i−j}K_{i} + G^{j}E_{n}X_{n}^{i−j}H_{i})) (I − H(X_{n}))^{−1} (16)

or, in equivalent form,

E_{n+1} = (E_{n} ∑_{i=1}^{+∞} X_{n}^{i−1}K_{i} + ∑_{j=2}^{+∞} G^{j−1}E_{n} ∑_{i=j}^{+∞} X_{n}^{i−j}A_{i}) (I − H(X_{n}))^{−1}. (17)

Proof. From (3) and (15) it follows that

E_{n+1} = G − X_{n+1} = K(G)(I − H(G))^{−1} − K(X_{n})(I − H(X_{n}))^{−1} =

∑_{i=0}^{+∞} G^{i}K_{i}(I − H(G))^{−1} − ∑_{i=0}^{+∞} X_{n}^{i}K_{i}(I − H(X_{n}))^{−1} =

∑_{i=0}^{+∞}(G^{i} − X_{n}^{i})K_{i}(I − H(X_{n}))^{−1} + K(G) ((I − H(G))^{−1} − (I − H(X_{n}))^{−1}) =

∑_{i=0}^{+∞}(G^{i} − X_{n}^{i})K_{i}(I − H(X_{n}))^{−1} + G (I − (I − H(G))(I − H(X_{n}))^{−1}) .

Whence, from lemma 2

\begin{matrix}∑_{i=0}^{+∞}(G^{i} − X_{n}^{i})K_{i}(I − H(X_{n}))^{−1} = \\ ∑_{i=1}^{+∞} ∑_{j=1}^{i} G^{j−1}(G − X_{n})X_{n}^{i−j}K_{i}(I − H(X_{n}))^{−1} \end{matrix}

and

I − (I − H(G))(I − H(X_{n}))^{−1} =

I − (I − ∑_{i=0}^{+∞} X_{n}^{i}H_{i} − ∑_{i=1}^{+∞} ∑_{j=1}^{i} G^{j−1}(G − X_{n})X_{n}^{i−j}H_{i}) (I − H(X_{n}))^{−1} =

(∑_{i=1}^{+∞} ∑_{j=1}^{i} G^{j−1}(G − X_{n})X_{n}^{i−j}H_{i}) (I − H(X_{n}))^{−1}.

Therefore

E_{n+1} = (∑_{i=1}^{+∞} ∑_{j=1}^{i}(G^{j−1}E_{n}X_{n}^{i−j}K_{i} + G^{j}E_{n}X_{n}^{i−j}H_{i})) (I − H(X_{n}))^{−1}.

From the above equation and from the relation A_{i} = H_{i−1} + K_{i}, i ≥ 1, we obtain formula (17). □

From the above proposition the generalization of proposition 3 follows:

Corollary 10 Let X_{n} = F (X_{n−1}), X_{0} = O, where F (·) is defined by (15). Then for every integer n ≥ 0 it holds that

e^{T} E_{n+1} = e^{T} E_{n}R_{n}  (18)

where, for n ≥ 0,

R_{n} = (∑_{j=1}^{+∞} ∑_{i=j}^{+∞} X_{n}^{i−j}A_{i} − H(X_{n})) (I − H(X_{n}))^{−1}.  (19)

Proof. For the stochasticity of the matrix G, from (16) and from the relation A_{i} = H_{i−1} + K_{i} we obtain

e^{T} E_{n+1} =
e^{T} E_{n} (∑_{i=1}^{+∞} ∑_{j=1}^{i} X_{n}^{i−j}(K_{i} + H_{i})) (I − H(X_{n}))^{−1} =

e^{T} E_{n} (∑_{j=1}^{+∞} ∑_{i=j}^{+∞} X_{n}^{i−j}A_{i} − H(X_{n})) (I − H(X_{n}))^{−1}.
□

The following result is a straightforward consequence of corollary 10 and a generalization of corollary 4 (we leave the proof to the reader):

Corollary 11 For every integer n ≥ 1 it holds that

∥E_{n}∥_{1} = ∥∏_{i=0}^{n−1} R_{i}∥_{1}  (20)

where the matrices R_{n} are defined in (19) and ∏_{i=0}^{n−1} R_{i} = R_{0}R_{1} … R_{n−1}.

Define the matrix

R_{F} = (∑_{i=1}^{+∞} Ā_{i} − H(G)) (I − H(G))^{−1}  (21)

such that lim_{n} R_{n} = R_{F}. The matrix R_{F} allows us to express the asymptotic rate of convergence of the sequence {X_{n}}, by means of the following

Theorem 12 Let

r = lim_{n} ^{n} \sqrt{∥E_{n}∥}

be the mean asymptotic rate of convergence of the sequence X_{n} = F (X_{n−1}), X_{0} = O, where F (·) is defined in (15) and ∥ · ∥ is any matrix norm. Then

r = ρ(R_{F})

where R_{F} is given in (21).

Proof. The theorem can easily be proved by using the same arguments of theorem 5. □

From the above theorem it follows that the speed of convergence of the sequence X_{n} = F (X_{n−1}), X_{0} = O, defined by (15) is related to the spectral radius of the matrix R_{F} of (21). Observe that the matrix R_{F} can be written in the form

R_{F} = CD^{−1},

where

C = ∑_{i=1}^{+∞} Ā_{i} − H(G), D = I − H(G)

are such that

C ≥ O, D^{−1} ≥ O, D − C = I − ∑_{i=1}^{+∞} Ā_{i},

that is, R_{F} is obtained by means of a regular splitting (compare [10]) of the M-matrix I − ∑_{i=1}^{+∞} Ā_{i}. Suppose that R_{1} = C_{1}D_{1}^{−1} and R_{2} = C_{2}D_{2}^{−1} are the matrices of (21) associated with the functions F_{1}(·), F_{2}(·) of (15), respectively, and suppose that C_{1} ≤ C_{2}. From the Perron-Frobenius theorem (compare [2]), it follows that

ρ(R_{1}) ≤ ρ(R_{2}).

On the other hand observe that, for any choice of F (·) defined by (15), since H_{i} ≥ 0 for any i, it holds that

H(G) ≤ ∑_{i=1}^{+∞} G^{i−1}A_{i} = Ā_{1}.

Whence

C = ∑_{i=1}^{+∞} Ā_{i} − H(G) ≥ ∑_{i=2}^{+∞} Ā_{i} = C^{(U)},

where R^{(U)} = C^{(U)}D^{(U)}^{−1} is the matrix associated with the sequence {X_{n}^{(U)}} defined by (6). Therefore, the functional iteration method defined by (6) is the fastest one in the class of method defined by (15), with X_{0} = O.

4 Functional iteration methods generating a sequence of stochastic matrices.

Consider functional iteration methods

X_{n} = F (X_{n−1}), X_{0} = Q  (22)

where F (·) is given by (15) and Q is a column stochastic matrix. The sequence {X_{n}} is a sequence of stochastic matrices that still converges to the matrix G (compare [5]). We prove that such sequence of matrices converges faster than the sequence of matrices obtained by starting with X_{0} = O.

Proposition 9 expresses recursively the error E_{n} of the sequence X_{n} = F (X_{n−1}), where F (·) is defined in (15) and X_{0} is any nonnegative matrix. In the case where X_{0} = O, since at each step n it holds X_{n} ≤ G, we easily obtain formula (20) by left multiplying (17) by e^{T}. In the case where X_{0} is not necessarily equal to zero (hence E_{n} is not necessarily a nonnegative matrix), in order to express the error at step n in a simpler form, we introduce the matrix tensor product ⊗ such that A ⊗ B is the block matrix defined by the blocks a_{i,j}B, where A = {a_{i,j}}_{i,j}. In fact, if we represent a matrix A by the vector a obtained by column-wise arranging the entries of A, we find that the relation AXB can be rewritten as (B^{T} ⊗ A)x, where A, X, B are matrices of compatible sizes. Hence, the recursive formula (17) can be rewritten as

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

where ê_{n}, for n ≥ 0, is the k^{2}-dimensional vector associated with the matrix E_{n} and

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

where I denotes the k × k identity 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}  (24)

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}

Let Y_{j} = lim_{n} Y_{j,n}, j ≥ 0, and let Y = ∑_{j=0}^{+∞} Y_{j}. Observe that

Y = ((∑_{j=1}^{+∞} ∑_{i=j}^{+∞} G^{i−j}A_{i} − H(G))(I − H(G))^{−1})^{T} = R_{F}^{T} ,

where R_{F} = lim_{n} R_{n} of (19) and is related to the speed of convergence of the sequence X_{n} = F (X_{n−1}), X_{0} = O. From the above equation it follows, in particular, that ρ(Y) = ρ(R_{F}). We prove also that

Proposition 13 The spectral radius of the matrix R̂_{F} = lim_{n} R̂_{n} is given by

ρ(R̂_{F}) = ρ(R_{F}).

Moreover, a nonnegative k^{2}-dimensional vector w such that

w^{T} R̂_{F} = λw^{T} ,

where λ = ρ(R_{F}), is given by w^{T} = (y^{T} ⊗ e^{T}), where y is a nonnegative k-dimensional vector such that R_{F} y = λy.

Proof. Let S be the Schur canonical form [10] of the matrix G. Then, the matrix T = ∑_{j=0}^{+∞}(Y_{j} ⊗ S^{j}) is similar to the matrix R̂_{F}. Moreover, due to the properties of the matrix S, the eigenvalues of the matrix T are given by the set

⋃_{α∈Λ}{μ|μ is eigenvalue of ∑_{j=0}^{+∞} Y_{j}α^{j}},

where Λ is the set of the eigenvalues of G; in particular, we deduce that the eigenvalues of the matrix R_{F} are eigenvalues of the matrix R̂_{F}. On the other hand, for the properties of nonnegative matrices [8], since for any α ∈ Λ it holds that

| ∑_{j=0}^{+∞} Y_{j}α^{j}| ≤ ∑_{j=0}^{+∞} Y_{j}|α^{j}| ≤ ∑_{j=0}^{+∞} Y_{j} = R_{F}^{T} ,

we find that ρ(R̂_{F}) ≤ ρ(R_{F}), whence ρ(R̂_{F}) = ρ(R_{F}). From the properties of tensor product and from the relation e^{T} G = e^{T}, we obtain

(y^{T} ⊗ e^{T})R̂_{F} = (y^{T} ⊗ e^{T}) ∑_{j=0}^{+∞}(Y_{j} ⊗ G^{j}) =
∑_{j=0}^{+∞}(y^{T} Y_{j}) ⊗ (e^{T} G^{j}) = ∑_{j=0}^{+∞}(y^{T} Y_{j}) ⊗ e^{T} =
(y^{T} ∑_{i=0}^{+∞} Y_{j}) ⊗ e^{T} = λ(y^{T} ⊗ e^{T}).

Therefore, the vector w^{T} = y^{T} ⊗ e^{T} is a left nonnegative eigenvector of the matrix R̂_{F}. □

By following the same argument of the proof of proposition 13, it is straightforward to prove the following

Proposition 14 For every n ≥ 0 the spectral radius of the matrix R̂_{n} is given by

ρ(R̂_{n}) = ρ (∑_{j=0}^{+∞} Y_{j,n}) .

Moreover, a nonnegative k^{2}-dimensional vector w_{n} such that

w_{n}^{T} R̂_{n} = λ_{n}w_{n}^{T} ,

where λ_{n} = ρ(R_{n}), is given by w_{n}^{T} = (y_{n}^{T} ⊗ e^{T}), where y_{n} is a nonnegative k-dimensional vector such that y_{n}^{T} ∑_{j=0}^{+∞} Y_{j,n} = λ_{n}y_{n}^{T}.

From proposition 14 it follows that, for any integer n ≥ 0, the left nonnegative eigenvector w_{n} of the matrix R̂_{n} associated with its spectral radius λ_{n} belongs to the linear space S_{1} 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),  (25)

where e_{i} is the k-dimensional vector having the i-th entry equal to 1 and the remaining entries equal to zero. Let S_{2} be the linear space S_{2} = S_{1}^{⊥}. It can easily be observed that, since the matrices X_{n} and G are column stochastic, the error vector ê_{n} belongs to the linear space S_{2}, for any n. In particular, if S_{i} is the matrix whose columns are an orthogonal basis of the linear space S_{i}, for i = 1, 2, then the k^{2} × k^{2} matrix S = (S_{1}|S_{2}) is such that

S^{T} S = I, S^{T} ê_{n} = ( \begin{matrix}0 \\ f_{n} \end{matrix} ) ,

where f_{n} is a k(k − 1)-dimensional vector. Whence the recursive equations (23) can be rewritten in the form

( \begin{matrix}0 \\ f_{n+1} \end{matrix} ) = (S^{T} R̂_{n}S)(S^{T} ê_{n}) = (S^{T} R̂_{n}S) ( \begin{matrix}0 \\ f_{n} \end{matrix} ) .  (26)

On the other hand, observe that, if w ∈ S_{1} then, from (24) and from the properties of tensor product, it follows that (w^{T} R̂_{n})^{T} ∈ S_{1}, whence S_{1}^{T} R̂_{n}S_{2} = O. Therefore the matrix S^{T} R̂_{n}S has the structure

S^{T} R̂_{n}S = ( \begin{matrix}S_{1}^{T} \\ − \\ S_{2}^{T} \end{matrix} ) R̂_{n} (S_{1}|S_{2}) = ( \begin{matrix}V_{n} & O \\ T_{n} & W_{n} \end{matrix} ) ,  (27)

where V_{n} = S_{1}^{T} R̂_{n}S_{1}, W_{n} = S_{2}^{T} R̂_{n}S_{2} and T_{n} = S_{2}^{T} R̂_{n}S_{1}. Hence the recurrences (26) can be rewritten in the form

f_{n+1} = W_{n}f_{n}.  (28)

The above properties allow us to prove the following convergence results:

Theorem 15 Let

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

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

r ≤ ρ(W_{F})

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

Proof. From (26) and (28) it follows that

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

Therefore, let us analyze the convergence of the sequence {f_{n}}. Let ϵ > 0 be fixed. Then (compare [10]) there exists a matrix norm ∥ · ∥_{ϵ} such that

∥W_{F} ∥_{ϵ} ≤ ρ(W_{F}) + ϵ.

On the other hand, since the sequence {W_{n}} converges to W_{F}, there exists an integer i_{0} such that

∥W_{i}∥_{ϵ} ≤ ∥W_{F} ∥_{ϵ} + ϵ

for any i ≥ i_{0}. Therefore

\begin{matrix}∥f_{n}∥_{ϵ} = ∥W_{n−1}W_{n−2} … W_{0}f_{0}∥_{ϵ} ≤ c_{i_{0}}∥W_{n−1}∥_{ϵ}∥W_{n−2}∥_{ϵ} … ∥W_{i_{0}}∥_{ϵ} ≤ \\ c_{i_{0}}(∥W_{F} ∥_{ϵ} + ϵ)^{n−i_{0}} ≤ c_{i_{0}}(ρ(W_{F}) + 2ϵ)^{n−i_{0}} \end{matrix}

where c_{i_{0}} is a positive constant. Hence, by taking the n-th root, we have, for n → +∞,

r = lim_{n} ^{n} \sqrt{∥f_{n}∥_{ϵ}} ≤ ρ(W_{F}) + 2ϵ.

Therefore, from the arbitrarity of ϵ and from (29), we conclude the proof. □

Observe that, from the structure (27) of R̂_{F}, it holds ρ(W_{F}) ≤ ρ(R̂_{F}) = ρ(R_{F}). In the case where the matrix R̂_{F} is irreducible the inequality is strict, as shown by the following:

Theorem 16 If the matrix R̂_{F} is irreducible, then the mean asymptotic rate of convergence

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

of the sequence (22), where F (·) is defined by (15) and ∥ · ∥ is any vector norm, is

r ≤ ρ_{2}(R̂_{F})

where ρ_{2}(R̂_{F}) denotes the second largest modulus eigenvalue of the matrix R̂_{F}.

Proof. From theorem 15 it is sufficient to prove that

ρ(W_{F}) ≤ ρ_{2}(R̂_{F}).

First observe that, from (27), the eigenvalues of the matrix R̂_{n} coincide with the eigenvalues of the matrices V_{n} and W_{n}. On the other hand, if the columns of the matrix S_{1} are given by the vectors v_{1}, …, v_{k} of (25), then the matrix V_{n} is nonnegative and such that, for the properties of tensor product and from the structure (24) of R̂_{n},

V_{n} = S_{1}^{T} R̂_{n}S_{1} = \frac{1}{k} (I ⊗ e^{T})R̂_{n}(I ⊗ e) = ∑_{j=0}^{+∞} Y_{j,n}.

Therefore, for proposition 14, the spectral radius of the matrix V_{n} is equal to the spectral radius of the matrix R̂_{n}, for any n. On the other hand, since the matrix R̂_{F} is irreducible, it follows that there exists an integer n_{0} such that the matrices R̂_{n} are irreducible for n ≥ n_{0}. Whence, for the Perron Frobenius theorem (compare [2]), for any n ≥ n_{0}, there exists a unique positive eigenvalue of the matrix R̂_{n} having modulus equal to the spectral radius and its algebraic multiplicity is equal to 1. Therefore, for any n ≥ n_{0}, the spectral radius of the matrix W_{n} is such that

ρ(W_{n}) ≤ ρ_{2}(R̂_{n}),

where ρ_{2}(R̂_{n}) is the second largest modulus eigenvalue of the matrix R̂_{n}. In particular, by taking limits for n → +∞, we obtain

ρ(W_{F}) ≤ ρ_{2}(R̂_{F}),  (29)

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

Let us associate with a function F (·) of (15) the matrix

R̂_{F} = ∑_{i=1}^{+∞} (G^{i−1}K_{i}(I − H(G))^{−1})^{T} ⊗ I+
∑_{j=2}^{+∞} ∑_{i=j}^{+∞} (G^{i−j}A_{i}(I − H(G))^{−1})^{T} ⊗ G^{j−1} =
((I − H(G))^{−1^{T}} ⊗ I) (∑_{j=1}^{+∞} Ā_{j}^{T} ⊗ G^{j−1} − H(G)^{T} ⊗ I) .

From the above theorems and from the results of the previous section, it follows that, in the cases where the matrix R̂_{F} is irreducible, then the mean asymptotic rate of convergence of the sequence X_{n} = F (X_{n−1}) is:

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

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

In many cases the modulus of the second largest eigenvalue of a matrix can be much smaller than its spectral radius. In the following figures we report, for the functional iterations defined by (4), (5), (6) the logarithm (to the base 10) of the residual ∥X_{n} − ∑_{i=0}^{+∞} X_{n}^{i}A_{i}∥_{1}, for the sequences obtained by starting with X_{0} = O and with X_{0} = I, for a problem arising from the modelling of a metropolitan network [1]. It is worth pointing out the increasing of the speed of convergence obtained by starting with X_{0} = I instead of X_{0} = O. Moreover it can be observed that the method defined by (6) (method based on the matrix U) converges more quickly than the methods based on (4) and (5) (natural and traditional algorithm, respectively). More precisely the spectral radius ρ_{1} of the matrix R̂_{F}, which gives the mean asymptotic rate of convergence of the sequence obtained by starting with X_{0} = O, is given by ρ_{1} = 0.998448 for the natural algorithm defined by (4), ρ_{1} = 0.998379 for the traditional algorithm defined by (5) and ρ_{1} = 0.998176 for the algorithm based on the matrix U defined by (6). On the other hand the second largest modulus eigenvalue of the matrix R̂_{F} is ρ_{2} = 0.883677 for the natural algorithm, ρ_{2} = 0.875044 for the traditional algorithm and ρ_{2} = 0.858737 for the algorithm based on the matrix U. Moreover, as it can easily be observed from the figures, for this particular example the rate of convergence of the sequence obtained by starting with X_{0} = I is equal to the second largest modulus eigenvalue of the matrix R̂_{F}.

References

[1] Anastasi, G., Lenzini L. (1995): Personal communication.

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

[3] Bini D., Meini B. (1996): On the solution of a nonlinear matrix equation arising in queueing problems. SIAM J. Matrix Anal. Appl., 17(4), pp. 906– 926.

[4] Çinlar E. (1975): Introduction to Stochastic Processes, Prentice-Hall, Englewood Cliffs, N.J..

[5] Latouche G. (1992): Algorithms for evaluating the matrix G in Markov chains of PH/G/1 type. Bellcore Technical Report.

Figure 1: Natural algorithm Figure 2: Traditional algorithm

Figure 3: Algorithm based on matrix U

[6] Latouche G. (1995): Personal communication.

[7] Latouche G., Stewart G.W. (1996): Numerical methods for M/G/1 type queues. in Computations with Markov Chains, W. J. Stewart Ed., pp. 571– 581, Kluwer Academic Publishers.

[8] Minc H. (1988): Nonnegative Matrices, John Wiley and Sons, Inc..

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

[10] Varga R.S. (1962): Matrix Iterative Analysis, New Jersey, Prentice Hall .