Scattered Data Interpolation and the Conditioning of Radial Basis Functions

Measurement does not arrive on a grid. Weather stations, sensor buoys, scattered probe points, and particle tracers produce positions and values with no adjacency relation between them, and the visualization needs a value everywhere, including at positions where nothing was measured. The reconstruction problem is therefore stated without topology: given $N$ pairs $(\mathbf{x}_i, f_i)$ in $\mathbb{R}^d$, produce a function that agrees with the data at those points and behaves sensibly between them.

Grid methods answer this by inventing a grid and interpolating on it, which moves the problem rather than solving it. Meshless methods keep the points and choose a basis whose coefficients are determined by the data. The design tension is stark: the methods that produce the smoothest, most pleasing surfaces are the ones whose linear systems are the most delicate to solve, and the methods that are numerically tame produce surfaces with visible artefacts.

The interpolation matrix is dense, and that is the whole cost

Radial basis function interpolation takes the form

$$s(\mathbf{x}) = \sum_{j=1}^{N} w_j , \varphi(\lVert \mathbf{x} - \mathbf{x}_j \rVert) + p(\mathbf{x}),$$

where $\varphi$ is a scalar function of distance and $p$ is a low-degree polynomial, present when the kernel is only conditionally positive definite. Imposing $s(\mathbf{x}_i) = f_i$ at every data point gives a linear system of size $N + m$, with $m$ the number of polynomial terms. The matrix’s entries depend only on pairwise distances, so it has no sparsity structure at all: every data point influences every basis function.

That single fact sets the cost. Forming the matrix requires $O(N^2)$ distance evaluations; solving it by a dense direct method requires $O(N^3)$ work and $O(N^2)$ storage. For $N = 10^4$ double-precision entries, the matrix alone is 800 MB, and a factorization is out of reach for interactive use. For $N = 10^6$, no direct approach is viable at all. The usual remedies are not improvements to the solve but changes to the problem: use a compactly supported kernel so the matrix becomes sparse, use a local formulation so the solve is many small systems, or use a fast multipole or hierarchical method to apply the matrix without forming it.

The kernel choice determines both the quality of the surface and the conditioning of that matrix. Four families dominate.

The Gaussian, $\varphi(r) = e^{-(\varepsilon r)^2}$, is infinitely smooth and produces very accurate reconstruction of smooth data. Its matrix becomes numerically singular quickly as $N$ grows or as the shape parameter $\varepsilon$ shrinks toward a flat kernel, because the columns of the matrix become nearly linearly dependent. This is not a defect of any particular solver; the interpolation problem itself becomes ill-posed in that limit, and the coefficients that solve it are enormous with alternating signs.

The multiquadric, $\varphi(r) = \sqrt{1 + (\varepsilon r)^2}$, behaves similarly: smooth, accurate, and increasingly ill-conditioned as the shape parameter decreases. The thin plate spline, $\varphi(r) = r^2 \log r$ in two dimensions and $\varphi(r) = r$ in three, is the fundamental solution of a biharmonic operator and minimizes a bending energy functional over all interpolants. It has no shape parameter to tune, which removes one failure mode, but its matrix is dense and grows in condition number with $N$ regardless.

Compactly supported kernels, of which Wendland’s are the standard construction, are positive definite on $\mathbb{R}^d$ for a stated smoothness and vanish beyond a radius. Truncating at that radius makes the matrix sparse, with a bandwidth set by the support radius relative to point spacing. The cost is a visible seam where the support ends unless the support radius spans enough neighbours, typically eight to twenty per dimension, and the surface quality is lower than a global smooth kernel’s for the same data.

Shape parameter, Runge phenomena, and the flat-kernel limit

The shape parameter controls how quickly the basis function decays with distance. Large $\varepsilon$ makes each basis function peaked and local in effect; small $\varepsilon$ makes it broad and nearly polynomial. Accuracy on smooth data typically improves as $\varepsilon$ decreases, while conditioning worsens, and the two curves cross at a point that depends on $N$, on the point distribution, and on the floating-point precision in use.

This is the same trade that appears in polynomial interpolation as the Runge phenomenon, and the connection is not superficial. In the flat limit, the RBF interpolant converges to a polynomial interpolant of high degree, and polynomial interpolation on equispaced points diverges near the interval ends. RBF interpolation on scattered points avoids the specific equispaced failure because the points are irregularly placed, but the underlying ill-conditioning remains, and it appears as high-frequency oscillation in the reconstructed field between data points, which a visualization renders as ripples that have no counterpart in the data.

In double precision, the practical consequence is that the shape parameter cannot be chosen purely for accuracy. RBF-specific stabilization methods exist: RBF-QR rewrites the basis in a numerically stable form for near-flat kernels; replacement of the ill-conditioned system by a regularized least-squares problem (RBF regression with smoothing) accepts a residual at the data points; and contour integral or arbitrary-precision evaluation extends the usable range at a cost of many times the arithmetic. All of them exist because the straightforward solve fails, not because the mathematics is wrong.

import numpy as np

def rbf_interpolate(x, f, query, eps, kernel="mq"):
    """Global RBF interpolation in d dimensions. x: (N, d), f: (N,)."""
    def phi(r):
        if kernel == "mq":
            return np.sqrt(1.0 + (eps * r) ** 2)
        if kernel == "gauss":
            return np.exp(-(eps * r) ** 2)
        if kernel == "tps":
            r = np.maximum(r, 1e-14)          # r^2 log r -> 0 as r -> 0
            return r ** 2 * np.log(r)
        raise ValueError(kernel)

    def dist(a, b):
        # pairwise squared distances via the identity, avoiding an (N,M,d) cube
        return np.sqrt(np.maximum(
            (a * a).sum(1)[:, None] + (b * b).sum(1)[None, :] - 2 * a @ b.T, 0.0))

    A = phi(dist(x, x))
    # TPS is only conditionally positive definite: add a linear polynomial block
    if kernel == "tps":
        P = np.hstack([np.ones((len(x), 1)), x])
        M = np.block([[A, P], [P.T, np.zeros((P.shape[1], P.shape[1]))]])
        rhs = np.concatenate([f, np.zeros(P.shape[1])])
        sol = np.linalg.solve(M, rhs)
        w, c = sol[:len(x)], sol[len(x):]
        Pq = np.hstack([np.ones((len(query), 1)), query])
        return phi(dist(query, x)) @ w + Pq @ c
    w = np.linalg.solve(A, f)
    return phi(dist(query, x)) @ w

The polynomial augmentation in the thin-plate branch is not decoration. Conditional positive definiteness guarantees invertibility only on the subspace of coefficient vectors that satisfy the polynomial side conditions, so the kernel alone does not determine a unique interpolant; the augmented block supplies those conditions and lets the reconstruction represent a linear field exactly.

Local stencils trade smoothness for tractability

The dense solve’s cost is a property of the global formulation, not of the kernel. A local formulation builds one small system per query point from its nearest neighbours, typically ten to fifty in three dimensions, and evaluates a low-degree polynomial or a local RBF fit there. The work per query is $O(k^3)$ with $k$ the stencil size, independent of $N$, and the whole reconstruction is $O(N k^3)$ for $N$ queries.

The cost of locality is smoothness. A global thin-plate spline is the minimizer of a bending energy over the whole domain, so the surface has continuous curvature by construction. Local fits agree with the data and with each other only approximately, and where stencils differ the surface can have small discontinuities in slope or value. The standard mitigation is to blend overlapping local fits with partition-of-unity weights, which restores continuity at the price of a weighting function that must sum to one everywhere and whose support must be large enough to overlap.

Local methods also have a well-behaved conditioning story that global ones lack, which is why they dominate practical meshless solvers. A stencil of twenty points in three dimensions, fitted by a quadratic polynomial with ten terms, has an overdetermined system whose normal equations are modest and well-scaled if the points are not nearly coplanar. Degenerate stencils, with all points on a plane or a line, are a real and common failure: a sensor array on a surface, or points sampled along a coastline, produce rank deficiency that a general solver reports as singularity and that a careful implementation must detect and handle by falling back to a lower-degree fit.

Delaunay interpolation: cheap, local, and C0

The alternative to a smooth basis is a piecewise-linear interpolant over a triangulation. In two dimensions, the Delaunay triangulation maximizes the minimum angle over all triangulations of the same point set, which makes it the preferred choice for interpolation because it avoids the thin slivers that a greedy triangulation produces. The interpolant on each triangle is the unique linear function through its three vertex values, so the reconstruction is continuous and its gradient is not.

Cost is $O(N \log N)$ for the construction and $O(\log N)$ per point location with a suitable structure, or $O(1)$ by walking from the previously located triangle when queries are spatially coherent, which they are for a raster image. The memory is linear in $N$. Against a global RBF solve, this is a different regime entirely: millions of points are routine.

The visual consequences are exactly what the continuity class predicts. Contours of a Delaunay-interpolated field are piecewise linear, so an iso-line has corners at cell boundaries; contours of an RBF field are smooth. The Delaunay surface does not overshoot, because a linear function on a triangle is bounded by its vertex values, and its extrema lie at data points. The smooth surface can overshoot between points, and on sparse, irregular data that overshoot appears as haloes around isolated high values.

In three dimensions the Delaunay triangulation becomes a tetrahedralization, and the favourable angle property does not survive: slivers with near-zero volume are possible and are not removed by the Delaunay criterion alone. Interpolation on a sliver is numerically unstable because the barycentric coordinates involve division by a volume near zero. This is the reason three-dimensional scattered interpolation often prefers local RBF or moving least squares over tetrahedralization.

Choosing, given what the picture needs

The decision is driven by the derivative continuity the visualization requires and by the number of points. A contour plot needs only a continuous field, so piecewise linear is adequate; a surface rendering shaded with normals needs a continuous gradient, or the shading shows the triangulation; a computation of curvature, of streamlines through the reconstructed field, or of a gradient-based feature needs a continuous second derivative.

Shepard’s method, with weights proportional to an inverse power of distance, sits at the cheap end: no linear system, $O(N)$ per query, exact at data points and continuous everywhere except at them, where the derivative is zero, producing flat spots at every sample. Modified Shepard methods restrict the sum to a local neighbourhood to reduce the cost and improve the shape, and they remain the simplest thing that produces a smooth-looking surface.

Once the field exists, its rendering raises a separate set of questions. The choice of color ramp for the reconstructed values is constrained as described in Perceptually Uniform Colormaps and the Failure of Rainbow Ramps; extracting a level set from it requires a grid or a tetrahedralization, as described in Marching Cubes and the Case Table; and integrating through a vector field reconstructed this way, where the interpolation error changes trajectories rather than values, is the subject of Streamline Integration in Vector Fields.