Scalable Gaussian Process Regression Using Deep Neural Networks

1 Introduction

Gaussian Processes (GPs) are widely used in regression [Rasmussen and Williams, 2006; Seeger, 2004] and classification [Nickisch and Rasmussen, 2008; Naish-Guzman and Holden, 2007] for their flexibility and competitive learnability. A GP regression model applies a feature-mapping function to project inputs into a feature space and then performs Bayesian linear regression [Rasmussen and Williams, 2006], as shown in Figure 1 (a). A kernel function is usually defined to replace the inner product of the feature-mapping function so as to hold the condition that the feature space is infinite-dimensional. However, such a kernel method leads the training of GPs to suffer from computational issue for large datasets due to the unfavourable inversion of kernel matrices that scales cubically with the number of data points.

There are currently two kinds of methods to address the computational issue of GPs. One is to construct low-rank approximations of kernel matrices, namely Sparse Gaussian Processes (SPGPs) [Smola and Bartlett, 2000; Lawrence et al., 2002; Snelson and Ghahramani, 2005]. The state-ofthe-art SPGP was proposed in [Snelson and Ghahramani, 2005], which was later renamed as Fully Independent Training Conditional (FITC) model in [Qui˜nonero-Candela and Rasmussen, 2005]. The performance of FITC is nearly comparable to that of full GP while its training complexity only scales linearly with the number of training samples. The other one is to learn feature-mapping functions explicitly by fixing the feature space to be finite-dimensional, such as the work of [L´azaro-Gredilla and Figueiras-Vidal, 2010]. In [L´azaroGredilla and Figueiras-Vidal, 2010], the authors used the hidden unit transfer functions of one-hidden-layer Neural Networks (NNs) as feature-mapping functions. The resulting model, Marginalized NN mixture (mixMNN 1 ), performs better than FITC at almost the same computational cost on several low-dimensional datasets.

∗This work was jointly supported by the National Science Funds of China for Major Instrument Project (Grant No: 61327809) and Major International Cooperation Project (Grant No: 61210013).

The idea of applying the hidden unit transfer functions as the feature-mapping functions of GPs has been explored in many literatures [Neal, 1995; Williams and Barber, 1998; Rasmussen and Williams, 2006]. The famous work was performed by Neal [Neal, 1995]. He proved that a one-hiddenlayer neural network equipped with infinite number of hidden units becomes a nonparametric GP model by placing independent zero-mean Gaussian priors to all the weights of the network and integrating them out. Nevertheless, this kind of GPs still suffer from the computational issue due to the assumption of the infinite hidden units. The mixMNN method alleviates this problem by fixing the number of the hidden units to be finite and only marginalizing out the output weights. The resulting model is a parametric GP, of which the kernel function is parameterized by the hidden unit transfer function.

However, in practical, the shallow architecture of the mixMNN makes it inefficient for the highly-nonlinear data, e.g. the image-based regression tasks [Zhou et al., 2005]. In this context, we propose a more flexible model by applying Deep Neural Networks (DNNs) as the feature-mapping functions of GPs (as shown in Figure 1 (c)). We perform a Bayesian linear regression on the top layer of the DNN, thus resulting in a new GP model, which is named as Deep-Neural-Network-based Gaussian Process (DNN-GP). DNNs such as Stacked Denoising Auto-Encoders (SDAEs) [Vincent et al., 2008; Erhan et al., 2010], deep belief networks [Hinton et al., 2006; Hinton and Salakhutdinov, 2006; Lee et al., 2008], and deep Boltzmann machines [Salakhutdinov and Hinton, 2009; Cho et al., 2013] have been widely adopted in machine learning and computer vision. They are capable of decomposing the data into regular patterns of multiple levels, in which higher level representations are the abstractions of lower level ones [Bengio, 2009]. Employing DNNs as the feature-mapping functions of GPs is expected to learn much more meaningful representation of the input although such feature map is finite-dimensional.

1 The method is denoted as MNNmix in [L´azaro-Gredilla and Figueiras-Vidal, 2010]. We rename it as mixMNN to be consistent with the denotations of our methods.

[Salakhutdinov and Hinton, 2007] also employed deep networks as feature extractors to improve the performance of GP regression. They applied the output of deep belief networks as the input of Gaussian kernels. which is substantially different from our algorithmic framework. Therefore, their model also needs to invert the kernel matrix, still suffering from the computational limitation as general GP models do. Our model takes advantage of deep networks by embedding them in the feature space, thus leading to the considerable reduction of computational complexity. [Snoek et al., 2011] and [Damianou and Lawrence, 2012] apply GP-based building blocks as the feature extractors for highly-nonlinear regression. Unlike our DNN-based models, such GP-based building blocks are not scalable to large data unless some extra techniques are added for computational reduction.

In sum, we make the following contributions: (1) We formulate a new kind of GPs, i.e. DNN-GPs, by applying DNNs as the feature-mapping functions. We fix the dimensionality of the feature space to be finite and learn the featuremapping functions explicitly. Therefore, our GPs are feasible for large datasets since the computational complexity only scales linearly with the size of the training set. (2) In DNN-GPs, we pre-train DNN-GPs with SDAEs in an unsupervised way, which not only provides better initializations to the parameters of DNNs, but also makes DNN-GPs perform competitively when the labeled data are limited. (3) We propose mixture models by combining the predictions of several DNN-GPs, dubbed mixDNN-GPs. The mixDNN-GPs algorithm significantly outperforms all compared models in our experiments.

2 Gaussian Process Regression

The most common way to interpret Gaussian processes regression is defining a kernel function as the covariance of the distribution over latent functions, which is known as the function-space view [Rasmussen and Williams, 2006]. In this paper, however, we start from the equivalent weightspace view [Rasmussen and Williams, 2006] that is more convenient to establish our model. Readers can resort to Mercer's theorem for the proof of the equivalence between the function-space view and the weight-space view in [Rasmussen and Williams, 2006]. Suppose that we have a dataset D of N observations {(xn, yn)}Nn=1 , where the dimensionality of each input is D, i.e. xn ∈ RD . We denote X =

[x1, · · · , xN] as the input matrix and y = [y1, · · · , yN] as the output vector. For a regression task, we assume that each output is generated from a probabilistic model

Specifying the latent function f(xn) to be the linear projection wT xn, we obtain the Bayesian regression model, where the vector of weights w is sampled from a Gaussian distribution N (0, Σp). A bias weight or offset can be derived by augmenting the input xn with an additional element. We do not explicitly include it in our notation for simplicity. If we make a further extension by first projecting the inputs into a feature space with a mapping function φ(x), we formulate a Gaussian process regression model, as illustrated in Figure 1 (a). It means that in a Gaussian process, the latent function is defined as

After marginalizing out the weight vector w, we obtain the prior

where the latent functions f = [f(x1), · · · , f(xN)] and K is the covariance matrix whose entries are specified by Kij = φ(xi)T Σpφ(xj). To improve the model representation ability in standard GP models, the feature space is always assumed to be infinite-dimensional [Rasmussen and Williams, 2006], thereby leading to the application of a kernel function k(xi, xj) as the entries of the covariance matrix to circumvent the difficulty of computing φ(xi)T Σpφ(xj). The kernel function depends on a small number parameters θ. For example, the Automatic Relevance Determination (ARD) kernel is defined as k(xi, xj) = α exp(- 1 2 ∑Dd=1 (x(d) i - x(d) j )2/β2 d ), where θ = {α, β1, · · · , βD}.

Integrating out latent functions, the negative marginal loglikelyhood takes the form

which is a minimization objective to train the parameters θ and σ. The prediction of a new input data x∗ can be made by conditioning on observed data and parameters. The distribution of a new target y∗ at the new input x∗ is

where k∗ is a vector with elements k(xi, x∗) and k∗∗ = k(x∗, x∗).

Equations (4-5) indicate that the key computation in GPs including the training and the prediction is to invert the covariance matrix K + σ2 I, which is of complexity O(N 3 ).

3 FITC and mixMNN

To reduce the computational time of GPs from cubic scale to linear scale with respect to N, one approach resorts to a small set of inducing points to construct a low-rank approximation to the kernel matrix K over the full dataset , which is known as Sparse GP [Smola and Bartlett, 2000; Lawrence et al., 2002; Snelson and Ghahramani, 2005]. The most successful Sparse GP is the method proposed in [Snelson and Ghahramani, 2005], and then renamed as FITC in [Qui˜nonero-Candela and Rasmussen, 2005]. FITC augments the training set D by adding a noiseless pseudo-data set ¯D of size M < N: pseudo-inputs ¯X = {¯xn}Mn=1 and pseudolatent-functions ¯f = {¯fn}Mn=1 . Assuming fns to be conditionally independent given the pseudo-data set, we attain the prior over the latent variables f by integrating out the pseudolatent-functions:

where Q = Kf¯f K-1 ¯f¯f K¯ff , Kf¯f , K¯f¯f and K¯ff are matrices with the elements k(xi, ¯xj), k(¯xi, ¯xj) and k(¯xi, xj), respectively. Compared to the original covariance matrix K in Equation (3), the covariance matrix here is a low-rank matrix Q plus a diagonal matrix diag(K - Q). With the new covariance matrix, Equations (4-5) can be computed in O(NM 2 ) time, using the Woodbury matrix identity [Rasmussen and Williams, 2006].

Another way to reduce the computational cost is to learn the feature-mapping function explicitly. Denoting φ(X) = [φ(x1), · · · ,φ(xN)], the covariance matrix K is equal to φ(X)T Σpφ(X). Setting the dimensionality of φ(X) to be M < N, the term φ(X)T Σpφ(X) becomes a lowrank construction of K. Explicitly substituting this term into Equations (4-5), the computational cost can be reduced in O(NM 2 ) as well. One example is the neural-networkbased models mixMNNs [L´azaro-Gredilla and FigueirasVidal, 2010], in which the hidden unit transfer functions of one-hidden-layer neural networks are employed as featuremapping functions. Specifically, mixMNNs assign φ(x) = g(Ux + u0), where U and u0 are the input weight matrix and input bias vector, respectively, and g(·) is the activation function.

4 Our Model: DNN-GP

Similar to mixMNNs, we also reduce the computational cost of GPs by explicitly learning the finite-dimensional featuremapping functions. Our models can be learned in two steps: the unsupervised pre-training (Section 4.1) and the supervised fine-tuning (Section 4.2).

4.1 Model Pre-training

Stacked Denoising Auto-Encoders (SDAEs) specify deep networks by applying a Denoising Auto-Encoder (DAE) to initialize each layer of deep networks [Vincent et al., 2008]. In a DAE, the original input data X are corrupted by some kind of stochastic noise. In this paper, we follow the corruption process used in [Vincent et al., 2008]. For the input X, a fixed proportion v of the dimensions of each sample are chosen randomly to be forced to 0, while other dimensions are left invariant, thus leading to the corrupted ones X′ . X′ are used for the input while X are used for the reconstruction targets. The parameters Wi, bi and ci are trained to minimize the squared error loss function

where ‖ · ‖F is the matrix Frobenius norm, Wi is a weight matrix , bi and ci are bias vectors, g(·) is a nonlinear active function, e.g. a sigmoid function.

Once an DAE is trained, its internal representation g(WiX′ +bi) can be adopted as the input for training a second DAE. Figure 1 (b) illustrates this training procedure. After such greedy layer-by-layer training, we will obtain a top layer as the highly-nonlinear representation of the input, i.e. ϕ(X), which generally captures high-order patterns in original data. After training the whole architecture, we initialize a DNN with the parameters of the learned SDAE. SDAEs can be learned efficiently on large datasets. Training SDAEs scales linearly with the number of training samples. In addition, Equation (7) can be written as the summation of data points, so we can use a stochastic optimization procedure on the mini-batches of data points for each iteration.

4.2 Model Fine-tuning

Substituting the new feature-mapping function into Equation (2), we attain a GP with the elements of the covariance matrix K(xi, xj) = ϕ(xi)T Σpϕ(xj), as illustrated in Figure1 (c).

For simplicity, we consider the isotropic covariance for the weights, namely Σp = s-2 p I, where s2p is the precision parameter. Thus, K(xi, xj) = s-2 p ϕ(xi)T ϕ(xj). Substituting the new covariance matrix into Equation (4) and applying the Woodbury matrix identity (A.9 and A.10 in [Rasmussen and Williams, 2006]) to the terms (K+ σ2 I)-1 and |K+ σ2 I|, we attain

where F = [ϕ(x1), · · · ,ϕ(xN)] ∈ RM×N and A = σ2 s2p I + FFT ∈ RM×M . M is the number of the units of the top layer of the SDAE. The computational complexity of Equation (8) is O(NM 2 ).

Substituting the new covariance matrix into the prediction Equation (5), we have

The computational complexity of the new prediction is only O(NM). 2 Denoting h(x∗) = L-1ϕ(x∗) as the new feature-mapping function and letting s = L-1Fy, we can compactly rewrite Equation (9) as

where the new predictive mean is exactly a linear combination of the new features. It means that in terms of the prediction of the expectation, we finally derive a feature space where the output is the linear function of the input, after a highly-nonlinear transformation, i.e. h(x∗). The vector s depends on the training samples but is independent of the testing input. It can measure the relevances of the features to the prediction.

Now we perform the discriminative fine-tuning using back propagation. We first compute the derivatives of the objection function (8) with respect to the explicit parameters including the mapping matrix F and parameters σ2 and s2p

2 Equation (8) and Equation (9) are similar to Equation (7) and Equation (6) in [L´azaro-Gredilla and Figueiras-Vidal, 2010] respectively. The difference is that ϕ(X) in [L´azaro-Gredilla and Figueiras-Vidal, 2010] is only a simple hidden unit transfer function of one-hidden-layer NN, while ϕ(X) here is the output of the learned SDAEs.

where the matrix C = 2 × ∂L ∂A , and tr(·) is the matrix trace. Then by the chain rule, we readily back propagate ∂L ∂F to the implicit parameters W and b. For W , we compute the gradients as∂L

where {a2, · · · , ad} are the active values of the hidden units of the corresponded layers, and d is the depth of DNN-GPs, i.e. the number of the layers of DNNs. It is worth noting that we do not include the output layer when counting d, as already shown in Figure 1 (c). The computation of the gradients of b is similar to Equation (14). One may refer to Section 5.3 in [Bishop, 2006] for details. Similar to pre-training, backpropagation in Equation (14) scales linearly with N [Hinton and Salakhutdinov, 2006]. The computational complexity of Equation (8) is O(NM 2 ). Thus, our model is scalable to the number of training samples .

Based on Equations (11)-(14), we apply the conjugate gradient method to fine-tune our models. Because the objective function (8) cannot be decomposed into a sum over data points, it encounters the difficulty when using a stochastic optimization procedure on the mini-batches of data points. Thus, we apply the full-batch optimization for fine-tuning. Even so, fine-tuning our models on full-batch performs with the competitive speed of convergence in our experiments.

4.3 Combining predictions of several DNN-GPs

To overcome the risk of overfitting, mixMNNs average individual predictions of several one-hidden-layer NNs of which parameters are initialized randomly. In DNN-GPs, we initialize the parameters of DNNs at random and then pre-train them with SDAEs. In this way, a single DNN-GP is capable of achieving satisfactory performance. However, we find that combining predictions of several DNN-GPs can further dramatically boost the performance. Following [L´azaro-Gredilla and Figueiras-Vidal, 2010], we derive the combined mean and variance of the distribution of y∗ at a new test sample x∗

where K is the number of combined GPs, µ∗k and σ2 ∗k are the predictive mean and variance of k-th GP, respectively. We denote the mixture version of DNN-GPs as mixDNN-GPs.

5 Experiments

We compare the performance of proposed models with full GPs, FITCs [Snelson and Ghahramani, 2005] and mixMNNs

[L´azaro-Gredilla and Figueiras-Vidal, 2010] on one synthetic image dataset Rectangles [Larochelle et al., 2007] and two real face datasets Olivetti Faces [Salakhutdinov and Hinton, 2007] and FGnet [Zhou et al., 2005]. Rectangles is used for predicting the area of rectangles, Olivetti Faces is for face orientation extraction and FGnet is for age estimation. We first enlarge each dataset to contain more than 10000 samples, and then pre-process all the images to the same size 32 × 32 and normalize the values of pixels in [0, 1]. For reader's convenience, we list the information of datasets in Table 1.

For mixMNNs and DNN-GPs, we use the pre-processed images as the inputs and the sigmoid functions as the activation functions. For full GPs and FITCs, PCA is applied to the images for feature extraction so as to improve the performance, meaning that the first p principal components are chosen as the input. To select the optimal p, the performance of a full GP will be evaluated for various values of p and the value resulted in best performance will be chosen. The setting of training SDAEs in DNN-GPs is fixed as follows: The corrupted proportion v in each DAE is set to be 0.5. Equation (7) is optimized with the gradient descent method whose learning rate is 0.1. The dataset is divided into mini-batches of size 100, and the parameters are updated after each minibatch iteration. Each DAE is trained for 30 passes through the entire training set. We train full GPs, FITCs, mixMNNs, and fine-tune DNN-GPs with the conjugate gradient method [Rasmussen and Williams, 2006] on the full-batch of training samples. If the objective function dose not decrease within 50 evaluations in each line search, or the epoch of the line search exceeds 100, the training process will halts.

Following [L´azaro-Gredilla and Figueiras-Vidal, 2010; Titsias and L´azaro-Gredilla, 2013], all the experiments here are measured with Normalized Mean Square Error (NMSE)

where yn and µn are the true label and predictive mean of the n-th testing sample, respectively, ¯ y is the mean of the labels of the training set and Ntest is the number of testing samples.

5.1 Notations: M and d

We use a unified symbol M to denote the number of the units at the top layer of SDAEs in DNN-GPs, the number of pseudo-inputs in FITCs and the number of hidden units in mixMNNs. In all of our experiments, we fix the number of the neural units of the other hidden layers to be 500 excluding the top layer of SDAEs. Such a constraint may limit the ability of DNN-GPs. In our practice, however, we find that the DNN-GPs of this architecture are sufficient to achieve desirable performance. We denote K as the number of the individual architectures used to construct the mixture models (i.e. mixMNNs and mixDNN-GPs). Increasing K leads to enhanced performance and higher computational cost. We set K = 4 for both mixMNNs and mixDNN-GPs.

5.2 Analysis on ProposedModels

We compare the performance of DNN-GP with its several variants to evaluate the influence of the depth d, the pretraining, the fine-tuning and the mixture setting. Here, two sub-experiments are implemented.

In the first one, we define a simplified version of DNNGPs as simDNN-GPs, in which no pre-training procedure is involved. The parameters of simDNN-GPs are initialized at random. We report both the training and testing NMSEs on Olivetti Face in Table 2. Here, the training NMSE and the testing NMSE mean the NMSE measured on the training set and the testing set, respectively. Table 2 shows that the testing NMSEs of all the three models will decrease when d increases. Building a deep architecture is helpful for improving the performance of GP regression. With the same d, simDNN-GPs provide lower training NMSEs but higher testing NMSEs than DNN-GPs. The unsupervised pre-training acts as a regularizer to enforce the parameters of DNN-GP to a constrained region, where a better generalization yields. The mixDNN-GPs approach provides both lower testing NMSEs and training NMSEs than DNN-GPs. Committing a mixture of DNN-GPs is able to attain an optimized solution, thereby resulting in an improved performance.

In the second one, we compare the performance of DNNGP with its two counterparts: LBR+DNN and GP+DNN. LBR+DNNs and GP+DNNs apply the output of SDAEs as the input of Linear-Bayesian-Regression models (LBRs) and full GPs, respectively. Unlike DNN-GPs, there is no global fine-tuning involved in this two models. The average testing NMSEs are reported in Tabel 3. Clearly, DNN-GPs achieve the best accuracy. LBR+DNNs perform grossly worse than DNN-GPs, which verifies that the global fine-tuning plays the essential role in our models.

5.3 Comparison on Full Training Sets

In this experiment, we perform the comparison between FITC+PCAs, mixMNNs, 5-layer DNN-GPs and 5-layer mixDNN-GPs. To evaluate the performance of above models at different computational complexities, we vary M from 25 to 800 and display the results of the testing NMSEs in Figure 2. We also implement fullGP+PCA on the full training set so as to provide a reference of the desirable performance. From Figure 2, we observe that a larger M results in a better performance for both DNN-GPs and mixDNN-GPs. DNN-GPs perform better than FITC+PCAs and mixMNNs on all three datasets when M > 50. Compared with fullGP+PCA, DNNGPs obtain lower testing errors on Rectangles and Olivetti Faces, and the comparable accuracy on FGnet. mixDNNGPs achieve the best performance among all the compared models when M = 800.

The training of DNN-GPs includes the layer-by-layer pretraining (Equation (7)), the computation of the derivatives of the explicit parameters (Equations (11) - (13)) and the back propagation of the weight parameters (Equation (14)). As we have demonstrated in Section 4, all these computations scale linearly with the number of training samples, which leads DNN-GPs to be feasible for large datasets. The computational cost for training the full GP and the DNN-GP is presented in Tabel 4. Obviously, the training of DNN-GPs is much faster than that of full GPs.

5.4 Comparison on Limited Labeled Sets

In this task, we evaluate the performance of FITC+PCAs, mixMNNs, 5-layer simDNN-GPs, 5-layer DNN-GPs and 5layer mixDNN-GPs on limited labeled training sets. M is fixed to be 800. We construct the labeled training set of size N by randomly selecting a subset from the original training set. Training FITC+PCAs, mixMNNs and simDNN-GPs are based on the labeled training set. For DNN-GPs and mixDNN-GPs, we first pre-train these two models with the original training set and then fine-tune them with the limited labeled set. Figure 3 illustrates the performance of the compared models on the labeled sets of different sizes.

On Olivetti Faces, DNN-GP achieves very low error even when the number of labeled data N is very small. Particularly, DNN-GP with only 2000 labeled training samples outperforms FIIC+PCAs and mixMNN which are trained on the full dataset, i.e. N = 16000. In contrast, simDNN-GP performs worse than FIIC+PCA and mixMNN when N < 8000. When N = 16000, simDNN-GP finally outperforms FIIC+PCA and mixMNN but still is inferior to DNN-GP and mixDNN-GP. On the other two datasets Rectangles and FGnet, we still observe the fact that DNN-GPs perform better than simDNN-GPs in all cases, thus verifying the importance of the unsupervised pre-training. As expected, mixDNN-GPs outperform all the compared models in all cases.

6 Conclusion

In this paper, we applied stacked denoising auto-encoders as the feature-mapping functions to reformulate GPs, which have been named as Deep-Neural-Network-based GPs (DNN-GPs). Our experiments verified the following conclusions: (1) DNN-GPs considerably outperform FITCs and mixMNNs, and perform better than full GPs as a whole, while the training of DNN-GPs is scalable to the size of the training set. (2) The unsupervised training of SDAEs is important to the performance of DNN-GPs, which not only provides better initializations to the parameters of DNNs, but also makes DNN-GPs perform competitively when the number of the labeled data is limited. (3) The mixture setting can significantly improve the performance of DNN-GPs. The mixDNN-GPs algorithm achieves the state-of-the-art performance in our experiments.

Overall, applying deep neural networks as the featuremapping functions is able to effectively facilitate the performance of GPs for regression problems with a large amount of training data. For the future work, we will extend DNN-GPs to multi-output regression problems.