Multiscale Polynomial Interpolation

• The paper introduces a multiscale polynomial interpolation framework that integrates dilation and translation operators with Chebyshev bases to achieve robust high-accuracy function approximation.
• It employs a structured matrix pencil approach to simultaneously recover sparse polynomial degrees and coefficients through generalized eigenvalue problems.
• The methodology extends to quantized tensor train approximations and finite element methods, offering adaptive, super-resolution, and efficient solutions for high-dimensional problems.

A multiscale polynomial interpolation perspective provides a unified framework for constructing, analyzing, and applying hierarchical polynomial approximations across various domains, exploiting the structure and properties of polynomial bases (notably Chebyshev) at multiple resolution scales. This perspective systematically integrates dilation and translation operations, multilevel hierarchies, and structured linear algebra (matrix pencils, tensor decompositions) to enable robust, efficient, and high-accuracy interpolation or representation of functions and data—especially in high-dimensional and multiscale settings.

1. Theoretical Foundations and Multiscale Operators

Multiscale polynomial interpolation extends classical sparse polynomial interpolation by introducing dilation and translation operators on polynomial bases, enabling the analysis and reconstruction of functions sampled at different scales and shifts. For the Chebyshev first-kind basis $T_m(t) = \cos(m\theta)$ with $t = \cos\theta$, the fundamental model takes the form

$f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$

where $\alpha_i$ are unknown coefficients and $\varphi_i$ are unknown (typically integer) degrees. Given samples $f(t_j)$ on a uniform grid $t_j = \cos(j\Delta)$, the two primary multiscale operators are:

• Dilation ($s$): Samples $f(t_{js})$ correspond to multiplying basis elements by $\cos(m s \Delta)$, acting as a "frequency decimation" analogous to Prony-type methods.
• Translation ($t = \cos\theta$0): Shifted samples $t = \cos\theta$1 encode additional spectral/degree information via $t = \cos\theta$2, breaking degeneracies and enabling alias resolution.

This operator framework supports function analysis and interpolation at different scales (i.e., "multiscale") and provides the algebraic structure necessary for stable numerical methods (Cuyt et al., 2020).

2. Structured Matrix Pencil Approach for Sparse Polynomial Interpolation

The multiscale matrix pencil method reformulates sparse polynomial interpolation as a structured generalized eigenvalue problem, allowing simultaneous recovery of the sparsity pattern (polynomial degrees) and coefficients. Given chosen integers $t = \cos\theta$3 (with $t = \cos\theta$4), construct $t = \cos\theta$5 Hankel-type and shifted-Hankel matrices $t = \cos\theta$6, whose entries are symmetrized combinations of function samples: $t = \cos\theta$7

$t = \cos\theta$8

where $t = \cos\theta$9. These matrices admit factorizations involving Vandermonde-type matrices and diagonal matrices of coefficients and cosines of the degrees. The key computational step is the solution of the generalized eigenvalue problem

$f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$0

with eigenvalues $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$1. The degrees $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$2 are uniquely determined by de-aliasing using both the dilation and translation eigenstructure, and coefficients $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$3 are obtained via structured Vandermonde-type linear solves. Optimal choices of $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$4 (large, co-prime to $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$5) improve conditioning when degrees are clustered, yielding super-resolution and robust performance (Cuyt et al., 2020).

3. Multiscale Interpolation in Quantized Tensor Train (QTT) Approximations

The QTT decomposition for functions on dyadic grids is rigorously connected to a multiscale polynomial interpolation hierarchy. Each point $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$6 is represented via a binary expansion $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$7. Multilevel decompositions correspond to nested partitioning of the interval into $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$8 dyadic subintervals at scale $f(t) = \sum_{i=1}^n \alpha_i T_{\varphi_i}(t)$9, on each of which local Chebyshev interpolation is applied: $\alpha_i$0 Here, $\alpha_i$1 are Chebyshev-Lobatto points and $\alpha_i$2 are cardinal basis polynomials. The resulting multidimensional tensor $\alpha_i$3 admits a quantized tensor train representation whose ranks at each scale are tightly controlled by the local interpolation error $\alpha_i$4.

This construction explains and quantifies why QTT ranks decay exponentially or superalgebraically with depth for analytic, smooth, or bandlimited functions, and also why sharp features can be efficiently represented via adaptive multiscale interpolation, with dyadic grid adaptivity focused on "dangerous" intervals where the function is non-smooth (Lindsey, 2023).

4. Error Estimates, Rank Bounds, and Adaptivity

The error analysis of multiscale polynomial interpolation rests on standard approximation theory in Chebyshev bases. For $\alpha_i$5, the interpolation error in subintervals of size $\alpha_i$6 is $\alpha_i$7 for degree $\alpha_i$8 interpolation. For analytic $\alpha_i$9 (e.g., with a Bernstein-ellipse parameter $\varphi_i$0) or $\varphi_i$1-bandlimited $\varphi_i$2, the error decays exponentially in $\varphi_i$3 or as $\varphi_i$4. The QTT core ranks $\varphi_i$5 at depth $\varphi_i$6 are bounded by $\varphi_i$7, and more refined rank bounds tie $\varphi_i$8 to local smoothness or bandlimit:

• $\varphi_i$9: $f(t_j)$0
• Bandlimited $f(t_j)$1: $f(t_j)$2

In practice, the rank at finer scales ($f(t_j)$3) decays, reflecting that $f(t_j)$4 is locally nearly polynomial or constant at small length scales. Adaptive algorithms monitor local smoothness or features and refine grid sizes or raise interpolation order only where necessary (e.g., on intervals containing singularities) (Lindsey, 2023).

5. Multiscale Polynomial Interpolation in Numerical Multiphysics and Electromagnetics

In complex multiphysics settings, multiscale polynomial interpolation underpins modern finite element and spectral-element methods for PDE-based simulation. The Hermite Finite Element Method (HFEM) leverages nodal Hermite interpolation polynomials—of high order and $f(t_j)$5 continuity—across all mesh elements. The degrees of freedom in each element are function values and derivatives up to second order at nodes, ensuring continuity and avoiding spurious modes common in vector field discretizations. The group-theoretic basis construction (e.g., by Kassebaum–Boucher–Ram-Mohan) enables systematic high-order, symmetry-adapted, $f(t_j)$6-continuous element bases.

Key properties include:

• Exact enforcement of divergence constraints (suppression of spurious solutions),
• Spectral accuracy for smooth problems with rapid $f(t_j)$7-convergence,
• Flexibility for complex geometries via arbitrary triangular/tetrahedral meshes,
• Strictly local mesh refinement without global accuracy degradation,
• Computational scaling that is highly favorable against traditional vector FEM or Fourier (plane-wave) methods (Pham et al., 2019).

6. Algorithmic Realizations and Computational Complexity

Multiscale polynomial interpolation admits several algorithmic instantiations:

• Matrix pencil eigendecomposition: Dominated by $f(t_j)$8 complexity for $f(t_j)$9-term interpolation, further reducible via structure to $t_j = \cos(j\Delta)$0 or $t_j = \cos(j\Delta)$1.
• QTT core construction: Requires $t_j = \cos(j\Delta)$2 function evaluations and $t_j = \cos(j\Delta)$3 per level for dense TTs, or $t_j = \cos(j\Delta)$4 for rank-adaptive constructions.
• Rank-revealing and sparse Chebyshev core methods: Truncated SVDs at each level, error propagation controlled by moderate Lebesgue constants, sparsity exploited via localized interpolation supports efficient computation and storage.
• HFEM assembly and solution: Yields sparse, banded Hermitian matrices with fill-factors $t_j = \cos(j\Delta)$5–$t_j = \cos(j\Delta)$6 and global solves scaling as $t_j = \cos(j\Delta)$7 compared to $t_j = \cos(j\Delta)$8–$t_j = \cos(j\Delta)$9 for plane-wave or boundary element methods.

These algorithmic frameworks are inherently black-box, stable, and rank-adaptive. For sparse polynomial interpolation, model order $s$0 can be found from matrix pencil rank, while for QTT and HFEM, adaptivity is achieved by error monitoring and local refinement (Cuyt et al., 2020, Lindsey, 2023, Pham et al., 2019).

7. Applications and Impact Across Computational Science

Multiscale polynomial interpolation perspectives have driven advances in:

• Separable nonlinear inverse problems: Enabling recovery of both nonlinear (e.g., frequency/degree) and linear coefficients via stable, explicitly structured solvers (Cuyt et al., 2020).
• High-dimensional function representation: QTT decompositions and tensor train methods tailored by local smoothness and features, essential for high-dimensional numerical analysis and uncertainty quantification (Lindsey, 2023).
• Numerical eigenvalue problems: HFEM achieves superior accuracy with minimal DoF for waveguides, resonators, photonic crystals, and multiphysics coupling, demonstrating the power of high-order, derivative-continuous, multiscale polynomial bases (Pham et al., 2019).

A plausible implication is that the multiscale polynomial interpolation framework is foundational to future developments in black-box, adaptive, and efficient algorithms for both classical and emerging scientific computing challenges involving structured high-dimensional data and PDEs.