domain: info math

text:
Keanu Sisouk

Julie Delon

Julien Tierny

Wasserstein Dictionaries of Persistence Diagrams

Keywords: Topological data analysis, ensemble data, persistence diagrams. ✦

This paper presents a computational framework for the concise encoding of an ensemble of persistence diagrams, in the form of weighted Wasserstein barycenters [100], [102] of a dictionary of atom diagrams. We introduce a multi-scale gradient descent approach for the efficient resolution of the corresponding minimization problem, which interleaves the optimization of the barycenter weights with the optimization of the atom diagrams. Our approach leverages the analytic expressions for the gradient of both sub-problems to ensure fast iterations and it additionally exploits shared-memory parallelism. Extensive experiments on public ensembles demonstrate the efficiency of our approach, with Wasserstein dictionary computations in the orders of minutes for the largest examples. We show the utility of our contributions in two applications. First, we apply Wassserstein dictionaries to data reduction and reliably compress persistence diagrams by concisely representing them with their weights in the dictionary. Second, we present a dimensionality reduction framework based on a Wasserstein dictionary defined with a small number of atoms (typically three) and encode the dictionary as a low dimensional simplex embedded in a visual space (typically in 2D). In both applications, quantitative experiments assess the relevance of our framework. Finally, we provide a C++ implementation that can be used to reproduce our results.

INTRODUCTION

A S measurement devices and numerical techniques are becom- ing more and more advanced, datasets are becoming more and more complex geometrically. This geometrical complexity makes interactive exploration and analysis difficult, which challenges the interpretation of the data by the users. This motivates the creation of expressive data abstractions, capable of encapsulating the main features of interest of the data into simple representations, visually conveying the main information to the user.

Topological Data Analysis (TDA) [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] is a family of techniques which precisely addresses this issue. It provides concise topological descriptors of the main structural features hidden in a dataset. The relevance of TDA for analyzing scalar data, its efficiency and robustness have been documented in a number of visualization tasks [START_REF] Heine | A survey of topology-based methods in visualization[END_REF]. Examples of successful applications include turbulent combustion [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF], [START_REF] Gyulassy | Stability of dissipation elements: A case study in combustion[END_REF], [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF], material sciences [START_REF] Favelier | Visualizing ensembles of viscous fingers[END_REF], [START_REF] Gyulassy | Topologically Clean Distance Fields[END_REF], [START_REF] Gyulassy | Interstitial and Interlayer Ion Diffusion Geometry Extraction in Graphitic Nanosphere Battery Materials[END_REF], [START_REF] Soler | Ranking Viscous Finger Simulations to an Acquired Ground Truth with Topology-Aware Matchings[END_REF], nuclear energy [START_REF] Maljovec | Topology-inspired partition-based sensitivity analysis and visualization of nuclear simulations[END_REF], fluid dynamics [START_REF] Kasten | Two-dimensional timedependent vortex regions based on the acceleration magnitude[END_REF], [START_REF] Nauleau | Topological Analysis of Ensembles of Hydrodynamic Turbulent Flows -An Experimental Study[END_REF], bioimaging [START_REF] Anderson | Topological data analysis of functional MRI connectivity in time and space domains[END_REF], [START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF], [START_REF] Carr | Simplifying Flexible Isosurfaces Using Local Geometric Measures[END_REF], chemistry [START_REF] Bhatia | Topoms: Comprehensive topological exploration for molecular and condensed-matter systems[END_REF], [START_REF] Guenther | Characterizing Molecular Interactions in Chemical Systems[END_REF], [START_REF] Olejniczak | A Topological Data Analysis Perspective on Non-Covalent Interactions in Relativistic Calculations[END_REF], [START_REF] Olejniczak | Topological Data Analysis of Vortices in the Magnetically-Induced Current Density in LiH Molecule[END_REF] or astrophysics [START_REF] Shivashankar | Felix: A topology based framework for visual exploration of cosmic filaments[END_REF], [START_REF] Sousbie | The Persistent Cosmic Web and its Filamentary Structure: Theory and Implementations[END_REF].

Among the different topological descriptors studied in TDA (such as the merge and contour trees [START_REF] Acharya | A parallel and memory efficient algorithm for constructing the contour tree[END_REF], [START_REF] Carr | Computing contour trees in all dimensions[END_REF], [START_REF] Carr | Parallel peak pruning for scalable SMP contour tree computation[END_REF], [START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF], [START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF], [START_REF] Tarasov | Construction of contour trees in 3D in O(n log n) steps[END_REF], the Reeb graph [START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF], [START_REF] Doraiswamy | Computing Reeb Graphs as a Union of Contour Trees[END_REF], [START_REF] Gueunet | Task-based Augmented Reeb Graphs with Dynamic ST-Trees[END_REF], [START_REF] Parsa | A deterministic o(m log m) time algorithm for the reeb graph[END_REF], [START_REF] Pascucci | Robust on-line computation of Reeb graphs: simplicity and speed[END_REF], [START_REF] Tierny | Loop surgery for volumetric meshes: Reeb graphs reduced to contour trees[END_REF] , or the Morse-Smale complex [START_REF] Bremer | A Multi-Resolution Data Structure for 2-Dimensional Morse Functions[END_REF], [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF], [START_REF] Edelsbrunner | Morse-Smale complexes for piecewise linear 3-manifolds[END_REF], [START_REF] Edelsbrunner | Hierarchical morse complexes for piecewise linear 2-manifolds[END_REF], [START_REF] Forman | A User's Guide to Discrete Morse Theory[END_REF], [START_REF] Gyulassy | Shared-Memory Parallel Computation of Morse-Smale Complexes with Improved Accuracy[END_REF], [START_REF] Robins | Theory and Algorithms for Constructing Discrete Morse Complexes from Grayscale Digital Images[END_REF], [START_REF] Shivashankar | Parallel Computation of 3D Morse-Smale Complexes[END_REF]), the Persistence Diagram (Fig. 1) is a particularly prominent example. As described in Sec. 2.1, it is a concise topological descriptor which captures the main structural features in a dataset and which assesses their individual importance.

In addition to the challenge of increased geometrical complexity (discussed above), a new difficulty has recently emerged in many applications, with the notion of ensemble dataset. These representations describe a given phenomenon not only with a single dataset, but with a collection of datasets, called ensemble members. In that context, the topological analysis of an ensemble dataset consequently results in an ensemble of corresponding topological descriptors (e.g. one persistence diagram per ensemble member).

Then, a major challenge consists in developing practical tools for such an ensemble of topological descriptors, to facilitate its processing, analysis and visualization. Such tools include compression approaches (to facilitate the manipulation of the ensemble of descriptors) or visualization methods (for instance, with planar layouts, where each point encodes a descriptor and the distance between a pair of points encodes the intrinsic differences between the corresponding descriptors).

To enable the above tools, a key research question deals with the definition of a concise, yet informative, encoding of the ensemble of descriptors. A promising research direction consists in defining a dictionary (i.e. a set of reference descriptors, or atoms), such that the topological descriptors of the ensemble can be concisely encoded by expressing them as a specific function of the atoms (e.g. a linear combination). At a technical level, this requires to accurately capture and model the implicit relations (i.e. the possible functions) which link the different descriptors of the ensemble.

A series of recent works started the exploration of this overall direction, in particular with the notion of average topological representation [START_REF] Lacombe | Large Scale computation of Means and Clusters for Persistence Diagrams using Optimal Transport[END_REF], [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF], [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF], [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF], [START_REF] Yan | A structural average of labeled merge trees for uncertainty visualization[END_REF]. These techniques can produce a topological descriptor which nicely summarizes the ensemble. However, they do not capture the implicit relations between the different topological descriptors.

This paper addresses this issue by introducing a simple and efficient approach for the estimation of linear relations between persistence diagrams on their associated Wasserstein metric space. Inspired by previous work on histograms [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning[END_REF], our approach provides a linear encoding of the input ensemble, where each diagram is represented as a weighted Wasserstein barycenter [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF], [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] of a dictionary of automatically optimized diagrams called atom diagrams. We introduce a novel multi-scale gradient descent algorithm (Sec. 4) for the efficient resolution of the corresponding minimization problem (Sec. [START_REF] Anderson | Topological data analysis of functional MRI connectivity in time and space domains[END_REF], for which we interleave the optimization of the barycenter weights (Sec. 3.2) with the optimization of the atom diagrams (Sec. 3.3). Extensive experiments (Sec. 6) on public ensembles demonstrate the efficiency of our approach, with Wasserstein dictionary computations in the orders of minutes for the largest examples. We illustrate the relevance of our contributions for the visual analysis of ensemble data with two applications, data reduction (Sec. 5.1) and dimensionality reduction (Sec. 5.2).

Related Work

The literature related to our work can be classified into three main classes: (i) uncertainty visualization, (ii) ensemble visualization, and (iii) topological methods for ensembles. (i) Uncertainty visualization: Data variability can be represented in the form of uncertain datasets, by considering the data at each point of the domain as a random variable, associated with an explicit probability density function (PDF). The analysis and visualization of uncertain data has been recognized as a major challenge in the visualization community [START_REF]2008 uncertainty of measurement -part 3: Guide to the expression of uncertainty in measurement (GUM)[END_REF], [START_REF] Bonneau | Overview and state-of-the-art of uncertainty visualization[END_REF], [START_REF] Johnson | A next step: Visualizing errors and uncertainty[END_REF], [START_REF] Maceachren | Visualizing geospatial information uncertainty: What we know and what we need to know[END_REF], [START_REF] Pang | Approaches to uncertainty visualization[END_REF], [START_REF] Potter | From quantification to visualization: A taxonomy of uncertainty visualization approaches[END_REF]. Several techniques have been proposed either dealing with the entropy of the random variables [START_REF] Potter | Visualization of uncertainty without a mean[END_REF], or their correlation [START_REF] Pfaffelmoser | Visualization of global correlation structures in uncertain 2D scalar fields[END_REF] or gradient variation [START_REF] Pfaffelmoser | Visualizing the variability of gradients in uncertain 2D scalar fields[END_REF]. The effect of data uncertainty on feature extraction has also been studied (for instance for level set extraction [START_REF] Athawale | Uncertainty quantification in linear interpolation for isosurface extraction[END_REF], [START_REF] Athawale | Isosurface visualization of data with nonparametric models for uncertainty[END_REF], [START_REF] Pfaffelmoser | Visualizing the positional and geometrical variability of isosurfaces in uncertain scalar fields[END_REF], [START_REF] Pöthkow | Positional Uncertainty of Isocontours: Condition Analysis and Probabilistic Measures[END_REF]- [START_REF] Pöthkow | Probabilistic Marching Cubes[END_REF], [START_REF] Schlegel | On the interpolation of data with normally distributed uncertainty for visualization[END_REF]), for various interpolation schemes and PDF models (e.g. Gaussian [START_REF] Liebmann | Critical Points of Gaussian-Distributed Scalar Fields on Simplicial Grids[END_REF], [START_REF] Otto | Uncertain 2D vector Field Topology[END_REF], [START_REF] Otto | Uncertain topology of 3D vector fields[END_REF], [START_REF] Petz | Probabilistic local features in uncertain vector fields with spatial correlation[END_REF] or uniform [START_REF] Bhatia | Flow visualization with quantified spatial and temporal errors using edge maps[END_REF], [START_REF] Günther | Mandatory critical points of 2D uncertain scalar fields[END_REF], [START_REF] Szymczak | Hierarchy of stable Morse decompositions[END_REF] distributions). In general, a central limitation of existing methods for uncertain data is their design dependence on a specific PDF model (Gaussian, uniform, etc). This challenges their usability for ensemble data, where the PDFs estimated from the ensemble can follow an arbitrary, unknown model. Moreover, most of these techniques do not consider multi-modal PDFs, which are however essential when multiple trends appear in the ensemble. (ii) Ensemble visualization: Another approach to model data variability consists in using ensemble datasets. In this context, the variability is encoded by a sequence of empirical observations (i.e the members of the ensemble). Established techniques typically compute geometrical objects, such as level sets or streamlines, thereby capturing the main features for each member of the ensemble. From there, a representative of the resulting ensemble of geometrical objects can be computed. For this task, a few methods have been introduced. For instance spaghetti plots [START_REF] Diggle | The Analysis of Longitudinal Data[END_REF] are used in the case of level-set variability, more particularly for weather data [START_REF] Potter | Ensemble-vis: A framework for the statistical visualization of ensemble data[END_REF], [START_REF] Sanyal | Noodles: A tool for visualization of numerical weather model ensemble uncertainty[END_REF], and box-plots [START_REF] Mirzargar | Curve boxplot: Generalization of boxplot for ensembles of curves[END_REF], [START_REF] Whitaker | Contour boxplots: A method for characterizing uncertainty in feature sets from simulation ensembles[END_REF] for the variability of contours and curves. In the case of trend variability, Hummel et al. [START_REF] Hummel | Comparative visual analysis of lagrangian transport in CFD ensembles[END_REF] conceived a Lagrangian framework for classification purposes in flow ensembles. More specifically, clustering techniques have been used to identify the main trends in ensemble of streamlines [START_REF] Ferstl | Streamline variability plots for characterizing the uncertainty in vector field ensembles[END_REF] and isocontours [START_REF] Ferstl | Visual analysis of spatial variability and global correlations in ensembles of iso-contours[END_REF]. However, only few techniques have applied this strategy to topological objects. Favelier et al. [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF] and Athawale et al. [START_REF] Athawale | Uncertainty Visualization of 2D Morse Complex Ensembles Using Statistical Summary Maps[END_REF] respectively introduced techniques to analyze the geometrical variability of critical points and gradient separatrices. Overlapbased heuristics have been studied for estimating a representative contour tree from an ensemble [START_REF] Kraus | Visualization of uncertain contour trees[END_REF], [START_REF] Wu | A contour tree based visualization for exploring data with uncertainty[END_REF]. In the context of ensembles of histograms, Schmitz et al. [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning[END_REF] introduced a dictionary encoding approach based on optimal transport [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. However, this method is not directly applicable to persistence diagrams. It focuses on a fundamentally different object (histograms). Thus, the employed distances, geodesics and barycenters are defined differently (in particular in an entropic form [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF]) and the algorithms for their computations are drastically different (based on Sinkhorn matrix scaling [START_REF] Sinkhorn | Diagonal equivalence to matrices with prescribed row and column sums[END_REF]). In contrast, our work focuses on Persistence diagrams (Sec. 2.1), whose associated metric space is also inspired from optimal transport, but with various formal and computational specificities (Sec. 2.2). Moreover, our approach is based on gradient descent which, from our experience, provides better practical convergence for this kind of problems than quasi-Newton techniques. Finally, we contribute a multi-scale progressive optimization algorithm, which provides improved solutions in comparison to a naive optimization. (iii) Topological methods for ensembles: To analyze the relations between the persistence diagrams of an ensemble, several key low level notions are required, such as the notion of distance and barycenters between diagrams, for which we review the literature here. Inspired by optimal transport [START_REF] Kantorovich | On the translocation of masses[END_REF], [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF], the Wasserstein distance between persistence diagrams [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] (Sec. 2.2) has been extensively studied [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF], [START_REF] Cohen-Steiner | Lipschitz Functions Have Lp-Stable Persistence[END_REF]. It relies on a bipartite assignment problem, for which exact [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] and approximate [START_REF] Bertsekas | A new algorithm for the assignment problem[END_REF], [START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF] implementations are available in open-source [START_REF] Tierny | The Topology ToolKit[END_REF]. Based on this distance, several approaches have explored the possibility to define a representative diagram of an ensemble of persistence diagrams, with the notion of Wasserstein barycenter. Turner et al. [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF] introduced the first approach for the computation of such a barycenter. Lacombe et al. [START_REF] Lacombe | Large Scale computation of Means and Clusters for Persistence Diagrams using Optimal Transport[END_REF] presented an approach based on entropic transport [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF]. However, it requires a pre-vectorization step which is subject to several parameters, and which is not conducive to visualization tasks (features can no longer be individually tracked beyond the pre-vectorization step). In contrast, Vidal et al. [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] introduced a vectorization-free approach which maintains the feature assignments explicitly. It is based on a progressive scheme, which greatly accelerates computation in practice. These concepts have been recently investigated for other topological descriptors, such as merge trees [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF], [START_REF] Yan | A structural average of labeled merge trees for uncertainty visualization[END_REF]. Recently, several authors have investigated another compact representation of ensembles of topological descriptors, via a basis of representative descriptors. For instance, Li et al. [START_REF] Li | Sketching merge trees for scientific visualization[END_REF] introduce a vectorization for merge trees, which was subsequently used by matrix sketching procedures [START_REF] Woodruff | Sketching as a Tool for Numerical Linear Algebra[END_REF] to create a basis of representative merge trees. In contrast, our work focuses on persistence diagrams (which can encode different features). Also, it directly operates on the Wasserstein metric space of persistence diagrams, thereby avoiding the typical technical difficulties associated with vectorizations (e.g. quantization and/or linearization artifacts, potential stability issues, possible inaccuracies in vectorization reversal, etc.). Pont et al. [START_REF] Pont | Principal geodesic analysis of merge trees (and persistence diagrams)[END_REF] introduced the notion of principal geodesic analysis of merge trees (and persistence diagrams), with the same overall goal of characterizing the relations between the topological descriptors of an ensemble. In this work, we introduce a different formulation of the problem, which is both simpler (based on the construction of weighted Wasserstein barycenters) and more flexible (our optimization is not subject to complicated constraints such as geodesic orthogonality). This results in a simpler implementation and slightly faster computations (Sec. 6).

Contributions

This paper makes the following new contributions:

1) A simple approach for the linear encoding of Persistence Diagrams: We formulate the linear encoding of an ensemble of persistence diagrams on their associated Wasserstein metric space as a dictionary optimization (Sec. 3), which simply optimizes, simultaneously, (i) the barycentric weights (Sec. 3.2) and (ii) the atoms of the dictionary (Sec. 3.3). 2) A multi-scale algorithm for the computation of a Wasserstein dictionary of Persistence Diagrams: We introduce a novel, efficient algorithm for the optimization of the above dictionary encoding problem. Our algorithm leverages the analytic expressions of the gradient of both of the above sub-problems, to ensure fast iterations. Moreover, in comparison to a naive optimization, our algorithm reaches solutions of improved energy thanks to a multi-scale strategy. Finally, we leverage shared-memory parallelism to further improve performances. 3) An application to data reduction: We present an application to data reduction (Sec. 5.1), where the persistence diagrams of the input ensemble are significantly compressed, by solely storing their barycentric weights as well as the atom diagrams. 4) An application to dimensionality reduction: We present an application to dimensionality reduction (Sec. 5.2), by embedding each input diagram as a point within a 2D view, based on its weights relative to a Wasserstein dictionary composed of three atoms (thereby defining a 2-simplex). 5) Implementation: We provide a C++ implementation of our algorithms that can be used for reproducibility purposes.

PRELIMINARIES

This section presents the theoretical elements needed for the formalization of our work. We introduce the topological data representation that we use -the persistence diagram (Sec. 2.1)and its associated metric (Sec. 2.2). Then we define the notion of Wasserstein barycenter of persistence diagrams (Sec. 2.3), which is a core component of our approach (Sec. 3).

Persistence diagrams

Each input ensemble member is given in the form of a piecewise linear (PL) scalar field f : M → R defined on a PL (d M )-manifold M , with d M ≤ 3 for our applications. Given an isovalue w ∈ R,

we denote f -1 -∞ (w) = f -1 (-∞, w
] the sub-level set of f at w. As w increases, the topology of f -1 -∞ changes at specific points of M , called "critical points". Critical points are classified by their index I : 0 for minima, 1 for 1-saddles, d M -1 for (d M -1)-saddles, and d M for maxima (in practice, f is enforced to contain only isolated, non-degenerate critical points [START_REF] Edelsbrunner | Hierarchical morse complexes for piecewise linear 2-manifolds[END_REF], [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF]). According to the Elder rule [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF], each topological feature of f -1 -∞ (w) (e.g. a connected component, a cycle, a void) can be associated with a pair of critical points (c, c ′ ) (with f (c) < f (c ′ ) and I c = I c ′ -1), corresponding to its birth and death during the sweep of the data by w (from -∞ to +∞). Such a pair (c, c ′ ) is called a persistence pair. For instance, when two connected components of f -1 -∞ (w) meet at a critical point c ′ , the younger one (created last) dies, letting the oldest one (created first) survive. Then, the critical points are represented visually as 2D bar codes where the horizontal axis encodes the birth of a feature noted b = f (c) and where the vertical axis encodes its lifespan and death noted d = f (c ′ ) . This representation is called the Persistence Diagram, noted X. In the diagram, salient features stand out from the diagonal and smallamplitude noise is typically located near the diagonal, as shown in Fig. 1. In the remainder, we enumerate the points of X with indices such that X = {x 1 , . . . , x K } and we note i X = {1, . . . , K} the set of indices (i.e. the set of all integers going from 1 to K) .

Wasserstein distance

To evaluate the distance between two diagrams, a typical preprocessing step consists in augmenting each diagram, by including the diagonal projection of all the off-diagonal points of the other diagram. To illustrate that, let us consider

X 1 = {x 1 1 , . . . , x K 1 1 } and X 2 = {x 1 2 , . . . , x K 2 2 }.
Given an off-diagonal point x (i.e. b < d), let ∆x be its diagonal projection, specifically: ∆x = ( b+d 2 , b+d 2 ). Let P 1 and P 2 be the sets of the diagonal projections of the points of X 1 and X 2 respectively. Then, X 1 and X 2 are augmented into X 1 ′ and X 2

′ by considering

X 1 ′ = X 1 ∪ P 2 and X 2 ′ = X 2 ∪ P 1 . This ensures that |X ′ 1 | = |X ′ 2 | = K (
which eases distance evaluation). Specifically, we consider in the remainder that the notations X 1 and X 2 refer to augmented diagrams (i.e.

|X 1 | = |X 2 | = K).
Then, given two persistence diagrams X 1 and X 2 the L 2 -Wasserstein distance between them is defined as:

W (X 1 , X 2 ) = min ψ:i X bi j -→i X K ∑ j=1 c(x j 1 , x ψ( j) 2 ), (1)
where ψ, the matching, is a bijection of the index set i X towards itself (i.e. ψ is a permutation of i X ). This bijection matches one persistence pair x 1 of X 1 (i.e. a pair of critical points of index I and I + 1 respectively) to one persistence pair x 2 in X 2 (another pair of critical points, with the same indices I and I + 1) whether x 1 and x 2 are on the diagonal or not (Fig. 2). Given the cost c in the definition of W , ψ is the optimal way to transport X 1 onto X 2 .

In our work, we consider the cost c(x, y) = 0 when x and y are both diagonal points, and c(x, y) = ∥x -y∥ 2 otherwise (∥x -y∥ denotes the Euclidean distance between x and y in the birth/death space).

Wasserstein barycenter

Given a set of persistence diagrams D = {a 1 , . . . , a m } (which we will call in the remainder dictionary), a Wasserstein barycenter (Fig. 3) -or Fréchet mean -of the dictionary D with barycentric weights λ λ λ = (λ 1 , . . . , λ m ) is a diagram, which we note Y (λ λ λ , D) in the following, which minimizes the Fréchet energy E F (B):

E F (B) = m ∑ i=1 λ i W 2 (a i , B)
.

λ λ λ is such that λ i ≥ 0 and m ∑ i=1 λ i = 1.
We denote Σ m the simplex of such vectors. Intuitively, Y (λ λ λ , D) is a diagram which minimizes the above linear combination, given λ λ λ , of its squared Wasserstein distances to the diagrams of the dictionary D.

The computation of the barycenter Y (λ λ λ , D) requires generalizing the pairwise augmentation described in Sec. 2.2. Specifically, each non-diagonal point of each dictionary diagram a i is projected to the diagonal of all the other dictionary diagrams a j (with i ̸ = j). After this first augmentation, each dictionary diagram a i contains K = ∑ m i=1 |a i | points (where |a i | is the number of nondiagonal points in a i ). Then Y (λ λ λ , D) is typically initialized on the dictionary diagram a * which initially minimizes the Fréchet energy E F . Let |Y (λ λ λ , D)| = |a * | be the number of non-diagonal points of a * . Then, all the non-diagonal points of all the atoms are projected on the diagonal of Y (λ λ λ , D), and reciprocally, all the non-diagonal points of Y (λ λ λ , D) are projected on the diagonal of each atom. Thus, at this stage, after this second augmentation, each dictionary diagram a i and the candidate barycenter

Y (λ λ λ , D) contains K = ∑ m i=1 |a i | + |Y (λ λ λ , D)
| points (mostly on the diagonal). Next, we optimize Y (λ λ λ , D) in practice with the approach by Vidal et al. [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF], which provides a time-efficient approximation of the original algorithm by Turner et al. [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF]. Similar to Turner et al., it is based on an iterative optimization, where each iteration includes an Assignment step, followed by an Update step. Specifically, the Assignment step computes the optimal assignments ψ i between the candidate Y (λ λ λ , D) and each dictionary diagram a i . Next, the Update step minimizes the Fréchet energy E F under the current assignments ψ i . Since the L 2 -Wasserstein distance considers the Euclidean distance as a cost function (Sec. 2.2), this minimization is achieved by simply placing each point of Y (λ λ λ , D) at the arithmetic mean in the birth/death space of its assigned points in the dictionary diagrams. This can be done since the arithmetic mean minimizes the Fréchet energy defined respectively to Euclidean distances (more sophisticated Update procedures, e.g. based on an optimization routine, would need to be derived for other distances in the birth/death space). After this Update step, the subsequent Assignment further improves the assignments ψ i , hence decreasing the Fréchet energy constructively at each iteration.

The approach by Vidal et al. [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] revisits this framework by integrating tailored approximations throughout the computation. Specifically, it approximates the optimal assignments ψ i with the fast Auction optimization [START_REF] Bertsekas | A new algorithm for the assignment problem[END_REF] (instead of the traditional, yet prohibitive, Munkres algorithm [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]). Further, it improves performance with a mechanism called price memorization, which enables the initialization of the Auction optimization with the assignments ψ i computed in the previous Assignment step. This allows the barycenter optimization to resume the assignment optimization instead of re-computing it from scratch at each iteration. This approach also includes a strategy for the adaptive increase of the accuracy parameter of the Auction optimization, allowing for fast assignments in the early iterations of the barycenter algorithm, and slower but more accurate assignments towards its convergence.

WASSERSTEIN DICTIONARY ENCODING

This section formalizes our approach for the Wasserstein dictionary encoding of an ensemble of persistence diagrams. Sec. 3.1 provides an overview of our approach, which interleaves barycentric weight optimization (λ λ λ ) with atom optimization (D). Finally, Secs. 3.2 and 3.3 detail the gradient estimation for both sub-problems.

Overview

Let {X 1 , . . . , X N } be the input ensemble of N persistence diagrams. The goal of our approach is to jointly optimize two sub-problems:

• Optimize a set D of m reference persistence diagrams, called the atoms of the Wasserstein dictionary D; • Optimize for each input diagram X n a vector of m barycentric weights λ λ λ n ∈ Σ m , in order to accurately approximate X n with a Wasserstein barycenter Y (λ λ λ n , D) (Sec. 2.3). This can be formalized as a joint optimization, where one wishes to find the optimal barycentric weights Λ * = λ λ λ * 1 , . . . , λ λ λ * N and the optimal Wasserstein dictionary D * = {a * 1 , . . . , a * m } (with m ≪ N), in order to minimize the following dictionary energy:

E D (Λ, D) = N ∑ n=1 W 2 Y (λ λ λ n , D), X n . (2)
Our overall strategy for optimizing Eq. 2 consists in iteratively interleaving two sub-optimizations:

1) For a fixed dictionary D, the set of barycentric weights Λ is optimized with one step of gradient descent (Sec. 3.2); 2) For a fixed set of barycentric weights Λ, the dictionary D is optimized with one step of gradient descent (Sec. 3.3). Then, this sequence of two sub-procedures is iterated until a predefined stopping condition is reached (Sec. 4.2).

Finally, the output of our approach is the optimized Wasserstein dictionary D * (a set of m atom diagrams) and, for each input diagram X n , a vector of weights λ λ λ * n ∈ Σ m , which can be interpreted as the barycentric coordinates of X n in D * (thus capturing linear relations between the input diagrams on the Wasserstein dictionary).

Weight optimization

This section details the optimization of the barycentric weights Λ = λ λ λ 1 , . . . , λ λ λ N . Let D = {a 1 , . . . , a m } be a fixed dictionary of atom diagrams, with m > 0. Let X be a diagram of the input ensemble. For a given set of weights λ λ λ = λ 1 , . . . , λ m , let Y (λ λ λ ) = y 1 (λ λ λ ), . . . , y K (λ λ λ ) be its barycentric approximation, relative to D (i.e. each point y j (λ λ λ ) of Y (λ λ λ ) approximates a point in X).

We recall that after augmentation (Sec. Then, in order to compare it to X, Y (λ λ λ ) is further augmented by projecting on its diagonal the |X| non-diagonal points of X. Then, at this stage, the size

K of Y (λ λ λ ) is given by K = ∑ m i=1 |a i | + |Y (λ λ λ )| + |X|.
We augment similarly X (i.e. by projecting the nondiagonal points of Y (λ λ λ ) to its diagonal) and the m atoms (i.e. by projecting the non-diagonal points of X to their diagonals). Then, at this point, Y (λ λ λ ), X, and the m atoms a i all have the same size

K = ∑ m i=1 |a i | + |Y (λ λ λ )| + |X|.
In this section, we describe a gradient descent on λ λ λ to minimize the weight energy:

E W (λ λ λ ) = W 2 Y (λ λ λ ), X . (3)
A step of the corresponding gradient descent is illustrated in Fig. 4.

Given the set of optimal matchings φ 1 , . . . , φ m between Y (λ λ λ ) and the m atoms, the j th point of Y (λ λ λ ), noted y j (λ λ λ ), is given by:

∀ j ∈ {1, . . . , K}, y j (λ λ λ ) = m ∑ i=1 λ i a φ i ( j) i . (4)
In other words, the j th point y j (λ λ λ ) of the diagram Y (λ λ λ ) is a linear combination (with the weights λ λ λ ) of the m points it matches to in the atoms (one point per atom a i ), under the optimal assignments φ i (i.e. minimizing Eq. 1).

For a fixed set of assignments φ 1 , . . . , φ m , the Wasserstein distance (Eq. 1) between X and its approximation Y (λ λ λ ) is then:

E W (λ λ λ ) = W 2 Y (λ λ λ ), X = K ∑ j=1 c y j (λ λ λ ), x ψ( j) ,
where ψ denotes the optimal assignment (Eq. 1) between X and its approximation Y (λ λ λ ). When y j (λ λ λ ) and x ψ( j) are not both diagonal points, the cost c y j (λ λ λ ), x ψ( j) is given by their squared Euclidean distance in the birth/death space (it is zero otherwise, see Sec. 2.2). Then, by exploiting Eq. 4, E W (λ λ λ ) can be re-written as:

E W (λ λ λ ) = W 2 Y (λ λ λ ), X = K ∑ j=1 ∥y j (λ λ λ ) -x ψ( j) ∥ 2 = K ∑ j=1 ∥ m ∑ i=1 λ i a φ i ( j) i -x ψ( j) ∥ 2 . Since m ∑ i=1 λ i = 1, E W (λ λ λ
) can finally be re-written as:

E W (λ λ λ ) = W 2 Y (λ λ λ ), X = K ∑ j=1 ∥ m ∑ i=1 λ i (a φ i ( j) i -x ψ( j) )∥ 2 .
(5) Intuitively, this energy measures the error (in terms of Wasserstein distance) induced by approximating the input diagram X with its barycentric approximation Y (λ λ λ ). In Eq. 5, it is computed for each j th point y j (λ λ λ ) of the diagram Y (λ λ λ ), by considering the birth/death distances between the points y j (λ λ λ ) maps to, in the atoms on one hand and in the input diagram X on the other.

Then, by applying the chain rule on Eq. 5, the gradient of the weight energy (Eq. 3) is given by:

∇E W (λ λ λ ) = 2 m ∑ i=1 K ∑ j=1     (a φ 1 ( j) 1 -x ψ( j) ) T . . . (a φ m ( j) m -x ψ( j) ) T     λ i (a φ i ( j) i -x ψ( j) ) . (6)
Now that the gradient of the weight energy is available (Eq. 6), we can proceed to gradient descent. Specifically, the barycentric weights at the iteration t + 1 (noted λ λ λ t+1 ) are obtained by a step ρ λ λ λ from the weights at the iteration t (noted λ λ λ t ) along the gradient:

λ λ λ t+1 = Π Σ m λ λ λ t -ρ λ λ λ ∇E W (λ λ λ t ) , (7)
where Π Σ m is the projection onto the simplex of admissible barycentric weights (i.e. positive and summing to 1, c.f. Sec. 2.3). Since ∇E W is L-Lipschitz (see the computation details in Appendix A), a gradient step will guarantee an energy decrease as long as:

ρ λ λ λ ≤     2 K ∑ j=1 (a φ 1 ( j) 1 -x ψ( j) ) T . . . (a φ m ( j) m -x ψ( j) ) T 2     -1 < 1 L .
Overall, for a given input diagram X, each iteration t of gradient descent for the optimization of E W consists in the following steps:

1) Computing the Wasserstein barycenter Y (λ λ λ t ) (Sec. 2.3); 2) Computing the Wasserstein distance W 2 Y (λ λ λ t ), X (Eq. 3); 3) Estimating the gradient ∇E W (λ λ λ ) (Eq. 6); 4) Applying one step ρ λ λ λ of gradient descent (Eq. 7).

Atom optimization

This section details the optimization of the atoms of the dictionary D = {a 1 , . . . , a m }. Similarly to Sec. 3.2, let X be a diagram of the input ensemble and let λ λ λ ∈ Σ m be its -fixed -vector of barycentric weights. For a given dictionary D, let Y (D) = y 1 (D), . . . , y K (D) be the barycentric approximation of X, relative to λ λ λ . In this section, we describe a step of gradient descent on D to minimize the following atom energy:

E A (D) = W 2 Y (D), X .
A step of the corresponding gradient descent is illustrated in Fig. 5.

Given the set of optimal matchings φ 1 , . . . , φ m between Y (D) and the m atoms, the j th point of Y (D), noted y j (D), is given by:

∀ j ∈ {1, . . . , K}, y j (D) = m ∑ i=1 λ i a φ i ( j) i .
This expression is identical to Eq. 4 (Sec. 3.2). However, y j now depends on D, which is the variable of the current optimization. Then, the gradient of y j (D) with regard to D is simply given by:

∇y j (D) = λ 1 • • • λ m T . (8)
For a fixed set of assignments φ 1 , . . . , φ m , the Wasserstein distance (Eq. 1) between X and its approximation Y (λ λ λ ) is then:

E A (D) = W 2 Y (D), X = K ∑ j=1 c y j (D), x ψ( j) ,
where ψ( j) denotes the optimal assignment between X and its barycentric approximation Y (D). Similarly to Eq. 5 (Sec. 3.2), the above equation can be re-written as:

E A (D) = W 2 Y (D), X = K ∑ j=1 ∥ m ∑ i=1 λ i (a φ i ( j) i -x ψ( j) )∥ 2 . Let D j = [a φ 1 ( j) 1 , . . . , a φ m ( j) m
] T be the (m × 2)-matrix formed by the atom points matching to a given point y j (D) of Y (D), via the fixed assignments φ 1 , . . . , φ m . Specifically, the i th line of this matrix refers to the point a

φ i ( j) i
in the atom a i where y j (D) maps to (via the optimal assignment φ i ). For this line, the two columns of the matrix encode the birth/death coordinates of the point a φ i ( j) i

. Then, the pointwise atom energy of y j (D), noted e A (D j ), is given by:

e A (D j ) = m ∑ i=1 λ i (a φ i ( j) i -x ψ( j) ) 2 . ( 9
)
Then, by applying the chain rule on Eq. 9 (using Eq. 8), the gradient of the pointwise atom energy is given by:

∇e A (D j ) = 2 λ 1 • • • λ m T m ∑ i=1 λ i (a φ i ( j) i -x ψ( j) ) T . (10)
Now that the gradient of the pointwise atom energy is available (Eq. 10), we can proceed to a step of gradient descent. Specifically, the matrix of atom points matched to y j (D) at the iteration t + 1 (noted D j t+1 ) is obtained by a step ρ D from the same matrix at the iteration t (noted D j t ) along the gradient:

D j t+1 = Π X D j t -ρ D ∇e A (D j t ) , (11)
where Π X projects each atom point to an admissible region of the 2D birth/death space (i.e. above the diagonal and within the global scalar field range). Since ∇e A is L-Lipschitz (see the computation details in Appendix B), a gradient step will guarantee an energy decrease as long as:

ρ D < (4m) -1 < L -1 .
Note that, in order to control the final size S m of the dictionary D, after each iteration of atom optimization, each atom a i is thresholded by removing its K = (mK -S m )/m least persistent points (at the subsequent optimization iteration, all diagrams will be re-augmented again in a pre-preprocess, as detailed in Sec. 3.2).

Overall, for a given input diagram X, each iteration t of gradient descent for the optimization of E A consists in the following steps: Fig. 6. Illustration of our initialization strategy on a toy 2D point set (top left). First, the entries of the distance matrix of the input (middle) are summed on a per-line basis. The line maximizing this sum (cyan), noted l 1 , identifies the first atom, noted a 1 , as the point which is the furthest away from all the others (cyan sphere, top right). Next, the atom a 2 (grey sphere, top right) is selected as the point which maximizes its distance to a 1 . At this point, the line l 2 (grey, corresponding to the point a 2 ) is added to the line l 1 , to encode the distances to these two atoms (a 1 and a 2 ). Then, the point a 3 is selected as the maximizer of l 1 + l 2 : it is the point which is the furthest away from all the previously selected atoms. Then, the corresponding line, l 3 , is added to l 1 + l 2 and the process is iterated until the target number of atoms has been achieved.

ALGORITHM

This section presents our overall algorithm for the resolution of the optimization formulated in Sec. 3. Sec. 4.1 details our initialization strategy. Our overall multi-scale scheme is presented in Sec. 4.2. Finally, shared-memory parallelism is discussed in Sec. 4.3.

Initialization

Our strategy for the initialization of the Wasserstein dictionary D, illustrated in Fig. 6, is inspired by the celebrated k-means++ strategy [START_REF] Celebi | A comparative study of efficient initialization methods for the k-means clustering algorithm[END_REF]. Specifically, we iteratively select the m atoms among the N input diagrams. At the first iteration, we select as first atom the diagram which maximizes the sum of its Wasserstein distances (Eq. 1) to all the input diagrams (cyan point in Fig. 6). Next, each iteration selects as the next atom the diagram which maximizes the sum of its Wasserstein distances to all the previously selected atoms. This process stops when the desired number of atoms, m, has been selected. As illustrated in Fig. 6 in the case of a toy 2D point set, this initialization strategy has the nice property that it tends to select atoms on the convex hull of the input point set, which ensures that the non-atom points can indeed be expressed as a convex combination of the atoms, hence leading to accurate initial barycentric approximations. As for the barycentric weights, these are uniformly initialized (i.e. to 1/m).

Multi-scale optimization algorithm

In real-life data, persistence diagrams tend to contain many lowpersistence features, which essentially encode the noise in the data (see Fig. 1, right). In this section, we present a multi-scale optimization strategy which addresses this issue by prioritizing the optimization on the most persistent pairs, which correspond to the most salient features of the data. As detailed in Sec. 6.2, this strategy leads the optimization to solutions of improved energy in comparison to a naive (non-multi-scale) approach.

Our multi-scale strategy consists in iterating our optimization procedure by progressively increasing the resolution (in terms of persistence) of the input diagrams. This is inspired by the progressive strategy by Vidal et al. [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] for the problem of Wasserstein barycenter optimization. Specifically, given an input diagram X, let ∆ f be the span in scalar values in the corresponding ensemble member f (i.e. ∆ f = max v∈M f (v)min v∈M f (v)). Given a threshold τ ∈ [0, 1], we note X τ = {x ∈ X| d xb x ≥ τ∆ f } the version of X at resolution τ. It is a subset of X which contains persistence pairs whose relative persistence is above τ. Note that the input diagrams are not normalized by persistence, which would prevent the capture of variability in data ranges within the ensemble. Instead, we normalize the above persistence threshold, by expressing it as a fraction τ ∈ [0, 1] of the scalar field range ∆ f . Then, our multi-scale optimization will first consider the input diagrams at a resolution τ 0 and then will progressively consider finer resolutions τ 1 , . . . , τ r until the full diagrams are considered at τ r = 0. This multi-resolution strategy, based on a per-diagram normalized persistence threshold (τ ∈ [0, 1]) prevents diagrams from being empty in the early resolutions in case of large variations in data range within the ensemble (which would occur for instance with a per-ensemble normalization). The multi-resolution is illustrated in Fig. 7. In our experiments, we set τ 0 = 0.2 and decrease τ by 0.05 at each resolution (i.e. τ 1 = 0.15, τ 2 = 0.10, τ 3 = 0.05, τ 4 = 0). At each resolution, the solution for the previous resolution is used as an initialization. Note that alternative strategies were considered for decreasing τ (for instance by dividing it by 2 at each resolution), but the best experimental results were obtained for the above decrease strategy. Alg. 1 summarizes our overall approach. For each suboptimization (i.e. weight and atom optimization), although each gradient step is guaranteed to decrease the corresponding energy (see the end of Secs. Sec. 3.2 and Sec. 3.3), this is only true for fixed assignments (between a diagram X and its barycentric approximation as well as between the barycentric approximation and the atoms). Since the assignments can change along the iterations of the optimization, the overall energy E D (Eq. 2) may increase between consecutive iterations. Hence, pragmatic stopping conditions need to be considered. In practice, if E D has not decreased for more than 10 iterations, we return the solutions λ λ λ * and D * reached by the optimization with the lowest energy E D .

Parallelism

Our approach can be trivially parallelized with shared-memory parallelism. First, its most computationally demanding task, the N barycentric approximations of the input diagrams can be computed independently. Thus, for each barycentric approximation, we use one parallel task per input diagram. Next, the estimation of the gradient of E W (Sec. 3.2) is done on a per input diagram basis, independently. Thus, we use one parallel task per input diagram. Regarding the estimation of the gradient of E A (Sec. 3.3), given a barycentric approximation Y (D) of an input diagram X, each of its points y j (D) defines independently a pointwise version of the gradient of the atom energy (see the last paragraph of Sec. 3.3). Thus, we use one parallel task per point y j (D) of a barycentric approximation Y (D) of an input diagram X.

APPLICATIONS

This section illustrates the utility of our approach in concrete visualization tasks: data reduction and dimensionality reduction.

Data reduction

Like any data representation, persistence diagrams can benefit from lossy compression. This can be useful in in-situ [START_REF] Ayachit | ParaView Catalyst: Enabling In Situ Data Analysis and Visualization[END_REF] use-cases, where time-steps are represented on permanent storage with topological signatures [START_REF] Brown | Utilising urgent computing to tackle the spread of mosquito-borne diseases[END_REF]. In such scenarios, lossy compression is useful to facilitate the manipulation (i.e. storage and transfer) of the resulting ensemble of persistence diagrams. We present now an application to data reduction where the input ensemble of persistence diagrams is compressed, by only storing to disk: (i) the Wasserstein dictionary of persistence diagrams D * and (ii) the N barycentric weights λ λ λ * 1 , . . . , λ λ λ * N . The compression quality can be controlled with two input parameters (i) the number of atoms m and (ii) the maximum S m Fig. 8. Visual comparison (left) between the input persistence diagrams (top insets, saddle-maximum persistence pairs only) and our compressed diagrams (bottom insets, Sec. 5.1, saddle-maximum persistence pairs only) for three members of the Isabel ensemble (one member per ground-truth class). For each member, the sphere color encodes the matching between the input and the compressed diagrams (for the meaningful persistence pairs, above 10% of the function range). This visual comparison shows that the main features of the diagrams (encoding the main hurricane wind gusts in the data) are well preserved by the data reduction, especially for the members coming from the cluster 2, for which a lower relative reconstruction error (Err ) can be observed. The planar overview of the ensemble (right) generated by our dimensionality reduction (Sec. 5.2) enables the visualization of the relations between the different diagrams of the ensemble. Specifically, this illustration shows a larger disparity for two clusters. Fig. 9. Visual comparison (left) between the input persistence diagrams (top insets) and our compressed diagrams (bottom insets, Sec. 5.1) for four members of the Ionization front (3D) ensemble (one member per ground-truth class). The color encoding is the same as in Fig. 8. This visual comparison shows that the main features of the diagrams (the extremities of the ionization front) are well preserved by the data reduction, especially for the members coming from the clusters 2 and 3, for which a lower reconstruction error (Err ) can be observed. The planar overview of the ensemble (right) generated by our dimensionality reduction (Sec. 5.2) enables the visualization of the relations between the diagrams of the ensemble. Specifically, it shows a larger disparity for the clusters 0 and 1 (spread out purple and pink spheres), which are also the most difficult to reconstruct.

of the total size of the atoms (i.e. ∑ m i=1 |a i |). The reconstruction error (given by the energy E D , Eq. 2.) will be minimized for large values of both parameters, while the compression factor will be maximized for low values. In our data reduction experiments, we set the number of atoms m to the number of ground-truth classes of each ensemble, as documented in the ensemble descriptions [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF]. Moreover, we set S m to c -1 f ∑ N i=1 |X n |, where c f is a target compression factor and |X n | is the number of non-diagonal points in the input diagram X n (see Sec. 6.2 for a quantitative evaluation).

Fig. 8 (left) provides a visual comparison between the diagram compressed with this strategy (bottom insets) and the original diagram (top insets), for three members of the Isabel ensemble. This experiment shows that diagrams can be significantly compressed (c f = 5.49), while still faithfully encoding the main features of the data. Fig. 9 (left) provides a similar visual comparison for the Ionization front (3D) ensemble (c f = 2.9).

We have applied our data reduction approach to topological clustering [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF], where the main trends within the ensemble are identified by clustering the ensemble members based on their persistence diagrams. For the large majority of our test ensembles, the outcome of the clustering algorithm [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] was identical when used with the input diagrams or our compressed diagrams (Sec. 6.3 documents a counter-example). This confirms the viability and utility of our data reduction scheme.

Dimensionality reduction

Our framework can also be used to generate low-dimensional layouts of the ensemble, for its global visual inspection. Specifically, we generate 2D planar layouts by using m = 3 atoms and by embedding our Wasserstein dictionary D * as a triangle in the plane, such that its edge lengths are equal to the Wasserstein distances between the corresponding atoms. Next, each diagram X of the input ensemble is embedded as a point in this triangle by using its barycentric weights λ λ λ * as barycentric coordinates.

As illustrated in Figs. 8 (right) and 9 (right), our dimensionality reduction provides a planar overview of the ensemble which groups together diagrams which are close in terms of Wasserstein distances. Specifically, in both examples, the ground-truth classification of the ensemble is visually respected: the points of a given class (same color) indeed form a distinct cluster in the planar view.

RESULTS

This section presents experimental results obtained on a computer with two Xeon CPUs (3.2 GHz, 2x10 cores, 96GB of RAM). The input persistence diagrams were computed with the Discrete Morse Sandwich algorithm [START_REF] Guillou | Discrete morse sandwich: Fast computation of persistence diagrams for scalar data-an algorithm and a benchmark[END_REF]. We implemented our approach in C++ (with OpenMP), as modules for TTK [START_REF] Tierny | The Topology ToolKit[END_REF], [START_REF] Bin Masood | An Overview of the Topology ToolKit[END_REF]. Experiments were ran on the benchmark of public ensembles [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees -Ensemble Benchmark[END_REF] described in [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF], which includes simulated and acquired 2D and 3D ensembles from previous work and past SciVis contests [START_REF] Organizers | The IEEE SciVis Contest[END_REF]. The considered type of persistence pairs (i.e. the index of the corresponding critical points, Sec. 2.1) was selected on a per-ensemble basis, depending on the features of interest present in the ensemble. All types of pairs (i.e. minimum-saddle pairs, saddle-saddle pairs and saddlemaximum pairs) were considered for the following ensembles: Cloud processes, Isabel, Starting Vortex, Sea Surface Height, Vortex Street. Only the persistence pairs including extrema were considered for the ensembles Ionization front (2D) and Ionization front (3D). Finally, only the persistence pairs containing maxima were considered for the remaining ensembles: Asteroid Impact, Dark Matter, Earthquake, Viscous Fingering, Volcanic Eruptions.

Time performance

The most computationally expensive part of our approach is the computation of the N Wasserstein barycenters, for which we use the algorithm by Vidal et al. [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF]. Each iteration of barycenter optimization approximatively requires O(mK 2 ) steps in practice (where K is the size of the augmented diagrams, cf. Sec. 2.2). As discussed in Sec. 4.3, each barycenter is computed in parallel. The evaluations of the gradient of the weight energy (Sec. 3.2) and the atom energy (Sec. 3.3) both require O(NmK) steps. As described in Secs. 3.2 and 3.3, both evaluations can be run in parallel.

Tab. 1 evaluates the practical time performance of our multiscale algorithm for the optimization of the Wasserstein dictionary. In sequential, the runtime is roughly a function of the number of input diagrams (N) as well as their average size (|X|). The parallelization of our algorithm (with 20 cores) induces a significant speedup (up to 18 for the largest ensembles), resulting in an average computation time below 5 minutes, which we consider to be an acceptable pre-processing time, prior to interactive exploration. In comparison to the principal geodesic analysis of persistence diagrams (Tab. 1 of [START_REF] Pont | Principal geodesic analysis of merge trees (and persistence diagrams)[END_REF]), on a per ensemble basis, our approach is 1.56 times faster on average (on the same hardware).

Framework quality

Tab. 2 reports compression factors and average relative reconstruction errors for our application to data reduction (Sec. 5.1). For each ensemble, the compression factor is the ratio between the storage size of the input diagrams and that of the Wasserstein dictionary D * (the m atoms, of average size |a|, plus the N sets of barycentric weights). The relative reconstruction error is obtained by considering the Wasserstein distance between an input diagram and its barycentric approximation, divided by the maximum pairwise Wasserstein distance observed in the input ensemble. Then this relative reconstruction error is averaged over all the diagrams of the ensemble. Tab. 2 compares a naive optimization (Sec. 3) to our multi-scale strategy (Sec. 4.2). Specifically, for a given ensemble, the same target compression factor was used for both approaches (by imposing the same upper boundary on the total size of the atoms, Sec. 5.1). Tab. 2 shows that our multi-scale strategy (Sec. 4.2) enables the optimization to progress towards better solutions, as assessed by the improvement in reconstruction error of 30% on average. In comparison to the principal geodesic analysis of persistence diagrams (Appendix D of [START_REF] Pont | Principal geodesic analysis of merge trees (and persistence diagrams)[END_REF]), for the same compression factors, the error induced by our approach is on average 1.79 times larger. However, our approach is simpler, more flexible (our optimization is not subject to restrictive constraints, such as geodesic orthogonality) and slightly faster (Sec. 6.1). Fig. 10 provides a visual comparison between the 2D layouts obtained with our approach on the Isabel ensemble and those obtained with two typical dimensionality reduction techniques, namely MDS [START_REF] Kruskal | Multidimensional Scaling[END_REF] and tSNE [START_REF] Van Der Maaten | Visualizing Data Using t-SNE[END_REF], directly applied on the distance matrix obtained by computing the Wasserstein distance between all the pairs of input diagrams. For a given technique, to quantify Fig. 10. Comparison between the 2D layouts obtained with our approach (W2-Dict) and these obtained with typical dimensionality reduction approaches (W2-MDS [START_REF] Kruskal | Multidimensional Scaling[END_REF], W2-tSNE [START_REF] Van Der Maaten | Visualizing Data Using t-SNE[END_REF]) on the Isabel ensemble (all persistence pairs are considered). Here, the three approaches preserve well the clusters of the ensemble (NMI/ARI). As expected, W2-MDS provides (by design) the best metric preservation (SIM, bold). Our approach constitutes a trade-off between W2-MDS and W2-tSNE. Fig. 11. Comparison between the 2D layouts obtained with our approach (W2-Dict) and these obtained with typical dimensionality reduction approaches (W2-MDS [START_REF] Kruskal | Multidimensional Scaling[END_REF], W2-tSNE [START_REF] Van Der Maaten | Visualizing Data Using t-SNE[END_REF]) on a challenging ensemble. In this example (Asteroid Impact), the presence of an outlier (time step of the actual impact, red entry in the distance matrix, left) challenges cluster preservation. While W2-tSNE provides the best cluster preservation scores (NMI/ARI), it fails at visually depicting the outlier (red circle) as being far away from the other ensemble members. In contrast, W2-MDS and W2-Dict do a better job at isolating this outlier (red circle), with W2-Dict providing slightly improved cluster preservation scores (NMI/ARI).

TABLE 3

Detailed layout quality scores (i.e. bold: best values). On average (bottom row), our approach (W2-Dict) provides a trade-off between W2-MDS and W2-tSNE: it preserves the clusters (NMI/ARI) slightly better than W2-MDS and the metric (SIM) clearly better than W2-tSNE. its ability to preserve the structure of the ensemble, we run kmeans in the 2D layouts and evaluate the quality of the resulting clustering (given the ground-truth [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF]) with the normalized mutual information (NMI) and adjusted rand index (ARI). To quantify its ability to preserve the geometry of the ensemble, we report the metric similarity indicator SIM [START_REF] Pont | Principal geodesic analysis of merge trees (and persistence diagrams)[END_REF], which evaluates the preservation of the Wasserstein metric in the 2D layout. All these scores vary between 0 and 1, with 1 being optimal. In Fig. 10, the three approaches preserve well the clusters of the ensemble (NMI/ARI) and our approach provides a trade-off between MDS and tSNE in terms of metric preservation (SIM). Fig. 11 provides another visual comparison on a challenging ensemble (Asteroid Impact). There, the presence of an outlier (time step of the actual impact) challenges cluster preservation. While tSNE provides the best cluster preservation (NMI/ARI), it fails at visually depicting the outlier (red circle) as being far away from the other ensemble members. In contrast, MDS and our approach do isolate this outlier (red circle), with our approach providing slightly improved cluster preservation (NMI/ARI) over MDS. This illustrates the viability of our dimensionality reductions for outlier detection. Appendix C extends this visual analysis to all our test ensembles.

NMI ARI SIM Dataset W2-MDS W2-tSNE W2-Dict W2-MDS W2-tSNE W2-Dict W2-MDS W2-tSNE W2-
Tab. 3 extends our quantitative analysis to all our ensembles. MDS preserves well the metric (high SIM), at the expense of mixing ground-truth classes (low NMI/ARI). tSNE behaves symmetrically (higher NMI/ARI, lower SIM). Our approach provides a trade-off between the extreme behaviors of MDS and tSNE, with a cluster preservation slightly improved over MDS (NMI/ARI), and a clearly improved metric preservation over tSNE (SIM). Fig. 12 reports the evolution of the normalized energy E D along the optimization for all test ensembles, for the naive optimization strategy (Sec. 3), by using a number of atoms equal to the number of ground-truth classes (cf. our application to data reduction, Sec. 5.1). In this figure, the energy is normalized on a per ensemble basis, based on its initial value. This figure shows that the energy does decrease for most ensembles, but still with large oscillations due to the non-convex nature of the dictionary energy E D . In contrast, the energy evolution with our multi-scale strategy (Fig. 13) results in much less oscillations, which indicates the ability of this strategy to help the optimization explore in a more stable manner the locally convex areas of the energy (Appendix D discusses a counterexample). Specifically, in Fig. 13, one can observe sequences of discontinuous decrease patterns, characterized by an abrupt drop followed by a plateau. Each of these patterns corresponds to one persistence scale of our multi-scale strategy (this is particularly apparent on the Cloud Processes ensemble).

Fig. 14 provides a closer comparison between the two strategies on a selection of four ensembles. The Cloud Processes ensemble is an example where the naive optimization reaches a solution of slightly lower energy. For the other ensembles, our multiscale strategy leads to solutions of much lower energy, visually confirming the conclusions of Tab. 2. In this figure, one can also observe the characteristic decrease patterns discussed above, particularly apparent on the Ionization Front (3D) ensemble, which correspond to the distinct scales of our multi-scale strategy.

Limitations

Similarly to other optimization problems based on topological descriptors [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF], [START_REF] Pont | Principal geodesic analysis of merge trees (and persistence diagrams)[END_REF], [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF], [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF], our energy is not convex. Additionally, as shown in Fig. 12, the interleaving of the weight optimization (Sec. 3.2) with atom optimization (Sec. 3.3) can even lead to oscillations in the energy. As discussed in Sec. 6.2, our multi-scale strategy (Sec. 4.2) greatly mitigates both issues, with a more stable optimization than a naive approach (Sec. 3), which leads to relevant solutions which are exploitable in the applications (Sec. 5). However, we have found one example in our test ensembles (the Sea Surface Height ensemble), where our multi-scale strategy reached solutions which were arguably worse than obtained with a naive solution, as described in details in Fig. 15. In this example, the most persistent features in the diagrams are not particularly discriminative for the separation of the ground-truth classes. On the contrary, the variations between these classes seem mostly encoded by the low persistence features: in Fig. 15 clear separations in the distance matrices between the ground-truth classes only start to occur in the latest persistence scales (dashed sub-matrices, top right inset). This counter-intuitive observation goes against the rule of thumb traditionally used in topological data analysis, which states that the most persistent pairs encode the most important features in the data. For this example, when applying our framework to dimensionality reduction, the nondiscriminative aspect of the early persistence scales eventually lead our multi-scale strategy towards a local minium which does not separate the ground-truth classes well (planar layout, bottom left) in comparison to the naive strategy (planar layout, bottom right). Thus, for this ensemble, we reported dimensionality reduction results (Tab. 3, Appendix C) obtained with the naive optimization. In general, this means that when users are confronted with ensembles where the most persistent pairs are not the most responsible for data variability (hence class separation), the naive optimization may need to be considered additionally as it might provide solutions which better encode the ground-truth classes. Finally, as detailed in Appendix D, the presence of clear outliers can also challenge our optimization, especially when the selected number of atoms equals the number of ground-truth classes. Then, in this case, the best dictionary encoding will consequently be obtained by increasing the number of atoms, specifically, by considering that each outlier forms a singleton class.

CONCLUSION

In this paper, we presented an approach for the encoding of linear relations between persistence diagrams, given the Wasserstein metric. Specifically, we introduced a dictionary based representation of an ensemble of persistence diagrams, inspired by previous work on histograms [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning[END_REF]. We first documented a naive optimization, which interleaves the optimization of the barycentric weights of the input diagrams with the optimization of the atoms of the dictionary (Sec. 3). Then, we presented a multi-scale strategy (Sec. 4.2) leading to more stable optimizations and relevant solutions (Sec. 6.2). We demonstrated the utility of our contributions in applications (Sec. 5) to data reduction and dimensionality reduction, where the visualizations generated by our framework enable the visual identification of the main trends in the ensembles (Figs. 8,9), and the quick identification of outliers (Fig. 11). In contrast to previous work on persistence diagram encoding [START_REF] Pont | Principal geodesic analysis of merge trees (and persistence diagrams)[END_REF], our framework is simpler, less constrained and slightly faster in practice.

A natural direction for future work is the extension of our framework to other topological descriptors such as Reeb graphs or Morse-Smale complexes. However, this requires the definition of key geometrical tools, such as geodesic or barycenter computation algorithms, which is still an active research problem. We believe our framework for the dictionary encoding of persistence diagrams is an interesting practical step for the analysis of collections of persistence diagrams. In the future, we will continue our investigation of the adaptation of tools from optimal transport to the analysis of ensembles of topological descriptors, as we believe it can become a key solution in the long term for the advanced analysis of large-scale ensembles.

APPENDIX A ∇E W IS L-LIPSCHITZ

Proposition 1 : Let X be a persistence diagram and D = (a 1 , . . . , a m ) a Wasserstein dictionary of persistence diagrams. If the optimal matchings are constant, then E

W (λ λ λ ) = W 2 Y (λ λ λ ), X is convex and ∇E W is L-Lipschitz on Σ m .
Proof. Let λ λ λ = (λ 1 , . . . , λ m ) ∈ Σ m , Y (λ λ λ ) = y 1 (λ λ λ ), . . . , y K (λ λ λ ) the barycenter computed and φ λ λ λ ,1 , . . . , φ λ λ λ ,m the matchings between Y (λ λ λ ) and each atom (a 1 , . . . , a m ):

∀ j ∈ {1, . . . , K}, y j (λ λ λ ) = m ∑ i=1 λ i a φ λ λ λ ,i ( j) i . ( 12
)
We suppose the optimal matchings to be constant, thus we write φ i = φ λ λ λ ,i . Like before we consider the following gradient:

∇y j (λ λ λ ) = a φ 1 ( j) 1 . . . a φ m ( j) m . ( 13
)
Now recall the following expression for:

W 2 Y (λ λ λ ), X = min ψ λ λ λ :i X bi j -→i X K ∑ j=1 ∥y j (λ λ λ ) -x ψ λ λ λ ( j) ∥ 2 . ( 14
)
This minimum is always attained, and with the hypothesis on the optimal matchings we write ψ = ψ λ λ λ . Thus we rewrite:

W 2 Y (λ λ λ ), X = K ∑ j=1 ∥y j (λ λ λ )-x ψ( j) ∥ 2 = K ∑ j=1 ∥ m ∑ i=1 λ i (a φ i ( j) i -x ψ( j) )∥ 2 .
(15) W 2 Y (λ λ λ ), X is convex with λ λ λ and the gradient follows naturally: For the following part we denote H j = a φ 1 ( j) 1 x ψ( j) . . . a φ m ( j) m x ψ( j) . The Hessian then writes as

H = H(λ λ λ ) = 2 K ∑ j=1
(H j ) T H j . This shows that λ λ λ → W 2 Y (λ λ λ ), X is convex. Indeed for u ∈ R m we have:

u T Hu = 2 K ∑ j=1 u T (H j ) T (H j )u = 2 K ∑ j=1 ∥H j u∥ 2 ≥ 0 ( 17
)
This also shows that ∇E W is L-Lipschitz with L = ∥H∥. For numerical reasons, we bound L as follows:

L = ∥H∥ = 2 K ∑ j=1 (H j ) T (H j ) ≤ 2 K ∑ j=1 ∥(H j ) T (H j )∥ = 2 K ∑ j=1 ∥H j ∥ 2 .
(18) Thus for our algorithm, we consider the following gradient step:

ρ ≤ 2 K ∑ j=1 ∥H j ∥ 2 -1 . ( 19
)
APPENDIX B ∇e A IS L-LIPSCHITZ Proposition 2 : Let X be a persistence diagrams and λ λ λ = (λ 1 , . . . , λ m ) ∈ Σ m . If the optimal matchings are constant, the functions e A are convex and ∇e A is L-Lipschitz.

Proof. Let U = (u 1 , . . . , u m ) ∈ (R 2 ) m , for j ∈ {1, . . . , K} we have:

e A (U) = m ∑ i=1 λ i (u i -x ψ( j) ) 2 (20)
The gradient follows naturally:

∇e A (U) = 2    λ 1 . . . λ m    u i -x ψ( j) T (21)
Immediately we have the Hessian H j = H g j (U) = 2λ λ λ λ λ λ T , giving us the convexity of e A and the L-Lipschitzianity of ∇e A with L = ∥H j ∥ ≤ 2∥λ λ λ ∥ 2 ≤ 2m. For numerical reasons, we consider the larger upper bound: L ≤ 4m. Thus for our algorithm, we consider the following gradient step ρ ≤ (4m) -1 .

APPENDIX C DIMENSIONALITY REDUCTION

Fig. 16 extends Figure 10 (main manuscript) to all our test ensembles and it confirms visually the conclusions of the table of quality scores (Table 3 of the main manuscript).

APPENDIX D VOLCANIC ERUPTION ENSEMBLE

This appendix discusses the special case of the Volcanic eruption ensemble (12 members), for which a consistent energy increase can be observed in the Figure 13 of the main manuscript (normalized energy of our multi-scale optimization as a function of computation time), beyond 70% of the completion time (the optimization reaches the stopping conditions at 100%). The ground-truth classification of this ensemble contains 3 classes [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF]. However, one of these classes contains a clear outlier (light purple entry in the bottom views of Fig. 17), corresponding to a peak of activity in the eruption (see the terrain views of 4 members, bottom left of Fig. 17, including the outlier, light purple frame). The corresponding persistence diagram (light purple diagram in the aggregated birth/death space, bottom middle of Fig. 17) contains features which are significantly more persistent than the other diagrams (taken from distinct ground-truth classes, one color per class). Then, this outlier exhibits an excessively high distance to the rest of the ensemble, as illustrated in the Wassertein distance matrix (bottom right of Fig. 17, light purple entry).

The presence of this outlier challenges our optimization when using a number of atoms equal to the number of ground-truth classes (which is the default strategy documented in the main manuscript). As shown in the energy plots (Fig. 17, top), a consistent energy increase can be observed when using only 3 atoms (1 per ground-truth class, black curve). When removing the outlier, the energy evolution exhibits a more characteristic oscillating behavior (green curve). Finally, when initializing the optimization with 4 atoms (1 per class, plus 1 for the outlier), the Fig. 16. Comparison of the planar layouts for typical dimensionality reduction techniques all our test ensembles. The color encodes the classification ground-truth [START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF]. For each quality score, the best value appears bold. For the Sea Surface Height ensemble, the naive optimization procedure has been used (cf. Sec. 6.3 of the main manuscript).

optimization results in few oscillations and a consistent energy decrease (yellow curve). This indicates that the outlier member (light purple) should be interpreted as a singleton class and that the best dictionary encoding will consequently be obtained with 4 atoms.

Fig. 1 .

 1 Fig. 1. Persistence diagrams of a clean (left) and noisy (right) terrain (dark blue spheres: minima, dark yellow: maxima, other: saddles). The three main hills are clearly represented with long bars in the persistence diagrams. In the noisy persistence diagram, small bars encode noise.

Fig. 2 .

 2 Fig. 2. Optimal matching (green dashes, right) with regard to W between the two persistence diagrams (center) of two terrains (left).

Fig. 3 .

 3 Fig. 3. Wasserstein barycenter (cyan, uniform weights) of 3 persistence diagrams (center) of 3 terrains (left). Each barycenter point (cyan sphere) is the barycenter of its matched points in the inputs (cyan triangle).

  2.3), Y (λ λ λ ) and the atoms contain ∑ m i=1 |a i |+|Y (λ λ λ )| points each, where |a i | and |Y (λ λ λ )| denote the number of non-diagonal points in a i and Y (λ λ λ ) respectively.

Fig. 4 .

 4 Fig. 4. Optimizing the weights of the barycenter Y (λ λ λ ) (cyan diagram) to improve its approximation of X (grey diagram), given a fixed Wasserstein dictionary D of 3 atoms (dark blue, yellow, green). At a given iteration t (center), a step ρ λ λ λ is made along the gradient of the weight energy E W (cyan arrows), resulting in an improved estimation at iteration t + 1 (right).

Fig. 5 .

 5 Fig. 5. Optimizing the atoms of the Wasserstein dictionary D (dark blue, yellow and green diagrams). At a given iteration t (center), a step ρ D is made along the gradient of the pointwise atom energy e A (arrows on each triangle), resulting in a dictionary (right) that enables an improved barycentric approximation (Y (D), cyan) of the input diagram X (grey).

1 )

 1 Computing the Wasserstein barycenter Y (D t ) (Sec. 2.3); 2) Computing the Wasserstein distance W 2 Y (D t ), X (Eq. 3); 3) For each point y j (D) of Y (D): a) Estimating the gradient ∇e A (D j ) (Eq. 10); b) Applying one step ρ D of gradient descent on D j (Eq. 11); 4) Remove the K least persistent points from each atom a i .

Fig. 7 .

 7 Fig. 7. Multi-resolution representation of an input persistence diagram (taken from the Isabel ensemble). At a given resolution (from left to right), only the persistence pairs above a given persistence threshold (red dash line) are considered in the optimization.

Algorithm 1 :

 1 Multi-scale Wasserstein Dictionary Optimization. Input: Set of persistence diagrams {X 1 , . . . , X N }; Output 1: Wasserstein Dictionary D * ; Output 2: Barycentric weights λ λ λ * 1 , . . . , λ λ λ * N ; for τ ∈ {τ 0 , . . . , τ r } do if τ == τ 0 then Initialization (Sec. 4.1); end while E D (Eq. 2) decreases do for n ∈ {1, . . . , N} do Perform a gradient step ρ λ λ λ n along ∇E W relative to X n (Sec. 3.2); end for n ∈ {1, . . . , N} do Perform a gradient step ρ D along ∇E A relative to X n (Sec. 3.3); end end end

Fig. 12 .

 12 Fig. 12. Evolution of the (normalized) energy E D along the optimization, with a naive optimization (Sec. 3), for all our test ensembles.

Fig. 13 .

 13 Fig. 13. Evolution of the (normalized) energy E D along the optimization, with our multi-scale strategy (Sec. 4.2), for all our test ensembles.

Fig. 14 .

 14 Fig.14. Comparison of the evolutions of the (normalized) energy E D between the naive optimization (Sec. 3, N, dashed curves) and our multiscale strategy (Sec. 4.2, MS, solid curves) for four ensembles. For this experiment, the energy has been normalized with regard to the initial energy of the naive optimization. The Cloud Processes ensemble is an example where the naive optimization reaches a solution of slightly lower energy. For the other three ensembles, our multi-scale strategy leads to solutions of much lower energy, through a sequence of characteristic, discontinuous decrease patterns (abrupt drop followed by a plateau) corresponding to the five persistence scales of our multi-scale strategy.

Fig. 15 .

 15 Fig.15. Counter-example for our multi-scale strategy (Sea Surface Height ensemble). Top: Wasserstein distance matrices for the first four persistence scales of our multi-scale strategy. The ground-truth classes only start to become visible in the distance matrix between the third and fourth scale (dashed sub-matrices in the fourth scale). As a result, our multi-scale strategy is attracted in the first scales towards a local minimum of the energy which does not encode well the ground-truth classes (dimensionality reduction, bottom left). In contrast, the naive optimization manages to reach a solution which separates well the ground-truth classes (dimensionality reduction, bottom right).

∇W 2 Y 1 -

 21 (λ λ λ ), X = 2 x ψ( j) ) T . . .

  (λ λ λ )x ψ( j) .[START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF]

Fig. 17 .

 17 Fig.17. Evolution of the (normalized) energy E D along the optimization (top curves), with our multi-scale strategy, for the Volcanic eruption ensemble, for distinct initializations. The ground-truth classification of this ensemble contains 3 classes[START_REF] Pont | Wasserstein Distances, Geodesics and Barycenters of Merge Trees[END_REF], including one outlier (light purple entry in the bottom views, from left to right: terrain view of the data, aggregated birth/death space, distance matrix). A clear energy increase can be observed when considering the entire ensemble (black curve), while a more characteristic oscillating behavior occurs when discarding the outlier (green curve). When initializing the optimization with 4 atoms (1 per class, plus 1 for the outlier), the optimization results in few oscillations and a consistent energy decrease (yellow curve).

TABLE 1

 1 Running times (in seconds) of our multi-scale algorithm (1 and 20 cores).

	Dataset	N	|X|	1 core	20 cores	Speedup
	Asteroid Impact (3D)	20	220	259	35	7.50
	Dark matter (3D)	40	216	1,323	188	7.04
	Earthquake (3D)	12	97	113	92	1.23
	Ionization front (3D)	16	757	4,230	595	7.11
	Isabel (3D)	12	1,310	1,609	270	5.96
	Viscous Fingering (3D)	15	158	252	49	5.14
	Cloud processes (2D)	12	1,176	914	64	14.28
	Ionization front (2D)	16	186	145	45	3.22
	Sea surface height (2D)	48	1,567	14,587	792	18.42
	Starting vortex (2D)	12	125	140	24	5.83
	Vortex street (2D)	45	43	1,061	241	4.40
	Volcanic eruptions (2D)	12	860	2,798	706	3.96

TABLE 2

 2 Comparison of the average relative reconstruction error (for a common target compression factor), between a naive optimization (Sec. 3) and our multi-scale strategy (Sec. 4.2). Our multi-scale algorithm improves the error by 30% on average over the naive approach.

	Dataset	N	|X|	m	|a|	Factor	Error (Naive)	Error (Multi-Scale)
	Asteroid Impact (3D)	20	220	4	493	2.20	0.09	0.06
	Dark matter (3D)	40	216	4	215	10.87	0.15	0.12
	Earthquake (3D)	12	98	3	120	3.05	0.16	0.04
	Ionization front (3D)	16	757	4	1,044	2.90	0.29	0.20
	Isabel (3D)	12	1,310	3	1,049	5.49	0.34	0.37
	Viscous Fingering (3D)	15	158	3	41	2.78	0.15	0.11
	Cloud processes (2D)	12	1,176	3	381	5.97	0.38	0.41
	Ionization front (2D)	16	186	4	300	2.68	0.38	0.17
	Sea surface height (2D)	48	1,567	4	534	20.98	0.54	0.61
	Starting vortex (2D)	12	125	2	379	1.98	0.22	0.09
	Vortex street (2D)	45	43	5	75	5.08	0.18	0.04
	Volcanic eruptions (2D)	12	860	3	345	9.97	0.20	0.20

ACKNOWLEDGMENTS

This work is partially supported by the European Commission grant ERC-2019-COG "TORI" (ref. 863464, https://erc-tori.github.io/).

Keanu Sisouk is a Ph.D. student at Sorbonne University. He received his master degree in Mathematics from Sorbonne University in 2021. His fields of interests lie on topological methods for data analysis, optimal transport, optimization methods, statistics and partial derivative equations.

Julie