Opened 12 years ago

Closed 12 years ago

Last modified 12 years ago

#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 mabshoff)

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)

903-charpoly_pari.patch (2.5 KB) - added by AlexGhitza 12 years ago.
903-charpoly_pari_doctests.patch (2.8 KB) - added by AlexGhitza 12 years ago.
apply after 903-charpoly_pari.patch
903-ncalexan-charpoly-over-number-field-1.patch (4.1 KB) - added by ncalexan 12 years ago.
903-small_fix.patch (1.3 KB) - added by AlexGhitza 12 years ago.
apply after Nick's patch
903-ncalexan-charpoly-over-number-field-2.patch (3.7 KB) - added by ncalexan 12 years ago.

Download all attachments as: .zip

Change History (23)

comment:1 Changed 12 years ago by mabshoff

  • Description modified (diff)

comment:2 Changed 12 years ago by ncalexan

The following IRC transcript is relevant.

ncalexan-827: wstein-1606: is there a good way to find degree 1 primes in a number field.
ncalexan-827: ?
wstein-1606: theoretically or in sage right now?
ncalexan-827: wstein-1606: right now, but also in theory.
wstein-1606: for all but finitely many primes you factor the defining poly of the field and see if it splits completely.
dmharvey left the chat room.
wstein-1606: you can do that more quickly by taking the gcd mod p with x^p - x
ncalexan-827: Right.
wstein-1606: or something like that.
wstein-1606: I bet you get the idea.
ncalexan-827: Doesn't this only work if your ring of integers is generated by a single element?
wstein-1606: nick -- that's why I said "all but finitely many".
wstein-1606: For all but a couple primes the ring of integers is monogenic.
wstein-1606: just avoid primes dividing the discriminant.
wstein-1606: for those do the usual factorization algorithm.
ncalexan-827: wstein-1606: sure.  It sounds like a degree_one_primes() generator would be handy for multi modular algorithms, then.
wstein-1606: yep!
wstein-1606: But beware -- it is very very very hard to find any degree one primes if the field isn't cyclotomic.
wstein-1606: For a galois S_n extension of degree n, only 1/n! of the primes split completely (by cebo.)
ncalexan-827: Hmm.  So it's not even worth it?
wstein-1606: Well it's very worth it in the cyclotomic case.
wstein-1606: more generally, in the abelian case.
wstein-1606: E.g., for quadratic fields.
wstein-1606: and even 1/n! isn't bad if n is small, which it often is.

comment:3 follow-up: Changed 12 years ago by cremona

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,xp-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 cremona

Replying to cremona:

Sorry about bad formatting.

comment:5 Changed 12 years ago by ncalexan

See #2835 for a possible primes of degree one iterator implementation.

Changed 12 years ago by AlexGhitza

comment:6 Changed 12 years ago by AlexGhitza

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

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

Changed 12 years ago by AlexGhitza

apply after 903-charpoly_pari.patch

comment:8 Changed 12 years ago by AlexGhitza

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

The patches apply and pass tests for me. Nick, what are your thoughts?

comment:10 Changed 12 years ago by ncalexan

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

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

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

Changed 12 years ago by AlexGhitza

apply after Nick's patch

comment:13 Changed 12 years ago by AlexGhitza

  • 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

comment:14 Changed 12 years ago by ncalexan

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

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

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 mabshoff

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

Yes, that patch was the only one that needed merging.

Note: See TracTickets for help on using tickets.