🔢 mathematics

Linear-cost Polyharmonic Spline Interpolation of Arbitrary Degree

This paper introduces a highly efficient method for polyharmonic spline interpolation of arbitrary degree that combines the fast multipole method with sparse inverse approximations and preconditioned conjugate gradients to achieve linear-cost computation and rapid convergence for large-scale datasets while maintaining the accuracy of traditional dense solvers.

Original authors: Christopher J. Geoga, Michael O'Neil

Published 2026-08-13

📖 8 min read🧠 Deep dive

CC BY 4.0

Original authors: Christopher J. Geoga, Michael O'Neil

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a cartographer trying to draw a perfect map of a bumpy, mountainous landscape, but you only have a handful of scattered weather stations reporting the height of the ground. Your goal is to guess the elevation of every single spot in between those stations so you can build a smooth, continuous surface. This is the heart of a field called "interpolation," a branch of mathematics used everywhere from weather forecasting to computer graphics. The tricky part is that the more data points you have, the harder the math gets. In fact, for many traditional methods, doubling your data doesn't just double the work; it multiplies it by a huge number, making it impossible to solve on a normal computer if you have millions of points.

To solve this, scientists often use a tool called a "polyharmonic spline." Think of this as a magical, stretchy rubber sheet that you pin down at your known data points. The sheet naturally settles into a shape that connects all the dots smoothly. The problem is that calculating exactly how this rubber sheet bends requires solving a massive, tangled web of equations. Usually, this takes so much computer power that it's like trying to count every grain of sand on a beach by hand. However, there are two clever tricks in the scientific toolbox that can speed things up. The first is the "Fast Multipole Method" (FMM), which is like a super-efficient way of grouping distant friends together so you don't have to talk to every single person individually to send a message. The second is the "Vecchia approximation," which is a way of guessing the answer by only looking at your closest neighbors, assuming that people far away don't influence you much.

This paper introduces a new, super-fast way to draw that rubber sheet map, even when you have over a million data points. The authors, Christopher J. Geoga and Michael O'Neil, combined those two clever tricks—the grouping method and the neighbor-guessing method—with a few new mathematical shortcuts. They found that by treating the problem like a physics puzzle involving electric charges and using a specific type of "pre-conditioner" (a mathematical warm-up exercise that helps the computer solve the puzzle faster), they could get the answer almost instantly. Their method is so efficient that it can handle a million points in under 15 seconds on a regular laptop, a task that would usually take hours or days. They also showed that this approach is incredibly accurate, matching the results of the slow, perfect methods almost exactly, without needing to tweak any settings. It's a bit like finding a shortcut through a dense forest that gets you to the same destination as the long, winding path, but in a fraction of the time.

The Magic of the Stretchy Sheet

At the core of this work is a problem that sounds simple but gets messy fast: how do you fill in the blanks between data points? The authors use a method called Polyharmonic Spline (PHS) interpolation. Imagine you have a sheet of rubber and you pin it down at specific locations where you know the height. The sheet naturally curves to connect them. The math behind this involves a "kernel matrix," which is just a giant spreadsheet showing how every single point talks to every other point.

The trouble is, this spreadsheet is "dense," meaning every cell has a number in it. If you have 1,000 points, you have a million cells to calculate. If you have a million points, you have a trillion trillion cells. Traditional computers would need to do a cubic amount of work ( $O(n^3)$ ) to solve this, which is why it's usually impossible for huge datasets.

The authors' first big insight is that they don't need to calculate every single cell directly. Instead, they realized that the math behind the rubber sheet can be broken down into two simpler parts. One part is a "core" kernel, which is like a basic building block (either a logarithm or a simple distance). The other part is a low-rank matrix, which is a fancy way of saying it has a lot of repeating patterns that can be simplified. By using a mathematical trick called the Hadamard product (which is just multiplying matrices element-by-element), they showed that they could compute the whole thing by just running a fast algorithm on that simple "core" building block.

The Fast Multipole Method: Grouping the Crowd

To speed up the calculation of that "core" building block, the authors use the Fast Multipole Method (FMM). Imagine you are at a massive concert and you need to shout a message to everyone in the crowd. If you shout to every single person one by one, it takes forever. But, if you group people into clusters, you can shout to the center of a cluster, and the sound carries to everyone in that group.

The FMM does exactly this for math. It organizes the data points into a tree-like structure (a quadtree). If a group of points is far away from the point you are calculating, the algorithm treats the whole group as a single "super-point" with a combined effect. This turns a problem that would take forever into one that scales linearly ( $O(n)$ ). If you double the number of points, the time only doubles, rather than exploding. The authors adapted this method, originally used for electrostatics (calculating how electric charges push and pull on each other), to handle the specific math of the rubber sheet.

The Preconditioner: Warming Up the Engine

Even with the fast grouping trick, the computer still needs to solve a system of equations to find the exact shape of the rubber sheet. This is where the "preconditioner" comes in. Think of the computer solver as a car trying to drive up a steep, winding hill. If the hill is too steep or twisty, the car might stall or take forever. A preconditioner is like a road crew that smooths out the path, making the hill easier to climb so the car can zoom to the top.

The authors propose a new, incredibly fast preconditioner based on the "Vecchia approximation." This method assumes that a point is mostly influenced by its nearest neighbors, not by points on the other side of the world. By using a statistical model called the Matérn covariance (which describes how things smooth out over distance), they can build a sparse matrix—a spreadsheet where most of the cells are zero. This sparse matrix is easy to compute and acts as a perfect warm-up for the solver.

The authors found that this specific combination works wonders. In their tests, the computer solver (a method called Preconditioned Conjugate Gradient) converged in fewer than 15 iterations, even for datasets with over one million points. This means the car didn't just climb the hill; it flew up it.

The Results: Speed Meets Accuracy

The paper puts this new method to the test with several experiments. First, they compared it to older methods. They found that while other approaches might work for small datasets, they often fail to control the number of steps needed as the data gets bigger. The new Vecchia-based preconditioner, however, kept the number of steps low and steady, regardless of size.

They also tested the accuracy. In one experiment, they tried to predict a complex function that had both smooth waves and a sharp, jagged spike. The new method produced errors that were virtually identical to the "exact" method (the slow, perfect one), proving that the shortcuts didn't sacrifice quality.

Perhaps the most impressive demonstration was a real-world test using sea surface temperature data from the Pacific Ocean. They had about 58,000 measurements with some missing due to "cloud cover" (simulated gaps). Using their method, they filled in the missing data in just 5 seconds with a very low error rate. In contrast, a traditional method using the same statistical model took over 400 seconds and actually performed worse. This highlights a key feature of their approach: because the polyharmonic spline is "scale-invariant," it doesn't need to be tuned or adjusted for different sizes of data, making it a "plug-and-play" solution that just works.

Why This Matters

The authors conclude that this approach offers a "truly end-to-end linear-cost" solution. This means that as your data grows, the time it takes to solve the problem grows at a manageable, steady pace. They have even released a software library that allows others to use this method for 2D data. While they focused on 2D and specific orders of the spline, they suggest that the same logic could work for 3D and other variations in the future.

In short, Geoga and O'Neil have taken a problem that was previously too heavy for most computers to lift and made it light enough to carry in a backpack. By combining the speed of grouping distant points with the efficiency of neighbor-based guessing, they've created a tool that can map the world, one million points at a time, in the blink of an eye.

Technical Summary: Linear-cost Polyharmonic Spline Interpolation of Arbitrary Degree

Problem Statement
Polyharmonic spline (PHS) interpolation is a widely used kernel interpolation method that decomposes an estimator into a polynomial trend and a dense kernel-based remainder. The core challenge in applying PHS to large datasets ( $n$ ) is the computational cost associated with the kernel matrix $M$ , which is dense and conditionally positive-definite. Direct solution of the augmented linear system required for PHS weights scales as $O(n^3)$ , rendering it prohibitively expensive for large $n$ . While iterative methods (e.g., GMRES, PCG) can reduce this to $O(n^2)$ by avoiding explicit factorization, they rely heavily on the existence of effective preconditioners to control iteration counts. Previous attempts to accelerate PHS using the Fast Multipole Method (FMM) or specific preconditioning strategies have struggled to maintain linear cost or control iteration counts under standard convergence metrics for arbitrary spline orders.

Methodology
The authors propose a unified approach combining the Fast Multipole Method (FMM) for matrix-vector products and the Vecchia approximation for preconditioning, achieving linear-cost interpolation for arbitrary PHS orders.

Linear-cost Matrix-Vector Products via FMM:
The authors observe that the PHS kernel matrix $M$ for any order $k$ can be decomposed using Hadamard products involving low-rank matrices and a "core" kernel. Specifically, for even $k$ , the term $\|x_j - x_k\|^k$ is separable, allowing $M$ to be expressed as a sum of Hadamard products between diagonal matrices and a core kernel matrix $M_{\log}$ (where entries are $\log \|x_j - x_k\|$ ). Since $\log r$ is the Green's function for the 2D Poisson equation, the matrix-vector product $Mv$ can be computed using a single FMM routine for the logarithmic kernel (and a similar approach for the linear kernel $\|r\|$ for odd orders). This reduces the complexity of applying $M$ to $O(n)$ (specifically $O(p^2 n)$ where $p$ is the expansion order).
Handling Scale Invariance and Subspaces:
PHS interpolation is scale-invariant. The authors utilize this property to rescale data points to a domain convenient for the FMM without affecting the interpolation weights. They employ a projection strategy using the orthogonal complement $Q_\perp$ of the polynomial basis $B$ . This transforms the indefinite augmented system into a strictly positive-definite system of size $(n-S) \times (n-S)$ , allowing the use of the Preconditioned Conjugate Gradient (PCG) method.
Vecchia Approximation Preconditioning:
To accelerate PCG convergence, the authors introduce a sparse preconditioner based on the Vecchia approximation of the Matérn covariance model. They leverage the observation that the thin-plate spline kernel (and higher-order PHS) shares asymptotic behavior with the Matérn kernel (specifically with smoothness parameter $\nu=1$ ) near the origin.
- The preconditioner is constructed as $P = Q_\perp^T \tilde{T}\tilde{T}^T Q_\perp$ , where $\tilde{T}\tilde{T}^T$ is a sparse approximation of the inverse Matérn covariance matrix derived via Vecchia's approximation (conditioning each point on a small set of neighbors).
- This approach allows the preconditioner to be assembled in $O(n)$ time with embarrassingly parallel operations.
- The authors note a trade-off: rescaling the domain and tuning the Matérn lengthscale parameter $\alpha$ balances the accuracy of the kernel approximation against the conditioning of the Vecchia approximation.

Key Contributions

Generalization of FMM for PHS: The paper demonstrates that FMM with just two kernels (logarithmic and distance) suffices for matrix-vector products of any PHS order $k$ , simplifying previous approaches that required specific kernels for specific degrees.
Novel Preconditioner: The introduction of a Vecchia-based sparse preconditioner using the Matérn model, which is shown to be significantly more effective than prior approximate cardinal function approaches.
End-to-End Linear Cost: The combination of these techniques yields a procedure that converges in fewer than 15 iterations (often <10) even for $n > 10^6$ , achieving true linear-cost interpolation.
Software Implementation: A high-performance Julia library (FastPolyharmonicSplines.jl) is provided for odd-order PHS interpolation in 2D.

Results
Numerical experiments on synthetic and real-world data demonstrate the efficacy of the method:

Iteration Control: Unlike previous methods where iteration counts grew with $n$ , the Vecchia preconditioner maintains iteration counts below 10 for thin-plate splines ( $k=2$ ) and around 15 for $k=4$ , regardless of dataset size up to $n=10^6$ .
Accuracy: The iterative solver produces prediction errors identical to exact dense linear algebra methods (using Bunch-Kaufman factorization) within machine precision when tolerances are tight.
Performance: For $n=10^6$ points, the method achieves interpolation in under 15 seconds on a standard laptop. The cost is approximately equivalent to 50–60 FMM evaluations.
Real Data Application: In a sea surface temperature anomaly interpolation task ( $n \approx 58,000$ ), the scale-invariant thin-plate spline method outperformed a tuned Matérn Gaussian process model (RMSE 0.287 vs. 0.398) while requiring significantly less runtime (5s vs. 418s), highlighting the benefit of PHS when parameter tuning is not desired.

Significance and Claims
The paper claims to provide a "simple and performant approach" that removes the primary computational bottleneck of PHS interpolation. By decoupling the matrix-vector product cost from the spline order and providing a robust, linear-cost preconditioner, the authors enable the use of high-order PHS interpolation for massive datasets without the need for parameter tuning (a key advantage of PHS). The work establishes that thin-plate spline interpolation, previously limited by $O(n^3)$ costs, can now be performed at a cost of roughly 50–60 FMMs while matching the accuracy of dense solvers. The authors modestly note that while the current library focuses on 2D odd-order PHS, the theoretical framework extends to arbitrary orders and higher dimensions, representing a significant step toward scalable, parameter-free kernel interpolation.