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 be an abs.

Write for the dimensional skeleton of , formed by restricting to those elements of with dimension .

If is a -simplex with vertices labelled we may write to denote this. To denote the -simplex consisting of the vertices of except the vertex we write .

For each dimension we write for the set of -dimensional simplices, and write for its cardinality.

For each dimension we write for the –*dimensional chain group* given by the group presentation . What this says is that the generators of the group are the -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 has basis and is over the field (mod 2 arithmetic). So every element can be represented as a binary vector of length , and addition is bitwise xor.

Now we define boundary maps . These are both group homomorphisms and linear maps. For a -simplex define:

, the sum of its -dimensional faces.

Since is a a basis for this determines a linear mapping.

Define called the –*dimensional cycle group* (subspace).

Define called the –*dimensional boundary group* (subspace).

Now we can define the –*dimensional homology group* .

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

Finally define the *Betti number* .

**Computation**

We are going to be interested in the Betti numbers as these can be thought of as counting the -dimensional holes in our space. For 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 , which we have called , we simply count the number of -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 will be have rows, one for each -simplex, and columns, one for each -simplex. Then the column of the matrix will contain the coefficients of the image of the simplex expressed in terms of our basis for .

So to find the 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 greater than the dimension of the simplicial complex, the trivial group, and similarly for negative . We don’t have to calculate anything in this case.

This also means that the kernel of is . This is handled by us having the empty set as a kind of -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 is given by then is a representation of a triangle. Then the zeroth Betti number as there is one connected component. However the zeorth Betti number of 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 , then this structure can be specified with the maximal simplices then restricting to the -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

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 the number of pivots is equal to and 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[0] n=matrix.shape[1] 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))[1] 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))[2] 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 since there is one connected component, from the -dimensional hole made by the arm and 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))[0] 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.