Opened 4 years ago
Closed 3 years ago
#22970 closed enhancement (fixed)
use flint fmpq_mat for rational dense matrices
Reported by:  vdelecroix  Owned by:  

Priority:  major  Milestone:  sage8.0 
Component:  linear algebra  Keywords:  thursdaysbdx 
Cc:  cpernet, mmasdeu  Merged in:  
Authors:  Vincent Delecroix  Reviewers:  Marc Masdeu 
Report Upstream:  N/A  Work issues:  
Branch:  cb0dae1 (Commits)  Commit:  cb0dae187cbbf8e8b594a64ef5c87553459b7f71 
Dependencies:  Stopgaps: 
Description (last modified by )
Attachments (2)
Change History (44)
comment:1 Changed 4 years ago by
 Branch set to u/vdelecroix/22970
 Commit set to 83e3a87687ee96e5a4b8cbcb91e825d8009bc4e2
 Component changed from PLEASE CHANGE to linear algebra
comment:2 Changed 4 years ago by
For a discussion about an annoying segfault related to this ticket, see this thread.
comment:3 Changed 4 years ago by
 Commit changed from 83e3a87687ee96e5a4b8cbcb91e825d8009bc4e2 to a254a43e3f4db5e850691bded1fbe7bb59ef0a03
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
3daa722  22970: more flint declarations

26822a2  22970: nonzero versions of gmp randomize functions

769d86c  22970: pari / fmpq_mat_t conversions

99cf141  22970: let rational dense matrix uses fmpq_mat_t

a254a43  22970: adapt codes using rational dense matrices

comment:4 Changed 4 years ago by
 Status changed from new to needs_review
comment:5 Changed 4 years ago by
 Description modified (diff)
comment:6 Changed 4 years ago by
 Keywords thursdaysbdx added
comment:7 Changed 4 years ago by
 Commit changed from a254a43e3f4db5e850691bded1fbe7bb59ef0a03 to 1621ac1e7945ee0df2aa92b47de5bf614bb9148d
Branch pushed to git repo; I updated commit sha1. New commits:
1621ac1  22970: use flint by default in det/inverse

comment:8 Changed 4 years ago by
New commits:
3daa722  22970: more flint declarations

26822a2  22970: nonzero versions of gmp randomize functions

769d86c  22970: pari / fmpq_mat_t conversions

99cf141  22970: let rational dense matrix uses fmpq_mat_t

a254a43  22970: adapt codes using rational dense matrices

1621ac1  22970: use flint by default in det/inverse

comment:9 Changed 4 years ago by
 Commit changed from 1621ac1e7945ee0df2aa92b47de5bf614bb9148d to 867ed5650aec68870282416055aaecaf0b735dfb
Branch pushed to git repo; I updated commit sha1. New commits:
867ed56  29970: echelonize/echelon_form using flint

comment:10 Changed 4 years ago by
 Commit changed from 867ed5650aec68870282416055aaecaf0b735dfb to d1c8d52b65b83268fc1ff6eca155cf85c7797850
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
d1c8d52  29970: echelonize/echelon_form using flint

comment:11 Changed 4 years ago by
Thanks for the work Vincent! I think that it's my duty to review this one :).
While I update my develop branch and check the ticket out, could you check that that the docstring in the determinant method is correct? In algorithm
, I see
None
missing.
Also, do you happen to have timings? that compare with the previous implementation? Otherwise I can produce some...
comment:12 Changed 4 years ago by
 Status changed from needs_review to needs_work
Hi Marc,
Thanks for proposing the review! I need some more work on my branch as I have strange failing doctests in sage/modular
that I don't understand right now.
I will provide proper timings for the following methods (tell me if you have more in mind):
determinant
multiplication
echelonize
charpoly
(not yet in flint, but can go through integer matrices)minpoly
(not yet in flint, but can go through integer matrices)
Changed 4 years ago by
comment:13 Changed 4 years ago by
 Commit changed from d1c8d52b65b83268fc1ff6eca155cf85c7797850 to edfa6f46b5bab2c2f9e555c2cd673a4795dd4983
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
f01157c  22970: more flint declarations

acce30c  22970: nonzero versions of gmp randomize functions

be5ccb4  22970: pari / fmpq_mat_t conversions

9ba2c4d  22970: let rational dense matrix uses fmpq_mat_t

150631d  22970: adapt codes using rational dense matrices

a3b6b37  22970: use flint by default in det/inverse

10f6d8c  29970: echelonize/echelon_form using flint

aef538e  22970: derandomize tests in free_module_integer.py

edfa6f4  22970: various cleaning in matrix_rational_dense

comment:14 followup: ↓ 15 Changed 4 years ago by
At commit edfa6f46b5
using test_file.sage I obtain (within different session because of caching)
sage: %runfile test_file.sage sage: test(seed=0) f = q + 1/2*a1*q^3 + (4*a1  26)*q^5 + O(q^6) K = Number Field in a1 with defining polynomial x^2 + 32*x  9984 a = 1/2*a1 True False
and
sage: %runfile test_file.sage sage: test(seed=1) f = q + 1/2*a1*q^3 + (4*a1  26)*q^5 + O(q^6) K = Number Field in a1 with defining polynomial x^2 + 32*x  9984 a = 1/2*a1 False True
And it looks like a bug to me that this result depends on the seed of the random generator...
comment:15 in reply to: ↑ 14 Changed 4 years ago by
Replying to comment:14: Discussed at this sagedevel thread.
Changed 4 years ago by
comment:16 Changed 4 years ago by
 echelonize: flint is better except when coefficients are huge in which case multimodular is best
nrows=5 ncols=5 rank=5, density=1.0000 nbits(height)=100 flint 0.0565 multimodular 0.3284 padic 0.5626 nrows=10 ncols=10 rank=5, density=1.0000 nbits(height)=37 flint 0.1063 multimodular 2.6420 padic 7.6511 nrows=100 ncols=100 rank=20, density=0.1799 nbits(height)=12 flint 9.1834 multimodular 12.2192 padic 74.8152 nrows=10 ncols=10 rank=10, density=1.0000 nbits(height)=200 flint 10.1606 multimodular 1.5586 padic 2.1540
 determinant
dim=3 density=0.7778 nbits(height)=3 flint 0.0092 integer 0.1334 pari 0.0017 dim=3 density=1.0000 nbits(height)=100 flint 0.0113 integer 0.1188 pari 0.0269 dim=30 density=1.0000 nbits(height)=5 flint 0.9006 integer 1.0754 pari 34.5452 dim=30 density=0.2544 nbits(height)=5 flint 0.4389 integer 0.9986 pari 179.9622
 charpoly
dim=3 density=0.7778 nbits(height)=2 flint 0.3875 linbox 1.3019 generic 0.0431 dim=4 density=0.5625 nbits(height)=2 flint 0.1776 linbox 1.2013 generic 0.0450 dim=5 density=0.6400 nbits(height)=2 flint 0.1737 linbox 1.2952 generic 0.0532 dim=6 density=0.6111 nbits(height)=2 flint 0.1850 linbox 1.2773 generic 0.0651 dim=10 density=0.4100 nbits(height)=2 flint 0.2343 linbox 1.4250 generic 0.1818 dim=10 density=1.0000 nbits(height)=2 flint 0.2521 linbox 1.4264 generic 0.4842 dim=30 density=1.0000 nbits(height)=2 flint 4.8531 linbox 3.4320 generic 50.0362
comment:17 Changed 4 years ago by
 Commit changed from edfa6f46b5bab2c2f9e555c2cd673a4795dd4983 to 37614675ea8c814e0af3e2fb6664f6c5c0ade0f1
comment:18 Changed 4 years ago by
 Status changed from needs_work to needs_review
Needs review again. I tried to make independent of the random seed the annoying doctests in sage/modular...
comment:19 Changed 4 years ago by
patchbot is green!! (with a weird ascii complaints)
comment:20 followup: ↓ 21 Changed 4 years ago by
On my computer, a call to
sage: CuspForms(101,12).hecke_matrix(3)
takes 765 seconds with the current develop branch, and 1157 seconds with the proposed patch. I think that most of it is explained by the fact that the method echelon_form is called 1394 times on quite large matrices (dense, over QQ). I was actually looking for speedups in existing highlevel functionality (like this one), so this is a bit disappointing.
Maybe the default algorithm should be changed if the entries are large?
comment:21 in reply to: ↑ 20 Changed 4 years ago by
Replying to mmasdeu:
On my computer, a call to
sage: CuspForms(101,12).hecke_matrix(3)takes 765 seconds with the current develop branch, and 1157 seconds with the proposed patch. I think that most of it is explained by the fact that the method echelon_form is called 1394 times on quite large matrices (dense, over QQ). I was actually looking for speedups in existing highlevel functionality (like this one), so this is a bit disappointing.
Maybe the default algorithm should be changed if the entries are large?
I did not do proper tunings (since I will redo it after #22924 anyway). However, I can try to do something better for echelon_form/echelonize
. Do you have other relevant examples in mind?
comment:22 Changed 4 years ago by
Not at the moment. I haven't done exhaustive tests either, but it seems that in many cases the performance is similar.
As for other examples, many operations with subquotients involve doing this kind of linear algebra.
comment:23 Changed 4 years ago by
Some remarks
 flint has several algorithms for echelon form and this is not finely tuned.
 On very specific instances, multimodular can be much faster (up to x100). This is the reason why your example
CuspForms(101,12).hecke_matrix(3)
gets slower. I have no idea how to detect such instances.  for
random_matrix(QQ, ...)
(ie independent entries) flint looks like to be the best option excepted in the range10 <= nrows=ncols <= 20
where flint tuning (over Z) is very bad
One way to go might be to have a global option that controls the default behavior of the echelon form. I am thinking of something that can be used like
def my_modular_form_algorithm(): with rational_echelon_form('multimodular'): # some relevant code here ....
What do you think?
comment:24 Changed 4 years ago by
Second option, allow the option algorithm=my_function
where my_function
takes a matrix as input and returns a valid algorithm for echelon form. (I do prefer this one better)
EDIT: this will not work since echelon_form
is indirectly called in modular function algorithms...
comment:25 Changed 4 years ago by
Implementing the idea from comment:23 I got
sage: %time m = CuspForms(31, 12).hecke_matrix(3) CPU times: user 7.52 s, sys: 63.3 ms, total: 7.59 s Wall time: 7.55 s sage: %time c = CuspForms(51, 12).hecke_matrix(3) CPU times: user 2min 13s, sys: 203 ms, total: 2min 13s Wall time: 2min 13s sage: %time c = CuspForms(101, 12).hecke_matrix(3) CPU times: user 5min 10s, sys: 663 ms, total: 5min 10s Wall time: 5min 10s
against the following in develop
sage: %time m = CuspForms(31, 12).hecke_matrix(3) CPU times: user 13.5 s, sys: 43.3 ms, total: 13.6 s Wall time: 13.5 s sage: %time c = CuspForms(51, 12).hecke_matrix(3) CPU times: user 3min 8s, sys: 250 ms, total: 3min 9s Wall time: 3min 9s %time c = CuspForms(101, 12).hecke_matrix(3) CPU times: user 8min 7s, sys: 527 ms, total: 8min 8s Wall time: 8min 8s
(so close to x2 improvement)
I propose to postpone this improvement to another ticket since this one is already a big change. If you agree I will open a new one and work on it as soon as this one is closed.
comment:26 Changed 4 years ago by
 Status changed from needs_review to needs_work
1) I agree that this can be dealt with in another ticket, but it is moderately urgent. ModularSymbols? / ModularForms? becoming twice slower than currently is not good news...
Another example of a reasonable calculation: with the Sage 7.5.1 the following calculation takes 22 seconds on my laptop. It takes 1min22s (!) with the proposed patch.
sage: %time A = ModularSymbols(11,60).hecke_matrix(3)
2) While reviewing, in sage/modular/hecke/module.py
there is a doctest that changed in a suspicious way. It seems that a number field was created with a different generator and that causes the output to look different. I propose that the line that causes trouble gets changed into
sage: [o.minpoly() for o in A.system_of_eigenvalues(3)] [x  1, x + 1, x^2  2*x  2]
The output that appears should not depend on the implementation.
3) There seems to be quite a mess on what are the acceptable values to the 'algorithm' parameter. Maybe it's a good time to change that! My suggestion would be that the default would always be None, and then the docstring would specify what is done in this case.
comment:27 followup: ↓ 28 Changed 4 years ago by
I might have an explanation for why Flint is slow for these modular symbol matrices. See:
https://github.com/wbhart/flint2/issues/335
It can be fixed in Flint (if my analysis is correct), but this requires some amount of coding.
comment:28 in reply to: ↑ 27 Changed 4 years ago by
Replying to fredrik.johansson:
I might have an explanation for why Flint is slow for these modular symbol matrices. See:
https://github.com/wbhart/flint2/issues/335
It can be fixed in Flint (if my analysis is correct), but this requires some amount of coding.
Would be tremendous!
comment:29 Changed 4 years ago by
@mmasdeu: for the sake of this ticket I will set the default to multimodular
that should actually speed up all modular symbol computations. Then #23026 will deal with making flint the default (that is much faster on random instances).
comment:30 Changed 4 years ago by
That sounds reasonable...at least there is no regression in speed (hopefully!).
comment:31 Changed 4 years ago by
 Commit changed from 37614675ea8c814e0af3e2fb6664f6c5c0ade0f1 to b5f5600ddc92113f3307b07a13c3d60e0bb5a736
comment:32 Changed 4 years ago by
 Commit changed from b5f5600ddc92113f3307b07a13c3d60e0bb5a736 to 5f2e1bac18839e95198e455897ab81cc35c8923f
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
14d31c5  22970: better doctest in modular/hecke/module.py

0c8b2ce  22970: improved doc + simpler echelonize/echelon_form logic

c172064  22970: _new_matrix/from_columns/add_to_entry

5f2e1ba  22970: small optimization in modsym/ambient.py

comment:33 Changed 4 years ago by
at commit b5f5600 I got
sage: %time m = CuspForms(31, 12).hecke_matrix(3) CPU times: user 11.7 s, sys: 40 ms, total: 11.7 s Wall time: 11.7 s sage: %time c = CuspForms(51, 12).hecke_matrix(3) CPU times: user 2min 52s, sys: 283 ms, total: 2min 52s Wall time: 2min 52s sage: %time c = CuspForms(101, 12).hecke_matrix(3) CPU times: user 7min 17s, sys: 620 ms, total: 7min 18s Wall time: 7min 18s
and at 5f2e1ba the even better
sage: %time m = CuspForms(31, 12).hecke_matrix(3) CPU times: user 7.5 s, sys: 36.7 ms, total: 7.53 s Wall time: 7.51 s sage: %time c = CuspForms(51, 12).hecke_matrix(3) CPU times: user 2min 10s, sys: 223 ms, total: 2min 10s Wall time: 2min 10s sage: %time c = CuspForms(101, 12).hecke_matrix(3) CPU times: user 5min 9s, sys: 603 ms, total: 5min 10s Wall time: 5min 10s
comment:34 Changed 4 years ago by
 Status changed from needs_work to needs_review
comment:35 Changed 4 years ago by
And tests also pass on beta7!
comment:36 Changed 4 years ago by
 Reviewers set to Marc Masdeu
 Status changed from needs_review to positive_review
This is definitely more than satisfying. The code looks very clean, and all tests pass. Thanks for the great work!
comment:38 Changed 4 years ago by
Too bad. The code is fine on beta7. Can I try to fix it by merging another ticket?
comment:39 Changed 4 years ago by
 Commit changed from 5f2e1bac18839e95198e455897ab81cc35c8923f to cb0dae187cbbf8e8b594a64ef5c87553459b7f71
Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:
86a305b  22970: adapt codes using rational dense matrices

b5320ef  22970: use flint by default in det/inverse

fad1b0d  29970: echelonize/echelon_form using flint

b68993f  22970: derandomize tests in free_module_integer.py

91c5835  22970: various cleaning in matrix_rational_dense

9f83912  22970: fix doctests in sage/modular/

bcf6b7d  22970: better doctest in modular/hecke/module.py

14f3fdb  22970: improved doc + simpler echelonize/echelon_form logic

fceef50  22970: _new_matrix/from_columns/add_to_entry

cb0dae1  22970: small optimization in modsym/ambient.py

comment:40 Changed 4 years ago by
 Status changed from needs_work to needs_review
rebased (I launched the patchbot)
comment:41 Changed 4 years ago by
 Status changed from needs_review to positive_review
patchbot still happy!
comment:42 Changed 3 years ago by
 Branch changed from u/vdelecroix/22970 to cb0dae187cbbf8e8b594a64ef5c87553459b7f71
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
22970: more flint declarations
22970: pari / fmpq_mat_t conversions
22970: implement libs/flint/randomize.pxd
22970: let rational dense matrix uses fmpq_mat_t
22970: adapt codes using rational dense matrices