Large-Scale Gaussian Processes via Alternating Projection

Kaiwen Wu1 Jonathan Wenger2 Haydn Jones1 Geoff Pleiss3,4 Jacob R. Gardner1
1University of Pennsylvania 2Columbia University 3University of British
Columbia 4Vector Institute

###### Abstract

Gaussian process (GP) hyperparameter optimization requires repeatedly solving
linear systems with $n\times n$ kernel matrices. To address the prohibitive
$\mathcal{O}(n^{3})$ time complexity, recent work has employed fast iterative
numerical methods, like conjugate gradients (CG). However, as datasets
increase in magnitude, the corresponding kernel matrices become increasingly
ill-conditioned and still require $\mathcal{O}(n^{2})$ space without
partitioning. Thus, while CG increases the size of datasets GPs can be trained
on, modern datasets reach scales beyond its applicability. In this work, we
propose an iterative method which only accesses subblocks of the kernel
matrix, effectively enabling _mini-batching_. Our algorithm, based on
alternating projection, has $\mathcal{O}(n)$ per-iteration time and space
complexity, solving many of the practical challenges of scaling GPs to very
large datasets. Theoretically, we prove our method enjoys linear convergence
and empirically we demonstrate its robustness to ill-conditioning. On large-
scale benchmark datasets up to four million datapoints our approach
accelerates training by a factor of 2$\times$ to 27$\times$ compared to CG.

##  1 INTRODUCTION

(a) 3droad

(b) house electric

Figure 1:  Comparison of the convergence of alternating projection and
(preconditioned) conjugate gradient. Both algorithms are _initialized at zero_
, but CG increases the residual after the first iteration. Left: While the
asymptotic convergence rate of CG can be faster than alternating projection,
CG does not find a better solution than alternating projection in the first
$1000$ iterations. Right: CG struggles with convergence due to ill-
conditioning and does not reach the tolerance $\delta$. In contrast,
alternating projection convergences. See §4 for more details.

Scaling Gaussian process (GP) models to large datasets has been a central
research topic in probabilistic machine learning for nearly two decades. The
primary challenge is the cubic complexity of computing both the marginal log
likelihood (MLL) during training and the predictive distribution at test time.
Over the years, this problem has been addressed both from a modeling
perspective (_e.g_., Hensman et al.,, 2013, 2015; Titsias,, 2009; Snelson and
Ghahramani,, 2005; Salimbeni et al.,, 2018; Jankowiak et al.,, 2020; Katzfuss
and Guinness,, 2021) and from a numerical methods perspective (_e.g_., Cutajar
et al.,, 2016; Pleiss et al.,, 2018; Gardner et al.,, 2018; Wang et al.,,
2019; Maddox et al.,, 2022), and contemporary work even unifies these
perspectives to a degree (Artemev et al.,, 2021; Wenger et al., 2022b, ). In
recent years, numerical methods have increasingly relied on matrix-free
iterative methods, which access the kernel matrix through matrix-vector
multiplications. These iterations are suitable for GPU acceleration (Gardner
et al.,, 2018) and have shown success on medium to moderately large datasets
(Wang et al.,, 2019), outperforming modeling-based approaches such as
stochastic variational GPs (SVGP) (Hensman et al.,, 2013).

Most GP training and inference approaches based on iterative methods use
classic general-purpose algorithms for matrix solves, such as conjugate
gradients (CG) (Cutajar et al.,, 2016; Gardner et al.,, 2018; Wang et al.,,
2019), MINRES (Pleiss et al.,, 2020), or (stochastic) gradient descent (Lin et
al.,, 2023). There is reason to believe that such algorithms are suboptimal
for modern hardware-accelerated Gaussian processes. For example, CG was
purpose-built for sparse linear systems that require high-precision solutions.
Neither of these properties applies to GP regression: the necessary solves
involve dense covariance matrices, and tasks such as hyperparameter
optimization can be performed with extremely coarse-grained solves (Wang et
al.,, 2019; Maddox et al.,, 2022). These characteristics of large-scale dense
operations and low precision amenability are in line with existing trends in
machine learning (Courbariaux et al.,, 2015; Micikevicius et al.,, 2018), but
ultimately place Gaussian processes at odds with much of the literature on
numerical methods.

Much in the way that deep learning has been revolutionized by purpose-built
optimizers that exploit properties of neural networks (Kingma and Ba,, 2015;
Loshchilov and Hutter,, 2019), this paper aims to accelerate GPs with a
purpose-built method leveraging (coarse-grained) covariance matrix solves on
modern hardware. We introduce an iterative method to compute gradients of the
marginal log-likelihood (MLL) and the posterior mean, that improves over CG in
the following ways: 1) It requires $\mathcal{O}(n)$ computation per iteration
(rather than CG’s $\mathcal{O}(n^{2})$); 2) It converges rapidly and
monotonically in its early stages (but does not necessarily obtain higher
precision than CG); and 3) It demonstrates improved numerical stability in
floating point arithmetic.

In summary, we make the following contributions:

  * •

We propose an iterative method to train Gaussian processes, which computes the
MLL derivatives and posterior mean via _alternating projection_. Each update
accesses only subblocks of the kernel matrix, has linear complexity, and
decreases the residual near-monotonically.

  * •

We prove that our algorithm converges linearly at a rate no slower than
gradient descent, despite never operating on the full kernel matrix.
Empirically, our method achieves a 2-27$\times$ speed-up over CG on a wide
range of datasets.

  * •

As a demonstration of its scalability and robustness to ill-conditioning, we
are able to train a GP on $4$ million data points, the largest dataset
reported in the literature to-date without using inducing points or similar
modeling approximations—to the best of our knowledge. We find that our method
outperforms SVGP by a significant margin at this scale.

##  2 SETUP AND BACKGROUND

Notation. Let $(\mathbf{X},\mathbf{y})$ be a training set of $n$ training
inputs
$\mathbf{X}=(\begin{matrix}\mathbf{x}_{1}&\cdots&\mathbf{x}_{n}\end{matrix})^{\top}\in\mathcal{X}\subseteq\mathbb{R}^{n\times
d}$ and labels
$\mathbf{y}=(\begin{matrix}y_{1}&\cdots&y_{n}\end{matrix})^{\top}\in\mathbb{R}^{n}$.
Let the set $\\{1,2,\dots,n\\}$ be denoted by $[n]$. Given a matrix
$\mathbf{A}\in\mathbb{R}^{n\times n}$ and an index set $I\subseteq[n]$,
$\mathbf{A}_{I}=\mathbf{A}_{I,:}$ is the $|I|\times n$ row-indexed submatrix,
$\mathbf{A}_{:,I}$ the $n\times|I|$ column-indexed submatrix, and
$\mathbf{A}_{I,I}$ is the $|I|\times|I|$ principal submatrix. We use similar
indexing notations for vectors.

Now, let $f:\mathcal{X}\to\mathbb{R}$ be a latent function, and let
$k_{\bm{\theta}}:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a (known)
positive definite kernel function with hyperparameters $\bm{\theta}$. We write
$\mathbf{f}=f(\mathbf{X})=(\begin{matrix}f(\mathbf{x}_{1})&\cdots&f(\mathbf{x}_{n})\end{matrix})^{\top}\in\mathbb{R}^{n}$.
Similarly, $k_{\bm{\theta}}(\mathbf{X},\cdot):\mathcal{X}\to\mathbb{R}^{n}$
denotes the vector-valued function given by
$(\begin{matrix}k(\mathbf{x}_{1},\cdot)&\cdots&k(\mathbf{x}_{n},\cdot)\end{matrix})^{\top}\in\mathbb{R}^{n}$,
and $\mathbf{K}_{\bm{\theta}}\in\mathbb{R}^{n\times n}$ is the Gram matrix
with
$[\mathbf{K}_{\bm{\theta}}]_{ij}=k_{\bm{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})$.
We omit the subscript $\bm{\theta}$ unless the context needs it.

Gaussian Process Regression. In supervised GP regression, we assume a
response-generating function $f$ that is Gaussian process distributed a
priori— _i.e_. $f\sim\mathcal{GP}\big{(}\mu,k_{\bm{\theta}})$. For simplicity
of presentation, we assume without loss of generality an exact observation
model— _i.e_. $\mathbf{y}=f(\mathbf{X})$.111Note that we can easily recover an
observational noise model by setting
$k_{\bm{\theta}}(\mathbf{x},\mathbf{x}^{\prime})=k_{\mathrm{base}}(\mathbf{x},\mathbf{x}^{\prime})+\sigma^{2}\mathbbm{1}[\mathbf{x}=\mathbf{x}^{\prime},\mathbf{x}\in\mathbf{X}]$
for some $k_{\mathrm{base}}$ and $\sigma>0$. Given a finite test dataset
$\mathbf{x}_{1}^{*},\ldots,\mathbf{x}_{M}^{*}$, we can obtain a posterior
distribution over $f(\mathbf{x}^{*}_{1}),\ldots,f(\mathbf{x}^{*}_{M})$ using
standard Gaussian conditioning rules with the posterior mean and covariance:

$\displaystyle\mathbb{E}[\mathbf{f}^{*}\mid\mathbf{f}]$
$\displaystyle=\bm{\mu}+\mathbf{K}_{*\mathbf{f}}\mathbf{K}^{-1}(\mathbf{y}-\bm{\mu}),$
$\displaystyle\mathbb{C}[\mathbf{f}^{*}\mid\mathbf{f}]$
$\displaystyle=\mathbf{K}_{**}-\mathbf{K}_{*\mathbf{f}}\mathbf{K}^{-1}\mathbf{K}_{\mathbf{f}*}.$

We refer the reader to Rasmussen and Williams, (2006, Ch. 2) for more details.

Hyperparameter Training. The hyperparameters $\bm{\theta}$ of the GP are
learned by minimizing the negative marginal log likelihood (MLL)
$\ell(\bm{\theta}):=-\log p(\mathbf{y};\bm{\theta})$. With a Gaussian process
prior on $f$, we have
$p(\mathbf{y};\bm{\theta})=\mathcal{N}(\mathbf{y};\bm{\mu},\mathbf{K}_{\bm{\theta}})$,
yielding the following minimization:

$\displaystyle\operatorname*{\mathrm{minimize}}_{\bm{\theta}}\ell(\bm{\theta})\overset{c}{=}{\textstyle\frac{1}{2}}\left(\mathbf{y}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{y}+\log\det(\mathbf{K}_{\bm{\theta}})\right)$
(1)

Equation 1 is commonly optimized with first-order methods, which require an
(unbiased) estimate of $\frac{\partial\ell(\bm{\theta})}{\partial\theta}$.
Unfortunately, as (1) cannot be written in the usual
$\sum_{i=1}^{n}\ell(\mathbf{x}_{i},y_{i})$ form common to many machine
learning algorithms, standard minibatching strategies are not readily
applicable. Following prior work (_e.g_. Cutajar et al.,, 2016; Gardner et
al.,, 2018; Wenger et al., 2022a, ), we use the following unbiased estimate:

$\displaystyle-\tfrac{1}{2}\mathbf{y}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\tfrac{\partial\mathbf{K}_{\bm{\theta}}}{\partial\theta}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{y}+{\textstyle\frac{1}{2l}\sum_{i=1}^{l}}\Big{(}\mathbf{z}_{i}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\Big{)}\tfrac{\partial\mathbf{K}_{\bm{\theta}}}{\partial\theta}\mathbf{z}_{i},$
(2)

where $\mathbf{z}_{i}$ are i.i.d. random vectors with
$\mathbb{E}\left[\mathbf{z}_{i}\right]=\bm{0}$ and
$\mathbb{E}\left[\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\right]=\mathbf{I}$. Note
that the first term is an unbiased approximation of
$\mathbf{tr}\left(\mathbf{K}_{\bm{\theta}}^{-1}\frac{\partial\mathbf{K}_{\bm{\theta}}}{\partial\theta}\right)$.
Crucially, computing (2) primarily involves computing solves with
$\mathbf{K}_{\bm{\theta}}$.

Linear Solves via Iterative Methods. When $\mathbf{K}$ is large, direct
methods for solving $\mathbf{K}\mathbf{w}=\mathbf{b}$ are prohibitively slow.
Iterative methods, such as conjugate gradients (CG), offer reduced asymptotic
complexity (Cutajar et al.,, 2016), significant GPU acceleration (Gardner et
al.,, 2018), and memory savings if $\mathbf{K}$ is accessed in a map-reduce
fashion (Wang et al.,, 2019; Charlier et al.,, 2021).

CG minimizes the quadratic objective
$\tfrac{1}{2}\mathbf{w}^{\top}\mathbf{K}\mathbf{w}-\mathbf{b}^{\top}\mathbf{w}$
by iteratively searching along conjugated directions. Each iteration requires
a $\mathcal{O}(n^{2})$ matrix-vector multiplication with $\mathbf{K}$. In
exact arithmetic, CG returns an exact solution after $n$ iterations. In
practice for ill-conditioned problems, CG is terminated once the residual
$\mathbf{r}=\mathbf{b}-\mathbf{K}\mathbf{w}$ is small enough, _e.g_.,
$\|\mathbf{r}\|\\!\leq\delta\|\mathbf{b}\|$ for some predefined tolerance
parameter $\delta$.

For GP hyperparameter learning often large values of the tolerance $\delta$
are used despite the potential for overfitting (Potapczynski et al.,, 2021),
for example $\delta=1$ is used in practice (Wang et al.,, 2019; Maddox et
al.,, 2022) and has been the default setting of CG during training in popular
GP software packages (_e.g_., GPyTorch222GPyTorch setting https://rb.gy/qi8er
and GPflow333GPflow setting https://rb.gy/mozif).

For hyperparameter training, each MLL derivative evaluation requires a batched
linear solve $\mathbf{K}\mathbf{W}=\mathbf{B}$, where
$\mathbf{B}=(\begin{matrix}\mathbf{y}&\mathbf{z}_{1}&\dots&\mathbf{z}_{l}\end{matrix})$
with $\mathbf{z}_{i}$ are random samples for stochastic MLL derivative
estimation in (2).

RKHS. Every kernel $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ induces a
space of functions
$\mathcal{H}:=\overline{\mathrm{span}}\\{k(\mathbf{x},\cdot):\mathbf{x}\in\mathcal{X}\\}\subset\mathbb{R}^{\mathcal{X}}$,
known as a reproducing kernel Hilbert space (RKHS) where the inner product
$\langle\cdot,\cdot\rangle_{\mathcal{H}}$ is defined as $\langle
k(\mathbf{x},\cdot),k(\mathbf{x}^{\prime},\cdot)\rangle_{\mathcal{H}}=k(\mathbf{x},\mathbf{x}^{\prime})$
for all $\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$.

RKHS Projection. Define the following finite dimensional linear subspaces of
$\mathcal{H}$ for indices $I\subseteq[n]$:

$\displaystyle V_{[n]}$
$\displaystyle=\mathrm{span}\\{k(\mathbf{x}_{i},\cdot):i=1,2,\cdots,n\\}\subset\mathcal{H},$
(3) $\displaystyle V_{I}$
$\displaystyle=\mathrm{span}\\{k(\mathbf{x}_{i},\cdot):i\in I\\}\subseteq
V_{[n]},$

By definition these subspaces contain functions of the form
$f(\cdot)=\sum_{i=1}^{n}\alpha_{i}k(\mathbf{x}_{i},\cdot)$ and
$f(\cdot)=\sum_{i\in I}\alpha_{i}k(\mathbf{x}_{i},\cdot)$ respectively. We can
map any $f\in\mathcal{H}$ onto these subspaces using the projection operator.

######  Definition 1 (Projection Operator).

Let $V\subseteq\mathcal{H}$ be a closed linear subspace. The projection of any
$f\in\mathcal{H}$ onto $V$ is given by the _projection operator_

$\displaystyle\operatorname{proj}_{V}(f)=\operatorname*{argmin}_{g\in
V}~{}~{}\tfrac{1}{2}\|f-g\|_{\mathcal{H}}^{2},$

which is well-defined, i.e. the unique minimizer exists.

Intuitively, the projection operator finds the best approximation of $f$ in
$V$, where approximation error is measured by the norm
$\lVert\cdot\rVert_{\mathcal{H}}$. For $V=V_{[n]}$ and $V=V_{I}$, the
projection operator has a simple form:

$\displaystyle\operatorname{proj}_{V_{[n]}}(f)$
$\displaystyle=f(\mathbf{X})^{\top}\mathbf{K}^{-1}k(\mathbf{X},\cdot),$ (4)
$\displaystyle\operatorname{proj}_{V_{I}}(f)$
$\displaystyle=f(\mathbf{X})^{\top}\mathbf{E}_{I}^{\top}\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}k(\mathbf{X},\cdot).$

Importantly, these projections only evaluate $f$ and the kernel $k$ on the
data $\mathbf{X}$ (or subset $\mathbf{X}_{I}$). In other words, it is
unnecessary to evaluate $f$ or $k$ outside of $\mathbf{X}$ (or
$\mathbf{X}_{I}$). The complexity of computing the projection
$\operatorname{proj}_{V}(f)$ depends on the dimension of $V$:
$\operatorname{proj}_{V_{[n]}}(f)$ takes $\mathcal{O}(n^{3})$ time and
$\operatorname{proj}_{V_{I}}(f)$ takes $\mathcal{O}(\lvert I\rvert^{3})$ time.

##  3 METHOD

In this section, we develop an iterative method for computing solves
$\mathbf{K}^{-1}\mathbf{b}$ by alternating projection. The method supports
batch linear solves with multiple right-hand sides, as required by estimating
the MLL derivative (2), and is amenable to GPU parallelism. We cast the linear
solve as a projection in the RHKS $\mathcal{H}$ and decompose the projection
into a sequence of small-scale subproblems. Each subproblem is solved in
$\mathcal{O}(n)$ time, allowing frequent updates. Alternating projection
typically makes rapid progress in the early stage and finds a medium-precision
solution quickly.

$g=r_{0}$$k(\mathbf{x}_{1},\cdot)$$k(\mathbf{x}_{2},\cdot)$$k(\mathbf{x}_{3},\cdot)$$s_{1}$$r_{1}$

Figure 2:  Left: Illustration of alternating projection. $s_{1}$ is the
projection of $g=r_{0}$ onto the subspace spanned by $k(\mathbf{x}_{1},\cdot)$
and $k(\mathbf{x}_{2},\cdot)$. The residual $r_{1}=g-s_{1}$ will be projected
to other coordinates in the next iteration. Right: Gauss-Southwell block
selection results in faster convergence than random/cyclic.

High Level Approach. Assume $k$ is strictly positive definite and there is no
duplicate data, then there exists $g\in\mathcal{H}$ interpolating
$\mathbf{b}$, i.e. $g(\mathbf{X})=\mathbf{b}$. The exact form of $g$ is not
important (or unique for that matter); rather, we are interested in its
projection onto the subspace $V_{[n]}$, which by (4) is

$\displaystyle\operatorname{proj}_{V_{[n]}}(g)=\mathbf{b}^{\top}\mathbf{K}^{-1}k(\mathbf{X},\cdot),$

Thus the linear solve can be obtained from the coefficients of the projection
$\operatorname{proj}_{V_{[n]}}(g)$.

Directly projecting $g$ onto $V_{[n]}$ is computationally infeasible, as the
time complexity is cubic in $n$. Instead, we partition $[n]$ into subsets
$\mathcal{P}=\\{I_{1},I_{2},\cdots,I_{m}\\}$. For each subset
$I\in\mathcal{P}$, the projection to the linear subspace $V_{I}\subseteq
V_{[n]}$ is cheap, provided that $\lvert I\rvert$ is small. Thus, we construct
the (full) projection $\operatorname{proj}_{V_{[n]}}(g)$ by iteratively
computing the projection onto the linear subspaces $V_{I}$ where
$I\in\mathcal{P}$.

Starting from $r_{0}=g$ and $s_{0}=0$, the $j$-th iteration selects an index
set $I\subseteq[n]$ and updates as follows

$\displaystyle s_{j+1}=s_{j}+\operatorname{proj}_{V_{I}}(r_{j})$ (5)
$\displaystyle r_{j+1}=r_{j}-\operatorname{proj}_{V_{I}}(r_{j})$ (6)

Intuitively, $s_{j}$ progressively approximates the true projection
$\operatorname{proj}_{V_{[n]}}(g)$, since (5) iteratively adds the projection
onto subspaces $V_{I}$ to the current approximation $s_{j}$. Meanwhile, (6)
consistently updates the residual. As $j\to\infty$, $s_{j}$ converges to the
true projection $\mathbf{b}^{\top}\mathbf{K}^{-1}k(\mathbf{X},\cdot)$
(Wendland,, 2004). See Figure 2 (left panel) for an illustration of
alternating projection.

Implicit Representation of $r_{j}(\cdot)$ Crucially, in the updates (5) and
(6), the function $r_{i}$ is only ever accessed through its evaluation on
$\mathbf{X}$ (recall the projection formula (4)). Therefore, we only need to
maintain the vector $\mathbf{r}_{i}=r_{i}(\mathbf{X})\in\mathbb{R}^{n}$
instead of the entire function. The update (6) thus reduces to

$\displaystyle\mathbf{r}_{j+1}$
$\displaystyle:=\mathbf{r}_{j}-\operatorname{proj}_{V_{I}}(r_{j})(\mathbf{X})$
$\displaystyle=\mathbf{r}_{j}-\mathbf{K}\mathbf{E}_{I}^{\top}\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}\mathbf{r}_{j}$
(7)
$\displaystyle=\mathbf{r}_{j}-\mathbf{K}_{:,I}\mathbf{K}_{I,I}^{-1}[\mathbf{r}_{j}]_{I},$
(8)

where $\mathbf{E}_{I}$ denotes the rows of the identity matrix corresponding
to $I$. The final line comes from the right multiplication
$\mathbf{K}\mathbf{E}_{I}^{\top}$ and left multiplication
$\mathbf{E}_{I}\mathbf{r}_{j}$.

Representing $s_{i}(\cdot)$ via Kernel Functions. Every $s_{i}$ is in
$V_{[n]}$ and can thus be written as a linear combination
$\mathbf{w}_{i}^{\top}k(\mathbf{X},\cdot)$ for some
$\mathbf{w}_{i}\in\mathbb{R}^{n}$, which is proved by induction. At the $0$-th
iteration, we see that $s_{0}(\cdot)$ is the zero function, which can be
written as $\bm{0}^{\top}k(\mathbf{X},\cdot)$. For the $j$-th iteration,
assuming $I\subseteq[n]$ is selected and
$s_{j}=\mathbf{w}_{j}^{\top}k(\mathbf{X},\cdot)$, then we have

$\displaystyle s_{j+1}$
$\displaystyle=s_{j}+\operatorname{proj}_{V_{I}}(r_{j})$
$\displaystyle=\underbrace{\left(\mathbf{w}_{j}^{\top}+r_{j}(\mathbf{X})^{\top}\mathbf{E}_{I}^{\top}\mathbf{K}^{-1}_{I,I}\mathbf{E}_{I}\right)}_{\mathbf{w}_{j+1}}k(\mathbf{X},\cdot),$

where the last line gives an explicit update on $\mathbf{w}_{j}$:

$\displaystyle\mathbf{w}_{j+1}=\mathbf{w}_{j}+\mathbf{E}_{I}^{\top}\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}\mathbf{r}_{j}.$
(9)

Recall that $\mathbf{E}_{I}$ simply selects rows/columns. Only entries in
$\mathbf{w}_{j}$ indexed by $I$ need to be updated, while keeping the entries
outside $I$ unchanged:

$\displaystyle\begin{split}[\mathbf{w}_{j+1}]_{I}&=[\mathbf{w}_{j}]_{I}+\mathbf{K}_{I,I}^{-1}[\mathbf{r}_{j}]_{I},\\\
[\mathbf{w}_{j+1}]_{[n]\setminus I}&=[\mathbf{w}_{j}]_{[n]\setminus
I}.\end{split}$ (10)

Summary. (8) and (10) yield an iteration on
$s_{i}(\cdot)\\!=\\!\mathbf{w}_{i}^{\top}k(\mathbf{X},\cdot)$ where the
$\mathbf{w}_{i}$ are obtained through simple matrix operations. Since the
$s_{i}$ are produced by alternating projections, we have
$s_{i}\to\operatorname{proj}_{V_{[n]}}(g)$ and thus
$\mathbf{w}_{i}\to\mathbf{K}^{-1}\mathbf{b}$. We summarize this approach in
Algorithm 1. Note that the algorithm can be adapted to perform multiple right-
hand solves in parallel by replacing
$\mathbf{w}_{i},\mathbf{r}_{i},\mathbf{b}_{i}$ vectors with matrices
$\mathbf{W},\mathbf{R},\mathbf{B}$.

Input: A batched linear system $\mathbf{K}\mathbf{W}\\!=\\!\mathbf{B}$

Output: The solution $\mathbf{W}^{*}=\mathbf{K}^{-1}\mathbf{B}$

1 Initialize $\mathbf{W}=\mathbf{O}$ and $\mathbf{R}=\mathbf{B}$

2 for  _$t=1,2,\cdots$_ do // epoch

3 for  _$j=1,2,\cdots,m$_ do // mini-batch

4 select a block $I\in\mathcal{P}$ from the partition

5 $\mathbf{W}_{I}=\mathbf{W}_{I}+\mathbf{K}_{I,I}^{-1}\mathbf{R}_{I}$

6 $\mathbf{R}=\mathbf{R}-\mathbf{K}_{:,I}\mathbf{K}_{I,I}^{-1}\mathbf{R}_{I}$

7 end for

8 if  _converged_ then return $\mathbf{W}$

9

10 end for

Algorithm 1 Alternating Projection

Block Selection. Selecting which block to update is crucial for fast
convergence. The simplest block selection rules are random selection (sample
$I$ uniformly from $\mathcal{P}$) and cyclic selection ($I=I_{j}$), which
usually converge slowly (see Figure 2). Instead, we select the block $I$ with
the largest residual norm

$\displaystyle
I=\operatorname*{argmax}_{I\in\mathcal{P}}~{}\lVert\mathbf{R}_{I,:}\rVert_{\mathrm{F}}^{2}.$
(11)

In the special case that $\mathbf{R}$ is an $n\times 1$ vector, (11) reduces
to the Gauss-Southwell (GS) rule (Nutini et al.,, 2015). (11) is a
modification adapted to our setting.

Cached Cholesky. Algorithm 2 requires solving a linear system with the
submatrix $\mathbf{K}_{I,I}$. To avoid repeatedly inverting the same matrices,
we compute and cache the Cholesky factors of all principal submatrices
$\\{\mathbf{K}_{I,I}:I\in\mathcal{P}\\}$ once whenever the GP hyperparameters
are updated (e.g., once per gradient computation). To facilitate parallelism,
we partition the blocks evenly so that every block has the same size $\lvert
I\rvert=b$ and factorize all matrices in a single batch Cholesky call, which
takes $\mathcal{O}(nb^{2})$ time and $\mathcal{O}(nb)$ memory.

Complexity. The block selection takes $\mathcal{O}(n)$ time. Updating the
weights $\mathbf{W}$ takes $\mathcal{O}(b^{2})$ time. Updating the residual
$\mathbf{R}$ takes $\mathcal{O}(nb)$ time. Each epoch runs $m=n/b$ inner loops
and thus takes $\mathcal{O}(nb+n^{2})$ time in total. Thus, the complexity of
each epoch has the same quadratic complexity as a single CG iteration. A more
fine-grained analysis in Appendix F shows that each epoch requires
$(2+\frac{3}{b})n^{2}+(2b+1)n$ FLOPs. Thus, for typical batch sizes $1\ll b\ll
n$, each epoch requires roughly $2n^{2}$ FLOPs, the same number as a single CG
iteration. We note that every update in Algorithm 1 has linear (in terms of
$n$) time and memory complexity.

Connection with Coordinate Descent. Interestingly, we can show that Algorithm
1 produces iterates equivalent to coordinate descent on the quadratic form
(see §A for details). We will exploit this connection to prove the rate of
convergence of Algorithm 1. We introduce this algorithm as alternating
projection for two reasons: (a) unlike in coordinate descent, the update rules
based on alternating projection maintain the residual $\mathbf{R}$, which
enables efficient block selection strategies like the GS rule without re-
evaluating the residual; (b) alternating projection can be easily extended to
different settings. For instance, a parallel coordinate descent algorithm was
discovered via the connection with (Dykstra’s) alternating projection (Boyle
and Dykstra,, 1986; Tibshirani,, 2017) in the setting of regularized least-
squares, which hints that Algorithm 1 may be distributed.

##  4 CONVERGENCE

Figure 3:  Convergence of alternating projection with different batch sizes
$b$ on 3droad. Left: Smaller batch sizes converge faster within the same
epochs. Right: However, smaller batch sizes result in more sequential updates
on the GPU and thus longer wall-clock time.

Let $\lambda_{\max}$ and $\lambda_{\min}$ be the largest and smallest
eigenvalues of $\mathbf{K}$, $\kappa=\lambda_{\max}/\lambda_{\min}$ its
condition number, and define
$\lambda_{\max}^{\prime}=\max_{I\in\mathcal{P}}\lambda_{\max}(\mathbf{K}_{I,I})$
as the maximum of the largest eigenvalues of the principal submatrices
$\\{K_{I,I}:I\in\mathcal{P}\\}$. By leveraging the connection with coordinate
descent (Nutini et al.,, 2022), we can prove an explicit convergence rate for
Algorithm 1 when applied to a linear system with multiple right-hand sides.

######  Theorem 1.

Let $\mathbf{W}^{*}$ be the (unique) solution of the linear system
$\mathbf{K}\mathbf{W}=\mathbf{B}$ and $\mathbf{W}^{(t)}$ its approximation
after $t$ epochs of Algorithm 1 using the modified GS rule (11). Then it holds
that

$\displaystyle\lVert\mathbf{W}^{(t)}-\mathbf{W}^{*}\rVert_{\mathbf{K}}^{2}\leq\exp\big{(}-t/\kappa^{\prime}\big{)}\lVert\mathbf{W}^{(0)}-\mathbf{W}^{*}\rVert_{\mathbf{K}}^{2},$

where $\kappa^{\prime}=\lambda_{\max}^{\prime}/\lambda_{\min}\leq\kappa$.

The rate in Theorem 1 improves over gradient descent despite only needing sub
matrices, for which the above holds with $\exp(-t/\kappa)$, since generally
$\kappa^{\prime}\leq\kappa$. For comparison, the convergence rate of (batched)
CG is $4\big{(}(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)\big{)}^{2t}\approx
4\exp\big{(}-4t/\sqrt{\kappa}\big{)}$ for a sufficiently large condition
number $\kappa\gg 1$. The convergence rate of alternating projection is
asymptotically faster than that of CG if
$\kappa^{\prime}\leq\frac{1}{4}\sqrt{\kappa}$. In general, we do not expect
this condition to hold. However, alternating projection has practical
advantages despite a slower asymptotic convergence rate. First, alternating
projection performs $m$ times more updates than CG with the same number of
FLOPs. Second, alternating projection generally decreases the residual in
every epoch, while the CG residual is not monotonic. We empirically observe
that CG often increases the residual dramatically in the early stage and it
takes time for CG to enter the “linear convergence phase”. In addition, the
dependency on $\kappa^{\prime}$ suggests that alternating projection
implicitly works on better-conditioned matrices, which may imply robustness
against ill-conditioning.

Figure 1 shows the above two points in practice. The figure is plotted on two
checkpoints at the $50$ epoch GP training on the 3droad and house electric
datasets respectively. The (batched) linear system $\mathbf{K}^{-1}\mathbf{B}$
has $16$ right-hand sides, where $\mathbf{b}_{0}=\mathbf{y}$ is the training
labels and $\\{\mathbf{b}_{i}\\}_{i=1}^{15}$ are _i.i.d_. samples from a
Gaussian. We can prove that the random selection strategy in Figure 2 (right
panel) achieves a similar rate in Theorem 1, but only in expectation. In
practice, the GS rule converges faster than random selection.

The batch size $b$ affects the rate in Theorem 1 through the condition number
$\kappa^{\prime}=\lambda_{\max}^{\prime}/\lambda_{\min}$. Note that the
largest eigenvalue of the principal submatrix is bounded by its trace
$\lambda_{\max}(\mathbf{K}_{I,I})\leq\mathbf{tr}\left(\mathbf{K}_{I,I}\right)$,
where the trace grows linearly in $\lvert I\rvert$. A small batch size
$b=\lvert I\rvert$ is likely to have a small $\lambda_{\max}^{\prime}$ and a
faster convergence rate. We compare the convergence of different batch sizes
in Figure 3. Although small batch sizes lead to faster convergence, they
generally have a longer running time due to more sequential updates.
Therefore, in practice, we recommend using the largest batch size possible
subject to memory constraints. In addition, we note that the convergence rate
in Theorem 1 is loose for large batch sizes $b$. In the extreme case where
$b=n$, Algorithm 1 is equivalent to the Cholesky decomposition on the entire
matrix $\mathbf{K}$ and thus converges to the exact solution in one update.
However, Theorem 1 does not reflect that. The convergence rate in practice may
be much faster than the theory predicts.

##  5 EXPERIMENTS

We evaluate the efficacy of our alternating projections solver in a GP
regression task. Our evaluation includes a training dataset of $n=4M$, which,
to the best of our knowledge, is considerably larger than any other dataset
where a GP has been applied without inducing points or employing modeling
approximations.

All experiments are performed on a single $24$ GB NVIDIA RTX A5000 GPUs with
single precision floating point, and all numerical algorithms/GP models are
implemented in PyTorch/GPyTorch (Gardner et al.,, 2018). We use the KeOps
library (Charlier et al.,, 2021) to implement all matrix-free numerical
algorithms in a map-reduce fashion, thus eliminating the need to store large
$n\times n$ kernel matrices in memory.

###  5.1 Main Result: GP Regression

We first evaluate our method on large-scale GP training tasks. We compare
against GPs trained with CG, which is the predominant matrix-free GP training
approach (Gardner et al.,, 2018; Wang et al.,, 2019; Maddox et al.,, 2022).

Metrics. Our primary desiderata for GPs are 1) low computational costs for
training and 2) generalization. Therefore, we compare the different training
methods using the following metrics: 1) the total number of floating point
operations (FLOPs) normalized by $2n^{2}$ (the FLOPs of a single matmul), 2)
the wall clock training time, and 3/4) the trained model’s RMSE and NLL
measured on the test set.

Datasets and Models. We conduct experiments on UCI regression datasets, whose
statistics are shown in Section E.6. Each dataset is split into $80\%$
training and $20\%$ test. The labels are normalized so that they have zero
mean and unit variance. Almost all experiments are averaged over $5$ runs.
Because of resource constraints, we limit the two largest datasets—House
Electic and Gas Sensors—to 3 and 1 run respectively.

We train GP regression models with $\nu=2.5$ Matérn kernels and a constant
prior mean. We optimize the following hyperparameters: a scalar constant for
the prior mean, a $d$-dimensional kernel lengthscale, a scalar outputscale,
and a scalar observational noise parameter $\sigma^{2}$. We include
experiments with $\nu=1.5$ Matérn kernels in Appendix E.

MLL Optimization. To compute the stochastic MLL gradient (2), we use
$l\\!=\\!15$ random samples $\mathbf{z}_{i}$. Thus, all matrix-free methods
solve a batched linear system with $16$ right-hand sides $\mathbf{y}$ and
$\\{\mathbf{z}_{i}\\}_{i=1}^{15}$ in each training iteration. On the first
five datasets, the GPs are trained by $50$ iterations of Adam with a step size
$0.1$. On house electric and gas sensors, the GPs are trained by $100$
iterations of Adam with a step size $0.1$.

Alternating Projection Details. As discussed in §4, a large batch size is
preferred empirically. We use the largest batch size that we can fit on a $24$
GB GPU. The batch sizes $b$ are set as: 6000 on SGEMM, air quality and 3droad;
4000 on song and buzz; 1000 on house electric; 500 on gas sensors. We use the
sequential partition $\mathcal{P}$: the data points from $(j-1)b+1$ to $jb$
belong to the $j$-th block $I_{j}$ for $j=1,2,\cdots n/b$.

The maximum CG iterations and the maximum number of alternating projection
epoch is set to $1000$. Following GPyTorch’s CG stopping criteria, we
terminate the alternating projection solves after (a) the average relative
residual norm is strictly smaller than the tolerance $\delta=1$ or (b) $1000$
total epochs, whichever comes first. However, we ensure that at least $11$
epochs of alternating projections have been run before termination (again
following GPyTorch). We define the average relative residual norm as
$\frac{1}{l+1}\sum_{i=0}^{l}\lVert\mathbf{r}_{i}\rVert/{\lVert\mathbf{b}_{i}\rVert}$
when there are $l+1$ right hand sides
$(\begin{matrix}\mathbf{b}_{0}&\mathbf{b}_{1}&\cdots&\mathbf{b}_{l}\end{matrix})$.

CG Details. We use GPyTorch’s implementation of CG, which uses the same
stopping criteria as our alternating projections implementation. Following
Wang et al., (2019); Wenger et al., 2022a , we use a pivoted Cholesky
preconditioner of size $500$ on all datasets except: house electric uses a
size $300$ and gas sensors uses a size $150$ due to GPU memory overflow.

Prediction. At test time, the predictive mean is computed by the same
iterative method used for training (e.g., CG for the CG trained GP,
alternating projection for the AP trained GP). A limitation of our method is
that it does not easily result in a cache for variances. Therefore, we use
1000 Lanczos iterations as in Pleiss et al., (2018); Wang et al., (2019).

Figure 4:  GP training on air quality dataset. Left: Because the likelihood
noise $\sigma^{2}$ decreases during training, the matrix $\mathbf{K}$ gets
more ill-conditioned. Right: CG is sensitive to this increased ill-
conditioning, while alternating projections is robust.  Table 1:  _Gaussian
process training on UCI benchmark datasets._ Metrics are computed across
multiple runs and reported with $\pm$ one standard deviation.

Dataset | Method | RMSE | NLL | FLOPs$/2n^{2}$ | Training time | Speed up
---|---|---|---|---|---|---
SGEMM $n=241,600$ $d=14$ | CG | $0.048\pm 0.000$ | $\bm{-1.037\pm 0.001}$ | $551\pm 1$ |  $9.1$m $\pm 0.0$ |
Alt. Proj. | $\bm{0.046\pm 0.000}$ | $-0.999\pm 0.001$ | $550\pm 0$ |  $12.2$m $\pm 0.2$ | $0.7\times$
2-7 | SVGP | $0.086\pm 0.000$ | $-0.934\pm 0.003$ | NA |  $14.8$m $\pm 0.1$ |
air quality $n=382,168$ $d=13$ | CG | $\bm{0.261\pm 0.001}$ | $0.143\pm 0.004$ | $2965\pm 19$ |  $33.5$m $\pm 1.5$ |
Alt. Proj. | $\bm{0.262\pm 0.001}$ | $\bm{0.137\pm 0.003}$ | $550\pm 0$ |  $16.9$m $\pm 0.5$ | $2.0\times$
2-7 | SVGP | $0.363\pm 0.003$ | $0.399\pm 0.006$ | NA |  $23.4$m $\pm 0.1$ |
3droad $n=434,874$ $d=3$ | CG | $\bm{0.069\pm 0.000}$ | $1.324\pm 0.002$ | $5128\pm 114$ |  $53.2$m $\pm 2.8$ |
Alt. Proj. | $0.076\pm 0.000$ | $1.203\pm 0.001$ | $676\pm 1$ |  $21.1$m $\pm 0.5$ | $2.5\times$
2-7 | SVGP | $0.327\pm 0.002$ | $\bm{0.320\pm 0.005}$ | NA |  $26.1$m $\pm 0.1$ |
song $n=515,345$ $d=90$ | CG | $\bm{0.747\pm 0.002}$ | $1.140\pm 0.003$ | $4431\pm 110$ |  $13.8$h $\pm 0.8$ |
Alt. Proj. | $\bm{0.749\pm 0.002}$ | $\bm{1.132\pm 0.002}$ | $550\pm 0$ |  $2.7$h $\pm 0.1$ | $5.1\times$
2-7 | SVGP | $0.790\pm 0.002$ | $1.184\pm 0.002$ | NA |  $0.5$h $\pm 0.0$ |
buzz $n=583,250$ $d=77$ | CG | $0.321^{*}\pm 0.144$ | $0.669^{*}\pm 1.152$ | $16726\pm 2724$ |  $31.1$h $\pm 5.4$ |
Alt. Proj. | $\bm{0.239\pm 0.001}$ | $\bm{0.018\pm 0.003}$ | $550\pm 0$ |  $2.0$h $\pm 0.1$ | $15.6\times$
2-7 | SVGP | $0.259\pm 0.002$ | $0.066\pm 0.006$ | NA |  $0.6$h $\pm 0.0$ |
house electric $n=2,049,280$ $d=11$ | CG | - | - | $\geqslant 50441$ |  $\geqslant 11$d |
Alt. Proj. | $\bm{0.030\pm 0.000}$ | $-1.148\pm 0.001$ | $1100\pm 0$ |  $9.8$h $\pm 0.4$ | $\geqslant 26.9\times$
2-7 | SVGP | $0.050\pm 0.000$ | $\bm{-1.549\pm 0.001}$ | NA |  $2.1$h $\pm 0.0$ |
gas sensors $n=4,178,504$ $d=17$ | CG | - | - | - | - |
Alt. Proj. | $\bm{0.203}$ | $\bm{0.070}^{\dagger}$ | 1100 |  $84.5$h |
2-7 | SVGP | $0.330\pm 0.001$ | $0.339\pm 0.003$ | NA |  $8.7$h $\pm 0.03$ |
* : At test time, CG does not reach the tolerance $\delta=0.01$ after $4000$ iterations on some checkpoints.
\- : CG does not finish GP training.
$\dagger$ : This predictive variance is calculated using only $500$ Lanczos
iterations to save time and avoid numerical instability.

Results on $10^{5}<n<10^{6}$ datasets. Section E.6 compares the predictive
performance and the training speed of CG-based versus alternating projection-
based GPs. Both training procedures produce GPs with similar RMSE and NLL. We
conjecture that this similarity occurs because both approaches solve linear
systems up to the same tolerance, and thus find similar hyperparameters. One
exception is the buzz dataset: CG struggles to converge while training on this
dataset, resulting in considerably worse RMSE and NLL.

The primary difference between the two methods is training time. Alternating
projection-based training is up to $27\times$ faster than CG. The only
exception is SGEMM GPU, which seems to be a well-conditioned dataset since CG
converges quickly.

For reference, we also report the training/test performance of stochastic
variational Gaussian processes (SVGP) (see Appendix E for experimental design
details). GPs trained by alternating projection achieve substantially lower
RMSE and comparable NLL compared with SVGP. We do note that SVGPs have lower
NLL on 3droad and house electric, which we suspect is a limitation of the
Lanczos predictive variance estimates used on the alternating projections
models. (Note that SVGP’s predictive variances can be computed exactly and do
not make use of the Lanczos estimator.) Indeed, in Appendix E we find that the
NLL gap shrinks as we increase the rank of the Lanczos variance estimator,
suggesting that this gap is not a fundamental limitation of the alternating
projections training methodology.

Results on $n\geq 10^{6}$ datasets. Previous attempts to train GPs using
iterative methods on datasets with $n\geq 10^{6}$ examples have used a large
noise constraint $\sigma^{2}\geq 0.1$ to improve the conditioning of the
kernel matrix (_e.g_., Wang et al.,, 2019; Maddox et al.,, 2022). Since
alternating projection is much less conditioning-sensitive than CG (see §5.2),
for the first time, we are able to train the model with a much smaller noise
constraint $\sigma^{2}\geq 10^{-4}$, _i.e_. the default in GPyTorch for the
Gaussian likelihood.444GPyTorch likelihood setting https://rb.gy/fv41w
Removing the noise constraint yields much better predictive performance: the
RMSE $0.030$ is significantly lower than what can be achieved with high-noise
constraint models (see Appendix E).

We additionally train a GP on the gas sensors dataset with $4$ million data
points. To the best of our knowledge, this is the largest dataset trained on
using GPs without the use of inducing point or other modeling approximations.
CG-based training appears to be intractable on such a large dataset, requiring
over a week to train. In contrast, the alternating projections method required
84.5 hours.

###  5.2 Effect of Kernel Matrix Conditioning

As implied by our theoretical dependence on $\lambda_{\text{max}}^{\prime}$
rather than $\lambda_{\text{max}}$, we observe that our alternating
projections method is less sensitive to ill-conditioning than CG. We
demonstrate this phenomenon in Figure 4, which depicts training on the
$n\approx 400K$ air quality dataset. Over the course of training, the noise
parameter $\sigma^{2}$ decreases for both methods, resulting in an
increasingly ill-conditioned kernel matrix (as
$\lambda_{\min}\approx\sigma^{2})$. At the end of training, when
$\sigma^{2}\approx 0.01$, CG requires over $120$ iterations to
converge—$10\times$ as many iterations as the beginning of training. In
contrast, alternating projection consistently converges in $11$ iterations
despite the decreasing noise and increasing condition number. See more
datasets in Appendix E.

###  5.3 Alternating Projection at Test Time

(a) air quality

(b) buzz

Figure 5:  Running CG and alternating projection on test-time solves
$\mathbf{K}^{-1}(\mathbf{y}-\bm{\mu})$. For alternating projection, the x-axis
is the number of epochs. Left: CG has faster convergence rate, but CG does not
reach the test-time tolerance $\delta=0.01$ much faster. Right: Alternating
projection reaches the tolerance $\delta=0.01$ faster despite its slower
asymptotic rate.

Any linear solver $\mathbf{K}^{-1}\mathbf{b}$ can be used to compute the
posterior mean on the test data. We explore alternating projection at test
time, as shown in Figure 5 and Table 3 in Appendix E. With a test-time
tolerance $\delta=0.01$, the posterior mean computed by alternating projection
is practically the same as CG: the RMSE of both methods are the same up to the
$3$rd digit after the decimal point. While alternating projection is slightly
slower on medium-size datasets such air quality and 3droad, we observe strong
speed up on large datasets such as buzz and house electric. Our method
computes the posterior mean $17.2\times$ faster in wall-clock time than CG on
buzz, and requires only $5$ min to compute the posterior mean on house
electric.

##  6 RELATED WORK

The early usage of conjugate gradients in GPs dates back at least to Yang et
al., (2004); Shen et al., (2005). They proposed methods speeding up CG by
approximate matrix-vector multiplications. More recently, CG has been
revisited by Davies, (2015); Cutajar et al., (2016). Then, a series of work
(Gardner et al.,, 2018; Wang et al.,, 2019; Artemev et al.,, 2021) and
software such as GPyTorch (Gardner et al.,, 2018) and GPflow (Matthews et
al.,, 2017) have popularized CG for GPs.

Alternating projection (Von Neumann,, 1949) is a general algorithm finding a
point in the intersection of convex sets, enjoying applications in convex
optimization (Agmon,, 1954) and scattered data approximation (Wendland,,
2004). An early work applying coordinate descent with greedy block selection
for GP inference is done by Bo and Sminchisescu, (2008). However, the
algorithm is not parallelizable on modern hardware like GPUs due to the
inherent sequential nature of the greedy selection, and lacks an explicit
convergence rate with explicit constants. Lin et al., (2023) recently have
applied stochastic gradient descent for approximate GP posterior sampling.
They also observe CG struggles with convergence in ill-conditioned settings.

##  7 CONCLUSION

In this work we proposed an alternating projection method with provable linear
convergence for solving dense kernel linear systems and applied it to GP
training and inference. Our method quickly reaches commonly used tolerances
faster than CG, requires only linear time per iteration, and is highly robust
to ill-conditioning. Experiments on several large-scale benchmark datasets
show that we achieve a 2-27$\times$ speed-up over CG-based training and a
2-17$\times$ speed up over CG-based inference with an _increase in predictive
performance_. This includes results on datasets as large as $4$ million data
points which is state-of-the-art for GPs trained with iterative methods
without artificially inflating observation noise for stability.

### Acknowledgements

JW was supported by the Gatsby Charitable Foundation (GAT3708), the Simons
Foundation (542963) and the Kavli Foundation.

## References

  * Agmon, (1954) Agmon, S. (1954).  The relaxation method for linear inequalities.  Canadian Journal of Mathematics, 6:382–392.
  * Artemev et al., (2021) Artemev, A., Burt, D. R., and van der Wilk, M. (2021).  Tighter bounds on the log marginal likelihood of Gaussian process regression using conjugate gradients.  In International Conference on Machine Learning (ICML), volume 139, pages 362–372.
  * Bertin-Mahieux, (2011) Bertin-Mahieux, T. (2011).  YearPredictionMSD.  UCI Machine Learning Repository.
  * Bo and Sminchisescu, (2008) Bo, L. and Sminchisescu, C. (2008).  Greedy block coordinate descent for large scale Gaussian process regression.  In Conference on Uncertainty in Artificial Intelligence (UAI).
  * Boyle and Dykstra, (1986) Boyle, J. P. and Dykstra, R. L. (1986).  A method for finding projections onto the intersection of convex sets in Hilbert spaces.  In Advances in Order Restricted Statistical Inference: Proceedings of the Symposium on Order Restricted Statistical Inference, pages 28–47.
  * Charlier et al., (2021) Charlier, B., Feydy, J., Glaunès, J. A., Collin, F.-D., and Durif, G. (2021).  Kernel operations on the GPU, with autodiff, without memory overflows.  Journal of Machine Learning Research, 22(74):1–6.
  * Chen, (2019) Chen, S. (2019).  Beijing Multi-Site Air-Quality Data.  UCI Machine Learning Repository.
  * Courbariaux et al., (2015) Courbariaux, M., Bengio, Y., and David, J.-P. (2015).  BinaryConnect: Training deep neural networks with binary weights during propagations.  In Advances in Neural Information Processing Systems (NeurIPS), volume 28.
  * Cutajar et al., (2016) Cutajar, K., Osborne, M., Cunningham, J., and Filippone, M. (2016).  Preconditioning kernel matrices.  In International Conference on Machine Learning (ICML), volume 48 of Proceedings of Machine Learning Research, pages 2529–2538.
  * Davies, (2015) Davies, A. J. (2015).  Effective implementation of Gaussian process regression for machine learning.  PhD thesis, University of Cambridge.
  * Fonollosa, (2015) Fonollosa, J. (2015).  Gas sensor array under dynamic gas mixtures.  UCI Machine Learning Repository.
  * Gardner et al., (2018) Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G. (2018).  GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration.  In Advances in Neural Information Processing Systems (NeurIPS), volume 31.
  * Hebrail and Berard, (2012) Hebrail, G. and Berard, A. (2012).  Individual household electric power consumption.  UCI Machine Learning Repository.
  * Hensman et al., (2013) Hensman, J., Fusi, N., and Lawrence, N. D. (2013).  Gaussian processes for big data.  In Conference on Uncertainty in Artificial Intelligence (UAI).
  * Hensman et al., (2015) Hensman, J., Matthews, A., and Ghahramani, Z. (2015).  Scalable variational Gaussian process classification.  In International Conference on Artificial Intelligence and Statistics (AISTATS), volume 38, pages 351–360.
  * Jankowiak et al., (2020) Jankowiak, M., Pleiss, G., and Gardner, J. (2020).  Parametric Gaussian process regressors.  In International Conference on Machine Learning (ICML), pages 4702–4712.
  * Katzfuss and Guinness, (2021) Katzfuss, M. and Guinness, J. (2021).  A general framework for Vecchia approximations of Gaussian processes.  Statistical Science, 36(1):124–141.
  * Kaul, (2013) Kaul, M. (2013).  3D Road Network (North Jutland, Denmark).  UCI Machine Learning Repository.
  * Kelly et al., (2023) Kelly, M., Longjohn, R., and Nottingham, K. (2023).  The UCI machine learning repository.
  * Kingma and Ba, (2015) Kingma, D. and Ba, J. (2015).  Adam: A method for stochastic optimization.  In International Conference on Learning Representations (ICLR).
  * Lin et al., (2023) Lin, J. A., Antorán, J., Padhy, S., Janz, D., Hernández-Lobato, J. M., and Terenin, A. (2023).  Sampling from Gaussian process posteriors using stochastic gradient descent.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Loshchilov and Hutter, (2019) Loshchilov, I. and Hutter, F. (2019).  Decoupled weight decay regularization.  In International Conference on Learning Representations (ICLR).
  * Maddox et al., (2022) Maddox, W. J., Potapcynski, A., and Wilson, A. G. (2022).  Low-precision arithmetic for fast Gaussian processes.  In Conference on Uncertainty in Artificial Intelligence (UAI), volume 180, pages 1306–1316.
  * Matthews et al., (2017) Matthews, A. G. d. G., van der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León-Villagrá, P., Ghahramani, Z., and Hensman, J. (2017).  GPflow: A Gaussian process library using TensorFlow.  Journal of Machine Learning Research, 18(40):1–6.
  * Micikevicius et al., (2018) Micikevicius, P., Narang, S., Alben, J., Diamos, G., Elsen, E., Garcia, D., Ginsburg, B., Houston, M., Kuchaiev, O., Venkatesh, G., and Wu, H. (2018).  Mixed precision training.  In International Conference on Learning Representations (ICLR).
  * Nutini et al., (2022) Nutini, J., Laradji, I., and Schmidt, M. (2022).  Let’s make block coordinate descent converge faster: Faster greedy rules, message-passing, active-set complexity, and superlinear convergence.  Journal of Machine Learning Research, 23(131):1–74.
  * Nutini et al., (2015) Nutini, J., Schmidt, M., Laradji, I., Friedlander, M., and Koepke, H. (2015).  Coordinate descent converges faster with the Gauss-Southwell rule than random selection.  In International Conference on Machine Learning (ICML), volume 37, pages 1632–1641.
  * Paredes and Ballester-Ripoll, (2018) Paredes, E. and Ballester-Ripoll, R. (2018).  SGEMM GPU kernel performance.  UCI Machine Learning Repository.
  * Pleiss et al., (2018) Pleiss, G., Gardner, J., Weinberger, K., and Wilson, A. G. (2018).  Constant-time predictive distributions for Gaussian processes.  In International Conference on Machine Learning (ICML), volume 80, pages 4114–4123.
  * Pleiss et al., (2020) Pleiss, G., Jankowiak, M., Eriksson, D., Damle, A., and Gardner, J. (2020).  Fast matrix square roots with applications to Gaussian processes and Bayesian optimization.  In Advances in Neural Information Processing Systems (NeurIPS), volume 33, pages 22268–22281.
  * Potapczynski et al., (2021) Potapczynski, A., Wu, L., Biderman, D., Pleiss, G., and Cunningham, J. P. (2021).  Bias-free scalable Gaussian processes via randomized truncations.  In International Conference on Machine Learning (ICML), volume 139, pages 8609–8619.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. (2006).  Gaussian processes for machine learning.  MIT Press.
  * Salimbeni et al., (2018) Salimbeni, H., Cheng, C.-A., Boots, B., and Deisenroth, M. (2018).  Orthogonally decoupled variational Gaussian processes.  In Advances in Neural Information Processing Systems (NeurIPS), volume 31.
  * Shen et al., (2005) Shen, Y., Seeger, M., and Ng, A. (2005).  Fast Gaussian process regression using kd-trees.  In Advances in Neural Information Processing Systems (NeurIPS), volume 18.
  * Snelson and Ghahramani, (2005) Snelson, E. and Ghahramani, Z. (2005).  Sparse Gaussian processes using pseudo-inputs.  In Advances in Neural Information Processing Systems (NeurIPS), volume 18.
  * Tibshirani, (2017) Tibshirani, R. J. (2017).  Dykstra's algorithm, admm, and coordinate descent: Connections, insights, and extensions.  In Advances in Neural Information Processing Systems (NeurIPS), volume 30.
  * Titsias, (2009) Titsias, M. (2009).  Variational learning of inducing variables in sparse Gaussian processes.  In International Conference on Artificial Intelligence and Statistics (AISTATS), volume 5, pages 567–574.
  * Von Neumann, (1949) Von Neumann, J. (1949).  On rings of operators. reduction theory.  Annals of Mathematics, pages 401–485.
  * Wang et al., (2019) Wang, K., Pleiss, G., Gardner, J., Tyree, S., Weinberger, K. Q., and Wilson, A. G. (2019).  Exact Gaussian processes on a million data points.  In Advances in Neural Information Processing Systems (NeurIPS), volume 32.
  * Wendland, (2004) Wendland, H. (2004).  Scattered Data Approximation, volume 17.  Cambridge University Press.
  * (41) Wenger, J., Pleiss, G., Hennig, P., Cunningham, J., and Gardner, J. (2022a).  Preconditioning for scalable Gaussian process hyperparameter optimization.  In International Conference on Machine Learning (ICML), volume 162, pages 23751–23780.
  * (42) Wenger, J., Pleiss, G., Pförtner, M., Hennig, P., and Cunningham, J. P. (2022b).  Posterior and computational uncertainty in Gaussian processes.  In Advances in Neural Information Processing Systems (NeurIPS), volume 35, pages 10876–10890.
  * Yang et al., (2004) Yang, C., Duraiswami, R., and Davis, L. S. (2004).  Efficient kernel machines using the improved fast Gauss transform.  In Advances in Neural Information Processing Systems (NeurIPS), volume 17.
  * Yang et al., (2015) Yang, Z., Wilson, A., Smola, A., and Song, L. (2015).  A la Carte – Learning Fast Kernels.  In Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, volume 38, pages 1098–1106.

Large-Scale Gaussian Processes via Alternating Projection: Supplementary
Material

[sections] [sections]l1

##  Appendix A Connection between Coordinate Descent and Alternating
Projection

This section shows the connection between Algorithm 1 and coordinate descent,
as shown in Algorithm 2.

Input: A kernel linear system $\mathbf{K}\mathbf{W}=\mathbf{B}$

Output: The solution $\mathbf{K}^{-1}\mathbf{B}$

1 Initialize $\mathbf{W}=\mathbf{O}$

2 for  _$i=1,2,\cdots$_ do // epoch

3 for  _$j=1,2,\cdots,m$_ do // mini-batch

4 select a block $I\in\\{I_{1},I_{2},\cdots,I_{m}\\}$

5
$\mathbf{W}_{I}=\mathbf{K}_{I,I}^{-1}\big{(}\mathbf{B}_{I}-\mathbf{K}_{I,\neg
I}\mathbf{W}_{\neg I}\big{)}$

6

7 end for

8 if  _converged_ then

9 return $\mathbf{W}$

10

11 end for

Algorithm 2 Block Coordinate Descent

Observe that the minimizer of the quadratic objective

$h(\mathbf{W})=\frac{1}{2}\mathbf{tr}\left(\mathbf{W}^{\top}\mathbf{K}\mathbf{W}\right)-\mathbf{tr}\left(\mathbf{B}^{\top}\mathbf{W}\right)$
(12)

is exactly the solution $\mathbf{K}^{-1}\mathbf{B}$ of the linear system
$\mathbf{K}\mathbf{W}=\mathbf{B}$.

Given a partition of indices $\\{I_{1},I_{2},\cdots,I_{m}\\}$ where $I_{i}\cap
I_{j}=\emptyset$ for all $i\neq j$ and $\cup_{i=1}^{m}I_{i}=[n]$, coordinate
descent minimizes (12) by minimizing over a subset of variables
$\mathbf{W}_{I,:}$ in each iteration. The derivative of (12) is

$\displaystyle\nabla h(\mathbf{W})=\mathbf{K}\mathbf{W}-\mathbf{B}.$

Thus, the derivative w.r.t. the subblock $\mathbf{W}_{I,:}$ is

$\displaystyle[\nabla h(\mathbf{W})]_{I}$
$\displaystyle=\mathbf{K}_{I}\mathbf{W}-\mathbf{B}_{I}$
$\displaystyle=\left(\begin{matrix}\mathbf{K}_{I,I}&\mathbf{K}_{I,\neg
I}\end{matrix}\right)\left(\begin{matrix}\mathbf{W}_{I}\\\ \mathbf{W}_{\neg
I}\end{matrix}\right)-\mathbf{B}_{I}$

where the second line splits $\mathbf{K}_{I}$ and $\mathbf{W}$ into two
blocks. The index $\neg I=[n]\setminus I$ denotes the complement of $I$.
Setting the derivative to zero gives the following update

$\displaystyle\mathbf{W}_{I}^{(j+1)}=\mathbf{K}_{I,I}^{-1}\big{(}\mathbf{B}-\mathbf{K}_{I,\neg
I}\mathbf{W}_{\neg I}^{(j)}\big{)}$

which minimizes (12) over $\mathbf{W}_{I,:}$ exactly. The full algorithm of
coordinate descent is shown in Algorithm 2.

The following lemma shows the $\mathbf{R}$ matrix in Algorithm 1 is indeed the
residual of the linear system. This lemma will be useful in proving the
equivalence between Algorithm 1 and Algorithm 2.

######  Lemma 1.

Let $\mathbf{R}^{(j)}$ and $\mathbf{W}^{(j)}$ be the residual and weight after
$j$ updates of Algorithm 1. Then we have

$\displaystyle\mathbf{R}^{(j)}=\mathbf{B}-\mathbf{K}\mathbf{W}^{(j)}.$

###### Proof.

The proof is an induction on the number of updates $j$. At the initialization
$j=0$, the equality holds trivially. Suppose after the $j$-th update we have
$\mathbf{R}^{(j)}=\mathbf{B}-\mathbf{K}\mathbf{W}^{(j)}$. All we need to do is
to verify this equality in the case of $j+1$ by direct calculation:

$\displaystyle\mathbf{B}-\mathbf{K}\mathbf{W}^{(j+1)}$
$\displaystyle=\mathbf{B}-\mathbf{K}\big{(}\mathbf{W}^{(j)}+\mathbf{E}_{I}^{\top}\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}\mathbf{R}^{(j)}\big{)}$
$\displaystyle=\mathbf{R}^{(j)}-\mathbf{K}\mathbf{E}_{I}^{\top}\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}\mathbf{R}^{(j)}$
$\displaystyle=\mathbf{R}^{(j+1)}$

where the first line uses the update rule (7) of $\mathbf{W}^{(j)}$ and the
last line uses the update rule (9) of $\mathbf{R}^{(j)}$. ∎

With Lemma 1, now we can show the equivalence between Algorithm 1 and
Algorithm 2.

######  Lemma 2.

Let $\mathbf{W}^{(j)}$ be the weight produced by Algorithm 1 after $j$
updates. Them, we have

$\displaystyle\mathbf{W}_{I}^{(j+1)}$
$\displaystyle=\mathbf{K}_{I,I}^{-1}\big{(}\mathbf{B}_{I}-\mathbf{K}_{I,\neg
I}\mathbf{W}_{\neg I}^{(j)}\big{)}$ $\displaystyle\mathbf{W}_{\neg I}^{(j+1)}$
$\displaystyle=\mathbf{W}_{\neg I}^{(j)}$

where $\neg I=[n]\setminus I$. Thus, Algorithm 1 produces the same iterates as
Algorithm 2.

###### Proof.

Recalling the update rule (9), we have

$\displaystyle\mathbf{W}^{(j+1)}=\mathbf{W}^{(j)}+\mathbf{E}_{I}^{\top}\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}\mathbf{R}^{(j)}.$

Thanks to $\mathbf{E}_{I}^{\top}$, entries outside $I$ are unchanged and thus
$\mathbf{W}_{\neg I}^{(j+1)}=\mathbf{W}_{\neg I}^{(j)}$. It remains to verify
the entries indexed by $I$. By Lemma 1, we can plug in
$\mathbf{R}^{(j)}=\mathbf{B}-\mathbf{K}\mathbf{W}^{(j)}$ and thus

$\displaystyle\mathbf{W}_{I}^{(j+1)}$
$\displaystyle=\mathbf{W}_{I}^{(j)}+\mathbf{K}_{I,I}^{-1}\mathbf{E}_{I}\big{(}\mathbf{B}-\mathbf{K}\mathbf{W}^{(j)}\big{)}$
$\displaystyle=\mathbf{W}_{I}^{(j)}+\mathbf{K}_{I,I}^{-1}\big{(}\mathbf{B}_{I}-\mathbf{K}_{I}\mathbf{W}^{(j)}\big{)}$
$\displaystyle=\mathbf{W}_{I}^{(j)}+\mathbf{K}_{I,I}^{-1}\big{(}\mathbf{B}_{I}-\mathbf{K}_{I,I}\mathbf{W}_{I}^{(j)}-\mathbf{K}_{I,\neg
I}\mathbf{W}_{\neg I}^{(j)}\big{)}$
$\displaystyle=\mathbf{K}_{I,I}^{-1}\big{(}\mathbf{B}_{I}-\mathbf{K}_{I,\neg
I}\mathbf{W}_{\neg I}^{(j)}\big{)}$

where the second line uses the definition of $\mathbf{E}_{I}$; the third line
split the matrix $\mathbf{K}_{I}$ into blocks
$\mathbf{K}_{I}=\big{(}\begin{matrix}\mathbf{K}_{I,I}&\mathbf{K}_{I,\neg
I}\end{matrix}\big{)}$; the last line is straightforward algebra. ∎

##  Appendix B Technical Lemmas

######  Lemma 3.

The quadratic objective function (12) satisfies the Polyak-Łojasiewicz (PL)
inequality

$\displaystyle\frac{1}{2}\|\nabla
h(\mathbf{W})\|_{\mathrm{F}}^{2}\geq\lambda_{\min}(h(\mathbf{W})-h(\mathbf{W}^{*})).$

###### Proof.

If $\mathbf{W}$ has only a single column this follows directly from the strong
convexity of the quadratic function. When $\mathbf{W}$ has multiple columns,
$h$ is a separable function across each column. Therefore, $h$ is also
$\lambda_{\min}$ strongly convex which implies the PL inequality. ∎

######  Lemma 4.

For $h(\mathbf{W})$ as in (12), it holds that

$h(\mathbf{W})-h(\mathbf{W}^{*})=\frac{1}{2}\|\mathbf{W}-\mathbf{W}^{*}\|_{\mathbf{K}}^{2}.$

###### Proof.

Using $\mathbf{B}=\mathbf{K}\mathbf{W}^{*}$, we have

$\displaystyle h(\mathbf{W})-h(\mathbf{W}^{*})$
$\displaystyle=\frac{1}{2}\langle\mathbf{W},\mathbf{K}\mathbf{W}\rangle-\langle\mathbf{B},\mathbf{W}\rangle-\frac{1}{2}\langle\mathbf{W}^{*},\mathbf{K}\mathbf{W}^{*}\rangle+\langle\mathbf{B},\mathbf{W}^{*}\rangle$
$\displaystyle=\frac{1}{2}\langle\mathbf{W},\mathbf{K}\mathbf{W}\rangle-\langle\mathbf{K}\mathbf{W}^{*},\mathbf{W}\rangle-\frac{1}{2}\langle\mathbf{W}^{*},\mathbf{K}\mathbf{W}^{*}\rangle+\langle\mathbf{K}\mathbf{W}^{*},\mathbf{W}^{*}\rangle$
$\displaystyle=\frac{1}{2}\langle\mathbf{W},\mathbf{K}\mathbf{W}\rangle-\langle\mathbf{K}\mathbf{W}^{*},\mathbf{W}\rangle+\frac{1}{2}\langle\mathbf{W}^{*},\mathbf{K}\mathbf{W}^{*}\rangle$
$\displaystyle=\frac{1}{2}\|\mathbf{W}-\mathbf{W}^{*}\|_{\mathbf{K}}^{2}.$

∎

##  Appendix C Proof of Theorem 1

See 1

###### Proof.

For any residual $\mathbf{R}$, note the following inequality

$\displaystyle\|\mathbf{R}\|_{\mathrm{F}}^{2}$
$\displaystyle=\sum_{I\in\mathcal{P}}\|\mathbf{R}_{I,:}\|_{\mathrm{F}}^{2}\leq\lvert\mathcal{P}\rvert\cdot\max_{I\in\mathcal{P}}\|\mathbf{R}_{I,:}\|_{\mathrm{F}}^{2}=m\cdot\max_{I\in\mathcal{P}}\|\mathbf{R}_{I,:}\|_{\mathrm{F}}^{2}.$
(13)

The improvement on the objective $h$ as in (12) after the update in the $j$-th
iteration is bounded by

$\displaystyle h(\mathbf{W}^{(j+1)})-h(\mathbf{W}^{(j)})$
$\displaystyle=-\frac{1}{2}\lVert\mathbf{R}_{I,:}^{(j)}\rVert_{\mathbf{K}_{I,I}^{-1}}^{2}$
$\displaystyle\leq-\frac{1}{2\lambda_{\max}^{\prime}}\|\mathbf{R}_{I,:}^{(j)}\|_{\mathrm{F}}^{2}$
$\displaystyle\leq-\frac{1}{2m\lambda_{\max}^{\prime}}\|\mathbf{R}^{(j)}\|_{\mathrm{F}}^{2}$

where the last inequality is due to the Gauss-Southwell selection rule and
(13). Subtract $h^{*}=h(\mathbf{W}^{*})$ from both sides. Then, we have

$\displaystyle h(\mathbf{W}^{(j+1)})-h^{*}$
$\displaystyle=h(\mathbf{W}^{(j)})-h^{*}-\frac{1}{2m\lambda_{\max}^{\prime}}\|\mathbf{R}^{(j)}\|_{\mathrm{F}}^{2}$
$\displaystyle\leq\big{(}1-\frac{\lambda_{\min}}{m\lambda_{\max}^{\prime}}\big{)}\big{(}h(\mathbf{W}^{(j)})-h^{*}\big{)}$
$\displaystyle\leq\big{(}1-\frac{1}{m\kappa^{\prime}}\big{)}\big{(}h(\mathbf{W}^{(j)})-h^{*}\big{)}$

where the second line uses
$\mathbf{R}^{(j)}=\mathbf{B}-\mathbf{K}\mathbf{W}^{(j)}=-\nabla
h(\mathbf{W}^{(j)})$ by Lemma 1 and the PL inequality by Lemma 3. Using the
inequality $(1-x)^{t}\leq\exp(-tx)$, we obtain a convergence rate in the
number of updates $j$

$\displaystyle
h(\mathbf{W}^{(j+1)})-h^{*}\leq\exp\Big{(}-\frac{j}{m\kappa^{\prime}}\Big{)}\big{(}h(\mathbf{W}^{(0)})-h^{*}\big{)}.$

Since each epoch has $m$ updates, the convergence rate in the number of epochs
$t$ is

$\displaystyle
h(\mathbf{W}^{(t+1)})-h^{*}\leq\exp\Big{(}-\frac{t}{\kappa^{\prime}}\Big{)}\big{(}h(\mathbf{W}^{(0)})-h^{*}\big{)}.$

By Lemma 4, the left and right hand sides can be written as
$\|\mathbf{W}^{(t)}-\mathbf{W}^{*}\|_{\mathbf{K}}^{2}$ and
$\|\mathbf{W}^{(0)}-\mathbf{W}^{*}\|_{\mathbf{K}}^{2}$ respectively, which
concludes the proof.

∎

##  Appendix D Descriptions of the UCI Datasets in the Experiments

This section lists the relevant information of the datasets with citations.
The datasets used in the papers are SGEMM GPU (Paredes and Ballester-Ripoll,,
2018), air quality (Chen,, 2019), 3droad (Kaul,, 2013), song (Bertin-Mahieux,,
2011), buzz (Yang et al.,, 2015), house electric (Hebrail and Berard,, 2012),
and gas sensors (Fonollosa,, 2015). All of them are downloaded from the UCI
machine learning repository (Kelly et al.,, 2023).

##  Appendix E Additional Experiments

This section presents more experimental details and additional experiments.

###  E.1 Further Experimental Details

GP Training. All Gaussian processes, including the stochastic variational
Gaussian processes, use an observation noise constraint $\sigma^{2}\geq
10^{-4}$, which is the default in GPyTorch. For the stochastic trace
estimation (2), we use $\ell=15$ random probe vectors. For CG, the probe
vectors are sampled from $\mathcal{N}(\mathbf{0},\mathbf{P})$, where
$\mathbf{P}$ is the pivoted Cholesky preconditioner. Again, these settings are
the default in GPyTorch. For alternating projection, the probe vectors are
sampled from the Rademacher distribution.

Preconditioning. CG uses the pivoted Cholesky preconditioner both in training
and test. During training, the preconditioner size is $500$ on SGEMM, air
quality, 3droad, song and buzz; $300$ on house electric; $150$ on gas sensors.
We decrease the preconditioner size on house electric and gas sensors due to
GPU memory overflow. During test, the preconditioner size is $500$ on SGEMM,
air quality, 3droad, song, buzz and house electric; $300$ on gas sensors.
Again, we decrese the preconditioner size on gas sensors due to GPU memory
flow.

SVGP Training. All SVGPs use $1024$ inducing points and a batch size of
$4096$. On the first six datasets, SVGPs are trained with $50$ iterations of
Adam with a step size $0.01$ and another $150$ iterations of Adam with a step
size $0.001$. On gas sensors, we train the SVGP with $50$ iterations of Adam
with a step size $0.01$ followed by $350$ iterations of Adam with a step size
$0.001$.

The right panel of Figure 2 is produced on with an alternating projection-
trained GP on air quality with batch size $1000$. The linear system solved in
the figure is $\mathbf{K}^{-1}\mathbf{y}$.

Figure 3 is plotted with an alternating projection-trained GP on 3droad. The
linear system in the figure is
$\mathbf{K}^{-1}(\begin{matrix}\mathbf{y}&\mathbf{z}_{1}&\mathbf{z}_{2}&\cdots&\mathbf{z}_{15}\end{matrix})$
where $\mathbf{z}_{i}$ are sampled from a standard Gaussian distribution.

###  E.2 GP Training on House Electric with Large Noise Constraint
$\sigma^{2}\geq 0.1$

We compare Gaussian processes on house electric trained with two different
noise constraints $\sigma^{2}\geq 0.1$ and $\sigma^{2}\geq 10^{-4}$, as shown
Table 2. We observe significant improvements on both RMSE and NLL when the
noise is smaller. In particular, the GP trained with small noise constraint
$\sigma^{2}\geq 10^{-4}$ has $40\%$ smaller RMSE and significantly smaller
NLL. This indicates that artificially inflating the observation noise
$\sigma^{2}$, while making the kernel matrix well-conditioned, ultimately
hurts the predictive performance.

With alternating projection, training the GP with small noise constraint
$\sigma^{2}\geq 10^{-4}$ is as fast as the GP with large noise constraint
$\sigma^{2}\geq 10^{-1}$.

Table 2: Comparison of GP training on the house electric dataset with large noise constraint $\sigma^{2}\geq 0.1$ and small noise constraint $\sigma^{2}\geq 10^{-4}$. Dataset | Method | RMSE | NLL | FLOPs / $2n^{2}$ | Time
---|---|---|---|---|---
house electric $n=2,049,280$ $d=11$ | CG ($\sigma^{2}\geq 10^{-1}$) | $0.050\pm 0.000$ | $-0.196\pm 0.000$ | $1200\pm 8$ |  $9.6$h $\pm 0.6$
Alt. Proj. ($\sigma^{2}\geq 10^{-1}$) | $0.053\pm 0.000$ | $-0.197\pm 0.000$ | $1100\pm 0$ |  $9.8$h $\pm 0.4$
Alt. Proj. ($\sigma^{2}\geq 10^{-4}$) | $\bm{0.030\pm 0.000}$ | $\bm{-1.148\pm 0.001}$ | $1100\pm 0$ |  $9.8$h $\pm 0.4$

###  E.3 CG Iterations During Training

Figure 4 in the main paper is produced on air quality. This section presents
figures on more datasets, as shown in Figure 6. We observe similar phenomenon:
as the noise decreases durining training, the number of CG iteration
increases; in contrast, alternating projection converges steadily.

(a) SGEMM GPU

(b) 3droad

(c) song

(d) buzz

Figure 6:  The observation noise $\sigma^{2}$ and the number of CG
iterations/alternating projection epochs during training. Top: The observation
noisea $\sigma^{2}$ decreases as the training goes. Bottom: CG takes more
iterations to converge as the observation noise decreases during training.
However, alternating projection is less sensitive to the decrease of
observation noise.

###  E.4 Increasing Lanczos Iterations Improves NLL

In the experiment, we use $1000$ Lanczos iterations to compute the predictive
variance and the test negative log likelihood (NLL). This section investigates
the relation between test NLL and the Lanczos iterations, as shown in Figure
7. We empirically observe that increasing the Lanczos iterations always
decreases the test NLL. This suggests that the true NLL of the GPs may be even
lower than what is reported in Section E.6.

(a) SGEMM

(b) air quality

(c) 3droad

(d) song

(e) buzz

(f) house electric

Figure 7:  Test negative log likelihood (NLL) _vs_. the number of Lanczos
iterations. Empirically, the test NLL decreases as the number of Lanczos
iterations increases on all datasets.

###  E.5 Alternating Projection in Test Time

This section compares alternating projection and CG in the test time. We use
CG and alternating projection to compute the predictive mean and the RMSE of
alternating projection-trained GPs, as shown Table 3. CG and alternating
projection report virtually the same RMSE (exactly the same up to the third
digit after the decimal point). However, we observe significant speed up on
large datasets. On the four largest datasets, alternating projection achieves
$2.3\times$ to $72.3\times$ speed up. In particular, the predictive mean on
house electric can be computed in $5$ min with alternating projection.

Table 3:  Compute the predictive mean and the RMSE of the GPs using CG and alternating projection.  Dataset | RMSE | Time |  Speed up
---|---|---|---
CG | Alt. Proj. | CG | Alt. Proj.
SGEMM | $0.046\pm 0.000$ | $0.046\pm 0.000$ |  $35.0$s $\pm 1.1$ |  $13.5$s $\pm 0.3$ | 0.4$\times$
air quality | $0.256\pm 0.001$ | $0.256\pm 0.001$ |  $2.8$s $\pm 0.3$ |  $3.6$s $\pm 0.8$ | 0.7$\times$
3droad | $0.076\pm 0.000$ | $0.076\pm 0.000$ |  $5.8$m $\pm 0.4$ |  $9.6$m $\pm 0.6$ | $0.6\times$
song | $0.749\pm 0.002$ | $0.749\pm 0.001$ |  $38.1$m $\pm 0.7$ |  $16.4$m $\pm 1.0$ | $2.3\times$
buzz | $0.240\pm 0.001$ | $0.239\pm 0.001$ |  $1.2$h $\pm 0.6$ |  $4.4$m $\pm 1.2$ | $17.2\times$
house electric | $0.032\pm 0.000$ | $0.030\pm 0.000$ |  $5.6$h $\pm 0.6$ |  $4.7$m $\pm 0.2$ | $72.3\times$
gas sensors | $0.203$ | $0.203$ |  $16.1$h |  $27.7$m | $34.9\times$

###  E.6 Training Gaussian Processes with Matérn $\nu=1.5$

Lastly, we report results using Matérn $\nu=1.5$. The experimental settings
are exactly the same as Matérn $\nu=2.5$ GPs. We observe similar phenomenon:
while CG-trained GPs and alternating projection-trained GPs have similar RMSE
and NLL, alternating projection achieves $1.4\times$ to $27.2\times$ speed up
against CG.

Table 4:  Gaussian process training on UCI benchmark datasets with Matérn
$\nu=1.5$. Metrics are computed across multiple runs and reported with $\pm$
one standard deviation.

Dataset | Method | RMSE | NLL | FLOPs$/2n^{2}$ | Training time | Speed up
---|---|---|---|---|---|---
SGEMM $n=241,600$ $d=14$ | CG | $\bm{0.048\pm 0.000}$ | $\bm{-1.071\pm 0.001}$ | $550\pm 0$ |  $8.9$m $\pm 0.2$ |
Alt. Proj. | $\bm{0.048\pm 0.000}$ | $-1.060\pm 0.001$ | $550\pm 0$ |  $12.1$m $\pm 0.2$ | $0.7\times$
2-7 | SVGP | $0.085\pm 0.000$ | $-0.932\pm 0.001$ | NA |  $18.3$m $\pm 0.1$ |
air quality $n=382,168$ $d=13$ | CG | $\bm{0.227\pm 0.002}$ | $0.131\pm 0.003$ | $1825\pm 26$ |  $22.5$m $\pm 1.2$ |
Alt. Proj. | $0.253\pm 0.001$ | $\bm{0.033\pm 0.002}$ | $550\pm 0$ |  $16.1$m $\pm 0.5$ | $1.4\times$
2-7 | SVGP | $0.358\pm 0.002$ | $0.387\pm 0.005$ | NA |  $28.8$m $\pm 0.1$ |
3droad $n=434,874$ $d=3$ | CG | $\bm{0.065\pm 0.001}$ | $1.062\pm 0.003$ | $6086\pm 142$ |  $44.4$m $\pm 2.2$ |
Alt. Proj. | $0.069\pm 0.001$ | $0.896\pm 0.002$ | $572\pm 1$ |  $16.5$m $\pm 0.3$ | $2.7\times$
2-7 | SVGP | $0.319\pm 0.002$ | $\bm{0.294\pm 0.007}$ | NA |  $32.4$m $\pm 0.1$ |
song $n=515,345$ $d=90$ | CG | $\bm{0.743\pm 0.001}$ | $1.135\pm 0.003$ | $4393\pm 159$ |  $13.7$h $\pm 0.6$ |
Alt. Proj. | $\bm{0.746\pm 0.002}$ | $\bm{1.129\pm 0.002}$ | $550\pm 0$ |  $2.6$h $\pm 0.0$ | $5.3\times$
2-7 | SVGP | $0.790\pm 0.002$ | $1.184\pm 0.002$ | NA |  $0.6$h $\pm 0.0$ |
buzz $n=583,250$ $d=77$ | CG | $\bm{0.238\pm 0.000}$ | $0.027\pm 0.002$ | $13608\pm 2299$ |  $25.4$h $\pm 4.7$ |
Alt. Proj. | $\bm{0.238\pm 0.001}$ | $\bm{0.002\pm 0.004}$ | $550\pm 0$ |  $1.9$h $\pm 0.1$ | $13.4\times$
2-7 | SVGP | $0.255\pm 0.002$ | $0.049\pm 0.009$ | NA |  $0.7$h $\pm 0.0$ |
house electric $n=2,049,280$ $d=11$ | CG | - | - | - |  $\geqslant 11$d |
Alt. Proj. | $\bm{0.029\pm 0.000}$ | $-1.321\pm 0.000$ | $1100\pm 0$ |  $9.7$h $\pm 0.1$ | $\geqslant 27.2\times$
2-7 | SVGP | $0.048\pm 0.000$ | $\bm{-1.580\pm 0.003}$ | NA |  $2.6$h $\pm 0.0$ |
gas sensors $n=4,178,504$ $d=17$ | CG | - | - | - | - |
Alt. Proj. | $\bm{0.201}$ | $\bm{0.245}^{\dagger}$ | $1100$ |  $42$h∗ |
2-7 | SVGP | $0.311\pm 0.002$ | $0.286\pm 0.004$ | NA |  $10.6$h $\pm 0.1$ |
$\dagger$ : This predictive variance is calculated using only $500$ Lanczos
iterations to save time and avoid numerical instability.
* : Time measured on a A100 GPU.

##  Appendix F FLOPs in Algorithm 1

The following table gives floating point operations (FLOPs) and memory
complexity of Algorithm 1. There is no hidden constant in the leading term.
Throughout, we assume $l\ll n$ and $1\ll b\ll n$.

Table 5: FLOPs Counting in Algorithm 1. Operation | FLOPs | Memory
---|---|---
Cache Cholesky decomposition of $\\{\mathbf{K}_{I,I}:I\in\mathcal{P}\\}$ | $\frac{1}{3}nb^{2}$ | $nb$
GS rule $I=\operatorname*{argmax}_{I\in\mathcal{P}}~{}\lVert\mathbf{R}_{I,:}\rVert_{\mathrm{F}}^{2}$ | $2nl$ | -
$\mathbf{W}_{I}=\mathbf{W}_{I}+\mathbf{K}_{I,I}^{-1}\mathbf{R}_{I}$ | $(b^{2}+b)l$ | -
$\mathbf{R}=\mathbf{R}-\mathbf{K}_{:,I}\mathbf{K}_{I,I}^{-1}\mathbf{R}_{I}$ | $(b^{2}+2nb+n)l$ | $nb$
total FLOPs of a single epoch | $\big{(}(2+\frac{3}{b})n^{2}+(2b+1)n\big{)}l$ | $2nb$