Opened 3 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: sage-8.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 vdelecroix)

The flint library has a builtin type for rational matrices: fmpq_mat_t. This ticket proposes to replace the current array of mpq_t wrapped by Matrix_rational_dense in favor of fmpq_mat_t.

(See also the related #22924 and #22872)

Attachments (2)

test_file.sage (730 bytes) - added by vdelecroix 3 years ago.
benchmark_matrices.py (6.6 KB) - added by vdelecroix 3 years ago.

Download all attachments as: .zip

Change History (44)

comment:1 Changed 3 years ago by vdelecroix

  • Branch set to u/vdelecroix/22970
  • Commit set to 83e3a87687ee96e5a4b8cbcb91e825d8009bc4e2
  • Component changed from PLEASE CHANGE to linear algebra

New commits:

f7d368f22970: more flint declarations
fce2da522970: pari / fmpq_mat_t conversions
b52e2c322970: implement libs/flint/randomize.pxd
d5ea72b22970: let rational dense matrix uses fmpq_mat_t
83e3a8722970: adapt codes using rational dense matrices

comment:2 Changed 3 years ago by vdelecroix

For a discussion about an annoying segfault related to this ticket, see this thread.

comment:3 Changed 3 years ago by git

  • Commit changed from 83e3a87687ee96e5a4b8cbcb91e825d8009bc4e2 to a254a43e3f4db5e850691bded1fbe7bb59ef0a03

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

3daa72222970: more flint declarations
26822a222970: nonzero versions of gmp randomize functions
769d86c22970: pari / fmpq_mat_t conversions
99cf14122970: let rational dense matrix uses fmpq_mat_t
a254a4322970: adapt codes using rational dense matrices

comment:4 Changed 3 years ago by vdelecroix

  • Status changed from new to needs_review

comment:5 Changed 3 years ago by vdelecroix

  • Description modified (diff)

comment:6 Changed 3 years ago by vdelecroix

  • Keywords thursdaysbdx added

comment:7 Changed 3 years ago by git

  • Commit changed from a254a43e3f4db5e850691bded1fbe7bb59ef0a03 to 1621ac1e7945ee0df2aa92b47de5bf614bb9148d

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

1621ac122970: use flint by default in det/inverse

comment:8 Changed 3 years ago by vdelecroix

New commits:

3daa72222970: more flint declarations
26822a222970: nonzero versions of gmp randomize functions
769d86c22970: pari / fmpq_mat_t conversions
99cf14122970: let rational dense matrix uses fmpq_mat_t
a254a4322970: adapt codes using rational dense matrices
1621ac122970: use flint by default in det/inverse

comment:9 Changed 3 years ago by git

  • Commit changed from 1621ac1e7945ee0df2aa92b47de5bf614bb9148d to 867ed5650aec68870282416055aaecaf0b735dfb

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

867ed5629970: echelonize/echelon_form using flint

comment:10 Changed 3 years ago by git

  • Commit changed from 867ed5650aec68870282416055aaecaf0b735dfb to d1c8d52b65b83268fc1ff6eca155cf85c7797850

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

d1c8d5229970: echelonize/echelon_form using flint

comment:11 Changed 3 years ago by mmasdeu

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 3 years ago by vdelecroix

  • 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)
Last edited 3 years ago by vdelecroix (previous) (diff)

Changed 3 years ago by vdelecroix

comment:13 Changed 3 years ago by git

  • Commit changed from d1c8d52b65b83268fc1ff6eca155cf85c7797850 to edfa6f46b5bab2c2f9e555c2cd673a4795dd4983

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

f01157c22970: more flint declarations
acce30c22970: nonzero versions of gmp randomize functions
be5ccb422970: pari / fmpq_mat_t conversions
9ba2c4d22970: let rational dense matrix uses fmpq_mat_t
150631d22970: adapt codes using rational dense matrices
a3b6b3722970: use flint by default in det/inverse
10f6d8c29970: echelonize/echelon_form using flint
aef538e22970: derandomize tests in free_module_integer.py
edfa6f422970: various cleaning in matrix_rational_dense

comment:14 follow-up: Changed 3 years ago by vdelecroix

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...

Last edited 3 years ago by vdelecroix (previous) (diff)

comment:15 in reply to: ↑ 14 Changed 3 years ago by vdelecroix

Replying to comment:14: Discussed at this sage-devel thread.

Changed 3 years ago by vdelecroix

comment:16 Changed 3 years ago by vdelecroix

From benchmark_matrices.py

  • 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 3 years ago by git

  • Commit changed from edfa6f46b5bab2c2f9e555c2cd673a4795dd4983 to 37614675ea8c814e0af3e2fb6664f6c5c0ade0f1

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

ba1787a22970: various cleaning in matrix_rational_dense
376146722970: fix doctests in sage/modular/

comment:18 Changed 3 years ago by vdelecroix

  • 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...

Last edited 3 years ago by vdelecroix (previous) (diff)

comment:19 Changed 3 years ago by vdelecroix

patchbot is green!! (with a weird ascii complaints)

comment:20 follow-up: Changed 3 years ago by 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 speed-ups in existing high-level 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 3 years ago by vdelecroix

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 speed-ups in existing high-level 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 3 years ago by mmasdeu

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 3 years ago by vdelecroix

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 range 10 <= 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 3 years ago by vdelecroix

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...

Last edited 3 years ago by vdelecroix (previous) (diff)

comment:25 Changed 3 years ago by vdelecroix

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 3 years ago by mmasdeu

  • 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 follow-up: Changed 3 years ago by 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.

comment:28 in reply to: ↑ 27 Changed 3 years ago by vdelecroix

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 3 years ago by vdelecroix

@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 3 years ago by mmasdeu

That sounds reasonable...at least there is no regression in speed (hopefully!).

comment:31 Changed 3 years ago by git

  • Commit changed from 37614675ea8c814e0af3e2fb6664f6c5c0ade0f1 to b5f5600ddc92113f3307b07a13c3d60e0bb5a736

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

948cab122970: better doctest in modular/hecke/module.py
b5f560022970: improved doc + simpler echelonize/echelon_form logic

comment:32 Changed 3 years ago by git

  • Commit changed from b5f5600ddc92113f3307b07a13c3d60e0bb5a736 to 5f2e1bac18839e95198e455897ab81cc35c8923f

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

14d31c522970: better doctest in modular/hecke/module.py
0c8b2ce22970: improved doc + simpler echelonize/echelon_form logic
c17206422970: _new_matrix/from_columns/add_to_entry
5f2e1ba22970: small optimization in modsym/ambient.py

comment:33 Changed 3 years ago by vdelecroix

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
Last edited 3 years ago by vdelecroix (previous) (diff)

comment:34 Changed 3 years ago by vdelecroix

  • Status changed from needs_work to needs_review

comment:35 Changed 3 years ago by vdelecroix

And tests also pass on beta7!

comment:36 Changed 3 years ago by mmasdeu

  • 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:37 Changed 3 years ago by vbraun

  • Status changed from positive_review to needs_work

Merge conflict

comment:38 Changed 3 years ago by vdelecroix

Too bad. The code is fine on beta7. Can I try to fix it by merging another ticket?

comment:39 Changed 3 years ago by git

  • Commit changed from 5f2e1bac18839e95198e455897ab81cc35c8923f to cb0dae187cbbf8e8b594a64ef5c87553459b7f71

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

86a305b22970: adapt codes using rational dense matrices
b5320ef22970: use flint by default in det/inverse
fad1b0d29970: echelonize/echelon_form using flint
b68993f22970: derandomize tests in free_module_integer.py
91c583522970: various cleaning in matrix_rational_dense
9f8391222970: fix doctests in sage/modular/
bcf6b7d22970: better doctest in modular/hecke/module.py
14f3fdb22970: improved doc + simpler echelonize/echelon_form logic
fceef5022970: _new_matrix/from_columns/add_to_entry
cb0dae122970: small optimization in modsym/ambient.py

comment:40 Changed 3 years ago by vdelecroix

  • Status changed from needs_work to needs_review

rebased (I launched the patchbot)

comment:41 Changed 3 years ago by vdelecroix

  • Status changed from needs_review to positive_review

patchbot still happy!

comment:42 Changed 3 years ago by vbraun

  • Branch changed from u/vdelecroix/22970 to cb0dae187cbbf8e8b594a64ef5c87553459b7f71
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.