possible memory bug in fmpz_mat_snf()

[John Cremona]

I have been using the Smith Normal Form function on some integer matrices with the function fmpz_mat_snf().  The matrices are not particularly large (perhaps 1-300 rows/cols) and with small entries, and usuallt fmpx_mat_snf() works find (and quickly), but sometimes it eats up more and more RAM until it quits.

I don't think there is anything wrong in the input matrix M and by code does

fmpz_mat_t S;

fmpz_mat_init_set(S, M); // to set to the right size

fmpz_mat_snf(S, M);

so please tell me if that is incorrect usage.  I will try to prepare a test code + minimal program to replicate the error.

FLINT is great, thanks!

John

[John Cremona]

Here is a test program snf_test.cc (which is C++ as that is what I use) and two sample data files, compressed. If I compile the program and run it,

./snf_test < A.txt

works fine, but

./snf_test < B.txt

blows up.

Incidentally, I think is a bit silly for the fmpz_mat_read functions to require the fmpz_mat_t object to be initialised with the correct number of rows and columns, when the input stream is required to include those two numbers. Why not have the read function do the init itself after reading them? That is why my input files have the number of rows and columns in twice.

John

[Fredrik Johansson]

Hi John,

Thanks for the bug report.

I put a printout in fmpz_mat_snf_kannan_bachem to see what is going on.

For matrix A the intermediate entries remain smallish most of the way through the algorithm, just ballooning to a few thousand bits in the last few rounds.

For matrix B they suddenly explode around iteration 600/908, growing to more than a million bits. I guess this is a pathological input for the algorithm and that we just need a better SNF algorithm here.

Fredrik

--

---
You received this message because you are subscribed to the Google Groups "flint-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAD0p0K6t5BOmD%3DQ4u_O6pQ8F%2BRYiRHkpTibrBmm3325ECcUE1w%40mail.gmail.com.

[Claus Fieker]

Just for fun:

fmpz_mat_hnf_classical
terminates in <2 sec....

Greetings
Claus

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAD0p0K6t5BOmD%3DQ4u_O6pQ8F%2BRYiRHkpTibrBmm3325ECcUE1w%40mail.gmail.com.

> #include <iostream>
> #include <flint/fmpz.h>
> #include <flint/fmpz_mat.h>
>
> int main ()
> {
> fmpz_mat_t M, S;
> int r, c;
> std::cin >> r >> c;
> fmpz_mat_init(M,r,c);
> fmpz_mat_read(M);
> std::cout << "Input matrix (" << r<<"x"<<c<<"):\n";
> fmpz_mat_print_pretty(M);
> std::cout << "\n";

> fmpz_mat_init_set(S, M); // to set to the right size
> fmpz_mat_snf(S, M);

> std::cout << "SNF matrix:\n";
> fmpz_mat_print_pretty(S);
> std::cout << "\n";
> fmpz_mat_clear(M);
> fmpz_mat_clear(S);
> }

[Oscar Benjamin]

The matrix has some kind of banded structure. I have also seen cases
where fmpq_mat_rref_classical is faster than fmpz_mat_rref for
banded-like integer matrices. In that case it seems the banded
structure makes it faster to control coefficient growth locally with
division if it prevents global coefficient growth. Maybe something
similar applies here.

--
Oscar

> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/ZdMyNExnIOP0zJAD%40mathematik.uni-kl.de.

[Edgar Costa]

If one doesn't care about the transformation matrices, then one can apply HNF repeatedly to compute the SNF.

(some transposes are necessary along the way).

I believe this is what Magma does when one does not ask for a change of basis:
> A := Matrix(Integers(),  StringToInteger(foo[1]), StringToInteger(foo[2]), [StringToInteger(x) : x in foo[3..#foo]]) where foo := Split(Split(Read("A.txt"),"\n")[2], " ");
> B := Matrix(Integers(),  StringToInteger(foo[1]), StringToInteger(foo[2]), [StringToInteger(x) : x in foo[3..#foo]]) where foo := Split(Split(Read("B.txt"),"\n")[2], " ");
> time S := SmithForm(A);
Time: 0.250
> time S := SmithForm(B);
Time: 0.060
> time S, P, Q := SmithForm(A);
Time: 0.260
> time S, P, Q := SmithForm(B);
Time: 0.110

> GetMaximumMemoryUsage()/1024.^2;
156.875000000000000000000000000

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAHVvXxSX3A0o%2BngYA-7uH5gf0v3s5tnBqNKi10zkdQoMOB-EEA%40mail.gmail.com.

[Claus Fieker]

If one wants the and with Trafo, one asks for the hnf with Trafo....

Needs some transpose and inverse and the final sorting of the diagonal

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CA%2BiQ7x425PF0kzMS0XtUGMtNaniS%3DabK6m63-8nwOcgD5EzwJg%40mail.gmail.com.

[Claus Fieker]

On Tue, Feb 20, 2024 at 06:19:00AM +0100, Claus Fieker wrote:
> If one wants the and with Trafo, one asks for the hnf with Trafo....
> Needs some transpose and inverse and the final sorting of the diagonal

PS.: in Hecke/Oscar in julia, I cannot quite reproduce Magma's times,
but I can get the HNF in <0.2 sec, using the "same" idea that, as far as
I recall, is used in Magma: if the matrix is sparse, use a sparse HNF.
In Magma, the sparse HNF is used until the remaining part is dense, then
the dense code is called. In Hecke we do a pure, vanilla, Kannan-Bachem,
but sparse.

The Pernet-Stein algo, that is used for the HNF in Nemo (for this kind
of matrix) is supposed to be "optimal" for random input. But as
observed, John's matrix is far from random, so it might also just be the
worst input possible...
The flint Kannan-Bachem runs in 2sec. this might also be a good start
for the snf.

> >> https://groups.google.com/d/msgid/flint-devel/CAHVvXxSX3A0o%2BngYA-7uH5gf0v3s5tnBqNKi10zkdQoMOB-EEA%40mail.gmail.com
> >> .

> >>
> > --
> >
> > ---
> > You received this message because you are subscribed to the Google Groups
> > "flint-devel" group.
> > To unsubscribe from this group and stop receiving emails from it, send an
> > email to flint-devel...@googlegroups.com.
> > To view this discussion on the web, visit

> > https://groups.google.com/d/msgid/flint-devel/CA%2BiQ7x425PF0kzMS0XtUGMtNaniS%3DabK6m63-8nwOcgD5EzwJg%40mail.gmail.com
> > <https://groups.google.com/d/msgid/flint-devel/CA%2BiQ7x425PF0kzMS0XtUGMtNaniS%3DabK6m63-8nwOcgD5EzwJg%40mail.gmail.com?utm_medium=email&utm_source=footer>
> > .

> >
>
> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.

> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAA9NTHK4yapu8YnSVBg_XtsxAa9KFujs%3DAnRbHh_dskk6V4DAA%40mail.gmail.com.

[John Cremona]

I'm just wondering if anyone is reading this mailing list?

John

[Edgar Costa]

Plenty of us are, you got several replies.

If you didn't get the emails, perhaps because they were marked as spam, you can use the Google groups website.

--

---
You received this message because you are subscribed to the Google Groups "flint-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAD0p0K4WtNgLK43ZLRK7NxGiva8Mi%3DA3fGnLv3EtJ2VHhSn%3D%3Dw%40mail.gmail.com.

[John Cremona]

Apologies -- although I had emailed the list, I had not realised that I had my email preferences set to "No email". How stupid of me.

Now I am readong all the very helpful andc onstructive responses to my query -- thanks all!

In case anyone is interested, I am computing the integral homology of Bianchi groups. I don't need the transformation matrices, only the invariants, i.e. the diagonal entries.  I did have Sage code which managed the matrices (using pari) OK.

I am amused by the suggestion that my matrices are the "worst possible"!

[John Cremona]

In my test program (which I attached here) I just added the line 

fmpz_mat_hnf(M, M);

before

fmpz_mat_hnf(S, M);
and it ran OK without blowing up.  So I'll do the same in my production version.  I may also try Edgar's suggestion of doing this repeatedly with transposing (which means alternately doing  a right- and left- HNF). I have never thought about whether or not that would stabilise at the SNF.

[Claus Fieker]

On Fri, Feb 23, 2024 at 06:42:11AM -0800, John Cremona wrote:
> In my test program (which I attached here) I just added the line
> fmpz_mat_hnf(M, M);
> before
> fmpz_mat_hnf(S, M);
> and it ran OK without blowing up. So I'll do the same in my production
> version. I may also try Edgar's suggestion of doing this repeatedly with
> transposing (which means alternately doing a right- and left- HNF). I have
> never thought about whether or not that would stabilise at the SNF.

It will not. It will stabilize at a diagonal matrix - which then has to
be transformed into the SNF.
But that is cheap

> >>> <https://groups.google.com/d/msgid/flint-devel/CAD0p0K4WtNgLK43ZLRK7NxGiva8Mi%3DA3fGnLv3EtJ2VHhSn%3D%3Dw%40mail.gmail.com?utm_medium=email&utm_source=footer>
> >>> .

> >>>
> >>
>
> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.

> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/915689da-65e6-4db2-b7b3-8ec002d47940n%40googlegroups.com.

[John Cremona]

On Friday 23 February 2024 at 15:45:19 UTC fie...@mathematik.uni-kl.de wrote:

On Fri, Feb 23, 2024 at 06:42:11AM -0800, John Cremona wrote:
> In my test program (which I attached here) I just added the line
> fmpz_mat_hnf(M, M);
> before
> fmpz_mat_hnf(S, M);
> and it ran OK without blowing up. So I'll do the same in my production
> version. I may also try Edgar's suggestion of doing this repeatedly with
> transposing (which means alternately doing a right- and left- HNF). I have
> never thought about whether or not that would stabilise at the SNF.
It will not. It will stabilize at a diagonal matrix - which then has to
be transformed into the SNF.

Ah yes of course -- thanks, Claus.  I think I will start to do this as after >100 good runs I am seeing memory use expanding alarmingly in the current one (size 1861x1019).

[Fredrik Johansson]

Should we be doing something like this systematically in fmpz_mat_snf itself? Patches would be welcome!

Fredrik

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/Zdi9im-9tk5F57ir%40mathematik.uni-kl.de.

[John Cremona]

Well, having coded it up and checked on small examples, I set it going on my most recent failing case (size 1861 rows and 1001 columns) and the initial HNF step failed, in the sense that I killed it when the meory reached 50% (of 16G).  So the blow-up can also happend with HNF.  I could try with different specific HNF algorithms (e.g. classical) but not now, as it is the weekend (and I am retired, after all).

[Edgar Costa]

I'm happy to implement something during the workshop week.

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/643df4ef-4f44-4041-b8d0-f5f79314fa8fn%40googlegroups.com.