Scattered-Data Interpolation

This documentation was generated with the assistance of AI. Please report any inaccuracies.

Two-Dimensional Grid Interpolation assumes function values are tabulated on a complete rectangular grid. When they are instead known only at an arbitrarily scattered set of points in -dimensional space, a different family of methods is needed. This package implements the two most widely used general-purpose approaches — radial basis function (RBF) interpolation and kriging — plus Shepard interpolation, a cheap variant of RBF. All three are described in Press et al., 2007, Section 3.7 (page 139), and implemented as RadialBasisFunctionInterpolator, ShepardInterpolator, and KrigingInterpolator in the com.irurueta.numerical.interpolation package.

Both RBF interpolation and kriging are comparatively expensive: they need work to digest the data points once, followed by work per interpolated value (kriging can additionally supply an error estimate, at a further cost per value). This effectively limits their practical use to . Shepard interpolation trades accuracy for avoiding the setup entirely.

Radial basis function (RBF) interpolation

The idea: every known point "influences" its surroundings identically in every direction, according to some radial basis function of the distance from that point. Approximate the interpolating function everywhere as a linear combination of these, centered on every known point:

The weights are found by requiring the interpolation to be exact at every known data point — a linear system of equations in unknowns, solved once via LU decomposition:

A normalized variant (NRBF) instead requires the basis functions to sum to unity, replacing the two equations above with:

There is no consistent evidence that either variant is generally superior; both are available via a constructor flag.

Choosing a radial basis function

This package provides four pluggable kernels, each with a scale factor that should be chosen larger than the typical spacing between points but smaller than the "outer scale" of the function being interpolated:

Class Formula Notes

MultiQuadricRadialBasisFunction

The most commonly used choice; comparatively insensitive to the choice of .

InverseMultiQuadricRadialBasisFunction

Comparable accuracy to the multiquadric; decays to zero away from the data, useful when that is desired.

ThinPlateRadialBasisFunction

,

Physically motivated by the energy of a warped thin elastic plate.

GaussianRadialBasisFunction

Can reach very high accuracy for smooth functions with a well-tuned , but accuracy is correspondingly sensitive to that choice; also decays to zero away from the data.

Press et al., 2007 suggests trying these roughly in the order listed above, starting with the multiquadric, and always experimenting with the scale parameter — the difference in accuracy between a good and a poor choice of can span several orders of magnitude.

Shepard interpolation

Shepard interpolation is the special case of normalized RBF interpolation where as (and is finite, decreasing, for ). In that case the weights collapse to simply the function values themselves, and no linear system needs solving at all:

using the power-law kernel (typically ) — a nearness-weighted average that favors nearby points over distant ones. ShepardInterpolator implements exactly this, with no setup cost at all, at the price of being generally less accurate than a well-tuned RBF.

Interpolation by kriging

Kriging — named for mining engineer D.G. Krige — is a form of linear prediction (Gauss-Markov estimation / Gaussian process regression) that additionally supplies an error estimate for every interpolated value. It requires a variogram , the expected mean-square variation of the function over a given offset distance:

Writing and \(v^*_j \equiv v(\|\mathbf{x}^* - \mathbf{x}_j\|)\) for a query point , and defining the augmented vectors and matrix:

the kriging estimate and its variance are:

The extra row/column of ones (and the corresponding zero) make this an unbiased estimator by automatically correcting for a weighted average of the data. Computing the LU decomposition of once makes each subsequent interpolated value cost ; computing its variance costs a further .

For interpolation (as opposed to smoothing/fitting with known measurement error), an adequate variogram can usually be estimated automatically, without deep domain knowledge, using the simple power-law model that this package’s nested KrigingInterpolator.Variogram class fits by unweighted least squares over all pairs of data points:

with fixed (default , must satisfy ) and fitted from the data.

flowchart TD A["Scattered points (x_i, y_i)"] --> B{"Which method?"} B -->|"RBF"| C["Solve N x N linear system for weights w_i\n(LU decomposition, once)"] C --> D["interpolate: sum w_i * phi(|x - x_i|)"] B -->|"Shepard"| E["No setup: weights are just y_i"] E --> F["interpolate: weighted average using phi(r) = r^-p"] B -->|"Kriging"| G["Fit variogram v(r) = alpha * r^beta;\nsolve (N+1)x(N+1) system (LU, once)"] G --> H["interpolate: V*^T * V^-1 * Y (+ optional variance)"]

Library implementation

RadialBasisFunctionInterpolator and ShepardInterpolator both extend BaseRadialBasisFunctionInterpolator, which provides the shared Euclidean-distance and dimension bookkeeping; KrigingInterpolator is standalone, since its variogram-based system has a different shape (the extra unbiasedness row/column).

Class Setup cost Notes

RadialBasisFunctionInterpolator

Constructor takes a RadialBasisFunction kernel and a boolean flag for the normalized (NRBF) variant.

ShepardInterpolator

None

Constructor takes the power p (default 2.0).

KrigingInterpolator

Constructor takes a Variogram (or builds a default power-law one automatically); interpolate(pt, esterr) additionally returns the estimated standard error.
Kriging’s interpolated values are quite insensitive to the exact variogram model used, but its error estimates are considerably more sensitive — treat them as order-of-magnitude guidance rather than precise figures.

Example

final var rbf = new MultiQuadricRadialBasisFunction(scale);
final var interpolator = new RadialBasisFunctionInterpolator(points, values, rbf);
final var value = interpolator.interpolate(target);

// Kriging, with an error estimate
final var kriging = new KrigingInterpolator(points, values);
final var esterr = new double[1];
final var krigedValue = kriging.interpolate(target, esterr);