Perfect sampling of Jackson queueing networks

Abstract

We consider open Jackson networks with losses with mixed finite and infinite queues and analyze the efficiency of sampling from their exact stationary distribution. We show that perfect sampling is possible, although the underlying Markov chain may have an infinite state space. The main idea is to use a Jackson network with infinite buffers (that has a product form stationary distribution) to bound the number of initial conditions to be considered in the coupling from the past scheme. We also provide bounds on the sampling time of this new perfect sampling algorithm for acyclic or hyper-stable networks. These bounds show that the new algorithm is considerably more efficient than existing perfect samplers even in the case where all queues are finite. We illustrate this efficiency through numerical experiments. We also extend our approach to variable service times and non-monotone networks such as queueing networks with negative customers.

Appendix: Solution of the traffic equations

Let us consider the traffic equations for a general Jackson network:

$$\begin{aligned} \left\{ \begin{array}{llcl} \forall i\not = 0,&{}\lambda _i = \sum _{j=1}^M p_{ji}(\lambda _j \wedge \mu _j), \\ \lambda _0 = \mu _0. \\ \end{array}\quad \right. \end{aligned}$$

(17)

Theorem 5

The system (17) admits exactly one solution.

Proof

Let $h$ be the function $h(\mathbf {x}){:=}\mathbf {x'}, x'_i = \sum \limits _{j \in \mathcal {J}} p_{ji} (x_j \wedge \mu _j)$. $P$ only has non negative coefficients, so $h$ is clearly monotone on $\mathbb {R}^\mathcal {J}$.

We know that P is stochastic, so we have

$$\begin{aligned} h(\mathbf {x})_i \!=\!\! \sum \limits _{j\in \mathcal {J}} p_{ji} (x_j \wedge \mu _j) \!\le \!\! \sum \limits _{j\in \mathcal {J}} p_{ji} \max (\varvec{\mu }) \le \max (\varvec{\mu }) \Rightarrow h(x) \le \max \varvec{\mu }\;\text {for}\;\mathbf {x} \in \mathbb {R}^\mathcal {J}, \end{aligned}$$

giving an upper bound of the image of $h$. We can also notice that for all inputs, the outputs of $h$ are positive linear combinations of the inputs and of positive values ($\varvec{\mu }$). Therefore, if the inputs are non negative, then so are the outputs.

Hence the restriction of $h$ to the complete lattice $[0, \max \varvec{\mu }]^\mathcal {J}$ is a monotone function on the same lattice. We can apply the Knaster-Tarski theorem: its fix points form a non empty lattice.

The minimum ${\mathbf {x}}^{-}$, maximum ${\mathbf {x}}^+$ of this lattice, and their difference, $\varvec{\delta } \mathop {=}\limits ^{\mathrm{def}}{\mathbf {x}}^+ - {\mathbf {x}}^ -$, satisfy:

$$\begin{aligned} \delta _i = \sum \limits _{j \in \mathcal {J}, x^ +_j \le \mu _j} p_{ji} \delta _j + \sum \limits _{j\in \mathcal {J}, x^ -_j < \mu _j < x^ +_j} p_{ji}(\mu _j - x^ -_j) \le \sum \limits _{j \in \mathcal {J}} \delta _j p_{ji}, \quad \text { for all } i \ne 0. \end{aligned}$$

As for queue $0$, $\delta _0 =x^+_0 -x^-_0 =\mu _0 - \mu _0 = 0$. Since $P$ is stochastic and irreducible, its maximum eigenvalue is 1 and the corresponding eigenspace is of dimension one, generated by an eigenvector without null coordinates.

So $\varvec{\delta }$ can only be the vector $0$ meaning that the lattice is reduced to a singleton. Hence the solution of the semi linear system exists and is unique. $\square$

1.1 Computation

One could iterate the function $h$ until reaching the fixed point, but this method would not give any guarantee on the number of iterations. Instead, we propose working on a pair $({\mathcal {L}}, \mathbf {x})$, where ${\mathcal {L}}$ is a set of queues and $x$ is the solution of the following system; both system and solution will be updated together with the modification of ${\mathcal {L}}$.

$$\begin{aligned} \left\{ \begin{array}{llll} x_0 =\mu _0, \\ x_i = \sum _{j\in {\mathcal {L}}} \mu _j p_{ji} + \sum _{j \notin {\mathcal {L}}} (x_j\wedge \mu _j) p_{ji}&{} for\ i\ne 0. \end{array}\right. \end{aligned}$$

(18)

Starting with ${\mathcal {L}}$ equal to the set of all queues: ${\mathcal {L}} = {\mathcal {J}}$, we repeat the following procedure as long as possible.

• Compute the vector ${\mathbf {x}}$, solution of (18).

• Find one queue $i^ *$ in ${\mathcal {L}}$ such that $i^ *\in {\mathcal {L}}, x_{i^ *} \le \mu _{i^ *}$.

• set the new set ${\mathcal {L}} = {\mathcal {L}} \backslash \{ i^*\}$

Lemma 6

The solutions ${\mathbf {x}}$ of the successive systems form a non-increasing sequence and are expressible as linear combinations of $\varvec{\mu }_{|{\mathcal {L}}}$.

Proof

This proof is in two parts: we first prove that the solutions ${\mathbf {x}}$ form a non-increasing sequence.

Consider the family of functions $h_{{\mathcal {L}}}$ similar to the function $h$ from the proof of Lemma 5:

$$\begin{aligned} h_{{\mathcal {L}}} ({\mathbf {x}}) \mathop {=}\limits ^{\mathrm{def}}\sum \limits _{j \not \in {\mathcal {L}}} p_{ji} (x_j \wedge \mu _j) + \sum \limits _{j \in {\mathcal {L}}} p_{ji} \mu _j. \end{aligned}$$

Every function of this family is monotone and admits only one fixed point on the lattice $[0,\max \varvec{\mu }]^\mathcal {J}$ for exactly the same reason as in the proof of Theorem 5.

We also have immediately (using $\forall a, b, (a\wedge b) \le a$)

$$\begin{aligned} \forall {\mathcal {L}}, {\mathcal {L}}',\; {\mathcal {L}}\subset {\mathcal {L}}' \Rightarrow h_{\mathcal {L}}\preceq h_{{\mathcal {L}}'}. \end{aligned}$$

Using the characterization of the minimum fixed point as the lower bound of the set where the output is lower than or equal to the input, we obtain:

$$\begin{aligned} \forall {\mathcal {L}}, {\mathcal {L}}',\; {\mathcal {L}}\subset {\mathcal {L}}' \Rightarrow \forall {\mathbf {x}}, h_{\mathcal {L}}\preceq h_{{\mathcal {L}}'}\Rightarrow h_{{\mathcal {L}}}(\text {FP} (h_{\mathcal {L}}')) \le \text {FP} (h_{\mathcal {L}}') \Rightarrow \text {FP} (h_{\mathcal {L}}) \le \text {FP} (h_{{\mathcal {L}}'}). \end{aligned}$$

Since the sequence ${\mathcal {L}}$ is decreasing, the fixed points ${\mathbf {x}}$ of $h_{\mathcal {L}}$ form a non-increasing sequence.

As a corollary, for any integer $n$ and for any $i\in \mathcal {J}\backslash {\mathcal {L}}^n$, let $n_0$ be the step when $i$ was chosen. We have $x^n_i \le x^{n_0}_i \le \mu _i$. Every term $(x_i \wedge \mu _i), i \notin {\mathcal {L}}^n$ in the system generated by ${\mathcal {L}}^n$ takes the value $x_i$. Therefore, the system is equivalent to the new system:

$$\begin{aligned} \left\{ \begin{array}{llll} \forall i \not = 0,\; &{}x_i = \sum \limits _{j\in {\mathcal {L}}} \mu _j p_{ji} + \sum \limits _{j \notin {\mathcal {L}}} x_j p_{ji}, \\ &{}x_0 =\mu _0.\\ \end{array} \right. \end{aligned}$$

(19)

We denote by $S({\mathcal {L}}, \varvec{\mu })$ this new system, by ${\mathbf {x}}^n$ the solution of $S({\mathcal {L}}= {\mathcal {L}}^n,\varvec{\mu })$ and by $i_n$ the element chosen before the $(n+1)$th step.

Using the sequence $({\mathcal {L}}^n)_n$, we build a sequence of matrices $(D^{(n)})_{n\in {\mathbb {N}}}$ such that for all indexes n, for all vector of real $\varvec{\mu }$, ${\mathbf {x}}$ is solution of $S({\mathcal {L}}^n,\varvec{\mu })$ if and only if ${\mathbf {x}}= \varvec{\mu }D^{(n)}$ with the additional property of partial sub-stochasticity:

$$\begin{aligned} D^n\text { is }{\mathcal {L}}^n\text {-sub-stochastic if }\forall n, \forall i \in \mathcal {J}, \sum \limits _{j\in {\mathcal {L}}^n} d^n_{ij} \le 1. \end{aligned}$$

We take $D^{(0)} = P$; indeed, the solution ${\mathbf {x}}$ of $S({\mathcal {L}}^0=\mathcal {J},\varvec{\mu })$ satisfies ${\mathbf {x}}= \varvec{\mu }P$, and $P$ is stochastic.

Computation of $D^{(n+1)}$ from $D^{(n)}$: We want to use the vector $(\mu _1,\ldots , \mu _{i_n-1}, x_{i_n},\mu _{i_n+1},\ldots ,\mu _{N})$, where coordinate $i_n$ in $\varvec{\mu }$ is replaced by the variable $x_{i_n}$.

We obtain a new system of variables ${\mathbf {x}}$ and ${\mathbf {m}}$ and equations:

$$\begin{aligned} \left\{ \begin{array}{ll} {\mathbf {x}}\text { solution of } S({\mathcal {L}}^n, \mathbf {m}),\\ {\mathbf {m}}=(\mu _1,\ldots , \mu _{i_n-1},x_{i_n},\mu _{i_n+1},\ldots ,\mu _{N}). \end{array} \right. \end{aligned}$$

With an immediate substitution on $\mathbf {m}$, we can see that this system is exactly $S({\mathcal {L}}^{n+1}, \varvec{\mu })$. By the definition of $D^n$, this system is also equivalent to

$$\begin{aligned} \left\{ \begin{array}{ll} {\mathbf {x}}=\mathbf {m} D^{(n)}, \\ {\mathbf {m}}= (\mu _1,\ldots , \mu _{i_n-1},x_{i_n},\mu _{i_n+1},\ldots ,\mu _{N}), \end{array}\right. \end{aligned}$$

in which we can make the same substitution, therefore the solution ${\mathbf {x}}^{n+1}$ of $S({\mathcal {L}}^{n+1}, \varvec{\mu })$ is also a solution of $S({\mathcal {L}}^{n}, \mathbf {{\mathbf {m}}} )$, so that it must satisfy ${\mathbf {x}}^{n+1} = (\mu _1,\ldots ,$ $\mu _{i_n-1},x^{n+1}_{i_n},\mu _{i_n+1},\ldots ,\mu _{N}) D^n$.

We notice that $x^{n+1}_{i_n}$ is the only variable appearing on the right-hand side of this system, and satisfies:

$$\begin{aligned} x^{n+1}_{i_n}= \sum \limits _{j\in \mathcal {J}\backslash \{i_n\}} d^n_{ji_n} \mu _j + d^n_{i_ni_n} x^{n+1}_{i_n}. \end{aligned}$$

Hence⁴

$$\begin{aligned} x^{n+1}_{i_n}= \sum \limits _{j\in \mathcal {J}\backslash \{i_n\}} \frac{d^n_{ji_n}}{1-d^n_{i_ni_n}} \mu _j. \end{aligned}$$

Now, we substitute in every equation of the system the variable $x^{n+1}_{i_n}$ by this combination. We obtain exactly the system $S({\mathcal {L}}^{n+1},\varvec{\mu })$ in the form of a linear combination ${\mathbf {x}}^{n+1}= \varvec{\mu } D^{n+1}$.

To summarize, we obtain $D^{n+1}$ from $D^n$ by applying the following successive elementary operations:

$$\begin{aligned} d^{n+1}_{j i_n}&:= \frac{d^n_{j i_n}}{1-d^n_{i_ni_n}} \quad \text{ for } \text{ all }\; j\ne i_n\nonumber \\ d^{n+1}_{ij}&= d^n_{ij} + \frac{d^n_{i i_n} d^ n_{i_nj}}{1-d_{i_ni_n}} = d^n_{ij} + d^{n+1}_{i i_n} d^{n}_{i_nj} \quad \text{ for } \text{ all }\;i,j\ne i_n\nonumber \\ d^{n+1}_{i_n j}&:= 0 \quad \text{ for } \text{ all }\;j. \end{aligned}$$

(20)

Finally, let us show that $D^{n+1}$ is ${\mathcal {L}}$-sub-stochastic:

Consider row $i\ne i_n$ (the $i_n$th row is null). Using the previous result and the fact that $D^{(n)}$ is ${\mathcal {L}}^n$-sub-stochastic:

$$\begin{aligned} \sum \limits _{j\in {\mathcal {L}}^{n+1}} d^{n+1}_{ij} = \sum \limits _{j\in {\mathcal {L}}^n} d^{n}_{ij} + d^n_{ii_n}\left( \sum \limits _{j\in {\mathcal {L}}^n} d^{n+1}_{i_nj} -1\right) \le \sum \limits _{j\in {\mathcal {L}}^n} d^{n}_{ij} \le 1 \end{aligned}$$

Also, by using (20) and the fact that all coefficient of $d_{ij}^{n}$ are non negative, all coefficients $d_{ij}^{n+1}$ are clearly non negative as well. Therefore, $D^{(n+1)}$ is ${\mathcal {L}}^{n+1}$-sub-stochastic. $\square$

Lemma 7

The generated sequence stops after at most $|\mathcal {J}|$ steps and can only end on a solution of the traffic equations.

Proof

At each step, the size of ${\mathcal {L}}$ decreases by one. This size begins at $|\mathcal {J}|$ and cannot be negative. Hence the sequence cannot have more than $|\mathcal {J}|$ steps.

An ending is a state where every element of ${\mathcal {L}}$ satisfies: ${\mathbf {x}}^\text {end}_i > \mu _i\; \text {for}\; i\in {\mathcal {L}}$ (stopping condition) and ${\mathbf {x}}^\text {end}_i \le \mu _i\; \text {for}\; i\in \mathcal {J}\backslash {\mathcal {L}}$ (part of recurrence property).

Hence the $\min$ operator $\wedge$ in (7) gives $\mu _i$ on exactly the set ${\mathcal {L}}$. Using the function $h$ defined in the proof of Theorem 5 in the Appendix, we have

$$\begin{aligned} \forall i \in \mathcal {J},\, h({\mathbf {x}}^\text {end})_i = \sum \limits _{j\in {\mathcal {L}}} p_{ji}\mu _i +\sum \limits _{j\notin {\mathcal {L}}} p_{ji}{\mathbf {x}}^\text {end}_i = {\mathbf {x}}^\text {end}_i \end{aligned}$$

by definition of ${\mathbf {x}}^\text {end}$ as the solution of exactly this system. Since ${\mathbf {x}}^\text {end}$ is a fixed point of $h$, it is also the only solution of the system (7). $\square$

Remark 3

We don’t need to implement the instruction $\{\forall j \in \mathcal {J},{d}_{ji^*}\leftarrow 0\}$ since the coefficients in columns outside ${\mathcal {L}}$ will never be accessed again.

Proposition 4

(complexity) Algorithm 3 computes the solution of the traffic equations in $O(M^3)$ operations.

Proof

There will be at most $M$ executions of the loop.

Cost of the different operations:

• 1. initialization: The first step is a product of a matrix and a vector: $O(M^2)$.

• 2. choice of $i_n$: At most $|{\mathcal {L}}|$ tests are done in one set: $O(M)$.

• 3. update of $D$: At each step less than $M^2$ coefficients are updated, and each update consists in a constant number of operations.

• 4. update of ${\mathcal {L}}$: With a direct pointer, this operation is done in constant time.

• 5. computation of $\mathbf {x}$ It is a matrix-vector product: done in $O(M^ 2)$.

Finally, the algorithm is made of at most $M$ repetitions of this $O(M^2)$ loop, so its execution is in $O(M^3)$. $\square$