Hostname: page-component-7b9c58cd5d-sk4tg Total loading time: 0 Render date: 2025-03-15T13:01:34.081Z Has data issue: false hasContentIssue false

GENERALIZED CLASS [Cscr ] MARKOV CHAINS AND COMPUTATION OF CLOSED-FORM BOUNDING DISTRIBUTIONS

Published online by Cambridge University Press:  27 February 2007

Mouad Ben Mamoun
Affiliation:
PRiSM, Université Versailles St Quentin, 78035 Versailles, France, and, Université Mohammed V, Rabat, Maroc
Ana Bušić
Affiliation:
PRiSM, Université Versailles St Quentin, 78035 Versailles, France
Nihal Pekergin
Affiliation:
PRiSM, Université Versailles St Quentin, 78035 Versailles, France, and, Centre Marin Mersenne, Université Paris 1, 75013 Paris, France, E-mail: {mobe,abusic,nih}@prism.uvsq.fr
Rights & Permissions [Opens in a new window]

Abstract

In this article we first give a characterization of a class of probability transition matrices having closed-form solutions for transient distributions and the steady-state distribution. We propose to apply the stochastic comparison approach to construct bounding chains belonging to this class. Therefore, bounding chains can be analyzed efficiently through closed-form solutions in order to provide bounds on the distributions of the considered Markov chain. We present algorithms to construct upper-bounding matrices in the sense of the ≤st and ≤icx order.

Type
Research Article
Copyright
© 2007 Cambridge University Press

1. INTRODUCTION

Markovian models have been largely used in computer systems performance evaluation as well as in different areas of applied probability to analyze quantitatively the performance measures of interest. Despite considerable works on the numerical analysis of Markov chains, the analysis of large-size models still remains a problem [17]. The state space reduction techniques, such as lumpable or quasilumpable Markov chains, and models having specific numerical solutions, such as product-form and matrix geometric solutions, have been largely studied to overcome the state space explosion problem [17]. It is also usual to provide approximations on the performance measures of interest by ignoring or simplifying the complicating aspects of the underlying models. Bounding methods have received attention among other approximative methods since they provide bounds on the performance measures of interest. These methods differ in concepts on which they are based and in types of obtained bound.

The bounding method proposed by Courtois and Semal [6] and then extended by other authors [12,14] provides bounds on steady-state rewards without the computation of the steady-state distribution. Since the analysis is performed in steady state, the problem can be considered as a linear algebra problem, and by using polyhedra theory, the lower and the upper bounds on the steady-state rewards are computed. Although this method provides quite accurate results for some classes of model, there are some constraints on the structures of the underlying Markov chains for its application.

The Markov reward approach proposed by van Dijk [18] and van Dijk and van der Wal [19] provides comparison results for expected cumulative and steady-state rewards of two Markov chains. This approach has been successfully applied to provide bounds for several queuing network models that do not have product-form solutions. The general theorems and their conditions as well as technical steps are given in the survey article [18]. The comparison between the rewards of the original model and the bounding model is established on an expectation basis without providing a comparison of the underlying distributions. Since a specific reward is considered in this approach, the comparison constraints can be, in some cases, relaxed compared to the comparison on a sample path basis. Although there are general comparison and error bounding theorems, the establishment of technical steps of proofs is model driven.

There are also bounding methods and comparison results based on the stochastic comparison that have been largely applied in different areas of applied probability: economics, reliability, and queuing networks. We state here two self-contained books as central references [13,16]. Several stochastic orders have been proposed and the best known is the strong stochastic order ≤st. Different comparison results have been established in this stochastic order sense by coupling constructions, by sample-path comparisons, or by means of the stochastic monotonicity. Among the large literature on ≤st comparison results, let us state here some articles that construct bounding Markov chains to overcome the state space explosion problem. In [15], the state space reduction by taking into account the dynamic of the underlying model to establish bounding models has been discussed. In [8,11], the underlying models have been modified to have bounding models that can be analyzed by means of matrix-geometric methods or by the stochastic complementation method. The comparison results are established by sample-path comparison techniques.

Sufficient conditions for the existence of stochastic comparison of two time-homogeneous Markov chains are given by the stochastic monotonicity and bounding properties of their one-step transition probability matrices [13]. In [1], this idea is used to construct a monotone bounding Markov chain in the sense of the ≤st order on a reduced-size state space. Furthermore, this method is extended to construct monotone, numerically efficient bounding chains [7]. In the sequel, this approach is called the algorithmic stochastic comparison. For a given Markov chain, we construct a monotone bounding matrix that is endowed with some properties that lead to simpler numerical computations or analytical expressions. This approach is not model driven, so it can be applied for any Markov chain; however, the bounds might be not accurate for some cases.

In [3,4], a particular class of Markov chains having a closed-form solution for the steady-state distribution has been introduced. These Markov chains are called class [Cscr ] Markov chains and the closed-form solution of their transient distributions has been given in [5]. In this article, we first give the extension of class [Cscr ] Markov chains to [Cscr ][Gscr ] Markov chains also having closed-form solutions for transient and steady-state distributions. Second, we present algorithms to construct [Cscr ][Gscr ] bounding matrices in the sense of the ≤st and ≤icx stochastic orders. Therefore, we have the stochastic comparison between the transient and the steady-state distribution of the underlying model and the closed-form bounding distributions.

The variability orders (≤icx and ≤icv) are especially used to compare models to emphasize the impact of model parameters or policies, but not to construct bounding Markov chains. Contrary to the ≤st order, these stochastic orders represent some difficulties in constructing monotone bounding matrices [1]. In [4], an algorithm to construct an ≤icx-monotone, class [Cscr ] bounding matrix for a given matrix is given. To our best knowledge, this was the first work to construct ≤icx-monotone bounding matrices. In fact, the construction of ≤icx bounding chains has been possible by means of the specific monotonicity characterization of class [Cscr ] Markov chains. The generalization of class [Cscr ] Markov chains will extend the chains that have closed-form solutions to compute transient and the steady-state distributions. Therefore, it would be possible to compute more accurate bounds by constructing monotone bounding chains that belong to this extended class.

The article is organized as follows. In Section 2 we give some preliminaries on the stochastic comparison approach. Section 3 is devoted to class [Cscr ][Gscr ] Markov chains. We first present main properties of this class and the closed-form solutions to compute transient distributions and the steady-state distribution. Then we give algorithms to construct ≤st and ≤icx bounding matrices belonging to this class. In Section 4 numerical results are presented in order to give some insights into how the accuracy of bounds might depend on the underlying matrix structure. We conclude in Section 5 and give some perspectives to extend this work.

2. PRELIMINARIES ON THE STOCHASTIC COMPARISON

Here we state the basic definitions and theorems on the stochastic comparison approach. We refer to [13,16] for further details. First, we give the definition to compare two random variables taking values on a totally ordered space [Escr ]. Let [Fscr ]st denote the class of all increasing real functions on [Escr ] and [Fscr ]icx denote the class of all increasing and convex real functions on [Escr ]. We denote by ≤[Fscr ] the stochastic order relation, where [Fscr ] can be replaced by st or icx to be associated respectively to the class of functions [Fscr ]st or [Fscr ]icx. In the sequel, all of the vectors and matrices are boldfaced, all of the vectors are row vectors, and xt denotes the transposed vector for x. When comparing two vectors or matrices, ≤ denotes the componentwise comparison.

Definition 1: Let X and Y be two random variables taking values on a totally ordered space [Escr ],

whenever the expectations exist.

Notice that for discrete random variables X and Y with probability vectors p and q, the notations p[Fscr ] q and X[Fscr ] Y are used interchangeably. Stochastic orders ≤st and ≤icx can be also defined by matrices (see [9,10]). Throughout this work, we assume that [Escr ] = {1,…, n}, but the following statements can be extended to the infinite case. We give Kst and Kicx matrices related respectively to the ≤st and ≤icx orders. In the sequel, K[Fscr ] denotes the matrix related to the ≤[Fscr ] order, [Fscr ] ∈ {st,icx}:

Proposition 1: Let X and Y be two random variables taking values on [Escr ] = {1,2,…, n} and let p = (pi)i=1n and q = (qi)i=1n, be probability vectors such that

which can be also given as follows:

It is shown in Theorem 5.2.11 of [13, p.186] that monotonicity and comparability of the probability transition matrices of time-homogeneous Markov chains yield sufficient conditions to compare the underlying chains stochastically. We first define the monotonicity and comparability of stochastic matrices and then state this theorem.

Definition 2: Let P be a stochastic matrix. P is said to be stochastically[Fscr ]-monotone if for any probability vectors p and q,

Definition 3: Let P and Q be two stochastic matrices. Q is said to be an upper-bounding matrix of P in the sense of the[Fscr ] order (P[Fscr ] Q) if

Let us remark that this is equivalent to saying that P[Fscr ] Q if

where Pi,* denotes the ith row of matrix P.

Theorem 1: Let P (resp. Q) be the probability transition matrix of the time-homogeneous Markov chain {Xn,n ≥ 0} (resp. {Yn,n ≥ 0}). If

  • X0[Fscr ] Y0,
  • at least one of the probability transition matrices is monotone (i.e., either P or Q is[Fscr ]-monotone),
  • the transition matrices are comparable (i.e., P[Fscr ] Q),

then

The following lemma [13] lets us compare the steady-state distributions of Markov chains, which can be considered as the limiting case (ΠP = limn→∞ X(n)).

Lemma 1: Let Q be a monotone, upper-bounding matrix for a finite stochastic matrix P in the sense of the[Fscr ] order. If the steady-state distributions (πP, πQ) exist, then

3. GENERALIZED CLASS [Cscr ] ([Cscr ][Gscr ]) MARKOV CHAINS

In this section we give the generalization of class [Cscr ] Markov chains that will be called class [Cscr ][Gscr ] Markov chains. First we present definitions and main properties of this class; then we give algorithms to construct bounding chains.

3.1. Class [Cscr ][Gscr ] Markov Chains and Closed-Form Solutions

Definition 4 (Class [Cscr ][Gscr ]): A stochastic matrix P belongs to [Cscr ][Gscr ] if there are three vectors v, c, and r in

such that v is a probability vector (vi ≥ 0, ∀i,

,

, and

where e = (ei)i=1n and ei = 1, ∀i.

Remark 1: Class [Cscr ] Markov chains [3] can be obtained by taking r = (ri)i=1n, ri = i − 1, ∀i. In this case, v equals P1,*.

In the following lemma we give the properties of power matrices for class [Cscr ][Gscr ] matrices.

Lemma 2: Let P be a stochastic matrix from class [Cscr ][Gscr ]: P = etv + rtc. Let α and β be the following constants:

The power matrices Pk,k ≥ 2, belong to class [Cscr ][Gscr ] and they can easily be computed as follows:

Proof: The proof is done by induction on k.

Basic step (k = 2): P2 = (etv + rtc)2 = etvetv + etvrtc + rtcetv + rtcrtc. Using the fact that vet = 1 and cet = 0 (Definition 4), crt = α and vrt = β (Eq. (4)), we have: P2 = etv + [αrt + βet] c.

Induction step: Suppose that the lemma holds for k and let us show that it holds also for k + 1. Let

. Hence, Pk = etv + ftc, so

On the other hand,

and

. Therefore,

. █

We can now state the theorem on the closed-form solution for transient distributions of this class of Markov chains.

Theorem 2: Let {X(k), k ≥ 0} be a time-homogeneous discrete-time Markov chain with probability transition matrix P and let πk = (πik)i=1n be the distribution of X(k).

If P belongs to class [Cscr ][Gscr ], then for all k ≥ 0,

where ak is the constant defined as

The constants α and β have been already defined in Eq. (4) and g is defined as

Proof: Using Lemma 2, we have:

The following corollary gives the conditions to have the limiting distribution of [Cscr ][Gscr ] matrices.

Corollary 1 (Steady-state distribution for class [Cscr ][Gscr ] matrices): Let {X(k), k ≥ 0} be a time-homogeneous discrete-time Markov chain with probability transition matrix P from class [Cscr ][Gscr ] and let α be defined as α = crt. {X(k), k ≥ 0} has a stationary distribution if and only if |α| < 1. In that case, the unique stationary distribution (steady-state distribution) π = limk→∞ πk is given by

where R is a constant given by R = vrt/(1 − crt).

Proof: From Theorem 2 we have πk = v + ak c, ∀k, where πk is the distribution of X(k) and

. For |α| < 1, limk→∞ ak = β/(1 − α), so we have π = v + (β/(1 − α))c = v + Rc, which is independent of π0. For |α| ≥ 1, ak diverges; thus, the stationary distribution does not exist. █

Let us now give the conditions to have the limiting distribution in terms of entries of P belonging to class [Cscr ][Gscr ].

Remark 2: Since

, we have α ≥ −1. The condition |α| < 1 is satisfied if and only if

.

It follows from Theorem 2 and Corollary 1 that the computation of transient distributions and the steady-state distribution are linear with the size of the model.

Remark 3: The complexity of computing πk is θ(max{n,k}) and it equals θ(n) for π.

3.2. Monotonicity of Class [Cscr ][Gscr ] Matrices

The following proposition gives sufficient conditions for the monotonicity of class [Cscr ][Gscr ] matrices in terms of vectors c and r. For a vector

, we will write x ∈ [Fscr ] if and only if vector x, seen as a real function on [Escr ] = {1,…, n}, is in class [Fscr ].

Proposition 2: A matrix P = etv + rtc ∈ [Cscr ][Gscr ] is[Fscr ]-monotone if

Proof: First, we define vectors c′ = cK[Fscr ] and v′ = vK[Fscr ]. Hence, PK[Fscr ] = etv′ + rtc′. To prove that P is ≤[Fscr ]-monotone, we consider two probability vectors p and q such that p[Fscr ] q. According to Definition 2 we must show that pP[Fscr ] qP, which is equivalent to pPK[Fscr ]qPK[Fscr ] (see Proposition 1). We have

On the other hand, if X and Y denote random variables corresponding respectively to probability vectors p and q, then using Definition 1 and the fact that r ∈ [Fscr ], we have E [r(X)] ≤ E [r(Y)]. Remark that

and, similarly, E [r(Y)] = qrt, so prtqrt. Using in addition the fact that c′ = cK[Fscr ] ≥ 0, we conclude that prtc′ ≤ qrtc′, and the proof is achieved. █

Remark 4: For the ≤st order, r ∈ [Fscr ]st if and only if vector r is increasing; that is,

For the ≤icx order, r ∈ [Fscr ]icx if and only if vector r is increasing and convex; that is,

In the general case, when the underlying matrix P does not belong to class [Cscr ][Gscr ], the monotonicity characterization is given as follows:

  • P is ≤st-monotone if and only if for each column j, the function
    is increasing with respect to i [9].
  • P is ≤icx-monotone if and only if for each column j, the function
    is increasing and convex with respect to i (see [10] for the infinite state space case and [2] for the finite state space case).

3.3. Bounding Algorithms

Obviously, having closed-form solutions makes [Cscr ][Gscr ] matrices attractive for computational issues. However, for a given matrix, the constraints to belong to [Cscr ][Gscr ] matrices might not be satisfied. We propose to construct class [Cscr ][Gscr ] bounding matrices by applying the stochastic comparison approach. In fact, we apply Theorem 1: For a given P, we construct Q ∈ [Cscr ][Gscr ] such that P[Fscr ] Q and Q is ≤[Fscr ]-monotone. Therefore, we can compute the transient distributions and the steady-state distribution of Q through its closed-form solutions to provide bounding distributions on P.

3.3.1. Strong stochastic bounds.

In this subsection we give the algorithm to construct bounding matrix in the sense of the ≤st order. As stated earlier, a [Cscr ][Gscr ] stochastic matrix is defined through three vectors: v, r, and c. For a given matrix P and a nonnegative, increasing vector r, Algorithm 1 returns vectors c and v such that Q = etv + rtc ∈ [Cscr ][Gscr ]; Pst Q and Q is ≤st-monotone. The nonnegativity constraint is used only to simplify the algorithm and its proof. Some possible choices for r will be discussed in Section 3.4.

Algorithm 1:st-Monotone [Cscr ][Gscr ] Upper Bound Q of P

Input: P: stochastic matrix; r: nonnegative, r ∈ [Fscr ]st (increasing)

Output: Q: ≤st-monotone, [Cscr ][Gscr ] upper bound of P, Q = etv + rtc

Notation: K := min{i | ri > 0},M := min{i | ri = rn}

We suppose that M ≠ 1. (In the other case, the algorithm is trivial: Taking the maximal value of each column of PKst, we obtain a vector x. We can then take v = xKst−1; Q = etv.)

For M ≠ 1, K is well defined and we have KM.

//Column n:

1.

;

2.

;

for j = n − 1 downto 2 do

3.

;

4.

;

end

//Column 1:

5.

;

6.

;

We give the following theorem to prove Algorithm 1.

Theorem 3: Let P be a stochastic matrix and Q be the matrix computed by Algorithm 1. The matrix Q is a stochastic matrix, Q ∈ [Cscr ][Gscr ], Q isst-monotone and Pst Q.

The proof of Theorem 3 is given in the Appendix.

3.3.2. Increasing convex bounds.

In Algorithm 2 we give the construction of ≤icx bounding matrices. Let us first define the operators that will be used in this algorithm.

and

are defined as follows:

Let us remark that φicx(x,j) = (xKicx)j and

.

Algorithm 2:icx-Monotone [Cscr ][Gscr ] Upper Bound Q of P

Input: P: stochastic matrix; r: nonnegative vector, r ∈ [Fscr ]icx (increasing, convex)

Output: Q: ≤icx-monotone [Cscr ][Gscr ] upper bound of P, Q = etv + rtc

Notation: K := min{i | ri > 0},M := min{i | ri = rn}

We suppose that M ≠ 1. (In the other case, the algorithm is trivial: Taking the maximal value of each column of PKicx, we obtain a vector x. We can then take v = xKicx−1; Q = etv.)

For M ≠ 1, K is well defined and we have KM.

//Column n:

1.

;

2.

;

for: j = n − 1 downto 2 do

3.

;

4.

;

end

//Column 1:

5.

;

6.

;

The following theorem is given to prove Algorithm 2.

Theorem 4: Let P be a stochastic matrix and Q be the matrix computed by Algorithm 2. The matrix Q is a stochastic matrix from class [Cscr ][Gscr ] and it is anicx-monotone upper bound for P.

The proof of Theorem 4 is given in the Appendix.

Remark 5: The computation complexity of Algorithm 1 and Algorithm 2 is θ(n2).

3.4. Choosing Vector r

For both algorithms, vector r of row increments, satisfying monotonicity conditions (see Proposition 2), is fixed. Moreover, vector r is also supposed to be nonnegative in order to simplify the algorithm and its proof. Thus, any vector satisfying the above conditions can be chosen as the vector of row increments. We present some of possible choices for vector r in the following.

If we take vector r = (ri)i=1n, ri = i − 1, we obtain the bounding matrix that belongs to class [Cscr ] (see Remark 1). It is important to notice that this choice is independent of the underlying matrix entries.

We propose some other possibilities for vector r in order to compute more accurate bounds. The main idea here is to improve bounding entries in the last column. In fact, since in both algorithms the construction is done recursively from the last column to the first column, the perturbations in the last column entries of the bounding matrix are propagated to other bounding matrix entries. So by diminishing the perturbations in the last column, it might be possible to improve the accuracy of the bounds. The other reason comes essentially from performance evaluation studies. Quite often, we are interested in bounding only the tail distribution. For instance, in reliability applications we are interested in bounding the probability that the studied system is not operational; in queuing applications we bound the overflow probability of the buffer. In general, these measures concern a few states and they can be ordered as last states. Thus, accurate bounds for last states would be sufficient for these cases.

3.4.1. Strong stochastic case.

For the strong stochastic order, there is a vector r that minimizes the differences between the exact and the bounding entries for the last column.

Proposition 3: Let the input vector r of Algorithm 1 be defined as

Then the bounding matrix Q computed by Algorithm 1 satisfies the following:

Proof: Pst U and Pst Q (see Theorem 3) yield that pi,nui,n and pi,nqi,n, ∀i, so Eq. (9) is equivalent to qi,nui,n, ∀i. Since U is ≤st-monotone, ui,n are increasing in i. Moreover, ui,npi,n, ∀i, so ui,n ≥ maxki pk,n = ri, ∀i.

Finally, we will show that Algorithm 1 with the input vector r defined by Eq. (9) yields vn = 0 and cn = 1. Therefore, qi,n = vn + ri cn = riui,n, ∀i, and the proof is concluded. Indeed, we have rnri > 0, 1 ≤ i < M, ri = maxki pk,ipi,n, ∀i, and rn ≤ 1; thus, it can be seen from the first line of Algorithm 1 that vn = 0. On the other hand, pi,nri, ∀i, and rK = pK,n so it follows from the second line of Algorithm 1 that cn = 1. █

In the sequel, vector r obtained by Eq. (8) will be called the St-LC vector. In the following example, we present class [Cscr ][Gscr ] matrices obtained with vector r corresponding to class [Cscr ] and with the St-LC vector computed from Eq. (8).

Example 1: Let P be the original matrix with the steady-state distribution

. By taking vector r that corresponds to class [Cscr ], r = (0, 1, 2), Algorithm 1 returns matrix A with vector v = (0.4,0.6,0) and vector c = (−0.2, −0.2, 0.4). The steady-state distribution of A is πA = (0.1, 0.3, 0.6).

The St-LC vector computed from Eq. (8) is r = (0, 0.4, 0.6). In this case, Algorithm 1 returns bounding matrix A′ with vector v = (0, 0.4, 0.6) and vector c = (−0.5, −0.5, 1). The steady-state distribution of this bounding matrix is πA = (0.2, 0.4, 0.4).

We note that A′ belongs to class [Cscr ][Gscr ] but not to class [Cscr ] and πPst πAst πA. However, we do not always have this relation (see Fig. 4).

Steady-state analysis (≤st order).

3.4.2. Increasing convex case.

In the case of the increasing convex order, for a given matrix P it is not always possible to find a monotone upper-bounding matrix Q satisfying the following conditions:

Let us show this for matrix P of Example 1. It can be shown that the following matrices B and B′ are ≤icx-monotone upper-bounding matrices for P. To satisfy Eq. (10) we must have pi,nqi,n ≤ min(bi,n,bi,n′), ∀i; thus, q1,3 = 0, q2,3 = 0.4, and q3,3 = 0.6. However, this vector r = (0, 0.4, 0.6) is not compatible with the monotonicity constraints (see Remark 4).

Let us recall that an ≤icx-monotone upper-bounding matrix Q of P is optimal in the sense of ≤icx order if

In fact, Eq. (10) represents the optimality conditions for the last column, so it is a necessary condition for the global optimality [Eq. (11)]. Hence, the previous example shows that it is not always possible to find the optimal bound in the ≤icx sense. Moreover, this is the case for any stochastic matrix, not only for class [Cscr ][Gscr ] matrices, since the properties of this class have not been considered.

As we cannot always satisfy Eq. (10), we propose to minimize the sum of distances between qi,n and pi,n:

Proposition 4: Let vector r be a solution of the following problem: Find vector x with a minimal sum of entries satisfying the monotonicity (x ∈ [Fscr ]icx) and the comparison constraints,

Then the upper-bounding matrix Q computed by Algorithm 2 with this input vector r minimizes Eq. (12); that is,

Proof: Let us remark first that Picx U and Picx Q (see Theorem 4) give pi,nui,n and pi,nqi,n, ∀i; so Eq. (14) is equivalent to

. Constraints (13) are fulfilled by all ≤icx-monotone upper bounds of P, so we have

, ∀Uicx-monotone upper bound of P.

Similarly to the proof of Proposition 3, it can be shown from Algorithm 2 that vn = 0 and cn ≤ 1. Moreover, if cn < 1, then

, which is in contradiction with the fact that r is a solution of (13). Thus, we have cn = 1. █

The Simplex algorithm can be used to find a solution of problem (13). Let us remark that solution of this problem is not necessarily unique. The bounds computed from Algorithm 2 with input vector r is obtained through the Simplex algorithm will be called Icx-Simplex bounds.

Although problem (13) can be efficiently solved by the Simplex algorithm, we cannot guarantee its complexity in the worst case. Therefore, we present an alternative heuristic choice for r. Unrolling the comparison [Eq. (2)] and the monotonicity (Remark 4) constraints for the last column increasingly in rows,

we obtain a nonnegative vector r such that r ∈ [Fscr ]icx. In the sequel, the bounds obtained from Algorithm 2 with this input vector will be called Icx-LC bounds.

Let us remark that a vector r obtained by relations (15) satisfy constraints (13) except the last one, rn ≤ 1. However, this will be compensated by the fact that we will have cn < 1 during the construction of the bounding matrix through Algorithm 2. Thus, the last column of the bounding matrix Q computed by Algorithm 2 satisfies all of the constraints (13). Although the Icx-LC bound does not necessarily minimize Eq. (12), since we take into consideration the information contained in the last column of the original matrix, this choice could yield better bounds than class [Cscr ] bounds for rewards that are mostly based on the last states. This is true in the case of the example studied in Section 4. Moreover, for this example, Icx-LC bounds are remarkably close to Icx-Simplex bounds.

Remark 6: The complexity of computing different vectors r is not important compared to that of Algorithms 1 and 2 (Remark 5). In fact, Eqs. (8) and (15) can be solved with linear complexity, and Eq. (13) can be solved efficiently with the Simplex algorithm.

4. NUMERICAL EXAMPLE

Let us first emphasize that the proposed bounding method can be applied to any finite discrete-time Markov chain, but the tightness of bounds depends largely on the underlying transition matrices. Our main goal in this section is to give some insights on the relationship between the accuracy of [Cscr ][Gscr ] bounds and the underlying matrix structures. We study an example that has been also studied in previous articles because of its simplicity in generating different structures of probability transition matrices.

We consider a finite-capacity buffer with multiple servers and batch arrivals. The time is supposed to be divided into slots and the analysis is provided in discrete time. The number of servers is denoted by S and the buffer capacity by B; thus, C = B + S is the maximum number of customers in the system (see Fig. 1). The size of batch arrivals is assumed to be distributed according to a truncated modified geometrical distribution with parameter pa. The probability that i customers arrive during a slot is ai = (1 − pa)pai, 1 ≤ iA, where A is the maximal batch size. The service time is distributed geometrically, and ps denotes the probability that the service is finished in one slot. Obviously, the values of S and A have a large impact on the sparseness of the matrix. The number of nonzero entries increases with the values of S and A and the probability transition matrix is full when A = S = B. The probabilities pa and ps for arrivals and services, respectively, determine where the probability mass is located in the matrix. In the following figures the curves of bounds are noted according to vector r used as input in bounding algorithms (see Section 3.4).

Considered system.

In Figure 2, the values of B, pa, and ps are fixed and we assume that S = A. We vary this last value and give the average number of customers in the steady state. We note that when S = A increases, the accuracy of bounds is improved for all class [Cscr ][Gscr ] bounds for both ≤st and ≤icx orders. In fact, whatever the structure of the original matrix, the bounding matrix is filled due to the structure of class [Cscr ][Gscr ] matrices. Thus, if the number of zero entries in the original matrix is large, the number of modified entries would be also large and this might cause a deterioration of the accuracy of bounds. However, we note that it is possible to have tight bounds even for relatively sparse matrices, especially with the ≤icx order. In fact, the presented results confirm the intuition which follows from the structure of [Cscr ][Gscr ] matrices that the proposed bounds are more adapted for full matrices. On the other hand, it is not possible to state that the bounds for relatively sparse matrices are always loose, since the quality of bounds depends also on the magnitude of the perturbations.

Impact of sparseness of the transition matrix on the quality of class [Cscr ][Gscr ] bounds (steady state).

In the case of the ≤st order, we can see in Figure 2 that the St-LC and the St-Class C curves cross each other: The class [Cscr ] bound is tighter for large values of A = S (full matrices) and the St-LC bound is tighter for sparser matrices. Thus, each of these ≤st bounds might be of interest. In the case of the ≤icx order, class [Cscr ] bounds are tighter than the other considered class [Cscr ][Gscr ] bounds and they are very close to the exact values for full matrices. Moreover, there is not a sensible difference between Icx-LC and Icx-Simplex bounds. This phenomena has been observed not only in Figure 2 but also for different sets of parameters that we considered. Hence, the heuristic proposed in Section 3.4.2 seems to be useful to avoid the complexity of the Simplex algorithm.

We also note in Figure 2 that ≤icx bounds are generally more accurate than the ≤st bounds. This is due to specific monotonicity properties of class [Cscr ][Gscr ] matrices. In fact, the ≤st comparison of random variables and stochastic matrices always implies the ≤icx comparison [13], as the class of increasing convex functions is contained in the class of increasing functions. For the monotonicity, this is not generally the case. However, for class [Cscr ] matrices, the monotonicity in the ≤st sense implies the monotonicity in the ≤icx sense [4]. We can also show that, for [Cscr ][Gscr ] matrices, if the chosen vector r is increasing and convex (r ∈ [Fscr ]icx), then cKst ≥ 0 implies the ≤icx monotonicity (see Proposition 2). On the other hand, ≤st bounds might be useful if one wants to compute bounds on increasing but not convex rewards. For instance, if we want to compare the tail distributions, we cannot use ≤icx bounds.

Let us remember that the proposed approach can be also applied to provide transient bounds. In Figure 3 we present transient behaviors of the average number of customers for the parameter values of Figure 2 with A = S = 10. We can see here that in the transient case we can make the same observations as in the steady-state case (Figure 2). More generally, we have noticed also through other experiments that the quality of transient bounds has the same tendency as in the case of steady-state bounds.

Transient [Cscr ][Gscr ] bounds for the case A = S = 10.

In Figures 4 and 5, we consider bounds on steady-state distributions. Let N be the number of customers in the steady-state case. In Figure 4, we present ≤st bounds for the tail distribution (Prob(Nx)). For ≤icx bounds, it is not possible to compare tail distributions. However, we have the following characterization of the ≤icx order (see [13]):

Steady-state analysis (≤icx order).

Informally speaking, this characterization says that, for a given level x, we compare the mean overflow from that level. In Figure 5, we compare E [(Nx + 1)+] for different ≤icx bounds. We note that for this set of parameters, the bounds obtained using vectors r derived from the last column are more accurate for both ≤icx and ≤st bounds. For instance, in Figures 4 and 5, for x ≥ 10 we obtain tighter bounds by taking into account the last column of the original matrix.

It is not surprising to obtain more accurate bounds for last states (i.e., bounds on the rewards of last states), by taking into account the information contained in the last column of the original matrix. In fact the information used to derive the input vector r is valid for the last column and also, to a lesser extent, for the columns that are relatively close to it. However, it is less intuitive that we obtain tighter results from class [Cscr ] ≤icx bounds when studying average behavior of the underlying model.

It can be seen from Figures 6 and 7 that different class [Cscr ][Gscr ] bounds can cross each other several times. Due to the low computational complexity of class [Cscr ][Gscr ] bounds, it is worth computing different class [Cscr ][Gscr ] bounds and taking the best one for each value. In Table 1 we give the best ≤st and ≤icx bounds that can be derived from the bounds of Figures 6 and 7. In Table 2 we present, for the same set of parameters, the best transient bounds for Prob(N ≥ 15) for the ≤st order and E [(N − 14)+] for the ≤icx order. We observe that the best transient bounds are obtained from different vectors r, as in the steady-state case.

Bounds can cross more than once (≤st case).

Bounds can cross more than once (≤icx case).

Best ≤st and ≤icx Bounds Derived from Bounds of Figures 6 and 7

Best ≤st and ≤icx Transient Bounds for the Same Set of Parameters as in Figures 6 and 7

5. CONCLUSION

In this article we present the generalization of class [Cscr ] matrices to class [Cscr ][Gscr ] matrices. This class of Markov chains has closed-form solutions to compute transient distributions and the steady-state distribution, which makes it attractive for computational issues of large Markov chains. We propose to apply the algorithmic stochastic comparison approach to construct upper-bounding matrices belonging to this class. For a given finite Markov chain we construct a bounding class [Cscr ][Gscr ] chain and compute its distributions to provide bounding distributions. By analyzing bounding class [Cscr ][Gscr ] chains, we significantly gain in terms of both computational and storage complexity, compared to the exact analysis.

We give algorithms to construct class [Cscr ][Gscr ], monotone, upper-bounding matrices in the sense of the ≤st and the ≤icx orders. Class [Cscr ][Gscr ] matrices are completely defined by means of three vectors: v, c, and r. Our method consists of two steps. First, we use the information contained in the last column of the original transition matrix to compute vector r of row increments. Once vector r is fixed, we apply Algorithm 1 (≤st bounds) or Algorithm 2 (≤icx bounds) to determine the remaining v and c vectors.

We present different choices of vector r to construct class [Cscr ][Gscr ] bounding matrices in the sense of both ≤st and ≤icx orders. By changing input vector r, it is possible to easily construct different bounding matrices in class [Cscr ][Gscr ]. These matrices can be analyzed through the closed-form solutions and the better bound can be taken for a given performance measure of interest.

We have performed numerical experiments to study the accuracy of bounds as a function of matrix structures and different choices of vector r. We observe that class [Cscr ][Gscr ] bounds become tighter when the sparseness of the underlying matrix decreases and that ≤icx bounds are tighter in general than ≤st bounds. We show that ≤st bounds can be improved by class [Cscr ][Gscr ] matrices compared to class [Cscr ] matrices. This is especially useful when we are interested in functionals of distributions that are increasing but not convex. For ≤icx bounds, we conclude that class [Cscr ] bounds are generally accurate. Other class [Cscr ][Gscr ] bounds could be interesting when computing bounds for rewards defined on last states.

Computing closed-form stochastic bounds seems to be a promising approach to overcome the state space explosion problem of Markov chains. We are considering extending our work to look for other matrix structures having closed-form solutions and to other stochastic orders.

APPENDIX

In this appendix we give a common proof for Theorems 3 and 4. For that, we will first give some definitions to have a common notation for Algorithms 1 and 2 and related properties. In the sequel, ≤[Fscr ] will denote either the ≤st order or the ≤icx order. It corresponds to the ≤st order in the context of Algorithm 1 and to the the ≤icx order in the context of Algorithm 2.

Let φst and

be the operators

and

, respectively, defined by

Note that lines 1, 2, 5, and 6 are the same in Algorithms 1 and 2. In addition, as Qi,* = v + ri c, lines 3 and 4 can be written in a common way as follows:

3.

;

4.

.

On the other hand, both conditions r ∈ [Fscr ]st and r ∈ [Fscr ]icx, corresponding respectively to ≤st and ≤icx orders, imply that r is increasing. Moreover, the input vector r for both algorithms is nonnegative. Note also that the operators φ[Fscr ] and

are linear in the first component and that

for both orders ≤st and ≤icx. However, we do not have the same relation between φ[Fscr ](x,j) and φ[Fscr ](x,j + 1), which is why we separate the cases ≤st and ≤icx for some points of the proof of Theorems 3 and 4. Let us remark that

For the sake of clarity of this proof, we give first the following two lemmas.

Lemma A1: Let P be a stochastic matrix. The matrix Q = etv + rtc computed by Algorithm 1 or by Algorithm 2 satisfies

Proof:

Column j = n: We have φ[Fscr ](Pi,*,n) = pi,n and φ[Fscr ](Qi,*,n) = qi,n. For 1 ≤ i < K, as ri = 0, from line 1 it follows that qi,n = vn + ri cn = vnpi,n. For Kin, line 2 and r nonnegative ⇒ qi,n = vn + ri cnpi,n.

Column j, 2 ≤ jn − 1: For i < K, ri = 0 and line 3 of the algorithm

. For Kin, line 4

.

Column j = 1: Lines 5 and 6 imply that

. Thus, the case of the ≤st order is trivial

. For the ≤icx order, using (A2) we have

. █

Lemma A2: Let P be a stochastic matrix and Q = etv + rtc the matrix computed by Algorithm 1 or by Algorithm 2. Then we have the following:

Proof: cet = 0 follows directly from line 6 of the algorithm. Line 2 of the algorithm assures cn ≥ 0 and line 4 gives

; hence, φ[Fscr ](c,j) ≥ 0, ∀j > 1. For j = 1, as

, φst(c,1) = 0 and φicx(c,1) = φicx(c,2) ≥ 0. We conclude that φ[Fscr ](c,j) ≥ 0, ∀j (i.e., cK[Fscr ] ≥ 0).

vet = 1 follows directly from line 5. By the construction of v (lines 1 and 3) we have vj ≥ 0, ∀j > 1. It remains us to show that

(i.e.,

). In order to show this inequality, we will show that

by induction on j, 2 ≤ jn. Since pi,n ≤ 1, 1 ≤ i < M, from line 1 of the algorithm it follows that vn ≤ 1.

Let us suppose that

. From line 3, if vj = 0, then

. Otherwise,

In the case of the ≤st order,

, since

and rn > 0.

In the case of the ≤icx order, as φicx(Pi,*,j + 1) − φicx(Qi,*,j + 1) ≤ 0 (Lemma A1),

Proof of Theorems 3 and 4: From (A4) and the definition of Q as Q = etv + rtc, we have Q ∈ [Cscr ][Gscr ]. Lemma A2 (relation A5) and Lemma A1 give that Q is monotone and an upper bound for P with respect to the order ≤[Fscr ]. It remains just to show that Q is a stochastic matrix.

From Lemma A2 we have Qi,* et = vet + ri cet = 1, ∀i. In order to show that Q is stochastic, we must show that qi,j ≥ 0, ∀i,j.

For the last column we have vn ≥ 0 and cn ≥ 0 (lines 1 and 2 of the algorithm), so since r is nonnegative, qi,n = vn + ri cnvn ≥ 0, ∀i.

From line 4 of the algorithm we have cj ≥ −(vj /rn) ⇒ qn,j = vj + cj rn ≥ 0, 2 ≤ jn − 1. As the vector r is increasing, all of the columns of Q are monotone (increasing for cj ≥ 0 and decreasing for cj < 0) and contained between two values that are both nonnegative (vj and qn,j); thus, qi,j ≥ 0, ∀i, 2 ≤ jn − 1.

For the first column we have

, ∀i, so we must show that

, ∀i. If

, then

, ∀i. In the other case (

, as r is increasing, we have

, ∀i, so it is sufficient to show that

. We will show it by induction on j. Let us remark that ri = rn for iM. For j = n from line 1 of the algorithm we have

Let us suppose that

. We deduce that

On the other hand, note that

and

(Lemma A2). Hence,

. As rn > 0 and

(from (A4)), we have

From line 3 of the algorithm it follows that

For Ki < M, ri > 0, so we have

For Min, ri = rn. In the case of the ≤icx order, using the fact that φicx(Pi,*,j + 1) ≤ φicx(Qi,*,j + 1) = φicx(v + ri c,j + 1) (Lemma A1) and relation (A2), we have

The last inequality is obviously satisfied in the case of the ≤st order, since

and

. We conclude that

From (A6)–(A9) it follows that

, so

. By induction on j it follows that

, so we proved that Q is a stochastic matrix. █

References

REFERENCES

Abu-amsha, O. & Vincent, J.M. (1998). An algorithm to bound functionals of Markov chains with large state space. In 4th INFORMS Conference on Telecommunications, Florida.
Ben Mamoun, M. (2002). Encadrements stochastiques et évaluation de performances des réseaux. PhD thesis, University of Versailles.
Ben Mamoun, M. & Pekergin, N. (2000). Computing closed-form stochastic bounds on the stationary distribution of Markov chains. In D. Gardy & A. Mokkadem (eds.), Mathematics and computer science: Algorithms, trees, combinatorics and probabilities. Basel: Birkhauser-Verlag, pp. 197209.
Ben Mamoun, M. & Pekergin, N. (2002). Closed-form stochastic bounds on the stationary distribution of Markov chains. Probability in the Engineering and Informational Sciences 16: 403426.Google Scholar
Ben Mamoun, M. & Pekergin, N. (2005). Computing closed-form stochastic bounds on transient distributions of Markov chains. In M. Papazoglu & K. Yamazaki (eds.), Applications and the Internet SAINT2005 workshops. New York: IEEE Computer Society, pp. 260263.
Courtois, P.J. & Semal, P. (1984). Bounds for the positive eigenvectors of nonnegative matrices and their approximation by decomposition. Journal of the ACM 31(4): 804825.Google Scholar
Fourneau, J.M. & Pekergin, N. (2002). An algorithmic approach to stochastic bounds. In M. Calzrossa & S. Tucci (eds.), Performance evaluation of complex systems: Techniques and tools. Berlin: Springer-Verlag, pp. 6489.
Golubchik, L. & Lui, J. (1997). Bounding of performance measures for a threshold-based queuing systems with hysteresis. In Proceeding of ACM SIGMETRICS'97, pp. 147157.
Keilson, J. & Kester, A. (1977). Monotone matrices and monotone Markov processes. Stochastic Processes and their Applications 5: 231241.Google Scholar
Li, H. & Shaked, M. (1994). Stochastic convexity and concavity of Markov processes. Mathematics in Operations Research 19(2): 477493.Google Scholar
Lui, J., Muntz, R., & Towsley, D. (1998). Computing performance bounds of Fork–Join parallel programs under a multiprocessing environment. IEEE Transactions on Parallel and Distributed Systems 9(3): 295311.Google Scholar
Mahevas, S. & Rubino, G. (2001). Bound computation of dependability and performance measures. IEEE Transactions on Computers 50(5): 399413.Google Scholar
Muller, A. & Stoyan, D. (2002). Comparison methods for stochastic models and risks. New York: Wiley.
Muntz, R.R. & Lui, J. (1994). Computing bounds on steady-state availability of repairable computer systems. Journal of the ACM 41(4): 676707.Google Scholar
Pekergin, N. (1999). Stochastic performance bounds by state space reduction. Performance Evaluation 36–37, 117.
Shaked, M. & Shantikumar, J.G. (1994). Stochastic orders and their applications. San Diego: Academic Press.
Stewart, W.J. (1994). Introduction to the numerical solution of Markov chains. Princeton, NJ: Princeton University Press.
van Dijk, N.M. (1998). Bounds and error bounds for queueuing networks. Annals of Operations Research 79: 295319.Google Scholar
van Dijk, N.M. & van der Wal, J. (1989). Simple bounds and monotonicity results for finite multi-server exponential tandem queues. Queueing Systems 4: 116.Google Scholar
Figure 0

Steady-state analysis (≤st order).

Figure 1

Considered system.

Figure 2

Impact of sparseness of the transition matrix on the quality of class [Cscr ][Gscr ] bounds (steady state).

Figure 3

Transient [Cscr ][Gscr ] bounds for the case A = S = 10.

Figure 4

Steady-state analysis (≤icx order).

Figure 5

Bounds can cross more than once (≤st case).

Figure 6

Bounds can cross more than once (≤icx case).

Figure 7

Best ≤st and ≤icx Bounds Derived from Bounds of Figures 6 and 7

Figure 8

Best ≤st and ≤icx Transient Bounds for the Same Set of Parameters as in Figures 6 and 7