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]))
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
.