A Data-Driven Interpolation Method on Smooth Manifolds via Diffusion Processes and Voronoi Tessellations

1. Problem Statement

The paper addresses the challenge of approximating smooth, real-valued functions $f: M \to \mathbb{R}$ defined on a smooth manifold $M \subset \mathbb{R}^d$, given only a finite set of pointwise evaluations $\{x_i, f(x_i)\}_{i=1}^N$.

Key Challenges:

• Intrinsic Geometry: The data lies on a low-dimensional manifold embedded in a high-dimensional space, requiring methods that respect the manifold's intrinsic geometry rather than just Euclidean distance.
• Computational Efficiency: Classical data-driven methods (e.g., Gaussian Process Regression, Radial Basis Functions, Neural Networks) often suffer from high computational costs. GPR scales as $O(N^3)$, and RBFs as $O(N^2)$, while Neural Networks require iterative training phases that can be unstable and time-consuming.
• Sparse Data & Noise: In applications like sparse-view Computed Tomography (CT), data is limited and noisy, leading to ill-posed reconstruction problems prone to spurious oscillations.

2. Methodology

The proposed approach is a fully data-driven, non-iterative interpolation method that combines Diffusion Maps concepts with Voronoi Tessellations.

Core Components:

1. Gaussian Kernel Approximation:
   The method utilizes a rescaled Gaussian kernel $G_\varepsilon(z) = \exp(-\|z\|^2/\varepsilon^2)$. Based on the Diffusion Maps framework, the approximation of a function $g$ at a point $x$ is defined via a normalized convolution:
   $$\bar{g}(x) = \frac{\sum_{i=1}^N G_{\varepsilon(x)}(x - x_i) g(x_i)}{\sum_{i=1}^N G_{\varepsilon(x)}(x - x_i)}$$

2. Adaptive Scale Parameter ($\varepsilon(x)$):
   A critical innovation is the selection of the scale parameter $\varepsilon$. Instead of a global constant, the method defines $\varepsilon(x)$ adaptively for each query point $x$ based on the Voronoi tessellation of the dataset:
   $$\varepsilon(x) = \min_{1 \le i \le N} \|x - x_i\|$$

   • Significance: This choice ensures the normalization factor in the denominator is uniformly bounded away from zero (specifically $\ge e^{-1}$), preventing numerical instability. It also allows the method to adapt to local data density.

3. Algorithm:

   • Input: Sample set $X$, function values, and a query point $x$.
   • Step 1: Compute $\varepsilon(x)$ as the distance to the nearest neighbor.
   • Step 2: Compute the weighted sum of Gaussian kernels using $\varepsilon(x)$.
   • Output: The interpolated value $\bar{g}(x)$.
   • Complexity: The inference step scales linearly with the number of samples $N$ ($O(N)$), as it requires no training phase.

3. Key Theoretical Contributions

The paper establishes several rigorous mathematical properties of the proposed interpolant:

• Vanishing Gradient at Samples: The gradient of the interpolant $\nabla \bar{g}(x)$ vanishes exactly at the sample points $x_i$ and decays exponentially in their neighborhoods.
• Total Variation (TV) Minimization: The method is shown to be an approximate solution to a Total Variation minimization problem. Specifically, it minimizes an energy functional of the form $\|f - g\| + \lambda \|\nabla f\|$. In the case where the forward operator is the identity, it provides a closed-form analytical solution to a compressed sensing problem.
• Spectral Attenuation (Low-Pass Filtering): The method acts as a data-driven low-pass filter. By penalizing the Dirichlet energy (gradient norm), it suppresses high-frequency components of the signal. As the number of samples $N \to \infty$, the gradient energy of the interpolant converges to zero with high probability, effectively attenuating noise and high-frequency oscillations.
• No Training Required: Unlike Neural Networks or iterative optimization methods, this approach requires no backpropagation or hyperparameter tuning for training; the approximation is derived directly from the data geometry.

4. Experimental Results

The authors validated the method through two primary experiments:

A. Synthetic Interpolation on a Torus

• Setup: Interpolation of a high-frequency function $f(x,y) = \sin(10\pi x)\cos(10\pi y)$ on a 2D domain.
• Comparison: Compared against Feedforward Neural Networks (FNN), Gaussian Process Regression (GPR), and Radial Basis Function Networks (RBFN).
• Results:
  • Accuracy: The proposed method achieved the lowest Mean Squared Error (MSE) across all sample sizes (100 to 2000 points).
  • Efficiency: It was significantly faster. For 2000 points, the proposed method took 0.53 seconds, whereas RBFN took 1436 seconds and GPR took 13.49 seconds. FNN required training time not included in the inference-only comparison but still lagged in total pipeline efficiency.

B. Sparse-View CT Reconstruction

• Setup: Reconstruction of Shepp-Logan phantoms, human heads, and human abdomens from sparse angular projections (100 and 300 angles vs. the standard ~720).
• Method: The proposed method was used to interpolate the sinogram (angular data) before applying Filtered Back Projection (FBP). It was also used as an initialization for iterative Total Variation (TV) regularization.
• Results:
  • Reconstruction Quality: The proposed method produced the most accurate reconstructions (lowest MSE) compared to cubic splines, GPR, FNN, and RBFN.
  • TV Regularization Synergy: When used as an initialization for iterative TV minimization, the proposed method often reached the optimal solution faster than standard FBP initialization. In some high-sampling cases, further TV iterations actually increased error when starting from the proposed method, suggesting the method already provides a near-optimal, regularized solution.
  • Computational Cost: The method drastically reduced the time required to achieve high-quality regularized reconstructions compared to iterative baselines.

5. Significance and Impact

• Scalability: The linear scaling ($O(N)$) and lack of a training phase make this method highly scalable for large datasets where traditional kernel methods (GPR) become computationally prohibitive.
• Stability: The adaptive Voronoi-based scale parameter ensures numerical stability, a common failure point in kernel methods with small bandwidths.
• Theoretical Bridge: The work bridges Manifold Learning, Compressed Sensing, and Total Variation Regularization, providing a closed-form solution that implicitly solves complex optimization problems.
• Practical Application: It offers a powerful preprocessing tool for medical imaging (CT), enabling high-quality reconstruction from sparse data without the computational burden of iterative optimization, potentially reducing radiation exposure and scan times.

In conclusion, the paper presents a robust, efficient, and theoretically grounded interpolation technique that outperforms state-of-the-art data-driven methods in both accuracy and computational speed, particularly in scenarios involving sparse data and geometric constraints.