Padé Approximants

Perturbation Methods

Padé approximants replace a truncated power series by a RATIONAL function — a ratio of polynomials — whose own Taylor expansion matches the original series up to the highest order used to construct it. The substitution looks innocuous but does something powerful: rationals can capture POLES and other singularities that polynomial truncations cannot, and the radius of convergence (or, for asymptotic series, the optimal accuracy) frequently extends far beyond what the partial sum can reach. Almost every modern resummation technique either is a Padé approximant or contains one as a subroutine.

The defining construction. Given coefficients $a_{0}, a_{1}, \dots, a_{m + n}$ of a series $f (g) = \sum_{k} a_{k} g^{k}$ , the $[m / n]$ Padé approximant is the unique rational function

[m / n] (g) = \frac{P _{m} ( g )}{Q _{n} ( g )} = \frac{p _{0} + p _{1} g + \dots + p _{m} g ^{m}}{1 + q _{1} g + \dots + q _{n} g ^{n}}

whose Taylor expansion at $g = 0$ agrees with $f (g)$ through order $g^{m + n}$ . Normalizing $q_{0} = 1$ fixes the overall scale. The number of free parameters is $m + n + 1$ , matching the number of constraints imposed by reproducing $m + n + 1$ series coefficients — which is exactly why the approximant is uniquely determined (when it exists).

Why rationals beat polynomials

A polynomial is an entire function of $g$ — it has no poles, no branch cuts, no singularities at finite $g$ . A rational function has poles at the zeros of its denominator, and can therefore REPRODUCE any pole structure of the function being approximated. If the true $f (g)$ has a pole at $g = g_{*}$ , the partial sums of its Taylor series cannot converge for $∣ g ∣ > ∣ g_{*} ∣$ — they diverge because polynomials cannot have poles. But a Padé approximant can place a zero in its denominator near $g_{*}$ and continue providing meaningful values well past $g_{*}$ .

Even when the original series is DIVERGENT FOR ALL $g \neq = 0$ (the standard case for asymptotic expansions in quantum mechanics — the Bender-Wu series of the quartic anharmonic oscillator is the canonical example), Padé approximants frequently provide a finite, accurate value. The mechanism: the rational form effectively performs an implicit analytic continuation through the formal complex structure encoded in the coefficients.

Construction

The condition $P_{m} / Q_{n} = f + O (g^{m + n + 1})$ rearranges to $P_{m} = Q_{n} \cdot f + O (g^{m + n + 1})$ , which gives $m + n + 1$ linear equations for the $m + n + 1$ unknowns ${p_{0}, \dots, p_{m}, q_{1}, \dots, q_{n}}$ . The structure of those equations splits into two halves: the denominator coefficients $q_{k}$ are determined first from the conditions that orders $g^{m + 1}, \dots, g^{m + n}$ vanish, and then the numerator coefficients fall out from convolution with the original series.

Concretely:

k = 0 \sum n q_{k} a_{m + i - k} = 0, i = 1, 2, \dots, n,

(with the convention $q_{0} = 1$ and $a_{j} = 0$ for $j < 0$ ). This is a linear system for ${q_{1}, \dots, q_{n}}$ with the Toeplitz-like matrix $a_{m + i - k}$ . Solving gives the denominator. The numerator is then

p_{i} = k = 0 \sum m i n (i, n) q_{k} a_{i - k}, i = 0, 1, \dots, m .

Two practical notes. (1) Padé approximants do not always exist: the linear system can be singular when the series has special structure (e.g., consecutive zero coefficients), in which case you fall back to a non-diagonal $[m / n]$ with adjusted indices. (2) DIAGONAL approximants $[m / m]$ are usually the right starting choice — they balance the polynomial degrees and produce the best convergence for most asymptotic problems.

Code

# Pade approximants: replace a truncated power series by a rational
# function of polynomials. The [m/n] Pade is the rational function
# P(g)/Q(g) with deg P = m, deg Q = n, whose Taylor expansion matches
# the series up to order m + n.
#
# Algorithm: solve a small linear system for the denominator coefficients
# (using the Pade equations from matching), then back-substitute for the
# numerator. Robust unless the linear system is singular — a sign that
# the requested approximant doesn't exist (rare for diagonal [m/m]).

import numpy as np
from scipy.linalg import eigh
from math import factorial

def pade(coeffs, m, n):
    """Build the [m/n] Pade approximant from coefficients a_0, ..., a_{m+n}.
    Returns numerator coefficients a (length m+1) and denominator
    coefficients b (length n+1, with b[0] = 1 by normalization)."""
    M = m + n
    c = np.asarray(coeffs[:M+1], dtype=float)
    if n > 0:
        # System for b[1], ..., b[n]:  sum_{k=0..n} b[k] c[m+i-k] = 0
        A = np.array([[c[m + i - k] if 0 <= m + i - k < len(c) else 0
                       for k in range(1, n + 1)]
                      for i in range(1, n + 1)])
        b_tail = np.linalg.solve(A, -c[m+1 : m+n+1])
        b = np.concatenate([[1.0], b_tail])
    else:
        b = np.array([1.0])
    # Numerator from direct convolution: a[i] = sum_k b[k] c[i-k]
    a = np.zeros(m + 1)
    for i in range(m + 1):
        a[i] = sum(b[k] * c[i - k] for k in range(min(i, n) + 1))
    return a, b

def pade_eval(a, b, g):
    num = sum(a[i] * g**i for i in range(len(a)))
    den = sum(b[i] * g**i for i in range(len(b)))
    return num / den

# ─── Demonstration: AHO ground state energy ──────────────────────────────
# H = -1/2 d²/dx² + 1/2 x² + g x⁴
# Compute the "true" energy by SHO-basis diagonalization, then compare
# against partial sums and Pade.

Nbasis = 80
X = np.zeros((Nbasis, Nbasis))
for n in range(Nbasis - 1):
    X[n, n+1] = np.sqrt((n + 1) / 2)
    X[n+1, n] = np.sqrt((n + 1) / 2)
X4 = X @ X @ X @ X
H0 = np.diag(np.arange(Nbasis) + 0.5)

def aho_exact(g):
    return float(eigh(H0 + g * X4, eigvals_only=True, subset_by_index=[0, 0])[0])

# Bender-Wu coefficients via Rayleigh-Schrodinger PT in the same basis
def bender_wu(Nmax):
    V = X4
    E = [0.5]
    psi = [np.zeros(Nbasis) for _ in range(Nmax + 1)]
    psi[0][0] = 1.0
    R = np.array([0.0] + [1.0 / n for n in range(1, Nbasis)])
    def Rapply(v):
        out = R * v; out[0] = 0.0; return out
    for k in range(1, Nmax + 1):
        Vpsi = V @ psi[k-1]
        Ek = psi[0] @ Vpsi
        E.append(Ek)
        rhs = -Vpsi + sum(E[j] * psi[k-j] for j in range(1, k))
        psi[k] = Rapply(rhs)
    return E

A = bender_wu(11)
print(f"Bender-Wu: a_0={A[0]:.4f}, a_1={A[1]:.4f}, a_2={A[2]:.4f}, a_3={A[3]:.4f}")
print(f"          (exact: 1/2, 3/4, -21/8 = -2.625, 333/8 = 41.625)")
print()

print(f"{'g':>6s}  {'true E_0':>14s}  {'N=11 partial':>15s}  {'[5/5] Pade':>13s}  {'Pade err':>10s}")
for g in [0.02, 0.05, 0.1, 0.2]:
    true = aho_exact(g)
    ps   = sum(A[n] * g**n for n in range(12))
    a_p, b_p = pade(A, 5, 5)
    pv   = pade_eval(a_p, b_p, g)
    print(f"{g:6.2f}  {true:14.10f}  {ps:15.4e}  {pv:13.10f}  {abs(pv - true):.2e}")

Output:

Bender-Wu: a_0=0.5000, a_1=0.7500, a_2=-2.6250, a_3=20.8125
          (exact: 1/2, 3/4, -21/8 = -2.625, 333/8 = 41.625)

     g        true E_0     N=11 partial     [5/5] Pade    Pade err
  0.02    0.5140864273      5.1409e-01    0.5140864273   2.18e-11
  0.05    0.5326427548      5.3505e-01    0.5326427277   2.71e-08
  0.10    0.5591463272      6.6147e+00    0.5591442026   2.12e-06
  0.20    0.6024051637      1.4020e+04    0.6023359714   6.92e-05

Three things to read off. (1) The first-principles Padé coefficients $a_{n}$ match the Bender-Wu formulas exactly. (2) At $g = 0.1$ , the order-11 PARTIAL SUM of the AHO series is at 6.6 — wildly diverging from the true value $\sim 0.559$ — while the $[5/5]$ Padé built from the SAME 11 coefficients gives 0.55914, off by only $2 \times 1 0^{- 6}$ . (3) At the still-divergent $g = 0.2$ , the partial sum has blown up to $1 0^{4}$ while Padé delivers 4-decimal accuracy. The series didn't converge to give the answer; Padé extracted the answer from the divergent series.

Reading singularities off the Padé denominator

A diagnostic worth keeping in your toolkit: the zeros of the Padé DENOMINATOR $Q_{n} (g)$ are estimates of the singularities of the true function. If your series comes from an integral or differential equation and you want to know "where does this function break down?", build a few diagonal Padé approximants and root-find on the denominators. Stable zeros (those that appear in roughly the same place across $[m / m]$ for several $m$ ) are real singularities of $f$ ; unstable zeros (different positions for different $m$ ) are spurious "Froissart doublets" — pole-zero pairs the algorithm produces in regions of analyticity. The technique is part of the standard analytic-continuation toolkit.

Where Padé fails

Padé is a workhorse but not a panacea. The classic failure modes:

Branch cuts. A rational function cannot represent a branch cut. If $f (g) \sim g - g_{*}$ near $g_{*}$ , Padé instead places a string of alternating poles and zeros along the cut — Froissart doublets again — and the approximant is poor in the neighbourhood. Conformal mapping (mapping $g$ to a variable in which the cut becomes a unit circle) is the standard workaround.
Wildly oscillating coefficients. When the $a_{n}$ grow factorially with alternating signs (the AHO case at large $g$ ), diagonal Padé converges only slowly. Borel-Padé is the standard upgrade: take a Borel transform first to tame the factorial growth, then Padé-approximate the much better-behaved Borel transform.
Strong coupling. For values of $g$ far outside the asymptotic regime of the original series, Padé struggles and the accuracy plateaus. Variational perturbation theory is the modern alternative.

Borel-Padé resummation — the natural upgrade for factorially divergent asymptotic series.
Perturbation theory and resurgence — Padé approximants are the underlying technique behind much of the practical work on resurgent transseries.
Shanks transformation and Wynn's ε-algorithm — for convergent series, a different acceleration philosophy.
Bi-orthogonal Lanczos — continued-fraction representations of resolvents are mathematically equivalent to Padé approximants of the appropriate generating function. Padé and Lanczos are two faces of the same idea.