\section{Installation}

\texttt{Lagrange2D} requires Mathematica 10.0 or higher. The source code is available on GitHub: \texttt{https://github.com/williamgilpin/lagrange2d}. After downloading the repository, the functions contained in the package can be accessed in a local notebook by running the following (with an appropriate absolute path substituted):

\begin{verbatim}
Import["/Users/me/path/to/folder/lagrange2d/lagrange2d.wl"]
\end{verbatim}

Alternatively, the file \texttt{lagrange2d.wl} can be manually moved into the working directory of the local notebook, in which case the following lines will import the package:

\begin{verbatim}
SetDirectory[NotebookDirectory[]]
<< lagrange2d.wl
\end{verbatim}

After successful installation, detailed and executable examples demonstrating the various functions in the package are available in the included \texttt{demos.nb} notebook.

\section{Included functions}

We summarize in Table \ref{funcs} the different functions included in Lagrange2D. Here, we highlight several functions and their underlying theoretical basis.

\begin{table}[htp]
\centering
\caption{Functions included in Lagrange2D. Full documentation and runnable examples are included in the notebook "demos.nb" on GitHub at \texttt{https://github.com/williamgilpin/lagrange2d/blob/master/demos.nb} }
\vspace{4 mm}
\label{funcs}
\begin{tabular} {|p{5cm}|p{6cm}|}
\hline
Function & Description\\
\hline
\texttt{advectPoints} & Numerically advect of a set of initial conditions\\
\hline
\texttt{findFTLEField}		& Find all finite-time Lyapunov exponents\\
\hline
\texttt{findMaxFTLEField}		& Find the maximum finite-time Lyapunov exponents\\
\hline
\texttt{findStretchlines}		& Find the vector fields corresponding to maximum and minimal stretching\\
\hline
\texttt{findKYDim}	& Find the Kaplan-Yorke fractal dimension\\
\hline
\texttt{flushingTimes}		& Find the flushing time\\
\hline
\texttt{pathPlot}		& Draw pathlines of fixed duration for a set of initial conditions\\
\hline
\texttt{animateFlow}		& Generate a movie file showing the temporal evolution of a set of initial conditions\\
\hline
\texttt{makeMesh} 	& Draw a uniform rectangular mesh\\
\hline
\end{tabular}
\end{table}

\subsection{Finite-time Lyapunov Exponents}

The Lyapunov exponents of a dynamical system quantify the degree to which trajectories originating at given points in the domain spread apart over time. While these formally may be calculated by computing trajectories of pairs of neighboring points on the attractor, recent work has shown that they are well-approximated by the finite-time Lyapunov exponent \cite{haller2015lagrangian,shadden2011lagrangian}, defined as
\begin{equation}
\lambda(\v r) \equiv \dfrac{1}{\tau}\max\left[\text{eig}\left(			\log\left(	\d{\phi_{0}^{\tau}(\v r)}{\v r}	\right)		\right)\right]
\end{equation}
where the flow map $\phi_{0}^{\tau}$ transforms a tracer particle at initial position $\v r$ to its final position $\v r'$  after a time $\tau$. In practice, this requires numerically integrating an ensemble of trajectories for a duration $\tau$, and then numerically interpolating the map $\phi_{0}^{\tau}$.

In the functions \texttt{findFTLEField} and \texttt{findMaxFTLEField}, a region of initial conditions is specified, and then advected forward in time for a specified duration. The displacement map $\phi_{0}^{\tau}$ is calculated using Mathematica's built-in \texttt{Interpolate} function.

\subsection{Kaplan-Yorke Fractal Dimension}

The Kaplan-Yorke fractal dimension quantifies the structure of a dynamical system's attractor \cite{kaplan1979functional}. Unlike the box-counting and information dimensions, the Kaplan-Yorke dimension has the advantage of being directly computable from particle trajectories, making it well-suited for the analysis of unsteady flows. For a general dynamical system, the Kaplan-Yorke dimension is defined as
\begin{equation}
D_{KY} \equiv j+{\frac  {\sum _{{i=1}}^{j}\lambda _{i}}{|\lambda _{{j+1}}|}},
\label{ky}
\end{equation}
where the Lyapunov exponents of the system, $\lambda_i$, have been ranked from largest to smallest, and $j$ denotes the first index that satisfies the relationship
\[
\sum _{i=1}^{j}\lambda _{i}\geqslant 0, \qquad \sum _{{i=1}}^{{j+1}}\lambda _{i}<0.
\]
In two dimensional, incompressible flows, $j = 1$. However, when using estimates of $\lambda$, such as those inferred from experimental velocity fields, or those resulting from differential methods like particle image velocimetry, strict incompressibility is not necessarily preserved, and so $j \in \{0,1,2\}$. In the function \texttt{findKYDim} in \texttt{Lagrange2D}, we use estimates of the Lyapunov exponents generated by the FTLE calculation in \eqref{ky}, in order to produce a timescale-dependent scalar field $\mathcal D_{KY}^t(x,y)$.

In fluid flows, the Kaplan-Yorke dimension quantifies the tendency of particles to cluster within subregions of the fluid's attractor: wherever the Kaplan-Yorke dimension equals the dimensionality of the dependent variables ($d=2$ here), then there is no particular tendency of trajectories to localize \cite{vilela2017map,tuval2004opening}. In contrast, in subregions where the Kaplan-Yorke dimension is much less than the dimension of the underlying dynamics, particles will tend to cluster and form filamentary structures. In periodic or statistically stationary flows, the spatially-resolved Kaplan-Yorke dimension can therefore be used to differentiate the stable and unstable manifolds of the flow.

\section{Example Usage}

The following example and images are available in executable form in the \texttt{demos.nb} notebook included in the installation directory. Plots produced by the examples below are shown in Figure \ref{fig_lagrange}, and full documentation of arguments and keyword arguments for each function can be accessed via the function's docstring. For example, to access the documentation for the function \texttt{pathPlot}, the command \texttt{?pathPlot} can be executed anywhere within an active notebook.

\begin{figure}
{
\centering
\includegraphics[width=\linewidth]{fig_lagrange.pdf}
\caption{
{\bf Example plots produced by Lagrange2D.} (A) Streamlines of the time-dependent "double gyre" flow, generated at times corresponding to times when the flow is maximally symmetric (blue) and maximally asymmetric (red). (B) Pathlines for particles advected by the double gyre flow, generated using the function \texttt{pathPlot}. (C) The motion of a blob of particles advected by the flow, with successive timepoints indicated in darkening shades of blue. Plot generated using the function \texttt{advectPoints}. (D) The maximal finite-time Lyapunov exponents, generated using the function \texttt{findFTLE}. (E) The Kaplan-Yorke fractal dimension, generated using the function \texttt{findKYDim}. (F) The flushing times for the shown domain, generated using the function \texttt{flushingTimes}. For these images, the double gyre flow was used, with parameters $A=0.1,\epsilon=0.25,\omega=\pi/5$.
}
\label{fig_lagrange}
}
\end{figure}

\subsubsection{Define a velocity field}
Figure \ref{fig_lagrange} shows applications of the functions described in the previous section to the study of the double gyre flow, which has the form
\begin{align*}
v_x(x,y,t) &= - \pi A \sin(\pi f(x,t)) \cos(\pi y), 			\\
v_y(x,y,t) &= \pi A \cos(\pi f(x,t)) \sin(\pi y) \big(2 a(t) x + b(t)\big)
\end{align*}
where
\begin{align*}
a(t) &= \epsilon \sin(\omega t), \\
b(t) &= 1-2\epsilon \sin(\omega t), \\
f(x,t) &= a(t) x^2 + b(t) x.
\end{align*}
The parameters $\omega$ and $A$ control the transport timescales of the flow, while $\epsilon$ determines the degree of temporal variation, and thus chaotic mixing present \cite{shadden2005definition}.

When using \texttt{Lagrange2D}, this flow may be defined as a standard Mathematica function,
\begin{verbatim}
a[t_] := e Sin[w t];
b[t_] := 1 - 2 e Sin[w t];
f[t_, x_] := a[t] x^2 + b[t] x;

g[t_, x_, y_] := Pi A Sin[Pi f[t, x]] Cos[Pi f[t, y]](*D[f[t,x],x]*)

(* double gyre flow *)
{vx[t_, x_, y_],  vy[t_, x_,  y_]} :=
     {-Pi A Sin[Pi f[t, x]] Cos[Pi y],
       Pi A Cos[Pi f[t, x]] Sin[Pi y] (2 a[t] x + b[t])};

params = {A -> 0.1, e -> 0.1, w -> Pi/5};
\end{verbatim}

\subsection{Visualize and simulate a time-varying flow}

Plotting the pathlines associated with a flow uses similar syntax to Mathematica's built-in \texttt{StreamPlot} function,
\begin{verbatim}
pathPlot[{vx[t, x, y], vy[t, x, y]} /. params,
  {x, 0, 2}, {y, 0, 1}, {t, 0, 15}]
\end{verbatim}
The density and spacing of pathlines can be manipulated using similar optional arguments to those used by \texttt{StreamPlot}, which are further described in the function's docstring.

Similar syntax is used to integrate the motion of a set of initial conditions subject to the flow. However, the full initial conditions must be defined by the user in a separate command,
\begin{verbatim}
blobPoints = makeMesh[{.2, .6}, {.1, .5}, 50, 50];
ptLocs = advectPoints[
   blobPoints,
   {vx[t, x, y], vy[t, x, y]} /. params,
   {t, 0, 20}, x, y];
\end{verbatim}

The motion of this collection of points may be visualized by exporting a video,
\begin{verbatim}
animateFlow[
  {vx[t, x, y], vy[t, x, y]} /. params,
  blobPoints, {t, 0, 20},{x, 0, 2}, {y, 0, 1}, 20, 1/10];
\end{verbatim}

\subsection{Find the maximal FTLE field}

Finding the FTLE field associated with a velocity field uses similar syntax to the commands described above. This function requires the exact mesh points at which to compute the FTLE field to be specified by the user.
\begin{verbatim}
domainPoints = makeMesh[{0, 2}, {0, 1}, 150, 150];
ftleField =
  findMaxFTLEField[{vx[t, x, y], vy[t, x, y]} /. params,
   domainPoints, {t, 0, 10}, x, y];
\end{verbatim}

\subsection{Find a Kaplan-Yorke fractal dimension field}

In order to calculate the Kaplan-Yorke fractal dimension, two commands are required: one to generate the finite-time Lyapunov exponents associated with each point in a user-specified mesh, and another to convert these exponents into an estimate of the Kaplan-Yorke fractal dimension,
\begin{verbatim}
domainPts = makeMesh[{0, 2}, {0, 1}, 150, 150];
ftleVals =
  findFTLEField[{vx[t, x, y], vy[t, x, y]} /. params,
   domainPoints, {t, 0, 10}, x, y];
kyVals = findKYDim[ftleVals];
\end{verbatim}

\subsection{Find the stretching vector field}

The stretching distribution corresponds to a vector field defined at each point in the domain. Similar to the above commands for computing FTLE fields, this function requires a mesh to be specified by the user. There are two orthogonal stretching fields, one associated with the rate of maximal stretching, and the other associated with the rate of minimal stretching,
\begin{verbatim}
blobPoints = makeMesh[{0, 2}, {0, 1}, 100, 100];
stretchVecsPositive =
  findStretchlines[{vx[t, x, y], vy[t, x, y]} /. params,
   blobPoints, {t, 0, 10}, x, y, "Positive"];
stretchVecsNegative =
  findStretchlines[{vx[t, x, y], vy[t, x, y]} /. params,
   blobPoints, {t, 0, 10}, x, y, "Negative"];
\end{verbatim}

\subsection{Find the flushing time field}

The flushing, or residence, time distribution has similar syntax to Mathematica's built-in \texttt{StreamPlot} function. However, because the boundaries of the domain determine the residence time, this field depends directly on the bounds of the specified mesh.
\begin{verbatim}
flushField = flushingTimes[
   {vx[t, x, y], vy[t, x, y]} /. params,
   {x, -1, 1},{y, 0, 1},{t, 0, 50},{n -> 20000}];
\end{verbatim}

The results of some of the functions described above are showing in Figure \ref{fig_lagrange}. Comparison of snapshots of advected particle positions, and the FTLE field, reveals that the FTLE field successfully reveals kinetic barriers to transport in the flow. Comparison of the FTLE field and advection of a group of particles shows that groups particles tend to deform around regions of high FTLE values. The Kaplan-Yorke dimension indicates regions of the flow that are kinetically accessible within the integration time, and it establishes regions near which trajectories tend to cluster. Consistent with the pathlines of the flow, the flushing time plot reveals asymmetry intrinsic to the map, indicated by strong differences in the residence time of particles depending on whether they originate in regions to the left or right of the unstable manifold of the flow.

\section{Conclusions}

We present a Mathematica package, \texttt{Lagrange2D}, for the analysis and visualization of transport in two-dimensional fluid flows. Our package has particular relevance to characterizing unsteady flows, for which standard velocity fields and streamline visualization can fail to capture the full structure of the flow. \texttt{Lagrange2D} includes utilities for integrating ensembles of trajectories, as well as functions for computing finite-time Lyapunov exponents, the fractal dimension, and the flushing time distribution. Additionally, our package contains visualization tools for viewing pathlines and for generating videos of time-varying flows. \texttt{Lagrange2D} is self-contained and simple to install, and we hope that it will simplify the characterization of two-dimensional flows, particularly those for which the velocity field has an analytic form.

\section{Acknowledgements}

W.G. thanks the U. S. Department of Defense for their support through the NDSEG Fellowship program. The author declares no competing interests.

\bibliographystyle{elsarticle-num}


\section{Introduction}
In the last decades, automatic analysis of human motion has been an active research topic, with applications that have been exploited in a number of different contexts, including video surveillance, semantic annotation of videos, entertainment, human computer interaction and home care rehabilitation, to say a few.
Differences in body proportion (size, height, corpulence), body stiffness and training, influence the way different people perform an action. Even one same person is not able to perform the same action twice, exactly replicating the same sequence of body poses in space and time. This variability makes the task of human motion analysis very challenging.

For years, the approaches could be distinguished in two main classes: those operating on pixel values extracted from the RGB stream (either stacking groups of consecutive frames or extracting motion vectors) and those building upon the higher level representation of body skeletons.
These latter approaches were supported by the diffusion of low-cost RGB-D cameras (such as the Microsoft Kinect) that can operate in real-time, while reliably extracting the 3D coordinates of body joints.
More recently, deep CNN architectures have demonstrated real-time and accurate extraction of the coordinates of body joints from RGB streams~\cite{cao2017realtime}.
These advances make it possible to use a skeleton-based body representation in a much broader range of domains and operative contexts than before, being not limited by the short operative range of RGB-D sensors that typically operate indoor and in the range of a few meters.
The design of the recognition/classification module on top of the body skeleton representation makes it possible to describe an action as a sequence of body poses, each one corresponding to a point in a feature space, whose dimension is proportional to the number of body joints.
By exploiting the geometric properties of the manifold where these pose descriptors lie, it is possible to define a similarity metric that is invariant under translation, scaling, rotation and also under variations of the speed of execution of the action.
Furthermore, the explicit representation of an action as a trajectory, \emph{i.e.}, a sequence of poses, on the manifold makes it possible to extract statistical summaries, such as mean and deviation from the mean, from a group of actions.
Through these summaries, one action can be better characterized for the purpose of detecting outliers corresponding to the anomalous execution of an action, that can be of particular relevance for action prediction.
In fact, when analyzing the skeleton sequences, there are four main aspects to challenge: (1) A shape representation invariant to undesirable transformations; (2) A temporal modeling of landmark sequences; (3) A suitable rate-invariant distance between arbitrary sequences, and (4) A solution for temporal sequence classification.

This paper lies in the continuity of recent works that model the comparison and classification of temporal sequences of landmarks on the Riemannian manifold of positive-semidefinite matrices. Building on the work~\cite{KacemPAMI2019}, our approach involves four different steps: 1) We build a trajectory on the Riemannian manifold from the body skeletons; 2) We apply a curve fitting algorithm on the trajectories to denoise the data points; 3) We perform a temporal alignment using a Global Alignment Kernel, defining a positive-semidefinite kernel; 4) Finally, we use this kernel with a classic SVM to classify the actions. An overview of the full approach is given in Fig.~\ref{Fig:AppOverview}.

\begin{figure*}[!ht]
  \centering
    \includegraphics[width=0.85\textwidth]{Overview-ICCV.png}
    \caption{Overview of the proposed approach - After automatic body skeletons detection for each frame of a sequence, the Gram matrices are computed to build the trajectory on the $\mathcal{S}^+(d,n)$ manifold. We apply a curve fitting algorithm on the trajectory to smooth the curve and reduce noise. Global Alignment Kernel (GAK) is then used to align the trajectories on the manifold. Finally, we use the kernel generated from GAK with SVM to classify the actions.}
    \label{Fig:AppOverview}
\end{figure*}

The novelties with respect to~\cite{KacemPAMI2019} are:
\begin{itemize}
\item The manifold of positive-semidefinite matrices is here endowed with a different metric;
\item A recent curve fitting method is used to smooth trajectories on the manifold;
\item We use Global Alignment Kernel for temporal alignment, instead of Dynamic Time Warping.
\end{itemize}

The rest of the paper is organized as follows: in Section~\ref{sect:2}, we summarize the closest works in the literature. In Section~\ref{sect:3}, we explain how we represent our data on the manifold of fixed-rank positive-semidefinite matrices, and present the Riemannian metric that we will use in this paper. In Section~\ref{sect:4}, we describe the curve fitting algorithm that we will use to denoise the trajectories. The curve alignment method and the classifier are presented in Section~\ref{sect5:classification}, while results and discussions are reported in Section~\ref{sect:exp-results}. Finally, in Section~\ref{sect:conclusion}, we conclude and draw some perspectives of the work.

\section{Related Work}\label{sect:2}
A detailed review of the many approaches to human action recognition and classification is out of the scope of this paper. The interested reader can refer to~\cite{kong-2018} for a detailed and updated survey.
In the following, we focus on approaches that use body skeletons as inputs to the recognition/classification module.

One of the first approaches to perform action recognition by the analysis of trajectories of tracked body interest points was presented in Matikainen \etal~\cite{Matikainen_2009_6466}.
Despite the promising results obtained, the authors did not take into account the geometric information of the trajectories.
More recently, in the case of human skeleton in RGB-D images, Devanne \etal~\cite{Devanne2015Cybernetics}
proposed to formulate the action recognition task as the problem of computing a distance between trajectories generated by the joints moving during the action. An action is then interpreted as a normalized parameterized curve in $\mathbb{R}^N$. However, this approach does not take into account the relationship between the joints.
In the same direction, Su \etal~\cite{Su2014} proposed a metric that considers time-warping on a Riemannian manifold, thus allowing the registration of trajectories and the computation of statistics on the trajectories. Su \etal~\cite{Su_2014_CVPR} applied this framework to the problem of visual speech recognition. Similar ideas have been developed by Ben Amor \etal{}~\cite{amor2016action} on the Kendall's shape space with application to action recognition using rate-invariant analysis of skeletal shape trajectories.

Anirudh \etal~\cite{AnirudhTSS17} started from the framework of Transported Square-Root Velocity Fields (TSRVF), which has desirable properties including a rate-invariant metric and vector space representation. Based on this framework, they proposed to learn an embedding such that each action trajectory is mapped to a single point in a low-dimensional Euclidean space, and the trajectories that differ only in the temporal rate map to the same point. The TSRVF representation and accompanying statistical summaries of Riemannian trajectories are used to extend existing coding methods such as PCA, KSVD, and Label Consistent KSVD to Riemannian trajectories. In the experiments, it is demonstrated that such coding efficiently captures distinguishing features of the trajectories, enabling action recognition, stroke rehabilitation, visual speech recognition, clustering, and diverse sequence sampling.

In~\cite{vemulapalli2014human}, Vemulapalli \etal proposed a Lie group trajectory representation of the skeletal data on a product space of special Euclidean ($SE$) groups. For each frame, this representation is obtained by computing the Euclidean transformation matrices encoding rotations and translations between different joint pairs. The temporal evolution of these matrices is seen as a trajectory on $SE(3) \times \cdots \times SE(3)$ and mapped to the tangent space of a reference point. A one-versus-all SVM, combined with Dynamic Time Warping and Fourier Temporal Pyramid (FTP) is used for classification. One limitation of this method is that mapping trajectories to a common tangent space using the logarithm map could result in significant approximation errors. Aware of this limitation, in~\cite{vemulapalli2016rolling} the same authors proposed a mapping combining the usual logarithm map with a rolling map that guarantees a better flattening of trajectories on Lie groups.

More recently, Kacem \etal~\cite{KacemPAMI2019} proposed a geometric approach for modeling and classifying dynamic 2D and 3D landmark sequences based on Gramian matrices derived from the static landmarks. This results in an affine-invariant representation of the data. Since Gramian matrices are positive-semidefinite, the authors relies on the geometry of the manifold of fixed-rank positive-semidefinite matrices, and more specifically, to the metric investigated in~\cite{Bonnabel2009SIAM}. However, this metric is parametrized, and the parameter should ideally be learned from the data. In addition, this paper adopts Dynamic Time Warping for sequence alignment. The resulting distance does not generally lead to a positive-definite kernel for classification.

All the approaches described above rely on the use of \textit{hand-crafted features} enabling representation of the action as a trajectory or point in some suitable manifold. Differently from these approaches, many neural network models have been proposed that rely on training for the extraction of \textit{deep learned features}.
Recurrent neural networks (RNNs) and particularly Long-Short-Term Memory Networks (LSTMs) have been used to perform action recognition by the analysis of sequences of skeleton poses~\cite{Zhu-2016}. However, these methods typically lose structural information when converting the skeleton data and joint connectivity into the vector-shaped input of the neural network. As an alternative, some approaches introduce Convolutional Neural Networks (CNNs)~\cite{Chao-2017} and Graph Convolutional Networks (GCNs)~\cite{Yan-2018,Li-2019} in the overall architecture so as to retain the structural information among joints of the skeleton.
Although these approaches result in state-of-the-art performance~\cite{Gao-2019} on public action recognition benchmarks, it is not possible to define a formal mathematical framework to compute a valid metric on the internal, learned feature representation so as to perform a statistical analysis of the learned actions.

\section{Our Approach}\label{sect:3}
\subsection{Shape Representation}
To represent body movement dynamics, we rely on the time series made of the coordinates of the $n$ tracked body points (\emph{i.e.}, $p_1 = (x_1,y_1), \ldots, p_{n} = (x_{n},y_{n})$ in 2D, or $p_1 = (x_1,y_1,z_1), \ldots, p_{n} = (x_{n},y_{n},z_n)$ in 3D), during each video sequence.
Each video sequence is thus characterized by a set of landmark configurations $\{Z_0,\ldots, Z_{\tau}\}$, where $\tau$ is the number of frames of the video sequence, and where each configuration matrix $Z_{i}$ $(1 \leq i\leq \tau) \in  \mathrm{R}^{n \times d}$ encodes the position of the $n$ landmarks in $d$ dimensions (with $d = 2$ or $d = 3$). We aim to measure the dynamic changes of the curves made of the landmark configurations, remaining invariant to rotation and translation.
Similarly as in~\cite{KacemPAMI2019}, this goal is achieved through a Gram matrix representation, where we compute the Gram matrices as:
\begin{equation}
\label{eq:gram}
G=Z Z^T\; .
\end{equation}

These Gram matrices are $n \times n$ positive-semidefinite matrices, of rank smaller than or equal to $d$ (always equal to $d$ in the datasets considered). Conveniently for us, the Riemannian geometry of the space ${\mathcal S}^+(d,n)$ of $n \times n$ positive-semidefinite matrices of rank $d$ has been studied in~\cite{Bonnabel2009SIAM,journee2010low,vandereycken2009embedded,vandereycken2013riemannian,Massart2018, massart2019curvature}, and used in, \emph{e.g.},~\cite{faraki2016image,meyer2011regression,gousenbourger2017piecewise, massart2019interpolation}.

A classical approach in the design of algorithms on manifolds consists in resorting to first order local approximations on the manifold, called tangent spaces. This requires two tools: the Riemannian exponential (that allows us to map tangent vectors from the tangent space to the manifold), and the Riemannian logarithm (mapping points from the manifold to the tangent space).

In~\cite{KacemPAMI2019}, the manifold $\psd{d}{n}$ is identified to the quotient manifold $(\mathrm{St}(d,n) \times \mathcal{P}_d)/\mathcal{O}_d$, where $\mathrm{St}(d,n) := \{  Y \in \mathbb{R}^{n \times d} \rvert Y^T Y = I_d \}$ is the Stiefel manifold, $\mathcal{P}_d$ is the manifold of $d \times d$ positive-definite matrices, and $\mathcal{O}_d$ is the orthogonal group in dimension $d$. We consider here another representation of the manifold $\psd{d}{n}$, that will result in different expressions for the distance between two points, the Riemannian exponential and logarithm.

\subsection{The Quotient Manifold $\frmat{n}{d} / \mathcal{O}_d$}\label{metrics}

We consider here the identification of $\psd{d}{n}$ to the quotient manifold $\frmat{n}{d} / \mathcal{O}_d$, where $\frmat{n}{d}$ is the set of full-rank $n \times d$ matrices. This geometry has been studied in~\cite{journee2010low,Massart2018,massart2019curvature}.

The identification of $\psd{d}{n}$ with the quotient $\frmat{n}{d} / \mathcal{O}_d$ comes from the following observation. Any PSD matrix $G \in \psd{d}{n}$ can be factorized as $G = ZZ^T$, with $Z \in \frmat{n}{d}$. However, this factorization is not unique, as any matrix $\tilde Z := ZQ$, with $Q \in \mathcal{O}_d$, satisfies $\tilde Z \tilde Z^T = Z QQ^T Z^T = G$. The two points $Z$ and $\tilde Z$ are thus \emph{equivalent} with respect to this factorization, and the set of equivalent points
\[
Z \mathcal{O}_d := \{  ZQ | Q \in \mathcal{O}_d \},
\]

\noindent
is called the equivalence class associated to $G$. The quotient manifold $\frmat{n}{d} / \mathcal{O}_d$ is defined as the set of equivalence classes. The mapping $\pi : \frmat{n}{d} \to \frmat{n}{d} / \mathcal{O}_d$, mapping points to their equivalence class, induces a Riemannian metric on the quotient manifold from the Euclidean metric in $\frmat{n}{d}$. This metric results in the following distance between PSD matrices:
\begin{equation}
\label{eq:PSD-norm}
    d(G_i, G_j) = \left[  \mathrm{tr}(G_i) + \mathrm{tr}(G_j) - 2 \mathrm{tr}\left( \left( G_i^{\frac{1}{2}}  G_j G_i^{\frac{1}{2}} \right)^{\frac{1}{2}} \right) \right]^{\frac{1}{2}}.
\end{equation}

This distance can be expressed in terms of the landmark variables $Z_i, Z_j \in \frmat{n}{d}$ as follows:
\begin{equation}
\label{eq:frob-norm-new-metrix}
d(G_i, G_j) = \min_{Q \in \mathcal{O}_d} \lVert Z_jQ-Z_i\rVert_F.
\end{equation}

\noindent
The optimal solution is $Q^* := VU^\top \!$, where $Z_i^\top \! Z_j = U \Sigma V^\top \!$ is a singular value decomposition.

As stated by the next theorem, when $d = 2$, the distance can also be formulated as follows:

\begin{theorem}
\label{theo1}
Let $G_i, G_j \in \psd{2}{n}$ be two Gram matrices, obtained from landmark matrices $Z_i, Z_j \in \mathbb{R}^{n \times 2}$. The Riemannian distance~\eqref{eq:PSD-norm} can be expressed as:
\begin{equation}
\label{eq:distp2}
\resizebox{1\hsize}{!}{$d(G_i, G_j)= \left( \Tr(G_i) - 2\sqrt{(a + d)^2 + (c - b)^2} + \Tr(G_j) \right)^{\frac{1}{2}}$},
\end{equation}
\end{theorem}

\noindent
where $Z_j^TZ_i$ = $\left( \begin{array}{cc} a & b \\ c & d \end{array} \right)$.
\begin{proof}
See Appendix.
\end{proof}

Expressions for the Riemannian exponential and logarithm are given in~\cite{Massart2018}. We used the implementations provided in the Manopt toolbox~\cite{manopt}.

\section{Trajectory Modeling}\label{sect:4}
The dynamic changes of body joints movement are characterized by trajectories on the Riemannian manifold of positive-semidefinite matrices of fixed rank (see Fig.~\ref{Fig:AppOverview}). More specifically, we fit a curve $\beta_G$ to a sequence  of landmark configurations $\{Z_0, \ldots, Z_\tau\}$ represented by their corresponding Gram matrices $\{G_0,\ldots,G_\tau\}$ in $\mathcal{S}^+(d,n)$. This curve will enable us to model the spatio-temporal evolution of the elements on $\mathcal{S}^+(d,n)$.

Modeling a sequence of landmarks as a piecewise-geodesic curve on $\mathcal{S}^+(d,n)$ showed very promising results when the data are well acquired, \emph{i.e.}, without tracking errors or missing data, see~\cite{KacemPAMI2019, KacemICCV17, OtberdoutBMVC2018}. To account for both missing data and tracking errors, we rely on a more recent curve fitting algorithm: fitting by composite cubic blended curves, proposed in~\cite[\S 5]{Gousenbourger2018}. Specifically, given a set of points $G_0,\ldots,G_\tau \in \mathcal{S}^+(d,n)$ associated to times $t_0, \dots, t_\tau$, with $t_i := i$, the curve $\beta_G$, defined on the interval $[0,\tau]$, is defined as:
\[
\beta_G(t) := \gamma_i(t-i), \qquad t \in [i, i+1],
\]

\noindent
where each curve $\gamma_i$ is obtained by blending together fitting cubic B\'ezier curves computed on the tangent spaces based on the data points $d_i$ and $d_{i+1}$.

These fitting cubic B\'ezier curves depend on a parameter $\lambda$, allowing us to balance two objectives: proximity to the data points at the associated time instants, and regularity of the curve (measured in terms of mean square acceleration). A high value of $\lambda$ results in a curve with possibly high acceleration, but that will almost interpolate the data, while taking $\lambda \rightarrow 0$ will result in a geodesic. The interested reader can refer to~\cite[\S 5]{Gousenbourger2018} for more information about the curve fitting procedure.

\section{Classification}\label{sect5:classification}
Now that we have defined how to represent a sequence and how to compare two distinctive landmark configurations, we present in this section how we compare two landmark sequences and how to classify the actions performed in these same sequences.

\subsection{Global Alignment}\label{global-alignment-kernel}
As we described in Section~\ref{sect:4}, we represent a sequence as a trajectory of Gram matrices in $\psd{d}{n}$. The sequences represented in this manifold can be of different length as the execution rate of the actions can vary from one person to another, meaning that we can not effectively compare them. A common method to do so is to use Dynamic Time Warping (DTW) as proposed in several works~\cite{amor2016action, KacemPAMI2019, gritai2009matching}. However, DTW does not define a proper metric and can not be used to derive a valid positive-definite kernel for the classification phase. To address the problem of non positive definiteness of the kernel defined by DTW, Cuturi \etal~\cite{CuturiVBM07} proposed the Global Alignment Kernel (GAK), which allows us to derive a valid positive-definite kernel when aligning two time series. More recently Otberdout \etal~\cite{Otberdout2019} have proposed to classifiy deep trajectories in SPD manifold using GAK. The generated kernel can be used directly with Support Vector Machine (SVM) for the classification phase, whereas it is not the case with kernels generated with DTW. In fact, the kernels built with DTW do not show favorable positive definiteness properties as they rely on the computation of an optimum rather than the construction of a feature map. Note that the computation of the kernels with GAK can be done in quadratic complexity, similarly to naive implementation of DTW. The next paragraph describes how to compute the similarity score between two sequences, using this Global Alignment Kernel.

Let us now consider $Z^1= \{Z^1_0,\cdots, Z^1_{\tau_1}\}$ and $Z^2=\{Z^2_0,\cdots, Z^2_{\tau_2}\}$, two sequences of landmark configuration matrices. Given a metric to compute the distance between two elements of each sequence, we propose to compute the matrix $D$ of size $\tau_1 \times \tau_2$, where each $D(i,j)$ is the distance between two elements of the sequences, with $1 \leq i \leq \tau_1$ and $1 \leq j \leq \tau_2$.
\begin{equation}
\label{eq:matrix-distance}
    D(i,j) = d(Z^1_i, Z^2_j) .
\end{equation}

The kernel $\tilde{k}$ can now be computed using the halved Gaussian Kernel on this same matrix $D$. Therefore, the kernel $\tilde{k}$ can be defined as:
\begin{equation}
\label{eq:k-tilde}
    \tilde{k}(i,j) = \frac{1}{2} * exp\left(-\frac{D(i,j)}{\sigma^2}\right).
\end{equation}
 As reported in~\cite{CuturiVBM07}, we can redefine our kernel such as:
\begin{equation}
\label{eq:kernel-PSD}
    k(i,j) = \frac{\tilde{k}(i,j)}{(1-\tilde{k}(i,j))} .
\end{equation}
This strategy assures us that the kernel is positive semi-definite and can be used in its own. Finally, we can compute the similarity score between the two sequences $Z^1$ and $Z^2$. Remember that this computation is performed in quadratic complexity, like DTW. To do so, we define a new matrix $M$ that will contain the path to the similarity between our two sequences. We define $M$ as a zeros matrix of size $(\tau_1+1) \times (\tau_2+1)$ and $M_{0,0} = 1$. Computing the terms of $M$ is done using Theorem 2 in~\cite[\S 2.3]{CuturiVBM07}:
\begin{equation}
\label{eq:matrix-similarity}
    M_{i,j} = (M_{i,j-1} + M_{i-1,j-1} + M_{i-1,j})*k(i,j) .
\end{equation}

The similarity score we seek is the value at $M_{(\tau_1+1),(\tau_2+1)}$. Algorithm~\ref{algo:similarity} describes all the steps to get the similarity score.

Finally, we build a new matrix $K$ of size $n_{seq} \times n_{seq}$, where $n_{seq}$ is the number of sequences in the dataset we test. This matrix is symmetric and contains all the similarity scores between all the sequences of the dataset and it is used as the kernel for the classification phase with SVM. As this matrix is built with values computed from positive semi-definite kernel, it is a positive semi-definite matrix itself.

\begin{algorithm}
\caption{Computing the similarity score between two sequences using Global Alignment Kernel~\cite{CuturiVBM07}
\label{algo:similarity}}
\DontPrintSemicolon
\SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
\Input{Two sequences of landmark configurations $Z^1= \{Z^1_0,\cdots, Z^1_{\tau_1}\}$, where $Z^1_{0 \leq i \leq \tau_1}$ and $Z^2=\{Z^2_0,\cdots, Z^2_{\tau_2}\}$, where $Z^2_{0 \leq j \leq \tau_2}$.}
\Output{The similarity score between two sequences $Z^1, Z^2$}
$\tilde{k} \longleftarrow \frac{1}{2}*exp\left(-\frac{D(Z^1, Z^2)}{\sigma^2} \right)$ Equations ~\eqref{eq:matrix-distance} and \eqref{eq:k-tilde} \\

\For{$i \leftarrow 0$ \KwTo $\tau_1$}{
    \For{$j \leftarrow 0$ \KwTo $\tau_2$}{
        $k(i,j) \longleftarrow \frac{\tilde{k}(i,j)}{(1 - \tilde{k}(i,j))}$ Equation~\eqref{eq:kernel-PSD}
    }
}

$M \longleftarrow zeros(\tau_1+1, \tau_2+1)$ \\
$M_{0,0} \longleftarrow 1$ \\

\For{$i \leftarrow 1$ \KwTo $\tau_1+1$}{
    \For{$j \leftarrow 1$ \KwTo $\tau_2+1$}{
        $M_{i,j} \longleftarrow (M_{i,j-1} + M_{i-1,j-1} + M_{i-1,j})*k(i,j)$ See Equation~\eqref{eq:matrix-similarity}
    }
}
$similarity \longleftarrow M_{\tau_1+1, \tau_2+1}$

\Return similarity, the similarity score between $Z^1$ and $Z^2$
\end{algorithm}

\subsection{Classification with SVM}
Our trajectory representation reduces the problem of landmark sequence classification to that of trajectory classification in $\mathcal{S}^+(d,n)$. Given that GAK provides a valid PSD kernel as demonstrated by Cuturi \etal{}~\cite{CuturiVBM07}, and given that our local kernel $K$ satisfies this condition as discussed before, we use the standard SVM with the $K$ kernel that represents the matrix containing the similarity scores between all the sequences of a dataset to classify the aligned trajectories with global alignment on $\mathcal{S}^+(d,n)$.

By contrast, DTW may define a non positive definite kernel. Hence, we adopt the pairwise proximity function SVM (ppfSVM), which assumes that instead of a valid kernel function, all that is available is a proximity function without any restriction.
That is, let us consider $\mathcal{T} = \{\beta_G~:[0,1]\rightarrow \mathcal{S}^+(d,n)\}$, the set of time-parameterized trajectories of the underlying manifold. Like in~\cite[\S 4.1]{KacemPAMI2019}, we define a matrix $D_{dtw}$ containing the similarity measure between two trajectories aligned with DTW.
In that case, given $m$ trajectories in $\mathcal{T}$, the proximity function $\mathcal{P} : \mathcal{T} \times \mathcal{T} \rightarrow \mathbb{R}^+$ between two trajectories $Z^{1}$ and $Z^{2}$ is defined by,
\begin{equation}
\mathcal{P}(Z^{1}, Z^{2})=D_{dtw}(Z^{1}, Z^{2}) \; .
\end{equation}

\noindent
Using this proximity function, the main idea of ppfSVM is to represent each training example $Z$ with a vector $[\mathcal{P}(Z,Z^1), \dots , \mathcal{P}(Z,Z^m)]^T$. The set of trajectories can be represented by a $m \times m$ matrix $P$, where $P(i,j) = \mathcal{P}(Z^{1}, Z^{2})$, with $1 \leq i,j \leq m$. From this matrix $P$ we can use a classical linear SVM.

\section{Experimental Results}
\label{sect:exp-results}
To validate the proposed approach, we have conducted experiments on three publicly available datasets with 3D and 2D actions: UTKinect-Action3D, KTH-Action and UAV-Gesture. Our experiments followed the experimental settings commonly used for these datasets.

\subsection{UTKinect-Action3D Dataset}
The UTKinect-Action3D dataset~\cite{xia2012view} is a widely used dataset for 3D action recognition. It contains 199 sequences, consisting of 10 actions, namely \textit{walk, sit down, stand up, pick up, carry, throw, push, pull, wave hands} and \textit{clap hands} performed by 10 different subjects. The videos and the skeletons were captured with a Microsoft Kinect and the skeletons are composed of 20 body joints. In our approach, we use the available skeletal joint locations, where each body joint is defined with its $x$, $y$ and $z$ coordinates. Following the same experimental settings of~\cite{TanfousDA18, LiuWDAK18, KeBASB18}, we performed the Leave-One-Out cross validation, meaning that we used one sequence for testing and the rest for training.
Our experimental results are summarized in Table~\ref{results-our-utk}. In particular, the columns are as follows: \emph{Curve Fitting} indicates if we performed the curve fitting algorithm described in Section~\ref{sect:4}; \emph{Lambda} indicates the value of the lambda parameter in curve fitting;
\emph{Alignment Method} indicates if we used the standard DTW to align sequences or GAK as described in Section~\ref{global-alignment-kernel}; \emph{Sigma} indicates the value of the sigma parameter for the Gaussian Kernel when using GAK; and \emph{Results} indicates our scores.

\begin{table}[!ht]
\centering
\caption{Our results on the UTKinect-Action3D dataset}
\begin{adjustbox}{width=0.47\textwidth}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Curve Fitting} & \textbf{Lambda} & \textbf{Alignment Method} & \textbf{Sigma} & \textbf{Results} \\ \hline
Yes & 0.5 & DTW & - & 97\% \\
No & - & DTW & - & 97.49\% \\ \hline
Yes & 0.5 & GAK & 0.3 & 97.49\% \\
No & - & GAK & 0.3 & \textbf{97.99\%} \\
Yes & 0.5 & GAK & 0.5 & \textbf{97.99\%} \\ \hline
\end{tabular}
\end{adjustbox}
\label{results-our-utk}
\end{table}

The best accuracy that we obtained on this dataset is 97.99\%. Overall, we can say that the application of curve fitting does not increase our results. Our assumption is that the data in this dataset are very clean, and we can loose some information with the application of smoothing on clean data. Note that we obtained better results when using the Global Alignment Kernel rather that DTW.

\begin{table}[H]
\centering
\caption{Comparison of our approach with state-of-the-art results for the UTKinect-Action3D dataset. *: Deep Learning approach}
\begin{adjustbox}{width=0.38\textwidth}
\begin{tabular}{|l|c|c|}
\hline
 & \multicolumn{2}{c|}{\textbf{Protocol}} \\
 \hline
\textbf{Methods} & \textbf{H-H} & \textbf{LOOCV} \\ \hline
Trajectory on ${\mathcal S}^+(d,n)$~\cite{KacemPAMI2019} (2019) & - & 96,48\% \\
SCK+DCK~\cite{KoniuszCP16} (2016) & 98.2\% & - \\
Bi-LSTM~\cite{TanfousDA18} (2018)* & - & 98.49\% \\
LM$^3$TL~\cite{YangDTZLG17} (2017) & - & 98.8\% \\
GCA-LSTM~\cite{LiuWDAK18} (2018)* & - & 99\% \\
MTCNN~\cite{KeBASB18} (2018)* & - & 99\% \\
Hankel \& Gram matrices~\cite{ZhangWGSC16} (2016) & - & \textbf{100\%} \\
\hline
Ours & - & 97.99\% \\
\hline
\end{tabular}
\end{adjustbox}
\label{sota-utk}
\end{table}

In Table~\ref{sota-utk}, we compare our method with recent state-of-the-art results. Overall, our approach achieves competitive results with respect to most recent approaches. We directly compare our results with~\cite{KacemPAMI2019} as we work on the same geometric space of ${\mathcal S}^+(d,n)$ manifold. The main differences between our method and the method in~\cite{KacemPAMI2019} is the use of a different metric and of the Global Alignment Kernel instead of DTW. Our metric is simpler that the metric in~\cite{KacemPAMI2019}, as we do not have to estimate the parameter $k$ used in Eq.~(7) in~\cite{KacemPAMI2019} for distance computation. Furthermore, the $k$ parameter in~\cite{KacemPAMI2019} is more of a constraint as they have to determine its best value for each dataset they test. The use of GAK is also an advantage for us as it defines a positive semi-definite kernel, which is not the case for DTW allowing us to use a classic SVM instead of ppfSVM.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.47\textwidth]{UTK-Comparison-Throw-Push.png}
    \caption{Comparison of two sequences that are confused in UTKinect-Action3D dataset (top: \textit{Throw} action, bottom: \textit{Push} action)}
    \label{qualitative-utk}
\end{figure}

The Figure~\ref{qualitative-utk} presents two sequences that are confused, leading to a misclassification for one of them. In that case, the top action (\emph{i.e.} Throw) is misclassified as the bottom action (\emph{i.e.} Push). One of the reasons can be the position of the arm at the end of the action, which is the same in the two sequences. The \textit{Throw} action is the most confused action in the dataset.

\subsection{KTH-Action Dataset}
The KTH-Action dataset~\cite{kth} is a 2D action recognition dataset. It consists of six actions, namely \textit{boxing, handclapping, handwaving, jogging, running} and \textit{walking} performed by 25 subjects in four different conditions, which are outdoor, outdoor with scale variations, outdoor with different clothes and indoor. The sequences were acquired with a static camera at a frame rate of 25 fps and a resolution of $160 \times 120$ pixels. The dataset contains a total of 599 clips, with 100 clips per actions (1 clip is missing for one action). As the sequences in the dataset are 2D videos, we have to extract the skeletons of the subjects performing the actions. To do so, we used the OpenPose framework~\cite{openpose} to extract the skeletons in the COCO format, with 18 body joints. Note that we clean the landmark sequences by removing the frames where the body joints where not effectively estimated. Keeping all the frames leads to worst results due to misdetected joints, meaning that we do not need all the frames available to recognize an action. Figure~\ref{coco-skeleton} shows the configuration of the body joints that we analyzed. For this dataset, we followed the Leave-One-Actor-Out cross validation protocol, meaning that we use one subject for testing and the rest for training. Table~\ref{results-ours-kth} summarizes our experimental results on this dataset.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.22\textwidth]{coco-body-joints.png}
    \caption{Skeleton with the COCO format.}
    \label{coco-skeleton}
\end{figure}

\begin{table}[H]
\centering
\caption{Our experimental results on the KTH-Action dataset}
\begin{adjustbox}{width=0.47\textwidth}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Curve Fitting} & \textbf{Lambda} & \textbf{Alignment Method} & \textbf{Sigma} & \textbf{Results} \\ \hline
No & - & DTW & - & 94.49\% \\
Yes & 10 & DTW & - & 94.66\% \\ \hline

No & - & GAK & 0.2 & 95.16\% \\
Yes & 10 & GAK & 0.2 & \textbf{96.16\%} \\ \hline
\end{tabular}
\end{adjustbox}
\label{results-ours-kth}
\end{table}

Here, again, we obtained better results when using the GAK, demonstrating superior performance over DTW. The results reported with DTW are the best accuracy over all the configurations we tested. Unlike the data in UTKinect-Action3D dataset, the data in KTH-Action are 2D and low resolution videos, with presence of noise in the background, leading to noisy skeleton data after extraction. In this regard, the application of the curve fitting algorithm improves our results by 1\%. We compare our approach with the state-of-the-art in Table~\ref{sota-kth}. Overall, our method achieves competitive results with recent approaches, while only using skeletal data.

\begin{table}[H]
\centering
\caption{Comparison of our approach with state-of-the-art results for the KTH-Action dataset. *: Deep Learning approach}
\begin{adjustbox}{width=0.41\textwidth}
\begin{tabular}{|l|c|c|c|}
\hline
\textbf{Methods} & \textbf{Input data} & \textbf{Protocol} & \textbf{Accuracy} \\
\hline
Schüldt \etal~\cite{kth} (2004) & RGB & Split & 71.7\% \\
Liu \etal~\cite{LiuLS09} (2009) & RGB & LOAO & 93.8\% \\
Yoon \etal\cite{YoonK10} (2010) & Skeleton & - & 89\% \\
Raptis \& Soatto~\cite{RaptisS10} (2010) & RGB & LOAO & 94.5\% \\
Wang \etal~\cite{WangKSL11} (2011) & RGB & Split & 94.2\% \\
Gilbert \etal~\cite{GilbertIB11} (2011) & RGB & LOAO & 95.7\% \\
Jiang \etal~\cite{JiangLD12} (2012) & RGB & LOAO & 95.77\% \\
Vrigkas \etal~\cite{VrigkasKNK14} (2014) & RGB & LOAO & \textbf{98.3\%} \\
Veeriah \etal~\cite{VeeriahZQ15} (2015)* & RGB & Split & 93.96\% \\
Liu \etal~\cite{LiuSLL16} (2016) & RGB & Split & 95\% \\
Almeida \etal~\cite{AlmeidaPG17} (2017) & RGB & LOAO & 98\%  \\
\hline
Our & Skeleton & LOAO & 96.16\% \\
\hline
\end{tabular}
\end{adjustbox}
\label{sota-kth}
\end{table}

\subsection{UAV-Gesture Dataset}
The UAV-Gesture dataset~\cite{UAV-Gesture} is a 2D videos dataset, consisting of 13 actions corresponding to UAV (\emph{i.e.}, Unmanned Aerial Vehicles) gesture signals. These actions are \textit{All Clear, Have Command, Hover, Land, Landing Direction, Move Ahead, Move Downward, Move To Left, Move To Right, Move Upward, Not Clear, Slow Down} and \textit{Wave Off}. The actions are performed by 11 different subjects in an outdoor scenario with slight camera movements. The dataset contains 119 high-quality clips consisting of 37151 frames. As reported in~\cite{UAV-Gesture}, this dataset is not primarily designed for action recognition, but it can be used for this specific task. The skeletons are available with the dataset and the OpenPose framework was also used to extract them in the COCO format. Table~\ref{results-ours-uav} compares our results with the baseline experiment reported in~\cite{UAV-Gesture}

\begin{table}[H]
\centering
\caption{Comparison of our approach with the baseline on the UAV-Gesture dataset. *: Deep Learning approach}
\begin{adjustbox}{width=0.478\textwidth}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Method} & \textbf{Curve Fitting} & \textbf{Lambda} & \textbf{Alignment Method} & \textbf{Results} \\ \hline
P-CNN~\cite{UAV-Gesture} (2018)* & - & - & - & 91.9\% \\ \hline
Ours & No & - & GAK & 91.6\% \\
Ours & Yes & 10 & GAK & \textbf{92.44\%} \\ \hline
\end{tabular}
\end{adjustbox}
\label{results-ours-uav}
\end{table}

This is a very recent dataset and its principal interest does not rely on action recognition, meaning a lack of results to compare our results with. However, the authors have tested their dataset for the case of action recognition based on skeletons with Pose-Based Convolutional Neural Network (P-CNN) descriptors, that gives us a baseline to compare our results.
The baseline achieves an accuracy of 91.9\% with a Deep Learning based approach, whereas our approach achieves an accuracy of 92.44\%, outperforming the state-of-the-art results when applying curve fitting and the GAK alignment method.

The Figure~\ref{qualitative-uav} presents two very similar actions in the dataset, that is \textit{All Clear} and \textit{Not Clear}. The only big difference between these two actions is the orientation of the hand on the raised arm.
This information is not captured when only retrieving the body skeleton from the OpenPose framework, even if it is possible to retrieve the hand skeleton. Note that with our method, the two actions are only confused three time on a total of 22 sequences, meaning that our method is capable of differentiate minimal changes in the actions.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.47\textwidth]{UAV-Comparison-All-Not-Clear.png}
    \caption{Comparison of two sequences that are confused in UAV-Gesture dataset (top: \textit{All Clear} action, bottom: \textit{Not Clear} action)}
    \label{qualitative-uav}
\end{figure}

\subsection{Computation Time Comparison}
In this analysis, we have computed the time for each step of our approach. Applying the curve fitting algorithm can be resource demanding on some manifolds, as it requires successive computations of Riemannian exponentials and logarithms (see~\cite{Gousenbourger2018} for more information on the computational cost of the method). The alignment method can also be resource demanding, regarding the size of the trajectory. With all these parameters in mind, we propose to compute the time that our method takes to compute specific tasks from pose extraction to action classification. The tests were conducted on a laptop equipped with an Intel Core i7-8750H CPU, 16G of RAM and a NVidia Quadro P1000 GPU. Table~\ref{computation-time-kth} and Table~\ref{computation-time-uav} summarize the execution time of each part of our method for the KTH-Action and UAV-Gesture datasets, respectively. For metric notation, $M_1$ refers to Eq.~\eqref{eq:frob-norm-new-metrix} and $M_2$ refers to Eq.~\eqref{eq:distp2}, in Section~\ref{metrics}.

\begin{table}[H]
\centering
\caption{Execution time (in seconds) obtained on the KTH-Action dataset for the different steps of the method for one sequence.}
\begin{adjustbox}{width=0.478\textwidth}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Pose extraction} & \textbf{Curve Fitting} & \textbf{Alignment method - Metric} & \textbf{Alignment} & \textbf{Classification} \\ \hline
147 & 0.069 & DTW - $M_1$ & 0.034 & 0.41 \\
147 & 0.069 & DTW - $M_2$ & \textbf{0.02} & 0.41 \\ \hline
147 & 0.069 & GAK - $M_1$ & 0.04 & 0.49 \\
147 & 0.069 & GAK - $M_2$ & \textbf{0.019} & 0.49 \\ \hline
\end{tabular}
\end{adjustbox}
\label{computation-time-kth}
\end{table}

\begin{table}[H]
\centering
\caption{Execution time (in seconds) obtained on the UAV-Gesture dataset for the different steps of our method for one sequence.}
\begin{adjustbox}{width=0.478\textwidth}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Pose extraction} & \textbf{Curve Fitting} & \textbf{Alignment method - Metric} & \textbf{Alignment} & \textbf{Classification} \\ \hline
- & 0.504 & DTW - $M_1$ & 0.128 & 0.038 \\
- & 0.504 & DTW - $M_2$ & \textbf{0.072} & 0.038 \\ \hline
- & 0.504 & GAK - $M_1$ & 0.138 & 0.53 \\
- & 0.504 & GAK - $M_2$ & \textbf{0.066} & 0.053 \\ \hline
\end{tabular}
\end{adjustbox}
\label{computation-time-uav}
\end{table}

For the KTH-Action dataset, we consider a sequence of 61 frames and a sequence of 192 frames for UAV-Gesture. First, we can observe that the pose extraction phase takes most of the execution time for the KTH-Action dataset. This is partially due to the fact that our GPU is not powerful enough (we get around 3.5fps with our Quadro P1000). The extraction time is not reported for the UAV-Gesture dataset as the skeletons are available with the dataset. The second thing we can observe is the low difference in computation time for the alignment part when switching from DTW to GAK. We can also note that when we use $M_2$, the computation time can be reduced by a factor 2 compared to the use of $M_1$, showing that the formula~\eqref{eq:distp2} is in our case cheaper to evaluate than~\eqref{eq:frob-norm-new-metrix}. If we only consider the execution time for the treatment of the skeletons, it takes around 0.499 seconds to classify an action of the KTH-Action dataset and around 0.614 seconds for an action of the UAV-Gesture dataset in the best case scenario.

\section{Conclusion and Future Work}\label{sect:conclusion}
In this paper, we have proposed a method for comparing and classifying temporal sequences of 2D/3D landmarks on the positive semi-definite manifold. Our approach involves three different steps: 1) We build a trajectory on the Riemannian manifold from the body skeletons; 2) we apply a curve fitting algorithm on the trajectories to denoise the data points; 3) we perform a temporal alignment using a Global Alignment Kernel. Our experiments on three publicly available datasets show that the proposed approach gives competitive results with respect to state-of-the-art methods.

\section{Acknowledgements}
We thank Dr. A. Kacem for fruitful discussions. We also thank Prof. J-C. Alvarez Paiva  for fruitful discussions on the formulation of the distance between $n \times 2$ landmarks configurations.

\section{Appendix}
\begin{proof} of Theorem~\ref{theo1}.
We can reformulate our metric introduced in Eq.~\eqref{eq:frob-norm-new-metrix} with:
\begin{align*}
d^2(G_i, G_j) & = \Tr \left[ (Z_jQ - Z_i)(Z_jQ - Z_i)^T \right] \\
  & = \Tr(G_i) - 2\Tr(Z_iQ^TZ_j^T) + \Tr(G_j) .
\end{align*}

To minimize our distance, we need to maximize the term $\Tr(Z_iQ^TZ_j^T)$. Let $Z_j^TZ_i$ be a $2 \times 2$ matrix with four unknown values $a, b, c, d$ and let $Q \in \mathcal{O}_p$, we  maximize:
\begin{equation}
\label{eq:maximize-term}
\begin{split}
    \max\;\Tr \left[ \left( \begin{array}{cc} a\cos{\Theta} - b\sin{\Theta} & - \\ - & c\sin{\Theta} + d\cos{\Theta} \end{array} \right) \right] .
\end{split}
\end{equation}

From Eq.~\eqref{eq:maximize-term} we now have to find the maximum of $(a + d)\cos{\Theta} + (c - b)\sin{\Theta}$, meaning that we have to maximize $\sqrt{(a + d)^2 + (c - b)^2}\cos{(O-O')}$. As we want to maximize this value, $O$ has to be equal to $O'$, so $\sqrt{(a + d)^2 + (c - b)^2}\cos{(O-O')} \leqslant \sqrt{(a + d)^2 + (c - b)^2}$. Therefore we can say that:
\begin{equation}
    \max \Tr(Z_iQ^TZ_j^T) = \sqrt{(a + d)^2 + (c - b)^2} .
\end{equation}
\end{proof}

{\small
\bibliographystyle{ieee}


\section{Introduction}

Probabilistic models have been extremely important for successful language processing, from the first accomplishments in natural language processing (NLP), with $n$-grams and their successors \cite[chapters 22 and 23]{RusselNorvigAIModernApproach3rd}, to the tools routinely used today for the annotation of biological sequences in bioinformatics, based on the profile hidden Markov models (pHMM) or covariance models (CM) \citep{durbin_eddy_krogh_mitchison_1998,Coste2016}. Yet, after its success in computer vision and pattern recognition, deep learning is now replacing probabilistic models to become the new key player in the NLP field, showing that neural networks based on dense vector representations are very powerful for performing a large variety of NLP tasks. Impact of deep learning is until now limited in biology, but promising advances are being made (see for instance  \citet{doi:10.1098/rsif.2017.0387}, notably section 3 for the applications to biological sequences) and we can expect progress in NLP to be beneficial for biological sequence analysis once again.\\

While the main appealing feature of neural networks was initially \emph{non-linearity}, the historical ``S''-shape activation functions (such as the $\tanh$ function) are more and more replaced by the Rectified Linear Unit (ReLU) function, reducing non-linearity to its simplest form: a function returning $0$ or a weighted linear combination of its input values. The most important feature of deep learning in NLP seems now to be its ability to learn intermediate representations ---such as the word vector representations (word embeddings) learned by word2vec \citep{word2vec} and GloVe \citep{pennington2014glove}, or the recent deeper contextual representations learned by ELMo \citep{Peters:2018}, ULMFiT \citep{howard2018universal}, and GPT \citep{gpt,gpt2}--- in a hierarchical manner up to the last layers of the neural networks in charge of their weighted combination for solving the task at hand. Using weighted combination of weight vectors (or matrices) optimized by gradient descent appears to be more and more the core of deep learning, as witnessed for instance by the growing importance of tensor-based libraries and of automatic differentiation \citep{DBLP:journals/jmlr/BaydinPRS17}.

Hierarchical representations are common in formal grammars and gradient descent is not an original approach. The key advantage of deep learning with respect to the inference of classical grammatical representations could come from using and combining \emph{weights} rather than \emph{probabilities}. We think that better understanding and evaluating the importance of this difference is essential for future research in this area. As a first contribution, we present here practical problems illustrating the fundamental limitations of current probabilistic modeling approaches for language learning. The question we raise is whether they could be properly handled with more elaborate probabilistic modeling, or whether they demonstrate the interest of the deep "Probxit" towards weighted models.

\section{Modeling with probabilistic grammars}
To cope with inherent "noise", "variations" and "errors", languages have long been modelled with probabilistic models.
In practice, these models rarely exceed the expressiveness of probabilistic context-free grammars, or even probabilistic automata. We introduce briefly the main related definitions and notations following \citet{Vidal:2005:PFM:1070615.1070803,Vidal:2005:PFM:1070615.1070798} to which we refer the reader for a more detailed presentation of the different probabilistic grammatical models of languages.

A probabilistic context-free grammar (PCFG) $G$ is defined as a tuple $\langle Q, \Sigma, S, R, P \rangle$ where $Q$ is a set of non-terminal symbols, $\Sigma$ is a finite alphabet, $S \in Q$ is the initial symbol, $R$ is a set of rules $A \rightarrow \alpha$ with $\alpha \in (Q \cup \Sigma)^*$ and $P : R \mapsto [0,1]$ define the probabilities of rules for each $A\in Q$, under the constraint:
\begin{equation}\label{eq:derivprobdistrib}
  \forall A \in Q \colon \sum_{(A\rightarrow\alpha) \in R} P(A\rightarrow\alpha) = 1.
\end{equation}
A successful (leftmost) derivation of a sequence $x$ from $S$ by a succession of rules $r_1\ldots r_l \in R^l$ is denoted by $S \overset{r_1\ldots r_l}{\Rightarrow} x$ and its probability is $P(r_1\ldots r_l)=\prod_{i=1}^{l} P(r_i)$. If a sequence $x$ can be successfully derived from $S$, its probability is $P(x)=\sum_{r\in R^+\colon S\overset{r}{\Rightarrow} x}P(r)$, otherwise its probability is $0$. Only consistent PCFGs, i.e.\ satisfying:
\begin{equation}\label{eq:seqprobdistrib}
  \sum_{x\in\Sigma^*}P(x)=1
\end{equation}
defining thus a probability distribution on $\Sigma^*$, are classically considered.

Probabilistic regular grammars, more often depicted as probabilistic automata (PA), are PCFGs with simpler rules of the form $A\rightarrow aB$, $A\rightarrow a$, or $A\rightarrow \lambda$  with $A,B \in Q$, $a\in \Sigma$ and $\lambda$ denoting the empty sequence.

\begin{figure*}[t]
  \centering
  \includegraphics[width=0.93\textwidth]{SH3align-to-PA.pdf}
  \caption{pHMM-like probabilistic automata (right) from a toy alignment of 15 SH3 sequences with 7 conserved columns (left). Height of amino-acid letters is proportional to the probability of the transition labeled by the letter. Here, probability of amino-acids in left-to-right transitions is their observed frequency in the modelled column while they are assumed equiprobable in the self-loop transitions introduced at beginning and end of the PA, to parse complete proteins, and in its middle, to allow insertions between columns 4 and 5 (as in protein O45009).}
  \label{fig:SH3aliPA}
\end{figure*}
To give intuition, we will illustrate some limitations of purely probabilistic approaches on a practical task: modeling a family of protein sequences. This task is classically achieved with pHMMs, which are discrete HMMs with a predefined left-to-right topology. Illustration will focus on the core of this approach, simplified and transposed to PA to remain in the grammatical formalism. As an example of a family of protein sequences, we consider the SH3 domain, a region fairly well conserved among several protein sequences that is known to be important for binding and interacting with other proteins. This example has been realistically covered in the pHMM tutorial by \citet{krogh1998introduction}. We will simply use here a toy ``profile PA'' built from an excerpt of the alignment of the SH3 regions from some protein sequences, shown in figure \ref{fig:SH3aliPA}. The choice of this example was obviously driven by our own domain of expertise, but also because it enables us to illustrate modeling problems with respect to ''objective'' physico-chemical considerations. The problems presented in the following sections should still be generic enough to be easily transposed to other applications and methods using probabilistic sequence models.

\section{Limitations}

We present here a list of practical limitations of classical probabilistic language modeling that should not arise when using weighted models. This does not mean that more elaborate probabilistic approaches could not solve these issues, nor that the best way to learn weighted models overcoming these limitations is known. We see these problems as critical testbeds for studying and comparing both approaches in order to better measure the contribution of the shift from probabilities to weights in the success of deep learning.
Ideally, we would like these examples to also help developing principled approaches for learning weighted grammatical models, such as those available for probabilistic models. A pragmatic first step in this direction could be the elaboration of well-founded schemes for parameter estimation of ``profile weighted automata'' outperforming the simple maximum likelihood estimation of pHMMs' probabilities (or, better, the sophisticated Bayesian inference methods used to boost their performances by the addition of pertinent pseudocounts from Dirichlet mixture priors, as initiated for instance by \cite{SjolanderKBHKMH96}).

\paragraph{L1: Probability scattering at each non-deterministic derivation step}

Comparing the probabilities given by probabilistic grammars to sequences of different lengths is a well-know practical problem. The profile PA in figure \ref{fig:SH3aliPA} will assign smaller probabilities to longer proteins even if they have exactly the same SH3 region, because of the loss of probability in self-loop transitions consuming the rest of the protein sequence. This problem is tackled in bioinformatics by looking at the ratio between the probability of the sequence and its probability according to a null model: these ratios are usually directly integrated in the profiles by associating log-odds scores to transitions rather than probabilities. Profile are thus already used under the form of weighted models which combine two probabilistic models.

More generally for PCFGs, equation \ref{eq:derivprobdistrib} implies a probability loss by a factor strictly smaller than one, at each derivation step involving a non-deterministic rule. Such derivations scatter more and more probabilities between alternative suffix languages, resulting in negligible probabilities for longer derivations. Moreover, they induce a dependence of probability values to the number of non-deterministic derivations, which has often to be corrected or compensated in practical tasks, while this could be avoided by not satisfying equation \ref{eq:derivprobdistrib}.

\paragraph{L2: Dependence to number of choices for mass probability distribution within alternative rules}

Another consequence of equation \ref{eq:derivprobdistrib} is that the probability mass has to be shared between all alternative rules. Its amount by rule depends then from the number of alternatives. For the characterization of SH3 regions, figure \ref{fig:SH3aliPA} shows the importance of having the amino-acid {\tt G} at the first position. But when we look at the second conserved column, we can see that it contains only hydrophobic amino-acids and it might be as important for SH3 function to have a hydrophobic residue at the second position as a {\tt G} at the first position. This cannot be modeled with PCFGs using the amino-acid alphabet: assuming that there are 14 hydrophobic acids and that they are equiprobable, the probability of transition with these amino-acids will be at most $\frac{1}{14}$ making a small difference with the completely equiprobable $\frac{1}{20}$ null model, in contrast to $\frac{14}{15}$, or any value close to $1$, that can be assigned at the first position to a transition with {\tt G}. We could imagine working on a larger alphabets, such as the amino-acids powerset, but the exponential growth of the alphabet size would then complicate training. Considering weights (e.g. reflecting adequacy of particular amino-acids at position, instead of occurrence) could be a more efficient solution in practice.

\paragraph{L3: Identical probability mass for all choice points}
Because of the normalisation to $1$ in equation \ref{eq:derivprobdistrib}, all choice points have the same total weight in final probability. Yet, some choices can be more important than others. For instance, prior expert knowledge could tell us that amino acid {\tt P} in the fourth position is more important for the function of SH3 than amino-acid {\tt G} in the first position, but this could not be modelled with a PCFG. Authorizing a normalisation to variable values (e.g. reflecting relevance of position modelled with a particular non-terminal) would help here. Note that it could also be an indirect way to overcome limitation L2.

\section{Towards language modeling Probxit?}

\paragraph{Analysis of limitations}
Limitations seen above are induced by equation \ref{eq:derivprobdistrib} which enables proper PCFGs to satisfy equation \ref{eq:seqprobdistrib} and thus define a probability distribution over all sequences. While discarding the constraint of equation \ref{eq:derivprobdistrib} is tempting, this constraint is intrinsically needed to distribute the remaining probability mass between all sequences that can be generated after each choice point, and it makes perfect sense to model the probability of sequence occurrence in the language with respect to the uncertainty brought by the non-deterministic rules.

Actually, what is illustrated here is more a misused feature of probabilistic grammars rather than a real limitation of these models: probability of sequences under constraint of equation  \ref{eq:derivprobdistrib} specifies the likelihood of their occurrence inside the language, not their probability of being in the language. Considering the SH3 example, probabilities tell us how often we can expect to see a given SH3 sequence comparatively to other SH3 sequences, not to estimate the uncertainty that the sequence can or not perform the SH3 function. If we are interested in membership to language, it is clearly given by the non-probabilistic part of the grammar: a sequence $w$ is in the language $L$ if it can be successfully derived by rules which have non-zero probability. If we were interested in the uncertainty of membership, we should consider learning a probability $P(x \in L)$, with $P(x \in L) + P(x \not\in L) = 1$ instead of equation \ref{eq:seqprobdistrib}.
Likelihood of occurrence in the language can thus be considered at most as an imperfect proxy for membership prediction.

\paragraph{Research directions}
The first research direction to overcome the enumerated limitations is to work on learning the topology of the grammars. Note that this could be completed by learning also a membership probability, but the meaning of this uncertainty, and how it could be learned, would have to be clarified. It could be a measure of the uncertainty of decision with respect to current knowledge (e.g.\ based on the number of available examples supporting the decision). In the example of SH3, we can also imagine, for instance, that the binding affinity of different SH3 regions could influence their probability of binding to another protein and to perform their function: this could be formalized as a kind of membership probability, but would require information about the probability of performing the function in the training set. We will not consider these extensions here and will continue to focus on membership problem.

Learning a good grammar topology can be hard. In contrast, adding probabilities (and a recognition threshold) to simpler topologies has been shown successful to improve sufficiently their expressiveness for many classification task
, while maintaining the number of parameters low enough to be estimable from training sets. This strategy has been shown efficient in practice, despite the limitations that we have identified. A second research direction is to study if more elaborate probabilistic approaches could be used more successfully here. The third research direction is to get rid of the constraint of equation \ref{eq:derivprobdistrib}, which is not really needed for membership prediction, to learn this way weights instead of probabilities of grammar rules.

In this third approach, the weights of sequences could be divided by a partition function (if it exists) to be considered as probabilities taking values in the interval $[0,1]$ and satisfy equation \ref{eq:seqprobdistrib}. Besides the problem of the existence and the computation cost of such a partition function, we think that the constraint of equation \ref{eq:seqprobdistrib} has also to be discarded, since it restricts the expressiveness of the approach. Indeed, (positive) weights are useful if a threshold greater than zero is chosen so as sequences accepted in the language are those whose weight is above this threshold (otherwise, successful derivation by transitions with non-zero weight is sufficient). By transforming weights into probabilities, the choice of a threshold $t$ would then only enable to define finite languages (with at most $\frac{1}{t}$ sequences since total probability mass is $1$), restricting thus unnecessarily the languages that can be represented.

\paragraph{From probabilities to weights}
Getting rid of both equations means shifting from probabilistic to weighted representations. This move is visible in research on learning weighted grammars (see for instance contrastive estimation from \citet{Smith:2005:CET:1219840.1219884}, still motivated by mass distribution constraint) and can be considered an intrinsic feature of deep learning approaches. While practical successes show the interest of these approaches, better understanding their fundamental contribution and being able to develop principled approaches for learning weighted grammars, such as those developed for probabilistic grammars, remains a challenge. A key we identified is not (mis)using occurrence probabilities anymore and we exhibited some advantages of using weighted rather than probabilistic models, but theoretical and practical studies are still needed. The problem of training weighted, instead of probabilistic, profile automata on proteins seems a good support for these studies since it is simple and enables objective evaluation of the approaches. To address these questions, we think that we will need first to find the equivalent for weighted grammars of maximum likelihood or Bayesian approaches (especially for the incorporation of prior knowledge, an essential feature of using pHMMs on proteins) and the nice compositional properties of probabilistic approaches, and would really like to discuss this during the workshop
\clearpage
\section*{Acknowledgments}
I am grateful to Hugo Talibart and Witold Dyrka for fruitful discussions that motivated this work submitted as a four pages opinion paper to ACL 2019 Workshop ``Deep Learning and Formal Languages: Building Bridges''.