Among many extensions of GPs for enhancing expressivity, e.g. warped GP in [1]}, Deep Gaussian Processes (DGPs) [2]} are a general hierarchical composition of GPs. The compositional structure enhances its expressive power, e.g. a GP with SE kernel can not fit a step function well but a DGP can. Consider for simplicity a two-layer function \(f({\bf x})=f_2({\bf f}_1({\bf x}))\)  where the input \({\bf x}\in \mathbb {R}^D\)  is mapped to the hidden output \({\bf h=f}_1({\bf x})\in \mathbb {R}^H\)  and then to a real output \(f_2({\bf h})\) . The hidden layer with finite \(H\)  is referred to as the bottleneck in [3]}, [4]}. DGP is defined by the joint density \(p(f_2({\bf f}_1({\bf x}_{1:N}))\) ,
\(\mathcal {N}(f_2({\bf h}_{1:N})|0,\Sigma _2({\bf H}))\prod _{i=1}^H\mathcal {N}(h_i({\bf x}_{1:N})|0,\Sigma _1({\bf X}))\:,\) 
