Fast Approximation of Persistence Diagrams with Guarantees

Jules Vidal and Julien Tierny * CNRS -Sorbonne Universit é

Figure 1: Approximations of persistence diagrams for a CT scan of a backpack, with different Bottleneck approximation errors. High persistence features correspond to high density objects present in the bag (leftmost, bottom). In this example, our controlled approximation reduces computation time by 70% for an error tolerance of 10%. Our algorithm also provides an approximation of the scalar field that is precise around persistent features (top, volume rendering of the approximated fields) and deteriorated elsewhere (bottom, isocontours capturing the cloth of the bag). For each approximation, the thirty most persistent maxima are represented by spheres (top views) which correctly capture the features of the data. The uncertainty resulting from our approximation is visualized in the diagram with (i) colored squares, which bound the correct location of certain features (which are guaranteed to be present in the exact result). These are located outside a (ii) red band (in the vicinity of the diagonal), which denotes features which might not exist in the exact diagram. Our approximations precisely capture the high persistence features of the data (out of the red band, black numbers) and collect significantly less noisy features (red numbers, red band) than the exact result (numbers in parentheses).

ABSTRACT

This paper presents an algorithm for the efficient approximation of the saddle-extremum persistence diagram of a scalar field. Vidal et al. introduced recently a fast algorithm for such an approximation (by interrupting a progressive computation framework [START_REF] Vidal | A progressive approach to scalar field topology[END_REF]). However, no theoretical guarantee was provided regarding its approximation quality. In this work, we revisit the progressive framework of Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] and we introduce in contrast a novel approximation algorithm, with a user controlled approximation error, specifically, on the Bottleneck distance to the exact solution. Our approach is based on a hierarchical representation of the input data, and relies on local simplifications of the scalar field to accelerate the computation, while maintaining a controlled bound on the output error. The locality of our approach enables further speedups thanks to shared memory parallelism. Experiments conducted on real life datasets show that for a mild error tolerance (5% relative Bottleneck distance), our approach improves runtime performance by 18% on average (and up to 48% on large, noisy datasets) in comparison to standard, exact, publicly available implementations. In addition to the strong guarantees on its approximation error, we show that our algorithm also provides in practice outputs which are on average 5 times more accurate (in terms of the L 2 -Wasserstein distance) than a naive approximation baseline. We illustrate the utility of our approach for interactive data exploration and we document visualization strategies for conveying the uncertainty related to our approximations.

INTRODUCTION

Modern datasets, acquired or simulated, are continuously gaining in geometrical complexity, thanks to the ever-increasing accuracy of acquisition devices or computing power of high performance systems. This geometrical complexity makes interactive exploration and analysis difficult, which challenges the interpretation of large datasets by end users. This motivates the definition of expressive data abstractions, capable of capturing the main features present in large datasets into concise representations, which visually convey the most important information to the users. In that context, Topological Data Analysis (TDA) [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] forms a family of generic, robust, and efficient techniques whose utility has been demonstrated in a number of visualization tasks [START_REF] Heine | A survey of topology-based methods in visualization[END_REF] for revealing the implicit structural patterns present in complex datasets. Examples of popular application fields 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], 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 time-dependent vortex regions based on the acceleration magnitude[END_REF], bioimaging [START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF][START_REF] Carr | Simplifying Flexible Isosurfaces Using Local Geometric Measures[END_REF], quantum 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] 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]. These applications rely on established topological data abstractions, such as contour trees [START_REF] Boyell | Hybrid techniques for real-time radar simulation[END_REF][START_REF] Carr | Computing contour trees in all dimensions[END_REF][START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF][START_REF] Smirnov | Triplet Merge Trees[END_REF], Reeb graphs [START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF][START_REF] Pascucci | Robust on-line computation of Reeb graphs: simplicity and speed[END_REF][START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF] or Morse-Smale complexes [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF][START_REF] Gyulassy | Shared-Memory Parallel Computation of Morse-Smale Complexes with Improved Accuracy[END_REF][START_REF] Gyulassy | A practical approach to Morse-Smale complex computation: Scalability and generality[END_REF][START_REF] Robins | Theory and Algorithms for Constructing Discrete Morse Complexes from Grayscale Digital Images[END_REF]. In particular, the Persistence diagram [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] is a concise data representation, which visually summarizes the population of features of interest present in large datasets, as a function of a measure of importance called topological persistence. Its conciseness made it increasingly popular in data visualization and analysis, where it quickly provides visual hints regarding the number and salience of the structural features present in large datasets.

While topological methods usually have an acceptable time complexity, the construction of the above data representations can still be time consuming in practice for large, real-life datasets. Thus, when they are integrated into larger analysis pipelines, TDA algorithms can become a serious time bottleneck. Several research directions are possible to improve the time performance of TDA algorithms. A natural avenue for performance improvement consists in revisiting these algorithms for a parallel computation [START_REF] Carr | Parallel peak pruning for scalable SMP contour tree computation[END_REF][START_REF] Gueunet | Contour forests: Fast multi-threaded augmented contour trees[END_REF][START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF][START_REF] Gueunet | Task-based Augmented Reeb Graphs with Dynamic ST-Trees[END_REF][START_REF] Gyulassy | Shared-Memory Parallel Computation of Morse-Smale Complexes with Improved Accuracy[END_REF][START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF][START_REF] Shivashankar | Parallel Computation of 3D Morse-Smale Complexes[END_REF], which is particularly relevant for high-performance hardware. Another direction consists in considering degraded computations, which trade accuracy for speed. This direction is particularly appealing for persistence diagrams, as they can contain, for noisy datasets, many features with low persistence which would gain in being only approximated, as they often have in practice only little relevance (important features typically have a high persistence, see Fig. 1).

In this work, we introduce an algorithm for the fast approximation of extremum persistence diagrams of large scalar datasets, with a user-controlled error bound. In particular, we re-visit the progressive computation framework of Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] and derive a multiresolution strategy for accelerating the computation, based on a controlled degradation of the input data (hence resulting in a controlled approximation of the result). Extensive experiments demonstrate the time performance gain provided by our algorithm (up to 48% on large, noisy datasets) for mild error tolerance (e.g. 5%), as well as its superior accuracy in comparison to a naive approximation baseline. We also document several visualization strategies for our approximated topological features, to easily enable the visual identification of certain features (which will be present in the exact diagram) along with visual clues of their positional uncertainty in the diagram.

Related Work

The literature related to our work can be classified into two main groups: (i) topological methods, and (ii) multiresolution methods. (i) Topological methods have gained a growing interest from the visualization community for the last two decades [START_REF] Heine | A survey of topology-based methods in visualization[END_REF]. Specifically, our work is closely related to the construction of topological data abstractions on scalar data, that provide a generic and robust description of the features in the data. When considering discretized data as piecewise-linear (PL) scalar fields defined on PL-manifolds, many concepts from the original Morse theory [START_REF] Milnor | Morse Theory[END_REF] can be adapted to the computational setting. As such, combinatorial algorithms have been developed to efficiently compute topological abstractions in the discrete setting. A local, combinatorial characterization of the critical points of a PL scalar field was introduced by Banchoff [START_REF] Banchoff | Critical points and curvature for embedded polyhedral surfaces[END_REF]. Critical points often correspond to features of interest in the data, but can be numerous in the presence of noise. Topological persistence [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF][START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] addresses this issue by introducing an importance measure on critical points. It can be visualized in the Persistence diagram (Sec. 2.3), a visual barcode that can be generally computed using matrix reduction operations [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF][START_REF] Edelsbrunner | Topological persistence and simplification[END_REF]. The Persistence diagram has been shown to be stable [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF][START_REF] Cohen-Steiner | Lipschitz functions have lp-stable persistence[END_REF] under well-established metrics such as the Bottleneck [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] and Wasserstein [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF] distances. However, some applications call for more discriminative topological abstractions. Merge trees and contour trees, for instance, track the merge and split events of connected components of sub-level sets and level sets that happen at the critical points. Combinatorial algorithms were first designed to compute these trees in low dimension [START_REF] Boyell | Hybrid techniques for real-time radar simulation[END_REF]. An efficient algorithm was later introduced by Carr et al. [START_REF] Carr | Computing contour trees in all dimensions[END_REF] for the computation of the contour tree in all dimensions, with optimal time-complexity. A lot of work also focused on the parallel computation of these structures [START_REF] Carr | Parallel peak pruning for scalable SMP contour tree computation[END_REF][START_REF] Gueunet | Task-based Augmented Merge Trees with Fibonacci Heaps[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] Morozov | Distributed contour trees[END_REF][START_REF] Pascucci | Parallel Computation of the Topology of Level Sets[END_REF][START_REF] Smirnov | Triplet Merge Trees[END_REF]. The Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF] is the generalization of contour trees to domains that are not simply connected, which motivates its use for shape analysis applications [START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF]. In order to track the connected components of level sets on such domains, original algorithms resorted to slicing strategies [START_REF] Biasotti | Extended Reeb Graphs for Surface Understanding and Description[END_REF]. Later work optimized the slicing approach to restrict it to iso-contours at critical points [START_REF] Doraiswamy | Output-Sensitive Construction of Reeb Graphs[END_REF][START_REF] Patanè | Reeb graph computation based on a minimal contouring[END_REF]. Algorithms with optimal time complexity were introduced for the 2D case [START_REF] Cole-Mclaughlin | Loops in Reeb graphs of 2-manifolds[END_REF] and in arbitrary dimension [START_REF] Parsa | A deterministic o(m log m) time algorithm for the reeb graph[END_REF]. Other work focused on its parallel computation [START_REF] Gueunet | Task-based Augmented Reeb Graphs with Dynamic ST-Trees[END_REF] or presented efficient algorithms for its computa-tion in specific settings [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]. An other widely used topological abstraction is the Morse-Smale complex [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF], which encodes the relations between critical points in terms of unique integral lines of the gradient field. These integral lines segment the domain into cells in which the gradient integrates to identical extremities. Algorithms for its computation were designed for PL manifolds [START_REF] Edelsbrunner | Morse-smale complexes for piecewise linear 3-manifolds[END_REF], or later [START_REF] Gyulassy | A practical approach to Morse-Smale complex computation: Scalability and generality[END_REF][START_REF] Robins | Theory and Algorithms for Constructing Discrete Morse Complexes from Grayscale Digital Images[END_REF] for the Discrete Morse Theory [START_REF] Forman | A User's Guide to Discrete Morse Theory[END_REF]. Parallel methods have been documented [START_REF] Gyulassy | Shared-Memory Parallel Computation of Morse-Smale Complexes with Improved Accuracy[END_REF][START_REF] Shivashankar | Parallel Computation of 3D Morse-Smale Complexes[END_REF]. For the case of multivariate scalar data, recent work [START_REF] Carr | Joint Contour Nets[END_REF][START_REF] Tierny | Jacobi Fiber Surfaces for Bivariate Reeb Space Computation[END_REF] investigated efficient algorithms for the Reeb space [START_REF] Edelsbrunner | Reeb spaces of piecewise linear mappings[END_REF] computation, a generalization of the Reeb graph. (ii) Multiresolution hierarchical representations have been largely used by the visualization community to compute and store data representations at different levels of details [START_REF] Weiss | Diamond hierarchies of arbitrary dimension[END_REF][START_REF] Weiss | Supercubes: A high-level primitive for diamond hierarchies[END_REF]. They have been notably applied to the efficient and robust extraction of isosurfaces [START_REF] Gerstner | Topology preserving and controlled topology simplifying multiresolution isosurface extraction[END_REF][START_REF] Gregorski | Interactive view-dependent rendering of large isosurfaces[END_REF][START_REF] Pascucci | Time critical isosurface refinement and smoothing[END_REF] and isocontours [START_REF] Lewiner | Hierarchical isocontours extraction and compression[END_REF]. Such representations have been also documented for the contour tree [START_REF] Pascucci | Multi-resolution computation and presentation of contour trees[END_REF], the Reeb graph [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF] and the Morse-Smale complex [START_REF] Bremer | A Multi-Resolution Data Structure for 2-Dimensional Morse Functions[END_REF][START_REF] Guenther | Efficient computation of a hierarchy of discrete 3d gradient vector fields[END_REF][START_REF] Iuricich | Hierarchical forman triangulation: A multiscale model for scalar field analysis[END_REF]. These representations allow to adaptively simplify these topological abstractions in a fineto-coarse manner. As such, only the output data structures are processed hierarchically, while the input data are processed with classical algorithms on the finest level of the hierarchy.

More closely related to our work, Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] recently introduced novel ideas for the progressive computation of critical points and saddle-extremum persistence diagrams. In their approach, the input data is processed in a coarse-to-fine manner with the help of a multiresolution hierarchy on the model of the classical red subdivision [START_REF]Simplizialzerlegungen von beschrankter Flachheit[END_REF][START_REF] Bey | Tetrahedral grid refinement[END_REF][START_REF] Loop | Smooth Subdivision Surfaces Based on Triangles[END_REF][START_REF] Bank | Refinement algorithms and data structures for regular local mesh refinement[END_REF][START_REF] Zhang | Successive subdivision of tetrahedra and multigrid methods on tetrahedral meshes[END_REF] of triangulated regular grids [START_REF] Bey | Simplicial grid refinement: on Freudenthal's algorithm and the optimal number of congruence classes[END_REF][START_REF] Kuhn | Some combinatorial lemmas in topology[END_REF]. The output of the approach is progressively and efficiently updated from one level of the hierarchy to the next. In practice, their algorithms are interruptible, which means that the processing of the hierarchy can be stopped before reaching the finest level and still deliver an exploitable result in a lower amount of computation time. Although the authors demonstrated the quality of the interrupted outputs experimentally, there are no theoretical guarantees on their approximation quality (specifically, no error bound). In contrast, we present in this work a method for the approximate computation of a persistence diagram, with strong guarantees on the approximation error.

Contributions

This paper makes the following new contributions:

1. A fast approximation algorithm for extremum persistence diagrams, with controlled error: We introduce a fast algorithm for the approximated computation of persistence diagrams of extremum/saddle pairs. Our approach provides strong, usercontrolled guarantees on the Bottleneck distance between the approximated diagram and the exact result. By construction, our approximation method processes a smaller number of topological events than exact methods, which reduces its run time. We observe an 18% reduction of the run time on average, for a mild tolerance of 5% on the relative Bottleneck error. Moreover, our approach provides much more accurate approximations (in terms of the L 2 Wasserstein distance) than a naive approximation baseline. 2. An approximated visualization of topological features with uncertainty assessment: We describe several strategies for the visualization of our topological approximations. As a by-product of our algorithm, we provide an approximated scalar field, for visualization purposes, which exactly matches our approximated diagrams and in which the approximated critical points can be reliably embedded. We augment the approximated diagrams with visual glyphs assessing the uncertainty of the approximation. These glyphs enable the visual identification of certain features (which will be present in the exact diagram) along with visual clues on their positional uncertainty in the diagram. 3. A reference implementation (additional material): Finally, we provide a reference C++ implementation of our algorithms that can be used to replicate our results as well as for future benchmarks.

PRELIMINARIES

This section presents the theoretical background of our work. It contains definitions adapted from the Topology ToolKit [START_REF] Tierny | The Topology ToolKit[END_REF]. We refer the reader to textbooks [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] for an introduction to topology.

Input Data

The input data is a piecewise linear (PL) scalar field f : M → R defined on a PL d-manifold M, with d ≤ 3 in our applications. The field f is characterized by its values at the vertices of M, as scalar values are linearly interpolated on the simplices of M of higher dimension. f is enforced to be injective on the vertices of M, using a symbolic perturbation inspired by Simulation of Simplicity [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF].

Critical Points

We define a sub-level set of f as the pre-image of (-∞, w) by

f , noted f -1 -∞ (w) = {p ∈ M | f (p) < w}.
When continuously increasing w, the number of topological features of f -1 -∞ (w) (its connected components, independant cycles and voids) changes only at specific vertices, called the critical points of f [START_REF] Milnor | Morse Theory[END_REF].

A local characterization of critical points, based on their link, was introduced by Banchoff [START_REF] Banchoff | Critical points and curvature for embedded polyhedral surfaces[END_REF]. For two simplices τ, σ ∈ M, we note τ < σ if τ is a face of σ , i.e. if τ is defined by a non-empty, strict subset of the vertices of σ . We call star of a vertex v ∈ M the set of its co-faces, noted

St(v) = {σ ∈ M | v < σ }. The link of v, noted Lk(v) = {τ ∈ M | τ < σ , σ ∈ St(v) , τ ∩ v = /
0}, is defined as the set of all the faces τ of the simplices of St(v) that do not intersect v. Intuitively, the link of v is the boundary of a small simplicial neighborhood St(v) of v. The lower link Lk -(v) is the set of simplices of Lk(v) that are lower than v in terms of f values:

Lk -(v) = {σ ∈ Lk(v) | ∀v ∈ σ , f (v ) < f (v)}. Symmetrically, the upper link of v is given by Lk + (v) = {σ ∈ Lk(v) | ∀v ∈ σ , f (v ) > f (v)}. If both Lk -(v)
and Lk + (v) are non-empty and simply connected, v is a regular point. In any other case, v is a critical point, of criticality or critical index I(v): 0 for a local minimum (Lk

-(v) = / 0), d for a local maximum (Lk + (v) = / 0), or i for a i-saddle (0 < i < d).

Persistence Diagrams

As an isovalue w continuously increases, the topology of f -1 -∞ (w) evolves. For a regular point v, the sub-level sets enter St(v) through Lk -(v) and exit through Lk + (v) without any local change in the topology, as Lk -(v) and Lk + (v) both have exactly one connected component. Otherwise, if v is critical, a topological feature (a connected component, a cycle or a void) of f -1 -∞ (w) is either created or destroyed at the value f (v). For instance, the connected components of f -1 -∞ (w) are created at the local minima of f , said to be their birth points. When two topological features meet at a critical point, we use the Elder rule [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] to state that the younger one (the one created last, in terms of f values) dies at the benefit of the oldest. In particular, the connected components of f -1 -∞ (w) die at 1-saddles. Each topological feature of f -1 -∞ (w) is thus characterized by a unique pair (c 0 , c 1 ) of critical points, corresponding to its birth and death, with f (c 0 ) < f (c 1 ) and I(c 0 ) = I(c 1 ) -1. The topological persistence [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] of this pair, noted p(c 0 , c 1 ) = f (c 1 )f (c 0 ), denotes the lifetime of the corresponding feature in f -1 -∞ (w) in terms of its scalar range. The persistence of the connected components in f -1 -∞ (w) is encoded by minimum/1-saddle pairs. In 3D, (d -1)saddle/maximum pairs characterize the persistence of the voids of f -1 -∞ (w) while 1-saddle/2-saddle pairs encode the lifetime of its independent cycles. In practice, the features of interest in the data are often characterized by the extrema of the field f . Thus we will only focus on extremum/saddle pairs in the following.

The topological persistence of each critical pair gives a measure of importance on the corresponding extremum of the field, that has been shown to be reliable to distinguish between noise and important features. The Persistence diagram [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] of f , noted D( f ), is a visual representation of the ensemble of features in the data, where each persistence pair (c 0 , c 1 ) is embedded as a point in the 2D plane, at coordinates ( f (c 0 ), f (c 1 )). The persistence of each pair can thus be read in the diagram as the height of the point to the diagonal. Consequently, each topological feature of f -1 -∞ (w) can be visualized in the diagram as a bar (Fig. 2), whose height indicates its importance in the data. Large bars corresponding to high persistence features stand out visually, while low persistence pairs, likely to be associated with noisy features, are represented by small bars in the vicinity of the diagonal. The persistence diagram is a concise visual depiction of the repartition of features in the data and has been shown to be a stable [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] and useful tool for data summarization tasks. As seen in Fig. 2, it encodes the number, ranges and salience of features of interest, and gives hints about the level of noise in the data.

The Wasserstein distance is a well-established metric between persistence diagrams. Originally from the field of Transportation theory, it is defined as the cost of an optimal assignment between two diagrams D( f ) and D(g), where the cost of matching a point

a = (x a , y a ) ∈ D( f ) to a point b = (x b , y b ) ∈ D(b)
is given by a pointwise distance d q in the 2D birth/death space of the diagrams:

d q (a, b) = (|x b -x a | q + |y b -y a | q ) 1/q = a -b q (1)
By convention, d q (a, b) = 0 if both a and b are on the diagonal (x a = y a and x b = y b ). The L q -Wasserstein distance between D( f ) and D(g) is then given by:

W q D( f ), D(g) = min φ ∈Φ ∑ a∈D( f ) d q a, φ (a) q 1/q (2)
where Φ is the set of all possible maps φ mapping each point a ∈ D( f ) to a point b ∈ D(g) or to its projection onto the diagonal. When q converges to infinity, the Wasserstein distance converges to the Bottleneck distance, another practical metric that measures the worst mismatch between D( f ) and D(g):

W ∞ D( f ), D(g) = min φ ∈Φ max a∈D( f ) a -φ (a) ∞ (3)
The persistence diagram is stable under these two metrics, which intuitively means that a small variation in the scalar field entails a small difference in the resulting distances. In particular, the Bottleneck distance between two diagrams is bounded by the L ∞ distance between the corresponding scalar fields [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF]:

W ∞ D( f ), D(g) ≤ f -g ∞ (4)

Edge-Nested Hierarchical Data Representation

To design progressive algorithms for topological data analysis, Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] used a progressive representation of the input data. This representation is based on an edge-nested multiresolution hierarchy of the input domain M. The hierarchy H = {M 0 , M 1 , . . . , M h } is a growing sequence of PL manifolds, such that all the vertices present in one level M i of H are present in all the following levels M j , j > i. This naturally yields a hierarchy of PL scalar fields { f 0 , . . . , f h } such as each vertex maintains the same scalar value from its level of insertion i forward:

∀v ∈ M i , ∀ j | i ≤ j ≤ h, f j (v) = f i (v) = f (v).
The edge-nested property states that new vertices inserted at a level i of the hierarchy are exclusively inserted on existing edges (connecting old vertices) at level i -1, one per edge, and that new edges only connect new vertices (see Fig. 3). We refer the reader to Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] for a detailed formalization. Edge-nested hierarchies can be easily generated from regular grids. Given a d-dimensional regular grid G 0 , a valid PL d-manifold M h can be obtained from G 0 by considering its Kuhn triangulation [START_REF] Kuhn | Some combinatorial lemmas in topology[END_REF] (Fig. 3). A sequence {G 0 , G 1 , . . . , G h } of decimated grids can be recursively defined from G 0 , where G i is obtained for i > 0 by sub-sampling only vertices with even coordinates in G i-1 . The edge-nested triangulation hierarchy H = {M 0 , M 1 , . . . , M h } is then defined by considering the triangulation M i of each level G h-i . Further implementations details can be found in Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF].

Topological Invariant Vertices

The advantage of using an edge-nested triangulation hierarchy to process the data resides in the fact that topological information can be computed on the vertices of M i and remain valid throughout the rest of the hierarchy. In particular, as discussed by Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF], when progressing from a level of the hierarchy to the next, the topological structure of the link is invariant for all old vertices: their number of neighbor vertices remains the same, and the adjacency relations between neighbors are preserved. This important property (Fig. 4) enables the efficient identification of vertices that do no change criticality when progressing to the next hierarchical level.

For each vertex v of M i , its link polarity (blue +/glyphs in Fig. 4) indicates for each one of its neighbors n ∈ Lk(v) i , whether n ∈ Lk + (v) i (positive polarity) or n ∈ Lk -(v) i (negative polarity). As the structure of Lk(v) i is invariant, the criticality of v is characterized by the polarity of its link. In particular, if the link polarity of v stays unchanged when progressing a level down the hierarchy, the criticality of v remains the same and does not need to be recomputed. Given an edge (v 0 , v 1 ) of M i that gets subdivided into two old edges

(v 0 , v n ) and (v n , v 1 ) along the new vertex v n of M i+1 , v n is said to be monotonic if f (v n ) ∈ ( f (v 0 ), f (v 1 )), assuming that f (v 0 ) < f (v 1 ).
Otherwise, v n is called a non-monotonic vertex, and the polarity of the link of v 0 or v 1 is necessarily changed, as shown in Fig. 4b. On the opposite, an old vertex keeps the same link polarity (and thus the same criticality) if all its new neighbors are monotonic (Fig. 4c). It is then called a Topologically Invariant old vertex. Potentially, depending on the regularity of the data, a large part of the old vertices at level i can be topologically invariant. Similarly, the criticality of a part of the new vertices can be deduced from the hierarchy with the identification of non-monotonic vertices. It can be shown [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] that if a new vertex v of M i is monotonic and if all its new neighbors are also monotonic, then v is necessarily a regular point. It is called a Topologically Invariant new vertex.

It has been shown [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] that Topologically Invariant vertices (TIs) represent a large part of the data in practice (around 70% on average on a diverse sample of real-life data sets). Thus, the edge-nested hierarchical representation of the data enables a highly efficient detection of these TI vertices. It proved itself useful to carry out topological data analysis tasks such as the extraction of critical points or the computation of a persistence diagram, as significant shortcuts can be made in the computation for TI vertices.

OVERVIEW

Fig. 5 provides an overview of our approach, which revisits the progressive framework of Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF], to derive a fast approximation algorithm with strong guarantees. First, we exploit their multiresolution hierarchy of the input data (Sec. 2.4) to quickly update, down to the finest hierarchy level, the polarity of each vertex (a local information used to identify critical points). This step is described in Sec. 4.1. During the hierarchy traversal, in contrast to their original approach, we artificially increase the number of topologically invariant vertices (Sec. 2.5) in order to significantly speedup the computation, through a procedure called vertex folding, which artificially degrades the input data. This step is described in Sec. 4.2. The data degradation induced by the vertex folding procedure is precisely controlled in the process, to provide strong guarantees on the approximation error of the output. This is described in Sec. 4.3. Finally, we describe in Sec. 4.4 how to handle degenerate configurations such as flat plateaus. Overall, our new approach involves three major differences to the progressive framework of Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF], which are detailed in the rest of this section:

1. Our new approximation algorithm is not progressive: it does not generate a sequence of progressively refined outputs. Instead, our traversal of the multiresolution hierarchy only updates a minimal amount of information (the vertex polarity). The criticality of each vertex is only evaluated after the hierarchy traversal is finished (while the criticality is updated at each hierarchy level in [START_REF] Vidal | A progressive approach to scalar field topology[END_REF]).

2. Our approximation approach is based on a multiresolution degradation of the input data, which accelerates the overall computation, while maintaining a controlled output error.

3. Overall, in contrast to the progressively interrupted results of Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF], our approximated output is provided with strong guarantees on its approximation error (in terms of relative Bottleneck distance) and with a significantly improved practical accuracy (in terms of the L 2 Wasserstein distance).

TOPOLOGY APPROXIMATION

This section presents our novel approach for the controlled approximation of an extremum/saddle persistence diagram. In the following, we focus on the case of minimum/1-saddle pairs ((d -1)saddle/maximum pairs being treated symmetrically). Our method is based on folding operations for non-monotonic vertices. As the hierarchy is processed, non-monotonic vertices are inserted at each level. When a non-monotonic vertex is inserted in the hierarchy on a new edge, we can purposely decide to reinterpolate this vertex to enforce its monotony, and hence accelerate the computation of its criticality. The resulting error is the difference between the real scalar value at this vertex and the interpolated value, which bounds the Bottleneck error on the diagram estimation, as detailed next.

Hierarchy Processing

This section presents our computation strategy in the case where the approximation error ε is set to 0 (i.e. exact computations). Given an input edge-nested triangulation hierarchy H = {M 0 , M 1 , . . . , M h }, persistence diagrams are evaluated in three steps. First, the hierarchy is completely traversed, from the coarsest level M 0 to the finest M h , to efficiently detect topologically invariant vertices, for which the criticality can be efficiently estimated in a second step, at the last level of the hierarchy only (M h ). Third, the persistence diagram is computed from the saddle points identified at the second step. 1) Hierarchy traversal: In order to identify topologically invariant vertices, we compute the link polarity (Sec. 2.5) of the vertices for each level of the hierarchy. The link polarity of a vertex v is encoded in our setting as an array of bits, one per neighbor n of v, denoting whether n is higher or lower than v. The size of the link polarity is the same for each level of H (Sec. 2.5): at most 6 bits in 2D and 14 bits in 3D. At each level i, the polarity of new vertices is initialized, while the polarity of old vertices is updated to account for the insertion of non-monotonic vertices. The detection of non-monotonic vertices enables the fast identification of regular points: new topologically invariant vertices are known to be regular points of f i , and old topologically invariant vertices keep the same criticality at level i than at level i -1 (Sec. 2.5). We leverage these informations to avoid the explicit computation of the criticality for some of the vertices. As a new topologically invariant vertex v is identified in M i , it is flagged as a regular vertex. If its link polarity is not changed in the remaining levels (e.g. v is topologically invariant through the rest of the hierarchy), v is guaranteed to be a regular vertex of M h and no further computation will be needed for it.

2) Critical points: At M h , we compute explicitly the criticality of the vertices which are not yet guaranteed to be regular (as described above). The criticality of a vertex v is computed by enumerating the connected components of Lk + (v) and Lk -(v) (Sec. 2.2).

3) Persistence diagram: The minimum/1-saddle persistence diagram is deduced from the critical points. In particular, for each saddle point s, a backward integral line is launched from each connected component of Lk -(s), if there are more than one. These integral lines end at a local minimum m. For each integral line, we back-propagate the vertex index of m back to s. This way, we build for each saddle s the list {m 0 , m 1 , . . . , m j } of found local minima, yielding at least a minimum per connected component of sub-level sets merging at s. A merge event can thus be represented by a triplet (s, m i , m j ). All found triplets are then sorted in ascending values of s, and the merge events are processed in this order with a UnionFind (UF) data structure [START_REF] Cormen | Introduction to Algorithms[END_REF]. A persistence pair is created when two different components are merged. We refer to Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] for more details, as this last step of our method is identical to theirs.

Vertex Folding

This section describes the variations of the above strategy for the case where ε = 0, to decrease the computational cost. Specifically, we present a strategy to artificially increase the number of topologically invariant vertices during the hierarchy traversal (step 1), which will consequently result in skipping the estimation of vertex criticality (step 2) for a larger number of vertices (hence the overall speedup).

A new vertex n, appearing at a level i of the hierarchy, is inserted at the center of an old edge

(o 0 , o 1 ) of M i-1 (Sec. 2.4). Assuming that f (o 0 ) < f (o 1 ) < f (n)
, n is a non-monotonic vertex and impacts the link polarity of o 1 . The apparition of such vertices reduces the overall performance as these will trigger explicit criticality computations in step 2. To reduce the number of non-monotonic vertices inserted in M i , we can choose to reinterpolate a non-monotonic vertex n between its two new neighbors to enforce its monotony. We note the resulting monotonic vertex n and say that this vertex is folded. We define its new approximated value f ( n) as the interpolation of the approximated values of its two old neighbors:

f ( n) = f (o 0 ) + f (o 1 ) /2.
The values f (o 0 ) and f (o 1 ) are themselves either the result of a linear interpolation if o 0 or o 1 have been previously folded. Otherwise, they are equal to f (o 0 ) and f (o 1 ). The linear interpolation is preferred here to alternatives as it provides a good balance between accuracy and efficiency.

Formally, we build a sequence { f 0 , f 1 , . . . , f h } of PL scalar fields defined on each hierarchy level. The sequence is defined recursively:

1.

f 0 = f 0 2. For each old vertex o of M i , f i (o) = f i-1 (o) 3. For each new vertex n of M i that is not folded, f i (n) = f i (n) 4.
For each folded new vertex n of M i on the edge (o 0 , o 1 ),

f i ( n) = f i-1 (o 0 ) + f i-1 (o 1 ) 2
We note F i the set of folded vertices at level i (with F i ⊂ M i 0 ). By construction, folded new vertices are monotonic. Fig. 6 illustrates how a higher amount of folded vertices at level i implies a higher number of topologically invariant vertices identified on M i . In particular, if all non-monotonic vertices are folded at level i, all vertices of M i are topologically invariant.

Bottleneck Error Control

Computing persistence diagrams with vertex folding results in approximations of the exact result D( f ), given by D( f ) (diagram of the approximated field f ). Then the resulting approximation error (in terms of Bottleneck distance) is given by [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF]: The numbers denote the values of the approximated scalar field f . Red squares indicate the non-monotonic new vertices, whose folding error δ ε is labelled in red. On the left, no vertex is folded and no approximation is made on the scalar field (i.e. f = f ). It results in a high number of non-monotonic vertices and a low number of topologically invariant old and new vertices (respectively blue and green squares). In contrast, a folding threshold ε = 30% is applied on the right. Every non-monotonic vertex n with a folding error δ (n) ≤ 8 gets folded (yellow hats). This reduces the number of non-monotonic vertices and more than doubles the number of TI vertices (blue and green squares) of the full precision approach.

W ∞ D( f ), D( f ) ≤ f -f ∞ ,
which is rather easy to estimate. For each new vertex n inserted at level i of the hierarchy on the edge (o 0 , o 1 ), we define its folding error δ (n) as the difference between its original scalar value and its reinterpolation value at its level of insertion:

δ (n) = f i (o 0 ) + f i (o 1 ) /2 -f (n)
. Then, we have:

f -f ∞ = max n∈F h δ ( n) = max n∈F h | f ( n) -f ( n)|
In the light of these observations, we use the following folding strategy as we process the hierarchy. Given a target approximation error ε, the hierarchy is processed as described in Sec. 4.1, except that vertex polarity is estimated from the approximated field f . For each level i, we choose to fold non-monotonic vertices n with an error δ (n) < ε. Monotonic vertices or non-monotonic vertices with a higher folding error get added into the hierarchy without being reinterpolated. Let f ε : M → R be the final field approximation (after the hierarchy traversal is completed), given the target approximation error ε. Then it is clear that ∀v ∈ M, δ (v) ≤ ε. Thus:

W ∞ D( f ε ), D( f ) ≤ f ε -f ∞ = max n∈F h δ ( n) ≤ ε
The computation must reach the last level h, to accurately capture max δ ( n) (required for guarantees), thus it is not progressive. In the remainder, ε is noted as a percentage of the data range (ε = 0% is the exact computation, ε = 100% folds all non-monotonic vertices).

Monotony offsets

We now discuss how to handle degenerate flat plateaus, which become more frequent given our vertex folding strategy. The input scalar field f is injective on the vertices of M (Sec. 3). This is enforced in practice with an offset field O : M 0 → N, typically corresponding for each vertex to its offset in memory (Fig. 7a). O is then used to disambiguate vertices with identical scalar values.

In theory, a folded vertex n is guaranteed to be monotonic. However in practice, if δ ( n) falls below the precision of the data type used to encode the scalar field, a flat plateau emerges. This occurs frequently for instance when the input data is expressed with integers. Then, given a new folded vertex n inserted on an edge (o 0 , o 1 ), we may have: f ( n) = f (o 0 ), which means n and o 0 will be disambiguated in the algorithm by their offset O. However, this can introduce undesired monotony changes (Fig. 7b, red squares).

To guarantee the monotony of folded vertices, we introduce a monotony offset on each vertex v, noted M(v), which is modified when the vertex gets folded. The purpose of the monotony offset M is to take over the regular offset O if it contradicts the monotony of newly folded vertices. Given a new folded vertex n inserted on an edge (o 0 , o 1 ) at level i, that is non-monotonic with respect to o 0 (i.e.

f ( n) < f (o 0 ) < f (o 1 ) or f ( n) > f (o 0 ) > f (o 1 )
), we set:

M( n) =            M(o 0 ) -2 h-i if f ( n) < f (o 0 ) < f (o 1 ) and O( n) > O(o 0 ) M(o 0 ) + 2 h-i if f ( n) < f (o 0 ) < f (o 1 ) and O( n) < O(o 0 ) M(o 0 ) else
The monotony offset is initially set to zero for all new vertices and modified only in case of vertex folding. Then, the field M explicitly encodes the monotony of newly folded vertices (Fig. 7). The monotony offsets are used to disambiguate the comparison of two vertices of identical approximated scalar value in the rest of the approach (criticality estimation, integral lines, etc).

Parallelism

Our approach can be easily parallelized using shared-memory parallelism. The first step of our approach, the traversal of the hierarchy, can be trivially parallelized on the vertices, as all operations are local to a vertex (Sec. 4.2). However the hierarchy must be processed in sequential, which implies synchronization between the different levels. The computation of the criticality for non topologically invariant vertices is also completely parallel. The computation of the persistence diagram from the critical points can be parallelized on the saddle points, however locks are necessary for the parallel computation of integral lines, in order to back-propagate the indices of found extrema. The saddle points are sorted in parallel, using the GNU implementation [START_REF] Singler | The GNU libstdc++ Parallel Mode: Software Engineering Considerations[END_REF]. Finally, processing the merge events represented by the triplets is a sequential task, but represents only a fraction of the total sequential computation time (less that 1% in practice).

Uncertainty

By construction, our approach induces a Bottleneck error of ε, which corresponds to a maximum mismatch of ε between the pairs of D( f ε ) and these of D( f ). This means that some approximated persistence pairs (with a persistence below 2ε) may not be present in the exact diagram (as they may be matched to the diagonal at a cost lower than ε). We call these pairs uncertain. In contrast, the approximated pairs with a persistence beyond 2ε will certainly be present in the exact diagram, and their exact location is bounded within a square of side 2ε. We call these pairs certain. Thus, we can visually convey the level of uncertainty of our approximations directly in the persistence diagrams (Fig. 8). For this, we use a red band to indicate uncertain pairs, and we draw the bounding squares for certain pairs (Sec. 5.3).

RESULTS

This section presents experimental results obtained on a variety of datasets, available on public repositories [START_REF] Klacansky | Open Scientific Visualization Data Sets[END_REF][START_REF]TTK Data[END_REF]. We implemented our approach in C++, as modules for the Topology ToolKit (TTK) [START_REF] Bin Masood | An Overview of the Topology ToolKit[END_REF][START_REF] Tierny | The Topology ToolKit[END_REF]. The experiments were carried out on a desktop computer with two Xeon CPUs (3.0 GHz, 2×4 cores) and 64 GB of RAM. 7: Monotony offsets on a toy example of an integer field f (a: black numbers) defined on a 2D grid. The injectivity of f is guaranteed by the offset field O (a: blue numbers). As the hierarchy is processed (b), some vertices get folded (red and green squares) according to a given error threshold. Due to the precision of the field (here integer precision), their interpolated value (numbers on the right) might be identical to that of an old neighbor. In some cases (red squares), the folding actually entails a monotony change as the offset field O (blue numbers) provides the wrong order relation between neighbor vertices. If such an event occurs, the monotony offset of the folded vertex is updated (red numbers) to enforce monotony. Green squares denote folding cases where O does not contradict the folding monotony. At the finest resolution (b, rightmost), this results in plateaus (bottom row of the grid) where the injectivity of f is guaranteed and where the monotony is correctly enforced. Table 1: Increase of the number of topologically invariant (TI) vertices for different levels of approximation. For real-world datasets (MinMax and Random excluded), the average proportion of TI vertices rises from 70% for the full precision approach (i.e. ε = 0%) to 94% for a relative Bottleneck error of 5% (ε = 5%).

Dataset ∑ h i=0 |M h 0 | % TI ε = 0% ε = 1% ε = 5% ε = 10% At 931

Time Performance

The time complexity of our approach is similar to the complexity of the non-interruped algorithm by Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF]. The first two steps of our approach, that amount to the hierarchy traversal and the computation of critical points, have a linear complexity in the number of vertices: O(∑ i=h i=0 |M i 0 |). The third step, to compute the saddle-extremum persistence from the critical points, is identical to the approach by Vidal, except that we compute the persistence diagram exclusively on the last, finest level of the hierarchy. It has a practical time complexity [START_REF] Vidal | A progressive approach to scalar field topology[END_REF]

of O |M h 1 | + n s log n s + n s α(n m )
, where α stands for the inverse of the Ackermann function, and n s and n m respectively denote the number of saddle points and extrema.

Table 1 reports the number of TI vertices for all datasets, for various approximation errors. The column ε = 0 corresponds to the numbers of TI vertices reported by Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] in the exact case. For the majority of datasets, we observe a large increase in Table 2: Sequential computation times (in seconds) of our approach for the approximation of persistence diagram, for different approximation errors. The column Default reports the run time of the default approach in the Topology ToolKit [START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF]. The last column indicates the speedup of our approach with an approximation error of 5%, against the fastest of the two exact methods (left). Bold numbers indicate the smallest ε providing a speedup over reference approaches.

Dataset

Default [START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF] Progressive [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] Ours the number of TI vertices, from a proportion of 70% on average on the real-life datasets to a proportion of 90% for a small error of 1%, and 94% for a mild tolerance of 5% on the approximation. The largest increase in the number of TI vertices is reported for Backpack (being a large and noisy dataset). This table confirms that our strategy of vertex folding (Sec. 4.2) indeed implies a sensible increase in the number of TI vertices, even for mild approximation errors. Regarding the criticality estimation (step 2, Sec. 4.1), as no computation is needed for the vertices which remained topologically invariant throughout the hierarchy (Sec. 4.1), a higher proportion of TI vertices in the data is thus likely to significantly decrease the computational workload, resulting in lower computation times. Table 2 details the sequential computation times of our approach for different approximation errors. They are compared with public implementations of exact algorithms, both available in TTK [START_REF] Tierny | The Topology ToolKit[END_REF]: the progressive approach by Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] (run up to the finest resolution, hence producing an exact result), and the default algorithm used in TTK [START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF] (run at the finest hierarchy level). The last column of Tab. 2 present the speedups obtained with a Bottleneck error tolerance of 5%, compared with the fastest of either reference approaches. We observe an average reduction of the run times of 18% on real-world datasets, which confirms that our strategy of maximizing the number of TI vertices effectively reduces the computation times. The observed speedups are consistent with the increases in the proportion of TI vertices reported in Tab. 1: a large increase in the number of TI vertices implies an important reduction of the computation time. Interestingly, we find our method most beneficial on datasets that initially present a low amount of TI vertices. This usually corresponds to a high level of noise, which impedes both reference methods. In particular, the highest speedup is achieved on Backpack, our largest, noisiest real-world dataset, with a reduction of computation time of nearly 50%. In contrast, our method fails to reduce the computation times on smooth datasets such as Min-Max, Boat or Aneurism, for which the progressive approach really shines. These datasets exhibit high initial proportions of TI vertices, which limits the increase in TI vertices enabled by our approach (Tab. 1). For Random, the interpolation cost of folded vertices seems to counterbalance the speedup induced by the increase in TI vertices. Table 3 lists the computation times obtained with the parallel version of our algorithm. We find an overall average parallel efficiency of 43% with an error level of 5%, which is on par with the progressive approach [START_REF] Vidal | A progressive approach to scalar field topology[END_REF]. Although the traversal of the hierarchy can be trivially parallelize over vertices, it is subject to synchronization steps between hierarchy levels. The last step of the approach, deducing the persistence diagram from the critical points, is less balanced. Indeed, the parallel computation of integral lines between saddle points and extrema (Sec. 4.1) necessitates locks. The 3D top views show the ten most persistent maxima (spheres), corresponding to the bones of the foot (isocontour, computed on the approximated field f ε ). Our approximated diagrams correctly capture the most persistent features at a reduced computational cost.

ε = 1% ε = 5% ε = 10%

Approximation Accuracy

By design, our approach produces approximated persistence diagrams with a guaranteed bound on the Bottleneck distance to the exact result (Sec. 4.3). In the following, we additionally evaluate the accuracy of our approximation by considering the L 2 -Wasserstein distance (Sec. 2.3, computed with the efficient progressive approximation [START_REF] Vidal | Progressive wasserstein barycenters of persistence diagrams[END_REF] of TTK [START_REF] Tierny | The Topology ToolKit[END_REF]) between our approximations and the exact result. We also evaluate ff 2 to quantify the pointwise error of the approximated field f . Results are given in Table 4.

For comparison, we perform the same evaluation for a naive baseline approximation which provides identical guarantees. This baseline consists in computing, with an exact algorithm, the diagram of a staircase function f , i.e. a quantized version of the input data f , with a quantization step of 2ε. By construction, the staircase function verifies ff ∞ < ε (Fig. 9) and its persistence diagram is then guaranteed not to exceed a Bottleneck error of ε [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF]. Table 4 compares the accuracy of this baseline approximation to our algorithm. In terms of L 2 distance, our approximations of the input scalar fields are around 2 times more accurate on average (on realworld datasets, Random and MinMax excluded). This difference could be expected, as our method performs local and adaptive linear interpolations, while the staircase approach systematically flattens the data. However, our approach is also significantly more accurate with regards to the L 2 -Wasserstein distance to the exact persistence diagrams: 4 times in average on our real-world datasets.

Figure 9 further illustrates the limitations of the staircase baseline. For a given approximation error (5%), our method gives an approximated diagram that is visually more similar to the exact result, and which better depicts the number of salient features, as well as the noise in the data. In contrast, the diagram produced by the staircase approximation is more difficult to interpret, as the positions of persistence pairs are quantized on a grid in the 2D birth/death space, resulting in several co-located pairs, which cannot be distinguished visually. Our approximation of the scalar field is also closer to the exact field, both visually and in terms of the L 2 norm, enabling more accurate critical point approximations in the domain (top).

Qualitative Analysis

This section discusses the utility of our approximations from a qualitative point of view, for data analysis and visualization purposes. Fig. 1 shows the result of our approach on the Backpack dataset. In this example, our approximated diagrams correctly capture the high persistence maxima of the scalar field, which correspond to high density objects inside the bag (bottles, wires, and metallic parts of the bag). The approximate field f ε resulting from the vertex folding (Sec. 4.2) is more precise in the vicinity of features of high persistence. Indeed, the volume rendering in the top views of Fig. 1 shows clearly the objects inside the bag (high persistence maxima), while isocontours capturing the cloth of the bag (Fig. 1, bottom) illustrate Our approximated persistence diagrams correctly capture these certain pairs in the order of their persistence, with higher persistence features being detected at higher tolerance on the Bottleneck error.

the deterioration of the field in this region. The same phenomenon can be noted for the Foot dataset (Fig. 10). More vertices are folded in the vicinity of low persistence features, typically the skin of the foot, and the level of deterioration of the field increases with the approximation error. Conversely, high persistence features (bones of the toes) are well captured.

The above observation suggests that our method provides a better approximation for persistent features, and a more degraded evaluation for less persistent structures, which was an original motivation for our approach (to focus the computational efforts on relevant structures). This is confirmed in Fig. 1 by the number of noisy pairs (of low persistence, within the red band) in the approximated diagrams, which is significantly lower than the amount of noisy pairs (of identical persistence) in the exact diagram (Fig. 1, in parenthesis).

Figure 12 compares our approximations to the progressive approach by Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF]. To generate an approximated result with the progressive approach, we interrupt its computation at the penultimate hierarchy level (the computation would become exact at the final level). As documented by its authors, such intermediate results can provide useful previews, but however, with no guarantee on the approximation error. This is illustrated in Fig. 12, where the progressive approximation fails at correctly capturing the maximum of largest persistence. In contrast, our method produces persistence diagrams that correctly convey the salience of the features, and are five times more accurate, in terms of the Wasserstein distance.

Indications about the approximation uncertainty (Sec. 4.6) can be displayed in the output diagrams. Figure 11 shows our approximations for the Isabel dataset. For each approximation error, the red band indicates uncertain pairs, which may not be part of the exact result. Certain pairs are represented with a square bounding their correct location. These glyphs give a good sense of the approximation uncertainty, and are useful to assess the reliability of the diagram. For instance, a large pair in the uncertain zone may indicate the presence of a medium persistence feature in the data. This can be confirmed with the computation of a slighlty better estimation, as illustrated in Fig. 12 for the fourth most persistent feature. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF] fails at correctly capturing the global maximum (accurately detected only at the last level), resulting in a diagram that is 5 times less accurate than our 5% approximation (L 2 -Wasserstein distance to the exact result). In contrast, our approximation correctly captures the high persistence features of the data. On the far right, our 4% approximated diagram detects as certain the fourth most persistent maxima, which was marked as uncertain with ε = 5%.

LIMITATIONS AND DISCUSSION

Our approximations tend to generate much less low-persistence features than exact algorithms (Fig. 1), which can be an issue if features of interested are hidden among noisy features near the diagonal. On the upside, this characteristic of our approximations make them well suited for subsequent analysis and processing (e.g. distances and clustering), as diagrams are often thresholded in practice prior to further computations, to remove low persistence pairs anyway.

An important limitation of our approach, compared to the work of Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF], is its lack of progressivity. Indeed, to provide strong approximation guarantees, the hierarchy has to be completely traversed in our work and no intermediate result can be provided.

Another limitation is that our approach only supports saddleextremum persistence pairs at the moment. However, from our experience, these correspond in practice to the key features users tend to be interested in.

Finally, our approach provides strong guarantees on the Bottleneck distance. Future work is needed for the theoretical study of the impact of our approximations on the L 2 Wasserstein metric.

CONCLUSION

This paper introduced a method for the approximation of the persistence diagram of a scalar field. Our work revisits the progressive approach by Vidal et al. [START_REF] Vidal | A progressive approach to scalar field topology[END_REF], that generated preview diagrams upon interruption of a progressive framework. We addressed the main drawback of their approach, namely the lack of guaranteed error bounds on the diagram estimations. In contrast, we presented a novel algorithm that efficiently computes the approximation of a persistence diagram within a user controlled approximation error on the Bottleneck distance to the exact result. We showed that the approximated persistence diagrams are relevant for visualization and data analysis tasks, as they correctly describe the high persistence features in the data (i.e. the number and salience of important features), and they are more concise in practice than the exact diagrams. The uncertainty related to our approximations can be effectively depicted visually inside the diagrams.

We believe that the development of approximative approaches (with guarantees) for data analysis and visualization is an important and exciting research direction. They are especially relevant in the field of Topological Data Analysis, as most of the computational workload is typically spent in practice on the capture of small scale features (Fig. 1), whereas they usually are less relevant in the applications, and post-process simplification techniques are often applied to eliminate them anyway. In this context, a logical avenue for future work would be the development of approximative methods to revisit existing topological data representations (merge trees, Morse-Smale complexes, Reeb graphs).

Figure 2 :

 2 Figure 2: Persistence diagrams of a clean (left) and noisy (right) 2D scalar field. Critical points are represented by spheres (light brown: minima, cyan: maxima, others: saddles). The persistence diagram highlights the three main features of the data (high persistence pairs). On the right, the small pairs near the diagonal indicate the presence of noisy features with low persistence.

Figure 3 :

 3 Figure 3: Edge-nested triangulation hierarchy of a regular grid. New vertices are only inserted on existing edges, one per edge, splitting each old edge in two. Old and new edges are shown (right, M 2 ) in grey and white respectively; old and new vertices are shown in black and white respectively. New edges only connect new vertices.

Figure 4 :

 4 Figure4: Evolution of the upper link (light green) and the lower link (dark green) of a vertex v (red sphere) when progressing down the hierarchy. The topological structure of the link at level i (a) remains unchanged at level i + 1, thanks to the edge-nested property of the hierarchy. In the example (b), neighbor 5 is non-monotonic: it switches its polarity (from blue minus to red plus), which modifies the criticality of v (saddle point in M i+1 ). If the link polarity of v is not changed (c), v is called a Topologically Invariant vertex and it maintains its criticality.

Figure 5 :

 5 Figure5: Overview of our approach. First, the traversal of the hierarchy (a) enables the efficient detection of the vertices (blue spheres) that are not topologically invariant (TI), and for which the criticality must be computed. Non-monotonic vertices (red spheres) can be folded (transparent red spheres) during the traversal, i.e. reinterpolated to artificially increase the number of TI vertices. Second (b), the criticality of all non-TI vertices detected in the first step is computed, to identify the critical points of the approximated field (spheres, cyan: maxima, green: saddles, beige: minimum). Third, the saddle points are used to seed integral lines ending at extrema (c), from which the persistence diagram is deduced (d).

Figure 6 :

 6 Figure 6: Traversal of the hierarchy and identification of topological invariants vertices for two different folding threshold (left: ε = 0%, right: ε = 30%).The numbers denote the values of the approximated scalar field f . Red squares indicate the non-monotonic new vertices, whose folding error δ ε is labelled in red. On the left, no vertex is folded and no approximation is made on the scalar field (i.e. f = f ). It results in a high number of non-monotonic vertices and a low number of topologically invariant old and new vertices (respectively blue and green squares). In contrast, a folding threshold ε = 30% is applied on the right. Every non-monotonic vertex n with a folding error δ (n) ≤ 8 gets folded (yellow hats). This reduces the number of non-monotonic vertices and more than doubles the number of TI vertices (blue and green squares) of the full precision approach.

Figure

  Figure7: Monotony offsets on a toy example of an integer field f (a: black numbers) defined on a 2D grid. The injectivity of f is guaranteed by the offset field O (a: blue numbers). As the hierarchy is processed (b), some vertices get folded (red and green squares) according to a given error threshold. Due to the precision of the field (here integer precision), their interpolated value (numbers on the right) might be identical to that of an old neighbor. In some cases (red squares), the folding actually entails a monotony change as the offset field O (blue numbers) provides the wrong order relation between neighbor vertices. If such an event occurs, the monotony offset of the folded vertex is updated (red numbers) to enforce monotony. Green squares denote folding cases where O does not contradict the folding monotony. At the finest resolution (b, rightmost), this results in plateaus (bottom row of the grid) where the injectivity of f is guaranteed and where the monotony is correctly enforced.

Figure 8 :

 8 Figure 8: Approximation uncertainty visualization. The diagram at ε = 10% (a., right) exhibits one uncertain pair within the red band, which is absent from the exact result (a., left). Squares bound the correct location of certain pairs (b., transparent: exact pairs).

Figure 9 :

 9 Figure9: Comparison of our approach to a naive baseline approximation (staircase field), for the same tolerance of 5% on the Bottleneck error. Our approach (right) provides an approximated persistence diagram that resembles the exact result (center), and is around 2 times closer in terms of the Wasserstein distance (W 2 ). The high persistence critical points (spheres, 2D domain) capture more precisely the features of the data in our case. Our approximated field is also two times closer (L 2 norm) to the exact field than the naive approximation.

Figure 10 :

 10 Figure 10: Approximated persistence diagrams on the Foot dataset.The 3D top views show the ten most persistent maxima (spheres), corresponding to the bones of the foot (isocontour, computed on the approximated field f ε ). Our approximated diagrams correctly capture the most persistent features at a reduced computational cost.

Figure 11 :

 11 Figure 11: Approximated diagrams for the Isabel hurricane dataset. The most persistent maxima of the field (the magnitude of the wind velocity) are represented in the 3D view as spheres (scaled by persistence). The approximation uncertainty is visualized in the diagram: in the red band indicates uncertain persistence pairs that may not exist in the exact persistence diagram, and colored squares indicate a bound on the location of the persistence pairs that are certain to exist in the exact diagram.Our approximated persistence diagrams correctly capture these certain pairs in the order of their persistence, with higher persistence features being detected at higher tolerance on the Bottleneck error.

Figure 12 :

 12 Figure 12: Approximation of persistence diagrams on the hydrogen dataset. The best interrupted result of Vidal et al.[START_REF] Vidal | A progressive approach to scalar field topology[END_REF] fails at correctly capturing the global maximum (accurately detected only at the last level), resulting in a diagram that is 5 times less accurate than our 5% approximation (L 2 -Wasserstein distance to the exact result). In contrast, our approximation correctly captures the high persistence features of the data. On the far right, our 4% approximated diagram detects as certain the fourth most persistent maxima, which was marked as uncertain with ε = 5%.

Table 3 :

 3 Parallel computation times (in seconds, on 8 physical cores) of our algorithm with different approximation errors. The presented speedups relate to the sequential run times.

	Dataset	Default [34] Progressive [78] ε = 1% speedup ε = 5% speedup ε = 10% speedup
	At	0.20	0.07	0.05	2.61	0.07	2.24	0.07	2.25
	SeaSurfaceHeight	0.21	0.13	0.08	3.65	0.07	3.91	0.07	3.77
	EthaneDiol	0.21	0.11	0.08	3.69	0.09	3.60	0.10	3.35
	Hydrogen	0.72	0.20	0.13	3.18	0.16	2.95	0.17	2.80
	Isabel	0.47	0.46	0.54	1.77	0.28	3.15	0.30	3.05
	Combustion	0.34	0.51	0.54	1.83	0.31	2.92	0.25	3.41
	Boat	0.29	0.29	0.33	2.94	0.37	2.85	0.34	2.98
	MinMax	0.80	0.61	0.69	3.11	0.54	4.23	0.54	3.97
	Aneurism	2.03	1.51	1.62	2.23	1.36	2.59	1.33	2.25
	Foot	3.12	2.34	2.65	3.91	1.91	4.25	2.04	3.01
	Heptane	2.44	2.35	2.08	2.60	1.33	3.88	1.51	3.41
	Random	27.04	8.47	7.15	4.05	7.77	3.74	8.89	3.43
	Backpack	30.36	24.88	12.63	4.92	8.50	4.72	7.02	4.53

Table 4 :

 4 Accuracy comparison between our approach and a naive baseline approximation based on the computation of a staircase function, for the same relative Bottleneck error of 5%. Our approximations are more accurate in average, both in terms of the L 2 distance of the approximated field (2 times more accurate) and in terms of the L 2 -Wasserstein distance to the exact persistence diagram (5 times).

	Dataset	Staircase L 2	Ours L 2	Ratio Staircase W 2	Ours W 2 Ratio
	At	276.66	75.31	3.67	3.31	1.01	3.29
	SeaSurfaceHeight	92.3	58.05	1.59	8.75	1.69	5.17
	EthaneDiol	337.55	73.2	4.61	1.78	0.61	2.9
	Hydrogen	11.0	13.0	0.85	25.69	12.73	2.02
	Isabel	3,591.99	1,569.98	2.29	42.26	8.08	5.23
	Combustion	38.04	17.59	2.16	0.87	0.21	4.1
	Boat	24.18	9.11	2.66	2.14	0.49	4.37
	MinMax	117.55	0.37 314.77	0.0	0.0	-
	Aneurism	6.0	9.0	0.67	1,198.14	128.61	9.32
	Foot	16,006.11 10,784.33	1.48	5,859.67 1,419.51	4.13
	Heptane	5.0	12.0	0.42	716.1	61.12 11.72
	Random	29,270.73	3,468.14	8.44	23,863.56	1,038.1 22.99
	Backpack	142.0	142.0	1.0	17,941.13 1,933.83	9.28

ACKNOWLEDGMENTS

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