Reimplement matrix_integer_dense using FLINT

[mmasdeu]

In this ticket we reimplement all of matrix_integer_dense using FLINT fmpz_mat_t.

The speed-up is substantial. With the new code:

sage: A = Matrix(ZZ,1000,1000,range(10^6))
sage: %time B = A*A
CPU times: user 883 ms, sys: 3.33 ms, total: 887 ms
Wall time: 888 ms

This takes more than 1 minute in the same computer using Sage 6.3.

CC: @pjbruin

Component: linear algebra

Keywords: flint, matrix

Author: Marc Masdeu

Branch: 812a509

Reviewer: William Stein, Jeroen Demeyer

Issue created by migration from https://trac.sagemath.org/ticket/16803

[mmasdeu]

Branch: u/mmasdeu/16803-matrix_integer_dense_flint

[mmasdeu]

Commit: a51a2fb

[mmasdeu]

comment:2

By the way, after rebasing it does not compile. I will fix it tomorrow hopefully.

[fredrik-johansson]

comment:3

A word of caution about algorithm defaults: FLINT should generally be faster for small input but IML/Linbox etc. probably have much better code for things like nullspace, charpoly... as soon as the matrices start to get large.

Replacing _rational_kernel_iml with _rational_kernel_flint, for example, is strange. Surely both algorithms should be available.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

f1862b3	Trac 16803: Fixed compilation problems.
d180ee2	Trac 16803: Finished re-adding legacy functionality.

[sagetrac-git]

Changed commit from a51a2fb to d180ee2

[mmasdeu]

comment:5

I rewrote the code to have all the previous functionality available. The code compiles and passes all the doctests in the matrix/ folder. I will have it run all the doctests during the weekend.

[sagetrac-git]

Changed commit from d180ee2 to e26c347

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

e26c347

Trac 16803: Modified to pass doctests from outside matrix folder.

[sagetrac-git]

Changed commit from e26c347 to f2ac03c

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

a6e9a3d	Moved ntl dependency inside the pyx file.
f2ac03c	Trac 16803: Added sig_on/off, changed some defaults after benchmarking.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

888451b	Trac 16803: Fixed a missing sig_off() and doctests.
166bae5	All doctests pass, documentation builds well also.

[sagetrac-git]

Changed commit from f2ac03c to 166bae5

[williamstein]

comment:9

Hi,

A first remark playing around with timings is that in sage-6.3 somebody or something has definitely completely massively screwed up Linbox, or how Linbox is built, or how we use linbox, which was the default representation and library for arithmetic. Consider this:

sage: n=200; set_random_seed(0); a = random_matrix(ZZ,n); b = random_matrix(ZZ,n)
sage: %time c=a._multiply_linbox(b)
CPU times: user 326 ms, sys: 3.63 ms, total: 330 ms
Wall time: 347 ms
sage: %time c=a._multiply_multi_modular(b)
CPU times: user 51.5 ms, sys: 1.6 ms, total: 53.1 ms
Wall time: 56.1 ms
sage: n=400; set_random_seed(0); a = random_matrix(ZZ,n); b = random_matrix(ZZ,n)
sage: %time c=a._multiply_linbox(b)
CPU times: user 5.28 s, sys: 22.8 ms, total: 5.3 s
Wall time: 5.37 s
sage: %time c=a._multiply_multi_modular(b)
CPU times: user 230 ms, sys: 6.08 ms, total: 236 ms
Wall time: 248 ms

Linbox should easily beat the multiply_multi_modular approach -- that's just basically a naive version of the same thing linbox would be doing if it were properly built and linked. Instead, Linbox is massively slower. In fact, almost exactly as much slower as I would expect if the wrong BLAS were being linked. Note that I tested against the new flint backend and now a._multiply_multi_modular(b) is asymptotically similar to flint, though asymptotically (for n=2000 already), they are pretty close to the same in speed (8.55s versus 11.4s). In any case, somebody should really look into why linbox is so screwed up -- we switched to it initially since it was much better than doing a._multiply_multi_modular(b) and also our own strassen implementation.

[williamstein]

comment:10

OK, I've been skeptically looking at and testing this patch much of today. I made a few tiny changes, but frankly I think this is overall AWESOME. FLINT's linear algebra seems quite good.

The following three methods are now still gone. Add them back so I can run some randomized consistency tests more easily...

     a._multiply_classical      a._multiply_linbox         a._multiply_multi_modular

When these are added back, and I run tests of correctness by comparing with them, then I'll give this an enthusiastic positive review.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

7f4f903	Merge branch 'develop' into matflint
5f21e3d	Minor changes during review. Put back old multiplication functionality.

[sagetrac-git]

Changed commit from 166bae5 to 5f21e3d

[sagetrac-git]

Changed commit from 5f21e3d to a13b79c

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

a13b79c

Fixed errors, all doctests pass now.

[mmasdeu]

comment:15

I have added the three functions back, and run the long tests in sage/matrix/.

[williamstein]

comment:16

MAJOR PROBLEM! Now that multiply_multi_modular is back, I wrote code to do some random testing of correctness. I did not find situations where things are wrong... however, I quickly ran into major trouble because this ticket introduces a huge memory leak. I just defined two matrices and multiplied them a few hundred times, and had major ram problems. In particular, in sage-6.3 (without this branch):

sage: get_memory_usage()
1047.1875
sage: for i in range(50):
....:         a=random_matrix(ZZ,300); del a
....: 
sage: get_memory_usage()
1047.1875

With this ticket:

sage: get_memory_usage()
2132.09765625
sage: for i in range(50):
    a=random_matrix(ZZ,300); del a
....: 
sage: get_memory_usage()
2231.3203125

Conclusion: making nontrivial use of the code on this ticket will very quickly use up all RAM on your computer.

So... needs work until that is fixed.

[mmasdeu]

comment:77

Done.

[vbraun]

comment:79

Linear algebra should really give the right answer to a few ulp, and not just 1e-10. But the test is about testing that the division operator works, so ok.

[vbraun]

Changed branch from u/mmasdeu/16803-matflint-4 to 812a509

[jdemeyer]

Changed commit from 812a509 to none

[jdemeyer]

comment:81

Follow-ups: #17080 and #17090.

[jdemeyer]

comment:82

Did any of the reviewers (Volker?) actually read the code on this ticket? I've started doing this and there are many things to be improved... perhaps this ticket was merged too fast.

[vbraun]

comment:83

Is it perfect? No. But better than before? IMHO yes. Faster? For sure.

[williamstein]

comment:84

Replying to @jdemeyer:

Did any of the reviewers (Volker?) actually read the code on this ticket? I've started doing this and there are many things to be improved...

I certainly read the version of the code 5 weeks ago, and had many, many, many suggestions for improvement, found numerous major bugs, etc.

perhaps this ticket was merged too fast.

One has to balance that with the fact that for whatever reason, without this ticket a vast range of linear algebra in Sage is so utterly orders of magnitude too slow, as to make it painful and useless for research work. This new code is dramatically faster than before.

[jdemeyer]

comment:85

Replying to @williamstein:

This new code is dramatically faster than before.

I only hope we're not trading correctness for speed.

[williamstein]

comment:86

Replying to @jdemeyer:

Replying to @williamstein:

This new code is dramatically faster than before.

I only hope we're not trading correctness for speed.

Well when I looked over the code 5 weeks ago there are numerous memory leaks and other causes for concern (e.g., unitialized memory assumed to be 0). But those were I think all addressed. Did you find anything that is actually incorrect? Or does the overall code just make you unhappy and nervous? Or are you not happy building on FLINT?

[williamstein]

comment:87

I just have to add jdemeyer that I'm really glad you're looking at this too, since you have an amazing eye for quality (similar to your recent remarks about the pexpect interfaces).

[jdemeyer]

comment:88

Replying to @williamstein:

Did you find anything that is actually incorrect?

The segmentation fault in _lmul_, which is the reason for creating #17090.

Or does the overall code just make you unhappy and nervous?

When looking at this code, there are some sloppy things (like why is the _det_pari() function duplicated?). Overall, I have less faith in sloppy code than clear code.

Or are you not happy building on FLINT?

That's not the problem. In general, I think it's good to rely on external packages if they do a good job (which seems to be the case here).

[jdemeyer]

comment:90

What's with the changes in padics in this ticket? Surely, yet another mistake...

[jdemeyer]

comment:91

Yet another follow-up: #17094

[mmasdeu]

comment:92

I must say that IMO some of the comments (82 and 90) made above are not very constructive. jdemeyer questions the thoroughness of the reviewing process, while clearly he himself didn't read the comments at the beginning of the ticket (at least those that refer to the convenience of leaving in the legacy functions).

After the experience with this ticket I feel less inclined to take on the analogous job for matrix_rational_dense. Sorry!

[jdemeyer]

comment:93

Replying to @mmasdeu:

I must say that IMO some of the comments (82 and 90) made above are not very constructive.

[comment:82] was made out of frustration by seeing so many things here which are dubious. So okay, that's not constructive.

However in [comment:90], I noticed a mistake, mentioned it in the comment and fixed it in #17090. How is that not constructive?

while clearly he himself didn't read the comments at the beginning of the ticket (at least those that refer to the convenience of leaving in the legacy functions).

That's true, but after William Stein's comment, I have fixed that back in #17090.

I feel less inclined to take on the analogous job for matrix_rational_dense.

Please don't. I think it's absolutely good to improve matrices like you did here. The fact that I am making so many comments and creating follow-up tickets means that I care. I make mistakes, you make mistakes, that's normal and that's why have this review process.

My only complaint (and I still stand behind this) is that this ticket was given positive_review too soon. Ideally, those 3 follow-up tickets should have been reviewer patches on this ticket.

[mmasdeu]

comment:94

Replying to @jdemeyer:

[comment:82] was made out of frustration by seeing so many things here which are dubious. So okay, that's not constructive.

I can understand your frustration. But the bar is high enough (python + cython + sage oddities + legacy code + git + trac, and I am probably missing something here) that makes it very hard for people to contribute, and these comments do not help. But I take your last comment in good faith :-).

Also, if there are other things that are dubious please share them because I (and possibly others!) would like to learn from it.

However in [comment:90], I noticed a mistake, mentioned it in the comment and fixed it in #17090. How is that not constructive?

Surely I did not go and touch p-adics for no reason. These changes weren't there in the first version of the my patch, and I had to add them because during the time it took to from the first version until this ticket got closed other changes were being merged, and these made the patch to stop compiling. Changing those lines in p-adics is how I managed to fix it, but the hackish idea was not mine, I saw it used elsewhere. I appreciate you having fixed it back.

[jdemeyer]

comment:95

Replying to @mmasdeu:

Surely I did not go and touch p-adics for no reason. These changes weren't there in the first version of the my patch, and I had to add them because during the time it took to from the first version until this ticket got closed other changes were being merged, and these made the patch to stop compiling.

I think you simply did something compiling Sage, because those changes to p-adics are totally unrelated to matrices. I cannot imagine how the change to matrices could break building p-adics or vice versa. And indeed, I simply undid the changes in #17090 without any bad consequences.

[williamstein]

comment:96

Replying to @mmasdeu:

Replying to @jdemeyer:

[comment:82] was made out of frustration by seeing so many things here which are dubious. So okay, that's not constructive.

I can understand your frustration. But the bar is high enough (python + cython + sage oddities + legacy code + git + trac, and I am probably missing something here) that makes it very hard for people to contribute, and these comments do not help. But I take your last comment in good faith :-).

Marc -- don't you want to be a bad ass computer programmer, in addition to mathematician? Just redo matrix_rational_dense, and you'll be way closer :-)

[jdemeyer]

comment:97

Replying to @mmasdeu:

Also, if there are other things that are dubious please share them because I (and possibly others!) would like to learn from it.

Well, you can look at the follow-up tickets #17090 and #17094 and see what I changed (and preferably review those tickets).

Some general things which haven't been mentioned yet:

1. clean-up your code when you change things: when you remove code, try to remove everything related to the old code. When you change code, the old structure might not be applicable anymore, so refactor if needed (e.g. the sepatate branch for n <= 4 in determinant() no longer serves any purpose)
2. also fix documentation: when you change something, also change the corresponding documentation. When you copy/paste something, be sure to update documentation and doctests (see the various set_unsafe methods).
3. don't add unneeded import statements or includes, those add needless extra dependencies. Like, why did you add include "sage/libs/ntl/decl.pxi" in several places when you don't use any NTL code?
4. try to use correct formatting: spacing, indentation, docstring formatting (this was mostly good, but a bit messed up in src/module_list.py). If anything, this makes your code look more professional and less sloppy.
5. don't use assert for checking user input: use assert for when you know that something is true (but want to check just in case), not when you need to check that something is true.

[mmasdeu]

comment:98

Replying to @jdemeyer:

Thanks! These comments are helpful indeed.

William -- I rather try to be a better mathematician...it has the advantage that although you get constantly evaluated, the evaluations are mostly anonymous :-). But that was in my to-do list. And it still is, despite my previous comments in this ticket!