Mathieu Carriere

Marco Cuturi

Steve Y Oudot

Steve Y Oudot Sliced Wasserstein

Mathieu Carrière

Sliced Wasserstein Kernel for Persistence Diagrams

come   L'archive ouverte pluridisciplinaire

Introduction

Topological Data Analysis (TDA) is an emerging trend in data science, grounded on topological methods to design descriptors for complex data-see e.g. [START_REF] Carlsson | Topology and data[END_REF] for an introduction to the subject. The descriptors of TDA can be used in various contexts, in particular statistical learning and geometric inference, where they provide useful insight into the structure of data. Applications of TDA can be found in a number of scientific areas, including computer vision [START_REF] Li | Persistence-Based Structural Recognition[END_REF], materials science [START_REF] Hiraoka | Hierarchical structures of amorphous solids characterized by persistent homology[END_REF], and brain science [START_REF] Singh | Topological analysis of population activity in visual cortex[END_REF], to name Proceedings of the 34 th International Conference on Machine Learning, Sydney, Australia, PMLR 70, 2017. Copyright 2017 by the author(s). a few. The tools developed in TDA are built upon persistent homology theory [START_REF] Edelsbrunner | Computational Topology: an introduction[END_REF][START_REF] Oudot | Persistence Theory: From Quiver Representations to Data Analysis[END_REF], and their main output is a descriptor called persistence diagram (PD), which encodes the topology of a space at all scales in the form of a point cloud with multiplicities in the plane R2 -see Section 2.1 for more details.

PDs as features. The main strength of PDs is their stability with respect to perturbations of the data (Chazal et al., 2009b;[START_REF] Chazal | Persistence stability for geometric complexes[END_REF]. On the downside, their use in learning tasks is not straightforward. Indeed, a large class of learning methods, such as SVM or PCA, requires a Hilbert structure on the descriptors space, which is not the case for the space of PDs. Actually, many simple operators of R n , such as addition, average or scalar product, have no analogues in that space. Mapping PDs to vectors in R n or in some infinite-dimensional Hilbert space is one possible approach to facilitate their use in discriminative settings.

Related work.

A series of recent contributions have proposed kernels for PDs, falling into two classes. The first class of methods builds explicit feature maps: One can, for instance, compute and sample functions extracted from PDs [START_REF] Bubenik | Statistical Topological Data Analysis using Persistence Landscapes[END_REF][START_REF] Adams | Persistence Images: A Stable Vector Representation of Persistent Homology[END_REF][START_REF] Robins | Principal Component Analysis of Persistent Homology Rank Functions with case studies of Spatial Point Patterns, Sphere Packing and Colloids[END_REF]; sort the entries of the distance matrices of the PDs [START_REF] Carrière | Stable Topological Signatures for Points on 3D Shapes[END_REF]; treat the PD points as roots of a complex polynomial, whose coefficients are concatenated [START_REF] Fabio | Comparing persistence diagrams through complex vectors[END_REF]. The second class of methods, which is more relevant to our work, defines implicitly feature maps by focusing instead on building kernels for PDs. For instance, [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF] use solutions of the heat differential equation in the plane and compare them with the usual L 2 (R 2 ) dot product. [START_REF] Kusano | Persistence Weighted Gaussian Kernel for Topological Data Analysis[END_REF] handle a PD as a discrete measure on the plane, and follow by using kernel mean embeddings with Gaussian kernels-see Section 4 for precise definitions. Both kernels are provably stable, in the sense that the metric they induce in their respective reproducing kernel Hilbert space (RKHS) is bounded above by the distance between PDs. Although these kernels are injective, there is no evidence that their induced RKHS distances are discriminative and therefore follow the geometry of the bottleneck distances, which are more widely accepted distances to compare PDs.

Contributions.

In this article, we use the sliced Wasserstein (SW) distance [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] to define a new ker-nel for PDs, which we prove to be both stable and discriminative. Specifically, we provide distortion bounds on the SW distance that quantify its ability to mimic bottleneck distances between PDs. This is in contrast to other kernels for PDs, which only focus on stability. We also propose a simple approximation algorithm to speed up the computation of that kernel, confirm experimentally its discriminative power and show that it outperforms experimentally both proposals of [START_REF] Kusano | Persistence Weighted Gaussian Kernel for Topological Data Analysis[END_REF] and [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF] in several supervised classification problems.

Background on TDA and Kernels

We briefly review in this section relevant material on TDA, notably persistence diagrams, and technical properties of positive and negative definite kernel functions.

Persistent Homology

Persistent homology [START_REF] Zomorodian | Computing persistent homology[END_REF][START_REF] Edelsbrunner | Persistent homology-a survey[END_REF][START_REF] Oudot | Persistence Theory: From Quiver Representations to Data Analysis[END_REF] is a technique inherited from algebraic topology for computing stable signatures on real-valued functions. Given f : X → R as input, persistent homology outputs a planar point set with multiplicities, called the persistence diagram of f and denoted by Dg f . See Figure 1 for an example. To understand the meaning of each point in this diagram, it suffices to know that, to compute Dg f , persistent homology considers the family of sublevel sets of f , i.e. the sets of the form f -1 ((-∞, t]) for t ∈ R, and it records the topological events (e.g. creation or merge of a connected component, creation or filling of a loop, void, etc.) that occur in f -1 ((-∞, t]) as t ranges from -∞ to +∞. Then, each point p ∈ Dg f represents the lifespan of a particular topological feature (connected component, loop, void, etc.), with its creation and destruction times as coordinates. See again Figure 1 for an illustration.

For the interested reader, we point out that the mathematical tool used by persistent homology to track the topological events in the family of sublevel sets is homological algebra, which turns the parametrized family of sublevel sets into a parametrized family of vector spaces and linear maps. Computing persistent homology then boils down to computing a family of bases for the vector spaces, which are compatible with the linear maps.

Distance between PDs. We now define the pth diagram distance between PDs. Let p ∈ N and Dg 1 , Dg 2 be two PDs. Let Γ : Dg 1 ⊇ A → B ⊆ Dg 2 be a partial bijection between Dg 1 and Dg 2 . Then, for any point x ∈ A, the cost of x is defined as c(x) := x -Γ(x) p ∞ , and for any point y ∈ (Dg 1 Dg 2 ) \ (A B), the cost of y is defined as c (y) := yπ ∆ (y) p ∞ , where π ∆ is the projection onto the diagonal ∆ = {(x, x) | x ∈ R}. The cost c(Γ) is defined as: c(Γ) := ( x c(x) + y c (y)) 1/p . We then define the pth diagram distance d p as the cost of the best partial bijection between the PDs:

d p (Dg 1 , Dg 2 ) = inf Γ c(Γ).
In the particular case p = +∞, the cost of Γ is defined as c(Γ) := max{max x δ(x) + max y δ (y)}. The corresponding distance d ∞ is often called the bottleneck distance. One can show that d p → d ∞ when p → +∞. A fundamental property of PDs is their stability with respect to (small) perturbations of their originating functions. Indeed, the stability theorem [START_REF] Bauer | Induced matchings and the algebraic stability of persistence barcodes[END_REF]Chazal et al., 2009a;2016;[START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] asserts that for any f, g : X → R, we have

d ∞ (Dg f, Dg g) ≤ f -g ∞ , (1)
See again Figure 1 for an illustration.

In practice, PDs can be used as descriptors for data via the choice of appropriate filtering functions f , e.g. distance to the data in the ambient space, eccentricity, curvature, etc.

The main strengths of the obtained descriptors are: (a) to be provably stable as mentioned previously; (b) to be invariant under reparametrization of the data; and (c) to encode information about the topology of the data, which is complementary and of an essentially different nature compared to geometric or statistical quantities. These properties have made persistence diagrams useful in a variety of contexts, including the ones mentioned in the introduction of the paper. For further details on persistent homology and on applications of PDs, the interested reader can refer e.g. to [START_REF] Oudot | Persistence Theory: From Quiver Representations to Data Analysis[END_REF] and the references therein.

Kernel Methods

Positive Definite Kernels. Given a set X, a function k : X × X → R is called a positive definite kernel if for all integers n, for all families x 1 , ..., x n of points in X, the matrix [k(x i , x j )] i,j is itself positive semi-definite.

For brevity we will refer to positive definite kernels as kernels in the rest of the paper. It is known that kernels generalize scalar products, in the sense that, given a kernel k, there exists a Reproducing Kernel Hilbert Space (RKHS) H k and a feature map φ :

X → H k such that k(x 1 , x 2 ) = φ(x 1 ), φ(x 2 ) H k .
A kernel k also induces a distance d k on X that can be computed as the Hilbert norm of the difference between two embeddings:

d 2 k (x 1 , x 2 ) def. = k(x 1 , x 1 ) + k(x 2 , x 2 ) -2 k(x 1 , x 2 ).
We will be particularly interested in this distance, since one of the goals we will aim for will be that of designing a kernel k for persistence diagrams such that d k has low distortion with respect to d 1 . Negative Definite and RBF Kernels. A standard way to construct a kernel is to exponentiate the negative of a Euclidean distance. Indeed, the Gaussian kernel for vectors with parameter σ > 0 does follow that template approach:

R X p q s (a) R X p q s (b) -∞ +∞ +∞ f (p) f (s) (c)
k σ (x, y) = exp -x-y 2 2σ 2
. An important theorem of [START_REF] Berg | Harmonic Analysis on Semigroups: Theory of Positive Definite and Related Functions[END_REF] (Theorem 3.2.2,p.74) states that such an approach to build kernels, namely setting

k σ (x, y) def. = exp - f (x, y) 2σ 2 ,
for an arbitrary function f can only yield a valid positive definite kernel for all σ > 0 if and only if f is a negative semi-definite function, namely that, for all integers n, ∀x 1 , ..., x n ∈ X, ∀a 1 , ..., a n ∈ R n such that i a i = 0, i,j a i a j f (x i , x j ) ≤ 0. Unfortunately, as observed in Appendix A of [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF], d 1 is not negative semi-definite (it only suffices to sample a family of point clouds to observe experimentally that more often than not the inequality above will be violated for a particular weight vector a). In this article, we use an approximation of d 1 with the Sliced Wasserstein distance, which is provably negative semi-definite, and we use it to define a RBF kernel that can be easily tuned thanks to its bandwidth parameter σ.

Wasserstein distance for unnormalized measures on R

The Wasserstein distance (Villani, 2009, §6) is a distance between probability measures. For reasons that will become clear in the next section, we will focus on a variant of that distance: the 1-Wasserstein distance for nonnegative, not necessarily normalized, measures on the real line (Santambrogio, 2015, §2). Let µ and ν be two nonnegative mea-sures on the real line such that |µ| = µ(R) and |ν| = ν(R) are equal to the same number r. We define the three following objects:

W(µ, ν) = inf P ∈Π(µ,ν) R×R |x -y|P (dx, dy) (2) Q r (µ, ν) = r R |M -1 (x) -N -1 (x)|dx (3) L(µ, ν) = inf f ∈1-Lipschitz R f (x)[µ(dx) -ν(dx)] (4)
where Π(µ, ν) is the set of measures on R 2 with marginals µ and ν, and M -1 and N -1 the generalized quantile functions of the probability measures µ/r and ν/r respectively.

Proposition 2.1. We have W = Q r = L. Additionally (i) Q r is negative definite on the space of measures of mass r;

(ii) for any three positive measures µ, ν, γ such that |µ| = |ν|, we have

L(µ + γ, ν + γ) = L(µ, ν).
Equation ( 2) is the generic Kantorovich formulation of optimal transport, which is easily generalized to other cost functions and spaces, the variant being that we consider an unnormalized mass by reflecting it directly in the set Π.

The equality between (2) and ( 3) is only valid for probability measures on the real line. Because the cost function | • | is homogeneous, we see that the scaling factor r can be removed when considering the quantile function and multiplied back. The equality between (2) and ( 4) is due to the well known Kantorovich duality for a distance cost (Villani, 2009, Particular case 5.4) which can also be trivially generalized to unnormalized measures, proving therefore the main statement of the proposition. The definition of Q r shows that the Wasserstein distance is the l 1 norm of rM -1 -rN -1 , and is therefore a negative definite kernel (as the l 1 distance between two direct representations of µ and ν as functions rM -1 and rN -1 ), proving point (i).

The second statement is immediate.

We conclude with an important practical remark: for two unnormalized uniform empirical measures µ = n i=1 δ xi and ν = n i=1 δ yi of the same size, with ordered

x 1 ≤ • • • ≤ x n and y 1 ≤ • • • ≤ y n , one has: L(µ, ν) = n i=1 |x i -y i | = X -Y 1 , where X = (x 1 , ..., x n ) ∈ R n and Y = (y 1 , ..., y n ) ∈ R n .

The Sliced Wasserstein Kernel

In this section we define a new kernel between PDs, called the Sliced Wasserstein (SW) kernel, based on the Sliced Wasserstein metric of [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF]. The idea underlying this metric is to slice the plane with lines passing through the origin, to project the measures onto these lines where W is computed, and to integrate those distances over all possible lines. Formally: Definition 3.1. Given θ ∈ R 2 with θ 2 = 1, let L(θ) denote the line {λ θ | λ ∈ R}, and let π θ : R 2 → L(θ) be the orthogonal projection onto L(θ). Let Dg 1 , Dg 2 be two PDs, and let µ θ 1 := p∈Dg 1 δ π θ (p) and µ θ 1∆ := p∈Dg 1 δ π θ •π∆(p) , and similarly for µ θ 2 , where π ∆ is the orthogonal projection onto the diagonal. Then, the Sliced Wasserstein distance is defined as:

SW(Dg 1 , Dg 2 ) def. = 1 2π S1 W(µ θ 1 + µ θ 2∆ , µ θ 2 + µ θ 1∆ )dθ.
Note that, by symmetry, one can restrict on the half-circle [-π 2 , π 2 ] and normalize by π instead of 2π. Since Q r is negative semi-definite, we can deduce that SW itself is negative semi-definite: Lemma 3.2. Let X be the set of bounded and finite PDs. Then, SW is negative semi-definite on X.

Proof. Let n ∈ N * , a 1 , ..., a n ∈ R such that i a i = 0 and Dg 1 , ..., Dg n ∈ X. Given 1 ≤ i ≤ n, we let μθ i := µ θ i + q∈Dg k ,k =i δ π θ •π∆(q) , μθ ij∆ := p∈Dg k ,k =i,j δ π θ •π∆(p) and d = i |Dg i |. Then: i,j a i a j W(µ θ i + µ θ j∆ , µ θ j + µ θ i∆ ) = i,j a i a j L(µ θ i + µ θ j∆ , µ θ j + µ θ i∆ ) = i,j a i a j L(µ θ i + µ θ j∆ + µ θ ij∆ , µ θ j + µ θ i∆ + µ θ ij∆ ) = i,j a i a j L(μ θ i , μθ j ) = i,j a i a j Q d (μ θ i , μθ j ) ≤ 0
The result follows by linearity of integration.

Hence, the theorem of [START_REF] Berg | Harmonic Analysis on Semigroups: Theory of Positive Definite and Related Functions[END_REF] allows us to define a valid kernel with:

k SW (Dg 1 , Dg 2 ) def. = exp - SW(Dg 1 , Dg 2 ) 2σ 2 . ( 5
)
Metric equivalence. We now give the main theoretical result of this article, which states that SW is equivalent to d 1 . This has to be compared with [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF] and [START_REF] Kusano | Persistence Weighted Gaussian Kernel for Topological Data Analysis[END_REF], which only prove stability and injectivity. Our equivalence result states that the k SW , in addition to be stable and injective, also preserves the metric between PDs, which should intuitively lead to an improvement of the classification power. This intuition is illustrated in Section 4 and Figure 4, where we show an improvement of classification accuracies on several benchmark applications. Theorem 3.3. Let X be the set of bounded PDs with cardinalities bounded by N ∈ N * . Let Dg 1 , Dg 2 ∈ X. Then, one has:

d 1 (Dg 1 , Dg 2 ) 2M ≤ SW(Dg 1 , Dg 2 ) ≤ 2 √ 2d 1 (Dg 1 , Dg 2 ),
where

M = 1 + 2N (2N -1). Proof. Let s θ : Dg 1 ∪ π ∆ (Dg 2 ) → Dg 2 ∪ π ∆ (Dg 1 ) be the one-to-one bijection between Dg 1 ∪ π ∆ (Dg 2 ) and Dg 2 ∪ π ∆ (Dg 1 ) induced by W(µ θ 1 + µ θ 2∆ , µ θ 2 + µ θ 1∆
), and let s be the one-to-one bijection between Dg 1 ∪ π ∆ (Dg 2 ) and Dg 2 ∪ π ∆ (Dg 1 ) induced by the partial bijection achieving d 1 (Dg 1 , Dg 2 ).

Upper bound. Recall that θ 2 = 1. We have:

W(µ θ 1 + µ θ 2∆ , µ θ 2 + µ θ 1∆ ) = | p -s θ (p), θ | ≤ | p -s(p), θ | ≤ √ 2 p -s(p) ∞ ≤ 2 √ 2d 1 (Dg 1 , Dg 2 ),
where the sum is taken over all p ∈ Dg 1 ∪ π ∆ (Dg 2 ). The upper bound follows by linearity.

Lower bound. The idea is to use the fact that s θ is a piecewise-constant function of θ, and that it has at most

2 + 2N (2N -1) critical values Θ 0 , ..., Θ M in [-π 2 , π 2 ].
Indeed, it suffices to look at all θ such that p 1 -p 2 , θ = 0 for some p 1 , p 2 in Dg 1 ∪ π ∆ (Dg 2 ) or Dg 2 ∪ π ∆ (Dg 1 ). Then:

Θi+1 Θi | p -s θ (p), θ |dθ = p -s Θi (p) 2 Θi+1 Θi |cos(∠(p -s Θi (p), θ))|dθ ≥ p -s Θi (p) 2 (Θ i+1 -Θ i ) 2 /2π ≥ (Θ i+1 -Θ i ) 2 d 1 (Dg 1 , Dg 2 )/2π,
where the sum is again taken over all p ∈ Dg 1 ∪ π ∆ (Dg 2 ), and where the inequality used to lower bound the integral of the cosine is obtained by concavity. The lower bound follows then from the Cauchy-Schwarz inequality.

Note that the lower bound depends on the cardinalities of the PDs, and it becomes close to 0 if the PDs have a large number of points. On the other hand, the upper bound is oblivious to the cardinality. A corollary of Theorem 3.3 is that d kSW , the distance induced by k SW in its RKHS, is also equivalent to d 1 in a broader sense: there exist continuous, positive and monotone functions g, h such that

g(0) = h(0) = 0 and g • d 1 ≤ d kSW ≤ h • d 1 .
When the condition on the cardinalities of PDs is relaxed, e.g. when we only assume the PDs to be finite and bounded, with no uniform bound, the feature map φ SW associated to k SW remains continuous and injective w.r.t. d 1 .

This means that k SW can be turned into a universal kernel by considering exp(k SW ) (cf Theorem 1 in [START_REF] Kwitt | Statistical Topological Data Analysis -A Kernel Perspective[END_REF]). This can be useful in a variety of tasks, including tests on distributions of PDs.

Computation. In practice, we propose to approximate k SW in O(N log(N )) time using Algorithm 1. This algorithm first samples M directions in the half-circle S + 1 ; it then computes, for each sample θ i and for each PD Dg, the scalar products between the points of Dg and θ i , to sort them next in a vector V θi (Dg). Finally, the 1 -norm between the vectors is averaged over the sampled directions:

SW M (Dg 1 , Dg 2 ) = 1 M M i=1 V θi (Dg 1 ) -V θi (Dg 2 ) 1 .
Note that one can easily adapt the proof of Lemma 3.2 to show that SW M is negative semi-definite by using the linearity of the sum. Hence, this approximation remains a kernel. If the two PDs have cardinalities bounded by N , then the running time of this procedure is O(M N log(N )). This approximation of k SW is useful since, as shown in Section 4, we have observed empirically that just a few directions are sufficient to get good classification accuracies. Note that the exact computation of k SW is also possible in O(N 2 log(N )) time using the algorithm described in [START_REF] Carrière | Sliced Wasserstein Kernel for Persistence Diagrams[END_REF].

Experiments

In this section, we compare k SW to k PSS and k PWG on several benchmark applications for which PDs have been proven useful. We compare these kernels in terms of classification accuracies and computational cost. We review first our experimental setting, and then all our tasks.

Experimental setting All kernels are handled with the LIBSVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] implementation of C-SVM, and results are averaged over 10 runs on a 2.4GHz Intel Xeon E5530 Quad Core.

The

Algorithm 1 Computation of SW M Input: Dg 1 = {p 1 1 ... p 1 N1 }, Dg 2 = {p 2 1 ... p 2 N2 }, M . Add π ∆ (Dg 1 )
to Dg 2 and vice-versa. Let SW M = 0; θ = -π/2; s = π/M ; for i = 1 ... M do Store the products p 1 k , θ in an array V 1 ; Store the products p 2 k , θ in an array V 2 ; Sort V 1 and V 2 in ascending order; cost factor C is cross-validated in the following grid: {0.001, 0.01, 0.1, 1, 10, 100, 1000}. Table 1 summarizes the number of labels, and the number of training and test instances for each task. Figure 2 illustrate how we use PDs to represent complex data. We first describe the two baselines we considered, along with their parameterization, followed by our proposal.

SW M = SW M + s V 1 -V 2 1 ; θ = θ
PSS. The Persistence Scale Space kernel k PSS [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF] is defined as the scalar product of the two solutions of the heat diffusion equation with initial Dirac sources located at the PD points. It has the following closed form expression:

k PSS (Dg 1 , Dg 2 ) = 1 8πt p∈Dg 1 q∈Dg 2 exp -p-q 2 8t -exp -p-q 2 8t
, where q = (y, x) is the symmetric of q = (x, y) along the diagonal. Since there is no clear heuristic on how to tune t, this parameter is chosen in the applications by ten-fold cross-validation with random 50%-50% training-test splits and with the following set of N PSS = 13 values: 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 500 and 1000.

PWG. Let K, p, ρ be positive parameters. Let k ρ be the Gaussian kernel with parameter ρ and associated RKHS H ρ . Let Dg 1 , Dg 2 be two PDs, and let

µ 1 := x∈Dg 1 arctan(Kd ∞ (x, ∆) p )k ρ (•, x
) ∈ H ρ be the kernel mean embedding of Dg 1 weigthed by the diagonal distances. Let µ 2 be defined similarly. Let τ > 0. The Persistence Weighted Gaussian kernel k PWG [START_REF] Kusano | Persistence Weighted Gaussian Kernel for Topological Data Analysis[END_REF]2017) is defined as

k PWG (Dg 1 , Dg 2 ) = exp - µ1-µ2 Hρ 2τ 2
, i.e. the Gaussian kernel with parameter τ on H ρ . The authors in [START_REF] Kusano | Persistence Weighted Gaussian Kernel for Topological Data Analysis[END_REF] provide heuristics to compute K, ρ and τ and give a rule of thumb to tune p. Hence, in the applications we select p according to the rule of thumb, and we use ten-fold cross-validation with random 50%-50% training-test splits to chose K, ρ and τ . The ranges of possible values is obtained by multiplying the values computed with the heuristics with the following range of 5 factors: 0.01, 0.1, 1, 10 and 100, leading to N PWG = 5 × 5 × 5 = 125 different sets of parameters.

Parameters for k SW . The kernel we propose has only one parameter, the bandwidth σ in Eq. ( 5), which we choose using ten-fold cross-validation with random 50%-50% training-test splits. The range of possible values is obtained by computing the squareroot of the median, the first and the last deciles of all SW(Dg i , Dg j ) in the training set, then by multiplying these values by the following range of 5 factors: 0.01, 0.1, 1, 10 and 100, leading to N SW = 5 × 3 = 15 possible values.

Parameter Tuning. The bandwidth of k SW is, in practice, easier to tune than the parameters of its two competitors when using grid search. Indeed, as is the case for all infinitely divisible kernels, the Gram matrix does not need to be recomputed for each choice of σ, since it only suffices to compute all the Sliced Wasserstein distances between PDs in the training set once. On the contrary, neither k PSS nor k PWG share this property, and require recomputations for each hyperparameter choice. Note however that this improvement may no longer hold if one uses other methods to tune parameters. For instance, using k PWG without cross-validation is possible with the heuristics given by the authors in [START_REF] Kusano | Persistence Weighted Gaussian Kernel for Topological Data Analysis[END_REF], and leads to smaller training times, but also to worse accuracies.

3D shape segmentation

Our first task, whose goal is to produce point classifiers for 3D shapes, follows that presented in [START_REF] Carrière | Stable Topological Signatures for Points on 3D Shapes[END_REF].

Data. We use some categories of the mesh segmentation benchmark of Chen et al. [START_REF] Chen | A Benchmark for 3D Mesh Segmentation[END_REF], which contains 3D shapes classified in several categories ("airplane", "human", "ant"...). For each category, our goal is to design a classifier that can assign, to each point in the shape, a label that describes the relative location of that point in the shape. For instance, possible labels are, for the human category, "head", "torso", "arm"... To train classifiers, we compute a PD per point using the geodesic distance function to this point-see [START_REF] Carrière | Stable Topological Signatures for Points on 3D Shapes[END_REF] for details. We use 1-dimensional persistent homology (0-dimensional would not be informative since the shapes are connected, leading to solely one point with coordinates (0, +∞) per PD). For each category, the training set contains one hundredth of the points of the first five 3D shapes, and the test set contains one hundredth of the points of the remaining shapes in that category. Points in training and test sets are evenly sampled. See Figure 2. Here, we focus on comparison between PDs, and not on achieving state-of-the-art results. It has been proven that PDs bring complementary information to classical descriptors in this task-see [START_REF] Carrière | Stable Topological Signatures for Points on 3D Shapes[END_REF], hence reinforcing their discriminative power with appropriate kernels is of great interest. Finally, since data points are in R 3 , we set the p parameter of k PWG to 5.

Results. Classification accuracies are given in Table 2. For most categories, k SW outperforms competing kernels by a significant margin. The variance of the results over the run is also less than that of its competitors. However, training times are not better in general. Hence, we also provide the results for an approximation of k SW with 10 directions. As one can see from Table 2 and from Figure 3, this approximation leaves the accuracies almost unchanged, while the training times become comparable with the ones of the other competitors. Moreover, according to Figure 3, using even less directions would slightly decrease the accuracies, but still outperform the competitors performances, while decreasing even more the training times.

Orbit recognition

In our second experiment, we use synthetized data. The goal is to retrieve parameters of dynamical system orbits, following an experiment proposed in [START_REF] Adams | Persistence Images: A Stable Vector Representation of Persistent Homology[END_REF].

Data. We study the linked twist map, a discrete dynamical system modeling fluid flow. It was used in [START_REF] Hertzsch | DNA microarrays: design principles for maximizing ergodic, chaotic mixing[END_REF] to model flows in DNA microarrays. Its orbits can be computed given a parameter r > 0 and initial positions (x 0 , y 0 ) ∈ [0, 1] × [0, 1] as follows:

x n+1 = x n + ry n (1 -y n ) mod 1 y n+1 = y n + rx n+1 (1 -x n+1 ) mod 1
Depending on the values of r, the orbits may exhibit very different behaviors. For instance, as one can see in Figure 2, when r is 2, there seems to be no interesting topological features in the orbit, while voids form when r is 1. Following [START_REF] Adams | Persistence Images: A Stable Vector Representation of Persistent Homology[END_REF], we use 5 different parameters r = 2.5, 3.5, 4, 4.1, 4.3, that act as labels. For each parameter, we generate 100 orbits with 1000 points and random initial positions. We then compute the PDs of the distance functions to the point clouds with the GUDHI library (The GUDHI Project, 2015) and we use them (in all homological dimensions) to produce an orbit classifier that predicts the parameter values, by training over a 70%-30% training-test split of the data. Since data points are in R 2 , we set the p parameter of k PWG to 4.

Results. Since the PDs contain thousands of points, we use kernel approximations to speed up the computation of the Gram matrices. In order for the approximation error to be bounded by 10 -3 , we use an approximation of k SW with 6 directions (as one can see from Figure 3, this has a small impact on the accuracy), we approximate k PWG with 1000 random Fourier features [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], and we approximate k PSS using Fast Gauss Transform [START_REF] Morariu | Automatic online tuning for fast Gaussian summation[END_REF] with a normalized error of 10 -10 . One can see from Table 2 that the accuracy is increased a lot with k SW . Concerning training times, there is also a large improvement since we tune the parameters with grid search. Indeed, each Gram matrix needs not be recomputed for each parameter when using k SW .

Texture classification

Our last experiment is inspired from [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF] and [START_REF] Li | Persistence-Based Structural Recognition[END_REF]. We use the OUTEX00000 data base [START_REF] Ojala | Outex -new framework for empirical evaluation of texture analysis algorithms[END_REF] for texture classification.

Data. PDs are obtained for each texture image by computing first the sign component of CLBP descriptors [START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF] with radius R = 1 and P = 8 neighbors for each image, and then compute the persistent homology of this descriptor using the GUDHI library (The GUDHI Project, 2015). See Figure 2. Note that, contrary to the experiment of [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF], we do not downsample the images to 32 × 32 images, but keep the original 128 × 128 images. Following [START_REF] Reininghaus | A Stable Multi-Scale Kernel for Topological Machine Learning[END_REF], we restrict the focus to 0-dimensional persistent homology. We also use the first 50%-50% training-test split given in the database to produce classifiers. Since data points are in R 2 , we set the p parameter of k PWG to 4.

Results

We use the same approximation procedure as in Section 4.2. According to Figure 3, even though the approximation of SW is rough, this has again a small impact on the accuracy, while reducing the training time by a significant margin. As one can see from Table 2, using k PSS leads to almost state-of-the-art results [START_REF] Ojala | Outex -new framework for empirical evaluation of texture analysis algorithms[END_REF][START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF], closely followed by the accuracies of k SW and k PWG . The best timing is given by k SW , again because we use grid search. Hence, k SW almost achieves the best result, and its training time is better than the ones of its competitors, due to the grid search parameter tuning.

Metric Distortion. To illustrate the equivalence theorem, we also show in Figure 4 a scatter plot where each point represents the comparison of two PDs taken from the Airplane segmentation data set. Similar plots can be obtained with the other datasets considered here. For all points, the x-axis quantifies the first diagram distance d 1 for that pair, while the y-axis is the logarithm of the RKHS distance induced by either k SW , k PSS , k PWG or a Gaussian kernel directly applied to d 1 , to obtain comparable quantities. We use the parameters given by the cross-validation procedure described above. One can see that the distances induced by k SW are less spread than the others, suggesting that the metric induced by k SW is more discriminative. Moreover the distances given by k SW and the Gaussian kernel on d 1 exhibit the same behavior, suggesting that k SW is the best natural equivalent of a Gaussian kernel for PDs.

Conclusion

In this article, we introduce the Sliced Wasserstein kernel, a new kernel for PDs that is provably equivalent to the first diagram distance between PDs. We provide fast algorithms to approximate it, and show on several datasets substantial improvements in accuracy and training times (when tuning parameters is done with grid search) over competing kernels. A particularly appealing property of that kernel is that it is infinitely divisible, substantially facilitating the tuning of parameters through cross validation.

Figure 1 .

 1 Figure 1. Sketch of persistent homology: (a) the horizontal lines are the boundaries of sublevel sets f ((-∞, t]), which are colored in decreasing shades of grey. The vertical dotted lines are the boundaries of their different connected components. For instance, a new connected component is created in the sublevel set f -1 ((-∞, t]) when t = f (p), and it is merged (destroyed) when t = f (s); its lifespan is represented by a copy of the point with coordinates (f (p), f (s)) in the persistence diagram of f (Figure (c)); (b) a piecewiselinear approximation g (blue) of the function f (red) from sampled values; (c) superposition of Dg f (red) and Dg g (blue), showing the partial matching of minimum cost (magenta) between the two persistence diagrams.

Figure 2 .

 2 Figure 2. Examples of PDs computed on orbits, texture images and 3D shapes.

Figure 3 .

 3 Figure3. The first row corresponds to the orbit recognition and the texture classification while the second row corresponds to 3D shape segmentation. On each row, the left plot shows the dependence of the accuracy on the number of directions, the middle plot shows the dependence of a single Gram matrix computation time, and the right plot shows the dependence of the ratio of the approximation of SW and the exact SW. Since the box plot of the ratio for orbit recognition is very similar to that of 3D shape segmentation, we only give the box plot of texture classification in the first row.

Figure 4 .

 4 Figure4. We show how the metric d1 is distorted. Each point represents a pair of PDs and its abscissae is the first diagram distance between them. Depending on the point color, its ordinate is the logarithm of the distance between PDs in the RKHS induced by either kPSS (blue points), kPWG (green points), kSW (red points) and a Gaussian kernel on d1 (black points).

Table 1 .

 1 Number of instances in the training set, the test set and number of labels of the different applications.

	+ s;
	end for
	Output: (1/π)SW M ;
	TASK	TRAINING TEST LABELS
	ORBIT	175	75	5
	TEXTURE	240	240	24
	HUMAN	415	1618	8
	AIRPLANE	300	980	4
	ANT	364	1141	5
	BIRD	257	832	4
	FOURLEG	438	1097	6
	OCTOPUS	334	1447	2
	FISH	304	905	3

INRIA Saclay

CREST, ENSAE, Université Paris Saclay.Correspondence to: Mathieu Carrière <mathieu.carriere@inria.fr>.

Acknowledgements. We thank the anonymous referees for their insightful comments. SO was supported by ERC grant Gudhi and by ANR project TopData. MC was supported by a chaire de l'IDEX Paris Saclay.