Deep kernel processes

Introduction

The deep learning revolution has shown us that effective performance on difficult tasks such as image classification (Krizhevsky et al., 2012) requires deep models with flexible lower-layers that learn task-dependent representations. Here, we consider whether these insights from the neural network literature can be applied to purely kernel-based methods. (Note that we do not consider deep Gaussian processes or DGPs to be “fully kernel-based” as they use a feature-based representation in intermediate layers).

Importantly, deep kernel methods (e.g. Cho & Saul, 2009) already exist. In these methods, which are closely related to infinite Bayesian neural networks (Lee et al., 2017; Matthews et al., 2018; Garriga-Alonso et al., 2018; Novak et al., 2018), we take an initial kernel (usually the dot product of the input features) and perform a series of deterministic, parameter-free transformations to obtain an output kernel that we use in e.g. a support vector machine or Gaussian process. However, the deterministic, parameter-free nature of the transformation from input to output kernel means that they lack the capability to learn a top-layer representation, which is believed to be crucial for the effectiveness of deep methods (Aitchison, 2019).

Contributions

We propose deep kernel processes (DKPs), which combine nonlinear transformations of the kernel, as in Cho & Saul (2009) with a flexible learned representation by exploiting a Wishart or inverse Wishart process (Dawid, 1981; Shah et al., 2014).

We show that models ranging from DGPs (Damianou & Lawrence, 2013; Salimbeni & Deisenroth, 2017) to Bayesian neural networks (BNNs; Blundell et al., 2015, App. C.1), infinite BNNs (App. C.2) and infinite BNNs with bottlenecks (App. C.3) can be written as DKPs.

We define a specific DKP, the deep inverse Wishart process (DIWP) which offers convenient variational approximate posteriors.

We develop a novel doubly-stochastic variational inducing-point inference scheme purely in the kernel domain (as opposed to Salimbeni & Deisenroth, 2017, who described DSVI for standard feature-based DGPs) for DIWPs.

We demonstrate improved performance of DIWPs on fully-connected benchmark datasets.

DKPs and specifically DIWPs offer two key advantages over feature-based methods such as DGPs and BNNs. First, DGPs and BNNs have complex approximate posteriors (Li et al., 2018), due in part to permuation/rotation symmetries in the posterior over weights/features (App. D.1 and D.2; MacKay, 1992; Moore, 2016; Pourzanjani et al., 2017). This complexity means that common variational approximate posteriors can give a very poor approximation to the true posterior. In contrast, the Gram matrices in DKPs are invariant to permutations/rotations of the weights/features and thus have much simpler true posteriors which are more easily captured by variational approximate posteriors. Second, in DIWPs the “width” parameter is learnable, and in the limit of infinite width gives a series of deterministic kernel transformations, as in an infinite neural network. This gives DIWPs the ability to learn on a layer-by-layer basis where a deterministic kernel transformation is appropriate, or where more flexibility in the kernel is needed.

Background

We briefly revise Wishart and inverse Wishart distributions. The Wishart distribution is a generalization of the gamma distribution that is defined over positive semidefinite matrices. Suppose that we have a collection of $P$-dimensional random variables $\mathbf{x}_{i}$ with $i\in\{1,\dotsc,N\}$ such that

has Wishart distribution with scale matrix $\mathbf{V}$ and $N$ degrees of freedom. When $N>P-1$, the density is,

where $\Gamma_{P}$ is the multivariate gamma function. Further, the inverse, $\mathbf{S}^{-1}$ has inverse Wishart distribution, $\mathcal{W}^{-1}\!\left(\mathbf{V}^{-1},N\right)$. The inverse Wishart is defined only for $N>P-1$ and also has closed-form density. Finally, we note that the Wishart distribution has mean $N\mathbf{V}$ while the inverse Wishart has mean $\mathbf{V}^{-1}/(N-P-1)$ (for $N>P+1$).

Deep kernel processes

We define a kernel process to be a set of distributions over positive definite matrices of different sizes, that are consistent under marginalisation (Dawid, 1981; Shah et al., 2014). The two most common kernel processes are the Wishart process and inverse Wishart process, which we write in a slightly unusual form to ensure their expectation is $\mathbf{K}$. We take $\mathbf{G}$ and $\mathbf{G}^{\prime}$ to be finite dimensional marginals of the underlying Wishart and inverse Wishart process,

and where we explicitly give the consistent marginal distributions over $\mathbf{K}^{*}$, $\mathbf{G}^{*}$ and $\mathbf{G}^{\prime*}$ which are $P^{*}\times P^{*}$ principal submatrices of the $P\times P$ matrices $\mathbf{K}$, $\mathbf{G}$ and $\mathbf{G}^{\prime}$ dropping the same rows and columns. In the inverse-Wishart distribution, $\delta$ is a positive parameter that can be understood as controlling the degree of variability, with larger values for $\delta$ implying smaller variability in $\mathbf{G}^{\prime}$.

We define a deep kernel process by analogy with a DGP, as a composition of kernel processes, and show in App. A that under sensible assumptions any such composition is itself a kernel process. Note that we leave the question of the full Kolmogorov extension theorem (Kolmogorov, 1933) for matrices to future work: for our purposes, it is sufficient to work with very large but ultimately finite input spaces as in practice, the input vectors are represented by elements of the finite set of 32-bit or 64-bit floating-point numbers (Sterbenz, 1974).

For DGP regression, the outputs, $\mathbf{Y}$, are most commonly given by a likelihood that can be written in terms of the output features, $\mathbf{F}_{L+1}$. For instance, for regression, the distribution of the $\lambda$th output feature column could be

alternatively, we could use a classification likelihood,

Importantly, our methods can be used with any likelihood with a known probability density function.

Except at the output, the model is phrased entirely in terms of positive-definite kernels and Gram matrices, and is consistent under marginalisation (assuming a valid kernel function) and is thus a DKP. At a high level, the model can be understood as alternatively sampling a Gram matrix (introducing flexibility in the representation), and nonlinearly transforming the Gram matrix using a kernel (Fig. 2).

This highlights a particularly simple interpretation of the DKP as an autoregressive process. In a standard autoregressive process, we might propagate the current vector, $\mathbf{x}_{t}$, through a deterministic function, $\mathbf{f}(\mathbf{x}_{t})$, and add zero-mean Gaussian noise, $\boldsymbol{\xi}$,

By analogy, the next Gram matrix has expectation centered on a deterministic transformation of the previous Gram matrix,

Remember that we are restricted to kernels that can be written as a function of the Gram matrix,

Superficially, this might seem very limiting — leaving us only with dot-product kernels (Williams & Rasmussen, 2006) such as,

However, in reality, a far broader range of kernels fit within this class. Importantly, isotropic or radial basis function kernels including the squared exponential and Matern depend only on the squared distance between points, $R$, (Williams & Rasmussen, 2006)

These kernels can be written as a function of $\mathbf{G}$, because the matrix of squared distances, $\mathbf{R}$, can be computed from $\mathbf{G}$,

Variational inference in deep kernel processes

A key part of the motivation for developing deep kernel processes was that the posteriors over weights in a BNN or over features in a deep GP are extremely complex and multimodal, with a large number of symmetries that are not captured by standard approximate posteriors (MacKay, 1992; Moore, 2016; Pourzanjani et al., 2017). For instance, in the Appendix we show that there are permutation symmetries in the prior and posteriors over weights in BNNs (App. D.1) and rotational symmetries in the prior and posterior over features in deep GPs with isotropic kernels (App. D.2). The inability to capture these symmetries in standard variational posteriors may introduce biases in the parameters inferred by variational inference, because the variational bound is not uniformly tight across the state-space (Turner & Sahani, 2011). Gram matrices are invariant to permutations or rotations of the features, so we can sidestep these complex posterior symmetries by working with the Gram matrices as the random variables in variational inference. However, variational inference in deep Wishart processes (equivalent to DGPs Sec. 4.1 and infinite NNs with bottlenecks App. C.3) is difficult because the approximate posterior we would like to use, the non-central Wishart (App. E), has a probability density function that is prohibitively costly and complex to evaluate in the inner loop of a deep learning model (Koev & Edelman, 2006). Instead, we consider an inverse Wishart process prior, for which the inverse Wishart itself makes a good choice of approximate posterior.

The resulting inverse Wishart process does not have a direct interpretation as e.g. a deep GP, but does have more appealing properties for variational inference, as it is always full-rank and allows independent control over the approximate posterior mean and variance. Finally, it is important to note that Wishart and inverse Wishart distributions do not differ as much as one might expect; the standard Wishart and standard inverse Wishart distributions have isotropic distributions over the eigenvectors so they only differ in terms of their distributions over eigenvalues, and these are often quite similar, especially if we consider a Wishart model with ResNet-like structure (App. H).

2 An approximate posterior for the deep inverse Wishart process

Inspired by this exact posterior that is available in simple models, we choose the approximate posterior in our model to be,

3 Doubly stochastic inducing-point variational inference in deep inverse Wishart processes

where the $\mathbf{X}$-dependence enters on the right because $\mathbf{G}_{1}=\tfrac{1}{N_{0}}\mathbf{X}\mathbf{\Omega}\mathbf{X}^{T}$. Taking inspiration from Salimbeni & Deisenroth (2017), the full approximate posterior is the product of an approximate posterior over inducing points and the conditional prior for train/test points,

and the prior can be written in the same form,

To obtain the full ELBO, we substitute Eqs. (25) and (26) into Eq.(23), the conditional prior terms cancel,

Importantly, the first term in the ELBO (Eq. 27 is a summation across test/train datapoints, and the second term depends only on the inducing points, so as in Salimbeni & Deisenroth (2017) we can compute unbiased estimates of the expectation by taking only a minibatch of datapoints, and we never need to compute the density of the conditional prior (Eq. 30), we only need to be able to sample it.

Finally, to sample the test/training points, conditioned on the inducing points, we need to sample,

and where $P_{\text{i}}$ is the number of inducing inputs, and $P_{t}$ is the number of train/test inputs. Defining the Schur complements,

The full algorithm is given in Alg. 1, where the $\operatorname{P}$ and $\operatorname{Q}$ distributions for $\mathbf{\Omega}$ and for inducing points are given by Eq. (18) and (21). We optimize using standard reparameterised variational inference (Kingma & Welling, 2013; Rezende et al., 2014) (for details on how to reparameterise samples from the Wishart, see Ober & Aitchison, 2020).

Computational complexity

As in non-deep GPs, the complexity is $\mathcal{O}(P^{3})$ for time and $\mathcal{O}(P^{2})$ for space for standard DKPs (the $\mathcal{O}(P^{3})$ time dependencies emerge e.g. because of inverses and determinants required for the inverse Wishart distributions). For DSVI, there is a $P_{\text{i}}^{3}$ time and $P_{\text{i}}^{2}$ space term for the inducing points, because the computations for inducing points are exactly the same as in the non-DSVI case. As we can treat each test/train point independently (App. I), the complexity for test/training points must scale linearly with $P_{\text{t}}$, and this term has $P_{\text{i}}^{2}$ time scaling, e.g. due to the matrix products in Eq. (32). Thus, the overall complexity for DSVI is $\mathcal{O}(P_{\text{i}}^{3}+P_{\text{i}}^{2}P_{\text{t}})$ for time and $\mathcal{O}(P_{\text{i}}^{2}+P_{\text{i}}P_{\text{t}})$ for space which is exactly the same as non-deep inducing GPs. Thus, and exactly as in non-deep inducing-GPs, by using a small number of inducing points, we are able to convert a cubic dependence on the number of input points into a linear dependence, which gives considerably better scaling.

Surprisingly, this is substantially better than standard DGPs. In standard DGPs, we allow the approximate posterior covariance for each feature to differ (Salimbeni & Deisenroth, 2017), in which case, we are in essence doing standard inducing-GP inference over $N$ hidden features, which gives complexity of $\mathcal{O}(NP_{\text{i}}^{3}+NP_{\text{i}}^{2}P_{\text{t}})$ for time and $\mathcal{O}(NP_{\text{i}}^{2}+NP_{\text{i}}P_{\text{t}})$ for space (Salimbeni & Deisenroth, 2017). It is possible to improve this complexity by restricting the approximate posterior to have the same covariance for each point (but this restriction harms performance).

Results

We began by comparing the performance of our deep inverse Wishart process (DIWP) against infinite Bayesian neural networks (known as the neural network Gaussian process or NNGP) and DGPs. To ensure sensible comparisons against the NNGP, we used a ReLU kernel in all models (Cho & Saul, 2009). For all models, we used three layers (two hidden layers and one output layer), with three applications of the kernel. In each case, we used a learned bias and scale for each input feature, and trained for 8000 gradient steps with the Adam optimizer with 100 inducing points, a learning rate of $10^{-2}$ for the first 4000 steps and $10^{-3}$ for the final 4000 steps. For evaluation, we used approximate posterior 100 samples, and for each training step we used 10 approximate posterior samples in the smaller datasets (boston, concrete, energy, wine, yacht), and 1 in the larger datasets.

Next, we considered fully-connected networks for small image classification datasets (MNIST and CIFAR-10; Table 2). We used the same models as in the previous section, with the omission of learned bias and scaling of the inputs. Note that we do not expect these methods to perform well relative to standard methods (e.g. CNNs) for these datasets, as we are using fully-connected networks with only 100 inducing points (whereas e.g. work in the NNGP literature uses the full $60,000\times 60,000$ covariance matrix). Nonetheless, as the architectures are carefully matched, it provides another opportunity to compare the performance of DIWPs, NNGPs and DGPs. Again, we found that DIWP usually gave statistically significant gains in predictive performance (except for CIFAR-10 test-log-likelihood, where DIWP lagged by only $0.01$). Importantly, DIWP gives very large improvements in the ELBO, with gains of $0.09$ against DGPs for MNIST and $0.08$ for CIFAR-10 (Table 2). For MNIST, remember that the ELBO must be negative (because both the log-likelihood for classification and the KL-divergence term give negative contributions), so the change from $-0.301$ to $-0.214$ represents a dramatic improvement.

Related work

Our first contribution was the observation that DGPs with isotropic kernels can be written as deep Wishart processes as the kernel depends only on the Gram matrix. We then gave similar observations for neural networks (App. C.1), infinite neural networks (App. C.2) and infinite network with bottlenecks (App. C.3, also see Aitchison, 2019). These observations motivated us to consider the deep inverse Wishart process prior, which is a novel combination of two pre-existing elements: nonlinear transformations of the kernel (e.g. Cho & Saul, 2009) and inverse Wishart priors over kernels (e.g. Shah et al., 2014). Deep nonlinear transformations of the kernel have been used in the infinite neural network literature (Lee et al., 2017; Matthews et al., 2018) where they form deterministic, parameter-free kernels that do not have any flexibility in the lower-layers (Aitchison, 2019). Likewise, inverse-Wishart distributions have been suggested as priors over covariance matrices (Shah et al., 2014), but they considered a model without nonlinear transformations of the kernel. Surprisingly, without these nonlinear transformations, the inverse Wishart prior becomes equivalent to simply scaling the covariance with a scalar random variable (App. L; Shah et al., 2014).

In addition, there are generalised Wishart processes (Wilson & Ghahramani, 2010, contrasting with our deep Wishart processes). While the term “generalised Wishart process” is not yet in widespread use, it allows us to make a distinction that is very useful in our context. In particular, a generalised Wishart process is a distribution over infinitely many finite-dimensional marginally Wishart matrices. For instance, these might represent the noise in a dynamical system. In that case, there would in principle be infinitely covariance matrices, one for each state-space location or time-point (Wilson & Ghahramani, 2010; Heaukulani & van der Wilk, 2019; Jorgensen et al., 2020). In contrast, kernel processes (Dawid, 1981; Bru, 1991) are distributions over a single infinite dimensional matrix. We stack these kernel process to form a (non-genearlised) deep kernel process. Importantly, generalised Wishart priors are actually quite inflexible. They are not capable of capturing a DKP prior because in a generalised Wishart process, the Wishart matrices are generated from underlying features, and these features are jointly multivariate Gaussian at all locations (Sec. 4 in Wilson & Ghahramani, 2010) and therefore lack the required nonlinearities between layers. In addition, inference is also very different. In particular, inference for the generalised Wishart is generally performed on the underlying multivariate Gaussian feature vectors (Eq. 1 e.g. Eq. 15-18 in Wilson & Ghahramani 2010, Eq. 12 in Heaukulani & van der Wilk 2019 or Eq. 24 in Jorgensen et al. 2020). Unfortunately, variational approximate posteriors defined over multivariate Gaussian feature vectors fail to capture symmetries in the true posterior (Eq. 14). In contrast, we define approximate posteriors directly over the symmetric positive semi-definite Gram matrices themselves, which required us to develop new, more flexible distributions over these matrices.

Further linear (inverse) Wishart processes have been used in the financial domain to model how the volatility of asset prices changes over time (Philipov & Glickman, 2006b, a; Asai & McAleer, 2009; Gourieroux & Sufana, 2010; Wilson & Ghahramani, 2010; Heaukulani & van der Wilk, 2019). Importantly, inference in these dynamical (inverse) Wishart processes is often performed by assuming fixed, integer degrees of freedom, and working with underlying Gaussian distributed features. This approach allows one to leverage standard GP techniques (e.g. Kandemir & Hamprecht, 2015; Heaukulani & van der Wilk, 2019), but it is not possible to optimize the degrees of freedom and the posterior over these features usually has rotational symmetries (App. D.2) that are not captured by standard variational posteriors. In contrast, we give a novel doubly-stochastic variational inducing point inference method that operates purely on Gram matrices and thus avoids needing to capture these symmetries.

Conclusions

We proposed deep kernel processes which combine nonlinear transformations of the Gram matrix with sampling from matrix-variate distributions such as the inverse Wishart. We showed that DGPs, BNNs (App. C.1), infinite BNNs (App. C.2) and infinite BNNs with bottlenecks (App. C.3) are all instances of DKPs. We defined a new family of deep inverse Wishart processes, and give a novel doubly-stochastic inducing point variational inference scheme that works purely in the space of Gram matrices. DIWP performed better than fully connected NNGPs and DGPs on UCI, MNIST and CIFAR-10 benchmarks.

References

Appendix A DKPs are kernel processes

where $N$ is a positive integer and $\delta$ is a positive real number. A kernel process is defined by consistency under marginalisation and row/column exchangeability. Consistency under marginalisation implies that if we define $\mathbf{K}^{*}$ and $\mathbf{G}^{*}$ as principle submatrices of $\mathbf{K}$ and $\mathbf{G}$, dropping the same rows and columns, then $\mathbf{G}$ being distributed according to a kernel process implies that $\mathbf{G}^{*}$ is distributed according to that same kernel process,

where $\sigma$ is a permutation of the rows/columns. Note that both the Wishart and inverse Wishart as defined in Eq. (35) are consistent under marginalisation and are row/column exchangeable.

A deep kernel process, $\mathcal{D}$, is the composition of two (or more) underlying kernel processes, $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$,

We define $\mathbf{K}^{*}$, $\mathbf{G}^{*}$ and $\mathbf{H}^{*}$ as principle submatrices of $\mathbf{K}$, $\mathbf{G}$ and $\mathbf{H}$ respectively, dropping the same rows and columns, and again, $\mathbf{K}_{\sigma}$, $\mathbf{G}_{\sigma}$ and $\mathbf{H}_{\sigma}$ are those matrices with the rows and columns permuted. To establish that $\mathcal{D}$ is consistent under marginalisation, we use the consistency under marginalisation of $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$

Likewise, to establish row/column exchangeability, we use row/column exchangeability of $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$,

The deep kernel process $\mathcal{D}$ is thus consistent under marginalisation and has exchangeable rows/columns, and hence a deep kernel process is indeed itself a kernel process.

Further, note that we can consider $\mathcal{K}$ to be a deterministic distribution that gives mass to only a single $\mathbf{G}$. In that case, $\mathcal{K}$ can be thought of as a deterministic function which must satisfy a corresponding consistency property,

and this is indeed satisfied by all deterministic transformations of kernels considered here. In practical terms, as long as $\mathbf{G}$ is always a valid kernel, it is sufficient for the elements of $G_{i\neq j}$ to depend only on $K_{ij}$, $K_{ii}$ and $K_{jj}$ and for $G_{ii}$ to depend only on $K_{jj}$, which is satisfied by e.g. the squared exponential kernel (Eq. 4.1) and by the ReLU kernel (Cho & Saul, 2009).

Appendix B The first layer of our deep GP as Bayesian inference over a generalised lengthscale

Under our approach, we sample $\mathbf{F}=\mathbf{F}_{1}$ from Eq. (5b), so $\mathbf{F}$ can be written as,

where $\mathbf{f}_{i}$ is the $i$th row of $\mathbf{F}$. Putting this into a squared exponential kernel without a lengthscale parameter,

We find that a parameter-free squared exponential kernel applied to $\mathbf{F}$ is equivalent to a squared-exponential kernel with generalised lengthscale hyperparameters applied to the input.

Appendix C BNNs as deep kernel processes

Here we show that standard, finite BNNs, infinite BNNs and infinite BNNs with bottlenecks can be understood as deep kernel processes.

Standard, finite BNNs are deep kernel processes, albeit ones which do not admit an analytic expression for the probability density. In particular, the prior for a standard Bayesian neural network (Fig. 3 top) is,

C.2 Infinite NNs

C.3 Infinite NNs with bottlenecks

This generative process is given graphically in Fig. 4 (top).

Appendix D Standard approximate posteriors over features and weights fail to capture symmetries

We have shown that it is possible to represent DGPs and a variety of NNs as deep kernel processes. Here, we argue that standard deep GP approximate posteriors are seriously flawed, and that working with deep kernel processes may alleviate these flaws.

In particular, we show that the true DGP posterior has rotational symmetries and that the true BNN posterior has permutation symmetries that are not captured by standard variational posteriors.

Permutation symmetries in neural network posteriors were known in classical work on Bayesian neural networks (e.g. MacKay, 1992). Here, we spell out the argument in full. Taking $\mathbf{P}$ to be a permutation matrix (i.e. a unitary matrix with $\mathbf{P}\mathbf{P}^{T}=\mathbf{I}$ with one $1$ in every row and column), we have,

i.e. permuting the input to a nonlinearity is equivalent to permuting its output. Expanding two steps of the recursion defined by Eq. (47b),

the output is the same under the original or permuted weights,

Introducing a different perturbation between every pair of layers we get a more general symmetry,

and as the prior over elements within a weight matrix is IID Gaussian (Eq. 47a), the prior probability density is equal under original and permuted weights,

Thus, the joint probability is invariant to permutations,

Due in part to these permutation symmetries, the posterior distribution over weights is extremely complex and multimodal. Importantly, it is not possible to capture these symmetries using standard variational posteriors over weights, such as factorised posteriors, but it is not necessary to capture these symmetries if we work with Gram matrices and kernels, which are invariant to permutations (and other unitary transformations; Eq. 14).

D.2 Rotational symmetries in deep GP posteriors

The same argument can straightforwardly be extended to the inputs, $\operatorname{P}\left(\mathbf{F}_{1}|\mathbf{X}\right)$,

and to the final probability density, for output activations, $\mathbf{F}_{L+1}$ which is not invariant to permutations,

Therefore, applying Bayes theorem the posterior is invariant to rotations,

Importantly, these posterior symmetries are not captured by standard variational posteriors with non-zero means (e.g. Salimbeni & Deisenroth, 2017).

D.3 The true posterior over features in a DGP has zero mean

Appendix E Difficulties with VI in deep Wishart processes

The deep Wishart generative process is well-defined as long as we admit nonsingular Wishart distributions (Uhlig, 1994; Srivastava et al., 2003). The issue comes when we try to form a variational approximate posterior over low-rank positive definite matrices. This is typically the case because the number of datapoints, $P$ is usually far larger than the number of features. In particular, the only convenient distribution over low-rank positive semidefinite matrices is the Wishart itself,

However, a key feature of most variational approximate posteriors is the ability to increase and decrease the variance, independent of other properties such as the mean, and in our case the rank of the matrix. For a Wishart, the mean and variance are given by,

To go about resolving this issue, we need to find a distribution over low-rank matrices with independent control of the mean and variance. The natural approach is to use a non-central Wishart, defined as the outer product of Gaussian-distributed vectors with non-zero means. While this distribution is easy to sample from and does give independent control over the rank, mean and variance, its probability density is prohibitively costly and complex to evaluate (Koev & Edelman, 2006).

Appendix F Singular (inverse) Wishart processes at the input layer

Now, we are able to use the full-rank matrix, $\mathbf{\Omega}$ rather than the low-rank matrix, $\mathbf{G}_{1}$ as the random variable for variational inference. For the approximate posterior over $\mathbf{\Omega}$, in a deep inverse Wishart process, we use

Note in the usual case where there are fewer inducing points than input features, then the matrix $\mathbf{K}_{0}$ will be full-rank, and we can work with $\mathbf{G}_{1}$ as the random variable as usual.

Appendix G Approximate posteriors over output features

To define approximate posteriors over inducing outputs, we are inspired by global inducing point methods (Ober & Aitchison, 2020). In particular, we take the approximate posterior to be the prior, multiplied by a “pseudo-likelihood”,

This is valid both for global inducing inputs and (for small datasets) training inputs, and the key thing to remember is that in either case, for any given input (e.g. an MNIST handwritten 2), there is a desired output (e.g. the class-label “2”), and the top-layer global inducing outputs, $\mathbf{v}_{\lambda}$, express these desired outcomes. Substituting for the prior,

and computing this value gives the approximate posterior in the main text (Eq. 21).

Appendix H Using eigenvalues to compare deep Wishart, deep residual Wishart and inverse Wishart priors

One might be concerned that the deep inverse Wishart processes in which we can easily perform inference are different to the deep Wishart processes corresponding to BNNs (Sec. C.1) and infinite NNs with bottlenecks (App. C.3). To address these concerns, we begin by noting that the (inverse) Wishart priors can be written in terms of samples from the standard (inverse) Wishart

Note that as the standard Wishart and inverse Wishart have uniform distributions over the eigenvectors (Shah et al., 2014), they differ only in the distribution over eigenvalues of $\mathbf{\Omega}$ and $\mathbf{\Omega}^{\prime}$. We plotted the eigenvalue histogram for samples from a Wishart distribution with $N=P=2000$ (Fig. 5 top left). This corresponds to an IID Gaussian prior over weights, with $2000$ features in the input and output layers. Notably, there are many very small eigenvalues, which are undesirable as they eliminate information present in the input. To eliminate these very small eigenvalues, a common approach is to use a ResNet-inspired architecture (which is done even in the deep GP literature, e.g. Salimbeni & Deisenroth, 2017). To understand the eigenvalues in a residual layer, we define a Res$\mathcal{W}$ distribution by taking the outer product of a weight matrix with itself,

where the weight matrix is IID Gaussian, plus the identity matrix, with the identity matrix weighted as $\alpha$,

With $\alpha=1$, there are still many very small eigenvalues, but these disappear as $\alpha$ increases. We compared these distributions to inverse Wishart distributions (Fig. 5 bottom) with varying degrees of freedom. For all degrees of freedom, we found that inverse Wishart distributions do not produce very small eigenvalues, which would eliminate information. As such, these eigenvalue distributions resemble those for Res$\mathcal{W}$ with $\alpha$ larger than $1$.

Appendix I Doubly stochastic variational inference in deep inverse Wishart processes

Due to the doubly stochastic results in Sec. 5.3, we only need to compute the conditional distribution over a single test/train point (we do not need the joint distribution over a number of test points). As such, we can decompose $\mathbf{G}$ and $\mathbf{\Psi}$ as,

Appendix J Samples from the 1D prior and approximate posterior

First, we drew samples from a one-layer (top) and two-layer (bottom) deep inverse Wishart process, with a squared-exponential kernel (Fig. 6). We found considerable differences in the function family corresponding to different prior samples of the top-layer Gram matrix, $\mathbf{G}_{L}$ (panels). While differences across function classes in a one-layer IW process can be understood as equivalent to doing inference over a prior on the lengthscale, this is not true of the two-layer process, and to emphasise this, the panels for two-layer samples all have the same first layer sample (equivalent to choosing a lengthscale), but different samples from the Gram matrix at the second layer. The two-layer deep IW process panels use the same, fixed input layer, so variability in the function class arises only from sampling $\mathbf{G}_{2}$.

Next, we exploited kernel flexibilities in IW processes by training a one-layer deep IW model with a fixed kernel bandwidth on data generated from various bandwidths. The first row in Figure 7 shows posterior samples from one-layer deep IW processes trained on different datasets. For each panel, we first sampled five full $\mathbf{G}_{1}$ matrices using Eq.(34a) and (34b). Then for each $\mathbf{G}_{1}$, we use Gaussian conditioning to get a posterior distribution on testing locations and drew one sample from the posterior plotted as a single line. Remarkably, these posterior samples exhibited wiggling behaviours that were consistent with training data even outside the training range, which highlighted the additional kernel flexibility in IW processes. On the other hand, when model bandwidth was fixed, samples from vanilla GPs with fixed bandwidth in the second row displayed almost identical shapes outside the training range across different sets of training data.

Appendix K Why we care about the ELBO

While we have shown that DIWP offers some benefits in predictive performance, it gives much more dramatic improvements in the ELBO. While we might think that predictive performance is the only goal, there are two reasons to believe that the ELBO itself is also an important metric. First, the ELBO is very closely related to PAC-Bayesian generalisation bounds (e.g. Germain et al., 2016). In particular, the bounds are generally written as the average training log-likelihood, plus the KL-divergence between the approximate posterior over parameters and the prior. This mirrors the standard form for the ELBO,

where $x$ is all the data (here, the inputs, $\mathbf{X}$ and outputs, $\mathbf{Y}$), and $z$ are all the latent variables. Remarkably, Germain et al. (e.g. 2016) present a bound on the test-log-likelihood that is exactly the ELBO per data point, up to additive constants. As such, in certain circumstances, optimizing the ELBO is equivalent to optimizing a PAC-Bayes bound on the test-log-likelihood. Similar results are available in Rivasplata et al. (2019). Second, we can write down an alternative form for the ELBO as the model evidence, minus the KL-divergence between the approximate and true posterior,

As such, for a fixed generative model, and hence a fixed value of the model evidence, $\log\operatorname{P}\left(x\right)$, the ELBO measures the closeness of the variational approximate posterior, $\operatorname{Q}\left(z\right)$ and the true posterior, $\operatorname{P}\left(z|x\right)$. As we are trying to perform Bayesian inference, our goal should be to make the approximate posterior as close as possible to the true posterior. If, for instance, we can set $\operatorname{Q}\left(z\right)$ to give better predictive performance, but be further from the true posterior, then that is fine in certain settings, but not when the goal is inference. Obviously, it is desirable for the true and approximate posterior to be as close as possible, which corresponds to larger values of $\mathcal{L}$ (indeed, when the approximate posterior equals the true posterior, the KL-divergence is zero, and $\mathcal{L}=\log\operatorname{P}\left(x\right)$ ).

Appendix L Differences with Shah et al. (2014)

For a one-layer deep inverse Wishart process, using our definition in Eq. (18)

Importantly, we do the nonlinear kernel transformation after sampling the inverse Wishart, so the inverse-Wishart sample acts as a generalised lengthscale hyperparameter (App. B), and hence dramatically changes the function family.

In contrast, for Shah et al. (2014), the nonlinear kernel is computed before, the inverse Wishart is sampled, and the inverse Wishart sample is used directly as the covariance for the Gaussian,

This difference in ordering, and in particular, the lack of a nonlinear kernel transformation between the inverse-Wishart and the output is why Shah et al. (2014) were able to find trivial results in their model (that it is equivalent to multiplying the covariance by a random scale).