We introduce CUTS-GPR, a new method for performing numerically exact GPR in high-dimensional settings. The key component of CUTS-GPR is an extremely fast kernel matrix-vector product, which exhibits near-linear or even linear scaling with the amount of training data, N, and low-order polynomial scaling with dimensionality, D. This is obtained by combining an additive kernel with an incomplete grid and exploiting the resulting structure of the kernel matrix. The scalability of the matrix-vector product is verified by benchmarks with billions of data points and thousands of dimensions. We demonstrate the end-to-end scalability of CUTS-GPR by running full GPR calculations, including hyperparameter optimization, on synthetic datasets with up to N = 4 494 001 and D = 500. As a realistic and challenging test, we finally apply CUTS-GPR to a set of ten potential energy surfaces (PESs) with N = 447 265 and D = 24. The calculations are completed in a matter of hours, showing that CUTS-GPR enables Bayesian modelling of high-dimensional PESs - a longstanding challenge in computational chemistry.