# Curvature-based rejection sampling

Isabella Costa Maia, Marco Congedo, Pedro Luiz Coelho Rodrigues, Salem Said

► To cite this version:

Isabella Costa Maia, Marco Congedo, Pedro Luiz Coelho Rodrigues, Salem Said. Curvature-based rejection sampling. 2025. hal-05332278v3

HAL Id: hal-05332278

https://hal.science/hal-05332278v3

Preprint submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Curvature-based rejection sampling

**Isabella Costa Maia**

GIPSA-lab, University Grenoble Alpes, CNRS, Grenoble-INP

isabella.costa-maia@grenoble-inp.fr

**Marco Congedo**

GIPSA-lab, University Grenoble Alpes, CNRS, Grenoble-INP

marco.congedo@gipsa-lab.grenoble-inp.fr

**Pedro L. C. Rodrigues**

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK

pedro.rodrigues@univ-grenoble-alpes.fr

**Salem Said**

Univ. Grenoble Alpes, CNRS, Grenoble INP, LJK

salem.said@univ-grenoble-alpes.fr

## Abstract

The present work introduces curvature-based rejection sampling (CURS). This is a method for sampling from a general class of probability densities defined on Riemannian manifolds. It can be used to sample from any probability density which “depends only on distance”. The idea is to combine the statistical principle of rejection sampling with the geometric principle of volume comparison. CURS is an exact sampling method and (assuming the underlying Riemannian manifold satisfies certain technical conditions) it has a particularly moderate computational cost. The aim of the present work is to show that there are many applications where CURS should be the user’s method of choice for dealing with relatively low-dimensional scenarios.

# 1 Introduction

The present work introduces a new method for sampling from probability densities defined on Riemannian manifolds. This method is called curvature-based rejection sampling (CURS). It addresses the problem of sampling from a probability density function which “depends only on distance”.

Specifically, let $M$ be a Riemannian manifold and let $x_o$ be some point in $M$. CURS can be used to sample from any probability density $p$ on $M$, of the form

$$p(x) \propto f(r(x)) \quad \text{where } r(x) = d(x_o, x) \qquad (1)$$

Here, $p \propto f$ means $p$ is proportional to $f$, and $d(x_o, x)$ denotes the Riemannian distance between $x_o$ and $x$ in $M$ (the precise definition of a probability density on $M$ is given in the following Section 2).

CURS is a geometric rejection sampling method which has the following advantages:

* it is a rather general method, since it applies to any density $p$ of the form (1) and to a large class of manifolds $M$. For example, it applies successfully whenever $M$ is a so-called Hadamard manifold, such as a manifold of covariance matrices (real, complex, Toeplitz, etc.) (Lee, 2018; Said et al., 2018), but is by no means limited to this special case of Hadamard manifolds.

* it is an exact method, which produces samples that immediately follow the density $p$, with no need for any kind of burn-in period. This is in contrast to existing Markov-Chain Monte Carlo methods.

* the computational cost of a single iteration of this method is quite moderate.

* it is often possible to theoretically evaluate the time complexity of this method, by computing its rejection probability. In this way, one has a theoretical performance guarantee.
The introduction of CURS was motivated by the specific problem of sampling from a Gaussian density

$$p(x) \propto \exp\left[-\frac{d^2(x_o, x)}{2\sigma^2}\right] \quad \text{where } \sigma > 0 \qquad (2)$$

which is defined on the manifold M of real covariance matrices. In the notation of (1), this corresponds to $f(r) = \exp(-r^2/2\sigma^2)$. Section 5, below, will return to this problem.

At present, existing methods for sampling on manifolds can be divided into exact and approximate methods. Exact methods tend to be specialized, applicable only to specific probability distributions and specific classes of manifolds. For example, there exist highly effective methods for sampling from uniform distributions on compact Lie groups, or from classical directional distributions on Stiefel and Grassmann manifolds (Meckes, 2019; Chikuse, 2003). Clearly, the major drawback of such methods is that they are too specialized, as they are of no immediate use beyond their originally intended goal.

Approximate methods are based on geometric Markov Chain Monte Carlo techniques (geometric MCMC). They are very general, as they can be used to sample from any positive density, on any Riemannian manifold which satisfies rather mild technical conditions. Famous examples of such methods include Hamiltonian Monte Carlo and geometric Langevin MCMC (Girolami & Calderhead, 2011; Bharath et al., 2025). Their main drawback is their computational complexity, due to the use of often complicated discretization schemes.

CURS follows a kind of middle ground which allows it to make the best of both worlds. Being neither too specialized nor completely general, it is able to perform exact sampling at a moderate computational cost. Its main drawback is that, being a rejection sampling method, its performance is significantly affected by the curse of dimensionality.

Still, users should find it useful to know about CURS, and retain it as part of their toolbox, due to its highly advantageous performance-complexity tradeoff in relatively low-dimensional scenarios.

The present work is organised as follows. Section 2 explains the general idea behind the CURS method. Section 3 formalises this idea into pseudo-code instructions. Section 4 discusses the geometric interpretation of CURS, while Section 5 applies it to the problem of sampling from a Gaussian density of the form (2). Finally, Section 6 is concerned with a more general version of CURS, extending its applicability to a larger class of Riemannian manifolds.

# 2 The idea behind CURS

## 2.1 Geodesic spherical coordinates

Consider the problem of sampling from a probability density p of the form (1). Since p(x) depends only on the distance r(x), one may hope to perform a separation of variables, by introducing geodesic spherical coordinates (Chavel, 2006; Lee, 2018).

These coordinates parameterise the point $x \in M$ using two quantities, the distance r and the direction s (in fact, s is a unit-length tangent vector to M at the point $x_o$). Then, x is obtained by moving for a distance r along a geodesic curve starting at $x_o$ and whose initial velocity is equal to the vector s. In terms of the Riemannian exponential mapping (Chavel, 2006; Lee, 2018)

$$x = \operatorname{Exp}_{x_o}(rs) \qquad (3)$$

which will be abbreviated to $x = x(r, s)$.

In general, geodesic spherical coordinates do not regularly cover the entire manifold M, but are limited to a subset $D_{x_o}$ called the domain of injectivity. For now (to simplify the presentation) it is assumed that $D_{x_o} = M$, so the mapping $(r, s) \mapsto x(r, s)$ in (3) is a diffeomorphism from $(0, \infty) \times S_{x_o} M$ to $M - \{x_o\}$ (here, $S_{x_o} M$ denotes the set of unit-length tangent vectors to M at $x_o$). This assumption is certainly verified when M is a Hadamard manifold, such as a manifold of covariance matrices (real, complex, Toeplitz, etc.) (Lee, 2018; Said et al., 2018). It will be dropped later on, in Section 6.
One hopes to sample $x$ by first separately sampling $r$ and $s$ and then replacing them into (3). The major stumbling block for this approach will be the difficulty in expressing the Riemannian volume measure of $M$. To say that $p$ is a probability density on $M$ means that the probability distribution of a random point $x$ (when $x$ is sampled from $p$) is given by

$$P(dx) = p(x) \mathrm{vol}(dx) \qquad (4)$$

where vol denotes the Riemannian volume measure of $M$. From (1), the joint distribution of $r$ and $s$ is then

$$P(dr \times ds) \propto f(r) \mathrm{vol}(dr \times ds) \qquad (5)$$

where vol(dr × ds) is just vol(dx), but expressed in terms of $r$ and $s$. This is given as follows.

Note that $S_{x_o} M$ is the unit sphere in the tangent space to $M$ at $x_o$. Let $\omega$ denote the surface area measure on this sphere. Then, the following formula holds (Chavel, 2006)

$$\mathrm{vol}(dr \times ds) = |\det A(r, s)| dr\omega(ds) \qquad (6)$$

where the matrix $A(r, s)$ is obtained by solving a second-order linear differential equation, known as the matrix Jacobi equation (Chavel, 2006) (Page 114). Replacing (6) into (5) yields

$$P(dr \times ds) \propto f(r) |\det A(r, s)| dr\omega(ds) \qquad (7)$$

Thus, in order to sample $r$ and $s$, one must sample from the joint distribution (7), and this is where the difficulties come in. Since $|\det A(r, s)|$ depends on both $r$ and $s$, these two cannot be sampled separately (in other words, independently). Rejection sampling comes to the rescue, as seen in the following subsection.

## 2.2 Introducing rejection sampling

Looking back to (7), the following observation can be made. In order to sample $r$ and $s$ independently, one should aim to have

$$|\det A(r, s)| = \text{a function of } r \times \text{a function of } s$$

Unfortunately, such a factorisation holds only quite rarely. However, in many situations, one can still replace the equality by an inequality. This is thanks to a classical theorem of Riemannian geometry, Bishop's volume comparison theorem (Chavel, 2006) (Page 130).

This theorem states that, if the sectional curvatures of $M$ are bounded below by some negative number $-\kappa^2$, then the following inequality holds

$$|\det A(r, s)| \le (\kappa^{-1} \sinh(\kappa r))^{d-1} \qquad (8)$$

where $d$ is the dimension of $M$. Now, the right-hand side is indeed the product of a function of $r$ by a function of $s$ (the latter being a constant function, equal to 1).

The condition required for (8) (lower bound on sectional curvatures) is verified in many, and even most, realistic situations. In particular, it holds for spaces of real or complex covariance matrices, with $\kappa = 1/\sqrt{2}$ (see Appendix B).

With (8) at hand, it becomes possible to apply a rejection sampling strategy, along the general lines explained in (Robert & Casella, 2004). This leads to the following formulation of the CURS method.

The target joint density of $r$ and $s$ in (7) reads

$$p(r, s) = \frac{1}{Z} \times f(r) |\det A(r, s)| \qquad (9)$$

where $Z$ is a normalising constant (the density is with respect to the reference measure $dr\omega(ds)$). Consider the proposal joint density (of course, with respect to the same reference measure)

$$g(r, s) = \frac{1}{Z_{\kappa}} \times f(r) (\kappa^{-1} \sinh(\kappa r))^{d-1} \qquad (10)$$
where $Z_κ$ is again a normalising constant. Then, note from (8)

$$p(r, s) \le T \times g(r, s) \text{ where } T = \frac{Z_κ}{Z} \qquad (11)$$

This is the key inequality required for rejection sampling. According to (Robert & Casella, 2004), to sample from $p(r, s)$, it is enough to sample from the proposal $g(r, s)$ and then reject any samples which do not verify

$$U \le \frac{1}{T} \times \frac{p(r, s)}{g(r, s)} \qquad (12)$$

where $U$ is an independently generated random variable, with uniform distribution in the unit interval $[0, 1]$. As explained in the previous Section 2, samples $x$ from the original density $p(x)$ are then obtained by replacing $r$ and $s$ into (3).

CURS is just a direct implementation of this idea. Everything works out because the proposal density $g(r, s)$ does not depend on $s$. This means that $r$ and $s$ can be sampled independently under this density. Moreover, it implies that $s$ follows a uniform distribution with respect to the surface area measure $\omega$ on the sphere $S_{x_o}M$. This makes it particularly easy to sample $s$.

## 3 CURS in pseudo-code

It is now possible to formulate the CURS method in pseudo-code. Its major steps turn out to be the following. A few comments are in order.

<caption>**Algorithm 1** CURS algorithm</caption>

```c
## Sample (r, s) from g(r, s) in (10)
1: Sample s uniformly on the sphere $S_{x_o}M$
2: Sample $r > 0$ from density
   g(r) ~ f(r) (κ⁻¹ sinh(κr))⁽⁻¹⁾
3: Apply the rejection procedure
4: Sample $U \sim U(0, 1)$.
5: if condition (14) holds then
6:   Reject (r, s) and return to step 2.
7: else
8:   Accept (r, s) and put $x = x(r, s)$ as in (3).
```

For step 1, sampling a uniform distribution on a sphere is a well-known procedure. First, $s$ is sampled from a standard normal distribution in the tangent space to $M$ at $x_o$. Then, $s$ is replaced with $s/||s||$ which clearly belongs to the unit sphere $S_{x_o}M$ ($||s||$ denotes the norm of $s$ with respect to the Riemannian metric of $M$). For this step, the vector $s$ should always be expressed in an orthonormal basis of the tangent space at $x_o$.

For step 2, if $f(r)$ is log-concave, then the density $g(r)$ in (13) will also be log-concave. It can then be efficiently sampled using the universal black-box algorithm of (Devroye, 1986).

For step 4, note that the condition in (14) is equivalent to (12), according to the definition of $M$ in (11). This step requires a single evaluation of $\det A(r, s)$. Oftentimes, this is known in closed form (see Section 5).

Eventually, if this is not the case, accurate numerical evaluation of $A(r, s)$ remains a straightforward option. This is not discussed here, since it is felt that situations where $\det A(r, s)$ is known in closed-form already cover an extensive range of applications.
The main indicator of the time complexity of CURS is the acceptance probability

$$ \Pi = \frac{1}{T} = \frac{Z}{Z_{\kappa}} \qquad (15) $$

In fact, $T$ is the expected number of iterations necessary to produce a single new sample. In order to have a theoretical performance guarantee for CURS, one should be capable of computing $\Pi$. While $Z_{\kappa}$ essentially reduces to a one-dimensional integral, the target density normalising constant $Z$ will be a multiple integral. The following expressions are completely general (see Appendix A)

$$ Z = \int_{S_{x_o} M} \int_0^\infty f(r) |\det A(r, s)| dr \omega(ds) \qquad (16) $$

$$ Z_{\kappa} = \Omega_{d-1} \int_{0}^{\infty} f(r) (\kappa^{-1} \sinh(\kappa r))^{d-1} dr \qquad (17) $$

where $\Omega_{d-1}$ is the surface area of a $(d-1)$-dimensional sphere (in other words, of $S_{x_o} M$). Clearly, the core problem related to computing $\Pi$ is the evaluation of $Z$.

# 4 Geometric interpretation

CURS is born from the marriage of the statistical principle of rejection sampling to the geometric principle of volume comparison. This confers a new geometric meaning onto the proposal density (10) and the rejection condition (14). In order to uncover this new meaning, recall Bishop's volume comparison inequality (8). This inequality becomes an equality in exactly one case: when $M$ is a hyperbolic space, a space of constant negative curvature equal to $-\kappa^2$.

If $M$ were in fact such a hyperbolic space, the target joint density (9) would be identical to the proposal density (10). In other words, the proposal density pretends that the space $M$ is a hyperbolic space with constant curvature $-\kappa^2$. The role of the rejection condition (14) is then to recover the true geometry of $M$.

The main difference between the hyperbolic “proposal geometry” and the true geometry of $M$ comes from anisotropy. In the hyperbolic case, a geodesic starting from $x_o$ in any initial direction $s$ only meets one value of the sectional curvature, $-\kappa^2$. In the general case, different directions $s$ can meet different values of the sectional curvature, and the true geometry is then anisotropic.

The rejection condition re-introduces the anisotropy which was left out from the proposal density, by favoring (accepting with higher probability) those samples which fall along certain directions. By assumption, the sectional curvatures of $M$ are bounded below by $-\kappa^2$. The rejection condition favors directions $s$ which meet curvatures closer to $-\kappa^2$ and tends to discard those directions with curvatures closer to positive values (incidentally, this justifies the name curvature-based rejection sampling).

Underlying this behavior is the following mathematical fact. If the curvatures along a certain direction $s$ are closer to $-\kappa^2$, then (along that direction) the volume comparison inequality (8) becomes sharper. In turn, the right-hand side of (14) becomes larger (closer to 1) and acceptance more probable. This is due to a complementary volume comparison inequality, called Günther's inequality (Chavel, 2006) (Page 128), which implies that, if the sectional curvatures along the direction $s$ are bounded above by some negative $-\delta^2$,

$$ |\det A(r, s)| \geq (\delta^{-1} \sinh(\delta r))^{d-1} \qquad (18) $$

It is precisely this lower bound which controls the sharpness of the upper bound in (8). The rejection condition favors directions which admit a $\delta$ closer to $\kappa$, because this makes (8) closer to being an equality, as it were.

In short, samples generated from the proposal density of CURS explore all directions with equal probability. However, directions which display “atypical curvatures” are then gotten rid of by the rejection condition. On the other hand, it should be clear that the probability of rejection increases with the distance $r$, for samples which fall along the same direction.

A concrete example of the above geometric interpretation will be provided in the following section (in particular, see the discussion after (27)).
# 5 Generalised Gaussian distributions

Recall the problem of sampling from the Gaussian density (2) on the space $M$ of real covariance matrices. Here, CURS is applied to the more general class of densities

$$p(x) \propto \exp\left[-\frac{d^{\alpha}(x_o, x)}{2\sigma^2}\right] \qquad (19)$$

defined on the space $M$ of $N \times N$ real covariance matrices, for $\alpha > 1$ and $\sigma > 0$. In fact, two versions of CURS will be considered. First, the general version outlined in the Section 3, and then an improved version, more specific to the problem at hand.

Everything in the present section extends immediately, from the particular case where $M$ is the space of real covariance matrices to the general case where $M$ is any negatively-curved symmetric space (see Appendix D). In particular, this general case includes various spaces of covariance matrices (real, complex, Toeplitz, etc.), which are useful for an extensive range of applications (Said et al., 2018).

With regard to the discussion at the end of Section 3, note that $\det A(r, s)$ is here known in closed form.

## 5.1 General CURS

The space $M$ is here equipped with its well-known affine-invariant Riemannian metric (Pennec et al., 2006),

$$\langle u, v \rangle_x = \operatorname{tr}(x^{-1}ux^{-1}v) \qquad (20)$$

for tangent vectors $u$ and $v$ to $M$ at $x$ (these are identified with real symmetric matrices, in the above formula). Therefore, one has the following Riemannian distance (Pennec et al., 2006)

$$d^2(x, y) = \operatorname{tr}(\log(x^{-1/2} y x^{-1/2}))^2 \qquad (21)$$

where $\operatorname{tr}$ denotes the trace and $\log$ the symmetric matrix logarithm. In order to spell out the main steps of Algorithm 1, an additional bit of information is needed. Namely, this is (see Appendix B)

$$\det A(r, s) = r^{N-1} \prod_{i<j} \left( \frac{\sinh(\kappa_{ij}(s)r)}{\kappa_{ij}(s)} \right) \qquad (22)$$

for the volume density in (6). Here,

$$\kappa_{ij}(s) = (\zeta_i - \zeta_j)/2 \quad 1 \le i < j \le N \qquad (23)$$

where $(\zeta_i)$ denote the eigenvalues of the symmetric matrix $s$. With this in mind, consider the following.

**step 1:** without any loss of generality, one may assume that $x_o = \text{id}$, the $N \times N$ identity matrix. Then, the metric (20) becomes the usual trace metric (i.e. scalar product)

$$\langle u, v \rangle_{x_o} = \operatorname{tr}(uv) \qquad (24)$$

An orthonormal basis for this metric is afforded by the matrices $u_{ij}$ ($1 \le i \le j \le N$),

$$\begin{aligned} u_{ij} &= (e_{ij} + e_{ji})/\sqrt{2} & \text{if } i < j \\ &= e_{jj} & \text{if } i = j \end{aligned}$$

where $e_{ij}$ is a matrix all of whose entries are zero, except the entry at row $i$ and column $j$, which is equal to 1. Accordingly, to carry out step 1, let $s_{ij}$ be independent standard normal real random variables and introduce

$$s = \sum_{i \le j} s_{ij} u_{ij} \qquad (25)$$

Then, replace $s$ with $s/\|s\|$, where $\|s\|^2 = \operatorname{tr}(s^2)$ since $\|s\|$ is the norm of $s$ with respect to (24).
Note here that it is possible to construct s as in (25) without explicitly dealing with the matrices ($u_{ij}$). If $t$ is an $N \times N$ matrix with independent elements, each having standard normal distribution, it is enough to set $s = (t+t^\dagger)/2$, where $\dagger$ denotes transposition (in the language of random matrix theory, $s$ is then distributed according to the so-called Gaussian orthogonal ensemble (Livan et al., 2018)).

**step 2:** $r$ must be sampled from the density

$$g(r) \propto \exp(-r^\alpha/2\sigma^2) (\kappa^{-1} \sinh(\kappa r))^{d-1} \quad (26)$$

where $\kappa = 1/\sqrt{2}$ and $d = N(N+1)/2$, the dimension of $M$. The fact that one may take $\kappa = 1/\sqrt{2}$ was mentioned in Section 2.2 and is proved in Appendix B.

Since $\alpha > 1$, the density (26) is a log-concave univariate density, and can be sampled using the universal black-box algorithm of (Devroye, 1986).

**step 4:** using (22) and also the fact that $\kappa = 1/\sqrt{2}$, it is straightforward to evaluate the condition (14). This becomes

$$U > (\kappa r / \sinh(\kappa r))^{N-1} \prod_{i<j} \frac{(\sinh(\kappa_{ij}(s)r) / \kappa_{ij}(s))}{(\sinh(\kappa r) / \kappa)} \quad (27)$$

Recalling (23), it can be seen that this condition tends to reject those $s$ whose eigenvalues ($\zeta_i$) are too close to one another, and that rejection becomes more likely as $r$ increases. This is in agreement with the geometric interpretation of Section 4 (rejection of directions $s$ with curvatures closer to positive values), as can be seen from the principal curvatures formula (49), which is discussed in Appendix B.

Some care is required with the right-hand side of (27), as this involves the product of a potentially large number (equal to $d-1$) of small quantities. In practice, it is better to evaluate the condition after taking logarithms of both sides.

**step 8:** when $x_o = \text{id}$, the Riemannian exponential mapping $\text{Exp}_{x_o}$ in (3) is just the matrix exponential. Finally, an additional step (step 9, then) is needed if $x_o \neq \text{id}$. This re-centres samples $x$ around the new $x_o$ using the transformation $x \mapsto x_o^{1/2} x x_o^{1/2}$.

When $\alpha = 2$, the acceptance probability $\Pi$ in (15) can be computed analytically, thanks to the expressions from (Said et al., 2023; Santilli & Tierz, 2021). Here, $Z$ and $Z_\kappa$ are denoted $Z(\sigma)$ and $Z_\kappa(\sigma)$, in order to highlight their dependence on $\sigma$. Thus,

$$Z(\sigma) = (\pi\sigma^2/2)^{N/2} 2^{N(N-1)/4} \prod_{j=1}^{N} \Omega_{j-1} \times \exp[-N(N-1)^2(\sigma^2/8)] \times \text{Pf}(\Lambda(\sigma)) \quad (28)$$

where $\Omega_{j-1}$ is the surface area of a $(j-1)$-dimensional sphere, and $\text{Pf}$ denotes the Pfaffian, while $\Lambda(\sigma)$ is the $N \times N$ matrix whose elements are given by

$$\Lambda_{ij}(\sigma) = \exp[(i^2 + j^2)(\sigma^2/2)] \text{erf}((j - i)(\sigma/2))$$

for $0 \le i, j \le N-1$, with $\text{erf}$ the error function (the definition of the Pfaffian can be recalled from (Livan et al., 2018)). On the other hand,

$$Z_{\kappa}(\sigma) = \left( \frac{\sigma \Omega_{d-1}}{(2\kappa)^{d-1}} \right) \sum_{j=0}^{d-1} (-1)^j C_j^{d-1} \frac{\Phi((d-1-2j)\kappa\sigma)}{\Phi'((d-1-2j)\kappa\sigma)} \quad (29)$$

where $C_j^{d-1}$ is a binomial coefficient and $\Phi$ is the standard normal cumulative distribution function. It is important to note that (28) only holds when $N$ is even (there is a similar expression that can be used for odd $N$ (Santilli & Tierz, 2021)).

Thanks to Expressions (28) and (29), it is possible to match the empirical acceptance probability $\hat{\Pi}$, observed from multiple runs of CURS, with its theoretical counterpart $\Pi$, computed from (15). Agreement between theory (II) and observation ($\hat{\Pi}$) is clear from the following Figures 1 and 2.
Figure 1 shows the plot of the theoretical acceptance probability $\Pi$ as a function of $\sigma$ (on the horizontal axis). This plot is overlaid with individual values of the empirical probability $\hat{\Pi}$, with each value estimated from $10^6$ iterations of the CURS algorithm (Algorithm 1). These are marked by dot points, which show no visible deviation away from the theoretical plot.

The acceptance probability clearly decreases with $\sigma$ and the rate of this decrease becomes more rapid as $N$ (the matrix dimension) increases. This is a manifestation of the curse of dimensionality. To put this effect into perspective, Figure 2 shows the dependence of $\Pi$ and $\hat{\Pi}$ on a new variable

$$ \delta = \int_M d^2(x_o, x) p(x) \mathrm{vol}(dx) \qquad (30) $$

the expected squared Riemannian distance away from $x_o$, with respect to the target density $p(x)$ in (19), with $\alpha = 2$. It is more intuitive to think in terms of $\delta$, since this new variable fits our usual idea of “variance”.

For example, Figure 2b shows that when $N = 4$, for a variance $\delta = 2$, the acceptance probability is $\Pi \approx 0.27$. Thus, on average, to obtain one new sample from the target density, 3.7 iterations of CURS are necessary. The actual time for a single iteration is of the order of $10^{-2}$ seconds (on a standard desktop computer).

By way of comparison, note that sampling from the same target distribution in (19) (as before, with $\alpha = 2$) was considered in (Bharath et al., 2025). Using a geometric Langevin algorithm, in the case where $N = 3$, this work reported having to compute upwards of $10^6$ trajectories of a discretised Langevin diffusion, in order to obtain a good approximation of the target density (see Table 3, in (Bharath et al., 2025)).

This requires a computational effort far greater than the one involved in CURS, since a single trajectory may contain hundreds of intermediate points, and computing each point requires evaluating both a matrix logarithm and a matrix exponential (see Algorithm 2, in (Bharath et al., 2025)). By contrast, one iteration of CURS evaluates only a matrix exponential, and this only once.

<figure><img src="image_2.png" /><figcaption>Figure 1: Theoretical (solid line) and empirical (dot points) acceptance probabilities, plotted against $\sigma$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Theoretical (solid line) and empirical (dot points) acceptance probabilities, plotted against $\delta$</figcaption></figure>
In all generality, one should also recall CURS is an exact sampling method, whereas the Langevin algorithm (being an MCMC method) is only approximate. To produce new samples, the Langevin algorithm either has to compute new trajectories, or to run a single trajectory for a sufficiently long time.

## 5.2 Sharp CURS

The acceptance probability observed in the previous subsection suffers from a significant drop as the matrix dimension $N$ increases. While this is to be expected, as rejection sampling methods tend to suffer from the curse of dimensionality, it can still be partially mitigated by employing a new, “sharper” version of CURS.

This new version applies the same steps as in the previous subsection, but replaces (27) with the condition

$$U > \prod_{i<j} \frac{(\sinh(\kappa_{ij}(s)r)/\kappa_{ij}(s))}{(\sinh(\kappa r)/\kappa)} \quad (31)$$

When looked at closely, this is the same as (27), but without the factor $(\kappa r/\sinh(\kappa r))^{N-1}$ before the product. Since this factor quickly becomes small as $N$ increases, removing it should improve the acceptance probability.

Both conditions (27) and (31) aim to put upper bounds on the volume density (22). This density contains $d-1$ factors ($d = N(N+1)/2$ being the dimension of $M$). The first $N-1$ factors are all equal to $r$ and the remaining $N(N-1)/2$ factors are $\sinh(\kappa_{ij}(s)r)/\kappa_{ij}(s)$ where the $\kappa_{ij}(s)$ are given by (23) for $1 \le i < j \le N$.

Condition (27) is obtained by upper bounding each one of these $d-1$ factors by the same quantity $\sinh(\kappa r)/\kappa$. This is indeed a correct upper bound because $\kappa_{ij}(s) \le \kappa$, as explained in Appendix B after Formula (49). However, this upper bound is an exaggerated overestimate for the first $N-1$ factors, because $\sinh(\kappa r)/\kappa$ increases much faster than $r$.

Condition (31) upper bounds the first $N-1$ factors (all equal to $r$) by $r$, and maintains the upper bound $\sinh(\kappa r)/\kappa$ for the remaining $N(N-1)/2$ factors. A sharper inequality is thus obtained, whose effect is to remove $(\kappa r/\sinh(\kappa r))^{N-1}$ before the product from (27).

As before, it is helpful to have a meaningful geometric interpretation of the method at hand. In Section 4, it was explained that the idea of the CURS method is to reject directions which display “atypical curvatures”. Now, the sectional curvatures encountered along a direction $s$ are determined by the principal curvatures (given by (49), Appendix B). Out of these $d-1$ principal curvatures, there are always $N-1$ equal to 0 (the others are negative). Since these $N-1$ principal curvatures equal to 0 occur along every direction $s$, they can hardly be considered “atypical”, even though they are different from the other principal curvatures which mostly have strictly negative values.

The main difference between general CURS (using Condition (27)) and sharp CURS (using Condition (31)), is that the former systematically discriminates against curvatures closer to positive values, whereas the latter (using specific information about the space of real covariance matrices) takes into account the existence of principal curvatures equal to 0.

Tables 1 and 2 show the improvement in acceptance probability obtained by using sharp CURS instead of CURS. These tables report empirical acceptance probabilities, each one estimated from $10^6$ iterations of either method. Note that standard deviations are not reported, since such a large number of iterations makes them negligible.

<table><caption>Table 1: ($N = 4$) Empirical acceptance probabilities for CURS ($\hat{\Pi}$) and sharp CURS ($\hat{\Pi}_s$)</caption><tr><td>σ</td><td>0.2</td><td>0.4</td><td>0.6</td><td>0.8</td><td>1.0</td><td>1.2</td><td>1.4</td></tr><tr><td>Π̂</td><td>0.7817</td><td>0.3430</td><td>0.0638</td><td>0.0031</td><td>0.0000</td><td>0</td><td>0</td></tr><tr><td>Π̂<sub>s</sub></td><td>0.8682</td><td>0.5510</td><td>0.2364</td><td>0.0606</td><td>0.0086</td><td>0.0006</td><td>0.0000</td></tr></table>
When $N = 4$ and $\sigma = 1$, Table 1 shows that the acceptance probability for CURS is less than $10^{-4}$, while sharp CURS has acceptance probability $86 * 10^{-4}$ — a rather large improvement. Thus, sharp CURS produces a new sample every $\approx 116$ iterations on average. Each iteration requires a time of the order of $10^{-2}$ seconds. For $N = 4$ and $\sigma = 1$ the variance $\delta$ in (30) is a pretty big 13.3 (compare to discussion in the previous subsection). This can be found from $\delta = \sigma^3 Z'(\sigma)/Z(\sigma)$ with $Z(\sigma)$ given by (28) (Said et al., 2023) (the prime denotes the derivative).

Table 2 reports a very similar situation when $N = 6$. For example, at $\sigma = 0.5$, the variance is $\delta = 7.47$. The acceptance probability for CURS is less than $10^{-4}$, while it is $11 * 10^{-4}$ for sharp CURS. The time required for a single iteration (of either method) is still of the order of $10^{-2}$ seconds.

# 6 Dealing with the cut locus

The present section discusses dropping the simplifying assumption made in Subsection 2.1. This was the assumption that geodesic spherical coordinates regularly cover the entire manifold $M$, so that the domain of injectivity $D_{x_o}$ satisfies $D_{x_o} = M$.

In general, the manifold $M$ is the disjoint union of $D_{x_o}$, which is an open subset of $M$, and of the so-called cut locus $\text{Cut}_{x_o}$. This cut locus is a closed subset of $M$ and has zero Riemannian volume (Chavel, 2006). Now, the target distribution (4) clearly assigns zero probability to any subset of $M$ that has zero volume. In particular, the probability that a sample $x$ should belong to $\text{Cut}_{x_o}$ is ( according to (4))

$$P(\text{Cut}_{x_o}) = \int_{\text{Cut}_{x_o}} P(dx) = \int_{\text{Cut}_{x_o}} p(x) \text{vol}(dx) = 0$$

where the last equality follows from the fact that $\text{Cut}_{x_o}$ has zero volume. For the sampling problem at hand, this has the following implication: it is possible to consider the target distribution (4) (with its density (1)) as a distribution on the injectivity domain $D_{x_o}$ and thereby to restrict sampling to this open subset of $M$.

This allows for the CURS algorithm to be extended to the fully general situation, where one has $D_{x_o} \neq M$, but at the cost of an additional modification. In fact, the final step of Algorithm 1 (step 8) produces a sample $x$ by replacing $(r, s)$ into the exponential map in (3). As of now, the couple $(r, s)$ should be restricted in such a way that (3) will only produce samples $x$ that belong to $D_{x_o}$. This restriction is indeed necessary. If it is not applied, the exponential map is not anymore a diffeomorphism (it may even fail to be bijective).

To introduce the correct restriction, one needs to know the so-called c-function (Chavel, 2006) (Page 115), $c: S_{x_o}M \to (0, \infty)$, where $c(s)$ is the least value of $r > 0$ for which $x(r, s)$ (given by (3)) does not belong to $D_{x_o}$. Using this function, the target joint density $p(r, s)$ from (9) is now replaced with the truncated version

$$p_c(r, s) = \frac{1}{Z} \times f(r) |\det A(r, s)| \times \mathbf{1}_{\{r<c(s)\}} \quad (32)$$

where $Z$ is a normalising constant, $\mathbf{1}_A$ denotes the indicator function of a set $A$ (equal to 1 on $A$ and to 0 on the complement of $A$), and the truncated density $p_c$ is with respect to the reference measure $dr\omega(ds)$.

In this new setting, the task of CURS is to sample $(r, s)$ from the target joint density (32) and then to obtain $x$ by replacing into (3). The general idea introduced in Subsection 2.2 can be employed once more. Namely, using a volume comparison inequality, obtain an upper bound on the target joint density (32), and then use this to perform rejection sampling.

<table><caption>Table 2: ($N = 6$) Empirical acceptance probabilities for CURS ($\hat{\Pi}$) and sharp CURS ($\hat{\Pi}_s$)</caption><tr><td>σ</td><td>0.1</td><td>0.2</td><td>0.3</td><td>0.4</td><td>0.5</td><td>0.6</td><td>0.7</td></tr><tr><td>Π̂</td><td>0.7377</td><td>0.2798</td><td>0.0449</td><td>0.0022</td><td>0.0000</td><td>0</td><td>0</td></tr><tr><td>Π̂<sub>s</sub></td><td>0.8067</td><td>0.4126</td><td>0.1224</td><td>0.0179</td><td>0.0011</td><td>0.0000</td><td>0</td></tr></table>
Here, this idea will not be formulated in general terms, but rather illustrated through a concrete example. For this example, $M$ is the group of $N \times N$ unitary matrices, $M = U(N)$, equipped with the trace metric

$$ \langle u, v \rangle_x = \text{tr}(uv^\dagger) \qquad (33) $$

where $\dagger$ denotes conjugation-transposition and the tangent vectors $u$ and $v$ at $x \in M$ are $N \times N$ matrices, such that $x^\dagger u$ and $x^\dagger v$ are skew-Hermitian.

Consider a target density of the general form (1), and assume without loss of generality that $x_o = \text{id}$, the $N \times N$ identity matrix. Then, to express the truncated density $p_c$ from (32), note the following facts (see Appendix C for background).

• the tangent space to $M$ at $x_o = \text{id}$ is the space of $N \times N$ skew-Hermitian matrices. The unit sphere $S_{x_o} M$ is then made up of those skew-Hermitian $s$ with $\|s\|^2 = 1$, where the norm is determined by the metric (33).

• for $r > 0$ and $s \in S_{x_o} M$, the volume density is given by

$$ \det A(r, s) = r^{N-1} \prod_{i<j} (\sin(\kappa_{ij}(s)r) / \kappa_{ij}(s))^2 \qquad (34) $$

where, similar to (23),

$$ \kappa_{ij}(s) = (\zeta_i - \zeta_j)/2 \quad 1 \le i < j \le N \qquad (35) $$

with the eigenvalues of $s$ being $j\zeta_i$ for $i = 1, \dots, N$ (here, $j = \sqrt{-1}$).

• the $c$-function is given in the following manner, for each $s \in S_{x_o} M$,

$$ c(s) = \frac{\pi}{\max_i |\zeta_i|} \qquad (36) $$

Formulas (34)–(36) completely determine the truncated density $p_c$ in (32), with the normalising constant

$$ Z = \int_{S_{x_o} M} \int_0^{c(s)} f(r) |\det A(r, s)| dr \omega(ds) \qquad (37) $$

As stated above, the job of CURS is to sample from this truncated density. In order to apply the idea of Subsection 2.2, one has to upper bound $p_c(r, s)$ using some other density $g(r, s)$ which does not depend on $s$. The chosen upper bound should reflect knowledge of the curvature of the manifold $M$ at hand. Specifically, this is a manifold with positive sectional curvatures (Chavel, 2006), so the volume comparison inequality (8) should be applied with $\kappa = 0$, which yields $|\det A(r, s)| \le r^{d-1}$ ($d = N^2$ is the dimension of $M$).

Replacing this last inequality into (32), and noting that the indicator function of $\{r < c(s)\}$ is always $\le 1$, one is immediately lead to

$$ p_c(r, s) \le T \times g(r, s) \qquad (38) $$

in terms of the following

$$ g(r, s) = \frac{1}{Z_0} \times f(r) r^{d-1} \times \mathbf{1}_{\{r < \Delta_{x_o}\}} \quad T = \frac{Z_0}{Z} \qquad (39) $$

where $Z_0$ is a normalising constant, and $\Delta_{x_o} = \max c(s)$ (where the maximum is taken over $s \in S_{x_o} M$). Note that $\Delta_{x_o}$ is the greatest possible Riemannian distance $d(x_o, x)$, between $x_o$ and any other point $x \in M$.

Having obtained (38) and (39), it is possible to spell out the steps of a CURS algorithm for sampling on the manifold $M = U(N)$. While specific to this manifold, these steps hopefully offer a clear idea of what a fully general CURS algorithm should look like.

**step 1:** as in Subsection 5.1, the first step is to sample $s$ from a uniform distribution on the sphere $S_{x_o} M$. To do so, let $s = jh$ where $h$ is distributed according to the Gaussian unitary ensemble, and then replace $s$ by $s/\|s\|$ (where the norm is determined by (33)). Note that $h = (t + t^\dagger)/2$ where $t$ has independent entries, $t_{ij} = a_{ij} + j b_{ij}$ with $a_{ij}$ and $b_{ij}$ independent, Gaussian, each with mean 0 and variance 1/2 (Forrester, 2010).
**step 2:** again as in Subsection 5.1, sample the distance $r$ from the proposal density $g(r) = g(r, s)$ in (39). However (unlike Subsection 5.1), then reject any couple $(r, s)$ such that $r \ge c(s)$, where $c(s)$ is given by (36). For (39), note that $\Delta_{x_o} = \sqrt{N\pi}$ (see Appendix C).

**step 3:** for the rejection procedure, sample $U$ from a uniform distribution in the unit interval (0, 1).

**step 4:** reject $(r, s)$ if

$$U > \prod_{i<j} (\sin(\kappa_{ij}(s)r)/\kappa_{ij}(s)r)^2 \qquad (40)$$

Once these four steps lead to acceptance of some couple $(r, s)$, it is enough to put $x = \exp(rs)$ where $\exp$ denotes the matrix exponential. This amounts to applying the general formula (3). Moreover, if $x_o \neq 0$, it is also necessary to re-centre the sample $x$ around $x_o$ using the transformation $x \mapsto x_o x$.

The above described steps closely mirror those in Subsection 5.1. The only difference appears in step 2, which has an additional rejection step that discards $(r, s)$ when $r \ge c(s)$. The manifold $M = U(N)$ is compact, and the Riemannian distance between any two points $x_o$ and $x$ cannot exceed $\sqrt{N\pi}$ (the diameter of $M$). If only for this reason, one cannot accept to sample a distance $r$ that is too large. Furthermore, without the restriction $r < c(s)$, the exponential map will fail to be a diffeomorphism, so that taking $x = \exp(rs)$ will not produce correct samples from the target distribution.

Steps 3 and 4 look just like the corresponding steps in Section 5.1 but the rejection condition (40) works in quite a different way from (27) (and also from (31)). Condition (40) favors directions $s$ which have eigenvalues $(\zeta_i)$ close to one another (these have small $\kappa_{ij}(s)$, by (35)), while the previous (27) tends to reject $s$ with close eigenvalues. This is a striking difference between how CURS works for spaces of positive curvature (the present case) and spaces of negative curvature (Section 5).

A complete analysis of the above CURS algorithm for sampling on $M = U(N)$ would require evaluation of the acceptance probability $\Pi = 1/T$ (where $T$ is given by (39)). This will not be carried out here, since the mathematical tools required are somewhat outside of the scope of the present work.

It is helpful to have a feeling for (or intuition about) the acceptance region $r < c(s)$. When $N = 2$, this can be obtained through visualisation in a two-dimensional plane, with the coordinates $x_1 = r\zeta_1$ and $x_2 = r\zeta_2$, as in Figure 3. In this figure, the interior square encloses the acceptance region $r < c(s)$, while the outer circle corresponds to $r = \Delta_{x_o}$ (this is equal to $\sqrt{2\pi}$, since $N = 2$). In step 2 (which was just described), the proposal density produces samples inside the circle, but only those that fall inside the square are retained.

The samples which fall inside the square are then submitted to the rejection condition (40), whose right-hand side is just $(\text{sinc}((x_1 - x_2)/2))^2$ (where $\text{sinc}(x) = \sin(x)/x$). This is equal to 1 when $x_1 = x_2$, so samples which fall along the main diagonal have a probability equal to 1 of being accepted. The probability of acceptance decreases along the antidiagonal, and becomes zero at the corners of the square, $(x_1, x_2) = (-\pi, \pi)$ or $(\pi, -\pi)$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Proposal and acceptance regions for CURS ($M = U(2)$)</figcaption></figure>
A Integration in spherical coordinates

The aim here is to explain how the integral

$$
\int_M f(x) \mathrm{vol}(dx) \qquad (41)
$$

where $f: M \to \mathbb{R}$ and vol denotes the Riemannian volume, can be expressed in geodesic spherical coordinates, introduced in Subsection 2.1. The following material is based on (Chavel, 2006) (Chapter III).

Let $M$ be a complete Riemannian manifold and $x_o \in M$. Then, $M$ is the disjoint union of the open subset $D_{x_o}$, called the domain of injectivity, and of the closed subset $\text{Cut}_{x_o}$, called the cut locus.

Now, the cut locus has zero Riemannian volume, so the integral (41) may be restricted to its complement

$$
\int_M f(x) \mathrm{vol}(dx) = \int_{D_{x_o}} f(x) \mathrm{vol}(dx) \quad (42)
$$

In addition, the mapping $x \mapsto x(r, s)$ from (3) becomes a diffeomorphism onto $D_{x_o}$ when the couple $(r, s)$ is restricted in a certain way. This restriction is given by the so-called $c$-function (Chavel, 2006) (Page 115). This is $c: S_{x_o}M \to (0, \infty)$, where $c(s)$ is the least value of $r > 0$ for which $x(r, s)$ does not belong to $D_{x_o}$.

Then, integration over $D_{x_o}$ can be replaced by integration over $(r, s)$ such that $r < c(s)$. From (42),

$$
\begin{align}
\int_M f(x) \mathrm{vol}(dx) &= \int_{D_{x_o}} f(x) \mathrm{vol}(dx) \nonumber \\
&= \int_{S_{x_o}M} \int_0^{c(s)} f(r, s) \mathrm{vol}(dr \times ds) \tag{43}
\end{align}
$$

where $f(r, s) = f(x(r, s))$ and $\mathrm{vol}(dr \times ds)$ is just $\mathrm{vol}(dx)$, but expressed in terms of $r$ and $s$.

This was given by (6) in Subsection 2.1, which reads

$$
\mathrm{vol}(dr \times ds) = |\det A(r, s)| \, dr \omega(ds)
$$

where $\omega$ is the surface area measure on the unit sphere $S_{x_o}M$ and the matrix $A(r,s)$ is obtained by solving
a second-order linear differential equation, known as the matrix Jacobi equation (Chavel, 2006) (Page 114).

Replacing into (43) gives the general integral formula in geodesic spherical coordinates

$$
\int_M f(x) \mathrm{vol}(dx) = \int_{S_{x_o}M} \int_0^{c(s)} f(r,s) |\det A(r,s)| \, dr \omega(ds) \quad (44)
$$

Formulas (16) and (17) in Section 3 are both special cases of (44). Indeed, (16) is just the same as (44), but with $c(s) = \infty$. This reflects the assumption, made in Section 2, that $D_{x_o} = M$. On the other hand, (17) is a consequence of (16), which follows by putting

$$
\det A(r, s) = (\kappa^{-1} \sinh(\kappa r))^{d-1} \quad (45)
$$

which does not depend on $s$, and then integrating over $s$, which leads to the factor $\Omega_{d-1}$. In (16) and (17), by comparison with (44), $f(r,s) = f(r)$ does not depend on $s$.

The key to using (44) is the matrix Jacobi equation. Note that $A(r,s)$ is a $(d-1) \times (d-1)$ matrix ($d = \dim M$), and the matrix Jacobi equation reads (the prime denotes differentiation with respect to $r$)

$$
A''(r, s) + R(r, s) A(r, s) = 0 \qquad A(0, s) = 0 \text{ and } A'(0, s) = I_{d-1} \tag{46}
$$

where $R(r, s)$ is a matrix given in terms of the Riemann curvature tensor of $M$ (Chavel, 2006) (Page 114), and $I_{d-1}$ is the $(d-1) \times (d-1)$ identity matrix. For example, if $M$ has constant negative curvature $-\kappa^2$, then $R(r, s) = -\kappa^2 I_{d-1}$, and the solution of (46) is $A(r, s) = \kappa^{-1} \sinh(\kappa r) I_{d-1}$. This immediately gives (45).
## B Background for Section 5

The present appendix provides some additional background on the Riemannian geometry of the space M, of $N \times N$ real covariance matrices, equipped with the affine-invariant metric (20). Essentially, the aim is to apply the integral formula (44) to this space, and this will require solving the matrix Jacobi equation (46).

The main ingredient in this equation is the Riemann curvature tensor. Recall that the tangent space to M at $x_o = id$ is identified with the vector space of $N \times N$ real symmetric matrices. Denote this space by $T_{x_o}M$. The Riemann curvature tensor is a map $R: T_{x_o}M \times T_{x_o}M \times T_{x_o}M \to T_{x_o}M$, which is linear in each one of its three arguments (Lee, 2018). In the present case, it is given by (Helgason, 2001),

$$R(u, v)w = [w, [u, v]] \qquad (47)$$

where $[a, b] = ab - ba$ for any matrices $a, b$.

The space M at hand is moreover a symmetric space. In particular, its curvature tensor is parallel (Helgason, 2001), and this implies the matrix $R(r, s)$ in (46) is constant (does not depend on r). Moreover, it is just the matrix (in any chosen basis of $T_{x_o}M$) of the radial curvature operator (Lee, 2018),

$$R_s(u) = R(s, u)s \qquad (48)$$

This operator $R_s: T_{x_o}M \to T_{x_o}M$ is diagonalisable, with the following eigenvalues: 0 with multiplicity $N$, and $-(\kappa_{ij}(s))^2$ (given by (23)) for $1 \le i < j \le N$, with multiplicity 1 (see (Said, 2021) for details). In addition, the corresponding eigenvectors form an orthonormal basis of $T_{x_o}M$.

Let $(u_\alpha; \alpha = 1, \dots, d)$ be such an orthonormal basis, where each $u_\alpha$ corresponds to the eigenvalue $\lambda_\alpha$ ($d = N(N+1)/2$ is the dimension of M). In this basis, the matrix $A(r, s)$ (the solution of (46)) is diagonal, and its diagonal elements are solutions to ordinary differential equations,

$$A''_{\alpha\alpha}(r, s) = -\lambda_{\alpha} A_{\alpha\alpha}(r, s) \qquad A_{\alpha\alpha}(0, s) = 0 \text{ and } A'_{\alpha\alpha}(0, s) = 1$$

For the $N$ eigenvalues $\lambda_\alpha = 0$, this yields $A_{\alpha\alpha}(r, s) = r$. For the remaining eigenvalues $\lambda_\alpha = -(\kappa_{ij}(s))^2$, one has $A_{\alpha\alpha}(r, s) = \sinh(\kappa_{ij}(s)r)/\kappa_{ij}(s)$. Multiplying these together yields

$$\det A(r, s) = r^N \prod_{i<j} (\sinh(\kappa_{ij}(s)r)/\kappa_{ij}(s))$$

The volume density (22) is obtained by replacing $r^N$ with $r^{N-1}$, since one of the eigenvectors with $\lambda_\alpha = 0$ is $u_\alpha = s$, which is normal to the sphere $S_{x_o}M$ and should not be included when integrating over this sphere.

Assume that $u_1 = s$, so each $u_\alpha$ with $\alpha > 1$ is orthogonal to $s$. Then, the vectors $s$ and $u_\alpha$ span a two-dimensional plane with sectional curvature (Said, 2021)

$$\sec(s, u_\alpha) = \lambda_\alpha = 0 \text{ or } -(\kappa_{ij}(s))^2 \qquad (49)$$

These are known as principal sectional curvatures, and the sectional curvatures of M at any point are included between the minimum and maximum of these $\lambda_\alpha$. It is clear from (23) that those $s$ which have eigenvalues $(\zeta_i)$ closer to one another encounter sectional curvatures closer to 0.

The sectional curvatures of M are always negative or zero, by (49). The least possible sectional curvature is

$$-\kappa^2 = \min - (\kappa_{ij}(s))^2$$

where the minimum is over $s \in S_{x_o}M$. A straightforward differentiation shows that this minimum is achieved when there are two eigenvalues $\zeta_i = -\zeta_j$ both with absolute value $|\zeta_i| = |\zeta_j| = 1/\sqrt{2}$. This implies that $-\kappa^2 = -(\kappa_{ij}(s))^2 = -1/2$, and thus justifies the value of $\kappa = 1/\sqrt{2}$ used in Section 5.

Note that the presentation in this appendix is completely specific to the space of real covariance matrices. The following appendix puts it in a more general context, by considering the application of CURS to Riemannian symmetric space of non-compact type.
## C Background for Section 6

The present appendix provides a background discussion for the facts used in Section 6 in order to express the truncated density (32) when $M = U(N)$.

The first one of these three facts states that the tangent space to $M = U(N)$ at $x_o = \text{id}$ is the space of $N \times N$ skew-Hermitian matrices. This is a well known statement in the elementary theory of Lie groups, and a short proof may be found (for example) in (Hall, 2015) (Chapter 3).

The second fact is the expression (34) for the volume density. This can be derived following a reasoning analogous to the one sketched in Appendix B, for the space of real covariance matrices. Indeed, $M = U(N)$ is also a Riemannian symmetric space (because it is a compact Lie group (Helgason, 2001)). Its curvature tensor is therefore parallel, and the matrix $R(r, s)$ in the matrix Jacobi equation (46) is constant.

The eigenvalues of this matrix, which are just the same as those of the radial curvature operator (48), are the following: 0 with multiplicity $N$, and $(\kappa_{ij}(s))^2$ (given by (35)) for $1 \le i < j \le N$, each with multiplicity 2 (details can be found in (Said, 2021)).

Knowledge of these eigenvalues leads to the solution of the matrix Jacobi equation (46), as in Appendix B. For the determinant $\det A(r, s)$, each zero eigenvalue contributes a factor $r$, while each eigenvalue $(\kappa_{ij}(s))^2$ contributes a factor $\sin(\kappa_{ij}(s)r)/\kappa_{ij}(s)$ — note that since these eigenvalues are positive, rather than negative as in Appendix B, there is a sine function instead of the hyperbolic sine function.

Upon multiplying these factors together,

$$ \det A(r, s) = r^N \prod_{i<j} \left( \frac{\sin(\kappa_{ij}(s)r)}{\kappa_{ij}(s)} \right)^2 $$

Then, (34) is recovered by replacing $r^N$ with $r^{N-1}$. This is for the same reason as in Appendix B. Namely, one of the eigenvectors corresponding to an eigenvalue equal to zero happens to be $s$, which should not counted in the volume density.

Finally, the third fact is the expression (36) for the c-function. This can be found in (Sakai, 1977), which computes $c(s)$ using the characterisation that $c(s)$ is the least value of $r > 0$ for which there exists $s' \in S_{x_o} M$, different from $s$, such that $x(r, s) = x(r, s')$ (in the notation of (3)). In other words, $c(s)$ marks the first time the geodesic curve starting at $x_o$ with initial velocity $s$ intersects another geodesic curve also starting at $x_o$.

The sectional curvatures of $M = U(N)$ belong to the interval between the minimum and maximum of the eigenvalues 0 and $(\kappa_{ij}(s))^2$ of the radial curvature operator. Clearly, then, this manifold $M$ has non-negative curvature, and the maximum sectional curvature is $\max(\kappa_{ij}(s))^2$ (with the maximum taken over $s \in S_{x_o} M$). It is possible to show this maximum is equal to $1/2$.

The diameter of $M = U(N)$ is the greatest possible Riemannian distance between two points $x_o$ and $x$ in $M$. This is $\Delta_{x_o} = \max c(s)$ with $c(s)$ given by (36). This maximum occurs when all $\zeta_i$ are equal, and therefore equal to $1/\sqrt{N}$. This yields the diameter $\sqrt{N\pi}$. In fact, if $x_o = \text{id}$ then the maximum Riemannian distance $d(x_o, x)$ occurs for $x = -\text{id}$, and only for this point.

## D CURS for symmetric spaces

The development in Section 5 and Appendix B directly generalises to a large class of Riemannian manifolds. This is the class of negatively-curved Riemannian symmetric spaces, of which the space of real covariance matrices is just one instance, but which also includes spaces of complex, quaternion, Toeplitz or block Toeplitz covariance matrices (Said et al., 2018).

The present appendix presents a general version of CURS, for negatively-curved Riemannian symmetric spaces. For now, the general treatment for positively-curved spaces (generalising Section 6 and Appendix C) will not be detailed. This is just to avoid making the present paper too long, and to focus on general ideas rather than technical details.
Now, let $M$ be a negatively-curved Riemannian symmetric space (Helgason, 2001). In particular, $M$ is a Hadamard manifold and therefore satisfies the simplifying assumption of Subsection 2.1 (that $D_{x_o} = M$). Therefore, application of CURS on this space reduces to the steps described in Algorithm 1 of Section 3.

Application of these steps requires knowledge of three things: (a) the lower bound $-\kappa^2$ on the sectional curvatures of $M$, (b) the volume density det $A(r, s)$ (e.g. as in (22)), (c) finally (although this is not part of the algorithm in itself) the normalising constants $Z$ and $Z_\kappa$ which are needed to evaluate the acceptance probability $\Pi$ in (15).

To describe these three things, the most important mathematical objects will be the so-called positive restricted roots of the Lie algebra of the group of isometries of $M$ (Helgason, 2001). These will first be presented in general form, and then discussed through basic examples.

To begin, recall that $M$ admits a transitive action of a group of isometries $G$. Each isometry $g \in G$ is a mapping $g: M \to M$ that preserves the Riemannian metric of $M$. The subgroup of $G$, made up of those $g \in G$ that leave invariant the point $x_o \in M$, in the sense that $g(x_o) = x_o$ (the notation is that of (1)), will be denoted $K$.

Without any loss of generality, it is possible to think of the group $G$ as a matrix Lie group (Hall, 2015). This means that $G$ is a closed subgroup of the group $GL(q, \mathbb{R})$ of invertible $q \times q$ matrices (for some $q \ge 2$). The Lie algebra $\mathfrak{g}$ of $G$ is then a certain vector space of $q \times q$ matrices (not necessarily invertible), which is closed under the bracket operation: $u, v \in \mathfrak{g}$ implies $[u, v] \in \mathfrak{g}$ (where $[u, v] = uv - vu$).

This Lie algebra $\mathfrak{g}$ admits a Cartan decomposition $\mathfrak{g} = \mathfrak{k} + \mathfrak{p}$. This is a direct sum of vector spaces where $\mathfrak{k}$ is the Lie algebra of $K$ and $\mathfrak{p}$ is a subspace of $\mathfrak{g}$, which can be identified with the tangent space to $M$ at $x_o$.

A subspace $\mathfrak{a}$ of $\mathfrak{p}$ is called abelian if $[u, v] = 0$ for any $u, v \in \mathfrak{a}$. If $\mathfrak{a}$ is such an abelian subspace, which has maximal dimension, then each $v \in \mathfrak{p}$ can be put under the form $v = k\mathfrak{a}k^{-1}$ for some $k \in K$ and $a \in \mathfrak{a}$ — note that this is just a matrix product. In fact, it should be thought of as a “diagonalisation” of $v$.

Formulas (47) and (48) from Appendix B remain true in general. The radial curvature operator is given by

$$R_s(u) = [s, [s, u]] \qquad s, u \in T_{x_o} M \simeq \mathfrak{p} \qquad (50)$$

where $s, u$ belong to $T_{x_o} M$ (the tangent space to $M$ at $x_o$), which is identified with the subspace $\mathfrak{p}$ of $\mathfrak{g}$. It is important to understand that (50) should only be applied after $s, u$ have been recast as elements of $\mathfrak{p}$.

With this in mind, if $s = k\zeta k^{-1}$, with $k \in K$ and $\zeta \in \mathfrak{a}$, the eigenvalues of the radial curvature operator $R_s$ are the following: 0 with multiplicity $N$ (known as the rank of $M$), and $-(\kappa(\zeta))^2$ with respective multiplicities $m_\kappa$, where each $\kappa: \mathfrak{a} \to \mathbb{R}$ is a linear function, knows as a positive restricted root of the Lie algebra $\mathfrak{g}$ with respect to the subspace $\mathfrak{a}$.

As in Appendix B, the eigenvalues of the radial curvature operator play the role of principal sectional curvatures. The sectional curvatures of $M$ are therefore included between the minimum and maximum of these eigenvalues. It follows that $M$ has negative sectional curvatures, and that these are bounded below by

$$-\kappa^2 = \min - (\kappa(\zeta))^2 \qquad (51)$$

where the minimum is taken over all positive restricted roots $\kappa$ and all $s \in S_{x_o} M$, where $s = k\zeta k^{-1}$.

The eigenvalues of the radial curvature operator also allow one to directly write down the volume density det $A(r, s)$. Each eigenvalue equal to 0 contributes a factor $r$, while each eigenvalue $-(\kappa(\zeta))^2$ contributes a factor $\sinh(\kappa(\zeta)r)/\kappa(\zeta)$. Just as in Appendix B (or also in Appendix C), this yields

$$\det A(r, s) = r^{N-1} \prod_{\kappa} (\sinh(\kappa(s)r)/\kappa(s))^{m_{\kappa}} \qquad (52)$$

where the product is over all positive restricted roots. Note here that $\kappa(\zeta)$ has been written as $\kappa(s)$, by analogy with (22). There is no ambiguity in using this notation, since one can prove that the product in (52) only depends on $s$ and not on the particular decomposition $s = k\zeta k^{-1}$.
Finally, for the acceptance probability $\Pi$ in (15), Formulas (16) and (52) should be used in computing $Z$, while Formula (17) remains the same for $Z_\kappa$. The issue of obtaining expressions similar to the one in (28) (in other words, determinantal or Pfaffian expressions for $Z$), in the general setting described by (16) and (52), is the subject of ongoing research(Said & Mostajeran, 2023).

To make the above general treatment of negatively-curved symmetric spaces more concrete, it is helpful to consider three basic examples. In these examples, $M = M_\beta$ is the space of $N \times N$ covariance matrices, whose entries are real numbers, complex numbers, or quaternions, according to $\beta = 1, 2$, or 4. The group $G$ of isometries is the group of transformations $g : M_\beta \to M_\beta$ of the form $g(x) = gxg^\dagger$, where $g$ is an $N \times N$ invertible (real, complex, or quaternion) matrix and $\dagger$ denotes the conjugate-transpose (Said, 2021).

Therefore, $G$ may be identified with the matrix Lie group $GL(N, \mathbb{K}_\beta)$ of invertible matrices $g$ with entries from $\mathbb{K}_\beta = \mathbb{R}, \mathbb{C}$ or $\mathbb{H}$. The Lie algebra $\mathfrak{g}$ of $G$ is typically denoted $gl(N, \mathbb{K}_\beta)$, and is just the vector space of $N \times N$ matrices $v$ (not necessarily invertible) with entries from $\mathbb{K}_\beta$.

The Cartan decomposition is $\mathfrak{g} = \mathfrak{t} + \mathfrak{p}$, where $\mathfrak{t}$ is the subspace of matrices $v \in \mathfrak{g}$ such that $v + v^\dagger = 0$, and $\mathfrak{p}$ is the subspace of matrices $v \in \mathfrak{g}$ such that $v - v^\dagger = 0$. Therefore, matrices in $\mathfrak{t}$ are skew-Hermitian and those in $\mathfrak{p}$ are Hermitian.

The subspace $\mathfrak{t}$ of $\mathfrak{g}$ is the Lie algebra of the subgroup $K$ of $G$, which is made up of $g \in G$ such that $gg^\dagger = \text{id}$ (Said, 2021). These are exactly those $g \in G$ for which $g(x_o) = x_o$ when $x_o$ is taken to be $\text{id}$ (the $N \times N$ identity matrix). On the other hand, the subspace $\mathfrak{p}$ of $\mathfrak{g}$ is isomorphic to the tangent space to $M_\beta$ at $x_o = \text{id}$. It is usual to identify $u$ in this tangent space with $\tilde{u} = \frac{1}{2}u$ in $\mathfrak{p}$.

A maximal abelian subspace $\mathfrak{a}$ of $\mathfrak{p}$ is (for example) the space of $N \times N$ real matrices (these are real even when $\mathbb{K}_\beta = \mathbb{C}$ or $\mathbb{H}$). This is clearly abelian because $[u, v] = 0$ whenever $u$ and $v$ are diagonal matrices. Then, any $v \in \mathfrak{p}$ can be written $v = k\mathfrak{a}k^{-1}$ where $k \in K$ and $\mathfrak{a} \in \mathfrak{a}$ — this is just a general statement of the fact that any Hermitian matrix is unitarily diagonalisable.

The unit sphere $S_{x_o} M_\beta$ is the set of Hermitian matrices $s$ with $\text{tr}(s^2) = 1$. The radial curvature operator $R_s$ is given by (50). Here, one should be careful to replace $s$ by $\tilde{s} = \frac{1}{2}s$ and $u$ by $\tilde{u} = \frac{1}{2}u$ into (50). Then,

$$R_s(u) = \frac{1}{4}[s, [s, u]] \qquad (53)$$

This has a factor $\frac{1}{4}$ instead of $\frac{1}{8}$ because after computing the right-hand side of (50), one has to multiply by 2 in order to return from $\mathfrak{p}$ to the tangent space at $x_o$ of $M_\beta$ (i.e. from $\tilde{v} = \frac{1}{2}v$ to $v$, where $v$ denotes $R_s(u)$).

If $s = k\zeta k^{-1}$ where $k \in K$ and $\zeta \in \mathfrak{a}$, so $\zeta$ is a real diagonal matrix, then the eigenvalues of $R_s$ are the following: 0 with multiplicity $N$, and $-(\kappa_{ij}(\zeta))^2$ for $1 \le i < j \le N$, each one with multiplicity $\beta$, where $\kappa_{ij}(\zeta) = (\zeta_i - \zeta_j)/2$ when $\zeta = \text{diag}(\zeta_i; i = 1, \dots, N)$. For the three examples at hand ($M = M_\beta$), these $\kappa_{ij}$ (reminiscent of (23)) are the positive restricted roots.

From (51), it is possible to find $\kappa = 1/\sqrt{2}$, while (52) gives the volume density

$$\det A(r, s) = r^{N-1} \prod_{i<j} (\sinh(\kappa_{ij}(s)r)/\kappa_{ij}(s))^\beta \qquad (54)$$

which reduces to (22) when $\beta = 1$.

The three examples just considered are rather immediate generalisations of the space of real covariance matrices from Section 5. However, an extensive variety of interesting matrix spaces can be understood using symmetric space techniques which were reviewed in the present appendix (see (Edelman & Jeong, 2023)). Using the information about positive restricted roots (for example in (Said et al., 2018)), Formula (52) can be written down for other spaces of covariance matrices, such as spaces of Toeplitz or block-Toeplitz covariance matrices, and this directly leads to the adequate implementation of CURS.
## References

K. Bharath, A. Lewis, A. Sharma, and Tretyakov M. V. Sampling and estimation on manifolds using the langevin diffusion. *Journal of Machine Learning Research*, 26:1–50, 2025.

I. Chavel. *Riemannian geometry: a modern introduction*. Cambridge University Press, 2006.

Y. Chikuse. *Statistics on special manifolds*. Springer Science, 2003.

L. Devroye. *Non-uniform random variate generation*. Springer-Verlag, 1986.

A. Edelman and S. Jeong. Fifty three matrix factorizations: A systematic approach. *SIAM Journal on Matrix Analysis and Applications*, 44(2):415–480, 2023.

P. J. Forrester. *Log-gases and random matrices*. Princeton University Press, 2010.

M. Girolami and B. Calderhead. Riemannian manifold Langevin and Hamiltonian Monte Carlo methods. *Journal of the Royal Statistical Society (series B)*, 73:123–214, 2011.

B. C. Hall. *Lie groups, Lie algebras, and representations*. Springer International Publishing, 2015.

S. Helgason. *Differential geometry, Lie groups, and symmetric spaces*. American Mathematical Society, 2001.

J. M. Lee. *Introduction to Riemannian geometry*. Springer International Publishing, 2018.

G. Livan, M. Novaes, and P. Vivo. *Introduction to random matrices: theory and practice*. Springer Cham, 2018.

E. S. Meckes. *The random matrix theory of the classical compact groups*. Cambridge University Press, 2019.

X. Pennec, P. Fillard, and N. Ayache. A Riemannian Framework for Tensor Computing. *International Journal of Computer Vision*, 66:41–66, 2006.

C. Robert and G. Casella. *Monte Carlo statistical methods*. Springer-Verlag, 2004.

S. Said. Statistical models and probabilistic methods on Riemannian manifolds. *Habilitation thesis (HDR), Université de Bordeaux*, 2021. URL (arxiv:2101.10855).

S. Said and C. Mostajeran. Determinantal expressions of certain integrals on symmetric spaces. In *Geometric Science of Information*. Springer Nature Switzerland, 2023.

S. Said, H. Hajri, L. Bombrun, and B. C. Vemuri. Gaussian distributions on Riemannian symmetric spaces: statistical learning with structured covariance matrices. *IEEE Transactions on Information Theory*, 64:752–772, 2018.

S. Said, S. Heuveline, and C. Mostajeran. Riemannian statistics meets random matrix theory: Toward learning from high-dimensional covariance matrices. *IEEE Transactions on Information Theory*, 69:472–481, 2023.

T. Sakai. On the cut loci of compact symmetric spaces. *Hokkaido mathematical journal*, 6:136–161, 1977.

L. Santilli and M. Tierz. Riemannian Gaussian distributions, random matrix ensembles, and diffusion kernels. *Nuclear Physics B.*, 973, 2021.