A Bayesian Framework for Persistent Homology

TL;DR

This paper introduces a Bayesian approach to analyze persistence diagrams in topological data analysis, modeling them as Poisson point processes and enabling statistical inference with conjugate priors.

Contribution

It develops a novel Bayesian framework for persistence diagrams using point process theory, providing closed-form posterior intensities with Gaussian mixture priors.

Findings

Effective classification in materials science using Bayes factors

Closed-form posterior intensities derived from Gaussian mixture priors

Framework extends statistical inference capabilities for topological data analysis

Abstract

Persistence diagrams offer a way to summarize topological and geometric properties latent in datasets. While several methods have been developed that utilize persistence diagrams in statistical inference, a full Bayesian treatment remains absent. This paper, relying on the theory of point processes, presents a Bayesian framework for inference with persistence diagrams relying on a substitution likelihood argument. In essence, we model persistence diagrams as Poisson point processes with prior intensities and compute posterior intensities by adopting techniques from the theory of marked point processes. We then propose a family of conjugate prior intensities via Gaussian mixtures to obtain a closed form of the posterior intensity. Finally we demonstrate the utility of this Bayesian framework with a classification problem in materials science using Bayes factors.

Tables

Table 1. Table 1: The parallels between the Bayesian framework for RVs and its counterpart for random PDs.

	Bayesian Framework for RVs	Bayesian Framework for Random PDs
Prior	Modeled by a prior density $f$	Modeled by a Poisson PP with prior intensity $\lambda$
Likelihood	Depends on observed data	Stochastic kernel that depends on observed PDs
Posterior	Compute the posterior density	A Poisson PP with posterior intensity

Table 2. Table 2: List of Gaussian mixture parameters of the prior intensities in Equation ( 5 ). The means μ i 𝒟 X superscript subscript 𝜇 𝑖 subscript 𝒟 𝑋 \mu_{i}^{\mathcal{D}_{X}} are 2 × 1 2 1 2\times 1 vectors and the rest are scalars

	${\mu }_i^{{𝒟}_X}$	${\sigma }_i^{{𝒟}_X}$	$c_i^{{𝒟}_X}$
Informative Prior	$(0.5,1.2)$	$0.01$	$1$
Weakly informative Prior	$(0.5,1.2)$	$0.2$	$1$
Unimodal Uninformative Prior	$(1,1)$	$1$	$1$
Bimodal Uninformative Prior	$(0.5,0.5)$ $(1.5,1.5)$	$0.2$ $0.2$	$1$ $2$

Table 3. Table 3:                                                                    TABLE 3 Parameters for (M3 ′ ) in Equation ( 6 ) and ( 7 ). We set the weight and mean of the Gaussian component, c 𝒟 Y S = 1 superscript 𝑐 subscript 𝒟 subscript 𝑌 𝑆 1 c^{\mathcal{D}_{Y_{S}}}=1 and μ 𝒟 Y S = ( 0.5 , 0 ) superscript 𝜇 subscript 𝒟 subscript 𝑌 𝑆 0.5 0 \mu^{\mathcal{D}_{Y_{S}}}=(0.5,0) respectively for all of the cases. The first row corresponds to parameters in the functions characterizing 𝒟 Y subscript 𝒟 𝑌 \mathcal{D}_{Y} that are used in computing the posterior depicted in the first column of Figure 8 . The second row corresponds to analogous parameters that are used in computing the posterior depicted in the second columns of Figures 5 – 8 . Similarly, the third row corresponds to parameters in the functions characterizing 𝒟 Y subscript 𝒟 𝑌 \mathcal{D}_{Y} used for computing the posterior presented in the third columns of Figures 5 – 8 .

Case-I	Case-II	Case-III	Case-IV
			${\sigma }^{{𝒟}_{Y_O}}=0.1$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$
${\sigma }^{{𝒟}_{Y_O}}=0.01$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$	${\sigma }^{{𝒟}_{Y_O}}=0.1$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$	${\sigma }^{{𝒟}_{Y_O}}=0.01$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$	${\sigma }^{{𝒟}_{Y_O}}=0.1$ ${\sigma }^{{𝒟}_{Y_S}}=1$
${\sigma }^{{𝒟}_{Y_O}}=0.1$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$	${\sigma }^{{𝒟}_{Y_O}}=0.1$ ${\sigma }^{{𝒟}_{Y_S}}=1$	${\sigma }^{{𝒟}_{Y_O}}=0.1$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$	${\sigma }^{{𝒟}_{Y_O}}=0.01$ ${\sigma }^{{𝒟}_{Y_S}}=0.1$

Table 4. Table 8: FIG. 8: Case-IV: The first, second and third columns match the parameters of observed persistence diagram 𝒟 Y subscript 𝒟 𝑌 \mathcal{D}_{Y} used in computing the posteriors depicted in the third column of Figure 5 and Figure 6 , and second column of Figure 7l respectively with α = 0.5 𝛼 0.5 \alpha=0.5 . The parameters are presented in Table 3 . The color maps are representing scaled intensities. A variation in the level of intensity is observed for all of them compared to their respective cases due to the added term in the posterior intensity. The posterior intensities in first and second columns exhibits the estimation of hole with higher variability as compared to the respective figures in Case-I and Case-II. The posteriors in the third column demonstrate dominance of the prior relative to their corresponding figures in Case-III, especially when one examines those for informative, weakly informative and bimodal uninformative priors.

	Case-I	Case-II	Case-III
Informative
Weakly Informative
Unimodal Uninformative
Bimodal Uninformative

Table 5. Table 4:                                                                    TABLE 4 Parameters for the prior intensities used in cross-validation of materials science data. Each prior λ 𝜆 \lambda is indexed by its corresponding class for Prior-1 or U 𝑈 U in the case of the Prior-2. The summary of AUCs across 10-folds for materials science data after scoring with Algorithm 1 is presented in the last three columns.

Priors	Parameters for Prior Intensities				Summary of AUC
		${\mu }_i^{𝒟}$	${\sigma }_i^{𝒟}$	$c_i^{𝒟}$	5th percentile	Mean	95th percentile
Prior-1	${\lambda }_{BCC}$	(0.5,0.24) (3.6,3.6) (3.7,0.65)	2 2 2	1 1 1	$0.931$	$0.941$	$0.958$
	${\lambda }_{FCC}$	(0.4,0.27) (2.8,1.2) (2.9,3)	2 2 2	1 1 1
Prior-2	${\lambda }_U$	(1,1)	20	1	$0.928$	0.94	$0.951$

Equations

$$j_{n|M}(x_1,\dots ,x_n)=\{\begin{array}{l}{\sum }_{\pi \in {\mathcal{S}}_n}{\prod }_{i=1}^{n}p(x_i|y_{\pi (i)}),n=m,\\ 0,\text{otherwise},\end{array}$$

$$G(h)=J_0+\sum _{n=1}^{\infty }\frac{1}{n!}{\int }_{{\mathbb{X}}^n}\left(\prod _{j=1}^{n}h(x_j)\right){\mathbb{J}}_n(d{x}_1\dots d{x}_n)$$

$${\lambda }_{{\mathcal{D}}_X|D_{Y^{1:m}}}(x)=\left(1-\alpha (x)\right){\lambda }_{{\mathcal{D}}_X}(x)+\frac{1}{m}\alpha (x)\sum _{i=1}^m\sum _{y\in D_{Y^i}}\frac{\ell (y|x){\lambda }_{{\mathcal{D}}_X}(x)}{{\lambda }_{{\mathcal{D}}_{Y_S}}(y)+{\int }_{\mathbb{W}}\ell (y|u)\alpha (u){\lambda }_{{\mathcal{D}}_X}(u)du}\text{a.s.}$$

$${\mathcal{N}}^{\ast }(z;\upsilon ,\sigma I):=\mathcal{N}(z;\upsilon ,\sigma I)\text{mathbbm}{1}_{\mathbb{W}}(z),$$

$${\lambda }_{{\mathcal{D}}_X}(x)=\sum _{j=1}^N{c}_j^{{\mathcal{D}}_X}{\mathcal{N}}^{\ast }(x;{\mu }_j^{{\mathcal{D}}_X},{\sigma }_j^{{\mathcal{D}}_X}I),\text{vspace}-0.1in$$

$$\ell (y|x)={\mathcal{N}}^{\ast }(y;x,{\sigma }^{{\mathcal{D}}_{Y_O}}I).$$

$${\lambda }_{{\mathcal{D}}_{Y_S}}(y)=\sum _{k=1}^M{c}_k^{{\mathcal{D}}_{Y_S}}{\mathcal{N}}^{\ast }(y;{\mu }_k^{{\mathcal{D}}_{Y_S}},{\sigma }_k^{{\mathcal{D}}_{Y_S}}I),$$

$${\lambda }_{{\mathcal{D}}_X|D_{Y^{1:m}}}(x)=(1-\alpha ){\lambda }_{{\mathcal{D}}_X}(x)+\frac{\alpha }{m}\sum _{i=1}^m\sum _{y\in {\mathcal{D}}_{Y^i}}\sum _{j=1}^N{C}_j^y{\mathcal{N}}^{\ast }(x;{\mu }_j^y,{\sigma }_j^{y}I),$$

$$\text{where}{C}_j^y$$

$$w_j^y$$

$$\text{and}{\mu }_j^y$$

$$\sum_{j=1}^{N}\Big{[}\frac{c_{j}^{\mathcal{D}_{X}}\mathcal{N}(y;\mu_{j}^{\mathcal{D}_{X}},(\sigma^{\mathcal{D}_{Y_{O}}}+\sigma_{j}^{\mathcal{D}_{X}})I)}{\lambda_{\mathcal{D}_{Y_{S}}}(y)+\alpha\sum_{j=1}^{N}c_{j}^{\mathcal{D}_{X}}\mathcal{N}(y;\mu_{j}^{\mathcal{D}_{X}},(\sigma^{\mathcal{D}_{Y_{O}}}+\sigma_{j}^{\mathcal{D}_{X}})I)\int_{\mathbb{W}}\mathcal{N}(u;\mu_{j}^{x|y},\sigma_{j}^{x|y}I)du}\Big{]}\mathcal{N}^{*}(x;\mu_{j}^{x|y},\sigma_{j}^{x|y}I),$$

$$p_{\mathcal{D}}(D)=\frac{e^{-\lambda }}{|D|!}\prod _{d\in D}{\lambda }_{\mathcal{D}}(d)=\frac{e^{-\lambda }}{|D|!}\prod _{d\in D}\sum _{i=1}^N{c}_i^{\mathcal{D}}{\mathcal{N}}^{\ast }(d;{\mu }_i^{\mathcal{D}},{\sigma }_i^{\mathcal{D}}I),$$

$$p_{\mathcal{D}|\mathcal{D}_{Y}}(D|T_{Y})=\frac{e^{-\lambda}}{|D|!}\prod_{d\in D}\lambda_{D|T_{Y}}(d)=\frac{e^{-\lambda}}{|D|!}\prod_{d\in D}\,\Big{[}(1-\alpha)\lambda_{\mathcal{D}}(d)+\frac{\alpha}{n}\sum_{y_{j}\in T_{Y}}\sum_{i=1}^{N}C_{i}^{d|y_{j}}\mathcal{N}(d;\mu_{i}^{d|y_{j}},\sigma_{i}^{d|y_{j}}I)\Big{]},$$

$$BF(D)=\frac{p_{D|{\mathcal{D}}_Y}(D|T_Y)}{p_{D|{\mathcal{D}}_{Y^{\prime }}}(D|T_{Y^{\prime }})}$$

$${\lambda }_{{\mathcal{D}}_X|D_{Y^{1:m}}}={\lambda }_{{\mathcal{D}}_{X_V}|D_{Y^{1:m}}}+{\lambda }_{{\mathcal{D}}_{X_O}|D_{Y^{1:m}}}=(1-\alpha (x)){\lambda }_{{\mathcal{D}}_X}+{\lambda }_{{\mathcal{D}}_{X_O}|D_{Y^{1:m}}},$$

$${\lambda }_{{\mathcal{D}}_{X_O}|D_{Y^{1:m}}}=\frac{1}{m}\sum _{i=1}^m{\lambda }_{{\mathcal{D}}_{X^i}|D_{Y^i}}.$$

$$\mathcal{H}:=\Big{\{}(x,y)\in(\mathcal{D}_{X_{O}},\mathcal{D}_{Y_{O}})\Big{\}}\bigcup\Big{\{}(\Delta,y)|y\in\mathcal{D}_{Y_{S}}\Big{\}}.$$

$${\lambda }_{\mathcal{H}}(x,y)=\alpha (x){\lambda }_{{\mathcal{D}}_X}(x)\ell (y|x)\text{mathbbm}{1}_{x\in \mathbb{W}}+{\lambda }_{D_{Y_S}}(y)\text{mathbbm}{1}_{x=\Delta }.$$

$${\lambda }_{\mathcal{H}}(x,y)={\lambda }_{{\mathcal{H}}_Y}(y)p(x|y).$$

$$p(x|y)=\frac{\alpha (x){\lambda }_{{\mathcal{D}}_X}(x)\ell (y|x)\text{mathbbm}{1}_{x\in \mathbb{W}}+{\lambda }_{D_{Y_S}}(y)\text{mathbbm}{1}_{x=\Delta }}{{\lambda }_{{\mathcal{H}}_Y}(y)},{\lambda }_{{\mathcal{H}}_Y}(y)\ne 0.$$

$${\lambda }_{{\mathcal{H}}_X|D_{Y^i}}(x)=\sum _{y\in D_{Y^i}}\frac{\alpha (x){\lambda }_{{\mathcal{D}}_X}(x)\ell (y|x)\text{mathbbm}{1}_{x\in \mathbb{W}}+{\lambda }_{D_{Y_S}}(y)\text{mathbbm}{1}_{x=\Delta }}{{\lambda }_{{\mathcal{H}}_Y}(y)},{\lambda }_{{\mathcal{H}}_Y}(y)\ne 0$$

$${\lambda }_{{\mathcal{D}}_{X_O}|D_{Y^i}}(x)=\sum _{y\in D_{Y^i}}\frac{\alpha (x){\lambda }_{{\mathcal{D}}_X}(x)\ell (y|x)}{{\lambda }_{{\mathcal{H}}_Y}(y)},a.s.$$

$${\lambda }_{{\mathcal{D}}_{X_O}|D_{Y^i}}(x)=\alpha (x)\sum _{y\in D_{Y^i}}\frac{\ell (y|x){\lambda }_{{\mathcal{D}}_X}(x)}{{\lambda }_{{\mathcal{D}}_{Y_S}}(y)+{\int }_{\mathbb{W}}\ell (y|u)\alpha (u){\lambda }_{{\mathcal{D}}_X}(u)du},$$

Full text

Keywords Bayesian inference and classification, intensity, marked Poisson point processes, topological data analysis, high entropy alloys, atom probe tomography.

1 Introduction

A crucial first step in understanding patterns and properties of a crystalline material is determining its crystal structure. For highly disordered metallic alloys, such as high entropy alloys (HEAs), atom probe tomography (APT) gives a snapshot of the local atomic environment; see Figure 1. However, APT has two main drawbacks: experimental noise and an abundance of missing data. Approximately 65% of the atoms in a sample are not registered in a typical experiment [50], and the spatial coordinates of those identified atoms are corrupted by experimental noise [42]. Understanding the atomic pattern within HEAs using an APT image requires observation of atomic cubic unit neighborhood cells under a microscope. This is problematic as APT may have a spatial resolution approximately the length of the unit cell under consideration [28, 42]. Hence, the process is unable to see the finer details of a material, rendering the determination of a lattice structure a challenging problem [55, 38]. Existing algorithms for detecting the crystal structure [13, 24, 25, 32, 43, 56] are not able to establish the crystal lattice of an APT dataset, as they rely on symmetry arguments based on identifying repeating parts of molecules. Consequently, the field of atom probe crystallography, i.e., determining the crystal structure from APT data, has emerged in recent years [21] and [43]. Algorithms in this field rely on knowing the global lattice structure a priori and aim to determine local small-scale structures within a larger sample. For some materials this information is readily known, while for others, such as HEAs, the global structure is unknown and must be inferred.

A recent work [62] proposes a machine-learning approach to classifying crystal structures of a noisy and sparse materials dataset without knowing the global structure a priori. The authors employ a convolutional neural network for classifying the crystal structure by examining a diffraction image, a computer-generated diffraction pattern. The authors suggest their method could be used to determine the crystal structure of APT data. However, the synthetic data considered in [62] is not a realistic representation of experimental APT data, where about 65% of the data is missing and furthermore corrupted by observational noise. Most importantly, their synthetic data is either sparse or noisy, not a combination of both. Herein, we consider a combination of noise and sparsity such as is the case in real APT data.

In this work, we specifically classify unit cells that are either body-centered cubic (BCC) or face-centered cubic (FCC). These lattice structures are the essential building blocks of HEAs [61] and have fundamental differences that set them apart in the case of noise-free, complete materials data. The BCC structure has a single atom in the center of the cube, while the FCC has a void in its center but has atoms on the center of the cubes’ faces, see Figure 1 (b-c). These two crystal structures are distinct when viewed through the lens of topology. Differentiating between the empty space and connectedness of these two lattice structures allows us to create an accurate classification rule. This fundamental distinction between BCC and FCC point clouds is captured well by topological methods and explains the high degree of accuracy in the classification scheme presented herein. Indeed, we offer a Bayesian classification framework for persistence homology.

Overall, topological data analysis (TDA) encompasses a broad set of techniques that explore topological structure in datasets [17, 22, 12, 59]. One of these techniques, persistent homology, associates shapes to data and summarizes salient features with persistence diagrams – multisets of points that represent homological features along with their appearance and disappearance scales [17]. Features of a persistence diagram that exhibit long persistence describe global topological properties in the underlying dataset, while those with shorter persistence encode information about local geometry and/or noise. Hence, persistence diagrams can be considered multiscale summaries of data’s shape. While there are several methods present in literature to compute persistence diagrams, we adopt geometric complexes that are typically used for applications of persistent homology to data analysis in various settings such as handwriting analysis [2], studying of brain arteries [4, 5], image analysis [7, 11, 10], neuroscience [14, 54, 3], sensor network [16, 52], protein structure [20, 31], biology [51, 39, 45], dynamical system [29], action recognition [58], signal analysis [35, 34, 47, 36], chemistry [60], genetics [26], object data [46], etc.

Researchers desire to utilize persistence diagrams for inference and classification problems. Several achieve this directly with persistence diagrams [40, 35, 6, 19, 41, 49, 9], while others elect to first map them into a Hilbert space [8, 48, 1, 57, 18]. The latter approach enables one to adopt traditional machine learning and statistical tools such as principal component analysis, random forests, support vector machines, and more general kernel-based learning schemes. Despite progress toward statistical inference, to the best of our knowledge, a full Bayesian treatment predicated upon creating posterior distributions of persistence diagrams is still absent in literature. The first Bayesian considerations in a TDA context take place in [41] where the authors discuss a conditional probability setting on persistence diagrams where the likelihood for the observed point cloud has been substituted by the likelihood for its associated topological summary.

The homological features in persistence diagrams have no intrinsic order implying they are random sets as opposed to random vectors. This viewpoint is embraced in [40] to construct a kernel density estimator for persistence diagrams. This kernel density estimator gives a sensible way to obtain priors for distributions of persistence diagrams; however, computing posteriors entirely through the random set analog of Bayes’ rule is computationally intractable in general settings [23]. Intuitively, this follows because evaluation of fixed sets in the likelihoods for parametric densities of random sets may involve a term for each possible association of input points to parameters, resulting in exponential scaling with respect to the number of parameters. To address this, we model random persistence diagrams as Poisson point processes. The defining feature of these point processes is that they are solely characterized by a single parameter known as the intensity. Utilizing the theory of marked point processes, we obtain a method for computing posterior intensities that does not require us to consider explicit maps between input diagrams and underlying parameters, alleviating the computational burden associated with deriving the posterior intensity from Bayes’ rule alone.

In particular, for a given collection of observed persistence diagrams, we consider the underlying stochastic phenomena generating persistence diagrams to be Poisson point processes with prior uncertainty captured in presupposed intensities. In applications, one may select an informative prior by choosing an intensity based on expert opinion, or alternatively choose an uninformative prior intensity when information is not available. The likelihood functions in our model represent the level of belief that observed diagrams are representative of the entire population. We build this analog using the theory of marked Poisson point processes [15]. A central idea of this paper is to utilize the topological summaries of point clouds in place of the actual point clouds. This provides a powerful tool with applications in wide ranging fields. The application considered in this paper is the classification of the crystal structure of materials, which allows scientists to predict the properties of a crystalline material. Our goal is to view point clouds through their topological descriptors as this can reveal essential shape peculiarities latent in the point clouds. Our Bayesian method adopts a substitution likelihood technique by Jeffreys in [27] instead of considering the full likelihood for the point cloud. A similar sort of discussion was considered in [41] for defining conditional probability on persistence diagrams.

Another key contribution of this paper is the derivation of a closed form of the posterior intensity, which relies on conjugate families of Gaussian mixtures. An advantage of this Gaussian mixture representation is that it allows us to perform Bayesian inference in an efficient and reliable manner. Indeed, this model can be viewed as an analog of the ubiquitous example in standard Bayesian inference where a Gaussian prior and likelihood yield a Gaussian posterior. We present a detailed example of our closed form implementation to demonstrate computational tractability and showcase its applicability by using it to build a Bayes factor classification algorithm; we test the latter in a classification problem for materials science data.

The contributions of this work are:

Theorem 3.1, which provides the Bayesian framework for computing the posterior distribution of persistence diagrams. 2. 2.

Proposition 3.1, which yields a conjugate family of priors based on a Gaussian mixture for the proposed Bayesian framework. 3. 3.

A classification scheme using Bayes factors considering the posteriors of persistence diagrams and its application to a materials science problem.

This paper is organized as follows. Section 2 provides a brief overview of persistence diagrams and general point processes. Our methods are presented in Section 3. In particular, Subsection 3.1 establishes the Bayesian framework for persistence diagrams, while Subsection 3.2 contains the derivation of a closed form for a posterior distribution based on a Gaussian mixture model. A classification algorithm with Bayes factors is discussed in Section 4. To assess the capability of our algorithm, we investigate its performance on materials data in Subsection 4.1. Finally, we end with discussions and conclusions in Section 5.

2 Background

We begin by discussing preliminary definitions essential for building our model. In Subsection 2.1, we briefly review simplicial complexes and provide a formal definition for persistence diagrams (PDs). Pertinent definitions and theorems from point processes (PPs) are discussed in Subsection 2.2 .

2.1 Persistence Diagrams

We start by discussing simplices and simplicial complexes, intermediary structures for constructing PDs.

Definition 2.1**.**

A $\mathbbm{d}$-dimensional collection of data $\{v_{0},\dots,v_{n}\}\subset\mathbb{R}^{\mathbbm{d}}\setminus\{0\}$ is said to be geometrically independent if for any set $t_{i}\in\mathbb{R}$ with $\sum_{i=0}^{n}t_{i}=0$, the equation $\sum_{i=0}^{n}t_{i}v_{i}=0$ implies that $t_{i}=0$ for all $i\in\{0,\dots,n\}.$

Definition 2.2**.**

A $k-$simplex, is a collection of $k+1$ geometrically independent elements along with their convex hull: [v_{0},\dots,v_{k}]=\Big{\{}\sum_{i=0}^{k}\alpha_{i}v_{i}|\sum_{i=0}^{k}\alpha_{i}=1\Big{\}}. We say that the vertices $v_{0},\dots,v_{n}$ span the $k-$dimensional simplex, $[v_{0},\dots,v_{k}]$. The faces of a $k-$simplex $[v_{0},\dots,v_{k}]$, are the $(k-1)-$simplices spanned by subsets of $\{v_{0},\dots,v_{k}\}$.

Definition 2.3**.**

A simplicial complex $S$ is a collection of simplices satisfying two conditions: (i) if $\xi\in S$, then all faces of $\xi$ are also in $S$, and (ii) the intersection of two simplices in $S$ is either empty or contained in $S$.

Given a point cloud $X$, our goal is to construct a sequence of simplicial complexes that reasonably approximates the underlying shape of the data. We accomplish this by using the Vietoris-Rips filtration.

Definition 2.4**.**

Let $X=\{x_{i}\}_{i=0}^{L}$ be a point cloud in $\mathbb{R}^{\mathbbm{d}}$ and $r>0$. The Vietoris-Rips complex of $X$ is defined to be the simplicial complex $\mathcal{V}_{r}(X)$ satisfying $[x_{i_{1}},\dots,x_{i_{l}}]\in\mathcal{V}_{r}(X)$ if and only if $\operatorname{diam}(x_{i_{1}},\dots,x_{i_{l}})<r$. Given a nondecreasing sequence $\{r_{n}\}\in\mathbb{R}^{+}\cup\{0\}$ with $r_{0}=0$, we denote its Vietoris-Rips filtration by $\{\mathcal{V}_{r_{n}}(X)\}_{n\in\mathbb{N}}$.

A persistence diagram $\mathcal{D}$ is a multiset of points in $\mathcal{W}:=\mathbb{W}\times\{0,1,\dots,\mathbbm{d}-1\}$, where $\mathbb{W}:=\{(b,d)\in\mathbb{R}^{2}|\,\,d\geq b\geq 0\}$ and each element $(b,d,k)$ represents a homological feature of dimension $k$ that appears at scale $b$ during a Vietoris-Rips filtration and disappears at scale $d$. Intuitively speaking, the feature $(b,d,k)$ is a $k-$dimensional hole lasting for duration $d-b$. Namely, features with $k=0$ correspond to connected components, $k=1$ to loops, and $k=2$ to voids. An example of a PD is shown in Figure 2.

2.2 Poisson Point Processes

This section contains basic definitions and fundamental theorems from PPs, primarily Poisson PPs. Detailed treatments of Poisson PPs can be found in [15] and references therein. For the remainder of this section, we take $\mathbb{X}$ and $\mathcal{X}$ to be a Polish space and its Borel $\sigma$-algebra, respectively.

Definition 2.5**.**

A finite point process $\mathcal{P}$ is a pair $\left(\left\{p_{n}\right\},\left\{\mathbb{P}_{n}\right\}\right)$ where $\sum_{n=0}^{\infty}p_{n}=1$ and $\mathbb{P}_{n}$ is a symmetric probability measure on $\mathcal{X}^{n}$, where $\mathcal{X}^{0}$ is understood to be the trivial $\sigma$-algebra.

The sequence $\{p_{n}\}$ defines a cardinality distribution and the measures $\{\mathbb{P}_{n}\}$ give spatial distributions of vectors $(x_{1},\dots,x_{n})$ for fixed $n$. Definition 2.5 naturally prescribes a method for sampling a finite PP: (i) determine the number of points $n$ by drawing from $\{p_{n}\}$ then, (ii) spatially distribute $(x_{1},\dots,x_{n})$ according to a draw from $\mathbb{P}_{n}$. As PPs model random collections of elements in $\{x_{1},\dots,x_{n}\}\subset\mathbb{X}$ whose order is irrelevant, any sensible construction relying on random vectors should assign equal weight to all permutations of $(x_{1},\dots,x_{n})$. This is ensured by the symmetry requirement in Definition 2.5. We abuse notation and write $\mathcal{P}$ for samples from $\mathcal{P}$ as well as their set representations. It proves useful to describe finite PPs by a set of measures that synthesize $p_{n}$ and $\mathbb{P}_{n}$ to simultaneously package cardinality and spatial distributions.

Definition 2.6**.**

Let $\left(\left\{p_{n}\right\},\left\{\mathbb{P}_{n}\right\}\right)$ be a finite PP. The Janossy measures $\left\{\mathbb{J}_{n}\right\}$ are defined as the set of measures satisfying $\mathbb{J}_{n}=n!p_{n}\mathbb{P}_{n},\hskip 5.69054pt\text{for all}\hskip 5.69054ptn\in\mathbb{N}.$

Given a collection of disjoint rectangles $A_{1},\dots,A_{n}\subset\mathbb{X}$, the value $\mathbb{J}_{n}(A_{1}\times\dots\times A_{n})$ is the probability of observing exactly one element in each of $A_{1},\dots,A_{n}$ and none in the complement of their union. For applications, we are primarily interested in Janossy measures $\mathbb{J}_{n}$ that admit densities $j_{n}$ with respect to a reference measure on $\mathbb{X}$. We are now ready to describe the class of finite PPs that model PDs.

Definition 2.7**.**

Let $\Lambda$ be a finite measure on $\mathbb{X}$ and $\mu:=\Lambda(\mathbb{X})$. 7The finite point process $\Pi$ is Poisson if, for all $n\in\mathbb{N}$ and measurable rectangles $A_{1}\times\dots\times A_{n}\in\mathcal{X}^{n}$, $p_{n}=e^{-\mu}\frac{\mu^{n}}{n!},$ and $\mathbb{P}_{n}(A_{1}\times\dots\times A_{n})=\prod_{i=1}^{n}\left(\frac{\Lambda(A_{i})}{\mu}\right).$We call $\Lambda$ an intensity measure.

Equivalently, a Poisson PP is a finite PP with Janossy measures $\mathbb{J}_{n}(A_{1}\times\dots\times A_{n})=e^{-\mu}\prod_{i=1}^{n}\Lambda(A_{i}).$ The intensity measure in Definition 2.7 admits a density, $\lambda$, with respect to some reference measure on $\mathbb{X}$. Notice that for all $A\in\mathcal{X}$, $\mathbb{E}(|\Pi\cap A|)=\sum_{n=0}^{\infty}p_{n}\mathbb{E}_{\mathbb{P}_{n}}\left(\sum_{k=0}^{n}k\binom{n}{k}\mathbbm{1}_{A^{k}\times(A^{c})^{n-k}}\right)$. Elementary calculations then show $\mathbb{E}(|\Pi\cap A|)=\Lambda(A)$. Thus, we interpret the intensity measure of a region $A$, $\Lambda(A)$ as the expected number of elements in $\Pi$ that land in $A$. The intensity measure serves as an analog to the first order moment for a random variable.

The next two definitions involve a joint PP wherein points from one space parameterize distributions for the points living in another. Consequently, we introduce another Polish space $\mathbb{M}$ along with its Borel $\sigma$-algebra $\mathcal{M}$ to serve as the mark space in a marked Poisson PP. These model scenarios in which points drawn from a Poisson PP provide a data likelihood model for Bayesian inference with PPs.

Definition 2.8**.**

Suppose $\ell:\mathbb{X}\times\mathbb{M}\rightarrow\mathbb{R}^{+}\cup\{0\}$ is a function satisfying: 1) for all $x\in\mathbb{X}$, $\ell(x,\bullet)$ is a probability measure on $\mathbb{M}$, and 2) for all $B\in\mathcal{M}$, $\ell(\bullet,B)$ is a measurable function on $\mathbb{X}$. Then, $\ell$ is a stochastic kernel from $\mathbb{X}$ to $\mathbb{M}$.

Definition 2.9**.**

A marked Poisson point process $\Pi_{M}$ is a finite point process on $\mathbb{X}\times\mathbb{M}$ such that: (i) $\left(\left\{p_{n}\right\},\left\{\mathbb{P}_{n}(\bullet\times\mathbb{M})\right\}\right)$ is a Poisson PP on $\mathbb{X}$, and (ii) for all $(x_{1},\dots,x_{n})\in\mathbb{X}^{n}$, measurable rectangles $B_{1}\times\dots\times B_{n}\in\mathcal{M}^{n}$ , $\mathbb{P}_{n}((x_{1},\dots,x_{n})\times B_{1}\times\dots\times B_{n})=\frac{1}{n!}\sum_{\pi\in\mathcal{S}_{n}}\prod_{i=1}^{n}\ell(x_{\pi(i)},B_{i})$, where $\mathcal{S}_{n}$ is the set of all permutations of $(1,\dots,n)$ and $\ell$ is a stochastic kernel.

Given a set of observed marks $M=\ \{y_{1},\dots,y_{m}\}$ it can be shown [53] that the Janossy densities for the PP induced by $\Pi_{M}$ on $\mathbb{X}$ given M are

[TABLE]

where $p$ is the stochastic kernel for $\Pi_{M}$ evaluated in $\mathbb{X}$ for a fixed value of $y\in\mathbb{M}$.

The following theorems allow us to construct new Poisson PPs from existing ones. Their proofs can be found in [30].

Theorem 2.1** (The Superposition Theorem).**

Let $\{\Pi_{n}\}_{n\in\mathbb{N}}$ be a collection of independent Poisson PPs each having intensity measure $\Lambda_{n}$. Then their superposition $\Pi$ given by $\Pi:=\bigcup_{n\in\mathbb{N}}\Pi_{n}$ is a Poisson PP with intensity measure $\Lambda=\sum_{n\in\mathbb{N}}\Lambda_{n}$.

Theorem 2.2** (The Mapping Theorem).**

Let $\Pi$ be a Poisson PP on $\mathbb{X}$ with $\sigma$-finite intensity measure $\Lambda$ and let $(\mathbb{T},\mathcal{T})$ be a $\sigma$-algebra. Suppose $f:\mathbb{X}\rightarrow\mathbb{T}$ is a measurable function. Write $\Lambda^{*}$ for the induced measure on $T$ given by $\Lambda^{*}(B):=\Lambda(f^{-1}(B))$ for all $B\in\mathcal{T}$. If $\Lambda^{*}$ has no atoms, then $f\circ\Pi$ is a Poisson PP on $T$ with intensity measure $\Lambda^{*}$.

Theorem 2.3** (The Marking Theorem).**

The marked Poisson PP in Definition 2.9 has the intensity measure given by $\Lambda_{M}(C)=\iint_{C}\Lambda(dx)\ell(x,dm)$, where $\Lambda$ is the intensity measure for the Poisson PP that $\Pi_{M}$ induces on $\mathbb{X}$, and $\ell$ is a stochastic kernel.

The final tool we need is the probability generating functional as it enables us to recover intensity measures using a notion of differentiation. The probability generating functional can be interpreted as the PP analog of the probability generating function.

Definition 2.10**.**

Let $\mathcal{P}$ be a finite PP on a Polish space $\mathbb{X}$. Denote by $\mathcal{B}(\mathbb{C})$ the set of all functions $h:\mathbb{X}\rightarrow\mathbb{C}$ with $||h||_{\infty}<1$. The probability generating functional of $\mathcal{P}$ denoted $G:\mathcal{B}(\mathbb{C})\rightarrow\mathbb{R}$ is given by

[TABLE]

Definition 2.11**.**

Let $G$ be the probability generating functional given in Equation (2). The functional derivative of $G$ in the direction of $\eta$ evaluated at $h$, when it exists, is given by $G^{\prime}(h;\eta)=\lim_{\epsilon\rightarrow 0}\frac{G(h+\epsilon\eta)-G(h)}{\epsilon}$.

It can be shown that the functional derivative satisfies the familiar product rule [33]. As is proved in [44], the intensity measure $\Lambda$ of the Poisson PP in Definition 2.7 can be obtained by differentiating $G$, i.e., $\Lambda(A)=G^{\prime}(1;\mathbbm{1}_{A})$, where $\mathbbm{1}_{A}$ is the indicator function for any $A\in\mathcal{X}$. Generally speaking, one obtains the intensity measure for a general point process through $\Lambda(A)=\lim_{h\rightarrow 1}G^{\prime}(h;\mathbbm{1}_{A})$, but the preceding identity suffices for our purposes since we only consider point processes for which Equation (2) is defined for all bounded $h$.

Corollary 2.1**.**

The intensity function for the PP whose Janossy densities are listed in Equation (1) is $\sum_{i=1}^{m}p(x|y_{i})$.

Proof.

This directly follows from writing the probability generating functional for the PP in question using its Janossy densities then applying $\Lambda(A)=G^{\prime}(1;\mathbbm{1}_{A})$. By Definition 2.10, linearity of the integral, and Fubini’s theorem, we have G(h)=\prod_{i=1}^{m}\Big{(}\int_{\mathbb{X}}h(x)p(x|y_{i})\operatorname{d}x\Big{)}=\prod_{i=1}^{m}G_{i}(h) where $G_{i}(h)$ is the probability generating functional for the PP with Janossy densities $j_{1}(x)=p(x|y_{i})$ and $j_{n}=0$ for $n\neq 1$. One arrives at the desired result by applying the product rule for functional derivatives and the intensity retrieval property of probability generating functionals. ∎

3 Bayesian Inference

In this section, we construct a framework for Bayesian inference with PDs by modeling them as Poisson PPs. First, we derive a closed form for the posterior intensity given a PD drawn from a finite PP, and then we present a family of conjugate priors followed by an example.

3.1 Model

Given a persistence diagram $\mathcal{D}$, the map $T:\mathbb{W}\rightarrow T(\mathbb{W})$ given by $T(b,d)=(b,d-b)$ defines tilted representation of $\mathcal{D}$ as $T(\mathcal{D})=\cup_{(b,d,k)\in\mathcal{D}}(T(b,d),k)$; see Figure 2. In the sequel, we assume all PDs are given in their tilted representations and, unless otherwise noted, abuse notation by writing $\mathbb{W}$ and $\mathcal{D}$ for $T(\mathbb{W})$ and $T(\mathcal{D})$, respectively. We also fix the homological dimension of features in a PD by defining $\mathcal{D}^{k}:=\{(b,d)\in\mathbb{W}\hskip 5.69054pt|\hskip 5.69054pt(b,d,k)\in\mathcal{D}\}$.

According to Bayes’ theorem, posterior density is proportional to the product of a likelihood function and a prior. To adopt Bayesian framework to PDs, we need to define two models. In particular, our Bayesian framework views a random PD as a Poisson PP equipped with a prior intensity while observed PDs $\mathcal{D}_{Y}$ are considered to be marks from a marked Poisson PP. This enables modification of the prior intensity by incorporating observed PDs, yielding a posterior intensity based on data. Some parallels between our Bayesian framework and that for random variables (RVs) are illustrated in Table 1.

Let $(\mathcal{D}_{X}^{k},\mathcal{D}_{Y}^{k})\in\mathbb{W}\times\mathbb{W}$ be a finite PP and consider the following:

• (M1)

For $k_{1}\neq k_{2}$, $(\mathcal{D}_{X}^{k_{1}},\mathcal{D}_{Y}^{k_{1}})$ and $(\mathcal{D}_{X}^{k_{2}},\mathcal{D}_{Y}^{k_{2}})$ are independent.

• (M2)

For $k$ fixed, $\mathcal{D}_{X}^{k}=\mathcal{D}_{X_{O}}^{k}\cup\mathcal{D}_{X_{V}}^{k}$ and some $\alpha:\mathbb{W}\rightarrow[0,1]$, $\mathcal{D}_{X_{O}}^{k}$ and $\mathcal{D}_{X_{V}}^{k}$ are independent Poisson PPs having intensity functions $\alpha(x)\lambda_{\mathcal{D}_{X}^{k}}(x)$ and $\left(1-\alpha(x)\right)\lambda_{\mathcal{D}_{X}^{k}}(x)$, respectively.

• (M3)

For $k$ fixed, $\mathcal{D}_{Y}^{k}=\mathcal{D}_{Y_{O}}^{k}\cup\mathcal{D}_{Y_{S}}^{k}$ where

• (i)

$(\mathcal{D}_{X_{O}}^{k},\mathcal{D}_{Y_{O}}^{k})$ is a marked Poisson PP with a stochastic kernel density $\ell(y|x)$.

• (ii)

$\mathcal{D}_{Y_{O}}^{k}$ and $\mathcal{D}_{Y_{S}}^{k}$ are independent finite Poisson PPs where $\mathcal{D}_{Y_{S}}^{k}$ has intensity function $\lambda_{\mathcal{D}_{Y_{S}}^{k}}$.

Hereafter we abuse notation by writing $\mathcal{D}_{X}$ for $\mathcal{D}_{X}^{k}$. The modeling assumption (M1) allows us to develop results independently for each homological dimension $k$ then combine them using independence. In (M2), the random persistence diagram $\mathcal{D}_{X}$ modeled as a Poisson PP with prior intensity $\lambda_{\mathcal{D}_{X}}$. There are two cases we may encounter for any point $x$ from the prior intensity due to the nature of persistence diagrams. We assign a probability function $\alpha(x)$ to accommodate these two possibilities. Depending upon the noise level in data, any feature $x$ in $\mathcal{D}_{X}$ may not be represented in observations and this scenario happens with probability $1-\alpha(x)$ and we denote this case as $\mathcal{D}_{X_{V}}$ in (M2). Otherwise a point $x$ observed with a probability of $\alpha(x)$ and this scenario is presented as $\mathcal{D}_{X_{O}}$ in (M2). Consequently, the intensities of $\mathcal{D}_{X_{O}}$ and $\mathcal{D}_{X_{V}}$ are proportional to the intensity $\lambda_{\mathcal{D}_{X}}$ weighted by $\alpha(x)$ and $1-\alpha(x)$ respectively and the total prior intensity for $\mathcal{D}_{X}$ is given by their sum. (M3) considers observed persistence diagram $\mathcal{D}_{Y}$ and decomposes it into two independent PDs, $\mathcal{D}_{Y_{O}}$ and $\mathcal{D}_{Y_{S}}$. $\mathcal{D}_{Y_{O}}$ is linked to $\mathcal{D}_{X_{O}}$ via a marked point process with likelihood $\ell(y|x)$ defined in Equation (1), whereas the component $\mathcal{D}_{Y_{S}}$ includes any point $y$ that arises from noise or unanticipated geometry. See Figure 3 for a graphical representation of these ideas.

Theorem 3.1** (Bayesian Theorem for Persistence Diagrams).**

Let $\mathcal{D}_{X}$ be a persistence diagram modeled by a Poisson PP as in (M2). Suppose $\mathcal{D}_{X_{O}}$ and $\mathcal{D}_{X_{V}}$ have prior intensities $\alpha(x)\lambda_{\mathcal{D}_{X}}$ and $(1-\alpha(x))\lambda_{\mathcal{D}_{X}}$, respectively. Consider $D_{Y^{1}},\dots,D_{Y^{m}}$ independent samples from the point process that characterizes the persistence diagram $\mathcal{D}_{Y}$ of (M3) and denote $D_{Y^{1:m}}:=\cup_{i=1}^{m}D_{Y^{i}}$ where $D_{Y^{i}}=D_{Y_{O}^{i}}\cup D_{Y_{S}^{i}}$ for all $i=1,\cdots,m$. Moreover, $\ell(y|x)$ is the likelihood associated for the stochastic kernel between $\mathcal{D}_{X_{O}}$ and $\mathcal{D}_{Y_{O}}$, and $\lambda_{\mathcal{D}_{Y_{S}}}$ is the intensity of $\mathcal{D}_{Y_{S}}$ as defined in (M3). Then, the posterior intensity of $\mathcal{D}_{X}$ given $D_{Y^{1:m}}$ is

[TABLE]

The proof of Theorem 3.1 can be found in the appendix. One important point about the above theorem is that, instead of relying on a likelihood function for the point cloud data, our Bayesian model considers the likleihood for the persistence diagram generated by the observed point cloud data at hand. This is analogous to the idea of substitution likelihood by Jeffreys in [27].

3.2 A Conjugate Family of Prior Intensities: Gaussian Mixtures

This section focuses on constructing a family of conjugate prior intensities, i.e., a collection of priors that yield posterior intensities of the same form when used in Equation (3). Exploiting Theorem 3.1 with Gaussian mixture prior intensities, we obtain Gaussian mixture posterior intensities. As PDs are stochastic point processes on the space $\mathbb{W}$, not $\mathbb{R}^{2}$, we consider a restricted Gaussian density restricted to $\mathbb{W}$. Namely, for a Gaussian density on $\mathbb{R}^{2}$, $\mathcal{N}(z;\upsilon,\sigma I)$, with mean $\upsilon$ and covariance matrix $\sigma I$ , we restrict the Gaussian density on $\mathbb{W}$ as

[TABLE]

where $\mathbbm{1}_{\mathbb{W}}$ is the indicator function of the wedge $\mathbb{W}$. ‘

Consider a random persistence diagram $\mathcal{D}_{X}$ as in (M2) and a collection of observed PDs $\{D_{Y^{1}},\cdots,D_{Y^{m}}\}$ that are independent samples from Poisson PP characterizing the PD $\mathcal{D}_{Y}$ in (M3). We denote $D_{Y^{1:m}}:=\cup_{i=1}^{m}D_{Y^{i}}$. Below we specialize (M2) and (M3) so that applying Theorem 3.1 to a mixed Gaussian prior intensity yields a mixed Gaussian posterior:

1. (M2*′*)

$\mathcal{D}_{X}=\mathcal{D}_{X_{O}}\cup\mathcal{D}_{X_{V}}$, where $\mathcal{D}_{X_{O}}$ and $\mathcal{D}_{X_{V}}$ are independent Poisson PPs with intensities $\alpha\lambda_{\mathcal{D}_{X}}(x)$ and $(1-\alpha)\lambda_{\mathcal{D}_{X}}(x)$, respectively, with

[TABLE]

where $N$ is the number of mixture components. 2. (M3*′*)

$\mathcal{D}_{Y}=\mathcal{D}_{Y_{O}}\cup\mathcal{D}_{Y_{S}}$ where

1. (i)

the marked Poisson PP $(\mathcal{D}_{X_{O}},\mathcal{D}_{Y_{O}})$ has density $\ell({y|x})$ given by

[TABLE] 2. (ii)

$\mathcal{D}_{Y_{O}}$ and $\mathcal{D}_{Y_{S}}$ are independent finite Poisson PPs and $\mathcal{D}_{Y_{S}}$ has intensity function given below.

[TABLE]

where $M$ is the number of mixture components.

Proposition 3.1**.**

Suppose that the assumptions (M1),(M2′), and (M3′) hold; then, the posterior intensity of Equation (3) in Theorem 3.1 is a Gaussian mixture of the form

[TABLE]

[TABLE]

The proof of Proposition 3.1 follows from well known results about products of Gaussian densities given below; for more details, the reader may refer to [37] and references therein.

Lemma 3.1**.**

For $p\times p$ matrices $H,R,P$, with $R$ and $P$ positive definite ,and a $p\times 1$ vector $s$,

$\mathcal{N}(y;Hx,R)\,\mathcal{N}(x;s,P)=q(y)\,\mathcal{N}(x;\hat{s},\hat{P})$, where $q(y)=\mathcal{N}(y;Hs,R+HPH^{T}),\,\,\hat{s}=s+K(y-Hs),\,\,\hat{P}=(I-KH)P$ and $K=PH^{T}(HPH^{T}+R)^{-1}$.

Proof of Proposition 3.1.

Using Lemma 3.1, we first derive $\ell(y|x)\,\lambda_{\mathcal{D}_{X}}(x)$ by observing that, in our model, $H=I,R=\sigma^{\mathcal{D}_{Y_{O}}}I,s={\mu}_{j}^{\mathcal{D}_{X}}$ and $P=\sigma_{j}^{\mathcal{D}_{X}}I$. By typical matrix operations we obtain, $K=\frac{\sigma_{j}^{\mathcal{D}_{X}}}{\sigma_{j}^{\mathcal{D}_{X}}+\sigma^{\mathcal{D}_{Y_{O}}}},\hat{s}=\frac{\sigma_{j}^{\mathcal{D}_{X}}y+\sigma^{\mathcal{D}_{Y_{O}}}\mu_{j}^{\mathcal{D}_{X}}}{\sigma_{j}^{X}+\sigma^{\mathcal{D}_{Y_{O}}}}$, and $\hat{P}=\frac{\sigma^{\mathcal{D}_{Y_{O}}}\sigma_{j}^{\mathcal{D}_{X}}}{\sigma_{j}^{\mathcal{D}_{X}}+\sigma^{\mathcal{D}_{Y_{O}}}}$. Hence the numerator and denominator of the second term in Equation (3), $\sum_{j=1}^{N}c_{j}^{\mathcal{D}_{X}}\mathcal{N}(y;\mu_{j}^{\mathcal{D}_{X}},(\sigma^{\mathcal{D}_{Y_{O}}}+\sigma_{j}^{\mathcal{D}_{X}})I)\,\mathcal{N}^{*}(x;\mu_{j}^{y},\sigma_{j}^{y}I)$, and $\lambda_{\mathcal{D}_{Y_{S}}}(y)+\alpha\sum_{j=1}^{N}c_{j}^{\mathcal{D}_{X}}\mathcal{N}(y;\mu_{j}^{\mathcal{D}_{X}},(\sigma^{\mathcal{D}_{Y_{O}}}+\sigma_{j}^{\mathcal{D}_{X}})I)\int_{\mathbb{W}}\,\mathcal{N}(u;mu_{j}^{y},\sigma_{j}^{y}I)du,$ respectively, yield

[TABLE]

where the bracketed expression is the definition of $C_{j}^{y}$. ∎

3.2.1 Example

Here, we present a detailed example of computing the posterior intensity according to Equation (8) for a range of parametric choices. Reproducing these results, the interested reader may download our R-package BayesTDA. We consider circular point clouds often associated with periodicity in signals [35] and focus on estimating homological features with $k=1$ as they correspond to 1-dimensional holes, which describe the prominent topological feature of a circle. Precisely our goals are to: (i) illustrate posterior intensities and draw analogies to standard Bayesian inference; (ii) determine the relative contributions of the prior and observed data to the posterior; and (iii) perform sensitivity analysis.

We start by considering a Poisson PP with prior intensity $\lambda_{\mathcal{D}_{X}}$ that has the Gaussian mixture form given in (M2*′*). We take into account four types of prior intensities: (i) informative, (ii) weakly informative, (iii) unimodal uninformative, and (iv) bimodal uninformative; see Figures 5–7l (a), (d), (g), (j), respectively. We use one Gaussian component in each of the first three priors as the underlying shape has single $1-$dimensional feature and two for the last one to include a case where we have no information about the cardinality of the underlying true diagram. The parameters of the Gaussian mixture density in Equation (5) used to compute these prior intensities are listed in Table 2. To present the intensity maps uniformly throughout this example while preserving their shapes, we divide the intensities by their corresponding maxima. This ensures all intensities are on a scale from [math] to $1$, and we call it the scaled intensity. The observed PDs are generated from point clouds sampled uniformly from the unit circle and then perturbed by varying levels of Gaussian noise; see Figure 4 wherein we present three point clouds sampled with Gaussian noise having variances $0.001I_{2}$, $0.01I_{2}$, and $0.1I_{2}$, respectively. Consequently, these point clouds provide persistence diagrams $D_{Y^{i}}$ for $i=1,2,3$, which are considered as independent samples from Poisson point process $\mathcal{D}_{Y}$, exhibiting distinctive characteristics such as only one prominent feature with high persistence and no spurious features (Case-I), one prominent feature with high persistence and very few spurious features (Case-II), and one prominent feature with medium persistence and more spurious features (Case-III).

For each observed PD, persistence features are presented as green dots overlaid on their corresponding posterior intensity plots. For Cases-I-III, we set the probability $\alpha$ of the event that a feature in $\mathcal{D}_{X}$ appears in $\mathcal{D}_{Y}$ to $1$, i.e., any feature in $\mathcal{D}_{X}$ is certainly observed through a mark in $\mathcal{D}_{Y}$, and later in Case-IV, we decrease $\alpha$ to $0.5$ while keeping all other parameters the same for the sake of comparison. The choice of $\alpha=0.5$ anticipates that any feature has equal probability to appear or disappear in the observation and in turn provides further intuition about the contribution of prior intensities to the estimated posteriors. We observe that in all cases, the posterior estimates the $1-$dimensional hole; however, with different uncertainty each time. For example, for the cases where the data are trustworthy, expressed by a likelihood with tight variance, or in the case of an informative prior, the posterior accurately estimates the $1-$dimensional hole. In contrast, when the data suffer from high uncertainty and the prior is uninformative, then the posterior offers a general idea that the true underlying shape is a circle, but the exact estimation of the 1-dimensional hole is not accurate. We examine the cases below.

Case-I: We consider informative, weakly informative, unimodal uninformative and bimodal uninformative prior intensities as presented in Figure 5 (a), (d), (g) and (j) respectively to compute corresponding posterior intensities. The prior intensities parameters are listed in Table 2. The observed PD is obtained from the point cloud in Figure 4 (left). The parameters associated to the observed PD are listed in Table 3. For the observed PD arising from data with very low noise, we observe that the posterior computed from any of the priors predicts the existence of a one dimensional hole accurately. Firstly, with a low variability in observed persistence diagram $(\sigma^{\mathcal{D}_{Y_{O}}}=0.01\,\,\,\text{and}\,\,\,\sigma^{\mathcal{D}_{Y_{S}}}=0.1)$, the posterior intensities estimate the hole with high certainty (Figure 5 (b), (e), (h) and (k) respectively). Next, to determine the effect of observed data on the posterior, we increase the variance of the observed PD component $\mathcal{D}_{Y_{O}}$, which consists of features in observed PDs that are associated to the underlying prior. Here, we observe that the posterior intensities still estimate the hole accurately due to the trustworthy data; this is evident in Figure 5 (c), (f), (i) and (l). In Figure 5, the posteriors in (b), (e), (h), and (k) have lower variance around the 1-dimensional feature in comparison to (c), (f), (i), and (l) respectively.

Case-II: Here, we consider all four priors as in Case-I (see Figure 6 (a), (d), (g) and (j)). The point cloud in Figure 4 (center) is more perturbed around the unit circle than that of Case-I (Gaussian noise with variance $0.01I_{2}$). Due to this, the associated PD exhibits spurious features. The parameters used for this case are listed in Table 3. We compute the posterior intensities for each type of prior. First, to illustrate the posterior intensity and check the capability of detecting the $1-$dimensional feature, we use moderate noise for the observed PD $(\sigma^{\mathcal{D}_{Y_{O}}}=0.1\,\,\,\text{and}\,\,\,\sigma^{\mathcal{D}_{Y_{S}}}=0.1)$. The results are presented in Figure 6 (b), (e), (h), and (k); overall, the posteriors estimate the prominent feature with different variances in their respective posteriors. Next, to illustrate the effect of observed data on the posterior, we increase the variance $\sigma^{\mathcal{D}_{Y_{S}}}$ of $\mathcal{D}_{Y_{S}}$. According to our Bayesian model, the persistence diagram component $\mathcal{D}_{Y_{S}}$ contains features that are not associated with $\mathcal{D}_{X}$, so increasing $\sigma^{\mathcal{D}_{Y_{S}}}$ yields that every observed point is linked to $\mathcal{D}_{X}$, and therefore one may expect to observe increased intensity skewed towards the spurious points that arise from noise. Indeed, posterior intensities with weakly informative, unimodal uninformative, and bimodal uninformative priors exhibit the skewness toward the spurious point in Figure 6 (f), (i) and (l) respectively, but this is not the case when an informative prior is used. In (f), we observe increased intensity skewing towards the spurious points, and in (i) and (l) the intensity appears to be bimodal with two modes – one at the prominent and other at the spurious point. For the bimodal uninformative prior since one mode is located close to the spurious point in the observed PD, we observe higher intensity for that mode in the posterior (Figure 6 (l)) with another mode estimating the prominent feature.

Case-III: We again consider the four types of priors here. The observed PD constructed from the point cloud in Figure 4 (right). The point cloud has Gaussian noise with variance $0.1I_{2}$ and due to the high noise level in sampling relative to the unit circle, the associated PD exhibits one prominent feature and several spurious features. We repeat the parameter choices as in Case-I for the variances of observed PD. For the choice of $\sigma^{\mathcal{D}_{Y_{O}}}=0.01$ and $\sigma^{\mathcal{D}_{Y_{S}}}=0.1$, the posteriors computed from all of the four priors are able to detect the difference between the one prominent and other spurious points. We increase the variance of $\mathcal{D}_{Y_{O}}$ to determine the effect of observed PD on the posterior and we observe that only the posterior intensity from informative prior has evidence of the hole (Figure 7l(c)). For the weakly informative and uninformative priors, while the posteriors in (f), (i) and (l) may not detect the hole clearly, in (f) we observe a mode with higher variance and in (i) and (l), a tail towards the high persistence point implying presence of a hole. It should be noted that with the informative prior the posterior intensity identifies the hole closer to the mode of the prior as we increase the variance in $\sigma^{\mathcal{D}_{Y_{O}}}$.

Case-IV: Lastly, in this case we concentrate on the effect of $\alpha$. The rest of the parameters used for this case remain the same and are listed in Table 3. We decrease $\alpha$ to $0.5$ to model the scenario that a feature in $\mathcal{D}_{X}$ has equal probability to appear or vanish in observed $\mathcal{D}_{Y}$. The columns of Figure 8 correspond to the parameters of the observed persistence diagram $\mathcal{D}_{Y}$ used in computing the posteriors depicted in the third column of Figure 5, third column of Figure 6, and second column of Figure 7l respectively. By comparing them with their respective cases, we notice a change in the intensity level in all of these due to the first term of the posterior intensity on the right hand side of Equation (8). Comparing with the respective figures in Case-I, we observe that the posterior intensities are estimating the hole with higher variability for the weakly informative and unimodal uninformative priors. For bimodal prior, we observe bimodality in the posterior. Next for Case-II, the existence of a hole is evident for informative and weakly informative priors with higher uncertainty when compared to their previous cases. The unimodal and bimodal uninformative priors lead to bimodal and trimodal posteriors, respectively. We observe that the posterior resembles the prior intensity more closely when we compare them to respective figures in Case-III. One can especially see this with the informative, weakly informative and bimodal uninformative priors, which have significantly increased intensities at the location of the modes of prior.

4 Classification

The Bayesian framework introduced in this paper allows us to explicitly compute the posterior intensity of a PD given data and prior knowledge. This lays the foundation for supervised statistical learning methods in classification. In this section, we build a Bayes factor classification algorithm based on notions discussed in Section 3 and then apply it on materials data, in particular, on measurements for spatial configurations of atoms.

We commence our classification scheme with a persistence diagram $D$ belonging to an unknown class. We assume that $D$ is sampled from a Poisson point process $\mathcal{D}$ in $\mathbb{W}$ with the prior intensity $\lambda_{\mathcal{D}}$ having the form in (M2*′*). Consequently, its probability density has the form

[TABLE]

where $\lambda=\int_{\mathbb{W}}\lambda_{\mathcal{D}}=\mathbb{E}(|\mathcal{D}|)$, with probability $\alpha$ as in (M2*′*). Next suppose we have two training sets $T_{Y}:=D_{Y_{1:n}}$ and $T_{Y^{\prime}}:=D_{Y^{\prime}_{1:m}}$ from two classes of random diagrams $\mathcal{D}_{Y}$ and $\mathcal{D}_{Y^{\prime}}$, respectively. The likelihood densities of respective classes take the form of Equation (6). We then follow Equation (8) to obtain the posterior intensities of $\mathcal{D}$ given the training sets $T_{Y}$ and $T_{Y^{\prime}}$ from the prior intensities and likelihood densities. In particular, the corresponding posterior probability density of $\mathcal{D}$ given the training set $T_{Y}$ is

[TABLE]

and the posterior probability density given $T_{Y}^{\prime}$ is given by an analogous expression. The Bayes factor defined by

[TABLE]

provides the decision criterion for assigning $D$ to either $\mathcal{D}_{Y}$ or $\mathcal{D}_{Y^{\prime}}$. More specifically, for a threshold $c$, $BF(D)>c$ implies that $D$ belongs to $\mathcal{D}_{Y}$ and $BF(D)<c$ implies otherwise. We summarize this scheme in Algorithm 1 .

4.1 Atom Probe Tomography Data

Our goal in this section is to use Algorithm 1 to classify the crystal lattice of a noisy and sparse materials dataset, where the unit cells are either Body centered cubic (BCC) or Face centered cubic (FCC); recall Figure 1. The BCC structure has a single atom in the center of the cube, while the FCC has a void in its center but has atoms on the centers of the cubes’ faces (Figure 1 (b-c)). However, sparsity and noise do not allow the crystal structure to be revealed. For high-entropy alloys, our object of interest, APT, provides the best atomic level characterization possible. Due to the sparsity and noise in the resulting data, there are only a few algorithms for successfully determining the crystal structure; see [21, 43]. These algorithms, designed for APT data, rely on knowing the global structure a priori (which is not the case for High Entropy Alloys (HEAs)) and seek to discover small-scale structure within a sample.

Bypassing this restriction, the neural network architecture of [62] provides a way to classify the crystal structure of a noisy or sparse dataset by looking at a diffraction image. However, the authors use a test data case with very low noise or sparsity, but not both cases together, which is the representative case of the APT data. The algorithm is also not publicly available, so a side by side comparison of our method with theirs using HEAs is not feasible.

It is natural to consider persistence diagrams in this setting because they distill salient information about the materials patterns with respect to connectedness and empty space (holes) within cubic unit cells, i.e we can differentiate between atomic unit cells by examining their homological features. In particular, after storing both types of spatial configurations as point clouds, we compute their Rips filtrations (see Section 2), collecting resultant 1-dimensional homological features into PDs; see Figure 9. The data set had 200 diagrams from each class. To perform classification with Algorithm 1, we started by specifying priors for each class, $\lambda_{\mathcal{D}_{BCC}}$ and $\lambda_{\mathcal{D}_{FCC}}$. Two scenarios were considered, namely using separate priors (Prior-1 in Table 4) and the same prior (Prior-2 in Table 4) for both BCC and FCC classes. In particular, for Prior-1 we superimpose 50 PDs from each class and find the highly clustered areas by using K-means clustering. The centers of the clusters from K-means are then used as the means in Gaussian mixture priors; see Eqn. (8). In this manner, we produce different priors for BCC and FCC classes. On the other hand for Prior-2 we choose a flat prior with higher variance level than that of Prior-1 for both of the classes. The parameters for these two prior intensities are in Table 4. For all cases, we set $\sigma^{\mathcal{D}_{Y_{O}}}=0.1$ and $\lambda_{\mathcal{D}_{Y_{S}}}(x)=5\mathcal{N}^{*}\left(x;(0,0),0.2I\right)$. We chose a relatively high weight for $\lambda_{\mathcal{D}_{Y_{S}}}$ because the nature of the data implied that extremely low persistence holes were rare events arising from noise. To perform 10-fold cross validation, we partitioned PDs from both classes into training and test sets. During each fold, we took the training sets from each class, $T_{BCC}$ and $T_{FCC}$, and input them into Algorithm 1 as $T_{Y}$ and $T_{Y^{\prime}}$, respectively. Next, we computed the Bayes factor $BF(D)=\frac{p_{\mathcal{D}_{BCC}}(D|T_{BCC})}{p_{\mathcal{D}_{FCC}}(D|T_{FCC})}$ for each diagram $D$ in the test sets. After this, we used the Bayes factors to construct receiver operating characteristic (ROC) curves and computed the resulting areas under the ROC curves (AUCs) . Finally, we used the AUCs from 10-fold cross validation to build a bootstrapped distribution by resampling 2000 times. Information about these bootstrapped distributions is summarized in Table 4, which shows our scoring method almost perfectly distinguishes between the BCC and FCC classes using the Bayesian framework of Section 3. Also, it exemplifies the robustness of our algorithm as two different types of priors produce near perfect accuracy.

5 Discussion and Conclusions

This work is the first approach to introduce a Bayesian framework for persistent homology. This toolbox will give the opportunity to an expert to incorporate their prior belief about the data as well as analyze the data using topological data analysis methods. To that end, we introduce point processes to model random persistence diagrams. Indeed, we incorporate the prior uncertainty by modeling persistence diagrams as Poisson point processes and noisy observations of persistence diagrams as marked Poisson PP to model the level of confidence that observations are representatives of the ground truth. Considering a Poisson point process, one needs to focus on the intensity of the random process. Adapting a prior intensity and a pertinent likelihood, we prove that a posterior intensity can be retrieved. It should be noted that our Bayesian model considers persistence diagrams, which are summaries of the data at hand, for defining a substitution likelihood rather than using the underlying point cloud data. This does not adhere to a strict Bayesian viewpoint, as we model the behavior of the persistence diagrams without considering the underlying data (materials data in our example) used to create it; however, our paradigm incorporates prior knowledge and observed data summaries to create posterior probabilities, analogous to the notion of substitution likelihood detailed in [27]. The general relationship between the likelihood models related to point cloud data and those of their corresponding persistence diagrams remains an important open problem. Furthermore we show that using Gaussian mixture conjugate family of priors. A detailed example is presented to demonstrate posterior intensities for several interesting instances resulted from varying parameters of the model. We establish evidence that our Bayesian framework offers update of prior uncertainty in the light of new evidence in a similar way as the standard Bayesian for random variables. Thus, the Bayesian inference developed here can be reliably used for machine learning and data analysis techniques directly on the space of PDs. Indeed a classification algorithm is derived and successfully applied on materials science data to assess the capability of our Bayesian framework.

Appendix A-Proof of Theorem 3.1

Proof.

By Theorem 2.1, we decompose $\lambda_{\mathcal{D}_{X}|D_{Y^{1:m}}}$ to write

[TABLE]

where the second equality follows because $\mathcal{D}_{X_{V}}$ is independent of $\mathcal{D}_{Y}$. Theorem 2.1 allows us to express $\lambda_{\mathcal{D}_{X_{O}}}$ as the average of intensity functions $\lambda_{\mathcal{D}_{X_{O}^{i}}}$ for $i=1,\cdots,m$, where the $\mathcal{D}_{X_{O}^{i}}$ are independent and equal in distribution to $\mathcal{D}_{X_{O}}$. That is, $\lambda_{\mathcal{D}_{X_{O}}}=\frac{1}{m}\sum_{i=1}^{m}\lambda_{\mathcal{D}_{X_{O}^{i}}}$, and by conditioning we have,

[TABLE]

So to expand Equation (12) it suffices to compute $\lambda_{\mathcal{D}_{X_{O}^{i}}|D_{Y^{i}}}$ for fixed $i$. First, we express the finite PP $(\mathcal{D}_{X},\mathcal{D}_{Y})$ as a marked Poisson PP. To this end, we adopt a construction from [53], the augmented space $\mathbb{W}^{\prime}:=\mathbb{W}\cup\{\Delta\}$, where $\Delta$ is a dummy set used for labeling points in $\mathcal{D}_{Y_{S}}$. Next, we define the random set, $\mathcal{H}=\mathcal{H}_{\mathbb{W}}\cup\mathcal{H}_{\Delta}$ such that

[TABLE]

One can observe that $\mathcal{H}$ is the superposition of two marked Poisson PPs $\mathcal{H}_{\mathbb{W}}$ and $\mathcal{H}_{\Delta}$, taking values in $\mathbb{W}\times\mathbb{W}$ and $\Delta\times\mathbb{W}$, respectively. Moreover, it directly follows from $(M2)$ and $(M3)(i)$ that $\mathcal{H}_{\mathbb{W}}$ has marginal intensity function $\alpha(x)\lambda_{\mathcal{D}_{X}}(x)$ on $\mathbb{W}$ and stochastic kernel density $\ell(y|x)$ while $(M3)(ii)$ shows that $\mathcal{H}_{\Delta}$ has marginal intensity function $\lambda_{D_{Y_{S}}}(\mathbb{W})$ on $\{\Delta\}$ with stochastic kernel density $\frac{\lambda_{D_{Y_{S}}}(y)}{\lambda_{D_{Y_{S}}}(\mathbb{W})}$. By Theorem 2.3, the intensity functions for $\mathcal{H}_{\mathbb{W}}$ and $\mathcal{H}_{\Delta}$ are $\alpha(x)\lambda_{\mathcal{D}_{X}}(x)\ell(y|x)$ and $\lambda_{D_{Y_{S}}}(y)$, respectively. Hence, applying Theorem 2.1 to Equation (14) reveals that the intensity function for $\mathcal{H}$, $\lambda_{\mathcal{H}}$, is given by

[TABLE]

Let $\mathcal{H}_{Y}:=\{y:(x,y)\in\mathcal{H}\}$, $\mathcal{H}_{X}:=\{x:(x,y)\in\mathcal{H}\}$ be the projections of $\mathcal{H}$ onto its first and second coordinates, respectively. It immediately follows from Theorem 2.2 that $\mathcal{H}_{Y}$ is a Poisson PP on $\mathbb{W}$ since it is the image of $\mathcal{H}$ under a projection. Therefore, by treating the first coordinates of $\mathcal{H}$ as marks, we may express $\mathcal{H}$ as a marked Poisson PP having intensity function $\lambda_{\mathcal{H}_{Y}}$ on $\mathbb{W}$ and stochastic kernel density $p(x|y)$ from $\mathbb{W}$ to $\mathbb{W}^{\prime}$. Another application of Theorem 2.3 then implies

[TABLE]

From Equations (15) and (16), we obtain the identity

[TABLE]

Equation (17) describes the probability density of $\mathcal{H}$ at $x\in\mathbb{W}^{\prime}$ for $y\in\mathbb{W}$ fixed. Substituting Equation (17) for the Janossy density in Equation (1) and applying Corollary 2.1 gives the intensity function for the point process $\mathcal{H}_{X}|D_{Y^{i}}$ whenever $\lambda_{\mathcal{H}_{Y}}(y)\neq 0$ for any $y\in D_{Y^{i}}$,

[TABLE]

Restricting Equations (15) and (16) to $\mathbb{W}\times\mathbb{W}$, we obtain

$p(x|y)\lambda_{\mathcal{H}_{Y}}(y)=\ell(y|x)\alpha(x)\lambda_{\mathcal{D}_{X}}(x)$. Thus, $\ell(y|x)\alpha(x)\lambda_{\mathcal{D}_{X}}(x)=0$

whenever $\lambda_{\mathcal{H}_{Y}}(y)=0$, from which we conclude $\lambda_{\mathcal{H}_{Y}}(y)\neq 0$ a.s. . Hence, restricting Equation (18) to $\mathbb{W}\times\mathbb{W}$ yields

[TABLE]

Notice that $\mathcal{H}_{Y}$ is the same PP as $\mathcal{D}_{Y_{O}}\cup\mathcal{D}_{Y_{S}}$. Theorems 2.2 and 2.3 imply that $\mathcal{D}_{Y_{O}}$ is a Poisson PP, and $\mathcal{D}_{Y_{S}}$ is a Poisson PP by (M3), so by Theorem 2.1, $\lambda_{\mathcal{H}_{Y}}=\lambda_{\mathcal{D}_{Y_{O}}}+\lambda_{\mathcal{D}_{Y_{S}}},$ where $\lambda_{\mathcal{D}_{Y_{O}}}(y)=\lambda_{(\mathcal{D}_{X_{O}},\mathcal{D}_{Y_{O}})}(\mathbb{W}\times y)=\int_{\mathbb{W}}\alpha(u)\lambda_{\mathcal{D}_{X_{O}}}(u)\ell(y|u)du$ by Theorem 2.3. Employing Equation (19) one gets that

[TABLE]

which proves Theorem 3.1 after substituting into Equation (12). ∎

Acknowledgments

Research has been partially funded by the Army Research Office, W911NF-17-1-0313, the National Science Foundation, MCB-1715794, and DMS-1821241, and Thor Industries/Army Research Lab, W911NF-17-2-0141.