Scalable Training of Inference Networks for Gaussian-Process Models
Appendix

A. Derivation using the Dual Representation of GP and Natural Gradient
We note that the stochastic mirror descent is equivalent to natural gradient updates for exponential family densities (Raskutti& Mukherjee, 2015; Khan & Lin, 2017). To show this, we derive the same adaptive Bayesian filter in eq. (6) using the dual
representation of GP as a Gaussian measure and natural gradient.
A Gaussian process GP(m(x), κ(x, x0 )) has a dual representation as a Gaussian measure ν on a separable Banach space
B (Cheng & Boots, 2017; Mallasto & Feragen, 2017). There is an RKHS H that corresponds to a positive definite kernel kdensely embedded in B. The measure ν is constructed as follows. First define a canonical Gaussian cylinder set measureν H on H, denoted as N (f |µ, Σ), where µ ∈ H is the mean function, Σ : H → H is a bounded positive semi-definite linear
operator. They satisfym(x) = µ> k(x, ·),κ(x, x0 ) = k(x, ·)> Σk(x0 , ·),where h> g denotes inner product in the RKHS: h> g = hh, giH . Let i be the inclusion map from H into B. Then themeasure ν is induced by ν H using this map. In the measure theory of infinite-dimensional space, ν is known as the abstract
Wiener measure. The RKHS H is sometimes called the Cameron-Martin space.
Remark 1. The intuition for this construction is that the canonical Gaussian cylinder set measure ν H is not a proper
measure. In fact, we can show that countable additivity does not hold for this ”measure” (Eldredge, 2016). The inclusion
map i here radonifies ν H into a true measure ν. One way to think about this is that functions drawn from the Gaussian
process fall outside of H with probability one (Kanagawa et al., 2018), but are contained in B. Despite this, as pointed out
in Cheng & Boots (2017), we can conveniently work with the canonical form N (f |µ, Σ) and get correct results as long as
the conclusion is independent of the dimension of f .

It is easy to check that the canonical Gaussian measure which corresponds to the GP prior GP(0, k(x, x0 )) is p(f ) =N (f |0, I). Assuming the canonical form of the variational distribution is q(f ) = N (f |µ, Σ), we have stochastic approximation of the lower bound as:L̂(q) = N Eq log p(yn |f ) − KL [q(f )kp(f )] .
Now that we can interpret GPs in the form of canonical Gaussian measures, we can then write q(f ) and p(f ) in the
exponential family form to simplify the derivation:q(f ) ∝ exp{λ> t(f ) − A(λ)},

p(f ) ∝ exp{λ>0 t(f ) − A(λ0 )}.
where t(f ) = {f, f f > } denotes the sufficient statistics, A(λ) = 12 µ> Σ−1 µ + 12 log |Σ| is the partition function. The
natural parameters of p(f ) and q(f ) are λ0 = {0, − 21 I} and λ = {Σ−1 µ, − 12 Σ−1 }, respectively. Let u denote the mean
parameters of q(f ): u = Eq [t(f )] = {µ, Σ}. There is a dual relationship between the natural parameter λ and the mean
parameter u: u = ∇λ A(λ). The mapping ∇A is one-to-one when the exponential family is minimal (Wainwright et al.,2008). The stochastic natural gradient of the lower bound with respect to λ is defined as
˜ λ L̂(q) = F (λ)−1 ∇λ L̂(q),∇where F (λ) = Eq ∇λ log q(f ) ∇λ log q(f )> is the Fisher information matrix. The natural gradient of the KL divergenceScalable Training of Inference Networks for Gaussian-Process Models

term is


q(f )
−1
˜
∇λ KL [q(f )kp(f )] = F (λ) ∇λ Eq log
p(f )


−1
= F (λ) ∇λ (λ − λ0 )> u − A(λ) + A(λ0 )= F (λ)−1 [u − ∇λ A(λ) + ∇λ u(λ − λ0 )]= F (λ)−1 ∇2λ A(λ)(λ − λ0 )= λ − λ0 ,where we have used the fact that ∇2λ A(λ) = F (λ). We can also derive a simplified form of the natural gradient of theconditional log likelihood term, by writing it as the gradient with respect to the mean parameter u:
˜ λ Eq [log p(yn |f )] = F (λ)−1 ∇λ Eq [log p(yn |f )]∇= F (λ)−1 ∇λ u · ∇u Eq [log p(yn |f )]= F (λ)−1 ∇2λ A(λ) · ∇u Eq [log p(yn |f )]= ∇u Eq [log p(yn |f )].
So the natural gradient update can be written as
λt+1 = λt + βt (N ∇u Eq [log p(yn |f )] − λt + λ0 )= (1 − βt )λt + βt (N ∇u Eq [log p(yn |f )] + λ0 ) .

Reinterpreting the above equation in the density space, we haveqt+1 (f ) ∝ qt (f )1−βt p(f )βt exp{h∇u Eq [log p(yn |f )], t(f )i}N βt .

(15)

The likelihood p(yn |f ) is said to be conjugate with the prior if it has a form as p(yn |f ) ∝ exp{λ(yn )> t(f )}. For example,
in GP regression, the likelihood is p(yn |f ) = N (yn |f (xn ), σ 2 ), which has an above form with the natural parameter
λ(yn ) = { σ12 yn k(xn , ·), − 2σ1 2 k(xn , ·)k(xn , ·)> }. By plugging in p(yn |f ) ∝ exp{λ(yn )> t(f )}, we can verify that eq. (15)is equivalent to
qt+1 (f ) ∝ qt (f )1−βt p(f )βt p(yn |f )N βt ,which turns out to be the same adaptive Bayesian filter we get in eq. (6). As for the non-conjugate case, we can view thenatural gradient update as the projection of the functional mirror descent update onto exponential families, by approximating
the likelihood term with the exponential family exp{h∇u Eq [log p(yn |f )], t(f )i}.

B. Experiment Details and Additional Results
B.1. Synthetic Data
The full figures including FBNN with M = 2, 5, 20 are shown in Fig. 4. We can see that FBNN’s problem of overestimating
uncertainty is consistent with different values of M . We set β0 = 1, ξ = 0.1 in this experiment.
B.2. Regression
For all regression experiments, we set the measurement points to be sampled from the empirical distribution of training data
convolved with the prior RBF kernel if the prior hyperparameters are initialized by optimizing GP marginal likelihood on arandom subset, otherwise we set c(x) simply to be the training data distribution. We fix c(x) and do not adapt it together
with the prior kernel parameters when the prior hyperparameters are updated during training.
Benchmarks The GP we use has a RBF kernel with dimension-wise lengthscales, also known as Automatic Relevance
Determination (ARD) (MacKay, 1996). The RFE inference network we use has 1000 features (hidden units). We initialize
the lengthscales in the network using the lengthscales of the prior kernel, and initialize the frequencies (first-layer weights)
with random samples from a standard Gaussian. We then train these frequencies and lengthscales as inference network
parameters. For FBNN, we keep all settings (including the inference network) the same as GPNet except the training

Scalable Training of Inference Networks for Gaussian-Process Models

Figure 4. Posterior process on the Snelson dataset, where shaded areas correspond to intervals of ±3 standard deviations, and dashed lines
denote the ground truths. (a) From top to bottom: SVGP with M ∈ {2, 5, 20} (left to right) inducing points; GPNet with M ∈ {2, 5, 20}
measurement points; FBNN with M ∈ {2, 5, 20} measurement points.

objective used. We use 20 random splits for each dataset, where we keep 90% of the dataset as the training set and use theremaining 10% for test. The inputs and outputs of all data points are scaled to have nearly zero mean and unit variance using
the mean and standard deviation calculated from training data. We use minibatch size 500, learning rate η = 0.003 for all
the experiments. With each random seed we ran 10K iterations. For Boston, Concrete, Kin8nm and Protein we initialize
the prior hyperparameters by maximizing GP marginal likelihood for 1K iterations on a randomly chosen subset of 1000
points. For Power and Wine we optimize the prior hyperparameters during training using the minibatch lower bound andthe same learning rate as η, as described in section 3.3. We set β0 = 1 and ξ = 1 except for Power and Protein we use
β0 = 0.01, ξ = 0.1 and β0 = 0.1, ξ = 0.1, respectively. In Table 3 and Table 4 we list the mean and standard errors of all
experiments.
Airline Delay We use the same type of GP prior and inference networks as in benchmark datasets above. We set
β0 = 0.1, ξ = 0.1. For all methods we train for 10K iterations with minibatch size 500 and learning rate η = 0.003. For
GPNet we optimize the prior hyperparameters during training with the same learning rate as η, for which the objective
is described in section 3.3, while we found that doing the same with FBNN seriously hurt its performance (leading to
RMSE 27.186 for M = 100), therefore we did not update hyperparameters during training for FBNN. We initialize theprior hyperparameters by maximizing GP marginal likelihood for 1K iterations on a randomly chosen subset of 1000 points.
We did this for both GPNet and FBNN, though we found that for GPNet this does not improve the performance.
B.3. Classification
CNN-GP Prior The prior is defined as follows. Let Z(`) (x) denote the pre-activation output of the `-th layer of the
ConvNet. The shape of Z(`) (x) is C (`) × (H (`) D(`) ). Each row of it represents the flattened feature map in a channel. A

Scalable Training of Inference Networks for Gaussian-Process Models
Table 3. Regression: RMSE.
DATA SETN

D

SVGP, 100

GPN ET, 100

SVGP, 500

GPN ET, 500

FBNN, 500

B OSTON
C ONCRETE
E NERGY
K IN 8 NM
P OWERP ROTEIN
W INE

506
1030
768
8192
9568
45730
1599

12
8
8
8
4
9
11

2.897±0.132
5.768±0.094
0.469±0.014
0.086±0.001
3.941±0.033
4.536±0.010
0.625±0.009

2.786±0.142
5.301±0.127
0.493±0.022
0.080±0.001
3.942±0.032
4.540±0.014
0.614±0.010

3.023±0.187
5.075±0.119
0.439±0.015
0.074±0.000
3.791±0.034
4.154±0.010
0.626±0.009

2.754±0.143
5.050±0.132
0.461±0.014
0.067±0.000
3.898±0.032
4.329±0.013
0.627±0.009

2.775±0.141
5.089±0.117
0.459±0.013
0.072±0.000
4.135±0.029
4.087±0.051
0.633±0.008

Table 4. Regression: Test log likelihood.
DATA SETN

D

SVGP, 100

GPN ET, 100

SVGP, 500

GPN ET, 500

FBNN, 500

B OSTON
C ONCRETE
E NERGY
K IN 8 NM
P OWERP ROTEIN
W INE

506
1030
768
8192
9568
45730
1599

12
8
8
8
4
9
11

-2.465±0.054-3.166±0.015-0.675±0.024
1.006±0.004-2.793±0.008-2.932±0.002-0.949±0.014-2.421±0.049-3.115±0.024-1.060±0.008
1.095±0.011-2.794±0.007-3.057±0.032-0.917±0.014-2.458±0.072-3.027±0.023-0.600±0.033
1.183±0.004-2.755±0.008-2.841±0.002-0.949±0.015-2.429±0.055-3.066±0.022-0.847±0.013
1.283±0.005-2.783±0.008-2.986±0.029-0.948±0.014-2.437±0.025-3.046±0.029-0.755±0.018
1.189±0.005-2.847±0.006-2.821±0.014-0.961±0.013

hidden layer in the network makes the transformation:
(`+1)
(`)
Zj,g (x) = bj +(`)(`) (`)
C
D
XHX

i=1

(`)

(`)

Wj,i,g,h a(Zi,h (x)),

h=1where Wji is the pseudo weight matrix that corresponds to the convolutional filter Uji . The elements of each row in Wji
are zero except where Uji applies. bj denotes the bias in the j-th channel. a is the ReLU activation function. Let x, y denote(`)(`)
the positions within a filter, independent Gaussian priors are placed over uj,i,x,y and bj to form a Bayesian ConvNet:(`)

2uj,i,x,y ∼ N (0, σw/C (`) ),bj ∼ N (0, σb2 ).

By carefully taking the limit of hidden-layer widths, one can prove that each row in Z(`) (x) form a multivariate Gaussian,and different rows (channels) are independent and identically distributed (i.i.d.), thus showing that the Bayesian ConvNet(`+1)
defines a GP (Garriga-Alonso et al., 2019). It is easy to show the prior mean function is zero: E[Zj,g (x)] = 0. To(`)

determine the prior covariance kernel of the output, we can follow a recursive procedure (For simplicity, we use vg (x, x0 )(`)(`)
to denote the covariance between Zj,g (x) and Zj,g (x0 )):X

2vg(`+1) (x, x0 ) = σb2 + σw(`)

sh (x, x0 ),h∈g-th patch(`)

sh (x, x0 ) = 1/2π

q

(`)

(`)

vg (x, x)vg (x0 , x0 )J1 (θg(`) ),

(`)
(`)
(`)
(`)
(`)
where J1 (θg ) = sin θg +(π −θg ) cos θg and θg = arccos



q

(`)
vg (x, x0 )/



(`)
(`)
vg (x, x)vg (x0 , x0 ). The prior convnet

we used is a deep convolutional neural network with 6 residual blocks, each two of them operates on a different size offeature maps, with the first two on feature maps with the same size as the original image. There are strided convolution(stride=2) between the three groups.
Inference Networks The inference network we used has the same structure as the prior ConvNet, except the number ofconvolutional filters are [64, 64, 128, 128, 256, 256]. On top of it we have a fully-connected layer of size 512 and neuraltangent kernels defined by a MLP with 100 hidden units for the output of each class.

Scalable Training of Inference Networks for Gaussian-Process Models

We use batch size 64 and M = 64 measurement points in this experiment. We set c(x) to be the empirical distribution of
the training data. In implementation this simply means that we use two different shuffles of the training dataset, and pick aminibatch from each of them. Then we use one of the two minibatches as training points, and the other as measurement
points. The learning rate is η = 0.0003. We set β0 = 0.01, ξ = 0.1, and ran for 10K iterations. We did not update priorhyperparameters in this experiment.