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 pdimensional chain group given by the group presentation C_p =<N_p |[x,y], x=x^{-1},\ \forall \, x,y \in N_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+1dimensional cycle group (subspace).

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

Now we can define the pdimensional 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):

jeremy_example

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[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 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))[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.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: