Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language
As of Version 14.3, PolynomialSmithDecomposition is a built-in symbol.

Function Repository Resource:

 PolynomialSmithDecomposition

Source Notebook

Compute the Smith decomposition of a matrix of univariate polynomials

Contributed by: Daniel Lichtblau

ResourceFunction["PolynomialSmithDecomposition"][mat]

computes the Smith decomposition of the matrix mat of univariate polynomials.

ResourceFunction["PolynomialSmithDecomposition"][mat,x]

computes the Smith decomposition for a matrix of polynomials in the variable x.

Details and Options

The result is given in the form {u,s,v}, where u and v are unimodular matrices, s is a diagonal matrix and u.mat.v⩵s. Moreover, each diagonal element of s divides all subsequent diagonal elements.
A unimodular matrix over a ring of univariate polynomials is a matrix with nonzero determinant lying in the coefficient field (that is, constant).
The Smith decomposition is computed using iterations of the Hermite decomposition, as implemented in the resource function PolynomialHermiteDecomposition.
Multivariate polynomials are regarded as univariate in the specified variable, with all others treated as symbolic coefficients.
ResourceFunction["PolynomialSmithDecomposition"] takes Method and Modulus as options. These are the same as for PolynomialHermiteDecomposition.
For specifics of the Method option, see documentation for PolynomialHermiteDecomposition. The only uses in ResourceFunction["PolynomialSmithDecomposition"] are to pass it to the Hermite computations.
ResourceFunction["PolynomialSmithDecomposition"] is intended for matrices of polynomials in a single variable, with all coefficients being exact numbers.

Examples

Basic Examples (4) 

Compute the Smith decomposition of a 2×3 matrix of low degree polynomials:

In[1]:=
mat = {{x^3 - 11, x^2 - 2 x + 5, x^3 - 2 x^2 + 3 x - 7}, {x^2 + 4 x + 9, x^2 - 19, x^3 + 5 x - 6}}; {u, s, v} = ResourceFunction["PolynomialSmithDecomposition"][mat]
Out[2]=

Check the matrix equation:

In[3]:=
Expand[u . mat . v - s]
Out[3]=

Check that u and v are unimodular; that is, they are constants:

In[4]:=
{Det[u], Det[v]}
Out[4]=

Check that the diagonal entries of s divide into subsequent entries:

In[5]:=
Apply[Divisible, Reverse[#]] & /@ Partition[Diagonal[s], 2, 1]
Out[5]=

Options (6) 

Modulus (6) 

Generate and compute the Smith decomposition for an 8×12 matrix of random degree-5 polynomials with coefficients between -10 and 10:

In[6]:=
randomPoly[deg_, max_, x_] := RandomInteger[{-max, max}, deg + 1] . x^Range[0, deg] randomMatrix[m_, n_, deg_, max_, x_] := Table[randomPoly[deg, max, x], {m}, {n}] SeedRandom[1111]; bigmat = randomMatrix[8, 12, 5, 10, x]; Timing[{u, s, v} = ResourceFunction["PolynomialSmithDecomposition"][bigmat]] // First
Out[10]=

As is often the case in symbolic computation, speed improves when one works over a prime field of modest size:

In[11]:=
Timing[{umod, smod, vmod} = ResourceFunction["PolynomialSmithDecomposition"][bigmat, Modulus -> 103]] // First
Out[11]=

Check that the Smith form is a diagonal matrix:

In[12]:=
DiagonalMatrixQ[smod]
Out[12]=

Check that the decomposition satisfies the necessary matrix identity:

In[13]:=
PolynomialMod[Expand[umod . bigmat . vmod - smod, Modulus -> 103], Modulus -> 103]
Out[13]=

Check that the left factor is unimodular:

In[14]:=
Det[umod]
Out[14]=

The determinant of the right factor cannot easily be computed due to expression swell. A probabilistic test that it is unimodular (that is, a constant) is to check that it gives the same value with multiple substitutions of the polynomial variable:

In[15]:=
Table[Mod[Det[vmod /. x -> j], 103], {j, 1, 25, 3}]
Out[15]=

Neat Examples (2) 

Use the Smith decomposition to compute the minimal polynomial of a matrix:

In[16]:=
mat = \!\(\* TagBox[ RowBox[{"(", "", GridBox[{ {"4", RowBox[{"-", "1"}], RowBox[{"-", "4"}], "9", RowBox[{"-", "2"}], RowBox[{"-", "3"}], "6", "1"}, { RowBox[{"-", "5"}], "1", "0", RowBox[{"-", "2"}], "2", "2", RowBox[{"-", "3"}], "0"}, {"8", RowBox[{"-", "2"}], RowBox[{"-", "5"}], "14", RowBox[{"-", "4"}], RowBox[{"-", "6"}], "10", "2"}, {"4", RowBox[{"-", "1"}], RowBox[{"-", "4"}], "8", RowBox[{"-", "2"}], RowBox[{"-", "3"}], "5", "1"}, {"7", RowBox[{"-", "1"}], RowBox[{"-", "4"}], "15", RowBox[{"-", "4"}], RowBox[{"-", "6"}], "10", "1"}, {"5", RowBox[{"-", "1"}], RowBox[{"-", "4"}], "7", RowBox[{"-", "3"}], RowBox[{"-", "3"}], "4", "1"}, { RowBox[{"-", "3"}], "1", "4", RowBox[{"-", "7"}], "1", "2", RowBox[{"-", "5"}], RowBox[{"-", "1"}]}, {"15", RowBox[{"-", "3"}], RowBox[{"-", "8"}], "27", RowBox[{"-", "8"}], RowBox[{"-", "11"}], "20", "4"} }, GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}}, GridBoxSpacings->{"Columns" -> { Offset[0.27999999999999997`], { Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> { Offset[0.2], { Offset[0.4]}, Offset[0.2]}}], "", ")"}], Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\);
In[17]:=
ResourceFunction["PolynomialSmithDecomposition"][ IdentityMatrix[Length[mat]] x - mat, x][[2, -1, -1]]
Out[17]=

The result from PolynomialSmithDecomposition is a scalar multiple of the result of the resource function MatrixMinimalPolynomial:

In[18]:=
ResourceFunction["MatrixMinimalPolynomial"][mat, x] // Expand
Out[18]=

For a square matrix mat, the nontrivial diagonal elements in the Smith normal form for mat-x ℑ are the invariant factors of the characteristic polynomial of mat (here ℑ is the identity matrix of same dimension as mat). We can use this to compute the rational canonical form.

We start with code to create a companion matrix from a polynomial:

In[19]:=
companionMatrix[poly_, x_] := Module[ {clist = CoefficientList[poly, x], len = Exponent[poly, x], cmat}, cmat = ConstantArray[0, {len, len}]; Do[cmat[[i, i - 1]] = 1, {i, 2, len}]; cmat[[All, -1]] = -Most[clist]/Last[clist]; cmat ]

Work with an upper triangular matrix having a nontrivial minimal polynomial and hence nontrivial rational canonical form:

In[20]:=
SeedRandom[1234]; mat = DiagonalMatrix[{1, 1, 2, 2, 2, 2, 3, 3, 3}]; Do[mat[[i, j]] = RandomInteger[], {i, 1, 4}, {j, i + 1, 5}]; MatrixForm[mat]
Out[23]=

Find the diagonal for the Smith normal form of mat-x:

In[24]:=
{uu, ss, vv} = ResourceFunction["PolynomialSmithDecomposition"][ mat - x*IdentityMatrix[Length[mat]], x]; diag = Diagonal[ss]
Out[25]=

Compute the submatrices for the diagonal block:

In[26]:=
compmats = Map[companionMatrix[#, x] &, diag]
Out[26]=

Use SparseArray and Band to reconstruct the companion matrix for mat as a diagonal block matrix:

In[27]:=
compmat = Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]]
Out[27]=

Version History

  • 1.0.5 – 30 January 2025
  • 1.0.4 – 17 January 2025
  • 1.0.3 – 18 October 2022
  • 1.0.2 – 19 January 2022

Related Resources

Author Notes

PolynomialSmithDecomposition inherits any weaknesses found in the resource function PolynomialHermiteDecomposition. Since the Smith decomposition uses iterations of the Hermite decomposition, such issues could be compounded. It might not give good results, for example, in the case where input polynomials have approximate coefficients.

License Information