wish list: smith normal form transformation

[argriffing]

Flint has fmpz_mat_hnf_transform and fmpz_mat_snf but I'm curious about a hypothetical fmpz_mat_snf_transform. According to http://www.informatik.uni-kiel.de/~gej/publ/transsmith.pdf "There are many algorithms for efficiently computing the SNF. Some of these algorithms also compute the transformation matrices."

[alexjbest]

Hi, I definitely agree that snf_transform would be a nice addition and I think I remember two ways of doing this that are not too complicated.

First let me explain why there is hnf_transform but no snf_transform: there is actually a very easy way computing the transform along with the HNF (and this is what is done if you look at hnf_transform) you put an identity matrix on the right hand side of your matrix and then use whatever HNF method you like. Your HNF function will automatically record what it is doing on this identity matrix, so that your resulting block matrix can be split into HNF and transform matrix. (Unfortunately as far as I can tell no such simple trick would work for SNF as sticking an identity matrix on anywhere makes your SNF the identity automatically.)

As for straightforwardish ways of actually finding the SNF + transform in Flint: The first is to make transform versions of the existing SNF algorithms (or possible only the most general one snf_kannan_bachem) which every time an operation is performed on the input matrix does the same thing to a transform matrix (one for column ops and one for row ops). One issue with this is having potentially two copies of broadly the same code, one with these operations in and one without which might lead to problems with maintaining both sets of code, though this could probably be solved by having a with_transform flag or some such with only a slight performance cost to the non transformed function (or maybe there is some clever compile time stuff that can be done here to produce both functions from the same code, this seems possible but I don't know how right now), anyway I think this is why I didn't add a transform version of snf_kannan_bachem. There are also probably issues with this method when using modular SNF algorithms (like snf_iliopoulos), but I haven't thought about it recently so couldn't really say.

The other way is kind of funny and I think probably very inefficient, but if the matrices involved are small or have some other structure that might make their SNFs quick to compute anyway then it could be workable and quite simple to implement. The method is to repeatedly call hnf_transform and then transpose the resulting HNF matrix to get a new matrix to hnf_transform again, repeating until the matrix ends up in SNF, the transform matrices found along the way can be multiplied to get the row transform (from the odd iteration transformations) and the column transform (from the even iterations).

There are undoubtedly far better ways of doing this (i.e. faster and/or returning smaller transformation matrices), but the two approaches above would be easy enough to write given the existing code in Flint I think.

Did you have a specific use case in mind (size of matrix, entries, special structure) and do you think either of these approaches could work there?

Alex

[argriffing]

I like the idea of implementing a fmpz_mat_snf_transform by copying the Kannan-Bachem code and adding some bookkeeping. I'm not interested in a modular variant, and I think Flint is OK with the 5th rule of good C programming.

Did you have a specific use case in mind (size of matrix, entries, special structure)

Honestly I don't know what I'm doing. I have C code for a combinatorial optimization problem which involves some recursions and log/exp, but with certain rational parameters the computer jams up because even with arbitrary precision math it doesn't understand the exact cancellations. My idea is to use a change of basis so that detecting the cancellations becomes easier or trivial, and I think that doing this efficiently somehow involves Smith or Hermite transforms.

[argriffing]

I ended up using only fmpz_mat_hnf_transform instead. But I'm explicitly inverting the U matrix afterwards which is slow and makes me feel stupid. I notice that GAP has options to build the inverse transform matrix while the normal form is being created. Maybe this could be added to flint? Or maybe the user should augment their input matrices by putting an identity matrix somewhere?

[alexjbest]

I finished testing the snf_transform code I wrote, and while it was a little trickier than I anticipated it all looks to be working (with the above caveats about speed and size of transformation matrices mentioned above), the code is on my github repo, specifically the functions fmpz_mat_snf_transform, fmpz_mat_snf_kannan_bachem_trans and fmpz_mat_snf_diagonal_trans.

That's interesting, which gap function is it? I can't see a particularly nice way of doing it as you go, really the only thing I can think of at the moment is using a similar method to as mentioned above and recording the inverse transformations, once again this is not too complicated to write for a straightforward hnf implementation, but considerably more difficult for a complicated one like Pernet-Stein and probably not very efficient. As your transformation U satisfies U^-1 H = A you could maybe also solve for U^-1 quickly afterwards (as H is upper triangular, but also possibly singular), but I don't know if any of this is really going to be faster than just inverting the U hnf_transform gives you.

[argriffing]

That's interesting, which gap function is it?

I did not try using it, but I found a function named HermiteNormalFormIntegerMatInverseTransforms mentioned in the docs at http://www.math.rwth-aachen.de/~Greg.Gamble/gap4r3/doc/htm/ref/CHAP025.htm#SSEC001.7

[alexjbest]

Hmm well I can't see that function in the doc on the gap site, so I guess they might have removed it? In any case I don't think there is a clever way of doing this for fast Hermite form algorithms so the best way in Flint is likely to simply invert the transformation matrix.

[antieau]

What is the status of this issue? It seems like this functionality now exists as the function fmpz_mat_snf_transform in src/acb_theta/sp2gz_decompose.c. I cannot find that in the documentation however. I would like to use this functionality (in Python-FLINT too).

[fredrik-johansson]

@j-kieffer Any reason for not making fmpz_mat_snf_transform public? How does the algorithm compare to our current fmpz_mat_snf?

[j-kieffer]

I didn't realize there was an open issue about this ! I just needed a quick hack as in the acb_theta module, the matrices will never be larger than 20x20. I'm using successive HNF transforms as above:

The other way is kind of funny and I think probably very inefficient, but if the matrices involved are small or have some other structure that might make their SNFs quick to compute anyway then it could be workable and quite simple to implement. The method is to repeatedly call hnf_transform and then transpose the resulting HNF matrix to get a new matrix to hnf_transform again, repeating until the matrix ends up in SNF, the transform matrices found along the way can be multiplied to get the row transform (from the odd iteration transformations) and the column transform (from the even iterations).

I would guess that the "right" way to make such a function public would be to adapt the code of fmpz_mat_snf instead, which should be much faster.

[antieau]

Thanks for the comments! I expect to have a couple of students programming for me over the summer. I will ask one of them to write a version of this adapting fmpz_mat_snf and have them submit a pull request.

[alexjbest]

The old code I wrote earlier in the thread is still available on a branch https://github.com/flintlib/flint/compare/main...alexjbest:flint2:alex/old?expand=1, I'd be curious if it still builds and how the runtime compares with Jean's. Your students should feel free to reuse that also if it is helpful of course!