Efficient and Robust Carathéodory-Steinitz Pruning of Positive Discrete Measures
Filip Bělík, Jesse Chan, Akil Narayan
TL;DR
The paper tackles efficient compression of positive discrete measures to achieve moment-preserving quadrature by extending Carathéodory-Steinitz pruning to streaming settings. It introduces SCSP and its Givens-rotations-based enhancement, GSCSP, achieving storage O($N^2$) and per-iteration complexity O($N^2$) with overall cost roughly linear in the original support size M. A formal stability analysis shows local Lipschitz continuity in total-variation distance under admissible perturbations, and numerical results demonstrate scalability to extremely large rules (up to M = 10^9) and practical use in cut-cell discontinuous Galerkin methods. The work provides competitive performance against NNLS and LP while offering improved stability and storage robustness, and it releases CaratheodoryPruning.jl for reproducibility.
Abstract
In many applications, one seeks to approximate integration against a positive measure of interest by a positive discrete measure: a numerical quadrature rule with positive weights. One common desired discretization property is moment preservation over a finite dimensional function space, e.g., bounded-degree polynomials. Carathéodory's theorem asserts that if there is any finitely supported quadrature rule with more nodes than the dimension of the given function space, one can form a smaller (and hence more efficient) positive, nested, quadrature rule that preserves the moments of the original rule. We describe an efficient streaming procedure for Carathéodory-Steinitz pruning, a numerical procedure that implements Carathéodory's theorem for this measure compression. The new algorithm makes use of Givens rotations and on-demand storage of arrays to successfully prune very large rules whose storage complexity only depends on the dimension of the function space. This approach improves on a naive implementation of Carathéodory-Steinitz pruning whose runtime and storage complexity are quadratic and linear, respectively, in the size of the original measure. We additionally prove mathematical stability properties of our method with respect to a set of admissible, total-variation perturbations of the original measure. Our method is compared to two alternate approaches with larger storage requirements: non-negative least squares and linear programming, and we demonstrate comparable runtimes, with improved stability and storage robustness. Finally, we demonstrate practical usage of this algorithm to generate quadrature for discontinous Galerkin finite element simulations on cut-cell meshes.
