Approximating Continuous Functions on Persistence Diagrams Using Template Functions

Abstract

The persistence diagram is an increasingly useful tool from Topological Data Analysis, but its use alongside typical machine learning techniques requires mathematical finesse. The most success to date has come from methods that map persistence diagrams into vector spaces, in a way which maximizes the structure preserved. This process is commonly referred to as featurization. In this paper, we describe a mathematical framework for featurization called template functions, and we show that it addresses the problem of approximating continuous functions on compact subsets of the space of persistence diagrams. Specifically, we begin by characterizing relative compactness with respect to the bottleneck distance, and then provide explicit theoretical methods for constructing compact-open dense subsets of continuous functions on persistence diagrams. These dense subsets—obtained via template functions—are leveraged for supervised learning tasks with persistence diagrams. Specifically, we test the method for classification and regression algorithms on several examples including shape data and dynamical systems.

Appendices

Appendix A: Implementation of the Interpolating Polynomials Algorithm

In this appendix, we give more details on the implementation of the interpolating polynomials described in Sect. 6.2. The barycentric formula for Lagrange interpolation described by [7] is given by

$$\begin{aligned}&f(x) := \sum \nolimits _{j=0}^m{\ell ^{\mathcal {A}}_j(x) c_j} = \frac{\sum \nolimits _{j=0}^m{\frac{w_j}{x-{\tilde{x}}_j}}c_j}{\sum \nolimits _{j=0}^k{\frac{w_j}{x-{\tilde{x}}_j}}}; \text { where } \quad \nonumber \\&w_j=\frac{1}{\ell '(a_j)}; \quad \ell '(a_j) = \prod \limits _{i=0, i\ne j}^{m}{(a_j-a_i)}, \end{aligned}$$

(A1)

while ${\mathcal {A}}= \{a_i\}_{i=0}^m \subset {\mathbb {R}}$ is a finite set of distinct mesh values, and $\{c_i \in {\mathbb {R}}\}$ is a collection of evaluation values. The function in Eq. (A1) has the property that $f(a_i) = c_i$ for all i, and it also satisfies the partition of unity condition $\sum \nolimits _{j=0}^{m}{f(x)}= 1, \,\, \forall \, x$.

Barycentric Lagrange interpolation is often used for approximating ${\mathbb {R}}$-valued functions and there are efficient algorithms for obtaining the weights associated with it. However, in our formulation we need to an interpolating polynomial over an ${\mathbb {R}}^2$-valued function. Therefore, we next describe how to expand the algorithm for interpolating a scalar valued function to interpolating a function on the plane. Note that the notation used here is self-contained from Sect. 6.2.

We assume that our planar mesh is the outer product of $m+1$ mesh points along the birth time x-axis, and $n+1$ points along the lifetime y-axis. We also assume that the persistence diagram has N pairs of (birth, lifetime) points.

• 1. Get $\tilde{\gamma }$ and $\phi$ which correspond to the interpolation matrices along the x-mesh and the y-mesh, respectively. These are the matrices that describe the linear transformation from the $m+1$ mesh points of birth times ($n+1$ mesh of lifetimes) to the corresponding interpolated values of the N query birth times (N query lifetimes) for a given diagram. This step is equivalent to separately obtaining the interpolation matrices for the birth times and the lifetimes.

• 2. Set $\gamma =\tilde{\gamma }^T$.

• 3.

  • (a) Replicate each column in $\gamma$ $n+1$ times to obtain $\Gamma$ whose dimensions are $(m+1)\times (N\times (n+1))$.

  • (b) Unravel $\phi$ row-wise into a row vector, then replicate each row $m+1$ times to obtain $\Phi$ whose dimensions are $(m+1)\times (N\times (n+1))$.

• 4. Use element-wise multiplication to obtain $\tilde{\Psi }=\Gamma \cdot \Phi$, where $\cdot$ means element-wise multiplication, and $\tilde{\Psi }$ has dimension $(m+1)\times (N\times (n+1))$.

• 5.

  • (a) Split $\tilde{\Psi }$ into N chunks of $(m+1)\times (n+1)$ matrices along the columns axis.

  • (b) Concatenate the split pieces row-wise to obtain an $(N\times (m+1))\times (n+1)$ matrix $\Psi$.

• 6. Reshape $\Psi$ by concatenating each $(m+1)\times (n+1)$ piece row-wise to obtain an $N \times ((m+1)\times (n+1))$ matrix $\Xi$.

• 7. Let the 2D base mesh be given as

  $$\begin{aligned} \begin{bmatrix} f_{00} &{} f_{01} &{} \ldots &{} f_{0n} \\ f_{10} &{} f_{11} &{} \ldots &{} f_{1n} \\ \vdots &{} &{} &{} \vdots \\ f_{m0} &{} f_{m1} &{} \ldots &{} f_{mn} \end{bmatrix}, \end{aligned}$$

  where $f_{ij} = f(x_i, y_j)$ and $(x_i, y_j)$ is a unique point in the 2D mesh. Define the vector $[f_{00} \, f_{01} \, \ldots \, f_{mn}]$ which is obtained by unraveling the 2D mesh row-wise.

• 8. We can interpolate the query points $(x_q, y_q)$ using

  $$\begin{aligned} p(x_q, y_q) = \begin{bmatrix} \ell _0(x_0) \ell _0(y_0) &{} \ldots &{} \ell _m(x_0) \ell _n(y_0) \\ \ell _0(x_1) \ell _0(y_1) &{} \ldots &{} \ell _m(x_1) \ell _n(y_1) \\ \vdots &{} &{} \vdots \\ \ell _0(x_{N-1}) \ell _0(y_{N-1}) &{} \ldots &{} \ell _m(x_{N-1}) \ell _n(y_{N-1}) \end{bmatrix} \begin{bmatrix} f_{00} \\ f_{01} \\ \vdots \\ f_{mn} \end{bmatrix}. \end{aligned}$$

Here is a sketch of the resulting matrices:

$$\begin{aligned} \tilde{\gamma }&= \begin{bmatrix} \ell _0(x_0) &{} \ell _1(x_0) &{} \ldots &{} \ell _m(x_0) \\ \vdots &{} &{} &{} \vdots \\ \ell _0(x_{N-1}) &{} \ell _1(x_{N-1}) &{} \ldots &{} \ell _m(x_{N-1}) \end{bmatrix}_{N\times (m+1)},\\ \phi&= \begin{bmatrix} \ell _0(y_0) &{} \ell _1(y_0) &{} \ldots &{} \ell _n(y_0) \\ \vdots &{} &{} &{} \vdots \\ \ell _0(y_{N-1}) &{} \ell _1(y_{N-1}) &{} \ldots &{} \ell _n(y_{N-1}) \end{bmatrix}_{N\times (n+1)},\\ \gamma = \tilde{\gamma }^T&= \begin{bmatrix} \ell _0(x_0) &{} \ell _0(x_1) &{} \ldots &{} \ell _0(x_{N-1}) \\ \ell _1(x_0) &{} \ell _1(x_1) &{} \ldots &{} \ell _1(x_{N-1})\\ \vdots &{} &{} &{} \vdots \\ \ell _m(x_0) &{} \ell _m(x_1) &{} \ldots &{} \ell _m(x_{N-1}) \end{bmatrix}_{(m+1)\times N}, \\ \Gamma&= \begin{bmatrix} \ell _0(x_0) &{} \ell _0(x_0) &{} \ldots &{} \ell _0(x_0) &{} \ldots &{} \ell _0(x_{N-1}) &{} \ell _0(x_{N-1}) &{} \ldots &{} \ell _0(x_{N-1})\\ \ell _1(x_0) &{} \ell _1(x_0) &{} \ldots &{} \ell _1(x_0) &{} \ldots &{} \ell _1(x_{N-1}) &{} \ell _1(x_{N-1}) &{} \ldots &{} \ell _1(x_{N-1})\\ \vdots &{} &{} \vdots &{} &{} \vdots &{} &{} \vdots &{} &{} \vdots \\ \ell _m(x_0) &{} \ell _m(x_0) &{} \ldots &{} \ell _m(x_0) &{} \ldots &{} \ell _m(x_{N-1}) &{} \ell _m(x_{N-1}) &{} \ldots &{} \ell _m(x_{N-1}) \end{bmatrix} \end{aligned}$$

where $\Gamma$ has dimension $(m+1)\times (N\times (n+1))$.

$$\begin{aligned} \Phi = \begin{bmatrix} \ell _0(y_0) &{} \ell _1(y_0) &{} \ldots &{} \ell _n(y_0) &{} \ldots &{} \ell _0(y_{N-1}) &{} \ell _1(y_{N-1}) &{} \ldots &{} \ell _n(y_{N-1}) \\ \ell _0(y_0) &{} \ell _1(y_0) &{} \ldots &{} \ell _n(y_0) &{} \ldots &{} \ell _0(y_{N-1}) &{} \ell _1(y_{N-1}) &{} \ldots &{} \ell _n(y_{N-1}) \\ \vdots &{} &{} \vdots &{} &{} \vdots &{} &{} \vdots &{} &{} \vdots \\ \ell _0(y_0) &{} \ell _1(y_0) &{} \ldots &{} \ell _n(y_0) &{} \ldots &{} \ell _0(y_{N-1}) &{} \ell _1(y_{N-1}) &{} \ldots &{} \ell _n(y_{N-1}) \end{bmatrix} \end{aligned}$$

where $\Phi$ has dimension $(m+1)\times (N\times (n+1))$.

We can now compute the elementwise product $\Psi = \Gamma \cdot \Phi$, which has the dimension $(m+1)\times (N\times (n+1))$.

We then need to apply the following operations: (i) reshaping $\Psi$ to obtain ${\hat{\Psi }}_1$ given by

$$\begin{aligned} {\hat{\Psi }}_1 = \begin{bmatrix} \ell _0(x_0)\ell _0(y_0) &{} \ell _0(x_0)\ell _1(y_0) &{} \ldots &{} \ell _0(x_0)\ell _n(y_0) \\ \ell _1(x_0)\ell _0(y_0) &{} \ell _1(x_0)\ell _1(y_0) &{} \ldots &{} \ell _1(x_0)\ell _n(y_0) \\ \vdots &{} &{} \vdots &{} \\ \ell _m(x_0)\ell _0(y_0) &{} \ell _m(x_0)\ell _1(y_0) &{} \ldots &{} \ell _m(x_0)\ell _n(y_0) \\ \vdots &{} &{} \vdots &{} \\ \ell _0(x_{N-1})\ell _0(y_{N-1}) &{} \ell _0(x_{N-1})\ell _1(y_{N-1}) &{} \ldots &{} \ell _0(x_{N-1})\ell _n(y_{N-1}) \\ \ell _1(x_{N-1})\ell _0(y_{N-1}) &{} \ell _1(x_{N-1})\ell _1(y_{N-1}) &{} \ldots &{} \ell _1(x_{N-1})\ell _n(y_{N-1}) \\ \vdots &{} &{} \vdots &{} \\ \ell _m(x_{N-1})\ell _0(y_{N-1}) &{} \ell _m(x_{N-1})\ell _1(y_{N-1}) &{} \ldots &{} \ell _m(x_{N-1})\ell _n(y_{N-1}) \end{bmatrix}. \end{aligned}$$

(ii) unraveling ${\hat{\Psi }}_1$ into an $N\times ((m+1)\times (n+1))$ matrix ${\hat{\Psi }}_2$ given by

$$\begin{aligned} {\hat{\Psi }}_2 = \begin{bmatrix} \ell _0(x_0) \ell _0(y_0) &{} \ldots &{} \ell _0(x_0) \ell _n(y_0) &{} \ldots &{} \ell _m(x_0) \ell _n(y_0) \\ \vdots &{} &{} \vdots &{} \\ \ell _0(x_k) \ell _0(y_k) &{} \ldots &{} \ell _0(x_k) \ell _n(y_k) &{} \ldots &{} \ell _m(x_k) \ell _n(y_k) \\ \vdots &{} &{} \vdots &{} \\ \ell _0(x_{N-1}) \ell _0(y_{N-1}) &{} \ldots &{} \ell _0(x_{N-1}) \ell _n(y_{N-1}) &{} \ldots &{} \ell _m(x_{N-1}) \ell _n(y_{N-1}) \end{bmatrix}.\nonumber \\ \end{aligned}$$

(A2)

The collection of all the scores constitutes the feature vector corresponding to the chosen base mesh point and to the query points where the latter are the persistence diagram points. In this study we summed the rows of ${\hat{\Psi }}_2$ after taking the absolute value of each entry. The resulting number represents the score at each base mesh point. If the persistence diagram contains the mesh points and we want to find the interpolated values at query points $p_{\mathrm{interp}}$, then we would compute $p_{\mathrm{interp.}}={\hat{\Psi }}_2\, f$.

The implementation of this algorithm can be found in the teaspoon package at teaspoon.ML.feature_functions.interp_polynomial.

Appendix B: Additional Shape Data Results

This appendix gives additional results for the SHREC data set described in Sect. 8.4 using tent functions instead of interpolating polynomials. Table 4 should be compared to the results of Table 3.

Table 4 Results of classification of shape data discussed in Sect. 8.4