\section{Introduction}

Matrix Factorization (MF) is a core machine learning technique for applications of
collaborative filtering, such as recommender systems or drug discovery, where a
data matrix $\mathbf{R}$ is factorized into a product of two matrices, such that $\mathbf{R}
\approx \mathbf{U} \mathbf{V}^\top$. The main task in such applications is to predict
unobserved elements of a partially observed data matrix. In recommender
systems, the elements of $\mathbf{R}$ are often ratings given by users to items, while
in drug discovery they typically represent bioactivities between chemical
compounds and protein targets or cell lines.

Bayesian Matrix Factorization (BMF) \cite{bhattacharya2011sparse,Salakhutdinov+Mnih:2008b},
formulates the matrix factorization task as
a probabilistic model, with Bayesian inference conducted on the unknown
matrices  $\mathbf{U}$ and $\mathbf{V}$.  Advantages often associated with BMF include
robustness to over-fitting and improved predictive accuracy, as well as
flexible utilization of prior knowledge and side-data.  Finally, for
application domains such as drug discovery, the ability of the Bayesian
approach to quantify uncertainty in predictions is of crucial importance
\cite{10.1093/bioinformatics/btz335}.

Despite the appeal and many advantages of BMF,
scaling up the posterior inference for industry-scale problems has proven
difficult. Scaling up to this level requires both data and
computations to be distributed over many workers, and so far only very few distributed
implementations of BMF have been presented in the literature.
In \cite{VanderAa+others:2017}, a high-performance computing implementation of BMF using Gibbs
sampling for distributed systems was proposed. The authors considered three different
distributed programming models: Message Passing Interface (MPI),  Global Address Space
Programming Interface (GASPI) and ExaSHARK.
In a different line of work, \cite{Ahn+others:2015} proposed to use a distributed version
of the minibatch-based Stochastic Gradient Langevin Dynamics algorithm for posterior inference
in BMF models.
While a key factor in devising efficient distributed solutions is to be able to minimize
communication between worker nodes, both of the above solutions require some degree of communication
in between iterations.

A recent promising proposal, aiming at minimizing communication and thus
reaching a solution in a faster way, is to use a hierarchical embarrassingly
parallel MCMC strategy \cite{Qin2019}. This technique, called BMF with
Posterior Propagation (BMF-PP),  enhances regular embarrassingly parallel MCMC
(e.g. \cite{Neiswanger:2014:AEE:3020751.3020816,Wang+others:2015}), which does
not work well for matrix factorization \cite{Qin2019} because of
identifiability issues. BMF-PP introduces communication at predetermined
limited phases in the algorithm to make the problem identifiable, effectively
building one model for all the parallelized data subsets, while in previous
works multiple independent solutions were found per subset.

The current paper is based on a realization that the approaches of
\cite{VanderAa+others:2017} and \cite{Qin2019} are compatible, and in fact
synergistic. BMF-PP will be able to parallelize a massive matrix but will be
the more accurate the larger the parallelized blocks are. Now replacing the
earlier serial processing of the blocks by the distributed BMF in
\cite{VanderAa+others:2017} will allow making the blocks larger up to the
scalability limit of the distributed BMF, and hence decrease the wall-clock
time by engaging more processors.

The main contributions of this work are:

\begin{itemize}
\item We combine both approaches, allowing for parallelization
both at the algorithmic and at the implementation level.
\item We analyze what is the best way to subdivide the original matrix into
subsets, taking into account both compute performance \emph{and} model
quality.
\item We examine several web-scale datasets and show that datasets with
different properties (like the number of non-zero items per row) require
different parallelization strategies.
\end{itemize}

The rest of this paper is organized as follows. In Section~\ref{sec:algo}
we present the existing Posterior Propagation algorithm and distributed
BMF implementation and we explain how to combine both. Section~\ref{sec:exp} is
the main section of this paper, where we document the experimental setup and
used dataset, we compare with related MF methods, and present the results both
from a machine learning point of view, as from a high-performance compute point
of view. In Section~\ref{sec:conclusions} we draw conclusions and propose future
work.

\section{Distributed Bayesian Matrix Factorization with Posterior Propagation}
\label{sec:algo}

In this section, we first briefly review the BMF model and then describe the individual
aspects of distributed computation and Posterior Propagation and how to combine them.

\subsection{Bayesian Matrix Factorization}

In matrix factorization, a (typically very sparsely observed) data matrix
$\mathbf{R}\in \R^{N\times D}$ is factorized into a product of two matrices $\mathbf{U} \in
\R^{N\times K}= (\mathbf{u}_1,\ldots,\mathbf{u}_N)^\top$ and $\mathbf{V} = (\mathbf{v}_1,\ldots,\mathbf{v}_D)^\top\in
\R^{D\times K}$. In the context of recommender systems, $\mathbf{R}$ is a rating
matrix, with $N$ the number of users and $D$ the number of rated items.

In Bayesian matrix factorization
\cite{bhattacharya2011sparse,Salakhutdinov+Mnih:2008b}, the data are
modelled as

\begin{equation}\label{eq:BMF_likelihood_normal}
p(\mathbf{R}|\mathbf{U},\mathbf{V}) = \prod_{n=1}^N \prod_{d=1}^D \left[\mathcal{N}\left(r_{nd}|\mathbf{u}_n^\top\mathbf{v}_d,\tau^{-1}\right) \right]^{I_{nd}}
\end{equation}
where $I_{nd}$ denotes an indicator which equals 1 if the element $r_{nd}$ is
observed and 0 otherwise, and $\tau$ denotes the residual noise precision. The
two parameter matrices $\mathbf{U}$ and $\mathbf{V}$ are assigned Gaussian priors.
Our goal is then to compute the joint posterior density $p(\mathbf{U},\mathbf{V}|\mathbf{R}) \propto p(\mathbf{U})p(\mathbf{V})p(\mathbf{R}|\mathbf{U},\mathbf{V})$,
conditional on the observed data.
Posterior inference is typically done using Gibbs sampling, see \cite{Salakhutdinov+Mnih:2008b} for details.

\subsection{Bayesian Matrix Factorization with Posterior Propagation}

In the \emph{Posterior Propagation} (PP) framework \cite{Qin2019}, we start by partitioning
$\mathbf{R}$ with respect to both rows and columns into $I\times J$ subsets
$\mathbf{R}^{(i,j)}$, $i=1,\ldots,I$, $j = 1,\ldots, J$.
The parameter matrices $\mathbf{U}$ and $\mathbf{V}$ are correspondingly partitioned into $I$ and $J$ submatrices, respectively.
The basic idea of PP is to
process each subset using a hierarchical embarrassingly parallel MCMC scheme in
three phases, where the posteriors from each phase are propagated forwards and
used as priors in the following phase, thus introducing dependencies between
the subsets. The approach proceeds as follows (for an illustration, see
Figure~\ref{fig:PP}):

\paragraph{Phase (a):}
    Joint inference for submatrices $(\mathbf{U}^{(1)}, \mathbf{V}^{(1)})$, conditional on data subset $\mathbf{R}^{(1,1)}$:
    \begin{align*}
    &p\left(\mathbf{U}^{(1)},\mathbf{V}^{(1)}|\mathbf{R}^{(1,1)}\right) \propto  p\left(\mathbf{U}^{(1)}\right)p\left(\mathbf{V}^{(1)}\right) p\left(\mathbf{R}^{(1,1)}|\mathbf{U}^{(1)},\mathbf{V}^{(1)}\right).
    \end{align*}
\paragraph{Phase (b):}
    Joint inference in parallel for submatrices $(\mathbf{U}^{(i)},\mathbf{V}^{(1)})$, $i=2,\ldots,I$, and $(\mathbf{U}^{(1)},\mathbf{V}^{(j)})$, $j=2,\ldots,J$, conditional on data subsets which share columns or rows with $\mathbf{R}^{(1,1)}$, and using posterior marginals from phase (a) as priors:
    \begin{align*}
    p\left(\mathbf{U}^{(i)},\mathbf{V}^{(1)}|\mathbf{R}^{(1,1)},\mathbf{R}^{(i,1)}\right) & \propto  p\left(\mathbf{V}^{(1)}|\mathbf{R}^{(1,1)}\right) p\left(\mathbf{U}^{(i)}\right) p\left(\mathbf{R}^{(i,1)}|\mathbf{U}^{(i)},\mathbf{V}^{(1)}\right),\\
    p\left(\mathbf{U}^{(1)},\mathbf{V}^{(j)}|\mathbf{R}^{(1,1)},\mathbf{R}^{(1,j)}\right) & \propto  p\left(\mathbf{U}^{(1)}|\mathbf{R}^{(1,1)}\right) p\left(\mathbf{V}^{(j)}\right) p\left(\mathbf{R}^{(1,j)}|\mathbf{U}^{(1)},\mathbf{V}^{(j)}\right).
    \end{align*}
\paragraph{Phase (c):}
    Joint inference in parallel for submatrices $(\mathbf{U}^{(i)},\mathbf{V}^{(j)})$, $i=2,\ldots,I$, $j=2,\ldots,J$, conditional on the remaining data subsets, and using posterior marginals propagated from phase (b) as priors:
    \begin{align*}
    &p\left(\mathbf{U}^{(i)},\mathbf{V}^{(j)}|\mathbf{R}^{(1,1)},\mathbf{R}^{(i,1)},\mathbf{R}^{(1,j)},\mathbf{R}^{(i,j)}\right)  \\
    &\propto p\left(\mathbf{U}^{(i)}|\mathbf{R}^{(1,1)},\mathbf{R}^{(i,1)}\right) p\left(\mathbf{V}^{(j)}|\mathbf{R}^{(1,1)},\mathbf{R}^{(1,j)}\right)
    p\left(\mathbf{R}^{(i,j)}|\mathbf{U}^{(i)},\mathbf{V}^{(j)}\right).
    \end{align*}

Finally, the aggregated posterior is obtained by combining the posteriors obtained in phases (a)-(c) and  dividing away the multiply-counted propagated posterior marginals; see \cite{Qin2019} for details.

\begin{figure}
	\begin{center}
		\begin{tikzpicture}

    		\draw[help lines] (2,0) grid (6,3);
    		\node[fit={(2,0) (6,3)}, inner sep=0pt, draw=black, ultra thick] (Y) {};
    		\node[fit={(2,0) (3,3)}, inner sep=0pt, draw=black, very thick] (Yi1) {};
    		\node[fit={(2,2) (6,3)}, inner sep=0pt, draw=black, very thick] (Y1j) {};
		\node[align=center] at (2.5+0.02,2.5) {(a)};
	    \node[align=center] at (3.5+0.02,2.5) {(b)};
		\node[align=center] at (4.5+0.02,2.5) {(b)};
	    \node[align=center] at (5.5+0.02,2.5) {(b)};
		\node[align=center] at (2.5+0.02,1.5) {(b)};
	    \node[align=center] at (3.5+0.02,1.5) {(c)};
		\node[align=center] at (4.5+0.02,1.5) {(c)};
	    \node[align=center] at (5.5+0.02,1.5) {(c)};
	    	\node[align=center] at (2.5+0.02,0.5) {(b)};
	    \node[align=center] at (3.5+0.02,0.5) {(c)};
		\node[align=center] at (4.5+0.02,0.5) {(c)};
	    \node[align=center] at (5.5+0.02,0.5) {(c)};

    		\draw[help lines] (0,0) grid (1,3);
    		\node[fit={(0,0) (1,3)}, inner sep=0pt, draw=black, ultra thick] (X) {};
	    \node[align=center,font=\small] at (0.5+0.02,2.5) {$\mathbf{U}^{(1)}$};
	    \node[align=center,font=\small] at (0.5+0.02,1.5) {$\mathbf{U}^{(2)}$};
	    	\node[align=center,font=\small] at (0.5+0.02,0.5) {$\mathbf{U}^{(3)}$};

    		\draw[help lines] (2,4) grid (6,5);
    		\node[fit={(2,4) (6,5)}, inner sep=0pt, draw=black, ultra thick] (W) {};
	    \node[align=center,font=\small] at (2.5+0.02,4.5) {$\mathbf{V}^{(1)}$};
	    \node[align=center,font=\small] at (3.5+0.02,4.5) {$\mathbf{V}^{(2)}$};
	    	\node[align=center,font=\small] at (4.5+0.02,4.5) {$\mathbf{V}^{(3)}$};
	    \node[align=center,font=\small] at (5.5+0.02,4.5) {$\mathbf{V}^{(4)}$};

    	    \node[align=center,font=\large] at (6.5+0.1,1.5) {$\mathbf{R}$};
    	    \node[align=center,font=\large] at (0.5-0.1,3.5) {$\mathbf{U}$};
    	    \node[align=center,font=\large] at (1.5-0.1,4.5) {$\mathbf{V}^\top$};
		\end{tikzpicture}
	\end{center}
	\caption{Illustration of Posterior Propagation (PP) for a data
		matrix $\mathbf{R}$ partitioned into $3 \times 4$ subsets. Subset inferences for the
		matrices $\mathbf{U}$ and $\mathbf{V}$ proceed
		in three successive phases, with posteriors obtained in one phase being
		propagated as priors to the next one. The letters (a), (b), (c) in the matrix $\mathbf{R}$ refer to the
		phase in which the particular data subset is processed. Within each phase, the
		subsets are processed in parallel with no communication. Figure adapted from \cite{Qin2019}.}\label{fig:PP}
\end{figure}

\subsection{Distributed Bayesian Matrix Factorization}

In \cite{VanderAa+others:2017} a distributed parallel implementation of BMF is
proposed. In that paper an implementation the BMF algorithm~\cite{Salakhutdinov+Mnih:2008b} is
proposed that distributes the rows and columns of $\mathbf{R}$ across
different nodes of a supercomputing system.

Since Gibbs sampling is used, rows of $\mathbf{U}$ and rows of $\mathbf{V}$ are independent
and can be sampled in parallel, on different nodes. However, there is a dependency
between samples of $\mathbf{U}$ and $\mathbf{V}$. The communication pattern
between $\mathbf{U}$ and $\mathbf{V}$ is shown in Figure~\ref{fig:comm_struct}.

\begin{figure}
  \centering
  \def.7\columnwidth{.7\columnwidth}
  \import{figs/}{comm_struct.pdf_tex}
  \caption{Communication pattern in the distributed implementation of BMF}
  \label{fig:comm_struct}
\end{figure}

The main contribution of this implementation is how to distribute $\mathbf{U}$ and $\mathbf{V}$
to make sure the computational load is distributed as equally as possible
and the amount of data communication is minimized.  The authors of
\cite{VanderAa+others:2017} optimize the distributions by analysing the sparsity structure
of $\mathbf{R}$.

As presented in the paper, distributed BMF provides a reasonable speed-up for
compute systems up to 128 nodes. After this, speed-up in the strong scaling
case is limited by the increase in communication, while at the same
computations for each node decrease.

\subsection{Distributed Bayesian Matrix Factorization with Posterior Propagation}

By combining the distributed BMF implementation with PP, we
can exploit the parallelism at both levels. The different subsets of $\mathbf{R}$ in
the parallel phases of PP are independent and can thus be computed on in
parallel on different resources.  Inside each subset we exploit parallelism by
using the distributed BMF implementation with MPI communication.

\section{Experiments}
\label{sec:exp}

In this section we evaluate the distributed BMF implementation with Posterior
Propagation (D-BMF+PP) and competitors
for benchmark
datasets. We first have a look at the datasets we will use for the experiments.
Next, we compare our proposed D-BMF+PP implementation with other related matrix
factorization methods in terms of accuracy and runtime, with a special
emphasis on the benefits of Bayesian methods.
Finally we look at the strong scaling of D-BMF+PP and what is a good way to split the
matrix in to blocks for Posterior Propagation.

\subsection{Datasets}

\begin{table}
\centering
\caption{Statistics about benchmark datasets. }
    \begin{tabular}{l|r|r|r|r|r}
    Name  & Movielens \cite{Harper:movielens} & Netflix \cite{Gomez-Uribe:2015:NRS:2869770.2843948}
& Yahoo \cite{Yahoo} & Amazon \cite{Amazon}\\
\hline
    Scale &  1-5  &  1-5  &  0-100 & 1-5 \\
    \# Rows & 138.5K & 480.2K & 1.0M  & 21.2M \\
    \# Columns & 27.3K & 17.8K & 625.0K & 9.7M \\
    \# Ratings & 20.0M & 100.5M & 262.8M & 82.5M \\
\hline
    Sparsity & 189   & 85    &  2.4K  & 2.4M \\
    Ratings/Row & 144   & 209   & 263   & 4 \\
    \#Rows/\#Cols & 5.1   & 27.0 & 1.6   & 2.2 \\
\hline
    K     & 10    & 100   & 100   & 10 \\
    rows/sec ($\times 1000$) & 416   & 15   & 27    & 911 \\
    ratings/sec ($\times 10^6$) & 70    & 5.5  & 5.2   & 3.8 \\
    \end{tabular}%
   \label{tab:datasets}%
\end{table}%

Table~\ref{tab:datasets} shows the different webscale datasets we used for
our experiments. The table shows we have a wide diversity of properties: from
the relatively small Movielens dataset with 20M ratings, to the 262M ratings
Yahoo dataset.  Sparsity is expressed as the ratio of the total number of
elements in the matrix (\#rows $\times$ \#columns) to the number of filled-in
ratings (\#ratings). This metric also varies: especially the Amazon dataset is
very sparsely filled.

The scale of the ratings (either 1 to 5 or 0 to 100) is only important when
looking at prediction error values, which we will do in the next section.

The final three lines of the table are not data properties, but rather properties
of the matrix factorization methods. K is the number of latent dimensions used, which
is different per dataset but chosen to be common to all matrix factorization methods.  The
last two lines are two compute performance metrics, which will be discussed in
Section~\ref{sec:scaling}.

\subsection{Related Methods}

In this section we compare the proposed BMF+PP method to other matrix
factorization (MF) methods, in terms of accuracy (Root-Mean-Square Error or
RMSE on the test set), and in terms of compute performance (wall-clock time).

Stochastic gradient decent (SGD \cite{Teflioudi2012:DSGD++}), Coordinate
Gradient Descent (CGD \cite{Yu2012:ccd++}) and Alternating Least Squares (ALS
\cite{Tan2016:FCP}) are the three most-used algorithms for non-Bayesian MF.

CGD- and ALS-based algorithms update along one dimension at a time while the other
dimension of the matrix remains fixed. Many variants of ALS and CGD exist that
improve the convergence \cite{hsieh2011fastcoord}, or the parallelization degree
\cite{Pilszy2010FastAM}, or are optimized for  non-sparse rating matrices
\cite{parallelSGD}.
In this paper we limit ourselves to comparing to methods, which divide the rating
matrix into blocks, which is necessary to support very large matrices.
Of the previous methods, this includes the SGD-based ones.

FPSGD~\cite{Teflioudi2012:DSGD++} is a very efficient SGD-based library for
matrix factorization on multi-cores. While FPSGD is a single-machine
implementation, which limits its capability to solve large-scale problems, it
has outperformed all other methods on a machine with up to 16 cores. NOMAD
\cite{Yun2014:NOMAD} extends the idea of block partitioning, adding the
capability to release a portion of a block to another thread before its full
completion. It performs similarly to FPSGD on a single machine, and can scale out
to a 64-node HPC cluster.

\begin{table}
\centering
\caption{RMSE of different matrix factorization methods on benchmark datasets.}
\label{tab:related_rmse}
\begin{tabular}{l|r|r|r}
\hline
 Dataset    & BMF+PP & NOMAD & FPSGD  \\
\hline
 Movielens  &  0.76 &  0.77  &  0.77  \\
 Netflix    &  0.90 &  0.91  &  0.92  \\
 Yahoo      & 21.79 & 21.91  & 21.78  \\
 Amazon     &  1.13 &  1.20  &  1.15  \\
\end{tabular}
\end{table}

Table~\ref{tab:related_rmse} compares the RMSE values on the test sets of the
four selected datasets.  For all methods, we used the $K$ provided in
Table~\ref{tab:datasets}. For competitor methods, we used the default
hyperparameters suggested by the authors for the different data sets.  As already
concluded in \cite{Qin2019}, BMF with Posterior Propagation results in
equally good RMSE compared to the original BMF. We also see from the table
that on average the Bayesian method produces only slightly better results,
in terms of RMSE. However we do know that Bayesian methods have significant
statistical advantages \cite{10.1093/bioinformatics/btz335} that stimulate
their use.

\begin{table}
\centering
    \caption{Wall-clock time (hh:mm) of different matrix
    factorization methods on benchmark datasets, running on
    single-node system with 16 cores.}
\label{tab:related_time}
\begin{tabular}{l|r|r|r|r}
\hline
Dataset & BMF+PP & BMF & NOMAD & FPSGD  \\
\hline
Movielens  & 0:07  &  0:14 & 0:08  & 0:09 \\
Netflix    & 2:02  &  4:39 & 0:08  & 1:04 \\
Yahoo      & 2:13  & 12:22 & 0:10  & 2:41 \\
Amazon     & 4:15  & 13:02 & 0:40  & 2:28 \\
\end{tabular}
\end{table}

For many applications, this advantage outweighs the much higher
computational cost. Indeed, as can be seem from Table~\ref{tab:related_time},
BMF is significantly more expensive, than NOMAD and FPSGD, even when taking in to
account the speed-ups thanks to using PP. NOMAD is the fastest method, thanks
to the aforementioned improvements compared to FPSGD.

\subsection{Block Size}

\begin{figure}
   \centering
   \includegraphics[width=0.9\columnwidth]{figs/rmse_vs_blocksize}
    \caption{Block size exploration: RMSE on test set and wall-clock time (hh:mm) for different block size on Netflix.
	Each bubble is the result for a different block size indicated as the number of blocks across the rows and the columns
        inside the bubble. The size of bubbles is an indication for the aspect ratio of the
	blocks inside PP. Smaller bubbles indicate the blocks are more square.}
\label{fig:rmse_vs_blocksize}
\end{figure}

The combined method will achieve some of the parallelization with the PP
algorithm, and some with the distributed BMF within each block. We next compare
the performance as the share between the two is varied by varying the block
size.  We find that blocks should be approximately square, meaning the number
of rows and columns inside the each block should be more or less equal. This
implies the number of blocks across the rows of the $\mathbf{R}$ matrix will be less if
the $\mathbf{R}$ matrix has fewer rows and vice versa.

Figure~\ref{fig:rmse_vs_blocksize} shows optimizing the block size is crucial to
achieve a good speed up with BMF+PP and avoid compromising the quality of the
model (as measured with the RMSE).  The block size in the figure, listed as
$I\times J$, means the $R$ matrix is split in $I$ blocks (with equal amount of rows)
in the vertical direction and $J$ blocks (with equal amount of columns) in the
horizontal direction.  The figure explores different block sizes for the Netflix
dataset, which has significantly more rows than columns ($27\times$) as can be seen
from Table~\ref{tab:datasets}. This is reflected by the fact that the data point
with the smallest bubble area ($20 \times 3$) provides the best trade-off between
wall-clock time and RMSE.

We stipulate that the reason is the underlying trade-off between the amount of
information in the block and the amount of compute per block. Both the amount of
information and the amount of compute can optimized (maximized and minimized
respectively) by making the blocks approximately squared, since both are proportionate
to the ratio of the area versus the circumference of the block.

We will come back to this trade-off when we look at
performance of BMF+PP when scaling to multiple nodes for the different datasets in
Section~\ref{sec:scaling}.

\subsection{Scaling}
\label{sec:scaling}

In this section we look at how the added parallelization of Posterior
Propagation increases the strong scaling behavior of BMF+PP. Strong scaling
means we look at the speed-up obtained by increasing the amount of compute
nodes, while keeping the dataset constant.

The graphs in this section are displayed with a logarithmic scale on both axes.
This has the effect that linear scaling (i.e. doubling of the amount of
resources results in a halving of the runtime) is shown as a straight line . Additionally,
we indicate Pareto optimal solutions with a blue dot. For these solutions one
cannot reduce the execution time without increasing the resources.

We performed experiments with different block sizes, indicated with different
colors in the graphs. For example, the yellow line labeled $16 \times 8$ in Figure
Figure~\ref{fig:netflix_yahoo} means we performed PP with $16$ blocks for the rows of $\mathbf{R}$
and $8$ blocks for the columns of $\mathbf{R}$.

We have explored scaling on a system with up to 128K nodes, and
use this multi-node parallelism in two ways:

\begin{enumerate}
    \item Parallelisation inside a block using the distributed version of
        BMF~\cite{VanderAa+others:2017}.

    \item Parallelism across blocks. In phase (b) we can exploit parallelism
        across the blocks in the first row and first column, using up to $I + J$ nodes
        where $I$ and $J$ is the number of blocks in the vertical, respectively horizontal
        direction of the matrix. In phase (c) we can use up to $I \times J$ nodes.
\end{enumerate}

\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{figs/netflix}
    \includegraphics[width=0.9\columnwidth]{figs/YahooMusic_kddcup2011}
    \caption{Strong scaling results for the Netflix (top) and Yahoo (bottom) datasets.
    X-axis indicates the amount of compute resources (\#nodes). Y-axis is wall-clock time.
    The different series correspond to different block sizes.}
    \label{fig:netflix_yahoo}
\end{figure}

\paragraph{General Trends}
When looking at the four graphs (Figure~\ref{fig:netflix_yahoo}
and~\ref{fig:movielens_amazon}), we observe that for the same amount of nodes
using more blocks for posterior propagation, increases the wall-clock time. The
main reason is that we do significantly more compute for this solution, because
we take the same amount of samples for each sub-block. This means for a
partitioning of $32 \times 32$ we take $1024 \times$ more samples than $1
\times 1$.

On the other hand, we do get significant speedup, with up to $68\times$
improvement for the Netflix dataset. However the amount of resources we need
for this is clearly prohibitively large (16K nodes). To reduce execution time,
we plan to investigate reducing the amount of Gibbs samples and the effect this
has on RMSE.

\paragraph{Netflix and Yahoo}
Runs on the Netflix and Yahoo datasets (Figure~\ref{fig:netflix_yahoo})
have been performed with 100 latent dimensions (K=100).
Since computational intensity (the amount of compute per row/column of $\mathbf{R}$),
is $O(K^3)$, the amount of compute versus compute for these experiments is
high, leading to good scalability for a single block ($1 \times 1$) for Netflix, and almost
linear scalability up to 16 or even 64 nodes. The Yahoo dataset is too large
to run with a $1 \times 1$ block size, but we see a similar trend for $2 \times 2$.

\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{figs/MovieLens20M}
    \includegraphics[width=0.9\columnwidth]{figs/AmazonProductReview}
    \caption{Strong scaling results for the Movielens (top) and Amazon (bottom) datasets.
    X-axis indicates the amount of compute resources (\#nodes). Y-axis is wall-clock time.
    The different series correspond to different block sizes.}
    \label{fig:movielens_amazon}
\end{figure}

\paragraph{Movielens and Amazon}
For Movielens and Amazon (Figure~\ref{fig:movielens_amazon}), we used $K=10$
for our experiments. This implies a much smaller amount of compute for the same
communication cost inside a block.  Hence scaling for $1 \times 1$ is mostly
flat. We do get significant performance improvements with more but smaller
blocks where the Amazon dataset running on 2048 nodes, with $32 \times 32$
blocks is $20 \times$ faster than the best block size ($1 \times 1$) on a
single node.

When the nodes in the experiment align with the parallelism across of blocks
(either $I + J$ or $I \times J$ as explained above), we see a significant drop
in the run time. For example for the Amazon dataset on $32 \times 32$ blocks going
from 1024 to 2048 nodes.

\section{Conclusions and Further Work}
\label{sec:conclusions}

In this paper we presented a scalable, distributed implementation of Bayesian
Probabilistic Matrix Factorization, using asynchronous communication.  We
evaluated both the machine learning and the compute performance on several
web-scale datasets.

While we do get significant speed-ups, resource requirements for these
speed-ups are extermely large.  Hence, in future work, we plan to investigate the effect
of the following measures on execution time, resource usage and RMSE:

\begin{itemize}
    \item Reduce the number of samples for sub-blocks in phase (b) and phase (c).
    \item Use the posterior from phase (b) blocks to make predictions for phase (c) blocks.
    \item Use the GASPI implementation of \cite{VanderAa+others:2017}, instead of the MPI implementation.
\end{itemize}

\section*{Acknowledgments}
The research leading to these results has received funding from the European
Union’s Horizon2020 research and innovation programme under the EPEEC project,
grant agreement No 801051. The work was also supported by the Academy of Finland (Flagship programme: Finnish Center for Artificial Intelligence, FCAI; grants 319264, 292334).
We acknowledge PRACE for awarding us access to Hazel Hen at GCS@HLRS, Germany.


\section{Introduction}
\label{sec:introduct}
While cloud computing has changed human's life by providing service oriented architecture (SOA) as a service (SaaS), Platform as a Service (PaaS), Infrastructure as a Service (IaaS) \cite{fox2009above}, the proliferation of the Internet of Things (IoT), Artificial Intelligence (AI) and edge computing technologies leads us enter the era of post-cloud computing \cite{shi2016promise}. According to the estimation of Cisco Global Cloud Index, the total amount of data created by any device will reach 847 Zettabytes (ZB) annually by 2021 while the Internet traffic will reach 19.5 ZB by then \cite{index2017forecast}. Some of the produced data might require timeliness, involve privacy or cause unpredictable impact on the network. For example, in applications like smart traffic lights, real-time urban surveillance, cloud computing is not able to serve the purpose well \cite{chen2018smart, chen2017smart}.

Rapid technological advances in cloud computing and Internet of Things (IoT) make Smart Cities feasible by integrating heterogeneous computing devices to collaboratively provide variant pervasively deployed smart services \cite{nikouei2019kerman, nikouei2018smart}, in which the capability of data collection and processing at the edge is the key \cite{nikouei2019toward}. Migrating lower-level data processing tasks to edge devices enables the system to meet the requirements for delay-sensitive, mission-critical applications, including smart urban surveillance \cite{chen2017enabling}, instant privacy protection \cite{fitwi2019lightweight, fitwi2019no}, real-time public safety \cite{xu2018real}, etc. However, the heterogeneity and resource constraint at edge necessitate a scalable, flexible and lightweight system architecture that supports fast development and easy deployment among multiple service providers.

Because of many attractive features, such as good scalability, fine granularity, loose coupling, continuous development, and low maintenance cost, the microservices architecture has emerged and gained a lot of interests both in industry and academic community \cite{nagothu2018microservice, nikouei2019decentralized}. Compared to traditional SOAs in which the system is a monolithic unit, the microservices architecture divides an monolithic application into multiple atomic microservices that run independently on distributed computing platforms. Each microservice performs one specific sub-task or service, which requires less computation resource and reduces the communication overhead. Such characteristics make the microservices architecture an ideal candidate to build a flexible platform, which is easy to be developed and maintained for cross-domain applications like smart public safety, smart traffic surveillance systems, etc \cite{chen2018smart, nikouei2018smart}.

Motivated by the observations obtained in our previous research \cite{nikouei2019decentralized, xu2019blendmas}, this work tries to answer three questions:
\begin{itemize}
    \item Is it suitable to run multiple microservices inside one container just considering the performance of the edge device? Which type of microservices with different resource consuming inclination could be put together if the answer is yes?
    \item Is the effect of interference that executes multiple microservices running on fog computing and edge computing scenarios different? and
    \item What is the trade-off between the ``one process per container'' rule and the existing limitation of resource at edge side, such as computation and memory?
\end{itemize}

The rest of this paper is organized as follows. Section \ref{sec:background} discusses the current arts of research on microservices technology and performance evaluation. The microservices deployment policy and performance matrix are explained in Section \ref{sec:systemdesign}. Given comprehensive experimental results, the performance evaluations are discussed in Section \ref{sec:experiment}. Finally, a summary is presented in Section \ref{sec:conclusions}.

\section{Related Work}
\label{sec:background}

\subsection{Microservices Architecture}

Microservices architectures (MSA) are extensions of a Service Oriented Architecture (SOA) \cite{wilde2016approaches}. The traditional SOA uses a monolithic architecture that constitutes different software features in a single interconnected application and database. Relaying on the tightly coupled dependence among functions and components, it is more challenging to adapt to new requirements in an IoT-enabled system than the SOA, which requires scalability, service extensibility, and cross-platform interoperability \cite{datta2018next}. The main difference between MSA and SOA is that a MSA is more decentralized and it distributes all the logic (such as routing, message parsing) in to smart end points \footnote{http://www.opengroup.org/soa/source-book/msawp/p3.htm}. MSA also adopts a light weight applications programming interface (API) gateway for managing services instead of using heavier and more sophisticated Enterprise Service Bus \cite{lander2019microservices}. Each microservice runs its own process and communicates with peers using light weight communication mechanisms (REST API or SOAP protocol \cite{malavalli2015scalable}, etc). The services are built around business capabilities and independently deployed by fully automated deployment tools \cite{pathania2017setting}.

Fine granularity and loose coupling are two significant features in MSA \cite{yu2018survey}. Fine granularity means that there is a bare minimum of centralized management of these services. Moreover, there are instances where the service is not micro \cite{sheng2019micro}. Loose coupling implies that each of microservices components has few dependencies on other separate components, which makes the deployment and development of Micro-services more independent \cite{pautasso2017microservices}.

Granularity should be considered from the standpoint of its eventual impact on service performance as the decision process is influenced by numerous environmental factors when deploying MSA\cite{shadija2017microservices}.

\subsection{Performance Evaluation}

Container technology offers a more lightweight method to abstract the applications from the system environment which allows the microservices to be deployed quickly but also consistently. Compared to VMs, containers not only provide all the libraries and other dependencies but also consume less resource and produce lower overhead \cite{wu2015pseudo, wu2017container}. Applications with different dependencies are wrapped into a container image and shared to various users. These features allows the applications to be run at the edge of the network. Ruchika et al. \cite{tech2016evaluation} evaluated the feasibility of using container technology like docker for IoT applications.

Recent studies compared the performance of containers to virtual machines and showed that in most occasions containers work better than or almost equal to virtual machines. The IBM Research Division measured the performance of Docker in terms of CPU, memory, disk I/O and compared the result with KVM \cite{felter2015updated}. The research showed that in all evaluated cases container works better than VMs. In a study researchers have compared the performance of Docker with VMs when running the Spark applications  \cite{bhimani2017accelerating}. They claimed that Spark works better with docker for calculation intensive applications. Similar work has been done in big data area with the interference caused by neighbor containers running big data microservices \cite{ye2017performance}. The performance of different single-board computer devices are evaluated by deploying containerized benchmarks \cite{morabito2017virtualization}, which focused on hardware comparison and the choice of one device rather than another, however, it did not pay much attention on container analysis. An optimized version of Docker called Docker-pi was proposed for container deployment on single-board machines to speed up docker deployment \cite{ahmed2018docker}.

Researchers have studied the performance of collocated microservices running on a common host and compared the result in container with VMs \cite{sharma2016containers}. Evaluation is also reported on how microservices are interfered with each other in cloud computing environment \cite{jha2018holistic}, in which the performances between processed in a container and the situation where they follow the ``one process per-container'' rule are compared. The conclusion was that putting microservices inside one container has a positive effect to the overall performance. However, the work is based on a cloud computing platform.

A survey of VM management lists the virtualization frameworks for edge computing and elaborates the techniques tailored for it \cite{tao2019survey}. Focusing on Docker, a container based technology as an edge computing platform, researchers evaluated four criteria and conclude that Docker is a viable candidate \cite{ismail2015evaluation}.

To the best of our knowledge, there is not a reported effort that evaluates the interference effect of containers in the edge computing environment. The trade-off between lower overheads and the ``one process per container'' rule seems more important considering the resource limitation at the edge.

\section{Microservices Benchmark Evaluation}
\label{sec:systemdesign}
For microservices development, we choose Docker container owing to popularity and familiarity in container community. Docker wraps up the application with all its dependencies into a container so that it can easily be executed on any Linux server (on-premise, bare metal, public or private cloud). It uses layered file system images along with the other Linux kernel features (namespace and cgroups) for the management of containers. This feature allows Docker to create any number of containers from a base image, which are copies of the base image wrapped up with additional features. This also reduces the overall memory and storage requirements that eases fast startup.

To measure the performance of the containerized microservices deployment policy, several widely used benchmarks are considered in designing the microservices benchmark based on typical computing paradigms used in the smart public safety system at the edge. The key metrics of the benchmark are defined as following:

\begin{enumerate}
    \item CPU performance: To evaluate the performance of the CPU, we chose the Linpack benchmarks which measure the system's floating point computing power by solving linear equations. And the capability of the CPU is measured in terms of FLOPS (floating point operations per second).
    \item Memory performance:We chose the STREAM benchmarks to measure the performance of the memory. STREAM is a simple synthetic benchmark program that measures sustainable memory bandwidth (in MB/s) and the corresponding computation rate for simple vector kernels.
    \item Disk performance: To measure the disk I/O capability, we chose the Bonnie++ benchmark which gives the results in terms of input, output, rewrite (in Kb/s) and seeks (in per second).
\end{enumerate}

\begin{figure}[b]
    \centering
    \includegraphics[height=0.42\textwidth]{JPG/Flow.jpg}
    \caption{Containerized benchmark evaluation workflow.}
    \label{fig:workflow}
\end{figure}

For each metric evaluation, the containerized microservices are deployed both on fog (desktop) side and edge (Raspberry Pi 4) side. Given above defined benchmarks, the performance and the interference effect of the containerized microservices are evaluated based on a set of microservice deployment policies. Figure \ref{fig:workflow} illustrates the whole containerized benchmark evaluation workflow, and the microservice deployment policies are described as following four cases:

\begin{itemize}
    \item \emph{Case 1}: One microservice is developed as single container and only one container is running on host machine. Since we would put multiple microservices in one container later, the system resource is constrained by using the control groups or Cgroups, a Linux kernel feature. Such case is considered as a baseline for the entire evaluation.
    \item \emph{Case 2}: Multiple microservices are developed as single container and only one container is running on host machine. In this case we do not employ Cgroups so the microservices can have access to all the system resources. This case is used to evaluate resource competition caused by wrapping multiple microservices into one container.
    \item \emph{Case 3}: Each container only holds single microservice and multiple containers are running on host machine without employing Cgroups. This case could evaluate resource competition caused by multiple containers deployed on single host platform.
    \item \emph{Case 4}: Each container only holds single microservice and multiple containers are running on host machine. Like case 1, this case employs Cgroups to limit the system resource during test.
\end{itemize}

\section{Experimental Analysis}
\label{sec:experiment}
\subsection{Experimental Setup}
A concept-proof prototype system has been implemented on a simulating SPS scenarios. Table \ref{tab:testbed} describes configurations of experimental testbed. The desk top simulates fog node while Raspberry Pi 4 acts as edge node.

\begin{table}[ht]
\caption{Configuration of Experimental Testbed.}
\label{tab:testbed}
\begin{center}
\begin{tabular}{|l|p{2.5cm}|p{3.5cm}|}
\hline
\rule[-1ex]{0pt}{3.5ex} \textbf{Device} & Desktop & Raspberry Pi 4 \\
\hline
\rule[-1ex]{0pt}{3.5ex} \textbf{CPU} & 4.5GHz Intel Core TM (4 cores) & Quad core Cortex-A72 (ARM v8), 1.5GHz \\
\hline
\rule[-1ex]{0pt}{3.5ex} \textbf{Memory} & 4GB DDR3 & 4GB LPDDR4 RAM \\
\hline
\rule[-1ex]{0pt}{3.5ex} \textbf{Storage} & 50G HHD & 32GB (microSD card) \\
\hline
\rule[-1ex]{0pt}{3.5ex} \textbf{OS} & Ubuntu 16.04 & Raspbian GNU/Linux (Jessie) \\
\hline
\end{tabular}
\end{center}
\end{table}
\vspace{-10pt}

For workload size and configuration on fog, we provide a specific configuration for each microservice. For Linpack, we consider the matrix of size N as 15000. We also considered the problem size (number of equations to solve) as 15000. We configure the STREAM by setting the array size as 60M and DNTIMES as 200. The total memory requirement for this configuration is 1373.3 MiB. For Bonnie++, we considered the file size as 8192 MB and set the uid to use as root.

For workload size and configuration on edge (Raspberry Pi 4) side, each microservice, we provide a specific configuration. For Linpack, we consider the matrix of size N as 15000. We also considered the problem size (number of equations to solve) as 15000. We configure the STREAM by setting the array size as 30M and DNTIMES as 100. The total memory requirement for this configuration is 1373.3 MiB. For Bonnie++, we considered the file size as 4096 MB and set the uid to use as root.

The experimental test is conducted based on performance of benchmarks and microservices deployment policies. For the ease of presentation of the results, we used the following abbreviations for the microservices used in the simulation of SPS system, they are namely Linpack (L), STREAM (S), and Bonnie++ (B). To get experimental results, for case 1, each microservice (L, S and B) is performed 30 iterations. For case 2, case 3 and case 4, since the microservices are in all possible combinations and corresponding isolated running time of the containers are not identical, it is not suitable to repeat the process for a certain number of iterations. Hence, for these cases, we keep the container running for a certain period of time (60 minutes) and compute the average performance.

\subsection{Evaluation on Fog Side}

In this paper we define a desktop as the fog side, the most usual physical machine in a LAN system and  the role of a manager of the other edge computing device in our IoT system.

\emph{1. CPU Performance and Analysis}:
According to the CPU architecture, we chose the Linpack version for Intel and implemented it in Docker. We ran the container(s) and recorded the performance of Linpack microservice as described in Section \ref{sec:systemdesign}. Figure \ref{fig:cpu_fog} shows the average CPU performance of running Linpack in terms of GFLOPS. Compared with case 1, all combination of microservices have more computation overhead on host performance. The L+B in case 2 has the least impact on CPU performance while L+L in case 3 is the worst result that has 64\% degradation compared to the baseline L in case 1. Owing to higher computing operations, like float calculation, used by Linpack (L), so L+L combination shows worst CPU performance in all deployment policies. Furthermore, in case 3, control groups (Cgroups) are disabled so that containers have to compete for the system resource. Therefore, test results in case 4 have better performance than those in case 3.

\begin{figure}[b]
    \centering
    \includegraphics[width=.45\textwidth]{JPG/Linpack-fog.jpg}
    \caption{CPU Performance on Fog.}
    \label{fig:cpu_fog}
\end{figure}

\noindent \emph{2. Memory Performance and Analysis}:
To evaluate the memory performance, we wrapped the STREAM benchmark inside a Docker container. The STREAM benchmark employs four vector operations and gives a result of bandwidth. The average result of all the four vector operations (COPY, SCALE, ADD and TRIAD) demonstrate similar shapes in bar chart, so we just present the result of COPY operation as memory performance. As shown in Fig. \ref{fig:mem_fog}, except for combination S+S in case 2, the others have limited performance deduction compared to the baseline in case 1. Given the results, wrapping memory type workloads microservices as one container introduce higher memory overhead on fog computing environment. Similar to the results in CPU performance, all the combinations show best performance in case 4 when Cgroups is enabled.

\begin{figure}[t]
    \centering
    \includegraphics[width=.45\textwidth]{JPG/Copy-fog.jpg}
    \caption{Bandwidth of COPY on Fog.}
    \label{fig:mem_fog}
    \vspace{-10pt}
\end{figure}

\begin{figure}[b]
    \centering
    \includegraphics[width=.48\textwidth]{JPG/IO_fog.png}
    \caption{a) Performance of Input on Fog; b) Performance of Output on Fog.}
    \label{fig:IO_fog}
\end{figure}

\noindent \emph{3. Disk I/O Performance and Analysis}:
The Bonnie++ benchmark evaluates various performance parameters in which only sequential block input and output scenarios are considered during test. The performance evaluation is based on I/O throughput in terms of MB/s. Figure \ref{fig:IO_fog} shows that the input performance is better than output performance on fog node. The results show the similar interference pattern that running multiple instances of Bonnie++ leads to lower I/O throughput, no matter cases that are deployed in a single container or separate containers. From case 2 to case 4, the performance is nearly equal to the baseline case 1 except for the combination of B+B.

In our IoT system, for instance, a local smart surveillance, a fog device plays the role of managing and monitoring the performance of the edge device like smart cameras. Consider that the nodes could be numerous, the fog device is more sensitive to the granularity of the containers. Based on the evaluation results, generally it seems more practical to reduce the number of containers which are CPU consuming inclined. Meanwhile we should avoid deploying same type microservices together. The evaluation gives a general direction of setting up the fog device, we should adjust the deployment of the specified microservices according to the actual situation of the system and capability of the hardware.

\subsection{Evaluation on Edge Side}
The experimental process on edge ( Raspberry Pi) is quite the same as fog side. However, given limited resource on edge node, we set the time interval to 120 minutes for the following evaluations.

\noindent \emph{1. CPU Performance and Analysis}:
Since the CPU architecture of Pi 4 is ARM, we picked another Linpack version and adjusted the codes to the system environment. Then We containerized the benchmark and ran experiments as per the various case described in Section \ref{sec:systemdesign}. Figure \ref{fig:cup-edge} shows the CPU performance in terms of MFLOPS.

Similar to the results on fog side, all the combinations have a significant impact on the computing performance of host except for the combination of L+B. The combination of L+S in case 2 has worst performance, so it causes 41\% degradation compared to the baseline of case 1. But all three combinations demonstrate the best performance in case 4, knowing the fact that the microservices are separated in two containers with the control group enabled. This is quite different from the case in physical machine where L+B in case 4 did not give a good performance.

The results also indicate that computing type workloads tend to have the worst performance when working with memory consuming ones when they are wrapped into a single container with cgroups employed. And it shows that CPU consuming microservice works well with disk I/O workloads and the performance is much closer to the baseline relatively compared to the fog environment. In addition, as Fig. \ref{fig:cup-edge} shows, case 4 has better performance than case 3 while running all combinations. The reason behind is that in case 3 Cgroups are disabled so the containers have to compete for the system resource, however, the containers have no interference with each other in case 4 when Cgroups are enabled.

\begin{figure}[t]
    \centering
    \includegraphics[width=.45\textwidth]{JPG/Linpack-edge.jpg}
    \caption{CPU Performance on Edge.}
    \label{fig:cup-edge}
    \vspace{-10pt}
\end{figure}

\noindent \emph{2. Memory Performance and Analysis}:
Similar to fog side, only the of COPY operation is used for evaluation for memory performance in terms of GB/s. As shown in Fig. \ref{fig:copy-edge}, all four vector operations generate almost the same shapes in charts. Compared to the baseline in case 1, the results shows that all combinations have limited performance reduction except for combination S+L in case 2. Hence, enclosing memory type workloads inside one container as microservice is not suggested on edge computing platform. Like the results in CPU performance, the combinations work best in case 4 that microservices are separated in containers with resource control.

\begin{figure}[t]
    \centering
    \includegraphics[width=.45\textwidth]{JPG/Copy-edge.jpg}
    \caption{Bandwidth of COPY on Edge.}
    \label{fig:copy-edge}
    \vspace{-10pt}
\end{figure}

\noindent \emph{3. Disk I/O Performance and Analysis}:
To evaluate Disk I/O performance on edge device, we used Bonnie++ microservice that creates a large file at least twice the size of inbuilt memory, and only sequential block input and output are considered during the test. The throughput performance of running microservices on edge is shown in Fig. \ref{fig:IO-edge}. The result shows the similar interference pattern indicating that running multiple instances of Bonnie++ creates higher resource contention when they are deployed either in a single container or separate containers. The results also indicate that the performance of combinations is nearly equal to the baseline performance except for the case of running multiple instances of Bonnie++ (B+B). All the combinations of deployment policies show very small deviation from the baseline of case 1 except for case of (B+B).

Basically in our IoT system, the edge device like Raspberry Pi suffer from the limited resource so we should carefully consider the combination of microservices or simply isolate some inside a single container. Based on the results, generally we should deploy the CPU consuming inclined microservice alone, in smart surveillance system for example, the Computer Vision(CV) type containerized services. And it is considerable to wrap memory and disk inclined ones together to reduce the overheads and spare more system resources.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=.48\textwidth]{JPG/IO_edge.png}
    \caption{a) Performance of Input on Edge; b) Performance of Output on Edge.}
    \label{fig:IO-edge}
    \vspace{-10pt}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}
In this paper, a set of microservice deployment policies are proposed and performance matrix are considered. To evaluate the microservice deployment policy, a simulating test is developed and tested both on fog and edge computing platform. Given the experimental results, we conclude as following:

\begin{itemize}
    \item Under edge computing environment, running multiple microservices in a single container is an option considering the comparable performance (except for certain cases) if we ignore the ``\emph{one process per container}'' rule. The story behind the rule is that when deploy multiple process such as microservices in a single container we may face operations problem. For instance, we may have to replace the whole container if one of the microservices need updates which could be a disaster in a large scale. Considering the resource limitation, however, sometimes we have no choices.
    \item Memory intensive microservices are not heavily effected except for competing with same type ones inside one container. Disk I/O microservices is quite the same situation, however, both of the two kinds have a significant effect to the CPU intensive process.
\end{itemize}

This is a preliminary study toward a microservices architecture based smart public safety as an edge service. Our ongoing efforts include building a larger scale edge surveillance system with 32 smart cameras, which will enable a more comprehensive study for deeper insights.

\bibliographystyle{IEEEtranS}