computePI {TDAvec}

A Vector Summary of the Persistence Surface

Description

For a given persistence diagram D=\{(b_i,p_i)\}_{i=1}^N, computePI() computes the persistence image (PI) - a vector summary of the persistence surface:

\rho(x,y)=\sum_{i=1}^N f(b_i,p_i)\phi_{(b_i,p_i)}(x,y),

where \phi_{(b_i,p_i)}(x,y) is the Gaussian distribution with mean (b_i,p_i) and covariance matrix \sigma^2 I_{2\times 2} and

f(b,p) = w(p)=\left\{ \begin{array}{ll} 0 & \quad y\leq 0 \\ p/p_{max} & \quad 0<p<p_{max}\\ 1& \quad y\geq p_{max} \end{array} \right.

is the weighting function with p_{max} being the maximum persistence value among all persistence diagrams considered in the experiment. Points of D with infinite persistence value are ignored

Usage

computePI(D, homDim, xSeq, ySeq, sigma)

Arguments

D

matrix with three columns containing the dimension, birth and persistence values respectively

homDim

homological dimension (0 for H_0, 1 for H_1, etc.)

xSeq

numeric vector of increasing x (birth) values used for vectorization

ySeq

numeric vector of increasing y (persistence) values used for vectorization

sigma

standard deviation of the Gaussian

Value

A numeric vector whose elements are the average values of the persistence surface computed over each cell of the two-dimensional grid constructred from xSeq=\{x_1,x_2,\ldots,x_n\} and ySeq=\{y_1,y_2,\ldots,y_m\}:

\Big(\frac{1}{\Delta x_1\Delta y_1}\int_{[x_1,x_2]\times [y_1,y_2]}\rho(x,y)dA,\ldots,\frac{1}{\Delta x_{n-1}\Delta y_{m-1}}\int_{[x_{n-1},x_n]\times [y_{m-1},y_m]}\rho(x,y)dA\Big),

where dA=dxdy, \Delta x_k=x_{k+1}-x_k and \Delta y_j=y_{j+1}-y_j

Author(s)

Umar Islambekov

References

1. Adams, H., Emerson, T., Kirby, M., Neville, R., Peterson, C., Shipman, P., ... & Ziegelmeier, L. (2017). Persistence images: A stable vector representation of persistent homology. Journal of Machine Learning Research, 18.

Examples

N <- 100 
set.seed(123)
# sample N points uniformly from unit circle and add Gaussian noise
X <- TDA::circleUnif(N,r=1) + rnorm(2*N,mean = 0,sd = 0.2)

# compute a persistence diagram using the Rips filtration built on top of X
D <- TDA::ripsDiag(X,maxdimension = 1,maxscale = 2)$diagram 

# switch from the birth-death to the birth-persistence coordinates
D[,3] <- D[,3] - D[,2] 
colnames(D)[3] <- "Persistence"

resB <- 5 # resolution (or grid size) along the birth axis
resP <- 5 # resolution (or grid size) along the persistence axis 

# compute PI for homological dimension H_0
minPH0 <- min(D[D[,1]==0,3]); maxPH0 <- max(D[D[,1]==0,3])
ySeqH0 <- seq(minPH0,maxPH0,length.out=resP+1)
sigma <- 0.5*(maxPH0-minPH0)/resP 
computePI(D,homDim=0,xSeq=NA,ySeqH0,sigma) 

# compute PI for homological dimension H_1
minBH1 <- min(D[D[,1]==1,2]); maxBH1 <- max(D[D[,1]==1,2])
minPH1 <- min(D[D[,1]==1,3]); maxPH1 <- max(D[D[,1]==1,3])
xSeqH1 <- seq(minBH1,maxBH1,length.out=resB+1)
ySeqH1 <- seq(minPH1,maxPH1,length.out=resP+1)
sigma <- 0.5*(maxPH1-minPH1)/resP
computePI(D,homDim=1,xSeqH1,ySeqH1,sigma)