Smith Normal Form?

Specific Domains Numerics
1

Does anyone have knowledge of packages for Smith Normal Form (e.g., SmithNormalForm.jl)?

I did a simple test of SmithNormalForm.jl:

Apart from the example in the GitHub page,… function smith seems to work:

  • rectangular matrices
  • matrices with integer and floating point numbers

Long time ago, I used Maple to find the Smith form of rectangular linear pencils , and am looking for this possibility in Julia. I specifically tried with:

  1. apply smith to a matrix of rational numbers – this did not work (a little surprising, since it worked with floats…)
  2. apply smith to a matrix s*I - A where s is defined as @variables s using the Symbolics.jl package

Neither 1 nor 2 work.

The Smith form of a linear matrix pencil is of some interest in control theory.

→ is there a way to update the algorithms such that both rational numbers and Symbolics.jl variables work?

2 
julia> using Nemo

Welcome to Nemo version 0.37.2

Nemo comes with absolutely no warranty whatsoever


julia> const nums = -50:50
-50:50

julia> const dens = ((-50:-1)..., (1:50)...);

julia> r() = rand(nums)//rand(dens)
r (generic function with 1 method)

julia> m = [r() for _ ∈ 1:3, _ ∈ 1:5]
3×5 Matrix{Rational{Int64}}:
  27//4   27//10   -3//5    -23      43//29
 -29//37  28//39  -36//49  -23//25  -13//22
  -9//4     8      -9//11    7//6   -35//29

julia> a = matrix(QQ, m)
[  27//4   27//10     -3//5       -23    43//29]
[-29//37   28//39   -36//49   -23//25   -13//22]
[  -9//4        8    -9//11      7//6   -35//29]

julia> snf(a)
[1   0   0   0   0]
[0   1   0   0   0]
[0   0   1   0   0]

julia> _, t, u = snf_with_transform(a);

julia> t * a * u
[1   0   0   0   0]
[0   1   0   0   0]
[0   0   1   0   0]

Nemo also has a sister package, AbstractAlgebra, which is supposed to have a lot of the same functionality as Nemo, but implemented in pure Julia. The above is supposed to work even if we replaced using Nemo with using AbstractAlgebra, however it seems there’s a bug or missing functionality: Smith normal form seems broken for matrices of rational numbers · Issue #1477 · Nemocas/AbstractAlgebra.jl · GitHub

3 Likes
3

For a matrix with rational entries (or entries in an arbitrary field), the Smith normal form is rather simple. It is just the matrix here: Matrix equivalence - Wikipedia. It is uniquely determined by the rank of the matrix.

For task 2):

julia> using Nemo

julia> Qx, x = QQ["x"];

julia> M = matrix(Qx, [1 2; 4 5])
[1   2]
[4   5]

julia> snf(x*identity_matrix(Qx, 2) - M)
[1               0]
[0   x^2 - 6*x - 3]

julia> elementary_divisors(x * identity_matrix(Qx, 2) - M) # also known as invariant factors sometimes
2-element Vector{QQPolyRingElem}:
 1
 x^2 - 6*x - 3
3 Likes
4

I have a package, NormalForms.jl, providing Smith and Hermite normal forms as subtypes of LinearAlgebra.Factorization, and it’s under active development. As of now, it is written specifically for integer matrix operations. However, I’m sure 1) could be implemented in the package pretty easily. As for 2), I’ve never used Symbolics.jl, but if I can find maintainers who are familiar with it, I’d also be happy to integrate this functionality too.

I will open issues relating to this on the repo and you are welcome to file a PR if you are interested in contributing.

1 Like