On the Wasserstein Geodesic Principal Component Analysis of probability measures
Authors: Nina Vesseron
Elsa Cazelles
Alice Le Brigant
Klein
Venue: ICLR 2026 Oral (2026)
Keywords: wasserstein PCA
optimal transport
deep learning
Primary area: unsupervised, self-supervised, semi-supervised, and supervised representation learning

Abstract: This paper focuses on Geodesic Principal Component Analysis (GPCA) on a collection of probability distributions using the Otto-Wasserstein geometry. The goal is to identify geodesic curves in the space of probability measures that best capture the modes of variation of the underlying dataset. We first address the case of a collection of Gaussian distributions, and show how to lift the computations in the space of invertible linear maps. For the more general setting of absolutely continuous probability measures, we leverage a novel approach to parameterizing geodesics in Wasserstein space with neural networks. Finally, we compare to classical tangent PCA through various examples and provide illustrations on real-world datasets.

Published as a conference paper at ICLR 2026
ON THE WASSERSTEIN GEODESIC PRINCIPAL COMPO- NENT ANALYSIS OF PROBABILITY MEASURES
Nina Vesseron CREST-ENSAE, IP Paris
Elsa Cazelles
CNRS, IRIT, Université de Toulouse
(nina. www.eressseeeeeeeoe Bensae. fr
{el caze 1eeellees les) fr
Alice Le Brigant
Thierry Klein
Université Paris Panthéon Sorbonne {alice. 1tw.:::gaatt/unnii-aaaaaaaee
ENAC, IMT, Université de Toulouse (thierry. .lein) @math. univ-toulouse. fr
ABSTRACT
This paper focuses on Geodesic Principal Component Analysis (GPCA) a collection of probability distributions using the Otto-Wassersterf geometry. The goal is to identify geodesic curves in the space of probability measures that best capture the modes variation of the underlying dataset. We first address the case of collection of Gaussian distributions, and show how to lift the computations to the space of invertible linear maps. For the more general setting of absolutely continuous probability measures, leverage novel approach parameterizing geodesics in Wasserstein space with neural networks. Finally, we compare to classical tangent PCA through various examples and provide illustrations on real- world datasets.
INTRODUCTION
In this paper, we are interested in computing the main modes of variation of dataset of absolutely continuous (a.c.) probability measures supported in RR, For data points living an arbitrary Hilbert space, the classical approach defined by Principal Component Analysis (PCA) consists in finding sequence of nested affine subspaces on which the projected data retain a maximal part of the variance of the original dataset, equivalently, yield best lower-dimensional approximations. When dealing with a set of a.c. probability distributions, a natural choice to identify the probability measures with their probability density functions and perform PCA on these using the L Hilbert metric. Unfortunately, as highlighted in Cazelles et al. (2018), the components computed in this manner fail capture the intrinsic structure of the dataset: the projections onto the components most likely result in non- non-pooitivee and un-rommale functions. Using the Wasserstein metric W2 instead has proven to overcome these limitations by taking into account the geometry of the space of distributions.
The Wasserstein metric endows the space of probability distributions with a Riemannian- structure, framing the problem as PCA on a (positively) curved Riemannian manifold. first approach to solve this task, known as Tangent PCA (TPCA), consists in embedding the data into the tangent space reference point, and applying classical PCA in this flat space, as in Fletcher al. (2003). TPCA is computationally advantageous but can generically induce distortion in the embedded data, depending on the curvature of the manifold at the reference point and the dispersion of the data. A more geometrically coherent approach is Geodesic PCA (GPCA) proposed for Riemannian manifolds in Huckemaan et al. (2010); Huckemann & Ziezold (2006), where principal modes of variations are geodesics that minimize the variance of the projection residuals. Following this approach, the first geodesic component of a set of probability measures V1, Vn in the Wasserstein space solves
inf (+++(() geodesic
Σ inf (a(
(1)
"These authors contributed equally this work.

---

Published as a conference paper at ICLR 2026
Interestingly, unlike in the Hilbert setting, this criterion is not equivalent to maximizing the variance of the projections, which leads to a different notion of PCA on Riemannian manifolds (see Sommer et al. 2014))
Related works TPCA in the Wasserstein space was considered by Wang et al. (2013) through the use of the linearized Wasserstein distance. In similar approach, [Boissrdd et al. (2013) restrict to distributions that can be obtained by deforming a single template measure. For probability measures, Bigot et al. (2017) have shown that GPCA and its linearized approximation TPCA coincide, as the embedding into a tangent space is then an isometry when constrained to convex set. An algorithm in this case has been proposed Cazelles et al. (2188), with an approximate extension in dimension 2. For higher-dimensional measures, (Seguy & Cuturi (2015) solve an approximate version of GPCA, replacing geodesics by generalized geodesics as defined in Ambrosio et al. (2008). Despite all this, method to solve the exact GPCA problem described in equation is still missing for -value probability measures. The goal of this paper is to fill this gap.
Main contributions In this paper, introduce two algorithms solve the exact GPCA problem in the Wasserstein space of (1) centered Gaussian distributions and (2) a.c. probability measures on iRd Our methods are exact in the sense that they do not rely linearization of the Wasserstein space, and the components are true geodesics that minimize the cost in equation I In the Gaussian case, we leverage the matrices. We show an example where GPCA and TPCA differ significantly, and relate this effect to curvature. In the general case of a.c. probability distributions, we lift the probability distributions to the space of (non necessarily optimal) maps that pushforward given reference measure, as described by vyto (200)). This approach independent of the chosen reference measure and yields convenient way to parametrize geodesic components and define orthogonality with respect to the Wasserstein metric. In practice, we parametrize geodesic components using multilayer perceptrons (MLPs), trained to minimize the cost in equation IT] We show illustrations images and 3D point clouds. Along the way, we prove that for univariate Gaussian distributions, GPCA yields the same results whether it is performed in the space of a.c. distributions restricted to the Gaussian submanifold.
geometry to lift the computations the flat space of invertible
Organization of the paper In Section 2 we present the Wasserstein metric and its restriction to Gaussian distributions, as well as the related Mures-Wasserstehi and Otto-Wasserster geometries. We present GPCA for centered Gaussian distributions in Section 3 and the general case of a.c. probability measures in Section Experiments are presented in Section and the paper ends with discussion in Section All the proofs and additional experiments are deferred to the appendices.
2 BACKGROUND
The Wasserstein distance Optimal transport is about finding the optimal way to transport mass from one distribution μ on Idd to another with respect to a ground cost, say the Euclidean squared distance. The total transport cost defines the Wasserstein distance W2 between a.c. measures μr, with moment of order 2, whose Monge (178I) formulation is given by
W (u, v/ = { | Th (x)
and where the map T is the f-a.s. unique gradient of a convex function verifying Tp # V. as proven by Brenier (199T). When the distributions ju and v are centered (non-degenerate) Gaussian distributions, they can be identified with their covariance matrices En. Σ, and the induced distance on the manifold S++ of symmetric positive definite (SPD) matrices is called the Bures- Wassersesssrttn distance BW₂ (see e.g. Modin (017)); Bhatia et al. 2019))
BW Ex)
Exp E,
Both distances can be induced by a Riemannian metric on their respective manifolds, i.c. the space of a.c. distributions and S++, we will see in the following. For more details, see Appendix B
Bures- Wasserstein geometry of centered Gaussian distributions The set of centered non- degenerate Gaussian distributions on Rdd is identified with the manifold S++ of SPD matrices.
2

---

Published as a conference paper at ICLR 2026
The Riemannian geometry of the Bures- Wasserstein metric in equation can be deseribed by con- sidering St as the quotient of the manifold L of invertible matrices by the right action of the orthogonal group Oa. In this geometry, GLa is decomposed into equivalence classes called fibers. The fiber over E S2* is defined to be the pre-image of Σ under the projection
T: A E GLa + AAT S
and can be obtained as the result of the action of O2 on representative, e.g. 11/2 the only SPD square root of 2: π. (22) {A E GLd, AAT 2 ///20aa
Tangent vectors to G are said to be horizonial if they are orthogonal to the fibers with respect to the Frobenius metric, i.c. if they belong to the space
GL
honaontal al grodes
A2
HOr (X E Raxd X A A 0}, (5)
1d
4
c
Hora
for given point A E GLg- Then the projection in equation 4 defines an isometry between the horizontal subspace Hor equipped with the Frobenius inner prod- uct (X = tr(XYT). and S equipped with a Riemannian metric that induces the Bures- Wassersissn distance (equation 3 as the geodesic distance. In partic- ular, this means that moving horizontally along straight lines in the top space GLd is equivalent to moving along geodesics in the bottom space ST (see Figure as recalled in the following proposition.
(ld)
Wassenseein goodesic
Σ1
The Bures- Wasserstein geom- etry of centered Gaussian distributions, inspired by Khesin al. (202T).
Figure 1:
Proposition (Takatsu (2011; Malagó al. (2018); Bhatia et al. (2019)). ). Any geodesic E(t) in
Sá S++ for the Bures- Wasserstein metric equation is the n-projection of horizontal line segment
in GLa that is
E(t) - (x) =
-A+
AE GLa, XE HOrA,
(6)
where is defined in certain time interval (tmins m Also, the Bures- Wasserstein distance between two covariance matrices Σ1,E, St++ is given by the minimal Euclidean distance between their fibers
BW Σ2)
inf Qu02
- 121-E1220 11/22211 -
inf QESOa
|||//2 1/2 =1/²
where || is the Frobenius norm and SO₂ is the special orthogonal group.
essential to note that the geodesic equation 6 cannot extended for all time (the only geodesic lines are those obtained by translation (Kloeckner) 2010, Proposition 3.6)). Therefore, equation is only defined on a time interval (Amin: /mmxx that depends on the eigenvalues of X A (see Appendix B.3). More details the Bures- Wasserstein geometry can be found in Appendix B.2
Otto- Wasserstein geometry of a.c. probability measures The Riemannian structure described for Gaussian distributions is a special case of Otto (2001)ss 's more general construction the bottom space becomes the space Prob([?) of a.c. distributions supported on a compact set f2 C IR while the top space is the space of diffeomorphisms Diff (i) endowed with the 7,2 metric with respect fixed reference measure (see Figure 19 in Appendix B}. The fibers of (1f are then defined to be the
pre-images under the projection
:: o E Diff + ((p) = P*p p Prob()).
(8)
In this setting, horizontal displacements Diff are along vector fields that are gradients of functions. The projection æ defines an isometry between the horizontal subspace equipped with the (p)))innee product and Prob(?)) equipped with a Riemannian metric that induces the Wasserstein distance as the geodesic distance. In particular, we have the following result.
Propasition (Otto (001I). Any geodesic t + jut)) for the Wasserstein metric given in equation 2 is the T-projection of a line segment DII going through diffeomorphism sp at horizontal speed V f sp. for some smooth function C(RR).. That is, for I defined in certain interval (Imins mmx)),
v(t) Tlpp + VV f () (id +1V f) ()P+P).
3

---

Published as a conference paper at ICLR 2026
Another geodesic ju(t) T(p + +Vf () is orthogonal to puf) at 0 for the Riemannian metric inducing the Wasserstein distance if and only if ( s,, Vf () [P(p) p((p) 0.
We emphasize that need not be convex in equation unlike in the more classical parametrization of geodesics due to McCann 1997) between two distributions Ho and μ1 = (u
(10) Note that equation parametrizes geodesics provided that id +/V f is diffeomorphism, and thus it is defined on time interval that depends on the cigenvalues of the Hessian of f. On the other hand, the convexity condition on the function u in the parametrization of equation [0 enenuress that time t is defined on [0, 1]. Both are completely equivalent (see Appendix B.3 for details).
((t) (id + Vu id))+Mo. with & E [0, 1| and
convex function.
3 GEODESIC PCA ON CENTERED GAUSSIAN DISTRIBUTIONS
In this section, we consider the exact GPCA problem for the Bures- Wasserstein metric in equation The data are centered Gaussian distributions identified with their covariance matrices E...... 2n
S++ Following Huckemann al. (2010), we define the first component the geodesic t + 21)
S that minimizes the sum of squared residuals of the 2W2--rojection of the data:
inf b++)(t)) geodesic
Σ inf BW2 Ex). t
(11)
The second principal component is defined to be the geodesic that minimizes the same cost function, with the constraint of intersecting the previous component orthogonally. The subsequent principal components have the additional constraint of going through the intersection of the first two principal geodesics. This definition does not impose that the geodesic components go through the Wasserstein barycenter (see (Agueh& & Carlier 20)1)). and in Section 5 we show an example where this is indeed not verified. This gives an observation of the phenomenon already described in Huckemann & Ziezold (2006) for spherical geometry. The proofs of this section are deferred to Appendix
Learning the geodesic components Following Propositon I we lift the GPCA problem equa- tion [I to the total space GLa of Otto's fiber bundle. This has several advantages: the Bures- Wasserstein distance in the cost function of equation is replaced by the Frobenius norm ||. the geodesic is replaced by horizontal line segment, and the projection times 4 become explicit. The price to pay is an optimization over variables (Q+) in SO,, needed to represent the covariance matrices Σ, by invertible matrices D./²Q, Qu in their respective fibers.
Proposition 3. Let T: GLd + S A + AAT and (A1. X1, (Q.)) be solution of
a
inf F(A1, (Q.) Σ ||A + ( -2 ///²Q.112.
(12)
subject to A GLd, X1 E HOTA,, ||x,|| = 1, QI...., On SOa.
Then there exist mins tmax R such that the geodesic E t [tmin: max]] T(A1 + tXx) in S
minimizes equation 71
Here the t, are projection times given by t -1, X1), and PA.X is a projection operator that clips any E IR onto closed interval [fmin mmax] depending A and X, such that A PA, (t) invertible for any t in this interval (see Appendix [B.3). Clipping the time parameter of the line segment is necessary ensure it remains within GT. and projects onto geodesic in St+.
The second component is thus defined by E2(t) = TAA + tX2), found by solving:
inf F(A2, X2, (Q:)
subject to
A = (Ai + R") R' E S0a, 1" [Imins (max)
(13)
X₂ E Hor
|X 2 (X2,XIR") = 0, Qi...., Qn E SOa-
Note that this step requires to find new rotation matrices (Q:) The first two components fix the intersection point π (Ag) through which all other geodesic components will pass, see Figure 2 For every higher order component, the fiber over (A2)) and orthogonal to the lifts of the velocity vectors of the previous components. Details on the implementation of these components given in Appendix D.2
search for velocity vector X₁ that is horizontal at some point in
4

---

Published as a conference paper at ICLR 2026
The second component is geodesic of St that orthogonally intersects the first component. Lifting again the problem to GLa. this boils down to searching for horizontal line + A2 IX2 where A2 = (Ai + 1" IR for a rotation ma- trix R2, time 1* mm /max] and horizontal vector X2 E Hor such that X2,)) 0. The equation for Ae ensures that the я p--proeectoons of the first two horizontal lines intersect, while the condition on X2 ensures that they intersect orthogonally (since Xi is horizontal at Ag easily be checked). See Figure
GLd
A
X1
A
Aa
X1R
X2
(TA,))
T: (A2))
dm d,
S
d (X2)
Figure 2: First (red) and second (blue) geodesic compo- nents of Gaussian GPCA, where d denotes the dif- ferential of the projection æ:: A AA GLd-
Quantifying the difference between TPCA and GPCA In the following, we quantify the distortion induced by linearization in the case of covariances matrices with same eigenvalues.
Proposition 4. Let E E + with eigenvalues a² i and E' = PaEP2 where Po is the rotation matrix of angle 0. Then, denoting = ((a + b)/2)² I, we have
BW3 BW2 2.8 (Σ
a b (a=b) cos² 0 + O((a (()-)))) (a+b)
(14)
where BW2,5 2,5, the linearized Cures-Wassersterr distance at E. recalled in equation 29
For given 0, equation 14 shows that the distorsion is most important for la+b close to 1, which corresponds to matrices that are close to the border of the cone, as illustrated in Section 5.1
[a-b]
On the restriction to the space of Gaussian distributions Geodesic PCA can also be defined in the more general space of a.c. probability distributions, as presented in Section 4 A natural question that arises is whether performing GPCA in the whole space of probability distributions gives the same result as restricting to the space of Gaussian distributions, which is totally geodesic. The answer is yes in dimension one, shown in Appendix D
Proposition 5. Let = Amm,, of) for i 1...., n. be n univariate Gaussian distributions. The first principal geodesic component & E [0, 1] + put) solving equation remains in the space of Gaussian distributions for all L E [0, 1].
Up to our knowledge, this remains an open question in higher dimension.
GEODESIC PCA ON A.C. PROBABILITY MEASURES: GPCAGEN
We now tackle the task of performing GPCA on a set a.c. probability measures M,,..., Un using the
geometry. We propose parameterization of the geodesic principal components based on Otto's 's formulation, leveraging neural networks. Additionally, we introduce a dedicated cost function to optimize the different geodesic components.
Parameterizing geodesics Following Proposition 2 and equation any geodesic 1. + (u(t) in the Wasserstein space (Prob(S2), W2) can be expressed (u(t) ( IV
(p) +) for t in some interval ['min mmax]. +p : [ + 1d diffeomorphism, ; R2 IR smooth function, and p fixed reference measure, taken to be the standard Gaussian distribution in this work. Using multilayer perceptrons (MLPs) to parametrize the functions & and f. denoted (P9 and fø, respectively, the curve
+ Ha (t) = (id + tV fw) (Po#P)
geodesic for E [mminn | provided that id + fV E Diff( for all t in this interval. Equivalently, this condition holds if the Hessian matrix La + tH (x) is positive definite for all Idd and & [mming (max], where H (I) denotes the Hessian of f at In practice, we enforce
5

---

Published as a conference paper at ICLR 2026
this constraint by monitoring the eigenvalues of I (z) (see Appendix B.3 and either clipping t or adjusting the interval [minn fmax] to ensure that all eigenvalues remain positive. This representation enables to sample from the distributions along the geodesic. Specifically, given the learned vector field se and function fa, one can sample from (t) by first drawing p and then applying the transformations bpy and + fv sequentially as (x) + (ypo ~ (t,
Learning the geodesic components The first principal component in GPCA minimizes the objec- tive in equation! [ The scalar variables & specify the projection time of each distribution vs onto the geodesic t ++ p(t). Leveraging the explicit form of Otto's geodesic, equation I can be reformulated
as:
inf FEC(R") (EDi((()) (1.... Ea tamxi
pff.ptt....... t.): = Σ W2((dd++VV))(<P))))
(15)
We jointly learn the parameters t, together with the neural networks sp and fo to minimize the the Sinkhorn divergence Sg that has been proven to be differentiable and computationally efficient approximation of the squared Wasserstein distance W see/ Frogner al. ): Genevay et al. (2018); Chizat et al. (2020), and represent the distributions p and P4 using batches of m samples parameters based on a single distribution v4 sampled each iteration, detailed Algorithm To compute min and mmax on line of Algorithm| [ we approximate the extremal eigenvalues of Hf by evaluating the largest and smallest eigenvalues over the finite set {H (Fx)) ) and substitute these estimates into the theoretical bounds from Appendix B.3
objective in equation| 13 In practice,
and V/5 ~ The optimization proceeds by updating the
Algorithm 1 Geodesic PCA algorithm for a.c. measures: GPCAGEN
1: Initialize po. fú and the & for 1 iis
2: while not converged do
3: for = 1 to n do
4: Draw i.i.d samples yj vs and draw i.i.d samples Fy ~ 1 1, S m
5: Estimate mins mmax with { [= and set E = min(max((ss, tmin), tmax)
2 (6) (id VVV s+) (FF) for k
6:
7: Ca,y. + St m Ek= 8 6_m,, m
ΣØ= 5 v3
8:
Update spe. S and the 4 with V VL.
9: end for
10:
end while
The second principal component minimizes the objective in equation [ subject to the constraint that it intersects the first component orthogonally. Similar to the first component, we use two MLPs, fea
and Pez, parameterize the geodesic ++ Me,.ve whe (t). along with n scalar variables tf. to optimize the objective in equation 15 We also introduce two additional scalar variables, e and 22 meerr which define the intersection times of the two geodesics, along with the regularization terms:
(g. 1g
I(Ei, Ez, Limer (macr) = ||E: (1mer)) (2 and Ogg,h) 110122
1 11:
where I enforces the two geodesics in Diff (1) (id +/V oyo and (2(1) (id +zV
SPe.. to intersect the respective times me and e while O(g, h) ensures orthogonality between the corresponding horizontal vector fields g = fo(())) and h V (Po,) in L²(p). The total objective used to optimize the second principal component incorporates these regularization terms and is given by:
[(Ser: Pog, .....,(e) + A Eeg M ( O (s,,(
with (t) = (tv o ie and (0,,v(()) (t) = (id +tV f fwz) Po, and where x and ho are the regularization parameters controlling the of between the intersection and orthogonality regularization terms, respectively. Note that in virtue of Proposition 2 the Lp) inner product in the regularization term O truly enforces orthogonality of the geodesic components with respect

---

Published as a conference paper at ICLR 2026
to the Riemannian metric associated the Wasserstein distance. Note also that I enforces the geodesics to intersect in (12 which means that, at the intersection time, the geodesics H and 10 in Prob([t) intersect and share the same representative. An alternative implementation would be to enforce the intersection the geodesics 14 and P02.v22 in Prob([2) and to impose the orthogonality V & (PP)) o Rt and V (Pe.) in L (2)),, where R = E0g.v/ (fmer) LB,,v, (f)mer)) This approach is the used in the Gaussian case. However, computing R is computationally expensive, and we therefore preferred to impose Eo,,,, (fmec) = 60 (fatte) which directly yields R id.
The training algorithm used to optimize the second principal component follows the same structure as Algorithm [ except for the seventh line, where the regularization terms, estimated using the minibatch Ik ~ are added to the loss function. Higher-order components can be computed similarly.
5
EXPERIMENTS
5.1 EXPERIMENTS ON CENTERED GAUSSIAN DISTRIBUTIONS
In this section, consider toy examples in -5t and compare GPCA to its widely used linearized approximation, TPCA (see Appendix C). We use coordinate systems for matrices in Sh the first comes from the spectral decomposition, and the second maps any SPD matrix to a point in the interior of the cone C {1 U, 2) E RRa,
2² x² +2):
= Pa
a²
0 6²
P
+
x - (a,b,0) ( X RR (I, (x) EC (16)
where Pe is the rotation matrix of angle 9. Generically, GPCA and TPCA yield very similar results: for sets of n = 50 covariance matrices randomly generated using a uniform distribution on the parameters (a, b, 0). GPCA reduces the objective in equation o1 of less than 1% w.r.t TPCA, average for 100 trials. This suggests that TPCA is generally a very good approximation of GPCA. Two extreme cases are described below: GPCA and TPCA are equivalent and (ii) GPCA and TPCA drastically differ.
Matrices with same orientation we consider set of covariance matrices that live in the subspace 0 constant in notations of equation I6 then both GPCA and TPCA yield exactly the same results,
namely that of linear PCA in the (a, b) boooordiisssss. This is because any such subspace has zero curva- ture for the Wasserstein metric, and geodesics are straight lines in the (a, b) coorrdinasss (Appendix D.I] FigureBJshows the geodesic compo- nents obtained for a of matrices in the subspace = 0 that form regular rectangular grid in the (a, b) coordinates, i.c. Σy = diag(a; (a? where the a, 's and by 's are equally spaced. They are indeed straight lines that capture the variations in and b respectively.
Figure 3: GPCA on a set of diagonal covariance matrices Enj with varying eigenvalues < af < 3. by < The matrices form a planar grid inside the cone C of SPD matrices in equation 16 (left), and correspond to ellipses of varying width and height (right). The first component (red) captures the variation in a, while the second component (blue) captures the variation in b.
Matrices with same eigenvalues Now we consider covariance matrices that all have the same cigenvalues but different orientations. Specifically, we choose Σ; = de.diaa((e) 6² PT, for positive reals > b, 0; = in/n h for i O,..., n and an even number n. the (x, y, z) coordinates (equation 16), the covariance matrices are displayed on circle of equation and x² + = cst (constant determinant), as shown in Figure (in practice, we choose a slightly open circle to break the symmetry). Then the Bures- Wasss-sesnnnttinn barycenter of the covariance matrices E....., Ex, is given by Σ (a + b)² /4I (see Proposition 15 in Appendix D.I.. When performing TPCA on Σ,, E1,..., En, the barycenter Σ, the radial distances between Σ. and E; preserved, but
cst (constant trace)
7

---

go through Published as a conference paper at ICLR 2026
Cost imprrveemmee (5G
10-
0.5
-0.5
0.0
0.5
-0.5
[e-b]/ [a+b]
Figure 4: Comparison between tangent and geodesic PCA on set of n = 20 covariance matrices with same eigenvalues a², 62 and different orientations 0. (left) They equally spaced on an (open) circle in a horizontal plane inside the cone of SPD matrices. The first component of TPTC (dashed red line) goes through the Fréchet mean Σ (magenta dot), a multiple of the identity, while the component of GPCA (solid red line) does not. Here la b + bl ~ 0.8. (middle) Representation of the left figure in the (z. y) coordinates. (right) Evolution of the first component cost improvement (in the sense of minimization) of GPCA with respect to TPCA, function of the ratio la - bl/ aa [a++a bl.
the pairwise distances between the Σ, Proposition Proposition evaluuaes the level of this distorsion. Note that since (a b)2 /(a + 6)2 (x² + /²)/22, the distorsion is most important when covariance matrices are close to the border of the cone, see Figure
Indeed, in that case, the results of GPCA can be very different from those of TPCA and the first component may not
the Wasserstein barycenter Σ, see Figure (middle) and Figure g in Appendix A In that case GPCA worse-behaved TPCA, some of the Gaussian distributions will project onto
may be scen the first geodesic component boundaries, yielding a poor separation. Figure 4 (right) shows the percentage of improvement of the cost in equation [T (in terms of minimization) of GPCA with respect to TPCA, in the setting previously described for different values of the ratio |a - bb/ ///a+ 6]. on average for 10 runs per value of the ratio. The blue strip indicates standard deviation.
Weather dataset In this paragraph, we use the Weather CORGIS Dataset to illustrate GPCA based on empirical covariance matrices. The dataset provides weekly measures of precipitation and wind speed recorded from March to January 2016 across the 50 U.S. states and the territory of Puerto Rico. From these measures, we construct two histograms for each state: one for precipitation and one for wind speed. We then compute the 50 empirical covariances from these histograms. We show in Figure 14 the projection of each state onto the two first GPCA components computed from the empirical covariance matrices. We can clearly identify clusters of different weather behavior among the states.
5.2 EXPERIMENTS ON ABSOLUTELY CONTINUOUS DISTRIBUTIONS
We conduct preliminary experiment on synthetic dataset with known geodesics verify that our algorithm, GPCAGEN (Section<) 4). accurately recovers the first two principal components. We then apply GPCAGEN 3D point clouds from the ModelNet40 dataset (Wu al. (2015)) and to color distributions of images from the Landscape Pictures dataset (Rougetet (020)). An additional experiment in Appendix(A. A.33 demonstrates how GPCA can be used for outlier detection. For these experiments, fo and j0g are ML.Ps with four hidden layers of size 128 and an output layer of size 1
8
8
8
8
1
8
8
3
22
2
Figure Densities of probability distributions uniformly sampled along the first and second principal geodesics components. GPCAGEN successfully recovers the two orthogonally intersecting geodesics constructed from MNIST data. The first component (left) captures variation color space, while the second component (right) recovers the interpolation from the digit "1" to the digit "2"
8

---

Published as a conference paper at ICLR 2026
1
3
-2 -1
1
-3
-2
-1
0
2
3
t=
t= 0.88
1=11
3.0
t=
-2.4
0.2
Figure 6: Empirical distributions uniformly sampled along the geodesics corresponding to the first (left) and second (right) principal components, computed by GPCAGEN in the 3D point cloud of chairs experiment (top row). the 3D point cloud of lamps experiment (middle row) and the Landscape images experiment (bottom row).
and d respectively. We found that setting the regularization coefficients A and λo to 1.0 ensures the algorithm works as expected in all experiments. A discussion of the regularization coefficients, along with details on the architecture and hyperparameters, provided in Appendix
MNIST geodesics. We represent each image from the MNIST dataset (LeCun et al. (2010)
probability measure over R4, The grayscale pixel intensities define a normalized density over spatial coordinates (x, yy IR2, and further assign each pixel two additional values corresponding red and blue color channels. We construct two orthogonal geodesics; the first one interpolates between digit "I" and digit "2", both assigned fixed purple by setting the color channels to 0.5. The second one is defined from the midpoint of the first, by linearly interpolating the color from red to blue. As shown Figures 5 and GPCAGEN successfully recovers the two geodesics intersecting orthogonally. A second experiment on the MNIST dataset is displayed in Appendix
3D point cloud. We use the ModelNet40 3D point cloud dataset Wu et al. 2015 and apply GPCA subset of 100 randomly selected lamp point clouds, Figure (middle row) and Figure (left) demonstrate that the first principal component captures the distinction between hanging lamps (chandeliers) and standing lamps (floor lamps), while the second component reflects variations in the thickness of the lamp structure. We conduct similar experiment ModelNet40 representing different chairs. As shown in Figure 6 (top row) and Figure [0, the first principal component distinguishes between chairs and armchairs, while the second component captures the height of the seat.
100 point clouds from
took 39 images from the Landscape Pictures dataset (Rougetet) (2020) and use GPCAGEN the corresponding point clouds, where each point cloud represents color distribution in the image. ( row) and Figure ( (right) show that the first component captures variations in overall brightness, ranging from dark to bright images, while the second component separates mostly blue images from mostly green ones.
Landscape images.
Baselines An obvious baseline for GPCAGEN is TPCA. Unlike GPCAGEN, which learns continu- geodesics from empirical distributions of absolutely continuous measures, TPCA acts on discrete measures. A direct numerical comparison between the two methods is therefore not meaningful. However, include in Appendix A.2 the two principal components returned by TPCA the 3D point cloud experiments. We observe in Figure I6 that the discrete nature of TPCA produces artifacts, including holes in certain regions and excessive concentration in others.

---

Published as a conference paper at ICLR 2026
Another natural baseline consists in embedding point clouds into latent space of dimension d then performing standard PCA on the resulting latent vectors. This approach, in addition to being computationally expensive, does not produce meaningful modes of variation, shown in Section A.2 of the appendices.
3
1
0
-1
-4
-2
-1
0 d))
1
12)
Figure 7: Each lamp point cloud (left) and each image (right) is embedded in the plane according to its projection times onto the first and second principal components computed by GPCAGEN.
DISCUSSION
GPCA a statistical approach learning the main modes of variations of a set of probability distributions. The first components capture meaningful structure for data lying on curved space, which then enables downstream tasks such as classification, clustering, and outlier detection. In this work, have proposed two methods for computing exact GPCA distributions and the other for the more general case of a.c. probability distributions. In the Gaussian case, experiments suggest that GPCA and TPCA generically yield very similar results, except for distributions with covariance matrices that are close to the boundary of the SPD cone, for which GPCA can yield undesirable effects as suggested by the pathological example of Figure In the general case of a.c. probability measures, a key advantage of our approach is that operates directly on continuous distributions, avoiding the need for empirical approximations sample sizes and can introduce discretization artifacts in the recovered components. Additionally, our method enables sampling from any point along the geodesic components- something not possible with discrete approximations commonly used in TPCA. Otto's parametrization also allowed us to avoid relying on input convex neural networks (ICNNs) by not requiring convex functions, with the trade-off being the need to estimate the eigenvalues of the Hessian of f. This perspective opens new directions for parametrizing convex functions without imposing hard architectural constraints.
tailored for Gaussian
the v4, which would require equal
ACKNOWLEDGEMENT
This work benefited from the support of the Agence nationale de la recherche, through the PEPR This work performed using HPC resources from GENCI-IDRIS (Grant 2023-103245). This work was partially supported by Hi! Paris through the PhD funding of Nina Vesseron.
PDE-AI project
10

---

Published as a conference paper at ICLR 2026
REPRODUCIBILITY STATEMENT
All implementation details our proposed method, including model architectures, training proce- dures, and hyperparameter settings, are provided Section the main paper and in Appendix D Dn2 Original theoretical results are presented with complete proofs in Appendix D The datasets used our experiments are publicly available. We will release the source code reproduce all experiments associated with this paper at a later stage.
REFERENCES
Martial Agueh and Guillaume Carlier. Barycenters in the Wasserstein space. SIAM Journal on Mathematical Analysis, 43(2):904-924, 2011.
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media, 2008.
Luigi Ambrosio, Alberto Bressan, Dirk Helbing, Axel Klar, Enrique Zuazua, Luigi Ambrosio, and Nicola Gigli. user's guide to optimal transport. Modelling and Optimisation of Flows on
Networks: Cetraro, Italy 2009, Editors: Benedetto Piccoli, Michel Rascle, pp. 1-155, 2013.
Rajendra Bhatia, Tanvi Jain, and Yongdo Lim. On the Cures-Wasserstee distance between positive definite matrices. Expositiones Mathematicae, 37(2):165-119) 191 2019.
Jérémie Bigot, Raúl Gouet, Thierry Klein, and Alfredo López. Geodesic PCA in the Wasserstein space by convex PCA. Annales de l'Institut Henri Poincaré Probabilités et Statistiques, 53(1): 1-26, 2017.
Emmanuel Boissard, Thibaut Le Gouic, and Jean- Loubes. Distribution's template estimate with Wasserstein metrics. Bernoulli, 21(2):740-759, 2015.
Nicolas Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, 2023.
Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. Commu- nications on pure and applied mathematics, 44(4):375 -47, 1991.
Anna Calissano, Aasa Feragen, and Simone Vantini. Populations of unlabelled networks: Graph space geometry and generalized geodesic principal components. Biometrika, pp. 147-170, 2024.
Elsa Cazelles, Vivien Seguy, Jérémie Bigot, Marco Cuturi, and Nicolas Papadakis. Geodesic PCA versus log-PCA of histograms in the Wasserstein space. SIAM Journal Scientific Computing, 40(2):B429-B433 B456, 2018.
Lenaic Chizat, Pierre Roussillon, Flavien Léger, François-Xaver Xavier Vialard, and Gabriel Peyré. Faster wasserstein distance estimation with the sinkhorn divergence. Advances in neural information processing systems, 33:2257-2269, 2020.
Jed A Duersch, Meiyue Shao, Chao Yang. and Ming Gu. A robust and efficient implementation of LOBPCG. SIAM Journal on Scientific Computing, (0(5):C655-C677 2018
David G Ebin and Jerrold Marsden. Groups of diffeomorphisms and the motion an incompressible fluid. Annals of Mathematics, 92(1):102-16) 163, 1970.
Thomas Fletcher, Conglin Lu, and Sarang Joshi. Statistics shape via principal geodesic analysis on Lie groups. In 2003 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003. Proceedings., volume 1, pp. I-I. IEEE, 2003.
Charlie Frogner, Chiyuan Zhang, Hossein Mohahi, Mauricio Araya, and Tomaso A Poggio. Learning with wasserstein loss. Advances in neural information processing systems, 28, 2015.
Aude Genevay, Gabriel Peyré, and Mareo Cuturi. Learning generative models with sinkhorn di- vergences. In International Conference on Artificial Intelligence and Statistics, pp. 1608-1617. PMLR, 2018.
11

---

Published as a conference paper at ICLR 2026
Philipp Harms, Elodie Maignant, and Stefan Schlager. Approximation of Riemannian distances and applications to distance- based learning on manifolds. arXiv preprint :1904.11860, 2019.
Wen Huang and Ke Wei. Riemannian proximal gradient methods. Mathematical Programming, 194 (1):371-413) 2022.
Stephan Huckemann and Herbert Ziezold. Principal component analysis for Riemannian manifolds, with an application to triangular shape spaces. Advances in Applied Probability, 38(2):299-319, 2006.
Stephan Huckemann, Thomas Hotz, and Axel Munk. Intrinsic shape analysis: Geodesic PCA for Riemannian manifolds modulo isometric Lie group actions. Statistica Sinica, pp. 1-58, 2010.
Boris Khesin, Gerard Misiotek, and Klas Modin. Geometric hydrodynamics and infinite-dimensional Newton's 's equations. Bulletin of the American Mathematical Sociery, 58(3):377-44 442, 2021.
Benoit Kloeckner. A geometric study of Wasserstein spaces: Euclidean spaces, Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 9(2):297-323, 2010.
Yann LeCun, Corinna Cortes, and CJ Burges. MNIST handwritten digit database. ATT Labs [Online].
Available: http://yann.lecu htw.://www.cc...o/////////mmmeae cow./exm/////isstt..ooooommaama. 2, 2010.
Cyril Letrouit and Quentin Mérigot. Gluing methods for quantitative stability of optimal transport maps. arXiv preprini arXiv:241 1.04908, 2024.
Chi-Kwong Li and Fuzhen Zhang. Eigenvalue continuity and Gersgorin's 's theorem. Electronic Journal of Linear Algebra, 2019.
Yingjie Ma, Xi Gao, Chao Liu, and Jie Li. Improved SQP and SLSQP algorithms for feasible path-based process optimisation. Computers & Chemical Engineering, pp. 108751, 2024.
Luigi Malagö, Luigi Montrucchio, and Giovanni Pistone. Wasserstein Riemannian geometry of
Gaussian densities. Information Geometry, 1:137-179, 2018.
Robert McCann. A convexity principle for interacting gases. Advances in mathematics, 128(1): 153-179, 1997.
Nina Miolane, Nicolas Guigui, Alice Le Brigant, Johan Mathe, Benjamin Hou, Yann Thanwerdas, Stefan Heyder, Olivier Peltre, Niklas Koep, Hadi Zaatiti, Hatem Hajri, Yann Cabanes, Thomas Gerald, Paul Chauchat, Christian Shewmake, Daniel Brooks, Bernhard Kainz, Claire Donnat, Susan Holmes, and Xavier Pennec. Geomstats: A Python package for Riemannian geometry in machine learning. Journal of Machine Learning Research, 21(223):1 2020. URL Stthttp
Klas Modin. Geometry of matrix decompositions seen through optimal transport and information geometry. Journal of Geometric Mechanics, 9(3):335 390, 2017.
Gaspard Monge. Mémoire sur la théorie des déblais et des remblais. Mem. Math. Phys. Acad. Royale Sci., pp. 666-704, 1781.
Felix Otto. The geometry of dissipative evolution equations: the porous medium equation. Communi- cations in Partial Differential Equations, 2001.
Michael J.D. Powell. On search directions for minimization algorithms. Matherattiaal programming, 4:193-201, 1973.
Charles R Qi. Hao Su, Kaichun Mo, and Leonidas J Guibas. Pointnet: Deep learning point sets for 3d classification and segmentation. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 652-660, 2017.
James O. Ramsay and Bernard W. Silverman. functionaal data analysis: methods and case
studies. Springer, 2002.
Araad Rougetet. Landscape pictures dataset.
httpp://ww.kagggee.co cttp://le..oo///aattttttt
2020.
12

---

Published as a conference paper at ICLR 2026
Klaus D Schmidt. On inequalities for moments and the covariance of monotone functions. Insurance: Mathematics and Economics, 55:91-95, 2014.
Vivien Seguy and Marco Cuturi. Principal geodesic analysis for probability measures under the
optimal transport metric. Advances in Neural Information Processing Systems, 28, 2015.
Stefan Sommer, François Lauze, and Mads Nielsen. The differential the exponential map, Jacobi fields and exact principal geodesic analysis. CoRR, abs/1008. / /902, 2010.
Stefan Sommer, François Lauze, and Mads Nielsen. Optimization geodesics for exact principal geodesic analysis. Advances Computational Mathematics, 40:283-313, 2014.
Asuka Takatsu. Wasserstein geometry of Gaussian measures. Osaka J. Math., 48: 1005-1026, 2011
Yann Thanwerdas. Riemannian and stratified geometries on covariance and correlation matrices. PhD thesis, Université Côte d' Azur, 2022.
Pauli Virtanen, Ralf Gommers, Travis Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, ilhan Polat, Yu Feng. Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272 2020. doi: 10903//41192-09--00066
Wei Wang, Dejan Slepčev, Saurav Basu, John Ozolek, and Gustavo K Rohde. A linear optimal transportation framework for quantifying and visualizing variations in sets of images. International journal of computer vision, 101:254-269, 2013.
Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao. 3d shapenets: A deep representation for volumetric shapes, 2015. URL https:://s/e
axq/abs/14055
Constantin Zalinescu. Convex analysis in general vector spaces. World scientific, 2002.
Juan Zhou, Kangkang Deng, Hongxia Wang, and Zheng Peng. Inexact Riemannian gradient descent method for nonconvex optimization, arXiv preprint arXiv:2409. 2024.
13

---

Published as a conference paper at ICLR 2026
A
ADDITIONAL EXPERIMENTS AND FIGURES
A.1
GEODESIC PCA
Here present additional figures to further explain the experiments described in the paper. Figure concerns the experiment on Gaussian distributions with diagonal covariances described Section 5.1 corresponding to Figure 4 It shows all three principal components found by tangent PCA (left) and geodesic PCA, in two equally optimal solutions (middle, right).
0.5
0.5
0.0
0.0 0.0
0.0
0.0
0.0
0.5
-0. 0.5
-0.5
Figure 8: Principal geodesic components of a set of Gaussian distributions whose covariance matrices have same eigenvalues and different orientations, as described in Section 5. Tangent PCA yields unique solution (left) where geodesic components cross at the barycenter, while geodesic PCA yields two equally optimal solutions (middle, right) where the geodesic components cross at another point. The first geodesic component is shown in red, the second in blue, the third in green.
Figure displays the plane the first geodesic components of the MNIST experiment of Section 5.2. while Figure I0 shows the planar representation of the 3D point cloud of chairs experiment given by the projection onto the first two geodesic components found by GPCAGEN algorithm and depicted in Figure 6 (top row).
4
2
2
0
-2
-2
0
-2
1
2
Figure 9: Each point cloud, corresponding to distribution along one of the artificially con- structed geodesics, is embedded in the plane ac- cording to its projection times onto the first and second geodesics returned by the GPCAGEN algorithm. We observe that GPCAGEN success- fully recovers the two orthogonally intersecting geodesics designed from MNIST-based interpo- lations of digit shape and color.
Figure 10: Each chair point cloud is embedded in the plane according its projection times onto the first and second geodesics returned by the GPCAGEN algorithm.
14

---

Published as a conference paper at ICLR 2026
Finally, present additional experiment the MNIST dataset. We use the same color construc- tion as in the experiment presented in Section we then apply GPCAGEN to a dataset of 20 red digits "1", 20 blue digits "1", 20 red digits "2" and 20 blue digits "2" (see Figure 12] As shown in Figures 11 and 12 GPCAGEN again identifies two orthogonal geodesics: the first primarily captures variation in color, while the second captures variation in shape-from from digit "2" to digit "1".
0
2
2
2
Figure 11: Densities of probability distributions uniformly sampled along the geodesics corresponding to the first and second principal components. The first component (left) returned by GPCAGEN captures variation in color space, while the second component (right) recovers the interpolation between digit "2" and digit "1".
2
2
-3
-2
2
3
4
Figure 12: Each MNIST digit is embedded in the plane (the arrows indicate the exact position of each digit) according to its projection times onto the first and second geodesics returned by the GPCAGEN algorithm. We observe that the first principal component recovered by GPCAGEN captures variation in color, while the second component reflects the transformation from digit "2" to digit "1".
A.2 COMPARISON OF GPCA TO RELATED METHODS
Other notions of PCA on Gaussian distributions There exist wide variety of metrics on the space of symmetric positive definite matrices, such as e.g. the log-Euclidean, Buclidean-Cholesky or affine-invarian metrics (see Thanwerdas (2022) for comprehensive overview), Each of these metrics could be used to perform PCA on centered Gaussian distributions. However, there is no obvious quantitative way to compare the results. Each method optimizes its own criterion, and any metric that one could think to compare the methods would rely on a choice of underlying metric on the space of SPD matrices. Comparison of PCA methods with two different metrics thus boils down to comparing the metrics themselves. We illustrate in Figure B the behavior of covariances matrices along geodesics for different metrics.
15

---

Published as a conference paper at ICLR 2026
0
0
0
0
Figure 13: Geodesics on the space of symmetric positive definite matrices from left to right, for (top)
the Bures-Was Wasserstein metric, (middle) the log-Euclidean metric and (bottom) the Euclidean metric on the Cholesky coefficients.
New Hampssize
Puerto Rico
HAwaii
0.04
Massachusetts
Alaska
0.02
Piorida
Califonnia
Musiissinnl
Vontana
Calorapo
-0.01
Vermar
-0.04
0.02
0.08
Figure 14: Each dot represents a state projected onto the first two GPCA components computed from the empirical covariance matrices, which are also shown in the figure.
TPCA on 3D Point Cloud Data Here we present the results returned by TPCA on the 3D point and compare them to from those obtained by GPCAGEN.
cloud experiments, Figures 16 and 15
For the lamps dataset, the first component is similar and captures the distinction between hanging and standing lamps. The second component focuses the object thickness, like the second GPCAGEN component, but also on whether mass is concentrated at the extremitics or the middle of the lamp structure.
For the chairs dataset, both geodesics obtained by TPCA resemble those returned by GPCA. However, the second TPCA component also appears to account for whether the mass is concentrated or not.
Finally, due to the discrete nature of the TPCA algorithm, we observe discretization artifacts in the TPCA components: holes in some parts of the space, mass concentration in others.
16

---

Published as a conference paper at ICLR 2026
0.3
0.2
0.2
0.1
h
0.0
-0.1
-0.1
-0.2
-0.2
-0.1
0.0
-0.4
-0.2
0.2
4/
Figure 15: For the chair and the lamp experiment, each point cloud is embedded in the plane according to its projection times onto the first and second principal components computed by TPCA.
-0.15
0.100
0.05
0.00
0.05
0.10
0.15
0.20
0.25
-0.1
0.0
0.1
0.2
0.3
-0.4
-0.3
-0.
-0.1
0.0
0.1
0.2
0.3
t
-0.2
-0.1
0.0
0.1
0.2
0.3
Figure 16: Empirical distributions uniformly sampled along the geodesics corresponding the first (first line) and second (second line) principal components, as computed by TPCA in the 3D point cloud of chairs experiment (top rows) and the 3D point cloud of lamps experiment (bottom rows),
PCA computed in the latent space of PointNet. For the 3D plint-cloud datasets, we evaluated the natural baseline that consists in embedding point clouds into latent space of dimension d
17

---

Published as a conference paper at ICLR 2026
and then performing standard PCA on the resulting latent vectors. We used a pretrained PointNet autoencoder (Qi et al. 2017) from the public repository https:: haas://githbb. wtt.//oii//////////i/ pe autoencoder trained ModelNet40, to encode each point cloud (chairs and lamps) into a d-dimensional latent representation, on which PCA was applied. Figure shows the resulting 21) projections. We observe some clustering of similar objects; for example, large lamps tend to group together in the lamp dataset, and chairs versus armchairs form distinguishable clusters. The second principal component for chairs appears correlate with the height of the seat. Beyond these observations, however, PCA provides limited separability (especially for lamps), and the recovered components are difficult interpret.
More generally. this approach presents several important limitations:
Training point-clous autoencoder requires large collection of distributions. our case (100 distributions), we need to rely pretrained autoencoder trained on related dataset.
PCA on autoencoder embeddings relies heavily on the geometry learned by the encoder. The learned geometry is not guaranteed align with the Wasserstein structure and the recovered principal components may not reflect meaningful modes variation (as observed in the experiments ahove). Moreover, for given autoencoder that wish train, different random seeds at initialization can lead to different learned geometries and thus different PCA components, which is not suitable.
1.5
1.5
1.0
1.0
0.5
0.5
0.0
-0.5
-1.0
-1.0
-1.5
-1.5
-2.0
-2.0
-3
2
Figure 17: For the lamp (left) and the chair (right) experiment, each point cloud is embedded in the plane according to its projection times onto the first and second principal components computed by the POINTNET PCA method.
A.3 APPLICATION OF GPCA TO OUTLIER DETECTION
In this section, demonstrate how GPCA can be used for outlier detection. The underlying intuition is that GPCA components capture the structure of the dataset on which they are trained, and samples from different dataset are expected lie far from the learned components Wasserstein distance. In this experiment, we use the ModelNet40 3D point cloud dataset (Wu et al., 2015) and apply GPCA to subset of 100 randomly selected chair point clouds compute the first two components. For a new point cloud X, we define its score as the sum of the Wasserstein distances between X and projections onto the first two learned GPCA components. To compute the Wasserstein distance between X and component, component, we perform a grid search over 20 equally spaced values of between fmin and max computing the Wasserstein distance between X and 2048 samples drawn from the component each f. and select the that minimizes this distance. We repeat the same procedure for the second component and sum the two minimal distances to obtain the final score.
ot emd from the POT library. Specifically, for each
We evaluate this approach on 120 point clouds: 60 new chairs (not used for training) and 60 point clouds of cars. The left histogram in Figure T8 shows the resulting scores. We observe that the scores
18

---

Published as a conference paper at ICLR 2026
of the chair point clouds (in blue) lower than those of the point clouds (in green), indicating that it is possible to detect whether point cloud is not a chair using this score. We also repeat the experiment with 60 point clouds of planes, shown in the right histogram of Figure T8 and observe that the separation between chair and plane scores is even more pronounced.
Chairs 16
Chairs
Cars
Frequency
10
0.02
0.04
0.06
0.08
0.10
0.12
0.025 0.050
0.075
0.100
0.125
0.150
0.175
0.200
Score
Score
Figure 18: GPCA scores obtained on 60 new point clouds of chairs (never seen during training) and 60 point clouds cars (left) planes (right). The separation the histograms indicates that GPCA can be used for outlier detection.
B THE WASS-RRSNNNN GEOMETRY
In this section, briefly describe the fiber bundle structure the Wasserstein space due Otto (2001), that is behind the Riemannian interpretation of the Wasserstein distance. We then present its restriction the space of centered non-degenerate Gaussian distributions, which coincides with the Bures- Wasssssesssrrtnin Riemannian geometry on SPD matrices. Finally, we relate Otto's parametrization of geodesics McCann's interpolation.
We present these well-known results without proofs and refer the interested reader to Otto 2001);
Khesin al. (2021) and (Ambrosio al. 2013 Section 6.1) 1) for more details the general setting and to Takatsu (201) Malago et al. (2018) Bhatia et al. (2099) for details and proofs in the Gaussian setting.
B. THE ОтTо- ATTO-WASSERSTEIN GEOMETRY OF A.C. DISTRIBUTIONS
Consider the space Prob(S2) of absolutely continuous probability measures with smooth densities with respect to the Lebesgue measure, and support included in compact set C IRd well as the space (If of diffeomorphisms on f. These spaces can be equipped with infiniteddimensional manifold structure, see e.g. Ebin & Marsden 1970), that we will not describe here. The tangent space of Diff (S2) at o Dif is given by
(2f ouo p, v: Ω Rd vector field).
We fix reference measure p E Prob(St) and equip Dif((S) with the 1--metric with respect to p. defined for any tangent vectors u o
E Dif as
p) f uudd
(u pp o p))
(uo (
where μ p. Then the space of diffeomorphisms can be decomposed into fibers, defined to be equivalence classes under the projection
n: Diff + Prob(!?), + P#P.
Specifically, the fiber over μ E Prob(S)) is given by π - (p) {sp E Dif((N), ppp= n1., see Figure 19 (right). The tangent space to the fiber π (1)) at y E Diff (f) and its orthogonal with respect to the L (p) -e)-metric are refered to the vertical and horizontal spaces respectively:
Verg:
ker dr
Hor = (V
19

---

Published as a conference paper at ICLR 2026
where d (I) - (r(e) denotes the differential of at yp Moving along vertical vectors in (St) means staying in the same fiber, i.e. projecting always to the same measure μ in the bottom space. On the contrary, moving along horizontal vectors means moving orthogonally to the fibers, ie., in the direction that gets fastest away from the fiber. The following proposition gives the form vertical and horizontal vectors.
Proposition 6. Let yp (if. Then
Verp {w 9, V (wy)) 0},
Hor = [Vf o p, f C (9t)
(*))
The following results state that line segments and -distances in Diff(I) can be used compute Wasserstein geodesics and distances in the space probability measures Prob([t), provided we restrict horizontal displacements.
Proposition 7. The projection :: Dif - Prob([2) is a Riemannian submersion, i.e. is an isometry for any ip E Diff(It..
d Hor,
This implies the following.
Proposition 8 (Proposition 2 in main). Any geodesic 1 + p(l) for the Wasserstein metric in equation 2 is the x-projection of a line segment in (if going through diffeomorphism ip at horizontal speed f up for some smooth function C(R.). That is, for defined in certain interval (min )
(tt) x(p+)Vfop) (id +VV (P+
(17)
Another geodesic ji(t) = (lp + fof 4) is orthogonal to p(t) at = 0 0.for the Riemannian metric inducing the Wasserstein distance if and only if fVf o, V (p))ppp))))))))))))))))) = 0.
We comment the link between this parametrization and McCann's interpolation in Section B.3
B.2 THE OTTO-WASSEREN GEOMETRY OF GAUSSIAN DISTRIBUTIONS
The Wasserstssnr distance in equation on the space S of symmetric positive definite (SPD) matrices is the geodesic distance induced by a Riemannian metric gBW, which can be written in different ways. Here we the expression from (Thanwerdas 2022
Table 4.7), defined for
= PDP E Ett+ and = PU'P" E S4. by
Σ 0110 IS4
BW g
(U,U)=
U dd
(18)
d,
where the d 's are the diagonal elements of D. The associated Riemannian geometry can be described by Otto's fiber bundle restricted to the space of centered Gaussian distributions, in the following way.
In this setting, diffeomorphisms are restricted to invertible linear maps yp: u + Au for some invertible matrix A. i.c. the space of diffeomorphisms is replaced by the Lic group of invertible matrices GLa- Tangent vectors are then given by linear maps u + Xu for any matrix standard normal distribution p (O, Id) as reference measure, the I--mmttic with respect p between u + Xu and u + Yu is then written, for any X, Y E Rdxd;
Rdxd, Fixing the
[ p(u) (e))(dp(()) = / (r(p(())(u)))ddp )dp(u) = -r
Xuu T
= tr
yielding the standard Frobenius inner product on (the tangent space of) GT... We obtain fibration of the top space GLa over the bottom space S by considering the following projection
#: GLa S
A+ AA
(19)
see Figure 19 (left). The fiber over Σ E S+++ is
( = (A E GLa. AAT = = =
(20)
where O₂ denotes the space of orthogonal matrices and 1// denotes the only SPD square root of the
SPD matrix Σ. The differential of the projection (AA) = AAT given by
d (X) XAT T AXT
(21)
20

---

Published as a conference paper at ICLR 2026
Therefore, vertical vectors, which those tangent to the fibers, equivalently, those belonging to the kernel of (X),, are given by
Vera: = {X E Bddd XAT + AX 0}
= = E Raxd XA is antisymmetric)
{XX K(AT)= K S S. (AT)-).
where S denotes the space of antisymmetric matrices of size d. Once again, moving along vertical vectors in GL means staying in the same fiber, i.e projecting always the same SPD matrix in the bottom space S7 Horizontal vectors are those that are orthogonal to all vertical vectors (for the Frobenius metric), i.e. matrices X such that for any antisymmetric matrix K:
(X, (X, tr(XA-1K))
which is equivalent to XA- symmetric (this can be seen by taking for K the basis elements of S in the above equation), yickding
Hor {xx E Raxd, -A
= =
= =(x E Rddd A x
(X KA, K E Sa) S&A
where Sa denotes the space of symmetric matrices.
Proposition 9. The projection GLd + ST+ A + AAT is Riemannian submersion, i.e. dr is an isometry from Aor a equipped with the Frobenius inner product Tx equipped with the inner product (y for any A E Gla-
Just like in the general case, this yields a way to lift the computation of geodesics and distances. Proposition 10 (Propositon in main). Any geodesic t ++ E{t) in St for the Bures- Wasserstein metric in equation] 3 is the T n-projection of a horizontal line segment in GLa, that is
E(t) x(A+tX) = (A+ + +tx) A E GLa. X E HOr
(22)
where is defined in certain time interval (Imin mmax)- Also, the Bures- Wrees-WWsserstein distance between two covariance matrices Σ,, Eg E S++ is given by the minimal distance between their fibers
B E2)
inf 1 -2/2 inf Qu
- 22/²Q111
(23)
QESO+
where || the Frobenius norm and SO is the special orthogonal group.
Formula in equation 22 and the first equality of equation 23 are direct consequences of the fact that a Riemannian submersion. To obtain the second equality of equation 23 we first notice that optimizing on Q1, Q2 Ou is equivalent to optimizing on a single Q O₂ thanks to the invariance of the Frobenius metric w.r.t. the right action of Og. And second, that the infimum attained at (see (Bhatia et al. 2019 Equations 3 and 35))
=z= 1-
where T-M "/,33/>)//
the Monge map from Σ1 En (see (Malagó et al. 2018 equation 8)), and Q" has positive determinant and belongs to SOa-
Thus the closest element of the fiber ग = 1/2 is given by / T i.e. by left multiplying E//² by the Monge map T. This is more generally true for any representative of Σ;: Proposition 11. Let Σ,, En S T the Monge map from Σn to Ez. A E ( Then A2 : TA, is said to be aligned with respect to A₁, that is, it is the closest point in ग ((2) to A More precisely, we have
1. A₂ A (T- 11A) E Hor
2. Logg, (E2): dIT., ((T T)A1) ( -T- TJE, ET
21

---

Published as a conference paper at ICLR 2026
GLu
D
honzontal gobdeste
honiontal goodesie
A2
Id
4
X
Hora
id
Vf
Il
(pp))
(Id)
x Mpp)
Wasterstein grodesie
Prob([))
Wasserstein geadesic
Id
M
Figure 19: The Otto-Wassent Wasserstein geometry of (left) centered Gaussian distributions and (right) a.c. probability distributions. Figures inspired by Khesin et al. (202I).
3. (W2 E2) = Logs, De 1E = II(TT - I)A1
where Log is the Riemannian logarithm map, || 1iW
BW (---) and || gE
the Frobenius norm.
This means that to compute the Bures- Wasserstein distance between two covariance matrices Σ₁ and Eg, one can consider any representative A₁ in the fiber over Σ, compute the representative A2 of Σ₂ aligned to A₁ (using the Monge map) and finally compute the Frobenius norm of A2 Ai.
B.3 GEODESIC PARAMETRIZATION
There are two classical parameterizations for Wasserstein geodesics in the space of a.c. probability measures.
McCann's interpolation The first one, due McCann (997). is given between two probability distributions Si and μu,, and depends on the optimal transport map in equation 2 gradient of a convex function u, that is T Vu and
obtained as the
He ((1 t) id +tV do (id -t(Vu id)) t [0, 1].
(24)
Otto's geodesic The second one, exploiting Otto's fiber bundle geometry in Otto (2001), consists in writing a geodesic in the Wasserstein space as the projection of a horizontal geodesic the total space of diffeomorphisms. Such horizontal geodesic is line segment going through diffeomorphism y with horizontal speed f o v, where f is any smooth function (not necessarily convex). Therefore we get
(25) In this second expression, the bounds on the time depends on the function f. Indeed, for μ., to be geodesic, id + f needs to remain is the space of diffcomorphisms for given s, which means that id +HHsss needs be positive definite. Therefore, we get the following conditions depending on the minimum Amin and maximum Amax eigenvalues of Hess f:
= (p + vV = (id VVV (++P), s E (So. $1).
8 E (- 00, -1/mmmn) if Amax < 0, E (- 1/ Am/ +00) E 1/A / -1 11/Amin) if Amin < 0 < Amax-
if Amin > O,
(26)
clear that equation 24 is particular case of equation 25| where we choose P Mo and f = Vu - id. Conversely, one can write equation 25 under the form of equation 24| For given diffeomorphism & and function f. consider the geodesic given by equation 25 and set A4 = P#P. Assume that we are in the case where all eigenvalues of Hess f are negative, then must be in
00, -1/Amin Consider EJO, 10-1// 11/Amin, and define Jr: = P.- = (id + 'V Setting t s/ s" have that the geodesic between do and μn is written
Pr (id +1." V.) (id +VV id)) EE [0,
22

---

Published as a conference paper at ICLR 2026
for u(z) f f 1|x11//2. Now for any eigenvalue X of H the Hessian of f, we have
A Amn > -1/5" i.c. by the interval of definition of s*. This means that the Hessian H, = H + id is positive definite, which means that u is necessarily convex. The other cases work similarly.
The Caussian case Gaussian distributions means that for E GLd and X E Hor A such that ||X]| = 1, the interval of definition of geodesic depends on the invertibility of A + sx. In turn, the maximal interval of definition of E (8o, 81) is defined from the eigenvalues of X XA-1, through the same formula in
Transposing Otto's formulation in equation 25 to the case of geodesic between
equation 26
C LINEARIZED OPTIMAL TRANSPORT AND TANGENT PCA
In this section, we provide the definition of linearized Wasserstein distance and details on how to perform tangent PCA for both Gaussian distributions and general a.c. distributions. Tangent PCA is a widely used approach to compute PCA on the Wasserstein space, that consists in embedding probability distributions into the tangent space some reference measure p. and performing PCA in the tangent space with respect to the linearized Wasserstein distance.
C.1 THE CASE OF CENTERED GAUSSIAN DISTRIBUTIONS
We consider n covariance matrices E....., En and their Bures- sarrs-Wasserstein barycenter (or Fréchet mean) Σ, that is, the SPD matrix verifying (see Agueh & Carlier (2011))
Σ arg min EESt
EE,
(27)
The idea behind tangent PCA is to represent each data point by the corresponding tangent vector, given by the Riemannian logarithm map, in the tangent space at the reference point Σ, i.e.
{Logs C. TsStt.
(28)
Now, can lift the computations from the tangent space at Σ to the horizontal space at point in the fiber over Σ, say A: E1/2 by aligning all representatives to A, see Proposition 11 The key point that the tangent space Σ equipped with the Bures Wasserstein Riemannian metric isometric to Hor S.A equipped with the Frobenius inner product where we recall that S₂ is the space of symmetric matrices. This means that instead of performing PCA for the Bures- warrs-WWssersttin inner instead perform linear PCA their pre-images
product the tangent vectors in equation 28 by A see Proposition 11
{(T:- ( 1) C HOrA, where 7,
Tj is the optimal transport map from Σ E,, see Section B.2 Now, noticing that
(24) TrKKAA K,
(t(KKEK)) AKG. Sa.
we see that the space A equipped with the Frobenius inner product is itself isometric to Sa equipped with the Frobenius inner product weighted by Σ. Therefore, tangent PCA is performed through Euclidean PCA on the (centered) vectors (T. - 1= to the Frobenius metric weighted by E. Another way to see this is by noticing that the linearized Wassersessn distance BW,, 2.2 with respect to Σ is given by
in the vector space Sa, with respect
BW2s(En, E2):
||Togs E₁ Logs, 11 Bw
| (TT - 1/2 dogiva ((T2 - 1) Bw [[ 1) (T2 1)
= [[ - -T>)E1/11
where Bl denotes the norm associated to the Bures Wasserstein Riemannian metric in equation 18 is Otto's projection in equation [9] and we have used Propositions 2 11 Finally,
BW Σz): [| Logs BW =| Tells, (29)
where ||
denotes the Frobenius norm weighted by Σ.
23

---

Published as a conference paper at ICLR 2026
C.2 THE CASE OF A.C. DISTRIBUTIONS
Similarly, one can embed a.c. probability distributions ,1...., Mn into the (p) space at some a.c. reference measure p through the optimal maps v 7% in the Monge problem in equation 2 Then, the Wasserstein distance can be approximated by the linearized Wasserstein distance in Wang et al. (2013) given by
W2,p(11, 122 1 - Th
(30)
Note that previously mentioned, this metric induces distortions while the radial distances from p to any b4 are preserved, that |lid To (2=11 l1T - 7 (مم * W2(U1, u2). recent paper by Letrouit & Mérigot (2024) proved however, that under some assumptions, 22,p is bi-Hölder equivalent to W2, which indicates that the distortion effect can be controlled.
W2(p, vi), other distances are not
Then, denoting Pn the Wasserstein barycenter in Agueh & Carlier 2011 of PI..... Mas that is the solution of
Un arg min
Σ (2,,,))
(31)
langent PCA consists in performing classical PCA, see e.g. (T id) in the Hilbert space (1,))
Ramsay & Silverman (2002), of
D GEODESIC PCA FOR GAUSSIAN DISTRIBUTIONS
In this section, we present the proofs related geodesic PCA for Gaussian distributions and the implementation of our algorithm in this case.
D.1 PROOFS RELATED TO GPCA FOR GAUSSIAN DISTRIBUTIONS
We first prove the existence of mimimizers for the GPCA problems lifted to Otto's fiber bundle. Lemma 1. The GPCA problem in equation 12 for the first component admits global minimum.
Proof. First, let us define the set of normalized matrices IB: {X gdxd |||||| 1}. By denoting Amin (resp. Amax) the smallest (resp. largest) eigenvalue of XA- extending the geodesic t + A + t.X as far as possible (see Section B.3) means that the closed interval [fmin [max] is defined
for some fixed > 0 by {
00, Amin - E] if Amax < 0, 1/2aax + 8, +00) if Amin 0, 1/ + -, -1.-1Ammin E] if Amia
(32)
Amax-
Let us now consider the function
F: GLa x x (Ruxdgn R
2 Σ |A (txx))(txXxX 2//22012 = - (A,
where t, = (2, and P(A,x) R R is the projection operator that clips point t into
mimn m. which depends on A and X. Then the function F' is continuous on GLa x B x (Rdxdyn as composition of linear and continuous functions. Note that the function XX. + P(),x) (t:) continuous by eigenvalue continuity, see [i Zhang (2019). Additionally, the function F is coercive (see e.g. Zalinescu )2002) GLa x IB [RaxaTa Indeed, a diagonal {A
2//2,, for (A, GLa Rdxd] for some i E {1.... m). we have t, O. and therefore we have either (A,A,, Q:) = 0 if (0) = 0, or (A X, Q:) = x|| 1|2
otherwise. This would imply that 9: (A, (.) doesn't 't go to infinity when the norm (I(A, + co. However, in this case, we have gj (A, Q2)
oo when [,A, X,), - oo for any * i. Moreover, as P(A.x) (f,) is clipping, it won't play role in the coercivity. We conclude by the fact that the function (A, + X A A X is continuous, implying that the set of constraint {(A,X) E GL, x Rdxdd X A A = 0} closed and and SO₂ are compact. The optimization problem in equation 12 thus admits global minimum.
24

---

Published as a conference paper at ICLR 2026
Note that this result also applies for the second component in equation 13 and the higher order components.
Proposition 12 (Proposition main), Let n. GLa + SJ A + AAT and (1.))
be a solution of
inf F(A1, (QQ))))) = Σ PA - 2./2
subject to A GLad, X1 E Hor || QI .0 ESAa
Then there exist bmin: Lmax E IK such that the geodesic E: [mmin: bmax) + a(A + 1X1) in 5
minimizes equation 11
Proof. A horizontal geodesic in GLa is straight line going through base point A E GLa in the direction of horizontal vector X A (that we consider normalized, ic. || 1), i.c. t + A tX E GLa. Denoting min maax] the interval constructed in equation (32 which depends on the eigenvalues of XA- 11 we have that (T (X) (E(tmmtt is geodesic in the Bures Wasserstein sense, see Proposition and
min (T LE frals
min inf tE(tun, Q,ESO,
||4 MA++X- tX -
inf ||4 + (()x - 21/22 Q.ESO,
where t, = (2/ is the (orthogonal) projection time of Σ 2 onto the line t + A A+tX.
We therefore deduce that a of solution (A, X, (6.x((Q.. (Q,) of equation 12 defines proper geodesic (T(A + tX) fmm)) solution of problem in equation
Proposition 13 (Proposition 5 in main). Let v4 = = of) for i = 1....n be n univariate Gaussian distributions. The first principal geodesic component E [0, 1] - p(l) solving equation / remains in the geodesic space of Gaussian distributions for all E [0, 1].
Proof. Let Prob2(R) be the set of a.c. probability measures on R that have finite second moment, and Q the set of corresponding quantile functions:
(F v e Proby(R))
Q is the of increasing, left-continuous functions q : (0, + I. and a convex cone in L2 1)). the set of square- inttrr-ittegrable functions on [0, 1]. The mapping
t: F,1
(33)
defines an isometry between Prob (R) equipped with the Wasserstein metric, and Q equipped with the L metric (see e.g. Bigot et al. (2017)), that is, for any μ, v E Proba(R),
W2 v) |1 - 11 [=([011)-
The map & in equation 33 also defines isometry from the set of (univariate) Gaussian distributions the set of all Gaussian quantile functions G. This space g is the upper-half of the plane F spanned by the constant function 1 and the quantile function F of the standard normal distribution:
G= R 1RR1+RR R Fo FF: span (1, 1.
Now, consider n normal distributions 1/,..... ... Vn. and (p(t))re(o.1) the first principal geodesic com- ponent found by minimizing equation [ the
of squared residuals in Probz(RR.. Since μ is
Wasserstein geodesic in Proby(R) and & is isometry, the curve t ++ (())(t) F-1 is
((0, geoodescc in Q, i.c. line segment
t E (0, + F-1 ((8)
(1 F1 + +FF
((0) m())
25

---

Published as a conference paper at ICLR 2026
Since {1, Fi forms orthonormal basis of F, the orthogonal projection of this line segment F is given by
1. E [0, 1] - (FTee> 1)1 (Fatey Fo
which lies in G. To see this, we need to show that the following value is positive:
M(0) Fo S/¹ F (at Fo) (v)dy- xFFo' o (z)) -B()))((((())=XX(((())))))))))))
where X plt) and T' Fo of is the Monge map from p(t) the standard normal distribution. Since T is increasing, indeed have B(XT(X)) (see e.g. the proof of Theorem 2.2 in Schmidt (2014).
Finally, since (())) orthogonally projects from Q w.r.t the 1 metric and T defines isometry, we get that the geodesic μ orthogonally projects to geodesic I(μ) in the space of Gaussian distributions, w.r.t the Wasserstein metric. By the distance minimizing property of orthogonal projections, we know that the cost function in equation [ evaluated at I (H)) is no larger than its value at μ. Since 4 optimal, we get that μ T()) and y belongs the space of Gaussian distributions.
Proposition 14. Let Σ1, D2 two SPD matrices that are diagonalizable in the same orthonormal basis, i.e.
E1 = (of
0) pT and Σg =
(n3
PT.
0
PT.
b7
63
where P is orthogonal. Then BW?(E,, (E1, Σz) = (a) (a1-a2)² -2) (b1 b2) and thus the Bures- Wasserstein geodesic between Σn and Σg given by
Σ(t) ( t)a1 + tb1)
((1-t)ag - th,)3 tbb))apPP
0
PT, 05ts1.
Proof. This is a straightforward computation using equation
Proposition 15. Let us consider n = 2p covariance matrices Σ, = E(a, 0,) as defined in equation 16 where 0, = in in/n for 0. .......2 1. Then, the Bures- Wasserstein barycenter in equation 27 of these covariance matrices is given by Σ = (a + /42
Proof. Each pair of covariance matrices
aa² 3 = Po. (o²
b² Po.. and Ditp = Po.to/2DP #:40/2 = Po. (88
0
PJ.
a²)
are diagonalizable in the same basis, and so by Proposition
14 the geodesic from Eg to E,+p is
M(t) = Ps. ((1 ((11 l)a + 1b)²
((1 -1-11 1)b+ to)2 PT,, 0
In particular, the Fréchet mean is given by Σ = E1/2) = ((a + b)/2)3I. Since each pair of covariance matrices has the same Fréchet mean, the Fréchet mean of the whole set E,...., E,, also given by Σ.
Proposition 16 (Proposition in main). Let Σ E S with eigenvalues a² b² and E' = P&EPI where Po is the rotation matrix of angle 0. Then, denoting Σ ((a 2)/2)² I we have
BW3 E') BW2 1,2,5 (E,E) = (a (a-b) (a=b)= cos² 0 + O((a () la+ bb
(34)
Proof. Recall that the linearized Bures-Wasserseen Wasserstein distance at between E and E' is given by the distance between their images by the Riemannian logarithm map U: = Logs and U': = Logs in the tangent space at Σ, i.c.
BW, 2,5. (E,
|IU U' BW
26

---

Published as a conference paper at ICLR 2026
where ||-|| RW denotes the norm associated the wessessnssr Riemannian metric in equation 18
As in any Riemannian manifold, the true geodesic distance can be approximated by this linearized distance in the tangent space, corrected by the curvature (see e.g. Lemma in Harms al. (019))
BW2 E' - - UII 3
Rg(UUU',,,,0 O aw + -10' (35)
where Rg is the curvature tensor.
Recall from equation 18 that the Bures- Wasserstein norm of vector U is expressed in an eigenvector basis of the base point, here Σ. Since any basis is an eigenvector basis of Σ, it is convenient to choose that 2. which we can assume without loss generality to be the canonical basis. Thus write Σ D where D diag(a", b²) and E' PoDP, and the norm associated to the Bures- Wasserstein Riemannian metric is given by
Σ 186552
U
di + dj
where the d, are the eigenvalues of Σ, given here by d2 = ((a + b)/2)2. From Proposition 11 we have
U: = Logs = = (T (T- - I),
U²: = Logg = (T" 1)E + E(T' I)),
where
T: E 5-11"(("/255//))133355
1/2
a+b 2
T:: = 5-112(5112//5/11113333333
Po 12 pp
a a b
and easily get
a² 12 J,
U' b²
cos 28
sin 29
P
where PoJP
sin 28
cos 20
and J = =diag(1, -1). Thus after some computations we obtain
| U'| BW =aa- bl/V2, BW
(36)
BW, 2.8 (Σ = U U'I
BW = v2((a - b) sin
To compute the curvature tensor, we use the following formula from Thanwerdas 2022 Table 4.7)
Re(U,U", = Σ d,dg UV, Un JJ
d, d
where [A, B] AB BA is the Lie bracket of matrices, Uo and U/ are the only symmetric matrices
verifying the Sylvester equations U = UQ + Vo and U Un EU respectively. Since Σ is multiple of the identity, we easily get
b J. Un b
b PoJP
Up
a a b
and straightforward computations yield
3 (a b)a
Rs(UUU',UU,Y)
sin 20.
(37)
(a + b)²
Finally, putting together equation
35 equation
36 and equation 37 and obtain
- (a + b)²
6)
BW = (2, (2,
sin 0 cos 2 + O((a - b)%)
and dividing by the squared linearized optimal transport distance yields the desired result.
27

---

Published as a conference paper at ICLR 2026
D.2 IMPLEMENTATION GPCA FOR GAUSSIAN DISTRIBUTIONS
As described in Section 3 the first and second components geodesic PCA are respectively found
by solving the minimization problems in equation and equation 13 given by
The geodesic components
E,(1) - (A; + + (A, + for i= 1,1,2
where A1 E GLd and X, E HO are minimizers of equation and Ag E GLd and X₂ E Hor minimizers of equation| 13 The matrix T(42) is the crossing point through which all geodesic components intersect, see Figure 2 The higher order components are found in analogous way: for the k-th component, we search for horizontal segment the fiber over the intersection point (we parametrize it w.r.t. the previous position in the fiber, i.c. Aş- 1 for certain Rx-1- SOa) and the horizontal velocity vector X is orthogonal to the lifts of the velocity vectors of the previous component. Thus, the k-th component, k 2 3, solves:
As LX where Ag belongs to
inf F(Au, Xx, (Q1)
subject to
A = A R&-1 E SOO, X2 E Hor ||all =1 (X X Rg-1) = 0, 1 s -, eI....,QU SOd-
(38)
Following Huckemann al. (2100) and Calissano al. (2024), propose iterative algorithm to implement these components, that, for each component, alternates two steps:
(Step 1) minimization of the objective function (see equation X),
12 with respect to (Q4);=1 for fixed
(Step 2) minimization of the objective function F' with respect to (A, for fixed (Q+)
In dimension d = 2, any rotation matrix Q can be parametrized by a scalar angle 0 and both steps are solved using the Sequential Least Squares Programming (SLSQP) algorithm (see e.g. Ma et al. (2024) available on the scipy python library and given by Virtanen et al. (2020). In higher dimension, each minimization with respect rotation matrix performed using Riemannian gradient descent on SOa, relying on the Riemannian geometry of SOa induced by the standard Frobenius metric of the ambient space ddd In particular we use the exponential map implemented in the Python library geomstats developed by Miolane et al. (2020). More details the Riemannian geometry of SO, and the Riemannian gradient descent procedure can be found e.g. in (Boumal 2023 Sections 7.4 and 4.3).
Unfortunately, we cannot ensure the convergence the iterates of the proposed block alternating algorithm, classical arguments require uniqueness of the minimizer at each iterations proven in Powell (1773). This is unachievable our problem: the line with base point A and direction X E Hor A and the line with base point AQ and direction XQ E Hor AQ for Q E Oa project onto the same geodesic in the bottom space. However, regarding (Step 1), and thanks to Theorem 3.7 in Huang Weil (ei22022), we have for fixed (A, X) that the cost function f: (QI,..., Qa) F(A, (Q.) has the Riemannian convergence of the iterates towards an accumulation point thanks Theorem 3.14 in Zhou al.) (2024). The three assumptions in this theorem are verified in our case Assumption (3.5) (I.- Retraction Smoothness) is obtained because grad is Lipschitz, and Corollary 10.54 Boumal (2023); Assumption (3.7) (bounded from below) directly holds because f > 0; Assumption (3.8) (ndividual Retraction Lipschitzness) is verified thanks to Corollary 10.47 Boumal 2023).
Lojjskawwoasiiiiizzz property at any point of (Oa)a. Finally, we have the
Scalability of the algorithm Surely, the computational time of our algorithm for Gaussian distri- butions will increase with the dimension. However, the algorithm can be made less sensitive to the number of input covariance matrices by parallelizing (Step 2) of our algorithm, which consists in updating the orthogonal matrices (Q.) This would significantly reduce the overall computational cost of the algorithm. Also, we currently use the scipy toolbox to solve (Step 1), which could also be accelerated using more powerful optimization toolbox.
28

---

Published as a conference paper at ICLR 2026
E HYPERPARAMETERS
E.1 HYPERPARAMETERS SETTING
Hyperparameter
Value
dense MLP
fo architecture
- 128 - 128 - 128 - 128 - 1
ELU activation functions
Adam
step size 0.0005
f optimizer
31 0.9
3g = 0.999
dense MLP
vPe architecture
+ 128 - 128 - 128 - 128 - d
RELU activation functions
Adam
step size 0.0005
Pe optimizer
31 0.9
32 = 0.999
Adam
step size 0.005
4 optimizer
31 0.9
82 = 0.999
batch size
1024
number of gradient steps first component
120,000
number of gradient steps second component
200,000
A
1.0
1.0
Table 1: Hyperparameters used across all experiments.
All experiments were conducted on single V100 GPU with 32GB of memory, using a shared set of hyperparameters detailed in Table The same hyperparameters are used for computing both the first and second geodesic components, except for the number of gradient steps (see Table I), which is increased for the second component. This is likely due to the additional complexity introduced by the intersection and orthogonality constraints enforced through regularization. Both f& and yo are implemented standard multilayer perceptrons (MLPs) with four hidden layers of width 1 128. We use ELU activation functions in fø because its gradient is used to parameterize a transport map in our formulation, and ELUs are commonly employed in such settings. The Sinkhorn divergence Se is used in the loss function as a surrogate for the squared Wasserstein distance to compute the geodesic components. The regularization parameter & must be adapted to the scale of the data; we set as = 0.01 Em,r">v; x'||2, where the expectation is approximated via Monte Carlo using the current minibatch samples. Note that setting & this way is the default configuration in the OTT-JAX library. For computing the second geodesic component, we fix the regularization coefficients Ao and Àz to 1.0. which we found to be robust across all experiments. While increasing them (e.g., to 10.0) typically yields similar results, excessively large values may degrade performance. Conversely, if these regularization terms are too small, the algorithm tends recover the first component the second, due to its lower cost. In practice, we monitor the regularization terms during optimization to ensure they decrease sufficiently relative to their initial values, confirming that the optimization effectively optimize the intersection and orthogonality constraints. To determine the hyperparameters in Table [
we performed a grid search over the optimizer learning rate for the t; in 5e 1e 5e -3 1c and over the regularization coefficients o and Az in 0.1, 1.0, 10.0, 100.0. We found that setting both regularization terms to 1.0 consistently yielded good performance across all experiments, see Section H.2
29

---

Published as a conference paper at ICLR 2026
Note yp parameterization. Note that although yp is theoretically required to be diffeomorphism in Otto's parameterization of geodesics (equation 9). we parameterize it using a simple MLP. Initially, experimented with normalizing flows to ensure invertibility, but observed that standard MLP yielded similar results. In Otto's ' geodesic framework, ip serves to modify the reference measure and define the measure at & O along the geodesic. If yo is not diffeomorphism and the pushforward P#P is not absolutely continuous, the resulting geodesic becomes degenerate, which may hinder optimization of the loss equation in equation In practice, however, we found that the MLP ye reliably produces absolutely continuous measures, which is sufficient for our method.
E.2 IMPACT OF THE REGULARIZATIONS GPCA
For the estimation of the second GPCA component, we introduce two regularization terms, S.0... +1 +2 (mter) and fe fe.((e))) with their associated regularization co- efficients λ and λo. The first term enforees that the two components intersect, while the second ensures that the components remain orthogonal. Experimentally, we observe that setting both co- efficients A5 Xo 1.0 robustly enforces these constraints across all experiments while still producing meaningful principal components. Conversely, if these regularization terms are too small, the algorithm tends recover the first component as the second, at it gives the lowest cost. In practice, we monitor the regularization terms during optimization to ensure they decrease sufficiently relative to their initial values. This permits confirm that the optimization effectively optimize the intersection and orthogonality constraints. This section aims at quantifying the impact of the two regularizing coefficients \ and λo on the computed geodesics. We focus on the 3D point-cloud experiments with lamps.
Vitter:
E.2. ORTHOGONALITY REGULARIZATION
In this part, we set the regularization term x to 1.0 and compute GPCA for different values of λo. The resulting second component is shown in Figure 20 The GPCA cost of this component, as defined in equation 15 together with the quantity measuring the orthogonality between components. fg V f (Po,)) )). are reported in Table 2 The quantities reported in Table 2 are estimated on batches of size 2048. The variance computed over 100 runs for the orthogonality measure and 5 runs for the GPCA cost. Note that each run of the orthogonality estimation already involves computing 100 Wasserstein distances, since we have 100 point clouds.
Note that the GPCA cost of the second component should be compared with that of the first component, which is around 3.0. Table 2 shows that for low values of λo (i.e., 0.001 and 0.01), the orthogonality quantity is large, and the recovered "second" component is in fact identical to the first component, as illustrated in Figure 20 This is also reflected in the GPCA cost (see Table ( which matches the one of the first component. For higher values of do (0.1, 1.0, 10.0, 100.0), the algorithm successfully recovers a distinct second component. For the highest value (i.e., Xo = 100.0), a loss of performance is observed.
Orthogonality: (VVS) (Vy(()))))
GPCA cost (second component)
0.001
0,909 + 0.005
2.96 + 0.01
0.01
0.811 + 0.008
3.00 + 0.01
0.1
2.1 x 10-3 + 3 x 10
5.75 + 0.02
1.0
3.1 x 10-4
4 + 6.5 x 10-5
5.76 + 0.02
10.0
2.2 x 1
5 x 10
5.89 + 0.02
100.0
1.1 x 10
+ x x 10
5.99 + 0.02
Table 2: Orthogonality regularization value and second-component loss for different values of λo
30

---

Published as a conference paper at ICLR 2026
0.001
Ào = 0.01
-2
-1
0
1
2
-3
-2
-1
1
2
40
Ag = 1.0
-4
-2
0
2
4
-3
-2
-1
0
1
3
10.0
2
100.0
-2.
-15
-1.0
-0.5
0.0
-2.
-1.5
-1.0
-0.5
0.0
Figure 20: Empirical distributions sampled uniformly along the geodesics associated with the second GPCA principal component for different values of the regularization coefficient λo In all experiments, the other regularization coefficient is fixed A 1.0.
E.2.2 REGULARIZATION ON THE INTERSECTION OF THE GEODESICS
In this part, we set the regularization term Åg to 1.0 and compute GPCA for different values of A The second component displayed in Figure 21 the GPCA cost of this component, well the quantity measuring the intersection the components, I[Ee, w E02. interr ) are reported in Table 3 The quantities reported Table are estimated on batches of 2048. The variance is computed over 100 runs for the intersection measure and 5 runs for the GPCA cost.
We observe from the recovered geodesics in Figure 21 that this regularization term plays less significant role than the orthogonality term. Moreover, Table shows that increasing A does not affect negatively the GPCA cost of the recovered component.
A 0.001
4, 0.01
-4
-2
0
2
33
-2
-1
0
1
2
4 = 0.1
1.0
-4
-2
0
4
-3
-2
-1
0
1
4, = 10.0
2, = 100.0
-3
-2
-1
1
2
-0.5
0.0
0.5
1.0
1.5
2.0
Figure 21: Empirical distributions sampled uniformly along the geodesics associated with the second GPCA principal component for different values of the regularization coefficient Å₁. In all experiments, the other regularization coefficient fixed at λo = 1.0.
31

---

Published as a conference paper at ICLR 2026
X
Intersection: [(5 gi S02. iiter: )
22
GPCA cost (second component)
0.001
6.5 x 10-2 + 1 x 10-3
5.76 + 0.02
0.01
2.3 x 10-2
4 x 10-4
5.77 + 0.02
0.1
2.7 10-3
x 10-4
5.77 + 0.01
1.0
1.0 10-3
+ 6 X 10-5
5.76 + 0.02
10.0
7.2 x 10-5
3 x 10-6
5.74 + 0.02
100.0
3.1 x 10-3
1 10-6
5.91 - 0.02
Table 3: Squared Euclidean distance between Ç1 (finter) and E2 (femee) and second-component loss for different values Xz-
E.2.3 SCALABILITY OUR GPCAGEN ALGORITHM
For general distributions, there are two types of "scaling" that can affect the algorithm:
1. Number of probability measures (n): The number of measures v; directly determines the iterations of the inner loop in Algorithm I (line 3). Consequently, the training time scales linearly with
2. Dimension of the space (d): As the dimension of the space in which the v; lies increases, the main challenge consists in accurately estimating the maximum and minimum eigenvalues that the Hessian of can take. As discussed with reviewer oUM,, in high dimensions, becomes necessary to use algorithms that avoid computing the full Hessian and instead rely on matrix -ectox vectorr products, such as the LOBPCG algorithm Duersch et al. (2018). Further- more, rather than relying solely on the samples in the training batch, an ad adversarial approach would be needed to track the eigenvectors corresponding to the worst-case eigenvalues.
USE OF LARGE LANGUAGE MODELS (LLMS)
LLMss were used only assist with polishing the writing; all research ideas, experiments, and analyses were conducted independently by the authors.
32


Reviews:
[{"review_id": "c340AgHtlP", "reviewer": "Reviewer_oUMT", "summary": "This paper introduces a method for exact PCA in the Wasserstein space of probability measures P_2(R^d) under the Wasserstein metric. Instead of performing tangent PCA, which has been an established linearization technique for performing PCA on curved spaces (most heavily used in Riemannian settings) or geodesic PCA (Huckemann et al. Sommer et al.) that have tackled this problem in general Riemannian settings, this paper proposes two methods, i) for geodesic PCA for Gaussian measures, and a ii) neural network approach for parameter optimization for absolutely continuous probability measures. The authors show experimental evaluations on point clouds and image color distributions (mostly qualitative).", "questions": "For experimental results, the authors compare their method against TPCA. How will their method compare with Bigot et al. 2017 (for 1D case)?", "limitations": "", "rating": 10, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The work shows a method to compute principal modes of variation in datasets of probability measures, specifically using the Wasserstein geometry. For Gaussian measures, the method leverages the Bures-Wasserstein geometry and lifts computations to the space of invertible matrices, providing exact geodesics as principal components.\n\nThis is a significant contribution over earlier methods which have used linearized Wasserstein distances (Wang et al. (2013) and Boissard et al. (2015)), have approximated Geodesic PCA in 2D (Cazelles et al. (2018)), and the approach by Seguy & Cuturi (2015), which relied on  generalized geodesics to approximate PCA. \n\n\nFor general absolutely continuous distributions, the reliance on neural network parameterizations of Wasserstein geodesics is clever. \n\nFor the case of centered Gaussians, the application Bures-Wasserstein geometry and the decomposition via invertible matrices is well-motivated and theoretically grounded.", "weaknesses": "The block alternating algorithm for Gaussian GPCA is not guaranteed to always converge to a unique minimum due to non-uniqueness in the problem geometry (the authors acknowledge this).\n\n\nIn the general case, one needs to verify the eigenvalues of the Hessian at each step during the Otto geodesic update. This may be computationally expensive. \n\nWhile the neural network implementation facilitates computational tractability, the construction of geodesics needs further tuning and learning from large examples and may have scalability issues. \n\nThe comparisons and experimental results are all qualitative and thus the practical application of this method is in question. However, this is an extremely minor point, not important for the publication of the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "zlYVkxvEOk", "reviewer": "Reviewer_2t9m", "summary": "This paper introduces an algorithm for geodesic PCA using optimal transport ideas. This is done by generalizing a certain formulation of Euclidean PCA, giving an optimization problem for the first principle component involving geodesics in Wasserstein space. Further principal components are defined using orthogonality. There are two versions presented: for Gaussians, and for more general distributions relying on neural network approximation. Evaluations are performed on synthetic Gaussian-type data, and point clouds.", "questions": "The overwhelming majority of my concerns with this paper concern the comprehensiveness of evaluations and specifically comparisons with various baselines like PCA on PointNet latent space. Please address those comparisons, or argue why they are not appropriate.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I liked the following:\n* **Interesting problem.** Generalizing PCA to spaces of probability measures seems to be a generically useful tool, since comparing distributions is a central task throughout machine learning which recurs in many situations.\n* **Technically sound - especially for Gaussian distributions.** The approach involves the Bures-Wasserstein geometry and relationships between certain matrix groups, and makes it easier to see\n* **Easily provides use cases beyond what the authors have illustrated.** One way one can think about PCA is that it provides a very crude form of manifold learning, where the manifold is a subset of the ambient Euclidean space. This can be used to define for instance an outlier score, by comparing distance between each datapoint and the PCA subspace. The authors do not consider this, but the fact that it is easily possible reveals the general value of the technique.", "weaknesses": "I am worried about the following:\n* **Use of regularization in neural network objectives.** In particular, using regularization to enforce geometric constraints is much weaker than incorporating them as a hard constraint via a clever parametrization. In practice, I suspect the different directions do not end up orthogonal, and it would be helpful to quantify how much this is a problem in practice, and how sensitive it is to hyperparameter tuning. I did not see an experiment directly addressing this.\n* **Experiments: Gaussian case should have at least *some* kind of non-synthetic evaluation.** Gaussian experiments are limited to toy datasets designed to highlight differences with ordinary PCA. While this is certainly necessary, it would have been better to also include some kind of non-synthetic datasets where the Gaussian mean and covariance is estimated empirically. For example, in continuous action imitation learning and RL, one often considers policies consisting of neural networks that learn a mean and covariance. In this domain, the authors' form of PCA could be used for things like outlier detection. I am not suggesting the authors necessarily need to add this specific example, more that it's important to have *something* non-synthetic for the Gaussian case since it's a central part of the paper.\n* **Experiments: why is there a theorem in the experiments section?** Theory should be presented separately from evaluation, this part should be moved to some part in main body.\n* **Experiments: Gaussian case is missing obvious baselines.** In 5.1, a natural baseline is to compare against is Euclidean PCA in the space of Cholesky coefficients of the covariance matrix. How would this produce different results compared to Euclidean PCA with respect to the SPD cone coordinates, or the proposed method?\n* **Experiments: non-Gaussian case is missing obvious baselines.** For the point cloud data presented, the most obvious way I can think of to do PCA is take some kind of off-the-shelf model like PointNet, and perform PCA in its latent space, and present the results. While this is certainly a much heavier method, I think some kind of visual comparison needs to be included in the main body, as this is the single most obvious performance question that readers might have.\n* **Experiments: no evaluation of downstream applications of PCA like outlier detection.** The authors do not evaluate how good the principal components produced are for any downstream purpose, instead they are just plotted to show their properties. I would like to see at least some kind of quantitative evaluation. For instance, this could be considered for outlier detection, to compare how well linear vs. the proposed geodesic PCA performs. Non-linear PCA has been studied for this specific purpose before, see for instance \"Quadric Hypersurface Intersection for Manifold Learning in Feature Space\" by Pavutnitskiy et al. which includes full experiments on this. Note that I am not suggesting the authors implement this specific downstream application of PCA, instead I think it is important to have at least some kind of evaluation of this kind.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0nxGe29lhb", "reviewer": "Reviewer_F5cW", "summary": "The paper introduces a framework for Geodesic Principal Component Analysis (GPCA) of probability measures under the Otto-Wasserstein geometry. It provides exact solutions for Gaussian distributions by lifting computations to the space of invertible matrices and extends the method to general absolutely continuous measures using a neural network-based parametrization of Wasserstein geodesics.  Experiments on Gaussian data, 3D point clouds, and images show that captures meaningful nonlinear modes of variation and highlight the difference to tangent PCA.", "questions": "- Intro: It would be nice to list some concrete examples on the shortcomings of using $L_2$ for probability densities\n  - It would be nice to show some concrete examples for prop. 6 in terms of measures and vector fields \n  - Problem 38: Why not a global optimization instead of the block optimization?\n  - In the paper, one component is computed after another. Could one find multiple components at once? Is it equivalent? Especially for the neural network part, I do not think this should be hard to implement?\n  - What happens for Gaussian distributions with varying centers? Is there still a nice reduced model?\n  - Typesetting: please use \\colon instead of : when introducing functions\n  - Can geodesics be extended in the Otto formulation beyond t_min/max? E.g. by choosing a new horizontal vector?\n  - Language: I believe the correct usage is \"to lift to \\[a space\\]\" instead of \"to lift in \\[a space\\]\"\n  - Fig. 3 and others: I believe it would be more instructional to show the cone as a volumetric object. This might be harder a render, but maybe the authors could give it a try.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "**Exposition:**\nI found the paper to be very well written. It has a common thread running through it that makes it easy to follow the story. Thus, I could read it in one go and understood all the core ideas. Furthermore, I think all the necessary information is included in the paper needed to reproduce the method and the experiments. The division of information between main text and appendix is also sensible.\n\n**Novelty:**\nI think the introduced method is novel and advances the state-of-the-art in measure-based (geometric) data analysis. Thus, I consider it to be a contribution worthy of being published at ICLR. The distinction to similar methods (esp. TPCA) is clearly made and well-explained.\n\n**Theory:**\nThe included theory is sound and reasonable. To me, all the included theoretical results clearly connect to the overall story of the paper and are well explained. \n\n**Reproducability:**\nWith the provided information, the method (esp. the ac. formulation) should be easy to implement. I liked especially that pseudo-code (Alg. 1) was provided, as this always makes things easier to understand. Furthermore, all the experiments are well documented to level where one could reproduce them.", "weaknesses": "**Scalability:** I am a bit worried about the scalability of the method. All examples are conducted on a small scale with at most two components. Thus, the paper leaves the gap what would happen for larger datasets and what kind of resources the method requires in such a scenario. It would be great if the authors could discuss this in the paper and also illuminate if it is, indeed, a problem.\n\n**Practical applications:**\nThis ties into the second weakness I see with the paper: a lack of practical applicability. As mentioned above, the examples are small and restricted to synthetic examples. Thus, I do not see on which applications the paper could have an (immediate) impact. As I see this as a more of a theoretical/modelling work, I do not consider this as a major issue. However, it would be nice if the authors could include a discussion what kind of applications they envision their method could be useful for (down the line) and how they image their ideas could impact other machine learning work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "k8fI7sz6hY", "reviewer": "Reviewer_Pw4V", "summary": "This paper concern with the problem of Geodesic Principal Component Analysis (GPCA), i.e., given $n$ probability measure $\\nu_1,\\ldots,\\nu_n$, we gradually construct a sequence of geodesics $\\mu_1,\\ldots,\\mu_k$, known as the principle components, such that   \n\\begin{equation*}\n    \\mu_i = \\inf_{t \\mapsto \\mu(t) \\text{ geodesics}} \\sum_{i=1}^n \\inf_{t_i} W_2^2(\\mu(t_i),\\nu_i), \n\\end{equation*}\nand the subsequent principle components satisfy the same cost function with additional constrains about orthogonality with previous principle components. To resolve this problem, the authors leverage tool from Otto-Wassernstein geometry. They propose two algorithms to tackle two different settings: centered Gaussian distributions and a.c. probability measures: for first setting, the problem reduces to the matrix quadratic optimisation, for the second setting, by choosing a Gaussian reference measure, the problem reduces to optimise a function with neural network. Then, they validate their algorithms using toy examples (for Gaussian distributions) and some real dataset (MNIST, 3D point cloud, etc. )", "questions": "1. As pointed out in the weakness section, I miss the point that how can we \nuse the proposition 6, 7, 8 in our algorithms. \n\n2. Can we use this algorithms for some LLMs model or other complicated image tasks? \n\n3. Can you point out an example that the global minimum point in Lemma 1 is not unique?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "$\\bullet$ Originality: The problem appears to be novel and applicable to the real dataset, as the authors demonstrate in their experience. \n\n$\\bullet$ Quality: Firstly, the authors provide both theoretical intuition and experimental verification for the algorithm. \n\n$\\bullet$ Clarity: The article is easy to follow in general. \n\n$\\bullet$ Significance: The algorithm provide a method to analysis the principle components of data. In the real data, they relate to some important aspects of data, such as brightness, shape, etc.", "weaknesses": "I think that the article does not provide a quite clear theoretical verification: despite of the fact that the authors provide a nice intuition from Otto-Wassernstein geometry, they do not provide the proof of convergence of these two algorithms. For example, in Lemma 1, they just point out that the loss function admit one global minimum and do not show that the algorithm converge to the, at least, stationary point. In addition, for the case of centered Gaussian distributions, as the author admit, it is unknown that the optimal solution is still Gaussian distribution in general distribution space. Another aspect is that some of the theoretical results seems to not useful in the analysis of these two algorithms. For example, I am not sure that how we use the proposition 6, 7, 8 in our algorithms. \n\nMinor suggestion: \n- Line 32: accound → account\n\n- Equation 26: inconsistent notation, the author use \"]\" to denote open bracket \"(\" and \"[\" to denote close bracket, while in other formula they use the \"(\" and \")\".\n\n- \"hessian\" should be \"Hessian\"\n\n- Line 180: GL should be italic \n\n- The proof for coercivity is not clear (among the notion of coercivity, which one that you use, and why the norm $\\|()\\|$ goes to infinity implies the coercivity).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)