Title: GAUSSIAN PROCESS NEURONS

Abstract: We propose a method to learn stochastic activation functions for use in probabilistic neural networks. First, we develop a framework to embed stochastic activation functions based on Gaussian processes in probabilistic neural networks. Second, we analytically derive expressions for the propagation of means and covariances in such a network, thus allowing for an efficient implementation and training without the need for sampling. Third, we show how to apply variational Bayesian inference to regularize and efficiently train this model. The resulting model can deal with uncertain inputs and implicitly provides an estimate of the confidence of its predictions. Like a conventional neural network it can scale to datasets of arbitrary size and be extended with convolutional and recurrent connections, if desired.

1. INTRODUCTION:

The popularity of deep learning and the implied race for better accuracy and performance has lead to new research of the fundamentals of neural networks. Finding an optimal architecture often focusses on a hyperparameter search over the network architecture, regularization parameters, and one of a few standard activation functions: tanh, ReLU (Glorot et al. (2011) ), maxout (Goodfellow et al. (2013) ) . . . Focussing on the latter, looking into activation functions has only taken off since Nair & Hinton (2010) introduced the rectified linear unit (ReLU), which were shown to produce significantly better results on image recognition tasks Krizhevsky et al. (2012) . Maas et al. (2013) then introduced the leaky ReLU, which has a very small, but non-zero, slope for negative values. He et al. (2015) proposed the parameterized ReLU, by making the slope of the negative part of the leaky ReLU adaptable. It was trained as an additional parameter for each neuron alongside the weights of the neural network using stochastic gradient descent. Thus, the activation function was not treated as a fixed hyper-parameter anymore but as adaptable to training data. While the parameterized ReLU only has one parameter, this was generalized in Agostinelli et al. (2014a) to piecewise linear activation functions that can have an arbitrary (but fixed) number of points where the function changes it slope. This can be interpreted as a different parameterization of a Maxout network (Goodfellow et al. (2013) ), in which each neuron takes the maximum over a set of different linear combinations of its inputs.

Instead of having a fixed parameter for the negative slope of the ReLU, Xu et al. (2015) introduced stochasticity into the activation function by sampling the value for the slope with each training iteration from a fixed uniform distribution. Clevert et al. (2015) and Klambauer et al. (2017) replaced the negative part of ReLUs with a scaled exponential function and showed that, under certain conditions, this leads to automatic renormalization of the inputs to the following layer and thereby simplifies the training of the neural networks, leading to an improvement in accuracy in various tasks.

Nearly fully adaptable activation functions have been proposed by Eisenach et al. (2016) . The authors use a Fourier basis expansion to represent the activation function; thus with enough coefficients any (periodic) activation function can be represented. The coefficients of this expansion are trained as network parameters using stochastic gradient descent or extensions thereof.

Promoting a more general approach, Agostinelli et al. (2014b) proposed to learn the activation functions alongside the layer weights. Their adaptive piecewise linear units consist of a sum of hinge-shaped functions with parameters to control the hinges and the slopes of the linear segments. However, by construction the derivative of these activation functions is not continuous at the joints between two linear segments, which often leads to non-optimal optimizer performance.

To our knowledge, previous research on learning activation functions took place in a fully deterministic setting, i.e. deterministic activation functions were parameterized and included in the optimization of a conventional neural network. Here instead, we explore the setting of probabilistic activation functions embedded in a graphical model of random variables resembling the structure of a neural network. We develop the theory of Gaussian-process neurons and subsequently derive a lower-bound approximation using variational inference, in order to develop a computationally efficient version of the Gaussian Process neuron.

NOTATION:

To define the model we will need to slice matrices along rows and columns. Given a matrix X, we will write X i to select all elements of the i-th row and X j to select all elements of the j-th column.

2. GAUSSIAN PROCESSES:

Gaussian Processes (GPs) are nonparametric models that provide flexible probabilistic approaches for function estimation. A Gaussian Process (Rasmussen & Williams, 2006) defines a distribution over a function f (x) ∼ GP(m(x), k(x, x )) where m is called the mean function and k is the covariance function. For S inputs x ∈ R S×N of N dimensions the corresponding function values

f with f i f (X i ) follow a multivariate normal distribution 1 f ∼ N (m, K(X, X))

with mean vector m i m(X i ) and K(X, X) is the covariance matrix defined by (K(X, X)) ij k(X i , X j ). In this work we use the zero mean function m(x) = 0 and the squared exponential (SE) covariance function with scalar inputs,

EQUATION

)

where ν is called the length-scale and determines how similar function values of nearby inputs are according to the GP distribution. Since this covariance function is infinitely differentiable, all function samples from a GP using it are smooth functions.

3. GAUSSIAN PROCESS NEURONS:

We will first describe the fundamental, non-parametric model, which will be approximated in the following sections for efficient training and inference. Let the input to the l-th layer of Gaussian Process neurons (GPNs) be denoted by X l-1 ∈ R S×N l-1 where S is the number of data points (samples) and N l-1 is the number of input dimensions. A layer l ∈ {1, . . . , L} of N l GPNs indexed by n ∈ {1, . . . N l } is defined by the joint probability

P(X l n , F l n | X l-1 ) = P(F l n | X l-1 ) P(X l n | F l n ) (2)

with the GP prior F l conditioned on the layer inputs X l-1 multiplied with the weights W l ,

EQUATION

and an additive Gaussian noise distribution,

EQUATION

This corresponds to a probabilistic activation function

F l sn = f l n (X l-1 s W l n ) with f l n (z) ∼ GP(0, k(z, z )) .

(5)

This GP has scalar inputs and uses the standard squared exponential covariance function. Analogous to standard neural networks, GPN layers can be stacked to form a multi-layer feed-forward network.

The joint probability of such a stack is

P(X {l} , F {l} ) = P(X 0 ) L l=1 P(X l n , F l n | X l-1 ) . (6) X l-1 sN . . . X l-1 s1 A l sn F l sn V l n U l n X l sn sample s W n σ 2 n inputs

activation response output:

Figure 1 : The auxiliary parametric representation of a GPN using virtual observation inducing points V and targets U .

All input samples X 0 s , s ∈ {1, . . . , S}, are assumed to be normally distributed with known mean and covariance,

EQUATION

To obtain predictions P(X L | X 0 ), all latent variables in eq. ( 6) would need to be marginalized out; unfortunately due to the occurrence of X l in the covariance matrix in eq. ( 3) analytic integration is intractable.2 

3.1. AUXILIARY PARAMETRIC REPRESENTATION:

The path to obtain a tractable training objective is to temporarily parameterize the activation function f l n (z) of each GPN using virtual observations (originally proposed by Quiñonero-Candela & Rasmussen (2005) for sparse approximations of GPs) of inputs and outputs of the function. These virtual observations are only introduced as an auxiliary device and will be marginalized out later. Each virtual observation r consists of a scalar inducing point V l rn and corresponding target U l rn . Under these assumptions on f l d , the GP prior P(F l n | X l-1 ) is replaced with

EQUATION

where the mean and variance are those obtained by using the virtual observations as "training" points for a GP regression evaluated at the "test" points X l W l d ,

EQUATION

)

where K l nrt = K(V l rn , V l tn ). Given enough inducing points that lie densely between the layer's activations A l n = X l-1 W l n , the shape of the activation function becomes predominantly determined by the corresponding targets of these inducing points. Consequently, the inter-sample correlation in eq. ( 8) becomes negligible, allowing us to further approximate this conditional by factorizing it over the samples; thus we have

EQUATION

We now marginalize eq. ( 4) over F l n and get P(

EQUATION

) (a) X 0 A 1 U 1 F 1 X 1 A 2 U 2 F 2 X 2 A 3 U 3 F 3 X 3 (b) X 0 A 1 U 1 F 1 X 1 A 2 U 2 F 2 X 2 A 3 U 3 F 3 X 3

Figure 2 : A GPN feed-forward network prior distribution (a) for three layers and the approximation of its posterior obtained using variational inference and the central limit theorem on the activations (b). Each node corresponds to all samples and GPN units within a layer. Dotted circles represent variational parameters.

Although we now have a distribution for X l that is conditionally normal given the values of the previous layer, the marginals P(X l n ), l ∈ {1, . . . , L}, will, in general, not be normally distributed, because the input from the previous layer X l-1 appears non-linearly through the kernel function in the mean eq. ( 9).

By putting a GP prior on the distribution of the virtual observation targets U l as shown in fig. 1 .,

EQUATION

it can easily be verified that the marginal distribution of the response,

EQUATION

)

recovers the original, non-parametric GPN response distribution given by (3). The first use of this prior-restoring technique was presented in Titsias ( 2009) for finding inducing points of sparse GP regression using variational methods.

3.2. APPROXIMATE BAYESIAN INFERENCE FOR GPN NETWORKS:

The joint distribution of a GPN feed-forward network is given by

EQUATION

)

where to notation

{•} L 1 should be read as {• 1 , • 2 , . . . , • L }.

A graphical model corresponding to that distribution for three layers is shown in fig. 2a . Since exact marginalization over the latent variables is infeasible we apply the technique of variational inference Wainwright et al. (2008) to approximate the posterior of the model given some training data by a variational distribution Q.

The information about the activation functions learned from the training data is mediated via the virtual observation targets U l , thus their variational posterior must be adaptable in order to store that information. Hence, we choose a normal distribution factorized over the GPN units within a layer with free mean and covariance for the approximative posterior of U l ,

EQUATION

This allows the inducing targets of a GPN to be correlated, but the covariance matrix can be constrained to be diagonal, if it is desired to reduce the number of variational parameters. We keep the with β l n K(V l n , V l n ) -1 µ U l n and, using the same method as above,

EQUATION

For off-diagonal elements, n = n , we observe that F l sn and F l sn are conditionally independent given A l because the activation functions of GPNs n and n are represented by two different GPs. Hence we have

EQUATION

where

Λ l srtnn E Q(A l s ) K(A l sn , V l tn ) K(A l sn , V l rn ) = ν 2 exp(A l srtnn /B l snn ) [ν 2 + Σ A l snn ] [ν 2 + Σ A l sn n ] -(Σ A l nn ) 2 with A l srtnn (V l rn -µ A l sn ) 2 [ν 2 + Σ A l snn ] + (V l tn -µ A l sn ) 2 [ν 2 + Σ A l sn n ] + 2 (V l tn -µ A l sn ) (µ A l sn -V l rn ) Σ A l snn B l snn 2 [ν 2 + Σ A l snn ] [ν 2 + Σ A l sn n -(Σ A l sn n ) 2 .

This concludes the calculation of the moments of F l . We can now state how the activation distribution propagates from a layer to the next. The marginal distribution of A l+1 is given by

EQUATION

)

Thus Q(A L ) in ( 31) can be calculated by iterating the application of eqs. ( 33), ( 35), ( 37), ( 39) and ( 40) over the layers l. To save computational power only the variances can be propagated by assuming that Σ F l s is diagonal and therefore ignoring (37).

Now Ł pred can be identified as the expected log-probability of the observations under the marginal distribution Q(F L ) and thus we can expand it as follows,

Ł pred = E Q(F L ) log P(X L | F L ) ∝ E Q(F L ) -S N l n=1 log σ L n - 1 2 S s=1 N L n=1 X L sn -F L sn 2 (σ L n ) 2 = -S N L n=1 log σ L n - 1 2 S s=1 N L n=1 (X L sn ) 2 -2 X L sn E Q(F L ) F L sn + E Q(F L ) (F L sn ) 2 (σ L n ) 2 (41)

where S is the number of training samples and X L are the training targets. The distribution Q(F L ) is of arbitrary form, but only its first and second moments are required to evaluate Ł pred . For the first moment we obtain E Q(F L ) F l sn = µ F L sn with µ F L sn given by ( 33) and the second moment evaluates to

EQUATION

with

β L n K(V L n , V L n ) -1 µ U L n

and Ω L from (36). This concludes the calculation of all terms of the variational lower bound (18). The resulting objective is a fully deterministic function of the parameters. Training of the model is performed by maximizing Ł w.r.t. to the variational parameters µ U l , Σ U l and the model parameters σ l , W l and V l . This can be performed using any gradient-descent based algorithm. The necessary derivatives are not derived here and it is assumed that this can be performed automatically using symbolic or automatic differentiation in an appropriate framework.

3.3. COMPUTATIONAL AND MODEL COMPLEXITY:

The activation function are represented using 2R variational parameters per GPN, where R is the number of inducing points and targets. It can be shown that R = 10 linearly spaced inducing points are enough to represent the most commonly used activation functions (sigmoid, tanh, soft ReLU) with very high accuracy. The number of required parameters can be reduced by sharing the same activation function within groups of neurons or even across whole layers of neurons. If the inducing points V l are fixed (for example by equally distributing them in the interval [-1, 1]), the kernel matrices K(V l , V l ) and their inverses can be precomputed since they are constant. The number of parameters and the computational complexity of propagating the means and covariances only depend on R and are therefore independent of the number of training samples. Thus, like a conventional neural network, a GPN network can inherently be trained on datasets of unlimited size.

4. CONCLUSION:

We have presented a non-parametric model based on GPs for learning of activation functions in a multi-layer neural network. We then successively applied variational to make fully Bayesian inference feasible and efficient while keeping its probabilistic nature and providing not only best guess predictions but also confidence estimations in our predictions. Although we employ GPs, our parametric approximation allows our model to scale to datasets of unlimited size like conventional neural networks do.

We have validated networks of Gaussian Process Neurons in a set of experiments, the details of which we submit in a subsequent publication. In those experiments, our model shows to be significantly less prone to overfitting than a traditional feed-forward network of same size, despite having more parameters.

:

We omit writing P (•).

The inverse of the covariance matrix appears in the PDF of a normal distribution. Thus, the dependency on X l is highly non-linear and analytic calculation of the integral is not feasible.

Figures and Tables:

Peer review (iclr2018):
 In Bayesian neural networks, a deterministic or parametric activation is typically used. In this work, activation functions are considered random functions with a GP prior and are inferred from data. - Unnecessary complexity The presentation of the paper is unnecessarily complex. It seems that authors spend extra space creating problems and then solving them. Although some of the derivations in Section 3.2.2 are a bit involved, most of the derivations up to that point (which is already in page 6) follow preexisting literature. For instance, eq. (3) proposes one model for p(F|X). Eq. (8) proposes a different model for p(F|X), which is an approximation to the previous one. Instead, the second model could have been proposed directly, with the appropriate citation from the literature, since it isn*t new. Eq. (13) is introduced as a "solution" to a non-existent problem, because the virtual observations are drawn from the same prior as the real ones, so it is not that we are "coming up" with a convenient GP prior that turns out to produce a computationally tractable solution, we are just using the prior on the observations consistently. In general, the authors seem to use "approximately equal" and "equal" interchangeably, which is incorrect. There should be a single definition for p(F|X). And there should be a single definition for L_pred. The expression for L_pred given in eq. (20) (exact) and eq. (41) (approximate) do not match and yet both are connected with an equality (or proportionality), which they shouldn*t. Q(A) is sometimes taken to mean the true posterior (i.e., eq. (31)), sometimes a Gaussian approximation (i.e., eq (32) inside the integral), and both are used interchangeably. - Incorrect references to the literature Page 3: "using virtual observations (originally proposed by Quiñonero-Candela & Rasmussen (2005) for sparse approximations of GPs)" The authors are citing as the origin of virtual observations a survey paper on the topic. Of course, that survey paper correctly attributes the origin to [1]. Page 4: "we apply the technique of variational inference Wainwright et al. (2008)". How can variational inference be attributed to (again) a survey paper on the topic from 2008, when for instance [2] appeared in 2003? - Correctness of the approach Can the authors guarantee that the variational bound that they are introducing (as defined in eqs. (19) and (41)) is actually a variational bound? It seems to me that the approximations made to Q(A) to propagate the uncertainty are breaking the bounding guarantee. If it is no longer a lower bound, what is the rationale behind maximizing it? The mathematical basis for this paper is actually introduced in [3] and a single-layer version of the current model is developed in [4]. However, in [4] the authors manage to avoid the additional Q(A) approximation that breaks the variational bound. The authors should contrast their approach with [4] and discuss if and why that additional central limit theorem application is necessary. - No experiments The use of a non-parametric definition for the activation function should be contrasted with the use of a parametric one. With enough data, both might produce similar results. And the parameter sharing in the parametric one might actually be beneficial. With no experiments at all showing the benefit of this proposal, this paper cannot be considered complete. - Minor errors: Eq. (4), for consistency, should use the identity matrix for the covariance matrix definition. Eq. (10) uses subscript d where it should be using subscript n Eq. (17) includes p(X^L|F^L) in the definition of Q(...), but it shouldn*t. That was particularly misleading, since if we take eq. (17) to be correct (which I did at first), then p(X^L|F^L) cancels out and should not appear in eq. (20). Eq. (23) uses Q(F|A) to mean the same as P(F|A) as far as I understand. Then why use Q? - References [1] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. [2] Beal, M.J. Variational Algorithms for Approximate Bayesian Inference. [3] M.K. Titsias and N.D. Lawrence. Bayesian Gaussian process latent variable model. [4] M. Lázaro-Gredilla. Bayesian warped Gaussian processes.