“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Padé via Lanczos (Feldmann & Freund 1995): faithful reproduction

Circuits

Padé via Lanczos (PVL) is the algorithm Feldmann and Freund published in 1995 to fix the numerical disaster that was Asymptotic Waveform Evaluation (AWE). AWE had been the workhorse for fast linear-circuit transfer-function approximation since Pillage and Rohrer 1990, but everyone in the room knew it broke down past about $q = 10$ poles — the explicit moment matrix became so ill-conditioned that any further Padé coefficients were noise. F&F's contribution was to realize that the same Padé approximant could be assembled without ever forming the moment matrix, by running a Lanczos iteration on the same operator and reading the approximant off directly from the tridiagonalization. This page reproduces the paper's two key diagnostic experiments — AWE moment-matrix condition numbers blowing up, PVL convergence to machine precision — on a small symmetric RC pencil, and walks through the F&F section 4 derivation that makes the Lanczos↔Padé identity precise.

This is the paper-faithful companion to /circuits/lanczos_mor, which builds the same Krylov-subspace machinery as a continued fraction in the symmetric case. The two pages use the same RC ladder; the difference is in framing. The MOR page emphasizes the physics (where do the Ritz values land in the spectrum, why do short chains miss the low-frequency mode); this page emphasizes the algorithm (why AWE failed, how PVL routes around the failure, and what the underlying algebra of the Padé identity is).

Setup: MNA descriptor form and the transfer function

Following F&F section 2: a lumped linear time-invariant circuit obeys the descriptor-form ODE

C \dot{x} = - G x + b u, y = l^{T} x .

$G$ collects memoryless elements (resistors, algebraic constraints, sources), $C$ collects memory elements (capacitors, inductors), $b$ is the input selector, and $l$ is the output projection. Laplace-transform and rearrange to get the transfer function as a matrix element of the resolvent of the pencil $(G, C)$ :

H (s) = l^{T} (G + s C)^{- 1} b .

For an interconnect mesh, the pencil is huge — millions of nodes — and you need $H (s)$ across many decades of frequency. Re-factoring $G + s C$ at every $s$ is out of the question; you want a reduced rational approximation valid across the operating band.

Change of variables: shift to an expansion point

Pick $s_{0} \in C$ in the operating band of interest (real, typically). Provided $G + s_{0} C$ is invertible, change variables to $s = s_{0} + σ$ and define

A := - (G + s_{0} C)^{- 1} C, r := (G + s_{0} C)^{- 1} b .

Algebra gives

H (s_{0} + σ) = l^{T} (I - σ A)^{- 1} r .

Taylor-expand the resolvent around $σ = 0$ :

H (s_{0} + σ) = k = 0 \sum \infty m_{k} σ^{k}, m_{k} := l^{T} A^{k} r .

The moments $m_{k}$ are what AWE tries to compute explicitly. They're available via the recursive solve $(G + s_{0} C) u_{k} = C u_{k - 1}$ with $u_{0} = r$ , then $m_{k} = l^{T} u_{k}$ . One LU factorization, then $2 q$ triangular solves to get $m_{0}$ through $m_{2 q - 1}$ . So far so good.

The Padé approximant — what it is, what AWE tries to do

The $q -th$ Padé approximant of $H (s_{0} + σ)$ in F&F's convention is the rational function

H_{q} (s_{0} + σ) = \frac{b _{q - 1} σ ^{q - 1} + \dots + b _{1} σ + b _{0}}{a _{q} σ ^{q} + \dots + a _{1} σ + 1}

whose Taylor expansion matches $H (s_{0} + σ)$ in the first $2 q$ coefficients. That's the $(q - 1, q)$ Padé entry, and it's the natural choice for a transfer function that decays to zero at high frequency.

AWE solves for the $a_{j}$ coefficients via the moment Hankel system

M_{q} a_{q} ⋮ a_{1} = - m_{q} ⋮ m_{2 q - 1}, M_{q} = [m_{j + k - 2}]_{j, k = 1, \dots, q}

then roots the denominator polynomial to get the $q$ poles, and reads off residues via another linear solve.

Why AWE fails: the moments collapse onto one eigenvalue

The vectors $u_{k} = A^{k} r$ AWE generates are exactly the right Krylov sequence for $A$ starting from $r$ . Power iteration: they converge to an eigenvector of $A$ corresponding to the eigenvalue of $A$ with largest magnitude. So the moments $m_{k} = l^{T} u_{k}$ rapidly collapse to information about that one eigenvalue — and you can't extract $q$ distinct poles from a sequence that only knows about one.

The symptom is in the condition number of the moment matrix $M_{q}$ . Build a small symmetric RC ladder (200 nodes, $R = 1 Ω$ , $C = 1 nF$ ), expand around $s_{0} = 1 0^{6}$ , and compute the AWE condition numbers up to $q = 16$ :

import numpy as np
import scipy.linalg as sla
import scipy.sparse as sp
import scipy.sparse.linalg as spla

# Symmetric RC ladder: the simplest pencil where the AWE-vs-PVL story is
# visible. F&F's algorithm extends to non-symmetric (G, C) via biorthogonal
# Lanczos with separate left and right starting vectors — the math
# specializes here because the pencil is symmetric.
N = 200
R = 1.0; C = 1e-9
g = 1.0 / R
diag = np.full(N, 2.0 * g); diag[0] = g; diag[-1] = g; diag += 1e-9
off = np.full(N - 1, -g)
G  = sp.diags([off, diag, off], [-1, 0, 1], format='csr')
Cm = sp.diags([C] * N, format='csr')
b  = np.zeros(N); b[0] = 1.0     # current source at node 0
l  = np.zeros(N); l[0] = 1.0     # observe voltage at node 0 (driving-point Z)

# AWE: compute moments  m_k = l^T A^k r  where  A = -(G + s0 C)^{-1} C,
# then form the Hankel moment matrix M_q[j,k] = m_{j+k} and try to solve
# for the Pade denominator coefficients.

def awe_moments(G, Cm, b, l, s0, K):
    M_pencil = (G + s0 * Cm).tocsc()
    LU = spla.splu(M_pencil)
    r = LU.solve(b)
    u = r.copy()
    m = [float(l @ u)]
    for _ in range(1, K):
        u = -LU.solve(Cm @ u)            # u_k = A^k r
        m.append(float(l @ u))
    return np.array(m)

s0 = 1e6                                  # expansion point, in the operating band
moments = awe_moments(G, Cm, b, l, s0, 2 * 16)
print("AWE moment-matrix condition numbers")
print(f"{'q':>4}  {'cond(M_q)':>16}")
for q in [2, 4, 6, 8, 10, 12, 14, 16]:
    M_q = np.array([[moments[j + k] for k in range(q)] for j in range(q)])
    print(f"{q:>4}  {np.linalg.cond(M_q):>16.3e}")

AWE moment-matrix condition numbers
   q         cond(M_q)
   2         8.123e+12
   4         1.681e+39
   6         1.267e+66
   8         3.492e+94
  10        1.842e+123
  12        1.478e+146
  14        2.692e+171
  16        3.058e+195

Same shape as F&F Table 1: $κ (M_{q}) \sim 1 0^{12}$ at $q = 2$ , climbing past $1 0^{195}$ by $q = 16$ . Each factor-of-10 increase in $κ (M_{q})$ is one more decimal digit of accuracy lost when solving the Hankel system. By $q = 6$ the $M_{q}$ system is solvable only in the formal sense — the numerical answer is dominated by round-off.

Verify the consequence: try to extract the poles by solving the Hankel system and rooting the denominator.

# AWE attempts to extract poles by solving the Hankel system for the
# denominator polynomial and rooting it.

print(f"{'q':>4}  {'cond(M_q)':>14}  {'pole extraction result':>40}")
for q in [4, 8, 12, 16]:
    M_q = np.array([[moments[j + k] for k in range(q)] for j in range(q)])
    rhs = -moments[q:2 * q]
    try:
        denom_coeffs = np.linalg.solve(M_q, rhs)
        poly = np.concatenate([denom_coeffs[::-1], [1.0]])
        poles = np.roots(poly)
        n_real = sum(1 for p in poles if abs(p.imag) < 1e-10 * abs(p.real))
        msg = f"{n_real}/{q} poles real; max|pole| = {np.max(np.abs(poles)):.2e}"
    except np.linalg.LinAlgError:
        msg = "SINGULAR Hankel system"
    print(f"{q:>4}  {np.linalg.cond(M_q):>14.2e}  {msg:>40}")

   q       cond(M_q)                  pole extraction result
   4        1.68e+39      0/4 poles real; max|pole| = 2.63e+01
   8        3.49e+94      0/8 poles real; max|pole| = 4.91e+00
  12       1.48e+146     0/12 poles real; max|pole| = 3.22e+00
  16       3.06e+195     0/16 poles real; max|pole| = 2.20e+00

For this pencil all true poles of $H (s)$ are real (the poles are at $s = - λ_{j}$ where $λ_{j}$ are eigenvalues of the symmetrized Cholesky-conjugated $K = L^{- 1} G L^{- T}$ , which is symmetric positive-definite). AWE returns zero real poles out of $q$ at every order tested. The Hankel-system solve is so corrupted that the resulting denominator polynomial doesn't have a single real root in the right ballpark. The numerical Padé approximant is garbage.

The PVL idea: route the same Padé approximant through Lanczos

F&F's observation: there's another way to compute the same Padé approximant that doesn't form $M_{q}$ at all. Run a Lanczos iteration on $A$ . In F&F's general non-Hermitian setting that means biorthogonal Lanczos (Algorithm 1 in the paper), with right starting vector $r$ and left starting vector $l$ :

A V_{q} = V_{q} T_{q} + ρ_{q + 1} v_{q + 1} e_{q}^{T}, A^{T} W_{q} = W_{q} \tilde{T}_{q} + η_{q + 1} w_{q + 1} e_{q}^{T} .

With biorthogonality $W_{q}^{T} V_{q} = D_{q}$ (diagonal). Then F&F show — and the derivation walk below reproduces step by step — that the $q -th$ Padé approximant of $H (s_{0} + σ)$ equals

H_{q} (s_{0} + σ) = (l^{T} r) e_{1}^{T} (I - σ T_{q})^{- 1} e_{1} .

No moment matrix. No Hankel system. The Padé coefficients are implicit in the Lanczos tridiagonal $T_{q}$ . Diagonalize $T_{q} = S_{q} Λ_{q} S_{q}^{- 1}$ and you have the poles and residues directly: $p_{j} = s_{0} + 1/ λ_{j}$ , $k_{j} = - (l^{T} r) μ_{j} ν_{j} / λ_{j}$ with $μ = S_{q}^{T} e_{1}$ and $ν = S_{q}^{- 1} e_{1}$ . F&F's Algorithm 2.

Derivation walk: from the Lanczos relations to the Padé identity

This is F&F section 4 spelled out in algebraic moves. Each step names the move so you can verify mentally instead of pattern-matching. Toggle "Every step" to see the micro algebra; use arrow keys to step.

What you'll see: the right and left Lanczos relations get pushed against the moment

m_{k} = l^{T} A^{k} r

until biorthogonality collapses the inner factors and the moment becomes

(l^{T} r) \cdot e_{1}^{T} T_{q}^{k} e_{1}

. Summing the Taylor series of

H (s_{0} + σ)

with this substitution gives the closed-form PVL formula, and an eigendecomp of

T_{q}

splits it into pole-residue form.

Step 1 of 15 · ← → keys to navigate

Step 1

A V_{q} = V_{q} T_{q} + ρ_{q + 1} v_{q + 1} e_{q}^{T}

The right Lanczos relation from Feldmann-Freund Algorithm 1.

V_{q} = [v_{1} v_{2} \dots v_{q}]

is the matrix of right Lanczos vectors built from starting vector

v_{1} = r /∥ r ∥

. The tridiagonal

T_{q}

carries the Lanczos

α

's on the diagonal and

β

's/

ρ

's on the off-diagonals. The trailing column

ρ_{q + 1} v_{q + 1} e_{q}^{T}

is the residual that hasn't been absorbed yet.

PVL demo: convergence to machine precision

The same RC ladder, same expansion point. PVL specialized to the symmetric pencil collapses to a symmetric Lanczos run (the biorthogonal vectors become equal up to scaling) on the Cholesky-symmetrized $K = L^{- 1} G L^{- T}$ . The full non-symmetric F&F algorithm — needed for general circuits with inductors that break symmetry — uses biorthogonal Lanczos with two starting vectors. The Lanczos↔Padé identity derived above doesn't care about symmetry; the simplification is purely algorithmic.

# PVL on the symmetric pencil (F&F section 4 specialized to G = G^T, C = C^T).
# Symmetrize:  C = L L^T  (diagonal, so L = sqrt(C))
#              K = L^{-1} G L^{-T}  (sparse, symmetric, positive-definite)
#              (G + s C) = L (K + s I) L^T
# Run symmetric Lanczos on  A_sym = -(K + s0 I)^{-1}  starting from
# y_0 = (K + s0 I)^{-1} bp,  bp = L^{-1} b.

def lanczos_symmetric(matvec, v0, n_iter):
    alphas, betas, Q = [], [], [v0.copy()]
    q_prev = np.zeros_like(v0); q = v0.copy(); beta_prev = 0.0
    for j in range(n_iter):
        w = matvec(q)
        alpha = float(q @ w)
        w = w - alpha * q - beta_prev * q_prev
        for qk in Q[:-1]:                    # full reorthogonalization
            w -= float(qk @ w) * qk
        beta = float(np.linalg.norm(w))
        alphas.append(alpha)
        if beta < 1e-14:
            break
        if j < n_iter - 1:
            betas.append(beta)
            q_prev = q; q = w / beta
            Q.append(q.copy())
    return np.array(alphas), np.array(betas), np.column_stack(Q)

def pvl(s_grid, G, Cm, b, l, s0, q):
    sqrtC = np.sqrt(Cm.diagonal()); Linv = 1.0 / sqrtC
    K = (sp.diags(Linv) @ G @ sp.diags(Linv)).tocsc()
    K_plus = (K + s0 * sp.diags(np.ones(K.shape[0]))).tocsc()
    LU = spla.splu(K_plus)
    bp = Linv * b; lp = Linv * l
    y0 = LU.solve(bp); norm_y0 = np.linalg.norm(y0)
    v0 = y0 / norm_y0
    alphas, betas, V = lanczos_symmetric(lambda x: -LU.solve(x), v0, q)
    T = np.diag(alphas) + np.diag(betas, 1) + np.diag(betas, -1)
    coeff_row = norm_y0 * (lp @ V)
    eigvals, S = sla.eigh(T)
    e1 = np.zeros(T.shape[0]); e1[0] = 1.0
    cS = coeff_row @ S
    Se1 = S.T @ e1
    H_q = np.zeros(len(s_grid), dtype=complex)
    for i, s in enumerate(s_grid):
        sigma = s - s0
        H_q[i] = np.sum(cS * Se1 / (1.0 - sigma * eigvals))
    return H_q

omegas = np.logspace(3, 9, 200)
H_ref = np.array([complex(l @ spla.spsolve((G + 1j * w * Cm).tocsc(), b))
                   for w in omegas])
print("PVL transfer-function convergence")
print(f"{'q':>4}  {'max rel error':>16}")
for q in [4, 8, 16, 32, 64, 128]:
    Hq = pvl(1j * omegas, G, Cm, b, l, s0, q)
    rel_err = np.max(np.abs(Hq - H_ref)) / np.max(np.abs(H_ref))
    print(f"{q:>4}  {rel_err:>16.3e}")

PVL transfer-function convergence
   q     max rel error
   4         9.798e-01
   8         1.560e-04
  16         1.297e-05
  32         4.115e-10
  64         1.623e-10
 128         1.625e-10

Compare with the AWE table from above. AWE's moment matrix at $q = 16$ has condition number $3 \times 1 0^{195}$ and produces zero real poles out of 16. PVL at $q = 32$ matches the direct sparse solve to ten significant digits. Same problem, same pencil, same expansion point — different algorithm.

Two things to notice about the convergence. First, the dip at $q = 4$ is wild (almost 100% error) because the Krylov subspace hasn't reached the low-eigenvalue mode that dominates low-frequency impedance (this is the failure mode the MOR page diagnoses in detail). Second, once $q$ is large enough to capture the relevant spectrum, convergence is essentially limited by round-off ( $\sim 1 0^{- 10}$ for double precision plus a moderately conditioned problem). That's the PVL story: monotone, predictable, controlled by problem geometry rather than by numerical degradation of the algorithm.

Pole-quality bound

PVL produces both the approximant and a quality measure on the poles. The poles are $1/ λ_{j}$ in the shifted variable, with $λ_{j}$ eigenvalues of $T_{q}$ . These are Ritz values of $A$ in the Krylov basis. Each Ritz pair $(λ_{j}, z_{j})$ with $z_{j} := V_{q} s_{j}$ has an explicit residual

A z_{j} - λ_{j} z_{j} = ρ_{q + 1} v_{q + 1} s_{q j}

where $s_{q j}$ is the last component of $s_{j}$ . Taking norms and rescaling by $n (A) \approx ∥ A ∥_{2}$ (easily estimated from the Lanczos coefficients themselves) gives the F&F pole-quality measure

Q_{j} := \frac{ρ _{q + 1} ∣ s _{q j} ∣}{n ( A ) \cdot ∥ s _{j} ∥ _{2}} \geq \frac{∥ A z _{j} - λ _{j} z _{j} ∥ _{2}}{∥ A ∥ _{2} \cdot ∥ s _{j} ∥ _{2}} .

Small $Q_{j}$ = converged pole; large $Q_{j}$ = spurious pole that's an artifact of the projection. This is the other thing AWE can't give you. AWE produces $q$ poles and you have no idea which ones to trust; PVL produces $q$ poles AND $q$ quality numbers, and you keep the ones with small $Q_{j}$ .

What PVL doesn't fix, and what came after

Two limitations of single-shift PVL, both addressed by follow-up work:

Wide-band coverage. PVL is local to $s_{0}$ . Accuracy degrades far from the expansion point. F&F cite complex frequency hopping (CFH: Chiprout-Heeb-Nakhla-Ruehli 1993) for multi-point Padé panels stitched across decades of frequency. CFH costs several matrix factorizations, one per shift; PVL with a single well-chosen $s_{0}$ often matches CFH's accuracy with one factorization (their PEEC example, 60 PVL iterations vs. CFH up to 5GHz).
Passivity preservation. The PVL reduced model doesn't generically remain passive (it can violate energy dissipation rules of the original network). Odabasioglu-Celik-Pileggi 1998 fix this with PRIMA: block Krylov projection followed by careful symmetrization that preserves the original pencil's positive-real property. Same Lanczos↔Padé math underneath; the addition is a structure-preserving twist on the projection.

There's also the issue of Lanczos breakdowns. The biorthogonal recurrence can hit $δ_{n} = w_{n}^{T} v_{n} = 0$ mid-chain even when neither vector is zero (Gutknecht's "serious breakdown," corresponding to singular blocks in the Padé table). F&F's implementation of PVL uses the look-ahead Lanczos variant of Freund-Gutknecht-Nachtigal 1993 that skips over these blocks. For symmetric pencils (the demo above) breakdowns don't occur, which is why the demo can use straight symmetric Lanczos.

Related on this site

Lanczos + Continued Fractions for Circuit Model-Order Reduction — same RC ladder, different framing: emphasizes where Ritz values land in the spectrum and why short chains miss low-frequency modes. Use that page first for the physics, then this one for the F&F algorithm and Padé identity.
Lanczos iteration — the underlying algorithm, with the explicit tridiagonalization derivation.
Bi-orthogonal Lanczos — the non-Hermitian variant that PVL actually uses for general circuits (the demo above specializes to symmetric because biorthogonal Lanczos collapses there).
Liouville-Lanczos — the physics analog: same Krylov+continued-fraction machinery applied to TDDFT linear response in solids.

Reference

P. Feldmann and R.W. Freund, "Efficient Linear Circuit Analysis by Padé Approximation via the Lanczos Process," 1995. The numerical experiments above reproduce Table 1 (moment-matrix condition numbers) and Figure 2 (PVL convergence) of that paper on a small symmetric RC ladder. The derivation walk reproduces the algebra of F&F section 4 (equations 19–27) in algebraic-move-named steps.