Xavier Pennec
email: xavier.pennec@sophia.inria.fr

Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements

In medical image analysis and high level computer vision, there is an intensive use of geometric features like orientations, lines, and geometric transformations ranging from simple ones (orientations, lines, rigid body or affine transformations, etc.) to very complex ones like curves, surfaces, or general diffeomorphic transformations. The measurement of such geometric primitives is generally noisy in real applications and we need to use statistics either to reduce the uncertainty (estimation), to compare observations, or to test hypotheses. Unfortunately, even simple geometric primitives often belong to manifolds that are not vector spaces. In previous works [1, 2], we investigated invariance requirements to build some statistical tools on transformation groups and homogeneous manifolds that avoids paradoxes. In this paper, we consider finite dimensional manifolds with a Riemannian metric as the basic structure. Based on this metric, we develop the notions of mean value and covariance matrix of a random element, normal law, Mahalanobis distance and χ 2 law. We provide a new proof of the characterization of Riemannian centers of mass and an original gradient descent algorithm to efficiently compute them. The notion of Normal law we propose is based on the maximization of the entropy knowing the mean and covariance of the distribution. The resulting family of pdfs spans the whole range from uniform (on compact manifolds) to the point mass distribution. Moreover, we were able to provide tractable approximations (with their limits) for small variances which show that we can effectively implement and work with these definitions.

1 In facts, geodesics are defined as the critical points of the energy functional E(γ) = 1 2 b a ∂γ 2 dt. It turns out that they also optimize the length functional but they are moreover parameterized proportionally to arc-length.

Introduction

To represent the results of a random experiment, one theoretically consider a probability measure on the space of all events. Although this probabilized space contains all the information about the random experiment, one often have access only to some measurements depending of the outcome of the experiment. The mathematical way to formalize this is to investigate random variables or observables which are maps from the probabilized space into R. One usually further simplify by restricting to random variables that have a probability density function (pdf).

However, from a computational point of view, the pdf is still too informative and we have to restrict the measurements to a few numeric characteristics of a random variable. Thus, one usually approximate a unimodal pdf by a central value and a dispersion value around it. The most used central value is the mean value or expectation of the random variable: x = E [ x ] = x dPr = R y p x (y) dy. The corresponding dispersion value is the variance σ 2

x = E (x -x) 2 . In real problems, we can have several simultaneous measurements of the same random experiment. If we arrange these n random variables x i into a vector x = (x 1 . . . x n ), we obtain a random vector. As the expectation is a linear operator, it is easily generalized to vector or matrix functions in order to define the mean value and the covariance matrix of a random vector: Σ xx = E [ (x -x) (x -x) T ] = (y -x) (y -x) T p x (y) dy. If one has to assume a probability distribution, the Gaussian distribution is usually well adapted, as it is completely determined by the mean and the covariance. It is moreover the entropy maximizing distribution knowing only these moments. Then, one can use a statistical distance between distributions such as the Mahalanobis distance and the associated statistical tests.

The problem we investigate in this article is to generalize this framework to measurements in finite dimensional Riemannian manifolds instead of measurements in a vector space. Examples of manifolds we routinely use in medical imaging applications are 3D rotations, 3D rigid transformations, frames (a 3D point and an orthonormal trihedron), semi-or non-oriented frames (where 2 (resp. 3) of the trihedron unit vectors are given up to their sign) [START_REF] Pennec | Landmark-based registration using features identified through differential geometry[END_REF][START_REF] Pennec | Feature-based registration of medical images: Estimation and validation of the pose accuracy[END_REF], oriented or directed points [START_REF] Granger | Rigid point-surface registration using an EM variant of ICP for computer guided oral implantology[END_REF][START_REF] Granger | Statistiques exactes et approchées sur les normales aléatoires[END_REF], positive definite symmetric matrices coming from diffusion tensor imaging [START_REF] Fletcher | Gaussian distributions on Lie groups and their application to statistical shape analysis[END_REF][START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors[END_REF][START_REF] Ch | Statistics on multivariate normal distributions: A geometric approach and its application to diffusion tensor MRI[END_REF][START_REF] Batchelor | A rigorous framework for diffusion tensor calculus[END_REF][START_REF] Pennec | A Riemannian framework for tensor computing[END_REF] or from variability measurements [START_REF] Fillard | Extrapolation of sparse tensor fields: Application to the modeling of brain variability[END_REF]. We have already shown in [START_REF] Pennec | A framework for uncertainty and validation of 3D registration methods based on points and frames[END_REF][START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF] that this is not an easy problem and that some paradoxes can arise. In particular, we cannot generalize the expectation to give a mean value since it would be an integral with value in the manifold: a new definition of mean is needed, which implies revisiting an important part of the theory of statistics.

Statistical analysis on manifolds is a relatively new domain at the confluent of several mathematical and application domains. Its goal is to study statistically geometric object living in differential manifolds. It is linked to the theory of statistical manifolds [START_REF] Amari | Differential-geometric methods in Statistics[END_REF][START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF], which aims at providing a Riemannian structure to the the space of parameters of statistical distribution. However, the targeted geometrical objects are usually different. Directional statistics [START_REF] Bingham | An antipodally symmetric distribution on the sphere[END_REF][START_REF] Jupp | A unified view of the theory of directional statistics, 1975-1988[END_REF][START_REF] Kent | The art of Statistical Science, chapter 10 : New Directions in Shape Analysis[END_REF][START_REF] Mardia | Directional statistics and shape analysis[END_REF] provide a first approach to statistics on manifold. As the manifolds considered here are spheres and projective spaces, the tools developed were mostly extrinsic, i.e. relying on the embedding of the manifold in the ambient Euclidean space. More complex objects are obtained when we consider the "shape" of a set of k points, i.e. what remains invariant under the action of a given group of transformation (usually rigid body ones or similarities). The statistics on these shape spaces [START_REF] Kendall | A survey of the statistical theory of shape (with discussion)[END_REF][START_REF] Dryden | Theoretical and distributional aspects of shape analysis[END_REF][START_REF] Le | The Riemannian structure of euclidean shape space: a novel environment for statistics[END_REF][START_REF] Small | The Statistical Theory of Shapes[END_REF] raised the need for intrinsic tools. However, the link between the tools developed in these works, the metric used and the space structure was not always very clear.

Another mathematical approach was provided by the study of stochastic processes on Lie groups. For instance, [START_REF] Grenander | Probabilities on Algebraic Structures[END_REF] derived central limit theorems on different families of groups and semi-groups with specific algebraic properties. Since then, several authors in the area of stochastic differential 1 geometry and stochastic calculus on manifolds proposed results related to mean values [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF][START_REF] Kendall | Convexity and the hemisphere[END_REF][START_REF] Emery | Sur le barycentre d'une probabilité dans une variété[END_REF][START_REF] Arnaudon | Barycentres convexes et approximations des martingales continues dans les variétés[END_REF][START_REF] Picard | Barycentres et martingales sur une variété[END_REF][START_REF] Darling | Martingales on non-compact manifolds: maximal inequalities and prescribed limits[END_REF]. On the applied mathematics and computer science side, people get interested in computing and optimizing in specific manifolds, like rotations and rigid body transformations [START_REF] Pennec | Feature-based registration of medical images: Estimation and validation of the pose accuracy[END_REF][START_REF] Grenander | Hilbert-schmidt lower bounds for estimators on matrix Lie groups for atr[END_REF][START_REF] Pennec | Computing the mean of geometric features -application to the mean rotation[END_REF][START_REF] Gramkow | On averaging rotations[END_REF][START_REF] Moakher | Means and averaging in the group of rotations[END_REF], Stiefel and Grassmann manifolds [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF], etc.

Over the last years, several groups attempted to federate some of the above approaches in a general statistical framework, with different objectives in mind. For instance, [START_REF] Hendricks | A Cramer-Rao type lower bound for estimators with values in a manifold[END_REF] and [START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF] aimed at characterizing the performances of statistical parametric estimators, like the bias and the mean square error. [START_REF] Hendricks | A Cramer-Rao type lower bound for estimators with values in a manifold[END_REF] considered extrinsic statistics, based on the Euclidean distance of the embedding space, while [START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF] considered the intrinsic Riemannian distance, and refined the Cramer-Rao lower bound using bounds on the sectional curvature of the manifold. In [START_REF] Bhattacharya | Nonparametric estimation of location and dispersion on Riemannian manifolds[END_REF][START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds, I[END_REF], the authors focused on the asymptotic consistency properties of the extrinsic and intrinsic means and variances for large sample sizes, and were able to propose a central limit theorem for flat manifolds. Here, in view of computer vision and medical image analysis applications, our concern is quite different: we aim at developing computational tools that can consistently deal with geometric features, or that provide at least good approximations. As we often have few measurements, we are interested in small sample sizes rather than large one, and we prefer to obtain approximations rather than bounds on the quality of the estimation. Thus, one of our special interest is to develop Taylor expansions with respect to the variance, in order to evaluate the quality of the computations with respect to the curvature of the manifold. In all cases, the chosen framework is the one of geodesically complete Riemannian manifolds, which appears to be powerful enough to support an interesting theory. To ensure a maximal consistency of the theory, we rely in this paper only on intrinsic properties of the Riemannian manifold, thus excluding methods based on the embedding of the manifold in an ambient Euclidean space.

We review in Section 2 some basic notions of differential and Riemannian geometry that will be needed afterward. This synthesis was inspired from [39, chap. 9], [START_REF] Klingenberg | Riemannian Geometry[END_REF][START_REF] Carmo | Riemannian Geometry[END_REF][START_REF] Gallot | Riemannian Geometry[END_REF], and the reader can refer to these books to find more details. In the remaining of the paper, we consider that our Riemannian manifold is connected and geodesically complete. We first detail in Section 3 the measure induced by the Riemannian metric on the manifold, which allows to define probability density functions, in particular the uniform one. Then, we turn in Section 4 to the expectation of a random point. We provide a quite comprehensive survey of the definitions that have been proposed. Among them, we focus on the Karcher and Fréchet means, defined as the set of points minimizing locally or globally the variance (the expected Riemannian distance). We provide a new proof of the barycentric characterization theorem and an original Gauss-Newton gradient descent algorithm to practically compute the mean. Once the mean value is determined, one can easily define the covariance matrix of a random point (and possibly higher order moments) using the exponential map at the mean point (Section 4). To generalize the Gaussian distribution, we propose in Section 6 a new family of distributions based on a maximum entropy approach. Under some reasonable hypotheses, we show that is amounts to take a truncated Gaussian distribution in the exponential map at the mean point. We illustrate the properties of this pdf family on the circle, and provide computationally tractable approximations for concentrated distributions. Last but not least, we investigate in Section 7 the generalization of the Mahalanobis distance and the χ 2 law. A careful analysis shows that, with our definition of the generalized Gaussian, the χ 2 law remains independent of the variance and of the manifold curvature, up to the order 3. This demonstrate that the whole framework is computationally sound and particularly consistent.

Differential geometry background 2.1 Riemannian metric, distance and geodesics

In the geometric framework, one specifies the structure of a manifold M by a Riemannian metric. This is a continuous collection of dot products . | . x on the tangent space T x M at each point x of the manifold. A local coordinate system x = (x 1 , . . . x n ) induces a basis ∂ ∂x = (∂ 1 , . . . ∂ n ) of the tangent spaces (∂ i is a shorter notation for ∂/∂x i ). Thus, we can express the metric in this basis by a symmetric positive definite matrix G(x) = [g ij (x)] where each element is given by the dot product of the tangent vector to the coordinate curves: g ij (x) = ∂ i | ∂ j . This matrix is called the local representation of the Riemannian metric in the chart x and the dot products of two vectors v and

w in T x M is now v | w x = v T G(x) w. The matrix G(x) is called the local representation of the Riemannian metric in the chart x.
If we consider a curve γ(t) on the manifold, we can compute at each point its instantaneous speed vector γ(t) and its norm, the instantaneous speed. To compute the length of the curve, we can proceed as usual by integrating this value along the curve:

L b a (γ) = b a γ(t) dt = b a γ(t) | γ(t) γ(t) 1 2 dt
The Riemannian metric is the intrinsic way of measuring length on a manifold. The extrinsic way is to consider the manifold as embedded in a larger vector space E (think for instance to the sphere S 2 in R 3 ) and compute the length of a curve in M as for any curve in E. In this case, the corresponding Riemannian metric is the restriction of the dot product of E onto the tangent space at each point of the manifold. By Whitney's theorem, there always exists such an embedding for a large enough vector space E (dim(E) ≤ 2dim(M) + 1).

To obtain a distance between two points of a connected Riemannian manifold, we simply have to take the minimum length among the smooth curves joining these points:

dist(x, y) = min γ L(γ) with γ(0) = x and γ(1) = y (1)
The curves realizing this minimum for any two points of the manifold are called geodesics 1) be the inverse of the metric matrix (in a given coordinate system x) and Γ i jk = 1 2 g im (∂ k g mj + ∂ j g mk -∂ m g jk ) the Christoffel symbols (using Einstein summation convention that implicit sum upon each index that appear up and down in the formula). The calculus of variations shows the geodesics are the curves satisfying the following second order differential system (in the chart

1 . Let [g ij ] = [g ij ] (-
x = (x 1 , . . . x n )): γi + Γ i jk γj γk = 0
The manifold is said to be geodesically complete if the definition domain of all geodesics can be extended to R. This means that the manifold has no boundary nor any singular point that we can reach in a finite time (for instance, R n -{0} with the usual metric is not geodesically complete, but R n or S n are). As an important consequence, the Hopf-Rinow-De Rham theorem state that such a manifold is complete for the induced distance (equation 1), and that there always exist at least one minimizing geodesic between any two points of the manifold (i.e. which length is the distance between the two points). From now on, we will assume that the manifold is geodesically complete.

Exponential map and cut locus

From the theory of second order differential equations, we know that there exists one and only one geodesic γ (x,∂v) going through the point x ∈ M at t = 0 with tangent vector ∂ v ∈ T x M. This geodesic is theoretically defined in a sufficiently small interval around zero but since we the manifold is geodesically complete, its definition domain can be extended to R. Thus, the point γ (x,∂v) (t) is defined for all vector ∂ v ∈ T x M and all parameter t. The exponential map maps each vector ∂ v to the point of the manifold reached in a unit time:

exp x : T x M -→ M ∂ v -→ exp x (∂ v ) = γ (x,∂v) (1)
This function realizes a local diffeomorphism from a sufficiently small neighborhood of 0 in T x M into a neighborhood of the point x ∈ M. We denote by log x = exp (-1) x the inverse map or simply -→ xy = log x (y). In this chart, the geodesics going through x are the represented by the lines going through the origin: log x γ (x, -→ xy) (t) = t -→ xy. Moreover, the distance with respect to the development point x is preserved:

dist(x, y) = -→ xy = ( -→ xy | -→ xy x ) 1/2
Thus, the exponential chart at x can be seen as the development of the manifold in the tangent space at a given point along the geodesics. This is also called a normal coordinate system if it is provided with an orthonormal basis. At the origin of such a chart, the metric reduces to the identity matrix and the Christoffel symbols vanish. Now, it is natural to search for the maximal domain where the exponential map is a diffeomorphism. If we follow a geodesic γ (x,∂v) (t) = exp x (t ∂ v ) from t = 0 to infinity, it is either always minimizing all along or it is minimizing up to a time t 0 < ∞ and not any more after (thanks to the geodesic completeness). In this last case, the point z = γ (x,∂v) (t 0 ) is called a cut point and the corresponding tangent vector t 0 ∂ v a tangential cut point. The set of of all cut points of all geodesics starting from x is the cut locus C(x) ∈ M and the set of corresponding vectors the tangential cut locus C(x) ∈ T x M. Thus, we have C(x) = exp x (C(x)), and the maximal definition domain for the exponential chart is the domain D(x) containing 0 and delimited by the tangential cut locus.

It is easy to see that this domain is connected and star-shaped with respect to the origin of T x M. Its image by the exponential map covers all the manifold except the cut locus and the segment [0, -→ xy] is transformed into the unique minimizing geodesic from x to y. Hence, the exponential chart at x is a chart centered at x with a connected and star-shaped definition domain that covers all the manifold except the cut locus C(x):

D(x) ∈ R n ←→ M -C(x) -→ xy = log x (y) ←→ y = exp x ( -→ xy)
From a computational point of view, it is often interesting to extend this representation to include the tangential cut locus. However, we have to take care of the multiple representations: points in the cut locus where several minimizing geodesics meet are represented by several points on the tangential cut locus as the geodesics are starting with different tangent vectors (e.g. rotation of π around the axis ±n for 3D rotations, antipodal point on the sphere). This multiplicity problem cannot be avoided as the set of such points is dense in the cut locus. The size of this definition domain is quantified by the injectivity radius i(M, x) = dist(x, C(x)), which is the maximal radius of centered balls in T x M on which the exponential map is one-to-one. The injectivity radius of the manifold i(M) is the infimum of the injectivity over the manifold. It may be zero, in which case the manifold somehow tends towards a singularity (think e.g. to the surface z = 1/ x 2 + y 2 as a sub-manifold of R 3 ).

Example: On the sphere S n (center 0 and radius 1) with the canonical Riemannian metric (induced by the ambient Euclidean space R n+1 ), the geodesics are the great circles and the cut locus of a points x is its antipodal point x = -x. The exponential chart is obtained by rolling the sphere onto its tangent space so that the great circles going through x become lines. The maximal definition domain is thus the open ball D = B n (π). On its boundary ∂D = C = S n-1 (π), all the points represent x. For the real projective space P n (obtained by identification of antipodal points of the sphere S n ), the geodesics are still the great circles, but the cut locus of the point {x, -x} is now the equator of the two points where antipodal points are still identified (thus the cut locus is P n-1 ). The definition domain of the exponential chart is the open ball D = B n ( π 2 ), and the tangential cut locus is the sphere ∂D = S n-1 ( π 2 ) where antipodal points are identified.

Taylor expansion of a real function

Let f be a smooth function from M to R (an observable). Its Gradient grad f at point x is the linear form on T x M corresponding to the directional derivatives ∂ v :

∀ v ∈ T x M grad f (v) = ∂ v f
Thanks to the dot product, we can identify the linear form dω in the tangent space T x M with the vector ω if dω(v) = ω | v x for all vector v ∈ T x M. All these notions can be extended to the whole manifold using vector fields: in a local chart x, we have

∂ v f = ∂f (x) ∂x v and ω | v x = ω T G(x) v.
Thus, the expression of the gradient in a chart is:

grad f = G (-1) (x) ∂f ∂x T = g ij ∂ j f
This definition corresponds to the classical gradient in R n even in the case of a non orthonormal basis. The second covariant derivative (the Hessian) is a little bit more complex and makes use of the connection ∇. We just give here its expression in a local coordinate system:

Hess f = ∇df = (∂ ij f -Γ k ij ∂ k f ) dx i dx j
Let now f x be the expression of f in a normal coordinate system at x. Its Taylor expansion around the origin is:

f x (v) = f x (0) + J fx v + 1 2 v T H fx v + O( v 3 )
where

J fx = [∂ i f ] and H fx = [∂ ij f ].
Since we are in a normal coordinate system, we have f x (v) = f (exp x (v)). Moreover, the metric at the origin reduces to the identity: J fx = grad f T , and the Christoffel symbols vanish so that the matrix of second derivatives H fx corresponds to the Hessian Hess f . Thus, The Taylor expansion can be written in any coordinate system:

f (exp x (v)) = f (x) + grad f (v) + 1 2 Hess f (v, v) + O( v 3 ) (2)
3 Random points on a Riemannian Manifold

In this paper, we are interested in measurements of elements of a Riemannian manifold that depend on the outcome of a random experiment. Particular cases are given by random transformation and random feature for the particular case of transformation groups and homogeneous manifolds.

Definition 1 (Random point on a Riemannian Manifold) Let (Ω, B(Ω), Pr) be a probability space, B(Ω) being the Borel σ-algebra of Ω (i.e. the smallest σ-algebra containing all the open subsets of Ω) and Pr a measure on that σ-algebra such that Pr(Ω) = 1. A random point in the Riemannian manifold M is a Borel measurable function x = x(ω) from Ω to M.

As in the real or vectorial case, we can now make abstraction of the original space Ω and directly work with the induced probability measure on M.

Riemannian measure or volume form

In a vector space with basis A = (a 1 , . . . a n ), the local representation of the metric is given by G = A T A where A = [a 1 , . . . a n ] is the matrix of coordinates change from A to an orthonormal basis. Similarly, the measure (or the infinitesimal volume element) is given by the volume of the parallelepipedon spanned by the basis vectors:

dV = | det(A)| dx = | det(G)| dx.
Assuming now a Riemannian manifold M, we can see that the Riemannian metric G(x) induces an infinitesimal volume element on each tangent space, and thus a measure on the manifold:

dM(x) = |G(x)| dx
One can show that the cut locus has a null measure. This means that we can integrate indifferently in M or in any exponential chart. If f is an integrable function of the manifold and

f x ( -→ xy) = f (exp x ( -→ xy)
) is its image in the exponential chart at x, we have: A simple example of a pdf is the uniform pdf in a bounded set X :

M f (x) dM = D(x) f x ( z) G x ( z) d z

Probability density function

p X (y) = 1 X dM 1 X (y) = 1 X (y) Vol(X )
One must be careful that this pdf is uniform with respect to the measure dM and is not uniform for another measure on the manifold. This problem is the basis of the Bertrand paradox for geometrical probabilities [START_REF] Poincaré | Calcul des probabilités[END_REF][START_REF] Kendall | Geometrical probability. Number 10 in Griffin's statistical monographs and courses[END_REF][START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF] and raise the problem of the measure to choose on the manifold. In our case, the measure is induced by the Riemannian metric but the problem is only lifted: which Riemannian metric do we have to choose ? For transformation groups and homogeneous manifolds, we showed in [START_REF] Pennec | L'incertitude dans les problèmes de reconnaissance et de recalage -Applications en imagerie médicale et biologie moléculaire[END_REF] that an invariant metric is a good geometric choice, even if such a metric does not always exist for homogeneous manifolds or if it leads in general to a partial consistency only between the geometric and statistical operations in non compact transformation groups [START_REF] Pennec | Probabilities and Statistics on Riemannian Manifolds: A Geometric approach[END_REF].

Expression of the density in a chart

Let x be a random point of pdf p x . If x = π(x) is a chart of the manifold defined almost everywhere, we obtain a random vector x = π(x) which pdf ρ x is defined with respect to the Lebesgue measure dx in R n instead of dM in M. Using the expression of the Riemannian measure, the two pdfs are related by

ρ x (y) = p x (y) |G(y)|
We shall note that the density ρ x depends on the chart used whereas the pdf p x is intrinsic to the manifold.

Expectation of an observable

Let ϕ(x) be a Borelian real valued function defined on M and x a random point of pdf p x . Then, ϕ(x) is a real random variable and we can compute its expectation:

E [ ϕ(x) ] = E x [ ϕ ] = M ϕ(y) p x (y) dM(y)
This notion of expectation corresponds to the one we defined on real random variables and vectors. However, we cannot directly extend it to define the mean value of the distribution since we have no way to generalize this integral in R into an integral with value in the manifold.

Expectation or Mean of a Random point

We focus in this section to the notion of central value of a distribution. We will preferably use the denomination mean value or mean point than expected point to stress the difference between this notion and the expectation of a real function.

Fréchet expectation or mean value

Let x be a random vector of R n . Fréchet observed in [START_REF] Fréchet | L'intégrale abstraite d'une fonction abstraite d'une variable abstraite et son application à la moyenne d'un élément aléatoire de nature quelconque[END_REF][START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF] that the variance σ 2 x (y) = E dist(x, y) 2 is minimized for the mean vector x = E [ x ]. The major point for the generalization is that the expectation of a real valued function is well defined for our connected and geodesically complete Riemannian manifold M.

Definition 3 Variance of a random point

Let x be a random point of pdf p x . The variance σ 2

x (y) is the expectation of the squared distance between the random point and the fixed point y:

σ 2 x (y) = E dist(y, x) 2 = M dist(y, z) 2 p x (z) dM(z) (3)
Definition 4 Fréchet expectation of a random point Let x be a random point. If the variance σ 2 x (y) is finite for all point y ∈ M (which is in particular true for a density with a compact support), every point x minimizing this variance is called an expected or mean point. Thus, the set of the mean points is:

E [ x ] = arg min y∈M E dist(y, x) 2
If there exists a least one mean point x, we call variance the minimal value σ 2

x = σ 2 x (x) and standard deviation its square-root.

Similarly, one defines the empirical or discrete mean point of a set of measures x 1 , . . .

x n : E [ {x i } ] = arg min y∈M E { dist(y, x i ) 2 } = arg min y∈M 1 n i dist(y, x i ) 2
If there exists a least a mean point x, one calls empirical variance the minimal value s 2 = 1 n i dist(x, x i ) 2 and empirical standard deviation or RMS (for Root Mean Square) its squareroot.

Following the same principle, one can define other types of central values. The mean deviation at order α is

σ x,α (y) = (E [ dist(y, x) α ]) 1/α = M dist(y, z) α p x (z) dM(z) 1/α
If this function is bounded on M, one call central point at order α every point xα minimizing it. For instance, the modes are obtained for α = 0. Exactly like in a vector space, they are the points where the density is maximal on the manifold (which is generally not a maximum for the density on the charts). The median point is obtained for α = 1. For α → ∞, we obtain the "barycenter" of the distribution support (which has to be compact).

The definition of these central values can be extended to the discrete case easily, except perhaps for the modes and for α → ∞. We note that the Fréchet expectation is defined for all metric space and not only for Riemannian manifolds.

Existence and uniqueness: Riemannian center of mass

As our mean point is the result of a minimization, its existence is not ensured (the global minimum could be unreachable) and anyway the result is a set and no longer a single element. This has to be compared with some central values in vector spaces, for instance the modes. However, the Fréchet expectation does not define all the modes even in vector spaces: one only keep the modes of maximal intensity.

To get rid of this constraint, Karcher [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] proposed to consider the local minima of the variance σ 2

x (y) (equation 3) instead of the global ones. We call these new set of means Riemannian centers of mass. As global minima are local minima, the Fréchet expected points are a subset of the Riemannian centers of mass. However, the use of local minima allows to characterize the Riemannian centers of mass using only local derivatives of order two. Using this extended definition, Karcher [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] and Kendall [START_REF] Kendall | Probability, convexity, and harmonic maps with small image I: uniqueness and fine existence[END_REF] established conditions on the manifold and the distribution to ensure the existence and uniqueness of the mean. We just recall here the results without the proofs.

Definition 5 (Regular geodesic balls) The ball B(x, r) = {y ∈ M/ dist(x, y) < r} is said geodesic if it does not meet the cut locus of its center. This means that there exists a unique minimizing geodesic from the center to any point of a geodesic ball. The ball is said regular if its radius verifies 2r √ κ < π, where κ is the maximum of the Riemannian curvature in this ball.

For instance, on the sphere S 2 with radius one, the curvature is constant and equal to 1. A geodesic ball is regular if r < π/2. Such a ball can almost cover an hemisphere, but not the equator. In a Riemannian manifold with non positive curvature, a regular geodesic ball can cover the whole manifold (according to the Cartan-Hadamard theorem, such a manifold is diffeomorphic to R n if it is simply connected and complete).

Theorem 1 (Existence and uniqueness of the Riemannian center of mass) Let x be a random point of pdf p x .

• Kendall [START_REF] Kendall | Probability, convexity, and harmonic maps with small image I: uniqueness and fine existence[END_REF] If the support of p x is included in a regular geodesic ball B(y, r), then there exists one and only one Riemannian center of mass x on this ball.

• Karcher [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] If the support of p x is included in a geodesic ball B(y, r) and if the ball of double radius B(y, 2 r) is still geodesic and regular, then the variance σ 2 x (z) is a convex function of z and has only one critical point on B(y, r), necessarily the Riemannian center of mass.

These conditions ensure a correct behavior of the mean for sufficiently localized distributions. However, they are quite restrictive as they only address pdfs with a compact support. Kendall's existence and uniqueness theorem was extended by [START_REF] Darling | Martingales on non-compact manifolds: maximal inequalities and prescribed limits[END_REF] to distributions with non-compact support in manifolds with Ψ-convexity. This notion, already introduced by Kendall in his original proof, is here extended to the whole manifold. Unfortunately, this type of argument can only be used for a restricted class of manifolds as a non-compact connected and geodesically complete Ψ-convex manifold is diffeomorphic to R m . It remains that this extension of the theorem applies to the important class of Hadamard manifolds (i.e. simply connected, complete and with non-positive sectional curvature), whose curvature is bounded from below.

Other possible definitions of the mean points

The Riemannian center of mass is perfectly adapted for our purpose, thanks to its good properties for optimization (see Section 4.6 below). However, there are other works proposing different ways to generalize the notion of mean value or barycenter of a distribution in a manifold. We review them for the sake of completeness and for their mathematical interest, but they do not seem to be practically applicable.

Doss [START_REF] Doss | Sur la moyenne d'un élément aléatoire dans un espace distancié[END_REF] used another property of the expectation as the starting point for the generalization: if x is a real random variable, the only real number x verifying:

∀y ∈ R |y -x| ≤ E [ |x -x| ]
is the mean value E [ x ]. Thus, in a metric space, the mean according to Doss is defined as the set of points x ∈ M verifying:

∀y ∈ M dist(y, x) ≤ E [ dist(x, x) ]
Herer shows in [START_REF] Herer | Espérance mathématique au sens de Doss d'une variable aléatoire à valeur dans un espace métrique[END_REF][START_REF] Herer | Espérance mathématique d'une variable aléatoire à valeur dans un espace métrique à courbure négative[END_REF] that this definition includes the classical expectation in a Banach space (with possibly other points) and develop on this basis a conditional expectation.

A similar definition that uses convex functions on the manifold instead of metric properties proposed by Emery [START_REF] Emery | Sur le barycentre d'une probabilité dans une variété[END_REF] and Arnaudon [START_REF] Arnaudon | Espérances conditionnelles et C-martingales dans les variétés[END_REF][START_REF] Arnaudon | Barycentres convexes et approximations des martingales continues dans les variétés[END_REF]. A function from M to R is convex if its restriction to all geodesic is convex (considered as a function from R to R). The convex barycenter of a random point x with density p x is the set B(x) of points y ∈ M such that α(y) ≤ E [ α(x) ] holds for every real bounded and convex function α on a neighborhood of the support of p x .

This definition seems to be of little interest in our case since for compact manifolds, such as the sphere or SO 3 (the manifold of 3D rotations), the geodesics are closed and the only convex functions on the manifold are the constant ones. Thus, every random point for which the support of the distribution is the whole manifold has the whole manifold as convex barycenter.

However, in the case where the support of the distribution is included in a strongly convex open set2 U, Emery [START_REF] Emery | Sur le barycentre d'une probabilité dans une variété[END_REF] showed that the exponential barycenters, defined as the critical points of the variance σ 2

x (y) are subset of the convex barycenter B(x). Local and global minima being particular critical points, the exponential barycenters include the Riemannian centers of mass that include themselves the Fréchet means.

Picard [START_REF] Picard | Barycentres et martingales sur une variété[END_REF] realized a good synthesis of most of these notions of mean value and show that the definition of a "barycenter" (i.e. a mean value) is linked to a connector, which determines itself a connection, and thus possibly a metric. An interesting property brought by this formulation is that the distance between two barycenters (with different definitions) is of the order of O(σ x ). Thus, for sufficiently centered random points, all these values are close.

Characterizing a Riemannian center of mass

To characterize a local minimum of a twice differentiable function, we just have to require a null gradient and a negative definite Hessian matrix. The problem with the variance function σ 2 (y) is that the integration domain (namely M\C(y)) depends on the derivation point y. Thus we cannot just use the Lebesgue theorem to differentiate under the sum, unless there is no cut locus, or the distribution has a sufficiently small compact support, which is the property used by Karcher Kendall and Emery for the previous existence, uniqueness and equivalence theorems. We were able to generalize in appendix A a differentiability proof of Pr Maillot [START_REF] Maillot | Différentielle de la variance et centrage de la plaque de coupure sur une variété riemannienne compacte[END_REF] originally designed for the uniform distribution on compact manifolds. The theorem we obtain is the following: Theorem 2 Gradient of the variance function Let P be a probability on the Riemannian manifold M. The variance σ 2 (y) = M dist(y, z) 2 dP (z) is differentiable at any point y ∈ M where it is finite and where the cut locus C(y) has a null probability measure:

P (C(y)) = C(y) dP (z) = 0 and σ 2 (y) = M dist(y, z) 2 dP (z) < ∞
At such a point, it has the following gradient:

(grad σ 2 )(y) = -2 M/C(y) -→ yz dP (z) = -2 E -→ yx
Now, we know that the variance is continuous but may not be differentiable at the points where the cut locus has a non-zero probability measure. At these points, the variance can have an extremum (think for instance to x in vector spaces). Thus, the extrema of σ 2 are characterized by (grad σ 2 )(y) = 0 if this is defined or P (C(y)) > 0.

Corollary 1 (Characterization of Riemannian centers of mass)

Assume that the random point x has a finite variance everywhere and let A be the set of points where the cut locus has a non-zero probability measure. A necessary condition for x to be a Riemannian center of mass is

E -→ xx = 0 if x ∈ A, or x ∈ A.
For discrete or empirical means, the characterization is the same but we can write explicitly the set

A = ∪ i C(x i ).
If the manifold does not have a cut locus (for instance in Hadamard manifolds), we have no differentiation problem. One could think of going one step further and computing the Hessian matrix. Indeed, we have in the vector case: Hess σ 2

x (y) = -2 Id everywhere, which proves that any extremum of the variance is a minimum. In Riemannian manifolds, one has to be more careful because the Hessian is modified by the curvature of the manifold. One solution is to compute the Hessian matrix using the theory of Jacobi fields, and then take estimates of its eigenvalues based on bounds on the curvature. This is essentially the idea exploited in [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] to show the uniqueness of the mean is small enough geodesic balls, and by [START_REF] Kendall | The propeller: a counterexample to a conjectured criterion for the existence of certain harmonic functions[END_REF] to exhibit an example of a manifold without cut-locus that is strongly convex (i.e. there is one and only one minimizing geodesic joining any two points), but that support finite mass measures that have non-unique Riemannian centers of mass. Thus, the absence of a cut locus is not enough: one should also have some constraints on the curvature of the manifold. In order to remain simple, we stick in this paper to the existence and uniqueness theorem provided by [START_REF] Darling | Martingales on non-compact manifolds: maximal inequalities and prescribed limits[END_REF] for simply connected and complete manifolds whose curvature is non-positive (i.e. Hadamard) and bounded from below.

Corollary 2 (Characterization of the Fréchet mean for Hadamard manifolds with a curvature bounded from below). Assume that the random point x has a finite variance. Then, there is one and only one Riemannian center of mass characterized by E -→ xx = 0. For discrete or empirical means, the characterization is similar.

Results similar to Theorem 2 and the above corollaries have been derived independently. [START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF] defined the mean values in manifolds as the exponential barycenters. To relate them with the Riemannian centers of mass, they determined the gradient of the variance. However, they only investigate the relation between the two notions when the probability is dominated by the Riemannian measure, which excludes explicitly point-mass distributions. In [START_REF] Bhattacharya | Nonparametric estimation of location and dispersion on Riemannian manifolds[END_REF][START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds, I[END_REF], the gradient of the variance is also determined and the existence of the mean is established for simply connected Riemannian manifolds with non-positive curvature.

Basically, the characterization of the Riemannian center of mass is the same as in Euclidean spaces if the curvature of manifold is non-positive (and bounded from below), in which case there is no cut-locus (we assumed that the manifold was complete and simply connected). If the sectional curvature becomes positive, a cut locus may appear, and a non-zero probability on this cut-locus induces some discontinuities in the first derivative of the variance. This corresponds to something like a Dirac measure on the second order derivative, which is an additional difficulty to compute the Hessian matrix of the variance on these manifolds.

Example on the circle

The easiest example of this difficulty is probably a symmetric distribution on the circle. Let p = cos(θ) 2 /π be the probability density function of our random point θ on the circle. For a circle with the canonical metric, the exponential chart centered at α is -→ αθ = θ -α for θ ∈]α -π; α + π[, the distance being obviously dist(α, θ) = |α -θ| within this domain.

Let us first compute the mean points by computing explicitly the variance and its derivatives. The variance is:

σ 2 (α) = α+π α-π dist(α, θ) 2 p(θ) dθ = π -π γ 2 cos(γ + α) 2 π dγ = π 2 3 - 1 2 + cos(α) 2 .
Its derivative is rather easy to compute: grad σ 2 (α) = -2 cos(α) sin(α), and the second order derivative is H(α) = 4 sin(α) 2 -2. Solving for grad σ 2 (α) = 0, we get four critical points:

• α = 0 and α = ±π with H(0) = H(±π) = -2, • α = ±π/2 with H(±π/2) = +2.
Thus, there are two relative (and here absolute) minima:

E [ θ ] = {0 , ±π}.
Let us use now the general framework developed on Riemannian manifolds. According to theorem 2, the gradient of the variance is

grad σ 2 (α) = -2E -→ αθ = -2 α+π α-π -→ αθ dθ = -2 α+π α-π (θ -α) cos(θ) 2 π dθ = -2 cos(α) sin(α),
which is in accordance with our previous computations. Now, differentiating once again under the sum, we get:

α+π α-π ∂ 2 dist(α, θ) ∂α 2 p(θ) dθ = -2 α+π α-π ∂ -→ αθ ∂α p(θ) dθ = 2 α+π α-π p(θ) dθ = 2,
which is clearly different from our direct calculation. One way to see the problem is the following: the vector field -→ αθ is continuous and differentiable on the circle except at the cut locus of α (i.e. at θ = α ± π) where it has a jump of 2π. Thus, the second order derivative of the squared distance should be -2(-1 + 2πδ (α±π) (θ)), where δ is the Dirac distribution, and the integral becomes:

H(α) = -2 α+π α-π -1 + 2πδ (α±π) (θ) p(θ) dθ = 2 -4π p(α ± π) = 2 -4 cos(θ) 2
which is this time in accordance with the direct calculation.

A gradient descent algorithm to obtain the mean

Gradient descent is a usual technique to compute a minimum. Moreover, as we have a canonical way to go from the tangent space to the manifold thanks to the exponential map, this iterative algorithm seems to be perfectly adapted. In this section, we assume that the conditions of theorem (2) are fulfilled.

Let y be an estimation of the mean of the random point x and f (y) = σ 2 x (y) the variance. A practical gradient descent algorithm is to minimize the second order approximation of the cost function at the current point. According to the Taylor expansion of equation ( 2), the second order approximation of f and y is:

f (exp y (v)) = f (y) + grad f (v) + 1 2 Hess f (v, v)
This is a function of the vector v ∈ T y M. Assuming that Hess f is positive definite, this function is convex and has thus a minimum characterized by a null gradient. Let H f (v) denote the linear form verifying H f (v)(w) = Hess f (v, w) for all w and H (-1) f denote the inverse map. The minimum is characterized by

grad v f y = 0 = grad f + H f (v) ⇔ v = -H (-1) f (grad f )
We saw in the previous section that grad f = -2 E -→ yx . Neglecting the "cut locus term" in the Hessian matrix gives us a perfect positive definite matrix Hess f 2 Id. Thus, the gradient descent algorithm is

y t+1 = exp yt E -→ y t x
This gradient descent algorithm can be seen as the discretization of the ordinary differential equation ẏ(t) = E ---→ y(t)x . Other discretization scheme are possible [START_REF] Hairer | Geometric numerical integration : structure preserving algorithm for ordinary differential equations[END_REF], sometimes with convergence theorems [START_REF] Dedieu | Newton method on Riemannian manifolds: Covariant alpha-theory[END_REF].

In the case of the discrete or empirical mean, which is much more interesting from a statistical point of view, we have exactly the same algorithm, but with the empirical expectation:

y t+1 = exp yt 1 n n i=1 --→ y t x i
We note that in the case of a vector space, these two formula simplify to y t+1 = E [ x ] and y t+1 = 1 n i x i , which are the definition of the mean value and the barycenter. Moreover, the algorithm converges in a single step.

An important point for this algorithm is to determine a good starting point. In the case on a set of observations {x i }, one can choose at random one of the observations as the starting point. Another solution is to map to each point x i its mean distance with respect to other points (or the median distance to be robust) and choose as the starting point the minimizing point. From a computer science point of view, the complexity is k 2 (where k is the number of observations) but the method can be randomized efficiently [START_REF] Huber | Robust Statistics[END_REF][START_REF] Rousseeuw | Robust Regression and Outliers Detection[END_REF].

To verify the uniqueness of the solution, we can repeat the algorithm from several starting points (for instance all the observations x i ). If we know the Riemannian curvature of the manifold (for instance if it is constant or if there is an upper bound κ), we can use theorem (1, Section 4.2). We just have to verify that the maximum distance between the observations and the mean value we have found is sufficiently small so that all observations fits into a regular geodesic ball of radius:

r = max i dist(x, x i ) < π 2 √ κ

Covariance matrix

With the mean value, we have a dispersion value: the variance. To go one step further, we observe that the covariance matrix of a random vector x with respect to a point y is the directional dispersion of the "difference" vector -→ yx = x -y:

Cov x (y) = E -→ yx -→ yx T = R n ( -→ yx) ( -→ yx) T p x (x) dx
This definition is readily extendible to a complete Riemannian manifold using the random vector -→ yx in T y M and the Riemannian measure. In fact, we are usually interested in the covariance relative to the mean value: Definition 6 (Covariance) Let x be a random point and x ∈ E [ x ] a mean value that we assume to be unique to simplify the notations (otherwise we have to keep a reference to the mean value). We define the covariance Σ xx by the expression:

Σ xx = Cov x (x) = E -→ xx -→ xx T = D(x) ( -→ xx) ( -→ xx) T p x (x) dM(x)
The empirical covariance is defined in the same way using the discrete version of the expectation operator.

We observe that the covariance depends on the basis used for the exponential chart if we see it as a matrix, but it does not depend on it if we consider it as a bilinear form over the tangent plane.

The covariance is related to the variance just as in the vector case:

Tr(Σ xx ) = E Tr( -→ xx -→ xx T ) = E dist(x, x) 2 = σ 2 x
This formula is still valid relatively to any fixed point: Tr(Cov x (y)) = σ 2 x (y). In fact, as soon as we have found a (or the) mean value and that the probability of its cut locus is null, everything appears to be similar to the case of a centered random vector by developing the manifold onto the tangent space at the mean value. Indeed, -→ xx is a random vector of pdf ρ x (y) = p x (y) |G(y)| with respect to the Lebesgue measure in the connected and star-shaped domain D(x) ⊂ T xM. We know that its expectation is E -→ xx = 0 and its covariance matrix is defined as usual. Thus, we could define higher order moments of the distribution by tensors on this tangent space, just as we have done for the covariance.

An information-based generalization of the Normal distribution

In this section, we turn to the generalization of the Gaussian distribution to Riemannian manifolds. Several generalizations have already be proposed so far. In the stochastic calculus community, one usually consider the heat kernel p(x, y, t), which is the transition density of the Brownian motion [START_REF] Grenander | Probabilities on Algebraic Structures[END_REF][START_REF] Emery | Stochastic Calculus in Manifolds[END_REF][START_REF] Grigor | Heat kernels on weighted manifolds and applications[END_REF]. This is the smallest positive fundamental solution to the heat equation ∂f ∂t -∆f =0, where ∆ is the Laplace-Beltrami operator (i.e. the standard Laplacian with corrections for the Riemannian metric). On compact manifolds, an explicit basis of the heat kernel is given by the spectrum of the manifold-Laplacian (eigenvalues λ i with associated eigenfunctions f i solutions of ∆f = λf ). The practical problem is that the computation of this spectrum is impossible but in very few cases [START_REF] Gallot | Riemannian Geometry[END_REF].

To obtain tractable formulas, several distributions have been proposed in directional statistics [START_REF] Bingham | An antipodally symmetric distribution on the sphere[END_REF][START_REF] Jupp | A unified view of the theory of directional statistics, 1975-1988[END_REF][START_REF] Kent | The art of Statistical Science, chapter 10 : New Directions in Shape Analysis[END_REF][START_REF] Mardia | Directional statistics and shape analysis[END_REF][START_REF] Mardia | Directional statistics[END_REF], in particular the wrapped Gaussian distributions. The basic idea is to take the image by the exponential of a Gaussian distribution on the tangent space centered at the mean value (see e.g. [START_REF] Mardia | Directional statistics[END_REF] for the circular and spherical case). It is easy to see that the wrapped Gaussian distribution tends towards the mass distribution if the variance goes to zero. In the circular case, one can also show that is tends toward the uniform distribution for a large variance. [START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF] extended this definition by considering non-centered Gaussian distributions on the tangent spaces of the manifold in order to tackle the asymptotic properties of estimators. In this case, the mean value is generally not any more simply linked to the Gaussian parameters. In view of a computational theory, the main problem is that the pdf of the wrapped distributions can only be expressed if there is a particularly simple geometrical shape of the cut-locus. For instance, considering an anisotropic covariance on the n-dimensional sphere leads to very complex calculations.

Thus, instead of keeping a Gaussian pdf in some tangent space, we propose in this paper a new variational approach based on global properties, consistent with the previous definitions of the mean and covariance. The property that we take for granted is the maximization of the entropy knowing the mean and the covariance matrix. For many aspects, this may not be the best generalization of the Gaussian. However, we demonstrate that it provides a family ranging from the point-mass distribution to the uniform measure (for compact manifolds) and that we can provide computationally tractable approximations for any manifold in case of small variances. In this section the symbols log and exp denote the standard logarithmic and exponential functions in R.

Entropy and uniform law

As we can integrate a real valued function, the extension of the entropy H [ x ] of a random point is straightforward:

H [ x ] = E [ -log(p x (x)) ] = - M log(p x (x)) p x (x) dM(x)
This definition is coherent with the measure inherited from our Riemannian metric since the pdf p U that maximizes the entropy when we only know that the measure is in a compact set U is the uniform density in this set:

p U (x) = 1 U (x) U dM(y)

Maximum entropy characterization

Now assume that we only know the mean (that we suppose to be unique) and the covariance of a random point: we denote it by x ∼ (x, Σ). If we need to assume a pdf for that random point, it seems reasonable to choose the one which is the least informative, while fitting the mean and the covariance. The pdf is maximizing in this case the conditional entropy

H [ x | x ∈ E [ x ] , Σ xx = Σ ]
In standard multivariate statistics, this maximum entropy principle is one characterization of the Gaussian distributions [62, p. 409]. In the Riemannian case, we can express all the constraints directly in the exponential chart at the mean value. Let ρ(y) = p(exp x(y)) be the density in the chart with respect to the induced Riemannian measure dM x(y) = |G x(y)| dy (we use here the Riemannian measure instead of the Lebesgue one to simplify equations below). The constraints are:

• the normalization: E [ 1 M ] = D(x) ρ(y) dM x(y) = 1 • a nul mean value: E -→ xx = D(x)
y ρ(y) dM x(y) = 0,

• and a fixed covariance Σ :

E -→ xx -→ xx T = D(x) y y T ρ(y) dM x(y) = Σ
To simplify the optimization, we won't consider any continuity or differentiability constraint on the cut locus C(x) (which would mean constraints on the border of the domain). This means that we can do the optimization in the exponential chart at the mean point as if we were in the open domain D(x) ∈ R n .

Using the convexity of the function -x log(x), one can show that the maximum entropy is attained by distributions of density ρ(y) = k exp -β T y -y T Γ y 2 , if there exists constants k, β and Γ such that our constrains are fulfilled [62, Theorem 13.2.1, p. 409]. Assuming that the definition domain D(x) is symmetric with respect to the origin, we find that β = 0 ensures a null mean. Substituting in the constraints gives the following equations.

Theorem 3 (Normal law)

We call Normal law on the manifold M the maximum entropy distribution knowing the mean value and covariance. Assuming no continuity nor differentiability constraint on the cut locus C(x) and a symmetric domain D(x), the pdf of the Normal law of mean x (the mean value) and concentration matrix Γ is given by:

N (x,Γ) (y) = k exp - -→ xy T Γ -→ xy 2
where the normalization constant and the covariance are related to the concentration matrix by:

k (-1) = M exp - -→ xy T Γ -→ xy 2 dM(y) and Σ = k M -→ xy -→ xy T exp - -→ xy T Γ -→ xy 2 dM(y)
From the concentration matrix, we can compute the covariance of the random point, at least numerically, but the reverse is more difficult.

The vectorial case

The integrals can be entirely computed, and we find k 1) . The Normal density is thus the usual Gaussian density:

(-1) = (2π) n 2 √ |Γ| and Γ = Σ (-
N (x,Γ) (x) = k exp - (x -x) T Γ (x -x) 2 = 1 (2 π) n/2 |Σ| exp - (x -x) T Σ (-1) (x -x) 2

Example on a simple manifold: the circle

The exponential chart for the circle of radius 1 with the canonical metric is the angle θ ∈ D =]-π; π[ with respect to the development point, and the measure is simply dθ. For a circle of radius r, the exponential chart becomes x = r θ. The domain is D =] -a; a[ (with a = π r) and the measure is dx = r dθ. Thus, the normalization factor of the Normal density is:

k (-1) = a -a exp - γ x 2 2 dx = 2π γ erf γ 2 a
where erf is the error function erf(x) = 2 √ π x 0 exp(-t 2 ) dt. The density is the truncated Gaussian

N (0,γ) (x) = k exp -γ x 2 2
. It is continuous but not differentiable on the cut locus π ≡ -π. The truncation introduces a bias in the relation between the variance and the concentration parameter:

σ 2 = a -a x 2 k exp - γ x 2 2 dx = 1 γ 1 -2 a k exp - γ a 2 2
It is interesting to have a look on limit properties: if the circle radius goes to infinity, the circle becomes the real line and we obtain the usual Gaussian with the relation σ 2 = 1/γ, as expected. Now, let us consider a circle with a fixed radius. As anticipated, the variance goes to zero and the density tends toward a point mass distribution (see figure 3) if the concentration γ goes to infinity. On the other hand, the variance cannot becomes infinite (as in the real case) when the concentration parameter γ goes to zero because the circle is compact: a Taylor expansion gives σ 2 = a 2 /3 + O(γ). Thus, the maximal variance on the circle is

σ 2 0 = lim γ→0 σ 2 = a 2 3 with the density N (0,0) (x) = 1 2 a
Interestingly, the Normal density of concentration 0 is the uniform density. In fact, this result can be generalized to all compact manifolds Figure 3: Variance σ 2 with respect to the concentration parameter γ on the circle of radius 1 and the real line. This variance tends toward σ 2 0 = π 2 /3 for the uniform distribution on the circle (γ = 0) whereas it tends to infinity for the uniform measure on R. For a strong concentration (γ > 1), the variance on the circle can be accurately approximated by σ 2 1/γ, as in the real case.

Small concentration matrix on compact manifolds: uniform distribution

Let n be the dimension of the manifold. Using bounds on the eigenvalues of Γ, we have:

Tr(Γ) -→ xy 2 /n ≤ -→ xy T Γ -→ xy ≤ Tr(Γ) -→ xy 2
. This means that we can bound the exponential by:

1 + +∞ k=1 (-1) k Tr(Γ) k -→ xy 2k k! 2 k ≤ exp(- -→ xy T Γ -→ xy/2) ≤ 1 + +∞ k=1 (-1) k Tr(Γ) k -→ xy 2k k! 2 k n k
As all the moments M -→ xy 2k dM are finite since the manifold is compact, one can conclude that

k (-1) = M exp(- -→ xy T Γ -→ xy/2) dM = V ol(M) + O(Tr(Γ))
It follows immediately that the limit of the normal distribution is the uniform one, and that the limit covariance is finite.

Theorem 4 Let M be a compact manifold. The limit of the normal distribution for small concentration matrices is the uniform density N (y) = 1/V ol(M) + O(Tr(Γ)) . Moreover, the covariance matrix tends towards a finite value:

Σ = 1 V ol(M) M -→ xy -→ xy T dM + O(Tr(Γ)) < +∞

Approximation for a large concentration matrix

If the pdf is sufficiently concentrated (a high concentration matrix Γ or a small covariance matrix Σ), then we can use a Taylor expansion of the metric in a normal coordinate system around the mean value to approximate the previous integrals and obtain a Taylor expansions of the normalization factor and the concentration matrix with respect to the covariance matrix.

The Taylor expansion of the metric is given by [63, p84]. We easily deduce the Taylor expansion of the measure around the origin (Ric is the Ricci (or scalar) curvature matrix in the considered normal coordinate system):

dM(y) = det(G(y)) dy = 1 - y T Ric y 6 + O( y 3 ) dy
Substituting this Taylor expansion in the integrals and manipulating the formulas (see appendix B) leads to the following theorem.

Theorem 5 (Approximate normal density)

In a complete Riemannian manifold, the normal density is

N (y) = k exp(- -→ xy T Γ -→ xy/2). Let r = i(M,
x) be the injectivity radius at the mean point (by convention r = +∞ if there is no cut-locus). Assuming a finite variance for any concentration matrix Γ, we have the following approximations of the normalization constant and concentration matrix for a covariance matrix Σ of small variance σ 2 = Tr(Σ):

k = 1 + O(σ 3 ) + ε σ r (2π) n det(Σ) and Γ = Σ (-1) - 1 3 Ric + O(σ) + ε σ r
Here, ε(x) is a function that is a O(x k ) for any positive k, with the convention that ε σ

+∞ = ε(0) = 0. More precisely, this is a function such that ∀k ∈ R + , lim 0+ x -k ε(x) = 0

Discussion

The maximum entropy approach to generalize the normal distribution to Riemannian manifolds is interesting since we obtain a whole family of densities going from the Dirac measure to the uniform distribution (or the uniform measure if the manifold is only locally compact). Unfortunately, this distribution is generally not differentiable at the cut locus, and often even not continuous. However, if the relation between the parameters and the moments of the distribution are not as simple as in the vector case (but can we expect something simpler in the general case of Riemannian manifolds?), the approximation for small covariances turns out to be rather simple. Thus, this approximate distribution can be handled quite easily for computational purposes. It is likely that similar approximations hold for wrapped Gaussian, but this remains to establish.

Mahalanobis distance and χ 2 law

The problem we are now tackling is to determine if an observation x could reasonably have come from a given probability distribution x on M or if it should be considered as an outlier. From a computational point of view, the pdf of the measurement process is often too rich an information to be estimated or handled. In practice, one often characterizes this pdf by its moments, and more particularly by the mean and the covariance. We denote it by x ∼ (x, Σ xx ). Based on these characteristics only, we have to decide if the observation x is compatible with this measurement process.

In the vector case, a well adapted tool is the Mahalanobis distance µ 2 = (x -x) T Σ (-1) xx (x -x), which measures the distance between the observation x and the mean value x according to the "metric" Σ (-1) xx . To test if the observation x is an outlier, the principle of the Mahalanobis D 2 test is to compute the tail probabilities of the resulting distribution (the so called p-value), i.e. the risk of error when saying that the observation does not come from the distribution. The distribution of x is usually assumed to be Gaussian, as this distribution minimizes the added information (i.e. minimizes the entropy) when we only know the mean and the covariance. In that case, we know that the Mahalanobis distance should be χ 2 n distributed if the observation is correct (n is the dimension of the vector space). If the probability of the current distance is too small (i.e. µ 2 is too large), the observation x can safely be considered as an outlier.

The definition of the Mahalanobis distance can be easily generalized to complete Riemannian manifolds with our tools. We note that it is well defined for any distribution of the random point and not only the normal one.

Definition 7 (Mahalanobis distance)

We call Mahalanobis distance between a random point x ∼ (x, Σ xx ) and a (deterministic) point y on the manifold the value

µ 2 x (y) = -→ xy T Σ (-1)
xx -→ xy.

Properties

Since µ 2

x is a function from M to R, µ 2 x (y) is a real random variable. The expectation of this random variable is well defined and turns out to be quite simple:

E µ 2 x (y) = M µ 2 x (z) p y (z) dM(z) = M -→ xz T Σ (-1) xx -→ xz p y (z) dM(z) = Tr Σ (-1) xx M -→ xz -→ xz T p y (z) dM(z) = Tr Σ (-1) xx Cov y (x)
The expectation of the Mahalanobis distance of a random point with itself is even simpler:

E µ 2 x (x) = Tr(Σ (-1) xx Σ xx ) = Tr( Id n ) = n Theorem 6
The expected Mahalanobis distance of a random point with itself is independent of the distribution and does only depend on the dimension of the manifold:

E µ 2 x (x) = n.
This identity can be used to verify with a posteriori observations that the covariance matrix has been correctly estimated. It can be compared with the expectation of the "normalized" squared distance, which is by definition: E dist(x, x) 2 /σ 2 x = 1.

A generalized χ 2 law

Assuming that the random point x ∼ (x, Σ xx ) is normal, we can go one step further and compute the probability that χ 2 = µ 2 x < α 2 (see appendix C). This generalization of the χ 2 law turns out to be still independent of the mean value and the covariance matrix of the random point (at least up to the order O(σ 3 )):

Theorem 7 (Approximate χ 2 law) With the same hypotheses as for the approximate normal law, the χ 2 probability is

Pr{χ 2 ≤ α 2 } = (2π) -n 2 x ≤α exp - x 2 2 dx + O(σ 3 ) + ε σ r
while the density is

p χ 2 (u) = 1 2 Γ n 2 u 2 n 2 -1 exp - u 2 + O(σ 3 ) + ε σ r
The χ 2 probability can be computed using the incomplete gamma function Pr{χ 2 ≤ α 2 } = P n 2 , α 2

2

(see for instance [START_REF] Press | Numerical Recipices in C[END_REF]). In practice, one often use this law to test if an observation x has been drawn from a random point x that we assume to be Gaussian: if the hypothesis is true, the value µ 2

x (x) will be less than α 2 with a probability γ = Pr{χ 2 ≤ α 2 }. Thus, one choose a confidence level γ (for instance 95% or 99%), then we find the value α(γ) such that γ = Pr{χ 2 ≤ α 2 } and accept the hypothesis if µ 2

x (x) ≤ α 2 .

Discussion

On a (geodesically complete) Riemannian manifold, it is easy to define probability density functions associated to random points, thanks to the availability of a metric. However, as soon as the expectation is concerned, we may only define the expectation of an observable (a real or vectorial function of the random point). Thus, the definition of a mean value for a random point is much more complex than for the vectorial case and it requires a distance-based variational formulation: the Riemannian center of mass basically minimize locally the variance. As the mean is now defined through a minimization procedure, its existence and uniqueness are not ensured any more (except for distributions with a sufficiently small compact support). In practice, one mean value almost always exists, and it is unique as soon as the distribution is sufficiently peaked. The properties of the mean are very similar to those of the modes (that can be defined as central values of order 0) in the vectorial case. We present here a new proof of the barycentric characterization theorem that is valid for distribution with any kind of support. The main difference with the vector case is that we have to ensure a null probability measure of the cut locus. To compute the mean value, we designed an original Gauss-Newton gradient descent algorithm that essentially alternates the computation of the barycenter in the exponential chart centered at the current estimation, and a re-centering step of the chart at the newly computed barycenter.

To define higher moments of the distribution, we used the exponential chart at the mean point (which may be seen as the development of the manifold onto its tangent space at this point along the geodesics): the random point is thus represented as a random vector with null mean in a star-shaped and symmetric domain. With this representation, there is no more problem to define the covariance matrix and potentially higher order moments. Based on this covariance matrix, we defined a Mahalanobis distance between a random and a deterministic point that basically weights the distance between the deterministic point and the mean point using the inverse of the covariance matrix. Interestingly, the expected Mahalanobis distance of a random point with itself is independent of the distribution and is equal to the dimension of the manifold, as in the vectorial case.

Like for the mean, we choose a variational approach to generalize the Normal law: we define it as the maximum entropy distribution knowing the mean and the covariance. Neglecting the cut-locus constraints, we show that this amounts to consider a truncated Gaussian distribution on the exponential chart centered at the mean point. However, the relation between the concentration matrix (the "metric" used in the exponential of the pdf) and the covariance matrix is slightly more complex that the simple inversion of the vectorial case, as it has to be corrected for the curvature of the manifold.

Last but not least, using the Mahalanobis distance of a Normally distributed random point, we can generalize the χ 2 law: we were able to show that is has the same density as in the vectorial case up to an order 3 in σ. This opens the way to the generalization of many other statistical tests, as we may expect similarly simple approximations for sufficiently centered distributions.

In this paper, we focused on purpose on the theoretical formulation of the statistical framework in geodesically complete Riemannian spaces, eluding applications examples. The interested reader will find practical applications in computer vision to compute the mean rotation [START_REF] Pennec | Computing the mean of geometric features -application to the mean rotation[END_REF][START_REF] Gramkow | On averaging rotations[END_REF] or for the generalization of matching algorithms to arbitrary geometric features [START_REF] Pennec | Toward a generic framework for recognition based on uncertain geometric features[END_REF]. In medical image analysis, selected applications cover the validation of the rigid registration accuracy [START_REF] Pennec | Feature-based registration of medical images: Estimation and validation of the pose accuracy[END_REF][START_REF] Roche | Rigid registration of 3D ultrasound with mr images: a new approach combining intensity and gradient information[END_REF][START_REF] Nicolau | Evaluation of a new 3D/2D registration criterion for liver radio-frequencies guided by augmented reality[END_REF], shape statistics [START_REF] Fletcher | Gaussian distributions on Lie groups and their application to statistical shape analysis[END_REF] and more recently tensor computing, either for processing and analyzing diffusion tensor images [START_REF] Batchelor | A rigorous framework for diffusion tensor calculus[END_REF][START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors[END_REF][START_REF] Ch | Statistics on multivariate normal distributions: A geometric approach and its application to diffusion tensor MRI[END_REF][START_REF] Pennec | A Riemannian framework for tensor computing[END_REF], or to model the brain variability [START_REF] Fillard | Extrapolation of sparse tensor fields: Application to the modeling of brain variability[END_REF]. One can even find applications in rock mechanics with the analysis of fracture geometry [START_REF] Rasouli | Application of Riemannian multivariate statistics to the analysis of rock fracture surface roughness[END_REF].

In the theory presented here, all definitions are derived from the Riemannian metric of the manifold. A natural question is how to chose this metric? Invariance properties requirements provide a partial answer for connected Lie groups and homogeneous manifolds [START_REF] Pennec | L'incertitude dans les problèmes de reconnaissance et de recalage -Applications en imagerie médicale et biologie moléculaire[END_REF][START_REF] Pennec | Probabilities and Statistics on Riemannian Manifolds: A Geometric approach[END_REF]. However, an invariant metric does not always exist on an homogeneous manifolds. Likewise, left and right invariant metric and generally different in non-compact Lie groups, so that we only have a partial consistency between the geometric and statistical operations. Another way to chose the metric could be to estimate it from empirical data.

More generally, we could design other definitions of the mean value using the notion of connector introduced in [START_REF] Picard | Barycentres et martingales sur une variété[END_REF]. This connector formalizes a relationship between the manifold and its tangent space at one point, exactly in the way we used the exponential map of the Riemannian metric. Thus, we could generalize easily the higher order moments and the other statistical operations we defined by replacing everywhere the exponential map with the connector. One important point is that these connectors can model extrinsic distances (like the Euclidean distance on unit vectors), and could lead to very efficient approximations of the mean value for sufficiently peaked distributions. For instance, the "barycenter / re-centering" algorithm we designed will most probably converge toward a first order approximation of the Riemannian mean if we use any chart that is consistent to the exponential chart at the first order (e.g. Euler's angle on re-centered 3D rotations). We believe that this research track may become one of the most productive from the practical point of view.

Hypotheses

Let P be a probability on the Riemannian manifold M. We assume that the cut locus C(y) of the derivation point y ∈ M has a null measure with respect to this probability (as it has with the Riemannian measure) and that the variance is finite at that point:

P (C(y)) = C(y) dP (z) = 0 and σ 2 (y) = M dist(y, z) 2 dP (z) < ∞

Goal and problem

Let now g(y) = M\C(y) -→ yz dP (z). As -→ yz = d(z, y) ≤ 1 + d(z, y) 2 , and using the null probability of the cut locus, we have:

g(y) ≤ M\C(y) -→ yz dP (z) ≤ M\C(y) 1 + d(z, y) 2 dP (z) = 1 + σ 2 (y) < ∞
Thus g(y) is well defined everywhere. Let h(y) = d(z, y) 2 = -→ yz 2 . For a fixed z ∈ C(y) (which is equivalent to y ∈ C(z)), we have (grad h)(y) = -2 -→ yz. Thus the proposition:

(grad σ 2 )(y) = -2 M/C(y) -→ yz dP (z)
corresponds to a derivation under the sum, but the usual conditions of the Lebesgue theorem are not fulfilled: the zero measure set C(y) varies with y. Thus, we have to come back to the original definition of the gradient.

Definition of the gradient

Let γ(t) be a curve with γ(0) = y and γ(0) = w. By definition, the function σ 2 : M → R is derivable if there exists a vector (grad σ 2 )(y) ∈ T y M (the gradient) such that:

∀w ∈ T y M (grad σ 2 )(y) | w = ∂ w σ 2 (y) = lim t→0 σ 2 (γ(t)) -σ 2 (y) t
Since tangent vectors are defined as equivalent classes, we can choose the geodesic curve γ(t) = exp y (t w). Using v = t w, the above condition can then be rewritten:

∀v ∈ T y M lim v →0 σ 2 (exp y (v)) -σ 2 (y) -(grad σ 2 )(y) | v v = 0
which can be rephrased as: for all η ∈ R + , there exists ε sufficiently small such that:

∀v ∈ T y M, v < ε σ 2 (exp y (v)) -σ 2 (y) -(grad σ 2 )(y) | v ≤ η v

General idea

Let ∆(z, v) be the integrated function (for z ∈ C(y)):

∆(z, v) = dist(exp y (v), z) 2 -dist(y, z) 2 + 2 -→ yz | v
and H(v) = M\C(y) ∆(z, v) dP (z) be the function to bound:

H(v) = σ 2 (exp y (v)) -σ 2 (y) -(grad σ 2 )(y) | v = M dist(exp y (v), z) 2 dP (z) - M dist(y, z) 2 dP (z) + 2 M\C(y) -→ yz | v dP (z)
The idea is to split this integral in two in order to bound ∆ on a small neighborhood W around the cut locus of y and to use the standard Lebesgue theorem to bound the integral of ∆ on M\W . Let us first reformulate the definition of W ε :

z ∈ W ε ⇔ ∃x ∈ B(y, ε) / z ∈ C(x) ⇔ ∃x ∈ B(y, ε) / x ∈ C(z) ⇔ C(z) ∩ B(y, ε) = ∅
Going to the limit, we have: z

∈ W 0 = lim ε→0 W ε ⇔ C(z) ∩ {y} = ∅ ⇔ z ∈ C(y)
. Thus, we have a continuous series of included and decreasing sets all containing C(y) and converging toward it. Now, let us prove that W ε is an open set.

Let U = {u = (x, v) ∈ M × T x M ; v = 1} be the unit tangent bundle of M and ρ : U → R+ = R + ∪ {+∞} be the cutting abscissa of the geodesic starting at x with the tangent vector v. Let now U = ρ (-1) (R + ) be the subset of the unit tangent bundle where ρ(u) < +∞. The function ρ being continuous on U , the subspace U = ρ (-1) ([0, +∞[) is open. Let π : U → M be the canonical projection along the fiber π((x, v)) = x. This is a continuous and open map. Let us denote by U x = π (-1) ({x}) the unit tangent bundle at point x and U x = π (-1) ({x}) ∩ U its subset that lead to a cutting point of x.

Consider u = (x, v) ∈ U x and the geodesic exp x (t ρ(u) v) for t ∈ [0; 1]: it is starting from x with tangent vector v and arriving at z = exp x (ρ(u) v) with the same tangent vector by parallel transportation. Reverting the time course (t → 1 -t), we have a geodesic starting at z with tangent vector -v and arriving at x with the same tangent vector. By definition of the cutting function, we have ρ(u) = ρ(u ) with u = (z, -v). Thus, the function q(u) = (exp x (ρ(u) v), -v) is a continuous bijection from U into itself with q (-1) = q.

Let z ∈ W ε . By definition of W ε , z is in the cut locus of a point x ∈ B(y, ε): there exists a unit tangent vector v at that point such that z = π(q(x, v)). Conversely, the cut locus of z intersects B(y, ε): there exists a unit tangent vector v at z such that x = π(q(z, v)) ∈ B(y, ε). Thus, we can rewrite W ε = π(U ε ) where U ε = q (-1) (π (-1) (B(y, ε)) ∩ U ). The functions π and q being continuous,

(i) For all x ∈ B(y, ε), C(x) ⊂ W ε , (ii) P (W ε ) = Wε dP (z) < α (iii) Wε dist(y, z) dP (z) < α
By hypothesis, the cut locus C(y) has a null measure for the measures dP (z). The distance being a measurable function, its measure is null on the cut locus: C(y) dist(y, z) dP (z) = 0. Thus, the functions Wε dP (z) and Wε dist(y, z) dP (z) are converging toward zero as ε goes to zero. By continuity, we can make both terms smaller than any positive α by choosing ε sufficiently small.

Bounding ∆ on W ε

Let W ε , ε and α verifying the conditions of lemma 2 and x, x ∈ B(y, ε). We have dist(z, x) ≤ dist(z, x ) + dist(x , x). Thus:

dist(z, x) 2 -dist(z, x ) 2 ≤ dist(x, x ) dist(x, x ) + 2 dist(z, x )
Using the symmetry of x and x and the inequalities dist(x, x ) ≤ 2ε and dist(z, x ) ≤ dist(z, y)+ε, we have: dist(z, x) 2 -dist(z, x ) 2 ≤ 2 dist(x, x ) 2ε + dist(z, y) . Applying this bound to x = exp y (v) and x = y, we obtain:

dist(z, exp y (v)) 2 -dist(z, y) 2 ≤ 2 2ε + dist(z, y) v
Now, the last term of ∆(z, v) is easily bounded by: -→ yz | v ≤ dist(y, z) v . Thus, we have:

|∆(z, v)| ≤ 4 ε + dist(z, y) v and its integral over W ε is bounded by: Wε ∆(z, v) dP (z) ≤ 4 v Wε ε + dist(z, y) dP (z) < 8α v
Bounding ∆ on M\W ε : Let x = exp y (v) ∈ B(y, ε). We know from lemma 2 that the cut locus C(x) of such a point belongs to W ε . Thus, the integration domain M\W ε is now independent of y and we can use the usual Lebesgue theorem to differentiate under the sum:

grad M\Wε dist(y, z) 2 dP (z) = -2 M\Wε -→ yz dP (z)
By definition, this means that for v small enough, we have:

M\Wε ∆(z, v) dP (z) < α v

Conclusion

Thus, for v small enough, we have M ∆(z, v) dP (z) < 9α v , which means that the variance has a derivative at the point y:

(grad σ 2 )(y) = -2 M/C(y) -→ yz dP (z)

B Approximation of the generalized Normal density

In this section, we only work with a normal coordinate system at the mean value of the considered normal law. The generalized Normal density is: N (y) = k exp (-y T Γ y/2), where the normalization constant, the covariance and concentration are related by:

k (-1) = M exp - y T Γ y 2 dM(y) and Σ = k M y y T exp - y T Γ y 2 dM(y)
We assume here that there two integrals converge toward a finite value for all concentration matrices Γ. In these expressions, the parameter is the concentration matrix Γ. The goal of this appendix is to invert these relations in order to obtain a Taylor expansion of the concentration matrix and the normalization coefficient k with respect to the variance. This can be realized thanks to the Taylor expansion of the Riemannian metric around the origin in a normal coordinate system [63, section 2.8, corollary 2.3, p.84]: det Let us investigate the normalization coefficient first. We have:

(g ij (exp v)) = 1 -Ric (v, v)/3 + O( v 3 ),
k (-1) = D exp - y T Γ y 2 1 - y T Ric y 6 + R dM (y) dy
where D is the definition domain of the exponential chart. Since the concentration matrix Γ is a symmetric positive definite matrix, we have a unique symmetric positive definite square root Γ 1/2 , which is obtained by changing the eigenvalues of Γ into their square root in a diagonalization. Using the change of variable z = Γ 1/2 y in the first two terms and the matrix T = Γ -1 2 Ric Γ -1 2 to simplify the notations, we have:

k (-1) = det(Γ) -1/2 D exp - z 2 2 1 - z T T z 6 dy + D exp - y T Γ y 2 R dM (y) dy
where the new definition domain is D = Γ 1/2 D. Likewise, we have for the covariance:

k (-1) Σ = 1 det(Γ) D Γ -1 2 (zz T )Γ -1 2 exp - z 2 2 1 - z T T z 6 dz+ D yy T exp - y T Γy 2 R dM (y)dy

B.1 Manifolds with an infinite injectivity radius at the mean point

At such a point, the cut locus is void and the definition domain of the exponential chart is D = D = R n . This simplifies the integrals so that we can compute the first two terms explicitly.

The normalization coefficient: The first term is

R n exp - z 2 2 dz = i R exp - z 2 i 2 dz i = (2π) n/2
By linearity, we also easily compute the second term:

R n z T T z exp - z 2 2 dz = Tr T R n z z T exp - z 2 2 dz = (2π) n/2 Tr(T )
Thus, we end up with:

k (-1) = (2π) n/2 det(Γ) 1 - Tr(Γ (-1) Ric ) 6 + R n exp - y T Γ y 2 R dM (y) dy
Unfortunately, the last term is not so simple to simplify as we only know that the remainder R dM (y) behaves as a O( y 3 ) for small values of y. The idea is to split the right integral into one part around zero, say for y < α, where we know how to bound the remainder, and to show that exp(-y 2 ) dominate the Remainder elsewhere. Let γ m be the smallest eigenvalue of Γ. As y T Γy ≥ γ m y 2 , we have:

R n exp - y T Γ y 2 R dM (y) dy ≤ R n exp -γ m y 2 2 R dM (y) dy
For any η > 0, there exists a constant α > 0 such that y < α implies that |R dM (y)| < η y 3 . Thus:

y <α exp -γ m y 2 2 R dM (y) dy < η y <α exp -γ m y 2 2 y 3 dy < η R n exp -γ m y 2 2 y 3 dy = 2 η V ol(S n-1 ) γ 2 m
Thus, we know that this part of the integral behaves as a O γ -2 m . For the other part, as exp is a monotonous function, one has: exp(-γ y 2 ) < exp(-y 2 )/γ 2 for y 2 > 2 log(γ)/(γ -1). Moreover, as the limit of log(γ)/(γ -1) is zero at γ = +∞, it can be made smaller than α/2 provided that γ is large enough. Thus, we have that (for γ m large enough):

y 2 >α exp - γ m y 2 2 R dM (y) dy < 1 γ 2 m y 2 >α exp - y 2 2 |R dM (y)| dy < 1 γ 2 m R n exp - y 2 2 |R dM (y)| dy
The last integral is a constant which is finite, since we assumed that k (-1) was finite for all Γ. Thus, this part of the integral behaves also as a O γ -2 m when γ m goes to infinity. We have obtained that:

R n exp(-y T Γ y/2) R dM (y) dy = O γ -2 m (4) so that k (-1) = (2π) n/2 det(Γ) 1 - Tr(Γ (-1) Ric ) 6 + O γ -2 m
As Tr(Γ (-1) Ric ) is a term in γ -1 m and det(Γ) a term in γ 1/2 m , we finally have:

k = det(Γ) (2π) n/2 1 + Tr(Γ (-1) Ric ) 6 + O γ -3/2 m (5)
The covariance matrix: The principle is the same as for the normalization coefficient, with slightly more complex integrals. We basically have to compute

k (-1) Σ = det(Γ) -1/2 Γ -1/2 I - 1 6 J Γ -1/2 + D y y T exp - y T Γ y 2 R dM (y) dy with I = R n z z T exp - z 2 2 dz and J = R n Tr(T z z T ) z z T exp - z 2 2 dz
These are classical calculations in multivariate statistics: for the first integral, the off diagonal elements I ij (with i = j) vanish because we integrate antisymmetric functions over R:

R n z i z j exp( z 2 /2) dz = 0 for i = j Summary of the approximate normal density: In a manifold without a cut locus at the mean point, the normal density in a normal coordinate system at the mean value is:

N (y) = k exp - y T Γ y 2
The normalization constant and the concentration matrices are approximated by the following expressions for a covariance matrix Σ of small variance σ 2 = Tr(Σ): B.2 Manifolds with a cut locus at the mean point Now, we have to integrate all the integrals over the definition domain D of the exponential chart at the mean point, which is only a subset of R n . For the integrals involving the remainder of the Taylor expansion of the metric, there is no problem as a bound on D is dominated by the bound on R n we already computed. For the other integrals, the idea is to bound each term from above and below to show that the same results still holds with just a slight modification of the Taylor expansion bounding term. Let r be the (finite) injectivity radius at the mean point. The open ball B(r) is the greatest geodesic ball included in the definition domain: B(r) ⊂ D ⊂ R n . The idea is to bound the integrals of positive terms by integrating over these three sets.

In fact, with the change in variable z = Γ 1/2 y, the ball is transformed into the ellipsoid Γ 1/2 B(r) = {z | z T Γ (-1) z < r 2 }. As z T Γ (-1) z ≤ γ (-1) m z 2 , we may use smaller ball B( √ γ m r) for the integrations in polar coordinates. For the integrals that are separable along axes, it will be simpler to use the maximal cube C = {z/ -√ γ m r/ √ n < z i < √ γ m r/ √ n}. As a summary, we have the inclusion chain: C ⊂ B( √ γ m r) ⊂ Γ 1/2 B(r) ⊂ D = Γ 1/2 D ⊂ R n .

In the following, we bound the different integrals used in the last section and summarize afterward the modifications it implies in the results of the last section. The first integral is bounded above by the previous value:

D exp - z 2 2 dz ≤ R n exp - z 2 2 dz = (2π) n/2
and below by an integration over C:

D exp - z 2 2 dz ≥ i √ γm r/ √ n - √ γm r/ √ n exp - z 2 i 2 dz i ≥ (2π) n/2 1 -erfc √ γ m r √ 2n n
Here, erfc (x) = 1 -2 √ π x 0 exp(-t 2 ) dt is the complement of the error function. An interesting property is that this function (like the exponential) tends toward zero faster than any fraction 1/x k as x goes to infinity. Putting things the other way, like for Taylor expansions, we will denote by ε(x) a function that is an O(x k ) for every positive k:

lim 0+ ε(x) x k = lim 0+ erfc (1/x) x k = lim 0+ exp(-1/x) x k = 0
We can summarize the value of the first integral by: D exp -z 2 2 dz = (2π) n/2 + ε γ m -1/2 r -1 appendix B. Let Σ 1/2 be the positive symmetric square root of Σ. With the change of variable y = Σ 1/2 x, this probability becomes:

Pr{χ 2 ≤ α 2 } = √ 2 k det(Σ) x ≤α exp - 1 2 x T Σ 1 2 Γ Σ 1 2 x 1 - 1 6 x T Σ 1 2 Ric Σ 1 2 x dx + O(σ 3 )
Using S = 1 3 Σ 1/2 Ric Σ 1/2 and the fact that k det(Σ) = (2 π) -n/2 (1 + O(σ 3 )) (from (Eq. 8)), we end up with

Pr{χ 2 ≤ α 2 } = 1 + O(σ 3 ) (2 π) n x ≤α exp - 1 2 x T Σ 1 2 Γ Σ 1 2 x 1 - 1 2 x T S x dx + O(σ 3 )
The goal is now to show that the above integrated term is basically exp(-x 2 /2) plus a remainder which integrates into a O(σ 3 ). In this process, we will have to compute the integrals:

I k (α) = x <α exp - x 2 2 x 2k dx
By changing into polar coordinates (r = x is the radius and u is the corresponding unit vector so that x = r u), we have dx = r n-1 dr du and thus:

I k (α) = S n-1 du α 0 r 2k+n-1 exp(-r 2 /2) dr = π n/2 Γ(n/2) α 2 0 t k-1+ n 2 exp - t 2 dt
with the change of variable t = r 2 . In this formula, Γ = +∞ 0 t x-1 exp(-t)dt is the Gamma function which can be recursively computed from Γ(x + 1) = x Γ(x) with Γ(1) = 1 and Γ(1/2) = √ π. The value I 0 (α) is in fact the standard cumulated probability function of the χ 2 law (up to a normalization factor). For α = +∞, the remaining integral can be computed using the Gamma function:

I k (+∞) = (2π) n 2 2 k Γ(k + n 2 ) Γ( n 2 )
Series expansion of the exponential From Eq. 7, we have Σ 1/2 Γ Σ 1/2 = Id -S -R Γ , where the remainder R Γ behaves as a O(σ 3 ). Thus, we have exp -

x T Σ 1/2 Γ Σ 1/2 x 2 = exp - 1 2 x 2 1 + x T S x 2 + R exp
where the remainder is given by the convergent series

R exp (x) = 1 2 x T R Γ x + +∞ k=2 (x T (S + R Γ ) x) k 2 k k! .
Substituting in the probability, we find that:

Pr{χ 2 ≤ α 2 } = 1 + O(σ 3 ) (2 π) n/2 x ≤α exp - x 2 2 (1 + R exp (x)) dx + O(σ 3 )
Integral of the Remainder R exp S = 1 3 Σ 1/2 Ric Σ 1/2 is a term of the order O(σ 2 ) and the remainder R Γ behaves as a O(σ 3 ). Thus, we can find positive constants C 1 and C 2 such that, for σ sufficiently small, we have: x T R Γ x ≤ C 1 σ 3 x 2 and x T (S + R Γ ) x ≤ C 2 σ 2 x 2 . This means that we can bound the integral of the remainder by:

x <α exp - x 2 2 R exp (x) dx ≤ C 1 σ 3 I 1 (+∞) + +∞ k=2 C k 2 σ 2k 2 k k! I k (+∞) =≤ C 1 σ 3 + (2 π) n 2 Γ n 2 ) +∞ k=2 a k
where a k = C k 2 σ 2k Γ(k + n/2)/k!. To show that the series converges, let us investigate the ratio of successive terms:

a k+1 a k = C 2 σ 2 k + 1 Γ(k + 1 + n/2) Γ(k + n/2) = C 2 σ 2 k + n/2 k + 1
The limit for k = +∞ is obviously C 2 σ 2 , which can be made smaller than 1 for σ sufficiently small. Thus, by d'Alembert's ratio test, the series converges. Moreover, the smaller term of the series is a σ 4 (for k=2), which shows that the integral of the remainder is a finally dominated by the first term, a O(σ 3 ).

Conclusion Finally, the researched probability is

Pr{χ 2 ≤ α 2 } = (2π) -n 2 x ≤α exp - x 2 2 dx+O(σ 3 ) = 1 2 Γ( n 2 ) α 2 0 t 2 n/2
exp(-t/2)dt+O(σ 3 ) Thus, the probability density function of a χ 2 is: and we have already seen that this integral is 1 + ε σ r . As α > √ γ m r, the same integral is itself of the same order, and we obtain in all cases the same result as with no cut locus where O(σ 3 ) is replaced by O(σ 3 ) + ε (σ/r).

p χ 2 (u) = 1 2 Γ n

Figure 1 :

 1 Figure 1: Exponential chart and cut locus for the sphere S 2 and the projective space P 2

Definition 2

 2 Let B(M) be the Borel σ-algebra of M. The random point x has a probability density function p x (real, positive and integrable function) if: ∀X ∈ B(M), Pr(x ∈ X ) = X p(y) dM(y) and Pr(M) = M p(y) dM(y) = 1

Figure 2 :

 2 Figure 2: The covariance is defined in the tangent plane at S 2 at the mean point as the classical covariance matrix of the random vector "deviation from the mean" Σ xx = E

Lemma 1 W

 1 ε = x∈B(y,ε) C(x) is a continuous series of included and decreasing open sets all containing C(y) and converging towards it.

Lemma 2

 2 U ε is open. Finally, π being an open map, we conclude that W ε is open. Let y ∈ M and α > 0. Then there exists an open neighborhood W ε of C(y) such that

  where Ric is the expression of the Ricci tensor of scalar curvatures in the exponential chart. Thus[42, p.144]:dM(y) = det(G(y)) dy = 1 -y T Ric y 6 + R dM (y)

C. 2 2 B

 22 Manifold with a cut locus at the mean pointIf there is a cut locus, we have to replace O(σ k ) by O(σ k ) + ε (σ/r) and we should only integrate on B Σ (α) ∩ D. After the change of coordinates, we should integrate on B(α) ∩ D . Thus, we have:Pr{χ 2 ≤ α 2 } = (2π) -n As we have B( √ γ m r) ⊂ D ⊂ R n , we can enclose the researched domain into: B min( √ γ m r, α) ⊂ B(α) ∩ D ⊂ B(α).For α ≤ √ γ m r, there is not problem but for α > √ γ m r, we have:Pr{χ 2 ≤ α 2 } ≥ (2π) -

Here, strongly convex means that for every two points of U there is a unique minimizing geodesic joining them that depend in a C ∞ of the two points.

Acknowledgments

The main part of this work was done at MIT, Artificial Intelligence Lab in 1997. If all the ideas presented in this paper were in place at that time [69], the formula for gradient of the variance (appendix A) remained a conjecture. The author would like to thank specially Pr Maillot for providing a first proof of the differentiability of the variance for the uniform distribution on compact manifolds [53]. Its generalization to the non compact case with arbitrary distributions considerably delayed the submission of this paper. In the meantime, other (and simpler) proofs were proposed [15,38]. The author would also like to thank two anonymous referees for very valuable advises and for the detection of a number of technical errors.

A Gradient of the variance

This proof is a generalization of a differentiability proof for the uniform distribution on compact manifolds by Pr Maillot [START_REF] Maillot | Différentielle de la variance et centrage de la plaque de coupure sur une variété riemannienne compacte[END_REF]. The main difficulties were to remove the compactness hypothesis, and to generalize to arbitrary distributions.

The diagonal elements are:

exp(-z 2 j /2) dz j = (2π) n/2 so that we have I = (2π) n/2 Id. Since Tr(T z z T ) = k,l T kl z k z l , the elements of the second integral are:

Let us investigate first off diagonal elements (i = j): for k = i and l = j or l = i and k = j, we integrate an antisymmetric function and the result is zero. Since the Ricci curvature matrix is symmetric, the matrix T is also symmetric and the sum is reduced to a single term:

The diagonal terms are:

For k = l, we integrate antisymmetric functions: the result is zero. Thus, we are left with:

Combining with off diagonal terms, we get: J = (2π) n/2 (Tr(T ) Id + 2 T ), so that:

Like for the normalization constant, we now have to bound the integral of the remainder. The principle is the same. Let γ m be the smallest eigenvalue of Γ. As y T Γy ≥ γ m y 2 , we have:

The last integral is a constant which is finite, since we assumed that Σ was finite for all Γ. Thus, this part of the integral behaves also as a O γ -3 m when γ m goes to infinity. As 1/ det(Γ) a term in γ -1/2 m , we finally have:

Simplifying this expression using the previous Taylor expansion of k (Eq. 5), we obtain thus the following relation between the covariance and the concentration matrices:

Inverting the variables of the Taylor expansions This relation can be inverted to obtain the Taylor expansion of the concentration matrix with respect to the covariance. First, we shall note that from the above equation, O(γ

where σ max is the square root of the largest eigenvalue of Σ. However, a global variable such as the variance

Then, one easily verifies that the Taylor expansion of Γ (-1) is: Γ (-1) = Σ + Σ Ric Σ/3 + O(σ 5 ), and finally:

To express k with respect to Σ instead of Γ, we have to compute Tr(Γ (-1) Ric ) and det(Γ). For the determinant, one verifies that if A is a differential map from R to the matrix group, one has det(A) = (det(A)) Tr(A A (-1) ), so that det( Id + η B) = 1 + η Tr(B) + O(η 2 ). Since Γ Σ = Id -1 3 Ric Σ + O(σ 3 ) and Σ is a term in O(σ 2 ), we have

and thus

Substituting this expression in equation 5, we obtain:

Which simplifies in:

Now, for the second integral, we have

We obtain thus

In fact, it is easy to see that every integral we computed over R n in the previous paragraph has the same value over D plus a term whose absolute value is of the order of ε(γ

). Thus, we can directly generalize the previous results by replacing O(σ k ) with O(σ k ) + ε (σ/r).

Summary of the approximate normal density: In a manifold with injectivity radius r at the mean point, the normal density in a normal coordinate system at this mean value is:

The normalization constant and the concentration matrices are approximated by the following expressions for a covariance matrix Σ of small variance σ 2 = Tr(Σ):

Assuming that the random point x ∼ (x, Σ) follows a normal law, we want to compute the probability that

Let B Σ (α) be the "elliptic ball" of covariance Σ and radius α: B Σ (α) = {y / y T Σ (-1) y ≤ α 2 }.

C.1 Manifold without a cut locus at the mean point

Assuming that there is not cut locus, this probability can be written in a normal coordinate system: