Computing the Smith Normal Form of a matrix, and homology groups?

A simple algorithm for computing the Smith Normal Form of a matrix in $\ZZ$

The proof and the algorithm are the same.

First, a few auxiliary functions. Given a matrix $M$, the follow two functions are self-explanatory.

def dims(M):
 num_righe=len(M)
 num_colonne=len(M[0])
 return (num_righe,num_colonne)

def MinAij(M,s):
 num_righe, num_colonne=dims(M)
 ijmin=[s,s]
 valmin=max( max([abs(x) for x in M[j][s:]]) for j in range(s,num_righe) )
 for i in (range(s,num_righe)):
  for j in (range(s,num_colonne)):
   if (M[i][j] != 0 ) and (abs(M[i][j]) <= valmin) :
    ijmin = [i,j]
    valmin = abs(M[i][j])
 return ijmin

def IdentityMatrix(n):
 res=[[0 for j in range(n)] for i in range(n)]
 for i in range(n):
  res[i][i] = 1
 return res

def display(M):
 r=""
 for x in M:
  r += "%s\n" % x
 return r +""

Then, one needs the elementary operations on rows and columns on the matrix $M$: swap (permute) two rows, add to a row an integer multiple of another row, and change sign of a row. The same for columns.

def swap_rows(M,i,j):
 tmp=M[i]
 M[i]=M[j]
 M[j]=tmp

def swap_columns(M,i,j):
 num_of_columns=len(M)
 for x in range(num_of_columns):
  tmp=M[x][i]
  M[x][i] = M[x][j]
  M[x][j] = tmp

def add_to_row(M,x,k,s):
 num_righe,num_colonne=dims(M)
 for tmpj in range(num_colonne):
  M[x][tmpj] += k * M[s][tmpj]

def add_to_column(M,x,k,s):
 num_righe,num_colonne=dims(M)
 for tmpj in range(num_righe):
  M[tmpj][x] += k * M[tmpj][s]

def change_sign_row(M,x):
 num_righe,num_colonne=dims(M)
 for tmpj in range(num_colonne):
  M[x][tmpj] = - M[x][tmpj]

def change_sign_column(M,x):
 num_righe,num_colonne=dims(M)
 for tmpj in range(num_righe):
  M[tmpj][x] = - M[tmpj][x]

Finally, two important functions. The first, checks that the entry at place $(s,s)$ in the matrix $M$ is the only one non-zero in its $s$-th column under $s$ and $s$-th row at the right of $s$. If this is true, then the algorithm will proceed.

The second, will scan in the submatrix at the lower-right of $M_ {s,s}$, to check if there is any entry which is not divisible by $M_ {s,s}$. If it cannot be found, then it returns an empty set.

def is_lone(M,s):
 num_righe,num_colonne=dims(M)
 if [M[s][x] for x in range(s+1,num_colonne) if M[s][x] != 0] + [ M[y][s] for y in range(s+1,num_righe) if M[y][s] != 0] == []:
  return True
 else:
  return False

def get_nextentry(M,s):
  # find and element which is not divisible by M[s][s]
  num_righe,num_colonne=dims(M)
  for x in range(s+1,num_righe):
   for y in range(s+1,num_colonne):
    if M[x][y] % M[s][s]  != 0:
     return (x,y)
  return None

Given all these functions, the algorithm/proof is rather simple, as written in the folliwing block. The main idea is to put at the UL corner $M_ {0,0}$ the GCD of all matrix entries, by row/column operations.

First, find the smallest non-zero entry, and put it in the UL corner. Then, as in the Euclidean Algorithm, sum to all rows and columns suitable multiples of the first row/column. After this, either all the elements in the $0$-th row and $0$-th column are zero except $M_ {0,0}$, or not. If all the elements in the $0$-th row and $0$-th column V except $M_ {0,0}$ are zero, then either $M_ {0,0}$ divides all the entries of $M$, or not.

If it does divide all, then proceed to the submatrix obtained by deleting the first row and the first column. It is important to note that since $M_ {0,0}$ divides all the entries of $M$, it will divide all the entries of the submatrix and also all the entries of the subsequente modifications of the submatrix. Hence, when at the end the matrix will be in diagonal form, the first term $M_ {0,0}$ will divide all the others.

If it does not divide all, then add to the first row the row containing the entry which is not divisible by $M_ {0,0}$. So, the new matrix is of the type with $M_ {0,0}$ at the UL corner, and not all the elements in the $0$-th row and $0$-th column except $M_ {0,0}$ are zero.

Now go back to the step of finding the new non-zero minimum, put it in $(0,0)$, and repeat from there.

Why in a finite number of steps $M_ {0,0}$ has to divide all the entries? At the UL corner, $M_ {0,0}$ was chosen to be the minimum of non-zero terms of $M$, and the new one will be $M'_ {0,0} \leq M_ {0,0}$. But, if there was an entry
not divisible by $M_ {0,0}$, then $M'_ {0,0} < M_ {0,0}$. We show this with an example. Assume $a$ does not divide $b$ (and hence $a\neq b$), both positive. First, add the second row to the first. $$ \begin{bmatrix} a & 0 \\\ 0 & b \end{bmatrix} \mapsto \begin{bmatrix} a & b \\\ 0 & b \end{bmatrix} $$ Then, there are two cases: either $a < b$, and we do nothing, or $a > b$, and we permute the two columns. $$ \mapsto \begin{cases} \begin{bmatrix} a & b \\\ 0 & b \end{bmatrix} & \text{if $a<b$} \\\ \begin{bmatrix} b & a \\\ b & 0 \end{bmatrix} & \text{if $b < a$} \end{cases} $$ After this, it is time to add suitable multiples. In the two cases, it happens (adding forst rows, then columns) the following. $$ \mapsto \begin{cases} \begin{bmatrix} a & b-qa \\\ 0 & b \end{bmatrix} & \text{if $a<b$} \\\ \begin{bmatrix} b & a-qb \\\ 0 & -a \end{bmatrix} & \text{if $b < a$} \end{cases} $$ If $a<b$, it is $|b-qa|<a$. If $b<a$, it is $|a-qb| < b < a$. Hence in both cases after the row/column operations, the minimum of the matrix entries will be strictly smaller than $a$. Since this cannot happen indefinitively, after a finite number of steps $M_ {0,0}$ will divide all the remaining entries in $M$.

def Smith(M):
 num_righe,num_colonne=dims(M)
 L = IdentityMatrix(num_righe)
 R = IdentityMatrix(num_colonne)
 maxs=min(num_righe,num_colonne)
 for s in range(maxs):
  print ("step %s/%s\n" % (s+1,maxs))
  print "M:", display(M)
  while not is_lone(M,s):
   i,j = MinAij(M,s) # the non-zero entry with min |.|
   swap_rows(M,s,i)
   swap_rows(L,s,i)
   swap_columns(M,s,j)
   swap_columns(R,s,j)
   for x in range(s+1,num_righe):
    if M[x][s] != 0:
     k = M[x][s] // M[s][s]
     add_to_row(M,x,-k,s)
     add_to_row(L,x,-k,s)
   for x in range(s+1,num_colonne):
    if M[s][x] != 0:
     k = M[s][x] // M[s][s]
     add_to_column(M,x,-k,s)
     add_to_column(R,x,-k,s)
   if is_lone(M,s):
    res=get_nextentry(M,s)
    if res:
     x,y=res
     add_to_row(M,s,1,x)
     add_to_row(L,s,1,x)
    else:
     if M[s][s]<0:
      change_sign_row(M,s)
      change_sign_row(L,s)
 return L,R

A minor step is the one for changing sign of the diagonal entry, just in the case it is negative. At the end, the matrix $M$ will be in Smith Normal Form, and the routine will return the two unimodular matrices $L$ and $R$ resulting from the row/columns operations respectively. They are obtained by applying the same row/column operations on two identity matrices.

Example: the triangulated sphere.

As an example, here we compute the Smith Normal Forms of the differentials of che chain complex of the triangulate $2$-sphere, from sage. First, the differential $\partial_ 1 \from C_ 1 \cong \ZZ^6 \to C_ 0 \cong \ZZ^4 $.

sage: K=simplicial_complexes.Sphere(2)
sage: CK=K.chain_complex()
sage: 
sage: M=matrix(ZZ,CK.differential(1)).numpy().tolist()
sage: Smith(M)
step 1/4

M: [-1, -1, -1, 0, 0, 0]
[1, 0, 0, -1, -1, 0]
[0, 1, 0, 1, 0, -1]
[0, 0, 1, 0, 1, 1]

step 2/4

M: [1, 0, 0, 0, 0, 0]
[0, 0, 0, -1, -1, 1]
[0, 1, 1, 1, 1, 0]
[0, -1, -1, 0, 0, -1]

step 3/4

M: [1, 0, 0, 0, 0, 0]
[0, 1, 0, 0, 0, 0]
[0, 0, 1, 1, 1, 1]
[0, 0, -1, -1, -1, -1]

step 4/4

M: [1, 0, 0, 0, 0, 0]
[0, 1, 0, 0, 0, 0]
[0, 0, 1, 0, 0, 0]
[0, 0, 0, 0, 0, 0]

([[0, 0, 0, 1], [-1, 0, 0, 0], [-1, -1, 0, 0], [1, 1, 1, 1]],
 [[0, 1, -1, 1, 1, 0],
  [0, 0, 1, -1, -1, -1],
  [0, 0, 0, 0, 0, 1],
  [0, 0, 0, 1, 0, 0],
  [0, 0, 0, 0, 1, 0],
  [1, 0, 0, 0, -1, -1]])
sage: print "\nsmith form of d_1:\n", display(M)

smith form of d_1:
[1, 0, 0, 0, 0, 0]
[0, 1, 0, 0, 0, 0]
[0, 0, 1, 0, 0, 0]
[0, 0, 0, 0, 0, 0]

Hence $\operatorname{coker} \partial_ 1 \cong \ZZ^1$ and $\ker \partial_ 1 \cong \ZZ^3$. What about the generators of $\ker \partial_ 1$? Let $L$ and $R$ denote the left and right matrices, and $M$ the matrix of $\partial_ 1$. Then let $S$ denote the Smith Normal Form of $\partial_ 1$, so that $LMR = S$. Since $\ker S$ is generated by the three vectors $$ \begin{bmatrix} 0\\\ 0\\\ 0\\\ 1\\\ 0\\\ 0 \end{bmatrix} \begin{bmatrix} 0\\\ 0\\\ 0\\\ 0\\\ 1\\\ 0 \end{bmatrix} \begin{bmatrix} 0\\\ 0\\\ 0\\\ 0\\\ 0\\\ 1 \end{bmatrix}, $$ then $\ker LM = \ker M$ is generated by their images by $R$, i.e. the last three columns of $R$: $$ \begin{bmatrix} 1\\\ -1\\\ 0\\\ 1\\\ 0\\\ 0 \end{bmatrix} \begin{bmatrix} 1\\\ -1\\\ 0\\\ 0\\\ 1\\\ -1 \end{bmatrix} \begin{bmatrix} 0\\\ -1\\\ 1\\\ 0\\\ 0\\\ -1 \end{bmatrix}, $$ These three cycles are, in the standard basis of $1$-simplices of $K$ $[(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]$: $$ [01]+[20]+[12] = [01] + [12] + [20]~, $$ $$ [01]+[20]+[13]+[32] = [01] + [13] + [32] + [21]~, $$ $$ [20] + [03] + [32]~. $$

Furthermore, about the differential $\partial_ 2$ from $C_ 2 \cong \ZZ^4 \to C_ 1 \cong \ZZ^6 $.

sage: M=matrix(ZZ,CK.differential(2)).numpy().tolist()
sage: Smith(M)
step 1/4

M: [1, 1, 0, 0]
[-1, 0, 1, 0]
[0, -1, -1, 0]
[1, 0, 0, 1]
[0, 1, 0, -1]
[0, 0, 1, 1]

step 2/4

M: [1, 0, 0, 0]
[0, 0, 1, -1]
[0, -1, -1, 0]
[0, 0, -1, 1]
[0, 1, 1, 0]
[0, 1, 0, 1]

step 3/4

M: [1, 0, 0, 0]
[0, 1, 0, 0]
[0, 0, -1, -1]
[0, 0, -1, -1]
[0, 0, 1, 1]
[0, 0, 1, 1]

step 4/4

M: [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, 1],
  [1, 0, 0, 0, 0, 0],
  [1, 1, 0, 0, 0, 0],
  [0, 1, 0, 1, 0, -1],
  [-1, -1, 0, 0, 1, 1],
  [1, 1, 1, 0, 0, 0]],
 [[0, 1, -1, 1], [0, 0, 1, -1], [0, 0, 0, 1], [1, 0, 0, -1]])
sage: print "\nsmith form of d_2:\n", display(M)

smith form of d_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]

Therefore, $\ker \partial_ 2 \cong \ZZ$, and $\operatorname{coker} \partial_ 2 \cong \ZZ^3$. The generator of $\ker \partial_ 2$ is, as above, the last column of $R$, which is $$ \begin{bmatrix} 1\\\ -1\\\ 1\\\ -1\\\ \end{bmatrix} $$ which corresponds to the cycle in the standard basis $ [ (0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3) ] $ $$ [ 012 ] - [ 013 ] + [ 023 ] - [ 123 ] = - ( [ 123 ] - [ 023 ] + [ 013 ] - [ 012 ] ) = $$ $$ = - \sum_ {j=0}^3 (-1)^j [ 0 \ldots \hat j \ldots 3 ]~. $$

Is it possible to compute the homology with these data?

Computing homology groups

Now, consider just $H_ 1(K) = \ker \partial_ 1 / \operatorname{Im} \partial_ 2$. Let $Z_ 1 = \ker \partial_ 1$ and $B_ 1 = \operatorname{Im} \partial_ 2$. Let $S_ 1$ and $S_ 2$ the SNF’s of $ \partial_ 1 $ and $ \partial_ 2 $, and $L_ 1$ , $L_ 2$, $R_ 1$ and $R_ 2$ the unimodular matrices of row/columns operations. Then $L_ 1 \partial_ 1 R_ 1 = S_ 1$ and $L_ 2 \partial_ 2 R_ 2 = S_ 2$, as in the following diagram.

Let $\alpha_ 1$, $\alpha_ 2$, $\ldots$, $\alpha_ q$ be the elementary divisors of $\partial_ 2$, i.e. the non-zero diagonal terms of $S_ 2$, and $\beta_ 1$, $\beta_ 2$, $\ldots$, $\beta_ r$ the elementary divisors of $\partial_ 1$, i.e. the non-zero diagonal terms of $S_ 1$. Let $\ve_ 1$, $\ve_ 2$, $\ldots$, $\ve_ {c_ 1}$ be the standard basis of $\ZZ^{c_ 1}$.

Since $\ker S_ 1$ is generated by the $c_ 1 - r$ elements $\ve_ {r+1}$, $\ldots$, $\ve_ {c_ 1}$, as show above, $Z_ 1$ is freely generated by the last $c_ 1 - r$ columns of $R_ 1$.

On the other hand, $\operatorname{Im} S_ 2$ is freely generated by $$ \alpha_ 1 \ve_ 1, \alpha_ 2 \ve_ 2, \ldots, \alpha_ q \ve_ q $$ hence $B_ 1\subset Z_ 1 \subset C_ 1$ is freely generated by $\alpha_ j$-multiples of the first $q$ columns of the inverse $L_ 2^{-1}$. Now, $C_ 1 = Q \oplus Q'$, where $Q$ is the subgroup generated by the first $q$ columns of $L_ 2^{-1}$, and $Q'$ the subgroup generated by the last $c_ 1 - q$ columns of $L_ 2^{-1}$.

By construction $B_ 1 \subset Q$, and since $C_ 1$ is torsion-free, $Q\subset Z_ 1$. Now, define $Q'' = Q' \cap Z_ 1$, which is a subgroup of $Z_ 1$. Since $Q'' \cap Q = 0$, it is defined the direct sum $Q\oplus Q'' \subset Z_ 1$. If $x \in Z_ 1 \subset C_ 1 = Q \oplus Q'$, then $x=q + q'$ with $q \in Q$ and $q' \in Q'$, and $0 = \partial_ 1 x = \partial_ 1 q + \partial_ 1 q' = \partial_ 1 q'$, hence $x=q + q'$ with $q \in Q$ and $q' \in Q' \cap Z_ 1 = Q''$, and hence $Z_ 1 = Q \oplus Q''$.

Since $Z_ 1=Q\oplus Q''$ is a free abelian group of rank $c_ 1 - r$, and $Q$ is a free abelian subgroup of rank $q$,

(1)

$Q''\subset Z_ 1$ is a free abelian group of rank $c_ 1 - r - q$.

Furthermore, since $B_ 1 \subset Q$, the quotient $Z_ 1 / B_ 1$ splits as $$ \dfrac{Z_ 1 }{ B_ 1 } = \dfrac{Q \oplus Q''}{ B_ 1\oplus 0} \cong \dfrac{Q}{ B_ 1 } \oplus Q''
\cong $$ $$ \cong \ZZ_ {\alpha_ 1} \oplus \ldots \oplus \ZZ_ {\alpha_ q} \oplus \ZZ^{c_ 1 - r - q}, $$ where $\ZZ_ \alpha$ denotes the cyclic group $\ZZ / \alpha \ZZ$ of order $\alpha$.

(2)

The homology group $H_ j$ is isomorphic to $ \ZZ_ {\alpha_ 1} \oplus \ldots \oplus \ZZ_ {\alpha_ q} \oplus \ZZ^{c_ j - r - q}$, where $\alpha_ 1$, $\ldots,$ $\alpha_ q$ are the elementary divisors of $S_ {j+1}$, $c_ j$ the rank of $C_ j$, and $r$ the number of elementary divisors of $S_ j$.

So, for the chain complex of the sphere: $c_ 2 = 4$, $c_ 1 = 6$, $c_ 0 = 4$; from $S_ 1$: $\beta_ 1 = \beta_ 2 = \beta_ 3 = 1$, $r= 3$; from $S_ 2$: $\alpha_ 1 = \alpha_ 2 = \alpha_ 3= 1$, $q=3$.

Hence, for the sphere, $$ \begin{aligned} H_ 0 & \cong \ZZ_ 1 \oplus \ZZ_ 1 \oplus \ZZ_ 1 \oplus \ZZ^{4-0-3} = \ZZ \\\ H_ 1 & \cong \ZZ_ 1 \oplus \ZZ_ 1 \oplus \ZZ_ 1 \oplus \ZZ^{6-3-3} = 0 \\\ H_ 2 & \cong \ZZ^{4-0-3} = \ZZ . \end{aligned} $$

Remark: We did not compute explicit isomorphisms of $H_ j \cong Z_ j / B_ j$. In fact, we even did not compute the inverse $L_ j^{-1}$. To have an explicit presentation (generators and relations) of $H_ j$ one needs to proceed as follows.

First, generators of $\ker S_ 1$ in $\ZZ^{c_ 1}$ are $\ve_ {q+1}$, $\ldots,$ $\ve_ {c_ 1}$. Next, generators of $R_ 1^{-1}Q$ are $\vw_ j = R_ 1^{-1} L_ 2^{-1} \ve_ j$ for $j=1,\ldots, q$. Since for each $j=1,\ldots, q$, $L_ 2^{-1} \ve_ j \in Q$, hence $$ S_ 1 \vw_ j = S_ 1 R_ 1^{-1} L_ 2^{-1} \ve_ j = L_ 1 \partial_ 1 L_ 2^{-1} \ve_ j = 0, $$ which means that the first $r$ entries of $\vw_ j$ are zero. Note that the projection $p\from Z_ 1 = Q\oplus Q'' \to Q $, after $R_ 1^{-1}$, is a projection $ \ZZ^{c_ 1 - r} \to \ZZ^{q}$. Its kernel is $Q''$.

So, let $T$ denote the $(c_ 1 -r ) \times q$ matrix whose $q$ columns are the last $c_ 1 - r$ entries of $\vw_ j$, for $j=1,\ldots, q$. It represents the monomorphism $T\from R_ 1^{-1} Q \to \ker S_ 1$, or $T \from \ZZ^q \to \ZZ^{c_ 1 - r}.$ Let $S_ T$ its Smith Normal Form, and $L_ T$ and $R_ T$ the unimodular matrices such that $S_ T = L_ T T R_ T$. By the arguments above, $S_ T $ has $q$ $1$’s as elementary divisors, and the rest are $0$. Furthermore, the unimodular matrix $L_ T ^{-1} \from \ZZ^{c_ 1 - r} \to \ZZ^{c_ 1 - r}$
sends the first $q$ elements of the standard basis $\ve_ j$, $j=1,\ldots, q$, to a set of generators for $Q$ (the image of $T$ in $\ker S_ 1$). The latter $c_ 1 - r - q$ columns of $L_ T^{-1}$ are the coefficients of a set of generators of $Q''$ in the standard basis of $\ker S_ 1$, which we can denote by $\vw_ {q+1}, \ldots, \vw_ {c_ 1 - r}$ (by adding $q$ zeros at the head of the $(c_ 1 - r)$-tuple).

So, going back to $C_ 1$, $Z_ 1$ is generated by the elements of $Q$, i.e. $$ R_ 1 \vw_ 1 = L_ 2^{-1} \ve_ 1, \ldots R_ 1 \vw_ q = L_ 2^{-1} \ve_ q, $$ and the elements of \( Q'' \), i.e.
$$ R_ 1 \vw_ {q+1} = R_ 1 L_ T^{-1} \ve_ {q+1} \ldots R_ 1 \vw_ {c_ 1 - r} = R_ 1 L_ T^{-1} \ve_ {c_ 1 - r} $$ where suitable zeros are added when necessary. More in Computational Homology, by Tomasz Kaczynski, Konstantin Mischaikow, Marian Mrozek, Springer, 2004.