SABER: Stable adaptive barycentric extension for reconstruction of high-dimensional functions

In this paper, we develop a sparse barycentric interpolation method for high-dimensional function approximation. The method constructs a sparse barycentric operator through tensor contraction on Chebyshev-Gauss-Lobatto grids. By extending Higham’s backward error analysis to multiple dimensions, we establish that floating-point errors are strictly bounded by the number of local tensor-product contributions and the Lebesgue constants of the component one-dimensional interpolants. To maintain numerical stability under adaptive refinement, we introduce a dynamic node insertion strategy governed by a Leja geometric repulsion penalty. The approach can balance residual maximization with geometric regularization, controlling the growth of the Lebesgue constant and the condition number of the barycentric weights. For discontinuous functions, we investigate an adaptive-fiber dimension-wise reconstruction approach by using the idea from WENO. Extensive numerical experiments, including highly anisotropic benchmarks and three-dimensional curved spherical shocks, demonstrate the efficiency of the proposed method.