Inverting a power series

Inverting a power series

As a student, one of the things that seemed curious to me about power series was that a function might have a simple power series, but its inverse could be much more complicated. For example, the coefficients in the series for arctangent have a very simple pattern

but the coefficients in the series for tangent are mysterious.

There’s no obvious pattern to the coefficients. It is possible to write the sum in closed form, but this requires introducing the Bernoulli numbers which are only slightly less mysterious than the power series for tangent.

To a calculus student, this is bad news: a simple, familiar function has a complicated power series. But this is good news for combinatorics. Reading the equation from right to left, it says a complicated sequence has a simple generating function!

Computing the coefficients

The example above suggests that inverting a power series, i.e. starting with the power series for a function and computing the power series for its inverse, is going to be complicated. We introduce exponential Bell polynomials to encapsulate some of the complexity, analogous to introducing Bernoulli numbers above.

Assume the function we want to invert, f(x), satisfies f(0) = 0 and its derivative satisfies f‘(0) ≠ 0. Assume f and its inverse g have the series representations

and

Note that this isn’t the power series per se. The coefficients here are k! times the ordinary power series coefficients. (Compare with ordinary and exponential generating functions explained here.)

Then you can compute the series coefficients for g from the coefficients for f as follows. We have g₁ = 1/f₁ (things start out easy!) and for n ≥ 2,

where the B‘s are exponential Bell polynomials,

and

is the kth rising power of n. (More on rising and falling powers here.)

Example

Let’s do an example in Python. Suppose we want a series for the inverse of the gamma function near 2. The equations above assume we’re working in a neighborhood of 0 and that our function is 0 at 0. So we will invert the series for f(x) = Γ(x+2) − 1 and then adjust the result to find the inverse of Γ(x).

import numpy as np
from scipy.special import gamma
from sympy import factorial, bell 

def rising_power(a, k):
    return gamma(a+k)/gamma(a)

# Taylor series coefficients for gamma centered at 2
gamma_taylor_coeff = [
    1.00000000000000,
    0.42278433509846,
    0.41184033042643,
    0.08157691924708,
    0.07424901075351,
    -0.0002669820687,
    0.01115404571813,
    -0.0028526458211,
    0.00210393334069,
    -0.0009195738388,
    0.00049038845082,    
]
N = len(gamma_taylor_coeff)

f_exp_coeff = np.zeros(N)
for i in range(1, N):
    f_exp_coeff[i] = gamma_taylor_coeff[i]*factorial(i)

# Verify the theoretical requirements on f
assert( f_exp_coeff[0] == 0 )
assert( f_exp_coeff[1] != 0 )

f_hat = np.zeros(N-1)
for k in range(1, N-1):
    f_hat[k] = f_exp_coeff[k+1] / ((k+1)*f_exp_coeff[1])
    
g_exp_coeff = np.zeros(N)
g_exp_coeff[1] = 1/f_exp_coeff[1]

for n in range(2, N):
    s = 0
    for k in range(1, n):
        # Note that the last argument to bell must be a list,
        # and not a NumPy array.
        b = bell(n-1, k, list(f_hat[1:(n-k+1)]))
        s += (-1)**k * rising_power(n, k) * b
    g_exp_coeff[n] = s / f_exp_coeff[1]**n

def g(x):
    return sum(
        x**i * g_exp_coeff[i] / factorial(i) for i in range(1, N)
    )

def gamma_inverse(x):
    return 2. + g(x-1.)

# Verify you end up where you started when you do a round trip
for x in [1.9, 2.0, 2.1]:
    print( gamma_inverse(gamma(x)) )

Output:

1.90000003424331
2.00000000000000
2.09999984671755