Home » maths » Computing Homology

# Computing Homology

Last time we introduced abstract simplicial complexes (abs) and how to represent them in Python. Today we will look at the homology groups of an abs.

For an excellent primer on homology I recommend checking out this post , I will just quickly recap the basics.

Homology: A quick reminder

Let $X$ be an abs.

Write $X_i$ for the $i^{th}$dimensional skeleton of $X$, formed by restricting to those elements of $X$ with dimension $\le i$.

If $\sigma \in X$ is a $p+1$-simplex with vertices labelled $0, \dots, p$ we may write $\sigma = [0, \dots, p]$ to denote this. To denote the $p$-simplex consisting of the vertices of $\sigma$ except the $i^{th}$ vertex we write $[0, \dots, \hat{i}, \dots p]$.

For each dimension $p$ we write $N_p$ for the set of $p$-dimensional simplices, and write $n_p = |N_p|$ for its cardinality.

For each dimension $p$ we write $C_p$ for the $p$dimensional chain group given by the group presentation $C_p =$. What this says is that the generators of the group are the $p$-simplices, the group is abelian and each element is its own inverse.

An alternative way to think of this is as a vector space, where $C_p$ has basis $N_p$ and is over the field $\mathbb{Z}_2$ (mod 2 arithmetic). So every element can be represented as a binary vector of length $n_p$, and addition is bitwise xor.

Now we define boundary maps $\delta_{p+1}: C_{p+1} \to C_p$. These are both group homomorphisms and linear maps. For a $p+1$-simplex $\sigma = [0,1, \dots, p]$ define: $\delta_{p+1}\sigma = \sum \limits _{i=0}^{p}[0, \dots, \hat{i}, \dots p]$, the sum of its $p$-dimensional faces.

Since $N_{p+1}$ is a a basis for $C_{p+1}$ this determines a linear mapping.

Define $Z_{p+1} := ker \, \delta_{p+1}$ called the $p+1$dimensional cycle group (subspace).

Define $B_p: = Im \, \delta_{p+1}$ called the $p$dimensional boundary group (subspace).

Now we can define the $p$dimensional homology group $H_p:= \frac{Z_p}{B_p}$.

This is both a quotient group and a quotient vector space.

Finally define the $p^{th}$ Betti number $B_p:= \text{rank }H_p = \text{rank } Z_p - \text{rank }B_p$.

Computation

We are going to be interested in the Betti numbers as these can be thought of as counting the $p$-dimensional holes in our space. For $B_0$ this just means the number of connected components.

In order to calculate the ranks of interest, we must construct the matrices representing the boundary mappings. In order to do this we must agree upon an ordering of the simplices. So we simply use the one generated by the data structure, which is why each simplex has an index attribute.

To find the rank of $C_p$, which we have called $n_p$, we simply count the number of $p$-simplices. So we can add that feature to our SimplicialComplex class easily.

    def chain_rank(self,p):
return len(self.n_simplex_dict[p])


Now the matrix of $\delta_p: C_p \to C_{p-1}$ will be have $n_{p-1}$ rows, one for each $p-1$-simplex, and $n_p$ columns, one for each $p$-simplex. Then the $i^{th}$ column of the matrix will contain the coefficients of the image of the $i^{th} \ p-$simplex expressed in terms of our basis for $C_{p-1}$.

So to find the $i^{th}$ column vector, we ask each of a simplex’s boundary faces for their index. These will be the positions of the 1s in the column.

    def boundary_column(self,simplex):
column_size = self.chain_rank(simplex.dimension-1)
faces = simplex.faces
column = [0 for i in range(column_size)]
for face in faces:
column[face.index]=1
return column


Now we simply combine these columns to form the boundary matrix. Note that for $p$ greater than the dimension of the simplicial complex, $C_p=0$ the trivial group, and similarly for negative $p$. We don’t have to calculate anything in this case.

This also means that the kernel of $\delta_0$ is $C_0$. This is handled by us having the empty set as a kind of $-1$-dimensional simplex.

    def get_boundary_matrix(self,p, skeleton=None):
if skeleton==None:
skeleton=self.dimension
if ((p>skeleton)or(p<0)):
return np.asarray([])
mat=[self.boundary_column(simplex) for simplex in self.n_simplex_dict[p]]
return np.asarray(mat, dtype=int).transpose()


The skeleton argument is used when you want to only consider simplices up to a certain dimension. For instance if $X$ is given by $[0,1,2]$ then $X$ is a representation of a triangle. Then the zeroth Betti number $B_0=1$ as there is one connected component. However the zeorth Betti number of $X_0$ is 3, one for each vertex.

Let’s test what we’ve got so far. In this post Jeremy Kun computes the homology of the tetrahedron with an ‘arm’, he provides a nice picture of it here. If we label the vertices of the tetrahedron $0,1,2,3$, then this structure can be specified with the maximal simplices $[0,1,2,3], [0,4], [2,4$ then restricting to the $2$-skeleton, since we want there to be a hole in the tetrahedron.

Let’s see what the boundary matrices are:

def jeremy_example():
maximal_simplices = [[0,1,2,3],[0,4],[2,4]]
S = SimplicialComplex(maximal_simplices)
S.plot(direction='up', label_color='blue')
for i in range(3):
print "Boundary matrix of dimension ", i
print S.get_boundary_matrix(i,skeleton=2)

jeremy_example()


Which produced the following plot (of the 3-skeleton): And outputted to the console:

Boundary matrix of dimension 0
[[1 1 1 1 1]]
Boundary matrix of dimension 1
[[1 1 1 0 0 0 1 0]
[1 0 0 1 1 0 0 0]
[0 1 0 1 0 1 0 1]
[0 0 1 0 1 1 0 0]
[0 0 0 0 0 0 1 1]]
Boundary matrix of dimension 2
[[1 1 0 0]
[1 0 1 0]
[0 1 1 0]
[1 0 0 1]
[0 1 0 1]
[0 0 1 1]
[0 0 0 0]
[0 0 0 0]]

If you check his post, you will see that the matrices are the same, up to a reordering of basis, and noticing that $1 \equiv -1 \ (mod \ 2)$

So far so good. Now all that remains is to calculate the rank and nullity of each matrix. But this is easily done, we simply diagonalize the matrix and count the pivots. For $\delta_{p+1}$ the number of pivots is equal to $\text{rank } B_p$ and $n_{p+1} - \text{rank }B_p = \text{ rank } Z_{p+1}$ by the rank-nullity theorem.

The form we will put the matrix in is called Smith-Normal form I won’t go through the details of how the algorithm works, they are on the wiki page. If you have seen Gaussian elimination before this is basically that.

    def reduce_matrix(self, matrix):
#Returns [reduced_matrix, rank, nullity]
if np.size(matrix)==0:
return [matrix,0,0]
m=matrix.shape
n=matrix.shape

def _reduce(x):
#We recurse through the diagonal entries.
#We move a 1 to the diagonal entry, then
#knock out any other 1s in the same  col/row.
#The rank is the number of nonzero pivots,
#so when we run out of nonzero diagonal entries, we will
#know the rank.
nonzero=False
#Searching for a nonzero entry then moving it to the diagonal.
for i in range(x,m):
for j in range(x,n):
if matrix[i,j]==1:
matrix[[x,i],:]=matrix[[i,x],:]
matrix[:,[x,j]]=matrix[:,[j,x]]
nonzero=True
break
if nonzero:
break
#Knocking out other nonzero entries.
if nonzero:
for i in range(x+1,m):
if matrix[i,x]==1:
matrix[i,:] = np.logical_xor(matrix[x,:], matrix[i,:])
for i in range(x+1,n):
if matrix[x,i]==1:
matrix[:,i] = np.logical_xor(matrix[:,x], matrix[:,i])
#Proceeding to next diagonal entry.
return _reduce(x+1)
else:
#Run out of nonzero entries so done.
return x
rank=_reduce(0)
return [matrix, rank, n-rank]

def boundary_rank(self,p, skeleton=None):
if skeleton==None:
skeleton=self.dimension
if (p >= skeleton) or (p<0):
#Trivial case.
return 0
else:
return self.reduce_matrix(self.get_boundary_matrix(p+1, skeleton))

def cycle_rank(self, p, skeleton=None):
if skeleton==None:
skeleton=self.dimension
if p==0:
#Mapping to trivial group so kernel is everything in this case.
return self.chain_rank(0)
elif (p>skeleton) or (p<0):
#Trivial case.
return 0
else:
return self.reduce_matrix(self.get_boundary_matrix(p, skeleton))

def betti_number(self, p, skeleton=None):
if skeleton==None:
skeleton=self.dimension
return self.cycle_rank(p,skeleton) - self.boundary_rank(p, skeleton)


Now let’s just test it out on the same example. We expect $B_0=1$ since there is one connected component, $B_1=1$ from the $1$-dimensional hole made by the arm and $B_2=1$ for the ball inside the tetrahedron’s shell.

def jeremy_example(skeleton=2):
maximal_simplices = [[0,1,2,3],[0,4],[2,4]]
S = SimplicialComplex(maximal_simplices)
S.plot(direction='up', label_color='blue')
for i in range(4):
print "Boundary matrix of dimension", i
print S.get_boundary_matrix(i,skeleton)
print "Reduced boundary matrix of dimension", i
print S.reduce_matrix(S.get_boundary_matrix(i,skeleton))
print "Betti number",i,"=", S.betti_number(i,skeleton)

jeremy_example()


The console outputs:

Boundary matrix of dimension 0
[[1 1 1 1 1]]
Reduced boundary matrix of dimension 0
[[1 0 0 0 0]]
Betti number 0 = 1
Boundary matrix of dimension 1
[[1 1 1 0 0 0 1 0]
[1 0 0 1 1 0 0 0]
[0 1 0 1 0 1 0 1]
[0 0 1 0 1 1 0 0]
[0 0 0 0 0 0 1 1]]
Reduced boundary matrix of dimension 1
[[1 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0]
[0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 0 0]]
Betti number 1 = 1
Boundary matrix of dimension 2
[[1 1 0 0]
[1 0 1 0]
[0 1 1 0]
[1 0 0 1]
[0 1 0 1]
[0 0 1 1]
[0 0 0 0]
[0 0 0 0]]
Reduced boundary matrix of dimension 2
[[1 0 0 0]
[0 1 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]]
Betti number 2 = 1
Boundary matrix of dimension 3
[]
Reduced boundary matrix of dimension 3
[]
Betti number 3 = 0

Just as expected. So we’re done. Get the code here.