Make the orthogonal polynomial eval_* functions ufuncs. #339

merged 15 commits into from Dec 1, 2012


None yet
2 participants

ewmoore commented Oct 17, 2012

As I proposed here and branching from Pauli's work in #333 and #334, I've moved the functions for evaluating orthogonal polynomials to cython and made them into ufuncs. At present there are no changes to how things are actually evaluated. Since this is based on #333 and #334, only ed50214 and on are of interest here.

At present none of the polynomials can be evaluated for a complex number. This probably needs to be fixed since it removes an existing feature. I'll add this soon.

Ticket 1435 still needs to be addressed, but this is easy, too (Need to change

Some other polynomials can also be evaluated using recurrences as long as care is taken, (c.f. This ought to be faster for small n at least, so I'd like to add it where it makes sense.

I'm submitting this as something of a work in progress because wanted to get something up before I ended up just letting it rot on my machine.


ewmoore commented Oct 26, 2012

So I'd like a little input before I spend more time here. I've started filling out various routines using the recurrances directly. However I see two possible issues here.

  1. We are now doing dispatch to different implementations based on the type of the array that is passed. This ends up meaning that for some functions you will possibly get a NaN or not depending on the type. For instance:
In [23]: s.eval_legendre(44000,np.linspace(-1,1,20),sig='ld->d')
array([  1.00000000e+00,   4.59024091e-03,  -1.30797158e-03,
         2.07979703e-03,   2.44349495e-03,  -8.74737533e-04,
        -3.29821075e-04,   1.85696376e-03,  -2.57296402e-03,
        -1.45712334e-04,  -1.45712334e-04,  -2.57296402e-03,
         1.85696376e-03,  -3.29821075e-04,  -8.74737533e-04,
         2.44349495e-03,   2.07979703e-03,  -1.30797158e-03,
         4.59024091e-03,   1.00000000e+00])

In [24]: s.eval_legendre(44000,np.linspace(-1,1,20),sig='dd->d')
/home/ewm/pyenvs/scipy/bin/ipython:1: RuntimeWarning: overflow encountered in eval_legendre
/home/ewm/pyenvs/scipy/bin/ipython:1: RuntimeWarning: invalid value encountered in eval_legendre
array([             nan,              nan,              nan,
                    nan,              nan,              nan,
                    nan,              nan,              nan,
                    nan,              nan,  -2.57296402e-03,
         1.85696376e-03,  -3.29821075e-04,  -8.74737533e-04,
         2.44349495e-03,   2.07979703e-03,  -1.30797158e-03,
         4.59024091e-03,   1.00000000e+00])
  1. The way currently works means that you cannot define a ld->d function and a dd->d function and have things work as expected. Either they end up being defined in an incorrect order (dd->d, ld->d, instead of ld->d, dd->d) or it selects the ld->d version for all of them and rounds the doubles to longs. Is this by design or does this just need to be extended? I'm also curious if I should just wait and see how the possible change from silently rounding doubles to longs to raising an error comes out? This was discussed in a recent pull request, although I'm not sure right off which one that was.




pv commented Oct 26, 2012

The dispatch mechanism indeed needs a tweak, so that the exact signatures (with "i>l" replacement) come first, and all up/downcast versions follow after that. This was the intent, but wasn't required by any of the existing functions, so the current behavior is incorrect.

The outcome from the double -> long downcast issue is probably that we will simply not generate the "l>d" versions of the functions so that an error is raised (so that the error is determined by input type), or we add data value checks in the downcast loops (so that the error is raised by input value). The "i>l" downcast probably needs checking by input value. The issue is that on 64-bit systems it can be that i != l, with l being the main integer type, and this can cause major inconvenience if only i is accepted. The backwards compatibility for the previous unsafe functions is done by adding manually the unsafe versions to the lists.

However, I don't see this issue overlapping very much with this PR.

It might be interesting to see, however, how fast the recurrence-based evaluation is. The advantage of using the existing special function routines is that they tend to have asymptotic expansions which work also when the evaluation would require an absurd number of cycles in the recurrence. One alternative is of course to use the recurrences when they make sense, and fall back to the special functions if not.


ewmoore commented Oct 30, 2012

I've investigated the speed of the recurrence for Legendre polynomials. I've looked at n between 0 and 10,000 for x vectors of 10, 100, 1000, 10,000 elements. The recurrence is almost always faster, except for small vectors at high order.

Some plots: (curves labeled ld use the recurrence, and dd call out to hyp2f1)
dropbox, where the image is hosted seems to be having issues with this one.

This is the same data as the above plot where each curve has been normalized by the length of x.


ewmoore commented Oct 30, 2012

Since Jacobi, Gegenbauer, and Chebyshev polynomials also result in calls to hyp2f1 they'll likely be faster using the recurrence directly as well. I'll check Laguerre and Hermite separately.

@ewmoore ewmoore closed this Oct 30, 2012

@ewmoore ewmoore reopened this Oct 30, 2012


pv commented Oct 31, 2012

Hi, thanks for the results. Now that I think of it again, I fixed our hyp2f1 implementation at one point by adding recurrence evaluation (hyp2f1.c:hyp2f1ra), which is probably equivalent to a generalized ortho-polynomial recurrence. Importantly, there is no code path in hyp2f1 for asymptotic expansions in large orders, so it will never be faster than the recurrence. Which means that orthogonal polynomials with only recurrences can only improve the situation.

Of course, adding asymptotic expansions to hyp2f1 would be possible, but this requires digging up literature, and there are more important things to fix in scipy.special...


pv commented Nov 11, 2012

Note that you can use now also fused types to get the complex-valued versions with the same amount of work. In the ufunc loop description, one can then specify somefunc[double]: d->d, somefunc[double complex]: D->D


pv commented Nov 25, 2012

Ok, this looks good. (Sorry for the delay, apparently Github doesn't send notifications for only pushes, but only for comments).


pv commented Nov 25, 2012

Rebased branch here:
I also converted the complex loops to use fused types.

To work on that:

git remote add pv git://
git fetch pv
git checkout ortho_eval
git branch ortho_eval_backup ortho_eval
git reset --hard pv/ortho_eval   # NB. resets also working tree

Apparently, the long-loops were not run before due to the loop selection issues in They are now after the rebase, and the tests currently seem to indicate that the Jacobi and Gegenbauer polynomial recurrences apparently give wrong results.

Also, some new tests need to be added to compare the double and long loops to each other.


ewmoore commented Nov 28, 2012

Ok, I've pulled down your changes and have started looking at the failures.

I'm using forms of the recurrences given here: . Which may or may not be necessary. One of things they do is to work with a rescaled recurrence. That rescaling for the Gegenbauer polynomials rescales the value at x = 1.0 to 1.0, using C(n, a, 1.0) = binom(n+2a-1, n). However, this requires that the binomial coefficient, binom(n,k) be computed for non-integer values of n. This doesn't seem like a big deal, but it gets tricky. We can extend the binomial coefficient using either of the gamma or beta function, but we are calculating it using gammaln and exp, which is giving us the wrong sign when the binomial coefficient should be negative. Consider:

In [1]: import numpy as np; import scipy.special as s;

In [2]: n = 0.264

In [3]: k = 2

In [4]: s.binom(n,k)
Out[4]: 0.097152000000000016

In [5]: np.exp(s.gammaln(n+1) - s.gammaln(k+1) - s.gammaln(n-k+1))
Out[5]: 0.097152000000000016

In [6]: s.gamma(n+1) / ( s.gamma(k+1) * s.gamma(n-k+1) )
Out[6]: -0.097152000000000002

In [7]: 1.0 / ((n+1) * s.beta(n-k+1, k+1))
Out[7]: -0.09715200000000003

Computing the binomal coefficient directly with gamma or beta fixes the Gegenbauer failure, I'll push something tomorrow.. It may well be something similar for Jacobi, I'll have to see.

It seems to me that the current binom implementation gives the wrong answer at times.


pv commented Nov 28, 2012

In [5]: np.exp(s.gammaln(n+1) - s.gammaln(k+1) - s.gammaln(n-k+1))
Out[5]: 0.097152000000000016

In [6]: s.gamma(n+1) / ( s.gamma(k+1) * s.gamma(n-k+1) )
Out[6]: -0.097152000000000002

Indeed, gammaln gives the logarithm of the absolute value. We'd probably
want an additional function in scipy.special giving the sign of the
Gamma function.

On the real axis,

sgn Gamma(x)  = (-1)**(floor(x)) , x < 0
sgn Gamma(x)  = 1, x > 0

The sign at non-positive integers is undefined, but I think it should be
chosen as a finite number rather than nan because we want

1 / Gamma(x) = gammasgn(x) * exp(-gammaln(x))

be an entire function. gammasgn(-n) = 1 is probably not a bad choice,
and gammasgn(-n) = 0 would have the advantage that for it,
gammasgn(-n) * exp(gammaln(-n)) = nan rather than inf.

I'd avoid the expression in terms of plain gamma functions --- that will
start spitting out nans for large n even when the binomial coefficient
is not large.

The expression with the Beta function is probably OK, although may need
an additional power series expansion (very) close to the limit k->0 (or
at least dealing with the point k=0).

Related: we'd also need a new implementation of the Pochhammer symbol,
Pochhammer(a, n) = Gamma(a + n) / Gamma(a).



ewmoore commented Nov 29, 2012

I can't seem to push my new changes.

(scipy)ewm@hollow:/scipy$ git push
Username for '': ewmoore
Password for '':
! [rejected] ortho_eval -> ortho_eval (non-fast-forward)
error: failed to push some refs to ''
hint: Updates were rejected because the tip of your current branch is behind
hint: its remote counterpart. Merge the remote changes (e.g. 'git pull')
hint: before pushing again.
hint: See the 'Note about fast-forwards' in 'git push --help' for details.

Any ideas?


pv commented Nov 29, 2012

git push -f -- rebased commits are not on top of the old ones, so it doesn't let you push them without confirmation.


ewmoore commented Nov 30, 2012

Okay, I've fixed up the issues with eval_gegenbauer and eval_jacobi, added a gammasgn function with test, fixed the bug in binom, add a test and add a test between the ld->d and dd->d version.

The only thing I know of that might still be outstanding here is that we don't have a dd->d version for either form of the Hermite polynomials. There are a bunch of relationships that might work, but at this point we haven't lost any abilities we had before I started.


ewmoore commented Nov 30, 2012

Oh and I see that lgam actually provides the sign in cephes, but it does so rather unpleasantly by placing it in a global variable named sgngam. This isn't AFAICT exposed anywhere. I don't think that we could actually use sgngam in binom without moving that to cephes. I also don't really see a backward compatible way to add the sign to the exposed gammaln function.

pv added a commit that referenced this pull request Dec 1, 2012

Merge pull request #339 from ewmoore/ortho_eval
ENH: special: Make the orthogonal polynomial eval_* functions ufuncs

Convert the orthogonal polynomial eval_* functions to real ufuncs, and
implement their evaluation via recurrences where applicable.  This makes e.g.
their out= arguments work properly, and makes the simple cases faster.

Also add new gammasgn() and binom() functions.

@pv pv merged commit e1ae0a0 into scipy:master Dec 1, 2012


pv commented Dec 1, 2012

Thanks, merged. The general Hermite polynomials can go in a separate PR if needed.

I opted in the end to go for the beta function relation in binom(), because the case n >> k needs special handling to avoid loss of precision. Unfortunately, this isn't yet implemented in the beta() function either, but that's a separate bug.


pv commented Dec 1, 2012

Hmm, lots of failures turn up on mingw-windows build:

Not immediately clear what's going on --- some routines seem to be returning nans.

Some syntax unsupported by earlier Python versions was also present:


pv commented Dec 2, 2012

facd638 fixed the nan issues.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment