11institutetext: Maryam Parvizi, Amirreza Khodadadian, Sven Beuchler, Thomas
Wick 22institutetext: Leibniz University Hannover, Institute of Applied
Mathematics, Welfengarten 1, 30167 Hannover, Germany 22email:
<EMAIL_ADDRESS>33institutetext: and Cluster of Excellence PhoenixD (Photonics, Optics, and
Engineering - Innovation Across Disciplines), Leibniz Universität Hannover,
Germany

# Hierarchical LU preconditioning for the time-harmonic Maxwell equations

Maryam Parvizi  Amirreza Khodadadian  Sven Beuchler  Thomas Wick

###### Abstract

The time-harmonic Maxwell equations are used to study the effect of electric
and magnetic fields on each other. Although the linear systems resulting from
solving this system using FEMs are sparse, direct solvers cannot reach the
linear complexity. In fact, due to the indefinite system matrix, iterative
solvers suffer from slow convergence. In this work, we study the effect of
using the inverse of $\mathcal{H}$-matrix approximations of the Galerkin
matrices arising from Nédélec’s edge FEM discretization to solve the linear
system directly. We also investigate the impact of applying an
$\mathcal{H}-LU$ factorization as a preconditioner and we study the number of
iterations to solve the linear system using iterative solvers.

##  1 Introduction

The Maxwell system describes the behaviour of electromagnetic fields.
Nédélec’s edge finite element method is an efficient discretization technique
to solve this equation which involves the curl-curl operator numerically
monk2003finite ; hiptmair2002finite . Although the matrix arising from the
linear system is sparse, direct solvers cannot solve the problem with linear
complexity. In the presence of quasi-uniform meshes, the memory requirement of
$O(N^{4/3})$ and computational time of $O(N^{2})$ are expected where the
problem size is $N$ liu1992multifrontal .

Matrices with full rank often can be approximated using low-rank matrices;
but, it is not always applicable. Thus, it is desirable to present a block-
wise partitioning of the matrix and approximate appropriately chosen (using an
admissibility condition) blocks by their low-rank decompositions. Hierarchical
matrices ($\mathcal{H}$-matrices), hackbusch2015hierarchical , are block wise
low-rank matrices that allow us to represent dense matrices with data sparse
approximations and the logarithmic-linear storage complexity, i.e.,
$\mathcal{O}(Nm\log(N))$ where $m$ is a parameter that controls the accuracy
of the approximation.

Besides the storage complexity reduction, another application of the
$\mathcal{H}$-matrix approximations is to use them as a preconditioner to
solve the system directly, or to use them as preconditioner to reduce the
number of iterations in Krylov-based iterative solvers based on matrix-vector
multiplication, e.g., GMRES saad2003iterative . In the time-harmonic case, the
system matrix may become indefinite and ill-conditioned, in particular for
high frequency cases. In this regime, the usual factorization methods such as
incomplete LU (iLU) do not lead to reliable results and converge to the exact
solution poorly. Then it is very difficult to design Galerkin discretizations
MeSau11  and efficient iterative solvers ErnstGander  (see also
HENNEKING202185  for recent studies).

In this paper, we study the effect of applying hierarchical $LU$
decompositions, i.e., ${\mathcal{H}}-LU$ decompositions, as a preconditioner
to solve the Maxwell equations using iterative solvers. The idea of using
$\mathcal{H}$ matrices for the curl-curl operators and magnetostatic model
problems was introduced in bebendorf2009parallel . A preconditioner based on
the $\mathcal{H}$ matrices was used in bebendorf2012hierarchical  and
ostrowski2010cal  for solving the Maxwell equations in the low-frequency
regime. In faustmann2022h , the authors showed that the inverse of the
Galerkin matrix corresponding to the FEM discretization of time-harmonic
Maxwell equations can be approximated by an $\mathcal{H}$-matrix and proved
the exponential convergence in the maximum block-rank. This approximation can
be used to prove the existence of $\mathcal{H}$-$LU$ factorization without
frequency restriction.

In this work, in addition to addressing the advantage of iterative
hierarchical preconditioning, we exploit the influence of applying an inverse
of $\mathcal{H}$-matrix approximations as a preconditioner to solve the linear
systems directly. Our numerical tests also include studies on the influence of
the wave number. As observed, for both low and high frequency materials using
${\mathcal{H}}-LU$ factorization will lead to a fast and accurate convergence
of the iterative solvers.

The paper is organized as follows. In Section 2, we briefly introduce the
Maxwell system and present the Nédélec’s finite element discretization. The
$\mathcal{H}$-matrices and how to compute the inverse of an
$\mathcal{H}$-matrix approximation for the Galerkin system matrix are
explained in Section 3. We also present an algorithm to compute the
$\mathcal{H}-LU$ approximation of a matrix, and use it as a preconditioner to
solve the system directly. In Section 4, numerical results are presented to
substantiate the efficiency of the hierarchical matrix as a direct and
iterative solver. Finally, the conclusions are given in Section 5.

##  2 The Maxwell equations

Denoting $\hat{\mathcal{H}}$ the magnetic field intensity, $\mathcal{E}$ the
electric field density for the domain $\Omega\subset\mathbb{R}^{d}\,(d=2,3)$,
we have the time-harmonic Maxwell equations as

$\displaystyle\nabla\cdot(\beta\hat{\mathcal{H}})$ $\displaystyle=0\quad\hskip
28.45274pt$ $\displaystyle\text{in}\,\,\Omega\times\mathcal{I},$ (1a)
$\displaystyle\nabla\cdot(\alpha\mathcal{E})$ $\displaystyle=\rho\quad\hskip
28.45274pt$ $\displaystyle\text{in}\,\,\Omega\times\mathcal{I},$ (1b)
$\displaystyle\left(\alpha\frac{\partial}{\partial
t}+{\color[rgb]{0,0,0}\chi}\right)\mathcal{E}{\color[rgb]{0,0,0}-}\nabla\times\hat{\mathcal{H}}$
$\displaystyle=\mathcal{F}\quad\hskip 28.45274pt$
$\displaystyle\text{in}\,\,\Omega\times\mathcal{I},$ (1c)
$\displaystyle\beta\frac{\partial}{\partial
t}\hat{\mathcal{H}}+\nabla\times\mathcal{E}$ $\displaystyle=0\quad\hskip
28.45274pt$ $\displaystyle\text{in}\,\,\Omega\times\mathcal{I},$ (1d)

where $\mathcal{F}$ is the applied electrical force, $\rho$ is the charge
density, and $\mathcal{I}=(0,T]$ is the time interval. In addition, $\alpha$
and $\beta$ are the dielectric and magnetic permeabilities, and
${\color[rgb]{0,0,0}\chi}$ is the conductivity constant. Considering an
arbitrary frequency $\omega$, with respect to time, the electric and magnetic
fields can be represented as follows

$\displaystyle\mathcal{E}(x,t)$ $\displaystyle=e^{-i\omega
t}\boldsymbol{E}(x),$ (2a) $\displaystyle\hat{\mathcal{H}}(x,t)$
$\displaystyle=e^{-i\omega t}\boldsymbol{H}(x),$ (2b)
$\displaystyle\mathcal{F}(x,t)$ $\displaystyle=e^{-i\omega
t}\boldsymbol{F}(x).$ (2c)

Replacing (2a) and (2b) into (1c) and (1d), we obtain

$\displaystyle-\nabla\times\boldsymbol{H}-i\omega\zeta\boldsymbol{E}$
$\displaystyle=\textbf{F}(x)\quad\hskip 28.45274pt$
$\displaystyle\text{in}\,\,\Omega,$ (3a)
$\displaystyle\nabla\times\boldsymbol{E}-i\omega\beta\boldsymbol{H}$
$\displaystyle=0\quad\hskip 28.45274pt$ $\displaystyle\text{in}\,\,\Omega,$
(3b)

where
$\zeta:={\color[rgb]{0,0,0}\alpha}+{i{\color[rgb]{0,0,0}\chi}}/{\omega}$.
Also, we consider a perfect conduction boundary condition (we surround
$\Omega$ by a perfect bounded material i.e., $E\times\boldsymbol{n}=0$).
Therefore, we have the following second-order operator for (3)

$\displaystyle{\boldsymbol{\mathcal{L}}}\boldsymbol{E}:=\nabla\times(\beta^{-1}\nabla\times\boldsymbol{E})-\kappa\boldsymbol{E}=\boldsymbol{J}_{S}\qquad\text{in}\,\,\Omega,$
(4)

where $\kappa:=\omega^{2}\zeta$ and
$\boldsymbol{J}_{S}:=-i\omega\boldsymbol{F}$.

###  2.1 Discretization by edge elements

To present a Galerkin weak formulation for (4),  we denote by
$\boldsymbol{L}^{2}(\Omega)$ as the space of vector field with three entries
from $L^{2}(\Omega)$, i.e.,

$\boldsymbol{L}^{2}(\Omega):=\left\\{\textbf{U}=(U_{1},U_{2},U_{3})\quad:\quad
U_{i}\in L^{2}(\Omega),\,\,i=1,2,3\right\\},$

with $\left<\cdot\,,\cdot\right>$ as the inner product on this space, and
continue with defining the following space

$\displaystyle\boldsymbol{H}(\operatorname*{curl},\Omega):=\left\\{\mathbf{U}\in\boldsymbol{L}^{2}(\Omega)\;\colon\;\nabla\times\mathbf{U}\in\boldsymbol{L}^{2}(\Omega)\right\\},$

equipped with the norm

$\left\|\textbf{U}\right\|^{2}_{\boldsymbol{H}(\operatorname*{curl},\Omega)}:=\left\|\textbf{U}\right\|^{2}_{\boldsymbol{L}^{2}(\Omega)}+\left\|\nabla\times\textbf{U}\right\|^{2}_{\boldsymbol{L}^{2}(\Omega)}.$

Considering homogeneous Dirichlet boundary conditions, the space
$\boldsymbol{H}_{0}(\operatorname*{curl},\Omega)\subset\boldsymbol{H}(\operatorname*{curl},\Omega)$
is introduced as follows

$\boldsymbol{H}_{0}(\operatorname*{curl},\Omega):=\\{\textbf{U}\in\boldsymbol{L}^{2}(\Omega)\;\colon\;\nabla\times\textbf{U}\in\boldsymbol{L}^{2}(\Omega),\;\textbf{U}\times\textbf{n}=0\;\text{on}\;\Gamma\\}.$

Then, the weak formulation for (4) can be written as : find
$\boldsymbol{E}\in\boldsymbol{H}_{0}(\operatorname*{curl},\Omega)$

$\displaystyle
a(\boldsymbol{E},\Phi):=\left<\nabla\times\boldsymbol{E},\nabla\times\Phi\right>-\kappa\left<\boldsymbol{E},\Phi\right>=\left<\boldsymbol{J}_{S},\Phi\right>\qquad\forall\Phi\in\boldsymbol{H}_{0}(\operatorname*{curl},\Omega).$
(5)

Here, we should note that $\kappa$ is not an eigenvalue of the operator
$\nabla\times\nabla\times$ (monk2003finite, , Corollary . 4.19), in addition
we have $\kappa\neq 0$, and we set $\beta=1$. The existence of the unique
solution for the variational formulation (5) is proven in (hiptmair2002finite,
, Thm. 5.2), and the following a priori estimate is obtained

$\left\|\boldsymbol{E}\right\|_{\boldsymbol{H}(\operatorname{curl},\Omega)}\leq
C^{*}\left\|\boldsymbol{J}_{S}\right\|_{L^{2}(\Omega)},$ (6)

where $C^{*}$ depends on $\Omega$ as well as $\kappa$.

For the discretization, we consider quasi-uniform mesh simplices
${\mathcal{T}}=\\{T_{1},\dots,T_{N_{\mathcal{T}}}\\}$ where
$T_{j}\in\mathcal{T}$ are open elements and denote
$h:=\max_{T_{j}\in\mathcal{T}}{\rm diam}(T_{j})$. We assume $\mathcal{T}$ is a
Ciarlet-regular mesh, i.e., it does not contain any hanging nodes. We also
assume there is $\gamma>0$ such that ${\rm
diam}(T_{\ell})\leq\gamma\,|T_{\ell}|^{1/3}$ for all $T_{\ell}\in\mathcal{T}$.
In order to present the Galerkin FEM for (5), we consider lowest order
Nédélec’s $\operatorname*{curl}$-conforming elements of the first kind, i.e.,

$\displaystyle\textbf{V}_{h}$
$\displaystyle:=\\{\mathbf{v}_{h}\in\boldsymbol{H}(\operatorname*{curl},\Omega)\;\colon\;\mathbf{v}_{h}|_{T}\in\mathcal{N}_{0}(T)\quad\forall
T\in\mathcal{T}\\},$ $\displaystyle\textbf{V}_{h,0}$
$\displaystyle:=\textbf{V}_{h}\cap\boldsymbol{H}_{0}(\operatorname{curl},\Omega),$

where for all $T\in\mathcal{T}_{h}$, the lowest order local Nédélec space of
the first kind $\mathcal{N}_{0}(T)$ is defined as monk2003finite

$\displaystyle{\color[rgb]{0,0,0}\mathcal{N}_{o}(T)=\\{\boldsymbol{a}+\boldsymbol{b}\times\mathbf{x}\;\colon\;\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{3},\quad\mathbf{x}\in
T\\}.}$

We denote $\mathcal{Y}_{h}:=\\{\Phi_{1},\ldots,\Phi_{N}\\}$ as a basis with
$N$ as the dimension of $\textbf{V}_{h,0}$. This basis is uniquely defined by
the property $\sigma(\Phi_{i},e_{j})=\delta_{ij}$, where $e_{j}$ denotes an
interior edge of the mesh and $\sigma(p,e)$ is the line integral of the
tangential component of $p$ along $e$. Then, the Galerkin FEM for (5) is given
as: Find $\boldsymbol{E}_{h}\in\textbf{V}_{h,0}$ such that

$\displaystyle
a(\boldsymbol{E}_{h},\Phi_{h})=\left<\boldsymbol{J}_{S},\Phi_{h}\right>\qquad\forall\Phi_{h}\in\textbf{V}_{h,0}.$
(7)

This is equivalent to solve the following system

$\displaystyle\textbf{A}x=b\quad\textrm{where}\quad\textbf{A}=[\textbf{A}_{ij}]_{i,j=1}^{N}\quad\textrm{with}\quad\textbf{A}_{ij}:=a(\Phi_{i},\Phi_{j}),\qquad\Phi_{j},\,\Phi_{i}\in\mathcal{Y}_{h},$
(8)

and the right hand side vector $b$ is defined as
$b_{j}:=\left<\boldsymbol{J}_{S},\Phi_{j}\right>$, $j\in\\{1,2,\ldots,N\\}$.

##  3 $\mathcal{H}$-matrices and $\mathcal{H}$-matrix arithmetic

$\mathcal{H}$-matrices are defined based on a partition $P$ generated from a
clustering algorithm that allows us to determine which blocks can be
approximated by low-rank matrices or are small hackbusch2015hierarchical .

Applying $\mathcal{H}$-matrix approximations allows us to store large matrices
in the format of low-rank block-wise matrices which could lead to logarithmic-
linear storage complexity provided that a proper method is used to define the
hierarchical structure that results in the final block-wise format of the
matrix. In the following lemma from faustmann2022h , it is shown that the
inverse of the Galerkin matrix A (8) can be approximated using an
$\mathcal{H}$-matrix, and proven that this approximation converges
exponentially with respect to the maximum block rank to A.
Definition 1 [$\mathcal{H}$-matrices] A matrix
$\mathbf{B}_{\mathcal{H}}\in\mathbb{C}^{N\times N}$ is called an
$\mathcal{H}$-matrix, if for every admissible block $(\tau,\sigma)$, we have
the following factorization

$\mathbf{B}_{\mathcal{H}}|_{\tau\times\sigma}=\mathbf{X}_{\tau\sigma}\mathbf{Y}^{H}_{\tau\sigma},$

of rank $r$ where $\mathbf{X}_{\tau\sigma}\in\mathbb{C}^{\tau\times r}$ and
$\mathbf{Y}_{\tau\sigma}\in\mathbb{C}^{\sigma\times r}$.
In order to use the inverse of the $\mathcal{H}$-matrix approximation of A as
a preconditioner, first, we need to find an $\mathcal{H}$-matrix approximation
for A, then we obtain the inverse using the iterative method of Schulz
hackbusch2000sparse .

###### Lemma 1

faustmann2022h  Let A be the Galerkin matrix defined in (8). Then, there
exists an $\mathcal{H}$-matrix approximation $\mathbf{B}_{\mathcal{H}}$ with
the maximum block rank $r$ (rank of all the blocks of
$\mathbf{B}_{\mathcal{H}}$ is either smaller than or equal to $r$) such that

$\displaystyle\left\|\mathbf{A}^{-1}-\mathbf{B}_{\mathcal{H}}\right\|_{2}\leq\bar{C}h^{-1}e^{-c(r^{1/4}/\ln
r)},$

where $\bar{C}$ and $c$ are constants depending only on material parameters,
and the properties of $\Omega$.

In the definition of $\mathcal{H}$-matrices, the low-rank blocks are
determined based on the concept of $\eta$-admissibility defined in
hackbusch2015hierarchical . In the following, the mathematical definition of
an $\mathcal{H}$-matrix is given.

Although computing the inverse of the $\mathcal{H}$-matrix approximation of
the Galerkin system matrix leads to logarithmic-linear complexity, the
computational cost to solve the linear system directly is still too high.
Thus, we need to use another alternative to reduce the numerical cost such as
the $\mathcal{H}-LU$ factorization. In the following, we present an algorithm
from (MR2451321, , Sec. 2.9), and use it as a preconditioner to solve the
linear systems.

1) Compute the $\mathcal{H}$-matrix approximation of A, i.e.,
$\textbf{A}_{\mathcal{H}}$
2) Compute the $\mathcal{H}$\- matrix $LU$ decomposition of
$\textbf{A}_{\mathcal{H}}$ as follows

1) for $j=1,\ldots,N$

for $k=1,\ldots,j-1$

Solve the system
$\displaystyle{\sum_{i=1}^{k}}L_{ji}\,U_{ik}={(\textbf{A}_{\mathcal{H}})}_{jk}$
to get $L_{jk}.$

2) Compute $L_{jj}$ and $U_{jj}$ by
$\displaystyle{L_{jj}U_{jj}=A_{jj}-\sum_{i=1}^{j-1}L_{ji}U_{ij}}$.

3) for $k=j+1,\ldots,N$

Solve the system
$\displaystyle{\sum_{i=1}^{j}}L_{ji}\,U_{ik}={(\textbf{A}_{\mathcal{H}})}_{jk}$
to get $U_{jk}.$

3) for $i=1,\ldots,\text{MaxIt}$

Compute $r=b-\textbf{A}\,x$ and
err=$\left\|r\right\|_{2}/\left\|b\right\|_{2}$

Compute $x=x+U^{-1}\left(L^{-1}\,r\right)$.

if $\texttt{err}<\texttt{TOL}$ break

4) Compute $\texttt{error}=\left\|\textbf{A}\,x-b\right\|_{2}$.

Algorithm 1 $\mathcal{H}-LU$ decomposition and application in preconditioning
of a simple iterative solver for (7).

##  4 Main algorithm and numerical experiments

In this section, we first present the main algorithm of this work, and then we
study three numerical examples to solve the linear system arising from the
Maxwell equations using an $\mathcal{H}$-matrix approximation as a
preconditioner. For this, we employ the geometrically balanced cluster tree
presented in grasedyck2004adaptive , and we set the admissibility parameter
$\eta=2$. We use a truncated singular value decomposition (SVD) with different
ranks $r$ to compute $\textbf{B}_{\mathcal{H}}$ as the inverse of
$\textbf{A}_{\mathcal{H}}$ obtained from the Schulz iteration. In other words,
for admissible blocks $(\tau,\sigma)$, we have
$\textbf{A}_{\mathcal{H}}|_{\tau\times\sigma}:=\textbf{U}_{r}\textbf{S}_{r}\textbf{V}^{T}_{r}$
where $\textbf{U}_{r}\in\mathbb{R}^{\tau\times r}$,
$\textbf{S}_{r}\in\mathbb{R}^{r\times r}$ and
$\textbf{V}_{r}\in\mathbb{R}^{\sigma\times r}$ are the first $r$ columns of U,
S and V, respectively. Then, we find the inverse $\textbf{B}_{\mathcal{H}}$
for the matrix $\textbf{A}_{\mathcal{H}}$.

###  4.1 Example 1: a unit box

The geometry is $\Omega=(0,1)^{3}$ and $\boldsymbol{J}_{S}=[0,0,1]^{\top}$.
The coarse mesh consists of 6 tetrahedrons. This mesh is uniformly refined $k$
times, $k=2,\ldots,7$. All computations are performed in MATLAB with 125
cores. Table 1 displays the iteration numbers of the preconditioned GMRES
method with the described $\mathcal{H}$-matrix preconditioner for different
values of $\kappa$. The GMRES method is stopped if a relative accuracy of
${\color[rgb]{0,0,0}TOL}=10^{-5}$ of the residual is reached. In all
experiments, the parameter $\beta=1$ is chosen. After 100 iterations, we
restart GMRES. The results show that the $\mathcal{H}$-matrices can be used as
an efficient preconditioner if $\kappa\leq 20^{2}$. For higher frequencies,
the iteration numbers grow in some, but not all levels $k$. This is due to the
computation of the $LU$ decomposition of the $\mathcal{H}$-matrix. The
approximation of the original matrix by the $\mathcal{H}$-matrix is still very
good, also in the case of $k=5$ and $\kappa=900$.

$\kappa~{}~{}$ | $N_{\text{dof}}$ | $~{}~{}25~{}~{}$ | $~{}~{}100~{}~{}$ | $~{}~{}225~{}~{}$ | $~{}~{}400~{}~{}$ | $~{}~{}625~{}~{}$ | $~{}~{}900~{}~{}$
---|---|---|---|---|---|---|---
$k=2$ | $98$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$
$k=3$ | $604$ | $2$ | $2$ | $2$ | $2$ | $2$ | $2$
$k=4$ | $3\,184$ | $2$ | $4$ | $5$ | $5$ | $4$ | $3$
$k=5$ | $41\,024$ | $3$ | $4$ | $8$ | $9$ | $11$ | $>3000$
$k=6$ | $238\,688$ | $3$ | $6$ | $20$ | $4$ | $48$ | $25$
$k=7$ | $1\,807\,264$ | $4$ | $5$ | $7$ | $19$ | $46$ | $93$
Table 1: GMRES iterations numbers for Example 1.

###  4.2 Example 2: two boxes

Here the geometry consists of two boxes, i.e.,
$\Omega:=(-1,1)\times(-1,1)\times[-1,-2)\cup(-2,2)\times[1,2)\times(-1,1)\cup(-2,2)\times[-1,1]\times[-1,1]\cup(-1,1)\times(-1,1)\times[1,2)\cup(-2,2)\times[-1,2)\times(-1,1)$.
We set the parameters $\kappa=1$, $\boldsymbol{J}_{S}=[1,1,1]^{\top}$, and
$\beta=1$. The computational domain with and the inverse of
$\mathcal{H}$-matrix approximation of the stiffness matrix is shown in Figure
1. We have 24 440 admissible blocks, 48 404 small blocks and the depth is 16.
The decay of the approximation error versus $r$ and the corresponding
allocated memory are shown in Figure 2. As shown, using higher $r$ gives rise
to a reliable inverse of the $\mathcal{H}$-matrix approximation. The computed
$\textbf{B}_{\mathcal{H}}$ can be used as a preconditioner to solve the linear
system (7) directly.

Figure 1: Example 2. The computational geometry (left) and the inverse of
$\mathcal{H}$-matrix approximation of the stiffness matrix (7)
$\mathcal{H}$-matrix clustering (right) for $N_{\text{dof}}=22\,001$. Figure
2: Approximation error of $B_{\mathcal{H}}$ in Example 2 for a stiffness
matrix with $N_{\text{dof}}=22\,001$ (left) and the relative allocated memory
(right).

###  4.3 Example 3: a magnet surrounded by air

We consider a magnet surrounded by air where the box is $1\times 1\times 1$
and the magnet dimension is $0.5\times 0.5\times 0.75$. We set $\kappa=10$,
$\beta=10$ and $\boldsymbol{J}_{S}=[10,10,10]^{\top}$. The geometry and the
$\mathcal{H}$-approximation of A for $N_{\text{dof}}=122\,202$ is shown in
Figure 3. In this approximation, we have 215 964 admissible blocks, 402 451
small blocks, and the depth is 15. The ${\mathcal{H}}-$LU decomposition of A
is given in Figure 4. We use Algorithm 1 for different $N_{\text{dof}}$ to
solve the linear system. Table 2 shows the results for different matrices
using $TOL=1\times 10^{-8}$. For all cases, the negligible error indicates the
accuracy the method, and the elapsed CPU time points out its efficiency. For
the last two examples, we used Netgen/ngsolve schoberl2017netgen  to produce
the meshes (denoting different $N_{\text{dof}}$) and assembling the matrices
of (7).

Figure 3: Example 3. The computational geometry (left) and $\mathcal{H}$-matrix approximation of A for $N_{\text{dof}}=122\,202$ (right). Figure 4: The $\mathcal{H}$-LU decomposition for the stiffness matrix resulting from Example 3 corresponding to Figure 3. Table 2: The iterative solver preconditioned by $\mathcal{H}-LU$ used to solve the Maxwell equations for different $N_{\text{dof}}$. $N_{\text{dof}}$ |  5 492 |  8 050 |  13 602 |  33 933 |  48 811 |  70 133 |  78 603 |  96 846 |  129 200 |  304 309
---|---|---|---|---|---|---|---|---|---|---
error |  7.22 e -9 |  1.07 e -9 |  5.80e-9 |  2.23e-8 |  5.11e-8 |  9.06e-8 |  5.98e-10 |  4.38e-10 |  5.04e-10 |  5.78e-9
time [s] | 11.52 | 35.24 | 60.52 | 151.71 | 256.18 | 546.9 | 602.73 | 481 | 820.9 | 2335
iterations | 3 | 3 | 3 | 3 | 4 | 4 | 5 | 5 | 6 | 7

##  5 Conclusion

In this work, the $\mathcal{H}$-matrix approximation was used to solve the
time-harmonic Maxwell equations. As a direct solver, the inverse of the
hierarchical matrix approximation of the linear system was employed as a
preconditioner, where an accurate approximation was achieved. Additionally, we
then employed an $\mathcal{H}-LU$ factorization as a preconditioner.

In both cases, the use of $\mathcal{H}$ matrix approximations could reduce the
computational cost and increase the accuracy of the solution. The
$\mathcal{H}$ matrices can be coupled with the domain decomposition to take
advantage of both approaches, i.e., to reduce the complexity and accelerate
the convergence of the iterative solvers. This possibility will be addressed
in the future papers.

###### Acknowledgements.

The authors acknowledge the Deutsche Forschungsgemeinschaft (DFG) under
Germany Excellence Strategy within the Cluster of Excellence PhoenixD (EXC
2122, Project ID 390833453). Maryam Parvizi is funded by Alexander von
Humboldt Foundation project named $\mathcal{H}$-matrix approximability of the
inverses for FEM, BEM and FEM-BEM coupling of the electromagnetic problems.
Finally, the authors thank Sebastian Kinnewig for fruitful discussions.

## References

  * [1] M. Bebendorf.  Hierarchical matrices.  Springer, 2008.
  * [2] M. Bebendorf and F. Kramer.  Hierarchical matrix preconditioning for low-frequency–full-maxwell simulations.  Proceedings of the IEEE, 101(2):423–433, 2012.
  * [3] M. Bebendorf and J. Ostrowski.  Parallel hierarchical matrix preconditioners for the curl-curl operator.  Journal of Computational Mathematics, pages 624–641, 2009.
  * [4] O. G. Ernst and M. J. Gander.  Why it is difficult to solve Helmholtz problems with classical iterative methods.  pages 325–363. Springer, 2012.
  * [5] M. Faustmann, J. M. Melenk, and M. Parvizi.  $\mathcal{H}$-matrix approximability of inverses of FEM matrices for the time-harmonic Maxwell equations.  Advances in Computational Mathematics, 48(5), 2022.
  * [6] L. Grasedyck, W. Hackbusch, and S. L. Borne.  Adaptive geometrically balanced clustering of $\mathcal{H}$-matrices.  Computing, 73(1):1–23, 2004.
  * [7] W. Hackbusch.  Hierarchical matrices: algorithms and analysis, volume 49.  Springer, 2015.
  * [8] W. Hackbusch and B. N. Khoromskij.  A sparse $\mathcal{H}-$matrix arithmetic: general complexity estimates.  Journal of Computational and Applied Mathematics, 125(1-2):479–501, 2000.
  * [9] S. Henneking and L. Demkowicz.  A numerical study of the pollution error and dpg adaptivity for long waveguide simulations.  Computers & Mathematics with Applications, 95:85–100, 2021.
  * [10] R. Hiptmair.  Finite elements in computational electromagnetism.  Acta Numerica, 11:237–339, 2002.
  * [11] J. W. Liu.  The multifrontal method for sparse matrix solution: Theory and practice.  SIAM Review, 34(1):82–109, 1992.
  * [12] J. M. Melenk and S. Sauter.  Wavenumber explicit convergence analysis for galerkin discretizations of the helmholtz equation.  SIAM Journal on Numerical Analysis, 49(3):1210–1243, 2011.
  * [13] P. Monk et al.  Finite element methods for Maxwell’s equations.  Oxford University Press, 2003.
  * [14] J. M. Ostrowski, M. Bebendorf, R. Hiptmair, and F. Krämer.  $\mathcal{H}$-matrix-based operator preconditioning for full maxwell at low frequencies.  IEEE Transactions on Magnetics, 46(8):3193–3196, 2010.
  * [15] Y. Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [16] J. Schöberl et al.  Netgen/ngsolve.  https://ngsolve.org, 2017.


# DRPC: Distributed Reinforcement Learning Approach for Scalable Resource
Provisioning in Container-based Clusters

Haoyu Bai , Masaki Owari M. Hayashi is with Graduate School of Mathematics,
Nagoya University, Nagoya, JapanM. Owari is with the Faculty of Informatics,
Shizuoka University, Hamamatsu, Shizuoka, Japan. Haoyu Bai, Minxian Xu,
Kejiang Ye, Rajkumar Buyya, Chengzhong Xu H. Bai and R. Buyya are with the
School of Computing and Information Systems, the University of Melbourne,
Melbourne, Australia. H. Bai contributed this work when he was a visiting
student at Shenzhen Institute of Advanced Technology, Chinese Academy of
Sciences. M. Xu and K. Ye are with Shenzhen Institute of Advanced Technology,
Chinese Academy of Sciences, Shenzhen, China. C. Xu is with State Key Lab of
IOTSC, University of Macau, Macau, China. This work is supported by the
National Key R&D Program of China (No.2021YFB3300200), the National Natural
Science Foundation of China (No. 62072451, 62102408), Shenzhen Industrial
Application Projects of undertaking the National key R & D Program of China
(No. CJGJZD20210408091600002), and Shenzhen Science and Technology Program
(No. RCBS20210609104609044). M. Xu is the corresponding author.

###### Abstract

Microservices have transformed monolithic applications into lightweight, self-
contained, and isolated application components, establishing themselves as a
dominant paradigm for application development and deployment in public clouds
such as Google and Alibaba. Autoscaling emerges as an efficient strategy for
managing resources allocated to microservices’ replicas. However, the dynamic
and intricate dependencies within microservice chains present challenges to
the effective management of scaled microservices. Additionally, the
centralized autoscaling approach can encounter scalability issues, especially
in the management of large-scale microservice-based clusters. To address these
challenges and enhance scalability, we propose an innovative distributed
resource provisioning approach for microservices based on the Twin Delayed
Deep Deterministic Policy Gradient algorithm. This approach enables effective
autoscaling decisions and decentralizes responsibilities from a central node
to distributed nodes. Comparative results with state-of-the-art approaches,
obtained from a realistic testbed and traces, indicate that our approach
reduces the average response time by 15% and the number of failed requests by
24%, validating improved scalability as the number of requests increases.

###### Index Terms:

Cloud computing, Kubernetes, microservice, reinforcement learning, distributed
resources management

##  I introduction

Microservice architecture has emerged as a transformative approach in the
field of software design and development. It is characterized by its modular
and decentralized structure, where complex applications are broken down into
smaller, independently deployable services[1][2]. Each microservice focuses on
a specific business capability, allowing teams to work on distinct components
simultaneously and enabling rapid development, deployment, and scalability.
This architecture promotes flexibility, resilience, and maintainability by
minimizing the impact of changes within one service on the overall system [3].
As microservices communicate through well-defined APIs, they facilitate
seamless integration and support heterogeneous technology stacks. The cloud
service providers, such as Amazon, Google, Microsoft, and Alibaba, have
embraced microservice to develop and deploy their applications in cloud
computing environment [4].

The rapid growth of the microservices paradigm has introduced challenges in
resource management [5]. As the number of microservices increases within a
system, ensuring optimal resource allocation and utilization becomes complex.
Microservices autoscaling has emerged as a viable solution to address these
challenges [6]. Autoscaling involves dynamically adjusting the number of
instances or resources allocated to microservices based on real-time demand,
ensuring efficient resource utilization while maintaining desired performance
levels [7]. This technique encompasses both horizontal scaling, which involves
adding or removing replicas of a service, and vertical scaling, which involves
adjusting the resources allocated to each instance. Autoscaling mechanisms
utilize various metrics and triggers, such as CPU usage, memory consumption,
and request rates, to make decisions about scaling actions [8]. Implementing
effective autoscaling strategies enhances system responsiveness, minimizes
operational costs, and optimizes resource allocation in dynamic and
unpredictable environments.

Moreover, the intricate interdependencies among microservices often give rise
to critical paths and key nodes that significantly impact performance.
Efficiently implementing autoscaling for microservices is confronted with its
own set of challenges and intricacies. The dynamic nature of these
dependencies necessitates a profound understanding of the application’s
behavior, workload patterns, and their implications for system performance.
Identifying the optimal scaling strategy, striking a balance in resource
allocation across interconnected services while maintaining overall system
stability, proves to be a non-trivial undertaking [9]. Additionally, the real-
time nature of scaling decisions and the imperative to minimize disruptions to
service quality and user experience further contribute to the complexity.

Reinforcement learning (RL) [10] has the potential to optimize the scaling of
microservices, taking into account their intricate interdependencies.
Nevertheless, the majority of RL solutions adopt a centralized decision-making
approach, which poses scalability challenges. This centralized structure
encounters difficulties with the increasing complexity and interdependencies
of microservices, resulting in adverse effects on scalability, response times,
and reliability. Centralized management approaches such as Borg [11] permit
users to over-provision resources while assigning jobs to machines, leading to
resource wastage and consequently diminishing overall performance in resource-
limited environments. In addition, centralized approaches can compromise
system performance as microservices scale up, given the continuous
communication between huge amount of services and the central node may
introduce bottlenecks and latency.

In this paper, we explore novel approaches to address these limitations and
propose a distributed reinforcement learning framework for microservices
automatic scaling. This framework aims to enhance scalability, improve
decision-making efficiency, and ensure the adaptability of microservices
scaling strategies to dynamic and evolving environments. By decentralizing the
decision-making process, we aim to mitigate the challenges associated with
centralized approaches and unlock the full potential of RL in optimizing
microservices resource allocation and scalability.

The main contributions of this paper can be summarized as follows:

  * •

We design a distributed reinforcement learning framework for resource
provisioning for container-based autoscaling (DRPC). This framework
facilitates precise modelling of system resources and their scalable
allocation to address the dynamic demands of microservices.

  * •

We propose a method for selecting optimization strategies based on
reinforcement learning and a distributed algorithm that incorporates domain
knowledge through deep imitation learning. This approach facilitates efficient
and adaptive decision-making, optimizing the choice of scaling strategies
across clusters of microservices.

  * •

We perform comprehensive testing and validation of our proposed model and
policies, utilizing real-world traces and a dedicated testing platform. These
experiments are conducted to thoroughly assess the effectiveness and
performance of our distributed reinforcement learning framework for
microservices autoscaling.

##  II Related Work

In this section, we discuss the current autoscaling methods for microservices,
categorizing them into three groups according to their key mechanisms:
threshold-based and heuristic approaches, machine learning (ML) based
approaches, and deep learning (DL) based approaches.

###  II-A Threshold-based and Heuristic Autoscaling

The threshold-based heuristic approach for microservice resource allocation
relies on predefined rules, scaling resources up when utilization surpasses a
predetermined threshold (e.g., 75%). This method proves efficient in scenarios
with abundant resources and stable request patterns. He et al. [12] leveraged
both genetic and heuristic algorithms to determine the optimized microservice
deployment location within an edge-cloud environment. Horizontal pod auto-
scaler (HPA) [13] is a scaling technique adopted in Kubernetes, primarily
focusing on horizontal scaling. It is capable of dynamically adjusting the
number of replicas based on resource variations, such as CPU and memory,
within the existing servers. The primary objective of HPA is to determine the
number of replicas to be added or removed based on system running status.
Kannan et al. [14] conceptualized multi-stage tasks using a Directed Acyclic
Graph (DAG). By leveraging the DAG, they could predict the task’s completion
duration. Their method continually adapts thresholds and algorithms for each
microservice. As a result, every microservice can manage its loads
autonomously, incurring minimal communication expenses.

###  II-B Machine Learning based Autoscaling

The ML-based autoscaling for microservices dynamically adjusts resources
through ML algorithms that analyze historical data and workload patterns. This
method optimizes resource utilization and ensures consistent performance by
scaling services as required. Hou et al. [15] proposed an autoscaling approach
that emphasizes both power and latency considerations for resource
provisioning at micro and macro levels. By employing decision trees and a
tagging methodology, they were able to expedite the resource-matching process.
Yu et al. [16] introduced a framework called Microscaler that utilizes a
service mesh to monitor resource usage patterns. By integrating online
learning and heuristic methods, the framework achieves near-optimal solutions
to satisfy resource needs and quality of service (QoS) requirements. Liu et
al. [17] undertook a detailed analysis of bottlenecks in prevalent
microservice applications and incorporated various machine learning models to
facilitate resource scheduling. Meanwhile, Gan et al. [18] harnessed
predictive methodologies to detect QoS violations. They leveraged a
combination of machine learning models and expansive historical data to
pinpoint the microservices responsible for these QoS infractions, enabling
better resource reallocation to mitigate QoS degradation.

TABLE 1: Comparison of related work

Approach | Scaling Techniques | Workload Prediction | Scheduling Mechanism | Decision-making Pattern
---|---|---|---|---
Vertical | Horizontal | Brownout | Linear | Non-Linear | Heuristic | Machine Learning | Deep Learning | Distributed | Centralised
He et al. [12] | $\checkmark$ | $\checkmark$ |  |  |  | $\checkmark$ |  |  | $\checkmark$ |
Burns et al. [13] |  | $\checkmark$ |  |  |  | $\checkmark$ |  |  |  | $\checkmark$
Kannan et al. [14] | $\checkmark$ |  |  | $\checkmark$ |  | $\checkmark$ |  |  |  | $\checkmark$
Hou et al. [15] | $\checkmark$ |  |  | $\checkmark$ |  |  | $\checkmark$ |  |  | $\checkmark$
Yu et al. [16] |  | $\checkmark$ |  |  | $\checkmark$ |  | $\checkmark$ |  |  | $\checkmark$
Liu et al. [17] | $\checkmark$ |  |  |  | $\checkmark$ |  | $\checkmark$ |  |  | $\checkmark$
Gan et al. [18] | $\checkmark$ |  |  |  | $\checkmark$ |  | $\checkmark$ |  |  | $\checkmark$
Rzadca et al. [19] | $\checkmark$ | $\checkmark$ |  | $\checkmark$ |  |  |  | $\checkmark$ |  | $\checkmark$
Xu et al. [20] | $\checkmark$ | $\checkmark$ | $\checkmark$ |  | $\checkmark$ |  |  | $\checkmark$ |  | $\checkmark$
Qiu et al. [21] | $\checkmark$ |  |  |  |  |  |  | $\checkmark$ |  | $\checkmark$
Zhang et al. [22] |  | $\checkmark$ |  |  |  |  |  | $\checkmark$ |  | $\checkmark$
Rossi et al. [23] | $\checkmark$ | $\checkmark$ |  | $\checkmark$ |  |  |  | $\checkmark$ |  | $\checkmark$
Our method (DRPC) | $\checkmark$ | $\checkmark$ | $\checkmark$ |  | $\checkmark$ |  |  | $\checkmark$ | $\checkmark$ |

###  II-C Deep Learning based Autoscaling

The emphasis has shifted towards employing DL for microservice autoscaling,
utilizing neural networks for decision-making and pattern analysis. This
approach dynamically adjusts resource allocation based on real-time demand and
optimizes provisioning by learning from historical and current data, thereby
minimizing resource wastage. Autopilot [19] accumulates historical server data
and then utilizes two scheduling techniques: the sliding window algorithm
based on past data, and the meta-algorithm inspired by RL. CoScal [20]
classifies resource usage into four representative levels. It then applies an
approximated Q-learning algorithm to these segments, establishing an estimated
policy. This policy encompasses horizontal and vertical scaling, as well as
brownout [24] capabilities that can dynamically activate and deactivate
application component. With gated recurrent unit (GRU), CoScal predicts
upcoming workloads and refers to the trained lookup table to determine the
suitable scaling action for each pod, ensuring optimization across all pods.
FIRM [21] employs a systematic approach to allocating resources in cloud
systems. It pinpoints the crucial path in the microservice dependency by
carefully examining the connections between components. Once this path is
identified, FIRM employs a specialized Support Vector Machine (SVM) tailored
to operate on both a per-critical-path and per-micro-service-instance basis.
This helps in identifying specific microservice instances that require
optimization. Zhang et al. [22] introduced a predictive RL algorithm for
horizontal container scaling, which combines the Autoregressive Integrated
Moving Average (ARIMA) model with a neural network model, ensures the
predictability and precision of the scaling procedure. Rossi et al. [23]
developed RL-based strategies to manage both horizontal and vertical scaling
for containers. This approach enhances system adaptability in the face of
fluctuating workloads and hastens the learning phase by harnessing varying
levels of environmental knowledge.

###  II-D Critical Analysis

Although the existing representative methods have brought valuable
contributions, our proposed method advances the relevant area in several key
points. First, unlike HPA and other threshold-based approaches that primarily
focus on horizontal scaling, our approach leverages multi-faceted scaling
approaches (horizontal, vertical and brownout that can dynamically activate or
deactivate optional microservices [24]), ensuring optimal resource allocation
even during non-overload-states. Second, compared with the machine learning-
based approach, we have applied deep learning to capture the features of
workloads and utilized RL to make scaling decisions to achieve more accurate
and efficient resource provisioning decisions. The compared differences with
the related work are highlighted in Table 1.

Our approach is most similar to CoScal [20] and FIRM [21] based on RL to auto-
scale microservices while having significant differences compared with them.
In contrast to CoScal’s large-grained segmentation of system states, our
method employs a deep neural network, offering nuanced decision-making and the
capability for simultaneous multiple scaling actions, crucial for sudden load
changes. Compared to FIRM, which focuses on optimizing resources of the
critical path and nodes, our method ensures comprehensive optimization, with
the advantage of supporting both horizontal and vertical scaling, and
brownout. Furthermore, our advanced resource allocation framework provides
precise solutions, overcoming the scalability limitations inherent in FIRM’s
centralized reinforcement learning approach, especially under a high volume of
requests.

##  III Motivation

In this section, we perform motivational experiments with a use case that
requests increase quickly to explore the system scalability under a
centralized design, such as the centralized database capturing the dynamic
change of microservices (FIRM) [21] and the centralized RL-based approach
(CoScal) [20]. The investigation focuses on response time as the number of
requests significantly increases. For these experiments, we utilize the
TrainTicket [25] microservice application comprising approximately 40
services, including user, station, price, and route. The experiments are
conducted on two nodes: one for the initial TrainTicket deployment and another
for scaling, with a configuration of up to 8 CPU cores and 8 GB memory.

Fig. 1 shows the performance of FIRM and CoScal when the number of requests
increases from 200 to 800 per second within a short time (e.g. 1 minute), and
we can notice apparent performance fluctuations. For instance, for the FIRM
approach, the response time increases from 160 ms to 290 ms when requests
increase from 200 to 400 per second. This resulted from the limitation of the
centralized design of FIRM in that its central node contains functionalities
including data collection, training, inference, and service provisioning. The
CoScal based on centralized RL design also suffers from the same issue when
the number of requests increases significantly within a short time, CoScal
might over-provision into the system to ensure response time but also incurs
resource wastage. For instance, several replicas are added to provision
resources when the number of requests increases and the response time is also
reduced.

Figure 1: Response time of centralized approaches when the number of requests
increases significantly

In this work, our objective is to provision resources efficiently via an RL-
based approach, and utilize a distributed framework to overcome the limitation
of centralized design, e.g. system performance bottleneck. A centralized
module responsible for data collection and inference can exhibit drawbacks
such as unreliability under a container-based environment, slow performance,
and inability to adjust system resources asynchronously. The system could also
potentially fail to adjust resources if the central nodes experience
performance degradation. Hence, we propose a distributed framework utilizing
multiple lightweight neural networks on distributed nodes to execute the
operations sent from the central node to relieve the burden of the central
node. This approach can potentially hasten resource allocation, and provide a
more accurate resource usage prediction of the cloud system’s behaviour.
Moreover, it could facilitate differential resource adjustment, accommodating
services with rapid changes more frequently than others. In the subsequent
sections, we will introduce our detailed solution.

##  IV System model

In this section, we will introduce the system model of DRPC, which adheres to
the Monitor-Analyze-Plan-Execute framework over a shared Knowledge pool
(MAPE-K), as depicted in Fig. 2. The model is composed of three key
components: (1) a Workload Processor and Predictor that preprocess the raw
workloads and predicts the future workloads, (2) a Central Teacher Network
that aims to learn the globally optimal policy, and (3) a Per-
Deployment111Please note that the term Per-Deployment is inherited the naming
conventions module of Kubernetes, which consists a set of pods to run an
application workload, usually does not need to maintain state. Distributed
Student Network responsible for data collection as well as asynchronous for
imitation learning.

Figure 2: System Model of DRPC

###  IV-A Workload Processor and Predictor

The Workload Processor is composed of key components including a Workload
Preprocessor, a Load Generator (e.g. Locust [26]), and a Historical Database
containing historical workloads. It functions to extract necessary workload
data and attributes from the Historical Database, preprocesses the dataset via
Workload Preprocessor, and handles realistic user interactions with the Load
Generator. It transforms these interactions into requests allocated to the
microservice-based cluster. The Historical Database stores the historical data
derived from realistic traces (e.g. Alibaba traces[27]). Detailed dataset
information such as timestamps, machine identifiers and different types of
resources are contained in these workloads.

The Workload Predictor is a critical component, responsible for forecasting
future loads coming into the system. It communicates to the auto-scaler about
the required amount of resources to be scaled-in or scaled-out. The Workload
Predictor obtains preprocessed workload data from the Workload Processor and
employs predictive models such as Long Short-Term Memory (LSTM)[28] or Gated
Recurrent Unit (GRU) to forecast future workloads. The responsibilities of the
Workload Predictor include managing the training process of the model,
updating the trained model when necessary, and deploying the trained models in
their final stages. This module interacts with the Workload Analyzer and
Workload Generator, receiving historical system data as input.

###  IV-B Central Teacher Network

The use of multi-agent RL often achieves sub-optimal solutions [29],
indicating a need for a central agent (“teacher”) module that suggests actions
to another agent (“student”) [30]. This module investigates the global state,
formulating an optimal policy function to enhance system performance. The
collected data is employed to train the decentralized student network via the
technique of imitation learning [31].

Existing performance modelling-based or heuristic-based strategies suffer
several limitations, including struggles with model reconstruction and
retraining due to their inability to adapt to dynamic system statuses [32].
These strategies require considerable expert knowledge [33], demanding
extensive efforts for the interpretation of microservice workloads and
infrastructure. Compared with them, RL is suitable for formulating resource
provisioning policies due to its capabilities of offering a close feedback
loop, exploring scaling actions, and formulating optimal policies independent
of erroneous assumptions. This allows for learning directly from actual
various workloads and operational conditions, providing deeper insights into
how changes in low-level resources influence application performance. We have
incorporated two techniques below to achieve our objective:

####  IV-B1 Twin Delayed Deep Deterministic Policy Gradient (TD3) [34]

To overcome the limitation of overestimation of $Q$ value due to the flexible
cloud environment, we utilize the TD3 algorithm, an advanced model-free,
actor-critic RL framework. Our RL formulation provides two distinct
advantages:

  * •

The model-free RL eliminates the need for precise modelling of the complete
distribution of states or the environment dynamics (transitions between
states). In the event of microservice updates, the simulations of state
transitions used in model-based RL become inefficient.

  * •

The actor-critic framework, merging policy-based and value-based methods, is
suitable for continuous stochastic environments. This accelerates convergence,
reduces variance, and provides a robust and efficient solution for managing
the dynamic nature of container-based environments.

####  IV-B2 Retraining Notifier

The Retraining Notifier signals when the system requires retraining. There are
instances where the imitation learner might not make efficient decisions,
triggering the Retraining Notifier and initiating a retraining state. This is
typically initiated by the following scenarios:

  * •

Insufficient exploration by the teacher network: If the teacher network
encounters a rarely seen state about which it lacks confidence, it might fail
to guide the student appropriately. This can potentially result in the cloud
service scaler misjudging the system’s scaling needs.

  * •

Sub-optimal learning state of the student deployment network: In this case,
despite the central network learning an optimal policy function, the
distributed deployment network fails to establish a correlation between the
optimal policy and the domain information it possesses. This can occur if the
information provided by the central teacher network is over-complex, requiring
an adjustment in the number of epochs that the distributed deployment net runs
to achieve convergence.

  * •

Outdated knowledge held by the student deployment network: For instance, if
the majority of users previously searched for tickets, the student net may
over-rely on this pattern. If the request composition changes, such as users
predominantly purchasing tickets, the deployment of student nets would require
retraining.

###  IV-C Distributed Student Deployment Network

The Distributed Student Deployment Network module operates in a distributed
fashion, enabling the parallel execution of scaling actions for microservices.
Rather than maintaining global states, it relies on local information to
function. This module consists of two key components: the Deployment Buffer,
which stores the policies communicated by the teacher node for training the
student network, and the Imitation Learner, which scales each deployment
independently.

####  IV-C1 Deployment Buffer

It maintains the policy transmitted by the Central Module, and stores the
state of the deployment when the policy is recorded. This policy-state pairing
is utilized to train the Deployment Student Network via imitation learning,
which aims to ensure replicas’ successful behaviours in a given context.

####  IV-C2 Imitation Learner

The Distributed Student Deployment Net is characterized by its lightweight
nature and lower resource requirements compared to the Central Teacher
Network. This characteristic facilitates more frequent and adaptable decision-
making. The system operates in two distinct modes: learning mode, where it
receives guidance from the Central Network, and acting mode, where it
autonomously makes decisions based on its state, including factors like
resource utilization. This dual-mode configuration enables continuous learning
and adaptation, leading to optimized resource allocation and enhanced system
performance.

##  V Distributed RL Algorithm for Scaling Microservice

In this section, we introduce the algorithm design of our proposed DRPC
approach, which can achieve efficient distributed decisions.

In DRPC, the RL-based resources provisioning is modelled as a Markov Decision
Process, it views the state $s_{t}\in S$ as the current microservices system
status and interprets the action $a\in A$ as a scaling operation modifying
system status and allocated resources. We use $M=(m^{1},m^{2},$ …, $m^{i})$ to
denote the physical machines provisioning resources for the microservices in
our system. For each physical machine, say $m^{i}$, the amount of resources
that can be allocated are denoted as
$R\textsuperscript{i}=(r^{i,1},r^{i,2},\ldots,r^{i,j})$, where the $j$
represents the type of resources, such as CPU and memory. Finally, the set of
actions that can be performed on each physical machine is denoted as
$A^{i}=(a^{i,1},a^{i,2},\ldots,a^{i,k})$. The actions indicate the amount of
resources that can be altered on each machine. The supported actions are
Horizontal_scaling, CPU_scaling, Memory_scaling, and Brownout. A positive or
negative sign before an action implies the addition or reduction of resources
to a specific machine, respectively. The term ’Brownout’ is a binary indicator
for whether a brownout can be triggered on the machine. If this value is set
to True, the minimum allowable replicas for horizontal scaling would be set to
0. The collective action space for the system,
$A=\prod_{i=1}^{I}\prod_{k=1}^{k}A^{i,k}$, is the product of the action spaces
for each machine.

###  V-A System-wide Action Execution

Input : A central module for controlling scaling actions

Output : None

1

2Initialize $trainingMode$ as True or False based on the system state;

3

4while  _True_ do

5 if  _$trainingMode$_ then

6 // Stage 1: The central module explores the state and chooses actions
$actions$ = central module.detection();

7 centralModule.Scale($actions$) based on Alg. 2;

8 foreach  _deploymentAction in getActionByDeployment( $actions$)_ do

9 // Training deployment networks deployment.train(deploymentAction);

10

11 end foreach

12

13 else

14 // Stage 2: Deployments make independent scaling decisions foreach
_deployment in allDeployments_ do

15 deployment.Scale() based on Alg. 2;

16

17 end foreach

18

19 end if

20 // Check if the retraining notifier is triggered and update $trainingMode$
$trainingMode$ = retrainingNotifier();

21

22 end while

Algorithm 1 DRPC: System-wide Action Execution

As illustrated in Algorithm 1, the execution of the system-wide action plan
progresses through two distinct stages: an exploration phase and a distributed
provision phase.

In Stage 1, the exploration phase, (lines 3-9), the central teacher node
explores the state space, choosing multi-dimensional scaling actions
(horizontal, CPU, memory scaling, or maintaining the current state) for each
deployment as described in the section V-C. Simultaneously, deployment-level
networks initiate training.

In Stage 2, the distributed provision phase, (lines 10-12), decision-making
transitions to deployment-level networks once the central node’s exploration
is deemed sufficient, with the same scaling options. This phase continually
updates the teacher network’s buffer with action-state pairs, as described in
section V-E. The system alternates between these stages based on the
retraining notifier (line 14), enabling the central node to re-explore states
or delegate decision-making based on its readiness.

Algorithm 2 illustrates the general scaling procedure of DRPC. The algorithm
seeks advice from the Scaling Agent to obtain Q-values that guide resource
adjustments (line 1). This includes CPU usage (lines 2-4), and memory
utilization (lines 5-7). Such guidance manifests as actions within the DRPC
framework. If the recommended scaling for CPU or memory is significant,
adjustments are made proportionally to ensure system stability or enhance
efficiency. Adding or removing replicas (lines 9-13) follows a similar process
as above, except if a microservice is pre-configured to support brownout, the
minimum number of replicas can be set to 0.

###  V-B Workload Processor and Predictor

The connection between different levels of resource usage and workloads is
determined using a deep neural network model for the workload analysis, as
shown in Fig. 2. To mimic realistic resource utilization, we apply data from
Alibaba traces, which contain workload traces from 4000 machines, encompassing
eight days of resource usage data. The performance profiling procedure is as
follows: we define a scheduling interval of 5 minutes, over which we gradually
increase the number of requests, with each test case comprising 200 requests
sent to the host. A three-layer Multi-Layer Perceptron (MLP) is applied to the
profiling data, thereby efficiently representing this relationship.
Consequently, we can translate the host utilization into the number of
requests to our system for any given level of utilization. Then, our workload
generator produces these requests in accordance with the current user type
composition.

We employ multivariate time series forecasting (MTFS), which converts MTFS
problem into a supervised learning task. The Gated Recurrent Unit (GRU) [35],
a type of recurrent neural network (RNN), is then used as the prediction
method, as it has been validated with good performance in [36] for MTFS. The
GRU solves the vanishing gradient problem encountered in conventional RNNs
through the utilization of gating mechanisms.

TABLE 2: Workloads Prediction Accuracy

Predicted length (mins) | 5 | 10 | 15 | 20 | 25
---|---|---|---|---|---
MSE | 0.0025 | 0.0027 | 0.0033 | 0.0049 | 0.0057

Table 2 shows the mean square errors (MSE) of actual utilization versus
predicted utilization for Alibaba workloads over different predicted length.
With MSE values ranging between 0.002 and 0.006, it shows that the workload
prediction algorithm can achieve good performance in resource utilization
prediction.

###  V-C Central Module

Input: Scaling Agent: $Agent$, CPU usage: $Cpu\\_usage$, Memory usage:
$memory\\_usage$, number of replicas: $replicas$, step size for CPU:
$Cpu\\_step$, step size for memory: $Memory\\_step$, Microservice state:
$thisMicroservice.Brownout$

Output: Updated $Cpu\\_usage$, $memory\\_usage$, $replicas$

1 [$Cpu\\_Scaling$, $Memory\\_Scaling$, $Horizontal\\_Scaling$] $\leftarrow
Agent.getQvalue()$

2

3if  _$|Cpu\\_Scaling| >0.5$_ then

4 $Cpu\\_usage\leftarrow Cpu\\_usage+Cpu\\_Scaling.Clip(-1,1)\times
Cpu\\_step$

5 end if

6if  _$|Memory\\_Scaling| >0.5$_ then

7 $memory\\_usage\leftarrow memory\\_usage+Memory\\_Scaling.Clip(-1,1)\times
Memory\\_step$

8 end if

9if  _$|Horizontal\\_Scaling| >0.5$_ then

10 if  _$replicas+\lceil Horizontal\\_Scaling\rceil <0$ and
$thisMicroservice.Brownout==False$_ then

11 exit

12 else

13 $replicas\leftarrow replicas+\lceil Horizontal\\_Scaling\rceil$

14 end if

15

16 end if

Algorithm 2 DRPC: General Scaling Procedure

Conventional RL models such as Sarsa [37] and Q-learning encounter challenges
when dealing with infinite or continuous states and actions, a common scenario
in microservices characterized by variable CPU and memory usage. To address
this issue, we employ the TD3 algorithm, a more sophisticated approach well-
suited for continuous, high-dimensional spaces. TD3 utilizes twin critic
networks to reduce bias, incorporates delayed policy updates for stable
learning, and employs target policy smoothing to prevent overfitting. These
attributes render TD3 robust and efficient, making it an ideal choice for
optimizing autoscaling in dynamic microservices environments.

####  V-C1 TD3 modelling

Our dual objectives are improving the Quality of Service (QoS) and maximizing
the utilization of physical machines. Our reward model, encompassing response
time $R_{qos}(rt)$ and resource utilization $R_{util}(u)$ are accordingly
formulated in Equations (1) and (2):

$R_{qos}(rt)=\begin{cases}e^{-}({\frac{rt-
RT_{max}}{RT_{max}})^{2}}&,rt>RT_{max}\\\ 1&,rt\leq RT_{max},\end{cases}$ (1)

where the maximum tolerant latency, $RT_{max}$, is pre-defined into the QoS
reward function. Normal system operation rewards 1, while performance that
exceeds $RT_{max}$ is penalized, gradually approaching 0, thus discouraging
SLO violation.

Resource utilization $R_{util}(u)$, measured using a devised model shown in
Equation (2),

$R_{util}(u)=\begin{cases}\frac{\sum_{k=1}^{K}\sum_{r=1}^{R}(U_{r,k}^{pred}-u_{r,k})^{3}}{K}+1&,u_{rk}\leq
U_{r,k}^{pred}\\\
\frac{\sum_{k=1}^{K}\sum_{r=1}^{R}(u_{r,k}-U_{r,k}^{pred})^{3}}{K}+1&,u_{r,k}>U_{r,k}^{pred},\end{cases}$
(2)

which guides the system towards resource conservation. Here, $K$ and $R$
represent the $kth$ physical machine and $rth$ resource type residual on the
$kth$ physical machine, respectively, for example, the $u_{1,1}$ indicate the
1st adjustable resources on the 1st machine. $U_{r,k}^{pred}$, the predefined
utility, represents the ideal utility for a certain resource type $r$ on
machine $k$. Proximity to this ideal utilization is rewarded; under-
provisioning or wastage is discouraged.

The final reward value as shown in Equation (3) is a combination of both
response time and resource utilization, The objective is:

$r(s_{t},a_{t})=\frac{R_{qos}(rt)}{R_{util}(u)}.$ (3)

The final objective is to decrease the response time while keeping the system
running in a stable state $U^{pred}$ as shown in Equation (4).

$\min_{u\in Total\\_Resources}|U^{pred}-R_{\text{util}}(u)|\land\min_{rt\in
RT}R_{qos}(rt).$ (4)

####  V-C2 TD3 Implementation Details

Result: TD3 (Twin Delayed Deep Deterministic policy gradient)

1

2Initialize actor network $\pi$ and critic networks $Q_{\theta_{1}}$,
$Q_{\theta_{2}}$ with random parameters $\theta,\theta_{1},\theta_{2}$

3 Initialize target networks $\pi^{\prime}$ and $Q^{\prime}$ with weights
$\pi^{\prime}\leftarrow\pi,Q^{\prime}_{\theta_{1^{\prime}}}\leftarrow
Q_{\theta_{1}},Q^{\prime}_{\theta_{2^{\prime}}}\leftarrow Q_{\theta_{2}}$

4 Initialize replay buffer $R$

5 for  _episode = 1 to M_ do

6 Observe state $s$

7 for  _t = 1 to T_ do

8 Select action $a=\text{clip}(\pi(s)+\epsilon,a_{High},a_{Low})\text{ where
}\epsilon\sim\mathcal{N}$ and execute it

9 Observe next state $s^{\prime}$, reward $r$ and done signal $d$ to end
episode

10 Store $(s,a,r,s^{\prime},d)$ in $R$

11 if  _it is time to update_ then

12 Sample a batch B of transitions $(s,a,r,s^{\prime},d)$ from $R$

13
$a^{\prime}\leftarrow\text{clip}(\pi^{\prime}(s^{\prime})+\text{clip}(\epsilon,-c,c),a_{High},a_{Low})$

14 $y\leftarrow
r+\gamma(1-d)\min_{i=1,2}Q^{\prime}_{\theta_{i^{\prime}}}(s^{\prime},a^{\prime})$

15 Update $Q_{\theta_{i}}$ by minimizing the loss:
$L(\theta_{i})=\frac{1}{|B|}\sum_{(s,a,r,s^{\prime},d)\in
B}(Q_{\theta_{i}}(s,a)-y)^{2}$

16 if  _$t\mod\text{policy\\_update}==0$_ then

17 Update $\pi$ by one step gradient ascent using the loss:
$L(\theta)=\frac{1}{|B|}\sum_{s\in B}-Q_{\theta_{1}}(s,\pi(s))$

18 Soft update the target networks:
$\theta_{i^{\prime}}\leftarrow\rho\theta_{i^{\prime}}+(1-\rho)\theta_{i}$ (for
both actor and critic)

19

20 end if

21

22 end if

23 $s\leftarrow s^{\prime}$

24

25 end for

26

27 end for

28

Algorithm 3 DRPC: TD3 Algorithm

As shown in Algorithm 3, TD3 begins by initializing the actor and critic
networks, as well as the replay buffer (lines 1-3). It then selects an action
from the actor network, executes the action, and observes the subsequent
state, reward, and episode termination status (lines 4-8). This tuple
$(s,a,r,s^{\prime},d)$ is stored in the replay buffer (line 9). Once enough
data has accumulated and it is time for an update, TD3 samples a batch of
transitions from the buffer, which includes data from both the central and
distributed modules (lines 10-11). In certain instances, scale-out may
considerably diminish the SLO violation, while scale-in, by circumventing the
communication time, could realize superior response times after several
scaling operations. The robust reward resulting from scale-out might engender
a fragile, or erroneous, ”sharp peak,” potentially obscuring our model’s
capacity to investigate the long-term impact of Scale-in. Furthermore, due to
the dynamic nature of cloud environments, where internet connections
fluctuate, some reallocation actions may reduce QoS violations more
effectively than anticipated. Ultimately, this can lead to an overestimation
of Q values. To alleviate this, a clipped (limits values to a range) colour
black noise, between bound $-c\text{ and }c$, $\epsilon$ is incorporated into
the target policy $\pi^{\prime}(s^{\prime})$ to get a final policy
$a^{\prime}$ and that policy is ensured to lie within $a_{Low}\text{ and
}a_{High}$ as shown in Equation (5) (line 12).

$a^{\prime}\leftarrow\text{clip}(\pi^{\prime}(s^{\prime})+\text{clip}(\epsilon,-c,c),a_{High},a_{Low}),$
(5)

where both Q-functions employ a singular target. Moreover, as microservices
are deployed in a distributed manner and communicate via cable, certain
scaling actions may lead to a significant reduction in communication time due
to internet fluctuations, ultimately generating a substantial reward. To
circumvent the overestimation of rewards, Clipped double-Q learning is
employed. Both Q-functions utilize a singular target $y(r,s^{\prime},d)$,
calculated by summing the immediate reward $r$ and the minimum value obtained
from the two Q-functions
$Q^{\prime}_{\theta_{i^{\prime}}}(s^{\prime},a^{\prime})$, multiplied by the
discount factor $\gamma$. If the subsequent state is a terminal state ($d=1$),
no future reward is considered (line 13) as shown in Equation (6):

$y(r,s^{\prime},d)\leftarrow
r+\gamma(1-d)\min_{i=1,2}Q^{\prime}_{\theta_{i^{\prime}}}(s^{\prime},a^{\prime}).$
(6)

The loss indicated as $L(\theta_{i})$ in line 14, for the $i_{th}$ critic
network, is calculated as the average mean of the sum of squared differences
between the target Q-value $y$, and the critic network predicted Q-value
$Q_{\theta_{i}}(s,a)$, across a mini-batch experiences $B$ containing a set of
$(s,a,r,s^{\prime},d)$ experience. The number of such experiences is denoted
by $|B|$ as formulated in Equation (7):

$L(\theta_{i})=\frac{1}{|B|}\sum_{(s,a,r,s^{\prime},d)\in
B}(Q_{\theta_{i}}(s,a)-y)^{2}.$ (7)

To minimize the variability in microservice communication and to learn a more
stable $Q$ function, the policy will only be updated by one step of gradient
ascent once every $policy\\_update$ (line 15) times using the loss function
(line 16) with Equation (8):

$L(\theta)=\frac{1}{|B|}\sum_{s\in B}-Q_{\theta_{1}}(s,\pi(s)).$ (8)

Finally, the target network will be updated (lines 17-20).

###  V-D Deployment Unit with Imitation Learning

Imitation learning involves training an agent to perform tasks by emulating
the actions of an expert. This approach has gained popularity due to its
capability to teach complex behaviors without the need for explicit
programming or heavy reliance on reward signals, which is common in
traditional RL. In our scenario, the deployment unit utilizes its information
to imitate the resource adjustment actions of the Central module. This
objective, as shown in Eq. (9), is accomplished by minimizing the mean square
error between the Q-values $Q_{\text{Ctr}}(s)$ generated by the central
network and those $Q_{\text{Dtr}}(s)$ produced by the distributed network
under the same state $s$.

$\min\left(Q_{\text{Ctr}}(s)-Q_{\text{Dtr}}(s)\right)^{2}$ (9)

When the retraining notifier signals deployment readiness, scaling actions,
guided by human knowledge-driven frequency, are executed simultaneously and
asynchronously. State-action-reward pairs are then gathered, combined with
data from other deployments, and utilized to fine-tune the policy of the
central network.

###  V-E Retraining Notifier

When the RL-based model needs to be re-trained based on the conditions as
introduced in Section IV.B, the Retraining Notifier as shown in Algorithm 4
informs the teacher network to repeat the training phase whenever the student
network faces insufficient information or the average reward from recent
actions falls below a predetermined threshold. It initializes an array,
$rewardHistory$, to record the last $npr$ rewards and the iterator number
$iter$; continuously, the algorithm calculates the current action’s reward,
storing it cyclically in $rewardHistory$ (lines 1-7). If collected rewards are
below the defined number $npr$, retraining will not be triggered due to
insufficient data (line 8). However, if the average of $rewardHistory$ drops
below the threshold $TH$, the retraining process will be triggered, indicating
suboptimal performance (lines 10-12). Otherwise, the system continues its
distributional operations (lines 13-15).

Input : Retraining Threshold $TH$, number of past rewards to monitor $npr$,
action sets $a_{t}=\\{a_{k}^{i}(t)|i\in\\{0,1,\ldots,I\\},$
$k\in\\{1,2,\ldots,K\\}\\}$

Output : Boolean indicating whether retraining mode is triggered

1

2Initialize an array $rewardHistory$ of size $npr$ to store the rewards from
past actions

3

4Initialize an Iterator $iter$ = 0

5

6while  _True_ do

7 Calculate the current action reward $currentReward$ using
get_reward($a_{t}$)

8

9 Store the current reward at the position $iter$ % $npr$ in $rewardHistory$

10

11 $rewardHistory[iter\%npr]$ = $currentReward$

12

13 $iter$ ++

14

15 if  _$iter <npr$_ then

16 // Not enough data gathered yet, do not trigger retraining return False

17

18 if  _average( $rewardHistory$) $<$ $TH$_ then

19 // The average reward is below the threshold, trigger retraining return
True

20

21 else

22 // The average reward is above the threshold, do not trigger retraining
return False

23

24 end if

25

26 end while

Algorithm 4 DRPC: Retraining Notifier

##  VI PERFORMANCE EVALUATIONS

(a) Number of success requests

(b) Number of failure requests

(c) Average response time

(d) Transactions per second

Figure 3: Performance comparison of KuScal, FIRM, CoScal, and DRPC

To assess the efficacy of DRPC in the autoscaling of microservices, we conduct
experiments on a realistic testbed based on Kubernetes. We detail the
experimental settings, benchmarks, and a comprehensive analysis of the
results. The primary objective is to verify that a distributed asynchronous
framework can improve performance, potentially offering potential directions
for future research.

###  VI-A Experimental Setup

We utilize TrainTicket as a testbed in our motivational example in Section
III. It provides functionalities such as ticket purchases, availability
checks, cancellations, and news browsing, simulating backend processes like
verification code generation, user login checks, and database ticket
availability searches. The platform establishes latency-measuring chains using
domain name service for service connectivity. Load balancers or service
routers based on cloud server tests are employed. Each TrainTicket
microservice has its own database, reducing wait times and facilitating
scheduling for developers.

We use a cluster with five machines (one master and four workers) for
microservice-based cluster. Each machine has 8 CPU and 8GB memory. This
prototype system is implemented using a suite of toolkits, including Cgroup
v2, Python, TensorFlow, PyTorch, Sklearn, and Locust, which support deep
learning and reinforcement learning environments.

We simulate four user types (Normal Query Users, Ticket Buyers, Cancel Ticket
Users, and Admin Users) with various ratios in our cluster using the Alibaba
work trace, each assigned with a unique ID for tracking. This emulates real-
world user behavior dynamics, where service requests fluctuate over time. Data
is collected every 200 milliseconds and averaged over 1000-minute intervals in
accordance with the existing work [20][21].

###  VI-B Baselines

Several state-of-the-art approaches from both industry and academia have been
selected as baselines.

KuScal [13] is a Kubernetes mechanism using horizontal scaling to adjust
replicas dynamically. It determines replica quantity based on resource
fluctuations, adjusting microservice numbers via continuous real-time resource
utilization tracking and predefined CPU utilization thresholds. The threshold
is set to be 75% as it is evaluated as the most effective in some existing
work surveyed in [38].

CoScal [20] divides resource usage into quarterlies, applies an approximated
Q-learning algorithm to devise a policy, and supports horizontal, vertical
scaling, and brownout. Using the GRU unit for workload prediction, it
forecasts workload, consulting a trained table to guide pod-level scaling
actions, and even optimizing non-critical path pods.

FIRM [21] systematically allocates cloud resources by identifying the critical
microservice path and using machine learning to target specific instances
needing optimization. The process strategically moves from the longest to the
shortest chain, prioritizing longer chains to reduce microservice time delays
and boost efficiency.

###  VI-C Experiment Analysis

(a) CDF of response time

(b) Scalability comparison

(c) Reward convergence of DRPC

Figure 4: CDF comparison and reward convergence

Several widely used metrics, including the number of success requests, failure
rate and response time have been utilized to evaluate the performance of our
proposed approach and the baselines.

####  VI-C1 Computation Efficiency Analysis

Based on the efficient design of our distributed student network, DRPC
necessitates considerably fewer computational resources compared to FIRM, with
trainable parameters being 643 for our model and 11,352 for FIRM. DRPC
utilizes only $6\times 10^{-5}$ of CPU time to determine the subsequent
action, compared to $4\times 10^{-3}$ taken by FIRM. Our model also showcases
superior efficiency, with a peak CPU utilization of 0.8 as opposed to FIRM’s
maximum usage of 2.0. This illustrates that our approach is markedly more
resource-efficient than FIRM. DRPC incorporates 40 services, and enables
simultaneous application of CPU, memory, and horizontal scaling for each
deployment, executing up to 120 actions per time step. KuScal ranks second,
when CPU utilization for all services surpasses the threshold, and all
services undergo horizontal scaling. However, FIRM can only scale one resource
type per step, as FIRM’s Actor-critic network estimates the Q-value for
scaling actions and applies the action with the highest Q-value.

The offline training phase required approximately 10 days to converge, with
the majority of this time dedicated to the data collection process within the
cluster. However, once the initial training is completed, the model can be
fine-tuned within minutes via online learning to adapt to changes in the
composition of user requests, which can handle the significant changes in
realistic system.

####  VI-C2 Comparison of the Number of Success Requests per Second

Fig. 3(a) presents the requests per second processed by four distinct
algorithms. We partitioned the 5000-minute data into five time periods, with
each period representing the average results over the corresponding period. We
observed that, across the five time periods, FIRM is capable of processing a
higher number of requests compared to our method during the initial periods.
However, our method outperformed FIRM over the subsequent four periods.
Likewise, CoScal also demonstrated the ability to adapt and improve over time,
but it did not outperform our method or FIRM consistently. Furthermore, the
difference between our method, CoScal, and FIRM was relatively small during
the initial periods, until our method started surpassing them. It can be also
reasoned that a scheduling algorithm solely reliant on the longest chain
becomes increasingly challenging to handle workloads under significant
variances in request rate over time. This might be caused by the longest chain
typically signifying more requests at the onset of the scenario when chain
variations are minimal. However, as the load changes over time, the chains
experiencing high loads, which may not necessarily be the longest chains,
become the actual targets for optimization. In summary, FIRM and Coscal have
achieved good performance in handling requests successfully. Compared to FIRM,
DRPC improves the success rate of requests by around 2%.

####  VI-C3 Comparison of Failure Rate

As shown in Fig. 3(b), DRPC has a failure rate of 4.5%, which is notably lower
than FIRM’s 5.9%, CoScal’s 8.2%, and Kuscal’s significantly high rate with
29%. When we evaluate the improvement of DRPC over the other methods, it
reduced the failure requests by 24% compared to FIRM, 44% in relation to
CoScal, and exhibited an impressive 85% reduction when compared to KuScal. In
summary, DRPC demonstrates a significant performance improvement over the
benchmarks. This enhancement is attributed to our design’s ability to adjust
deployments more frequently, thereby reducing the likelihood of failures
compared to other methods.

####  VI-C4 Comparison of Response Time

Fig. 3(c) demonstrates the average response time for the four algorithms under
the same configurations. KuScal shows the longest average response time across
all five periods, being three to four times longer than others, as its
inability to predict forthcoming requests inhibits it from conducting pre-
emptive horizontal scaling to accommodate these incoming requests. In
contrast, CoScal performs better than KuScal but trails behind FIRM and DRPC.
The CoScal model demonstrates an average response time that is generally
similar to that of FIRM, except for the last two time periods, where FIRM
shows a slight gap. We observed that the DRPC achieves the shortest response
time, typically ranging between 40 ms and 70 ms, for most of the time
intervals, thereby reducing the average response time by around 15% compared
to FIRM and 24% compared to Coscal. Remarkably, it reduces the response time
by 72.4% compared to Kuscal.

####  VI-C5 Scalability Analysis

In comparing the transactions per second (TPS) of distributed networks with a
centralized teacher network (as shown in Figure 3(d)), it was observed that
asynchronous, frequently updated distributed models consistently surpass the
teacher model (centralized) in performance. These models not only emulate the
teacher network’s patterns but also demonstrate reduced volatility. In Figure
4(b), an analysis shows that as request numbers significantly increase, our
method outperforms FIRM with 30% improvement in terms of response time in
average. Additionally, DRPC’s lower quadratic coefficient in fitted line than
FIRM confirms its effectiveness, scalability, and stability.

####  VI-C6 Cumulative Distribution Function (CDF) of Response Time Analysis

TABLE 4: Response time at different percentile | 50th | 66th | 75th | 80th | 90th | 95th | 99th | 99.99th
---|---|---|---|---|---|---|---|---
FIRM | 9 | 10 | 12 | 15 | 45 | 120 | 360 | 2900
CoScal | 22 | 36 | 51 | 63 | 95 | 180 | 510 | 1800
DRPC | 11 | 17 | 23 | 28 | 44 | 84 | 210 | 900
Kuscal | 240 | 660 | 960 | 1100 | 1700 | 2300 | 3900 | 7900

We also utilize the CDF in Table 4 to highlight the differences between
different approaches, as well as the results in Figure 4(a). Until the 90th
percentile, the CDF curve of FIRM consistently outperforms DRPC. This
difference appears because FIRM allocates requests to specific nodes rather
than central ones. Our approach, which releases resources on non-critical
chains or nodes, introduces delays on these nodes. Since FIRM does not
allocate resources to these nodes, resulting in over-provisioning, it handles
certain simple requests more quickly than our method. However, beyond the 90th
percentile, our strategy efficiently redistributes the released resources to
nodes critical for QoS enhancement. In summary, DRPC significantly reduces
tail latency by reallocating resources to nodes crucial for enhancing QoS,
outperforming baselines consistently.

####  VI-C7 Reinforcement Learning Convergence Analysis

As illustrated in Figure 4(c), the system requires approximately 300 episodes
to reach convergence. Beyond this point, the reward per action experiences no
significant increase, prompting us to halt the operation. The reward coverage
settles at approximately 0.7, as seen from the figure. This indicates that
CPU, memory usage, and latency are approaching an optimal configuration by
analyzing the reward function. In terms of actual performance, memory
utilization is around 75%, while CPU utilization is approximately 60%. The QoS
is guaranteed, ensuring that 99% of requests can be handled within 210 ms.

##  VII Conclusions and Future Work

In this work, we introduced a novel framework for distributed resource
provisioning, named DRPC. By adopting an asynchronous, parallel, and
differential approach, DRPC facilitates the global optimal allocation of
resources. This accommodates the dynamic nature of microservice-based clusters
while ensuring QoS. Notably, DRPC incorporates DL methodologies for workload
prediction, achieving a higher level of accuracy compared to conventional
gradient-based methods. Additionally, it leverages a distributed RL algorithm
to make informed decisions on scaling strategies, effectively managing the
infinite action-states space associated with microservices. The results, based
on realistic testing and comparisons with state-of-the-art algorithms,
demonstrate that DRPC outperforms the baselines in terms of successful request
rate and average response time, particularly under significantly increased
requests.

The limitation of this approach is that it increases network usage due to the
essential communication of distributed RL. Future research will focus on
automating the asynchronous updating process, which currently requires manual
setting of timing intervals for specific microservices, to improve system
efficiency. In addition, we would like to explore anomaly-aware workloads
management under container-based environment, and incorporate our approach
into large-scale and production environment (e.g. Alibaba Cloud) with further
validations.

## SOFTWARE AVAILABILITY

The codes have been open-sourced to https://github.com/vincent-haoy/DRPC for
research usage.

## References

  * [1] S. Newman, _Building microservices: Designing fine-grained systems_. O’Reilly, 2021.
  * [2] G. Quattrocchi, D. Cocco, S. Staffa, A. Margara, and G. Cugola, “Cromlech: Semi-automated monolith decomposition into microservices,” _IEEE Transactions on Services Computing_ , vol. 17, no. 2, pp. 466–481, 2024.
  * [3] O. M. A. Khan and A. Chandaka, _Developing Microservices Architecture on Microsoft Azure with Open Source Technologies_. Microsoft Press, 2021.
  * [4] M. Xu, L. Yang, Y. Wang, C. Gao, L. Wen, G. Xu, L. Zhang, K. Ye, and C. Xu, “Practice of alibaba cloud on elastic resource provisioning for large-scale microservices cluster,” _Software: Practice and Experience_ , vol. 54, no. 1, pp. 39–57, 2024.
  * [5] S. Pallewatta, V. Kostakos, and R. Buyya, “Placement of microservices-based iot applications in fog computing: A taxonomy and future directions,” _ACM Computing Surveys_ , 2023.
  * [6] S. N. A. Jawaddi, M. H. Johari, and A. Ismail, “A review of microservices autoscaling with formal verification perspective,” _Software: Practice and Experience_ , vol. 52, no. 11, pp. 2476–2495, 2022.
  * [7] S. Xie, J. Wang, B. Li, Z. Zhang, D. Li, and P. C. K. Hung, “Pbscaler: A bottleneck-aware autoscaling framework for microservice-based applications,” _IEEE Transactions on Services Computing_ , vol. 17, no. 2, pp. 604–616, 2024\.
  * [8] V. Kalavri, J. Liagouris, M. Hoffmann, D. Dimitrova, M. Forshaw, and T. Roscoe, “Three steps is all you need: fast, accurate, automatic scaling decisions for distributed streaming dataflows,” in _13th USENIX Symposium on Operating Systems Design and Implementation (OSDI 18)_ , 2018, pp. 783–798.
  * [9] S. Luo, C. Lin, K. Ye, G. Xu, L. Zhang, G. Yang, H. Xu, and C. Xu, “Optimizing resource management for shared microservices: A scalable system design,” _ACM Transactions on Computer Systems_ , 2023.
  * [10] S. Padakandla, “A survey of reinforcement learning algorithms for dynamically varying environments,” _ACM Computing Surveys (CSUR)_ , vol. 54, no. 6, pp. 1–25, 2021.
  * [11] A. Verma, L. Pedrosa, M. Korupolu, D. Oppenheimer, E. Tune, and J. Wilkes, “Large-scale cluster management at google with borg,” in _Proceedings of the Tenth European Conference on Computer Systems_. ACM, 2015.
  * [12] X. He, Z. Tu, X. Xu, and Z. Wang, “Re-deploying microservices in edge and cloud environment for the optimization of user-perceived service quality,” in _International Conference on Service-Oriented Computing_. Springer, 2019, pp. 555–560.
  * [13] B. Burns, J. Beda, and K. Hightower, _Kubernetes: up and running: dive into the future of infrastructure_. O’Reilly Media, 2019.
  * [14] R. S. Kannan, L. Subramanian, A. Raju, J. Ahn, J. Mars, and L. Tang, “Grandslam: Guaranteeing slas for jobs in microservices execution frameworks,” in _Proceedings of the Fourteenth EuroSys Conference 2019_ , ser. EuroSys ’19. New York, NY, USA: Association for Computing Machinery, 2019.
  * [15] X. Hou, C. Li, J. Liu, L. Zhang, Y. Hu, and M. Guo, “Ant-man: towards agile power management in the microservice era,” in _2020 SC20: International Conference for High Performance Computing, Networking, Storage and Analysis (SC)_. IEEE Computer Society, 2020, pp. 1098–1111.
  * [16] G. Yu, P. Chen, and Z. Zheng, “Microscaler: Automatic scaling for microservices with an online learning approach,” in _2019 IEEE International Conference on Web Services (ICWS)_ , July 2019, pp. 68–75.
  * [17] L. Liu, “Qos-aware machine learning-based multiple resources scheduling for microservices in cloud environment,” _arXiv preprint arXiv:1911.13208_ , 2019\.
  * [18] Y. Gan, Y. Zhang, K. Hu, D. Cheng, Y. He, M. Pancholi, and C. Delimitrou, “Seer: Leveraging big data to navigate the complexity of performance debugging in cloud microservices,” in _Proceedings of the twenty-fourth international conference on architectural support for programming languages and operating systems_ , 2019, pp. 19–33.
  * [19] K. Rzadca, P. Findeisen, J. Swiderski, P. Zych, P. Broniek, J. Kusmierek, P. Nowak, B. Strack, P. Witusowski, S. Hand _et al._ , “Autopilot: workload autoscaling at google,” in _Proceedings of the Fifteenth European Conference on Computer Systems_ , 2020, pp. 1–16.
  * [20] M. Xu, C. Song, S. Ilager, S. S. Gill, J. Zhao, K. Ye, and C. Xu, “Coscal: Multifaceted scaling of microservices with reinforcement learning,” _IEEE Transactions on Network and Service Management_ , vol. 19, no. 4, pp. 3995–4009, 2022.
  * [21] H. Qiu, S. S. Banerjee, S. Jha, Z. T. Kalbarczyk, and R. K. Iyer, “$\\{$FIRM$\\}$: An intelligent fine-grained resource management framework for $\\{$SLO-Oriented$\\}$ microservices,” in _14th USENIX Symposium on Operating Systems Design and Implementation (OSDI 20)_ , 2020, pp. 805–825.
  * [22] S. Zhang, T. Wu, M. Pan, C. Zhang, and Y. Yu, “A-sarsa: A predictive container auto-scaling algorithm based on reinforcement learning,” in _2020 IEEE International Conference on Web Services (ICWS)_ , 2020, pp. 489–497.
  * [23] F. Rossi, M. Nardelli, and V. Cardellini, “Horizontal and vertical scaling of container-based applications using reinforcement learning,” in _2019 IEEE 12th International Conference on Cloud Computing (CLOUD)_ , 2019, pp. 329–338.
  * [24] M. Xu and R. Buyya, “Brownout approach for adaptive management of resources and applications in cloud computing systems: A taxonomy and future directions,” _ACM Comput. Surv._ , vol. 52, no. 1, jan 2019.
  * [25] X. Zhou, X. Peng, T. Xie, J. Sun, C. Xu, C. Ji, and W. Zhao, “Benchmarking microservice systems for software engineering research,” in _Proceedings of the 40th International Conference on Software Engineering: Companion Proceeedings, ICSE 2018, Gothenburg, Sweden, May 27 \- June 03, 2018_ , M. Chaudron, I. Crnkovic, M. Chechik, and M. Harman, Eds. ACM, 2018, pp. 323–324.
  * [26] Locust Contributors, “Locust,” 2023, accessed: 2023-06-05. [Online]. Available: https://locust.io/
  * [27] J. Guo, Z. Chang, S. Wang, H. Ding, Y. Feng, L. Mao, and Y. Bao, “Who limits the resource efficiency of my datacenter: An analysis of alibaba datacenter traces,” in _Proceedings of the international symposium on quality of service_ , 2019, pp. 1–10.
  * [28] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [29] L. Buşoniu, R. Babuška, and B. De Schutter, “Multi-agent reinforcement learning: An overview,” _Innovations in multi-agent systems and applications-1_ , pp. 183–221, 2010.
  * [30] M. Zimmer, P. Viappiani, and P. Weng, “Teacher-student framework: a reinforcement learning approach,” in _AAMAS Workshop Autonomous Robots and Multirobot Systems_ , 2014.
  * [31] A. Hussein, M. M. Gaber, E. Elyan, and C. Jayne, “Imitation learning: A survey of learning methods,” _ACM Computing Surveys (CSUR)_ , vol. 50, no. 2, pp. 1–35, 2017.
  * [32] Y. Cao, Y. Zhao, J. Li, R. Lin, J. Zhang, and J. Chen, “Multi-tenant provisioning for quantum key distribution networks with heuristics and reinforcement learning: A comparative study,” _IEEE Transactions on Network and Service Management_ , vol. 17, no. 2, pp. 946–957, 2020.
  * [33] J. Rahman and P. Lama, “Predicting the end-to-end tail latency of containerized microservices in the cloud,” in _2019 IEEE International Conference on Cloud Engineering (IC2E)_. IEEE, 2019, pp. 200–210.
  * [34] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _International conference on machine learning_. PMLR, 2018, pp. 1587–1596.
  * [35] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” in _NIPS 2014 Workshop on Deep Learning, December 2014_ , 2014.
  * [36] M. Xu, C. Song, H. Wu, S. S. Gill, K. Ye, and C. Xu, “Esdnn: Deep neural network based multivariate workload prediction in cloud computing environments,” _ACM Trans. Internet Technol._ , vol. 22, no. 3, aug 2022\.
  * [37] C. Song, M. Xu, K. Ye, H. Wu, S. S. Gill, R. Buyya, and C. Xu, “Chainsformer: A chain latency-aware resource provisioning approach for microservices cluster,” in _Service-Oriented Computing_ , F. Monti, S. Rinderle-Ma, A. Ruiz Cortés, Z. Zheng, and M. Mecella, Eds. Cham: Springer Nature Switzerland, 2023, pp. 197–211.
  * [38] C. Carrión, “Kubernetes scheduling: Taxonomy, ongoing issues and challenges,” _ACM Computing Surveys_ , vol. 55, no. 7, pp. 1–37, 2022.

|  Haoyu Bai received his BSc degree from the University of Melbourne. Now he
is a PhD student at the University of Melbourne, he is also a visiting student
at Shenzhen Institute of Advanced Technology, Chinese Academy of Science. His
main research interest includes efficient microservice-based application and
system management.
---|---
|  Minxian Xu (Senior Member, IEEE) is currently an Associate Professor at
Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences. his
PhD degree from the University of Melbourne in 2019. His research interests
include resource scheduling and optimization in cloud computing. He has co-
authored 60+ peer-reviewed papers published in prominent international
journals and conferences. His PhD thesis was awarded the 2019 IEEE TCSC
Outstanding Ph.D. Dissertation Award. He was also awarded the 2023 IEEE TCSC
Award for Excellence (Early Career Award).
---|---
|  Kejiang Ye (Senior Member, IEEE) received the BSc and PhD degrees from
Zhejiang University, in 2008 and 2013, respectively. He was also a joint PhD
student with the University of Sydney from 2012 to 2013. After graduation, he
works as post-doc researcher with Carnegie Mellon University from 2014 to 2015
and Wayne State University from 2015 to 2016. He is currently a professor with
the Shenzhen Institute of Advanced Technology, Chinese Academy of Science. His
research interests focus on the performance, energy, and reliability of cloud
computing and network systems.
---|---
|  Rajkumar Buyya (Fellow, IEEE) is currently a Redmond Barry distinguished
professor and director with Cloud Computing and Distributed Systems (CLOUDS)
Laboratory, University of Melbourne, Australia. He has authored more than 625
publications and seven textbooks including ”Mastering Cloud Computing”
published by McGraw Hill, China Machine Press, and Morgan Kaufmann for Indian,
Chinese and international markets, respectively. He is one of the highly cited
authors in computer science and software engineering worldwide (h-index=168,
g-index=322, 149,000+citations).
---|---
|  Chengzhong Xu (Fellow, IEEE) received the Ph.D. degree in computer science
and engineering from the University of Hong Kong in 1993. He is the Dean of
Faculty of Science and Technology and the Interim Director of Institute of
Collaborative Innovation, University of Macau. He published two research
monographs and more than 300 peer-reviewed papers in journals and conference
proceedings; his papers received about 17K citations with an H-index of 72.
His main research interests lie in parallel and distributed computing and
cloud computing.
---|---