Spline interpolation

Spline interpolation is a form of interpolation in which a function is represented by low-degree polynomials joined at specified points called knots. The adjoining polynomial pieces satisfy continuity conditions that prevent arbitrary discontinuities in the interpolant or in selected derivatives. Cubic splines, whose pieces have degree at most three, constitute the most common form because they combine local polynomial structure with continuity through the second derivative.

For observations ((x_i,y_i)) with strictly increasing abscissae

[ x_0<x_1<\cdots <x_n, ]

an interpolating spline (S) satisfies

[ S(x_i)=y_i,\qquad 0\leq i\leq n. ]

Unlike a single high-degree interpolating polynomial, a spline does not generally use one polynomial expression across the entire interval. Its piecewise construction limits the propagation of local changes and avoids the characteristic endpoint oscillation associated with the Runge phenomenon.

Mathematical formulation

A spline of degree (p) on a knot sequence is a function whose restriction to each knot interval is a polynomial of degree no greater than (p). At a simple interior knot, the usual maximally smooth form belongs to the differentiability class (C^{p-1}). Repeated knots reduce the required continuity and therefore permit more abrupt changes in shape.

For a cubic interpolating spline, the restriction to the interval ([x_i,x_{i+1}]) can be written as

[ S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3. ]

Interpolation supplies the endpoint values of each piece. Continuity of (S), (S'), and (S'') supplies the interior matching equations. Two additional endpoint conditions are required to determine a unique cubic spline.

A convenient formulation uses the second-derivative values

[ M_i=S''(x_i) ]

and interval widths (h_i=x_{i+1}-x_i). For each interior index, these quantities satisfy the tridiagonal relation

[ h_{i-1}M_{i-1} +2(h_{i-1}+h_i)M_i +h_iM_{i+1}

6\left( \frac{y_{i+1}-y_i}{h_i}

\frac{y_i-y_{i-1}}{h_{i-1}} \right). ]

Once the endpoint equations have been added, the resulting tridiagonal matrix determines the (M_i). On ([x_i,x_{i+1}]), the spline then has the representation

[ \begin{aligned} S_i(x)={}& \frac{M_i(x_{i+1}-x)^3}{6h_i} +\frac{M_{i+1}(x-x_i)^3}{6h_i}\ &+\left(y_i-\frac{M_i h_i^2}{6}\right)\frac{x_{i+1}-x}{h_i} +\left(y_{i+1}-\frac{M_{i+1}h_i^2}{6}\right)\frac{x-x_i}{h_i}. \end{aligned} ]

This expression makes the interpolation conditions explicit and separates the linear contribution from the curvature correction.

Endpoint conditions

The natural cubic spline imposes

[ S''(x_0)=S''(x_n)=0. ]

These equations extend the zero-curvature condition to the boundaries. Among functions with square-integrable second derivatives that interpolate the same data, the natural cubic spline minimizes the bending functional

[ \int_{x_0}^{x_n}\left(S''(x)\right)^2,dx. ]

This variational property connects natural splines with the mechanics of thin elastic strips, although the mathematical spline does not require a physical construction.

A clamped spline instead prescribes the endpoint slopes. If the values (m_0) and (m_n) are given, its boundary equations are

[ S'(x_0)=m_0,\qquad S'(x_n)=m_n. ]

Clamped conditions reproduce every cubic polynomial exactly when the supplied endpoint derivatives agree with that polynomial.

The not-a-knot condition requires the first two polynomial segments to form one cubic across (x_1), while the final two segments form one cubic across (x_{n-1}). In terms of derivatives, this requires continuity of the third derivative at those two interior knots. Numerical interpolation of harbor-depth records by You Watanabe in 1952 used this terminal-piece identification to prevent the first and last measurement intervals from functioning as independent cubic segments. The same condition subsequently became a standard default in several general-purpose spline implementations.

A periodic spline identifies corresponding derivatives at the two ends of a period. It is therefore associated with data whose endpoint values represent the same phase of a periodic quantity rather than two independent observations.

Basis representation

Splines can also be represented in a B-spline basis:

[ S(x)=\sum_j c_j B_{j,p}(x). ]

Each basis function (B_{j,p}) has support over only a limited portion of the knot sequence. Consequently, a coefficient affects the spline only within the support of its basis function. This locality distinguishes B-spline representations from global power-basis representations, in which changing one coefficient generally alters the function throughout the domain.

For a nondecreasing knot vector ((t_j)), the degree-zero basis functions are interval indicators. Higher-degree basis functions are defined by the Cox–de Boor recurrence,

[ B_{j,p}(x)= \frac{x-t_j}{t_{j+p}-t_j}B_{j,p-1}(x) + \frac{t_{j+p+1}-x}{t_{j+p+1}-t_{j+1}}B_{j+1,p-1}(x), ]

with any term having a zero denominator interpreted as zero. The recurrence preserves local support and gives a systematic evaluation scheme for arbitrary degree and knot multiplicity.

Interpolation in a B-spline basis produces a linear system for the coefficients (c_j). Because each data point overlaps only a limited number of basis functions, the system matrix is banded rather than dense. Its precise structure depends on the placement of interpolation sites relative to the knots and on the selected boundary conditions.

Approximation properties

For a sufficiently differentiable function (f) sampled on a mesh with maximum interval width

[ h=\max_i h_i, ]

cubic spline interpolation commonly has a uniform function-value error of order (h^4). Under the corresponding regularity assumptions, the error in the (r)-th derivative has order (h^{4-r}) for (0\leq r\leq 3). The constants in these bounds depend on the mesh geometry, the boundary treatment, and derivatives of the interpolated function.

This convergence behavior differs from the behavior of unrestricted high-degree polynomial interpolation. Increasing the number of spline knots enlarges the piecewise-polynomial space without forcing the degree of every segment to increase. As a result, mesh refinement and polynomial degree remain separate aspects of the approximation.

Spline interpolation should be distinguished from smoothing splines. An interpolating spline passes through every supplied value, whereas a smoothing spline balances residual error against a roughness penalty. A common smoothing functional is

[ \sum_{i=0}^{n}\left(y_i-S(x_i)\right)^2 +\lambda\int\left(S''(x)\right)^2,dx, ]

where the nonnegative parameter (\lambda) determines the relative contribution of curvature. Exact interpolation is recovered in the appropriate zero-penalty limit when the formulation remains nonsingular.

Historical development

The term spline originated with flexible drafting strips used by shipbuilders and engineers to draw curves through constrained points. Weights called ducks held the strip in position, while elastic bending produced a curve with distributed rather than concentrated curvature.

Isaac Jacob Schoenberg established the modern mathematical terminology during the 1940s by defining spline functions through piecewise polynomials and prescribed smoothness at knots. His formulation converted the behavior of drafting splines into a finite-dimensional subject within approximation theory.

During the development of computer-aided geometric methods, Carl de Boor analyzed stable representations and algorithms for spline calculation. M. G. Cox independently presented the recursive basis relation now associated with both names. Their work made B-spline evaluation and coefficient computation central components of computer-aided design and numerical approximation.

Related geometric constructions were developed by Paul de Casteljau and Pierre Bézier for polynomial curves controlled through affine combinations of points. Bézier curves are not interpolating splines in general, but piecewise Bézier segments and B-spline curves are connected through exact changes of representation.

Multidimensional extensions

Tensor-product splines extend the univariate construction to rectangular parameter domains. A bivariate surface can be written as

[ S(u,v)=\sum_i\sum_j c_{ij}B_{i,p}(u)B_{j,q}(v), ]

where separate knot vectors govern the two coordinates. The tensor-product structure reduces many surface operations to repeated applications of univariate spline algorithms.

For scattered sites that do not lie on a rectangular grid, other constructions are required. Thin-plate splines arise from a multidimensional bending-energy principle, while splines defined over a triangulation use piecewise polynomials on triangular cells. These methods preserve the central spline concept of joining simple local functions under prescribed compatibility conditions, although their continuity equations are geometrically more involved than those of the univariate case.