New Julia package offering Smith and Hermite normal forms

[brainandforce]

I'm currently working on a package, NormalForms.jl, which intends to make calculations of the Hermite and Smith normal forms available to any Julia user. Honestly, I wish this functionality was in the standard library - but since it isn't, I figured it would be good to combine the functionality of two existing packages (SmithNormalForm.jl and HermiteNormalForm.jl) as these matrix decompositions are both used heavily in crystallography.

The code isn't completely correct yet for the Smith normal form (in particular, I haven't implemented the divisibility checks for the diagonal elements, and I'm dealing with a bug that gets these calculations stuck in a while loop), but I hope having a few extra eyes on it can sort out the issues.

(Disclaimer: I have almost no experience with implementing algorithms, especially for matrix decompositions...I'm in uncharted territory)

[brainandforce]

Quick update: I just solved the issues with the Smith normal form calculations. Barring any major issues, I should have the package registered in the next few days.

[thchr]

This would be really nice: I'd be glad to drop the weird "copy-vendoring" that I currently resort to with SmithNormalForm.jl.

How does this compare performance-wise with SmithNormalForm.jl? This is not critical, but would be nice if it could match that.

[brainandforce]

Well, it looks like performance needs work:

julia> M = [-1 1 1; 1 -1 1; 1 1 -1] * diagm([1, 1, 2])
3×3 Matrix{Int64}:
 -1   1   2
  1  -1   2
  1   1  -2

julia> @benchmark SmithNormalForm.smith(M)
BenchmarkTools.Trial: 10000 samples with 11 evaluations.
 Range (min … max):  990.455 ns … 721.662 μs  ┊ GC (min … max): 0.00% … 99.06%
 Time  (median):       1.163 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.527 μs ±   8.443 μs  ┊ GC (mean ± σ):  7.48% ±  1.40%

  █▆▅▆▄▄▃▂▁▂▃▃▄▂▁ ▁ ▂▂▁▁▁                                       ▂
  ████████████████████████▇▇▆▇▆▆▆▆▅▆▅▆▅▆▆▅▅▅▅▅▄▅▄▄▃▄▅▅▄▄▅▄▄▅▄▅▅ █
  990 ns        Histogram: log(frequency) by time       4.94 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 13.

julia> @benchmark NormalForms.snf(M)
BenchmarkTools.Trial: 10000 samples with 3 evaluations.
 Range (min … max):   8.084 μs …  2.214 ms  ┊ GC (min … max):  0.00% … 97.50%
 Time  (median):     10.111 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):   14.288 μs ± 51.645 μs  ┊ GC (mean ± σ):  12.07% ±  3.41%

  ▇█▅▃▄▂▃▆▅▄▅▄▃▄▂▄▃▂▁▁▁▂▁▂▁▁▁                                 ▂
  ████████████████████████████████▇█▆█▇▇▇▇▇▆▆▆▆▆▄▆▆▆▅▄▅▅▄▅▅▄▅ █
  8.08 μs      Histogram: log(frequency) by time      38.5 μs <

 Memory estimate: 22.53 KiB, allocs estimate: 286.

[brainandforce]

Quick update: I have registered the package and solved the performance issues listed above. However, I think the algorithms that are being used may have some issues. They seem to work for 3×3 matrices I've been testing, but some other shapes/sizes have run into problems.

[brainandforce]

Another update: I've sorted out issues that have come up in my own use cases (generating atomic positions when building a supercell), which primarily involve 3x3 matrices. There are still some matrices which show failures in unimodularity when calculating the Hermite or Smith normal forms, but overall I've found the package to work in the 3D case, which is what I assume you're using it for.

[brainandforce]

One more update: I believe I've solved all outstanding issues regarding failed Smith/Hermite normal form factorizations for 3x3 matrices (thanks @glwhart). The unimodularity failures I've run into seem to be restricted to cases where the matrices have large enough (pseudo)determinants to overflow.

[thchr]

Hi @brainandforce; thanks for your work on a proper, maintained Smith/Hermite normal form package.

Unfortunately, I've come to realize that some of the work I'm doing that uses a Smith normal form depends on the the precise form that is obtained; basically, to keep my head sane, I need to not change what smith returns, even if it is another valid decomposition, at least for a while.

The bits that I use the Smith normal form for here are strictly for quite big matrices; maybe on the order of 10x10 and bigger, generally.

[brainandforce]

By precise form, do you mean the forms of the unimodular matrix factors? I've actually noticed this issue in testing, but the Smith normal forms themselves are correct. It appears the factors aren't uniquely determined.

[thchr]

Indeed, yeah: to stick with whatever conventions I've implicitly been using, I need the $T$ and $S$ matrices to keep having the same form they did with the old Smith normal form code. They are indeed not uniquely determined.