## Has Mathematica a function to compute the Smith Normal Form?

13

5

The Smith normal form is a matrix that can be calculated for any matrix (not necessarily square) with integer entries. See Wikipedia for a more elaborate description. Has Mathematica a function to compute the Smith Normal Form ? If not, is there an easy way to compute it in Mathematica?

12

I don't think Mathematica has a built-in function to compute the Smith Normal Form, but here are a couple packages by David Jabon that can do this for you: Smith Normal Form

Usage

Needs@"IntegerSmithNormalForm"


Here is an integer matrix mat

mat = {{1, 2, 3}, {-2, 3, 1}, {3, 2, 1}}


We compute it's SmithForm

smith = SmithForm[mat]


One can even compute the ExtendedSmithForm

(* "ExtendedSmithForm[A] gives, for an integral
matrix A, {D, {P,Q}} where
D is in Smith normal form and P and Q are
matrices such that P A Q = D. "  *)

extsmith = ExtendedSmithForm[mat]


Gives:

{{{1, 0, 0}, {0, 1, 0}, {0, 0, 28}},
{{{1, 0, 0}, {4, -1, -2}, {13, -4, -7}}, {{1, -2, 15}, {0, 1, -9}, {0, 0, 1}}}}


We can check that indeed P A Q == D

extsmith[[2, 1]].mat.extsmith[[2, 2]] == smith


True

1Coded in 1994 and it still works in Mathematica 9, that's how it should be. – nilo de roock – 2014-01-22T18:59:58.663

@ndroock1. I agree. – RunnyKine – 2014-01-22T19:01:35.950

Is there a package for Mathematica that computes the SNF (actually just the invariant factors) of a sparse matrix (Sparse Array)? – Leo – 2016-04-07T18:38:30.667

@Leon. Have you checked if the new built-in function gives you that? – RunnyKine – 2016-04-07T18:44:16.910

I have Mathematica 10.0 which doesn't have this built in yet. In http://reference.wolfram.com/language/ref/SmithDecomposition.html there is no mention of sparse matrices, and there is no InvariantFactors command.

– Leo – 2016-04-07T19:54:26.297

11

Just to separate this from a package-based answer. In Mathematica 10.2, you can now do this with the built-in function SmithDecomposition. So using the same matrix from my previous answer:

mat = {{1, 2, 3}, {-2, 3, 1}, {3, 2, 1}};
MatrixForm /@ SmithDecomposition[mat]


Where the second element is the Smith normal form.

In case you want to use it in a script, SmithDecomposition[mat][[2]] selects just the SNF. – 黄雨伞 – 2018-04-10T17:26:58.107

8

Other than indenting the code blocks, what follows comes verbatim from an old file of mine. I have not tested it in the current Mathematica, though I have no reason to think it might fail.

diagonalQ[mat_?MatrixQ] := With[
{posns=Flatten[Map[Position[#, _?(#!=0&)]&, mat]]},
Length[Union[posns]] == Length[posns]
]

diagonalize[mat_?MatrixQ] := Module[
{hnf=mat, umat=IdentityMatrix[Length[mat]],
vmat=IdentityMatrix[Length[mat[[1]]]], tmpu, tmpv},
While[Not[diagonalQ[hnf]],
{tmpu, hnf} = HermiteDecomposition[hnf];
umat = tmpu . umat;
{tmpv, hnf} = HermiteDecomposition[Transpose[hnf]];
vmat = vmat . Transpose[tmpv];
hnf = Transpose[hnf];
];
{umat, hnf, vmat}
]

divides[a_, b_] := GCD[a,b]===a

smithNormalForm[mat_?MatrixQ] := Module[
{uu, dd, vv, diags, gcd, col=0, dim, tmpu, tmpv},
{uu, dd, vv} = diagonalize[mat];
diags = Select[Flatten[dd], #!=0&];
dim = Length[diags];
While[col+1<dim,
col++;
If [divides[diags[[col]],GCD[Apply[Sequence,Drop[diags,col]]]],
Continue[]];
vv = Transpose[vv];
Do [
dd[[j,col]] = diags[[j]];
vv[[col]] += vv[[j]],
{j,col+1,dim}];
vv = Transpose[vv];
{tmpu, dd, tmpv} = diagonalize[dd];
uu = tmpu  . uu;
vv = vv . tmpv;
diags = Select[Flatten[dd], #!=0&];
];
{uu, dd, vv}
]


For examples:

SeedRandom[1111];
m = 4; n = 6;
uppermat[m_] :=
Table[If[i==j,1,If[i>j,0,Random[Integer,{-10,10}]]], {i,m}, {j,m}]
umat[m_] := uppermat[m] . Transpose[uppermat[m]]
uu = umat[m];
vv = umat[n];
hh = {{1,0,0,0,0,0}, {0,3,0,0,0,0}, {0,0,0,18,0,0}, {0,0,0,0,0,90}};
testmat = uu . hh . vv;


To change this just rig hh so it is appropriate for the dimensions.
To get nontrivial examples you might want to have some rows and/or
columns zeroed, and you'd want nontrivial gcds for the nonzero
elements.

{newu, newh, newv} = smithNormalForm[testmat]
Max[Abs[newu . testmat . newv - newh]] (* want zero *)

{newu, newh, newv} = smithNormalForm[Transpose[testmat]]
Max[Abs[newu . Transpose[testmat] . newv - newh]] (* want zero *)
`

1I really wonder why Mma does not have this builtin. – Mariano Suárez-Álvarez – 2014-06-28T19:18:02.027

2@Mariano Suárez-Alvarez About all I can say is I wonder about it myself. – Daniel Lichtblau – 2014-06-28T20:05:58.123