SmithNormalForm (empty) → 0.1.0.0

raw patch · 6 files changed

+274/−0 lines, 6 filesdep +basedep +matrixdep +vector

Dependencies added: base, matrix, vector

Files

ChangeLog.md +4/−0
LICENSE +8/−0
README.md +60/−0
SmithNormalForm.cabal +40/−0
src/Data/Matrix/SmithNormalForm.hs +29/−0
src/Data/Matrix/SmithNormalForm/Internal.hs +133/−0

+ ChangeLog.md view

@@ -0,0 +1,4 @@+# Changelog for SmithNormalForm++## 1.0.0+Initial release.

+ LICENSE view

@@ -0,0 +1,8 @@+MIT License+Copyright (c) 2021 Brian Hwang++Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

+ README.md view

@@ -0,0 +1,60 @@+# SmithNormalForm++This is a simple Haskell implementation of Smith normal form for matrices over the *integers*, built upon the [matrix](https://hackage.haskell.org/package/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.++## What is Smith Normal Form?++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.)++++## Some Applications of Smith Normal Form+ +* 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](https://arxiv.org/abs/1602.00166).

+ SmithNormalForm.cabal view

@@ -0,0 +1,40 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.34.4.+--+-- see: https://github.com/sol/hpack++name:           SmithNormalForm+version:        0.1.0.0+synopsis:       A lightweight Haskell implementation of Smith normal form over the integers.+description:    This package contains an implementation of Smith normal form for arbitrary matrices over the integers. In particular, it does not assume that the matrix is invertible or square.+category:       Math, LinearAlgebra+homepage:       https://notabug.org/bwh/SmithNormalForm+bug-reports:    https://notabug.org/bwh/SmithNormalForm/issues+author:         Brian Hwang+maintainer:     bhwang@math.cornell.edu+copyright:      2021 Brian Hwang+license:        MIT+license-file:   LICENSE+build-type:     Simple+extra-source-files:+    README.md+    ChangeLog.md++source-repository head+  type:                git+  location:            https://notabug.org/bwh/SmithNormalForm.git++library+  exposed-modules:+      Data.Matrix.SmithNormalForm+      Data.Matrix.SmithNormalForm.Internal+  other-modules:+      Paths_SmithNormalForm+  hs-source-dirs:+      src+  build-depends:+      base >=4.7 && <5+    , matrix >=0.3.6 && <0.4+    , vector >=0.10 && <0.13+  default-language: Haskell2010

+ src/Data/Matrix/SmithNormalForm.hs view

@@ -0,0 +1,29 @@+{-+Module      : Data.Matrix.SmithNormalForm +Description : A simple implementation of Smith normal form over the integers.+Copyright   : (c) Brian Hwang 2021+License     : MIT+Stability   : stable+-}+module Data.Matrix.SmithNormalForm+    ( smithNormalForm+    , invariantFactors+    ) where++import qualified Data.Matrix as M+import qualified Data.Vector as V++import qualified Data.Matrix.SmithNormalForm.Internal++-- | Returns the Smith normal form of an matrix, i.e. a diagonal matrix+-- obtained by applying elementary row and column operations+-- whose diagonal entries \([d_1,..,d_n]\) are such that \(d_k | d_{k+1}\)+-- (Does not assume that the matrix is square.)+smithNormalForm :: Integral a => M.Matrix a -> M.Matrix a+smithNormalForm = Data.Matrix.SmithNormalForm.Internal.smithNF +++-- | Given a matrix M, returns the invariant factors of M, i.e. the list of+-- diagonal entries of the Smith normal form of M.+invariantFactors :: Integral a => M.Matrix a -> [a]+invariantFactors = V.toList . M.getDiag . smithNormalForm

+ src/Data/Matrix/SmithNormalForm/Internal.hs view

@@ -0,0 +1,133 @@+module Data.Matrix.SmithNormalForm.Internal+    ( smithNF+    , rectifyDiagonal+    , diagonalize+    , isDiagonalMatrix+    , divides+    ) where++import qualified Data.Matrix as M+import qualified Data.Vector as V+import Data.Maybe (fromJust)++-- | Main method that returns the Smith normal form of a given matrix.+smithNF :: Integral a => M.Matrix a -> M.Matrix a+smithNF m = (\diags -> extraZeros (M.diagonalList (length diags) 0 diags)) $ map abs $ rectifyDiagonal $ diagonalize m+  where extraZeros d = if M.nrows m > M.ncols m+                       then (M.<->) d (M.zero ((M.nrows m) - (M.ncols m)) (M.ncols m))+                       else (M.<|>) d (M.zero (M.nrows m) ((M.ncols m) - (M.nrows m)))++-- | Given a diagonal matrix, outputs a list \([d_1,..,d_n]\) that satisfies +--  \(d_k \mid d_{k+1}\) and represents the diagonal entries.+-- Assumes input is a diagonal matrix (not checked).+rectifyDiagonal :: Integral a => M.Matrix a -> [a]+rectifyDiagonal diagonalMatrix+  | length diags <= 1 = diags+  | allDivisible = diags+  | otherwise = rectifyDiagonal $ diagonalize $ (\op -> op (M.diagonalList (length diags) 0 diags)) $ head $ map modifier $ filter (\(b, _) -> b == False) $ zip divisibles divIndices +  where diags = V.toList (M.getDiag diagonalMatrix)+        divIndices = zip [0..(length diags) - 2] [1..(length diags) - 1]+        divPairs = map (\(k, k') -> (diags !! k, diags !! k')) $ divIndices+        divisibles = map (\(d, d') -> d `divides` d') divPairs+        allDivisible = and $ divisibles+        modifier (b, (i, j)) = if b then id else M.setElem (diags !! j) (i+2, i+1)+++-- | Given a matrix, returns a diagonal matrix obtained by applying +-- elementary row and column operations, but which does not necessarily satisfy the divisibility property+diagonalize :: Integral a => M.Matrix a -> M.Matrix a+diagonalize = diagonalizer 1++diagonalizer :: Integral a => Int -> M.Matrix a -> M.Matrix a+diagonalizer rowIndex m +  | rowIndex > M.nrows m = m+  | isDiagonalMatrix m = m+  | pivotPosition == (-1, -1) = m -- means there's no more cols+  | hasNonzeroAmongZeros areZeroCols = diagonalizer rowIndex $ (\(zeroColIndex, _) -> M.switchCols zeroColIndex (zeroColIndex+1) m) $ head $ dropWhile (\(_, b) -> not b) $ zip [1..] areZeroCols +  | isZero restOfRow && isZero restOfCol = diagonalizer (rowIndex + 1) m+  | otherwise = diagonalizer rowIndex (clearPivotRow pivotPosition $ improvePivot pivotPosition (i + 1) mat)+  where areZeroCols = map isZero $ cols m+        (pivotPosition, mat) = choosePivot rowIndex m +        (i, j) = pivotPosition+        restOfRow = M.rowVector $ V.drop j $ M.getRow i m+        restOfCol = M.colVector $ V.drop i $ M.getCol j m++hasNonzeroAmongZeros :: [Bool] -> Bool+hasNonzeroAmongZeros (a:b:xs) = (a && (not b)) || (hasNonzeroAmongZeros (b:xs))+hasNonzeroAmongZeros _ = False++-- | Returns whether a matrix (not necessarily square) is diagonal+isDiagonalMatrix :: (Num a, Eq a) => M.Matrix a -> Bool+isDiagonalMatrix m = and $ map (== 0) [M.getElem i j m | i <- [1..M.nrows m], j <- [1..M.ncols m], i /= j]++clearPivotRow :: Integral a => (Int, Int) -> M.Matrix a -> M.Matrix a+clearPivotRow (t, jt) m = M.transpose $ improvePivot (jt, t) (jt+1) (M.transpose m)++-- INVARIANT: all operations do not change the absolute value of the determinant+-- i.e. (I) scale a row by a unit (which are just 1 and -1 over the integers)+--     (II) switch two rows +--    (III) add a multiple of one row to another+-- This does most of the real work.+improvePivot :: Integral a => (Int, Int) -> Int -> M.Matrix a -> M.Matrix a+improvePivot (t, jt) rowIndex m+  | pivot == 0 = error "Zero pivot entry."  +  | rowIndex > M.nrows m = m+  | pivot `divides` nextEntry = improvePivot (t, jt) (rowIndex + 1) (M.mapRow (\j elt -> elt + ((M.getElem 1 j pivotRow)*(-1)*(nextEntry `div` pivot))) rowIndex m)+  | otherwise = improvePivot (t, jt) rowIndex (M.switchRows t rowIndex (M.mapRow (\j elt -> elt - (nextEntry `div` pivot)*(M.getElem 1 j pivotRow)) rowIndex m))        +  -- idea: nextEntry = pivot * q + r+  where pivot = M.getElem t jt m+        nextEntry = M.getElem rowIndex jt m   +        pivotRow = row t m++-- | Returns whether a `divides` b := \(a \mid b\) +-- and handles the special case of \(a=0\)+divides :: Integral a => a -> a -> Bool+divides 0 0 = True+divides 0 _ = False+divides a b = b `mod` a == 0++---------------------+-- PIVOT-SELECTION --+---------------------+choosePivot :: Integral a => Int -> M.Matrix a -> ((Int, Int), M.Matrix a)+choosePivot rowIndex m+  | pivotPosition == Nothing = ((-1, -1), m)    +  | (M.getElem t jt m) /= 0 = ((t, jt), m)+  | otherwise = ((t, jt), makePivotNonzero (t, jt) m)+  where pivotPosition = nextPivotPosition rowIndex m+        (t, jt) = fromJust pivotPosition++-- ASSUME: there is a nonzero entry in column jt, not at (t, jt), but at jt' > jt+makePivotNonzero :: Integral a => (Int, Int) -> M.Matrix a -> M.Matrix a+makePivotNonzero (t, jt) m = M.switchRows t nonzeroIndex m+  where nonzeroIndex = head [i | (i, entry) <- zip [1..] (M.toList (col jt m)), entry /= 0, i > t]++-- ASSUMES: the rows with index < rowIndex only have a single nonzero entry, +--   which occurs in columnIndex < rowIndex +nextPivotPosition :: Integral a => Int -> M.Matrix a -> Maybe (Int, Int)+nextPivotPosition rowIndex m +  | null nonzeroColIndices = Nothing+  | otherwise = Just (rowIndex, (head nonzeroColIndices))+  where nonzeroColIndices = [j | (j, column) <- zip [1..] (cols m), not (isZero column), j >= rowIndex]++isZero :: (Num a, Eq a) => M.Matrix a -> Bool+isZero m = m == M.zero (M.nrows m) (M.ncols m)++---------------------------+-- MATRIX HELPER METHODS --+---------------------------+-- return a list of rows, as 1 x n matrices+rows :: M.Matrix a -> [M.Matrix a]+rows m = map (\k -> row k m) [1..(M.nrows m)]++-- return a list of cols, as n x 1 matrices+cols :: M.Matrix a -> [M.Matrix a]+cols m = map (\k -> col k m) [1..(M.ncols m)]++-- get a row, represented as an 1 x n matrix+row :: Int -> M.Matrix a -> M.Matrix a+row k = (M.rowVector . M.getRow k) ++-- get a column, represented as a n x 1 matrix+col :: Int -> M.Matrix a -> M.Matrix a+col k = (M.colVector . M.getCol k)