I'm trying to sample uniformly from a ball around the identity matrix in the manifold of symmetric positive-definite matrices (SPD), i.e., $\mathcal{B}(R) = \{X \in \mathcal{P}(n) : d_\mathcal{P}(I_n, X) \le R\}$, where $\mathcal{P}(n)$ is endowed with the usual metric, $ds = \lVert X^{-1} dX \rVert_F$. I'm getting some unexpected results and the question is what exactly is the issue with the approach described below (if any)?
I've recently learned (from a previous question) that this can be achieved by sampling from a density $p(X) \propto \sqrt{\det g(X)}$ where $g$ is the Riemannian metric in some normal coordinates. (I'm using the suggested rejection sampling approach for this.)
Additionally, I've learned about the following identity concerning the Riemannian metric of the SPD manifold (see, e.g., [ref, Prop.1]),
$$\det(g(X)) = 2^{n (n - 1) / 2} \det(X)^{n + 1}.$$
Now, given that this is a homogeneous space, I'm restricting to sampling around the identity matrix. Then, each matrix $X$ can be represented in the exponential chart at $I_n$ as $P = \log X$, with $P$ a symmetric $n \times n$ matrix. The density that I want to sample from is then,
$$p(X) \propto \exp\Big(\frac{n+1}{2} \mathop{\textrm{Tr}} \log X\Big).$$
The surprising result that makes me check the correctness of my approach is that the pairwise distances between the sampled matrices are mostly less than the radius $R$:
The distance is given by $d_\mathcal{P}^2(X, Y) = \lVert \log X^{-1/2} Y X^{-1/2} \rVert_F^2 = \sum_{i=1}^n \log^2 \lambda_i(X^{-1} Y)$, where $\lambda_i(P)$ denotes the $i$th eigenvalue of the matrix $P$ (here and below).
From the above formulas, I can see why this happens: being at most $R$-away from $I_n$ means $$ d_\mathcal{P}^2(I_n, X) = \sum_{i=1}^n \log^2 \lambda_i(X) = \sum_{i=1}^n \log^2 \exp(\lambda_i(\log X)) = \sum_{i=1}^n \lambda_i^2(\log X) \le R, $$ while large density mass is assigned to matrices with large values of $$ \mathop{\textrm{Tr}} \log X = \sum_{i=1}^n \lambda_i(\log X). $$
This seems analogous to sampling uniformly from the unit disk/ball and then weighing the positive quadrant/octant (exponentially) more -- only that in the spectrum space. Is this really the case or am I misapplying or misinterpreting some of the references from above?
Thank you.
[ref]: Moakher, Maher, and Mourad Zéraï. "The Riemannian geometry of the space of positive-definite matrices and its application to the regularization of positive-definite matrix-valued data." Journal of Mathematical Imaging and Vision 40.2 (2011): 171-187.
