\section{Introduction}
The growing diffusion of Web-based services allows an
increasingly large number of users to access and consume online
content. People access web pages,
purchase items in e-commerce sites, watch online content through
streaming services, or interact within social networks and media.
Understanding the factors that characterize user preferences
and shape their future behavior is crucial in order to provide
users with a better experience through the recommendation of new
content and data that users are likely to appreciate.

Collaborative filtering approaches to recommendation were extensively
investigated by the current literature
\cite{Aggarwal:2016:RST:2931100}. Among these, latent variable models
\cite{Hofmann:2004,Ruslan:2008,Wang:2011,Kabbur:2013}
gained substantial attention, due to their capabilities in modeling
the hidden causal relationships that ultimately influence user
preferences.  Recently, however, new approaches based on neural
architectures (see \cite{ZhangYS17aa} for a comprehensive survey) were
proposed, achieving competitive performance with respect to the
current state of the art. Also, new paradigms based on the combination
of deep learning and
probabilistic latent variable modeling
\cite{kingma14,RezendeMW14} were proven quite
successful in domains such as computer vision and speech
processing. However, their adoption for modeling user preferences is
still unexplored, although recently it is starting to gain attention
\cite{Liang:2018:VAC:3178876.3186150,Li:2017}

The aforementioned approaches rely on the ``bag-of-word''
assumption: when considering a user and her preferences, the order of
such preferences can be neglected and all preferences are
exchangeable. This assumption works with general user trends which
reflect a long-term behavior. However, it fails to capture the
short-term preferences that are specific of several application
scenarios, especially in the context of the Web. Sequential
data can express causalities and dependencies
that require ad-hoc modeling and algorithms. And in fact, efforts to
capture this notion of causality have been made, both in the context of latent variable modeling
\cite{Rendle:2010,Barbieri:2013,Tavakol:2014,He:2017} and deep learning
\cite{Hidasi:2016,Devooght:2017,Wu:2017,Tang:2018}.

In this paper, we consider the task of sequence recommendation from the
perspective of combining deep learning and latent variable
modeling. Inspired by the approach in
\cite{Liang:2018:VAC:3178876.3186150}, we assume that at a given
timestamp the choice of a given item is influenced by a latent
factor that models user trends and preferences. However, the latent
factor itself can be influenced by user history and modeled to capture
both long-term preferences and short-term behavior. As a matter of
fact, the recent studies in the recurrent neural network literature
\cite{ChungGCB14,Greff:2017,Vaswani:2017} demonstrate the capabilities of these
architectures in capturing both long-term and short-term
relationships. It is hence natural to exploit them in a variational
setting.

Our contribution can be summarized as follows.
\begin{itemize}
\item We review the framework of Variational Autoencoders and discuss
  its adoption for the problem of modeling implicit preferences in a
  collaborative filtering scenario.
\item We extend the framework to the case of sequential
  recommendation, where user's preferences exhibit temporal
  dependencies. In particular, we discuss different modeling
  alternatives for tackling this task through the adoption of
  recurrent neural network architectures that model latent
  dependencies at different abstraction levels.
\item We evaluate the proposed framework on standard benchmark
  datasets, by showing that \emph{(a)} approaches not considering
  temporal dynamics are not totally adequate to model user
  preferences, and \emph{(b)} the combination of latent variable
  and temporal dependency modeling produces a substantial
  gain, even with regard to other approaches that only focus on
  temporal dependencies through recurrent relationships.
\end{itemize}

The rest of the paper is organized as follows. Section
\ref{sec:related} discusses the recent contributions in the current
literature and provides a systematic review of the approaches related
to our task of interest. Sections~\ref{sec:back} and~\ref{sec:vae} propose
the modeling of user preferences in a variational setting, by
illustrating how the general framework can be adapted to the case of
temporally ordered dependencies. The effectiveness of the proposed
modeling is illustrated in section~\ref{sec:exp}, and pointers to
future developments are discussed in section~\ref{sec:conc}.

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

Recommender systems have been extensively studied over the last two
decades. Within the collaborative filtering framework, recommendation is
essentially modeled as a prediction problem: Given a user and an item,
we would like to predict the preference of the user for that item,
exploiting the user's past choices, i.e. her history.

Most approaches disregard the temporal
order of the preferences in a user's history. Among these,
latent variable models
\cite{Hofmann:2004,Ruslan:2008,Salakhutdinov:2008,Rendle:2009,flda,Ning:2011,Wang:2011,Rendle:2012,Kabbur:2013,2014Barbieri,BM2011}
were proven extremely effective
in modeling user preferences and providing reliable recommendations.
Essentially, these approaches embed users and items into latent spaces
that translate relatedness into geometrical closeness. The latent
embeddings can be used to decompose
the large sparse \emph{preference matrix}
\cite{Salakhutdinov:2008,Ruslan:2008,flda}, to devise item
similarity~\cite{Ning:2011, Kabbur:2013}, or more generally to
parameterize probability distributions for item preference \cite{Rendle:2009,
  Hofmann:2004,Wang:2011,Rendle:2012} and sharpen the prediction
quality by means of meaningful priors.

The recent literature is currently focusing on deep learning, which
shows substantial advantages over traditional approaches.
For example, \emph{Neural Collaborative Filtering} (NCF) \cite{He17}
generalizes matrix factorization to a non-linear setting, where users,
items and preferences are modeled through a simple multilayer
perceptron network that exploits latent factor transformations.

Notably, prominent deep learning approaches to collaborative
filtering are based on the idea of autoencoding the features from the
preference matrix.  \emph{AutoRec} \cite{Sedhain15} exploits
autoencoders to encode preference histories. Unseen preferences can
be devised by looking at the reconstructed decoding, which is shaped to
include scores for all possible items of interest.
Autoencoders are also amenable to consider side information
\cite{strub15} to mitigate the sparsity of the data and to tackle the
cold start problem.

Hybrid approaches that integrate latent variable modeling and deep
learning have also gained attention. \emph{Collaborative Deep
  Learning} \cite{Wang15} embeds a \emph{stacked denoising
  autoencoder} \cite{Vincent10} into a Bayesian matrix factorization
setting.  Similarly, \emph{AutoSVD++} \cite{Zhang17} exploits
the notion of \emph{contractive autoencoder} \cite{Rifai:2011} to
learn latent item representations that are integrated into the
\emph{SVD++} model \cite{Koren08}.

The introduction of the variational autoencoding framework
\cite{kingma14,RezendeMW14} has suggested a tighter coupling
between deep learning and latent variable modeling.
\emph{Collaborative Variational Autoencoder} (CVA) \cite{Li:2017} and
\emph{Hybrid Variational Autoencoder} (HVAE)
\cite{2018arXiv180801006G} exploit side information to feed a
variational autoencoder whose goal is to produce a latent
representation of the items. In CVA, the preference matrix is hence
modeled by combining user and item embeddings with the item
latent representations, while HVAE uses another variational
autoencoder to reproduce the whole users' preference history.
By contrast,
\cite{Liang:2018:VAC:3178876.3186150} proposes a neural generative
model where a user's history is modeled through a multinomial
likelihood conditioned to a latent user representation which in turn
is modeled through a variational autoencoder.

Within the context of collaborative filtering, a strong effort has
also been made to model temporal dynamics within the history of user preferences
\cite{Quadrana:2018}.
The \emph{Factorizing Personalized Markov Chain} model (FPMC)
\cite{Rendle:2010}, for example, proposes a combination of matrix
factorization and Markov chains. FPMC considers personalized
first-order transition probabilities between items, which are modeled
by decomposing the underlying tensor through user and item
embeddings.
Transition probabilities can also be measured by exploiting more
sophisticated modeling \cite{HeM16,He:2017}, where users are mapped into
translation (latent) vectors operating on item sequences and
consequently a transition corresponds to a geometric affinity of these
latent vectors.
Orthogonally, Markov dependencies can also be exploited to model
dependencies between latent variables~\cite{Barbieri:2013}, thus
resulting in richer formalizations and more accurate recommendations.

Recently, a revamped interest for sequence-based recommendation has
taken place, motivated by both the success of recurrent neural
networks \cite{Cho14,Greff:2017} in domains such as language modeling,
and the need to focus on \emph{session-based}
recommendations~\cite{Hidasi:2016, Tan16}, i.e. recommendation that do not rely
on a user model and instead can cope with single anonymous preference
sessions.

\emph{GRU4Rec} \cite{Hidasi:2016} proposes a recurrent neural network
model based on \emph{Gated Recurrent Units} to predict the next item
in a user session, based on the history seen so far. Since its
introduction, this model has witnessed several evolutions, and similar
architectures were proposed in
\cite{Devooght:2017,Liu:2016,Twardowski:2016, Wu16,Tan16,Jannach:2017,Quadrana:2017}.
Recurrent networks were also exploited to strengthen matrix
factorization, by producing history-aware embeddings of users and
items. The \emph{RRN} approach proposed in \cite{Wu:2017} combines two
recurrent networks, whose output at any time-step, relative to a specific
user and item, can be hence exploited to predict the current
preference.

Finally, the CASER (\emph{Convolutional Sequence
  Embedding Recommendation}) model \cite{Tang:2018} proposes an
approach that departs from RNN modeling and instead exploits a
\emph{convolutional neural network}, by transforming a sequence into
a matrix built from the concatenation of the embeddings of the items
appearing in the sequence. The matrix can hence feed convolutional
layers that can extract relevant useful features for predicting
the next items.

\section{Background}\label{sec:back}

The framework of variational autoencoders draws from the idea of
latent variable modeling \cite{Murphy:2012:MLP:2380985}. Essentially,
we can assume a $K$-dimensional
latent variable space $\mathcal{Z}$, upon which we can devise a probability
density function $P_{\theta}(\bs{z})$ for
$\bs{z}\in\mathcal{Z}$, where $\theta$ represents a set of density
parameters. The datapoints
$\bs{X} = \left\{\bs{x}_{1}, \ldots, \bs{x}_{M}\right\}$ we observe can be
modeled through a dependency $P_\phi(\bs{x}|\bs{z})$, so
that the overall likelihood of $\bs{X}$ can be specified by marginalizing
over $\mathcal{Z}$:
\begin{align*}
P(\bs{X}) = & \prod_i P(\bs{x}_i) = \prod_i\int P_\phi\left(\bs{x}_i|\bs{z}\right) P_\theta(\bs{z}) \diff \bs{z} \, .
\end{align*}

Within the classical VAE framework \cite{kingma14,RezendeMW14}, \bs{z} is assumed to
be distributed according to a standard normal distribution, that is
$\theta = \{\mathbf{0},\mathbf{I}_K\}$ and consequently
$\mathbf{z} \sim \mathcal{N}(\mathbf{0},\mathbf{I}_K)$. The
key intuition here is that even complex dependencies can be generated
starting from normally distributed variables. Thus, $P_\phi(\bs{x}_i|\bs{z})$ is
parameterized by the function $f_\phi(\bs{z})$, representing any
(unknown) transformation of $z$ that expresses such a dependency. We
can model this transformation through a neural network, so that
$\phi$ represents the set of network parameters to be optimized.
Hence, we can model the inference problem for $P(\bs{X})$ as
optimization problem, where we aim at finding the optimal parameter
set $\phi$ that maximize $P(\bs{X})$.

However, $P(\bs{X})$ is typically intractable and
we need to find an approximation. Variational inference
\cite{Blei2017} usually
tackles this approximation by introducing a proposal distribution
$Q(\bs{z}|\bs{x})$, which approximates the true posterior
$P(\bs{z}|\bs{x})$. The relationship between the likelihood and the
proposal distribution is given by the following equation:
\begin{align*}
  \log P(\mathbf{x}) & = \int Q(\bs{z}|\bs{x}) \log P(\bs{x}) \diff \bs{z}\\
  & =  \int Q(\bs{z}|\bs{x}) \log \frac{P(\bs{x}|\bs{z})P(\bs{z})}{P(\bs{z}|\bs{x})} \diff \bs{z}\\
  & =  \int Q(\bs{z}|\bs{x}) \log P(\bs{X}|\bs{z}) \diff \bs{z}
   +\int Q(\bs{z}|\bs{x}) \log \frac{Q(\bs{z}|\bs{x})}{P(\bs{z}|\bs{x})} \diff \bs{z} \\
& - \int Q(\bs{z}|\bs{x}) \log \frac{Q(\bs{z}|\bs{x})}{P(\bs{z})} \diff \bs{z} \\
   & =  E_{z\sim Q}[\log P(\bs{x}|\bs{z})]
     +\mathit{KL}\left(Q(\bs{z}|\bs{x}) \|P(\bs{z}|\bs{x})\right) \\
& - \mathit{KL}\left(Q(\bs{z}|\bs{x})  \| P(\bs{z})\right) \, .
\end{align*}
By rearranging, we obtain
\begin{equation}
\begin{split} \label{eq:elbo}
\log P(\mathbf{x}) &- \mathit{KL}\left[Q(\bs{z}|\bs{x})
\|P(\bs{z}|\bs{x})\right] \\
& = E_{z\sim Q}[\log P(\bs{x}|\bs{z})] - \mathit{KL}\left[Q(\bs{z}|\bs{x})  \| P(\bs{z})\right] \, .
\end{split}
\end{equation}
This equation is crucial to the development variational autoencoders. The
left-hand side represents the term to maximize, plus an error term due
to approximating the true posterior $P(\bs{z}|\bs{x})$ with
$Q(\bs{z}|\bs{x})$. A  good choice of $Q$ gives a small
(hopefully zero) approximation error, and
consequently allows to directly optimize the likelihood. The
right-hand side of the equation (called Evidence Lower Bound, ELBO)
represents the equivalent of the
left-hand side, but with the added bonus that, for a given choice of
$Q$ it becomes tractable and amenable to optimization. In particular,
by noticing that
\begin{equation*}
\begin{split}
\mathit{KL}&\left[\mathcal{N}(\bs{\mu}, \bs{\Sigma}) \|
  \mathcal{N}(\bs{m}, \bs{S})\right]=\\
  & \frac{1}{2}\left(log\frac{|\bs{S}|}{|\bs{\Sigma}|} -K
+ \mathit{tr}\left(\bs{S}^{-1}\bs{\Sigma}\right)
+ \left(\bs{m} - \bs{\mu}\right)^T\bs{S}^{-1}\left(\bs{m} -
  \bs{\mu}\right)
\right) \, ,
\end{split}
\end{equation*}
we can choose to model  $Q(\bs{z}|\bs{x}) = \mathcal{N}(\bs{z}|\mu_\lambda(\bs{x}),
\Sigma_\lambda(\bs{x}))$, so that the term $\mathit{KL}\left[Q(\bs{z}|\bs{x})  \|
  P(\bs{z})\right]$ has a closed form.
Again, we can model the parameters $\mu_\lambda(\bs{x})$
and $\Sigma_\lambda(\bs{x})$ through a neural
network trained on $\bs{x}$ and parameterized by $\lambda$. A
particular case that we shall use throughout this paper is given by
$\bs{\Sigma}_\lambda (\bs{x}) = \mathit{diag}\left(\sigma_{\lambda,1}
  (\bs{x}), \ldots, \sigma_{\lambda,K} (\bs{x})\right)$.

As a consequence, the modeling can resort to two neural networks: the
first one for ``encoding'' an input $\bs{x}$ into a latent variable
$\bs{z}$ by means of $Q_\lambda(\bs{z} |\bs{x})$. The second one
``decodes'' the latent representation $\bs{z}$ into the corresponding
$\bs{x}$ by means of $P_\phi(\bs{x}|\bs{z})$. The learning process is
governed by the loss function given by the right-hand side of
eq.~\ref{eq:elbo}, computed on the training data \bs{X}.

A problem with this loss is that term $E_{\bs{z}\sim Q}[\log
P(\bs{x}|\bs{z})]$ is typically approximated by
sampling the values $\bs{z}$ according to $Q(\bs{z} |\bs{x})$. However, the
sampling is a nondeterministic function that depends on the parameters
$\bs{\mu}_\lambda$ and $\bs{\Sigma}_\lambda$ (which in turn depend on $\lambda$) and
is not differentiable. To overcome
this, \cite{kingma14} propose a \emph{reparametrization trick}: instead
of sampling \bs{z}, we can sample an auxiliary noise variable
\bs{\epsilon} according to a fixed distribution $P(\bs{\epsilon})$,
and obtain \bs{z} by means of a differentiable transformation
depending on $\lambda, \bs{\epsilon}$ and $\bs{x}$. Specifically, we can sample
gaussian noise $\bs{\epsilon} \sim \mathcal{N}(\bs{0}, \bs{I})$, and
obtain $\bs{z}_\lambda(\bs{\epsilon},\bs{x}_i) =
\bs{\mu}_\lambda(\bs{x}_i) +
\bs{\Sigma}_\lambda(\bs{x}_i)\cdot\bs{\epsilon}$,   normally
distributed according to parameters $\bs{\mu}_\lambda$ and $\bs{\Sigma}_\lambda$.

To summarize, given the loss function
\begin{equation}\label{eq:vaeloss}
\begin{split}
\mathcal{L}(\phi, \lambda; \bs{X})  =  \sum_i &\left\{
\frac{1}{2}
    \sum_k \left(\sigma_{\lambda,k} (\bs{x}_i)  - 1 -
      \log\sigma_{\lambda,k}(\bs{x}_i)
    + \mu_{\lambda,k}(\bs{x}_i)^2
\right)\right.\\
& - \left.
\vphantom{\sum_{k} k}
  E_{\bs{\epsilon} \sim \mathcal{N}(\bs{0}, \bs{I})} \left[
    \log P_\phi\left(\bs{x}_i | \bs{z}_\lambda(\bs{\epsilon}, \bs{x}_i)\right) \right]
\right\} \, ,
\end{split}
\end{equation}
we can learn the parameters $\lambda$ and $\phi$, and consequently
the encoder $Q_\lambda(\bs{z} | \bs{x})$ and the decoder
$P_\phi(\bs{x} | \bs{z})$.

\section{Variational Autoencoders for User Preferences}\label{sec:vae}

The general framework described in the previous section can be
instantiated to the case of collaborative filtering by specifying the
\bs{x} variables, and consequently the probability density
$P_\phi(\bs{x} | \bs{z})$. We analyze some alternative models here. We
shall use the following shared notation: $u \in U =\{1,
  \ldots, M\}$ indexes a user and $i \in I = \{1, \ldots, N\}$ indexes an
item for which the user can express a preference. We model implicit
feedback, thus assuming a preference matrix $\bs{X} \in
\{0,1\}^{N\times M}$, so that $\bs{x}_u$ represents the (binary) row
with all the item preferences for user $u$. Given $\bs{x}_u$, we define
$I_u= \{i\in I | x_{u,i} = 1 \}$ (with $N_u = |I_u|$). Analogously, $U_i= |\{u\in U |
x_{u,i} = 1 \}|$ and $M_i = |U_i|$.

We also consider a precedence and temporal relationships
within \bs{X}. First of all, the preference matrix induces a natural ordering
relationship between items: $i \prec_u j$ has the meaning that
$x_{u,i} > x_{u,j}$ in the rating matrix.  Also, we assume the
existence of timing information $\bs{T}\in
{\rm I\!R}_{+}^{M\times N}\cup \{\emptyset\}$, where the term $t_{u,i}$
represents the time when $i$ was chosen by $u$ (with $t_{u,i} =
\emptyset$ if $x_{u,i} = 0$). Then, $i <_u j$ denotes that $t_{u,i} <
t_{u,j}$, With an abuse of notation, we also introduce a temporal mark
in the elements of $\bs{x}_u$: the term $\bs{x}_{u(t)}$ (with $1 \leq t
\leq N_u$) represents the $t$-th item in $I_u$ in the sorting induced
by $<_u$, whereas $\bs{x}_{u(1:t)}$ represents the sequence
$\bs{x}_{u(1)}, \ldots, \bs{x}_{u(t)}$.

\subsection{Multinomial model}
\label{sec:mvae}

The reference model is the Multinomial variational autoencoder (MVAE)
proposed in \cite{Liang:2018:VAC:3178876.3186150}. Within this
framework, for a given user $u$ it is possible to devise $\bs{x}_u$
according to the generative setting:
\begin{equation}
\begin{gathered}
\bs{z}_u \sim \mathcal{N}(\bs{0}, \bs{I}_K) \, , \qquad  \pi(\bs{z}_u)
\sim \exp\left\{f_\phi(\bs{z}_u)\right\} \, , \\
 \bs{x}_u \sim \mathit{Multi}\left(N_u,\pi(\bs{z}_u)\right) \, .
\end{gathered}
\end{equation}
The underlying ``decoder'' is modeled by
$$
\log P_\phi(\bs{x}_u |\bs{z}_u) = \sum_i x_{u,i} \log \pi_i(\bs{z}_u) \, ,
$$
thus enabling a complete specification of the overall variational
framework. Prediction for new items is accomplished by resorting to
the learned functions $f_\phi$ and $Q_\lambda$: given a user history
$\bs{x}$, we compute $\bs{z} = \bs{\mu}_\lambda(\bs{x})$ and then
devise the probabilities for the whole item set through
$\pi(\bs{z})$. Unseen items can then be ranked according to
their associated probabilities.

\subsection{Ranking Models}
\label{sec:bpr}

A different formulation is inspired by the Bayesian Personalized Ranking
(BPR) model introduced in
\cite{Rendle:2009}. Here, we explicitly model an
order among observations, so that for each triplet $u,i,j$ we can
devise whether $i \prec_u j$. In classical BPR, this is modeled
by associating a factorization rank $\bs{p}^T_u \bs{q}_i$ for each
pair $(u,i)$, and requiring that $\bs{p}^T_u \bs{q}_i < \bs{p}^T_u
\bs{q}_j$ whenever $i \prec_u j$. In a VAE setting, $\bs{p}^T_u
\bs{q}_i$ is replaced by a score based on the latent gaussian variable
$\bs{z}_{u,i}$:
\begin{equation}
\begin{gathered}
\bs{z}_{u,i} \sim \mathcal{N}(\bs{0}, \bs{I}_K) \, ,\\
i \prec_u j \sim \mathit{Bernoulli}\left(\pi(\bs{z}_{u,i},
   \bs{z}_{u,j})\right) \, ,
\end{gathered}
\end{equation}
Two different specifications for $\pi(\bs{z}_{u,i}, \bs{z}_{u,j})$ are
possible: either we let the neural network implicitly models the
comparison between $i$ and $j$
\begin{equation}\label{eq:rvaep}
\pi(\bs{z}_{u,i}, \bs{z}_{u,j})
\sim \sigma\left(f_\phi(\bs{z}_{u,i}, \bs{z}_{u,j})\right) \, ,
\end{equation}
or alternatively, we can assume that the network simply provides the
score upon which to compare:
\begin{equation}\label{eq:rvaed}
\pi(\bs{z}_{u,i}, \bs{z}_{u,j})
\sim \sigma\left(f_\phi(\bs{z}_{u,j}) -  f_\phi(\bs{z}_{u,i})\right) \, .
\end{equation}
Here, $\sigma$ represents the standard logistic function $\sigma(z) =
1/(1+\exp\{-z\})$.

These alternatives represent two different instantiations of the
Ranking Variational Autoencoder. The
difference between them is substantial at prediction time: given $u$
and a set $J \subset I$ of unseen items, in both models we need to
devise $\bs{z}_{u,i} = \bs{\mu}_\lambda(\bs{x}_u | x_{u,i} = 1)$ for
each $i\in J$. However, the instance of eq. \ref{eq:rvaep} requires to
organize the data as triplets $u,i,j$  and to
compute $P_\phi(i \prec_u j)$ for each triplet. These probabilities
can be then exploited to devise an order among the items in
$U$. Inconsistencies can arise in this setting, which can in principle
prevent to produce a global rank order.

Conversely, the instance of eq. \ref{eq:rvaed} can
directly sort items according to $f_\phi(\bs{z}_{u,i})$ and does not
produce inconsistencies, thus being better suited for ranking
items. In the following we shall only focus on this
instance, that we shall refer to as RVAE.

\subsection{Sequential Models}\label{seq:sequential}

\begin{figure*}[th!]
 \centering
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
          \includegraphics[height=1.2in]{mvae-arch}
        \caption{MVAE}
        \label{fig:mvae}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
          \includegraphics[height=1.2in]{rvae-arch}
        \caption{RVAE}
        \label{fig:rvae}
    \end{subfigure}
   ~
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
          \includegraphics[height=1.2in]{svae-arch}
        \caption{SVAE}
        \label{fig:svae}
    \end{subfigure}
    \caption{Variational architectures. Terms with
    hat represent decoding reconstruction. Dotted boxes represent neural
    layers. }
    \label{fig:architecture}
\end{figure*}

The basic framework proposed in section \ref{sec:back} can also
be exploited to model time-aware user preferences. Ideally, latent
variable modeling should be able to express temporal dynamics and hence
causalities and dependencies among preferences in a user's
history. In the following we elaborate on this intuition.

\subsubsection{Modeling History-aware User Preferences}\label{sec:svae-simple}
Within a probabilistic framework, we can model temporal dependencies by
conditioning each event to the previous events: given a sequence
$\bs{x}_{(1:T)}$, we have
$$
P\left(\bs{x}_{(1:T)}\right) = \prod_{t=0}^{T-1}
P\left(\bs{x}_{(t+1)}| \bs{x}_{(1:t)}\right) \, .
$$
This specification suggests two key aspects:
\begin{itemize}
\item There is a recurrent relationship between $\bs{x}_{(t+1)}$
  and $\bs{x}_{(1:t)}$ devised by
  $P\left(\bs{x}_{(t+1)}|\bs{x}_{(1:t)}\right)$, upon which the
  modeling can take advantage,   and
\item each time-step can be handled separately, and in particular it
  can be modeled through a conditional VAE.
\end{itemize}

Let us consider the following (simple) generative model
\begin{equation}\label{eq:svae}
\begin{gathered}
\bs{z}_{u(t)} \sim \mathcal{N}\left(\bs{0},\bs{I}_K\right) \, ,
\qquad  \pi\left(\bs{z}_{u(t)}\right)
\sim \exp\left\{f_\phi\left(\bs{z}_{u(t)}\right)\right\} \, , \\[2mm]
 \bs{x}_{u(t)} \sim \mathit{Multi}\left(1,\pi\left(\bs{z}_{u(t)}\right)\right) \, ,
\end{gathered}
\end{equation}
which results in the joint likelihood
$$
P(\bs{x}_{u(1:T)}, \bs{z}_{u(1:T)}) = \prod_t
P(\bs{x}_{u(t)} |\bs{z}_{u(t)})P(\bs{z}_{u(t)}) \, .
$$
Here, we can approximate the posterior $P(\bs{z}_{u(1:T)}|
\bs{x}_{u(1:T)})$ with the factorized proposal distribution
$$
Q_\lambda(\bs{z}_{u(1:T)}|\bs{x}_{(1:T)}) = \prod_t q_\lambda(\bs{z}_{u(t)}|\bs{x}_{(1:t-1)}) \, ,
$$
where $q_\lambda(\bs{z}_{u(t)}|\bs{x}_{(1:t-1)})$ is a gaussian
distribution whose parameters $\bs{\mu}_\lambda(t)$ and
$\bs{\sigma}_\lambda(t)$ depend upon the
current history $\bs{x}_{u(1:t-1)}$, by means of a recurrent layer $\bs{h}_t$:
\begin{equation}\label{eq:recurrent}
\begin{split}
\bs{\mu}_\lambda(t), \bs{\sigma}_\lambda(t) & =
                                                               \varphi_\lambda(\bs{h}_t)\\
\bs{h}_t & = \mathit{RNN}_\lambda\left(\bs{h}_{t-1},   \bs{x}_{u(t-1)} \right) \, .
\end{split}
\end{equation}

The resulting loss function follows directly from eq.~\ref{eq:vaeloss}:
\begin{equation*}
\begin{split}
\mathcal{L}(\phi, \lambda,\bs{X})  =  \sum_u \sum_{t=1}^{N_u}&\left\{
\frac{1}{2}
    \sum_k \left(\sigma_{\lambda,k} (t)  - 1 -
      \log\sigma_{\lambda,k}(t) + \mu_{\lambda,k}(t)^2
\right)\right.\\
& - \left.
\vphantom{\sum_{k} k}
  E_{\bs{\epsilon} \sim \mathcal{N}(\bs{0}, \bs{I})} \left[
    \log P_\phi\left(\bs{x}_{u(t)} | \bs{z}_\lambda(\bs{\epsilon}, t)\right) \right]
\right\} \, .
\end{split}
\end{equation*}

The proposal distribution introduces a
dependency of the latent variable from a recurrent layer, which allows
to recover the information from the previous history. We call this
model SVAE. Figure \ref{fig:architecture} shows the main
architectural difference with respect to the models proposed so
far. In SVAE, we can observe the recurrent relationship occurring in the layer
upon which $\bs{z}_{u(t)}$ depends.

Notably, the prediction step can be easily accomplished in a similar
way as for MVAE: given a user history $\bs{x}_{u(1:t-1)}$, we can resort
to eq.~\ref{eq:recurrent} and set $\bs{z} =
\bs{\mu}_\lambda(t)$, upon which we can devise the probability for the
$\bs{x}_{u(t)}$ by means of $\pi(\bs{z})$.

\subsubsection{A taxonomy of sequential variational autoencoders}
\label{sec:svaetaxonomy}
We already discussed that within a sequence modeling framework, the
core of the approach is on modeling $P\left(\bs{x}_{(t+1)}|
  \bs{x}_{(1:t)}\right)$ through a conditional variational autoencoder. SVAE describes
just one of several possible modeling choices.
\begin{figure*}[ht!]
 \centering
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \begin{tabular}{c}
          \includegraphics[height=1.2in]{rvae1}
          \\
          $            \scriptscriptstyle
            P\left(\bs{x}_{(1:T)}\right) = \int P\left(\bs{z}_{(1:T)}\right)\prod_{t=1}^{T-1} P_\phi\left(\bs{x}_{(t+1)}|
  \bs{z}_{(t+1)}, \bs{x}_{(1:t)}\right) \diff \bs{z}_{(1:T)}
          $
        \end{tabular}
        \caption{Single latent dependency}
        \label{fig:SLD}
    \end{subfigure}%
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \begin{tabular}{c}
          \includegraphics[height=1.2in]{rvae2}
          \\
          $
            \scriptscriptstyle
            P\left(\bs{x}_{(1:T)}\right) = \int P\left(\bs{z}_{(1:T)}\right)\prod_{t=1}^{T-1} P_\phi\left(\bs{x}_{(t+1)}|
              \bs{z}_{(1:t)},\bs{x}_{(1:t)}\right) \diff
            \bs{z}_{(1:T)}
          $
        \end{tabular}
        \caption{Multiple latent dependencies}
        \label{fig:MLD}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \begin{tabular}{c}
          \includegraphics[height=1.2in]{rvae3}
          \\
          $
            \scriptscriptstyle
            P\left(\bs{x}_{(1:T)}\right) = \int\prod_{t=1}^{T-1} P_\phi\left(\bs{x}_{(t+1)}|
              \bs{z}_{(t+1)}\right)
            P\left(\bs{z}_{(t+1)}|\bs{x}_{(1:t)}\right)\diff
            \bs{z}_{(1:T)}
          $
        \end{tabular}
        \caption{Recurrent latent dependency}
        \label{fig:RLD}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \begin{tabular}{c}
          \includegraphics[height=1.2in]{rvae4}
          \\
          $
            \scriptscriptstyle
            P\left(\bs{x}_{(1:T)}\right) = \int\prod_{t=1}^{T-1} P_\phi\left(\bs{x}_{(t+1)}|
              \bs{z}_{(1:t+1)}, \bs{x}_{(1:t)}\right)
            P\left(\bs{z}_{(t+1)}|\bs{z}_{(1:t)}, \bs{x}_{(1:t)}\right)\diff
            \bs{z}_{(1:T)}
          $
        \end{tabular}
        \caption{Global recurrent dependency}
        \label{fig:GLD}
    \end{subfigure}
    \caption{Recurrent Variational Autoencoders. Diamond boxes
      represent deterministic variables.}
    \label{fig:recvae}
\end{figure*}
In fact, four alternate formalizations can take place, as illustrated in
fig.~\ref{fig:recvae}:
\begin{itemize}
\item The simplest modeling, considers a single gaussian variable
  \bs{z} and a parameterization of the conditional distribution
$P_\phi\left(\bs{x}_{(t+1)}| \bs{z},\bs{x}_{(1:t)}\right)$ with a function
$f_\phi(\bs{z}, \bs{h}_{t+1})$, where $\bs{h}_t$ represents the hidden
state of a recurrent neural function
\begin{equation}\label{rnn-x}
  \bs{h}_t = \mathit{RNN}\left(\bs{h}_{t-1}, \bs{x}_{(t)}\right) \, .
\end{equation}
Figure \ref{fig:SLD} illustrates the graphical model and the sequence
likelihood.
\item Following \cite{BayerO14a}, we can alternatively introduce $t$
  independent gaussian variables and parameterize
  the conditional likelihood
$P_\phi\left(\bs{x}_{(t+1)}|
  \bs{z}_{(1:t)},\bs{x}_{(1:t)}\right)
$
by means of a function
$f_\phi(\bs{h}_{t+1})$, where again $\bs{h}_t$ represents the hidden
state of a recurrent neural function
\begin{equation}\label{rnn-x-z}
  \bs{h}_t = \mathit{RNN}\left(\bs{h}_{t-1}, \bs{x}_{(t)}, \bs{z}_{(t)}\right) \, .
\end{equation}
Figure \ref{fig:MLD} illustrates the graphical model and the sequence
likelihood.
\item So far, the modeling combines history and latent variables to
  define the conditional distribution. An alternative consists in
  assuming that history affects the latent variable \bs{z} instead. In
  practice, the conditional
  likelihood $P_\phi\left(\bs{x}_{(t+1)}| \bs{z}_{(t+1)}\right)$ would
  only depend on the latent variable, which exhibits a prior distribution
  $P\left(\bs{z}_{(t+1)}|\bs{x}_{(1:t)}\right)$, modeled as a gaussian
  with parameters depending on the current state $\bs{h}_t$ of the
  network, devised as in eq.~\ref{rnn-x}. The graphical model and
  sequence likelihood are shown in fig. \ref{fig:RLD} and it
  resembles the SVAE model discussed above.
\item Finally, \cite{Chung:2015:RLV:2969442.2969572} propose a comprehensive
  model, where the gaussian latent variables $t+1$ also depend on each
  other through the Markovian dependency $P\left(\bs{z}_{(t+1)}|\bs{z}_{(1:t)},
    \bs{x}_{(1:t)}\right)$. Both this dependency and
  $P_\phi\left(\bs{x}_{(t+1)}| \bs{z}_{(1:t+1)}, \bs{x}_{(1:t)}\right)$ can be specified by
  the hidden state $\bs{h}_{t}$ of a recurrent network, devised as in
  eq. \ref{rnn-x-z}.  The graphical model and
  sequence likelihood are shown in fig. \ref{fig:GLD}.
\end{itemize}

\subsubsection{Extending SVAE}\label{sec:svae-ext}
The generative model of eq. \ref{eq:svae} only focuses on the
next item in the sequence. The base model however is flexible enough
to extend its focus on the next $k$ items, regardless of the time:
$$
 \bs{x}_{u(t:t+k-1)} \sim \mathit{Multi}\left(k,\pi\left(\bs{z}_{u(t)}\right)\right),
$$
Again, the resulting joint likelihood can be modeled in different
ways.
The simplest way consists in considering $\bs{x}$ as a time-ordered multi-set,
\begin{equation}\label{eq:svae-ext-std}
  P(\bs{x}_{u(1:T)}, \bs{z}_{u(1:T)}) = \prod_t
  P(\bs{x}_{u(t:t+ k -1)} |\bs{z}_{u(t)})P(\bs{z}_{u(t)}) \, .
\end{equation}
Alternatively, we can consider the probability of an item as a mixture
relative to all the time-steps where it is considered:
\begin{equation}\label{eq:attention}
P(\bs{x}_{u(1:T)}, \bs{z}_{u(1:T)}) = \prod_t P(\bs{z}_{u(t)})
\sum_{t-k +1\leq j \leq t}P(\bs{x}_{u(t)} |\bs{z}_{u(j)}),
\end{equation}
where $P(\bs{x}_{u(t)} |\bs{z}_{u(j)})$ is the probability of
observing $\bs{x}_{u(t)}$ according to $\pi(\bs{z}_{u(j)})$.
In both cases, the variational approximation is modeled exactly as
shown above, and the only difference lies in the second component of
the loss function, which has to be adapted according to the above
equations.

There is an interesting analogy between eq. \ref{eq:attention} and the
attention mechanism \cite{Vaswani:2017}. In fact, it can be noticed in
the equation that the prediction of $\bs{x}_{u(t)}$ depends on the
latent status of the $k$ previous steps in the sequence. In practice,
this enables to capture short-term dependencies and to encapsulate
them in the same probabilistic framework by weighting the likelihood
based on $\bs{z}_{u(t-k+1)}, \ldots, \bs{z}_{u(t)}$.

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

We evaluate SVAE on some benchmark datasets, by comparing with various
baselines and the current state-of-the-art competitors, in order to
assess its capabilities in modeling preference data. Additionally, we
provide a sensitivity analysis relative to the configurations/contour
conditions upon which SVAE is better suited. The main highlight from
our experiments is that SVAE provides a huge edge over the current
state-of-the-art for the task of top-N recommendation across various
metrics.

\vspace*{-5mm}
\subsection{Datasets}
We evaluate our model along with the competitors on two popular
publicly available datasets, namely \textit{Movielens-1M} and
\textit{Netflix}.
\textrm{Movielens-1M} is a time series dataset containing
user-item ratings pairs along with the corresponding
timestamp. Since we work on implicit feedback, we binarize the
data, by considering only the user-item pairs where the rating provided
by the user was strictly greater than 3 on a range 1:5.
\textrm{Netflix} has the same data format as \textrm{Movielens-1M} and
the same technique is used to binarize the ratings. We use a subset
of the full dataset that matches the user-distribution with the
full dataset. The subset is built by stratifying users according to
their history length, and then sampling a subset of size inversely
proportional to the size of the
strata. Figure~\ref{fig:netflix_dist} compares the distributions of
the full and the sampled dataset: We can notice that the
distributions share the same shape, but in the sample users with
small history length are undersampled whereas users with large
histories are kept.

\begin{table}
\begin{center}
  \includegraphics[width=0.4\textwidth]{table-datasets}
\end{center}
\caption{Basic statistics of the datasets used in the experiments.}
\label{dataStats}
\end{table}

Table~\ref{dataStats} shows the basic statistics of the data. For
illustration purposes, we also show the basic statistics of the full
Netflix dataset. We can see that the average length of sequences in
the Netflix subset is significantly increased, as a result of
downsampling users with small history length. Also, notice that the
sampling procedure does not affect the number of items.
To preprocess the data, we
first group the interacted items for each user, and ignore the users
who have interacted with less than five items. After preprocessing, we
split the remaining users into train, validation and test sets.

\begin{figure}[h!]
\begin{minipage}[t]{0.45\linewidth}
    \includegraphics[height=1.3in]{Netflix-Full_count_dist_user}

\end{minipage}%
    \hfill%
\begin{minipage}[t]{0.45\linewidth}
    \includegraphics[height=1.3in]{Netflix-Subset_count_dist_user}

  \end{minipage}
  \vspace{-5mm}
  \caption{Comparison of the user distributions of the full Netflix
    dataset and the subsample produced by progressive stratified
    sampling.}
    \label{fig:netflix_dist}
\end{figure}

\vspace*{-5mm}
\subsection{Evaluation Metrics and protocol}
Since we are considering implicit preferences, the evaluation is
done on top-n recommendation, and it relies on the following
metrics.
\begin{description}
\item[Normalized Discounted Cumulative Gain.] Also abbreviated as
  \textit{NDCG@n}, the metric gives more weight to the relevance of
  items on top of the recommender list and is defined as
$$
  \textit{NDCG@n} = \frac{\textit{DCG@n}}{\textit{IDCG@n}}
$$
where
$$
\textit{DCG@n} =\sum_{i=1}^{n} \frac{r_{i}}{log_{2}(i+1)}  \quad\text{and}\quad
\textit{IDCG@n} = \sum_{i=1}^{|\textit{R}|} \frac{1}{log_{2}(i+1)}.
$$
Here, $r_i$ is the relevance (either 1 or 0 within the implicit
feedback scenario) of the $i$-th recommended items in the
recommendation list, and $R$ is the set of relevant items.
\item[Precision.] By defining $\textit{Hits} = \sum_i r_i$ as the number of items
  occurring in the recommendation list that were actually
  preferred by the user, we have
  $$
  \vspace*{-1mm}
  \textit{Precision}@n = \frac{\textit{Hits}@n}{n}
  $$
\item[Recall,] defined as the percentage of items actually preferred by the user
  that were present in the recommendation list:
$$
Recall@n  = \frac{\textit{Hits}@n}{|R|}
$$
\end{description}
In our experiments, we use the above metrics with two values of $n$,
respectively $10$ and $100$.

The evaluation protocol works as follows.
We partition users into training, validation and test set. The
model is trained using the full histories of the users in the training
set. During evaluation, for each user in the validation/test
set we split the \textbf{time-sorted} user history into two parts,
\textit{fold-in} and \textit{fold-out} splits. The \textit{fold-in}
split is used as a basis to learn the necessary representations and provide a
recommendation list which is then evaluated with the \textit{fold-out} split
of the user history using the metrics defined above.
We believe that this strategy is more robust when compared to
other methodologies wherein the same user can be in both the training
as well as testing sets. Table~\ref{dataStats} shows the number of
heldout users for each datasets.

It is worth noticing that Liang et. al
\cite{Liang:2018:VAC:3178876.3186150} follow a similar strategy but
with a major difference: they do not consider any sorting for user
histories. That is, for the validation/test users, the \textit{fold-in}
set doesn't precede the \textit{fold-out} with respect to time. By
contrast, we keep the
\textit{fold-in} set to be the first \textit{80\%} of the time-sorted
user history, and the last $20\%$ represents the \textit{fold-out} set. We
shall see in the following sections that this difference is substantial in the
evaluation.

\subsection{Competitors}
We compare our model with various baselines and current
state-of-the-art models including recently published neural
architectures and we now present a brief summary about our
competitors to provide a better understanding of these models.

\begin{itemize}
\item \textbf{POP} is a simple baseline where users are recommended
  the most popular items in the training set.

\item \textbf{BPR}, already mentioned in section \ref{sec:bpr}, is a
  state of the art model based on Matrix
  Factorization, which ranks items
  differently for each user \cite{Rendle:2009}.

  There is a subtle
  issue concerning BPR: by separating users on
  training/validation/test as discussed above, the latent
  representation of users in the validation/test is not
  meaningful. That, is, BPR is only capable of providing meaningful
  predictions for users that were already exploited in the training
  phase. To solve this, we extended the training set to include
  the partial history in \textit{fold-in} for each user in the
  validation/test. The evaluation still takes place on their
  corresponding \textit{fold-out} sets.

\item \textbf{FPMC} \cite{Rendle:2010} is a model which
  clubs both Matrix Factorization and Markov Chains together using
  personalized transition graphs over underlying Markov
  chains.

\item \textbf{CASER} \cite{Tang:2018}, already discussed in section
  \ref{sec:related}, is a convolutional model that uses vertical and
  horizontal convolutional layers to embed a sequence of recent items
  thereby learning sequential patterns for next-item
  recommendation. The authors have shown that this model outperforms
  other approaches based on recurrent neural network modeling, such as
  GRU4Rec. We use the implementation provided by the authors and tune
  the network by keeping the number of horizontal filters to be 16,
  and vertical filters to be 4.

\item \textbf{MVAE}, discussed in section \ref{sec:mvae}, from which
  the SVAE model draws heavily. We use the implementation provided by
  the authors, with the default hyperparameter settings.

\end{itemize}

We also include the \textbf{RVAE} model proposed in section
\ref{sec:bpr}, that we consider a baseline here. Notably, despite
being considered a simple extension of the BPR
model, RVAE relies on a neural layer for embedding users: as a
consequence,  it does a better job in ranking items for the users which the
model has never seen before, contrary to BPR. In practice, RVAE
upgrades BPR to session-based recommendations.

\subsection{Training Details}
The experiments only consider the SVAE model illustrated in subsection
\ref{sec:svae-simple} and the extensions of subsection~\ref{sec:svae-ext}. We
reserve a more detailed analysis of the extensions discussed in
subsection~\ref{sec:svaetaxonomy} to future work.
The model is trained end-to-end on the full histories of the training
users.
Model hyperparameters are set using
the evaluation metrics obtained on validation users.

The SVAE architecture includes an embedding layer of size 256, a recurrent layer
realized as a GRU with 200 cells, and two encoding layers (of size 150
and 64) and finally two decoding layers (again, of size 64 and 150).
We set the number $K$ of latent factors for the
variational autoencoder to be 64. Adam \cite{KingmaB14} was used to
optimize the loss function coupled with a weight decay of $0.01$.
As for RVAE, the architecture includes user/item embedding layers (of
size 128), two encoding layers (size 100 and 64), and a final layer
that produces the score $f_\phi(\bs{z}_{u,i})$.
Both SVAE and RVAE were implemented in PyTorch \cite{paszke2017automatic}
and trained on a single GTX 1080Ti GPU.
The source code is available
on GitHub\footnote{\url{https://github.com/noveens/svae_cf}.}.

\subsection{Results}

\begin{table*}
\begin{center}
  \includegraphics[width=0.8\linewidth]{table-results.pdf}
\caption{Results of the evaluation (in percentage). MVAE$^\ast$ considers random
  splits that disregards the temporal order of user history.
  BPR$^\ast$ relies on including the \textit{fold-in} subsequences in
  the training phase.}
\label{table:resultsTable}
\end{center}
\end{table*}

In a first set of experiments, we compare SVAE with all competitors
described above. Table~\ref{table:resultsTable} shows the results of
the comparison. SVAE consistently outperforms the competitors on both
datasets with a significant gain on all the metrics. It is important
here to highlight how the temporal \textit{fold-in}/\textit{fold-out}
split is crucial for a fair evaluation of the predictive capabilities:
MVAE was evaluated both on temporal and random split, exhibiting
totally different performances. Our interpretation is that, with
random splits, the prediction for an item is easier if the encoding
phase is aware of forthcoming items in the same user history. This
severely affects the performance and overrates the predictive
capabilities of a model: In fact, the accuracy of MVAE drops
substantially when a temporal split is considered.

By contrast, SVAE is trained to capture the actual temporal
dependencies and ultimately results in better predictive accuracy.
This is also shown in fig.~\ref{fig:seq_len_vs_ndcg}, where we show
that SVAE consistently outperforms the competitors irrespective of the
size of \textit{fold-in}. The only exception is with very short
sequences (less than 10 items), where MVAE gets better results with
respect to the sequential models. It is also worth noticing how the
performance of both sequential models tend to degrade with increasing
sequences, but SVAE maintains its advantage over CASER.

\begin{figure}[th!]
  \vspace*{-2mm}
  \centering
  \includegraphics[width=\columnwidth]{bars_dist.pdf}
  \vspace*{-8mm}
  \caption{Average NDCG@100 for MVAE, CASER \& SVAE across
    various history lengths.}
  \label{fig:seq_len_vs_ndcg}
\end{figure}

We discussed in section \ref{sec:svae-ext} how the basic SVAE framework can
be extended to focus on predicting the next $k$ items, rather then
just the next item. We analyse this capability in
fig.~\ref{fig:next_k_vs_ndcg}, where the accuracy for different values
of $k$ is considered according to the modeling in
\ref{eq:svae-ext-std}. On Movielens, the best value is achieved for
$k=4$, and acceptable values range within the interval $1-10$.

\begin{figure}[th!]
  \centering
  \includegraphics[height=1.2in]{next_k_vs_ndcg.pdf}
        \vspace*{-4mm}
    \caption{NDCG@100 values for SVAE on Movielens, across different sizes for the
      number of items forward in time to predict on.
}
    \label{fig:next_k_vs_ndcg}
\end{figure}

Finally, in fig.~\ref{fig:learningcurve} we analyse the convergence
rate of SVAE, in terms of NDCG (left y-axis and blue line) and loss
function values (right y-axis and red line). The learning phase
converges quickly and does not exhibit
overfitting: on Movielens, a stable model is reached within 8
epochs, whereas Netflix requires 13 epochs. The average runtime per
epoch is 197 seconds on Movielens and 2 hours on Netflix: in practice,
the learning scales linearly with the number of
interactions in the dataset.

\begin{figure}[th!]
  \centering
  \includegraphics[height=1.4in,width=1.6in]{learning_curve_ml.pdf}
  ~
  \includegraphics[height=1.4in,width=1.6in]{learning_curve_netflix.pdf}
  \vspace*{-4mm}
  \caption{Learning curves on validation data.}
  \label{fig:learningcurve}
\end{figure}

\section{Conclusions and future work}\label{sec:conc}

Combining the representation power of latent spaces, provided by
variational autoencoders, with the sequence modeling capabilities of
recurrent neural networks is an effective strategy to sequence
recommendation. To prove this, we devised SVAE, a simple yet robust
mathematical framework capable of modeling temporal dynamics
upon different perspectives, within the fundamentals of variational
autoencoders. The experimental evaluation highlights the
capability of SVAE to consistently outperform state-of-the-art models.

The framework proposed here is worth further extensions
that we plan to accomplish in a future work. From a conceptual point
of view, we need to perform a thorough analysis of the taxonomy defined in section
\ref{sec:svaetaxonomy}. From an architectural point of view,
the attention mechanism, outlined in
section \ref{eq:attention}, requires a better understanding and a more
detailed analysis of its possible impact in view of the recent
developments \cite{Bahdanau14} in the literature. Also, the SVAE
framework relies on recurrent networks. However, different
architectures (e.g. based on convolution \cite{Tang:2018} or translation
invariance \cite{He:2017}) are worth being investigated within a
probabilistic variational setting.


\section{1 Introduction}
\label{Introduction}
In high energy nuclear physics, productions of prompt photon and photon associated jet
with high transverse momentum are valuable observations
of the short-distance dynamics of quarks and gluons~\cite{Catani:2002ny}.
Since prompt photon is precisely calculable by perturbative quantum chromodynamics (pQCD)
at higher order corrections and carries no color charge like other gauge bosons,
it has been widely regarded as an optimal probe of the initial state of the collisions\cite{Janus:2018jto,Ru:2014yma,Ru:2015pfa,Ru:2016ifu,Ru:2016wfx},
as well as an excellent tag of inversive parton (jet) to quantify the mechanisms of
jet quenching in ultra-relativistic heavy ion collisions~\cite{Chatrchyan:2012vq,Wang:1996pe,Wang:1996yh,Albacete:2017qng,Dai:2013xca,Helenius:2014qla,Goharipour:2017uic,Goharipour:2018sip}.
In the past few years, ATLAS~\cite{Aaboud:2016sdm,Aaboud:2017cbm,ATLAS:2017ojy} and CMS
\cite{Sirunyan:2018gro,Sirunyan:2018qec,Sirunyan:2017qhf} collaborations have made lots of measurements on
isolated prompt photon and $\gamma$+jet productions in proton-proton,
proton-nucleus and nucleus-nucleus collisions.
The isolated photon's nuclear modification factor depending on photon transverse energy $\rm E_T^{\gamma}$
measured in most central PbPb collisions has large uncertainty, it also shows centrality dependence
at limited $\rm E_T^{\gamma}$ intervals~\cite{Chatrchyan:2012vq}.
And $\gamma$+jet events have been discussed to extend the study on the tomography of
quark-gluon plasma(QGP) created in heavy ion collisions.
For instance, the distribution of photon plus jet transverse momentum imbalance $(\rm x_{j\gamma}=p_T^{jet}/p_T^{\gamma})$ indicates
that a larger part of jet's initial energy would be damped and deposited in most central Pb$+$Pb collisions\cite{Sirunyan:2017qhf}.   In this study we may employ the productions of isolated photon and
photon+jet in p+A collisions to probe the initial-state cold nuclear matter (CNM) effects.

In elementary hadron-hadron collisions, with the perturbative QCD (pQCD), the cross section of leading particle (and jet) in general could be expressed as an convolution of
the parton distribution functions (PDFs), and the hard partonic cross section, and the fragmentation functions (FFs) if applicable.
The parton distribution function (PDF) for a parton i from the free proton $\rm (f_i^p(x,Q^2))$ is of nonperturbative property in the frame of
the QCD collinear factorization theorem~\cite{Collins:1989gx,Brock:1994er} and its evolution in the scale $\rm Q^2$
can be depicted as the DGLAP equations \cite{Dokshitzer:1977sg,Gribov:1972ri,Altarelli:1977zs}.
In p+A collisions, PDFs in nuclear environment should be modified due to different CNM effects, such as
shadowing, anti-shadowing, EMC effect and Fermi motion {\it etc.} \cite{deFlorian:2003qf}.  It is expected that the QCD factorization theorem may hold for nuclei as a good approximation, and we can replace
the PDF in a free proton  $\rm (f_i^p(x,Q^2))$ with the nuclear PDF
$\rm (f_i^A(x,Q^2))$ to effectively include different CNM effects to study hard processes in p+A collisions.

In past three decades, our understanding on the global fits of nuclear PDFs (nPDFs)
have been regularly enriched by the growing experimental results of
the fixed-target deeply inelastic scattering (DIS) and low-mass Drell-Yan (DY) dilepton measurements
\cite{Eskola:1998df,Eskola:2012rg,Paukkunen:2014nqa} and
the theoretical predictions from leading order (LO) up to next-to-next-to-leading order (NNLO)
\cite{Pumplin:2002vw,Gao:2013xoa,Ball:2014uwa,Harland-Lang:2014zoa}.
In the DSSZ framework~\cite{deFlorian:2011fp}, the global analysis for the nPDFs is presented as the ratio of
parton distributions in a proton of a nucleus and in the free proton,
$\rm R_i^A(x,Q^2)=f_i^A(x,Q^2)/f_i^p(x,Q^2)$ evolving at initial scale $\rm Q_0 = 1~GeV$.
The DSSZ analysis not only uses the $l^\pm/\mu-$DIS data sets and p$+$A DY data sets,
but also firstly includes the inclusive pion production in the deutron and gold collisions at PHENIX
to constrain the nuclear gluon PDF.
The EPPS16~\cite{Eskola:2016oht} is the extension of the previous EPS09~\cite{Eskola:2009uj} with
the additional experimental data from proton-lead collisions at LHC~\cite{Chatrchyan:2014hqa,Khachatryan:2015hha,Aad:2015gta} for the first time.
It offers a less biased and flavor-dependent fitting analysis for nuclear PDFs.
Following the CTEQ global PDF fitting framework~\cite{Gao:2013xoa,Kovarik:2013sya,Dulat:2016rzo},
nCTEQ15~\cite{Kusina:2015vfa} describes the nuclear dependence of nPDFs at NLO on different nuclei, including Pb.
The nIMParton~\cite{Wang:2016mzo} is a global analysis based on two data sets of nuclear DIS data,
which either only contains isospin-scalar nuclei or all nuclear data. The difference of the fitting
functions obtained by these two data sets is on the shadowing (small x region) effects.
In addition, the Fermi motion and off-shell effect, nucleon swelling, and parton-parton recombination
are taken into account together in the nIMParton framework. In the market, some other parameterizations
of nPDFs have also been proposed~\cite{AtashbarTehrani:2012xh,Hirai:2007sx,Kulagin:2004ie,Kulagin:2014vsa,AtashbarTehrani:2017rwo}, and so far
due to the lack of enough experimental data, limited constraints and large uncertainties
appear in nuclear gluon distribution, especially at small x and large $Q^2$ region, and nuclear quark
distribution at large x for all sets of nPDFs.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.35in,height=2.4in,angle=0]{atlas_pp8TeV_gamma_pt.eps}
\hspace*{-0.1in}
\caption{The cross section of isolated photon plus jet as function of $\rm p_T^{\gamma}$ in $\rm{p+p}$ collisions at
$8$~TeV and the NLO pQCD theoretical calculations(JETHOX).
}
\label{fig:ppgammapt}
\end{center}
\end{figure}
\hspace*{-1.5in}

In this work, we study the isolated prompt photon and $\gamma+$jet productions
in proton-lead collisions at LHC energy $\sqrt{\rm s_{NN}}=8.16$TeV with a NLO
pQCD program JETPHOX~\cite{Catani:2002ny,Aurenche:2006vj,Belghobsi:2009hx} with updated
proton's PDFs $-$ CT14 parametrization~\cite{Dulat:2016rzo}.
The nuclear parton distribution functions (nPDFs) parametrizations (DSSZ, EPPS16, nCTEQ15)
are performed at next-to-leading order accuracy and the nIMParton is based on Leading Order (LO)
calculation with parton-parton recombination. These four sets of nPDFs parametrizations have been
utilized in obtaining the cross sections of photon and photon associated jet in p+A collisions.
We calculate the nuclear modification factors $\rm R_{pPb}$ for isolated photon production as function of
transverse momentum $\rm p_T^{\gamma}$ and pseudo-rapidity $\eta^\gamma$ at both forward
and backward rapidity region, the forward-backward asymmetry $\rm Y_{p\rm Pb}^{asym}$
for isolated photon production as function of $\rm p_T^{\gamma}$, and the nuclear modification
factors $\rm R_{pPb}$ for $\gamma+$jet production as function of $\eta_{\gamma \rm Jet}$
at limited $\rm p_T^{\rm avg}$ intervals.
We address the quantification of dominant Bjorken x regions detected under different
specific rapidity and transverse momentum ranges.

This paper is organized as follows: in Section 2, we describe our
perturbative QCD predictions of prompt photon associated jet inclusive
cross section in proton-proton collisions at $8$TeV.
In Section 3, we discuss the nuclear modification of isolated
prompt photon productions at both forward and backward region at $8.16$TeV.
In Section 4, the cold nuclear matter effects on photon$+$jet productions are studied
at different transverse momentum intervals at $8.16$TeV.
And we give the summary in the last Section.

\section{2 Photon and Photon+Jet Productions in p+p}
\label{sec:pp}
High-$\rm p_T$ prompt photons mainly arise from two possible
mechanisms in hadronic collisions, produced directly in the hard sub-processes
referred to as "direct" photons or fragmented from energetic parton.
We consider that the next-to-leading order(NLO) inclusive cross section for the production of
prompt photon
with transverse momentum $\rm p_T^{\gamma}$ is given by the sum of the fragmentation and direct contributions,
written as~\cite{Catani:2002ny,Belghobsi:2009hx},
\begin{eqnarray}
\begin{aligned}
\sigma(p_T^{\gamma})=&\hat{\sigma}^D(p_T^{\gamma};\mu;M;M_F) \\
&+\sum_k\int_0^1\frac{dz}{z}\hat{\sigma}^F(p_T^{\gamma}/z;\mu;M;M_F)D_k^{\gamma}(z;M_F)
\end{aligned}
\label{eq:crosssection}
\end{eqnarray}
Where $\mu$ is the renormalization scale, $M$ is the initial state factorisation scale
and $M_F$ is an arbitrary final state fragmentation scale. The contribution $\hat{\sigma}^F$ denotes
the partonic cross section for producing a parton convoluted with the PDF of the incoming proton, and $D_k^{\gamma}$
is the fragmentation function of a parton $k$ (quarks,anti-quarks~and~gluon) into a photon.
$\hat{\sigma}^D$ includes the partonic cross section for producing a direct photon and the corresponding PDFs.
Experimentally, there are also secondary photons originated from hadron decay during the collisions,
therefore an isolation cut would be applied for the substantial production of photons.
A photon is isolated if the amount of deposited hadronic transverse energy
$\rm E_T$ is not more than an specific upper limit $\rm E_T^{iso}$
in a fixed radius $\rm R_{iso}=\sqrt{(\eta-\eta_\gamma)^2+(\phi-\phi_\gamma)^2}$
in pseudo-rapidity and azimuthal angle around the photon direction.
This restriction on the yields of isolated photons could not only reject the secondary decay photons,
but also reduce the contribution from fragmentation processes.
In the following we focus on the production of isolated photon and
isolated photon tagged jets in hadronic collisions.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=3.2in,angle=0]{gammaRaF34CMS01.eps}
\hspace*{-0.1in}
\caption{(Upper) A comparison between the nuclear modification ratios $\rm R_{pPb}$
for p-Pb collisions at$\sqrt{s}$ = 8.16 TeV and $\rm 3 < \eta^{\gamma} < 4$
using the nCTEQ15 , EPPS16, DSSZ and nIMParton nuclear modifications and the
CT14 free-proton PDFs. (Bottom) The corresponding average Bjorken $\rm\langle x_{Pb}\rangle$
as function of $\rm p_T^{\gamma}$.
}
\label{fig:RaF34}
\end{center}
\end{figure}
\hspace*{-1.5in}

We calculate isolated photon and jet productions in
proton-proton collisions at $\rm 8TeV$ with a next-to-leading order
pQCD program JETPHOX~\cite{Catani:2002ny,Aurenche:2006vj,Belghobsi:2009hx}
with CT14 nucleon parton distribution functions~\cite{Dulat:2016rzo}
in accordance with the ATLAS experiment~\cite{Aaboud:2016sdm}.
The isolated energy cut for a photon has been set as $\rm E_T^{iso}<6GeV$, and
the isolated cone of radius in the pseudo-rapidity and azimuthal angle plane is $\rm R_{cone}=0.4$.
Moreover photons are selected if its transverse momentum $\rm p_T^{\gamma}>130 GeV$ and $|\eta^{\gamma}|<2.37$,
except $\rm 1.37<|\eta^{\gamma}|<1.56$.
Jets are reconstructed by anti-$\rm k_t$ algorithm with cone size $\rm R=0.6$ with $\rm p_T^{jet}>100 GeV$
and $\rm|\eta^{jet}|<4.4$.
In Fig.\ref{fig:ppgammapt}, we calculate the differential cross section
$\rm d\sigma/dp_T^{\gamma}$ up to $\rm p_T^{\gamma}=1$ TeV
in proton-proton collisions at $\rm 8$ TeV, and our theoretical
predication shows good agreement with the ATLAS experimental results.

\section{3 Isolated Photon in p+Pb collisions at very forward and backward rapidity}
\label{Nuclear modification ratio for Isolated Photon in proton$+$lead collisions at very forward and backward rapidity region}
The inclusive cross section for the isolated photon production in proton-nucleus collisions
could be evaluated by using nuclear PDFs (nPDFs) as substitutes for the free-nucleon PDFs in the collinear
factorization framework as stated above, which could effectively include different CNM effects.

In our calculations, we obtain the nPDFs $\rm f_i^A(x,Q^2)$ by multiplying the CT14 parton distribution functions (PDFs)\cite{Dulat:2016rzo}
with a flavor and scale dependent factor $\rm R_i^A(x,Q^2)$ taken from four different parametrizations
DSSZ~\cite{deFlorian:2011fp}, EPPS16~\cite{Eskola:2016oht}, nCTEQ15~\cite{Kusina:2015vfa}, nIMParton~\cite{Wang:2016mzo}.
These four parametrizations for nPDFs are similar that they categorize CNM effects with Bjorken x region into shadowing,
anti-shadowing, EMC effect and so on, but differ in the specific formalisms and parameters for describing CNM effects and the input
experimental data used in global fits. DSSZ,~nCTEQ15,~EPPS16 could be convoluted in the expression for calculating the photon
production at next-to-leading order(NLO) since they are also quantitated in the NLO pQCD framework. Whereas the leading order
results for photon production are applied with nIMParton parametrization to maintain the consistency of the analysis.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=3.2in,angle=0]{gammaRaB34CMS01.eps}
\hspace*{-0.1in}
\caption{The same as Fig.\ref{fig:RaF34}, except at backward rapidity $\rm -4 < \eta^{\gamma} < -3$.
}
\label{fig:RaB34}
\end{center}
\end{figure}
\hspace*{-1.5in}
\hspace{0.7in}
\begin{figure*}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=2.2in,angle=0]{F_40_300.eps}
\includegraphics[width=3.2in,height=2.2in,angle=0]{B_40_300.eps}
\hspace*{-0.1in}
\caption{(Left)NLO fluctuations at Forward rapidity $\rm 3<\eta^{\gamma}<4$;
(Right)NLO fluctuations at Backward rapidity $-4<\eta^{\gamma}<-3$.
}
\label{fig:fbcorr}
\end{center}
\end{figure*}
\hspace*{-1.5in}

The nuclear modification factors in proton+lead collisions are defined as:
\begin{eqnarray}
\rm R_{pPb}=\frac{d\sigma^{pPb}/dp_T}{\langle N_{coll}\rangle d\sigma^{pp}/dp_T}
\label{eq:Rab}
\end{eqnarray}
with $\rm\langle N_{coll}\rangle$ representing the number of binary nucleon-nucleon collisions by the glauber model~\cite{dEnterria:2003xac}.

Now we can make our theoretical predictions for the isolated prompt photon
production in p$+$p and p$+$Pb collisions at very forward rapidity region
$\rm 3<\eta^{\gamma}<4$ at $\rm \sqrt{s_{NN}}=8.16$TeV with ATLAS isolated cuts for photons~\cite{Aaboud:2016sdm},
along with photon's transverse momentum constrained in $\rm 40 GeV<p_T^{\gamma}<300 GeV$.
We display the nuclear modification ratio $\rm R_{pPb}^{\gamma}$ as function of $\rm p_T^{\gamma}$ in upper panel of Fig.\ref{fig:RaF34}.
And the momentum fraction carried by initial parton from the incoming particle can be roughly estimated at LO as
$\rm x_{1,2}=\frac{p_T}{\sqrt{s_{NN}}}(e^{\pm y_1}+e^{\pm y_2})$, which $\rm x_1(x_p)$ is the initial parton
coming from proton on the $\rm +z$ direction, $\rm x_2(x_{pb})$ is the initial parton coming from lead on the
$\rm -z$ direction in $\rm p+Pb$ collisions and $y_{1,2}$ is the rapidity of $\gamma$ and the associated jet respectively.
The estimated average Bjorken $\rm\langle x_{Pb}\rangle$ has been defined as the events average value of Bjorken $\rm x_{Pb}$
in JetPhox simulation.
In the bottom panel of Fig.\ref{fig:RaF34}, we show the estimation of the parton's average momentum fraction off nucleus based on NLO results in JETPHOX. We have checked that $\rm\langle x_{Pb}\rangle$ for nPDFs parametrizations vary slightly from each other,
which $\rm\langle x_{Pb}\rangle$ for nIMParton can be calculated directly at LO.
We can see that the average Bjorken $\rm\langle x_{Pb}\rangle$ is lower than 0.055 at very forward rapidity region,
which represents the shadowing effect dominating the CNM effects. Moreover the average Bjorken $\rm\langle x_{Pb}\rangle$ has a linear
positive dependence with $\rm p_T^{\gamma}$ expected in its LO estimation.

In the Fig.\ref{fig:RaF34}, we see that the DSSZ's shadowing effect is unremarkable on the suppression
of isolation photon production in p$+$Pb collisions when the average Bjorken $\rm\langle x_{Pb}\rangle<0.055$.
We also notice that DSSZ's $\rm R_{pPb}^{\gamma}(p_T^{\gamma})$ shows a very weak $\rm p_T^{\gamma}$ dependence,
which means its shadowing effect is nearly independent on photon's transverse momentum in DSSZ at forward rapidity region.
Meanwhile the other three parametrizations' shadowing decrease with $\rm p_T^{\gamma}$ increasing upon $\rm 3<\eta^{\gamma}<4$.
Also we can find that the brand new nPDF parametrization nIMParton's shadowing is only weaker than nCTEQ15's in our discussion.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=2.4in,angle=0]{gammaPtasym34CMS01.eps}
\hspace*{-0.1in}
\caption{A comparison between the forward-to-backward asymmetry $\rm Y_{asym}$
for p-Pb collisions at$\sqrt{s_{NN}}$ = 8.16 TeV and $\rm 3 < \eta^{\gamma} < 4$
using the nCTEQ15 , EPPS16, DSSZ and nIMParton nuclear modifications and the
CT14 free-proton PDFs.
}
\label{fig:asym34}
\end{center}
\end{figure}
\hspace*{-1.5in}

In Fig.\ref{fig:RaB34}, similar phenomenon that the positively linear correlation
between $\rm\langle x_{Pb}\rangle$ and $\rm p_T^{\gamma}$ has been shown at backward rapidity $-4<\eta^{\gamma}<-3$.
Whereas the estimated average Bjorken $\rm\langle x_{Pb}\rangle$ ranges from 0.25 to 0.8,
which mostly correspond to EMC effect. We could go a little further to distinguish
four different parametrizations' EMC maximum from each nuclear modification factor's extreme point,
such as the nIMParton's EMC minimum appears in $\rm p_T^{\gamma}=150GeV$ and the corresponding
average Bjorken locates around $\rm\langle x_{\rm Pb}\rangle=0.5$, which is the lowest in our results.
We further investigate the correlations between $\rm x_{Pb}$ and $\rm p_T^{\gamma}$ at both forward and backward
rapidities at NLO, shown in Fig.~\ref{fig:fbcorr}.
We observe the broadening of Bjorken $\rm x_{Pb}$ at specific $\rm p_T^{\gamma}$ interval due to higher corrections,
and the spreading of  $\rm x_{Pb}$ at small $\rm p_T^{\gamma}$ is rather wider.
Also we can see a very dense statistics cluster around the low $\rm p_T^{\gamma}$ in our Monte-Carlo simulation,
because the possibility distribution of the hard sub-processes for the photon production following
double-logarithmic declining with photon's transverse momentum.

Note the nuclear modification factor is sensitive to the nucleon
PDF(nPDF) and p+p baseline~\cite{CMS:2013cka}, we may calculate the ratio
of the photon production at forward and backward rapidity, which could get rid of
the large uncertainty in free nucleon PDFs, which could be used to probe the CNM effects with less
arbitrariness\cite{Helenius:2014qla,Goharipour:2017uic,Goharipour:2018sip}.
We define the forward-backward yield asymmetry as:
\begin{eqnarray}
\rm Y_{p\rm Pb}^{asym}=\frac{d\sigma/dp_T(p+\rm Pb \rightarrow \gamma+X)|_{\eta\in[\eta_1,\eta_2]}}{d\sigma/dp_T(p+\rm Pb \rightarrow \gamma+X)|_{\eta\in[-\eta_2,-\eta_1]}}
\label{eq:asym}
\end{eqnarray}

Our predictions of the forward-to-backward yield asymmetries $\rm Y_{p\rm Pb}^{asym}$
for the isolated prompt photon production in p$+$Pb collisions at $\rm \sqrt{s}=8.16TeV$
and $\rm 3 < |\eta^{\gamma}| < 4$ are shown in Fig.\ref{fig:asym34}.
As a result of the symmetry of the colliding system, there is nearly no forward-to-backward
yield asymmetry observed in proton-proton collisions.
$\rm Y_{p\rm Pb}^{asym}$ is larger than one in p$+$Pb collisions,
which means the photon production suffering more suppression in the backward rapidity
region.  And the EMC effect reduces the photon production more effectively
than the shadowing effect does on the whole.
Besides, we notice that the value of $\rm Y_{p\rm Pb}^{asym}$ starts going down to one
with all parametrizations due to the decreasing of EMC effect at relatively high $\rm p_T^{\gamma}$.
This manifestation seems less obvious in nCTEQ15 as its EMC maximum is close to $\rm p_T^{\gamma}$'s highest boundary
and our approximate curve fitting.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=3.2in,angle=0]{gammaEtaF34CMS01.eps}
\hspace*{-0.1in}
\caption{The same as Fig.\ref{fig:RaF34}, but as a function of photon's rapidity $\eta^{\gamma}$.
}
\label{fig:ReF34}
\end{center}
\end{figure}
\hspace*{-1.5in}
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=3.2in,angle=0]{gammaEtaB34CMS01.eps}
\hspace*{-0.1in}
\caption{The same as Fig.\ref{fig:RaB34}, but as a function of photon's rapidity $\eta^{\gamma}$.
}
\label{fig:ReB34}
\end{center}
\end{figure}
\hspace*{-1.5in}

In order to further explore the impact of input nuclear modifications
on the cross section of isolated prompt photon productions in proton-nucleus collision.
We further discuss the isolated photon's nuclear modification factor $\rm R_{pPb}^{\gamma}(\eta^{\gamma})$
as a function of photon's rapidity $\eta^{\gamma}$ at both forward and
backward rapidities. The Fig.\ref{fig:ReF34} tells us that a growing suppression
on the photon productions from DSSZ, EPPS16, nIMParton and nCTEQ15 at forward pseudo-rapidity,
which quantitatively appears in accordance with the $\rm R_{pPb}^{\gamma}(p_T^{\gamma})$ at
$p_T^{\gamma}=50 GeV$ due to the highest statistics at lowest $p_T^{\gamma}$ region in Monte-Carlo
simulations exhibited above in Fig.\ref{fig:RaF34}. On the other hand, $\rm R_{pPb}^{\gamma}(\eta^{\gamma})$
shows very weak $\eta^{\gamma}$ dependence, because the variation on Bjorken $\rm x_{Pb}$ is at the magnitude
of $10^{-3}$ at region $3<\eta^{\gamma}<4$, shown in the bottom of Fig.\ref{fig:ReF34}.
Combining the results of $\rm R_{pPb}$ evolved with $p_T^{\gamma}$ and $\eta^{\gamma}$, the suppression pattern of isolated
photon could be quantitatively analyzed through $\langle x_{Pb}\rangle$ at both forward and backward rapidities.
In Fig.\ref{fig:ReB34}, the nuclear modification factors using four different nPDFs all show the nearly positive linear relation
with $\eta^{\gamma}$, which the values could match with $\rm R_{pPb}^{\gamma}(p_T^{\gamma})$'s through the average Bjorken
$\langle x_{Pb}\rangle$, shown in Fig.\ref{fig:ReB34} and Fig.\ref{fig:RaB34} respectively.
And the nCTEQ15 parametrization gives a stronger suppression than the others', which could be
confirmed in the prediction of $\rm R_{pPb}^{\gamma}(p_T^{\gamma})$ at low $\rm p_T^{\gamma}$ region in Fig.\ref{fig:RaB34}.

\hspace{0.7in}
\begin{figure*}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=2.25in,height=1.85in,angle=0]{115_4000.eps}
\includegraphics[width=2.25in,height=1.85in,angle=0]{xp_115_4000.eps}
\includegraphics[width=2.25in,height=1.85in,angle=0]{xpb_115_4000.eps}
\hspace*{-0.1in}
\caption{Correlation between $\rm x_{Pb}/x_p$(Left), $x_p$(Middle), $x_{Pb}$(Right) and $\gamma+$jet pseudo-rapidity $\rm \eta_{\gamma jet}$.
}
\label{fig:correlation1}
\end{center}
\end{figure*}
\hspace*{-1.0in}

\section{4 Isolated Photon+Jet in p+Pb}
\label{Isolated Photon associated with Jet Productions in proton$+$lead collisions}

As compared to isolated photon productions, the isolated photon associated jet production in p+A reactions has more leverage power to access CNM effects in a wider kinematic regions due to its exclusive property.
To understand the nuclear modifications for isolated prompt photon associated jet productions,
one usually estimates the momentum fractions of the initial-state partons at leading order(LO)
to evaluate the CNM effects' contribution by the final-state kinematics,
namely $\rm x_{1,2}$ defined above. In the following,
we are enlightened by the work on dijets productions in the CMS collaboration~\cite{Sirunyan:2018qel},
and provide the $\gamma$+jet pseudo-rapidity $\rm \eta_{\gamma jet}=\frac{1}{2}(\eta_\gamma+\eta_{Jet})$ distributions
of a photon tagged jet at specific range of their average transverse momentum $\rm p_T^{avg}=\frac{1}{2}(p_T^{\gamma}+p_T^{Jet})$.
As $\rm \eta_{\gamma jet}$ would be equal to $\rm \frac{1}{2}\ln(x_p/x_{Pb})$ in
the center-of-mass frame when two partons collide with each other at LO.
The NLO simulation would give rise to complicated correlations between $\rm x_{Pb}$ ($\rm x_{p}$) and $\rm \eta_{\gamma jet}$  describing the nuclear matter's influence,
shown in  Fig.~\ref{fig:correlation1}.
We also compute the correlations between $\rm x_{Pb}/x_p$ and $\rm \eta_{\gamma Jet}$
at $\rm 115GeV<p_T^{avg}<4000GeV$ interval in Fig.~\ref{fig:correlation1}. It is shown though at NLO accuracy, both distributions of $\rm x_{Pb}$ and
$\rm x_{p}$ over $\rm \eta_{\gamma jet}$ are rather wider, when the ratio $\rm x_{Pb}/x_p$ at NLO is rather narrow and centered at values at LO with very high statistics.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=2.2in,angle=0]{diff3.eps}
\hspace*{-0.1in}
\caption{The ratio of $\langle x_2\rangle/\langle x_1\rangle$ in p$+$Pb and p$+$p collisions
as a function of $\rm\eta_{\gamma jet}$, considering different CNM effects descriptions.
}
\label{fig:diff}
\end{center}
\end{figure}
\hspace*{-1.0in}

In Fig.\ref{fig:diff}, we can notice that there are tiny shifts on $\frac{<x_2>}{<x_1>}$ as a function of $\rm\eta_{\gamma jet}$
caused by nuclear matter. Based on the relation between $\rm \eta_{\gamma Jet}$ and $\rm\langle x_{Pb}\rangle$,
we could reach an assessment about different CNM effects predominate region at $\rm \eta_{\gamma Jet}$,
which the $\gamma+$jet production is sensitive to shadowing ($\rm \eta_{\gamma Jet}>1.6$),
anti-shadowing (\rm $-0.2<\eta_{\gamma Jet}<1.6$), and EMC effects ($\rm \eta_{\gamma Jet}<-0.2$)
at $\rm 150<p_T^{avg}<200$ interval, shown as the black line in Fig.\ref{fig:correlation2}.
The average Bjorken $\rm\langle x_{Pb}\rangle$ shows a nearly negative log-linear relation with
$\rm \eta_{\gamma Jet}$, and becomes globally higher when the $\rm p_T^{avg}$ interval increasing.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.2in,height=2.2in,angle=0]{xpb4.eps}
\hspace*{-0.1in}
\caption{Mean Bjorken x of the parton from the lead ion $\rm x_{Pb}$ obtained from JetPhox
as a function of $\rm\eta_{\gamma jet}$ in different $\gamma+$jet events' $\rm p_T^{avg}$ intervals.
}
\label{fig:correlation2}
\end{center}
\end{figure}
\hspace*{-1.0in}

Furthermore, we discuss the nuclear modification factor for $\gamma+$jet production
as a function of $\rm \eta_{\gamma jet}$ at $\rm 115GeV<p_T^{avg}<4000GeV$ and
$\rm 200GeV<p_T^{avg}<250GeV$. We have roughly estimated that the nuclear modification
factor is sensitive to shadowing when $\rm \eta_{\gamma jet}>1.6$
and is sensitive to EMC effect when $\rm \eta_{\gamma jet}<-0.2$.
In the total average transverse momentum interval ($\rm 115GeV<p_T^{avg}<4000GeV$),
we could find that nIMParton offers the strongest shadowing effect, and EPPS16 has the most
predominant anti-shadowing effect and EMC effect, along with DSSZ's is in between these two.
At the meantime, nCTEQ15 provides a almost constant suppression when $\rm \eta_{\gamma jet}$ varies.
When the average transverse momentum interval limited to $\rm 200GeV<p_T^{avg}<250GeV$,
nCTEQ15's shadowing starts to restore, as well as it remains flattened at small $\rm \eta_{\gamma jet}$
and large average Bjorken $\rm \langle x_{Pb}\rangle$. nIMParton's and EPPS16's shadowing becomes more damped.
And EPPS16 has more clear anti-shadowing peak as well as nIMParton and DSSZ begin to reveal anti-shadowing peak.
DSSZ provides the strongest EMC suppression at this situation. It is emphasized by leveraging the rapidity and
transverse momentum of both photon and jet with measured $\rm p_T^{avg}$ and $\rm \eta_{\gamma jet}$,
we can get access to CNM effects in a wide regime and also allocate different kinematics precisely, where
differences between varieties of nPDFs sets may be investigated more effectively.
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.3in,height=2.8in,angle=0]{ratio0.eps}
\hspace*{-0.1in}
\caption{(Upper) The nuclear modification factor for $\gamma +$ jet production as a function of $\rm\eta_{\gamma jet}$
at $\rm 115 GeV<p_T^{avg}<4000 GeV$. The NLO pQCD calculation is based on JetPhox with nCTEQ15 , EPPS16, DSSZ and nIMParton
as the nPDFs. (Bottom) The corresponding average Bjorken $\rm\langle x_{Pb}\rangle$
as function of $\rm \eta^{\gamma Jet}$.
}
\label{fig:RpA0}
\end{center}
\end{figure}
\hspace*{-1.5in}
\hspace{0.7in}
\begin{figure}[!t]
\begin{center}
\hspace*{-0.1in}
\includegraphics[width=3.3in,height=2.8in,angle=0]{ratio2.eps}
\hspace*{-0.1in}
\caption{The same as Fig.\ref{fig:RpA0},except for at $\rm 200 GeV<p_T^{avg}<250 GeV$
}
\label{fig:RpA2}
\end{center}
\end{figure}
\hspace*{-1.5in}
\section{5 Summary}
\label{Summary and Conclusions}
We calculate the productions of isolated prompt photon and $\gamma +$ jet
 in p+A with CNM effects from four sets of nuclear parton
distribution functions (nPDFs) parametrizations, {\it i.e.} DSSZ, EPPS16, nCTEQ15, nIMParton,
by utilizing the NLO pQCD approach(nIMParton at LO) at LHC $8.16$ TeV. We present the nuclear modification ratio of
isolated prompt photon $\rm R_{pA}^{\gamma}(p_T^{\gamma})$ and $\rm R_{pA}^{\gamma}(\eta^{\gamma})$,
and find shadowing effect and EMC effects dominating at very forward and backward rapidity region respectively.
And the rapidity dependence of prompt photon's nuclear modification ratio
shows weak rapidity dependence at forward region and varies linearly at backward rapidity region.
The production of the isolated photon associated with jet gives us the leveraged power to study
the tomography of cold nuclear matter.
And the CNM effects of $\gamma$+jet productions could be intuitively presented at
constraint rapidity $\rm\eta_{\gamma jet}$ and average transverse momentum $\rm p_T^{avg}$ region in our discussion.
In terms of sensitivity, comparisons of four different nPDFs parametrization's cold nuclear matter
contribution have been exhibited, and nCTEQ15 shows peculiar results than the others,
which could extend our understanding on the constraints of different nPDFs descriptions.
It is noted in our work the nPDF parametric form (nIMParton) proposed by Institute of Modern Physics in China
has been applied for the first time to investigate hard processes in ultra-relativistic heavy-ion collisions,
with a comparison again other three mainstream nPDF groups' predictions.

{\bf Acknowledgments:} This research is supported by the NSFC of China with Project
Nos. 11435004, 11322546, 11805167
and partly supported by China University of Geosciences
(Wuhan) (No. 162301182691). .

\vspace*{-.6cm}