\section{Introduction}

The barycenter in the space of probability measures equipped with the Wasserstein distance, first
introduced by \cite{agueh_carlier_wbary} gained a lot of popularity in recent years.
The most striking property of the Wasserstein space is probably how its geodesics can be interpreted as a displacement of
particles in $\RR^d$, allowing the Wasserstein barycenter of measures to take into account the geometry of the
underlying space $\RR^d$.
Thanks to its geometric interpretation, it proved useful in a variety of applications.
For instance, it has been applied in image processing:
\cite{julien_texture_mixing} developed an algorithm for
texture synthesis built around this measure, \cite{barre_atm_wbary} used the
Wasserstein barycenter of images to improve the precision of gas emission sourcing,
and it played a central role in \cite{gramfort_foptand} to compare brain image data.
In the field of fairness, \cite{barrio_fairness} and \cite{legouic_fairness} leveraged
the distribution to tackle the problem of fairness in regression and classification problems.
Barycenters also found applications in Bayesian inference; \cite{srivastava2015wasp,bayes_inference_example}
proposed a method to accelerate the computation of posterior distributions for Bayesian inference
using the Wasserstein barycenter.

More formally, let us denote by
$(P_2(\RR^d),W_2)$ the set of all measures defined on $\RR^d$ with finite second order moment, endowed with the Wasserstein distance
\begin{equation*}
    W_2 : (\mu, \nu) \longmapsto
        \sqrt{\inf_{\gamma} \int_{\RR^d \times \RR^d} \nabs{x-y}^2 \Di\gamma(x,y)} \,,
\end{equation*}
where the infimum is taken over the set of couplings of $\mu, \nu$.
A Wasserstein barycenter $b$ of $\mu_1, \dots, \mu_n$ with weights $\lambda_1, \dots, \lambda_n$
is a minimizer of the map $\nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$,
\ie the measure $b$ corresponds to a Fréchet mean of $\mu_1, \dots, \mu_n$ in the Wasserstein
space $(P_2(\RR^d), W_2)$.
Equivalently, the multimarginal formulation of the barycenter
problem asserts that $b$ is obtained by pushing forward the minimizer $\gamma^\star$ of
a functional $G$ (see equation \eqref{eqn:G_functional}) defined over the set of couplings of the
$\mu_i$'s. Under mild conditions on the $\mu_i$'s,
both problems admit unique solutions, and yield the same measure.
Section \ref{sect:notions} provides more details on the Wasserstein barycenter.

Most methods in the literature propose to estimate the Wasserstein barycenter with a discrete measure \eg
\cite{cuturi2014fast,solomon2015convolutional,benamou2015iterative}.
However, this approach does not scale well with the dimension as noted by \cite{altschuler2021wasserstein}.
Indeed, they showed that computing the infimum of $\nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$ when the
$\mu_i$ are discrete measures is a NP-hard problem in the dimension, even
when only approximate solutions are acceptable.
We can, however, avoid estimating the density altogether
and focus on generating samples distributed according to the barycenter of known measures.
Given the broad applicability of the Wasserstein barycenter and of sampling techniques in general,
we believe that such sampling procedures should be part of the statistician's toolbox.

We are motivated by the success of sampling methods in high dimensional settings, which is
due to the possibility of integrating functions against the target measure without resorting to discretization on a large grid.
Markov chain methods have become popular tools to this end.
The most well known among them are probably Hamiltonian Monte Carlo methods, variants of the
Metropolis-Hastings algorithm relying on simulations to generate diverse samples from a distribution
(\cite[Chapter 11]{bishop_prml}), and Langevin diffusion
(\cite[Chapter 4]{pavliotis_stochastic_processes_and_apps_2015}), which transports
points along random trajectories to redistribute them according to the target measure. Such transportation
methods are computationally attractive since one only needs to store minimal information about how to
move the particles at each iteration.

In their celebrated work, \cite{jko} showed that the marginals of Langevin
diffusion are distributed according to gradient flows in the Wasserstein space of the Kullback-Leibler divergence with
respect to the stationary measure. This insight, studied rigorously in \cite{ags}, brought to light a
connection between sampling and optimization which added a new perspective to the study of Monte Carlo
algorithms (\eg \cite{vempala_rates_ula} analyze the Unadjusted Langevin Algorithm and
\cite{chewi2020exponential} analyze the Mirror Langevin diffusion in this framework).
Inspired by this insight, we aim to minimize the multimarginal
formulation $G$ with gradient descent in the Wasserstein space. This procedure, like
Langevin diffusion, iteratively redistributes randomly initialized points to produce a sample from
the barycenter.
We therefore recover an approximation of the optimal
coupling through the samples, which can be desirable in applications.

However since the barycenter is defined by a constrained optimization problem,
one cannot expect that the constraints will remain satisfied along the gradient flow.
This issue of constraints has been tackled in a variety of ways in the literature. A natural
solution consists in restricting the domain of possible directions at each iteration,
which hints to the Frank-Wolfe algorithm, \ie to choose the steepest descent direction
in the admissible set.
This approach was studied in \cite{luise2019sinkhorn} within the
context of regularized optimal transport, where the regularized Wasserstein barycenter (also known
as the Sinkhorn barycenter) minimizes a sum of Sinkhorn divergences.
When the admissible set of directions forms a Reproducing Kernel Hilbert Space (RKHS) with suitable kernel,
\cite{sinkhorn_bary} propose to generate samples distributed
according to the Sinkhorn barycenter by iterative pushforward of an initial measure with the map
$\id_{\RR^d} - h \cdot d$, where $h$ is a step size and $d$ is the direction of steepest descent in the RKHS.
Both methods operate on discrete measures and the authors prove that the
continuous measure is recovered as a (weak) limit when the number of samples increases to ensure consistency.
We refer to \cite{peyre_cot} for details about regularized optimal transport.
We note that in some settings the issue of constraints can be addressed more easily, for example in the Bures-Wasserstein
manifold (\ie the subspace of Gaussian measures in the Wasserstein space)
where the barycenter problem reduces to a finite dimensional optimization problem.
\cite{chewi_bures} propose a
gradient descent algorithm for the original formulation of the barycenter problem.
However, we obviously cannot hope that similar properties hold for arbitrary continuous measures.

In this paper, we perform gradient descent on a penalized functional
$F^\al : P_2((\RR^d)^n) \to \RR^+$ obtained from $G$ by adding a penalization term to control the distance
between the coupling marginals and the $\mu_i$'s.
We weigh the penalization with a coefficient $\al$ and control the induced error with $\al$. Taking advantage of
the differential calculus on the Wasserstein space, we can define a gradient flow for $F^\al$.
To implement this procedure, we focus on the popular \ssf{SVGD} algorithm introduced in \cite{liu_wang_svgd}
to approximate the gradient of the penalization; technical details are given in Section \ref{sect:simulations}.

\begin{figure}[ht]

    \centering
    \vspace*{-1em}
    \includegraphics[ width=\textwidth]{Images/final_final_arctan.png}
    \vspace{-2em}
    \caption{\small Top plot: 200 samples drawn from $\bary{\prN(0,1), \arctan \pf \prN(0,1)}$ using
    \ssf{BARYGD} with a Gaussian kernel and initial step size $h = 10^{-3}$, and the resulting kernel
    density estimate. Bottom plot: Trajectories of the samples evolving over 220 iterations, from green to
    red samples. We can see the $\arctan$ transport map is approximately reconstructed.}
    \label{fig:norm_arctan}
\end{figure}

\paragraph{Contributions.}

Inspired by the work of Jordan, Kinderlehrer and Otto cited above, we introduce a new sampling
algorithm, called \ssf{BARYGD},
that performs kernelized gradient descent on a well chosen functional $F^\al$ built as a penalized version of
the now classical multimarginal formulation of the barycenter problem introduced by \cite{agueh_carlier_wbary}.
We show that our method is consistent: with fixed penalization strength $\al$, if gradient descent yields
a coupling $\geal$ such that $F^\al(\geal) - \inf F^\al \leq \epsilon$, then we get a quantitative bound on
the Wasserstein distance between the approximate barycenter obtained from $\geal$ and the true barycenter
obtained from the optimal coupling $\gamma^\star$.
This bound vanishes when $\epsilon \to 0$ and $\al \to \infty$.
As a consequence, we show that the minimizers of $F^\al$ converge to
$\gamma^\star$ when $\al \to \infty$ and quantify the rate of convergence. Furthermore, we
perform numerical experiments with the algorithm in several settings (see Figures \ref{fig:norm_arctan} and \ref{fig:norm}).

\paragraph{Organization.}

In Section \ref{sect:notions}, we describe the Wasserstein barycenter in more details.
We present and analyze our method in Section \ref{sect:algorithm}.
Section \ref{sect:simulations} is devoted to implementation details. In Section
\ref{sect:discussion}, we discuss some open questions.

\paragraph{Notation.}

\emph{All probability measures are assumed to be absolutely
continuous with respect to the Lebesgue measure and, with a slight abuse of notation, we
will identify measures with their densities.}
We write $X\sim \nu$ when random variable $X$ has distribution $\nu$.
For $\mu, \nu \in P_2(\RR^d)$ and a measurable map $f : \RR^d \to \RR^d$,
we write $\nu = f \pf \mu$ if $f(X) \sim \nu$ when $X \sim \mu$.
The set of couplings of $\mu$ and $\nu$, \ie probability measures on $\RR^d \times \RR^d$ with marginals
$\mu$ and $\nu$, is denoted by $\Pi(\mu, \nu)$.
We denote marginals by indices so $\gamma_1 = \mu$ and $\gamma_2 = \nu$ for $\gamma \in \Pi(\mu, \nu)$;
the notation extends naturally to couplings of $n$ measures.
We denote Euclidean gradients by $\nabla$
and write $\nabla_W$ for gradients in the Wasserstein space. Finally, we write indifferently $c_t = c(t)$
and $c' = \dot{c} = \partial_t c$ when $c : [a,b] \to X$ is a curve in some space $X$.
\section{Wasserstein barycenters}
\label{sect:notions}

This section briefly presents background notions on the Wasserstein space and its barycenters that were first
introduced in the seminal paper of \cite{agueh_carlier_wbary}.
We refer to \cite{legouic_existence} for questions of existence and stability of the barycenter.

The Wasserstein space over $\RR^d$ is defined as the set $P_2(\RR^d)$ of probability measures over $\RR^d$
with finite second order moment, endowed with the distance $W_2$ defined by
\[
    W_2^2(\mu_0,\mu_1)
        = \inf_{\gamma \in \Pi(\mu_0, \mu_1)} \int \nabs{x-y}^2\Di\gamma(x,y), \quad \forall \mu_0, \mu_1 \in P_2(\RR^d)
\]
where the infimum is taken over the set $\Pi(\mu_0, \mu_1)$ of all couplings $\gamma$ of $\mu_0$ and $\mu_1$.
The Wasserstein space is a geodesic space: for every $\mu_0$ and $\mu_1$ in $P_2(\RR^d)$ there exists
a path $(\mu_t)_{t\in [0,1]}$ such that
\[
    W_2(\mu_s,\mu_t) = |s-t| \, W_2(\mu_0,\mu_1), \quad \forall s,t \in [0,1].
\]
Such a path $(\mu_t)_{t\in[0,1]}$ is called a constant-speed \emph{geodesic} between its end points $\mu_0$ and $\mu_1$
--- we will often shorten it to \emph{geodesic}.
A functional $G$ defined over $P_2(\RR^d)$ is said to be \emph{geodesically convex} if it is convex along every geodesic.

Let $\mu_1,\dots,\mu_n\in P_2(\RR^d)$ and let $\lambda_1,\dots,\lambda_n$ be strictly positive weights
such that $\sum \lambda_i = 1$.
A Wasserstein barycenter of $(\mu_i)$ with weights $(\lambda_i)$ is defined as any measure
\begin{equation}
    \label{eqn:bary_def}
    b \in \argmin_{\nu \in P_2(\RR^d)} \sum \lambda_i W_2^2(\mu_i, \nu),
\end{equation}
\ie a barycenter corresponds to a Fréchet mean of the $\mu_i$ in the metric space $(P_2(\RR^d), W_2)$.
Whenever one of $\mu_1, \dots, \mu_n$ has density with respect to the Lebesgue measure, the barycenter is
unique and always defined. We assume throughout that the $\mu_i$ are absolutely continuous with
respect to the Lebesgue measure on $\RR^d$.

The minimization problem \eqref{eqn:bary_def} is equivalent to the following multimarginal problem.
Let
$T(x) = \sum \lambda_i x_i$ for any $x \in (\RR^d)^n$ and recall that $\Pi(\mu_1, \dots, \mu_n)$ denotes the set of all couplings of
the $\mu_i$.
The infimum of the functional
\begin{equation}
    \label{eqn:G_functional}
    G : \gamma \longmapsto \int \sum \lambda_i \nabs{x_i - T(x)}^2 \Di\gamma(x)
\end{equation}
over $\Pi(\mu_1, \dots, \mu_n)$ is achieved for a coupling $\gamma^\star$ such that
\begin{subequations}
    \begin{align}
        & b = \bary{\gamma^\star} \label{eqn:bary_pf}, \\
        & G(\gamma^\star) = \sum \lambda_i W_2^2(\mu_i, \bary{\gamma^\star}) \label{eqn:multmargform}.
    \end{align}
\end{subequations}
Moreover, if one of the $\mu_i$'s is absolutely continuous with respect to the Lebesgue measure, then $\gamma^\star$ is unique.
Note that by introducing the probability measure $P = \sum \lambda_i \delta_{\mu_i}$ on $P_2(\RR^d)$, the right hand side
of \eqref{eqn:multmargform} corresponds to the variance of $P$.
With this notations, $b$ minimizes
the variance functional $\nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$ of $P$ and
$b$ corresponds to the barycenter of $P$.

In the next section, we leverage the smooth structure of the Wasserstein space in order to develop an optimization scheme
for $G$ that will be the base of our sampling algorithm.

\section{Sampling with gradient flows}
\label{sect:algorithm}
This section defines the flow gradient of a functional on the Wasserstein space.
Some useful details of this technique introduced in \cite{jko} are gathered
in Appendix~\ref{sect:sosws}.

\subsection{Definition of the problem}

We aim to optimize the functional $G$ defined in Section \ref{sect:notions}
with a gradient flow on the Wasserstein space. For simplicity, we first define the cost function $c$ by
\begin{equation*}
    c : x \mapsto \sum \lambda_i \nabs{x_i - T(x)}^2
\end{equation*}
then, for any $\gamma \in P_2((\RR^d)^n)$, we have
\begin{equation*}
    G(\gamma) = \int c \Di\gamma.
\end{equation*}

The set of marginal constraints $\Pi(\mu_1, \dots, \mu_n)$ is convex
in $L^2((\RR^d)^n)$ but is not geodesically convex in the Wasserstein space, which makes the multimarginal
problem a difficult non-convex optimization problem on $P_2(\RR^d)$.
Since a gradient descent scheme on $G$ will inevitably leave the constraint set $\Pi(\mu_1, \dots, \mu_n)$, we modify the problem by enforcing
the constraints with a penalization term that ensures that the marginals are close to the $\mu_i$'s in $\Penalization$ divergence.
We recall that the $\Penalization$ divergence between two measures $\mu, \nu$ is given by
\[
    \Penalization_\nu(\mu) = \int \bp{\frac{\mu}{\nu} - 1}^2 \Di\nu
        \textrm{  if } \mu \ll \nu \quad\textrm{and}\quad
    \Penalization_\nu(\mu) = +\infty \textrm{  otherwise}.
\]
The problem thus becomes to minimize the functional $F^\al$ defined by
\begin{equation*}
    \gamma \mapsto F^\al(\gamma) := \int c \Di\gamma + \al \sum \lambda_i \Penalization_{\mu_i}(\gamma_i),
\end{equation*}
for $\al\in[0,\infty]$ (we use the convention $0\times\infty=0$).

For $\al$ large enough, the minimizer of $F^\al$ should be close to the minimizer of $G$ with the marginal
constraints --- which corresponds also to the minimizer of $F^\infty$.
The following proposition ensures that the minimizer exists.
\begin{proposition}
    \label{prop:barygd_falminexist}
    Suppose at least one of the $\mu_i$'s is absolutely continuous.
    Then, for any $\al > 0$, the functional $F^\al$ admits at least one
    minimizer in $P_2((\RR^d)^n)$.
    Moreover, this minimizer is absolutely continuous with respect to the Lebesgue measure.
\end{proposition}

Given a minimizer $\gamma^\al$ of $F^\al$, the measure $\bary{\gamma^\al}$, dubbed the associated barycenter to $\gamma^\al$,
should be close to the barycenter $b$ of $P$ thanks to \eqref{eqn:bary_pf} and \eqref{eqn:multmargform}.
Our next result quantifies this proximity with respect to $\al$ under extra assumptions on $\mu_1, \dots, \mu_n$.

We first assume that $\mu_1, \dots, \mu_n$ satisfy the Poincaré inequality with a strictly positive constant $\Ciso$, \ie that
for any $i \in \eint{1,n}$ and any Lipschitz $f \in L^2(\mu_i)$, we have
\begin{equation}
    \label{def:logsobolev}
    \nabs{f}_{L^2(\mu_i)}^2 \leq \Ciso \nabs{\nabla f}_{L^2(\mu_i)}^2
\end{equation}
where $\nabla f$ is defined Lebesgue-almost everywhere.
Such inequalities are very common in the sampling setting.
We also assume  that $P = \sum \lambda_i \delta_{\mu_i}$ satisfies a variance inequality with constant $k$,
\ie that there exists a strictly positive constant $k$ such that for any
$\nu \in P_2(\RR^d)$, it holds
\begin{equation*}
    k W_2^2(\nu, b_P) \leq B(\nu) - B(b_P),
\end{equation*}
where $b$ is the barycenter of $\mu_1, \dots, \mu_n$, and $B : \nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$
is the functional appearing in the original formulation of the barycenter problem.

Such inequalities were introduced by \cite{sturm_varineq} to describe curvature properties of geodesic spaces and have played
a central role in understanding the behavior of the empirical Wasserstein barycenter
(see \cite{ahidar_coutrix_rates,legouicFastConvergenceEmpirical2019}). We refer to Appendix \ref{sect:proofs} for
details on variance inequalities.
We can now state the following result.

\begin{proposition}
    \label{prop:barygd_barstab}
    Let $\al \geq 1$ and $\varepsilon \leq 1$.
    Suppose $\mu_1, \dots, \mu_n$ satisfy a Poincar\'e inequality with constant $\Ciso$ and that
    $\sum \lambda_i \delta_{\mu_i}$ satisfies a variance
    inequality with constant $\Cvar$. Then there exists a strictly positive constant $C$, only depending
    on the variance $\sigma^2 = \int c \Di\gamma^\star$ and $\Ciso$, such that for any $\geal \in P_2((\RR^d)^n)$
    satisfying $F^\al(\geal) - \inf F^\al \leq \epsilon$, it holds
    \begin{equation}
        \label{ieqn:barygd_stab_wbars_bound}
        \frac{\Cvar}{4} W_2^2(\bary{\geal}, \bary{\gamma^\star}) \leq \epsilon + \frac{C}{\sqrt{\al}}.
    \end{equation}
\end{proposition}

Proposition \ref{prop:barygd_barstab} shows that any approximate minimizer of $F^\al$ produces a barycenter
which distance to the desired barycenter $b$ is of order $O(1 / \sqrt{\al})$ and can therefore be controlled
to arbitrary precision with $\al$.

The fact that the minimization of $F^\al$ is carried out without constraints on the marginals of $\gamma$ allows
to use classical optimization techniques.
The next section is devoted to the introduction of the gradient flow on the Wasserstein space of $F^\al$ that will be the
backbone of our algorithm.
\subsection{Minimizing scheme}

We first briefly recall that the Wasserstein gradient of a functional $F : P_2(\RR^d) \to \RR$ at
a point $\nu \in P_2(\RR^d)$ is a vector field on $\RR^d$ given by the Euclidean gradient of its first variation $f$, \ie
\begin{equation}
    \label{eqn:ngrad_W}
    \nabla_W F(\nu) := x\mapsto\nabla f(x) \in L^2(\RR^d,\RR^d;\nu)
\end{equation}
where $f$ is such that for any signed measure $\delta$ with $\int \!\Di\delta = 0$, we have
\begin{equation*}
    F(\nu + \epsilon \delta) = F(\nu) + \epsilon \int f \Di\delta + o(\epsilon).
\end{equation*}
For instance, the Wasserstein gradient of $\mu \mapsto \Penalization_{\nu}(\mu)$ is given by
\[
    \nabla_W \Penalization_{\nu}(\mu) = 2 \nabla \frac{\mu}{\nu}.
\]
Intuitively, when a probability measure is seen as a large number of particles, the Wasserstein
gradient of a function defines a vector field along which each particle should be moved in order to locally maximize the function.
More details about the differential calculus on the Wasserstein space are provided in Appendix \ref{sect:sosws}.

To sample from the barycenter of $P=\sum \lambda_i\, \delta_{\mu_i}$, we consider the flow described by
\begin{equation}
    \label{eq:algo_cont}
    \begin{cases}
        \partial_t X_t = - \nabla_W F^\al(\gamma(t))(X_t) \\
        X(0) = X_0 \sim \gamma(0),
    \end{cases}
\end{equation}
where $\gamma(t)$ denotes the distribution of $X_t$ at time $t > 0$.
Recall that for a measure $\gamma \in P_2((\RR^d)^n)$, we write $\gamma_i$ for its $i$-th marginal.
Using \eqref{eqn:ngrad_W}, the Wasserstein gradient of $F^\al$ is given by
\begin{equation*}
    \nabla_W F^\al(\gamma)(x)
        = \nabla c(x) + 2 \al \sum \lambda_i \nabla \frac{\gamma_i}{\mu_i}(x),
\end{equation*}
for any $\gamma\in P_2((\RR^d)^n)$ and $x\in\RR^d$.

The explicit Euler scheme for \eqref{eq:algo_cont} provides an approximation of $(X_{t_m})_m$ for $t_{m+1}-t_m=h_m$, with
\begin{equation}
    \label{eqn:bad_euler}
    X_{m+1} = X_m - h_m \left[\nabla c(X_m) + 2 \al \sum \lambda_i \nabla \frac{\gamma_i(t_m)}{\mu_i}(X_m)\right],
\end{equation}
where $\gamma(t_m)$ is the distribution of $X_m$ and $h_m$ is a given step size.

However, note that since $\gamma(t)$ is unknown, this discretization scheme is not implementable as it is.
In order to implement it, we use a kernel to approximate the Wasserstein gradient
of the penalization term in $F^\al$ as is done in \cite{liu_wang_svgd, chewi_legouic_lu_maunu_rigollet_lawgd}.
The next section provides more details.

\section{Implementation}
\label{sect:simulations}

As noted in Section~\ref{sect:algorithm}, Algorithm \eqref{eqn:bad_euler} cannot be implemented directly
since there is no canonical way to recover the unknown distribution $\gamma(t_m)$ of the current
state $X_{t_m}$ from the mere knowledge of $X_m$. In this section, we present a concrete implementation of
the scheme \eqref{eqn:bad_euler}.

\paragraph{The \ssf{SVGD} algorithm.}

We consider the popular \ssf{SVGD}
algorithm introduced in \cite{liu_wang_svgd} which uses an explicit kernel $\ksvgd$ independent
of the target measure.
The kernel $\ksvgd$ can be chosen to have a convenient
analytical form (\eg a Gaussian kernel) amenable to direct evaluation. Let
$\pi \propto e^{-V}$, with $V \in \prC^1(\RR^d)$, be the target measure.
\ssf{SVGD} aims to minimize the $\Penalization_\pi$ divergence by transporting randomly initialized points
along the flow described by
\begin{equation}
    \label{eqn:svgd_flow}
    \partial_t X_t = - \iosvgd_{\pi} \nabla \Penalization_{\pi}(\mu_t)(X_t),
\end{equation}
where $X_t \sim \mu_t$ and $\iosvgd_{\pi} : f \mapsto \int \ksvgd(\scdot, x) f(x) \Di\pi(x)$ is the integral operator
associated with $\ksvgd$. Note that integration by parts yields
\[
    \iosvgd_{\pi} \nabla \Penalization_{\pi}(\mu_t)(X_t)
        = \int \bb{\ksvgd(X_t, x) \, \nabla V(x)  - \nabla_2 \ksvgd(X_t, x)} \Di\mu_t(x).
\]
Approximating the integrals by averaging with samples $(X^i_t)_{i=1,\dots,N}$ distributed according to $\mu_t$,
one gets the \ssf{SVGD} algorithm
\begin{equation}
    \label{eqn:svgd_iter}
    X_{t+1}^i = X_t^i +
        \frac{h_t}{N} \sum_{j=1}^N  \Big\{ -\ksvgd(X_t^i, X_t^j)\, \nabla V(X_t^j)
            + \nabla_2 \ksvgd(X_t^i, X_t^j) \Big\}
\end{equation}
where the initial points $X_0^1, \dots, X_0^N$ are chosen at random.

\paragraph{Proposed algorithm.}

To minimize $F^\al$, we implement \eqref{eqn:bad_euler}, replacing each gradient of
the penalization sum, $\sum_{i=1}^n\lambda_i\Penalization_{\mu_i}$, with the \ssf{SVGD} iteration for each marginal.
Thus, we couple $n$ batches of $N$ points and generate $N$ samples approximately distributed according to the
barycenter $\bary{\gamma^\al}$ of a coupling $\gamma^\al$ minimizing $F^\al$. Our iteration takes the form
\begin{equation*}

    X_{t+1}^{i,j}
        = X_t^{i,j} - h_t \nabla_j c(X_t)

            + \frac{\al h_t \lambda_j}{N} \sum_{\ell=1}^N
                \bb{ \ksvgd(X_t^{i,j}, X_t^{\ell,j}) \, \nabla \log \mu_j(X_t^{\ell,j}) + \nabla_2 \ksvgd(X_t^{i,j}, X_t^{\ell,j}) }
\end{equation*}
where $j \in \eint{1,n}$ is the index of marginal $\mu_j$ and $X_t^{i,j}$, $i \in \eint{1,N}$, stands for the
$i$-th sample associated with the $j$-th marginal.

\paragraph{Numerical details.}

All the simulations are carried out with a Gaussian kernel
\begin{equation*}
    \ksvgd(x,y) = \exp\bp{-\nabs{x-y}^2} \quad \forall x, y \in \RR.
\end{equation*}
The step size $h_t$ as well as $\al = \al_t$ vary over iterations, starting from $h_0=0.1$
and $\al_0=1000$. We first keep the parameters fixed and let gradient descent concentrate the particles on
the graph of an increasing function before doubling the penalization strength while keeping $h_t \cdot \al_t$
constant. Thus, we enforce the marginal constraints once the coupling corresponds to an optimal transport.
In Figure \ref{fig:norm}, we the step-size evolves according to AdaGrad (see \cite{duchi_adagrad}).
We give further experimental results in Appendix \ref{sect:app_simuls}.

\section{Conclusion}
\label{sect:discussion}

We have presented a new algorithm to sample from the Wasserstein barycenter of a finite set of measures.
Our numerical experiments illustrate its convergence in several examples.
We also proved a theoretical bound for samples drawn from almost minimizers of $F^\al$, the gradient
flow of which inspired our algorithm.
It remains to formally prove that this algorithm, a kernelized gradient descent of $F^\al$, converges to a
minimum of $F^\al$ and bound the rate of convergence.
On another note, it would be interesting to import other common tools of optimization to the Wasserstein setting
in order to derive a sampling algorithm by optimizing $F^\al$.

\begin{figure}[ht!]

    \centering
    \hspace*{-2em}
    \includegraphics[width=0.8\textwidth]{Images/svgd_norms.png}
    \includegraphics[width=0.8\textwidth]{Images/distances_norms.png}

    \caption{\small \textit{Top figure}: Top plot: 150 samples drawn from the barycenter of two unit variance normal
        distributions (shown in light gray) using our algorithm with Gaussian kernel with $\al = 10^3$ and $h_m = 10^{-4}$,
        and the resulting kernel density estimator.
        Bottom plot: Trajectories of the samples evolving over 1000 iterations, starting from green to red sample.
        The transport map is almost perfectly reconstructed.
        \textit{Bottom figure}: the Wasserstein distance to the true barycenter over iterations plotted in $\log-\log$ scale.
        From the graph, we read that \ssf{BARYGD} performs as $O(1/\sqrt{m})$ in the number of iterations $m$.
    }
    \label{fig:norm}
\end{figure}

\titleformat{\section}[block]{\sffamily\textbf{\MakeUppercase{#1}}}{}{1em}{}

\section*{Acknowledgments}

We thank Philippe Rigollet for useful discussions.
Thibaut Le Gouic was supported by NSF award IIS-1838071.

\clearpage
\printbibliography

\newpage
\onecolumn


\section{Introduction}

With recent advances of vision technologies on smartphones and Head Mounted Devices such as Facebook Oculus~\cite{fb_oculus}, Samsung Gear VR~\cite{samsung_gear}, Microsoft Hololens~\cite{MS_hololens}, VR/\ang{360} video streaming is becoming increasingly popular. The major service providers such as YouTube and Facebook have already started streaming \ang{360} videos through optimized spherical video players~\cite{fb_360,yt_360}.

\ang{360} video frames are typically 4-6 times larger than normal videos ~\cite{qian2018flare} and it can grow up to 80 times~\cite{guan_pano_2019}, requiring more resources at every entity involved in the process, i.e. servers, networks and end-user devices. Viewport-adaptive streaming, \textit{a.k.a.} Field-of-View (FoV) aware streaming, has emerged as an efficient approach for \ang{360} video streaming \cite{qian_optimizing_2016, qian2018flare,he2018rubiks}. FoV-aware streaming is expected to significantly reduce the amount of data that has to be transferred by predicting the next viewport or FoV (visible area) of the user and sending only the selected portion of the panoramic frame to the client. Despite the great promise and research, viewport-adaptive streaming has not yet been adopted by streaming service providers like YouTube and Facebook. This is primarily due to lack of end-device software support in viewport-adaptive streaming and the lack of holistic understanding of how user viewports change from frame to frame and from video to video.

There have been scattered efforts to understand viewport characteristics of \ang{360} videos, but with small-scale experiments with limited number of videos~\cite{xie_cls_2018,rossi_spherical_2019,petrangeli_trajectory-based_2018}. This paper presents a comprehensive analysis of viewport trajectories of \ang{360} video users by analysing a large unified dataset of over 3700 viewport traces. \emph{To the best of our knowledge, this is the largest \ang{360} video user behavior analysis conducted to date.} We propose a viewport clustering mechanism to understand similarities of user behaviors that takes spatial, motion and other behavioral patterns into account. We particularly highlight the importance of adding temporal features such as head movement speed to the viewport clustering to ensure that subtle differences in highly dynamic viewport patterns can be captured.

Video categorizations are widely studied for normal video streaming~\cite{sun2016exploiting,jiang2017exploiting}, which
can support video indexing, storage management, popularity prediction, and efficient resource provisioning at both content and network service providers~\cite{sun2016exploiting, mcclanahan2017interplay}.
However, these traditional video categorizations will not be effective anymore for \ang{360} videos as they categorize entire video into one category or cluster~\cite{afzal_characterization_2017, carlsson_had_2020, nasrabadi_taxonomy_2019}. Due to the freedom of users to request different portions of the same frame and tile-based streaming of \ang{360} videos, user behavior can vary significantly even within a given video. In this paper, we propose a novel video segment categorization mechanism based on the dynamic viewport characteristics of users that cluster chunks of the same video into different groups.

The key contributions of this paper are;\vspace{-1mm}
\begin{itemize}
    \item An aggregated dataset with 88 videos with over 3700 viewport traces containing over 142 hours of viewing logs. We consolidate different formats of individual datasets, e.g. head orientation, data sampling rates, trace duration, into unified framework, which has been made publicly available\footnote{\url{https://github.com/theamaya/Viewport-Aware-Dynamic-360-Video-Segment-Categorization}}.
    \item A novel feature based viewport clustering algorithm that takes spatial, motional and other behavioral features of viewport patterns into account. We validate the effectiveness of the proposed clustering mechanism using the aggregated dataset compared to previously proposed clustering methods. The results show an average of 81.17\% of viewport overlap and limits the pairwise difference in angular speed of 0.47 rad/s for the traces within the resulting clusters.
    \item A novel dynamic video segment categorization algorithm which is based on different user behavior types induced by the corresponding video. The results show 28.32\% improvement in similarity in viewport distribution within the clusters when compared to the existing static genre-based video categorizations.
\end{itemize}

\section{Related Work}

\textbf{Viewport clustering}.
Majority of the existing algorithms for viewport clustering, are based on the hypothesis that users tend to exhibit similar viewing patterns for a given video. One notable clustering approach, proposed by \cite{rossi_spherical_2019} leverages the spatial arrangement of user viewpoints at a particular play-point or during a specific time window of the video. The clustering is done in the spherical domain using geodesic distance as the closeness metric. Another approach proposed by \cite{xie_cls_2018} extracts user fixations from viewport traces and performs density-based clustering using a spatial distance metric. The algorithm proposed by \cite{petrangeli_trajectory-based_2018} utilizes viewport trajectory data in the long term and clusters the viewports with a time and distance based metric.

\textbf{Video categorization}.
The literature surrounding video categorization is limited, with most studies focusing on video titles and genres. Afzal \textit{et al.}~\cite{afzal_characterization_2017} categorized 2285 \ang{360} videos based on their genre which is determined by keywords and found 14 different categories, ranging from documentaries to sports and driving. Another categorization approach, proposed by \cite{almquist_prefetch_2018} and \cite{carlsson_had_2020} used manual inspection of video content and expected behaviors of user FoV trajectories, resulting in categories such as \textit{Exploration}, \textit{Static focus}, \textit{Rides} etc. An approach proposed by \cite{nasrabadi_taxonomy_2019} categorizes videos based on the extent of camera motion and the number of moving targets.

\textit{In contrast to these methods, we aim to leverage a viewport-aware approach to the problem of \ang{360} video categorization. Our proposed solution aims to go beyond spatial features by considering temporal and behavioral patterns as metrics for viewport similarity modeling.}

\section{Dataset}\label{sec:dataset}
We created a unified dataset  combining 6 different \ang{360} video datasets from the literature. The selected datasets contain \ang{360} videos
and their corresponding user head orientation data. Table~\ref{tab:dataset} summarises the
datasets used. The aggregated dataset contains 88 videos with an average of 45 user traces per video. Each video is at least 30s long. We have made the aggregated dataset publicly available \href{https://github.com/theamaya/Viewport-Aware-Dynamic-360-Video-Segment-Categorization/}{\underline{here}}.

First,
we converted all head orientation traces to the same representation; Yaw and Pitch angles in radians (Roll angle variation was negligible, therefore was ignored in this study). Next,
to tackle the differences in sampling rates, we re-sampled the entire dataset at 10Hz.
For this study, we only considered the datasets collected using
HMDs, excluding
mobile and PC users to eliminate the impact of device on viewing patterns~\cite{xu_analyzing_2019}\footnote{We keep viewport patterns from different rendering devices and its impact on our categorization task as a future work}.

\vspace{-2mm}
\begin{figure}[h]
\centering
    \subfloat[Viewport fixations for 4 sample videos]{\includegraphics[trim=0 0 0 0, clip,width=0.7\columnwidth]{Motivation_FOV_added.pdf}%
    \label{fig:motivation_2_diff_videos}}

    \subfloat[Mean distribution: Yaw]{\includegraphics[width=0.35\columnwidth]{Motivation_feature_3.png}%
    \label{fig:motivation_f3_dist}}\quad

    \hspace{0.5mm}
    \subfloat[Max yaw movement]{\includegraphics[width=0.35\columnwidth]{Motivation_feature_11.png}
    \label{fig:motivation_f11_dist}}

    \vspace{-3mm}
    \caption{Viewport and two feature distribution between selected videos in two categories}
    \vspace{-4mm}
    \label{fig:motivation}
\end{figure}

\vspace{-2mm}
\subsection{Relationship of viewport and video types}

In our preliminary analysis, we studied how viewport distributions differ between different videos. Fig.~\ref{fig:motivation} shows four example videos analyzed in this study. In both \textit{i)} V1 (Rollercoster) \& V2 (Riding on a Hog): the camera is moving forward with a high speed whereas in \textit{ii)} V3 (Landscape) \& V4 (Panel): Camera movement is not significant, and contents are dispersed in the frame.

Viewport distributions in Fig.~\ref{fig:motivation_2_diff_videos} shows that, content of the video has an impact on the users' viewport distribution. We see that the viewports of V1 \& V2 are centered in the frame whereas, viewports are comparably dispersed in V3 \& V4. This indicates that,
viewport distributions can be used as a proxy to explain certain content patterns in \ang{360} videos
indicating that a viewport based categorization of videos can be effective. Fig.~\ref{fig:motivation_f3_dist} plots the spatial distribution of 2s viewport chunks in the yaw axis, and~\ref{fig:motivation_f11_dist} plots the maximum yaw angle moved by a viewport trace within a 2s chunk. We observe that patterns of raw viewport distributions in \ref{fig:motivation_2_diff_videos} are reflected in these features, indicating that carefully extracted features from viewport traces can result in an effective video categorization.

\section{Conclusion and Future work}

In this study, we conducted an in-depth analysis of user behaviors in \ang{360} video streaming using a large dataset of over 3700 viewport traces.
By extracting both spatial and temporal features from viewport traces, the proposed method can identify viewport similarities across different videos whilst achieving
81.17\% \textit{within-cluster} viewport overlap and outperforms the existing viewport clustering algorithms.
We extended the above solution to develop a novel \ang{360} video categorization which successfully captures temporal changes of content and user behaviors by independently processing 2s video chunks. The proposed categorization shows 28.32\% improvement in clustering videos with similar combinations of user behaviors over existing genre-based \ang{360} video categorizations.

In future work, we aim to extend our video categorization method for to develop a novel caching algorithm to enhance the \ang{360} video processing at edge servers. We plan to explore new approaches for edge resource optimization further analysing the content of the categorized videos.

\section{feature-based viewport clustering}\label{section 4}

In this section, we propose a feature representation for user viewport traces, which is used for a novel viewport clustering method. In Section \ref{section 5}, we extend this feature based viewport clustering to develop a novel dynamic video categorization.
We divided the duration of a video/ viewport trace into bins (chunks) of 2 seconds.
Table~\ref{tab:basic notation} briefly describes the notations used in our further explanations.

\begin{table}[h]
  \small
  \vspace{-2mm}
  \caption{Feature selection for viewport clustering}
  \vspace{-4mm}
  \begin{tabular}{l|l}
    \toprule
    Notation & Description\\
    \midrule
    $i$, $k$ &  video id: $i \in [ 0,87]$, chunk id: $k \in [ 0,13]$,  \\
    $n_i$ & num. of users for $i^{th}$ video\\
    $j$ & user id, $j\in [ 0,n_{j})$\\
    $\{V_{i}\}$ & Set of all videos $\forall i \in [ 0,87]$ \\
    $\{c_{k}\}$ & Set of all time chunks $\forall k \in [ 0,13]$ \\
    $\{u_{i,j}\}$ & Set of all users $\forall i \in [ 0,87], \forall j\in [ 0,n_{i})$ \\
    $\{v_{i,k}\}$ & Set of all video chunks $\forall i \in [ 0,87], \forall k \in [ 0,13]$\\
    $\{t_{i,k,j}\}$ & Set of all viewport trace chunks  $\forall i \in [ 0,87], \forall k \in [ 0,13],$\\
    & $\forall j\in [ 0,n_{i})$\\
    \bottomrule
\end{tabular}
\label{tab:basic notation}
\vspace{-2mm}
\end{table}

\vspace{-2mm}
\subsection{Feature extraction}\label{Section 4- feature extraction}

We observed that the
\textit{users who consume \ang{360} video content differ to each other in multiple aspects such as positioning, speed of head movement, extent of the sphere explored and the maximum angle moved within a defined window}.

We extracted features for every viewport trace chunk in $\{t_{i,k,j}\}$ as summarized in Table \ref{tab:Feature_selection_for_viewport_clustering}. Let $F^{( 1)}$ represent the feature extraction function for a viewport trace chunk. $\underline{f}^{( 1)}_{i,k,j}$ is a 15 dimensional feature vector that represent a viewport trace chunk in the dataset.

\vspace{-2mm}
\begin{equation}\label{eq:feature_extraction_func_vp_clust}
 F^{( 1)}( t_{i,k,j}) =\underline{f}^{( 1)}_{i,k,j}
\end{equation}
\vspace{-3mm}

\begin{table}[h]
  \small
  \caption{Feature selection for viewport clustering}
  \vspace{-4mm}
  \begin{tabular}{p{2.0cm}| p{5.5cm}}
    \toprule
    Feature & Measurements y--yaw, p--pitch\\
    \midrule
    Position & mean, 25\textsuperscript{th} and 75\textsuperscript{th} percentiles: y \& p \\
    Speed &  mean, 25\textsuperscript{th} and 75\textsuperscript{th} percentiles: y \& p \\
    Maximum angle &

    angular position relative to initial point: y \& p\\
    Area explored & \% coverage on the sphere by \ang{100}x\ang{100}  viewport\\
  \bottomrule
\end{tabular}
\label{tab:Feature_selection_for_viewport_clustering}
\end{table}

\vspace{-2mm}
\subsection{Clustering the viewport trace chunks}

Fig.~\ref{fig:view_port_based_cat} summarizes the complete procedure of generating feature vector representations for viewport traces
which are then used as input to  a K-Means algorithm that clusters the viewport trace chunks into \textit{M} clusters.
We extract the cluster to which a viewport trace chunk belongs to using Eq.~\ref{eq:cluster_alloction_func} (further used in Section \ref{section 5})

\begin{equation}\label{eq:cluster_alloction_func}
    C^{( 1)}( t_{i,k,j}) =m,\ m\in [ 0,M)
\end{equation}

We defined three pairwise metrics to quantify the similarity of two viewport trace chunks for this viewport cluster analysis.
\begin{itemize}
\item\textit{Pairwise viewport overlap- VPO-} \cite{rossi_spherical_2019} measures the spatial coherence of two traces. It calculates the geodesic distance between each two corresponding sample points of 2 viewport traces and uses it as a proxy for pairwise viewport overlap. In our representation, a viewport trace chunk has 20 points (10Hz, 2 seconds).
The average of geodesic distance between each pair of corresponding points is taken as pairwise VPO.

\item\textit{Pairwise difference in head movement speed-}The average angular speed is calculated for each viewport trace segment, and their pairwise absolute difference is measured.

\item\textit{Pairwise difference in percentage of sphere explored-}within the 2s duration from the features extracted as mentioned in \ref{tab:Feature_selection_for_viewport_clustering}.

\end{itemize}
\textit{maximization} of the Pairwise VPO, and \textit{minimization} of Pairwise difference in head movement speed and Pairwise difference in percentage of sphere explored are desirable within a cluster.

\vspace{-1mm}
\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{Blcok_Diagram_VP_2.pdf}%
\label{fig:viewport_clustering_based_categorization}
\vspace{-4mm}
\caption{Overview of viewport clustering}
\label{fig:view_port_based_cat}
\vspace{-4mm}
\end{figure}

\subsubsection{The effectiveness of the proposed feature representation}

Clustering similar user viewport traces is often used in viewport prediction for \ang{360} video streaming. In this experiment we prove that the feature representation introduced in Eq.~\ref{eq:feature_extraction_func_vp_clust} can be used to accurately cluster viewports in such an application.
The feature vectors $\underline{f}^{( 1)}_{i,k,j},\ \forall j\in [0,n_{i}]$ as input to the K-Means clustering to yield \textit{M} different clusters separately for different videos $i$ and different time bins $k$.
The results were compared with 3 other algorithms,

\textbf{Spherical Clustering\cite{rossi_spherical_2019}}- Models viewport traces as nodes of a unit distance graph. Edges are connected if the two viewport traces have high viewport overlap. Finds maximal cliques as clusters

\textbf{Trajectory based\cite{petrangeli_trajectory-based_2018}}- Represents viewport traces as graph nodes with edges weighted by their spatial similarity. Uses spectral clustering while minimizing the distortion score.

\textbf{DBSCAN\cite{xie_cls_2018}}- Maximizes spatial density of clusters by taking the mean positions of trajectories in yaw, pitch in the 2s interval

\begin{figure}[h]
\vspace{-2mm}
\centering
\subfloat[Pairwise viewport overlap]{\includegraphics[width=0.18\textwidth]{BOX_PLOT_VPalgoComparison-VPoverlap-_2_.pdf}%
    \label{fig:Pairwise viewport overlap}}

    \hspace{2mm}
\subfloat[Pairwise difference in head movement speed]{\includegraphics[width=0.18\textwidth]{BOX_PLOT_VPalgoComparison-Speeddif-_2_.pdf}%
    \label{fig:Pairwise difference in head movement speed}}
\vspace{-3mm}
\subfloat[Pairwise difference in percentage of sphere explored ]{\includegraphics[width=0.18\textwidth]{BOX_PLOT_VPalgoComparison-_sphere.pdf}%
    \label{fig:Pairwise difference in percentage of sphere explored }}

    \hspace{2mm}
\subfloat[Number of resulting clusters ]{\includegraphics[width=0.115\textwidth,angle=90 ]{BOX_PLOT_VPalgoComparison-nclusters-_2_.pdf}%
    \label{fig:Number of resulting clusters }}
\vspace{-4mm}
\caption{Within-cluster measurement comparison of T-Trajectory based\cite{petrangeli_trajectory-based_2018}, S- Spherical\cite{rossi_spherical_2019}, D- DBSCAN\cite{xie_cls_2018} with the P- proposed for the aggregated dataset.}
\vspace{-2mm}
\label{fig:VP clustering algorithm comparison}
\end{figure}

Fig.~\ref{fig:VP clustering algorithm comparison} contains the aggregated comparison of the proposed algorithm with \cite{rossi_spherical_2019, petrangeli_trajectory-based_2018, xie_cls_2018} for the task of viewport clustering. The analysis was done for all $\{v_{i,k}\}, \forall i \in [ 0,87], \forall k \in [ 0,13]$.
For a clear comparison, at each video and time chunk we used the resulting number of clusters from Spherical clustering\cite{rossi_spherical_2019} as the number of clusters for K-means in the proposed algorithm. It is clear that the proposed algorithm outperforms all others in speed and percentage exploration based metrics, and shows on par performance with Spherical clustering algorithm in \textit{VPO} metric.

Fig.~\ref{fig:Temporal VP clustering vs proposed} contains a temporal comparison of the 4 viewport clustering algorithms for two videos picked from the dataset; Rollercoaster and Timelapse.
Majority of \textit{Roller-coaster} users fixate on the center of the frame while \textit{Time-lapse} users explore the entire frame freely without a significant correlation. The number of clusters generated by each algorithm at each time point is indicated on the corresponding plot on Fig.~\ref{fig:VP overlap rollercoaster}~and~\ref{fig:VP overlap timelapse}. Since Spherical clustering\cite{rossi_spherical_2019} shows the best performance out of the three algorithms used for comparison, we used the same number of clusters yielded by Spherical clustering for our algorithm at each time bin.

\begin{figure}[h]
\centering
\subfloat[Viewport overlap - rollercoaster]{\includegraphics[width=0.244\textwidth]{VPAlgoCompareTemporal_VPOverlapRollercoaster.pdf}%
    \label{fig:VP overlap rollercoaster}}

\subfloat[Viewport overlap - timelapse]{\includegraphics[width=0.244\textwidth]{VPAlgoCompareTemporal_VPOverlapTimelapse-_2_.pdf}%
    \label{fig:VP overlap timelapse}}
\hfill
\subfloat[Pairwise Speed diff. -rollercoaster ]{\includegraphics[width=0.244\textwidth]{VPAlgoCompareTemporal_SpeeddiffRollercoaster.pdf}%
    \label{fig:Speed difference rollercoaster}}

\subfloat[Pairwise Speed diff. -timelapse ]{\includegraphics[width=0.244\textwidth]{VPAlgoCompareTemporal_SpeeddiffTimelapse.pdf}%
    \label{fig:Speed difference timelapse}}
\hfill
\subfloat[Pairwise \% sphere diff. -rollercoaster ]{\includegraphics[width=0.244\textwidth]{VPAlgoCompareTemporal__sphereRollercoaster.pdf}%
    \label{fig:percent sphere rollercoaster}}

\subfloat[Pairwise \% sphere diff. -timelapse ]{\includegraphics[width=0.244\textwidth]{VPAlgoCompareTemporal__sphereTimelapse.pdf}%
    \label{fig:percent sphere timelapse}}
\vspace{-4mm}
\caption{Comparison of the viewport (VP) clustering algorithms for Rollercoaster and Timelapse videos. }
\vspace{-2mm}
\label{fig:Temporal VP clustering vs proposed}
\end{figure}

Spherical clustering algorithm produces significantly high number of clusters as it is based on maximal clique clustering.
Nevertheless, for the same number of clusters, the proposed algorithm yields almost on par performance in terms of viewport overlap.
Furthermore, our proposed algorithm significantly outperforms all others in pairwise speed and pairwise \% sphere explored metrics by explicitly using those features to define the viewprot trace chunk. The literature suggests that having viewport traces with similar speeds, orientation and spatial distribution within a cluster facilitates more accurate and long term viewport prediction. \textit{Viewport prediction followed by clustering is not in the scope of this paper.}

\textit{\textbf{Takeaway-} Unlike other algorithms that rely primarily on the spatial details of the viewports, the proposed feature extraction and clustering method successfully identifies clusters of viewport trace chunks that not only attend to the same content, but are also similar in behavioral aspects such as the speed of head of movement, maximum angular displacement within the window, and the extent of the sphere viewed. Moreover, the proposed algorithm is fast and lightweight in comparison to Trajectory based\cite{petrangeli_trajectory-based_2018} and Spherical clustering\cite{rossi_spherical_2019} approaches that utilize algorithms with a high computational overhead.}

\vspace{-2mm}
\subsection{Viewport clustering for aggregated dataset}\label{subsection:Viewport clustering for the entire dataset}

As the first step of our dynamic video categorization, we generated clusters inside the set of all viewport trace chunks $\{t_{i,k,j}\}$ as introduced in Table\ref{tab:basic notation}. The same feature vectors derived in \ref{Section 4- feature extraction}, $\underline{f}^{( 1)}_{i,k,j}$ were used as inputs to the K-Means clustering algorithm.

The optimum number of clusters for the set $\{t_{i,k,j}\}$ was determined by a DB (Davies- Bouldin) Score analysis as indicated in Fig.~\ref{fig:DBScore analysis VP}. The graph gives 2 minimums; \textit{M}=3 and \textit{M}=10. According to the plot, \textit{M}=3 is not a stable minimum and as shown in Fig.~\ref{fig:Similarity measurement within clusters for M=3}~and~\ref{fig:Similarity measurement within clusters for M=10}, upon clustering the dataset to 3 clusters, 2.29\% of the datapoints were identified as outliers. Where as for 10 clusters the outlier percentage was only 1.12\%. \textbf{Therefore, we selected \textit{M}=10}.

\vspace{-2mm}
\begin{figure}[h]
\centering
\subfloat[.]{\includegraphics[width=0.235\textwidth]{DBScoreVPClusters_reduced_size.pdf}
    \label{fig:DBScore analysis VP}}

\subfloat[.]{\includegraphics[width=0.245\textwidth]{BehavioursPerVideoChunk_reduced_size.pdf}
    \label{fig:Behaviours per video chunk}}\\

\vspace{-4mm}
\caption{a)-DB-score analysis to determine the optimum number of clusters for viewport trace chunk clustering. b)-Number of different user behaviors available in video chunks and the CDF of their distribution}\vspace{-5mm}
\label{fig:Similarity measurement VP m=10}
\end{figure}

\begin{figure}[h]
\centering
\subfloat[M=3 ]{\includegraphics[width=0.16\textwidth, , angle=90]{3VPClustersWithinBOXPLOT2.pdf}
    \label{fig:Similarity measurement within clusters for M=3}}

\subfloat[M=10 ]{\includegraphics[width=0.34\columnwidth, , angle=90]{10VPClustersWithinBOXPLOT2.pdf}
    \label{fig:Similarity measurement within clusters for M=10}}\\

\vspace{-4mm}
\caption{Pairwise VPO and Pairwise difference in head movement speed within and between clusters for \textit{M=3, 10}}
\vspace{-4mm}
\label{fig:Similarity measurement VP m=10}
\end{figure}

Fig.~\ref{fig:Similarity measurement VP m=10} illustrates the results for similarity measurement within and between viewport clusters generated by the proposed algorithm.
We observe that for each cluster, the resulting percentage pairwise viewport overlap between the traces of a cluster is significantly higher than that between traces from different clusters.
The averaged pairwise viewport overlap within viewport clusters is 81.17\%, and that between different clusters is 47.80\%.
Moreover, the pairwise difference of viewport trace speeds within a cluster is significantly less than that between different clusters. While the averaged pairwise speed difference within clusters is 0.47 rad/s, that between different clusters is as high as 2.47 rad/s
It is evident that the proposed viewport clustering algorithm captures similarities in terms of both spatial and motion features.

\vspace{-2mm}
\begin{figure}[h!]
\centering
\subfloat[Viewport trace cluster 7]{\includegraphics[width=0.22\textwidth]{VPClusterHeatmap7.pdf}
    \label{fig:Viewport trace cluster 7}}

\subfloat[Viewport trace cluster 9]{\includegraphics[width=0.22\textwidth]{VPClusterHeatmap9.pdf}
    \label{fig:Viewport trace cluster 9}}
    \hfill

\subfloat[Cluster 7 features]{\includegraphics[width=0.18\textwidth]{VP_clusters_combined_reduced_fig_size_fig_1.pdf}
    \label{fig:features cluster 7}}

\subfloat[Cluster 9 features]{\includegraphics[width=0.26\textwidth]{VP_clusters_combined_reduced_fig_size_fig_2.pdf}
    \label{fig:features cluster 9}}
    \hfill
\vspace{-3mm}
\caption{Viewport distributions and corresponding features of clusters 7 and 9}
\vspace{-3mm}

\label{fig:Viewport distributions of different clusters}
\end{figure}

Fig.~\ref{fig:Viewport trace cluster 7} and \ref{fig:Viewport trace cluster 9} show example intermediate results of viewport clustering indicating viewport trace distribution on Equirectangular frame for clusters 7 and 9. Each colored dot indicates one sample of viewport centers of a given user in the cluster. A viewport trace for 2s is composed of 20 such points.
It is evident how different user behaviors are captured by different clusters. Cluster 7 contains user traces that are slowly moving or are stationary around the center of the frame while cluster 9 contains users that are moving fast in the equatorial axis biased to the right of the frame.
Each cluster can contain traces with different combinations of features describing position and motion as shown in in Fig.~\ref{fig:Viewport distributions of different clusters}.
Fig.~\ref{fig:Behaviours per video chunk} shows the number of different viewport trace clusters (different behaviors) to which the users of a particular video chunk belongs. To eliminate the effects of outliers, a behavior is counted only if there are two or more users belonging to that particular cluster. According to the CCDF, More than 75\% of all video chunks have users displaying more than 4 different user behavior types (out of 10).

\textit{\textbf{Takeaway-} Different users consuming different \ang{360} videos show similarities in terms of positioning, speed, and \% of sphere explored by their viewports. The novel feature definition successfully identifies such similar behaviors. Majority of the video chunks have users displaying a variety of user trace patterns. Hence, for a successful characterization of videos, all different user behaviors should be used }

\section{Dynamic Video chunk categorization} \label{section 5}

In this section, we present a novel algorithm to dynamically categorize video chunks. Inspired by the results of Section \ref{subsection:Viewport clustering for the entire dataset}, we statistically characterise a video chunk by the distribution of how the users perceive its content by using different viewport trace patterns. We start off with the hypothesis that \textit{for a large enough dataset, the aggregate set of all user viewport trace chunks can reasonably represent all practical types of short-term user behaviors in perceiving 360°videos.}
In the following subsections we define a feature representation for video chunks based on our results from section \ref{subsection:Viewport clustering for the entire dataset}, categorize the video chunks into optimum number of clusters, and evaluate the performance of the proposed categorization.

\vspace{-2mm}
\subsection{Feature engineering for video chunks}
Let ${F^{( 2)}( v_{i,k}) =\underline{f}^{( 2)}_{i,k}}$ be the feature extraction function for video chunks. .
For a particular video chunk $v_{i,k}$ it generates an \textit{M} dimensional feature vector $\underline{f}^{( 2)}_{i,k}$ as follows. (\textit{M} is the number of viewport trace clusters as explained in section \ref{subsection:Viewport clustering for the entire dataset}. We selected \textit{M}=10.)

\vspace{-2mm}
\begin{equation}\label{eq:vp_traces_vector_4_vid_clustering}
{\underline{f}^{( 2)}_{i,k}} =<\frac{|A^{0}_{i,k} |}{n_{i}},\frac{\ |A^{1}_{i,k} |}{n_{i}},\ ...,\frac{\ |A^{M-1}_{i,k} |}{n_{i}}>
\end{equation}
\vspace{-2mm}

such that $A^{m}_{i,k} =\left\{t_{i,k,j} \ |\ C^{( 1)}( t_{i,k,j}) =m,\ \forall j\in [ 0,\ n_{i})\right\}$

In a gist, the feature representation for a particular video chunk indicates the population distribution of its users that belongs to each of the  different viewport trace clusters (representing different user trace behaviors) resulting in section \ref{subsection:Viewport clustering for the entire dataset}. For this analysis we used 1232 unique \ang{360} video chunks. (88 videos, 14 time bins)

\vspace{-2mm}
\subsection{Video chunk categorization}\label{Video chunk categorization}

Next we generate clusters of similar video chunks in the defined \textit{M} dimensional feature space as shown in Fig.~\ref{fig:overview_dynamic_video_cat}. $\underline{f}^{( 2)}_{i,k} \ \forall i\in [ 0,87] ,\ \forall k\in [ 0,14)$ are taken as inputs to a K-Means clustering algorithm that yields \textit{Q} clusters. Similar to Eq.~\ref{eq:cluster_alloction_func}, corresponding cluster/category for a video chunk is given by $C^{( 2)}( v_{i,k}) =q,\ q\in [ 0,Q)$.

We define a category as a collection of video chunks that are similarly perceived by the users. i.e, videos with similar distributions of different user trace behaviors. The following distance metric for 2 video chunks $v_{i,k}$ and $v_{i',k'}$ quantifies the above definition.

\vspace{-2mm}
\begin{equation}\label{eq:similarity_metric2}
    S( v_{i,k} ,v_{i',k'}) =\sqrt{\sum\limits ^{M}_{m=0}\left( |A^{m}_{i,k} |-|A^{m}_{i',k'} |\right)^{2}}
\end{equation}
\vspace{-2mm}

It calculates the euclidean distance between the feature vectors of two video chunks, quantifying the difference in how the users of two videos are distributed among the 10 identified user trace patterns/behaviors.

We conducted a Davies-Bouldin score analysis to determine the optimum number of categories as shown in Fig.~\ref{fig:DB-Score_analysis_for_optimum_Q}, which yielded \textit{Q}=6.
Fig.~\ref{fig:overview_dynamic_video_cat} summarizes the feature extraction and categorization of video chunks.
Fig.~\ref{fig:Similarity measurement within and between video categories New Metric} represent the \textit{within-category} and \textit{cross-category} evaluation of the six dynamic video chunk categories. On average, the within-cluster pairwise feature distance is 0.27 and the cross-cluster pairwise feature distance is 0.56.

\vspace{-2mm}
\begin{figure}[h!]
\centering
\subfloat[DB-Score analysis for optimum \textit{Q} ]{\includegraphics[width=0.24\textwidth]{DBScoreVideoClusters_reduced_size.pdf}
    \label{fig:DB-Score_analysis_for_optimum_Q}}

\subfloat[Pairwise similarity]{\includegraphics[width=0.24\textwidth]{6VideoCategoryWithinBOXPLOTNewMetric_reduced_size.pdf}
    \label{fig:Similarity measurement within and between video categories New Metric}}
\vspace{-2mm}
\caption{(a)- DB-score analysis to determine the optimum number of categories, (b)- \textit{within-category} and \textit{cross-category} analysis of pairwise distance measurement calculated according to Eq.~\ref{eq:similarity_metric2} }\vspace{-4mm}
\label{fig:Optimal number of video categories}
\end{figure}

Fig.~\ref{fig:Video categories heatmaps and feature hist} shows heatmaps of equirectangular spatial viewport distribution for all six dynamic video categories obtained by combining the viewport traces of all of the users corresponding to its video chunks. A clear difference is observed between different spatial viewport distributions. In Fig.~\ref{fig:Video categories heatmaps and feature hist}, a histogram representation of the feature vectors of video chunks in each category is indicated below each corresponding heatmap. It is evident that the proposed dynamic video categorization algorithm successfully identifies different combinations of user trace types for video chunk clustering.

\textit{\textbf{Takeaway- }In viewing any \ang{360} video, there are different users who display different behaviors. It is advantageous to take all of these different behaviors into account, rather than generalizing it to the majority behavior when characterizing a video. The proposed video chunk categorization clusters the videos that induce similar types of behaviors in their users, irrespective of the video content.}

\vspace{-1mm}
\subsection{Static vs Dynamic video chunk clustering}

Majority of the existing categorizations of \ang{360} videos are static; the entire video belongs to a particular category. These static methods often fail since significant temporal variations in video content and corresponding user behavior can occur in long videos.
To the best of our knowledge, ours is the first proposed algorithm that dynamically categorizes \ang{360} videos by taking 2s video chunks. Fig \ref{fig:DynamicClustersPerVideo} shows how different chunks of the same video can belong to different categories according to the proposed algorithm.

We compared the performance of the proposed categorization with genre-based static categorization by \cite{afzal_characterization_2017}. The 88 videos in the dataset were assigned to 10 genres such as \textit{Sports, scenety etc.}
Fig.~\ref{fig:static vs dynamic new metric} shows the comparison results using \textit{within-category} and \textit{cross-category} feature distances as introduced in Eq.~\ref{eq:similarity_metric2}.

\vspace{-2mm}
\begin{figure}[h!]
\centering
\subfloat[\# of dynamic clusters per video  ]{\includegraphics[width=0.22\textwidth]{nDynamicClustersPerVideo.pdf}
    \label{fig:DynamicClustersPerVideo}}

    \hspace{0.5mm}
\subfloat[Static vs dynamic cat. ]{\includegraphics[width=0.22\textwidth]{GenreVsDynamicNewMetric.pdf}
    \label{fig:static vs dynamic new metric}}
    \vspace{-2mm}
\caption{a)- \# of dynamic clusters per a particular video b)- Comparison static vs dynamic categorizations using Eq \ref{eq:similarity_metric2}}\vspace{-2mm}
\label{fig:dynamic cat plots}
\end{figure}

\textit{\textbf{Takeaway- }Static categorizations of a long \ang{360} videos fail to take temporal variations of contents and user behaviors into account. Breaking the videos down to independant chunks enables more accurate categorizations. Quantifying the categorization in terms of user viewport patterns eliminates ambiguities and the requirement of semantic information for categorization.}


\section{Introduction}

Most of the modern astronomical sky surveys are characterized by fast pace data ingestion, data intensive science cases or automatic reduction pipelines \citep[e.g.][]{feigelson2012big}, which often lie on the verge of technological developments and analysis capabilities.
This unprecedented availability of observations challenges the traditional approaches for data analysis, leading to a shift in the paradigm for knowledge discovery \citep{bell2009beyond}, notably dominated by machine learning (ML) techniques \citep{ball2010data}.
Albeit the difficult mathematical and statistical foundations, a complex terminology driven by the confluence of several sciences and the arduous interpretation of the results, the training of intelligent agents has become an every-day practice in astronomy. The accessibility of easy-to-use free software resources mostly written for \textit{R} \citep{team2000r} or \textit{Python} \citep{van2003python} languages was fundamental to this step.

In most cases ML methods can be separated into two basic steps: first, raw data is converted into a set of useful features that are relevant to the task at hand (e.g. periods or intensities) and then these features are fed to a classifier or a statistical method \citep[see e.g. ][]{mitchell1997machine}.

ML methods have a number of limitations. For instance, they are highly susceptible to errors produced by the limitations of the datasets \citep{cai2015challenges}.
The results can also be hampered by the not fully understood role of the features \citep{duboue2020art} or by  the biases introduced by improperly defined experiments \citep{domingos2012few}.
These facts are well known and have not been ignored in the astronomical community \citep{2020MNRAS.492.5377L}.

Here we are interested in the role of some sources of noise that are present in commonly used features in astronomical research, and their impact on the results of ML methods in this context. We use data from the synoptic survey ``\textit{Vista Variables in Via L\'{a}ctea}'' (VVV, \citealp{2010NewA...15..433M}), observed with the Vista telescope \citep{sutherland2015visible}, which pursues, among its main objectives, to produce a three-dimensional map of a large part of the galactic center (Bulge) of the Milky Way and a fraction of the internal Galactic Disk. The VVV data is presented in units called ``tiles'', which are rectangular areas of the sky surveyed over time. For each tile, the VVV data reduction pipeline \citep{emerson_vista_2004} provides a pre-processed image and a database of files with the values of position, magnitude and color indices of the light sources present in the image, which comprises the ``photometric catalog''. These catalogs are the main subject of this study.

The images are subject to two noise sources, namely, experimental errors and observation conditions. The derived catalogs are also affected, since the noise permeates all the survey information which is comprised on a set of features or observables. Atmospheric conditions, moon phases, maintenance of the camera and telescope or modifications to the software, among many others, can influence the observation or recording of the data. As a consequence, the derived measurements that are used as features for a ML analysis can be also prone, at different levels, to these errors and conditions.

Random measurement errors are present in every experimental or observational science
They are unavoidable, but each error typically affects only a single or a reduced set of observations\rev{, in particular wide-field surveys images can be affected by issues such as weather, astronomical conditions, software updates, etc.} ML methods can cope efficiently with this kind of errors. For a large survey as VVV, observational conditions can change slightly (but not randomly) over long periods of time or for different regions in the sky.
This problem is more difficult for ML methods.
In many situations we want to train an intelligent agent using a well--known portion of the survey, and then use it to predict other less--known zones searching for a given astronomical phenomenon.
Given the ML methodology, the agent will work efficiently on training data, but will probably fail to generalize to other zones due to this slight change in observational conditions.
Due to the diverse nature of the features extracted from the data (intensity, periods, colours, etc.), possibly they will reflect this effect in different proportions.
It is then interesting to ask whether it is possible to automatically detect which of the extracted features are more sensitive to these changes in observational conditions.
Hereafter we call ``Drifting Features'' to the ones in a dataset that are sensitive to the observational conditions. We aim at evaluating their influence over a large scale ML experiment.

As a working example, we focus on the problem of detecting RR-Lyrae (RRLs) variable objects over VVV data. That is, we train classifiers using data from some VVV tiles, and evaluate them in the task of identifying RRLs on other tiles.

Drifting Features should be consistent within a limited zone of the sky (for instance, a tile or two consecutive tiles) but should show slight changes, almost undetectable by most simple statistics, between tiles that are separated apart\footnote{\rev{Whether two consecutive tiles or, in general, two regions on the sky are similar depends on the survey strategy. For example, regions that are repeatedly observed in quick succession in a survey, or that share important observational parameters such as the exposure time, could show this property even if they have long angular distances. The opposite is also valid, regions that are close at angular distances can be observed under very diverse conditions.}}. Those changes could potentially alter the capabilities of the classifier. To detect these features and their effect on automatic classification, we propose to use, again, ML methods. If we face a ML method with the task of discriminating data from two tiles, it will be forced to learn differences among the tiles that are present in the features. We can then use Feature Selection methods \citep{guyon2002gene} to evaluate the importance of each feature for this classifier that discriminates tiles, and mark highly relevant features as candidates of been Drifting Features. \rev{In other words, we propose to learn a separate task (the tile of origin of a source) not because it is unknown or difficult, but as a method to detect which are the features that are more useful to this task, the features that contain information that change with the tile of origin.}

This work is divided into the following sections: in Section~\ref{section:vvv}
we explain our experimental setup (data, feature extraction, model selection,etc.). In Section \ref{section:drift} we introduce our procedure for the identification of Drifting Features and in Section \ref{section:evaluation} we evaluate the effect of these features in the task of RRL identification. Finally in Section \ref{section:conclusions} we discuss our results and draw our conclusions.

\section{Experimental setup}
\label{section:vvv}

\subsection{Data}

One of the main objectives of the VVV is the creation of a 3D map of the bulge and the galactic center \citep{2010NewA...15..433M} for which the search for variable stars in general, and RRLs in particular, is important due to their use as standard candles \citep{bailey1902discussion}.
To this end, the survey relies on data from the VIRCAM infrared camera, mounted on ESO's VISTA survey telescope \citep{sutherland2015visible}, which at the time of its construction was the largest NIR camera with 16 non-contiguous \mbox{2k x 2k} detectors.
To complete a contiguous tile, VIRCAM simultaneously exposes its detectors 6 times with a suitable offset. Each one of the exposures is called a \textit{pawprint} and the combination of the six overlapping pawprints is a \textit{tile}. For this reason each pixel is observed in at least 2 pawprints and also the edges are shared with the observations of the adjacent tiles.
The survey observation plan was organized in two stages: during the first year the tile is observed in five astronomical filters Z, Y, J, H and Ks separated by a few hours; then, in subsequent years, it is re-observed using the Ks band for variability studies. Only some tiles are observed in multi-band after the first year.

The dataset used in this work is the one presented in
\citet{cabral2020automatic}, which consists of {62} features extracted with \texttt{feets} \citep{cabral_fats_2018}
from light curves that were reconstructed from the photometric catalogs provided by the Cambridge Astronomical Survey Unit (CASU).

From the original dataset we selected 8 tiles located at different zones of the Bulge, as shown in Figure \ref{fig:used_tiles}. For each tile we extracted all the RRLs plus a uniform sample of $2000$ unclassified, unsaturated and non-faint sources (average magnitude between $12$ and $16.5$).
From these selections, sources with invalid values were removed, leaving the final dataset for this work, described in  Table~\ref{tab:sample}.

We choose to use a reduced dataset with around $2000$ sources for each tile in order to decrease dramatically the computational burden of our experiments. As shown in \citet{cabral2020automatic}, the use of a reduced dataset can lead to optimistic estimations of the accuracy of the detections, but our main objective is to find and characterize the features that best represent the differences between the tiles and not the accuracy of the detection of the RRLs.

\begin{figure*}
\centering
    \includegraphics[width=.7\textwidth]{figures/Fig_1.pdf}
    \caption{\label{fig:used_tiles}
    Map of the Bulge tiles of the  \textsc{VVV} survey over an extinction map (Extinction map adapted from and \citet{gonzalez2012reddening}). We highlight the tiles used in this work with red borders.}
\end{figure*}

\begin{table}
\caption{Total number of sources, RRL and sample taken in each tile used in this work.
}
    \centering
\begin{tabular}{l|rrr}
\toprule
 \textbf{Tile} &    \textbf{Total} &  \textbf{RRL} &  \textbf{Sample} \\
 b206 &   407720 &   47 &    2047 \\
 b214 &   376822 &   35 &    2034 \\
 b216 &   334773 &   43 &    2043 \\
 b261 &   735838 &  253 &    2252 \\
 b277 &   831323 &  430 &    2429 \\
 b278 &   857887 &  437 &    2436 \\
 b360 &  1029149 &  679 &    2669 \\
 b396 &   729671 &   15 &    2015 \\
\bottomrule
\end{tabular}
\label{tab:sample}
\end{table}

\subsection{Error measures}

We face two different binary classification problems in this work.
First we try to separate sources between two tiles, this instance is made in order to construct a Tile Classifier (hereafter TC) that allows to asses the relevance of the features. \rev{As we stated in the Introduction, we use the TC as an auxiliary method that allow us to detect which features are candidates to Drifting Features.}
Then, we build a Source Classifier (hereafter SC) that seeks to discriminate RRL sources from unknown sources.
In the first problem both classes are nearly balanced in all cases. On the other hand, as discussed in \citet{cabral2020automatic}, the identification of a few variable stars within a large set of unknown sources is usually a highly imbalanced problem, which generates several inconveniences such as those discussed in the recent work by \citet{hosenie2020imbalance} and requires specific error measures.

In the RRL detection problem (SC) we will call RRL samples as the positive class and the other sources as the negative class.
In the tile identification problem (TC) both classes (the two tiles) are equivalent, so we will arbitrarily call one of them as positive and the other as negative.
All positive samples (in this case, either a source or a tile) that are correctly identified by the classifier are called true positives (TP), otherwise if they are missed by the classifier they are called false negatives (FN).
Negative samples that are wrongly classified are called false positives (FP), and those correctly identified are called true negatives (TN).
Using a combination of these four outcomes, we can define two complementary performance measures, called Precision and Recall, which are adequate to deal with unbalanced problems. The Precision is defined as $\displaystyle{TP / (TP + FP)}$. It measures, for example, the fraction of real RRLs detected over all those retrieved by the classifier. The Recall, on the other hand, is defined as $\displaystyle{TP / (TP + FN)}$. It measures, in the same example, the fraction of all RRLs that are detected by the classifier.

Many classifiers can change their decisions outputs by adjusting the probability threshold that considers an observation to be positive or negative.
A high threshold increases the Precision and decreases the Recall since fewer cases are classified as positive, while a low threshold generates the opposite effect.
To evaluate Precision and Recall together, we consider the Precision-Recall curves, where we plot a set of pairs of values corresponding to different thresholds.
A curve that approaches the top-right corner is, in general, considered to represent a better classifier.

For balanced classification problems it is common to find more traditional metrics in the literature. According to that, for the tile identification problem we also use Accuracy $\displaystyle{(TP + TN) / (TP + FP + TN + FN)}$ and the area under the ROC curve (ROC-AUC) measures. The ROC curve is equivalent in concept to the Precision-Recall curve described before, and the area under it is a global measure of the performance of the classifier. The only difference between both curves is that a ROC curve that approaches the top-left corner represents a better classifier.

\subsection{Model Selection}
\label{section:model_select}

For the TC problem we evaluated four classifiers with diverse foundations: SVM (Support Vector Machine) with linear kernel \citep{vapnik2013nature}, SVM with Radial Basis Function (RBF) kernel, K-Nearest Neighbors \citep[KNN, ][]{mitchell1997machine} and Random Forest \citep[RF, ][]{breiman2001random}, \rev{all implementations from the Scikit-Learn Python package \citep{pedregosa2011scikit}}.

To determine the best hyper-parameters for every model, we executed a grid-search of all possible combinations of values for each hyper-parameter over a fixed list. We used a \rev{5 k-fold} setup on a dataset with tiles b278 and b261, using precision as performance measure. These tiles were chosen  because they are not extreme in location or in their balance, such as b396 or b220.

With this setup, we selected the following hyper-parameter values:
\begin{description}
    \item[SVM-Linear:] $C = 100$.
    \item[SVM-RBF:] $C = 100$ and $\gamma = 0.003$.
    \item[KNN:] $K = 56$ with a $Manhattan$ metric; also, the importance of the neighbor class was not weighted by distance.
    \item[RF:] We created $500$ decision trees with Information-Gain as metric, the maximum number of random selected features for each tree is the $0.5$ of the total number of features, and the minimum number of observations in each leaf is $5$.
\end{description}

\begin{table}
    \caption{
    Classification metrics of the SVM models (with kernel linear and RBF), RF and KNN  on the tiles b278 and b261 sources.}
    \centering
\begin{tabular}{l|rrr}
\toprule
      \textbf{Model} &  \textbf{Precision} &  \textbf{Recall} &    \textbf{AUC} \\
\midrule
 SVM-Linear &     0.8511 &  0.8511 & 0.9286 \\
    SVM-RBF &     0.8003 &  0.8003 & 0.8680 \\
         RF &     0.7707 &  0.7707 & 0.8548 \\
        KNN &     0.6973 &  0.6973 & 0.7685 \\
\bottomrule
\end{tabular}

\label{tab:model_select}
\end{table}

\begin{figure*}
    \centering
    \includegraphics[width=0.8\textwidth]{figures/Fig_2.pdf}
    \caption{\textsc{ROC} (left) and Precision-Recall (right) curves of the  SVM models (with linear and RBF kernels), RF and KNN, for the prediction of the tile of a given source, using 10--fold CV with tiles b278 and b261.}
    \label{fig:model_select}
\end{figure*}

Using the optimal values for the hyper-parameters we compared the four models on the same dataset using a \rev{10--fold cross--validation} setup. Table~\ref{tab:model_select} shows the corresponding results using the default threshold (0.5) for all models. For all three metrics considered (Precision, Recall and AUC) the SVM-Linear classifier clearly outperformed all the other classifiers. More importantly, Figure \ref{fig:model_select} shows the corresponding ROC and Precision-Recall curves, which show that SVM-Linear also outperforms the other methods for all possible thresholds.
Given these results, we selected SVM-Linear as the classifier for the tile identification problem.

For the SC problem \citet{cabral2020automatic} already determined that RF is the classifier with the best performance for our dataset and general experimental setup.
\subsection{Feature Selection}
\label{subsection:rfe}

Feature selection \citep{guyon2013introduction} is the process of extracting some subsets of features from the entire set in order to optimize the classification performance and/or the computational complexity of the problem.
We choose for this work the Recursive Feature Elimination algorithm \citep[RFE, ][]{guyon2002gene}. The method is widely adopted, characterized by its good performance and simplicity. As a backward selection method, RFE starts with all the features and sequentially eliminates the unimportant features using a recursive process.

RFE is integrated with a classification method, which provides at each step of the recursion the importance score of the features. RFE iteratively executes the underlying classifier and extracts the score for each variable, then the variable (or group of variables) with worse performance (according to the score) is eliminated.

The method typically ends when the desired (fixed) number of features to select is reached. Another possibility is to monitor a performance metric for the subsets (for example the accuracy on an independent validation set) and stop the recursion when the metric is optimal.

In this work we rely on the RFE implementation with k-fold Cross-Validation (RFECV) as stopping criteria \rev{also} provided by the Scikit-Learn package \citep{pedregosa2011scikit}.
RFECV produces $k$ replicated experiments ($k=5$ in our work), each one selecting features over ($k-1$) folds and monitoring the classification error ($1-accuracy$) over the remaining fold. Then it determines the number of features to select looking for the least average error throughout all the folds. In the last step, RFECV produces a final selection using all the dataset to select the features and stopping at the previously selected point.

It is worth mentioning that the classifier embedded in the RFE in the feature selection stage may be different from the eventual method in the final classification stage.

\section{Finding Drifting Features}
\label{section:drift}

As we stated in the Introduction, we propose to use ML methods to detect Drifting Features, looking for features that are useful to know the tile of origin of a given source (exclusively from features derived from the pawprint stack photometry, without any other Header Keyword Data). With this goal, in this first experiment we consider all the sources in each tile together (RRLs \& unknowns) and train classifiers to learn the tile of origin of each source and not its astronomic type.

We apply the RFECV method, as described in the previous section, on 28 binary classification problems, consisting each one in separating a different pair of tiles from the set of 8 tiles in our dataset. Thus, for each SC problem (for example, separating tiles b206 and b214) we obtain from RFECV a subset of selected features for that problem. Each subset has potentially a different length, as discussed before.

Figure \ref{fig:lless_heatmap} shows the number of features selected for each problem. It is evident that there are two different behaviours. In some cases RFE selects just a few features, as for example tile b216 with any other except b206; opposite, in some other cases the selected subset contains a high number of features (tile b206 against b216, b396 against b277 or b278, for example). But the number of selected features by itself is not relevant, what is more important to identify Drifting Features is how well they separate the two tiles.

We can arbitrarily divide the problem into ``few features" (4 or less selected features) and ``high number" (more than 4 selected features). Figure~\ref{fig:nl_acurves_1} shows ROC curves for the 12 ``high number" problems, and their relative locations in space. The figure shows curves for three classifiers: one trained with all the features in the dataset (All features), a second trained using only those selected by RFE, i.e. our candidates to Drifting Features, and a third one trained with those not selected by RFE (we call them the Stable subset). All the ``few features" subset cases (b216-b278 for example) produce trivial ROC curves saturated on the top-left corner for the three subsets, with $AUC>0.99$, which we do not show.

Analyzing the results, the first observation is that, as expected, the classifier trained with the Drifting Features (those selected by feature selection) is always very similar in performance to the one trained with all the features. This result is a confirmation that RFECV does its work, selecting a subset of features that is responsible for the separation of the classes. The second result is that the performance of the models for the ``few features" problems is clearly superior to the ones shown in the figure, i.e. the ``high number" problems. This implies that the two behaviours noted in Figure  \ref{fig:lless_heatmap} correspond to problems that are easy to solve, where the separation is almost perfect and can be done with a few features, and problems that are harder, where the tiles cannot be fully separated and RFE selects bigger subsets.

It is interesting to note the different response of the classifiers trained on the Stable subset on the easy and hard problems. In the hard problems RFECV selects a high number of features, which means that there are features with no considerable information about the tile of origin of the source. After several features were selected by RFECV, the rest of the features, the Stable subset, contains much less information about the origin and produces a classifier with low performance. For the easy problems, on the other hand, there are several features with much information about the tile of origin. Using just 2 to 4 features selected by RFECV is enough to produce an almost perfect classifier. The Stable subset in this case contains plenty of features with good information about the origin, producing a classifier with also almost perfect performance. If we take into account the relative positions of each pair of tiles for the Hard and Easy problems, there is not a definite pattern emerging. Most problems involving neighbour tiles, as for example b277, b278 and b261 are hard, and most problems involving tiles in the low-right region are easy.

\begin{figure}
    \centering
    \includegraphics[width=.85\columnwidth]{figures/Fig_3.pdf}
    \caption{
    Number of features selected by RFE for each binary TC problem. Each cell corresponds to the dataset including Tile A (Rows) and Tile B (Columns).}
    \label{fig:lless_heatmap}
\end{figure}

\begin{figure*}[tbh!]
    \centering

    \includegraphics[page=1, width=.7\textwidth]{figures/Fig_4.pdf}
    \caption{
        ROC curves for the combinations of two tiles selecting the more important features. Each panel has three curves: one using all the features (thick grey lines), the Drifting Features selected by RFE (dashed lines), and those considered stable, not selected by RFE (dotted lines).
    }
    \label{fig:nl_acurves_1}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=.8\textwidth]{figures/Fig_5.pdf}
    \caption{
        Total times that each feature was selected by RFECV over the 28 SC problems considered in this work. The different colors identify which of the three groups each feature belongs to: orange for those based on color, blue for those based on period, and white for those based only on magnitudes.
    }
    \label{fig:lless_freq}
\end{figure*}

Another relevant analysis for a feature selection method is which are the features that are selected in each case. The upper half of the Figure~\ref{fig:lless_freq} shows the number of times that each feature was selected by RFECV over the 28 TC problems, for all features that were selected at least two times (Table \ref{tab:rfe_no_limit_alt} in Appendix \ref{appendix:rfe} shows the list of features selected on each problem). Two particular features (\texttt{c89\_hk\_color} and \texttt{c09\_hk\_color}) were selected in almost all cases, while all the features related to pseudo-color \citep{cabral2020automatic} were the most frequently selected, appearing in at least half of the cases. This information suggests that color-related information in general is the most important characteristic to distinguish tiles.

\begin{figure}
    \centering
    \includegraphics[width=.41\textwidth]{figures/Fig_6.pdf}
    \caption{Same as Fig. \ref{fig:nl_acurves_1} for a fixed RFE selection of 10 features for the Drifting subset. Top row shows two Easy datasets while bottom row shows two Hard datasets.}
    \label{fig:l10_curve}
\end{figure}

The two very different behaviours of Hard and Easy problems will make harder the evaluation of the real influence of Drifting Features on the TC problems, because deleting only two features or half of the features will lead to diverse scenarios. To make an easier and fairer comparison we changed the feature selection method, using RFE with a fixed number of ten selected features.

The list of these selected features and their frequency can be observed in the bottom half of Figure~\ref{fig:lless_freq}.
Only 15 features were selected in total over the 28 TC problems, from whom the 11 more relevant are related to color, probably with a high dependence with the location of the tile.

The overall performance of the classifiers trained with the Full, Drifting and Stable subsets, for some exemplars of Easy and Hard datasets, can be seen in Figure~\ref{fig:l10_curve}. For the Easy problems (bottom row) we are now using less features in the Stable subset, leading to a lower performance. On the opposite side, for the hard problems (top row) we are now using more features in the Stable subset, leading to a clear increase in its performance. The rest of the TC problems show the same type of result (data not shown).

Using a fixed selection of features with RFE we obtain in all cases a subset of 10 features (the Drifting subset) that can discriminate the tile of origin with high accuracy, and another subset (Stable) with much less information about the tile of origin of the sources. We analyze in the following the impact of these subsets on the SC problems.

\section{Evaluation of the influence of Drifting Features}
\label{section:evaluation}

In this section we evaluate the influence of the Drifting subsets selected in the previous step on the SC problems, i.e., to discriminate unknown sources and RRL variable stars.

For each pair of tiles we have three datasets, one with all the features (Full), a second with only 10 Drifting Features selected by RFE and finally one with the remaining features, the Stable subset. Different from the previous problem we have now, for the SC problem, two possibilities for each pair of tiles: first we train our classifiers in one of the tiles and search for RRLs in the other, and second we invert the tiles, training classifiers in the second tile of the pair and looking for RRLs in the first tile. Thus, for each one of the 56 SC problems we train corresponding RF classifiers and obtain three PR curves for the Full, Drifting and Stable classifiers.

The complete results are presented in Appendix~\ref{appendix:evaluation}, while a summary of some representative cases can be seen in Figure~\ref{fig:nl_acurves_1_PR}. A first result is that, clearly, the Drifting subset shows lower performance in all cases. This was expected as most Drifting Features are related to color and \citet{cabral2020automatic} demonstrated that color alone cannot clearly identify RRLs. More interesting, if we compare the performance of the Full datasets with the Stable datasets, we can see that there is no clear advantage in eliminating the Drifting Features from the datasets. Full and Stable curves are almost similar in all cases. Differences are small and there is no clear pattern of when eliminating Drifting Features will improve the performance of the ML methods.

\begin{figure*}
    \centering
    \includegraphics[width=.95\textwidth]{figures/Fig_7.pdf}
    \caption{Precision-Recall curves for the SC problems. We show results for six combinations of train-test tiles
    using three classifiers trained with the Full, Drifting and Stable subsets of features.}
    \label{fig:nl_acurves_1_PR}
\end{figure*}

\section{Discussion}
\label{section:conclusions}
In this work we introduced and discussed the concept of Drifting Features, related with small changes in the properties measured by the features, which can potentially harm the result of ML methods in astronomy. Using the identification of RRLs in VVV as a working example, we introduced a method for detecting Drifting Features, using an indirect ML method. We forced a classifier to learn the tile of origin of diverse sources, and select the features more relevant to this task as candidates to Drifting Features. We showed that this method can efficiently identify a reduced set of features that contains useful information about the tile of origin of the sources. We also showed that, for our particular example of detecting RRLs in VVV, Drifting Features are mostly related to color. On the other hand, we showed in the Section~\ref{section:evaluation} that, even if we have a clear set of Drifting Features in our problem, they are almost harmless for the identification of RRLs.

As a future work we will explore the influence of Drifting Features on the detection of other types of variable sources and other large scale ML experiments. We will also explore a different way of setting the number of selected features by RFE, considering all features that are relevant to the problem and not only the subset that shows the best performance for some metric or a fixed-length subset.

\begin{acknowledgements}
The authors would like to thank to their families and friends, and also IATE astronomers for useful comments and suggestions.
This work was partially supported by the Consejo Nacional
de Investigaciones Cient\'ificas y T\'ecnicas (CONICET, Argentina) and the Secretar\'ia de Ciencia y Tecnolog\'ia de la Universidad Nacional de C\'ordoba (SeCyT-UNC, Argentina).
J.B.C, are supported by a fellowship from CONICET.
Some processing was achieved with Argentine VO (NOVA) infrastructure, for which the authors express their gratitude.
We gratefully acknowledge data from the ESO Public Survey program ID 179.B-2002 taken with the VISTA telescope and products from the Cambridge Astronomical Survey Unit (CASU).
J.B.C. thanks to Maren Hempel by creating the template for the creation of the template on which the figure~\ref{fig:used_tiles} is based, and finally Bruno Sánchez and Martín Beroiz for the continuous support and friendship.
This research has made use of the
\url{http://adsabs.harvard.edu/}, Cornell University xxx.arxiv.org repository, adstex (\url{https://github.com/yymao/adstex}), astropy and the Python programming language.
\end{acknowledgements}

\bibliographystyle{aa}