  For regression tasks, standard Gaussian processes (GPs) provide natural
uncertainty quantification, while deep neural networks (DNNs) excel at
representation learning. We propose to synergistically combine these two
approaches in a hybrid method consisting of an ensemble of GPs built on the
output of hidden layers of a DNN. GP scalability is achieved via Vecchia
approximations that exploit nearest-neighbor conditional independence. The
resulting deep Vecchia ensemble not only imbues the DNN with uncertainty
quantification but can also provide more accurate and robust predictions. We
demonstrate the utility of our model on several datasets and carry out
experiments to understand the inner workings of the proposed method.
