Persistence Images: A Stable Vector Representation of Persistent Homology
Abstract
Many datasets can be viewed as a noisy sampling of an underlying space, and tools from topological data analysis can characterize this structure for the purpose of knowledge discovery. One such tool is persistent homology, which provides a multiscale description of the homological features within a dataset. A useful representation of this homological information is a persistence diagram (PD). Efforts have been made to map PDs into spaces with additional structure valuable to machine learning tasks. We convert a PD to a finitedimensional vector representation which we call a persistence image (PI), and prove the stability of this transformation with respect to small perturbations in the inputs. The discriminatory power of PIs is compared against existing methods, showing significant performance gains. We explore the use of PIs with vectorbased machine learning tools, such as linear sparse support vector machines, which identify features containing discriminating topological information. Finally, high accuracy inference of parameter values from the dynamic output of a discrete dynamical system (the linked twist map) and a partial differential equation (the anisotropic KuramotoSivashinsky equation) provide a novel application of the discriminatory power of PIs. Keywords: topological data analysis, persistent homology, persistence images, machine learning, dynamical systems
1 Introduction
In recent years, the field of topology has grown to include a large set of computational tools (Edelsbrunner and Harer, 2010). One of the fundamental tools is persistent homology, which tracks how topological features appear and disappear in a nested sequence of topological spaces (Edelsbrunner and Harer, 2008; Zomorodian and Carlsson, 2005). This multiscale information can be represented as a persistence diagram (PD), a collection of points in the plane where each point corresponds to a topological feature that appears at scale and disappears at scale . We say the feature has a persistence value of . This compact summary of topological characteristics by finite multisets of points in the plane is responsible, in part, for the surge of interest in applying persistent homology to the analysis of complex, often highdimensional data. Computational topology has been successfully applied to a broad range of datadriven disciplines (Perea and Harer, 2013; Dabaghian et al., 2012; Chung et al., 2009; Heath et al., 2010; Singh et al., 2008; Topaz et al., 2015; Pearson et al., 2015).
Concurrent with this revolution in computational topology, a growing general interest in data analysis has driven advances in data mining, pattern recognition, and machine learning (ML). Since the space of PDs can be equipped with a metric structure (bottleneck or Wasserstein (Mileyko et al., 2011; Turner et al., 2014)), and since these metrics reveal the stability of PDs under small perturbations of the data they summarize (CohenSteiner et al., 2007, 2010; Chazal et al., 2014), it is possible to perform a variety of ML techniques using PDs as a statistic for clustering data sets. However, many other useful ML tools and techniques (e.g. support vector machines (SVM), decision tree classification, neural networks, feature selection, and dimension reduction methods) require more than a metric structure. In addition, the cost of computing the bottleneck or Wasserstein distance grows quickly as the number of offdiagonal points in the diagrams increases (Di Fabio and Ferri, 2015). To resolve these issues, considerable effort has been made to map PDs into spaces which are suitable for other ML tools (Bubenik, 2015; Reininghaus et al., 2015; Bendich et al., 2014; Adcock et al., 2013; Donatini et al., 1998; Ferri et al., 1997; Chung et al., 2009; Pachauri et al., 2011; Bendich et al., 2015; Chen et al., 2015; Carrière et al., 2015; Di Fabio and Ferri, 2015). Each approach has benefits and drawbacks, which we review in §2. With these in mind, we are led to pose the following question:
Problem Statement: How can we represent a persistence diagram so that:

the output of the representation is a vector in ,

the representation is stable with respect to input noise,

the representation is efficient to compute,

the representation maintains an interpretable connection to the original PD, and

the representation allows one to adjust the relative importance of points in different regions of the PD?
The main contribution of this paper is to provide a finitedimensionalvector representation of a PD called a persistence image (PI). We first map a persistence diagram to an integrable function called a persistence surface. The surface is defined as a weighted sum of Gaussian functions^{1}^{1}1Or more generally, a weighted sum of probability density functions, one centered at each point in the PD. The idea of persistence surfaces has appeared even prior to the development of persistent homology, in Donatini et al. (1998) and Ferri et al. (1997). Taking a discretization of a subdomain of defines a grid. A persistence image, i.e. a matrix of pixel values, can be created by computing the integral of on each grid box. This PI is a “vectorization” of the PD, and provides a solution to the problem statement above.
Criteria (i) is our primary motivation for developing PIs. A large suite of ML techniques and statistical tools (means and variances) already exist to work with data in . Additionally, such a representation allows for the use of various distance metrics (norms and angle based metrics) and other measures of (dis)similarity. The remaining criteria of the problem statement (iiv) further ensure the usefulness of this representation.
The desired flexibility of 5 is accomplished by allowing one to build a PI as a weighted sum of Gaussians, where the weightings may be chosen from a broad class of weighting functions.^{2}^{2}2Weighting functions are restricted only to the extent necessary for our stability results in §5. For example, a typical interpretation is that points in a PD of high persistence are more important than points of low persistence (which may correspond to noise). One may therefore build a PI as a weighted sum of Gaussians where the weighting function is nondecreasing with respect to the persistence value of each PD point. However, there are situations in which one may prefer different measures of importance. Indeed, Bendich et al. (2015) find that, in their regression task of identifying a human brain’s age from its arterial geometry, it is the points of medium persistence (not high persistence) that best distinguish the data. In such a setting, one may choose a weighting function with largest values for the points of medium persistence. In addition, the Homology Inference Theorem (CohenSteiner et al., 2007) states that when given a sufficiently dense finite sample from a space , it is the points in the PD with sufficiently small birth times (and sufficiently high persistence) which recover the homology groups of the space; hence one may choose a weighting function that emphasizes points near the deathaxis and away from the diagonal, as indicated in the leftmost yellow rectangle of Bendich (2009, Figure 2.4). A potential disadvantage of the flexibility in 5 is that it requires a choice; however, prior knowledge of one’s particular problem may inform that choice. Moreoever, our examples illustrate the effectiveness of a standard choice of weighting function that is nondecreasing with the persistence value.
The remainder of this article is organized as follows. Related work connecting topological data analysis and ML is reviewed in §2, and §3 gives a brief introduction to persistent homology, PDs from point cloud data, PDs from functions, and the bottleneck and Wasserstein metrics. PIs are defined in §4 and their stability with respect to the 1Wasserstein distance between PDs is proved in §5. Lastly, §6 contains examples of ML techniques applied to PIs generated from samples of common topological spaces, an applied dynamical system modeling turbulent mixing, and a partial differential equation describing pattern formation in extended systems driven far from equilibrium. Our code for producing PIs is publicly available at https://github.com/CSUTDA/PersistenceImages.
2 Related Work
The space of PDs can be equipped with the bottleneck or Wasserstein metric (defined in §3), and one reason for the popularity of PDs is that these metrics are stable with respect to small deviations in the inputs (CohenSteiner et al., 2007, 2010; Chazal et al., 2014). Furthermore, the bottleneck metric allows one to define Fréchet means and variances for a collection of PDs (Mileyko et al., 2011; Turner et al., 2014). However, the structure of a metric space alone is insufficient for many ML techniques, and a recent area of interest in the topological data analysis community has been encoding PDs in ways that broaden the applicability of persistence. For example, Adcock et al. (2013) study a ring of algebraic functions on the space of persistence diagrams, and Verovšek (2016) identifies tropical coordinates on the space of diagrams. Ferri and Landi (1999) and Di Fabio and Ferri (2015) encode a PD using the coefficients of a complex polynomial that has the points of the PD as its roots.
Bubenik (2015) develops the notion of a persistence landscape, a stable functional representation of a PD that lies in a Banach space. A persistence landscape (PL) is a function , which can equivalently be thought of as a sequence of functions . For the landscape distance between two landscapes and is defined as ; the landscape distance is stable with respect to the bottleneck distance on PDs, and the landscape distance is continuous with respect to the Wasserstein distance on PDs. One of the motivations for defining persistence landscapes is that even though Fréchet means of PDs are not necessarily unique (Mileyko et al., 2011), a set of persistence landscapes does have a unique mean. Unique means are also a feature of PIs as they are vector representations. An advantage of PLs over PIs is that the map from a PD to a PL is easily invertible; an advantage of PIs over PLs is that PIs live in Euclidean space and hence are amenable to a broader range of ML techniques. In §6, we compare PDs, PLs, and PIs in a classification task on synthetic data sampled from common topological spaces. We find that PIs behave comparably or better than PDs when using ML techniques available to both representations, but PIs are significantly more efficient to compute. Also, PIs outperform PLs in the majority of the classification tasks and are of comparable computational efficiency.
A vector representation of a PD, due to Carrière et al. (2015), can be obtained by rearranging the entries of the distance matrix between points in a PD. In their Theorem 3.2, they prove that both the and norms between their resulting vectors are stable with respect to the bottleneck distance on PDs. They remark that while the norm is useful for nearestneighbor classifiers, the norm allows for more elaborate algorithms such as SVM. However, though their stability result for the norm is wellbehaved, their constant for the norm scales undesirably with the number of points in the PD. We provide this as motivation for our Theorem 4, in which we prove the , , and norms for PI vectors are stable with respect to the 1Wasserstein distance between PDs, and in which none of the constants depend on the number of points in the PD.
By superimposing a grid over a PD and counting the number of topological features in each bin, Bendich et al. (2014) create a feature vector representation. An advantage of this approach is that the output is easier to interpret than other more complicated representations, but a disadvantage is that the vectors are not stable for two reasons:

an arbitrarily small movement of a point in a PD may move it to another bin, and

a PD point emerging from the diagonal creates a discontinuous change.
Source (i) of instability can be improved by first smoothing a PD into a surface. This idea has appeared multiple times in various forms — even prior to the development of persistent homology, Donatini et al. (1998) and Ferri et al. (1997) convert size functions (closely related to 0dimensional PDs) into surfaces by taking a sum of Gaussians centered on each point in the diagram. This conversion is not stable due to (ii), and we view our work as a continued study of these surfaces, now also in higher homological dimensions, in which we introduce a weighting function^{3}^{3}3Our weighting function is continuous and zero for points of zero persistence, i.e. points along the diagonal. to address (ii) and obtain stability. Chung et al. (2009) produce a surface by convolving a PD with the characteristic function of a disk, and Pachauri et al. (2011) produce a surface by centering a Gaussian on each point, but both of these methods lack stability again due to (ii). Surfaces produced from random PDs are related to the empirical intensity plots of Edelsbrunner et al. (2012).
Reininghaus et al. (2015) produce a stable surface from a PD by taking the sum of a positive Gaussian centered on each PD point together with a negative Gaussian centered on its reflection below the diagonal; the resulting surface is zero along the diagonal. This approach is similar to ours, and indeed we use (Reininghaus et al., 2015, Theorem 3) to show that our persistence surfaces are stable only with respect to the 1Wasserstein distance (Remark 1). Nevertheless, we propose our independentlydeveloped surfaces as an alternative stable representation of PDs with the following potential advantages. First, our sum of nonnegatively weighted Gaussians may be easier to interpret than a sum including negative Gaussians. Second, we produce vectors from our surfaces with wellbehaved stability bounds, allowing one to use vectorbased learning methods such as linear SVM. Indeed, Zeppelzauer et al. (2016) report that while the kernel of Reininghaus et al. (2015) can be used with nonlinear SVMs, in practice this becomes inefficient for a large number of training vectors because the entire kernel matrix must be computed. Third, while the surface of Reininghaus et al. (2015) weights persistence points further from the diagonal more heavily, there are situations in which one may prefer different weightings, as discussed in §1 and item 5 of our Problem Statement. Hence, one may want weightings on PD points that are nonincreasing or even decreasing when moving away from the diagonal, an option available in our approach.
We produce a persistence surface from a PD by taking a weighted sum of Gaussians centered at each point. We create vectors, or PIs, by integrating our surfaces over a grid, allowing ML techniques for finitedimensional vector spaces to be applied to PDs. Our PIs are stable, and distinct homology dimensions may be concatenated together into a single vector to be analyzed simultaneously. Our surfaces are studied from the statistical point of view by Chen et al. (2015), who cite a preprint version of our work; their applications in Section 4 use the norm between these surfaces, which can be justified as a reasonable notion of distance due to our Theorem 3 that proves the distance between such surfaces is stable.
Zeppelzauer et al. (2016) apply persistent images to 3D surface analysis for archeological data, in which the machine learning task is to distinguish scans of natural rock surfaces from those containing ancient humanmade engravings. The authors state they select PIs over other topological methods because PIs are computationally efficient and can be used with a broader set of ML techniques. PIs are compared to an aggregate topological descriptor for a persistence diagram: the first entry of this vector is the number of points in the diagram, and the remaining entries are the minimum, maximum, mean, standard deviation, variance, 1stquartile, median, 3rdquartile, sum of square roots, sum, and sum of squares of all the persistence values. In their three experiments, the authors find the following.

When classifying natural rock surfaces from engravings using persistent diagrams produced from the sublevel set filtration, PIs outperform the aggregate descriptor.

When the natural rock and engraved surfaces are first preprocessed using the completed local binary pattern (CLBP) operator for texture classificiation (Guo et al., 2010), PIs outperform the aggregate descriptor.

The authors added PIs and the aggregate descriptor to eleven different nontopological baseline descriptors, and found that the classification accuracy of the baseline descriptor was improved more by the addition of PIs than by the addition of the aggregate descriptor.
Furthermore, Table 1 of Zeppelzauer et al. (2016) demonstrates that for their machine learning task, PIs have low sensitivity to the parameter choices of resolution and variance (§4).
3 Background on Persistent Homology
Homology is an algebraic topological invariant that, roughly speaking, describes the holes in a space. The dimensional holes (connected components, loops, trapped volumes, etc.) of a topological space are encoded in an algebraic structure called the th homology group of , denoted . The rank of this group is referred to as the th Betti number, , and counts the number of independent dimensional holes. For a comprehensive study of homology, see Hatcher (2002).
Given a nested sequence of topological spaces , the inclusion for induces a linear map on the corresponding th homology for all . The idea of persistent homology is to track elements of as the scale (or “time”) parameter increases (Edelsbrunner and Harer, 2008; Zomorodian and Carlsson, 2005; Edelsbrunner and Harer, 2010). A standard way to represent persistent homology information is a persistence diagram (PD)^{4}^{4}4Another standard representation is the barcode (Ghrist, 2008)., which is a multiset of points in the Cartesian plane . For a fixed choice of homological dimension , each homological feature is represented by a point , whose birth and death indices and are the scale parameters at which that feature first appears and disappears, respectively. Since all topological features die after they are born, necessarily each point appears on or above the diagonal line . A PD is a multiset of such points, as distinct topological features may have the same birth and death coordinates.^{5}^{5}5By convention, all points on the diagonal are taken with infinite multiplicity. This facilitates the definitions of the Wasserstein and bottleneck distances below. Points near the diagonal are often considered to be noise while those further from the diagonal represent more robust topological features.
In this paper, we produce PDs from two different types of input data:

When our data is a a point cloud, i.e. a finite set of points in some space, then we produce PDs using Vietoris–Rips filtration.

When our data is a realvalued function, then we produce PDs using the sublevel set filtration.^{6}^{6}6As explained in §A.3, (i) can be viewed as a special case of (ii).
For setting (i), point cloud data often comes equipped with a metric or a measure of internal similarity and is rich with latent geometric content. One approach to identifying geometric shapes in data is to consider the dataset as the vertices of a simplicial complex and to add edges, triangles, tetrahedra, and higherdimensional simplices whenever their diameter is less than a fixed choice of scale. This topological space is called the Vietoris–Rips simplicial complex, which we introduce in more detail in §A.2. The homology of the Vietoris–Rips complex depends crucially on the choice of scale, but persistent homology eliminates the need for this choice by computing homology over a range of scales (Carlsson, 2009; Ghrist, 2008). In §6.1–6.4.1, we obtain PDs from point cloud data using the Vietoris–Rips filtered simplicial complex, and we use ML techniques to classify the point clouds by their topological features.
In setting (ii), our input is a real valued function defined on some domain . One way to understand the behavior of map is to understand the topology of its sublevel sets . By letting increase, we obtain an increasing sequence of topological spaces, called the sublevel set filtration, which we introduce in more detail in §A.3. In §6.4.2, we obtain PDs from surfaces produced from the KuramotoSivashinsky equation, and we use ML techniques to perform parameter classification.
In both settings the output of the persistent homology computation is a collection of PDs encoding homological features of the data across a range of scales. Let denote the set of all PDs. The space can be endowed with metrics as studied by CohenSteiner et al. (2007) and Mileyko et al. (2011). The Wasserstein distance defined between two PDs and is given by
where and ranges over bijections between and . Another standard choice of distance between diagrams is referred to as the bottleneck distance. These metrics allow us to measure the (dis)similarity between the homological characteristics of two datasets.
4 Persistence Images
We propose a method for converting a PD into a vector while maintaining an interpretable connection to the original PD. Figure 1 illustrates the pipeline from data to PI starting with spectral and spatial information in from an immunofluorescent image of a circulating tumor cell (Emerson et al., 2015).
Precisely, let be a PD in birthdeath coordinates^{7}^{7}7We omit points that correspond to features with infinite persistence, e.g. the feature corresponding to the connectedness of the complete simplicial complex.. Let be the linear transformation , and let be the transformed multiset in birthpersistence coordinates^{8}^{8}8Instead of birthpersistence coordinates, one could also use other choices such as birthdeath or (average size)persistence coordinates. Our stability results (§5) still hold with only a slight modification to the constants., where each point corresponds to a point . Let be a differentiable probability distribution with mean . In all of our applications, we choose this distribution to be the normalized symmetric Gaussian with mean and variance defined as
Fix a nonnegative weighting function that is zero along the horizontal axis, continuous, and piecewise differentiable. With these ingredients we transform the PD into a scalar function over the plane.
Definition 1.
For a PD, the corresponding persistence surface is the function
The weighting function is critical to ensure the transformation from a PD to a persistence surface is stable, which we prove in §5.
Finally, the surface is reduced to a finitedimensional vector by discretizing a relevant subdomain and integrating over each region in the discretization. In particular, we fix a grid in the plane with boxes (pixels) and assign to each the integral of over that region.
Definition 2.
For a PD, its persistence image is the collection of pixels .
PIs provide a convenient way to combine PDs of different homological dimensions into a single object. Indeed, suppose in an experiment the PDs for , , …, are computed. One can concatenate the PI vectors for , , …, into a single vector representing all homological dimensions simultaneously, and then use this concatenated vector as input into ML algorithms.
There are three choices the user makes when generating a PI: the resolution, the distribution (and its associated parameters), and the weighting function.
Resolution of the image: The resolution of the PI corresponds to the grid being overlaid on the PD. The classification accuracy in the PI framework appears to be fairly robust to choice of resolution, as discussed in §6.2 and (Zeppelzauer et al., 2016).
The Distribution: Our method requires the choice of a probability distribution which is associated to each of the points in the PD. The examples in this paper use a Gaussian centered at each point, but other distributions may be used. The Gaussian distribution depends on a choice of variance: we leave this choice as an open problem, though the experiments in §6.2 and (Zeppelzauer et al., 2016) show a low sensitivity to the choice of variance.
The Weighting Function: In order for our stability results in §5 to hold, our weighting function must be zero along the horizontal axis (the analogue of the diagonal in birthpersistence coordinates), continuous, and piecewise differentiable. A simple choice is a weighting function that depends only on the vertical persistence coordinate . In order to weight points of higher persistence more heavily, functions which are nondecreasing in , such as sigmoidal functions, are a natural choice. However, in certain applications such as Bendich et al. (2015) it may be points of small or medium persistence that perform best for ML tasks, and hence, one may choose to use more general weighting functions. In our experiments in §6, we use a piecewise linear weighting function which only depends on the persistence coordinate . Given , define via
We use , where is the persistence value of the most persistent feature in all trials of the experiment.
In the event that the birth coordinate is zero for all points in the PD, as is often the case for , it is possible to generate a 1dimensional (instead of 2dimensional) PI using 1dimensional distributions. This is the approach we adopt. Appendix B displays examples of PIs for the common topological spaces of a circle and a torus with various parameter choices.
5 Stability of Persistence Surfaces and Images
Due to the unavoidable presence of noise or measurement error, tools for data analysis ought to be stable with respect to small perturbations of the inputs. Indeed, one reason for the popularity of PDs in topological data analysis is that the transformation of a data set to a PD is stable (Lipschitz) with respect to the bottleneck metric and – given some mild assumptions about the underlying data – is also stable with respect to the Wasserstein metrics (Edelsbrunner and Harer (2010)). In §5.1, we show that persistence surfaces and images are stable with respect to the 1Wasserstein distance between PDs. In §5.2, we prove stability with improved constants when the PI is constructed using the Gaussian distribution.
5.1 Stability for general distributions
For differentiable, define to be the maximal norm of the gradient vector of , i.e. the largest directional derivative of . It follows by the fundamental theorem of calculus for line integrals that for all , we have
(1) 
We may safely denote by and by since the maximal directional derivative and supremum of a fixed differentiable probability distribution are invariant under translation. Note that
(2) 
since for any we have .
Recall that our nonnegative weighting function is defined to be zero along the horizontal axis, continuous, and piecewise differentiable.
Lemma 1.
For , we have .
Theorem 1.
The persistence surface is stable with respect to the 1Wasserstein distance between diagrams: for we have
Proof.
Since we assume and consist of finitely many points, there exists a matching that achieves the infimum in the Wasserstein distance. Then
by Lemma 1.  
∎
It follows that persistence images are also stable.
Theorem 2.
The persistence image is stable with respect to the 1Wasserstein distance between diagrams. More precisely, if is the maximum area of any pixel in the image, is the total area of the image, and is the number of pixels in the image, then
The constant for the norm bound containing goes to infinity as the resolution of the image increases. For this reason, in Theorem 4 we provide bounds with better constants in the specific case of Gaussian distributions.
Proof.
Remark 1.
Recall is the set of all PDs. The kernel defined by is nontrivial and additive, and hence Theorem 3 of Reininghaus et al. (2015) implies that is not stable with respect to for any . That is, when there is no constant such that for all we have .
5.2 Stability for Gaussian distributions
In this section, we provide stability results with better constants in the case of Gaussian distributions. With Gaussian distributions we can control not only the distance but also the distance between two persistence surfaces.
Our results for 2dimensional Gaussians will rely on the following lemma for 1dimensional Gaussians.
Lemma 2.
For , let be the normalized 1dimensional Gaussians, defined via . If , then
Proof.
Let . We show in Appendix C that
(3) 
where is defined by
Furthermore, function is differentiable for . By searching for real roots of , one can show
The result follows by letting . ∎
Lemma 3.
For , let be normalized 2dimensional Gaussians. Then
The proof of Lemma 3 is shown in Appendix C and uses a similar construction to that of Lemma 2. We are prepared to prove the stability of persistence surfaces with Gaussian distributions.
Theorem 3.
The persistence surface with Gaussian distributions is stable with respect to the 1Wasserstein distance between diagrams: for we have
Proof.
Since we assume and consist of finitely many offdiagonal points, there exists a matching that achieves the infimum in the Wasserstein distance. Then
∎
It follows that persistence images are also stable.
Theorem 4.
The persistence image with Gaussian distributions is stable with respect to the 1Wasserstein distance between diagrams. More precisely,
Proof.
6 Experiments
In order to assess the addded value of our vector representation of PDs, we compare the performance of PDs, PLs, and PIs in §6.1 in a classification task for a synthetic data set consisting of point clouds sampled from six different topological spaces using medoids, which utilizes only the internal dissimilarities of each of the topological data descriptions to classify. We find that PIs produce consistently high classification accuracy and, furthermore, the computation time for PIs is significantly faster than computing bottleneck or Wasserstein distances between PDs. In §6.2, we explore the impact that the choices of parameters determining our PIs have on classification accuracy. We find that the accuracy is insensitive to the particular choices of PI resolution and distribution variance. In §6.3, we combine PIs with a sparse support vector machine classifier to identify the most strongly differentiating pixels for classification; this is an example of a ML task which is facilitated by the fact that PIs are finite vectors. Finally, as a novel machine learning application, we illustrate the utility of PIs to infer dynamical parameter values in both continuous and discrete dynamical systems: a discrete time system called the linked twist map in §6.4.1, and a partial differential equation called the anisotropic KuramotoSivashinsky equation in §6.4.2.
6.1 Comparison of PDs, PLs, and PIs using medoids Classification
Our synthetic dataset consists of six shape classes: a solid cube, a circle, a sphere, three clusters, three clusters within three clusters, and a torus. Given a level of Gaussian random noise, we produce 25 point clouds of 500 randomly sampled noisy points from each of the six shapes; giving 150 point clouds in total. We then compute the and PDs for the Vietoris–Rips filtration (§A.2) built from each point cloud which have been endowed with the ambient Euclidean metric on
Our goal is to compare various methods for transforming PDs into distance matrices to be used to establish proximity of topological features extracted from data. We create distance matrices of size , using three choices of representation (PDs, PLs, PIs), three choices of metric (, , )^{9}^{9}9The , , distances on PDs are more commonly known as the 1Wasserstein, 2Wasserstein, and bottleneck distances., two choices of Gaussian noise (, ), and two homological dimensions (, ). For example, the PD, , , , distance matrix contains the 2Wasserstein distances between the PDs for the random point clouds with noise level . By contrast, the PI, , , distance matrix contains all pairwise distances between the PIs^{10}^{10}10For PIs in this experiment, we use variance , resolution , and the weighting function defined in §4. produced from the PDs with noise level .
We first compare these distance matrices based on how well they classify the random point clouds into shape classes via medoids clustering (Kaufman and Rousseeuw, 1987; Park and Jun, 2009). medoids produces a partition of a metric space into clusters by choosing points from the dataset called medoids and assigning each metric space point to its closest medoid. The score of such a clustering is the sum of the distances from each point to its closest medoid. The desired output of medoids is the clustering with the minimal clustering score. Unfortunately, an exhaustive search for the global minimum is often prohibitively expensive. A typical approach to search for this global minimum is to choose a large selection of random initial medoids, improve each selection of medoids iteratively in rounds until the clustering score stabilizes and then return the identified final clustering with the lowest score for each initialization. In our experiments, we choose 1,000 random initial selections of medoids (as there are six shape classes) for each distance matrix, improve each selection of medoids using the Voronoi iteration method (Park and Jun, 2009), and return the clustering with the lowest classification score. To each medoids clustering we assign an accuracy which is equal to the percentage of random point clouds identifed with a medoid of the same shape class. In Table 1, we report the classification accuracy of the medoids clustering with the lowest clustering score, for each distance matrix.
Our second criterion for comparing methods to produce distance matrices is computational efficiency. In Table 1, we report the time required to produce each distance matrix, starting with 150 precomputed PDs as input. In the case of PLs and PIs, this time includes the intermediate step of transforming each PD into the alternate representation, as well as computing the pairwise distance matrix. All timings are computed on a laptop with a 1.3 GHz Intel Core i5 processor and 4 GB of memory. We compute bottleneck, 1Wasserstein, and 2Wasserstein distance matrices using the software of Kerber et al. (2016). For PL computations, we use the Persistence Landscapes Toolbox by Bubenik and Dlotko (2016). Our MATLAB code for producing PIs is publically available at https://github.com/CSUTDA/PersistenceImages.
Distance Matrix  Accuracy (Noise 0.05)  Time (Noise 0.05)  Accuracy (Noise 0.1)  Time (Noise 0.1) 

PD, ,  96.0%  37346s  96.0%  42613s 
PD, ,  91.3%  24656s  91.3%  25138s 
PD, ,  60.7%  1133s  63.3%  1149s 
PD, ,  100%  657s  96.0%  703s 
PD, ,  100%  984s  97.3%  1042s 
PD, ,  81.3%  527s  66.7%  564s 
PL, ,  92.7%  29s  96.7%  33s 
PL, ,  77.3%  29s  82.0%  34s 
PL, ,  60.7%  2s  63.3%  2s 
PL, ,  83.3%  36s  80.7%  48s 
PL, ,  83.3%  50s  66.7%  69s 
PL, ,  74.7%  8s  66.7%  9s 
PI, ,  93.3%  9s  95.3%  9s 
PI, ,  92.7%  9s  95.3%  9s 
PI, ,  94.0%  9s  96.0%  s 
PI, ,  100%  17s  95.3%  18s 
PI, ,  100%  17s  96.0%  18s 
PI, ,  100%  17s  96.0%  18s 
We see in Table 1 that PI distance matrices have higher classification accuracy than nearly every PL distance matrix, and higher classification accuracy than PDs in half of the trials. Furthermore, the computation times for PI distance matrices are significantly lower than the time required to produce distance matrices from PDs using the bottleneck or Wasserstein metrics. In this experiment, persistent images provide a representation of persistent diagrams which is both useful for the classification task and also computationally efficient.
6.2 Effect of PI Parameter Choice
In any system that relies on multiple parameters, it is important to understand the effect of parameter values on the system. As such, we complete a search of the parameter space used to generate PIs on the shape dataset described in §6.1 and measure medoids classification accuracy as a function of the parameters. We explore 20 different resolutions (ranging from to in increments of 5), use a Gaussian function with 20 different choices of variance (ranging from 0.01 to 0.2 in increments of 0.01), and the weighting function described in §4.
For each set of parameters, we compute the classification accuracy of the medoids clustering with the minimum clustering score on the two sets of noise levels for the homology dimensions and . We observe that the classification accuracy is insensitive to the choice of resolution and variance. The plots in Figure 2 are characteristic of the 2dimensional accuracy surface over all combinations of parameters in the ranges of variances and resolutions we tested. In an application to archeology, Zeppelzauer et al. (2016) find a similar robustness of PIs to the choices of resolution and variance.
6.3 Differentiating Homological Features by Sparse Support Vector Machine
The 1norm regularized linear support vector machine (SVM), a.k.a. sparse SVM (SSVM) classifies data by generating a separating hyperplane that depends on very few input space features (Bradley and Mangasarian, 1998; Zhu et al., 2004; Zhang and Zhou, 2010). Such a model can be used for reducing data dimension or selecting discriminatory features. Note that linear SSVM feature selection is implemented on vectors and, therefore, can be used on our PIs to select discriminatory pixels during classification. Other PD representations in the literature (Reininghaus et al., 2015; Pachauri et al., 2011) are designed to use kernel ML methods, such as kernel (nonlinear) SVMs. However, constructing kernel SVM classifiers using the 1norm results in minimizing the number of kernel functions, not the number of features in the input space (i.e. pixels in our application) (Fung and Mangasarian, 2004). Hence, for the purpose of feature selection or, more precisely, PI pixel selection, we employ the linear SSVM.
We adopt the oneagainstall (OAA) SSVM on the sets of and PIs from the six class shape data. In a oneagainstall SSVM, there is one binary SSVM for each class to separate members of that class from members of all other classes. The PIs were generated using resolution , variance , and noise level . Note that because of the resolution parameter choice of , each PI is a dimensional vector, and the selected features will be a subset of indices corresponding to pixels within the PI. Using 5fold crossvalidated SSVM resulted in accuracy comparing six sparse models with indications of the discriminatory features. Feature selection is performed by retaining the features (again, in this application, pixels) with nonzero SSVM weights, determined by magnitude comparison using weight ratios; for details see Chepushtanova et al. (2014). Figure 3 provides two examples, indicating the pixels of PIs that discriminate circles and tori from the other classes in the synthetic data set.
Feature selection produces highly interpretable results. The discriminatory pixels in the PIs that separate circles from the other classes correspond to the region where highly persistent topological features exist across all samples of a noisy circle (highlighted in Figure 3a). Alternatively, the discriminatory pixels in PIs that separate tori from the other classes correspond to points of short to moderate persistence (see Figure 3b). In this way, Figure 3b reiterates an observation of Bendich et al. (2015) that points of short to moderate persistence can contain important discriminatory information. Similar conclusions can be drawn from the discriminatory pixels of others classes (Appendix D). Our classification accuracy of 100% is obtained using only those pixels selected by SSVM (a cumulative set of only 10 distinct pixels).
6.4 Application: Determination of Dynamical System Parameters
Models of dynamic physical phenomenon rarely agree perfectly with the reality they represent. This is often due to the presence of poorlyresolved (or poorlyunderstood) processes which are parameterized rather than treated explicitly. As such, determination of the influence of a model parameter — which may itself be an incompletelydescribed conglomeration of several physical parameters — on model dynamics is a mainstay of dynamical system analysis. In the case of fitting a dynamic model to data, i.e. explicit determination of optimal model parameters, a variety of techniques exist for searching through parameter space, which often necessitate costly simulations. Furthermore, such approaches struggle when applied to models exhibiting sensitivity to initial conditions. We recast this problem as a machinelearning exercise based on the hypotheses that model parameters will be reflected directly in dynamic data in a way made accessible by persistent homology.
6.4.1 A discrete dynamical model
We approach a classification problem with data arising from the linked twist map, a discrete dynamical system modeling fluid flow. Hertzsch et al. (2007) use the linked twist map to model flows in DNA microarrays with a particular interest in understanding turbulent mixing. This demonstrates a primary mechanism giving rise to chaotic advection. The linked twist map is a Poincaré section of eggbeatertype flow (Hertzsch et al., 2007) in continuous dynamical systems. The Poincaré section captures the behavior of the flow by viewing a particle’s location at discrete time intervals. The linked twist map is given by the discrete dynamical system
where is a positive parameter. For some values of , the orbits are dense in the domain. However, for other parameter values, voids form. In either case, the truncated orbits exhibit complex structure.
For this experiment, we choose a set of parameter values, 2.5, 3.5, 4.0, 4.1 and 4.3, which produce a variety of orbit patterns. For each parameter value, 50 randomlychosen initial conditions are selected, and 1000 iterations of the linked twist map are used to generate point clouds in . Figure 4 shows examples of typical orbits generated for each parameter value. The goal is to classify the trials by parameter value using PIs to capitalize on distinguishing topological features of the data. We use resolution and a Gaussian with variance to generate the PIs. These parameters were chosen after a preliminary parameter search and classification effort. Similar results hold for a range of PI parameter values.
For a fixed parameter value and a large number of points (many thousands), the patterns in the distributions of iterates show only small visible variations for different choices of the initial condition . However, with few points, such as in Figure 5, there are more significant variations in the patterns for different choices of initial conditions, making classification more difficult.
We perform classification and crossvalidation with a discriminant subspace ensemble. This ML algorithm trains many “weak” learners on randomly chosen subspaces of the data (of a fixed dimension), and classifies and assigns a score to each point based on the current subspace. The final classification arises from an average of the scores of each data point over all learners (Ho, 1998). We perform 10 trials and average the classification accuracies. For the concatenated and PIs, this method achieves a classification accuracy of 82.5%; compared to 49.8% when using only PIs, and 65.7% using PIs. This experiment highlights two strengths of PIs: they offer flexibility in choosing a ML algorithm that is well suited to the data under consideration, and homological information from multiple dimensions may be leveraged simultaneously for greater discriminatory power.
This application is a brief example of the utility of PIs in classification of data from dynamical systems and modeling realworld phenomena, which provides a promising direction for further applications of PIs.
6.4.2 A partial differential equation
The KuramotoSivashinsky (KS) equation is a partial differential equation for a function of spatial variables and time that has been independently derived in a variety of problems involving pattern formation in extended systems driven far from equilibrium. Applications involving surface dynamics include surface nanopatterning by ionbeam erosion (Cuerno and Barabási, 1995; Motta et al., 2012), epitaxial growth (Villain, 1991; Wolf, 1991; Rost and Krug, 1995), and solidification from a melt (Golovin and Davis, 1998). In these applications, the nonlinear term in the KS equation may be anisotropic, resulting in the anisotropic KuramotoSivashinsky (aKS) equation
(4) 
where , and the real parameter controls the degree of anisotropy. At a fixed time , is a patterned surface (periodic in both and ) defined over the plane. Visibly, the anisotropy appears as a slight tendency for the pattern to be elongated in the vertical or horizontal direction.
Numerical simulations of the aKS equation for a range of parameter values (columns) and simulation times (rows) are shown in Figure 6. For all simulations, the initial conditions were lowamplitude white noise. We employed a Fourier spectral method with periodic boundary conditions on a spatial grid, with a fourthorder exponential time differencing RungeKutta method for the time stepping. Five values for the parameter were chosen, namely , 1.25, 1.5, 1.75 and 2, and thirty trials were performed for each parameter value. Figure 7 shows the similarity between surfaces associated to two parameter values and at an early time.
We aim to identify the anisotropy parameter for each simulation using snapshots of surfaces as they evolve in time. Inference of the parameter using the surface alone proves difficult for several reasons. First, Equation (4) exhibits sensitivity to initial conditions: initially nearby solutions diverge quickly. Second, although the surface at a fixed time is an approximation due to the finite discretization of its domain, the spatial resolution is still very large: in fact, these surfaces may be thought of as points in . We were unable to perform standard classification techniques in this space. It was therefore necessary to perform some sort of dimension reduction. One such method is to simply ‘resize’ the surface by coarsening the discretization of the spatial domain after computing the simulation at a high resolution by replacing a block of grid elements with their average surface height. The surfaces were resized in this way to a resolution of and a subspace discriminant ensemble was used to perform classification. Unsurprisingly, this method performs very poorly at all times (first row of Table 2).
The anisotropy parameter also influences the mean and amplitude of the surface pattern. We eliminate differences in the mean by meancentering each surface after the simulation. To assess the impact of the variance of surface height on our task, classification was performed using a normal distributionbased classifier built on the variances of the surface heights. In this classifier, a normal distribution was fit to a training set of 2/3 of the variances for each parameter value, and the testing data was classified based on a test for each of the different models. That is, a value for each new variance was computed for membership to the five normal distributions (corresponding to the five parameter choices of ), and the surface was classified based on the model yielding the highest value. After the pattern has more fully emerged (by, say, time ) this method of classification yields 75 accuracy^{11}^{11}11Accuracy reported is averaged over 100 different training and testing partitions., as shown in Table 2. However, early on in the formation of the pattern, this classifier performs very poorly because height variance is not yet a discriminating feature. Figure 8 shows the normal distribution fit to the variance of the surfaces for each parameter value at times and 5, and illustrates why the variance of surface height is informative only after a surface is allowed to evolve for a sufficiently long time.
Variance of a surface is reflected in its sublevel set filtration (see §A.3 for more details) PD. Yet, the PD and the subsequent PI contain additional topological structure, which may reveal other influences of the anisotropy parameter on the evolution of the surface. Persistence diagrams were computed using the sublevel set filtration, and PIs were generated with resolution and a Gaussian with variance . We think of our pipeline to a PI as a dimensionality reduction in this case, taking a surface which in actuality is a very highdimensional point and producing a much lower dimensional one that retains meaningful characteristics of the original surface.
We again use a subspace discriminant ensemble to classify PIs by parameter. Table 2 compares these results to the same technique applied to low dimensional approximations of the raw surfaces and the normal distributionbased classifier built from surface variance alone. At each time in the system evolution, the best classification accuracy results from using PIs, improving accuracies over using either low resolution approximations of the surfaces or variance of surface height alone by at least 20, including at early times in the evolution of the surface when pattern amplitudes are not visibly differentiated (see Figure 7). We postulate that PIs capture more subtle topological information that is useful for identifying the parameter used to generate each surface.
As we observed in §6.4.1, concatenating and PIs can notably improve the classification accuracy over either feature vector individually. We again note that classification accuracy appears insensitive to the PI parameters. For example, when the variance of the Guassians used to generate the PIs was varied from 0.0001 to 0.1, the classification accuracy of the H PIs, changed by less than one percentage point. The classification accuracy for H fluctuated in a range of approximately five points. For a fixed variance, when the resolution of the image was varied from 5 to 20, the H accuracy varied by little more than three points until the accuracy dropped by six points for a resolution of 25.
PIs performed remarkably well in this classification task, allowing one to capitalize on subtle structural differences in the patterns and significantly reduce the dimension of the data for classification. There is more to be explored in the realm of pattern formation and persistence that is outside the scope of this paper.
Classification Approach  Time t=3  Time t=5  Time t=10 

Subspace Discriminant Ensemble, Resized Surfaces  26.0 %  19.3%  19.3 % 
Variance Normal Distribution Classifier  20.74%  75.2%  77.62 % 
Subspace Discriminant Ensemble, PIs  58.3 %  96.0 %  94.7 % 
Subspace Discriminant Ensemble, PIs  67.7 %  87.3 %  93.3% 
Subspace Discriminant Ensemble, and PIs  72.7 %  95.3 %  97.3 % 
7 Conclusion
PIs offer a stable representation of the topological characteristics captured by a PD. Through this vectorization, we open the door to a myriad of ML tools. This serves as a vital bridge between the fields of ML and topological data analysis and enables one to capitalize on topological structure (even in multiple homological dimensions) in the classification of data.
We have shown PIs yield improved classification accuracy over PLs and PDs on sampled data of common topological spaces at multiple noise levels using medoids. Additionally, computing distances between PIs requires significantly less computation time compared to computing distances between PDs, and comparable computation times with PLs. Through PIs, we have gained access to a wide variety of ML tools, such as SSVM which can be used for feature selection. Features (pixels) selected as discriminatory in a PI are interpretable because they correspond to regions of a PD. We have explored datasets derived from dynamical systems and illustrated that topological information of solutions can be used for inference of parameters since PIs encapsulate this information in a form amenable to ML tools, resulting in high accuracy rates for data that is difficult to classify.
The classification accuracy is robust to the choice of parameters for building PIs, providing evidence that it is not necessary to perform largescale parameter searches to achieve reasonable classification accuracy. This indicates the utility of PIs even when there is not prior knowledge of the underlying data (i.e. high noise level, expected holes, etc.). The flexibility of PIs allows for customization tailored to a wide variety of realworld data sets.
Acknowledgments: We would like to acknowledge the research group of Paul Bendich at Duke University for allowing us access to a persistent homology package which greatly reduced computational time and made analysis of large point clouds feasible. This code can be accessed via GitLab after submitting a request to Paul Bendich. This research is partially supported by the National Science Foundation under Grants No. DMS1228308, DMS1322508, NSF DMS1115668, NSF DMS1412674, and DMR1305449 as well as the DODUSAF under Award Number FA95501210408.
References
 Adcock et al. (2013) Aaron Adcock, Erik Carlsson, and Gunnar Carlsson. The ring of algebraic functions on persistence bar codes. arXiv preprint arXiv:1304.0530, 2013.
 Bendich (2009) Paul Bendich. Analyzing Stratified Spaces Using Persistent Versions of Intersection and Local Homology. PhD thesis, Duke University, 2009.
 Bendich et al. (2014) Paul Bendich, Sang Chin, Jesse Clarke, Jonathan deSena, John Harer, Elizabeth Munch, Andrew Newman, David Porter, David Rouse, Nate Strawn, and Adam Watkins. Topological and statistical behavior classifiers for tracking applications. arXiv preprint arXiv:1406.0214, 2014.
 Bendich et al. (2015) Paul Bendich, JS Marron, Ezra Miller, Alex Pieloch, and Sean Skwerer. Persistent homology analysis of brain artery trees. Annals of Applied Statistics, 2015. To appear.
 Bradley and Mangasarian (1998) Paul S Bradley and Olvi L Mangasarian. Feature selection via concave minimization and support vector machines. In Machine Learning Proceedings of the Fifteenth International Conference, ICML 1998, pages 82–90, 1998.
 Bubenik (2015) Peter Bubenik. Statistical topological data analysis using persistence landscapes. The Journal of Machine Learning Research, 16(1):77–102, 2015.
 Bubenik and Dlotko (2016) Peter Bubenik and Pawel Dlotko. A persistence landscapes toolbox for topological statistics. Jounral of Symbolic Computations, 2016. Accepted.
 Carlsson (2009) Gunnar Carlsson. Topology and data. Bulletin of the American Mathematical Society, 46(2):255–308, 2009.
 Carrière et al. (2015) Mathieu Carrière, Steve Y Oudot, and Maks Ovsjanikov. Stable topological signatures for points on 3d shapes. In Computer Graphics Forum, volume 34, pages 1–12, 2015.
 Chazal et al. (2014) Frédéric Chazal, Vin de Silva, and Steve Oudot. Persistence stability for geometric complexes. Geometriae Dedicata, 173(1):193–214, 2014.
 Chen et al. (2015) YenChi Chen, Daren Wang, Alessandro Rinaldo, and Larry Wasserman. Statistical analysis of persistence intensity functions. arXiv preprint arXiv:1510.02502, 2015.
 Chepushtanova et al. (2014) Sofya Chepushtanova, Christopher Gittins, and Michael Kirby. Band selection in hyperspectral imagery using sparse support vector machines. In Proceedings SPIE DSS 2014, volume 9088, pages 90881F–90881F15, 2014.
 Chung et al. (2009) Moo K Chung, Peter Bubenik, and Peter T Kim. Persistence diagrams of cortical surface data. In Information Processing in Medical Imaging, pages 386–397. Springer, 2009.
 CohenSteiner et al. (2007) David CohenSteiner, Herbert Edelsbrunner, and John Harer. Stability of persistence diagrams. Discrete & Computational Geometry, 37(1):103–120, 2007.
 CohenSteiner et al. (2010) David CohenSteiner, Herbert Edelsbrunner, John Harer, and Yuriy Mileyko. Lipschitz functions have stable persistence. Foundations of computational mathematics, 10(2):127–139, 2010.
 Cuerno and Barabási (1995) Rodolfo Cuerno and AlbertLásló Barabási. Dynamic scaling of ionsputtered surfaces. Physical Review Letters, 74:4746, 1995.
 Dabaghian et al. (2012) Yu Dabaghian, Facundo Memoli, L Frank, and Gunnar Carlsson. A topological paradigm for hippocampal spatial map formation using persistent homology. PLoS Computational Biology, 8(8):e1002581, 2012.
 Di Fabio and Ferri (2015) Barbara Di Fabio and Massimo Ferri. Comparing persistence diagrams through complex vectors. In International Conference on Image Analysis and Processing 2015 Part I; Editors V. Murino, E. Puppo, LNCS 9279, pages 294–305, 2015.
 Donatini et al. (1998) Pietro Donatini, Patrizio Frosini, and Alberto Lovato. Size functions for signature recognition. In SPIE’s International Symposium on Optical Science, Engineering, and Instrumentation, pages 178–183, 1998.
 Edelsbrunner and Harer (2010) Herbert Edelsbrunner and John Harer. Computational topology: An introduction. American Mathematical Society, 2010.
 Edelsbrunner and Harer (2008) Herbert Edelsbrunner and John Harer. Persistent homology – a survey. Contemporary Mathematics, 453:257–282, 2008.
 Edelsbrunner et al. (2012) Herbert Edelsbrunner, A Ivanov, and R Karasev. Current open problems in discrete and computational geometry. Modelirovanie i Analiz Informats. Sistem, 19(5):5–17, 2012.
 Emerson et al. (2015) Tegan Emerson, Michael Kirby, Kelly Bethel, Anand Kolatkar, Madelyn Luttgen, Stephen O’Hara, Paul Newton, and Peter Kuhn. Fourierring descriptor to characterize rare circulating cells from images generated using immunofluorescence microscopy. Computerized Medical Imaging and Graphics, 40:70–87, 2015.
 Ferri and Landi (1999) Massimo Ferri and Claudia Landi. Representing size functions by complex polynomials. Proc. Math. Met. in Pattern Recognition, 9:16–19, 1999.
 Ferri et al. (1997) Massimo Ferri, Patrizio Frosini, Alberto Lovato, and Chiara Zambelli. Point selection: A new comparison scheme for size functions (with an application to monogram recognition). In Computer Vision ACCV’98, pages 329–337. Springer, 1997.
 Fung and Mangasarian (2004) Glenn M. Fung and O.L. Mangasarian. A feature selection newton method for support vector machine classification. Computational Optimization and Applications, 28(2):185–202, 2004.
 Ghrist (2008) Robert Ghrist. Barcodes: The persistent topology of data. Bulletin of the American Mathematical Society, 45(1):61–75, 2008.
 Golovin and Davis (1998) Alexander A. Golovin and Stephen H Davis. Effect of anisotropy on morphological instability in the freezing of a hypercooled melt. Physica D: Nonlinear Phenomena, 116:363–391, 1998.
 Guo et al. (2010) Zhenhua Guo, Lei Zhang, and David Zhang. A completed modeling of local binary pattern operator for texture classification. IEEE Transactions on Image Processing, 19(6):1657–1663, 2010.
 Hatcher (2002) Allen Hatcher. Algebraic Topology. Cambridge University Press, 2002.
 Heath et al. (2010) Kyle Heath, Natasha Gelfand, Maks Ovsjanikov, Mridul Aanjaneya, and Leonidas J Guibas. Image webs: Computing and exploiting connectivity in image collections. In Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, pages 3432–3439. IEEE, 2010.
 Hertzsch et al. (2007) JanMartin Hertzsch, Rob Sturman, and Stephen Wiggins. DNA microarrays: Design principles for maximizing ergodic, chaotic mixing. Small, 3(2):202–218, 2007.
 Ho (1998) Tin Kam Ho. The random subspace method for constructing decision forests. IEEE Transactions on Pattern Analysis and Machine Intelligence, 20(8):832–844, 1998.
 Kaufman and Rousseeuw (1987) Leonard Kaufman and Peter Rousseeuw. Clustering by means of medoids. NorthHolland, 1987.
 Kerber et al. (2016) Michael Kerber, Dmitriy Morozov, and Arnur Nigmetov. Geometry helps to compare persistence diagrams. Proceedings of the Workshop on Algorithm Engineering and Experiments, 2016. Accepted.
 Mileyko et al. (2011) Yuriy Mileyko, Sayan Mukherjee, and John Harer. Probability measures on the space of persistence diagrams. Inverse Problems, 27(12):124007, 2011.
 Motta et al. (2012) Francis C Motta, Patrick D Shipman, and R Mark Bradley. Highly ordered nanoscale surface ripples produced by ion bombardment of binary compounds. Journal of Physics D: Applied Physics, 45(12):122001, 2012.
 Pachauri et al. (2011) Deepti Pachauri, Christian Hinrichs, Moo K Chung, Sterling C Johnson, and Vikas Singh. Topologybased kernels with application to inference problems in Alzheimer’s disease. IEEE Transactions on Medical Imaging, 30(10):1760–1770, 2011.
 Park and Jun (2009) HaeSang Park and ChiHyuck Jun. A simple and fast algorithm for medoids clustering. Expert Systems with Applications, 36(2):3336–3341, 2009.
 Pearson et al. (2015) Daniel A. Pearson, R. Mark Bradley, Francis C. Motta, and Patrick D. Shipman. Producing nanodot arrays with improved hexagonal order by patterning surfaces before ion sputtering. Phys. Rev. E, 92:062401, Dec 2015. doi: 10.1103/PhysRevE.92.062401. URL http://link.aps.org/doi/10.1103/PhysRevE.92.062401.
 Perea and Harer (2013) Jose A Perea and John Harer. Sliding windows and persistence: An application of topological methods to signal analysis. Foundations of Computational Mathematics, pages 1–40, 2013.
 Reininghaus et al. (2015) Jan Reininghaus, Stefan Huber, Ulrich Bauer, and Roland Kwitt. A stable multiscale kernel for topological machine learning. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4741–4748, 2015.
 Rost and Krug (1995) Martin Rost and Joachim Krug. Anisotropic kuramotosivashinsky equation for surface growth and erosion. Physical Review Letters, 75:3894, 1995.
 Singh et al. (2008) Gurjeet Singh, Facundo Memoli, Tigran Ishkhanov, Guillermo Sapiro, Gunnar Carlsson, and Dario L Ringach. Topological analysis of population activity in visual cortex. Journal of Vision, 8(8):11, 2008.
 Topaz et al. (2015) Chad M Topaz, Lori Ziegelmeier, and Tom Halverson. Topological data analysis of biological aggregation models. PloS One, 10(5):e0126383, 2015.
 Turner et al. (2014) Katharine Turner, Yuriy Mileyko, Sayan Mukherjee, and John Harer. Fréchet means for distributions of persistence diagrams. Discrete & Computational Geometry, 52(1):44–70, 2014.
 Verovšek (2016) Sara Kališnik Verovšek. Tropical coordinates on the space of persistence barcodes. arXiv preprint arXiv:1604.00113, 2016.
 Villain (1991) J Villain. Continuum models of crystal growth from atomic beams with and without desorption. J. Phys. I France, 1:19–42, 1991.
 Wolf (1991) Dietrich E Wolf. Kinetic roughening of vicinal surfaces. Physical Review Letters, 67:1783, 1991.
 Zeppelzauer et al. (2016) Matthias Zeppelzauer, Bartosz Zieliński, Mateusz Juda, and Markus Seidl. Topological descriptors for 3d surface analysis. arXiv preprint arXiv:1601.06057, 2016.
 Zhang and Zhou (2010) Li Zhang and Weida Zhou. On the sparseness of 1norm support vector machines. Neural Networks, 23(3):373–385, 2010.
 Zhu et al. (2004) Ji Zhu, Saharon Rosset, Trevor Hastie, and Rob Tibshirani. 1norm support vector machines. Advances in neural information processing systems, 16(1):49–56, 2004.
 Zomorodian and Carlsson (2005) Afra Zomorodian and Gunnar Carlsson. Computing persistent homology. Discrete & Computational Geometry, 33(2):249–274, 2005.
Appendix A Homology and Data
Homology is an invariant that characterizes the topological properties of a topological space . In particular, homology measures the number of connected components, loops, trapped volumes, and so on of a topological space, and can be used to distinguish distinct spaces from one another. More explicitly, the dimensional holes of a space generate a homology group, . The rank of this group is referred to as the th Betti number, , and counts the number of dimensional holes of . For a comprehensive study of homology, see Hatcher (2002).
a.1 Simplicial Complexes and Homology
Simplicial complexes are one way to define topological spaces combinatorially. More precisely, a simplicial complex consists of vertices (0simplices), edges (1simplices), triangles (2simplices), tetrahedra (3simplices), and higherdimensional simplices (containing vertices), such that

if is a simplex in then contains all lowerdimensional simplices of , and

the nonempty intersection of any two simplices in is a simplex in .
The following setup is necessary for a rigorous definition of (simplicial) homology. To a simplicial complex, one can associate a chain complex of vector spaces over a field (often a finite field for a small prime),
Here, vector space consists of all linear combinations of the simplices of , and has as a basis the set of all simplices. The linear map , known as the boundary operator, maps a simplex to its boundary, a sum of its faces. More formally, the boundary map acts on a simplex by
where is the simplex obtained from by removing vertex . We define two subspaces of : subspace is known as the cycles, and subspace is known as the boundaries. The boundary operator satisfies the property , which implies the inclusion .
Homology seeks to uncover an equivalence class of cycles that enclose a dimensional hole—that is, cycles which are not also boundaries of simplices. To this end, the th order homology is defined as , a quotient of vector spaces. The th Betti number is the dimension of this vector space, and counts the number of independent holes of dimension . More explicitly, counts the number of connected components, the number of loops, the number of trapped volumes, and so on. Betti numbers are a topological invariant, meaning that topologically equivalent spaces have the same Betti number.
a.2 Persistence Diagrams from Point Cloud Data
One way to approximate the topological characteristics of a point cloud dataset is to build a simplicial complex on top of it. Though there are a variety of methods to do so, we restrict attention to the Vietoris–Rips simplicial complex due to its computational tractability (Ghrist, 2008). Given a data set (equipped with a metric) and a scale parameter , the Vietoris–Rips complex has as its set of vertices and has a simplex for every collection of vertices whose pairwise distance is at most . However, it is often not apparent how to choose scale . Selecting too small results in a topological space with a large number of connected components, and selecting too large results in a topological space that is contractible (equivalent to a single point).
The idea of persistent homology is to compute homology at many scales and observe which topological features persist across those scales (Ghrist, 2008; Carlsson, 2009; Edelsbrunner and Harer, 2008). Indeed, if is an increasing sequence of scales, then the corresponding Vietoris–Rips simplicial complexes form a filtered sequence . As varies, so does the homology of , and for any homological dimension we get a sequence of linear maps . Persistent homology tracks the homological features over a range of values of . Those features which persist over a larger range are considered to be true topological characteristics, while shortlived features are often considered as noise.
For each choice of homological dimension , the information measured by persistent homology can be presented as a persistence diagram (PD), a multiset of points in the plane. Each point corresponds to a topological feature that appears (is ‘born’) at scale parameter and which no longer remains (‘dies’) at scale . Since all topological features die after they are born, this is an embedding into the upper half plane, above the diagonal line . Points near the diagonal are considered to be noise while those further from the diagonal represent more robust topological features.
a.3 Persistence Diagrams from Functions
Let be a topological space and let be a realvalued function. One way to understand the behavior of map is to understand the topology of its sublevel sets , where . Indeed, given , one can study map using the persistent homology of the resulting filtration of topological spaces, known as the sublevel set filtration:
If is a simplicial complex, then one can produce an increasing sequence of simplicial complexes using a modification of this procedure called the lower star filtration (Edelsbrunner and Harer, 2010). Similarly, if is a cubical complex (an analogue of a simplicial complex that is instead a union of vertices, edges, squares, cubes, and higherdimensional cubes), then one can produce an increasing sequence of cubical complexes.
In §6.4.2, we study surfaces produced from the KuramotoSivashinsky equation. The domain is discretized into a grid of vertices, i.e. a 2dimensional cubical complex with vertices, horizontal edges, vertical edges, and squares. We produce an increasing sequence of cubical complexes as follows:

A vertex is included at scale if .

An edge is included at scale if both of its vertices are present.

A square is included at scale if all four of its vertices are present.
Our PDs are obtained by taking the persistent homology of this cubical complex sublevel set filtration.
We remark that PDs from point cloud data in §A.2 can be viewed as a specific case of PDs from functions. Indeed, given a data set in some metric space , let be the distance function to set , defined by for all . Note that is the union of the metric balls of radius centered at each point in . For , the persistent homology of
is identical to the persistent homology of a simplicial complex filtration called the Čech complex. Furthermore, the persistent homology of the Vietoris–Rips complex is an approximation of the persistent homology of the Čech complex (Edelsbrunner and Harer, 2010, Section III.2).