Over a million developers have joined DZone.
{{announcement.body}}
{{announcement.title}}

How to Compute the Power Series for an Inverse Function

DZone's Guide to

How to Compute the Power Series for an Inverse Function

To a calculus student, a simple, familiar function has a complicated power series is bad news. But this is good news for combinatorics.

· Big Data Zone ·
Free Resource

The open source HPCC Systems platform is a proven, easy to use solution for managing data at scale. Visit our Easy Guide to learn more about this completely free platform, test drive some code in the online Playground, and get started today.

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 following series representations:

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 g1 = 1/ f1 (things start out easy!) and for n ≥ 2:

...where the Bs 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. 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

Managing data at scale doesn’t have to be hard. Find out how the completely free, open source HPCC Systems platform makes it easier to update, easier to program, easier to integrate data, and easier to manage clusters. Download and get started today.

Topics:
big data ,mathematics ,data science

Published at DZone with permission of

Opinions expressed by DZone contributors are their own.

{{ parent.title || parent.header.title}}

{{ parent.tldr }}

{{ parent.urlSource.name }}