Bayesian probabilistic numerical methods for numerical integration offer
significant advantages over their non-Bayesian counterparts: they can encode
prior information about the integrand, and can quantify uncertainty over
estimates of an integral. However, the most popular algorithm in this class,
Bayesian quadrature, is based on Gaussian process models and is therefore
associated with a high computational cost. To improve scalability, we propose
an alternative approach based on Bayesian neural networks which we call
Bayesian Stein networks. The key ingredients are a neural network architecture
based on Stein operators, and an approximation of the Bayesian posterior based
on the Laplace approximation. We show that this leads to orders of magnitude
speed-ups on the popular Genz functions benchmark, and on challenging problems
arising in the Bayesian analysis of dynamical systems, and the prediction of
energy production for a large-scale wind farm