[with patch, positive review] speed up elementary_divisors for sparse integer matrices?

[jhpalmieri]

It seems to me that if mat is a sparse integer matrix, then

mat.dense_matrix().elementary_divisors()

is much faster than

mat.elementary_divisors()

Is this correct? I've checked this on certain families of matrices, but probably not extensively enough.

If so, we should change how elementary divisors for sparse integer matrices are computed. I've patched this, pretty naively, by sticking a new method in matrix_integer_sparse.pyx which just contains the above code. I would appreciate any comments or corrections.

[jhpalmieri]

This is a simple patch, but may raise complicated issues:

Here's the situation, as I see it: the code for sparse matrices needs a lot of work, and the particular problem on this ticket is one symptom of this. Right now, if you call elementary_divisors on a sparse matrix, it will likely take a long time, and might even crash. The patch proposes a fix to this -- converting to a dense matrix first. This will crash, too, for some matrices; I'm guessing that it will crash for strictly fewer matrices than elementary_divisors will. Therefore this patch is not a terrible idea.

It is rather makeshift, though. The right thing to do is to rewrite the sparse matrix code so that it works and is faster (and also less buggy -- see #5099, for instance).

This patch would be a temporary fix. Is it a good idea to include this now, with the idea to fix it the right way later? Or is it better to leave less than ideal code there now, first as an encouragement to fix it, and second so that we don't leave a patch like this in as a longterm solution?

[mabshoff]

Another simple fix that ought to be in 3.3.

Cheers,

Michael

[mabshoff]

Hi John,

this patch does not apply to my 3.3.rc0 merge tree. Please try to rebase it against 3.3.alpha6:

sage-3.3.rc0/devel/sage$ patch -p1 < trac_5142.patch 
patching file sage/matrix/matrix_integer_sparse.pyx
Hunk #1 FAILED at 296.
1 out of 1 hunk FAILED -- saving rejects to file sage/matrix/matrix_integer_sparse.pyx.rej

Cheers,

Michael

[jhpalmieri]

rebased against 3.3.alpha6

[jhpalmieri]

The new patch is just a rebase against 3.3.alpha6.

[AlexGhitza]

The patch looks good; I have three complaints regarding the docstring:

• in the description of the algorithm 'pari', you presumably mean "works robustly", since "works robustless" doesn't mean anything
• this is a one-line method, and it is pretty self-explanatory, so I don't think it needs a description of the implementation in the docstring
• the OUTPUT description claims that the method returns a list of int's; this is not true, since the output is a list of Integer's

I'll very happily give this a positive review once these issues are resolved. This patch is a good idea and, as John points out, similar things should be done for other methods for sparse matrices (determinant is another example). For the record, trying

sage: A = random_matrix(ZZ, 100, 100, sparse=True)
sage: time e = A.elementary_divisors()

simply fails with a mysterious TypeError in the current code, whereas with the patch applied it works in 1.44 seconds.

[jhpalmieri]

only apply this patch

[jhpalmieri]

Thanks for the comments. Here is a new patch which addresses them.

It also makes another change: my docstring for elementary_divisors was just copied from the docstring for the version in matrix_integer_dense, so since you pointed out the issues with my new docstring, I thought I should change that one, too. I made one change to the old docstring in addition to the ones you suggested: I deleted the paragraph which suggests that we use linbox. (We haven't used linbox for elementary_divisors for almost two years.)

[AlexGhitza]

Good stuff.

I'm upgrading this to a critical bug, due to the fact that it's not just a question of speed, but it also fails (see the example above).

[mabshoff]

Merged in Sage 3.3.rc0.

Cheers,

Michael

[mabshoff]

As John pointed out in  http://groups.google.com/group/sage-devel/t/312339a77bf58908: I merged 5142-rebased.patch instead of 5142-new.patch. So I reverted 5142-rebased.patch and merged 5142-new.patch in 3.3.rc1. Sorry for the screwup.

Cheers,

Michael