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?

nilo de roock

Posted 2014-01-22T16:39:50.207

Reputation: 8 701

Answers

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]

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

RunnyKine

Posted 2014-01-22T16:39:50.207

Reputation: 32 260

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]

Mathematica graphics

Where the second element is the Smith normal form.

RunnyKine

Posted 2014-01-22T16:39:50.207

Reputation: 32 260

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 *)

Daniel Lichtblau

Posted 2014-01-22T16:39:50.207

Reputation: 52 368

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