Differences ψₙ(x + m) − ψₙ(x) for integer m: subtraction loses up to 9 digits for large x, while the exact finite sum is accurate and about 20× cheaper

Open
  • Is something wrong?
#569 2 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 1 day

Nobody has claimed this yet.

Assessment

Difficulty
5/5
Estimated time
Over a week
Newbie friendliness
40/100
Issue type
Feature
Clarity
Mostly clear
Activity status
Active
Tech stack
julia

Research direction

Start with src/gamma.jl, which the issue identifies as containing the current special-function implementation, and read the existing loggammadiv API and issue #93 for related design. The proposal leaves the public API shape open, so first determine how integer-shift differences should fit the package. Done means an agreed API that computes the stated difference accurately and has tests covering the examples in the reproducer.

Written by the indexing model from the issue text.

Description

polygamma(n, x + m) - polygamma(n, x) for a nonnegative integer m comes up whenever the score or the curvature of a count likelihood is taken with respect to a shape parameter: the negative binomial's ψ(r + k) − ψ(r) and ψ′(r + k) − ψ′(r) for a count k, and the same shapes in beta-binomial and Dirichlet-multinomial models (the loggamma(x + α) − loggamma(α) analogue is what Distributions.jl writes out for the Dirichlet-multinomial). Taking the difference by subtraction cancels badly once x is large (the Poisson limit of the negative binomial, large concentration parameters), while the identity

ψₙ(x + m) − ψₙ(x) = (−1)ⁿ n! Σₖ₌₀^{m−1} (x + k)^−(n+1)

gives the difference exactly for small m, with no cancellation and no special-function call.

Reproducer

using SpecialFunctions
# ψₙ(x + m) − ψₙ(x) for integer m ≥ 0, by subtraction and as the exact finite sum
# ψₙ(x + m) − ψₙ(x) = (−1)ⁿ n! Σₖ₌₀^{m−1} (x + k)^−(n+1)
subtract(n, x, m) = polygamma(n, x + m) - polygamma(n, x)
finitesum(n, x, m) = (-1)^n * factorial(n) * sum(k -> (x + k)^-(n + 1), 0:m-1; init=zero(x))
reference(n, x, m) = setprecision(BigFloat, 512) do
    finitesum(n, big(x), m)
end
relerr(a, ref) = Float64(abs((big(a) - ref) / ref))
for x in (0.5, 10.0, 1e3, 1e6, 1e9), m in (1, 3, 50)
    print("x = ", lpad(x, 6), "  m = ", lpad(m, 2))
    for n in 0:2
        ref = reference(n, x, m)
        print("   n=$n: subtraction ", round(relerr(subtract(n, x, m), ref); sigdigits=2),
            "  sum ", round(relerr(finitesum(n, x, m), ref); sigdigits=2))
    end
    println()
end

Output on Julia 1.10.12 with SpecialFunctions 2.9.1; master at d8ffcc1b456171884112a157b7be1bc4707f77be has a byte-identical src/gamma.jl, so it behaves the same:

x =    0.5  m =  1   n=0: subtraction 2.2e-16  sum 0.0   n=1: subtraction 4.4e-16  sum 0.0   n=2: subtraction 2.2e-16  sum 0.0
x =    0.5  m =  3   n=0: subtraction 6.8e-17  sum 7.7e-17   n=1: subtraction 2.7e-16  sum 7.4e-17   n=2: subtraction 2.6e-16  sum 4.8e-17
x =    0.5  m = 50   n=0: subtraction 1.4e-17  sum 2.9e-16   n=1: subtraction 3.6e-16  sum 1.8e-16   n=2: subtraction 2.0e-16  sum 1.4e-17
x =   10.0  m =  1   n=0: subtraction 8.9e-16  sum 5.6e-17   n=1: subtraction 5.0e-16  sum 1.9e-16   n=2: subtraction 2.6e-15  sum 2.1e-17
x =   10.0  m =  3   n=0: subtraction 1.1e-15  sum 3.8e-17   n=1: subtraction 6.7e-17  sum 7.1e-17   n=2: subtraction 8.5e-16  sum 1.1e-16
x =   10.0  m = 50   n=0: subtraction 5.0e-16  sum 3.8e-16   n=1: subtraction 6.9e-17  sum 2.3e-16   n=2: subtraction 4.3e-16  sum 2.7e-16
x = 1000.0  m =  1   n=0: subtraction 3.3e-13  sum 2.1e-17   n=1: subtraction 3.5e-13  sum 4.5e-17   n=2: subtraction 1.4e-13  sum 6.2e-17
x = 1000.0  m =  3   n=0: subtraction 3.3e-13  sum 2.1e-17   n=1: subtraction 5.3e-14  sum 1.2e-16   n=2: subtraction 3.9e-14  sum 7.8e-18
x = 1000.0  m = 50   n=0: subtraction 4.9e-15  sum 4.0e-17   n=1: subtraction 4.9e-15  sum 2.4e-16   n=2: subtraction 3.6e-15  sum 1.4e-16
x =  1.0e6  m =  1   n=0: subtraction 1.0e-9  sum 4.5e-17   n=1: subtraction 1.8e-10  sum 2.0e-17   n=2: subtraction 1.3e-10  sum 7.2e-17
x =  1.0e6  m =  3   n=0: subtraction 3.5e-10  sum 8.3e-17   n=1: subtraction 3.4e-11  sum 8.9e-17   n=2: subtraction 4.9e-11  sum 3.3e-17
x =  1.0e6  m = 50   n=0: subtraction 4.6e-12  sum 1.5e-16   n=1: subtraction 5.4e-12  sum 2.5e-16   n=2: subtraction 1.9e-12  sum 4.4e-17
x =  1.0e9  m =  1   n=0: subtraction 8.3e-8  sum 6.2e-17   n=1: subtraction 1.5e-7  sum 7.2e-17   n=2: subtraction 1.2e-7  sum 3.8e-17
x =  1.0e9  m =  3   n=0: subtraction 8.4e-8  sum 2.0e-17   n=1: subtraction 1.3e-8  sum 1.2e-16   n=2: subtraction 1.9e-9  sum 4.9e-17
x =  1.0e9  m = 50   n=0: subtraction 3.6e-8  sum 4.3e-17   n=1: subtraction 3.8e-9  sum 7.6e-17   n=2: subtraction 6.0e-10  sum 1.9e-16

These are relative errors against a 512-bit reference. Subtraction loses about 3 digits at x = 10³, 7 at x = 10⁶ and 9 at x = 10⁹ (for m = 1: 8.3e-8, 1.5e-7 and 1.2e-7 at n = 0, 1, 2), while the finite sum stays within a few ulp everywhere. This is not a defect of digamma/trigamma/polygamma themselves, which are accurate to a few ulp (#450, #451): it is the subtraction of two nearly equal values.

Cost

For m up to 8, the three sums together cost about 7 ns per evaluation against about 170 ns for the six special-function calls they replace (Julia 1.10.12, x in [0.5, 3.5], measured on a busy machine, so take the ratio as rough: roughly 20×).

Question

Would a function for this difference be in scope? For example

polygamma_shift(n, x, m)   # ψₙ(x + m) − ψₙ(x), the exact finite sum for integer m ≥ 0

possibly with digamma_shift / trigamma_shift conveniences, or a two-point form ψₙ(y) − ψₙ(x) that uses the finite sum when y − x is a small integer and the difference of the asymptotic series otherwise. The package already has loggammadiv(a, b) for log(Γ(b) / Γ(a + b)), which exists for the same reason, and the rising factorial of #93 is the loggamma counterpart of this sum.

If a different shape fits the package better, I'm happy to adjust the proposal.

Dominant language
Julia
Stars
384
Forks
115
Avg merge
1d 13h
Merged PRs (30d)
5

Getting set up

This project ships no dev container, Dockerfile or contributing guide, so setting up is up to you: start from its README, and see our first-contribution guide for the general steps.

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from JuliaMath/SpecialFunctions.jl

All issues in JuliaMath/SpecialFunctions.jl

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.