An Analysis of Posterior Collapse, Parameterization and Initialization in Variational Deep Gaussian Processes Francisco Javier Sáez-Maldonado∗
Department of Computer Science and Artificial Intelligence Universidad de Granada Granada, Spain
arXiv:2606.25882v1 [cs.LG] 24 Jun 2026
Juan Maroñas∗
Department of Quantitative Methods CUNEF Universidad Madrid, Spain
Daniel Hernández-Lobato
Department of Computer Science and Centro de Investigación Avanzada en Física Fundamental Universidad Autónoma de Madrid Madrid, Spain.
Abstract Deep Gaussian Processes (DGPs) are probabilistic models with remarkable prediction performance that concatenate Gaussian Processes (GPs) across several layers. Exact inference in DGPs is, however, intractable, and variational inference (VI) is often used to approximate the posterior with a parametric distribution tuned by minimizing the KullbackLeibler divergence. Moreover, finding a good VI approximation is challenging due to the model’s strong nonlinearities and multiple posterior modes. In particular, a problem of VI is posterior collapse, where VI converges to a variational posterior that matches the prior. In variational DGPs, this implies explaining the data as noise. In this work, we study posterior collapse in DGPs and identify its connection to the Double Stochastic VI algorithm and the widely used linear prior mean function employed in all but the last layer of a DGP. We show that the benefit of the linear prior mean does not arise from avoiding the non-injective pathology in very deep DGPs, as previously believed, but from improving the conditioning of the optimization problem at initialization. Thus, we propose an alternative initialization of a zero prior mean DGP that mimics a DGP with a linear prior mean at initialization. This enables successful training of DGPs without imposing optimization-driven constraints on the prior, allowing to choose the prior based on modeling assumptions rather than optimization convenience. Our analysis considers three common parameterizations of DGPs and shows that not all of them benefit from a linear prior mean. We also explain why a whitened parameterization of the DGP provides more stable convergence, something often assumed from experience, but lacking a rigorous analysis. Furthermore, we show that this stability is also beneficial to avoid the posterior collapse problem. Extensive experiments on synthetic and UCI datasets validate our findings. Namely, the proposed initialization prevents posterior collapse, improves stability, and achieves performance comparable to (and sometimes better than) DGPs with a linear prior mean.
∗. Equal Contribution. Juan Maroñas is also with the Machine Learning Group at Universidad Autónoma de Madrid.
©2026 Francisco Javier Sáez-Maldonado, Juan Maroñas and Daniel Hernández-Lobato. License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/.
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Keywords: Deep Gaussian Processes, Variational Inference, Posterior Collapse, Initialization, Mean Function, Reparameterization, Whitening, Optimization.
Contents 1 Introduction
3
2 Related Work
7
3 Deep Gaussian Processes 11 3.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Parameterizing DGPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.1 Non-Whitened DGPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.2 Whitened DGPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3.3 Non-Whitened DGPs with a Reparameterized Variational Mean . . . . 16 3.4 Prior Mean Functions for DGPs . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Model Initialization and Specification . . . . . . . . . . . . . . . . . . . . . . . 19 4 Posterior Collapse in the Context of DGPs 21 4.1 Impact of Inner Layer Prior Mean Functions . . . . . . . . . . . . . . . . . . . 22 4.2 An Analysis of the Posterior Collapse Problem on a Toy Problem . . . . . . . 23 4.2.1 ELL Term and Initial Predictive Distribution . . . . . . . . . . . . . . 24 4.2.2 Explaining the Differences in the Initial Predictive Variance . . . . . . 26 4.2.3 Virtual Dataset Seen by a ZERO Prior Mean DGP at Initialization . . 28 4.2.4 Impact of Initial Predictive Variances in the ELL . . . . . . . . . . . . 28 4.2.5 Whitened Parameterization . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2.6 Initializing S = Sv = I at the Output Layer . . . . . . . . . . . . . . . 30 4.2.7 Impact of the Kernel Hyper-parameters . . . . . . . . . . . . . . . . . 34 4.2.8 The Effect of More Inducing Points . . . . . . . . . . . . . . . . . . . . 35 4.2.9 The Standard Non-whitened Parameterization . . . . . . . . . . . . . . 35 4.2.10 Mixing Matrices and Warping Functions . . . . . . . . . . . . . . . . . 38 4.3 An Approximate Explanation for the Posterior Collapse Problem . . . . . . . 38 4.3.1 Approximating at Initialization a DGP Using a SVGP . . . . . . . . . 38 4.3.2 Simulating the Coordinate Updates . . . . . . . . . . . . . . . . . . . . 41 4.4 Impact of the Parameterization on the Optimization Process . . . . . . . . . . 45 4.4.1 The KLD at Initialization . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.4.2 Optimization Instabilities in the KLD Term . . . . . . . . . . . . . . . 47 4.5 Multi-output Models, Other Likelihoods and Alternative Assumptions . . . . 48 4.6 Conclusions about how to Initialize DGPs . . . . . . . . . . . . . . . . . . . . 49 5 An Initialization to Alleviate Posterior Collapse 49 5.1 The Predictive Mean at Initialization . . . . . . . . . . . . . . . . . . . . . . . 50 5.2 Solving the System of Equations . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.3 Arbitrary Prior Mean Functions . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.4 Fixing Xxz and Initializing the Output Layer To Predict the Targets . . . . . 52 2
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
5.5 5.6 5.7 5.8
Initialization of the Inducing Points . . . . . . . . . . . . . . . . . . . . . . . . Selecting the Initial Kernel Length-scale ℓ . . . . . . . . . . . . . . . . . . . . Summary of the Proposed Initialization Strategy . . . . . . . . . . . . . . . . Dependent and Non-Gaussian Processes at Each Layer . . . . . . . . . . . . .
53 54 56 56
6 Experiments 57 6.1 Toy Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.1.1 An Illustration of the Posterior Collapse Problems of DSVI . . . . . . 57 6.1.2 Evaluating a Very Deep DGP Model . . . . . . . . . . . . . . . . . . . 62 6.1.3 Results for Each DGP Model and Initialization Considered . . . . . . . 62 6.2 Experiments on the UCI Datasets . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2.1 Main Results For Both Parameterizations . . . . . . . . . . . . . . . . 66 6.2.2 Analyzing the Posterior Collapse Problem on the UCI Datasets . . . . 69 7 Conclusions and Future Work
73
A Coordinate Updates for the Noise Parameter
81
B Additional Details of the Proposed Initialization 83 B.1 Equations and Parameters from the System of Equations . . . . . . . . . . . . 83 B.2 Fixing the Kernel Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 B.3 Option 1: 2 dim(Z) Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 B.4 Option 2: dim(Z) Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 B.5 Selection of the Initial Length-scale . . . . . . . . . . . . . . . . . . . . . . . . 85 C Further Experimental Results 87 C.1 Additional Figures for Section 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 87 C.2 Coordinate Updates for 64 Inducing Points . . . . . . . . . . . . . . . . . . . 87 C.3 Additional Toy Experiments Illustrating other Parameter Initializations . . . 94 C.4 Additional Toy Experiments Illustrating Optimization Instabilities Depending on the Parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 C.5 Additional UCI Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 C.5.1 RMSE Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 C.5.2 Comparing Whiten and Non-Whiten Models . . . . . . . . . . . . . . . 99 C.5.3 Split-wise Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
1 Introduction Gaussian Processes (GPs) (Rasmussen and Williams, 2005) are flexible non-parametric models that have shown promising results across many applications, such as molecule optimization (Gómez-Bombarelli et al., 2018; Tripp et al., 2023) or uncertainty estimation in Neural Networks (Ortega et al., 2024; Liu et al., 2023). They are used as prior distributions over some target function f , i.e. f ∼ GP ν (µ(·), K(·, ·)), where µ(·) and K(·, ·) are, respectively, the mean and covariance functions of the GP, parameterized by the set of parameters ν. Deep Gaussian Processes (DGPs) (Damianou and Lawrence, 2013) extend GPs by concatenating multiple GPs in a layered network, thereby defining a hierarchical structure that leads to 3
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
a more flexible probabilistic model. This concatenation enlarges the class of functions that can be modeled, allowing to capture complex patterns within the data. Fig. 1 (a) shows a graphical illustration of a DGP. Typically, the inner layers have as many units as there are input dimensions (Salimbeni and Deisenroth, 2017). In spite of the aforementioned advantages of DGPs, these models have the drawback that exact Bayesian inference is intractable. Thus, the community has been exploring different ways of computing posterior approximations, most of them relying on sparse GP representations based on inducing points and variational methods (Damianou and Lawrence, 2013; Salimbeni and Deisenroth, 2017; Yu et al., 2019, 2021), on Expectation Propagation (Bui et al., 2016), on Markov Chain Monte Carlo (MCMC) (Havasi et al., 2018a), and, recently, on Denoising Diffusion Variational Inference (Xu et al., 2024). One of the most popular variational inference approximations is the Doubly Stochastic Variational Inference (DSVI) algorithm (Salimbeni and Deisenroth, 2017), due to its implementation simplicity, performance, scalability, and its widespread adoption into various DGP inference frameworks and alternative formulations of DGP priors (Havasi et al., 2018a; Jankowiak and Gardner, 2019; Hamelijnck et al., 2019; Sáez-Maldonado et al., 2023; Salimbeni et al., 2019; Rudner et al., 2020). Two design choices have gained popularity over the last years when training DGP models using DSVI. The first one involves using a whitened representation of the process values at the inducing points (Hensman et al., 2015). The second one consists of using a nonlearnable linear-mean function in the inner layers of a DGP, firstly proposed by Salimbeni and Deisenroth (2017). The motivation for the first design choice is that it is well known that parameterizations that remove correlations enhance mixing when fitting probabilistic models using MCMC (Papaspiliopoulos et al., 2007; Betancourt and Girolami, 2015). The use of this parameterization is nowadays the de facto way of fitting these models when using Variational Inference (VI), as implemented by many popular GP/DGP libraries (Matthews et al., 2017; Gardner et al., 2018; Dutordoir et al., 2021). Beyond libraries, the community has also adopted the whitened representation in different research directions involving GPs or DGPs (Salimbeni et al., 2019; Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023). However, to the best of our knowledge, the reason why this parameterization is beneficial for VI algorithms has not been adequately studied. In fact, these modern libraries do not explicitly argue why using a whitened representation is a better choice, beyond being recommended in its source code as seen, e.g., in the library documentation by Gardner et al. (2018)1 . Some works also outline that whitening can facilitate optimization (Leibfried et al., 2022), but provide no theoretical nor empirical evidence for it. The second design choice was motivated to avoid the non-injective pathologies of very deep DGPs with a ZERO prior mean function (Duvenaud et al., 2014). To alleviate these pathologies, Salimbeni and Deisenroth (2017) introduced a prior mean function that is simply equal to the identity function, or a PCA projection of the inputs, when the number of units in each hidden layer is smaller than the number of input dimensions. For these reasons, it is also known as the PCA prior mean function in the literature. Notwithstanding, such a mean function is used even in 2-layer DGPs, which are far from being close to the models 1. https://github.com/cornellius-gp/gpytorch/blob/44993efcc180bdbdeaaf2107c7cc1ba532b2da9b/ gpytorch/variational/unwhitened_variational_strategy.py#L30
4
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
studied by Duvenaud et al. (2014), with up to 11 layers. Moreover, Damianou and Lawrence (2013) successfully trained 5 layers ZERO prior mean DGPs2 which suggests that the original motivation to change the prior mean function of the inner layers is not completely justified. Similar to the whitening reparameterization, the PCA prior mean function has become the default choice in the inner layers of a DGP in modern software (Dutordoir et al., 2021; Gardner et al., 2018) and research (Salimbeni et al., 2019; Rudner et al., 2020; Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023; Sáez-Maldonado et al., 2023), with all attributing its success to the non-injective pathological problem of DGPs (Rudner et al., 2020; Salimbeni et al., 2019). In this work, we found that the success of the PCA prior mean function is linked to optimization issues that arise from the selection of the initial variational parameters within the DSVI inference framework. In particular, issues that arise when initializing the variational mean to m = 0, a common choice in the literature (Salimbeni and Deisenroth, 2017; Salimbeni et al., 2019; Matthews et al., 2017; Dutordoir et al., 2021; Gardner et al., 2018). Specifically, we observed that the DSVI algorithm applied to a ZERO prior mean DGP with this initialization can result in a bad local optimum. The variational learning objective collapses the approximate posterior to the prior, and the model learns to explain the data as noise. This is a well-known effect studied in Variational Autoencoders equipped with an expressive decoder, known as posterior collapse (He et al., 2019; Chen et al., 2017; Bowman et al., 2016). However, the sources of the collapse problem differ from those we found in the context of DGPs. One may find this behavior unsurprising since previous works have shown that sparse GPs fitted via VI tend to learn observation uncertainty through the noise parameter (Bauer et al., 2016). Moreover, some works suggest training sparse GPs by first freezing model hyper-parameters (Leibfried et al., 2022; Hensman et al., 2013), or just the noise parameter (Hamelijnck et al., 2019), before fitting everything end-to-end, including variational parameters. Here, we show that the posterior collapse problem in DGPs is alleviated by using the PCA prior mean function in the inner layers of the DGP model. Since the architectural choice of the prior should be dictated by modeling considerations, and not by limitations due to an optimization issue (Knoblauch et al., 2022), we propose a straightforward initialization of the variational parameters that allows successful training of DGPs using the DSVI inference algorithm with independence of the specified prior mean function. For this, we propose to initialize the parameters of the model so that the output of each GP within the inner layers of the DGP mimics the output of a linear-mean prior GP at initialization. This can be efficiently achieved by initializing the variational means so that the inner layer predictions just at the inducing locations are similar to the output of the PCA prior mean function, i.e., the identity function, when the number of hidden units is equal to the number of input dimensions, or a PCA projection of the inputs, when this number is smaller than the input dimensions. Furthermore, we extend this idea to initialize the last layer of the DGP so that the initial predictive mean of the DGP explains the observed targets well, which is expected to facilitate fitting the DGP model. Our analysis considers both a whitened and a non-whitened parameterization of a DGP. We show that posterior collapse can occur under both parameterizations. Thus, the parameterization itself is not the root cause of collapse. However, the whitened pa2. Neil Lawrence’s blog also showed successful optimization of 5-layers zero mean DGPs https:// inverseprobability.com/talks/notes/deep-gaussian-processes.html.
5
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
rameterization consistently provides much more stable optimization, because whitening removes prior-induced correlations and yields better-conditioned gradients. By contrast, the non-whitened parameterization is significantly more unstable, exhibiting noisy updates, large fluctuations in the KL divergence, and a higher chance of converging to poor local optima. Regarding the non-whitened reparameterization, we also investigate the impact of a reparameterization of the variational mean (i.e., the variational mean is simply the variational mean minus the prior mean). This reparameterization is used, in libraries such as GPFLOW, and is originally motivated in the context of inter-domain GPs (Leibfried et al., 2022). Our initial analysis shows that such a non-whitened reparameterization benefits further from the PCA prior mean function than the standard non-whitened parameterization implemented, e.g., in GPYTORCH or GPJAX. However, both non-whitened parameterizations suffer from optimization instabilities, even if the PCA prior mean function is used, which may result in poorly fitted solutions that suffer from the posterior collapse problem. An initial study of optimization instabilities, in terms of Jacobians and Hessians of the Kullback-Leibler divergence regularizer of VI, is carried out in a separate technical report (Maroñas, 2025). The experimental evaluation carried out includes a series of studies on both a controlled synthetic dataset and eight real-world UCI benchmarks, designed to validate the theoretical analysis of posterior collapse and the proposed initialization strategy. These experiments clearly show that ZERO prior mean DGPs often collapse under standard initialization, while PCA prior mean DGPs avoid collapse only because of their favorable initialization behavior. Furthermore, the results obtained show that the proposed initialization consistently stabilizes optimization, avoids collapse in cases where standard DGPs fail, and achieves predictive performance comparable to (and sometimes surpassing) PCA prior mean DGPs models, confirming its effectiveness in practical scenarios. Summing up, this paper extends our unpublished preliminary results (Sáez-Maldonado et al., 2024) through the following contributions: • We provide a comprehensive analysis of posterior collapse in DGPs trained with DSVI, showing that collapse arises from the interaction between standard variational parameter initialization and the inner-layer mean function, rather than from the non-injective pathologies traditionally cited in the literature. • We reinterpret the role of the widespread PCA prior mean function, demonstrating that its effectiveness stems from how it conditions the optimization landscape at initialization, rather than from improved modeling properties. • We propose a simple and computationally efficient initialization scheme for ZERO prior mean DGPs that mimics the beneficial initialization behavior of PCA prior mean DGPs. This enables successful training of DGPs under principled ZERO mean priors, avoiding optimization-driven architectural constraints. • We analyze three commonly used parameterizations of DGPs. Namely, the standard non-whitened, the whitened, and the GPFLOW reparameterized non-whitened forms, and show that they exhibit fundamentally different optimization behaviors, with only some benefiting from PCA prior mean functions. 6
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
• We provide the first detailed explanation of why the whitened parameterization yields more stable convergence, clarifying an assumption that has been widely accepted in practice but lacked rigorous justification. • We identify conditions under which even whitened DGPs with PCA prior mean functions can suffer posterior collapse, particularly when kernel output scale or inducing-point configurations introduce excessive variance in the inner layers. • Through extensive experiments on synthetic and UCI benchmark datasets, we show that the proposed initialization prevents posterior collapse, stabilizes optimization, and achieves predictive performance comparable to or better than that of DGPs with PCA prior mean functions. The rest of the paper is as follows: Sec. 2 describes important related work from the literature that is relevant for the current investigation. Sec. 3 describes the DGP models, the role of prior mean functions, and the typical initializations and parameterizations that are used for these models. Sec. 4 analyzes the problem of posterior collapse in the context of DGP and investigates when it is more likely to happen, for each prior mean function considered and initialization. Sec. 5 proposes novel initializations of a ZERO prior mean DGP to alleviate the posterior collapse problem. Sec. 6 includes synthetic and real-world experiments to validate our findings and the proposed initializations. Finally, Sec. 7 gives the conclusions of the paper.
2 Related Work This section reviews important work from the literature that is relevant to the current investigation. First, inference sub-optimality in variational methods has been widely studied and is associated with a poorly expressive approximate posterior distribution, which is only able to fit some modes from the true posterior (Bishop, 2007). This effect is known as mode collapse. The problem arises from the way in which the Kullback-Leibler Divergence (KLD) is minimized, which forces the variational distribution to avoid regions of low posterior density. Consequently, another issue derived from this minimization problem is that the variational posterior tends to underestimate the posterior variance. In modern machine learning, the mode collapse problem has mainly been studied in the context of deep generative models based on variational inference (VI), such as the variational autoencoder (Kingma and Welling, 2014). For a long time, the community has focused on creating expressive posterior approximations that can circumvent the limitations of factorized Gaussian distributions (Rezende and Mohamed, 2015), showing improved quality and diversity of generated samples. While not explicitly stated as a solution to mode collapse, improving upon a unimodal Gaussian distribution directly targets this issue. That is, better approximations to the posterior necessarily lead either to unimodal non-Gaussian densities or to multimodal distributions. However, other research directions have also attributed this collapse problem to the decoder (Cremer et al., 2018). Beyond deep generative models, the mode collapse problem is still studied in probabilistic graphical models such as the Gaussian mixture model. For instance, Soletskyi et al. (2025) provides a theoretical explanation of this problem using the gradient flow associated with the optimization procedure. 7
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
On the other hand, in the deep generative models literature, one often finds the term mode collapse used to define the issue in which a generative model is only able to generate samples from a subset of the data distribution, a problem initially observed in generative adversarial networks (Goodfellow et al., 2014). However, this term has also been adopted and studied in many different generative models, such as variational autoencoders (Javaloy et al., 2022), or diffusion models (Barceló et al., 2024), where the collapse problem is attributed to the presence of conflicting gradients during training. Another issue associated with a collapse problem is what is known as posterior collapse. Again adopted by the deep generative community, this term refers to the approximate posterior collapsing to the prior in generative models trained via VI. This effect is favored because that is precisely what minimizes the KLD term in the objective. Some works attribute this issue to variational autoencoders equipped with expressive decoders (Chen et al., 2017; He et al., 2019; Bowman et al., 2016) showing that optimization results in a model working as a standard autoregressive model. This is because the decoder model is so expressive that it is able to map the noise distribution parameterized by the prior and transform it into a sample from the data distribution, which is the fitting mechanism used in many of such models (van den Oord et al., 2016a,b). This results in a model that does not need to encode meaningful information in the inference network (the encoder) to reconstruct the data samples, which makes the objective function set the approximate posterior equal to the prior. Other works show that posterior collapse arises due to spurious local maxima in the training objective (Lucas et al., 2019) or attribute it to a non-identifiability problem (Wang et al., 2021). In our particular case, the issues we have found in the context of DGPs can be associated with the problem studied by Lucas et al. (2019), since under certain initialization, there are solutions in the training objective of the DGP where the model learns everything as noise and collapses the variational distribution to the prior. Both the mode collapse and posterior collapse problems naturally appear in Bayesian neural networks fitted via VI. In the modern era of machine learning, Graves (2011) and Blundell et al. (2015) studied unimodal posterior approximations to the true posterior, which were then improved in follow-up works to yield more expressive approximate posteriors, directly targeting the mode collapse problem (Louizos and Welling, 2016, 2017). On the other hand, Wu et al. (2019) studied a hierarchical prior with parameters fitted via empirical Bayes. The authors state that by allowing too many degrees of freedom in the prior distribution, the empirical Bayes update would always lead to a KLD of zero because the prior parameters would be set to the variational parameters. This results in a posterior collapse problem derived from an architectural design and optimization procedure. Our work studies posterior collapse in the context of DGPs where the problem arises from a different optimization issue. In the context of sparse Gaussian processes, the standard regression likelihood (homoscedastic Gaussian) results in a posterior distribution that is Gaussian. Thus, variational sparse GPs were originally designed to scale computations (Titsias, 2009). Here, the choice of a Gaussian variational distribution was justified, and no collapse problem is expected. This idea opened the possibility of using variational distributions to fit sparse GPs with non-Gaussian likelihoods such as classification (Hensman et al., 2015) or heteroscedastic observation models (Lázaro-Gredilla and Titsias, 2011). Naturally, this approximation suffers from both the collapse and variance underestimation problems. In consequence, inference 8
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
sub-optimality has also been addressed using MCMC (Hensman et al., 2015) to accommodate non-Gaussian posteriors arising from the choice of the likelihood. In DGPs, the posterior distribution is non-Gaussian. The initial variational approach by Damianou and Lawrence (2013) considers a joint variational Gaussian distribution between all the latent functions in the hierarchy. Later, Salimbeni and Deisenroth (2017) proposed a variational joint posterior where the distributions are conditionally Gaussian, given a sample from the previous layer, since they observed that the original inference distribution in Damianou and Lawrence (2013) tends to underestimate the posterior variance. This is because the approximate posterior is not able to model correlations between the latent functions within layers. While the joint distribution over the latent functions suggested by Salimbeni and Deisenroth (2017) is non-Gaussian, and thus aims at directly tackling more complex posteriors, Havasi et al. (2018a) shows further that the posterior over just the inducing points is both non-Gaussian and multimodal. Therefore, those authors propose using MCMC to draw samples from it. Following this observation Yu et al. (2019, 2021) draw inspiration from the deep generative community and propose improving the variational distribution to make it non-Gaussian and multimodal using normalizing flows. Furthermore, there are recent alternatives based on approximating the posterior distribution through denoising diffusion models(Xu et al., 2024, 2026). Although none of these works explicitly mention the problem of mode collapse, they, similarly to deep generative models, target both the mode collapse and the variance underestimation problem that characterizes variational inference by considering more flexible posteriors. Notwithstanding, none of these works explicitly studied or reported inference suboptimality due to a posterior collapse, as we do. In the context of DGPs, there is also a line of research work studying variance collapse in the inner layers of the model. In particular, Ustyuzhaninov et al. (2020) found that, in noise-free settings, DGPs fitted via variational inference with factorized distributions result in the variational distributions of the inner layers collapsing to a deterministic function. In principle, with respect to any parameterization of the model. Later, Popescu et al. (2022) performed a similar study and showed that when the number of inducing points increases, the variance collapses under a DGP with a ZERO prior mean function in the inner layer. They also show that the PCA prior mean function solves this problem. While this is probably the closest line of work related to ours, we do not attribute the posterior collapse nor optimization issues found to the variance tending to zero, but to the actual signal seen by the objective function to be optimized at initialization. Regarding the parameterization of the process values at the inducing points, the use of a whitened representation of random variables goes beyond the field of optimization in the context of DGPs. For instance, in latent Gaussian models, Murray and Adams (2010) discussed the use of a whitened representation of the random variables to make those variables independent under the prior, simplifying posterior sampling using MCMC. This is also used in single-layer sparse GPs (Hensman et al., 2015). In general, a reparameterization of probabilistic graphical models is often used to remove correlations in the posterior that enhance mixing (Papaspiliopoulos et al., 2007; Betancourt and Girolami, 2015). Moreover, in the context of GPs, using a whitened representation of the process is commonly assumed. For instance, Zhu et al. (2023) decouples the kernel hyper-parameters in the mean and covariance conditionals of the SVGP. Then, these authors highlight that, empirically, they found it useful to use whitening in the mean kernel matrix. Also, Shi et al. (2020) proposed 9
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
to decompose the SVGP as a sum of two independent processes, pointing out the importance of using a whitened representation to make optimization easier by reducing the correlation in the posterior distributions. The justifications provided in these previous works for using a whitening representation are, however, scarce and rely on experimental arguments, as also indicated in GP libraries (Gardner et al., 2018), GP tutorials (Leibfried et al., 2022) or other research works (Salimbeni et al., 2019). Here, we provide explanations on why using the whitening reparameterization helps the optimization process, finding the prior covariance as the main cause for this. Notwithstanding, note that the non-whitened parameterization has also been used as the default parameterization in other DGP research works (Salimbeni and Deisenroth, 2017; Rudner et al., 2020). Our work also allows us to conclude that, when using a non-whitening parameterization, a reparameterization performed of the variational mean (the prior mean is subtracted from the variational mean), as done in GPFLOW, does not provide significant additional convergence stability. However, it contributes to the effectiveness of the PCA prior mean function in avoiding the posterior collapse problem. In the literature, this extra reparameterization is simply introduced as a way to unlock the versatility of interdomain GP frameworks (Leibfried et al., 2022). Furthermore, its adoption into DGPs or GP libraries (Salimbeni and Deisenroth, 2017; Dutordoir et al., 2021; Salimbeni et al., 2019; Matthews et al., 2017) has not been properly justified, and is probably just due to a library implementation simplicity. We would like to remark that, to our knowledge, none of these previous works study the relation between optimization results and the different parameterizations considered in the context of DGPs. Regarding the architectural designs of DGPs, Havasi et al. (2018b) proposes a decoupled parameterization of the GP based on the duality between GPs and Gaussian measures, which was first considered by Cheng and Boots (2017). To optimize the proposed model, the authors use the PCA prior mean function in the inner layers, arguing that it avoids degenerate solutions for the covariance matrix and that it also helps with the initialization of the inducing points. Furthermore, the proposed parameterization is compared to the standard SVGP and to the whitened reparameterization of the variational mean parameter, mentioning that both of them are more stable than the one from Cheng and Boots (2017). Also, in deep inter-domain GPs, Rudner et al. (2020) state that the use of a PCA prior mean function in the inner layers is required to avoid the pathologies presented by Duvenaud et al. (2014), similar to Salimbeni and Deisenroth (2017) and Salimbeni et al. (2019). In the three papers, even though the decision of using a PCA prior mean function is supported by practical evidence, no strong theoretical or empirical arguments are provided to justify that decision. Additionally, De Ath et al. (2020) provides a purely experimental study on the role of the prior mean function in GPs, limiting the scope of the study to Bayesian optimization. By contrast, our work analyzes the role of the prior mean function in DGPs trained with the DSVI algorithm from an optimization viewpoint, something that, to our knowledge, has not been done yet. Finally, the initialization of the SVGP parameters is studied in Zhu et al. (2022), where a novel method for initializing the inducing points, based on kernel-based least squares, is proposed. While this initialization helps to improve the performance of the model, it also requires solving an optimization problem iteratively before training the model. Furthermore, it is only applied to standard GPs and the study does not consider DGPs. In our work, by 10
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
contrast, the proposed initialization requires solving a simpler optimization problem and naturally considers DGPs.
3 Deep Gaussian Processes In this section, we will detail the theoretical background on DGPs needed to follow the rest of this work. Throughout the exposition, we will assume the reader is familiarized enough with sparse variational GPs (Titsias, 2009), their training using stochastic variational inference (Hensman et al., 2013), and their use in the context of the Doubly Stochastic Variational Inference algorithm for DGPs (Salimbeni and Deisenroth, 2017). 3.1 Problem Formulation Consider the observed data D = {(xn , yn ) | xn ∈ X ⊆ RD , yn ∈ Y ⊆ RC , n = 1, . . . , N }, where D and C are the input and output dimensionality, respectively. We write it in matrix form as X = (x1 , · · · , xN )T and Y = (y1 , · · · , yN )T . Our goal is to learn a predictive distribution for y given x and the observed data, i.e., p(y | x, D). We can model this problem by placing a prior distribution over a function that maps x to a distribution over y, and then use the posterior distribution to integrate over all possible functions to make predictions on unseen data instances. We consider a DGP as the prior distribution of this function. DGPs define a prior distribution over vector-valued functions that surpass some of the limitations of standard GP models, such as misspecification from assumed joint Gaussianity between points in the function, or smoothness. The joint prior probability over these functions is obtained by a hierarchical concatenation of GPs. The idea is to organize GPs into layers, as in a fully connected neural network, and feed the GP output of one layer as input to the next one. In its original formulation (Damianou and Lawrence, 2013), within the same layers, the GPs are conditionally independent given a sample realization of the previous layer. Fig. 1 (a) shows a graphical depiction of a 4-layered DGP. Fig. 1 (b) shows the probabilistic graphical model of a L-layered DGP, where each layer has Dl GPs (units). In our problem formulation, fdl is a vector random variable with the values of the d-th latent function in the layer l evaluated at a particular set of points X. We have a GP as the prior of that latent function. Such a GP is characterized by a kernel K(·, ·) with parameters νdl . A particular sample s of such a latent variable, as in a Monte Carlo method, l , and for a single point xn , we might also denote with f l is denoted as fd,s d,n,s the sample of l the corresponding random variable. We assume each layer has D processes, and denote with l ) the set of vector random variables belonging to layer l. Then, we can write f l = (f1l , . . . , fD l the joint distribution between the observed data and the DGP latent variables as: p(Y, X, f 1 , . . . , f L ) = p(Y | f L )
L Y l=1
l
p(f l | f l−1 ) = p(Y | f L )
L Y D Y l=1 d=1
p(fdl | f l−1 ) ,
(1)
with f 0 := X and for simplicity we have removed the dependence at each layer with respect l ). The form of p(Y | f ) will depend on the task considered. For example, to ν l = (ν1l , . . . , νD L for multiclass classification, we can use a categorical likelihood or, for multivariate regression, a multivariate Gaussian distribution with some assumed noise model that may include 11
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Inputs
GP
GP
GP
GP
f (1) (x)
f (2) (x)
f (3) (x)
y
x1
x2
x3 x
(a) Neural Network inspired graphical representation of a 4-layer DGP. D1
x
f1
D2
f2
DL−1
...
f L−1
DL
fL
y
N ν1
ν2
ν L−1
νL
(b) Probabilistic Graphical Model of the Deep Gaussian Process as originally formulated by Damianou and Lawrence (2013).
Figure 1: Graphical representations of a Deep Gaussian Process. dependencies between the regressed dimensions. Each of the factors p(fdl | f l−1 ) describes the prior distribution given by a GP at inputs f l−1 . Furthermore, one can easily model dependencies within the processes in the same layer by aggregating the processes using a mixing matrix (Álvarez et al., 2012; Jankowiak and Gardner, 2019), a warping function (Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023; Sáez-Maldonado et al., 2023), or a combination of both (Lin et al., 2023), and then, at inference, make some assumptions to achieve scalability, such as removing dependencies. This is a practice implemented in many common GP libraries (Matthews et al., 2017), as already noted by Maroñas and Hernández-Lobato (2023). 3.2 Inference Inference in a DGP, i.e., computing p(f 1 , . . . , f L | X, Y), is intractable. One of the most popular ways of performing approximate inference in this model is the Doubly Stochastic Variational Inference (DSVI) algorithm (Salimbeni and Deisenroth, 2017). Since this is our object of study, we present the formulation involved precisely. 12
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
The DSVI formulation first introduces a set of M inducing points Zld and associated process values uld for each latent function in the DGP. Let ul = (ul1 , . . . , ulDl ). The joint distribution of the observed data and all latent variables then becomes: l
L p(Y, X, {f l }L l=1 , {ul }l=1 ) = p(Y | f L )
L Y D Y l=1 d=1
p(fdl | f l−1 , uld )p(uld ) ,
(2)
with p(fdl | f l−1 , uld ) a GP predictive distribution at f l−1 given observed targets uld at Zld , and with p(uld ) a GP prior distribution for uld . Again, for simplicity, we have removed the l ) and with respect to Zl . dependence at each layer with respect to ν l = (ν1l , . . . , νD d L Then, the intractable posterior distribution p({f l }l=1 , {ul }L l=1 | Y, X) is approximated via variational inference, by minimizing the KLD between an approximate family and the true posterior. The novel contribution from Salimbeni and Deisenroth (2017) is to define a variational family that maintains the conditional structure between the layers, as in the prior distribution, while assuming independence within the layers. To yield an efficient algorithm, they also used the uncollapsed version (Hensman et al., 2013) of the variational sparse GP formulation (Titsias, 2009). More precisely, the approximate posterior distribution is l
L q({f l }L l=1 , {ul }l=1 | Y, X) =
L Y D Y l=1 d=1
p(fdl | f l−1 , uld )q(uld ) ,
(3)
with q(uld ) a Gaussian variational distribution with mean and covariances given by mld and Sld , respectively. Thus, the approximate predictive distribution for the d-th GP process values in layer l is:
q(fdl | f l−1 ) =
Z
p(fdl | uld , f l−1 )q(uld )duld ,
(4)
which is Gaussian with mean and covariances given by: −1 (m − µZ ) + µX , µqf = KXZ KZZ
(5)
−1 −1 −1 KZX , SKZZ KZX + KXZ KZZ Kqf = KXX − KXZ KZZ
(6)
where we have removed the dependence with respect to the layer, l, and the hidden dimension, d, to simplify the notation. We denote with µX = µ(f l−1 ) and KXX = K(f l−1 , f l−1 ), the prior mean and corresponding kernel matrix (prior covariances) evaluated at a set of input points f l−1 , respectively. Similarly, we denote with the kernel matrix KZZ = K(Zld , Zld ) the prior covariances associated with the inducing points. The kernel matrix KXZ = K(f l−1 , Zld ) denotes the prior cross-covariances between the process values at the inducing points and the input locations. Note that at the first layer f l−1 = X. This notation will simplify exposing the ideas described later on in our paper. When using the variational posterior in Eq. (3), the learning objective results in the Evidence Lower Bound (ELBO) given by: ELBO = Eq(f ,...,f ) 1 L
log p(Y | f 1 , . . . , f L ) −
L X l=1
KLD [q(ul ) || p(ud )]
S L X Dl h i X 1X l l ≈ Eq(f |f log p(Y | f L ) − KLD q(ud ) || p(ud ) , ) L L−1,s S s=1
l=1 d=1
13
(7)
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
where the first expectation is partially approximated using samples obtained at each level in the hierarchy. These samples are generated recursively by sampling from the Gaussian distribution in Eq. (4). Namely, f l−1,s ∼ q(f l−1 | f l−2,s ). The expectation over the last layer can be computed either via quadrature methods, Monte Carlo, or exact integration, depending on the likelihood. Each of the KLDs between the prior and the approximate posterior in Eq. (7) are given by: KLD =
i −1 1h −1 log | KZZ | − log | S | −M + Tr KZZ S + (m − µZ )T KZZ (m − µZ ) , 2
(8)
where, again, we omit subindices for simplicity. By using likelihoods that factorize across training instances and output dimensions, this training objective in Eq. (7) can be computed using minibatches of data and expectations over univariate variational posterior distributions (Salimbeni and Deisenroth, 2017). Moreover, in some regression problems, the last expectation in Eq. (7) admits a closed-form solution. Gradients are typically computed efficiently using the reparameterization trick (Kingma and Welling, 2014) and automatic differentiation. When the likelihood does not factorize across target dimensions, expectations can still be computed by sampling from a multivariate Gaussian with at most C dimensions per datapoint, as in multiclass classification, with C the number of output dimensions. If a mixing matrix is used at any level of the hierarchy, one often ignores dependencies so that the joint multivariate Gaussian distribution results in independent sampling per GP, alleviating the quadratic memory cost of keeping the full covariance per datapoint (Maroñas and Hernández-Lobato, 2023). Overall, this results in an efficient algorithm that justifies its wide adoption in many settings. When training the model, a single sample per observed instance, i.e., S = 1, is usually enough. A predictive distribution for test instances can be computed in a similar way to that of training instances. In this case, more samples S can be used to increase the accuracy of the Monte Carlo approximation of the predictive distribution. 3.3 Parameterizing DGPs This section discusses how DGP are parameterized in practice. Model parameterization is relevant to multiple parts of the modeling procedure. Here, we focus on how to parameterize the hierarchical prior defined by a DGP. In particular, we will discuss the non-whitened and whitened parameterizations of the prior and approximate posterior over the process values associated with the inducing points. We also describe a reparameterized version of the non-whitened model used by Salimbeni and Deisenroth (2017); Matthews et al. (2017) and Dutordoir et al. (2021). Throughout this subsection, we will omit the layer l and hidden dimension d indices for brevity. 3.3.1 Non-Whitened DGPs In the previous section, we implicitly described the standard non-whitened parameterization. The prior and approximate posterior are given by the following distributions: p(u) = N (u | µZ , KZZ ) , q(u) = N (u | m, S) , 14
(9) (10)
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
which results in a marginal variational posterior distribution with parameters given by Eqs. (5) and (6), and a KLD in the loss function given by Eq. (8). We can find this parameterization in libraries such as GPYTORCH (Gardner et al., 2018) or GPJAX (Pinder and Dodd, 2022). We have also found the inter-domain DGP of Rudner et al. (2020) to use this parameterization. 3.3.2 Whitened DGPs The term whitening is often referred to as a reparameterization of the distributions in a graphical model that removes correlations from the random variables being sampled using MCMC algorithms. This reparameterization allows an efficient exploration of the target distribution (Papaspiliopoulos et al., 2007), something also known in the community as mixing. Intuitively, this reparameterization unstucks the sampling algorithm from regions of the distribution with high geometrical changes. More precisely, it often results in samples that are more representative of the target distribution and less correlated. The funnel distribution is a good example of this phenomenon (Neal, 2003). This reparameterization accounts for transforming the prior process values associated with the inducing points to have 0 mean and unit variance, and, in the context of sparse GPs, was initially introduce to facilitate sampling via MCMC (Hensman et al., 2015)). More precisely, instead of performing inference on u, one performs inference on v, whose prior is defined as a whitened Gaussian distribution. Namely, (11)
p(v) = N (v | 0, I) , with the inducing values u being recovered by the following linear transformation:
(12)
u = LZZ · v + µZ ,
where LZZ is the Cholesky factor of KZZ , i.e., KZZ = LZZ LTZZ . When used in the context of variational inference, the variational distribution is also defined in the whitened space: (13)
q(v) = N (v | mv , Sv ) .
With the reparameterization from Eq. (12), the prior and variational posterior on the unwhitened space takes the form: p(u) = N (u | µZ , KZZ ) , q(u) = N
u | LZZ mv + µZ , LZZ Sv LTZZ
(14)
,
(15)
which results in a predictive Gaussian distribution in Eq. (4) with mean and covariance given by: T µqfw = KXZ L−1 mv + µX , ZZ
(16)
T −1 Kqfw = KXX − KXZ KZZ KZX + KXZ L−1 Sv L−1 ZZ ZZ KZX .
(17)
Note that this reparameterization leaves the GP prior invariant, and so the prior remains the same as in the non-whitened case, meaning that both parameterizations lead to the same 15
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
statistical model. The KLD in the ELBO can be computed from the prior and posterior distributions in the whitened space. That is, KLD[q(v) || p(v)] =
1 − log | Sv | −M + Tr [Sv ] + mTv mv , 2
(18)
since the KLD is invariant under variable transformations. Note that now, the KLD does not depend on the location of the inducing points Z, nor kernel parameters. Only the data-dependent term of VI depends on them. This is expected to facilitate optimization. While originally proposed for MCMC algorithms, this reparameterization has become a de facto standard in many VI algorithms applied to GPs and DGPs (Salimbeni et al., 2019; Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023; Sáez-Maldonado et al., 2023) or directly recommended by GP libraries (Gardner et al., 2018) and sparse GP tutorials (Leibfried et al., 2022). The main reason is the association of this parameterization with better convergence results, but without a rigorous experimental or theoretical study. Our previous unpublished work showed that whitening results in more stable optimization (Sáez-Maldonado et al., 2024). We discuss this later on in this work. 3.3.3 Non-Whitened DGPs with a Reparameterized Variational Mean Another option to implement sparse non-whitened GPs considers a reparameterization of the variational parameters. This is observed in the GPFLUX DGP library (Dutordoir et al., 2021), in the work by Salimbeni and Deisenroth (2017), and also in the implementation of the sparse GPs in the GPFLOW library (Matthews et al., 2017), which is the framework considered in this work. In particular, the GPFLOW non-whitened parameterization reparameterizes the variational mean as: mr = m − µZ .
(19)
Plugging this reparameterization into the variational mean and covariance yields a predictive Gaussian distribution in Eq. (4) with parameters: −1 µqfr = KXZ KZZ mr + µX ,
(20)
−1 −1 −1 Kqf = KXX − KXZ KZZ KZX + KXZ KZZ SKZZ KZX ,
(21)
and the resulting KLD in the ELBO results in: KLD =
−1 1 −1 log | KZZ | − log | S | −M + Tr KZZ S + mTr KZZ mr . 2
(22)
By inspection, one can see that this corresponds to the following prior and variational posterior: p(u) = N (u | 0, KZZ ) , q(u) = N (u | mr , S) .
(23) (24)
There are some concerns about the use of this reparameterization in DGP models. First, the implications behind how this reparameterization modifies the statistical model defined 16
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
by the DGP prior remain underexplored. Note that, now, the prior mean over the values associated with the inducing points is zero, and different from the prior mean associated with any other input. Moreover, we found its use in the context of DGPs not justified by Salimbeni and Deisenroth (2017). Similarly, we did not find a justification for it in the GPFLOW and GPFLUX library’s documentation for both GPs and DGPs (Matthews et al., 2017; Dutordoir et al., 2021). The only reference to the use of this parameterization is found in Leibfried et al. (2022) and is related to inter-domain sparse GP models. It turns out that this reparameterization unlocks the use of any prior mean function in these models, since it removes the convolution integral needed to compute the variational mean posterior, which is often analytically intractable. In spite of the aforementioned connection with inter-domain sparse GP models, we have not found references to the advantages of this reparameterization in the seminal work of inter-domain GPs (Lázaro-Gredilla and Figueiras-Vidal, 2009), nor in modern frameworks for these models (van der Wilk et al., 2020), nor in their practical implementation (Dutordoir et al., 2020). We think this is because single-layer GP models often assume a ZERO prior mean function, resulting in a marginal variational posterior mean common for both non-whitened parameterizations, which does not require the convolution integral. In inter-domain DGPs (Rudner et al., 2020), this reparameterization is not used since the convolutions involving the mean function in the framework do not need to be solved. Overall, this might explain why its adoption in other frameworks such as GPYTORCH or general DGPs models has not been popularized. Thus, in principle, for general GP or DGP models, this reparameterization just implies implementation simplicity, since the routine that evaluates the KLD can be used for both the whitened and non-whitened parameterizations (i.e., by setting KZZ = I). In any case, our analysis shows that this reparameterization is important to precondition the VI algorithm at initialization, such that posterior collapse is avoided. More precisely, we found that the standard non-whitened parameterization with a big number of inducing points and the PCA prior mean function can present a similar behavior to that of DGPs with the ZERO-mean prior function at initialization, locking the potential preconditioning that this mean function introduces. On the other hand, from an optimization viewpoint, this parameterization does not improve stability in the way the whitened parameterization does. Table 1 summarizes the different parameterizations described in this section.
Table 1: A summary of the different DGP prior parameterizations considered in this work. Parameterization
p(u) and q(u)
Non-whitened
p(u) = N (u | µZ , KZZ )
Whitened
GPFLOW Non-whitened
Marginal Variational Posterior −1 µqf = KXZ KZZ (m − µZ ) + µX
−1 −1 −1 Kqf = KXX − KXZ KZZ (KZZ − S)KZZ KZX
q(u) = N (u | m, S)
p(u) = N (u | µZ , KZZ ) q(u) = N
u | LZZ mv + µZ , LZZ Sv LTZZ
p(u) = N (u | 0, KZZ ) q(u) = N (u | mr , S)
T µqfw = KXZ L−1 mv + µX ZZ T Kqfw = KXX − KXZ L−1 [I − Sv ] L−1 ZZ ZZ KZX −1 µqfr = KXZ KZZ mr + µ X −1 −1 −1 Kqf = KXX − KXZ KZZ (KZZ − S)KZZ KZX
17
KLD 1 [log | KZZ | − log | S | −M + 2 −1 −1 Tr KZZ S + (m − µZ )T KZZ (m − µZ )] 1 − log | S | −M + Tr [S] + mT m 2
1 [log | KZZ | − log | S | −M + 2 −1 −1 Tr KZZ S + mTr KZZ mr ]
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
3.4 Prior Mean Functions for DGPs We give details in this section about the evaluation of the prior mean function in the context of DGPs. We specifically describe the PCA prior mean function suggested by Salimbeni and Deisenroth (2017), which is a non-learnable linear mean function. The goal is to provide a solid basis to understand the proposed initialization of a ZERO prior mean DGP described in the following sections, designed to reduce the impact of the posterior collapse problem. We consider the standard DGP architecture, suggested by Salimbeni and Deisenroth (2017), which introduces hidden layers of dimensionality D in the DGP, with D the input data dimension. Therefore, at each hidden layer there are D GPs that compute a non-linear transformation of the input. Consider now a D dimensional input xn which is fed into a hidden DGP layer. If the prior mean is zero, then each GP prior mean outputs a scalar zero, i.e., µ(xn ) = 0. However, if the mean function is a linear mean function parameterized by w, then each GP prior mean outputs the dot product between the parameter and the input, i.e., µ(xn ) = wT xn . This means that, for an arbitrary number of GPs per layer, the whole prior mean value for the layer is evaluated as Wxn , with W a D × D matrix with rows being the parameters of the linear prior mean for each GP unit within the layer. Typically, W is not learnable and has a predefined structure. For example, if the input dimension to a layer and the number of units in the layer are the same, as described above, then, we have W = I, as suggested by Salimbeni and Deisenroth (2017). For this choice, the output of the prior mean of each GP is the corresponding dimension of xn . In other words, if D = 3 then µ1 (xn ) = xn1 , µ2 (xn ) = xn2 and µ3 (xn ) = xn3 , where µd (·) denotes the prior mean function of the d-th GP. In this case, we can think of each mean function being parameterized by a one-hot-encoding vector w with 1 in the d-th dimension and all the other dimensions equal to 0. When the input dimension D is large, it is no longer feasible to have D units in each DGP hidden layer. In this case, a smaller number of units Dl < D is considered at each layer l. Moreover, W is set to a PCA projection matrix from D dimensions to Dl dimensions in the first layer, as suggested by Salimbeni and Deisenroth (2017). The parameters of W, in this case, are found using PCA on the training data. In the other layers, we simply set W = I, as described before. If the number of GPs in a hidden layer is larger than the input dimensionality of that layer D, then W is set equal to the identity matrix with zero-padded column vectors. Thus, the remaining outputs of the layer are set to have a prior mean equal to 0 (Salimbeni and Deisenroth, 2017). Summing up, when the DGP has as many units in each layer as the input dimensionality D, which is the typical scenario, the PCA prior mean function of Salimbeni and Deisenroth (2017) simply sets each hidden layer to model deviations from the identity function that maps each input vector xn to the same input vector xn . Importantly, using this mean function implies that the statistical model defined by the
DGP in the whitened parameterization is different from the non-whitened with a repa-
rameterized variational mean. This contrasts with the zero mean function, where both parameterizations lead to the same statistical model. 18
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
3.5 Model Initialization and Specification In variational DGPs, one usually finds the prior design choices and initialization of the model and variational parameters described in this section. Thus, our analysis will consider them, but also other alternatives. Prior design choices: 1. The covariance functions of the DGP are chosen to be stationary since, by construction, the DGP already models non-stationary processes. Stationary covariance functions depend on the distance between the points. It is common to find square exponential kernels combined with an output scale parameter (Salimbeni and Deisenroth, 2017; Rudner et al., 2020; Salimbeni et al., 2019). 2. It is assumed independence between the latent functions at a layer. This is the original prior formulation by Damianou and Lawrence (2013). However, we also discuss typical mixing procedures (such as linear combinations using a mixing matrix), showing that our analysis and proposed solutions are equivalent, under a correct initialization. 3. We focus our experiments on single-output regression problems, but provide a small discussion in both coupled and independent multi-output Gaussian observation models, and other likelihoods such as those for classification. For regression problems, one usually considers a homoscedastic Gaussian noise model at the output of the DGP. Namely, p(y n |fnL ) = N y n | fnL , σ 2 , with fnL the DGP output associated to xn . The target value is y n , and the noise variance is σ 2 . Thus, for N training instances, we have: p(Y | f L ) =
N Y n=1
N y n | fnL , σ 2 ,
(25)
which results in an expected log-likelihood (ELL) term in the ELBO given by: ELL =
N S X 1X n=1
S
n L Eq(f L |fn,s L−1 ) log p(y | fn ) . n
s=1
(26)
Note that for each training point, we are required to evaluate an expectation approximated by S Monte Carlo samples. Importantly, however, the generated samples are independent for each different point in the training set or mini-batch X (Salimbeni and Deisenroth, 2017). Thus, we only need to generate separate samples from the predictive distribution q(f 1,n , . . . , f L−1,n ) for each individual point xn in X. These independent samples are represented in Fig. 2 using green dotted lines. They are used to compute the ELL term in the ELBO. The red dashed lines represent dependent samples from the predictive distribution at N specific locations. Note that these last samples result in smooth functions sampled from the posterior, unlike the independent samples mentioned before. 19
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Model Initialization: 1. The inducing point locations are initialized following Salimbeni and Deisenroth (2017). At the input layer, each Zld is initialized by running the k-means algorithm over the training points with as many clusters as inducing points. For the layers after the first one, the inducing points are initialized in the following way: • When the layer input dimension is equal to the data dimension, the inducing points are initialized with the same values as those in the input layer. • When the layer input dimension is smaller than the data dimension, a PCA projection of the inducing points is performed. • When the layer input dimension is bigger than the data dimension, the inducing points are copied, and the new dimensions are zero-padded. 2. The noise variance parameter σ 2 is initialized to values typically between 0.01 and 1.0. For example, Salimbeni and Deisenroth (2017) mentions an initial value of 0.01, while their implementation uses an initial value of 1.0 in practice. Hamelijnck et al. (2019) and Sáez-Maldonado et al. (2023) use 0.01 for this parameter, while Maroñas et al. (2021) consider 0.05, and Salimbeni et al. (2019) use 0.1. GPFLOW implements a lower bound of 10−6 in the noise variance, for numerical stability. We use an intermediate noise variance value of 0.5 to illustrate the optimization difficulties of DGPs in the next section. However, in our experiments, we set the noise variance to 0.05. The reason is that one of our conclusions is that smaller values of this parameter at initialization tend to reduce the posterior collapse problem. This is because larger initial noise variances force the model to represent the observed targets just as noise. 3. The variational parameters of each GP within the DGP are often initialized to m = 0, S = 10−5 I for the hidden layers. The output layer is initialized with S = I in the works by Salimbeni and Deisenroth (2017) and Rudner et al. (2020)3 . Other works that report results using the DSVI algorithm consider S = 10−5 I in the output layer (Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023; Sáez-Maldonado et al., 2023). Hamelijnck et al. (2019) mentions a random initialization of the variational mean and variances between 0 and 1. In any case, note that the variational variances should always be initialized to a value smaller than or equal to the prior variances. Moreover, regarding the hidden layers, it seems reasonable to initialize the variational covariance to be small for two reasons. First, we start by being certain about the process values associated with the inducing points, enforcing their adaptation in input space, and we let the model learn the necessary uncertainty. Second, this provides samples from the DGP hierarchy with less variance, improving convergence of the stochastic optimization algorithm. Such an initialization is also recommended in works such as that of Leibfried et al. (2022). We have, however, observed that Salimbeni et al. (2019) uses I to initialize the hidden layer variational covariance. In our work, we will analyze several settings of the initialization of the variational distribution at each layer, following the aforementioned observations and the initial values described. 3. The implementation released by Salimbeni and Deisenroth (2017) initializes S = 10−5 KZZ (hidden layer) and S = KZZ (output layer) but the work mentions S = 10−5 I and S = I. In our analysis, we will show that these two initializations are valid.
20
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
4. The kernel parameters depend upon the specific kernel. We use a Radial Basis Function kernel, given by ! ′ ∥x − x ∥2 ′ K(x, x ) = σo exp − . (27) 2ℓ2 Salimbeni and Deisenroth (2017) initialize ℓ = 2, which results in smooth functions (assuming the data have been standardized). The output scale parameter σo takes a value of 1.0 at initialization and controls the variance of the covariance function, as in the demos given in their source code. In their paper, they use an output scale of 2.0. For this reason, we use an output scale of 2.0 in the UCI benchmark experiments. In our toy problem, we set the output scale to 1.0. We also analyze the influence of this parameter at initialization. As we will justify and analyze later, in some of our models, we propose to adapt the length-scale to be more representative of the data. In particular, the length-scale is selected so that it minimizes the RMSE of the initialized model on the train set. Thus, we follow the common practice of choosing the length-scale adapted to the type of data we are trying to model (Hamelijnck et al., 2019; Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023; Rudner et al., 2020; Leibfried et al., 2022). As with the output scale, we also analyze the influence of the initial length-scale in the next section. The variational covariance S is parameterized via its Cholesky factor, which results in an unconstrained parameterization. Parameters requiring some restriction, such as the variance being positive, are usually parameterized on an unconstrained domain and then mapped through an invertible mapping to the restricted domain (Matthews et al., 2017).
4 Posterior Collapse in the Context of DGPs In this section, we present and analyze the optimization problems that often arise when training DGPs with the DSVI algorithm. Furthermore, we show some modeling choices that help overcome these issues. In particular, we focus our analysis on the choice of the inner layer mean function and on the three parameterizations previously introduced. We also analyze and discuss the use of different numbers of inducing points, or the initial variational and model parameters outlined in Sec. 3.5. On one side, we show that a constant prior mean function, with the typical initialization discussed in the previous section (in particular m = mv = mr = 0), is likely to lead to the posterior collapse problem. That is, the model minimizes the KLD in the ELBO by setting the variational parameters to those of the prior. The ELL term is then maximized by setting the variance of the noise equal to the variance of the targets, assuming they have zero mean. That is, the data are explained as pure noise. Our analysis starts with illustrative and intuitive descriptions of the problem to get a general idea. The problem is then formalized by noticing that a DGP at initialization can be approximated via a SVGP. With this observation, we then use the coordinate updates of the variational parameters of a SVGP (Hensman et al., 2013; Titsias, 2009) and derive the updates for the noise parameter, showing that a constant prior mean function enforces the model to learn a large noise variance while collapsing the variational parameters to the prior. This gives an idea of the type of local minima that the gradient ascent method often finds during the first iterations of optimization. 21
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Furthermore, we show that the posterior collapse problem is present in both parameterizations (whitening and non-whitening), justifying the use of the PCA prior mean function proposed by Salimbeni and Deisenroth (2017) in the non-whitened parameterization and subsequent usage in the whitened parameterization. We also show that, in the non-whitened case, the variational mean reparameterization of Sec. 3.3.3 benefits further from the PCA prior mean function than the standard non-whitened parameterization. Notwithstanding, we observe that posterior collapse in non-whitened models can additionally arise from training instabilities resulting from this parameterization (Sáez-Maldonado et al., 2024). Additionally, we provide insights into the source of these instabilities. Following this observation, we realize that any form of noise in the optimization procedure might result in posterior collapse, regardless of the parameterization. This additional noise might come from a small number of inducing points, a big variational variance S, or a big kernel output scale parameter σo . Based on this observation, the experiment section shows that a whitened DGP with 10 layers and the PCA prior mean function, can also present posterior collapse when the variance coming from the inner layers is big at initialization, and a small number of inducing points is used. This result questions the single original motivation behind using the PCA prior mean function, i.e., simply a large depth of the model (Salimbeni and Deisenroth, 2017). 4.1 Impact of Inner Layer Prior Mean Functions While GPs usually have ZERO prior mean functions, Salimbeni and Deisenroth (2017) suggested using a non-learnable linear mean function in the inner layers of a DGPs, as described in Sec. 3.4. Whenever the input and output dimensionality of the layer is the same (the most common scenario), the mean function copies the input into the output X = µ(X). When the input dimensionality is higher, a PCA projection of the inputs is performed. By contrast, when the output dimensionality is higher, the input is copied to the output and padded with zeros. Attending to Salimbeni and Deisenroth (2017), this design choice is motivated to avoid the pathologies discovered by Duvenaud et al. (2014) in very deep DGPs with ZERO prior mean functions. They mentioned that with this inner mean function, it is effective to initialize the variational mean mld = 0 for any GP. However, the following two observations raise concerns regarding the validity of this statement: 1) Salimbeni and Deisenroth (2017) use this mean function even in 2−layer DGPs, which are far from the deepness of the models (around 11 hidden layers) studied by Duvenaud et al. (2014); and 2) Damianou and Lawrence (2013) reported good predictive performance using ZERO prior mean DGPs with 5 hidden layers, being the inference algorithm the main difference w.r.t. Salimbeni and Deisenroth (2017). Building upon these two observations, we conduct a deeper analysis and illustrate that the poor performance of a DGP with a constant prior mean function and the number of layers considered by Salimbeni and Deisenroth (2017) is a consequence of posterior collapse, which originates from the inference algorithm with the initialization outlined in Sec. 3.5. In particular, the choice m = mr = mv = 0. This behavior is therefore not attributable to the non-injective mapping problem discussed by Duvenaud et al. (2014). We show that setting the variational mean equal to zero conditions the output layer of a ZERO prior mean DGP model to see, just after initialization, a virtual dataset consisting of pairs of {0, yn }N n=1 , for any training point xn . That is, the output of the last inner layer, which is the input to the 22
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
last DGP layer, tends to predict only zeros (just after initialization) for any input point xn . By contrast, a PCA prior mean DGP does not have this problem. At initialization, the last inner layer will tend to predict xn given xn as an input. The consequence is that the PCA prior mean function, when m = mr = mv = 0, conditions the optimization algorithm to run away from a minimum that occurs when the ELL term in the ELBO is minimized by learning a huge observation noise. The reason is that, on expectation, it makes the model output a predictive variance that does not agree with the observed data. This forces the model to adjust the variational parameters to compensate for this behavior. By contrast, a constant ZERO prior mean function with the aforementioned initialization outputs the same predictive variance at any point (recall that each input point xn is mapped to 0). This implies that a homoscedastic noise model will favor learning a huge observation noise, because the last GP model cannot disaggregate information from training points when the input it receives is always 0 for each training point. This could be seen as a non-injective problem attributable to the inference algorithm, but not to the DGP prior, which is the one studied by Duvenaud et al. (2014). Summing up, the origin of the posterior collapse problem is not the initial predictive mean at the training points (Sáez-Maldonado et al., 2024), which is equal to zero at initialization for each prior mean function considered in the inner layers (i.e., ZERO or PCA), but to the initial predictive mean at the inner layers, and the final output predictive variance, which are different at initialization, for each inner-layer prior mean function. 4.2 An Analysis of the Posterior Collapse Problem on a Toy Problem To validate our claims, we first analyze the training signals that are propagated through the DGP hierarchy and used to evaluate the ELL term in the ELBO. To do so, we design a toy experiment which will allow us to analyze how the different mean functions precondition the DSVI algorithm at initialization, under different settings, including the number of inducing points, different initialization parameters, or the different parameterizations described in Sec. 3.3. We depict how the PCA prior mean function consistently shapes the predictive variance in a way such that the ELL can be very low at points where data are present, and the only way to increase it is by modifying the variational parameters. Furthermore, we show that m = mr = mv = 0 is the primary factor contributing to the bad performance of any constant prior mean DGP. This analysis reveals that the shape of the predictive variance is determined by the inputs received at the output layer from the inner layers at initialization. This observation allows us to formally extend our explanation by approximating a DGP at initialization using a SVGP. This enables formalizing a coordinate optimization algorithm that can be used to understand the parameter updates carried out in a DGP during the initial training iterations. This section considers DGPs with one single GP per layer, for illustration purposes. However, the exact same analysis and behavior is expected to hold when an arbitrary number of GPs is used per layer. The rest of the section is organized as follows: First, in Sec. 4.2.1, we describe the example setup, introduce our toy problem, and then proceed by examining different initializations and the whitened and non-whitened parameterizations. We also explain why the two different prior mean functions generate different predictive distributions at initialization in Sec. 4.2.2. After this, in Sec. 4.2.3 and Sec. 4.2.4, we analyze how the initial predictive distribution affects 23
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
model fitting and why it may lead to posterior collapse. This is done under a non-whitened with a reparameterized variational mean under the initialization S = 10−5 I in all layers. Then, we analyze different initialization settings and discuss the whitened parameterization. After this, in Sec. 4.2.6, we study the configuration of Salimbeni and Deisenroth (2017), where the inner layer’s variational covariance are initialized to Sv = S = 10−5 I, but the output layer variational covariances are initialized to S = Sv = I. Later, we discuss the role of the kernel parameters in Sec. 4.2.7, and the number of inducing points in Sec. 4.2.8. Finally, we briefly discuss the standard non-whitened parameterization implemented in GPYTORCH and GPJAX, the effect of a mixing matrix and warping functions. 4.2.1 ELL Term and Initial Predictive Distribution In a DGP, the ELL term is defined in Eq. (26). When training single-output DGPs with a homoscedastic Gaussian observation model, as outlined in Sec. 3.5, the ELL can be estimated via Monte Carlo sampling as: ELL =
N S X 1X n=1
S
s=1
L log N yn | fn,s , σ2 ,
(28)
L ∼ q(f L | f L−1 ), . . . , f 1 ∼ q(f 1 | xn ). While the expectation over the last layer, with fn,s n,s n,s n,s n,s for this particular noise model, can be integrated out analytically, our explanation remains simpler considering this Monte Carlo approximation. Recall that the required samples are obtained from the predictive distribution at each training point xn . Thus, we only need to analyze the propagated samples, which are independently generated for each xn . To understand the effect of the generated samples on the training objective, we consider a one-dimensional single-output regression problem, where the goal is to model a latent function f : R → R. We consider a DGP with 2 layers. That is, a hidden layer and an output layer. The number of GPs (units) in the hidden layer is set equal to one. The problem considered is displayed in Fig. 2. This is the toy problem introduced by Rudner et al. (2020). Blue points denote the training data. We observe that the target latent function is close to a step function, with sharp transitions that are challenging to model using the squared exponential kernel and standard GPs. Each sub-figure in Fig. 2 corresponds to different DGPs using the non-whitened parameterization implemented by GPFLOW, where the prior mean or the initialization of the variational parameters is different. In this figure, the DGP modes have not been trained, and we simply report how they behave at initialization. The first column of the figure shows the predictive distribution at initialization (mean and variance) of the GP in the hidden layer q(f 1 | X). See Eq. (20) and Eq. (21). Moreover, we also display joint samples from such a predictive distribution (using full covariances) Kqf (red dashed line) and from the marginal distribution, i.e. with a diagonal covariance given by diag(Kqf ) (green dotted line). Note that the samples from the marginal distribution (green dotted line) are not smooth, unlike the samples from the joint distribution (red dashed line). Recall that the marginal samples are the ones used to evaluate the ELL term, as mentioned previously. We also report the initial location of the inducing points as yellow points. The second column of Fig. 2 displays similar information, at initialization, but for the GP in the last layer of the DGP model (second layer). This is the output layer. Importantly, however, we consider here the
24
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
(a) ZERO prior mean DGP. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
0.4
0.6
0.8
1.0
(b) PCA prior mean DGP (only for the hidden layer). 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
(c) ZERO prior mean DGP with the proposed variational initialization.
Figure 2: From top to bottom, we display the predictive distribution of two layers DGPs with ZERO prior mean function (first row), with PCA prior mean functions in the hidden layer (second row), and with ZERO prior mean functions with the proposed variational initialization (third row). The hidden layer and output layer’s S is set to 10−5 I and the non-whitened parameterization implemented by GPFLOW is considered. The left column shows samples and predictive distribution of the hidden layer GP. The center column shows samples and predictive distribution of the output layer GP, if we condition on the hidden layer’s posterior mean, i.e. by ignoring the uncertainty in the predictive distribution of the hidden layer. The right column shows the predictive distribution of the full DGP obtained using the Monte Carlo approximation in Eq. (29). predictive distribution for each input point after computing the predictive mean of the first layer and letting that predictive mean go through the second layer. This ignores prediction uncertainty resulting from the GP in the first layer (only considers the nonlinearities given by the predictive mean), but gives a general idea about how the predictive distribution at the output layer changes with respect to the input of the DGP in a tractable form. In other 25
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
words, it is the predictive distribution of the output layer if we condition on the inner layer posterior mean (black line). Last, the third column shows the full predictive distribution of the DGP at initialization. This predictive distribution is obtained by letting the input data go through the two layers of the DGP. Since this distribution is intractable, we approximate it using S = 1000 Monte Carlo samples as follows: S 2 1X 2 E fn = fn,s , S
1 VAR[fn2 ] =
S
s=1
S X s=1
2 fn,s
2
2 − E fn2 ,
(29)
2 ∼ q(f 2 | f 1 ); f 1 ∼ q(f 1 | xn ). with fn,s n,s n,s n,s n,s When we look at the figure row-wise instead of column-wise, we observe similar information, but when the prior mean or the variational initialization is changed. Specifically, the first row of Fig. 2 shows the predictive distributions, at initialization of a DGP with the ZERO prior mean function in each GP. The second row corresponds to the PCA prior mean function (only in the hidden layer, in the output layer, the ZERO prior mean function is employed). In each row, the variational initialization considered sets S = 10−5 I in the inner layer and the output layer. We also set a length-scale of 0.1, and an output scale of 1.0. The noise variance σ 2 is set to a small value. Finally, the third row corresponds to a ZERO prior mean DGP with our proposed initialization, which we describe in detail later in Sec. 5. In all cases, the number of inducing points is set to 5, initialized via K-means (for both the inner and output layer). Fig. 2 shows that a different prior mean function, or a different variational initialization, has a strong effect on the resulting initial DGP predictive distribution (third column), but only in terms of the predictive standard deviation. In all cases, the DGP predictive mean is equal to zero at initialization. The predictive variance is small at the inducing points when the PCA prior mean function is used in the hidden layer GP, and when the proposed variational initialization is used. By contrast, when the ZERO prior mean function is used, the predictive variance is constant for each input point.
4.2.2 Explaining the Differences in the Initial Predictive Variance The most striking difference in Fig. 2 is found in the DGP initial predictive variance when ZERO and PCA prior mean functions are used (first row and second row, right column). The initial DGP predictive means are, however, similar and equal to zero. Eq. (20) indicates that when mr = 0, the predictive mean µqfr at the inner layer is 0 for the ZERO mean prior function, and X, i.e., the DGP inputs, for the PCA prior mean function. Thus, samples from the predictive distribution of the inner layer will be around zero (the top-left sub-figure) and around values of X (middle-left sub-figure), respectively. However, since the output layer prior mean is zero for both models, the output layer predictive mean is zero for each model at initialization (top-middle and center sub-figures). Thus, when we consider the full DGP predictive distribution, we obtain zero as the initial predictive mean (top-right and middle-right sub-figures). Regarding the initial predictive variance, there is a clear difference in the output layer. In the ZERO mean prior DGP, it is constant. By contrast, in the PCA prior mean DGP, we observe a wavy shape, with zero variance at the inducing points. For the input layer, the initial predictive variance next to an inducing point zi will take the value corresponding to 26
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
the i-th entry of diag(S). In particular, in Eq. (21), when X = Z, the predictive covariances are equal to S. Since the variational covariance is initialized to S = 10−5 I, the predictive variance is almost zero at the inducing point locations. This is something that happens both in PCA and ZERO mean prior DGP models. The reason is that the predictive covariances do not depend on the prior mean in Eq. (21). For input points far from the inducing points, one should expect that Kxi Z ≈ 0. Thus, the predictive covariance in Eq. (21) will approximately be given by Kxi xi . This implies that any point far from the inducing points will get a predictive variance given by the kernel output scale parameter, σo , which is initialized to 1.0. This is visible in Fig. 2 for the inner layer (left column) in each DGP model. We do not observe a value exactly at σo = 1.0 since Kxi Z is not exactly zero, due to the value of the length-scale. Regarding the predictive variance of the GP at the output layer (middle-column in Fig. 2), we observe that in the ZERO prior mean DGP it is constant. Recall that, for any layer after the first one, the predictive covariance depends on the inducing locations Z (which, at initialization, are the same across layers, see Sec. 3.5), and the samples of the previous layer. Importantly, however, the samples from the previous layer will be similar to X for the PCA prior mean DGP, and to a vector of 0 for the ZERO prior mean DGP (ignoring the prediction uncertainty from the first layer). Consequently, for the ZERO mean prior function, the predictive covariances are: −1 −1 −1 Kqf = K00 − K0Z KZZ KZ0 + K0Z KZZ SKZZ KZ0 ,
(30)
for any point from the input space, even at the inducing points locations. This results in similar predictive variances for each point in the input space. Furthermore, the generated samples at this layer are expected to be constant in terms of the DGP inputs, as illustrated by the red dashed curve in the top-middle sub-figure in Fig. 2. The reason is that the predictive mean of the first layer will map each input to a value equal to zero. This implies a posterior covariance Kqf with the same values in all its entries. Since the marginal variances are constant, we also observe a constant variance across independent samples, as illustrated by the green dotted curve. Whether the output layer predictive variance is close to the output scale kernel parameter or the initial values specified in S depends on how far 0 is from the inducing points. A greater number of inducing points will result in a posterior more concentrated at the initial S values. The reason is that there is a bigger chance for an inducing point being close to 0. We note, however, that for the PCA prior mean DGP, things are different. In this case, at initialization, the hidden layer predictive mean is equal to the identity function. Namely, it outputs X given X as an input. Thus, the predictive variance for the output layer has exactly the same behavior as the input layer for the ZERO prior mean DGP. Specifically, rather than always evaluating on 0, we evaluate on X when X is the DGP input. The consequence is that the output layer predictive variance has the same wavy shape as shown in the input layer, but around a mean value of zero, since the output layer prior mean is the ZERO mean function. The last row in Fig. 2 shows that our proposed initialization results in a ZERO mean DGP with a similar behavior as that of the PCA prior mean DGP at initialization. While our reasoning regarding the shape of the output layer posterior variance relies on ignoring the predictive variance from the hidden layer, we expect a similar behavior in the full DGP predictive distribution. This is indeed the case, as illustrated by Fig. 2 (right 27
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
column). The reason is that the initial predictive variances are fairly small at initialization. ∼
Thus, the ZERO prior mean DGP will output values 0 close to 0 in the hidden layer, for each input in X. Similarly, the PCA prior mean DGP will output x̃n close to xn , for each input xn in X, in the hidden layer. Overall, this implies that, at initialization, the ELL for each training sample in the ZERO mean DGP is evaluated using a sampling procedure similar to ∼
∼
1 | xn ) and f 2 ∼ q(f 2 | 0). By contrast, the PCA mean DGP is evaluated using a 0 ∼ q(fn,s n,s n,s 1 | xn ) and f 2 ∼ q(f 2 | x̃n ). sampling procedure similar to x̃n ∼ q(fn,s n,s n,s
4.2.3 Virtual Dataset Seen by a ZERO Prior Mean DGP at Initialization An intuitive argument to understand why the ZERO prior mean DGP is likely to have a posterior collapse problem is related to the virtual dataset that it sees during training, just after initialization. Specifically, we have seen that the output layer of the ZERO prior mean ∼
DGP receives input values 0 close to 0, for any point in the training set. This means that the ∼
virtual dataset that is being used to optimize the ELL is made up from pairs {(0, y n )}N n=1 . Thus, from the model’s perspective, every input vector is the same (all of them being close to the zero vector). Therefore, in the output layer, there is no variability in xn that could explain the variability in y n . As a consequence, for the model, it is easier to represent y n purely as noise and set the KLD equal to zero. 4.2.4 Impact of Initial Predictive Variances in the ELL How do the differences in the initial predictive variance influence the ELL? An intuitive explanation starts by noting that the ELL term involves the expected value of the score function log N y n | fnL , σ 2 for each training point, under the marginal distribution q(fnL ) over the last layer. Assuming q(fnL ) has finite mean and variance, then the ELL can be obtained analytically. That is, ELL =
=
=
N X n=1 N X n=1 N X n=1
Eq [log N y n | fnL , σ 2 ]
(31)
1 log N y n | Eq [fnL ], σ 2 − 2 VAR(fnL ) 2σ
(32)
√ L 1 1 2 − log σ 2π − 2 (y n − µL Tr Kqf . qf n ) − 2σ 2σ 2
(33)
Here, while we use a similar notation to that of the GP predictive mean and variance, L represent the mean and covariance of a non-Gaussian distribution. µqf n = Eq [fnL ] and Kqf Therefore, the learning algorithm finds parameter values so that this score function is high, which, broadly speaking, implies that the resulting predictive mean is close to the targets, and the predictive variance VAR(fnL ) is as small as possible, ideally 0. The optimal parameter σ 2 is the expected value of the average squared distance between the target y n and the DGP output fnL , see appendix A. In other words, it is the sample variance given optimal posterior mean, i.e. how much uncertainty remains under-explained, plus a contribution that comes from the trace term. 28
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
The right column of Fig. 2 shows that the PCA prior mean DGP provides a predictive distribution at initialization that can highly disagree with the observed data. For example, at nearly xn = 0.5 (which coincides with an inducing point location i), the variance is almost zero while the corresponding y n is around 1.2. This will incur a very low ELL term, and so the model will be forced to update the parameters to improve the prediction of y n . At this point, xn , the kernel parameters (ν ) and the inducing point locations (Z) do not influence the ELL term. This is so because the predictive variance is given by S (since we are computing predictions at an inducing point) and the predictive mean is zero due to the initialization. More precisely, as shown by the two equations below, any hyper-parameter gradient resulting from the predictive mean will be multiplied by 0 (due to the predictive mean being linearly related to mr = 0) and any hyper-parameter gradient from the predictive variance will be zero, as initial predictive variances are given just by S. In particular, note that if we approximate the predictive mean of the DGP through Monte Carlo, then using Eq. (20) and Eq. (21), the gradients of the ELL w.r.t. the hyper-parameters and inducing points at this location xn at the last layer are given by:
1 ∂ ELLn = ∂ν S
S X s=1
h iT T −1 L−1 K Z ZZ 1 n 1 ∂ diag(S)i ∂ Kfn,s (y − µL =0 mr − 2 qf n,s ) 2 σ ∂ν T 2σ ∂ν
h iT T −1 L−1 K Z ZZ ∂ ELLn 1 1 n 1 ∂ diag(S)i ∂ Kfn,s L = (y − µ =0 n,s ) mr − 2 qf T ∂vec Z S σ2 2σ ∂vec Z ∂ [vec Z] s=1 S X
(34)
(35)
where we have omitted the L superscript from some of the variables for clarity. Note that by the chain rule, the gradients at each layer will all be multiplied by mL r . Thus, such a 2 low ELL term can only be increased by moving m, S and σ to compensate and explain the data, since they are the only gradients taking a value different from zero. More precisely, since the predictive variance is very small (due to the initial value of S) , the noise variance σ 2 does not influence the ELL term for the input xn via − 2σ1 2 VAR(fnL ). Thus, the ELL is 2 . From this point, making the maximized by setting this parameter to 1/N || Y − µL || qf 2 √ variance greater will necessarily make the ELL small, since limσ→∞ − log σ 2π = −∞. On the other hand, for an initial small σ 2 at initialization, gradient ascent pushes S to stay at low values. This is so because the gradient w.r.t. S at the i − th inducing location is −1/2σ2 and zero at other point. Thus, it is expected to be large, and because S cannot be negative, this keeps S to small values. This is to be expected since maximizing the ELL w.r.t. S is achieved by a small trace term, which is achieved by keeping the diagonal terms in S small. Thus, the only way we can make the ELL higher is by making the predictive mean to be closer to the observed targets, which necessarily implies moving m away from its initial value. Therefore, it is expected to be more efficient to update m to improve the ELL term. Once m is different from 0, the kernel hyper-parameters and the inducing points start being optimized. This avoids having big variational variances S (note that the ELL forces S to be small) or big noise values σ 2 during the first optimization iterations, reducing the possibility of explaining everything as noise and, consequently, the posterior collapse problem. 29
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
In contrast, for the ZERO mean prior DGP, at initialization, we have a constant predictive variance. This implies that there are no predictions that are in greater disagreement with the observed data than others. Thus, since the variational mean m and the variational covariance S in the last layer are initialized close to the prior mean and covariances, it is easier to keep them to these values. Moreover, in this case, the predictive variance will just depend on values coming from the previous layer, which are close to 0 for each input (see the previous section). Thus, making the term 2σ1 2 VAR(fnL ) small is more easily done by increasing the noise variance than by moving the inducing points or the kernel parameters so that KXZ ≈ 0, KXX ≈ 0 is satisfied. This results in a low KLD term, and the algorithm is expected to maximize the ELL simply by increasing the noise parameter, indicating that everything is noise. That will maximize the ELL term for all points without incurring a strong disagreement in terms of the KLD. The consequence is that it is more likely to have a posterior collapse. From this analysis, we can understand why more inducing points or a small output scale parameter σo are likely to induce less collapse. These design choices contribute to making the trace in the ELL small, forcing the model to move m from its initial value. 4.2.5 Whitened Parameterization So far, we have analyzed the initial predictive variance of a non-whitened reparameterized DGP implemented in GPFLOW. In this section, we show that the posterior collapse problem may also appear under the whitened parameterization, for similar reasons. Specifically, Fig. 3 shows the equivalent predictive distributions to those displayed Fig. 2, for different DGPs, using the whitened parameterization. All variational parameters and hyper-parameters are initialized as before. We observe that the initial predictive distributions resulting from both parameterizations look very similar. Therefore, a similar behavior should be expected under the whitened parameterization. Namely, the ZERO prior mean DGP is more likely to have a posterior collapse problem. We also observe how the proposed initialization works with this parameterization as well. In general, however, the predictive distribution under non-whitened and whitened parameterizations will not always be that similar. Initialization of the output layer variational covariances to I reveals some differences, as illustrated later on. The results obtained, however, show that the origin of the problem is independent of the parameterization, and that the ZERO prior mean DGP is consistently more likely to have a posterior collapse problem. This analysis also motivates a set of ablation studies carried out in the experimental section and provides guidance for more principled initializations of both variational parameters and model hyper-parameters. 4.2.6 Initializing S = Sv = I at the Output Layer First, we consider the non-whitened case. Fig. 4 shows the predictive distributions that are obtained in this case. This figure illustrates that in the output GP the predictive variance at the inducing points, when ignoring the prediction uncertainty from the inner layer, is now larger and equal to 1, as expected from S = I (middle column). We observe, however, that in the PCA prior mean DGP, the variance is slightly smaller than the prior variance when the input lies between two inducing points. This seems counterintuitive, as we could expect that far from the inducing points the uncertainty is close to that of the prior, i.e., 30
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
(a) ZERO mean DGP 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(b) PCA mean DGP 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(c) ZERO mean DGP with the proposed variational initialization
Figure 3: Same predictive distributions as the ones displayed in Fig. 2, but when the the whitened parameterization is used. KXX , with possibly added uncertainty coming from S. Nevertheless, it is a consequence of the variational posterior S providing more uncertainty than what the prior kernel and the location of the inducing points reduce. Specifically, we initialize S = I, while the prior covariances are given by KZZ . Therefore, −1 −1 Kqf = KXX − KXZ KZZ KZX + KXZ KZZ SK −1 KZX , | {z } | {z } | {z ZZ } (1)
(2)
(3)
−1 −1 = KXX − KXZ KZZ [KZZ − S] KZZ KZX ,
(36)
may lead to a smaller predictive variance than the prior, i.e., KXX . To see so, the predictive variance has three terms. The first term (1) is the prior variance. The second term (2) is how much variance is reduced by the prior distribution and the position of the inducing points (which has no uncertainty associated). The third term (3) reintroduces uncertainty coming 31
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
(a) ZERO mean DGP 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
(b) PCA mean DGP
Figure 4: ZERO (top) and PCA (bottom) DGP models with the output layer variational covariance initialized to S = I, in the GPFLOW’s unwhitened parameterization. from our actual lack of knowledge about the inducing point function evaluations. Thus, despite uncertainty arising from the inducing points (through S) and from the prior (through KXX ), the information provided by the inducing point locations together with the prior −1 kernel (through KXZ KZZ KZX ) might be stronger, consequently contributing to reducing the predictive variance at points far from the inducing points. This can be formally stated through the Loewner order of positive definite matrices. Formally, whenever (KZZ − S) is −1 positive definite, then by Loewner order KZZ ≻ S. Since for any point vT = Kxn Z KZZ is T T a vector, then v KZZ v > v Sv which implies that the uncertainty will never exceed the prior variance far from the inducing points. This will be true whenever diag(S) ≤ σo . On the other hand, the ZERO prior mean DGP does not show this behavior, as in the first layer each input is mapped to a vector of zeros. Finally, independently of the prior mean used, the full DGP predictive variance also increases, as expected, since now S = I at the output layer. The third column in Fig. 4 shows the corresponding predictive distribution. Because the full DGP predictive variances are bigger in the output layer when we set S = I, one should expect that it is easier for the DGP model to suffer from a posterior collapse problem. This is confirmed by our experiments later on. In addition, in this setting, the ZERO prior mean DGP is also expected to be more likely to suffer from a posterior collapse problem than the PCA prior mean DGP. The reason is that the first layer remains unchanged. ∼
Therefore, the ZERO prior mean DGP will still observe the virtual dataset {(0, y n )}N n=1 at initialization, which will make difficult explaining y n from xn . When we consider a whitened parameterization, and we set Sv = I, the initial predictive distributions that are obtained are very similar to those obtained in the unwhitened case, 32
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
(a) ZERO mean DGP 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
(b) PCA mean DGP
Figure 5: ZERO (top) and PCA (bottom) DGP models with the output layer variational covariance initialized to Sv = I, in the whitened parameterization. as illustrated by Fig. 5. In this case, however, the predictive variance of the output GP is constant and equal to 1 in the PCA prior mean DGP, when the prediction uncertainty from the inner layer is ignored (bottom-middle figure). The reason is that now Sv = I equals the prior whitened covariances associated with the inducing points and, consequently, at any point, at initialization, the predictive covariances are given by those of the prior, i.e., KXX . Note that: T −1 −1 Kqfw = KXX − KXZ KZZ KZX + KXZ L−1 ILZZ KZX ZZ −1 −1 = KXX − KXZ KZZ KZX + KXZ KZZ KZX
= KXX
Note that once we start optimizing the variational parameters, and Sv gets a different value from I, the same wavy shape appears, as we illustrate in Fig. 23 in the appendix. Thus, the full DGP predictive distributions look very similar to the non-whitened case in this setting, and one should expect a similar behavior of each DGP with respect to the posterior collapse problem. One may ask what may happen when we set S = Sv = I also at the inner layer. In this case, we should expect an increment in the initial predictive variance at that layer, for each parameterization, see Fig. 24 in the appendix. This bigger variance will result in an even bigger initial predictive variance at the output layer. Therefore, such an initialization is expected to lead to a higher probability of having a posterior collapse problem, and is hence not recommended. This initialization provides samples in the inner layers, for the ZERO prior 33
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
mean DGP model, which might be more different from 0 than when S = 10−5 I. Thus, at first, it might seem to help to alleviate the posterior collapse problem since now it is less likely ∼
that the output layer receives as input a dataset of virtual pairs {(0, y n )}N n=1 . However, in general, we have observed that this additional noise in the inner layer during training usually leads to sub-optimal models. In our experiments, we show that these instabilities, with independence of the parameterization, might lead to posterior collapse. More precisely, in the experimental section, we show that a very deep PCA prior mean DGP with the whitened parameterization can lead to posterior collapse when extra noise is introduced in the inner layers, and that reducing this noise alleviates the problem. We note that some works also consider random initialization of S (Hamelijnck et al., 2019). In this setting, we shall expect a similar behavior to the one already described, for both the non-whitened and the whitened parameterizations. Namely, that the ZERO prior mean DGP is more likely to suffer from a posterior collapse problem. The reason for this will be again the difficulty to model the dependence of y n with respect to xn , at initialization, as a consequence of the first layer predicting values close to 0 for any input. However, one should be careful when using a random initialization of S, since that may result in inconsistent variational covariances that make the predictive variance go beyond the prior, as illustrated by Fig. 26 in the appendix. This happens when random S results in KZZ − S not being positive definite. From this section, we conclude that spherical initializations S = Sv = αI with small α ≤ σo , e.g., 10−5 , yields an initialization that remains conceptually well-motivated, keeps the model within a plausible region of the parameter space, and avoids convergence issues associated with pathological starting points. This is because KZZ − S results in positive definite matrices in both parameterizations, resulting in a reduction of predictive variance beyond the prior. The suitability of this initialization is confirmed later on in the experimental section. 4.2.7 Impact of the Kernel Hyper-parameters The impact of the kernel hyper-parameters in the initial predictive distribution can be easily understood. In the RBF kernel, the output scale parameter, σo , is directly related to the diagonal of KXX , which is the initial prior variance. Reducing σo is expected to reduce the initial predictive variance when we move away from the inducing points. At the inducing points, however, the predictive variance will still be given by S. The length-scale, on the other hand, affects the relative distance between the points, making the prior covariances KXX more or less isotropic. It also has an impact on how quickly the predictive mean and variance go back to those of the prior as we move away from the inducing points. This will be the case for both the whitened and non-whitened parameterizations. Therefore, changing the output scale or the length-scale will have a small impact on the posterior collapse problem of the ZERO prior mean DGP. More precisely, the influence of xn on the DGP predictive distribution will still be rather small, making the posterior collapse likely to occur. By contrast, in the PCA prior mean DGP there will still be a significant dependence of the predictive distribution of the output GP w.r.t. xn , reducing the possibility of a posterior collapse. However, if posterior collapse appears, it might be beneficial to reduce σo . This will reduce the variance in the inner layers at locations far from the inducing points, reducing overall noise in the 34
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
optimization process. At the output layer, it will also reduce the predictive variance, and thus the trace term of the ELL, which results in a larger ELL. Summing up, changing the kernel hyper-parameters will affect mostly the predictive distribution in regions of the input space far from the inducing points or, similarly, far from the observed data, which is expected to have a smaller impact on the evaluation of the ELL term during the first iterations of the optimization problem, as compared with the initial choice of the variational parameters. 4.2.8 The Effect of More Inducing Points We now study the effect of increasing the number of inducing points when S = 10−5 I in all layers. Fig. 6 displays the initial predictive distributions of the ZERO prior mean DGP and the PCA prior mean DGP when 7 and 20 inducing points are considered, instead of 5. We only consider here the non-whitened parameterization. We observe that when a large number of inducing points are used, the predictive variance becomes closer to zero, which is the posterior variance specified by S associated to the inducing points. In this setting, we can expect the ELL term to be very low, because of the small variance. As we saw in Sec. 4.2.4, this favors not making the noise variance σ 2 large in order to explain the observed data, but updating the variational mean parameter. Note that the more inducing points, the smaller will be the trace term in the ELL if S is initialized with a small α. This is confirmed by our coordinate update approach in the next section. For completeness, appendix C.1 provides figures with more inducing points and other initializations for S in the inner and output layers. In those figures, we observe that with more inducing points, the same behaviors we have been analyzing are observed. For example, S with smaller α is also preferable. We also observe a wiggling shape appearing in the ZERO DGP. This is because with more inducing points, it is likely that function evaluations around 0 are closer to the inducing points, making KXZ = ̸ 0. In consequence, it is expected that more inducing points help in avoiding the posterior collapse problem in the ZERO DGP, which we confirm in the experimental section. Therefore, increasing the number of inducing points should favor avoiding the collapse problem. Moreover, it is also expected to make the dependence of the predictive variance on the output scale parameter σo less relevant since now the predictive variance is likely to depend mostly on S. We confirm these arguments in the experimental section. 4.2.9 The Standard Non-whitened Parameterization So far, we have analyzed what happens under the non-whitened parameterization by GPFLOW, in which the tuned parameter of the posterior approximation is the posterior mean minus the prior mean (see Sec. 3.3.3). However, as we mentioned earlier, both GPYTORCH and GPJAX implement the standard non-whitened parameterization. For a PCA prior mean DGP, this parameterization results in the following predictive mean at initialization: −1 µqf = X + KXZ KZZ (0 − Z) .
(37)
Thus, in the first layer GP of the DGP we will not observe a straight increasing line, as we have seen previously (see Fig. 2 middle row). In particular, right at the inducing points, 35
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
(a) ZERO mean DGP with 7 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(b) PCA mean DGP with 7 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(c) ZERO mean DGP with 20 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(d) PCA mean DGP with 20 inducing points.
Figure 6: Initial predictive distributions of a ZERO prior mean and a PCA prior mean DGP with 7 and 20 inducing points using the GPFLOW non-whitened parameterization.
the predictive mean will be µqf = 0, as indicated by Eq. (37). Fig. 7 shows the predictive mean of the inner layer GP of a DGP with the standard non-whitened parameterization and the PCA prior mean, for different numbers of inducing points. We observe that, as 36
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
5 inducing points
10 inducing points
1.0
1.0
0.5
0.5
0.0
0.0
−0.5
data Posterior mean Inducing points
−1.0 0.00
0.25
0.50
0.75
−0.5 −1.0
1.00
0.00
20 inducing points 1.0
0.5
0.5
0.0
0.0
−0.5
−0.5
−1.0
−1.0 0.25
0.50
0.75
0.50
0.75
1.00
100 inducing points
1.0
0.00
0.25
1.00
0.00
0.25
0.50
0.75
1.00
Figure 7: Predictive mean of the inner layer GP of a DGP with the standard non-whitened parameterization for different numbers of inducing points and the PCA prior mean function. the number of inducing points increases, the predictive mean becomes nearly constant and equal to zero. It is only equal to X when we are far away from the inducing points. This indicates that under this parameterization, a PCA prior mean DGP will behave similarly to a ZERO prior mean DGP in the inner layers, if the number of inducing points is big. Therefore, it is expected to be more likely to suffer from the posterior collapse problem than the non-whitened parameterization of GPFLOW. Even though Rudner et al. (2020) reports successfully training DGPs with the standard non-whitened parameterization, our analysis shows later on that such a parameterization is of less practical interest since it results in unstable optimization. 37
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
4.2.10 Mixing Matrices and Warping Functions While the original formulation of DGPs considered independence between processes within a layer, we can easily introduce dependencies either using mixing matrices (Jankowiak and Gardner, 2019) or warping functions (Sáez-Maldonado et al., 2023). Note that random initialization of these transformations might impact both the PCA and ZERO prior mean DGP predictive distribution at initialization. Thus, a practitioner might need to carefully choose the initial values. For mixing matrices, it seems reasonable to assume that the processes are (at the beginning) independent by initializing the mixing matrix to the identity, and let the marginal likelihood estimate any possible dependency. Also, warping functions are usually initialized to an identity mapping (Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023). In those cases, at initialization, the DGP models will behave in the same way as we have analyzed. 4.3 An Approximate Explanation for the Posterior Collapse Problem So far, we have used an intuitive explanation to understand the source of the posterior collapse problem, under different initializations, on a toy dataset. We have found that the main reason for the posterior collapse problem is that, at initialization, the output GP in the DGP sees a virtual dataset of input points very similar to zero. Even if most computations involving DGPs with the DSVI algorithm are not analytically tractable, we show in this section that we can make some approximations that yield tractable analytical conclusions. More precisely, at initialization, we approximate the output layer of a DGP with an SVGP that receives as an input either 0 or xn , and then analyze the corresponding coordinate gradient ascent updates of the variational and noise parameters. This gives extra insights into the reasons for posterior collapse. In this section, we focus on the standard non-whitened parameterization. 4.3.1 Approximating at Initialization a DGP Using a SVGP While in Sec. 4.2.4 we showed via the ELL term how a DGP is more or less likely to have a posterior collapse problem, we ignored the role of the KLD term. To take it into account, we note that many of the initializations considered in the previous sections lead to a small predictive variance in the inner layer of a DGP, making it nearly zero. This included a small variational variance S = 10−5 I, a big number of inducing points, or a small output scale kernel parameter. This means that, at initialization, we will mainly propagate the predictive means µqf until the last layer. Given this, we can see that, at initialization, the output layer resembles a sparse Variational GP (Titsias, 2009; Hensman et al., 2013). As mentioned earlier, for the PCA prior mean, the SVGP will be trained on {(xn , y n )}N n=1 , while the ZERO prior mean DGP will be trained on {(0, y n )}N . Thus, by studying the SVGP’s objective at n=1 initialization and the gradient updates, we can draw conclusions on what early local minima a DGP might find during optimization, revealing interesting properties. A SVGP can be trained by coordinate ascent on the variational and model parameters, where we mix exact and gradient-based updates. In particular, Titsias (2009) showed that 38
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
the optimal variational updates on m and S are: −1 1 KZZ , (38) Sopt ← KZZ KZZ + 2 KZX KXZ σ 1 −1 mopt ← 2 Sopt KZZ KZX y , (39) σ for the non-whitened parameterization. Furthermore, the optimal noise parameter σ 2 is: 2 σopt ←
1 || Y − µqf ||22 + N | {z } aleatoric
Tr Kqf | {z }
,
(40)
aleatoric + approximation error
where µqf and Kqf are the predictive mean and the predictive covariance, respectively. Appendix A provides the derivation. The coordinate update on σ 2 is the same for the whitened and non-whitened parameterizations. We just need to plug in the corresponding predictive mean and covariance. The coordinate update on the noise in Eq. (40) reveals that optimality is achieved as the expected value of the squared distances between the targets and the expected mean (the sample variance), plus another term which has an aleatoric component, plus some uncertainty that comes from representing the problem through the inducing points. Note that attending to Eq. (33), the optimal ELL is achieved when this trace term is 0. This is only achieved when S = 0 and X = Z. Attending to the coordinate update for S, this term is zero when either σ 2 → 0 or the kernel scale parameter σo is zero. σ 2 contributes as a source of aleatoric uncertainty, while σo is only driven to zero by the marginal likelihood if the data comes from a constant function. Thus, this is a source of aleatoric uncertainty. When S is zero, we require X = Z for the trace to be zero. This is variance coming from an approximation error. It is not aleatoric nor epistemic, because it is not made zero by collecting more data, nor is noise implicit in the data to be modeled. The immediate consequence of this observation in the ZERO prior mean DGP is that the h i contribution from Tr Kqf is larger since the predictive variance is constant across training data, as illustrated in Fig. 2 and Fig. 3. By contrast, the PCA prior mean DGP shrinks the predictive variance at some points. Clearly, the variational or hyper-parameters initialization will also contribute to this term, but given similar initializations, the ZERO prior mean DGP will increase the predictive variance more. Furthermore, in the limit of σ 2 being arbitrarily big, the coordinate updates on the variational parameters will result exactly in those that make the KLD be zero, i.e. m = 0, S = KZZ , resulting in a posterior collapse problem. See Eq. (38) and Eq. (39). Thus, a higher σ 2 is expected to drive the variational parameters towards the prior. Coordinate Updates for the ZERO prior mean DGP: For a ZERO mean DGP, we can see the output layer as a SVGP that receives 0 as every input. Using the initial variational parameters m = 0, S = αI, which make µqf = 0, the coordinate update is given by: 2 σopt ←
N
N
n=1
n=1
1 X n 2 1 X −1 −1 −1 (y ) + K00 − K0Z KZZ KZ0 + αK0Z KZZ KZZ KZ0 . N N 39
(41)
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Note that if we condition on 0 for each training location xn , the diagonals from the marginal variational posterior covariance evaluated at all pairs X are exactly the same, which results in N equal terms being added up. As we already saw, this posterior covariance can be above or below the prior variance, depending on how much uncertainty from the prior we are able to subtract due to the inducing point locations. For S = αI, with α ≤ σo , is below due to Lowener’s order. Assuming the initial output scale kernel parameter is σo = 1, we have: # N X 1 −1 −1 −1 2 σopt ← (yn )2 + 1 − K0Z KZZ KZ0 + αK0Z KZZ KZZ KZ0 , N "
(42)
n=1
which implies that initializing with smaller α values should favor smaller noise variances, reducing the posterior collapse problem, a similar conclusion as the one from Sec. 4.2.6. Coordinate update for the PCA prior mean DGP: For the PCA prior mean DGP model, the coordinate update is different, since we now observe xn as an input. Thus, we have that when m = 0, S = αI 2 σopt ←
N
N
n=1
n=1
1 X n 2 1 X −1 −1 −1 KZZ KZxn , (43) KZxn + αKxn Z KZZ (y ) + Kxn xn − Kxn Z KZZ N N
which for Kxn xn = 1 we have: # # " N N X X 1 1 −1 −1 −1 2 σopt ← KZxn − αKxn Z KZZ KZZ KZxn . (44) (yn )2 + 1 − Kxn Z KZZ N N "
n=1
n=1
The difference with respect to the ZERO prior mean DGP is clear. Here, the resulting noise variance parameter does not receive a constant contribution but an average value at each training point. Recall that some of them will have a smaller predictive variance when α = 10−5 , due to the input points being closer to the inducing points. The value of α will have a similar effect as the one described previously. Namely, a larger α will increase the noise variance and is expected to make the posterior collapse problem more likely. However, as we observed previously, the PCA prior mean DGP with α = 1 slightly shrinks the predictive variance compared to the ZERO model, and so it is expected to suffer less from posterior collapse under this initialization as well. In conclusion, the PCA prior mean DGP is expected to result in a smaller noise variance estimate than the ZERO prior mean DGP and to reduce the posterior collapse problem, when the same initialization is used for both models. Our results also highlight that a larger number of inducing points is more likely to make the estimate of σ 2 smaller, even more so when α = 10−5 , and to reduce the likelihood of a posterior collapse problem. More precisely, more inducing points will result in smaller predictive variances for each xn . This analysis reveals that an initial small variance kernel output scale parameter σo also favors small values for σ 2 . Only in the particular case of the inducing points being close to 0, the variance of the ZERO DGP will be smaller. However, the initialization of inducing points very often tends to cover the whole training domain and, hence, this pathological case is highly unlikely in practice. 40
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Coordinate updates for the variational parameters: We have shown that the ZERO prior mean DGP is expected to consistently learn higher variances than the PCA prior mean DGP because, as we mentioned, the noise variance parameter receives a constant contribution rather than an average value at each training point. However, the resulting noise variances will not be the only cause of the posterior collapse. By looking at the coordinate updates for the variational parameters in Eq. (38) and Eq. (39), we observe that Sopt will tend to KZZ also when KZX KXZ is a matrix of zeros, with independence of the value of the noise variance σ 2 . A careful look at how this matrix is computed reveals that it will be closer to a matrix of zeros for the ZERO prior mean DGP. Specifically, each position in the matrix is given by: (KZX KXZ )ij =
N X
K(zi , xn )K(xn , zj ) .
(45)
n=1
Unless an inducing point is close to 0 (unlikely for a small number of inducing points), most of the entries of the matrix will be close to zero for the ZERO prior mean DGP, for a small enough length-scale. This will make the ZERO prior mean DGP more likely to collapse the variational covariance to the prior KZZ . Consider now that the optimal value for S is approximatelyKZZ , the optimal update for the variational mean is: mopt ←
1 KZX Y . σ2
(46)
Again, given that KZX takes smaller values for the ZERO prior mean DGP and, in conjunction with σ 2 taking higher values, the variational mean will also tend to the prior, i.e., m = 0, more likely for the ZERO prior mean DGP than for the PCA prior mean DGP. When this collapse happens, the coordinate updates on the noise variance parameter are given by: N
2 σopt ← σo +
1 X n 2 (y ) N
(47)
n=1
for each model (when Kxn xn = K00 = σo ). See Eqs. (41),(43). Thus, both models learn the same noise variance in this scenario. These results also show why a small number of inducing points is more likely to result in a posterior collapse. With fewer inducing points, KXZ is more likely to be close to 0 for many training points. Common sense also supports this reasoning. With very few inducing points, most of the targets will have a high predictive variance given by the prior (i.e. a bigger trace term), and the model will tend to learn everything as noise. Overall, coordinate updates on the variational parameters also support why the ZERO DGP is more likely to suffer from posterior collapse. 4.3.2 Simulating the Coordinate Updates We have shown that coordinate updates make the ZERO prior mean DGP learn, in general, higher noise variances than the PCA prior mean DGP. This, in conjunction with the coordinate update on the variational parameters, where KXZ → 0 in some settings of the ZERO prior mean DGP, indicates that the model is more likely to collapse. However, showing that 41
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
this actually leads to a posterior collapse after some optimization steps can only be done experimentally because: 1) the limit σ → ∞ is never achieved since it results in ELL → −∞, 2) achieving KXZ = 0 depends on the kernel parameters, number of inducing points and its location, 3) analyzing the convergence of the iterative maps that result from the coordinate updates cannot be done analytically. Therefore, we perform several simulations under different conditions to validate our claims, fixing the kernel parameters and including points to the same initial values we have been using in this section. Since the order in which coordinates are updated matters, we perform two different ordering schemes. The first ordering scheme implies updating S − m − σ 2 , and the second one uses the order σ 2 − S − m. We will see that regardless of the order, the ZERO prior mean DGP consistently learns higher noise variances and collapses the variational parameters to the prior. Before proceeding with the simulation, we analyze the two ordering schemes to validate what we should expect from the simulations. Order update S − m − σ 2 : This order can consistently lead to smaller noise variances. Specifically, for the initial value m = 0, changes in m necessarily maximize the ELL. This is because for m = 0 the KLD in the ELBO takes its minimum possible value w.r.t. this parameter. Thus, maximizing the ELBO w.r.t. m implies that the ELL is increased more than the KLD, after the update. Note that maximizing the ELL w.r.t. m can only be achieved by making the term 1/2(y n − µL n,s )2 smaller. This is exactly the sample variance term appearing in the update qf for σ 2 in Eq. (40). Since this term is necessarily smaller than that of m = 0, because we have shown that the ELL is increased, updating first m and then σ 2 consistently provides smaller variance estimates than updating first σ 2 and then m. By contrast, it is not possible to guarantee that updating S first will result in a trace term in Eq. (40) that is smaller than in the previous step, and hence in smaller noise variances, unless S is initialized to the prior. However, we can show that updating S first results in a variance that never exceeds the prior variance. Specifically, if we replace Sopt into the predictive covariances we have: −1 Kqf = KXX − KXZ KZZ KZX + KXZ
−1 1 KZZ + 2 KZX KXZ KZX . σ
(48)
To show that any update on S does not increase the variance beyond the prior, we need −1 −1 to show that KXZ KZZ + σ12 KZX KXZ KZX . KZX is always smaller than KXZ KZZ −1 T For a single point x, v = KXZ is a column vector. Thus, showing that v KZZ v ≥ −1 vT KZZ + σ12 KZX KXZ v implies showing that KZZ + σ12 KZX KXZ ⪰ KZZ . Since 1 KZZ + σ2 KZX KXZ − KZZ = σ12 KZX KXZ is positive definite4 , we can use Lowener’s order to show that the inequality holds. A similar procedure can be followed to show that Sopt never exceeds the prior covariance KZZ . One just needs to apply the Woodbury matrix identity. This shows that the predictive variances will never exceed the prior variances given by KXX . Thus, updating σ 2 first may consider bigger Kqf than the prior variance, if S is initialized in a pathological way (i.e., with higher uncertainty than the prior, KZZ − S is not positive definite). Furthermore, this proof shows that Sopt in Eq. (38) will never result 4. Note that KZX KXZ is the product of a matrix by its transpose.
42
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
in a value that makes the KLD exactly zero, i.e., S = KZZ , unless KXZ = 0 is zero. This is because the variance can never grow to ∞ and the update never exceeds the prior since KZZ ⪰ Sopt with equality only when KXZ = 0. In summary, this update order leads to more controlled predictive variances (which never go beyond the prior), independently of the initial variational parameters, and should be preferred. Order update σ 2 − S − m: This order scheme is sensitive to the initial variational parameters. Attending to our previous results, for badly selected initial S, we could add variance beyond that of the prior, or be exactly the prior if S = KZZ . Also, initializations with S = αI and α = 1 will provide higher variances when compared to α = 10−5 . Posterior collapse can happen here earlier since higher variances will contribute to making the matrix 1/σ 2 KZX KXZ be closer to all zeros. As in the other order, total collapse can only happen if KXZ = 0. Simulations: We simulate different scenarios on the toy problem described in Sec. 4.2. The inputs for the PCA prior mean sparse GP are xn . In the case of the ZERO mean sparse GP, they are a vector of zeros. Fig. 8 report the optimization of the ELBO using coordinate updates for 5 inducing points (see Fig. 27 in the appendix for simulations using 64 inducing points). In the figure, we find 3 rows. The first and second rows consider the update order σ 2 − S − m for S = 10−5 I and S = I, respectively. The third row shows coordinate updates for the order S − m − σ 2 . In this case, the initial value of S is irrelevant. The left plots show the estimated noise variance for both ZERO and PCA prior mean models, for five iterations. Iteration 0 corresponds to the σ 2 value at initialization, which is 0.5. The right plots show the optimization, for both models, of the variational mean and S − KZZ using images. When m and S − KZZ take values around 0, we expect a posterior collapse to the prior. That is, m = 0 and S = KZZ , which are the KLD minimizers. Fig. 8 shows that the model with the ZERO prior mean learns higher noise variances, and mostly collapses the variational parameters to the prior (the variational mean is slightly different from zero), while the PCA prior mean model does not. Total collapse is not produced in practice unless KXZ = 0, because σ 2 cannot grow to ∞. This supports our previous analysis, where we have shown that the predictive variances are bigger for the ZERO prior mean, and also the probability that KXZ = 0, effectively yielding higher noise variance parameter updates. Regarding the order of the updates, S − m − σ 2 consistently provides better convergence, as we expected. Importantly, the PCA prior mean model updates significantly reduce the noise variance parameter. In each model, both ordering updates converge to the same maximum of the ELBO. We also observe in Fig. 8 that smaller values of α yield smaller values of σ 2 at the first iteration (order σ 2 − S − m) since they result in a smaller initial predictive variance. These coordinate updates also show that an initialization of S = I in the output layer is more prone to the posterior collapse problem, as already discussed. Finally, an experiment with more inducing points supports these conclusions and is analyzed in appendix C.2. The coordinate updates provide an intuitive idea behind the local minima that we may find early at the training stage of our DGP. Since we have closed-form updates, we are removing any noise coming from stochastic gradients or from steps coming from different learning rates, which might help to jump out of the pathological local optima. Nevertheless, 43
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
2.0 1.5 1.0 0.5 0.0 0
1
2
3
Iterations
4
5
m
PCA Mean Zero Mean
ZERO MEAN PCA MEAN
PCA Mean Zero Mean
Noise Variance
2.5
1.0 0.5 0.0 −0.5 −1.0
S − KZZ
0.00 −0.25 −0.50 −0.75
0
1
2
3
Iterations
4
5
(a) Order of updates: σ 2 -S-m, with initialization S = 10−5 I, m = 0
2.0 1.5
m
PCA Mean Zero Mean
ZERO MEAN PCA MEAN
1.0 0.5 0.0 −0.5 −1.0
1.0 0.5 0.0 0
1
2
3
Iterations
4
5
PCA Mean Zero Mean
Noise Variance
2.5
S − KZZ
0.0 −0.2 −0.4 −0.6 −0.8
0
1
2
3
Iterations
4
5
2
(b) Order of updates: σ -S-m, with initialization S = I, m = 0
2.0 1.5
PCA Mean Zero Mean
ZERO MEAN PCA MEAN
1.0 0.5 0.0 −0.5 −1.0
1.0 0.5 0.0 0
1
2
3
Iterations
4
5
PCA Mean Zero Mean
Noise Variance
2.5
m
S − KZZ
0.0 −0.2 −0.4 −0.6 −0.8
0
1
2
3
Iterations
4
5
2
(c) Order of updates: S-m-σ . Initialization does not matter since coordinate updates on S do not depend on previous variational parameter values.
Figure 8: Each row shows coordinate updates, different initializations, and update order (indicated within each sub-caption). The left figures show the noise variance parameter. The right figures top show the evolution of the variational parameter m and right-bottom figures show the value of S − KZZ . This allows us to illustrate whether the parameters are close to a posterior collapse whenever m and S − KZZ are close to zero. Simulations are shown for both PCA and ZERO prior mean sparse GPs using 5 inducing points. 44
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
overall our analysis shows that early in the training, it is likely that a ZERO prior mean DGP will end up learning the targets as just noise and collapsing the variational parameters, unless previous layers move out quickly from providing a zero input vector signal at the last layer. Finally, to complete our analysis, we need to consider what happens when all the parameters get optimized. In this case, there are no exact coordinate updates for the inducing points nor the kernel parameters. These parameters get optimized by gradientbased coordinate updates. In some cases, even if no exact coordinate updates are given, it is easy to understand what happens to some of these parameters. For instance, we have previously shown in Eq. (47) that, when collapse happens, the likelihood noise variance is given by the sample variance and the kernel output scale variance σo . However, when σo gets optimized, this is not true. In this case, it is tempting to think that σo will grow 2 that attains arbitrarily large values and, without bound. This leads to a noise parameter σopt consequently, favors overestimating the variance and explaining the data as noise. However, in the whitening case, when collapse happens, S = I. Thus, attending to the objective in Eq. 33, the trace only depends on KXX for any point, and takes its smallest value by making σo = 0. This is because in the whitening case, σo does not influence the KLD. This is observed in one of our experiments. 4.4 Impact of the Parameterization on the Optimization Process In Sec. 4.2.5, we observed that the posterior collapse problem may happen independently of the parameterization used. Intuitively, however, the whitened parameterization should yield a more stable optimization procedure than the non-whitened parameterization. This verified empirically in Fig. 9. This figure shows, in the toy problem, the optimization of the ELBO and its decomposition in the ELL term and the KLD term, when the PCA prior mean function is used, for the whitened parameterization (solid orange line) and the non-whitened parameterization of GPFLOW (dashed blue line). Here, we initialize m = 0 and Sv = S = 10−5 I. In this section, we show that the reason behind these instabilities for the non-whitened parameterization is a poorly conditioned Hessian of the KLD term. Our analysis also shows that, in the non-whitened model, the reparameterization of the variational mean of GPFLOW does not provide better convergence guarantees. In general, the training instabilities of the non-whitened parameterization can make the model arrive at the local minimum in which the variational parameters are set equal to the prior, making it more likely for the posterior collapse problem. 4.4.1 The KLD at Initialization First, in the non-whitened representation, the KLD at initialization, for m = 0, is expected to be smaller for the GPFLOW reparameterization of the variational mean (see Table 1). −1 The reson is the impact of the quadratic term (m − µZ )T KZZ (m − µZ ), which makes the KLD greater for mean functions such as the PCA. Regarding the non-whitened vs. whitened parameterization, it will depend on the initial choice of the variational parameters. Initially, we may think the whitened representation will have smaller KLD on initialization, since the variational and prior mean are zero and the covariance matrix S = αI is diagonal at initialization. For α = 1, this is clearly true, and KLD = 0. However, a careful inspection of 45
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
ELBO
−2000 −3000 −4000 −5000 0
2000
4000
6000
8000
2500
0 −500 −1000 −1500 −2000
PCA W PCA NWR
−2500
10000
0
2000
Epoch
4000
6000
8000
KL Divergence
Expected Log Likelihood
0 −1000
2000 1500 1000 500
10000
0
2000
Epoch
4000
6000
8000
10000
8000
10000
Epoch
ELBO
−1000 −2000 −3000 −4000 −5000 0
2000
4000
6000
8000
10000
2500
0 −500 −1000 −1500 −2000
PCA W PCA NWR
−2500 0
2000
Epoch
4000
6000
8000
KL Divergence
Expected Log Likelihood
(a) Optimization curves for unnormalized data
2000 1500 1000 500
10000
Epoch
0
2000
4000
6000
Epoch
(b) Optimization curves for normalized data
Figure 9: Training curves reporting the ELBO, the ELL term, and the KLD term of a PCA prior mean DGP on the toy problem, for both the whitened and the non-whitened parameterization of GPFLOW. Figures show that the non-whitened parameterization is more unstable. the KLD at initialization reveals opposite observations for other values of α. In particular, −1 1 , log | KZZ | −M log α − M + αTr KZZ 2 1 KLDwhitened = [−M log α − M + M α] . 2
KLDnon-whitened =
(49) (50)
Ignoring common terms −M log α and −M , and expressing computations in terms of the eigenvalues of KZZ , λk , we have: ! M M X 1 X 1 KLDnon-whitened = log λk + α + const. , (51) 2 λk k=1
k=1
1 KLDwhitened = M α + const. 2
(52)
Thus, there is no reason to establish that any of the parameterizations will have a higher KLD at initialization. For small α, the KLD in the non-whitened case is dominated by the P term M k=1 log λk . This implies that only when the length-scale is small relative to the data, we will have KZZ ≈ I (for an output scale parameter of 1), and this term will vanish. If the matrix has a large number of large eigenvalues, we can expect this term P to be above 0, 1 and KLD to be greater. For moderately big α, the term is dominated by M k=1 /λk . In any case, any degenerated direction will incur a high KLD. However, for the usual length-scales used at initialization (from 0.1 to 2.0) and using the fact that smaller values for α are generally used, we observe that the KLD of the non-whitened parameterization will be smaller 46
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
at initialization, as reported in Fig. 9. A smaller initial KLD is expected to increase the probability of a posterior collapse problem. 4.4.2 Optimization Instabilities in the KLD Term Fig. 9 shows that the KLD has a great impact on the instabilities of the non-whitened parameterization of GPFLOW. Even when the data is normalized, the non-whitened parameterization presents unstable convergence. Therefore, we look at the Jacobian of the KLD, for the non-whitened parameterization of GPFLOW, reported in Table 1. At initialization, such a Jacobian is: ∂ KLD = 0, ∂m ∂ KLD 1 −1 = KZZ − S −1 , ∂S 2
∂ KLD = 0, ∂µZ ∂ KLD 1 −1 −1 −1 = KZZ − KZZ SKZZ . ∂KZZ 2
(53) (54)
By contrast, in the case of the whitened parameterization, the Jacobian is: ∂ KLD = 0, ∂mv ∂ KLD 1 −1 = I − Sv , ∂Sv 2
∂ KLD = 0, ∂µZ ∂ KLD = 00T . ∂KZZ
(55) (56)
These equations show that the bigger the KLD, (i.e., the bigger the difference between prior and variational covariances), the bigger the Jacobians are. This would imply using an adaptive learning rate algorithm to take more time to adjust and produce stable convergence. Moreover, since, in our toy problem, the whitened model has higher KLD, we should expect worse-behaved convergence. However, this is not the case. We believe this may be because in the whitened parameterization, the KLD does not depend on KZZ and hence the KLD Jacobian is zero w.r.t. the inducing points, which should lead to more stable gradient updates. One way to analyze the optimization stability of a loss function is via its Hessian, which measures the curvature of the loss in all parameter-space directions. If the Hessian is poorly conditioned, i.e., its eigenvalues differ widely across directions, first–order gradient optimizers struggle, and convergence slows down due to exceedingly small steps in some directions and excessively large steps in others. This produces updates that oscillate around the minima, which is the noisy behavior observed in the non-whitened reparameterized model. We have compiled the derivations needed to analyze the Hessian of both parameterizations in a separate work (Maroñas, 2025). This work shows that the Hessian is highly dependent on KZZ , in the case of the non-whitened parameterization, and that the reparameterization implemented by GPFLOW (see Sec. 3.3.3) does not provide any convergence improvement. Overall, this justifies why the non-whitened parameterization provides unstable convergence, and we show that these instabilities may result in poor solutions where collapse appears. However, we note that these observations open the door to a broader and more detailed theoretical analysis of optimization stability that falls outside the scope of this paper. A comprehensive treatment of these aspects will be presented in a separate work. 47
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
4.5 Multi-output Models, Other Likelihoods and Alternative Assumptions Regarding other likelihoods, the first natural question is what we can expect from a C multi-output DGP regression model with a Multivariate Gaussian likelihood. This likelihood has a noise covariance parameter Σ, which might not be diagonal. One can show that coordinate updates in this model can be written down in a similar way to the coordinate update of the 1-dimensional regression setting (see appendix A). Namely, N
Σopt =
T 1 X n n y − µnqf yn − µnqf + Kqf , N
(57)
n=1
where now yn is a C-dimensional target, µnqf is the C-dimensional predictive mean vector n is the C × C predictive covariances at xn . From this coordinate update, we can and Kqf n are expected to identify similar aleatoric and epistemic terms. The predictive covariances Kqf maintain the same pathological properties we have observed so far, both if the processes are independent or dependent (assuming we have initialized the mixing matrix to I). Specifically, n is obtained via similar expressions as those given in Eqs. (6) and Eqs. (17) using block Kqf (sometimes diagonal) matrices. Thus, the problem of the output layer receiving as input a vector of zeros is exactly the same. In general, any likelihood that has a noise parameter to model aleatoric uncertainty, such as the Student-t for robust regression or the Laplace, may be more likely to suffer from posterior collapse. This includes, obviously, heteroscedastic regression, where the model can learn high variances per datapoint. Other likelihoods for other tasks, such as classification (categorical or Bernoulli likelihoods) or regression over count data (Poisson Likelihood), might need more time to achieve convergence, but do not have a parameter to learn everything as noise. For example, in classification problems, aleatoric and epistemic uncertainty are obtained from the average entropy and mutual information after sampling from the posterior distribution, so uncertainty estimation depends purely on the latent functions (Depeweg et al., 2018). This means that posterior collapse is expected to be less relevant and is likely to depend only on the number of datapoints and processes. For a categorical predictive distribution, the training objective in a C classification problem can be written down (see, for instance, Maroñas and Hernández-Lobato (2023) for a formal derivation) as: S
ELBO =
N
1 XX S
s=1 n=1
−
C X c=1
Z
L L L L q(f1,s . . . fC,s ) log p(yn | Softmax(f1,s . . . fC,s )) c
c
KLD [q(u ) || p(u )] .
(58)
Now, if we consider the model at initialization, we know that the inputs to the last layer will be around 0, for the ZERO prior mean DGP. This implies that the ELL evaluates to log(1/C ) for each training sample at initialization. For large C this term becomes very small, the KLD is not going to dominate, and the model is forced to push the parameters to stop predicting latent functions around 0. For small C, it might be more beneficial to output constant predictive distributions and minimize the KLD term. Therefore, for classification problems, we might expect either posterior collapse or a substantial time to converge. The PCA prior mean DGP will not have this problem. 48
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
We have seen that most of the issues related to the posterior collapse problem come from the distance between 0 and the inducing points. This suggests that the ZERO prior mean DGP might benefit, at initialization, from using non-stationary kernels. However, previous work has demonstrated that non-stationary kernels can incur severe misspecification, and that well non-stationary inductive biases need to be encoded for successful performance. See the work of Maroñas and Hernández-Lobato (2023) for further discussion and comparison. Regarding the use of more Monte Carlo samples to estimate the ELL in the ELBO, we shall expect an exacerbated effect if more samples are used. We believe that fewer samples relate to higher variance of the ELBO gradients, which may be useful to escape bad local optima such as the one associated with the posterior collapse problem. On the other hand, as the number of Monte Carlo samples increases, the estimator of the ELL converges to its expected value. Since this expectation is taken under a distribution with zero mean, this convergence amplifies the degeneracy we have been illustrating, bringing the optimization process closer to the pathological regime associated with posterior collapse. Since in practice one sample is enough for the model to achieve good convergence, and due to computational cost, we have not carried out extra experiments exploring this direction. 4.6 Conclusions about how to Initialize DGPs In this section, we have consciously analyzed the posterior collapse problem under different parameterizations and initializations of DGPs. The takeaways from the analysis are mainly concerned with removing most of the noise from the optimization procedure. This is, in turn, achieved by initializing the likelihood noise parameter σ 2 to a small value, and the variational variance Sv = S = αI with α ≤ σo . Other valid options are S = αKZZ with α ≤ σo . In general, select S so that KZZ − S is positive definite since it implies KZZ ⪰ S. This will result in a predictive variance that does not go beyond the prior at initialization, resulting in a smaller trace term in the ELL. Setting Sv = S = αI, for α ≤ σo , always satisfies this requirement in both parameterizations. A careful look at some DGP implementations has revealed that these two options are the most common in many works under any of the three parameterizations. Moreover, more inducing points help to alleviate posterior collapse. Finally, if posterior collapse is still present, one can reduce the kernel scale parameter σo , which might result in smaller predictive variances at input locations far from the inducing points. Also, using a whitened parameterization results in better-behaved gradients, reducing the noise in the optimization and reducing the posterior collapse problem.
5 An Initialization to Alleviate Posterior Collapse In the previous section, we have described the optimization problems suffered by ZERO prior mean DGPs, and explained why the PCA prior mean function is able to reduce their impact. Since Bayesian theory dictates that the choice of the prior must be driven by modeling assumptions, and not by limitations of the optimization procedure, this section introduces an alternative way of initializing ZERO prior mean DGPs such that the possibility of a posterior collapse problem is reduced. The idea relies on initializing the hyper-parameters and the variational parameters such that the predictive distribution of the ZERO prior mean DGP mimics the PCA prior mean DGP at initialization. By doing this, we can avoid the local optimum in which the posterior 49
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
collapse occurs. We achieve this behavior by setting adequate initial values for the model and variational parameters so that the predictive mean µqf , µqfr , µqfw at initialization in the inner layers resembles that of a DGP with the PCA prior mean function. We present a general framework and impose restrictions to yield an efficient implementation. As we will show in the experiments section, this approach allows effective training of DGPs with ZERO prior mean functions. However, the proposed solution should be valid for a DGP with any prior mean function. This general case is presented in Sec. 5.3, but has not been validated experimentally. In this section, we will consider three parameterizations: whitened (W), standard nonwhitened (NW), and non-whitened with the reparameterization of the variational mean of GPFLOW (NWR). Besides this, we will consider two mean functions PCA and ZERO. Therefore, we define the following acronyms to refer to each of these models throughout the rest of the work, see Table 2. Table 2: Summary of the acronyms used to refer to the different models considered in this work, attending to their mean function and the parameterization of the inducing points. Parameterization Whitened Non-whitened Non-whitened reparameterized Whitened Non-whitened Non-whitened reparameterized
Mean Function
Acronym
ZERO
ZERO-W
ZERO
ZERO-NW
ZERO
ZERO-NWR
PCA
PCA-W
PCA
PCA-NW
PCA
PCA-NWR
5.1 The Predictive Mean at Initialization Our goal is to mimic the PCA prior mean DGP at initialization using a ZERO prior mean function. Therefore, we start by inspecting the predictive mean of this model, for the initialization outlined in Sec. 3.5 (i.e., m = mv = mr = 0). To lighten notation, we will be presenting our approach as if we had a single GP per layer. Following Eqs. (6), (17), (21), the predictive mean of a PCA prior mean DGPs at initialization take values: −1 µqf = Xs − KXs Z KZZ Z,
µqfw = Xs ,
µqfr = Xs ,
(59)
for each parameterization, W, NW and NWR, respectively, where Xs is the set of points in which we are evaluating our posterior. In general, Xs will represent the input points where we want our ZERO prior mean DGP to mimic the PCA prior mean DGP at initialization. This may include, e.g., the train set, a mini-batch of points, or the whole input domain X . Ideally, we would like to be able to mimic the PCA prior mean DGP at initialization on X , yet we will show in Sec. 5.2 that this is not possible. Remark 1 Following Sec. 3.4, for a general DGP with an arbitrary number of GPs per layer, Eq. (59) is extended by noting that Xs represents the d−th coordinate of Xs in the evaluation of the predictive mean of the PCA prior mean DGP. Also, when the number of GPs at the level of the hierarchy is smaller than at the previous level, the predictive mean 50
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
function at that layer is given by a PCA projection of the inputs Xs . Thus, we just need to change Xs by the corresponding PCA projection. When the model is a ZERO prior mean DGP, the predictive mean is: T −1 −1 µqf = KXZ KZZ m, µqfw = KXZ L−1 mv , µqfr = KXZ KZZ mr , ZZ
(60)
for each parameterization. Clearly, at initialization, these predictive means evaluate to 0. Thus, to select initial parameters so that, at initialization, a ZERO prior mean DGP has the same predictive mean at initialization as that of a PCA prior mean DGP, we need to match Eq. (59) with Eq. (60), at each inner layer. This involves solving a system of linear equations. The learnable parameters are m, mv , mr , Z, X, Xs and ν . 5.2 Solving the System of Equations The three systems of equations that result from the above idea are: −1 −1 KXZ KZZ m = Xs − KXs Z KZZ Z, −1 T KXZ LZZ mv = Xs , −1 KXZ KZZ mr = Xs ,
(NW) (W) (NWR)
for each parameterization. We note that we do not require the set of points Xs to be the same as X in order to define the system, but it is required that dim(Xs ) = dim(X), where dim(·) represents the dimensions (rows and columns) of its argument, for the three systems to be well-defined. Since it is counterintuitive to require the ZERO prior mean DGP to mimic the PCA prior mean DGP at different points in the domain, we simply set X = Xs . Ideally, we would like to mimic the PCA prior mean DGP at the full domain, i.e. Xs = X . However, this implies having an infinite number of equations. With this in mind, under some assumptions, these systems can be turned into linear systems that can be solved without requiring numerical gradient-based optimization. However, this is not a restriction but a decision we have made. Thus, there is freedom to go through other approaches that might be interesting for other reasons (for example, incrementing the support Xs ), and still yield the desired initialization. These additional details are presented in Appendix B. The set of assumptions we have made includes fixing the values of the kernel hyperparameters ν , Xs , and Z and solving for m. In this case, the optimal values for the variational means, which would resemble the predictions of the PCA prior mean DGP at initialization, at some selected Xs and Z, are given by −1 −1 −1 m = KXs Z KZZ Xs − KXs Z KZZ Z , (NW) h −1 T i−1 mv = KXs Z LZZ Xs , (W) −1 −1 mr = KXs Z KZZ Xs , (NWR) T −1 for the three parameterizations considered, assuming KXs Z L−1 and KXs Z KZZ are ZZ invertible. If we additionally impose that the inducing points and Xs are exactly the same 51
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
points, i.e. Xs = Z := Xxz , for some set of inputs points Xxz , the solutions become: m = 0,
mv = L−1 Xxz Xxz Xxz ,
mr = Xxz .
(61)
This last assumption is used to remove many of the operations with cubic cost in the solutions. Remark 2 For the NW model, the optimal value for a ZERO prior mean DGP to mimic a PCA prior mean DGP at initialization is m = 0. However, this is exactly the initialization value we want to avoid. In other words, with the assumptions taken up to this point, the ZERO prior mean DGP with the standard initialization already mimics the PCA prior mean DGP at the inducing points values. However, in this work, we consider the non-whitened parameterization implemented by GPFLOW (NWR) and also give the necessary steps to develop initializations for the NW model implemented by GPYTORCH and GPJAX. These can be found in Appendix B. 5.3 Arbitrary Prior Mean Functions The presented derivations consider a ZERO prior mean DGP that mimics a PCA prior mean DGP. However, our results can easily be extended for DGPs with arbitrary mean functions to mimic a PCA prior mean DGP. In particular, for an arbitrary prior mean function µ, the means of the predictive distribution posteriors are: −1 µqf = KXZ KZZ (m − µZ ) + µX , −1 T µqfw = KXZ LZZ mv + µX ,
(62)
−1 µqfr = KXZ KZZ mr + µ X ,
(64)
(63)
for each parameterization, respectively. Therefore, the systems of equations to solve are: −1 −1 KXZ KZZ (m − µZ ) + µX = Xs − KXs Z KZZ Z, −1 T KXZ LZZ mv + µX = Xs , −1 KXZ KZZ mr + µX = Xs .
(65) (66) (67)
Given the same assumptions we have taken so far (fix values for ν , Xs , Z and setting Xs = Z = Xxz ), the solutions are: m = 0,
mv = L−1 Xxz Xxz (Xxz − µXxz ) ,
mr = Xxz − µXxz .
(68)
Remark 3 As a sanity check, note that if the prior PCA mean function is used, we obtain m = mv = m = 0, as expected, since µXxz = Xxz . Since both models are the same, they should have exactly the same parameters at initialization if one wants to mimic the other. 5.4 Fixing Xxz and Initializing the Output Layer To Predict the Targets Above, we assumed that Xs = Z := Xxz , for some set of input points Xxz , to reduce the computational cost of the approach. There are two ways in which we can initialize Xxz : 52
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1. Initialize the inducing points and let Xxz = Z := Xz . In this way, our initialization will resemble the predictive mean of the PCA prior mean DGP at the inducing locations Z, which are usually initialized via K-means. The inducing points obtained will be representative of the training data, and will guarantee that we match the initial predictions on a representative set. 2. Select Xxz = Z ⊂ X := Xx to be representative of the training data. This means that the inducing points are a subset of the training inputs. Note that the solutions to the initialization problem for both options are given by: m = 0,
mv = L−1 ZZ Z ,
mr = Z ,
(69)
where Z is obtained either through k-means or as a subset of the training data. The second option, however, has the advantage that it allows us to initialize the output GP of the DGP so that it predicts values that are closer to the corresponding targets associated with each selected point. Usually, a ZERO prior mean function is used at the output layer, with m = mr = mv = 0. However, by using a similar reasoning to that used in Sec. 5.1, we can initialize the output GP so that it has a predictive mean equal to the corresponding target values at some of the training points. This is expected to start the optimization process at a good initial solution and to help avoid the posterior collapse problem. We can think of this approach as mimicking a DGP model whose predictive mean function takes as input X and outputs the corresponding y. Initializing the output layer like this is similar to what we have described so far. One only has to replace Xs by the corresponding ys . This results in the following systems of equations: −1 −1 KXZ KZZ m = ys − KXs Z KZZ yz , −1 T KXZ LZZ mv = ys ,
(70) (71)
−1 KXZ KZZ mr = ys ,
(72)
for each parameterization. Note that we require the targets, ys , associated with Xs . This is the reason why selecting Z ⊂ X (option 2 above) becomes useful. By selecting the inducing points from the training data, we have the associated targets, something not possible in the first option. Using the same assumptions as before, namely, the inducing points and training points are the same locations, the solutions for this particular case are: m = 0,
mv = L−1 Xx Xx yx ,
mr = yx .
(73)
The fact that the inducing and training points are the same, and that the inducing points are selected from the training set, implies that ys = yz = yx and Xxz = Xx . 5.5 Initialization of the Inducing Points As mentioned earlier, our initialization requires a subset from the training data that is as representative as possible of the full support over where we are making predictions. Selecting a random subset of M instances, X, from the training set may concentrate the inducing points in a specific region of the domain. Thus, to alleviate this problem, we propose the following initialization of the inducing points Z: 53
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1. We compute M centroids C = {c1 , · · · , cM } from the training data, Xtrain , using k-means. 2. We select the inducing locations zj from Xtrain using the cosine distance to each centroid: zj = argmin d(cj , xj ) xj ∈Xtrain
= argmin cj · xj /(∥cj ∥∥xj ∥), xj ∈Xtrain
∀cj ∈ C .
(74)
When a training point is selected as the closest one to a centroid, we remove it from the training set so that it is not selected again. In this way, all the inducing points differ from each other, but they are the training instances closest to the generated centroids. Fig. 10 shows the predictive distribution of a DGP with 2-layers with the ZERO prior mean. The whitened parameterization is used. In all cases, the variational covariances are initialized at 10−5 I. We consider three scenarios: (a) when the variational mean is initialized to zero; (b) when the proposed initialization is used for the variational mean at each layer, but the inducing points are randomly chosen from the training set; and (c) when the proposed initialization is used, but the inducing points are selected using the algorithm described above. We observe that the last scenario is the one producing a more accurate initial predictive distribution that best matches the observed data.
(a) Standard initialization: ml = 0, Sl = 10−5 I.
(b) Proposed initialization using random Z ⊂ X.
(c) Proposed initialization selecting Z ⊂ X using Eq. (74).
Figure 10: Initial predictions of whitened 2-layer DGPs using the ZERO prior mean. Each sub-figure considers a different initialization of the variational mean of the output layer. In the following sections, we will label as ZERO-Points-M0-NWR and ZERO-Points-M0-W the models that use m = 0 in the output layer, but use the proposed initialization of the inner layers given in Eq. (69), for the non-whitened and the whitened parameterizations, respectively. Similarly, we will label with ZERO-Points-MY-NWR and ZERO-Points-MY-W to the model that initializes the predictive mean of the output layer so that the corresponding targets are predicted well, i.e. using Eq. (69) in the inner layers and Eq. (73) in the output layer. In all cases, we consider the non-whitened parameterization of GPFLOW. 5.6 Selecting the Initial Kernel Length-scale ℓ The proposed initialization for the variational mean aims at predicting the observed target values. Therefore, it becomes relevant to carefully set the initial kernel length-scale, ℓ, for 54
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Algorithm 1 Lengthscale ℓ Selection via RMSE Minimization Require: Train Dataset D, number of trials Ntrials , number of steps S, step size δ 1: pairs← [ ] 2: for s = 0 to Ntrials − 1 do 3: Dst , Dsv , ← random_partition(D) ▷ Randomly select 10% of D as validation 4: d ← median_pairwise_distance(Dst ) 5: G ← {d ± i · δ | i = 1, · · · , S} ∪ {d} ▷ Generate candidate lengthscales 6: G ← {ℓ ∈ G | ℓ > 0} 7: for each l ∈ G do 8: Mℓ ← initialize_model(Dst , ℓ) ▷ Init model using lengthscale ℓ 9: rℓ ← evaluate_rmse(Mℓ , Dsv ) ▷ Compute RMSE on Ds 10: Append (ℓ, rℓ ) to pairs 11: end for 12: end for 13: return ℓ⋆ with lowest rℓ⋆ in pairs
the particular given dataset in the case of using the initialization in the output layer (MY models). In particular, in such a case, the quality of the initial predictive distribution may strongly depend on the length-scale, as illustrated by Fig. 11. A heuristic used to initialize ℓ may guarantee that the initial predictive distribution on the observed data is good. With this goal, given a training set, we split it into training and validation, and compute the median Euclidean distance, md , between the datapoints in the training set. Then, we create a grid G of evenly-spaced length-scales with the middle point equal to md . For each length-scale in G, we initialize the model using that length-scale and evaluate RMSE of the corresponding initial distribution on the validation set. We finally select the length-scale ℓ⋆ ∈ G with the lowest RMSE value. In practice, however, we have found that on some datasets, md may vary a lot depending on the train-validation split. To avoid this issue, we repeat the selection in multiple train-validation splits and select the length-scale that provided the lowest RMSE across different splits. The full algorithm is described in Alg. 1. Note that the proposed heuristic focuses on selecting ℓ with the lowest RMSE, without considering the resulting KLD at initialization. We noted that adjusting the length-scale to ℓ = 0.01
ℓ = 0.1005
ℓ = 1.0
Figure 11: Initial predictive distributions of the ZERO-Points-MY-W model using different initial length-scale ℓ values. 55
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
improve RMSE hardness the KLD. This comes from replacing mr , mv by L−1 Xx ,Xx yx and yx in −1 T T the quadratic terms mr KZZ mr and mv mv . Fig. 22 in the appendix provides a comparison between the KLD value and the length-scale. Future work may consider improving this selection strategy by designing a new heuristic that also takes into account the initial KLD, which may increase, depending on the chosen initial length-scale. 5.7 Summary of the Proposed Initialization Strategy As a summary, the set of steps taken to yield the initial variational parameters are: 1. Find Z via the procedure described by Eq. (74) for the MY models. The M0 models can also be initialized through k-means. 2. Set Xs = Z = Xx . 3. In the case of initializing the output layer to predict the targets, heuristically determine the length-scale. Otherwise, use the default length-scale. 4. Set the optimal variational mean by Eq. (61). The proposed initialization has the limitation of being able to mimic the predictive mean of the PCA prior mean DGP, at the inner layers, only at a subset of the training points, given by the inducing points. Increasing the support of the initialization can be done with the procedure described in Appendix B, but results in more complex systems of equations. Importantly, our initialization should only improve the results substantially, in those ZERO mean models that suffer from posterior collapse. Our initialization does not change the statistical model, nor the approximate inference distribution. This implies that the space of possible solutions is the same across all the ZERO models. Thus, we shall only expect either solving the posterior collapse problem (through either M0 or MY initializations) or, given the same number of training iterations, faster convergence to a solution of the optimization problem (MY initialization). As a consequence, the MY initialized model may surpass the rest of the models. 5.8 Dependent and Non-Gaussian Processes at Each Layer This work considers DGP where the GPs at each layer are independent a priori. However, the proposed initialization strategy can work as well with models in which the GPs at each layer are dependent or where the latent processes are transformed in some way, yielding non-Gaussian processes as building blocks. Such modifications of the standard DGP can be achieved either by using a mixing matrix (Álvarez et al., 2012) or by considering warping functions (Maroñas et al., 2021; Maroñas and Hernández-Lobato, 2023; Sáez-Maldonado et al., 2023; Lin et al., 2023). Since these models use GPs as building blocks, and act directly on the samples or distributions over the function evaluations, our initialization and analysis remain valid. We just need to proceed in the same way to mimic the GP at initialization, and then possibly initialize the next elements of the modeling processes in a way that at initialization the behavior is the same as the standard DGP. For this, we just need the mixing matrix to be initialized to the identity matrix, and the warping functions to be an identity function. 56
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
6 Experiments In this section, we will detail the experiments performed to validate our theoretical observations and to test the proposed initialization. We have considered two scenarios for this validation: the first one is the toy dataset used in Rudner et al. (2020), which has already been used to analyze the optimization issues of DGPs in Sec. 4. The experiments performed on this dataset are useful to exhibit the pathologies of DGPs supported by the visualizations of the predictive distribution. The second scenario involves real datasets from the UCI repository (Dua and Graff, 2019). In these experiments, we will check that the conclusions obtained in the toy dataset are extrapolated to real-world datasets. The baseline models considered are DGPs with the ZERO prior mean, and the whitened and non-whitened parameterization of GPFLOW. We refer to them as ZERO-W, ZERONWR. Similarly, we consider equivalent DGP models with the PCA prior mean. We refer to them as PCA-W and PCA-NWR. See Table 2. Besides this, we also consider our proposal initializations, labeled as ZERO-Points-M0-W, ZERO-Points-M0-NWR, ZERO-Points-MY-W , ZERO-Points-MY-NWR, where M0 indicates that the output layer’s variational mean is L mL r = m = 0 at initialization, and MY indicates that we use our proposed initialization L for m in the output layer to predict the observed targets. We note that the standard nonwhitened parameterization is not implemented in GPFLOW. Therefore, we do not evaluate it. The covariance function is, for all the models, parameterized by the RBF kernel defined in Eq. (27). 6.1 Toy Dataset As seen in Sec. 4, the toy dataset of Rudner et al. (2020) describes a 1-dimensional regression problem with a multi-step function. This problem can be hard for DGPs since they need to learn to use a very small length-scale in order to model the non-smooth steps of the target function. We utilize this dataset to analyze the optimization issues that arise within the DSVI algorithm and different initialization parameters. We divide this section into two fundamental parts. The first aims at supporting the arguments provided in Sec. 4 for the posterior collapse problem under the DSVI algorithm. The second one shows a comparison of all the models in this dataset. 6.1.1 An Illustration of the Posterior Collapse Problems of DSVI We begin by validating the statements made in Sec. 4.2.6 and Sec. 4.2.8. We focus on the case of ZERO-W and PCA-W, standard DGP models, as the whitened parameterization is expected to facilitate the training of the models. As a common setup, in this section, we train all the models for 10.000 epochs using Adam and a learning rate of 10−3 . We use 2, 3 and 5 layer DGPs. We analyze the effect of the different initializations for S and its relation with the number of inducing points. Sl will denote the variational covariance of the inner layers, and SL of the output layer. • Sl = 10−5 I and SL = 10−5 I: Fig. 12 shows the predictive distribution of the ZERO-W and PCA-W models for each configuration. We observe that as the number of inducing points decreases and the number of layers increases, the posterior collapse problem of the ZERO-W model becomes more evident. Specifically, 5 inducing points show 57
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
total collapse with 3 and 5 layers. With 20 inducing points, the model learns a poor predictive distribution when the DGP has 3 layers, and results in posterior collapse when 5 layers are considered. This effect is reduced when the number of inducing points is increased to 100. In this case, there is no posterior collapse, but the predictive distribution is very poor in the deeper model. The PCA-W model is able to escape posterior collapse. These results corroborate the ideas presented in Sec. 4.2.8, indicating that more inducing points should alleviate the posterior collapse problem. Furthermore, they also illustrate that the ZERO mean prior DGP is more likely to have posterior collapse than the PCA prior mean DGP. • Sl = 10−5 I and SL = I: Fig. 13 shows the results obtained when the output layer covariances are initialized to those of the prior. We observe that, in this case, the ZERO-W model with 5 layers and 100 inducing points suffers from posterior collapse. When SL = 10−5 I, the model does not collapse, but the predictive distribution is poor. This confirms that setting SL = I increases the probability of a posterior collapse since 4 out of 9 models have collapsed. Again, we observe that using a higher number of inducing points leads to better predictive distributions and reduces the probability of a posterior collapse. Specifically, the predictive distribution of the ZERO-W DGP with 3 layers improves substantially when 100 inducing points are considered, similar to the previous setting considered in Fig. 12. Again, the PCA-W model does not result in posterior collapse, unlike the ZERO-W model. However, in this model we can observe a poorer fit with 5 inducing points, especially with 3 layers. This confirms that, in general, setting SL = I results in suboptimal fitting compared to SL = 10−5 I. • Sl = I and SL = 10−5 I: Fig. 14 shows the predictive distributions in this case, where the inner layer variance is higher. Again, we observe that the ZERO-W model suffers from posterior collapse in four out of the nine cases. This confirms that setting Sl = I in the inner layers increases the probability of posterior collapse. Compared to the previous setting, we observe that the ZERO prior mean model with 3 layers and 5 inducing points slightly escapes from collapse. As we discussed at the end of Sec. 4.2.6, this might be because, since the variability in the samples from the inner layer is higher, the chance of samples being far from 0 increases, which avoids the pathological behavior of the ZERO DGP at initialization. Nevertheless, the solution is poor. In Appendix C we provide further insights on this configuration. Furthermore, in this setting, the predictive distribution of the PCA-W model is also slightly worse than when Sl = 10−5 I for 20 inducing points. These results confirm that one should avoid setting high initial uncertainties in the inner layers of the DGP. • Sl = I and SL = I: For completeness, the pictures from this configuration are provided in Appendix C.3. Overall, this subsection corroborates our analysis carried out in Sec. 4, confirming that using smaller initial variational variances S and a bigger number of inducing points favors avoiding posterior collapse. Moreover, injected noise in the optimization results in suboptimal fitted models. In this case, the injected noise came either from a bigger initial Sl or a smaller number of inducing points, which makes the inner layer predictive variance take the value of the kernel output scale parameter σo at locations far from the inducing points. 58
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
3 Layers
5 Layers
100 Inducing
20 Inducing
5 Inducing
2 Layers
100 Inducing
20 Inducing
5 Inducing
ZERO-W
PCA-W
Figure 12: Predictive distribution (mean and standard deviation) of ZERO-W (top 3 rows) and PCA-W (bottom 3 rows) DGP models initialized with Sv = 10−5 in all layers. Blue points show training data. Columns correspond to different number of layers. Rows correspond to different number of inducing points: 5, 20 and 100 inducing points. 59
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
3 Layers
5 Layers
100 Inducing
20 Inducing
5 Inducing
2 Layers
100 Inducing
20 Inducing
5 Inducing
ZERO-W
PCA-W
Figure 13: Predictive distribution (mean and standard deviation) of ZERO-W (top 3 rows) and PCA-W (bottom 3 rows) DGP models initialized with Slv = 10−5 in the hidden layers and SL v = I in the output layer. Blue points show training data. Higher variance in the output layer results in poorer solutions. Columns correspond to different numbers of layers. Rows correspond to different numbers of inducing points: 5, 20, and 100 inducing points. 60
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
3 Layers
5 Layers
100 Inducing
20 Inducing
5 Inducing
2 Layers
100 Inducing
20 Inducing
5 Inducing
ZERO-W
PCA-W
Figure 14: Predictive distribution of both ZERO-W prior mean DGP (top 3 rows) and the PCA-W prior mean DGP (bottom 3 rows) initialized with Slv = I in the inner layers and −5 SL v = 10 I in the output layer. Blue points show training data. Using a high variance on the inner layers leads to a noisy optimization which yields poor solutions. 61
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
6.1.2 Evaluating a Very Deep DGP Model The previous results show that the PCA-W model is robust to posterior collapse in the evaluated configurations. However, we also observed that the models evaluated are more likely to suffer from posterior collapse when the number of layers increases and the number of inducing points is reduced. Therefore, in this section we train a 10-layer PCA-W DGP −5 model with 20 inducing points on the toy problem using Slv = SL v = 10 I. The predictive distribution obtained is shown in the top-left corner of Fig. 15. We observe that the model is close to posterior collapse. This is supported by the zero KLD shown in some of its hidden layers during the last training epochs, as displayed in the top-middle sub-figure. Moreover, the noise variance, displayed in the top-right sub-figure, is also increasing with each training epoch, suggesting that the model is simply explaining the observed data as noise. We hypothesize that the observed behavior may be due to a big initial prior variance, as specified by the kernel output scale parameter, i.e., σo in Eq. (27). The bottom row of Fig. 15 shows that when this initial parameter is reduced, the model does not suffer posterior collapse and learns to perfectly fit the observed data. Since there are few inducing points in this setting, far from them, the model will output the prior variance, which is big when σo = 1.0. This introduces noise in the optimization procedure, resulting in suboptimal fitting. By reducing the kernel output scale, we reduce this noise, which results in more effective training of the model. Summing up, these experiments show that the PCA prior mean function with a stable optimization procedure coming from the whitened parameterization might not necessarily avoid the posterior collapse problem, nor the non-injective pathologies (Duvenaud et al., 2014), and also that setting smaller initial prior variances may be beneficial to avoid the posterior collapse problem. This further confirms the result from the previous section, where injected noise in the optimization procedure drives the model to collapse. 6.1.3 Results for Each DGP Model and Initialization Considered We now proceed with a detailed performance comparison of all models in this dataset. The results obtained here will determine the model configurations that are run in the real-world datasets. For this comparison, we use 5 layers in all the experiments and train the models for 10.000 epochs, which is enough to achieve convergence of the ELBO. We use 100 inducing points. We initialize the inducing locations Z using K-means for the ZERO and PCA prior mean function models, and we use Z = Xx for the M0 and MY models. This means that the inducing points between the ZERO and PCA models and the M0 and MY differ. We use S = 10−5 I in all layers since attending to our first experiments, this configuration is more likely to produce the best results. Here, we consider the non-whitened parameterization of GPFLOW, i.e., NWR, and the whitened parameterization, i.e., W. First, we compare M0 and MY at initialization. We illustrate the initial predictive distribution of the DGP models using the standard value for the variational mean in the last layer, mL = 0, and our proposed initialization MY in Fig. 16. The initial predictive distribution of all the models using mL = 0 is the same, and is displayed in the left sub-figure. In the case of the MY models, we observe small differences in the W and NWR initializations. However, both are able to explain in a reasonable way the observed targets (middle and right sub-figures). Table 3 shows the initial ELL, KLD and RMSE of each model and initialization. 62
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Predictive distribution
Layer-wise KLD
Likelihood Variance
60
KL Divergence
σ=1
Layer 1 Layer 2 Layer 3 Layer 4 Layer 5 Layer 6 Layer 7 Layer 8 Layer 9 Layer 10
40 30 20 10
Likelihood Variance
0.8
50
0.6
0.4
0.2
0 0
9991 9992 9993 9994 9995 9996 9997 9998 9999 10000
2000
4000
6000
8000
10000
6000
8000
10000
Epoch
Epoch
30
KL Divergence
σ = 0.1
20 15 10
Likelihood Variance
0.08 Layer 1 Layer 2 Layer 3 Layer 4 Layer 5 Layer 6 Layer 7 Layer 8 Layer 9 Layer 10
25
0.06 0.04 0.02 0.00
9991 9992 9993 9994 9995 9996 9997 9998 9999 10000
Epoch
0
2000
4000
Epoch
10 layer PCA-W Figure 15: (top-row) Predictive distribution, layer-wise KLD and noise variance during training of the 10-layer PCA-W model with 20 inducing points. The kernel output scale parameter σ is initialized to 1.0. (bottom-row) Same information, but when σ is initialized to 0.1. We observe that σ = 0.1 results in a better model. By contrast, σ = 1 results in a model with near posterior collapse. As expected, all the models with a zero variational parameter in the output layer have the same ELL and RMSE since they output the same predictive distribution at initialization (the kernel hyper-parameters are the same). The models with the proposed initialization, i.e., ZERO-Points-MY-W and ZERO-Points-MY-NWR, obtain much better initial ELL and RMSE values. The reason is that the output layer aims to predict the targets at the inducing locations Z, which were selected as a subset of X. The initial KLD is, however, different. The big differences appear in the MY models, in which the variational mean is strongly modified in the last layer to predict the targets at the inducing points Z. This results in a big difference between the variational distribution and the prior, which yields a high KLD. In the case of the M0 models, there is an increment in the KLD due to the modified variational distribution in the inner layers. The large initial KLD value of the MY models results in a significantly larger ELBO, which can affect its optimization. We now analyze the final predictive distribution of each method after training. For this, in this toy problem, we use the learning rate λ = 10−4 . Table 4 shows the final results obtained by each method in terms of the RMSE, the KLD, and the estimated noise variance. We observe that a non-whitened model, i.e., ZERO-NWR, is the worst performing one in terms of RMSE, exhibiting a much lower KLD and a higher noise variance than the rest of the models. These are indicators of posterior collapse. The KLD does not go to zero due to noise in the optimization process. NWR models have, in general, worse RMSE, higher noise variances, and smaller KLD than their corresponding whitened versions. This indicates that the NWR parameterization results in models that tend to present posterior collapse, due to the noise injected coming from the parameterization. We observe how the ZERO-W model also provides suboptimal results with considerably worse results in terms of RMSE 63
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
ZERO output mean, mL = 0
ZERO-Points-MY-W
ZERO-Points-MY-NWR
Figure 16: (left) Initial predictive distributions of all the models that have mL = 0 (ZERO-W, ZERO-NWR, PCA-W, PCA-NWR, ZERO-Points-M0-W, ZERO-Points-M0-NWR), and for the proposed initialization in both W (middle) and NWR cases (right). The proposed initialization shows a very accurate initial predictive distribution. Table 3: Initial values for RMSE, ELL and KLD for each model on the toy dataset. Values are reported for the optimal length-scale for MY models (obtained using Alg. 1), i.e., ℓ = 0.1, and for the standard initial value ℓ = 2.0 for the other models.
ZERO-W ZERO-NWR PCA-W PCA-NWR ZERO-Points-M0-W ZERO-Points-M0-NWR ZERO-Points-MY-W ZERO-Points-MY-NWR
ℓ = 2.0
ℓ = 0.1
Model Name ELL
KLD
RMSE
ELL
KLD
RMSE
-2816.35 -2816.34 -2816.40 -2817.12 -2816.35 -2816.50 79.12 83.15
2628.23 1763.54 2628.23 1763.54 2665.76 1802.41 772874.04 772010.58
1.00 1.00 1.00 1.00 1.00 1.00 0.16 0.16
-2816.35 -2816.32 -2816.35 -2816.32 -2816.35 -2816.32 -339.79 -339.21
2628.23 1691.91 2628.23 1691.91 2672.95 1736.62 8797481.98 8796545.65
0.99 0.99 0.99 0.99 0.99 0.99 0.41 0.41
and a higher noise variance. The ZERO-NWR model presents nearly a collapsed distribution. Importantly, our proposed initialization obtains the best RMSE. However, we observe that the MY initialization typically results in a larger KLD term. Table 4: Test metrics by all the models on the toy dataset using λ = 10−4 . Model ZERO-W ZERO-NWR PCA-W PCA-NWR ZERO-Points-M0-W ZERO-Points-M0-NWR ZERO-Points-MY-W ZERO-Points-MY-NWR
KLD
241.4252 39.7501 259.3910 229.2737 281.4474 228.2416 316030.0797 483.9357 64
Noise Var. 0.1098 0.1205 0.0329 0.0344 0.0340 0.0331 0.0190 0.0199
RMSE
0.4591 0.9896 0.0321 0.0423 0.0294 0.0406 0.0055 0.0166
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Fig. 17 shows the predictive distribution of all the models, after training. We observe that, in general, standard ZERO prior mean DGPs with the non-whitened parameterizations suffer from a posterior collapse (bottom-left figure). Furthermore, a whitened parameterization alleviates such a problem, but results in poor predictive distributions (top-left figure). That is, the model tends to predict mostly noise. The PCA prior mean function is very effective in alleviating these problems, resulting in better predictive distributions. Similarly, our proposed initialization is also successful in avoiding the posterior collapse problem for all the evaluated configurations. It generates predictive distributions that look very similar to those of the PCA prior mean DGPs. In Appendix C.4, we extend with experiments using a higher learning rate. ZERO-W
PCA-W
ZERO-Points-M0-W
ZERO-Points-MY-W
ZERO-NWR
PCA-NWR
ZERO-Points-M0-NWR
ZERO-Points-MY-NWR
Figure 17: Predictive distributions of all trained models with λ = 10−4 . From the experiments carried out so far, we can conclude that: • In this dataset, a ZERO prior mean function with the standard initialization often leads to a posterior collapse problem or a bad predictive distribution, independently of the parameterization. • A whitened parameterization is to be preferred over a non-whitened one. Whitening often reduces the posterior collapse problem or improves the predictive distribution. • The proposed initialization of ml in the inner layers (that is, the M0 model), offers a similar behavior to the PCA prior mean DGP model and alleviates the posterior collapse problem, independently of the parameterization considered. • Using the proposed initialization in the last layer of the model (i.e., the MY model) results in slightly more accurate predictive distributions. However, it may induce high KLD values before and after training, which may complicate ELBO optimization. 6.2 Experiments on the UCI Datasets In this section, following observations made in the toy problem experiments, we train several DGP models on 8 real-world datasets from the UCI repository (Dua and Graff, 2019). These 65
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
are Boston, Concrete, Energy, Kin8nm, Power, Protein, Redwine, and Yacht. These datasets have been used in other works involving DGP (Salimbeni and Deisenroth, 2017). The goal of these experiments is two-fold. First, comparing the performance of the proposed initializations M0 and MY of the ZERO prior mean DGP model with the standard initializations used in the ZERO and PCA prior mean DGPs. Second, studying the posterior collapse problem and optimization stability in real-world datasets. In these experiments, we consider the experimental setup considered in Salimbeni and Deisenroth (2017). For optimization, we use Adam (Kingma and Ba, 2015). We train each model for 20, 000 iterations using a learning rate of λ = 0.01 and a batch size of at most 10, 000 instances. Each model is evaluated on 20 different train/test splits, and we report average results across splits. We train 2, 3, 4, and 5-layer DGP models, using 100 inducing points and min{D, 30} different GPs per layer, where D is the dimensionality of the data. The kernels are shared between the GPs in a layer, and the kernel parameters are initialized as ℓ = σ = 2.0. The noise variance is initialized to 0.05. For the PCA and ZERO mean prior models, and also for the M0 model, inducing points are initialized with K-means in the real-world experiments, which means Z = Xz . With this, the differences with respect to the PCA and ZERO prior mean DGP models just rely on the initialization of mr , mv . The MY model can only use Z = Xx , which means that inducing points also differ from those used by PCA and ZERO prior mean DGP since now the inducing points must be selected from the training data using the proposal in Sec. 5.5. For every model, we initialize the variational covariance to S = 10−5 I in all the layers. We estimate the test metrics using S = 100 Monte Carlo samples. In this section, we will show the results in terms of the test log-likelihood. The results in terms of RMSE draw the same conclusions, and we place them in the Appendix C.5.1 for brevity. Although the NWR parameterization showed poor performance in the toy example, we include results for both the whitened and non-whitened models for two main reasons. First, to the best of our knowledge, there has been no previous report in the literature comparing the performance of both parameterizations. Second, we have observed that, despite exhibiting greater optimization instability, the NWR parameterization achieves better results on several datasets. A detailed comparison between the two parameterizations is provided in Appendix C.5.2. 6.2.1 Main Results For Both Parameterizations Fig. 18a and Fig. 18b show the test log-likelihood obtained by each whitened and nonwhitened model, respectively. The figure shows that this metric is, in most of the datasets, very similar for all the models when we take into account the standard deviation of the results. In general, there are no significant differences when varying the number of layers. Importantly, the models using the proposed initializations, i.e., M0 and MY, show competitive results when compared with the PCA and ZERO models. Regarding the whitened parameterization, shown in Fig. 18a, looking first at the M0 model, we observe how this model is either better or equal to the ZERO-W model. As we mentioned in Sec. 5.7, we shall not expect the M0 model to surpass the ZERO-W model in predictive performance. Examples of clear improvements are given in Power, Kin8nm, and Yacht datasets. Among these, the most notable case is the Yacht dataset, where the ZERO-W 66
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
boston
concrete
energy
kin8nm
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 2.5
2.4
power
2.3
3.15
3.10
3.05
protein
0.8
0.7
0.6
redwine
1.36
1.38
1.40
yacht
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 2.76
2.74
2.72
2.01
2.00
ZERO-POINTS-M0 W
1.99
1.98
0.99
0.98
0.97
ZERO-POINTS-MY W
2.0
PCA W
(a) Whitened models. concrete energy
boston
1.5
1.0
0.5
ZERO W
kin8nm
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 3.25 3.00 2.75 2.50 2.25
3.15
power
3.10
3.05
protein
0.7
0.6
0.5
redwine
0.4
1.34
1.36
yacht
1.38
1.40
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 3.4
3.2
3.0
2.8
2.00
ZERO-POINTS-M0 NWR
1.99
1.98
1.97
0.99
ZERO-POINTS-MY NWR
0.98
0.97
PCA NWR
4
3
2
1
ZERO NWR
(b) Non-whitened models.
Figure 18: Test log-likelihood (right is better) in all UCI datasets. 67
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
model suffers from posterior collapse (as we will demonstrate later in Sec. 6.2.2), while the proposed M0 initialization mitigates this issue and achieves performance comparable to that of the PCA-W model. This case illustrates the beneficial properties of the proposed initialization. Only in the Concrete dataset, with 4-layers, and the Energy dataset, with 3-layers, M0 models work slightly worse. In the rest of the datasets, we observe slightly better or equal results. Regarding the MY model, we observe that this model can surpass the PCA-W model (Kin8nm). This model also performs similarly to the ZERO-W and M0 models (for instance, in Yacht or Boston). Nevertheless, we observe that the MY model achieves worse results than the M0 or even the ZERO-W models in cases like the Power and Energy datasets. We have further investigated and found out that this occurs due to the initial length-scale selection. In the Energy dataset, for instance, we have observed that the length-scale selected using the proposed algorithm (see Alg. 1) is around 6, which results in a big KLD. Thus, for the same number of epochs, the model might not have achieved a maximum in the ELBO; or the maximum achieved is poor. Looking at the RMSE result in Fig. 33a, we observed that the MY model in this dataset is not degraded. According to the Generalized Variational Inference theory (Knoblauch et al., 2022), the KLD is the uncertainty quantifier. Therefore, it is reasonable to expect that if optimization fails to achieve an appropriate KLD value, the model will struggle to quantify uncertainty accurately. This primarily affects the test LL, which is sensitive to uncertainty estimation, while leaving point-estimate metrics, like RMSE, largely unchanged. Finally, we observed that the Boston dataset exhibits degraded performance when using the PCA prior mean function. This finding suggests that architectural choices should be motivated by modeling considerations rather than optimization-related issues. Although the PCA-W variant consistently performs better, the results on the Boston, Kin8nm, Protein, and Power datasets, where either the M0 or MY models match or outperform the PCA model, indicate that the improvements associated with the PCA prior mean model do not necessarily arise from a change in the underlying statistical model. Consequently, our results support the use of the ZERO prior mean, provided that the model is initialized appropriately. Regarding the non-whitened parameterization, shown in Fig. 18b, we observe similar results. The most notable difference is that, here, the ZERO-NWR yields consistently worse results in a greater number of datasets than the whitened case. While the whitened model showed degraded performance in Kin8nm and Yacht, the NWR model also degrades performance in Boston and Power, and the 4th layer Yacht DGP. The most important conclusion drawn from this figure is that the proposed M0 and MY models effectively correct the poor solutions of the ZERO-NWR model by simply changing the initialization of the underlying shared statistical model. This highlights the benefits of the proposed initialization also for non-whitened models. We further highlight that in Concrete 2 layers, the MY model clearly outperforms the rest of the models. This contrasts with the whitened parameterization where the PCA clearly outperformed the rest of models in this dataset. As mentioned previously, these parameterizations are compared in the appendix for each of the models. Interestingly, we observe that models with the ZERO mean and NWR parameterization show degraded performance compared to the whitened parameterization (e.g., Kin8nm with 4 layers or Power with 5 layers). However, when using our M0 initialization, NWR slightly surpasses the whitened version. 68
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
To provide additional insights on these figures, we summarize the obtained results using whitened models in two different ways: Average model rank, and standard deviation analysis. In the following tables, the best result is highlighted in bold, and the second best is underlined. Table 5 shows the average rank of each model per dataset, in terms of the LL. The ranks were computed for each split. The table shows that the PCA-W model achieves the best results overall. Notwithstanding, the tables also show that the proposed initializations ZERO-PointsM0-W and ZERO-Points-MY-W perform better than the ZERO-W model. Furthermore, we observe that both ZERO-Points-M0-W and ZERO-Points-MY-W achieve results closer to those of the PCA-W model. This validates our design approach: we specifically initialized the inner layers of the ZERO DGP to approximate the PCA mean function, which is expected to yield more similar performance. This confirms further that this initialization successfully narrows the gap between PCA-W DGPs and ZERO-W DGPs, highlighting the beneficial properties of the proposed initialization. To assess the beneficial properties of the proposed initialization with respect to the stability of the resulting model, in Table 6 we display the ranking of the models based on the standard deviation of the results across the splits (lower is better). The table shows that the proposed initializations achieve more stable results than both PCA-W and ZERO-W models. This indicates that the initialization induces a smoother optimization process than those of the ZERO-W and PCA-W models, which is a highly beneficial property. For non-whitened models, we draw the same conclusions as those obtained for whitened ones. In Table 7 show the split-wise average rank in terms of the LL for each model. Table 8 shows the ranks in terms of standard deviations of the test LL. Table 5: Average rank of each whitened model in terms of the test LL for each dataset.
PCA ZERO ZERO-Points-M0 ZERO-Points-MY
boston
concrete
energy
kin8nm
power
protein
redwine
yacht
Overall
3.66 2.33 2.46 1.55
1.64 2.89 3.30 2.17
1.46 2.62 2.90 3.01
2.23 3.34 2.99 1.45
1.30 3.30 2.04 3.36
2.21 3.01 2.56 2.21
2.70 2.64 2.69 1.98
1.40 3.14 3.06 2.40
2.07 2.91 2.75 2.27
Table 6: Average rank of each whitened model in terms of the standard deviations (lower is better) of the LL (horizontal bars) in Fig. 18a.
PCA ZERO ZERO-Points-M0 ZERO-Points-MY
boston
concrete
energy
kin8nm
power
protein
redwine
yacht
Overall
2.00 4.00 3.00 1.00
4.00 2.00 1.00 3.00
2.00 3.00 1.00 4.00
3.00 4.00 1.00 2.00
1.00 4.00 3.00 2.00
4.00 3.00 2.00 1.00
4.00 2.00 3.00 1.00
3.00 2.00 1.00 3.00
2.88 3.00 1.88 2.12
6.2.2 Analyzing the Posterior Collapse Problem on the UCI Datasets We now analyze the previous results to determine whether the models that exhibited poor performance were affected by posterior collapse. Posterior collapse occurs when the model effectively reverts to the prior when making predictions. This phenomenon is typically 69
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Table 7: Average rank of each NWR model in terms of the test LL for each dataset.
PCA ZERO ZERO DatapointsM0 ZERO DatapointsMY
boston
concrete
energy
kin8nm
power
protein
redwine
yacht
Overall
3.51 2.41 2.49 1.59
1.81 2.77 3.04 2.38
1.19 2.46 2.74 3.61
2.30 3.19 2.75 1.76
1.35 3.66 2.31 2.67
2.17 2.83 2.51 2.49
2.20 2.81 2.56 2.42
1.44 3.40 2.90 2.26
2.00 2.94 2.66 2.40
Table 8: Average rank of each NWR model in terms of the standard deviations (lower is better) of the LL (horizontal bars) in Fig. 18a.
PCA ZERO ZERO DatapointsM0 ZERO DatapointsMY
boston
concrete
energy
kin8nm
power
protein
redwine
yacht
Overall
3.00 4.00 2.00 1.00
4.00 2.00 1.00 3.00
1.00 2.00 3.00 4.00
3.00 4.00 2.00 1.00
3.00 2.00 1.00 4.00
4.00 2.00 3.00 1.00
4.00 3.00 1.00 2.00
3.00 4.00 1.00 2.00
3.12 2.88 1.75 2.25
characterized by vanishing KL divergence values together with an increase in the estimated noise variance. To assess whether this behavior is present, we examine the evolution of these quantities during training. We first consider the whitened parameterization. In this case, Fig. 18a shows poor performance of the 5-layer ZERO-W model on the Kin8nm and Yacht datasets, especially on the latter. To further analyze these results, we plot the KLD and estimated noise variance of the ZERO-W and ZERO-Points-M0-W models over the last 50 epochs of the optimization process in Fig. 19, for each of the 20 training–test splits. The figure shows a horizontal line equal to 0 for the Yacht dataset. More precisely, in this dataset, the results show that in 7 out of the 20 splits (35% of them) there is posterior collapse. This causes the model to estimate the prior distribution as its posterior, effectively providing a poor solution. This also explains the high standard deviation of the ZERO-WPprior mean DGP model in Yacht. The value of 2 5 the likelihood variance coincides with 1/N N n=1 (yn ) (horizontal line around 1.0) , which is the optimal coordinate update of the noise parameter. Since we are under a whitening parameterization, the model has learned σ0 = 0, which is coherent with our discussion in Sec. 4.3.2. Importantly, our initialization solves the posterior collapse problem, as Fig. 18a shows. Fig. 19 shows how the M0 model (green lines) presents, in general, higher KLD and smaller likelihood variance than the ZERO model (blue line). This behavior is consistent with a tendency toward posterior collapse, even when it does not occur completely. Furthermore, we compare all the models in Fig. 20, which shows for a representative train/test split, the ELL and the KLD of each method for each iteration epoch, on the Yacht dataset. We observe that the ZERO-W model results in a posterior collapse as the KLD is zero. The ZERO-Points-M0-W and ZERO-Points-MY-W models, which also consider a ZERO prior mean function, but use the proposed initialization, avoid this pathology, achieving a solution that is closer to that of the DGP with the PCA prior mean function, i.e., PCA-W, as indicated by Fig. 18a. 5. The data is standardized by subtracting the mean and dividing by the variance, hence the expected second moment coincides with the variance.
70
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Likelihood Variance
0.062 0.060 0.058 0.056 0.054 0.052 19950
KL Divergence
1400 1300 1200 19950
19960
19970
19980
Epoch
ZERO W
19990
20000
Likelihood Variance
KL Divergence
100 19970
ZERO W
19980
Epoch
19990
20000
19980
Epoch
19990
20000
ZERO-POINTS-M0 W
(b) Kin8nm – Noise variance
1.0 0.8 0.6 0.4 0.2 0.0 19950
200
19960
19970
ZERO W
ZERO-POINTS-M0 W
(a) Kin8nm – KLD
0 19950
19960
19960
19970
ZERO W
ZERO-POINTS-M0 W
(c) Yacht – KLD
19980
Epoch
19990
20000
ZERO-POINTS-M0 W
(d) Yacht – Noise variance
Figure 19: KLD (left) and likelihood variance (right) obtained by the ZERO-W and ZEROPoints-M0-W models in the Kin8nm and Yacht datasets, for each train/test split. No mode collapse is shown in Kin8nm. Seven out of 20 splits suffer from posterior collapse in Yacht. 17500
KL Divergence
12500
100
10000
0 19940 19960 19980
Expected Log Likelihood
200
15000
7500 5000 2500 0 0
0 500 1000 1500 2000 2500
2500 5000 7500 10000 12500 15000 17500 20000
Epoch
ZERO DatapointsM0 W ZERO DatapointsMY W
500
PCA W ZERO W
0
2500 5000 7500 10000 12500 15000 17500 20000
Epoch
ZERO DatapointsM0 W ZERO DatapointsMY W
(a) KLD. The figure zooms in on the last 50 epochs.
PCA W ZERO W
(b) ELL.
Figure 20: KLD and ELL obtained by the four models during training on the Yacht dataset, when using 5 layers, for a representative train/test split. The posterior collapse of the ZERO-W model is shown by the KLD, which becomes zero, and induces a poor ELL term in the ELBO.
71
0.0625
KL Divergence
1400 1350 1300 1250 1200 1150 19950
Likelihood Variance
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
0.0600 0.0575 0.0550
19960
19970
ZERO NWR
19980
Epoch
19990
20000
0.0525 19950
ZERO-POINTS-M0 NWR
Likelihood Variance
KL Divergence
20000
ZERO-POINTS-M0 NWR
0.6 0.4
50
0.2
19960
19970
ZERO NWR
19980
Epoch
19990
20000
19950
Likelihood Variance
KL Divergence
19970
19980
Epoch
19990
20000
ZERO-POINTS-M0 NWR
(d) Boston – Noise variance
1.0
800
19960
ZERO NWR
ZERO-POINTS-M0 NWR
(c) Boston – KLD
0.8
600
0.6
400
0.4
200
0.2
19960
19970
ZERO NWR
19980
Epoch
19990
20000
0.0 19950
19970
ZERO NWR
19980
Epoch
19990
20000
1.0 0.8 0.6 0.4 0.2 0.0 19950
Likelihood Variance
KL Divergence
19960
19960
19970
ZERO NWR
ZERO-POINTS-M0 NWR
(e) Power – KLD
250 200 150 100 50 0 19950
19990
0.8
100
0 19950
19980
Epoch
(b) Kin8nm – Noise variance
1.0
150
0 19950
19970
ZERO NWR
(a) Kin8nm – KLD
200
19960
(g) Yacht – 4 Layers – KLD
19990
20000
ZERO-POINTS-M0 NWR
(f) Power – Noise variance
19960
19970
ZERO NWR
ZERO-POINTS-M0 NWR
19980
Epoch
19980
Epoch
19990
20000
ZERO-POINTS-M0 NWR
(h) Yacht – 4 layers– Noise variance
Figure 21: KLD (left) and likelihood variance (right) obtained by the ZERO-NWR and ZEROPoints-M0-NWR models in the Boston, Power, and Yacht (with 4 layer model) datasets, for each train/test split.
72
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
In the Kin8nm dataset, however, the poorer performance of the ZERO-W DGP model is not accompanied by posterior collapse. As shown in Fig. 19a, the KLD remains significantly above zero throughout training, providing no evidence of posterior collapse. Nevertheless, the ZERO-W model consistently exhibits a lower KLD and a higher estimated noise variance than the ZERO-Points-M0 model (see Fig. 19a and 19b). This suggests a tendency towards collapse, which results in explaining part of the data variability as noise, and accounts for the performance gap observed in Fig. 18a. Notice that the ZERO-W model shows a test loglikelihood of ≈ 1.352, compared to ≈ 1.38 for the ZERO-Points-M0-W model. To test whether this difference is significant, we perform a Wilcoxon signed-rank test which determines if the differences are centered around zero. The test outputs a p−value of 4.7 × 10−5 , which is very close to zero, indicating that the differences are significant. Motivated by this observation, we further study the previous low KLD–high observation noise, which we observe in more scenarios, observation in Appendix C.5.3. We confirm that many of the performance gaps between the ZERO and M0 models come from this behavior, reaffirming the tendency of the ZERO model towards collapsed solutions. Lastly, we show that the poor solutions of the ZERO-NWR models presented in Fig. 18b often correspond to posterior collapse. First, Fig. 21a and 21b indicate that the poor results on the Kin8nm dataset exhibit behavior similar to that observed in the whitened case: although ZERO-NWR does not suffer from posterior collapse, it still shows a significant difference in test log likelihood (p = 3.81 × 10−6 ). In contrast, ZERO-NWR suffers from posterior collapse in 13 out of 20 splits (65%) on the Boston dataset and in 6 out of 20 splits (30%) on the Power dataset. The corresponding KLD and likelihood variance plots are shown in Fig. 21c and 21d for Boston, and in Fig. 21e and 21f for Power. Finally, the model exhibits complete posterior collapse in all 20 splits of the Yacht dataset when using four layers, making it the only configuration with fewer than five layers that displays this pathology. This behavior is illustrated in Fig. 21g and 21h. Most importantly, both proposed initializations, M0 and MY, consistently avoid the posterior collapse suffered by the ZERO model, further highlighting the benefits of our initialization strategy.
7 Conclusions and Future Work In this work, we presented a detailed analysis of the DSVI algorithm used to train DGPs and characterized the frequently encountered problem of posterior collapse. Specifically, the optimization process reaches a local optimum that minimizes the KLD by simply setting the variational parameters equal to those of the prior. Consequently, the model’s predictions collapse into noise, yielding a poor and uninformative solution. We observed that posterior collapse is more likely to happen when using a ZERO prior mean function in the inner layers of the DGP. By contrast, using a PCA prior mean function may sometimes solve this issue. However, this imposes an important constraint on the prior mean functions that may not be desirable in some problems or extensions of the DGP model. An example is the use of transformed GPs instead of standard GPs (Maroñas et al., 2021; Sáez-Maldonado et al., 2023). Moreover, the choice of the model’s prior should be done in such a way that reflects particular characteristics of the learning problem, not simply as a consequence of the optimization problems that arise when fitting the model (Knoblauch et al., 2022). In addition, our analysis indicates that avoiding non-injectivity pathologies 73
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
is not the fundamental reason for the effectiveness of the PCA prior mean function, in the context of DGPs, as previously thought in the literature. In fact, we have shown that a very deep DGP with a whitened parameterization and this prior mean function can lead to posterior collapse if enough noise is introduced in the inner layers. Besides this, we have also shown that the whitened parameterization yields more stable optimization for the ZERO prior mean function model, which is often translated into better predictive distributions. Yet, our UCI results show that the non-whitened parameterization can result in better predictive performance. Since the PCA mean function under this parameterization imposes a different statistical model, this represents a line of research. Based on our observations, we proposed a novel initialization of the DGP model that, by modifying the initial value of the variational mean parameter, mimics the behavior of the PCA prior mean in the inner layers of the model, at initialization. This enables maintaining the typical ZERO prior mean function used in the context of GPs. Furthermore, using a minimal modification of the proposed initialization in the output layer produces an initial predictive distribution that is accurate at the inducing points (selected from a subset of the training set). We then experimentally showed that, thanks to the proposed initialization, the ZERO prior mean DGP can successfully avoid the posterior collapse phenomenon. In addition, we demonstrate that the results obtained, when using our proposed initialization, outperform the results of standard initialized ZERO prior mean DGPs, and that can more closely match those of the PCA prior mean DGP. Even though our proposed initialization alleviates the posterior collapse problem and outperforms the standard initialization of ZERO prior mean DGPs, it presents two main limitations. Firstly, it does not match exactly the predictive performance of the PCA mean function model, which is still better in most of the cases. Secondly, the proposed solution requires solving a system of linear equations. The numerical stability of the solution depends on the number of inducing points and the initial parameters of the kernel, which could compromise the quality of the proposed initialization. Additionally, in the case of initializing the output layer to accurately predict the targets at the inducing points, we also observed a big initial and final KLD term, which can also complicate the fitting of the model. Specifically, a big initial KLD results in a smaller initial ELBO value. Improving the support where our initialization matches the PCA mean function, and improving the numerical stability of the initialization are lines of future research. The first line can be suitably approached through Deep Transformed GPs (Sáez-Maldonado et al., 2023). In any case, by avoiding posterior collapse without constraining the prior mean, the method provides a practical strategy for challenging problems modeled by DGPs and opens a new direction for research on initialization strategies in such models. Moreover, it also enables considering more flexible models in the context of DGP, such as transformed GPs (Maroñas et al., 2021; Sáez-Maldonado et al., 2023), in which using the PCA prior mean function is challenging.
Acknowledgments and Disclosure of Funding The authors acknowledge financial support by Grant PID2022-140189OB-C22 funded by MICIU AEI/10.13039/501100011033 and by NextGenerationEU/PRTR. The authors also ac74
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
knowledge support from the project PID2022-139856NB-I00, funded by MCIN/AEI/10.13039/ 501100011033/FEDER, UE; from project IDEA-CM (TEC-2024/COM-89), funded by the Autonomous Community of Madrid; and from the ELLIS Unit Madrid. The authors also acknowledge computational support from the Centro de Computación Científica-Universidad Autónoma de Madrid (CCC-UAM).
References R. Barceló, C. Alcázar, and F. Tobar. Avoiding mode collapse in diffusion models fine-tuned with reinforcement learning, 2024. URL https://arxiv.org/abs/2410.08315. M. Bauer, M. van der Wilk, and C. E. Rasmussen. Understanding Probabilistic Sparse Gaussian Process Approximations. In Advances in Neural Information Processing Systems, pages 1533 – 1541, 2016. M. Betancourt and M. Girolami. Hamiltonian Monte Carlo for Hierarchical Models. Current Trends in Bayesian Methodology with Applications, pages 79–101, 2015. C. M. Bishop. Pattern Recognition and Machine Learning (Information Science and Statistics). Springer, 2007. ISBN 0387310738. C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra. Weight Uncertainty in Neural Network. In International Conference on Machine Learning, pages 1613–1622, 2015. S. R. Bowman, L. Vilnis, O. Vinyals, A. Dai, R. Jozefowicz, and S. Bengio. Generating Sentences from a Continuous Space. In SIGNLL Conference on Computational Natural Language Learning, pages 10–21, Berlin, Germany, Aug. 2016. T. Bui, D. Hernandez-Lobato, J. Hernandez-Lobato, Y. Li, and R. Turner. Deep Gaussian Processes for Regression using Approximate Expectation Propagation. In International Conference on Machine Learning, pages 1472–1481, 2016. X. Chen, D. P. Kingma, T. Salimans, Y. Duan, P. Dhariwal, J. Schulman, I. Sutskever, and P. Abbeel. Variational Lossy Autoencoder. In International Conference on Learning Representations, 2017. C.-A. Cheng and B. Boots. Variational Inference for Gaussian Process Models with Linear Complexity. In Advances in Neural Information Processing Systems, pages 5184–5194, 2017. C. Cremer, X. Li, and D. Duvenaud. Inference Suboptimality in Variational Autoencoders. In International Conference on Machine Learning, pages 1078–1086, 2018. A. Damianou and N. D. Lawrence. Deep Gaussian Processes. In International Conference on Artificial Intelligence and Statistics, pages 207–215, 2013. G. De Ath, J. E. Fieldsend, and R. M. Everson. What do you mean? the role of the mean function in bayesian optimisation. In Genetic and Evolutionary Computation Conference Companion, pages 1623–1631, 2020. 75
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
S. Depeweg, J.-M. Hernandez-Lobato, F. Doshi-Velez, and S. Udluft. Decomposition of Uncertainty in Bayesian Deep Learning for Efficient and Risk-sensitive Learning. In International Conference on Machine Learning, pages 1184–1193, 2018. D. Dua and C. Graff. UCI Machine Learning Repository, 2019. V. Dutordoir, N. Durrande, and J. Hensman. Sparse Gaussian Processes with Spherical Harmonic Features. In International Conference on Machine Learning, pages 2793–2802, 2020. V. Dutordoir, H. Salimbeni, E. Hambro, J. McLeod, F. Leibfried, A. Artemev, M. van der Wilk, J. Hensman, M. P. Deisenroth, and S. John. GPflux: A Library for Deep Gaussian Processes, 2021. URL https://arxiv.org/abs/2104.05674. D. Duvenaud, O. Rippel, R. Adams, and Z. Ghahramani. Avoiding pathologies in very deep networks. In International Conference on Artificial Intelligence and Statistics, pages 202–210, 2014. J. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson. GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. In Advances in Neural Information Processing Systems, pages 7587 – 7597, 2018. R. Gómez-Bombarelli, J. N. Wei, D. Duvenaud, J. M. Hernández-Lobato, B. SánchezLengeling, D. Sheberla, J. Aguilera-Iparraguirre, T. D. Hirzel, R. P. Adams, and A. AspuruGuzik. Automatic chemical design using a data-driven continuous representation of molecules. ACS central science, 4:268–276, 2018. I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio. Generative Adversarial Nets. In Advances in Neural Information Processing Systems, pages 2672 – 2680, 2014. A. Graves. Practical Variational Inference for Neural Networks. In Advances in Neural Information Processing Systems, pages 2348 – 2356, 2011. O. Hamelijnck, T. Damoulas, K. Wang, and M. Girolami. Multi-resolution Multi-task Gaussian Processes. In Advances in Neural Information Processing Systems, 2019. M. Havasi, J. M. Hernández-Lobato, and J. J. Murillo-Fuentes. Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonian Monte Carlo. In Advances in Neural Information Processing Systems, pages 7517 – 7527, 2018a. M. Havasi, J. M. Hernández-Lobato, and J. J. Murillo-Fuentes. Deep Gaussian Processes with Decoupled Inducing Inputs, 2018b. URL https://arxiv.org/abs/1801.02939. J. He, D. Spokoyny, G. Neubig, and T. Berg-Kirkpatrick. Lagging inference networks and posterior collapse in variational autoencoders. In International Conference on Learning Representations., 2019. J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian Processes for Big data. In Uncertainty in Artificial Intelligence, page 282–290, Arlington, Virginia, USA, 2013. 76
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
J. Hensman, A. G. Matthews, M. Filippone, and Z. Ghahramani. MCMC for Variationally Sparse Gaussian Processes. In Advances in Neural Information Processing Systems, pages 1648 – 1656, 2015. M. Jankowiak and J. Gardner. Neural Likelihoods for Multi-Output Gaussian Processes, 2019. URL https://arxiv.org/abs/1905.13697. A. Javaloy, M. Meghdadi, and I. Valera. Mitigating Modality Collapse in Multimodal VAEs via Impartial Optimization. In International Conference on Machine Learning, pages 9938–9964, 2022. D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. D. P. Kingma and M. Welling. Auto-Encoding Variational Bayes. In International Conference on Learning Representations, 2014. J. Knoblauch, J. Jewson, and T. Damoulas. An Optimization-centric View on Bayes’ Rule: Reviewing and Generalizing Variational Inference. Journal of Machine Learning Research, 23:1–109, 2022. M. Lázaro-Gredilla and A. Figueiras-Vidal. Inter-domain Gaussian Processes for Sparse Inference using Inducing Features. In Advances in Neural Information Processing Systems, pages 1087 – 1095, 2009. M. Lázaro-Gredilla and M. K. Titsias. Variational heteroscedastic Gaussian Process regression. In International Conference on Machine Learning, pages 841–848, 2011. F. Leibfried, V. Dutordoir, S. John, and N. Durrande. A Tutorial on Sparse Gaussian Processes and Variational Inference, 2022. URL https://arxiv.org/abs/2012.13962. Z. Lin, F. Yin, and J. Maroñas. Towards Flexibility and Interpretability of Gaussian Process State-Space Model, 2023. URL https://arxiv.org/abs/2301.08843. J. Z. Liu, S. Padhy, J. Ren, Z. Lin, Y. Wen, G. Jerfel, Z. Nado, J. Snoek, D. Tran, and B. Lakshminarayanan. A Simple Approach to Improve Single-Model Deep Uncertainty via Distance-Awareness. Journal of Machine Learning Research, 24:1–63, 2023. C. Louizos and M. Welling. Structured and Efficient Variational Deep Learning with Matrix Gaussian Posteriors. In International Conference on Machine Learning, pages 1708–1716, 2016. C. Louizos and M. Welling. Multiplicative Normalizing Flows for Variational Bayesian Neural Networks. In International Conference on Machine Learning, pages 2218–2227, 2017. J. Lucas, G. Tucker, R. Grosse, and M. Norouzi. Understanding Posterior Collapse in Generative Latent Variable Models. In Workshop on Deep Generative Models for Highly Structured Data. International Conference on Learning Representations, 2019. 77
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
J. Maroñas and D. Hernández-Lobato. Efficient Transformed Gaussian Processes for NonStationary Dependent Multi-class Classification. In International Conference on Machine Learning, pages 24045–24081, 2023. J. Maroñas, O. Hamelijnck, J. Knoblauch, and T. Damoulas. Transforming Gaussian Processes With Normalizing Flows . In International Conference on Artificial Intelligence and Statistics, pages 1081–1089, 2021. J. Maroñas. The Jacobian and Hessian of the Kullback-Leibler Divergence between Multivariate Gaussian Distributions (Technical Report), 2025. URL https://arxiv.org/abs/ 2506.23996. A. G. d. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. LeónVillagrá, Z. Ghahramani, and J. Hensman. GPflow: A Gaussian Process library using TensorFlow. Journal of Machine Learning Research, 18:1–6, apr 2017. T. Minka. Old and New Matrix Algebra Useful for Statistics. December 2000. I. Murray and R. P. Adams. Slice sampling covariance hyperparameters of latent Gaussian models. In Advances in Neural Information Processing Systems, page 1732–1740, 2010. R. M. Neal. Slice sampling. The annals of statistics, 31:705–767, 2003. L. A. Ortega, S. Rodriguez Santana, and D. Hernández-Lobato. Variational Linearized Laplace Approximation for Bayesian Deep Learning. In International Conference on Machine Learning, pages 38815–38836, 2024. O. Papaspiliopoulos, G. O. Roberts, and M. Sköld. A general framework for the parametrization of hierarchical models. Statistical Science, pages 59–73, 2007. T. Pinder and D. Dodd. GPJax: A Gaussian Process Framework in JAX. Journal of Open Source Software, 7:4455, 2022. S. Popescu, D. Sharp, J. Cole, and B. Glocker. Hierarchical Gaussian Processes with Wasserstein-2 Kernels, 2022. URL https://arxiv.org/abs/2010.14877. C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2005. D. Rezende and S. Mohamed. Variational Inference with Normalizing Flows. In International Conference on Machine Learning, pages 1530–1538, 2015. T. G. J. Rudner, D. Sejdinovic, and Y. Gal. Inter-domain Deep Gaussian Processes. In International Conference on Machine Learning, pages 8286–8294, 2020. F. J. Sáez-Maldonado, J. Maroñas, and D. Hernández-Lobato. Mode Collapse in Variational Deep Gaussian Processes. In NeurIPS 2024 Workshop on Bayesian Decision-making and Uncertainty, 2024. 78
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
H. Salimbeni and M. Deisenroth. Doubly Stochastic Variational Inference for Deep Gaussian Processes. In Advances in Neural Information Processing Systems, pages 4591 – 4602, 2017. H. Salimbeni, V. Dutordoir, J. Hensman, and M. Deisenroth. Deep Gaussian Processes with Importance-Weighted Variational Inference. In International Conference on Machine Learning, pages 5589–5598, 2019. J. Shi, M. Titsias, and A. Mnih. Sparse Orthogonal Variational Inference for Gaussian Processes. In International Conference on Artificial Intelligence and Statistics, pages 1932–1942, 2020. R. Soletskyi, M. Gabrié, and B. Loureiro. A theoretical perspective on mode collapse in variational inference. Machine Learning: Science and Technology, 6:025056, 2025. F. J. Sáez-Maldonado, J. Maroñas, and D. Hernández-Lobato. Deep Transformed Gaussian Processes, 2023. URL https://arxiv.org/abs/2310.18230. M. Titsias. Variational Learning of Inducing Variables in Sparse Gaussian Processes. In International Conference on Artificial Intelligence and Statistics, pages 567–574, 2009. A. Tripp, S. Bacallado, S. Singh, and J. M. Hernández-Lobato. Tanimoto Random Features for Scalable Molecular Machine Learning. In Advances in Neural Information Processing Systems, pages 33656–33686, 2023. I. Ustyuzhaninov, I. Kazlauskaite, M. Kaiser, E. Bodin, N. Campbell, and C. Henrik Ek. Compositional uncertainty in Deep Gaussian Processes. In Uncertainty in Artificial Intelligence, pages 480–489, 2020. A. van den Oord, N. Kalchbrenner, L. Espeholt, k. kavukcuoglu, O. Vinyals, and A. Graves. Conditional Image Generation with PixelCNN Decoders. In Advances in Neural Information Processing Systems, pages 4797 – 4805, 2016a. A. van den Oord, N. Kalchbrenner, and K. Kavukcuoglu. Pixel Recurrent Neural Networks. In International Conference on Machine Learning, pages 1747–1756, 2016b. M. van der Wilk, V. Dutordoir, S. John, A. Artemev, V. Adam, and J. Hensman. A Framework for Interdomain and Multioutput Gaussian Processes, 2020. URL https: //arxiv.org/abs/2003.01115. Y. Wang, D. Blei, and J. P. Cunningham. Posterior Collapse and Latent Variable Nonidentifiability. In Advances in Neural Information Processing Systems, pages 5443–5455, 2021. A. Wu, S. Nowozin, E. Meeds, R. E. Turner, J. M. Hernandez-Lobato, and A. L. Gaunt. Deterministic Variational Inference for Robust Bayesian Neural Networks. In International Conference on Learning Representations, 2019. J. Xu, D. Zeng, and J. Paisley. Sparse Inducing Points in Deep Gaussian Processes: Enhancing Modeling with Denoising Diffusion Variational Inference. In International Conference on Machine Learning, pages 55490–55500, 2024. 79
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
J. Xu, Q. Zhao, J. Paisley, and D. Zeng. Diffusion Bridge Variational Inference for Deep Gaussian Processes. In International Conference on Learning Representations, 2026. H. Yu, Y. Chen, B. K. H. Low, P. Jaillet, and Z. Dai. Implicit Posterior Variational Inference for Deep Gaussian Processes. In Advances in Neural Information Processing Systems, 2019. H. Yu, D. Liu, B. K. H. Low, and P. Jaillet. Convolutional normalizing flows for deep Gaussian Processes. In International Joint Conference on Neural Networks, pages 1–6. IEEE, 2021. X. Zhu, J. R. Gardner, and D. Bindel. Efficient variational Gaussian processes initialization via kernel-based least squares fitting. In NeurIPS Workshop on Gaussian Processes, Spatiotemporal Modeling, Decision-making Systems, pages 1–4, 2022. X. Zhu, K. Wu, N. Maus, J. Gardner, and D. Bindel. Variational Gaussian Processes with Decoupled Conditionals. In Advances in Neural Information Processing Systems, pages 46191–46211, 2023. M. A. Álvarez, L. Rosasco, and N. D. Lawrence. Kernels for Vector-Valued Functions: A Review. Foundations and Trends® in Machine Learning, 4:195–266, 2012. ISSN 1935-8237.
80
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Appendix A. Coordinate Updates for the Noise Parameter By noting that a one-dimensional Gaussian likelihood function can be compactly expressed through: N Y N y n | fn , σ 2 = N Y | f, σ 2 I (75) n=1
The objective function of a one-dimensional SVGP can be written as: Z ELBO = q(f) log N Y | f, σ 2 I df − KLD [q(u) || p(u)] i 1 h −1 = log N y | µqf , σ 2 I + Kqf − Tr σ 2 I Kqf − KLD [q(u) || p(u)] 2 Getting the differential over σ 2 : T i 1 h N Y − µqf + Tr Kqf dσ 2 dσ2 ELBO = − + 2 Y − µqf 2 2σ Since this is a scalar-valued, scalar argument function, the differential directly identifies the gradient. Thus: T i N 1 h + 2 Y − µqf Y − µqf + Tr Kqf = 0 2 2σ T i 1 h 2 Y − µqf Y − µqf + Tr Kqf σ = N 1 || Y − µqf ||22 +Tr Kqf σ2 = N −
This is equivalent to the expected value under q(fn ) of the average squared distance between the target y n and the DGP output fn . For the ELL term in Eq. 33 for the DGP, the derivation is exactly the same. Note that in the derivation we have done, the KLD does not take part in the differential. Now for any distribution q(f) with finite mean and covariance it holds: Z
q(f) log N Y | f, σ 2 I df
i 1 h −1 = log N y | µqf , σ 2 I + Kqf − Tr σ 2 I Kqf 2
(76)
Thus, if we substitute q(f L ) by the distribution over the last layer of a DGP, then the ELL takes the above expression and so the coordinate update for σ 2 is exactly the same as that of a SVGP using the predictive mean and predictive covariance of the non-Gaussian distribution q(f L ). For a C Multi-output SVGP one can follow a similar reasoning. Here, however, while we can express the likelihood function in a compact way using Kronecker products, the final expression does not match the same ordering as the latent processes of the variational posterior. Thus, here we explicitly write the likelihood function as a sum over samples. The 81
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
noise parameter of the likelihood function is expressed through Σ, assuming the most general case, which is a dense matrix. This gives: ELBO =
N X n=1
1 log p(yn | µqf n , Σ + K qf n ) − Tr Σ−1 K qf n − KLD [q(u) || p(u)] 2
Here now yn ∈ RC and is a vector µqf n containing the variational mean of each process at point n. In contrast, K qf n is a C × C matrix containing dependencies between all the processes at point n. The first order differential w.r.t. the noise matrix: N h iT T 1X vec Σ−1 yn − µqf n yn − µqf n Σ−1 − Σ−1 + Σ−1 K qf n Σ−1 dvec Σ dΣ ELBO = 2 n=1
From here, we can identify the gradient and undo the vec operator to retrieve the optimal update in terms of a matrix. Setting the expression to zero and operating yields: N X n=1
Σ−1 yn − µqf n
yn − µqf n
T
Σ−1 − Σ−1 + Σ−1 K qf n Σ−1 = 0
N
Σ=
T 1 X n y − µqf n yn − µqf n + K qf N n=1
If one now assumes that all the outputs share the same noise parameter, or that we want independent outputs each with its own noise, then either the differential must be recomputed or just use the differential of a diagonal matrix, which is easier than using Lagrange multipliers (Minka, 2000).
82
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Appendix B. Additional Details of the Proposed Initialization This appendix provides additional insights into how to solve the systems of equations that result from the proposed initialization. While we have only evaluated the most efficient version, for completeness, we provide the full analysis and derivations. This might be used to initialize the non-whitened models implemented in GPYTORCH or GPJAX or to increase the support of points Xs where the ZERO model mimics the PCA model at initialization. From the main paper, we know that the three systems of equations that result from our idea are: −1 −1 KXZ KZZ m = Xs − KXs Z KZZ Z −1 T KXZ LZZ mv = Xs −1 KXZ KZZ mr = Xs
(NW) (W) (NWR)
First, we noted that we do not require the set of points Xs to be the same as X in order to define the system, but it is required that dim(Xs ) = dim(X) for the three systems to be well-defined. However, since it is counterintuitive to require the ZERO DGP to mimic the PCA DGP at different points in the domain, our first assumption is X = Xs . With this in mind, we will now see that these systems can be solved in many ways. We will analyze all the options and show how we can yield a linear solvable system that can be solved without requiring numerical gradient-based optimization. As we mentioned, this is not a restriction of our approach but a decision we have made for our experiments. Thus, there is freedom to go through other approaches that might be interesting for other reasons (for example, incrementing the support Xs where the desired initialization is satisfied). B.1 Equations and Parameters from the System of Equations The resulting systems of equations have a total of dim(Xs ) equations, which correspond to the points where we aim to mimic the PCA DGP. Ideally, we would like Xs = X . However, since this would imply having an infinite number of equations, we start by fixing a finite number of points Xs , yet to be determined, where the equation will be satisfied. On the other hand we have a total of dim(Xs ) + 2 dim(Z) + dim(ν) variables: the Xs where the posterior is evaluated, the value of the inducing points Z, the variational mean mv , mr , m and the kernel hyper-parameters ν. Note that dim(m) = dim(Z). Since there are more variables than equations, solving the system for all the variables results in an infinite number of solutions. Thus, we need to fix some values for the variables and solve for the others. B.2 Fixing the Kernel Parameters Having fixed the number of equations to dim(Xs ), we need to fix a number of variables equal to the number of equations for the system to be compatible. If dim(Xs ) is greater than the number of parameters, then the system is overestimated, which results in a system without a solution unless some of the equations are linear combinations of the others. The immediate consequence is that fixing values for Xs , Z and m and solving for ν is not a good choice, because the number of kernel hyper-parameters is usually small, and that 83
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
would imply that our initialization only resembles the PCA model at very few points. For instance, for an RBF kernel, we would only be able to mimic the PCA DGP at two points in the domain (assuming the kernel does not implement automatic relevance determination). Moreover, this case entails solving a nonlinear system of equations using gradient-based optimization, as the kernel parameters appear nonlinearly in the equation. Consequently, by fixing ν, we now have dim(Xs ) equations and dim(Xs ) + 2 dim(Z) variables, and so we need to fix more variables. We describe two alternatives. B.3 Option 1: 2 dim(Z) Equations The first option we can think about is a system in which we have 2 dim(Z) equations. In other words, the number of points at which we want to mimic the PCA DGP is twice the number of inducing points. Since the number of variables dim(Xs ) + 2 dim(Z) is greater than the number of equations, we have three possibilities for fixing values in order to solve the system: • Fix values for m and Z and solve the system for Xs . • Fix either m or Z and a total of dim(Z) points from Xs and solve the system for the rest of Xs locations and either m or Z, depending on what we have fixed. • Fix Xs and solve the system for m and Z. However, this approach presents several difficulties. The first one is that when fixing a value for m we need to choose one far from 0, since these are the values that can make the model collapse and because that would imply the solution Xs = 0 for the NWR and W model, with independence of the value selected for Z. Second, solving the system for Xs may not be useful since optimization can solve for a Xs that is far from our training data (since Xs can be any point in the domain), and we would like our model to be initialized at points next to the training data or the points where we are likely to make predictions. Third, fixing Xs or Z can result in a system with no closed-form solution, due to the nonlinear relation in the equation from Xs and Z. In this case, the system would be solved by numerical gradient-based optimization, minimizing least squares. For example fixing Xs and solving for mr and Z can be done through: −1 argmin||KXs Z KZZ mr − Xs ||p
(77)
mr ,Z
which shows that this approach also has a cubic cost per gradient update. Thus, none of the options outlined in this section are good due to either an efficiency bottleneck or the plausible set of solutions achieved. B.4 Option 2: dim(Z) Equations From the previous section, we have learned that it is interesting to solve the system at points in which Xs is representative of the places where we make predictions, and we would like to avoid cubic cost and numerical gradient-based optimization. 84
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
As we now see, if we sacrifice the total number of points where we mimic the PCA DGP to the half, i.e. from 2 dim(Z) to dim(Z), we can yield an easier system to solve. Again, we have three options. • Fix a value for m and Z and solve for Xs , which is not a good option since it requires numerical gradient-based optimization and does not ensure initializing the model in a desired region of the space X . • Fix a representative value for Xs (for example a subset of the training data or some iterations of k-means) and m and solve the system for Z, which is not a good option because we require non-linear optimization of the least square function, and also require m to be different from 0. • Fix values for Xs and Z and solve for m, which results in a linear system of equations that can be solved, at most, in two cubic operations, which makes it the most computationally efficient option for solving the systems. With this last option, we have that the solution for the optimal value m, which would resemble the PCA mean function at initialization, at some selected Xs and Z is given by (for the three parameterizations considered): −1 −1 −1 m = KXs Z KZZ Xs − KXs Z KZZ Z h T i−1 mv = KXs Z L−1 Xs ZZ −1 −1 mr = KXs Z KZZ Xs
(NW) (W) (NWR)
Thanks to fixing Xs and Z, we yield linear systems with exact solutions, but require two cubic operations in the number of inducing points, one per matrix inversion. Nevertheless, choosing Xs and Z to be exactly the same points, we can yield a much more efficient solution. There are two options here: either select Xs to be the inducing points Z (obtained through k-means at initialization), or select a subset from Xs to be the inducing points. In any case since now Xs = Z := Xxz , the solutions to the systems are: m = 0;
−1 using KXxz Xxz KX =I xz Xxz
T mv = L−1 Xxz Xxz Xxz ; using KXxz Xxz = LXxz Xxz LXxz Xxz
mr = Xxz ;
(78)
−1 using KXxz Xxz KX =I xz Xxz
B.5 Selection of the Initial Length-scale In Sec. 5.6, we highlighted the importance of selecting an appropriate initial length-scale for the ZERO-Points-MY model to obtain good initial predictive performance. Nevertheless, this choice can entail a trade-off between the achieved RMSE and the resulting KLD. Fig. 22 presents the variation of the initial KLD and RMSE of the ZERO-Points-MY-W model with respect to the initial length-scale.
85
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
0.16
0.6
0.12 0.10
RMSE
KL Divergence (x1e-7)
0.7 0.14
0.08
0.5
0.4
0.06
0.3
0.04 0.02
0.2
0.00
0.1 0.000
0.025
0.050
0.075
0.100
0.125
0.150
0.175
0.000
0.200
0.025
0.050
0.075
0.100
0.125
0.150
0.175
0.200
Lengthscale
Lengthscale
(a) KLD.
(b) RMSE.
Figure 22: KLD and RMSE obtained by the two layer ZERO-Points-MY-W model when varying the length-scale ℓ at initialization in the steps dataset.
86
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
Appendix C. Further Experimental Results This section completes some of the experiments or simulations of the main paper. More precisely, we first provide additional figures for the theoretical analysis performed in Sec. 4. Then, we provide further observations from the results presented in Sec. 6, both for the toy dataset and UCI real-world datasets. C.1 Additional Figures for Section 4 In this section, we provide figures that complete the analysis of section 4. Fig. 23 complements the analysis of the initial predictive distribution according to a different initialization of Sv . Fig. 24 shows the wiggle effect when Sl = I and SL = 10−5 I in non-whitened models, for different number of inducing points. In Fig. 25, we swap the values of S with respect to the last figure, to show the wiggle effect when Sl = 10−5 I and SL = I in non-whitened models. To create Fig. 26, where predictive variance goes beyond the prior, we have initialized S in a way such that the prior information coming from the kernel is broken. Remember, we saw that predictive variance depends on KZZ − S. For S = I being diagonal, this removes the information contained in the diagonal of KZZ but keeps covariances, effectively reducing predictive variance far from the inducing points. To yield S which breaks prior information, we need to initialize S such that i) is positive semi-definite and ii) KZZ − S also breaks the information contained in the off-diagonal. This can be achieved by a diagonal matrix D with entries taking the values of 1 and −1 and setting: S = DT KZZ D
(79)
If D is an orthogonal matrix, then DT KZZ D results in a positive semi-definite matrix. By this initialization S takes values in the off-diagonal such that some covariances from KZZ are removed, and some are duplicated. This breaks information contained in the prior and inducing point location and effectively exacerbates posterior variance. Thus, one must take care with random initialization of S. C.2 Coordinate Updates for 64 Inducing Points Fig. 27 shows coordinate updates for 64 inducing points. In general, we observe the same behavior as with the 5 inducing points. As a difference, we observe that more inducing points yield lower noise variance as expected, since the uncertainty coming from approximation error is reduced. Particularly interesting is the fact that for S = I, 64 inducing points and the order σ 2 − S − m, the variance explodes in the first iteration to a value around 17. This shows that early in the training, gradient-based optimization can be directed to a zone where posterior collapse is likely.
87
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1.5
Predictive distribution of the 2nd Layer GP Predictive Distribution of the 1st Layer GP Ignoring the Predictive Variance of the First Layer 1.5 1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0 data posterior mean posterior variance Sample from posterior Sample from diagonal posterior Inducing points
−0.5 −1.0 −1.5 0.0
0.2
0.4
0.6
0.8
0.0
0.0
−0.5
−0.5
−1.0
−1.0
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
Full DGP Predictive Distribution
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(a) ZERO mean DGP with Sv = 0.8I 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(b) ZERO mean DGP with Sv = 1.3I 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(c) PCA mean DGP with Sv = 0.8I 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(d) PCA mean DGP with Sv = 1.3I
Figure 23: ZERO (top) and PCA (bottom) DGP models with the output layer variational covariance initialized to Sv = 1.3I and Sv = 0.8I , in the whitened parameterization.
88
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1.5
Predictive distribution of the 2nd Layer GP Predictive Distribution of the 1st Layer GP Ignoring the Predictive Variance of the First Layer 1.5 1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
−0.5
−0.5
−1.0
−1.0
0.0 data posterior mean posterior variance Sample from posterior Sample from diagonal posterior Inducing points
−0.5 −1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
Full DGP Predictive Distribution
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(a) ZERO DGP with 5 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(b) PCA DGP with 5 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(c) ZERO DGP with 20 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(d) PCA DGP with 20 inducing points.
89
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(e) ZERO DGP with 100 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(f) PCA DGP with 100 inducing points.
Figure 24: ZERO and PCA mean DGPs with the variational covariance S = I for the inner layer and S = 10−5 I in the output layer, varying the number of inducing points in the non-whitened parameterization.
90
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1.5
Predictive distribution of the 2nd Layer GP Predictive Distribution of the 1st Layer GP Ignoring the Predictive Variance of the First Layer 1.5 1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
−0.5
−0.5
−1.0
−1.0
0.0 data posterior mean posterior variance Sample from posterior Sample from diagonal posterior Inducing points
−0.5 −1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
Full DGP Predictive Distribution
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(a) ZERO mean DGP with 10 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(b) PCA mean DGP with 10 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(c) ZERO mean DGP with 100 inducing points. 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(d) PCA mean DGP with 100 inducing points.
Figure 25: ZERO and PCA mean DGPs with 10 and 100 inducing points using the GPFLOW non-whitened parameterization, with S = I in the output layer and S = 10−5 I in the inner layers.
91
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1.5
Predictive distribution of the 2nd Layer GP Predictive Distribution of the 1st Layer GP Ignoring the Predictive Variance of the First Layer 1.5 1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
−0.5
−0.5
−1.0
−1.0
0.0 data posterior mean posterior variance Sample from posterior Sample from diagonal posterior Inducing points
−0.5 −1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
Full DGP Predictive Distribution
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(a) ZERO mean DGP 1.5
1.5
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
0.2
0.4
0.6
0.8
1.0 −1.5 0.0
(b) PCA mean DGP
Figure 26: ZERO (top) and PCA (bottom) DGP models with the output layer variational covariance initialized to break prior information in the GPFLOW’s unwhitened parameterization. A horizontal line at 1 shows how variance far from the inducing point goes beyond the prior.
92
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
1
2
3
Iterations
4
5
PCA Mean Zero Mean
ZERO MEAN PCA MEAN
PCA Mean Zero Mean
Noise Variance
1.4
m 1
0
−1
S − KZZ
0.00 −0.25 −0.50 −0.75
0
1
2
3
Iterations
4
5
(a) Order of updates: σ 2 -S-m, with initialization S = 10−5 I, m = 0
1.2 1.0
17.31
0.8 0.6 0.4 0.2 0.0 0
1
2
3
Iterations
4
5
PCA Mean Zero Mean
ZERO MEAN PCA MEAN
PCA Mean Zero Mean
Noise Variance
1.4
m 1
0
−1
S − KZZ
0.0 −0.2 −0.4 −0.6 −0.8
0
1
2
3
Iterations
4
5
(b) Order of updates: σ 2 -S-m, with initialization S = I, m = 0
1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
1
2
3
Iterations
4
5
PCA Mean Zero Mean
ZERO MEAN PCA MEAN
PCA Mean Zero Mean
Noise Variance
1.4
m 1
0
−1
S − KZZ
0.0 −0.2 −0.4 −0.6 −0.8
0
1
2
3
Iterations
4
5
2
(c) Order of updates: S-m-σ . Initialization does not matter since coordinate updates on S does not depend on previous variational parameter values.
Figure 27: Coordinate updates experiment with 64 inducing points.
93
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
C.3 Additional Toy Experiments Illustrating other Parameter Initializations Fig. 28 shows toy results using Sv = I in all layers. Following our previous analysis, we can observe how this configuration also yields poorer solutions in the ZERO model and slightly higher predictive variance in the PCA model. For the experiment with Sl = I and SL = 10−5 I, Fig. 29 shows the layer-wise KLD of the model during training. We observe that in the cases with 20 and 100 inducing points and 3 − 5 layers, the inner layers of the models essentially do not modify their KLD, effectively not learning, which results in the final posterior collapse. Using a lower number of inducing points (in this figure, 5) and 3 layers, we observe how at the beginning of the training the model increases its KLD to model the data, thus escaping from the posterior collapse effect. Nevertheless, the final solution is still suboptimal due to noisy optimization. With 5 layers and 5 inducing points, we observe how even if the model escapes from collapse, at the end all layers are tending to collapse. The PCA-W model never falls into the zero KLD local minimum.
94
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
3 Layers
5 Layers
100 Inducing
20 Inducing
5 Inducing
2 Layers
100 Inducing
20 Inducing
5 Inducing
ZERO-W
PCA-W
Figure 28: Predictive distributions of the ZERO-W and PCA-W DGP models when initialized with Sv = I in all layers. 95
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
5 Inducing
2 Layers 30 20
20 Inducing
5 Layers
KLD Layer 0 KLD Layer 1 KLD Layer 2 Total KLD
30
KLD Layer 0 KLD Layer 1 KLD Layer 2 KLD Layer 3 KLD Layer 4 Total KLD
20
20
10 KLD Layer 0 KLD Layer 1 Total KLD
0 0
100 Inducing
3 Layers
2500
5000
10
10 0
0
7500 10000
0
2500
5000
7500 10000
KLD Layer 0 KLD Layer 1 KLD Layer 2 Total KLD
0
100
100
75
75
50
50
50
25
25
25
KLD Layer 0 KLD Layer 1 Total KLD
0 0
KLD Layer 0 KLD Layer 1 Total KLD
400
200
KLD Layer 0 KLD Layer 1 KLD Layer 2 Total KLD
400
0
2500 5000 7500 10000
KLD Layer 0 KLD Layer 1 KLD Layer 2 KLD Layer 3 KLD Layer 4 Total KLD
400
200
0 2500 5000 7500 10000
KLD Layer 0 KLD Layer 1 KLD Layer 2 KLD Layer 3 KLD Layer 4 Total KLD
75
2500 5000 7500 10000
200
0
7500 10000
0 0
0
5000
100
0
2500 5000 7500 10000
2500
0 0
2500 5000 7500 10000
0
2500 5000 7500 10000
5 Inducing
ZERO-W KLD Layer 0 KLD Layer 1 Total KLD
40
20 Inducing
KLD Layer 0 KLD Layer 1 KLD Layer 2 KLD Layer 3 KLD Layer 4 Total KLD
60
40 40
20 20 0
20
0 0
100 Inducing
80
KLD Layer 0 KLD Layer 1 KLD Layer 2 Total KLD
60
2500
5000
0
7500 10000
0
100
100
75
75
50
50
25
KLD Layer 0 KLD Layer 1 Total KLD
0 0
400
200
5000
7500 10000
KLD Layer 0 KLD Layer 1 KLD Layer 2 Total KLD
2500 5000 7500 10000
5000
7500 10000
50
0
0 0
2500 5000 7500 10000
KLD Layer 0 KLD Layer 1 KLD Layer 2 Total KLD
400
0
2500 5000 7500 10000
KLD Layer 0 KLD Layer 1 KLD Layer 2 KLD Layer 3 KLD Layer 4 Total KLD
400
200
0 0
2500
KLD Layer 0 KLD Layer 1 KLD Layer 2 KLD Layer 3 KLD Layer 4 Total KLD
100
200
0
0
25
2500 5000 7500 10000
KLD Layer 0 KLD Layer 1 Total KLD
2500
0 0
2500 5000 7500 10000
0
2500 5000 7500 10000
PCA-W
Figure 29: Layer-wise KLD during training of both ZERO-W and PCA-W models initialized −5 with Slv = I in the inner layers and SL v = 10 I in the output layer. In some cases, the ZERO-W model can not escape from the local optimum of zero KLD in some layers, especially with a big number of inducing points. With a lower number of inducing points, it is more 96 likely that the ZERO-W escapes from the posterior collapse, as in the case of 2 and 3 layers and 5 inducing points.
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
C.4 Additional Toy Experiments Illustrating Optimization Instabilities Depending on the Parameterization In the experimental section, we run the models using a learning rate of 10−4 . We have also run the same experiments using 10−2 . With this learning rate, as we see in Table 9, the ZERO-W model suffers posterior collapse, while the other models avoid it, see Fig. 30. In the case of NWR models, all of them (including the PCA) suffer from optimization instabilities that result in posterior collapse. These NWR models clearly obtain worse results than their whitened version. While we do not observe exactly a KLD of zero in the PCA-NWR model, we do observe a high likelihood variance. By inspecting the KLD at each layer, see Fig. 31, we observe how the PCA NWR model collapses two of the layers and compensates this behavior by increasing the variance. This is an indicator of how even the PCA model might result in a collapsed distribution when noise is injected in the optimization. Note that, in this work, we have presented two examples of a PCA model suffering collapse: one coming from the NWR parameterization, and another in a whitened parameterization where a small number of inducing points is used. Here, however, the KLD from the NWR models is higher than the whitened counterparts. This comes from the peaks in the optimization algorithm, see Fig. 32. Model ZERO-W ZERO-NWR PCA-W PCA-NWR ZERO-Points-M0-W ZERO-Points-M0-NWR ZERO-Points-MY-W ZERO-Points-MY-NWR
KLD
0.0010 869.0811 144.5761 3310.1553 184.2085 92017.6629 144253.7899 2164.8801
Lik. Var. 0.9967 1.1559 0.0053 1.1707 0.0058 1.0156 0.0321 1.5375
RMSE
0.9896 0.9895 0.0420 0.9895 0.0580 0.9903 0.0946 0.9901
Table 9: Test metrics achieved by all the methods in the toy dataset using a higher learning rate (λ = 10−2 ).
97
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
↑
↑
↑
1.0
Likelihood Variance
KL Divergence
2500
2000
1500
1000
500
0
0.8
0.6
0.4
0.2
0.0
0
2000
4000
6000
8000
10000
0
2000
4000
Epoch ZERO W
PCA W
6000
8000
10000
Epoch
ZERO DatapointsM0 W
ZERO DatapointsMY W
ZERO W
(a) KLD per epoch.
PCA W
ZERO DatapointsM0 W
ZERO DatapointsMY W
(b) Likelihood variance per epoch.
Figure 30: KLD and likelihood variance during training of the 5-layer whitened models in the toy dataset using λ = 10−2 .
600
400
Likelihood Variance
KL Divergence
2.5
Layer 1 Layer 2 Layer 3 Layer 4 Layer 5
800
200
0
2.0 1.5 1.0 0.5 0.0
9991 9992 9993 9994 9995 9996 9997 9998 9999 10000
0
2000
4000
Epoch
6000
8000
10000
Epoch
(a) Layer-wise KLD of the PCA-NWR in the last (b) Likelihood variance per epoch of the PCA10 epochs of training. NWR model.
Figure 31: Layer-wise KLD in the last 10 epochs of training and likelihood variance value during training of the 5-layer PCA-NWR trained using λ = 10−2 . 14 12
Log KL Divergence
Log KL Divergence
13 12 11 10 9 8
11 10 9 8 7
7 0
2
4
6
0
8
2
ZERO W ZERO NWR
PCA W PCA NWR
ZERO DatapointsM0 W ZERO DatapointsM0 NWR
4
6
8
Epoch
Epoch ZERO DatapointsMY W ZERO DatapointsMY NWR
ZERO W ZERO NWR
(a) Learning rate: 10−2
PCA W PCA NWR
ZERO DatapointsM0 W ZERO DatapointsM0 NWR
ZERO DatapointsMY W ZERO DatapointsMY NWR
(b) Learning rate: 10−4
Figure 32: Log KLD during the first 10 epochs in the toy dataset. Non-visible curves overlap with visible ones. This figure shows how the NWR model with λ = 10−2 results in a very unstable optimization, manifested through the peaks in the learning curve. 98
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
C.5 Additional UCI Results In this section, we present additional experimental results on UCI datasets that complement the analysis performed in Sec. 6.2. We include the results of the models in terms of the RMSE, compare the results between whitened and non-whitened parameterizations, and analyze the behavior of the models in some particular, interesting splits. C.5.1 RMSE Results Fig. 33a and Fig. 33b shows RMSE results for UCI experiments on the whitened and nonwhitened parameterizations. Furthermore, Table 10 and Table 11 show the average split-wise ranks of each model in terms of the RMSE, with both parameterizations respectively. Table 10: Average rank of each W model in terms of the test RMSE for each dataset. boston
concrete
energy
kin8nm
power
protein
redwine
yacht
Overall
3.55 2.35 2.21 1.89
1.74 2.94 3.05 2.27
1.57 2.64 2.85 2.94
2.21 3.34 2.99 1.46
1.27 3.34 2.04 3.35
2.34 2.70 2.45 2.51
2.55 2.61 2.65 2.19
1.48 3.09 2.99 2.45
2.09 2.88 2.65 2.38
PCA ZERO ZERO-Points-M0 ZERO-Points-MY
Table 11: Average rank of each NWR model in terms of the test RMSE for each dataset.
PCA ZERO ZERO DatapointsM0 ZERO DatapointsMY
boston
concrete
energy
kin8nm
power
protein
redwine
yacht
Overall
3.36 2.62 2.23 1.79
1.91 2.77 2.89 2.42
1.31 2.70 2.74 3.25
2.25 3.19 2.77 1.79
1.35 3.64 2.31 2.70
2.48 2.42 2.64 2.46
2.21 2.54 2.64 2.61
1.82 3.42 2.65 2.10
2.09 2.91 2.61 2.39
C.5.2 Comparing Whiten and Non-Whiten Models The results obtained in the toy dataset (see Sec. 6.1) showed that the NWR in all the models achieved the worst performance. Moreover, these results are confirmed in the UCI experiments (see Sec. 6.2) where we have observed that the ZERO NWR collapses more than its whitened version. We now compare both parameterizations for the rest of the models, namely PCA, M0, and MY. We show the comparison of both parameterizations for each model in Fig. 34. The only conclusion we might draw is that the ZERO-NWR is worse than the ZERO-W, with clear collapse in some cases (Power, Boston, Yacht, and Kin8nm). Only Protein benefits from the NWR parameterization with this mean function, where it works better in all depths except 5. In the rest of the parameterizations, there is no clear winner. In the PCA model, we see that the Yacht dataset is better fitted using the whitened parameterization, in contrast to the Energy dataset. In the rest of the datasets, both parameterizations perform similarly, with one being the winner depending on the dataset. It must be noticed that, within a dataset, different DGP depths perform differently depending on the parameterization, as in Concrete, where 2, 3 layers work better with the whitened, but 4 works better with non-whitened. The 99
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs boston concrete energy kin8nm
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 2.4
2.6
2.8
3.0
power
5.4
5.6
protein
5.8
0.425 0.450 0.475 0.500 0.059 0.060 0.061 0.062 0.063
redwine
yacht
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2
3.6
3.7
3.8
1.82
1.84
ZERO-POINTS-M0 W
1.86
0.640
0.645
0.650
ZERO-POINTS-MY W
0.655
0
PCA W
(a) Whitened models. concrete energy
boston
2
4
6
ZERO W
kin8nm
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 4
power
6
5.2
5.4
5.6
protein
5.8
0.35
0.40
0.45
redwine
0.50
0.060
0.062
yacht
0.064
L=5
L=5
L=4
L=4
L=3
L=3
L=2
L=2 4
6
8
10
1.80
ZERO-POINTS-M0 NWR
1.82
1.84
1.86
0.640 0.645 0.650 0.655
ZERO-POINTS-MY NWR
PCA NWR
0
5
10
ZERO NWR
(b) Non-whitened models.
Figure 33: Test RMSE (left is better) in all UCI datasets. 100
15
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
reason behind performance differences for this mean function could rely on the different statistical models being implemented by each parameterization. This is a research direction for future work. For the M0 and MY models, the tendency is similar. The MY model performs better using the NWR parameterization in Power or Protein. The whitened works better in Kin8nm. In the M0 both parameterizations work on pair. To confirm that this performance gap between parameterizations does not come from a poor solution coming from a noisy algorithm, we select (by observing Fig. 34) a model configuration where the NWR parameterization outperforms the whitened one and plot its learning curves. In Fig. 35, we represent the PCA and ZERO-Points-M0 models in a 5 layer DGP trained on the Power dataset. The figure shows that the NWR parameterization presents a very noisy behavior, confirming that the improved results do not come from a stable optimization procedure. This confirms that, although the NWR can achieve better results, it still yields to a poor optimization procedure. This suggests that ZERO NWR is more prone to posterior collapse due not only to its initialization, but also to the noisy gradients induced by the NWR parameterization. boston
L=5
concrete
energy
kin8nm
boston
concrete
energy
kin8nm
L=5
L=4
L=4
L=3
L=3
L=2
L=2 2.5
L=5
2.4
3.15
power
3.10
protein
3.05
0.6
redwine
0.4 1.36
1.38
yacht
1.40
3.0
2.5
power
3.150 3.125
0.70
protein
redwine
0.65
1.350 1.375
yacht
L=5
L=4
L=4
L=3
L=3
L=2
L=2 2.72
L=5
2.70 2.00
1.98
PCA W
0.99 0.98 0.97
PCA NWR
0.6
(a) PCA DGP
boston
concrete
energy
0.4
0.2
3.5
3.0
boston
kin8nm
2.00
1.98
0.99 0.98 0.97
ZERO W
ZERO NWR
(b) ZERO DGP
concrete
energy
4
2
kin8nm
L=5
L=4
L=4
L=3
L=3 L=2
L=2 2.35 2.30 2.25
L=5
power
3.1753.1503.125 0.725 0.700 0.675 1.36
protein
redwine
2.35 2.30 2.25
1.38
power
yacht
3.15
3.10
protein
0.8
0.7
redwine
1.36 1.38 1.40
yacht
L=5
L=4
L=4
L=3
L=3 L=2
L=2 2.74
2.72
2.00
1.98
ZERO-POINTS-M0 W
0.99 0.98 0.97 0.8
0.6
ZERO-POINTS-M0 NWR
0.4
(c) M0 ZERO DGP
2.76 2.74
2.00
1.98
ZERO-POINTS-MY W
0.99 0.98 0.97
0.8
0.6
ZERO-POINTS-MY NWR
0.4
(d) MY ZERO DGP
Figure 34: Comparison of parameterization across different models and parameter initialization 101
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
ELL (x1e-6)
KLD (x1e-6)
ELL (x1e-6)
0.0 0.15
−0.1
KLD (x1e-6)
0.150
0.0
0.125
−0.1
0.100
−0.2 −0.2
0.10
−0.3 −0.4
0.075
−0.3
0.050
−0.4
0.025
0.05
−0.5 −0.6 0
50
100
150
200
0.00
ZERO-POINTS-M0 W
0
50
100
150
200
0
50
ZERO-POINTS-M0 NWR
100
150
200
PCA W
0.000
0
50
100
150
200
PCA NWR
Figure 35: Learning curves of both parameterizations, for each split, for the 5 depth layer DGPs trained in the Power dataset. Models displayed represent the PCA and the M0 models. Table 12: Average split-wise difference (in absolute value) between the performance of the whitened and NWR parameterizations of each of the models. Lower values indicate higher similarity between the performance of both parameterizations.
PCA ZERO ZERO-Points-M0 ZERO-Points-MY
Boston
Concrete
Energy
Kin8nm
Power
Protein
Redwine
Yacht
Overall
0.053 0.863 0.031 0.051
0.041 0.023 0.020 0.026
0.213 0.025 0.026 0.057
0.011 0.024 0.007 0.020
0.021 0.607 0.012 0.021
0.016 0.017 0.017 0.015
0.009 0.008 0.004 0.006
0.227 2.280 0.146 0.127
0.074 0.481 0.033 0.040
The proposed initializations, M0 and MY, do not mitigate the noisy optimization process itself, as expected, but they nevertheless help avoid convergence to poor solutions. Lastly, in Table 12, we show the average difference between the performance of the W and NWR versions of each model. This table reveals that the proposed initializations yield the lowest average difference between the two different parameterizations of the model. This supports the previously indicated idea that our proposal induces a much more stable optimization process, which critically leads the model to avoid the posterior collapse pathology.
C.5.3 Split-wise Differences In the experimental section, we observed that the proposed initialization strategy sometimes outperforms or simply matches the performance of the ZERO model. As we mentioned, since the statistical model or the inference family is not changed, we just shall expect improved performance of the proposed strategy when the ZERO model with the initialization m = 0 cannot be correctly optimized. This section confirms this is the case most of the time, and that the reason behind suboptimal performance comes from a tendency towards posterior collapse. We now investigate the performance gap. Average results in Fig.18a,18b might hide particular splits in which the proposed solution solves the collapse problem. Only datasets with a large number of train-test splits collapsing, as in Yacht, will show an average degraded performance. This is because we are not expecting an improved performance of our proposal when the ZERO-W model optimizes correctly. Thus, for some datasets, we take the split with the lowest and the highest test log likelihood difference between the PCA and ZERO models and plot the corresponding KLD and likelihood variance in Fig. 36. We can observe the same tendency discussed. Models with the highest 102
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
performance difference (left column) usually present a higher likelihood variance, and a smaller KLD for the ZERO-W and ZERO-NWR model’s (blue line), showing this tendency to collapse. We observe how the proposed initialization (green line) solves the problem, presenting a behavior similar to the PCA (orange model). When the performance gap is the lowest (right column), we observe how the ZERO model matches the PCA and M0. Interestingly, we observe how, at the beginning, the variance of the ZERO model rapidly increases. At the end, however, the optimization algorithm is able to reduce it to a level comparable with the M0 and PCA models. Importantly, in the Concrete dataset, using 3 layer models, the ZERO-W works better than the proposed solution, see Fig.18a. This confirms that, effectively, we shall only expect improved performance of our approach when the model presents collapse. This is because both statistical models are the same. Regarding the rest of the datasets and models, we have widely inspected most of the splits. We observe that not all of them present this behavior, and that sometimes our proposed solution cannot correct the problem or directly works worse. This might be because we are only initializing the ZERO model close to the PCA model at the inducing point locations. This is a limitation of our proposed solution, with further research going in this direction. 325
KL Divergence
KL Divergence
Lik. Variance 0.10
Lik. Variance 0.10
280
300
0.09 270
275
0.08
0.08
250 0.06
225 19960
19980
PCA W (-3.16)
20000
0
ZERO W (-3.30)
5000
260
0.07
250
0.06
10000 15000 20000
19960
ZERO-POINTS-M0 W (-3.25)
19980
PCA W (-3.07)
0.05
20000
ZERO W (-3.07)
0
5000
10000 15000 20000
ZERO-POINTS-M0 W (-3.06)
Concrete dataset. 3 layers DGP. Whitened parameterization. KL Divergence
KL Divergence
Lik. Variance
Lik. Variance
0.065
800
500
0.060
700
0.065 0.060
0.055 600
480 0.055
0.050 460
500
0.050
0.045 19960
19980
PCA W (-2.61)
20000
0
ZERO W (-2.67)
5000
10000 15000 20000
19960
ZERO-POINTS-M0 W (-2.65)
19980
PCA W (-2.77)
20000
0
ZERO W (-2.77)
5000
10000 15000 20000
ZERO-POINTS-M0 W (-2.79)
Power dataset. 3 layers DGP. Whitened parameterization. KL Divergence
KL Divergence
Lik. Variance
1400
0.10
1350
0.08
19960
19980
PCA W (1.38)
20000
ZERO W (1.30)
0.14
1375
0.06
1300
0.12
1350
0.10
1325
0.08
1300 0
5000
10000 15000 20000
ZERO-POINTS-M0 W (1.35)
Lik. Variance
0.16
1400
0.06 19960
19980
PCA W (1.36)
20000
ZERO W (1.36)
0
5000
Kin8nm dataset. 5 layers DGP. Whitened parameterization.
103
10000 15000 20000
ZERO-POINTS-M0 W (1.36)
An Analysis of Posterior Collapse, Parameterization and Initialization in Variational DGPs
KL Divergence
KL Divergence
Lik. Variance
1000
0.150
800
Lik. Variance 0.09
580
0.125
0.08 560
0.100
600
0.075 400
0.07 0.06
540
0.050 19960
PCA NWR (-2.69)
19980
20000
0.05 0
ZERO NWR (-2.77)
5000
10000 15000 20000
ZERO-POINTS-M0 NWR (-2.73)
19960
PCA NWR (-2.77)
19980
20000
0
ZERO NWR (-2.77)
5000
10000 15000 20000
ZERO-POINTS-M0 NWR (-2.76)
Power dataset. 3 layers DGP. Non-whitened parameterization. 200
KL Divergence
KL Divergence
Lik. Variance 1.0
0.5
180
0.8
150
0.4
0.6
100
0.3
160
0.4
0.2
50 0.2 0 19950 19960 19970 19980 19990 20000
PCA NWR (-2.28)
Lik. Variance
0.6
140 0
ZERO NWR (-3.75)
5000
10000 15000 20000
ZERO-POINTS-M0 NWR (-2.22)
0.1
19950 19960 19970 19980 19990 20000
PCA NWR (-2.25)
0
ZERO NWR (-2.22)
5000
10000 15000 20000
ZERO-POINTS-M0 NWR (-2.30)
Boston dataset. 5 layers DGP. Non-whitened parameterization. KL Divergence
KL Divergence
Lik. Variance 1.0
600
0.8 400 200
0.8
650
0.6
625
0.6
0.4
600
0.4
0.2
575
0.2
0 0.0 19950 19960 19970 19980 19990 20000 0
PCA NWR (-2.66)
Lik. Variance 1.0
675
ZERO NWR (-4.28)
5000
10000 15000 20000
ZERO-POINTS-M0 NWR (-2.69)
19950 19960 19970 19980 19990 20000
PCA NWR (-2.74)
0.0
ZERO NWR (-2.76)
0
5000
10000 15000 20000
ZERO-POINTS-M0 NWR (-2.77)
Power dataset. 5 layers DGP. Non-whitened parameterization. (a) Greatest split performance difference.
(b) Lowest split performance difference.
Figure 36: Final KLD and likelihood variance parameter across several datasets, model deepness, and parameterization. The left column represents the curves for models in which the difference between the PCA and ZERO model performance is the highest across all the splits. The right column represents the same, but when the model performance is the lowest.
104