  We present a model and variance reduction method for the fast and reliable
computation of statistical outputs of stochastic elliptic partial differential
equations. Our method consists of three main ingredients: (1) the hybridizable
discontinuous Galerkin (HDG) discretization of elliptic partial differential
equations (PDEs), which allows us to obtain high-order accurate solutions of
the governing PDE; (2) the reduced basis method for a new HDG discretization of
the underlying PDE to enable real-time solution of the parameterized PDE in the
presence of stochastic parameters; and (3) a multilevel variance reduction
method that exploits the statistical correlation among the different reduced
basis approximations and the high-fidelity HDG discretization to accelerate the
convergence of the Monte Carlo simulations. The multilevel variance reduction
method provides efficient computation of the statistical outputs by shifting
most of the computational burden from the high-fidelity
  HDG approximation to the reduced basis approximations. Furthermore, we
develop \textit{a posteriori} error estimates for our approximations of the
statistical outputs. Based on these error estimates, we propose an algorithm
for optimally choosing both the dimensions of the reduced basis approximations
and the sizes of Monte Carlo samples to achieve a given error tolerance. We
provide numerical examples to demonstrate the performance of the proposed
method.
