Multivariate Newton interpolation in downward closed spaces reaches the optimal Bernstein–Walsh approximation rate

Abstract Recent advances in Bernstein—Walsh theory have extended Bernstein’s Theorem to multiple dimensions, stating that a multivariate function can be approximated with a geometric rate in a downward-closed polynomial space if and only if it is analytic in a generalized Bernstein polyellipse. To compute approximations of this class of functions—which we term Bos–Levenberg–Trefethen–(BLT) functions—we extend the classic univariate Newton interpolation algorithm to arbitrary multivariate downward-closed polynomial spaces, while maintaining its quadratic runtime and linear storage complexity. The present generalization supports any choice of (nontensorial) unisolvent interpolation nodes, whose number coincides with the dimension of the chosen downward-closed space. We prove that by selecting Leja nodes, the optimal geometric approximation rates for BLT-functions are achieved and that these rates extend to the derivatives of the interpolants. Choosing a Euclidean degree results in downward-closed spaces whose dimension only grows sub-exponentially with spatial dimension, while delivering approximation rates close to or matching those of the tensorial maximum-degree case, hence mitigating the curse of dimensionality. Importantly, our constructive proof directly inspires an algorithm for multivariate polynomial interpolation. We implemented this algorithm as a Python package and use it here to validate our theoretical findings in numerical experiments. These experiments corroborate the superiority of multivariate Newton interpolation over state-of-the-art alternatives, and they suggest that Leja-ordered Chebyshev–Lobatto nodes offer the same approximation power as Leja nodes.