Opened 8 years ago

use FLINT to speed up Chebyshev T polynomial creation

Reported by: Owned by: rws major sage-8.2 symbolics flint, speedup Ralf Stephan N/A u/rws/16812 d7a055df9b360f1f4941046bc99020455f4da09d #24554, #24668

Description

In #16670 the superiority of FLINT for creation of big Chebyshev-T-polynomials was confirmed. At T_10000 the speedup is already about 50x versus the present implementation. The issue is outsourced to this ticket.

comment:1 Changed 8 years ago by rws

• Branch set to u/rws/use_flint_to_speed_up_chebyshev_t_polynomial_creation

comment:2 Changed 8 years ago by rws

• Commit set to e8dd233a1d09f90bfa53ed7542c76879e616e755
• Status changed from new to needs_review

New commits:

 ​ba4e34d 16670: use FLINT for chebyshev_T polynomial ​df781ab 16670: fix previous patch: use argument parent gen ​e8dd233 16812: final fix and doctest

comment:3 Changed 8 years ago by jdemeyer

• Status changed from needs_review to needs_work

Documentation is missing:

def chebyshev_T(unsigned long n, var='x'):
"""
"""

comment:4 Changed 8 years ago by jdemeyer

Also the obvious question: can this be made to work also for chebyshev_U?

comment:5 follow-up: ↓ 12 Changed 8 years ago by 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:

1. It should work for any n, at least negative n should be supported.
1. It should work for any ring of characteristic 0, not just ZZ.
1. Ideally, it should work for any kind of non-constant 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 follow-up: ↓ 8 Changed 8 years ago by fredrik.johansson

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

Last edited 8 years ago by fredrik.johansson (previous) (diff)

comment:7 Changed 8 years ago by git

• Commit changed from e8dd233a1d09f90bfa53ed7542c76879e616e755 to d89d9ad736fe284182a8fd998dcf0c1e2077657b

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

comment:8 in reply to: ↑ 6 Changed 8 years ago by rws

...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 cimporting 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 8 years ago by fredrik.johansson

Putting a method on Polynomial_integer_dense_flint sounds like the easiest way to do it.

comment:10 Changed 8 years ago by git

• Commit changed from d89d9ad736fe284182a8fd998dcf0c1e2077657b to 70214f5a36157e3758fcd620709a59bc3f767a09

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

 ​85a141b 16812: move chebyshev_t into class Polynomial_integer_dense_flint ​70214f5 16812: remove restriction on n

comment:11 Changed 8 years ago by git

• 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 ; follow-ups: ↓ 15 ↓ 17 Changed 8 years ago by rws

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:

1. It should work for any n, at least negative n should be supported.
1. It should work for any ring of characteristic 0, not just ZZ.
1. Ideally, it should work for any kind of non-constant 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 8 years ago by rws

• Status changed from needs_work to needs_review

comment:14 Changed 8 years ago by git

• Commit changed from 100622a45e0de7a891ac2cb579c5935571edf49a to 61972b5574632e914e178fe211062d0dd0e35b87

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

 ​b032291 Merge branch 'develop' into t/16812/use_flint_to_speed_up_chebyshev_t_polynomial_creation ​61972b5 16812: fix doctest

comment:15 in reply to: ↑ 12 Changed 8 years ago by jdemeyer

• Status changed from needs_review to needs_work

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 follow-up: ↓ 18 Changed 8 years ago by jdemeyer

For the p-adics, how about this: compute the leading term of the resulting polynomial, which is 2^(n-1). 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 8 years ago by jdemeyer

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 ; follow-up: ↓ 19 Changed 8 years ago by rws

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 p-adic one. Actually I was aware that the p-adic 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 8 years ago by jdemeyer

I think you haven't noticed that I fixed a different doctest, not the p-adic 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 p-adic 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 temporarily-removed doctests pass, a ticket should remain needs_work.

comment:20 follow-ups: ↓ 21 ↓ 27 Changed 8 years ago by 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 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 8 years ago by jdemeyer

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 8 years ago by rws

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.

Last edited 8 years ago by rws (previous) (diff)

comment:23 follow-up: ↓ 26 Changed 8 years ago by rws

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

• 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 8 years ago by rws

This still needs some doctest cleanup.

comment:26 in reply to: ↑ 23 Changed 8 years ago by fredrik.johansson

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:

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 ; follow-up: ↓ 28 Changed 8 years ago by maldun

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?

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.

Last edited 8 years ago by maldun (previous) (diff)

comment:28 in reply to: ↑ 27 ; follow-up: ↓ 29 Changed 8 years ago by rws

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

• 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 8 years ago by rws

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

comment:32 follow-up: ↓ 35 Changed 8 years ago by 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. 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 5000-10000: It is reasonable from view of speed/memory and an explicit representation would be too cluttered anyway

Last edited 8 years ago by maldun (previous) (diff)

comment:33 Changed 8 years ago by maldun

• 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 8 years ago by maldun

I proposed a patch for that

comment:35 in reply to: ↑ 32 ; follow-up: ↓ 36 Changed 8 years ago by rws

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 5000-10000: 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 8 years ago by maldun

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

Last edited 8 years ago by maldun (previous) (diff)

comment:37 Changed 8 years ago by maldun

sage: T = chebyshev_T(100000000000,x)
Exception (FLINT memory_manager). Unable to allocate memory.
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-1-bef10d0acafb> in <module>()
----> 1 T = chebyshev_T(Integer(100000000000),x)

/home/maldun/sage/sage-6.3/local/lib/python2.7/site-packages/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/sage-6.3/local/lib/python2.7/site-packages/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/sage-6.3/local/lib/python2.7/site-packages/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/sage-6.3/local/lib/python2.7/site-packages/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/sage-6.3/local/lib/python2.7/site-packages/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
Last edited 8 years ago by maldun (previous) (diff)

comment:38 Changed 8 years ago by git

• Commit changed from 59d04f5fd05174cc0197424363db69c420e3e8f1 to 6165cbbe1ac982bab546c87558d77d7aa0c26b3d

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

 ​e5c0a83 Trac 16812: Make switch in default evaluation method ​9e6001c Merge branch 'cheby3' into cheby ​f02d699 Trac 16812: removed non default switch mechanism ​6165cbb Trac 16812: Make another proposal for a switch mechanism

comment:39 Changed 8 years ago by maldun

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 8 years ago by rws

• 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 8 years ago by rws

• Commit changed from 6165cbbe1ac982bab546c87558d77d7aa0c26b3d to 2ee34da7b00e784b8d3be8d9e1f743860d806277
• Status changed from needs_work to needs_review

I fixed your doctest. Let's see what the reviewer/s say.

New commits:

 ​2ee34da 16812: fix doctest

comment:42 Changed 7 years ago by rws

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

• Dependencies set to #17531

comment:44 Changed 6 years ago by rws

It maybe better to just make the Chebys GinacFunctions and use Flint from within Pynac.

comment:45 Changed 5 years ago by rws

While other orthogonal polynomial functions (as well as almost all other symbolic functions) are simply BuiltinFunctions or GinacFunctions 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 5 years ago by rws

• 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 sage-6.4 to sage-8.2

We will convert the Chebyshevs to GinacFunctions and call Flint from Pynac, which also will be the fastest way of creating the expression polynomial.

comment:47 Changed 5 years ago by rws

• Branch set to u/rws/16812

comment:48 Changed 5 years ago by rws

• Commit set to d7a055df9b360f1f4941046bc99020455f4da09d
• Dependencies changed from #24554 to #24554, pynac-0.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 4 years ago by rws

• Dependencies changed from #24554, pynac-0.7.16 to #24554, #24668
• Status changed from needs_work to needs_review

comment:50 Changed 4 years ago by rws

• Status changed from needs_review to needs_work

Branch sems corrupt.

Note: See TracTickets for help on using tickets.