Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
Authors: Thibault de Surrel
Fabien Lotte
Sylvain Chevallier
Florian Yger
Venue: ICML 2025 poster (2025)
Keywords: Gaussian distribution
Wrapped distributions
Symmetric positive definite matrices
estimation
classification
Riemannian geometry
density estimation
Primary area: probabilistic_methods
TL;DR: We propose and study the properties and the estimation of wrapped Gaussian distributions on the manifold of SPD matrices

Abstract: Circular and non-flat data distribution are prevalent across diverse domains of data science, yet their specific geometric structures often remain underutilized in machine learning frameworks.
A principled approach to accounting for the underlying geometry of such data is pivotal, particularly when extending statistical models, like the pervasive Gaussian distribution.
In this work, we tackle those issue by focusing on the manifold of symmetric positive definite matrices, a key focus in information geometry.
We introduced a non-isotropic wrapped Gaussian by leveraging the exponential map, we derive theoretical properties of this distribution and propose a maximum likelihood framework for parameter estimation. Furthermore, we reinterpret established classifiers on SPD through a probabilistic lens and introduce new classifiers based on the wrapped Gaussian model.
Experiments on synthetic and real-world datasets demonstrate the robustness and flexibility of this geometry-aware distribution, underscoring its potential to advance manifold-based data analysis.
This work lays the groundwork for extending classical machine learning and statistical methods to more complex and structured data.

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
Thibault de Surrel Fabien Lotte?
Sylvain Chevallier Florian Yger
Abstract
solution by leveraging the exponential map to wrap proba- bility distributions defined in the Euclidean tangent space onto the manifold itself. This approach yields wrapped dis- tribution, intrinsically aligned with the manifold's s geometry. Wrapped distributions have found applications across many domains, such as embedding single-cell RNA data (Ding & Regev, 2021), analyzing wave patterns (Jona-Lasinio et al., 2012), recognizing video and image features (Turaga et al., 2011), modeling Gaussian processes manifolds (Mal- lasto & Feragen, 2018 88 Liu et al., 2024) or doing statistics measurements of orientations
Circular and non-flat data distributions are preva- lent across diverse domains of data science, yet their specific geometric structures often remain underutilized in machine learning frameworks. A principled approach to accounting for the underly- ing geometry of such data is pivotal, particularly when extending statistical models, like the per- vasive Gaussian distribution. In this work, we tackle those issue by focusing on the manifold of symmetric positive definite (SPD) matrices, a key focus in information geometry. We introduce non-isotropic wrapped Gaussian by leveraging the exponential map, we derive theoretical prop- erties of this distribution and propose maximum likelihood framework for parameter estimation. Furthermore, we reinterpret established classifiers on SPD through probabilistic lens and introduce new classifiers based on the wrapped Gaussian model. Experiments on synthetic and real-world datasets demonstrate the robustness and flexibility of this geometry-aware distribution, underscor- ing its potential to advance manifold-based data analysis. This work lays the groundwork for ex- tending classical machine learning and statistical methods to more complex and structured data.
2024).
In this paper, we focus on wrapping Gaussian distributions on the Riemannian manifold of Symmetric Positive Definite (SPD) matrices equipped with the Affine Invariant Rieman- nian Metric (Pennec, 2020). Gaussian distributions are a comnessonn of machine learning and statistics (see p.102 of Casella & Berger 2001) due to their ubiquity and theoretical underpinnings, such as the Central Limit Theorem (CLT) (Section 5.4 of Wasserman 2004), which ensures that Gaus- sian distributions naturally arise in many scenarios. We will see that we are able to extend the CLT to wrapped Gaus- sians, providing a theoretical justification for their study, SPD matrices, which form the Riemannian manifold Pa. are pivotal in numerous applications, including Diffusion Tensor Imaging (Pennec, 2020), Brain-Computer Interfaces (BCI) (Lotte et al., 2018), 8)) process control (Willjuice Irutha- yarajan & Baskar, 2010), and video processing (Tuzel et al., 2008). The inherent Riemannian geometry of Pa necessi- tates adopting methods that respect its manifold structure when analyzing SPD data.
1. Introduction
When dealing with complex data, modeling them as lying on a manifold often provides a powerful solution (Sanborn et al., 2024; Jo & Hwang, 2024). However, classical Eu- clidean probability distributions fail to capture the intrinsic geometry of the underlying manifold. This limitation ne- cessitates adapting the choice of probability distributions to respect the manifold's s structure. In this work, we propose a
This paper is organized as follows: in Section 3, we in- troduce the Riemannian geometry of Pa. Next, we define wrapped Gaussians on Pa and explore their theoretical prop- erties in Section 4. In Section 5, we develop a Maximum Likelihood Estimator for parameter estimation and validate it with synthetic experiments. Building on this founda- tion, Section 6 revisits existing classifiers on Pa through probabilistic lens and introduces novel classifiers based on wrapped Gaussians. Finally, we evaluate these classifiers with synthetic and real-world datasets in Section 6.4.
'LAMSADE, CNRS, PSL Univ. Paris- Dauphine, France Inria center at the University of Bordeaux / LaBRI, France 'TAU, LISN, University Paris-Saclay, France 'LITIS, INSA Rouen-Normandy, France. Correspondence to: Thibault de Surrel <thibault.de- surrel@ Lamsade daw.hiie..d.pii...f.......
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
2. Related works
framework that unifies and generalizes classification on SPD matrices, and propose new classifiers. This application shows the potential of our wrapped Gaussian to become generic, flexible and powerful tool for manifold- based data analysis.
Other works have already tried to extend the Gaussian distri- bution a Riemannian manifold. Said et al. (2018) propose an isotropic Gaussian on Riemannian Symmetric Space M defined using center of mass ÿ E .M and scaling fac- tor o > 0. In our work, we are looking for more complex model, requiring a non-isotropic distribution with some pre- ferred directions. A non-isotropic Gaussian distribution on a manifold has been proposed in Pennec (2006), in which the authors use the characterization of the Gaussian distribution as the distribution that maximizes entropy given mean and covariance matrix (theorem 13.2.2 of Kagan et al. 1973). However, sampling this distribution leads to computational difficulties, the normalizing constant cannot be computed explicitly and in the case of a full covariance matrix, the estimator of the parameters becomes problematic.
3. How to deal with SPD matrices ?
3.1. The Riemannian geometry of SPD matrices
The set of d d symmetric, positive definite (SPD) matrices, denoted Pa is defined as follows:
Pa = {p E IRdxd p² = and Vo E Rd (0), x pr > 0}.
This set is convex and open in the set of d x d symmetric matrices Sq and thus, it is a manifold of dimension d(d+1)/2 For all p Pa, the tangent space TyPa at p can be identified with Sa. Moreover, for p E Pa one can define an inner product on the tangent space T.Pa at by:
Wrapped distributions have first been studied in directional statisties (Mardia & Jupp. 2000), on a circle (Collett & a sphere (Hauberg, 2018). Wrapped Gaussians have also been instantiated on hyperbolic spaces, first by Nagano et al. (2019) and then by Mathieu et al. (2019) and Cho et al. (2022). They mainly use it as the distribution of the latent space of a Variational Autoencoder which trained to learn the distribution. Apart from the manifold, another difference with our approach is that they wrap the Gaussian using a composition of the exponential map with parallel transport where we will only use the expo- nential map. Wrapped distributions have also been studied on more general classes of Riemannian manifolds. For example, Gaalza-Grriia et al. (2022) define wrapped distri- butions on homogeneous Riemannian manifolds. A major difference with our work is that they use a volume preserv- ing map to fush-forward the density from the tangent space to the manifold, leading to simpler expression of the den- sity. without any volume change term. In Chevallier et al. (2022) and Chevallier& & Guigui (2020), the authors work on general symmetric spaces and mainly study the Jacobian determinant of the exponential map, first a general setting and then on different examples (Grassmannians,, pseudohy- perboloids and special Euclidean group). Unlike our work, they consider the distribution the tangent space to always be centered, where we consider more general setting by al- lowing/ μA f O. To estimate the parameters of their wrapped Gaussians, they use moment estimation. Finally, in Troshin & Niculae (2023), they define a more general wrapped Gaus- sian, the 3 Gaussian that has compact support.
Lewis, 1981)
Vu, U E TpPa, (u, u), tr(p- up)
(1)
This inner product is called the Affine Invariant Rieman- nian Metric (AIRM) (Pennec, 2020) as, if a is an invert- ible matrix, one has (ama", T, .ava"' ) Yapa" = {u, u)p- Once endowed with this metric, Pa is a complete connected Riemannian manifold of non-positive curvature (see Ap- pendix I of Criseitiello & Boumal 2023). It is there- fore a Hadamard manifold (Shiga, 1984). Using the Car- tan-Hadamard theorem (theorem 12.8 of Lee 2018) one has that Pa is diffeomorphic to IR Rd through the ex- ponential map Expp: TpPa ~ R 1d + Pa. Another consequence of the completeness of Pa is that each pair of points (p, 4) E P2 can be connected by a unique minimiz- ing geodesic whose length defines distance on Pa. This AIRM distance, between p and q is given by:
8(p. g) = log( 11////111///////
(2)
where log the matrix logarithm and || ||r the Frobe- nius norm. Other useful tools of Riemannian geometry that will be used in the following are the exponential map Expp: TpPa - Pa and its inverse, the logarithm map Logp: : Pa + TpPa. For p, q E Pa and u E TpPa. those maps are given by the following expressions (see chapter 6 of (Bhatia, 2007)):
(xp,(u) = 122 exp(p 1/2 up 1/3)p//2 Log, (q) = 1/2 log(p 1/2 up
(3)
1/2ppp1ppppppppp 11///pp1/221p1 1//2
In the following, we propose a wrapped Gaussian on the manifold of SPD matrices that is not centered on the tangent space. After deriving some theoretical properties, we show that our wrapped Gaussian can be used in practice, showcas- ing the estimation of the parameters from a finite number of samples. Finally, we use our wrapped Gaussian to build
Finally, when Pa is equipped with the AIRM metric given at Equation (1), the Riemannian volume element at p [pn]]]] Pa is given by (see Scction 4.1.3 of Terras 1988):
dvol(p) = dct(p) - (d+1)/2
dp
(4)
2

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
(M, Σ)
Algorithm Sampling from Wrapped Gaussian
0g
WG μ. 2)
Require: p E Pa. Raddiiy DE Pataaa/a MM,, 2)
TpPa
1: Sample
Logp
ExPp
2: Compute X + Exp, (Vect, (t))
3: Return X ~ (,,
*p
WG(p:H,
Pd
given wrapped Gaussian can be parametrized by different sets of parameters, leading to an equivalence relation. the following, we will denote θ = Pa x 12dd++1// Paaan1/2 the space of parameters of the wrapped Gaussian which is product manifold.
Figure 1. An illustration of a wrapped Gaussian WG p4, E).
where dpi is the Lebesgue measure on R. The volume dvol is also invariant under congruence by invertible matrices. Remark 3.1. Other metrics could be used on Pd. such as the Log-Euclidean metric (Arsigny et al., 2005) or the s- divergence metric (Sra, 2015). However, the AIRM metric widely used in the literature and has many properties that make it suitable for our work (such as the affine invariance). For more details on the different metrics on Pa, we refer the reader to Chevallier et al. (2021) and Thanwerdas (2022).
4.1. The definition
To define a wrapped Gaussian, start with classical Euclidean Gaussian random variable in R 1g<4+1/2 and push it on the manifold Pa using the exponential map as follows:
Definition 4.1 (Wrapped Gaussian). Let p E Pa. P E R Rdd=+1//2 and Σ E Pararaya A random variable X on Pa follows a wrapped Gaussian denoted WG(p; 2) if
X = (Vect, (t)), to (NH,
3.2. The vectorization of the tangent spaces
A wrapped Gaussian is illustrated in Figure 1. This gives us very simple algorithm to sample a wrapped Gaussian WG μ4, E), since it is simply based on the sampling of (H, Σ) in TaPa and the computation of Vect, and Exp, that has closed-form formulas. We provide a sampling al- gorithm at Algorithm 1. Moreover, definition using push-forward (see definition 2.1 of Peyré & Cuturi (2020), Section 3.6 of Bogachev (2007) or Ap- pendix B):
As described in Section 3.1, the tangent space Ta Pa at point p Pa is identified with the space of d d symmetric matrices SJ which is isomorphic ddd1///,, We define such isomorphism called the vectorization as follows: Definition 3.2 (Vectorization). We start by defining the vectorization at identity for symmetric matrix [[,,]]]:
can rewrite this
Vecti,: : u Ti,Pa (W11, 22112, U22, V2u1a, V2u23, U33,
V2 uda) C 1y/
WG(P:H, = (Exp, o (oott"))N((4 Σ)
(5)
Then, for p E Pa, we define the vectorization at p:
where # denotes the push-forward operator.
Vectp: u E TpPa + Vectr, (1-1/2 up2222
1/2
Let us comment on the different parameters: p gives us a tan- gent plan from which the Gaussian is wrapped, so it locates the Gaussian on the manifold. yu is the mean of the Gaussian in the tangent space T,Pa. As Pa and TaPaa are in bijection through the exponential map. p and μ play a symmetric role that will be unveiled in Section 4.4. Finally, Σ is the covari- ance matrix between the entries of the SPD matrices. In the special case where the SPD matrices are covariance matri- ces, Σ models the covariance of the covariances, therefore, it can be seen as a fourth order moment.
One of the important property of Vect, that it is an isom- etry between (T,Pa, (p) and (Ra(d+1)/2, (2)) More information this vectorization can be found in Section 3.3.3.3. of Pennec (2020) in Appendix A.
4. Wrapped Gaussian on the manifold of SPD matrices
In this work, our objective is to define a non-isotropic Gaus- sian on the manifold of SDP matrices Pa. In this section, we define a wrapped Gaussian through the way it will be sam- pled. We will also give the density of a wrapped Gaussian and, we will show that, unlike usual probability distributions,
4.2. The density of a wrapped Gaussian
Now that we have defined the wrapped Gaussian, we give its density using the push- push-forward definition of WG μ. E) given in Equation (5) (see Theorem B.2 in Appendix B)
3

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
Theorem 4.2. The density E of the wrapped Gaussian
One can find a proof of this result, as well as an explicit expression for Va in Appendix D.
WG nb, E) exists and is:
Vx E Pd. fpux(()= Iu. (Legg
Then, we give a wrapped version of the multivariate Central Limit Theorem (CLT) (Theorem 5.12 of Wasserman 2004) for the manifold Pa. For this, we define the logarithmic product introduced in Arsigny et al. (2006):
(6)
11 (Logg (x))]
where 9, is the density of the multivariate Gaussian N(, Σ) and Jp(+) = det(d Ex),()) the Jacobian deter- minant of the exponential map Expp-
Definition 4.6 (Logarithmic product). Let p, E P.. The logarithmic product of p and q is defined as:
The Jacobian determinant of the exponential map describes the volume change induced by the identification of Pu and Pa by the exponential map Expp. To compute the density explicitly, one needs to compute the Jacobian Ju(), which has a closed form expression for the manifold Pa (see Section 5.3 of Chevallier Guigui (2020) or Appendix C).
p g exp (logp + log q)
Equipped with this logarithmic product, (Pa. )) forms a commutative group that is isomorphic to (Sa. +)- One can generalize this notation to the sum of SPD matrices Pl.-..: Pn Oi Pi- = P @... Pn- Using this logarith- mic product, we can state the following theorem:
Proposition 4.3. The Jacobian determinant of the exponen- tial map at the identity ExP1, is:
Vue E Pd. (u)(u) = 2// sinh Ac(())-,((u)
Theorem 4.7 (Wrapped CLT). Let (X;) be sequence of i.i.d. random variables on Pa- We suppose that the sequence (Vectz, (Togı, (X:)))en- of random variables Rd(d+1)/2 admits finite second order moment. We denote by l the mean and by E the covariance matrix of V (Log1, (X1)). Then,
2
(uuu) Az (
where the (u) are the eigenvalues of u. Then, one can use the previous formula to compute the Jacobian determinant of the exponential map at any point p Pa:
Vu E TaPa, Jp(u)= = JI. (p-2 1/22).
n O o m)=)) 3/2
WG(1a;
It should be noted that, unlike the wrapped Gaussians de- fined in Galaz-Garcia et al. (2022), Troshin & Niculae (2023) and Nagano et al. (2019), we do not restrict our- selves to a centered multivariate Gaussian (0, Σ) on the tangent space TaPa. In our case, we allow the wrapped Gaussian have an extra parameter μ. thus extending the flexibility and applicability of the model. Having non- centered distribution the tangent space Pa leads to new considerations that will be discussed in Section 4.4. Remark 4.4. In this work, we focus on extending the mul- tivariate Gaussian to a Riemannian setting. With no extra difficulty, one could wrap an Elliptically Contoured distri- bution (EC) on Pa (chapter 6 of Johnson (1987) or Delmas et al. 2024). We give more detail on this in Appendix F.
2-000
i=1
where d
denotes the convergence in distribution and -0
where m Expr, Expr, (ffctt,(()).
This theorem shows the interest of wrapped Gaussians, as they naturally appear in the limit of a logarithmic product of random SPD matrices. The proof of this theorem can be found in Appendix E. One can note that by generalizing the logarithmic product defined in Definition 4.6 to another tan- gent space T Pa. one can extend the limit to WG(p; 0, E).
We can also give information on the mean of WG μa, E) in the special case of μ. = 0 i.e. when the distribution on the tangent space TpPa is centered. We recall from definition 3 of Pennec (2019) that a mean, or exponential barycenter, of a probability distribution a on Pa is defined as a point P E Pa satisfying JP, (z))(=)dd(22))
4.3. Some properties of wrapped Gaussians
Let us now give some properties of the wrapped Gaussians. We start by rescaling property.
O. For wrapped Gaussians, one has the following result:
Proposition 4.5. Let X ~ WG Oaie+1/2s La())s)))) and let (p. U, Σ) E O. There exists a transformation of X de- noted V such that
Proposition 4.8. A mean of WG 0, E) is
The proof is straightforward from the definition of the wrapped Gaussian WG(p; E).
v(X) WG(55,.))
Thus, the wrapped Gaussian WG Ia; O2s Ia) is a building block of the wrapped Gaussians.
More information on the properties of this logarithmic product can be found in Arsigny et al. (2006).

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
We denote by [0x] the equivalence class of 0a:
Nuu - D)
Pa
[0a] = (0''0',,(,)) [0
Pa
0,
N(pEE
Using Proposition 4.9, one has the immediate corollary:
Corollary 4.11. Let la = (Pa.Pa Pa, Ea) E ④ and 0g (pa, Ma, Σx) E O. If there exists t E R such that pa e'Pas P2 = Pa + t Vectp, (Pa) and Σg = Ea. then Pa = 0g.
Expetp
ExPp
c'y
WG Σ)
(rr d'p
Remark 4.12. All equivalence classes do not contain a wrapped Gaussian that is centered on the tangent space Let us consider a wrapped Gaussian WG(p; μ. ))) Then, the equivalent wrapped Gaussians are of the form WG(e'p; p; μ Lv, Σ) for E R. If μ and v are aligned i.e. there exists / E R such that f fv, then the equivalence class contains wrapped Gaussian with μ 0. However, when μl and V = (1, ample, take = v + (1..... 0) = (2, 1, ,,,0, O)), then there exists no t E IR such that μ = tv and the equiva- Ience class does not contain a wrapped Gaussian with jr = O. Therefore, allowing μ * increases the expressiveness of the model.
Pd
Figure 2. Illustration of the equivalence between two wrapped
Gaussians given by Proposition 4.9.
1, 0, .0) are not aligned (for ex-
4.4. An equivalence relation between wrapped Gaussians
A key property of the Euclidean Gaussian (H, Σ) is its identifiability, meaning that the map (N,, Σ) / (n, Σ) is one-to-one (definition 1.5.2 of Lehmann & Casella 1998). The notion of identifiability is important when one wants to learn the parameters of a distribution from a finite set of samples as we will do in Section 5. If we consider a model based on the wrapped Gaussian, we lose this property as we have the following proposition, illustrated in Figure 2:
Once we have defined an equivalence class [0] of parameters that define the same wrapped Gaussian, it is natural to define a representative of [0). We define it as follows:
Definition 4.13 (Representative of an equivalence class). We choose as representative of the class [0], the tuple of pa- rameters gmin = (p"iin, ,1 ,umin, Emin) such that μ m is minimal in the sens of || 1|2. We call it the minimal representative.
Proposition 4.9. Let (p, Σ) E and E R. One has that WG(p; p. Σ) and WG(ee'r p3 Lu, Σ) are equal where (ec بلاا ,,.... Ruddiny 1//2
One is able, given tuple of parameters 0 = (p, μ, 2), to compute the minimal representative of the class [0] of equiv- alent tuples of parameters using the following proposition:
Lnnes
d(d-1)/2 zevos
All the proofs of this section can be found Appendix G. One can verify using Equation (3), that t + e'p is the geodesic y: & + Exp, (fp) in Pa starting at point p and with initial velocity p (as p symmetrie matrix, it also belongs to TpPPa Sa). Moreover, the map L + J Iv j t (p) is also the geodesic in ddd+1) with initial point and initial velocity -v = (p) Therefore, when p is pushed in one "direction" " (initial velocity p). μ is pushed in the opposite "direction" (initial velocity - (p).
Proposition 4.14. Let 0 (p. H,, Σ) E ④ be parameters. Then, the minimal representative of the class [0] as defined at Definition 4. 1.3 is 0 (pmin, mimin E.iin) where
nain ed Eim Hip,
min 1
Σ 1=1
Smin E
where we recall that v = (1, 1,0, ..... E Rid++1//2
A wrapped Gaussian can thus be represented by several sets define an equivalence relation between sets of parameters that define the same wrapped Gaussian:
This minimal representative will be used in the following.
of parameters,
5. Estimation of the parameters of a wrapped Gaussian distribution
Definition 4.10. Let 0,, = (Pa, Ea) E θ and 03- (ps. Ap, -n)) E θ be two sets of parameters. Then, Oa, and f are equivalent, which we denote by 0_, 83. they define the same wrapped Gaussian i.c.
In this section, we tackle the parameter estimation problem of a wrapped Gaussian given samples. After introducing the Maximum Likelihood Estimator we will use, we lead some synthetic experiments to assess its performance.
WG(pai Pa; E.) = WG M43 Eg).
5

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
Distance to the real p
Distance the real Pe Dimensioo
Distance to the real X
Dimension
Dimension
C1
10
30-1
3:30
Ntaboo
Nazioo potzza
Naaterr potzts
Figure 3. Results of the synthetic experiment on the estimation of parameters of a wrapped Caussian using an MLE.
5.1. A Maximum Likelihood Estimator
gate Gradient algorithm (Boumal, 2023) on the prod- uct manifold o Pa x Rd(d+1)/2 x Paiataya to com- pute the optimal parameters (ĀN, juN, Σn). We imple- mented this MLE in Python using the toolbox Pymanopt (Townsend et al., 2016). The codes for the different experiments is available at https://github. com/ htwp:///dsuurree//pppeepeeeeeeee
Let .i...., NN be N SPD matrices of size d d indepen- dently sampled from a wrapped Gaussian with unknown parameters 6* = (p", ,E) supposed to be minimal in the sens of Definition 4. 13. Our goal is to estimate 0* given the samples (II, .... IN using a Maximum Likelihood Esti- mator (MLE) (see Section 9.3 of Wasserman 2004). Let us introduce the likelihood CN of the model:
This estimation problem can become challenging as the dimension d of the considered SPD matrices increases. In fact, the number of coefficients to estimate is:
CNlPi Σ)
II mmi)((s))
d(d + 1)
2 d(d + d2 (d + 1)2 + 2ddd +1)
2
where Fp is the density of the wrapped Gaussian WG(p; ne, Σ) as given in Equation (6). We will also consider the log-likelihood yn defined as (x(p: p. Σ) = log (p: p, Σ). Then, we define the classical MLE 0,- = (fv. ΣN) ) as the parameter that maximizes LN (or equivalently n)). In the Euclidean setting, one has closed form of the MLE of μ and Σ, obtained by computing the gra- dient of CN. When dealing with wrapped Gaussians on Pa. were not able to derive closed form for the MLE pvv of the parameter p. Moreover, it is unlikely that there exists such a closed form, as for example, there is no closed form for the Riemannian mean on Pa (Moakher, 2005), Never- theless, the MLE of μ and Σ are analogous to the Euclidean setting but depend on p":
P
which grows at a rate of O(d)). For example, if d 10. there are 1, 650 coefficients to estimate and if dl = 30, the number jumps to 109, 275. One should thus make sure that the number of samples N is sufficiently big for results of the MLE to make sense. To get around this issue when the number of samples is small, one can assume that the co- variance matrix is diagonal, which reduces the number of coefficients to (dd)). This assumption implies independent entries of the SPD matrices and will be used in Section 6.4.
In other works, such as in Chevallier et al. (2022) in Chevallier & Guigui (2020), the authors use the method of moments (see Section 9.2 of Wasserman 2004) to estimate the parameter. In our case, we can use the method of mo- ment only when we know a priori that μ = Then, as given in Proposition 4.8, mean of the wrapped Gaussian is p", therefore, it can be estimated using the Riemannian mean PN = 6 ..... 2N) (Moakher, 2005), and Σ can be estimated using Proposition 5.1. However, in more general case of μ² f O, estimating p" using the Riemannian mean does not lead to the correct estimation of the true parameters. We give more details on why in Appendix H.
Proposition 5.1. The MLE jun and ΣN of the parameters I and E of the wrapped Gaussian are:
N
AN =
VLogp= (I;),
N
N
ΣN= =
(VLogp: (x:) - ) (VLog, (x:) - (uN)"
N
where VL.ogp = Vectp o Logg..
5.2. Synthetic experiments
One can note that in proposition 4.7 of Galaz-Garcia et al. (2022), they also have a closed from for the MLE of Σ that depends on p" without any closed form for Pv. In practice, we used a Riemannian Conju-
We led some synthetic experiments to evaluate the MLE's performances. For this, we sampled N points from a wrapped Gaussian in Pa whose minimal parameters
6

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
,p',"" ,,,x,,)) are known. Here, Σx is always chosen to be full SPD matrix. Then, we optimize the MLE find (DN,UN, nx), and compare the true and minimal estimated parameters: for and Σ, we compare them using the AIRM distance and for jr, we use || ||2- In this experiment, we looked at how the estimation error evolves when the number of samples N grows from 100 to 10, 000. We also compare different dimensions d E {2, 5, 103. More details on the experimental setup are given in Appendix I. The results of this experiment can be found in Figure 3. We can see that, as one would expect, as the number N of points sampled grows, the estimation error decreases. Moreover, we can remark that the dimension d does not affect the results of the estimation of p and Jr, but really affects the estimation of Σ. The higher the dimension d, the higher the error 8(Σ,, 2**) is. This is coherent as the number of parameters of Σ grows as Odd)), ), so even small increase of the dimension d leads to an important increment in the estimation error of Σ. We also led some experiments in the case where the covariance matrix Σ is diagonal. We show, in Appendix J, that in this case, one needs fewer samples to have a good estimation of E when the dimension d rises.
class is modeled by Mu Ex), ), the decision boundaries are quadratic. This classifier is called Quadratic Discriminant Analysis (QDA). One can restrict the covariance matrices to be diagonal, which leads to the Diagonal LDA and Diagonal QDA classifiers (Dudoit et al., 2002). Then, the Diagonal QDA classifier is equivalent to the Gaussian Naive Bayes classifier (see chapter 8 sec 3.3 of Bishop 2007).
A possible extension of LDA (or QDA) to the manifold of SPD matrices Pa is call Tangent Space LDA (or Tangent Spuce QDA) and is described in part IV B. of Barachant et al. (2012). The Riemannian mean T of the training set is computed, and all training points are sent to the tangent space Te Pa via the logarithm map Loge- Then, a classical I.DA (or QDA) can be used this Euclidean space.
Other classifiers from SPD matrices Other classifiers that have been developed for SPD matrices. For example, Multinomial Logistics Regression has been extended for SPD matrices in Chen al. (2024). They rely on metric that are pulled back from the Euclidean space which is not the case of the AIRM metric use in our work. Several deep learning approaches have been proposed to classify SPD matrices (Huang & Van Gool, 2017; Brooks et al., 2019; Nguyen, 2021). However, most of these approaches distort the geometry of the manifold and are out of the scope of this work, as our approach does not rely on deep learning.
6. Classification using wrapped Gaussians
In this section, we demonstrate that widely used geometry- aware classifiers on Pa can be integrated into a probabilistic framework. We also introduce new classifiers based on the wrapped Gaussian and conduct experiments on real-world data from various applications.
6.2. A general probabilistic framework
Our goal is to show that the MDM, Tangent Space LDA and Tangent Space QDA can seen as part of a proba- bilistic framework on the manifold of SPD matrices. More precisely, we will show that the previous classifiers can be rewritten as Maximum Likelihood based classifiers (like the classical LDA or QDA) where the different classes are mod- eled using distributions on the manifold Pa. Let us consider K classes of labeled SPD matrices and let us denote ag the modeled distribution of class k:
6.1. Classifiers used for SPD matrices
MDM The Minimum Distance to Mean (MDM) algorithm described in Barachant et al. (2010) a popular classifier for SPD matrices. Given training set of labeled SPD matrices, the MDM computes the Riemannian mean (Moakher, 2005) G₂ of each class i E f1...., K}. Then, given a new SPD matrix p. the predicted class ie is the class for which the
distance between p and G₂ is the smallest.
MDM For the MDM classifier, we first need to recall the isotropic Gaussians on Pa introduced in Said et al. (2018). Let ÿ E Pa and a > O, then, the isotropic Gaussian denoted G 0)) is defined by the following density:
LDA and QDA In Euclidean setting, the Linear Dis- criminant Analysis (LDA) (section 4.3 of Hastie et al. 2009) is classifier that assumes that each class k E {l....,K} is modeled by a multivariate Gaussian N N(uu, Σ) where the covariance matrix Σ is shared among all the classes. First, the parameters of each class are learned using an MLE on the training data. Then, to classify a new point 2, LDA compares the ligeliill of according each class and chooses the class that has the highest log- likelioel As the covariance matrix is shared among the classes, the decision boundaries are linear, which led to its name: linear discriminant analysis. When one assumes that the covari- ance matrices are not equal among the classes, i.e. each
1 ((
8(y. 7)²
Vy E Pa. f3.0())
exp
20
where ((o) is a normalizing constant. This normalizing constant depends only on a, and not ÿ shown in propo- sition 1 of Said et al. (2018). If one supposes that each class is modeled by an isotropic Gaussian with shared o among all the classes i.e. ag = G( σ) then, the MDM is equivalent to a maximum likelihood classifier. Here, the spread a of the isotropic Gaussians does not play any role
7

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
Dataset
Application
Dimension
matrices
classes
Reference
BNC12014066 4004
BCI BCT
3 x 3
720 x 9 subjects
(Leeb et al., 2007)
Zhou2016
5 5
320 4 subjects
(Zhou al., 2016)
AirQuality
Atmospheric data
102
(Smith et al., 2022)
Indiana Pines
Hyperspectral imaging
5
14,641 641
12
(Raumgardner al., 2015)
Pavia Univ.
Hyperspectral imaging
5 5
185, 176
Salinas
Hyperspectral imaging
5.x 5
94, 94,184 184
17
Textile
Image Analysis
10 x 10
16, 000
(Bergmann et al, 2021)
BreizhCrops
Multispectral imaging
13 13
178,658
(RuBwurm et al.. 2020)
Table 1. Summary of the datasets used for the experiments.
Dataset
||
Acc. MDM
Acc. TS-LDA
Acc. TS-QDA
Acc. Ho- WD-
Acc. Hc-WDA
BNC12014004
78.71 (=14.53)
78.73 (+14.52)
76.07 (+13.94)
75.38 (+14.40)
74.37 (+14.73)
Zhou2016
91.18 (15.51)
91.21 (15.50)
89.45 (17.43)
85.92 (19.20)
82.86 (=11.65
Air Quality
94.05 (16.53)
94.05 (+6.53)
97.05 (=4.45)
96.05 (#4.17)
97.00 (+4.47)
Indiana
58,01 (=0.72)
67.07 (+0.53)
73.38 (=0.45)
73.74 (+0.51)
74.30 (+0.83)
Pavia Uni.
72.32 (+00.9))
81.61 (4-0.6)
87.16 (+0.09)
87.36 (+00.7)
85.54 (+0.15)
Salinas
36.42 (+0.12)
46.30 (40017)
69.87 (+00.21)
71.20
(+0033)
62.39 (+0.22)
Textile
83.08 (+0062)
83.12 (1.0633
86.03 (+0066)
86.26 (10059)
85.93 (+077
BreizhCrops
45.48 (+.2)
47.67 (+.)3
50.72 (=0.28)
54.66 (=0.45)
51.33 (0058
Table 2. Accuracy of the different classifiers on the different datasets consider.
in the classification process. So, during training, one only has to estimate yk for each class, which is the center of mass and can be estimated using the Riemannian mean (see proposition 7 of Said et al. 2018).
implementation the same as in Section 5. We call it the Homogeneous Wrapped Diseriminant Analysis (Ho-Wi) WDA).
As for the QDA, we propose another version of this classifier where each class has its own covariance matrix E3: ag WG(pk; Pk, Σx). We call this classifier the Heterogeneous Wrapped Discriminant Analysis (D-). In that case, an MLE is optimized on each class individually, as in Section
Tangent Space LDA or QDA For the Tangent Space LDA, we will leverage the Wrapped Gaussians introduced in Section 4. Suppose that each class is modeled by a wrapped Gaussian centered at G, the Riemannian mean of the training set, and with shared covariance matrix Σ. for the Tangent Space LDA i.e. Org WG(G; M4c, Σ) or with one covariance matrix Σg per class for the Tangent Space QDA i.c. a₃ = WG (6; pk: Σg). Then, the Tangent Space LDA (or Tangent Space QDA) is a maximum likelihood classifier based on those distributions.
Wk: E (1, K}, (pk. Pk: Ex) E &argmax II fp
P.K.E
6.4. Experiments
In this section, we want to compare the Ho- WDA and He- WDA to the other classifiers (MDM, Tangent Space LDA denoted TS-LDA and Tangent Space QDA denoted TS- QDA) used on the manifold of SPD matrices Pa and detailed in Section 6.1. For this, we lead some experiments on 8 different real datasets coming from several applications. We give summary of the datasets used in Table 1 and more details on each one of them in Appendix K. For this experiment, we restricted ourselves to the case where the covariance matrices are diagonal. Therefore, one has fewer coefficients to estimate and the MLEE needs less points to converge. We give the accuracy of the classifiers we study the different datasets at Table 2.
6.3. Wrapped Discriminant Analysis
Having placed the various classifiers that are used on the manifold of SPD matrices in a probabilistic framework, we propose a new maximum likelihood classifier based on the wrapped Gaussians introduced Section 4. First, let model each class by wrapped Gaussian with shared covariance matrix Σ among the classes: og WG(pk: Mk: Σ). To learn the parameters of each class, we optimize an MLE on the whole model to find the parameters:
K N
We can see two different behaviors. First, on the datasets with a lot of matrices (Textile, Salinas, Indiana Pines, Pavia Univ., BreizhCrops), the Ho- WDA and He- WDA perform the best. The number of parameters to estimate is high, so the more samples one has, the better the estimation will be
(pi...,,Re (1, -.. UK, Σ) E argmax
II k=1==
l xx)))
P.
where (I:))= are the training points of class The
8

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
as we illustrated in the synthetic experiments Section 5. For BreizhCrops, even if the matrices are of size 13 13, we have a lot of points (177, 658) so the MLE is able to correctly estimate all the parameters. Secondly, on the BCI datasets (BNC12014004 and Zhou2016), we have signif- icantly less points (less than 1, 1,00) so the estimation of the parameters of the underlying wrapped Gaussians is less precise. In this case, one can see that the MDM and the TS-LDA perform the best and the Ho- WD- and He-WDA perform less well. However, the AirQuality dataset, the Ho-WDA and He-WDA perform the best. This is interest- ing as the number of matrices available in the dataset is small (102). An explanation could be that the underlying distribution of the data is not very complex, so a few points are enough to correctly estimate them. In BCI datasets, the distribution is more complex, so one needs more points to correctly estimate the distribution. Finally, we do not ob- serve clear dominance of the He-WDA over the Ho-W WDA. This is similar to the difference between the LDA and QDA where, often, an LDA can correctly classify data. In the- ory, the Ho- WDA should be a special case of the He-WDA where the covariance matrices for each class are the same. However, in practice, as the number of points per class can be small and as the number of parameters to estimate is the estimation of the covariance
classical machine learning models that rely on Gaussian distributions the manifold of SPD matrices. We aware that computing exponential and logarithmic maps remains bottleneck in SPD matrix geometry. However, trade- off may exist between computational cost and performance gains. Theoretically, Euclidean Gaussian- methods ex- tend to Pa via our wrapped Gaussian, though practical chal- lenges will arise in applications, requiring careful choices. For example, one could develop Gaussian Mixture Models, Hidden Markov Models or Variational Autoencoders on Pa using wrapped Gaussians. These models could then be ap- plied to various tasks such as clustering, sequence modeling, or anomaly detection on manifold-valuee data. It could also be possible to explore the use of wrapped Gaussians in the context of Bayesian inference, which could open new avenues for probabilistic modeling and uncertainty quantifi- cation in manifold mandfold-based data analysis.
A Acknowledgments
This work was funded by the French National Research Agency for project PROTEUS (grant ANR-22-CE33-02 01). Part of this work was carried out while Florian Yger was a member of PSI.-Dauphine University, and he acknowl- edges the support of the ANR as part of the "Investisse- ments d'avenir" program, reference ANR-11--3IA-00 (PRAIRIE 31A Institute). Sylvain Chevallier is supported by DATAIA (ANR-IT-CONV-0003. The authors would also like to thank Bastien Cavarretta that helped inspire Theorem 4.7 and Lucas Gnecco Heredia for his help in proofreading the paper.
higher for the matrices per class can be noisy and the performance of the He-WDA can be worse than that of the Ho- HD- For ex- ample, for the dataset Salinas, some classes have only a few hundred samples, which is not enough to estimate the covariance matrix accurately. The Ho- ID- on the other hand, estimates single covariance matrix for all the classes and is less sensitive to the number of samples per class.
Impact statement
7. Conclusion
This work introduces a flexible and theoretically grounded extension of Gaussian distributions to the Riemannian man- ifold of Symmetric Positive Definite (SPD) matrices, which appear in numerous application domains such as neuroimag- ing (e.g.. EEG/BCI), remote sensing, atmospheric modeling, and computer vision. By leveraging a non-isotropic wrapped Gaussian model, our approach respects the intrinsic geome- try of SPD matrices, offering a more principled statistical modeling framework for manifold valued data.
In this work, we present a generalization of non-isotropic multivariate Gaussians on the manifold of SPD matrices: Wrapped Gaussians, and we give some theoretical prop- erties. We solved the non-identifiability of our model by defining an equivalence relation between the set of parame- ters that define the same wrapped Gaussian. We also give all the tools needed to use the distribution in practice, such as an easy-to-use sampling algorithm or an MLE that correctly estimates the parameters of wrapped Gaussian. Finally, we showed that the MDM, TS-LDA and TS-QDA classi- fiers can be seen as part of a probabilistic framework using wrapped Gaussians. We introduced two new classifiers based on the wrapped Gaussian: Ho- WDA and He-WDA. We showed that the Ho-WDA and He- HD-WDA perform well on real data when the number of samples is sufficient. In future work, we plan to investigate the use of wrapped Gaussians to perform data augmentation or transfer learning. More- over, as we have developed in the paper geometry-aware Gaussian distribution, it becomes possible to extend all the
The proposed distribution, along with the associated maxi- mum likelihood estimator and probabilistic classifiers, can enhance the interpretability. robustness, and effectiveness of machine learning models in applications where geometric constraints are crucial. For example, in Brain- Brain-Computer Interfaces (BCI), this work may contribute to more accurate and stable classification of neural signals, potentially bene- fiting assistive technologies. In environmental sciences, the method can aid in more accurate statistical modeling of air quality data represented as covariance matrices.
9

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
However, as with any advancement in data modeling, espe- cially in sensitive domains such as neuroscience, there is possibility that improved interpretabilityy or classification performance could be misused. For instance, fine aineggrined neural decoding could be exploited for persuasive technolo- gies or behavioral profiling, raising ethical concerns around privacy and consent. We thus emphasize the importance of applying these techniques within responsible and ethically guided frameworks.
Baumgardner, M. F,, Biehl, L. L.,, and Landgrebe, D. A. 220 band aviris hyperspectral image data set: June 12, 1992 indian pine test site 3, Sep 2015. URL https:
/ //purr. wwwppueedd/ppuuuuueepepppt wad.//ubblccaiiosss /1947/1.
Bergmann, P., Batzner, K., Fauser, M., Sattlegger, D., and Steger, C. The mvtec anomaly detection dataset: A comprehensive real-world dataset for unsupervised anomaly detection. Int. J. Comput. Vision, 129(4): 1038-1059, April 2021. ISSN 0920-5691, doi: 10.1007/ s1 11263-022001400000 URL https: /tl:://doooo aooiiorrg/l
References
Aristimunha, B., Carrara, I,, Guetschel, P., Sedlar, S., Ro- drigues, P., Sosulski, J., Narayanan, D., Bjareholt, E., Quentin, B., Schirrmeister, R. T., Kalunga, E., Darmet, L.. Gregoire, C., Abdul Hussain, A., Gatti, R., Gon- charenko, V., Thielen, J., Moreau, T., Roy, Y., Jayaram, V., Barachant, A., and Chevallier, S. Mother of all BCI Benchmarks, 2023. URL https:/ / github. com/ NeuroTechX/moss
Bhatia, R. Positive Definite Matrices. Princeton Series in Applied Mathematics. Princeton University Press, Prince- ton, N.J, 2007. ISBN
Bishop, C. M. Pattern Recognition and Machine Learning (Information Science and Statistics). Springer, edition, 2007. ISBN 0387310738.
Bogachev, V. I. Measure Theory. Springer, Berlin, Heidel- berg. 2007. ISBN 978-3-540-34513-8 978-3-540-34514- 5. doi: 10. 10.1007/978-3-540-445111
Arsigny, V., Fillard, P., Pennec, X., and Ayache, N. Fast and Simple Calculus on Tensors in the Log-Euclidean Framework. In Duncan, Medical Image Computing and Computer-Asissed In- tervention MICCAI 2005, pp. 115-122, Berlin, Hei- delberg, 2005. Springer. ISBN 978-3-540-32094-4 doi: 10. 15661651166.5555
and Gerig, G. (eds.),
Bouchard, F., Mian, A., Tiomoko, M., Ginolhac, G., and Pas- cal, F. Random matrix theory improved fr) 'echet mean of symmetric positive definite matrices. arXiv preprint S.Kiv:2405.065555 2024.
Boumal, N. An Introduction to Optimization on Smooth Manifolds. Cambridge University Press, edition, March 2023. ISBN 1------000--6666--1--011--1-- 978-1-00-916617-1 978- 1-00-916615-7. doi: 10.101779781009166611
Arsigny, V., Fillard, P., Pennec, X., and Ayache, N. Ge- ometric Means in Novel Vector Space Structure on Sysmetric Positive-Definite Matrices. SIAM J. Matrix Analysis Applications, 29:328-344 January 2006. doi: 10.1 13//0006666
Brooks, D., Schwander, O., Barbaresco, F,, Schneider, J.-Y., and Cord, M. Riemannian batch normalization for spd neural networks. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL https: .ww.:eddddiiiinnnnneeeeeeeeeeeee cc ae/paper__Bae co//apper_filles/ppprr///////// pdf.
Barachant, A., Bonnet, S., Congedo, M., and Jutten, C. Riemannian geometry applied BCI classification. In LVA/ICA 2010 9th International Conference on La- tent Variable Analysis and Signal Separation, volume 6365, pp. 629. Springer, September 2010. doi: 10.1007/ 978-3-648 11995-977777
Barachant, A., Bonnet, S., Congedo, M., and Jutten, C. Mul- ticlass Brain CramutCrmputer Interface Classification by Rie- mannian Geometry. IEEE Transactions Biomedical Engineering, 59(4):2009928 April 2012. ISSN 0018- 9294, 1558-2531. doi: 10.1109/BM..22111122
Casella, G. and Berger, R. Statistical Inference. Duxbury
Resource Center, June 2001. ISBN 0534243126.
Chen, Y., Wiesel, A., Eldar, Y. C., and III A. O. H. Shrinkage algorithms for mmse covariance estima- tion. IEEE Trans. Signal Process, 58(10)5016-55029 2010. URL ht / /twpp.lpppppco/cp www-trier........ wem.iiiirrrrdddabb.. fournals/tsp/tee http://pppeppppp.tttnceeeecaeeee htw.1
Barachant, A., Barthelemy, Q.. King, J.-R., Gramfort, A., Chevallier, S., Rodrigues, P. L. C., Olivetti, E., Gon- charenko, V., vom Berg, G. W., Reguig, G., Lebeurrier, A., Bjäreholt, E., Yamamoto, M. S., Clisson, P., Corsi, M.-C., Carrara, I., Mellot, A., Lopes, B. J., Gaisford, B., Mian, A., Andreev, A., Cattan, G., and Lebeurrier, A. pyriemann, October 2024. URL https:t//// /wt.dooooorooooo
Chen, Z., Song, Y., Liu, G., Kompella, R. R., Wu, X., and Sebe, N. Riemannian multinomial logistics regression for SPD neural networks. In Conference on Computer Vision and Pattern Recognition 2024, 2024.
10. 5281/211/// wwwp//eeooo.ooo 593.m 593816.
10

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
Chevallier, E. and Guigui, N. Wrapped statistical models manifolds: Motivations, the case SE(n), and generaliza- tion to symmetric spaces. In Joint Structures and Com- mon Foundations of Statistical Physics, Information Ge- ometry and Inference for Learning, Les Houches, France, July 2020.
Dudoit, S., Fridlyand, J., and Speed, T. P. Comparison of discrimination methods for the classification of tu- mors using gene expression data. Journal of the Amer- ican Statistical Association, 97(457):77 87, 2002. doi: 00.11988016224452275347777922999 URL tttps::///..... /doi. org/10.
Chevallier, E., D., Lu, Y., and Dunson, D. Exponential- Wrapped Distributions on Symmetric Spaces. SIAM Journal on Mathematics of Data Science, 4(4):1347- 1368, December 2022. ISSN 2577-0187. doi: 10.1137/ 21M1461551.
Galaz F., Papamichalis, M., Turnbull, K., Lu- nagome7, S., and Airoldi, E. Wrapped distributions on homogeneous riemannian manifolds. arXiv preprint arXiv:2204.09790, 2022.
Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning. Springer Series in Statistics. Springer, New York, NY, 2009. ISBN 978-0-38--888877- 978-0-387-84858-7. doi: 10.100079788---38--8488555
Chevallier, S., Kalunga, E. K., Barthélemy, Q. and Mona- celli, E. Review of Riemannian Distances and Diver- gences, Applied to SSVEP-based BCI. Newroinfornat- ics, 1911);93-006, January 2021. ISSN 1559-0089, doi: 10.1007/s1 00007//12221-00009947777777
Hauberg, S. Directional Statistics with the Spherical Normal Distribution. In 2018 21 International Conference Information Fusion (FUSION), pp. 704-711, 114471,,July 2018. doi: 10.23919/CCIF 8.49999999
Cho, S., Lee, J., Park, J., and Kim, D. A rotated hy- perbolic wrapped normal distribution for hierarchical representation learning. In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K. (eds.), Advances in Neural Information Processing Systems, 2022. URL https: //openreview net /otfffoore forumm www.rummidddnniiiie idxxHnbaaqzxx..
Hua, J., Zhang, Y., de Foy, B., Mei, X., Shang. J., and Feng. C. Competing PM2.5 and NO2 holiday effects in the Beijing area vary locally due to differences in residential coal burning and traffic patterns. Science of The Total Environment, 750: 141575 January 2021. ISSN 0048- 9697. doi: 10.101 2wc.itotenw.220.1444
Collas, A., Bouchard, F., Breloy, A., Ginolhac, G., Ren, C.,, and Ovarlez, J.-P. Probabilistic pea from heteroscedastic signals: Geometric framework and application to clus- tering. IEEE Transactions on Signal Processing, 69: 6546-6560, 2021. doi: 10.1 10.1109//SP.2221..111 1339997
Huang. Z. and Van Gool, L. A riemannian network for spd matrix learning. In Proceedings of the AAAI conference on artificial intelligence, volume 31, 2017.
Collett, D. and Lewis, T. Discriminating Between the Von Mises and Wrapped Normal Distributions. Australian Journal of Statistics, 23(1):73-79, 1981. ISSN 1467- 842X. doi: 10.111 1467-8427.78899 1981.,900763...
Jo, and Hwang, S. J. Generative modeling on manifolds through mixture of riemannian diffusion processes. In International Conference on Machine Learning, 2024.
Criscitiello, C. and Boumal, N. An Accelerated First- Order Method for Non-convex Optimization on Mani- folds. Foundations of Computational Mathematics, 23 (4):1433-1509, August 2023, ISSN 1615-3383, doi: 1110077/11020--022-09577777
Johnson, M. Multivariate Statistical Simulation: Guide Selecting and Generating Continuous Multivariate Distributions. Wiley Series in Probability and Statistics. Wiley, 1987. ISBN
Daletskii, L. and Krein, G. Integration and differentia- tion of functions of hermitian operators and applications to the theory of perturbations. AMS Translations (2), 47 10-10):001100, 1965.
Jona-Lasinio, G., Gelfand, A., and Jona-Lasinio, M. Spatial analysis of wave direction data using wrapped Gaussian processes. The Annals of Applied Statistics, 6(4), Decem- ber 2012. ISSN 1932-6157. doi: 10.1214/1-AAAAS766
Delmas, J.-P., El Korso, M. N., Pascal, F., and Fortunati, S. Elliptically symmetric distributions signal processing and machine learning. Springer Nature, 2024.
Kagan, A., Linnik, I., and Rao, C. Characterization Prob- lems in Mathematical Statistics. Probability and Statistics Series. Wiley, 1973. ISBN 978-0-471-44522-22
Ding, J. and Regev, A. Deep generative model embed- ding of single-cell RNA-Seg profiles on hyperspheres and hyperbolic spaces. Nature Communications, 12 (1):2554, May 2021. ISSN 2041-1723, doi: 10.1038/
Ledoit, O. and Wolf, M. A well-conditioned estimator for Large-dimensiona covariance matrices. Journal of Multivariate Analysis, 88(2):365 February 2004. ISSN 0047259X. doi: 10. 121101690044-255XX333333
-41467-001-22851-4

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
Lee, J. M. Introduction to Riemannian Manifolds, volume 176 of Graduate Texts in Mathematics. Springer Interna- tional Publishing, Cham, 2018. ISBN 978-3-319-91755-2 978-3-319-9175599 doi:
Nguyen, X. S. GeomNet: A Neural Network Based on Rie-
mannian Geometries of SPD Matrix Space and Cholesky Space for 3D Skeleton-Based Interaction Recognition. In 2021 IEEE/CVF International Conference on Com- puler Vision (ICCV), Montreal, Canada, October 2021. doi: 10. 10.1099/CCV449222220..0112 / //hal hhhallscienneehaa-000022222
Leeb, R., Lee, F., Keinrath, C., Scherer, R., Bischof, H., and Pfurtscheller, G. Brain-computer communication: Motivation, aim, and impact of exploring a virtual apart- ment. IEEE Transactions on Neural Systems and Re- habilitation Engineering, 15(4):473-482 2007. doi: 10.1109/TNSREE2207..66666
URL https:
Pedregosa, F, Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cour- napeau, D., Brucher, M., Perrot, M., and Duchesnay, E. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
Lehmann, E. L. and Casella, G. Theory of Point Estimation. Springer- Verlag, New York, NY, USA, second edition, 1998.
Liu, J., Liu, C., Shi, J. Q. and Nye, T. Wrapped gaussian process functional regression model for batch data on riemannian manifolds. arXiv preprint 13X1,:2409.03111,, 2024.
Pennec, X. Intrinsic Statistics on Riemannian Manifolds:
Basic Tools for Geometric Measurements. Journal of Mathematical Imaging and Vision, 25(1): 127-14,,Jull July 2006. ISSN 0924-9907, 1573-7683. s 11085--0066228--
doi: 10.1007/
Lopez-Custodio, P. A cheat sheet for probability distributions of orientational data. arXiv preprint arXiv:2412.08934, 2024.
Pennec, X.
Curvature effects on the empirical mean in rie- mannian and affine manifolds: non-asympiotic high con- centration expansion in the small-sample regime. arXiv preprint arXiv: 2019.
Lotte, F., Bougrain, I.., Cichocki, A., Clere, M., Congedo, M., Rakotomamonjy, A., and Yger, F. A review of clas- sification algorithms for EEG-based crain-computer in- terfaces: A 10 year update. Journal of Neural Engineer- ing. 15(3):031005, April 2018. ISSN 1741-2552. doi: 10.1088//741122552abb222
Pennec, X. Manifold-valued image processing with SPD trices. In Pennec, X., Sommer, S., and Fletcher, T. (eds.), Riemannian Geometric Statistics in Medical Image Anal- ysis, pp. 75-134, Academic Press, January 2020. ISBN 978-0-12-814725.1 doi: 00.0016/97780011--1147777 00010-8.
Mallasto, A. and Feragen, A. Wrapped Gaussian Pro- cess Regression on Riemannian Manifolds. In 2018 IEEE/CVF Conference on Computer Vision and Pat- tern Recognition, pp. 5580-5588, Salt Lake City, UT, June 2018. IEEE. ISBN 978-1-5386-6620-9 doi:
Peyré, G. and Cuturi, M. Computational Optimal Transport, March 2020.
10.1109/CVPR.2116 005055
Pfurtscheller, G. and Neuper, C. Motor imagery and direct brain-computer communication. Proceedings of the IEEE, 89(7)11233-1134 2001. doi: (0.1109//.99829)
Mardia, K. V. and Jupp. P. Directional Statistics. Wiley Series in Probability and Statisties. Wiley, Chichester, new ed. edition, 2000. ISBN 978-0-471-95533--..
Rußwurm, M., Pelletier, C., Zollner, M., Lefèvre, S., and Körner, M. Breizhcrops: A time series dataset for crop type mapping. International Archives of the Photogram- metry, Remote Sensing and Spatial Information Sciences ISPRS (2020), 2020.
Mathieu, E., Le Lan, C., Maddison, C. J., Tomioka, R., and Teh, Y. W. Continuous Hierarchical Representations with Poincaré Variational Auto-Encoders, In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Said, S., Hajri, H., Bombrun, L., and Vemuri, B. C. Gaus- sian Distributions on Riemannian Symmetric Spaces: Statistical Learning With Structured Covariance Matri- ces. IEEE Transactions on Information Theory, 64 (2):552-772, February 2018. ISSN 1557-9654, doi: 10. 0..110997T22017221333999
Moakher, M. A Differential Geometric Approach to the Ge- ometric Mean of Symmetric Desitite-Definite Matrices.
SIAM Journal Matrix Analysis and Applications, 26 (0):735-447, January 2005. ISSN 0895-4798, 1095-7162. doi: 10. 10.11377508854779933369777
Nagano, Y., Yamaguchi, S., Fujita, Y., and Koyama, M. A wrapped normal distribution on hyperbolic space for gradient-based learning. In International Conference on Machine Learning, pp. 4693 4702. PMLR, 2019.
Sanborn, S., Mathe. J., Papillon, M., Buracas, D., Lille- mark, H. J., Shewmake, C., Bertics, A., Pennec, X., and Miolane, N. Beyond euclid: An illustrated guide to mod- ern machine learning with geometric, topological, and
12

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
algebraic structures. arXiv preprini arXiv:2407.09468. 2024.
Wasserman, L. All of Statistics: Concise Course in Sta- tistical Inference. Springer Texts in Statistics. Springer, New York, NY, 2004. ISBN 978-1-4419-2322-6 978-0- 387-21736-9, doi: 1007/978 0-387-21736--
Shiga, K. Hadamard Manifolds. In Geometry of Geodesics and Related Topics, volume 3, pp. 239-282, Mathemati- cal Society of Japan, January 1984. doi: 10.2969/aspm/ 00310239.
Willjuice Iruthayarajan, M. and Baskar, S. Covariance ma- trix adaptation evolution strategy based design of central- ized PID controller. Expert Systems with Applications, 37(8):5775-578) August 2010. ISSN 0957 4174. doi: 10. 10.1016ffeeea.221000....33
Smith, A., Hua, J., de Foy, B., Schauer, J., and Zavala, V. Multi-Site, Multi-Pollutant Atmospheric Data Analysis Using Riemannian Geometry, December 2022.
Zhou, B., Wu, X., Lv, Z., Zhang, L., and Guo, X. A fully automated trial selection method for optimiza- tion of motor imagery based brain-computer interface. PLOS ONE, 1 09 2016. doi: 10. 10. pone. 0162671 URL ht https: :tta:///oooo htt//ddi..rtt 1371/ journal. 0162657.
Smith, A., Hua, J., de Foy, B., Schauer, J. J., and Zavala, V. M. Multi-site, multi pollltaollut atmospheric data anal- ysis using riemannian geometry. Science of The Total Environment, 892:164064 2023.
Sra, S. Positive definite matrices and the S-divergence. Proceedings of the American Mathematical Society, 144 (7):2787-2797, October 2015. ISSN 0002-9939, 1088- 6826. doi: 10.1090/proc/12953.
Sra, S. and Hosseini, R. Conic geometric optimisation the manifold of positive definite matrices. SIAM Journal on Optimization, 25(1):713-733 January 2015. ISSN 1052-6234, 1095-7189. doi: 10.1137/1449988168
Terras, A. Harmonic Analysis on Symmetric Spaces and Applications II. Springer, New York, NY, 1988. ISBN doi: 10.1007/
978-0-387-96663-2 978-1-4612-3820- 978- 18811-1612-82-11.
Thanwerdas, Riemannian and Stratified Geometries on Covariance and Correlation Matrices. PhD thesis, Uni- versité Côte d' Azur, 2022.
Townsend, J., Koep, N., and Weichwald, Pymanopt: A python toolbox for optimization manifolds using automatic differentiation. Journal of Machine Learning Research, 17(137):1-5, 2016. URL http:// jmlr. org/papers /vl.s/ss/////// 477.html.
Troshin, S. and Niculae, V. Wrapped $ SbbetaS-gaussinae with compact support for exact probabilistic model- ing on manifolds. Transactions on Machine Leam- ing Research, 2023. ISSN 2835-8856. URL. https: //opereeviiw..c. wtwp::/wwwww.....m.m...... Kreekrur
Turaga, P., Veeraraghavan, A., Srivastava, A., and Chel- lappa, R. Statistical Computations on Grassmann and Stiefel Manifolds for Image and Video-Bil Based Recogni- tion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(1 33(11))2273-22) 0162-8828, 2160-9292,
November 2011. ISSN doi: 10. 10.1109/TPMM..2011..0
Tuzel, O., Porikli, F., and Meer, P. Pedestrian detection via classification on riemannian manifolds. IEEE Transac- tions on Pattern Analysis and Machine Intelligence, 1177 2008. doi: 10.1 1099/PAMM..200.755
13

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
A. The vectorization used
In Section 3.2, we defined an isomorphism between the tangent space TaPa at point p E Pa and R This isomorphism called the vectorization and is denoted Vectp. This isomorphism is not the only one that exists between the two spaces, so in this section we will motivate our choice of the vectorization and give some of its properties.
Let E Pa- We recall that the tangent space Pu at point p is Euclidean space once equipped with the inner product defined at Equation (1). Therefore, one can unveil an orthonormal basis of Pa. One can see that the tangent space at the identity I., Pa is the classical Euclidean space Sa equipped with the Frobenius inner product. Therefore, one can easily build an orthonormal basis of Ti,7 Pa and then, transport it to the other tangent spaces.
Proposition A.1 (Orthonormal basis of the tangent spaces), Let Eij be the d x d matrix with a I at position (i.j) and zeros everywhere else. Then
An orthonormal basis of PTi.Pa (())1)) is (Etast)is) defined as follows:
1 (et + tgx)
for i < j,
Fiadj
v2
Csi
for = j.
An orthonormal basis of (TpPa, )p) is (Epss)es) where Ep,is 1p//E///saa Elaiigp
Proof. One has that Pa = Sa and that (t,) is the Frobenius inner product, so one can use the classical basis of Sa to build an orthonormal basis of Pa, Then, by transporting the basis of Pa to Pa using the isometry x + p xy 1/2 one has a basis of TuP It is still orthonormal as x + p'/22pp' xp3/2/ 12 is an isometry.
Let us give another intuition on the choice of this basis for (TyPa. (1 Proposition A.2. The basis (Ep.g)is)) of (TpPd, (p) given at Proposition A. is the parallel transport of the basis (Erass)iss of Pa. (1) from Ti, to T,Pa.
Proof. According to Equation 22 of Sra & Hosseini (2015), in the case of Pa. the parallel transport I-p from PipPa to TpPa is:
Yu E Pap (u))p(u)= 11/2up
The result follows from the definition of (Epsi)iss.
Now that we have orthonormal basis of the tangent space TyPaa, we give the link between this basis and the vectorization Vectp:
Proposition A.3. Let (Pe.ii)sss be the orthonormal basis of the tangent space TuPa a described Proposition A. I. Let PsP Then,
(w) ((u, Ep. 1)pp (u, Ep, p2) (u, Ep.22) ... (w, (hna-1d)m (M, p))
Proof. We start with the case where p Ia. Let a [lusy]l E Ti, Sg. We simply need to show that, for i S j, one
has
- { Uni
if = j, if <
(u, Elasshis -
- (22mg
One has, when i = j:
(u, Elatt)aa = (u, eviia = tr(uess) = uit
And when i < j:
(u, Elatshas = (u. v2 (ess + egillia = v2
(tr(uess) + tr(uesi)) = V2u...
Therefore, one has the results for Vects
14

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
Now, in the general case of p Pa. one has that, for i S
(u, Enss) (u.p' F1atsp/))) (p upp 1/2 Eaa
By using the definition of (uct,(u) Vect,, (p up-1/12p11/3) and the result for p Las one has the result.
1/2)
The previous proposition helps motivate the choice of the vectorization. Indeed, this vectorization is simply the coordinates of the tangent vector in the orthonormal basis (Fg. of the tangent space Pa. Now that we are more convinced on this choice of vectorization, we give some of its important properties. Proposition A.4. Let E TyPa. then
I| (eccp(u))l Vectp(u)" Vectp(u) = 1/ull v u),
Therefore, Vect, is not only an isomorphism, it an isometry between (I,Pa: (:s 1)) and (IR (Ra4, || 112).
Let u E PyPd. then || Vectp(Log, u) 12 = (eog,(Log, u)T Veclp(Log, u) = 6),,u)²
Proof. Let us prove the two points of the proposition.
Let u E TaPa. One has, using Proposition A.3,
|| Vectp (u)1 = (ectp(u) Vectp(u) = Σ
(u, Ep.isl
(5
As (Epay)isz is an orthonormal basis of the Euclidean space (InPa, (p). one has that
uu Ep.s) 12 = |10/|p-
which proves the first point.
Let u E PaP One has, using the previous point, the definition of || Ipp- the expression of the Riemannian logarithm given at Equation (3) and the expression of the AIRM distance given at Equation (2):
|| (eot u)ll2 = Log, u I|p Log, u p 1/2 1-11//11 = || log(p 1/2 log(p-1/2upp w//2up-1/3)11 = 6(p, u)
Finally, let us give a consequence of the previous property on the Jacobian of the vectorization: as Vect, is an isometry, there is no volume change via the vectorization.
B. The push-forward
In this section, we give the definition and an important result on the push-forward measure. One can find more information on the push-forward measure in Section 3.6 of Bogachev (2007).
Definition B.1 (Pushforward measure). Given two measurable spaces (X', f2x) and (J,, Qy), a measurable map f: Sx + By and measure p:: 2x + [0, +oo], the pushforward of y is defined be the measure f #u Gy [0, +oo] given by
vB E Sy, (f#r)(B) (B)-
where f [B] is the preimage of B by f-
We now give the most important result on pushforward measures:
the change of variables
Theorem B.2 (Change of variables). Let μ be a non-negative measure. An mey-maasurable function g on y integrable with respect to the pushforward measure # u if and only if the function g o f integrable with respect to the measure p. In this case, has:
/ g(0)d((HH)(0)
J (g f)())d((x).
15

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
C. The Jacobian determinant of the exponential map on the manifold of SPD matrices
In this section, we give a proof of the result stated at Proposition 4.3, More specifically, we want to show that the Jacobian determinant of the exponential map Exp,, at identity is:
sinh A1(w)->,(u)
Jra 2dh
1/22
dc(u) (T(
i<j
where the (A)( are the eigenvalues of u.
Let us start by recalling that, when the tangent plan of interest is at the identity Ia, the Riemannian exponential map Exp1, is simply the matrix exponential exp:
Expr, u E Pa + exp(u) E Pa.
One can see this result using the expression of the Riemannian exponential given in Equation (3).
To prove the relation Equation (7), we will start by the case where u is a diagonal matrix. We will then extend the result to the general case.
Case 1: E Pu diagonal Let us consider u Tı,Po Pa diagonal, u diag(A1..... Ag). In the following, we will denote by V the differential of the Riemannian exponential in V = d Expr, (u) = d exp(u). We therefore want to compute the determinant of Vr det W One has that v Pa Ta, To compute the determinant, one need to choose adequate bases in both tangent spaces Ti,Pa and Texp(a)PA, Pa. By "adequate", we mean that the transformation between the two bases does not imply any volume change. For this, we consider for Pu the basis (Elaaj)ig and for Texp(u) Pa the basis (Eexp(n),&j) as defined at Proposition A.1. According to Proposition A.2, the transformation from the first to the second basis is the parallel transport, which does not imply any volume changes, as the parallel transport is an isometry (see Prop 10.36 of Boumal 2023).
Pa where we have identified T, with
Now that have two basis, want compute the matrix of V in those two bases. For this, we need to compute T Ery, and express it in the basis ( (Eexp(u).s) )isj. (y)sj-- As u is diagonal, we can use the Duletskii-Kred formula (see Daletskii & Krein (1965) or Equation 2.40 of Bhatia 2007) that states the following in our case: for h E TpPd.
(h)
[[ep 11 (w)
(8)
where
for i = j,
expli(())
for + j-
Using the previous formula, one can compute Ti (EL for S j. Now, one needs to compute the coefficients of V (EL,,5)) in the basis (Eexp(u),kt)esst As the basis (Eexp(w). kss orthonormal, one simply needs to compute the dot product between V (E1a,))) ) and one element of the basis to get the corresponding coefficient. For S one has:
(VD (E1ass)))
(EIg,ss) exp(u/2))r,,..k ext.s/xxxpp2222 (2))()) exp(-u/22 (cxp(-u/2) (E1,, (()ss)) exp( u/2), E1
using the definition of Exxpts) using the definition of the AIRM metric
tr (V)
Therefore, it is the coefficient (ks 1) of the matrix need to compute this matrix. As u is diagonal, u = diag (A1... Aa), one has that exp( (op(-uu2) diag(e >1/2 ..... e Az/2). Therefore, and using Equation (8), the coefficient (k, 1) of the matrix exp(- (op(-u/2)) u/2) ( Fig exp( 2/2) is:
(Erg,is) exp(-u (up to factor v2 when k / l). We now
e >1/2 Xr/2gdie -xx/12
for i = = k =
sinh ( Ai/2 c -ots c- -2/222 v2 A; A
for i # j and (k:, (i, if (k. # (i.j)-
v A
Therefore, the matrix of v in the bases (Fsss and (Fexp(s).a)est is diagonal with diagonal coefficients and 2 sinh Ar-A²) Ao A3
Thus, the determinant of the product of these coefficients, which gives the result for the diagonal case.
16

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
Case 2: E Pa general Let us now consider the general case where u E Ti,Pa is not diagonal. One has that, as u is symmetric, one can diagonalize it: u gdg" where g is an orthogonal matrix and d is diagonal. As exp(u) g exp(d)g), " one has that d exp(u) h = 9 (dexp(d) (g hg)) g Therefore, of d (w), and as g are orthogonal, one has that det d exp(u) h det (d edp(d)) (g hg)). One thus need compute exp((d) (g) hg), and using Kreiin formula, as in the diagonal case, one has:
are only interested in the determinant
exxp(d) (g) hg) (g) [exp D (w)
where (u) is defined as above and g² gh modify the basis used in Pa in Equation (9), the coefficients of g hy appear (rather than directly the coefficients of h as in the previous case). Therefore, we choose as basis for Ti. the basis ( where E(a) g T. One can easily check that this basis is orthonormal and does not imply any volume changes. One can now use the same proof as for
[[hay]]- |- In order to do the same proof as for the diagonal case, one needs to
the diagonal case compute the determinant of d exp(u), which gives the result for the general case.
D. The building block of the wrapped Gaussians
In this section, we give a proof of Proposition 4.5. For this, we will actually show a more precise proposition: Proposition D.1. Ler (p, 14 Σ) E O and X ~ H, Then,
1. p- 1p-1//2pp11/ 1/2 ~ WG(
2. Exp,(Togp X "ect)"(1) ~ WG(p; Ou(+1)/2) E). (22 ( Vect, X])) NG
3. Exp,,(Vect,
Proof. In this proof, we will only show the first two points of the above definition, the third one being similar to them.
1. Let Y p p-1/2Xp- 1/2, We want to show that Y WG pt, E). For this, let pp: Pa + R be a continuously bounded function. One has
Bip(Y)) f Jr. 1pp-112200 1/2)
= .fa p0p plp-112200 112
1
exp 1 Vectp(Log, (x)) - - 1 /Vectp(Log, 1i) dvol(ar). |JI, (log(p 1/22p 1/2)]]
)2 d Ee
Let us now define vpp: I + p-1221/2 1/2, vp is -'ddiffeomorphism between Pa and Pa. Moreover, as the volume clement dvol is invariant by congruence of GL,d, R), the transformation top does not imply any volume change. Therefore, by change of variables y = (x). one has:
Bip(Y)) = S p(v)
1
exp (- 1 (Vect, (Log, 22202 11) T Σ E (Vect, (p2/2 1)) dvol(y) (11glo((p11/2 2/2p2222pp2/////// 22pp1p11ppppp1111
Pa
V(2 d
= 1 pp p(y)
exp (- 1 (Vect,, (p 1/2 1og(v))112 11) Σ (Vect (1) (22 log log(u)///2) 1n) dvol(y). 11 (log(u))
(2T) det
Finally, using that Vect, (2 10) = (p 1/2 1/ log(u)p 1/2) = (logl()) and Logi, log, we have,
E[((Y)) = / ((v) (2 det
exp (- 1 (Vect1, Veect,(Logg,,eeo (v)) - 1) I (Logi, - -R)) dvol(y). |Jr, (log(y))
This shows us that Y WG(Ta: 11,
17

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
2. Let now Y (Logp Vect, (w)). We want show that Y (((++)22))) E). Let g:: Pa + R be bounded continuous function. One has
E(p(Y\)
((Expp Vect, (w)))
pp
S Pa
p(Exp, (Logp
Vectp -1 (1)))
exp (- 1 (Vectp (Log, (x)) - T (Veett(VV, - -1)) dvol(r).
x
v (2 det E
[Jg(Log((z))
Let us now define Up: + (Logp Vect, (k)). (u)). vp Cl difffoomorpphim between Pa and Pa and its inverse is : + Expp (Logp I + (m)). By change of variables (,)(2), one has:
/ ((v) Pa
exp (- (Log, T -
dvol(y) det dupp (w--(1)))
E|b(Y)]
((2x)d det Σ
(p(L (1) + Vect, (n))
We need to compute the change of volume term det ( ( For this, let us start by saying that d)p(()) d Expp, x - Now using the fact that (LLogp()) = (dExp, Theorem 4.2), we have
o d Log, x therefore, det (2) = det d Exp, (Log, - Vect, (H4) det d Log, x.
and the definition of Jp(u) = det d ( (see
det d = (Log, x - Voct, (1)) Jp(Logp
and thus, plugging (2) into the equation:
det dtt (Logp v)
Jp(Log y + Vect, (1))
Therefore,
E [[((()] = / p(v)
1
exp (-1 Voct, (Log, T Veet, Voetp(Log, " (v))) dvol(y).
(2 det
1J (Log (v))
This shows us that ~ WG(p; 02(41)/2; E).
3. For the third point, can prove it similarly as the two previous one, having in mind that the vectorization Vect, is an isometry (see Proposition A.4) therefore, neither Vectp nor Vect, -1 implies any volume changes.
Therefore, Proposition 4.5 is a direct corollary of the previous result:
Corollary D.2. Let X ~ WG(La; 024411/2 la(+))))) and let (p, n, E) O. Let us define
W: N Pa 1/2 (Vect,, (5)/2 (Veet, o Log, 11)) 12/2.
Then, (X) W.)
E. The wrapped Central Limit Theorem
In this section, we give a proof of the wrapped Central Limit Theorem stated at Theorem 4.7. For this, let (Xt)en be sequence of independent and identically distributed random variables on the Riemannian manifold Pa. We suppose that the sequence (,ect,, (Logı, (X,)))eeR~ of random variables on IB 40d+1/// satisfies the classical Central Limit Theorem. We want show that the sequence (X;) satisfies the wrapped Central Limit Theorem.
18

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
Remark E.1. Let start by saying that, the map + V, Log1, x is a diffeomorphism between Pa and IR Raddt//, the sequence (Vectz, (Togı, (X,))) is also independent and identically distributed.
As the sequence (Vecty, (Togr, (X, xe)) satisfies the classical Central Limit Theorem in 1/ has that,
Veotz - p)) n--00
N/ Σ).
By defining m Expr, (Vect (m)) Pa. and using the linearity of Vectp:
1 vn
Va
(Logr, (X,) Logr,(m))
n-+00
Thus, by applying the continuous map Exp,, Vect to the previous equation, by considering the fact that the convergence in distribution is stable by continuous maps (theorem 5.5 of Wasserman 2004) and the definition of the wrapped Gaussian (Definition 4.1), one has that
a
ExPra
Σ
(LogI, (X.)
Logis (m))
WG([s;0,,))
We can now simplify the left-hand side of the previous equation:
ExPιa
(Log,, (X,) - Logi, (m))
exp
(log X; - llog m)
using the expression of Expr, and Logi, of Equation (3)
exp
(log X; log logmmmm!! m-
using log log and exp(aar) = exp(xx) for R
i=1
log(X, 2 mm))
exp
using delinition logarithmic product
1
vn X 3 m
Therefore, one has the final result of the wrapped Central Limit Theorem:
Ö X, m00mm))
W,(
The generalized wrapped CLT identity matrix la. However, one can generalize this theorem to any point p Pa- For this, we need to introduce generalized logarithmic product @p between two points 41,92 E Pa:
The previous version of the wrapped Central Limit Theorem is centered around the
Q1 Op 92 Exp, (Log, a + Logp 92).
In the same way as before, one can show the following generalized wrapped CLT
V d O X, Op --11
WG(p; 0, E).
22000
19

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
F. An extension to wrapped Elliptically Contoured Distributions
As one can see from the expression of the density of wrapped Gaussian (G μa, 2) given at Theorem 4.2, is intrinsically linked to the density of the multivariate Gaussian N (p, 2). This suggests a possible extension to Elliptically Contoured Distributions (chapter 6 of Johnson (1987) or Delmas et al. 2024). We recall the definition of an Elliptically Contoured distribution:
Definition F.1 (Elliptically Contoured Distribution). random vector X E IR follows an Elliptically Contoured distribution if there exists y R Σ E Pa and function g such that X has density
fx(x) ls det ((x - u) " 4))
where & normalizing factor. We denote X EC E)
For example, the multivariate Gaussian (u, Σ) is an Elliptically Contoured distribution with g: + exp(-1/2). Another example is the multivariate t-distributions for which g: Ê (1+) 4/1/2 (see V.B of chapter of Delmas et al. 2024).
One can then extend what has been done previously on the wrapped Gaussian to define Wrapped Elliptically Contoured Distributions just like above:
Definition F.2 (Wrapped Elliptically Contoured). Let p E Pd, E Rddd+//, ,E Pac+n/2 and g be a function. Then, a random vector on Pa follows a Wrapped Elliptically Contoured denoted WEC μ4, E,g) if
X = Exp, (t)), ~~ EC(u, E,
One can then compute the density of WEC μl Σ, g) similarly as in Theorem 4.2. Moreover, all the work done on the equivalence relation for wrapped Gaussians stays valid for wrapped elliptically contoured distributions.
G. The proofs on the equivalence relation
In this section, want give proofs of the different results of Section 4.4. We recall the propositions and give their proofs. Let us start by Proposition 4.9. Proposition G.1. Let (p, pt Σ) θ and I E TR. One has that WG(p; pb E) and WG(e'p:: Lu, E) are equal where (p) .I, .. 1,0, ... .0) 1//////
Proof. In this following, we denote by y the function y: t ト+ e'p. Let us denote by f the density of WG(Y(f); tv, E) and by f the density of WG u, E). We want to show that f = Let x E Pa. by Theorem 4.2, one has:
- (Veet (Log + tv) T E-1 (Vectrco)Log - 1 + tvv)
exp
1
f(z)
(( )d det Σ
Jn(()) (L0g)(1((x))
One has, that
Log (x) 12 log(:(()=1/222)))11111))))
= =e'p' logle 1/2 1/2))////
using that 7(1) = e¹
e'p 12 log(e *Ta) 1/2 2 Log,(x) using that e- --1 Ta and p xpp 1/2 commute
= -te'pt + Log,
Furthermore, one has:
(Log (x))
te Vocty (p) Vo (Log)) using the linearity of Vect,(t)
-te' Vect (occt),(((t)1111)) (1(t) 1/2) + (x(t) 1/2 (z) 1/2 Log, (x)p 1/2) -t Vecti pp 1/2) + c Vecti, (p-1/2
-t (ect,(p) + Vect, (Log, (Log,
20

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
Therefore, the numerator of the density / can be rewritten as:
(-1 (Vec))(Log,(())) - --
'VVecl (Logpp())) - - 11) n))
exp
which is the same numerator as f-
Let us now focus on the denominator. One has:
(z)g = JI,(7(t) 1/2 Log Lo) (z))(()))
= J1, Aalee'p-1//0 1/2 1o0 p (x) 1/2)
= J1,(-tla + p--pp-p00 Log,
Log 10)p-11))
using the computation Log, (I))
We recall that the Jacobian determinant of the exponential map the identity is:
sinh A
J1 =
2 II (x 22(2)
i
where (Ar(s)): are the eigenvalues of u. Moreover, the eigenvalues of : == - attl + 11/2pp111111/
1//1/////g Log,(x))
1/2 11/2 are
+ X1 + Au(u)= = -at
1/2 (1-111// (z)p
Thus, for all < j. one has:
(a(u) A5(1)= = di (1-11//T.0 (x)) 112) A3 (p-1/2 Log, (z)p 1/2)
and therefore, this leads to:
(1)( = J1s (p-1/210m) (x)p 1112) = Jp (Logp(x))
So the denominator the density f is the same as the denominator of the density f and therefore, the two densitics are equal.
Remark G.2. The function 7: + etp is actually the geodesic with initial point p and initial velocity p. Indeed, the expression of the geodesic L, (t) with initial point q and initial velocity V E TT is (see Pennec 2020):
Vt E IR, (t) = -2 exp(tt 1VV
1
Therefore, the geodesic 7 with initial point p and initial velocity p (which is symmetric matrix, therefore an element of TuPa Sq) is:
1A E R, (1) = =epp=(1))
We now want to show Proposition 4.14 that we recall underneath:
Proposition G.3. Let 0 (p: μ, Σ) E θ be a tuple of parameters. Then, the minimal representative of the class [0] as defined at Definition 13 is Omnn (y/nin y min Mimiay where
min
ed 2im1 Hip,
d
min
uav.
1=1
Symin
Σ.
Proof. We want find the smallest mi in the sens of ll2 and the corresponding y such that (15 E) ~ (ymin, min 52min).. ). As all the l in the equivalence class of [0] are of the form μ tu for IR, to find the smallest umin one needs to minimize the following function:
: || tol3 11 20 v) + 11111112
21

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
One thus has:
<'(1) 2( 1) 20111112-
The minimum reached tmin
with ||0112 d and (H.1) Ei= Mi. Therefore, one has:
p min
ed
Ei
Hip
d
1
11 min
1
H&u
[1
N1: .. pd
Ni: pa+1," -Mal+11)/2
d
H. Why does estimating p using the Riemannian mean fails in the general case?
We said in Section that when μ² * O, using the Riemannian mean ( does not work to estimate the parameters (p* ,=" Σ²). Let us explain why. For this, we suppose in the following that μ² * 0. Let Pvv be the Riemannian mean: PN B (211 ZN). Then, we can use Proposition 5.1 compute the MLE of μ and Σ:
AN
VLoggy (Iz),
ΣN
(Loog (ax) (NN) (VLogan (2t) - pNN)
N
where we recall that VLogav is the vectorization at pv of Logø i.e. VLoggN Using the linearity of Vect,x we can write that
VecteN o Logan Let us focus on fx.
pN
Vecl
Σ
Loypw (1::)
N
According to proposition 3.4 of Moakher (2005), as P. the Riemannian mean of the points (I1,..., RN), we have the following:
Σ Logpy (xx) 0.
Therefore, AN = 0. It therefore not a good estimator of 1²' + O. That is why we do not use the Riemannian mean as an estimator of a general setting when we do not know a priori that p = 0.
More details on the MLE experiments
In this section, we give more details on the setup of the synthetic experiments lead in Section 5 to assess the quality of the estimation of parameters of a wrapped Gaussian using an MLE. Upon acceptation, we will release the code used to perform these experiments. To obtain the results plotted at Figure 3, we repeated 5 times the experiment with different true parameters 8* = (p", ,²= Σ*) randomly generated. Here are details on how we generated the true parameters:
For p", we use the function generate www.errttterandmmsspppepeeee www.dmmsppdmattri. from the library PyRiemann (Barachant et al., 2024). This function generates a random SPD matrix by generating a random matrix and then computing exp((X + 8 A) (X + s * A)) where X and s are parameters chosen by the user. We set = 0.14g and s =
For ju,, we generate random vector of size did+1)/2 with values in [0, 0.1..
For * we generate random SPD matrix using the same function as for p= with X 0.01 and s 0.02.
We chose relatively small values for and s because otherwise, when the dimension is large, the generated parameters are very far from identity leading to numerical instability.
22

---

Wrapped Gaussian on the manifold of Symmetrie Positive Definite Matrices
Distanne to the roal
Distamee to the roal
Distance to the roil X
Conritann modere
Conariana matri St
1P-
Chnertance notris X:
Fil
Daped
2M46
3011
Nunber poiats
Nasiaa
Namlen
Figure 4. Comparison of the estimation of the parameters of wrapped Gaussian when the covariance matrix Σ is diagonal or full. The dimension of the SPD matrices in this experiment is d = 10.
J. Estimating the parameters of a wrapped Gaussian when the covariance matrix Σ is diagonal
In this section, we detail the experiments on the estimation of the parameters of wrapped Gaussian when the covariance matrix E is diagonal. We used the same setup as in the previous experiments detailed in Appendix I except that we generated E' as diagonal matrix. The diagonal was uniformly sampled in [0, 1J. We repeated the experiment 5 times with different true parameters. The goal for this experiment was to assess the impact of the structure of the covariance matrix on the estimation of the parameters. We wanted to show that when the covariance matrix is diagonal, the estimation of the parameters requires fewer samples. We plotted the results in the case of dimension = 10 at Figure 4. One can see that, while the estimation of p and ju is not affected by the structure of Σ, the estimation of Σ is better when E is diagonal. With a lot less samples, one can achieve significantly better estimation of Σ when it is diagonal. This is coherent as the number of coefficients to estimate is reduced when Σ is diagonal. Therefore, this setting of Σ diagonal can be a good choice when the number of samples is limited.
K. More details on the real data experiments
Let us start by giving more details on the different datasets used in the experiments and the preprocessing we used.
BCT Datasets: We considered 2 datasets from Brain Computer Interfaces (BCT) for our experiments: BNC1201400/ (Leeb et al., 2007) and Zhou2016 (Zhou et al., 2016). They consist of several subjects and several sessions per subjects doing Motor Imagery task (Pfurtscheller & Neuper, 2001). We used the library MOABB (Aristimunha et al., 2023) to load and preprocess the data. For each EEG, We start by applying a standard band-pass filter with range (7: 35] Hz. Then, we used the Ledoit-Woo shrunk covariance matrix (Ledoit & Wolf, 2004) to in order to compute the covariance matrices and to avoid ill-conditioned matrices. The experiment we lead was cross-subject: each classifier was trained on all subject except one and tested on this last subject.
AirQuality: This dataset is from the Beijing Municipal Monitoring Center. It is dataset of air quality monitored from 34 different sites in Beijing, China (Hua et al., 2021). For each site, six atmospheric pollutants where recorded every hour: CO, NO2, O3, PM PM2.5 and SO2. We used the same preprocessing as in Smith et al. (2022) to get a point cloud of 102 covariance matrices of size 6 x 6. Each covariance matrix has a label depending on which period it weekdays, weekends or holidays.
represents:
Indiana, Pavia Uni, Salinas: These three dataset of hyperspectral remote sensing datasets are all publicly available https:/ /wwww ehu. eus htts//cwintc//// index. phpppppppppa// http///ypeeppccteeeeee _er____emoe________ Rewote_eemsiiieeeeeeeeeeeeeeeeee Sensing Scenes. Each dataset contains one hyperspectral image of a certain region containing a unique number of re- flectance bands. We applied the same preprocessing that in Collas et al. (2021) (Bouchard et al., 2024) that consists in four main steps. First, we normalize the data by subtracting the image global mean. Then, we apply PCA to reduce the dimension of the data to 5. A sliding window with no overlap is then used around each pixel for data sampling and then vectorized. In our experiments, we used a window of size 25 x 25. Finally, we compute the covariance matrix of each vectorized window using the Sample Covariance Matrix to get a point cloud of covariance matrices. For each covariance matrix, its class was computed by taking the majority class of the pixels in the window.
23

---

Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices
This dataset is made of of real images from textile manufacturing that contain non-defective and defective woven textiles. These images come from the public MVTec Anomaly Detection dataset (Bergmann al., 2021). The same preprocessing as in Smith et al. (2023) was applied to get the covariance matrices. The two classes correspond to defective and non-defective textiles.
Textile:
BreizhCrops: This dataset is intended for supervised classification of field crops from satellite time series. We used the dataset FRH01 that composed of satellite time series from the French region Finistère. As they multivariate time-series, we simply converted them to covariance matrices using the Oracle Approximating Shrinkage estimator (Chen et al., 2010). The classes correspond to different types of crops. More details on this dataset can be found in the original paper (RuBwurm et al., 2020).
For the non-BCI datasets (AirQuality, Indiana, Pavia Uni, Salinas, Textile and BreizhCrops), used 5-fold cross volldaaaii to evaluate the performance of the classifiers.
Let us also give some details on the implementation of the different classifiers used in the experiments.
The MDM is implemented using the library PyRiemann (Barachant et al., 2024).
The TS-LDA uses the Tangent class from PyRiemann (Barachant et al., 2024) and the LDA from Scikit-learn (Pedregosa et al., 201
The TS-QDA uses the Tangent Space class from PyRiemann (Barachant et al., 2024) and the Naive Bayes classifier from Scikit-learn (Pedregosa et al., 2011).
For the Ho-WDA and He-WDA, we implemented them using our MLE estimate the parameters of the wrapped Gaussians. To optimize the MLE, we used in practice the Riemannian Conjugate Gradient method (Boumal, 2023) with a maximum of 1, 000, 000 iterations and a max time set to 2 hours.
24


Reviews:
[{"review_id": "dl8okbwdpA", "reviewer": "Reviewer_9rYP", "summary": "The authors proposed a wrapped Gaussian formalism and give ML estimator for mean and covariance. The authors showed usefulness of proposed formulation in LDA and QDA on several datasets.", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Regarding the “Essential References Not Discussed”\nThe reference [1] you mention is actually the first reference given in Section 2 devoted to related works, on the first page of our paper. It was a key reference for the development of our theory as the authors proposed an __isotropic__ Gaussian distribution Riemannian Symmetric Spaces (e.g. SPD). Our goal was to study a __non-isotropic__ Gaussian distribution on SPD matrices. This is clearly stated in the first 6 lines of Section 2. \n\n## Regarding the \"Other Strengths and weaknesses\":\n### Minor comments:\nRegarding the choice of geometry, not all metrics behave the same on the manifold of SPD matrices. We could have chosen a flat metric leading to a Euclidean space.\n### Major comments:\n1. Although we only used the word \"intrinsic\" twice in our paper, it may have confused our reader. There is a trade off between working with an anisotropic distribution and having a fully intrinsic distribution. Indeed, as we detail in Section 2, Pennec has defined in [2] a purely anisotropic and intrinsic Gaussian distribution on a Riemannian manifold. However, as detailed in section 2:\n\t- The normalizing constant is expressed by an integral over the whole manifold $P_d$ (intractable in practice).\n\t- It uses a concatenation matrix rather than a covariance matrix, and the relation between those two matrices requires the computation of an integral over $P_d$. \n\t- No sampling method has been studied for it. \n\tHence, we tried to take the best of both world in our work by tackling the anisotropy and to provide an easy-to-sample distribution, and we completed these objectives.\n2. The Gaussian distribution defined in [1] also relys (indirectly) on the choice of a particular tangent space. Indeed, in the p.d.f. of the Gaussian, the difference between a given point $y$ and the mean $\\bar{x}$ is computed using $Log_{\\bar{x}}y$ (differenc in the tangent space). However, this distribution is called “intrinsic”. \n3. Allowing the distribution to be non-centered in the tangent space raises several theoretical and practical issues that are not present when $\\mu = 0$:\n\t- The non-identifiability issue raise by the choice of $\\mu \\neq 0$ is solved in section 4.4. \n\t- As stated in remark 4.12, the equivalence classes of a given wrapped Gaussian (WG) does not necessarily contain a centered WG. Therefore, the choice of $\\mu$ has a real impact on the expressivity of our distribution. \n\t- The estimation of the parameters is another issue: if one assumes $\\mu = 0$, then the estimation of the parameters is straightforward (by the methods of moments as it has been done by Chevallier et al. (2022)). However, if $\\mu \\neq 0$, then the estimation of the parameters is more complex and requires the use of a full MLE.\n\t- We would like to recall that, unlike the Euclidean case, $\\mu$ does not model the mean of the distribution, as we introduce a new parameter $p$. Maybe you got confused between those two parameters.\n\n\tGiven these precisions, could you argue what is specifically trivial and why ? \n4. The CLT on the tangent space is not novel to prove. However, the CLT on the manifold was (to the best of our knowledge) never stated in the literature. We believe that this result is interesting as it shows the interest of considering WG on SPD matrices as a WG appears as the asymptotic distribution of a logarithmic product of i.i.d. random variables on the manifold. If this result is known, could you give a reference ?\n5. The estimation of the parameters of a WG is not a trivial extension of the Euclidean Gaussian case. As we explain in proposition 5.1, __in the special case where $p^\\star$ is known a priori__, the MLE of $\\mu$ and $\\Sigma$ are the same as the Euclidean case. However, in the general case, the MLE of $p$ does not have a closed-form solution and is dependent on the MLE of ($\\mu$, $\\Sigma$). Thus, the estimation of the three parameters $p$, $\\mu$ and $\\Sigma$ must rely on the full maximum likelihood optimization problem stated in section 5.1. Once again, could you specify what is trivial here ?\n6. Using WG for classification with LDA and QDA primarily served to illustrate its applications. Our goal was to demonstrate how the theory developed earlier enables the construction of machine learning tools directly on the manifold of SPD matrices. The main value of WG lies in offering a novel way to model data on SPD matrices. Additionally, in Section 6.2, we showed how classical SPD classifiers can be reinterpreted through the WG framework, providing a unifying probabilistic perspective. We believe that introducing such a perspective is valuable, and this paper represents a first step in that direction.\n\n[1] Said, S. et al, Gaussian Distributions on Riemannian Symmetric Spaces: Statistical Learning With Structured Covariance Matrices. IEEE Trans. Inf. Theory, 2018. \n\n[2] Pennec, X. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements. J. Math. Imaging Vis, 2006.", "claims_and_evidence": "1. My biggest concern regarding the intrinsic nature of the formalism as claimed, the formalism is essentially a extrinsic formulation on tangent space, define Gaussian on the tangent space and push-forward on the manifold", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "Yes, the claims hold.", "experimental_designs_or_analyses": "The experiments are sound, although not convincing how much benefit we get using wrapped Gaussian", "relation_to_broader_scientific_literature": "The ideas are relevant for the community.", "essential_references_not_discussed": "1. Gaussian distributions on Riemannian symmetric spaces", "other_strengths_and_weaknesses": "Minor comments:\n\nThe nations are uncommon, for example why use both vector and matrix using Lower case?!\nIn all metrics on SPD, it has non-positive curvature, not just affine-invariant metric: “Once\nendowed with this metric, Pd is a complete connected Riemannian manifold of non-positive curvature” — this statement can be clarified.Also all Riemannian metrics give same sign of curvature, as the metrics are equivalent in that sense.\n3. Section 3.1 reads dense and a laundry list of things, better to clarify saying why (if possible) (and where) you need this to use this. All of these are well-known and well-studied across geometers so the non-familiar reader should know in which section you will use these tools.\n4. In section 4.2, “pus-forward” should be “push-forward”.\n5. In section 4.2, “The Jacobian determinant” better call it determinant of the Jacobian.\n\nMajor Comments:\n\n1. Whenever one defines something on Eucliden space and pull over manifold the intrinsic nature is gone. The authors argued in Introduction how classical tools are non-intrinsic: “However, classical Eu-\nclidean probability distributions fail to capture the intrinsic geometry of the underlying manifold.” So I am confused why authors definition of wrapped Gaussian is non-intrinsic as in section 4.1.\n\n2. The entire argument of being intrinsic doesn’t hold as everything defined in relate of wrapped Gaussian is extrinsic including Therorem 4.2, proposition 4.3\n\n3. In section 4.2., the authors said they allow usage of \\mu, which is very trivial extension of centered Gaussian.\n\n4. Why in section 4.3, the CLT will be something novel to prove? As it depends on CLT on the tangent space. Am I missing something here?\n\n5. The estimation of ML estimators of mean and std follows from Euclidean Gaussians, so again section 5 follows naturally from Euclidean MLEs.\n\n6. The performance of LDA and QDA using wrapped Gaussian is not very superior over the competition as in Table 2. So not sure why we will use such a formalism even if we ignore the extrinsic nature of it.", "other_comments_or_suggestions": "Please see weakness section", "questions_for_authors": "Please see weakness.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "CkHTbiihPu", "reviewer": "Reviewer_nsaj", "summary": "The authors propose a new version of Gaussian on the SPD manifold (with the affine metric) by wrapping a Gaussian from a tangent space onto the manifold. Their method has two main differences than previously proposed methods 1. the distribution need not be isotropic and 2. the footpoint of the distribution on the manifold need not be the mean of the distribution. The authors further extend the ideas of statistical classification using their wrapped Gaussian distribution.\n\nThis review has been updated after the rebuttal period.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "First, we would like to thank the reviewer for their work, their valuable comments and interesting questions. \n\n## Regarding the Claims and Evidence:\nYou say that you have doubts about the claims on the MLE made in section 6.1. The only claim on the MLE made in this section is that the LDA uses an MLE on the training data to learn the parameters of each class. Are you referring to this claim? If it is the case, we refer to section 4.3 of [1] (page 109) where the formulas for the estimated parameters are given for the Euclidean LDA. On can check that they coincide with the formula of the MLE in the Euclidean Gaussian case.  \n## Regarding the Comments Or Suggestions:\n1. We will mention earlier in the introduction that we focus on the affine-invariant metric on the manifold of SPD matrices and mention that other choices could be made (for example Log-Euclidean metric).\n2. We will check and correct the different typos and errors throughout the paper. \n## Regarding the Questions For Authors: \n1. In section 5, $\\Sigma$ is never assumed to be diagonal. The experiments lead in the section were always led with a full $\\Sigma$. However, we mention that one can assume that $\\Sigma$ is diagonal to simply the estimation problem (reducing the number of parameters to estimate from $O(d^4)$ to $O(d^2)$ where $d$ is size of the SPD matrices). This assumption of a diagonal $\\Sigma$ is done in the experiments led in Section 6.1 and in Appendix J. We will clarify this point in the final version. \n2. The normalizing constant in the MDM does not depend on the mean of the distributions. Indeed, as shown in Proposition 1 of [2], the normalizing constant of the isotropic Gaussian distribution depends only on $\\sigma$, not on the mean.\n3. Absolutely, in the description of the Tangent Space LDA and QDA, one should read “logarithm map” instead of “exponential map”. We will correct this error. \n4. In section 5.1, the $p^2$ is not meant to be $p$ squared but the $^2$ refers to a footnote. This can indeed be confusing and will be corrected.  \n\nFinally, you mention paper [3] that completely fits in the scope of our paper, and we will make sure to add this reference to our section 2 on related works. \n\n[1] Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning. Springer Series in Statistics. Springer, New York, NY, 2009. \n\n[2] Said, S., Hajri, H., Bombrun, L., and Vemuri, B. C. Gaussian Distributions on Riemannian Symmetric Spaces: Statistical Learning With Structured Covariance Matrices. IEEE Transactions on Information Theory, 64 (2):752–772, February 2018.\n\n[3]  Hatem Hajri, Ioana Ilea, Salem Said, Lionel Bombrun, Yannick Berthoumieu. Riemannian Laplace distribution on the space of symmetric positive definite matrices. 2015", "claims_and_evidence": "Almost all their claims are supported. \n\nThe only claims which I am not entirely sure are on the section 6.1 claims of MLE. These claims do seem intuitive, but I think a more thorough justification is necessary.", "methods_and_evaluation_criteria": "Yes, the authors test their methods on both simulated data and real data. The simulated data is somewhat low dimensional with a max d of 10 however the real data is on the same order so it seems appropriate.", "theoretical_claims": "Only glanced at most proofs but took a longer look at the proof in section H of the appendix.", "experimental_designs_or_analyses": "Yes, all were checked.", "relation_to_broader_scientific_literature": "The SPD manifold is often used but the structure is often ignored as it can be difficult to handle, the authors propose a method on the manifold leveraging the structure. Many authors have considered other distributions on SPD such as Hajri et al. who proposed the Laplace and this paper seems like a direct comparison in terms of impact.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "See comments below.", "other_comments_or_suggestions": "In the introduction, I think the authors should mention they are considering the affine metric earlier as there are two natural options for a metric on this space.\n\nSecond page I think there is a hanging sentence \"the authors work on symmetric spaces.\" the sentence seems out of place or incomplete.\n\n\"In the sequel\" seems odd choice of wording.\n\n\n\nTypos and errors:\n4.1 \"tangent plan\" -> \"tangent plane\"\n4.2 starts with \"pus-forward\" \n5.1 the theta star parameters missing comma\n6.1 \"They rely on metric that arex pulled\" -> \"They rely on metrics that are pulled\"\nAll paper: Inconsistent notation of $WG(p,\\mu,\\Sigma$ v $WG(p;\\mu,\\Sigma)$\nFigure 3, the third panel has the star in the wrong location", "questions_for_authors": "In section 5 it's unclear when it is assumed $\\Sigma$ is diagonal. It seem in 5.1 that it will be always but then in 5.2 it seems as if only sometimes.\n\nIs the normalizing constant in MDM dependent on the mean on the distribution?\n\nWhen describing LDA and QDA the authors mention that \"all training points are sent to the tangent space via the exponential map\" should this not be the log map?\n\n\nIn section 5.1 is the $p^2$ a typo?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ugy0dPsmQ0", "reviewer": "Reviewer_PgG4", "summary": "This paper studies the non-isotropic wrapped Gaussian distribution on the manifold of positive definite (PD) matrices. Specifically, the authors derive theoretical properties of the non-isotropic wrapped Gaussian distribution and propose maximum likelihood estimators for its parameters. They also define an equivalence relation between the set of parameters of two wrapped Gaussian and resolve the non-identifiability issue of the wrapped Gaussian for PD matrices. Finally, the authors provide new interpretations to several known classifiers on PD matrices through the lens of wrapped Gaussian distributions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "First, we would like to thank the reviewer for their work, their valuable comments and interesting questions. \n\n## Regarding the \"Claims And Evidence\":\n\nWe agree that the He-WDA should be a special case of the Ho-WDA when the covariance matrices for each classe are the same. In order to evaluate the performance of the He-WDA and Ho-WDA, we used a 5-fold cross-validation on each dataset. This is classicly done in ML to evaluate pipelines. We were limited by the amount of data in some datasets. Moreover, as the He-WDA estimates one covariance matrix per class, it has more parameters to estimate and then requires a lot of data per class. If the number of samples per class is low, the estimation of the covariance matrix can be very noisy and the performance of the He-WDA can be worse than that of the Ho-WDA. For example, for the dataset Salinas, some classes have only a few hundred samples, which is not enough to estimate the covariance matrix accurately. The Ho-WDA, on the other hand, estimates a single covariance matrix for all the classes and is less sensitive to the number of samples per class. We will add this explanation in the final version of the paper.  \n## Regarding the \"Methods And Evaluation Criteria\":\n\nThe sampling procedure of a wrapped Gaussian is very simple and does not rely on any sophisticated rejection sampling algorithm. Indeed, to sample from $WG(p;\\mu,\\Sigma)$, one can simply sample a point $x$ from the Euclidean Gaussian $\\mathcal{N}(\\mu,\\Sigma)$ and then project the samples on the manifold $P_d$ by $Exp_p(Vect^{-1}_p(x))$. The exponential map $Exp_p$ as well as the vectorization $Vect_p$ are both simple to compute. We will add this algorithm in the final version of the paper. \n\n## Regarding the \"Other Strengths And Weaknesses\": \n\nWe would like to emphasize the fact that allowing the distribution to be non-centered in the tangent space is not just a “slightly more general setting''. For more details on this point, please refer to our answer to point 3 of reviewer 9rYP. \n\n## Regarding the \"Other comments or suggestions\":\nWe will correct the different typos highlighted in this review. \n\n3. If the SPD matrices are of size $d \\times d$ (like we consider in the paper), then $\\nu$ is the concatenation of $d$ ones and $d(d+1)/2-d = d(d-1)/2$ zeros. We will clarify this point in the final version.\n4. This is an important point. Let us reformulate remark 4.12: Let us consider a wrapped Gaussian $WG(p;\\mu, \\Sigma)$. Then, the equivalent wrapped Gaussians are of the form $WG(e^t p; \\mu + tVect_p(p), \\Sigma)$ for $t \\in R$. If $\\mu$ and $Vect_p(p)$ are aligned i.e., there exists $\\tilde{t}$ such that $\\mu = -\\tilde{t}Vect_p(p)$, then the equivalence class contains a wrapped Gaussian with $\\mu = 0$. However, as $\\mu$ and $Vect_p(p) = (1,\\cdots,1,0,\\cdots,0)$ (the concatenation of $d$ ones and $d(d-1)/2$ zeros, see previous point) are not aligned (for example, take $\\mu = \\nu + (1,\\cdots,0) = (2,1,\\cdots,1,0,\\cdots,0)$), then there exists no $t$ such that $\\mu = -tVect_p(p)$ and the equivalence class does not contain a wrapped Gaussian with $\\mu = 0$. We will clarify this point in the final version.\n5. We agree with you that we could expect that the problem with $d=10$ is more challenging than the problem with $d=5$. According to this, we would expect the error on the estimation of $p^\\star$ to be higher for $d=10$ than for $d=5$. However, the results of Figure 3 show the contrary. Maybe we need to repeat more than 5 times the experiment to have a more accurate estimation of the error. We will try to increase the number of repetitions in the final version of the paper.\n6. In order to center the data on the tangent space (to have $\\mu^\\star = 0$), one needs to know $p^\\star$ and $\\mu^\\star$ (more details on how to center the data are given in appendix D).  Indeed, we need to know the tangent plan in which we want to center the data. As here our goal is to estimate the parameters, we do not know $p^\\star$ and $\\mu^\\star$ and thus cannot center the data. Moreover, if one simply “centers'' the data using the Riemannian mean $\\mathfrak{G}(x_1,...,x_N)$ , one does not have any guarantees that the “centered” data will have $\\mu = 0$. \n8. Computing exponential and logarithmic maps remains a bottleneck in SPD matrix geometry. However, a trade-off may exist between computational cost and performance gains. Theoretically, Euclidean Gaussian-based methods extend to $P_d$ via our wrapped Gaussian, though practical challenges will arise in applications, requiring careful choices.\n9. We will modify the font size of the proofs to make them more readable.\n12. Indeed, the random variable $X$ is missing in the second point of Proposition D.1. We will correct this error.\n14. The $n$ in $||\\nu||_2^2 = nd$ is a typo. Since $\\nu$ consists of $d$ ones and $d(d-1)/2$ zeros, its squared norm is simply $d$. Therefore, the expression for $p_min$ is correct. This error will be fixed.", "claims_and_evidence": "Overall, all the claims and results in this paper are supported by rigorous proofs and/or simulation results. However, I am not quite convinced by the claim on the second column of Page 8 (Line 403) that \"we do not observe a clear dominance of the He-WDA over the Ho-WDA\". The major issue of this claim is that the Monte Carlo experiments are only repeated for 5 times, which is clearly not enough. It should be at least 100 times. Moreover, it is not intuitive about why the He-WDA behaves worse than the Ho-WDA on many of data examples. Shouldn't Ho-WDA be a special case of He-WDA when all the covariance matrices for different classes are the same.\n\nThere are also some minor issues in the paper that I pointed out below.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria basically makes sense for the wrapped Gaussian distribution problem and its related classification problems. As a side note, since the authors mentioned that they can sample from the wrapped Gaussian distribution, it would be better to outline the sampling procedure in the main paper or Appendix. Specifically, a procedure without rejection sampling is expected.", "theoretical_claims": "I have checked all the proofs and results in both the main paper and supplementary materials.", "experimental_designs_or_analyses": "Yes, I have checked the validity of experimental analyses. The only concern is the Monte Carlo repetition times that I mentioned above.", "relation_to_broader_scientific_literature": "To my knowledge, the wrapped distributions have been studied in directional statistics dating back to at least 1970s. This paper extends this wrapped technique for the Gaussian distribution, or more generally, elliptically contoured distributions, to the manifold of positive definite matrices. As the authors pointed out, wrapped distributions have been studied on homogeneous Riemannian manifolds, but using a different techniques as the one proposed in this paper. More relatedly, some exponential-wrapped distributions on symmetric spaces were studied in the literature as well, but as the authors in this paper pointed out, these related works consider the distribution on the tangent space to always be centered. In this paper, the authors consider a slightly more general setting with the wrapped Gaussian distribution on the manifold of positive definite matrices not necessarily being centered.", "essential_references_not_discussed": "To the best of my knowledge, the paper did a good job in discussing the related works. The only main concern is the novelty of this paper when compared with the prior works that consider the centered distribution on the tangent space. Intuitively, if we know the tangent space of a manifold, it does not seem to be very difficult to center the data and/or distribution with respect to the origin of the tangent space. I encourage the authors to address this concern in more detail. Appendix H could be one example, but still not convincing enough.", "other_strengths_and_weaknesses": "The writing and proofs of this paper are of good standard.\n\nThe only main concern is the novelty of this paper when compared with the prior works that consider the centered distribution on the tangent space. Intuitively, if we know the tangent space of a manifold, it does not seem to be very difficult to center the data and/or distribution with respect to the origin of the tangent space. I encourage the authors to address this concern in more detail. Appendix H could be one example, but still not convincing enough.", "other_comments_or_suggestions": "1. In the Abstract (Line 21-25), the sentence \"We introduced a non-isotropic wrapped Gaussian by leveraging the exponential map, we derive theoretical properties of this distribution and propose a maximum likelihood framework for parameter estimation.\" is not correct in grammar.\n\n2. Second line in Section 4.2: a typo \"pus-forward\".\n\n3. Proposition 4.9: it should be clearly stated how many ones are there in $\\nu=(1,...,1,0,...,0)$. The same applies to Proposition 4.14.\n\n4. I am confused about Remark 4.12: can we just set $\\mu_{\\alpha}= -t  \\mathrm{Vect}_{p{\\alpha}}(p{\\alpha})$ so that the equivalence class can contain $\\mu=0$?\n\n5. Figure 3: For the case $d=10$, why is it almost consistently better than the case when $d=5$ for MLE to estimate $p^*$? It seems that $d=10$ is a more challenging problem. Does it mean that the Riemannian conjugate gradient algorithm doesn't fit to this problem?\n\n6. The first column of Page 6 (line 308-310): can we center the data and then apply the method of moment? what difficulties are there for preventing this straightforward adaptation?\n\n7. The first column of Page 7 (Line 366): typo \"arex pulled back..\"\n\n8. In the conclusion part, it seems to me that it may not be feasible to extend all the classical machine learning models that rely on Gaussian distributions to the manifold of SPD matrices. The computational issues incurred by calculating the exponential and logarithmic maps are huge barriers.\n\n9. Why did the authors present all the proofs in tiny font? It is hard for readers to read them.\n\n10. Line 793 on Page 15: typo: \"does dot imply...\".\n\n11. Line 812: \"$Id$\" at the end of the equation should be $I_d$.\n\n12. Second point in Proposition D.1: the random variable $X$ is missing. It should be $\\mathrm{Log}_p X$.\n\n13. Line 860 on Page 16: typo \"diffeomorphisme\".\n\n14. Line 1061 on Page 20: If $||\\nu||_2^2 = n d$, \n\nthen where does the \"n\" go in the expression of $p_{\\min}$?", "questions_for_authors": "Please see my concerns and comments above. I am happy to increase my scores if the authors can carefully address my comments and concerns.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper is borderline, with a split between reviewers, all of whom are seasoned experts on this topic. 

Those that liked the paper praise its results, particularly that the SPD manifold occurs often and is difficult to handle, indicating a need to develop better theoretical tools for working with it. Reviewers also felt the paper covered the literature well, and there were no complaints about presentational issues in spite of the fact that such complaints are extremely common.

The reviewer who was more critical, to my understanding, was primarily concerned with whether it is valid or not to describe the proposed method as intrinsic. I am not convinced this criticism matters to a decisive degree: whether a method is formulated in an extrinsic or intrinsic manner is secondary to whether or not it is a useful tool for practitioners to use on their problems. This reviewer also did not respond to a request for comment in the discussion phase.

On basis of this comparison, I believe this work's merits outweigh its downsides, and recommend a weak accept.
Decision: Accept (poster)