Persistent  Homology  Template

A Point Cloud

Here is a sample point cloud of ten points in (ℝ^2):

points={{0.335371,0.729228},{0.93472,0.644295},{0.493764,0.921443},{0.831404,0.055325},{0.340745,0.198575},{0.764133,0.076333},{0.251056,0.760276},{0.0211607,0.495884},{0.0969048,0.284578},{0.0585317,0.16929}};
pointPlot=ListPlot[points->Range[Length[points]],PlotRange->{{0,1},{0,1}},AspectRatio->1]

Vietoris-Rips Filtration

The following function computes the Vietoris-Rips filtration of a given point cloud. The optional parameter maxDist specifices the maximum length of edges to include, which limits the size of the simplicial complexes in the filtration.

computeVRfiltration[points_,maxDist_:∞]:=Module[
{n=Length[points],
dists,
filtration,
dist,triangles
},

(* compute all edge lengths *) 
dists=Table[0,n,n];
Do[
dists[[i,j]]=EuclideanDistance[points[[i]],points[[j]]]
,{i,1,n},{j,i+1,n}];

(* build the simplicial complex *)
filtration=Table[{i}->0,{i,1,n}];
Do[
If[dists[[i,j]]<=maxDist,
AppendTo[filtration,{i,j}->dists[[i,j]]];
Do[
If[dists[[i,k]]<=maxDist&&dists[[j,k]]<=maxDist,
AppendTo[filtration,{i,j,k}->Max[dists[[i,j]],dists[[i,k]],dists[[j,k]]]]
]
,{k,j+1,n}]
]
,{i,1,n},{j,i+1,n}];

(* return the simplicial complex *) 
SortBy[filtration,Last]
]

Compute the filtration for the given point cloud:

filtration=computeVRfiltration[points,0.8]

The next function draws a sequence of simplicial complexes to visualize the filtration.

visualizeVRFiltration[filtration_,points_]:=Module[
{vals=DeleteDuplicates[Values[filtration]],plots={},complex, verts, edges, faces, plot},

(* loop over filtration values in increasing order *)
Do[
complex=Keys[Select[filtration,Values[#]<=v&]];
verts={};edges={};faces={};

(* make lists of simplices of each dimension *) 
Do[
If[Length[simplex]==1,AppendTo[verts,points[[simplex[[1]]]]]];
If[Length[simplex]==2,AppendTo[edges,points[[simplex]]]];
If[Length[simplex]==3,AppendTo[faces,points[[simplex]]]]
,{simplex,complex}];

(* make the drawing *) 
plot=Show[
ListPlot[points->Range[Length[points]],PlotRange->{{0,1},{0,1}},AspectRatio->1,PlotStyle->{PointSize[Medium],Black},Axes->False],
Graphics[{Blue,Opacity[0.4],Triangle[faces],Black,Opacity[1],Line[edges],Text["ϵ = "<>ToString[v],{0.01,0.99},{Left,Top}]}] 
];
AppendTo[plots,plot]
,{v,vals}];

(* return the sequence of plots *) 
plots
]

Let's visualize our sample filtration:

plotSequence=visualizeVRFiltration[filtration,points];
ListAnimate[plotSequence]

Boundary Matrices

The following function computes the boundary matrix of a filtration of a specified dimension. If the dimension is j, then the matrix specifies the boundaries of all j-dimensional simplices in the filtration. That is, the columns of the matrix correspond to j-dimensional simplices, arranged in order of increasing filtration value.

computeBoundaryMx[filtration_,dim_]:=Module[
{rowSimplices,colSimplices,matrix,faces,r},

(* get the lists of simplices in adjacent dimensions *) 
rowSimplices=Keys[Select[filtration,Length[First[#]]==dim&]];
colSimplices=Keys[Select[filtration,Length[First[#]]==dim+1&]];
Print["column simplices: ",colSimplices];
Print["row simplies: ",rowSimplices];

(* handle empty matrix *)
If[Length[rowSimplices]==0||Length[colSimplices]==0,Return[{{}}]];

(* initialize a matrix *) 
matrix=Table[0,Length[rowSimplices],Length[colSimplices]];

(* loop over all columns *)
Do[
(* loop over all codim 1 faces of this simplex *)
faces=Subsets[colSimplices[[c]],{dim}];
Do[
(* set the correct entry to 1 *) 
r=Position[rowSimplices,f,1][[1,1]];
matrix[[r,c]]=1;
,{f,faces}];
,{c,1,Length[colSimplices]}];

(* return the matrix*) 
matrix
]

Now we compute the boundary matrices of our filtration. The matrix (∂_1) gives the boundaries of all edges in the filtration:

d1=computeBoundaryMx[filtration,1];
MatrixForm[d1]

The matrix (∂_2) gives the boundaries of all triangular faces in the filtration:

d2=computeBoundaryMx[filtration,2];
MatrixForm[d2]

Reduce the Boundary Matrices

To finish the computation of persistent homology, write the following module that reduces a boundary matrix. Your module should accept a matrix as input and return the reduced matrix.

reduceMatrix[matrix_]:=Module[
{},

]

Test your function by reducing the boundary matrix (∂_1):

MatrixForm[d1]
MatrixForm[reduceMatrix[d1]]

This is the reduced form of matrix (∂_1):

(({{0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}))

Compute the Persistence Pairs

The following function computes the barcode of a given filtration. Note that this function calls the previously-defined functions computeBoundaryMx and reduceMatrix, defined above. The barcode is returned as a nested list of persistence pairs.

computeBarcode[filtration_]:=Module[
{d1, d2,numPoints,edges,faces,
low, pairs,val1,val2, h0bars={},h1bars={}
},

(* compute the boundary matrices *) 
d1=computeBoundaryMx[filtration,1];
d2=computeBoundaryMx[filtration,2];

(* reduce the boundary matrices *) 
d1=reduceMatrix[d1];
d2=reduceMatrix[d2];

(* get the (Subscript[H, 0 ])persistence pairs *)
(*Print[MatrixForm[d1]];*)
numPoints=Dimensions[d1][[1]];
edges=Select[filtration,Length[First[#]]==2&];
pairs=Table[0,numPoints];
Do[
low=getLow[d1,c];(* low entry in the current column *)
If[low>0,
pairs[[low]]=c;
(*Print["row ", low," paired with column ",c]*)
]
,{c,1,Length[edges]}];
Do[
If[pairs[[r]]==0,
AppendTo[h0bars,{0,∞}],
AppendTo[h0bars,{0,Values[edges][[pairs[[r]]]]}]]
,{r, 1, numPoints}] ;
(*Print["(Subscript[H, 0]) barcode: ",h0bars];*)

(* get the (Subscript[H, 1 ])persistence pairs *)
(*Print[MatrixForm[d2]];*)
faces=Select[filtration,Length[First[#]]==3&];
pairs=Table[0,Length[edges]];
Do[
low=getLow[d2,c];(* low entry in the current column *)
If[low>0,
pairs[[low]]=c;
(*Print["row ", low," paired with column ",c]*)
]
,{c,1,Length[faces]}];
Do[
If[pairs[[r]]>0,
val1=Values[edges][[r]];
val2=Values[faces][[pairs[[r]]]];
If[val2>val1,AppendTo[h1bars,{val1,val2}]]
(* WHAT ABOUT EDGES THAT CREATE ESSENTIAL CYCLES? *)
]
,{r, 1, Length[edges]}] ;

(* return the barcode *)
{ReverseSortBy[h0bars,Last],ReverseSortBy[h1bars,Last]}
]

Compute the barcode of our sample filtration:

barcode=computeBarcode[filtration]

Visualizing  the  Barcode

The following module draws a barcode. The input is a nested list of persistence pairs, and the output is a graphics object.

visualizeBarcode[barcode_]:=Module[
{vpos=0.05,offset=0.05,h0bars={},h1bars={},h0labelPos,h1labelPos,h0top},
(* make graphics objects for (Subscript[H, 0]) bars *)
Do[
If[int[[2]]<∞,
AppendTo[h0bars,Line[{{int[[1]],vpos},{int[[2]],vpos}}]],
AppendTo[h0bars,Arrow[{{int[[1]],vpos},{1,vpos}}]]
];
vpos+=offset
,{int,barcode[[1]]}];
h0top=vpos;
h0labelPos=h0top/2;
(* make graphics objects for (Subscript[H, 1]) bars *)
Do[
If[int[[2]]<∞,
AppendTo[h1bars,Line[{{int[[1]],vpos},{int[[2]],vpos}}]],
AppendTo[h1bars,Arrow[{{int[[1]],vpos},{1,vpos}}]]
];
vpos+=offset
,{int,barcode[[2]]}];
h1labelPos=(h0top+vpos-offset)/2;

(* return the graphics object *)
Graphics[{Blue, Thickness[0.03],CapForm["Butt"],Arrowheads[0.11],h0bars,Text["(Subscript[H, 0])",{-0.03,h0labelPos},{Right,Center}],Orange,h1bars,Text["(Subscript[H, 1])",{-0.03,h1labelPos},{Right,Center}],AxisObject[{"Horizontal",0}]},PlotRange->{{-0.1,1},{0,vpos}}]
]

Now we can finally draw the barcode for our sample filtration:

visualizeBarcode[barcode]