Distance Between Persistence Diagrams

As the persistence diagram is a descriptor of the underlying shape structure of a dataset, it can be useful to quantify the differences between pairs of persistence diagrams. There are several ways to compute distances between persistence diagrams in the same homological dimension (like dimension 0 for clusters, dimension 1 for loops, etc.). The most common two are called the 2-wasserstein and bottleneck distances (Kerber, Morozov, and Nigmetov 2017). These techniques find an optimal matching of the 2D points in their input two diagrams, and compute a cost of that optimal matching. A point from one diagram is allowed either to be paired (matched) with a point in the other diagram or its diagonal projection, i.e. the nearest point on the diagonal line $y=x$ (matching a point to its diagonal projection is essentially saying that feature is likely topological noise because it died very soon after it was born).

Allowing points to be paired with their diagonal projections both allows for matchings of persistence diagrams with different numbers of points (which is almost always the case in practice) and also formalizes the idea that some points in a persistence diagram represent noise. The “cost” value associated with a matching is given by either (i) the maximum of infinity-norm distances between paired points, or (ii) the square-root of the sum of squared infinity-norm between matched points. The cost of the optimal matching under loss (i) is called the bottleneck distance of persistence diagrams, and the cost of the optimal matching of cost (ii) is called the 2-wasserstein metric of persistence diagrams. Both distance metrics have been used in a number of applications, but the 2-wasserstein metric is able to find more fine-scale differences in persistence diagrams compared to the bottleneck distance. The problem of finding an optimal matching can be solved with the Hungarian algorithm, which is implemented in the R package clue (Hornik 2005).

In the picture we can see that there is a “better” matching between D1 and D2 compared to D1 and D3, so the (wasserstein/bottleneck) distance value between D1 and D2 would be smaller than that of D1 and D3.

The wasserstein and bottleneck distances have been implemented in the TDApplied function diagram_distance. We can confirm that the distance between D1 and D2 is smaller than D1 and D3 for both distances:

# calculate 2-wasserstein distance between D1 and D2
diagram_distance(D1,D2,dim = 0,p = 2,distance = "wasserstein")
#> [1] 0.3905125

# calculate 2-wasserstein distance between D1 and D3
diagram_distance(D1,D3,dim = 0,p = 2,distance = "wasserstein")
#> [1] 0.559017

# calculate bottleneck distance between D1 and D2
diagram_distance(D1,D2,dim = 0,p = Inf,distance = "wasserstein")
#> [1] 0.3

# calculate bottleneck distance between D1 and D3
diagram_distance(D1,D3,dim = 0,p = Inf,distance = "wasserstein")
#> [1] 0.5

There is a generalization of the 2-wasserstein distance for any $p \geq 1$, the $p$-wasserstein distance, which can also be computed using the diagram_distance function by varying the parameter p. To see the computational details and proof of correctness for wasserstein and bottleneck distance calculations, see the appendix.

Another distance metric between persistence diagrams, which will be useful for kernel calculations, is called the Fisher information metric, $d_{FIM}(D_1,D_2,\sigma)$ (details can be found in (Le and Yamada 2018)). The idea is to represent the two persistence diagrams as probability density functions, with a 2D-Gaussian point mass centered at each point in both diagrams (including the diagonal projections of the points in the opposite diagram), all of variance $\sigma^2 > 0$, and calculate how much those distributions agree on their pdf value at each point in the plane (called their Fisher information metric).

Points in the rightmost plot which are close to white in color have the most similar pdf values in the two distributions, and would not contribute to a large distance value, however having more points with a red color would contribute to a larger distance value.

The diagram_distance function can also calculate the Fisher information metric between persistence diagrams:

# Fisher information metric calculation between D1 and D2 for sigma = 1
diagram_distance(D1,D2,dim = 0,distance = "fisher",sigma = 1)
#> [1] 0.02354779

# Fisher information metric calculation between D1 and D3 for sigma = 1
diagram_distance(D1,D3,dim = 0,distance = "fisher",sigma = 1)
#> [1] 0.08821907

Again, D1 and D2 are less different than D1 and D3 using the Fisher information metric.

Multidimensional Scaling of Persistence Diagrams

Dimension reduction is a task in machine learning which is commonly used for data visualization, removing noise in data, and decreasing the number of covariates in a model (which can be helpful in reducing overfitting). One common dimension reduction technique in machine learning is called multidimensional scaling (MDS) (Cox and Cox 2008). MDS takes as input an $n$ by $n$ distance (or dissimilarity) matrix $D$, computed from $n$ points in a dataset, and outputs an embedding of those points into a Euclidean space of chosen dimension $k$ which best preserves the inter-point distances. MDS is often used for visualizing data in exploratory analyses, and can be particularly useful when the input data points do not live in a shared Euclidean space (as is the case for persistence diagrams). Using the R function cmdscale from the package stats (R Core Team (2021)) we can compute the optimal embedding of a set of persistence diagrams using any of the three distance metrics using the function diagram_mds. Here is an example of the diagram_mds function projecting nine persistence diagrams, three noisy copies sampled from each of D1, D2 and D3, into 2D space:

# create 9 diagrams based on D1, D2 and D3
g <- generate_TDApplied_vignette_data(3,3,3)

# calculate their 2D MDS embedding in dimension 0 with the bottleneck distance
mds <- diagram_mds(diagrams = g,dim = 0,p = Inf,k = 2,num_workers = 2)

# plot
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(mds[,1],mds[,2],xlab = "Embedding coordinate 1",ylab = "Embedding coordinate 2",
     main = "MDS plot",col = as.factor(rep(c("D1","D2","D3"),each = 3)),bty = "L")
legend("topright", inset=c(-0.2,0), 
       legend=levels(as.factor(c("D1","D2","D3"))), 
       pch=16, col=unique(as.factor(c("D1","D2","D3"))))

The MDS plot shows the clear separation between the three generating diagrams (D1, D2 and D3), and the embedded coordinates could be used for further downstream analyses.

Testing for Group Differences of Persistence Diagrams

One of the most important inference procedures in classical statistics is the analysis of variance (ANOVA), which can find differences in the means of groups of normally-distributed measurements (Casella, Berger, and Company 2002). Distributions of persistence diagrams and their means can be complicated (see (Turner 2013) and (Turner et al. 2014)). Therefore, a non-parametric permutation test has been proposed which can find differences in groups of persistence diagrams. Such a test was first proposed in (Robinson and Turner 2017), and some variations have been suggested in later publications. In (Robinson and Turner 2017), two groups of persistence diagrams would be compared. The null hypothesis, $H_0$, is that the diagrams from the two groups are generated from shapes with the same type and scale of topological features, i.e. they “come” from the same “shape.” The alternative hypothesis, $H_A$, is that the underlying type or scale of the features are different between the two groups. In each dimension a p-value is computed, finding evidence against $H_0$ in that dimension. A measure of within-group distances (a “loss function”) is calculated for the two groups, and that measure is compared to a null distribution for when the group labels are permuted.

This inference procedure is implemented in the permutation_test function, with several speedups and additional functionalities. Firstly, the loss function is computed in parallel for scalability since distance computations can be expensive. Secondly, we store distance calculations as we compute them because these calculations are often repeated. Additional functionality includes allowing for any number groups (not just two) and allowing for a pairing between groups of the same size as described in (Abdallah (2021)). When a natural pairing exists between the groups (like if the groups represent persistence diagrams from the same subject of a study in different conditions) we can simulate a more realistic null distribution by restricting the way in which we permute group labels, achieving higher statistical power.

In order to demonstrate the utility of the permutation test we will detect differences between noisy copies of D1, D2, D3:

# permutation test between three diagrams
g1 <- generate_TDApplied_vignette_data(3,0,0)
g2 <- generate_TDApplied_vignette_data(0,3,0)
g3 <- generate_TDApplied_vignette_data(0,0,3)
perm_test <- permutation_test(g1,g2,g3,
                              num_workers = 2,
                              dims = c(0))
perm_test$p_values
#>          0 
#> 0.04761905

As expected, a difference was found (at the $\alpha = 0.05$ significance level) between the three groups.

The package TDAstats also has a function called permutation_test which is based on the same test procedure, however it uses an unpublished distance metric between persistence diagrams and does not use parallelization for scalability. As such, care must be taken if both TDApplied and TDAstats are attached in an R script to use the desired permutation_test function.

Finding Real Topological Features via Bootstrapping

While the focus of TDApplied is applied topological data analysis, and as such this vignette focuses on machinery for carrying out such analyses, a fundamental question in topological data analysis which must not be overlooked is how to tell if a point in a persistence diagram should be considered “real” (signal) or not (noise). This question has been addressed via a bootstrapping approach in (B. Fasy et al. 2014), which has been implemented in the TDA package. Due to the speed increase of TDApplied’s distance function (discussed in a later section), TDApplied has its own implementation of this bootstrap procedure.

The idea of the bootstrap procedure is as follows, where $X$ is the input data set and $\alpha$ is the desired threshold for type 1 error:

1. We first compute the diagram of $X$, $D$.
2. Then we repeatedly sample, with replacement, the original data to obtain $\{X_1,\dots,X_n\}$and compute new persistence diagrams $\{D_1,\dots,D_n\}$.
3. We calculate the bottleneck distance of each new diagram with the original, $d_{\infty}(D_i,D)$, in each dimension separately.
4. Finally we compute the $1-\alpha$ percentile of these values in each dimension.

These thresholds form a square-shaped “confidence interval” around each point in $D$. In particular, if $t$ was the threshold found for dimension $k$ then the confidence interval around a point $(x,y) \in D$ (of dimension $k$) is the set of points $\{(x',y'):\mbox{max}(|x-x'|,|y-y'|) < t\}$. For example, if we sampled 50 points on the unit circle, calculated the bottleneck-threshold-based confidence interval around the single 1-dimensional point, we would get something like this:

In this setup, “real” points will be those whose (open) confidence intervals do not overlap the diagonal line, where birth and death is the same. Note that the persistence threshold is twice the bottleneck distance threshold calculated above: the (Euclidean) distance from the point to the bottom right corner of the box is $\sqrt{2}t$, which must be greater than or equal to the (Euclidean) distance of the point to its diagonal projection, which is $\frac{y-x}{\sqrt{2}}$. Therefore, for the point to be considered real, $\sqrt{2}t \leq \frac{y-x}{\sqrt{2}}$, implying that the persistence of the point, $y-x$, must be no less than twice the bottleneck threshold $t$.

The function bootstrap_persistence_thresholds can be used to determine these persistence thresholds. Here is an example for a circle dataset (which is not the same as the one from above), and the results can be plotted with the plot_diagram function:

# sample 50 points from the unit circle
circle <- TDA::circleUnif(n = 50)

# calculate the bootstrapped persistence thresholds using 2 cores
# and 20 iterations
thresh <- bootstrap_persistence_thresholds(X = circle,FUN = "ripsDiag",
                                 maxdim = 1,thresh = 2,num_workers = 2,
                                 num_samples = 30)
diag <- thresh$diag

# plot original diagram and thresholded diagram side-by-side:
par(mfrow = c(1,2))

plot_diagram(diag,title = "Circle diagram")

plot_diagram(diag,title = "Circle diagram with thresholds",
             thresholds = thresh$thresholds)

The bootstrap procedure can be done in parallel or sequentially, depending on which function is specified to calculate the persistence diagrams. There are three possible such functions - TDAstats’ calculate_homology, TDA’s ripsDiag and TDApplied’s PyH (described in later sections). Based on our simulations, the order of fastest to slowest of these functions are PyH, calculate_homology and ripsDiag. Moreover, both ripsDiag and PyH allow for the calculation of representative (co)cycles (i.e. the approximate location in the data of each topological feature), whereas calculate_homology does not. Therefore, our recommendation is to pick the function according to the following criteria: if a user can use the PyH function, then it should be used in all cases except for when the input data set is small, the machine has many available cores and the number of desired bootstrap iterations is large. Otherwise, use calculate_homology for speed or ripsDiag if the other two functions are not available.

Kernels Between Persistence Diagrams

A kernel function is a special (positive semi-definite) symmetric similarity measure between objects in some complicated space which can be used to project data into a space suitable for machine learning (Murphy 2012). Some examples of machine learning techniques which can be “kernelized” when dealing with complicated data are k-means (kernel k-means), principal components analysis (kernel PCA), and support vector machines (SVM) which are inherently based on kernel calculations.

There have been, to date, four main kernels proposed for persistence diagrams. In TDApplied the persistence Fisher kernel (Le and Yamada 2018) has been implemented because of its practical advantages over the other kernels – smaller cross-validation SVM error on a number of test data sets and a faster method for cross validation. For information on the other three kernels see (Kusano, Fukumizu, and Hiraoka 2018), (Carrière, Cuturi, and Oudot 2017), and (Reininghaus et al. 2014).

The persistence Fisher kernel is computed directly from the Fisher information metric between two persistence diagrams: let $\sigma > 0$ be the parameter for $d_{FIM}$, and let $t > 0$. Then the persistence Fisher kernel is defined as $k_{PF}(D_1,D_2) = \mbox{exp}(-t*d_{FIM}(D_1,D_2,\sigma))$. Computing the persistence Fisher kernel can be achieved with the diagram_kernel function in TDApplied:

# calculate the kernel value between D1 and D2 with sigma = 2, t = 2
diagram_kernel(D1,D2,dim = 0,sigma = 2,t = 2)
#> [1] 0.9872455
# calculate the kernel value between D1 and D3 with sigma = 2, t = 2
diagram_kernel(D1,D3,dim = 0,sigma = 2,t = 2)
#> [1] 0.9707209

As before, D1 and D2 are more similar than D1 and D3.

Kernel k-means of Persistence Diagrams

Kernel k-means (Dhillon, Guan, and Kulis 2004) is a method which can find hidden groups in complex data, extending regular k-means clustering (Murphy 2012) via a kernel. A “kernel distance” is calculated between a persistence diagram and a cluster center using only the kernel function, and the algorithm converges like regular k-means. This algorithm is implemented in the function diagram_kkmeans as a wrapper of the kernlab function kkmeans. Moreover, a prediction function predict_diagram_kkmeans can be used to find the nearest cluster labels for a new set of diagrams. Here is an example clustering three groups of noisy copies from D1, D2 and D3:

# create noisy copies of D1, D2 and D3
g <- generate_TDApplied_vignette_data(3,3,3)
                              
# calculate kmeans clusters with centers = 3, and sigma = t = 2
clust <- diagram_kkmeans(diagrams = g,centers = 3,dim = 0,t = 2,sigma = 2,num_workers = 2)

# display cluster labels
clust$clustering@.Data
#> [1] 3 3 3 1 1 1 2 2 2

As we can see, the diagram_kkmeans function was able to correctly separate the three generating diagrams D1, D2 and D3 (the cluster labels are arbitrary and therefore may not be 1,1,1,2,2,2,3,3,3, however the induce they correct partition).

If we wish to predict the cluster label for new persistence diagrams (computed via the largest kernel value to any cluster center), we can use the predict_diagram_kkmeans function as follows:

# create nine new diagrams
g_new <- generate_TDApplied_vignette_data(3,3,3)

# predict cluster labels
predict_diagram_kkmeans(new_diagrams = g_new,clustering = clust,num_workers = 2)
#> [1] 3 3 3 1 1 1 2 2 2

This function correctly predicted the cluster for each new diagram (assigning each diagram to the cluster label by D1, D2 or D3, depending on which diagram it was generated from).

Kernel Principal Components Analysis of Persistence Diagrams

PCA is another dimension reduction technique in machine learning, but can be preferable compared to MDS in certain situations because it allows for the projection of new data points onto an old embedding model (Murphy 2012). For example, this can be important if PCA is used as a pre-processing step in model fitting. Kernel PCA (kPCA) (Schölkopf, Smola, and Müller 1998) is an extension of regular PCA which uses a kernel to project complex data into a high-dimensional Euclidean space and then uses PCA to project that data into a low-dimensional space. The diagram_kpca method computes the kPCA embedding of a set of persistence diagrams, and the predict_diagram_kpca function can be used to project new diagrams using a pre-trained kPCA model. Here is an example using a group of noisy copies of D1, D2 and D3:

# create noisy copies of D1, D2 and D3
g <- generate_TDApplied_vignette_data(3,3,3)

# calculate their 2D PCA embedding with sigma = t = 2
pca <- diagram_kpca(diagrams = g,dim = 0,t = 2,sigma = 2,features = 2,num_workers = 2)

# plot
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(pca$pca@rotated[,1],pca$pca@rotated[,2],xlab = "Embedding coordinate 1",
     ylab = "Embedding coordinate 2",main = "PCA plot",
     col = as.factor(rep(c("D1","D2","D3"),each = 3)))
legend("topright",inset = c(-0.2,0), 
       legend=levels(as.factor(c("D1","D2","D3"))), pch=16, 
       col=unique(as.factor(c("D1","D2","D3"))))

The function was able to recognize the three groups, and the embedding coordinates can be used for further downstream analysis. However, an important advantage of kPCA over MDS is that in kPCA we can project new points onto an old embedding using the predict_diagram_kpca function:

# create nine new diagrams
g_new <- generate_TDApplied_vignette_data(3,3,3)

# project new diagrams onto old model
new_pca <- predict_diagram_kpca(new_diagrams = g_new,embedding = pca,num_workers = 2)

# plot
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(new_pca[,1],new_pca[,2],xlab = "Embedding coordinate 1",
     ylab = "Embedding coordinate 2",main = "PCA prediction plot",
     col = as.factor(rep(c("D1","D2","D3"),each = 3)))
legend("topright",inset = c(-0.2,0), 
       legend=levels(as.factor(c("D1","D2","D3"))), pch=16, 
       col=unique(as.factor(c("D1","D2","D3"))))

As we can see, the original three groups, and their approximate location in 2D space, is preserved during prediction.

Kernel Support Vector Machines of Persistence Diagrams

SVMs (Murphy 2012) are one of the most popular machine learning techniques for regression and classification tasks. SVMs use a kernel function to project complex data into a high-dimensional space and then find a sparse set of training examples, called “support vectors,” which maximally linearly separate the outcome variable classes (or yield the highest explained variance in the case of regression).

SVMs have been implemented in the function diagram_ksvm, tailored for input datasets which contain pairs of persistence diagrams and their outcome variable labels. A prediction method is supplied called predict_diagram_ksvm which can be used to predict the label value of a set of new persistence diagrams given a pre-trained model. A parallelized implementation of cross-validation model-fitting is used based on the remarks in (Le and Yamada 2018) for scalability (which avoids needlessly recomputing persistence Fisher information metric values). Here is an example of fitting an SVM model on a list of persistence diagrams for a classification task (guessing whether the diagram comes from D1, D2 or D3):

# create thirty noisy copies of D1, D2 and D3
g <- generate_TDApplied_vignette_data(10,10,10)

# create response vector
y <- as.factor(rep(c("D1","D2","D3"),each = 10))

# fit model with cross validation
model_svm <- diagram_ksvm(diagrams = g,cv = 2,dim = c(0),
                          y = y,sigma = c(1,0.1),t = c(1,2),
                          num_workers = 2)

We can use the function predict_diagram_ksvm to predict new diagrams like so:

# create nine new diagrams
g_new <- generate_TDApplied_vignette_data(3,3,3)

# predict
predict_diagram_ksvm(new_diagrams = g_new,model = model_svm,num_workers = 2)
#> [1] D1 D1 D1 D2 D2 D2 D3 D3 D3
#> Levels: D1 D2 D3

As we can see the best SVM model was able to separate the three diagrams We can gain more information about the best model found during model fitting and the CV results by accessing different list elements of model_svm.

Testing for Independence Between Two Groups of Paired Persistence Diagrams

An important question when presented with two groups of paired persistence diagrams is determining if the pairings are independent or not. A procedure was described in (Gretton et al. 2007) which can be used to answer this question using kernel computations, and importantly uses a parametric null distribution. The null hypothesis for this test is that the groups are independent, and the alternative hypothesis is that the groups are not independent. A test statistic called the Hilbert-Schmidt independence criteria is calculated, and its value is compared to a gamma distribution with certain parameters which can be estimated from the data.

This inference procedure has been implemented in the independence_test function, and returns the p-value of the test in each desired dimension of the diagrams (among other additional information). We would expect to find no dependence between noisy copies of D1, D2 and D3, since each copy is generated randomly:

# create 10 noisy copies of D1 and D2
g1 <- generate_TDApplied_vignette_data(10,0,0)
g2 <- generate_TDApplied_vignette_data(0,10,0)

# do independence test with sigma = t = 1
indep_test <- independence_test(g1,g2,dims = c(0),num_workers = 2)
indep_test$p_values
#>         0 
#> 0.1354271

The p-value of this test would not be significant at any typical significance threshold, reflecting the fact that there is no real (i.e. non-spurious) dependence between the two groups, as expected.

PyH Requirements and Usage

There are three prerequisites that must be satisfied in order to use the PyH function:

1. The reticulate package must be installed.
2. Python must be installed and configured to work with reticulate.
3. The ripser python module must be installed, via reticulate::py_install("ripser").

Note that the installation status of python is checked using the function reticulate::py_available(), which according to several online forums does not always behave as expected. If error messages occur using TDApplied functions regarding python not being installed then we recommend consulting online resources to ensure that the py_available function returns TRUE on your system. Due to the complicated dependencies required to use the PyH function, it is only an optional function in the package and therefore the reticulate package is only suggested in the TDApplied namespace.

For an example use of PyH we will use the following pre-loaded dataset called “circ” (which is stored as a data frame in this vignette):

We would then calculate the persistence diagram as follows:

# import the ripser module
ripser <- import_ripser()

# calculate the persistence diagram
diag <- PyH(X = circ,maxdim = 1,thresh = 1,ripser = ripser)

The ripser module is imported outside of the PyH function using the import_ripser function, because this operation can be slow and if many persistence diagrams are to be computed then repeating this operation would be highly prohibitive. Additionally, the import_ripser function checks the three prerequisites for PyH listed earlier in this section, to avoid repeating checks. The import_ripser function is the same as reticulate::import("ripser"), but with the extra checks. In the TDApplied testing it was verified that the PyH output matches the TDAstats calculate_homology output for the same input data and parameters, but with one small difference: in the python ripser calculation there is always one connected component (i.e. dimension 0 topological feature) with infinite death value, whereas this feature is omitted in the TDAstats calculation. Whether or not to include this feature in the outputted diagram can be decided by the user via the ignore_infinite_cluster parameter being set to TRUE (the default value) or FALSE respectively.

Representative Cocycles

One of the advantages of the R package TDA over TDAstats is its ability to calculate representative cycles in the data, i.e. localizing the persistence diagram topological features in the input data set. This can permit deep analyses of the original data set by finding particular types of features spanned by certain subsets of data points. For example, a representative cycle of a 1-dimensional topological feature would be a set of edges between data points which lie along that feature (a loop). The PyH function can also find representative (co)cycles in its input data, which are returned if the calculate_representatives parameter is set to TRUE. In that case, the PyH function returns a list with a data frame called “diagram,” containing the persistence diagram, and a list called “representatives.” The representatives element has one element for each dimension in the persistence diagram, with one matrix/array for each point in the persistence diagram of that dimension (except for dimension 0, where the list is always empty). The matrix for a $d$-dimensional feature (1 for loops, 2 for voids, etc.) has $d + 1$ columns, where row $i$ contains the row indices in the data set of the data points in the $i$-th substructure in the representative (a substructure of a loop would be an edge, a substructure of a void would be a triangle, etc.). Here is an example where we calculate the representative cocycles of our circ dataset:

# ripser has already been imported, so calculate diagram with representatives
diag <- PyH(circ,maxdim = 1,thresh = 2,ripser = ripser,calculate_representatives = T)

# identify the loops in the diagram
diag$diagram[which(diag$diagram$dimension == 1),]

#>    dimension     birth    death
#> 50         1 0.6008424 1.738894

# show the representative for the loop, just the first five rows
diag$representatives[[2]][[1]][1:5,]

#>      [,1] [,2]
#> [1,]   10    7
#> [2,]   36    7
#> [3,]   10    4
#> [4,]   17    7
#> [5,]   32    7

The representative of the one loop contains the edges found to be present in the loop. We could iterate over the representative for the loop to find all the data point in that representative:

unique(c(diag$representatives[[2]][[1]][,1],diag$representatives[[2]][[1]][,2]))

#>  [1] 10 36 17 32 44 30 29 40 14 27 39 28 11 31  9 25  7 41  4 16  8 49 21 34  3
#> [26]  2

Since the circ dataset is two-dimensional, we can actually plot the loop representative according to the following process:

1. Pick a threshold $\epsilon$ less than the death radius of the cocycle.
2. Plot all edges between pairs of points in circ of distance no more than $\epsilon$.
3. Highlight all edges in the representative.

Since the death radius of the main cocycle is 1.738894, we can use the following code to plot the cocycle at thresholds value 1.7:

plot(x = circ$x,y = circ$y,xlab = "x",ylab = "y",main = "circ with representative")
for(i in 1:nrow(circ))
{
  for(j in 1:nrow(circ))
  {
    pt1 <- circ[i,]
    pt2 <- circ[j,]
    if(sqrt((pt1[[1]] - pt2[[1]])^2 + (pt1[[2]] - pt2[[2]])^2) <= 1.7)
    {
      graphics::lines(x = c(pt1[[1]],pt2[[1]]),y = c(pt1[[2]],pt2[[2]]))
    }
  }
}

for(i in 1:nrow(diag$representatives[[2]][[1]]))
{
  pt1 <- circ[diag$representatives[[2]][[1]][i,1],]
  pt2 <- circ[diag$representatives[[2]][[1]][i,2],]
  if(sqrt((pt1[[1]] - pt2[[1]])^2 + (pt1[[2]] - pt2[[2]])^2) <= 1.7)
  {
    graphics::lines(x = c(pt1[[1]],pt2[[1]]),y = c(pt1[[2]],pt2[[2]]),col = "red")
  }
}

This plot shows the main loop that was found via persistent cohomology, and the representative is a set of edges (in this 2D case) whose removal would destroy the loop. A more intuitive notion of a “representative loop” can be found with persistent homology, for instance using the TDA ripsDiag function with the cycleLocation parameter set to TRUE. Another example of the representative cocycle can be found using another threshold value, for instance the (rounded up) birth value of 0.6009:

Since we know that the data points in the representatives help form a loop in the original data set, we could perform a further exploratory analysis in circ to explore the periodic nature of the feature. While in this example we know that a loop will be present, this type of analysis could help find hidden latent structure in data sets. In this way PyH is as flexible as the TDA persistent homology function ripsDiag, although representative cocycles are less intuitive than representative cycles.

Benchmarking PyH Against TDAstats’ calculate_homology Function

The long calculation time of persistence diagrams is likely a large contributing factor to the slow adoption of topological data analysis for applied data science. Much research has been carried out in order to speed up these calculations, but the current state-of-the-art is the persistent cohomology algorithm (Silva and Vejdemo-Johansson 2011a). In R, the TDAstats’ calculate_homology function is the fastest option for persistence diagram calculations (Somasundaram et al. 2021), being a wrapper for the ripser persistent cohomology engine. TDApplied’s PyH function is a wrapper for the same engine, so we benchmarked their run time on circles, spheres and tori. The results were as follows:

As we can see the run time of PyH was significantly less than that of calculate_homology. We then used a linear model to analyze how the two functions scale compared to each other for all three shapes, regressing the ratio of TDAstats’ runtime divided by TDApplied’s runtime onto the number of points in the data set. For each of the three datasets, the model intercept was highly significant ($p < 0.001$) and the coefficient estimate for number of points was not. This suggests that the runtime of the two functions scale similarly, but there was a constant speed increase of PyH which differed by shape (about 3 times faster for circles, 7 times faster for tori and 8 times faster for spheres). Overall, if python is available to a TDApplied user then the PyH function can provide a very fast means to calculate persistence diagrams in R making analyses containing many large persistent homology calculations much more feasible.

Benchmarking the TDApplied diagram_distance and TDA wasserstein Functions

Computing wasserstein (or bottleneck) distances between persistence diagrams is a key feature of some of the main topological data analysis software packages in R and python. However, these calculations can be very expensive, rendering practical applications of topological data analysis nearly unfeasible. Since TDAstats has implemented an unconventional distance calculation, we will benchmark TDApplied’s diagram_distance function against the TDA wasserstein function on spheres and tori, calculating their distance in dimensions 0, 1 and 2 and recording the total time. The results were as follows:

A linear model, regressing the ratio of TDA’s runtime divided by TDApplied’s runtime onto the number of points in the data set, found a significant positive coefficient of the number of points. This suggests that diagram_distance scales better than wasserstein, and the model estimated a 95x speed up for 1000 data points. These results suggest that distance calculations with TDApplied are faster and more scalable, making the applications of statistics and machine learning with persistence diagrams more feasible in R. This is why the TDA distance calculation was not used in the TDApplied package.

Benchmarking the TDApplied diagram_distance Function against persim’s wasserstein Function

While the functionality of python packages for topological data analysis packages are out of the scope for an R package, in order to fully situate TDApplied in the landscape of topological data analysis software we will benchmark the diagram_distance function against its counterpart from the sci-kit TDA collection of libraries, namely the wasserstein function from the persim python module. The R package reticulate (Ushey, Allaire, and Tang 2022) was used to carry out this benchmarking, via installing, importing and using the persim module. This benchmarking procedure also used spheres and tori, calculating distances in dimensions 0, 1 and 2, and the results were as follows:

The runtime of the persim wasserstein function was significantly faster than TDApplied’s diagram_distance function. However, a linear model of the runtime ratio of TDApplied vs. persim against the number of points in the shape finds evidence that the two functions scale similarly, since the estimated coefficient for number of points was not significant but the intercept (15) was highly significant. Nevertheless, the raw speed increase in python could be the basis for a very fast python counterpart to the TDApplied package in the future.

Visualizing mental state topologies

What does a mental state topology, or a collection of them, represent? Loops in a mental state topology are sets of spatially periodic patterns of spatial activity, and these sets vary along one neurological dimension of activity. Collections of mental state topologies, like combining all of the persistence diagrams from the resting state 1 task across all 100 subjects, could encode both a characteristic number and representative neurological dimensions of loops across the whole collection. In this section we visualize combined diagrams across subjects and tasks, and construct representative loops from all the resting state 1 mental state topologies.

We visualized the persistence diagrams of subjects and tasks by combining all the diagrams for that subjects/task and plotting a probability distribution induced by adding a Gaussian point mass centered at each point (like in the Fisher information metric) with $\sigma$ values 0.001, 0.01 and 0.1. Plots with $\sigma = 0.1$ (and $\sigma = 0.01$ to a lesser degree) were least informative because $\sigma$ was comparable to the scale of the data. Here are the distribution plots for three subjects with $\sigma = 0.001$:

As we can see the distributions were very different between these three subjects. Next we visualized the nine tasks (again with $\sigma = 0.001$):

Again we found that the tasks had visually different distributions, although the two resting state scans are fairly similar, as we would expect, with two nearby clusters of loops.

But what do these mental state loops represent in terms of brain activity? As an illustration we considered the resting state 1 persistence diagram and used kmeans clustering to find two “basis” resting state 1 loops which described all of the loops well:

These two loops were loops from two subject’s resting state 1 data, but note that topologically similar loops (in terms of birth and death values) may lie in different parts of the data (i.e. represent different loops of spatial activity patterns). We calculated the representative (co)cycles of these two loops, used cmdscale to project the representatives into 2D (needing two dimensions for the first loop and four dimensions for the second) and selected a top, bottom, left and right points along the loops (in red):

Since the highlighted points from the plots each represent a spatial pattern of neural activity at some time point, we then plotted these spatial patterns around the circumference of two circles (however, note that the loops did not form perfect circles in the data). Note that when the two brain hemispheres (i.e. sides) are plotted side-by-side, the left hemisphere is on the right and the right hemisphere is on the left. Prior to plotting, the fMRI data was first normalized to have temporal mean 0 and standard deviation 1 for each surface node. The spatial patterns were then smoothed by averaging spatial signal within a 5mm radius from each surface node, and the smoothed signal was thresholded by absolute value at least 0.75. We plotted the spatial patterns with boxes around areas which appeared to exhibit periodic behavior around the loops:

Loop 1 exhibited more drastic neural activity (both high and low) compared to loop 2, which showed more localized areas of low activity. We compared the locations of periodic activity (in the boxes) to a well-known parcellation of the brain into 7 resting-state networks (Yeo et al. 2011):

The periodic behavior in loop 1 appeared in the activation of the left hemisphere dorsal attention network (green box, high activity on top and bottom, lower activity on the sides) and in the deactivation of the right hemisphere default mode and ventral attention networks (red box, high activity on the right, medium on the top and bottom, low on the left). The periodic behavior in loop 2 appeared in the deactivation of the right hemisphere default mode and ventral attention networks (pink box, low activity on top and bottom, higher activity on the sides). Therefore only three of the seven resting-state networks may be sufficient to describe the group-level resting-state data.

Coming up with techniques for determining what the one circular dimension of one such loop represents in terms of neural activity, and determining if this dimension is correlated with time/study design (i.e. the temporal sequence of events in the task) for resting-state/task-based fMRI is left to future work.

Comparing subjects and tasks with inference

Some major functions of healthy human brains are conserved across people, but there is still much inter-subject variability (different people can respond differently to viewing the same image). But are task-induced mental state topologies conserved across subjects or tasks? To answer this question we first analyzed distributional properties of the persistence diagrams between subjects and between tasks, and compared the diagrams between subjects and between tasks using TDApplied’s permutation and independence tests.

We first analyzed basic properties of the persistence diagrams - the number of loops and the mean persistence of loops, both across subjects and tasks. The overall distribution of the number of loops across all diagrams was:

The bootstrap thresholding was evidently very conservative, with about 73% of the diagrams containing no loops. This is not to say that those diagrams were topologically trivial - on the contrary they may have contained non-trivial higher-dimensional topological structure, or it is possible that task-correlated loops may exist at smaller (i.e. less persistent) scales. To assess the validity of either explanation we would require significantly more computational resources than what was available at the time of writing, and we leave this to future work.

We may draw two important conclusions from this graphic. Firstly, the fact that many diagrams did not have any loops provides evidence that topological structure was not driven by physiological sources, such as heart-rate and breathing, which is known to drive signal in fMRI data (Caballero-Gaudes and Reynolds 2017). Secondly, resting-state analyses are often concerned with the activity of (about) seven major “resting-state networks” found in (Yeo et al. 2011), and task analyses often study the activity of a similar numbers of task-specific networks. However, our loop analysis found that fewer loop dimensions (perhaps up to 2 loops) were needed to describe the HCP data, suggesting high-dimensional topological features may be a useful low-dimensional representation of spatial fMRI data.

We then compared the distribution of the number of loops across subjects and tasks using paired t-tests (with two-sided alternative hypotheses and at the $\alpha = 0.05$ threshold for type 1 error). About 19% of the subject-subject pairs had statistically different numbers of loops, and about 22% of the task-task pairs had statistically significant numbers of loops. Interestingly, most of the differences between tasks were found between resting state scans and task scans, with less loops in the resting state scans, suggesting that tasks may induce large-scale and periodic cycles of mental states which are used to accomplish the task.

Next we performed the same analysis for the mean persistence of loops in the diagrams. About 12% of the subject-subject pairs had differences in mean persistence values. We found 14% statistically different persistence between tasks, mainly being lower persistence in resting state scans compared to task-based scans.

In order to better understand subject-subject and task-task differences we used TDApplied’s permutation test to compare subjects (matched by task) and tasks (matched by subjects) with both the 2-wasserstein and bottleneck distances, again judging significance at the $\alpha = 0.05$ level. For the subjects, both bottleneck and wasserstein distances found about 20% of the subject-subject pairs to have different topologies, with about 96% of the differences conserved between the two metrics. For the tasks, in the bottleneck tests there were 4 pairs of different tasks: both resting state scans with the emotion task, the first resting scan with the social task and the emotional task with the relational task. For the wasserstein tests there were 8 pairs of different tasks: the first resting state scan with the emotion, gambling, motor and relational tasks, the second resting state scan with the emotion, motor and working memory tasks, and the relational task with the working memory task. Therefore, the resting state scans were more different from task scans than the task scans were with each other.

Finally, we analyzed whether any of the subject-subject or task-task pairs were independent (again at the $\alpha = 0.05$ significance level). Only about 18% of the subject-subject pairs were significantly non-independent - an interesting future question would be to determine whether these non-independent pairs were blood-related. On the other hand, over 86% of the task-task pairs were found to be significantly non-independent in at least half of the kernel parameter values. Therefore, task-induced mental state topology may be strongly related between tasks - i.e. computational strategies in different neural domains may be related. In order to compare the results of the two inference procedures, we removed kernel parameters which resulted in no non-independent subject-subject pairs, and then we used a two-sided t-test to compare the percentage of those kernel values which found significant non-independence between subject-subject pairs who were significantly different in the permutation tests compared to pairs of subjects who were not significantly different in the tests (for both distance metrics). The results were very similar for both distance metrics, with significantly lower percentage of non-independence in non-different subject-subject pairs.

Machine learning

Our inference analyses revealed some of the underlying topological differences of rest and task-induced neural function in the HCP dataset. However, what is the overall structure of the space of all the diagrams? In this section we will first cluster the diagrams to find distinct mental-state topological “modes,” and then visualize the clusters using 2D MDS computed with both the bottleneck and wasserstein metrics.

Using all of the kernel parameters deduced earlier, we performed kernel kmeans clustering on all the persistence diagrams over a range of values of k - 2,3,4,…,18,25,50,100. These values were chosen because there were 18 scans per subject and 100 subjects. The within-cluster sum of squared (withinss) distances were calculated for each value of k and for each kernel parameter combination, and the value of k which minimized the withinss (without resulting in empty clusters) was selected as being optimal for each kernel parameter. The frequency of each number of optimal clusters was:

From the graph it seemed like the optimal number of clusters was likely either 3 or 4. We therefore performed clustering with both values of k (3 and 4) for all kernel parameters, and calculated a membership percentage matrix for each value of k - for what percentage of kernel parameters was each pair of subjects in the same cluster. When k was 3 we found that about 54% of the pairs were in the same cluster at least 90% of the time, and that 31% of the pairs were in the same cluster less than 5% of the time. Similarly for four clusters, 53% of the pairs were in the same cluster at least 90% of the time, and 34% percent of the pairs were in the same cluster less than 5% of the time. These results indicate that the clusterings for k being 3 and 4 were highly consistent across all kernel parameters. We then used the kernlab function kkmeans to cluster the two membership matrices, and only for k = 3 was the clustering able to be resolved (potentially due to the membership matrix not being a real kernel matrix). We therefore visualized the group of diagrams in each of the three found clusters:

The three clusters are quite distinct in terms of what topological information they likely represent. For cluster one this is mainly two loops but sometimes a small or very large loop (with mean persistence in the cluster of about 0.0018), for cluster two this is one loop of variable size (mean persistence in the cluster of about 0.0011), and cluster 3 is one main loop and sometimes one additional small loop (this group only had 16 loops in it, with mean persistence about 0.0005).

Next we analyzed the number of rows in the diagrams of each cluster. The results were as follows:

#> [1] "Cluster 1:"
#>   1   2   3   4   5   6   7 
#> 112  38  19   6   2   1   1
#> [1] "Cluster 2:"
#>   1   2 
#> 292   7
#> [1] "Cluster 3:"
#>    0    1 
#> 1306   16

This confirmed that cluster 1 had more loops than the other clusters on average, that cluster 2 mainly had one loop, and that cluster three was only a few loops of low persistence (and all of the no loop diagrams).

Next we analyzed the task composition of clusters, and found the following percentages (rounded to the hundredths position and ordered by most to least rounded variance):

#> [1] "Cluster 1:"
#>      rest1      rest2    emotion   gambling     social      motor relational 
#>       0.20       0.17       0.05       0.07       0.07       0.09       0.11 
#>   language 
#>       0.12
#> [1] "Cluster 2:"
#>      rest1      rest2    emotion   gambling     social      motor relational 
#>       0.09       0.08       0.12       0.15       0.14       0.09       0.12 
#>   language 
#>       0.11
#> [1] "Cluster 3:"
#>      rest1      rest2    emotion   gambling     social      motor relational 
#>       0.10       0.11       0.12       0.11       0.11       0.12       0.10 
#>   language 
#>       0.11

Cluster 1 had much higher percentages of resting state scans compared to the other clusters. However, it also had lower percentages of emotion, gambling and social tasks than the other clusters. Clusters 2 and 3 were moreso separated by percentages of resting state 2 and motor (more in cluster 3), and gambling and social (more in cluster 2).

We then calculated 2D MDS embeddings of all the diagrams using both the wasserstein and bottleneck distances. No external factor like task, subject ID, subject gender or subject age clearly separated the embedded points, however the cluster labels did a reasonably good job:

Behavioral measures may correlate with the mental state topologies, and therefore may be explaining some of the variance within the plots (like the 0-dimensional topological features in (Anderson et al. 2018)), and this could be an interesting future analysis. Also cluster 1 seemed to be more 2 dimensional in both embeddings whereas cluster 2 seemed to be more one dimensional, and determining what factors could explain these dimensions would also be an interesting future analysis.

Conclusions of HCP Analysis

We performed an exploratory analysis of the loops in HCP mental state topologies, and found interesting differences across subjects and tasks. Visual and statistical differences between subjects (across tasks) were found in just their loops, although higher-dimensional or more local-scale topology may need to be analyzed to be able to identify subjects based on their mental state topologies. Inter-related loops emerged across tasks, and fewer loops, compared to functional networks, were needed to describe the fMRI data. Three main signatures of loop structures existed across all subjects and tasks, suggesting largely shared functional strategies across subjects and domains, and these signatures were low dimensional in their variation. Overall our analyses point towards the utility of high (above 0) dimensional topology in finding task-related and subject-specific patterns in fMRI spatial patterns.

Nevertheless, important future directions of this exploratory analysis include

1. seeing if higher-dimensional topological features can improve subject or task discrimination,
2. checking if low-persistence features can distinguish between subjects or tasks by decreasing how conservative the bootstrap procedure is,
3. developing techniques for finding representative topological features when combining multiple diagrams (like in the combined resting state 1 diagram across all subjects),
4. creating methods for finding the neurological dimension which varies periodically across a loop, and determining if temporal patterns exist in this dimension,
5. determining if HCP subjects whose mental state topologies were non-independent are blood related, and
6. correlating HCP behavioral measures with the MDS embedding coordinates of the mental state topologies to see if behavior can explain variance in mental state topologies.