In this dissertation, we consider the approximation of unknown or intractable integrals using quadrature, especially when the evaluation of these integrals is considered very costly. This is a central problem both within and without machine learning, including model averaging, (hyper-)parameter marginalization, and computing posterior predictive distributions. Recently Batch Bayesian Quadrature has successfully combined the probabilistic integration techniques of Bayesian Quadrature with the parallelization techniques of Batch Bayesian Optimization, resulting in improved performance when compared to state-of-the-art Markov Chain Monte Carlo techniques, especially when parallelization is increased. While the selection of batches in Batch Bayesian Quadrature mitigates costs associated with individual point selection, every point within every batch is nevertheless chosen serially, which impedes the realization of the full potential of batch selection. We resolve this shortcoming. We have developed a novel Batch Bayesian Quadrature method which allows us to update points within a batch without incurring the costs traditionally associated with non-serial point selection. To implement this, we also devise a novel dynamic domain decomposition. Combining these, we show that this efficiently reduces uncertainty, leads to lower error estimates of the integrand, and therefore results in more numerically robust estimates of the integral. Furthermore, we close an open question in the Batch Bayesian Quadrature literature about the cessation criterion, which we determine and support using numerical methods on commonly used test functions in one and two dimensions. We present our findings within the context of the history of quadrature, show how our novel methods significantly improve what currently exists in the literature, and provide recommendations for future improvements.