  Gaussian process hyperparameter optimization requires linear solves with, and
log-determinants of, large kernel matrices. Iterative numerical techniques are
becoming popular to scale to larger datasets, relying on the conjugate gradient
method (CG) for the linear solves and stochastic trace estimation for the
log-determinant. This work introduces new algorithmic and theoretical insights
for preconditioning these computations. While preconditioning is well
understood in the context of CG, we demonstrate that it can also accelerate
convergence and reduce variance of the estimates for the log-determinant and
its derivative. We prove general probabilistic error bounds for the
preconditioned computation of the log-determinant, log-marginal likelihood and
its derivatives. Additionally, we derive specific rates for a range of
kernel-preconditioner combinations, showing that up to exponential convergence
can be achieved. Our theoretical results enable provably efficient optimization
of kernel hyperparameters, which we validate empirically on large-scale
benchmark problems. There our approach accelerates training by up to an order
of magnitude.
