# Modelling multi-cell edge video analytics ††thanks: This work has been
submitted to IEEE for publication. Copyright may be transferred without
notice, after which this version may no longer be accessible.

Jaume Anguera Peris and Viktoria Fodor  School of Electrical Engineering and
Computer Science
KTH Royal Institute of Technology, Stockholm, Sweden
Email<EMAIL_ADDRESS>
###### Abstract

Edge intelligence is a scalable solution for analyzing distributed data, but
it cannot provide reliable services in large-scale cellular networks unless
the inherent aspects of fading and interference are also taken into
consideration. In this paper, we present the first mathematical framework for
modelling edge video analytics in multi-cell cellular systems. We derive the
expressions for the coverage probability, the ergodic capacity, the
probability of successfully completing the video analytics within a target
delay requirement, and the effective frame rate. We also analyze the effect of
the system parameters on the accuracy of the detection algorithm, the
supported frame rate at the edge server, and the system fairness.

###### Index Terms:

Edge intelligence, stochastic geometry, queuing theory, edge computing

##  I Introduction

The evolution of mobile devices has introduced tremendous innovation
opportunities for developing new services beyond personal communications.
These advances have brought new mobile applications, and with them, the
proliferation of large amounts of app-level data. At the same time, we have
experienced major breakthroughs in deep learning, which have made it possible
to manage and leverage the rise in data volumes for solving complex problems,
such as image classification, language processing, or intrusion detection [1].

However, as we strive for more accurate and reliable deep learning models, the
computationally and energy-constrained mobile devices make real-time
processing intractable. One prominent solution to extend the capabilities of
mobile devices is to offload the computationally intensive tasks to an edge
computing server. Edge computing offers a distributed solution that brings the
computational resources closer to the mobile devices, and provides lower
latency, more energy efficiency, better privacy protection, and lesser
bandwidth consumption than its counterpart cloud-based infrastructure [2].
Edge computing also makes it easier to gather richer data and learn from the
behaviour of nearby mobile devices, making it possible to explore deep
learning models for interactive cloud gaming, cognitive assistance, and real-
time video analytics [2].

Indeed, the growing interest in both edge computing and artificial
intelligence has coined the term edge intelligence. This area of research is
however in its early stages, and most of the proposed solutions are optimized
to be deployed on local networks, such that it can be assumed that all users
have similar network conditions and similar probabilities of successfully
completing their task. In large-scale cellular networks, aspects of noise,
fading, or interference can significantly impact the wireless communication
[3], so those assumptions no longer hold, and analyzing the communication
between the mobile users and the edge server becomes as relevant as analyzing
the implementation of the deep learning models. Hence, there needs to be a
tractable approach to understand how the network, the computing resources, and
the deep learning models affect altogether the task offloading problem.

In this work, we present a mathematical framework for analyzing edge-
intelligent, multi-cell cellular systems supporting video analytics. We
specifically leverage the tools of stochastic geometry, queuing theory, and
deep learning for modelling the entire offloading process under different
accuracy-latency constraints. We provide a thorough analysis of the
relationship between the parameters of the video analytics, the uplink
transmission, the edge server, and the deep learning algorithm, and evaluate
their trade-offs via numerical results. Overall, our discussions and results
provide the first steps towards designing efficient resource-allocation
strategies and traffic-control protocols for edge-intelligent, multi-cell
cellular systems.

The remainder of this paper is organized as follows. Section II reviews the
related work. Section III describes the system model for the entire offloading
process. In Section IV, we derive the expressions for the performance metrics
that evaluate the trade-offs between the different parameters in the system.
Then, we provide the numerical results in Section V, followed by the
concluding remarks in Section VI.

##  II Related work

Edge intelligence is a promising solution to assist mobile devices with
computationally-intensive and latency-stringent applications, and it is
recently receiving a lot of attention for visual-aid services, e.g., augmented
reality [4] and video edge analytics [5]. As a response to these emerging
interests, recent works have built proof-of-concept testbeds [4, 5, 6]. DARE
[4] focuses on the edge server and analyzes the effect of the frame size, the
frame rate, and the deep learning model on the task offloading problem.
DeepDecision [5] considers the same parameters but instead focuses on the
mobile users and the wireless network. The work in [6] uses its own
measurements to design data-driven statistical models for real-time inference.

The results from the above works highlight the suitability of edge
intelligence for improving the accuracy and latency of visual-aid services,
but their simplified models for the communication between the mobile users and
the edge servers do not capture the inherent aspects of large-scale cellular
networks. In this context, [7] models the performance of edge-intelligent
systems with parallel processing for video analytics, and [8] models the
performance of edge computing in multi-cell networks, but does not address the
performance of the application itself.

Our work extends the above results and presents a mathematical framework for
analyzing how the conditions of the network, the frame size, and the latency
requirements affect the transmission rate, the achievable frame rate, the
probability of successfully completing the video analytics, and the accuracy
of the detection algorithm.

##  III System model

###  III-A Network

Consider a multi-cell wireless network where all the base stations employ
frequency-division multiple access with reuse factor $\delta$, such that
mobile users experience inter-cell interference but no intra-cell
interference. Each mobile user connects to the base station that provides the
maximum received power averaged over fading, and each base station has only
one active user in a given frequency band, which is arbitrarily chosen from
all the mobile users located in its Voronoi cell. Motivated by [9], the
locations of the base stations are modelled according to a two-dimensional
homogeneous Poisson point process (PPP) with intensity $\lambda_{b}$.
Similarly, the locations of the mobile users using the same frequency
resources are considered to be independent and modelled according to some
other stationary PPP with intensity $\lambda_{u}$. Considering this, it can be
shown that the distance between any user and its serving base station, denoted
by $r$, is Rayleigh distributed. The proof follows from the null probability
of a two-dimensional PPP [10],

$\mathcal{P}[R>r]=\mathcal{P}\left[\text{no BS in a circle of area }\pi
r^{2}\right]=e^{-\lambda_{b}\pi r^{2}},$

from which the probability density function (PDF) is

$f_{R}(r)=2\lambda_{b}\pi r\,e^{-\lambda_{b}\pi r^{2}},\quad r\geq 0.$ (1)

For the rest of the paper, we consider a target user $k$ and its serving base
station, and refer to the set of interfering mobile users as $\mathcal{Z}$. We
denote the distances from the interfering users $z\in\mathcal{Z}$ to the base
station of interest by $d_{z}$, and the distances from the interfering users
to their serving base stations by $r_{z}$. Note that the distances
$\\{r_{z}\\}_{z\in\mathcal{Z}}$ are identically distributed but not
necessarily independent. The dependence comes from the structure of the
Poisson-Voronoi tessellation and the restriction that there is only one active
user in each frequency band. However, [11] shows that this dependence is weak
and the distances $\\{r_{z}\\}_{z\in\mathcal{Z}}$ can be assumed to be
independent and identically distributed (i.i.d.). Furthermore, the marginal
distribution of $r_{z}$ can be assumed to be Rayleigh distributed as in (1),
which is motivated by the irregular deployment of base stations and the same
null-probability argument as for $r$.

###  III-B Uplink transmission

In edge video analytics, the amount of data transmitted in the uplink, i.e.,
from the users to the base stations, is much larger than the amount of data
transmitted in the downlink. Therefore, we focus on modelling the uplink
transmission.

We consider that the transmitted signals attenuate with path-loss exponent
$\alpha>2$ and experience i.i.d. Rayleigh fading. All the users utilize a
maximum-ratio transmission precoding vector, such that the effect of the
channel can be modelled at the receiver by the multiplicative Exponential
random variable $|g|^{2}\sim\exp(1)$. We define the largest possible time-
frequency space within which the effect of the channel can be modelled as a
scalar multiplication by $|g|^{2}$ as the coherence interval. We denote the
coherence interval by $\tau_{c}=T_{coh}\,B_{coh}$, where $T_{coh}$ represents
the coherence time, and $B_{coh}$ represents coherence bandwidth. Considering
this, we define the number of coherence intervals that fit into a time-
frequency block of $t$ seconds and $B$ Hertz as
$N_{\textit{coh}}=t\,B/\tau_{c}$.

Moreover, we consider that all mobile users utilize distance-proportional
fractional power control of the form $r^{\alpha\epsilon}$, where
$\epsilon\in[0,1]$ is the power control factor. We define a reference power at
a distance of $1$ kilometer from the base station as $P$ and assume it to be
the same for all the users. To model the maximum transmission power of the
mobile users, we consider the transmit powers $Pr^{\alpha\epsilon}$ to be
limited by the peak power $\bar{P}$. Finally, we consider the noise power at
the receiver to be $\sigma^{2}$. Under this transmission model, the signal-to-
interference-noise ratio (SINR) for any randomly selected user $k$ at its
serving base station is

$\text{SINR}_{k}=\frac{|g_{k}|^{2}\,\ell(r,\alpha,\epsilon)\,r^{-\alpha}}{\sigma^{2}_{k}+\sum_{z\in\mathcal{Z}}|g_{z}|^{2}\,\ell(r_{z},\alpha,\epsilon)\,d_{z}^{-\alpha}},$
(2)

where

$\ell(r,\alpha,\epsilon)=\min\left(Pr^{\alpha\epsilon},\bar{P}\right),$ (3)

and the sub-indices $k$ and $z$ indicate if a variable depends on the location
and the resources assigned to the target user $k$ or to the interfering users
$z\in\mathcal{Z}$, respectively.

For the uplink transmission, we consider that each base station assigns a
fixed bandwidth $B_{k}$ to every active user within its Voronoi cell. Once a
user $k$ has been assigned a frequency band, it estimates the state of the
channel in each coherence interval from pilot signals and adjusts its coding
and modulation according to its instantaneous SINRk. As a result, the
transmission rate $R_{k}$ is a random variable that depends on the user’s
location, the conditions of the channel, and the inter-cell interference.

###  III-C Mobile edge computing server

All the base stations in the network are equipped with an independent mobile
edge computing server that is responsible for serving all the video analytic
users within the Voronoi cell of its corresponding base station. For our base
station of interest, we consider that there is a set of $K$ active video
analytic users, denoted by $k$, for $k=1,\dots,K$, that are capable of
continuously recording images and sending them over the wireless channel. Once
the images of these users are completely received at the base station, they
are delivered to the server without delay and queued at an infinite buffer for
being processed on a first-come-first-served basis.

From the point of view of the server, each user $k$ is seen as an independent
entity that offloads images to the server at an arrival rate of $\lambda_{k}$
frames per second. Altogether, since the communication system employs
frequency-division multiple access, the queue at the server observes the
superposition of $K$ independent arrival processes with an aggregate arrival
rate of $\lambda=\sum_{k}\lambda_{k}$ frames per second. Under these
assumptions, the arrival process behaves as a Poisson process [12]. Once in
queue, the images are processed at the server in a deterministic amount of
time $T_{s}$, which is determined by the resolution of the transmitted images
and the neural network considered for the video analytics.

In this paper, we consider for simplicity that all the users within the same
cell transmit images of the same resolution and perform the video analytics
utilizing the same neural network. Hence, it yields that the queuing system at
the server can be modelled as an M/D/1 queue with load $\rho=\lambda T_{s}$.

###  III-D Video analytics

All the servers in the network incorporate an object detection algorithm
consisting of a pre-trained YOLOv5 neural network. We specifically consider
YOLOv5 [13] for the video analytics because it is capable of handling images
of different resolutions without changing its associated learnable parameters,
thus making it possible to parametrize the inference time and the accuracy of
the detection algorithm.

Following the results in [14], the number of instructions that are required to
process an image with a resolution of $s_{k}^{2}$ pixels can be expressed as a
convex function of $s_{k}$,

$f(s_{k})=c_{1}s_{k}^{3}+c_{2},$ (4)

for some positive constants $c_{1}$ and $c_{2}$, where the units of $f(s_{k})$
are number of floating-point operations (in trillions). Similarly, the
accuracy of the detection algorithm can be expressed as a concave function of
$s_{k}$,

$A(s_{k})=c_{3}-c_{4}e^{-c_{5}s_{k}},$ (5)

for some other positive constants $c_{3}$, $c_{4}$, and $c_{5}$, where the
accuracy of the detection is measured as the mean average precision of the
object detection algorithm for a predefined threshold of the intersection over
union [14]. For more information about the effect of the learning and
inference processes on the constants $c_{1},\dots,c_{5}$, please refer to
[15].

Moreover, all the computing servers in the network are equipped with a CPU of
$h$ TFLOPs. With that, it follows from (4) that the time to process an image
of resolution $s_{k}^{2}$ is

$T_{s,k}=\frac{f(s_{k})}{h}=\frac{c_{1}s_{k}^{3}+c_{2}}{h}.$ (6)

Note that images of high resolution take longer to process but they result in
higher detection accuracies, with a nonlinear trade-off between the processing
time and the accuracy of the detection algorithm.

###  III-E Offloading process

Altogether, the offloading process from beginning to end can be described as
follows. All active video analytics users are assigned some frequency
resources and are considered to transmit with high priority, that is, they
always get access to the wireless channel as soon as they want to offload
their images. The transmitted images are considered to be squared with a
resolution of $s_{k}^{2}$ pixels, encoded at $\theta$ bits per pixel and
compressed at rate $\xi$:1, so an image of $s_{k}^{2}$ pixels is represented
by $\theta s_{k}^{2}/\xi$ bits. Once the images are completely received at the
serving base station, the images are queued at the edge server and processed
one by one by a convolutional neural network. For each user $k$, we denote the
time to transmit an image as $T_{ul,k}$, the time spent in the queue as
$T_{w,k}$, and the time to process an image as $T_{s,k}$. We refer to the sum
of these three times as the delay of the entire offloading process,
disregarding the little time of eventual feedback to the user, and let each
user $k$ specify a maximum value $t_{k}$ for that delay. For simplicity, we
assume that the maximum delay $t_{k}$ is the same for all the users within the
same cell.

##  IV System analysis

###  IV-A Coverage probability

The coverage probability for any user $k$ in the network is the probability
that its uplink SINRk is greater than some threshold $\gamma_{k}$. Considering
this, the coverage probability conditioned to user $k$ being at a fixed
distance $r$ from its serving base station can be derived from [11, Theorem
1], for any reuse factor $\delta$ [16, Theorem 4], and power-limit constraint
(3) as

$\displaystyle\mathcal{P}(\text{SINR}_{k}\geq\gamma_{k}|r)$
$\displaystyle\;=\mathcal{P}\left(|g_{k}|^{2}\geq\frac{\gamma_{k}\,(\sigma_{k}^{2}+\sum_{z\in\mathcal{Z}}|g_{z}|^{2}\,\ell(r_{z},\alpha,\epsilon)\,d_{z}^{-\alpha})}{\ell(r,\alpha,\epsilon)\,r^{-\alpha}}\bigg{|}r\right)$
$\displaystyle\;\overset{(a)}{=}e^{-\frac{\gamma_{k}\,r^{\alpha}}{\ell(r,\alpha,\epsilon)}\,\sigma^{2}_{k}}\;\mathbb{E}\left[\exp\left(-\frac{\gamma_{k}\,\sum_{z\in\mathcal{Z}}|g_{z}|^{2}\,\ell(r_{z},\alpha,\epsilon)\,d_{z}^{-\alpha}}{\ell(r,\alpha,\epsilon)\,r^{-\alpha}}\right)\right]$
$\displaystyle\;\overset{(b)}{=}e^{-\frac{\gamma_{k}\,r^{\alpha}}{\ell(r,\alpha,\epsilon)}\,\sigma^{2}_{k}}\mathcal{L}_{I_{z}}\left(\frac{\gamma_{k}\,r^{\alpha}}{\ell(r,\alpha,\epsilon)}\right),$
(7)

where $(a)$ follows from the exponential distribution of $|g_{k}|^{2}$, and
$(b)$ follows from the definition of the Laplace transform of the interference
at the serving base station. This Laplace transform can be calculated from the
i.i.d. distributions of $r_{z}$ and $d_{z}$, and the Probability Generating
Function [10] of the PPP that models the location of the mobile users, such
that

$\mathcal{L}_{I_{z}}(s)=\exp\left(-2\pi\lambda_{u}\int_{r}^{\infty}\beta(x,s)\,x\,dx\right),$
(8)

with

$\beta(x,s)=1-\int_{0}^{\infty}\frac{\pi\lambda_{u}\,e^{-\pi\lambda_{u}u}}{1+s\,\ell\left(u^{\frac{1}{2}},\alpha,\epsilon\right)x^{-\alpha}}\,du,$

where we considered that the closest interferer to the serving base station of
user $k$ is at least at distance $d_{z}\geq r$.

Note that the density of users using the same frequency resources always
satisfies $\lambda_{u}\leq\lambda_{b}/\delta$ because there cannot be more
active users than available frequency resources at the base stations.

###  IV-B Ergodic capacity

From the estimates of the state of the channel, any active user $k$ can adjust
its coding and modulation according to its instantaneous SINRk to achieve a
transmission rate close to the Shannon capacity,
$B_{k}\log_{2}(1+\text{SINR}_{k})$. Considering this, and following the same
steps as the ones considered to find (7), we can calculate the ergodic
capacity for any user $k$ at distance $r$ from its serving base station as
[16, Theorem 3]

$\displaystyle\bar{R}_{k}|r$
$\displaystyle=\mathbb{E}[B_{k}\log_{2}(1+\text{SINR}_{k})|r]$
$\displaystyle=\frac{B_{k}}{\log(2)}\int_{0}^{\infty}\mathcal{P}(\text{SINR}_{k}\geq
e^{t}-1|r)dt$
$\displaystyle=\frac{B_{k}}{\log(2)}\int_{0}^{\infty}e^{-\frac{(e^{t}-1)\,r^{\alpha}}{\ell(r,\alpha,\epsilon)}\sigma^{2}_{k}}\mathcal{L}_{I_{z}}\left(\frac{(e^{t}-1)\,r^{\alpha}}{\ell(r,\alpha,\epsilon)}\right)dt,$
(9)

where the second equality follows from the fact that the expectation of any
positive random variable $X$ can be calculated as
$\mathbb{E}[X]=\int_{0}^{\infty}\mathcal{P}\left(X\geq t\right)dt$, and the
last equality follows from the result in (7) by substituting $\gamma_{k}$ for
$e^{t}-1$.

###  IV-C Probability of successful computation

The probability of successful computation measures the probability that the
entire offloading process is completed within the maximum delay requirement
$t_{k}$. We consider that the offloading process starts when a user begins
offloading an image and ends when the server finishes processing that image.
Hence, for a user $k$ at distance $r$ from its serving base station and for a
maximum delay requirement $t_{k}$, the probability of successful computation
can be expressed as

$\displaystyle p_{succ}$
$\displaystyle\left(t_{k},s_{k},\lambda,\lambda_{b},\lambda_{u},\delta,P,\alpha,\epsilon,\bar{P}\,|\,r\right)=$
$\displaystyle\mathcal{P}\left(T_{ul,k}+T_{w,k}+T_{s,k}\leq
t_{k}\,|\,r\right).$ (10)

Note that $T_{ul,k}$ is a random variable that depends on the sources of
randomness in SINRk, $T_{w,k}$ is a random variable that depends on the number
of images that are already waiting in queue, and $T_{s,k}$ is a deterministic
variable that depends on the image resolution $s_{k}^{2}$ according to (6).
Considering this, we can calculate (10) by applying the Bayes’ rule of
conditional probability and finding the PDF of $T_{ul,k}|r$ and $T_{w,k}$.

The PDF of $T_{ul,k}|r$ can be derived by considering that an image is
completely offloaded at the base station when all its bits arrive at the base
station. That is, when the sum of all the bits transmitted over $\left\lceil
N_{\textit{coh}}\right\rceil$ coherence intervals of length $\tau_{c}$ is
equal to (or larger than) the size of the image111For simplicity, we neglect
the pilot signals, headers, and cyclic prefixes.,

$\tau_{c}\sum_{n=1}^{\left\lceil
N_{\textit{coh}}\right\rceil}\log_{2}(1+\text{SINR}_{k}[n])|r\geq\frac{\theta
s^{2}_{k}}{\xi},$ (11)

where $\left\lceil\cdot\right\rceil$ is the least integer function. While
$\tau_{c}$ depends on the propagation environment, sending images with
resolution $s_{k}^{2}\geq 40$ kilopixels require sufficiently large values of
$N_{\textit{coh}}$ for the channel to experience a high variety of fading
realizations during the entire uplink transmission [17]. Considering this, the
left hand side of (11) can be approximated as

$\displaystyle\tau_{c}\sum_{n=1}^{\left\lceil
N_{\textit{coh}}\right\rceil}\log_{2}(1+\text{SINR}_{k}[n])|r\xrightarrow{N_{\textit{coh}}\rightarrow\infty}$
$\displaystyle\tau_{c}\left\lceil
N_{\textit{coh}}\right\rceil\,\mathbb{E}[\log_{2}(1+\text{SINR}_{k})|r]=\frac{\tau_{c}\left\lceil
N_{\textit{coh}}\right\rceil}{B_{k}}\,\bar{R}_{k}|r,$

where $\bar{R}_{k}|r$ is the conditional ergodic capacity given in (9). Since
$N_{\textit{coh}}=T_{ul,k}B_{k}/\tau_{c}$, we get

$T_{ul,k}|r\geq\frac{\theta s^{2}_{k}}{\xi}\frac{1}{\bar{R}_{k}|r}.$

Note that the minimum time to transmit an image of $\theta s_{k}^{2}/\xi$ bits
becomes a deterministic variable as $N_{\textit{coh}}\rightarrow\infty$. Thus,
in the limit, we can express the PDF of $T_{ul,k}|r$ for a user $k$ located at
distance $r$ from it serving base station as

$f_{T_{ul,k}|r}(\omega)\xrightarrow{N_{\textit{coh}}\rightarrow\infty}\Delta\left(\omega-\frac{\theta
s^{2}_{k}}{\xi\,\bar{R}_{k}|r}\right),$

where $\Delta(\cdot)$ represents the Dirac delta function.

The probability density function of $T_{w,k}$ can be derived from the state
probabilities of an M/D/1 queuing system using the Erlang’s principle of
statistical equilibrium [18, Sections 10.4.2 and 10.4.4]. For a detailed
explanation on how to derive the PDF of $T_{w,k}$, please refer to [18] and
the equations therein.

Putting all together, we can express the conditional probability of successful
computation as

$\displaystyle
p_{succ}\left(t_{k},s_{k},\lambda,\lambda_{b},\lambda_{u},\delta,P,\alpha,\epsilon,\bar{P}\;|\;r\right)$
$\displaystyle\;=\int_{0}^{\infty}\mathcal{P}(T_{w,k}+T_{s,k}\leq
t_{k}-\omega|T_{ul,k}=\omega,r)f_{T_{ul,k}|r}(\omega)d\omega$
$\displaystyle\;=\mathcal{P}\left(T_{w,k}\leq t_{k}-\frac{\theta
s^{2}_{k}}{\xi\,\bar{R}_{k}|r}-\frac{c_{1}s_{k}^{3}+c_{2}}{h}\right),$ (12)

which can be efficiently and accurately calculated by using a recursive
formula based on Fry’s equations of state [19].

###  IV-D Effective arrival rate

Depending on the system parameters, the number of images offloaded to the
server do not always correspond to the number of images that are successfully
transmitted or processed within the delay requirement $t_{k}$. Considering
this, we calculate the effective arrival rate of all the users connected to
the same server as

$\lambda_{\textit{eff}}=\lambda\int_{0}^{\infty}p_{succ}\left(t_{k},s_{k},\lambda,\lambda_{b},\lambda_{u},\delta,P,\alpha,\epsilon,\bar{P}\,|\,r\right)f_{R}(r)dr,$
(13)

and the maximum effective arrival rate with respect to $\lambda$ as

$\lambda^{*}_{\textit{eff}}\left(t_{k},s_{k},\lambda_{b},\lambda_{u},\delta,P,\alpha,\epsilon,\bar{P}\right)=\max_{\lambda}\lambda_{\textit{eff}}.$
(14)

##  V Numerical results

TABLE I: System parameters Path-loss exponent $\alpha$ | 3.7
---|---
Peak transmission power $\bar{P}$ | 23 dBm (200 mW)
Noise power density $\sigma^{2}$ | -174 dBm/Hz
Bandwidth assigned to each user $B_{k}$ | 2.1 MHz
Image height/width $s_{k}$ |  $[200,600]$ pixels
Image encoding rate $\theta$ | 24 bits/pixel
Image compression rate $\xi$ | 2
Maximum load at the server $\rho_{\max}$ |  $0.99$ Erlangs
Processing capabilities of the server $h$ | 10 TFLOPs

This section evaluates all the performance metrics presented in Section IV in
the context of realistic parameters for multi-cell, edge-intelligent systems
supporting video analytics. We also analyze and discuss the effect of the
system parameters on the accuracy, effective frame rate, and system fairness.

In all our numerical evaluations, we consider the system parameters from Table
I, corresponding to a typical cellular network [20]. Following the same
reasoning as [16], we set $\delta=1$ to maximize the average transmission rate
in the network. Moreover, we set the density of base stations to
$\lambda_{b}=0.25$ BS/km2, and consider the density of users using the same
frequency resources to be $\lambda_{u}=\lambda_{b}$, corresponding to a highly
loaded network. For the detection algorithm, we set the parameters that define
(4) and (5) to $c_{1}=7\cdot 10^{-10}$, $c_{2}=0.083$, $c_{3}=1$,
$c_{4}=1.578$, and $c_{5}=6.5\cdot 10^{-3}$, which have been shown to fit the
ground truth with a root mean square error less than $0.03$ [14]. Recall that
we assumed for simplicity that the image resolutions $s_{k}^{2}$ and the
maximum delay requirements $t_{k}$ are the same for all the users located
within the same cell.

Figure 1 shows the conditional transmission rate (9) as a function of the
distance $r$ for different fractional power control parameters $\epsilon$, and
reference powers $P$. For all cases, the transmission rate decays rapidly with
increasing $r$, and increasing $\epsilon$ leads to lower transmission rates
for the users close to the base station but to slightly higher transmission
rates for the users at the cell edge. Moreover, since we have an interference-
limited system, $P$ makes little difference. The figure also shows the
cumulative distribution function of the locations of the users within the
coverage area of a serving base station, demonstrating that most of the users
are located within the region from $0.5$ to $1.5$ km. To balance between
supporting the users at the cell edge and providing a high transmission rate
for the typical users, the rest of the numerical results are calculated with
$\epsilon=0.25$ and $P=10$ mW.

Figures 1 and 1 show the conditional probability of successful computation
(12) as a function of the distance $r$ for different image heights/widths
$s_{k}$, aggregate arrival rates $\lambda$, and delay requirements $t_{k}$.
Specifically, the image sizes $s_{k}=\\{280,430,600\\}$ correspond to the
maximum accuracy $A(s_{k})=\\{0.74,0.90,0.97\\}$ that can be supported for the
arrival rates $\lambda=\\{100,70,40\\}$, respectively, without overloading the
servers. Figure 1 shows that the server load and the delay limit $t_{k}$
determines the success probability for the users close to the base station.
Users farther away from the base station are instead limited by the uplink
transmission time, and therefore the success probability is strongly affected
by the image size. Figure 1 shows that the image size and the delay limit
determines the maximum distance of possibly successful image analytics. Within
this distance, the probability of success can be increased by decreasing the
aggregate arrival rate, but it does not help the users farther away from the
base station.

Figure 1 shows the effective aggregate arrival rate (13) as a function of the
aggregate arrival rate for different $s_{k}$ and $t_{k}$. In all cases, we
observe two distinct regions. In the first region, $\lambda_{\textit{eff}}$
grows linearly with $\lambda$, with shallower slopes for larger images sizes
and more stringent delay requirements. In this case, the loss is due to the
long transmission times of the users farther from the base station. In the
second region, $\lambda_{\textit{eff}}$ reaches its peak and then falls
rapidly, as the server becomes overloaded.

Figures 1–1 reflect that the system can be bandwidth-limited or
computationally-limited, and it is important to know in which region the
system operates to effectively decrease the losses. Now we continue with the
analysis of the performance and the fairness of multi-cell edge video
analytics.

Figure 1 shows the maximum effective arrival rate (14) – detection accuracy
(5) region of the system. We consider $B_{k}=2.1$ MHz and $B_{k}=10$ MHz to
distinguish between the bandwidth-limited and the computationally-limited
cases. The results are compared to the ideal case with maximum server load and
no losses. We see that the rate-accuracy trade-off of the computationally-
limited systems follows the shape of the ideal case, with decreasing values
under stricter deadlines $t_{k}$. However, the bandwidth-limited systems
behave differently, with an early decrease of the achievable rate as the
accuracy is increased. Thus, the bandwidth-limited systems have little
possibility to trade frame rate for increased accuracy.

Finally, to evaluate the fairness of the system, Figure 1 shows the Lorenz
curve of the probability of successful computation (12) with respect to the
location of the users, for different $s_{k}$, $\lambda$, $B_{k}$, and
$t_{k}=0.3$ seconds. The Lorenz curve gap is given in the legend. For all the
configurations, we select the values of $s_{k}$ and $\lambda$ such that the
load at server is close to $\rho_{\max}$. We see that the computationally-
limited systems (large bandwidth, smaller image size, and high arrival rate)
remain fair, even though the probability of success is not necessarily high.
On the other hand, there are bandwidth-limited systems (small bandwidth,
larger image size, and low arrival rate) that are highly unfair, with a Lorenz
curve gap close to one. In the middle we find configurations that are equally
limited by the transmission and the computation capacities, where the success
probability gradually decreases for increasing distances to the serving base
station. Overall, the system fairness can be improved by decreasing the frame
size.

Figure 1: (a)-(c) Conditional transmission rate and probability of successful
computation as a function of the distance between any randomly selected user
and its serving base station. (d) Effective aggregate arrival rate as a
function of the arrival rate. (e) Maximum effective aggregate arrival rate as
a function of the accuracy of the detection algorithm. (f) Lorenz curve with
Lorenz curve gap given in the legend. Figures (a)-(c) evaluate the
relationship between the different parameters of the video analytics, of the
uplink transmission, and of the edge server. Figures (d) and (e) analyze the
effect of the system parameters on the effective arrival rate and on the
accuracy of the detection algorithm. Finally, Figure (f) measures the
disparity of the probability that an image is successfully processed,
depending on the location of the users within the coverage area of their
serving base station.

##  VI Discussion

Motivated by the recent advances in edge computing and artificial
intelligence, we proposed a mathematical framework for modelling multi-cell,
edge-intelligent systems supporting video analytics. We expressed the
probability of successful computation and the effective arrival rate as a
function of the accuracy of the video analytics, as well as the network and
computing parameters. Our work provides the first steps towards designing
efficient resource-allocation strategies and traffic-control protocols. We
showed, for example, that different actions are necessary to improve the
effective arrival rate, depending on whether the system is operating in
bandwidth or computationally limited regimes. We also showed that the system
can be highly unfair, and thus centralized solutions may be necessary to
protect the weaker users. Further work should consider the modelling of more
flexible bandwidth allocation and scheduling solutions, as well as protocol
design.

## References

  * [1] J. Chen and X. Ran, “Deep learning with edge computing: A review.” _Proc. IEEE_ , vol. 107, no. 8, pp. 1655–1674, 2019.
  * [2] Z. Zhou, X. Chen, E. Li, L. Zeng, K. Luo, and J. Zhang, “Edge intelligence: Paving the last mile of artificial intelligence with edge computing,” _Proceedings of the IEEE_ , vol. 107, no. 8, pp. 1738–1762, 2019.
  * [3] F. Saeik, M. Avgeris, D. Spatharakis, N. Santi, D. Dechouniotis, J. Violos, A. Leivadeas, N. Athanasopoulos, N. Mitton, and S. Papavassiliou, “Task offloading in edge and cloud computing: A survey on mathematical, artificial intelligence and control theory solutions,” _Computer Networks_ , vol. 195, p. 108177, 2021.
  * [4] Q. Liu and T. Han, “Dare: Dynamic adaptive mobile augmented reality with edge computing,” in _2018 IEEE 26th International Conference on Network Protocols (ICNP)_. IEEE, 2018, pp. 1–11.
  * [5] X. Ran, H. Chen, X. Zhu, Z. Liu, and J. Chen, “Deepdecision: A mobile deep learning framework for edge video analytics,” in _IEEE INFOCOM 2018-IEEE Conference on Computer Communications_. IEEE, 2018, pp. 1421–1429.
  * [6] A. Galanopoulos, V. Valls, G. Iosifidis, and D. J. Leith, “Measurement-driven analysis of an edge-assisted object recognition system,” in _ICC 2020-2020 IEEE International Conference on Communications (ICC)_. IEEE, 2020, pp. 1–7.
  * [7] C. Wang, S. Zhang, Y. Chen, Z. Qian, J. Wu, and M. Xiao, “Joint configuration adaptation and bandwidth allocation for edge-based real-time video analytics,” in _IEEE INFOCOM 2020-IEEE Conference on Computer Communications_. IEEE, 2020, pp. 257–266.
  * [8] S. Ko, K. Han, and K. Huang, “Wireless networks for mobile edge computing: Spatial modeling and latency analysis,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 8, pp. 5225–5240, 2018.
  * [9] S. Wang, W. Guo, and M. D. McDonnell, “Distance distributions for real cellular networks,” in _2014 IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS)_. IEEE, 2014, pp. 181–182.
  * [10] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, _Stochastic geometry and its applications_. John Wiley & Sons, 2013.
  * [11] T. D. Novlan, H. S. Dhillon, and J. G. Andrews, “Analytical modeling of uplink cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 12, no. 6, pp. 2669–2679, 2013.
  * [12] K. Sriram and W. Whitt, “Characterizing superposition arrival processes in packet multiplexers for voice and data,” _IEEE journal on selected areas in communications_ , vol. 4, no. 6, pp. 833–846, 1986.
  * [13] J. Redmon and G. Jocher, _YOLOv5_ , 2021. [Online]. Available: https://github.com/ultralytics/yolov5
  * [14] Q. Liu, S. Huang, J. Opadere, and T. Han, “An edge network orchestrator for mobile augmented reality,” in _IEEE INFOCOM 2018-IEEE Conference on Computer Communications_. IEEE, 2018, pp. 756–764.
  * [15] A. Bochkovskiy, C.-Y. Wang, and H.-Y. M. Liao, “Yolov4: Optimal speed and accuracy of object detection,” _arXiv preprint arXiv:2004.10934_ , 2020.
  * [16] J. G. Andrews, F. Baccelli, and R. K. Ganti, “A tractable approach to coverage and rate in cellular networks,” _IEEE Transactions on communications_ , vol. 59, no. 11, pp. 3122–3134, 2011.
  * [17] R. P. Torres and J. R. Pérez, “A lower bound for the coherence block length in mobile radio channels,” _Electronics_ , vol. 10, no. 4, p. 398, 2021.
  * [18] V. B. Iversen, “Teletraffic engineering and network planning,” _Technical University of Denmark_ , p. 270, 2010.
  * [19] V. B. Iversen and L. Staalhagen, “Waiting time distribution in M/D/1 queueing systems,” _Electronics Letters_ , vol. 35, no. 25, pp. 2184–2185, 1999.
  * [20] A. Ghosh, J. Zhang, J. G. Andrews, and R. Muhamed, _Fundamentals of LTE_. Pearson Education, 2010.


# Machine-Learning-Driven Runtime Optimization of BLAS Level 3 on Modern
Multi-Core Systems

Yufan Xia  The Chinese University of Hong Kong
Hong Kong SAR, China
<EMAIL_ADDRESS>Giuseppe Maria Junior Barca
The University of Melbourne
Melbourne, Australia
<EMAIL_ADDRESS>
###### Abstract

BLAS Level 3 operations are essential for scientific computing, but finding
the optimal number of threads for multi-threaded implementations on modern
multi-core systems is challenging. We present an extension to the Architecture
and Data-Structure Aware Linear Algebra (ADSALA) library that uses machine
learning to optimize the runtime of all BLAS Level 3 operations. Our method
predicts the best number of threads for each operation based on the matrix
dimensions and the system architecture. We test our method on two HPC
platforms with Intel and AMD processors, using MKL and BLIS as baseline BLAS
implementations. We achieve speedups of 1.5 to 3.0 for all operations,
compared to using the maximum number of threads. We also analyze the runtime
patterns of different BLAS operations and explain the sources of speedup. Our
work shows the effectiveness and generality of the ADSALA approach for
optimizing BLAS routines on modern multi-core systems.

##  I Introduction

The linear algebra subroutines in the Basic Linear Algebra Subprograms (BLAS)
[1] form the backbone of scientific computing. Due to the critical role of
BLAS, great effort has been devoted to improving the performance its
subroutines. The automatically tuned linear algebra software (ATLAS) and its
improvements are the first batch of auto-tuning efforts on optimising linear
algebra operations; they are able to auto-tune the linear algebra operation
codes by searching over parameters like blocking factor, loop order, and
partial storage location on each specific computer architecture [2, 3, 4].
Later in the 2000s, the self-optimised linear algebra routine (SOLAR)
attempted the analytical modelling of the timing of multi-process linear
algebra operations [5, 6, 7, 8, 9, 10].

Within BLAS, Level 3 (L3) operations, which are concerned with matrix-matrix
linear algebra, are the most complex and computationally demanding [1]. In
addition to the General Matrix Multiplication (GEMM), BLAS L3 includes the
Symmetric Rank-k Update (SYR), Symmetric Rank-2k Update (SYR2), Triangular
Matrix-Matrix Multiplication (TRMM), and Triangular Rank-k Update (TRSM)
subroutines.

While there has been significant work in optimizing single-threaded CPU
performance of BLAS L3 by fine-tuning the size of matrix blocks for various
system architectures, similar efforts for optimizing these routines
specifically for multi-core CPU have been less frequent due to the high
complexity and large variety of the underlying computer architectures.

Work has been carried out by using established optimization techniques such as
parameter tuning and blocking[11, 12]. The Parallel Linear Algebra for
Scalable Multi-Core Architectures (PLASMA) [13] library uses polynomial
regression to optimise the number of threads or block size used based on
empirical features, obtaining performance comparable to the Intel® MKL for
Cholesky decomposition. Peise _et al._ adopted a polynomial regression to
model the dense linear algebra run times, and they applied additional
techniques to boost the performance of this simple model [14].

However, choosing the number of threads that minimises the execution time of a
given BLAS remains challenging and largely unsolved due to the underlying
diversity and complexity of modern shared-memory computer architectures.

To tackle this challenge, recent research conducted by Xia _et al._ [15],
which was integrated in the open-source Architecture and Data-Structure Aware
Linear Algebra library, has employed a systematic machine learning (ML)
methodology to significantly reduce the runtime of Single-precision General
Matrix Multiplication (SGEMM) operations. This approach leverages an ML model
to dynamically select the optimal number of threads for minimizing the
execution time of GEMM for a specific input configuration and computer system
architecture. The ML model itself undergoes training during the installation
phase, and this training process is customized to the particular computer
system architecture in use. Subsequently, during runtime, the ML model is
employed to make predictions regarding the most efficient number of threads
for a given task. An important advantage of this approach is its adoption of
existing GEMM implementations, treating them as black boxes. This work has
inspired a novel auto-tuning method in Graph Neural Network training [16].

In this study, we have expanded the runtime optimization capabilities of the
ADSALA library to encompass all single- and double-precision BLAS L3
operations, and we present the performance results achieved. Furthermore, we
have incorporated an automatic model selection feature during the installation
process for each distinct BLAS subroutine. This addition aids in identifying
and using the most appropriate machine learning algorithm for each subroutine
on every machine where the library is installed. We test our method on two HPC
platforms with Intel and AMD processors, using multi-threaded MKL and BLIS as
baseline BLAS implementations. Our software implementation enables us to speed
up all BLAS L3 subroutines, notwithstanding the runtime overhead of ML
evaluations.

On the Gadi supercomputer located at the National Computational Infrastructure
(NCI) and on the Setonix supercomputer located at the Pawsey Supercomputing
Centre (please refer to section V for details about NCI and Pawsey), we
achieve speedups of 1.5 $\times$ to 3.0 $\times$ compared to using the maximum
number of threads with hyperthreading enabled or disabled. For comparison, the
most recent method relevant to our study employs pure polynomial regression to
model and determine the number of threads for the PLASMA QR routine[13]. This
method reported average speed improvements of 1%, 13%, and 27% across three
different platforms. The sizes of the matrices used in these tests ranged from
2000 to 7000, which aligns closely with the range we tested in our study.

We also analyze the runtime patterns of different ADSALA BLAS L3 runs and
discuss the sources of speedup. The ADSALA library is provided as open-source
implementation for the community to use and extend.

The remainder of this article is structured as follows. We provide background
information concerning BLAS L3 operations and ML in Section II. We review the
improved software design in Section III, and discuss the adopted ML methods in
Section IV. Section V details the experimentation platform and settings. We
present and analyze the performance of the ML models and software speedup in
Section VI. Section VII concludes.

TABLE I: Specifications of BLAS level III subroutines. |  |  |  |  |  |  |
---|---|---|---|---|---|---|---
|  |  Matrix A |  |  Matrix B |  |  Matrix C |
|  dims |  shape |  type |  shape |  type |  shape |  type
GEMM |  3 |  m $\times$ k |  regular |  k $\times$ n |  regular |  m $\times$ n |  regular
SYMM |  2 |  m $\times$ m |  symmetric |  m $\times$ n |  regular |  m $\times$ n |  regular
SYRK |  2 |  n $\times$ k |  regular |  n $\times$ k |  regular |  n $\times$ n |  symmetric
SYR2K |  2 |  n $\times$ k |  regular |  n $\times$ k |  regular |  n $\times$ n |  symmetric
TRMM |  2 |  m $\times$ m |  triangular |  m $\times$ n |  regular |  — |  —
TRSM |  2 |  m $\times$ m |  triangular |  m $\times$ n |  regular |  — |  —
TABLE II: Comparisons of ML model characteristics.

|  |  |  |
---|---|---|---|---
Model Catagories |  Models |  Parametric |  Good with Data Imbalance |  Data Size Requirement
Linear Models |  Linear Regression |  Yes |  No |  Medium
|  ElasticNet |  |  |  Medium
|  Bayesian Regression |  |  |  Small
Tree Based Models |  Decision Tree |  No |  Yes |  Medium
|  XGBoost |  |  |
|  AdaBoost |  |  |
|  Random Forest |  |  |
|  LightGBM |  |  |
Other Models |  SVM Regressor |  No |  No |  Small
|  KNN Regressor |  |  |  Medium

##  II Background

In this Section, we introduce BLAS L3 subroutines, ML algorithms, and related
data preprocessing techniques.

###  II-A BLAS Level III Subroutines

BLAS L3 includes six matrix-matrix operations, including GEMM, SYMM, SYR,
SYR2, TRMM, and TRSM [1]. These operations have different number of matrix
operands, matrix shapes, and some matrix operands are required to be symmetric
or triangular. These requirements are listed in Table I. These subroutines
provide a unified interfaces for different BLAS implementations, thus enables
users to switch between different BLAS implementations without changing their
code.

The performance of BLAS L3 subroutines is highly dependent on the underlying
implementation and the system architecture. For example, the GEMM routine is
the most commonly used BLAS Level III subroutine, and it is usually the most
optimised one for BLAS implementations. However, as reported, even GEMM has
the issue where the maximum thread number, which we defined to be equal to the
number of available CPU cores times the hyper-threading level, may not provide
the best performance [15]. While this issue has not been reported to be
present in other BLAS subroutines, in later sections, we will show that the
other BLAS L3 subroutines show a similar behaviour.

###  II-B Machine Learning Algorithms

To build a high-performance ML model, we need to sift through suitable
potential candidates. We provide a brief overview of the dataset we will use
(for more details on data collection, please refer to subsection IV-B) to help
justify our selection of candidate models that best fit the data
characteristics. The dataset for a given BLAS L3 operation typically comprises
approximately $10^{3}$ data points spanning 4-15 dimensions, with most data
features exhibiting a skewed distribution. Despite being a relatively small
dataset with low-dimensional data, the relationship between the features and
the label can be quite complex due to the expected non-linear relation
stemming from the polynomial time complexity of most BLAS L3 operations and
the intricacies of multi-core computer architectures. Given these attributes,
we present several ML algorithms as potential candidates in Table I.

Linear algorithms produce simple but high speed ML models [17]. They usually
suffer from low predictive accuracy for non-linear mappings. Due to the nature
of our task, evaluation speed is an essential factor, so linear ML models may
enable pragmatic BLAS L3 speedups. We include linear regression, ElasticNet
and Bayesian regression as our candidates [18, 19].

We also include both Decision Tree and tree ensemble models as candidates. The
Decision Tree algorithm is non-parametric, it uses a set of rules to map each
data instance to a discrete class or a continuous value [20]. Random Forest,
AdaBoost, XGBoost and LightGBM are more complex algorithms who are built upon
an ensemble of Decision Trees. These tree ensemble algorithms generally can
reduce bias and variances compared to a single Decision Tree and thus can
produce a more promising result but with higher evaluation time.

Support Vector Machine (SVN) and k-Nearest-Neighbors (kNN) are non-parametric
algorithms but are designed to work well for high dimensions and unknown data
structures, respectively [21]. Since we have low data dimensionality, they
might not demonstrate their advantages when learning from our dataset; kNN is
also known to be slow in evaluation [22]. We will evaluate their performance
for completeness.

(a) The installation workflow of ADSALA GEMM. Upon ADSALA installation, the
library performs two sub-parts shown in the diagram. In the end, two files
containing the configurations together with the production-ready ML model will
be saved for later use at runtime.

(b) The runtime workflow of ADSALA for BLAS level III subroutines.
Configuration files and trained ML models outputted during installation (see
Fig. 0(a)) are used by this runtime library.

Figure 1: The software design of ADSALA.

###  II-C Data Preprocessing Techniques

Data outliers can affect the predictive performance of ML models. While
statistical methods can effectively remove global outliers, they often fail to
detect local outliers. the Local Outlier Factor (LOF), a density-based method,
overcomes this limitation by assigning each data point a degree of isolation
score based on its surrounding points [23, 24]. We use LOF to identify and
remove these outliers from our dataset.

The Yeo-Johnson transformation is a data transformation technique that remaps
each feature value such that the feature distribution is near-Gaussian, thus
improving the predictive performance of ML algorithms that assume normality,
such as linear regression, ElasticNet, and Bayesian regression. Unlike the
original Box-Cox transformation, which requires all feature values to be
positive, Yeo-Johnson accepts non-positive values and provides stable
parameter estimation [25, 26]. Yeo-Johnson transformation uses a parameter
$\lambda$ to control the strength of transformation. We apply Yeo-Johnson with
maximum likelihood estimation (MLE) for parameter estimation for the data
transformation, thereby automating the ML workflow.

##  III Software Workflow

Our software, upon installation, collects runtime data for each of BLAS L3
subroutine and different combinations of thread numbers, input matrices
dimensions, and BLAS packages availability. This data trains an ML model to
predict the optimal thread number for given dimensions of the GEMM input. The
final product is an ADSALA library that dynamically selects the optimal thread
number at runtime, adapting to different HPC platforms, BLAS packages, and
subroutines.

Figure 0(a) and 0(b) illustrate our software design procedures. The software
comprises data gathering, model training, and a runtime library. The first two
parts are executed at installation time, while the third part is typically
used at runtime when the user program links to our library.

###  III-A Installation Workflow

During installation, the software gathers training data through
experimentation. It quasi-randomly (_vide infra_) samples from the domain of
the specific subroutine matrices dimensions, passes them to a timing program
that runs and times the corresponding BLAS operations. The timings are stored
as training data for the ML model. The timing data is then preprocessed and
utilised to tune the hyper-parameters of the ML model.

###  III-B Runtime Workflow

Upon instantiation, the ML models and its configuration are loaded into
memory. At runtime, when the user program calls a BLAS function, the ML model
predicts the optimal thread number on-the-fly and runs the implementation with
that thread number. The ML model and configuration are wrapped in a C++ class
that can be destroyed after the last BLAS call to free memory.

To avoid redundant ML model evaluations for two consecutive BLAS call with
identical input dimensions, our software remembers the input to the last BLAS
call and its correlation ML prediction. If the current BLAS input dimensions
are the same as the previous, the software will simply read and apply the
predictions from the responsible class attributes without re-evaluation.

##  IV Machine Learning Methods

In this section, we build upon the original ADSALA proof-of-concept library
[15] and extend it to support all BLAS L3 subroutines. We first introduce the
mechanism of ML predictions, then we detail the methods for data collection,
data preprocessing, model selection, and model training process in our
software. Since our goal is a library that works on different HPC
architectures and BLAS packages, the following methodologies can be applied to
any HPC architectures and BLAS packages.

###  IV-A Mechanism for Predictions

For a given BLAS subroutine and a given combination of input dimensions, the
software selects the optimal number of threads by first predicting the
duration of BLAS subroutine associated with the possible number of threads. It
then selects the thread number with the shortest predicted runtime for the
ensuing BLAS subroutine execution.

###  IV-B Data Gathering

Two paragraphs that were in this section are combined to eliminate repetitions
mentioned by R1

In order to ensure that our method is effective for BLAS L3 subroutines across
a variety of matrix shapes and sizes within memory constraints, including
slim/square and big/small matrices, we need to sample domains that are evenly
distributed across the space. To achieve this, we employ a scrambled Halton
sequence to produce a quasi-random sequence with low discrepancy for data
sampling [27]. The upper bound of the sum size of matrices to 500 MB111Some
BLAS subroutines have their output matrix overwrite the input matrix (TRMM and
TRSM). We calculate the sum memory size as the sum of input/output matrices,
ignoring the overwritten matrices.. Given that the samples can have two or
three domains, we opt for the scrambled Halton sequence over the regular
Halton sequence to reduce the correlation between dimensions [27]. We utilize
bases 2, 3, and 4 for generating the sequence for dimensions $m$, $k$, and $n$
and bases 2, 3 for dimensions $m$, $n$ (or $n$, $k$). These samples are then
input into the timing program to gather their runtime data.

###  IV-C Feature Engineering and Data Preprocessing

Table III shows the features used for the ML model. We created two set of
features for subroutines with three matrix size parameter and two matrix size
parameter with respect to the computational complexity, their memory
footprint, and the multi-thread speedup.

The execution time for BLAS subroutines is a function of the matrix dimensions
and the number of threads ($nt$), and it varies depending on the architecture.
For three-dimension subroutines, the terms $m*k$, $k*n$, $m*n$, and
$m*k+k*n+m*n$ correspond to the sizes of matrices $A$, $B$, $C$, and the total
memory size in single-/double-precision words, respectively. Similar terms are
constructed for two-dimension subroutines as shown in the second group in
Table III. These memory terms are directly related to the number of memory
operations and hence, influence the runtime of various BLAS subroutines.
Specifically, these terms are dominant in the serial runtime for smaller
matrix dimensions. In three-dimension subroutines, the cubic term $m*k*n$ is
proportional to the number of floating-point operations performed and tends to
dominate the runtime in serial execution for larger matrix dimensions. In
parallel execution, the FLOP workloads are distributed across threads,
resulting in terms like ${m*n*k}/nt$.

For feature transformation, outlier removal, feature selection, and hyper-
parameter tuning stages, we use the same methodologies as in the original
ADSALA library that helps improving the speed and also predictive performance
of the model[15]. We employ the Yeo-Johnson transformation to approximate the
distribution of features to a near-Gaussian distribution, which is beneficial
for model learning.222In our tests, we observed an enhancement in the
performance of these models when the Yeo-Johnson Transformation was applied,
resulting in a 10-20% decrease in the Root Mean Square Error (RMSE) for Linear
Regression. The application of it does not affect the performance of other
candidate models much. This transformation identifies the most suitable
parameter values for $\lambda$ and adjusts the transformation’s impact on each
feature using the MLE method [26]. Following this transformation, we carry out
a standardisation process on features to ensure they all operate on a similar
scale [28]. Subsequently, we eliminate features that have correlation
coefficients with other features exceeding a threshold of $80\%$ to eliminate
the potentially redundant features among the candidates. For each correlation
feature pair, we remove the feature with the larger total correlation with the
other features. Following this, the hyper-parameter tuning is performed for
all models to compare model performance.

TABLE III: List of available features for BLAS subroutines with two or three matrix dimensions. nt stands for number of threads. | Three Dimensions $m$, $k$, $n$ |  | Two Dimensions $m$, $n$ |
---|---|---|---|---
1 | m | ${\text{m}}/{\text{nt}}$ | m |
2 | k | ${\text{k}}/{\text{nt}}$ | n |
3 | n | ${\text{n}}/{\text{nt}}$ | nt |
4 | nt | ${\text{m*k}}/{\text{nt}}$ | m*n |
5 | m*k | ${\text{m*n}}/{\text{nt}}$ | memory_footprint |
6 | m*n | ${\text{k*n}}/{\text{nt}}$ | m/nt |
7 | k*n | ${\text{m*k*n}}/{\text{nt}}$ | n/nt |
8 | m*k*n |  memory_footprint / nt | m*n/nt |
9 | memory_footprint |  | memory_footprint / nt |

###  IV-D Model Selection

Given that our method’s ultimate goal is to achieve optimal runtime speedup,
the model selection process incorporates both the predictive performance of
the model and the speed of model evaluation with a potential trade-off between
these two. This approach is akin to the proof-of-concept work conducted with
ADSALA [15] where the speedup is estimated using the formula: s =
toriginaltADSALA+teval In this context, $t_{\text{ADSALA}}$ is the runtime of
the BLAS subroutine with predicted threads, $t_{\text{eval}}$ is the ML model
evaluation time, and $t_{\text{original}}$ is the runtime with maximum
threads. The evaluation time $t_{\text{eval}}$ is measured by averaging
multiple runs. The ML model with the highest average estimated speedup $s$
across all BLAS subroutines is selected.

##  V Experimentation Information

This Section presents the supercomputing platforms used for experimentation
and the setup of the experiments on those platforms.

###  V-A Experimentation Platforms

We conducted our experiments on two supercomputing platforms.

####  V-A1 Setonix

Setonix is a supercomputer housed at the Pawsey Supercomputing Research Centre
in Australia. Each of its compute nodes has specific features, as depicted in
Fig. 2.

  * •

Each compute node is equipped with two CPU sockets, each housing an AMD® EPYC
64-Core Milan CPU (2.55 GHz). The CPUs support hyper-threading, enabling up to
256 threads to run concurrently per compute node. Each Milan CPU comprises
eight modules, with each module containing eight Zen 3 cores and a dedicated
32 MB level three cache that is shared among these cores.

  * •

The system has a memory capacity of 256GB, organized into eight NUMA domains,
with four NUMA domains per socket. Each socket supports eight memory channels.

Figure 2: A schematic diagram for the 2-socket EPYC CPU configuration on
Setonix.

####  V-A2 Gadi

Gadi is a supercomputer located at the Australian National Computational
Infrastructure. Figure 3 illustrates the features of each of Gadi’s compute
nodes:

  * •

Each compute node has two CPU sockets, each containing an Intel® Xeon 24-Core
Cascade Lake CPU (Platinum 8274, 3.2 GHz). The CPUs support hyper-threading,
enabling up to 96 threads to run concurrently per compute node.

  * •

The system has a memory capacity of 192GB, organized into four NUMA domains,
with two NUMA domains per socket. Each socket supports six memory channels.

Figure 3: A schematic diagram of the 2-socket Cascade Lake CPU configuration;
sockets are connected using Intel® UPI (Ultra Path Interconnect).

###  V-B Experimentation Setup

For the Setonix supercomputer, which uses AMD® processors, we will use the
performance of BLIS333 https://developer.amd.com/amd-aocl/blas-library/ as the
baseline for measuring performance improvements. For the Gadi supercomputer,
which uses Intel® processors, we will use the performance of MKL444
https://www.intel.com/content/www/us/en/develop/documentation/get-started-
with-mkl-for-dpcpp/top.html as the baseline for measuring performance
improvements. For convenience in comparion of performance with the proof-of-
concept ADSALA work, we will use the same experimental settings which can be
found in [15].

TABLE IV: Model performance and estimated speedups for ML models on Setonix. subroutine | best_model
---|---
dgemm | LinearRegression
dsymm | XGBRegressor
dsyr2k | XGBRegressor
dsyrk | XGBRegressor
dtrmm | XGBRegressor
dtrsm | XGBRegressor
sgemm | XGBRegressor
ssymm | XGBRegressor
ssyr2k | XGBRegressor
ssyrk | XGBRegressor
strmm | XGBRegressor
strsm | XGBRegressor

TABLE V: Model performance and estimated speedups for ML models on Gadi. subroutine | best_model
---|---
dgemm | BayesianRidge
dsymm | RandomForestRegressor
dsyr2k | XGBRegressor
dsyrk | LinearRegression
dtrmm | XGBRegressor
dtrsm | LinearRegression
sgemm | XGBRegressor
ssymm | RandomForestRegressor
ssyr2k | XGBRegressor
ssyrk | XGBRegressor
strmm | XGBRegressor
strsm | LinearRegression

TABLE VI: Detailed statistics for model performance and estimated speedups for
ML models on Gadi. The bold values of each column represents the best entry in
that column.

dgemm Gadi
---
| Model |  Normalised Test RMSE |  Ideal mean speedup |  Ideal aggregate speedup |  Model evaluation time in $\mu$s |  Estimated mean speedup |  Estimated aggregate speedup
| Linear Regression |  0.93 |  1.13 |  1.01 |  15.3 |  1.12 |  1.01
| ElasticNet |  1.00 |  1.07 |  0.94 |  10.63 |  1.07 |  0.94
| Bayes Regression |  0.93 |  1.13 |  1.01 |  8.11 |  1.13 |  1.01
| Decision Tree |  0.32 |  0.82 |  0.36 |  8.02 |  0.82 |  0.36
| Random Forest |  0.24 |  1.08 |  0.97 |  983.09 |  1.01 |  0.94
| AdaBoost |  0.52 |  0.84 |  0.46 |  89.77 |  0.84 |  0.46
| KNN |  0.28 |  1.07 |  0.95 |  6449.32 |  0.78 |  0.78
| XGBoost |  0.23 |  1.1 |  0.98 |  290.76 |  1.08 |  0.97

dsymm Gadi
---
| Model |  Normalised Test RMSE |  Ideal mean speedup |  Ideal aggregate speedup |  Model evaluation time in $\mu$s |  Estimated mean speedup |  Estimated aggregate speedup
| Linear Regression |  1.00 |  1.49 |  1.86 |  6.49 |  1.49 |  1.86
| ElasticNet |  1.00 |  1.45 |  1.8 |  4.78 |  1.45 |  1.8
| Bayes Regression |  1.00 |  1.49 |  1.86 |  4.65 |  1.49 |  1.86
| Decision Tree |  0.46 |  1.07 |  0.54 |  4.56 |  1.07 |  0.54
| Random Forest |  0.15 |  1.72 |  1.89 |  563.59 |  1.62 |  1.86
| AdaBoost |  0.61 |  0.76 |  0.36 |  69.88 |  0.75 |  0.36
| KNN |  0.18 |  1.76 |  1.95 |  2775.13 |  1.42 |  1.79
| XGBoost |  0.18 |  1.65 |  1.82 |  1046.46 |  1.5 |  1.77

ssyrk Gadi
---
| Model |  Normalised Test RMSE |  Ideal mean speedup |  Ideal aggregate speedup |  Model evaluation time in $\mu$s |  Estimated mean speedup |  Estimated aggregate speedup
| Linear Regression |  1.00 |  0.95 |  0.89 |  6.85 |  0.95 |  0.89
| ElasticNet |  1.00 |  0.99 |  0.96 |  5.07 |  0.99 |  0.96
| Bayes Regression |  1.00 |  0.95 |  0.89 |  4.97 |  0.95 |  0.89
| Decision Tree |  0.33 |  0.56 |  0.28 |  4.89 |  0.56 |  0.28
| Random Forest |  0.09 |  1.2 |  0.81 |  2324.21 |  0.96 |  0.77
| AdaBoost |  0.48 |  0.53 |  0.26 |  62.94 |  0.52 |  0.26
| KNN |  0.11 |  1.2 |  0.85 |  1760.56 |  1 |  0.81
| XGBoost |  0.09 |  1.15 |  0.82 |  446.43 |  1.08 |  0.81

strsm Gadi
---
| Model |  Normalised Test RMSE |  Ideal mean speedup |  Ideal aggregate speedup |  Model evaluation time in $\mu$s |  Estimated mean speedup |  Estimated aggregate speedup
| Linear Regression |  1.00 |  1.04 |  0.96 |  8.84 |  1.04 |  0.96
| ElasticNet |  1.00 |  1.02 |  0.95 |  6.39 |  1.02 |  0.95
| Bayes Regression |  1.00 |  1.04 |  0.96 |  5.72 |  1.04 |  0.96
| Decision Tree |  0.31 |  0.64 |  0.23 |  4.94 |  0.64 |  0.23
| Random Forest |  0.06 |  1.2 |  1.1 |  2191.65 |  0.96 |  1.01
| AdaBoost |  0.46 |  0.58 |  0.27 |  119.84 |  0.57 |  0.27
| KNN |  0.08 |  1.2 |  1.09 |  1688.15 |  1 |  1.02
| XGBoost |  0.07 |  1.21 |  1.05 |  1354.14 |  1.04 |  0.99

(a) Setonix

(b) Gadi

Figure 4: Heatmap of the optimal number of threads on Setonix and Gadi,
concerning all BLAS level III subroutines except GEMM. The horizontal and
vertical axes use a square root scale.

(a) Setonix

(b) Gadi

Figure 5: Heatmap of the optimal number of threads on Setonix and Gadi. The
horizontal and vertical axes use a square root scale. The dashed lines on each
sub-graph are contour lines of the sampling domain with each label showing the
value of the third dimension.

(a) Setonix

(b) Gadi

Figure 6: Heatmap of the testing speedup with respect to matrix sizes,
concerning all BLAS level III subroutines except GEMM. The horizontal and
vertical axes use a square root scale. TABLE VII: Speedup statistics on
Setonix and Gadi with hyper-threading

(a) Setonix |  Mean |  Std |  Min |  25% |  50% |  75% |  Max
---|---|---|---|---|---|---|---
dgemm |  1.54 |  0.66 |  0.87 |  1.16 |  1.29 |  1.74 |  4.79
sgemm |  1.32 |  0.41 |  0.76 |  1.05 |  1.18 |  1.37 |  9.05
dsymm |  2.89 |  1.80 |  0.82 |  1.38 |  2.36 |  4.09 |  8.46
ssymm |  2.22 |  1.69 |  0.34 |  1.04 |  1.60 |  2.88 |  7.42
dsyr2k |  1.48 |  0.46 |  0.94 |  1.17 |  1.36 |  1.56 |  3.59
ssyr2k |  1.53 |  0.49 |  0.75 |  1.20 |  1.43 |  1.72 |  3.79
dsyrk |  1.46 |  0.46 |  0.92 |  1.15 |  1.28 |  1.69 |  3.54
ssyrk |  1.55 |  0.42 |  0.98 |  1.22 |  1.46 |  1.82 |  2.80
dtrmm |  1.61 |  0.82 |  0.87 |  1.13 |  1.27 |  1.78 |  6.51
strmm |  1.67 |  0.92 |  0.94 |  1.10 |  1.34 |  1.81 |  7.13
dtrsm |  1.71 |  1.07 |  0.99 |  1.13 |  1.27 |  1.77 |  7.43
strsm |  1.68 |  1.21 |  0.98 |  1.12 |  1.34 |  1.85 |  12.38

(b) Gadi |  Mean |  Std |  Min |  25% |  50% |  75% |  Max
---|---|---|---|---|---|---|---
dgemm |  1.27 |  0.55 |  0.41 |  1.01 |  1.13 |  1.29 |  4.25
sgemm |  1.07 |  0.70 |  0.88 |  1.00 |  1.00 |  1.02 |  3.01
dsymm |  2.28 |  1.89 |  0.88 |  1.18 |  1.43 |  2.50 |  12.08
ssymm |  2.16 |  1.98 |  0.99 |  1.18 |  1.30 |  2.22 |  11.05
dsyr2k |  1.28 |  0.38 |  0.90 |  1.08 |  1.20 |  1.37 |  4.27
ssyr2k |  1.47 |  0.60 |  0.41 |  1.14 |  1.26 |  1.61 |  4.52
dsyrk |  1.40 |  0.36 |  0.86 |  1.20 |  1.30 |  1.47 |  3.98
ssyrk |  1.65 |  0.47 |  1.11 |  1.29 |  1.45 |  1.96 |  3.03
dtrmm |  1.30 |  0.24 |  1.03 |  1.16 |  1.25 |  1.34 |  2.75
strmm |  1.35 |  0.26 |  0.98 |  1.17 |  1.30 |  1.47 |  2.78
dtrsm |  1.31 |  0.25 |  0.75 |  1.17 |  1.29 |  1.41 |  2.86
strsm |  1.40 |  0.40 |  0.74 |  1.19 |  1.33 |  1.47 |  4.40

(a) Setonix

(b) Gadi

Figure 7: Heatmap for speedups with respect to matrix dimensions on Setonix and Gadi. The horizontal and vertical axes use a square root scale. The dashed lines on each sub-graph are contour lines of the sampling domain with each label showing the value of the third dimension. TABLE VIII: Profiling results on Gadi with GEMM, SYMM, and SYRK. m, k, n |  Thread Count |  Total Time (s) |  Thread Sync (s) |  Kernel Call (s) |  Data Copy (s)
---|---|---|---|---|---
dgemm 64,2048,64 no ML |  96 |  11.001 |  4.806 |  0.010 |  2.060
dgemm 64,2048,64 with ML |  16 |  4.430 |  1.453 |  0.001 |  0.887
sgemm 64,2048,64 no ML |  96 |  7.698 |  2.716 |  0.005 |  1.133
sgemm 64,2048,64 with ML |  5 |  2.701 |  0.155 |  0.001 |  0.301
dsymm 248,39944 no ML |  96 |  29.063 |  11.627 |  0.496 |  9.783
dsymm 248,39944 with ML |  25 |  22.067 |  4.962 |  0.441 |  9.698
ssymm 2759,41681 no ML |  96 |  27.299 |  8.986 |  1.208 |  10.605
ssymm 2759,41681 with ML |  12 |  15.597 |  2.085 |  0.241 |  8.254
dsyrk 124,160163 no ML |  96 |  36.060 |  8.545 |  0.005 |  13.672
dsyrk 124,160163 with ML |  43 |  34.371 |  7.543 |  0.001 |  12.529
ssyrk 175,15095 no ML |  96 |  65.844 |  62.666 |  1.639 |  1.103
ssyrk 175,15095 with ML |  46 |  4.581 |  1.323 |  1.744 |  1.263

##  VI Performance Analysis

In this section, we first provide a visual representation and succinct
explanation of our datasets. We then present the results of our machine
learning model selection process. Finally, we evaluate the performance of our
method and its software implementation using our test datasets. All analysis
are performed on both Setonix and Gadi supercomputers and for all BLAS L3
subroutines with double and single precision.

###  VI-A Datasets, Training Time, and Data Visualization

The size of the training dataset is 1000-1200 for each of the BLAS III
subroutines on each HPC platform. It is believed that the size of the dataset
is sufficient, since it is generally observed that the validation performance
does not improve significantly with more training set data.

The data gathering on Setonix and Gadi required approximately 100 node hours
for each subroutine. It was observed that, given the same upper limit of
matrix size (specifically 500 MB), the double precision subroutines took more
time to collect data compared to their single precision counterparts.

Figure 5 contains two subfigures that depict the distribution of the optimal
number of threads for GEMM on Setonix and Gadi, respectively. The patterns for
sGEMM and dGEMM appear to be quite similar on Setonix, but show more variation
on Gadi. A common trend observed is that irregular GEMM calls, where at least
one of the dimensions $m$, $k$, or $n$ is small, often lead to suboptimal
performance. Also, we can observe patches of abnormal area where choices of
the optimal number of threads is drastically different from the surrounding
area.

The optimal number of threads for SYMM, SYRK, SYR2K, TRMM, and TRSM on Setonix
and Gadi is displayed in Fig. 4. The patterns for double precision and single
precision subroutines are mostly similar, but there are some exceptions. Also,
the distribution of the optimal number of threads is more uneven than for
GEMM. For SYRK, TRMM, and TRSM on Setonix, many calls have optimal number of
threads higher than the number of physical cores (hyperthreads). On the
contrary, for SYRK, SYR2K, and TRMM on Gadi, almost all calls have optimal
number of threads lower than the number of physical cores. There are also
abnormal areas where the optimal number of threads is drastically different
from the surrounding area, similar to GEMM.

The complex performance patterns of BLAS level III subroutines are shown by
these results; our ML models aim to learn these patterns and select the best
number of threads.

As mentioned in Section IV-C, we use stratified sampling to split the data set
for model training and testing, with 15 $\%$ of the data set as the test set.

###  VI-B Model Performance and Selection

The best algorithms chosen by our algorithm of selecting the ML model with the
highest average estimated speedup s for all BLAS subroutines are shown in
Table V and V. We can see that only four algorithms are chosen at least once,
with XGBoost still being the most common option. The results also indicate
that the best algorithms for Setonix and Gadi are sometimes different, which
is expected since the two platforms have different architectures. We have also
included in Table VI detailed statistics supporting the selection of most
suitable models for four subroutines on Gadi.

We regret that we have not completed integrating all four models in our code
base. For now, we only use XGBoost as the ML model for both Setonix and Gadi.
Since the estimated speedup only differ by less than 10 %, we expect the
speedup difference will not be significant. We will integrate the other three
models in the future.

We used separate datasets for each subroutines from the training and test
datasets test ADSALA on BLAS L3 subroutines. These datasets consist of 100-120
data points sampled with a scrambled Halton sequence within the same domain as
how its training dataset is sampled. This is expected to ensure a more
uniform, low discrepancy sampling of the data used for the performance
analysis of our software. We compared the subroutine runtime with our ADSALA
software to the runtime of original using maximum number of threads. Please
note that the speedup results include the model evaluation time during
runtime.

The testing shows that our ADSALA software can reliably improve the
performance of BLAS L3 subroutines on both Setonix and Gadi. Table VII shows
the statistics of the speedup results. The mean speedup on Setonix is
generally higher with also higher standard deviation than on Gadi. Precision-
wise, double precision subroutines generally have higher speedups than their
single precision counterparts. Subroutine-wise, SYMM has the highest mean
speedup on both platforms, while GEMM and SYR2K has the lowest mean speedup on
Setonix and Gadi, respectively.

Figure 7 visualises the GEMM speedup distribution on Setonix and Gadi as a 3D
heatmap from three angles, where red indicates acceleration and blue indicates
deceleration. We can see that the speedup distribution pattern resembles the
pattern of room for improvements, and the speedup generally decreases as three
dimensions get larger.

Figure 6 shows the SYMM, SYRK, SYR2K, TRMM, and TRSM speedup distribution with
respect to matrix dimensions on Setonix and Gadi. We can observe that the
speedup patterns also resembles that in Fig. 4. The double/single subroutine
pairs generally show similar speedup patter. However, speedup shows various
patterns for different subroutines. Also, it does not always follows the rule
that the speedup decreases as the matrix dimensions get larger.

###  VI-C Performance Assessment using Software Implementation

We use Table VIII to explain the large speedups of ADSALA over the original
BLAS package in some cases. These cases have considerable speedup and are
selected from the test set for software testing.

We use Intel${}^{\circledR}AdvisorandIntel^{\circledR}$ Vtune to profile these
two GEMMs on Gadi, repeating each matrix multiplication 100 times with
different input values. Table VIII shows the time breakdown of GEMM calls. The
SGEMM wall-time consists of three main components:

  1. 1.

Data copies. The BLAS runtime uses a buffer for each thread as a workspace and
copies blocks of matrix operands into it before computation. The copy time
depends on the matrix sizes, memory placement, and thread number.

  2. 2.

Thread synchronization.

  3. 3.

BLAS L3 kernel calls, where the FLOPs are performed. The kernels depend on the
thread number and the matrix block sizes. They can be compute-bound for large
blocks.

Generally speaking, the speedup of ADSALA is due to the reduction of all three
parts. Because of the relatively small problem size, the most significant time
consumption is on thread synchronization most of the time. The thread sync
time can be improved from around 30 % to more than 50 times. Data copy has the
second largest time consumption, and the speedup on it also contributes to the
reduction of total runtime. Although the kernal call is consumes minimal time,
our number threads selected by ML can also reduce the time consumption on it.

##  VII Conclusions

We presented an extension to the ADSALA library where only GEMM was optimized.
We extend its approach to fit other BLAS L3 routines with different input
matrices. While keeping the same method for selecting the ML algorithm, we
found the optimal algorithm is architecture and subroutine dependent.

Out of the six BLAS L3 subroutines, we improved on all of them, with speedups
ranging from 1.1 to 3.0. We also analyzed visualise and analyse the speedup
patterns of different BLAS subroutines and HPC platforms. By profiling, we
explained the large speedups in some cases and showed that the source of
speedup is the reduction of all three main components of the multi-thread
subroutine runtime.

Our method is applicable to runtime parameter predictions, including multi-
thread BLAS I, II, and LAPACK operations, which are sensitive to prediction
duration. Future work will investigate automatic processor selection, such as
the feasibility of running BLAS/LAPACK operations on either a GPU or CPU. At
present, our approach is applicable to problems that satisfy certain
conditions: a well-defined objective function, a finite and discrete search
space, a set of relevant features, and sufficient data for training and
validating ML models.

Our work shows the effectiveness and generality of the ADSALA approach for
optimizing BLAS routines on modern multi-core systems. For future works, we
plan to extend our ML-driven runtime thread selection approach to other BLAS
operations and to a more diverse set of computing hardware, including GPU
accelerators and heterogeneous systems.

## References

  * [1] L. S. Blackford, A. Petitet, R. Pozo, K. Remington, R. C. Whaley, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry _et al._ , “An updated set of basic linear algebra subprograms (blas),” _ACM Transactions on Mathematical Software_ , vol. 28, no. 2, pp. 135–151, 2002.
  * [2] R. C. Whaley and J. J. Dongarra, “Automatically tuned linear algebra software,” in _SC’98: Proceedings of the 1998 ACM/IEEE conference on Supercomputing_. IEEE, 1998, pp. 38–38.
  * [3] J. A. Gunnels, G. M. Henry, and R. A. Geijn, “A family of high-performance matrix multiplication algorithms,” in _International Conference on Computational Science_. Springer, 2001, pp. 51–60.
  * [4] R. Vuduc, J. W. Demmel, and J. Bilmes, “Statistical models for automatic performance tuning,” in _International Conference on Computational Science_. Springer, 2001, pp. 117–126.
  * [5] J. A. Gunnels, _A systematic approach to the design and analysis of linear algebra algorithms_. The University of Texas at Austin, 2001.
  * [6] V. Valsalam and A. Skjellum, “A framework for high-performance matrix multiplication based on hierarchical abstractions, algorithms and optimized low-level kernels,” _Concurrency and Computation: Practice and Experience_ , vol. 14, no. 10, pp. 805–839, 2002.
  * [7] J. Demmel, J. Dongarra, V. Eijkhout, E. Fuentes, A. Petitet, R. Vuduc, R. C. Whaley, and K. Yelick, “Self-adapting linear algebra algorithms and software,” _Proceedings of the IEEE_ , vol. 93, no. 2, pp. 293–312, 2005.
  * [8] J. R. Herrero Zaragoza, _A framework for efficient execution of matrix computations_. Universitat Politècnica de Catalunya, 2006.
  * [9] J. Cuenca, D. Giménez, and J. González, “Towards the design of an automatically tuned linear algebra library,” in _10th Euromicro Workshop on Parallel, Distributed and Network-based Processing (EUROMICRO-PDP 2002)_. IEEE Computer Society, 2002, pp. 0201–0201.
  * [10] J. Cuenca, D. Giménez, and J. González, “Architecture of an automatically tuned linear algebra library,” _Parallel Computing_ , vol. 30, no. 2, pp. 187–210, feb 2004.
  * [11] T. Katagiri, S. Ohshima, and M. Matsumoto, “Auto-tuning on numa and many-core environments with an fdm code,” in _2017 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)_ , 2017, pp. 1399–1407.
  * [12] S. Catalán, J. R. Herrero, E. S. Quintana-Ortí, R. Rodríguez-Sánchez, and R. Van De Geijn, “A case for malleable thread-level linear algebra libraries: The lu factorization with partial pivoting,” _IEEE access_ , vol. 7, pp. 17 617–17 633, 2019.
  * [13] J. Cámara, J. Cuenca, L.-P. García, and D. Giménez, “Empirical modelling of linear algebra shared-memory routines,” _Procedia Computer Science_ , vol. 18, pp. 110–119, 2013.
  * [14] E. Peise and P. Bientinesi, “Performance modeling for dense linear algebra,” in _2012 SC Companion: High Performance Computing, Networking Storage and Analysis_. IEEE, 2012, pp. 406–416.
  * [15] Y. Xia, M. De La Pierre, A. S. Barnard, and G. M. J. Barca, “A machine learning approach towards runtime optimisation of matrix multiplication,” in _2023 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_ , 2023, pp. 524–534.
  * [16] Y.-C. Lin, Y. Chen, S. Gobriel, N. Jain, G. K. Jha, and V. Prasanna, “Argo: An auto-tuning runtime system for scalable gnn training on multi-core processor,” _arXiv preprint arXiv:2402.03671_ , 2024.
  * [17] M. P. Deisenroth, A. A. Faisal, and C. S. Ong, _Mathematics for machine learning_. Cambridge University Press, 2020.
  * [18] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [19] C. M. Bishop and M. E. Tipping, “Bayesian regression and classification,” _Nato Science Series sub Series III Computer And Systems Sciences_ , vol. 190, pp. 267–288, 2003.
  * [20] C. M. Bishop and N. M. Nasrabadi, _Pattern recognition and machine learning_. Springer, 2006, vol. 4, no. 4.
  * [21] M. Awad and R. Khanna, “Support vector regression,” in _Efficient learning machines_. Springer, 2015, pp. 67–80.
  * [22] N. Bhatia _et al._ , “Survey of nearest neighbor techniques,” _arXiv preprint arXiv:1007.0085_ , 2010.
  * [23] J. Han, J. Pei, and H. Tong, _Data mining: concepts and techniques_. Morgan kaufmann, 2022.
  * [24] M. M. Breunig, H.-P. Kriegel, R. T. Ng, and J. Sander, “Lof: Identifying density-based local outliers,” _SIGMOD Rec._ , vol. 29, no. 2, p. 93–104, may 2000.
  * [25] R. M. Sakia, “The box-cox transformation technique: a review,” _Journal of the Royal Statistical Society: Series D (The Statistician)_ , vol. 41, no. 2, pp. 169–178, 1992.
  * [26] S. Weisberg, “Yeo-johnson power transformations,” _Department of Applied Statistics, University of Minnesota. Retrieved June_ , vol. 1, p. 2003, 2001.
  * [27] M. Mascagni and H. Chi, “On the scrambled halton sequence,” vol. 10, no. 3-4, pp. 435–442, 2004.
  * [28] A. Géron, _Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow: Concepts, tools, and techniques to build intelligent systems_. ” O’Reilly Media, Inc.”, 2019.


# A Demonstration of Benchmarking Time Series Management Systems in the Cloud

Prabhav Arora  New York University Abu DhabiP.O. Box 1212Abu DhabiUAE

###### Abstract.

Time Series Management Systems (TSMS) are Database Management Systems that
have been configured with the primary objective of processing and storing time
series data. With the IoT expanding at exponential rates and there becoming
increasingly more time series data to process and analyze, several TSMS have
been proposed and are used in practice. Each system has its own architecture
and storage mechanisms and factors such as the dimensionality of the dataset
or the nature of the operators a user wishes to execute can cause differences
in system performance. This makes it highly challenging for practitioners to
determine the most optimal TSMS for their use case. To remedy this several
TSMS benchmarks have been proposed, yet these benchmarks focus primary on
simple and supported operators, largely disregarding the advanced analytical
operators (ie. Normalization, Clustering, etc) that constitute a large part of
the use cases in practice. In this demo, we introduce a new benchmark that
enables users to evaluate the performance of four prominent TSMS (TimescaleDB,
MonetDB, ExtremeDB, Kairos-H2) in their handling of over 13 advanced
analytical operators. In a simple and interactive manner, users can specify
the TSMS(s) to compare, the advanced analytical operator(s) to execute, and
the dataset(s) to utilize for the comparison. Users can choose from over eight
real-world datasets with varying dimensions or upload their own dataset. The
tool then provides a report and recommendation of the most optimal TSMS for
the parameters chosen.

PVLDB Reference Format:
PVLDB, 14(1): XXX-XXX, 2020.
doi:XX.XX/XXX.XX ††This work is licensed under the Creative Commons BY-NC-ND
4.0 International License. Visit https://creativecommons.org/licenses/by-nc-
nd/4.0/ to view a copy of this license. For any use beyond those covered by
this license, obtain permission by emailing<EMAIL_ADDRESS>Copyright is held
by the owner/author(s). Publication rights licensed to the VLDB Endowment.
Proceedings of the VLDB Endowment, Vol. 14, No. 1 ISSN 2150-8097.
doi:XX.XX/XXX.XX

PVLDB Artifact Availability:
The source code, data, and/or other artifacts have been made available at
https://github.com/prabhav2302/ABench-IoT-Demo.

##  1\. Introduction

Time series data is essentially a collection of measurements of a certain
quantity obtained at different moments in time. With the rapid expansion of
the Internet of Things (IoT), such data has become increasingly more available
and prevalent in many domains of industry and research. A Time Series
Management System (TSMS) is a Database Management System that has been
specifically formulated and optimized for processing, storing, and interacting
with time series data. Each TSMS has its own architecture and storage
mechanism which determines how data is represented and stored internally
within the system. TSMS architects make several design decisions that in turn
impact the overall performance of their system.

One primary decision is whether to have a row-oriented (storing entries of a
row in contiguous memory locations) or column-oriented (storing entries of a
column in contiguous memory locations) system. Another decision designers must
make is how to structure related time series. Some TSMS store each time series
independently whereas other group subsets or all related time series together
in storage. Though these decisions typically create differences between the
architecture of the systems, one common property of most TSMS is that they
partition data into chunks along the time axis. This representation results in
faster computation times as the data can be pruned more efficiently since the
work is parallelized across the time intervals.

In addition to the system design, each TSMS must also provide support for
various operators that allow users to interact with the data. The complexity
of these operators is highly dependent on the system architecture and design
as they dictate how each operator must be run. Certain architectures, such as
row-store systems, favor certain operations, such as inserting and reading,
whereas other architectures, such as column-store systems, favor other
operations, such as aggregation calculations. Aside from the operators a user
wishes to utilize, the dataset a user wishes to interact with also impacts the
performance of the TSMS. Once again certain TSMS architectures are more
tailored towards datasets of certain properties (granularity, anomalies,
cardinality).

These permutations and combinations of desired user operators and datasets all
result in different TSMS performance and thereby mitigate the existence of a
blanket, silver bullet, “best” TSMS. How then does a practitioner navigate the
performance of these TSMS and choose the most ideal TSMS for their use case?

Multiple benchmarks have been devised that contrast the runtime, memory usage,
architecture, or quality of various TSMS. These evaluations typically begin
with inserting a certain time series dataset into the respective TSMS and then
then querying the system using a fixed set of operations. These benchmarks,
however, largely concern themselves with applying simple and supported
operators, and fail to take into consideration the advanced analytical
operators which constitute a large portion of the use cases in time series
applications. The performance of TSMS executing operators that do Reading and
Inserting have been widely studied whereas the performance of executing
operators that do Normalization, Clustering, and Mapping has been scarcely
documented.

AIoT-Bench aims to be an interactive and easy-to-use tool that enables users
to examine the performance of several TSMS whilst executing various analytical
operators on a diverse set of datasets. Users can evaluate the performance of
eight prominent TSMS (InfluxDB, TimescaleDB, MonetDB, ExtremeDB, Graphite,
Kairos-Cassandra, Kairos-H2, Druid) in their handling of over 9 advanced
analytical operators (Normalize, Map, Decompose, Recover, Repair, Cluster,
Classify, Pattern, Similarity) and 4 simple analytical operators (Records,
SumRecords, MovingAverage, Distance). Users have access to 6 real-world
datasets, each with its own unique and distinct properties, that have been
specifically selected to emulate and encompass the properties of data observed
in practice. Users may also upload their own dataset to compare TSMS
performance for their specific use case. After a user specifies exactly which
TSMS to compare, which analytical operator to execute, and which dataset to
utilize, AIoT-Bench provides a report and recommendation of the most optimal
TSMS for the parameters chosen.

##  2\. Analytical Queries

Our benchmark supports over 13 analytical operators, which can be further
subdivided into 5 overarching types of operators. Below is an explanation of
each operator and its potential use cases.

###  2.1. Advanced Transformation Operators

Transformation operators, are those operators which when given data in one
state perform some transformations and return the data in another, more useful
or insightful state. Our benchmark supports Centroid Decomposition,
SaxRepresentation, ZNormalization, and ZNormalization-Operators.
Centroid Decomposition is an iterative algorithm that decomposes an input
matrix X into a product of two or more matrices. CD can be used to extract the
principal features of a time series.
SaxRepresentation first normalizes the data using ZNormalization, then reduces
its dimensionality using Piece-wise Aggregate Approximation. Then, values are
assigned to a set of pre-defined buckets and mapped to a lower-case letter
based on which bucket they belong to.
ZNormalization first calculates the standard deviation and mean of the data.
Then for each point the Z-Score is calculated. Z-Score for $x$ = ($x$ \- Mean)
/ Standard Deviation. Z-score normalization is used in most time series
analytical tasks as a preprocessing step. Our benchmark tests Z-normalization
using both UDF and the in-built operators of each database.

###  2.2. Machine Learning Operators

Machine Learning operators, are operators which perform a ML algorithm on a
dataset. Our benchmark supports KMeans and KNN.
KMeans is an unsupervised learning algorithm which iteratively clusters the
dataset into $n$ groups and retrieves $n$ centroids.
KNN is a supervised learning algorithm which classifies each point in the
dataset to a category from a finite set based on its euclidean distance to the
nearest neighbour from each category.

###  2.3. Erroneous Data Handling Operators

Erroneous Data Handling Operators are operators which deal with missing or
anomalous data. Our benchmark supports Screen, Recovery, and HotSax operators.
Screen is a technique that detects anomalies and repairs them.
Recovery is implemented using RecovDB(Arous et al., 2019). The input data has
missing values and output data estimates the missing values and adds them to
the initial set.
HotSax(Keogh et al., 2005) processes each column at a time and looks for
anomalies by first normalizing the data and then applying euclidean distance
and a SAX representation. It returns a number $x$ of anomalies for time series
within a given time interval

###  2.4. Simple Operators

In order to provide a comprehensive and an all-inclusive benchmark, Our
benchmark also supports basic operators such as Select and Sum.
Select the datapoints (observations) whose timestamps lie in a given window.
Provides a way for the clients to retrieve the data within the desired time-
range.
Sum computes an aggregate over the entire time series using the SUM operator.

###  2.5. Similarity Detection Operators

Similarity Detection operators are operators which identify or quantify the
degree of similarity between time series. Our benchmark supports Distance and
DSTree.
Distance computes the Eucledian distance between the time corresponding points
in time series $ts$1 and $ts$2\. Is the building block for more complex
operators like KNN, and KMeans.
DSTree(Leung and Khan, 2006) retrieves for each time series in set $ts$2, the
closest time series in set $ts$1\. (Closeness defined by Euclidean distance
between two time series). The algorithm is split into two parts: the indexing
(using DS-tree) and the actual search based on the index.

Figure 1. Architecture Diagram

##  3\. Benchmark GUI

Our benchmark has a graphical user interface implemented with the PyQt5
framework for the front end, and Flask server framework for the back end. The
architecture of Our benchmark can be visualized in Figure 1. In step 1, a
client interacts with the GUI and selects the parameters and clicks run, which
prompts the GUI to send a HTTP GET request to the Flask server. The remote
machine which the Flask server is hosted on has all the prerequisite databases
installed and running. In step 2 the Flask server handles the request and
executes the relevant scripts on the Databases as specified by the client’s
parameters. In step 3, each TSMS sends its runtime results back to the Flask
server. The Flask server then sends the GUI a CSV file containing the results
in step 4. Finally, the GUI parses the CSV file and plots a graph with the
results.

Figure 2. Our benchmark PyQt5-based GUI

###  3.1. Supported Operators

Our benchmark implements all the 13 operators discussed above, enabling users
to test the performance of the various databases on a wide spectrum of use
cases. Each operator script contains a UDF to execute the specified operator.
For Our benchmark all UDFs were coded using Python to ensure consistency and
fairness in experiments.

###  3.2. Datasets

In addition to the user being able to upload custom datasets, Our benchmark
comes preloaded with 8 datasets. Each preloaded dataset exhibits different
characteristics of data (dimensionality, anomalies, regularity, etc.) from the
other. The decision to include datasets with distinct differences has been
made consciously in order to account for the various use cases of data one can
expect to encounter in real life.

Name | Number of TS | Length of TS | Features
---|---|---|---
Alabama | 46 | 3.500.000 | anomalies, regular
Sports | 360 | 142.500 | regular, anomalies
MexData | 512 | 7.000 | irregular, anomalies
Hydraulic | 43.680 | 2.205 | regular,anomalies

The Alabama weather dataset(for Environmental Information, 2021) has been
included due its large magnitude and large length of time series (3.5m). The
Sports dataset(Dua and Graff, 2017) was chosen as a middle-way dataset, with a
medium number of time series (360) and a medium-large length (140k). Next, the
MexData(Wijekoon et al., 2019) is selected as an irregular dataset of a medium
number of time series and small-medium length. Finally, the Hydraulic
dataset(Chawathe, 2019) was chosen due to its very large number of time
series, small-medium length, and large total volume of data. The other four
preloaded datasets, also display varying data characteristics to those
mentioned above.

Figure 3. Results Plot

##  4\. Demonstration Scenarios

Scenario 1: Choosing A Preloaded Dataset. In this scenario a user must first
select which databases he/she wishes to execute queries on. Then, the user
must decide which operator they wish to use. Next, the user must choose from
one of the 8 pre-loaded datasets available on the benchmark. Then, they must
specify exactly how many rows and columns of the dataset, they want the
benchmark to take into consideration. As soon as the user hits the ”Run”
button, a HTTP GET request is sent to the Flask server and after the queries
are executed on the remote machine, the results are returned and plotted for
the user. Figure 3 shows the results of Our benchmark for Centroid
Decomposition on 1000 rows and 10 columns of the Alabama dataset. Subsequent
to retrieving the results of a query, the user may tinker with the parameters
and further execute more tests.
Scenario 2: Uploading a Custom Dataset. The workflow for this scenario is
identical to Scenario 1, except in the part where the user has to select a
Dataset. Instead of choosing one of the 8 preloaded Datasets, a user must
instead choose ”User Dataset” on the Dataset drop down menu. Then, the user
should click the button with 3 dots and choose which dataset file from their
machine they want to use. Subsequent to retrieving the results of a query, the
user may tinker with the parameters and further execute more tests.

###### Acknowledgements.

This work was supported by the eXascale Infolab in Switzerland. We also thank
Professor Djellel Difallah for his mentorship and oversight of this research
project.

## References

  * (1)
  * Arous et al. (2019) Ines Arous, Mourad Khayati, Philippe Cudré-Mauroux, Ying Zhang, Martin Kersten, and Svetlin Stalinlov. 2019.  RecovDB: Accurate and Efficient Missing Blocks Recovery for Large Time Series. In _2019 IEEE 35th International Conference on Data Engineering (ICDE)_. 1976–1979.  https://doi.org/10.1109/ICDE.2019.00218
  * Chawathe (2019) Sudarshan S. Chawathe. 2019\.  Condition Monitoring of Hydraulic Systems by Classifying Sensor Data Streams. In _2019 IEEE 9th Annual Computing and Communication Workshop and Conference (CCWC)_. 0898–0904.  https://doi.org/10.1109/CCWC.2019.8666564
  * Dua and Graff (2017) Dheeru Dua and Casey Graff. 2017.  UCI Machine Learning Repository.  http://archive.ics.uci.edu/ml
  * for Environmental Information (2021) National Centers for Environmental Information. 2021.  Alabama Near Coastal Meteorological & Hydrographic Continuous Data Sampling from 2003 to present.  https://doi.org/10.25921/3edp-9d76
  * Keogh et al. (2005) E. Keogh, J. Lin, and A. Fu. 2005.  HOT SAX: efficiently finding the most unusual time series subsequence. In _Fifth IEEE International Conference on Data Mining (ICDM’05)_. 8 pp.–.  https://doi.org/10.1109/ICDM.2005.79
  * Leung and Khan (2006) Carson Kai-sang Leung and Quamrul I. Khan. 2006.  DSTree: A Tree Structure for the Mining of Frequent Sets from Data Streams. In _Sixth International Conference on Data Mining (ICDM’06)_. 928–932.  https://doi.org/10.1109/ICDM.2006.62
  * Wijekoon et al. (2019) Anjana Wijekoon, Nirmalie Wiratunga, and Kay Cooper. 2019\.  MEx: Multi-modal Exercises Dataset for Human Activity Recognition.  https://arxiv.org/abs/1908.08992