On Fast Dropout and its Applicability to Recurrent Networks

Introduction

Recurrent Neural Networks are among the most powerful models for sequential data. The capabilty of representing any measurable sequence to sequence mapping to arbitrary accuracy (Hammer, 2000) makes them universal approximators. Nevertheless they were given only little attention in the last two decades due to the problems of vanishing and exploding gradients (Hochreiter, 1991; Bengio et al., 1994; Pascanu et al., 2012). Error signals either blowing up or decaying exponentially for events many time steps apart rendered them largely impractical for the exact problems they were supposed to solve. This made successful training impossible on many tasks up until recently without resorting to special architectures or abandoning gradient-based optimization. Successful application on tasks with long-range dependencies has thus relied on one of those two paradigms. The former ist to make use of long short-term memory (LSTM) (Hochreiter and Schmidhuber, 1997). These approaches are among the best methods for the modelling of speech and handwriting (Graves et al., 2013; Graves et al., 2008; Graves, 2013). The latter is to to rely on sensible initializations leading to echo-state networks (Jäger et al., 2003).

The publication of (Martens and Sutskever, 2011) can nowadays be considered a landmark, since it was shown that even standard RNNs can be trained with the right optimization method. While a sophisticated Hessian-free optimizer was employed initially, further research (Sutskever et al., 2013; Bengio et al., 2012) has shown that carefully designed first-order methods can find optima of similar quality.

After all, the problem of underfitting standard RNNs can be dealt with to the extent that RNNs are practical in many areas, e.g., language modelling (Sutskever et al., 2011; Mikolov et al., 2010). In contrast, the problem of overfitting in standard RNNs has (due to the lack of necessity) been tackled only by few. As noted in (Pascanu et al., 2012), using priors with a single optima on the parameters may have detrimental effects on the representation capability of RNNs: a global attractor is constructed in parameter space. In the case of a prior with a mode at zero (e.g. an $L2$ regularizer) this biases the network towards solutions which lets information die out exponentially fast in time, making it impossible to memorize events for an indefinite amount of time.

Graves, 2011 proposes to stochastically and adaptively distort the weights of LSTM-based RNNs, which is justified from the perspective of variational Bayes and the minimum description length principle. Overfitting is practically non-existent in the experiments conducted. It is untested whether this approach works well for standard RNNs–along the lines of the observations of Pachitariu and Sahani, 2013 one might hypothesize that the injected noise disturbs the dynamics of RNNs too much and leads to divergence during training.

The deep neural network community has recently embraced a regularization method called dropout (Hinton et al., 2012). The gist is to randomly discard units from the network during training, leading to less interdependent feature detectors in the intermediary layers. Here, dropping out merely means to set the output of that unit to zero. An equivalent view is to set the complete outgoing weight vector to zero of which it is questionable whether a straight transfer of dropout to RNNs is possible. The resulting changes to the dynamics of an RNN during every forward pass are quite dramatic. This is the reason why Pachitariu and Sahani, 2013 only use dropout on those parts of the RNN which are not dynamic., i.e. the connections feeding from the hidden into the output layer.

Our contribution is to show that using a recent smooth approximation to dropout (Wang and Manning, 2013) regularizes RNNs effectively. Since the approximation is deterministic, we may assert that all dynamic parts of the network operate in reasonable regimes. We show that fast dropout does not keep RNNs from reaching rich dynamics during training, which is not obvious due to the relation of classic dropout to L2 regularization (Wager et al., 2013).

The structure of the paper is as follows. We will first review RNNs and fast dropout (FD) (Wang and Manning, 2013). A novel analysis of the derivatives of fast dropout leads to an interpretation where we can perform a decomposition into a loss based on the average output of a network’s units and a regularizer based on its variance. We will discuss why this is a form that is well suited to RNNs and consequently conduct experiments that confirm our hypothesis.

Methods

In this section we will first review RNNs and fast dropout. We will then introduce a novel interpretation of what the fast dropout loss constitutes in section 2.2.2 and show relationships to two other regularizers.

We will define RNNs in terms of two components. For one, we are ultimately interested in an output $\mathbf{y}$, which we can calculate given the parameters $\theta$ of a network and some input $\mathbf{x}$. Secondly, we want to learn the parameters, which is done by the design and optimization of a function of the parameters $\mathcal{L}(\theta)$, commonly dubbed loss, cost, or error function.

where $\{W_{\text{in}},W_{\text{out}},W_{\text{rec}}\}$ are weight matrices and $\{b_{h},b_{y}\}$ bias vectors. These form the set of parameters $\theta$ together with initial hidden state $h_{0}$. The dimensionalities of all weight matrices, bias vectors, and initial hidden states are determined by the dimensionalities of the input sequences as well as desired hidden layer and output layer sizes. The functions $f_{h}$ and $f_{y}$ are so-called transfer functions and mostly coordinate-wise applied nonlinearities. We will call the activations of units pre-synaptic before the application of $f$ and post-synaptic afterwards. Typical choices include the logistic sigmoid $f(\xi)={{1}\over{1+\exp(-\xi)}}$, tangent hyperbolicus and, more recently, the rectified linear $f(\xi)=\max(\xi,0)$ (Zeiler et al., 2013). If we set the recurrent weight matrix $W_{\text{rec}}$ to zero, we recover a standard neural network Bishop, 1995 applied independently to each time step.

2 Fast Dropout

In fast dropout (Wang and Manning, 2013), each unit in the network is assumed to be a random variable. To assure tractability, only the first and second moments of those random variables are kept, which suffices for a very good approximation. Since the pre-synaptic activation of each unit is a weighted sum of its incoming units (of which each is dropped out with a certain probability) we can safely assume Gaussianity for those inputs due to the central limit theorem. As we will see, this is sufficient to find efficient ways to propagate the mean and variance through the non-linearity $f$.

We will now inspect the forward propagation for a layer into a single unit, that is

where $\circ$ denotes the element-wise product and $f$ is a non-linear transfer function as before. Let the input layer $\mathbf{x}$ to the unit be Gaussian distributed with diagonal covariance by assumption: $\mathbf{x}\sim\mathcal{N}(\mu_{\mathbf{x}},\sigma^{2}_{\mathbf{x}})$. Furthermore, we have Bernoulli distributed variables indicating whether an incoming unit is not being dropped out organized in a vector $\mathbf{d}$ with $d_{i}\sim\mathcal{B}(p)$, $p$ being the complementary drop out rate. The weight vector $\mathbf{w}$ is assumed to be constant.

A neural network will in practice consist of many such nodes, with some of them, the output units, directly contributing to the loss function $\mathcal{L}$. Others, the input units, will not stem from calculation but come from the data set. Each component of $\mathbf{x}$ represents an incoming unit, which might be an external input to the network or a hidden unit. In general, $y$ will have a complex distribution depending highly on the nature of $f$. Given that the input to a function $f$ is Gaussian distributed, we obtain the mean and variance of the output as follows:

Forward propagation through the non-linearity $f$ for calculation of the post-synaptic activation can be approximated very well in the case of the logistic sigmoid and the tangent hyperbolicus and done exactly in case of the rectifier (for details, see Wang and Manning, 2013). While the rectifier has been previously reported to be a useful ingredient in RNNs (Bengio et al., 2012) we found that it leads to unstable learning behaviour in preliminary experiments and thus neglected it in this study, solely focusing on the tangent hyperbolicus. Other popular transfer functions, such as the softmax, need to be approximated either via sampling or an unscented transform (Julier and Uhlmann, 1997).

To obtain a Gaussian approximation for $a=(\mathbf{d}\circ\mathbf{x})^{T}\mathbf{w}$, we will use $\hat{a}\sim\mathcal{N}(\text{E}[a],\text{V}[a])$. The mean and variance of $a$ can be obtained as follows. Since $\mathbf{d}$ and $\mathbf{x}$ are independent it follows that

For independent random variables $A$ and $B$, $\text{V}[AB]=\text{V}[A]\text{E}[B]^{2}+\text{E}[A]^{2}\text{V}[B]+\text{V}[A]\text{V}[B]$. If we assume the components of $\mathbf{x}$ to be independent, we can write In contrast to Wang and Manning, 2013, we do not drop the variance not related to dropout. In their case, $\text{V}[a]=p(1-p)\mu_{\mathbf{x}}^{2T}\mathbf{w}^{2}$ was used instead of Equation (2).

Furthermore the independency assumption is necessary such that the Lyapunov condition is satisfied (Lehmann, 1999) for the the central limit theorem to hold, ensuring that $a$ is approximately Gaussian.

Propagating the mean and the variance through $f$ via $f_{\mu}$ and $f_{\sigma}$ suffices for determining the pre-synaptic moments of the outgoing units. At the output $\mathbf{y}$ of the whole model, we will simplify matters and ignore the variance. Some loss functions take the variance into account (e.g., a Gaussian log-likelihood as done in (Bayer et al., 2013)). Sampling can be a viable alternative as well.

The extension of fast dropout to recurrent networks is straightforward from a technical perspective. First, we note that we can concatenate the input vector at time step $t$, $x_{t}$ and the hidden state at the previous layer $h_{t-1}$ into a single vector: $c_{t}=[h_{t-1},x_{t}]$. We obtain a corresponding weight matrix by concatenation of the input to hidden and recurrent weight matrices $W_{\text{in}}$ and $W_{\text{rec}}$: $W_{c}=[W_{\text{rec}},W_{\text{in}}]$. We can thus reduce the computation to the step from above.

2.2 Beyond the Backward Pass: A Regularization Term

Given the forward pass, we used automatic differentiation with Theano (Bergstra et al., 2010) to calculate the gradients. Nevertheless, we will contribute a close inspection of the derivatives. This will prove useful since it makes it possible to interpret fast dropout as an additional regularization term independent of the exact choice of loss function.

Consider a loss $\mathcal{L}(\mathcal{D};\theta)$ which is a function of the data $\mathcal{D}$ and parameters $\theta$ We will frequently omit the explicit dependency on $\mathcal{D}$ and $\theta$ where clear from context.. In machine learning, we wish this loss to be minimal under unseen data $\mathcal{D}$ although we only have access to a training set $\mathcal{D}_{\text{train}}$. A typical approach is to optimize another loss $\mathcal{J}(\mathcal{\mathcal{D}_{\text{train}}};\theta)$ as a proxy in the hope that a good minimum of it will correspond to a good minimum of $\mathcal{L}$ for unseen data. Learning is often done by the optimization of $\mathcal{J}=\mathcal{L}+\mathcal{R}$, where $\mathcal{R}$ is called a regularizer. A common example of a regularizer is to place a prior on the parameters, in which case it is a function of $\theta$ and corresponds to the log-likelihood of the parameters. For weight decay, this is a spherical Gaussian with inverse scale $\lambda$, i.e. $\mathcal{R}_{\text{wd}}(\theta)=\lambda||\theta||_{2}^{2}$. Regularizers can be more sophisticated, e.g. Rifai et al., 2011 determine directions in input space to which a model’s outputs should be invariant. More recently, dropout (i.e. non-fast dropout) for generalized linear models has been intepreted as a semi-supervised regularization term encouraging confident predictions by Wager et al., 2013.

While it seems difficult to bring the objective function $\mathcal{J}_{\text{fd}}(\mathcal{D};\theta)$ of fast dropout into the form of $\mathcal{L}+\mathcal{R}$, it is possible with the derivatives of each node. For this, we perform back-propagation like calculations.

Let $a=(\mathbf{d}\circ\mathbf{x})^{T}\mathbf{w}$ and $y=f(a)$ be the pre- and post-synaptic activations of a component of a layer in the network. First note that ${\partial\mathcal{J}/\partial w_{i}}={\partial\mathcal{J}/\partial a}\cdot{\partial a/\partial w_{i}}$ according to the chain rule. Since $a$ is a random variable, it will be described in one of two forms. In the case of a Gaussian approximation, we will summarize it in terms of its mean and variance; this approach is used if propagation through $f$ is possible in closed form. In the case of sampling, we will have a single instantiation $\hat{a}$ of the random variable, which we can propagate through $f$. An analysis of both cases is as follows.

We find the derivative of $\mathcal{J}$ with respect to one of its incoming weights $w_{i}$ to be

We know that $\text{E}[a]=(\mathbf{x}\circ\mathbf{d})^{T}\mathbf{w}$ and thus $\partial\text{E}[a]/\partial w_{i}=x_{i}d_{i}$. This can be recognized as the standard back-propagation term if we consider the dropout variable $d_{i}$ as fixed. We will thus define

and subsequently refer to it as the local derivative of the training loss. The second term can be analysed similarly. We apply the chain-rule once more which yields

for which any further simplification of $\delta^{a}$ depends on the exact form of $\mathcal{J}$. The remaining two factors can be written down explicitly, i.e.

In alignment with Equation (3) this lets us arrive at

and offers an interpretation of fast dropout as an additive regularization term. An important and limiting aspect of this decomposition is that it only holds locally at $a$.

We note that depending on the sign of the error signal $\delta^{a}$, fast dropout can take on three different behaviours:

$\delta^{a}=0$ The error signal is zero and thus the variance of the unit considered to be optimal for the loss. The fast dropout term vanishes; this is especially true at optima of the overall loss.

$\delta^{a}<0$ The unit should increase its variance. The exact interpretation of this depends on the loss, but in many cases this is related to the expectation of the unit being quite erroneous and leads to an increase of scatter of the output. The fast dropout term encourages a quadratic growth of the weights.

$\delta^{a}>0$ The unit should decrease its variance. As before, this depends on the exact loss function but will mostly be related to the expectation of the unit being quite right which makes a reduction of scatter desirable. The fast dropout term encourages a quadratic shrinkage of the weights.

This behaviour can be illustrated for output units by numerically inspecting the values and gradients of the pre-synaptic moments given a loss. For that we consider a single unit $y=f(a)$ and a loss $d(y,z)$ measuring the divergence of its output to a target value $z$. The pre-synaptic variance $\text{V}[a]$ can enter the loss not at all or in one of two ways, respected by either the loss (see (Bayer et al., 2013)) or the transfer function. Three examples for this are

Squared loss on the mean, i.e. $d(y,z)=(\text{E}[y]-t)^{2}$ with $y=a$,

Gaussian log-likelihood on the moments, i.e. $d(y,z)\propto{(\text{E}[y]-t)^{2}\over 2\text{V}[y]}+\log\sqrt{2\pi\text{V}[y]}$ with $y=a$,

Negative Bernoulli cross entropy, i.e. $d(y,z)=z\log\text{E}[y]+(1-z)\log(1-\text{E}[y])$ with $y={{1}\over{1+\exp(-a)}}$.

We visualize the pre-synaptic mean and variance, their gradients and their respective loss values in Figure 1. For the two latter cases, erroneous units first increase the variance, then move towards the correct mean and subsequently reduce the variance.

An already mentioned alternative to the forward propagation of $\text{E}[a]$ and $\text{V}[a]$ through $f$ is to incarnate $\hat{a}$ via sampling and calculate $f(\hat{a})$. Let $s\sim\mathcal{N}(0,1)$ and $\hat{a}=\text{E}[a]+s\sqrt{\text{V}[a]}$. We can then use $\hat{a}$ explicitly and it follows that

Again, we recognize ${\partial\mathcal{J}/\partial\hat{a}}\cdot{\partial\hat{a}/\partial\text{E}[a]}\cdot{\partial\text{E}[a]/\partial{w_{i}}}={\partial\mathcal{J}/\partial\hat{a}}~x_{i}d_{i}$ as the standard back-propagation formula with dropout variables.

which, making use of results from earlier in the section is equivalent to

The value of this is a zero-centred Gaussian random variable, since $s$ is Gaussian. The scale is independent of the current weight value and only determined by the post-synaptic moments of the incoming unit, the dropout rate and the error signal. We conclude, that also in this case, we can write

where ${\partial\mathcal{R}_{\text{draw}}^{a}/\partial w_{i}}$ is defined as in Equation (4) and essentially an adaptive noise term.

We want to stress the fact that in the approximation as well as in the sampling case the regularization term vanishes at any optima of the training loss. A consequence of this is that no global attractor is formed, which makes the method theoretically useful for RNNs. One might argue that fast dropout should not have a regularizing effect all. Yet, regularization is not only influencing the final solution but also the optimization, leading to different optima.

As already mentioned, imposing a Gaussian distribution centred at the origin with precision $\lambda$ as a prior on the weights leads to a method called weight decay. It is not only probabilistically sound but also works well empirically, see e.g. (Bishop et al., 2006). Recalling that the derivative of weight decay is of the form $\partial\mathcal{R}_{\text{wd}}/\partial w_{i}=2\lambda w_{i}$ we can reinterpret $\mathcal{R}_{\text{approx}}$ as a weight decay term where the coefficient $\lambda$ is weight-wise, dependent on the current activations and possibly negative. Weight decay will always be slightly wrong on the training set, since the derivative of the weight decay term has to match the one of the unregularized loss. In order for $\partial\mathcal{L}/\partial\theta+\partial\mathcal{R}_{\text{wd}}/\partial\theta$ to be minimal, $\mathcal{L}$ cannot be minimal unless so is $\mathcal{R}$.

In this method (Graves, 2011) not the units but the weights are stochastic, which is in pratice implemented by performing Monte Carlo sampling. We can use a similar technique to FD to find a closed form approximation. In this case, a layer is $y=f(\mathbf{x}^{T}\mathbf{w})$, where we have no dropout variables and the weights are Gaussian distributed with $\mathbf{w}\sim\mathcal{N}(\mu_{\mathbf{w}},\sigma^{2}_{\mathbf{w}})$, with covariance diagonal and organized into a vector. We assume Gaussian density for $a=\mathbf{x}^{T}\mathbf{w}$. Using similar algebra as above, we find that

It is desirable to determine whether fast dropout and “fast adaptive weight noise” are special cases of each other. Showing that the aproaches are different can be done by equating equations (1) and (5) and solving for $\mu_{\mathbf{w}}$. This shows that rescaling by $1-p$ suffices in the case of the expectation. It is however not as simple as that for the variance, i.e. for equations (2) and (6), where the solution depends on $\mu_{\mathbf{x}}$ and $\sigma^{2}_{\mathbf{x}}$ and thus is not independent of the input to the network. Yet, both methods share the property that no global attractor is present in the loss: the “prior” is part of the optimization and not fixed.

3 Bag of Tricks

Throughout the experiments we will resort to several “tricks” that have been introduced recently for more stable and efficient optimization of neural networks and RNNs especially. First, we make use of rmsprop (Tieleman and Hinton, 2012), an optimizer which divides the gradient by an exponential moving average of its squares. This approach is similar to Adagrad (Duchi et al., 2011), which uses a window based average. We found that enhancing rmsprop with Nesterov’s accelerated gradient (Sutskever, 2013) greatly reduces the training time in preliminary experiments.

To initialize the RNNs to stable dynamics we followed the initialization protocol of (Sutskever et al., 2013) of setting the spectral radius $\rho$ to a specific value and the maximum amount of incoming connections of a unit to $\nu$; we did not find it necessary to centre the inputs and outputs. The effect of not only using the recurrent weight matrix for propagating the states through time but also its element-wise square for advancing the variances can be quantified. The stability of a network is coupled to the spectral radius of the recurrent weight matrix $\rho(W_{\text{rec}})$; thus, the stability of forward propagating the variance is related to the spectral radius of its element-wise square $\rho(W_{\text{rec}}^{2})$. Since $\rho(AB)\leq\rho(A)\rho(B)$ for non-negative matrices and non-singular matrices $A$ and $B$ (Horn and Johnson, 2012), setting $W_{\text{rec}}$ to full rank and its spectral radius to $\rho(W_{\text{rec}})$ assures $\rho(W_{\text{rec}}^{2})=\rho(|W_{\text{rec}}|^{2})\leq\rho(W_{\text{rec}})^{2}$, where $|\cdot|$ denotes taking the absolute value element-wise. We also use the gradient clipping method introduced in (Pascanu et al., 2012), with a fixed threshold of 225.

Since the hidden-to-hidden connections and the hidden-to-output connections in an RNN can make use of hidden units in quite distinct ways, we found it beneficial to separate the dropout rates. Specifically, a hidden unit may have a different probability to be dropped out when feeding into the hidden layer at the next time step than when feeding into the output layer. Taking this one step further, we also consider networks in which we completely neglect fast dropout for the hidden-to-output connections; an ordinary forward pass is used instead. Note that this is not the same as setting the dropout rate to zero, since the variance of the incoming units is completely neglected. Whether this is done is treated as another hyperparameter for the experiment.

Experiments and Results

All experiments were done by performing a random search (Bergstra and Bengio, 2012) over the hyper parameters (see Table 2 in the Appendix for an overview), where 32 runs were performed for each data set. We report the test loss of the model with the lowest validation error over all training runs, using the same split as in (Bengio et al., 2012). To improve speed, we organize sequences into minibatches by first splitting all sequences of the training and validation set into chunks of length of 100. Zeros are prepended to those sequences which have less than 100 time steps. The test error is reported on the unsplitted sequences.

Training RNNs to generatively model polyphonic music is a valuable benchmark for RNNs due to its high dimensionality and the presense of long as well as short term dependencies. This data set has been evaluated previously by Bengio et al., 2012 where the model achieving the best results, RNN-NADE (Boulanger-Lewandowski et al., 2013), makes specific assumptions about the data (i.e. binary observables). RNNs do not attach any assumptions to the inputs.

The data consists of four distinct data sets, namely Piano-midi.de (classical piano music), Nottingham (folk music), MuseData (orchestral) and JSBChorales (chorales by Johann Sebastian Bach). Each has a dimensionality of 88 per time step organized into different piano rolls which are sequences of binary vectors; each component of these vectors indicates whether a note is occuring at the given time step. We use the RNN’s output to model the sufficient statistics of a Bernoulli random variable, i.e.

which describes the probability that note $i$ is present at time step $t$. The output non-linearity $f_{y}$ of the network is a sigmoid which projects the points to the interval $(0,1)$. We perform learning by the minimization of the average negative log-likelihood (NLL); in this case, this is the average binary cross-entropy

where $k$ indices the training sample, $i$ the component of the target and $t$ the time step.

1.2 Results

Although a common metric for evaluating the performance of such benchmarks is that of accuracy (Bay et al., 2009) we restrict ourselves to that of the NLL–the measure of accuracy is not what is optimized and to which the NLL is merely a proxy. We present the results of FD-RNNs compared with the various other methods in Table 1. Our method is only surpassed by methods which either incorporate more specific assumptions of the data or employ various forms of depth (Boulanger-Lewandowski et al., 2013; Pascanu et al., 2013). We want to stress that we performed only 32 runs for each data set once more. This shows the relative ease to obtain good results despite of the huge space of potential hyper parameters.

One additional observation is the range of the Eigenvalues of the recurrent weight matrix $W_{\text{rec}}$ during training. We performed an additional experiment on JSBChorales where we inspected the Eigenvalues and and the test loss. We found that the spectral radius first increases sharply to a rather high value and then decreases slowly to settle to a specific value. We tried to replicate this behaviour in plain RNNs, but found that RNNs never exceeded a certain spectral radius at which they stuck. This stands in line with the observation from Section 2.2.2 that weights are encouraged to grow when the error is high and shrink during convergence to the optimum. See Figure 2 for a plot of the spectral radius $\rho(W_{\text{rec}})$ over the training process stages.

Conclusion

We have contributed to the field of neural networks in two ways. First, we have analysed a fast approximation of the dropout regularization method by bringing its derivative into the same form as that of a loss regularized with an additive term. We have used this form to gain further insights upon the behaviour of fast dropout for neural networks in general and shown that this objective function does not bias the solutions to those which perform suboptimal on the unreguarized loss. Second, we have hypothesized that this is beneficial especially for RNNs We confirmed this hypothesis by conducting quantitative experiments on an already established benchmark used in the context of learning recurrent networks.

References

Appendix

We show the hyper parameters ranges for the musical data in Table 2. The ones from which the numbers in Table 1 resulted are given in Table 3.