Large matrix inversions have often been cited as a major impediment to scaling Gaussian process (GP) models. With the use of GPs as building blocks for ever more sophisticated Bayesian deep learning models, removing these impediments is a necessary step for achieving large scale results. We present a variational approximation for a wide range of GP models that does not require a matrix inverse to be performed at each optimisation step. Our bound instead directly parameterises a free matrix, which is an additional variational parameter. At the local maxima of the bound, this matrix is equal to the matrix inverse. We prove that our bound gives the same guarantees as earlier variational approximations. We demonstrate some beneficial properties of the bound experimentally, although significant wall clock time speed improvements will require future improvements in optimisation and implementation. One major obstacle to the wider adoption of Gaussian Process (GP) (Rasmussen and Williams, 2006) based models is their computational cost, which is mainly caused by matrix inverses and determinants. Advances in variational approximate inference methods have reduced the size of the matrices on which expensive operations need to be performed, leading to O N M 2 time costs instead of O N 3 (Titsias, 2009), with approximations arbitrarily good with M N (Burt et al., 2019) . Minibatches of size B N can be used for training at a cost of O BM 2 + M 3 per iteration (Hensman et al., 2013) . The usefulness of training with small minibatches is hampered by the iteration cost being dominated by O M 3 , which again comes from an inverse and determinant. The computation is usually done using the Cholesky decomposition, which requires serial operations and high-precision arithmetic. So in addition to being an asymptotically expensive operation, it is also poorly suited to modern deep learning hardware. Removing these per-iteration matrix operations therefore seems necessary to speed up training. In this work, we provide a variational lower bound that can be computed without expensive matrix operations like inversion. Our bound can be used as a drop-in replacement to the existing variational method of Hensman et al. (2013 Hensman et al. ( , 2015 , and can therefore directly be applied in a wide variety of models, such as deep GPs (Damianou and Lawrence, 2013) . We focus on the theoretical properties of this new bound, and show some initial experimental results for optimising this bound. We hope to realise the full promise in scalability that this new bound has in future work. We presented new variational bounds for GP models that function as drop-in replacements to those developed by Hensman et al. (2013 Hensman et al. ( , 2015 , but without needing to compute expensive matrix operations each iteration. We prove their properties and show that they behave as expected in simple experiments using a single layer and deep GP. We believe this method to be promising, as it removes the most frequently cited impediment against the scaling of GP models. However, more improvements are needed to obtain the full practical benefits.
We present a variational lower bound for GP models that can be optimised without computing expensive matrix operations like inverses, while providing the same guarantees as existing variational approximations.