\section{Introduction}
Human motion prediction aims to forecast the sequence of future poses of a person given past observations of such poses. To achieve this, existing methods typically rely on recurrent neural networks (RNNs) that encode the person's motion~\cite{martinez2017human,gui2018adversarial,walker2017pose,kundu2018bihmp,barsoum2018hp,pavllo2019modeling,pavllo2018quaternet}. While they predict reasonable motions, RNNs are deterministic models and thus cannot account for the highly stochastic nature of human motion; given the beginning of a sequence, multiple, diverse futures are plausible. To correctly model this, it is therefore critical to develop algorithms that can learn the \emph{multiple modes} of human motion, even when presented with only deterministic training samples.

Recently, several attempts have been made at modeling the stochastic nature of human motion~\cite{yan2018mt,barsoum2018hp,walker2017pose,kundu2018bihmp,lin2018human}. These methods rely on sampling a random vector that is then combined with an encoding of the observed pose sequence. In essence, this combination is similar to the conditioning of generative networks; the resulting models aim to generate an output from a random vector while taking into account additional information about the content.

\begin{figure}
    \centering
    \small
        \input{ignorance.tex}
    \caption{Diversity of $K$ RNN decoder inputs,
    generated with $K=50$ different random vectors. We report the mean diversity over $N=50$ samples and the corresponding standard deviation.
    }

    \label{fig:baseline_RHP}
\end{figure}

Here, we argue that, while standard conditioning strategies may be
effective for many tasks, as in~\cite{yan2016attribute2image,kulkarni2015deep,esser2018variational,engel2017latent,bao2017cvae,larsen2015autoencoding}, they are ill-suited for motion prediction. The reason is the following: In other tasks, the conditioning variable only provides auxiliary information about the output to produce,
such as the fact that a generated face should be smiling. By contrast, in motion prediction, it typically contains the core signal to produce the output, i.e., the information about the previous poses. Since the prediction model is trained using deterministic samples, it can then simply learn to ignore the random vector and still produce a meaningful output based on the conditioning variable only. In other words, the model can ignore the root of variations, and thus essentially become deterministic.

This problem was discussed in~\cite{bowman2015generating} in the context of text generation, and we identified it in our own motion prediction experiments. As evidence, we plot in Fig.~\ref{fig:baseline_RHP} the diversity of the representations used as input to the RNN decoders of~\cite{yan2018mt} (\texttt{LHP}) and~\cite{barsoum2018hp} (\texttt{RHP}), two state-of-the-art methods that are closest in spirit to our approach. Here, diversity is measured as the average pairwise distance across the $K$ representations produced for a single series of observations. We report the mean diversity over 50 samples and the corresponding standard deviation. As can be seen from the figure, the diversity of \texttt{LHP} and \texttt{RHP} decreases as training progresses, thus supporting our observation that the models learn to ignore the
perturbations.

In this paper, we introduce a simple yet effective approach to counteracting this loss of diversity and thus to generating truly diverse future pose sequences. At the heart of our approach lies the idea of \emph{Mix-and-Match} perturbations: Instead of combining a noise vector with the conditioning variables in a deterministic manner,
we randomly select and perturb a subset of these variables. By randomly changing this subset at every iteration, our strategy prevents training from identifying the root of variations and forces the model to take it into account in the generation process. As a consequence, and as evidenced by the black curve in Fig.~\ref{fig:baseline_RHP}, which shows an increasing diversity as training progresses, our approach produces not only high-quality predictions but also truly diverse ones.

In short, our contributions are \textbf{(1)} a novel way of imposing diversity into conditional VAEs, called \emph{Mix-and-Match perturbations}; \textbf{(2)} a new motion prediction model capable of generating multiple likely future pose sequences from an observed motion; and \textbf{(3)} a new evaluation metric for quantitatively measuring the quality and the diversity of generated motions, thus facilitating the comparison of different stochastic approaches.

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

Most motion prediction approaches are based on \emph{deterministic} models~\cite{pavllo2018quaternet,pavllo2019modeling,gui2018adversarial,jain2016structural,martinez2017human,gui2018few,fragkiadaki2015recurrent,ghosh2017learning}, casting motion prediction as a regression task where only one outcome is possible given the observations. While this may produce accurate predictions, it fails to reflect the stochastic nature of human motion, where multiple plausible outcomes can be highly likely for a single given series of observations. Modeling this diversity is the topic of this paper, and we therefore focus the discussion below on the other methods that have attempted to do so.

The general trend to incorporate variations in the predicted motions consists of combining information about the observed pose sequence with a random vector. In this context, two types of approaches have been studied: The techniques that directly incorporate the random vector into the RNN decoder and those that make use of an additional Conditional Variational Autoencoder (CVAE)~\cite{sohn2015learning}.

In the first class of methods,~\cite{lin2018human} samples a random vector $z_t\sim\mathcal{N}(0,I)$ at each time step and adds it to the pose input to the RNN decoder. By relying on different random vectors at each time step, however, this strategy is prone to generating discontinuous motions. To overcome this,~\cite{kundu2018bihmp} makes use of a single random vector to generate the entire sequence. This vector is both employed to alter the initialization of the decoder and concatenated with a pose embedding at each iteration of the RNN. By relying on concatenation, these two methods contain parameters that are specific to the random vector, and thus give the model the flexibility to ignore this information.
In~\cite{barsoum2018hp}, instead of using concatenation, the random vector is added to the hidden state produced by the RNN encoder. While addition prevents having parameters that are specific to the random vector, this vector is first transformed by multiplication with a parameter matrix, and thus can again be zeroed out so as to remove the source of diversity, as we observe empirically in Section~\ref{sec:experiments_QD}. In our experiments, we will refer to this method as \texttt{RHP}, for random hidden state perturbation.

The second category of stochastic methods introduce an additional CVAE between the RNN encoder and decoder. This allows them to learn a more meaningful transformation of the noise, combined with the conditioning variables, before passing the resulting information to the RNN decoder. In this context,~\cite{walker2017pose} proposes to directly use the pose as conditioning variable. As will be shown in our experiments, while this approach is able to maintain some degree of diversity, albeit less than ours, it yields motions of lower quality because of its use of independent random vectors at each time step. In our experiments, we will refer to this method as \texttt{LPP}, for learned pose perturbation.  In~\cite{butepage2018anticipating}, an approach similar to that of~\cite{walker2017pose} is proposed, but with one CVAE per limb. As such, this method suffers from the same discontinuity problem as~\cite{walker2017pose,lin2018human}. Finally, instead of perturbing the pose, the recent work~\cite{yan2018mt} uses the RNN decoder hidden state as conditioning variable in the CVAE, concatenating it with the random vector. While this approach generates high-quality motions, it suffers from the fact that the CVAE decoder gives the model the flexibility to ignore the random vector. In the remainder of the paper, we will refer to this method as \texttt{LHP}, for learned hidden state perturbation.

Ultimately, both classes of methods suffer from the fact that they allow the model to ignore the random vector, thus relying entirely on the conditioning information to generate future poses. Here, we introduce an effective way to maintain the root of diversity by randomizing the combination of the random vector with the conditioning variable.

\section{Proposed Method}
In this section, we first introduce our \emph{Mix-and-Match} approach to introducing diversity in CVAE-based motion prediction. We then describe the motion prediction architecture we used in our experiments and propose a novel evaluation metric to quantitatively measure the diversity and  quality of generated motions.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{mixandmatch_new.pdf}
    \caption{Mix-and-Match perturbation. \textbf{(Top)} Illustration of the \textit{Sampling} operation (left) and of the \textit{Resampling} one (right). Given a sampling rate $\alpha$ and a vector length $L$, the Sampling operation samples $\left\lceil\alpha L\right\rceil$ indices, say $\mathcal{I}$. The complementary, unsampled, indices are denoted by $\mathcal{\bar{I}}$. Then, given two $L$-dimensional vectors and the corresponding  $\left\lceil\alpha L\right\rceil$ and $\left\lfloor(1-\alpha)L\right\rfloor$ indices, the Resampling operation mixes the two vectors to form a new $L$-dimensional one. \textbf{(Middle)} Example of Mix-and-Match perturbation. \textbf{(Bottom)} Example of perturbation by concatenation, as in~\cite{yan2018mt}. Note that, in Mix-and-Match perturbations, sampling is stochastic; the indices are sampled uniformly randomly for each mini-batch. By contrast, in~\cite{yan2018mt}, sampling is deterministic, and the indices in $\mathcal{I}$ are fixed and correspond to $\mathcal{I}=\{1,\dots,\frac{L}{2}\}$.
    }
    \label{fig:mix_match}
\end{figure}

\subsection{Mix-and-Match Perturbation}
\label{sec:mix_match}

The main limitation of prior work in the area of stochastic motion modeling, such as~\cite{walker2017pose,barsoum2018hp,yan2018mt}, lies in the way they fuse the random vector with the conditioning variable, i.e., RNN hidden state or pose, which causes the model to learn to ignore the  randomness and solely exploit the deterministic conditioning information to generate motion~\cite{walker2017pose,barsoum2018hp,yan2018mt}.
To overcome this, we propose to make it harder for the model to decouple the random variable from the deterministic information. Specifically, we observe that the way the random variable and the conditioning one are combined in existing methods is deterministic. We therefore propose to make this process stochastic.

Similarly to~\cite{yan2018mt}, we propose to make use of the hidden state as the conditioning variable and generate a perturbed hidden state by combining a part of the original hidden state with the random vector. However, as illustrated in Fig.~\ref{fig:mix_match}, instead of assigning predefined, deterministic indices to each piece of information, such as the first half for the hidden state and the second one for the random vector, we assign the values of hidden state to \textit{random} indices and the random vector to the complementary ones.

More specifically, as depicted in Fig.~\ref{fig:mix_match}, a mix-and-match perturbation takes two vectors of size $L$ as input, say $h_t$ and $z$, and combines them in a stochastic manner. To this end, it relies on two operations. The first one, called \textit{Sampling}, chooses $\left\lceil\alpha L\right\rceil$ indices uniformly at random among the $L$ possible values, given a sampling rate $0\leq \alpha \leq 1$. Let us denote by ${\cal I} \subseteq\{1,\ldots,L\}$, the resulting set of indices and by ${\cal \bar{I}}$ the complementary set. The second operation, called \textit{Resampling}, then creates a new $L$-dimensional vector whose values at indices in $\mathcal{I}$ are taken as those at corresponding indices in the first input vector and the others at the complementary indices in the second input vector. Note that, the second vector can also have dimension $\left\lfloor(1-\alpha)L\right\rfloor$, and its values be divided among the remaining indices of the output vector.

\begin{figure*}[t]
    \centering
    \includegraphics[width=\textwidth]{train_inference_new.pdf}
    \caption{Overview of our approach. \textbf{(Top): Overview of the model during inference}. During inference, given past information and a random vector sampled from a Normal distribution, the model generate new motions.
    \textbf{(Bottom): Overview of the model during training.} During training, we use a future pose autoencoder with a CVAE  between the encoder and the decoder. The RNN encoder-decoder network mapping the past to the future then aims to generate good conditioning variables for the CVAE. }
    \label{fig:stochastic}
\end{figure*}

\subsection{M\&M Perturbation for Motion Prediction}
Let us now describe the way we use our mix-and-match perturbation strategy for motion prediction. To this end, we first discuss the network we rely on during inference, and then explain our training strategy.

\paragraph{Inference.} The high-level architecture we use at inference time is depicted by Fig.~\ref{fig:stochastic} (Top). It consists of an RNN encoder that takes $t$ poses $x_{1:t}$ as input and outputs an $L$-dimensional hidden vector $h_t$. A random $\left\lceil\alpha L\right\rceil$-dimensional portion of this hidden vector, $h_t^\mathcal{I}$, is then combined with an $\left\lfloor(1-\alpha)L\right\rfloor$-dimensional random vector $z\sim\mathcal{N}(0,I)$ via our mix-and-match perturbation strategy.
The resulting $L$-dimensional output is passed through a small neural network (i.e., \textit{ResBlock2} in Fig.~\ref{fig:stochastic}) that reduces its size to $\left\lceil\alpha L\right\rceil$, and then fused with the remaining $\left\lfloor(1-\alpha)L\right\rfloor$-dimensional portion of the hidden state, $h_t^\mathcal{\bar{I}}$. This, in turn, is passed through the VAE decoder to produce the final hidden state $h_z$, from which the future poses $x_{t+1:T}$ are obtained via the RNN decoder.

\paragraph{Training.} During training, we aim to learn both the RNN parameters and the CVAE ones. Because the CVAE is an \emph{auto}encoder, it needs to take as input information about future poses. To this end, we complement our inference architecture with an additional RNN future encoder, yielding the training architecture depicted in Fig.~\ref{fig:stochastic} (Bottom). Note that, in this architecture, we incorporate an additional mix-and-match perturbation that fuses the hidden state of the RNN past encoder $h_t$ with that of the RNN future encoder $h_T$ and forms $h_{tT}^p$. This allows us to condition the VAE encoder in a manner similar to the decoder.
Note that, for each mini batch, we use the same set of sampled indices for all mix-and-match perturbation steps throughout the network. Furthermore, following the standard CVAE strategy, during training, the random vector $z$ is sampled from a distribution $\mathcal{N}(\mu_\theta(x), \Sigma_\theta(x))$, whose mean $\mu_\theta(x)$ and covariance matrix $\Sigma_\theta(x)$ are produced by the CVAE encoder with parameters $\theta$. This is done by the technique of~\cite{kingma2013auto}, which computes $z_p$ as,

\begin{align}
z_p=\mu_\theta(x) + \Sigma_\theta(x)^{\frac{1}{2}} \epsilon\;,
\label{eq:reparameterization}
\end{align}
where $\epsilon\sim\mathcal{N}(0,I)$. Note that, during inference, $z_p=\epsilon$ since we do not have access to $x$, hence to $\mu_\theta(x)$ and $\Sigma_\theta(x)$.

To learn the parameters of our model, we rely on the availability of a dataset $D=\{X_1, X_2, ..., X_N\}$ containing $N$ videos $X_i$ depicting a human performing an action. Each video consists of a sequence of $T$ poses, $X_i=\{x_i^1, x_i^2, ..., x_i^T\}$, and each pose comprises $J$ joints forming a skeleton, $x_i^t=\{x_{i, 1}^t, x_{i, 2}^t, ..., x_{i, J}^t\}$. The pose of each joint is represented as a quaternion. Given this data, we train our model by minimizing a loss function of the form
\begin{align}
\mathcal{L} = \frac{1}{N} \sum_{i=1}^N\Big(\mathcal{L}_{rot}(X_i) + \mathcal{L}_{skl}(X_i)\Big) + \lambda\mathcal{L}_{prior}\;.
\label{eq:stochastic_loss}
\end{align}

The first term in this loss compares the output of the network with the ground-truth motion using the squared loss. That is,
\begin{align}
\mathcal{L}_{rot}(X_i) = -\frac{1}{T-t}\frac{1}{4J}\sum_{k=t+1}^T\sum_{j=1}^{J}{\|\hat{x}^{k}_{i,j} - x^{k}_{i,j}\|}^2\;,
\label{eq:loss_quat}
\end{align}
where $\hat{x}^{k}_{i,j}$ is the predicted 4D quaternion for the $j^{th}$ joint at time $k$ in sample $i$, and $x^{k}_{i,j}$ the corresponding ground-truth one. The main weakness of this loss is that it treats all joints equally. However, when working with angles, some joints have a much larger influence on the pose than others. For example, because of the kinematic chain, the pose of the shoulder affects that of the rest of the arm, whereas the pose of the wrists has only a minor effect.

To take this into account, we define our second loss term as the error in 3D space. That is,
\begin{align}
\mathcal{L}_{skl}(X_i) = -\frac{1}{T-t}\frac{1}{J\times 3}\sum_{k=t+1}^{T}\sum_{j=1}^{J} \|\hat{p}^{k}_{i,j} - p^{k}_{i,j}\|^2\;,
\label{eq:loss_skl}
\end{align}
where $\hat{p}^{k}_{i,j}$ is the predicted 3D position of joint $j$ at time $k$ in sample $i$ and $p^{k}_{i,j}$ the corresponding ground-truth one. These 3D positions can be computed using forward kinematics, as in~\cite{pavllo2018quaternet,pavllo2019modeling}. Note that, to compute this loss, we first perform a global alignment of the predicted pose and the ground-truth one by rotating the root joint to face [0, 0, 0].

Finally, following standard practice in training VAEs, we define our third loss term as the KL divergence
\begin{align}
    \mathcal{L}_{prior} = -KL\Big(\mathcal{N}(\mu_\theta(x), \Sigma_\theta(x))) \Big\| \mathcal{N}(0,I)\Big)\;.
    \label{eq:kl}
\end{align}

In practice, since our VAE appears within a recurrent model, we weigh $\mathcal{L}_{prior}$ by a function $\lambda$ corresponding to the KL annealing weight of~\cite{bowman2015generating}. We start from $\lambda=0$, forcing the model to encode as much information in $z$ as possible, and gradually increase it to $\lambda=1$, following a logistic curve.

\paragraph{Curriculum Learning of Variation.} The parameter $\alpha$ in our mix-and-match perturbation scheme determines a trade-off between stochasticity and motion quality. The larger $\alpha$, the larger the portion of the original hidden state that will be perturbed. Thus, the model incorporates more randomness and less information from the original hidden state. As such, given a large $\alpha$, it becomes harder for the model to deliver motion information from the observation to the future representation since a large portion of the hidden state is changing randomly. In particular, we observed that training becomes unstable if we use a large $\alpha$ from the beginning, with the motion-related loss terms fluctuating while the prior loss $\mathcal{L}_{prior}$ quickly converges to zero.

To overcome this while still enabling the use of sufficiently large values of $\alpha$ to achieve high diversity, we introduce the curriculum learning strategy depicted by Fig.~\ref{fig:curriculum}. In essence, we initially select $\left\lceil\alpha L\right\rceil$ indices in a deterministic manner and gradually increase the randomness of these indices as training progresses. More specifically, given a set of $\left\lceil\alpha L\right\rceil$ indices, we replace $c$ indices from the sampled ones with the corresponding ones from the remaining $\left\lfloor(1-\alpha)L\right\rfloor$ indices. Starting from $c=0$, we gradually increase $c$ to the point where all $\left\lceil\alpha L\right\rceil$ indices are sampled uniformly randomly. More details, including the pseudo-code of this approach, are provided in the supplementary material. This strategy helps the motion decoder to initially learn and incorporate information about the observations (as in~\cite{yan2018mt}), yet, in the long run, still prevents it from ignoring the random vector.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{curriculum.pdf}
    \caption{Example of curriculum perturbation of the hidden state. At the beginning of training, the perturbation occurs in a deterministic portion of the hidden state. As training progresses, we gradually, and randomly, spread the perturbation to the rest of the hidden state. This continues until the indices to perturb are uniformly randomly sampled.
    }
    \label{fig:curriculum}
\end{figure}

\subsection{Quality and Diversity Metrics}
\label{sec:eval}
When dealing with multiple plausible motions, or in general diverse solutions to a problem, evaluation is a challenge. The standard metrics used for deterministic motion prediction models are ill-suited to this task, because they typically compare the predictions to the ground truth, thus inherently penalizing diversity. For multiple motions, two aspects are important: the \emph{diversity} and the \emph{quality}, or realism, of each individual motion.
Prior work typically evaluates these aspects via human judgement. While human evaluation is highly valuable, and we will also report human results, it is very costly and time-consuming. Here, we therefore introduce two metrics that facilitate the quantitative evaluation of both quality and diversity.

To measure the quality of generated motions, we propose to rely on a binary classifier
trained to discriminate real (ground-truth) samples from fake (generated) ones. The accuracy of this classifier on the test set is thus inversely proportional to the quality of the generated motions. In other words, high-quality motions are those that are not distinguishable from real ones. Note that we do not advocate for adversarial training of our approach. That is, we do not define a loss based on this classifier when training our model.

To measure the diversity of the generated motions, a naive approach would consist of relying on the distance between the generated motion and a reference one. However, generating identical motions that are all far from the reference one would therefore yield a high value, while not reflecting diversity. To prevent this, we therefore propose to make use of the average distance between all pairs of generated motions.

\section{Experiments}

Let us now evaluate the effectiveness of our approach at generating multiple plausible motions. To this end, we use Human3.6M~\cite{h36m_pami}, the largest publicly available motion capture dataset.  Below, we first exploit the metrics introduced in Section~\ref{sec:eval} to compare the quality and diversity of the results of our approach with those obtained by state-of-the-art methods that produce multiple motions~\cite{yan2018mt, walker2017pose,barsoum2018hp}. We then compare our results to the state-of-the-art deterministic motion prediction techniques for long-term motion prediction using standard metrics.

\paragraph{Implementation details.} The motion encoders and decoders in our model are single layer GRU~\cite{cho2014learning} networks, comprising 1024 hidden units each. For the decoders, we use a teacher forcing technique~\cite{williams1989learning} to decode motion. At each time-step, the network chooses with probability $P_{tf}$ whether to use its own output at the previous time-step or the ground-truth pose as input. We initialize $P_{tf}=1$, and decrease it linearly at each training epoch such that, after a certain number of epochs, the model becomes completely autoregressive, i.e., uses only its own output as input to the next time-step.
We train our model on a single GPU with the Adam optimizer~\cite{kingma2014adam} for 100K iterations. We use a learning rate of 0.001 and a mini-batch size of 64. To avoid exploding gradients, we use the gradient-clipping technique of~\cite{pascanu2013difficulty} for all layers in the network. We implemented our model using the Pytorch framework of~\cite{paszke2017automatic}.

\subsection{Evaluating Quality and Diversity}
\label{sec:experiments_QD}
Quantitative evaluation of a qualitative task is very challenging. While the ideal case is reporting the (log-)likelihood on a held-out set of samples, in (non-probabilistic) decoder-based generative models this is not possible. An alternative is using non-parametric kernel density estimates (KDE), only via samples, however, KDE is only well-suited for very low dimensional data space.
Evaluating against one GT motion (i.e., one sample from multi-modal distribution) can lead to a high score for one sample while penalizing other plausible modes. This behaviour is undesirable since it cannot differentiate a multi-modal solution with a good, but uni-modal one. Note, there exist some metrics~\cite{yan2018mt} to evaluate motions, however, they do not reflect the quality of a prediction, but how likely ground-truth future motions are with the given model. Moreover, as discussed, the metrics in~\cite{yan2018mt} only evaluate quality given one single groundtruth. While the groundtruth has high quality, there exist multiple high quality continuations of an observation, which our metric accounts for.
As discussed in Section~\ref{sec:eval}, we evaluate both the quality and diversity of the predicted motions. Note, these two metrics should be considered together, since each one taken separately does not provide a complete picture of how well a model can predict \emph{multiple plausible} future motions. For example, a model can generate diverse but unnatural motions, or, conversely, realistic but identical motions.

We compare our \emph{Mix-and-Match} approach with the different means of imposing variation in motion prediction discussed in Section~\ref{sec:related_work}, i.e., concatenating the hidden state to a learned latent variable (\texttt{LHP})~\cite{yan2018mt}, concatenating the pose to a learned latent variable at each time-step (\texttt{LPP})~\cite{walker2017pose}, and adding a (transformed) random noise to the hidden state (\texttt{RHP})~\cite{barsoum2018hp}. For the comparison to be fair, we use 16 frames (i.e., 640ms) as observation to generate the next 60 frames (i.e., 2.4sec) for all baselines. All models are trained with the same motion representation, backbone network, and  losses, except for \texttt{RHP} which cannot make use of $\mathcal{L}_{prior}$.

\begin{figure}
    \centering
    \includegraphics[width=.4\textwidth]{QualityClassifier.pdf}
    \caption{Architecture of the quality binary classifier.}
    \label{fig:quality_classifier}
\end{figure}

\begin{figure}
    \centering
    \begin{tabular}{c}
    \input{QD_plot_div.tex} \\
    \input{QD_plot_qual.tex}

    \end{tabular}
    \caption{Quality and diversity evaluation. Our approach outperforms the baselines in terms of diversity while preserving a high quality, especially late in the training progress.}
    \label{fig:QD_plot}
\end{figure}

To evaluate quality, as discussed in Section~\ref{sec:eval}, we use a recurrent binary classifier whose task is to determine whether a sample comes from the ground-truth data or was generated by the model.
As depicted by Fig.~\ref{fig:quality_classifier}, the model is based on a single layer GRU network with 1024 hidden units to process the motion, followed by a three-layer fully connected network (with 512, 128 and 1 units, respectively) with ReLU non-linearity in between and a sigmoid non-linearity for binary classification. We train such a classifier for each method, using 25K samples generated at different training steps together with 25K real samples, forming a binary dataset of 50K motions for each method. We use stochastic gradient descent for 5K iterations, with a mini-batch size of 256, a learning rate of 0.01 and a momentum of 0.9. To evaluate diversity, as discussed in Section~\ref{sec:eval}, we compute the mean Euclidean distance from each motion to all other $K-1$ motions when generating $K=50$ motions. Furthermore, we also performed a human evaluation to measure the quality of the motions generated by each method. To this end, we asked eight users to rate the quality of 50 motions generated by each method, for a total of 200 motions. The ratings were defined on a scale of 1-5, 1 representing a low-quality motion and 5 a high-quality, realistic one. We then scaled the values to the range 0-50 to make them comparable with those of the binary classifier.

The results of the metrics of Section~\ref{sec:eval} are provided in  Fig.~\ref{fig:QD_plot} and those of the human evaluation in Fig.~\ref{fig:human_eval}. Below, we analyze the results of the different models.

\paragraph{\texttt{LHP}.} As can be seen from Fig.~\ref{fig:QD_plot}, \texttt{LHP} tends to ignore the random variable $z$, thus ignoring the root of variation. As a consequence, it achieves a low diversity, much lower than ours, but produces samples of high quality, albeit almost identical. Note that this decrease in diversity occurs after 16K iterations, indicating that the model takes time to identify the part of the hidden state that contains the randomness. Nevertheless, at iteration 16K, prediction quality is low, and thus one could not simply stop training at this stage. Note that the lack of diversity of \texttt{LHP} is also evidenced by Fig.~\ref{fig:baseline_RHP}. To further confirm it, we performed an additional experiment where, at test time, we sampled each element of the random vector independently from $\mathcal{N}(50, 50)$ instead of from the prior $\mathcal{N}(0, I)$. This led to neither loss of quality nor increase of diversity of the generated motions. As can be verified in Fig.~\ref{fig:human_eval}, where \texttt{LHP} appears in a region of high quality but low diversity, the results of human evaluation match those of our classifier-based quality metric.

\begin{figure}
    \centering
    \input{qd_ch.tex}
    \caption{Human (H) and classifier-based (C) evaluation of quality for different methods. We plot diversity vs quality. A good method should fall into the top-right part of the plot, i.e., have high quality and diversity. Only our approach, for both human and classifier-based evaluation, satisfies this criterion.
    Real motions (blue circle) are deterministic, i.e., one future per observation, and have 0 diversity. However, their quality is optimal, i.e., 50\%.
    Note that the human and classifier-based results depict the same behavior.
    }
    \label{fig:human_eval}
\end{figure}

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{poses.png}
    \caption{Qualitative visualization of diverse motions generated by our model and by the baselines. Each block of columns shows the results for one observation (the first three poses of each sequence). The first row corresponds to the ground-truth motion and the other rows illustrate  multiple motions generated by each method (better seen when zoomed in).}
    \label{fig:qaulitative}
\end{figure*}

\begin{figure}
    \centering
    \begin{tabular}{@{ }c@{ }c}
        \input{alpha_q.tex} & \input{alpha_d.tex}
    \end{tabular}

    \caption{Quality and the diversity of the motions generated with our approach as a function of $\alpha$.
    Note that with $\alpha>0.7$, diversity increases significantly, but this diversity is the result of poor-quality motions.}
    \label{fig:tradeoff}
\end{figure}

\begin{figure*}[t]
    \centering
    \begin{tabular}{@{ }c@{ }c@{ }c@{ }c}
         \input{walking_long.tex} &
         \input{eating_long.tex} &
         \input{smoking_long.tex} &
         \input{discussion_long.tex} \\

         \input{baselines_walking.tex} &
         \input{baseline_eating.tex} &
         \input{baseline_smoking.tex} &
         \input{baseline_discussion.tex} \\

    \end{tabular}
    \caption{Mean angle error (MAE) for the Human 3.6M actions commonly used to report long-term prediction results. \textbf{(Top)} We compare the best of $K$ motions generated by our approach with the state-of-the-art deterministic baselines. Note that, while our approach does exploits knowledge of the action during neither training nor inference (unlike some of the baselines), it performs on par with the state-of-the-art deterministic baselines. \textbf{(Bottom)} We compare our approach with the state-of-the-art stochastic baselines. Note that the results for each stochastic baseline were obtained from the best of $K$ generated motions.
    }
    \label{fig:long_term}
\end{figure*}

\begin{figure*}[!t]
    \centering
    \includegraphics[width=0.95\textwidth]{effectK.png}
    \caption{Effect of $K$ in the MAE of our actions of the Human3.6M dataset. Note, the bold black one is the best of $K=50$ motions, and the shaded area indicates the region between best of $K=1$ and $K=500$.}
    \label{fig:k_effect}
\end{figure*}

\paragraph{\texttt{RHP}.} As for \texttt{LHP}, Fig.~\ref{fig:QD_plot} evidences the limited diversity of the motions produced by \texttt{RHP} despite its use of random noise during inference. Note that the authors of~\cite{barsoum2018hp}  mentioned in their paper that the random noise was added to the hidden state. Only by studying their publicly available code\footnote{\texttt{https://github.com/ebarsoum/hpgan}} did we understand the precise way this combination was done. In fact, the addition relies on a parametric, linear transformation of the noise vector. That is, the perturbed hidden state is obtained as
\begin{align}
    h_{perturbed} = h_{original} + W^{z\rightarrow h} z\;.
\end{align}
Because the parameters $W^{z\rightarrow h}$ are \emph{learnt}, the model has the flexibility to ignore $z$, which causes the behavior observed in Figs.~\ref{fig:QD_plot} and~\ref{fig:baseline_RHP}. Note that the authors of~\cite{barsoum2018hp} acknowledged that, despite their best efforts, they noticed very little variations between predictions obtained with different $z$ values. Since the perturbation is ignored, however, the quality of the generated motions is high. By depicting \texttt{RHP} in a region of high quality but low diversity, the human evaluation results in Fig.~\ref{fig:human_eval} again match those of our classifier-based quality metric.

\paragraph{\texttt{LPP}.} As can be seen in Fig.~\ref{fig:QD_plot}, \texttt{LPP} produces motions with higher diversity than \texttt{LHP} and \texttt{RHP}, but of much lower quality. The main reason behind this is that the random vectors that are concatenated to the poses at each time-step are sampled independently of each other, which translates to discontinuities in the generated motions. This problem might be mitigated by sampling the noise in a time-dependent, autoregressive manner, as in~\cite{videoFlow} for video generation. Doing so, however, goes beyond the scope of our analysis. When it comes to human evaluation, Fig.\ref{fig:human_eval} further confirms that \texttt{LPP}'s results lie in a low-quality, medium-diversity region.

\paragraph{\texttt{Ours}.} The goal of our mix-and-match perturbations was to make it hard for the model to decouple the random vector from the deterministic hidden state information. The success of our approach is confirmed by Fig.~\ref{fig:QD_plot}. Our model generates diverse motions, even after a long training time, and the quality of these motions is high. While this quality is slightly lower than that of \texttt{LHP} and \texttt{RHP} when looking at our classifier-based metric, it is rated higher by humans, as can be verified from Fig.~\ref{fig:human_eval}. We believe that this discrepancy is related to the binary classifier memorizing the ground-truth motions and thus not generalizing to the large diversity of motions generated by our model. As such, human evaluation still nicely complements our less expensive automatic one. Altogether, these results confirm the ability of our approach to generate highly diverse yet realistic motions. In Fig.~\ref{fig:qaulitative}, we further evidence this qualitatively by providing samples obtained by our approach for four different input sequences, as well as  samples from the baselines. Note that, for each input sequence, we produce large, yet natural variations of future poses.

Note that our approach depends on the parameter $\alpha$, which defines the amount of randomness used in our mix-and-match perturbations. In Fig.~\ref{fig:tradeoff}, we report the quality and diversity of our results when varying $\alpha\in\{0.1, \dots, 0.9\}$. Note that these plots show a trade-off between quality and diversity. This is to be expected since, by aiming to increase diversity, the resulting motion will become unrealistic. Nevertheless, our results can be seen to be highly diverse and of high quality for a wide range of values, i.e., by setting $\alpha\in[0.3,0.7]$. Note that, $\alpha$ is the only model-related hyper-parameter of the Mix-and-Match. The quality and diversity metrics are monotonic functions of $\alpha$, thus, one can choose a proper $\alpha$ given a task. Note that, using $\alpha=0.2$, our method still achieves a SoTA diversity of 2.25 with a higher quality of 45.0\%, however, for the sake of fair comparison, we use the default value of $\alpha=0.5$.

\subsection{Comparison with the State of the Art}
\label{sec:sota}
We now compare the results of our approach with those obtained using the state-of-the-art deterministic motion prediction methods~\cite{martinez2017human,jain2016structural,gui2018few,pavllo2018quaternet,pavllo2019modeling,fragkiadaki2015recurrent,gui2018adversarial} for long-term motion prediction, i.e., up to 1000ms. For this experiment, following previous work~\cite{fragkiadaki2015recurrent,pavllo2019modeling,pavllo2018quaternet,martinez2017human,gui2018adversarial}, we model velocity instead of pose, and do the same for the stochastic baselines. This is achieved by adding a residual connection to the motion decoder.
We then report the standard metric, i.e., the Euclidean distance between the predicted and ground-truth Euler angles. To evaluate this metric for our method which generates multiple, diverse predictions, we make use of the best sample among the $K$ generated ones with $K=50$ for the stochastic baselines and for our approach (i.e., the S-MSE metric~\cite{yan2018mt}).
In other words, we aim to show that, among the $K$ motions we generate, at least one is close to the ground truth. As can be seen in the top portion of Fig.~\ref{fig:long_term}, our approach yields errors comparable to the best-performing baselines, despite their use of more complex architectures and strong losses, such as the adversarial loss used in~\cite{gui2018adversarial}. Note that, unlike some of the baselines~\cite{gui2018adversarial,jain2016structural,martinez2017human}, our model requires knowing the action class during neither training nor inference.

In the bottom portion of Fig.~\ref{fig:long_term}, we further compare our best estimate with that of the other stochastic baselines, \texttt{LPH}, \texttt{RHP}, and \texttt{LPP}, using the best of $K$ motions in all cases. Note that, by providing better diversity, our approach outperforms these baselines.

\section{Evaluating the Effect of $K$}
In the main paper, we used $K=50$ to compare our approach with the state-of-the-art deterministic and stochastic baselines. Here, we provide an ablation study on the effect of $K$. To this end, we provide results when using $K=1$ to $K=500$. In Fig.~\ref{fig:k_effect}, we plot the results with $K=50$ as bold black lines, and the shaded area covers the results obtained with $K=1$ to $K=500$. While smaller values of $K$ yield large errors, the difference between $K=50$ and $K=500$ is very small (barely visible in most cases).

 \section{Conclusion}
In this paper, we have proposed an effective way of perturbing the hidden state of an RNN such that it becomes capable of learning the multiple modes of human motions. Our evaluation of quality and diversity, based on both new quantitative metrics and human judgment, have evidenced that our approach outperforms the state-of-the-art stochastic methods.
Generating diverse plausible motions given a short sequence of observations has many applications, especially when the motions are generated in an action-agnostic manner. For instance, our model can be used for human action forecasting~\cite{aliakbarian2016deep,rodriguez2018action,shi2018action,Aliakbarian_2017_ICCV,aliakbarian2018viena}, where one seeks to anticipate the action as early as possible.
It can also be employed for motion inpainting, where, given partial observations, one aims to generate multiple in-between solutions.
In the future, we will therefore investigate the use of our approach in such applications.

{\small
\bibliographystyle{ieee}

\section{Curriculum Learning of Perturbation}
As discussed in the main paper, our approach benefits from curriculum perturbation of hidden states. The parameter $\alpha$ determines the trade-off between the level of determinism and the quality of motions. In Algorithm 1, we provide the pseudo-code for our Mix-and-Match index sampling.

\begin{algorithm}[H]
\scriptsize
\SetAlgoLined
\textbf{Inputs:} current epoch \texttt{e}, perturbation ratio $\alpha$, sampling step \texttt{step}

\KwResult{Sampled indices for perturbation at epoch \texttt{e}, \texttt{indices}}

$s$ = int($\frac{\texttt{e}}{\texttt{step}}$)

$th$ = $min\{\frac{\alpha L}{2}, \frac{L-\alpha L}{2}\}$

\texttt{indices} = []

\If{$s<th$}
{
\texttt{indices}.append(sample $\alpha L - s$ indices from known part)

\texttt{indices}.append(sample $s$ indices from the rest)
}
\Else
{
\texttt{indices}.append(sample $\alpha L$ indices from hidden state)
}

\Return \texttt{indices}

 \caption{Curriculum Index Sampling}
 \label{alg:curriculum}
\end{algorithm}

\section{Human Evaluation of Quality}
As discussed in the main paper, when evaluating quality with a classifier, the classifier may \emph{memorize} the training set and thus become biased toward real data. As such, this would penalize the motions generated by our model, which are diverse enough to contain motions that were not seen during training. Therefore, human evaluation nicely complements the classifier-based quality assessment. In Table~\ref{tab:huamn}, we provide statistics related to the classifier-based and human evaluation of the quality of motions.

\begin{table}[!t]
    \centering
    \begin{tabular}{l c c}
         Baseline & $\texttt{Q}_\texttt{Classifier}$ & $\texttt{Q}_\texttt{Human}$\\
         \hline
         \texttt{LHP}~\cite{yan2018mt} & 45.0 & 37.0\\
         \texttt{RHP}~\cite{barsoum2018hp} & 47.0 & 38.1\\
         \texttt{LPP}~\cite{walker2017pose} & 13.0 & 15.0\\
         \hline
         \texttt{Ours} & 42.0 & 40.9\\
    \end{tabular}
    \caption{Comparing $\texttt{Q}_\texttt{Classifier}$ and human evaluation of quality for our approach and the baselines. These statistics correspond to evaluation after the models are fully trained.}
    \label{tab:huamn}
\end{table}

\section{Quality and Diversity of Motions}
In Table~\ref{tab:QD}, we provide the quality and corresponding diversity measures for the motions generated by our approach as well as the stochastic baselines compared with.

\begin{table*}[h]
\scriptsize
    \centering
    \begin{tabular}{l c c c c c c c c c c}
    \hline
     & \multicolumn{10}{c}{Training Progress}\\
    \hline
    Baseline & 10\% & 20\% & 30\% & 40\% & 50\% & 60\% & 70\% & 80\% & 90\% & 100\%\\
    \hline
        \texttt{LHP}~\cite{yan2018mt} & 1.19/29.0 & 1.13/37.0 & 0.28/45.0 & 0.24/46.0 & 0.21/46.0 & 0.47/45.0 & 0.28/44.0 & 0.26/45.0 & 0.24/46.0 & 0.26/45.0\\
        \texttt{RHP}~\cite{barsoum2018hp} & 0.76/27.0 & 0.65/40.0 & 0.62/43.0 & 0.58/45.0 & 0.54/44.0 & 0.52/44.0 & 0.51/45.0 & 0.50/45.0 & 0.48/46.0 & 0.48/47.0 \\
        \texttt{LPP}~\cite{walker2017pose}& 2.21/03.0 & 1.35/11.0 & 1.23/12.0 & 1.28/12.0 & 1.43/14.0 & 1.58/11.0 & 1.69/12.0 & 1.74/14.0 & 1.75/14.0 & 1.70/13.0 \\
        \hline
        \texttt{Ours} &
         2.02/22.0 &
         2.34/27.0 &
         2.73/30.0 &
         2.95/34.0 &
         3.17/39.0 &
         3.26/40.0 &
         3.49/41.0 &
         3.63/42.0 &
         3.53/42.0 &
         3.52/42.0 \\
        \hline
    \end{tabular}
    \caption{Quality and diversity of our approach and the stochastic baselines. The results are provided in the form of $Diversity/Quality$ at different training epochs.}
    \label{tab:QD}
\end{table*}

\section{Evaluating the Effect of $K$}
In the main paper, we used $K=50$ to compare our approach with the state-of-the-art deterministic and stochastic baselines. Here, we provide an ablation study on the effect of $K$. To this end, we provide results when using $K=1$ to $K=500$. In Fig.~\ref{fig:k_effect}, we plot the results with $K=50$ as bold black lines, and the shaded area covers the results obtained with $K=1$ to $K=500$. While smaller values of $K$ yield large errors, the difference between $K=50$ and $K=500$ is very small (barely visible in most cases).
\begin{figure*}[!t]
    \centering
    \begin{tabular}{c@{ }c@{ }c@{ }c}

    \input{alpha_walking.tex}  &
    \input{alpha_eating.tex}  &
    \input{alpha_smoking.tex}  &
    \input{alpha_discussion.tex} \\
    \end{tabular}
    \caption{Effect of $K$ in the MAE of our actions of the Human3.6M dataset. Note, the bold black one is the best of $K=50$ motions, and the shaded area indicates the region between best of $K=1$ and $K=500$.}
    \label{fig:k_effect}
\end{figure*}

\section{Long-term Motion Prediction}

In the main paper, we provide the plots of human motion predictions up to 1000ms for our approach and the recent deterministic and stochastic baselines on four actions of the Human3.6M dataset.  Note that, as discussed in the main paper, for our approach and the stochastic baselines, we report the best of $K=50$ motions. The numbers associated with this experiment are provided in Table~\ref{tab:det}.

\begin{table*}[!t]
    \centering
    \scriptsize
    \begin{tabular}{c c}
        \begin{tabular}{l| c c c c c c}
    \hline
    \multicolumn{7}{c}{Walking}\\
    \hline
    Method & 80ms  & 160ms & 320ms & 400ms & 560ms & 1000ms \\
    \hline
    Zero Velocity &
    0.39 & 0.86 & 0.99 & 1.15 & 1.35 & 1.32 \\
    AGED~\cite{gui2018adversarial} &
    0.22 & 0.36 & 0.55 & 0.67 & 0.78 & 0.91 \\
    LSTM-3LR~\cite{fragkiadaki2015recurrent} &
    1.18 & 1.50 & 1.67 & 1.76 & 1.81 & 2.20 \\
    SRNN~\cite{jain2016structural} &
    1.08 & 1.34 & 1.60 & 1.80 & 1.90 & 2.13 \\
    DAE-LSTM~\cite{ghosh2017learning} &
    1.00 & 1.11 & 1.39 & 1.48 & 1.55 & 1.39 \\
    GRU~\cite{martinez2017human} &
    0.28 & 0.49 & 0.72 & 0.81 & 0.93 & 1.03 \\
        \hline
        \hline
    \texttt{LHP}~\cite{yan2018mt} &
    0.73 & 0.79 & 0.90 & 0.93 & 0.95 & 1.05 \\
    \texttt{RHP}~\cite{barsoum2018hp} &
    0.61 & 0.62 & 0.71 & 0.79 & 0.83 & 1.07 \\
    \texttt{LPP}~\cite{walker2017pose} &
    0.56 & 0.66 & 0.98 & 1.05 & 1.28 & 1.60 \\
    \hline
    Ours &
    0.33 & 0.48 & 0.56 & 0.58 & 0.64 & 0.68 \\
    \hline
    \end{tabular} &
    \begin{tabular}{l| c c c c c c}
    \hline
    \multicolumn{7}{c}{Eating}\\
    \hline
    Method & 80ms  & 160ms & 320ms & 400ms & 560ms & 1000ms \\
    \hline
    Zero Velocity &
    0.27 & 0.48 & 0.73 & 0.86 & 1.04 & 1.38 \\
    AGED~\cite{gui2018adversarial} &
    0.17 & 1.28 & 0.51 & 0.64 & 0.86 & 0.93 \\
    LSTM-3LR~\cite{fragkiadaki2015recurrent} &
    1.36 & 1.79 & 2.29 & 2.42 & 2.49 & 2.82 \\
    SRNN~\cite{jain2016structural} &
    1.35 & 1.71 & 2.12 & 2.21 & 2.28 & 2.58 \\
    DAE-LSTM~\cite{ghosh2017learning} &
    1.31 & 1.49 & 1.86 & 1.89 & 1.76 & 2.01 \\
    GRU~\cite{martinez2017human} &
    0.23 & 0.39 & 0.62 & 0.76 & 0.95 & 1.08 \\
    \hline
    \hline
    \texttt{LHP}~\cite{yan2018mt} &
    0.68 & 0.74 & 0.95 & 1.00 & 1.03 & 1.38 \\
    \texttt{RHP}~\cite{barsoum2018hp} &
    0.53 & 0.67 & 0.79 & 0.88 & 0.97 & 1.12 \\
    \texttt{LPP}~\cite{walker2017pose} &
    0.44 & 0.60 & 0.71 & 0.84 & 1.05 & 1.54 \\
    \hline
    Ours &
    0.23 & 0.34 & 0.41 & 0.50 & 0.61 & 0.91 \\
    \hline
    \end{tabular}
    \\
    \begin{tabular}{l| c c c c c c}
    \hline
    \multicolumn{7}{c}{Smoking}\\
    \hline
    Method & 80ms  & 160ms & 320ms & 400ms & 560ms & 1000ms \\
    \hline
    Zero Velocity &
    0.26 & 0.48 & 0.97 & 0.95 & 1.02 & 1.69 \\
    AGED~\cite{gui2018adversarial} &
    0.27 & 0.43 & 0.82 & 0.84 & 1.06 & 1.21 \\
    LSTM-3LR~\cite{fragkiadaki2015recurrent} &
    2.05 & 2.34 & 3.10 & 3.18 & 3.24 & 3.42 \\
    SRNN~\cite{jain2016structural} &
    1.90 & 2.30 & 2.90 & 3.10 & 3.21 & 3.23 \\
    DAE-LSTM~\cite{ghosh2017learning} &
    0.92 & 1.03 & 1.15 & 1.25 & 1.38 & 1.77 \\
    GRU~\cite{martinez2017human} &
    0.33 & 0.61 & 1.05 & 1.15 & 1.25 & 1.50 \\
    \hline
    \hline
    \texttt{LHP}~\cite{yan2018mt} &
    1.00 & 1.14 & 1.43 & 1.44 & 1.68 & 1.99 \\
    \texttt{RHP}~\cite{barsoum2018hp} &
    0.64 & 0.78 & 1.05 & 1.12 & 1.64 & 1.84 \\
    \texttt{LPP}~\cite{walker2017pose} &
    0.59 & 0.83 & 1.25 & 1.36 & 1.67 & 2.03 \\
    \hline
    Ours &
    0.23 & 0.42 & 0.79 & 0.77 & 0.82 & 1.25 \\
\hline
    \end{tabular}     &
    \begin{tabular}{l| c c c c c c}
    \hline
    \multicolumn{7}{c}{Discussion}\\
    \hline
   Method & 80ms  & 160ms & 320ms & 400ms & 560ms & 1000ms \\
    \hline
    Zero Velocity &
    0.31 & 0.67 & 0.94 & 1.04 & 1.41 & 1.96 \\
    AGED~\cite{gui2018adversarial} &
    0.27 & 0.56 & 0.76 & 0.83 & 1.25 & 1.30 \\
    LSTM-3LR~\cite{fragkiadaki2015recurrent} &
    2.25 & 2.33 & 2.45 & 2.46 & 2.48 & 2.93 \\
    SRNN~\cite{jain2016structural} &
    1.67 & 2.03 & 2.20 & 2.31 & 2.39 & 2.43 \\
    DAE-LSTM~\cite{ghosh2017learning} &
    1.11 & 1.20 & 1.38 & 1.42 & 1.53 & 1.73 \\
    GRU~\cite{martinez2017human} &
    0.31 & 0.68 & 1.01 & 1.09 & 1.43 & 1.69 \\
    \hline
    \hline
    \texttt{LHP}~\cite{yan2018mt} &
    0.80 & 1.01 & 1.22 & 1.35 & 1.56 & 1.69 \\
    \texttt{RHP}~\cite{barsoum2018hp} &
    0.79 & 1.00 & 1.12 & 1.29 & 1.43 & 1.71 \\
    \texttt{LPP}~\cite{walker2017pose} &
    0.73 & 1.10 & 1.33 & 1.34 & 1.45 & 1.85 \\
    \hline
    Ours &
    0.25 & 0.60 & 0.83 & 0.89 & 1.12 & 1.30 \\
\hline
    \end{tabular}

    \end{tabular}
    \caption{Comparison against the state of the art for four actions of the Human3.6M dataset.}
    \label{tab:det}
\end{table*}

{\small
\bibliographystyle{ieee}


\section{Introduction}

Windowing computations are imperative for visualizing anatomies with a suitable contrast. This involves the determination of two preset values called window width (WW) and window level (WL) and using them to transform the intensity distribution of the image. However, the presence of background noise in MR images biases the calculation of default WW/WL values since there has to be a trade-off between enhancing contrast of foreground/anatomy of interest vs suppressing background/ outside the anatomy of interest. A sub-optimal WW/WL would lead to poor visualization making human intervention necessary more often than needed to make the image readable. Also, as mentioned earlier, it is challenging to manually adjust WW/WL to achieve necessary contrast in anatomy and hiding the background pixels at the same time. All this may lead to visual fatigue for the viewer and consequently a delayed diagnosis.

This paper proposes a deep learning (DL) based method that improves the WW/WL for original as well as derived images (e.g., functional maps) by intelligently eliminating the unwanted background pixels. This also provides better control for user defined windowing in post processing as background pixels are already eliminated.

\section{Literature review}
Most of the literature for simple adaptive windowing focuses on suppressing background/ denoising via simple global thresholding \cite{felmlee1999automatic}. A slightly more advanced method uses thresholding followed by morphological operations to extract RoI \cite{belykh2007method}. Also, thresholding creates holes within anatomy leading to loss of information and hence sub-optimal windowing. There are other methods which uses segmentation as an alternative. For example, the adaptive windowing described in \cite{kaushikadaptive} is specific to contrast enhancement of the structure of interest within the anatomy and uses a level-set based segmentation algorithm to extract the RoI. A review of traditional segmentation methods\cite{sharma2010automated} shows that there is no universal algorithm for segmentation of every medical image and most of these methods require hand crafted tuning. The ideal segmentation technique for an image depends on the application and the body part that is imaged.

Alternatively, a deep learning (DL) based approach has the advantage that it is a generic method which performs equally well for all applications, anatomies and protocols. The existing methods using the same \cite{lai2005adaptive} are complex and intensive and needs user defined WW/WL values for the whole training set.

The proposed method overcomes the inherent limitations of the thresholding based segmentation and WW/WL calculation by using DL based background suppression and the same also provides additional benefit of enhanced control in manual adjustment of WW/WL which is not the case in \cite{lai2005adaptive} since it just computes the WW/WL through a neural network.

\section{Methodology}
The steps used to achieve the revised WW/WL are as follows.
\begin{itemize}[noitemsep]
\itemsep0em
\item[1]Obtain background suppressed image.
\item[2]Use the foreground pixels to compute the WW/WL.
\item[3]Display the segmented image with the newly computed WW/WL
\end{itemize}

In step 1 above, our method uses U-net for implementing background suppression.It was trained with 2700 images covering brain and abdomen and the corresponding masks representing the anatomy.
We propose to use background suppressed image with auto computed WW/WL for all the further reading and processing. This will ensure that not only the image (original/ post processed) is displayed with improved window settings but also provide better control on dynamically adjusting the WW/WL afterwards.

\section{Results}
The accuracy of auto-windowing computation depends on the accuracy of segmentation, our evaluation strategy was twofold- evaluating the accuracy of U-net for segmentation followed by image review from clinical application specialists.
\begin{itemize}
\item[1]U-net was used to perform segmentation on 1130 test images covering brain and abdomen and an average DICE score of 0.94 (min = 0.713; max = 0.982) was achieved.
\item[2]The background suppressed images were presented with a better WW/WL than the original images(refer \figureref{example})
\item[3]It was also confirmed that the background suppressed image provided better control while dynamically adjusting WW/WL.
\end{itemize}

\begin{figure}[htbp]

\floatconts
  {fig:example}
  {\caption{WW/WL before and after background suppression} \label{example}}
  {\includegraphics[width=0.7\linewidth]{figure.JPG}}
\end{figure}

\section{Conclusion}

The proposed solution addresses the problem of poor image quality and suboptimal WW/WL due to presence of background without compromising the anatomical details. The DL based segmentation demonstrates effective background suppression with high DICE scores. The WW/WL of both original as well as derived images can be significantly improved by preprocessing using proposed method.

\midlacknowledgments{We would like to thank Ananthakrishna Madhyastha, MR Viz and Recon Manager and Venkatesh Seshachar, MR Software Engg Director from GE Healthcare for all the support and encouragement that helped make this paper.}


\section{Methodology}

Customers experience in fashion e-commerce consists in selecting an article in a desired size, trying the article, forming an opinion on its size and returning or keeping it. To simulate the customers behavior towards sizing, we model the joint probability of a customer to pick a size and the resulting return status. Return status is described by three possible events: the customer keeps the article, the customers returns the article because it's too small and the customer returns the article because it is too big. We ignore the cases where the customer returns the article for any other reason.

\subsection{Notation}

Let us denote $\mathcal{C}$ the set of customers and $\mathcal{A}$ the set of articles. The size $\mathcal{S}_i$ is a continuous random variable. The variable $\mathcal{R}$ indicates the return status described above. Orders $\mathcal{O}$ are defined by a customer, an article, the purchased size and the return status. Both approaches introduced in the paper model the joint probability $p(\mathcal{S}_o, \mathcal{R}_o \mid \mathcal{C}_o, \mathcal{A}_o)$ as detailed in the following.
\subsection{Baseline Model}
\label{sec:baseline}
The baseline model makes a simplifying assumption that the size the customer chose and the return status are two independent events. Thus, the joint probability is defined as the product of the probability over sizes and the probability of return status:

\begin{equation}
p(\mathcal{S}, \mathcal{R} \mid \mathcal{C}, \mathcal{A}) = p(\mathcal{S}  \mid \mathcal{C}, \mathcal{A}) p(\mathcal{R} \mid\mathcal{C}, \mathcal{A})
\end{equation}

\paragraph{Probability over sizes}

We assume that multiple persons can use a single account. The probability distribution over sizes is obtained by Gaussian Kernel Density Estimation. To avoid getting degenerate distributions when customers always purchase the same size, we set a lower limit on the variances. Let denote $\mathcal{O}_j$ the set of $n_j$ orders and $S_j = \{s_i, i=1..n_j\}$ the set of $n_j$ sizes purchased by the customer $c_j$. The related probability density function is defined as:

\begin{equation}
p(s \mid c_j) = \dfrac{1}{n_jh_j} \displaystyle\sum_{i=1}^{n_j} \phi\left(\dfrac{s-s_i}{h_j}\right)
\end{equation}

where $\phi$ is the normal density function and $h_j$ is the bandwidth parameter for that specific customer $c_j$. The latter is obtained by minimizing the mean integrated squared error \cite{Silverman1986}.

\paragraph{Probability over return status}

Customers have different return behaviors. However, impact of customers on returns is assumed negligible compared to potential sizing issues of the article. The probability of each return status is  consequently marginalized over customers: $p(\mathcal{R} \mid\mathcal{C}, \mathcal{A}) = p(\mathcal{R} \mid \mathcal{A})$. %
The probability over return status is the empirical distribution over the three possible events: article is kept, too big or too small. The case of one of the events being not observed in the training data may lead to a null probability in validation. To avoid that problem, we add one to the counts of each event. For an article $a_i$, sold $n_i$ times, the probability of a return event $r$, observed $n_{i,r}$ times in the data, is defined as:
\begin{equation}
p(r \mid a_i) = \dfrac{n_{i,r} + 1}{n_i + 3}
\end{equation}

Though, this method seems inelegant, it has a Bayesian grounding. Indeed, this is equivalent to taking the maximum a posteriori of a categorical distribution with a Dirichlet conjugate prior, which concentration parameter is equal to one, i.e the uniform distribution.
In case of a cold start, i.e if a customer (resp. an article) is new, the marginal distribution over sizes of all the customers (resp. over return status of all articles) is used.

\subsection{Hierarchical Bayesian Model}
\label{sec:BHM}
The baseline described above has a risk of specious parameter estimation and overfitting. Bayesian approaches conversely aim at providing a probability of the estimated parameters given a set of observed data, supporting the decision process when offering a recommendation to the customer. Therefore, using a Bayesian approach, we aim at modeling the joint probability of a size to be purchased and a return status to be observed without the simplifying hypothesis from the base approach. For each pair of customer and article, orders $\mathcal{O}$ are drawn following a categorical distribution.

\begin{equation}
\mathcal{O} \sim Cat(p(\mathcal{S}, \mathcal{R} \mid \mathcal{C}, \mathcal{A}))
\end{equation}

Contrary to the baseline, both $\mathcal{S}$ and $\mathcal{R}$ are not assumed independent, instead the joint probability is factorized as:

\begin{equation}
p(\mathcal{S}, \mathcal{R} \mid \mathcal{C}, \mathcal{A}) = p(\mathcal{S} \mid \mathcal{R}, \mathcal{C}, \mathcal{A}) \times p(\mathcal{R} \mid \mathcal{C}, \mathcal{A})
\label{eq:factorization}
\end{equation}

It is worth noting that methods described in \cite{Sembium2018,abdulla2017} aim at modeling $p(\mathcal{R} \mid \mathcal{S}, \mathcal{C}, \mathcal{A})$. Doing so requires discretizing the continuous variable $\mathcal{S}$ leading to an increase in the number of parameters to be inferred and making the model more susceptible to the sparsity of the data. The factorization we chose in Equation \ref{eq:factorization} allows us to model $p(\mathcal{S} \mid \mathcal{R}, \mathcal{C}, \mathcal{A})$ as a continuous distribution. This enables us to learn a smaller set of parameters specifying the distribution over all sizes, which helps to alleviate part of the sparsity problem.

\paragraph{Probability over return status}

For the same reasons as explained in the Section \ref{sec:baseline}, the probability of return status is marginalized over customers: $p(\mathcal{R} \mid \mathcal{C}, \mathcal{A}) = p(\mathcal{R} \mid \mathcal{A})$. Returns are assumed independent from one another, allowing us to model them using a categorical distribution. As the number of purchases might be low for some articles, a Dirichlet prior is used. The concentration parameter of the prior is based on the counts at the brand level and at the category level (e.g. dresses, t-shirts, sneakers, etc.). Let $n_K$, $n_S$, and $n_B$ denote the counts of kept articles, returned articles for being too small and too big at the article level respectively. In a similar way, $m_K$, $m_S$ and $m_B$ indicates the counts at the brand level and $m'_K$, $m'_S$ and $m’_B$, at the category level.

\begin{equation}
p(\mathcal{R} \mid \mathcal{A}) \sim Dirichlet(\alpha)
\end{equation}
with $\alpha = w \cdot [m_K, m_S, m_B] + w’ \cdot [m'_K, m'_S, m'_B]$. The weights $w$ and $w'$ are learned under the assumption that they follow Beta distribution with a low first shape parameter and the second shape parameter equal to 1, in order to favor low weight values.
\begin{equation}
w \sim Beta(0.5, 1) \mbox{ and } w' \sim Beta(0.1, 1)
\end{equation}
Since the Dirichlet prior is the conjugate of the categorical distribution, the posterior probability can be analytically computed, easing the inference of the parameters of the model.
\begin{equation}
P(\mathcal{R}=r \mid \mathcal{A}, \mathcal{O} ) = \dfrac{n_r + \alpha_r}{\displaystyle\sum_{\substack{i \in \\ \{K,S,B\}}} n_i + \displaystyle\sum_{\substack{i \in \\ \{K,S,B\}}}\alpha_i}
\end{equation}

\paragraph{Probability over sizes conditionally on the return status}

For a given customer and article, the probability distribution of the customer buying a size is a mixture of Gaussians. Since the number of users of an account is unknown, we decide to use an infinite mixture model. However, we assume that the number of distinct persons using a single account is low. That's why we opt for a Dirichlet process with a truncation level fixed to four. In order to ease the inference, we use a truncated stick-breaking process \cite{sethuraman1994}.
\begin{equation}
\begin{split}
&p(s \mid r, c, a) = \displaystyle\sum_{i=1}^{4} \pi_i \phi(s \mid \mu_i, \sigma_i^{2}) \\
&\pi_i = b_i \displaystyle\prod_{j=1}^{i-1} (1-b_j) \\
&b_i \sim Beta(1, \alpha) \mbox{ for }  i=1..3 \mbox{ and } b_4 = 1
\end{split}
\label{eq:mix_gaussian}
\end{equation}
The parameter $\pi$ is the mixing proportion, that can be interpreted as the probability of person $i$ using the account of the customer $c$. The shape parameter $\alpha$ of the $Beta$ distribution is acting as a concentration parameter of the Dirichlet process: the case $\alpha = 1$ is equivalent to the uniform distribution, thus favoring the scenario of multiple persons sharing a single account; conversely, the case $\alpha \rightarrow 0$ produces high density around 1 which favors the scenario of a single person placing all the orders. In the context of size recommendation, $\alpha$ is fixed at $0.5$.

The function $\phi$ in equation \ref{eq:mix_gaussian} is the normal probability density function over sizes for the person $i$ using the account of the customer $c$, buying an article $a$, resulting in a return status $r$. The parameter $\mu$ is a combination of three parameters $\mu = \mu_C + \mu_A + \eta_R$:
\begin{itemize}
\item the average size of a person $\mu_C$,\begin{equation}
\mu_C \sim \mathcal{N}(\mu_0, \sigma_0^2)
\end{equation}
where hyperparameters $\mu_0$ and $\sigma_0^2$ depend on the category, the gender of the article and the size system;
\item the average offset of the article $\mu_A$,
\begin{equation}
\mu_A \sim \mathcal{N}(0, 1)
\end{equation}
where assumption is made that most articles have an accurate size, i.e. an offset of zero;
\item a shifting parameter $\eta = \{\eta_K, \eta_S, \eta_B\}$ for each return status: resp. article is kept, returned too small and returned too big,
\begin{equation}
\eta_K = 0  \mbox{ ; } \eta_S \sim \mathcal{N}(-1, 1) \mbox{ ; } \eta_B \sim \mathcal{N}(1, 1)
\end{equation}
shifting parameter is fixed at 0 for the case the article is kept, while it's sampled using a Gaussian distribution centered to 1 (resp. -1) when the customer has returned the article because it is too big (resp. too small).
\end{itemize}
We assume the parameter $\sigma$ in equation \ref{eq:mix_gaussian} depends on the customer only. It is sampled following an Inverse Gamma distribution, with the shape parameter $\gamma_1$ and the scale parameter $\gamma_2$.
\begin{equation}
\sigma_C ^2 \sim \Gamma^{-1}(\gamma_1,\gamma_2)
\end{equation}
We fix the parameters of the distribution to $\gamma_1 = 1$ and $\gamma_2 = 2$, so that the mode of the Inverse Gamma distribution is equal to 1. Figure \ref{fig:bayesian_graphe} represents the graphical model of the approach.
\paragraph{Inference}
Monte-Carlo Markov Chains are popular sampling methods for Bayesian inference. But those approaches are often slow to converge. Variational inference methods run faster than sampling based methods but also introduce an approximation bias that may lead to a bad estimation of the parameters. However, those approaches have demonstrated reasonable performances on Dirichlet processes \cite{Blei2006} and are well suited for problems involving large amount of data. The inference is consequently done using mean-field approximation.

\input{graphical_model.tex}

\subsection{Providing a Size Recommendation}

The set of sizes of an article is a finite set and as a consequence the probability density function needs to be discretized. For a customer $c$ and an article $a$ with a set of $k$ sizes $\mathcal{S} = \{s_i, i=1..k\}$, the probability over sizes is discretized as follows:
\begin{equation}
p(s, r) = p(r) \dfrac{\displaystyle\int_{s-\frac{1}{2}\epsilon}^{s+\frac{1}{2}\epsilon} f(x) dx}{\displaystyle\int_{s_1-\frac{1}{2}\epsilon}^{s_k+\frac{1}{2}\epsilon} f(x) dx}
\end{equation}
where $f$ is the probability over sizes, marginal for the baseline and conditionally to the return status for the Bayesian model ; $\epsilon$ is equal to the step between two sizes.
To provide a size recommendation, we choose the size having the highest probability to be kept by the customer.

\section{Conclusion}
A hierarchical Bayesian approach was proposed to tackle the challenging problem of size recommendation in e-commerce fashion. The size purchased by a customer and its possible return events were jointly modeled thanks to a Bayesian approach. Experimental results were presented on real (anonymized) data from millions of customers along with a detailed discussion and comparison with a baseline approach with simplified hypothesis. It was shown that the Bayesian approach outperforms the baseline approach while providing better theoretical framework for gaining deeper understanding of the predictions. Future work consists in exploring different approaches to learn the joint probability, making use of additional article related data from fashion industry, along with deeper dives in the segmentation of customers and articles.
\section{Experiments}
For this experiment, anonymized purchase data is collected for adult shoes. The data consists of 14.5 million purchases, 3 million distinct customers and 73,000 distinct articles. As the data has a strong temporal component, cross validation is performed under the following conditions \cite{arlot2010}: 1. validation data occurs later than training data, 2. a period of three weeks - corresponding to the time to collect most returns from customers - is ignored between train set and validation set and 3. validation sets must not overlap. Both models are trained and cross-validated on the same data. We reported the average logarithm of the joint probability $p(s,r)$ over all observations in the validation set in Table \ref{tab:likelihood}. Higher numbers show better performances of the model.

\begin{table}[!h]
  \caption{Average log joint probability}
  \vspace{-0.3cm}
  \begin{tabular}{ | l | l | l |}
    \hline
     & Baseline & Bayesian \\ \hline
    incl. unknown customers & -2.85 (± 0.15) & -2.35 (± 0.11)\\ \hline
    excl. unknown customers & -3.32 (± 0.26)& -1.83 (± 0.19)\\
    \hline
  \end{tabular}
\label{tab:likelihood}
\end{table}

Table \ref{tab:likelihood} compares likelihood for both approaches, including and excluding customers not observed in the training data. Bayesian approach shows better results in both cases, where results are the best when excluding unknown customers. Conversely, the baseline performs better when unknown customers are included. This is mainly due to the fact that the baseline overfits by putting high probability density on the events observed in the training set, and very low density on unseen events.

In the context of size recommendation, two indicators play a key role in the decision process: coverage and accuracy. The coverage is the percentage of purchases for which the algorithm is confident making a decision. The accuracy is the number of correctly predicted sizes and return status, over all predictions. Figure \ref{fig:cov_vs_acc} (Top) shows the accuracy versus the coverage for several values of a sliding threshold on the joint probability for both models. On Figure \ref{fig:cov_vs_acc} (Bottom), the results are presented for the Bayesian model where we also include a threshold on the posterior probability of parameters.

\begin{figure}[!ht]
\centering
\subfloat{\includegraphics[width=0.47\textwidth]{coverage_vs_accuracy.pdf}}
\quad
\subfloat{\includegraphics[width=0.47\textwidth]{coverage_vs_accuracy_post.pdf}}
\caption{Accuracy versus coverage. (Top) Baseline and Bayesian model. (Bottom) Bayesian model with and without threshold on the posterior probability of the parameters.}
\label{fig:cov_vs_acc}
\end{figure}
Figure \ref{fig:cov_vs_acc} (Top) shows that accuracy decreases as coverage increases, when changing the threshold on the joint probability. The performances of the Baseline and the Bayesian model are similar, when the decision is based on the value of the joint probability. In Figure \ref{fig:cov_vs_acc} (Bottom), by putting a threshold on the posterior probability of parameters, we prevent the model from recommending article sizes to customers, where the parameters are poorly estimated. Performances are slightly better, however only 13\% of the purchases can be covered.

It is worth mentioning that both approaches need to filter out a lot of purchases before starting to show reasonable accuracy levels. From the computational complexity point of view, inference of the Bayesian approach is more costly compared to the baseline model. Results are encouraging and demonstrate the complexity of the size recommender topic, motivating a deeper research work in the field.

\section{Introduction}

Fashion is a way to express identity, moods and opinions. Customers also tend to use fashion to either emphasize certain parts of their body or hide others. In that context, size and fit have been shown to be among factors influencing the most the overall satisfaction \cite{Pisut2017}. Online customers have to buy before trying their clothes on. The sensory feedback phase about how the article fits via touch and visual cues is then delayed. Because of these uncertainties, a lot of consumers are still reluctant to engage in the purchase process.

To make matters worse, fashion articles including shoes and apparel have important sizing variations primarily due to: 1. different definitions of respective sizes from brands: the size systems used for specific categories are limited (e.g. S, M, L, etc.), however the sizes themselves represent different physical measurements from one brand to another; 2. different ways of converting a local size system to another: in Europe, garment sizes are not standardized and brands don't always use the same conversion logic.

A way to circumvent the confusion created by these variations is to use size tables which map physical body measurements to the article size system, requiring customers to have accurate measurements of their body. However, size tables themselves might suffer from a large variance, up to one inch within a single size. These differences stem from either different datasets used for size tables (e.g. German vs. UK population) or are due to vanity sizing, i.e deliberate size inconsistencies in brands targeting a specific focus group based on age, sportiness, etc. which represent major influences on the body measurements \cite{ujevic2005, shin2007,faust2014}.
The combination of the above factors leaves the customers alone to face a highly challenging problem of determining the right size and fit during their purchase journey. In recent years, there has been a lot of interest in building recommender systems in fashion e-commerce with major focus on modeling style preferences based on customers past interactions, taste and affinities \cite{hu2015,arora2016,bracher2016}. However, few research work have been conducted to tackle the size recommendation problem.

The recommendation of size and fit has been recently studied in \cite{abdulla2017} where sparsity in purchased data is mentioned as a major issue, especially considering that articles have a limited stock. To minimize that problem, the authors propose to represent articles as a combination of brand, usage, size, and fit. A neural network is then trained to learn a latent vector describing each article defined as the combination of features mentioned before. Customer vector representation is obtained by aggregating over purchased articles and, finally, a gradient boosted classifier predicts the fit of an article to a customer.

Following a different approach, the authors of \cite{sembium2017} propose a solution for determining if an article of a certain size would be fit, large, or small for a certain customer, using the purchase history. This is achieved by iteratively deducing the true sizes for customers and products, fitting a linear function based on the difference in sizes, and performing ordinal regression on the output of the function to get the loss. Extra features are simply included by addition to the linear function. To handle multiple persons behind a single account, hierarchical clustering is performed on each customer account before doing the above. An extension of that work has been very recently published proposing a Bayesian approach on a similar model \cite{Sembium2018}. Instead of learning the parameters in an iterative process, the updates are done with mean-field variational inference with Polya-Gamma augmentation. This method therefore naturally benefits from the nice advantages of Bayesian modeling - the uncertainty outputs, and the use of priors.
In this paper, we present two approaches: a baseline algorithm, launched on shoes in 2016 and on garments in 2017, which consists in inferring the size that a customer intends to buy and, independently, the article's sizing characteristics ; and a hierarchical Bayesian approach which aims at jointly modeling the purchases of one or multiple sizes of an article along with their possible return events: 1. no return (article is kept), 2. returned too small 3. returned too big.

In the size recommendation context, data sparsity is severe since it affects both articles and customers. To tackle this, we propose two design choices: a) building a hierarchy on top of parameters, exploiting prior knowledge on articles and customers; b) virtually treating the size as a continuous variable in the training phase.