Salem Said

Gaussian distributions on Riemannian symmetric spaces of non-positive curvature

Keywords: Gaussian distribution, Symmetric space, Random matrix theory, Markov-Chain Monte Carlo, Convex optimisation, Bayesian inference

This article aims to give a coherent presentation of the theory of Gaussian distributions on Riemannian symmetric spaces and also to report on recent original developments of this theory. The initial goal is to define a family of probability distributions, on any suitable Riemannian manifold, for which maximum-likelihood estimation, based on a finite sequence of observations, is equivalent to computation of the Riemannian barycentre of these observations. As it turns out, this goal is achievable whenever the underlying Riemannian manifold is a Riemannian symmetric space of non-positive curvature. In this case, the required Gaussian distributions are exactly the maximum-entropy distributions, for fixed barycentre and dispersion. The second step is the search for efficient means of computing the normalising factors associated with these distributions. This leads to a fascinating connection with random matrix theory, and even with theoretical physics (Chern-Simons theory), which yields a series of original results that provide exact expressions, as well as high-dimensional asymptotic expansions of the normalising factors. Another outcome of this connection with random matrix theory is the new idea of duality, between Gaussian distributions on Riemannian symetric spaces of opposite curvatures. The present article also investigates Bayesian inference for Gaussian distributions on symmetric spaces. This investigation motivates original results regarding Markov-Chain Monte Carlo and convex optimisation on Riemannian manifolds. It also reveals a new open problem (roughly, this concerns the equality of a posteriori mode with a posteriori barycentre), which should be the focus of future developments.

Introduction

The realisation that an essentially new approach, beyond that of classical statistics, is needed in order to learn from data that live in non-Euclidean spaces, can be credited to Fréchet, who invented what we today call the Fréchet mean, back in 1948 [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF].

The Fréchet mean generalises the concept of the mean (average or expectation) of a sequence of observations (x 1 , . . . , x N ), from the classical case where these observations lie in a Euclidean space, to the general setting where they belong to a non-Euclidean space.

Let us call our sample space M (the observations belong to M ). If M is a Euclidean space, it has a vector space structure, and the mean of (x 1 , . . . , x N ) is just the arithmetic mean (x 1 + . . . + x N )/N . If M is a non-Euclidean space, it will have no vector space structure, and this definition will lose all meaning.

To salvage the concept of mean, Fréchet suggested looking at the set of global minima of the sum of squared distances (the factor 1/2 is included for later convenience)

E(x) = 1 2 N n=1 d 2 (x n , x) for x ∈ M
He noted that any global minimum of E deserves to be called a mean of (x 1 , . . . , x N ). In this way, the mean of a sequence of observations in a non-Euclidean space is well-defined, at the cost of eventually failing to be unique. Fast-forward to the present, learning from data that live in Riemannian manifolds (a particular class of non-Euclidean spaces) has become central to many applications, ranging from radar signal processing to neuroscience [START_REF] Cabanes | Multidimensional complex stationary centered gaussian regressive time series classification: Application for audio and dar clutter machine learning in hyperbolic and siegel spaces[END_REF][START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF], and the Fréchet mean (more descriptively called the Riemannian barycentre) a very popular tool in this respect [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness and convexity[END_REF][START_REF] Said | Riemannian barycentres of Gibbs distributions : new results on concentration and convexity[END_REF].

Naturally, it also became important to provide a statistical (specifically, inferential) foundation for the use of this tool. One way or another, this lead to the quest for a suitable definition of a Gaussian distribution on a Riemannian manifold. This appeared inevitable, already because of the intimate connection between arithmetic means and Gaussian distributions, in the classical Eucidean case (see Paragraph 2.1, for discussion).

Gaussian distributons, defined as maximum entropy distributions on a Riemannian manifold, for a given Fréchet mean and dispersion, where first introduced by Pennec [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF]. For a while, it remained difficult to study these distributions, as there was no practical means of computing the associated normalising factors. However, a first breakthrough came when these factors where expressed as multiple integrals, in the case of Gaussian distributions on the space of real positive-definite matrices [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF].

In [START_REF] Bombrun | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF][9], the approach of [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF] was generalised to Gaussian distributions on Riemannian symmetric spaces of non-positive curvature, which include hyperbolic spaces, as well as spaces of real, complex and quaternion positive-definite matrices, and spaces of structured (Toeplitz or block-Toeplitz) positive-definite matrices. This opened the way to rigorous learning algorithms for data that live in these spaces (this is partially discussed in [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF]).

The introduction of Riemannian symmetric spaces reduced normalising factors of Gaussian distributions to multiple integrals, which could be computed using Monte Carlo techniques [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF]. Only very recently, it was realised that the techniques of random matrix theory made it possible to write down both analytic expressions and high-dimensional asymptotic expansions of these multiple integrals. This was studied by the theoretical physics community [START_REF] Santilli | Riemannian gaussian distributions, random matrix ensembles and diffusion kernels[END_REF] (see also our paper, currently under review [START_REF] Heuveline | Gaussian distributions on Riemannian symmetric spaces, random matrices, and planar feynman diagrams[END_REF]).

The aim of the present article is to give a coherent presentation of the theory of Gaussian distributions on Riemannian symmetric spaces of non-positive curvature, and report on recent original developments of this theory, including (but not limited to) the ones just mentioned. Its main body (Sections 2 and 3) relies on a variety of new results, mostly contained in the habilitation thesis [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF] -one advantage of this situation is that the flow of results is not interrupted by their sometimes lengthy proofs, given in [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF].

In the following, Section 2 introduced Gaussian distributions and their connection with random matrix theory. Section 3 investigates Bayesian inference of these distributions. Each one of these sections opens with a description of the original results which it contains.

Another, more modest, contribution of this article is Appendix B (not based on [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF]). This appendix provides new results on the convergence rates for Riemannian gradient descent, applied to strictly convex and strongly convex functions, defined on a convex subset of a Riemannian manifold. The main results are Propositions B.8 and B.9.

Gaussian distributions and RMT

The starting point of this section is a historical discussion of the concept of a Gaussian distribution. This leads up to the definition of Gaussian distributions, adopted in Paragraph 2.2, as a family of distributions P (x, σ) on a Riemannian manifold M , with parameters x ∈ M and σ > 0, for which maximum-likelihood estimation of x is equivalent to computation of the Riemannian barycentre. It turns out that this definition can be pursued whenever M is a Riemannian symmetric space of non-positive curvature.

Paragraph 2.3 then gives a general expression of the normalising factor Z(σ) of the Gaussian distribution P (x, σ), in the form of a multiple integral [START_REF] Bombrun | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF]. When M is a space of positive-definite matrices, or when M is the so-called Siegel domain, [START_REF] Bombrun | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF] is further reduced to a kind of integral familiar in random matrix theory ( [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF] and [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF], respectively).

Paragraph 2.4 states the existence and uniqueness of maximum-likelihood estimates of the parameters x and σ. It also states the maximum-entropy property of the Gaussian distribution P (x, σ), in Proposition 2.5. Paragraph 2.5 provides expressions of the barycentre (shown to be equal to x) and the covariance tensor of P (x, σ).

Paragraph 2.6 begins the series of results based on random matrix theory (RMT). These concern Gaussian distributions on the space H(N ) of complex positive-definite matrices. First, the analytic expression of Z(σ) is given in Proposition 2.8. Then, an asymptotic expansion of this expression, in the limit where N goes to infinity while t = N σ 2 remains constant, is given in Proposition 2.9.

Paragraph 2.7 describes the asymptotic distribution of eigenvalues of a random positivedefinite matrix in H(N ), drawn from the Gaussian distribution P (I N , σ) (I N denotes the N × N identity matrix). This asymptotic distribution has a probability density function, whose explicit expression is provided in Proposition 2.10.

Paragraph 2.8 introduces Θ distributions. These are classical normal distributions, wrapped around the unitary group U (N ), which is the dual symmetric space of H(N ). Proposition 2.11 uncovers an unexpected relationship between Θ distributions on U (N ) and Gaussian distributions on H(N ) : the normalising factors of these distributions are connected by a simple identity [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF].

Proofs of the above-mentioned results may be found in Chapter 3 of [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF].

From Gauss to Shannon

The story of Gaussian distributions is a story of discovery and re-discovery. Different scientists, at different times, were repeatedly lead to these distributions, through different routes. It seems the story began in 1801, on New Year's day, when Giuseppe Piazzi sighted a heavenly body (in fact, the asteroid Ceres), which he thought to be a new planet. Less than six weeks later, this "new planet" disappeared behind the sun. Using a method of least squares, Gauss predicted the area in the sky, where it re-appeared one year later. His justification of this method of least squares (cast in modern language) is that measurement errors follow a family of distributions, which satisfies Property 1 : maximum-likelihood estimation is equivalent to the least-squares problem.

In his Theoria motus corporum coelestium (1809), he used this property to show that the distribution of measurement errors is (again, in modern language) a Gaussian distribution.

In 1810, Laplace studied the distribution of a quantity, which is the aggregate of a great number of elementary observations. He was lead in this (completely different) way, to the same distribution discovered by Gauss. Laplace was among the first scientists to show Property 2 : the distribution of the sum of a large number of elementary observations is (asymptotically) a Gaussian distribution.

Around 1860, Maxwell rediscovered Gaussian distributions, through his investigation of the velocity distribution of particles in an ideal gas (which he viewed as freely-colliding perfect elastic spheres). Essentially, he showed that Property 3 : the distribution of a rotationally-invariant random vector, which has independent components, is a Gaussian distribution.

Kinetic theory lead to another fascinating development, related to Gaussian distributions. Around 1905, Einstein (and, independently, Smoluchowsky) showed that Property 4 : the distribution of the position of a particle, which is undergoing a Brownian motion, is a Gaussian distribution.

In addition to kinetic theory, alternative routes to Gaussian distributions have been found in quantum mechanics, information theory, and other fields. In quantum mechanics, a Gaussian distribution is a position distribution with minimum uncertainty. That is, it achieves equality in Heisenberg's inequality. In information theory, one may attribute to Shannon the following maximum-entropy characterisation Property 5 : a probability distribution with maximum entropy, among all distributions with a given mean and variance, is a Gaussian distribution.

The above list of re-discoveries of Gaussian distributions may be extended much longer. However, the main point is the following. In a Euclidean space, identified with R d , any one of the above five properties leads to the same famous expression of a Gaussian distribution,

P (dx| x, σ) = 2πσ 2 -d 2 exp - (x -x) 2 2σ 2 dx
as a probability distribution on R d , with mean vector x ∈ R d and variance parameter σ > 0 (here dx denotes the Lebesgue measure on R d ).

In non-Euclidean space, each one of these properties may lead to a different distribution, which may then be called a Gaussian distribution, but only from a restricted point of view. People interested in Brownian motion may call the heat kernel of a Riemannian manifold a Gaussian distribution on that manifold. However, statisticians will not like this definition, since it will (in general) fail to have a straightforward connection to maximum-likelihood estimation.

The "right" Gaussian

As of now, the following definition of Gaussian distributions is chosen. Gaussian distributions, on a Riemannian manifold M , are a family of distributions P (x, σ), parameterised by x ∈ M and σ > 0, such that : a maximum-likelihood estimate xN of x, based on samples (x n ; n = 1, . . . , N ) from P (x, σ), is a solution of the least-squares problem

minimise over x ∈ M E N (x) = 1 2 N n=1 d 2 (x n , x) (1)
This means that xN is an empirical barycentre of the samples (x n ). In order to construct probability distributions P (x, σ), which satisfy this definition, consider the density profile

f (x|x, σ) = exp - d 2 (x, x) 2σ 2 (2)
and the normalising factor,

Z(x, σ) = M f (x|x, σ) vol(dx) (3)
where vol denotes Riemannian volume. If this is finite, then

P (dx| x, σ) = (Z(x, σ)) -1 f (x|x, σ)vol(dx) (4)
is a well-defined probability distribution on M . In 2.4, below, it will be seen that P (x, σ), as defined by (4), is indeed a Gaussian distribution, if M is a Hadamard manifold and also a homogeneous space. The following propositions will then be helpful.

Proposition 2.1. Let M be a Hadamard manifold, whose sectional curvatures lie in [κ, 0], where κ = -c 2 . Then, for any x ∈ M and σ > 0, if Z(x, σ) is given by (3),

Z 0 (σ) ≤ Z(x, σ) ≤ Z c (σ) (5)
where

Z 0 (σ) = (2πσ 2 ) d 2 and Z c (σ) is positive, given by (d denotes the dimension of M ) Z c (σ) = ω d-1 σ (2c) d-1 d-1 k=0 (-1) k d -1 k Φ ((d -1 -2k)σc) Φ ((d -1 -2k)σc) ( 6
)
where ω d-1 denotes the area of the unit sphere in R d , and Φ denotes the standard normal distribution function.

Proposition 2.2. If M is a Riemannian homogeneous space, and Z(x, σ) is given by (3), then Z(x, σ) does not depend on x. In other words, Z(x, σ) = Z(σ).

If M is a Hadamard manifold and also a homogeneous space, then both Propositions 2.1 and 2.2 apply to M . Indeed, if M is a Riemannian homogeneous space, then its sectional curvatures lie within a bounded subset of the real line. Therefore, Proposition 2.1 implies Z(x, σ) is finite for all x ∈ M and σ > 0. On the other hand, Proposition 2.2 implies that Z(x, σ) = Z(σ). Then, (4) reduces to

P (dx| x, σ) = (Z(σ)) -1 exp - d 2 (x, x) 2σ 2 vol(dx) (7)
and yields a well-defined probability distribution P (x, σ) on M . This will be the main focus, throughout the following. Remark : here, readers may wish to recall the concept of a Hadamard manifold, or of a homogeneous space, from [START_REF] Petersen | Riemannian geometry[END_REF] (or any other good Riemannian geometry textbook).

The point of appealing to these concepts is the following. The assumption that M is a Hadamard manifold implies that geodesic spherical coordinates, which cover all of M , can be introduced at any point x ∈ M . Proposition 2.1 is obtained by writing the integral (3) in terms of these spherical coordinates, and then applying Riemannian volume comparison theorems, that state, very roughly speaking, that manifolds with more positive curvature have less volume. On the other hand, to say that M is a homogeneous space means that all points x ∈ M are equivalent, so changing the point of origin x does not change the integral (3). This is the key to Proposition 2.2.

The normalising factor Z(σ)

Assume now M = G/K is a Riemannian symmetric space which belongs to the noncompact case, described in Appendix A.1. In particular, M is a Hadamard manifold and also a homogeneous space. Thus, for each x ∈ M and σ > 0, there is a well-defined probability distribution P (x, σ) on M , given by [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF]. Here, the normalising factor Z(σ) can be expressed as a multiple integral, using the integral formula (74) of Proposition A.1, from Appendix A.1. Applying this proposition (with o = x), it is enough to note

f (ϕ(s, a)|x, σ) = exp - a 2 B 2σ 2
where a 2 B = B(a, a), in terms of the Ad(G)-invariant symmetric bilinear form B (see Appendix A.1). Since this expression only depends on a, (74) yields the following formula

Z(σ) = ω(S) |W | a exp - a 2 B 2σ 2 λ∈∆ + |sinh λ(a)| m λ da (8)
This formula expresses Z(σ) as a multiple integral on the vector space a. Recall that the dimension of a is known as the rank of M [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF]. Example 1 : the easiest instance of (8) arises when M is a hyperbolic space of dimension d, and constant sectional curvature equal to -1. Then, M has rank equal to 1, so that a = Râ for some unit vector â ∈ a. Since the sectional curvature is equal to -1, there is only one positive root λ, say λ(â) = 1, with multiplicity m λ = d -1. In addition, |W | = 2 because there are two Weyl chambers, C + = {râ ; r > 0} and C -= {râ ; r < 0}. Accordingly, [START_REF] Bombrun | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF] reads

Z(σ) = ω d-1 2 +∞ -∞ exp - r 2 2σ 2 |sinh(r)| d-1 dr = ω d-1 +∞ 0 exp - r 2 2σ 2 sinh d-1 (r)dr
In general, if all distances are divided by c > 0, the sectional curvature -1 is replaced by -c 2 . Thus, when M is a hyperbolic space of dimension d, and sectional curvature -c 2 ,

Z(σ) = ω d-1 +∞ 0 exp - r 2 2σ 2 (c -1 sinh(cr)) d-1 dr
This is exactly Z c (σ), expressed analytically in [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF]. Example 2 : another example, also susceptible of analytic expression, is when M is a space of positive-definite matrices with real, complex, or quaternion coefficients. Then, M = G/K with G = GL(N, K), where K = R, C or H (real numbers, complex numbers, or quaternions), and K ⊂ G a maximal compact subgroup, K = O(N ), U (N ) or Sp(N ).

In each of these three cases, a is the space of N × N real diagonal matrices, and the positive roots are the linear maps λ(a) = a ii -a jj where i < j, each one having its multiplicity m λ = β, (β = 1, 2 or 4, for K = R, C or H). In addition, a 2 B = 4tr(a 2 ). The Weyl group W is the groupe of permutation matrices in K, so |W | = N !, while S = K/T N where T N is the subgroup of all diagonal matrices in K. Replacing all of this into (8), it follows that

Z(σ) = ω β (N ) N ! a N i=1 exp - 2a 2 ii σ 2 i<j |sinh(a ii -a jj )| β da (9)
where ω β (N ) stands for ω(S), and da = da 11 . . . da N N . Introducing x i = exp(2a ii ),

Z(σ) = ω β (N ) 2 N N β N ! R N + |V (x)| β N i=1 ρ(x i , 2σ 2 )x -N β i dx i
where N β = (β/2)(N -1) + 1, ρ(x, k) = exp(-log 2 (x)/k) and V (x) = i<j (x j -x i ) is the Vandermonde determinant. Finally, using the elementary identity

ρ(x, k)x α = exp k 4 α 2 ρ e -k 2 α x, k it is immediately found that Z(σ) = ω β (N ) 2 N N β N ! × exp -N N 2 β (σ 2 /2) × R N + |V (u)| β N i=1 ρ(u i , 2σ 2 )du i (10)
For the case β = 2, the integral in (10) will be expressed analytically in 2.6, below.

Remark : curious readers will want to compute ω β (N ). For example, ω 2 (N ) can be found using the Weyl integral formula on U (N ) [START_REF] Knapp | Lie groups, beyond an introduction[END_REF]. This yields ω 2 (N ) = vol(U(N ))/(2π) N . The volume of the unitary group can be found by looking at the normalising factor of a Gaussian unitary ensemble [START_REF] Mehta | Random matrices[END_REF]. Specifically, vol(U(N )) = (2π) (N 2 +N )/2 /G(N ), in terms of G(N ) = Γ(1) × Γ(2) × . . . × Γ(N ) (Γ denotes the Euler Gamma function).

Example 3 : for this last example, let M = D N be the Siegel domain [START_REF] Siegel | Symplectic geometry[END_REF]. This is the set of N × N symmetric complex matrices z, such that I N -z † z is positive-definite. Here, M = G/K, where G Sp(N, R) (real symplectic group) and K U (N ) (unitary group). Precisely, G is the group of 2N × 2N complex matrices g, with g t Ωg = Ω and g † Γg = Γ, where t denotes the transpose, and where Ω and Γ are the matrices

Ω = I N -I N ; Γ = I N -I N
In addition, K is the group of block-diagonal matrices k = diag(U, U * ) where U ∈ U (N ), and * denotes the conjugate. The action of G on M is given by Möbius transformations,

g • z = (Az + B)(C z + D) -1 g = A B C D (11)
This action preserves the Siegel metric, which is defined by

v,v z = (I N -zz † ) -1 v 2 B v 2 B = 1 2 tr(vv † ) ( 12
)
where each tangent vector v is identified with a symmetric complex matrix. Now [START_REF] Terras | Harmonic analysis on symmetric spaces and applications[END_REF],

a = a a ; a = diag(a 11 , . . . , a N N ) (13)
The positive roots are λ(a) = a ii -a jj for i < j, and λ(a) = a ii + a jj for i ≤ j, all with m λ = 1. The order of the Weyl group is |W | = 2 N N !, and ω(S) = vol(U (N ))/2 N . Replacing into [START_REF] Bombrun | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF], it follows that

Z(σ) = vol(U (N )) 2 2N N ! a N i=1 exp - a 2 ii 2σ 2 i<j sinh |a ii -a jj | i≤j sinh |a ii + a jj | da (14)
or, after introducing u i = cosh(2a ii ),

Z(σ) = 2 -2N vol(U (N )) × C N V (u) N n=1 w(u i , 8σ 2 )du i ( 15
)
where

C N = {u ∈ R N + : u 1 ≤ u 2 ≤ u 1 ≤ . . . ≤ u N } and w(u, k) = exp(-acosh 2 (u)/k), while V (u)
is the Vandermonde determinant, as in [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF].

MLE and maximum entropy

Let M be a Hadamard manifold, which is also a homogeneous space. Consider the family of distributions P (x, σ) on M , given by [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF] for x ∈ M and σ > 0. This family of distributions fits the definition of Gaussian distributions, stated at the beginning of 2.2. Proposition 2.3. Let P (x, σ) be given by [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF], for x ∈ M and σ > 0. The maximumlikelihood estimate of the parameter x, based on samples (x n ; n = 1, . . . , N ) from P (x, σ), is unique and equal to the empirical barycentre xN of the samples (x n ). This proposition is almost immediate. From [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF], one has the log-likelihood function

(x, σ) = -N log Z(σ) - 1 2σ 2 N n=1 d 2 (x n , x) (16)
Since the first term does not depend on x, one may maximise (x, σ), first over x and then over σ. Clearly, maximising over x is equivalent to minimising the sum of squared distances d 2 (x n , x). This is just the least-squares problem (1), whose solution is the empirical barycentre xN . Moreover, xN is unique, since M is a Hadamard manifold [20][4]. Consider now maximum-likelihood estimation of σ. This is better carried out in terms of the natural parameter η = (-2σ 2 ) -1 , or in terms of the moment parameter δ = ψ (η), where ψ(η) = log Z(σ) and the prime denotes the derivative.

Proposition 2.4. The function ψ(η), just defined, is a strictly convex function, which maps the half-line (-∞, 0) onto R. The maximum-likelihood estimates of the parameters η and δ are ηN = (ψ ) -1 ( δN ) and δN = 1

N N n=1 d 2 (x n , xN ) (17)
where (ψ ) -1 denotes the reciprocal function.

Remark : ηN in ( 17) is well-defined, since the range of ψ is equal to (0, ∞). Indeed, one has the following inequalities, analogous to (5),

ψ 0 (η) ≤ ψ (η) ≤ ψ c (η) (18)
where ψ 0 (η) = log Z 0 (σ), and ψ c (η) = log Z c (σ). Now, ψ 0 (η) = nσ 2 , which increases to +∞ when σ increases to +∞. On the other hand, since η = (-2σ 2 ) -1 ,

ψ c (η) = σ 3 d dσ (log Z c (σ)) (19)
which, from [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF], is = 0 when σ = 0. Thus, it follows from (18) that ψ maps the half-line (-∞, 0) onto the half-line (0, +∞). An alternative definition of Gaussian distributions is provided by their maximum-entropy property, stated in the following proposition. Here, entropy specifically means Shannon's differential entropy. If P is a probability distribution on M , with probability density function p, this entropy is equal to

S(P ) = M (log p(x)) p(x)vol(dx)
Proposition 2.5. The Gaussian distribution P (x, σ) is the unique distribution on M , having maximum Shannon entropy, among all distributions P with given barycentre x and dispersion δ = E x∼P [d 2 (x, x)]. Its entropy is equal to ψ * (δ) where ψ * is the Legendre transform of ψ.

Barycentre and covariance

Let M be a Hadamard manifold, which is also a homogeneous space. Consider the barycentre and covariance of the Gaussian distribution P (x, σ) on M , given by [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF].

First, it should be noted P (x, σ) has a well-defined Riemannian barycentre, since it has finite second-order moments. To see that this is true, it is enough to note that [START_REF] Siegel | Symplectic geometry[END_REF]. Proposition 2.6. Let P (x, σ) be given by ( 7), for x ∈ M and σ > 0. The Riemannian barycentre of P (x, σ) is equal to x.

M d 2 (x, x)P (dx| x, σ) < ∞ Ineded, this integral is just ψ (η) in
The proof of this proposition relies on the fact that the so-called variance function

E(x) = 1 2 M d 2 (x, y)P (dy| x, σ) (20)
is strongly convex (see [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF], Paragraph 2.2.3). Thus, if gradE(x) = 0, then x is the global minimiser of E, and therefore the barycentre of P (x, σ). However, gradE(x) = 0 follows by a direct application of the following "Fisher's identity",

M (grad x log p(x|x, σ)) P (dx| x, σ) = 0
where grad x denotes the gradient with respect to x, defined according to the Riemannian metric of M , and p(x|x, σ) is the probability density function, appearing in [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF].

The covariance form of P (x, σ) is the symmetric bilinear form C x on T xM ,

C x (u, v) = M u, Exp -1 x (x) Exp -1 x (x), v p(x|x, σ)vol(dx) u , v ∈ T xM (21)
where Exp denotes the Riemannian exponential map (Exp -1 is well-defined, since M is a Hadamard manifold). With σ > 0 fixed, the map which assigns to x ∈ M the covariance form C x is a (0,2)-tensor field on M , here called the covariance tensor of P (x, σ). In order to compute this tensor field, consider the following situation.

Assume M = G/K is a Riemannian symmetric space. Here,

K = K o , the stabiliser in G of o ∈ M . For k ∈ K and u ∈ T o M , it is clear k•u ∈ T o M .
This defines a representation of K in the tangent space T o M , called the isotropy representation. One says that M is an irreducible symmetric space, if this isotropy representation is irreducible.

If M is not irreducible, then it is a product of irreducible Riemannian symmetric spaces M = M 1 ×. . .×M s [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF] (Proposition 5.5, Chapter VIII. This is the de Rham decomposition of M ). Accordingly, for x ∈ M and u ∈ T x M , one may write x = (x 1 , . . . , x s ) and u = (u 1 , . . . , u s ), where x r ∈ M r and u r ∈ T xr M r . Now, looking back at [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF], it may be seen that

p(x|x, σ) = s r=1 p(x r |x r , σ) p(x r |x r , σ) = (Z r (σ)) -1 exp - d 2 (x r , xr ) 2σ 2 (22)
For the following proposition, let η = (-2σ 2 ) -1 and ψ r (η) = log Z r (σ).

Proposition 2.7. Assume that M is a product of irreducible Riemannian symmetric spaces, M = M 1 × . . . × M s . The covariance tensor C in ( 21) is given by

C x(u, u) = s r=1 ψ r (η) dim M r u r 2 xr (23)
for u ∈ T xM where x = (x 1 , . . . , xs ) and u = (u 1 , . . . , u s ), with xr ∈ M r and u r ∈ T xr M r .

Example : let M = H(N ), the space of N × N Hermitian positive-definite matrices, so M = GL(N, C)/U (N ), with U (N ) the stabiliser of o = I N (N × N identity matrix). The de Rham decomposition of

M is M = M 1 × M 2 , where M 1 = R and M 2 is the submanifold whose elements are those x ∈ M such that det(x) = 1. Accordingly, each x ∈ M is identified with the couple (x 1 , x2 ), x1 = 1 N log det(x) x2 = (det(x)) -1/N x and each u ∈ T xM is written u = u 1 x + u 2 u 1 = 1 N tr(x -1 u) u 2 = u - 1 N tr(x -1 u) x
These may be replaced into expression [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF],

C x(u, u) = ψ 1 (η)u 2 1 + ψ 2 (η) N 2 -1 u 2 2 x2 (24)
where

ψ 1 (η) = log (2π σ 2 ) 1 2
, and ψ 2 (η) = log Z(σ) -ψ 1 (η) (Z(σ) is given by ( 26) in 2.6, below). After a direct calculation, this can be brought under the form

C x(u, u) = g 1 (σ)tr 2 (x -1 u) + g 2 (σ)tr(x -1 u) 2 (25)
where g 1 (σ) and g 2 (σ) are certain functions of σ.

Remark : as a corollary of Proposition 2.7, the covariance tensor C is a G-invariant Riemannian metric on M . This is clear, for example, in the special case of ( 25), which coincides with the general expression of a GL(N, C)-invariant metric.

Z(σ) from RMT

Random matrix theory is very helpful in the calculation of integrals such as ( 10) and ( 15), leading both to exact expressions and to asymptotic expansions of these integrals.

Here, this is illustrated for the integral [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF], with β = 2. This corresponds to M = H(N ), the space of N × N Hermitian positive-definite matrices. In this case, it is possible to provide an analytic formula for the normalising factor Z(σ).

Proposition 2.8. When M = H(N ), the normalising factor Z(σ), given by ( 10) with β = 2, admits of the following analytic expression

Z(σ) = ω 2 (N ) 2 N 2 2π σ 2 N 2 exp N 3 -N 6 σ 2 N -1 n=1 1 -e -nσ 2 N -n (26)
The proof of this proposition is a direct application of a well-known formula from random matrix theory [START_REF] Mehta | Random matrices[END_REF]. The integral in [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF] reads

I N (σ) = R N + |V (u)| β N i=1 ρ(u i , 2σ 2 )du i
According to [START_REF] Mehta | Random matrices[END_REF] (Chapter 5, Page 79), if (p n ; n = 0, 1, . . .) are orthonormal polynomials with respect to the weight function ρ(u, 2σ 2 ) on R + , then I N (σ) is given by

I N (σ) = N ! N -1 n=0 p -2 nn (27)
where p nn is the leading coefficient in p n . The required orthonormal polynomials p n are given by p n = (2πσ 2 ) -1 4 s n , where s n are the Stieltjes-Wigert polynomials [START_REF] Szegö | Orthogonal Polynomials[END_REF] (Page 33). Looking up the expression of these polynomials, it is easy to find

p -2 nn = 2π σ 2 1 2 exp (2n + 1) 2 2 σ 2 n m=1 1 -e -mσ 2
Then, working out the product ( 27) and replacing into [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF], one is lead to [START_REF] Udriste | Convex functions and optimization methods on Riemannian manifolds[END_REF]. Moving on, it is possible to derive an asymptotic expression of Z(σ), valid in the limit where N goes to infinity while the product t = N σ 2 remains constant. Proposition 2.9. Let Z(σ) be given by [START_REF] Udriste | Convex functions and optimization methods on Riemannian manifolds[END_REF]. If N → ∞, while t = N σ 2 remains constant, then the following equivalence holds,

1 N 2 log Z(σ) ∼ - 1 2 log 2N π + 3 4 + t 6 - Li 3 (e -t ) -ζ(3) t 2 (28)
where Li 3 (x) = ∞ k=1 x k /k 3 for |x| < 1 (the trilogarithm), and ζ is the Riemann Zeta function.

The main idea behind [START_REF] Meyn | Markov chains and stochastic stability[END_REF] is that, taking the logarithm in [START_REF] Udriste | Convex functions and optimization methods on Riemannian manifolds[END_REF], the product on the right-hand side turns into a Riemann sum for the improper integral

1 0 (1 -x) log 1 -e -tx dx = -(Li 3 (e -t ) -ζ(3))/t 2
where the equality follows by integrating term-by-term the power series of the logarithm.

The asymptotic distribution

From the point of view of random matrix theory, a Gaussian distribution P (I N , σ) on M = H(N ) defines a unitary matrix ensemble. If x is a random matrix, drawn from this ensemble, and (x i ; i = 1, . . . , N ) are its eigenvalues, which all belong to (0, ∞), then the empirical distribution ν N , which is given by (as usual,

δ x i is the Dirac distribution at x i ) ν N (B) = E 1 N N i=1 δ x i (B) (29)
for measurable B ⊂ (0, ∞), converges to an absolutely continuous distribution ν t , when N goes to infinity, while the product t = N σ 2 remains constant.

Proposition 2.10. Let c = e -t and a(t

) = c(1 + √ 1 -c) -2 while b(t) = c(1 - √ 1 -c) -2 .
When N goes to infinity, while the product t = N σ 2 remains constant, the empirical distribution ν N converges weakly to the distribution ν t with probability density function Remark : as one should expect, when t = 0 (so

dν t dx (x) = 1 πtx arctan 4e t x -(x + 1) 2 x + 1 1 [a(t),b(t)] (x) (30)
σ 2 = 0), a(t) = b(t) = 1.
The proof of Proposition 2.10 is a relatively direct application of a result in [START_REF] Kuijlaars | The asymptotic zero distribution of orthogonal polynomials with varying recurrence coefficients[END_REF] (Page 191). In fact, the intergration variables in [START_REF] Zanini | Parameters estimate of Riemannian gaussian distribution in the manifold of covariance matrices[END_REF] are u i = e t x i . Let νN be the empirical distribution of the u i (this is the same as [START_REF] Jarner | Geometric ergodicity of Metropolis algorithms[END_REF], but with u i instead of x i ). By applying [START_REF] Mehta | Random matrices[END_REF] (Chapter 5, Page 81),

νN (B) = 1 N B R (1) N (u)(du) (31)
for measurable B ⊂ (0, ∞), where the one-point correlation function R (1) N (u) is given by

R (1) N (u) = ρ(u, 2σ 2 ) N -1 n=0 p 2 n (u) (32)
in the notation of 2.6 (p n are orthonormal polynomials, with respect to the weight ρ(u, 2σ 2 )). According to [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF] (Page 133), νN given by ( 31) converges weakly to the so-called equilibrium distribution νt , which minimises the electrostatic energy functional

E(ν) = 1 t ∞ 0 1 2 log 2 (u)ν(du) - ∞ 0 ∞ 0 log |u -v|ν(du)ν(dv) (33)
over probability distributions ν on (0, ∞). Also according to [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF] (Page 133), this equilibrium distribution is the asymptotic distribution of the zeros of the polynomial p N (in the limit N → ∞ while N σ 2 = t). Fortunately, p N is just a constant multiple of the Stieltjes-Wigert polynomial s N [START_REF] Szegö | Orthogonal Polynomials[END_REF] (Page 33). Therefore, the required asymptotic distribution of zeros can be read from [START_REF] Kuijlaars | The asymptotic zero distribution of orthogonal polynomials with varying recurrence coefficients[END_REF] (Page 191). Finally, [START_REF] Chavel | Riemannian geometry, a modern introduction[END_REF] follows by introducing the change of variables x = e -t u. Remark : in [START_REF] Mariño | Chern-Simons theory, matrix models, and topological strings[END_REF], the equilibrium distribution νt is derived directly, by searching for stationary distributions of the energy functional [START_REF] Bogachev | Measure Theory[END_REF]. This leads to a singular integral equation, whose solution reduces to a Riemann-Hilbert problem. Astoundingly, the Gaussian distributions on H(N ), as introduced in the present chapter, provide a matrix model for Chern-Simons quantum field theory (a detailed account is given in [START_REF] Mariño | Chern-Simons theory, matrix models, and topological strings[END_REF]).

Duality : the Θ distributions

Recall the Riemannian symmetric space M = H(N ) of 2.6. Its dual space is the unitary group M * = U (N ) (the definition of duality may be found in Appendix A.2). Consider now a family of distributions on M * , which will be called Θ distributions, and which display an interesting connection with Gaussian distributions on M , studied in 2.6. Recall Jacobi's ϑ function 1 ,

ϑ(e iφ |σ 2 ) = +∞ m=-∞ exp(-m 2 σ 2 + 2miφ)
As a function of φ, up to some minor modifications, this is a wrapped normal distribution (in other words, the heat kernel of the unit circle),

1 2π ϑ(e iφ | σ 2 2 ) = ∞ m=-∞ exp - (2φ -2mπ) 2 2σ 2
Each x ∈ M * can be written x = k•e iθ where k ∈ U (N ) and e iθ = diag(e iθ j ; j = 1, . . . , N ).

Here, k • y = kyk † for y ∈ M * . With this notation, define the following matrix ϑ function,

Θ x σ 2 = k • ϑ(e iθ | σ 2 2 ) ( 34
)
which is obtained from x by applying Jacobi's ϑ function to each eigenvalue of x. Further, consider the positive function,

f * (x|x, σ) = det 2π σ 2 1 2 Θ xx † σ 2 (35)
where x ∈ M * . This is also equal to

det 2π σ 2 1 2 Θ x † x σ 2
since the matrices xx † and x † x are similar. Then, let Z M * (σ) denote the normalising constant

Z M * (σ) = M * f * (x|x, σ) vol(dx) (36)
which does not depend on x, as can be seen, by introducing the new variable of integration z = xx † , and using the invariance of vol(dx). Now, define a Θ distribution Θ(x, σ) as the probability distribution on M * , whose probability density function, with respect to vol(dx), is given by

p * (x|x, σ) = (Z M * (σ)) -1 f * (x|x, σ) (37)
Proposition 2.11. Let Z M (σ) = Z(σ) be given by [START_REF] Udriste | Convex functions and optimization methods on Riemannian manifolds[END_REF], and Z M * (σ) be given by [START_REF] Lee | Introduction to smooth manifolds[END_REF]. Then, the following equality holds

Z M (σ) Z M * (σ) = exp N 3 -N 6 σ 2 (38)
Remark : the Gaussian density (7) on M , and the Θ distribution density (37) on M * are apparently unrelated. Therefore, it is interesting to note their normalising constants Z M (σ) and Z M * (σ) scale together according to the simple relation [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. The connection between the two distributions is due to the duality between M and M * .

Gaussian distriutions and Bayesian inference

This section aims to investigate Bayesian inference for Gaussian distributions. Precisely, it aims to study Bayesian estimation of the parameter x of a Gaussian distribution P (x, σ), when this parameter is assigned a prior density which is also Gaussian, say P (z , τ ). Both the prior density P (z , τ ) and the likelihood density P (x, σ) are defined on a Riemannian symmetric space of non-positive curvature M while the minimum mean square error estimator xMMS , classically understood as the mean (expectation) of the posterior density, is here the Riemanian barycentre of π(x).

It is seen that xMAP can be computed directly, being a geodesic convex combination of the prior barycentre z and a new observation y, with respective weights 1 -ρ and ρ, where ρ = τ2 /(σ 2 + τ 2 ). On the other hand, xMMS seems much harder to compute.

However, Proposition 3.1 states that xMMS = xMAP if ρ = 1/2, and Proposition 3.2 states that, in the special case where M is a hyperbolic space, xMMS is a geodesic convex combination of z and y, just like xMAP , but with different weights, say (1 -t * ) and t * .

Paragraph 3.2 reports on numerical experiments which show, again in the special case where M is a hyperbolic space, that xMMS and xMAP lie very close to each other, and that they even appear to be equal (this would mean t * = ρ). At present, I am unaware of any mathematical explanation of this phenomenon.

Paragraph 3.3 describes the computational tools employed in calculating xMMS . First, Proposition 3.3 provides easy-to-verify sufficient conditions, for the geometric ergodicity of an isotropic Metropolis-Hastings Markov chain, in a Riemannian symmetric space M . These conditions are shown to apply in the case of the posterior density π(x), making it possible to generate geometrically ergodic samples (x n ; n ≥ 1) from this density.

The next Proposition 3.4 states that the empirical barycentre xN of the samples (x 1 , . . . , x N ) converges almost-surely to xMMS , so xN may be used to approximate xMMS to any required accuracy.

Concretely, computing the empirical barycentre xN requires solving a strongly convex optimisation problem on the Riemannian manifold M (here, convexity is with respect to the Riemannian connection of M [START_REF] Udriste | Convex functions and optimization methods on Riemannian manifolds[END_REF]. This has come to be called "geodesic convexity"). Appendix B is devoted to a brief but systematic study of convex optimisation on Riemannian manifolds. Specifically, it establishes the rate of convergence of Riemannian gradient descent schemes, applied to strictly convex or strongly convex cost functions.

The gradient descent schemes under consideration are retraction schemes (not limited to the Riemannian exponential) with a constant step-size. The problem is then to find the largest possible step-size which guarantees a certain rate of convergence. Proposition B.8 addresses this problem for strictly convex functions, and Proposition B.9 for strongly convex functions. For any strictly convex cost function, and suitable retraction, Proposition B.8 gives the largest possible step-size which guarantees a rate of convergence at least as fast as O(1/t) (t is the number of iterations). Proposition B.9 does the same for strongly convex functions, but with an exponential rate of convergence. In fact, with regard to the original motivation of computing xN , this ensures that a gradient descent scheme, using the Riemannian exponential, converges after only a few iterations.

MAP versus MMS

Assume M = G/K is a Riemannian symmetric space which belongs to the non-compact case, described in Appendix A.1. Recall the Gaussian distribution P (x, σ) on M given by its probability density function ( 7)

p(y|x, σ) = (Z(σ)) -1 exp - d 2 (y, x) 2σ 2 (39)
In 2.4, it was seen that maximum-likelihood estimation of the parameter x, based on samples (y n ; n = 1, . . . , N ), amounts to computing the empirical barycentre of these samples. The one-sample maximum-likelihood estimate, given a single observation y, is therefore xML = y. Instead of maximum-likelihood estimation, consider Bayesian estimation of x, based on the observation y. To do so, assign to x a prior density, which is also Gaussian,

p(x|z, τ ) = (Z(τ )) -1 exp - d 2 (x, z) 2τ 2 (40)
Upon observation of y, Bayesian inference concerning x is carried out using the posterior density

π(x) ∝ exp - d 2 (y, x) 2σ 2 - d 2 (x, z) 2τ 2 (41)
where ∝ indicates a missing (unknown) normalising factor.

In particular, the maximum a posteriori estimator xMAP of x is equal to the mode of the posterior density π(x). In other words, xMAP minimises the weighted sum of squared distances d 2 (y, x)/σ 2 + d 2 (x, z)/τ 2 . This is expressed in the following notation 3

, xMAP = z # ρ y where ρ = τ 2 σ 2 + τ 2 (42)
Thus, xMAP is a geodesic convex combination of the prior barycentre z and the observation y, with respective weights σ 2 /(σ 2 + τ 2 ) and τ 2 /(σ 2 + τ 2 ). On the other hand, the minimum mean square error estimator xMMS is the barycentre of the posterior density π(x). That is, xMMS is the unique global minimiser of

E π (w) = 1 2 M d 2 (w, x)π(x)vol(dx) (43)
3 If p, q ∈ M and c : [0, 1] → M is a geodesic curve with c(0) = p and c(1) = q, then p # t q = c(t), for t ∈ [0, 1]. In other words, p # t q is a geodesic convex combination of p and q, with respective weights (1 -t) and t.

While it is easy to compute xMAP from (42), it is much harder to find xMMS , as this requires minimising the integral (43), where the density π(x) is known only up to normalisation.

Still, there is one special case where these two estimators are equal.

Proposition 3.1. In the above notation, if

σ 2 = τ 2 (that is ρ = 1/2), then xMMS = xMAP .
When M is a Euclidean space, it is well-known that xMMS = xMAP for any value of ρ. When M is a space of constant negative curvature, the following proposition indicates xMMS and xMAP cannot be too far away from one another. Proposition 3.2. In the above notation, if M is a space of constant negative curvature (hyperbolic space), then xMMS = z # t * y for some t * ∈ (0, 1).

Bounding the distance

In general, one expects xMMS and xMAP to be different from one another, when ρ = 1/2. However, when M is a space of constant negative curvature, Proposition 3.2 shows the distace between these two estimators is always less than the distance between z and y.

Surprisingly (again when M is a space of constant negative curvature), numerical experiments show that xMMS and xMAP lie very close to each other, and that they even appear to be equal. I am unaware of any mathematical explanation of this phenomenon.

It is possible to bound the distance between xMMS and xMAP , using the fundamental contraction property [START_REF] Sturm | Probability measures on metric spaces of nonpositive curvature[END_REF] (this is an immediate application of Jensen's inequality, as explained in the proof of Theorem 6.3 in [START_REF] Sturm | Probability measures on metric spaces of nonpositive curvature[END_REF]).

d(x MMS , xMAP ) ≤ W (π, δ xMAP ) (44)
where W denotes the Kantorovich (L 1 -Wasserstein) distance, and δ xMAP denotes the Dirac probability distribution concentrated at xMAP . Now, the right-hand side of (44) is equal to the first-order moment

m 1 (x MAP ) = M d(x MAP , x)π(x)vol(dx) (45)
Of course, the upper bound in (44) is not tight, since it is strictly positive, even when ρ = 1/2, as one may see from (45).

It will be shown below that a Metropolis-Hastings algorithm, with Gaussian proposals, can be used to generate geometrically ergodic samples (x n ; n ≥ 1) from the posterior density π. It is therefore possible to approximate (45) by an empirical average

m1 (x MAP ) = 1 N N n=1 d(x MAP , x n ) (46)
In addition, the samples (x n ) can be used to compute a convergent approximation of xMMS . Precisely, the empirical barycentre xMMS of the samples (x 1 , . . . , x N ) converges almost-surely to xMMS (this is a result of Proposition 3.4).

Numerical experiments were conducted in the case where M is a hyperbolic space of curvature equal to -1 and of dimension d. The following table was obtained for the values σ 2 = τ 2 = 0.1, using samples (x 1 , . . . , x N ) where N = 2 × 10 Other values of σ 2 and τ 2 lead to similar orders of magnitude for m1 (x MAP ) and d(x MMS , xMAP ). While m1 (x MAP ) increases with the dimension d, d(x MMS , xMAP ) does not appear sensitive to increasing dimension.

Based on these experimental results, one is tempted to conjecture that xMMS = xMAP , even when ρ = 1/2. Of course, numerical experiments do not equate to a mathematical proof.

Computing the MMS

Metropolis-Hastings algorithm

A crucial step in Bayesian inference is sampling from the posterior density. Here, this is π(x) given by (41). Since π(x) is known only up to normalisation, a suitable sampling method is afforded by the Metropolis-Hastings algorithm. This algorithm generates a Markov chain (x n ; n ≥ 1), with transition kernel [START_REF] Roberts | General state-space Markov chains and MCMC algorithms[END_REF] P f (x) = M α(x, y)q(x, y)f (y)vol(dy) + ρ(x)f (x) (47)

for any bounded measurable function f : M → R, where α(x, y) is the probability of accepting a transition from x to dy, and ρ(x) is the probability of staying at x, and where q(x, y) is the proposed transition density q(x, y) ≥ 0 and M q(x, y)vol(dy) = 1 for x ∈ M (48)

In the following, (x n ) will always be an isotropic Metropolis-Hastings chain, in the sense that q(x, y) = q(d(x, y)), so q(x, y) only depends on the distance d(x, y). In this case, the acceptance probability α(x, y) is given by α(x, y) = min {1, π(y)/π(x)}.

The aim of the Metropolis-Hastings algorithm is to produce a Markov chain (x n ) which is geometrically ergodic. Geometric ergodicity means the distribution π n of x n converges to π, with a geometric rate, in the sense that there exist β ∈ (0, 1) and R(x 1 ) ∈ (0, ∞), as well as a function V : M → R, such that (in the following, π(dx) = π(x)vol(dx))

V (x) ≥ max 1, d 2 (x, x * ) for some x * ∈ M (49) M f (x)(π n (dx) -π(dx)) ≤ R(x 1 )β n (50)
for any function f : M → R with |f | ≤ V . If the chain (x n ) is geometrically ergodic, then it satisfies the strong law of large numbers [START_REF] Meyn | Markov chains and stochastic stability[END_REF] 1

N N n=1 f (x n ) -→ M f (x)π(dx) (almost-surely) (51)
as well as a corresponding central limit theorem (see Theorem 17.0.1, in [START_REF] Meyn | Markov chains and stochastic stability[END_REF]). Then, in practice, the Metropolis-Hastings algorithm can be used to generate samples (x n ) from the posterior density π(x).

The following general statement can be proved, concerning the geometric ergodicity of isotropic Metropolis-Hastings chains. The proof (see [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF], Section 4.6) is a generalisation of the one carried out in the special case where M is a Euclidean space [START_REF] Jarner | Geometric ergodicity of Metropolis algorithms[END_REF]. Proposition 3.3. Let M be a Riemannian symmetric space, which belongs to the noncompact case. Assume (x n ; n ≥ 1) is a Markov chain in M , with transition kernel given by (47), with proposed transition density q(x, y) = q(d(x, y)), and with strictly positive invariant density π.

The chain (x n ) satisfies ( 49) and (50), if the following assumptions hold, (a1) there exists x * ∈ M , such that r(x) = d(x * , x) and (x) = log π(x) satisfy lim sup r(x)→∞ grad r, grad x r(x) < 0 (a2) if n(x) = grad (x)/ grad (x) , then n(x) satisfies lim sup r(x)→∞ grad r, n x < 0 (a3) there exist δ q > 0 and ε q > 0 such that d(x, y) < δ q implies q(x, y) > ε q

Remark : the posterior density π in (41) verifies Assumptions (a1) and (a2). To see this, let x * = z, and write grad (x) = -1 τ 2 r(x)grad r(x) -

1 σ 2 grad f y (x)
where f y (x) = d 2 (y, x)/2. Then, taking the scalar product with grad r, grad r, grad x = -1 τ 2 r(x) -1 σ 2 grad r, grad f y x (52) since grad r(x) is a unit vector, for all x ∈ M . Now, grad f y (x) = -Exp -1 x (y) (see [START_REF] Chavel | Riemannian geometry, a modern introduction[END_REF]). But, since r(x) is a convex function of x, it follows, by (85) in Appendix B.1, that grad r, Exp -1

x (y) ≤ r(y) -r(x)

for any y ∈ M . Thus, the right-hand side of (52) is strictly negative, as soon as r(x) > r(y), and Assumption (a1) is indeed verified. That Assumption (a2) is also verified can be proved by a similar reasoning. Remark : on the other hand, Assumption (a3) holds, if the proposed transition density q(x, y) is a Gaussian density, q(x, y) = p(y|x, τ q ). With this choice of q(x, y), all the assumptions of Proposition 3.3 are verified, for the posterior density π in (41). Therefore, Proposition 3.3 implies that the Metropolis-Hastings algorithm generates geometrically ergodic samples (x n ; n ≥ 1), from this posterior density.

The empirical barycentre

Let (x n ; n ≥ 1) be a Metropolis-Hastings Markov chain in M , with its transition kernel (47), and invariant density π. Assume the chain (x n ) is geometrically ergodic, so it satisfies the strong law of large numbers (51). Then, let xN denote the empirical barycentre of the first N samples (x 1 , . . . , x N ). This is the unique global minimum of the variance function The proof of Proposition 3.4 is nearly a word-for-word repetition of the proof in [START_REF] Bhattacharya | Large sample theory of instrinsic and extrinsic sample means on manifolds I[END_REF] (that of Theorem 2.3).

E N (w) = 1 2N N n=1 d 2 (w, x n ) ( 53
According to the remarks after Proposition 3.3, the Metropolis-Hastings Markov chain (x n ), whose invariant density is the posterior density π(x), given by (41), is geometrically ergodic. Therefore, by Proposition 3.4, the empirical barycentre xMMS , of the samples (x 1 , . . . , x N ), converges almost-surely to the minimum mean square error estimator xMMS (since this is just the barycentre of the posterior density π). This provides a practical strategy for approximating xMMS . Indeed, xMMS can be computed using the Riemannian gradient descent method (this method is discussed in Appendix B.4). Remark : this strategy for approxmating xMMS provided the numerical results discussed in Paragraph 3.2. For an additional, visual illustration, consider (as in Paragraph 3.2) the case where M is a space of constant negative curvature -1, and of dimension d = 2. Figure 1 Neither this proposition, nor its proof, appeared in [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF]. The proof is here given in a series of lemmas. Recall that M is now a hyperbolic space (simply-connected space of constant negative curvature -1).

Lemma 3.1. Let γ : R → M denote the geodesic curve with γ(0) = z and γ(1) = y. Then, xMMS lies on this geodesic curve γ.

Proof : recall from [START_REF] Alekseevskij | Geometry of spaces of constant curvature[END_REF] (Section 3) that there exists an isometry σ : M → M , such that σ • σ is the identity (σ is an involution), and the set of fixed points of σ is exactly the geodesic curve γ. The key point in the following, which can be seen from ( 41), is that

π(σ(x)) = π(x) for x ∈ M (54)
In other words, σ leaves invariant the posterior density π. Let E π be the function in (43).

Then, note that

(E π • σ)(w) = 1 2 M d 2 (w, σ(x))π(x)vol(dx) = 1 2 M d 2 (w, x)π(σ(x))σ * (vol)(dx)
where the first equality follows from (41), because σ is an isometry and an involution, and the second equality by a change of variables (σ * (vol) denotes the pullback of the volume form vol by σ). Using (54) and the fact that σ preserves the volume, it now follows that

(E π • σ)(w) = E π (w) for w ∈ M (55)

A Riemannian symmetric spaces

A Riemannian symmetric space is a Riemannian manifold M , such that, for each x ∈ M , there exists an isometry s x : M → M , with s x (x) = x and ds x (x) = -Id x . This isometry s x is called the geodesic symmetry at x [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF]. Let G denote the identity component of the isometry goup of M , and K = K o be the stabiliser in G of some point o ∈ M . Then, M = G/K is a Riemannian homogeneous space. The mapping θ :

G → G, θ(g) = s o • g • s o is an involutive isomorphism of G.
Let g denote the Lie algebra of G, and consider the Cartan decomposition, g = k + p, where k is the +1 eigenspace of dθ and p is the -1 eigenspace of dθ. One clearly has the commutation relations,

[k, k] ⊂ k ; [k, p] ⊂ p ; [p, p] ⊂ k (62)
In addition, it turns out that k is the Lie algebra of K, and that p may be identified with T o M , in a natural way.

The Riemannian metric of M may always be expressed in terms of an Ad(K)-invariant scalar product Q on g. If x ∈ M is given by x = g • o for some g ∈ G (where

g • o = g(o)), then u,v x = Q(g -1 • u, g -1 • v) ( 63
)
where the vectors g -1 • u and g -1 • v, which belong to T o M , are identified with elements of p. Here, by an abuse of notation, dg

-1 • u is denoted g -1 • u. Let exp : g → G denote the Lie group exponential. If v ∈ T o M , then the Riemannian exponential Exp o (v) is given by Exp o (v) = exp(v) • o (64)
Moreover, if Π t 0 denotes parallel transport along the geodesic c(t) = Exp o (tv), then

Π t 0 (u) = exp(tv) • u (65)
for any u ∈ T o M (note that the identification T o M p is always made, implicitly). Using (65), one can derive the following expression for the Riemann curvature tensor at o,

R o (v, u)w = -[[v , u], w] v, u, w ∈ T o M (66)
A fundamental property of symmetric spaces is that the curvature tensor is parallel :∇ R = 0. This is often used to solve the Jacobi equation [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF][14], and then express the derivative of the Riemannian exponential,

dExp x (v)(u) = exp(v) • sh(R v )(u) ( 67
)
where sh

(R v ) = ∞ n=0 (R v ) n /(2n + 1)! for the self-adjoint curvature operator R v (u) = [v , [v, u]].
As a result of (67), since exp(v) is an isometry, the following expression of the Riemannian volume is immediate

Exp * o (vol) = |det(sh(R v ))| dv ( 68
)
where dv denotes the volume form on T o M , associated with the restriction of the scalar product Q to p.

Expression (68) yields applicable integral formulae, when g is a reductive Lie algebra (g = z + g ss : z the centre of g and g ss semisimple). If a is a maximal Abelian subspace of p, any v ∈ p is of the form v = Ad(k) a for some k ∈ K and a ∈ a (see [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF], Lemma 6.3, Chapter V). Moreover, using the fact that Ad(k) is an isomorphism of g,

Ad(k -1 ) • R v • Ad(k) = R a = λ∈∆ + (λ(a)) 2 Π λ ( 69
)
where each λ ∈ ∆ + is a linear form λ : a → R, and Π λ is the orthogonal projector onto the corresponding eigenspace of R a . Here, ∆ + is the set of positive roots of g with respect to a [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF] (see Lemma 2.9, Chapter VII). It is possible to use the diagonalisation (69), in order to evaluate the determinant (68). To obtain a regular parameterisation, let S = K/K a , where K a is the centraliser of a in K. Then, let ϕ : S ×a → M be given by ϕ(s, a) = Exp o (β(s, a)) where β(s, a) = Ad(s) a. Now, by ( 68) and (69),

ϕ * (vol) = λ∈∆ + sinh λ(a) λ(a) m λ β * (dv)
where m λ is the multiplicity of λ (the rank of Π λ ). On the other hand, one may show

β * (dv) = λ∈∆ + |λ(a)| m λ da ω(ds) ( 70
)
where da is the volume form on a, and ω is the invariant volume induced onto S from K.

Finally, the Riemannian volume, in terms of the parameterisation ϕ, can be expressed in the following way ϕ * (vol) =

λ∈∆ + |sinh λ(a)| m λ da ω(ds) (71)
Using (71), it will be possible to write down integral formulae for Riemannian symmetric spaces, either non-compact or compact.

A.1 The non-compact case

This is the case were g admits an Ad(G)-invariant, non-degenerate, symmetric bilinear form B, such that Q(u, z) = -B(u, dθ(z)) is an Ad(K)-invariant scalar product on g. In this case, B is negative-definite on k and positive-definite on p. Moreover, the linear map ad(z) = [z, •] is skew-symmmetric or symmetric (with respect to Q), according to whether z ∈ k or z ∈ p.

If u 1 , u 2 ∈ p are orthonormal, the sectional curvature of Span(u 1 , u 2 ) is found from (66

), κ(u 1 , u 2 ) = -[u 1 , u 2 ] 2 o ≤ 0.
Therefore, M has non-positive sectional curvatures. In fact, M is a Hadamard manifold. It is geodesically complete by (64). It is moreover simply connected, because Exp o : p → M is a diffeomorphism [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF] (Theorem 1.1, Chapter VI). Thus, (68) yields a first integral formula,

M f (x) vol(dx) = p f (Exp o (v)) |det(sh(R v ))| dv (72)
To obtain an integral formula from (71), one should first note that β : S × a → p is not regular, nor one-to-one. Recall the following :

• the hyperplanes λ(a) = 0, where λ ∈ ∆ + , divide a into finitely many connected components, which are open and convex sets, known as Weyl chambers. From (70), β is regular on each Weyl chamber.

• let K a denote the normaliser of a in K. Then, W = K a /K a is a finite group of automorphisms of a, called the Weyl group, which acts freely transitively on the set of Weyl chambers [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF] (Theorem 2.12, Chapter VII). Then, for each Weyl chamber C, β is regular and one-to-one, from S ×C onto its image in p. Moreover, if a r is the union of the Weyl chambers (a ∈ a r if and only if λ(a) = 0 for any λ ∈ ∆ + ), then β is regular and |W |-to-one from S × a r onto its image in p. To obtain the desired integral formula, it only remains to note that ϕ is a diffeomorphism from S × C onto its image in M . However, this image is the set M r of regular values of ϕ. By Sard's lemma, its complement is negligible [START_REF] Bogachev | Measure Theory[END_REF].

Proposition A.1. Let M = G/K be a Riemannian symmetric space, which belongs to the "non-compact case", just described. Then, for any bounded continuous function

f : M → R, M f (x) vol(dx) = C + S f (ϕ(s, a)) λ∈∆ + (sinh λ(a)) m λ da ω(ds) (73) = 1 |W | a S f (ϕ(s, a)) λ∈∆
+ |sinh λ(a)| m λ da ω(ds) (74)
Here, C + is the Weyl chamber C + = {a ∈ a : λ ∈ ∆ + ⇒ λ(a) > 0}.

A.2 The compact case

In this case, g admits an Ad(G)-invariant scalar product Q. Therefore, ad(z) is skewsymmmetric, with respect to Q, for each z ∈ g. Using (66), it follows that M is compact, with non-negative sectional curvature.

In fact, the compact case may be obtained from the previous non-compact case by duality. Denote g C the complexification of g, and let g * = k + p * where p * = ip. Then, g * is a compact real form of g C (that is, g * is a compact Lie algebra, and its complexification is equal to g C ). Denote G * the connected Lie group with Lie algebra g * .

If M = G/K is a Riemannian symmetric space which belongs to the non-compact case, then M * = G * /K is a Riemannian symmetric space which belongs to the compact case. Formally, to pass from the non-compact case to the compact case, all one has to do is replace a by ia. Applying this recipe to (71), one obtains

ϕ * (vol) = λ∈∆ + |sin λ(a)| m λ da ω(ds) ( 75
)
where da is the volume form on a * = ia, and ω is the invariant volume induced onto S from K. Note that the image under Exp o of a * is the torus T * = a * /a K , where a K is the lattice given by a K = {a ∈ a * : Exp o (a) = o}. Recall the following :

• ϕ(s, a) only depends on t = Exp o (a). Thus, ϕ may be considered as a map from S × T * to M .

• if a ∈ a K then exp(2a) = e (the identity element in G * ). Thus, λ(a) ∈ iπ Z for all λ ∈ ∆ + [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF] (Page 383). Therefore, there exists a function D : T → R, such that Proposition A.2. Let M = G * /K be a Riemannian symmetric space, which belongs to the "compact case", just described. For any bounded continuous function f :

D(t) =
M → R, M f (x) vol(dx) = 1 |W | T * S f (ϕ(t, a))D(t)dt ω(ds) (77)
A.3 Example of Propositions A.1 and A.2

consider M = H(N ) the space of N × N Hermitian positive-definite matrices. Here, G = GL(N, C) and K = U (N ). Moreover, B(u,z) = Re(tr(uz)) and dθ(z) = -z † . Thus, p is the space of N × N Hermitian matrices, and one may choose a the space of N × N real diagonal matrices. The positive roots are the linear maps λ(a) = a ii -a jj where i < j, and each one has its multiplicity equal to 2. The Weyl group W is the group of permutation matrices in U (N ) (so |W | = N !). Finally, S = U (N )/T N ≡ S N , where T N is the torus of diagonal unitary matrices. By (74),

H(N ) f (x) vol(dx) = 1 N ! a S N f s exp(2a)s † i<j sinh 2 (a ii -a jj )da ω(ds) (78)
where da = da 11 . . . da N N . Now, assume f is a class function :

f (k • x) = f (x) for k ∈ K and x ∈ H(N ).
That is, f (x) depends only on the eigenvalues x i = e r i of x. By (78),

H(N ) f (x) vol(dx) = ω(S N ) 2 N N ! R N f (exp(r)) i<j sinh 2 ((r i -r j )/2)dr (79)
or, by introducing the eigenvalues x i as integration variables,

H(N ) f (x) vol(dx) = ω(S N ) 2 N 2 N ! R N + f (x 1 , . . . , x N ) |V (x)| 2 N i=1 x -N i dx i (80)
where

V (x) = i<j (x j -x i ) is the Vandermonde determinant.
The dual of H(N ) is the unitary group U (N ). Here, G * = U (N ) × U (N ) and K U (N ), is the diagonal group K = {(x, x) ; x ∈ U (N )}. The Riemannian metric is given by the trace scalar product Q(u,z) = -tr(uz). Moreover, T * = T N and S = S N (this is U (N )/T N ). The positive roots are λ(ia) = a ii -a jj where i < j and where a is N × N , real and diagonal 4 . By writing the integral over T N as a multiple integral, (77) reads,

U (N ) f (x) vol(dx) = 1 N ! [0,2π] N S N f s exp(2ia)s † i<j sin 2 (a ii -a jj )ω(ds) da (81)
where da = da 11 . . . da N N . Now, assume f is a class function. That is, f (x) depends only on eigenvalues e iθ i of x. Integrating out s, from (81), it follows,

U (N ) f (x) vol(dx) = ω(S N ) 2 N N ! [0,2π] N f (exp(iθ)) i<j sin 2 ((θ i -θ j )/2)dθ (82)
or, after an elementary manipulation,

U (N ) f (x) vol(dx) = ω(S N ) 2 N 2 N ! [0,2π] N f (θ 1 , . . . , θ N ) |V (e iθ )| 2 dθ 1 . . . θ N (83)
where V (e iθ ) = i<j (e iθ j -e iθ i ) is the Vandermonde determinant. Integrals such as (80) and ( 83) are familiar in random matrix theory [17][34]. The resemblance between these integrals (for example, in the rôle played by the Vandermonde determinant) is at the origin of the sort of "duality" described in Paragraph 2.8.

B Convex optimisation B.1 Convex sets and functions

In Euclidean geometry, a convex set A is any set which satisfies the definition : if points x and y belong to A, then the straight line segment between x and y lies entirely in A.

One hopes to extend this definition to Riemannian geometry, by letting geodesics play the rôle of straight lines. However, this does not lead to one, but to multiple definitions of a convex set. The present article will focus on the following [START_REF] Chavel | Riemannian geometry, a modern introduction[END_REF].

Definition B.1. A subset A of a complete Riemannian manifold M is called strongly convex if, whenever points x and y belong to A, there exists a unique length-minimising geodesic γ x,y connecting x and y, and γ x,y lies entirely in A. Remark : as an example of a different way of defining a convex set, consider the following. A subset A of M is called weakly convex if, whenever points x and y belong to A, there exists a unique geodesic γ in M , such that γ connects x and y, and γ lies entirely in A (this coincides with the definition in [START_REF] Chavel | Riemannian geometry, a modern introduction[END_REF], because γ is then the unique length-minimising curve, among all curves that connect x and y and lie entirely in A). In Euclidean geometry, a ball of any radius is convex. In a Riemannian manifold, a ball may fail to be strongly (or even weakly) convex, if its radius is too large. On the other hand, a ball with sufficiently small radius is always strongly convex [START_REF] Chavel | Riemannian geometry, a modern introduction[END_REF] [START_REF] Petersen | Riemannian geometry[END_REF]. Remark : if 1 2 inj(M ) is replaced by inj(M ) in (84), then R ≤ R c (M ) implies B(x, R) is weakly convex [START_REF] Chavel | Riemannian geometry, a modern introduction[END_REF] [START_REF] Kendall | Probability, convexity, and harmonic maps with small image I : uniqueness and fine existence[END_REF]. There is a certain class of Riemannian manifolds, where balls of any radius are strongly convex. Namely, these are Hadamard manifolds. Recall that a Hadamard manifold is a simply connected, complete Riemannian manifold with non-positive sectional curvatures.

In particular [START_REF] Lee | Introduction to smooth manifolds[END_REF], this implies inj(M ) = ∞ and κ max = 0, so that R c (M ) = ∞.

The following definition of a convex function on a Riemannian manifold directly extends the usual, well-known definition of a convex function on a Euclidean space. Definition B.2. Let A be a strongly convex subset of a complete Riemannian manifold M , and f : A → R. Then, f is called convex (respectively, strictly convex) if f (γ x,y (t)) is a convex (respectively, strictly convex) function of the time parameter t, for all x, y ∈ A. Further, if there exists α > 0 such that f (γ x,y (t)) is an α-strongly convex function of t, for all x, y ∈ A, then f is called α-strongly convex.

For differentiable functions, it is possible to write down first-order and second-order characterisations of convexity [START_REF] Udriste | Convex functions and optimization methods on Riemannian manifolds[END_REF]. Recall that γx,y (0) = Exp -1 x (y), for x, y ∈ A, where the dot denotes the time derivative and Exp the Riemannian exponential map [START_REF] Lee | Introduction to smooth manifolds[END_REF]. In addition, let gradf and Hessf denote the gradient and Hessian of a function f on M (defined with respect to the Riemannian metric and Levi-Civita connection of M ).

Proposition B.2. Let A be a strongly convex subset of a complete Riemannian manifold M , and f : A → R.

(i) assume f is differentiable. Then, f is convex, if and only if

f (y) -f (x) ≥ gradf (x), Exp -1 x (y) x for all x, y ∈ A (85)
Moreover, f is strictly convex if and only if the above inequality is strict whenever y = x.

(ii) assume f is differentiable. Then f is α-strongly convex, if and only if,

f (y) -f (x) ≥ gradf (x), Exp -1 x (y) x + (α/2)d 2 (y, x) for all x, y ∈ A (86)
(iii) assume f is twice differentiable. Then f is convex if and only if Hessf (x) 0, and strictly convex if and only if Hessf (x) 0, for all x ∈ A. Moreover, f is α-strongly convex if and only if Hessf (x) αg(x), for all x ∈ A.

Here, •, • and d(•, •) denote the Riemannian scalar product and distance, associated with the Riemannian metric tensor g of M . Moreover, stands for the Loewner order. A straightforward consequence of (ii) in Proposition B.2 is the so-called PL inequality (PL stands for Polyak-Lojasiewicz [START_REF] Karimi | Linear convergence of gradient and proximalgradient methods under the Polyak-Lojasiewicz condition[END_REF]). This inequality will be used in Paragraph B.4.2.

Proposition B.3. Let f : A → R be a twice differentiable, α-strongly convex function. If f has its minimum at x * ∈ A, then gradf (x) 2 x ≥ 2α(f (x) -f (x * )) for all x ∈ A (87)

B.2 Second-order Taylor formula

Consider the second-order Taylor formula, for a twice-differentiable function f : M → R (as usual, M is a complete Riemannian manifold). For x ∈ M and v ∈ T x M ,

f (Exp x (v)) = f (x) + gradf (x), v x + 1 2 Hessf γ(t * ) ( γ , γ) ( 88
)
where γ is the geodesic curve γ(t) = Exp x (tv) and t * ∈ (0, 1). Formula (88) is the firstorder Taylor expansion, with Lagrange remainder, of the function f (γ(t)), at t = 0 [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF]. This formula will be the starting point for the study of Riemannian gradient descent in Paragraph B.4. There, it will be applied with v = -µgradf (x) and µ ∈ (0, 1].

To apply (88), it is quite helpful to control the second-order term in its right-hand side. One says that f is L-smooth on B ⊂ M , if there exists L ≥ 0 with |Hessf y (u, u)| ≤ L u 2 y for all y ∈ B and u ∈ T y M . Then, if γ(t) = Exp x (tv) belongs to B for all t ∈ (0, 1), Remark : as a consequence of (90), if x * ∈ B c is such that f (x * ) is the minimum of f (x), taken over x ∈ B c , then

f (Exp x (v)) ≤ f (x) + gradf (x), v x + (L/2) v 2
2L c (f (x) -f (x * )) ≥ gradf (x) 2 x for all x ∈ B c (91)
which is complementary to (87).

B.3 Taylor with retractions

It is customary, in practical applications, to approximate the Riemannian exponential map by another so-called retraction map, which is easier to compute [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. In this context, it is helpful to derive new versions of Formula (88) and of Proposition B.4, which apply when the exponential map Exp is replaced with a retraction Ret.

Recall that a retraction is a smooth map Ret :

T M → M (denoted Ret(x, v) = Ret x (v) for x ∈ M and v ∈ T x M ), such that Ret x (0 x ) = x and dRet x (0 x ) = Id x (92)
for all x ∈ M . Here, 0 x is the zero element in T x M and Id x is the identity map of T x M . Most retractions, encountered in practical applications, are regular retractions, in the following sense [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF].

Definition B.3. A retraction Ret : T M → M is regular, if there exists a smooth bundle map Φ : T M → T M , such that

Ret x (v) = Exp x (Φ x (v)) for all x ∈ M and v ∈ T x M (93) Here, Φ is denoted Φ(x, v) = Φ x (v) ("bundle map" means Φ x (v) ∈ T x M for all v ∈ T x M ).
If Ret : T M → M is a regular retraction and f : M → R is a twice-differentiable function, then (88) and (93) directly imply

f (Ret x (v)) = f (x) + gradf, Φ x (v) x + 1 2 Hessf γ(t * ) ( γ, γ) ( 94
)
where γ is the geodesic curve γ(t) = Exp x (tΦ x (v)) and t * ∈ (0, 1). Formula (94) is the required new version of (88). The Riemannian exponential Exp is a regular retraction, with Φ x = Id x for x ∈ M . For a general regular retraction Ret, each map Φ x : T x M → T x M still agrees with Id x up to second-order terms [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF].

Proposition B.5. Let Ret : T M → M be a regular retraction, with Φ : T M → T M given by (93). Then, for each x ∈ M , the map Φ x : T x M → T x M verifies (a) Φ x (0 x ) = 0 x and Φ x (0 x ) = Id x (the prime denotes the Fréchet derivative). (b) Φ x (0 x )(v, v) = c(0), where the curve c(t) is given by c(t) = Ret x (tv).

The retraction Ret is called geodesic when Φ x (0 x )(v, v) = 0 for x ∈ M and v ∈ T x M . In this case, Φ x agrees with Id x up to third order terms. For geodesic regular retractions, the following Proposition B.6 provides a new, general version of Proposition B.4. Proposition B.6. Let f : M → R be a C 2 function, with B c = {x : f (x) ≤ c} compact (and not empty), and let Ret : T M → M be a geodesic regular retraction. There exist constants β c , δ c , H c ≥ 0, which depend on f and Ret, such that, for all x ∈ B c ,

f (y) ≤ f (x) -µ 1 -(β c H c /2)µ -(δ c gradf (x) 2 x )µ 2 gradf (x) 2 x ( 95
)
whenever y = Ret x (-µgradf (x)) for some µ ∈ (0, 1]. In particular,

1 2 -(β c H c /2)µ -(δ c gradf (x) 2 x )µ 2 ≥ 0 =⇒ f (y) ≤ f (x) -(µ/2) gradf (x) 2 x (96) Therefore, x ∈ B c implies y ∈ B c .
Remark : the application of Proposition B.6 is somewhat simplified when the retraction Ret, in addition to being regular and geodesic, is contractive and uniformly geodesic.

Here, contractive means that

Φ x (v) x ≤ v x for x ∈ M and v ∈ T x M (97)
In this case, it is always possible to put β c = 1 and H c = L c , where L c is the same constant as in Proposition B.4. Uniformly geodesic means there exists δ ≥ 0, such that

Φ x (v) -v x ≤ δ v 3 x for x ∈ M and v ∈ T x M (98)
In this case, it is always possible to put δ c = δ (independent of c and even of f ). The widely-used projection retractions for spheres, unitary groups and Grassmann manifolds, are examples of contractive, uniformly geodesic regular retractions [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF] (see Sections 1.5 and 1.6). Here is another example, for positive-definite matrices.

Example : let M = P(N ), the space of symmetric positive-definite N × N matrices, equipped with its usual affine-invariant metric [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF]. For x ∈ P(N ), the tangent space T x P(N ) is identified with the space S(N ) of symmetric N × N matrices. Then, recall the Riemannian exponential map Exp x (v) = x exp(x -1 v) (where exp denotes the matrix exponential), and consider the retraction Ret x (v) = x + v + (1/2)vx -1 v. The point of using this retraction is that the eigenvalues of Ret x (v) will always be greater than 1/2. In addition, it is a contractive, uniformly geodesic regular retraction.

B.4 Riemannian gradient descent

Let f : M → R be a C 2 function, and Ret : T M → M a retraction. Together, these yield the Riemannian gradient descent scheme, where µ ∈ (0, 1] is called the step-size,

x t+1 = Ret x t (-µgradf (x t )) t = 0, 1, . . . (99)
Assume f has a compact (non-empty) sublevel set B c , and choose a constant L c ≥ 0 as in Proposition B.4. Also, assume Ret is a contractive, uniformly geodesic regular retraction, as in the remark after Proposition B.6. Specifically, Ret verifies (97) and (98), for some constant δ ≥ 0. This allows for Ret = Exp, the Riemannian exponential, in which case δ = 0. Now, the following lemma is a direct consequence of (96) in Proposition B.6.

Lemma B.1. Under the two assumptions just described, on f and Ret, if x 0 ∈ B c , then

1 2 -(L c /2)µ -(δ gradf 2 Bc )µ 2 ≥ 0 =⇒ f (x t+1 ) ≤ f (x t ) -(µ/2) gradf (x t ) 2 x t (100)
for all t ≥ 0, so that x t ∈ B c for all t ≥ 0. Here, gradf Bc is the maximum of gradf (x) x , taken over x ∈ B c .

This lemma immediately yields the convergence of the Riemannian gradient descent scheme (99), to the stationary points of f in B c . Here, µ * c is the infimum of µ ∈ (0, 1] such that 1 2 -(L c /2)µ -(δ gradf 2 Bc )µ 2 < 0. Proposition B.7. Under the assumptions on f and Ret, made in Lemma B.1, if µ ≤ µ * c , then x 0 ∈ B c implies the sequence (x t ) generated by (99) converges to the set of stationary points of f , in the sublevel set B c .

The proof of this proposition is straightforward.

From Lemma B.1, if µ ≤ µ * c , then x 0 ∈ B c implies x t ∈ B c and (µ/2) gradf (x t ) 2 x t ≤ f (x t ) -f (x t+1
) for all t ≥ 0. Adding these inequalities, for t = 0, . . . , T , (µ/2)

T t=0 gradf (x t ) 2 x t ≤ f (x t ) -f (x T +1 )
Then, since f is bounded below on the compact set B c , the series ∞ t=0 gradf (x t ) 2

x t must converge. Finally, compactness of B c ensures every subsequence of (x t ) has a further subsequence that converges to a stationary point of f in B c . Proposition B.7 holds without any convexity assumptions, made on the function f . Consider now the case of a strictly convex, and then of a strongly convex f .

B.4.1 Strictly convex case

Now, assume the function f is strictly convex on some strongly convex subset A of M . Moreove, assume f has a compact (non-empty) sublevel set B c ⊂ A. Then, choose a constant L c ≥ 0 as in Proposition B.4, and let the retraction Ret be as in Lemma B.1.

Assume that f has a unique minimum at x * ∈ B c . Let R be the radius of the smallest ball B(x * , R) such that B c ⊂ B(x * , R), D the maximum of gradf (x) x for x ∈ B(x * , R), and R c = R + D. The following Lemma ensures that (99) "contracts the distance to x * ". Lemma B.2. Assume that the sectional curvatures of M lie in the interval [-κ 2 min , κ 2 max ], and that R c < inj(x * ). If x 0 ∈ B c and µ ≤ µ * c (with µ * c as in Proposition B.7), then

1 -κ(R c )L c µ -(2δR c gradf Bc )L c µ 2 ≥ 0 =⇒ d(x t+1 , x * ) ≤ d(x t , x * ) (101) where κ(R c ) = max {κ min R c coth(κ min R c ), |κ max R c cot(κ max R c )|} and gradf Bc is the maximum of gradf (x) x for x ∈ B c .
Remark : inj(x * ) denotes the injectivity radius of M at x * [START_REF] Petersen | Riemannian geometry[END_REF]. The condition that R c < inj(x * ) guarantees the x t stay away from the cut locus Cut(x * ). This condition is introduced because the distance function x → d(x, x * ) is not differentiable on Cut(x * ), where its Hessian may even diverge to -∞. Let µ * d denote the infimum of µ such that 1 -κ(R c )L c µ -(2δR c gradf Bc )L c µ 2 < 0.

Proposition B.9. Under the assumptions just described, if µ ≤ µ * c and x 0 ∈ B c , then f (x t ) -f (x * ) ≤ (1 -µα) t (f (x 0 ) -f (x * )) (105) for all t ≥ 0. In particular, the sequence (x t ) converges to x * .

The proof of this proposition follows by replacing Inequality (87) into Lemma B.1. Indeed, if µ ≤ µ * c and x 0 ∈ B c , then (100) in Lemma B.1 immediately implies f (x t+1 ) -f (x * ) ≤ f (x t ) -f (x * ) -(µ/2) gradf (x t ) 2

x t for t ≥ 0 Thus, replacing (87) into the right-hand side,

f (x t+1 ) -f (x * ) ≤ (1 -µα)(f (x t ) -f (x * ))
and ( 105) can be obtained by induction.

Remark : (105) shows that f (x t ) converges to the minimum f (x * ), at an exponential rate. In practice, this can still be quite slow, if µα is very small. Indeed, one should attempt to use µ as large as possible, in order to benefit from the exponential rate (105).

From the definition of µ * c , one cannot have µ any larger than 1/L c , and µ = 1/L c is only possible if δ = 0, which corresponds to using Ret = Exp. Example : let π be a probability distribution on a complete Riemannian manifold M , and define

E π (x) = 1 2 M d 2 (x, y) π(dy) for x ∈ M (106)
If the support of π is contained in a ball B(z, R), where R < R c (M ) given by (84), then E π is C 2 on B(z, R), and has a unique global minimum x * ∈ M , such that x * ∈ B(z, R) [START_REF] Afsari | Riemannian L p center of mass : existence, uniqueness and convexity[END_REF] (x * is the Riemannian barycentre of π). In addition, if R < R c (M )/2, then E π is αstrongly convex on B(z, R), with α equal to 2κ max R cot(2κ max R) (= 1 if κ max = 0). In this case, it is possible to apply Proposition B.9 to the present example (with f = E π ). If M has positive sectional curvatures, it is always possible to choose L c = 1. On the other hand, if M has negative sectional curvatures L c = 1 + 4κ min R always works. Now, in order to prove (95), note that γ(t) = Exp x (tΦ x (-µgradf (x))) belongs to B for all t ∈ (0, 1). It follows from (94) that (similarly to (89)), f (y) ≤ f (x) + gradf, Φ x (-µgradf (x)) x + (H c /2) Φ x (-µgradf (x)) 2

x Then, using (109) and (110), f (y) ≤ f (x) -µ gradf (x) 2

x + (β c H c /2)µ 2 gradf (x) 2

x + δ c µ 3 gradf (x) 4

x which is the same as (95). Finally, (96) is an immediate consequence of (95).

Proof of Lemma B.1 : (100) can be obtained immediately, upon replacing β c = 1, δ c = δ and H c = L c into (96).

Proof of Proposition B.7 : the proof has already been summarised, right after the proposition.

Proof of Lemma B.2 : let L(x) = d 2 (x, x * )/2. If R c < inj(x * ), then L(x) is κ(R c )smooth on B(x * , R c ) [START_REF] Petersen | Riemannian geometry[END_REF]. Note that x t ∈ B c for all t ≥ 0, because µ ≤ µ * c as in Proposition B.7. Then, from the Taylor expansion (94) of L, and since Ret is contractive ((109) holds with β c = 1), L(x t+1 ) ≤ L(x t ) + gradL(x t ), Φ x t (-µgradf (x t )) x t + (κ(R c )/2)µ 2 gradf (x t ) 2 xt However, applying (91) to the third term on the right-hand side, this implies L(x t+1 ) ≤ L(x t ) + gradL(x t ), Φ x t (-µgradf (x t )) x t + κ(R c )L c µ 2 (f (x t ) -f (x * )) (111) Now, consider the second term on the right-hand side, since gradL(x t ) = -Exp -1

x t (x * ), this second term is equal to µ Exp -1

x t (x * ), gradf (x t ) x t -Exp -1 x t (x * ), Φ x t (-µgradf (x t )) + µgradf (x t ) x t (112)

Applying ( 85) and (110) (with δ c = δ, since Ret is uniformly geodesic), (112) ≤ -µ(f (x t ) -f (x * )) + (δµ 3 ) Exp -1 x t (x * ) x t gradf (x t ) 3

x t

Using (91) once again, along with Exp -1 x t (x * ) x t ≤ R c and gradf (x t ) x t ≤ gradf Bc , Since f (x t ) ≥ f (x * ), whenever the expression in square brackets is positive, one has L(x t+1 ) ≤ L(x t ). However, this directly yields (101). Proof of Proposition B.8 : note from Lemma B.1 that µ ≤ µ * c implies f (x t+1 ) -f (x * ) ≤ f (x t ) -f (x * ) -(µ/2) gradf (x t ) 2

x t

On the other hand, note that gradf (x t )

x t ≥ f (x t ) -f (x * ) d(x t , x * ) ≥ f (x t ) -f (x * ) d(x 0 , x * )
where the first inequality follows by applying Cauchy-Schwarz to (85), and the second one from Lemma B.2, since µ ≤ µ * d . Letting ε(t) = f (x t ) -f (x * ), it is now clear that ε(t + 1) ≤ ε(t) -(µ/2) ε(t) d(x 0 , x * ) 2 so that (102) can be proved by a straightforward induction. Proof of Proposition B.9 : the proof was summarised after the proposition.

where 1

 1 [a(t),b(t)] denotes the indicator function of the interval [a(t), b(t)].

2 .

 2 Paragraph 3.1 begins by expressing the posterior density π(x), based on the general definition π(x) ∝ prior density × likelihood density Here, π(x) will remain partially unknown, as the missing normalising factor cannot be determined. Then, two Bayesian estimators are studied. The maximum a posteriori xMAP is the mode of π(x), xMAP = argmax x∈M π(x)

)

  Let x denote the Riemannian barycentre of the invariant density π. It turns out that xN converges almost-surely to x. Proposition 3.4. Let (x n ) be any Markov chain in a Hadamard manifold M , with invariant distribution π. Denote xN the empirical barycentre of (x 1 , . . . , x N ), and x the Riemannian barycentre of π. If (x n ) satisfies the strong law of large numbers (51), then xN converges to x, almost-surely.

  represents M in the shape of the Poincaré disc. The prior barycentre z is designated by a square and the observation y by a circle •. Grey crosses × mark the last 1000 out of N = 100000 samples x n generated using the Metropolis-Hastings kernel (47), and the empirical barycentre xMMS is designated by a black circle •. In both of the Subfigures 1a and 1b, xMMS is seen to lie on the geodesic connecting z and y, here the dashed circle arc. Note that 1a corresponds to Proposition 3.1 and 1b to Proposition 3.2.

(a) σ 2 =Figure 1 :

 21 Figure 1: Poincaré disc with z = ; y = • ; xMMS = •

  λ∈∆ + |sin λ(a)| m λ whenever t = Exp o (a) Now, T * is a totally flat submanifold of M . Therefore, Exp * (dt) = da, where dt denotes the invariant volume induced onto T * from M . With a slight abuse of notation, (75) now reads, ϕ * (vol) = D(t)dt ω(ds) (76) Denote (T * ) r the set of t ∈ T * such that D(t) = 0. By the same arguments as in the non-compact case, ϕ is a regular |W |-to-one map from S × (T * ) r onto M r , the set of regular values of ϕ.

Proposition B. 1 .

 1 Assume the sectional curvatures of M are bounded above by κ 2 max ≥ 0. Then, denoting inj(M ) the injectivity radius of M , let R c (M ) For any x ∈ M and R < R c (M ), the open ball B(x, R) is strongly convex (if κ max = 0, it should be understood that division by zero yields infinity).

  yields the following Proposition B.4. To state this proposition, consider a C 2 function f : M → R, and assume the sublevel set B c = {x : f (x) ≤ c} is compact (and not empty), for some real c. Of course, B c is contained in some closed ball B = B(z, R). Let G be the maximum of gradf (x) x , taken over x ∈ B, and B = B(z, R + G). Now, by compactness of B , f is L c -smooth on B , for some L c ≥ 0. Proposition B.4. Let f : M → R be a C 2 function, with B c and L c defined as above, and let y = Exp x (-µgradf (x)) for some µ ∈ (0, 1]. If µ ≤ 1/L c , then f (y) ≤ f (x) -(µ/2) gradf (x) 2 x for all x ∈ B c (90) In particlar, x ∈ B c implies y ∈ B c .

C 2 x 2 x 2 x 2 c 2 x 3 x

 222223 Proofs for Section BProof of Proposition B.3 : write inequality (86) under the equivalent form f (y) -f (x) ≥ gradf (x), Exp-1 x (y) x + (α/2) Exp -1 x (y)With y = x * , this becomesf (x) -f (x * ) ≤ -gradf (x), Exp -1 x (x * ) x -(α/2) Exp -1 x (x * ) 2x or, by completing the square on the righ-hand side,f (x) -f (x * ) ≤ -Then, (87) follows immediately, by noting the first term on the righ-hand side is negative. Proof of Proposition B.4 : the proof employs the notation introduced before the proposition. Letx ∈ B c and v = -µgradf (x). Then, note that v x ≤ gradf (x) x ≤ G.This implies γ(t) = Exp x (tv) belongs to B for all t ∈ (0, 1). From the definition of L c , it now follows by (89) thatf (y) ≤ f (x) + gradf (x), v x + (L c /2) v 2x and, by recalling v = -µgradf (x),f (y) ≤ f (x) -µ(1 -(L c /2)µ) gradf (x)Then, (90) follows because µ ≤ 1/L c implies the expression in parentheses is ≥ 1/2. Proof of Proposition B.5 : this is given in[START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF], Section 1.5. Proof of Proposition B.6 : assume Ret is a regular geodesic retraction, and let Φ be the corresponding map in (93). Since B c is compact, there exist β c , δ c ≥ 0 such thatsup { Φ x (u) op ; x ∈ B c and u ∈ T x M, u x ≤ gradf (x) x } ≤ β 1/(107) sup { Φ x (u) op ; x ∈ B c and u ∈ T x M, u x ≤ gradf (x) x } ≤ δ c(108)where • op denotes the operator norm of the linear map Φ x (u) : T x M → T x M , or of the tri-linear map Φ x (u) :T x M × T x M × T x M → R.In terms of these constants β c and δ c ,Φ x (-µgradf (x)) 2 x ≤ (β c µ 2 ) gradf (x) for x ∈ B c (109) Φ x (-µgradf (x)) + µgradf (x) x ≤ (δ c µ 3 ) gradf (x) for x ∈ B c(110)Furthermore, let B c be contained in a closed geodesic ball B = B(z, R). Denote G the maximum of gradf (x) x taken over x ∈ B, and B = B(z, R + β 1/2 c G). By compactness of B , there exists H c ≥ 0 such that f is H c -smooth on B .

( 112 )

 112 ≤ -µ(f (x t ) -f (x * )) + (2δR c gradf Bc )L c µ 3 (f (x t ) -f (x * )) (113)Finally, from (111) and (113),L(x t+1 ) ≤ L(x t ) -µ 1 -κ(R c )L c µ -(2δR c gradf Bc )L c µ 2 (f (x t ) -f (x * ))

  5 .The first table confirms Proposition 3.1. The second table, more surprisingly, shows that xMMS and xMAP can be quite close to each other, even when ρ = 1/2.In both of these tables, d(x MMS , xMAP ) is an approximation of d(x MMS , xMAP ), based on using the empirical barycentre xMMS instead of xMMS . The main source of error affecting this approximation is the fact that the samples (x 1 , . . . , x N ) follow from a Metropolis-Hastings algorithm, and not directly from the posterior density π.

	dimension d 2	3	4	5	6	7	8	9	10
	m1 (x MAP ) 0.28 0.35 0.41 0.47 0.50 0.57 0.60 0.66 0.70
	d(x MMS , xMAP ) 0.00 0.00 0.00 0.01 0.01 0.02 0.02 0.02 0.03
	and the following table for σ 2 = 1 and τ 2 = 0.5, again using N = 2 × 10 5 .
	dimension d 2	3	4	5	6	7	8	9	10
	m1 (x MAP ) 0.75 1.00 1.12 1.44 1.73 1.97 2.15 2.54 2.91
	d(x MMS , xMAP ) 0.00 0.00 0.03 0.02 0.02 0.03 0.04 0.03 0.12

To follow the original notation of Jacobi[START_REF] Whittaker | A course of modern analysis[END_REF], this should be written ϑ(e iφ |q) where q = e -σ

. In other popular notations, this function is called ϑ 00 or ϑ

.

Proofs of the results stated in this section can be found in Chapter 4 of[START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF].

Please do no confuse the imaginary number i with the subscript i.

Finally, taking w = xMMS and recalling that xMMS is the unique global minimiser of E π , it follows that σ(x MMS ) = xMMS , so that xMMS indeed lies on γ. Lemma 3.2. There exists a continuous function c : R → R such that grad E π (γ(t)) = c(t) γ(t)

for t ∈ R (56)

Remark : taking covariant derivatives in (56),

where c (t) = dc(t)/dt. Because E π is 1-strongly convex (see Item (iii) of Proposition B.2, Appendix B.1), it follows that c (t) ≥ 1. In particular, c(t) is strictly increasing.

Proof : let w = γ(t) and take the gradient of (55). This yields

However, the derivative dσ : T γ(t) M → T γ(t) M is equal to 1 on vectors parallel to γ(t) and to -1 on vectors orthogonal to γ(t). Thus, (58) implies that grad E π (γ(t)) should be parallel to γ(t). This is equivalent to (56). The next step in the proof will be to compute c(0) and c(1). This will show that c(0) is negative and c(1) is positive. Computing c(0) and c(1) requires taking a closer look at the posterior density π.

Lemma 3.3. Let Z(z, y, τ, σ) denote the missing normalising factor in (41). Then, Z(z, y, τ, σ) = Z(δ, τ, σ), where δ = d 2 (z, y).

Proof : a hyperbolic space is a two-point homogeneous space [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF] (Page 355). This means if z , y ∈ M have d(z , y ) = d(z, y), then there exists an isometry g : M → M such g(z) = z and g(y) = y . Now, since g is an isometry,

Thus, introducing the change of variables w = g(x),

In other words, Z(z, y, τ, σ) only depends on the distance between z and y.

It is now possible to compute c(0) and c(1).

Lemma 3.4.

There exists a positive constant ψ(δ, τ, σ), such that

In the notation of (56).

Proof : for any value of the parameters (z, y, τ, σ),

where δ = d 2 (z, y). Taking the gradient of this identity with respect to z, and using

where grad z denotes the gradient with respect to z, it follows that

where ψ(δ, τ, σ) = -2 × ∂ log Z(δ, τ, σ)/∂δ. However, here one has,

Thus, replacing (61) into (60), it follows that

which is the first part of (59). The second part can be proved in the same way, taking the gradient with respect to y rather than z. The fact that ψ(δ, τ, σ) > 0 follows because c(t) is strictly increasing (see the remark after Lemma 3.2), and has at most one zero (because E π has exactly one stationary point).

It is now possible to complete the proof of Proposition 3.2. Lemma 3.4, shows that c(0) is negative and c(1) is positive. Therefore, c(t * ) = 0 for some t * ∈ (0, 1

for all t ≥ 0. In particular, the sequence (x t ) converges to x * .

Remark : the quality of the convergence in (102) depends above all on the step-size µ.

The smaller this is, the slower the convergence. From the definitions of µ * c and µ * d , it is clear that there are two reasons why µ would be smaller : a larger constant δ, and a larger curvature (in absolute value) κ min . In theory, one can alway make δ = 0 by using the retraction Ret = Exp, but this requires the ability to compute the Riemannian exponential Exp with sufficient accuracy. Remark : the rate of convergence stated in (102) is a partial generalisation of the rate found in [START_REF] Nesterov | Lectures on convex optimization[END_REF], for gradient descent in a Euclidean space. In the Euclidean setting, δ = 0 and κ min = κ max = 0. It then follows from Proposition B.8, that (102) obtains whenever µ ≤ 1/L c . Essentially, this is Corollary 2.1.2 (Page 81) in [START_REF] Nesterov | Lectures on convex optimization[END_REF]. However, note the restriction µ ∈ (0, 1], which is necessary in a curved Riemannian manifold.

Here is an optimisation problem, which falls under the scope of Proposition B.8. Example : let M be a Hadamard manifold, with sectional curvatures bounded below by -κ 2 min ≤ 0. Fix a cutoff parameter q > 0, and define

In [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF], it was proved that V y : M → R is strictly convex, but not strongly convex, and that it is (1 + qκ min )-smooth on M . Now, let π be a probability distribution on M and consider the problem of minimising

Note that V π is strictly convex (but not strongly convex), and (1 + qκ min )-smooth on M , because the same is true of each function V y . In fact, V π has compact sublevel sets whenever the distribution π has finite first-order moments [START_REF] Said | Statistical models and probabilistic methods on riemannian manifolds[END_REF]. In this case, V π (x) is guaranteed to achieve its minimum at some x * ∈ M . This x * is called the robust Riemannian barycentre of π (the adjective "robust" comes from the field of robust statistics [START_REF] Huber | Robust statistics[END_REF]). When applying Lemma B.2 and Proposition B.8 to the present example (with f = V π ), note that inj(x * ) = ∞, since M is a Hadamard manifold, and L c = (1 + qκ min ) does not depend on c.

B.4.2 Strongly convex case

Here, assume the function f is α-strongly convex on some strongly convex subset A ⊂ M . Let B c ⊂ A be a sublevel set of f (where c > inf x f (x)). Because f is strongly convex, B c is compact, and it is possible to choose a constant L c ≥ 0, as in Proposition B.4. Then, let µ * c be given as in Proposition B.7. As usual, f has a unique minimum at x * ∈ B c .