Gaussian processes are popular tools for prediction that have been shown to be orders of magnitude more accurate than modern competitors on a host of prediction tasks. However, the computational cost of fitting them can be daunting. Inspired by the recent deployments of large-scale optimization in deep learning, this talk illustrates how carefully written optimization problems can be used to replace the usual matrix decomposition used to fit Gaussian process predictors. This can be used throughout hyperparameter tuning and is extremely scalable with modern linear algebra libraries. A computational deployment of the methodology has been used on large-scale machines and produces significant improvements to accuracy over existing methods while being quick to query.