pith. sign in

arxiv: 2505.12282 · v3 · submitted 2025-05-18 · 🧮 math.NA · cs.NA

Kernel interpolation on generalized sparse grids

Michael Griebel , Helmut Harbrecht , Michael Multerer This is my paper

Pith reviewed 2026-05-22 14:51 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords kernel interpolationsparse gridsscattered datasamplet compressioncombination techniqueproduct regionserror estimateslarge-scale problems
0
0 comments X

The pith

Generalized sparse grids with samplet compression enable kernel interpolation on up to a billion scattered points.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops a method for kernel interpolation of scattered data on product regions whose component domains may have equal or unequal dimensions. It constructs optimal sparse grids for these settings and obtains sharpened error bounds through duality arguments. The solution procedure rests on the sparse-grid combination technique, where each low-dimensional anisotropic tensor-product subproblem is solved after samplet compression turns the univariate kernel matrices into essentially sparse form. The resulting algorithm scales to interpolation problems with billions of points, especially when the underlying reproducing kernels are nonlocal. Readers would care because conventional dense kernel methods quickly exhaust memory and time on large or high-dimensional scattered datasets.

Core claim

We consider scattered data approximation on product regions of equal and different dimensionality. On each of these regions, we assume quasi-uniform but unstructured data sites and construct optimal sparse grids for scattered data interpolation on the product region. For this, we derive new improved error estimates for the respective kernel interpolation error by invoking duality arguments. An efficient algorithm to solve the underlying linear system of equations is proposed. The algorithm is based on the sparse grid combination technique, where a sparse direct solver is used for the elementary anisotropic tensor product kernel interpolation problems. The application of the sparse direct sol

What carries the argument

Sparse-grid combination technique applied to anisotropic tensor-product kernel interpolants after samplet compression renders each univariate kernel matrix essentially sparse.

If this is right

  • Optimal error rates are recovered for the generalized sparse-grid setting even when the component domains differ in dimension.
  • Dense nonlocal kernel matrices no longer prevent solution of systems with 10^9 points.
  • The combination technique plus compression yields a practical direct solver for each elementary tensor-product subproblem.
  • The same framework applies uniformly to both equal-dimensional and mixed-dimensional product regions.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The compression step may transfer to other nonlocal integral operators that produce dense matrices in high dimensions.
  • One could replace samplets by other hierarchical bases and measure the effect on the overall complexity.
  • Product-region problems of this type appear in uncertainty quantification, suggesting a route to scalable surrogate models there.

Load-bearing premise

The data sites on each sub-region of the product domain are quasi-uniform though unstructured.

What would settle it

A run of the algorithm on a ten-dimensional product domain with a nonlocal kernel and exactly one billion scattered points, checking whether memory and runtime remain practical while the observed error matches the derived duality bound.

Figures

Figures reproduced from arXiv: 2505.12282 by Helmut Harbrecht, Michael Griebel, Michael Multerer.

Figure 1
Figure 1. Figure 1: Visualization of the subsampling procedure starting from a set of 1000 uniformly chosen random points on [0, 1]2 . Different sophisticated algorithms for the construction of nested subsets from a given set of data sites have been proposed in the literature, see, e.g., [6, 9, 40]. Nonetheless, for our our purposes, the simple algorithm described below is sufficient. To construct a multilevel sequence (2.8) … view at source ↗
Figure 2
Figure 2. Figure 2: Sparsity patterns of the samplet compressed exponential ker￾nel on the unit square (left) for 300 000 data sites, the nested dissection reordering (middle), and the Cholesky factor (right). Each dot represents a matrix block of size 300 × 300. The number of entries per block is color coded, where lighter blocks have less entries. The sparse direct solver mitigates to some extent the computational cost for … view at source ↗
Figure 3
Figure 3. Figure 3: Computation times for the canonical sparse grid on the unit hypercube (0, 1)m and m = 3, 6, 9, 12, 15, 18 [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Convergence of the kernel interpolant on the canonical sparse grid in (0, 1)m. In [PITH_FULL_IMAGE:figures/full_fig_p016_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Sketch of the regular grid points (blue) and the evaluation points (red) on the unit interval, the unit square, and the unit cube. 0 1 2 3 4 5 6 100 10−1 10−2 10−3 10−4 10−5 10−6 10−7 level j relative error d = 1 d = 2 d = 3 h −3.125 j [PITH_FULL_IMAGE:figures/full_fig_p017_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Convergence of the d-variate kernel approximation in case of the hypercube [0, 1]d for d = 1, 2, 3. 100 101 102 103 104 105 106 107 108 100 10−1 10−2 10−3 10−4 10−5 number N of degrees of freedom relative error accuracy equilibrated cost-benefit equilibrated DoF equilibrated N −1.04 [PITH_FULL_IMAGE:figures/full_fig_p017_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Convergence rates of the kernel approximant with respect to different sparse grids on (0, 1) × (0, 1)2 × (0, 1)3 . Ωi is shown in [PITH_FULL_IMAGE:figures/full_fig_p017_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Sketch of the quasi-uniform points (blue) and evaluation points (red) on the unit interval, the unit sphere, and the Stanford bunny. We next consider the kernel interpolation of the respective sparse grid. For the present setting, we can summarize the parameters as (5.2) d1 = 1, d2 = 2, d3 = 3, s1 = s2 = s3 = 25 16 . Therefore, choosing the weights (5.3) w1 = w2 = w3 = 1 for the sparse grid construction eq… view at source ↗
Figure 9
Figure 9. Figure 9: Convergence of the kernel interpolant on the interval (d = 1), the sphere (d = 2), and the rabbit (d = 3). 100 101 102 103 104 105 106 107 100 10−1 10−2 10−3 Number N of degrees of freedom Relative error accuracy equilibrated cost-benefit equilibrated DoF equilibrated N −1.04 [PITH_FULL_IMAGE:figures/full_fig_p020_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Convergence rates of the kernel approximant with respect to different sparse grids on the product of general subregions in 1 + 2 + 3 dimensions. the cost-benefit-equilibrated sparse grid is given by the weights in (5.5). As can be inferred from [PITH_FULL_IMAGE:figures/full_fig_p020_10.png] view at source ↗
read the original abstract

We consider scattered data approximation on product regions of equal and different dimensionality. On each of these regions, we assume quasi-uniform but unstructured data sites and construct optimal sparse grids for scattered data interpolation on the product region. For this, we derive new improved error estimates for the respective kernel interpolation error by invoking duality arguments. An efficient algorithm to solve the underlying linear system of equations is proposed. The algorithm is based on the sparse grid combination technique, where a sparse direct solver is used for the elementary anisotropic tensor product kernel interpolation problems. The application of the sparse direct solver is facilitated by applying a samplet matrix compression to each univariate kernel matrix, resulting in an essentially sparse representation of the latter. In this way, we obtain a method that is able to deal with large problems up to billions of interpolation points, especially in case of reproducing kernels of nonlocal nature. Numerical results are presented to qualify and quantify the approach.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 3 minor

Summary. The manuscript develops kernel interpolation on generalized sparse grids for scattered data approximation on product regions of equal or different dimensionality. Assuming quasi-uniform but unstructured data sites on each factor, it constructs optimal sparse grids and derives new error estimates for the kernel interpolant via duality arguments. An efficient solver is proposed that applies the sparse-grid combination technique to anisotropic tensor-product subproblems, each solved by a sparse direct method after samplet compression of the univariate kernel matrices; this is claimed to enable treatment of up to billions of points, especially for nonlocal kernels. Supporting numerical results are presented.

Significance. If the error analysis extends rigorously to the compressed subproblems and the combination step, the work would provide a practical route to high-dimensional kernel interpolation at scales far beyond current direct methods, particularly for nonlocal reproducing kernels. The combination of duality-based a-priori estimates with the combination technique and samplet compression constitutes a clear algorithmic and analytic advance; the paper also supplies concrete numerical evidence of scalability.

major comments (2)
  1. [§4] §4 (error analysis): The duality argument establishing the interpolation error bound is stated for the uncompressed kernel matrices. It is not shown that the perturbation introduced by samplet compression of the univariate factors remains controlled under the combination technique, nor whether an additional consistency term must be added to the product-region error bound to preserve the claimed rate for nonlocal kernels.
  2. [§5.2] §5.2 (algorithm and implementation): The choice of compression tolerance is described only heuristically. No explicit relation is given between the tolerance, the univariate fill distance, and the final sparse-grid error, leaving open whether the headline claim of “no accuracy loss” at billion-point scale holds uniformly.
minor comments (3)
  1. [Abstract] The abstract asserts “new improved error estimates” without indicating the precise improvement (e.g., removal of a logarithmic factor or extension to anisotropic quasi-uniformity).
  2. [§2] Notation for the product domains and the associated sparse-grid index sets would benefit from an explicit diagram or table in §2.
  3. [§3] A few typographical inconsistencies appear in the displayed equations of §3 (e.g., missing parentheses around inner products).

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading of the manuscript and the constructive comments. We address the major points below and will revise the manuscript to incorporate additional analysis where needed.

read point-by-point responses

  1. Referee: [§4] §4 (error analysis): The duality argument establishing the interpolation error bound is stated for the uncompressed kernel matrices. It is not shown that the perturbation introduced by samplet compression of the univariate factors remains controlled under the combination technique, nor whether an additional consistency term must be added to the product-region error bound to preserve the claimed rate for nonlocal kernels.

    Authors: We agree that the duality-based error estimates presented in Section 4 apply directly to the exact (uncompressed) kernel interpolant. The samplet compression is introduced later as part of the solver to achieve scalability for large problems. Our numerical results indicate that the chosen compression tolerances preserve the observed accuracy, but we acknowledge that a rigorous perturbation analysis through the combination technique is not provided. In the revised manuscript we will add a brief analysis (or remark) bounding the compression-induced consistency error and showing that it can be made smaller than the interpolation error term by appropriate choice of tolerance, thereby preserving the claimed rates. revision: yes

  2. Referee: [§5.2] §5.2 (algorithm and implementation): The choice of compression tolerance is described only heuristically. No explicit relation is given between the tolerance, the univariate fill distance, and the final sparse-grid error, leaving open whether the headline claim of “no accuracy loss” at billion-point scale holds uniformly.

    Authors: The referee is correct that the selection of the compression tolerance is currently described in heuristic terms. In practice the tolerance is chosen relative to the univariate fill-distance error so that compression does not dominate the overall error. We will revise Section 5.2 to state an explicit guideline (or sufficient condition) relating the tolerance to the univariate fill distance and the target sparse-grid error, making the regime in which “no accuracy loss” holds more precise and uniform. revision: yes

Circularity Check

0 steps flagged

No significant circularity; derivation relies on independent duality arguments and combination technique.

full rationale

The paper derives improved error estimates via duality arguments for kernel interpolation on quasi-uniform sites, then applies the sparse grid combination technique with samplet compression for efficiency on large-scale problems. These steps are presented as independent mathematical constructions and algorithmic choices, not as redefinitions or fits of the target scalability claim. No quoted reduction shows a prediction or result equaling its own input by construction, and self-citations (if present) are not load-bearing for the central claims. The approach remains self-contained against external benchmarks for error control and compression.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Abstract-only review prevents identification of specific free parameters, axioms, or invented entities; none are explicitly introduced in the summary text.

pith-pipeline@v0.9.0 · 5680 in / 1016 out tokens · 47456 ms · 2026-05-22T14:51:03.629476+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

47 extracted references · 47 canonical work pages

  1. [1]

    Adams and J

    R. Adams and J. Fournier. Sobolev Spaces. Academic Press, Amsterdam, 2003

    work page 2003

  2. [2]

    A vesani, R

    S. A vesani, R. Kempf, M. Multerer, and H. Wendland. Multiscale scattered data analysis in samplet coordinates. SIAM J. Sci. Comput. , 47(5):A3038–A3063, 2025

    work page 2025

  3. [3]

    Bartel, M

    F. Bartel, M. Schäfer, and T. Ullrich. Constructive subsampling of finite frames with applications in optimal function recovery. Appl. Comput. Harmon. Anal. , 65:209–248, 2023

    work page 2023

  4. [4]

    Bungartz and M

    H.-J. Bungartz and M. Griebel. Sparse grids. Acta Numer. , 13:147–269, 2004

    work page 2004

  5. [5]

    Y. Chen, T. Davis, W. Hager, and S. Rajamanickam. Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Trans. Math. Software , 35(3):22:1–22:14, 2008

    work page 2008

  6. [6]

    De Floriani

    L. De Floriani. A pyramidal data structure for triangle-based surface description. IEEE Comput. Graph. Appl., 9(2):67–78, 1989

    work page 1989

  7. [7]

    F.-J. Delvos. d-variate Boolean interpolation. J. Approx. Theory , 34(2):99–114, 1982

    work page 1982

  8. [8]

    Delvos and W

    F.-J. Delvos and W. Schempp. Boolean Methods in Interpolation and Approximation . Longman Sci- entific & Technical, Harlow, 1989

    work page 1989

  9. [9]

    Floater and A

    M. Floater and A. Iske. Multistep scattered data interpolation using compactly supported radial basis functions. J. Comput. Appl. Math. , 73(1):65–78, 1996

    work page 1996

  10. [10]

    Dolbeault, D

    M. Dolbeault, D. Krieg, and M. Ullrich. A sharp upper bound for sampling numbers in L2. Appl. Comput. Harmon. Anal. , 63:113–134, 2023

    work page 2023

  11. [11]

    Dũng , V

    D. Dũng , V. Nguyen, C. Schwab, and J. Zech. Analyticity and Sparsity in Uncertainty Quantification for PDEs with Gaussian Random Field Inputs , volume 2334 of Lecture Notes in Mathematics , Springer Nature Switzerland AG, 2023

    work page 2023

  12. [12]

    D. Dũng, V. Temlyakov, and T. Ullrich. Hyperbolic Cross Approximation. Advanced Courses in Math- ematics. CRM Barcelona. Birkhäuser/Springer Basel, 2018

    work page 2018

  13. [13]

    Fasshauer

    G. Fasshauer. Meshfree Approximation Methods with MATLAB . World Scientific, River Edge, NJ, 2007

    work page 2007

  14. [14]

    Feldstein, D

    A. Feldstein, D. Lazzara, N. Princen, and K. Willcox. Multifidelity data fusion: Application to blended-wing-body multidisciplinary analysis under uncertainty. AIAA J. , 58(2):889–906, 2020

    work page 2020

  15. [15]

    A. George. Nested dissection of a regular finite element mesh. SIAM J. Numer. Anal. , 10(2):345–363, 1973

    work page 1973

  16. [16]

    Giebermann

    K. Giebermann. Multilevel approximation of boundary integral operators. Computing, 67(3):183–-207, 2001. KERNEL INTERPOLATION ON GENERALIZED SPARSE GRIDS 23

    work page 2001

  17. [17]

    Greengard and V

    L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comput. Phys. , 73(2):325– 348, 1987

    work page 1987

  18. [18]

    Griebel and H

    M. Griebel and H. Harbrecht. On the construction of sparse tensor product spaces. Math. Comput. , 82(282):975–994, 2013

    work page 2013

  19. [19]

    Griebel and H

    M. Griebel and H. Harbrecht. A note on the construction of L-fold sparse tensor product spaces. Constr. Approx., 38(2):235–251, 2013

    work page 2013

  20. [20]

    Griebel, M

    M. Griebel, M. Schneider and C. Zenger. A combination technique for the solution of sparse grid problems. In P. de Groen and R. Beauwens, editors, Iterative Methods in Linear Algebra , pages 263–

    work page

  21. [21]

    IMACS, Elsevier, North Holland, 1992

    work page 1992

  22. [22]

    Hackbusch

    W. Hackbusch. Hierarchische Matrizen: Algorithmen und Analysis . Springer, Berlin-Heidelberg, 2009

    work page 2009

  23. [23]

    Hackbusch

    W. Hackbusch. Tensor Spaces and Numerical Tensor Calculus . 2nd ed., Springer, Berlin-Heidelberg, 2019

    work page 2019

  24. [24]

    Haji-Ali, H

    A.-L. Haji-Ali, H. Harbrecht, M. Peters, and M. Siebenmorgen. Novel results for the anisotropic sparse grid quadrature. J. Complexity , 47:62–-85, 2018

    work page 2018

  25. [25]

    Harbrecht and M

    H. Harbrecht and M. Multerer. Samplets: Construction and scattered data compression. J. Comput. Phys., 471:111616, 2022

    work page 2022

  26. [26]

    Harbrecht and M

    H. Harbrecht and M. Multerer. Samplets: Wavelet concepts for scattered data. In R. DeVore and A. Kunoth, editors, Multiscale, Nonlinear and Adaptive Approximation II , pages 299–326. Springer, Berlin-Heidelberg, 2024

    work page 2024

  27. [27]

    Harbrecht, M

    H. Harbrecht, M. Multerer, and J. Quizi. The dimension weighted fast multipole method for scattered data approximation. J. Comput. Phys. , 532:113956, 2025

    work page 2025

  28. [28]

    Harbrecht, M

    H. Harbrecht, M. Multerer, O. Schenk, and C. Schwab. Multiresolution kernel matrix algebra. Numer. Math., 156:1085–1114, 2024

    work page 2024

  29. [29]

    Harbrecht, M

    H. Harbrecht, M. Peters, and M. Siebenmorgen. Combination technique based k-th moment analysis of elliptic problems with random diffusion. J. Comput. Phys. , 252:128–141, 2013

    work page 2013

  30. [30]

    Hertrich

    J. Hertrich. Fast kernel summation in high dimensions via slicing and Fourier transforms. SIAM J. Math. Data Sci. , 6(2):440–461, 2024

    work page 2024

  31. [31]

    R. Kempf. The Tensor Product Multilevel Method for High-dimensional Meshfree Approximation . PhD thesis, Universität Bayreuth, 2023

    work page 2023

  32. [32]

    Kempf and H

    R. Kempf and H. Wendland. High-dimensional approximation with kernel-based multilevel methods on sparse grids. Numer. Math. , 154:485–519, 2023

    work page 2023

  33. [33]

    B. Matérn. Spatial Variation, volume 36 of Lecture Notes in Statistics. Springer, Berlin, second edition, 1986

    work page 1986

  34. [34]

    Opschoor, C

    J. Opschoor, C. Schwab, and J. Zech. Deep learning in high dimension: ReLU neural network expres- sion for Bayesian PDE inversion De Gruyter, Berlin, 2022

    work page 2022

  35. [35]

    Peherstorfer, K

    B. Peherstorfer, K. Willcox, and M. Gunzburger. Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Rev. , 60(3):550–591, 2018

    work page 2018

  36. [36]

    Potts and G

    D. Potts and G. Steidl. Fast summation at non-equispaced knots by NFFT. SIAM J. Sci. Comput. , 24(6):2013–2037, 2003

    work page 2013

  37. [37]

    Rahimi and B

    A. Rahimi and B. Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems (NeurIPS) , pages 1177–1184. 2007

    work page 2007

  38. [38]

    Rasmussen and C

    C. Rasmussen and C. Williams. Gaussian Processes for Machine Learning . MIT Press, Cambridge, MA, 2006

    work page 2006

  39. [39]

    V. Rychkov. On restrictions and extension of the Besov and Triebel-Lizorkin spaces with respect to Lipschitz domains. J. London Math. Soc. 60(1):237–257, 1999

    work page 1999

  40. [40]

    Schaback

    R. Schaback. Superconvergence of kernel-based interpolation. J. Approx. Theory , 235:1–19, 2018

    work page 2018

  41. [41]

    Schlömer, D

    T. Schlömer, D. Heck, and O. Deussen. Farthest-point optimized point sets with maximized minimum distance. In Proc. ACM SIGGRAPH Symp. High Performance Graphics (HPG) , pages 135–142. ACM, Vancouver, 2011

    work page 2011

  42. [42]

    Schölkopf and A

    B. Schölkopf and A. Smola. Learning with Kernels: Support Vector Machines, Regularization, Opti- mization, and Beyond . MIT Press, Cambridge, MA, 2002

    work page 2002

  43. [43]

    Sloan and V

    I. Sloan and V. Kaarnioja. Doubling the rate: Improved error bounds for orthogonal projection with application to numerical analysis. BIT Numer. Math. , 65:10, 2025

    work page 2025

  44. [44]

    S. Smolyak. Quadrature and interpolation formulas for tensor products of certain classes of functions. Doklady Akademii Nauk SSSR , 4:240–243, 1963

    work page 1963

  45. [45]

    Theodoridis and K

    S. Theodoridis and K. Koutroumbas. Pattern Recognition. Academic Press, Waltham, MA, fourth edition, 2009

    work page 2009

  46. [46]

    Wendland

    H. Wendland. Scattered Data Approximation. Cambridge University Press, Cambridge, 2004

    work page 2004

  47. [47]

    C. Zenger. Sparse grids. In Parallel algorithms for partial differential equations , Proceedings of the 6th GAMM-Seminar, Kiel/Germany, 1990, Notes Numer. Fluid Mech. 31, pages 241–251, Vieweg, Braunschweig, 1991. 24 MICHAEL GRIEBEL, HELMUT HARBRECHT, AND MICHAEL MULTERER Michael Griebel, Institut für Numerische Simulation, Universität Bonn, Friedrich-Hir...

    work page 1990