Inverting the Pade approximation (going from $P_m(x)/Q_n(x)$ to $f_{m+n}(x)$)

Question

Is there a way to easily invert the Pade approximation and get back the Taylor series it was derived from? In other words, if I have two polynomials $P_m(x)/Q_n(x)$ derived from an $(m+n)$th degree polynomial $f_{m+n}(x)$, can I get that polynomial back easily?

(specifically I'm looking for the inverse of the scipy function scipy.misc.pade)

I found the Wikipedia link on power series of a rational function, but I'm having trouble figuring out how to convert it to a matrix equation to solve numerically... any help?

Answer 1

oh wait... the fog is beginning to clear here... if $\mathbf P$ and $\mathbf Q$ are the coefficients in column vector form e.g. $\mathbf P = \begin{bmatrix}p_0\\p_1\\ \vdots \\ p_m\end{bmatrix}$, then we have

$$\sum\limits_{k=0}^{m} p_kx^k = \left(\sum\limits_{k=0}^{m+n} a_kx^k\right)\left(\sum\limits_{k=0}^{n} q_kx^k\right)$$

which translates into matrix form as

$$ \begin{bmatrix}{\mathbf P} \\{\mathbf 0} \end{bmatrix} = {\mathbf K \mathbf A} = \begin{bmatrix} q_0 \\ q_1 & q_0 \\ q_2 & q_1 & q_0\\ q_3 & q_2 & q_1 & q_0\\ \vdots & \vdots & \vdots & \vdots & \ddots \\ q_n & q_{n-1} & q_{n-2} & q_{n-3} & \cdots & q_0 & \\ & q_{n} & q_{n-1} & q_{n-2} & \cdots & q_1 & q_0\\ && q_{n} & q_{n-1} & \cdots & q_2 & q_1 & q_0 \end{bmatrix} {\mathbf A}$$

where $\mathbf K$ is a Toeplitz matrix with $m+n+1$ rows and columns (looks like I can use scipy.linalg.toeplitz ), and then I can solve for $\mathbf A$.

Not sure I'm going to get back the same polynomial as the original power series, though.

a trial in Python:

import scipy.linalg
import numpy as np

def epoly(n):
    ''' returns the degree n polynomial Taylor series for e^-x '''
    def helper():
        a = 1.0
        yield a
        for i in xrange(1,n+1):
            a = -a / i
            yield a
    return np.poly1d(list(helper())[::-1])

def poly2array(p):
    if isinstance(p, np.poly1d):
        p = p.coeffs[::-1]
    else:
        p = np.array(p)
    assert len(p) > 1, "p must be a 1-D array"    
    return p

def rat2power(p,q,verbose=False):
    p = poly2array(p)
    q = poly2array(q)
    m = p.shape[0] - 1
    n = q.shape[0] - 1
    QQ = scipy.linalg.toeplitz(np.hstack((q,np.zeros(m))),np.hstack((1,np.zeros(m+n))))
    pv = np.array([np.hstack((p, np.zeros(n)))]).transpose()
    if verbose:
        print 'QQ=',QQ
        print 'pv=',pv
    return np.linalg.solve(QQ,pv).transpose()[0]

####

epoly5 = epoly(5)
epoly5coeffs = epoly5.coeffs[::-1]
print epoly5coeffs
p1,q1=scipy.misc.pade(epoly5coeffs, 3)

print 'p1 = ', p1
print 'q1 = ', q1
a=rat2power(p1,q1,verbose=True)
print a

which outputs

[ 1.         -1.          0.5        -0.16666667  0.04166667 -0.00833333]
p1 =        2
0.05 x - 0.4 x + 1
q1 =           3        2
0.01667 x + 0.15 x + 0.6 x + 1
QQ= [[ 1.          0.          0.          0.          0.          0.        ]
 [ 0.6         1.          0.          0.          0.          0.        ]
 [ 0.15        0.6         1.          0.          0.          0.        ]
 [ 0.01666667  0.15        0.6         1.          0.          0.        ]
 [ 0.          0.01666667  0.15        0.6         1.          0.        ]
 [ 0.          0.          0.01666667  0.15        0.6         1.        ]]
pv= [[ 1.  ]
 [-0.4 ]
 [ 0.05]
 [ 0.  ]
 [ 0.  ]
 [ 0.  ]]
[ 1.         -1.          0.5        -0.16666667  0.04166667 -0.00833333]

and this appears to work.