In 15th International Conference on Foundations of Software
Science and Computational Structures, FOSSACS 2012, volume 7213 of Lecture Notes in Computer Science, pages 437–451. Springer, 2012.

[14]
Yuxin Deng, Tom Chothia, Catuscia Palamidessi, and Jun Pang.
Metrics for Action-labelled Quantitative Transition Systems.
Electr. Notes Theor. Comput. Sci., 153(2):79–96, 2006.

[15]
Josee Desharnais, Vineet Gupta, Radha Jagadeesan, and Prakash Panangaden.
Metrics for Labelled Markov Processes.
Theor. Comput. Sci., 318(3):323–354, 2004.

[16]
Josée Desharnais, François Laviolette, and Mathieu Tracol.
Approximate Analysis of Probabilistic Processes: Logic, Simulation
and Games.
In 5th International Conference on the Quantitative Evaluaiton
of Systems, QEST 2008, pages 264–273. IEEE Computer Society, 2008.

[17]
Nathanaël Fijalkow, Bartek Klin, and Prakash Panangaden.
Expressiveness of Probabilistic Modal Logics, Revisited.
In 44th International Colloquium on Automata, Languages, and
Programming, ICALP 2017, volume 80 of LIPIcs, pages 105:1–105:12.
Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2017.

[18]
Hongfei Fu.
Computing Game Metrics on Markov Decision Processes.
In 39th International Colloquium on Automata, Languages, and
Programming, ICALP 2012, volume 7392 of Lecture Notes in Computer
Science, pages 227–238. Springer, 2012.

[19]
Hongfei Fu.
Verifying Probabilistic Systems: New Algorithms and Complexity
PhD thesis, RWTH Aachen, Aachen, Germany, Nov 2014.

[20]
Alessandro Giacalone, Chi-Chang Jou, and Scott A. Smolka.
Algebraic Reasoning for Probabilistic Concurrent Systems.
In Proceedings of the IFIP WG 2.2/2.3 Working Conference on
Programming Concepts and Methods, pages 443–458. North-Holland, 1990.

[21]
Felix Hausdorff.
Grundzüge der Mengenlehre.
Verlag Von Veit & Comp, Leipzig, 1914.

[22]
Bengt Jonsson and Kim G. Larsen.
Specification and Refinement of Probabilistic Processes.
In 6th Annual Symposium on Logic in Computer Science, LICS
1991, pages 266–277. IEEE Computer Society, 1991.

[23]
Brendan Juba.
On the Hardness of Simple Stochastic Games.
Master's thesis, Carnegie Mellon University, Pittsburgh, PA, USA, May

[24]
Leonid Vitalevich Kantorovich.
On the transfer of masses (in Russian).
Doklady Akademii Nauk, 5(5-6):1–4, 1942.
Translated in Management Science, 1958.

[25]
Barbara König and Christina Mika-Michalski.
(Metric) Bisimulation Games and Real-Valued Modal Logics for
In 29th International Conference on Concurrency Theory, CONCUR
2018, volume 118 of LIPIcs, pages 37:1–37:17. Schloss Dagstuhl -
Leibniz-Zentrum fuer Informatik, 2018.

[26]
Peter Kleinschmidt and Heinz Schannath.
A Strongly Polynomial Algorithm for the Transportation Problem.
Math. Program., 68:1–13, 1995.

[27]
Viktor Klee and Christoph Witzgall.
Facets and Vertices of Transportation Polytopes.
In 5th Summer Seminar on the Mathematics of the Decision
Sciences, volume 11 of Lectures in Applied Mathematics, pages
257–282. AMS, 1967.

[28]
Thomas Liggett and Steven A. Lippman.
Stochastic Games with Perfect Information and Time Average Payoff.
SIAM Review, 11(4):604–607, 1969.

[29]
Facundo Mémoli.
Gromov-Wasserstein Distances and the Metric Approach to Object
Foundations of Computational Mathematics, 11(4):417–487, 2011.

[30]
Matteo Mio.
On the Equivalence of Game and Denotational Semantics for the
Probabilistic $\mu$-calculus.
Logical Methods in Computer Science, 8(2), 2012.

[31]
James B. Orlin.
On the Simplex Algorithm for Networks and Generalized Networks.
In Mathematical Programming Essays in Honor of George B. Dantzig
Part I, pages 166–178. Springer Berlin Heidelberg, 1985.

[32]
Gabriel Peyré and Marco Cuturi.
Computational Optimal Transport.
Foundations and Trends in Machine Learning, 11(5-6):355–607,

[33]
Martin L. Puterman.
Markov Decision Processes: Discrete Stochastic Dynamic
John Wiley & Sons, Inc., New York, NY, USA, 1st edition, 1994.

[34]
Alexander Schrijver.
Theory of Linear and Integer Programming.
Wiley-Interscience series in Discrete Mathematics and Optimization.
Wiley, 1999.

[35]
Roberto Segala.
Modeling and Verification of Randomized Distributed Real-time
PhD thesis, Massachusetts Institute of Technology, Cambridge, MA,
USA, 1995.

[36]
Lloyd S. Shapley.
Stochastic Games.
Proceedings of the National Academy of Sciences,
39(10):1095–1100, 1953.

[37]
Roberto Segala and Nancy A. Lynch.
Probabilistic Simulations for Probabilistic Processes.
In 5th International Conference on Concurrency Theory, CONCUR
1994, volume 836 of Lecture Notes in Computer Science, pages 481–496.
Springer, 1994.

[38]
James K. Strayer.
Linear Programming and its Applications.
Undergraduate Texts in Mathematics. Springer-Verlag, New York, NY,
USA, 1989.

[39]
Qiyi Tang.
Computing Probabilistic Bisimilarity Distances.
PhD thesis, York University, Toronto, Canada, August 2018.

[40]
Qiyi Tang and Franck van Breugel.
Computing Probabilistic Bisimilarity Distances via Policy
In 27th International Conference on Concurrency Theory, CONCUR
2016, volume 59 of LIPIcs, pages 22:1–22:15. Schloss Dagstuhl -
Leibniz-Zentrum fuer Informatik, 2016.

[41]
Qiyi Tang and Franck van Breugel.
Deciding Probabilistic Bisimilarity Distance One for Labelled Markov
In 30th International Conference on Computer Aided Verification,
CAV 2018, volume 10981 of Lecture Notes in Computer Science, pages
681–699. Springer, 2018.

[42]
Qiyi Tang and Franck van Breugel.
Deciding Probabilistic Bisimilarity Distance One for Probabilistic
In Sven Schewe and Lijun Zhang, editors, 29th International
Conference on Concurrency Theory, CONCUR 2018, volume 118 of Leibniz
International Proceedings in Informatics (LIPIcs), pages 9:1–9:17. Schloss
Dagstuhl–Leibniz-Zentrum fuer Informatik, 2018.

[43]
Franck van Breugel.
On Behavioural Pseudometrics and Closure Ordinals.
Inf. Process. Lett., 112(19):715–718, 2012.

[44]
Franck van Breugel, Babita Sharma, and James Worrell.
Approximating a Behavioural Pseudometric without Discount for
Probabilistic Systems.
Logical Methods in Computer Science, 4(2), 2008.

[45]
Franck van Breugel and James Worrell.
The Complexity of Computing a Bisimilarity Pseudometric on
Probabilistic Automata.
In Horizons of the Mind. A Tribute to Prakash Panangaden
—Essays Dedicated to Prakash Panangaden on the Occasion of His 60th
Birthday, volume 8464 of Lecture Notes in Computer Science, pages
191–213. Springer, 2014.

[46]
Cédric Villani.
Optimal Transport: Old and New.
Grundlehren der mathematischen Wissenschaften. Springer, 2008.

[47]
Uri Zwick and Mike Paterson.
The complexity of mean payoff games on graphs.
Theoretical Computer Science, 158(1/2):343–359, May 1996.


# Provably Reliable Large-Scale Sampling from Gaussian Processes

Anthony Stephenson
Department of Mathematics
Bristol University
&Robert Allison
Department of Mathematics
Bristol University
&Edward Pyzer-Knapp
IBM Research

###### Abstract

When comparing approximate Gaussian process (GP) models, it can be helpful to
be able to generate data from any GP. If we are interested in how approximate
methods perform at scale, we may wish to generate very large synthetic
datasets to evaluate them. Naïvely doing so would cost
$\mathcal{O}\left(n^{3}\right)$ flops and $\mathcal{O}\left(n^{2}\right)$
memory to generate a size $n$ sample. We demonstrate how to scale such data
generation to large $n$ whilst still providing guarantees that, with high
probability, the sample is indistinguishable from a sample from the desired
GP.

##  1 Introduction

###  1.1 Motivation

In the GP literature, even when an approximate model designed to be highly
scalable is introduced, evaluation is usually done on a small-scale toy
synthetic dataset of low dimension and on a selection of real datasets. When
large synthetic data are used, they are often noisy versions of deterministic
functions rather than genuine samples from a Gaussian process prior. We
believe this leaves a gap in the analysis: careful assessment of performance
at scale under controlled, model-consistent conditions. As such, in this paper
we show how one can generate datasets that grant this option whilst defining
the criteria necessary for these datasets to function as benchmarks. We
acknowledge the work in [1] on sampling from the GP _posterior_ but point out
that it is unhelpful for our purposes.

###  1.2 Sampling from Gaussian Processes

To test GP approximations in a controlled environment we want to be able to
generate reliable synthetic datasets from arbitrary GPs. Since a draw from a
GP evaluated at a finite set of points is distributed according to a
multivariate normal distribution with some mean and covariance $m$ and $C$, we
can form a sample by generating points from a standard Normal distribution and
performing a linear transformation using some decomposition of the covariance
$C=AA^{T}$ i.e. $u\sim\mathcal{N}_{u}\left(0,I_{n}\right)$, $y=Au+m$
$\implies$ $y\sim\mathcal{N}_{y}\left(m,C\right)$. Although we can generate
$u$ efficiently, the decomposition $C=AA^{T}$ is expensive
($\mathcal{O}\left(n^{3}\right)$) using standard methods (SVD or Cholesky).
This quickly becomes infeasible at large $n$, much like GP regression. As a
result, we look to approximate methods to generate synthetic data at scale.

###  1.3 Outline

In this paper we seek to determine how to affordably generate large
(approximate) samples from GPs such that any given sample is
_indistinguishable_ (in some sense) from a sample drawn from the exact GP. We
do so by first reviewing GP approximation techniques (excluding conjugate
gradient based methods, which are considered in D) from the literature that
can be leveraged to generate samples from a GP prior and then computing bounds
on relevant parameters to constrain the error between the approximate and
exact samples (sections 2 and 3). We define this error in terms of the total
variation distance. We go on to define our notion of “indistinguishability” in
section 4 before running experiments to support our claims in section 5.

##  2 Random Fourier Features (RFF)

RFF were introduced as a method of approximating kernels at large scales in
Support Vector Machines and Kernel Ridge Regression problems in [2]. The
method has since been studied extensively, for example in [3, 4, 5, 6, 7, 8].
One of the appealing features of the RFF approximation for sampling from a GP
is that we don’t need to form the full Gram matrix (given by $ZZ^{T}$ with
$Z\in\mathbb{R}^{n\times D}$) to generate samples. To generate samples we need
only construct a single $Z$ matrix and transform a variable
$w\sim\mathcal{N}\left(0,I_{D}\right)$ to get $\hat{f}=Zw$. We use the
formulation suggested in [9] to minimise the variance of our sampling and
adapt their proof technique for our purposes in B.

######  Lemma 2.1.

To generate a sample of size $n$ whose marginal distribution differs from the
true marginal distribution from a given GP by a total variation distance
($\operatorname{\mathcal{TV}}$) of at most $\epsilon$, with probability
$1-\delta$ it is sufficient to use $D$ RFFs, where $D\geq
8\log\left(\frac{n}{\sqrt{\delta}}\right)\frac{n^{2}}{8\epsilon^{2}\sigma_{\xi}^{4}}$
for some $\delta>0$.

This leads to an overall sampling complexity of
$\mathcal{O}\left(nD\right)=\mathcal{O}\left(n^{3}\log n\right)$ which is
_worse_ than what we would get using Cholesky decomposition. However, it is
worth noting that, in terms of memory usage and ease of parallelism, this
method is still competitive since we need only generate a single sample at a
time, computing a single vector inner product per sample. With careful
implementation, memory cost can be as low as $\mathcal{O}\left(1\right)$ (see
B for details).

##  3 Contour Integral Quadrature (CIQ)

There is some literature dedicated to the computation of functions of square
matrices via the approximation of the Cauchy integral formula ([10, 11, 12]).
An algorithm for the function of interest for us, $A^{\frac{1}{2}}$, is
derived in [11] and subsequently built upon by [12] where the authors derive
an efficient quadrature algorithm in this and the $A^{-\frac{1}{2}}$ case,
specifically citing sampling as a potential usage. We make use of their
algorithm in this vein to estimate matrix-vector products of the form
$K^{\frac{1}{2}}u$. We refer the reader to these sources for a thorough
explanation of the method involved, but include a brief summary in C.

In addition to the superior time complexity we show below, this algorithm also
has a modest memory overhead ($\mathcal{O}\left(Qn\right)$ with $Q$ the number
of quadrature points) and general application to _any_ kernel, unlike the RFF
method which is only applicable to stationary kernels and necessitates non-
trivial derivations of Fourier features for non-RBF kernels.

Define parameters $d$ ($x$-dimension), $\theta$ (kernel parameters),
$\sigma_{\xi}^{2}$ (noise-variance) and $\eta$ (weight of noise-variance at
CIQ approximation stage).

  1. 1.

Sample $x$ data from some distribution, e.g.
$x\sim\mathcal{N}_{x}\left(0,\frac{1}{d}I_{d}\right)$.

  2. 2.

Construct partially noisy kernel
$K_{\xi,ij}=k(x_{i},x_{j})+\eta\sigma_{\xi}^{2}\delta_{ij}$

  3. 3.

Sample $u\sim\mathcal{N}_{u}\left(0,I_{n}\right)$.

  4. 4.

Use CIQ to approximate $f\approx\hat{f}=Mu$ where $M\approx
K_{\eta\xi}^{\frac{1}{2}}$.

  5. 5.

Add noise to the sample to get $\hat{y}=\hat{f}+\xi$ with
$\xi\sim\mathcal{N}\left(0,(1-\eta)\sigma_{\xi}^{2}I_{n}\right)$.

Algorithm 1 CIQ Sampling

######  Lemma 3.1.

To sample approximate draws from a Gaussian Process with
$\operatorname{\mathcal{TV}}<\epsilon$ when compared to a draw from the exact
Gaussian Process, $Q$ quadrature points and $J$ Lanczos iterations will be
sufficient provided we use the CIQ procedure from algorithm 1 to generate our
draw, where $Q$ and $J$ satisfy
$Q\geq\mathcal{O}\left(\log\left(\frac{n}{\eta\sigma_{\xi}^{2}}\right)(-\log\delta_{Q})\right)$
and
$J\geq\tilde{\mathcal{O}}\left(\frac{\sqrt{n}}{\sqrt{\eta}\sigma_{\xi}}\log\frac{n}{\sigma_{\xi}(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})}\right)$
with $0<\delta_{Q}<\epsilon\sigma_{\xi}\sqrt{1-\eta}$.

###  3.1 Preconditioning

It is shown in the appendix (C.1) that the number of iterations $J$ depends
primarily on the condition number of the kernel, which implies that we can
reduce $J$ using preconditioning.

######  Lemma 3.2.

Using a rank-$k$ ($k=\sqrt{n}$) Nyström preconditioner on an ($n\times n$)
kernel matrix with noise variance $\eta\sigma_{\xi}^{2}$ and some constant
$\tilde{C}^{\prime}>0$ means that setting
$J\geq
1+\frac{\sqrt{\lambda_{k+1}}n^{3/8}}{\sqrt{\eta}\sigma_{\xi}}\left(\frac{5}{4}\log
n-\log(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})+\tilde{C}^{\prime}\right)$
Lanczos iterations in the CIQ algorithm will satisfy our requirements.

To make 3.2 more useful, we rely on a result from 3.3 that shows that for a
certain class of kernels we can relate the $k+1^{th}$ eigenvalue to $n$, to
get a workable bound on $J$:

######  Lemma 3.3.

For a sufficiently smooth radial (C.2.3) kernel function $k\in L^{2}_{\mu}$,
some constant $c_{1}>0$ and Nyström preconditioner of rank
$\lfloor{\sqrt{n}}\rfloor$ we define a variable $\gamma=\frac{7}{8}\log
n-\frac{c_{1}}{2}n^{1/d}$ to obtain sufficient conditions for $J$ under three
possible scenarios:

  1. i

Moderate $n$ s.t. $\gamma>1$: $J\geq\mathcal{O}\left(n^{7/8}\log n\right)$

  2. ii

Larger $n$ s.t. $\gamma\in(0,1)$: $J\geq\mathcal{O}\left((\log n)^{2}\right)$

  3. iii

$n\rightarrow\infty$ s.t. $\gamma<0$: $J\geq\mathcal{O}\left(1\right)$

##  4 “Indistinguishable” distributions

We now define what we mean by ‘indistinguishable’. Assume that samples are
provided either from the true GP $P_{0}$ with $p=\frac{1}{2}$ or from the
distribution of the approximating method $P_{1}$ with $p=\frac{1}{2}$ (i.e. a
uniform prior on models). Our decision process is to select the model with the
largest (exact) posterior probability. This can be shown to produce the
smallest achievable error rate (A). If the models were completely
indistinguishable then the error rate would be $\frac{1}{2}$. Perfect
indistinguishability is an unachievable goal due to limited compute-resource
so we instead require $\Pr(\mbox{error})$ to be within $\epsilon$ of
$\frac{1}{2}$ for suitably small $\epsilon$.

######  Definition 4.1 ($\epsilon$-indistinguishable).

$P_{0}$ and $P_{1}$ are $\epsilon$-indistinguishable if the above optimal
Bayesian decision process has $\Pr(\mbox{error})\geq\frac{1}{2}-\epsilon$.

######  Lemma 4.2.

$P_{0}$ and $P_{1}$ are $\epsilon$-indistinguishable if
$\operatorname{\mathcal{TV}}(P_{0},P_{1})\leq 2\epsilon$.

When combining 4.2 with 2.1 and 3.1, 3.2 and 3.3 and by setting $\epsilon$, we
can obtain rigorous and stringent guarantees that synthetic data will behave
like exact-GP data during subsequent analysis - in particular for the purpose
of evaluating approximate-GP regression performance. A further justification
of this notion of indistinguishability and its relation to hypothesis testing
is given in A.

##  5 Experiments

The results above provide bounds on fidelity parameters ($D,J$) of the
sampling approximations. We run suboptimal hypothesis tests to demonstrate
where choices of $D,J$ are definitely insufficient to reach
‘indistinguishability’ and to enable a like-for-like comparison between CIQ
and RFF. The experiments we run generate data from a known GP (with an
isotropic RBF kernel) using the approximate sampling procedure being tested.
The data are then “whitened” using the true kernel matrix such that, when the
data is sufficiently close to the true generating GP distribution, the output
will be a vector of $\mathcal{N}\left(0,1\right)$ distributed points. We
therefore test the hypothesis that the data is from the true GP by running a
Cramér-von Mises test for normality on the output at a significance level
$\alpha$ (more detail on the specifics is provided in G).

We generate a series of M datasets of sizes $(2^{m})_{m=1}^{M}$ over a range
of what we consider to be the data hyperparameters; that is the kernel-scale
$\sigma_{f}^{2}$, noise variance $\sigma_{\xi}^{2}$, lengthscale $l$ and
dataset dimension $d$. For each of these setups we run experiments with
varying fidelity parameter to determine the value required for the rejection
rate from the experiments to converge on the type I error rate, $\alpha$,
implying that the data is indistinguishable from the true GP using the
(suboptimal) hypothesis test.

###  5.1 Results

Figure 1a shows the results of our experiments using the RFF method. We see
that for each choice of $n$ and for both lengthscales tested, the rate appears
to have converged to the significance level before the number of RFFs
predicted by 2.1. We suspect that the requirement that all elements of the
difference matrix are bounded (see B) is too stringent and could be relaxed,
on the grounds that there are only $n$, not $n^{2}$, independent elements.

1b and 1c show the results when we use the CIQ method without and with
preconditioning (respectively). As with the RFF experiments, we see
convergence before the theorised bounds, as we should hope. It is clear that
preconditioning improves the rate of convergence, as expected.

(a)

(b)

(c)

Figure 1: Rejection rate convergence with size of fidelity parameter.
Significance level ($\alpha$) is shown by a blue dashed line and the 95% CI
around $\alpha$ (for converged results) is in green. The range of results
obtained from running a Cholesky benchmark is shown by the grey bar. The
fidelity parameter is rescaled on the $x$-axis by the upper bound derived in
the relevant section of this paper. Vertical black dashed line is at $1.0$
indicating where we reach that bound. (a) shows the RFF case with no. RFF, D
and $\bar{D}(n)=n^{2}\log n$. (b) shows the convergence with Lanczos
iterations $J$ and $\bar{J}(n)=\sqrt{n}\log n$. (c) is preconditioned CIQ with
$\bar{J}(n)=n^{3/8}\log n$.

##  6 Discussion and conclusion

We show how to generate approximate samples from _any_ Gaussian Process that,
with high probability, cannot be distinguished from a draw from the assumed
GP. Bounds are derived to ensure that relevant approximation parameters are
chosen to satisfy the requirements on the fidelity of the sample for arbitrary
probabilistic bounds at a cost cheaper than a standard approach. We believe
this work to be of use to researchers aiming to develop GP approximations for
use on large datasets. For practical use, we generally suggest the use of CIQ
over RFF or other common approaches on the basis of the strong theoretical
guarantees we can provide, given some computational budget. We do note,
however, that when memory is the main bottleneck RFF may be a preferable
choice. We provide a table in E summarising the performance of different
algorithms and discuss future work in H.

## Acknowledgments

We would like to thank Hamza Alawiye for his original implementation of the
RFF method; Nick Baskerville and Adam Lee for their advice on aspects of the
GPytorch and CIQ code; IBM Research for supplying iCase funding for Anthony
Stephenson; NCSC for contributing toward Robert Allison’s funding and finally
the reviewers for their constructive comments.

## References

  * [1] James T. Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth.  Efficiently sampling functions from gaussian process posteriors.  CoRR, abs/2002.09309, 2020.
  * [2] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007.
  * [3] Zhu Li, Jean-Francois Ton, Dino Oglic, and Dino Sejdinovic.  Towards a unified analysis of random Fourier features.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 3905–3914. PMLR, 09–15 Jun 2019.
  * [4] Tianbao Yang, Yu-feng Li, Mehrdad Mahdavi, Rong Jin, and Zhi-Hua Zhou.  Nyström method vs random fourier features: A theoretical and empirical comparison.  In F. Pereira, C.J. Burges, L. Bottou, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
  * [5] Bharath Sriperumbudur and Zoltan Szabo.  Optimal rates for random fourier features.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
  * [6] Fanghui Liu, Xiaolin Huang, Yudong Chen, and Johan A. K. Suykens.  Random features for kernel approximation: A survey on algorithms, theory, and beyond.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44:7128–7148, 2022.
  * [7] Francis Bach.  On the equivalence between kernel quadrature rules & random feature expansions.  Journal of Machine Learning Research, 18:1–38, 2017.
  * [8] Krzysztof Choromanski, Mark Rowland, Tamas Sarlos, Vikas Sindhwani, Richard Turner, and Adrian Weller.  The geometry of random features.  In Amos Storkey and Fernando Perez-Cruz, editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1–9. PMLR, 09–11 Apr 2018.
  * [9] Danica J. Sutherland and Jeff Schneider.  On the error of random fourier features.  In Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence, UAI’15, page 862–871, Arlington, Virginia, USA, 2015\. AUAI Press.
  * [10] Philip I Davies and Nicholas J Higham.  Computing $f(A)b$ for matrix functions $f$.  In QCD and numerical analysis III, pages 15–24. Springer, 2005\.
  * [11] Nicholas Hale, Nicholas J Higham, and Lloyd N Trefethen.  Computing $A^{\alpha}$, $\log A$ and related matrix functions by contour integrals.  SIAM Journal on Numerical Analysis, 46(5):2505–2523, 2008.
  * [12] Geoff Pleiss, Martin Jankowiak, David Eriksson, Anil Damle, and Jacob Gardner.  Fast matrix square roots with applications to gaussian processes and bayesian optimization.  Advances in Neural Information Processing Systems, 33:22268–22281, 2020.
  * [13] Christopher K Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [14] J.P. Romano and E. L. Lehmann.  Testing Statistical Hypotheses, volume 47.  3rd edition, 2005.
  * [15] Luc Devroye, Abbas Mehrabian, and Tommy Reddad.  The total variation distance between high-dimensional gaussians with the same mean.  https://arxiv.org/abs/1810.08693, 2018.
  * [16] Eric W. Weisstein.  "chi distribution.".  https://mathworld.wolfram.com/ChiDistribution.html.  MathWorld–A Wolfram Web Resource.
  * [17] A. Laforgia and P. Natalini.  On the asymptotic expansion of a ratio of gamma functions.  Journal of Mathematical Analysis and Applications, 389(2):833–837, 2012.
  * [18] Gil Shabat, Era Choshen, Dvir Ben Or, and Nadav Carmel.  Fast and Accurate Gaussian Kernel Ridge Regression Using Matrix Decompositions for Preconditioning.  "http://arxiv.org/abs/1905.10587", 2019.
  * [19] Mikhail Belkin.  Approximation beats concentration? An approximation view on inference with smooth radial kernels.  https://arxiv.org/pdf/1801.03437, 2018.
  * [20] Colin Parker, Albert Fox.  SAMPLING GAUSSIAN DISTRIBUTIONS IN KRYLOV SPACES WITH CONJUGATE GRADIENTS.  SIAM J. SCI. COMPUT., 34(3):312–334, 2012.
  * [21] Jacob R. Gardner, Geoff Pleiss, David Bindel, Kilian Q. Weinberger, and Andrew Gordon Wilson.  Gpytorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration.  Advances in Neural Information Processing Systems, 2018-Decem(NeurIPS):7576–7586, 2018.

##  Appendix A “Indistinguishable” distributions

We adopt the definition of _total variation distance_ used in [14]:

######  Definition A.1.

(Total Variation Distance (TV))

$\displaystyle\operatorname{\mathcal{TV}}(P_{0},P_{1})$
$\displaystyle=\frac{1}{2}\left\lVert
P_{1}-P_{0}\right\rVert_{1}=\frac{1}{2}\int|p_{1}-p_{0}|d\mu$

where $p_{i}$ is the density of $P_{i}$ with respect to any measure $\mu$
dominating both $P_{0}$ and $P_{1}$, i.e. $p_{i}=\frac{dP_{i}}{d\mu}$.

Using this definition, $\operatorname{\mathcal{TV}}(P_{0},P_{1})\in[0,1]$ for
any measures $P_{0},P_{1}$. Note that we use notation such that
$\operatorname{\mathcal{TV}}(P_{0},P_{1})=\operatorname{\mathcal{TV}}(p_{0},p_{1})$.

The two proofs below draw heavily from the proof of Theorem 13.1.1 from [14]:

###### Proof that the decision process of section 4 achieves the minimum
possible error rate.

A general decision process for our problem is represented by a function
$\phi:\mathbb{R}^{n}\rightarrow\\{0,1\\}$ whereby, if we are presented with a
sample vector $y\ \in\mathbb{R}^{n}$, we select model $P_{\phi(y)}$. Under
this decision process $\Pr(\mbox{error})=\Pr(\mbox{select
}P_{1}\>|\>P_{0})\Pr(P_{0})+\Pr(\mbox{select }P_{0}\>|\>P_{1})\Pr(P_{1})$.
With an assumed uniform prior on models, we have

$\displaystyle 2\Pr(\mbox{error})$ $\displaystyle=\Pr(\mbox{select
}P_{1}\>|\>P_{0})+\Pr(\mbox{select }P_{0}\>|\>P_{1})$
$\displaystyle=\int\phi(x)p_{0}d\mu(x)+\int(1-\phi(x))p_{1}d\mu(x)$
$\displaystyle=1+\int\phi(x)(p_{0}(x)-p_{1}(x))d\mu(x)$
$\displaystyle=1+\int_{R_{+}}\phi(x)f(x)d\mu(x)+\int_{R_{-}}\phi(x)f(x)d\mu(x)+\int_{R_{0}}\phi(x)f(x)d\mu(x)$

where $f(x)=p_{0}(x)-p_{1}(x)$, $R_{+}=\\{x\in\mathbb{R}^{n}:f(x)>0\\}$,
$R_{-}=\\{x\in\mathbb{R}^{n}:f(x)<0\\}$ and
$R_{0}=\\{x\in\mathbb{R}^{n}:f(x)=0\\}$. From this it follows that setting
$\phi^{\star}$ to be 0 on $R_{+}$, 1 on $R_{-}$ and either 1 or 0 on $R_{0}$
minimises the probability of error in the decision process. This choice of
$\phi$ precisely agrees with the decision process described in [14].

∎

###### Proof of 4.2.

From the above proof we see that the probability of error achieved by the
optimal decision process is
$p^{\prime}=\frac{1}{2}\left[1+\int_{R_{-}}(p_{0}(x)-p_{1}(x))d\mu(x)\right]$.
If we interchange the roles of $p_{0}$ and $p_{1}$ we obtain the alternative
representation
$p^{\prime\prime}=\frac{1}{2}\left[1+\int_{R_{+}}(p_{1}(x)-p_{0}(x))d\mu(x)\right]$.
Since the probability of error, $p$, satisfies $p=p^{\prime}=p^{\prime\prime}$
we have $p=\frac{1}{2}[p^{\prime}+p^{\prime\prime}]$ which can be rearranged
to give

$\displaystyle p$
$\displaystyle=\frac{1}{2}\left[1-\frac{1}{2}\int\left|p_{1}(x)-p_{0}(x)\right|d\mu(x)\right]=\frac{1}{2}-\frac{1}{2}\operatorname{\mathcal{TV}}(P_{0},P_{1}).$

Hence if we can set $\operatorname{\mathcal{TV}}(P_{0},P_{1})\leq 2\epsilon$
we obtain a maximin error rate of at least $\frac{1}{2}-\epsilon$ as claimed.
∎

######  Remark A.2.

Note that we can consider the _worst_ decision rule, $\phi^{\dagger}$, to do
the exact opposite of $\phi^{\star}$ and obtain that for any $\phi$:
$\frac{1}{2}-\frac{1}{2}\tau\leq p\leq\frac{1}{2}+\frac{1}{2}\tau$ where
$\tau=\operatorname{\mathcal{TV}}(P_{0},P_{1})$. $p=\frac{1}{2}$ only when
$\tau=0$. Additionally, a best-case rate of 0 (under $\phi^{\star}$) is
achieved when $P_{0}$ and $P_{1}$ are maximally different in the sense that
$\tau=1$. In this regime we correspondingly obtain the worst-case rate (under
$\phi^{\dagger}$) of 1.

###  A.1 Bounding the total variation distance

We could try to directly bound the total variation distance, for example, by
using the work in [15]. However we make use of Pinsker’s inequality

$\operatorname{\mathcal{TV}}(P_{0},P_{1})\leq\sqrt{\frac{1}{2}\operatorname{\mathcal{KL}}(P_{0},P_{1})}$

and bound the $\operatorname{\mathcal{KL}}$ instead.

Rather than bound the KL-divergence between the marginal distributions (A.3),
we seek to bound the _conditional_ KL (A.4), since for Gaussian distributions
this collapses to the 2-norm of the difference of the sample vectors,
$f,\hat{f}$ (A.3). Details justifying this are provided in A.2.

We define the _marginal_ KL divergence as the usual KL divergence, to clearly
distinguish from the _conditional_ KL we define below. Note that we assume
$\operatorname{\mathcal{KL}}(Q,P)=\operatorname{\mathcal{KL}}(q,p)$ where
$q,p$ are the densities of $Q,P$ respectively w.r.t some dominating measure
$\mu$.

######  Definition A.3 (Marginal Kullback-Leibler Divergence).

$\displaystyle\operatorname{\mathcal{KL}}(q,p)$ $\displaystyle=\int
q(y)\log\frac{q(y)}{p(y)}dy.$

###  A.2 Conditional $\operatorname{\mathcal{KL}}$

First, note that $f$ and $\hat{f}$ are correlated RVs that actually depend on
an underlying RV $u\sim\mathcal{N}\left(0,I_{n}\right)$. We then condition
both the approximate and true distributions ($q$ and $p$) on $u$ rather than
$f$ and $\hat{f}$. In particular, $f=K^{\frac{1}{2}}u$,
$\hat{f}=\hat{K}^{\frac{1}{2}}u$.

######  Definition A.4 (Conditional Kullback-Leibler Divergence).

$\displaystyle\operatorname{\mathcal{KL}}(q,p\>|\>u)$ $\displaystyle=\int
q(y\>|\>u)\log\frac{q(y\>|\>u)}{p(y\>|\>u)}dy.$

######  Lemma A.5.

Using A.4,

$\displaystyle\operatorname{E}_{u}[\operatorname{\mathcal{KL}}(q,p\>|\>u)]$
$\displaystyle\geq\operatorname{\mathcal{KL}}(q,p)$

###### Proof of A.5.

$\displaystyle\operatorname{E}_{u}[\operatorname{\mathcal{KL}}(q,p\>|\>u)]$
$\displaystyle=\int\pi(u)q(y\>|\>u)\log\frac{q(y\>|\>u)}{p(y\>|\>u)}dydu$
$\displaystyle=\int q(y,u)(\log q(y\>|\>u)-\log p(y\>|\>u))dydu$
$\displaystyle=\int
q(y,u)\left[\log\frac{q(y,u)}{\pi(u)}-\log\frac{p(y,u)}{\pi(u)}\right]dydu$
$\displaystyle=\underbrace{\int
q(y,u)\log\frac{q(y,u)}{p(y,u)}dydu}_{\operatorname{\mathcal{KL}}(q(y,u),p(y,u))}$
$\displaystyle=\int
q(u\>|\>y)q(y)\left[\log\frac{q(y)}{p(y)}+\log\frac{q(u\>|\>y)}{p(u\>|\>y)}\right]dydu$
$\displaystyle=\underbrace{\int
q(y)\log\frac{q(y)}{p(y)}dy}_{\operatorname{\mathcal{KL}}(q,p)}+\int
q(y)\underbrace{\int
q(u\>|\>y)\log\frac{q(u\>|\>y)}{p(u\>|\>y)}du}_{\operatorname{\mathcal{KL}}(q(u\>|\>y),p(u\>|\>y))}dy$
$\displaystyle=\operatorname{\mathcal{KL}}(q,p)+\underbrace{q(y^{\prime})\operatorname{\mathcal{KL}}(q(u\>|\>y),p(u\>|\>y))}_{\geq
0}.$

where we have used Tonelli’s theorem to reorder the integral, the non-
negativity of the KL divergence and subsequently the Mean Value Theorem to
complete the proof. ∎

###  A.3 Gaussians

Since we are specifically interested in samples from GPs, we are fortunate
enough to be dealing with Gaussian distributions and, as such, we can write
down a closed form for the KL-divergence between two Gaussian measures.

For the marginal case, with $p(y)=\mathcal{N}_{y}\left(0,K_{\xi}\right)$ and
$q(y)=\mathcal{N}_{y}\left(0,\hat{K}_{\xi}\right)$, we have

$\displaystyle\operatorname{\mathcal{KL}}(q,p)$
$\displaystyle=\frac{1}{2}\left\\{\operatorname{Tr}K_{\xi}^{-1}\hat{K}_{\xi}-n+\log|K_{\xi}|-\log|\hat{K}_{\xi}|\right\\}.$

We now define $E=\hat{K}_{\xi}-K_{\xi}$ and set $\Delta=K_{\xi}^{-1}E$. Note
that $E$ is symmetric but not p.s.d. So we have:

######  Lemma A.6.

Using the definitions above,

$\displaystyle\operatorname{\mathcal{KL}}(q,p)$
$\displaystyle\leq\frac{\left\lVert
E\right\rVert_{F}^{2}}{4\sigma_{\xi}^{4}}.$

###### Proof.

$\displaystyle\operatorname{\mathcal{KL}}(q,p)$
$\displaystyle=\frac{1}{2}\left\\{\operatorname{Tr}\Delta-\log|I+\Delta|\right\\}$
$\displaystyle\leq\frac{1}{4}\operatorname{Tr}\Delta^{2}$
$\displaystyle\leq\frac{1}{4}\left\lVert
K_{\xi}^{-1}\right\rVert_{2}^{2}\left\lVert E\right\rVert_{F}^{2}$
$\displaystyle=\frac{1}{4}\lambda_{n}(K_{\xi})^{-2}\left\lVert
E\right\rVert_{F}^{2}$ $\displaystyle\leq\frac{\left\lVert
E\right\rVert_{F}^{2}}{4\sigma_{\xi}^{4}}$

where we have used that

$\displaystyle\log|I+\Delta|$
$\displaystyle=\sum_{i}\log(1+\lambda_{i}(\Delta))\geq\sum_{i}\lambda_{i}(\Delta)-\frac{1}{2}\lambda_{i}(\Delta)^{2}$
$\displaystyle\geq\operatorname{Tr}\Delta-\frac{1}{2}\operatorname{Tr}\Delta^{2},$

and

$\displaystyle\operatorname{Tr}\Delta^{2}$
$\displaystyle=\operatorname{Tr}[(K_{\xi}^{-\frac{1}{2}}EK_{\xi}^{-\frac{1}{2}})^{2}]$
$\displaystyle=\operatorname{Tr}[(\Lambda_{\xi}^{-\frac{1}{2}}U^{T}EU\Lambda_{\xi}^{-\frac{1}{2}})^{2}]$
$\displaystyle\leq\lambda_{1}(K_{\xi}^{-\frac{1}{2}})^{2}\operatorname{Tr}[EU\Lambda_{\xi}^{-1}U^{T}E]$
$\displaystyle=\lambda_{1}(K_{\xi}^{-1})\operatorname{Tr}[\Lambda_{\xi}^{-1}U^{T}E^{2}U]$
$\displaystyle\leq\lambda_{1}(K_{\xi}^{-1})^{2}\operatorname{Tr}(E^{2})$
$\displaystyle=\left\lVert K_{\xi}^{-1}\right\rVert_{2}^{2}\left\lVert
E\right\rVert_{F}^{2}$

with the eigendecomposition $K_{\xi}=U\Lambda_{\xi}U^{T}$.

∎

For the conditional distributions
$p(y\>|\>f)=\mathcal{N}_{y}\left(f,(1-\eta)\sigma_{\xi}^{2}I_{n}\right)$ and
$q(y\>|\>\hat{f})=\mathcal{N}_{y}\left(\hat{f},(1-\eta)\sigma_{\xi}^{2}I_{n}\right)$,
we get

$\displaystyle\operatorname{\mathcal{KL}}(q,p\>|\>f,\hat{f})$
$\displaystyle=\frac{1}{2}\left[\operatorname{Tr}I_{n}+(\hat{f}-f)^{T}\sigma_{\xi}^{-2}(1-\eta)^{-1}I_{n}(\hat{f}-f)-n\right]$
$\displaystyle=\frac{1}{2(1-\eta)\sigma_{\xi}^{2}}\left\lVert\hat{f}-f\right\rVert_{2}^{2}.$
(*)

We can then apply lemma A.5 to demonstrate that, if we can find an upper bound
on the 2-norm between true and approximate function evaluations on some $x$
data, we will be able to correspondingly bound the KL-divergence between the
marginal distributions of the noise-corrupted samples and finally the TV via
the inequality (A.1).

##  Appendix B Random Fourier Features

Let $M$ denote memory usage for each line.

Define a rule to set a random seed for each $j$ in $1..D$ to ensure the same
vector $\omega_{j}$ is used for each row of $X$.

for  _$i:n$_ do

Sample $(x_{i}^{(1)},\dots,x_{i}^{(d)})\sim\mathbb{P}_{X}$ $\triangleright$
$M=\mathcal{O}\left(d\right)$

$\hat{f}_{i}\leftarrow 0$

for  _$j:D$_ do

Sample $(\omega_{j}^{(1)},\dots,\omega_{j}^{(d)})\sim\mathbb{P}_{\Omega}$
$\triangleright$ $\mathcal{O}\left(1\right)\leq
M\leq\mathcal{O}\left(d^{2}\right)$

Compute $z_{j}(x_{i})=g(x_{i}^{T}\omega_{j})$ $\triangleright$
$M=\mathcal{O}\left(1\right)$

Sample $w_{j}\sim\mathcal{N}\left(0,1\right)$ $\triangleright$
$M=\mathcal{O}\left(1\right)$

$\hat{f}_{i}\leftarrow\hat{f}_{i}+z_{j}(x_{i})w_{j}$ $\triangleright$
$M=\mathcal{O}\left(1\right)$

Algorithm 2 Memory-efficient procedure to generate RFF samples.

Algorithm 2 represents an extreme example of how we can trade-off sequential
time complexity in the RFF procedure to produce an exceptionally memory-
efficient method of at worst $\mathcal{O}\left(d^{2}\right)$ (in the most
general case where a $d\times d$ matrix is required to sample from
$\mathbb{P}_{\Omega}$). We would not recommend this algorithm for a practical
implementation, but rather as an illustration of how far the principle can be
pushed: massively distributing work amongst many processors and opting to
write each sample to disk to avoid storing the full length $n$ vector output.

###### Proof of 2.1.

Define the matrix of differences $E_{ij}=z(x_{i})^{T}z(x_{j})-k(x_{i},x_{j})$.

$\displaystyle|E_{ij}|<\varepsilon/n$
$\displaystyle\implies\underbrace{\sum_{ij}E_{ij}^{2}}_{=\left\lVert
E\right\rVert_{F}^{2}}<\varepsilon^{2}\implies\left\lVert
E\right\rVert_{F}<\varepsilon$

which implies

$\displaystyle\Pr\left(\left|E_{ij}\right|<\frac{\varepsilon}{n}\quad\forall
i,j\right)$
$\displaystyle=\Pr\left(\left|E_{ij}\right|^{2}<\frac{\varepsilon^{2}}{n^{2}}\quad\forall
i,j\right)$
$\displaystyle\leq\Pr\left(\sum_{ij}\left|E_{ij}\right|^{2}<\sum_{ij}\frac{\varepsilon^{2}}{n^{2}}\right)$
$\displaystyle=\Pr\left(\sum_{ij}\left|E_{ij}\right|^{2}<\varepsilon^{2}\right).$

At a particular pair of locations $(x,x^{\prime})$ which correspond to an
element of the kernel matrix, we can use the bounds on the (vector) functions
$z_{j}(x)=\sqrt{\frac{2}{D}}(\sin(\omega_{j}^{T}x),\cos(\omega_{j}^{T}x))^{T}$
(as suggested in [9]) and Hoeffding’s inequality to get a probabilistic tail
bound on the error of an element of the approximate kernel matrix:

$\displaystyle S_{D/2}$
$\displaystyle=\sum_{i=1}^{D/2}z_{i}(x)^{T}z_{i}(x^{\prime})-k(x,x^{\prime})$
$\displaystyle-2/D$ $\displaystyle\leq z_{i}(x)^{T}z_{i}(x^{\prime})\leq 2/D$
$\displaystyle p=\Pr\left(|S_{D/2}|\geq\frac{\varepsilon}{n}\right)$
$\displaystyle\leq
2\exp\left(\frac{-2\varepsilon^{2}}{n^{2}\sum_{i=1}^{D/2}(2/D--
2/D)^{2}}\right)$ $\displaystyle=2\exp(-D\varepsilon^{2}/4n^{2}).$

Now we apply the union bound, assuming that the locations are relatively
uncorrelated. Note that $E$ is symmetric and we can ensure the diagonal
elements are 0 so that we only need to bound $\frac{1}{2}n(n-1)$ elements:

$\displaystyle
q=\Pr\left(\bigcup_{ij}\left\\{\left|E_{ij}\right|\geq\frac{\varepsilon}{n}\right\\}\right)$
$\displaystyle\leq
p\cdot\frac{1}{2}n(n-1)=n(n-1)\exp(-D\varepsilon^{2}/4n^{2})$
$\displaystyle\leq n^{2}\exp\left(-\frac{D\varepsilon^{2}}{4n^{2}}\right)$
$\displaystyle\Pr\left(\left|E_{ij}\right|<\frac{\varepsilon}{n}\quad\forall
i,j\right)$ $\displaystyle=1-q.$

If we now state that we wish to choose a number of RFF $D$ s.t. all elements
of the error matrix are less than $\varepsilon/n$ with probability $1-\delta$
then we can rearrange the final expression above with $q=\delta$ to find:

$\displaystyle D$ $\displaystyle\geq
8\log\left(\frac{n}{\sqrt{\delta}}\right)\frac{n^{2}}{\varepsilon^{2}}$
$\displaystyle\implies\left\lVert E\right\rVert_{F}^{2}<\varepsilon^{2}$
$\displaystyle\implies\operatorname{\mathcal{KL}}<\frac{\varepsilon^{2}}{4\sigma_{\xi}^{4}}\quad(\ref{lem:kl_frob})$
$\displaystyle\implies\operatorname{\mathcal{TV}}<\frac{\varepsilon}{\sqrt{8}\sigma_{\xi}^{2}}\quad\eqref{pinsker}.$

To complete the proof we set $\epsilon=\sqrt{8}\sigma_{\xi}^{2}\varepsilon$ so
that $\operatorname{\mathcal{TV}}<\epsilon$. ∎

##  Appendix C Contour Integral Quadrature

###  C.1 Summary of the CIQ method

Here we give a brief summary of the CIQ method, as discussed in [10], [11] and
[12]. We neglect much of the intricacies, which can be found in the
aforementioned sources. The CIQ method relies on a numerical quadrature
approximation of the matrix version of Cauchy’s integral theorem, for some
square matrix $A$:

$f(A)=\frac{1}{2\pi\it{i}}\int_{\Gamma}f(z)(zI-A)^{-1}dz.$

As usual in complex analysis $\Gamma$ is a closed anticlockwise contour over
which $f$ is analytic.

In our case, we want to use $f(z)=z^{\frac{1}{2}}$ so that

$A^{\frac{1}{2}}=\frac{1}{2\pi\it{i}}\int_{\Gamma}z^{\frac{1}{2}}(zI-A)^{-1}dz.$

Section 4 of [11] pays particular attention to this case and makes a change of
variables $w=z^{\frac{1}{2}}$ to find

$\displaystyle A^{\frac{1}{2}}$
$\displaystyle=\frac{A}{\pi\it{i}}\int_{\Gamma_{w}}(w^{2}I-A)^{-1}dw$
$\displaystyle=\frac{\it{i}A}{\pi}\int_{-\it{i}\infty}^{\it{i}\infty}(w^{2}I-A)^{-1}dw$

This expression is then approximated using a trapezoid rule with $Q$
quadrature points:

$\hat{A}_{Q}^{\frac{1}{2}}=\frac{-2K^{\prime}m^{\frac{1}{2}}A}{\pi
Q}\sum_{q=1}^{Q}(w(t_{q})^{2}I-A)^{-1}\mathrm{cn}(t_{q})\mathrm{dn}(t_{q})$

where $\mathrm{cn},\mathrm{dn}$ are Jacobi elliptic functions in standard
notation.

To compute matrix-vector products, as we intend to, note that Cauchy’s
integral formula can be adapted straightforwardly to

$f(A)b=\frac{1}{2\pi\it{i}}\int_{\Gamma}f(z)(zI-A)^{-1}bdz.$

Although these expressions are sufficient to calculate the desired
approximations, we refer the reader to [12] for significantly more detail on
the practicalities of an efficient implementation.

###  C.2 Proof of bounds for CIQ parameters (3.1)

###### Proof.

From [12] we get the following expression for the error when using CIQ to
approximately sample from a multivariate normal:

$\displaystyle\varepsilon_{J}=\left\lVert
a_{J}-K^{\frac{1}{2}}u\right\rVert_{2}$
$\displaystyle\leq\underbrace{\mathcal{O}\left(\exp\left(-\frac{2Q\pi^{2}}{\log\kappa+3}\right)\right)}_{\varepsilon_{Q}}+\underbrace{\frac{2Q\log(5\sqrt{\kappa})\kappa\sqrt{\lambda_{n}}}{\pi}\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{J-1}}_{B(Q,J)}\left\lVert
u\right\rVert_{2}.$ (C.1)

Clearly the most important factor in the number of iterations, $J$, required
for the algorithm to satisfy some prespecified degree of accuracy
$\varepsilon_{J}$ for a sample of size $n$ depends on the condition number
$\kappa$ of the kernel matrix. Thus we need a bound on the condition number,
which we can find by the following argument.

Using the spectral condition number $\kappa=\frac{\lambda_{1}}{\lambda_{n}}$,
where we have ordered the eigenvalues such that
$\lambda_{1}\geq\lambda_{2}\geq\dots\lambda_{n}$, we first bound the largest
eigenvalue by noting that for a kernel $K$ with kernel-scale
$\sigma_{f}^{2}=1$, $\operatorname{Tr}K=n$ and where
$\operatorname{Tr}K=\sum_{i=1}^{n}\lambda_{i}$, we see that $\lambda_{1}\leq
n$. For a kernel with added uncorrelated noise of the form
$K_{\eta\xi}=K+\eta\sigma_{\xi}^{2}I_{n}$, we similarly see that
$\lambda_{1}(K_{\eta\xi})\leq n+\eta\sigma_{\xi}^{2}$. We include this latter
case as it also allows us to bound the minimum eigenvalue as
$\lambda_{n}(K_{\eta\xi})\geq\eta\sigma_{\xi}^{2}$ and so

$\displaystyle\kappa(K_{\eta\xi})$
$\displaystyle\leq\frac{n}{\eta\sigma_{\xi}^{2}}+1.$

Now, starting from the required bound on the TV-distance we can obtain lower
bounds for the CIQ fidelity parameters $Q,J$ to ensure that $Q$ is such that
$\varepsilon_{Q}<\delta_{Q}$ for some free parameter $\delta_{Q}$ and
$\operatorname{\mathcal{TV}}(q,p)\leq\epsilon$:

$\displaystyle\operatorname{\mathcal{TV}}(q,p)$
$\displaystyle\leq\sqrt{\frac{1}{2}\operatorname{\mathcal{KL}}(q,p)}\leq\frac{1}{\sqrt{2}}\sqrt{\operatorname{E}_{u}\left[\operatorname{\mathcal{KL}}(q,p\>|\>u)\right]}$
$\displaystyle=\frac{1}{2}\sqrt{\operatorname{E}_{u}\left[\frac{1}{(1-\eta)\sigma_{\xi}^{2}}\left\lVert
f-\hat{f}\right\rVert^{2}_{2}\right]}$
$\displaystyle\leq\frac{1}{2\sqrt{1-\eta}\sigma_{\xi}}\sqrt{\operatorname{E}_{u}\left[\varepsilon_{Q}+B(Q,J)\left\lVert
u\right\rVert_{2}^{2}\right]}$
$\displaystyle=\frac{1}{2\sqrt{1-\eta}\sigma_{\xi}}\sqrt{\varepsilon_{Q}^{2}+B(Q,J)^{2}\operatorname{E}_{u}\left\lVert
u\right\rVert_{2}^{2}+\varepsilon_{Q}B(Q,J)\operatorname{E}_{u}\left\lVert
u\right\rVert_{2}}$
$\displaystyle\leq\frac{1}{2\sqrt{1-\eta}\sigma_{\xi}}(\delta_{Q}+\sqrt{n}B(Q,J))\leq\frac{1}{2}\epsilon.$

Where we have used Pinsker’s inequality (A.1), lemma A.5, (A.3), (C.1) and
that $u\sim\mathcal{N}\left(0,I_{n}\right)$ and hence $\left\lVert
u\right\rVert_{2}\sim\chi_{n}$ (note not $\chi^{2}$). Thus we can explicitly
find the mean (and variance) ([16]) and hence a bound on
$\operatorname{E}_{u}\left\lVert u\right\rVert_{2}$. Here we will assume $n$
is large and use the asymptotic expansion in [17] to simplify the result:

$\displaystyle\operatorname{E}[\left\lVert u\right\rVert_{2}]$
$\displaystyle=\sqrt{2}\frac{\Gamma(\frac{1}{2}(n+1))}{\Gamma(n/2)}$
$\displaystyle=\sqrt{n}(1-\frac{1}{4n}+\mathcal{O}\left(n^{-2}\right))\leq\sqrt{n}.$

From here we can write

$\displaystyle\frac{\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q}}{\sqrt{n}}$
$\displaystyle\geq B(Q,J)$ (C.2)

and rearrange for $J$ (once we have bounded $Q$ in terms of our free parameter
$\delta_{Q}$). Note that (C.2) determines an upper bound on $\delta_{Q}$ to
guarantee that the LHS is positive.

In the next sections we make extensive use of the following relationships
(where $\kappa\leq 1+\zeta$ with $\zeta\gg 1$):

$\displaystyle\sqrt{\kappa}\leq\sqrt{1+\zeta}=\sqrt{\zeta}(1+\zeta^{-1})^{\frac{1}{2}}\leq\sqrt{\zeta}+\frac{1}{2\sqrt{\zeta}}$
(C.3)

and

$\displaystyle\log\kappa=\log\zeta+\log(1+\zeta^{-1})\leq\log\zeta+\zeta^{-1}$
(C.4)

to see that

$\displaystyle\sqrt{\kappa}\log\kappa\leq\sqrt{\zeta}\log\zeta+\frac{1}{2\sqrt{\zeta}}\log\zeta+\frac{1}{\sqrt{\zeta}}+\frac{1}{2\zeta^{3/2}}.$
(C.5)

Similarly,

$\displaystyle\log\log\frac{x}{y}$ $\displaystyle=\log(\log x-\log y)=\log\log
x+\log\left(1-\frac{\log y}{\log x}\right)$ $\displaystyle=\log\log
x-\frac{\log y}{\log x}+\mathcal{O}\left((\log x)^{-2}\right)$

giving

$\displaystyle\log\log\frac{x}{y}$ $\displaystyle=\log\log
x+\mathcal{O}\left((\log x)^{-1}\right)$ (C.6)

provided $x\gg y$, and

$\displaystyle\log(\log x+y)$ $\displaystyle=\log\log x+\mathcal{O}\left((\log
x)^{-1}\right)$ (C.7)

if additionally $y<\log x$.

####  C.2.1 Bounding the number of Quadrature Points

To ensure $\varepsilon_{Q}$ does not exceed $\delta_{Q}$ we will constrain $Q$
as follows:

$\displaystyle Q$
$\displaystyle\geq\left(\log\kappa+3\right)(-\log\delta_{Q})\frac{1}{2\pi^{2}}.$
(C.8)

Using (C.4) with $\zeta=\frac{n}{\eta\sigma_{\xi}^{2}}$ we can achieve a
sufficient $Q$ by requiring that

$\displaystyle Q$
$\displaystyle\geq\frac{1}{2\pi^{2}}\left(\log\frac{n}{\eta\sigma_{\xi}^{2}}+3+\mathcal{O}\left(n^{-1}\right)\right)(-\log\delta_{Q}).$
(C.9)

####  C.2.2 Bounding the number of msMINRES Iterations

Taking (C.1) and (C.2) we rearrange in terms of $J$ to find

$\displaystyle J$ $\displaystyle\geq
1+\frac{1}{\log(\sqrt{\kappa}-1)-\log(\sqrt{\kappa}+1)}\log\left\\{\frac{\pi(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})}{2Q\sqrt{\lambda_{n}}\kappa\sqrt{n}(\log(5\sqrt{\kappa}))}\right\\}.$
(C.10)

We start by simplifying the prefactor (making use of Taylor expansions):

$\displaystyle\log(\sqrt{\kappa}-1)-\log(\sqrt{\kappa}+1)$
$\displaystyle=\log(1-1/\sqrt{\kappa})-\log(1+1/\sqrt{\kappa})$
$\displaystyle=-\frac{2}{\sqrt{\kappa}}-\mathcal{O}\left(\kappa^{-3/2}\right)$
$\displaystyle(\log(\sqrt{\kappa}-1)-\log(\sqrt{\kappa}+1))^{-1}$
$\displaystyle=-\frac{\sqrt{\kappa}}{2}\left(1-\mathcal{O}\left(\kappa^{-1}\right)\right).$

Before we substitute our bound for the condition number, we first give a more
general expression. To obtain it we define the RHS of (C.8) to be $\bar{Q}$
and that
$\log\bar{Q}\leq\log\log\kappa+\log(-\log\delta_{Q})+\mathcal{O}\left((\log\kappa)^{-1}\right)$
using (C.7). Hence,

$\displaystyle J$ $\displaystyle\geq
1+\frac{\sqrt{\kappa}}{2}\left[\log(\kappa\sigma_{\xi}\sqrt{n})+2\log\log\kappa-\log(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})+C\right]$
(C.11)

where $C$ is a pseudo-constant that contains constants, decaying functions of
$\kappa$ and similarly ‘negligible’ terms (such as $\log(-\log\delta_{Q})$ and
the small error term
$\left|\sqrt{\lambda_{n}}-\sqrt{\eta}\sigma_{\xi}\right|$). Note that the RHS
here is larger than the RHS of (C.10) so that it is a more conservative bound.

Now using our bounds for $\kappa$, $\sqrt{\kappa}$ (C.3) and $\log\kappa$
(C.4) (with $\zeta=\frac{n}{\eta\sigma_{\xi}^{2}}$) we see that (denoting
$\tilde{J}$ as the RHS of (C.11))

$\displaystyle\tilde{J}$ $\displaystyle\leq
1+\frac{\sqrt{n}}{2\sqrt{\eta}\sigma_{\xi}}\left\\{\log\left(\left(\frac{n}{\eta\sigma_{\xi}^{2}}+1\right)\sigma_{\xi}\sqrt{n}\right)+2\log\log\left(\frac{n}{\eta\sigma_{\xi}^{2}}+1\right)-\log(\pi(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q}))+C^{\prime}\right\\}$
$\displaystyle\leq
1+\frac{\sqrt{n}}{2\sqrt{\eta}\sigma_{\xi}}\left\\{\log\left(\frac{n^{3/2}}{\eta\sigma_{\xi}}\right)+2\log\log\frac{n}{\eta\sigma_{\xi}^{2}}-\log(\pi(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q}))+C^{\prime\prime}\right\\}$
$\displaystyle\leq 1+\frac{\sqrt{n}}{2\sqrt{\eta}\sigma_{\xi}}\left\\{\log
n^{3/2}-\log(\pi(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q}))+2\log\log
n+C^{\prime\prime\prime}\right\\}$
$\displaystyle=\tilde{\mathcal{O}}\left(\frac{\sqrt{n}}{\sqrt{\eta}\sigma_{\xi}}\log\frac{n}{\sigma_{\xi}(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})}\right)$

where we have again used (C.4) and (C.6) and absorbed the constant and
decaying terms into the sequence of ‘constants’
$C^{\prime},C^{\prime\prime},C^{\prime\prime\prime}$.

Having upper bounded $\tilde{J}$, if we now use this as a lower bound for $J$
then we have a sufficient condition to satisfy our TV requirement, i.e.

$\displaystyle J$
$\displaystyle\geq\tilde{\mathcal{O}}\left(\frac{\sqrt{n}}{\sqrt{\eta}\sigma_{\xi}}\log\frac{n}{\sigma_{\xi}(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})}\right)$
(C.12)

where we use $\tilde{\mathcal{O}}\left(\cdot\right)$ to mean ignoring
$\log\log$ terms.

∎

####  C.2.3 Preconditioning

###### Proof of 3.2.

If we take the rank-$k$ Nyström approximation ($\tilde{K}$) as our
preconditioner, with cost $\mathcal{O}\left(Nk^{2}\right)$, then Corollary
4.10 of [18] shows:

$\displaystyle\tilde{\kappa}=\operatorname{cond}[(\tilde{K}+\eta\sigma_{\xi}^{2}I)^{-1}(K+\eta\sigma_{\xi}^{2}I)]$
$\displaystyle\leq
1+\frac{2\lambda_{k+1}(K)\sqrt{4k(n-k)+1}}{\eta\sigma_{\xi}^{2}}.$ (C.13)

(Henceforth we use $\lambda_{k}$ to denote $\lambda_{k}(K)$.) To satisfy our
cost requirement we set $k=\lfloor{\sqrt{n}}\rfloor$ and since

$\displaystyle(4\sqrt{n}(n-\sqrt{n})+1)^{\frac{1}{2}}$
$\displaystyle=(1+4n^{3/2}(1-1/\sqrt{n}))^{\frac{1}{2}}$
$\displaystyle=2n^{3/4}(1-1/\sqrt{n})^{\frac{1}{2}}\left(1+\frac{1}{8}n^{-3/2}(1-1/\sqrt{n})^{-1}+\mathcal{O}\left(n^{-3}\right)\right)$
$\displaystyle=2n^{3/4}\left(1-\frac{1}{2\sqrt{n}}-\frac{1}{8n}-\frac{1}{16n^{3/2}}+\mathcal{O}\left(n^{-2}\right)\right)\left(1+\frac{1}{8n^{3/2}}+\mathcal{O}\left(n^{-2}\right)\right)$
$\displaystyle\leq 2n^{3/4}$

we have

$\displaystyle\tilde{\kappa}$ $\displaystyle\leq
1+\frac{4\lambda_{k+1}}{\eta\sigma_{\xi}^{2}}n^{3/4}$ (C.14)

which we write as $\tilde{\kappa}\leq 1+\zeta$, as before, but now with
$\zeta=\frac{4\lambda_{k+1}}{\eta\sigma_{\xi}^{2}}n^{3/4}$.

With this value of $\zeta$ we again make use of (C.3), (C.4) and (C.5) to show
that

$\displaystyle\sqrt{\tilde{\kappa}}\log\tilde{\kappa}$
$\displaystyle\leq\frac{2\sqrt{\lambda_{k+1}}}{\sqrt{\eta}\sigma_{\xi}}n^{3/8}\log\left(\frac{4\lambda_{k+1}}{\eta\sigma_{\xi}^{2}}n^{3/4}\right)+\mathcal{O}\left(n^{-3/8}\log
n\right).$ (C.15)

Similarly, we have

$\displaystyle\sqrt{\tilde{\kappa}}$
$\displaystyle\leq\sqrt{\zeta}+\frac{1}{2}\zeta^{-\frac{1}{2}}=\frac{2\sqrt{\lambda_{k+1}}}{\sqrt{\eta}\sigma_{\xi}}n^{3/8}+\mathcal{O}\left(n^{-3/8}\right).$
(C.16)

We finish the proof by rewriting (C.11) in the form

$\displaystyle J$ $\displaystyle\geq
1+\frac{\sqrt{\tilde{\kappa}}}{2}\left(\log(\tilde{\kappa}\sqrt{n}\sigma_{\xi})-\log(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})+\tilde{C}\right)$

and inserting (C.15) (updating $\tilde{C}\rightarrow\tilde{C}^{\prime}$ to
absorb additional approximately negligible terms). ∎

###### Proof of 3.3.

[19] shows us that for sufficiently smooth radial kernels, the $k^{th}$ matrix
eigenvalue is given by (e.g. Gaussian, Cauchy)

$\displaystyle\lambda_{k}$ $\displaystyle\lesssim
n\sqrt{\varphi}c_{2}\mathrm{e}^{-c_{1}k^{1/d}}$

for $c_{1},c_{2}>0$, $x\in\mathbb{R}^{d}$, $\varphi=\sup_{x\in\Omega}k(x,x)$.
For us, $\varphi=\sigma_{f}^{2}$ and we set $k=\lfloor{\sqrt{n}}\rfloor$.

We can use this to see

$\displaystyle\sqrt{\lambda_{k+1}}n^{3/8}$
$\displaystyle\leq\sqrt{\lambda_{k}}n^{3/8}\leq\sqrt{c_{2}\sigma_{f}}n^{7/8}\mathrm{e}^{-\frac{c_{1}}{2}n^{1/d}}\leq\sqrt{c_{2}\sigma_{f}}n^{7/8}$

and hence use this with lemma 3.2 (for i)

$\displaystyle J$ $\displaystyle\geq
1+\frac{\sqrt{c_{2}\sigma_{f}}}{\sqrt{\eta}\sigma_{\xi}}n^{7/8}\left(\frac{5}{4}\log
n-\log(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})+\tilde{C}^{\prime}\right)$
$\displaystyle=1+\mathcal{O}\left(n^{7/8}\log n\right).$

For modest $n$ such that $0<\frac{3}{8}\log
n-\frac{c_{1}}{2}n^{\frac{1}{d}}<1$ (ii),

$\displaystyle\exp\left[\frac{7}{8}\log
n-\frac{c_{1}}{2}n^{\frac{1}{d}}\right]$ $\displaystyle\leq 1+\frac{7}{8}\log
n$

which we combine with 3.2 to obtain

$\displaystyle J$ $\displaystyle\geq
1+\frac{\sqrt{c_{2}\sigma_{f}}}{\sqrt{\eta}\sigma_{\xi}}\left[\frac{31}{32}(\log
n)^{2}+\frac{5}{4}\log
n-\log(\epsilon\sigma_{\xi}\sqrt{1-\eta}-\delta_{Q})+\tilde{C}^{\prime\prime}\right]$
$\displaystyle=\mathcal{O}\left((\log n)^{2}\right)$

where we have again absorbed small constant, decaying and $\log\log$ terms
into $\tilde{C}^{\prime\prime}$.

At sufficiently large $n$ (iii), we can write

$\displaystyle\sqrt{\lambda_{k}}n^{3/8}\log n$
$\displaystyle\leq\sqrt{c_{2}\sigma_{f}}\exp\left[\frac{7}{8}\log n+\log\log
n-\frac{c_{1}}{2}n^{\frac{1}{d}}\right]$
$\displaystyle\leq\sqrt{c_{2}\sigma_{f}}$

since for all finite $d$ the $n^{\frac{1}{d}}$ term grows faster in $n$ than
$\log n$ (and definitely $\log\log n$), the exponent will become negative at
large $n$. This gives us

$J\geq 1+\mathcal{O}\left(1\right).$

∎

######  Remark C.1.

Note that the constant terms hidden in the $\mathcal{O}\left(\cdot\right)$
notation are largely consistent for all $\gamma$ and generally less than
$\mathcal{O}\left(\sigma_{\xi}^{-1}\right)$.

######  Remark C.2.

For the special case of the RBF kernel we can exploit the precise expressions
for the eigenvalues to obtain more specific bounds on $J$, but we skip these
details here.

##  Appendix D Conjugate Gradients (CG)

In addition to the sampling methods described in this paper we would also like
to acknowledge that the conjugate gradients algorithm can also be adapted to
facilitate approximate sampling as, for example, in [20]. Since the
computational complexity of such a procedure is in general
$\mathcal{O}\left(n^{2}k\right)$ where $k$ is the number of iterations
employed by the conjugate gradient algorithm, in order for this to be
competitive, we require $k\leq\mathcal{O}\left(\sqrt{n}\log n\right)$. Within
the time-frame of this paper, we have been unable to investigate this option
fully, but we believe that in cases where the Gram matrix can be well-
approximated by a low-rank matrix (e.g. at large lengthscales), CG sampling is
likely to be a promising approach.

##  Appendix E Summary of algorithms

Method | Time | Space
---|---|---
Cholesky | $\mathcal{O}\left(n^{3}\right)$ | $\mathcal{O}\left(n^{2}\right)$
RFF | $\mathcal{O}\left(n^{3}\log n\right)$ | $\mathcal{O}\left(n\right)$
CIQ | $\mathcal{O}\left(n^{5/2}\log n\right)$ | $\mathcal{O}\left(n\log n\right)$
PCIQ | $\mathcal{O}\left(n^{2.375}\log n\right)$ | $\mathcal{O}\left(n\log n\right)$
Table 1: Time and space complexity of competing methods of generating draws
from a GP. P=with preconditioning.

In addition to the algorithms in table 1 we point out that the RFF method is
highly parallelisable and the ‘wall-clock’ time could be significantly reduced
from that listed, with a large enough supply of processors; but it is beyond
the scope of this paper to delve into the details of such an implementation.
Finally, the PCIQ entry does not include the gains observed as $n$ becomes
very large (as outlined in 3.3).

##  Appendix F Implementation

We implemented the RFF sampling procedure in NumPy and made use of the
GPyTorch111github.com/cornellius-gp/gpytorch/ ([21]) library to facilitate
CIQ. We wrapped a SciPy implementation of interpolative decomposition to
integrate with GPyTorch for preconditioning with the CIQ method. These will be
made available on our GitHub (github.com/ant-stephenson/gpsampler) and can be
installed as a Python library.

To run our experiments we made use of the HPC system BluePebble at the
University of Bristol, using nodes with 12 CPUs with 15GB of memory each and
allocating a maximum of 200 hours per experiment. This was sufficient for our
purposes and this much parallel compute was only necessary to facillitate the
running of 1000 repeats per experiment.

##  Appendix G Experiments

We present results from some empirical experiments running hypothesis tests in
section 5. We chose to use hypothesis testing rather than directly implement
the Bayesian decision procedure to avoid considerable coding effort and
compute resource whilst still demonstrating our main points.

A more complete description of the method used is the following:

Input: A set of adjustable experiment parameters $\theta$; $N$ the number of
repeat experiments per parameter set.

Output: Hypothesis test rejection rate, $r$.

$r\leftarrow 0$

for  _$i:N$_ do

1:Generate sample $\hat{y}$ of length $n$ using method $M(\theta)$.

2:Whiten the sample using a Cholesky decomposition of the true kernel matrix,
i.e. $\hat{z}=L^{-1}\hat{y}$ for $K_{\xi}=LL^{T}$.

3:Run a (Cramér-von Mises) hypothesis test to determine whether the whitened
sample $\hat{z}$ is consistent with an i.i.d draw from a standard normal
distribution. Record the test result $t\in\\{0,1\\}$ at a predetermined
significance level $\alpha$.

4:$r\leftarrow r+t$

$r\leftarrow r/N$

Algorithm 3 Procedure used to test sample quality.

In addition to figure 1, figure 2 shows the results of further experiments run
at larger lengthscales to assess how performance of the CIQ method degrades as
condition number becomes more extreme. We note that apart from an increase in
variance (which is expected), the results appear to be quite stable,
oscillating around the blue line and mostly contained between the green lines.
The blue line represents the rejection rate we expect under the null
hypothesis (that the distributions are the same) and is thus the rate we
expect to achieve at convergence. The green lines are given by the 95%
confidence intervals for a large-sample of Bernoulli trials at the converged
rate.

(a)

(b)

Figure 2: Rejection rate convergence with size of fidelity parameter as
before, with additional plots at more extreme lengthscales.

On the results themselves, we note that the PCIQ method appears to converge
more slowly for $l=0.1$ than for $l=1.0$, an initially surprising result. At
first consideration, we expect the kernel matrix for the former to be closer
to the identity and thus have a smaller condition number. This is true, but
conversely, the effectiveness of the preconditioning step is hindered by the
fact that at small lengthscales the matrix will also be full rank and thus to
well-approximate the inverse we are likely to need a higher rank
approximation. To test this hypothesis we ran a simulation of our
(rank-$\sqrt{n}$) preconditioner acting on a series of random kernel matrices
with varying lengthscales. Figure 3 shows the results from this simulation
from which can be seen an apparent peak in the vicinity of $l=0.1$, in
particular for the $n\in\\{2000,5000\\}$ cases, which replicates our previous
findings. It can be shown that the worst-case lengthscale gets smaller as $n$
increases, tending to 0 in the limit as $n\rightarrow\infty$.

Figure 3: Preconditioner effectiveness as a function of lengthscale. Different
colour lines represent different sample sizes. For this simulation we used
parameter values consistent with our previous experiments:
$(d,\sigma_{\xi}^{2},\sigma_{f}^{2})=(2,0.001,1.0)$. The (first) black dashed
line is at $l=0.1$ and the second dot-dash line is at $l=1.0$. $\hat{P}$
represents the preconditioner approximation to the kernel matrix inverse.

##  Appendix H Further work

As noted in 5.1, we believe the bounds on at least the RFF method can be
improved upon. Additionally, more extensive experiments over a wider
hyperparameter space and more general kernel functions, as well as utilising
the Bayesian decision process outlined in the paper, would provide more
thorough support to the arguments. We also acknowledge the vast amount of
literature dedicated to efficient implementations of various linear algebra
routines (e.g. SVD and Cholesky) that could be utilised for our purposes,
albeit with considerable effort to derive similar theoretical guarantees.

In running the experiments we made extensive use of the GPyTorch machinery,
pushing it beyond its intended use; we wish to acknowledge that our
application of CIQ is a ‘misuse’ of the GPyTorch implementation as it was
originally designed. As a result, we believe it would be beneficial to adjust
it to be more compatible with this application, in order to facilitate further
adoption of synthetic data for GP evaluation.