We introduce a simple and performant approach for rapidly and accurately performing polyharmonic spline (PHS) interpolation using a combination of the fast multipole method (FMM) from computational electrostatics and the Vecchia approximation from the Gaussian process and sparse approximate inverse literatures. Using basic properties about Hadamard products and low-rank matrices, we demonstrate that an FMM with two kernels, the logarithmic and distance kernels, results in fast PHS interpolation for all orders. Furthermore, we demonstrate the exceptional performance of sparse inverse approximation methods with the Matérn covariance model for preconditioning. Combined with careful management of disallowed subspaces, we describe a procedure for obtaining prediction weights using preconditioned conjugate gradient that converges in less than $15$ iterations, even for problem sizes with over one million points. As a result, thin-plate spline interpolation---a particularly popular method that does not require parameter tuning---that matches the fully dense $\mathcal{O}(n^3)$ computation in accuracy can be done at the cost of approximately $50-60$ FMMs. A high-performance software library for odd-order PHS interpolation in two dimensions is made available as a companion to this work.
🌊
