✓✓ Beats tuned baseline 2026

FMM-Accelerated Polyharmonic Neural Field Head

Implementation & benchmark of arXiv:2608.11462 — Linear-cost Polyharmonic Spline Interpolation of Arbitrary Degree

Usefulness7/10

Difficulty6/10

Novelty7/10

Architecture Scalability Speedup Memory Accuracy Mlp Embedding Inference speedup Memory Harmonic analysis Approximation theory Linear algebra Numerical analysis

Source paper: Linear-cost Polyharmonic Spline Interpolation of Arbitrary Degree arXiv:2608.11462 ⓘ · analyzed Sep 1, 2026

AI-generated research hypothesis, automatically tested. Not peer-reviewed.

Idea description

Attach a polyharmonic spline decoder to a coordinate MLP or use it as a standalone neural-field output head over a large set of spatial anchors. The decoder represents the output as a low-degree polynomial trend plus a PHS kernel expansion, while FMM evaluates all anchor-to-query interactions in approximately linear or near-linear cost. When coefficients must be fitted or periodically recalibrated, solve the constrained interpolation system with projected conjugate gradients and a sparse Vecchia/Matérn preconditioner instead of dense factorization.

Formulas

$$\hat f(\mathbf{x})=\sum_{s=1}^{S}\beta_sP_s(\mathbf{x})+\sum_{j=1}^{n}w_j\varphi(\|\mathbf{x}_j-\mathbf{x}\|),\qquad \varphi(r)=\begin{cases}r^k,&k\in\mathbb{N}\text{ odd},\\r^k\log r,&k\in\mathbb{N}\text{ even},\end{cases}$$

$$\begin{bmatrix}\mathbf{M}&\mathbf{B}\\\mathbf{B}^{T}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{w}\\\boldsymbol\beta\end{bmatrix}=\begin{bmatrix}\mathbf{y}\\\mathbf{0}\end{bmatrix},\qquad M_{ij}=\varphi(\|\mathbf{x}_i-\mathbf{x}_j\|),\quad B_{is}=P_s(\mathbf{x}_i),\quad B^T\mathbf{w}=0$$

$$\mathbf{w}=Q_{\perp}\mathbf{z},\qquad (Q_{\perp}^{T}\mathbf{M}Q_{\perp})\mathbf{z}=Q_{\perp}^{T}\mathbf{y},\qquad Q_{\perp}^{T}Q_{\perp}=I,\quad \operatorname{range}(Q_{\perp})=\ker(B^T)$$

$$P_{\mathrm{exact}}=Q_{\perp}^{T}\Sigma^{-1}Q_{\perp},\qquad P_{\mathrm{Vecchia}}=Q_{\perp}^{T}\tilde{\mathbf{T}}\tilde{\mathbf{T}}^{T}Q_{\perp}$$

Mathematical statement

The paper represents an interpolant by $\hat f(\mathbf{x})=\sum_{s=1}^{S}\beta_sP_s(\mathbf{x})+\sum_{j=1}^{n}w_j\varphi(\|\mathbf{x}_j-\mathbf{x}\|)$, where $\mathbf{x}_j\in\mathbb{R}^d$ are anchor locations, $P_s$ are low-degree polynomial basis functions, $\beta_s$ are trend coefficients, and $w_j$ are kernel coefficients. The polyharmonic kernel is $\varphi(r)=r^k$ for odd $k\in\mathbb{N}$, and $\varphi(r)=r^k\log r$ for even $k\in\mathbb{N}$. With $M_{ij}=\varphi(\|\mathbf{x}_i-\mathbf{x}_j\|)$, $B_{is}=P_s(\mathbf{x}_i)$, and observed values $\mathbf{y}$, the coefficients satisfy the augmented system $\begin{bmatrix}M&B\\B^T&0\end{bmatrix}\begin{bmatrix}\mathbf{w}\\\boldsymbol\beta\end{bmatrix}=\begin{bmatrix}\mathbf{y}\\0\end{bmatrix}$, equivalently $B^T\mathbf{w}=0$. Let $Q_\perp$ have orthonormal columns spanning $\ker(B^T)$; then $\mathbf{w}=Q_\perp\mathbf{z}$ and the reduced system is $(Q_\perp^TMQ_\perp)\mathbf{z}=Q_\perp^T\mathbf{y}$. The paper accelerates multiplication by $M$ using FMM implementations for logarithmic and distance kernels. It preconditions the projected system using $P_{\mathrm{exact}}=Q_\perp^T\Sigma^{-1}Q_\perp$, approximated by $P_{\mathrm{Vecchia}}=Q_\perp^T\tilde T\tilde T^TQ_\perp$, where $\Sigma$ is a Matérn covariance matrix and $\tilde T$ is a sparse reverse-Cholesky/precision factor built from a fixed-size Vecchia conditioning set. The transferable property is that dense geometric interactions are retained while only matrix-vector products and sparse triangular operations are required.

Implementation notes

1. Integration point: implement the head after a coordinate encoder or directly after normalized spatial coordinates. Store $n$ anchor locations $X=[\mathbf{x}_1,\ldots,\mathbf{x}_n]$ and one or more coefficient vectors $w^{(c)}$, one per output channel. For a query batch $X_q$, output the polynomial trend plus the PHS interaction $K(X_q,X)W$; do not materialize the $|X_q|\times n$ kernel matrix. Use $k=2$ for the thin-plate kernel in two dimensions, or compare $k=1,2,3,4$ in the ablation. 2. Pseudocode: construct $B_{is}=P_s(x_i)$, compute or implicitly apply a null-space basis $Q_\perp$, initialize $z$; repeat PCG: r = Q.T @ (y - M_fmm(Q @ z)); p = preconditioner(r); update z += alpha*p, r -= alpha*Q.T@M_fmm(Q@p); finally w=Q@z and solve $B^T\beta=B^T(y-Mw)$. In the forward pass compute yq = Bq@beta + fmm_phs(Xq, X, w). 3. Mathematical computations: FMM supplies the approximation to $Mv$ and $K(X_q,X)v$; the polynomial constraint is exact. Build the Vecchia preconditioner by ordering anchors, selecting $m$ nearest earlier anchors for each point, fitting local Matérn conditional regressions, and assembling a sparse factor $\tilde T$ with at most $m$ off-diagonal entries per column; estimate the Matérn length scale from anchor spacing or tune it on a small validation subset. Estimate PCG quality using residual norm and iteration count, not by forming eigenvalues. 4. First experiment: use a 2D SIREN or 4-layer ReLU coordinate MLP on the synthetic function in the paper, with 10k, 100k, and 1M Sobol anchors and dense held-out grid queries. Compare (a) dense RBF/PHS solve, (b) Adam-trained learned RBF coefficients, (c) FMM-PHS with no preconditioner, and (d) FMM-PHS with Vecchia-PCG. Measure wall-clock time, peak memory, relative interpolation error, FMM tolerance, and PCG iterations. Success means subquadratic memory, fewer than 20 PCG iterations, and near-dense-solve accuracy; a strong result is at least 5x lower runtime and 10x lower peak memory at 100k anchors, with no increase in held-out error.

Verification

Beats tuned baseline

Stage 1 · Toy mechanism gate: Passed ✓

Stage 2 · Mechanism transferred to benchmark: Confirmed ✓

Stage 2 · Practical benchmark result: Significant win ✓

Methodology: 8 paired seeds, permutation test (p<0.05) against a learning-rate-tuned baseline. How verification works

Stage 1 — Mechanism check agent confidence 8/10

Built a reproducible PHS interpolation head with exact polynomial null-space projection, projected PCG, a Vecchia-style Matérn preconditioner, and chunked query evaluation. The core predictions were confirmed: projected k=2 PHS interactions scaled as rho^2 with maximum error 3.6e-15, projected squared-distance polynomial terms vanished at 3.1e-15 relative magnitude, and computed weights satisfied the polynomial constraint near 1e-14. Vecchia preconditioning reduced PCG iterations from 292 to 196, but did not reach the paper-scale fewer-than-20-iteration target; this is a promising small-scale signal rather than a demonstrated FMM speedup.

Agent confidence: 8/10

Baseline: Unpreconditioned projected PCG: 292 iterations, relative residual 7.65e-08; dense query evaluation: 0.00288 seconds.
Idea: Vecchia-preconditioned projected PCG: 196 iterations with m=16 or 32, relative residual 3.83e-08 to 6.32e-08; exact-vs-chunked query error 0.0. Scaling sweep observed [0.25, 1, 4, 16] versus predicted rho^2 [0.25, 1, 4, 16].

Limitations: This MVP uses dense scipy kernel construction and multiplication rather than a true FMM, tests only n=180 in 2D and one synthetic function, and implements a simplified Vecchia precision rather than the paper's full reverse-Cholesky construction. It does not test 10k-1M anchors, GPU execution, peak-memory scaling, neural coordinate encoders, or a rigorous wall-clock FMM comparison.

How to run: OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 /home/maxwelhelp/main/bin/python3 phs_experiment.py

Stage 2 — Benchmark vs. tuned baseline (arXiv:2608.11462) Sep 2, 2026

Beats the tuned baseline by 48.6% (p=0.0081, wins 8 / 8 paired seeds; metric: lower is better).

Benchmark: Spatial_anchor_regression
Model: coordinate_mlp
Paired seeds: 8
Baseline mean: 0.0127
Idea mean: 0.0065
Effect (Δ): -0.0062 (−48.6%; negative = idea better)
Wins: 8 / 8 paired seeds
p-value: 0.0081 (permutation test, 20 000 shuffles)
Smallest detectable effect: ±22.2%
Mechanism: Confirmed ✓
Practical effect: Helps
Baseline tuning: swept over 3 configs

Limitations:

This stage-2 MVP uses a local 400-sample custom spatial regression track rather than the absent built-in spatial/PDE track. It uses 96 anchors and a dense torch.cdist interaction per batch, not a true FMM; it does not test 10k-1M anchors, wall-clock or peak-memory scaling, projected null-space coefficient constraints, Vecchia-PCG recalibration, coordinate-encoder integration beyond the small MLP, or the built-in tracks. The custom track is registered in-memory by the driver and is not written into the read-only shared bench directory.

How to run: OMP_NUM_THREADS=1 /home/maxwelhelp/main/bin/python3 run_bench.py

Verdict computed by deterministic test code from paired-seed statistics — not by the language model.

Stage 2 — Benchmark Sep 1, 2026 latest Worked ✓

Agent confidence: 8/10

Built a structurally matched custom spatial neural-field track and trained two end-to-end coordinate-MLP systems differing only by the added chunked k=2 PHS anchor head. Across 8 paired seeds, the idea reduced MSE from 0.012713 to 0.006530, delta=-0.006183, with all 8 wins and permutation p=0.0081; the benchmark therefore reports a significant win. The trained PHS residual was nonzero and aligned with held-out baseline error correction (mean correlation 0.704), confirming the tested mechanism at NN scale.

Baseline: Baseline coordinate MLP, lr=0.01 selected by sweep over {0.001, 0.003, 0.01}: mean test MSE 0.0127128037, std 0.0031444912, per-seed [0.0102508850, 0.0145432632, 0.0139772529, 0.0109899705, 0.0170229226, 0.0169589557, 0.0095467241, 0.0084124552].
Delta mean: -0.0061826404
Idea: Coordinate MLP plus chunked differentiable k=2 PHS anchor head, lr=0.01 selected from parity sweep: mean test MSE 0.0065301632, std 0.0013944949, per-seed [0.0083841076, 0.0057844375, 0.0060189855, 0.0073700463, 0.0054297075, 0.0089141643, 0.0052020513, 0.0051378058].
Idea wins: 8
Mechanism signature: map[confirmed:true observed_mean_correction_correlation:0.7042187005 observed_mean_residual_rms:0.1142530991 prediction:The trained PHS residual should be nonzero and align with held-out error correction.]
Permutation p value: 0.0081
Relative change percent: -48.64

Benchmark result