12
$\begingroup$

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?

$\endgroup$
11
$\begingroup$

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]

Mathematica graphics

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

$\endgroup$
  • 1
    $\begingroup$ Coded in 1994 and it still works in Mathematica 9, that's how it should be. $\endgroup$ – nilo de roock Jan 22 '14 at 18:59
  • $\begingroup$ @ndroock1. I agree. $\endgroup$ – RunnyKine Jan 22 '14 at 19:01
  • $\begingroup$ Is there a package for Mathematica that computes the SNF (actually just the invariant factors) of a sparse matrix (Sparse Array)? $\endgroup$ – Leon Apr 7 '16 at 18:38
  • $\begingroup$ @Leon. Have you checked if the new built-in function gives you that? $\endgroup$ – RunnyKine Apr 7 '16 at 18:44
  • $\begingroup$ I have Mathematica 10.0 which doesn't have this built in yet. In reference.wolfram.com/language/ref/SmithDecomposition.html there is no mention of sparse matrices, and there is no InvariantFactors command. $\endgroup$ – Leon Apr 7 '16 at 19:54
10
$\begingroup$

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]

Mathematica graphics

Where the second element is the Smith normal form.

$\endgroup$
  • $\begingroup$ In case you want to use it in a script, SmithDecomposition[mat][[2]] selects just the SNF. $\endgroup$ – 黄雨伞 Apr 10 '18 at 17:26
8
$\begingroup$

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 *)
$\endgroup$
  • 1
    $\begingroup$ I really wonder why Mma does not have this builtin. $\endgroup$ – Mariano Suárez-Álvarez Jun 28 '14 at 19:18
  • 1
    $\begingroup$ @Mariano Suárez-Alvarez About all I can say is I wonder about it myself. $\endgroup$ – Daniel Lichtblau Jun 28 '14 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.