Solving certain queueing problems by means of regular splittings

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

Abstract

We analyze the problem of the computation of the solution of the nonlinear matrix equation X = ∑_{i=0}^{+∞} X^{i}A_{i}, arising in queueing models. We propose a technique based on regular splittings, that one hand leads to a new method for computing the solution, on the other hand it may be used to construct nonlinear matrix equations equivalent to starting one, that can be possibly solved by applying different algorithms.

Key Words. Regular splitting, Markov chain, 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 [5] and arise in a wide variety of queueing problems modeled by a Markov chain 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)

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 [5]) to the computation of the minimal nonnegative solution G of the nonlinear matrix equation

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

favati@imc.pi.cnr.it

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

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

where X is a k × k matrix. If the matrix P is irreducible and positive recurrent [2], equation (3) has a unique nonnegative solution, which is column stochastic [5].

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 [6], which involves the block entries A_{i}, B_{i} of the matrix P in (1), and G.

In this paper we derive a new method for solving the matrix equation (3), that consists in rewriting equation (3) in terms of a linear system and in applying an iterative method, based on regular splittings, for its solution. More precisely, the equation (3) can be rewritten as the following block Toeplitz block Hessenberg infinite system:

[G, G^{2}, G^{3}, …] [ \begin{matrix}I − A_{1} & −A_{0} & & ◯ \\ −A_{2} & I − A_{1} & −A_{0} \\ −A_{3} & −A_{2} & I − A_{1} & −A_{0} \\ ⋮ & ... & ... & ... & ... \end{matrix} ] = [A_{0}, O, O, …] .  (4)

In order to solve the above system we generate, by means of regular splittings, a sequence of equivalent systems having block Toeplitz block Hessenberg matrices. We prove that the sequence of transformed systems converges to the system:

[G, G^{2}, G^{3}, …] [ \begin{matrix}I & −A_{0}^{∗} & & ◯ \\ & I & −A_{0}^{∗} \\ ◯ & & ... & ... \end{matrix} ] = [A_{0}^{∗}, O, O, …] ,  (5)

from which we obtain G = A_{0}^{∗}. In this way we generate a sequence of nonlinear matrix equations X = ∑_{i=0}^{∞} X^{i}A_{i}^{(n)}, n ≥ 0, whose solution is still G, that converges to the linear matrix equation X = A_{0}^{∗}, that is immediately solved.

This approach one hand leads to a new method for computing the solution G of (3), on the other hand it may be used to construct nonlinear matrix equations equivalent to (3), that can be possibly solved by applying different algorithms. Indeed, we prove that functional iterations applied to the matrix equation obtained after few steps of the method, converge faster than functional iterations applied to (3), when the initial approximation is the null matrix.

The paper is organized as follows. In Section 2 we analyze the convergence properties of the proposed method, considering also the case where the matrix power series ∑_{i=0}^{∞} A_{i}z^{i} is rational. In Section 3 we relate the regular splittings method with functional iterations, and show some numerical results.

2 The regular splitting method and its convergence properties.

In this section we explain the idea which the method is based on and analyze its convergence properties.

Denoted by H, C and X the coefficient matrix, the right hand side and the unknown of system (4) respectively, we consider the following splitting: H = M − N where

M = [ \begin{matrix}I − A_{1} & & ◯ \\ −A_{2} & I − A_{1} \\ −A_{3} & −A_{2} & I − A_{1} \\ ⋮ & ... & ... & ... \end{matrix} ] , N = [ \begin{matrix}O & A_{0} & & ◯ \\ & O & A_{0} \\ ◯ & & ... & ... \end{matrix} ]

obtaining the equivalent system

X(I − NM^{−1}) = CM^{−1}.  (6)

M is a nonsingular matrix such that M^{−1} is nonnegative and N is a nonnegative matrix. Splittings with these properties are called regular splittings in [8] and are encountered in the numerical solution of finite Markov chains by means iterative methods (see [7]).

Since M and N are block triangular block Toeplitz matrices, the coefficient matrix H^{(1)} and the right hand side C^{(1)} of the resulting system have the same structure of the initial ones; namely

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

and

C^{(1)} = CM^{−1} = [A_{0}^{(1)}, O, O, …] ,

where A_{i}^{(1)} = A_{0}B_{i+1}, i = 0, 1, …, and B_{1} = (I − A_{1})^{−1}, B_{i} = (I − A_{1})^{−1} ∑_{h=1}^{i−1} A_{i+1−h}B_{h} , i ≥ 2, are the block entries of the first block column of M^{−1}. Moreover, it can be easily verified that A_{i}^{(1)} ≥ 0, i = 0, 1, … and that ∑_{i=0}^{+∞} A_{i}^{(1)} is a column stochastic matrix.

So we may iterate this transformation process obtaining a sequence of equivalent systems

[G, G^{2}, G^{3}, …] H^{(j)} = [A_{0}^{(j)}, O, O, …] , j ≥ 1,  (7)

with

H^{(j)} = [ \begin{matrix}I − A_{1}^{(j)} & −A_{0}^{(j)} & & ◯ \\ −A_{2}^{(j)} & I − A_{1}^{(j)} & −A_{0}^{(j)} \\ −A_{3}^{(j)} & −A_{2}^{(j)} & I − A_{1}^{(j)} & −A_{0}^{(j)} \\ ⋮ & ... & ... & ... & ... \end{matrix} ] ,  (8)

where the blocks A_{i}^{(j)} are defined by the recursions:

A_{i}^{(j+1)} = A_{0}^{(j)}B_{i+1}^{(j)}, i = 0, 1, …, j ≥ 1,  (9)

and

B_{1}^{(j)} = (I − A_{1}^{(j)})^{−1}, B_{i}^{(j)} = (I − A_{1}^{(j)})^{−1} ∑_{h=1}^{i−1} A_{i+1−h}^{(j)}B_{h}^{(j)}, i ≥ 2, j ≥ 1,  (10)

From (7) and (8) it follows that G = ∑_{i=0}^{∞} G^{i}A_{i}^{(j)}; in this way at each step of our method we construct a new nonlinear matrix equation having the same nonnegative solution G.

The formulae (9–10) relating the blocks {A_{i}^{(j)}}_{i} obtained at two subsequent steps can be expressed in functional form in terms of formal matrix power series. Indeed, if we associate with the sequence {A_{i}^{(j)}}_{i}, at step j the formal matrix power series φ^{(j)}(z) = ∑_{i=0}^{∞} A_{i}^{(j)}z^{i} we obtain that, for j ≥ 0,:

φ^{(j+1)}(z) = φ^{(j)}(0)(I − ∑_{i=1}^{∞} A_{i}^{(j)}z^{i−1})^{−1} = φ^{(j)}(0) (I − \frac{φ^{(j)}(z) − φ^{(j)}(0)}{z})^{−1} .  (11)

We prove that, if the matrix A_{0} is nonsingular, the generated sequence of systems converges to a system that is easy to solve. More precisely, we prove that the sequences {H^{(j)}}_{j} converges to the matrix

[ \begin{matrix}I & −A_{0}^{∗} & & ◯ \\ & I & −A_{0}^{∗} \\ ◯ & & ... & ... \end{matrix} ]

where A_{0}^{∗} = lim_{j→∞} A_{0}^{(j)}. From (7) we conclude that G = lim_{j→∞} A_{0}^{(j)}. This result is expressed in functional form by the following theorem.

Theorem 1 Let {A_{i}}_{i} be a sequence of k × k nonnegative matrices, such that ∑_{i=0}^{+∞} A_{i} is column stochastic and det A_{0} ≠ 0. Let φ^{(j)}(z) = ∑_{i=0}^{∞} A_{i}^{(j)}z^{i}, j = 1, 2, … be defined according to (11),where φ^{(0)}(z) = ∑_{i=0}^{∞} A_{i}z^{i}, then lim_{j→∞} φ^{(j)}(z) = lim_{j→∞} A_{0}^{(j)} = A_{0}^{∗} and A_{0}^{∗} = G.

Proof. First we prove that the lim_{j→∞} A_{0}^{(j)} exists. Since the following inequalities hold

A_{0}^{(j+1)} = A_{0}^{(j)}(I − A_{1}^{(j)})^{−1} ≥ A_{0}^{(j)}, e^{T} A_{0}^{(j)} ≤ e^{T} , ∀j

the sequence {A_{0}^{j}} is nondecreasing, bounded and hence convergent. Denote with A_{0}^{∗} the limit matrix.

Then we prove by induction on i that lim_{j→∞} A_{i}^{(j)} = 0, i ≥ 1:

• Step 1: lim_{j→∞} A_{1}^{(j)} = 0. Indeed

A_{0}^{(j)} = A_{0}^{(0)}(I − A_{1}^{(0)})^{−1}(I − A_{1}^{(1)})^{−1} … (I − A_{1}^{(j−1)})^{−1}

implies

A_{0}^{∗} = A_{0} ∏_{j=0}^{∞}(I − A_{1}^{(j)})^{−1}.

If det A_{0} ≠ 0 then ∏_{j=0}^{∞}(I − A_{1}^{(j)})^{−1} = A_{0}^{−1}A_{0}^{∗} is a convergent product. On the other hand, it is easy to see that

∀n ≥ 0, ∏_{j=0}^{n}(I − A_{1}^{(j)})^{−1} ≥ I + ∑_{j=0}^{n} A_{1}^{(j)},

hence also the series ∑_{j=0}^{∞} A_{1}^{(j)} is convergent, implying lim_{j→∞} A_{1}^{(j)} = 0.

• Inductive step: assuming that lim_{j→∞} A_{k}^{(j)} = 0, k = 1, …, i we prove that lim_{j→∞} A_{i+1}^{(j)} = 0. From relations (9) and (10) we have

A_{i}^{(j+1)} = A_{0}^{(j)}B_{i+1}^{(j)} = A_{0}^{(j)}(I − A_{1}^{(j)})^{−1} ∑_{h=1}^{i} A_{i+2−h}^{(j)}B_{h}^{(j)}.

We isolate in the sum the term containing the matrix A_{i+1}^{(j)}, obtaining

A_{i}^{(j+1)} = A_{0}^{(j+1)}A_{i+1}^{(j)}B_{1}^{(j)} + A_{0}^{(j+1)} ∑_{h=2}^{i} A_{i+2−h}^{(j)}B_{h}^{(j)}.  (12)

For j going to ∞, the left hand side of (12) goes to 0, by inductive hypotheses; since A_{0}^{(j+1)} is bounded and ∑_{h=2}^{i} A_{i+2−h}^{(j)}B_{h}^{(j)} contains matrices A_{k}^{(j)} with k ≤ i the second term in the right hand side of (12) goes to 0, too. Hence

lim_{j→∞} A_{0}^{(j+1)}A_{i+1}^{(j)}B_{1}^{(j)} = 0.  (13)

If, by contradiction, A_{i+1}^{(j)} does not converge to the null matrix, a subsequence {A_{i+1}^{(j_{h})}}_{h} and a nonnegative matrix R, R ≠ 0 exist such that A_{i+1}^{(j_{h})} ≥ R, for any h. Then A_{0}^{(j_{h}+1)}A_{i+1}^{(j_{h})}B_{1}^{(j_{h})} ≥ A_{0}A_{i+1}^{(j_{h})} ≥ A_{0}R, ∀h; that is the subsequence {A_{0}^{(j_{h}+1)}A_{i+1}^{(j_{h})}B_{1}^{(j_{h})}}_{h} is bounded away from zero, being det A_{0} ≠ 0, and this contradicts (13). □

In the rational case, that is when φ^{(0)}(z) is written as a fraction of two polynomial matrices

φ^{(0)}(z) = P (z)Q(z)^{−1},

the sequence {φ^{(j)}(z)}_{j} is expressed in the same way

φ^{(j)}(z) = P^{(j)}(z)Q^{(j)}(z)^{−1}.

It is easy to prove that the polynomial matrices associated with φ^{(j+1)}(z) and with φ^{(j)}(z) are related by the following formulae:

\begin{matrix}P^{(j+1)}(z) = P^{(j)}(0)Q^{(j)}(0)^{−1}Q^{(j)}(z), \\ Q^{(j+1)}(z) = Q^{(j)}(z) − \frac{P^{(j)}(z)−P^{(j+1)}(z)}{z}, j ≥ 0. \end{matrix}  (14)

Moreover if P (z) and Q(z) are polynomials of degree p_{0} and q_{0} respectively, the degrees of the two sequences {P^{(j)}(z)}_{j} and {Q^{(j)}(z)}_{j} do not increase:

p_{j} = deg(P^{(j)}(z)) ≤ max(p_{0}, q_{0}),
q_{j} = deg(Q^{(j)}(z)) ≤ max(p_{0}, q_{0}).

Hence in the rational case, it is more convenient to use formulae (14), rather than (11), since the degree of φ^{(j)}(z) could increase.

3 Regular splittings and functional iteration methods.

In this section we point out some relations between regular splitting method and functional iteration methods. These considerations, together with the results of numerical experimentations, suggest us to apply a few steps of the method, and then to solve the nonlinear matrix equation obtained in this way by functional iteration methods.

The most commonly used functional iteration method are based on the recursion:

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

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

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

or

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

or

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

In [3] and [4] it is proved that in the case where X_{0} = 0 the method based on (18) is faster than the method based on (17) and that the method based on (17) is faster than the method based on (16). This property results experimentally true also when X_{0} is a stochastic matrix.

We observe that method based on (17) coincides with method based on (16), applied to solve the matrix equation X = ∑_{i=0,i≠1}^{+∞} X^{i}Ā_{i} where Ā_{i} = A_{i}(I − A_{1})^{−1}, i = 0, 2, 3, …, are obtained by writing (4) as

[G, G^{2}, G^{3}, …] (I − N̄M̄^{−1}) = [A_{0}, O, O, …] M̄^{−1},

where

M̄ = [ \begin{matrix}I − A_{1} & & ◯ \\ & I − A_{1} \\ ◯ & & ... \end{matrix} ] , N̄ = [ \begin{matrix}O & A_{0} & & ◯ \\ A_{2} & O & A_{0} \\ A_{3} & A_{2} & O & ... \\ ⋮ & ... & ... & ... \end{matrix} ] .

This splitting cannot be iterated since the blocks on the main diagonal of H̄^{(1)} = I − N̄M̄^{−1} are identity matrices.

Consider now the matrix equation

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

where the blocks A_{i}^{(1)} are obtained after one step of the regular splitting method described in the previous section. From (11), it follows that the method based on (16), applied to solve (19) is the functional iteration method based on (18), that is the fastest one. Thus, if we apply functional iteration (16) to the matrix equation X = ∑_{i=0}^{∞} X^{i}A_{i}^{(2)}, obtained by applying one more step of regular splitting, we have a further improvement of the rate of convergence. This observation suggests us to apply a few steps of regular splittings, and then to apply functional iteration method (16).

We have tested this idea to solve a problem arising in telecommunication modeling, where the size of the blocks A_{i} is 16, and the number of nonzero blocks is 241 (we refer to [1] for more details). We have applied h = 1, 5, 10 regular splittings, and then functional iteration (16), with X_{0} = 0 and X_{0} = I. In the following figures we report the logarithm of the residual error ||X_{n}−∑_{i=0}^{∞} X_{n}^{i}A_{i}^{(h)}||_{1} of X_{n}, versus the number of iterations n, for the different values of h.

From the figures we observe that the asymptotic convergence is improved, as h grows, also in the case X_{0} = I.

References

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

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

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

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

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

Figure 1: Functional iterations after h regular splittings, X_{0} = 0 Figure 2: Functional iterations after h regular splittings, X_{0} = I [6] 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.

[7] D.J. Rose, “Convergent regular splittings for singular M-matrices”, SIAM J. Alg. Discr. Methods, \underline{5} (1984), pp. 133–144.

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