This is a simple Haskell implementation of Smith normal form for matrices over the integers, built upon the matrix library, and so emphasizes elegance and minimal dependencies over pure speed. In particular, it has no dependencies beyond those needed for the matrix library.
Linear algebra packages usually assume that you are working over a field, such as the rational numbers or the real numbers. Matrices over fields also have a Smith normal form, but as scaling a row by an arbitrary nonzero constant corresponds to an elementary matrix in linear algebra over a field, the resulting diagonal matrix is equivalent to one with r 1's and all other entries 0, where r is the rank of the matrix. More interesting information can be gleaned if we restrict ourselves to "linear algebra over the integers."
Example Usage
In ghci:
> import Data.Matrix
Data.Matrix> import Data.Matrix.SmithNormalForm
Data.Matrix Data.Matrix.SmithNormalForm> a = fromList 4 4 [-2, 1, 1, 0, 1, -2, 1, 0, 1, 1, -3, 1, 0, 0, 1, -1] :: Matrix Integer
Data.Matrix Data.Matrix.SmithNormalForm> a
┌ ┐
│ -2 1 1 0 │
│ 1 -2 1 0 │
│ 1 1 -3 1 │
│ 0 0 1 -1 │
└ ┘
Data.Matrix Data.Matrix.SmithNormalForm> smithNormalForm a
┌ ┐
│ 1 0 0 0 │
│ 0 1 0 0 │
│ 0 0 3 0 │
│ 0 0 0 0 │
└ ┘
Warning: For dense matrices beyond a small size, calculating the Smith normal form of a matrix can result in "coefficient explosion," that is, during the computation, the entries of the matrices involved in the computation can become very large, even when those of the initial matrix are small. So if your matrix has entries whose type consists of bounded integers (e.g. Int) as opposed to unbounded integers (e.g. Integer), you are likely to encounter integer overflow errors.
We recommend that you make sure that the matrices to which you're applying smithNormalForm have type Matrix Integer, especially if you found the remarks above baffling or if you're not sure what to do here.
In general, the main theorem regarding the existence and uniqueness of Smith normal form holds over any principal ideal domain (PID). The present implementation immediately extends to any Euclidean domain, but would require modifications for the general case.
The Smith normal form a matrix A can be thought of as the simplest matrix that is equivalent to A, in the sense that it can obtained by applying the kinds of transformations used in Gauss–Jordan elimination to both rows (probably the first thing you learn in linear algebra, in getting to reduced row echelon form) and columns (i.e. doing the same thing to the transposed matrix).
More precisely, the Smith normal form of a matrix A with entries in the integers is the unique matrix D (also with entries in the integers) with the following properties:
-
We have A = RDC where R and C are products of elementary matrices; here R and C can be taken to correspond to sequences of row and column operations, respectively. (Recall that the elementary matrices over the integers are matrices of determinant +1 or -1 that correspond to the operations of (i) multiplying a row or column by -1, (ii) switching two rows or two columns, or (iii) adding a scalar integer multiple of one row to another row or one column to another column.)
-
The entries of D are zero outside of the diagonal entries (i.e. D_{i,j} = 0 for i not equal to j).
-
The diagonal entries [d(1) = D_{1,1}, d(2) = D_{2,2},..,d(n) = D_{n,n}] of D are nonnegative and d(i) divides d(i+1) for i = 1, 2,..., n-1. (These numbers are called the invariant factors of the matrix A.)
- Compute the homology of a chain complex of finitely generated abelian groups.
- Determine the structure of a finitely generated module over the integers.
- Calculate the (abelian) sandpile group (a.k.a. critcal group, Jacobian group, Picard group) of a graph.
A number of combinatorial applications and related questions can be found in this 2016 survey by Stanley.
