[DOC] Levinson–Durbin recursion for computing Pade approximants

[arvoelke]

This may be useful for analysis. A work-in-progress:

import numpy as np
from scipy.misc import factorial, pade

# Taylor series
c = 1. / factorial(np.arange(10))

# Array representing the Toeplitz matrix
a = c[-2::-1]
assert len(a) % 2 == 1  # odd

# Form the Toeplitz matrix explicitly
N = len(a) // 2 + 1  # even
T = np.empty((N, N))
for i in range(N):
    T[i, :] = a[N-i-1:2*N-i-1]

# Solution vector
y = -c[N:]
assert len(y) == N

# Levinson algorithm using the underlying a vector
f = np.zeros(N)
b = np.zeros(N)

f[0] = b[0] = 1. / a[N-1]
x = y[0] * b

for n in range(1, N):
    ef = a[N-n-1:N-1].dot(f[:n])
    eb = a[N:N+n].dot(b[:n])
    ex = a[N-n-1:N-1].dot(x[:n])

    fn = np.append(f[:n], [0])
    bn = np.append([0], b[:n])

    f[:n+1] = 1. / (1 - eb*ef) * fn - ef / (1 - eb*ef) * bn
    b[:n+1] = -eb / (1 - eb*ef) * fn + 1. / (1 - eb*ef) * bn
    x[:n+1] += (y[n] - ex) * b[:n+1]

assert np.allclose(np.dot(T, x), y)

# Relationship to Pade approximants
assert np.allclose(pade(c, N)[1], np.append(x[::-1], [1]))

References:

• A new algorithm for computing Pade approximants (2011)
• https://en.wikipedia.org/wiki/Levinson_recursion

Also note that computing d = D^{m+1}[q(x)*f(x) - p(x)] eliminates p and reveals the higher-order error O(x^n) in terms of q and f:

from scipy.special import binom

m = 2
n = 3

# Exponential taylor series
a = 1. / factorial(np.arange(m + n + 1))

assert m + n + 1 == len(a)
p, q = pade(a, n)
assert len(p) == m
assert len(q) == n

a = np.append(a, [0]*(m+1))  # padding for convenience
def x_pow(n): return np.poly1d([1., 0]) ** n

d = np.poly1d([0.])
for l in range(m+n+1):
    for i in range(l, l+m+2):
        d += a[i] * binom(m+1, i-l) * factorial(i) / factorial(l) * np.polymul(
            x_pow(l), np.polyder(q, m+1-(i-l)))
print d

and the derivatives can be grouped together into a single linear differential operator for each index l.

Derived using:

• https://en.wikipedia.org/wiki/General_Leibniz_rule
• https://www-sop.inria.fr/apics/anap03/PadeTalk.pdf (generalizing the first part of page 14)

[arvoelke]

Related open questions:

• http://math.stackexchange.com/questions/2158451/closed-form-solution-to-a-toeplitz-system
• http://math.stackexchange.com/questions/2158233/proof-of-closed-form-pad%C3%A9-approximants-for-ex
• Recurrence relation from lower-orders to higher-orders?
• Relationship to special case of m=1? (which is just a discrete convolution; I think because the Toeplitz matrix becomes lower triangular which makes Gaussian elimination easy)

[arvoelke]

See appendix of The Holomorphic Embedding Loadflow Method for DC Power Systems and Nonlinear DC Circuits (2015) for a review of methods for computing Pade approximants. Among these, the text Practical Extrapolation Methods (2003) is a wealth of relevant information (see section 17.3 pages 327-329, for example).