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

Hortonworks Sandbox for HDP and HDF is your chance to get started on learning, developing, testing and trying out new features. Each download comes preconfigured with interactive tutorials, sample data and developments from the Apache community.

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

Hortonworks Community Connection (HCC) is an online collaboration destination for developers, DevOps, customers and partners to get answers to questions, collaborate on technical articles and share code examples from GitHub.  Join the discussion.

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 }}