Opened 8 years ago
Last modified 8 years ago
#17066 needs_work defect
always simplify hypergeometric() when it's a polynomial
Reported by:  Ralf Stephan  Owned by:  

Priority:  major  Milestone:  sage6.4 
Component:  symbolics  Keywords:  special, functions, pFq, evaluation 
Cc:  Merged in:  
Authors:  Ralf Stephan  Reviewers:  
Report Upstream:  N/A  Work issues:  
Branch:  u/rws/work_around_mpmath_problem_with_hypergeometric___zeroes (Commits, GitHub, GitLab)  Commit:  0e0fef65a021a9b0a9db30a2bc6d300ede68bf14 
Dependencies:  Stopgaps: 
Description
See https://groups.google.com/forum/?hl=en#!topic/sagedevel/8z2HyH5Y8Dc with Fredrik's assessment:
"For Sage, fixing the problem is actually trivial: when the hypergeometric function is a polynomial (and at least when the inputs are exact), don't call mpmath; just evaluate the polynomial directly and then call .n() on the result."
Change History (26)
comment:1 Changed 8 years ago by
Branch:  → u/rws/work_around_mpmath_problem_with_hypergeometric___zeroes 

comment:2 Changed 8 years ago by
Authors:  → Ralf Stephan 

Commit:  → b6f0fbc3983b520c28330e49ec38d79aadf1796f 
Status:  new → needs_review 
Summary:  work around mpmath problem with hypergeometric() zeroes → always try to simplify hypergeometric() 
comment:3 Changed 8 years ago by
Status:  needs_review → needs_work 

Maybe document the problem at ask.sagemath is fixed, too?
hypergeometric([2,1],[2],1).n(100)
was the example Dima posted on sagedevel...
comment:4 Changed 8 years ago by
The "solution" is tongueincheek, anyway, at least its application should be restricted in order to not burn CPU needlessly with higher orders of the function.
comment:6 followup: 8 Changed 8 years ago by
This is related to the following ask.sagemath question as well.
comment:7 Changed 8 years ago by
Commit:  b6f0fbc3983b520c28330e49ec38d79aadf1796f → dff348c1c81579de07e88485e05d5e40b0d31942 

comment:8 Changed 8 years ago by
Replying to kcrisman:
This is related to the following ask.sagemath question as well.
Indeed. While what you're saying there about decidability is certainly right, we can bit off a good chunk using known closed forms. I have now fixed a bug in closed_form()
and it looks like using that Sage function (provided by Fredrik with the hypergeometric code) instead of Maxima is not only faster but gives also better results:
sage: hypergeometric([3/2],[1/2,1/2],5) # using closed_form() 2*sqrt(5)*sinh(2*sqrt(5))  19*cosh(2*sqrt(5)) sage: from sage.interfaces.maxima_lib import maxima_lib, max_to_sr sage: max_to_sr(maxima_lib.hgfred([3/2],[1/2,1/2],5).ecl()) 4*(5*5^(1/4))*sqrt(1/5)*sqrt(pi)*sqrt(sqrt(5)/pi)*cosh(2*sqrt(5))  2*5^(3/4)*sqrt(pi)*bessel_I(3/2, 2*sqrt(5))
The latter, while numerically correct, is difficult to simplify further.
There are several doctests failing where I need to decide how to fix or rewrite them.
comment:9 Changed 8 years ago by
Commit:  dff348c1c81579de07e88485e05d5e40b0d31942 → a2253833d3edd9cd6fa9ccfc9509592ad2b9fad4 

Branch pushed to git repo; I updated commit sha1. New commits:
a225383  17066: restrict to nonpositive first arg; fix doctests

comment:10 Changed 8 years ago by
Commit:  a2253833d3edd9cd6fa9ccfc9509592ad2b9fad4 → 66ed0c0586f2ec1b78a28e76a1f34f258b6f4f96 

Branch pushed to git repo; I updated commit sha1. New commits:
66ed0c0  17066: more doctests

comment:11 Changed 8 years ago by
Status:  needs_work → needs_review 

Had to restrict the simplification, after all, as advised.
comment:12 Changed 8 years ago by
Yeah, make sure nothing is wrong!
I know that I won't be able to finish the review of this immediately. But I do notice that all the hold=True
will be quite confusing in the doctests. I recommend that we add some more doctests that test the methods in question with nonsimplifying stuff, and then make it clear to the reader why we are all of a sudden using hold=True
. On the plus side, hold=True
seems to have worked great!
comment:13 Changed 8 years ago by
Commit:  66ed0c0586f2ec1b78a28e76a1f34f258b6f4f96 → 55daf81080f46c54429e629b65e8049acfdb6d00 

Branch pushed to git repo; I updated commit sha1. New commits:
55daf81  17066: explain why we add hold=True to some doctests

comment:14 Changed 8 years ago by
Summary:  always try to simplify hypergeometric() → always simplify hypergeometric() when it's a polynomial 

comment:16 followups: 17 23 Changed 8 years ago by
Expanding hypergeometric functions in polynomial closed form before evaluating them can lead to numerical cancellation. Compare for example
sage: hypergeometric([1,100],[3], 0.7) 0.0278923450568346 sage: hypergeometric([1,100],[3], Reals(1000)(0.7)).n() 0.0278923450568346 sage: sage.functions.hypergeometric.closed_form(hypergeometric([1,100],[3],x)).subs(x=0.7) 376.859920391941
Won't that be a problem with this ticket? (I didn't read beyond the ticket's description, so please ignore my comment if it is beside the point.)
comment:17 followup: 19 Changed 8 years ago by
Replying to mmezzarobba:
Won't that be a problem with this ticket?
It appears the only solutions are A) restrict the range of parameters where ticket is applied B) do not simplify; just fix the evaluation problem C) rely on the user to supply needed precision; add warning
EDIT: nonsense about n()
removed
comment:18 Changed 8 years ago by
When it's a polynomial of low degree and the inputs are not too large, you could do everything with rational arithmetic and convert to a correctly rounded floatingpoint number in the end. Or you could use interval arithmetic.
comment:19 Changed 8 years ago by
Replying to rws:
It appears the only solutions are A) restrict the range of parameters where ticket is applied B) do not simplify; just fix the evaluation problem C) rely on the user to supply needed precision; add warning
If (as you seem to say in your post to sagedevel
) the result of the ticket is that
sage: hypergeometric([1,100],[3],x)
returns a polynomial while
sage: hypergeometric([1,100],[3],0.7)
still calls a numerically stable evaluation algorithm, then I don't think there is a problem. (Ideally, hypergeometric([1,100],[3],x,hold=True).subs(x=0.7)
should also call the numerically stable evaluation algorithm, though.)
What I was concerned about was that the implementation of n()
for an "atomiclooking" expression with dedicated numerical evaluation code would end up using a numerically unstable formula internally without adapting the precision or warning the user. I'm not saying such a behaviour is an absolute nogo, but I do think it is something we should avoid when possible.
comment:20 Changed 8 years ago by
It's not that I don't see your point. The real fix for this cannot be provided by hypergeometric
though but must be done in subs
, if at all.
comment:21 Changed 8 years ago by
Since mpmath uses a numerically stable algorithm but fails on zero, you could potentially also try calling mpmath, catch the exception, and then check symbolically if it's exactly zero.
comment:22 Changed 8 years ago by
Commit:  55daf81080f46c54429e629b65e8049acfdb6d00 → 0e0fef65a021a9b0a9db30a2bc6d300ede68bf14 

comment:23 followup: 24 Changed 8 years ago by
Expanding hypergeometric functions in polynomial closed form before evaluating them can lead to numerical cancellation. Compare for example
sage: hypergeometric([1,100],[3], 0.7) 0.0278923450568346 sage: hypergeometric([1,100],[3], Reals(1000)(0.7)).n() 0.0278923450568346 sage: sage.functions.hypergeometric.closed_form(hypergeometric([1,100],[3],x)).subs(x=0.7) 376.859920391941Won't that be a problem with this ticket? (I didn't read beyond the ticket's description, so please ignore my comment if it is beside the point.)
Just so I follow, what you are saying is that the polynomial (rational function, I guess) given will be correct, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision? I'm torn here because on the one hand, presumably that could happen with any rational function, there are standard examples given in a lot of calculus books, but on the other hand we "expect" hg to return a hg function and not a polynomial.
The real fix for this cannot be provided by hypergeometric though but must be done in subs, if at all.
Well, it could provide a good reason not to "always simplify hypergeometric when it's a polynomial"... I didn't realize that people meant rational function when they said polynomial. Is it possible to really only simplify when it's a true polynomial, or is that a very unlikely event?
comment:24 followup: 25 Changed 8 years ago by
Replying to kcrisman:
Just so I follow, what you are saying is that the polynomial (rational function, I guess)
No, it is really a polynomial—but one of degree 100, with large coefficients that alternate in sign.
given will be correct, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision?
Precisely. (Well, not necessarily a very high precision, but probably something like 80100 guard bits in my example.)
comment:25 Changed 8 years ago by
Just so I follow, what you are saying is that the polynomial (rational function, I guess)
No, it is really a polynomial—but one of degree 100, with large coefficients that alternate in sign.
Oh, I misread some of the output in the example I tested  hard to see where the dividing signs go :) So it's the alternation that is the cancellation problem, got it.
given will be correct, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision?
Precisely. (Well, not necessarily a very high precision, but probably something like 80100 guard bits in my example.)
Anyway, 53 default yields madness, as we see in comment:16. Not good, but I'm not sure what the "right" answer is; I'm leaning toward saying we should not do this automatically, then, but making the simplification a lot easier to find (which is probably already done in this ticket). There are likely as many people using hg functions for numerics as for symbolics, right?
comment:26 Changed 8 years ago by
Status:  needs_review → needs_work 

sage t long src/sage/symbolic/expression.pyx # 4 doctests failed sage t long src/sage/functions/bessel.py # 1 doctest failed sage t long src/sage/functions/hypergeometric.py # 102 doctests failed
This patch just changes
hypergeometric()
to always invokeMaxima_lib.hgfred()
.New commits:
17066: always try to simplify