Interpolation-Regression Approximation

• Interpolation-Regression Approximation is a continuum of methods that unifies exact interpolation with regression through RKHS and kernel-based formulations.
• It involves constrained, penalized, and reconstruction approaches which balance pointwise fidelity with regularization to enhance stability and accuracy.
• The framework applies to Gaussian process regression, polynomial reconstruction, and geometric approximations on manifolds, demonstrating both theoretical and practical benefits.

Across the cited works, interpolation-regression approximation appears not as a single algorithm but as a family of constructions that place exact interpolation, approximate interpolation, least-squares fitting, and regularization on a common continuum. In the most explicit RKHS formulation, exact kernel interpolation is the minimum-native-space-norm extension satisfying pointwise constraints, while regression replaces the hard constraint $u(\mathbb X)=\mathbb Y$ by a quadratic fidelity term and recovers interpolation in the limit $\alpha\searrow 0$ (Guidotti, 6 Feb 2026). In nonparametric regression, the same bridge is realized by representing the unknown function through an interpolator evaluated at finitely many knots and estimating only the knot values by least squares or penalized least squares (Xiong, 2018). A related Mercer-kernel viewpoint shows that Gaussian process regression, kernel interpolation, and pseudospectral polynomial approximation can coincide under specific choices of kernel, noise level, and design (Gorodetsky et al., 2015).

1. Conceptual landscape

A convenient way to organize the subject is to distinguish three recurring mechanisms. First, there are constrained formulations, where one interpolates exactly on a selected subset and fits the remaining data in a least-squares sense. Second, there are penalized formulations, where interpolation is replaced by a soft data-fidelity term weighted by a regularization parameter. Third, there are reconstruction formulations, where interpolation is used as an internal parametrization of a regression problem rather than as the final estimator itself. In a broader but explicitly qualified sense, some works also use interpolation as a deterministic surrogate inside a larger learning or optimization pipeline, even when no statistical regression objective is present (Guidotti, 6 Feb 2026, Xiong, 2018, Daulbaev et al., 2020).

Paradigm	Core construction	Representative sources
Constrained interpolation-regression	Exact interpolation on a subset, least-squares fit on the full sample or remainder	(Marchi et al., 2014, Dell'Accio et al., 9 Aug 2025, Dell'Accio et al., 13 May 2026)
Penalized kernel approximation	$\|u\|_{\mathcal H_K}^2$ regularization plus quadratic data mismatch	(Guidotti, 6 Feb 2026)
Reconstruction-based regression	Estimate knot values, then reconstruct by interpolation	(Xiong, 2018)
Internal interpolatory surrogate	Sample $\to$ interpolate $\to$ substitute in downstream equations	(Daulbaev et al., 2020)

A common misconception is that interpolation and regression are disjoint categories. The cited literature repeatedly treats them as adjacent regimes of the same approximation problem. Exact interpolation may arise as the zero-regularization or zero-noise endpoint of regression, while regression may be understood as approximate interpolation with controlled residuals. Conversely, interpolation can also serve as a structural map inside a regression estimator, so that the unknown function is reconstructed from a lower-dimensional parameter vector of knot values rather than fitted directly in an ambient infinite-dimensional space (Guidotti, 6 Feb 2026, Xiong, 2018).

2. RKHS, kernel, and Gaussian-process formulations

The cleanest unification is variational. For a positive definite kernel $K$ with native space $\mathcal H_K$, the exact interpolant at data $(\mathbb X,\mathbb Y)$ is

$$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$

and it is characterized as the minimizer of $\frac12\|u\|_{\mathcal H_K}^2$ under the hard constraint $\alpha\searrow 0$0. Regression or approximate interpolation replaces the constraint by

$\alpha\searrow 0$1

which yields

$\alpha\searrow 0$2

with strong convergence $\alpha\searrow 0$3 in $\alpha\searrow 0$4 as $\alpha\searrow 0$5. The same paper interprets exact interpolation as a Fredholm integral equation of the first kind and approximate interpolation as a second-kind regularization, thereby making the stabilization role of regression explicit (Guidotti, 6 Feb 2026).

The Gaussian-process interpretation is formally identical. With prior $\alpha\searrow 0$6 and observation noise variance $\alpha\searrow 0$7, the posterior mean is

$\alpha\searrow 0$8

Setting $\alpha\searrow 0$9 gives exact interpolation, while $\|u\|_{\mathcal H_K}^2$0 reproduces the RKHS regression formula. A Mercer-kernel expansion further shows that GP regression can coincide with pseudospectral polynomial approximation when the kernel eigenfunctions match the polynomial basis, the design orthogonalizes those eigenfunctions, and the nugget vanishes; departures from that coincidence are attributable to nugget regularization, kernel truncation, or loss of discrete orthogonality (Gorodetsky et al., 2015).

Scalable GP work preserves the interpolation viewpoint but alters the interpolation geometry. Soft Kernel Interpolation replaces SKI’s dense rectilinear lattice by a learned set of interpolation points $\|u\|_{\mathcal H_K}^2$1 and a softmax interpolation matrix

$\|u\|_{\mathcal H_K}^2$2

so that

$\|u\|_{\mathcal H_K}^2$3

The method has $\|u\|_{\mathcal H_K}^2$4 time and $\|u\|_{\mathcal H_K}^2$5 space, and is designed to remove the lattice-induced $\|u\|_{\mathcal H_K}^2$6 bottleneck of SKI while retaining an interpolation-based kernel approximation (Camaño et al., 2024).

3. Finite-dimensional polynomial and reconstruction methods

In nonparametric regression, the reconstruction approach makes interpolation the parametrization layer. One chooses knots $\|u\|_{\mathcal H_K}^2$7, an interpolator $\|u\|_{\mathcal H_K}^2$8 with $\|u\|_{\mathcal H_K}^2$9, and estimates the knot-value vector $\to$0 by

$\to$1

If the interpolation error is

$\to$2

and the knot-value estimator has mean-square accuracy $\to$3, then the pointwise MSE is $\to$4. This supplies the basic interpolation-regression decomposition: interpolation bias plus statistical estimation error. The same framework reinterprets polynomial regression and kernel ridge regression as special cases, and yields Gaussian process reconstruction regression with complexity $\to$5 instead of the $\to$6 cost of full GP/KRR when $\to$7 (Xiong, 2018).

Polynomial kernels furnish a finite-feature version of the same bridge. For

$\to$8

the feature map is finite-dimensional, the kernel matrix factorizes as $\to$9, and invertibility is equivalent to full row rank of the Vandermonde-type matrix $\to$0. When $\to$1 and the nodes are polynomially unisolvent, polynomial-kernel interpolation collapses exactly to ordinary polynomial interpolation in $\to$2 for $\to$3 or $\to$4 for $\to$5. The native spaces are explicit finite-dimensional polynomial spaces, and the paper notes that increasing $\to$6 may connect to the overparameterized regime and ridgeless regression, with $\to$7 acting as an implicit regularizer through the RKHS norm (Elefante et al., 2022).

On equispaced intervals, constrained mock-Chebyshev least squares is a canonical simultaneous interpolation-regression construction. One interpolates on a mock-Chebyshev subset $\to$8, writes every constrained polynomial as

$\to$9

and chooses $K$0 by weighted least squares on the remaining nodes. The final approximant

$K$1

is strictly better than pure mock-Chebyshev interpolation in the discrete $K$2-norm on the leftover set and, under a sufficient condition, improves it in the uniform norm as well. The method is designed precisely to use all equispaced samples without reverting to unstable full equispaced interpolation and the Runge phenomenon (Marchi et al., 2014).

4. Geometric and domain-specific constructions

On manifolds and point clouds, interpolation-regression approximation acquires a geometric role. Kernel interpolation and regression are used to build approximate level-set or defining functions for sampled hypersurfaces,

$K$3

from which normals, tangent planes, and principal curvatures can be computed analytically. The paper’s central message is that exact interpolation is the zero-regularization limit of regression, but for smooth kernels it may be extremely ill-posed, whereas regression offers a natural regularization for noisy or imperfect geometric data (Guidotti, 6 Feb 2026).

On the unit disk, the same philosophy is implemented with Zernike polynomials and a Bos-array-inspired interpolation subset. The approximation space is $K$4, represented in practice through Zernike polynomials, while the interpolation subset $K$5 is chosen as a mock-optimal approximation to the optimal Bos array $K$6. The coefficient vector $K$7 is obtained from the constrained least-squares problem

$K$8

and the resulting operator is a projection onto $K$9. The same construction is then used to derive cubature formulas by integrating the reconstructed polynomial over the disk (Dell'Accio et al., 9 Aug 2025).

On $\mathcal H_K$0, the hybrid object is a spherical polynomial of degree $\mathcal H_K$1 that interpolates exactly on a prescribed subset $\mathcal H_K$2 and minimizes the discrete residual over the full sample set. In coefficient form,

$\mathcal H_K$3

with solution characterized by a KKT system and by discrete orthogonality of the residual against the feasible variation space. In antipodally symmetric settings, the problem splits into independent even and odd components in a real spherical-harmonic basis, while for spherical designs the normal matrix becomes scalar and the spectral condition number of the KKT matrix can be written explicitly (Dell'Accio et al., 13 May 2026).

5. Approximation theory, conditioning, and basis choice

The literature repeatedly separates approximation quality from numerical stability. Exact solvability does not imply stable computation, and the decisive objects are often the node geometry, the basis, and the normalization. A rescaled RBF method makes this explicit by defining

$\mathcal H_K$4

where $\mathcal H_K$5 is the standard kernel interpolant of $\mathcal H_K$6 and $\mathcal H_K$7 is the interpolant of the constant function $\mathcal H_K$8. This quotient form is equivalent to interpolation with the rescaled kernel

$\mathcal H_K$9

and also to a Shepard-type approximation with normalized cardinal functions

$(\mathbb X,\mathbb Y)$0

Hence the method remains interpolatory, reproduces constants exactly, and forms a partition of unity, yet often reduces oscillations and improves empirical stability (Marchi et al., 2016).

For polynomial kernels, stability issues arise even when interpolation is theoretically admissible. The kernel matrix may be singular unless the Vandermonde-type matrix has full row rank, and direct solution of $(\mathbb X,\mathbb Y)$1 can be highly ill-conditioned. The paper therefore adapts the RBF-QR methodology to build a numerically stable basis for the interpolation space $(\mathbb X,\mathbb Y)$2, rather than working in the unstable translate basis $(\mathbb X,\mathbb Y)$3. The numerical experiments show that direct kernel-matrix inversion may fail to converge even for modest $(\mathbb X,\mathbb Y)$4, while RBF-QR recovers stable convergence and accuracy comparable to polynomial interpolation (Elefante et al., 2022).

Exact interpolation also provides benchmark approximation rates. For Gaussian interpolation at nonuniform samples $(\mathbb X,\mathbb Y)$5 generated by a Riesz-basis sequence, one has

$(\mathbb X,\mathbb Y)$6

together with derivative estimates $(\mathbb X,\mathbb Y)$7. In a different direction, piecewise affine nodal interpolation on suitably chosen simplicial triangulations approximates Sobolev functions directly in mixed $(\mathbb X,\mathbb Y)$8-/$(\mathbb X,\mathbb Y)$9-Sobolev topologies without first passing through smooth-density arguments (Hamm, 2013, Schaftingen, 2013).

A related misconception is that regression becomes theoretically relevant only when noise is present. In convex-constrained derivative-free optimization, feasible-only sampling near the boundary makes the unconstrained interpolation picture unavailable, and the relevant notion is $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$0-full linearity on $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$1. The paper proves that both linear regression models and underdetermined quadratic interpolation models of minimum Frobenius norm achieve the same $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$2 function-error and $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$3 directional-gradient-error orders required by a convex-constrained trust-region method, provided the feasible sample set is $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$4-poised on the feasible region (Roberts, 2024).

6. Extensions beyond classical surface fitting

Some works broaden interpolation-regression approximation beyond direct function reconstruction. In neural ODE training, the Interpolated Reverse Dynamic Method stores exact forward states at Chebyshev nodes, reconstructs the trajectory by barycentric Lagrange interpolation,

$$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$5

and then substitutes $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$6 into the adjoint equations. The interpolated quantity is the latent state trajectory, not the gradient itself and not a learned regression surrogate. The paper explicitly describes this as purely interpolation-based, but also notes that, in a broader sense, it instantiates the pipeline “sample $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$7 interpolate $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$8 use in downstream gradient computation” (Daulbaev et al., 2020).

Local deterministic reconstruction methods offer another extension. One paper proposes a gradient-based local hyperplane approximation and a higher-accuracy smooth approximating-function correction for multidimensional regression, interpolation, reconstruction, and imputation. The construction uses only nearby points, solves local linear systems for directional derivatives, and in the smoother variant corrects those derivatives through coordinatewise cross-sectional curves. The reported experiments reach predictor dimension $$u_{\mathbb X,\mathbb Y}=K(\cdot,\mathbb X)K(\mathbb X,\mathbb X)^{-1}\mathbb Y,$$9, and the method is presented as a local geometric alternative to global fitting strategies for high-dimensional data (Shestopaloff et al., 2017).

On the exact-interpolation side, fractal constructions show that approximation families need not be polynomial or kernel-linear. Lidstone fractal interpolation constructs a $\frac12\|u\|_{\mathcal H_K}^2$0 self-referential interpolant $\frac12\|u\|_{\mathcal H_K}^2$1 satisfying

$\frac12\|u\|_{\mathcal H_K}^2$2

with $\frac12\|u\|_{\mathcal H_K}^2$3-error of order $\frac12\|u\|_{\mathcal H_K}^2$4 for the function and $\frac12\|u\|_{\mathcal H_K}^2$5 for the $\frac12\|u\|_{\mathcal H_K}^2$6-th derivative. A multivariate $\frac12\|u\|_{\mathcal H_K}^2$7-fractal operator on hyper-rectangles similarly maps a continuous function to a self-referential analogue $\frac12\|u\|_{\mathcal H_K}^2$8, with perturbation bound

$\frac12\|u\|_{\mathcal H_K}^2$9

These works remain interpolatory rather than regression-based, but they enlarge the approximation side of the interpolation-regression spectrum by introducing parameterized self-referential approximants (Kapoor et al., 2014, Pandey et al., 2021).

Taken together, the cited literature supports a precise synthesis. Interpolation-regression approximation is best understood as a continuum of methods that trade hard constraints, soft constraints, subset interpolation, least-squares residual minimization, and internal interpolatory surrogates against one another. Exact interpolation is often a limiting case of regression, regression can be implemented through interpolation-based parametrizations, and many of the decisive practical questions—node geometry, basis design, regularization, and conditioning—are shared by both regimes rather than belonging exclusively to one.