Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices

First, we would like to thank the reviewer for their work, their valuable comments and interesting questions.

Regarding the "Claims And Evidence":

We 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.

Regarding the "Methods And Evaluation Criteria":

The 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. 

Regarding the "Other Strengths And Weaknesses": 

We 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. 

Regarding the "Other comments or suggestions":

We will correct the different typos highlighted in this review.

3. 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.
4. 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.
5. 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.
6. 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$. 
7. 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.
8. We will modify the font size of the proofs to make them more readable.
9. Indeed, the random variable $X$ is missing in the second point of Proposition D.1. We will correct this error.
10. 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.

First, we would like to thank the reviewer for their work, their valuable comments and interesting questions.

Regarding the Claims and Evidence:

You 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.

Regarding the Comments Or Suggestions:

1. 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).
2. We will check and correct the different typos and errors throughout the paper.

Regarding the Questions For Authors:

1. 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.
2. 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.
3. 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.
4. 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.

Finally, 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.

[1] Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning. Springer Series in Statistics. Springer, New York, NY, 2009.

[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.

[3] Hatem Hajri, Ioana Ilea, Salem Said, Lionel Bombrun, Yannick Berthoumieu. Riemannian Laplace distribution on the space of symmetric positive definite matrices. 2015

Regarding the “Essential References Not Discussed”

The 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.

Regarding the "Other Strengths and weaknesses":

Minor comments:

Regarding 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.

Major comments:

1. 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:

   • The normalizing constant is expressed by an integral over the whole manifold $P_d$ (intractable in practice).
   • 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$.
   • No sampling method has been studied for it. Hence, 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.

2. 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”.

3. Allowing the distribution to be non-centered in the tangent space raises several theoretical and practical issues that are not present when $\mu = 0$:

   • The non-identifiability issue raise by the choice of $\mu \neq 0$ is solved in section 4.4.
   • 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.
   • 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.
   • 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.

   Given these precisions, could you argue what is specifically trivial and why ?

4. 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 ?

5. 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 ?

6. 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.

[1] Said, S. et al, Gaussian Distributions on Riemannian Symmetric Spaces: Statistical Learning With Structured Covariance Matrices. IEEE Trans. Inf. Theory, 2018.

[2] Pennec, X. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements. J. Math. Imaging Vis, 2006.