Approximation: Padé & Rational Interpolation

Overview

Rational approximation approximates a target function by a ratio of polynomials $R(x) = P(x)/Q(x)$. It shines exactly where polynomial approximation (Taylor expansion, Chebyshev approximation) struggles:

• Representing poles: as the denominator $Q(x)$ approaches zero, the approximation can reproduce divergences and singularities of the function. A polynomial only ever takes finite values, so it cannot represent a pole.
• Asymptotic behavior: a function that saturates to a constant as $x \to \infty$ (such as $\tanh$) is naturally approximated as the ratio of the leading terms of numerator and denominator.
• Reaching beyond the radius of convergence: even at points where the original Taylor series diverges, a diagonal Padé approximant often returns a meaningful value (an analytic-continuation-like behavior).

Given the same number of degrees of freedom (number of coefficients), rational approximation often has smaller error over a wider range than polynomial approximation. The price is that the zeros of the approximate denominator can create spurious poles (and, paired with a numerator zero, spurious zeros = Froissart doublets), so it takes discipline to inspect where the poles land.

sangi's Approximation module offers three families of rational approximation depending on the kind of input:

• From the Taylor coefficients at a single point → padeApprox / padeTable (Padé approximation)
• From Taylor coefficients to a continued-fraction representation → taylorToSFraction / taylorToJFraction (qd algorithm)
• From values at distinct sample points → rationalInterpolation (multipoint Padé / Cauchy interpolation)
• Matrix exponential $e^A$ → matrixExpPade (scaling & squaring)

Related API: Approximation. For polynomial and orthogonal-polynomial approximation, see the previous page Least Squares, Chebyshev & Minimax.

Padé Approximation

Given the Taylor coefficients $a_0, a_1, \ldots$ of a function $f$ around $x = 0$, the $[M/N]$ Padé approximant is the following rational function:

$$R_{M/N}(x) = \frac{P_M(x)}{Q_N(x)} = \frac{p_0 + p_1 x + \cdots + p_M x^M}{1 + q_1 x + \cdots + q_N x^N}$$

The denominator is normalized as $Q_N(0) = 1$. The coefficients are determined by requiring that the Taylor expansion of $R_{M/N}$ match the Taylor series of $f$ up to degree $M+N$:

$$f(x) - \frac{P_M(x)}{Q_N(x)} = O(x^{M+N+1})$$

Coefficient-matching conditions and the linear system

Multiplying both sides by $Q_N(x)$ and comparing coefficients at each power of $x$ yields the following conditions (here $a_j = 0\ (j < 0)$):

$$\sum_{l=0}^{N} q_l\, a_{m-l} = p_m \quad (0 \le m \le M), \qquad \sum_{l=0}^{N} q_l\, a_{m-l} = 0 \quad (M < m \le M+N)$$

Since $q_0 = 1$, the latter $N$ equations form a linear system in the denominator coefficients $q_1, \ldots, q_N$:

$$\begin{pmatrix} a_M & a_{M-1} & \cdots & a_{M-N+1} \\ a_{M+1} & a_M & \cdots & a_{M-N+2} \\ \vdots & & \ddots & \vdots \\ a_{M+N-1} & a_{M+N-2} & \cdots & a_M \end{pmatrix} \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_N \end{pmatrix} = -\begin{pmatrix} a_{M+1} \\ a_{M+2} \\ \vdots \\ a_{M+N} \end{pmatrix}$$

Once the denominator is solved, the numerator follows explicitly from the first set of relations as a convolution:

$$p_m = a_m + \sum_{l=1}^{\min(m,N)} q_l\, a_{m-l} \quad (0 \le m \le M)$$

sangi's padeApprox(taylor, M, N) follows exactly this procedure. It solves the $N \times N$ Toeplitz-structured system by LU decomposition (algorithms::solve), constructs the numerator by convolution, and returns a PadeResult<T> (numerator / denominator / valid). When the denominator system is singular (a degenerate Padé) and cannot be solved, it returns valid = false rather than silently returning a wrong approximation. For evaluation, evaluatePade computes the numerator and denominator each by Horner's method and takes their ratio.

The Padé table

Arranging $[M/N]$ into a two-dimensional array with row $M$ and column $N$ gives the Padé table:

The 0th column ($N=0$) is the Taylor polynomial itself, and the diagonal $[M/M]$ and the antidiagonal $M+N=\text{const}$ are each important as convergence sequences. In practice the diagonal Padé approximants are often the most stable and the most accurate. sangi's padeTable(taylor, Mmax, Nmax) generates every cell with $0 \le M \le M_{\max},\ 0 \le N \le N_{\max}$ at once and returns it as a vector<vector<PadeResult<T>>>. Degenerate cells are included with valid = false. Since each cell solves an $N \times N$ system, the cost is $O(M_{\max} \cdot N_{\max} \cdot N^2)$, so if only the diagonal is needed it is cheaper to call padeApprox directly.

Approximation beyond the radius of convergence

Even for functions with a finite radius of convergence such as $f(x) = \log(1+x)$ or $f(x) = \tan x$, a Padé approximant can return a meaningful value outside the circle of convergence. This is because the zeros of the denominator approach the locations of the true poles of $f$, so the rational function plays the role of analytic continuation. Note, however, that unrelated nearby zero-pole pairs (Froissart doublets) can appear, so for high-order Padé approximants you should always inspect the zeros of the denominator.

Continued Fractions and the qd Algorithm

The continued-fraction representation is two sides of the same coin as Padé approximation. Converting a Taylor series into a continued fraction makes evaluation robust against cancellation, and truncating partway automatically yields a sequence of Padé approximants.

Stieltjes continued fraction (S-fraction)

An S-fraction has a first-degree term at each level:

$$f(z) = \cfrac{c_0}{1 + \cfrac{a_1 z}{1 + \cfrac{a_2 z}{1 + \cfrac{a_3 z}{1 + \cdots}}}}$$

Its $k$-th approximant (truncated partway) corresponds to a zigzag sequence in the Padé table. In sangi, StieltjesFraction<T> holds the leading constant b0 $= c_0$, the coefficient list a $= [a_1, a_2, \ldots]$, and terminated indicating whether it closed off at finite length.

Jacobi continued fraction (J-fraction)

A J-fraction groups the S-fraction two levels at a time, giving second-degree terms:

$$f(z) = \cfrac{c_0}{1 - \beta_0 z - \cfrac{\alpha_1^2 z^2}{1 - \beta_1 z - \cfrac{\alpha_2^2 z^2}{1 - \beta_2 z - \cdots}}}$$

This corresponds to the three-term recurrence of orthogonal polynomials, and the J-fraction up to level $n$ coincides with the diagonal $[n/n]$ Padé approximant. sangi's JacobiFraction<T> holds c0, the coefficient list beta ($\beta_0, \beta_1, \ldots$), and alpha2 ($\alpha_1^2, \alpha_2^2, \ldots$).

Taylor series → continued-fraction conversion (qd algorithm)

The core of the conversion is Rutishauser's qd (quotient-difference) algorithm. It generates the $q$-column and the $e$-column alternately by a recurrence from the Taylor coefficients. The first row begins with

$$q_1^{(n)} = \frac{c_{n+1}}{c_n}, \qquad e_0^{(n)} = 0$$

and the table is filled in by the rhombus rule:

$$e_{k-1}^{(n)} = q_{k-1}^{(n+1)} - q_{k-1}^{(n)} + e_{k-2}^{(n+1)}, \qquad q_k^{(n)} = q_{k-1}^{(n+1)}\, \frac{e_{k-1}^{(n+1)}}{e_{k-1}^{(n)}}$$

The continued-fraction coefficients are determined from the top-row values. The S-fraction comes from flipping the signs:

$$a_{2k-1} = -q_k^{(0)}, \qquad a_{2k} = -e_k^{(0)},$$

and the J-fraction comes from grouping two levels:

$$\beta_0 = q_1^{(0)}, \quad \beta_n = q_{n+1}^{(0)} + e_n^{(0)}\ (n \ge 1), \quad \alpha_n^2 = q_n^{(0)}\, e_n^{(0)}$$

sangi's taylorToSFraction(taylor, kmax) / taylorToJFraction(taylor, kmax) build this qd table internally. If $q_k$ or $e_k$ becomes $0$ partway (a degenerate point of the Padé), the process truncates there and sets terminated = true. When the leading coefficient $c_0 = 0$, $q_1 = c_1/c_0$ cannot be defined either, so it truncates immediately.

Stability of continued-fraction evaluation

Continued-fraction evaluation is stable with the backward (bottom-up) method: start from the deepest level and, working upward, fold up

$$\text{acc} \leftarrow 1 + \frac{a_i z}{\text{acc}}$$

and finally set $f(z) = b_0 / \text{acc}$ (the J-fraction uses an analogous recurrence including the second-degree term). Compared with the naive approach of evaluating the numerator and denominator polynomials separately to large values and then taking their ratio, this method keeps intermediate values near $1$ and is robust against cancellation. sangi's evalSFraction / evalJFraction implement this backward evaluation.

Rational Interpolation (Thiele / Cauchy)

Padé approximation matched the Taylor coefficients at a single point. By contrast, rational interpolation seeks a rational function that matches the values $y_i = f(x_i)$ at several distinct sample points $x_0, x_1, \ldots, x_{M+N}$. It is also called multipoint Padé, Newton-Padé, or Cauchy interpolation.

Thiele form (continued-fraction interpolation by inverse differences)

The classical construction is Thiele's continued-fraction interpolation, which uses the rational analogue of Newton's divided differences, the inverse (reciprocal) differences $\varphi[\cdot]$, to write

$$R(x) = y_0 + \cfrac{x - x_0}{\varphi[x_0,x_1] + \cfrac{x - x_1}{\varphi[x_0,x_1,x_2] + \cfrac{x - x_2}{\ddots}}}$$

in continued-fraction form. The inverse differences are computed by the recurrence

$$\varphi[x_i, x_{i+1}] = \frac{x_i - x_{i+1}}{y_i - y_{i+1}}, \qquad \varphi[x_0,\ldots,x_k] = \frac{x_{k-1} - x_k}{\varphi[x_0,\ldots,x_{k-2},x_{k-1}] - \varphi[x_0,\ldots,x_{k-2},x_k]}$$

and play, at multiple points, the role the S-fraction and J-fraction play in Padé. This rational extrapolation is used in the Bulirsch-Stoer ODE solver and in extrapolation for numerical integration (refining the samples and extrapolating the $h \to 0$ limit by a rational function).

sangi's implementation — the Cauchy-type linear system

Instead of explicitly building Thiele's inverse differences, sangi's rationalInterpolation(x, y, M, N) takes the Cauchy form, solving the interpolation conditions directly as a linear system. At each sample point $P(x_i) - y_i\, Q(x_i) = 0$, that is,

$$\sum_{k=0}^{M} p_k\, x_i^k \; - \; y_i \sum_{k=1}^{N} q_k\, x_i^k \; = \; y_i \qquad (i = 0, 1, \ldots, M+N)$$

is solved as a system of $M+N+1$ equations under the $q_0 = 1$ normalization by LU decomposition (exactly $M+N+1$ sample points are required). It returns the coefficient lists numerator / denominator as a RationalInterpResult<T>, and evalRational evaluates by Horner's method.

Differences from (single-point) Padé and things to watch

• Input: Padé takes Taylor coefficients at one point; rational interpolation takes sample values at distinct multiple points.
• Repeated samples: if the $x_i$ have repetitions, the linear system becomes ill-conditioned (a Hermite-type formulation is really needed then). Take the samples at distinct points. When you want the Taylor-coefficient version collapsed at $x_i = 0$, use padeApprox.
• Unattainable points: in general, rational interpolation has degenerate configurations (unattainable points) that no rational function can interpolate. If the denominator becomes zero at a sample point (a pole), the solution is invalid unless the numerator is simultaneously zero there.

The construction is based on the rational-interpolation framework of Stoer & Bulirsch (§2.2.4).

Padé for the Matrix Exponential (Scaling & Squaring)

Padé approximation applies not only to scalar functions but also to matrix-valued functions. The representative example is the matrix exponential

$$e^A = \sum_{k=0}^{\infty} \frac{A^k}{k!}$$

The naive approach of directly truncating the series loses accuracy when $\|A\|$ is large, because the terms first grow enormous and then cancel (cancellation). The practical, standard solution is scaling & squaring:

1. Scaling: look at $\|A\|$, choose $s$, and make the norm of $A/2^s$ sufficiently small.
2. Padé: apply a diagonal Padé approximant $e^{A/2^s} \approx R_{m/m}(A/2^s)$ to the scaled-down $A/2^s$. The diagonal Padé of the exponential has explicitly given numerator and denominator coefficients, and can be evaluated with a single matrix solve (one linear system).
3. Squaring: recover $e^A = \left(e^{A/2^s}\right)^{2^s}$ by squaring $s$ times.

sangi's matrixExpPade(A) is a thin wrapper around the library's matrix-exponential routine algorithms::expm. expm uses $[13/13]$ diagonal Padé as the core of scaling & squaring, and has a higher-accuracy Schur-Parlett path as needed (the Higham 2005 scheme). The argument $m$ is kept for source compatibility, but the Padé degree is chosen inside the routine, so it is ignored. A non-square input throws an exception.

Comparison Table

Rough guidance: if you have Taylor coefficients and want to raise the accuracy, use Padé approximation; if you want to stabilize that evaluation or obtain a Padé sequence by truncation, use a continued fraction; if you have discrete samples rather than Taylor coefficients, use rational interpolation; if you need the matrix exponential, use scaling & squaring.

References

• Baker, G. A., & Graves-Morris, P. (1996). Padé Approximants. 2nd ed. Encyclopedia of Mathematics and Its Applications, Cambridge University Press.
• Cuyt, A., & Wuytack, L. (1987). Nonlinear Methods in Numerical Analysis. North-Holland.
• Rutishauser, H. (1957). Der Quotienten-Differenzen-Algorithmus. Birkhäuser.
• Stoer, J., & Bulirsch, R. (2002). Introduction to Numerical Analysis. 3rd ed. Springer. §2.2 (rational interpolation, Bulirsch-Stoer extrapolation).
• Higham, N. J. (2005). "The scaling and squaring method for the matrix exponential revisited". SIAM Journal on Matrix Analysis and Applications, 26(4), 1179–1193.