{'fraction_non_alphanumeric': 0.047933703271028034, 'fraction_numerical': 0.01796144859813084, 'mean_word_length': 4.0963720930232554, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GeomNet: A Neural Network Based on Riemannian Geometries of SPD Matrix Space and Cholesky Space for 3D Skeleton-Based Interaction Recognition

In this paper, we propose a novel method for representation and classification of two-person interactions from 3D skeleton sequences. The key idea of our approach is to use Gaussian distributions to capture statistics on R n and those on the space of symmetric positive definite (SPD) matrices. The main challenge is how to parametrize those distributions. Towards this end, we develop methods for embedding Gaussian distributions in matrix groups based on the theory of Lie groups and Riemannian symmetric spaces. Our method relies on the Riemannian geometry of the underlying manifolds and has the advantage of encoding high-order statistics from 3D joint positions. We show that the proposed method achieves competitive results in two-person interaction recognition on three benchmarks for 3D human activity understanding.


Introduction
3D skeleton-based action recognition has been an active research topic in recent years with many potential applications. In this work, we focus on 3D skeleton-based twoperson interaction recognition (3DTPIR). Compared to a large number of general-purpose methods for 3D skeletonbased action recognition, methods for 3DTPIR are much less studied in the literature. Recent works [46,63] have shown that state-of-the-art action recognition methods do not always perform well on human interaction recognition. This is because they lack an effective mechanism for capturing intra-person and inter-person joint relationships [46].
In two-person interactions, arm and leg movements are highly correlated. However, these correlations are simply encoded by distances between joints in most existing works for 3DTPIR [17,18,42,65]. This motivates us to use highorder statistics, i.e. covariance matrices to better capture these interactions. It has been known that n × n covariance matrices lie on a special type of Riemannian manifolds, i.e. SPD manifolds (denoted by Sym + n ). A large body of works has been developed for classification of SPD-valued data.
Recently, SPD neural networks have demonstrated impressive results [13]. One of the core issues that remains open is the finding of effective and efficient methods for modeling probability distributions on Sym + n . Since Gaussian distributions (abbreviated as Gaussians) on R n are the most popular probability distributions used in statistics, existing works mainly focused on generalizing them to Sym + n . Such a generalization was first given in [43] in a more general context of Riemannian manifolds. However, the asymptotic formulae of Riemannian Gaussian distributions (abbreviated as Riemannian Gaussians) proposed in this work make them hard to evaluate and apply in practice. Some works aim to address this shortcoming by introducing notions of Riemannian Gaussians in symmetric spaces [49,50] and homogeneous spaces [6]. These have been successfully applied to classification problems. In this work, we also interested in Riemannian Gaussians for classification. However, differently from the above works, we seek methods for embedding Riemannian Gaussians in matrix groups. This allows us to perform classification of Riemannian Gaussians without having to resort to an exact expression of their probability density function as in [6,49,50].
In summary, the main contributions of this work are: • We propose an embedding method for Gaussians by mapping them diffeomorphically to Riemannian symmetric spaces.
• We consider representing a 3D skeleton sequence by a set of SPD matrices that leads us to the study of statistics on Sym + n . We show that the product space of mean and covariance on Sym + n can be viewed as a Lie group with an appropriate group product. Moreover, we point out a connection between this space and the group of lower triangular matrices with positive diagonal entries.
• Based on the theory described above, we introduce a neural network for learning a geometric representation from a 3D skeleton sequence.
• Experiments on three benchmarks for 3D human activity understanding demonstrate the competitiveness of our method with state-of-the-art methods.

Related Works
We will briefly discuss representative works for 3DTPIR (Section 2.1), embeddings of Gaussians (Section 2.2), and probability distributions on Sym + n (Section 2.3).
Approaches for 3DTPIR are much less studied. Handcrafted feature based methods mainly rely on distances [17,18,42,65] or moving similarity [27] between joints of two persons. Li and Leung [23] applied a multiple kernel learning method to an interaction graph constructed from the relative variance of joint relative distances. Two-stream RNNs are proposed in [37,57] where interactions between two persons are modeled by concatenating the 3D coordinates of their corresponding joints, or by augmenting the input sequence with distances between their joints. In [46], Relational Network [51] is extended to automatically infer intra-person and inter-person joint relationships. The recent work [63] deals with graph construction in graph convolutional networks for 3DTPIR.

Embedding of Gaussians
Methods for embedding Gaussians are widely used in statistics, e.g. for measuring the distance between probability distributions. The work of [47] first proposed a distance function based on the Fisher information as a Riemannian metric. However, in the general case of multivariate Gaussians, an exact formula for the distance function is difficult to obtain. In computer vision, one of the most widely used embedding is derived from [34]. The key idea is to identify Gaussians with SPD matrices by parametrizing the space of Gaussians as a Riemannian symmetric space. The work of [5] shares a similar idea of identifying Gaussians with SPD matrices. However, it is based on embedding Gaussians into the Siegel group. In [12], a connection is established between Gaussians and a subspace of affine matrices. The method of [24] relies on the Log-Euclidean metrics [1] for embedding Gaussians in linear spaces.

Probability Distributions on Sym +
n Existing works mainly focused on generalizing Gaussians to Sym + n due to their popularity in statistics. Generalizations of Gaussians are proposed in Riemannian manifolds [43,66], symmetric spaces [49,50], and homogeneous spaces [6]. In [2,4], Riemannian Gaussians are derived from the definition of maximum entropy on exponential families. Family of Alpha-Divergences and other related divergences have also been extensively studied [8].

Statistics on Riemannian Manifolds
The theory presented here is based on [43]. In this framework, the structure of a manifold M is specified by a Riemannian metric. Let x be a point of M as a local reference and T x M be the tangent space at x. A Riemannian metric is a continuous collection of dot products < .|. > x on T x M. The distance between two points of a connected Riemannian manifold is the minimum length among the smooth curves joining them. The curves realizing this minimum for any two points are called geodesics. Let v ∈ T x M be a tangent vector at x. We define the exponential map at x as the function that maps v to the point y ∈ M that is reached after a unit time by the geodesic starting at x with this tangent vector. This map is defined in the whole tangent space T x M but it is generally one-to-one only locally around 0 in the tangent space (i.e. around x in the manifold). Let − → xy = log x (y) be the inverse of the exponential map that is the smallest vector as measured by the Riemannian metric such that y = exp x ( − → xy). The exponential chart at x can be seen as the development of M in the tangent space at a given point along the geodesics.
The definitions of mean and covariance on a Riemannian manifold are given below.
Definition 1 Let x be a random point of probability density function p x . Denote by dist R (y, x) the distance between x, y ∈ M induced by the Riemannian metric of M. The set of expected or mean values is: where dM(z) is the volume measure induced by the Riemannian metric of M, and we assume that the integral is finite for all points y ∈ M (which is true for a density with a compact support).
Definition 2 Let x be a random point,x be a mean value that we assume to be unique. The covariance is defined as: where D(x) is the maximal definition domain for the exponential chart atx.

Proposed Approach
In what follows, we are interested in matrices over the field of real numbers, unless otherwise stated.

Embedding of Gaussians
Let N (n) be the space of n−variate Gaussians, Sym +,1 n be the space of SPD n × n matrices with determinant 1. Lemma 1 gives an embedding of Gaussians in N (n).

Lemma 1
We can identify the Gaussian (Σ Σ Σ, µ µ µ) ∈ N (n) with the following matrix in Sym +,1 n+k : where µ µ µ and Σ Σ Σ are the mean and covariance of the Gaussian, I k is the k × k identity matrix, µ µ µ(k) is a matrix with k identical column vectors µ µ µ.
Proof. See supplementary material. When k = 1, the embedding in (3) becomes the one introduced in [34]. The natural symmetric Riemannian metric resulting from the above embedding is given in Lemma 2.

Lemma 2
The Riemannian metric is given by: where A 1 and A 2 are two tangent vectors at P. Proof. See supplementary material. It turns out that the Riemannian metric given in (4) belongs to the family of affine-invariant metrics proposed in [44]. Consequently, the exponential map at a point can be obtained [45] as: where P is a SPD matrix, A is a tangent vector at P, and exp(.) is the matrix exponential. By inverting the exponential map, we obtain the logarithmic map: where P and Q are two SPD matrices, and log(.) is the matrix logarithm.

Statistics on Sym +
n Suppose that we are given a set of matrices P 1 , P 2 , . . . , P L ∈ Sym + n . From (1), one can define the empirical or discrete mean value of P 1 , P 2 , . . . , P L as: The mean can be computed by an iterative procedure consisting in: (1) projecting the SPD matrices in the tangent space at the current mean; (2) estimating the arithmetic mean in that space; (3) projecting the mean back in Sym + n . These steps are iterated until convergence [39].
Similarly, the empirical covariance of a set of L SPD matrices of mean P m is defined using the discrete version of the expectation operator in (2): where ⊗ denotes the tensor product. We propose to learn a transformation of −−−→ P m P i from T P m Sym + n to another tangent space so that the covariance computed in this space is more discriminative for classification. The transformation is performed by parallel transport (PT). We need Lemma 3 for our transformation.
Lemma 3 Let P, Q ∈ Sym + n . Let the Riemannian metric be the one given in (4). The PT from Q to P along geodesics connecting Q and P of a tangent vector A ∈ T Q Sym + n is given by: Proof. See supplementary material. The formula of PT in (9) is the same as those in [11,55,61] which are all based on the Riemannian metric of the following form: where A 1 and A 2 are two tangent vectors at P.
In [4], the authors also use PT for designing Riemannian batch normalization (RBN) layers. Our method differs from theirs in three main aspects. First, their method learns the parameters of RBN layers from the statistics of minibatches, while our method deals with the statistics within each sequence. Note that a RBN layer can also be designed in our framework and can potentially improve the accuracy of our network. Second, their formulation of Riemannian Gaussians involves only a Riemannian mean without notion of variance. Third, their method does not aim to leverage the second-order statistics (covariance) on SPD manifolds. Now suppose that −−−→ P m P i , i = 1, . . . , L are transported to another tangent space atP (the target point of PT), then the covariance can be estimated as: Let f v (.) be a mapping that vectorizes a symmetric matrix by taking its lower triangular part and applying a √ 2 coefficient on its off-diagonal entries in order to preserve the norm [45]. Then, the covariance is given by: (12) If P c is a matrix of size n × n , then we use a point (P m , P c ) that lies on the product manifold Sym + n ×Sym + n to parametrize the distribution of the given set of SPD matrices. Next, we propose an embedding of this point based on the Lie group theory.

Embedding of Riemannian Gaussians
We first define an appropriate group product on the product manifold Sym + n × Sym + n .
and P c i = L i L T i be the Cholesky decomposition of P c i . Denote by ϕ : Sym + n → M k ×n a smooth bijective mapping with a smooth inverse where M k ×n is a subset of the set of k × n matrices. The group product between two elements of M(n, n ) is defined as: Theorem 1 shows that M(n, n ) forms a Lie group.
Proof. See supplementary material. Based on Theorem 1, we can establish a Lie group isomorphism between M(n, n ) and a subgroup of a group of lower triangular matrices with positive diagonal entries.
Theorem 2 Denote by LT + (n ) the group of lower triangular n ×n matrices with positive diagonal entries, 0 n ×k the n × k matrix with all elements equal to zero. Let where H ∈ LT + (n ), and φ : K + (n + k ) → M(n, n ), φ(K P m ,L ) = (P m , P c ), (15) where P c = LL T , L ∈ LT + (n ). Then φ is a Lie group isomorphism.
Proof. See supplementary material. We now can give the embedding matrix of a point (P m , P c ) ∈ M(n, n ) where P c = LL T as follows: The embedding matrix in (16) depends on the choice of function ϕ(.). In this work, we set where P = UZU T is the eigenvalue decomposition of P, and log(Z) is the diagonal matrix of eigenvalue logarithms.

A Neural Network for 3DTPIR
We are now ready to introduce a neural network (Ge-omNet) for 3DTPIR based on the theory developed in the previous sections. Let N j and N f be the number of joints and that of frames in a given sequence, respectively, Let x in t,i ∈ R 3 , t = 1, . . . , N f , i = 1, . . . , N j be the feature vector (3D coordinates) of joint i at frame t. Two joints i and j are neighbors if they are connected by a bone. Denote by S i the set of neighbors of joint i. Let i 1,r , i 2,r be the two joints selected as the roots of the first and second skeleton, respectively (see Fig. 1). For any two joints i and i r ∈ {i 1,r , i 2,r } that belong to the same skeleton, the distance dist J (i, i r ) between them is defined as the number of bones connecting them (see Fig. 1). The first layer of Ge-omNet is a convolutional layer written as: where x out t,i ∈ R d is the output feature vector of joint i at frame t, andW t ,i is defined as: Here, the set of weights {W u,v }, u, v = 1, 2, 3 completely defines the convolution filters in Eq. (18). Let N j,1 and N j,2 be the numbers of joints belonging to the arms and legs of two skeletons, respectively (see Fig. 1). Let X out,1 and X out,2 respectively of size N j,1 × N f × d and N j,2 × N f × d be the data associated with the arms and legs of two persons. The motivation behind this partition is that the interaction between two persons often involve those among their arms and those among their legs. For each b ∈ {1, 2}, the set of N j,b N f d−dim feature vectors from X out,b is partitioned into L subsets using K-means clustering. Let Gaussian (Σ Σ Σ b l , µ µ µ b l ) whose parameters can be estimated as: Based on the theory developed in Section 4.1, the Gaussian (Σ Σ Σ b l , µ µ µ b l ) can be identified with the following matrix: The above computations can be performed by a layer as: The next layer is designed to compute statistics on SPD manifolds and can be written by: where W b pt , b = 1, 2 are the parameters corresponding to the target points of PT (see Section 4.2). Specifically, P b,m is the mean of P b l , l = 1, . . . , L, and P b,c is given by: The next layer computes the embeddings of statistics P b,m , P b,c , b = 1, 2 and can be written as: where B b is the embedding matrix of (P b,m , P b,c ) given in the right-hand side of (16).
The next layer transforms B b , b = 1, 2 to some matrices in LT + (n + k ) as: where are the parameters that are required to be in LT + (n + k ) so that the outputs D b are also in LT + (n + k ). The network then performs a projection: where Finally, a fullyconnected (FC) layer and a softmax layer are used to obtain class probabilities: where W f c are the parameters of the FC layer, the operator concat(V 1 , V 2 ) concatenates the two column vectors V 1 and V 2 vertically, and C out are the output class probabilities. We use the cross-entropy loss for training GeomNet.

Geometry Aware Constrained Optimization
Some layers of GeomNet rely on the eigenvalue decomposition. To derive the backpropagation updates for these layers, we follow the framework of [16] for computation of the involved partial derivatives. The optimization procedure for the parameters W b pt , W b lw , b = 1, 2 is based on the Adam algorithm in Riemannian manifolds [3]. The Riemannian Adam update rule is given by: where x t and x t+1 are respectively the parameters updated at timesteps t and t + 1, xt is an adaptivity term, g t is the gradient evaluated at timestep t, α, , β 1 , β 2 are constant values. The squared Riemannian norm g t 2 xt =< g t |g t > xt corresponds to the squared gradient value in Riemannian settings. Here, < .|. > xt is the dot product for the Riemannian metric of the manifold in consideration, as discussed in Section 3.1. After updating x t+1 in Eq. (30), we update τ t as the PT of m t along geodesics connecting x t and x t+1 , i.e. τ t = T xt,xt+1 (m t ).
The update rule in Eq. (30) requires the computation of the exponential map and the PT. For SPD manifolds, these operations are given in Eqs. (5) and (9). It remains to define these operations for the update of the parameters W b lw , b = 1, 2. To this aim, we rely on the Riemannian geometry of LT + (n) studied in the recent work [26]. By considering the following metric: where K ∈ LT + (n), U, V ∈ T K LT + (n), U ij is the element on the i th row and j th column of U, Lin has shown [26] that the space LT + (n) (referred to as Cholesky space) equipped with the above metric forms a Riemannian manifold. On this manifold, the exponential map at a point can be computed as: where K ∈ LT + (n), U ∈ T K LT + (n), K is a matrix of the same size as K whose (i, j) element is K ij if i > j and is zero otherwise, D(K) is a diagonal matrix whose (i, i) element is K ii . Also, the PT of a tangent vector U ∈ T K LT + (n) to a tangent vector at H ∈ LT + (n) is given by: where K, H ∈ LT + (n).

Experiments
Our network was implemented with Tensorflow deep learning framework and the experiments were conducted using two NVIDIA GeForce GTX 1080 GPUs. We used GeomStats library [38] for geometric computations. The dimension d of output vectors at the convolutional layer, the number of clusters L, and the learning rate were set to 9, 180, and 10 −2 , respectively. The batch sizes were set respectively to 30 and 256 for the experiments on SBU Interaction dataset and those on NTU datasets. The values of the pair (k, k ) (see (3) and (16)) were set to (2, 3) and (2, 1) for the experiments on SBU Interaction and NTU datasets, respectively. The values of α, , β 1 , and β 2 in the Riemannian Adam algorithm 1 were set to 10 −3 , 10 −8 , 0.9, and 0.999, respectively [20]. In our experiments, GeomNet converged well after 600 epochs. For more details on our experiments, we refer the interested reader to the supplementary material.
NTU RGB+D 60 dataset. This dataset [52] contains 56,880 sequences created from 40 subjects with three cameras views and categorized into 60 classes. We followed the two experimental protocols cross-subject (X-subject) and cross-view (X-view) [52].
NTU RGB+D 120 dataset. This dataset [28] contains 114,480 sequences in 120 action classes, captured by 106 subjects with three cameras views. We followed the two experimental protocols cross-subject (X-subject) and crosssetup (X-setup) [28].

Ablation Study
In this section, we study the impact of different components of GeomNet on its accuracy 2 on SBU Interaction and NTU RGB+D 60 datasets.
Embedding dimensions. Here we investigate the impact of the parameters k and k (see (3) and (16)). Fig. 2 shows the accuracies of GeomNet on SBU Interaction dataset with different settings of (k, k ), i.e. k = 0, 1, 2 and k = 0, . . . , 10. Note that when k = 0, the layer f gaussemb relies only on the covariance information. Also, when k = 0, the outputs B b , b = 1, 2 of the layer f spdstatsemb are simply obtained by the Cholesky decomposition of P b,c , i.e. B b (B b ) T = P b,c . It is interesting to note that Geom-Net achieves the best accuracy with (k, k ) = (2, 3), i.e. none of k and k is equal to 1. This is opposed to previous works [12,25,41,59], where n-variate Gaussians are always identified with SPD (n + 1) × (n + 1) matrices. To the best of our knowledge, this is the first work that shows the benefits of identifying n-variate Gaussians with SPD (n + k) × (n + k) matrices where k > 1. The results also reveal that the setting of (k, k ) has a non-negligible impact on the accuracy of GeomNet. Namely, the performance gap between two settings (k, k ) = (1, 1) (94.54%)   and (k, k ) = (2, 3) (96.33%) is 1.79%. We can also notice that when k is fixed, GeomNet always performs best with k > 1. This shows the effectiveness of our parameterization of Riemannian Gaussians in (16).
To investigate the effectiveness of our proposed embedding of Gaussians outside of our framework, we used it to improve the state-of-the-art neural network on SPD manifolds SPDNet [13]. In [13], the authors performed action recognition experiments by representing each sequence by a joint covariance descriptor. The covariance descriptor is computed from the second order statistics of the 3D coordinates of all body joints in each frame. For SBU Interaction dataset, the size of the covariance matrix is 90×90 (30 body joints in each frame). In our experiment, we combined the covariance matrix and the mean vector using the proposed embedding of Gaussians to represent each sequence. Each sequence is then represented by a SPD (90 + k) × (90 + k) matrix. We used the code of SPDNet 3 published by the authors. Fig. 3 shows the accuracies of SPDNet 4 on SBU Interaction dataset with different settings of k. As can be observed, SPDNet gives the best accuracy with the setting k = 10. The performance gap between two settings k = 1 (90.5%) and k = 10 (92.38%) is 1.88%. The accuracy of SPDNet when using only the covariance (k = 0) is 79.48%, which is significantly worse than its accuracy with the setting k = 10. The results confirm that our proposed embedding of Gaussians is effective in the framework of SPDNet and that it is advantageous over the one of [34]. This suggests that our method could also be beneficial to previous works that rely on Gaussians to capture local feature distribution, e.g. [12,24,25,36,41,59].   Parallel transport. Tab. 1 gives the accuracies of Geom-Net without using PT on SBU Interaction and NTU RGB+D 60 datasets. The accuracies of GeomNet are also shown for comparison purposes. When PT is not used, the covariance in Eq. (12) is computed as: It can be seen that the use of PT is crucial for obtaining high accuracy. Specifically, on NTU RGB+D 60 dataset, computing the covariance without PT results in a loss of 31.44% on X-Subject protocol and a loss of 29.49% on X-View protocol. On SBU Interaction dataset, a significant reduction in accuracy (24.82%) can also be observed when PT is not used. These results highlight the importance of learning the parameters W b pt , b = 1, 2 in GeomNet. Lower triangular matrix learning. Tab. 2 gives the accuracies of GeomNet without using the layer f trilmap on SBU Interaction and NTU RGB+D 60 datasets. Again, the accuracies of GeomNet are also shown for comparison purposes. We can note that the introduction of the layer f trilmap brings performance improvement, i.e. 1.43% on SBU Interaction dataset, and 1.32% on X-Subject protocol and 1.27% on X-View protocol on NTU RGB+D 60 dataset.

Results on SBU Interaction Dataset
Results of GeomNet and state-of-the-art methods on SBU Interaction dataset are given in Tab. 3. For SPDNet, we report its best accuracy using the embedding in (3) with k = 10. We can remark that the accuracies of most of the hand-crafted feature based methods [18,56] are lower than 90%. The state-of-the-art method [9] for skeleton-based action recognition only gives a modest accuracy of 80.35%,

Results on NTU RGB+D 60 Dataset
Tab. 4 shows the results of GeomNet and state-of-theart methods on NTU RGB+D 60 dataset. For ST-GCN and AS-GCN, we used the codes 5,6 published by the authors. For SPDNet, we report its best accuracy using the embedding in (3) with k = 3. We can observe that Geom-Net gives the best results on this dataset. Since ST-GCN is based on fixed skeleton graphs which might miss implicit joint correlations, AS-GCN improves it by learning actional links to capture the latent dependencies between joints. AS-GCN also extends the skeleton graphs to represent structural links. However, AS-GCN does not achieve significant improvements over ST-GCN. This indicates that actional and structural links in AS-GCN are still not able to cope with complex patterns in 3DTPIR. As can be seen, GeomNet outperforms ST-GCN and AS-GCN by large margins. We can also note a large performance gap between GeomNet and SPDNet. This can probably be explained by the fact that: (1) GeomNet aims to learn inter-person joint relationships; (2) GeomNet leverages the covariance information on SPD manifolds.  Table 5: Recognition accuracy (%) of GeomNet and stateof-the-art methods on NTU RGB+D 120 dataset.

Results on NTU RGB+D 120 Dataset
Results of GeomNet and state-of-the-art methods on NTU RGB+D 120 dataset are given in Tab. 5. For SPDNet, we report its best accuracy using the embedding in (3) with k = 3. As can be observed, GeomNet performs best on this dataset. Note that LSTM-IRN performs significantly worse than GeomNet on this most challenging dataset. By adapting the graph structure in ST-GCN to involve connections between two skeletons, ST-GCN-PAM achieves significant improvements. However, ST-GCN-PAM is still outperformed by GeomNet by 3.21% on X-Subject protocol 7 .
The results indicate that: (1) without any prior knowledge, automatic inference of intra-person and inter-person joint relationships is difficult; (2) even with prior knowledge, the state-of-the-art ST-GCN performs worse than GeomNet. Compared to the results on NTU RGB+D 60 dataset, the performance gap between GeomNet and SPDNet is more pronounced on this dataset. Notice that our method is based only on the assumption that the joints of the arms of two persons and those of their legs are highly correlated during their interaction. Therefore, no explicit assumption in pairwise joint connections is required for interaction recognition.

Conclusion
We have presented GeomNet, a neural network based on embeddings of Gaussians and Riemannian Gaussians for 3DTPIR. To improve the accuracy of GeomNet, we have proposed the use of PT and a layer that learns lower triangular matrices with positive diagonal entries. Finally, we have provided experimental results on three benchmarks showing the effectiveness of GeomNet.