Elsa Cazelles

Vivien Seguy

Jérémie Bigot

Marco Cuturi

Nicolas Papadakis

Log-PCA versus Geodesic PCA of histograms in the Wasserstein space

Keywords: Geodesic Principal Componant Analysis, Wasserstein Space, Non-convex optimization AMS classifications: 62-07, 68R10, 62H25

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Most datasets describe multivariate data, namely vectors of relevant features that can be modeled as random elements sampled from an unknown distribution. In that setting, Principal Component Analysis (PCA) is certainly the simplest and most widely used approach to reduce the dimension of such datasets. We consider in this work the statistical analysis of data sets whose elements are histograms supported on the real line, and also discuss extensions to the general case of probability measures supported on the d-dimensional Euclidean space. Just as with PCA, our main goal in that setting is to compute the principal modes of variation of histograms around their mean element and therefore facilitate the visualization of such datasets. (bottom-left) Geodesic PCA with respect to the Wasserstein metric using the iterative geodesic algorithm (4.1). The black curve represents the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

However, since the number, size or locations of significant bins in the histograms of interest may vary from one histogram to another, using standard PCA on histograms (with respect to the Euclidean metric) is bound to fail (see for instance Figure 1). In this paper, we propose to use the 2-Wasserstein metric [Vil03, §7.1] to measure the distance between histograms, and to compute their modes of variation accordingly. In our approach, histograms are seen as piecewise constant probability density functions (pdf) supported on a given interval Ω of the real line. In this setting, the variability in a set of histograms can be analyzed via the notion of Geodesic PCA (GPCA) of probability measures in the Wasserstein space W 2 (Ω) admitting these histograms as pdf. That approach has been recently proposed in the statistics and machine learning literature in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] for probability measures on the real line, and in [SC15, WSB + 13] for discrete probability measures on R d . However, implementing GPCA remains a challenging computational task even in the simplest case of pdf's supported on R. The purpose of this paper is to provide a fast algorithm to perform GPCA of probability measures supported on the real line, and to compare its performances with log-PCA, namely standard PCA in the tangent space at the Wasserstein barycenter of the data [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF][START_REF] Petersen | Functional data analysis for density functions by transformation to a hilbert space[END_REF].

Related results

Foundations of Geodesic PCA in the Wasserstein space The space of probability measures (with finite second moment) endowed with the 2-Wasserstein distance is not a Hilbert space. Therefore, standard PCA, which involves computing a covariance matrix, cannot be applied directly to compute principal mode of variations in a Wasserstein sense. Nevertheless, a meaningful notion of PCA can still be defined by relying on the pseudo-Riemannian structure of the Wasserstein space, which was extensively studied in [START_REF] Ambrosio | Gradient flows with metric and differentiable structures, and applications to the Wasserstein space[END_REF] and [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF]. Following this principle, a framework for GPCA of probability measures supported on a interval Ω ⊂ R was introduced in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF]. GPCA is defined as the problem of estimating a principal geodesic subspace (of a given dimension) which maximizes the variance of the projection of the data to that subspace. In that approach the base point of that subspace is the Wasserstein barycenter of the data as introduced in [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF], which is also known as a Fréchet mean in the general context of metric spaces. Existence, consistency and a detailed characterization of GPCA in W 2 (Ω) were studied in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF]. In particular, the authors have shown that this approach is equivalent to map the data in the tangent space of W 2 (Ω) at the Fréchet mean, and then to perform a PCA in this Hilbert space that is constrained to lie in a convex and closed subset of functions. Mapping the data to this tangent space is not difficult in the one-dimensional case as it amounts to computing a set of optimal maps between the data and their Wasserstein barycenter, for which a closed form is available using their quantile functions (see for example [START_REF] Villani | Topics in Optimal Transportation[END_REF] §2.2]). To perform PCA on the mapped data, [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] fell short of proposing an algorithm to minimize that problem, which has a non-convex and non-differentiable objective function as well as involved constraints. Only a numerical approximation to the computation of GPCA was proposed in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF], which amounts to applying log-PCA, namely a standard PCA of the dataset mapped beforehand to the tangent space of W 2 (Ω) at its Fréchet mean.

Previous work in the one-dimensional case PCA of histograms with respect to the Wasserstein metric has also been proposed in [START_REF] Verde | Dimension reduction techniques for distributional symbolic data[END_REF] in the context of symbolic data analysis. Their approach consists in computing a standard PCA in the Hilbert space L 2 ([0, 1]) of the quantile functions associated to the histograms. Therefore, the algorithm in [START_REF] Verde | Dimension reduction techniques for distributional symbolic data[END_REF] corresponds to log-PCA of probability measures as suggested in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF], but it does not solve the problem of convex-constrained PCA in a Hilbert space associated to an exact GPCA in W 2 (Ω). A related problem, which can be referred to as geodesic regression (considered in [START_REF] Fletcher | Geodesic regression on riemannian manifolds[END_REF][START_REF] Thomas | Geodesic regression and the theory of least squares on riemannian manifolds[END_REF] for data on a Riemannian manifold), has been considered by Jiang et al. in [START_REF] Jiang | Geometric methods for spectral analysis[END_REF] where the authors fit a single geodesic g t to indexed histograms in order to model nonstationary time series. In the problem of finding principal geodesics, we do not assume that the dataset is indexed.

PGA and log-PCA on Riemannian manifolds The method of GPCA proposed in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] clearly shares similarities with analogs of PCA for data belonging to a Riemannian manifold M of finite dimension. These methods, generally referred to as Principal Geodesic Analysis (PGA), extend the notion of classical PCA in Euclidean spaces for the purpose of analyzing data belonging to curved Riemannian manifolds (see e.g. [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF][START_REF] Sommer | Manifold valued statistics, exact principal geodesic analysis and the effect of linear approximations[END_REF]). This generalization of PCA proceeds by replacing Euclidean concepts of vector means, lines and orthogonality by the more general notions in Riemannian manifolds of Fréchet mean, geodesics, and orthogonality in tangent spaces.

In [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF], linearized PGA, which we refer to as log-PCA, is defined as follows. In a first step, data are mapped to the tangent space T xM at their Fréchet mean x by applying the logarithmic map log x to each data point. Then, in a second step, standard PCA in the Euclidean space T xM can be applied. This provides a family of orthonormal tangent vectors. Principal components of variation in M can then be defined by back-projection of these tangent vectors on M by using the exponential map at x, that is known to parameterize geodesics at least locally. Log-PCA has low computational cost, but this comes at the expense of two simplifications and drawbacks:

(1) First, log-PCA amounts to substituting geodesic distances between data points by the linearized distance in T xM, which may not always be a good approximation because of the curvature of M, see e.g. [START_REF] Sommer | Manifold valued statistics, exact principal geodesic analysis and the effect of linear approximations[END_REF].

(2) Secondly, the exponential map at the Fréchet mean parameterizes geodesics only locally, which implies that principal components in M obtained with log-PCA may not be geodesic along the typical range of the dataset.

Numerical approaches to GPCA and log-PCA in the Wasserstein space Computational methods have been introduced in [SC15, WSB + 13] to extend the concepts of PGA on Riemannian manifolds to that of the space W 2 (R d ) of probability measures supported on R d endowed with the Wasserstein metric. [WSB + 13] propose to compute a notion of template measure (using k-means clustering) of a set of discrete probability measures, and to consider then the optimal transport plans from that template measure to each measure in the data set. Computation of the barycentric projection of each optimal transport plan leads to a set of Monge maps over which a standard PCA can be applied, resulting in an orthonormal family of tangent vectors defined on the support of the template measure. Principal components of variation in R d can then be obtained through the push-forward operator, namely by moving the mass along these tangent vectors. This approach, analog to log-PCA on Riemannian manifolds, suffers from the main drawbacks mentioned above: for d > 1, the linearized Wasserstein distance may be a crude approximation of the Wasserstein distance, and there is no guarantee that the computed tangent vectors parameterize geodesics of sufficient length to summarize most of the variability in the dataset. Losing geodesicity means that the principal components are curves in W 2 (R d ) along which the mass may not be transported optimally, which may significantly reduce the interpretability of these principal components. A different approach was proposed in [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF], in which the notion of generalized geodesics in W 2 (R d ) (see e.g. Chapter 9 in [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF]) is used to define a notion of PGA of discrete probability measures. In [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF], generalized geodesics are parameterized using two velocity fields defined on the support of the Wasserstein barycenter. The authors proposed to minimize directly the distances from the measures in the dataset to these generalized geodesics, by updating these velocity fields which are constrained to be in opposite directions. This approach is more involved computationally than log-PCA, but it avoids some of the drawbacks highlighted above. Indeed, the resulting principal components yield curves in W 2 (R d ) that are approximately geodesics. Nevertheless, the computational method in [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] uses a heuristic projection on the set of optimal velocity fields, which results in a algorithm which has no convergence guarantees. Moreover, by optimizing over generalized geodesics rather than geodesics, it does not solve exactly the problem of computing geodesic PCA in W 2 (R d ).

Main contributions

In this paper, we focus on computing an exact GPCA on probability measures supported on Ω ⊂ R d . We mainly focus on the case d = 1 (discussing extensions in the last section), which has the advantage that the linearized Wasserstein distance in the tangent space is equal to the Wasserstein distance in the space W 2 (Ω). The main challenge is thus to obtain principal curves which are geodesics along the range of the dataset.

The first contribution of this paper is to propose two fast algorithms for GPCA in W 2 (Ω). The first algorithm finds iteratively geodesics such that the Wasserstein distance between the dataset and the parameterized geodesic is minimized with respect to W 2 (Ω). This approach is thus somewhat similar to the one in [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF]. However, a heuristic barycentric projection is used in [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] to remain in the feasible set of constraints during the optimization process. In our approach, we rely on proximal operators of both the objective function and the constraints to obtain an algorithm which is guaranteed to converge to a critical point of the objective function. Moreover, we show that the global minimum of our objective function for the first principal geodesic curve corresponds indeed to the solution of the exact GPCA problem defined in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF]. While this algorithm is able to find iteratively orthogonal principal geodesics, there is no guarantee that several principal geodesics parameterize a surface which is also geodesic. This is the reason we also propose a second algorithm which computes all the principal geodesics at once by parameterizing a geodesic surface as a convex combination of optimal velocity fields and relaxing the orthogonality constraint between principal geodesics. Both algorithms are a variant of the proximal Forward-Backward algorithm. They converge to a stationary point of the objective function, as shown by recent results in non-convex optimization based on proximal methods [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems[END_REF][START_REF] Ochs | ipiano: Inertial proximal algorithm for nonconvex optimization[END_REF]. Our second contribution is a numerical comparison of log-PCA in W 2 (Ω), as done in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] (for d = 1) or [WSB + 13], with our approach which solves the exact Wasserstein GPCA problem. Finally, we discuss extensions to the case of probability measures supported on the d-dimensional Euclidean space, providing detailed calculations in the two-dimensional case, and perform computation of GPCA on a two-dimensional example, comparing results with the ones obtained with the log-PCA approach.

In all our experiments, data are normalized in order to have a suitable representation as probability measures. We believe this preprocessing does not affect any useful properties of the histogram datasets considered in the present article, in the same way as centering or whitening are often used as a preprocessing step in many data-analysis tasks. Yet, if the total mass of a given histogram matters for some application, we could consider the use of unbalanced optimal transport [CPSV15, [START_REF] Liero | Optimal entropy-transport problems and a new hellinger-kantorovich distance between positive measures[END_REF][START_REF] Chizat | An interpolating distance between optimal transport and fisher-rao metrics[END_REF] which provides a distance between unnormalized measures. This generalization is out of the scope of this work and may be an interesting line of research in the future.

The code to reproduce the results of this paper is available online1 .

Structure of the paper

In Section 2, we provide some background on GPCA in the Wasserstein space W 2 (Ω), borrowing material from previous work in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF]. Section 3 describes log-PCA in W 2 (Ω), and some of its limitations are discussed. Section 4 contains the main results of the paper, namely two algorithms for computing GPCA. In Section 5, we provide a comparison between GPCA and log-PCA using statistical analysis of real datasets of histograms. In the last Section 6 we discuss extensions of our algorithms to the case d > 1, and perform GPCA computation on a twodimensional example, comparing again results with the log-PCA approach. Some perspectives on this work are also given. Finally, various details on the implementation of the algorithms are deferred to technical Appendices.

2 Background on Geodesic PCA in the Wasserstein space

Definitions and notations

Let Ω be a (possibly unbounded) interval in R. Let ν be a probability measure (also called distribution) over (Ω, B(Ω)) where B(Ω) is the σ-algebra of Borel subsets of Ω. For a mapping T : Ω → Ω, the push-forward measure T #ν is a probability measure on Ω defined by (T #ν)(A) = ν{x ∈ Ω|T (x) ∈ A}, for any A ∈ B(Ω). The cumulative distribution function (cdf) and the (generalized) quantile function of ν are denoted respectively by F ν and F - ν . The Wasserstein space W 2 (Ω) is the set of probability measures with support included in Ω and having a finite second moment, that is endowed with the quadratic Wasserstein distance d W defined by

d 2 W (µ, ν) := 1 0 (F - µ (α) -F - ν (α)) 2 dα, µ, ν ∈ W 2 (Ω). (2.1)
We also denote by W ac 2 (Ω) the set of measures ν ∈ W 2 (Ω) that are absolutely continuous with respect to the Lebesgue measure dx on R. If µ ∈ W ac 2 (Ω) then T * = F - ν •F µ will be referred to as the optimal mapping to push-forward µ onto ν and in this case d 2 W (µ, ν) = Ω (T * (x)-x) 2 dµ(x). For a detailed analysis of W 2 (Ω) and its connection with optimal transport theory, we refer to [START_REF] Villani | Topics in Optimal Transportation[END_REF].

The pseudo Riemannian structure of the Wasserstein space

In what follows, µ r denotes a reference measure in W ac 2 (Ω), whose choice will be discussed later on. The space W 2 (Ω) has a formal Riemannian structure described, for example, in [START_REF] Ambrosio | Gradient flows with metric and differentiable structures, and applications to the Wasserstein space[END_REF]. The tangent space at µ r is defined as the Hilbert space L 2 µr (Ω) of real-valued, µ r -squareintegrable functions on Ω, equipped with the inner product ⟨•, •⟩ µr defined by ⟨u,

v⟩ µr = Ω u(x)v(x)dµ r (x), u, v ∈ L 2
µr (Ω), and associated norm ∥ • ∥ µr . We define the exponential and the logarithmic maps at µ r , as follows.

Definition 2.1. Let id : Ω → Ω be the identity mapping. The exponential exp µr : L 2 µr (Ω) → W 2 (R) and logarithmic log µr : W 2 (Ω) → L 2 µr (Ω) maps are defined respectively as

exp µr (v) = (id + v)#µ r and log µr (ν) = F - ν • F µr -id. (2.2)
Contrary to the setting of Riemannian manifolds, the "exponential map" exp µr defined above is not a local homeomorphism from a neighborhood of the origin in the "tangent space" L 2 µr (Ω) to the space W 2 (Ω), see e.g. [START_REF] Ambrosio | Gradient flows with metric and differentiable structures, and applications to the Wasserstein space[END_REF]. Nevertheless, it is shown in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] that exp µr is an isometry when restricted to the following specific set of functions

V µr (Ω) := log µr (W 2 (Ω)) = log µr (ν) ; ν ∈ W 2 (Ω) ⊂ L 2 µr (Ω),
and that the following results hold (see [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF]).

Proposition 2.1. The subspace V µr (Ω) satisfies the following properties :

(P1) the exponential map exp µr restricted to V µr (Ω) is an isometric homeomorphism, with inverse log µr . We have hence

d W (ν, η) = || log µr (ν) -log µr (η)|| L 2 µr (Ω) . (P2) the set V µr (Ω) := log µ (W 2 (Ω)
) is closed and convex in L 2 µr (Ω).

(P3) the space V µr (Ω) is the set of functions v ∈ L 2 µr (Ω) such that T := id + v is µ r -almost everywhere non decreasing and that T (x) ∈ Ω, for x ∈ Ω.

Moreover, it follows, from [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF], that geodesics in W 2 (Ω) are exactly the image under exp µr of straight lines in V µr (Ω). This property is stated in the following lemma.

Lemma 2.1. Let γ : [0, 1] → W 2 (Ω) be a curve and let v 0 := log µr (γ(0)), v 1 := log µr (γ(1)). Then γ = (γ t ) t∈[0,1] is a geodesic if and only if γ t = exp µr ((1 -t)v 0 + tv 1 ), for all t ∈ [0, 1].

GPCA for probability measures

Let ν 1 , . . . , ν n be a set of probability measures in W ac 2 (Ω). Assuming that each ν i is absolutely continuous simplifies the following presentation, and it is in line with the purpose of statistical analysis of histograms. We define now the notion of (empirical) GPCA of this set of probability measures by following the approach in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF]. The first step is to choose the reference measure µ r . To this end, let us introduce the Wasserstein barycenter [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF] or Fréchet mean of the ν i 's, that is defined as the probability measure ν,

ν = argmin µ∈W 2 (Ω) 1 n n i=1 d 2 W (ν i , µ).
Note that it immediately follows from results in [AC11] that ν ∈ W ac 2 (Ω), and that its cdf satisfies

F - ν = 1 n n i=1 F - ν i . (2.3)
A typical choice for the reference measure is to take µ r = ν which represents an average location in the data around which can be computed the principal sources of geodesic variability. To introduce the notion of a principal geodesic subspace of the measures ν 1 , . . . , ν n , we need to introduce further notation and definitions. Let G be a subset of W 2 (Ω). The distance between

µ ∈ W 2 (Ω) and the set G is d W (ν, G) = inf λ∈G d W (ν, λ)
, and the average distance between the data and G is taken as

D W (G) := 1 n n i=1 d 2 W (ν i , G). (2.4) Definition 2.2. Let K be some positive integer. A subset G ⊂ W 2 (Ω) is said to be a geodesic set of dimension dim(G) = K if log µr (G) is a convex set such that the dimension of the smallest affine subspace of L 2 µr (Ω) containing log µr (G) is of dimension K.
The notion of principal geodesic subspace (PGS) with respect to the reference measure µ r = ν can now be presented below.

Definition 2.3. Let CL(W ) be the metric space of nonempty, closed subsets of W 2 (Ω), endowed with the Hausdorff distance, and

CG ν,K (W ) = {G ∈ CL(W ) | ν ∈ G, G is a geodesic set and dim(G) ≤ K} , K ≥ 1. A principal geodesic subspace (PGS) of ν = (ν 1 , . . . , ν n ) of dimension K with respect to ν is a set G K ∈ argmin G∈CG ν,K (W ) D W (G) = argmin G∈CG ν,K (W ) 1 n n i=1 d 2 W (ν i , G).
(2.5)

When K = 1 , searching for the first PGS of ν simply amounts to search for a geodesic curve γ (1) that is a solution of the following optimization problem:

γ(1) := argmin γ 1 n n i=1 d 2 W (ν i , γ) | γ is a geodesic in W 2 (Ω) passing through µ r = ν. .
We remark that this definition of γ(1) as the first principal geodesic curve of variation in W 2 (Ω) is consistent with the usual concept of PCA in a Hilbert space in which geodesic are straight lines.

For a given dimension k, the GPCA problem consists in finding a nonempty closed geodesic subset of dimension k which contains the reference measure µ r and minimizes Eq. (2.4). We describe in the next section how we can parameterize such sets G.

Geodesic PCA parameterization

GPCA can be formulated as an optimization problem in the Hilbert space L 2 ν (Ω). To this end, let us define the functions ω i = log ν (ν i ) for 1 ≤ i ≤ n that corresponds to the data mapped in the tangent space. It can be easily checked that this set of functions is centered in the sense that 1 n n i=1 ω i = 0. Note that, in a one-dimensional setting, computing ω i (mapping of the data to the tangent space) is straightforward since the optimal maps T * i = F - ν i • F ν between the data and their Fréchet mean are available in a simple and closed form.

For U = {u 1 , . . . , u K } a collection of K ≥ 1 functions belonging to L 2 ν (Ω), we denote by Sp(U) the subspace spanned by u 1 , . . . , u K . Defining Π Sp(U ) v as the projection of v ∈ L 2 ν (Ω) onto Sp(U), and Π Sp(U )∩Vν (Ω) v as the projection of v onto the closed convex set Sp(U) ∩ V ν (Ω), then we have

Proposition 2.2. Let ω i = log ν (ν i ) for 1 ≤ i ≤ n, and U * = {u * 1 , . . . , u * k } be a minimizer of 1 n n i=1 ∥ω i -Π Sp(U )∩Vν (Ω) ω i ∥ 2 ν , (2.6
)

over orthonormal sets U = {u 1 , . . . , u K } of functions in L 2 ν (Ω) of dimension K (namely such that ⟨u j , u j ′ ⟩ ν = 0 if j ̸ = j ′ and ∥u j ∥ ν = 1). If we let G U * := exp ν (Sp(U * ) ∩ V ν (Ω)),
then G U * is a principal geodesic subset (PGS) of dimension k of the measures ν 1 , . . . , ν n , meaning that G U * belongs to the set of minimizers of the optimization problem (2.5).

Proof. For v ∈ L 2 ν (Ω) and a subset C ∈ L 2 ν (Ω), we define d ν (v, C) = inf u∈C ∥v -u∥ ν . Remark that i ω i = 0. Hence by Proposition 3.3 in [BGKL15], if U * minimizes 1 n n i=1 d 2 ν (ω i , Sp(U * ) ∩ V ν (Ω)) = 1 n n i=1 ∥ω i -Π Sp(U * )∩Vν (Ω) ω i ∥ 2 ν , then Sp(U * ) ∩ V ν (Ω) ∈ argmin C 1 n n i=1 d 2 ν (ω i , C
), where C is taken over all nonempty, closed, convex set of V ν (Ω) such that dim(C) ≤ K and 0 ∈ C. By Proposition 4.3 in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF], and since log ν (ν) = 0, we can conclude that G * is a geodesic subset of dimension K which minimizes (2.4).

Thanks to Proposition 2.2, it follows that GPCA in W 2 (Ω) corresponds to a mapping of the data into the Hilbert space L 2 ν (Ω) which is followed by a PCA in L 2 ν (Ω) that is constrained to lie in the convex and closed subset V ν (Ω). This has to be interpreted as a geodesicity constraint coming from the definition of a PGS in W 2 (Ω). Because this geodesicity constraint is nontrivial to implement, recent works about geodesic PCA in W 2 (Ω) relied on a heuristic projection on the set of optimal maps [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF], or relaxed the geodesicity constraint by solving a linearized PGA [WSB + 13, BGKL15]. We describe the latter approach in the following section.

The log-PCA approach

For data in a Riemannian manifold, we recall that log-PCA consists in solving a linearized version of the PGA problem by mapping the whole data set to the tangent space at the Fréchet mean through the logarithmic map [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF]. This approach is computationally attractive since it boils down to computing a standard PCA. [WSB + 13] used this idea to define a linearized PGA in the Wasserstein space W 2 (R d ), by defining the logarithmic map of a probability measure as the barycentric projection of an optimal transport plan with respect to a template measure. This approach has the two drawbacks (1) and (2) of log-PCA mentioned in Section 1.1. A third limitation inherent to the Wasserstein space is that when this template probability measure is discrete, the logarithmic map cannot be defined straightforwardly as there is no guarantee about the existence of an optimal map solution of the optimal transport problem. This is why the authors of [WSB + 13] had to compute the barycentric projection of each optimal transport plan, which is obtained by simply averaging the locations of the split mass defined by this plan. This averaging process is however lossy as distinct probability measures can have the same barycentric projection.

We consider as usual a subset Ω ⊂ R. In this setting, W 2 (Ω) is a flat space as shown by the isometry property (P1) of Proposition 2.1. Moreover, if the Wasserstein barycenter ν is assumed to be absolutely continuous, then Definition 2.1 shows that the logarithmic map at ν is well defined everywhere. Under such an assumption, log-PCA in W 2 (Ω) corresponds to the following steps:

1. compute the log-maps (see Definition 2.1)

ω i = log ν (ν i ), i = 1, . . . , n, 2. perform the PCA of the projected data ω 1 , • • • , ω n in the Hilbert space L 2 ν (Ω) to obtain K orthogonal directions ũ1 , . . . , ũK in L 2 ν (Ω) of principal variations,
3. recover a principal subspace of variation in W 2 (Ω) with the exponential map exp ν (Sp( Ũ)) of the principal eigenspace Sp( Ũ) in L 2 ν (Ω) spanned by ũ1 , . . . , ũK . For specific datasets, log-PCA in W 2 (Ω) may be equivalent to GPCA, in the sense that the set exp ν (Sp( Ũ) ∩ V ν (Ω)) is a principal geodesic subset of dimension K of the measures ν 1 , . . . , ν n , as defined by (2.5). Informally, this case corresponds to the setting where the data are sufficiently concentrated around their Wasserstein barycenter ν (we refer to Remark 3.5 in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] for further details). However, carrying out a PCA in the tangent space of W 2 (R) at ν is a relaxation of the convex-constrained GPCA problem (2.6), where the elements of the sought principal subspace do not need to be in V ν . Indeed, standard PCA in the Hilbert space L 2 ν (Ω) amounts to finding Ũ = {ũ 1 , . . . , ũK } minimizing,

1 n n i=1 ∥ω i -Π Sp(U ) ω i ∥ 2 ν , , (3.1)
over orthonormal sets U = {u 1 , . . . , u k } of functions in L 2 ν (Ω). It is worth noting that the three steps of log-PCA in W 2 (Ω) are simple to implement and fast to compute, but that performing log-PCA or GPCA (2.6) in W 2 (Ω) are not necessarily equivalent.

Log-PCA is generally used for two main purposes. The first one is to obtain a low dimensional representation of each data measure ν i = exp ν (ω i ) through the coefficients ⟨ω i , ũk ⟩ L 2 ν . From this low dimensional representation, the measure ν i ∈ W 2 (Ω) can be approximated through the exponential mapping exp ν (Π Sp(U ) ω i ). The second one is to visualize each mode of variation in the dataset, by considering the evolution of the curve t → exp ν (tũ k ) for each ũk ∈ Ũ.

However, relaxing the convex-constrained GPCA problem (2.6) when using log-PCA results in several issues. Indeed, as shown in the following paragraphs, not taking into account this geodesicity constraint makes difficult the computation and interpretation of exp ν (Sp( Ũ)) as a principal subspace of variation, which may limit its use for data analysis.

Numerical implementation of pushforward operators A first downside to the log-PCA approach is the difficulty of the numerical implementation of the pushforward operator in the exponential map exp ν (v) = (id + v)#ν when the mapping id + v is not a strictly increasing function for a given vector v ∈ Sp( Ũ). This can be shown with the following proposition, which provides a formula for computing the density of a pushforward operator.

Proposition 3.1. (Density of the pushforward) Let µ ∈ W 2 (R) be an absolutely continuous measure with density ρ (that is possibly supported on an interval Ω ⊂ R). Let T : R → R be a differentiable function such that |T ′ (x)| > 0 for almost every x ∈ R, and define ν = T #µ. Then, ν admits a density g given by,

g(y) = x∈T -1 (y) ρ(x) |T ′ (x)| , y ∈ R. (3.2)
When T is injective, this simplifies to,

g(y) = ρ(T -1 (y)) |T ′ (T -1 (y))| . ( 3
h(x)|T ′ (x)|dx = R x∈T -1 (y) h(x)dy. (3.4) Let B a Borel set and choose h(x) = ρ(x) |T ′ | 1 T -1 (B) ,
x. Hence, using (3.4), one obtains that

T -1 (B) ρ(x)dx = R x∈T -1 (y) ρ(x) |T ′ (x)| 1 T -1 (B) (x)dy = B x∈T -1 (y) ρ(x) |T ′ (x)|
dy.

The definition of the pushforward ν(B) = µ(T -1 (B)) then completes the proof.

The numerical computation of formula (3.2) or (3.3) is not straightforward. When T is not injective, computation of the formula (3.2) must be done carefully by partitioning the domain of T in sets on which T is injective. Such a partitioning depends on the method of interpolation for estimating a continuous density ρ from a finite set of its values on a grid of reals. More importantly, when T ′ (x) is very small, ρ(x)

|T ′ (x)| may become very irregular and the density of ν = T #µ may exhibit large peaks, see Figure 2 for an illustrative example.

Pushforward of the barycenter outside the support Ω A second downside of log-PCA in W 2 (Ω) is that the range of the mapping Ti = id + Π Sp( Ũ ) ω i may be larger than the interval Ω. This implies that the density of the pushforward of the Wasserstein barycenter ν by this mapping, namely exp ν (Π Sp( Ũ ) ω i ), may have a support which is not included in Ω. This issue may be critical when trying to estimate the measure ν i = exp ν (ω i ) by its projected measure exp ν (Π Sp( Ũ ) ω i ). For example, in a dataset of histograms with bins necessarily containing only nonnegative reals, a projected distribution with positive mass on negative reals would be hard to interpret. Mapping T i = id+Π Sp(ũ2) ω i obtained from the projections of these 5 distributions onto the second eigenvector ũ2 given by log-PCA of the whole dataset. (Right) Pushforward exp ν (Π Sp(ũ2) ω i ) = T i #ν of the Wasserstein barycenter ν for each 1 ≤ i ≤ 5. As the derivative T ′ i take very small values, the densities of the pushforward barycenter T i #ν for 1 ≤ i ≤ 5 exhibit large peaks (between 0.4 and 0.9) whose amplitude is beyond the largest values in the original data set (between 0.08 and 0.12).

A higher Wasserstein reconstruction error Finally, relaxing the geodesicity constraint (2.6) may actually increase the Wasserstein reconstruction error with respect to the Wasserstein distance. To state this issue more clearly, we define the reconstruction error of log-PCA as

r( Ũ) = 1 n n i=1 d 2 W ν i , exp ν (Π Sp( Ũ ) ω i ) . (3.5)
and the reconstruction error of GPCA as

r(U * ) = 1 n n i=1 d 2 W ν i , exp ν (Π Sp(U * )∩Vν (Ω) ω i ) . (3.6)
where U * is a minimizer of (2.6). Note that in (3.5), the projected measures exp ν (Π Sp( Ũ ) ω i ) might have a support that lie outside Ω. Hence, the Wasserstein distance d W in (3.5) has to be understood for measures supported on R (with the obvious extension to zero of ν i outside Ω).

The Wasserstein reconstruction error r( Ũ) of log-PCA is the sum of the Wasserstein distances of each data point ν i to a point on the surface exp ν (Sp( Ũ)) which is given by the decomposition of ω i on the orthonormal basis Ũ. However, by Proposition 2.1, the isometry property (P1) only holds between W 2 (R) and the convex subset V ν ⊂ L 2 ν (R). Therefore, we may not have

d 2 W ν i , exp ν Π Sp( Ũ ) ω i = ∥ω i -Π Sp( Ũ ) ω i ∥ 2
ν as Π Sp( Ũ ) ω i is a function belonging to L 2 ν (R) which may not necessarily be in V ν . In this case, the minimal Wasserstein distance between ν i and the surface exp ν (Sp(U * )) is not equal to ∥ω i -Π Sp(U ) ω i ∥ ν , and this leads to situations where r( Ũ) > r(U * ) as illustrated in Figure 3. 4 Two algorithmic approaches for GPCA in W 2 (Ω), for Ω ⊂ R

In this section, we introduce two algorithms which solve some of the issues of log-PCA that have been raised in Section 3. First, the output of the proposed algorithms guarantees that the computation of mappings to pushforward the Wassertein barycenter to approximate elements in the data set are strictly increasing (that is they are optimal). As a consequence, the resulting pushforward density behaves numerically much better. Secondly, the geodesic curve or surface are constrained to lie in W 2 (Ω), implying that the projections of the data are distributions whose supports do not lie outside Ω.

Iterative geodesic approach

In this section, we propose an algorithm to solve a variant of the convex-constrained GPCA problem (2.6). Rather than looking for a geodesic subset of a given dimension which fits well the data, we find iteratively orthogonal principal geodesics (i.e. geodesic set of dimension one).

Assuming that that we already know a subset

U k-1 ⊂ L 2 ν (Ω) containing k -1 orthogonal principal directions {u l } k-1 l=1 (with U 0 = ∅), our goal is to find a new direction u k ∈ L 2 ν (Ω) of
principal variation by solving the optimization problem:

u k ∈ argmin v⊥U k-1 1 n n i=1 ∥ω i -Π Sp(v)∩Vν (Ω) ω i ∥ 2 ν , (4.1)
where the infimum above is taken over all v ∈ L 2 ν (Ω) belonging to the orthogonal of U k-1 . This iterative process is not equivalent to the GPCA problem (2.6), with the exception of the first principal geodesic (k = 1). Nevertheless, it computes principal subsets U k of dimension k such that the projections of the data onto every direction of principal variation lie in the convex set V ν .

The following proposition is the key result to derive an algorithm to solve (4.1) on real data.

Proposition 4.1. Introducing the characteristic function of the convex set V ν (Ω) as:

χ Vν (Ω) (v) = 0 if v ∈ V ν (Ω) +∞ otherwise the optimization problem (4.1) is equivalent to u k = argmin v⊥U k-1 min t 0 ∈[-1;1] H(t 0 , v), (4.2)
where

H(t 0 , v) = 1 n n i=1 min t i ∈[-1;1] ∥ω i -(t 0 + t i )v∥ 2 ν + χ Vν (Ω) ((t 0 -1)v) + χ Vν (Ω) ((t 0 + 1)v). (4.3)
Proof. We first observe that Π Sp(u)∩Vν (Ω) ω i = β i u, with β i ∈ R and β i u ∈ V ν (Ω). Hence, for u k solution of (4.1), we have:

1 n n i=1 ∥ω i -Π Sp(u k )∩Vν (Ω) ω i ∥ 2 ν = 1 n n i=1 ∥ω i -β i u k ∥ 2 ν .
such that β i ∈ R and β i u k ∈ V ν (Ω) for all i ∈ {1, . . . , n}. We take M ∈ argmax 1≤i≤n β i and m ∈ argmin 1≤i≤n β i . Without loss of generality, we can assume that β M > 0 and β m < 0. We

then define v = (β M -β m )u k /2 and t 0 = (β M + β m )/(β M -β m ), that checks |t 0 | < 1.
Hence, for all i = 1, . . . , n, there exists t i ∈ [-1; 1] such that:

β i u k = (t 0 + t i )v ∈ V ν . In particular, one has t M = 1 and t m = -1, which means that (t 0 ± 1)v ∈ V ν (Ω). Reciprocally, (t 0 ± 1)v ∈ V ν (Ω) ensures us by convexity of V ν (Ω) that for all t i ∈ [-1; 1], (t 0 + t i )v ∈ V ν (Ω).
Proposition 4.1 may be interpreted as follows. For a given

t 0 ∈ [-1; 1], let v ∈⊥ U k-1 satisfying (t 0 -1)v ∈ V ν and (t 0 + 1)v ∈ V ν . Then, if one defines the curve g t (t 0 , v) = (id + (t 0 + t)v)#ν for t ∈ [-1; 1], (4.4)
it follows, from Lemma 2.1, that (g t (t 0 , v)) t∈[-1;1] is a geodesic since it can be written as g t (t 0 , v) = exp ν ((1 -u)w 0 + uw 1 ), u ∈ [0, 1] with w 0 = (t 0 -1)v, w 1 = (t 0 + 1)v, u = (t + 1)/2, and with w 0 and w 1 belonging to V ν for |t 0 | < 1. From the isometry property (P1) in Proposition 2.1, one has min

t i ∈[-1;1] ∥ω i -(t 0 + t i )v∥ 2 ν = min t i ∈[-1;1] d 2 W (ν i , g t i (v)), (4.5)
and thus the objective function H(t 0 , v) in (4.2) is equal to the sum of the squared Wasserstein distances between the data set and the geodesic curve (g t (t 0 , v)) t∈[-1;1] . The choice of the parameter t 0 corresponds to the location of the mid-point of the geodesic g t (t 0 , v), and it plays a crucial role. Indeed, the minimization of H(t 0 , v) over t 0 ∈ [-1; 1] in (4.2) cannot be avoided to obtain an optimal Wasserstein reconstruction error. This is illustrated by the Figure 3, where the Wasserstein reconstruction error r( Ũ) of log-PCA (see equation (3.5)) is compared with the ones of GPCA, for different t 0 , obtained for k = 1 as

t 0 ∈ [-1; 1] → H(t 0 , u t 0 1 ) with u t 0 1 = argmin v H(t 0 , v)
. This shows that GPCA can lead to a better low dimensional data representation than log-PCA in term of Wasserstein residual errors.

Geodesic surface approach

Once a family of vectors (v 1 , • • • , v k ) has been found through the minimization of problem (4.1), one can recover a geodesic subset of dimension k by considering all convex combinations of the vectors ((t

1 0 + 1)v 1 , (t 1 0 -1)v 1 , • • • , (t k 0 + 1)v k , (t k 0 -1)v k ).
However, this subset may not be a solution of (2.6) since we have no guarantee that a data point ν i is actually close to this geodesic subset. This discussion suggests that we may consider solving the GPCA problem (2.6) over geodesic set parameterized as in Proposition 4.1. In order to find principal geodesic subsets which are close to the data set, we consider a family

V K = (v 1 , • • • , v K ) of linearly independant vectors and t K 0 = (t 1 0 , • • • , t K 0 ) ∈ [-1, 1] K such that (t 1 0 -1)v 1 , (t 1 0 + 1)v 1 , • • • , (t K 0 -1)v K , (t K 0 + 1)v K are all in V ν .
Convex combinations of the latter family provide a parameterization of a geodesic set of dimension K by taking the exponential map exp ν of

Vν (V K , t K 0 ) = { K k=1 (α + k (t k 0 + 1) + α - k (t k 0 -1))v k , α ± ∈ A} (4.6)
where A is a simplex constraint:

α ± ∈ A ⇔ α + k , α - k ≥ 0 and K k=1 (α + k + α - k ) ≤ 1.
We hence substitute the general sets Sp(U) ∩ V ν (Ω) in the definition of the GPCA problem (2.6) to obtain,

(u 1 , • • • , u K ) = argmin V K ,t K 0 1 n n i=1 ∥ω i -Π Vν (V K ,t K ) ω i ∥ 2 ν , = argmin v 1 ,••• ,v K min t K 0 ∈[-1,1] K 1 n n i=1 min α ± i ∈A ∥ω i - K k=1 (α + ik (t k 0 + 1) + α - ik (t k 0 -1))v k ∥ 2 ν(4.7) + K k=1 χ Vν (Ω) ((t k 0 + 1)v k ) + χ Vν (Ω) ((t k 0 -1)v k ) + n i=1 χ A (α ± i ).

Discretization and Optimization

In this section we follow the framework of the iterative geodesic algorithm. We provide additional details when the optimization procedure of the geodesic surface approach differs from the iterative one.

Discrete optimization problem

Let Ω = [a; b] be a compact interval, and consider its discretization over

N points a = x 1 < x 2 < • • • < x N = b, ∆ j = x j+1 -x j , j = 1, . . . , N -1.
We recall that the functions ω i = log ν (ν i ) for 1 ≤ i ≤ n are elements of L 2 ν (Ω) which correspond to the mapping of the data to the tangent space at the Wasserstein barycenter ν. In what follows, for each 1 ≤ i ≤ n, the discretization of the function ω i over the grid reads w i = (w j i ) N j=1 ∈ R N . We also recall that χ A (u) is the characteristic function of a given set A, namely χ A (u) = 0 if u ∈ A and +∞ otherwise. Finally, the space R N is understood to be endowed with the following inner product and norm ⟨u, v⟩ ν = N j=1 f (x j )u j v j and ∥v∥ 2 ν = ⟨v, v⟩ ν for u, v ∈ R N , where f denotes the density of the measure ν. Let us now suppose that we have already computed k -1 orthogonal (in the sense ⟨u, v⟩ ν = 0) vectors u 1 , • • • u k-1 in R N which stand for the discretization of orthonormal functions u 1 , . . . , u k-1 in L 2 ν (Ω) over the grid (x j ) N j=1 . Discretizing problem (4.2) for a fixed t 0 ∈] -1; 1[, our goal is to find a new direction u k ∈ R N of principal variations by solving the following problem over all v = {v j } N j=1 ∈ R N :

u k ∈ argmin v∈R N 1 n n i=1 min t i ∈[-1;1] ∥w i -(t 0 + t i )v∥ 2 ν +χ S (v)+χ V ((t 0 -1)v)+χ V ((t 0 +1)v), (4.8)
where S = {v ∈ R N s.t. ⟨v, u l ⟩ ν = 0, l = 1 • • • k -1} is a convex set that deals with the orthogonality constraint v ⊥ U k-1 and V corresponds to the discretization of the constraints contained in V ν (Ω). From Proposition 2.1 (P3), we have that ∀v ∈ V ν (Ω), T := id + v is non decreasing and T (x) ∈ Ω for all x ∈ Ω. Hence the discrete convex set V is defined as

V = {v ∈ R N s.t. x j+1 + v j+1 ≥ x j + v j , j = 1 • • • N -1 and x j + v j ∈ [a; b], j = 1 • • • N }
and can be rewritten as the intersection of two convex sets dealing with each constraint separately.

Proposition 4.2. One has

χ V ((t 0 -1)v) + χ V ((t 0 + 1)v) = χ D (v) + χ E (Kv),
where the convex sets D and E respectively deal with the domain constraints

x j +(t 0 +1)v j ∈ [a; b] and x j + (t 0 -1)v j ∈ [a; b], i.e.: D = {v ∈ R N , s.t. m j ≤ v j ≤ M j }, (4.9) with m j = max a-x j t 0 +1 , b-x j t 0 -1 and M j = min a-x j t 0 -1 , b-x j t 0 +1
, and the non decreasing constraint of id + (t 0 ± 1)v:

E = {z ∈ R N s.t. -1/(t 0 + 1) ≤ z j ≤ 1/(1 -t 0 )}.
(4.10) with the differential operator K : R N → R N computing the discrete derivative of v ∈ R n as

(Kv) j = (v j+1 -v j )/(x j+1 -x j ) if 1 ≤ j < N 0 if j = N, (4.11)
Having D and E both depending on t 0 is not an issue since problem (4.8) is solved for fixed t 0 .

Introducing t = {t i } n i=1 ∈ R n , problem (4.8) can be reformulated as:

min v∈R N min t∈R n J(v, t) := n i=1 ∥w i -(t 0 + t i )v∥ 2 ν F (v,t) + χ S (v) + χ D (v) + χ E (Kv) + χ B n 1 (t) G(v,t)
. (4.12)

where B n 1 is the L ∞ ball of R n with radius 1 dealing with the constraint t i ∈ [-1; 1]. Notice that F is differentiable but non-convex in (v, t) and G is non-smooth and convex.

Geodesic surface approach For fixed (t 1 0 , . . . ,

t K 0 ) ∈ R K and ff ± = {α + k , α - k } K k=1 , the discretized version of (4.7) is then min v 1 ,...,v K ∈R N min ff ± 1 ,...,ff ± n ∈R 2K F ′ (v, t) + G ′ (v, t), (4.13)
where

F ′ (v, t) = n i=1 ∥w i -K k=1 (α + ik (t k 0 + 1) + α - ik (t k 0 -1))v k ∥ 2 ν is still non-convex and differ- entiable, G ′ (v, t) = K k=1 χ E (Kv k ) + χ D k (v k ) + n i=1 χ A (ff ± i )
2 is convex and non smooth, A is the simplex of R 2K and D k is defined as in (4.9), depending on t k 0 . We recall that the orthogonality between vectors v k is not taken into account in the geodesic surface approach.

Optimization through the Forward-Backward Algorithm

Following [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems[END_REF], in order to compute a critical point of problem (4.12), one can consider the Forward-Backward algorithm (see also [START_REF] Ochs | ipiano: Inertial proximal algorithm for nonconvex optimization[END_REF] for an acceleration using inertial terms). Denoting as X = (v, t) ∈ R N +n , taking τ > 0 and X (0) ∈ R N +n , it reads:

X (ℓ+1) = Prox τ G (X (ℓ) -τ ∇F (X (ℓ) )), (4.14) where Prox τ G ( X) = argmin X 1 2τ ||X -X|| 2 + G(X) with the Euclidean norm || • ||.
In order to guarantee the convergence of this algorithm, the gradient of F has to be Lipschitz continuous with parameter M > 0 and the time step should be taken as τ < 1/M . The details of computation of ∇F and Prox τ G for the two algorithms are given in Appendix A.

5 Statistical comparison between log-PCA and GPCA on synthetic and real data

Synthetic example -Iterative versus geodesic surface approaches

First, for the synthetic example displayed in Figure 1, we compare the two algorithms (iterative and geodesic surface approaches) described in Section 4. The results are reported in Figure 4 by comparing the projection of the data onto the first and second geodesics computed with each approach. We also display the projection of the data onto the two-dimensional surface generated by each method. It should be recalled that the principal surface for the iterative geodesic algorithm is not necessarily a geodesic surface but each g t (t k 0 , u k ) t∈[-1;1] defined by (4.4) for k = 1, 2 is a geodesic curve for U = {u 1 , u 2 }. For data generated from a location-scale family of Gaussian distributions, it appears that each algorithm provides a satisfactory reconstruction of the data set. The main divergence concerns the first and second principal geodesic. Indeed enforcing the orthogonality between components in the iterative approach enables to clearly separate the modes of variation in location and scaling, whereas searching directly a geodesic surface in the second algorithm implies a mixing of these two types of variation.

Note that the barycenter of Gaussian distributions N (m i , σ 2 i ) can be shown to be Gaussian with mean m i and variance ( σ i ) 2 .

Population pyramids

As a first real example, we consider a real dataset whose elements are histograms representing the population pyramids of n = 217 countries for the year 2000 (this dataset is produced by the International Programs Center, US Census Bureau (IPC, 2000), available at https: //www.census.gov/programs-surveys/international-programs.html). Each histogram in the database represents the relative frequency by age, of people living in a given country. Each bin in a histogram is an interval of one year, and the last interval corresponds to people older than 85 years. The histograms are normalized so that their area is equal to one, and thus they represent a set of pdf. In Figure 5, we display the population pyramids of 4 countries, and the whole dataset. Along the interval Ω = [0, 84], the variability in this dataset can be considered as being small. For K = 2, log-PCA and the iterative GPCA algorithm lead to the same principal orthogonal directions in L 2 ν (Ω), namely that ũ1 = u * 1 and ũ2 = u * 2 where (ũ 1 , ũ2 ) minimizes (3.1) and (u * 1 , u * 2 ) are minimizers of (4.2). In this case, all projections of data ω i = log ν (ν i ) for i = 1, . . . , n onto Sp({ũ 1 , ũ2 }) lie in V ν (Ω), which means that log-PCA and the iterative geodesic algorithm lead exactly the same principal geodesics. Therefore, population pyramids is an example of data that are sufficiently concentrated around their Wasserstein barycenter so that log-PCA and GPCA are equivalent approaches (see Remark 3.5 in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF] for further details). Hence, we only display in Figure 6 the results of the iterative and geodesic surface algorithms.

In the iterative case, the projection onto the first geodesic exhibits the difference between less developed countries (where the population is mostly young) and more developed countries (with an older population structure). The second geodesic captures more subtle divergences concentrated on the middle age population. It can be observed that the geodesic surface algorithm gives different results since the orthogonality constraint on the two principal geodesics is not required. In particular, the principal surface mainly exhibit differences between countries with a young population with countries having an older population structure, but the difference between its first and second principal geodesic is less contrasted. The first row is the GPCA of the data set obtained with the iterative geodesic approach. The second row is the GPCA through the geodesic surface approach. The first (resp. second) column is the projection of the data into the first (resp. second) principal direction. The black curve is the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

Children's first name at birth

In a second example, we consider a dataset of histograms which represent, for a list of n = 1060 first names, the distribution of children born with that name per year in France between years 1900 and 2013. In Figure 7, we display the histograms of four different names, as well as the whole dataset. Along the interval Ω = [1900,2013], the variability in this dataset is much larger than the one observed for population pyramids. This dataset has been provided by the INSEE (French Institute of Statistics and Economic Studies) [INSur].

This is an example of real data where log-PCA and GPCA are not equivalent procedures for K = 2 principal components. We recall that log-PCA leads to the computation of principal orthogonal directions ũ1 , ũ2 in L 2 ν (Ω) minimizing (3.1). First observe that in the left column of Figure 8, for some data ω i = log ν (ν i ), the mappings Ti = id + Π Sp({ũ 1 }) ω i are decreasing, and their range is larger than the interval Ω (that is, for some x ∈ Ω, one has that Ti (x) / ∈ Ω). Hence, such Ti are not optimal mappings. Therefore, the condition Π Sp( Ũ ) ω i ∈ V ν (Ω) for all 1 ≤ i ≤ n (with Ũ = {ũ 1 , ũ2 }) is not satisfied, implying that log-PCA does not lead to a solution of GPCA thanks to Proposition 3.5 in [START_REF] Bigot | Geodesic PCA in the Wasserstein space by Convex PCA[END_REF].

Hence, for log-PCA, the corresponding histograms displayed in the right column of Figure 8 are such that Π Sp({ũ 1 }) ω i / ∈ V ν (Ω). This implies that the densities of the projected measures The dashed red curves represent the mapping Ti = id + Π Sp({ũ1}) ω i where ω i = log ν (ν i ), and ũ1 is the first principal direction in L 2 ν (Ω) obtained via log-PCA. The blue curves are the mapping T i = id + Π Sp({u * 1 }) ω i , where u * 1 is the first principal direction in L 2 ν (Ω) obtained via the iterative algorithm. (Right) The histogram stands for the pdf of measures ν i that have a large Wasserstein distance with respect to the barycenter ν. The red curves are the pdf of the projection exp ν (Π Sp(ũ1) ω i ) with log-PCA, while the blue curves are the pdf of the projection exp ν (Π Sp(u * 1 ) ω i ) with GPCA.

exp ν (Π Sp(ũ 1 ) ω i ) have a support outside Ω = [1900,2013]. Hence, the estimation of the measure ν i = exp ν (ω i ) by its projection onto the first mode of variation obtained with log-PCA is not satisfactory.

In Figure 8, we also display the results given by the iterative geodesic algorithm, leading to orthogonal directions u * 1 , u * 2 in L 2 ν (Ω) that are minimizers of (4.2). Contrary to the results obtained with log-PCA, one observes in Figure 8 that all the mapping T i = id + Π Sp({u * 1 }) ω i are non-decreasing, and such that T i (x) ∈ Ω for all x ∈ Ω. Nevertheless, by enforcing these two conditions, one has that a good estimation of the measure ν i = exp ν (ω i ) by its projection exp ν (Π Sp(u * 1 ) ω i ) is made difficult as most of the mass of ν i is located at either the right or left side of the interval Ω which is not the case for its projection. The histograms displayed in the right column of Figure 8 correspond to the elements in the dataset that have a large Wasserstein distance with respect to the barycenter ν. This explains why it is difficult to have good projected measures with GPCA. For elements in the dataset that are closest to ν, the projected measures exp ν (Π Sp(ũ 1 ) ω i ) and exp ν (Π Sp(u * 1 ) ω i ) are much closer to ν i and for such elements, log-PCA and the iterative geodesic algorithm lead to similar results in terms of data projection.

To better estimate the extremal data in Figure 8, a solution is to increase the support of the data to the interval Ω 0 = [1850, 2050], and to perform log-PCA and GPCA in the Wasserstein space W 2 (Ω 0 ). The results are reported in Figure 9. In that case, it can be observed that both algorithms lead to similar results, and that a better projection is obtained for the extremal data. Notice that with this extended support, all the mappings Ti = id + Π Sp({ũ 1 }) ω i obtained with log-PCA are optimal in the sense that they are non-decreasing with a range inside Ω 0 .

Finally, we display in Figure 10 and Figure 11 the results of the iterative and geodesic surface algorithms with either Ω = [1900,2013] or with data supported on the extended support Ω 0 = [1850,2050]. The projection of the data onto the first principal geodesic suggests that the distribution of a name is deeply dependent on the part of the century. The second geodesic expresses a popular trend through a spike effect. In Figure 10, the artefacts in the principal surface that are obtained with the iterative algorithm at the end of the century, correspond to the fact that the projection of the data ω i onto the surface spanned by the first two components is not ensured to belong to the set V ν (Ω).

Extensions beyond d > 1 and some perspectives

We now briefly show that our iterative algorithm for finding principal geodesics can be adapted to the general case d > 1. This requires to take into account two differences with the onedimensional case. First, the definition of the space V µr (Ω) in 2.1 relies on the explicit close-form solution (2.2) of the optimal transport problem which is specific to the one-dimensional case. We must hence provide a more general definition of V µr (Ω). Second, the isometry property (P 1) does not hold for d > 1, so that Wasserstein distances cannot be replaced by the L 2 ν norm between log-maps as in (4.5) and must be explicitly computed and differentiated.

Definition of V µr (Ω) in the general case In the one dimensional case, V µr (Ω) is characterized in Proposition 2.1 (P3) as the set of functions v ∈ L 2 µr (Ω) such that T := id + v is µ r -almost everywhere non decreasing. A important result by Brenier [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] is that, in any dimension, if µ r does not give mass to small sets, there exists an optimal mapping T ∈ L 2 µr (Ω) between µ r and 22 any probability measure ν, and T is equal to the gradient of a convex function u, i.e. T = ∇u. Therefore we define the set V µr (Ω) as the set of functions v ∈ L 2 µr (Ω) such that id + v = ∇u for an arbitrary convex function u.

In order to deal with the latter constraint, we note it implies that div(v) ≥ -1. Indeed, assuming that id + v = ∇u, then u being a convex potential involves div(∇u) ≥ 0 which is equivalent to div(id + v) = div(v) + 1 ≥ 0. We therefore choose to substitute this constraint by the constraint div(v) ≥ -1.

General objective function Without the isometry property (P1), the objective function H(t 0 , v) in (4.3) must be written with the explicit Wasserstein distance d W ,

H(t 0 , v) = 1 n n i=1 min t i ∈[-1;1] d 2 W (ν i , g t i (t 0 , v)) + χ Vν (Ω) ((t 0 -1)v) + χ Vν (Ω) ((t 0 + 1)v), ( 6
d 2 W (ν i , g t i (t 0 , v)) F (v,t) + χ S (v) + χ E (Kv) + χ D (v) + χ B n 1 (t) G(v,t) . (6.2)
where K is a discretized divergence operator, and E = {z ∈ R N : -1 t 0 +1 ⩽ z ⩽ 1 1-t 0 }, D = {v : id+(t 0 ±1)v ∈ Ω} deals with the domain constraint and S deals with the orthogonality constraint w.r.t. to the preceding principal components. As for the one-dimensional case, we minimize J through the Forward-Backward algorithmas detailed in the appendix B.

Extension to higher dimensions is straightforward. However, considering that we have to discretize the support of the Wasserstein mean ν, the approach becomes intractable for d > 3.

Application to grayscale images

We consider the MNIST dataset [START_REF] Lecun | The mnist database of handwritten digits[END_REF] which contains grayscale images of handwritten digits. All the images have identical size 28 × 28 pixels. Each grayscale image, once normalized so that the sum of pixel grayscale values sum to one, can be interpreted as a discrete probability measure, which is supported on the 2D grid of size 28 × 28. The ground metric for the Wasserstein distance is then the 2D squared Euclidean distance between the locations of the pixels of the two-dimensional grid. We compute the first principal components on 1000 images of each digit. Wasserstein barycenters, which are required as input to our algorithm, are approximated The second row is the GPCA through the geodesic surface approach. The first (resp. second) column is the projection of the data into the first (resp. second) principal direction. The black curve is the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

efficiently through iterative Bregman projections as proposed in [BCC + 15]. We use the network simplex algorithm2 to compute Wasserstein distances. Figure 12 displays the results obtained with our proposed forward-backward algorithm (with t 0 set to 0 for simplicity), and the ones given by Log-PCA as described in section 3. These two figures are obtained by sampling the first principal components. We then use kernel smoothing to display the discrete probability measures back to the original grid and present the resulting grayscale image with an appropriate colormap.

Visually, both the Log-PCA and GPCA approaches capture well the main source of variability of each set of grayscale images. We observe variations in the slant of the handwritten digits for all digits, the most obvious case being digit '1'. As a principal component is parameterized by a whole velocity field on the support of the Wasserstein mean of the data, single principal components can capture more interesting patterns, such as changes in the shape of the '0' or the presence or absence of the lower loop of the '2'. From purely visual inspection, it is difficult to tell which approach, Log-PCA or GPCA, provides a "better" principal component. For this purpose we compute the reconstruction error of each digit. This reconstruction error is computed in the same way for both Log-PCA and GPCA principal components: We sample the principal components at many times t and find for each image in a given dataset, the time at which the geodesic is the closest to the image sample. This provides an approximation of GPCA log-PCA min t∈[-1;1] d 2 W (ν i , g t (v)) for each image i = 1, . . . , n, where (g t ) t∈[-1,1] is the principal component. For the Log-PCA principal component, we take gt = (id + t1.25λv)#ν, where λ is the eigenvalue corresponding to the first principal component. The 1.25 factor is useful to consider a principal curve which goes through the whole range of the dataset. For the GPCA principal geodesic, we have g * t = (id + tv)#ν. The reconstruction errors are shown in Table 1. We see that, for each digit, we obtain a better, i.e. smaller, reconstruction error when using the proposed forward-backward algorithm. This result is not surprising, since the reconstruction error is explicitly minimized through the Forward-Backward algorithm. As previously mentioned, Log-PCA rather computes linearized Wasserstein distances. In one-dimension, the isometry property (P1) states that these quantities are equal. In dimension two or larger, that property does not hold.

Discussion

The proposed forward-backward algorithm minimizes the same objective function as defined in [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF]. The first difference with the algorithm provided [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] is that we take gradient steps with respect to both v and t, while the latter first attempts to find the optimal t (by sampling the geodesics at many time t), before taking a gradient step of v. Our approach reduces the cost of computing a gradient step by one order of magnitude. Secondly, [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] relied on barycentric projections of optimal plans to preserve the geodesicity of the principal curves in between gradient steps. That heuristic does not guarantee a decrease in the objective after a gradient step. Moreover, the method in [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] considered two velocity fields v 1 , v 2 rather than a single v since the optimality of both v and -v could not be preserved through the barycentric projection. When considering probability measures over high dimensional space (d > 3), our algorithm becomes intractable since we need to discretize the support of the Wasserstein mean of the data with a regular grid, while the approach of [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] is still tractable since an arbitrary support for the Wasserstein mean is used. A remaining challenge for computing principal geodesics in the Wasserstein space is then to propose an algorithm for GPCA which is still tractable in higher dimensions while not relying on barycentric projections of optimal transport plans.

Acknowledgement : The authors would like to thank the INSEE for providing the first name at birth dataset [INSur].

A Dimension d = 1

We here detail the application of Algorithm (4.14) to the iterative GPCA procedure that consists in solving the problem (4.12):

min v∈R N min t∈R n J(v, t) := n i=1 ∥w i -(t 0 + t i )v∥ 2 ν F (v,t) + χ S (v) + χ D (v) + χ E (Kv) + χ B n 1 (t) G(v,t)
.

A.1 Lispschitz constant of ∇F

Let us now look at the Lipschitz constant of ∇F (v, t) on the restricted acceptable set D × B n 1 . We first denote as H the hessian matrix (of size (N + n) × (N + n)) of the C 2 function F (X). We know that if the spectral radius of H is bounded by a scalar value M , i.e. ρ(H) ≤ M , then ∇F is a Lipschitz continuous function with constant M . Hence, we look at the eigenvalues of the Hessian matrix of F = n i=1 N j=1 f n (x j )(w j i -(t 0 + t i )v j ) 2 that is

∂ 2 F ∂t 2 i = N j=1 2v 2 j f n (x j ), ∂ 2 F ∂v 2 j = n i=1
2(t 0 + t i ) 2 f n (x j ), ∂ 2 F ∂t i ∂v j = 2 f n (x j )(2(t 0 + t i )v j -w j i ) and ∂ 2 F ∂t i ∂t i ′ = ∂ 2 F ∂v j ∂v j ′ = 0, for all i ̸ = i ′ or j ̸ = j ′ . Being This problem can be solved independently on v and t. For t, it can be done pointwise as where z = {z j } N j=1 ∈ R N is a dual variable and χ * E = sup v ⟨v, z⟩ -χ E (v) is the convex conjugate of χ E that reads: (χ * E (z)) j = -z j /(1 + t 0 ) if z j ≤ 0, z j /(1 -t 0 ) if z j > 0.

t * i = argmin t i 1 2τ ||t i -ti || 2 + χ B 1 1 (t i ) = Proj
everywhere. Gradients can be computed in the same fashion as [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF] to obtain, where z = {z j } N j=1 ∈ R N is a dual variable and χ * E = sup v ⟨v, z⟩ -χ E (v) is the convex conjugate of χ E . This can be solve with the same iterative steps as described in A.2,

∇ v F = 2 n i=1 (t 0 +t i )(Z t i -XP * T diag(1/ f n )), ∇ t i F = 2⟨Z t i diag( f n ), v ⟩-2⟨P * , v T X
   z (m+1) = Prox σχ * E (z (m) + σK v(m) ) v (m+1) = Prox θ(χ D +χ S ) (v (m) -θ(K * z (m+1) + 1 τ (v (m) -ṽ)) v(m+1) = 2v (m+1) -v (m) (B.4)
Here the definition of the divergence operator K and the transpose of the divergence operator K * are specific to the dimension. For d = 2, with a regular grid discretizing Ω in M × N points, we take

K T z = -∇z = - ∂ + x z ∂ + y z
, with ∂ + x z(i, j) = z(i + 1, j) -z(i, j) if i < M 0 otherwise, ∂ + y z(i, j) = z(i, j + 1) -z(i, j) if j < N 0 otherwise so that

Ku = K u x u y = ∂ - x u x + ∂ - y u y , with ∂ - x u(i, j) =    u(i, j) -u(i -1, j) if 1 < i < M u(i, j) if i = 1 -u(i -1, j) if i = M.
To ensure convergence of B.4, one can take 1/σ.(1/θ -1/τ ) = ||K|| 2 . See [START_REF] Chambolle | On the ergodic convergence rates of a first-order primal, Äìdual algorithm[END_REF][START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF] for more details. Since we have ||K|| 2 = 8, the parameters can be taken as σ = 1/4 and θ = τ /(1 + 2τ ).

Figure 1 :

 1 Figure 1: Synthetic example. (Right) A data set of n = 100 Gaussian histograms randomly translated and scaled. (Top-left) Standard PCA of this data set with respect to the Euclidean metric. The Euclidean barycenter of the data set is depicted in blue. (bottom-left) Geodesic PCA with respect to the Wasserstein metric using the iterative geodesic algorithm (4.1). The black curve represents the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

Figure 2 :

 2 Figure 2: (Left) Distribution of the total precipitation (mm) collected in a year in 1 ≤ i ≤ 5 stations among 60 in China -Source : Climate Data Bases of the People's Republic of China 1841-1988 downloaded from http://cdiac.ornl.gov/ndps/tr055.html. The black curve is the density of the Wasserstein barycenter of the 60 stations. (Middle)Mapping T i = id+Π Sp(ũ2) ω i obtained from the projections of these 5 distributions onto the second eigenvector ũ2 given by log-PCA of the whole dataset. (Right) Pushforward exp ν (Π Sp(ũ2) ω i ) = T i #ν of the Wasserstein barycenter ν for each 1 ≤ i ≤ 5. As the derivative T ′ i take very small values, the densities of the pushforward barycenter T i #ν for 1 ≤ i ≤ 5 exhibit large peaks (between 0.4 and 0.9) whose amplitude is beyond the largest values in the original data set (between 0.08 and 0.12).

Figure 3 :

 3 Figure 3: Comparison of the Wasserstein reconstruction error between GPCA and log-PCA on the synthetic dataset displayed in Figure 1 for the first component, with an illustration of the role of the parameter t 0 in (4.2).

Figure 4 :

 4 Figure 4: Synthetic example -Data sampled from a location-scale family of Gaussian distributions. The first row is the GPCA of the data set obtained with the iterative geodesic approach. The second row is the GPCA through the geodesic surface approach. The black curve is the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

Figure 5 :

 5 Figure 5: Population pyramids. A subset of population pyramids for 4 countries (left) for the year 2000, and the whole dataset of n = 217 population pyramids (right) displayed as pdf over the interval [0, 84].

Figure 6 :

 6 Figure6: Population pyramids. The first row is the GPCA of the data set obtained with the iterative geodesic approach. The second row is the GPCA through the geodesic surface approach. The first (resp. second) column is the projection of the data into the first (resp. second) principal direction. The black curve is the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

ROBERTFigure 7 :Figure 8 :

 78 Figure 7: Children's first name at birth. An subet of 4 histograms representing the distribution of children born with that name per year in France, and the whole dataset of n = 1060 histograms (right), displayed as pdf over the interval[1900, 2013]

Figure 9 :

 9 Figure 9: Children's first name at birth with extended support Ω 0 = [1850, 2050].(Left) The dashed red curves represent the mapping Ti = id + Π Sp({ũ1}) ω i where ω i = log ν (ν i ), and ũ1 is the first principal direction in L 2 ν (Ω) obtained via log-PCA. The blue curves are the mapping T i = id + Π Sp({u * 1 }) ω i , where u * 1 is the first principal direction in L 2 ν (Ω) obtained via the iterative algorithm. (Right) The histogram stands for the pdf of measures ν i that have a large Wasserstein distance with respect to the barycenter ν. The red curves are the pdf of the projection exp ν (Π Sp(ũ1) ω i ) with log-PCA, while the blue curves are the pdf of the projection exp ν (Π Sp(u * 1 ) ω i ) with GPCA.

Figure 11 :

 11 Figure 11: Children's first name at birth with extended support Ω 0 = [1850, 2050]. The first row is the GPCA of the data set obtained with the iterative geodesic approach.The second row is the GPCA through the geodesic surface approach. The first (resp. second) column is the projection of the data into the first (resp. second) principal direction. The black curve is the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

Figure 12 :

 12 Figure 12: First principal geodesics for 1000 images of each digit from the MNIST dataset, computed through the proposed Forward-Backward algorithm (left) and log-PCA (right).

  {µ k } n+N k=1 the eigenvalues of H, we have ρ(H) = max k |µ k | ≤ max k l |H kl |. We denote as f ∞ = max j | f n (x j )| and likewise w ∞ = max i,j |w j i |. Since |t 0 | < 1, t 2 i ≤ 1, ∀t ∈ B n 1 and v 2 j ≤ α 2 = (b -a) 2 , ∀v ∈ D, by defining γ = 2(1 + |t 0 |)α + w ∞ , we thus have ρ(H) ≤ 2f ∞ max nα 2 + N γ, nγ + N (1 + |t 0 |) 2 := M.(A.1)A.2 Computing Prox τ GIn order to implement the algorithm (4.14), we finally need to compute the proximity operator of G defined as:(v * , t * ) = Prox τ G (ṽ, t) = argmin v,t 1 2τ (||v -ṽ|| 2 + ||t -t|| 2 ) + χ S (v) + χ D (v) + χ E (Kv) + χ B n 1 (t).

  [-1;1] ( ti ). Unfortunately, there is no closed form expression of the proximity operator for the component v. It requires to solve the following intern optimization problem at each extern iteration (ℓ) of the algorithm (4.14):v * = argmin v 1 2τ ||v -ṽ|| 2 + χ S (v) + χ D (v) + χ E (Kv), (A.2)where, to avoid confusions, we denote by v the variable that is optimized within the intern optimization problem (A.2).Remark A.1. The Lipschitz constant of ∇F (v, t) in (A.1) relies independantly on v and |t 0 |, thus we can choose the optimal gradient descent step τ for v * and t * .Primal-Dual reformulation Using duality (through Fenchel transform), one has:min v∈R N 1 2τ ||v -ṽ|| 2 + χ S (v) + χ D (v) + χ E (Kv) 2 + χ S (v) + χ D (v) + ⟨Kv, z⟩ -χ * E (z), (A.3)

  ⟩, (B.1)Proximal operator of G The only difference between the one-dimensional case and the twodimensional case considered here concerns the projection step of v,v * = argmin v 1 2τ ||v -ṽ|| 2 + χ S (v) + χ D (v) + χ E (Kv), (B.2)Primal-Dual reformulation As for the on-dimensional case, one has, minv∈R N 1 2τ ||v -ṽ|| 2 + χ S (v) + χ D (v) + χ E (Kv) = min v∈R N max z∈R N 1 2τ ||v -ṽ|| 2 + χ S (v) + χ D (v) + ⟨Kv, z⟩ -χ * E (z), (B.3)

  Children's first name at birth with support Ω =[1900, 2013]. The first row is the GPCA of the data set obtained with the iterative geodesic approach. The second row is the GPCA through the geodesic surface approach. The first (resp. second) column is the projection of the data into the first (resp. second) principal direction. The black curve is the density of the Wasserstein barycenter. Colors encode the progression of the pdf of principal geodesic components in W 2 (Ω).

		0.06	First PG	0.06	Second PG	0.06	Principal Surface	Data set
	Iterative Geodesic Algorithm	0.01 0.02 0.03 0.04 0.05			0.01 0.02 0.03 0.04 0.05		0.01 0.02 0.03 0.04 0.05	0.2 0.25 0.3
		0			0		0
		1900 1950 2000	1900 1950 2000	1900 1950 2000
								0.15
		0.06			0.06		0.06
	Geodesic Surface Algorithm	0.01 0.02 0.03 0.04 0.05			0.01 0.02 0.03 0.04 0.05		0.01 0.02 0.03 0.04 0.05	0.05 0.1
		0			0		0	0
		1900 1950 2000	1900 1950 2000	1900 1950 2000	1900 1950 2000
	Figure 10: min v∈R N min t∈R n J(v, t) :=	n i=1
								.1)
	where g t (t 0 , v) = (id + (t 0 + t)v)#ν for t ∈ [-1; 1] as defined in (4.4). Optimizing over both the
	functions v ∈ (R d ) N and the projection times t, the discretized objective function to minimize
	is,

Table 1 :

 1 MNIST digit Log-PCA RE (•10 3 ) GPCA RE (•10 3 ) Reconstruction Errors (RE) computed on 1000 sample images of each digit of the MNIST dataset. (center) Reconstruction error w.r.t. the first principal component computed with the Log-PCA algorithm. (right) Reconstruction error w.r.t. the first principal geodesic computed with the proposed Forward-Backward algorithm.

	0	2.0355	1.9414
	1	3.1426	1.0289
	2	3.4221	3.3575
	3	2.6528	2.5869
	4	2.8792	2.8204
	5	2.9391	2.9076
	6	2.1311	1.9864
	7	4.7471	2.8205
	8	2.0741	2.0222
	9	1.9303	1.8728

https://github.com/ecazelles/2017-GPCA-vs-LogPCA-Wasserstein

http://liris.cnrs.fr/~nbonneel/FastTransport/

Hence, one can use the Primal-Dual algorithm proposed in [START_REF] Chambolle | On the ergodic convergence rates of a first-order primaldual algorithm[END_REF] to solve the problem (A.3). For two parameters σ, θ > 0 such that ||K|| 2 ≤ 1 σ ( 1 θ -1 τ ) and given v 0 , v0 , z 0 ∈ R N , the algorithm is:

where K * is defined as ⟨Kv, z⟩ = ⟨v, K * z⟩. Using the operator K defined in (4.11), we thus have:

where ∆ j = x j+1 -x j . We have that ||K|| 2 = ρ(K * K), the largest eigenvalue of K * K. With the discrete operators (4.11) and (A.5), ρ(K * K) can be bounded by

One can therefore for instance take σ = 1 δ and θ = τ /(1 + δτ ).

Proximity operators in (A.4) The proximity operator of χ D + χ S is obtained as:

since projecting onto D ∩ S is equivalent to first project onto the orthogonal of Sp(U k-1 ) and then onto D. One can finally show that the proximity operator of χ * E can be computed pointwise as:

(A.8)

A.3 Algorithms for GPCA

Gathering all the previous elements, we can finally find a critical point of the non-convex problem (4.12) using the Forward-Backward (FB) framework (4.14), as detailed in Algorithm 1.

Algorithm 1 Resolution with FB of problem (4.12): min v,t F (v, t) + G(v, t)

Require:

)) and M > 0 (defined in (A.1)). Set (v (0) , t (0) ) ∈ D × B n 1 Set τ < 1/M , σ = 1/δ and θ = τ /(1 + δτ ). %Extern loop: while ||v (ℓ) -v (ℓ-1) ||/||v (ℓ-1) || > η do % FB on t with t (ℓ+1) = P rox τ G (t (ℓ) -τ ∇F (v (ℓ) , t (ℓ) )):

Geodesic surface approach In order to solve the problem (4.13), we follow the same steps as in the section A.1-A.2. First we obtain the Lipchitz constant of the function F by the same computations performed for the iterative algorithm. Then, since the constraints' problem in G ′ are separable, we can compute each component v k and each ff ± i independantly. The only difference with the iterative algorithm concerns the proximal operator of the function χ A , which is the projection into the simplex of R 2K .

B Dimension d = 2

We now show how to generalize the algorithm to the two-dimensional case.

Gradients of F We write X = (x 1 , • • • , x N ) ∈ (R 2 ) N the discretized support of ν, Z t = (x 1 + (t 0 + t)v 1 , • • • , x N + (t 0 + t)v N ) the support g t (t 0 , v), the geodesic sampled at time t. Let P * be an optimal transport plan between ν and g t (t 0 , v). The function F (v, t) is differentiable almost