# Fast Matrix Square Roots with Applications to Gaussian Processes and Bayesian Optimization

arXiv: https://arxiv.org/abs/2006.11267 · HTML: https://mlsanity.com/paper/2006.11267/

## Key facts

- Cited by 26 papers in the map (data through 14 Oct 2026).
- First posted on arXiv 19 Jun 2020.

## Built on by (most cited first, 26 of 26)

- [Parallel Bayesian Optimization of Multiple Noisy Objectives with Expected Hypervolume Improvement](https://mlsanity.com/paper/2105.08195/) (2021)
- [Multi-Objective Bayesian Optimization over High-Dimensional Search Spaces](https://mlsanity.com/paper/2109.10964/) (2021)
- [Forward-backward Gaussian variational inference via JKO in the Bures-Wasserstein Space](https://mlsanity.com/paper/2304.05398/) (2023)
- [Bayesian Algorithm Execution: Estimating Computable Properties of Black-box Functions Using Mutual Information](https://mlsanity.com/paper/2104.09460/) (2021)
- [Bayesian Optimization with High-Dimensional Outputs](https://mlsanity.com/paper/2106.12997/) (2021)
- [Scaling Gaussian Processes with Derivative Information Using Variational Inference](https://mlsanity.com/paper/2107.04061/) (2021)
- [Surrogate Likelihoods for Variational Annealed Importance Sampling](https://mlsanity.com/paper/2112.12194/) (2021)
- [When are Iterative Gaussian Processes Reliably Accurate?](https://mlsanity.com/paper/2112.15246/) (2021)
- [Computation-Aware Kalman Filtering and Smoothing](https://mlsanity.com/paper/2405.08971/) (2024)
- [Large-Scale Gaussian Processes via Alternating Projection](https://mlsanity.com/paper/2310.17137/) (2023)
- [Low-Precision Arithmetic for Fast Gaussian Processes](https://mlsanity.com/paper/2207.06856/) (2022)
- [Universal Functional Regression with Neural Operator Flows](https://mlsanity.com/paper/2404.02986/) (2024)
- [Computation-Aware Gaussian Processes: Model Selection And Linear-Time Inference](https://mlsanity.com/paper/2411.01036/) (2024)
- [Variational Nearest Neighbor Gaussian Process](https://mlsanity.com/paper/2202.01694/) (2022)
- [Hierarchical Inducing Point Gaussian Process for Inter-domain Observations](https://mlsanity.com/paper/2103.00393/) (2021)
- [Actually Sparse Variational Gaussian Processes](https://mlsanity.com/paper/2304.05091/) (2023)
- [PRISM: Distribution-free Adaptive Computation of Matrix Functions for Accelerating Neural Network Training](https://mlsanity.com/paper/2601.22137/) (2026)
- [Low-Rank Updates of Matrix Square Roots](https://mlsanity.com/paper/2201.13156/) (2022)
- [Sparse Kernel Gaussian Processes through Iterative Charted Refinement (ICR)](https://mlsanity.com/paper/2206.10634/) (2022)
- [Practical Bayesian Algorithm Execution via Posterior Sampling](https://mlsanity.com/paper/2410.20596/) (2024)
- [Provably Reliable Large-Scale Sampling from Gaussian Processes](https://mlsanity.com/paper/2211.08036/) (2022)
- [Gradients of Functions of Large Matrices](https://mlsanity.com/paper/2405.17277/) (2024)
- [Adaptive Candidate Point Thompson Sampling for High-Dimensional Bayesian Optimization](https://mlsanity.com/paper/2604.08891/) (2026)
- [Search at the Cost of Sampling: Nearly-Instant Latent Space Bayesian Optimization](https://mlsanity.com/paper/2609.19476/) (2026)
- [Variance Reduction for Matrix Computations with Applications to Gaussian Processes](https://mlsanity.com/paper/2106.14565/) (2021)
- [Privacy-aware Gaussian Process Regression](https://mlsanity.com/paper/2305.16541/) (2023)

## Abstract (arXiv)

Matrix square roots and their inverses arise frequently in machine learning, e.g., when sampling from high-dimensional Gaussians N(mathbf 0, mathbf K) or whitening a vector mathbf b against covariance matrix mathbf K. While existing methods typically require O(N3) computation, we introduce a highly-efficient quadratic-time algorithm for computing mathbf K1/2 mathbf b, mathbf K-1/2 mathbf b, and their derivatives through matrix-vector multiplication (MVMs). Our method combines Krylov subspace methods with a rational approximation and typically achieves 4 decimal places of accuracy with fewer than 100 MVMs. Moreover, the backward pass requires little additional computation. We demonstrate our method's applicability on matrices as large as 50,!000 × 50,!000 - well beyond traditional methods - with little approximation error. Applying this increased scalability to variational Gaussian processes, Bayesian optimization, and Gibbs sampling results in more powerful models with higher accuracy.