{'fraction_non_alphanumeric': 0.03981982845368729, 'fraction_numerical': 0.010398198284536872, 'mean_word_length': 4.011887608069165, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Persistence Diagrams and the Heat Equation Homotopy

Persistence homology is a tool used to measure topological features that are present in data sets and functions. Persistence pairs births and deaths of these features as we iterate through the sublevel sets of the data or function of interest. I am concerned with using persistence to characterize the difference between two functions f, g : M ->R, where M is a topological space. Furthermore, I formulate a homotopy from g to f by applying the heat equation to the difference function g-f. By stacking the persistence diagrams associated with this homotopy, we create a vineyard of curves that connect the points in the diagram for f with the points in the diagram for g. I look at the diagrams where M is a square, a sphere, a torus, and a Klein bottle. Looking at these four topologies, we notice trends (and differences) as the persistence diagrams change with respect to time.


Introduction
Given two functions, we are interested in comparing the respective persistence diagrams. The persistence diagram is a set of points in the Cartesian plane used to describe births and deaths of homology groups as we iterate through sublevel sets of a function. There are many ways that we can compare two diagrams, including matching points between the diagrams. However, finding a meaningful matching (see Section 3) is a difficult task, especially if we are interested in capturing the relationship between the underlying functions.
As stated in the proposal, the goal of this Research Initiation Project is to gain an intuition for persistence and homology, as well as to understand the current state of research in these fields. I accomplished this goal by investigating the following problem: For a 2-manifold M , suppose we have two continuous functions f, g : M → R. We can create the persistence diagrams for f and for g. If we know the relationship between f and g, we can make informed decisions when matching points in the persistence diagrams. In particular, we are interested in the case where f and g are homotopic.
Given a homotopy between f and g, we can create a vineyard of the persistence diagrams. Then, we use the vines in the vineyard to help make an informed matching of the points in the persistence diagrams.
In this paper, we discuss two known methods of matching persistence diagrams, by measuring the bottleneck and the Wasserstein distances. Although stability results exist for matching persistence diagrams by minimizing either the bottleneck or the Wasserstein distance, these matchings are made without consideration of the underlying functions f and g. If we are able to create a continuous deformation of the function f to g, then we can use this additional information to aid in the matching of the points in the persistence diagrams. As a result, the matching obtained will be based on the underlying functions. We look into alternate way of measuring the distance between the persistence diagrams for the functions by assuming that there exists a homotopy between the functions. We create a homotopy that we call the heat equation homotopy and measure distances between the persistence diagrams for f and g by using these homotopies to aid in the pairing of points in the persistence diagrams of f and g. Then, we turn to analyzing an example of the heat equation homotopy and discuss various interesting patterns.

Computational Topology Preliminaries
Observing patterns and features in data sets is a common goal in many disciplines, including biology. Extracting the key features from a noisy data set can be an ambiguous task, and often involves simplifying and finding the best view of the data. Computational topology, and more specifically persistent homology, is a tool used for data analysis. Here, we give a brief review of the necessary background of computational topology, but refer you to [10], [11] and [13] for more details.
2.1. Homology. Let X be a simplicial complex of dimension d. For p ∈ N and p ≤ d, the symbol X p will denote the power set of all p-simplices in X. Each set of X p is called a p-chain. The chain group C p is defined by the set X p under the disjoint union, or symmetric difference, operation. This operation can be interpreted as addition modulo two. The group C p is therefore isomorphic to Z 2 to some non-negative integer power. All algebriac groups in this paper are vector spaces over Z 2 . Consider the boundary homomorphism : that maps the p-chain α ∈ C p to the boundary of α, a chain in C p−1 [13].
The p th homology group of X, denoted H p (X), is defined as the kernel of ∂ p modulo the image of ∂ p+1 : H p (X) = Ker(∂ p )/Im(∂ p+1 ). The kernel of the homomorphism ∂ p is the set of elements in the domain that are evaluated to zero (the empty set) and the image of ∂ p+1 is the set of elements of the form ∂ p+1 (x), where x is in the domain C p+1 : Ker(∂ p ) = {α ∈ C p |∂ p (α) = ∅} and Im(∂ p+1 ) = {α ∈ C p |∃α ′ ∈ C p+1 ∋ ∂ p+1 (α ′ ) = α}.
The p th Betti number, β p , is the rank of the p th homology group of X. By definition, the rank of a group is the (smallest) number of generators needed to define the group up to isomorphism. Since we are concerned with groups with Z 2 coefficients, the rank uniquely defines the group up to isomorphism. For example, the group with three generators is Z 3 2 = Z 2 ⊕ Z 2 ⊕ Z 2 and the group with n generators is Z n 2 . 2.2. Persistent Homology. Now, we define persistent homology for functions from R to R. The complete discussion of the extension of these ideas to higher dimensions is found in [10] and [11]. We present a simplified setting to focus on the relevant concepts, while avoiding the complications that arise in the general setting.
Suppose C is the graph of f : R → R. We can think of f as the height function on C. Now, we characterize the topology of the sublevel set R f s = f −1 ((−∞, s]), and we monitor how the homology groups change as s goes from negative infinity to infinity. The zeroth persistent homology group, denoted H 0 (R f s ), will change whenever s is a local maximum or a local minimum of the function f . The maxima and the minima are where the Betti numbers as well as the homology classes change for functions from R to R. A critical point is defined as the values r ∈ R where the derivative is zero: d dx f (r) = 0. Then, s = f (r) is called the critical value at r [18]. A Morse function defined over a subspace of R is a smooth function, such that no two critical values share a function value, and the second derivative at each critical value is non-zero. For a Morse function f , the critical points are the set of r ∈ R with s = f (r), such that a Betti number changes by only one from R f s−ǫ to R f s+ǫ for every sufficiently small value of ǫ > 0. If the sum of the Betti numbers increases, we call r a positive critical point. If the sum decreases, then r is a negative critical point.  Figure 1. On the left, we see the graph of a function f in R 2 . On the right is the corresponding persistence diagram, Dgm 0 (f ). Each point is drawn with multiplicity one. The birth at the point (7, 1) remains unpaired.
As s increases from negative infinity, we label each new component with the positive critical value that introduces the component. We pair each negative critical value s with the most recently discovered unpaired positive critical value representing the components joined at s. This will lead to the diagram Dgm 0 (f ) as demonstrated in Figure 1. Consider one pair: a positive critical value that was introduced at time s = s 1 and a negative critical value that was introduced at time s = s 2 , where s 1 < s 2 . This pair is represented in Dgm 0 (f ) as the point (s 1 , s 2 ). The persistence of that pair is equal to the difference in function values: s 2 − s 1 . In Figure 1, there are three pairs of points and one positive critical value that remains unmatched. This value represents an essential homology class. It would be paired if we were to consider the extended persistence diagram as presented in [6].

Relationships between Diagrams.
We now turn to looking at two functions. Two functions are called homotopic if there exists a continuous deformation of the first function into the second. We are interested in creating a homotopy between f and g in order to observe how the corresponding persistence diagram changes through time. More importantly, we use the homotopy when matching points in the persistence diagrams. Before proceeding, let us formally define a homotopy and give an example that we will use later. Definition 2.1 (Homotopy). We say that f, g : M → R are homotopic if there exists a continuous function F : M×[0, 1] → R such that F (x, 0) = f (x) and F (x, 1) = g(x), for all x ∈ M. We will denote the homotopy F (x, t) by f t (x).
For f, g : M → R, we have a diagram for each integer p ∈ {0, 1, 2}. Now, choose a value of p and assume that we have a homotopy (not necessarily the straight line homotopy) from g to f . Choose τ time-steps of the homotopy, 0 < t 1 < . . . < t τ = 1. Let t 0 = 0, so that Dgm p (f 0 ) is the initial persistence diagram. At each time t j , we have a persistence diagram Dgm p (f t j ). Moreover, given Dgm p (f t j ), we can compute Dgm p (f t j+1 ) in time linear in the number of simplices of the filtration by using a straight line homotopy between f t j and f t j+1 , as described in [8]. In Section 5.2, we describe how to use information obtained from this computation in order to pair points in the persistence diagrams. Then, we stack the diagrams so that Dgm p (f t ) is drawn at height t in R 3 and connect the points in the diagrams by curves formed by the line segments connecting matched points in consecutive diagrams. The result is a piecewise linear path between points in Dgm p (f 0 ) and points in Dgm p (f 1 ). If we let the time difference between any two consecutive diagrams approach zero, the piecewise linear path becomes a set of continuous curves by the stability result for the straight line homotopy (see Section 3.4). Each curve that traces the path of an off-diagonal points through time is called a vine. The collection of vines is referred to as a vineyard [8,10]. We pair the endpoints of each vine to obtain a matching of the persistence points in Dgm p (g) with the points in Dgm p (f ).

Matching Persistence Diagrams
To match the points in Dgm p (g) and Dgm p (f ) without considering the homotopy, we look at methods of finding matchings on bipartite planar graphs. We begin this section with a few definitions. A matching is a bipartite graph P = (A ∪ B, E) where the vertex sets A and B are disjoint, the edges are between one vertex in A and one vertex in B, and each vertex is incident on at most one edge. A matching is maximal if the addition of any edge would result in a graph that is no longer a matching. A matching is perfect if every vertex is incident upon exactly one edge. In other words, it is a matching where there does not exist an unmatched vertex.
We consider two different methods for measuring the distance between the persistence diagrams, Dgm p (f ) and Dgm p (g). The goal is to match every point in A =Dgm p (f ) to a point in B =Dgm p (g) in order to minimize the cost of the matching. One method for determining this cost is to consider the bipartite matching problem, where for every a ∈ A and b ∈ B, the edge {a, b} has a cost C(a, b). Then, the cost of a perfect matching is the sum (or the maximum) value of the edge costs. The cost for a matching that is not perfect is infinite. To resolve the issue where the number of off-diagonal points in both diagrams is not equal or the diagrams are dissimilar, we allow an off-diagonal point to be matched to a point on the line y = x. The use of the bottleneck and Wasserstein matchings for this purpose is presented in Chapter VIII of [10]. As we will show, each of these methods entails some notion of stability for the persistence diagrams. That is, we can bound the bottleneck and Wasserstein distance between persistence diagrams by the distance between f and g, for some class of functions.
3.1. Two Cost Metrics. Consider the matching problem where we must match the elements of set A ⊆ R 2 with elements of a second set B ⊆ R 2 , and where there is a cost associated with each pair (a, b) ∈ A × B. We will define both the bottleneck and the Wasserstein costs of a matching. Later, we will use these distance metrics as the costs to find the bottleneck and the Wasserstein matchings. Let C(a, b) be the L ∞ distance between the points a and b; that is,  The use of the notation W ∞ to denote the bottleneck distance will become clear in the next section. If |A| = |B| = n, then a maximal matching can be found in O(n 5/2 ) using the Hopcroft-Karp algorithm [15]. If we mimic the thresholding approach of the hungarian method [16], then the bottleneck solution can be found in O(n 5/2 log n). Since A and B are sets of points in the plane, we can improve the computational complexity of determining the matching under the bottleneck distance. Efrat, Itai, and Katz developed a geometric improvement to the Hopcroft-Karp algorithm with a running time of O(n 1.5 log 2 n) [12].

The Wasserstein Matching Criterion.
In the Wasserstein Matching, we seek to minimize the maximum degree q Wasserstein cost over all perfect matchings: Although for small values of q, the bottleneck and the Wasserstein criteria may produce different matchings, if we take the limit as q → ∞, we see that the Wasserstein criterion approaches the bottleneck criterion.
The Hungarian method computes the Wasserstein matching in O(n 4 ) computational complexity [16]. In [19], Vaidya maintains weighted Voronoi diagrams for an O(n 2.5 log n) computation of this matching. Further improvements were made by Agarwal, Efrat, and Sharir [1]. They utilize a data structure that improves the running time to O(n 2+ǫ ) for the min-weight Euclidean matching.

Stability Theorems. A distance metric (and corresponding matching) is stable if a small change in the input sets
A and B produces a small change in the measured distance between the sets. The property of stability is not obvious and sometimes not true. If a matching is stable, however, we can use it to create a vineyard from a smooth homotopy.
Let P ∞ be the matching obtained using the bottleneck criteria; that is, the matching P ∞ is the matching that minimizes the bottleneck distance. Let f, g be tame functions. This means that the homology groups of M f s and of M g s have finite rank for all s. In addition, only a finite number of homology groups are realized as H p (M f s ) or H p (M g s ) [11]. The function difference ||f − g|| ∞ is the maximum difference between the function values: The Stability Theorem for Tame Functions, which gives us that the bottleneck distance W ∞ (Dgm p (f ), Dgm p (g)) is bounded above by ||f − g|| ∞ [5]. The proof of this theorem presented in [8] uses the following result: Stability Result of the Straight Line Homotopy. Given f t (x), the straight line homotopy from g to f , we know that there exists a perfect matching P of the persistence diagrams for f and g such that the bottleneck cost of P is upper bounded by the distance between f and g: We will explain how to obtain this matching P in Section 5.3.
The Wasserstein distance is stable for Lipschitz functions with bounded degree k total persistence. This is proven in [7]. If we relax either of these two conditions, then the Wasserstein distance becomes unstable for two functions f and g where ||f − g|| ∞ ≤ ǫ as ǫ approaches zero.

The Heat Equation
The heat equation is a mathematical description of the dispersion of heat through a region in space. We start with the initial reading (or a guess) of the temperature throughout an enclosed space, say an empty room. After infinite time and given no external changes, the temperature at each point in the space will converge to the average initial temperature.

4.1.
Dispersing the Difference. Let the manifold M be a closed square subset of R 2 . Then, we can think of the functions f, g : M → R as surfaces in R 3 . Now, let u 0 be the difference g − f . If u 0 (x) = 0 for all x ∈ M , then we have f = g. Otherwise, define the average of a function as the integral divided by the area.
Then, we can calculate the average value of u 0 over the domain M by subtraction: avg(u 0 ) = avg(g) − avg(f ). We apply the heat equation to u 0 and obtain u(x, t) where u(x, 0) = u 0 (x) and lim t→∞ u(x, t) = c. For sake of simplicity, we will assume that the average of u 0 vanishes. If this is not the case, we can impose this condition by setting g = g − avg(u 0 ). Now, we can observe the difference u disperse through time until u becomes the zero function. Similarly, we know that f (x) + u(x, t) will go from g to f . Although the value of u(x, t) approaches zero for all x as t increases lim t→∞ u(x, t) = 0, we will stop at a time T when u(x, T ) ∈ (−ǫ, +ǫ) for all x ∈ M and for some ǫ > 0. Then, the function f + u goes from g to a function close to f . Furthermore, the manner in which the function f + u changes is dictated by the heat equation.

The Continuous Heat Equation.
Here, we describe the heat equation as it applies to a continuous function. For more details, please refer to [2,4]. In Section 4.3, we will modify these equations to approximate the solution in the discrete case and we introduce the heat equation homotopy.
Let {b 1 , b 2 } be an orthonormal basis for M. The general form of the heat equation satisfies the following conditions: and the initial condition: If we can solve this partial differential equation (PDE), we obtain u(x, t) defined for all x ∈ M and t ≥ 0. Typically, the heat equation has the additional constraint that u(x, t) is constant with respect to t for all x ∈ ∂M. However, we are interested in the case where u(x, t) is heat conserving. That is, we would like avg(u t 1 ) =avg(u t 2 ) for all t 1 , t 2 . As we will show, in order to obtain this goal, the values on the boundary will reflect the values interior to the boundary. The equation in (1) describes how the temperature changes with respect to time and space. We note that if we impose the condition ∂u ∂t (x, t) = 0, then the heat equation will not change with respect to time, and Equation (1) becomes Laplace's equation, ∆u(x, t) = 0. This is known as the steady-state heat equation and will have a unique solution. The iterative methods that we look at in Section 4.3 aim at finding an approximation of u(x, t) for this problem. The final solution will be constant with respect to time, and so we say it is approaching the steady-state. We are interested in following the behavior of heat equation as it approaches the solution to the steady-state heat equation.

The Discrete Heat Equation Homotopy
. Solving a partial differential equation is not a simple task. Thus, we must defer to numerical methods to estimate this solution, which require spatial and temporal discretization [4]. In the following computations, we use a regular grid decomposition of M = [0, 1] 2 , writing x i = (i − 1)h and y j = (j − 1)h, where h = 1/(n − 1) for some fixed integer n ≥ 2. In the next section, we explain how to apply the heat equation on different topologies.
The first step in creating the heat equation homotopy f t (x) is to compute u t (x) = u(x, t) over the discretized domain. There are three issues that can arise when using the continuous formulation of the heat equation described in Section 4.2.
1. We need to solve the partial differential equation presented in Equations (1) and (2). 2. The partial derivative ∂ 2 ut is defined for all non-negative t, but a homotopy must to be defined for t ∈ [0, 1]. In order to resolve these issues, we apply temporal and spatial discretization as well as scaling. The goal is to obtain a homotopy f t (x) from g to f using the heat equation solution u(x, t). Below, we describe one such resolution; note, however, that other approaches may be taken.

Mathematical Description of the Heat Equation.
Let us recall the steady-state heat equation over R 2 : In this equation, we are using {b 1 , b 2 } as the standard basis vectors for R 2 . At each mesh point x = (i, j), we employ the Taylor polynomial in the variable b 1 to obtain an approximation of the second derivative with respect to b 1 : where h is the spatial step size. Similarly, we also have an approximation for the second derivative with respect to b 2 : To simplify notation, we will now use x i,j to denote (i, j). If we plug (4) and (5) into (3), then we obtain the following equation: Thus, the approximation to the heat equation is made by looking at local neighborhoods for each point x i,j . Figure 2 highlights the four neighbors of the mesh needed to compute an approximation to the heat equation. The neighborhood of a point, Nhd(i, j), is defined to be the set of local neighbors of the point x i,j .
We have n 2 equations of the form presented in (6), one for each point in the n × n mesh. We relabel the mesh points in column-major order in order to use one index instead of two: v {(j−1)n+i} := u(x i,j , t). Then, we may express the n 2 linear equations in matrix-vector form, Av = 0.

4.3.2.
An Iterative Algorithm for Linear Systems of Equations. As we have shown above, solving the discrete heat equation finds a solution v to the linear system of equations Av = 0. Above, we have described how to construct the matrix A, as it is the coefficient matrix for the system of linear equations. In the above description, A is equal to L n , the Poisson matrix of order n. We note here that this matrix is n 2 ×n 2 . We can write L n = D−N , where D is the diagonal matrix 4 · I and N is a matrix with 0's on the diagonal and with only 1 as the non-zero entries of the matrix. Sometimes we refer to D as the valency matrix, since it expresses the degree of each mesh point. The matrix N is symmetric and we call it the neighborhood matrix since the non-zero entries in row i correspond to the neighbors of the mesh point v i [3]. That is, N (i, j) = 1 iff v i and v j are adjacent. The iterative algorithm can be defined by these matrices. We want a solution of the form Av = 0, which means (D − N )v = 0. We can re-write this so that Dv = N v. And, the iterative algorithm can immediately be seen: In the original formulation, this translates to: where the neighborhood Nhd(x) is the neighborhood of x defined by Equation (7). This iterative method is known as Jacobi iteration.

4.3.3.
Creating the Homotopy. We continue the process of obtaining u t+1 from u t until we have reached a halting point, where |u t+1 (x) − u t (x)| ≤ ǫ for all x ∈ M and for some predetermined value of ǫ. We will let T be the maximum time computed in the iterative method. Thus, we have a function u : M × [0, T ] → R. We reparameterize u(x, t) with respect to t in order to obtainũ(x, t) = u(x, t · T ) defined over the domain M × [0, 1]. Then, the heat equation homotopy f t can be defined by the equation: .
By using this homotopy, the initial difference between f and g disperses.
Althoughũ(x, t) approaches a constant function as t approaches T , interesting things can happen along the way. For example, critical values can be created.
provided that these elements are within the domain.
Each vertex can have two, three, or four neighbors, resulting in a system of equations that does not conserve heat. Let V be the set of vertices in M, and define the total heat to be the sum of the heat over the entire domain: x∈V u(x). Then, the total heat is not preserved between iterations. To fix this, let x be its own neighbor for every neighbor that x is missing. Hence, each x ∈ V is the neighbor of four other vertices and has four neighbors itself. For the second step, we have The four edges of the square create the boundary on the square topology. If we identify the boundary edges in pairs, we can then create a surface without boundary. One way to do this is to create a torus from the square. Formally, we define the torus topology as follows:  (Torus Topology). Given an n × n mesh, the neighborhood of (i, j) is: where addition and subtraction are calculated modulo n.
The torus topology can also be described by  two a edges together and the two b edges together, we obtain a manifold without boundary. In this manifold, each vertex v is the neighbor of four other vertices, and has four neighbors that contribute to the new value for that vertex. Similarly, we can define the Klein Bottle and the Spherical topologies as depicted in Figure 5.

4.5.
Convergence. When using iterative algorithms, we worry about how long finding the solution (or something close to the solution) will take. In this section, we discuss one metric that measures the convergence of iterative algorithms. If x (k) = u(x, t k ) is the solution at stage k, then the solution at stage k + 1 is We now see that and for the solution x ∞ : If we subtract (10) from (11), then we obtain The action of M = D −1 N on the initial error determines whether the error of the solution x will increase or decrease. Now, if M has n 2 nondegenerate eigenvalues and n 2 linearly independent eigenvectors, we can write M = V ΣV −1 , where Σ is the diagonal matrix of eigenvalues. Since M k = V Σ k V −1 , we use the spectral norm of the matrix M to determine if the error is compounding or decreasing. The spectral norm, ρ, of a matrix the largest absolute eigenvalue of that matrix. If ρ < 1, then the iterative algorithm converges. If ρ > 1, then the iterative algorithm does not converge. For example, if we are using Jacobi iteration as described in Section 4.3.2, then ρ is very close to one, and thus convergence is very slow [14]. Although this is undesirable behavior for finding a solution to the heat equation, for our purposes, it allows us to more closely examine the behavior of the heat equation.

Vineyards of the Heat Equation Homotopy
In Section 2.3, we defined homotopy and stated that we can stack the persistence diagrams associated with the homotopy to create a vineyard. In this section, we estimate the underlying vineyard for the heat equation homotopy by monitoring the transpositions in the filter.

Turning a Mesh into a Filter.
We begin with an n × n matrix of values. Although the heat equation computes neighbors based on a grid, we will be computing persistence using simplicial complexes. Thus, we triangulate the domain by adding in horizontal, vertical and diagonal edges, as well as the triangles formed by the voids. An n × n mesh gives us n 2 vertices, (n − 1)(3n − 1) edges, and 2(n − 1) 2 faces. In Figure 6, we see the triangulation of a four-by-four mesh. Each simplex σ is defined by the vertices that create it. In addition, the value of a simplex, f (σ), is the maximum function value of those vertices. Now, we order the simplices by the following two rules: 1. If f (σ 1 ) < f (σ 2 ). then σ 1 appears before simplex σ 2 .
2. If τ is a subsimplex of σ, denoted τ ≤ σ, then simplex τ appears before simplex σ. We observe here that the second rule does not contradict the first, because f (τ ) ≤ f (σ) whenever τ is a face of σ. The resulting ordering of the simplices is called a filter. The two rules imply that every initial subsequence of the filter defines a subcomplex of the mesh. Growing this initial subsequence until it equals the entire filter gives a sequence of simplicial complexes called the induced filtration. As with any sorting algorithm, to create a filter of m simplices will take O(m log m) time.
Create a filtration for f 0 (x) as prescribed above. We use the filtration to compute the persistence diagrams, just as we used the sublevel sets in Section 2.2. We progress from one complex in the filtration to another by adding one simplex. The addition of this simplex can be the birth of a new homology class or the death of an existing homology class. After we have iterated through the entire filter, we have completed the computation of the persistence diagrams for time 0.
To compute Dgm p (f t 1 ), we could repeat the same process. However, if we use Dgm p (f t 0 ), we can compute the new diagram in linear time and we can match points in the two diagrams [8].

5.2.1.
Transpositions. Suppose we have two filters on m simplices, such that the filters are identical except that two adjacent simplices have swapped order. Then, we can write the first filter: F 1 : σ 1 , σ 2 , . . . , σ i , σ i+1 , . . . , σ m and the second filter: If we have the persistence diagram for the filter F 1, then we can compute the persistence diagram for the the filter F 2 by performing one transposition. A transposition updates the persistence diagram by changing the pairings of two consecutive simplices, if necessary. In the case that dim(σ i ) = dim(σ i+1 ), the transposition does not affect the pairs in the persistence diagram. Thus, only the transposition of simplices with the same number of vertices can result in a pairing swap. One transposition can swap the births and deaths of at most two points in the persistence diagram. We must distinguish between nested and unnested persistence pairings. We say that two persistence points, (b 1 , d 1 ) and (b 2 , d 2 ), are nested if the birth at b 2 and the death at d 2 occur after b 1 and before d 2 . Assuming all events happen at distinct moments of time, this is equivalent to b 1 < b 2 < d 2 < d 1 . If we transpose σ i and σ i+1 , then the three types of pair-swapping transpositions are: 1. The births of two nested pairs are transposed. In this case, σ i is associated with b 1 in Dgm p (F 1), and with b 2 in Dgm p (F 2). 2. The deaths of two nested pairs are transposed. In this case, σ i is associated with d 2 in Dgm p (F 1), and with d 1 in Dgm p (F 2). 3. The birth of one pair and the death of another, unnested, pair are transposed. In this case, the addition of σ i to the filtration created a death in Dgm p (F 1), but a birth in Dgm p (F 2). Figure 7 illustrates the three types of transpositions that can be made. Although pair swaps of types 1 and 2 involve two persistence points in the same diagram, pair swaps of type 3 involve persistence points in diagrams of two consecutive dimensions. In each of the cases, the transposition results in a swap only if changing the order in which the simplices are added changes the persistence pairing that is made. For a complete algorithm to compute the transpositions, please refer to [8].
From these transpositions, we can create the matching referred to in Stability Result of the Straight Line Homotopy of 3.4. Suppose the transposition of σ i and of σ i+1 resulted in a pairing swap. Then, every persistence point is paired with itself in the matching, except for (b 1 , d 1 ) and (b 2 , d 2 ), whose pairings are swapped.

Sweep Algorithm.
Changing one filtration into another may require more than one transposition. In order to compute the persistence diagram, we first create a topological arrangement and then use a sweep algorithm. In the Cartesian plane, write the filtration of f t j horizontally. Below that, σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 σ 7 σ 4 σ 3 σ 5 σ 6 σ 7 σ 1 σ 2 Figure 8. We create a topological arrangement by connecting drawing a line between the vertices that represent the same simplex. By doing this, we find a finite number of crossings. Each crossing represents one transposition in the filter.
write the filtration of f t j+1 using the same simplex names. Then, we connect like simplices with a single curve (does not need to be straight). An example of this process is given in Figure 8. After this arrangement has been created, an ordering on the transpositions can be found by topologically sweeping the arrangement, as presented in [9].
We compute Dgm p (f t j+1 ) and keep track of the matching by progressing one transposition at a time in the order dictated by the sweep algorithm.

5.3.
Measuring Distance between Persistence Diagrams. Now, we assume that there exists a homotopy between functions f and g. Then, we have a vineyard that connects each point a in A =Dgm p (f ) with a point b in B =Dgm p (g). We will use this pairing as our matching and define a distance metric for it.
Assume a and b are connected by the vine s : [0, 1] → R 3 , as described in Section 5.2. Since the vine was created from a homotopy, we will use the index t to emphasize that s(t) is a persistence point for the function f t (x). The velocity of the vine ∂s ∂t will be integrated on [0, 1] in order to measure the distance traveled between a and b : In order to obtain a distance between persistence diagrams we sum these distances over all vines in the vineyard V : If s(t) is only defined for a discrete set of times t i with 0 ≤ i ≤ T , then we obtain an alternate definition for D f g : When a point a ∈ A enters the diagonal at t < 1, we pair a with the corresponding diagonal point, since the diagonal points can only occur as endpoints of a vine by definition of vine in Section 5.2. Then, we only measure the distance over the interval in which the vine is defined, [0, t). Symmetrically, we can have a diagonal point in A paired with an off-diagonal point in B.
We are interested in understanding how this distance metric compares with the bottleneck and Wasserstein distances, as well as to understand the properties of this distance metric and homotopy matching.

An Example
In this section, we present the results from one example. Although the results from one example cannot be generalized, we are able to capture the different behaviors of the homotopy under the four topologies (square, sphere, torus, and Klein bottle). In the subsequent examples, M is a square closed region of R 2 . The functions f and g are approximated by a 101 × 101 mesh. We obtain the function values from two gray-scale images of a bird and a flower respectively, as shown in Figure 9. The difference u 0 = g − f is computed. From here, we will find a homotopy from u 0 to the zero function. We note that this homotopy differs from the one previously defined, since we do not add f to the heat equation solution. In practice, we found that this homotopy more clearly displays the behavior of the heat equation.
In Figures 10-14, we see several stages of the heat equation homotopy using various topologies of the square. The persistence diagrams shown are combined diagrams for all dimensions. In Figure 15, we look at the diagram for the first step of the homotopy, separated for dimensions p = 0 and p = 1, under the Klein bottle and sphere topologies. We note that the points with the highest persistence appear in the dimension one persistence diagram. This is a property that remains true as the homotopy progresses.
6.1. Analyzing Different Topologies. In the topologies without a boundary (torus, Klein bottle, and sphere), a new feature is created with a relatively high persistence. For example, if we start with a triangular region of high values against a border, as in Figure 16, then there exists a 1-cycle in the sublevel sets. Note, however, that there is not a 1−cycle in any sublevel set of the square topology. Since we have created the four different topologies by gluing the edges of the square together in various ways, different behaviors along these edges can be expected. We keep this difference in mind as we continue to look for commonalities and other differences caused by the adopted topology. The persistence diagrams for the torus and the Klein bottle topologies behave similarly. In most of the graphs in this section, the curves of the torus and the Klein bottle are usually parallel. This is an indication that orientability has little effect on the heat equation homotopy. This does not come as a not a surprise, as we did not use the orientation when computing the heat equation.
6.2. Duration of Vines. In Figure 17, we see the distribution of the vine lengths in the vineyard Dgm 0 (u t ). The histogram is skewed right, since the mean is greater than the median. Table 1 confirms this observation. The same pattern is in fact observed under different topologies. When comparing the vine lengths of any two topologies, the Kolmogorov-Smirnov test for statistical difference with α = .05 fails to reject the null hypothesis that the distributions are the same. On the other hand, if we remove the shortlived vines, then we start to see that comparing the Klein bottle and sphere topologies results in the rejection of the null hypothesis. However, this is not a strong enough indication that these distributions are different. The details of this statistical method are out of the scope of this paper, but can be found in [17].  From these observations we can conclude that all topologies display a similar distribution of the length of vines, with many of the vines being short-lived. In addition, under the torus and Klein bottle topologies, a large number of vines span the entire vineyard.
6.3. Monitoring Total Persistence. The degree q total persistence is the sum of the q th powers of persistence over all points in the persistence diagram. In Figure 18, we have the graph of the degree one total persistence   Figure 19. The declining total persistence of degrees two and three. In the graph for the total persistence of degree three, we omit the sphere topology. Relative to the other topologies, the values were very high. versus the iteration. Although the total persistence rapidly deteriorates initially, the decay slows down around iteration 100. --We notice here   Figure 20 shows that the mean absolute value of the change between steps of the homotopy decreases rapidly at first, then slowly. Recall that the values of the mesh points range from −255 to 255. At the first step, the mean absolute change is between 10 and 15 (not shown in Figure 20), which is only a 2% − 3% initial change. The value, however, remains above 0.1 in all cases except for the square topology. Given the nature of the heat equation, we expect the values to decrease with respect to time. The small values for mean absolute change are in part due to the initial values of u 0 . The values of one vertex does not differ by a large amount from the values of its neighbors. We know that the iterative method chosen is slow to converge; however, this graph allows us to gauge how little change is occurring at each step of the homotopy. 6.5. Counting Transpositions. The total number transpositions between steps of the heat equation versus time is shown in Figure 21(a). All four topologies follow a decreasing pattern that levels off, with the Klein bottle and the torus behaving distinctively different from the square and the sphere.
We remark on several notable observations from these diagrams. The total number of transpositions is still significant when the heat equation algorithm reaches the halting condition. At the last step, the square topology makes over 400,000 total transpositions. The other three topologies make even more transpositions. In Figure 21 (b),(c), and (d) to see a different pattern for the number of the switches resulting from the transpositions. Under the square topology, the number of pair swaps of types 1, 2, and 3 is down to the double digits after 500 iterations of the heat equation. We notice that Figures (b) and (c) are almost identical. This is not a surprising observation, since pairing swaps 1 and 2 (described in Section 5.2.1) are symmetric cases of two births of the same dimension or two deaths of the same dimension being transposed, resulting in a pair swap in the diagram.
Three other patterns are worth noting in the graphs of Figure 21. First, we see that the number of transpositions (of all kinds) rapidly decreases until iteration 25. At this point, the four topologies show different behaviors. Second, after iteration 200, the torus and the Klein bottle topologies have both leveled off to a constant function. Finally, we notice that adding graphs (b), (c), and (d) does not result in a graph that looks like (a). Thus, a significant number of the transpositions made are those that do not result in a pairing swap. Proportionally, this occurs more often in the torus and Klein bottle topologies than in the square and the sphere topologies. In Figure 22, we restrict our counts to the number transpositions of vertices only. The values at the vertices dictates the values of the edges and the faces. The pattern that this graph follows is similar to Figure 21(a), the graph for the number of transpositions generated from all of the simplices Thus, is seems that the behavior of the vertex-vertex transpositions is proportional to the behavior of simplex-simplex transpositions.

Conclusion
The objective of this project was to find a new way of measuring the distance between two functions. We develop a homotopy of functions based on the heat equation and investigated the behavior of the vineyard resulting from this homotopy. We computed this homotopy using four topologies, and found some trends that differed from our expectation. For example, the heat equation on elliptic surfaces is known to converge more quickly than in Euclidean space. Yet, we found the behavior of the sphere topology to contradict this fact. We believe that this contradiction arises from the discretization of the sphere.
We began an investigation of the behavior of the heat equation homotopy; however, many directions still remain for where we can continue. We would like to explore if the observations made in Section 6 were specific to initial function u 0 , or if we are observing behaviors that arise from the topologies used in the heat equation. Moreover, we would like to formally compare the matching we obtain with the bottleneck and Wasserstein matchings.