Frequentist coverage and sup-norm convergence rate in Gaussian process
regression
Yun Yang
1

∗1

, Anirban Bhattacharya†2 , and Debdeep Pati‡2

Department of Statistics, Florida State University
Department of Statistics, Texas A&M University

arXiv:1708.04753v1 [] 16 Aug 2017

2

Abstract
Gaussian process (GP) regression is a powerful interpolation technique due to its flexibility in capturing non-linearity. In this paper, we provide a general framework for understanding the frequentist coverage of point-wise and simultaneous Bayesian credible sets
in GP regression. As an intermediate result, we develop a Bernstein von-Mises type result
under supremum norm in random design GP regression. Identifying both the mean and
covariance function of the posterior distribution of the Gaussian process as regularized M estimators, we show that the sampling distribution of the posterior mean function and the
centered posterior distribution can be respectively approximated by two population level
GPs. By developing a comparison inequality between two GPs, we provide exact characterization of frequentist coverage probabilities of Bayesian point-wise credible intervals and
simultaneous credible bands of the regression function. Our results show that inference
based on GP regression tends to be conservative; when the prior is under-smoothed, the
resulting credible intervals and bands have minimax-optimal sizes, with their frequentist
coverage converging to a non-degenerate value between their nominal level and one. As a
byproduct of our theory, we show that the GP regression also yields minimax-optimal posterior contraction rate relative to the supremum norm, which provides a positive evidence
to the long standing problem on optimal supremum norm contraction rate in GP regression.

Key words: Gaussian process regression; Bernstein-von Mises theorem; nonparametric regression; Gaussian comparison theorem; kernel estimator.

1

Introduction and preliminaries

Gaussian process (GP) regression is a popular Bayesian procedure for learning an infinitedimensional function f ∈ F in the nonparametric regression model
ε ∼ N (0, σ 2 ),

y = f (x) + ε,

where x ∈ X ∈ [0, 1]d is the covariate variable and y ∈ Rd the response variable. Through
specifying a GP as the prior measure over the infinite dimensional parameter space F, Bayes
rule yields a posterior measure for f that can be used to either construct a point estimator as the
posterior mean, or characterize estimation uncertainties through the corresponding point-wise
credible intervals and simultaneous credible bands. Examples of the wide usages of GP include
computer experiment emulations [18, 25, 35], spatial data modeling [1, 17], geostatistical kriging
[28, 40], and machine learning [32].
∗

yyang@stat.fsu.edu
anirbanb@stat.tamu.edu
‡
debdeep@stat.tamu.edu
†

1

Despite its long-standing popularity, formal investigations into theoretical properties of GP
regression from a frequentist perspective assuming the existence of a true data-generating function f ∗ is a much more recent activity. A majority of this recent line of work has been directed
towards understanding large sample performance of point estimation in the form of proving posterior consistency and deriving posterior contraction rates. As an incomplete survey, [44, 46]
provide a general framework for studying the rate of posterior contraction of GP regression,
and derive the contraction rates for several commonly used covariance kernels. In [45] and
[4], the authors show that by putting hyper priors over inverse bandwidth parameters in the
square-exponential covariance kernel, GP regression can adapt to the unknown (anisotropic)
smoothness levels of the target function. [49] show that a class of GP priors with Euclideanmetric based covariance kernels can additionally adapt to unknown low-dimensional manifold
structure even with a moderate-to-high dimensional ambient space X in which the covariate
lives.
There is a comparatively limited literature on the validity of conducting statistical inference
in GP regression, or more generally, in Bayesian nonparametric procedures, from a frequentist
perspective. Uncertainty quantification for GPs plays an important role—even more important
than point estimation itself—in some applications such as design of experiments [6] and risk
management [36], and hence it is of interest to investigate the frequentist validity of such.
However, unlike finite dimensional parametric models, where a Bernstein von Mises (BvM)
theorem under the total variation metric holds fairly generally and guarantees asymptotically
nominal coverage of Bayesian credible intervals, the story of the frequentist validity of credible
intervals/bands in infinite dimensional models is more complicated and delicate [16, 20, 24,
27]. For the Gaussian white noise model, [10] showed that a Bernstein von Mises (BvM)
theorem holds in weaker topologies for some natural priors, which yields the correct asymptotic
coverage of credible sets based on the weaker topology; however their result is not applicable
for understanding the coverage of credible intervals/bands. A similar result for the stronger
L∞ -norm can be found in [11]. In the context of linear inverse problems with Gaussian priors,
[26] showed that asymptotic coverage of credible sets can be guaranteed by under-smoothing
the prior compared to the truth. [43] investigated the frequentist coverage of Bayesian credible
sets in the Gaussian white noise model, and showed that depending on whether the smoothness
parameter in the prior distribution matches that in the truth, the coverage of the corresponding
credible sets can either be significantly below its nominal confidence level, or converge to one
as the sample size increases, even though the nominal level is fixed at a constant value. [43]
also investigated the adaptivity of the credible sets to unknown smoothness using an empirical
Bayes approach; see also [33, 37].
A majority of the work discussed in the previous paragraph focuses on the Gaussian white
noise model or its equivalent infinite dimensional normal mean formulation, where appropriate
conjugate Gaussian priors lead to an analytically tractable posterior. A major difficulty of
calculating the nominal coverage probabilities of Bayesian point-wise credible intervals and
simultaneous credible bands in GP regression lies in the fact that the covariance structure
in the posterior GP of f involves unwieldy matrix inverses that are complicated for direct
analysis. Moreover, when the design is random, the randomness in the covariance structure
further complicates the analysis. A relevant work in this regard is [53], where the authors derive
the posterior contraction rates under the supremum norm for GP regression induced by tensor
products of B-splines, and identify conditions under which the coverage of the simultaneous
credible band tends to one as the sample size tends to infinity. However, their result requires
the nominal level of the credible bands to also tend to one, and is based on a key property of
their prior distribution that the resulting posterior covariance function can be sandwiched by
two identity matrices with similar scalings, preventing it to be applicable to a broader class
of GP priors. Also relevant to the present discussion is [39], who obtained similar results for
point-wise credible intervals using a scaled Brownian motion prior. The authors exploit explicit
2

formulas for the eigenvalues and eigenfunctions of the covariance kernel of a Brownian motion
when the design points are on a regular grid, which along with other properties specific to
Brownian motion can be used to linearize the posterior mean and variance, aiding a direct
analysis.
The goal of this article is to provide a general framework for understanding the frequentist
coverage of Bayesian credible intervals and bands in GP regression by proving a BvM type
theorem under the supremum norm in random design nonparametric regression. Towards this
goal, we first show that the sampling distribution of the posterior mean function can be wellapproximated by a GP over the covariate space X with an explicit expression for its covariance
function. Second, we find a tractable population level GP approximation to the centered posterior GP whose covariance function is non-random and also admits an explicit expression. The
frequentist coverage of the Bayesian credible intervals and bands are derived from an interplay
between these two population level GPs. A salient feature of our technique is that it provides an explicit expression of the coverage along with explicit finite sample error bounds—it is
non-asymptotic and applies to any nominal level. Interestingly, we find that when the prior is
under-smoothed, the Bayesian credible intervals and bands are always moderately conservative
in the sense that its frequentist coverage is always higher than its nominal level, and converges
to a fixed number (with explicit expression) between its nominal level and one as the sample size
grows to infinity. For example, when the covariate is one-dimensional and uniformly distributed
over the unit interval, and the prior smoothness level is 2(3), the asymptotic frequentist coverage
of a 95% credible interval is 0.976(0.969). This phenomenon is radically different from existing
results, where the frequentist coverage is either tending to zero or one but never converges to
a non-degenerate value. As a byproduct of our theory, we show that the GP regression also
yields minimax-optimal posterior contraction rate relative to the supremum norm, which provides positive evidence to the long standing problem about supremum norm contraction rate of
general Bayesian nonparametric procedures.
In our proofs, we employ the equivalent kernel representation [32, 38] of the kernel ridge
regression estimator to establish a novel local linear expansion of the posterior mean function
relative to the supremum norm, which is of independent interest and builds a link between
GP regression and frequentist kernel type estimators [21]. This local linear expansion can be
utilized to show the limiting GP approximation to the sampling distribution of the posterior
mean function. Towards the proof of approximating the posterior GP with a population level
GP, we develop a new Gaussian comparison inequality that provides explicit bounds on the
Kolmogorov distance between two GPs through the supremum norm difference between their
respective covariance functions.
Overall, our results reveals delicate interplay between frequentist coverage of Bayesian credible sets and proper characteristics of the prior measure in infinite-dimensional Bayesian procedures. We validate GP regression for conducting statistical inference by showing that as long
as the prior measure is not over-smoothed, Bayesian credible sets always provide moderately
conservative uncertainty quantification with minimax-optimal sizes.
We begin the technical development by introducing notation used throughout the paper. A
summary of all the major notations are provided in Table 2 in §5 for the reader’s convenience.

1.1

Notation

Let H1 , H2 be normed linear spaces. The Fréchet derivative of an operator A at the point
a ∈ H1 is the bounded linear operator denoted DA(a) : H1 → H2 which satisfies
kA(a + h) − A(a) − DA(a)hk
= 0.
h→0
khk
lim

3

(1)

In particular, when H1 = Rn , H2 = Rm , the Fréchet derivative DA(a) is the Jacobian of A, a
linear operator which is represented by an m × n matrix (∂/∂xj Ai ).
Throughout C, C 0 , C1 , C2 , . . . are generically used to denote positive constants whose values
might change from one line to another, but are independent from everything else. We use .
and & denote inequalities up to a constant multiple; we write a  b when both a . b and a & b
hold. For α > 0, let bαc denote the largest integer strictly smaller than α.

1.2

Review of RKHS

We recall some key facts related to reproducing kernel Hilbert spaces (RKHS); further details
and proofs can be found in Chapter 1 of [48]. Let X denote a general index set. A symmetric
function K : X × X → R is positive definite (p.d.) if for any n ∈ N, a1 , . . . , an ∈ R, and
t1 , . . . , tn ∈ X ,
n X
n
X

ai aj K(ti , tj ) > 0.

i=1 j=1

A (real) RKHS H is a Hilbert space of real-valued functions on X such that for any t ∈ X , the
evaluation function Lt : H → R; f →
7 f (t), is a bounded linear functional, i.e., there exists a
constant Mt > 0 such that
|Lt f | = |f (t)| ≤ Mt kf kH , f ∈ H.
p
In the above display, kf kH = hf, f iH is the Hilbert space norm. Since the evaluation maps
are bounded linear, it follows from the Riesz representation theorem that for each t ∈ X , there
exists an element Kt ∈ H such that f (t) = hf, Kt iH . Kt is called the representer of evaluation
at t, and the kernel K(s, t) = Kt (s) is easily shown to be p.d. Conversely, given a p.d. function
K on X × X , one can construct a unique RKHS on X with K as its reproducing kernel. Given
a kernel K, we shall henceforth let Kt (·) = K(·, t).
R
Let L2 (X ) denote the Rspace of square integrable functions f : X → R with X f 2 (x)dx < ∞.
We denote hf, giL2 (X ) := X f (x)g(x)dx the usual inner product on X . If a p.d. kernel K(s, t) is
R R
continuous and X X K 2 (s, t)dsdt < ∞, then by Mercer’s theorem, there exists an orthonormal
2
sequence of continuous eigenfunctions {ψj }∞
j=1 in L (X ) with eigenvalues µ1 ≥ µ2 . . . ≥ 0, and
Z
K(s, t)ψj (s)ds = µj ψj (t), j = 1, 2, . . . ,
K(s, t) =

∞
X

Z Z
µj ψj (s)ψj (t),
X

j=1

K 2 (s, t)dsdt =

X

∞
X

µ2j < ∞.

j=1

The RKHS H determined by K consists of functions f ∈ L2 (X ) satisfying
where fj = hf, ψj iL2 (X ) . Further, for f, g ∈ H,
kf k2H

=

∞
X
fj2
j=1

µj

,

hf, giH =

∞
X
fj gj
j=1

µj

.

P∞

2
j=1 fj /µj

< ∞,

(2)

For stationary kernel K with X = [0, 1] and µ being the Lebesgue measure, we can always
choose ψ2j−1 (s) = sin(πjs) and ψ2j (s) = cos(πjs), j = 1, 2, . . ., by expanding K(s, t) into a
Fourier series over [−1, 1] and exploiting the identity cos(x − y) = cos(x) cos(y) + sin(x) sin(y).
Under these choices, we also have µ2j−1 = µ2j for j = 1, 2, . . . (more details can be found in
Appendix A of [50]).

4

2

Framework

To begin with, we introduce Gaussian process regression and draw its connection (for both
the posterior mean function and posterior covariance function) with kernel ridge regression
(KRR) in §2.1. In §2.2, we introduce the key mathematical tool in our proofs—equivalent
kernel representation of the kernel ridge regression. In §2.3, we present our first result on the
local linear expansion of the KRR estimator relative to the supremum norm, indicating the
asymptotic equivalence between the KRR estimator with a carefully constructed kernel type
estimator.

2.1

GP regression

For easy presentation, we focus on the univariate regression problem where X ⊂ R, and
our results can be straightforwardly extend to multivariate cases. Let (Yi , Xi ), i = 1, . . . , n, be
i.i.d., with Xi ∈ X and Yi ∈ R, with joint density ρY,X (y, x) ∝ φσ (y − f ∗ (x)) 1X (x), where φσ
denotes the N (0, σ 2 ) density, and f ∗ : X → R is the unknown function of interest. Our goal is
to estimate and perform inference on f ∗ based on the data Dn = {(Yi , Xi ), i = 1, . . . , n}. We
assume σ 2 is known throughout the paper.
We consider a nonparametric regression model
Yi = f (Xi ) + εi ,

εi ∼ N (0, σ 2 ),

(3)

and assume a mean-zero GP prior on the regression function f , f ∼ GP(0, σ 2 (nλ)−1 K), where
K is a positive definite kernel and λ > 0 is a tuning parameter to be chosen later.
By conjugacy, it is easy to check that the posterior distribution of f is also a GP, f | Dn ∼
enB (·, ·),
enB ), with mean function fbn (·) and covariance function C
GP(fbn , C
fbn (x) = E[f (x) | Dn ],

enB (x, x0 ) = Cov[f (x), f (x0 ) | Dn ],
C

x, x0 ∈ X .

(4)

e B . These
Since the posterior distribution is a GP, it is completely characterized by fbn and C
n
quantities can be explicitly calculated; we introduce some notation to express these succinctly.
For vectors E = (e1 , . . . , em )T and E0 = (e01 , . . . , e0r )T , let K(E, E0 ) denote the m × r matrix
(K(ei , e0j )). Also, let X = (X1 , . . . , Xn )T and Y = (Y1 , . . . , Yn )T . With these notations,
fbn (x) = K(x, X) [K(X, X) + nλ In ]−1 Y,


B
0
2
−1
0
−1
0
e
Cn (x, x ) = σ (nλ)
K(x, x ) − K(x, X) [K(X, X) + nλ In ] K(X, x ) .

(5)
(6)

In particular, the posterior variance function is


2
−1
−1
Var[f (x) | Dn ] = σ (nλ)
K(x, x) − K(x, X) [K(X, X) + nλ In ] K(X, x) .
The presence of the inverse kernel matrix in equations (5) and (6) renders analysis of the
GP posterior unwieldy. A contribution of this article is to recognize both the mean function and
covariance function as regularized M -estimators and use a equivalent kernel trick to linearize
the solutions that avoids having to deal with matrix inverses. It is well-known (see, e.g. Chapter
6 of [32]) that the posterior mean fbn coincides with the kernel ridge regression (KRR) estimator
fbn,λ = argmin `n,λ (f ),
f ∈H

 X

n
1
2
2
(Yi − f (Xi )) + λ kf kH ,
`n,λ (f ) :=
n
i=1

5

(7)

when the RKHS H corresponds to the prior covariance kernel K. The objective function `n,λ
in (7) combines the average squared-error loss with a squared RKHS norm penalty weighted by
the prior precision parameter λ. It follows from the representer theorem for RKHSs [48] that
the solution to (7) belongs to the linear span of the kernel functions {K(·, Xi )}ni=1 . Given this
fact, solving (7) only amounts to solving a quadratic program, and the solution coincides with
fbn in (5).
A novel observation aiding our subsequent development is that the posterior covariance
enB can be related to the bias of a noise-free KRR estimator. Write σ −2 nλ C
enB (x, x0 )
function C
as
e B (x, x0 ) = Kx (x0 ) − K
b x (x0 ),
σ −2 nλ C
n
b x (·) = K(·, X) [K(X, X) + nλ In ]−1 K(X, x).
where Kx (·) = K(x, ·) as defined earlier, and K
b x is the solution to the following KRR
Comparing with (5) and (7), it becomes apparent that K
problem with noiseless observations of Kx at the random design points {Xi }ni=1 ,

 X
n
1
2
x
2
b x = argmin
(Zi − g(Xi )) + λ kgkH ,
K
n
g∈H

(8)

i=1

where Zix = K(x, Xi ) = Kx (Xi ).
To summarize, the posterior mean corresponds to the usual KRR estimator, and an appropriately scaled version of the posterior covariance function can be recognized as the bias of a
noiseless KRR estimator. This motivates us to study the performance of KRR estimators in
the supremum norm, which to best of our knowledge, hasn’t been carried out before. For past
work on risk bounds for the KRR estimator in other norms, refer to [19, 29, 41, 52, 54].

2.2

Equivalent kernel representation of the KRR estimator

We first introduce an equivalent-kernel formulation that allows us to linearize the bias of a
KRR estimate. Let H ⊂ L2 (X ) be an RKHS, with reproducing kernel K having eigenfunctions
{ψj } and eigenvalues {µj } with respect to L2 (X ). Fix λ > 0. Define a new inner product on H
as
hf, giλ := hf, giL2 (X ) + λ hf, giH ,

(9)

p
and let kf kλ = hf, f iλ . Observe that (H, h·, ·iλ ) is again an RKHS, since for any x ∈ X ,
|f (x)| ≤ C kf kH ≤ C 0 kf kλ for some constant C 0 > 0 depending on C and λ, proving the
boundednessPof the evaluation maps.
P∞
2
Let f = ∞
j=1 fj ψj and g =
j=1 gj ψj be elements of L (X ). Then,
hf, giL2 (X ) + λ hf, giH =

∞
X

fj gj + λ

j=1

∞
X
fj gj
j=1

µj

=

∞
X
fj gj
j=1

νj

,

where
νj =

µj
1
=
,
λ
λ
+
µj
1 + µj

6

j = 1, 2, . . .

(10)


P
P∞ 2
2
Hence, (H, h·, ·iλ ) consists of f = ∞
j=1 fj ψj ∈ L (X ) :
j=1 fj /νj < ∞ , with
hf, giλ =

∞
X
fj gj
j=1

νj

.

(11)

The reproducing kernel associated with the new RKHS is thus
e t) :=
K(s,

∞
X

νj ψj (s)ψj (t).

(12)

j=1

e s (·) denote the representer of the evaluation map, i.e., K
e s (·) = K(s,
e ·), so
As before, we let K
e s iλ . The kernel K
e is known as the equivalent kernel (see,
that for any g ∈ H, g(s) = hg, K
e.g., Chapter 7 of [32]), motivated by the notion of equivalent kernels in the spline smoothing
literature [38]. We comment more on connections with the equivalent kernel literature once we
represent the KRR problem in terms of the equivalent kernel.
Before proceeding further, we introduce two operators which are routinely used subsequently.
First, define a linear operator Fλ : H → H given by
Z
e t)ds.
Fλ g(t) = g(s)K(s,
R
e s ds. The operator Fλ is easily recognized
We shall often use the abbreviation Fλ g = g(s)K
P∞
e
as a convolution
P∞with the equivalent kernel K. If g = j=1 gj ψj , a straightforward calculation
yields Fλ g = j=1 νj gj ψj . Thus, for functions f, g ∈ H, it follows from (11) that
hf, Fλ giλ =

∞
X

fj gj = hf, giL2 (X ) .

j=1

The above display also immediately tells us that Fλ is a self-adjoint operator, i.e., hf, Fλ giλ =
hFλ f, giλ for all f, g ∈ H. Define another self-adjoint operator Pλ : H → H given by Pλ = id−Fλ ,
where id is the identity operator on H. Then, it follows from the previous display and (9) that
hf, Pλ giλ = λ hf, giH ,

f, g ∈ H.

(13)

Having developed the necessary groundwork, let us turn our attention back to the KRR estimate.
e X iλ , and kf k2 =
Recall the objective function `n,λ defined in (7). Writing f (Xi ) = hf, K
i
H
hf, Pλ f iλ using (13), we can express

 X
n
1
2
e X iλ ) + hf, Pλ f i .
(Yi − hf, K
`n,λ (f ) =
i
λ
n
i=1

Viewing `n,λ : (H, h·, ·iλ ) → R and performing a Fréchet differentiation with respect to f , one
obtains a score equation for the KRR estimate fbn,λ as,
Sn,λ (fbn,λ ) = 0,
where Sn,λ : H → H is given by
n

Sn,λ (f ) =

1X
e X − Pλ f.
(Yi − f (Xi ))K
i
n
i=1

7

(14)

Define Sλ (f ) := Ef ∗ [Sn,λ (f )] to be the population version of the score equation, where the
expectation is assumed with respect to the true joint density ρY,X . Recall the convolution
operator Fλ . We then have,
Z
e x dx − Pλ f = Fλ (f ∗ − f ) − Pλ f = Fλ f ∗ − f,
{f ∗ (x) − f (x)} K
(15)
Sλ (f ) =
X

since Fλ f + Pλ f = f by definition. It is immediate that Sλ (Fλ f ∗ ) = 0, which therefore implies
that the function Fλ f ∗ is the solution to the population level score equation. We shall henceforth
refer to Fλ f ∗ as the population-level KRR estimator.
In the above treatment, we first differentiated the objective function and then took an
expectation with respect to the true distribution ρY,X to arrive at the population-level KRR
estimator Fλ f ∗ . One arrives at an identical conclusion if the expectation of the objective
function `n,λ is minimized, which is equivalent to minimizing
kf ∗ − f k2L2 (X ) + λ kf k2H =

∞
X

(fj∗ − fj )2 + λ

j=1

∞
X

fj2 /λj .

j=1

P∞

Solving for fj , one obtains fj = νj fj∗ , and hence f = j=1 νj fj∗ ψj = Fλ f ∗ . This approach thus
also leads to the equivalent kernel; see Chapter 7 of [32] for a detailed exposition along these
lines.

2.3

Sup-norm bounds for the KRR estimator

We now use the equivalent kernel representation to derive error bounds in the supremum
norm between a KRR estimator and its target function. We first lay down two kinds of parameter space considered for the true function f ∗ . Recall the orthonormal basis {ψj }∞
j=1 of
L2 (X ).
For any α > 1/2 and B > 0, define
∞
∞
X
X

ΘαS (B) = f =
fj ψj ∈ L2 (X ) :
j 2α fj2 ≤ B 2 .
j=1

(16)

j=1

For integer α and the Fourier basis, ΘαS corresponds to the α-smooth Sobolev functions with
absolutely continuous (α − 1) derivatives and whose αth derivative has uniformly bounded L2
norm.
Next, for any α > 0 and B > 0, define
∞
∞
X
X

2
ΘH (B) = f =
fj ψj ∈ L (X ) :
j α |fj | ≤ B .
α

j=1

(17)

j=1

The H-subscript is used to indicate a correspondence of this class of functions with α-Hölder
functions. Indeed, under the Fourier basis, if f ∈ ΘαH , then f has bαc continuous derivatives,
and
∞
X
(bαc)
(bαc) 0
f
(x) − f
(x ) .
j bαc |fj | j α−bαc |x − x0 |α−bαc . |x − x0 |α−bαc ,
j=1

which implies that the bαcth derivative of f is Lipschitz continuous of order α − bαc.
We next lay down some standard technical conditions on the eigenfunctions and eigenvalues
of the kernel K.

8

Assumption (B): There exists global constants Cψ , Lψ > 0 such that the eigenfunctions
{ψj }∞
j=1 of K satisfy |ψj (t)| ≤ Cψ for all j ≥ 1, t ∈ X , and |ψj (t) − ψj (s)| ≤ Lψ j|t − s| for all
j ≥ 1 and t, s ∈ X .
−2α for some α > 1/2.
Assumption (E): The eigenvalues {µj }∞
j=1 of K satisfy µj  j

As a motivating example, the Matérn kernel with smoothness index (α − 1/2) satisfies (B)
and (E) when expanded with respect to the Fourier basis
ψ2j−1 (x) = sin(πjx), ψ2j (x) = cos(πjx), j = 1, 2, . . .
on X = [0, 1]; see [3, 50] for more details. Observe that the RKHS associated with a Gaussian
process associated with Matérn kernel with smoothness index (α − 1/2) (in the multivariate
case, α − d/2, with d the dimensionality of X ) is ΘαS . As a passing comment, this space does
not contain functions with smoothness less than α, which includes functions with smoothness
α − 1/2. Although for concreteness we focus on kernels with polynomially decaying eigenvalues
in the paper, our theory is also applicable to other kernel classes, such as squared exponential
kernels and polynomial kernels.
Observe that νj = µj /(µj + λ)  1/(1 + λj 2α ) under Assumption (E). Bounding the sums
by integrals, the following facts are easily observed and used repeatedly in the sequel,
∞
X
j=1

1
. λ−1/(2α) ,
1 + λj 2α

∞
X
j=1

1
. λ−1/(2α) .
(1 + λj 2α )2

(18)

Recall the KRR estimator fbn,λ from (7). Let wi = Yi − f ∗ (Xi ), so that wi ∼ N(0, σ 2 ) are
independent for i = 1, . . . , n
also independent of (X1 , . . . , Xn )T . Recall the operators Fλ
Pand
∞
and Pλ from §2.2. For f = j=1 fj ψj ∈ L2 (X ), we continue to denote
Fλ f =

∞
X

νj fj ψj ,

Pλ f =

j=1

∞
X

(1 − νj )fj ψj .

(19)

j=1

Since νj ∈ (0, 1) for all j ≥ 1, Fλ f and Pλ f are elements of L2 (X ) for any f ∈ L2 (X ).
We now state a theorem which bounds the sup-norm distance between fbn,λ and the true
function f ∗ . To state the theorem in its most general form, we don’t make any smoothness
assumptions of f ∗ yet and state a high probability bound on kfbn,λ − f ∗ k∞ by only assuming
f ∗ ∈ L2 (X ). Reductions of the bound when f ∗ is in either of the smoothness classes ΘαS or ΘαH
are discussed subsequently.
2.1 (Sup-norm bounds for KRR estimator). Assume the kernel K satisfies assumptions (B) and (E). Define h via the equation λ = h2α (α > 1/2). Then, with probability at least
1 − n−10 with respect to the randomness in {(Xi , wi )}ni=1 , the following error estimate holds,
r
log n
kfbn,λ − f ∗ k∞ ≤ 2kPλ f ∗ k∞ + C σ
.
(20)
nh

Theorem

Moreover, with the same probability, the following higher-order expansion holds,
r


n
X
1
log n
∗
0
∗
b
e
fn,λ − Fλ f −
wi KXi
≤ C γn 2kPλ f k∞ + C σ
,
n
nh
∞
i=1

9

(21)

where

n
o r log n
p
−1+1/(2α) −1/(2α)
γn = max 1, n
h
log n
.
(22)
nh
e is the equivalent kernel of K
Here, C and C 0 are constants independent of (n, λ, h, σ), and K
defined in §2.2.
Remark 2.1 (Sub-Gaussian errors). An inspection of the proof of Theorem 2.1 reveals that we
haven’t made explicit use of the normality of the wi s. Indeed, the conclusions of the theorem,
and all subsequent results, extend to sub-Gaussian errors.

The first-order bound in (20) has two components: a bias term, kPλ f ∗ k∞ , and a variance
term, σ (log n)1/2 (nh)−1/2 . Since Pλ f ∗ = f ∗ − Fλ f ∗ , kPλ f ∗ k∞ measures the closeness (in terms
of the supremum norm) between the true function f ∗ and its convolution with the equivalent
e Fλ f ∗ . Recall also that Fλ f ∗ is the solution to the population-level score equation. The
kernel K,
smoothing parameter λ provides the trade-off between the bias and variance; larger λ (stronger
regularization) reduces variance at the cost of increasing the bias, and vice versa. Notice that
a typical analysis of the KRR estimator using basic inequality (such as [51]) requires a lower
bound on the regularization parameter λ, while our analysis is free of such an assumption.
Under additional assumptions on f ∗ , more explicit bounds can be obtained for the bias term
kPλ f ∗ k∞ . For example, if f ∗ ∈ ΘαS (B), then one can show that kPλ f ∗ k∞ . hα−1/2 . Similarly,
for f ∗ ∈ ΘαH (B), one obtains the bound kPλ f ∗ k∞ . hα . Choosing h optimally in either situation
lead to the following explicit bounds.
Corollary

2.1 (Minimax rates for Sobolev and Hölder classes). Suppose f ∗ ∈ ΘαS (B) defined

in (16). Set

h=

B2n
σ 2 log n

−1/(2α)
.

Then, with probability at least 1 − n−10 with respect to the randomness in {(Xi , wi )}ni=1 ,
kfbn,λ − f ∗ k∞ . B 1/(2α)



σ 2 log n
n

 (α−1/2)
2α

.

(23)

Next, suppose f ∗ ∈ ΘαH (B) defined in (17). Set

h=

B2n
σ 2 log n

−1/(2α+1)
.

Then, with probability at least 1 − n−10 with respect to the randomness in {(Xi , wi )}ni=1 ,
kfbn,λ − f ∗ k∞ . B 1/(2α+1)



σ 2 log n
n



α
(2α+1)

.

(24)

Corollary 2.1 implies that the rate of convergence (up to logarithmic terms) of the KRR
estimator to the truth in supremum norm is n−(α−1/2)/(2α) and n−α/(2α+1) for α-smooth Sobolev
and Hölder functions respectively. For functions in the Sobolev class, the KRR estimator
concentrates at the usual minimax rate of n−α/(2α+1) under the L2 norm [8, 54]. However, it
is known [5, 7] that under the point-wise and/or supremum norm, the minimax rate for αsmooth Sobolev functions deteriorates to n−(α−1/2)/(2α) . Hence, the KRR estimator achieves
the minimax rate under the supremum norm as well. For the Hölder class, the minimax rate
remains the same under the L2 and L∞ norms, and the KRR estimator achieves the minimax
rate.
10

The higher-order expansion in the display (21) provides a finer insight
into the distributional
p
P
e X with cobehavior of fbn,λ . Let U (·) denote the zero-mean random process h/n ni=1 wi K
i
0
0
b
b
variance function Cn , so that Cn (x, x ) = EU (x)U (x ). For any fixed x, the law of U (x) can be
bn (x)) distribution by the central limit theorem, and (21) can
approximated by the law of a N(0, C
√
b (x)) distribution.
be used to establish that the law of nh(fbn,λ (x)−Fλ f ∗ (x)) is close to a N(0, C
√ n
Indeed, we shall establish a stronger result that the law of the process nh(fbn,λ (·) − Fλ f ∗ (·))
bn ) in the Kolmogorov distance. The pointcan be approximated by a Gaussian process GP(0, C
wise and uniform approximation results will be crucial to prove asymptotic validity of Bayesian
point-wise and simultaneous credible bands in §3.
Theorem 2.1 specializes to the noiseless KRR problem in a straightforward fashion upon
setting σ = 0 in the bounds (20) and (21). As noted in §2.1, the posterior variance function
can be expressed as the solution to a noiseless KRR problem, motivating our interest in such
situations. The following corollary states a general result for the noiseless case which is used
subsequently to analyze the posterior variance function.
2.2 (Sup-norm bounds for KRR estimator: noiseless case). Consider a noiseless
version of the KRR problem in (7) where Yi = f ∗ (Xi ). Assume the kernel K satisfies assumptions (B) and (E). Define h via the equation λ = h2α . Then, with probability at least 1 − n−10
with respect to the randomness in {Xi }ni=1 , the following error estimate holds,
Corollary

kfbn,λ − f ∗ k∞ ≤ 2kPλ f ∗ k∞ .

(25)

Moreover, with the same probability,
kfbn,λ − Fλ f ∗ k∞ = k(f ∗ − fbn,λ ) − Pλ f ∗ k∞ ≤ C 0 γn kPλ f ∗ k∞ ,

(26)

e is
where γn is as in (22). As before, C and C 0 are constants independent of (n, λ, h), and K
the equivalent kernel of K defined in §2.2.
The proof of Corollary 2.2 follows directly from tracking the proof of Theorem 2.1 with
σ = 0, and hence omitted. Implications of Corollary 2.2 for the posterior variance function are
discussed in the next section.

3

Convergence limit of Bayesian posterior

We now use the sup-norm KRR bounds developed in §2.3 to analyze the GP posterior
enB ) defined in (4).
f | Dn ∼ GP(fbn , C
Implications for the posterior mean function: As noted in §2.1, the posterior mean
fbn under the GP(0, σ 2 (nλ)−1 K) prior coincides with the KRR estimator fbn,λ . Hence, the
conclusions of Theorem 2.1 apply to fbn , that is, with probability at least 1 − n−10 ,
r
log n
∗
∗
,
and
kfbn − f k∞ ≤ 2kPλ f k∞ + C σ
nh
r
 X



(27)
n
1
log n
∗
∗
∗
b
e
fn − f −
wi KXi − Pλ f
≤ C γn kPλ f k∞ + σ
: = δn ,
n
nh
∞
i=1

In particular, for optimal choices of the prior precision parameter λ as in Corollary 2.1,
the posterior mean achieves the minimax rates for the Sobolev and Hölder classes under the
supremum norm.
Implications for the posterior variance function: Recall from § 2.1 that the posterior
11

b x (x0 ), where
enB admits the representation σ −2 nλ C
enB (x, x0 ) = Kx (x0 ) − K
covariance function C
b x is the solution to the noiseless KRR problem in (8). From (31) in Corollary 2.2, it follows
K
with at least 1 − n−10 ,
e B (x, ·)
σ −2 nλ C
n

∞

≤ CkPλ Kx k∞ .

A simple calculation yields
Pλ Kx (·) =

∞
∞
X
X
e x (·).
(1 − νj )µj ψj (x)ψj (·) = λ
νj ψj (x)ψj (·) = λK
j=1

j=1

Since {ψj } are uniformly bounded, we have
∞
X

νj ψj (x)ψj (x0 ) - λ−1/(2α) = h−1 .

(28)

j=1

From (26) in Corollary 2.2, we obtain that with probability at least 1 − n−10 ,
e B (x, ·) − Pλ Kx
σ −2 nλ C
n

∞

≤ C γn kPλ Kx k∞ .

Combining with the previous display,
enB (x, ·) − K
ex
σ −2 n C

∞

e x k∞ .
≤ C γn kK

(29)

e B and
In particular, this relationship between the rescaled posterior covariance function σ −2 n C
n
e leads to a practically useful way to numerically approximate
the equivalent kernel function K
the equivalent kernel function without explicitly conducting the eigen-decomposition to the
original kernel function K.
e B ). Then the conditional distribution of fbn + feB given Dn is the posterior
Let feB ∼ GP(0, C
n
distribution of the mean function f | Dn , or in other words, the law of feB given Dn is that of
the centered posterior distribution. When studying second-order properties
of the posterior, it
√
will be useful to work with a√scaled version
of the posterior with a nh scaling, which has the
√
enB ). Write
same covariance function as nh feB ; nh feB | Dn ∼ GP(0, nh C

e B = σ 2 h σ −2 n C
eB .
nh C
n
n
b B ),
The approximation error bound in (29) motivates us to define a related GP as W B ∼ GP(0, C
n
with
bnB (x, x0 ) : = σ 2 h K(x,
e
C
x0 ) = σ 2 h

∞
X

νj ψj (x) ψj (x0 ).

j=1

Notice that
√ inequality (29) implies that following sup-norm difference between covariance functions of nhfeB and W B ,
e B (x, x0 ) − C
b B (x, x0 ) ≤ C σ 2 h γn sup K(x,
e
sup nh C
x 0 ) . γn ,
n
n
x,x0

e
since supx,x0 |K(x,
x0 )| .

(30)

x,x0

P∞

bnB is a fixed function (de. h−1 from (28). Importantly, C
enB which involves the random design
pending only on the eigenbasis and the scaling h) unlike C
j=1 νj

12

{Xi }ni=1 . We shall make repeated use of this approximation-error bound in (30) in the sequel
for studying sup-norm posterior convergence rate and frequentist converge of Bayesian credible
intervals/bands.

3.1

Point-wise and sup norm posterior convergence rate

We first state a result for the posterior rate of contraction around f ∗ in terms of the pointwise and supremum norm. The Hellinger and total variation norms are by far the most common
metrics for establishing posterior contraction rates and there is only a recent small literature
for stronger norms [9, 22, 23, 53]. The following theorem establishes that the GP posterior
concentrates around the true function in both the point-wise and supremum norms at the
respective minimax rates for the Sobolev and Hölder classes.
3.1. Assume the kernel K satisfies assumptions (B) and (E). Define h via the equa

−1/(2α+1)
−1/(2α)
tion λ =
If h  B 2 n/(σ 2 log n)
for f ∗ ∈ ΘαH (B) and h  B 2 n/(σ 2 log n)
for f ∗ ∈ ΘαS (B), then, with probability at least 1 − n−10 with respect to the randomness in
{Xi , wi }ni=1 , the following error estimates for the squared point-wise risk hold for any x∗ in X ,

 2α


 2/(2α+1) σ2 log n (2α+1)

if f ∗ ∈ ΘαH (B),
B
n
∗
∗ ∗ 2
(31)
E[|f (x ) − f (x )| | Dn ] 
 2α−1

2α
2 log n

σ

B 1/α
if f ∗ ∈ ΘαS (B).
n
Theorem

h2α .

With the same probability and the same choice of h, the following error estimates for the squared
supremum norm hold,

 2α

(2α+1)


σ 2 log n
2/(2α+1)

if f ∗ ∈ ΘαH (B),
B
n
∗ 2
(32)
E[kf − f k∞ | Dn ] 
 2α−1

2α

σ 2 log n
1/α
∗
α

B
if f ∈ ΘS (B).
n
Two remarks are in order. First, we can also apply similar techniques to show error estimates
regarding the derivatives f (k) of f with order k up to bαc by identifying the posterior covariance
function of f (k) (whose posterior is also a GP) with certain noiseless KRR estimate, which also
applies to results in the following subsections. Due to the space constraint, we omit these
results. Second, since the sup-norm is stronger than the L2 norm, our result shows that by
introducing the scaling (nλ)−1 in the covariance kernel of the prior GP, we can eliminate the
mismatch [44, 46] between the smoothness level (which is α + 1/2 for the univariate case and
α + d/2 for the d-variate case) of the RKHS of the α-Matérn kernel and the best smoothness
level of the truth (which is α without this scaling) that the posterior can adapt to.

3.2

Frequentist coverage of point-wise posterior credible intervals

In the following, we leverage on the (higher order) expansions for both the posterior mean
(27) and the variance (30) to derive the frequentist coverage properties of point-wise Bayesian
credible intervals.
Let Φ denote the N(0, 1) c.d.f. and for any γ ∈ (0, 1), let zγ denote the γth standard normal
enB ), we consider for
quantile with Φ(zγ ) = γ. Since the posterior distribution of f is GP(fbn , C
any x ∈ X a point-wise credible interval centered at fbn (x) with level β ∈ (0, 1) as


CIn (x; β) = fbn (x) − ln (x; β), fbn (x) + ln (x; β) ,

13

where the half length
q
ln (x; β) = z(1+β)/2

enB (x, x) ,
C

is chosen so that the posterior probability of f (x) falling into the credible interval is β, or


P f (x) ∈ CIn (x; β) Dn = β.
P
−1 implies
e
Note that a combination of (30) and the fact that supx,x0 |K(x,
x0 )| . ∞
j=1 νj . h
that the size of the credible interval CIn (x; β) is O((nh)−1/2 ), which is minimax-optimal under
choices of h in Theorem 3.1. Our goal below is to investigate the frequentist coverage of
CIn (x; β), and in particular, identify situations when


Pρ f ∗ (x) ∈ CIn (x; β) ≥ β.
Note that CIn (x; β) is a random interval in the above display and Pρ denotes the probability
under the true data generating distribution that CIn (x; β) contains the true function value
f ∗ (x). Let γ = (1 + β)/2. Write


Pρ f ∗ (x) ∈ CIn (x; β)


q
q
√ 

∗
B
B
b
e
e
= Pρ − zγ nh Cn (x, x) ≤ nh fn (x) − f (x) ≤ zγ nh Cn (x, x)


q
q
√ 

(i)
∗
B
B
b
b
b
≈ Pρ − zγ Cn (x, x) ≤ nh fn (x) − f (x) ≤ zγ Cn (x, x)

q
i
√ 
 h√
∗
∗
B
b
b
= Pρ nh fn (x) − Fλ f (x) ∈
(33)
nh Pλ f (x) ± zγ Cn (x, x) .
The approximation (i) in the above display follows from (30) andpthe approximation bound
P
e X (·) from
will be made concrete inside the proof. Recall the process U (·) = h/n ni=1 wi K
i
c
b
the discussion after Corollary 2.1. Consider a GP W ∼ GP(0, Cn ) whose covariance function
bn matches the covariance function of the process U (·) over X , that is, for any (x, x0 ) ∈ X 2 ,
C
∞
X


ψj (x) ψj (x0 )
0
0
2
0
2
b
e
e
Cn (x, x ) = E[U (x)U (x )] = σ h E K(X, x) K(X, x ) = σ h
.
(1 + h2α /µj )2
j=1

√
The second display in (27) can be used to establish that the law of the random process nh (fbn −
c . Substituting in (33) and tracking the approximation
Fλ f ∗ ) is (point-wise) close to that of W
errors leads to the following theorem.
Theorem 3.2 (Frequentist coverage of posterior credible intervals). There exists some constant
C independent of (n, h) such that the frequentist coverage of CIn (x; β) satisfies that for any
x ∈ X,
 1


 


P f ∗ (x) ∈ CIn (x; β) − Φ un (x; β) + bn (x) − Φ − un (x; β) + bn (x) ≤ C √
+ γn + δ n ,
nh
q
bnB (x, x)/C
bn (x, x) z(1+β)/2 is the inflated quantile and
where un (x; β) = C
√

bn (x, x) −1/2 nh Pλ f ∗ (x) is the bias at x.
bn (x) = C

We briefly comment on the source of each of the three approximation error terms on the
right hand side of the previous display; details can be found inside the proof. The distributional

bn (x, x) −1/2 U (x) with a N(0, 1) random variable incurs an error of (nh)−1/2
approximation of C
14

from the Berry–Essen theorem. The approximation of the appropriately standardized posterior
√


bn (x, x) −1/2 U (x) results in an error of δn
bn (x, x) −1/2 nh(fˆn (x) − Fλ f ∗ (x)) with C
mean C
from (27). The approximation (30) contributes an error of γn .
Come back to the concrete case when [0, 1] is the unit interval and K is the Matérn kernel
with smoothness index α > 1/2, where ψ2j (x) = cos(πjx), ψ2j−1 = sin(πjx), and µ2j = µ2j+1 
j −2α for j = 0, 1, . . .. Applying the identity cos(x − y) = cos(x) cos(y) + sin(x) sin(y), we can
bnB and C
bn into
simplify the two covariance functions C
bnB (x, x0 ) = σ 2 h
C

∞
X
cos(πj(x − x0 ))
j=1

1 + h2α /µ2j

bn (x, x0 ) = σ 2 h
and C

∞
X
cos(πj(x − x0 ))
j=1

(1 + h2α /µ2j )2

.

(34)

Define
CIR

b B (x, x)/C
bn (x, x) =
:= C
n

∞
X
j=1

∞
X

1
1+

h2α /µ2j

j=1

1
(1 +

h2α /µ2j )2

,

where we have applied identities in (34) so that the CIR is independent of x. Clearly CIR > 1.
For simplicity, we consider f ∗ ∈ ΘαH0 (B). All the results naturally extends to f ∗ ∈ ΘαS 0 (B).
3.1 (Matérn kernel class). (Under-smooth) Suppose f ∗ ∈ ΘαH0 (B) for α0 > α,
then for any x ∈ X , bn (x) → 0 as n → ∞, and


P f ∗ (x) ∈ CIn (x; β) → 2 Φ(CIR z(1+β)/2 ) − 1 ∈ (β, 1), as n → ∞.

Corollary

(Smooth-match) For any βe ∈ (0, 1), there is some sufficiently large constant B > 0 and a
sequence of functions {fn∗ }, each belongs to ΘαH0 (B), such that


e
P fn∗ (x) ∈ CIn (x; β) → β,

as n → ∞.

∗ ∈ Θα0 (B), such that for
(Over-smooth) For any α0 < α, there always exists a function fbad
S
any x ∈ X , bn (x) → ∞ as n → ∞, and

 ∗
(x) ∈ CIn (x; β) → 0, as n → ∞.
P fbad

3.3

Frequentist coverage of simultaneous posterior credible bands

In this subsection, we study the frequentist coverage of the following posterior credible band
centered at the posterior mean fbn with level β ∈ (0, 1),
n
o
CBn (β) = f ∈ L2 (X ) : f − fbn ∞ ≤ rn (β) ,
where the half length rn (β) is chosen so that posterior probability of f falling into the credible
band is β, i.e.,


P f ∈ CBn (β) Dn = β.
Whereas the point-wise intervals CIn (x; β) permitted an explicit description using the Gaussianity of f (x) | Dn for any x, the band CBn (β) needs to be defined implicitly due to the lack of
a similar exact distributional result. Establishing the frequentist validity of the band CBn (β)
follows a conceptually similar √
route as before; (i) approximate the sampling distribution of the
c ∼ GP(0, C
bn ) define previously,
standardized posterior mean nh (fbn − Fλ f ∗ ) by the GP W
√
and (ii) approximate the centered and scaled posterior measure nh (f − fbn ) | Dn by the GP
15

c B ∼ GP(0, C
bnB ). However, substantial work is necessary to obtain uniform counterparts of
W
the point-wise approximations obtained previously.
√
We first discuss approximation of the sampling distribution of nh (fbn −Fλ f ∗ ) in supremum
c is defined as a GP with law GP(0, C
bn ), where C
bn is the covariance
norm. Recall that W
p
Pn
e X (·) over X , which is the leading term
function of the random process U (·) = h/n i=1 wi K
i
√

∗
bnM = kW
c k∞ and
in the expansion of n h fbn − Fλ f in the second display in (27). Let Z
M
Zn = U ∞ . We will first show by applying the results in [15] on Gaussian approximation to
bnM and ZnM are close with respect
the suprema of empirical processes that the distributions of Z
to the Kolmogorov distance.
√
We also define ZenM = n h fbn −Fλ f ∗ ∞ as the rescaled sup-norm deviation of the posterior
mean function from its population counterpart Fλ f ∗ = (I − Pλ )f ∗ . By applying an anticoncentration bound [14] for GP, we obtain by combining the second display in (27) that the
eM can be well-approximated by that of Z M . Finally, a combining of these two
distribution of Z
n
n
enM can be well-approximated by ZbnM ,
approximation results implies that the distribution of Z
c , as summarized in the following theorem.
the supremum of a Gaussian process W
3.3 (Gaussian process approximation of the posterior mean function). There exists
some constant C independent of (n, h) such that for any t ≥ 0,
√




log n
.
(35)
P ZbnM ≤ t − P ZenM ≤ t ≤ C
(nh)1/8
Theorem

Next, we show that the posterior√measure can also be uniformly approximated by another
enB ), where feB has the law of (f −
specially designed GP. Recall that nh feB ∼ GP(0, nh C
enB depends on the random design {Xi }n and is hard
fbn ) | Dn . The covariance function C
i=1
bB )
to directly work with. For this reason, we defined a population level GP W B ∼ GP(0, C
n
B
B
b
e
whose √
covariance function Cn provides good approximation (30) to nh Cn in sup-norm. Let
eB = n h kfeB k∞ and Z B = kW B k∞ . In the next theorem, we show that the distributions of
Z
n
n
eB and its population level counterpart Z B are close with respect to the Kolmogorov distance.
Z
n
n
To prove Theorem 3.4, we develop a new Gaussian comparison inequality; see Theorem 5.1 in
§5.8; that explicitly bounds the Kolmogorov distance between two GPs in terms of the supremum
norm difference between their covariance functions. The comparison inequality extends the
Sudakov–Fernique inequality for finite-dimensional Gaussians as stated in [12] to GPs.
3.4 (Gaussian process approximation of the centered posterior measure). There
exists some constant C independent of (n, h) such that for any t ≥ 0,
Theorem

 B



e ≤ t Dn − P Z B ≤ t ≤ C γ 1/3 log n.
P Z
n
n
n

(36)

c ∼ GP(0, C
bn ) and W B ∼
From Theorem 3.3 and Theorem 3.4, we only need to compare W
B
b ) in order to study the frequentist coverage of the Bayesian credible band CBn (β).
GP(0, C
n
B
Let qn (β) denote the β-th quantile of the population level random variable ZnB = kW B k∞ .
3.5 (Frequentist coverage of posterior credible bands). There exists some constant
C independent of (n, h), such that for any β ∈ (0, 1),

Theorem

1
C
rn (β) − √ qnB (β) ≤ √ γn1/3 log n.
nh
nh

16

(37)

√
nh kPλ f ∗ k∞ → 0 as n → ∞, then
√



log n √
∗
+
nh
kP
f
k
≤ qnB (β) ≤ C γn1/3 log n +
∞ → 0,
λ
(nh)1/8
(38)

Moreover, if the bias term kPλ f ∗ k∞ satisfies



c k∞
P f ∗ ∈ CBn (β) − P kW



c k∞ ≤ q B (β) ∈ (β, 1). On the other hand, if the bias term
as n → ∞ and h → 0, where P kW
n
√
kPλ f ∗ k∞ satisfies nh kPλ f ∗ k∞ → ∞ as n → ∞, then


P f ∗ ∈ CBn (β) → 0, as n → ∞ and h → 0.
(39)
Since (nh)−1/2 qnB (β) = O((nh)−1/2 ), the width of the simultaneous credible band
√ is again
minimax-optimal under choices of h in Theorem 3.1. Unlike the point-wise case where nh (fbn (x)−
f ∗ (x)) weakly converges to a nondegenerate distribution, in the simultaneous case, the two nc generally do not weakly converge to non-degenerate
dependent population level GPs W B and W
laws (tight GPs) over [0, 1] as n → ∞ and h → 0, due to the kernel type estimator form of the
process U with an (n, h) dependent kernel (see, for example, [30, 42]).
Similar to Corollary 3.1, we have the following corollary when specialized to the Matérn
kernel with smoothness index α > 1/2 over [0, 1], and omit the proof.
3.2 (Matérn
kernel class). (Under-smooth) Suppose f ∗ ∈ ΘαH0 (B) for α0 > α,
√
then the bias term nh kPλ f ∗ k∞ → 0 as n → ∞, and


P f ∗ ∈ CBn (β) ∈ (β, 1), as n → ∞.

Corollary

(Smooth-match) For any βe ∈ (0, 1), there is some sufficiently large constant B > 0 and a
sequence of functions {fn∗ }, each belongs to ΘαH0 (B), such that


e
P f ∗ ∈ CBn (β) → β,

as n → ∞.

∗ ∈ Θα0 such that the bias
(Over-smooth)
For any α0 < α, there always exists a function fbad
H
√
term nh kPλ f ∗ k∞ → ∞ as n → ∞, and


P f ∗ ∈ CBn (β) → 0, as n → ∞.

4

Simulation study

In the following, we numerically investigate the behavior of the point-wise and simultaneous
credible intervals for a certain f ∗ when a Gaussian process prior with Matérn covariance kernel
given by
K(x, y) =

21−ν √
2ν |x − y|
Γ(ν)

ν

Bν

√

2ν |x − y| ,

x, y ∈ [0, 1],

(40)

where Bν is the modified Bessel function of the second kind for 0 < ν < ∞. We recall that
K with ν = (α − 1/2) satisfies (B) and (E) when expanded with respect to the Fourier
basis ψ2j−1 (x) = sin(πjx), ψ2j (x) = cos(πjx), j = 1, 2, . . .. Also the eigen-values satisfy
λ2j−1 = λ2j  j −2α .P
−1.7 sin j} cos(π(j − 0.5)x). Training data of size n are drawn from
We let f ∗ (x) = ∞
j=1 {j
∗
the model yi = f (xi ) + i with i ∼ N(0, 0.12 ) and xi ∼ Unif(0, 1). 200 equally spaced points
on [0, 1] are chosen as the test x values. For Figure 1, the prior has smoothness 0.6 where
f ∗ ∈ ΘαH0 (B) for α0 < 0.7. For Figure 2, the prior has smoothness 1.7 and f ∗ satisfies the
17

counter example in the proof of Corollary 3.1 in the over-smooth case. The posterior mean and
covariance functions are directly obtained using (5) and (6), from where the point-wise 95%
credible intervals at the 200 test points are constructed. To obtain the simultaneous credible
bands, 1000 random samples are drawn from the posterior distribution of the function evaluated
at the 200 test points. Then rn (0.95) is estimated from the 95% quantile of the posterior samples
of kf − fˆn,λ k∞ . In the over-smooth and the under-smooth case, it is suggestive from Figures 1
and 2 that the simultaneous and point-wise coverage tends to a non-zero fraction and zero in
the respective cases as the sample size increases to ∞, which is consistent with the prediction
from our theory.
The coverage of the simultaneous β%-credible intervals is further investigated through a
replicated simulation study. We consider β = 0.80 and 0.90. Using the same f ∗ , we consider
three smoothness parameters (0.1, 0.15, 1.2) of the Matérn covariance kernel. 1000 datasets were
drawn from the same model and the proportion of datasets for which f ∗ lies in the simultaneous
interval is recorded in Table 1. In the over-smooth case, the coverage tends to 0 and in the
under-smooth case, the coverage tends to a number between β and 1. It may appear from the
simulation study that in the under-smooth case, the coverage probability is tending to 1 as
the sample size increases. This is primarily because ν is close to zero and the actual coverage
probability becomes very close to 1 in these cases. As a result when the sample size is very large,
one needs a huge number of replicates to actually sample from the tail event that a dataset for
which the simultaneous credible band does not cover the true function.

n
Prior
Matérn(ν = 0.1)
Matérn(ν = 0.15)
Matérn(ν = 1.2)

200
0.80
0.977
0.875
0

500

0.90
0.993
0.939
0

0.80
0.995
0.926
0

0.90
0.999
0.953
0

2000
0.80
0.998
0.978
0

0.90
0.999
0.993
0

Table 1: Simultaneous coverage probability over 1000 replicates.

18

(a) n = 200

(b) n = 500

(c) n = 2000

Figure 1: Credible intervals with Matérn with ν = 0.1. The red and the blue lines represent f ∗
and the point-wise posterior mean respectively. The upper and lower brackets of the dark-grey
and the light-grey shaded regions represent the point-wise and the simultaneous 95% credible
intervals.

19

(a) n = 200

(b) n = 500

(c) n = 2000

Figure 2: Credible intervals with Matérn with ν = 1.2. The red and the blue lines represent f ∗
and the point-wise posterior mean respectively. The upper and lower brackets of the dark-grey
and the light-grey shaded regions represent the point-wise and the simultaneous 95% credible
intervals.

20

5

Proofs

In this section, we provide proofs of the results in the paper. For notational brevity, we drop
λ from the subscript of k · kλ within the proofs; all instances of k · k refer to the RKHS norm
(for the equivalent kernel) k · kλ .

5.1

Summary of notations

For the reader’s convenience, Table 2 provides a summary of some of the key notations
introduced earlier that are heavily used inside the proofs.
Table 2: Notation reference
Symbol
Dn
f∗
wi
fbn
enB
C
fbn,λ
e
K
{νj }∞
j=1
Fλ
Pλ
Sn,λ
Sλ

j=1

eigenvalues of equivalent kernel, νj = µj /(µj + λ)
R
e t)ds
convolution with equivalent kernel, Fλ g(t) = g(s)K(s,
Pλ = id − Fλ
P
e X − Pλ f
sample score function, Sn,λ f = n−1 n (Yi − f (Xi ))K
i=1

i

population score function, Sλ f = Fλ f ∗ − f
p
P
e X (·)
second-order scaled bias approximation, U (·) = h/n ni=1 wi K
i
bn (x, x0 ) = EU (x)U (x0 )
covariance function of U (·), C

c
W
eM
Z

c ∼ GP(0, C
bn )
Gaussian approximation to U , W
√
eM = nhkfbn − Fλ f ∗ k∞
rescaled sup-norm bias of posterior mean, Z

ZnM
bnM
Z

sup-norm of U , ZnM = kU k∞
c, Z
bnM = kW
c k∞
sup-norm of W

feB
bB
C

e B ), law of centered GP posterior
feB ∼ GP(0, C
n
eB
deterministic approximation to nh C
n
√
B
B
b
W ∼ GP(0, Cn ), approximation to law of nhfeB
√
enB = nh kfeB k∞ , sup-norm of centered and scaled posterior GP
Z
ZnB = kW B k∞ , sup-norm of W B

n
WB

enB
Z
ZnB

5.2.1

posterior mean function, E[f | Dn ]
enB (x, x0 ) = Cov(f (x), f (x0 ) | Dn )
posterior covariance function, C
 −1 Pn
2 
2
KRR solution, argmin f ∈H n
i=1 (Yi − f (Xi )) + λ kf kH
e t) = P∞ νj ψj (s)ψj (t)
equivalent kernel, K(s,

U
bn
C

n

5.2

definition
data, {(Yi , Xi ), i = 1, . . . , n}
true function
Gaussian error, wi = Yi − f ∗ (Xi ) ∼ N(0, σ 2 )

n

Proof of Theorem 2.1
A First order bound for fbn,λ

Our analysis proceeds in two steps: we first prove a rough estimate for the error bound
relative to the RKHS norm k · k using a sub-Gaussian type inequality, based on which we
21

can control the covering entropy and further obtain a sharper sup-norm error bound using a
Bernstein type inequality. We note that the sup-norm bound cannot be directly obtained by
applying the Sobolev inequality that relates k · k with the sup norm.
A rough error bound relative to the k · k norm: We make use of the following subGaussian type concentration inequality under the k · k norm.
5.1 (Lemma 6.1 in [52]). Suppose α > 1/2. There exists some constant CK > 0
depending only on the kernel K, such that for any t > 0,
Lemma




P sup
g∈H



n

 
1 X
2
−β
e
e
g(Xi ) KXi − E g(X) KX
kgk ≥ CK n t ≤ e−t ,
n

where the constant β =

i=1

(2α − 1)2
is strictly positive.
4α (2α + 1)

Let An denote an event such that under An , it holds for all g ∈ H,
n


1 X
e X ≤ CK n−β/2 kgk.
e X − E g(X) K
g(Xi ) K
i
n

(41)

i=1

According to Lemma 5.1, we may choose An such that P(An ) ≥ 1 − n−10 for any sufficiently
large n.
Notice that we have the following two identities regarding the finite-sample level and the
population level score functions,
Sn,λ (fbn,λ ) = 0

and

Sλ (Fλ f ∗ ) = 0.

(42)

Let ∆f = fbn,λ − Fλ f ∗ denote the difference between the finite-sample minimizer and the population level minimizer of the KRR. From the identity (15), note also that ∆f = −Sn,λ (Fλ f ∗ ).
Now we will obtain a bound on ∆f by repeatedly applying inequality (41).
By the definition of operators Sn,λ and Sλ , we have


n


 
 1 X
∗
∗
b
b
e X − E ∆f (X) K
eX .
∆f (Xi ) K
Sn,λ (fn,λ ) − Sλ (fn,λ ) − Sn,λ (Fλ f ) − Sλ (Fλ f ) =
i
n
i=1

(43)
Therefore, we obtain by applying inequality (41) that under An ,


 

Sn,λ (fbn,λ ) − Sλ (fbn,λ ) − Sn,λ (Fλ f ∗ ) − Sλ (Fλ f ∗ ) ≤ CK n−β/2 k∆f k.

On the other hand, by using the identities in (42) along with the identity ∆f = −Sn,λ (Fλ f ∗ ),
we get

 

Sn,λ (fbn,λ ) − Sλ (fbn,λ ) − Sn,λ (Fλ f ∗ ) − Sλ (Fλ f ∗ ) = ∆f − Sn,λ (Fλ f ∗ ).
(44)

22

Let us know attempt to bound kSn,λ (Fλ f ∗ )k. By applying the triangle inequality, we get
kSn,λ (Fλ f ∗ )k = Sn,λ (Fλ f ∗ ) − Sλ (Fλ f ∗ )
n


1X ∗
e X − E {f ∗ − Fλ f ∗ }(X) K
eX
≤
{f − Fλ f ∗ }(Xi ) K
i
n
i=1

n

1X
eX
{Yi − f ∗ (Xi )}K
i
n

+
1
=
n

(i)

i=1
n
X
i=1

(45)

n


1X
∗
eX
e
e
Pλ f (Xi ) KXi − E Pλ f (X) KX +
wi K
i
n
∗

i=1

:= T1 + T2 ,
where in step (i) we used the decomposition id = Pλ + Fλ , and recall that wi = Yi − f ∗ (Xi ) is
the i.i.d. N (0, σ 2 ) noise. Applying inequality (41), the first term T1 can be bounded as
T1 ≤ CK n−β/2 kPλ f ∗ k.
e i , Xj )]1≤i,j≤n and w = (w1 , . . . , wn )T . Then we
To bound the second term T2 , let Σ = [K(X
have T22 = n2 wT Σw. By the Hanson-Wright inequality [34], we have
h
p
i
2
P wT Σw ≥ σ 2 tr(Σ) + 2 tr(Σ2 ) t2 + 2kΣkF t2 ≤ e−t , ∀t > 0,
where k · kF denotes the matrix Frobenius norm. For some constant DK only depending on K,
tr(Σ) =
tr(Σ2 ) =

n
X

e i , Xi ) ≤ DK n h−1 ,
K(X

i=1
n
X

2
e i , Xj )2 ≤ DK
K(X
n2 h−2 ,

i,j=1

kΣkF =
we obtain
P

p
tr(Σ2 ) ≤ DK n h−1 ,

!
n
1X
e X ≥ DK σ (nh)−1/2 (1 + 2t + 2t2 ) ≤ e−t2 .
wi K
i
n
i=1

Therefore, there exists some event Bn with P(Bn ) ≥ 1 − n−10 such that on Bn we have
T2 ≤ Dk σ (nh)−1/2 log n,
where we slightly abuse the notation by using DK to mean a different constant depending only
on K.
By putting pieces together, we obtain that on An ∩ Bn ,
k∆f k ≤ k∆f − Sn,λ (Fλ f ∗ )k + kSn,λ (Fλ f ∗ )k
≤ CK n−β/2 k∆f k + CK n−β/2 kPλ f ∗ k + Dk σ (nh)−1/2 log n.
For n sufficiently large such that CK n−β/2 ≤ 1/2, the above implies
k∆f k ≤ 2 CK n−β/2 kPλ f ∗ k + 2 Dk σ (nh)−1/2 log n.

23

(46)

This further implies
kfbn,λ − f ∗ k ≤ kfbn,λ − Fλ f ∗ k + kf ∗ − Fλ f ∗ k ≤ 2 kPλ f ∗ k + 2 Dk σ (nh)−1/2 log n.

(47)

A sharp error bound relative to the sup norm: From inequality (46) and the definition
of k · k, we obtain the following control on the Sobolev norm k · kH


en .
k∆f kH = kfbn,λ − Fλ f ∗ kH ≤ 2λ−1/2 kPλ f ∗ k + Dk σ (nh)−1/2 log n : = A
Moreover, we also have k∆f k∞ = kfbn,λ − Fλ f ∗ k∞ ≤ C h−1/2 kfbn,λ − Fλ f ∗ k ≤ C 0 n for some
sufficiently large constant C 0 > 0.
For any An > 0 and Bn > 0, let Gn = {g ∈ H : kgk∞ ≤ Bn , kgkH ≤ An }. We will make
use of the following Bernstein type concentration inequality for suprema of empirical processes
under the sup norm. A proof of Lemma 5.2 can be found in §5.12.
Lemma


P

5.2. For any α > 1/2, it holds for some constant C independent of (n, h, An ) that
n

1X
e i , t) − E[g(X1 )K(X
e 1 , t)]
g(Xi )K(X
t∈T ;g∈Gn n i=1
s
hp
i
√
Bn
log n + x + log
+ A1/(2α)
kgk∞
> (nh)−1/2
n
kgk∞

h
i
Bn
1/α
1/(2α)−1/2
−1
kgk∞ ≤ 2e−x ,
log n + x + log
+ An max 1, (n/h)
+ (nh)
kgk∞
sup

x > 0.

en , Bn = n and t = c log n. Therefore,
Let Bn denote the event in this lemma with An = A
−10
for sufficiently large c we have P[Bn ] ≥ 1 − n , and under event Bn we have that for all g ∈ H
en ,
satisfying kgk∞ ≤ n and kgkH ≤ A
n

1X
e X − E[g(X)K
eX ]
g(Xi )K
i
n
∞
i=1
r
i kgk
h
 1/α
n
n
en log
√ ∞
≤C
log
+ (nh)−1/2 max 1, (n/h)1/(2α)−1/2 A
min(1, kgk∞ )
min(1, kgk∞ )
nh

r
r
n
n
kgk∞
√ ,
≤ C 0 max 1, n−1+1/(2α) h−1/(2α) log
log
min(1, kgk∞ )
min(1, kgk∞ )
nh
(48)
for any α > 1/2 and λ = h2α  n−2α/(2α+1) .
Let us divide the proof into two cases:
If k∆f k∞ ≤ kPλ f ∗ k∞ :

The claimed bound follows since kfbn,λ − f ∗ k∞ ≤ k∆f k∞ + kPλ f ∗ k∞ .

n
If k∆f k∞ ≥ kPλ f ∗ k∞ : In this case, we have log min(1, k∆f
k∞ )  log n. Similar to the previous
analysis, by applying (48) with g = ∆f , and comparing identities (43) and (44), we can get

∆f − Sn,λ (Fλ f ∗ )

∞

≤ C 0 max 1, n−1+1/(2α) h−1/(2α)

p
p
k∆f k∞
log n
log n √
.
nh

(49)

In order to bound Sn,λ (Fλ f ∗ ) ∞ , we will make use of the following concentration inequality
for bounding T2 relative to the sup norm.

24

Lemma

5.3. There exists some constant C independent of (n, h) such that for any x > 0,
r


n
1+x
1X
e
P
≤Cσ
≤ e−x .
(50)
wi KXi
n
nh
∞
i=1

By conducting a similar analysis as the steps in (45) with k · k being replaced by k · k∞ , and
using Lemma 5.3 and inequality (48), we can obtain that under intersection of event Bn and
the event in Lemma 5.3 with x = c log n,
n
o r log n
p

−1+1/(2α) −1/(2α)
∗
∗
h
Sn,λ (Fλ f ) ∞ ≤ C max 1, n
log n kPλ f k∞ + σ
.
(51)
nh
Now combining inequalities (49) and (51), we get
∆f

∞

≤C

0

∗

Sn,λ (Fλ f )

∞

≤C

n

max 1, n

−1+1/(2α)

h

−1/(2α)

o
p

log n kPλ f ∗ k∞ + σ

r

log n
,
nh

which implies the claimed bound as
r
∗

∗

∗

kfbn,λ − f k∞ ≤ k∆f k∞ + kPλ f k∞ ≤ 2kPλ f k∞ + C σ
5.2.2

log n
.
nh

(52)

A second order bound for fbn,λ

Define gn := fbn,λ − f ∗ . By the definition of operators Sn,λ and Sλ , we have
n

 
 1 X


eX .
e X − E gn (X) K
Sn,λ (fbn,λ ) − Sλ (fbn,λ ) − Sn,λ (f ∗ ) − Sλ (f ∗ ) =
gn (Xi ) K
i
n

(53)

i=1

On the other hand, equation (42) leads to the identity

 

Sn,λ (fbn,λ ) − Sλ (fbn,λ ) − Sn,λ (f ∗ ) − Sλ (f ∗ ) = −Sλ (fbn,λ ) + Sλ (f ∗ ) − Sn,λ (f ∗ )
(i)

= gn − Sn,λ (f ∗ ),

where in step (i) we used (15).
Let us divide the proof into two cases:
If kgn k∞ ≤ h1/2 kPλ f ∗ k∞ : The claimed bound follows by using the fact that


eX
E gn (X) K
≤ C h−1/2 kgn k∞ .

1
n

Pn

i=1 gn (Xi ) KXi −

e

∞

If k∆f k∞ ≥ h1/2 kPλ f ∗ k∞ : In this case, we have log min(1, nkgn k )  log n. Therefore, inequal∞
ity (48) and identity (53) together imply that on Bn ,

 

Sn,λ (fbn,λ ) − Sλ (fbn,λ ) − Sn,λ (f ∗ ) − Sλ (f ∗ )
0

≤ C max 1, n

−1+1/(2α)

−1/(2α)

h

∞

p
p
kgn k∞
log n
log n √
.
nh

Combining the last two displays and the first order bound (52) for kgn k∞ , we obtain that

25

for sufficiently large n, on event An ∩ Bn ,
fbn,λ − f ∗ −




n
1 X
∗
e
wi KXi − Pλ f
n
i=1

∞

(54)
r

p

≤ C 0 max 1, n−1+1/(2α) h−1/(2α) log n

log n 
2kPλ f ∗ k∞ + C σ
nh

r

log n 
.
nh

Moreover, since the randomness in An ∩ Bn is with respect to the noise {wi : i = 1, . . . , n}
and the random design {Xi : i = 1, . . . , n}, this inequality also holds uniformly for all f ∗ as
well, meaning that it holds for all f ∗ where fbn,λ is defined as a random function depending on
{Xi : i = 1, . . . , n} and {Yi = f ∗ (Xi ) + wi : i = 1, . . . , n}.

5.3

Proof of Corollary 2.1

We prove the bounds on kPλ f ∗ k∞ stated before the Corollary; the proof then is immediate
upon optimally choosing h.
e x k. First, kK
e x k2 =
First suppose f ∗ ∈ ΘαS (B). Bound kPλ f ∗ k∞ ≤ kPλ f ∗ k supx∈X kK
P
−1 by (18). Hence, sup
−1/(4α) . Next,
e
e
K(x,
x) . ∞
x∈X kKx k . λ
j=1 νj . h
kPλ f ∗ k2 = λ hPλ f ∗ , f ∗ iH = λ

∞
X
j=1

∞
X
(fj∗ )2
λ
≤λ
j 2α (fj∗ )2 .
µj + λ µj
j=1

This leads to kPλ f ∗ k∞ ≤ λ1/2−1/(4α) = h(α−1/2) by using the definition of ΘαS (B).
Now suppose f ∗ ∈ ΘαH (B). Recall
∗

Pλ f =

∞
X

(1 −

νj )fj∗ ψj

j=1

=

∞
X
j=1

λ
f ∗ ψj .
µj + λ j

Bound
∗

kPλ f k∞ ≤

∞
X
j=1

∞ p
∞
√ X
√ X
√
λµj fj∗
λ
∗
fj = λ
j α fj∗ . λ = hα ,
√ . λ
µj + λ
µj + λ µj
j=1

j=1

where we used the definition of ΘαH (B) and the AM-GM inequality to bound
1/2.

5.4

p
λµj /(µj + λ) ≤

Proof of Theorem 3.1

We prove the the first part. From the discussion in §3, we obtain a high probability bound
for the posterior expected point-wise loss. With large probability
E[|f (x∗ ) − f ∗ (x∗ )|2 | Dn ] = |fˆn,λ (x∗ ) − f ∗ (x∗ )|2 + Var[f (x∗ ) | Dn ]
log n
+ kPλ f ∗ k2∞ + kPλ Kx∗ k2∞ .
- σ2
nh
Since kPλ Kx∗ k2∞ - σ 2 /(nh), choosing h optimally, we have (31).
To prove (32), observe that
E[kf − f ∗ k2∞ | Dn ] ≤ 2kfˆn,λ − f ∗ k2∞ + 2Ekf − fˆn,λ k2∞ .

26

enB ), we use a version of Dudley’s entropy integral (TheoSince f − fˆn,λ | Dn := Z ∼ GP(0, C
rem
p 3.2 of [31]) for bounding moments of supremum of Gaussian processes. Define ρ(s, t) :=
Var{Z(s) − Z(t)} as the intrinsic psuedo-metric. From the proof of Theorem 4.4., we have
ρ(s, t) ≤ h−1 |s − t|1/2 . Hence the entropy of [0, 1]d with respect to ρ is of the order log n. The
d
d
width
√ of the process {Z(t) : t ∈ [0, 1] } with respect to t0 ∈ [0, 1] is given by supt ρ(t, t0 ) σ/ nh. Thus from Theorem 3.2 of [31]
σ2
σ 2 log n
Ekf − fˆn,λ k2∞ ≤
+
.
nh
nh
Again choosing h optimally, we have (32).

5.5

Proof of Theorem 3.2

√
Let us first find a normal approximation to the distribution of nh (fbn (x) − Fλ f ∗ (x)) for
any fixed x ∈ X . Using the Berry Esseen theorem, we obtain the following Kolmogorov disp

bn (x, x) −1/2 h/n Pn wi K
e X (x) and a
tance bound between the sampling distribution of C
i
i=1
standard normal random variable
h
bn (x, x)
sup P C

−1/2 p

u∈R

h/n

n
X

i
e X (x) ≤ u − Φ(u) ≤ √C ,
wi K
i
nh
i=1

∀x ∈ X ,

where Φ denotes the cdf of a standard normal random variable, and we have used the fact that
e X (x) is bounded by h−1/2 up to a constant. Since function Φ has
the third moment of h1/2 wi K
i
a bounded derivative, we can combine the above with the second display in (27) to obtain
h
bn (x, x)
sup P C

−1/2 √

u∈R

i

 1
nh (fbn (x) − Fλ f ∗ (x)) ≤ u − Φ(u) ≤ C √
+ δn ,
nh

∀x ∈ X .

Now let us turn to the posterior credible interval CIn (x; β). By combining the fact that

e B (x, x) and inequality (30), we obtain
f (x) ∼ N fbn (x), C
n
sup P
u∈R

h
bnB (x, x)
C

−1/2 √

i

nh f (x) − fbn (x) ≤ u Dn − Φ(u) ≤ C γn ,

∀x ∈ X .

This implies that the half length of the credible interval CIn (x; β) satisfies
s
s
B
b
bnB (x, x)
Cn (x, x)
C
ln (x; β) −
z(1+β)/2 ≤ C
γn ,
nh
nh
where we used the fact the inverse function Φ−1 of Φ has a bounded derivative, and zu denotes
the uth quantile of a standard normal distribution for u ∈ (0, 1).
We can express the frequentist coverage as


P f ∗ (x) ∈ CIn (x; β)


= P − ln (x; β) ≤ fbn (x) − f ∗ (x) ≤ ln (x; β)
h 
√
√

bn (x, x) −1/2 nh ln (x; β) + C
bn (x, x) −1/2 nh Pλ f ∗ (x)
=P − C
√

bn (x, x) −1/2 nh (fbn (x) − Fλ f ∗ (x))
≤ C
i
√
√


bn (x, x) −1/2 nh ln (x; β) + C
bn (x, x) −1/2 nh Pλ f ∗ (x) .
≤ C

27

Combining this with the previous approximation bounds, we can obtain

s b B
√

 ∗

Cn (x, x)
bn (x, x) −1/2 nh Pλ f ∗ (x)
z(1+β)/2 + C
P f (x) ∈ CIn (x; β) − Φ
bn (x, x)
C
s



 1
bnB (x, x)

C
−1/2 √
∗
b
nh Pλ f (x) ≤ C √
+ γn + δn ,
+Φ −
z(1+β)/2 + Cn (x, x)
bn (x, x)
nh
C
which implies the desired result.

5.6

Proof of Corollary 3.1.

We define a sequence of integers {an } to be of constant separation if {an } is strictly increasing
and |an − an−1 | is constant for all n ≥ 2.
Under-smooth: For any x ∈ X , there exists a constant Cx depending on x, and a subsequence
kj , j = 1, . . . , ∞ with constant separation such that
bn (x, x) ≥ Cx h
C

∞
X
j=1

1
1+

h2α /µkj

= Cx .

Also
√

∗

1/2

nh |Pλ f (x)| ≤ λ

∞
X
j=1

j −α0 |fj∗ |
.
λ + j −2α j −α0

−1/(2α) .
For 2α > α0 , the function x 7→ x−α0 /(λ + x−2α ) is maximized at x = {α0 λ/(2α
√ − α0 )}∗
For 2α ≤ α0 , the function is monotonically decreasing.
In the first case nh |Pλ f (x)| ≤
√
∗
0.5(α
/α−1)
0
Cλ
which goes to 0. In the later case, nh |Pλ f (x)| ≤ Cλ1/2 which also goes to 0.

Smooth-match: Let (fn∗ )j = Cβ j −α0 for j = jn = bλ−1/(2α) c, and (fn∗ )j = 0 for j 6= jn , where
Cβ is some tuning parameter to be determined later. It is obvious that f ∗ ∈ ΘαH0 (B) for any
B ≥ Cβ . Moreover, we have
√

nhPλ (fn∗ )(x) = λ1/2

∞
X
j=1

j −α0 (fn∗ )
bλ−1/(2α) c−α
1/2
ψ
(x)
=
C
λ
ψj (x)
j
β
λ + j −2α j −α0
λ + bλ−1/(2α) c−2α n

: = cλ Cβ ψjn (x) ∈ [0.5 Cβ ψjn (x), Cβ ψjn (x)].
Let ebn be the solution of
e
Φ(un (x; β) + ebn ) − Φ((un (x; β) − ebn ) = β.

en (x, x)−1/2 cλ Cβ ψjn (x), then Theorem 3.2 implies P f ∗ (x) ∈
If we choose Cβ such that ebn = C
n

CIn (x; β) → βe as n → ∞.
∗
Over-smooth: Let fbad,j
= (µ∗j )1/2 aj where aj = j −(1+δ) for δ < α − α0 . Let β := 1 + δ + α0 .
There exists a subsequence kj with constant separation such that

√

∗
nhPλ fbad
(x) = λ1/2

∞
X
(µ∗j )1/2
j=1

λ + µj

aj ψj (x) & λ1/2

∞
X

1

j=1

λkjβ + kjβ−2α

1

λ

28

1

β

 λ 2 − 2α + 2α −1
δ+α0 −α
2α

Z

Z

dt
tβ + tβ−2α

dt
.
tβ + tβ−2α

(55)

Since integral to the right of (55) is finite and noting that δ + α0 − α < 0, the conclusion follows
immediately.

5.7

Proof of Theorem 3.3

By applying Theorem 2.1 in [15] about a Gaussian approximation to the suprema of empirical
processes (similar calculations as Corollary 2.2 therein), there exists a random variable Bn that
bM = kW
c k∞ , the suprema of the GP W
c ∼ GP(0, C
bn ), such that
has the same distribution as Z
n
for any γ ∈ (0, 1),
h
b1/3 v 2/3 log2/3 n i
≤ C2 γ,
P ZnM − Bn ≥ C1
γ 1/3 n1/6
e
e 2 (x)] . 1 and C1 , C2 are constants
where b = h1/2 supx,x0 |K(x,
x0 )| . h−1/2 , v 2 = supx E[h K
X
independent of n. Combining this with Lemma 2.3 in [15] that converts this coupling bound to
convergence in Kolmogorov distance, and the anti-concentration bound for GP in Corollary 2.1
in [14] that provides an upper bound to the pdf of Bn or equivalently the pdf of ZbnM , we can
reach


 M

 M

log2/3 n
b
c |),
P Z n ≤ t − P Z n ≤ t ≤ C3
+ γ A(|W
γ 1/3 (nh)1/6
bn (x, x) of the GP W
c is
c |) = E[ZnM ] < ∞ is a constant since the variance function C
where A(|W
uniformly bounded by some constant independent of (n, h). Consequently, we obtain that for
any t ≥ 0,


 M

en ≤ t
P ZnM ≤ t − P Z
(i) 
√



≤ P ZnM ≤ t − P ZbnM ≤ t + nh δn + n−10
√
√
√






 M

en ≤ t + nh δn + n−10
≤ P ZnM ≤ t − P ZnM ≤ t + nh δn + P ZnM ≤ t + nh δn − P Z


(ii) √

log2/3 n
≤ 4 nh δn A(|W |) + 1 + C3
+ γ A(|W |) + n−10
γ 1/3 (nh)1/6


p
√
log2/3 n
2α
+ γ + nh γn h + γn σ log n + n−10 ,
≤ C4
γ 1/3 (nh)1/6
√
bM − Z
eM ≤ nh δn with
where step (i) follows since the second display in (27) implies Z
n
n
probability at least 1 − n−10 , and step (ii) follows by applying the anti-concentration bound
c . Similarly, we can obtain the same
in Corollary 2.1 in [14] for the suprema ZnM of the GP W
 M

 M

e
bound for P Zn ≤ t − P Zn ≤ t , which combined with above implies that for any t ≥ 0,
√



 M

en ≤ t .
P ZnM ≤ t − P Z
by choosing γ =

5.8

√

log n
,
(nh)1/8

(56)

log n/(nh)1/8 .

Proof of Theorem 3.4

We need to apply the following Gaussian process comparison inequality, which is of independent interest. A proof is provided in § 5.10.

29

5.1 (Comparison inequality for Gaussian processes). Let kXkF = supf ∈F Xf . Consider two centered Gaussian processes GP(0, C X ) and GP(0, C Y ) over an index set F. Suppose
|C X (f, g) − C Y (f, g)| ≤ γ for any (f, g) ∈ F 2 . Let ρX and ρY denote their respective intrinsic
pseudometrics. If we let ZX : = kXkF and ZX : = kY kF , then for any t ∈ R and κ ∈ (0, 1),
n
o




P ZX ≤ t − P ZY ≤ t ≤ νX + νY + log N (ε, F, ρX ) + log N (ε, F, ρY ) γ 1/3
p
+ νX δX + νY δY + 2ε 2 log(4/κ) + κ,
Theorem

where
νX

Z
q
X
= K C (f, f ) + K

ρX (F ) q

log N (u, F, ρX ) du + 1

0
Z ρY (F ) q

q
C Y (f, f ) + K
log N (u, F, ρY ) du + 1
0
Z εq
p
δX = K
log N (u, F, ρX ) du + ε 2 log(4/κ)
Z0 ε q
p
δY = K
log N (u, F, ρY ) du + ε 2 log(4/κ).
νY = K

0

As a remark, for a GP over T = [0, 1], we can expand the index set to TE = T ∪ [−2, −1]
and define Xt = −Xt+2 for any t ∈ [−2, −1]. For this expanded index set, the conditions
of Theorem 5.1 still hold. In addition, this simple modification to the original GPs implies a
Kolmogorov distance bound between their sup norms, that is, supt∈T |Xt | = supt∈TE Xt and
supt∈T |Yt | = supt∈TE Yt .
e B and C Y = C
b B for
We will use this Gaussian process comparison inequality with C X = C
n
n
proving the desired approximation error bound. Due to the sup norm error bound (30), we can
choose γ  γn in Theorem 5.1. Before applying the Gaussian process comparison
inequality, we
√ q 

0
B
e
also need study properties of intrinsic pseudometric ρbn (x, x ) = nh E (f (x) − feB (x0 ))2
√
of the GP nh feB by comparing it with the Euclidean metric over X = [0, 1] (so that we can
control the covering entropy of X relative to the intrinsic pseudometric). It is easy to see that
we can express ρbn as

 

σ −2 λ h−1 ρbn (x, x0 )2 = gx,x0 (x) − gx,x0 (x0 ) − gbx,x0 (x) − gbx,x0 (x0 ) ,
(57)
where gx,x0 (·) = K(x, ·) − K(x0 , ·) and gbx,x0 is the solution of the following KRR with noiseless
observations
 X

n
1
x,x0
2
2
gbx,x0 = argmin
(Zi − g(Xi )) + λ kgkH .
n
g∈H
i=1

0

where Zix,x = K(x, Xi ) − K(x0 , Xi ) = gx,x0 (Xi ). We can similarly apply the expansion (54) to
obtain that with probability at least 1 − n−10 ,
gbx,x0 − gx,x0 + Pλ gx,x0

∞

30

≤ C γn kPλ gx,x0 k∞ .

This implies that for all sufficiently large n,
gbx,x0 − gx,x0

∞

≤ 2 kPλ gx,x0 k∞
∞ 
X
ψj (x) − ψj (x0 ) ψj (x̃)
≤ 2λ sup
1 + h2α /µj
x̃
j=1

0

≤ C λ Lψ Cψ |x − x |

∞
X
j=1

−2

. λh

j
1 + (jh)2α

0

|x − x |.

Combining this with identity (57), we obtain
σ −2 λ h−1 ρbn (x, x0 )2 . 2λ h−2 |x − x0 |,
or

ρbn (x, x0 ) . h−1 |x − x0 |1/2 ,

∀(x, x0 ) ∈ X 2 .

This implies that the covering entropy of X = [0, 1] relative to metric ρX = ρbn in Theorem 5.1
satisfies log N (u, F, ρX ) . log(n/u) for any u > 0. Similarly, it is straightforward to verify that
the same bound holds for ρY as the intrinsic pseudometric associated with the population level
GP W B . Moreover, it is easy to see that we can bound ρX (F) = supx,x0 ρbn (x, x0 ) as
σ −2 λ h−1 ρX (F) ≤ 4 sup gbx,x0 − gx,x0
x,x0

∞

∞
X
ψj (x) − ψj (x0 ) ψj (x̃)
1 + h2α /µj



≤ 4λ sup
x̃

≤ C λ Cψ2

j=1
∞
X
j=1

−1

. λh

1
1 + (jh)2α

,

implying ρX (F) . 1. Similarly, it can be verified that ρY (F) . 1. Now we can apply the
1/3
Gaussian comparison Theorem 5.1 with γ  γn , κ  n−1 and ε  γn to obtain that for any
t ≥ 0,
 B



e ≤ t Dn − P Z B ≤ t . γ 1/3 log n.
P Z
n
n
n

5.9

(58)

Proof of Theorem 3.5

enB , we obtain by
By the definition
of rn (β) as the β-th posterior quantile of (nh)−1/2 Z
√
choosing t = nh rn (β) in inequality (36) that
√




P ZnB ≤ nh rn (β) − P ZnB ≤ qnB (β) . γn1/3 log n,


where we used the definition of qnB (β) that implies P ZnB ≤ qnB (β) = β. Now the anticoncentration inequality (Corollary 2.1 in [14]) for the supremum ZnB of the GP W B implies that
the pdf of ZnB is bounded by some constant independent of (n, h) because A(|W B |) = E[ZnB ] < ∞
is a constant, we obtain that
√

nh rn (β) − qnB (β) . γn1/3 log n,

31

implying the first desired bound comparing the quantiles.
Now we proceed to prove the second and third displayed inequalities in the theorem. By
the definition, we have






P f ∗ ∈ CBn (β) = P kfbn − f ∗ k∞ ≤ rn (β) = P kfbn − Fλ f ∗ − Pλ f ∗ k∞ ≤ rn (β)
√
If the bias term kPλ f ∗ k∞ satisfies nh kPλ f ∗ k∞ → ∞ as n → ∞, then for arbitrarily large
fixed constant M > 0, as long as n is sufficiently large, we always have the bound
√






P f ∗ ∈ CBn (β) ≤ P kfbn − Fλ f ∗ k∞ ≥ kPλ f ∗ k∞ − rn (β) ≤ P ZenM ≥ M nh rn (β) .
Combining this with inequality (37), Theorem 3.4, and the Borell’s inequality for GP, we obtain
that for n large enough,


 M

e ≥ M q B (β)/2 ≤ e−C M 2 ,
P f ∗ ∈ CBn (β) ≤ P Z
n
n
for some constant C independent of (n, h). This implies the
√last desired bound (39).
On the other hand, if the bias term kPλ f ∗ k∞ satisfies nh kPλ f ∗ k∞ → 0 as n → ∞, then
we have





 M √
e ≥ nh (rn (β) + kPλ f ∗ k∞ ) .
P f ∗ ∈ CBn (β) ≤ P kfbn − Fλ f ∗ k∞ ≤ rn (β) + kPλ f ∗ k∞ = P Z
n
Now using Theorem 3.4 and the anti-concentration inequality (Corollary 2.1 in [14]) for the
supremum ZnB of the GP W B lead to
√
√



 √
log n
P f ∗ ∈ CBn (β) − P ZbnM ≥ nh rn (β) . nh kPλ f ∗ k∞ +
.
(nh)1/8
√
 M



bn ≥ nh rn (β) − P f ∗ ∈ CBn (β) .
A similar analysis leads to the same upper bound for P Z
These two together imply
√
√

 ∗

 M √
log n
∗
b
P f ∈ CBn (β) − P Zn ≥ nh rn (β) . nh kPλ f k∞ +
.
(nh)1/8
Then the desired inequality (38) is a direct consequence of the preceding display, inequality (37)
bnM = kW
c k∞
and the anti-concentration inequality (Corollary 2.1 in [14]) for the supremum Z
c |) = E[ZbM ] < ∞ is a constant. Last but not
that implies a bounded pdf of ZbnM since A(|W
n
bn of the GP W
c is uniformly strictly smaller
least, since by construction the covariance function C
B
B
b
than the covariance function Cn of the GP W , we can conclude by a classical GP comparison
c k∞ is stochastically smaller than kW
c B k∞ . Since q B (β) is by definition the
inequality that kW
n




c B k∞ , we must have P kW
c k∞ ≤ q B (β) > P kW
c B k∞ ≤ q B (β) = β.
β-th quantile of kW
n
n

5.10

Proof of Theorem 5.1

Let kXkFε = supf, g∈F , ρ(f, g)≤ε |Xf − Xg |.
Step 1: Fix ε > 0. Let N = N (ε, F, ρX ) + N (ε, F, ρY ), and {f1 , . . . , fN } ⊂ F be the
union of an ε-net of F relative to ρX and an ε-net of F relative to ρY . It is easy to see that
{f1 , . . . , fN } serves as an ε-net of F relative to both ρX and ρY . To simplify the notation, we
denote Xfj and Yfj by Xj and Yj respectively, for j = 1, 2, . . . , N , and denote CfXj ,fk and CfYj ,fk
X and C Y for j, k = 1, . . . , N . Define
by Cjk
jk
ε
ZX
= max Xj
1≤j≤N

and ZYε = max Yj ,
1≤j≤N

32

ε − Z | ≤ kXk
ε
then |ZX
X
Fε and |ZY − ZY | ≤ kY kFε . By applying Borell-Tsirelson-IbragimovSudakov inequality [47, Proposition A.2.1], we obtain that with probability at least 1 − κ/2,
p
kXkFε ≤ E[kXkFε ] + ε 2 log(4/κ), and
p
kY kFε ≤ E[kY kFε ] + ε 2 log(4/κ).

By the maximal inequality for GP [47, Corollary 2.2.8], we obtain
Z εq
log N (u, F, ρX ) du, and
E[kXkFε ] ≤ K
0
Z εq
E[kY kFε ] ≤ K
log N (u, F, ρY ) du.
0

Putting pieces together, we obtain that with probability at least 1 − κ/2,
Z εq
p
ε
|ZX − ZX | ≤ K
log N (u, F, ρX ) du + ε 2 log(4/κ) = δX ,
Z 0ε q
p
|ZYε − ZY | ≤ K
log N (u, F, ρY ) du + ε 2 log(4/κ) = δY .

and
(59)

0

Step 2: For any t > 0, we approximate the non-smooth map x 7→ 1(max1≤j≤N xj ≤ t)
by a smooth function. Following [13], we first approximate the map x 7→ max1≤j≤N xj by
PN βxj 
. A straightforward calculation gives
Fβ : RN → R, defined as Fβ (x) = β −1 log
j=1 e
max xj ≤ Fβ (x) ≤ max xj + β −1 log N.
1≤j≤N

1≤j≤N

(60)


Then we approximate the step function x 7→ 1(x ≤ t) by g(x) = v ψ(x − t) for some (large)
ψ > 0 and any smooth non-increasing function v ∈ C 3 (R) satisfying v(t) = 1 for t ≤ 0,
v(t) ∈ [0, 1] for t ∈ [0, 1], v(t) = 0 for t ≥ 1 and max{kv 0 k∞ , kv 00 k∞ k} ≤ C for some universal
constant C > 0. It is straightforward to verify that
1(x ≤ t) ≤ g(x) ≤ 1(x ≤ t + ψ −1 ),

∀x ∈ R,

(61)

and kg 0 k∞ ≤ C ψ, kg 00 k∞ ≤ C ψ 2 . Combining these two, we may approximate x 7→ 1(max1≤j≤N xj ≤
t) by the smooth function f = g ◦ Fβ . Our construction is closely related to the construction
in [13] (they consider approximating the indicator function 1A (x) for any measurable subset A,
therefore their construction is more complicated), and using their Lemma 4.3 we obtain
N X
N
X


|∂j ∂k (g ◦ Fβ )(x)| ≤ kg 00 k∞ + 2kg 0 k∞ β ≤ C ψ 2 + 2ψβ .

(62)

j=1 k=1

Step 3: We need the following multivariate version of Stein’s lemma, which can be proved
by applying integration by parts.
5.4. If F : RN → R is a C 1 function with at most polynomial growth at infinity, and
W = (W1 , . . . , Wn ) is a centered Gaussian random vector, then for any 1 ≤ i ≤ N ,

Lemma

N

 X


E Wi F (W ) =
E[Wi Wj ] E ∂i F (W ) .
i=1

Let X =√
(X1 , . . .√, XN ) and Y = (Y1 , . . . , YN ). We consider
smart path interpola the Slepian

tion Z(t) = t X + 1 − t Y for t ∈ [0, 1], and let φ(t) = E g ◦ Fβ (Zt ) . Then φ is differentiable,
33

and

X
N
 Y
Xi 
j
.
φ (t) = E
∂j (g ◦ Fβ )(Zt ) √ − √
2 t 2 1−t
j=1
0

For each j, Lemma 5.4 implies
N

 √ X


Y
E ∂j (g ◦ Fβ )(Zt )Yj = t
Cjk
E ∂jk (g ◦ Fβ )(Zt ) ,

and

k=1
N
X

 √


X
E ∂j (g ◦ Fβ )(Zt )Xj = 1 − t
Cjk
E ∂jk (g ◦ Fβ )(Zt ) .
k=1

Combining the last three displays, we obtain that for any t ∈ [0, 1],
N N
 X
1 XX 
Y
φ (t) =
E ∂jk (g ◦ Fβ )(Zt ) (Cjk
− Cjk
).
2
0

j=1 k=1

Combining this and inequality (62), we finally reach




E g ◦ Fβ (Y ) − E g ◦ Fβ (X)
Z 1
C
C
Y
X
≤ (ψ 2 + 2ψβ) γ,
− Cjk
≤
|φ0 (t)| dt ≤ (ψ 2 + 2ψβ) Cjk
2
2
0

(63)

X − C Y ≤ γ.
where in the last step we used the condition Cjk
jk

Step 4: Combining pieces together, we obtain that for any t ∈ R (notice that g is nonincreasing)
P[ZX ≤ t − δX − β −1 log N ] − P[ZY ≤ t + ψ −1 + δY ]
p
ε
≤ P[ZX
+ β −1 log N ≤ t] − P[ZYε ≤ t + ψ −1 ] + 2ε 2 log(4/κ) + κ
p




≤ E g(X + β −1 log N ) − E g(Y ) + 2ε 2 log(4/κ) + κ
p




≤ E g ◦ Fβ (X) − E g ◦ Fβ (Y ) + 2ε 2 log(4/κ) + κ
p
C
≤ (ψ 2 + 2ψβ) γ + 2ε 2 log(4/κ) + κ.
2

(by (59))
(by (61))

(64)

(by (60))
(by (63))

We have the anti-concentration bound for suprema of separable GP [14, Theorem 2.1],


P |ZX − t| ≤  ≤ 4 E[ZX ] + 4, ∀ > 0,
and the maximal inequality for GP [47, Corollary 2.2.8],
ρX (F ) q

Z
q
X
E[ZX ] ≤ K C (f, f ) + K

log N (u, F, ρX ) du,

0

where ρX (F) = supf,g∈F ρX (f, g) denotes the diameter of F relative to ρX . Combining the two
preceding displays, we obtain
n q


P |ZX − t| ≤  ≤ 4 K C X (f, f ) + K

Z
0

34

ρX (F ) q

o
log N (u, F, ρX ) du + 1 = νX .

A similar anti-concentration bound holds for ZY as
n q


P |ZY − t| ≤  ≤ 4 K C Y (f, f ) + K

Z

ρY (F ) q

o
log N (u, F, ρY ) du + 1 = νY .

0

Putting pieces together, we finally reach that for any t ∈ R,
P[ZX ≤ t] − P[ZY ≤ t]
≤ P[ZX ≤ t − δX − (2β)−1 log N − (2ψ)−1 ] − P[ZY ≤ t + δY + (2β)−1 log N + (2ψ)−1 ]
(i)


≤ νX δX + (2β)−1 log N + (2ψ)−1 + νY δY + (2β)−1 log N + (2ψ)−1
p
C
+ (ψ 2 + 2ψβ) γ + 2ε 2 log(4/κ) + κ,
2
where in step (i) we applied inequality (64) with t = t + (2β)−1 log N − (2ψ)−1 . Now, we can
take β = ψ log N and ψ = γ −1/3 to obtain
p
P[ZX ≤ t] − P[ZY ≤ t] ≤ νX δX + νY δY + (νX + νY + log N ) γ 1/3 + 2ε 2 log(4/κ) + κ.

5.11

Proof of Lemma 5.3

To begin with, we collect a version of the classical Bernstein inequality and a resulting
expectation bound.
Bernstein’s inequality: Let
random variables. Let ν and c be
PnX1 , . . .2, Xn be independent
Pn
q
−1
q−2 for q ≥ 3. Let
positive
numbers
such
that
EX
≤
ν,
and
E|X
i | ≤ 2 q! νc
i
i=1
i=1
Pn
S = i=1 (Xi − EXi ). Then, for any λ ∈ (0, 1/c),


νλ2
λS
Ee ≤ exp
,
2(1 − cλ)
which in particular implies, for any x > 0,
√
P[|S| ≥ 2νx + cx] ≤ 2e−x .

(65)

Further, if S a random variable satisfying (65), then Then, for any q ≥ 1,
E[S 2q ] ≤ q!(8ν)q + (2q!)(4c)2q .
In particular, E(S 2 ) . C1 ν + C2 c2 , where C1 , C2 are absolute constants.
Our proof makes use of the following tail bound for supremum of empirical processes with
sub-exponential increments from [2].
Bernstein-type inequality for suprema of random processes (Theorem 2.1, [2]) Let
(Xt )t∈T be a centered family with T ⊂ RD for some finite D. Fix some t0 ∈ T and let
Z̄ = supt |Xt − Xt0 |. Consider norms d(s, t) = d(s − t) and δ(s, t) = δ(s − t) on RD , and assume
there exist v, b > 0, c ≥ 0 such that
T ⊂ {t ∈ RD : d(t, t0 ) ≤ v, cδ(t, t0 ) ≤ b}.
Further, assume that for all s 6= t ∈ T ,

 λ(Xt −Xs ) 
E e
≤ exp




λ2 d2 (t, s)
1
, ∀ λ ∈ 0,
.
2(1 − λcδ(t, s))
cδ(t, s)

35

(66)

Then, with C = 18,
p


v 2 (1 + x) + b(1 + x) ≤ 2e−x , ∀ x > 0.
P Z̄ ≥ C

(67)

P
e i , t), where wi ∼ N (0, 1),
Without loss of generality, we assume σ = 1. Let Ut = ni=1 wi K(X
Xi ∼ U (0, 1), and wi ’s are independent of Xi ’s. It suffices to find a tail bound for kU k∞ =
Pn
t,s
t,s
e i , t) −
= K(X
supt∈[0,1] |Ut |. For t, s ∈ [0, 1], we write Ut − Us =
i=1 wi ∆i , where ∆i
P
∞
e i , s) =
K(X
j=1 νj {ψj (t) − ψj (s)}ψj (Xi ). We suppress the dependence on s and t, and write
∆i subsequently.
P
We now proceed to the proof. Recall Ut − Us = ni=1 wi ∆i . We first define the two norms
d and δ. Let d(t, s) = [E(Ut − Us )2 ]1/2 be the intrinsic semi-metric. For any κ ∈ (0, α ∧ 1),
let δ(t, s) = |t − s|κ . Now we estimate the quantities v, b and c appearing in Bernstein’s
inequality (67). First, we notice
d2 (t, s) = E

n
n
X
X
(wi ∆i )2 =
E∆2i .
i=1

i=1

Under Assumption B, we can obtain by using orthonormality of the eigenfunctions that
E∆2i

=

∞
X

νj2 |ψj (t)

2

− ψj (s)| ≤

2Cψ2

j=1

∞
X

νj2 ≤ Ch−1 .

j=1

Combining the two previous displays, we conclude
√
sup d(t, s) ≤ nh−1 .
t,s

In order to verify
condition (66) through characterizing the growth of the moments, we need
P
to bound E ni=1 |wi ∆i |q for q ≥ 3. To that end, we first bound |∆i | as follows. For any
i = 1, . . . , n and any fixed κ ∈ (0, 1), we use Assumption B to get
|∆i | ≤

∞
X

νj |ψj (t) − ψj (s)||ψj (Xi )|

j=1

≤ Cψ

∞
X

νj |ψj (t) − ψj (s)|κ [2Cψ1−κ ]

j=1

≤

2Cψ2−κ Lκψ

|t − s|

κ

∞
X

j κ νj ≤ Ch−(1+κ) δ(t, s).

j=1

This implies that for any q ≥ 3,
E

n
X
i=1

q

|wi ∆i | ≤

n
X



q−2
q q/2 E |∆i |q−2 ∆2i ] ≤ 2−1 q! d2 (t, s) h−(1+κ) δ(t, s)
,

i=1

q/2 and q! ≥ (q/e)q . For the |∆ |q−2 term, we used the
where we used the fact that E|wi |q . qP
i
global bound on |∆i |, and finally used ni=1 E∆2i = d2 (t, s).
Let c = h−(1+κ) . Then using the classical Bernstein’s inequality described in the beginning
of this subsection, we obtain




λ2 d2 (t, s)
1
λ(Ut −Us )
Ee
≤ exp
, ∀λ ∈ 0,
.
2(1 − λcδ(t, s))
cδ(t, s)

36

Thus the sub-exponential increment condition (66) is satisfied, and quantities v 2 and b in Bernstein’s inequality (67) are nh−1 and h−(1+κ) , respectively. Thus,
p


P kU k∞ ≥ 18 nh−1 (1 + x) + h−(1+κ) (1 + x) ≤ 2e−x , ∀ x > 0,
which implies the claimed concentration inequality by dividing both sides inside the probability
by n.

5.12

Proof of Lemma 5.2

Define
n

U (t, g) =

1X
e i , t) − E[g(X1 )K(X
e 1 , t)].
g(Xi )K(X
n
i=1

To prove the result, we will apply the peeling technique and the following lemma. A proof
of the following lemma is deferred to the next subsection. In this subsection, the meaning of
constant C can be change from line to line.
Lemma

5.5. Let Gn (ρ) = {g ∈ H : kgk∞ ≤ ρ, kgkH ≤ An }. Then


−1/2 √
−1
P
sup |U (t, g)| > Hn + (nh)
ρ x + (nh) ρ x ≤ e−x ,

x > 0,

t∈T ;g∈Gn


√
1/(2α)
1/α
where Hn = (nh)−1/2 ρ [1+ log n+An
]+(nh)−1 ρ [1+log n+An max 1, (n/h)1/(2α)−1/2 ]
for any α > 1/2.
Now we S
proceed to the proof. To apply the peeling technique, we decompose the range
−j−1 B , 2−j B ]. For any g ∈ G , there always exists some j ≥ 0 such that
(0, Bn ] into ∞
n
n
n
0
j=0 (2
q
B
−j
−1
−j
kgk∞ ∈ (2 0 Bn , 2 0 Bn ], which implies j0 ≤ c1 log kgkn for some constant c1 > 0. Using
∞

Lemma 5.5 with ρ = 2−j0 Bn and x = x + j0 , we obtain


hp
√
log n + x +

s

i
Bn
+ An1/(2α) kgk∞
kgk∞
t∈T ;g∈Gn , kgk∞ ∈(2−j0 −1 Bn , 2−j0 Bn ]

h
i
Bn
1
−1
1/α
1/(2α)−1/2
kgk∞ ≤ j0 e−x ,
log n + x + log
+ An max 1, (n/h)
+ C (nh)
kgk∞
2
√
√
√
where we used the inequality that a + b ≤ a+ b and the fact that kgk∞ > 2−j0 −1 Bn = ρ/2.
P

sup

−1/2

|U (t, g)| > C (nh)

37

log

Now by combining this with a union bound over j0 , we obtain
s

i
hp
√
Bn
log n + x + log
+ An1/(2α) kgk∞
P
sup |U (t, g)| > C (nh)−1/2
kgk∞
t∈T ;g∈Gn

h
i
Bn
1/(2α)−1/2
+ C (nh)−1 log n + x + log
+ A1/α
max
1,
(n/h)
kgk
n
∞
kgk∞
s

∞
i
hp
X
√
Bn
≤
log n + x + log
+ A1/(2α)
kgk∞
P
sup
|U (t, g)| > C (nh)−1/2
n
kgk∞
t∈T ;g∈Gn , kgk∞ ∈(2−j0 −1 Bn , 2−j0 Bn ]
j0 =0

h
i
Bn
−1
1/α
1/(2α)−1/2
+ C (nh)
log n + x + log
+ An max 1, (n/h)
kgk∞
kgk∞
∞
X
1 −x
e = 2e−x ,
≤
2 j0
j0 =0

which is the desired result.

5.13

Proof of Lemma 5.5

Let η = (t, g) ∈ T ⊗ H. The proof applies an improved version of Bernstein’s inequality (67)
over the product space T ⊗ H by truncating the chaining in the proof of Theorem 5.1 in [2] at
a finite level. We will use the following simple inequality multiple times:
e 1 , t)| ≤ C kgk h−1/2 ,
E|g(X1 )K(X
∞

for any g ∈ H.

e 1 , t)|]2 ≤ E[g 2 (X1 )K
e 2 (X1 , t)] ≤ kgk2 EK
e 2 (X1 , t) ≤
In fact, this follows since [E|g(X1 )K(X
∞
2
−1
C kgk∞ h .
First, set d2 (η, η 0 ) = E[U (t, g) − U (t0 , g 0 )]2 to be the intrinsic semi-metric in Bernstein’s
inequality (67) as in the proof of Lemma 5.3. Let’s try to get a global bound for d first for
finding v,
d2 (η, η 0 ) = Var[U (t, g) − U (t0 , g 0 )]


≤ 2 VarU (t, g) + VarU (t0 , g 0 )

2 
e 2 (X1 , t) + g 02 (X1 )K
e 2 (X1 , t0 )
≤ E g 2 (X1 )K
n
2h−1 (kgk2∞ + kg 0 k2∞ )
≤
.
n
As before, we need verify condition (66) by establishing
P a Bernstein type inequality for the
difference U (t, g) − U (t0 , g 0 ). Write U (t, g) − U (t0 , g 0 ) = ni=1 Wi , where




−1
0
0
0
0
e
e
e
e
Wi = n
g(Xi )K(Xi , t) − g (Xi )K(Xi , t ) − E g(X1 )K(X1 , t) − g (X1 )K(X1 , t ) .
P
Clearly, ni=1 EWi2 = d2 (η, η 0 ). We first obtain a bound on |Wi | in order to find quantity c and
e i , t)−Eg(X1 )K(X
e 1 , t).
metric δ in (66). To simplify the notation, we write Vg,t (Xi ) = g(Xi )K(X
n|Wi | = |Vg,t (Xi ) − Vg0 ,t0 (Xi )| ≤ |Vg,t (Xi ) − Vg,t0 (Xi )| + |Vg,t0 (Xi ) − Vg0 ,t0 (Xi )|.
We bound the two terms on the right hand side of the above display separately. Under As-

38

sumption B, we have for any fixed κ ∈ (0, 1),
e i , t) − K(X
e i , t0 )} + E g(X1 ){K(X
e 1 , t) − K(X
e 1 , t0 )}
|Vg,t (Xi ) − Vg,t0 (Xi )| ≤ g(Xi ){K(X
1/2

2
−(1+2κ)
0 2κ
−(1+κ)
0 κ
|t − t |
≤ kgk∞ h
|t − t | + kgk∞ h
. kgk∞ h−(1+κ) |t − t0 |κ .
In addition, under the same assumption we have
e i , t0 ) + E (g − g 0 )(X1 )K(X
e 1 , t0 )}
|Vg,t0 (Xi ) − Vg0 ,t0 (Xi )| ≤ (g − g 0 )(Xi )K(X
1/2

−1
−1
0 2
0
≤ g−g ∞h + g−g ∞h
. g − g0

∞

h−1 .

Putting these together, we obtain
|Wi | . (nh)−1



g − g0

∞


+ h−κ kgk∞ |t − s|κ .

This means we can choose
δ(η, η 0 ) =



g − g0

∞


+ h−κ kgk∞ |t − s|κ ,

(68)

and c = (nh)−1 . Under this choice of c and δ, we can easily verify that for q ≥ 3,
E

n
X

|Wi |q ≤ [(nh)−1 δ(η, η 0 )]q−2 d2 (η, η 0 ),

i=1

which implies the sub-exponential increment condition (66). To calculate covering numbers in
Theorem 5.1 of [2], we need to bound the d metric by the Euclidean metric. To that end,


d2 (η, η 0 ) ≤ 2 Var{U (t, g) − U (t0 , g)} + Var{(U (t0 , g) − U (t0 , g 0 )}
2  −(1+2κ)
2 
≤
h
kgk2∞ |t − s|2κ + h−1 g − g 0 ∞ .
n
Therefore, on the product space T ⊗G, we can choose v = (nh)−1/2 ρ, c = (nh)−1 and b = (nh)−1 ρ
in Theorem 5.1 of [2]. Moreover, using their notations, we have the following bound on the level
k set Ak for telescoping,

|Ak | ≤ max N (T ⊗ G, d, 2−k v), N (T ⊗ G, cδ, 2−k b) ,
where N (F, d, u) denotes the u-covering number of space F relative to metric d. In our case,
since by definition G is the αth order Sobolev space with radius An , we have
N (T ⊗ G, d, 2−k v) ≤ exp{(2k An )1/α }(2k /hκ )1/α ,

N (T ⊗ G, cδ, 2−k b) ≤ exp{(2k An )1/α }(2k /hκ )1/α .

Now we describe an improved version of Theorem 5.1 in [2] that we will use. We still apply
the chaining technique, but applying an improved telescoping identity by truncating at some
finite level k0 (please refer to the first display in the proof of Theorem 5.1 in §5.1 in [2], where

39

we have adopted our notation by identifying X· with U (·) and t with η = (t, g)),
U (η) − U (η0 ) =

kX
0 −1 h

i 

U (πk+1 (η)) − U (πk (η)) + U (η) − U (πk0 (η)) ,

∀η ∈ T ⊗ H,

k=0


where πk maps any point η ∈ T ⊗ H to some point πk (η) in Ak such that d η, πk (η) ≤ 2−k v

and cδ η, πk (η) ≤ 2−k b, and satisfies π0 (η0 ) = η0 . We choose k0 so that
sup
η∈T ⊗H

U (η) − U (πk0 (η)) ≤ v = (nh)−1/2 ρ,

p
which is satisfied if 2k0 = b n/hc, since
U (η) − U (η 0 ) ≤

n
X

|Wi | ≤ h−1



g − g0

∞

 p
+ h−κ kgk∞ |t − s|κ ≤ h/n cδ(η, η 0 ) ≤ v,

i=1

for η 0 = πk0 (η) by the definition of Ak0 in their Theorem 5.1 as an 2−k0 v-net of T ⊗H relative to
the cδ metric. Under such a choice of k0 , the proof of Theorem 5.1 therein leads to the following
concentration inequality for any α > 1/2 (this truncated version is critical if α ∈ (1/2, 1], since
otherwise the bound of Hn below diverges to +∞ in the original version where k0 = +∞),


√
P
sup |U (t, g)| > Hn + v x + b x ≤ e−x , x > 0,
t∈T ; g∈Gn

where
k0
k0
X
X
p

1/α
1/(2α)
−k+k/(2α)
+ b 1 + log n + An
2−k+k/α
Hn . v 1 + log n + An
2



k=0

k=0

p



1/(2α)−1/2
≤ C1 v 1 + log n + A1/(2α)
+ C2 b 1 + log n + A1/α
,
n
n max 1, (n/h)
which proves the desired inequality.

References
[1] Sudipto Banerjee, Alan E Gelfand, Andrew O Finley, and Huiyan Sang. Gaussian predictive
process models for large spatial data sets. Journal of the Royal Statistical Society: Series
B (Statistical Methodology), 70(4):825–848, 2008.
[2] Yannick Baraud. A bernstein-type inequality for suprema of random processes with applications to model selection in non-Gaussian regression. Bernoulli, 16(4):1064–1085, 2010.
[3] Anirban Bhattacharya and Debdeep Pati. Posterior contraction in gaussian process regression using wasserstein approximations. Information and Inference, page To Appear,
2017.
[4] Anirban Bhattacharya, Debdeep Pati, and David Dunson. Anisotropic function estimation
using multi-bandwidth gaussian processes. Annals of statistics, 42(1):352, 2014.
[5] Lawrence D Brown and Mark G Low. A constrained risk inequality with applications to
nonparametric functional estimation. The Annals of Statistics, 24(6):2524–2535, 1996.
[6] Evgeny Burnaev and Maxim Panov. Adaptive design of experiments based on gaussian
processes. 2015.
40

[7] Cristina Butucea. Exact adaptive pointwise estimation on sobolev classes of densities.
ESAIM: Probability and Statistics, 5:1–31, 2001.
[8] Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares
algorithm. Foundations of Computational Mathematics, 7(3):331–368, 2007.
[9] Ismaël Castillo. On Bayesian supremum norm contraction rates. The Annals of Statistics,
42(5):2058–2091, 2014.
[10] Ismaël Castillo and Richard Nickl. Nonparametric Bernstein–von Mises theorems in Gaussian white noise. The Annals of Statistics, 41(4):1999–2028, 2013.
[11] Ismaël Castillo, Richard Nickl, et al. On the Bernstein–von Mises phenomenon for nonparametric Bayes procedures. The Annals of Statistics, 42(5):1941–1969, 2014.
[12] Sourav Chatterjee. An error bound in the sudakov-fernique inequality. arXiv preprint
math/0510424, 2005.
[13] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato. Gaussian approximations and
multiplier bootstrap for maxima of sums of high-dimensional random vectors. The Annals
of Statistics, 41(6):2786–2819, 2013.
[14] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato. Anti-concentration and honest,
adaptive confidence bands. The Annals of Statistics, 42(5):1787–1818, 2014.
[15] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato. Gaussian approximation of
suprema of empirical processes. The Annals of Statistics, 42(4):1564–1597, 2014.
[16] Dennis D Cox. An analysis of Bayesian inference for nonparametric regression. The Annals
of Statistics, pages 903–923, 1993.
[17] Noel Cressie. Statistics for spatial data. John Wiley & Sons, 2015.
[18] Carla Currin, Toby Mitchell, Max Morris, and Don Ylvisaker. Bayesian prediction of deterministic functions, with applications to the design and analysis of computer experiments.
Journal of the American Statistical Association, 86(416):953–963, 1991.
[19] Lee H Dicker, Dean P Foster, Daniel Hsu, et al. Kernel ridge vs. principal component
regression: Minimax bounds and the qualification of regularization operators. Electronic
Journal of Statistics, 11(1):1022–1047, 2017.
[20] David Freedman. Wald lecture: On the Bernstein-von Mises theorem with infinitedimensional parameters. The Annals of Statistics, 27(4):1119–1141, 1999.
[21] Evarist Giné and Richard Nickl. Uniform central limit theorems for kernel density estimators. Probability Theory and Related Fields, 141(3-4):333–387, 2008.
[22] Evarist Giné and Richard Nickl. Rates of contraction for posterior distributions in lr metrics, 1 ≤ r ≤ ∞. The Annals of Statistics, pages 2883–2911, 2011.
[23] Marc Hoffmann, Judith Rousseau, and Johannes Schmidt-Hieber. On adaptive posterior
concentration rates. The Annals of Statistics, 43(5):2259–2295, 2015.
[24] Iain M Johnstone. High dimensional Bernstein-von Mises: simple examples. Institute of
Mathematical Statistics collections, 6:87, 2010.
[25] Marc C Kennedy and Anthony O’Hagan. Bayesian calibration of computer models. Journal
of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001.
41

[26] Bartek T Knapik, Aad W van der Vaart, and J Harry van Zanten. Bayesian inverse
problems with Gaussian priors. The Annals of Statistics, 39(5):2626–2657, 2011.
[27] Haralambie Leahu. On the Bernstein-von Mises phenomenon in the Gaussian white noise
model. Electronic Journal of Statistics, 5:373–404, 2011.
[28] Georges Matheron. The intrinsic random functions and their applications. Advances in
applied probability, 5(3):439–468, 1973.
[29] Shahar Mendelson. Geometric parameters of kernel machines. Lecture notes in computer
science, pages 29–43, 2002.
[30] Yoichi Nishiyama. Impossibility of weak convergence of kernel density estimators to a
non-degenerate law in l2 (rd ). Journal of Nonparametric Statistics, 23(1):129–135, 2011.
[31] David Pollard. Asymptotics via empirical processes. Statistical science, pages 341–354,
1989.
[32] Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine
learning. MIT Press, 2006.
[33] Kolyan Ray. Adaptive Bernstein-von Mises theorems in Gaussian white noise. arXiv
preprint arXiv:1407.3397, 2014.
[34] Mark Rudelson and Roman Vershynin. Hanson-Wright inequality and sub-Gaussian concentration. Electronic Communications in Probability, 18, 2013.
[35] Jerome Sacks, William J Welch, Toby J Mitchell, and Henry P Wynn. Design and analysis
of computer experiments. Statistical science, pages 409–423, 1989.
[36] Daniel Jonathan Scansaroli. Stochastic Modeling with Temporally Dependent Gaussian
Processes Applications to Financial Engineering, Pricing and Risk Management. Lehigh
University, 2012.
[37] Paulo Serra and Tatyana Krivobokova. Adaptive empirical Bayesian smoothing splines.
Bayesian Analysis, 12(1):219–238, 2017.
[38] Bernard W Silverman. Spline smoothing: the equivalent variable kernel method. The
Annals of Statistics, pages 898–916, 1984.
[39] Suzanne Sniekers and Aad van der Vaart. Credible sets in the fixed design model with
Brownian motion prior. Journal of Statistical Planning and Inference, 166:78–86, 2015.
[40] Michael L Stein. Interpolation of spatial data: some theory for kriging. Springer Science &
Business Media, 2012.
[41] Ingo Steinwart, Don R Hush, and Clint Scovel. Optimal rates for regularized least squares
regression. In COLT, 2009.
[42] Gilles Stupfler et al. On the weak convergence of the kernel density estimator in the uniform
topology. Electronic Communications in Probability, 21, 2016.
[43] Botond Szabó, AW van der Vaart, JH van Zanten, et al. Frequentist coverage of adaptive
nonparametric bayesian credible sets. The Annals of Statistics, 43(4):1391–1428, 2015.
[44] Aad van der Vaart and Harry van Zanten. Information rates of nonparametric gaussian
process methods. Journal of Machine Learning Research, 12(Jun):2095–2119, 2011.
42

[45] A. W. van der Vaart and J. H. van Zanten. Adaptive bayesian estimation using a gaussian
random field with inverse gamma bandwidth. Ann. Statist., 37(5B):2655–2675, 2009.
[46] Aad W van der Vaart and J Harry van Zanten. Rates of contraction of posterior distributions based on gaussian process priors. The Annals of Statistics, pages 1435–1463,
2008.
[47] Aad W Van Der Vaart and Jon A Wellner. Weak convergence. In Weak Convergence and
Empirical Processes, pages 16–28. Springer, 1996.
[48] Grace Wahba. Spline models for observational data. SIAM, 1990.
[49] Yun Yang, David B Dunson, et al. Bayesian manifold regression. The Annals of Statistics,
44(2):876–905, 2016.
[50] Yun Yang and Debdeep Pati. Bayesian model selection consistency and oracle inequality
with intractable marginal likelihood. arXiv preprint arXiv:1701.00311, 2017.
[51] Yun Yang, Mert Pilanci, Martin J Wainwright, et al. Randomized sketches for kernels: Fast
and optimal nonparametric regression. The Annals of Statistics, 45(3):991–1023, 2017.
[52] Yun Yang, Zuofeng Shang, and Guang Cheng. Non-asymptotic theory for nonparametric
testing. arXiv preprint arXiv:1702.01330, 2017.
[53] William Weimin Yoo, Subhashis Ghosal, et al. Supremum norm posterior contraction
and credible sets for nonparametric multivariate regression. The Annals of Statistics,
44(3):1069–1102, 2016.
[54] Tong Zhang. Learning bounds for kernel regression using effective data dimensionality.
Neural Computation, 17(9):2077–2098, 2005.

43

