VKSR: Scalable Kernel Surface Reconstruction Using Vecchia's Approximation¶

Conference: ECCV 2026
Paper: ECCV Official
Code: https://mweiherer.github.io/vksr/
Area: 3D Vision
Keywords: Surface Reconstruction, Implicit Surfaces, Kernel Methods, Gaussian Processes, Vecchia Approximation

TL;DR¶

Addressing the excessive computational cost and oversmoothing of traditional kernel surface reconstruction on dense or large-scale point clouds caused by global low-rank approximations, VKSR introduces Vecchia's approximation from Gaussian Processes to shift computation from a global solve to fully localized response-process inference, reconstructing 14M+ points in minutes with up to 180× speedup over existing kernel baselines.

Background & Motivation¶

Recovering a continuous implicit surface (e.g., a Signed Distance Function, SDF) from an oriented point cloud is a foundational problem in 3D computer vision and computer graphics. In recent years, kernel surface reconstruction methods such as Neural Splines and Matérn kernels have emerged as a compelling paradigm. By formalizing implicit surface reconstruction as a Kernel Ridge Regression (KRR) problem, they offer closed-form solutions and enable easy injection of geometric inductive biases (such as smoothness priors) through kernel selection. Consequently, they achieve superior robustness over classical radial basis functions (RBFs), widely used Screened Poisson Surface Reconstruction (SPSR), and coordinate-based neural fields when reconstructing sparse, incomplete, or noisy inputs.

However, standard kernel ridge regression requires solving an $m \times m$ linear system with cubic time $O(m^3)$ and quadratic space $O(m^2)$ complexity in the number of input points $m$, rendering it computationally intractable for dense or large-scale point clouds containing millions of points. To address this limitation, existing methods rely on low-rank approximations of the kernel matrix using $n \ll m$ Nyström samples. This design implicitly assumes that the entire underlying surface can be represented by a small set of data-dependent global basis functions. This global smoothness assumption severely degrades high-frequency geometric details in dense inputs. To recover fine features, practitioners must dramatically increase the number of Nyström samples—causing runtimes to explode into hours—or fall back to spatial subdivision (chunking and blending), which is prone to boundary artifacts and fragile to chunk size selection.

To achieve scalable, artifact-free reconstruction without global low-rank assumptions or spatial partitioning, the authors exploit the theoretical equivalence between the KRR closed-form solution and the posterior predictive mean of a zero-mean Gaussian Process (GP). Inspiringly, they introduce Vecchia's approximation—a powerful technique from spatial statistics and GP literature—to implicit surface reconstruction. Core idea: replace global low-rank kernel approximations with Vecchia's conditional independence approximation formulated directly on the response process, turning an intractable global solve into embarrassingly parallel, highly conditioned local neighborhood linear systems to achieve minute-level high-fidelity reconstruction of 14M+ points on a single GPU.

Method¶

Overall Architecture¶

VKSR accepts an oriented point cloud $X_0 \subset \mathbb{R}^3$ with per-point normals, offsets points by $\pm \epsilon$ along normals to construct an augmented point cloud $\mathcal{X}$ of $m = 3m_0$ points with SDF target values $y$. Unlike existing approaches, VKSR does not compute or invert any global $m \times m$ or $n \times n$ matrix. Instead, the continuous surface is reconstructed on demand across spatial evaluation queries. For any spatial query point $x$, VKSR first retrieves its $k$-nearest neighbors $N_k(x) \subset \mathcal{X}$ using an accelerated partition-based index. It then forms an independent, well-conditioned $k \times k$ local kernel matrix to analytically solve for the query-dependent weight vector via a direct single-precision Cholesky solver. Finally, it evaluates the local kernel sum to compute the SDF value $\hat{f}^{(k)}(x)$ and extracts the final triangular mesh across an adaptive sparse multiresolution octree grid via Marching Cubes.

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    A["Input Oriented Point Cloud<br/>Generate Augmented Points with Offsets"] --> B["Response-Process Vecchia Approximation<br/>Decouple Global Matrix to Local $k$-NN Inference"]
    B --> C["Partition-Based Nearest Neighbor Search<br/>Lloyd Quantization & Multi-Probing"]
    C --> D["Local Kernel Solve<br/>Batched FP32 Cholesky Inversion"]
    D --> E["Sparse Multiresolution Surface Extraction<br/>Marching Cubes on Octree Voxel Grid"]

Key Designs¶

1. Response-Process Vecchia Approximation: Shifting from Global Low-Rank to Local Neighborhoods Existing kernel methods rely on Nyström low-rank matrix approximation $\tilde{K} = K_{\mathcal{X}\mathcal{Z}}K_{\mathcal{Z}}^{-1}K_{\mathcal{X}\mathcal{Z}}^\top$, which imposes a global smoothness prior that blurs sharp surface details and demands substantial Nyström samples. VKSR leverages the insight that the optimal KRR solution matches the posterior predictive mean of a zero-mean Gaussian Process $f \sim \mathcal{GP}(0, k)$ with observational noise variance $\sigma^2 = \lambda$. Crucially, the authors apply Vecchia's approximation directly to the noisy response process $y \sim \mathcal{N}(0, K_{\mathcal{X}} + \sigma^2 I)$ rather than the latent process. Through the chain rule of probabilities, Vecchia assumes each observation is conditionally independent given its $k$-nearest spatial neighbors, inducing a sparse precision matrix. For any spatial evaluation point $x \in \mathbb{R}^3$, the posterior predictive mean simplifies to a purely local formulation: $$ \hat{f}^{(k)}(x) = k_{N_k(x)}(x)^\top (K_{N_k(x)} + \lambda I)^{-1} y_{N_k(x)} = \sum_{x_i \in N_k(x)} \hat{\alpha}_i(x) k(x, x_i) $$ where $N_k(x)$ contains the $k$-nearest neighbors of $x$ in $\mathcal{X}$, and $\hat{\alpha}(x) = (K_{N_k(x)} + \lambda I)^{-1} y_{N_k(x)} \in \mathbb{R}^k$ represents local kernel weights. This mathematically reduces the global $O(m^3)$ solve to decoupled, miniature $O(k^3)$ linear systems where $k \in [16, 128]$ is sufficient for millimeter precision. Crucially, this imposes a local smoothness assumption that aligns naturally with physical surfaces, eliminating global over-smoothing while strictly converging to standard KRR as $k \to m$.

2. Partition-Based Approximate Nearest Neighbor Search: Eliminating the $O(m)$ Query Bottleneck While localizing the kernel solve lowers matrix inversion to $O(k^3)$, naive linear nearest neighbor scanning still scales as $O(m)$ per query point, which becomes the dominant bottleneck when evaluating $t$ points over millions of input observations. To break this dependency, VKSR integrates a partition-based nearest neighbor search using multi-probing (powered by FAISS). During a one-time preprocessing step, the augmented dataset $\mathcal{X}$ is partitioned into $C \ll m$ compact clusters via Lloyd's algorithm in $O(Cmd)$ time ($d=3$). At inference time, for each query point $x$, the system locates the $P \ll C$ nearest centroids in $O(C(d+P))$ time and restricts the $k$-NN search strictly within these $P$ candidate clusters in $O((Pm/C)(d+k))$ time. Furthermore, data arrays remain in system RAM while query batches stream through the GPU, allowing millions of points to be indexed and processed with a minimal GPU memory footprint of a few hundred megabytes.

3. Sparse Multiresolution Voxel Evaluation and Batched FP32 Cholesky Solves: Throughput Optimization The overall runtime of VKSR is governed by the total number of query evaluations $t$. Evaluating a dense $1024^3$ voxel grid would require an excessive number of localized solves in empty space. Instead, VKSR computes $\hat{f}^{(k)}$ on a sparse voxel grid constructed via adaptive, conditional subdivision from an initial coarse grid, evaluating only active cells near the zero-level surface. In addition, because local $k \times k$ matrices are inherently well-conditioned compared to ill-conditioned global kernel systems, VKSR bypasses expensive double-precision conjugate gradient solvers. It uses single-precision (FP32) direct Cholesky decomposition batched over tens of thousands of queries concurrently on the GPU, maximizing hardware utilization.

Loss & Training¶

As a non-parametric kernel regression method, VKSR does not undergo iterative backpropagation. The closed-form analytical solution directly minimizes the local empirical regularized risk at query point $x$: $$ \min_{f \in \mathcal{H}k} \left{ \sum^2 \right} $$ Following established kernel surface reconstruction practices, the default covariance function is chosen from the Matérn family (e.g., Matérn } (f(x_i) - y_i)^2 + \lambda |f|_{\mathcal{H}_k$\nu=1/2$, corresponding to the Laplace kernel $k(x, x') = \exp(-\frac{\|x - x'\|_2}{\ell})$, or Matérn $\nu=3/2$), where $\ell$ is the kernel length scale. Because queries are solved independently, inference can be partitioned into arbitrary query batch chunks, guaranteeing deterministic execution within any prescribed VRAM budget (e.g., under 1 GB).

Key Experimental Results¶

Main Results¶

On ShapeNet (256 objects across 13 classes), VKSR is evaluated under dense ($m_0=100\text{K}$ points, $k=32$) and sparse ($m_0=1\text{K}$ points, $k=128$) settings. Metrics include F-Score (FS ↑, threshold 0.01), Chamfer Distance (CD ↓, scaled by $10^3$), Normal Consistency (NC ↑, %), runtime (seconds), and peak VRAM (MB).

Dataset / Setting	Method	FS (%) ↑	CD ($\times 10^{-3}$) ↓	NC (%) ↑	Time (s) ↓	VRAM (MB) ↓
ShapeNet Dense (100K pts)	RIMLS	99.5	2.33	97.9	115.11	N.A.
	SAP	97.4	3.64	95.0	303.75	456.2
	SPSR	99.9	2.04	96.5	0.82	N.A.
	Neural Splines (NS)	99.5	2.34	96.5	131.95	1703.0
	Matérn ($\nu=1/2$)	99.7	2.25	96.5	47.91	1703.0
	Matérn ($\nu=3/2$)	99.6	2.28	96.7	47.35	1703.0
	VKSR (Ours)	99.8	2.19	96.8	0.77	342.2
ShapeNet Sparse (1K pts)	RIMLS	78.4	9.27	87.2	0.83	N.A.
	SAP	68.4	13.47	73.7	87.84	163.7
	SPSR	91.5	4.36	88.6	0.45	N.A.
	Neural Splines (NS)	94.1	3.91	93.0	2.69	274.7
	Matérn ($\nu=1/2$)	93.5	3.90	92.7	0.66	274.7
	Matérn ($\nu=3/2$)	94.0	3.89	92.9	0.75	274.7
	VKSR (Ours)	93.6	3.95	92.7	2.35	522.0

On the Stanford 3D Scanning Repository (dense point clouds up to 14M points, extracted at $1024^3$ grid resolution), VKSR demonstrates orders-of-magnitude faster execution over prior kernel baselines:

Object (Points)	Metric	SPSR	NS ($n=15\text{K}$)	NS (Chunked)	Matérn 1/2 (Chnk.)	VKSR (Exact)	VKSR (Approx.)
Armadillo (173K)	CD ↓ / Time (s) ↓	0.10 / 7.7s	1.27 / 14521s	0.10 / 3999s	0.10 / 1000s	0.10 / 90.3s	0.10 / 57.2s
Asian Dragon (3.6M)	CD ↓ / Time (s) ↓	0.10 / 18.5s	0.48 / 7192s	0.10 / 9182s	0.10 / 3496s	0.10 / 562.9s	0.10 / 126.6s
Lucy (14M)	CD ↓ / Time (s) ↓	0.94 / 25.7s	3.13 / 7169s	9.13 / 11314s	4.06 / 4231s	0.81 / 1809.8s	0.80 / 80.2s
Thai Statue (5.1M)	CD ↓ / Time (s) ↓	0.25 / 24.2s	0.71 / 7749s	2.89 / 12636s	0.29 / 4519s	0.24 / 840.3s	0.24 / 61.5s

Ablation Study¶

On the Surface Reconstruction Benchmark (evaluating resilience against range scan noise and missing data at $256^3$ resolution, reporting CD, Hausdorff Distance HD, and runtime):

Object	Metric	RIMLS	SAP	SPSR	Neural Splines	Matérn 3/2	VKSR (Ours)
Anchor (85K)	CD ↓ / HD ↓ / Time (s)	0.27 / 7.66 / 34s	0.35 / 8.91 / 350s	0.30 / 7.22 / 2.7s	0.26 / 5.33 / 308s	0.26 / 19.98 / 75s	0.24 / 6.33 / 13.8s
Daratech (61K)	CD ↓ / HD ↓ / Time (s)	0.21 / 2.86 / 7.3s	0.23 / 3.08 / 303s	0.23 / 6.04 / 1.5s	0.22 / 4.71 / 231s	0.23 / 4.66 / 61s	0.23 / 6.22 / 5.7s
DC (71K)	CD ↓ / HD ↓ / Time (s)	0.17 / 2.96 / 27s	0.19 / 3.32 / 354s	0.16 / 2.65 / 2.1s	0.15 / 1.26 / 313s	0.15 / 1.34 / 71s	0.16 / 2.08 / 4.9s
Gargoyle (95K)	CD ↓ / HD ↓ / Time (s)	0.18 / 4.07 / 38s	0.19 / 5.71 / 395s	0.18 / 4.46 / 2.5s	0.18 / 3.24 / 461s	0.17 / 3.13 / 105s	0.17 / 2.42 / 15.1s
Lord Quas (57K)	CD ↓ / HD ↓ / Time (s)	0.13 / 2.09 / 14s	0.14 / 3.36 / 302s	0.13 / 1.51 / 1.4s	0.12 / 0.96 / 207s	0.12 / 1.10 / 58s	0.13 / 2.79 / 24.9s

Ablation investigations into $k$ and query chunking reveal: 1. Neighborhood Size $k$: Denser point clouds achieve excellent reconstruction with smaller neighborhoods (e.g., $k=16$ or $32$ yields CD of 2.19 on ShapeNet Dense), whereas sparse inputs require slightly larger contexts ($k=128$) to maintain continuity. 2. Deterministic VRAM Cap: Because individual query evaluations are independent, adjusting query chunk sizes from 1K to 10K allows VKSR to strictly constrain its VRAM footprint to 1GB–2GB across massive datasets without risking out-of-memory errors.

Key Findings¶

Eliminating the Low-Rank and Chunking Bottleneck: On the 14M point Lucy model, global Nyström approximation takes hours and washes out fine geometric structures, while space-chunking introduces severe seam artifacts (CD degrades to 4.06–9.13). VKSR reconstructs Lucy in just 80.2 seconds with a state-of-the-art CD of 0.80, delivering over 50× to 140× speedups over chunked kernel methods.
The Ideal Compromise between SPSR and Kernel Methods: SPSR is fast on dense clouds but fails on sparse or incomplete inputs (CD degrades to 4.36 on ShapeNet Sparse). Global kernel methods shine on sparse data but freeze on dense data. VKSR bridges this divide: it matches SPSR in dense runtime (0.77s vs 0.82s) while preserving kernel methods' inductive bias on sparse data (CD 3.95 vs SPSR's 4.36).
Near-Lossless Approximate $k$-NN Acceleration: Multi-probing FAISS approximate nearest neighbor search yields a ~9× speedup over exact nearest neighbor search on large point clouds with negligible accuracy difference (e.g., Lucy runtime drops from 1809.8s to 80.2s while CD remains 0.80).

Highlights & Insights¶

Formulating Vecchia on the Response Process: Applying Vecchia's approximation to the noisy response process rather than the noise-free latent process is the pivotal algorithmic breakthrough, decoupling the posterior mean into fully local calculations and unlocking massive GPU parallelism.
Local Smoothness Aligns with Physical Reality: Real-world 3D geometries exhibit local spatial continuity rather than long-range covariance coupling. Localized kernels naturally preserve sharp local features, bypass ill-conditioned global matrix inversions, and allow stable single-precision Cholesky factorization.
Bridging Classical Non-Parametric Theory and Modern Hardware: VKSR seamlessly integrates analytical Gaussian Process theory with streaming RAM indices and GPU batched tensor operations, proving that classical closed-form kernel methods can rival modern deep learning models in scalability and throughput.

Limitations & Future Work¶

Performance in Ultra-Sparse Regions: On extremely sparse or irregularly sampled point clouds, fixed Euclidean $k$-NN neighborhoods may bridge disconnected topological surfaces. Implementing hierarchical multi-scale neighbor search represents a promising refinement.
Unexploited Uncertainty Quantification: While VKSR natively possesses a closed-form posterior predictive variance formulation from GP theory, the current work focuses exclusively on mean SDF extraction. Leveraging predicted variance for active robotic scanning and grasping remains an exciting future direction.
Integration with Learnable Neural Kernel Fields: Preliminary investigations in the supplementary material indicate that VKSR can be embedded as a differentiable layer within learnable frameworks (e.g., Neural Kernel Fields), combining learned priors with non-parametric surface guarantees.

vs Neural Splines & Matérn Global Kernels: Prior kernel surface reconstruction methods rely on Nyström low-rank matrix approximations or spatial chunking, suffering cubic/quadratic complexity explosions and boundary seams; VKSR replaces global low-rank solving with local Vecchia inference, achieving up to 180× speedups without chunking artifacts.
vs Screened Poisson Surface Reconstruction (SPSR): SPSR solves an octree Poisson system quickly on dense scans but deteriorates on sparse or partial inputs; VKSR retains continuous kernel interpolation priors, dramatically outperforming SPSR on sparse/noisy data while matching its dense runtime on GPU.
vs Implicit Moving Least Squares (IMLS / RIMLS): IMLS performs moving local regressions using constant or low-degree polynomials, often failing to close holes in incomplete scans; VKSR performs non-parametric kernel ridge regression in an RKHS, with theoretical proofs demonstrating that classical IMLS is a restricted special case of VKSR.

Rating¶

Novelty: ⭐⭐⭐⭐⭐ Pioneering introduction of Vecchia's response-process approximation into implicit surface reconstruction, resolving the long-standing scalability bottleneck of kernel methods.
Experimental Thoroughness: ⭐⭐⭐⭐⭐ Rigorous validation across ShapeNet, Stanford 3D, ScanNet, and noisy benchmark scans spanning 1K to 14M points with extensive ablation of $k$, chunk size, and precision.
Writing Quality: ⭐⭐⭐⭐⭐ Elegant mathematical derivations, thorough asymptotic complexity comparisons, and crystal-clear presentation of hardware engineering optimizations.
Value: ⭐⭐⭐⭐⭐ Establishes a fast, analytically elegant, and highly scalable baseline for large-scale 3D surface reconstruction.