2.1. A hierarchical formulation of queueing systems

Within a hierarchical multilevel formulation, rates in $\boldsymbol{\lambda}$ have a distribution (or image) $\mathbb{P}_{\boldsymbol{\lambda}}\equiv \boldsymbol{\lambda}_*\mathbb{P}$ under a reference measure $\mathbb{P}$ on $(\Omega,\mathcal{F})$ . We assume this to admit a density $f_{\boldsymbol{\lambda}}$ with respect to a base measure that will further induce (by properties of exponential transitions) distributions over the service rates $\boldsymbol{\mu}$ and routing topology. Next, note that a network trajectory over a finite interval is a piecewise deterministic jump process, such that $X\equiv (\textbf{t},\textbf{x})$ is represented by a sequence of transition times $\textbf{t}$ along with states $\textbf{x}$ . Each pair $(\textbf{t},\textbf{x})$ is furthermore a random variable on a measurable space $(\mathcal{X},\Sigma_\mathcal{X})$ supporting finite $\mathcal{S}$ -valued trajectories, and a conditional density $f_{X|\boldsymbol{\lambda}}$ may be defined with respect to a dominating base measure $\mu_\mathcal{X}$ , such that the regular conditional probability $\mathbb{P}(A|\boldsymbol{\lambda})$ , $A\in\mathcal{F}$ , satisfies

\begin{equation*}\mathbb{P}(X^{-1}(B)|\boldsymbol{\lambda}) = \int_B f_{X|\boldsymbol{\lambda}}(\textbf{t},\textbf{x}) \,\mu_\mathcal{X}(d\textbf{t},d\textbf{x})\end{equation*}

for all $B\in\Sigma_{\mathcal{X}}$ (see Appendix A for details). In this case,

(2) \begin{align}f_{X|\boldsymbol{\lambda}}(\textbf{t},\textbf{x}) = \pi(x_0)\, e^{Q_{x_{I}}(T-t_I)}\, \prod_{i=1}^I Q_{x_{i-1},x_{i}}\, e^{Q_{x_{i-1}}(t_i-t_{i-1})} \end{align}

for every pair of ordered times $\textbf{t}=\{0,t_1,\dots,t_I\}$ in [0, T] and states $\textbf{x}=\{x_0,\dots,x_I\}$ . Here, $\pi(\cdot)$ denotes an arbitrary distribution over initial states, and $Q\equiv Q(\boldsymbol{\lambda})$ is the matrix of infinitesimal rates associated with fixed values in $\boldsymbol{\lambda}$ . The QN model is thus fully parametrized by a collection of hyperparameters. Analogous modelling choices for continuous-time Markov chains or Markov jump processes can be found in [Reference Baele, van de Peer and Vansteelandt2, Reference Huelsenbeck, Bollback and Levine13] or [Reference Zhao32], to name a few.

Figure 2: Closed QN with a single FCFS service station and a delay.

2.2. Network evaluation and problem statement

Let $T>0$ denote some arbitrary terminal time and $\textbf{x}_0\in\mathcal{S}$ an initial state in X. For simplicity, this is assumed to be a 0-valued vector, where no jobs populate the system. Now, let $0\leq t_1<\dots<t_K\leq T$ denote some fixed network monitoring times along with observation variables $\{O_k\in\mathcal{O},k=1,\dots,K\}$ , for some arbitrary support set $\mathcal{O}$ , such that

(3) \begin{align}\mathbb{P}\left(\bigcap_{k=1}^K O_k^{-1}(\textbf{o}_k)\big|X\right) = \prod_{k=1}^K \mathbb{P}(O_k^{-1}(\textbf{o}_k)|X) = \prod_{k=1}^K f_{O|X_{t_k}}(\textbf{o}_k)\end{align}

for any sequence of elements $\textbf{o}_1,\dots,\textbf{o}_K$ , where $\textbf{o}_k$ denotes the time- $t_k$ network observation across all stations. Hence, observations are conditionally mutually independent given the network states. The term $f_{O|X_{t_k}}$ stands for a conditional mass function assigned to measurements across the M stations, defined with respect to a counting measure $\mu_{\mathcal{O}}$ . In this paper it is assumed that $f_{O|x}>0$ (everywhere) for all $x\in\mathcal{S}$ ; however, deterministic observations such as queue lengths can be easily approximated by means of regularized indicator functions. Extensions to continuous settings are straightforward.

Example 2. Consider a vanilla closed network as pictured in Figure 2. It includes a single FCFS service station, $K=1$ processing unit, and a delay routing a population of N jobs in a closed loop. The evolution of $X=(X_t)_{t\geq 0}$ monitors the total number of jobs within the service station, with $X_0=0$ . Variables $O_k=(O_{0,k},O_{1,k})$ at times $t_k$ approximate a system with deterministic observations of queue lengths (in both station and delay) and are supported on $\mathcal{O}=\{0,\dots,N\}^2$ . The observation model may thus factor across the network components, so that, for example, $f_{O|x}(\textbf{o}) = \tilde{f}_{O|N-x}(o_0) \cdot \tilde{f}_{O|x}(o_1)$ with $\tilde{f}_{O|x}(o)=\frac{\epsilon}{N} + \mathbb{I}(o=x) \cdot (1-\frac{N+1}{N}\epsilon)$ , and

(4) $$\matrix{ \mathbb{P}{(O_k^{ - 1}({\bf{o}})|X) = \left\{ {\matrix{ {{{(1 - {\epsilon})}^2},} \hfill & {{o_0} = N - {X_{{t_k}}},{\mkern 1mu} {o_1} = {X_{{t_k}}},} \hfill \cr {{{({\epsilon}/N)}^2},} \hfill & {{o_0} \ne N - {X_{{t_k}}},{\mkern 1mu} {o_1} \ne {X_{{t_k}}},} \hfill \cr {(1 -{\epsilon})\cdot{\epsilon}/N} \hfill & {{\rm{otherwise,}}} \hfill \cr } } \right.} \hfill \cr } $$

for a slack regularizing parameter $\epsilon>0$ and all $k=1,\dots,K$ . This will account for some $\epsilon \cdot 100 \%$ faulty measurements. Finally, note that the framework can be generalized to accommodate different observation types. For simplicity of presentation, in this paper we focus on queue lengths.

Now, let $\mathbb{P}(A|\textbf{o}_1,\dots,\textbf{o}_K)$ , $A\in\mathcal{F}$ , $(\textbf{o}_1,\dots,\textbf{o}_K)\in \mathcal{O}^K$ , denote the regular conditional probability across global events, conditioned on observations. Our primary interest lies in its induced distribution over the intensity rates (which we denote by $\mathbb{P}_{\boldsymbol{\lambda}|\textbf{o}_1,\dots,\textbf{o}_K}$ ) and trajectories. Within (Bayesian) inferential settings, these induced distributions are referred to as posteriors. For simplicity, we will restrict the problem statement to intensity rate posteriors; however, the reader will note that an analogous presentation is easily deduced for network trajectories. The posterior distribution exists and admits a density carried by its corresponding prior $\mathbb{P}_{\boldsymbol{\lambda}}$ (see Appendix A); moreover, the transformation is proportional to a weighted product of network paths, and is defined by the Radon–Nikodym derivative

(5) \begin{align}\frac{d\mathbb{P}_{\boldsymbol{\lambda}|\textbf{o}_1,\dots,\textbf{o}_K}}{d\mathbb{P}_{\boldsymbol{\lambda}}} = \frac{\int_{\mathcal{X}}\big[\prod_{k=1}^K \mathbb{P}(O_k^{-1}(\textbf{o}_k)|\textbf{t},\textbf{x})\big]\, f_{X|\boldsymbol{\lambda}}(\textbf{t},\textbf{x})\, \mu_{\mathcal{X}}(d\textbf{t},d\textbf{x})}{\mathbb{P}(O_1^{-1}(\textbf{o}_1)\cap\dots\cap O_K^{-1}(\textbf{o}_K)) } , \end{align}

which corresponds to Bayes’ equation. There, the denominator denotes a normalizing constant that integrates over trajectories and rates. This transformation will often induce a density representation $f_{\boldsymbol{\lambda}|\textbf{o}_1,\dots,\textbf{o}_K}$ for the posterior distribution with respect to a suitable (Lebesgue) base. In these cases, we may think of the above derivative as a likelihood ratio. However, this ratio poses a tractability problem; that is, the integral over trajectories cannot be computed analytically and must be approximated. This is a common problem in inferential tasks with jump processes (cf. [Reference Hobolth and Stone12, Reference Rao and Teh24, Reference Perez and Kypraios23]), and proposed solutions often rely on intensive MCMC procedures that iterate between trajectories and parameters, including direct sampling, rejection sampling, and uniformization-based methods. Yet algorithms may be hard to implement, computationally demanding, or applicable only to reduced classes of problems. In the case of QNs, strong temporal dependencies in the stochastic trajectories X impose hard coupling properties amongst rates and paths [Reference Sutton and Jordan27], which limits the applicability of state-of-the-art numerical solutions to the simplest types of network evaluation problems [Reference Perez, Hodge and Kypraios22].

Alternatively, the complex integrals in (5) may be approximated through a variational approach, where we suitably parametrize an alternative approximating measure to $\mathbb{P}$ that will decompose the integrand into multiple independent and analytically tractable parts. However, as yet there exists no approximating measure design to achieve this goal, in light of the complexity of jump processes induced by networks of queues. In the following, we present theoretical results leading to a new variational inferential design for use with jump processes, which significantly deviates from and overcomes multiple limitations found in standard methods that resort to full-independence assumptions [Reference Opper and Sanguinetti18, Reference Opper and Sanguinetti19, Reference Cohn, El-Hay, Friedman and Kupferman7, Reference Wang, Blei and Heckerman28, Reference Zechner30, Reference Zhang, Pan and Rao31].

4.1. A variational decomposition

The likelihood for observation events in (3) readily transfers to counts Y by means of (7); we thus may write $ f_{O|Y_{t_k}}(\textbf{o}_k)\equiv f_{O|X_{t_k}}(\textbf{o}_k) $ . Under the measure $\tilde{\mathbb{P}}$ , network states can have negative values; the likelihood is undefined in such instances. Now, note that piecewise $\mathcal{S}^Y$ -valued trajectories also represent elements $(\textbf{t},\textbf{y})$ in a space $(\mathcal{Y},\Sigma_{\mathcal{Y}})$ , similar to network paths in X. Let $f_{Y|\boldsymbol{\lambda},\textbf{o}_1,\dots,\textbf{o}_K}$ be a density function, with respect to some base measure $\mu_\mathcal{Y}$ , where for all $B\in\Sigma_{\mathcal{Y}}$ ,

\begin{equation*}\mathbb{P}(Y^{-1}(B)|\boldsymbol{\lambda},\textbf{o}_1,\dots,\textbf{o}_K) = \int_B f_{Y|\boldsymbol{\lambda},\textbf{o}_1,\dots,\textbf{o}_K} \,d\mu_\mathcal{Y}.\end{equation*}

It may be shown by properties of conditional distributions that, conditioned on observations, Y is a non-homogeneous semi-Markov process with hazard functions

(10) \begin{equation}\Lambda_{y,y'}(t)=\Xi_{y,y'}\cdot\frac{\mathbb{P}(\bigcap_{k : t_k>t}O_k^{-1}(\textbf{o}_{\textbf{k}})|Y_t=y')}{\mathbb{P}(\bigcap_{k : t_k>t}O_k^{-1}(\textbf{o}_{\textbf{k}})|Y_t=y)} \end{equation}

for $y'\neq y$ , and $\Lambda_{y}(t)=-\sum_{y'\neq y}\Lambda_{y,y'}(t)$ , so that

\begin{align*}f_{Y|\boldsymbol{\lambda},\textbf{o}_1,\dots,\textbf{o}_K}(\textbf{t},\textbf{y}) = \pi(y_0)\, e^{\int_{t_I}^{T} \Lambda_{y_{I}}(u)du} \, \prod_{i=1}^I \Lambda_{y_{i-1},y_{i}}(t_i)\, e^{\int_{t_{i-1}}^{t_i} \Lambda_{y_{i-1}}(u)du}\, .\end{align*}

Here, $\Xi\equiv \Xi(\delta,\boldsymbol{\lambda})$ denotes the generator matrix associated with fixed values in $\boldsymbol{\lambda}$ . For a deeper look at conditional jump processes we refer the reader to [Reference Serfozo25, Reference Daley and Vere-Jones9]. This conditional counting process is of key importance; however, the structure of rates in (10) poses a trivial analytical impediment. In our approximating effort, we assume the existence of an alternative measure $\mathbb{Q}$ on $(\Omega,\mathcal{F})$ . Under this measure, count trajectories in Y are subject to a full decomposition; that is, the $\mathbb{Q}$ -law of Y is that of a family of $|\mathcal{T}|$ independent non-homogeneous Poisson counting processes with state-dependent intensity functions $\nu^{\boldsymbol{\eta}} = (\nu^{\boldsymbol{\eta}}_t(\cdot))_{t\geq 0}$ , for all ${\boldsymbol{\eta}}\in\mathcal{T}$ . Here, the following hold:

• • Intensity rates for jumps $y\xrightarrow{t,{\boldsymbol{\eta}}} y'$ , $t\geq 0$ , are independent of $\boldsymbol{\lambda}$ , change over time, and are given by $\nu^{\boldsymbol{\eta}}_t(y_{\boldsymbol{\eta}})$ .

• • Holding rates in Y evolve according to $|\nu_t (Y_t)|$ , with $\nu_t (Y_t)=-\sum_{\boldsymbol{\eta}\in\mathcal{T}}\nu^{\boldsymbol{\eta}}_t(Y^{\boldsymbol{\eta}}_t)$ .

• • The state probability of the multivariate process Y factors across the job transition directions, so that

  \begin{equation*}\mathbb{Q}(Y_t=y) = \prod_{\boldsymbol{\eta}\in\mathcal{T}}\mathbb{Q}(Y^{\boldsymbol{\eta}}_t=y_{\boldsymbol{\eta}})\end{equation*}
  for every $y\in\mathcal{S}^Y$ .

Similar full-independence decompositions are often referred to as a variational mean-field approximations within Bayesian inferential settings [Reference Cohn, El-Hay, Friedman and Kupferman7]. In order to ensure computational tractability within forthcoming procedures, the intensity functions $\nu$ must be bounded from above by some arbitrary constant, so that $\nu^{\boldsymbol{\eta}}_t(y_{\boldsymbol{\eta}})\leq \bar{\nu}$ for all $t>0,\boldsymbol{\eta}\in\mathcal{T}$ , and $y\in\mathcal{S}^Y$ . We finally note that $\mathbb{Q}$ and $\tilde{\mathbb{P}}$ are equivalent on $\mathcal{F}$ , as both assign a positive measure to every marginally-increasing sequence of $\mathbb{N}_0^{|\mathcal{T}|}$ -valued counts.

4.2. Recapitulation and the variational lower bound

It is important to observe that the intensity functions $\nu^{\boldsymbol{\eta}}$ completely define our alternative measure. In the rest of the paper, we inspect how to fine-tune these functions to ensure that $\mathbb{Q}$ offers a good approximation to the measure $\mathbb{P}$ , conditioned on observations in $\textbf{O}$ . First, within the diagram in Figure 3 we offer an outline of the various measures, stochastic processes, probabilistic dependencies, and infinitesimal generators introduced up to this point. There, we also summarize the primary steps that help construct a variational framework to enable the inferential task. In short, these are as follows:

• • The reference measure $\mathbb{P}$ , along with the generator $Q=Q(\boldsymbol{\lambda})$ , describes the real unknown stochastic network trajectories in X.

• • $\tilde{\mathbb{P}}$ is a minor departure from $\mathbb{P}$ and expands the support of X so that queue lengths can become negative. It is associated with the generators $\tilde{Q}=\tilde{Q}(\delta,\boldsymbol{\lambda})$ and $\Xi=\Xi(\delta,\boldsymbol{\lambda})$ , for trajectories in X and job transition counts in Y, respectively.

• • $\mathbb{Q}$ is an alternative measure under which Y is a set of independent non-homogeneous counting processes. It is associated with intensity functions $\nu^{\boldsymbol{\eta}} = (\nu^{\boldsymbol{\eta}}_t(\cdot))_{t\geq 0}$ for each count process in the directions $\boldsymbol{\eta}\in\mathcal{T}$ .

Figure 3: Overview of the approximating variational framework presented in this paper, summarizing the various measures used and the primary goals that each step will accomplish.

In our next result, a variational lower bound will interrelate the transition rates $\boldsymbol{\lambda}$ for jumps under $\mathbb{P}$ and $\tilde{\mathbb{P}}$ with the intensity functions $\nu^{\boldsymbol{\eta}}$ , in light of the available observations/measurements in $\textbf{O}$ . From now on, transition rates $\boldsymbol{\lambda}$ are assumed mutually independent under $\mathbb{Q}$ and admit undetermined densities $d\mathbb{Q}_{\boldsymbol{\lambda}} = \{d\mathbb{Q}_{\lambda,\boldsymbol{\eta}}, \ \boldsymbol{\eta}\in\mathcal{T}\}$ that must integrate to 1 on $(\mathbb{R}_+,\mathcal{B}(\mathbb{R}_+))$ .

Lemma 1. (Variational lower bound.) Define $\textbf{O}=\bigcap_{k=1}^K O_k^{-1}(\textbf{o}_k)$ , and let $\tilde{\mathbb{P}}$ and $\mathbb{Q}$ be the probability measures on $(\Omega,\mathcal{F})$ as defined above. Recall the notation $\textstyle\Xi_{y,y'}\equiv\textstyle\Xi_{y,{\boldsymbol{\eta}}}$ for jumps $y\xrightarrow{{\boldsymbol{\eta}}} y'$ with direction ${\boldsymbol{\eta}}$ . Then

(11) \begin{align}\log\tilde{\mathbb{P}}(\textbf{O}) &\geq \sum_{k=1}^K\mathbb{E}^{\mathbb{Q}}_{Y_{t_k}} \big[ \log f_{O|Y_{t_k}}(\textbf{o}_k) \big] -\mathbb{E}^{\mathbb{Q}}_{\boldsymbol{\lambda}} \bigg[ \log\frac{ d\mathbb{Q}_{\boldsymbol{\lambda}} }{d\mathbb{P}_{\boldsymbol{\lambda}}} \bigg] \nonumber \\[3pt]&\quad - \int_0^T \mathbb{E}^{\mathbb{Q}}_{Y_t,\boldsymbol{\lambda}} \bigg[ \sum_{{\boldsymbol{\eta}}\in\mathcal{T}} \nu^{\boldsymbol{\eta}}_t (Y^{\boldsymbol{\eta}}_t) \log \frac{ \nu^{\boldsymbol{\eta}}_t(Y^{\boldsymbol{\eta}}_t)}{\Xi_{Y_t,{\boldsymbol{\eta}}}} - \Xi_{Y_t} + \nu_t(Y_t) \bigg] dt \end{align}

offers a lower bound on the $\tilde{\mathbb{P}}$ -probability of retrieved observation events.

Proof. Note that

(12) \begin{align}\log\tilde{\mathbb{P}}(\textbf{O})&=\log\int_{\mathcal{Y}\times\mathbb{R}_+^{|\mathcal{T}|}} \mathbb{P}(\textbf{O}|Y)\, d(Y,\boldsymbol{\lambda})_*\tilde{\mathbb{P}} = \log \mathbb{E}^{\mathbb{Q}}_{Y,\boldsymbol{\lambda}} \bigg[ \mathbb{P}(\textbf{O}|Y)\, \frac{d(Y,\boldsymbol{\lambda})_*\tilde{\mathbb{P}}}{d(Y,\boldsymbol{\lambda})_*\mathbb{Q}} \bigg] \nonumber \\[3pt] &\geq \mathbb{E}^{\mathbb{Q}}_{Y} \big[ \log \mathbb{P}(\textbf{O}|Y)\big] - \mathbb{E}^{\mathbb{Q}}_{Y,\boldsymbol{\lambda}} \bigg[ \log \frac{d(Y,\boldsymbol{\lambda})_*\mathbb{Q}}{d(Y,\boldsymbol{\lambda})_*\tilde{\mathbb{P}}} \bigg], \end{align}

where we use Jensen’s inequality for finite measures. This is known as a variational lower bound on the log-likelihood, and

\begin{align*}\mathbb{E}^{\mathbb{Q}}_{Y} \big[ \log \mathbb{P}(\textbf{O}|Y)\big] &= \mathbb{E}^{\mathbb{Q}}_{Y} \big[ \log \prod_{k=1}^K f_{O|Y_{t_k}}(\textbf{o}_k) \big] = \sum_{k=1}^K\mathbb{E}^{\mathbb{Q}}_{Y_{t_k}} \big[ \log f_{O|Y_{t_k}}(\textbf{o}_k) \big]\end{align*}

follows directly from (3). The negative part in (12) is the Kullback–Leibler (KL) divergence between image measures of $\mathbb{Q}$ and $\tilde{\mathbb{P}}$ . By noting that these share base measures, and $Y,\boldsymbol{\lambda}$ are independent under $\mathbb{Q}$ , we have

(13) \begin{align}\mathbb{E}^{\mathbb{Q}}_{Y,\boldsymbol{\lambda}} \bigg[ \log \frac{d(Y,\boldsymbol{\lambda})_*\mathbb{Q}}{d(Y,\boldsymbol{\lambda})_*\tilde{\mathbb{P}}} \bigg] &= \mathbb{E}^{\mathbb{Q}}_{\boldsymbol{\lambda}} \bigg[ \log\frac{ d\mathbb{Q}_{\boldsymbol{\lambda}} }{d\mathbb{P}_{\boldsymbol{\lambda}}} \bigg] \, + \, \mathbb{E}^{\mathbb{Q}}_{\boldsymbol{\lambda}} \bigg[ \mathbb{E}^{\mathbb{Q}}_{Y} \bigg[ \log \frac{g_Y}{f_{Y|\boldsymbol{\lambda}}} \bigg] \bigg] , \end{align}

where $g_Y$ and $f_{Y|\boldsymbol{\lambda}}$ denote the Y-trajectory densities associated with rates $\nu^{\boldsymbol{\eta}}$ and $\Xi(\boldsymbol{\lambda})$ , respectively. The last term in (13) is a $\mathbb{Q}$ -average of the KL divergence on Y, where the mean is taken across the infinitesimal transition rates. For a fixed starting $Y_0\in\mathcal{S}^Y$ , the inner expectation is shown in [Reference Opper and Sanguinetti18] to take the equivalent form

\begin{align*}\mathbb{E}^{\mathbb{Q}}_{Y} \bigg[ \log \frac{g_Y}{f_{Y|\boldsymbol{\lambda}}} \bigg] &= \lim_{R\rightarrow\infty} \sum_{r=0}^{R-1} \mathbb{E}^{\mathbb{Q}}_{Y_{\frac{Tr}{R}}} \Bigg[ \sum_{y\in\mathcal{S}^Y} \mathbb{Q}(Y_{\frac{T(r+1)}{R}}=y|Y_{\frac{Tr}{R}}) \log\frac{\mathbb{Q}(Y_{\frac{T(r+1)}{R}}=y|Y_{\frac{Tr}{R}})}{\tilde{\mathbb{P}}(Y_{\frac{T(r+1)}{R}}=y|Y_{\frac{Tr}{R}},\boldsymbol{\lambda})} \Bigg] .\end{align*}

Note that within an infinitesimal time interval a jump in Y may only happen in one direction within $\mathcal{T}$ . With this in mind, we retrieve the limit of a Riemann sum in the interval [0, T], i.e.

\begin{align*}\mathbb{E}^{\mathbb{Q}}_{Y} \bigg[ \log \frac{g_Y}{f_{Y|\boldsymbol{\lambda}}} \bigg] = &\lim_{R\rightarrow\infty} \frac{T}{R }\sum_{r=0}^{R-1} \mathbb{E}^{\mathbb{Q}}_{Y_{\frac{Tr}{R}}} \Bigg[ \sum_{{\boldsymbol{\eta}}\in\mathcal{T}} \nu^{\boldsymbol{\eta}}_{\frac{Tr}{R}}(Y^{\boldsymbol{\eta}}_{\frac{Tr}{R}}) \log \frac{ \nu^{\boldsymbol{\eta}}_{\frac{Tr}{R}}(Y^{\boldsymbol{\eta}}_{\frac{Tr}{R}})}{\Xi_{Y_{\frac{Tr}{R}},{\boldsymbol{\eta}}}} + \frac{R}{T} \log\frac{1 + \frac{T}{R} \nu_{\frac{Tr}{R}}(Y_{\frac{Tr}{R}})}{1+ \frac{T}{R} \Xi_{Y_{\frac{Tr}{R}}}} \Bigg] \\[3pt] = &\int_0^T \mathbb{E}^{\mathbb{Q}}_{Y_t} \bigg[ \sum_{{\boldsymbol{\eta}}\in\mathcal{T}} \nu^{\boldsymbol{\eta}}_t (Y^{\boldsymbol{\eta}}_t) \log \frac{ \nu^{\boldsymbol{\eta}}_t(Y^{\boldsymbol{\eta}}_t)}{\Xi_{Y_t,{\boldsymbol{\eta}}}} - \Xi_{Y_t} + \nu_t(Y_t) \bigg] dt ,\end{align*}

and

\begin{align*}\mathbb{E}^{\mathbb{Q}}_{Y,\boldsymbol{\lambda}} \bigg[ \log \frac{d(Y,\boldsymbol{\lambda})_*\mathbb{Q}}{d(Y,\boldsymbol{\lambda})_*\tilde{\mathbb{P}}} \bigg] = \ &\mathbb{E}^{\mathbb{Q}}_{\boldsymbol{\lambda}} \bigg[ \log\frac{ d\mathbb{Q}_{\boldsymbol{\lambda}} }{d\mathbb{P}_{\boldsymbol{\lambda}}} \bigg]\\[3pt] &+ \int_0^T \mathbb{E}^{\mathbb{Q}}_{Y_t,\boldsymbol{\lambda}} \bigg[ \sum_{{\boldsymbol{\eta}}\in\mathcal{T}} \nu^{\boldsymbol{\eta}}_t (Y^{\boldsymbol{\eta}}_t) \log \frac{ \nu^{\boldsymbol{\eta}}_t(Y^{\boldsymbol{\eta}}_t)}{\tilde{\Xi}_{Y_t,{\boldsymbol{\eta}}}} - \tilde{\Xi}_{Y_t} + \nu_t(Y_t) \bigg] dt\end{align*}

completes the proof.

Thus, the lower bound in (11) depends on both the latent variables Y and $\boldsymbol{\lambda}$ , accounting for the various counts and rates. On a basic level, this is built from three distinct components: the expected log-observations, the KL divergence across service rate densities, and a $\mathbb{Q}$ -weighted divergence across hazard functions and rates, further integrated along the entire network trajectory.

6.1. Single-class closed network

We begin by analysing the small closed network example previously shown in Figure 2. We recall that this includes one FCFS service station, with $K_1=1$ processing unit, along with a delay, together processing a population of N jobs cyclically in a closed loop. All jobs belong to the same class and have equal service rates; we denote by $\mu_1$ the job processing rate within the service station. On completion, a job proceeds to the delay station where it awaits for an exponentially distributed time before being routed back to the queue. We use $\mu_0$ to denote the delay rate, and we note that the arrival rate to the queue is directly proportional to the number of jobs at the delay.

Both stations are independent and $\mu_0$ is fixed in order to ensure model identifiability within the service station. In this instance, the network topology is deterministic and trivial, and the evolution of $X=(X_t)_{t\geq 0}$ monitors the total number of jobs within the service station, with $X_0=0$ . The generator Q of X is finite and satisfies

\begin{equation*}Q= \begin{array}{c} \\ \\[-7pt] 0 \\ \\[-7pt] 1\\ \\[-7pt] \vdots\\ \\[-7pt] N-1 \\ \\[-7pt] N \end{array} \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} & 0 & 1 & 2 & \dots & N-2 & N-1 & N \\ \\[-7pt] & -N\mu_0 & N\mu_0 & 0 & \dots & 0 & 0 & 0\\ \\[-7pt] & \mu_1 & -N\mu_0+\mu_0-\mu_1 & (N-1)\mu_0 & \dots & 0 & 0 & 0 \\ \\[-7pt] & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ \\[-7pt] & 0 & 0 & 0 & \dots & \mu_1 & -(\mu_0+\mu_1) & \mu_0 \\ \\[-7pt] & 0 & 0 & 0 & \dots & 0 & \mu_1 & -\mu_1 \end{array}\right],\end{equation*}

where row and column labels denote the number of jobs in the queueing station. Since $\mu_0$ is fixed, our interest lies in $\lambda\equiv \lambda_{1,0} = \mu_1 \cdot p_{1,0} = \mu_1$ . We assign to this rate a distribution $\mathbb{P}_\lambda\equiv \lambda_\star\mathbb{P}$ with (gamma) density $f_\lambda$ such that its hyperparameters fix some reasonably uninformative prior knowledge on the system. We monitor the delay and FCFS service station at fixed and equally spaced times $t_1<\dots<t_K$ in an interval [0, T]. Here, variables $O_k$ are supported on $\mathcal{O}=\{0,\dots,N\}^2$ , and the observation model factors across the network components so that $f_{O|x}(\textbf{o}) = \tilde{f}_{O|N-x}(o_0) \cdot \tilde{f}_{O|x}(o_1)$ with $\tilde{f}_{O|x}(o)=\frac{\epsilon}{N} + \mathbb{I}(o=x) \cdot (1-\frac{N+1}{N}\epsilon)$ and

(18) $$\matrix{ \mathbb{P}{(O_k^{ - 1}({\bf{o}})|X) = \left\{ {\matrix{ {{{(1 -{\epsilon})}^2},} \hfill & {{o_0} = N - {X_{{t_k}}},{\mkern 1mu} {o_1} = {X_{{t_k}}},} \hfill \cr {{{({\epsilon}/N)}^2},} \hfill & {{o_0} \ne N - {X_{{t_k}}},{\mkern 1mu} {o_1} \ne {X_{{t_k}}},} \hfill \cr {(1 -{\epsilon})\cdot{\epsilon}/N} \hfill & {{\rm{otherwise}}} \hfill \cr } } \right.} \hfill \cr } $$

for $\epsilon>0$ and all $k=1,\dots,K$ . This accounts for some $\epsilon \cdot 100 \%$ faulty measurements; also, we note that a system with discrete observations is approximated as $\epsilon\rightarrow 0$ . Now, assume there exist some sample observations $\textbf{o}_1,\dots,\textbf{o}_K$ from a model realization in this closed network. These can easily be produced from (18) given a trajectory $(X_t)_{t\in[0,T]}$ . In order to produce the trajectory from (2), given the service rates, we may employ Gillespie’s algorithm [Reference Gillespie10] or faster uniformization-based alternatives [Reference Rao and Teh24].

In the following, we analyse simulated data ( $N=50,\,\mu_0=0.1,\,\epsilon=0.2,\,T=100,\,K=50$ ) by assigning a conjugate gamma density to $\lambda\equiv\mu_1$ , so that $\lambda\sim\Gamma(\alpha,\beta)$ with $\alpha=5$ and $\beta=2$ under the reference measure $\mathbb{P}$ . Recall that $\mu_0$ is fixed to ensure model identifiability, and $X_t\in\{0,\dots,N\}$ denotes the number of jobs in the service station at any time $t\geq0$ . For later reference, the stationary distribution of the system is given by

\begin{align*}\pi_{\mathbb{P}}(x|\lambda) = \lim_{t\rightarrow\infty} \mathbb{P}(X_t=x|\lambda) = \left(\frac{\mu_0}{\lambda}\right)^x\frac{1}{(N-x)!}\bigg/\sum_{x=0}^N\left(\frac{\mu_0}{\lambda}\right)^x\frac{1}{(N-x)!},\end{align*}

so that, assuming the observations are sufficiently spaced and the system has reached stationarity, it holds that

(19) \begin{align}\mathbb{P}\Big(\bigcap_{k=1}^K O_k^{-1}(\textbf{o}_k)\big|\lambda\Big) \approx \prod_{k=1}^K \sum_{x=0}^N f_{O|x}(\textbf{o}_k) \pi_{\mathbb{P}}(x|\lambda) = \frac{\prod_{k=1}^K \sum_{x=0}^N f_{O|x}(\textbf{o}_k) \left(\frac{\mu_0}{\lambda}\right)^x\frac{1}{(N-x)!}}{\Big(\sum_{x=0}^N\left(\frac{\mu_0}{\lambda}\right)^x\frac{1}{(N-x)!}\Big)^K}, \end{align}

where $f_{O|x}(\textbf{o}_k)$ is as defined in (18). Note that here $|\mathcal{T}|=2$ , and the process $Y=(Y^{0,1}_t,Y^{1,0}_t)_{t\geq 0}$ monitors transitions between the delay and service station, in both the directions $0\rightarrow 1$ and $1\rightarrow 0$ . The lower bound to the log-likelihood in (11) reduces to

(20) \begin{align}\log\tilde{\mathbb{P}}(\textbf{O}) \geq \ & \sum_{k=1}^K\mathbb{E}^{\mathbb{Q}}_{Y_{t_k}} \big[ \log \mathbb{P}(O^{-1}_k(\textbf{o}_k) |Y^{0,1}_{t_k} - Y^{1,0}_{t_k}) \big]\nonumber\\[3pt] &-\mathbb{E}^{\mathbb{Q}}_{\lambda} \bigg[ \log\frac{ g_{\lambda}^{1,0} }{ d\mathbb{P}_{\lambda} } \bigg] + \int_0^T \mathbb{E}^{\mathbb{Q}}_{Y_t,\lambda} \big[ \Psi[Y_t,\nu_t,\lambda] \big] dt \end{align}

with

\begin{align*}\Psi[Y_t,\nu_t,\lambda] = \ & \nu^{1,0}_t (Y^{1,0}_t) + \nu^{0,1}_t (Y^{0,1}_t) - 2\delta - \lambda\cdot\mathbb{I}(Y^{0,1}_{t_k} - Y^{1,0}_{t_k} > 0) \\[3pt] &- \mu_0\cdot (N + Y^{1,0}_{t_k} - Y^{0,1}_{t_k}) - \nu^{1,0}_t (Y^{1,0}_t) \log \frac{ \nu^{1,0}_t (Y^{1,0}_t)}{\delta + \lambda\cdot\mathbb{I}(Y^{0,1}_{t_k} - Y^{1,0}_{t_k} > 0)} \\[3pt] &- \nu^{0,1}_t (Y^{0,1}_t) \log \frac{ \nu^{0,1}_t (Y^{0,1}_t) }{\delta + \mu_0\cdot (N + Y^{1,0}_{t_k} - Y^{0,1}_{t_k}) },\end{align*}

so that it contains only two hazard functions in the approximating measure $\mathbb{Q}$ , namely $\nu^{0,1}$ and $\nu^{1,0}$ . In (20), we again notice that the lower bound is dominated by three distinguishable components: (i) the expected log-observations, (ii) the KL divergence across the service rate density, and (iii) a weighted $\mathbb{P}$ -to- $\mathbb{Q}$ divergence in the expected path likelihood, further integrated along the entire network trajectory. The differential equations for functionals in (16) reduce to

\begin{align*}\frac{dr_t^{0,1}(y)}{dt} &= r_t^{0,1}(y) \big(\delta + \mu_0 \cdot \mathbb{E}^{\mathbb{Q}}_{Y^{1,0}_t} \big[ (Y^{1,0}_t - y)\vee 0) \big] \big) \\[3pt]&- \frac{1+ k_t^{0,1}(y)/\mathbb{Q}(Y_t^{0,1}=y)}{e^{ k_t^{0,1}(y)/\mathbb{Q}(Y_t^{0,1}=y)}} r_t^{0,1}(y+1) e^{\mathbb{E}^{\mathbb{Q}}_{Y^{1,0}_t} \big[ \log(\delta + \mu_0\cdot[(Y^{1,0}_t-y)\vee 0]) \big]}\end{align*}

and

\begin{align*}\frac{dr_t^{1,0}(y)}{dt} &= r_t^{1,0}(y) \big(\delta + \mathbb{E}^{\mathbb{Q}}_{\lambda}[\lambda]\cdot \mathbb{Q}(Y^{0,1}_t>y) \big) \\&\quad - \frac{1+k_t^{1,0}(y)/\mathbb{Q}(Y_t^{1,0}=y)}{e^{ k_t^{1,0}(y)/\mathbb{Q}(Y_t^{1,0}=y)}} r_t^{1,0}(y+1) e^{\mathbb{E}^{\mathbb{Q}}_{Y^{0,1}_t,\boldsymbol{\lambda}} \big[ \log(\delta + \lambda \cdot \mathbb{I}(Y^{0,1}_t>y)) \big]} .\end{align*}

In Figure 4 (left) we observe the evolution of the lower bound (20) during the iterative inferential procedure, for a sufficiently small and negligible value of $\delta$ . There, we notice that the procedure has converged to a (local) optimum within approximately 13 iterations. On the right-hand side of the figure, we further observe credible intervals and point estimates across iterations for $\lambda$ , under the approximating measure $\mathbb{Q}$ . Displayed in the back (in grey) are the estimates obtained using the proposed method in this paper; in the front (in blue), we find approximations extracted by adapting variational procedures in [Reference Opper and Sanguinetti18], to allow for prior knowledge and conjugacy properties.

Figure 4: Left, evolution of lower bound to the log-likelihood during the inferential procedure. Right, 95% credible intervals and point estimates for the service rate $\lambda$ under $\mathbb{Q}$ ; the back (grey) corresponds to the proposed method, while the front (blue) is for an (adapted) traditional variational procedure.

Next, in Figure 5 (left) we find the $\mathbb{P}$ -prior density for $\lambda$ (flat density in the back of the figure) with the final-iteration $\mathbb{Q}$ -posterior superimposed (grey density in the front); the additional red and blue densities represent

• • the posterior density obtained through Metropolis–Hastings MCMC, by means of strong stationarity assumptions leading to the likelihood function shown in (19), and

• • the approximate posterior obtained by adapting variational procedures in [Reference Opper and Sanguinetti18].

Figure 5: Left, prior density (light grey flat density in the back) and posterior density (dark grey density in the front) for $\lambda$ , along with MCMC (red) and traditional variational (blue) density estimates. The black dot on the horizontal axis represents the original value in the network simulation. Right, network observations along with mean-average network trajectory and 95% credible interval for job counts in the service station; in grey, our proposed method, in blue, existing variational alernative method.

On the right-hand side we have the network observations on both the service station and the delay. Delay observations are displayed by subtracting their value from the job population N (which represents a second measurement on the service station). Whenever both observations match, these are displayed with a large dot. Along with it, we find the following:

• • in grey, a mean-average network trajectory and 95% credible interval for job counts on the service station $X_t=(Y^{0,1}_t-Y^{1,0}_t)_{t\geq 0}$ , under the approximating measure $\mathbb{Q}$ and with methods introduced in this paper;

• • in blue, a similar confidence interval and mean-average path obtained using benchmark methods in [Reference Opper and Sanguinetti18].

Noticeably, the average trajectory under our proposed variational technique flows through the most informative observations (thick dots), and the credible interval widens up to account for some faulty measurements within either network station. On the other hand, traditional variational approaches quickly converge to a local optimum and restrict mean-average dynamics, further compressing confidence intervals in regions with noisy data. In the next example, we notice how this poses a problem for traditional methods; that is, within complex and synchronized stochastic processes, we will fail to obtain sensible estimates for network parameters. Moreover, in our next example, inference by MCMC/forward–backward methods is virtually intractable (cf. [Reference Perez, Hodge and Kypraios22]).

6.2. Multi-class parallel tandems with bottleneck and service priorities

We analyse an open multi-class QN as pictured in Figure 6. In this network, there exists two classes ( $c=1,2$ ) of jobs that simultaneously transit the system. The first class consists of high-priority jobs with low arrival and service intensity rates. The second class consists of low-priority jobs with high arrival and service rates. Once a job enters the system, a probabilistic routing junction (pictured as a square within the figure) sends this job through either a PS or a priority-FCFS tandem; later, it will be serviced within an infinite station before leaving the network. In the top processor-sharing tandem, each station has five processing units. These will fraction their working capacity as seen in (8), in order to simultaneously service all jobs regardless of their class and priority level; however, service rates will differ depending on the job class. By contrast, the bottom tandem includes two FCFS stations with a single processing unit and priority scheduling. Here, low-priority jobs are serviced only if each station is entirely empty of any high-priority jobs. Consequently, the station loads in (9) are rewritten so that

\begin{align*}\textstyle\textstyle\Upsilon(y,i,1) = 1 \wedge x_{i,1} \quad \text{and} \quad \textstyle\textstyle\Upsilon(y,i,2) = (1 \wedge x_{i,2}) \cdot \mathbb{I}(x_{i,1} <1)\end{align*}

at stations $i\in\{2,4\}$ and for any $y\in\mathcal{S}^Y$ , where we recall that

\begin{equation*}x_{i,c} = \sum_{{\boldsymbol{\eta}}\in\mathcal{T}^{\leftarrow}_{i,c}} y_{{\boldsymbol{\eta}}} - \sum_{{\boldsymbol{\eta}}\in\mathcal{T}^{\rightarrow}_{i,c}} y_{{\boldsymbol{\eta}}} \end{equation*}

and thus $x_{2,c}=y_{0,2,c}-y_{2,4,c}, x_{4,c}=y_{2,4,c}-y_{4,5,c}$ , for $c=1,2$ . Due to the presence of service priorities, the ordering of jobs within the queue is irrelevant (this is also the case with random order disciplines); hence, our inferential framework allows for different service rates assigned to jobs in each class. Finally, the last service station includes an infinite number of processing units, and processing rates also differ depending on the job class.

Figure 6: Open QN with one routing junction (pictured as a square) and five service stations with varied disciplines and processing rates.

We analyse synthetic data created during a time interval [0, T] ( $T=100$ ), with arrival intensities $\mu_0^1=0.5$ , $\mu_0^2=3$ , routing probabilities $p^c_{0,i}=0.5$ , $i,c\in\{1,2\}$ , and service rates as shown in Table 1. We collect a reduced set of noiseless and equally spaced observations with $K=50$ ; these are essentially snapshots of the full system state across its service stations and job classes, so that $\mathcal{O}=\mathbb{N}^{10}$ , and the observation density in (3) is defined with

\begin{equation*}f_{O|x}(\textbf{o}) = \prod_{i=1}^5 \prod_{c=1}^2 \mathbb{I}(x_{i,c} = o_{i,c})\end{equation*}

for $x\in\mathcal{S}$ , where $o_{i,c}$ is an indexed observation in the element $\textbf{o}$ denoting the class-c queue length at station $i>0$ . (Source code for the data simulation process may be found at https://github.com/IkerPerez/variationalQueues.) Within the inferential procedure, this observation likelihood must be approximated with some regularized variant similar to (18), while taking $\epsilon\rightarrow 0$ . Next, we assign conjugate gamma priors to the various service intensities; in order to ensure identifiability in the problem, arrival rates and routing probabilities are fixed, and we focus this inferential task on the various service stations. Hence $\boldsymbol{\lambda}\equiv \{\mu^c_i\, :\, c=1,2 \; \text{and} \; i=1,\cdots,5\}$ , and we set $\lambda_{\boldsymbol{\eta}}\sim\Gamma(1,0.3)$ under the reference measure $\mathbb{P}$ , for all $\boldsymbol{\eta}\in\mathcal{T}$ .

Table 1: Summary statistics for posterior service rates in the QN in Figure 6.

In the following, we omit the cumbersome mathematical details related to this complex model formulation, and we focus on discussing prior choices, calibration of the algorithm, results, and method comparisons following the inferential procedure.

Remarks on using MCMC data-augmentation for inference. Transient inference in a stochastic system with priorities is especially challenging, because of the strong dependencies this generates on the queue lengths across the stations and classes. Specifically,

• • data-augmentation methods relying on MCMC techniques do not scale (cf. [Reference Sutton and Jordan27; Reference Perez, Hodge and Kypraios22]), as dependences yield very autocorrelated output chains;

• • there exist no analytic product-form distributions to enable approximate inferential methods under assumptions of system stationarity, as discussed in the previous example; and

• • generator-matrix exponentiations with a forward–backward algorithm [Reference Zhao32] are simply unscalable to such large multivariate systems.

Remarks on using benchmark variational methods for inference. Note that traditional variational methods (cf. [Reference Opper and Sanguinetti18, Reference Cohn, El-Hay, Friedman and Kupferman7]) are centred around populations or lengths in the individual queues. In this example, populations may not be factorized under an approximating measure $\mathbb{Q}$ , since system jumps are synchronized; i.e., a jump down in one queue corresponds to a jump up in another. As a consequence, pairs of approximating rates under $\mathbb{Q}$ will be interlinked with the same real transition rate under $\mathbb{P}$ , and derivations such as the lower bound in (11) or Equations (14)–(16) are unattainable. For the sake of completeness and comparisons, we adapt existing variational algorithms to the current task; however, we must make the following changes:

• • We allow a factorization $\mathbb{Q}(X_t=\textbf{x})=\prod_{i}\mathbb{Q}(X_{i,t}=x_i)$ , so that jobs may be virtually created and removed in any queue; i.e., jobs do not transit a network, but rather reach and depart servers individually. The full population of jobs in the network is not preserved.

• • We duplicate intensities for transitions in the real model; that is, we have a rate for (i) a job departing a queue and (ii) the job arriving at another. Technically, a job could arrive at a new server before it departs the previous one; synchronization is lost; point estimates for parameters are averaged across pairs and weighted for network load.

As we will see, this leads to drastic performance issues that render the method unusable.

Figure 7: 95% credible intervals for queue lengths across the service stations. Dark (light) grey corresponds to high-priority (low-priority) jobs. Bottom right panel shows expected jump intensity and station load in the direction $\eta=(0,1,1)$ .

6.2.2 Results

In the plots in Figure 7, with the exception of the bottom right one, we observe 95% credible intervals for the queue length processes $X^{i,c}_t$ over time, across the various service stations and job classes. In the figure, intervals in dark grey relate to high-priority jobs, and their corresponding queue length observations are represented by black circles. This information is superimposed over its analogue for low-priority jobs, where intervals are coloured in light grey and observations represented by small diamonds. These interval approximations ignore small positive densities that are sometimes assigned to negative queue lengths. Note that this is a consequence of employing counts across job transitions in Y as a basis for inference on X; however, we recall that this is necessary to overcome the coupling challenges described in Sections 2 and 4. Overall, we note that the variational flow captures the collected observations well, with some few exceptions in the stations with priority scheduling; hence, it offers a good basis for building approximate estimates for parameters and the likelihood ratio (5).

The bottom right plot in Figure 7 shows an overview of the expected jump intensity

\begin{equation*}\mathbb{E}^{\mathbb{Q}}_{Y^{\boldsymbol{\eta}}_t} [ \nu^{\boldsymbol{\eta}}_t (Y^{\boldsymbol{\eta}}_t) ]\end{equation*}

and station load

\begin{equation*}\mathbb{E}^{\mathbb{Q}}_{Y_t} [ \Upsilon(Y_t,\eta_1,\eta_3)\vee 0 ]\end{equation*}

in the direction $\boldsymbol{\eta}=(0,1,1)$ at times $t\in[0,T]$ . The sharp peaks in the intensities come at observation times, ensuring that the process density transits through the observations. Finally, we notice that the expected station load differs from the estimate of the high-priority queue length in station 1, as this process combines and weights the queue length across the two priorities according to (8).

Next, we find in Table 1 summary statistics for the posterior service rates under the approximating variational measure $\mathbb{Q}$ , along with point estimates obtained by adapting benchmark variational techniques in [Reference Opper and Sanguinetti18]. We observe that the proposed framework allows us to gain a good overview of the system properties and variability in the processing speed across the various stations; by contrast, existing methods are far from offering reasonable approximations to system behaviour (they instead seem to construct an averaged estimate of network flow). Noticeably, there exist a few significant deviations from real values, within the posterior estimates for high-priority service rates in PS stations. This is likely due to a combination of sampling variance, high model complexity, and the limitations of such approximate variational procedures for transient analyses of stochastic processes.