# Smith Normal Form?

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?

``````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

2 Likes

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.

``````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

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