Persistence weighted Gaussian kernel for topological data analysis

TL;DR

This paper introduces a new kernel method for persistence diagrams in topological data analysis, offering stability, explicit persistence control, and fast approximation, demonstrated on protein and glass data.

Contribution

It proposes a novel kernel for persistence diagrams that enhances stability, control, and computational efficiency in topological data analysis.

Findings

The kernel satisfies stability properties.

It provides explicit control over persistence effects.

Applied to protein and glass data with improved results.

Abstract

Topological data analysis (TDA) is an emerging mathematical concept for characterizing shapes in complex data. In TDA, persistence diagrams are widely recognized as a useful descriptor of data, and can distinguish robust and noisy topological properties. This paper proposes a kernel method on persistence diagrams to develop a statistical framework in TDA. The proposed kernel satisfies the stability property and provides explicit control on the effect of persistence. Furthermore, the method allows a fast approximation technique. The method is applied into practical data on proteins and oxide glasses, and the results show the advantage of our method compared to other relevant methods on persistence diagrams.

Figures

Tables

Table 1. Table 1: Results of SVM with PWGK, PSSK, and Gaussian. Average classification rates ( % percent \% ) for 99 test data sets are shown.

	RKHS-Linear	RKHS-Gauss
PWGK	60.0	83.0
PSSK	49.5	54.5
Gauss	57.6	69.7

Table 2. Table 2: CV classification rates ( % percent \% ) of SVM with PWGK and MTF (cited from [ CMW + 15 ] ).

	Protein-Drug	Hemoglobin
PWGK	100	88.90
MTF-SVM	(nbd) 93.91 / (bd) 98.31	84.50

Equations

$$d_B(D,E):=\underset{\gamma }{\inf }\underset{x\in D}{\sup }{\Vert x-\gamma (x)\Vert }_{\infty },$$

$$d_H(X,Y):=\max \left\{\underset{\mathbf{x}\in X}{\sup }\underset{\mathbf{y}\in Y}{\inf }{d}_M(\mathbf{x},\mathbf{y}),\underset{\mathbf{y}\in Y}{\sup }\underset{\mathbf{x}\in X}{\inf }{d}_M(\mathbf{x},\mathbf{y})\right\}.$$

$$d_B(D_q(X),D_q(Y))\le d_H(X,Y).$$

$$E_k:M_b(\Omega )\to {\mathcal{H}}_k,\mu \mapsto \int k(\cdot ,x)d\mu (x).$$

$$d_k(\mu ,\nu )={\Vert E_k(\mu )-E_k(\nu )\Vert }_{{\mathcal{H}}_k}$$

$${\mu }_D^w\mapsto E_k({\mu }_D^w):=\sum _{x\in D}w(x)k(\cdot ,x).$$

$$k^w(x,y):=w(x)w(y)k(x,y)$$

$$E_{k^w}:{\mu }_D\mapsto \sum _{x\in D}w(x)w(\cdot )k(\cdot ,x)\in {\mathcal{H}}_{k^w}.$$

$${\mathcal{H}}_k\to {\mathcal{H}}_{k^w},f\mapsto wf$$

$$⟨E_k({\mu }_D^w),E_k({\mu }_E^w){⟩}_{{\mathcal{H}}_k}=⟨E_{k^w}({\mu }_D),E_{k^w}({\mu }_E){⟩}_{{\mathcal{H}}_{k^w}},$$

$$k_{PWG}(x,y)=w_{\mathrm{arc}}(x)w_{\mathrm{arc}}(y)e^{-\frac{{\Vert x-y\Vert }^2}{2{\sigma }^2}}.$$

$$d_{k_G}^{w_{\mathrm{arc}}}(D_q(X),D_q(Y))\le L(M,d;C,p,\sigma )d_H(X,Y),$$

$${\mathcal{P}}_{\mathrm{finite}}(M)\to {\mathcal{H}}_{k_G},X\mapsto E_{k_G}({\mu }_{D_q(X)}^{w_{\mathrm{arc}}})$$

$$K_L(D,E)$$

$$=\sum _{x\in \underline{D}}\sum _{y\in \underline{E}}{w}_{\mathrm{arc}}(x)w_{\mathrm{arc}}(y)k_G(x,y)$$

$$K_G(D,E)=\exp \left(-\frac{d_{k_G}^{w_{\mathrm{arc}}}(D,E{)}^2}{2{\tau }^2}\right),$$

$$d_{k_G}^{w_{\mathrm{arc}}}(D,E{)}^2$$

$$=\sum _{x\in \underline{D}}\sum _{x^{\prime }\in \underline{D}}{w}_{\mathrm{arc}}(x)w_{\mathrm{arc}}(x^{\prime })k_G(x,x^{\prime })$$

$$+\sum _{y\in \underline{E}}\sum _{y^{\prime }\in \underline{E}}{w}_{\mathrm{arc}}(y)w_{\mathrm{arc}}(y^{\prime })k_G(y,y^{\prime })$$

$$-2\sum _{x\in \underline{D}}\sum _{y\in \underline{E}}{w}_{\mathrm{arc}}(x)w_{\mathrm{arc}}(y)k_G(x,y).$$

$$K_{PSS}(D,E)$$

$$=\frac{1}{8\pi t}\sum _{x\in \underline{D}}\sum _{y\in \underline{E}}{e}^{-\frac{{\Vert x-y\Vert }^2}{8t}}-e^{-\frac{{\Vert x-\bar{y}\Vert }^2}{8t}},$$

$$S=\{\{1\},\{2\},\{3\},\{1,2\},\{1,3\},\{2,3\},\{1,2,3\}\},$$

$$T=\{\{1\},\{2\},\{3\},\{1,2\},\{1,3\},\{2,3\}\}.$$

$$C_q(S)={\mathrm{Span}}_K\{⟨\sigma ⟩\mid \sigma \in S_q\}$$

$${\partial }_q⟨i_0,\dots ,i_q⟩=\sum _{\ell =0}^q(-1{)}^{\ell }⟨i_0,\dots ,\widehat{i_{\ell }},\dots ,i_q⟩,$$

$$\displaystyle M_{2}=\left[\begin{array}[]{r}1\\ 1\\ -1\end{array}\right],\quad M_{1}=\left[\begin{array}[]{rrr}-1&0&-1\\ 1&-1&0\\ 0&1&1\end{array}\right]$$

$$H_q(S)=Z_q(S)/B_q(S).$$

$$Z_{1}(S)={\rm Span}_{K}[\begin{array}[]{ccc}1&1&-1\end{array}]^{T}=B_{1}(S)$$

$${\Delta }_q=\left\{\sum _{i=0}^q{\lambda }_i{e}_i|\sum _{i=0}^q{\lambda }_i=1,{\lambda }_i\ge 0\right\},$$

$${\Delta }_q^{\ell }=\left\{\sum _{i=0}^q{\lambda }_i{e}_i|\sum _{i=0}^q{\lambda }_i=1,{\lambda }_i\ge 0,{\lambda }_{\ell }=0\right\}.$$

$$C_q(X)={\mathrm{Span}}_K\{\sigma \mid \sigma \in X_q\}.$$

$${\partial }_q\sigma =\sum _{\ell =0}^q(-1{)}^{\ell }\sigma \circ {\iota }_q^{\ell }.$$

$$H_q(X)=Z_q(X)/B_q(X).$$

Full text

1 Introduction

Recent years have witnessed an increasing interest in utilizing methods of algebraic topology for statistical data analysis. This line of research is called topological data analysis (TDA) [Car09], which has been successfully applied to various areas including information science [CIdSZ08, dSG07], biology [KZP*+*07, XW14], brain science [LCK*+*11, PET*+*14, SMI*+*08], biochemistry [GHI*+*13], and material science [NHH*+*15a, NHH*+*15b]. In many of these applications, it is not straightforward to provide feature vectors or descriptors of data from their complicated geometric configurations. The aim of TDA is to detect informative topological properties (e.g., connected components, rings, and cavities) from such data, and use them as descriptors.

A key mathematical apparatus in TDA is persistent homology, which is an algebraic method for extracting robust topological information from data. To provide some intuition for the persistent homology, let us consider a typical way of constructing persistent homology from data points in a Euclidean space, assuming that the data lie on a submanifold. The aim is to make inference on the topology of the underlying manifold from finite data. We consider the $r$-balls (balls with radius $r$) to recover the topology of the manifold, as popularly employed in constructing an $r$-neighbor graph in many manifold learning algorithms. While it is expected that, with an appropriate choice of $r$, the $r$-ball model can represent the underlying topological structures of the manifold, it is also known that the result is sensitive to the choice of $r$. If $r$ is too small, the union of $r$-balls consists simply of the disjoint $r$-balls. On the other hand, if $r$ is too large, the union becomes a contractible space. Persistent homology [ELZ02] can consider all $r$ simultaneously, and provides an algebraic expression of topological properties together with their persistence over $r$. We give a brief explanation of persistent homology in Supplementary material A.3.

The persistent homology can be visualized in a compact form called a persistence diagram $D=\{(b_{i},d_{i})\in\mathbb{R}^{2}\mid i\in I,\ b_{i}\leq d_{i}\}$, and this paper focuses on persistence diagrams, since the contributions of this paper can be fully explained in terms of persistence diagrams. Every point $(b_{i},d_{i})\in D$, called a generator of the persistent homology, represents a topological property (e.g., connected components, rings, and cavities) which appears at $X_{b_{i}}$ and disappears at $X_{d_{i}}$ in the $r$-ball model. Then, the persistence $d_{i}-b_{i}$ of the generator shows the robustness of the topological property under the radius parameter. As an example shown in Figure 1, the rings $\alpha_{1},\alpha_{2}$ and other tiny ones are expressed as $x_{1},x_{2}$, and the other points in the persistence diagram shown in Figure 10. A topological property with large persistence can be regarded as a reliable structure, while that with small persistence (points close to the diagonal) is likely to be noise. In this way, persistence diagrams encode topological and geometric information of data points.

While persistence diagrams nowadays start to be applied to various problems such as the ones listed in the beginning of this section, statistical or machine learning methods for analysis on persistence diagrams are still limited. In TDA, analysts often elaborate only one persistence diagram and, in particular, methods for handling many persistence diagrams, which should be supposed to contain randomness from the data, are at the beginning stage (see the end of this section for related works). Hence, developing a statistical framework on persistence diagrams is a significant issue for further success of TDA.

To this aim, this paper discusses kernel methods for persistence diagrams (see Figure 3). Since a persistence diagram is a point set of variable size, it is not straightforward to apply standard methods of statistical data analysis, which typically assume vectorial data. Here, to vectorize persistence diagrams, we employ the framework of kernel embedding of (probability and more general) measures into reproducing kernel Hilbert spaces (RKHS). This framework has recently been developed, leading various new methods for nonparametric inference [MFDS12, SGSS07, SFG13]. It is known [SFL11] that, with an appropriate choice of kernels, a signed measure can be uniquely represented by the Bochner integral of the feature vectors with respect to the measure. Since a persistence diagram can be regarded as a non-negative measure, it can be embedded into an RKHS by the Bochner integral. Once such a vector representation is obtained, we can introduce any kernel methods for persistence diagrams systematically.

For embedding persistence diagrams in an RKHS, we propose a useful class of positive definite kernels, called persistence weighted Gaussian kernel (PWGK). It is important that the PWGK can discount the contributions of generators close to the diagonal (small persistence), since in many applications those generators are likely to be noise. The advantages of this kernel are as follows. (i) We can explicitly control the effect of persistence, and hence, discount the noisy generators appropriately in statistical analysis. (ii) As a theoretical contribution, the distance defined by the RKHS norm for the PWGK satisfies the stability property, which ensures the continuity from data to the vector representation of the persistence diagram. (iii) The PWGK allows efficient computation by using the random Fourier features [RR07], and thus it is applicable to persistence diagrams with a large number of generators, which are seen in practical examples (Section 4).

We demonstrate the performance of the proposed kernel method with synthesized and real-world data, including protein datasets (taken by NMR and X-ray crystallography experiments) and oxide glasses (taken by molecular dynamics simulations). We remark that these real-world problems have biochemical and physical significance in their own right, as detailed in Section 4.

There are already some relevant works on statistical approaches to persistence diagrams. Some studies discuss how to transform a persistence diagram to a vector [Bub15, CMW*+*15, COO15, RT15]. In these methods, a transformed vector is typically expressed in a Euclidean space ${\mathbb{R}}^{k}$ or a function space $L^{p}$, and simple and ad-hoc summary statistics like means and variances are used for data analysis such as principal component analysis and support vector machines. The most relevant to our method is [RHBK15] (see also [KHN*+*15]), where they vectorize a persistence diagram by using the difference of two Gaussian kernels evaluated at symmetric points with respect to the diagonal so that it vanishes on the diagonal. We will show detailed comparisons between this method and ours. Additionally, there are some works discussing statistical properties of persistence diagrams for random data points: [CGLM14] show convergence rates of persistence diagram estimation, and [FLR*+*14] discuss confidence sets in a persistence diagram. These works consider a different but important direction to the statistical methods for persistence diagrams.

The remaining of this paper is organized as follows. In Section 2, we review some basics on persistence diagrams and kernel embedding methods. In Section 3, the PWGK is proposed, and some theoretical and computational issues are discussed. Section 4 shows experimental results, and compares the proposed kernel method with other methods.

2 Background

We review the concepts of persistence diagrams and kernel methods. For readers who are not familiar with algebraic topology and homology, we give a brief summary in Supplementary material. See also [Hat01] as an accessible introduction to algebraic topology.

2.1 Persistence diagram

Let $X=\{\bm{x}_{1},\ldots,\bm{x}_{n}\}$ be a finite subset in a metric space $(M,d_{M})$. To analyze topological properties of $X$, let us consider a fattened ball model $X_{r}=\bigcup_{i=1}^{n}B(\bm{x}_{i};r)$ consisting of balls $B(\bm{x}_{i};r)=\{\bm{x}\in M\mid d_{M}(\bm{x}_{i},\bm{x})\leq r\}$ with radius $r$, and use the homology $H_{q}(X_{r})$ to describe the topology of $X_{r}$. Here, for a topological space $S$, its $q$-th homology $H_{q}(S)\ (q=0,1,\ldots)$ is defined as a vector space, and its dimension $\dim H_{q}(S)$ counts the number of connected components $(q=0)$, rings $(q=1)$, cavities $(q=2)$, and so on111Throughout this paper we use a field coefficient for homology.. For the precise definition of homology, see Supplementary material. For example, $X_{0.21}$ in Figure 1 consists of one connected component and two rings, and hence $\dim H_{0}(H_{0.21})=1$ and $\dim H_{1}(X_{0.21})=2$.

Because of $X_{r}\subset X_{s}$ for $r\leq s$, the set $\mathbb{X}=\{X_{r}\mid r\geq 0\}$ becomes a filtration222A filtration is a family of subsets $\{X_{a}\mid a\in A\}$ indexed by a totally ordered set $A$ such that $X_{a}\subset X_{b}$ for $a\leq b$.. When the radius changes as in Figure 1, a new generator $\alpha_{i}\in H_{q}(X_{r})$ appears at some radius $r=b_{i}$ and disappears at a radius $r=d_{i}$ larger than $b_{i}$ (called birth and death, respectively). By gathering all generators $\alpha_{i}\ (i\in I)$ in the filtration $\mathbb{X}$, we obtain the collection of these birth-death pairs $\underline{D}_{q}(X)=\{(b_{i},d_{i})\in\mathbb{R}^{2}\mid i\in I\}$ as a multi-set333A multi-set is a set with multiplicity of each point. We regard a persistence diagram as a multi-set, since several generators can have the same birth-death pairs.. The persistence diagram $D_{q}(X)$ is defined by the disjoint union of $\underline{D}_{q}(X)$ and the diagonal set ${\Delta}=\{(a,a)\mid a\in\mathbb{R}\}$ counted with infinite multiplicity. A point $x=(b,d)\in D_{q}(X)$ is also called a generator of the persistence diagram. The persistence ${\rm pers}(x):=d-b$ of $x$ is its lifetime and measures the robustness of $x$ in the filtration. We will see shortly that the diagonal set is included in a persistence diagram to simplify the definition of a distance on persistence diagrams.

Figure 10 shows the persistence diagram $D_{1}(X)$ of $X$ given in Figure 1. The generators $x_{1}$ and $x_{2}$ correspond to the rings $\alpha_{1}$ and $\alpha_{2}$ in Figure 1, respectively. The persistence of $x_{2}$ is the longest, while the other generators including $x_{1}$ have small persistences, implying that they can be seen as noisy rings. Although there are no topological rings in $X$ itself, the persistence diagram $D_{1}(X)$ shows that there is a robust ring $\alpha_{2}$ and several noisy rings in $\mathbb{X}$. In this way, the persistence diagram provides an informative topological summary of $X_{r}$ over all $r$.

We remark that, in the finite fattened ball model, there is only one generator in $D_{0}(X)$ which does not disappear in the filtration; its lifetime is $\infty$. Thus, from now on, we deal with $D_{0}(X)$ by removing this infinite lifetime generator in order to simplify the notation444This is called the reduced persistence diagram.. We also note that the cardinality of $\underline{D}_{q}(X)$ obtained from the finite fattened ball model is finite.

2.2 Stability with respect to $d_{B}$

Any statistical data involve noise or stochasticity, and thus it is desired that the persistence diagrams are stable under perturbation of data. A popular measure to study the similarity between two persistence diagrams $D$ and $E$ is the bottleneck distance

[TABLE]

where $\gamma$ ranges over all multi-bijections555A multi-bijection is a bijective map between two multi-sets counted with their multiplicity. from $D$ to $E$666For $z=(z_{1},z_{2})\in\mathbb{R}^{2}$, $\|z\|_{\infty}$ denotes $\max(|z_{1}|,|z_{2}|)$. . Note that the cardinalities of $D$ and $E$ are equal by considering the diagonal set $\Delta$ with infinite multiplicity. As a distance between finite sets $X,Y$ in a metric space $M$, let us recall the Hausdorff distance given by

[TABLE]

Then, we have the following stability property (for more general settings, see [CdSO14]).

Proposition 2.1**.**

[CdSO14, CSEH07*]**

Let $X$ and $Y$ be finite subsets in a metric space $(M,d_{M})$. Then the persistence diagrams satisfy*

[TABLE]

Proposition 2.1 provides a geometric intuition of the stability of persistence diagrams. Assume that $X$ is the true location of points and $Y$ is a data obtained from skewed measurement with ${\varepsilon}=d_{H}(X,Y)$ (Figure 4). If there is a point $(b,d)\in D_{q}(Y)$, then we can find at least one generator in $X$ which is born in $(b-{\varepsilon},b+{\varepsilon})$ and dies in $(d-{\varepsilon},d+{\varepsilon})$. Thus, the stability guarantees the similarity of two persistence diagrams, and hence we can infer the true topological features from one persistence diagram.

2.3 Kernel methods for representing signed measures

Let $\Omega$ be a set and $k:\Omega\times\Omega{\rightarrow}\mathbb{R}$ be a positive definite kernel on $\Omega$, i.e., $k$ is symmetric, and for any number of points $x_{1},\ldots,x_{n}$ in $\Omega$, the Gram matrix $\left(k(x_{i},x_{j})\right)_{i,j=1,\ldots,n}$ is nonnegative definite. A popular example of positive definite kernel on $\mathbb{R}^{d}$ is the Gaussian kernel $k_{G}(x,y)=e^{-\frac{\left\lVert x-y\right\rVert^{2}}{2\sigma^{2}}}\ (\sigma>0)$, where $\left\lVert\cdot\right\rVert$ is the Euclidean norm in $\mathbb{R}^{d}$. It is also known that every positive definite kernel $k$ on $\Omega$ is uniquely associated with a reproducing kernel Hilbert space $\mathcal{H}_{k}$ (RKHS).

We use a positive definite kernel to represent persistence diagrams by following the idea of the kernel mean embedding of distributions [SGSS07, SFL11]. Let $\Omega$ be a locally compact Hausdorff space, $M_{b}(\Omega)$ be the space of all finite signed Radon measures on $\Omega$, and $k$ be a bounded measurable kernel on $\Omega$. Then we define a mapping from $M_{b}(\Omega)$ to $\mathcal{H}_{k}$ by

[TABLE]

The integral should be understood as the Bochner integral [DU77], which exists here, since $\int\left\lVert k(\cdot,x)\right\rVert_{\mathcal{H}_{k}}d\mu(x)$ is finite.

For a locally compact Hausdorff space $\Omega$, let $C_{0}(\Omega)$ denote the space of continuous functions vanishing at infinity777A function $f$ is said to vanish at infinity if for any ${\varepsilon}>0$ there is a compact set $K\subset\Omega$ such that $\sup_{x\in K^{c}}|f(x)|\leq{\varepsilon}$.. A kernel $k$ on $\Omega$ is said to be $C_{0}$-kernel if $k(x,x)$ is of $C_{0}(\Omega)$ as a function of $x$. If $k$ is $C_{0}$-kernel, the associated RKHS $\mathcal{H}_{k}$ is a subspace of $C_{0}(\Omega)$. A $C_{0}$-kernel $k$ is called $C_{0}$-universal if $\mathcal{H}_{k}$ is dense in $C_{0}(\Omega)$. It is known that the Gaussian kernel $k_{G}$ is $C_{0}$-universal on $\mathbb{R}^{d}$ [SFL11]. When $k$ is $C_{0}$-universal, by the mapping (1), the vector $E_{k}(\mu)$ in the RKHS uniquely determines the finite signed measure $\mu$, and thus serves as a representation of $\mu$.

Proposition 2.2** ([SFL11]).**

If $k$ is $C_{0}$-universal, the mapping $E_{k}$ is injective. Thus,

[TABLE]

defines a distance on $M_{b}(\Omega)$.

3 Kernel methods for persistence diagrams

We propose a kernel for persistence diagrams, called the Persistence Weighted Gaussian Kernel (PWGK), to embed the diagrams into an RKHS. This vectorization of persistence diagrams enables us to apply any kernel methods to persistence diagrams. We show the stability theorem with respect to the distance defined by the embedding, and discuss efficient computation of the PWGK.

3.1 Persistence weighted Gaussian kernel

We propose a method for vectorizing persistence diagrams using the kernel embedding (1) by regarding a persistence diagram as a discrete measure. In vectorizing persistence diagrams, it is important to discount the effect of generators located near the diagonal, since they tend to be caused by noise. To this end, we explain slightly different two ways of embeddings, which turn out to introduce the same inner products for two persistence diagrams.

First, for a persistence diagram $D$, we introduce a weighted measure $\mu^{w}_{D}:=\sum_{x\in\underline{D}}w(x){\delta}_{x}$ with a weight $w(x)>0$ for each generator $x\in\underline{D}$ (Figure 5), where ${\delta}_{x}$ is the Dirac delta measure at $x$. The weight function $w(x)$ discounts the effect of generators close to the diagonal, and a concrete choice will be discussed later.

As discussed in Section 2.3, given a $C_{0}$-universal kernel $k$ on $\mathbb{R}^{2}_{ul}:=\{(b,d)\in\mathbb{R}^{2}\mid b<d\}$, the measure $\mu^{w}_{D}$ can be embedded as an element of the RKHS $\mathcal{H}_{k}$ via

[TABLE]

From Proposition 2.2, this mapping does not lose any information about persistence diagrams, and $E_{k}(\mu^{w}_{D})\in\mathcal{H}_{k}$ serves as a representation of the persistence diagram.

As the second construction, let

[TABLE]

be the weighted kernel with the same weight function as above, and consider the mapping

[TABLE]

This also defines vectorization of persistence diagrams, and it is essentially equivalent to the first one, as seen from the next proposition (See Supplementary material for the proof.).

Proposition 3.1**.**

The following mapping

[TABLE]

defines an isomorphism between the RKHSs. Under this isomorphism, $E_{k}(\mu_{D}^{w})$ and $E_{k^{w}}(\mu_{D})$ are identified.

Note that under the identification of Proposition 3.1, we have

[TABLE]

and thus the two constructions introduce the same similarity (and hence distance) among persistence diagrams. We apply methods of data analysis to vector representations $E_{k}(\mu_{D}^{w})$ or $E_{k^{w}}(\mu_{D})$. The first construction may be more intuitive by the direct weighting of a measure, while the second one is also practically useful since all the parameter tuning is reduced to kernel choice.

For a practical purpose, we propose to use the Gaussian kernel $k_{G}(x,y)=e^{-\frac{\left\lVert x-y\right\rVert^{2}}{2\sigma^{2}}}\ (\sigma>0)$ for $k$, and $w_{{\rm arc}}(x)=\arctan(C{\rm pers}(x)^{p})\ (C,p>0)$ for a weight function. The corresponding positive definite kernel is

[TABLE]

We call it Persistence Weighted Gaussian Kernel (PWGK). Since the Gaussian kernel is $C_{0}$-universal and $w>0$ on $\mathbb{R}^{2}_{ul}$, $d^{w_{{\rm arc}}}_{k_{G}}(D,E):=\|E_{k_{G}}(\mu_{D}^{w_{{\rm arc}}})-E_{k_{G}}(\mu_{E}^{w_{{\rm arc}}})\|_{\mathcal{H}_{k_{G}}}$ defines a distance on the persistence diagrams. We also note that $w_{{\rm arc}}$ is an increasing function with respect to persistence. Hence, a noisy (resp. essential) generator $x$ gives a small (resp. large) value $w_{{\rm arc}}(x)$. By adjusting the parameters $C$ and $p$, we can control the effect of the persistence.

3.2 Stability with respect to $d^{w_{\rm arc}}_{k_{G}}$

Given a data $X$, we vectorize the persistence diagram $D_{q}(X)$ as an element $E_{k_{G}}(\mu_{D_{q}(X)}^{w_{{\rm arc}}})$ of the RKHS. Then, for practical applications, this map $X\mapsto E_{k_{G}}(\mu_{D_{q}(X)}^{w_{{\rm arc}}})$ should be stable with respect to perturbations to the data as discussed in Section 2.2. The following theorem shows that the map has the desired property (See Supplementary material for the proof.).

Theorem 3.2**.**

Let $M$ be a compact subset in $\mathbb{R}^{d}$, $X,Y\subset M$ be finite subsets and $p>d+1$. Then

[TABLE]

where $L(M,d;C,p,\sigma)$ is a constant depending on $M,d,C,p,\sigma$.

Let $\mathcal{P}_{{\rm finite}}(M)$ be the set of finite subsets in a compact subset $M\subset\mathbb{R}^{d}$. Since the constant $L(M,d;C,p,\sigma)$ is independent of $X$ and $Y$, Theorem 3.2 concludes that the map

[TABLE]

is Lipschitz continuous. To the best of our knowledge, a similar stability result has not been obtained for the other Gaussian type kernels (e.g., [RHBK15] does not deal with the Hausdorff distance.). Our stability result is achieved by incorporating the weight function $w_{\rm arc}$ with appropriate choice of $p$.

3.3 Kernel methods on RKHS

Once persistence diagrams are represented by the vectors in an RKHS, we can apply any kernel methods to those vectors. The simplest choice is to consider the linear kernel

[TABLE]

on the RKHS. We can also consider a nonlinear kernel on the RKHS, such as the Gaussian kernel:

[TABLE]

where $\tau$ is a positive parameter and

[TABLE]

Note that we can observe better performance with nonlinear kernels for some complex tasks [MFDS12]. In this paper, we mainly apply the RKHS Gaussian kernel $K_{G}$ in the experimental section. In Section 4, we apply SVM, kernel PCA, and kernel change point detection.

3.4 Computation of Gram matrix

Let $\mathcal{D}=\{D_{\ell}\mid\ell=1,\ldots,n\}$ be a collection of persistence diagrams. In many practical applications, the number of generators in a persistence diagram can be large, while $n$ is often relatively small: in Section 4.3, for example, the number of generators is 30000, while $n=80$.

If the persistence diagrams contain at most $m$ points, each element of the Gram matrix $(K_{G}(D_{i},D_{j}))_{i,j=1,\ldots,n}$ involves $O(m^{2})$ evaluation of $e^{-\frac{\left\lVert x-y\right\rVert^{2}}{2\sigma^{2}}}$, resulting the complexity $O(m^{2}n^{2})$ for obtaining the Gram matrix. Hence, reducing computational cost with respect to $m$ is an important issue, since in many applications $n$ is relatively small

We solve this computational issue by using the random Fourier features [RR07]. To be more precise, let $z_{1},\ldots,z_{M}$ be random variables from the $2$-dimensional normal distribution $N((0,0),\sigma^{-2}I)$ where $I$ is the identity matrix. This method approximates $e^{-\frac{\left\lVert x-y\right\rVert^{2}}{2\sigma^{2}}}$ by $\frac{1}{M}\sum_{a=1}^{M}e^{\sqrt{-1}z_{a}x}(e^{\sqrt{-1}z_{a}y})^{*}$, where $*$ denotes the complex conjugation. Then, $\sum_{x\in\underline{D}_{i}}\sum_{y\in\underline{D}_{j}}w(x)w(y)k_{G}(x,y)$ is approximated by $\frac{1}{M}\sum_{a=1}^{M}B^{a}_{i}(B^{a}_{j})^{*}$, where $B^{a}_{\ell}=\sum_{x\in\underline{D}_{\ell}}w(x)e^{\sqrt{-1}z_{a}x}$. As a result, the computational complexity of the approximated Gram matrix is $O(mnM+n^{2}M)$.

We note that approximation by the random Fourier features can be sensitive to the choice of $\sigma$. If $\sigma$ is much smaller than $\left\lVert x-y\right\rVert$, the relative error can be large. For example, in the case of $x=(1,2),y=(1,2.1)$ and $\sigma=0.01$, $e^{-\frac{\left\lVert x-y\right\rVert^{2}}{2\sigma^{2}}}$ is about $10^{-22}$ while we observed the approximated value can be about $10^{-3}$ with $M=10^{3}$. As a whole, these $m^{2}$ errors may cause a critical error to the approximation. Moreover, if $\sigma$ is largely deviated from the ensemble $\left\lVert x-y\right\rVert$ for $x\in\underline{D}_{i},y\in\underline{D}_{j}$, then most values $e^{-\frac{\left\lVert x-y\right\rVert^{2}}{2\sigma^{2}}}$ become close to [math] or $1$.

In order to obtain a good approximation and extract meaningful values, choice of parameter is important. For supervised learning such as SVM, we use the cross-validation (CV) approach. For unsupervised case, we follow the heuristics proposed in [GFT*+*07]. In Section 4.3, we set $\sigma={\rm median}\{\sigma(D_{\ell})\mid\ell=1,\ldots,n\}$, where $\sigma(D)={\rm median}\{\left\lVert x_{i}-x_{j}\right\rVert\mid x_{i},x_{j}\in\underline{D},\ i<j\}$, so that $\sigma$ takes close values to many $\left\lVert x-y\right\rVert$. For the parameter $C$, we also set $C=({\rm median}\{{\rm pers}(D_{\ell})\mid\ell=1,\ldots,n\})^{-p}$, where ${\rm pers}(D)={\rm median}\{{\rm pers}(x_{i})\mid x_{i}\in\underline{D}\}$. Similarly, $\tau$ is defined by ${\rm median}\{d^{w_{{\rm arc}}}_{k_{G}}(D_{i},D_{j})\mid 1\leq i<j\leq n\}$. In this paper, since all points of data are in $\mathbb{R}^{3}$, we set $p=5$ from the assumption $p>d+1$ in Theorem 3.2.

4 Experiments

We demonstrate the performance of the PWGK using synthesized and real data. In this section, all persistence diagrams are $1$-dimensional (i.e., rings) and computed by CGAL [DLY15] and PHAT [BKRW14].

4.1 Comparison to the persistence scale space kernel

The most relevant work to our method is [RHBK15]. They propose a positive definite kernel called persistence scale space kernel (PSSK for short) $K_{PSS}$ on the persistence diagrams:

[TABLE]

where $\Phi_{t}(D)(x)=\frac{1}{4\pi t}\sum_{y\in\underline{D}}e^{-\frac{\left\lVert x-y\right\rVert^{2}}{4t}}-e^{-\frac{\left\lVert x-\bar{y}\right\rVert^{2}}{4t}}$, and $\bar{y}=(y^{2},y^{1})$ for $y=(y^{1},y^{2})$. Note that the PSSK also takes zero on the diagonal by subtracting the Gaussian kernels for $y$ and $\bar{y}$.

While both methods discount noisy generators, the PWGK has the following advantages over the PSSK. (i) The PWGK can control the effect of the persistence by $C$ and $p$ in $w_{{\rm arc}}$ independently of the bandwidth parameter $\sigma$ in the Gaussian factor, while in the PSSK only one parameter $t$ must control the global bandwidth and the discounting effect. (ii) The approximation by the random Fourier features is not applicable to the PSSK, since it is not shift-invariant in total. We also note that, in [RHBK15], only the linear kernel is considered on the RKHS, while our approach involves a nonlinear kernel on the RKHS.

Regarding the approximation of the PSSK, Nyström method [WS01] or incomplete Cholesky factorization [FS01] can be applied. In evaluating the kernels, we need to calculate $(k(x,y))_{x\in\underline{D}_{i},y\in\underline{D}_{j}}$, which is not symmetric. We then need to apply Nyström or incomplete Cholesky to the symmetric but big positive definite matrix $(k(x,y))_{\begin{subarray}{c}x\in\underline{D}_{i},y\in\underline{D}_{j}\\ i,j=1,\dots,n\end{subarray}}$ of size $O(nm)$. Either, we need to apply incomplete Cholesky to the non-symmtric matrix $(k(x,y))_{x\in\underline{D}_{i},y\in\underline{D}_{j}}$ for all the combination of $(i,j)$, which requires considerable computational cost for large $n$. In contrast, the random Fourier features can be applied to the kernel function irrespective to evaluation points, and the same Fourier expansion can be applied to any $(i,j)$. This guarantees the efficient computational cost.

The detailed comparisons will be experimentally verified in Sections 4.2 and 4.3. With respect to the parameter $t$ in the PSSK, since $e^{-\frac{\left\lVert x-\bar{y}\right\rVert^{2}}{8t}}$ is only used to vanish the value of the feature map $\Phi_{t}(D)$ on the diagonal, we set $t=\frac{\sigma^{2}}{4}$ by using the same $\sigma$ defined in Section 3.4.

4.2 Classification with synthesized data

We first use the proposed method for a classification task with SVM, and compare the performance with the PSSK. The synthesized data are generated as follows. Each data set assumes one or two circles, and data points are located at even spaces along the circle(s). It always contains one larger circle $S_{1}$ of radius $r_{1}^{o}$ raging from 1 to 10, and it may have a smaller circle $S_{2}$ of radius 0.2 (10 points) with probability 1/2. Roughly speaking, the class label $Y$ is made by ${\bf XOR}(z_{0},z_{1})$, where $z_{i}$ ($i=0,1)$ is a binary variable: $z_{0}=1$ if the smaller $S_{2}$ exists, and $z_{1}=1$ if the birth and death $(b^{o},d^{o})$ of the generator corresponding to $S_{1}$ satisfies $b^{o}<A_{B}$ and $d^{o}>A_{D}$ for fixed thresholds $A_{B},A_{D}$. We can control $b^{o}$ and $d^{o}$ by choosing the number of points $N_{1}^{o}$ along $S_{1}$ (for birth) and the radius $r_{1}^{o}$ (for death). To generate the data points, we add noise effects to make the classification harder: the radius $r_{1}$ and sample size $N_{1}$ are in fact given by adding noise to $r_{1}^{0}$ and $N_{1}^{o}$, and $S_{1}$ are made according to the shifted $r_{1}$ and $N_{1}$, while the class label is given by the non-shifted $r_{1}^{o}$ and $N_{1}^{o}$. For the precise description of the data generation procedure, see Supplementary material. By this construction, the classifier needs to look at both of the location of the generator and the existence of the generator for the smaller one around the diagonal.

SVMs are trained with persistence diagrams given by 100 data sets, and evaluated with 99 independent test data sets. For the kernel on RKHS, we used both of the linear and Gaussian kernels. The hyper-parameters $(\sigma,C)$ in the PWGK and $t$ in the PSSK are chosen by the 10-fold cross-validation, and the degree $p$ in the weight of the PWGK is set to be $5$. The variance parameter in the RKHS-Gaussian kernel is set by the median heuristics. We also apply the Gaussian kernel (without any weights) for embedding persistence diagrams to RKHS.

In Table 1, we can see that the PSSK does not work well for this problem, even worse than the Gaussian kernel, and the classification rate by the linear RKHS kernel used originally in [RHBK15] is almost the chance level. This must be caused by the difficulty in handling the global location of generators and close look around the diagonal simultaneously. This classification task involves strong nonlinearity on the RKHS, as seen in the large improvement by PWGK+Gauss kernel.

4.3 Analysis of ${\rm{\bf SiO_{2}}}$

In this experiment, we compare the PWGK and the PSSK to the non-trivial problem of glass transition on ${\rm SiO}_{2}$, focusing also on their computational efficiency.

When we rapidly cool down the liquid state of ${\rm SiO}_{2}$, it avoids the usual crystallization and changes into a glass state. Understanding the liquid-glass transition is an important issue for the current physics and industrial applications [GS07]. For estimating the glass transition temperature by simulations, we first prepare atomic configurations of ${\rm SiO_{2}}$ for a certain range of temperatures, and then draw the temperature-enthalpy graph. The graph consists of two lines in high and low temperatures with slightly different slopes which correspond to the liquid and the glass states, respectively, and the glass transition temperature is conventionally estimated as an interval of the transient region combining these two lines (e.g., see [Ell90]). However, since the slopes of the two lines are close to each other, determining the interval is a subtle problem, and usually the rough estimate of the interval is only available. Hence, it is desired to develop a mathematical framework to detect the glass transition temperature.

Our strategy is to regard the glass transition temperature as the change point and detect it from a collection $\mathcal{D}=\{D_{\ell}\mid\ell=1,\ldots,n\}$ of persistence diagrams made by atomic configurations of ${\rm SiO}_{2}$, where $\ell$ is the index of the temperatures listed in the decreasing order. We use the kernel Fisher discriminant ratio ${\rm KFDR}_{n,\ell,\gamma}(\mathcal{D})$ [HMB09] as a statistical quantity for the change point detection. Here, we set $\gamma=10^{-3}$ in this paper, and the index $\ell$ achieving the maximum of ${\rm KFDR}_{n,\ell,\gamma}(\mathcal{D})$ corresponds to the estimated change point. The ${\rm KFDR}$ is calculated by the Gram matrix $(K(D_{i},D_{j}))_{i,j=1,\ldots,n}$ with respect to the kernel $K$.

We compute $\mathcal{D}$ with $n=80$ from the data used in [NHH*+*15a, NHH*+*15b]. Since the persistence diagrams of ${\rm SiO_{2}}$ contain huge amount of points, we apply the random Fourier features and the Nyström methods [DM05] for the approximations of the PWGK with the Gaussian RKHS and the PSSK, respectively. The sample sizes used in both approximations are denoted by $M$ and $c$, where $c$ is the number of chosen columns. Figrue 6 summarizes the plots of the change points for several sample sizes and the computational time.

The interval of the glass transition temperature $T$ estimated by the conventional method explained above is $2000K\leq T\leq 3500K$, which corresponds to $35\leq\ell\leq 50$.

The computational complexity of the random Fourier features with respect to the sample size is $O(M)$, while that of the Nyström method involves matrix inversion of $O(c^{3})$. For this reason, the PSSK with $c>500$ cannot be performed in reasonable time, and hence we cannot check the convergence of the change points with respect to the sample size as shown in Figure 6. On the other hand, the PWGK plot shows the convergence to $\ell=37$, implying that $\ell=37$ is the true change point. We here emphasize that the computation to obtain $\ell=37$ by the PWGK is much faster than the PSSK.

Figure 7 shows the normalized plots of ${\rm KFDR}_{n,\ell,\gamma}(\mathcal{D})$ and the $2$-dimensional plot given by KPCA (the color is given by the result of the change point detection by the PWGK). As we see from the figure, the KPCA plot shows the clear phase change between before (red) and after (blue) the change point. This strongly suggests that the glass transition occurs at the detected change point.

4.4 Protein classification

We apply the PWGK to two classification tasks studied in [CMW*+*15]. They use the molecular topological fingerprint (MTF) as a feature vector for the input to the SVM. The MTF is given by the 13 dimensional vector whose elements consist of the persistences of some specific generators (e.g., the longest, second longest, etc.) in persistence diagrams. We compare the performance of the PWGK with the Gaussian RKHS kernel and the MTF method under the same setting of the SVM reported in [CMW*+*15].

The first task is a protein-drug binding problem, and we classify the binding and non-binding of drug to the M2 channel protein of the influenza A virus. For each form, 15 data were obtained by NMR experiments, in which 10 data are used for training and the remaining for testing. We randomly generated 100 ways of partitions, and calculated the classification rates.

In the second problem, the taut and relaxed forms of hemoglobin are to be classified. For each form, 9 data were collected by the X-ray crystallography. We select one data from each class for testing, and use the remaining for training. All the 81 combinations are performed to calculate the CV classification rates.

The results of the two problems are shown in Table 2. We can see that the PWGK achieves better performance than the MTF in both problems.

5 Conclusion

In this paper, we have proposed a kernel framework for analysis with persistence diagrams, and the persistence weighted Gaussian kernel as a useful kernel for the framework. As a significant advantage, our kernel enables one to control the effect of persistence in data analysis. We have also proven the stability result with respect to the kernel distance. Furthermore, we have analyzed the synthesized and real data by using the proposed kernel. The change point detection, the principal component analysis, and the support vector machine using the PWGK derived meaningful results in physics and biochemistry. From the viewpoint of computations, our kernel provides an accurate and efficient approximation to compute the Gram matrix, suitable for practical applications of TDA.

Supplementary Material

This supplementary material provides a brief introduction of topological tools used in the paper, the proof of Theorem 3.1, the proof of Proposition 3.1, and explanations of synthesized data used in Section 4.2. In oder to prove Theorem 3.1, we introduce sub-level sets in Section B and total persistence in Section C, and we will prove a generalization of Theorem 3.1 in Section D.

Appendix A Topological tools

This section summarizes some topological tools used in the paper. In general, topology concerns geometric properties invariant to continuous deformations. To study topological properties algebraically, simplicial complexes are often considered as basic objects. We start with a brief explanation of simplicial complexes, and gradually increase the generality from simplicial homology to singular and persistent homology. For more details, see [Hat01].

A.1 Simplicial complex

We first introduce a combinatorial geometric model called simplicial complex to define homology. Let $P=\{1,\dots,n\}$ be a finite set (not necessarily points in a metric space). A simplicial complex with the vertex set $P$ is defined by a collection $S$ of subsets in $P$ satisfying the following properties:

$\{i\}\in S$ for $i=1,\dots,n$, and 2. 2.

if $\sigma\in S$ and $\tau\subset\sigma$, then $\tau\in S$.

Each subset $\sigma$ with $q+1$ vertices is called a $q$-simplex. We denote the set of $q$-simplices by $S_{q}$. A subcollection $T\subset S$ which also becomes a simplicial complex (with possibly less vertices) is called a subcomplex of $S$.

We can visually deal with a simplicial complex $S$ as a polyhedron by pasting simplices in $S$ into a Euclidean space. The simplicial complex obtained in this way is called a geometric realization, and its polyhedron is denoted by $|S|$. In this context, the simplices with small $q$ correspond to points ($q=0$), edges ($q=1$), triangles ($q=2$), and tetrahedra ($q=3$).

Example A.1**.**

Figure 8 shows two polyhedra of simplicial complexes

[TABLE]

A.2 Homology

A.2.1 Simplicial homology

The procedure to define homology is summarized as follows:

Given a simplicial complex $S$, build a chain complex $C_{*}(S)$. This is an algebraization of $S$ characterizing the boundary. 2. 2.

Define homology by quotienting out certain subspaces in $C_{*}(S)$ characterized by the boundary.

We begin with the procedure 1 by assigning orderings on simplices. When we deal with a $q$-simplex $\sigma=\{i_{0},\dots,i_{q}\}$ as an ordered set, there are $(q+1)!$ orderings on $\sigma$. For $q>0$, we define an equivalence relation $i_{j_{0}},\dots,i_{j_{q}}\sim i_{\ell_{0}},\dots,i_{\ell_{q}}$ on two orderings of $\sigma$ such that they are mapped to each other by even permutations. By definition, two equivalence classes exist, and each of them is called an oriented simplex. An oriented simplex is denoted by $\langle i_{j_{0}},\dots,i_{j_{q}}\rangle$, and its opposite orientation is expressed by adding the minus $-\langle i_{j_{0}},\dots,i_{j_{q}}\rangle$. We write $\langle\sigma\rangle=\langle i_{j_{0}},\dots,i_{j_{q}}\rangle$ for the equivalence class including $i_{j_{0}}<\dots<i_{j_{q}}$. For $q=0$, we suppose that we have only one orientation for each vertex.

Let $K$ be a field. We construct a $K$-vector space $C_{q}(S)$ as

[TABLE]

for $S_{q}\neq\emptyset$ and $C_{q}(S)=0$ for $S_{q}=\emptyset$. Here, ${\rm Span}_{K}(A)$ for a set $A$ is a vector space over $K$ such that the elements of $A$ formally form a basis of the vector space. Furthermore, we define a linear map called the boundary map $\partial_{q}:C_{q}(S)\rightarrow C_{q-1}(S)$ by the linear extension of

[TABLE]

where $\widehat{i_{\ell}}$ means the removal of the vertex $i_{\ell}$. We can regard the linear map $\partial_{q}$ as algebraically capturing the $(q-1)$-dimensional boundary of a $q$-dimensional object.

For example, the image of the $2$-simplex $\langle\sigma\rangle=\langle 1,2,3\rangle$ is given by $\partial_{2}\langle\sigma\rangle=\langle 2,3\rangle-\langle 1,3\rangle+\langle 1,2\rangle$, which is the boundary of $\sigma$ (see Figure 8).

In practice, by arranging some orderings of the oriented $q$- and $(q-1)$- simplices, we can represent the boundary map as a matrix $M_{q}=(M_{\sigma,\tau})_{\sigma\in S_{q-1},\tau\in S_{q}}$ with the entry $M_{\sigma,\tau}=0,\pm 1$ given by the coefficient in (6). For the simplicial complex $S$ in Example A.1, the matrix representations $M_{1}$ and $M_{2}$ of the boundary maps are given by

[TABLE]

Here the $1$-simplices (resp. [math]-simplices) are ordered by $\langle 1,2\rangle,\langle 2,3\rangle,\langle 1,3\rangle$ (resp. $\langle 1\rangle$, $\langle 2\rangle$, $\langle 3\rangle$).

We call a sequence of the vector spaces and linear maps

[TABLE]

the chain complex $C_{*}(S)$ of $S$. As an easy exercise, we can show $\partial_{q}\circ\partial_{q+1}=0$. Hence, the subspaces $Z_{q}(S)={\rm ker}\partial_{q}$ and $B_{q}(S)={\rm im}\partial_{q+1}$ satisfy $B_{q}(S)\subset Z_{q}(S)$. Then, the $q$-th (simplicial) homology is defined by taking the quotient space

[TABLE]

Intuitively, the dimension of $H_{q}(S)$ counts the number of $q$-dimensional holes in $S$ and each generator of the vector space $H_{q}(S)$ corresponds to these holes. We remark that the homology as a vector space is independent of the orientations of simplices.

For a subcomplex $T$ of $S$, the inclusion map $\rho:T\hookrightarrow S$ naturally induces a linear map in homology $\rho_{q}:H_{q}(T)\rightarrow H_{q}(S)$. Namely, an element $[c]\in H_{q}(T)$ is mapped to $[c]\in H_{q}(S)$, where the equivalence class $[c]$ is taken in each vector space.

For example, the simplicial complex $S$ in Example A.1 has

[TABLE]

from (13). Hence $H_{1}(S)=0$, meaning that there are no $1$-dimensional hole (ring) in $S$. On the other hand, since $Z_{1}(T)=Z_{1}(S)$ and $B_{1}(T)=0$, we have $H_{1}(T)\simeq K$, meaning that $T$ consists of one ring. Hence, the induced linear map $\rho_{1}:H_{1}(T)\rightarrow H_{1}(S)$ means that the ring in $T$ disappears in $S$ under $T\hookrightarrow S$.

A topological space $X$ is called triangulable if there exists a geometric realization of a simplicial complex $S$ whose polyhedron is homeomorphic888A continuous map $f:X{\rightarrow}Y$ is said to be homeomorphic if $f:X{\rightarrow}Y$ is bijective and the inverse $f^{-1}:Y{\rightarrow}X$ is also continuous. to $X$. For such a triangulable topological space, the homology is defined by $H_{q}(X)=H_{q}(S)$. This is well-defined, since a different geometric realization provides an isomorphic homology.

A.2.2 Singular homology

We here extend the homology to general topological spaces. Let $e_{0},\dots,e_{q}$ be the standard basis of ${\mathbb{R}}^{q+1}$ (i.e., $e_{i}=(0,\dots,0,1,0,\dots,0)$, 1 at $(i+1)$-th position, and 0 otherwise), and set

[TABLE]

We also denote the inclusion by $\iota^{\ell}_{q}:\Delta^{\ell}_{q}\hookrightarrow\Delta_{q}$.

For a topological space $X$, a continuous map $\sigma:\Delta_{q}\rightarrow X$ is called a singular $q$-simplex, and let $X_{q}$ be the set of $q$-simplices. We construct a $K$-vector space $C_{q}(X)$ as

[TABLE]

The boundary map $\partial_{q}:C_{q}(X)\rightarrow C_{q-1}(X)$ is defined by the linear extension of

[TABLE]

Even in this setting, we can show that $\partial_{q}\circ\partial_{q+1}=0$, and hence the subspaces $Z_{q}(X)={\rm ker}\partial_{q}$ and $B_{q}(X)={\rm im}\partial_{q+1}$ satisfy $B_{q}(X)\subset Z_{q}(X)$. Then, the $q$-th (singular) homology is similarly defined by

[TABLE]

It is known that, for a triangulable topological space, the homology of this definition is isomorphic to that defined in A.2.1. From this reason, we hereafter identify simplicial and singular homology.

The induced linear map in homology for an inclusion pair of topological space $Y\subset X$ is similarly defined as in A.2.1.

A.3 Persistent homology

Let $\mathbb{X}=\{X_{a}\mid a\in\mathbb{R}\}$ be a (right continuous) filtration of topological spaces, i.e., $X_{a}\subset X_{b}$ for $a\leq b$ and $X_{a}=\bigcap_{a<b}X_{b}$. For $a\leq b$, we denote the linear map induced from $X_{a}\hookrightarrow X_{b}$ by $\rho^{b}_{a}:H_{q}(X_{a})\rightarrow H_{q}(X_{b})$. The persistent homology $H_{q}(\mathbb{X})=(H_{q}(X_{a}),\rho^{b}_{a})$ of $\mathbb{X}$ is defined by the family of homology $\{H_{q}(X_{a})\mid a\in\mathbb{R}\}$ and the induced linear maps $\rho^{b}_{a}$ for all $a\leq b$.

A homological critical value of $H_{q}(\mathbb{X})$ is the number $a\in\mathbb{R}$ such that the linear map $\rho^{a+{\varepsilon}}_{a-{\varepsilon}}:H_{q}(X_{a-{\varepsilon}})\rightarrow H_{q}(X_{a+{\varepsilon}})$ is not isomorphic for any ${\varepsilon}>0$. The persistent homology $H_{q}(\mathbb{X})$ is called tame, if $\dim H_{q}(X_{a})<\infty$ for any $a\in\mathbb{R}$ and the number of homological critical values is finite. A tame persistent homology $H_{q}(\mathbb{X})$ has a nice decomposition property:

Theorem A.2** ([ZC05]).**

A tame persistent homology can be uniquely expressed by

[TABLE]

where $I[b_{i},d_{i}]=(U_{a},\iota^{b}_{a})$ consists of a family of vector spaces

[TABLE]

and the identity map $\iota^{b}_{a}={\rm id}_{K}$ for $b_{i}\leq a\leq b<d_{i}$.

Each summand $I[b_{i},d_{i}]$ is called a generator of the persistent homology and $(b_{i},d_{i})$ is called its birth-death pair. We note that, when $\dim H_{q}(X_{a})\neq 0$ for any $a<0$ (or for any $a>0$ resp.), the decomposition (14) should be understood in the sense that some $b_{i}$ takes the value $-\infty$ (or $d_{i}=\infty$, resp.), where $-\infty,\infty$ are the elements in the extended real $\overline{\mathbb{R}}=\mathbb{R}\cup\{-\infty,\infty\}$.

From the decomposition in Theorem A.2, we define a multiset

[TABLE]

The persistence diagram $D_{q}(\mathbb{X})$ is defined by the disjoint union of $\underline{D}_{q}(\mathbb{X})$ and the diagonal set $\Delta=\{(a,a)\mid a\in\overline{\mathbb{R}}\}$ with infinite multiplicity.

By definition, a generator $(b,d)\in D_{q}(\mathbb{X})$ close to $\Delta$ possesses a short lifetime, implying a noisy topological feature under parameter changes. On the other hand, a generator far away from $\Delta$ can be regarded as a robust feature.

As an example, we show a filtration of simplicial complexes $\mathbb{X}$ in Figure 9 and its persistence diagram $D_{1}(\mathbb{X})$ in Figure 10.

In the paper, we considered a filtration $\mathbb{X}=\{X_{r}\mid r\in\mathbb{R}\}$ generated by a finite subset $X=\{\bm{x}_{1},\ldots,\bm{x}_{n}\}$ in some metric space $(M,d_{M})$ by $X_{r}=\bigcup_{i=1}^{n}B(\bm{x}_{i};r)$, where $B(\bm{x};r)=\{\bm{x}^{\prime}\in M\mid d_{M}(\bm{x},\bm{x}^{\prime})\leq r\}$. For notational simplicity, we denote the persistence diagram of this filtration model by $D_{q}(X)$.

Appendix B Sub-level sets

A popular way of constructing a filtration of topological spaces and thus persistent homology is to use the sub-level sets of a function.

Given a function $f:M{\rightarrow}\mathbb{R}$, a sub-level set $F_{a}=\{\bm{x}\in M\mid f(\bm{x})\leq a\}$ defines a filtration $\mathbb{F}=\{F_{a}\mid a\in\mathbb{R}\}$ (Figure 11) and its persistent homology $H_{q}(\mathbb{F})$. A function $f$ is said to be tame if the persistent homology $H_{q}(\mathbb{F})$ is tame for all $q$. For a tame function $f:M{\rightarrow}\mathbb{R}$, the persistence diagram can be defined and denoted by $D_{q}(f)$. Note that, for a finite set $X=\{\bm{x}_{1},\ldots,\bm{x}_{n}\}$ in $M\subset\mathbb{R}^{d}$ and the function $f_{X}:M{\rightarrow}\mathbb{R}$ defined by $f_{X}(\bm{x})=\min_{\bm{x}_{i}\in X}d_{M}(\bm{x},\bm{x}_{i})$, $f_{X}$ is tame and the persistence diagram $D_{q}(f_{X})$ of the sub-level set is the same as $D_{q}(X)$.

In case of sub-level sets, the stability of the persistence diagram is shown in [CSEH07]. Here, we measure the difference between two functions $f,g:M{\rightarrow}\mathbb{R}$ by $\left\lVert f-g\right\rVert_{\infty}=\sup_{\bm{x}\in M}|f(\bm{x})-g(\bm{x})|$.

Proposition B.1** ([CSEH07]).**

Let $M$ be a triangulable compact metric space with continuous tame functions $f,g:M{\rightarrow}\mathbb{R}$. Then the persistence diagrams satisfy

[TABLE]

For finite subsets $X$ and $Y$ in a metric space $M$, since the difference of norm $\left\lVert f_{X}-f_{Y}\right\rVert_{\infty}$ is nothing but the Hausdorff distance $d_{H}(X,Y)$, we obtain Proposition 2.1 as a corollary.

Appendix C Total persistence

In order to estimate constants $L$ appearing in Theorem 3.1, we will review several properties of persistence.

Let $(M,d_{M})$ be a triangulable compact metric space. For a Lipschitz function $f:M{\rightarrow}\mathbb{R}$, the degree-$p$ total persistence is defined by

[TABLE]

for $0\leq t\leq{\rm Amp}(f)$. Here, ${\rm Amp}(f):=\max_{\bm{x}\in M}f(\bm{x})-\min_{\bm{x}\in M}f(\bm{x})$ is the amplitude of $f$. Let $K$ be a triangulated simplicial complex of $M$ by a homeomorphism $\vartheta:|K|{\rightarrow}M$. The diameter of a simplex $\sigma\in K$ and the mesh of the triangulation $K$ are defined by ${\rm diam}(\sigma)=\max_{\bm{x},\bm{y}\in\sigma}d_{M}(\vartheta(\bm{x}),\vartheta(\bm{y}))$ and ${\rm mesh}(K)=\max_{\sigma\in K}{\rm diam}(\sigma)$, respectively. Furthermore, let us set $N(r)=\min_{{\rm mesh}(K)\leq r}{\rm card}\left(K\right)$. The degree-$p$ total persistence is bounded as follows:

Lemma C.1** ([CSEHM10]).**

Let $M$ be a triangulable compact metric space and $f:M{\rightarrow}\mathbb{R}$ be a tame Lipschitz function. Then ${\rm Pers}_{p}(f,t)$ is bounded from above by

[TABLE]

where ${\rm Lip}(f)$ is the Lipschitz constant of $f$.

For a compact triangulable subspace $M$ in $\mathbb{R}^{d}$, the number of $d$-cubes with length $r>0$ covering $M$ is bounded for every $r>0$. The minimum of these numbers is $N(r)$ and bounded from above by $C_{M}/r^{d}$ for some constant $C_{M}$ depending only on $M$.

For $p>d$, we can find upper bounds for the both terms as follows:

[TABLE]

and

[TABLE]

We note $\lim_{t{\rightarrow}0}t^{p}N(t/{\rm Lip}(f))=0$ for $p>d$. Then, an upper bound of the total persistence ${\rm Pers}(f):={\rm Pers}_{p}(f,0)$ is given as follows:

Lemma C.2**.**

Let $M$ be a triangulable compact subspace in $\mathbb{R}^{d}$ and $p>d$. For any Lipschitz function $f:M{\rightarrow}\mathbb{R}$,

[TABLE]

where $C_{M}$ is a constant depending only on $M$.

In case of a finite subset $X\subset\mathbb{R}^{d}$, there always exists an $R$-ball $M$ containing $X$ for some $R>0$, which is a triangulable compact subspace in $\mathbb{R}^{d}$, and the total persistence of $f_{X}$ is bounded as follows:

Lemma C.3**.**

Let $M$ be a triangulable compact subspace in $\mathbb{R}^{d}$, $X=\{\bm{x}_{1},\ldots,\bm{x}_{n}\}$ be a finite subset of $M$, and $p>d$. Then

[TABLE]

where $C_{M}$ is a constant depending only on $M$.

Proof.

The Lipschitz constant of $f_{X}$ is $1$, because, for any $\bm{x},\bm{y}\in M$,

[TABLE]

Moreover, the amplitude ${\rm Amp}(f_{X})$ is less than or equal to ${\rm diam}(M)$, because $\min_{\bm{x}\in M}f_{X}(\bm{x})=0$ and $\max_{\bm{x}\in M}f_{X}(\bm{x})\leq{\rm diam}(M)$. ∎

For more general, for a persistence diagram $D$, we define $p$-degree total persistence of $D$ by ${\rm Pers}_{p}(D):=\sum_{x\in D}{\rm pers}(x)^{p}$. Let $x_{1},\ldots,x_{n}$ be points of $\underline{D}$. Then we can consider the $n$-dimensional vector

[TABLE]

from $D$. Since each ${\rm pers}(x)\ (x\in\underline{D})$ is always positive, by using the $\ell^{p}$-norm of $v(D)$, we have ${\rm Pers}_{p}(D)=\left\lVert v(D)\right\rVert^{p}_{p}$. In general, since $\left\lVert v\right\rVert_{q}\leq\left\lVert v\right\rVert_{p}\ (v\in\mathbb{R}^{n},\ 1\leq p\leq q)$, we have

[TABLE]

Proposition C.4**.**

If $1\leq p\leq q<\infty$ and ${\rm Pers}_{p}(D)$ is bounded, ${\rm Pers}_{q}(D)$ is also bounded.

Appendix D Proof of Theorem 3.1 and its generalization

We can have the following generalized stability result:

Theorem D.1**.**

Let $D$ and $E$ be finite persistence diagrams. Then

[TABLE]

where $L(D,E;C,p,\sigma)$ is a constant depending on $D,E,C,p,\sigma$.

In fact, we can calculate the constant $L(D,E;C,p,\sigma)$ in Theorem D.1 such as

[TABLE]

Actually, this constant is dependent on $D$ and $E$, and hence we cannot say that the map $D\mapsto E_{k_{G}}(\mu^{w_{\rm arc}}_{D})$ is continuous. However, in the case of persistence diagrams obtained from finite sets, this constant becomes independent of $D$ and $E$. From now on, to obtain the constant $L(D,E;C,p,\sigma)$, we show several lemmas.

Lemma D.2**.**

For any $x,y\in\mathbb{R}^{2}$, $\left\lVert k_{G}(\cdot,x)-k_{G}(\cdot,y)\right\rVert_{\mathcal{H}_{k_{G}}}\leq\frac{\sqrt{2}}{\sigma}\left\lVert x-y\right\rVert_{\infty}$.

Proof.

[TABLE]

We have used the fact $1-e^{-t}\leq t\ (t\in\mathbb{R})$ in (15) and $\left\lVert x\right\rVert^{2}\leq 2\left\lVert x\right\rVert^{2}_{\infty}\ (x\in\mathbb{R}^{2})$ in (16). ∎

Lemma D.3**.**

For any $x,y\in\mathbb{R}^{2}$, the difference of persistences $\left\lvert{\rm pers}(x)-{\rm pers}(y)\right\rvert$ is less than or equal to $2\left\lVert x-y\right\rVert_{\infty}$.

Proof.

For $x=(x_{1},x_{2}),y=(y_{1},y_{2})$, we have

[TABLE]

∎

Lemma D.4**.**

For any $x,y\in\mathbb{R}^{2}$, we have

[TABLE]

Proof.

[TABLE]

We have used the fact that the Lipschitz constant of $\arctan$ is $1$ in (17) and for any $s,t\in\mathbb{R}$,

[TABLE]

in (18). ∎

Proof of Theorem D.1.

Here, let a multi-bijection $\gamma:D{\rightarrow}E$ such that $\sup_{x\in D}\left\lVert x-\gamma(x)\right\rVert_{\infty}\leq{\varepsilon}$ and $E^{\prime}:=\underline{E}-(\gamma(\underline{D})-{\Delta})$. By the definition of $\gamma$, each point in $E^{\prime}$ is ${\varepsilon}$-close to the diagonal. Then, we evaluate $d^{w_{{\rm arc}}}_{k_{G}}(D,E)$ as follows:

[TABLE]

Form Lemma D.2, we have

[TABLE]

from Lemma D.4,

[TABLE]

Moreover, for any $x\in\mathbb{R}^{2}$ we have $\left\lVert k_{G}(\cdot,x)\right\rVert_{\mathcal{H}_{k_{G}}}=\sqrt{k_{G}(x,x)}=1$. Therefore, the continuation of (19) is

[TABLE]

We have used the fact $w_{\rm arc}(x)\leq C{\rm pers}(x)^{p}$ in (20), ${\rm pers}(y)\leq{\varepsilon}\ (y\in E^{\prime})$ in (21) and

[TABLE]

in (23). Thus, if both $(p-1)$-degree total persistence of $D$ and that of $E$ are bounded, since $p$-degree total persistence of $D$ is also bounded from Proposition C.4, the coefficient of ${\varepsilon}$ appearing in (23) is bounded. ∎

Proof of Theorem 3.1.

For any finite set $X\subset M$, from Lemma C.3, there exists a constant $C_{M}>0$ such that

[TABLE]

By replacing $D$ and $E$ with $D_{q}(X)$ and $D_{q}(Y)$ in (23), respectively, we have

[TABLE]

This value is a constant dependent on $M,d,C,p,\sigma$ but independent of $X,Y$. ∎

Appendix E Proof of Proposition 3.1

Proof.

Let

[TABLE]

and define its inner product by

[TABLE]

Using $w>0$, it is easy to see that $\tilde{\mathcal{H}}$ is a Hilbert space, and the mapping $f\mapsto wf$ is an isomorphism of the Hilbert spaces. Moreover, we can see $\tilde{\mathcal{H}}$ is in fact $\mathcal{H}_{k^{w}}$. To see this, it is sufficient to see that $k^{w}(\cdot,x)=w(\cdot)w(x)k(\cdot,x)$ is a reproducing kernel of $\tilde{\mathcal{H}}$. Then, the uniqueness of a reproducing kernel for an RKHS completes the proof. The reproducing property is in fact proved from

[TABLE]

The second assertion is obvious from Equations (2) and (3) in Section 3.1. ∎

Appendix F Synthesized data

We describe the details of synthesized data used in Section 4.2.

Each data set contains a variable number of 2 dimensional points, which are located on one or two circles. A data set always contains points along a larger circle $S_{1}$ of random center $(x_{1},y_{1})$ and radius $r_{1}$ raging from 1 to 10. $N_{1}$ points are located along $S_{1}$ at even spaces. Another circle $S_{2}$ of radius $r_{2}=0.2$ at the origin exists with probability 0.5, and 10 points are located at even spaces, if exists. The variables $r_{1}$ and $N_{1}$ are in fact noisy version of $r_{1}^{o}$ and $N_{1}^{o}$, which are also randomly given. The generative model of these variables are shown later.

Consider the persistence diagram for the $r$-ball model with the above points, and let $(b,d)$ ($b\leq d$) be the birth-death coordinate for the generator corresponding to $S_{1}$. Note that there are no noise to the points on the circle, once the circle is given, and $(b,d)$ is essentially determined by $N_{1}$ and $r_{1}$: if $N_{1}$ is sufficiently large, the birth time $b$ is approximately $\pi r_{1}/N_{1}$ (here we use $\sin\theta\approx\theta$ for a small $\theta$), and the death time $d$ is equal to $r_{1}$, if no intersection occurs with balls around the points on $S_{2}$. See Figure 12.

For the binary classification, the class label is assigned depending on whether or not $S_{2}$ exists, and whether or not the hypothetical (or true) circle and virtual data points given by $r_{1}^{o}$ and $N_{1}^{o}$ has long persistence. Let $A_{B}$ and $A_{D}$ be constants (we set $A_{B}=1$ and $A_{D}=4$ in the experiments). We introduce binary variables $z_{0}$ and $z_{1}$; $z_{0}=1$ if $S_{2}$ exists, and $z_{1}=1$ if $\pi r_{1}^{o}/N_{1}^{o}\leq A_{B}$ and $r_{1}^{o}\geq A_{D}$. Note that $z_{1}=1$ if the hypothetical (or true) generator $(b^{o},d^{o})$ approximately satisfy $b^{o}\leq A_{B}$ and $d^{o}\geq A_{D}$. The class label $Y$ of the data set is then given by

[TABLE]

The generative models of $r_{1},N_{1},r_{1}^{o}$ and $N_{1}^{o}$ are given as follows. The radius $r_{1}$ is generated by

[TABLE]

and the number of points $N_{1}$ is given by

[TABLE]

The center point is given by $(x_{1},y_{1})=1.5(r_{1},r_{1})+(V_{1}^{2},V_{2}^{2})$ with independent $V_{1},V_{2}\sim N(0,2)$. The “true” radius $r_{1}^{o}$ is generated by

[TABLE]

and the number of points $N_{1}^{o}$ is a random integer with equal probability in $[\lceil N_{*}/2\rceil,4N_{*}]$, where $N_{*}:=\pi r_{1}^{o}/A_{B}$. Note that $N_{*}$ points on the circle of radius $r_{1}^{o}$ approximately give a birth time $A_{B}$. If $N_{1}^{o}>N_{*}$ the true birth $b^{o}$ is smaller than $A_{B}$ approximately.

Note that the class label $Y$ is given based on $r_{1}^{o}$ and $N_{1}^{o}$, while the data points and persistence diagram are based on the noisy version $r_{1}$ and $N_{1}$. By this randomness, the best classification boundary does not give 100$\%$ classification rate. Figure 13 shows two examples of data set.