On the convergence of graph Laplacians with a symmetric divergence Liane Xu∗
arXiv:2607.05892v1 [stat.ML] 7 Jul 2026
Abstract When analyzing a manifold learning algorithm for data lying on a smooth, compact, connected Riemannian submanifold (M, g) of Rd , a key estimate for the geodesic distance dg is that there exists K > 0 such that 0 ≤ dg (p, q)2 − ∥p − q∥2 ≤ Kdg (p, q)4 for all p, q ∈ M. We observe that more generally, when M is equipped with a smooth symmetric divergence D satisfying a non-degeneracy condition and g is given by gp := 21 Hessp (D(p, ·)) for all p ∈ M, there exists K > 0 such that D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4 for all p, q ∈ M. We demonstrate that this is sufficient for the pointwise convergence of graph Laplacians constructed with D and discuss examples where D is given by the Sinkhorn divergence on a family of probability measures parametrized by a manifold.
Contents 1 Introduction 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 2 3
2 Preliminaries 2.1 Construction of the graph Laplacian . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Sinkhorn divergences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 5 7 8
3 Approximations of the geodesic distance and convergence of the graph Laplacian 10 4 Application to Sinkhorn divergences 11 4.1 Regularity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.1.1 Particle system parametrized by a manifold . . . . . . . . . . . . . . . . . . . 12 4.1.2 Smooth deformations of a measure . . . . . . . . . . . . . . . . . . . . . . . . 13 4.2 Non-degeneracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5 Examples 19 5.1 Rotating particle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5.2 Dilation of two particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 ∗ Program in Applied and Computational Mathematics, Princeton University
Key words and phrases: manifold learning, graph Laplacian, Sinkhorn divergence
1
6 Conclusion and future work
28
A Proofs for Section 3 34 A.1 Pointwise convergence of the graph Laplacian (Theorem 3.1) . . . . . . . . . . . . . 36 B Proofs for Section 4
41
C Additional details and plots for Section 5
46
1
Introduction
1.1
Motivation
A popular assumption in data science is that high-dimensional data lie on or near a lower-dimensional manifold M [11, 33, 79]. However, in certain applications — for example, when trying to understand the different conformations that a molecule can take [87] — we may want to model the underlying data generation process using a family of probability measures {µx : x ∈ M} ⊂ P(X) supported on some compact X ⊂ Rd , parametrized by an unknown manifold M. Specifically, suppose we are given µx1 , . . . , µxN , where x1 , . . . , xN drawn i.i.d. from a probability distribution on M. Under what assumptions can we learn something about M, and how can we learn it? One option is to try to embed µx1 , . . . , µxN into a low-dimensional Euclidean space using a manifold learning algorithm such as Isomap [79], Laplacian eigenmaps [10, 11, 12] or diffusion maps [27]. More concretely, each of the aforementioned algorithms only requires a measure of dissimilarity to compare data points. For data that lie in a Euclidean space, the Euclidean distance is often used. For our problem, a natural approach would be to use one of these algorithms, but to compare µxi , µxj for each i, j ∈ [N ] using a metric or divergence defined on probability measures. Depending on the application (e.g., [53, 87]), one may want to use a metric arising from optimal transport, such as the Wasserstein-2 distance. Manifold learning with the Wasserstein-2 distance has also received interest [42, 61] because the Wasserstein-2 space has a formal Riemannian structure [62]. However, since this structure is only formal, it is much more difficult to work with than a bona fide finite-dimensional smooth Riemannian manifold. Perhaps unsurprisingly, regularity assumptions can simplify the analysis of manifold learning algorithms. A key estimate for a smooth, connected, compact Riemannian submanifold (M, g) of Rd is that there exists K > 0 such that ∥p − q∥2 ≤ dg (p, q)2 ≤ ∥p − q∥2 + Kdg (p, q)4
(1)
for all p, q ∈ M, where dg denotes the geodesic distance on (M, g) [14, 17]. Now suppose that we are instead only given a smooth, compact, connected, orientable manifold M, not necessarily embedded in Euclidean space. In previous work, we found that if M is embedded in a metric space (X, d), d2 is smooth on a neighborhood of the diagonal ΞM of M × M and gp :=
1 Hessp d2 (p, ·) 2
is non-degenerate for all p ∈ M, then a similar inequality to (1) holds, with the Euclidean distance replaced by d [86]. Using this, we can prove the pointwise convergence of graph Laplacians — a key object in the Laplacian eigenmap and diffusion map algorithms — following a similar argument as 2
the Euclidean setting [12, 14]. Separately, Arias-Castro and Qiao obtained results for Isomap under slightly different conditions [4]; although the analysis of Isomap does not require the fourth-order error in (1) [17], as noted in [4], assuming some regularity of d2 around the diagonal ΞM does simplify the analysis. This current work is motivated by the problem of understanding manifold learning algorithms when the Sinkhorn divergence is used to compare µxi , µxj for each i, j ∈ [N ]. The Sinkhorn divergence, a debiased version of entropy-regularized optimal transport [34, 39], is popular for several reasons. For one, it can be approximated quickly via Sinkhorn’s algorithm [29, 74]. In addition, for compact X, it is differentiable on P(X) × P(X) and convex with respect to each of its inputs [34], and its sample complexity is given by the parametric rate (as opposed to that of the squared Wasserstein-2 distance, which suffers from the curse of dimensionality) [38, 55]. However, since the Sinkhorn divergence itself is not the square of a metric [49], we cannot apply our previous results from [86].
1.2
Contributions
The contributions of this paper are twofold. First, we observe that for a smooth, symmetric divergence D : M×M → [0, ∞) under a non-degeneracy condition (see Assumption 2.2 for details), we can use the symmetry of D to obtain the estimate D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4 ,
(2)
for all p, q ∈ M, where K > 0 is a constant uniform across all p, q ∈ M and dg is the geodesic distance on (M, g) with gp := 12 Hessp (D(p, ·)). From here, the pointwise convergence of graph Laplacians constructed with D can be proven mutatis mutandis from the case when (M, g) is a compact Riemannian submanifold of Euclidean space and the graph Laplacians are constructed with the squared Euclidean distance [12, 14, 27]. Secondly, since the Sinkhorn divergence has nice regularity properties in several situations of interest, we can apply our results to these situations. Specifically, in Section 4, we will study the regularity of the Sinkhorn divergence and the non-degeneracy condition under the following two scenarios: Assumption 4.1 (Particle system parametrized by a manifold). Let M be a smooth, connected, compact orientable m-dimensional manifold with m ≥ 1. Assume that we have a smooth embedding ι = (ι1 , · · · , ιn ) : M → P (Rd )n such that ιj (x) ̸= ιk (x) for all x ∈ M and j ̸= k. Fix some n p1 , . . . , pn > 0 such that j=1 pj = 1, and define F : M → P(Rd ) by F : M → P(Rd ) n X x 7→ pj διj (x) . j=1
Moreover, assume that F (x) ̸= F (y) for all x ̸= y. Assumption 4.2 (Deformations of a measure). Fix any k ∈ Z>0 . Let X = Ω be the closure of a bounded open convex subset Ω ⊂ Rd , and let U ⊂ Rm be an open set. Let µ ∈ P(X). Suppose Ψ : U × X → X is a continuous function and that its first k partial derivatives with respect to the first variable ∂kΨ : U × X → Rd (29) ∂pi1 ∂pi2 · · · ∂pik 3
exist and are continuous on U × X for all i1 , · · · , ik ∈ [m]. Define F : U → P(X) p 7→ Ψ(p, ·)# µ. We will sometimes also write µp for F (p) = Ψ(p, ·)# µ. Assume that µp ̸= µq for all p ̸= q. Both assumptions are motivated in part by examples in existing literature. For example, one can construct an idealized model of a protein by treating each atom as a Dirac delta or point mass (e.g., [30, 72] and references therein), upon which Assumption 4.1 can be viewed as an idealized model for a protein whose conformation space is parametrized by a manifold. Although two structures are often considered equivalent if one can be transformed into the other via rotations and translations (see [30] and references therein for ways to handle this), we will not consider this issue in this paper. Assumption 4.2 is similar to [42, Section 2.2], the "horizontal perturbations" in [49, Remark 3.2] and the assumptions of [67] — each considers a family of measures defined by deforming a template measure — though we use different regularity conditions. We start with Assumption 4.1, as under this assumption, we can illustrate the main ideas without functional analysis. The argument under Assumption 4.2 is similar but involves more analysis. The rest of this paper is organized as follows. In Section 2, we provide some background on graph Laplacians and Sinkhorn divergences and discuss related work. In Section 3, we discuss (2), i.e., that one can approximate the induced squared geodesic distance with the given symmetric divergence with a fourth-order error, and the pointwise convergence of graph Laplacians. In Section 4, we focus more specifically on the regularity of the Sinkhorn divergence and the non-degeneracy condition under Assumptions 4.1 and 4.2. Lastly, we study two examples in Section 5 and conclude in Section 6.
2
Preliminaries
Throughout this paper, unless otherwise specified, smooth means C ∞ , and we will follow the notation and assumptions in Table 1. When working with a smooth Riemannian manifold (M, g), we assume that it is equipped with its Levi-Civita connection and that the Laplace-Beltrami operator is given by ∆g f := −divg (gradg f ) for f ∈ C 2 (M). Table 1: Notation and assumptions for this paper. M m h D ε N L(ε,N ) β OTβ Sβ P(X)
smooth, connected, compact orientable m-dimensional manifold, m ≥ 1 dimension of M function satisfying Assumption 2.1, used in the construction of the graph Laplacian divergence satisfying Assumption 2.2 bandwidth parameter for the graph Laplacian; see Section 2.1 number of samples discrete Laplacian, defined in Section 2.1 regularization parameter for the Sinkhorn divergence; see Section 2.2 optimal transport with entropy regularization β > 0; see (9) Sinkhorn divergence with regularization β > 0; see (10) the space of Borel probability measures on X ⊂ Rd
4
2.1
Construction of the graph Laplacian
We begin by describing the construction of the discrete Laplacians that we will use throughout this paper. Fix any ε > 0. We assume that we have functions h : [0, ∞) → [0, ∞) and D : M × M → R satisfying the following conditions: Assumption 2.1 (Assumptions on h (following [24, 25])). Let h : [0, ∞) → [0, ∞) be a continuous function such that h ∈ C 2 ((0, ∞)) and that h, h′ , h′′ are exponentially decaying, i.e., there exist c, c0 , c1 , c2 > 0 such that h(k) (x) ≤ ck e−cx for all x > 0 and k = 0, 1, 2. Assume also that h is not identically 0. Assumption 2.2 (Assumptions for the divergence D). Let D : M × M → R be a continuous function satisfying the following conditions: 1. D(p, q) ≥ 0 for all p, q ∈ M; 2. D(p, q) = 0 if and only if p = q; 3. D(p, q) = D(q, p) for all p, q ∈ M; 4. D is smooth on an open neighborhood of the diagonal ΞM := {(p, p) : p ∈ M} ⊂ M × M; 5. For all p ∈ M,
1 Hessp (D(p, ·)) (3) 2 is positive definite, so (M, g) is a Riemannian manifold. Although the Hessian typically depends on a Riemannian metric, for each p ∈ M, f (x) := D(p, x) attains a minimum at p, so for all V, W ∈ Tp M, Hessp f (V, W ) = V (W f ) gp :=
is well-defined independent of any Riemannian metric on M. Remark 2.1. In information geometry, for a smooth manifold M, a smooth function D : M×M → R satisfying 1., 2., and 5. in Assumption 2.2 above is known as a divergence; symmetry (condition 4) is not required [2, Definition 1.1]. Now suppose that x1 , . . . , xN are drawn i.i.d. from a Borel probability measure P on M. Define Kε : M × M → [0, ∞) via D(x, y) Kε (x, y) := h (4) ε for each x, y ∈ M. We now define the graph Laplacians that we will use: Definition 2.1 (Graph Laplacian). Construct a graph with {xi }i∈[N ] as vertices, and set the edge weight from xi to xj to be (ε,N ) := Kε (xi , xj ) Wij (5)
5
for each i, j ∈ [N ]. The degree matrix D(ε,N ) ∈ RN ×N is defined as (P (ε,N ) i=j (ε,N ) k Wik Dij = . 0 i ̸= j The (unnormalized) graph Laplacian is then given by (6)
L(ε,N ) := D(ε,N ) − W (ε,N ) . We also define a closely related operator L(ε,N ) via L(ε,N ) f (x) :=
N X
Kε (x, xj )(f (x) − f (xj ))
(7)
j=1
for any function f : M → R. Indeed, observe that for any f : M → R, (ε,N ) L f (x1 ) f (x1 ) .. L(ε,N ) ... = . . L(ε,N ) f (xN )
f (xN )
As an example, if there is a smooth embedding ι : M → Rd , then the squared Euclidean distance D(p, q) := ∥ι(p) − ι(q)∥2 satisfies the conditions in Assumption 2.2, and g in (3) is the Riemannian metric on M from restricting the canonical Riemannian structure on Rd to ι(M). We will sometimes refer to this as the Euclidean setting. This setting is a standard assumption for manifold learning: that our data lie on a compact Riemannian submanifold of a (potentially high-dimensional) Euclidean space. In this case, one well-known dimensionality reduction algorithm is Laplacian eigenmaps [10, 11, −1 (ε,N ) 12], which uses the spectral decomposition of the normalized graph Laplacian D(ε,N ) L . (ε,N ) −1 (ε,N ) (ε,N ) −1/2 (ε,N ) (ε,N ) −1/2 Since D L is similar to the symmetric matrix D L D , it has eigenvectors −1/2 w(j) v (j) = D(ε,N ) for j = 0, . . . , N − 1, where w(0) , . . . , w(N −1) forms an orthonormal basis of eigenvectors for −1/2 (ε,N ) (ε,N ) −1/2 D(ε,N ) L D . For each j = 0, . . . , N − 1, let λj denote the eigenvalue cor(j) responding to v . Since −1 −1 D(ε,N ) L(ε,N ) = I − D(ε,N ) W (ε,N ) , without loss of generality, we can order the eigenvalues 0 = λ0 ≤ λ1 ≤ · · · ≤ λN −1 (ε)
As in [12], but in contrast to [10, 11], we do not use a δ-neighborhood rule (Wij = 0 if D(xi , xj ) ≥ δ) or a kNN (ε)
rule (Wij = 0 if xj is not one of the k nearest neighbors of xi ).
6
and take v (0) to be a multiple of the all ones vector. For a fixed l ≤ N − 1, Laplacian eigenmaps maps each data point xi into Rl using the eigenvectors v (1) , . . . , v (l) : (1) vi . (8) xi 7→ .. . (l)
vi
Even when we are only given h as in Assumption 2.1 and D as in Assumption 2.2 (not necessarily the squared Euclidean distance), we can still apply the Laplacian eigenmap algorithm, in the sense −1 (ε,N ) that we can construct the matrix D(ε,N ) L and compute its eigenvectors.
2.2
Sinkhorn divergences
We now recall some facts about entropy-regularized optimal transport, which can be found in [26, 60, 63]. Let X ⊂ Rd be compact. For µ, ν ∈ P(X), define for each β > 0 Z := OTβ (µ, ν) inf ∥x − y∥2 dπ(x, y) + βKL(π|µ ⊗ ν), (9) π∈Π(µ,ν)
where Π(µ, ν) denotes the set of couplings between µ and ν and KL denotes the Kullback–Leibler divergence. (For this paper, we will only consider the squared Euclidean distance as the cost function.) Then for each µ, ν ∈ P(X), as β → 0, OTβ (µ, ν) converges to the square of the Wasserstein-2 distance 1/2 Z W2 (µ, ν) := inf ∥x − y∥2 dπ(x, y) π∈Π(µ,ν)
between µ and ν [60, Section 5]. However, for β > 0 and µ ∈ P(X), OTβ (µ, µ) may be non-zero, so [39] proposed to instead use the Sinkhorn divergence Sβ (µ, ν) := OTβ (µ, ν) −
OTβ (µ, µ) OTβ (ν, ν) − 2 2
(10)
as a loss function for machine learning. In [34], it is shown that for each β > 0, Sβ is non-negative and Sβ (µ, ν) = 0 if and only if µ = ν. From now on, we assume that β > 0 is fixed. The dual for (9) is given by Z Z Z f (x) + g(y) − ∥x − y∥2 OTβ (µ, ν) = sup f dµ + gdν − β exp − 1 d(µ ⊗ ν)(x, y) β f,g∈C(X) (11) for µ, ν ∈ P(X), and there exist fµ,ν , gµ,ν ∈ C(X), also known as Schrödinger potentials, which attain the maximum in (11). In addition, we have that πµ,ν given by fµ,ν (x) + gµ,ν (y) − ∥x − y∥2 d(µ ⊗ ν)(x, y) (12) dπµ,ν (x, y) = exp β is optimal for the primal (9). Since πµ,ν is a coupling between µ and ν, fµ,ν , gµ,ν satisfy Z fµ,ν (x) + gµ,ν (y) − ∥x − y∥2 exp dν(y) = 1 for µ-a.e. x β 7
(13)
Z exp
fµ,ν (x) + gµ,ν (y) − ∥x − y∥2 β
dµ(x) = 1 for ν-a.e. y
(14)
which can also be understood as the Euler-Lagrange equations or first-order optimality conditions for (11) [60, Remark 3.4]. Observe that maximizers of (11) are not unique, since if fµ,ν , gµ,ν is a maximizer, so is fµ,ν + λ, gµ,ν − λ for any λ ∈ R. Nonetheless, the Schrödinger potentials fµ,ν , gµ,ν are unique up to an additive constant upon requiring that (13) and (14) hold for all x, y ∈ X (rather than µ-a.e. x and ν-a.e. y), i.e., Z gµ,ν (y) − ∥x − y∥2 fµ,ν (x) = −β log exp dν(y) (15) β Z fµ,ν (x) − ∥x − y∥2 dµ(x) (16) gµ,ν (y) = −β log exp β for all x, y ∈ X.
2.3
Related work
The Laplacian eigenmap [10, 11, 12] and diffusion map [27, 28, 59] algorithms can be viewed as discrete analogues to embeddings of (M, g) using the eigenfunctions of the Laplace-Beltrami operator [9] and heat kernel [16], respectively. As such, under the setting where (M, g) is a compact Riemannian submanifold of Rd and D is the squared Euclidean distance, several notions of the convergence of graph Laplacians to the Laplace-Beltrami operator (or a weighted version) on (M, g) have been studied; see, for example, [12, 13, 14, 19, 20, 25, 27, 32, 37, 44, 71, 80, 83] and references therein. However, in several manifold learning algorithms (e.g., Isomap [79] in addition to the Laplacian eigenmap and diffusion map algorithms), one can replace the Euclidean distance with another measure of dissimilarity between data points. Indeed, depending on the application, other measures of dissimilarity may be more appropriate than the Euclidean distance; for examples, see [22, 46, 53, 57, 79, 87]. The particular application which motivates our Assumption 4.1 is that of understanding the conformation space of a molecule [30, 87]. When comparing two conformations of a molecule, an optimal transport-related distance is better for measuring displacements in space, especially when compared to the L2 distance between the densities. Therefore, it is not surprising that Zelesko et al. found using a wavelet approximation of the Wasserstein-1 metric [68] in the diffusion maps algorithm beneficial when analyzing the conformation space of a molecule [87]. Motivated by [87], [46] studied the pointwise convergence of the graph Laplacian when the typical Euclidean norm is replaced by an arbitrary norm on Rd , using tools from convex geometry. In [81], the spectral convergence of graph Laplacians constructed with Fermat distances is studied, though we note that Fermat distances depend on the data and thus our analysis will not apply. In [61], Oliver et al. study the consistency of Laplace learning on a submanifold of Wasserstein-2 space (see [61] and [42] for the precise assumptions) via Γ-convergence of the discrete Dirichlet energies. Though we focus on graph Laplacian-based methods, several authors have also studied adapting other dimensionality reduction and/or manifold learning methods to a space of probability measures; see, e.g., [4, 18, 41, 42, 66, 84]. The appeal of the Wasserstein-2 space is that it has a formal Riemannian structure, introduced by Otto in [62]. For example, the Benamou-Brenier dynamical formulation of the Wasserstein-2
8
distance [15] is similar to how length-minimizing geodesics can be characterized as energy minimizers in Riemannian geometry. We refer to [3, 26, 82] for more details on this Riemannian structure. However, the Wasserstein-2 space is not a true finite-dimensional smooth Riemannian manifold, so much care needs to be taken even when constructing a "Riemannian submanifold" of the Wasserstein-2 space [42] and defining the Laplace-Beltrami operator on such a submanifold [61]. Additional regularity constraints can bring us a bit closer to standard Riemannian geometry [52], but the space that [52] works with is still infinite-dimensional. In this work, we instead adapt our approach from [86], so that we only need to work with a finite-dimensional smooth Riemannian manifold. In [86], we observed that when a smooth compact orientable manifold M is embedded into a metric space (X, d), a sufficient condition for the pointwise convergence of graph Laplacians constructed with d2 is for d2 to be smooth on a neighborhood of the diagonal ΞM of M × M and for 21 Hessp d2 (p, ·) to be nondegenerate for all p ∈ M. The main idea is to equip M with the Riemannian metric gp :=
1 Hessp d2 (p, ·) 2
(17)
and establish that there exists K, ε0 > 0 such that d(p, q)2 ≤ dg (p, q)2 ≤ d(p, q)2 + Kdg (p, q)4 for all p, q ∈ M such that dg (p, q) < ε0 . From here, the pointwise convergence of the graph Laplacians follows mutatis mutandis from the Euclidean setting [12]; see Section 3 and Appendix A.1 for more details. More generally, we can consider gp :=
1 Hessp D(p, ·) 2
(18)
as a Riemannian metric so long as D : M × M → R≥0 is smooth on a neighborhood of ΞM , D(p, p) = 0 for all p ∈ M and (18) is nondegenerate for all p ∈ M. This idea dates back to Rao [64] and is central to information geometry [2, 8], which uses the Fisher information as a Riemannian metric. Namely, when D is given by, e.g., the Kullback–Leibler divergence or squared Hellinger distance, under some assumptions, one can show that (18) is a scalar multiple of the Fisher information [26, 47]. More concretely, in statistics, one often considers a family of probability measures {µx : x ∈ Θ} parametrized by a subset Θ ⊂ Rm ; under appropriate assumptions, we can define D(x, y) := KL(µx ||µy ), upon which (18) is half of the Fisher information at x [47]. (Note, however, that in Assumption 2.2, we also require D to be symmetric.) This can then be used in the natural gradient algorithm [1], a gradient descent algorithm where the gradient is taken with respect to g. Combining information geometry and the Riemannian structure of Wasserstein-2 space, more recent work has studied when the Kullback-Leibler divergence is replaced with (a variant of) the squared Wasserstein-2 distance [23, 50, 51]. (17) and/or (18) have also been used to study Isomap on a family of measures parametrized by an open set of Rm [4] and to approximate geodesics on a Riemannian manifold [30, 65]. Given the prevalence of the Sinkhorn divergence in machine learning and related areas, a natural direction then is to study (18) on a parametrized family of probability measures when D is given by the Sinkhorn divergence. This is done in [67], which is more focused on the natural gradient algorithm, and [49], which explores the properties of the Riemannian structure given by (18) in more detail. For this paper, we will follow [49]. We also mention the follow-up work [43], which 9
uses [49] to understand a variant of the JKO scheme with the Sinkhorn divergence. To the best of our knowledge, there has not been prior work studying graph Laplacian-based manifold learning algorithms when the Riemannian structure is given by the Hessian of the Sinkhorn divergence.
3
Approximations of the geodesic distance and convergence of the graph Laplacian
It is well-known that for a smooth, connected, compact Riemannian submanifold (M, g) of Rd , there exists K > 0 such that ∥p − q∥2 ≤ dg (p, q)2 ≤ ∥p − q∥2 + Kdg (p, q)4
(19)
for all p, q ∈ M, where dg denotes the geodesic distance on (M, g); see, e.g., [14, Lemma 4.3] or [17, Lemma 3]. Indeed, (19) is often helpful when analyzing manifold learning algorithms for data lying on M ⊂ Rd . Now consider the more general case, where we are only given a manifold M satisfying the assumptions in Table 1 and a function D satisfying Assumption 2.2. Endow M with the Riemannian metric g as defined in (3). For the rest of this section and Appendix A, let exp and dg denote the associated exponential map and geodesic distance, respectively. Fix any p ∈ M and v ∈ Tp M such that gp (v, v) = 1. A priori, from Taylor expanding around p, we may expect D to approximate the geodesic distance dg up to third order, i.e., for t ∈ R with |t| small enough, D(p, expp (tv)) = t2 + O(|t|3 ) = dg (p, expp (tv))2 + O(|t|3 )
(20)
as t → 0. If D = d2 for a metric d on M, then one can obtain that there exists K > 0 such that d(p, q)2 ≤ dg (p, q)2 ≤ d(p, q)2 + Kdg (p, q)4
(21)
for all p, q ∈ M using the triangle inequality and compactness of M [86]. However, a divergence is not necessarily the square of a distance, and in particular, the Sinkhorn divergence is not the square of a distance [49]. Instead, we will use the symmetry of D in Assumption 2.2 to improve upon the third-order error in (20): Proposition 3.1. With D satisfying Assumption 2.2, there exists K > 0 such that D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4 for all p, q ∈ M. The main idea is to apply the following lemma to D along the unit-speed geodesics of (M, g): Lemma 3.1. Suppose f ∈ C 4 ([a, b] × [a, b]) for some a < b and that f is symmetric, i.e., f (x, y) = f (y, x) for all x, y ∈ [a, b], that f ≥ 0 and f (x, x) = 0 for all x ∈ [a, b]. Let fx := f (x, ·) for each x ∈ [a, b]. Then ∂3f ∂3f (x, x) = −3 (x, x) (22) ∂y 3 ∂x∂y 2 10
for all x ∈ (a, b). In particular, if additionally we have fx′′ (x) = 2 for all x ∈ (a, b), then fx′′′ (x) = 0 for all x ∈ (a, b) and so ! 1 (4) 2 (23) f (x, y) − (x − y) ≤ sup f (z) (x − y)4 24 z∈(a,b) x for all x, y ∈ (a, b). Lemma 3.1 follows from Taylor expanding each fx and using the symmetry condition; see Appendix A. We also leave the proof of Proposition 3.1 for Appendix A, as it is similar to [86, Lemma 3.5]. As we cannot use the triangle inequality, in contrast to (21), we may not have D(p, q) ≤ dg (p, q)2 for all p, q ∈ M. Nonetheless, Proposition 3.1 is enough for us to prove, for example, the pointwise convergence of the graph Laplacian, following similar arguments as the Euclidean setting [14, 27]. Theorem 3.1 (c.f. [14, 27]). 2.1 and 2.2. Fix any f ∈ C 3 (M), x ∈ M and Assume Assumptions 1
α > 0, and suppose εN = Ω N − m+2+α
and εN → 0 as N → ∞. Define (as in [27]) Z m2 := v12 h(∥v∥2 )dv.
(24)
Rm
If x1 , . . . , xN are drawn i.i.d. from a probability distribution on M with a density P ∈ C 3 (M) with respect to the Riemannian volume form dVg (where g is as given in Assumption 2.2), then L(εN ,N ) f (x) m/2+1 N εN
→
m2 P (x)∆g f (x) − 2gx (gradg f (x), gradg P (x)) 2
almost surely as N → ∞. In particular, if x1 , . . . , xN are drawn i.i.d. from the uniform distribution on (M, g), then L(εN ,N ) f (x) m2 ∆g f (x) → m/2+1 2 volg (M) Nε N
almost surely as N → ∞. As it is very similar to the arguments in [14, 27], we leave the proof of Theorem 3.1 to Appendix A.1. We expect other results, such as the spectral convergence of graph Laplacians (see [20, 25] and references therein for the Euclidean setting), to hold as well, but we have not checked this for this work.
4
Application to Sinkhorn divergences
In this section, we discuss two models where a family of measures is parametrized by a manifold M and the Sinkhorn divergence is smooth around the diagonal of M × M. We then discuss when the non-degeneracy condition (5. in Assumption 2.2) holds. Throughout this section, β > 0 is fixed.
4.1
Regularity
We start with showing that the smoothness around the diagonal assumption (4. in Assumption 2.2) holds in two models of interest. 11
4.1.1
Particle system parametrized by a manifold
We first consider a particle system parametrized by a manifold: Assumption 4.1 (Particle system parametrized by a manifold). Let M be a smooth, connected, compact orientable m-dimensional manifold with m ≥ 1. Assume that we have a smooth embedding ι = (ι1 , · · · , ιn ) : M → P (Rd )n such that ιj (x) ̸= ιk (x) for all x ∈ M and j ̸= k. Fix some n p1 , . . . , pn > 0 such that j=1 pj = 1, and define F : M → P(Rd ) by F : M → P(Rd ) n X x 7→ pj διj (x) . j=1
Moreover, assume that F (x) ̸= F (y) for all x ̸= y. We claim that D(x, y) := Sβ (F (x), F (y)) is smooth on M × M under Assumption 4.1. This follows immediately from the following lemma. Lemma 4.1. Under Assumption 4.1, Gβ (x, y) := OTβ (F (x), F (y)) is smooth on M × M. The main idea is to use the implicit function theorem with the first-order optimality conditions (15) and (16). The proof is similar to the proof of [54, Theorem 2], but in contrast to Lemma 4.1, [54] assumes that the locations of the Dirac delta functions are fixed, whereas the coefficients pj are allowed to vary. The implicit function theorem approach has also been taken in [21, 49]. However, in these works, the measures are not necessarily assumed to be discrete with finite support, and so the proofs become more technical. We will discuss this more in Section 4.1.2, but here we first consider the particle system assumption (Assumption 4.1) to illustrate the main ideas. Proof of Lemma 4.1. For convenience, define Fβ : Rn × Rn × M × M → R as n X X fj + gk − ∥ιj (x) − ιk (y)∥2 − 1 pj pk , (fj + gj )pj − β Fβ (f, g, x, y) := exp β j=1
(25)
j,k
so that by (11), we have Gβ (x, y) = sup Fβ (f, g, x, y). f,g∈Rn
(26)
As discussed in Section 2.2, this problem admits maximizers (the Schrödinger potentials), though we also have Fβ (f + λ1, g − λ1, x, y) = Fβ (f, g, x, y) for all λ ∈ R and x, y ∈ M, where 1 denotes the n-dimensional vector of all ones. However, if we consider Rn × Rn / ∼, where (f, g) ∼ (fe, ge) if fe = f + λ1, ge = g − λ1 for some λ ∈ R, then for each x, y ∈ M, Fβ (·, ·, x, y) is strictly concave on Rn × Rn / ∼. 12
For concreteness, let us identify Rn × Rn / ∼ with the subspace P = {(f, g) ∈ Rn × Rn |f1 = 0}, in which case the strict concavity of Fβ (·, ·, x, y) on P implies that the maximizer f ∗ (x, y), g ∗ (x, y) of supf,g∈P Fβ (f, g, x, y) is unique in P for each x, y ∈ M. Define T (β) (f, g, x, y) := ∇f2 ,...,fn ,g1 ,...,gn Fβ (f, g, x, y),
(27)
where ∇f2 ,...,fn ,g1 ,...,gn denotes the gradient in the f2 , . . . , fn , g1 , . . . , gn directions. For each x, y ∈ M, using the first-order optimality conditions for (26) and the strict concavity of Fβ (·, ·, x, y), f ∗ (x, y), g ∗ (x, y) ∈ P are characterized by T (β) (f ∗ (x, y), g ∗ (x, y), x, y) = 0.
(28)
Observe that we can write the first-order optimality conditions for (26) more explicitly as X fj + gk − ∥ιj (x) − ιk (y)∥2 pk = 0 pj − pj exp β k X fk + gj − ∥ιk (x) − ιj (y)∥2 pj − pj exp pk = 0, β k
for each j = 1, . . . , n, which, upon rearranging, are precisely the Schrödinger system (13) and (14). For each x, y ∈ M, when restricted to a bounded convex subset of P , Fβ (·, ·, x, y) is strongly concave with respect to the norm 1/2 n X |||(f, g)||| := ∥fj + gk ∥2 , j,k=1
hence with respect to the Euclidean norm as well. (Here ||| · ||| is only a norm on P — not on Rn × Rn — since we require f1 = 0 for all (f, g) ∈ P .) Therefore for each x, y ∈ M, the Jacobian of T (β) at f ∗ (x, y), g ∗ (x, y), x, y with respect to the f2 , . . . , fn , g1 , . . . , gn variables is invertible, so the implicit function theorem implies that f ∗ and g ∗ are smooth functions on a neighborhood of (x, y). Hence Gβ is smooth on a neighborhood of (x, y) for any x, y ∈ M. 4.1.2
Smooth deformations of a measure
In this subsection, we consider measures parametrized by an open set U ⊂ Rm instead of M, but our results can be extended to a manifold by passing to charts. Assumption 4.2 (Deformations of a measure). Fix any k ∈ Z>0 . Let X = Ω be the closure of a bounded open convex subset Ω ⊂ Rd , and let U ⊂ Rm be an open set. Let µ ∈ P(X). Suppose Ψ : U × X → X is a continuous function and that its first k partial derivatives with respect to the first variable ∂kΨ : U × X → Rd (29) ∂pi1 ∂pi2 · · · ∂pik exist and are continuous on U × X for all i1 , · · · , ik ∈ [m]. Define F : U → P(X) p 7→ Ψ(p, ·)# µ. We will sometimes also write µp for F (p) = Ψ(p, ·)# µ. Assume that µp ̸= µq for all p ̸= q. 13
Proposition 4.1. Under Assumption 4.2, Gβ (p, q) := OTβ (F (p), F (q)) is C k on U × U . Remark 4.1. In Assumption 2.2, we assumed that D is smooth (i.e., infinitely differentiable) on a neighborhood of the diagonal ΞM ⊂ M × M, but with a careful reading one can track the precise number of derivatives required for Proposition 3.1 and Theorem 3.1. Alternatively, one can assume that Assumption 4.2 holds for all k ∈ Z>0 , upon which Proposition 4.1 implies that D(p, q) := Sβ (F (p), F (q)) is infinitely differentiable on U × U . We will prove Proposition 4.1 at the end of this subsection, after first establishing a few lemmata. Define (as in [21, Definition 2.1]) the map T (β) : C(X)2 × P(X)2 → C(X)2 by Z f (x) + g(y) − ∥x − y∥2 (β) dν(y) (30) exp T1 (f, g, µ, ν)(x) := log β Z f (y) + g(x) − ∥x − y∥2 (β) T2 (f, g, µ, ν)(x) := log exp dµ(y) (31) β for f, g ∈ C(X) and µ, ν ∈ P(X). Lemma 4.2. Under Assumption 4.2, the map (f, g, p, q) 7→ T (β) (f, g, F (p), F (q)) is a C k map from C k (X) × C k (X) × U × U to C k (X) × C k (X). Remark 4.2. Here we use C k (X) to refer to the space of functions in C k (Ω) whose derivatives of order ≤ k have continuous extensions to X = Ω, equipped with the norm ∥f ∥C k (X) :=
sup
|∂ α f (x)| ,
|α|≤k,x∈Ω
which makes it a Banach space. If one only assumes that X is compact and convex, then [21] (which we will use in the proof of Lemma 4.3) instead defines C k (X) as the space of functions on X with a C k extension on Rd with the norm |||f |||C k (X) := inf
sup
∂ α fe(x) ,
fe |α|≤k,x∈Rd
where the infimum is taken over all C k extensions of f to Rd . In our setting, where X = Ω with Ω ⊂ Rd open, bounded and convex, these two definitions are equivalent (up to equivalence of norms) by Whitney’s extension theorem [85] (see also [45, Theorem 2.3.6 and proof of Theorem 2.3.10]). In particular, since any f ∈ C k (X) has a C k extension fe on an open neighborhood of X, under Assumption 4.2 we can apply the chain rule to p 7→ f (Ψ(p, x)) for each x ∈ X and f ∈ C k (X). We leave the proof of Lemma 4.2 for Appendix B, though we note that this lemma is similar to [21, Lemma 3.4] and [49, Lemma 3.13]. [21, Lemma 3.4] is primarily concerned with families of measures arising from interpolations, whereas [49, Lemma 3.13] assumes different conditions, which 14
are satisfied, for instance, by µt = Ψ(t, ·)# µ for t ∈ (−δ, δ), for some δ > 0 and Ψ with only one continuous derivative with respect to the first variable (what they call horizontal perturbations) [49, Remark 3.2]. Therefore [49, Lemma 3.13] only proves the continuous differentiability of T (β) . We assume more regularity of Ψ with respect to the first variable (when k > 1) to obtain additional regularity of T (β) . Following the approach in Lemma 4.1, we now apply the implicit function theorem. Lemma 4.3. Assume Assumption 4.2. For each p, q ∈ U , let fp,q and gp,q denote the Schrödinger potentials between F (p) = Ψ(p, ·)# µ and F (q) = Ψ(q, ·)# µ. Then p, q 7→ fp,q , gp,q is a C k function e k (X) := C k (X) × C k (X)/ ∼, where from U × U to C (f, g) ∼ (fe, ge) ⇔ fe = f + λ, ge = g − λ for some λ ∈ R. e k (X) with the quotient norm, so that it is a Banach space. Here we equip C Proof. For each p, q ∈ U , the Schrödinger system (15) and (16) for F (p) = Ψ(p, ·)# µ and F (q) = Ψ(q, ·)# µ can be rewritten as T (β) (fp,q , gp,q , F (p), F (q)) = 0.
(32)
Using the definition of T (β) and Lemma 4.2, (f, g, p, q) 7→ T (β) (f, g, F (p), F (q)) is well-defined as e k (X) × U × U to C e k (X), and moreover it is C k . In addition, by [21, Lemma 3.2] a map from C (see also [40]), for any µ, ν ∈ P(X), Df,g T (β) (fµ,ν , gµ,ν , µ, ν) is a Banach space isomorphism from e k (X) to itself. Here Df,g denotes the derivative with respect to the C e k (X) variable. Although we C k k stated (32) as an equivalence of functions in C (X) × C (X), if for some f, g ∈ C k (X), µ, ν ∈ P(X) and λ ∈ R Z f (x) + g(y) − ∥x − y∥2 (β) T1 (f, g, µ, ν)(x) = log exp dν(y) = λ β Z f (y) + g(x) − ∥x − y∥2 (β) dµ(y) = −λ T2 (f, g, µ, ν)(x) = log exp β for all x ∈ X, then λ
e =
Z Z exp
f (x) + g(y) − ∥x − y∥2 β
dν(y)dµ(x) = e−λ ,
e k (X), it so λ must be 0. Therefore, even if we take (32) as an equivalence of functions in C characterizes the Schrödinger potentials for F (p) = Ψ(p, ·)# µ and F (q) = Ψ(q, ·)# µ. The lemma statement then follows from the implicit function theorem for Banach spaces [48, Theorem 5.9]. By Lemma 4.3 and the regularity of Ψ in Assumption 4.2, we obtain the following: Lemma 4.4. Following the notation from Lemma 4.3, under Assumption 4.2, p, q 7→ fp,q (Ψ(p, ·)), gp,q (Ψ(q, ·)) e e k (X), the space C(X) e := C(X) × C(X)/ ∼. As with C is a C k map from U × U to C(X) is equipped with the quotient norm, making it a Banach space. 15
As with Lemma 4.2, the proof involves checking that partial derivatives exist and are continuous with respect to appropriate norms; we leave the proof for Appendix B. Using Lemma 4.4, we can now prove Proposition 4.1. Proof of Proposition 4.1. Under Assumption 4.2, we have Z Z Z Gβ (p, q) = OTβ (F (p), F (q)) = fp,q dµp + gp,q dµq = fp,q (Ψ(p, x)) + gp,q (Ψ(q, x))dµ(x). Although fp,q , gp,q is technically an equivalence class (as an element of C(X) × C(X)/ ∼), the equation above holds for any representative. Since p, q 7→ fp,q (Ψ(p, ·)), gp,q (Ψ(q, ·)) is C k from U × U to C(X) × C(X)/ ∼ by Lemma 4.4, Gβ is C k on U × U .
4.2
Non-degeneracy
In this section, we discuss sufficient conditions for the non-degeneracy condition (5. in Assumption 2.2) to be satisfied. For this section, unless otherwise stated, X = Ω is the closure of a bounded open set Ω ⊂ Rd , µ ∈ P(X) and j ∈ Z≥0 . To check whether the non-degeneracy condition is satisfied, we use results on the Hessian of the Sinkhorn divergence from [49], which involve reproducing kernel Hilbert spaces (RKHS). We recall some facts from RKHS theory but refer to other references, e.g., [58, 78], for a more thorough exposition. In what follows, Hk denotes the RKHS with the symmetric, positive definite kernel k : X × X → R, where the input set X will be X or Rn , depending on context. For this paper, we will just work with RKHSs of real-valued functions. Recall that by the Moore-Aronzajn theorem [5], whenever k is symmetric and positive definite, there exists a RKHS Hk with kernel k and it is unique. A popular method in statistics and machine learning is to embed a probability measure µ into Hk via the map Z µ 7→ k(x, ·)dµ(x), known as the kernel mean embedding R [58, 75]. By the Riesz representation theorem, the right-hand side is indeed in Hk whenever f 7→ f dµ is a bounded linear functional on Hk . More generally, we follow [69] and define the kernel mean embedding for any element of Hk∗ : Definition 4.1 (Kernel mean embedding). For a symmetric, positive definite kernel k : X × X → R and its corresponding RKHS Hk , define the kernel mean embedding Hk : Hk∗ → Hk by Hk [ν](y) := ⟨ν, k(·, y)⟩ for each y ∈ X and ν ∈ Hk∗ . One can check that Hk [ν] ∈ Hk for each ν ∈ Hk∗ by observing that it is precisely the representation ψν ∈ Hk of ν ∈ Hk∗ given by the Riesz representation theorem: ψν (y) = ⟨ψν , k(·, y)⟩Hk = ⟨ν, k(·, y)⟩ = Hk [ν](y). As in [49], we are interested in RKHSs with the self-transport kernel: Definition 4.2 (Self-transport kernel). For µ ∈ P(X), the self-transport kernel kµ : X × X → R is given by fµ,µ (x) + gµ,µ (y) − ∥x − y∥2 kµ (x, y) := exp , β 16
where fµ,µ , gµ,µ are Schrödinger potentials for (µ, µ) (with regularization parameter β > 0). Following [49], for the rest of this section, we use the convention fµ,µ = gµ,µ , so that we can write fµ,µ (x) + fµ,µ (y) − ∥x − y∥2 kµ (x, y) = exp . β As the Gaussian kernel is a positive definite kernel, so is kµ for each µ ∈ P(X). As in [49], we will use Hµ (instead of Hkµ ) to denote the RKHS with kernel kµ , and we will use Hµ to denote the kernel mean embedding associated to Hµ . Since the Schrödinger potential fµ,µ is smooth, we have kµ ∈ C ∞ (X × X), which implies that Hµ injects continuously into C j (X) (see [69] or [49, Appendix B]) and Hµ is also well-defined as an operator from C j (X)∗ to C j (X). [49] additionally define Kµ : C(X) → C j (X) by Z Kµ [φ](y) := φ(x)kµ (x, y)dµ(x) for each φ ∈ C(X) [49, Definition 3.3]. We can now state the results from [49] that we will use: as long as a curve of measures {µt }t∈(−δ,δ) is continuously differentiable in C j (X)∗ with the weak-* topology (what they call a C j -perturbation), lim
t→0
β β Sβ (µ0 , µt ) 2 = ⟨µ̇, (id − Kµ2 )−1 Hµ [µ̇]⟩ ≥ ∥Hµ [µ̇]∥Hµ , t2 2 2
(33)
where µ = µ0 and µ̇ ∈ C j (X)∗ acts on φ ∈ C j (X) via ⟨µ̇, φ⟩ =
d ⟨µt , φ⟩ t=0 . dt
(34)
Observe that for any φ ∈ C j (X) and constant c ∈ R, ⟨µ̇, φ⟩ = ⟨µ̇, φ + c⟩, and so we can view µ̇ as acting on C j (X)/R, the quotient space obtained by quotienting out the constant functions in C j (X). In particular, [49] show that id − Kµ2 : C j (X)/R → C j (X)/R is well-defined and has a bounded inverse. Hence ⟨µ̇, (id − Kµ2 )−1 Hµ [µ̇]⟩ in (33) is well-defined. It is now much more straightforward to check whether the non-degeneracy condition (5. in Assumption 2.2) holds under Assumption 4.1 and Assumption 4.2. Proposition 4.2. Under Assumption 4.1, the non-degeneracy condition (5. in Assumption 2.2) holds. Under Assumption 4.2, if for all p ∈ U and i ∈ [m], there exists φ ∈ C 1 (X) such that Z ∂Ψ (p, x) · ∇φ(Ψ(p, x))dµ(x) ̸= 0, (35) ∂pi the non-degeneracy condition also holds. Proof. Let us work under Assumption 4.1 first; the arguments for Assumption 4.2 are analogous. Let γ : (−δ, δ) → M be a smooth curve with γ(0) = x ∈ M and γ ′ (0) = v ∈ Tx M with v ̸= 0. Set µt = F (γ(t)) =
n X j=1
17
pj διj (γ(t))
for t ∈ (−δ, δ), and let X = BR (0) for some R > 0 large enough that Im(ιj ◦ γ) ⊂ X for all j ∈ [n]. We first check that t 7→ µt is continuously differentiable in C 1 (X)∗ with the weak-* topology, i.e, that µt+s − µt (36) µ̇t := lim s→0 s exists for each t ∈ (−δ, δ) (where the limit is taken in the weak-* topology) and that t 7→ µ̇t is continuous. This is similar to the argument in [49, Remark 3.2]: for any t ∈ (−δ, δ) and φ ∈ C 1 (X), lim
s→0
n X µt+s − µt φ(ιj (γ(t + s))) − φ(ιj (γ(t))) , φ = lim pj s→0 s s j=1 =
n X
pj ∇φ(ιj (γ(t))) · (ιj ◦ γ)′ (t),
j=1
so µ̇t ∈ C 1 (X)∗ exists and is given by ⟨µ̇t , φ⟩ =
n X
pj ∇φ(ιj (γ(t))) · (ιj ◦ γ)′ (t)
j=1
for φ ∈ C 1 (X). For any φ ∈ C 1 (X), t 7→ ⟨µ̇t , φ⟩ is continuous, so t 7→ µ̇t is continuous from (−δ, δ) to C 1 (X)∗ with the weak-* topology. Recall that in Assumption 4.1 we assumed that ι is a smooth embedding, so (ι ◦ γ)′ (t) is nonzero whenever γ ′ (t) is non-zero. In particular, since we took γ : (−δ, δ) → M to be such that γ ′ (0) = v ̸= 0, there exists j ∈ [n] such that (ιj ◦ γ)′ (0) is non-zero, so µ̇0 is non-trivial in C 1 (X)∗ . We claim that µ̇0 is in fact non-trivial in Hµ∗ 0 , so that using (33) we can conclude that the nondegeneracy condition holds. Indeed, it suffices to observe that Hµ is dense in C 1 (X) for any µ ∈ P(X), as this would imply that any non-trivial element in C 1 (X)∗ is also non-trivial in Hµ∗ . The argument is similar to an argument in [49, Appendix B], where they show that Hµ is dense in C(X) (i.e., kµ is universal). Let kβ (x, y) := exp −∥x − y∥2 /β for each x, y ∈ X. Then for any φ ∈ C 1 (X), c1 , c2 , . . . , cl ∈ R and y1 , . . . , yl ∈ X, we have l X
ci kµ (yi , ·) −φ exp (fµ,µ (yi )/β) i=1
≤ 2 ∥exp (fµ,µ /β)∥C 1 (X)
l X i=1
C 1 (X)
ci kβ (yi , ·) −
φ exp (fµ,µ /β)
, C 1 (X)
and the right hand side can be made arbitrarily small, given that the RKHS for the Gaussian kernel kβ is dense in C 1 (X) [69] (see Remark 4.3 below). The proof of the statement for Assumption 4.2 is similar. Let {ei }i∈[m] denote the standard coordinate basis for Rm . With U, Ψ, X as in Assumption 4.2, fix any p ∈ U, i ∈ [m] and δ > 0 such that Bδ (p) ⊂ U , and let µt := Ψ(p + tei , ·)# µ for t ∈ (−δ, δ). Then, as in [49, Remark 3.2], by the chain rule and dominated convergence theorem, µ̇t = lims→0 (µt+s − µt )/s ∈ C 1 (X)∗ exists and is given by Z ∂Ψ ⟨µ̇t , φ⟩ = (p + tei , x) · ∇φ(Ψ(p + tei , x))dµ(x) ∂pi 18
for each φ ∈ C 1 (X), and so t 7→ ⟨µ̇t , φ⟩ is continuous for any φ ∈ C 1 (X). Therefore t 7→ µt is continuously differentiable in C 1 (X)∗ with the weak-* topology. The condition that (35) holds for some φ ∈ C 1 (X) ensures that µ̇0 is non-trivial in C 1 (X)∗ . The density of Hµ0 in C 1 (X) implies that µ̇0 is non-trivial in Hµ∗ 0 and we can conclude. Remark 4.3 (On the universality of the Gaussian kernel). It is well-known that the RKHS Hkβ ,X with the Gaussian kernel kβ (x, y) := exp −∥x − y∥2 /β is dense in C(X) whenever X ⊂ Rd is compact [77]. We include X in the subscript of Hkβ ,X to emphasize that each element of Hkβ ,X is a function over X. When a RKHS Hk is dense in C(X), its kernel k is called universal [56, 77] or, more specifically, c-universal [69, 76]. The terminology c-universal is used to emphasize that the RKHS Hk is dense in C(X). One can also study whether a RKHS is dense in other function spaces. In particular, [69] show that the RKHS Hkβ ,Rd (i.e., each element of Hkβ ,Rd is a function on Rd ) with the Gaussian kernel kβ is dense in C01 (Rd ). Using the terminology in [69], the Gaussian kernel is c10 -universal when the input set is Rd . Here the space C01 (Rd ) is the space of C 1 functions that vanish at infinity along with its first derivatives, equipped with the topology generated by the semi-norms ∥f ∥α := supx∈Rd |∂ α f (x)|, |α| ≤ 1. We refer to [70] for more details. What we need for the proof of Proposition 4.2 is for Hkβ ,X to be dense in C 1 (X) whenever X = Ω for a bounded open convex set Ω ⊂ Rd . This follows from the density of Hkβ ,Rd in C01 (Rd ) by using [69, Theorem 6], which states that if a RKHS Hk with kernel k injects continuously into a locally convex topological vector space of functions F, Hk is dense in F if and only if k is strictly positive definite over F ∗ , i.e., 2 ∥Hk [ν]∥Hk = ⟨ν, Hk [ν]⟩ > 0 for all non-trivial ν ∈ F ∗ . (We used the forward direction in the proof of Proposition 4.2; the reverse direction follows from Hahn-Banach.) Thus the density of Hkβ ,Rd in C01 (Rd ) implies that ⟨ν, Hkβ [ν]⟩ > 0
(37)
for all non-trivial ν ∈ C01 (Rd )∗ , where we recall that Hkβ [ν](y) = ⟨ν, kβ (·, y)⟩ for y ∈ Rd . Since any f ∈ C 1 (X) can be extended to some fe ∈ C01 (Rd ) (see Remark 4.2), (37) is true for all non-trivial ν ∈ C 1 (X)∗ , and so Hkβ ,X is dense in C 1 (X).
5
Examples
In this section, we discuss two examples where D is given by the Sinkhorn divergence. Unless otherwise specified, β > 0 is fixed. For the figures in this section, we set β = 1. The code for all of the figures in this paper can be found at https://github.com/lzx23/graph_laplacians_w_sdiv.
5.1
Rotating particle
Fix some R > 0. Consider F : S1 → P(R2 ) given by F (θ) :=
1 1 δ(0,0) + δ(R cos θ,R sin θ) . 2 2 19
(38)
For convenience, let µθ := F (θ) for θ ∈ [0, 2π). For each θ ∈ [−π, π), we have a closed-form formula for OTβ (µ0 , µθ ). Claim 5.1. For each θ ∈ [−π, π), the optimal coupling for OTβ (µ0 , µθ ) =
min
2 2 X θ +β Πij log (4Πij ) , (Π12 + Π21 )R + Π22 2R sin 2 i,j=1 2
Π∈R2×2 , 1 Π=[1/2,1/2], Π1=[1/2,1/2]⊤ ⊤
is given by Π∗θ =
1
p(θ) := We then have
p(θ) , 1 2 − p(θ)
(40)
1 . 2 1 + eR2 cos(θ)/β
(41)
2 − p(θ)
where p is given by
(39)
p(θ)
OTβ (µ0 , µθ ) = β log(2 − 4p(θ)) + R2 (1 − cos(θ)),
(42)
from which one can compute Sβ (µ0 , µθ ) = OTβ (µ0 , µθ ) −
1 (OTβ (µ0 , µ0 ) + OTβ (µθ , µθ )) = OTβ (µ0 , µθ ) − OTβ (µ0 , µ0 ) 2 d d 1 (β) gθ , = R2 (1 − 2p(0)) dθ dθ 2
for each θ ∈ [−π, π). We leave the proof of Claim 5.1 for Appendix C. For the plots in this subsection, we will use the closed-form formulae for Sβ and g (β) , and, unless otherwise specified, R = 2 and β = 1. For the rest of the subsection, we will write g for g (1) . Remark 5.1. In practice, Sinkhorn’s algorithm (or variants thereof ) are used to approximately solve the entropy-regularized optimal transport problem [29]. We did not include the error from Sinkhorn’s algorithm in our analysis in Section 4 (i.e., we have assumed that we can compute the Sinkhorn divergences exactly) and leave this for future work. Nonetheless, for this example and the parameters that we use, we find empirically that the default implementation of Sinkhorn’s algorithm in the Python Optimal Transport toolbox [35, 36] provides a good enough approximation of the Sinkhorn divergence that the resulting approximation error for the discrete Laplacian is very small; see Appendix C for more details. Fourth-order error. We first discuss Proposition 3.1. For this example, a more straightforward way of seeing that the fourth-order error in Proposition 3.1 holds is to observe that θ 7→ Sβ (µ0 , µθ ) is smooth and even, so that there is no third-order term in its Taylor expansion around 0. We can also observe that Proposition 3.1 holds numerically in Fig. 1. Convergence of the discrete Laplacian. Next, in Fig. 2 we test the pointwise convergence of the discrete Laplacian constructed with D(θ, φ) := S1 (µθ , µφ ) for θ, φ ∈ S1 . We draw N samples θ1 , . . . , θN i.i.d. uniformly from [0, 2π), take h(x) = e−x/2 as the kernel function for constructing 20
Figure 1: Visualization of the approximation error in Proposition 3.1 for the rotating particle example in Section 5.1.
Figure 2: Example of pointwise convergence of the discrete Laplacian for the example in Section 5.1. 2volg (M) L(εN ,N ) f (θ) We plot θ 7→ for N = 250, 500, 1000, 2000, 4000 and compare with ∆g f . m/2+1 m2 N εN
21
Figure 3: Embedding of samples into R2 using the first two non-constant eigenvectors v (1) , v (2) of (D(εN ,N ) )−1 L(εN ,N ) for the example in Section 5.1. Each point is colored according to the value of θi . the weight matrix, set εN = 4N −1/3.01 and test using the function f (θ) = sin(θ)+ 12 cos(2θ). For this √ choice of h, m2 as defined in (24) can be computed to be m2 = (2π)m/2 = 2π. For convenience, we plot θ 7→
2volg (M) L(εN ,N ) f (θ) 1 m/2+1 , so that by Theorem 3.1, for each θ ∈ S , the limit as N m2 Nε
(almost surely) should be
→∞
N
∆g f (θ) = −
1 g0
d
d dθ , dθ
d2 f (θ). dθ2
(43)
Forqour numerical experiments, we use the closed-form formula of g to compute (43) and volg (M) = d d 2π g0 dθ , dθ . (1)
(2)
In Fig. 3, à la the Laplacian eigenmap algorithm, we map each sample θi to vi , vi , where v (1) , v (2) are the first two non-constant eigenvectors of (D(εN ,N ) )−1 L(εN ,N ) , as introduced in Section 2.1. We color each point according to the value of θi . As expected, we can see the underlying geometry of (S1 , g) in Fig. 3.
5.2
Dilation of two particles
In this subsection, we discuss the following example from [49, Section 7.2]: for x > 0, set µx :=
1 1 δx + δ−x . 2 2
Then, as shown in [49, Lemma 7.5], the Sinkhorn divergence between µx and µy for x, y > 0 is given by the closed-form formula 1 Sβ (µx , µy ) = β − log (kβ (x, y) + kβ (x, −y)) + (log (1 + kβ (x, −x)) + log (1 + kβ (y, −y))) , 2 (44) where kβ (x, y) := exp −(x − y)2 /β is the Gaussian kernel. From (44), one can show (see [49, Eq. 7.3]) that the associated Riemannian metric on (0, ∞) is given by d d S (µ , µ ) k(x) 8x2 := lim β x 2y = 1 + gx , −2 (45) y→x (y − x) dx dx 1 + k(x) β(1 + k(x))
22
for x ∈ (0, ∞), where k(x) :=kβ (x, −x). For convenience, we will use Sβ (x, y) to denote Sβ (µx , µy ) d d and g(x) to denote gx dx , dx for x, y ∈ (0, ∞). We also observe that the squared geodesic distance on (0, ∞) with respect to the Riemannian metric g in (45) is given by dg (x, y)2 =
Z y
p g(z)dz
2
(46)
x
for x, y ∈ (0, ∞). Fourth-order error. We now discuss how Lemma 3.1 can be applied to this example. Fix any a ∈ (0, ∞). Consider the Taylor expansions of Sβ (a, ·) and d2g (a, ·) around a: Sβ (a, a + h) = g(a)h2 + d2g (a, a + h) = g(a)h2 + Lemma 3.1 implies that
1 ∂ 3 Sβ (a, a)h3 + O(h4 ) 6 ∂y 3
1 ∂ 3 [d2g ] (a, a)h3 + O(h4 ). 6 ∂y 3
Sβ (a, a + h) − d2g (a, a + h) = O(h4 )
(47)
∂ 3 [d2g ] ∂ 3 Sβ (a, a) = (a, a). ∂y 3 ∂y 3
(48)
as h → 0, or equivalently,
One way to see this is to take the approach of Proposition 3.1: consider some small enough T > 0 and a unit-speed, length-minimizing geodesic segment (with respect to g) γ : [−T, T ] → (0, ∞) such that γ(0) = a. Now apply (23) in Lemma 3.1 to f (s, t) := Sβ (γ(s), γ(t)) to show that there exists K > 0 such that Sβ (a, γ(t)) − t2 ≤ Kt4 for all t ∈ (−T, T ), or equivalently, Sβ (a, γ(t)) − d2g (a, γ(t)) ≤ Kd4g (a, γ(t)) for all t ∈ (−T, T ). For this example, upon fixing any a1 , a2 ∈ (0, ∞) such that a1 < a2 , we have ! inf z∈[a1 ,a2 ]
g(z) (x − y)2 ≤ d2g (x, y) ≤
sup
g(z) (x − y)2
z∈[a1 ,a2 ]
for all x, y ∈ [a1 , a2 ], so (47) holds. Alternatively, we can directly compute the derivatives of d2g with respect to the second variable: Z y p p ∂[d2g ] (a, y) = 2 g(z)dz g(y) ∂y a ! Z y ′ p ∂ 2 [d2g ] g (y) p (a, y) = 2 g(z)dz + g(y) ∂y 2 2 g(y) a 23
for y ∈ (0, ∞), so that ∂ 3 [d2g ] (a, a) = 3g ′ (a). ∂y 3 ∂2S
Recall that by construction g(x) = 12 ∂y2β (x, x), so we have ∂ 3 [d2g ] 3 (a, a) = ∂y 3 2
∂ 3 Sβ ∂ 3 Sβ ∂ 3 Sβ (a, a) + (a, a) = (a, a), ∂x∂y 2 ∂y 3 ∂y 3
where the second equality follows from Lemma 3.1. We check (47) numerically for β = 1 in Fig. 4. In particular, we set δ = .005 and plot S1 (1, x) − dbg (1, x)2 S1 (1, x) − g(1)(x − 1)2 S1 (1, x) − dbg (1, x)2 S1 (1, x) − g(1)(x − 1)2 , , , (x − 1)3 (x − 1)3 (x − 1)4 (x − 1)4
at x = 1 + kδ for k ∈ Z \ {0} such that 1 + kδ ∈ [.9, 1.1]. Here dbg denotes a numerical approximation for the geodesic distance dg which we will describe below. In Fig. 4, both S1 and g(1) are computed using their closed-form formulae, though we obtain similar results if we compute S1 using Sinkhorn’s algorithm as implemented in the Python Optimal Transport library [35, 36]; see Appendix C. As x 7→ g(1)(x − 1)2 is the second-order Taylor expansion of x 7→ S1 (1, x) around 1, we expect |S1 (1,x)−g(1)(x−1)2 | S1 (1,x)−g(1)(x−1)2 to converge as x → 1, though may go to infinity as x → 1. (x−1)3 (x−1)4 This is indeed what we see in Fig. 4. In contrast, if (47) holds, we expect
S1 (1,x)−dg (1,x)2 (x−1)4
to stay
S (1,x)−db (1,x)2 bounded as x → 1. This agrees with the behavior of 1 (x−1)g4
in Fig. 4. b Last but not least, we describe how we compute dg : for each k ∈ N such that kδ ∈ [.9, 1.1), dg (kδ, (k + 1)δ) is approximated by applying Simpson’s rule to (46), and these estimates are then summed up to approximate the geodesic distance dg (1, 1 + kδ) for k ∈ Z such that kδ ∈ [−.1, .1]. Remark 5.2 (On the approximation error for the geodesic distance). Per the standard error bound for Simpson’s rule (see, e.g., [6]), for each k ∈ N such that kδ ∈ [.9, 1.1), our approximation dbg (kδ, (k + 1)δ) of dg (kδ, (k + 1)δ) satisfies 5
Cδ , dbg (kδ, (k + 1)δ) − dg (kδ, (k + 1)δ) ≤ 2880 where C is an upper bound on the absolute value of the fourth derivative of S1 (1, 1 + kδ) − dbg (1, 1 + kδ)2 S1 (1, 1 + kδ) − dg (1, 1 + kδ)2 C|k|δ 5 − ≤ (kδ)3 (kδ)3 2880|kδ|3
√
g on [.9, 1.1]. Hence
C|k|δ 5 2dg (1, 1 + kδ) + 2880
and similarly S1 (1, 1 + kδ) − dbg (1, 1 + kδ)2 S1 (1, 1 + kδ) − dg (1, 1 + kδ)2 C|k|δ 5 C|k|δ 5 − ≤ 2d (1, 1 + kδ) + g (kδ)4 (kδ)4 2880|kδ|4 2880
for k ∈ Z \ {0} such that kδ ∈ [−.1, .1].
24
Figure 4: Visualization of the approximation error in (47) for the dilation example in Section 5.2. Both S1 and g(1) are computed using their closed-form formulae, whereas the geodesic distance is approximated using Simpson’s rule. 25
Figure 5: Examples of pointwise convergence of the discrete Laplacian for the example in Sec2volg (M) L(εN ,N ) f (x) tion 5.2. As in Fig. 2, we plot x 7→ , but in the top plot, L(εN ,N ) f is constructed m/2+1 m2 N εN
using D(x, y) = S1 (x, y) for x, y ∈ M, whereas in the bottom plot, L(εN ,N ) f is constructed using D(x, y) = (x − y)2 for x, y ∈ M.
Figure 6: Embedding of samples into R2 using the first two non-constant eigenvectors v (1) , v (2) of (D(εN ,N ) )−1 L(εN ,N ) (constructed with D = S1 ) for the example in Section 5.2. Each point is colored according to the value of xi . 26
Convergence of the discrete Laplacian. We also provide an example of the pointwise convergence of the discrete Laplacian in Fig. 5. For the rest of the subsection, we take β = 1. Although in previous sections we assumed that M does not have a boundary, for a manifold with boundary, the pointwise convergence of discrete Laplacians can be shown at points away from the boundary, following a similar proof. We use rejection sampling to sample x1 , . . . , xN i.i.d. from the uniform distribution with respect to g (as in (45), with β = 1) on M = [.8, 1.2], which has density p g(x) (49) P (x) = volg (M) with respect to the Lebesgue measure on [.8, 1.2]. We set h(x) = e−x/2 and εN = .15N −1/3.01 and test using the function f (x) := x4 − 8x.√ As noted previously, with this choice of h, the constant m2 in (24) is given by m2 = (2π)m/2 = 2π. Similar to Fig. 2, we plot x 7→
2volg (M) L(εN ,N ) f (x) , m/2+1 m2 Nε N
so that when L(εN ,N ) f is constructed using D(x, y) := S1 (x, y), we expect convergence to 1 d 1 df 1 d2 f 1 dg df =− + 2 ∆g f := − √ √ 2 g dx g dx g dx 2g dx dx dg and approximate volg (M) away from the boundary. We use the closed-form formulae for g and dx using our approximation for dg (.8, 1.2) described above. 2volg (M) L(εN ,N ) f (x) where L(εN ,N ) f is constructed using the m/2+1 m2 N εN squared Euclidean distance D(x, y) := (x − y)2 instead of S1 (x, y). Since we are sampling from the
For comparison, we also plot x 7→
probability measure on [0.8, 1.2] with density P as in (49), we expect that the limiting function is instead 2volg (M) m2 d2 f df dP 1 df dg √ d2 f −P 2 − 2 ) =− g 2 −√ . m2 2 dx dx dx dx g dx dx
This is confirmed in the second subplot of Fig. 5. (1) (2) In Fig. 6, we map the ith sample to (vi , vi ), where v (1) , v (2) are again the first two nonconstant eigenvectors of the normalized graph Laplacian (D(εN ,N ) )−1 L(εN ,N ) constructed with D = S1 . Each point is colored according to the value of xi . For manifolds with boundary, we expect that for x1 , . . . , xN sampled i.i.d. uniformly from (M, g), with high probability the eigenvectors of the graph Laplacian (appropriately scaled) converge to the Neumann eigenfunctions of the LaplaceBeltrami operator as N → ∞, as is the case when (M, g) is a Riemannian submanifold of Rd and the squared Euclidean distance is used in the construction of the graph Laplacian [27, 73]. We also observe in Fig. 6 that the eigenvector v (1) gives us an ordering of the samples. Indeed, this is very similar to the spectral algorithm for seriation proposed by Atkins et al. [7], which uses a Fiedler eigenvector (an eigenvector with the smallest eigenvalue among those orthogonal to the all ones vector) of the unnormalized graph Laplacian to order samples. See [84] for an example of seriation using a graph Laplacian constructed from the squared Wasserstein-2 distance.
27
6
Conclusion and future work
In this work, we observe that a smooth symmetric divergence D : M × M → R≥0 satisfying Assumption 2.2 gives us an approximation of the squared geodesic distance d2g (for g as in (3)) up to fourth order, in the sense that there exists K > 0 such that D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4 for all p, q ∈ M. This is a similar error bound to when (M, g) is a Riemannian submanifold of Rd and D is given by the squared Euclidean distance, which is the typical setting for manifold learning. As an example of the kinds of results we expect to carry over from the Euclidean setting, we prove the pointwise convergence of graph Laplacians. One motivation for this work was to better understand manifold learning with Sinkhorn divergences. To this end, we describe two models (Assumptions 4.1 and 4.2) and discuss how to check for the smoothness and non-degeneracy required for Assumption 2.2 under these models. We also study two concrete examples in Section 5. However, in our analysis, we fix the regularization parameter β > 0, and we assume that we can compute the Sinkhorn divergence exactly. One direction for future work would be to study the dependence of various quantities (e.g., the Riemannian metric, K from Proposition 3.1) on β; see [49] for an informal argument, for measures with density, on the convergence as β → 0 of the Riemannian metric induced by Sβ to the Riemannian structure of Wasserstein-2 space. Another direction would be to do a more detailed analysis of manifold learning algorithms when the Sinkhorn divergences are computed via Sinkhorn’s algorithm (or a variant thereof), which is used often in practice.
Acknowledgments This work was supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-2444107 and the Simons Foundation Math+X Investigator Award to Amit Singer. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation. We thank Amit Singer for comments on an earlier draft and Gilles Mordant for helpful discussions.
References [1] S.-I. Amari. Natural gradient works efficiently in learning. Neural Computation, 10(2):251–276, 1998. [2] S.-I. Amari. Information geometry and its applications, volume 194. Springer, 2016. [3] L. Ambrosio, N. Gigli, and G. Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media, 2008. [4] E. Arias-Castro and W. Qiao. Embedding distributional data. The Annals of Statistics, 53(2):615–646, 2025. [5] N. Aronszajn. Theory of Reproducing Kernels. Transactions of the American Mathematical Society, 68(3):337–404, 1950. 28
[6] U. M. Ascher and C. Greif. A First Course in Numerical Methods. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2011. [7] J. E. Atkins, E. G. Boman, and B. Hendrickson. A spectral algorithm for seriation and the consecutive ones problem. SIAM Journal on Computing, 28(1):297–310, 1998. [8] N. Ay, J. Jost, H. Vân Lê, and L. Schwachhöfer. Information geometry, volume 64. Springer, 2017. [9] J. Bates. The embedding dimension of Laplacian eigenfunction maps. Applied and Computational Harmonic Analysis, 37(3):516–530, 2014. [10] M. Belkin and P. Niyogi. Laplacian eigenmaps and spectral techniques for embedding and clustering. In T. Dietterich, S. Becker, and Z. Ghahramani, editors, Advances in Neural Information Processing Systems, volume 14. MIT Press, 2001. [11] M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373–1396, 2003. [12] M. Belkin and P. Niyogi. Towards a theoretical foundation for Laplacian-based manifold methods. In International Conference on Computational Learning Theory, pages 486–500. Springer, 2005. [13] M. Belkin and P. Niyogi. Convergence of Laplacian eigenmaps. Advances in Neural Information Processing Systems, 19, 2006. [14] M. Belkin and P. Niyogi. Towards a theoretical foundation for Laplacian-based manifold methods. Journal of Computer and System Sciences, 74(8):1289–1308, 2008. [15] J.-D. Benamou and Y. Brenier. A computational fluid mechanics solution to the MongeKantorovich mass transfer problem. Numerische Mathematik, 84(3):375–393, 2000. [16] P. Bérard, G. Besson, and S. Gallot. Embedding Riemannian manifolds by their heat kernel. Geometric & Functional Analysis GAFA, 4:373–398, 1994. [17] M. Bernstein, V. De Silva, J. C. Langford, and J. B. Tenenbaum. Graph approximations to geodesics on embedded manifolds. Technical report, Citeseer, 2000. [18] J. Bigot, R. Gouet, T. Klein, and A. Lopez. Geodesic PCA in the Wasserstein space by convex PCA. In Annales de l’Institut Henri Poincaré (B) Probabilités et Statistiques, volume 53, pages 1–26, 2017. [19] J. Calder, N. García Trillos, and M. Lewicka. Lipschitz regularity of graph Laplacians on random data clouds. SIAM Journal on Mathematical Analysis, 54(1):1169–1222, 2022. [20] J. Calder and N. G. Trillos. Improved spectral convergence rates for graph Laplacians on ε-graphs and k-NN graphs. Applied and Computational Harmonic Analysis, 60:123–175, 2022. [21] G. Carlier, L. Chizat, and M. Laborde. Displacement smoothness of entropic optimal transport. ESAIM: Control, Optimisation and Calculus of Variations, 30:25, 2024.
29
[22] K. M. Carter, R. Raich, W. G. Finn, and A. O. Hero III. FINE: Fisher information nonparametric embedding. IEEE Transactions on Pattern Analysis and Machine Intelligence, 31(11):2093–2098, 2009. [23] Y. Chen and W. Li. Optimal transport natural gradient for statistical manifolds with continuous sample space. Information Geometry, 3(1):1–32, 2020. [24] X. Cheng and H.-T. Wu. Convergence of graph Laplacian with kNN self-tuned kernels. Information and Inference: A Journal of the IMA, 11(3):889–957, 2022. [25] X. Cheng and N. Wu. Eigen-convergence of Gaussian kernelized graph Laplacian by manifold heat interpolation. Applied and Computational Harmonic Analysis, 61:132–190, 2022. [26] S. Chewi, J. Niles-Weed, and P. Rigollet. Statistical optimal transport, volume 2364 of Lecture Notes in Mathematics. Springer, Cham, [2025] ©2025. École d’Été de Probabilités de SaintFlour XLIX – 2019. [27] R. R. Coifman and S. Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5–30, 2006. [28] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker. Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps. Proceedings of the National Academy of Sciences, 102(21):7426–7431, 2005. [29] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013. [30] W. Diepeveen, C. Esteve-Yagüe, J. Lellmann, O. Öktem, and C.-B. Schönlieb. Riemannian geometry for efficient analysis of protein dynamics data. Proceedings of the National Academy of Sciences, 121(33):e2318951121, 2024. [31] M. P. Do Carmo and J. Flaherty Francis. Riemannian geometry, volume 2. Springer, 1992. [32] D. B. Dunson, H.-T. Wu, and N. Wu. Spectral convergence of graph Laplacian and heat kernel reconstruction in L∞ from random samples. Applied and Computational Harmonic Analysis, 55:282–336, 2021. [33] C. Fefferman, S. Mitter, and H. Narayanan. Testing the manifold hypothesis. Journal of the American Mathematical Society, 29(4):983–1049, 2016. [34] J. Feydy, T. Séjourné, F.-X. Vialard, S.-i. Amari, A. Trouvé, and G. Peyré. Interpolating between Optimal Transport and MMD using Sinkhorn Divergences. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2681–2690. PMLR, 2019. [35] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer. POT: Python Optimal Transport. Journal of Machine Learning Research, 22(78):1–8, 2021.
30
[36] R. Flamary, C. Vincent-Cuaz, N. Courty, A. Gramfort, O. Kachaiev, H. Quang Tran, L. David, C. Bonet, N. Cassereau, T. Gnassounou, E. Tanguy, J. Delon, A. Collas, S. Mazelet, L. Chapel, T. Kerdoncuff, X. Yu, M. Feickert, P. Krzakala, T. Liu, and E. Fernandes Montesuma. POT Python Optimal Transport (version 0.9.5), 2024. [37] N. García Trillos, M. Gerlach, M. Hein, and D. Slepčev. Error estimates for spectral convergence of the graph Laplacian on random geometric graphs toward the Laplace–Beltrami operator. Foundations of Computational Mathematics, 20(4):827–887, 2020. [38] A. Genevay, L. Chizat, F. Bach, M. Cuturi, and G. Peyré. Sample Complexity of Sinkhorn Divergences. In The 22nd international conference on artificial intelligence and statistics, pages 1574–1583. PMLR, 2019. [39] A. Genevay, G. Peyré, and M. Cuturi. Learning Generative Models with Sinkhorn Divergences. In International Conference on Artificial Intelligence and Statistics, pages 1608–1617. PMLR, 2018. [40] A. Gonzalez-Sanz, J.-M. Loubes, and J. Niles-Weed. Weak limits of entropy regularized optimal transport; potentials, plans and divergences. arXiv preprint arXiv:2207.07427, 2024. [41] K. Hamm, N. Henscheid, and S. Kang. Wassmap: Wasserstein isometric mapping for image manifold learning. SIAM Journal on Mathematics of Data Science, 5(2):475–501, 2023. [42] K. Hamm, C. Moosmüller, B. Schmitzer, and M. Thorpe. Manifold Learning in Wasserstein Space. SIAM Journal on Mathematical Analysis, 57(3):2983–3029, 2025. [43] M. Hardion and H. Lavenant. Gradient flows of potential energies in the geometry of sinkhorn divergences. arXiv preprint arXiv:2511.14278, 2025. [44] M. Hein, J.-Y. Audibert, and U. Von Luxburg. From graphs to manifolds–weak and strong pointwise consistency of graph Laplacians. In International Conference on Computational Learning Theory, pages 470–485. Springer, 2005. [45] L. Hörmander. The analysis of linear partial differential operators. I. Classics in Mathematics. Springer-Verlag, Berlin, 2003. Distribution theory and Fourier analysis, Reprint of the second (1990) edition [Springer, Berlin; MR1065993 (91m:35001a)]. [46] J. Kileel, A. Moscovich, N. Zelesko, and A. Singer. Manifold learning with arbitrary norms. Journal of Fourier Analysis and Applications, 27(5):82, 2021. [47] S. Kullback. Information theory and statistics. Courier Corporation, 1997. [48] S. Lang. Fundamentals of differential geometry, volume 191 of Graduate Texts in Mathematics. Springer-Verlag, New York, 1999. [49] H. Lavenant, J. Luckhardt, G. Mordant, B. Schmitzer, and L. Tamanini. The Riemannian geometry of Sinkhorn divergences. Annales de l’Institut Henri Poincaré C, Analyse non linéaire, Oct. 2025. [50] W. Li and G. Montúfar. Natural gradient via optimal transport. Information Geometry, 1(2):181–214, 2018. 31
[51] W. Li and J. Zhao. Wasserstein information matrix. Information Geometry, 6(1):203–255, 2023. [52] J. Lott. Some Geometric Calculations on Wasserstein Space. Commun. Math. Phys, 277:423– 437, 2008. [53] P. Y. Lu, R. Dangovski, and M. Soljačić. Discovering conservation laws using optimal transport and manifold learning. Nature Communications, 14(1):4744, 2023. [54] G. Luise, A. Rudi, M. Pontil, and C. Ciliberto. Differential properties of Sinkhorn approximation for learning with Wasserstein distance. Advances in Neural Information Processing Systems, 31, 2018. [55] G. Mena and J. Niles-Weed. Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem. Advances in Neural Information Processing Systems, 32, 2019. [56] C. A. Micchelli, Y. Xu, and H. Zhang. Universal kernels. Journal of Machine Learning Research, 7(12), 2006. [57] G. Mishne, R. Talmon, R. Meir, J. Schiller, M. Lavzin, U. Dubin, and R. R. Coifman. Hierarchical coupled-geometry analysis for neuronal structure and activity pattern discovery. IEEE Journal of Selected Topics in Signal Processing, 10(7):1238–1253, 2016. [58] K. Muandet, K. Fukumizu, B. Sriperumbudur, and B. Schölkopf. Kernel mean embedding of distributions: A review and beyond. Foundations and Trends® in Machine Learning, 10(12):1–141, 2017. [59] B. Nadler, S. Lafon, R. R. Coifman, and I. G. Kevrekidis. Diffusion maps, spectral clustering and reaction coordinates of dynamical systems. Applied and Computational Harmonic Analysis, 21(1):113–127, 2006. [60] M. Nutz. Introduction to entropic optimal transport. Lecture notes, Columbia University, 2022. [61] M. C. A. Oliver, M. Roberts, C.-B. Schönlieb, and M. Thorpe. Laplace learning in Wasserstein space. arXiv preprint arXiv:2511.13229, 2025. [62] F. Otto. The geometry of dissipative evolution equations: the porous medium equation. Communications in Partial Differential Equations, 26(1-2):101–174, 2001. [63] G. Peyré and M. Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019. [64] C. R. Rao. Differential metrics in probability spaces. Differential geometry in statistical inference, 10:217–240, 1987. [65] M. Rumpf and B. Wirth. Variational time discretization of geodesic calculus. IMA Journal of Numerical Analysis, 35(3):1011–1046, 2015. [66] V. Seguy and M. Cuturi. Principal geodesic analysis for probability measures under the optimal transport metric. Advances in Neural Information Processing Systems, 28, 2015. 32
[67] Z. Shen, Z. Wang, A. Ribeiro, and H. Hassani. Sinkhorn natural gradient for generative models. Advances in Neural Information Processing Systems, 33:1646–1656, 2020. [68] S. Shirdhonkar and D. W. Jacobs. Approximate earth mover’s distance in linear time. In 2008 IEEE Conference on Computer Vision and Pattern Recognition, pages 1–8. IEEE, 2008. [69] C.-J. Simon-Gabriel and B. Schölkopf. Kernel distribution embeddings: Universal kernels, characteristic kernels and kernel metrics on distributions. Journal of Machine Learning Research, 19(44):1–29, 2018. [70] C.-J. Simon-Gabriel and B. Schölkopf. Kernel distribution embeddings: Universal kernels, characteristic kernels and kernel metrics on distributions. arXiv preprint arXiv:1604.05251, 2019. [71] A. Singer. From graph to manifold Laplacian: The convergence rate. Applied and Computational Harmonic Analysis, 21(1):128–134, 2006. [72] A. Singer. Wilson statistics: derivation, generalization and applications to electron cryomicroscopy. Foundations of Crystallography, 77(5):472–479, 2021. [73] A. Singer and H.-T. Wu. Spectral convergence of the connection Laplacian from random samples. Information and Inference: A Journal of the IMA, 6(1):58–123, 2017. [74] R. Sinkhorn. A relationship between arbitrary positive matrices and doubly stochastic matrices. The Annals of Mathematical Statistics, 35(2):876–879, 1964. [75] A. Smola, A. Gretton, L. Song, and B. Schölkopf. A Hilbert space embedding for distributions. In International Conference on Algorithmic Learning Theory, pages 13–31. Springer, 2007. [76] B. Sriperumbudur, K. Fukumizu, and G. Lanckriet. On the relation between universality, characteristic kernels and RKHS embedding of measures. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 773–780. JMLR Workshop and Conference Proceedings, 2010. [77] I. Steinwart. On the influence of the kernel on the consistency of support vector machines. Journal of Machine Learning Research, 2(Nov):67–93, 2001. [78] I. Steinwart and A. Christmann. Support vector machines. Springer Science & Business Media, 2008. [79] J. B. Tenenbaum, V. d. Silva, and J. C. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500):2319–2323, 2000. [80] N. G. Trillos, C. Li, and R. Venkatraman. Minimax Rates for the Estimation of Eigenpairs of Weighted Laplace-Beltrami Operators on Manifolds. arXiv preprint arXiv:2506.00171, 2025. [81] N. G. Trillos, A. Little, D. McKenzie, and J. M. Murphy. Fermat distances: Metric approximation, spectral convergence, and clustering algorithms. Journal of Machine Learning Research, 25(176):1–65, 2024. [82] C. Villani et al. Optimal transport: old and new, volume 338. Springer, 2009. 33
[83] M. Wahl. A kernel-based analysis of Laplacian Eigenmaps. arXiv preprint arXiv:2402.16481, 2024. [84] A. Warren, A. Afanassiev, F. Kobayashi, Y.-H. Kim, and G. Schiebinger. Principal curves in metric spaces and the space of probability measures. arXiv preprint arXiv:2505.04168, 2025. [85] H. Whitney. Functions differentiable on the boundaries of regions. Annals of Mathematics, 35(3):482–485, 1934. [86] L. Xu and A. Singer. Manifold learning in metric spaces. Applied and Computational Harmonic Analysis, page 101813, 2025. [87] N. Zelesko, A. Moscovich, J. Kileel, and A. Singer. Earthmover-based manifold learning for analyzing molecular conformation spaces. In 2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI), pages 1715–1719. IEEE, 2020.
A
Proofs for Section 3
Proof of Proposition 3.1. It suffices to prove that there exists K, r > 0 such that (50)
D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4
for all p, q ∈ M such that dg (p, q) < r. Indeed, since we assumed that D is continuous and M is compact, D(p, q) − dg (p, q)2 sup < ∞. dg (p, q)4 p,q∈M,dg (p,q)≥r First consider any unit-speed, length-minimizing geodesic segment γ : [−δ, δ] → M on (M, g). Let f (s, t) := D(γ(s), γ(t)) for s, t ∈ [−δ, δ], and for each s ∈ [−δ, δ], let fs := f (s, ·). Assume δ is small enough so that f ∈ C 4 ([−δ, δ] × [−δ, δ]). By Assumption 2.2, f is symmetric, f ≥ 0 and f (s, s) = 0 for all s ∈ [−δ, δ]. Additionally, for any fixed s ∈ (−δ, δ), lim
f (s, t)
t→s (t − s)2
= lim
t→s
D(γ(s), γ(t)) = gγ(s) (γ ′ (s), γ ′ (s)) = 1 (t − s)2
since γ is a unit-speed geodesic segment and we defined g via gp =
1 Hessp (D(p, ·)) 2
for all p ∈ M. Thus fs′′ (s) = 2 for all s ∈ (−δ, δ), so Lemma 3.1 implies that (4)
2
D(γ(0), γ(t)) − dg (γ(0), γ(t))
= D(γ(0), γ(t)) − t
2
≤
sups∈(−δ,δ) f0 (s) 24
dg (γ(0), γ(t))4 (51)
for all t ∈ (−δ, δ). Taking q = γ(0) and v = γ ′ (0), we have f0 (t) = D(γ(0), γ(t)) = D(q, expq (tv)) for t ∈ (−δ, δ). Defining Gq,v (t) := D(q, expq (tv)) for (q, v) ∈ T M, it thus suffices to show that there exists K, r > 0 such that for all q ∈ M and v ∈ Tq M such that gq (v, v) = 1, 34
1. s, t 7→ D(expq (sv), expq (tv)) belongs to C 4 ([−r, r] × [−r, r]), and (4)
2. Gq,v (t) < K for all t ∈ (−r, r). This follows from the compactness of M and is almost the same as [86, Lemma 3.5]. We provide the details for completeness. First fix any p ∈ M. Define F : T M → M × M by F (q, v) := (q, expq (v)). (As M is compact, F is well-defined on all of T M by Hopf-Rinow; if M were not compact, F is still well-defined on a neighborhood of (p, 0).) Then by the inverse function theorem, there exists an open neighborhood U of (p, 0) and an open neighborhood V of (p, p) such that F |U : U → V is a diffeomorphism [31, Theorem 3.7]. Making U smaller if necessary, we can also assume D is smooth on V , so D ◦ F is smooth on U . For small enough r(p) > 0, the geodesic ball B(p, r(p)) := {q ∈ M : dg (p, q) < r(p)} is strongly convex normal neighborhood of p (see [31, Ch. 3]) and B(p, r(p)) × B(p, r(p)) ⊂ V . h i r(p) r(p) For each q ∈ B p, r(p) , v ∈ T M such that g (v, v) = 1 and t ∈ − , , recall that by q q 2 2 2 definition we have Gq,v (t) := D(F (q, tv)) = D(q, expq (tv)). As B(p, r(p)) is strongly convex, expq (tv) ∈ B(p, r(p)) for all q ∈ B p, r(p) , v ∈ Tq M such that 2 h i r(p) r(p) gq (v, v) = 1 and t ∈ − 2 , 2 , and so for all such q, v, t, we have (q, tv) ∈ F −1 B(p, r(p)) × B(p, r(p)) ⊂ U. As D ◦ F is smooth on U and F −1 B(p, r(p)) × B(p, r(p)) is compact, there exists Kp > 0 such that the fourth derivative G(4) q,v (t) < Kp r(p) r(p) for all q ∈ B p, r(p) , v ∈ T M such that g (v, v) = 1 and t ∈ − . , q q 2 2 2 As M is compact, one can cover M with finitely many balls of the form B p, r(p) , and so 2 there exist K, r > 0 such that s, t 7→ D(expq (sv), expq (tv)) belongs to C 4 ([−r, r] × [−r, r]) and (4)
Gq,v (t) < K for all t ∈ (−r, r) for all (q, v) in the unit tangent bundle of M. Proof of Lemma 3.1. For any x, y ∈ (a, b), we have the Taylor expansions f (x, y) =
1 ∂2f 1 ∂3f 1 ∂4f 2 3 (x, x)(y − x) + (x, x)(y − x) + (x, s)(y − x)4 2 ∂y 2 6 ∂y 3 24 ∂y 4
(52)
f (y, x) =
1 ∂2f 1 ∂3f 1 ∂4f (y, y)(x − y)2 + (y, y)(x − y)3 + (y, t)(x − y)4 2 3 2 ∂y 6 ∂y 24 ∂y 4
(53)
35
for some s, t between x and y. Since f is symmetric, we can equate (52) and (53) to obtain ∂2f ∂2f ∂y 2 (x, x) − ∂y 2 (y, y)
2(y − x)
1 + 6
∂3f 1 ∂4f ∂4f ∂3f (x, x) + 3 (y, y) + (x, s) − 4 (y, t) (y − x) = 0 ∂y 3 ∂y 24 ∂y 4 ∂y
for any x, y ∈ (a, b) such that x ̸= y. Fixing x ∈ (a, b) and taking y → x, we then have 3 1 ∂3f 1 ∂3f ∂ f − (x, x) + (x, x) + (x, x) = 0, 2 ∂x∂y 2 ∂y 3 3 ∂y 3 which implies (22) upon rearranging. To prove (23), observe that if we have fx′′ (x) = for all x ∈ (a, b), then
∂2f (x, x) = 2 ∂y 2
∂3f ∂3f (x, x) + (x, x) = 0 ∂x∂y 2 ∂y 3
for all x ∈ (a, b). Using (22), we obtain −
1 ∂3f ∂3f (x, x) + (x, x) = 0 3 ∂y 3 ∂y 3
for all x ∈ (a, b), so fx′′′ (x) = 0 for all x ∈ (a, b) and (23) follows from Taylor’s theorem with the Lagrange remainder.
A.1
Pointwise convergence of the graph Laplacian (Theorem 3.1)
Throughout this section, assume that D, g are given by Assumption 2.2 and that h satisfies Assumption 2.1. For each ε > 0 and continuous f : M → R, we follow [27] and define Gε f by Z 1 Gε f (p) := m/2 Kε (p, q)f (q)dVg (q) ε M for each p ∈ M, where dVg denotes the Riemannian volume form for g and we recall from (4) that Kε : M × M → R is given by D(p, q) Kε (p, q) := h . ε Using the notation from [27, Lemma 8], we have the following proposition: Proposition A.1. Fix any f ∈ C 3 (M). Then for each p ∈ M, Gε f (p) = m0 f (p) +
εm2 (ω(p)f (p) − ∆g f (p)) + O(ε3/2 ), 2
as ε → 0+ , where ω : M → R depends on g and D, and Z m0 := h(∥v∥2 )dv Rm
36
(54)
Z m2 := Rm
v12 h(∥v∥2 )dv.
If f ∈ C 4 (M), the error can be improved to O(ε2 ). Moreover, the multiplicative constant for the O(ε3/2 ) (O(ε2 ) if f ∈ C 4 (M)) error can be made uniform across p ∈ M. This proposition is well-known when (M, g) is a compact Riemannian submanifold of Rd and D is the squared Euclidean distance [27, Lemma 8] (see also [24, Lemma A.5]). Proposition A.1 allows us to control the bias of the discrete Laplacian L(εN ,N ) f , and so Theorem 3.1 follows from standard concentration inequalities. Proof of Theorem 3.1. We first show that L(εN ,N ) f concentrates around its mean. This step is the same as if (M, g) were a compact Riemannian submanifold of Euclidean space and D were the squared Euclidean distance, as it does not depend on any particular properties of the Euclidean distance (see [46, Section 4.1] or [86, Section 3]). For completeness, we briefly review the argument: for (ε) (ε) each ε > 0 and j ∈ N, upon defining the random variables Xj := Kε (x, xj ), Yj := Kε (x, xj )f (xj ), we have N X (ε) (ε) L(ε,N ) f (x) = f (x)Xj − Yj , j=1 (ε)
and we can apply Hoeffding’s inequality since f (x)X1 (ε) (ε) f (x)Xj − Yj
(ε)
(ε)
(ε)
− Y1 , ..., f (x)XN − YN
are i.i.d. and
≤ 2c0 ∥f ∥∞ almost surely for any j ∈ N. This gives us P
2 h i − 2N δ 2 L(ε,N ) f (x) (ε) (ε) > δ ≤ 2e 8c0 ∥f ∥∞ − E f (x)X1 − Y1 N 1
for any δ > 0. As εN = Ω(N − m+2+α ) for some α > 0, we have X
(εN ,N )
1 P
N ∈N
L
m/2+1 εN
N
h i f (x) (ε ) (ε ) >δ − E f (x)X1 N − Y1 N
! ≤
X
2e
−
m/2+1 2 N δεN 8c20 ∥f ∥2∞
< ∞,
N ∈N
for any δ > 0, so one can conclude 1 P
m/2+1
εN
! h i L(εN ,N ) f (x) n→∞ (εN ) (εN ) − E f (x)X1 − Y1 −−−−→ 0 = 1 N
by using Borel-Cantelli and observing that δ > 0 was arbitrary. It thus remains to prove that h i m 2 (εN ) (εN ) E f (x)X − Y = P (x)∆g f (x) − 2gx (gradg f (x), gradg P (x)) . 1 1 N →∞ εm/2+1 2 1
lim
N
37
h i (ε ) (ε ) This follows from expanding the definition of E f (x)X1 N − Y1 N and using Proposition A.1: 1
1
εN
εN
h i (εN ) (εN ) E f (x)X − Y = 1 1 m/2
m/2
Z f (x)
Z KεN (x, y)P (y)dVg (y) −
f (y)KεN (x, y)P (y)dVg (y)
= f (x)GεN P (x) − GεN [f P ](x) εN m2 3/2 = f (x) m0 P (x) + (ω(x)P (x) − ∆g P (x)) + O εN − 2 εN m2 3/2 m0 f (x)P (x) + (ω(x)f (x)P (x) − ∆g (f P )(x)) + O εN 2 εN m 2 3/2 P (x)∆g f (x) − 2gx (gradg f (x), gradg P (x)) + O εN . = 2 It remains to prove Proposition A.1. It is almost the same as the arguments in [14, 27]; again we provide the details for completeness. Proof of Proposition A.1. Let U be an open neighborhood of the diagonal ΞM ⊂ M × M such that D is smooth on U . Let inj(M, g) denote the injectivity radius of (M, g), and let K be a positive number satisfying D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4 (55) for all p, q ∈ M. Such a K > 0 exists by Proposition 3.1. Recall that we assumed that M is compact, so we also have inj(M, g) > 0. Choose !! r 1 r ∈ 0, min inj(M, g), inf c dg (p, q), 2K (p,q)∈U if the complement U c of U is non-empty, and r r∈
0, min inj(M, g),
1 2K
!!
if U = M × M. 1. Approximation of Gε f with an integral over a geodesic ball. For any p ∈ M, ε ∈ (0, r) and f ∈ C 3 (M), we can approximate Gε f (p) with an integral over the geodesic ball B(p, r) of radius r: Gε f (p) −
1 εm/2
Z Kε (p, q)f (q)dVg (q) ≤
volg (M \ B(p, r)) supq∈M\B(p,r) |Kε (p, q)f (q)| εm/2
B(p,r)
Since we assumed that D is continuous and D(p, q) = 0 if and only if p = q, D∗ :=
inf
D(p, q)
p,q∈M : dg (p,q)≥r
is strictly greater than 0. Per our assumptions on h (Assumption 2.1), cD(p,q) cD ∗ D(p, q) 0 ≤ Kε (p, q) = h ≤ c0 e− ε ≤ c0 e− ε ε 38
.
for all p, q ∈ M such that dg (p, q) ≥ r. Hence Gε f (p) −
cD ∗
c0 volg (M)e− ε ∥f ∥∞ Kε (p, q)f (q)dVg (q) ≤ εm/2 B(p,r)
Z
1 εm/2
(56)
for any p ∈ M. 2. Taylor expansions and integration in normal coordinates. Fix any p ∈ M. We R 2 now approximate ε−m/2 B(p,r) Kε (p, q)f (q)dVg (q) with m0 f (p) + εm 2 (ω(p)f (p) − ∆g f (p)), where ω : M → R is to be determined. For this, we use normal coordinates around p and Taylor expansions. For concreteness, we can fix an orthonormal basis {e1 , . . . , em } of Tp M to identify Tp M with Rm . For any q = expp (v), q ̸= p with ∥v∥ = dg (p, q) < r, h
D(p, q) ε
=h
∥v∥2 ε
+
D(p, q) − ∥v∥2 ε
h′
∥v∥2 ε
+
1 2
D(p, q) − ∥v∥2 ε
2
h′′ (s)
(57)
for some s between D(p, q)/ε and ∥v∥2 /ε, and moreover D(p, q) − ∥v∥2 = D(p, q) − dg (p, q)2 ≤ Kdg (p, q)4 = K∥v∥4 . 1 By the choice of r, we have r2 ≤ 2K , so the right-hand side can be bounded by
K∥v∥4 ≤ Kr2 ∥v∥2 ≤
∥v∥2 2
for all v ∈ Tp M such that ∥v∥ < r. By the exponential decay of the second derivative of h in Assumption 2.1, for v, s as in (57), |h′′ (s)| ≤ c2 e−cs ≤ c2 e− Since
Z
∥v∥k e−
(
c ∥v∥2 −K∥v∥4
α∥v∥2 ε
ε
)
c∥v∥2
≤ c2 e− 2ε .
m+k dv = O ε 2
(58)
Rm
√ for any fixed α > 0 and k ∈ N (use a change of variables v 7→ v/ ε), we have
Z ∥v∥<r
D(p, q) − ∥v∥2 ε
2
c∥v∥2
e− 2ε dv ≤
K2 ε2
Z
m+4 c∥v∥2 ∥v∥8 e− 2ε dv = O ε 2 ,
Rm
so that approximating h with its first-order Taylor expansion results in an O(ε2 )-error: 1 εm/2
Z Kε (p, q)f (q)dVg (q) B(p,r)
r h i D(p, expp (v)) f expp (v) det gexpp (v) dv = m/2 h ε ε ∥v∥<r ! r ! Z h i D(p, expp (v)) − ∥v∥2 ∥v∥2 ∥v∥2 1 ′ = m/2 h + h f (expp (v)) det gexpp (v) dv + O(ε2 ). ε ε ε ε ∥v∥<r 1
Z
39
h i Here we use the shorthand gexpp (v) for the matrix h i ∂ ∂ gexpp (v) := gexpp (v) , ∂vi ∂vj i,j for each v ∈ Tp M, ∥v∥ ≤ r. We can also Taylor expand v 7→ f (expp (v)), v 7→
r
h i det gexpp (v) and
v 7→ D(p, expp (v)) around 0. Following the notation in [27, Appendix B], we use Qp,k to denote a generic homogeneous polynomial of degree k whose coefficients can depend on p (i.e., the particular polynomial may change from line to line). Then as v → 0, we have 1 f (expp (v)) = f (p) + ⟨gradg f (p), v⟩ + Hessp f (v, v) + O(∥v∥3 ) 2 r h i Ricp (v, v) det gexpp (v) = 1 − + Qp,3 (v) + O(∥v∥4 ) 6 D(p, expp (v)) = ∥v∥2 + Qp,4 (v) + Qp,5 (v) + O(∥v∥6 )
(59) (60) (61)
for v ∈ Tp M such that ∥v∥ < r. Here Ricp denotes the Ricci curvature tensor at p. Notably, v 7→ D(p, expp (v)) has no third-order term in its Taylor expansion because of Proposition 3.1. If we additionally have f ∈ C 4 (M), then we have 1 f (expp (v)) = f (p) + ⟨gradg f (p), v⟩ + Hessp f (v, v) + Qp,3 (v) + O(∥v∥4 ). 2
(62)
Observing that odd functions integrate to 0 and using the exponential decay of h, h′ from Assumption 2.1 with (58), we obtain
Z
1
h
εm/2
∥v∥<r
=
∥v∥2 ε
f (expp (v))
Z
1 εm/2
r
h ∥v∥<r
∥v∥2 ε
h i det gexpp (v) dv
Ricp (v, v) 1 f (p) 1 − + Hessp f (v, v) dv + O(ε3/2 ), 6 2
where the error can be improved to O(ε2 ) if f ∈ C 4 (M) and we can use (62), and also 1 εm/2
! r h i D(p, expp (v)) − ∥v∥2 ∥v∥2 ′ h f (expp (v)) det gexpp (v) dv ε ε ∥v∥<r Z Qp,4 (v) ∥v∥2 1 = m/2 h′ f (p)dv + O(ε2 ). ε ε ε ∥v∥<r
Z
Again using the exponential decay of h and h′ , we can approximate the above integrals over {v ∈ Tp M : ∥v∥ < r} with their respective integrals over Rm with an error smaller than O(ε2 ), so that altogether we obtain 1 εm/2
Z Kε (p, q)f (q)dVg (q) B(p,r)
f (p)Qp,4 (v)h′ ∥v∥2 ε ∥v∥ Ricp (v, v) 1 1 = m/2 h f (p) 1 − + Hessp f (v, v) + dv + O(ε3/2 ) ε 6 2 ε ε Rm Z Z −h(∥v∥2 )Ricp (v, v) ε ′ = f (p) m0 + ε + Qp,4 (v)h (∥v∥)dv + h(∥v∥2 )Hessp f (v, v)dv + O(ε3/2 ) 6 2 m m R R Z
2
40
√ upon the change of variables v 7→ v/ ε, where again the error can be improved to O(ε2 ) if f ∈ C 4 (M). (54) now follows from the equality Z h(∥v∥2 )Hessp f (v, v)dv = −m2 ∆g f (p), Rm
which holds from diagonalizing Hessp f and using ∆g f (p) = −Tr(Hessp f ). The multiplicative constant for the error can be made uniform over p ∈ M by observing that the dependence of the error on p is in terms of finitely many derivatives of r h i v 7→ f (expp (v)),
det gexpp (v) , D(p, expp (v)),
which can be uniformly bounded on the compact set {(p, v) ∈ T M : ∥v∥ ≤ r}.
B
Proofs for Section 4
For this section, we use ε to denote an infinitesimal, not the bandwidth parameter for the graph Laplacian in Section 2.1. Proof of Lemma 4.2. For convenience, let kβ (x, y) := exp −∥x − y∥2 /β . Similar to [49, Lemma 3.13], it suffices to prove that Aβ : C k (X) × U → C k (X) given by Z [Aβ (g, p)](y) := g(x)kβ (x, y)dµp (x) is C k . The proof of [49, Lemma 3.13] establishes the continuous differentiability of Aβ when U = (−δ, δ) for some δ > 0 and {µt }t∈(−δ,δ) is what they call a C k -perturbation (continuously differentiable in C k (X)∗ with the weak-* topology). Although we could use similar arguments to establish the continuous differentiability of Aβ in our setting, under Assumption 4.2, a more direct proof is possible. Moreover, since this approach generalizes to higher order derivatives and elucidates the dependence of the regularity of Aβ on the regularity of Ψ with respect to the first variable, we will take this approach. To start, let us define Hg (p, x, y) := g(Ψ(p, x))kβ (Ψ(p, x), y) (63) Hg,α (p, x, y) := g(Ψ(p, x))∂yα kβ (Ψ(p, x), y) = ∂yα Hg (p, x, y)
(64)
for each p ∈ U , x, y ∈ X, g ∈ C k (X) and multi-index α with |α| ≤ k, so that Z [Aβ (g, p)](y) = Hg (p, x, y)dµ(x) for each g ∈ C k (X) and p ∈ U . One can then check via the dominated convergence theorem that Aβ (g, p) is indeed in C k (X) for each g ∈ C k (X) and p ∈ U , with Z Z α α ∂ [Aβ (g, p)](y) = g(Ψ(p, x))∂y kβ (Ψ(p, x), y)dµ(x) = Hg,α (p, x, y)dµ(x) (65) 41
for each multi-index α with |α| ≤ k. Continuity of Aβ . We first show that Aβ is continuous from C k (X) × U to C k (X). Suppose gn → g in C k (X) and pn → p in U . For any multi-index α with |α| ≤ k, since Hg,α is continuous on U × X × X, it is uniformly continuous on Br (p) × X × X for some small enough r > 0. Hence for any ε > 0, there exists δ > 0 such that whenever |pn − p| < δ, ∥Aβ (g, pn ) − Aβ (g, p)∥C k (X) ≤ ε. Moreover, for any n ∈ N, one can bound ∥Aβ (gn , pn ) − Aβ (g, pn )∥C k (X) ≤ ∥gn − g∥C 0 (X)
X |α|≤k
sup ∂yα kβ (x, y) ,
x,y∈X
so ∥Aβ (gn , pn ) − Aβ (g, p)∥C k (X) → 0 as n → ∞. If k ≥ 1, by [48, Proposition 3.5], to show that Aβ is C k , it suffices to check that the partial derivatives D1 Aβ : C k (X) × U → B(C k (X), C k (X)) D2 Aβ : C k (X) × U → B(Rm , C k (X)) exist and are C k−1 . Existence and continuity of D1 Aβ when k ≥ 1. We discuss D1 Aβ first. Upon fixing any p ∈ U , the map g 7→ Aβ (g, p) is a bounded linear operator from C k (X) to itself, so for any g ∈ C k (X), D1 Aβ (g, p) = (f 7→ Aβ (f, p)) . To show that D1 Aβ is continuous, it suffices to check that if pn → p in U , then the map f 7→ Aβ (f, pn ) converges to f 7→ Aβ (f, p) in the B(C k (X), C k (X)) norm. The continuity of Aβ gives us that for each fixed f ∈ C k (X), we have Aβ (f, pn ) → Aβ (f, p) in C k (X); we would now like to show that this convergence is uniform across f ∈ C k (X) such that ∥f ∥C k (X) = 1. Using (65) and the definition of Hf,α from (64), we have ∂ α [Aβ (f, pn )](y) − ∂ α [Aβ (f, p)](y) Z = f (Ψ(pn , x))∂yα kβ (Ψ(pn , x), y) − f (Ψ(p, x))∂yα kβ (Ψ(p, x), y)dµ(x) Z = Hf,α (pn , x, y) − Hf,α (p, x, y)dµ(x), for any f ∈ C k (X), n ∈ N and multi-index α such that |α| ≤ k. Let r > 0 be such that Br (p) ⊂ U . By the mean value theorem and the chain rule, for each multi-index α such that |α| ≤ k, there exists C > 0 (which can depend on the derivatives of kβ and the first derivatives of Ψ with respect to the first variable) such that for large enough n, sup |Hf,α (pn , x, y) − Hf,α (p, x, y)| ≤ ∥pn − p∥ x,y∈X
sup q∈Br (p), x,y∈X
42
∥∇q Hf,α (q, x, y)∥ ≤ C ∥pn − p∥ ∥f ∥C 1 (X) .
It follows that there exists C ′ > 0 such that ∥Aβ (f, pn ) − Aβ (f, p)∥C k (X) ≤ C ′ ∥pn − p∥ ∥f ∥C 1 (X) for large enough n, so the map f 7→ Aβ (f, pn ) does indeed converge to f 7→ Aβ (f, p) in the B(C k (X), C k (X)) norm for k ≥ 1. Thus D1 Aβ is continuous. Differentiability of D1 Aβ when k ≥ 2. Now suppose k ≥ 2. We want to show that D1 Aβ is C k−1 . Again we look at its partial derivatives: we can see that D12 Aβ is identically 0 (since D1 Aβ (g, p) does not depend on g), whereas the existence and continuity of D2 D1 Aβ : C k (X) × U → B Rm , B(C k (X), C k (X)) follows from a Taylor expansion. For simplicity, let us consider when m = 1 (i.e., U ⊂ R); in this case, we can identify B(R, B(C k (X), C k (X))) with B(C k (X), C k (X)). The case m > 1 is similar. Fix any p ∈ U ⊂ R, and define the operator Λp : C k (X) → C k (X) by Z Λp (f ) := y 7→ ∂p Hf (p, x, y)dµ(x) (66) for each f ∈ C k (X), where we recall that Hf is given by (63). One can use the dominated convergence theorem to check that we indeed have Λp (f ) ∈ C k (X) for each f ∈ C k (X) and Λp ∈ B(C k (X), C k (X)). Moreover, there exists C > 0 such that ∥Aβ (f, p + ε) − Aβ (f, p) − εΛp (f )∥C k (X) ≤ C∥f ∥C 2 (X) |ε|2 for any f ∈ C k (X) and ε ∈ R with small enough |ε|. Indeed, since Ψ is twice differentiable with respect to the first variable, we can apply Taylor’s theorem with the Lagrange remainder to Hf (·, x, y) and Hf,α (·, x, y) for each x, y ∈ X and multi-index |α| ≤ k. The existence of such a C > 0 then follows from the first two derivatives of Ψ with respect to p being continuous and using that ∂p ∂yα Hf = ∂yα ∂p Hf for each multi-index |α| ≤ k. Thus lim
∥D1 Aβ (g, p + ε) − D1 Aβ (g, p) − εΛp ∥B(C k (X),C k (X)) |ε|
ε→0
=0
for any g ∈ C k (X) and p ∈ U , so D2 D1 Aβ (g, p) exists and is given by Λp , upon identifying B(R, B(C k (X), C k (X))) with B(C k (X), C k (X)). That D2 D1 Aβ is continuous follows from similar arguments as the continuity of D1 Aβ . The arguments for higher-order derivatives are also similar. Existence and regularity of D2 Aβ when k ≥ 1. Again, for simplicity, consider m = 1, so U ⊂ R. Fixing any g ∈ C k (X) and p ∈ U , we have lim
ε→0
∥Aβ (g, p + ε) − Aβ (g, p) − εΛp (g)∥C k (X) |ε|
= lim
ε→0
Aβ (g, p + ε) − Aβ (g, p) − Λp (g) =0 ε C k (X)
by the mean value theorem, the equality ∂p ∂yα Hg = ∂yα ∂p Hg for each multi-index |α| ≤ k and a uniform continuity argument. Therefore D2 Aβ : C k (X) × U → B(R, C k (X)) exists and is given by D2 Aβ (g, p) = (ε 7→ εΛp (g)) 43
for each g ∈ C k (X) and p ∈ U . Upon making the identification B(R, C k (X)) ∼ = C k (X), we see that k D2 Aβ is continuous: if gn → g in C (X) and pn → p in U , then there exists C > 0 (which can depend on ∂p Ψ and the derivatives of kβ , but not n) such that ∥Λpn (gn ) − Λpn (g)∥C k (X) ≤ C∥gn − g∥C 1 (X) for all n ∈ N, and also ∥Λpn (g) − Λp (g)∥C k (X) → 0 as n → ∞ by a uniform continuity argument. Hence ∥Λpn (gn ) − Λp (g)∥C k (X) → 0 as n → ∞, i.e., g, p 7→ Λp (g) is continuous from C k (X) × U to C k (X). The arguments for establishing that D2 Aβ is in fact C k−1 when k ≥ 2 are similar to those already discussed. Proof of Lemma 4.4. It suffices to check that e H : U × U × U × U → C(X) p, q, r, s 7→ fp,q (Ψ(r, ·)), gp,q (Ψ(s, ·)) is a C k map. We first establish the continuity of H. Suppose (pn , qn , rn , sn ) → (p, q, r, s) in U × U × U × U . Then for any n ∈ N, ∥H(pn , qn , rn , sn ) − H(p, q, rn , sn )∥C(X) ≤ ∥(fpn ,qn , gpn ,qn ) − (fp,q , gp,q )∥C(X) , e e and the right-hand side goes to 0 as n → ∞ by Lemma 4.3. Moreover, ∥H(p, q, rn , sn ) − H(p, q, r, s)∥C(X) e = ∥(fp,q (Ψ(rn , ·)), gp,q (Ψ(sn , ·))) − (fp,q (Ψ(r, ·)), gp,q (Ψ(s, ·)))∥C(X) e also goes to 0 as n → ∞, as (upon fixing a representative of (fp,q , gp,q )) r′ , s′ , x 7→ fp,q (Ψ(r′ , x)), gp,q (Ψ(s′ , x)) is uniformly continuous on the compact set BR (r) × BR (s) × X for some small enough R > 0. Thus ∥H(pn , qn , rn , sn ) − H(p, q, r, s)∥C(X) →0 e as n → ∞. We now check that the first partial derivatives of H exist and are continuous. For the partial e q) := (fp,q , gp,q ) is a C k map from U × U derivatives of H with respect to p, q, by Lemma 4.3, H(p, k e e e q) ∈ B(Rm × Rm , C e k (X)). to C (X) (and hence C(X)). For each p, q ∈ U × U , let Λp,q := DH(p, Thus for any r, s ∈ U , ∥H(p + εp , q + εq , r, s) − H(p, q, r, s) − Λp,q (εp , εq ) ◦ (Ψ(r, ·), Ψ(s, ·))∥C(X) e ∥(εp , εq )∥ 44
→0
as (εp , εq ) → 0. Hence for each p, q, r, s ∈ U, the partial derivative Dp,q H(p, q, r, s) ∈ B(Rm × e Rm , C(X)) exists and is given by Dp,q H(p, q, r, s) = (εp , εq , 7→ Λp,q (εp , εq ) ◦ (Ψ(r, ·), Ψ(s, ·))) . e k (X)), the continuity of As p, q 7→ Λp,q is a C k−1 map from U × U to B(Rm × Rm , C 2m ∼ e e Dp,q H : U × U × U × U → B(Rm × Rm , C(X)) = C(X) follows from a similar argument as the continuity of H. For the partial derivative of H with respect to r, let fep,q,x (r) := fp,q (Ψ(r, x)) for each p, q ∈ U and x ∈ X. (By a slight abuse of notation, we have chosen an arbitrary representative for fp,q .) Since one partial derivative of Ψ with respect to the first variable exists and is continuous on U × X, fep,q,x is C 1 on U for each p, q ∈ U and x ∈ X, and r, x 7→ ∇fep,q,x (r) is continuous on U × X for each p, q ∈ U . Fix any p, q, r ∈ U . By the mean value theorem, for any ε ∈ Rm with small enough norm, fep,q,x (r + ε) − fep,q,x (r) = ∇fep,q,x (r + αε) · ε for some α ∈ (0, 1). As r′ , x 7→ ∇fep,q,x (r′ ) is uniformly continuous on BR (r) × X for some small enough R > 0, we have lim
supx∈X fep,q,x (r + ε) − fep,q,x (r) − ∇fep,q,x (r) · ε ∥ε∥
ε→0
= 0.
e Hence for each p, q, r, s ∈ U , Dr H(p, q, r, s) ∈ B(Rm , C(X)) exists and is given by h i [Dr H(p, q, r, s)](ε) = x 7→ ∇fep,q,x (r) · ε, 0 . To show that
m ∼ e e Dr H : U × U × U × U → B(Rm , C(X)) = C(X)
is continuous, observe that Lemma 4.3 implies that p, q 7→ ∇fp,q is a C k map from U × U to d C k−1 (X) . Using the assumptions on the regularity of Ψ, another uniform continuity argument then implies that p, q, r 7→ x 7→ ∇fep,q,x (r) d
is continuous from U × U × U to (C(X)) , and the continuity of Dr H follows. The arguments for Ds H are completely analogous. If k ≥ 2, the arguments for higher-order derivatives are similar. Indeed, for each p, q, r, s ∈ U , Dp,q H(p, q, r, s) is given by pre-composing Λp,q with (Ψ(r, ·), Ψ(s, ·)), where by a slight abuse of e k (X))2m ; now recall that p, q 7→ Λp,q is a C k−1 map notation we consider Λp,q as an element of (C m m ek k e (X))2m by Lemma 4.3. Similarly, for Dr H, if k ≥ 2 one from U × U to B(R × R , C (X)) ∼ = (C can prove the existence and continuity of additional derivatives with respect to r (again via the chain rule, assumptions on the regularity of Ψ and uniform continuity arguments) and with respect to p, q using continuity of the partial derivatives (29) and Lemma 4.3. 45
C
Additional details and plots for Section 5
We first prove Claim 5.1. Proof of Claim 5.1. Fix any θ ∈ [−π, π). By symmetry, we have that the optimal coupling is of the form 1 − p(θ) p(θ) ∗ 2 Πθ = , (67) 1 p(θ) 2 − p(θ) and it remains to prove that p(θ) is given by (41). Since µ0 and µθ are both discrete measures with finite support, we can write (12) as # " 2 a1 b1 a1 b2 e−R /β 1 ∗ 2 Πθ = (68) θ 2 4 a2 b1 e−R /β a2 b2 e−(2R sin( 2 )) /β where a1 = ef1 /β , a2 = ef2 /β , b1 = eg1 /β , b2 = eg2 /β for some f1 , f2 , g1 , g2 ∈ R. Without loss of generality, we can take a1 = b1 , and since 2
2
a1 b2 e−R /β = a2 b1 e−R /β we also have a2 = b2 . Therefore 1 − p(θ) Π∗θ = 2 p(θ)
" a21 1 p(θ) = 2 1 4 a1 a2 e−R /β 2 − p(θ)
2
a1 a2 e−R /β 2 θ a2 e−(2R sin( 2 )) /β
#
2
for some a1 , a2 > 0, so 2 2 1 θ 2 2 θ 2 a21 a22 e−(2R sin( 2 )) /β 2 − p(θ) = e2R (1−2 sin ( 2 ))/β = e2R cos(θ)/β . = 2 2 p(θ) a1 a2 e−R /β Taking the square root of both sides gives us 2 1 − 1 = eR cos(θ)/β , 2p(θ)
which, upon rearranging, gives us (41): 1 . 2 1 + eR2 cos(θ)/β
p(θ) = (42) is obtained by observing that e
f1 +g1 β
(
=e
2
( θ2 ))
f2 +g2 − 2R sin β
=4
1 − p(θ) 2
= 2 − 4p(θ)
and the dual representation (11) of OTβ gives us 1 (f1 + f2 + g1 + g2 ) 2 2 ! 1 θ = β log(2 − 4p(θ)) + β log(2 − 4p(θ)) + 2R sin 2 2
OTβ (µ0 , µθ ) =
= β log(2 − 4p(θ)) + R2 (1 − cos(θ)) . 46
(β)
The expression for Sβ (µ0 , µθ ) follows from the definition of Sβ . Now g0 can be computed directly by considering the second derivative of H(θ) := Sβ (µ0 , µθ ) evaluated at θ = 0. Indeed, 2 R2 sin(θ) 1 R cos(θ)/β p′ (θ) = − e − 2 β 2 1 + eR2 cos(θ)/β ! 2 eR cos(θ)/β R2 sin(θ) = p(θ) 2 β 1 + eR cos(θ)/β 2 R sin(θ) = p(θ)(1 − 2p(θ)) β for θ ∈ (−π, π), and so βp′ (θ) H ′ (θ) = − 1 + R2 sin(θ) − p(θ) 2 = −2p(θ) R2 sin(θ) + R2 sin(θ) = R2 sin(θ) (1 − 2p(θ)) for θ ∈ (−π, π). Taking one more derivative gives us d d 1 1 (β) g0 , = H ′′ (0) = R2 (1 − 2p(0)). dθ dθ 2 2 (β)
By rotational symmetry, we have gθ
(β)
= g0
for all θ ∈ [−π, π).
As noted in Remark 5.1, in practice we often do not have a closed-form formula for the Sinkhorn divergence between two probability measures, and instead we have to approximate it. For Figures 1 and 4 (demonstrating the approximation error on the squared geodesic distance arising from Lemma 3.1), we find that the difference between using the closed-form formula and the Sinkhorn divergence computed using the default implementation of Sinkhorn’s algorithm in the Python Optimal Transport (POT) library [35, 36] is very small; see Figures 7 and 8, respectively. One therefore expects the difference between the resulting discrete Laplacians to also be very small. As an example, in Figure 9, we recreate Figure 2 but with the Sinkhorn divergence computed using the default implementation of Sinkhorn’s algorithm in POT and different values of β (0.2, 0.4, 0.6, 0.8, 1). The difference between using the closed-form formula and the POT version is virtually unnoticeable unless we specifically plot the difference between the resulting discrete Laplacians (Figure 10). Observe that for each β ∈ {0.2, 0.4, 0.6, 0.8, 1}, the scale on the y-axis in Figure 10 is much smaller than that on the corresponding plot in Figure 9.
47
b (µ ,µ )−S (µ ,µ ) S Figure 7: Plots of θ 7→ β 0 θ θ4 β 0 θ , where Sbβ is the Sinkhorn divergence approximated with Sinkhorn’s algorithm and Sβ (µ0 , µθ ) is given by the closed-form formula for the example in Section 5.1.
48
b (µ ,µx )−Sβ (µ1 ,µx ) S Figure 8: Plots of x 7→ β 1 (x−1) , where Sbβ is the Sinkhorn divergence approximated 4 with Sinkhorn’s algorithm and Sβ (µ1 , µx ) is given by the closed-form formula for the example in Section 5.2. 49
Figure 9: Same as Fig. 2, but with different values of β (0.2, 0.4, 0.6, 0.8, 1) and using the Sinkhorn divergence Sbβ approximated via Sinkhorn’s algorithm.
50
Figure 10: Subtracting the right-hand column from the left-hand column of Fig. 9. 51