GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration

Jacob R. GardnerGeoff PleissDavid BindelKilian Q. WeinbergerAndrew Gordon Wilson

article2018NeurIPS1,505 citations

Introduces GPyTorch, a PyTorch framework that reduces exact Gaussian process inference complexity from cubic to quadratic time on GPUs using preconditioned blackbox matrix-matrix multiplication.

Listen

Gaussian processes are powerful, flexible machine learning models that provide principled uncertainty estimates, yet their computational demands have historically limited their practical adoption on large datasets. Conventional exact inference methods rely heavily on the Cholesky matrix decomposition, which scales cubically with sample size and fails to fully exploit the parallel compute capabilities of modern graphics processing units (GPUs). Furthermore, existing Gaussian process software packages often tightly couple model definitions with specialized inference algorithms, making it difficult and labor-intensive to prototype advanced architectures or scale to larger problems.

The article evaluates a unified framework called Blackbox Matrix-Matrix (BBMM) inference, introduced alongside GPyTorch, a dedicated software platform built on PyTorch. BBMM reframes the core computational bottlenecks of Gaussian processes into parallel matrix multiplications, demonstrating how exact inference and popular scalable approximations can achieve significant speedups on modern hardware.

To evaluate this framework, the authors implemented a modified batched conjugate gradients algorithm that simultaneously computes model predictions, loss estimates, and parameter gradients within a single unified call, while employing a low-rank pivoted Cholesky preconditioner to accelerate numerical convergence. They benchmarked this approach across exact models and approximation schemes using multiple standard benchmark datasets with sizes ranging from hundreds up to over 500,000 data points.

The analysis yields several key findings. First, BBMM lowers the computational complexity of exact inference from cubic to quadratic time, achieving up to 20 to 32 times faster execution on GPUs compared to traditional CPU-based inference and roughly 4 to 8 times faster performance than GPU-accelerated Cholesky baselines. Second, for scalable approximation techniques like structured kernel interpolation, BBMM accelerates computation by up to 15 times over existing iterative methods. Third, the method matches or slightly improves upon the predictive accuracy of exact models by avoiding the numerical instabilities and artificial noise adjustments frequently required by Cholesky solvers. Finally, using a low-rank preconditioner dramatically accelerates solver convergence with negligible computational overhead.

These results demonstrate that organizations can train Gaussian process models on significantly larger datasets in less time and at lower compute costs without sacrificing mathematical exactness or predictive accuracy. By requiring only standard matrix multiplication routines, the framework decouples model design from inference, allowing engineering and research teams to implement complex or structured approximations in under 50 lines of code.

The authors recommend that teams deploying Gaussian process workflows adopt the open-source GPyTorch library and default to using the pivoted Cholesky preconditioner. Next steps should focus on extending theoretical convergence guarantees beyond one-dimensional cases to multivariate and non-standard kernels, as well as applying the framework more broadly to variational classification tasks. While the results provide high confidence for regression tasks up to GPU memory limits, users applying Gaussian processes to non-Gaussian likelihoods or extreme-scale problems should evaluate specific variational approximations suited to their target domain.

Cornell University gaussian processes gaussians gps GPyTorch blackbox matrix-matrix Gaussian process inference Gaussian process inference

arXiv: 1809.11165

Paper: Variational Learning of Inducing Variables in Sparse Gaussian Processes, Michalis K. Titsias (2009). Titsias establishes the foundational variational framework for inducing variables in sparse Gaussian processes, which GPyTorch incorporates and accelerates using structured matrix-vector and matrix-matrix multiplications.
Paper: A Unifying View of Sparse Approximate Gaussian Process Regression, Joaquin Quiñonero-Candela et al. (2005). Quiñonero-Candela and Rasmussen provide the unified theoretical formulation of sparse Gaussian process approximations whose computational bottlenecks and kernel structures GPyTorch is designed to generalize and scale.
Paper: Sparse Gaussian Processes using Pseudo-inputs, Edward Snelson et al. (2005). Snelson and Ghahramani introduce continuous pseudo-inputs for sparse Gaussian processes, establishing the core inducing-point parameterization that scalable frameworks like GPyTorch optimize.
Paper: Using the Nyström Method to Speed Up Kernel Machines, Christopher K. I. Williams et al. (2000). Williams and Seeger demonstrate low-rank Nyström approximations for large-scale kernel systems, providing essential background on matrix approximation methods relevant to GP scaling.
Paper: Random Features for Large-Scale Kernel Machines, Ali Rahimi et al. (2007). Rahimi and Recht develop random Fourier feature mappings for scaling kernel machines, offering a classic baseline approach to scalable kernel computations.
Paper: PyTorch: An Imperative Style, High-Performance Deep Learning Library, Adam Paszke et al. (2019). Paszke et al. present the PyTorch architecture and automatic differentiation engine upon which the GPyTorch library and its GPU-accelerated operations are built.

Paper: A Tutorial on Bayesian Optimization, Peter I. Frazier (2018). Frazier outlines the principles of Bayesian optimization, providing a direct downstream application area that routinely employs scalable GP models implemented in frameworks like GPyTorch.
Paper: Triton: an intermediate language and compiler for tiled neural network computations, Philippe Tillet et al. (2019). Tillet et al. present Triton for authoring custom high-performance GPU kernels, offering advanced compilation techniques to further accelerate structured matrix-matrix operations in PyTorch-based frameworks.
Paper: PyTorch Distributed: Experiences on Accelerating Data Parallel Training, Shen Li et al. (2020). Li et al. detail PyTorch's distributed data-parallel system, showing how GPU-accelerated PyTorch models can scale out across multi-node and multi-GPU clusters.

Cover for GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration

Abstract

Despite advances in scalable models, the inference tools used for Gaussian processes (GPs) have yet to fully capitalize on developments in computing hardware. We present an efficient and general approach to GP inference based on Blackbox Matrix-Matrix multiplication (BBMM). BBMM inference uses a modified batched version of the conjugate gradients algorithm to derive all terms for training and inference in a single call. BBMM reduces the asymptotic complexity of exact GP inference from $O(n^3)$ to $O(n^2)$ . Adapting this algorithm to scalable approximations and complex GP models simply requires a routine for efficient matrix-matrix multiplication with the kernel and its derivative. In addition, BBMM uses a specialized preconditioner to substantially speed up convergence. In experiments we show that BBMM effectively uses GPU hardware to dramatically accelerate both exact GP inference and scalable approximations. Additionally, we provide GPyTorch, a software platform for scalable GP inference via BBMM, built on PyTorch.

Topics

gaussian processes gaussians gps inference multiplication preconditioners probabilistic inference gradient algorithm gpus

Keywords

GPyTorch blackbox matrix-matrix Gaussian process inference Gaussian process inference GPU acceleration structured kernel interpolation sparse Gaussian processes trace estimation numerical linear algebra

1 Introduction
2 Related Work
3 Background
4 Gaussian process inference through blackbox matrix multiplication
4.1 Preconditioning
5 Programmability with BBMM
6 Results
7 Discussion
References
A Analysis of the modified CG algorithm, mBCG
A.1 Adaptation to multiple right hand sides.
A.2 Obtaining Lanczos tridiagonal matrices from mBCG.
B Runtime analysis of computing inference terms with mBCG
C The Pivoted Cholesky Decomposition
C.1 Running time of the pivoted Cholesky decomposition.
D Convergence Analysis of Pivoted Cholesky Preconditioned CG
E Proofs of Lemmas
E.1 Proof of
E.2 Proof of

Knowls

Knowl 1 — Blackbox Matrix-Matrix (BBMM) Gaussian Process Inference Framework
model/method
Blackbox Matrix-Matrix (BBMM) inference is a general, hardware-accelerated framework for Gaussian process (GP) training and prediction that reduces all core computational steps to matrix-matrix multiplications (MMMs). Given training inputs $X = [x_1, \dots, x_n]^\top \in \mathbb{R}^{n \times d}$ , targets $y \in \mathbb{R}^n$ , kernel covariance matrix $K_{XX} \in \mathbb{R}^{n \times n}$ , and noise variance $\sigma^2$ defining $\tilde{K}_{XX} = K_{XX} + \sigma^2 I_n$ , GP inference requires three operations:
1. The linear solve $\tilde{K}_{XX}^{-1} y$ for the predictive mean and marginal likelihood data-fit term.
2. The log determinant $\log |\tilde{K}_{XX}|$ for the marginal log likelihood complexity penalty.
3. The derivative trace term $\operatorname{Tr}\left(\tilde{K}_{XX}^{-1} \frac{d\tilde{K}_{XX}}{d\theta}\right)$ for hyperparameter optimization with respect to $\theta$ .
BBMM derives all three quantities simultaneously in a single call to a modified Batched Conjugate Gradients (mBCG) routine. Rather than computing an explicit $O(n^3)$ Cholesky factorization, BBMM requires only blackbox functions that compute matrix-matrix multiplications with the kernel matrix $\tilde{K}_{XX} M$ and its derivative $\frac{d\tilde{K}_{XX}}{d\theta} M$ for an $n \times t$ matrix $M$ . For exact GPs, BBMM reduces the asymptotic time complexity from $O(n^3)$ to $O(p n^2 t)$ for $p$ iterations and $t$ probe vectors, and reduces auxiliary space complexity from $O(n^2)$ to $O(nt)$ .
Knowl 2 — Modified Batched Conjugate Gradients (mBCG)
algorithm
The modified Batched Conjugate Gradients (mBCG) algorithm solves multiple linear systems simultaneously for an $n \times (t+1)$ matrix $B = [y, z_1, \dots, z_t]$ using matrix-matrix products while concurrently computing $t$ partial Lanczos tridiagonalization matrices $\tilde{T}_1, \dots, \tilde{T}_t \in \mathbb{R}^{p \times p}$ directly from the scalar recurrence coefficients generated during conjugate gradient iterations.
Input: mmm_A() - function computing matrix-matrix product with matrix $A \in \mathbb{R}^{n \times n}$
        $B \in \mathbb{R}^{n \times (t+1)}$ - matrix of right-hand sides $[y, z_1, \dots, z_t]$
       P_inv() - function computing preconditioner solve $P^{-1} R$
        $p$ - maximum number of iterations
        $\text{tolerance}$ - residual convergence threshold
Output: $U_p \approx A^{-1} B$ , tridiagonal matrices $\tilde{T}_1, \dots, \tilde{T}_t \in \mathbb{R}^{p \times p}$
$U_0 \leftarrow 0$
$R_0 \leftarrow \text{mmm\_A}(U_0) - B$
$Z_0 \leftarrow \text{P\_inv}(R_0)$
$D_0 \leftarrow Z_0$
$\tilde{T}_1, \dots, \tilde{T}_t \leftarrow 0$
for $j \leftarrow 1$ to $p$ do
     $V_j \leftarrow \text{mmm\_A}(D_{j-1})$
     $\vec{\alpha}_j \leftarrow ((R_{j-1} \circ Z_{j-1})^\top \mathbf{1}) \oslash ((D_{j-1} \circ V_j)^\top \mathbf{1})$
     $U_j \leftarrow U_{j-1} + D_{j-1} \operatorname{diag}(\vec{\alpha}_j)$
     $R_j \leftarrow R_{j-1} - V_j \operatorname{diag}(\vec{\alpha}_j)$
    if $\max_{1 \le i \le t+1} \|[R_j]_{:, i}\|_2 < \text{tolerance}$ then
        return $U_j, \tilde{T}_1, \dots, \tilde{T}_t$
     $Z_j \leftarrow \text{P\_inv}(R_j)$
     $\vec{\beta}_j \leftarrow ((Z_j \circ Z_j)^\top \mathbf{1}) \oslash ((Z_{j-1} \circ Z_{j-1})^\top \mathbf{1})$
     $D_j \leftarrow Z_j - D_{j-1} \operatorname{diag}(\vec{\beta}_j)$
    for $i \leftarrow 1$ to $t$ do
         $[\tilde{T}_i]_{j, j} \leftarrow \frac{1}{[\vec{\alpha}_j]_{i+1}} + \frac{[\vec{\beta}_{j-1}]_{i+1}}{[\vec{\alpha}_{j-1}]_{i+1}}$ (with $[\vec{\beta}_0]_{i+1} = 0, [\vec{\alpha}_0]_{i+1} = \infty$ )
        if $j > 1$ then
             $[\tilde{T}_i]_{j-1, j} \leftarrow \frac{\sqrt{[\vec{\beta}_{j-1}]_{i+1}}}{[\vec{\alpha}_{j-1}]_{i+1}}$
             $[\tilde{T}_i]_{j, j-1} \leftarrow \frac{\sqrt{[\vec{\beta}_{j-1}]_{i+1}}}{[\vec{\alpha}_{j-1}]_{i+1}}$
end
return $U_p, \tilde{T}_1, \dots, \tilde{T}_t$
In the algorithm, $\circ$ denotes the elementwise Hadamard product, $\oslash$ denotes elementwise division, and $\mathbf{1} \in \mathbb{R}^n$ is the all-ones vector. Running $p$ iterations of mBCG takes $O(p \, \Xi(A) + pnt)$ time and $O(nt + p^2 t)$ space, where $\Xi(A)$ represents the computational cost of multiplying $A$ by an $n \times (t+1)$ matrix.
Knowl 3 — Estimating Marginal Log-Likelihood Derivatives and Log Determinants from mBCG
model/method
Using the outputs of mBCG $[u_0, u_1, \dots, u_t] = \tilde{K}_{XX}^{-1} [y, z_1, \dots, z_t]$ and the partial Lanczos tridiagonal matrices $\tilde{T}_1, \dots, \tilde{T}_t \in \mathbb{R}^{p \times p}$ obtained from random probe vectors $z_1, \dots, z_t \overset{\text{iid}}{\sim} \mathcal{N}(0, I_n)$ , Gaussian process training and prediction terms are estimated as follows:
1. Linear Solves and Predictions: The solve $\tilde{K}_{XX}^{-1} y$ equals $u_0$ . Given test input $x^*$ and cross-covariance vector $k_{Xx^*}$ , the predictive posterior mean and covariance are: $\mu_{f \mid \mathcal{D}}(x^*) = \mu(x^*) + k_{Xx^*}^\top u_0, \quad k_{f \mid \mathcal{D}}(x^*, x^{*\prime}) = k(x^*, x^{*\prime}) - k_{Xx^*}^\top \tilde{K}_{XX}^{-1} k_{Xx^{*\prime}}$
2. Log Determinant Estimation via Stochastic Lanczos Quadrature: Let $\tilde{T}_i = V_i \Lambda_i V_i^\top$ be the eigendecomposition of the $p \times p$ tridiagonal matrix $\tilde{T}_i$ , and let $e_1 = [1, 0, \dots, 0]^\top \in \mathbb{R}^p$ . The log determinant is estimated by: $\log |\tilde{K}_{XX}| \approx \frac{1}{t} \sum_{i=1}^t e_1^\top (\log \tilde{T}_i) e_1 = \frac{1}{t} \sum_{i=1}^t e_1^\top V_i (\log \Lambda_i) V_i^\top e_1$ Eigendecomposing each tridiagonal matrix requires $O(p^2)$ time, costing $O(tp^2)$ post-processing work total.
3. Marginal Log-Likelihood Derivative via Stochastic Trace Estimation: For negative log marginal likelihood $\mathcal{L}(\theta \mid X, y) \propto \log |\tilde{K}_{XX}| - y^\top \tilde{K}_{XX}^{-1} y$ , the trace component of the derivative $\frac{d\mathcal{L}}{d\theta} = y^\top \tilde{K}_{XX}^{-1} \frac{d\tilde{K}_{XX}}{d\theta} \tilde{K}_{XX}^{-1} y + \operatorname{Tr}\left(\tilde{K}_{XX}^{-1} \frac{d\tilde{K}_{XX}}{d\theta}\right)$ is unbiasedly estimated as: $\operatorname{Tr}\left(\tilde{K}_{XX}^{-1} \frac{d\tilde{K}_{XX}}{d\theta}\right) \approx \frac{1}{t} \sum_{i=1}^t u_i^\top \left( \frac{d\tilde{K}_{XX}}{d\theta} z_i \right)$ This requires a single matrix multiply $\frac{d\tilde{K}_{XX}}{d\theta} [z_1, \dots, z_t]$ and $t$ inner products costing $O(nt)$ additional operations.
Knowl 4 — Pivoted Cholesky Preconditioner for Gaussian Process Covariance Matrices
model/method
To accelerate mBCG convergence, BBMM preconditions $\tilde{K}_{XX} = K_{XX} + \sigma^2 I_n$ using a rank- $k$ pivoted Cholesky decomposition of the kernel matrix, $K_{XX} \approx L_k L_k^\top$ with $L_k \in \mathbb{R}^{n \times k}$ and $k \ll n$ , forming the positive definite preconditioner: $P_k = L_k L_k^\top + \sigma^2 I_n$

$P_k$ provides four essential computational properties for BBMM:
1. Efficient Construction: Computing $L_k$ requires accessing the diagonal of $K_{XX}$ and greedily selecting $k$ rows based on maximal diagonal Schur complement entries. This takes $O(\rho(K_{XX}) k^2)$ time, where $\rho(K_{XX})$ is the time required to evaluate a single row of $K_{XX}$ ( $O(n k^2)$ for standard dense kernels).
2. Linear Solves in $O(n k^2)$ Time: By the Woodbury matrix inversion identity: $P_k^{-1} v = \frac{1}{\sigma^2} v - \frac{1}{\sigma^4} L_k \left(I_k + \frac{1}{\sigma^2} L_k^\top L_k\right)^{-1} L_k^\top v$ Computing $L_k^\top v$ takes $O(nk)$ operations, forming and factorizing the $k \times k$ matrix $(I_k + \sigma^{-2} L_k^\top L_k)$ takes $O(nk^2 + k^3)$ operations, and back-solving takes $O(nk)$ operations.
3. Exact Log Determinant in $O(n k^2)$ Time: By the matrix determinant lemma: $\log |P_k| = \log \left| I_k + \frac{1}{\sigma^2} L_k^\top L_k \right| + 2n \log \sigma$ which is computed directly from the Cholesky factor of the $k \times k$ matrix in $O(nk^2 + k^3)$ time.
4. Exact Sampling in $O(nk)$ Time: Sampling $z \sim \mathcal{N}(0, P_k)$ is performed via the reparameterization trick using standard Gaussian vectors $\epsilon_1 \sim \mathcal{N}(0, I_k)$ and $\epsilon_2 \sim \mathcal{N}(0, I_n)$ : $z = L_k \epsilon_1 + \sigma \epsilon_2$
Knowl 5 — Preconditioned Stochastic Estimators for Log Determinants and Trace Derivatives
model/method
When preconditioning mBCG with symmetric positive definite preconditioner $P$ , probe vectors are drawn from $z_i \sim \mathcal{N}(0, P)$ rather than $\mathcal{N}(0, I_n)$ . Because $P^{-1/2} z_i \sim \mathcal{N}(0, I_n)$ , the preconditioned system allows unbiased recovery of the unpreconditioned log determinant and trace derivative:
1. Preconditioned Log Determinant Estimator: mBCG preconditioned with $P$ constructs Lanczos tridiagonal matrices $T_i \in \mathbb{R}^{p \times p}$ that approximate the spectrum of the preconditioned matrix $P^{-1/2} \tilde{K}_{XX} P^{-1/2}$ . The log determinant of $\tilde{K}_{XX}$ is recovered via the identity $\log |\tilde{K}_{XX}| = \log |P^{-1/2} \tilde{K}_{XX} P^{-1/2}| + \log |P|$ as: $\log |\tilde{K}_{XX}| \approx \left( \frac{1}{t} \sum_{i=1}^t e_1^\top (\log T_i) e_1 \right) + \log |P|$ where $e_1 = [1, 0, \dots, 0]^\top \in \mathbb{R}^p$ . As the preconditioner rank $k \to n$ , $P^{-1/2} \tilde{K}_{XX} P^{-1/2} \to I_n$ , $\log |P^{-1/2} \tilde{K}_{XX} P^{-1/2}| \to 0$ , and the estimate converges to the exact value $\log |P|$ .
2. Preconditioned Trace Derivative Estimator: Using the preconditioned solves $u_i = \tilde{K}_{XX}^{-1} z_i$ returned by mBCG for $z_i \sim \mathcal{N}(0, P)$ , the derivative trace term is estimated by applying the derivative to the preconditioned vectors $P^{-1} z_i$ : $\operatorname{Tr}\left(\tilde{K}_{XX}^{-1} \frac{d\tilde{K}_{XX}}{d\theta}\right) = \operatorname{Tr}\left(\tilde{K}_{XX}^{-1} \frac{d\tilde{K}_{XX}}{d\theta} \mathbb{E}_{z_i \sim \mathcal{N}(0, P)}[P^{-1} z_i z_i^\top]\right) \approx \frac{1}{t} \sum_{i=1}^t u_i^\top \left( \frac{d\tilde{K}_{XX}}{d\theta} P^{-1} z_i \right)$
Knowl 6 — Condition Number Bound for Pivoted Cholesky Preconditioned RBF Kernels
theoretical result
Let $K_{XX} \in \mathbb{R}^{n \times n}$ be a univariate RBF kernel matrix, let $\tilde{K}_{XX} = K_{XX} + \sigma^2 I_n$ , let $L_k L_k^\top$ be the rank- $k$ pivoted Cholesky decomposition of $K_{XX}$ , and let $P_k = L_k L_k^\top + \sigma^2 I_n$ .

There exists a constant $b > 0$ such that the condition number $\kappa$ of the preconditioned matrix satisfies: $\kappa\left(P_k^{-1/2} \tilde{K}_{XX} P_k^{-1/2}\right) = \kappa\left(P_k^{-1} \tilde{K}_{XX}\right) \le \|P_k^{-1} \tilde{K}_{XX}\|_2 \|\tilde{K}_{XX}^{-1} P_k\|_2 \le \left(1 + \mathcal{O}\left(n \exp(-bk)\right)\right)^2$

This bound establishes that the condition number of the preconditioned system decays exponentially toward $1$ as the rank $k$ of the pivoted Cholesky preconditioner increases.
Knowl 7 — Convergence Rate of Pivoted Cholesky Preconditioned Conjugate Gradients
theoretical result
Let $K_{XX} \in \mathbb{R}^{n \times n}$ be a univariate RBF kernel matrix, $\tilde{K}_{XX} = K_{XX} + \sigma^2 I_n$ , and let $L_k L_k^\top$ be the rank- $k$ pivoted Cholesky decomposition of $K_{XX}$ forming preconditioner $P_k = L_k L_k^\top + \sigma^2 I_n$ .

When preconditioned conjugate gradients is used to solve the linear system $\tilde{K}_{XX} u = y$ , the error of the $p$ -th iterate $u_p$ in the $\tilde{K}_{XX}$ -norm relative to the exact solution $u^* = \tilde{K}_{XX}^{-1} y$ is bounded by: $\|u^* - u_p\|_{\tilde{K}_{XX}} \le 2 \left(\frac{1}{1 + \mathcal{O}(\exp(kb)/n)}\right)^p \|u^* - u_0\|_{\tilde{K}_{XX}}$ for some constant $b > 0$ , where $\|v\|_{\tilde{K}_{XX}} = (v^\top \tilde{K}_{XX} v)^{1/2}$ and $u_0$ is the initial iterate.

This result demonstrates that the convergence rate of conjugate gradient linear solves accelerates exponentially as the rank $k$ of the pivoted Cholesky preconditioner is increased.
Knowl 8 — Eigenvalue Bound for Univariate RBF Kernel Matrices
theoretical result
Let $x_1, \dots, x_n \in [0, 1]$ and let $K_{XX} \in \mathbb{R}^{n \times n}$ be the univariate RBF kernel matrix with entries $K_{ij} = \exp(-\gamma (x_i - x_j)^2)$ for lengthscale parameter $\gamma > 0$ .

The $(2l+1)$ -th eigenvalue $\lambda_{2l+1}(K_{XX})$ (ordered in non-increasing magnitude) satisfies the upper bound: $\lambda_{2l+1}(K_{XX}) \le 2n e^{-\gamma/4} I_{l+1}(\gamma/4) \sim \frac{2n e^{-\gamma/4}}{\sqrt{\pi \gamma}} \left( \frac{e\gamma}{8(l+1)} \right)^{l+1}$ where $I_{l+1}$ denotes the modified Bessel function of the first kind with order $l+1$ .

This bound establishes that the eigenvalues of univariate RBF kernel matrices decay super-exponentially with index $l$ .
Knowl 9 — Computational Complexity of BBMM for Structured and Scalable GP Models
model/method
BBMM requires only a routine for matrix-matrix multiplication (MMM) with the kernel $\tilde{K}_{XX} M$ and its derivative $\frac{d\tilde{K}_{XX}}{d\theta} M$ for an $n \times t$ matrix $M$ , adapting directly to various structured GP families with the following time complexities for $p$ iterations of mBCG:
1. Exact Gaussian Process: $\tilde{K}_{XX} = K_{XX} + \sigma^2 I_n$ . Dense kernel multiplication requires $O(n^2 t)$ per iteration, yielding total runtime $O(p n^2 t)$ and space $O(nt)$ , compared to $O(n^3)$ time and $O(n^2)$ space for Cholesky factorization.
2. Bayesian Linear Regression: $\tilde{K}_{XX} = X X^\top + \sigma^2 I_n$ for $X \in \mathbb{R}^{n \times d}$ . Associating $(X (X^\top M) + \sigma^2 M)$ requires $O(tnd)$ operations per iteration, giving $O(ptnd)$ total time and achieving exact inference in $O(tnd^2)$ time without bespoke derivations.
3. Sparse Gaussian Process Regression (SGPR) / Subset of Regressors (SoR): $\tilde{K}_{XX} \approx K_{XU} K_{UU}^{-1} K_{UX} + \sigma^2 I_n$ with $m$ inducing points ( $m \ll n$ ). Computing $K_{XU} (K_{UU}^{-1} (K_{UX} M)) + \sigma^2 M$ requires $O(tnm + tm^3)$ time per iteration, which is asymptotically faster than the $O(nm^2 + m^3)$ time required by Cholesky-based SGPR inference engines.
4. Structured Kernel Interpolation (SKI / KISS-GP): $\tilde{K}_{XX} \approx W K_{UU} W^\top + \sigma^2 I_n$ with sparse interpolation weight matrix $W$ and Toeplitz or Kronecker structured $K_{UU}$ . Toeplitz matrix multiplies take $O(m \log m)$ , yielding $O(tn + tm \log m)$ per iteration.
5. Kernel Compositions: Additive kernels $(K_1 + K_2) M = K_1 M + K_2 M$ and product kernels execute compositionally via primitive matrix-matrix multiplications.
Knowl 10 — Empirical Speedup and Accuracy of GPU BBMM vs Cholesky and Lanczos Baselines
empirical result
Experiments evaluating GPU-accelerated BBMM (using $p = 20$ iterations of CG, $t = 10$ Rademacher probe vectors, and rank $k = 5$ pivoted Cholesky preconditioners on an NVIDIA Titan Xp GPU against CPU and GPU baselines with an Intel Xeon E5-2650 CPU) show:
1. Inference Speedups:
- Exact GPs ( $n \le 3500$ on UCI datasets: Skillcraft, Gas, Wine, Airfoil, Autompg): GPU BBMM achieves up to $32\times$ speedup over CPU Cholesky (GPFlow) and up to $20\times$ speedup over GPU Cholesky, because Cholesky's triangular solve operations parallelize poorly on GPUs (GPU Cholesky only achieves $\approx 4\times$ speedup over CPU).
- SGPR ( $n$ up to $50,000$ , $m = 300$ inducing points on KEGG, Protein, Kin40k, Elevators, PoleTele): GPU BBMM achieves up to $10\times$ speedup over CPU Cholesky and $4\times$ over GPU Cholesky.
- SKI + Deep Kernel Learning ( $n$ up to $515,000$ , $m = 10,000$ on Song, Buzz, KEGG, Protein, Kin40k): GPU BBMM achieves up to $32\times$ speedup over CPU Lanczos/CG and up to $15\times$ speedup over GPU Lanczos/CG (Dong et al.), because mBCG computes all probe solves in parallel rather than sequentially.
1. Accuracy and Regularization: Across all benchmark datasets, BBMM achieves test mean absolute error (MAE) comparable to or lower than Cholesky-based inference for both RBF and Matérn-5/2 kernels. In single precision, conjugate gradients exhibits a regularizing effect that avoids numerical instabilities without requiring the artificial diagonal jitter or double precision computation necessary in Cholesky factorization.
2. Preconditioning Efficiency: Pivoted Cholesky preconditioners with rank $k \in \{2, 5, 9\}$ substantially decrease the relative linear solve residual error $\|\tilde{K}_{XX} u - y\|_2 / \|y\|_2$ within 20 iterations, achieving lower test MAE in fewer wall-clock seconds than unpreconditioned conjugate gradients while adding negligible per-iteration time overhead.

Coverage note — Intermediate mathematical lemmata (Lemmas 4, 5, and 6) supporting the RBF eigenvalue decay proof, Theorem 2 restated from Ubaru et al. (2017), and the brief discussion on non-Gaussian variational ELBO extensions were deliberately omitted in accordance with the rules on non-redundancy and standalone significance.

References

1.M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al. Tensorflow: A system for large-scale machine learning. In OSDI, volume 16, pages 265–283, 2016.
2.A. Asuncion and D. Newman. Uci machine learning repository. https://archive.ics.uci.edu/ml/, 2007. Last accessed: 2018-05-18.
3.H. Avron and S. Toledo. Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. Journal of the ACM (JACM), 58(2):8, 2011.
4.F. Bach. Sharp analysis of low-rank kernel matrix approximations. In COLT, 2013.
5.E. V. Bonilla, K. M. Chai, and C. Williams. Multi-task Gaussian process prediction. In NIPS, 2008.
6.L. Bottou. Large-scale machine learning with stochastic gradient descent. In COMPSTAT, pages 177–186. Springer, 2010.
7.P. Chaudhari, A. Choromanska, S. Soatto, Y. LeCun, C. Baldassi, C. Borgs, J. Chayes, L. Sagun, and R. Zecchina. Entropy-sgd: Biasing gradient descent into wide valleys. arXiv preprint arXiv:1611.01838, 2016.
8.T. Chen, M. Li, Y. Li, M. Lin, N. Wang, M. Wang, T. Xiao, B. Xu, C. Zhang, and Z. Zhang. Mxnet: A flexible and efficient machine learning library for heterogeneous distributed systems. arXiv preprint arXiv:1512.01274, 2015.
9.J. P. Cunningham, K. V. Shenoy, and M. Sahani. Fast Gaussian process methods for point process intensity estimation. In ICML, 2008.
10.K. Cutajar, M. Osborne, J. Cunningham, and M. Filippone. Preconditioning kernel matrices. In ICML, 2016.
11.B. N. Datta. Numerical linear algebra and applications, volume 116. Siam, 2010.
12.J. W. Demmel. Applied numerical linear algebra, volume 56. Siam, 1997.
13.K. Dong, D. Eriksson, H. Nickisch, D. Bindel, and A. G. Wilson. Scalable log determinants for Gaussian process kernel learning. In NIPS, 2017.
14.J. K. Fitzsimons, M. A. Osborne, S. J. Roberts, and J. F. Fitzsimons. Improved stochastic trace estimation using mutually unbiased bases. arXiv preprint arXiv:1608.00117, 2016.
15.J. R. Gardner, G. Pleiss, R. Wu, K. Q. Weinberger, and A. G. Wilson. Product kernel interpolation for scalable Gaussian processes. In AISTATS, 2018.
16.G. H. Golub and G. Meurant. Matrices, moments and quadrature with applications. Princeton University Press, 2009.
17.G. H. Golub and C. F. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
18.R. H. Hahnloser, R. Sarpeshkar, M. A. Mahowald, R. J. Douglas, and H. S. Seung. Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit. Nature, 405(6789):947, 2000.
19.H. Harbrecht, M. Peters, and R. Schneider. On the low-rank approximation by the pivoted cholesky decomposition. Applied numerical mathematics, 62(4):428–440, 2012.
20.K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In CVPR, 2016.
21.J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In UAI, 2013.
22.J. Hensman, A. G. d. G. Matthews, and Z. Ghahramani. Scalable variational Gaussian process classification. In ICML, 2015.
23.S. Hochreiter and J. Schmidhuber. Flat minima. Neural Computation, 9(1):1–42, 1997.
24.G. Huang, Z. Liu, K. Q. Weinberger, and L. van der Maaten. Densely connected convolutional networks. In CVPR, 2017.
25.M. F. Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statistics-Simulation and Computation, 19(2):433–450, 1990.
26.S. Ioffe and C. Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In ICML, 2015.
27.P. Izmailov, D. Podoprikhin, T. Garipov, D. Vetrov, and A. G. Wilson. Averaging weights leads to wider optima and better generalization. In Uncertainty in Artificial Intelligence (UAI), 2018.
28.Y. Jia, E. Shelhamer, J. Donahue, S. Karayev, J. Long, R. Girshick, S. Guadarrama, and T. Darrell. Caffe: Convolutional architecture for fast feature embedding. In ACMMM, pages 675–678. ACM, 2014.
29.N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, and P. T. P. Tang. On large-batch training for deep learning: Generalization gap and sharp minima. arXiv preprint arXiv:1609.04836, 2016.
30.D. P. Kingma and M. Welling. Auto-encoding variational Bayes. In ICLR, 2014.
31.A. Krizhevsky, I. Sutskever, and G. E. Hinton. Imagenet classification with deep convolutional neural networks. In NIPS, 2012.
32.A. Krizhevsky, I. Sutskever, and G. E. Hinton. Imagenet classification with deep convolutional neural networks. In Advances in neural information processing systems, pages 1097–1105, 2012.
33.C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. United States Governm. Press Office Los Angeles, CA, 1950.
34.A. G. d. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. León-Villagrá, Z. Ghahramani, and J. Hensman. Gpflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40):1–6, 2017.
35.I. Murray. Gaussian processes and fast matrix-vector multiplies. In ICML Workshop on Numerical Mathematics in Machine Learning, 2009.
36.D. P. O’Leary. The block conjugate gradient algorithm and related methods. Linear algebra and its applications, 29:293–322, 1980.
37.C. Paige. Practical use of the symmetric Lanczos process with re-orthogonalization. BIT Numerical Mathematics, 10(2):183–195, 1970.
38.B. N. Parlett. A new look at the Lanczos algorithm for solving symmetric systems of linear equations. Linear algebra and its applications, 29:323–346, 1980.
39.A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer. Automatic differentiation in PyTorch. 2017.
40.G. Pleiss, J. R. Gardner, K. Q. Weinberger, and A. G. Wilson. Constant-time predictive distributions for Gaussian processes. In ICML, 2018.
41.J. Quiñonero-Candela and C. E. Rasmussen. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6(Dec):1939–1959, 2005.
42.C. E. Rasmussen and C. K. Williams. Gaussian processes for machine learning, volume 1. MIT press Cambridge, 2006.
43.Y. Saad. Iterative methods for sparse linear systems, volume 82. siam, 2003.
44.Y. Saatçi. Scalable inference for structured Gaussian process models. PhD thesis, University of Cambridge, 2012.
45.E. Snelson and Z. Ghahramani. Sparse Gaussian processes using pseudo-inputs. In NIPS, 2006.
46.M. K. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In AISTATS, pages 567–574, 2009.
47.S. Ubaru, J. Chen, and Y. Saad. Fast estimation of tr (f (a)) via stochastic Lanczos quadrature. SIAM Journal on Matrix Analysis and Applications, 38(4):1075–1099, 2017.
48.H. A. Van der Vorst. Iterative Krylov methods for large linear systems, volume 13. Cambridge University Press, 2003.
49.A. J. Wathen and S. Zhu. On spectral distribution of kernel matrices related to radial basis functions. Numerical Algorithms, 70(4):709–726, 2015.
50.A. G. Wilson. Covariance kernels for fast automatic pattern discovery and extrapolation with Gaussian processes. PhD thesis, University of Cambridge, 2014.
51.A. G. Wilson and H. Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In ICML, 2015.
52.A. G. Wilson, C. Dann, and H. Nickisch. Thoughts on massively scalable Gaussian processes. arXiv preprint arXiv:1511.01870, 2015.
53.A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing. Deep kernel learning. In AISTATS, 2016.
54.A. G. Wilson, Z. Hu, R. R. Salakhutdinov, and E. P. Xing. Stochastic variational deep kernel learning. In NIPS, 2016.

Citation

MLA

Gardner, J. R., et al. “GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration”. arXiv, 2018, http://arxiv.org/abs/1809.11165v6.

APA

Gardner, J. R., Pleiss, G., Bindel, D., Weinberger, K. Q., & Wilson, A. G. (2018). GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. arXiv. http://arxiv.org/abs/1809.11165v6

Chicago

Gardner, J. R., G. Pleiss, D. Bindel, K. Q. Weinberger, and A. G. Wilson. 2018. “GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration”. arXiv. http://arxiv.org/abs/1809.11165v6.

Harvard

Gardner, J.R. et al. (2018) “GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration”, arXiv [Preprint]. Available at: http://arxiv.org/abs/1809.11165v6.

Vancouver

1. Gardner JR, Pleiss G, Bindel D, Weinberger KQ, Wilson AG (2018) GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. arXiv

BibTeX

@article{gardner2018gpytorch,
  title = {GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration},
  author = {Gardner, Jacob R. and Pleiss, Geoff and Bindel, David and Weinberger, Kilian Q. and Wilson, Andrew Gordon},
  year = {2018},
  journal = {arXiv},
  url = {http://arxiv.org/abs/1809.11165v6},
  eprint = {1809.11165}
}

Metadata:arXiv

Access the Paper

This paper is available from its original source. Click below to access the PDF.

Open PDF

License: Published with permission

GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration

Abstract

Topics

Keywords

Table of Contents

Knowls

Knowl 1 — Blackbox Matrix-Matrix (BBMM) Gaussian Process Inference Framework

Knowl 2 — Modified Batched Conjugate Gradients (mBCG)

Knowl 3 — Estimating Marginal Log-Likelihood Derivatives and Log Determinants from mBCG

Knowl 4 — Pivoted Cholesky Preconditioner for Gaussian Process Covariance Matrices

Knowl 5 — Preconditioned Stochastic Estimators for Log Determinants and Trace Derivatives

Knowl 6 — Condition Number Bound for Pivoted Cholesky Preconditioned RBF Kernels

Knowl 7 — Convergence Rate of Pivoted Cholesky Preconditioned Conjugate Gradients

Knowl 8 — Eigenvalue Bound for Univariate RBF Kernel Matrices

Knowl 9 — Computational Complexity of BBMM for Structured and Scalable GP Models

Knowl 10 — Empirical Speedup and Accuracy of GPU BBMM vs Cholesky and Lanczos Baselines

References

Citation

BibTeX

Access the Paper