Opened 6 years ago
Last modified 2 years ago
#16812 needs_work enhancement
use FLINT to speed up Chebyshev T polynomial creation
Reported by:  rws  Owned by:  

Priority:  major  Milestone:  sage8.2 
Component:  symbolics  Keywords:  flint, speedup 
Cc:  Merged in:  
Authors:  Ralf Stephan  Reviewers:  
Report Upstream:  N/A  Work issues:  
Branch:  u/rws/16812 (Commits)  Commit:  d7a055df9b360f1f4941046bc99020455f4da09d 
Dependencies:  #24554, #24668  Stopgaps: 
Description
In #16670 the superiority of FLINT for creation of big ChebyshevTpolynomials was confirmed. At T_10000 the speedup is already about 50x versus the present implementation. The issue is outsourced to this ticket.
Change History (50)
comment:1 Changed 6 years ago by
 Branch set to u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation
comment:2 Changed 6 years ago by
 Commit set to e8dd233a1d09f90bfa53ed7542c76879e616e755
 Keywords flint speedup added
 Status changed from new to needs_review
comment:3 Changed 6 years ago by
 Status changed from needs_review to needs_work
Documentation is missing:
def chebyshev_T(unsigned long n, var='x'): """ """
comment:4 Changed 6 years ago by
Also the obvious question: can this be made to work also for chebyshev_U
?
comment:5 followup: ↓ 12 Changed 6 years ago by
I think the conditions n>10 and is_PolynomialRing(P) and P.base() == ZZ and P.ngens() == 1 and x == P.gen()
should be generalized:
 It should work for any
n
, at least negativen
should be supported.
 It should work for any ring of characteristic 0, not just
ZZ
.
 Ideally, it should work for any kind of nonconstant polynomial, not just generators of univariate polynomial rings.
Point 2 and 3 could be supported by first computing the polynomial in ZZ[x]
and then substituting and/or changing the base ring as needed. I'm not saying this is needed on this ticket, but it would be the right thing to do.
comment:6 followup: ↓ 8 Changed 6 years ago by
Why does flint_arith.chebyshev_T create an fmpz_poly, convert it to an array of ZZs, and then convert it back to a polynomial? Over ZZ[x], the roundtrip is probably more expensive than the actual computation. It should be enough to create a new Polynomial_integer_dense_flint instance and apply arith_chebyshev_t to its .__poly
(the n > 10 condition can certainly be dropped).
comment:7 Changed 6 years ago by
 Commit changed from e8dd233a1d09f90bfa53ed7542c76879e616e755 to d89d9ad736fe284182a8fd998dcf0c1e2077657b
comment:8 in reply to: ↑ 6 Changed 6 years ago by
Replying to fredrik.johansson:
...It should be enough to create a new Polynomial_integer_dense_flint instance and apply arith_chebyshev_t to its
.__poly
Well, I tried that but needed to cast polynomial.__poly
to fmpz_poly_t
which is not possible in Cython because the type is an array. OTOH cimport
ing the whole class Polynomial_integer_dense_flint
declaration needs instances of the member functions defined. If I don't do that and just declare
cdef class Polynomial_integer_dense_flint(Polynomial): cdef fmpz_poly_t __poly
the function libs.flint.arith.chebyshev_T
will finally compile but the result is not a rings.polynomial.Polynomial_integer_dense_flint
and returning it and assigning it as such will bomb Sage.
I will now add a member function chebyshev_T
to rings...Polynomial_integer_dense_flint
to avoid all that hassle. My first forage into Cython which certainly presents itself as a fickle customer.
comment:9 Changed 6 years ago by
Putting a method on Polynomial_integer_dense_flint sounds like the easiest way to do it.
comment:10 Changed 6 years ago by
 Commit changed from d89d9ad736fe284182a8fd998dcf0c1e2077657b to 70214f5a36157e3758fcd620709a59bc3f767a09
comment:11 Changed 6 years ago by
 Commit changed from 70214f5a36157e3758fcd620709a59bc3f767a09 to 100622a45e0de7a891ac2cb579c5935571edf49a
Branch pushed to git repo; I updated commit sha1. New commits:
f5837ef  16812: generalize chebyshev_T() by substituting the argument and changing the base ring

e5993da  16812: fix "except Exception: pass"

6f3786b  16812: _eval_numpy_() was not called on numpy objects

100622a  16812: padics doctests sent Sage into coma

comment:12 in reply to: ↑ 5 ; followups: ↓ 15 ↓ 17 Changed 6 years ago by
Replying to jdemeyer:
I think the conditions
n>10 and is_PolynomialRing(P) and P.base() == ZZ and P.ngens() == 1 and x == P.gen()
should be generalized:
 It should work for any
n
, at least negativen
should be supported.
 It should work for any ring of characteristic 0, not just
ZZ
.
 Ideally, it should work for any kind of nonconstant polynomial, not just generators of univariate polynomial rings.
Point 2 and 3 could be supported by first computing the polynomial in
ZZ[x]
and then substituting and/or changing the base ring as needed. I'm not saying this is needed on this ticket, but it would be the right thing to do.
I have implemented this now and see no reason for the restriction in (2) to characteristic 0, as it seems to me that changing the base ring to, e.g., a finite field works nicely too.
In the process I uncovered at least two bugs.
I have removed the padics doctest for the moment because it was evaluated using the new FLINT code and the subsequent substitution was a bit too much. In this case the recursive evaulation is actually faster. But I have no idea at the moment when to select which algorithm, in order to not to fall in this trap.
comment:13 Changed 6 years ago by
 Status changed from needs_work to needs_review
comment:14 Changed 6 years ago by
 Commit changed from 100622a45e0de7a891ac2cb579c5935571edf49a to 61972b5574632e914e178fe211062d0dd0e35b87
comment:15 in reply to: ↑ 12 Changed 6 years ago by
 Status changed from needs_review to needs_work
Replying to rws:
I have removed the padics doctest for the moment because it was evaluated using the new FLINT code and the subsequent substitution was a bit too much. In this case the recursive evaulation is actually faster. But I have no idea at the moment when to select which algorithm, in order to not to fall in this trap.
needs_work for exactly this reason.
comment:16 followup: ↓ 18 Changed 6 years ago by
For the padics, how about this: compute the leading term of the resulting polynomial, which is 2^(n1)
. If that result is zero, use the recursive algorithm. Otherwise, use the FLINT algorithm.
As a general rule, one should never remove doctests to "fix" something, unless you can show that there is a problem with the doctest itself.
comment:17 in reply to: ↑ 12 Changed 6 years ago by
Replying to rws:
I have implemented this now and see no reason for the restriction in (2) to characteristic 0, as it seems to me that changing the base ring to, e.g., a finite field works nicely too.
It would certainly work in theory, but I just don't know if the FLINT algorithm is also faster in this case (I have not tried it).
comment:18 in reply to: ↑ 16 ; followup: ↓ 19 Changed 6 years ago by
Replying to jdemeyer:
As a general rule, one should never remove doctests to "fix" something, unless you can show that there is a problem with the doctest itself.
I think you haven't noticed that I fixed a different doctest, not the padic one. Actually I was aware that the padic doctest should go in again, that's why I wrote "for the moment" earlier. I just wanted to have a version that passes all doctests (I failed to note that another doctest in a different file was bad).
comment:19 in reply to: ↑ 18 Changed 6 years ago by
Replying to rws:
I think you haven't noticed that I fixed a different doctest, not the padic one.
I don't understand what you're trying to say. Fixing one doctest is not an excuse for breaking a different doctest.
Actually I was aware that the padic doctest should go in again, that's why I wrote "for the moment" earlier. I just wanted to have a version that passes all doctests.
OK fine, I understand that. But as long as not all those temporarilyremoved doctests pass, a ticket should remain needs_work.
comment:20 followups: ↓ 21 ↓ 27 Changed 6 years ago by
On another note, the behaviour of the function up to this is unexpectedly complicated and not documented. Given integer n, if the second argument is symbolic then if n<32
a formula is applied else the recursive algorithm. With this ticket, additionally polynomials are handled with FLINT (modulo the refinements discussed).
Now, I have a patch ready for an algorithm
keyword (flint
/recursive
). Also, I want to remove that n<32
bit together with the eval_formula
method that really is no longer necessary now we have FLINT. It also leads to such inconsistencies as
sage: var('n,x') (n, x) sage: chebyshev_T(5,x) 16*x^5  20*x^3 + 5*x sage: chebyshev_T(64, x) 2*(2*(2*(2*(2*(2*x^2  1)^2  1)^2  1)^2  1)^2  1)^2  1
which I would like to have either so or so, not both.
I think I propose to handle a symbolic argument exactly as a polynomial. No one expects such a nested result except if explicitly requested. It's also only faster because it's not expanded.
Do you have any objections?
comment:21 in reply to: ↑ 20 Changed 6 years ago by
Replying to rws:
I think I propose to handle a symbolic argument exactly as a polynomial.
Do you have any objections?
Not at all, your suggestion sounds like a good idea.
comment:22 Changed 6 years ago by
Another bug: special values apply only to integer n
, as T(1/2, 0)
shows, where the formula for x==0
in _eval_special_values_()
gives complex values, while the correct result is 1/sqrt(2)
. So, assumptions have to be checked or the function restricted to n in ZZ
.
Update: changed n==0
to x==0
in the above.
comment:23 followup: ↓ 26 Changed 6 years ago by
I would also say there is a bug in mpmath because it cannot evaluate chebyshev_T(10^6, 0.1)
without throwing a mpmath.NoConvergence
while Sage's recursive algorithm is happy:
sage: chebyshev_T(10^6, 0.1, algorithm="recursive") 0.636384327171504 sage: timeit('chebyshev_T(10^6, 0.1, algorithm="recursive")') 625 loops, best of 3: 229 µs per loop
The result agrees with Wolfram to the 9th place which has 0.636384327247008190135495327481582817644932547700768253402752...
comment:24 Changed 6 years ago by
 Commit changed from 61972b5574632e914e178fe211062d0dd0e35b87 to cc2e34f198961ba4af7b15a925cf893e77898b24
Branch pushed to git repo; I updated commit sha1. New commits:
d21e70c  16812: implement algorithm keyword for eval()

c4b8c17  16812: move examples to top class; doc adaptations

a37b50e  16812: add links from top to classes/methods

e5b19e4  Merge branch 'develop' into t/16812/use_flint_to_speed_up_chebyshev_t_polynomial_creation

9fc5592  16812: we do not want _evalf_() called >100 times...

830ff87  16812: we do not like "except Exception:"

cc2e34f  16812: revamp Chebyshev polys

comment:25 Changed 6 years ago by
This still needs some doctest cleanup.
comment:26 in reply to: ↑ 23 Changed 6 years ago by
Replying to rws:
I would also say there is a bug in mpmath because it cannot evaluate
chebyshev_T(10^6, 0.1)
without throwing ampmath.NoConvergence
while Sage's recursive algorithm is happy:
Yep, mpmath currently doesn't optimize Chebyshev polynomials in any way for integer n (it always just calls 2F1). So there's no reason for Sage to use it for numerical evaluation.
comment:27 in reply to: ↑ 20 ; followup: ↓ 28 Changed 6 years ago by
Replying to rws:
On another note, the behaviour of the function up to this is unexpectedly complicated and not documented. Given integer n, if the second argument is symbolic then if
n<32
a formula is applied else the recursive algorithm. With this ticket, additionally polynomials are handled with FLINT (modulo the refinements discussed).Now, I have a patch ready for an
algorithm
keyword (flint
/recursive
). Also, I want to remove thatn<32
bit together with theeval_formula
method that really is no longer necessary now we have FLINT. It also leads to such inconsistencies assage: var('n,x') (n, x) sage: chebyshev_T(5,x) 16*x^5  20*x^3 + 5*x sage: chebyshev_T(64, x) 2*(2*(2*(2*(2*(2*x^2  1)^2  1)^2  1)^2  1)^2  1)^2  1which I would like to have either so or so, not both.
I think I propose to handle a symbolic argument exactly as a polynomial. No one expects such a nested result except if explicitly requested. It's also only faster because it's not expanded.
Do you have any objections?
Yes I have. There are users (like me) who want the possibility to switch between this representations, especially for large values of n. So please still allow to call the method explicitely.
comment:28 in reply to: ↑ 27 ; followup: ↓ 29 Changed 6 years ago by
Replying to maldun
Yes I have. There are users (like me) who want the possibility to switch between this representations, especially for large values of n. So please still allow to call the method explicitely.
Exactly this is implemented. You can call the method, or you can say algorithm=recursive
.
comment:29 in reply to: ↑ 28 Changed 6 years ago by
Replying to rws:
Replying to maldun
Yes I have. There are users (like me) who want the possibility to switch between this representations, especially for large values of n. So please still allow to call the method explicitely.
Exactly this is implemented. You can call the method, or you can say
algorithm=recursive
.
Good to hear! And thank you for your work! As mentioned in the other ticket: Don't hesitate to ask for help.
comment:30 Changed 6 years ago by
 Commit changed from cc2e34f198961ba4af7b15a925cf893e77898b24 to 90090fd774b972595ce6dccde1ecbb56462fe99d
Branch pushed to git repo; I updated commit sha1. New commits:
90090fd  16812: doctest every execution path in chebyshev_T

comment:31 Changed 6 years ago by
 Status changed from needs_work to needs_review
The FLINT routines are now called on polynomials with exact base ring and used for symbolic input. Everything else gets the recursive treatment which is really fast if you don't have to expand afterwards.
Please review.
comment:32 followup: ↓ 35 Changed 6 years ago by
Typing
chebyshev_T(1000000,x)
in the old version gives a proper result in 7.52 The exact same input gives me full memory (8GB on my machine) and about 1GB swap + a frozen system in the new version. Higher input gives a an error of flint (Memory allocation error)
Some switch would be nice. I guess the best measure would be the break even point on memory consumption.
Edit: I think a good point would be about order 500010000: It is reasonable from view of speed/memory and an explicit representation would be too cluttered anyway
comment:33 Changed 6 years ago by
 Branch changed from u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation to u/maldun/use_flint_to_speed_up_chebyshev_t_polynomial_creation
 Commit changed from 90090fd774b972595ce6dccde1ecbb56462fe99d to 59d04f5fd05174cc0197424363db69c420e3e8f1
 Status changed from needs_review to needs_work
New commits:
59d04f5  Trac 16812: added a switch for evaluation between recursive and flint

comment:34 Changed 6 years ago by
I proposed a patch for that
comment:35 in reply to: ↑ 32 ; followup: ↓ 36 Changed 6 years ago by
Replying to maldun:
Typing
chebyshev_T(1000000,x)in the old version gives a proper result in 7.52 The exact same input gives me full memory (8GB on my machine) and about 1GB swap + a frozen system in the new version.
That result may be proper for a symbolic x
but the times are completely different for ring elements because expansion, which is avoided with symbolic x
, happens automatically. This speed difference for polynomial ring elements actually is the reason for this ticket (see description).
Higher input gives a an error of flint (Memory allocation error)
Some switch would be nice. I guess the best measure would be the break even point on memory consumption.
Edit: I think a good point would be about order 500010000: It is reasonable from view of speed/memory and an explicit representation would be too cluttered anyway
I don't think Sage has to anticipate user memory size, I can go happily to 150.000 with 8GB:
sage: R.<x> = ZZ[] sage: timeit('chebyshev_T(10^5,x)') 5 loops, best of 3: 454 ms per loop sage: timeit('chebyshev_T(150000,x)') 5 loops, best of 3: 1.01 s per loop
Who knows what will be in ten years? Note also that noone displays such results, and that's good because the output process takes near all memory and CPU, not FLINT.
Now, it might have been a bad idea by me to use FLINT on symbolic input but I thought users might appreciate unconvoluted polynomials as output. OTOH, they might get another motivation from this to use rings if they want straight output of polys. So I'll roll back the symbolic involvement, no limit is needed, and Bob's your uncle.
comment:36 in reply to: ↑ 35 Changed 6 years ago by
Replying to rws:
Who knows what will be in ten years? Note also that noone displays such results, and that's good because the output process takes near all memory and CPU, not FLINT.
Now, it might have been a bad idea by me to use FLINT on symbolic input but I thought users might appreciate unconvoluted polynomials as output. OTOH, they might get another motivation from this to use rings if they want straight output of polys. So I'll roll back the symbolic involvement, no limit is needed, and Bob's your uncle.
Than the solution is clear: Make a type check on the input type: If symbolic + higher order use recursion, if Polynomial Ring use flint
It's not about availabilty of RAM: no one wants that a Chebyshev polynomal consumes 1GB of RAM, and no one wants a freezing system, especially if there is a compact represention which only needs the logarithmical amount of time and memory. And also: You don't have much use with an expanded symbolic Cheby poly of order 1e+5, only a cluttered screen, whereas the unexpanded style still is human readable.
Always consider: The main input method should provide some more or less smart default behaviour
On another note: Encourage people to use polynomial ring by putting proper examples in the docstring
comment:37 Changed 6 years ago by
sage: T = chebyshev_T(100000000000,x) Exception (FLINT memory_manager). Unable to allocate memory.  RuntimeError Traceback (most recent call last) <ipythoninput1bef10d0acafb> in <module>() > 1 T = chebyshev_T(Integer(100000000000),x) /home/maldun/sage/sage6.3/local/lib/python2.7/sitepackages/sage/functions/orthogonal_polys.pyc in __call__(self, n, *args, **kwds) 624 # If n is an integer: consider the polynomial as an algebraic (not symbolic) object 625 if n in ZZ and not kwds.get('hold', False): > 626 ret = self._eval_(n, *args, **kwds) 627 if ret is not None: 628 return ret /home/maldun/sage/sage6.3/local/lib/python2.7/sitepackages/sage/functions/orthogonal_polys.pyc in _eval_(self, n, x, *args, **kwds) 679 if ((is_Polynomial(x) and x.parent().base_ring().is_exact()) or 680 (is_Expression(x) and (not x.is_numeric()) and ((n > 10000) or (n < 10000)))): > 681 return self.eval_algebraic(n, x) 682 else: 683 return self.eval_recursive(n, x) /home/maldun/sage/sage6.3/local/lib/python2.7/sitepackages/sage/functions/orthogonal_polys.pyc in eval_algebraic(self, n, arg) 918 n = n 919 R = PolynomialRing(ZZ, 'x') > 920 pol = Polynomial_integer_dense_flint.chebyshev_T(n, R, R.gen()) 921 if is_PolynomialRing(P): 922 pol = pol.subs({pol.parent().gen():arg}) /home/maldun/sage/sage6.3/local/lib/python2.7/sitepackages/sage/rings/polynomial/polynomial_integer_dense_flint.so in sage.rings.polynomial.polynomial_integer_dense_flint.Polynomial_integer_dense_flint.chebyshev_T (build/cythonized/sage/rings/polynomial/polynomial_integer_dense_flint.cpp:14398)() /home/maldun/sage/sage6.3/local/lib/python2.7/sitepackages/sage/ext/c_lib.so in sage.ext.c_lib.sig_raise_exception (build/cythonized/sage/ext/c_lib.c:1040)() RuntimeError: Aborted
Recursive version works without abort
This kills my system:
sage: R.<x> = ZZ[] sage: timeit('chebyshev_T(10^7,x)') so it is not only the use of symbolic expressions
comment:38 Changed 6 years ago by
 Commit changed from 59d04f5fd05174cc0197424363db69c420e3e8f1 to 6165cbbe1ac982bab546c87558d77d7aa0c26b3d
comment:39 Changed 6 years ago by
I propose the following switch mechanism: If the input is symbolic and deg > 10000 return the recursive result else use the regular algorithm
For other input the user has to take responsibility Of course polynomials don't work good with the python recursion, flint does a better job. For symbolic expressions for very high degree the recursion algorithm does a quite good job
comment:40 Changed 6 years ago by
 Branch changed from u/maldun/use_flint_to_speed_up_chebyshev_t_polynomial_creation to u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation
comment:41 Changed 6 years ago by
 Commit changed from 6165cbbe1ac982bab546c87558d77d7aa0c26b3d to 2ee34da7b00e784b8d3be8d9e1f743860d806277
 Status changed from needs_work to needs_review
comment:42 Changed 6 years ago by
 Status changed from needs_review to needs_work
sage t long src/sage/graphs/bipartite_graph.py # 1 doctest failed sage t long src/sage/functions/orthogonal_polys.py # 11 doctests failed
comment:43 Changed 6 years ago by
 Dependencies set to #17531
comment:44 Changed 4 years ago by
It maybe better to just make the Chebys GinacFunctions
and use Flint from within Pynac.
comment:45 Changed 3 years ago by
While other orthogonal polynomial functions (as well as almost all other symbolic functions) are simply BuiltinFunction
s or GinacFunction
s this ticket suffers from the fact that the Cheby functions in Sage allow the algorithm
keyword (by using a __call__
method in the superclass). The way such polymorphisms are resolved in other functions is by having the interface (the callable global function chebyshev_T()
) dispatching on the keyword. This refactoring step (#24554) should be done before all items of this ticket can be implemented.
comment:46 Changed 3 years ago by
 Branch u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation deleted
 Commit 2ee34da7b00e784b8d3be8d9e1f743860d806277 deleted
 Dependencies changed from #17531 to #24554
 Milestone changed from sage6.4 to sage8.2
We will convert the Chebyshevs to GinacFunction
s and call Flint from Pynac, which also will be the fastest way of creating the expression polynomial.
comment:47 Changed 3 years ago by
 Branch set to u/rws/16812
comment:48 Changed 3 years ago by
 Commit set to d7a055df9b360f1f4941046bc99020455f4da09d
 Dependencies changed from #24554 to #24554, pynac0.7.16
New commits:
907aff6  24497: pkg version/chksum

d085587  24497: doctest fixes

810f238  Merge commit 'd085587a0f4319f03a830c39eb5b3d83836f18dd' of git://trac.sagemath.org/sage into HEAD

5c343ab  24554: refactor chebyshev_T

23190a6  24554: refactor chebyshev_U

2582330  24554: remove ChebyshevFunction

96e4e39  Merge branch 'u/rws/remove___call_____usage_in_chebyshev_functions' of git://trac.sagemath.org/sage into t16812

d7a055d  16812: make Chebyshev T/U GinacFunctions

comment:49 Changed 3 years ago by
 Dependencies changed from #24554, pynac0.7.16 to #24554, #24668
 Status changed from needs_work to needs_review
comment:50 Changed 2 years ago by
 Status changed from needs_review to needs_work
Branch sems corrupt.
New commits:
16670: use FLINT for chebyshev_T polynomial
16670: fix previous patch: use argument parent gen
16812: final fix and doctest