# Utilizing Radiomic Feature Analysis for Automated MRI Keypoint Detection:
Enhancing Graph Applications

###### Abstract

Graph neural networks (GNNs) present a promising alternative to CNNs and
transformers in certain image processing applications due to their parameter-
efficiency in modeling spatial relationships. Currently, a major area of
research involves the converting non-graph input data for GNN-based models,
notably in scenarios where the data originates from images. One approach
involves converting images into nodes by identifying significant keypoints
within them. SuperRetina, a semi-supervised technique, has been utilized for
detecting keypoints in retinal images. However, its limitations lie in the
dependency on a small initial set of ground truth keypoints, which is
progressively expanded to detect more keypoints. Having encountered
difficulties in detecting consistent initial keypoints in brain images using
SIFT and LoFTR, we proposed a new approach: radiomic feature-based keypoint
detection. Demonstrating the anatomical significance of the detected keypoints
was achieved by showcasing their efficacy in improving registration processes
guided by these keypoints. Subsequently, these keypoints were employed as the
ground truth for the keypoint detection method (LK-SuperRetina). Furthermore,
the study showcases the application of GNNs in image matching, highlighting
their superior performance in terms of both the number of good matches and
confidence scores. This research sets the stage for expanding GNN applications
into various other applications, including but not limited to image
classification, segmentation, and registration.

Index Terms—  Image Matching, Image Registration, Kepoint Detection, Radiomic
Features, Brain MRI, GNN

##  1 Introduction

Graph neural networks (GNNs) have shown promising results for reducing the
computational requirements for certain image processing tasks as shown in [1].
However, converting images into graphs is an active area of research. Many
papers have tried breaking down images into patches and treating each patch as
a node in a graph [1]. Our new method relies on detecting important keypoints
in the images along with their features and making graphs from them.

Detecting important keypoints in certain types of images, such as magnetic
resonance images of the brain is not straightforward due to the lack of well-
defined landmarks. In various registration competitions, such as the BraTSReg
challenge [2], experts have had to mark specific landmark points on 3D brain
images to help the participants align images and evaluate how well their
registration algorithms worked. But this marking is neither easy nor fast.
Also, only a few points, around 6 to 50 per set of images, get marked.

Traditional keypoint detection algorithms such as SIFT [3] fall behind deep
learning-based keypoint detection algorithms which have different types
include: supervised, unsupervised, and semi-supervised methods. Some examples
are UnsuperPoint [4], SuperPoint [5], GLAMpoints [6], and SuperRetina [7].

After trying and failing to find reliable keypoints in brain images using
methods such as SIFT [3] and LoFTR [8], we proposed a new algorithm. Our
proposed method finds keypoints using radiomic features in brain images and
their feature vectors. We proved these keypoints are important landmarks in
brain images by showing how they help in registering brain images. Then we
used this dataset we made, which consists of MRI images from OASIS dataset [9]
and the keypoints we detected, to train the LK-SuperRetina algorithm [10] to
detect new keypoints. We showed that using a GNN like SuperGlue [11] for image
matching improves the features of detected keypoints. This method of detecting
keypoints opens doors for using graph-based neural networks for different
tasks like brain classification, segmentation, and registration.

##  2 Proposed Method

We introduced an approach to keypoint detection based on radiomic features.
Given an image alongside its segmentation labels, our method identifies
radiomic keypoints as the centers of radiomic segmentation labels within the
image. These radiomic features encompass a range of intensity and shape
characteristics specific to the regions defined by the segmentation map. To
compute these features, we utilized the Pyradiomic library [12].

Radiomic keypoints are closely tied to segmentation regions predicted by
neural networks. We trained Swin UNetR [13] to predict the segmentation maps
of brain images and used them as masks to extract the keypoints. Each
segmentation mask yields a keypoint location, accompanied by 53 descriptive
radiomic features, see figure 1 for more details. Radiomic keypoints exhibit
repeatability across various brain samples even in the presence of varying
intensity and non-rigid deformations which make our proposed method a more
robust alternative to methods like SIFT [3]. These detected keypoints can
serve as a initial keypoints for detecting additional keypoints through deep
learning techniques as we will show in our results, in which we trained LK-
SuperRetina [10] to automatically detect these ground truth keypoints and
extra keypoints. This underscores the significance of these keypoints as
essential landmarks within brain images.

Fig. 1: (a) The image matching pipeline, encompassing keypoint detection using
neural networks, graph formation from detected keypoints, graph neural network
(GNN) processing to enhance keypoint features, and a dedicated head for
keypoint matching. (b) Radiomic features-based Keypoint detection method.

##  3 Experiments and Results

This section begins with an introduction to the dataset used in all our
experiments. Following that, we explore various approaches applied for
keypoint detection, explaining how the detected keypoints were utilized to
improve registration and identify additional keypoints. Finally, we will
conclude with our results on GNN-based image matching.

###  3.1 OASIS dataset

The OASIS dataset, as described in [9], contains MRI data obtained from 414
subjects. This dataset has been divided into three separate subsets for
training, validation, and testing, following a ratio of 314:50:50,
respectively. Each subject in the dataset has T1-weighted scan, as well as
segmentation masks of various regions of the brain. This dataset incorporates
three distinct types of brain segmentation: a four-label mask, a thirty-five
label mask, and a twenty-four label mask.The 3D T1-weighted scans and their
corresponding masks are of resolution (160, 192, 224). Additionally, the 2D
T1-weighted scans and their corresponding masks have a resolution of (160,
192).

###  3.2 SIFT

SIFT [3] was initially chosen for experimentation due to its reputation in
keypoint detection. Utilizing gradients, SIFT excels in identifying scale and
rotation invariant keypoints. The method calculates orientation gradients
across scales, creating a distinctive 128-dimensional vector for each
keypoint.

Figure 2 shows significant issues we encountered with keypoint detection using
SIFT. Notably, there was inconsistency in keypoint locations across different
MRI slices within the brain, reflecting a lack of repeatability. Our
experiments highlighted SIFT’s ineffectiveness under conditions involving
large deformations.

###  3.3 LoFTR

Due to SIFT’s inability to establish correspondences, we opted to explore a
deep learning method namely LoFTR [8], which stands for detector-free local
feature matching transformer.

In the LoFTR approach, a convolutional neural network (CNN) with an encoder-
decoder architecture extracts both low-resolution and high-level features. The
LoFTR module, incorporating self-attention and cross-attention blocks,
transforms these features, and a differential matching layer offers two
methods: optimal transport and dual softmax [14] . Our implementation used a
pretrained model from the Aachen Day Night Dataset [15]. LoFTR succeeded in
generating keypoint matches for similar intensity profiles but faced
challenges in other cases similar to SIFT.These findings emphasize the need
for keypoints exhibiting consistency, accurate matching, and intensity profile
invariance across the dataset.

Fig. 2: SIFT [3] and LoFTR [8] performances in keypoint detection. The left
pair of images depicts results for SIFT, while the right one illustrates the
performance of LoFTR. It is evident that LoFTR outperforms SIFT in keypoint
detection and matching. Nevertheless, LoFTR encounters challenges when the
source and target images exhibit varying intensity distributions.

###  3.4 Radiomic features-based keypoint detection

Having faced challenges with keypoint detection using SIFT and LoFTR, which
proved sensitive to intensity variations, we shifted to our proposed radiomic
features-based method. In this section, we prove the importance of the
detected keypoints by showing our results on registering brain images among
different subjects. we also showcase the applications of the dataset
containing the original OASIS scans and the corresponding detected radiomic
keypoints on automatic keypoint detection. And finally train a GNN-based image
matcher on the same dataset.

####  3.4.1 Image registration

To assess the significance of radiomics keypoints as landmarks, we integrated
them into the loss function of a registration network (TransMorph [16]). Our
findings revealed that incorporating the keypoints’ loss led to a notable 3%
enhancement in the registration performance.

Vision transformers, excelling in capturing long-range spatial relationships,
prove effective in medical image tasks due to their large receptive fields.
TransMorph [16], a hybrid Transformer-ConvNet model, utilizes these advantages
for volumetric medical image registration. The encoder divides input volumes
into 3D patches, projecting them to feature representations through linear
layers. Sequential patch merging and Swin Transformer blocks follow. The
decoder, with upsampling and convolutional layers, connects to the encoder
stages via skip connections, producing the deformation field. We contributed
by designing a customized loss function for keypoints, utilizing Gaussian-
blurred keypoints to create a ground truth heatmap. Combining Dice and
inverted Dice losses addressed imbalanced masks, resulting in a 3% Dice score
improvement over the OASIS test dataset. TransMorph achieved a dice score of
0.89 with keypoint loss, compared to 0.86 without, underscoring keypoints’
role in enhancing registration performance. A potential avenue for future
research involves developing a loss function that considers both the feature
descriptors of keypoints and the disparity between the locations of registered
keypoints and their counterparts in the target image.

####  3.4.2 Automated keypoint detection

SuperRetina, introduced in [7] is an apdaptive version of SuperPoint model [5]
for identifying important keypoints in retinal images. Utilizing a semi-
supervised learning framework, SuperRetina maximizes the utility of limited
labeled retinal image data by combining both supervised and unsupervised
techniques. Yet, its utilization requires an initial set of ground truth
keypoints to initiate the process, subsequently increasing the detected
keypoints iteratively. In our approach, we use our radiomic keypoints as the
initial sets for OASIS images.

LK-SuperRetina [10] which is a modified version of SuperRetina consists of an
encoder for downsampling, along with two decoders—one for keypoint detection
and another for descriptor generation. Keypoint detection utilizes a mix of
labeled and unlabeled data, while descriptor training employs self-supervised
learning.

Following the U-Net [17] design, LK-SuperRetina’s shallow encoder begins with
a single convolutional layer, followed by three blocks containing two
convolutional layers, a $2\times 2$ max-pooling layer, and ReLU activation.
The keypoint decoder has three blocks with two convolutional layers, ReLU
activation, and concatenation block. The detection map ($P$) is generated
through a convolutional block with three convolutional layers and a sigmoid
activation.

The loss function combines the detector and the descriptor losses as shown in
Equation 1.

$l_{det}=l_{clf}+l_{geo}$ (1)

The classification loss component ($l_{clf}$) is defined in Equation 2, where
$\tilde{Y}$ represents the smoothed version of the binary ground truth labels
$Y$ of the keypoints after blurring them with a 2D Gaussian.

$l_{clf}(I;Y)=1-\frac{2.\sum_{i,j}(P\circ\tilde{Y})_{i,j}}{\sum_{i,j}(P\circ
P)_{i,j}+\sum_{i,j}(\tilde{Y}\circ\tilde{Y})_{i,j}}$ (2)

When feeding both the image $I$ and its augmented version $I^{\prime}$ to the
network, two tensors for the descriptors $D$ and $D^{\prime}$ are obtained.
For each keypoint $(i,j)$ in the non-maximum suppressed keypoint set
$\tilde{P}$, two distances are computed: $\Phi_{i,j}^{rand}$ between the
descriptors of $(i,j)$ in the set $\tilde{P}$ and a random point from the
registered heatmap $H(\tilde{P})$, and $\Phi_{i,j}^{hard}$ representing the
minimal distance, as depicted in Equation 3.

$l_{des}(I,H)=\sum_{(i,j)\in\tilde{P}}max(0,m+\Phi_{i,j}-\frac{1}{2}(\Phi_{i,j}^{rand}+\Phi_{i,j}^{hard}))$
(3)

For more in-depth information on the loss function, please refer to the
SuperRetina paper [7]. Figure 3 shows the results obtained from LK-
SuperRetina. As demonstrated, the number of additionally identified keypoints
meets expectations, showcasing the network’s proficiency in capturing good new
keypoints. The network successfully detects both the ground truth keypoints
and extra keypoints during the testing phase.

Figure 3 presents two instances demonstrating the resilience of the keypoint
detection model against deformations. The images were deformed randomly by an
affine deformation. We passed the original image (target) and the deformed
image (reference) seperately to LK-SuperRetina. The model successfully
identified corresponding keypoints in both images, as indicated by the number
of good matches. We adjusted the thresholds of LK-SuperRetina to detect a
smaller set of keypoints for the clarity of the visualization, but in
practice, the model can detect over 300 good keypoints.

Fig. 3: Two examples show the robustness of the keypoint detection model. In
both rows, the sequence of images, from left to right, includes: the reference
image, the target image, and the registration output. The good matches
contribute to aligning the images effectively.

####  3.4.3 Image matching

Following the detection of keypoints within the brain images, we proceed to
construct graphs to be used as inputs of the GNN, for accomplishing specific
tasks such as matching in our study. Within this paper, we demonstrate our
success in training a GNN-based matcher (SuperGlue [11]) using the graphs
formed from the detected keypoints.

SuperGlue designed for matching two sets of local features by identifying
correspondences and filtering non-matchable points. Using attention-based
graph neural networks, it integrates context aggregation, matching, and
filtering within a unified architecture. SuperGlue employs self-attention to
enhance the receptive field of local descriptors and cross-attention for
cross-image communication. The network handles partial assignments and
occluded points by solving an optimal transport problem. With superior
performance over other learned approaches, SuperGlue achieves state-of-the-art
results in pose estimation for challenging real-world indoor and outdoor
environments.

Figure 4 and Table 1 show a performance comparison between the brute force
matcher and SuperGlue across the test dataset. SuperGlue enhances the features
of the detected keypoints which improves the matching performance. The
identification of keypoints and subsequent graph creation broadens the
application of GNNs to various tasks in brain image analysis, extending beyond
image matching.

Fig. 4: A comparison of the SuperGlue and brute force matcher performance in matching detected keypoints on brain images. Method | Avg. No. Good Matches Confidence Score
---|---
BF [18] | 7 $0.449\pm 0.007$
SuperGlue |  20 $\textbf{0.988}\pm\textbf{0.010}$
Table 1: A comparison between brute force matcher and SuperGlue. SuperGlue
outperforms the brute force matcher in terms of both evaluation metrics: the
average number of good matches and the average confidence score across the
entire test dataset.

##  4 Conclusion

To sum up, our radiomic keypoint detection algorithm provides a solution for
automated keypoint detection in MRI scans, overcoming challenges encountered
by traditional and other deep learning methods. The limited set of radiomic
keypoints facilitates training SuperRetina for increased keypoint detection.
These keypoints are consistent and deformation resilient.

Our approach paves the way for the application of GNN-based models on brain
images, offering a faster and more parameter-efficient alternative compared to
CNNs and transformers. Moreover, the detection of keypoints contributes to
various tasks, including registration as justified in this work.

## References

  * [1] Xiaoxin He, Bryan Hooi, Thomas Laurent, Adam Perold, Yann LeCun, and Xavier Bresson,  “A generalization of vit/mlp-mixer to graphs,”  in International Conference on Machine Learning. PMLR, 2023, pp. 12724–12745.
  * [2] Bhakti Baheti, Diana Waldmannstetter, Satrajit Chakrabarty, Hamed Akbari, Michel Bilello, Benedikt Wiestler, Julian Schwarting, Evan Calabrese, Jeffrey Rudie, Syed Abidi, et al.,  “The brain tumor sequence registration challenge: establishing correspondence between pre-operative and follow-up mri scans of diffuse glioma patients,”  arXiv preprint arXiv:2112.06979, 2021.
  * [3] Tony Lindeberg,  “Scale invariant feature transform,”  2012\.
  * [4] Peter Hviid Christiansen, Mikkel Fly Kragh, Yury Brodskiy, and Henrik Karstoft,  “Unsuperpoint: End-to-end unsupervised interest point detector and descriptor,”  arXiv preprint arXiv:1907.04011, 2019.
  * [5] Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich,  “Superpoint: Self-supervised interest point detection and description,”  in Proceedings of the IEEE conference on computer vision and pattern recognition workshops, 2018, pp. 224–236.
  * [6] Prune Truong, Stefanos Apostolopoulos, Agata Mosinska, Samuel Stucky, Carlos Ciller, and Sandro De Zanet,  “Glampoints: Greedily learned accurate match points,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 10732–10741.
  * [7] Jiazhen Liu, Xirong Li, Qijie Wei, Jie Xu, and Dayong Ding,  “Semi-supervised keypoint detector and descriptor for retinal image matching,”  in Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXI. Springer, 2022, pp. 593–609.
  * [8] Jiaming Sun, Zehong Shen, Yuang Wang, Hujun Bao, and Xiaowei Zhou,  “Loftr: Detector-free local feature matching with transformers,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2021, pp. 8922–8931.
  * [9] Daniel S Marcus, Tracy H Wang, Jamie Parker, John G Csernansky, John C Morris, and Randy L Buckner,  “Open access series of imaging studies (oasis): cross-sectional mri data in young, middle aged, nondemented, and demented older adults,”  Journal of cognitive neuroscience, vol. 19, no. 9, pp. 1498–1507, 2007.
  * [10] Sahar Almahfouz Nasser, Nihar Gupte, and Amit Sethi,  “Reverse knowledge distillation: Training a large model using a small one for retinal image matching on limited data,”  arXiv e-prints, pp. arXiv–2307, 2023.
  * [11] Paul-Edouard Sarlin, Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich,  “Superglue: Learning feature matching with graph neural networks,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2020, pp. 4938–4947.
  * [12] Joost JM Van Griethuysen, Andriy Fedorov, Chintan Parmar, Ahmed Hosny, Nicole Aucoin, Vivek Narayan, Regina GH Beets-Tan, Jean-Christophe Fillion-Robin, Steve Pieper, and Hugo JWL Aerts,  “Computational radiomics system to decode the radiographic phenotype,”  Cancer research, vol. 77, no. 21, pp. e104–e107, 2017.
  * [13] Ali Hatamizadeh, Vishwesh Nath, Yucheng Tang, Dong Yang, Holger R Roth, and Daguang Xu,  “Swin unetr: Swin transformers for semantic segmentation of brain tumors in mri images,”  in International MICCAI Brainlesion Workshop. Springer, 2021, pp. 272–284.
  * [14] John Bridle,  “Training stochastic model recognition algorithms as networks can lead to maximum mutual information estimation of parameters,”  Advances in neural information processing systems, vol. 2, 1989.
  * [15] Torsten Sattler, Will Maddern, Carl Toft, Akihiko Torii, Lars Hammarstrand, Erik Stenborg, Daniel Safari, Masatoshi Okutomi, Marc Pollefeys, Josef Sivic, et al.,  “Benchmarking 6dof outdoor visual localization in changing conditions,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 8601–8610.
  * [16] Junyu Chen, Eric C Frey, Yufan He, William P Segars, Ye Li, and Yong Du,  “Transmorph: Transformer for unsupervised medical image registration,”  Medical image analysis, vol. 82, pp. 102615, 2022.
  * [17] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18. Springer, 2015, pp. 234–241.
  * [18] David G Lowe,  “Distinctive image features from scale-invariant keypoints,”  International journal of computer vision, vol. 60, pp. 91–110, 2004.


# Can machines solve general queueing systems?

Eliran Sherzer, Arik Senderovich, Opher Baron and Dmitry Krass

(2022; Date: February 2022)

###### Abstract.

Computing performance measures in _general_ queueing systems is a challenging
task: most complex systems are intractable, while simpler systems often
require involved mathematical analysis to derive exact solutions. For more
complex systems only approximation formulas are available, which may lead to
large errors. While a common practical alternative is simulation modeling, a
simulation model can only be designed for a specific, not a general system.
Moreover, such simulation models may take a long time to converge to an
equilibrium solution.

In this paper, we analyze how well a machine can solve a general problem in
queueing theory. To answer this question, we use a _deep learning_ model to
predict the stationary queue-length distribution of an $M/G/1$ queue (Poisson
arrivals, general service times, one server). To the best of our knowledge,
this is the first time a machine learning model is applied to a general
queueing theory problem. We chose $M/G/1$ queue for this paper because it lies
“on the cusp” of the analytical frontier: on the one hand exact solution for
this model is available, which is both computationally and mathematically
complex. On the other hand, the problem (specifically the service time
distribution) is general. This allows us to compare the accuracy and
efficiency of the deep learning approach to the analytical solutions.

The two key challenges in applying machine learning to this problem are (1)
generating a diverse set of training examples that provide a good
representation of a “generic” positive-valued distribution, and (2)
representations of the continuous distribution of service times as an input.
To deal with the first challenge we propose a novel method, based on phase-
type distributions, to sample service distributions such that our data will
“widely cover” the space of non-negative distributions. We use the results
from queuing theory to obtain exact solutions of $M/PH/1$ models (phase-type
service distribution).

To overcome the second challenge, we use the arrival rate and the first $n$
moments of the service times as inputs to our deep learning model. The output
(that we learn) is the stationary queue-length distribution. To test our
method, we conduct an extensive empirical evaluation. We use the sum of
absolute errors between our predictions and the “ground truth” as our primary
metric for evaluating the accuracy.

Our results show that our model is indeed able to predict the stationary
behavior of the $M/G/1$ queue extremely accurately: the average value of our
metric over the entire test set is $0.0009$; other metrics we employ also show
an excellent predictive accuracy. Moreover, our machine learning model is very
efficient, computing very accurate stationary distributions in a fraction of a
second (an approach based on simulation modeling would take much longer to
converge). We also present a case-study that mimics a real-life setting and
shows that our approach is more robust and provides more accurate solutions
compared to the existing methods. This shows the promise of extending our
approach beyond the analytically solvable systems (e.g., $G/G/1$ or $G/G/c$).

One interesting insight yielded by our approach is that the number of moments
of the service time distribution required for accurately capturing the
dynamics of an $M/G/1$ queue is quite small; high accuracy is achieved with
$n=5$, while the higher service time moments have negligible effect on the
accuracy.

††copyright: acmcopyright††journalyear: 2022††doi:
XXXXXXX.XXXXXXX††conference: ACM SIGMETRICS / IFIP PERFORMANCE 2022; June
06–10, 2022; Mumbai, India

##  1\. Introduction

In this paper we attempt to “teach” a machine to compute performance measures
in a general queueing system. In effect, we ask whether a machine can replace
complicated mathematical methods and provide an accurate solution for problems
in queueing theory. As we show in this paper the answer is likely in the
_affirmative_.

In order to investigate how well a machine can solve a general problem from
the queueing domain, we revisit the computation of the stationary distribution
of an $M/G/1$ (Poisson arrivals, general service, single server) queue, a
well-studied problem in queueing theory, and solve it using machine learning
(ML). We focus on this problem, because it lies “on the edge” of the
analytical frontier: while exact solutions for this model are available, they
are both computationally and mathematically complex. On the other hand, the
problem (specifically the service time distribution) is general. This allows
us to compare the accuracy and efficiency of the deep learning approach to the
analytical solutions.

Motivation Typically, queueing systems are complex, hard to analyze, and
analytical methods often struggle even with simple stylized models. For
example, a single-server queue with general inter-arrival and service
distributions (the $G/G/1$ queue) cannot be analyzed in an exact fashion, and
robust approximation methods are not available. Furthermore, real-life
queueing systems introduce complexities that often make queueing models
analytically intractable, e.g., customer abandonment, correlations between
arrivals and service times, and many more. Thus, in many cases, since no
analytic solution is available, the system can only be analyzed via
simulation. Yet, simulation has its own downsides, ranging from the
requirement for a detailed specification of all underlying distributions to,
often, a long convergence time, especially for multi-server systems.

In this work, instead of using the traditional queueing theory methods (e.g.,
embedded Markov-chains at departure epochs or supplementary variables), we
leverage the power of ML, using deep neural networks as our solution
framework. This is motivated by the major successes of this approach in a
variety of fields ranging from computer vision, to natural language
processing, to traffic prediction in transportation systems (see (LeCun2015, ;
726791, ; 5537907, )). Moreover, within the Operations Research (OR) domain,
recent papers have shown the usefulness of deep learning in well-known
problems such as scheduling, resource allocation, and others (see (222611, ;
8943940, )).

Technical challenges Applying a deep neural model to solve111Throughout this
paper, “solve” refers to computing the stationary queue-length distribution an
$M/G/1$ system is a non-trivial task. One faces several major challenges:

  1. (1)

Generating a set of training instances that provide a good representation of a
“generic” positive-valued distribution. More specifically, we need to obtain a
sufficiently diverse set of input pairs of arrival rates and service-time
distributions that would be representative of the space of all non-negative
distributions, while maintaining a “conformable” parametric representation.

  2. (2)

Representation of service time distributions as an input to the deep learning
model. While traditional steady-state computation methods use arrival rates
and service time distributions as direct inputs, an input to a deep learning
model must be represented by a tensor; with possibly some information being
lost as a result. What is the right distribution-to-tensor mapping?

To deal with the first challenge, we leverage several key results for Phase-
Type (PH) family of distributions. The power of PH distributions lies in the
fact that they are dense in the class of all non-negative distribution
functions. Practically speaking, this means that a PH distribution with a
sufficient number of phases can, theoretically, approximate any non-negative
distribution to an arbitrary precision (see (Asmussen2003, ), Theorem 4.2).
Finally, an $M/PH/1$ system (PH-type service time distribution) allows for a
much more computationally efficient solution methodology than a general
$M/G/1$ queue. Sampling from the universe of PH distributions is not a
straightforward task, as its parameters are not independent and have multiple
constraints that must be enforced. To overcome this challenge, we propose a
novel and versatile sampling method that covers a “wide” range of
distributions. We show the latter claim empirically in Section 4.

To address the second challenge, we develop a representation of the service
time component of the input as the first $n$ moments of the service time
distribution; these moments are computed from the empirical data set of
observed service times. The number of moments $n$ is a hyper-parameter of our
approach, which is determined after examining the accuracy of the learned
model as function of $n$.

While representing the service time distribution using a finite number of
moments may seem counter-intuitive (since different distributions can have
identical set of first $n$ moments) this representation is natural in view of
a number of existing results in stationary queueing analysis which often
depend only on the first two moments of both the service time and the inter-
arrival time distributions; for example, the mean queue length for the $M/G/1$
queue given by the Khintchine–Pollaczek formula, or the bound of the mean
queue length in the $G/G/1$ queue (see Chapter 2 in
(https://doi.org/10.1002/net.3230070308, )). The first two moments are also
widely used in more sophisticated queueing models such as queueing networks
(see (Whitt1994, ; https://doi.org/10.1111/j.1937-5956.1993.tb00094.x, )).
While it is clear that the moments of the distribution play a key role in
queueing analysis, we do not know how many moments are actually needed to have
a sufficient representation of the entire queue length distribution. This
study allows us to shed some light on this important question.

Results Our results show that a deep neural network can almost perfectly
approximate the stationary queue-length distribution for $M/G/1$ queues, while
using a relatively small net (roughly 15,000 parameters). Hence, potentially,
there is room for using deep learning techniques for more complex queueing
systems. Further, we show that the only the first $5$ moments of service time
suffice to determine the distribution of the $M/G/1$ stationary queue length;
there is no substantial gain in accuracy from adding higher-order service time
moments. This result is very intriguing as it provides further understanding
of the dynamics of an $M/G/1$ queue, and perhaps more complex systems as well.
It also greatly simplifies the task of fitting a particular parametric
distribution to the observed service time data as only the first $5$ moments
have to be matched, rather than using more complex methods such as EM
(Expectation Maximization algorithm)222Papers like (TELEK20071153, ) are doing
just that..

Note that by applying deep learning, we are able to provide a new data-driven
solution to a fundamental problem in queueing theory. To use the traditional
analytical $M/G/1$ results in practice, when one is provided with data, rather
than a closed-form distribution of service times, one must carefully choose a
parametric family of distributions, apply a statistical or an algorithmic
method in order to find the parameters that optimize some notion of goodness-
of-fit, and only then apply the analytical technique to obtain steady-state
queue-length distribution. These procedures are seldom straightforward, are
often time consuming, and do not guarantee an accurate fit in the end. Our
approach, where only the first $5$ moments are computed and then the solution
is directly produced by a pre-trained neural network is far more simple,
computationally efficient, and accurate, as demonstrated in Section 7.

Contributions The main contributions of this paper are threefold:

  * •

We provide a novel data-driven technique based on deep learning to solve a
fundamental problem in queueing theory. The network takes as input a compact
representation of the problem by using only the arrival rate and the first $n$
moments of the service time distribution.

  * •

We introduce a new sampling algorithm for PH service time distributions and
show how it can be used to generate a set of training samples. The algorithm
guarantees various important properties for the sampled PH distribution,
including absorption and stability.

  * •

We provide an extensive empirical evaluation of the approach by conducting a
set of experiments as well as analyzing a case-study that mimics a real-world
scenario.

Organization The remainder of the paper is organized as follows. In Section 2
we present an overview of the related literature. Section 3 provides a high-
level end-to-end overview of our solution. We provide our PH sampling
technique in Section 4, describe our neural architecture, input pre-
processing, and loss function in Section 5. The main findings from our
empirical evaluation are provided in Section 6, while a realistic case-study
is provided in Section 7. We conclude the paper with a discussion of the
limitations of our approach (Section 8), and some concluding remarks and
directions for future work in Section 9.

##  2\. Related Work

###  2.1. Computing the Stationary Queue-Length Distribution

There are two well-known analytic methods for computing the stationary queue
length of the $M/G/1$ queue. One is via the probability generating function
given by the Pollaczek–Khinchine transform equation (see (Harchol2013, )). The
other is via the supplementary variable method by Cox and Kendall (cox_1955,
). While both methods are exact, they can be computationally burdensome: in
order to obtain the probability of having $n$ customers in the steady-state
system, the former method requires the $n^{th}$ derivative to the generating
function, while for the latter we need to solve a set of $n$ non-linear
equations. Even if service times are uniformly distributed both methods can be
computationally heavy.

Since our training sample is obtained using PH distributions, a simpler option
is available for computing the stationary queue-length distribution: the Quasi
Birth-and-Death method (QBD) developed by Neuts (Neuts1981, ). As mentioned in
Section 1, a PH distribution with a sufficient number of phases can
approximate any non-negative distribution arbitrarily closely. However, as we
show later in the paper, it is not always easy to fit a PH representation that
provides a good approximation for a given data trace.

###  2.2. Deep Learning for Queueing Systems

There have been only a few prior studies that utilize deep learning neural
networks to analyze queues. The most similar paper to ours is by Nii et. al.
(Nii20, ) where the authors propose a deep learning model to analyze a
$GI/G/s$ queue. The input includes only the first two moments of both the
inter-arrival and the service times distribution. However, they are evaluating
only the average waiting time, rather than estimating the whole distribution,
and consider only a small range of possible values using a very small NN
(roughly 100 variable and three hidden layers). Kyritsis and Deriaz
(Kyritsis19, ) also utilized a neural network to predict waiting times.
However, unlike our model or the method proposed in (Nii20, ), they do not
refer to a steady-state distribution, but rather focus on an online prediction
setting conditioned on the current state of a network. Their results are not
sufficiently robust to make predictions on new data sets. In a similar work to
(Kyritsis19, ), Olawoyin and Hijry (Hijry_2020, ) evaluate the expected
waiting times in an emergency room using deep learning. This is quite
different from our goal of studying the ability of a machine learning approach
to analyze a general queuing problem.

Other papers that combined queueing systems and machine-learning did not
employ deep learning for evaluating queueing systems, focusing instead on
optimizing and controlling the queue. Efrosinin and Stepanova (Efrosinin, )
used deep learning for estimating of the optimal threshold policy in queues
with heterogeneous servers. Dai and Gluzman (dai2021queueing, ) used a deep
reinforcement learning model for controling queueing networks by effectively
solving a Markov decision process. For additional papers that employ machine-
learning for optimization and control of queueing networks we refer the reader
to (CHEN20113323, ; 8919665, ; 9502440, ). These papers utilize machine
learning methods to solve an optimization problem, while we are using them for
capturing the stochastic dynamics of the system.

##  3\. Solution Overview

The main phases of our solution are presented in Figure 1, where green ovals
correspond to processing/computational steps and red rectangles to data sets.

Steps 1-4 represent the generation of training input and output data; they are
described in Section 4 below. We start by generating input samples (step 1,
“Sampling Procedure” on Figure 1); see Section 4.1 for details. The sample is
generated in the form of pairs, with each pair comprising an arrival rate and
a service distribution333we only consider pairs that result in a stable
system, i.e., the arrival rate must be smaller than the service rate.. This
results in the “Training Input” data set in (step 2 on the Figure). Next we
compute the stationary queue-length distribution by applying a well-known QBD
method (step 3, “Analytic Computation” on the Figure; details can be found in
Section 4.2). This results in the training output (step 4), used as the target
when learning from the training data. The output consists of the first $l$
values of the stationary queue-length distribution; the value of $l$ is
described below.

Steps 2, 6, 7 and 8 on Figure 1 correspond to the training of the Deep
Learning model, and steps 5,6,7,8 to application of the model to new test
data. These steps are detailed in Section 5.

Prior to entering the input into the neural network some pre-processing is
required (step 6, “Pre-processing” in Figure 1). First, we compute the first
$n$ moments of the service time distribution, where the value of $n$ is a
hyper-parameter that requires tuning. Next, the computed moments undergo some
additional pre-processing (standardization). These steps are detailed in
Section 5.1. The processed input and the training output are fed into the deep
learning model for training (step 7, “Deep learning Model” on Figure 1); see
Section 5.2 for details. The model outputs the predicted stationary queue-
length distribution (step 8) in the same format as the training output.

Once the training is completes, we use test data (step 5), apply the required
pre-processing, and use our deep learning model to obtain the predicted
values.

Figure 1. Work-flow diagram of our learning procedure.

##  4\. Data Generation

In this section we outline the input sampling procedure (Section 4.1), and the
training output generation step (Section 4.2).

###  4.1. Input Sampling Procedure

To ensure the stability of the queueing system we must make sure the arrival
rate is smaller than the corresponding service rate. To achieve this, it is
simpler to sample the service time distribution first, and subsequently obtain
the arrival rate. The latter is drawn uniformly over the stability range of
the sampled service rate. The main challenge here is to sample many service
distributions so that they provide a “wide” coverage of the state-space of all
positive distributions. To this end, we employ phase-type (PH) distributions;
as mentioned earlier, PH distribution can approximate any non-negative
distribution (Neuts2013, ).

Phase-Type Distribution A PH distribution can be represented by a random
variable describing the time until absorption of a continuous-time Markov
chain (CTMC) with a single absorbing state; all other states are transient.
Each of the transient states of the CTMC corresponds to one of the states of
the PH distribution, where a state sojourn time follows an exponential
distribution. For a full representation of PH with $m\geq 1$ states, we need
to specify two parameters, namely $\alpha$, the initial probability $m$-vector
of starting the process in any of the $m$ states, and $S$, the corresponding
generating $m\times m$ matrix of the CTMC. We restrict our method to PH
distributed random variables such that $\sum_{i}|\alpha_{i}|=1$444Implying
that we do not allow the system to start in the absorbing state.. In order for
the pair $(\alpha,S)$ to be a valid pair for the CTMC, we need to meet the
following constraints:

  * •

$S[i,i]<0,\forall i\leq m$.

  * •

$\sum_{j=1}^{m}S[i,j]\leq 0$.

Challenges in Sampling from PH For a general input, we wish to control two
degrees of freedom hidden in $\alpha$ and $S$: (1) the size (number of states)
of the PH distribution, and (2) the transition function that governs the
switches between the different states (to be explained below). The former is
relatively straightforward to deal with as we uniformly sample a value between
$1$ and $max_{ph}$, where $max_{ph}$ is the maximum PH size that we use for
our experiments (this is a hyper-parameter of our approach and can be tuned by
the user). The second consideration is not trivial and requires some effort.

When generating $\alpha$ and $S$ we wish to allow all possible transition
structures from state to state. The transition structure is determined by
three elements: (i) the initial probabilities vector $\alpha$, (ii) the
locations of the non-diagonal cells in $S$ in which the value is $0$, and
(iii) the sum of each row in $S$ . Note that iff $S[i,j]=0$, for $i\neq j$,
then it is impossible to directly reach state $j$ from state $i$. Clearly, if
all non-diagonal cells in a single row $i\leq m$ are $0$, then the underlying
CTMC will reach its absorption state once the system exits state $i$. We refer
to such states as fully-absorbing. These are controlled by element (ii).
Letting $s_{i}$ be the sum of row $i$ in $S$, i.e.,
$\sum_{j=1}^{m}S[i,j]=s_{i}$, if $s_{i}<0$, then with probability
$s_{i}/S[i,i]$ the CTMC will reach its absorption state after exiting state
$i$. We refer to such states as partially absorbing. However, if $s_{i}=0$,
then the CTMC cannot reach the absorption state once the system exits state
$i$. We refer to such states as non-absorbing. The composition of partially
absorbing and non-absorbing states is controlled by element (iii). Thus, in
order to control the structure of the PH distribution we need to control for
elements (i), (ii), and (iii).

Note that we can have different values of $\alpha$ and $S$ and yet still end
up having the same service time distribution. Moreover, in queueing systems,
the time scale does not matter. What matters is only the ratio between inter-
arrival and service times. Thus, when having two different distributions, they
may still be very similar once we put the two on the same time scale. In order
to avoid such cases and sample distributions as diverse as possible, we first
split states into “classes”.

Classes of States in PH We define a class as a set of states such that if we
reach one of them from the initial state, it is impossible to reach any state
outside this set before the underlying CTMC reaches its absorption state. Each
class, in fact, is a PH distribution on its own. Thus, we can sample each
class separately and then group them together for a singe PH representation.
The main reason for doing so is that if we have a single large class, even
with different structures of transition matrices, we will end up in many cases
with similar distributions (Neuts2013, ). Since class structures are crucial
for the generality of our result we shall first sample the number of classes
and their sizes and only then sample the specific states within each class.

The total number of classes may vary from 1 to $m$ as it cannot exceed the
total number of states. We note that the PH random variable (that corresponds
to the time until the CTMC is absorbed) has, by its definition, a finite
expected value. In order to avoid infinite expected values, each class must
adhere to two conditions: (1) having at least a single partially or fully
absorbing state, and (2) any non-absorbing state must have a positive
probability to reach one of the partially or fully absorbing states.
Otherwise, the system may enter a “livelock”, i.e., a loop in which it never
reaches the absorbing state. We enforce condition (1) in our sampling
methodology; however, enforcing condition (2) may increase the computational
complexity of our sampling approach. Instead, once the procedure is complete,
we check for finite expected value. If it the value is infinite, we reject the
sample. We note that these rejections occurred in our extensive experiment
only 5% of the times (see Section 6 for details on the experiment).

Algorithm 1 Sample service time distribution

INPUT: $max_{ph}$

OUTPUT: $\alpha$, $S$

1:procedure Sampling PH

2: $m\leftarrow U(1,max_{ph})$

3: $NumClasses\leftarrow U(1,m)$ $\triangleright$ Sample number of different
classes.

4: $P\leftarrow Dirichlet(a_{1},a_{2},\ldots,a_{NumClasses})$ where
$a_{j}\leftarrow U(0,1),\forall j\leq NumClasses$

5: $AssignedStates\leftarrow 0$

6: for $i=1,i\leq NumClasses,i++$ do

7: $Size_{i}\leftarrow U(1,m-AssignedStates-(NumClasses-i))$

8: $AssignedStates\leftarrow AssignedStates+Size_{i}$

9: $\alpha^{i},S^{i}\leftarrow$ Sampling Class $i$ ($Size_{i}$)
$\triangleright$ Sample PH distribution for each class

10: $i_{low}=1+\sum_{j=1}^{i-1}Size_{j}$ $\triangleright$ $i_{low}$ is the
lower index of the $i^{th}$ class in $\alpha$ and $S$

11: $i_{up}=1+\sum_{j=1}^{i}Size_{j}$ $\triangleright$ $i_{up}$ is the upper
index of the $i^{th}$ class in $\alpha$ and $S$

12: $S[i_{low}:i_{up},i_{low}:i_{up}]=S^{i}$

13: $\alpha[i_{low}:i_{up}]=\alpha^{i}P[i]$

14: end for

15: if $E[X]<\infty$, where $X\sim PH(\alpha,S)$ then

16: Return $\alpha$, $S$

17: end if

18:end procedure

Algorithm 2 Sample PH of a single class

INPUT: $Size_{i}$

OUTPUT: $\alpha^{i}$, $S^{i}$

1:procedure Sampling class $i$($Size_{i}$)

2: $\alpha^{i}\leftarrow Zeros(1,Size_{i}),S^{i}\leftarrow
Zeros(Size_{i},Size_{i})$, $Trans\leftarrow Zeros(Size_{i},Size_{i})$

3: $FullA_{i},PartA_{i},NonA_{i}\leftarrow$ Sampling state types $(Size_{i})$
$\triangleright$ Assigning states to different types

4: $Trans\leftarrow$ Sample Trans($Trans,Size_{i},PartA_{i},NonA_{i}$)

5: $S^{i}[j,j]\leftarrow-U(1,1000),\forall j\\{j\leq Size_{i}\\}$

6: $S^{i}[j,k]\leftarrow U(1,1000),\forall j,k\\{j\neq k,Trans[j,k]=1\\}$

7: $p_{j}\leftarrow U(0,1)$, $\forall j\in PartA_{i}$ $\triangleright$ $p_{j}$
is the absorbing probability of state $j$

8: $S^{i}[j,k]=S^{i}[j,k](1-p_{j})\frac{-S^{i}[j,j]}{\sum_{l=1,l\neq
j}^{Size_{i}}S^{i}[j,l]}$ ,$\forall j\in PartA_{i},\forall k\neq j$

9: $S^{i}[j,k]=S^{i}[j,k]\frac{-S^{i}[j,j]}{\sum_{l=1,j\neq
j}^{Size_{i}}S^{i}[j,l]}$, $\forall j\in PartA_{i}\forall,k\neq j$

10: $\alpha^{i}\leftarrow Dirichlet(a_{1},\ldots,a_{Size_{i}})$, $a_{j}\sim
U(0,1),\forall j\leq Size_{i}$

11: Return $\alpha^{i}$, $S^{i}$

12:end procedure

13:procedure Sample Trans($Trans,Size_{i},PartA_{i},NonA_{i}$)

14: $NumTrans_{j}\leftarrow 1+Binomial(Size_{i}-2,U(0,1))$, $\forall j\in
PartA_{i},NonA_{i}$

15: $Seleced_{j}\leftarrow Choose(\\{k\leq Size_{i},k\neq j\\},NumTrans_{j})$
, $\forall j\leq Size_{j}$

16: $Trans[j,k]=1,$ $\forall j\leq Size_{i},k\in Seleced_{j}$

17: Return $Trans$

18:end procedure

Algorithm 3 Assign states into different classes

INPUT: $max_{ph}$

OUTPUT: $FullA_{i}$, $PartA_{i}$, $NonA_{i}$

1:procedure Sampling state types($Size_{i}$)

2: $FullA_{i}\leftarrow\\{\\}$, $PartA_{i}\leftarrow\\{\\}$,
$NonA_{i}\leftarrow\\{\\}$

3: $p_{1},p_{2},p_{3}\leftarrow Dirichlet(a_{1},a_{2},a_{3})$, where
$a_{j}\leftarrow U(0,1),j=1,2,3$

4: for $j=1,j\leq Size_{i},j++$ do

5: if $j=1$ then $\triangleright$ To ensure at least one state in $FullA_{i}$
or $PArtA_{i}$

6: $state\leftarrow Categorical(p_{1}/(p_{1}+p_{2}),p_{2}/(p_{1}+p_{2}))$

7: else

8: $state\leftarrow Categorical(p_{1},p_{2},p_{3})$

9: end if

10: Assign state to $FullA_{i}$ if $state=1$, $PartA_{i}$ if $state=2$,
$NonA_{i}$ if $state=3$,

11: end for

12: Return $FullA_{i}$, $PartA_{i}$, $NonA_{i}$

13:end procedure

Sampling Algorithms The sampling procedure is described in our main algorithm,
Algorithm 1, while Algorithms 2 and 3 are procedures called from within
Algorithm 1.

The idea of Algorithm 1 is partitioning the PH into classes, sample each class
independently (which takes place in Algorithm 2) and then aggregate all
classes into a unified PH representation $(\alpha,S)$. We first sample the
number of states $m$555Choosing $1000$ was a design choice. In principle, this
parameter can be chosen by the user. of the PH, the number of classes
($NumClasses$) and the initial probability of each class ($P$). Once those are
obtained, we sample the PH representation of class $i$ (i.e., $\alpha^{i}$ and
$S^{i}$), $\forall i\leq NumClasses$, given its size $(Size_{i})$. This
sampling procedure takes place in Sampling Class $i$ (i.e., the procedure in
Algorithm 2). We note that sampling $Size_{i}$, should be done carefully, as
the constraint $\sum_{i=1}^{NumClasses}Size_{i}=m$ must be met. For this
purpose, $Size_{i}$ is sampled uniformly over 1 to
$m-AssignedStates-(NumClasses-i)$. The upper-bound ensures that we assign only
from the remaining states, (i.e., that are not assigned by the previous
classes: $1,\ldots,i-1$), and that there will be at least one state for all
future classes (i.e,. states $i+1,\ldots,NumClasses$).

Once $\alpha^{i}$ and $S^{i}$ are sampled, we assign their values according to
their designated locations in $\alpha$ and $S$, as described in lines 12 and
13. We note that for $S$ the assignment are straightforward, while in $\alpha$
we multiply them by $P[i]$. This is due to the fact that “Sampling class $i$”
returns $\alpha^{i}$, which is summed up to $1$. However, the total
probability of going into class $i$ is $P[i]$ and hence $\alpha^{i}$ is
multiplied by $P[i]$. Finally, as mentioned above it is possible that the
underlying CTMC will reach a livelock, hence we return $\alpha$ and $S$ only
if the expected value of the sampled service time is finite.

###### Remark 1.

Note that we sample from the Dirichlet distribution with random weights
throughout our sampling procedure (see also Algorithm 2 and 3). This is a
useful way to sample a vector of probabilities (i.e., a vector that sums to 1)
since it allows us to get unbalanced, as well as balanced, realizations. For
example, if we sample $(a_{1},a_{2},a_{3})$ from $U(0,1)$ and then normalize
the weights to sum to 1, it is very unlikely to get an unbalanced realization
like $(0.98,0.01,0.01)$. However, we want realizations that are as diverse as
possible. When using a Dirichlet distribution, realization such as these are
not that scarce.

Next, we provide details of the two sampling procedures used by Algorithm 1,
first with Algorithms 2 and then 3. Algorithm 2 samples the $i^{th}$ class in
our CTMC. The logic of the sampling procedures is as follows. We first
classify each state and then determine all possible direct transitions.
Naturally, the next step should sampling the rates where the transitions are
possible. However, within any row of $S^{i}$, the transition rates are not
independent. For example, if state $j$ is non-absorbing, then the sum of the
$j^{th}$ row must be 0. In order to avoid the complexity of the conditional
sampling we take a different approach, which is first sample without
conditioning and then normalize the rates accordingly such that all
constraints are met. Once $S^{i}$ is determined we sample $\alpha^{i}$, which
is done according to Dirichlet distribution.

For the sake of clarity, we would like to point out that each state $j\leq
Size_{i}$ is classified to be either full-absorbing, partial-absorbing and
non-absorbing and assigned to $FullA_{i}$, $PartA_{i}$ and $NonA_{i}$,
respectively. $Trans$ is a $(Size_{i},Size_{i})$ matrix where $Trans[j,k]=1$
only if a direct transition from $j$ to $k$ is possible. $Trans$ is sampled in
‘Sample Trans’, where the number of potential states that can be reached
directly from state $j$ are $Size_{i}-1$. However, if $j\in
PartA_{i},NonA_{i}$ then there has to be at least one state $k\neq j$ to
transfer to. This is taking care of in line 14 of Algorithm 2. Once $Trans$ is
fixed, we sample both diagonal and non-diagonal states, and then normalize the
rates such that they will meet the constraints according to their type. For
full-absorbing states nothing is requires as there is only one non-zero value,
which is on the diagonal. For normalizing the rates of partial-absorbing and
non-absorbing states, the idea is similar in both cases: $S^{i}[j,k]$ is being
multiplied by a constant $\frac{-S^{i}[j,j]}{\sum_{j=1,j\neq
k}^{Size_{i}}S^{i}[j,k]}$, which makes sure that the sum of all non-diagonal
values will sum to the diagonal value. For partial-absorbing states we also
multiply by $1-p_{j}$, whereas, $p_{j}$ is the probability that state $j$ will
reach the absorbing state upon departure. This is because, we wish that the
ratio between the sum of non-diagonal entries and the absolute value of the
diagonal entry will be the probability that the process does not reach its
absorbing state once it departs state $j$.

We next provide details of Algorithm 3. Here, we split all states
$(1,\ldots,Size_{i})$ into $FullA_{i}$, $PartA_{i}$ and $NonA_{i}$, which
represent the full-absorbing, partial-absorbing and non-absorbing sets,
respectively. It is important that at least one of the sets $FullA_{i}$ and
$PartA_{i}$ will be non-empty (otherwise, it is impossible to reach the
absorbing state). We commence by sampling the probabilities $p_{1}$, $p_{2}$
and $p_{3}$, which represent the probability that a state assigned to
$FullA_{i}$, $PartA_{i}$ and $NonA_{i}$, respectively (line 3). Then we
iterate over all states, where the first state must be in either $FullA_{i}$
or $PartA_{i}$. This is done using the Categorical distribution. That is, a
value is $i$ with probability $p_{i}$, for $i=1,2,3$. In line 6, the
probabilities of the categorical distribution are normalized to sum up to 1.
Once the variable $state$ is sampled, we assign state $j$ to its designated
set in line 10.

Illustration of the Results. We empirically show that provided enough samples,
our approach produces a set of phase-type distributions that provides a “wide
cover” of the positive-valued distributions. For the sake of illustration, and
without the loss of generality, we scale all service time distributions such
that the mean value is 1. As shown in Figure 2, the PH distributions generated
by the Algorithms 1-3 above can take on various forms and capture complex
structures.

Figure 2. An example of 15 service time distributions sampled using Algorithm
1.

For further empirical evidence, we generated a data set containing 60,000
samples, which was later used as our test set. The mean value was set to $1$.
We present the range of the service times for 6 different percentiles (see
Table 1). The purpose of this table is to demonstrate that we are indeed
generating a very diverse set of distributions with our sampling method. For
each percentile level we provide the minimum and the maximum service time in
our data set. For example the “Lower” value of $8.511e^{-6}$ and the “Upper”
value of $1.0715$ in the $25\%$ percentile column shows the very wide range of
service times corresponding to this percentile in our sample of 60,000 service
time distributions. Since we know that all distributions have a mean value of
1, the very low and high values shows a large variability in where the
probability mass is concentrated among the distributions in our sample. This
indicates a wide range of the types of distributions included in our sample,
which was the goal of our sample generating process.

Percentile | 25% | 50% | 75% | 90% | 99% | 99.9 %
---|---|---|---|---|---|---
Lower | $8.5111e^{-6}$ | 0.0003 | 0.0027 | 0.0308 | 0.8391 | 1.1684
Upper | 1.0715 | 1.5812 | 3.0376 | 7.7300 | 29.9446 | 29.9871
Table 1. Range of percentiles.

Once the service time distribution (a PH) is sampled, we scale it such that
its expected value is $1$. To conclude input generation procedure, we sample
arrival rates from $U(0,\rho_{max})$ (in our experiments, we chose
$\rho_{max}=0.95$ to ensure that we do not cover very heavily loaded systems).
Note that our scaling of mean service times to $1$ and arrival rates to
$(0,0.95)$ does not result in any loss of generality with respect to the
steady state probabilities for systems with the utilization ratio in this
range.

###### Remark 2.

High utilization values, e.g., larger than 0.95, change the queue dynamics
significantly, increasing the probability of very high queue lengths. If
prediction for systems with very high utilization rates is desired, the model
should be trained only on samples corresponding to very high utilization
values of the queue.

###  4.2. Generating Output

To generate the output (stationary queue-length distribution), we use the
arrival rate and PH representation of service time distribution (per input) to
compute the first $l$ steady-state probabilities using the QBD method
(Neuts1981, ), which is significantly faster than the other analytical methods
for the case of PH-distributed service times.

Since we need to feed the computed distribution to the deep learning model in
the form of a fixed vector, we truncate our computation at $l$ such that the
total probability of having more than $l$ customers is smaller then $\epsilon$
(these can be viewed as meta-parameters). In our empirical evaluation we used
$\epsilon=10^{-9}$ for which we achieved the desired performance with $l=70$
in all generated samples. Since we also consider the probability of having an
empty system, this corresponds to covering the probability of up to 69
customers in the system.

###### Remark 3.

In term of run-times, sampling a single service time distribution and
computing the output takes on average 0.56 seconds. For the sizes of training,
validation, and test data sets provided in Table 2, sampling the entire data
took less than 7 hours using a cluster with 30 nodes. Each node uses Intel
Xeon Gold 5115 Tray Processor, with RAM memory of 128 GB.

##  5\. Deep Learning

In this section, we provide details into our deep learning model including our
choice of network architecture and the loss function. We start by describing
the input pre-processing step in Section 5.1 and then the network architecture
in Section 5.2.

###  5.1. Input pre-processing

Once the sampling procedure in Algorithm 1 is completed we end up with pairs
of arrival rate and service time distributions (in the form of $\alpha$ and
$S$). Hypothetically, we could enter $\alpha$ and $S$ along with their
corresponding arrival rates directly into the model. However, we have made a
design choice not to do so as it leads to very sparse and inefficient
representations that can have more than $10,000$ values. Instead, we decided
to compute the first $n$ moments of the service time distribution, and use it
as an input into the deep learning model together with the sampled arrival
rate.

Note that it is common practice in deep learning to normalize the input so
that all input values will be within the same range. We wish to enforce range
similarity (1) between samples and (2) within a sample (i.e., arrival rate and
service moments). While scaling all service time distributions so that the
expected value is 1, achieves between-sample similarity, as we increase the
moments order they increase exponentially, violating (2). Thus we apply a
$log$ transform to all service time moments. This resolves (2). Since the
first moment in all cases is equal to $1$, we omit it.

###  5.2. Network Architecture

We deploy a Fully-Connected (FC) feed-forward neural network as depicted in
Figure 3. While other architecture types, such as Convolution Neural Nets
(CNN) and Recurrent Neural Nets (RNN) can certainly be considered, we felt
neither one was appropriate for this case: CNN is not well suited to our needs
as the input size is only a small vector, while RNN is more suited for time-
dependant, rather than stationary, systems.

We used the Rectified Linear Unit (ReLU) (DBLP:journals/corr/abs-1803-08375, )
as the activation function for all hidden layers. Furthermore, we applied the
Softmax function to guarantee that the output-layer weights sum up to 1 (while
this is not a common use of the Softmax function, it is quite convenient in
our case as our output must sum to 1). An exhaustive tuning of all the related
hyper-parameters is out of the scope of this paper, and commonly used default
values for the architecture of the neural network were used throughout our
experiments.

Figure 3. Diagram of our neural network.

Loss function We now present the loss function for our network. For the
training batch size $B$ (a hyper-parameter), let $Y$ and $\hat{Y}$ present a
matrix of size $(B,l)$ of the true values of the stationary queue length
distribution and the network prediction, respectively. Our loss function
$Loss(\cdot)$ is given by

(1)
$Loss(Y,\hat{Y})=\frac{1}{B}\sum_{i=1}^{B}\sum_{j=0}^{l-1}|Y_{i,j}-\hat{Y}_{i,j}|+\frac{1}{B}\sum_{i=1}^{B}max_{j}(|Y_{i,j}-\hat{Y}_{i,j}|),$

where the values $Y_{i,j}$ and $\hat{Y}_{i,j}$ refers to the $i^{th}$ sample
in the batch and the probability of having $j$ customers in the system for the
respective distributions, and, as discussed in the previous section, $l$ is
the output length, which is fixed to 70. From equation (1) we observe that the
loss function is comprised of two terms. The first term aims to minimize the
overall absolute distances between the true and predicted distributions, while
the second term seeks to minimize the maximum distance. This functional form
is motivated by typical queueing applications where we prefer a few smaller
errors to a a single large one (particularly since large errors tend to occur
in the tail, which may be particularly important for providing service level
guarantees). While the are many other measures of distance between two
distributions, we felt they were less suitable for our purposes. For instance,
cross-entropy is more suitable for zero-one classification tasks, while KL-
divergence is less suitable for distributions having numerous zero values,
which is a common scenario in stationary queue-length distributions.

The loss function refers to a single batch size denoted by $B$, a hyper-
parameter that is be tuned during training. All selected hyper-parameters
including $B$ are described in Section 6.1.4 below.

##  6\. Model Tuning and Experimental Evaluation

In this section we describe training and validation of our model, and provide
an extensive empirical evaluation of our approach on test data.

###  6.1. Experimental Setting

In this section we discuss the different data-sets used to train and validate
our model, the training procedure, the metrics used to evaluate our model, and
the architecture and values of hyperparameters for our final model.

####  6.1.1. Data generation

We used our sampling procedure described in Section 4.1 above to produce three
different data sets: the training, validation, and test data sets. Their sizes
are provided in Table 2.

Data set | Training | Validation | Test
---|---|---|---
Size | 1,200,000 | 60,000 | 60,000
Table 2. Dataset sizes.

Training the deep learning model (namely, its weights) is performed on the
training set. The validation set is used during training for the purpose of
tuning hyper-parameters of the network. The model performance is assessed
based on the test data set.

####  6.1.2. Model Training and Tuning

The main logical flow of our process is as follows. Once all three data-sets
are sampled and pre-processed, including the computation of the output, we
train our model for every value of the number of input moments $n=2,..,20$.
For each $n$, we fine-tune our model, i.e., choose the best hyper-parameters
(e.g., number of hidden layers, learning-rate schedule, optimizer), based on
our validation set. Different hyper-parameters settings are compared via a
simple metric, $Metric1$, given by the following equation:

(2) $\displaystyle
Metric1(Y,\hat{Y})=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=0}^{l-1}|Y_{i,j}-\hat{Y}_{i,j}|,$

where $Y$ and $\hat{Y}$ correspond to two matrices of size $(N,l)$ of the true
values of the stationary queue-length distribution and of the network’s
predicted values, respectively. Here $N$ is the size of the entire validation
data set. Note that $Metric1$ is the absolute sum of errors (i.e., absolute
distances between the true and predictied distributions) computed for each
sample and then averaged over the entire validation data set. The main
advantages of this function is that it has a straightforward probabilistic
meaning and it covers the entire PDF of the stationary queue length. Also, for
each sample it provides an upper bound on the maximum difference between the
two CDF functions, which has statistical justification when examining
distances between two distributions, e.g., in the Kolmogorov–Smirnov test.

For each value of $n$ we choose settings that minimizes the value of
$Metric1$. Once the optimal settings for each $n$ are determined, we search
for the ‘best’ value of $n$, i.e., the value that achieved the minimal value
of $Metric1$ over the validation set. However, if we reach convergence, i.e.,
for different values of $n$ the value of $Metric1$ does not change by more
than $\epsilon=10^{-6}$, we select the smallest value of $n$ for which this
convergence is achieved. The intuitive idea is that if adding more moments
will give only a negligible added benefit value, then it is not worth doing.

####  6.1.3. Analyzing the number of service time moments

In Figure 4 we display our results of the service moment analysis. We compute
the value of $Metric1$ from Equation (2) as a function of the number of
moments. As shown in the figure, convergence is achieved once the number of
moments reaches 5. This implies that only the first 5 service time moments
appear to determine the dynamics of the $M/G/1$ queue. This result is somewhat
surprising as two distributions sharing the same first 5 moments may differ
significantly, yet will, apparently, result in the same stationary behavior.
On the other hand, as discussed earlier, the average queue length is
determined by only the first two service time moments, which makes this result
more sensible.

Figure 4. Service moments analysis.

####  6.1.4. Deep learning model architecture for $n=5$

Once the best value of $n$ has beem determined to be $5$, we fixed the values
of all other hyper-parameters. The resulting network is described below.

We have trained a fully connected neural network described in Section 5. In
total, there are 5 hidden layers with 14,399 parameters. We have used the Adam
optimization algorithm (Kingma14, ) for an iterative update of the network
weights based on the training data. The corresponding weight decay value for
avoiding over-fitting was $10^{-5}$. The selected training batch size was
$128$. The learning rate decayed exponentially, starting from 0.01. The number
of trained epochs was 300. The model was trained on an NVidia A100 GPU. For
further details about the chosen deep learning model see Appendix 6.

In terms of computational complexity this is considered to be a small network
(a large network typically contains several millions parameters666For example
the ResNet101, which is a commonly used architecture in computer vision has
44.5 million parameters.). This leaves a lot of room for adding network
complexity as machine learning methodology is extended to more complex
queueing systems.

###  6.2. Results

In this part, we present the evaluation of our machine learning system on the
test dataset. Throughout this section we use the “best” network with $n=5$,
tuned as described above.

For each sample in the test set we compare the true stationary queue length
distribution to its counterpart predicted by our model. In addition to
$Metric1$ (absolute distance) given by equation (2) above, where $N$ now
represents the size of the test data set, we also define the following
additional metric which focuses on the inverse of the CDF of the number of
customers in the system. Specifically, we compute the average relative error
for each of the $6$ different percentile values, namely,
$25\%,50\%,75\%,90\%,99\%,99.9\%$, as follows:

(3) $\displaystyle
Metric2(Y,\hat{Y},percentile)=\frac{1}{N}\sum_{i=1}^{N}\frac{F_{Y_{i}}^{-1}(percentile)-F_{\hat{{Y_{i}}}}^{-1}(percentile)}{F_{Y_{i}}^{-1}(percentile)}.$

This metric is particularly effective in detecting differences in the tails of
the two distribution functions; as noted earlier accurate estimates of tail
probabilities are particularly important in many practical applications of
queuing systems.

####  6.2.1. Performance analysis for the “best” model

In this section we demonstrate how well our deep learning model approximates
the stationary probabilities. The mean value of the error of our first metric
(i.e., Equation (2)) under $n=5$ is $0.0009$. In Figure 5 we present the error
histogram of $Metric1$, to provide a better picture of our performance. In
Table 3 we present the performance of our model under both metrics.

Note that the interpretation of values for $Metric1$ and $Metric2$ are quite
different. For $Metric1$ (the seond row of the table) the values represent the
exact percentiles of the absolute difference between true and predicted
distributions over the test data set. For example, the $99.9$-th percentile
was only $.02$, indicating an excellent agreement between the two
distributions. For $Metric2$ (last row of the table) the values represent
relative absolute differences between the true and predicted distribution at
each percentile level. For example, the largest relative difference was only
$0.412\%$ achieved at the $99.9$-th percentile.

Percentile | 25% | 50% | 75% | 90% | 99% | 99.9 %
---|---|---|---|---|---|---
Metric 1 | 0.000220 | 0.0004272 | 0.00087 | 0.002001 | 0.008757 | 0.02001
Metric 2 | 0.000342 | 0.000718 | 0.00094 | 0.000897 | 0.00148 | 0.00412
Table 3. Deep learning model performance on Test Data. Figure 5. An histogram
of the error metric presented in Equation (2).

As demonstrated in Figure 5 and in Table 3 the errors are very small, where in
99% of the cases the error is smaller than 0.01. This implies that our
predictions are extremely reliable. Furthermore, as the result of $Metric2$
suggest, the predicted distribution produced by our model remains very
accurate in the tail as well, where the average error for all percentiles is
very low.

##  7\. Case Study: Deep Learning vs. Analytical Queueing Methods

The results presented in the preceding section show that our deep learning
network is able to achieve excellent predictive accuracy on the test data set.
However, the experimental conditions used in the current section could be
criticized as not very reflective of “real-life” settings: after all, for each
test sample the network was supplied with the exact values of the first 5
moments of the service distribution, as well as the exact arrival rate. In
more realistic settings, only an empirical set of data observations for
service times and arrivals is likely to be available. How can our network be
expected to perform under these conditions? We aim to provide an answer in the
current section.

We shall assume that instead of having a closed-form expressions for the
inter-arrival and service time distributions, we are only provided with the
timestamps of arrivals, service commencements, and service completions -
essentially the system event log. This is the case in most real-life
situations.

From the managerial point of view, given past data one may wish to ask “what-
if?” questions. For example, what if the arrival rate will change next month,
while the service process remains the same? How will it effect system waiting
times? In this scenario, we would wish to analyze an $M/G/1$ queue where the
arrival rate is given, but the service times are known only via the data
sample.

Given such data, and assuming the system is known to operate as an $M/G/1$
queue, there are two approaches for obtaining the stationary distribution: the
”analytic approach” and the ”deep learning” one. Under the analytical
approach, we first use the data to fit some parametric distribution to service
times (we will assume a PH distribution is used), and then we apply one of the
available analytical expressions to obtain the stationary distributions (we
will assume the QBD method is used). The machine learning approach (assuming a
trained network is available) is simpler: we only need to compute the first
five service time moments from the data, and then we apply the deep learning
network to obtain the predicted stationary distribution.

Both methods will incur some errors. In the analytic approach the main
accuracy loss can be expected when fitting the service time distribution;
numerical errors from the application of the QBD calculation should be very
small (and controllable). In the deep learning method there are two potential
sources of error: due to the estimation of the service moments from the data,
and due to the prediction error of the deep learning model.

The case-study was executed as follows:

  1. (1)

We sample a service distribution using our Algorithm 1 from Section 4.1, scale
the mean service time to $1$, set the service rate to $0.85$ and then compute
the corresponding stationary queue length distribution. This gives us the
“Ground Truth”.

  2. (2)

Next, we sample 50,000 data points from the service time distribution
producing our “observed” service time sample.

  3. (3)

We apply the analytic approach to the data. This consists of using the EM
algorithm to fit the distribution (details are provided below), followed by
the QBD method. This gives us the “Analytic” estimate.

  4. (4)

We apply the deep learning approach producing the “Deep Learning” estimate.

Note that while the arrival rate is fixed at $0.85$ (even if arrivals were
provided as observational data, the estimation of the arrival rate is
straightforward by using the inverse of the average), fitting a service time
distribution is not easy. The true and empirical distributions of service
times along with a histogram of the data are presented in Figure 7.

Figure 6. Service data sample

Figure 7. PDF of number of customers in the system.

The best practice is to use the EM algorithm to fit the closest PH
distribution (see (Fackrell098, )). The weakness of the approach is its
computational complexity. The runtime of the EM algorithm depends on the size
of the PH, and thus we get a clear trade-off between accuracy and runtime. In
our experience, even for PH of size 20 (this refers to the dimension of the
transition matrix) it takes over an hour to run the EM algorithm in most
cases777This depends of course on the number of samples in the data set..
Moreover, the runtime grows exponentially with the size of the PH. In our
case-study we used maximum PH size of 20.

For the deep learning approach the first five moments were fitted by employing
an unbiased estimator directly from the data. That is,
$\bar{X^{i}}=\sum_{i=1}^{N}\frac{X^{i}}{N}$, where $X$ is the service time,
$N$ is the sample size, and $i=1,2,3,4,5$.

In Table 4 we present the first $5$ moments of service times estimated by each
method (for the analytic approach these were computed from the PH distribution
obtained by the EM method). We see that (not surprisingly) the direct
estimates used by the Deep Learning approach are far more accurate,
particularly with respect to higher-order moments. This points to substantial
differences between the service time distribution estimated by the EM method
and the ground truth distribution from which the data was generated. This can
also be observed from the comparison of actual and fitted distributions on
Figure 7.

Table 4 also presents pre-processing running times (for the analytic approach
this is the running time of the EM method, for the deep learning time this is
the time to compute the $5$ moments) and the algorithm running times (QBD
computation times for the analytic and ground truth approaches; network
computation time for the deep learning approach). The largest difference is in
the very long pre-processing time for the analytic approach compared to near-
instantaneous time for the deep learning one888The runtimes presented in Table
4 corresponding to the ‘EM’ procedure computed on the same server detailed in
Remark 3, i.e., Intel Xeon Gold 5115 Tray Processor, with RAM memory of 128
GB. .

| Analytic (EM) | Deep Learning | Ground Truth
---|---|---|---
First moment | 0.994 | 0.994 | 1
Second moment | 2.705 | 2.398 | 2.418
Third moment | 10.093 | 7.037 | 7.106
Fourth momnet | 44.131 | 21.961 | 22.19
Fifth moment | 217.936 | 71.227 | 71.97
Pre-procsseing time (sec) | 4320 | 0.3 | 0
Running time (sec) | 0.16 | 0.12 | 0.44
Table 4. Service moments and runtimes

Next, we compare the accuracy of the two approaches. The histogram of
stationary queue-length distributions under the three approaches is presented
on Figure 7 and the values of $Metric2$ for different percentiles are shown on
Table 5. We also computed $Metric1$ for both techniques.

While both the analytical and the deep learning approach do a fairly good job
of estimating the stationary queue length distribution, the deep learning
approach is the clear winner. The absolute difference between the estimated
and the true distributions according to Equation (2) (i.e., $Metric1$) for the
analytic method is 0.094, while it is only 0.002 for the deep learning method.
Figure 7 and Table 5 show that the difference in performance is mostly due to
the worse fit in the tail of the distribution for the analytic method. While
the deep learning model incurs minimal errors at every percentile on Table 5,
the errors under the analytic approach are quite substantial starting with
$75$-th percentile.

To conclude, this case-study demonstrates the practical capability of our
model, especially in answering “what-if?” questions. Whereas, moments are easy
to estimate from the data, fitting a service time distribution is both
computationally burdensome, and can lead to significant accuracy loss.

##  8\. Discussion & Limitations

In this section, we discuss some of the main limitations of our approach,
propose several solutions to overcome them, and compare our deep learning
model to employing simulation in queueing systems for the task of performance
prediction.

###  8.1. Limitations of our Approach

The main challenge is extending the deep learning approach to more complex and
realistic queuing systems for which, unlike the $M/G/1$ case analyzed above,
no analytical solutions are available. The issue is obtaining large and
reliable data sets on which the deep learning network can be trained. We
distinguish between two types of data sets.

##### Real data.

Real data that describes the queueing dynamics of the system can be extracted
from databases of many organizations. The advantage of the real data is that
“ground truth” of actual system performance is directly available. The main
problem is that one cannot typically generalize it to other systems, or even
to different time periods for the same system, as the data essentially
captures a single realization for a single queuing system. This issue can
potentially be addressed by combining different data sets together, and by
applying principles from transfer learning. Transfer learning is a machine
learning method where a model developed for a task is reused as the starting
point for a model on a second task . This is common practice in deep learning,
and it is empirically proven to be efficient (see
(10.1007/978-3-030-01424-7_27, )).

Percentile | Analytic | Deep learning
---|---|---
25 | 0.5 | 0.5
50 | 0 | 0
75 | 0.222 | 0
90 | 0.125 | 0
99 | 0.187 | 0
99.9 | 0.163 | 0.020
Table 5.  Comparison of $Metric2$

##### Synthetic queueing data.

Such data is sampled from a queueing system with a known structure, as we have
done in this paper. However, producing training output is not straightforward
since we are typically interested in applying machine-learning models to
systems for which closed-form solutions do not exist. This situation imposes a
catch: we can use our method to learn “queueing theory”, but only for queueing
models that the theory already knows how to solve.

We discuss several approaches for dealing with this issue. The most direct
extension is using our PH sampling algorithm to produce synthetic data for
systems where arrival and/or service processes are of PH form (i.e., $PH/M/1$,
$PH/PH/1$ systems) - such systems often have analytical solutions. As our
results indicate, PH-type distribution is a good proxy for the general case,
thus extensions to, for example, $G/G/1$ queue should be relatively easy.

For more complex systems, we may be able to perform transfer learning from
tractable queueing systems. An example for such an approach would be the task
of deriving stationary queue-length distribution of a $G/G/c$ queue. We can
use data from $M/M/c$ and $PH/PH/1$ queues, which both have analytical
solutions.

Another option is to use simulation models to generate training samples. In
principle, a simulation model can deal with very high level of complexity of
the underlying system and produce stationary distributions. The practical
limitation is that the runtime until convergence tends to be long, even for
relatively simple systems (particularly in the multi-server case). Thus, while
we were able to generate over a million training samples in a matter of hours
using our sampling algorithm, this may not be possible with simulation models.

While constructing a large training data set with simulation alone is not
practical, it may be possible to train a simpler model (e.g., $G/G/1$ with
distinct customer classes) and then (again) apply ’transfer learning’ with
fewer data points on a more complex model (e.g., the $G/G/c$ counterpart).

###  8.2. Discussion: Machine Learning vs. Simulation

Simulation is often the tool of choice when the underlying system is
intractable, i.e., a closed-form solution of the steady-state distribution is
unavailable. If the model is correct, we are guaranteed that if we simulate
for long enough our results will be very accurate. However, there are several
advantages we would like to note when using our approach compared to
simulation. The first is runtime: as noted earlier, it takes a long time for
simulations to converge, while a deep learning model provide predictions
immediately after it is trained. The second benefit of our approach lies in
the deep learning model itself. Although these models are often considered
black-boxes system abstractions, similarly to simulation, they may provide
insights beyond what can be achieved with simulation models. For example, the
moments analysis we performed done in Section 6.1.3 would be impossible when
using simulation, as the latter require a full specification of the service
time distribution.

More generally, the ability of the deep learning model to work with limited
data (first $5$ moments in our case) rather than a full specification of the
distribution gives it potentially great advantages over the simulation
approach.

###  8.3. Deep Learning for Fitting Distributions to Data

Before leaving this section we should mention one more interesting application
of our methodology. Fitting a distribution to data is of general interest in
many contexts, both within and beyond queuing applications. Although the PH
distribution family is very versatile, fitting it to data is not a simple
task. As we saw in Section 7, the standard approach of using EM algorithm is
neither accurate, nor computationally efficient.

In our model we used moments to represent a given PH distribution. However, we
could replace the roles, using moments computed from the data as inputs, while
the PH representation (i.e., the initial state distribution vector and the
transition matrix) would become the output. Our sampling algorithm can easily
generate training samples, and - as we saw in Section 7 \- computing moments
directly from the observational data is both accurate and efficient.

In fact, there are several studies that deal with fitting PH given the first
$n$ moments. For example in (doi:10.1081/STM-200056210, ), the authors present
an explicit method to compose minimal order continuous-time acyclic phase type
(APH) distributions given only the first three moments. Yet their method does
not always find a solution and it is limited to the first three moments. We
feel our deep learning approach can be re-purposed to provide more efficient
solutions to this important problem.

##  9\. Conclusion

In this paper we empirically demonstrated the ability of machine to analyze a
general queueing model, namely the stationary queue-length distribution of an
$M/G/1$ queue. Our results show that machines can accurately learn these
distributions using only a relatively compact deep learning architecture. This
gives hope for future studies where more complex queueing models (and thus
more complex deep learning models) will be incorporated.

The unique method of evaluating the stationary queue length distribution using
only service time moments instead of the full analytical representation,
allows us to (empirically) evaluate the number of service time moments
required for determining the stationary queue-length distribution.
Surprisingly, the accuracy does not significantly improve when we go beyond
the $5^{th}$ moment. We also provided a case study that describes a real-life
setting when data samples are available. We demonstrated that the deep
learning model can outperform analytical (data-driven) methods both in terms
of accuracy and runtime.

Finally, we believe that this approach will be very useful in future studies,
as our deep learning approach excelled in capturing the stochastic dynamics of
a stationary queueing model. As discussed in the preceding section, extensions
to $G/M/1$, $G/G/1$ and similar models should be rather direct. Applications
to more complex settings, as well as to non-queuing settings such as fitting
PH distribution to data, are also promising.

## References

  * [1] Abien Fred Agarap.  Deep learning using rectified linear units (relu).  arXiv preprint arXiv:1803.08375, abs/1803.08375, 2018.
  * [2] Asmussen, Soren.  Applied Probability and Queues.  Springer New York, New York, NY, 2003.
  * [3] A. Bobbio, A. Horváth, and M. Telek.  Matching three moments with minimal acyclic phase type distributions.  Stochastic Models, 21(2-3):303–326, 2005.
  * [4] Mu-Song Chen and Hao-Wei Yen.  Applications of machine learning approach on multi-queue message scheduling.  Expert Systems with Applications, 38(4):3323–3335, 2011.
  * [5] D. R. Cox.  The analysis of non-markovian stochastic processes by the inclusion of supplementary variables.  Mathematical Proceedings of the Cambridge Philosophical Society, 51(3):433–441, 1955.
  * [6] J. G. Dai and Mark Gluzman.  Queueing network controls via deep reinforcement learning, 2021.
  * [7] Dmitry Efrosinin and Natalia Stepanova.  Estimation of the optimal threshold policy in a queue with heterogeneous servers using a heuristic solution and artificial neural networks.  Mathematics, 9(11), 2021.
  * [8] Mark Fackrell.  Modelling healthcare systems with phase-type distributions.  Health Care Manag Sci, (12):11–26, 2009.
  * [9] Harchol-Balter, Mor.  Performance Modeling and Design of Computer Systems: Queueing Theory in Action.  Cambridge University Press, 2013.
  * [10] Olawoyin R. Hijry H.  Predicting patient waiting time in the queue system using deep learning algorithms in the emergency room.  International Journal of Industrial Engineering and Operations Management, 3(1):33–45, 2021.
  * [11] Leonard. Kleinrock.  Queueing systems, volume 2: Computer applications. by leonard kleinrock john wiley & sons, inc., new york 1976, 549 pages,.  Networks, 7(3):285–286, 1977.
  * [12] Athanasios I. Kyritsis and Michel Deriaz.  A machine learning approach to waiting time prediction in queueing scenarios.  In 2019 Second International Conference on Artificial Intelligence for Industries (AI4I), pages 17–21, 2019.
  * [13] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [14] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521(7553):436–444, May 2015.
  * [15] Yann LeCun, Koray Kavukcuoglu, and Clement Farabet.  Convolutional networks and applications in vision.  In Proceedings of 2010 IEEE International Symposium on Circuits and Systems, pages 253–256, 2010.
  * [16] Le Liang, Hao Ye, Guanding Yu, and Geoffrey Ye Li.  Deep-learning-based wireless resource allocation with application to vehicular networks.  Proceedings of the IEEE, 108(2):341–356, 2020.
  * [17] Bai Liu, Qiaomin Xie, and Eytan Modiano.  Reinforcement learning for optimal control of queueing systems.  In 2019 57th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 663–670, 2019.
  * [18] Yuta Makino, Tuan Phung-Duc, and Fumio Machida.  A queueing analysis of multi-model multi-input machine learning systems.  In 2021 51st Annual IEEE/IFIP International Conference on Dependable Systems and Networks Workshops (DSN-W), pages 141–149, 2021.
  * [19] Marcel F. Neuts.  Phase-type Probability Distributions, pages 1132–1134.  Springer US, Boston, MA, 2013.
  * [20] Neuts, Marcel F.  Matrix-Geometric Solutions in Stochastic Models. The Johns Hopkins University Press, 1981.
  * [21] Suguru Nii, Takashi Okudal, and Takuya Wakita.  A performance evaluation of queueing systems by machine learning.  In 2020 IEEE International Conference on Consumer Electronics - Taiwan (ICCE-Taiwan), pages 1–2, 2020.
  * [22] Kingma D. P. and Ba J.  Adam: A method for stochastic optimization.  ArXiv preprint:1412.6980, abs/1803.08375, 2014.
  * [23] Chuanqi Tan, Fuchun Sun, Tao Kong, Wenchang Zhang, Chao Yang, and Chunfang Liu.  A survey on deep transfer learning.  In Věra Kůrková, Yannis Manolopoulos, Barbara Hammer, Lazaros Iliadis, and Ilias Maglogiannis, editors, Artificial Neural Networks and Machine Learning – ICANN 2018, pages 270–279, Cham, 2018. Springer International Publishing.
  * [24] M. Telek and G. Horváth.  A minimal representation of markov arrival processes and a moments matching method.  Performance Evaluation, 64(9):1153–1168, 2007.  Performance 2007.
  * [25] WARD WHITT.  Approximations for the gi/g/m queue.  Production and Operations Management, 2(2):114–161, 1993.
  * [26] Ward Whitt.  Towards better multi-class parametric-decomposition approximations for open queueing networks.  Annals of Operations Research, 48(3):221–248, Jun 1994.
  * [27] Wencong Xiao, Romil Bhardwaj, Ramachandran Ramjee, Muthian Sivathanu, Nipun Kwatra, Zhenhua Han, Pratyush Patel, Xuan Peng, Hanyu Zhao, Quanlu Zhang, Fan Yang, and Lidong Zhou.  Gandiva: Introspective cluster scheduling for deep learning.  In 13th USENIX Symposium on Operating Systems Design and Implementation (OSDI 18), pages 595–610, Carlsbad, CA, October 2018. USENIX Association.

##  Appendix A Deep Learning model

In this section we provide additional details on the chosen deep learning
model. In Figure 8, we present the training loss and the value of $Metric1$ on
the validation set as a function of the training epoch. In Table 6 we specify
the network specific architecture.

(a)

(b)

Figure 8. Training Phase Hidden layer | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
Nodes | 30 | 40 | 50 | 60 | 60
Table 6. Network nodes.