Hyperparameter Optimization for Randomized Algorithms: A Case Study on Random Features

We begin by prescribing the following statistical model for the data $D_{N}$:

$\displaystyle y_{n}=F(x_{n})+\eta_{n},\quad\text{where}\quad\eta_{n}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathsf{N}(0,\sigma^{2})\quad\text{for all}\quad n\in[N].$

(1)

Here $E\coloneqq\{\eta_{n}\}_{n=1}^{N}$ represents noisy observations in the form of independent and identically distributed (i.i.d.) real-valued zero mean Gaussian random variables with common variance $\sigma^{2}$, and $F$ determines the input-output relationship. The ideal situation is the well-specified setting¹¹1The observed data $D_{N}$ may not necessarily follow the likelihood model (1) in reality. For example the noise may not be i.i.d. Gaussian or may appear multiplicatively. From this perspective, $\sigma>0$ may be viewed as a tunable hyperparameter of the GPR algorithm and not fixed in advance by the data. in which the observed $D_{N}$ is actually generated according to (1) for some ground truth function $F=F^{\star}$.

GPR proceeds by imposing additional structure on (1) in the form of a prior probability distribution over the function $F$. That is, we model $F$ as a GP, independent of $E$, given by

$\displaystyle F\sim\mathsf{GP}\bigl{(}\overline{F},K\bigr{)}.$

(2)

This means that $F\colon\mathcal{X}\to\mathbb{R}$ is a random function with mean function $\overline{F}\colon\mathcal{X}\to\mathbb{R}$ and symmetric covariance function $K\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ such that for any $T\in\mathbb{N}$ and any distinct points $Z\coloneqq\{z_{t}\}_{t=1}^{T}\subset\mathcal{X}$, the finite-dimensional marginals $F(Z)\coloneqq(F(z_{1}),F(z_{2}),\ldots,F(z_{T}))^{\top}\in\mathbb{R}^{T}$ are distributed according to

$\displaystyle F(Z)\sim\mathsf{N}\bigl{(}\overline{F}(Z),K(Z,Z)\bigr{)},$

(3)

where the mean vector and positive semidefinite covariance matrix of the multivariate Gaussian distribution in (3) are given by

$\displaystyle\overline{F}(Z)\coloneqq\bigl{(}\overline{F}(z_{1}),\ldots,\overline{F}(z_{T})\bigr{)}^{\top}\in\mathbb{R}^{T}\quad\text{and}\quad K(Z,Z)\coloneqq\bigl{[}K(z_{i},z_{j})\bigr{]}_{i,j\in[N]}\in\mathbb{R}^{T\times T}.$

(4)

Typically the prior mean $\overline{F}$ is set to zero. The choice of covariance function $K$ greatly influences prediction accuracy and is informed by prior knowledge about the problem, such as smoothness, sparsity, and lengthscales.

To solve the original regression problem, GPR simply applies Bayes’ rule by conditioning the prior function $F$ on the observed data $D_{N}$ to obtain the posterior distribution

$\displaystyle F\,|\,D_{N}\sim\mathsf{GP}\bigl{(}\overline{F}^{(N)},K^{(N)}\bigr{)},$

(5)

where the posterior mean function and posterior covariance function are given by²²2The posterior mean formula here is simply kernel ridge regression (Kanagawa et al., 2018) with a specific ridge penalty parameter depending on $\sigma^{2}$ and $N$.

$\displaystyle\begin{split}\overline{F}^{(N)}(x)&=\overline{F}(x)+K(x,X)\bigl{(}K(X,X)+\sigma^{2}I_{N}\bigr{)}^{-1}\bigl{(}Y-\overline{F}(X)\bigr{)}\quad\text{and}\quad\\ K^{(N)}(x,x^{\prime})&=K(x,x^{\prime})-K(x,X)\bigl{(}K(X,X)+\sigma^{2}I_{N}\bigr{)}^{-1}K(X,x^{\prime})\end{split}$

(6)

for any $x\in\mathcal{X}$ and $x^{\prime}\in\mathcal{X}$ (Williams and Rasmussen, 2006). Here $I_{N}\in\mathbb{R}^{N\times N}$ represents the $N$-dimensional identity and

$\displaystyle K(X,x)\coloneqq\bigl{(}K(x_{1},x),\ldots,K(x_{N},x)\bigr{)}^{\top}\in\mathbb{R}^{N}$

(7)

is the vector of cross covariances between the training inputs and point $x\in\mathcal{X}$. We identify $K(x,X)\coloneqq K(X,x)^{\top}\in\mathbb{R}^{1\times N}$ as a row vector. The whole posterior distribution (5) itself is the GPR estimator. The posterior mean is a standard point predictor and can be related to deterministic kernel methods. The posterior covariance enables the quantification of uncertainty in the prediction and beyond. The underlying regression problem is actually regularized by the prescription of the prior distribution over $F$ in the sense that the Bayesian inversion $D_{N}\mapsto F\,|\,D_{N}$ is a stable mapping (Stuart, 2010). The posterior also being Gaussian is due to the fact that both the likelihood (1) and prior (2) were assumed Gaussian.³³3One can also perform GPR with non-Gaussian likelihoods. However, the posterior will no longer be Gaussian in this case and the closed form expressions (6) for the mean and covariance will no longer be valid. More sophisticated approaches to compute the posterior, such as MCMC, would be required.

The closed form expressions for the posterior mean and covariance functions in (6) are what standard GPR software packages implement in practice. For $x\in\mathcal{X}$, these implementations involve solving the linear systems

$$\bigl{(}K(X,X)+\sigma^{2}I_{N}\bigr{)}\alpha=\bigl{(}Y-\overline{F}(X)\bigr{)}\quad\text{and}\quad\bigl{(}K(X,X)+\sigma^{2}I_{N}\bigr{)}\alpha(x)=K(X,x)$$

(8)

for coefficients $\alpha\in\mathbb{R}^{N}$ and $\alpha(x)\in\mathbb{R}^{N}$, respectively. We observe here that the linear algebraic complexity falls into three stages: space (i.e., memory storage), offline (i.e., operations on stored objects), and online (i.e., operations for evaluation at new $x\in\mathcal{X}$). GPR is a nonparametric method that requires access to the full input data $X$ to evaluate the posterior mean and covariance at online cost $\mathcal{O}(N)$. Storing $X$ requires $\mathcal{O}(dN)$ memory. Since $(K(X,X)+\sigma^{2}I_{N})$ is a dense matrix in general, it requires $\mathcal{O}(N^{2})$ storage in space. Often applications require the solution of the linear systems (8) many times with the same system matrix. In this case, we use a Cholesky factorization that has an offline cost of $\mathcal{O}(N^{3})$ operations. Then the cost of obtaining $\alpha$ in (8) is reduced to $\mathcal{O}(N^{2})$. Similarly, for a new input $x\in\mathcal{X}$, we can use the precomputed Cholesky factor to solve for $\alpha(x)$ with cost $\mathcal{O}(N^{2})$. Nevertheless, for problems with large $N$, such costs make GPR infeasible as presented. This motivates randomized approximations that reduce computational complexity.

Several scalable approximations to GPR have been proposed in the literature. These include, but are not limited to, sparse variational approaches based on inducing points (Hensman et al., 2015, 2018; Titsias, 2009), Nyström subsampling (Sun et al., 2015), and RFs (Rahimi and Recht, 2007). We focus on RFs and in particular the method of RFR. We interpret RFR as an approximation to GPR obtained by first replacing the GP prior covariance $K$ with a low-rank approximation $K_{M}$ and then performing exact GPR inference (5) with this low-rank covariance function. We now describe this procedure in detail.

To begin, let $\Theta$ be a set and $\varphi\colon\mathcal{X}\times\Theta\to\mathbb{R}$ be a feature map. Let $\mu$ be a probability distribution over $\Theta$. The choice of pair $(\varphi,\mu)$ defines the particular RF method. From this pair, define a GP prior covariance function $K\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ for $x$ and $x^{\prime}$ in $\mathcal{X}$ by

$$K(x,x^{\prime})\coloneqq\operatorname{\mathbb{E}}_{\theta\sim\mu}\bigl{[}\varphi(x;\theta)\varphi(x^{\prime};\theta)\bigr{]}.$$

(9)

A covariance defined this way is always symmetric and positive semidefinite by construction. Conversely, any given GP covariance function $K$ can be written in the form (9) for some set $\Theta$ and RF pair $(\varphi,\mu)$.⁴⁴4Indeed, by definition $\mathsf{GP}(0,K)$ satisfies $K(x,x^{\prime})=\operatorname{\mathbb{E}}_{F\sim\mathsf{GP}(0,K)}[F(x)F(x^{\prime})]$. To satisfy (9), take $\Theta$ to be a sufficiently smooth function space, $\mu=\mathsf{GP}(0,K)$, and $\varphi(x;\theta)=F(x)$ with $\theta=F\sim\mu$. In practice, many popular covariance functions fit such a form (Rudi and Rosasco, 2017, Supplement E). The functions $\varphi({\,\cdot\,};\theta)$ with $\theta\sim\mu$ are also called RFs and play a role similar to that of feature maps in the kernel methods literature.

So far, the kernel $K$ is completely deterministic. RF methods apply a Monte Carlo approximation to the integral in (9), replacing the full expectation with the empirical average

$$K_{M}(x,x^{\prime})\coloneqq\frac{1}{M}\sum_{m=1}^{M}\varphi(x;\theta_{m})\varphi(x^{\prime};\theta_{m}),\quad\text{where}\quad\theta_{m}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mu,$$

(10)

for $x\in\mathcal{X}$ and $x^{\prime}\in\mathcal{X}$ and some $M\in\mathbb{N}$. This is a random finite rank prior covariance function. RFR is then simply the process of applying Bayes’ rule to the RF prior $\mathsf{GP}(\overline{F},K_{M})$ under the likelihood model (1) to obtain the RF posterior $\mathsf{GP}(\overline{F}_{M}^{(N)},K_{M}^{(N)})$. Here, the RF posterior mean $\overline{F}_{M}^{(N)}$ and RF posterior covariance $K_{M}^{(N)}$ are given by (6) except with all instances of $K$ replaced by $K_{M}$. In terms of sample paths, we write $F_{M}\sim\mathsf{GP}(\overline{F},K_{M})$ for the RF prior and $F_{M}\,|\,D_{N}\sim\mathsf{GP}(\overline{F}_{M}^{(N)},K_{M}^{(N)})$ for the RF posterior. Although it appears that RFR is just GPR with a special choice of prior, from an implementation and practical point of view it is provably beneficial to view RFR as an entirely distinct methodology (Rahimi and Recht, 2008a, b). Indeed, RFR can be interpreted as a finite rank Bayesian linear model in weight space (Bishop and Nasrabadi, 2006, Section 3.3), which simplifies sampling of the posterior distribution. Moreover, in what follows we show that RFR involves different system matrices that are cheaper to invert and coefficients that have different interpretations than those in appearing in GPR. We leave a rigorous analysis of the convergence of RFR to GPR to future work.

We now follow through the linear algebra to define the system matrices that characterize the posterior mean and covariance of RFR. To this end, we adopt the notation

$\displaystyle\begin{split}\Phi_{M}(x)&\coloneqq\bigl{(}\varphi(x;\theta_{1}),\ldots,\varphi(x;\theta_{M})\bigr{)}\in\mathbb{R}^{1\times M}\quad\text{and}\quad\\ \Phi_{M}(X)&\coloneqq\bigl{[}\varphi(x_{n};\theta_{m})\bigr{]}_{n\in[N],m\in[M]}\in\mathbb{R}^{N\times M}\end{split}$

(11)

for $x\in\mathcal{X}$. Thus $K_{M}(x,x^{\prime})=\Phi_{M}(x)\Phi_{M}(x^{\prime})^{\top}/M$. We observe from (6) and (8) that $\overline{F}_{M}^{(N)}(x)=\overline{F}(x)+K_{M}(x,X)\alpha=\overline{F}(x)+M^{-1}\Phi_{M}(x)\beta$, where

$\displaystyle\beta\coloneqq\Phi_{M}(X)^{\top}\alpha\in\mathbb{R}^{M}.$

(12)

This shows that the difference between the posterior mean and the prior mean belongs to the linear span of the RFs $\{\varphi({\,\cdot\,};\theta_{m})\}_{m=1}^{M}$. We see that the random neural network example from Section 1 holds whenever the RFs $\varphi({\,\cdot\,};\theta)$ take the form of a hidden neuron with sampled weights $\theta\sim\mu$ (or combinations thereof). The coefficients $\beta$ represent the final linear output layer weights in this analogy. It remains to derive a computationally tractable equation for $\beta$. To do this, left multiply the first equation in (8) by $\Phi_{M}(X)^{\top}$ to obtain

$\displaystyle\biggl{(}\frac{1}{M}\Phi_{M}(X)^{\top}\Phi_{M}(X)+\sigma^{2}I_{M}\biggr{)}\beta=\Phi_{M}(X)^{\top}\bigl{(}Y-\overline{F}(X)\bigr{)}.$

(13)

Thus, obtaining $\beta$ only requires the inversion of an $M\times M$ matrix instead of an $N\times N$ one. This is advantageous whenever $M\ll N$.

We similarly seek to write the RF posterior covariance function in the span of the RFs. First, we project the second equation in (8) to obtain the parametrized $M\times M$ system

$$\biggl{(}\frac{1}{M}\Phi_{M}(X)^{\top}\Phi_{M}(X)+\sigma^{2}I_{M}\biggr{)}\beta(x^{\prime})=\biggl{(}\frac{1}{M}\Phi_{M}(X)^{\top}\Phi_{M}(X)\biggr{)}\Phi_{M}(x^{\prime})^{\top}.$$

(14)

for $\beta(x^{\prime})\coloneqq\Phi_{M}(X)^{\top}\alpha(x^{\prime})\in\mathbb{R}^{M}$. From the second equation in (6), we deduce that $K_{M}^{(N)}(x,x^{\prime})=K_{M}(x,x^{\prime})-M^{-1}\Phi_{M}(x)\beta(x^{\prime})=M^{-1}\Phi_{M}(x)(\Phi_{M}(x^{\prime})^{\top}-\beta(x^{\prime}))$. However, realizing that the system matrix on the left hand side of (14) is actually the precision matrix of the RF weight space posterior (Bishop and Nasrabadi, 2006, Section 3.3) reveals an even more efficient form of $K_{M}^{(N)}$. This form requires the solution $\widehat{\beta}(x^{\prime})$ of

$$\biggl{(}\frac{1}{M}\Phi_{M}(X)^{\top}\Phi_{M}(X)+\sigma^{2}I_{M}\biggr{)}\widehat{\beta}(x^{\prime})=\Phi_{M}(x^{\prime})^{\top}$$

(15)

and use of the matrix identity $I-(A+\sigma^{2}I)^{-1}A=\sigma^{2}(A+\sigma^{2}I)^{-1}$ for symmetric positive semidefinite $A$. We summarize the closed form formulas for the RF posterior mean and covariance functions in terms of the coefficients $\beta$ from (13) and $\widehat{\beta}(x^{\prime})$ from (15) as

$\displaystyle\begin{split}\overline{F}^{(N)}_{M}(x)&=\overline{F}(x)+\frac{1}{M}\Phi_{M}(x)\beta=\overline{F}(x)+\frac{1}{M}\sum_{m=1}^{M}\beta_{m}\varphi(x;\theta_{m})\quad\text{and}\quad\\ K_{M}^{(N)}(x,x^{\prime})&=\frac{\sigma^{2}}{M}\Phi_{M}(x)\widehat{\beta}(x^{\prime})=\frac{\sigma^{2}}{M}\sum_{m=1}^{M}\bigl{(}\widehat{\beta}(x^{\prime})\bigr{)}_{m}\,\varphi(x;\theta_{m})\end{split}$

(16)

for any $x$ and $x^{\prime}$ in $\mathcal{X}$.

We now assess the computational complexity for RFR. Since RFR is a parametric method of estimation, we can discard the data $X$ from memory once the dense matrix $\Phi_{M}(X)\in\mathbb{R}^{N\times M}$ is stored with cost $\mathcal{O}(NM)$. The new system matrix in (13) leads to offline complexity $\mathcal{O}(M^{3})$ for a Cholesky factorization. Subsequent solves for $\beta$ in (13) offline and $\widehat{\beta}(x)$ in (14) online both have complexity $\mathcal{O}(M^{2})$. The online cost to evaluate the RF posterior mean and covariance functions is $\mathcal{O}(M)$. Thus, RFR is computationally advantageous to GPR if the RF regression error is comparable to that of GPR even with $M\ll N$. This is indeed true. Theoretical analysis establishes that RFR delivers the same asymptotic error rate as GPR with only $M=\mathcal{O}(\sqrt{N})$ RFs (Lanthaler and Nelsen, 2023).

We conclude this subsection with an example of a commonly used RF pair $(\varphi,\mu)$.

More details about random Fourier features are provided in Section 3.4.

We model $F$ with an $\mathbb{R}^{p}$-valued GP prior $F\sim\mathsf{GP}(\overline{F},K)$. The mean function $\overline{F}\colon\mathcal{X}\to\mathbb{R}^{p}$ is vector-valued. More importantly, the covariance function $K\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}^{p\times p}$ is matrix-valued in order to model correlations in the output space. Several examples of matrix- or operator-valued covariance kernels may be found in the work of Brault et al. (2016); Kadri et al. (2016); Micchelli and Pontil (2004, 2005); Nelsen and Stuart (2021).

Under this $\mathbb{R}^{p}$-valued GP prior, the resulting posterior $F\,|\,D_{N}\sim\mathsf{GP}(\overline{F}^{(N)},K^{(N)})$ under the likelihood (19) is also an $\mathbb{R}^{p}$-valued GP with mean and covariance functions

$\displaystyle\begin{split}\overline{F}^{(N)}(x)&=\overline{F}(x)+K(x,X)\bigl{(}K(X,X)+B_{\Sigma})^{-1}\bigl{(}Y-\overline{F}(X)\bigr{)}\quad\text{and}\quad\\ K^{(N)}(x,x^{\prime})&=K(x,x^{\prime})-K(x,X)\bigl{(}K(X,X)+B_{\Sigma}\bigr{)}^{-1}K(X,x^{\prime})\end{split}$

(20)

for any $x\in\mathcal{X}$ and $x^{\prime}\in\mathcal{X}$. These equations are interpreted in block form. In particular, $B_{\Sigma}\coloneqq\mathrm{blockdiag}(\Sigma,\ldots,\Sigma)\in(\mathbb{R}^{p\times p})^{N\times N}$, $K(X,X)\coloneqq\bigl{[}K(x_{i},x_{j})\bigr{]}_{i,j\in[N]}\in(\mathbb{R}^{p\times p})^{N\times N}$, and

$\displaystyle K(X,x^{\prime})\coloneqq\bigl{(}K(x_{1},x^{\prime}),\ldots,K(x_{N},x^{\prime})\bigr{)}^{\top}\in\bigl{(}\mathbb{R}^{p\times p}\bigr{)}^{N}.$

(21)

The factor $K(x,X)=K(X,x)^{\top}\in(\mathbb{R}^{p\times p})^{1\times N}$ is interpreted as a “row vector” with matrix entries. The concatenated data $Y=\{y_{n}\}_{n=1}^{N}$ and prior mean entries $\overline{F}(X)$ are viewed as length $N$ vectors with each entry taking value in $\mathbb{R}^{p}$. The formulas (20) follow from whitening (19) by pre-multiplying by $\Sigma^{-1/2}$ and then applying an existing white noise result (Owhadi, 2022, Theorem 5.3).

Although (20) visually mimics its scalar-valued counterpart (6), the computational effort is greatly magnified in the vector-valued setting. Effectively, the sample size $N$ is replaced by $Np$ in all scalar-valued output space GPR complexity estimates. Indeed, (20) requires solving $Np$ by $Np$ linear systems of equations instead of $N$ by $N$ ones. Table 1 summarizes upper bounds on the space, offline, and online costs of the approach. This poor scaling with output space dimension $p$ has severely limited the practical utility of vector-valued GPR to date. Indeed, existing studies usually only work with separable (or even scalar) matrix-valued kernels, which suffer from a simplistic rank-one tensor structure (Kadri et al., 2016; Owhadi, 2022), or diagonal kernels, which ignore correlations in the output space (Cleary et al., 2021). Vector-valued RFs enable scalable approximate GP inference without neglecting correlations in the output space.

Vector-valued RFR greatly alleviates the cost of learning with matrix-valued covariance kernels. The method is characterized by a vector-valued feature map $\varphi\colon\mathcal{X}\times\Theta\to\mathbb{R}^{p}$ and a probability distribution $\mu$ over $\Theta$. We work with RF kernels of the form

$$K(x,x^{\prime})\coloneqq\operatorname{\mathbb{E}}_{\theta\sim\mu}\Bigl{[}\varphi(x;\theta)\varphi(x^{\prime};\theta)^{\top}\Bigr{]}\quad\text{and}\quad K_{M}(x,x^{\prime})\coloneqq\frac{1}{M}\sum_{m=1}^{M}\varphi(x;\theta_{m})\varphi(x^{\prime};\theta_{m})^{\top},$$

(22)

where $\theta_{m}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mu$ for $m\in[M]$ and $x$ and $x^{\prime}$ belong to $\mathcal{X}$. The vector-valued RFR method proceeds by assigning a RF prior distribution $F_{M}\sim\mathsf{GP}(\overline{F},K_{M})$ to $F$ in (19). The RF posterior is then given by the vector-valued Gaussian process $\mathsf{GP}(\overline{F}_{M}^{(N)},K_{M}^{(N)})$, where

$\displaystyle\begin{split}\overline{F}^{(N)}_{M}(x)&=\overline{F}(x)+\frac{1}{M}\sum_{m=1}^{M}\beta_{m}\varphi(x;\theta_{m})\quad\text{and}\quad\\ K_{M}^{(N)}(x,x^{\prime})&=\frac{1}{M}\sum_{m=1}^{M}\varphi(x;\theta_{m})\bigl{(}\widehat{\beta}(x^{\prime})\bigr{)}_{m}\end{split}$

(23)

for any $x$ and $x^{\prime}$ in $\mathcal{X}$. The coefficients $\beta\in\mathbb{R}^{M}$ and $\widehat{\beta}(x^{\prime})\in\mathbb{R}^{M\times p}$ solve the linear systems

$\displaystyle\begin{split}\biggl{(}\frac{1}{M}\Phi_{M}(X)^{\top}B_{\Sigma}^{-1}\Phi_{M}(X)+I_{M}\biggr{)}\beta&=\Phi_{M}(X)^{\top}B_{\Sigma}^{-1}\bigl{(}Y-\overline{F}(X)\bigr{)}\quad\text{and}\quad\\ \biggl{(}\frac{1}{M}\Phi_{M}(X)^{\top}B_{\Sigma}^{-1}\Phi_{M}(X)+I_{M}\biggr{)}\widehat{\beta}(x^{\prime})&=\Phi_{M}(x^{\prime})^{\top},\end{split}$

(24)

respectively. The term $\Phi_{M}(x^{\prime})\in\mathbb{R}^{p\times M}$ is analogous to the first line of (11). The formulas in (24) are derived using an approach similar to that in Subsection 2.1.2, except now with an additional left multiplication by $B_{\Sigma}^{-1/2}$. More precisely, we define (24) by

$\displaystyle\begin{split}&\sum_{m=1}^{M}\biggl{(}\frac{1}{M}\sum_{n=1}^{N}\varphi(x_{n};\theta_{\ell})^{\top}\Sigma^{-1}\varphi(x_{n};\theta_{m})+\delta_{\ell m}\biggr{)}\beta_{m}=\sum_{n=1}^{N}\varphi(x_{n};\theta_{\ell})^{\top}\Sigma^{-1}\bigl{(}y_{n}-\overline{F}(x_{n})\bigr{)}\quad\text{and}\quad\\ &\sum_{m=1}^{M}\biggl{(}\frac{1}{M}\sum_{n=1}^{N}\varphi(x_{n};\theta_{\ell})^{\top}\Sigma^{-1}\varphi(x_{n};\theta_{m})+\delta_{\ell m}\biggr{)}\bigl{(}\widehat{\beta}(x^{\prime})\bigr{)}_{m}=\varphi(x^{\prime};\theta_{\ell})^{\top},\end{split}$

(25)

where $\ell\in[M]$ and $\delta_{\ell m}=1$ if $\ell=m$ and equals zero otherwise. The $\Sigma$-weighted system Gram matrices in (25) are still only of size $M$ by $M$. These are much cheaper to invert than the vector-valued GPR system matrices whenever $M\ll Np$; this is often the case in practice. Table 1 summarizes the computational tradeoffs. We observe that the complexity of vector-valued GPR depends directly on $Np$, with cubic scaling in $Np$ offline and quadratic scaling online for storage. On the other hand, the complexity of RFR depends instead on $M$ instead of $Np$, with similar scaling in this variable. Thus, RFR is computationally advantageous whenever $M\ll Np$, while still maintaining the same average prediction accuracy as vector-valued GPR (Lanthaler and Nelsen, 2023).