#903 closed defect (fixed)
[with patch, very positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
Reported by: | was | Owned by: | was |
---|---|---|---|
Priority: | major | Milestone: | sage-3.0 |
Component: | linear algebra | Keywords: | |
Cc: | Merged in: | ||
Authors: | Reviewers: | ||
Report Upstream: | Work issues: | ||
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
See the below examples that illustrate that charpoly(A), where A is over a number field, currently totally sucks. This has very bad implications for modular forms computations. Even switching to use PARI for charpoly would give a significant speedup (which is a lot faster than Magma, interestingly).
K.<a> = NumberField(x^2 + 17) a^2 /// -17
n = 40 m = matrix(K, n, [(a+1)^randint(0,3) for _ in xrange(n^2)]) m /// 40 x 40 dense matrix over Number Field in a with defining polynomial x^2 + 17
m[0] # first row /// (2*a - 16, 1, -14*a - 50, a + 1, 2*a - 16, 2*a - 16, a + 1, -14*a - 50, 1, -14*a - 50, 2*a - 16, a + 1, 1, a + 1, a + 1, a + 1, a + 1, -14*a - 50, -14*a - 50, 2*a - 16, 1, a + 1, 1, a + 1, 2*a - 16, -14*a - 50, 1, 2*a - 16, 2*a - 16, 1, a + 1, a + 1, 1, 2*a - 16, -14*a - 50, 2*a - 16, 2*a - 16, -14*a - 50, 2*a - 16, -14*a - 50)
time k = m*m /// CPU time: 0.14 s, Wall time: 0.27 s
time f=m.charpoly() /// CPU time: 23.93 s, Wall time: 26.22 s
NOTE: Sage should use PARI for charpoly's over number fields, but currently it doesn't. Notice how much faster PARI is at the same computation!
m._clear_cache() g = pari(m) time h = g.charpoly() /// Time: CPU 2.52 s, Wall: 2.76 s
But a multimodular algorithm should blow away all of them. Student project? Basically, charpoly mod p is extremely fast in Sage, and one could reduce modulo primes of a number field, do the computation mod p, then lift, and get the correct result. I know of no implementations of this. A student project could be to implement this and tune it to be the fastest program in the world for charpoly over number fields.
time m.echelonize() /// CPU time: 0.35 s, Wall time: 0.35 s
%magma R<x> := PolynomialRing(RationalField()); K<a> := NumberField(x^2 + 17); n := 40; m := MatrixAlgebra(K, n)![(1+a)^Random(0, 3) : i in [1..n^2]]; time k := m*m; time f := CharacteristicPolynomial(m); time e := EchelonForm(m); /// Time: 0.000 Time: 15.220 Time: 0.150
Attachments (5)
Change History (23)
comment:1 Changed 12 years ago by
- Description modified (diff)
comment:2 Changed 12 years ago by
comment:3 follow-up: ↓ 4 Changed 12 years ago by
I needed exactly this for my work on Unimodular Integer Circulants (To appear in: Mathematics of Computation) or see my web page. At the time I used pari and have been meaning to try it in Sage too. Rather than take gcd(f,x^{p-x) -- all mof p of course -- , I took Mod(x,f*Mod(1,p))}p and that worked pretty well. My worst example: deg(f) = 18, the first degree 1 prime was about 250 million and did not work for me so I sed the second one which was around 2 billion. This poly has a Galois group of order only 2*(n/2)! ; if it been the generic n! I would probably still be waiting.
This suggests to me that for anything other than small degrees the implementation will need to use primes of higher degree.
comment:4 in reply to: ↑ 3 Changed 12 years ago by
Replying to cremona:
Sorry about bad formatting.
comment:5 Changed 12 years ago by
See #2835 for a possible primes of degree one iterator implementation.
Changed 12 years ago by
comment:6 Changed 12 years ago by
- Summary changed from charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
I suggest that we implement the easy fix and get that merged, and open an enhancement ticket along the lines "implement multimodular algorithm for charpoly over number fields".
I'm attaching a patch with the easy fix. On my machine, William's example runs 15 times faster after the patch than before.
comment:7 Changed 12 years ago by
- Summary changed from [with patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, with negative review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
I am not happy with this patch. I want to see an explicit doctest (not random) that has a polynomial defined over the number field, not over QQ.
I also want to see doctests for stacked relative fields -- I.e.
x = QQ['x'].gen() K.<a> = NumberField(x^2 - 2) L.<b> = K.extension(x^3 - a) M.<c> = L.extension(x^2 - a*x + b)
and see it work for matrices over all those fields.
This is hard -- we need polynomials over arbitrary pari number fields, and it doesn't work for #2329 yet, so it shouldn't work here yet.
comment:8 Changed 12 years ago by
- Summary changed from [with patch, with negative review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with additional patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
The point about the doctests is correct. I put up a new patch that adds the appropriate doctests and fixes a small bug (due to some unrelated weirdness that I don't want to deal with right now).
The doctests pass and I checked the results. It looks to me like it's working, so I'm not sure I agree with the last part of the last comment.
comment:9 Changed 12 years ago by
The patches apply and pass tests for me. Nick, what are your thoughts?
Changed 12 years ago by
comment:10 Changed 12 years ago by
- Summary changed from [with additional patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
I prefer the attached patch because it does not require string parsing, at least not in the code. It will be more robust to changes to PARI and uses the Sage library. It's shorter, too, and I prefer the factored-out function.
comment:11 Changed 12 years ago by
- Summary changed from [with patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
Nick's patch applies and passes the tests. I think it is preferable since it avoids string parsing. I'm going to give it a positive review. Alex, if you disagree, feel free to comment and change it.
comment:12 Changed 12 years ago by
- Summary changed from [with patch, positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, positive review pending fix of doctest failures] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
Yes, I much prefer Nick's code over mine. I was struggling to get things from Pari into Sage and the way I was doing things was obscure and rather fragile.
There is one problem I run into with this: 2 doctests in number_field_element.pyx fail. Probably an easy fix?
comment:13 Changed 12 years ago by
- Summary changed from [with patch, positive review pending fix of doctest failures] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
It was indeed an easy fix. Pari treats the t_POLMOD 0 in a different way than the others, so the code chocked on that. I've put up a patch that fixes this and adds a direct doctest for this behavior.
We should apply 903-ncalexan-charpoly-over-number-field-1.patch and 903-small_fix.patch. When this ticket gets closed we should make a new enhancement ticket containing the above discussion of a multimodular algorithm.
BTW, in 3.0.alpha5, before the patches:
sage: K.<a> = NumberField(x^2 + 17) sage: n = 40 sage: m = matrix(K, n, [(a+1)^randint(0,3) for _ in xrange(n^2)]) sage: time f = m.charpoly('Z') CPU times: user 59.69 s, sys: 0.01 s, total: 59.70 s Wall time: 59.70
and after the patches:
sage: K.<a> = NumberField(x^2 + 17) sage: n = 40 sage: m = matrix(K, n, [(a+1)^randint(0,3) for _ in xrange(n^2)]) sage: time f = m.charpoly('Z') CPU times: user 3.96 s, sys: 0.02 s, total: 3.98 s Wall time: 3.98
Changed 12 years ago by
comment:14 Changed 12 years ago by
- Summary changed from [with patch, positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
903-ncalexan-charpoly-over-number-field-2.patch
stands alone, and doesn't muck about with PARI at all: it turns out that polynomials over number fields work to the extent that special parsing is not necessary.
BTW, great team effort guys!
comment:15 Changed 12 years ago by
- Summary changed from [with patch, needs review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist) to [with patch, very positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)
I *love* the end result, it's so clean and simple (I wish there was a way to delete my old patches from this ticket, I'm somewhat embarrassed :)
Anyway, it works great and fast and it should be merged.
comment:16 Changed 12 years ago by
Nick, Alex,
am I correct to understand that I must only apply 903-ncalexan-charpoly-over-number-field-2.patch?
it applies by itself and does build, and doctests fine.
Cheers,
Michael
comment:17 Changed 12 years ago by
- Resolution set to fixed
- Status changed from new to closed
Merged 903-ncalexan-charpoly-over-number-field-2.patch in Sage 3.0.alpha6 only. In case that is not what is intended please comment here and I will merge those patches.
Cheers,
Michael
comment:18 Changed 12 years ago by
Yes, that patch was the only one that needed merging.
The following IRC transcript is relevant.