The design of interpolatory quadrature rules with positive weights is of great interest in approximation theory and scientific computing: Such quadrature rules achieve near-optimal approximation of integrals and are associated with well-behaved Lebesgue constants. Such quadrature rules are used heavily in applications like uncertainty quantification and computational finance.
Computing such quadrature rules in more than one dimension is an arduous task, frequently attempted with non-convex optimization schemes. The success or failure of such approaches typically depends on the type of domain, dimension, or approximation space. We present a new procedure whose implementation is a probabilistic algorithm based on a novel constructive proof of Tchakaloff's theorem. In particular, we can successfully compute size-N quadrature rules in high dimensions for general approximation spaces. The main feature of our algorithm is a complexity that depends algebraically on N, and in some cases is strictly linear in the dimension. We illustrate that our procedure is effective even for quadrature on complex and unbounded multivariate domains.