{'fraction_non_alphanumeric': 0.040475509815132454, 'fraction_numerical': 0.008552506194015628, 'mean_word_length': 3.9718109676655216, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perfect Sampling of Markov Chains with Piecewise Homogeneous Events

Perfect sampling is a technique that uses coupling arguments to provide a sample from the stationary distribution of a Markov chain in a finite time without ever computing the distribution. This technique is very efficient if all the events in the system have monotonicity property. However, in the general (non-monotone) case, this technique needs to consider the whole state space, which limits its application only to chains with a state space of small cardinality. We propose here a new approach for the general case that only needs to consider two trajectories. Instead of the original chain, we use two bounding processes (envelopes) and we show that, whenever they couple, one obtains a sample under the stationary distribution of the original chain. We show that this new approach is particularly effective when the state space can be partitioned into pieces where envelopes can be easily computed. We further show that most Markovian queueing networks have this property and we propose efficient algorithms for some of them.


Introduction
Simulation approaches can be efficient ways to estimate the stationary behavior of Markov chains by providing independent samples distributed according to their stationary distributions, even when it is impossible to compute this distribution numerically.
Propp and Wilson used a backward coupling scheme [14] to derive a simulation algorithm -called PSA (Perfect Sampling Algorithm) in the following -providing perfect samples (i.e. whose distribution is stationary) of the state of a discrete time Markov chain over a finite state space S. The main idea is to simulate trajectories starting from all states x ∈ S at some time in the past until time t = 0. If the state at time 0 is the same for all trajectories, then the chain has coalesced and the final state has the stationary distribution of the Markov chain. Otherwise, the simulations are started further in the past. Perfect sampling procedures have been developed in various contexts (for more information, see the annotated bibliography [16]).
Although the Perfect Sampling Algorithm provides perfect samples for all irreducible and aperiodic finite Markov chains, its complexity can be high. It suffers from two drawbacks that jeopardize its applicability for very large chains.
• The first one is the fact that the coupling time can be very large. Some recent work focused on the estimation of the coupling time for certain classes of Markov chains. For example, it was shown in [4] that Markov chains, modeling a class of networks of queues with finite capacities, have a quadratic coupling time with respect to capacities.
• The second factor in the complexity of PSA is the fact that one needs to run one simulation per state in S, which is prohibitive for most applications. Various techniques have been developed to reduce the number of trajectories that need to be considered in the coupling from the past procedure. A first crucial observation already appears in [14]: for a monotone Markov chain, one has to consider only extremal initial conditions. For anti-monotone systems, an analogous technique that also considers only extremal initial conditions has been developed by Kendall [11] and Häggström and Nelander [8].
To cope with non-monotonicity, Kendall and Møller introduced in [12] the general idea of two bounding processes sandwiching all the normal trajectories, reducing the number of processes having to be computed to two. They applied the sandwiching technique to perfect sampling of point processes, together with the idea of a dominating chain that allowed to handle the infinite state space. The idea of sandwiching processes was also applied in communication networks [1,7]. Huber [9] also introduced a similar idea of bounding chains for determining when coupling has occurred. In his model, the bounding chains are evolving in some bigger state space than the target Markov chain. However, both Kendall-Møller's and Huber's constructions are modeldependent and do not have an algorithmic counterpart so that they are not straightforward to apply in general.
In this paper we introduce a new perfect sampling algorithm (EPSA: Envelope Perfect Sampling Algorithm) that uses the lattice structure of the state space to design an effective and automatic construction of a bounding interval chain that can be expressed using only two trajectories. This is done by computing from the past, a lower and an upper bound of all trajectories of the Markov chain. Whenever the two bounds meet, the coupling state is distributed according to the stationary distribution of the chain. The idea of bounding interval chain is close to sandwiching processes in Kendall and Møller [12]. However, their approach was not constructive and the sandwiching processes were computed ad-hoc for each application.
Our aim is to propose, under a constructive approach, several algorithms to compute bounding interval chains that are easily computable and as tight as possible. This is done in two cases. The first (easy) case is when events are almost state space homogeneous (defined in Section 4). In that case, we provide an algorithm that computes the update rule of the bounding interval chain. The second case is more difficult since events are assumed to be piecewise defined (Section 5) over arbitrary polytopes (Section 6). In that case we show how to obtain an interval containing all trajectories by solving a set of linear programs. Since the worst case complexity of linear programs may be too high, we show how to avoid solving these, at the expense of obtaining bounding intervals that may be less tight.
In the last part of the paper (Section 7), we consider applications to queueing networks. We exhibit several events (such as negative customers, fork and join, batch routing) that satisfy almost state space homogeneity. We also provide numerical evidence of the very good performance of EPSA in these cases, by computing the coupling time experimentally.
We also show examples (such as generalized join the shortest queue) that are piecewise homogeneous, where the smallest interval containing all trajectories can also be computed in linear time in the number of queues. Our numerical experiments show that EPSA behaves just as well as in the easy case of almost space homogeneous events.

Definitions and Basic Properties
If not stated otherwise, all real (or integer) vectors are row vectors and v t denotes the transposed vector. For x, y ∈ R d , we denote x ∨ y def = (max{x 1 , y 1 }, . . . , To avoid multiple brackets, for x, y, z ∈ R d , we will write: Note that (x ∨ y) ∧ z = x ∨ (y ∧ z); the relations involving operators ∨ and ∧ will be evaluated from left to right. For any x ∈ R n , ||x|| 1 = n i=1 |x i |, and for any set U , |U | denotes its cardinality.
Throughout the paper, we consider discrete time homogeneous Markov chains (DTMC) with a finite state space S. Since we only consider finite state spaces, the continuous time models can be easily handled using the same techniques, after uniformization.

Markov Automaton
The evolution of a finite DTMC can always be obtained using a finite number of discrete events (or actions). To formalize this, we introduce the notion of Markov automaton. Markov automatons are similar to Generalized Semi Markov Processes [3] with a focus on state changes rather than on time.
S is a finite state space, A is an alphabet called the set of events, D is a probability distribution on A, and · is a right action by the letters of A on the states in S and is called the transition function : Equivalently, this action can be most naturally extended to a right action by words u 1→n def = u 1 . . . u n in A n on S for any n ∈ N: A Markov automaton A = (S, A, D, ·) naturally induces a Markov chain with the following construction: Let (u 1 , . . . , u n , . . . ) be an infinite sequence of random letter of A i.i.d. distributed according to D. We consider this sequence as an infinite word u def = u 1 . . . u n . . . on A. Then for any x 0 ∈ S, the random process (X n def = x 0 · u 1→n ) n∈N is a Markov chain issued from x 0 with probability transition matrix P given by for all x, y in S, P (x, y) = a∈A x·a=y P D (a). (2.1) We will say that the Markov chain (X n ) is generated by A and u. Moreover, with this construction, we can build a family of Markov chains {(X n (x) = x · u 1→n ) n≥0 | x ∈ S} starting from each state x ∈ S. We will refer to this family of Markov chain as the grand coupling generated by A and u. We will say that the grand coupling has coalesced at time n if all the Markov chains of the family has reached the same state, i.e. when S · u 1→n is reduced to a singleton. Conversely, for any probability transition matrix P on a finite space state S, it is easy to see that there exists a Markov automaton A = (S, A, D, ·) such that (2.1) holds, i.e. such that A generates a Markov chain on S with same transition matrix P , but that automaton is in general not unique. However as we will see in our examples, such representation of Markovian systems with Markov automatons naturally arises for many systems such as queueing networks.

Perfect Sampling
Let {X n } n∈N be an irreducible and aperiodic DTMC with a finite state space S and a transition matrix P . Let π denote the steady state distribution of the chain: π = πP .
The Perfect Sampling Algorithm (PSA) was first used by Propp and Wilson. It allows one to draw in finite time steady state distributed variable using coupling from the past.
The grand coupling generated by A and u couples if = 1. In that case, let τ def = inf n : S · u −n+1→0 = 1 be the coupling time of the chain.
Then E(τ ) < ∞ and S · u −τ +1→0 is steady state distributed. The algorithmic counterpart (PSA) of Theorem 2.2 is given in Algorithm 1. The main drawback of PSA algorithm is the fact that one needs to simulate one Markov chain starting from each state in S, that could be too large for a practical use of the algorithm.
Several approaches have been used to overcome this problem. The main one for a partially ordered state space (S, ) and monotone events is already present in [14].
3. An event a ∈ A is said to be monotone if x y =⇒ x · a y · a. An event a ∈ A is said to be anti-monotone if x y =⇒ x · a y · a.
If all events are monotone, then it is possible to generate a steady state by starting Algorithm 1 with maximal and minimal states only [14]. This technique has been successfully used in [15] to construct PSA for network of queues. Similarly, an efficient perfect sampling algorithm for the chains with anti-monotone events was given in [8].
When events are neither monotone nor anti-monotone, one can still use monotone bounds, as in [5], but this technique provides a sample from a bounding chain and not the initial one. We give in this paper an algorithm for perfect sampling of general (non-monotone) chains.

Envelopes
The approach proposed here generalizes what has been done in [6] and [2] to sample non-monotone Markov chains. It uses a bounding chain method similar to the one employed by Kendall and Møller in [12]. The sandwiching chains in [12] are in general difficult to compute, except for some models for which effective ad-hoc rules can be derived. Our approach is more algorithmic. We introduce a general framework that provides efficient algorithms for a wide variety of events, including many queueing events.

Bounding Interval Chains
We assume from now on that the state space S is equipped with a lattice order relation . We define bounding interval chains as follows. In particular we notice that when m n = M n , the grand coupling inside has necessary coalesced. Thus, by constructing efficiently such bounding intervals, we can detect easily if a grand coupling has coalesced and apply it to PSA.
We now construct the tightest bounding interval chain computable step by step for a given grand coupling, that we will call envelope chain.
inf U sup U  The construction of is illustrated on a simple example in Figure 2. As with '·', we extend the operator to a finite word of events u 1→n = u 1 . . . u n ∈ A n : It is also the tightest bounding chain computable step by step, in the sense that if, at time n, we only know that S · u 1→n is in [m n , M n ], then the tightest interval containing S · u 1→n+1 we can compute at time n + 1 is [m n , M n ] u n+1 . However, it is not tight in the stronger sense of [m n , M n ] being equal to S · u 1→n , that does not hold in general.
Notice that the lower envelope (m n ) n∈N alone is not a Markov chain, neither is the upper one (M n ) n∈N , since they depend on each other.
We can then use the envelope process to detect coupling in PSA. Proof. The proof simply uses the fact that the set S · u −n+1→0 is included in [⊥, ] u −n+1→0 . Consequently, if the latter is reduced to one point, so is the former and the grand coupling has coalesced.
Algorithm 1 can be adapted to the envelope simulation: start the simulation with only two states, ⊥ and , and iterate using and an i.i.d. sequence of events. This new algorithm, called EPSA (Envelope Perfect Sampling Algorithm), is given in Algorithm 2. The reason why we double n at each loop of the algorithm is that we need to compute S u −n+1→0 in each loop, which corresponds to n iterations of . While increasing n by 1 at each loop would lead to a quadratic cost in n, doubling it keeps the complexity linear. EPSA is a generalization of both monotone and anti-monotone sampling algorithm (see [14,8]): • Monotone events: if an event a ∈ A is monotone, then for any m M , (3.2) • Anti-monotone events: if an event a ∈ A is anti-monotone, then for In the next section we discuss the complexity of EPSA, that depends on the coupling time of the envelope chain, but also on the complexity of the computation of the envelope transition. The efficiency of the (anti-) monotone case is due to the fact that when ' · ' is (anti-)monotone, the complexity of is the same as using ' · ' twice.
Note however that the construction of the envelopes depends on the representation of the Markov chain by a Markov automaton that is not unique. Different automatons representing the same chain would lead to different envelope chains with different coupling properties (one may coalesce almost surely and the other not, or if they both coalesce their coupling times may be different). The complexity of the envelope transition may also differ depending on the representation chosen.

Complexity of EPSA
The envelope approach may not gain over the general PSA because of two drawbacks: The assumption that [m n , M n ] hits the set of single point intervals X in Theorem 3.2 may not hold. By doing so, we would still obtain a bounding interval chain, however its coupling time might be longer than for the envelope chain, or the bounding interval chain might not couple at all. There is then an obvious tradeoff between computing simple large bounds quickly with an increased coupling time, or computing tight bounds slowly with a potentially shorter coupling time.
Remark 3.4. Actually, even if the envelope chain does not couple, it is still possible to obtain perfect samples in finite time by using the splitting algorithm [2]. The splitting algorithm is hybrid: it first runs EPSA when the envelopes are too far apart and switches to the usual PSA algorithm as soon as the number of states inside the envelopes becomes manageable.
The rest of this paper focuses on improving the cost C e of computing envelopes under several assumptions except for Section 4.2 and numerical examples in Section 7 that discuss the coupling time issue.
Indeed, if C e is not too large compared to C (as it is the case for the various transitions proposed in this article), then the comparison of the two methods is essentially a comparison between the coupling time of EPSA (E(τ e )) and the coupling time of PSA multiplied by the cardinality of the state space (|S| × E(τ )).

Almost Space Homogeneous Events (ASHEs)
We now present a family of simple events called Almost Space Homogeneous Events (ASHEs) for which the corresponding envelope transitions are easy to compute. These events will be the first brick to build more elaborate ones. Intuitively, these events consist in moving homogeneously every point of the space in the same direction, with some simple conditions on the boundaries of the state space.

Computing the envelopes for ASHEs
We assume from now on that S is an interval of Definition 4.1. An event a is an almost space homogeneous event (ASHE) if there exists a vector v ∈ Z n and a binary relation R (called blocking relation) on the set of components {1, . . . , d} such that for all x in S we have the following property: (4.1) We call v the direction vector of event a.
Remark 4.2. The blocking relation R needs only to be specified on Act def = {i | v i = 0} the set of active components of the transition, since the event has an effect only along these components.
The cases B(x) = ∅ or B(x) = d are interesting particular cases that correspond respectively to the case where, if x + v ∈ S, then the event is truncated at the boundary of the state space: Here is a more complex case in the context of queues. Let us consider a Markovian queueing network made of 3 finite queues and exponential service times and a fork node: one departing packet from queue 1 is duplicated into two packets joining queues 2 and 3. This event is an ASHE, corresponds to adding a vector v = (−1, +1, +1) to all states except on the boundaries. The critical cases are when the first queue is empty or when one of the next queues is full. If the first queue is empty, then the event is blocked so that (1, 2)    The proof is given in Appendix A.

Nonexpansive properties of ASHEs
In the previous Section, we have shown that computing the envelopes for ASHEs can be done efficiently, so that the term C e in the complexity is in O(d 2 ). This section gives some insights on the coupling time. Since our approach is general, no tight bound on the coupling can be derived. The most immediate property is that, as mentioned in Section 3, the coupling time of envelopes is larger than the coupling time of the original chain. However, as seen in the experimental section 7, the difference between two coupling times is not large. This can be explained in part by several nonexpansive properties of ASHEs. In this section we show that the size of the envelope begin interval cannot increase more than linearly in the worse case, and does not increase (or even contracts) under some conditions on the blocking relation.
In the latter case, the contracting property also gives some insight of the splitting condition. Assume that a ∈ A is an ASHE given by (4.1). We first give a useful lemma about the set of blocked components of an ASHE. Proof. We first prove the inclusion ⊆ : Let x ∈ U and j ∈ B. Then Proof. Since there are no interactions between components, we can consider them separately. Thus, for i ∈ B: M i − m i ≤ M i − m i , as in Proposition 4.6. For i ∈ B, we have three different cases: There are two cases: The above result tells us that an event without blocking interactions between different components (i.e. (i, j) ∈ R =⇒ i = j) is nonexpansive in components that are not auto-blocking. For components that are autoblocking such event is nonexpansive when M i − m i ≥ |v i | − 1.
Suppose that for at least one component i ∈ A, we have two states x and y such that 5 Partition of the state space 5

.1 Zone partition
In this section we focus on a single event a. The state space is partitioned into zones, Z 1 , · · · , Z K , where K is assumed to be small compared to the size of the state space. The action a may have a different effect a k on each zone Z k . Formally, the action '·' by a is defined for all x ∈ S as : x·a = x·a k where k is such that x ∈ Z k . To compute the envelopes of the action of a over an arbitrary interval [m, M ], one can use the partition into zones.
The construction of these envelopes is illustrated in Figure 4. The first and second equalities are direct consequence of the definition of the partition; the third one is by definition of the action on a set; the fourth one simply uses the fact that A ∪ B = A ∪ B .
The last inclusion (5.2) does not hold with equality in general. It uses the fact that for any set A, A · a ⊂ A a. It gives us an easy way to compute an over-approximation of . Moreover, in some cases the inclusion can be replaced by an equality, such that the over-approximation is exact (or tight). We give some simple cases where this equality holds in Section 5.2.
The computation of an interval containing [m, M ] a using (5.2) reduces to two problems: (P 1 ) The computation of [m, M ] ∩ Z k . This can be done by integer linear programming (see next section).
(P 2 ) The computation of [m, M ] ∩ Z k a k . This is easy when a k is an ASHE (using Algorithm 3) or when a k is monotone (see Equation (3.2)).
In Section 6, we will give methods to solve (P 1 ) by computing [m, M ] ∩ Z k or an over-approximation. We can now prove that U a ⊆ U · a . Let [m , M ] = U a and [q, Q] = U · a . To prove the inclusion we just need to show that m ≥ q and M ≤ Q.
Corollary 5.2. Let a be a piecewise ASHE over a partition Z 1 , . . . , Z K such that each ASHE a k is non-blocking on its zone Z k . Then the inclusion (5.2) can be replaced by an equality.

Remark 5.3.
It is easy to see that any ASHE a can be partitioned into zones such that a is non-blocking on each zone. Although this is not useful for computing the envelope of a simple ASHE (since we already know how to compute it), it can however be used to obtain tight envelopes on piecewise ASHEs.
Proof. The inclusion U · a ⊆ U a is always true, and was already proved. We keep the same notation as in the previous proof: On component i, we need to consider several cases (we will show that m i ≥ q i and M i ≤ Q i ): , then the transition is non-blocking on U , and the previous Proposition applies. We now prove that • If i ∈ B(M ) \ B(m) a similar reasoning shows the same result.
We finally obtain the reverse inclusion U a ⊆ U · a .

Polytopic Zones and Linear Programming
In this section we consider a class of piecewise events for which the computation of envelopes can be done in a reasonable amount of time and memory space. We will focus here on one event a assumed to be defined piecewise over polytopic zones. Again, the state space is a product S = [0, C 1 ] × · · · × [0, C d ], and is here partitioned into zones defined by linear constraints.
Let us consider H hyperplanes of R d in general position and assume they do not intersect S (by translating them by a small ε if needed). These hyperplanes define a partition of R d in K convex polytopes (P 1 , ..., P K ), and these polytopes induce a partition of the state space S in K zones Z i = S ∩ P i . Without loss of generality, we consider only the nonempty zones and assume we have K of them.
Finally, over each zone Z k , the event a k is assumed to be easily computable (e.g. an ASHE).
There exists a relation between the total number of hyperplanes H needed to define the zones and the maximal number of zones that can be defined using them. The number K of zones is smaller than 1 + H + H 2 + · · ·+ H d . This is because the state space S may intersect all zones (bounded or not) formed by H hyperplanes in a general position. The number of such zones is given by the value of the characteristic polynomial of the hyperplanes at point -1 (see [17]). Example 6.1. A queueing example is displayed in Figure 5. We consider a Join the Shortest Waiting time (JSW) event, that is an arrival that allocates the new customer to the queue with the shortest expected waiting time (after arrival). The zones for the JSW event are displayed in Figure 6. If (x 2 + 1)/µ 2 > (x 1 + 1)/(2µ 1 ) + ε (with ε very small) then the new customer is sent to queue 1 as long as the queue is not full, in which case it is sent to queue 2. If (x 2 + 1)/µ 2 ≤ (x 1 + 1)/(2µ 1 ) + ε, it is sent to queue 2 as long as the queue is not full, in which case it is sent to queue 1. If both queues are full, the client is rejected. Here 3 hyperplanes divide the two-dimensional state space into 4 zones that are polytopes: (x 2 + 1)/µ 2 = (x 1 + 1)/(2µ 1 ) + ε, Note that in this case, the partition makes the event a non-blocking ASHE in each zone, so that Equation (5.2) holds with equality (by Proposition 5.1). However, the two zones in the upper right corner could be merged in one single zone with a blocking ASHE, and the equality would still hold (see Proposition 5.4).

Zone intersection with Linear Programming
Since we assumed the event a to be easily computable on each zone Z k , the only problem left to solve is (P 1 ) the computation of [m , M ] = [m, M ] ∩ Z k . When the zones Z k are induced by polytopes P k (Z k = S ∩ P k ), this computation amounts to solving 2d integer linear programs, one for each component m i and M i : Figure 5: Example of a queueing event that can be defined with ASHEs on polytopic zones. The first queue has two identical servers with rate µ 1 while the second has only one server with rate µ 2 . ∀i ∈ {1, . . . , d}, Unfortunately, integer linear programs are well known to be NP-hard in general. These linear programs can be solved in simple cases that we will not detail in this paper. A general approach to overcome this difficulty is to relax the integrity condition and to solve the corresponding rational linear programs, by computing q k and Q k in Z d defined by It is easy to see that the discrete interval [q k , Q k ] contains [m, M ] ∩ Z k . We can thus compute an over-approximation of for a :  Figure 7 gives an illustration on a simple example.

Complexity and data structure
In order to deal with the set of linear programs in an efficient way, they can be represented as a binary acyclic graph. Each node corresponds to an additional linear constraint. The left child corresponds to the half plane where the inequality is negative and the right child to the one where the inequality is positive. When an intersection can be reached through several paths, the nodes are merged.
The leaves of the graph are all the polytopes and the number of nodes in the tree is less than twice the number of leaves. This construction is illustrated in Figure 8.
This representation has a better amortized complexity to compute [m, M ]·a than an exhaustive resolution for each polytopic zone since the intersection of [m, M ] with a zone can be found empty early in the tree (this can even happen at the root of the tree if [m, M ] does not intersect one half plane. Sorting the constraints (from the root of the tree to the leaves) to get the best possible expected cost of computing [m, M ] · a can be done for a given distribution of m and M (for example uniform over S × S) by using dynamic programming. This is not in the scope of this paper.
The Algorithm 4 defines the procedure Bound(V, m, M, a) such that, when V is the graph G, Bound (G, m, M, a) = K k=1 [q k , Q k ] a k which is an over-approximation of [m, M ] a. It uses the procedure LP-Solve (not detailed) that solves linear programs over rational numbers in polynomial time. It also contains the test V ∩ [m, M ] Q = ∅ that can be done much faster (in linear time, see the next section) and helps to cut the cost: The LP-solve procedure is not used when the intersection is empty.
By solving linear programs with the interior point method [10], the complexity of each main call to Algorithm 4 is O(d 4.5 K).

Fast bounds using Minkowski sums
The previous cost to compute over-approximations of envelopes by solving linear programs may not be sustainable for chains with a large dimension d and when the number of calls to Bound is large.
In this section we will derive a way to compute looser bounds on the envelopes, with a complexity O(dHK). This technique is based on the computation of Minkowski sums of all K zones with intervals [m, M ] (see [13]).
Each polytope P k can be expressed as the intersection of h k hyperplanes: where A k is a d×h k matrix and b k a row vector of Algorithm 4: Recursive procedure computing the bounds on the envelopes using linear programming Data: Partition of the state space into Z 1 , . . . , Z K arranged in a graph structure G using linear constraints;  Figure 9 shows an illustration of the Minkowski sum in a simple case. It is direct to see that and that the set W k is a polytopic zone defined by h k + 2d hyperplanes: ∀i ∈ {1, . . . , d}, Indeed, since P k ⊂ [ k , L k ] Q , we have where k (resp. L k ) is the component-wise ceiling (resp. floor) of k (resp. L k ) and is pre-computed. The intervals [m ∨ k , M ∧ L k ] can replace the intervals [q k , Q k ] in Algorithm 4, and are easier to compute. The complexity of Algorithm 4 would then be O(dHK). However, these intervals provide looser bounds on the envelopes and might thus increase the coupling time of the bounding processes, at worst to infinity. As we pointed out earlier, this tradeoff between coupling time and complexity of computations needs to be adjusted depending on the Markov chain to simulate.

Queueing Network Examples
In this section, we illustrate the usefulness and the expressive power of piecewise ASHEs by showing examples of queueing networks that can be described using Markov automatons with piecewise ASHEs over polytopic zones.
The rule of thumb is that many classical (and some non-classical) queueing networks are such chains so that sampling their stationary distribution is made possible using ASHES over polytopic zones, while the classical PSA would be intractable in practice.

Jackson Network with finite capacities
We consider open queueing systems similar to Jackson networks except that queues have finite capacities. The network is formed by d queues numbered from 1 to d: Each queue i is a M/M/1/C i queue with finite capacity C i and a FIFO server. Entering customers arrive in each server i according to a Poisson process of parameter λ i , and each server is assumed to have exponential service time µ i . After being served at queue i, a customer is directed to queue j with probability P i,j or leaves the network with probability We will consider two types of blocking policy: • Client Loss (CL): If a client served by queue i is sent to a full queue, it is lost.
• Restart Service (RS): If a client served by queue i is sent to a full queue, it restarts its service at queue i.
We assume that new clients arriving in the network from outside are lost when sent to a full queue. The state of the network is entirely described by the number of customers in each queue, and the corresponding Markov chain can be generated with the following Markov automaton A = (S, A, D, ·) where S = [0, where λ a the frequency rate of each event a is defined as: and where the a i,j are ASHEs with direction vector v a i,j = e j − e i for all i, j ∈ {0, . . . , d} (with the convention e 0 = (0, . . . , 0)), and with a blocking relation R a i,j depending on the blocking policy used: R CL a i,j = {(i, j)} for the Client Loss policy and R RS a i,j = {(i, j), (j, i)} for the Service Restart policy. Note that when i = 0 or j = 0, the choice of R has no incidence. The blocking policy may also differ for each couple i, j.
It is easy to check that the events of this Markov automaton are all monotone for the product order on S = [0, C 1 ] × [0, C d ]. Therefore for this type of networks, the coupling times for PSA and EPSA are equal.

Examples of non-monotone ASHEs
In this section, we consider other types of routing events that are not monotone nor anti-monotone and can be described with ASHEs. These events can easily be added to any Jackson networks with finite capacity, and they can easily be combined with each other to give other types of ASHEs.
Envelopes for these events were already computed in [2] using ad-hoc methods, while the ASHE framework now embraces them all. The numerical results for these events are taken from [2].
Fork and Join A common type of non-monotone event used in queueing network is fork and join. A fork queue is such that, at service completion, several customers are sent simultaneously to several disjoint output queues. A join is a queue with several input buffers. The server is only active when no buffer is empty. A service completion removes one customer in each input buffer and sends one (merged) customer to an output buffer.
Fork events were already introduced in Example 4.3 and shown to be ASHEs. A join event with two input buffers i and j and output buffer k is described by an ASHE with direction vector v = e k − e i − e j and blocking relation depending on the blocking policy: Join events are not monotone for the natural product order: Consider for instance the join event a with i = 1, j = 2 and k = 3 in a 3-server network, and consider the two states x = (0, 1, 0) and y = (1, 1, 0). Then x < y while x · a = (0, 1, 0) and y · a = (0, 0, 1) are incomparable. Similarly, it is easy to show that fork is in general not monotone.
Negative customers We consider the additional events corresponding to negative customers. After queue i finishes a service, a negative customer is sent to queue j and kills one customer in that queue, if there is any. This event can be represented by an ASHE of direction vector v = −e i − e j and blocking relation R = {(i, j)}. Such events are not monotone for the natural product order.
We ran simulations for the small network displayed in Figure 11. The average coupling time of PSA and EPSA are reported in Figure 10 (in the following all data are given with 95 % confidence intervals).
In the example, the coupling time of EPSA is around 3 times larger than PSA. Since the number of states is 15 2 = 225, this means that EPSA is about 225/6 = 40 times faster than PSA in this case. Batch routing Here, we consider batch routing: After the end of service in queue i, a batch of K customers leaves queue i and creates L new customers in queue j. In general, K and L do not need to be equal: for instance 10 products that are assembled into one (e.g. put into a box), and then sent to next phase of a production line can be modeled by taking K = 10 and L = 1. This event can be represented by an ASHE of direction vector v = −Ke i + Le j and blocking relation depending on the blocking policy: As before, the exogenous arrivals (resp. services) can be modeled by setting i = 0 (resp. j = 0). We allow partial arrivals to queue j (in the case of R CL ), while the service must be complete (i.e. at least K customers must be present in queue i). However, partial service can also be modeled by piecewise ASHEs, using similar ideas as for multiple servers in Section 7.3. Figure 12 reports the coupling time of EPSA for a single queue with batch arrivals of size K = 1, 2, 3 or 5 with rate λ, departures of size 1 and rate 1, and blocking policy R = {(1, 1)}. Unlike what happens for negative customers and fork-join nodes, the coupling time of EPSA increases exponentially fast with the load of the queue (the scale being logarithmic). Such increase is due to the fact that envelopes can only couple after they both hit 0. The time for the upper envelope starting in C to hit 0 when the load ρ is larger than 1 becomes exponentially large (larger than e ρC ).
One way to prevent this is to use the Split algorithm [2, Algorithm 3], that is an hybrid between PSA and EPSA: It starts by computing bounding intervals as EPSA until the size of the interval becomes smaller than the size of the largest batch. When this happens, the bounding interval is then split into all the single trajectories contained in it, and these normal trajectories are then computed as in PSA. The advantage of the Split algorithm is that these normal trajectories as not as numerous as the |S| trajectories needed for the classical PSA, and it is guaranteed to terminate when PSA does unlike EPSA. Figure 13 reports the coupling times of EPSA, PSA and Split for a queue with capacity 20, arrival batches +2 and +3 with probabilities 0.49 and 0.51, and departures of size one (with rate 1). In this figure, the coupling time of PSA and Split are similar while EPSA's coupling time grows very fast when λ/µ is larger than one. In this case Split has a clear advantage over the other two algorithms: It couples exponentially faster than EPSA and deals with two envelopes instead of C trajectories for PSA. The efficiency of Split can be partially explained by the choice of the splitting point: as long as the size of the bounding interval is larger than K (the size of a batch), the batch arrival event is nonexpansive. This is a direct consequence of Proposition 4.8.

Examples of Piecewise ASHEs
We now consider even more complex queueing events that will prove to be piecewise ASHEs over polytopic zones.
Multiple Servers Let us assume now that queue i is a M/M/C/C queue with finite capacity C and C identical servers with serving rate µ, and that customers are sent to queue j after being served. Such type of service can be described using C piecewise ASHEs a 1 , . . . , a C as illustrated in Figure 14.
The state space of queue i is [0, C] and the transition rate matrix Q associated to this event is such that Q(x, x − 1) = xµ, the other coordinates being null. We add C events a 1 , . . . , a C with each a k defined as the following piecewise ASHE: on zone [0, k − 1], a k is the identity (v = 0, R = ∅) and on [k, C], a k is the ASHE with direction vector (v = e j − e i ) and blocking relation R = ∅. Finally, we set λ a k = µ, and take D(a k ) = λ a k / a∈A λ a .
This construction can be easily adapted to any queue with n < C servers, which can be represented with n piecewise ASHEs. It can also be easily combined with any other type of events discussed previously.
Join the Shortest Waiting Time (JSW) and other routing policies JSW is described in Example 6.1 and is shown to be a piecewise ASHE over polytopic zones.
More generally, consider the routing event r i where a client, after being served at queue i, is sent to queue j = R(x) for some function R : S → [0, d] (if j = 0 the client leaves the system), and the client is either lost (CL) or restarts its service (RS) if the destination queue is full. Such event can

Application to a network
Here is an example of a network where events are ASHEs or piecewise ASHEs over polytopic zones. The following model permits us to compare the routing policies JSW to the Random (Rand) routing.
In this model, we consider two parallel queues 1, 2 with same capacity C 1 = C 2 = 20 and one server with service rate µ 1 , µ 2 , respectively. We introduce µ = µ 1 + µ 2 and a = µ 1 /µ such that µ 1 = aµ and µ 2 = (1 − a)µ. Clients arrive in the system at a rate λ = 1 and go in one of the two queues, according to the chosen routing policies.
To compare the two policies, the 2-queue system is duplicated to give a 4-queue-network, with two 2-queue subsystems with coupled arrivals and departures (see Figure 15). Incoming clients are duplicated via a fork and sent to each subsystem. If one of the subsystem rejects its copy of the client, the other must reject it as well. This is to ensure that both systems always receive the same quantity of clients.  Figure 15: Network to compare the JSW and the Random routing policies The first system then routes its clients using the JSW policy, while the second system sends its clients randomly to queue 1 with probability √ µ 1 √ µ 1 + √ µ 2 and to queue 2 with remaining probability. If the packet arrives in a full queue, it is rejected (and the client is rejected by the other subsystem as well). The departures are also coupled: with rate µ 1 [resp. µ 2 ], a packet leaves queue 1 [resp. 2] (if the queue is not already empty), simultaneously for each subsystem. When the queue is empty in one subsystem, the corresponding queue in the other subsystem can still serve its clients. Figure 16 shows the results of the simulation of that network using the envelope technique. We plotted the average difference of clients in the queue between the two subsystems, when the ratio a varies between 0.1 and 0.9, and for several values of the network load ρ = λ/µ.
As we can see, the JSW routing strategy performs better than the random routing in all cases, and particularly when the serving rates are not close from each other or when the workload is high.
The coupling time is much lower for low load. In that case the dominant events are non-blocking services. If the load is high, the dominant events are arrivals that are blocking. Further, the coupling can only be achieved by hitting the lower boundary of the state space (by emptying the system), so the coupling time grows exponentially with the load.

Conclusion
Perfect sampling can be a very efficient method to study the stationary behavior of a Markov chain. This efficiency relies on three conditions. The first one concerns the small number of simulated trajectories needed for unbiased coupling. The second condition is the requirement that a small time should be needed to compute those trajectories and the third condition is to have a small coupling time. When these three conditions are satisfied, the perfect sampling technique is hard to beat to compute the stationary distributed with a given confidence interval.
In this paper we show that the first condition is easy to meet by using envelopes. We also show that the second condition can be met for piecewise space homogeneous chains by using a special type of events (ASHEs) and by over approximating the envelopes. The third condition is more chain dependent. We provide several examples where coupling using the envelope technique happens fast and others where the coupling time is very large. In those cases, the splitting technique can help reduce the coupling time without compromising with the other two points.