summaries: The author analyze the convergence properties of batch normalization for the ordinary least square (OLS) objective. They also provide experimental results on the OLS objective as well as small scale neural networks. First of all, understanding the properties of batch normalization is an important topic in the machine learning community so in that sense, contributions that tackle this problem are of interest for the community. However, this paper has a significant number of problems that need to be addressed before publication, perhaps the most important one being the overlap with prior work. Please address this point clearly in your rebuttal.

summaries: This paper provides a theoretical analysis for batch normalization with gradient descent (GDBN) under a simplified scenario, i.e., solving an ordinary least squares problem. The analysis shows that GDBN converges to a stationary point when the learning rate is less than or equal to 1, regardless of the condition number of the problem. Some practical experiments are carried out to justify their theoretical insights. The paper is in general easy to follow.

abstractText: Despite its empirical success, the theoretical underpinnings of the stability, convergence and acceleration properties of batch normalization (BN) remain elusive. In this paper, we attack this problem from a modeling approach, where we perform a thorough theoretical analysis on BN applied to a simplified model: ordinary least squares (OLS). We discover that gradient descent on OLS with BN has interesting properties, including a scaling law, convergence for arbitrary learning rates for the weights, acceleration effects, as well as insensitivity to the choice of learning rates. We then demonstrate numerically that these findings are not specific to the OLS problem and hold qualitatively for more complex supervised learning problems. This points to a new direction towards uncovering the mathematical principles that underlies batch normalization.

authors:
Sanjeev Arora
Nadav Cohen
Elad Hazan

title: On the optimization of deep networks: Implicit acceleration by overparameterization

venue: Proceedings of the 35th International Conference on Machine Learning,

authors:
J. Bjorck
C. Gomes
B. Selman

title: Understanding Batch Normalization

venue: ArXiv e-prints,

authors:
Léon Bottou
Frank E Curtis
Jorge Nocedal

title: Optimization methods for large-scale machine learning

venue: SIAM Review,

authors:
Tim Cooijmans
Nicolas Ballas
César Laurent
Aaron C. Courville

title: Recurrent batch normalization

authors:
Xavier Glorot
Yoshua Bengio

title: Understanding the difficulty of training deep feedforward neural networks

venue: In Proceedings of the thirteenth international conference on artificial intelligence and statistics,

authors:
Kaiming He
Xiangyu Zhang
Shaoqing Ren
Jian Sun

title: Deep residual learning for image recognition

venue: In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR),

title: Batch renormalization: Towards reducing minibatch dependence in batch-normalized models

authors:
Sergey Ioffe
Christian Szegedy

title: Batch normalization: Accelerating deep network training by reducing internal covariate shift

authors:
J. Kohler
H. Daneshmand
A. Lucchi
M. Zhou
K. Neymeyr
T. Hofmann

title: Towards a Theoretical Understanding of Batch Normalization

authors:
Steven G Krantz
Harold R Parks

title: A primer of real analytic functions

venue: Springer Science & Business Media,

authors:
Alex Krizhevsky
Geoffrey Hinton

title: Learning multiple layers of features from tiny images

venue: Technical report, Citeseer,

authors:
Yann LeCun
Léon Bottou
Yoshua Bengio
Patrick Haffner

title: Gradient-based learning applied to document recognition

venue: Proceedings of the IEEE,

authors:
J.D. Lee
M. Simchowitz
M.I. Jordan
B. Recht

title: Gradient Descent Converges to Minimizers

authors:
Yintai Ma
Diego Klabjan

title: Convergence analysis of batch normalization for deep neural nets

authors:
Ioannis Panageas
Georgios Piliouras

title: Gradient Descent Only Converges to Minimizers: NonIsolated Critical Points and Invariant Regions

venue: 8th Innovations in Theoretical Computer Science Conference (ITCS 2017),

title: Submersions and preimages of sets of measure zero

venue: Siberian Mathematical Journal,

title: Iterative methods for sparse linear systems, volume 82

venue: siam,

authors:
S. Santurkar
D. Tsipras
A. Ilyas
A. Madry

title: How Does Batch Normalization Help Optimization? (No, It Is Not About Internal Covariate Shift)

title: Global stability of dynamical systems

authors:
Han Xiao
Kashif Rasul
Roland Vollgraf

title: Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017

heading: 1 INTRODUCTION

text:
Batch normalization (Ioffe & Szegedy, 2015) (BN) is one of the most important techniques for training deep neural networks and has proven extremely effective in avoiding gradient blowups during back-propagation and speeding up convergence. In its original introduction (Ioffe & Szegedy, 2015), the desirable effects of BN are attributed to the so-called “reduction of covariate shift”. However, it is unclear what this statement means in precise mathematical terms. To date, there lacks a comprehensive theoretical analysis of the effect of batch normalization.
In this paper, we study the convergence and stability of gradient descent with batch normalization (BNGD) via a modeling approach. More concretely, we consider a simplified supervised learning problem: ordinary least squares regression, and analyze precisely the effect of BNGD when applied to this problem. Much akin to the mathematical modeling of physical processes, the least-squares problem serves as an idealized “model” of the effect of BN for general supervised learning tasks. A key reason for this choice is that the dynamics of GD without BN (hereafter called GD for simplicity) in least-squares regression is completely understood, thus allowing us to isolate and contrast the additional effects of batch normalization.
The modeling approach proceeds in the following steps. First, we derive precise mathematical results on the convergence and stability of BNGD applied to the least-squares problem. In particular, we show that BNGD converges for any constant learning rate ε ∈ (0, 1], regardless of the conditioning of the regression problem. This is in stark contrast with GD, where the condition number of the problem adversely affect stability and convergence. Many insights can be distilled from the analysis of the OLS model. For instance, we may attribute the stability of BNGD to an interesting scaling law governing ε and the initial condition; This scaling law is not present in GD. The preceding analysis also implies that if we are allowed to use different learning rates for the BN rescaling variables (εa) and the remaining trainable variables (ε), we may conclude that BNGD on our model converges for any ε > 0 as long as εa ∈ (0, 1]. Furthermore, we discover an acceleration effect of BNGD and moreover, there exist regions of ε such that the performance of BNGD is insensitive to changes in ε, which help to explain the robustness of BNGD to the choice of learning rates. We reiterate that contrary to many previous works, all the preceding statements are precise mathematical results that we derive for our simplified model.
The last step in our modeling approach is also the most important: we need to demonstrate that these insights are not specific features of our idealized model. Indeed, they should be true characteristics, at least in an approximate sense, of BNGD for general supervised learning problems. We do this
by numerically investigating the convergence, stability and scaling behaviors of BNGD on various datasets and model architectures. We find that the key insights derived from our idealized analysis indeed correspond to practical scenarios.

heading: 1.1 RELATED WORK

text:
Batch normalization was originally introduced in (Ioffe & Szegedy, 2015) and subsequently studied in further detail in (Ioffe, 2017). Since its introduction, it has become an important practical tool to improve stability and efficiency of training deep neural networks (He et al., 2016; Bottou et al., 2018). Initial heuristic arguments attribute the desirable features of BN to concepts such as “covariate shift”, which lacks a concrete mathematical interpretation and alternative explanations have been given (Santurkar et al., 2018). Recent theoretical studies of BN includes (Ma & Klabjan, 2017), where the authors proposed a variant of BN, the diminishing batch normalization (DBN) algorithm and analyzed the convergence of the DBN algorithm, showing that it converges to a stationary point of the loss function. More recently, (Bjorck et al., 2018) demonstrated that the higher learning rates of batch normalization induce a regularizing effect. Another related work is (Kohler et al., 2018), where the authors also considered the convergence properties of BNGD on linear networks (similar to the least-squares problem), as well as other special problems, such as learning halfspaces and extensions. In the OLS case, the authors showed that for a particularly adaptive choice of dynamic learning rate schedule, which can be seen as a fixed effective step size in our terminology (see equation (11) and the discussion that immediately follows), BNGD converges linearly if λmax is known. Moreover, the analysis also requires setting the rescaling parameter a every step to satisfy a stationarity condition, instead of simply performing gradient descent on a, as is done in the original BNGD.
The present research differs from these previous analysis in an important way - we study the BNGD algorithm itself, and not a special variant. More specifically, we consider constant learning rates (without knowledge of properties of the OLS loss function) and we perform gradient descent on rescaling parameters. We prove that the convergence occurs for even in this case (and in fact, for arbitrarily large learning rates for ε, as long as 0 < εa ≤ 1). This poses more challenges in the analysis and contrasts our work with previous analysis on modified versions of BNGD. This is an important distinction; While a decaying or dynamic learning rate is sometimes used in practice, in the case of BN it is critical to analyze the non-asymptotic, constant learning rate case, precisely because one of the key practical advantages of BN is that a bigger learning rate can be used than that in GD. Hence, it is desirable, as in the results presented in this work, to perform our analysis in this regime.
Finally, through the lens of the least-squares example, BN can be viewed as a type of overparameterization, where additional parameters, which do not increase model expressivity, are introduced to improve algorithm convergence and stability. In this sense, this is related in effect to the recent analysis of the implicit acceleration effects of over-parameterization on gradient descent (Arora et al., 2018).

heading: 1.2 ORGANIZATION

text: Our paper is organized as follows. In Section 2, we outline the ordinary least squares (OLS) problem and present GD and BNGD as alternative means to solve this problem. In Section 3, we demonstrate and analyze the convergence of the BNGD for the OLS model, and in particular contrast the results with the behavior of GD, which is completely known for this model. We also discuss the important insights to BNGD that these results provide us with. We then validate these findings on more general supervised learning problems in Section 4. Finally, we conclude in Section 5.

heading: 2 BACKGROUND

text:
Consider the simple linear regression model where x ∈ Rd is a random input column vector and y is the corresponding output variable. Since batch normalization is applied for each feature separately, in order to gain key insights it is sufficient to the case of y ∈ R. A noisy linear relationship is assumed between the dependent variable y and the independent variables x, i.e. y = xTw + noise
where w ∈ Rd is the parameters. Denote the following moments: H := E[xxT ], g := E[xy], c := E[y2]. (1)
To simplify the analysis, we assume the covariance matrix H of x is positive definite and the mean E[x] of x is zero. The eigenvalues of H are denoted as λi(H), i = 1, 2, ...d,. Particularly, the maximum and minimum eigenvalue of H is denoted by λmax and λmin respectively. The condition number of H is defined as κ := λmaxλmin . Note that the positive definiteness of H allows us to define the vector norms ‖.‖H and ‖.‖H−1 by ‖x‖2H = xTHx and ‖x‖2H−1 = x TH−1x respectively.

heading: 2.1 ORDINARY LEAST SQUARES

text:
The ordinary least squares (OLS) method for estimating the unknown parameters w leads to the following optimization problem
min w∈Rd
J0(w) := 1 2Ex,y[(y − x Tw)2] = c2 − g Tw + 12w THw. (2)
The gradient of J0 with respect to w is ∇wJ0(w) = Hw − g, and the unique minimizer is w = u := H−1g. The gradient descent (GD) method (with step size or learning rate ε) for solving the optimization problem (2) is given by the iterating sequence,
wk+1 = wk − ε∇wJ0(wk) = (I − εH)wk + εg, (3) which converges if 0 < ε < 2λmax =: εmax, and the convergence rate is determined by the spectral radius ρε := ρ(I − εH) = maxi{|1− ελi(H)|} with ‖u− wk+1‖ ≤ ρε‖u− wk‖. (4) It is well known (for example see Chapter 4 of (Saad, 2003)) that the optimal learning rate is εopt = 2 λmax+λmin , where the convergence estimate is related to the condition number κ(H):
‖u− wk+1‖ ≤ κ−1κ+1‖u− wk‖. (5)

heading: 2.2 BATCH NORMALIZATION

text:
Batch normalization is a feature-wise normalization procedure typically applied to the output, which in this case is simply z = xTw. The normalization transform is defined as follows:
NBN (z) := z−E[z]√
Var[z] = x Tw σ , (6)
where σ := √ wTHw. After this rescaling, NBN (z) will be order 1, and hence in order to reintroduce the scale (Ioffe & Szegedy, 2015), we multiply NBN (z) with a rescaling parameter a (Note that the shift parameter can be set zero since E[wTx|w] = 0). Hence, we get the BN version of the OLS problem (2):
min w∈Rd,a∈R
J(a,w) : = 12Ex,y [( y − aNBN (xTw) )2] = c2 − wT g σ a+ 1 2a 2. (7)
The objective function J(a,w) is no longer convex. In fact, it has trivial critical points, {(a∗, w∗)|a∗ = 0, w∗T g = 0}, which are saddle points of J(a,w). We are interested in the nontrivial critical points which satisfy the relations,
a∗ = sign(s) √ uTHu,w∗ = su, for some s ∈ R \ {0}. (8)
It is easy to check that the nontrivial critical points are global minimizers, and the Hessian matrix at each critical point is degenerate. Nevertheless, the saddle points are strict (Details can be found in Appendix), which typically simplifies the analysis of gradient descent on non-convex objectives (Lee et al., 2016; Panageas & Piliouras, 2017).
Consider the gradient descent method to solve the problem (7), which we hereafter call batch normalization gradient descent (BNGD). We set the learning rates for a and w to be εa and ε respectively. These may be different, for reasons which will become clear in the subsequent analysis. We thus have the following discrete-time dynamical system
ak+1 = ak + εa ( wTk g σk − ak ) , (9)
wk+1 = wk + ε ak σk
( g − w T k g
σ2k Hwk
) . (10)
We now begin a concrete mathematical analysis of the above iteration sequence.

heading: 3 MATHEMATICAL ANALYSIS OF BNGD ON OLS

text: In this section, we discuss several mathematical results one can derive concretely for BNGD on the OLS problem (7). First, we establish a simple but useful scaling property, which an important ingredient in allowing us to prove a linear convergence result for arbitrary constant learning rates. We also derive the asymptotic properties of the “effective” learning rate of BNGD (to be precisely defined subsequently), which shows some interesting sensitivity behavior of BNGD on the chosen learning rates. Detailed proofs of all results presented here can be found in the Appendix.

heading: 3.1 SCALING PROPERTY

text:
In this section, we discuss a straightforward, but useful scaling property that the BNGD iterations possess. Note that the dynamical properties of the BNGD iteration are governed by a set of numbers, or a configuration {H,u, a0, w0, εa, ε}. Definition 3.1 (Equivalent configuration). Two configurations, {H,u, a0, w0, εa, ε} and {H ′, u′, a′0, w′0, ε′a, ε′}, are said to be equivalent if for iterates {wk}, {w′k} following these configurations respectively, there is an invertible linear transformation T and a nonzero constant t such that w′k = Twk, a ′ k = tak for all k.
The scaling property ensures that equivalent configurations must converge or diverge together, with the same rate up to a constant multiple. Now, it is easy to check the system has the following scaling law.
Proposition 3.2 (Scaling property). Suppose µ 6= 0, γ 6= 0, r 6= 0, QTQ = I , then
(1) The configurations {µQTHQ, γ√µQu, γa0, γQw0, εa, ε} and {H,u, a0, w0, εa, ε} are equivalent.
(2) The configurations {H,u, a0, w0, εa, ε} and {H,u, a0, rw0, εa, r2ε} are equivalent.
It is worth noting that the scaling property (2) in Proposition 3.2 originates from the batchnormalization procedure and is independent of the specific structure of the loss function. Hence, it is valid for general problems where BN is used (Lemma A.9). Despite being a simple result, the scaling property is important in determining the dynamics of BNGD, and is useful in our subsequent analysis of its convergence and stability properties (see the sketch of the proof of Theorem 3.3).

heading: 3.2 BATCH NORMALIZATION CONVERGES FOR ARBITRARY STEP SIZE

text:
We have the following convergence result.
Theorem 3.3 (Convergence for BNGD). The iteration sequence (ak, wk) in equation (9)-(10) converges to a stationary point for any initial value (a0, w0) and any ε > 0, as long as εa ∈ (0, 1]. Particularly, we have the following sufficient conditions of converging to global minimizers.
(1) If a0wT0 g > 0 (or a0 = 0, w T 0 g 6= 0), εa ∈ (0, 1] and ε is sufficiently small (the smallness
is quantified by Lemma A.13), then (ak, wk) converges to a global minimizer.
(2) If εa = 1 and ε > 0, then (ak, wk) converges to global minimizers for almost all initial values (a0, w0).

heading: Sketch of Proof.

text:
We first prove that the algorithm converges for any εa ∈ (0, 1] and small enough ε, with any initial value (a0, w0) such that ‖w0‖ ≥ 1 (Lemma A.13). Next, we observe that the sequence {‖wk‖} is monotone increasing, and thus either converges to a finite limit or diverges. The scaling property is then used to exclude the divergent case if {‖wk‖} diverges, then at some k the norm ‖wk‖ should be large enough, and by the scaling property, it is equivalent to a case where ‖wk‖=1 and ε is small, which we have proved converges. This shows that ‖wk‖ converges to a finite limit, from which the convergence of wk and the loss function value can be established, after some work. The proof is fully presented in Theorem A.17 and preceding Lemmas.
In addition, using the ’strict saddle point’ arguments in (Lee et al., 2016; Panageas & Piliouras, 2017), we can prove the set of initial value for which (ak, wk) converges to saddle points has Lebesgue measure 0, provided some conditions, such as when εa = 1, ε > 0 (Lemma A.20). It is important to note that BNGD converges for all step size ε > 0 of wk, independent of the spectral properties of H . This is a significant advantage and is in stark contrast with GD, where the step size is limited by λmax(H), and the condition number of H intimately controls the stability and convergence rate. Although we only prove the almost sure convergence to global minimizer for the case of εa = 1, we have not encountered convergence to saddles in the OLS experiments even for εa ∈ (0, 2) with initial values (a0, w0) drawn from typical distributions.

heading: 3.3 CONVERGENCE RATE, ACCELERATION AND ASYMPTOTIC SENSITIVITY

text:
Now, let us consider the convergence rate of BNGD when it converges to a minimizer. Compared with GD, the update coefficient before Hwk in equation (10) changed from ε to a complicated term which we named as the effective step size or learning rate ε̂k
ε̂k := ε ak σk
wTk g
σ2k , (11)
and the recurrence relation in place of u− wk is
u− w T k g
σ2k wk+1 = (I − ε̂kH)
( u− w T k g
σ2k wk
) . (12)
Consider the dynamics of the residual ek := u− (wTk g/σ2k)wk, which equals 0 if and only if wk is a global minimizer. Using the property of H-norm (see section A.1), we observe that the effective learning rate ε̂k determines the convergence rate of ek via
‖ek+1‖H ≤ ∥∥∥u− wTk gσ2k wk+1∥∥∥H ≤ ρ(I − ε̂kH)‖ek‖H , (13)
where ρ(I−ε̂kH) is spectral radius of the matrix I−ε̂kH . The inequality (13) shows that the convergence of ek (and hence the loss function, see Lemma A.23) is linear provided ε̂k ∈ (δ, 2/λmax− δ) for some positive number δ. In fact, if we enforce ε̂k = 1/λmax for each k, which is done in the analysis in Kohler et al. (2018), then one immediately obtains the same linear convergence rate. But this requires knowledge of λmax (problem-dependent) and a modification the BNGD algorithm. We instead focus our analysis on the original BNGD algorithm.
Next, let us discuss below an acceleration effect of BNGD over GD. When (ak, wk) is close to a minimizer, we can approximate the iteration (9)-(10) by a linearized system. The Hessian matrix for BNGD at a minimizer (a∗, w∗) is diag(1, H∗/‖w∗‖2), where the matrix H∗ is
H∗ = H − Huu TH
uTHu . (14)
The matrix H∗ is positive semi-definite (H∗u = 0) and has better spectral properties than H , such as a lower pseudo-condition number κ∗ = λ ∗ max
λ∗min ≤ κ, where λ∗max and λ∗min are the maximal
and minimal nonzero eigenvalues of H∗ respectively. Particularly, κ∗ < κ for almost all u (see section A.1 ). This property leads to acceleration effects of BNGD: When ‖ek‖H is small, the contraction coefficient ρ in (13) can be improved to a lower coefficient. More precisely, we have the following result: Proposition 3.4. For any positive number δ ∈ (0, 1), if (ak, wk) is close to a minimizer, such that λmaxε|ak|
σ2k ‖ek‖H ≤ δ, then we have
‖ek+1‖H ≤ min{ρ ∗(I−ε̂kH∗)+δ
1−δ , ρ(I − ε̂kH)}‖ek‖H , (15) where ρ∗(I − ε̂kH) = max{|1− ε̂kλ∗min|, |1− ε̂kλ∗max|}.
Generally, we have ρ∗(I − ε̂kH∗) ≤ ρ(I − ε̂kH) provided ε̂k > 0, and the optimal rate is ρ∗opt := κ∗−1 κ∗+1 ≤ κ−1 κ+1 =: ρopt, where the inequality is strict for almost all u. Hence, the estimate (15) indicates that the optimal BNGD could have a faster convergence rate than the optimal GD, especially when κ∗ is much smaller than κ.
Finally, we discuss the dependence of the effective learning rate ε̂k (and by extension, the effective convergence rate (13) or (15)) on ε. This is in essence a sensitivity analysis on the performance of BNGD with respect to the choice of learning rate. The explicit dependence of ε̂k on ε is quite complex, but we can nevertheless give the following asymptotic estimates.
Proposition 3.5. Suppose εa ∈ (0, 1], a0wT0 g > 0, and ||g||2 ≥ wT0 g
σ20 gTHw0, then
(1) When ε is small enough, ε 1, the effective step size has a same order with ε, i.e. there are two positive constants, C1, C2, independent on ε and k, such that C1 ≤ ε̂k/ε ≤ C2.
(2) When ε is large enough, ε 1, the effective step size has order O(ε−1), i.e. there are two positive constants, C1, C2, independent on ε and k, such that C1 ≤ ε̂kε ≤ C2.
Observe that for finite k, ε̂k is a differentiable function of ε. Therefore, the above result implies, via the mean value theorem, the existence of some ε0 > 0 such that dε̂k/dε|ε=ε0 = 0. Consequently, there is at least some small interval of the choice of learning rates εwhere the performance of BNGD is insensitive to this choice. In fact, empirically this is one commonly observed advantage of BNGD over GD, where the former typically allows for a variety of (large) learning rates to be used without adversely affecting performance. The same is not true for GD, where the convergence rate depends sensitively on the choice of learning rate. We will see later in Section 4 that although we only have a local insensitivity result above, the interval of this insensitivity is actually quite large in practice.

heading: 4 EXPERIMENTS

text:
Let us first summarize our key findings and insights from the analysis of BNGD on the OLS problem.
1. A scaling law governs BNGD, where certain configurations can be deemed equivalent 2. BNGD converges for any learning rate ε > 0, provided that εa ∈ (0, 1]. In particular,
different learning rates can be used for the BN variables (a) compared with the remaining trainable variables (w)
3. There exists intervals of ε for which the performance of BNGD is not sensitive to the choice of ε
In the subsequent sections, we first validate numerically these claims on the OLS model, and then show that these insights go beyond the simple OLS model we considered in the theoretical framework. In fact, much of the uncovered properties are observed in general applications of BNGD in deep learning.

heading: 4.1 EXPERIMENTS ON OLS

text:
Here we test the convergence and stability of BNGD for the OLS model. Consider a diagonal matrix H = diag(h) where h = (1, ..., κ) is a increasing sequence. The scaling property (Proposition 3.2) allows us to set the initial value w0 having same 2-norm with u, ‖w0‖ = ‖u‖ = 1. Of course, one can verify that the scaling property holds strictly in this case.
Figure 1 gives examples of H with different condition numbers κ. We tested the loss function of BNGD, compared with the optimal GD (i.e. GD with the optimal step size εopt), in a large range of step sizes εa and ε, and with different initial values of a0. Another quantity we observe is the effective step size ε̂k of BN. The results are encoded by four different colors: whether ε̂k is close to the optimal step size εopt, and whether loss of BNGD is less than the optimal GD. The results indicate that the optimal convergence rate of BNGD can be better than GD in some configurations. This acceleration phenomenon is ascribed to the pseudo-condition number of H∗ (discard the only zero eigenvalue) being less than κ(H). This advantage of BNGD is significant when the (pseudo)condition number discrepancy between H and H∗ is large. However, if this difference is small, the acceleration is imperceptible. This is consistent with our analysis in section 3.3.
Another important observation is a region such that ε̂ is close to εopt, in other words, BNGD significantly extends the range of ‘optimal’ step sizes. Consequently, we can choose step sizes in BNGD at greater liberty to obtain almost the same or better convergence rate than the optimal GD. However, the size of this region is inversely dependent on the initial condition a0. Hence, this suggests that small a0 at first steps may improve robustness. On the other hand, small εa will weaken the performance of BN. The phenomenon suggests that improper initialization of the BN parameters weakens the power of BN. This experience is encountered in practice, such as (Cooijmans et al., 2016), where higher initial values of BN parameter are detrimental to the optimization of RNN models.

heading: 4.2 EXPERIMENTS ON PRACTICAL DEEP LEARNING PROBLEMS

text:
We conduct experiments on deep learning applied to standard classification datasets: MNIST (LeCun et al., 1998), Fashion MNIST (Xiao et al., 2017) and CIFAR-10 (Krizhevsky & Hinton, 2009). The goal is to explore if the key findings outlined at the beginning of this section continue to hold for more general settings. For the MNIST and Fashion MNIST dataset, we use two different networks: (1) a one-layer fully connected network (784 × 10) with softmax mean-square loss; (2) a fourlayer convolution network (Conv-MaxPool-Conv-MaxPool-FC-FC) with ReLU activation function and cross-entropy loss. For the CIFAR-10 dataset, we use a five-layer convolution network (ConvMaxPool-Conv-MaxPool-FC-FC-FC). All the trainable parameters are randomly initialized by the Glorot scheme (Glorot & Bengio, 2010) before training. For all three datasets, we use a minibatch size of 100 for computing stochastic gradients. In the BNGD experiments, batch normalization is performed on all layers, the BN parameters are initialized to transform the input to zero mean/unit variance distributions, and a small regularization parameter =1e-3 is added to variance √ σ2 + to avoid division by zero.
Scaling property Theoretically, the scaling property 3.2 holds for any layer using BN. However, it may be slightly biased by the regularization parameter . Here, we test the scaling property in practical settings. Figure 2 gives the loss of network-(2) (2CNN+2FC) at epoch=1 with different learning rate. The norm of all weights and biases are rescaled by a common factor η. We observe that the scaling property remains true for relatively large η. However, when η is small, the norm of weights are small. Therefore, the effect of the -regularization in √ σ2 + becomes significant, causing the curves to be shifted.
Stability for large learning rates We use the loss value at the end of the first epoch to characterize the performance of BNGD and GD methods. Although the training of models have generally not converged at this point, it is enough to extract some relative rate information. Figure 3 shows the loss value of the networks on the three datasets. It is observed that GD and BNGD with identical learning rates for weights and BN parameters exhibit a maximum allowed learning rate, beyond which the iterations becomes unstable. On the other hand, BNGD with separate learning rates exhibits a much larger range of stability over learning rate for non-BN parameters, consistent with our theoretical results in Theorem 3.3.
Insensitivity of performance to learning rates Observe that BN accelerates convergence more significantly for deep networks, whereas for one-layer networks, the best performance of BNGD and
GD are similar. Furthermore, in most cases, the range of optimal learning rates in BNGD is quite large, which is in agreement with the OLS analysis (Proposition 3.5). This phenomenon is potentially crucial for understanding the acceleration of BNGD in deep neural networks. Heuristically, the “optimal” learning rates of GD in distinct layers (depending on some effective notion of “condition number”) may be vastly different. Hence, GD with a shared learning rate across all layers may not achieve the best convergence rates for all layers at the same time. In this case, it is plausible that the acceleration of BNGD is a result of the decreased sensitivity of its convergence rate on the learning rate parameter over a large range of its choice.

heading: 5 CONCLUSION AND OUTLOOK

text:
In this paper, we adopted a modeling approach to investigate the dynamical properties of batch normalization. The OLS problem is chosen as a point of reference, because of its simplicity and the availability of convergence results for gradient descent. Even in such a simple setting, we saw that BNGD exhibits interesting non-trivial behavior, including scaling laws, robust convergence properties, acceleration, as well as the insensitivity of performance to the choice of learning rates. Although these results are derived only for the OLS model, we show via experiments that these are qualitatively valid for general scenarios encountered in deep learning, and points to a concrete way in uncovering the reasons behind the effectiveness of batch normalization.
Interesting future directions include the extension of the results for the OLS model to more general settings of BNGD, where we believe the scaling law (Proposition 3.2) should play a significant role. In addition, we have not touched upon another empirically observed advantage of batch normalization, which is better generalization errors. It will be interesting to see how far the current approach takes us in investigating such probabilistic aspects of BNGD.

heading: A PROOF OF THEOREMS

heading: A.1 GRADIENTS AND HESSIAN MATRIX

text:
The objective function in problem (7) has an equivalent form:
J(a,w) = 12 (u− a σw) TH(u− aσw) = 1 2‖u‖ 2 H − wT g σ a+ 1 2a 2, (16)
where u = H−1g.
The gradients are: ∂J ∂a = − 1 σ (w THu− aσw THw) = − 1σw
T g + a, (17) ∂J ∂w = − a σ (Hu− a σHw) + a σ3 (w THu− aσw THw)Hw = − aσ g + a σ3 (w T g)Hw. (18)
The Hessian matrix is ( ∂2J ∂a2
∂2J ∂a∂w
∂2J ∂w∂a ∂2J ∂w2
) = ( 1 AT21 A21 A22 ) (19)
where
A22 = a σ3 (w
T g) [ H + 1
wT g
( (Hw)gT + g(Hw)T ) − 3σ2 (Hw)(Hw) T ] , (20)
A21 = − 1σ ( g − 1σ2 (w T g)Hw ) . (21)
The objective function J(a,w) has trivial critical points, {(a∗, w∗)|a∗ = 0, w∗T g = 0}. It is obvious that a∗ is the minimizer of J(a,w∗), but (a∗, w∗) is not a local minimizer of J(a,w) unless g = 0, hence (a∗, w∗) are saddle points of J(a,w). The Hessian matrix at those saddle points has at least a negative eigenvalue, i.e. the saddle points are strict. In fact, the eigenvalues at the saddle
point (a∗, w∗) are {
1 2 (1±
√ 1 + 4 ‖g‖ 2
w∗THw∗ ), 0, ..., 0
} which contains d−2 repeated zero, a positive
and a negative eigenvalue.
On the other hand, the nontrivial critical points satisfies the relations,
a∗ = ± √ uTHu,w∗//u, (22)
where the sign of a∗ depends on the direction of u,w∗, i.e. sign(a∗) = sign(uTw∗). It is easy to check that the nontrivial critical points are global minimizers. The Hessian matrix at those minimizers is diag(1, H∗/‖w∗‖2) where the matrix H∗ is
H∗ = H − Huu TH
uTHu (23)
which is positive semi-definite and has a zero eigenvalue corresponding to the eigenvector u, i.e.H∗u = 0. The following lemma, similar to the well known Cauchy interlacing theorem, gives an estimate of eigenvalues of H∗.
Lemma A.1. IfH is positive definite andH∗ is defined asH∗ = H− Huu TH
uTHu , then the eigenvalues
of H and H∗ satisfy the following inequalities:
0 = λ1(H ∗) < λ1(H) ≤ λ2(H∗) ≤ λ2(H) ≤ ... ≤ λd(H∗) ≤ λd(H). (24)
Here λi(H) means the i-th smallest eigenvalue of H .
Proof. (1) According to the definition, we have H∗u = 0, and for any x ∈ Rd,
xTH∗x = xTHx− (x THu)2
uTHu ∈ [0, xTHx], (25)
which implies H∗ is semi-positive definite, and λi(H∗) ≥ λ1(H∗) = 0. Furthermore, we have the following equality:
xTH∗x = min t∈R ‖x− tu‖2H . (26)
(2) We will prove λi(H∗) ≤ λi(H) for all i, 1 ≤ i ≤ d. In fact, using the Min-Max Theorem, we have
λi(H ∗) = min
dimV=i max x∈V xTH∗x ‖x‖2 ≤ mindimV=imaxx∈V xTHx ‖x‖2 = λi(H).
(3) We will prove λi(H∗) ≥ λi−1(H) for all i, 2 ≤ i ≤ d. In fact, using the Max-Min Theorem, we have
λi(H ∗) = max
dimV=n−i+1 min x∈V xTH∗x ‖x‖2 = maxdimV=n−i+1,u⊥V min x∈V min t∈R ‖x−tu‖2H ‖x‖2
≥ max dimV=n−i+1,u⊥V min x∈V min t∈R
‖x−tu‖2H ‖x−tu‖2
= max dimV=n−i+1 min y∈span{V,u}
‖y‖2H ‖y‖2 , y = x− tu
≥ max dimV=n−(i−1)+1 min y∈V
yTHy ‖y‖2 = λi−1(H),
where we have used the fact that x ⊥ u, ‖x− tu‖2 = ‖x‖2 + t2‖u‖2 ≥ ‖x‖2.
There are several corollaries related to the spectral property of H∗. We first give some definitions. Since H∗ is positive semi-definite, we can define the H∗-seminorm.
Definition A.2. The H∗-seminorm of a vector x is defined as ‖x‖H∗ := xTH∗x. ‖x‖H∗ = 0 if and only if x is parallel to u.
Definition A.3. The pseudo-condition number of H∗ is defined as κ∗(H∗) := λd(H ∗)
λ2(H∗) .
Definition A.4. For any real number ε, the pseudo-spectral radius of the matrix I − εH∗ is defined as ρ∗(I − εH∗) := max
2≤i≤d |1− ελi(H∗)|.
The following corollaries are direct consequences of Lemma A.1, hence we omit the proofs.
Corollary A.5. The pseudo-condition number of H∗ is less than or equal to the condition number of H :
κ∗(H∗) := λd(H ∗)
λ2(H∗) ≤ λd(H)λ1(H) =: κ(H), (27)
where the equality holds up if and only if u ⊥ span{v1, vd}, vi is the eigenvector of H corresponding to eigenvalue λi(H). Corollary A.6. For any vector x ∈ Rd and any real number ε, we have ‖(I − εH∗)x‖H∗ ≤ ρ∗(I − εH∗)‖x‖H∗ . Corollary A.7. For any positive number ε > 0, we have
ρ∗(I − εH∗) ≤ ρ(I − εH), (28)
where the inequality is strict if uT vi 6= 0 for i = 1, d.
It is obvious that the inequality in (27) and (28) is strict for almost all u.

heading: A.2 SCALING PROPERTY

text:
The dynamical system defined in equation (9)-(10) is completely determined by a set of configurations {H,u, a0, w0, εa, ε}. It is easy to check the system has the following scaling property: Lemma A.8 (Scaling property). Suppose µ 6= 0, γ 6= 0, r 6= 0, QTQ = I , then
(1) The configurations {µQTHQ, γ√µQu, γa0, γQw0, εa, ε} and {H,u, a0, w0, εa, ε} are equivalent.
(2) The configurations {H,u, a0, w0, εa, ε} and {H,u, a0, rw0, εa, r2ε} are equivalent.
The scaling property is valid for general loss functions provided batch normalization is used. Consider a general problem
min w∈Rd
J0(w) := Ex,y[f(y, x Tw)], (29)
and its BN version
min w∈Rd,a∈R
J(a,w) : = Ex,y [ f ( y, aNBN (x Tw) )] . (30)
Then the gradient descent method gives the following iteration,
ak+1 = ak + εa wTk h̃ σk , (31)
wk+1 = wk + ε ak σk
( h̃− w T k h̃
σ2k Hwk
) , (32)
where h̃ = h(akwk/σk), and h is the gradient of original problem:
h(w) := Ex,y[xf ′ 2(y, x Tw)]. (33)
It is easy to check the general BNGD has the following property: Lemma A.9 (General scaling property). Suppose r 6= 0, then the configurations {w0, ε, ∗} and {rw0, r2ε, ∗} are equivalent. Here the sign * means other parameters.

heading: A.3 PROOF OF THEOREM 3.3

text:
Recall the BNGD iterations
ak+1 = ak + εa ( wTk g σk − ak ) ,
wk+1 = wk + ε ak σk
( g − w T k g
σ2k Hwk
) .
The scaling property simplify our analysis by allowing us to set, for example, ‖u‖ = 1 and ‖w0‖ = 1. In the rest of this section, we only set ‖u‖ = 1. For the step size of a, it is easy to check that ak tends to infinity with εa > 2 and initial value a0 = 1, w0 = u. Hence we only consider 0 < εa < 2, which make the iteration of ak bounded by some constant Ca. Lemma A.10 (Boundedness of ak). If the step size 0 < εa < 2, then the sequence ak is bounded for any ε > 0 and any initial value (a0, w0).
Proof. Define αk := wTk g σk , which is bounded by |αk| ≤ √ uTHu =: C, then
ak+1 = (1− εa)ak + εaαk = (1− εa)k+1a0 + (1− εa)kεaα0 + ...+ (1− εa)εaαk−1 + εaαk.
Since |1− εa| < 1, we have |ak+1| ≤ |a0|+ 2C ∑k i=0 |1− εa|i ≤ |a0|+ 2C 1 1−|1−εa| .
According to the iterations (34), we have
u− w T k g
σ2k wk+1 = ( I − ε akσk wTk g σ2k H )( u− w T k g σ2k wk ) . (34)
Define
ek := u− w T k g
σ2k wk, (35)
qk := u THu− (w
T k g) 2
σ2k = ‖ek‖2H ≥ 0, (36)
ε̂k := ε ak σk
wTk g
σ2k , (37)
and using the property w T g σ2k = argmin t ‖u− tw‖H , and the property of H-norm, we have
qk+1 ≤ ∣∣∣∣∣∣u− wTk gσ2k wk+1∣∣∣∣∣∣2H = ‖(I − ε̂kH)ek‖2H ≤ ρ(I − ε̂kH)2qk. (38)
Therefore we have the following lemma to make sure the iteration converge:
Lemma A.11. Let 0 < εa < 2. If there are two positive numbers ε− and ε̂+, and the effective step size ε̂k satisfies
0 < ε −
‖wk‖2 ≤ ε̂k ≤ ε̂ + < 2λmax (39)
for all k large enough, then the iterations (34) converge to a minimizer.
Proof. Without loss of generality, we assume ε −
‖wk‖2 < 1 λmax and the inequality (39) is satisfied for
all k ≥ 0. We will prove ‖wk‖ converges and the direction of wk converges to the direction of u. (1) Since ‖wk‖ is always increasing, we only need to prove it is bounded. We have,
‖wk+1‖2 = ‖wk‖2 + ε2 a 2 k
σ2k ‖Hek‖2 (40)
= ‖w0‖2 + ε2 k∑ i=0 a2i σ2i ‖Hei‖2 (41)
≤ ‖w0‖2 + ε2λmax k∑ i=0 a2i σ2i qi (42)
≤ ‖w0‖2 + ε2 λmaxC 2 a
λmin k∑ i=0 qi ‖wi‖2 . (43)
The inequality in last lines are based on the fact that ‖Hei‖2 ≤ λmax‖ei‖2H , and |ak| are bounded by a constant Ca. Next, we will prove ∑∞ i=0 qi ‖wi‖2 <∞, which implies ‖wk‖ are bounded.
According to the estimate (38), we have
qk+1 ≤ max i {|1− ε̂+λi|2, |1− ε −λi ‖wk‖2 | 2}qk (44)
≤ max{1− γ+, 1− ε −λmin ‖wk‖2 }qk, (45)
where 1− γ+ = maxi{|1− ε̂+λi|2} ∈ (0, 1). Using the definition of qk, we have
qk − qk+1 ≥ min{γ +‖w0‖2,ε−λmin} ‖wk‖2 qk =: Cqk ‖wk‖2 ≥ 0. (46)
Since qk is bounded in [0, uTHu], summing both side of the inequality, we get the bound of the infinite series ∑ k qk ‖wk‖2 ≤ uTHu C <∞.
(2) Since ‖wk‖ is bounded, we denote ε̂− := ε − ‖w∞‖2 , and define ρ := maxi {|1 − ε̂±λi|} ∈ (0, 1), then the inequality (38) implies qk+1 ≤ ρ2qk. As a consequence, qk tends to zero, which implies the direction of wk converges to the direction of u.
(3) The convergence of ak is a consequence of wk converging.
Since ak is bounded, we assume |ak| < C̃a √ uTHu, C̃a ≥ 1, and define ε0 := 12C̃aκλmax . The following lemma gives the convergence for small step size.
Lemma A.12. If the initial values (a0, w0) satisfies a0wT0 g > 0, and step size satisfies εa ∈ (0, 1], ε/‖w0‖2 < ε0, then the sequence (ak, wk) converges to a global minimizer.
Remark 1: If we set a0 = 0, then we have w1 = w0, a1 = εa wT0 g σ0
, hence a1wT1 g > 0 provided wT0 g 6= 0. Remark 2: For the case of εa ∈ (1, 2), if the initial value satisfies an additional condition 0 < |a0| ≤ εa |w T 0 g| σ0 , then we have (ak, wk) converges to a global minimizer as well.
Proof. Without loss of generality, we only consider the case of a0 > 0, wT0 g > 0, ‖w0‖ ≥ 1.
(1) We will prove ak > 0, wTk g > 0 for all k. Denote yk := w T k g, δ = ‖g‖ 4κ .
On the one hand, if ak > 0, 0 < yk < 2δ, then
yk+1 ≥ yk + ε akσk ‖g‖2 2 ≥ yk. (47)
On the other hand, when ak > 0, yk > 0, ε < ε0, we have
yk+1 ≥ εak‖g‖ 2
σk + yk
( 1− ε ak
σ2k
√ gTHg ) ≥ 12yk, (48)
ak+1 ≥ min{ak, yk/σk}. (49)
As a consequence, we have ak > 0, yk ≥ δy := min{y0, δ} for all k by induction. (2) We will prove the effective step size ε̂k satisfying the condition in Lemma A.11.
Since ak is bounded, ε < ε0, we have
ε̂k := ε ak σk
wTk g
σ2k ≤ εC̃aλmaxλmin‖wk‖2 ≤ εC̃aκ =: ε̂ + < 12λmax , (50)
and
qk+1 ≤ (1− ε̂kλmin)2qk ≤ (1− ε̂kλmin)qk < qk. (51)
which implies w T k+1g
σk+1 ≥ w T k g σk ≥ w T 0 g σ0 . Furthermore, we have ak ≥ min{a0, w T 0 g σ0 }, and there is a
positive constant ε− > 0 such that
ε̂k ≥ ε akλmax‖wk‖2 wTk g σk ≥ ε − ‖wk‖2 . (52)
(3) Employing the Lemma A.11, we conclude that (ak, wk) converges to a global minimizer.
Lemma A.13. If step size satisfies εa ∈ (0, 1], ε/‖w0‖2 < ε0, then the sequence (ak, wk) converges.
Proof. Thanks to Lemma A.12, we only need to consider the case of akwTk g ≤ 0 for all k, and we will prove the iteration converges to a saddle point in this case. Since the case of ak = 0 or wTk g = 0 is trivial, we assume akwTk g < 0 below. More specifically , we will prove |ak+1| < r|ak| for some constant r ∈ (0, 1), which implies convergence to a saddle point. (1) If ak and ak+1 have same sign, hence different sign with wTk g, then we have |ak+1| = |1 − εa‖ak| − εa|wTk g|/σk ≤ |1− εa‖ak|. (2) If ak and ak+1 have different signs, then we have
|wTk g| |akσk| ≤ ε 1 σ2k
( ‖g‖2 − w T k g
σ2k gTHwk
) ≤ 2εκλmax < 1. (53)
Consequently, we get
|ak+1| |ak| = εa |wTk g| |akσk| − (1− εa) ≤ 2εεaκλmax − (1− εa) < εa ≤ 1. (54)
(3) Setting r := max(|1− εa|, 2εεaκλmax − (1− εa)), we finish the proof.
To simplify our proofs for Theorem 3.3, we give two lemmas which are obvious but useful.
Lemma A.14. If positive series fk, hk satisfy fk+1 ≤ rfk + hk, r ∈ (0, 1) and lim k→∞ hk = 0, then
lim k→∞ fk = 0.
Proof. It is obvious, because the series bk defined by bk+1 = rbk + hk, b0 > 0, tends to zeros.
Lemma A.15 (Separation property). For δ0 small enough, the set S := {w|y2q < δ0, ‖w‖ ≥ 1} is composed by two separated parts: S1 and S2, dist(S1, S2) > 0, where in the set S1 one has y2 < δ1, q > δ2, and in S2 one has q < δ2, y2 > δ1 for some δ1 > 0, δ2 > 0. Here y := wT g, q := uTHu− (w THu)2
wTHw = uTHu− y
2
wTHw .
Proof. The proof is based on H being positive. The geometric meaning is illustrated in Figure 4.
Corollary A.16. If lim k→∞ ‖wk+1 −wk‖ = 0, and lim k→∞ (wTk g) 2qk = 0, then either lim k→∞ (wTk g) 2 = 0 or lim k→∞ qk = 0.
Proof. Denote yk := wTk g. According to the separation property (Lemma A.15), we can chose a δ0 > 0 small enough such that the separated parts of the set S := {w|y2q < δ0, ‖w‖ ≥ 1}, S1 and S2, have dist(S1, S2) > 0.
Because y2kqk tends to zero, we have wk belongs to S for k large enough, for instance k > k1. On the other hand, because ‖wk+1 − wk‖ tends to zero, we have ‖wk+1 − wk‖ < dist(S1, S2) for k large enough, for instance k > k2. Then consider k > k3 := max(k1, k2), we have all wk belongs to the same part S1 or S2.
If wk ∈ S1, (qk > δ2), for all k > k3, then we have lim k→∞ (wTk g) 2 = 0.
On the other hand, if wk ∈ S2, (y2k > δ1), for all k > k3, then we have lim k→∞ qk = 0.
Theorem A.17. Let εa ∈ (0, 1] and ε > 0. The sequence (ak, wk) converges for any initial value (a0, w0).
Proof. We will prove ‖wk‖ converges, then prove (ak, wk) converges as well. (1) We will prove that ‖wk‖ is bounded and hence converges. In fact, according to the Lemma A.13, once ‖wk‖2 ≥ ε/ε0 for some k, the rest of the iteration will converge, hence ‖wk‖ is bounded. (2) We will prove lim
k→∞ ‖wk+1 − wk‖ = 0, and lim k→∞ (wTk g) 2qk = 0.
The convergence of ‖wk‖ implies ∑ k a 2 kqk is summable. As a consequence,
lim k→∞ a2kpk = 0, lim k→∞ akek = 0, (55)
and lim k→∞ ‖wk+1 − wk‖ = 0. In fact, we have
‖wk+1 − wk‖2 = ε2 a 2 k σ2 ‖Hek‖ 2 ≤ λmaxε 2 λ2min a2kqk → 0. (56)
Consider the iteration of series |ak − wTk g/σk|,∣∣∣ak+1 − wTk+1gσk+1 ∣∣∣ ≤ ∣∣∣ak+1 − wTk+1gσk ∣∣∣+ ∣∣∣wTk+1gσk − wTk+1gσk+1 ∣∣∣ ≤ (1− εa)
∣∣∣ak − wTk gσk ∣∣∣+ ε |akgTHek|σ2k + |wTk+1g|(σkσk+1) |σk+1 − σk| ≤ (1− εa)
∣∣∣ak − wTk gσk ∣∣∣+ ε‖g‖H‖akek‖Hσ2k + |wTk+1g|(σkσk+1)ελmaxσk ‖akek‖H ≤ (1− εa)
∣∣∣ak − wTk gσk ∣∣∣+ 2C‖akek‖H . (57) The constant C in (57) can be chosen as C = ελmax‖u‖Hλmin‖w0‖2 . Since ‖akek‖H tends to zero, we can use Lemma A.14 to get lim
k→∞ |ak − wTk g/σk| = 0. Combine the equation (55), then we have
lim k→∞
(wTk g) 2pk = 0.
(3) According to the Corollary A.16, we have either lim k→∞ y2k = 0, or lim k→∞ qk = 0. In the former case, the iteration of (ak, wk) converges to a saddle point. However, in the latter case, (ak, wk) converges to a global minimizer. In both cases we have (ak, wk) converges.
To finish the proof of Theorem 3.3, we have to demonstrate the special case of εa = 1 where the set of initial values such that BN iteration converges to saddle points is Lebeguse measure zero. We leave this demonstration in next section where we consider the case of εa ≥ 1.
A.4 IMPOSSIBILITY OF CONVERGING TO STRICT SADDLE POINTS
In this section, we will prove the set of initial values such that BN iteration converges to saddle points is (Lebeguse) measure zero, as long as εa ≥ 1. The tools in our proof is similar to the analysis of gradient descent on non-convex objectives (Lee et al., 2016; Panageas & Piliouras, 2017). In addition, we used the real analytic property of the BN loss function (16).
For brevity, here we denote x := (a,w) and let εa = ε, then the BN iteration can be rewrote as
xn+1 = T (xn) := xn − ε∇J(xn). Lemma A.18. If A ⊂ T (Rd/{0}) is a measure zero set, then the preimage T−1(A) is of measure zero as well.
Proof. Since T is smooth enough, according to Theorem 3 of (Ponomarev, 1987), we only need to prove the Jacobian of T (x) is nonzero for almost all x ∈ Rd. In other words, the set {x : det(I − ε∇2J(x)) = 0} is of measure zero. This is true because the function det(I − ε∇2J(x)) is a real analytic function of x ∈ Rd/{0}. (Details of properties of real analytic functions can be found in (Krantz & Parks, 2002) for instance).
Lemma A.19. Let f : X → R be twice continuously differentiable in an open set X ⊂ Rd and x∗ ∈ X be a stationary point of f . If ε > 0, det(I − ε∇2f(x∗)) 6= 0 and the matrix ∇2f(x∗) has at least a negative eigenvalue, then there exist a neighborhood U of x∗ such that the following set B has measure zero,
B := {x0 ∈ U : xn+1 = xn − ε∇f(xn) ∈ U,∀n ≥ 0}. (58)
Proof. The detailed proof is similar to (Lee et al., 2016; Panageas & Piliouras, 2017).
Define the transform function as F (x) := x − ε∇f(x). Since det(I − ε∇2f(x∗)) 6= 0, accorded to the inverse function theorem, there exist a neighborhood U of x∗ such that T has differentiable inverse. Hence T is a local C1 diffeomorphism, which allow us to use the central-stable manifold theorem (Shub, 2013). The negative eigenvalues of ∇2f(x∗) indicates λmax(I − ε∇2f(x∗)) > 1 and the dimension of the unstable manifold is at least one, which implies the set B is on a lower dimension manifold hence B is of measure zero.
Lemma A.20. If εa = ε ≥ 1, then the set of initial values such that BN iteration converges to saddle points is of Lebeguse measure zero.
Proof. We will prove this argument using Lemma A.18 and Lemma A.19. Denote the saddle points set asW := {(a∗, w∗) : a∗ = 0, w∗T g = 0}. The basic point is that the saddle point x∗ := (a∗, w∗)
of the BN loss function (16) has eigenvalues {
1 2 (1 ±
√ 1 + 4 ‖g‖ 2
w∗THw∗ ), 0, ..., 0
} of the Hessian
matrix.
(1) For each saddle point x∗ := (a∗, w∗) of BN loss function, ε ≥ 1 is enough to allow us to use Lemma A.19. Hence there exist a neighborhood Ux∗ of x∗ such that the following set Bx∗ is of measure zero,
Bx∗ := {x0 ∈ Ux∗ : xn ∈ Ux∗ ,∀n ≥ 1}. (59)
(2) The neighborhoods Ux∗ of all x∗ ∈ W forms a cover of W , hence, accorded to Lindelöf’s open cover lemma, there are countable neighborhoods {Ui : i = 1, 2, ...} cover W , i.e. U := ∪iUi ⊇W . As a consequence, the following set A0 is of measure zero,
A0 := ∪iBi = ∪i{x0 ∈ Ui : xn ∈ Ui,∀n ≥ 1}. (60)
(3) Define Am+1 := T−1(Am) = {x ∈ Rd : T (x) ∈ Am},m ≥ 0. According to Lemma A.18, we have all Am and ∪mAm are of measure zero. (4) Since each initial value x0 such that the iteration converges to a saddle point must be contained in some set Am, we finish the proof.
Combine the results of Lemma A.20, scaling property 3.2 and the convergence theorem A.17, we have the following theorem directly.
Theorem A.21. If εa = 1, ε ≥ 0, then the BN iteration (9)-(10) converges to global minimizers for almost all initial values.

heading: A.5 CONVERGENCE RATE

text:
In the last section, we encountered the following estimate for ek = u− w T k g
σ2k wk
‖ek+1‖H ≤ ρ(I − ε̂kH)‖ek‖H . (61)
We can improve the convergence rate of the above if H∗ has better spectral property. This is the content of Proposition 3.4 and the following lemma is enough to prove it.
Lemma A.22. The following inequality holds, (1− δk)‖ek+1‖H ≤ ( ρ∗(I − ε̂kH∗) + δk ) ‖ek‖H , (62)
where δk := λmaxε|ak|
σ2k ‖ek‖H .
Proof. The case of wTk g = 0 is trivial, hence we assume w T k g 6= 0 in the following proof. Rewrite the iteration on wk as the following equality,
u− w T k g σ2k wk+1 = (I − ε̂kH)ek = (I − ε̂kH∗)ek − ε̂k
( 1− (w T k g) 2
uTHuσ2k
) Hu. (63)
Then we will use the properties of H∗-seminorm to prove our argument.
(1) Estimate the H∗-seminorm on the right hand of equation (63).
‖right‖H∗ ≤ ‖(I − ε̂kH∗)ek‖H∗ + |ε̂k| ( 1− (w T k g) 2
uTHuσ2k
) ‖Hu‖H∗ (64)
≤ ρ∗(I − ε̂kH∗)‖ek‖H∗ + λmax|ε̂k|√ uTHu ‖ek‖2H (65) = ρ∗(I − ε̂kH∗) |w T k g|√
uTHuσk ‖ek‖H + λmaxε|akw T k g|√ uTHuσ3k ‖ek‖2H (66)
= |wTk g|√ uTHuσk
( ρ∗(I − ε̂kH∗) + δk ) ‖ek‖H . (67)
(2) Estimate the H∗-seminorm on the left hand of equation (63). Using the H-norm on the iteration of wk, we have
σk+1 = ‖wk + ε akσkHek‖H ≥ σk − ε λmax|ak| σk ‖ek‖H . (68)
Consequently, we have
‖left‖H∗ = |w T k g|√
uTHuσk
σk+1 σk ‖ek+1‖H ≥ |w T k g|√ uTHuσk (1− δk)‖ek+1‖H . (69)
(3) Combining (1) and (2), we finish the proof.
Now, we turn to the convergence of the loss function which can be rewritten as Jk = 12‖ẽk‖ 2 H with ẽk = u− akσkwk. There is an useful equality between ‖ẽk‖ 2 H and ‖ek‖2H :
‖ẽk‖2H = ‖ek‖2H + ( ak − w T k g σk )2 . (70)
Recalling the inequality (57) and the boundedness of ak, we have a constant C0 such that∣∣∣ak+1 − wTk+1gσk+1 ∣∣∣ ≤ |1− εa|∣∣∣ak − wTk gσk ∣∣∣+ C0‖ek‖H , (71) which indicates that we can use the convergence of ek to estimate the convergence of the loss value Jk. In fact we have the following lemma.
Lemma A.23. If ‖ek‖H ≤ Cρk for some constant C and ρ ∈ (0, 1), εa ∈ (0, 1], then we have
‖ẽk‖2H ≤ C2ρ2k + ( C1(1− εa)k + C2kγk )2 , (72)
where γ = max(ρ, 1− εa), C1 = |a0 − wT0 g/σ0| and C2 = CC0.
Proof. According to the inequality (71), we have
∣∣∣ak − wTk gσk ∣∣∣ ≤ C1(1− εa)k + C2 k−1∑ i=0 (1− εa)iρk−i ≤ C1(1− εa)k + C2kγk. (73)
Put it in the equality (70), then we finish the proof.

heading: A.6 ESTIMATING THE EFFECTIVE STEP SIZE

text:
Accorded to Lemma A.12, the effective step size ε̂k has same order with ε‖wk‖2 provided a0w T 0 g > 0, ε/||w0|| < ε0. In fact, we have
C1ε ‖wk‖2 :=
a0w T 0 g
σ0 ε λmax‖wk‖2 ≤ ε̂k ≤ √ uTHu Caελmin‖wk‖2 =: C2ε ‖wk‖2 . (74)
Hence, to prove the Proposition 3.5, we only need to estimate the norm of wk.
proof of Proposition 3.5. According to the BNGD iteration, we have (see the proof of Lemma A.11)
‖wk+1‖2 ≤ ‖w0‖2 + ε2λmax k∑ i=0 a2i σ2i qi. (75)
(1) When ε‖w0‖2 < ε0 (ε0 is defined in Lemma A.12), the sequence qk satisfies qk+1 ≤ (1 − ε̂kλmin)qk. Hence the norm of wk is bounded by
‖wk‖2 ≤ ‖w0‖2 + εκCa σ0wT0 g ∞∑ i=0 (qi − qi+1) ≤ ‖w0‖2 + Cε, (76)
for some constant C. As a consequence,
C̃1ε := C1ε ‖w0‖2(1+Cε0) ≤ ε̂k ≤ C2ε ‖w0‖2 =: C̃2ε. (77)
(2) When ε is large enough, the increment of the norm ‖wk‖ at the first step is large as well. In fact, we have
‖w1‖2 − ‖w0‖2 = ε2 a 2 0
σ20 ‖He0‖2 = C3ε2. (78)
Since ||g||2 ≥ w T 0 g
σ20 gTHw0, we have a1wT1 g > a1w T 0 g > 0. Choose ε to be larger than some value
ε1 such that ε‖w1‖2 < ε0, then we can use the argument in (1) on (a1, w1). More precisely, there are two constants, C1, C2, such that
C1ε ‖w1‖2 ≤ ε̂k ≤ C2ε ‖w1‖2 . (79)
Plugging the equation (78) into it, we have
C1ε 2 1
‖w0‖2+C3ε21 ≤ C1ε
2
‖w0‖2+C3ε2 ≤ ε̂kε ≤ C2ε
2
‖w0‖2+C3ε2 ≤ C2 C3 . (80)

heading: B MODIFIED BNGD

text:
Through our analysis, we discovered that a modification of the BNGD, which we call MBNGD, becomes much easier to analyze and possesses better convergence properties. Note that the results in the main paper do not depend on the results in this section.
The modification is simply to enforce ak = wTk g σk
at every iteration, which yields the modified iterations:
wk+1 = wk + ε wTk g
σ2k
( g − w T k g
σ2k Hwk
) . (81)
In a sense, one can view the above as a limiting version of BNGD where the BN rescaling variable a is adjusted every step to the optimal value based on the current value of the weightsw. The MBNGD iterations is governed by the variables: H,u, a0, w0, ε, where the scaling properties (Lemma 3.2, omit the parameter εa now) remains. More importantly, we find the iteration will converge to a saddle point if and only if it exactly meets the saddle point at a finite step. More precisely, we have the convergence theorem: Theorem B.1 (Convergence for modified BNGD). The iteration sequence wk in equation (81) converges for any initial value w0 and any step size ε > 0. It converges to a global minimizer almost sure, in the sense that the set of initial values such that wk converges to a saddle point is of Lebesgue measure zero. Furthermore, It converges to a saddle point if and only if wTk g = 0 for some k.
Particularly, if ε < 2‖w0‖ 2
λmaxκ‖u‖2 , w T 0 g 6= 0, then wk converges to a global minimizer.
In the following, we assume ‖u‖ = 1.

heading: B.1 PROOF

text:
Lemma B.2. If wT0 g 6= 0 and ε‖w0‖2 < ε0 := 2 κλmax , then the sequence wk converges to a global minimizer.
Proof. Similar to the proof of Lemma A.12, but here the effective step size is always nonnegative which is defined as
ε̂k := ε (wTk g σ2k )2 ≤ εκ‖wk‖2 ≤ εκ‖w0‖2 =: ε̂+ < 2λmax . (82) The inequality (38) immediately gives qk+1 ≤ qk, which implies (wTk+1g) 2
σ2k+1 ≥ (w
T k g) 2
σ2k ≥ (w
T 0 g) 2 σ20 . As
a consequence, the effective step size has a lower bound
ε̂k ≥ ε (w T 0 g) 2
σ20 1 λmax‖wk‖2 =:
ε−
‖wk‖2 . (83)
Employing the Lemma A.11, we conclude that wk converges to a global minimizer. Lemma B.3. If wTk g 6= 0 for all k, then ∑∞ k=0 |wTk g| = ∑∞ k=0(w T k g) 2 =∞.
Proof. Without loss of generality, we assume ‖w0‖ ≥ 1, denote yk := wTk g and set δ = ‖g‖ 4κ . From the iteration of wk, we have
yk+1 = yk + yk σ2k
( ‖g‖2 − yk
σ2k gTHwk
) . (84)
If 0 < |wTk g| < 2δ, then we have the inequality:
‖g‖2 − yk σ2k gTHwk ≥ ‖g‖2 − κ‖g‖‖wk‖ |yk| ≥ 1 2‖g‖ 2, (85)
then
|yk+1| ≥ ( 1 + λmax‖wk‖2 ‖g‖2 2 ) |yk| > |yk| > 0. (86)
As a consequence, limk→∞ wTk g = 0 is not possible unless w T k g = 0 for some k, which implies the results we want.
Theorem B.4. The iteration sequence wk in equation (81) converges for any initial value w0 and any step size ε > 0. Furthermore, wk will converge to a global minimizer unless wTk g = 0 for some k.
Proof. Obviously, if wTk g = 0 for some k = k0, then wk = wk0 for all k ≥ k0, hence wk converges to wk0 . Without losing generality, we consider w T k g 6= 0 for all k and ‖w0‖ ≥ 1 below.
(1) Firstly, we will prove that ‖wk‖ is bounded and hence converges. In fact, according to the Lemma B.2, once ‖wk‖2 ≥ ε/ε0 for some k, the rest of the iteration will converge, hence ‖wk‖ is bounded. (2) Secondly, we will prove wk converges to a vector parallel to u.
Denote yk := wTk g, zk := wTk g σ2 . The convergence of ‖wk‖ indicates that ∑∞ k=0 z 2 kqk is summable,
and then ∑∞ k=0 y 2 kqk is summable as well. Therefore we have
‖wk+1 − wk‖2 = ε2 (w T k g) 2
σ4 ∣∣∣∣∣∣g − wTk gσ2k Hwk∣∣∣∣∣∣2 ≤ λmaxε2z2kqk, (87) and the above tends to zero, i.e. limk→∞ ‖wk+1 − wk‖ = 0. According to the separation property (Lemma A.15), we can chose a δ0 > 0 small enough such that the separated parts of the set S := {w|y2q < δ0, ‖w‖ ≥ 1}, S1 and S2, have dist(S1, S2) > 0. Because y2kqk tends to zero, we have wk belongs to S for k large enough, for instance k > k1. On the other hand, because ‖wk+1 − wk‖ tends to zero, we have ‖wk+1 − wk‖ < dist(S1, S2) for k large enough, for instance k > k2. Then consider k > k3 := max(k1, k2), we have all wk belongs to the same part S1 or S2.
However, Lemma B.3 says ∑∞ k=0 y 2 k = ∞, hence wk ∈ S1 (qk > δ2) for all k > k3 is not true.
Thereforewk ∈ S2 (y2k > δ1) for all k > k3. Consequently, we can claim that ∑∞ k=0 qk is summable and wk converges to a vector parallel to u.

heading: B.2 EXPERIMENT

text:
Here we test the convergence and stability of MBNGD for OLS model. Consider the diagonal matrix H = diag(h), where h = (1, ..., κ) is an increasing sequence. The scaling property allows us to set the initial value w0 having same 2-norm with u, ‖w0‖ = ‖u‖ = 1. Figure 5 gives an example of a 5-dimensional H with condition number κ = 2000. The GD and MBNGD iteration are executed k = 5000 times where u and w0 are randomly chosen from the unit sphere. The values of effective step size, loss ‖ek‖2H and error ‖ek‖ are plotted. Furthermore, to explore the performance of GD and MBNGD, the mean values over 300 random tests are given. It is worth to note that, the geometric mean (G-mean) is more reliable than the arithmetic mean (Amean), where the geometric mean of x can be defined as exp(E(lnx)). Here the reliability means that the G-mean converges quickly when the number of tests increase, however the A-mean does not converge as quickly. In this example, the optimal convergence rate of MBNGD is observably better than GD. This acceleration phenomenon is ascribed to the pseudo-condition number of κ∗(H∗) being less than κ(H). However, if the difference between (pseudo-)condition number of H and H∗ is small, the acceleration is imperceptible.
Another important observation is that the BN significantly extends the range of ‘optimal’ step size, which is embodied by the effective step size ε̂k having a large constant C in ε̂ = O(Cε−1). This means we can chose step size in BN at a large interval to get almost same or better convergence rate than that of the best choice for GD.
Figure 6 gives an example of 100-dimension H with condition number κ = 2000. Similar results as those in the 5-dimensional case are obtained. However, the best optimal convergence rate of MBNGD here has not noticeably improved compared with GD with the optimal learning rate, which is due to the fact that large d decrease the difference between eigenvalues of H and H∗.
Additional tests indicate that:
(1) larger dimensions leads to larger intervals of ‘optimal’ step size, (Figure 7)
(2) the effect of condition number on the ‘optimal’ interval is small (Figure 8).


summaries: This submission closely builds upon an earlier work (Xie et al., 2017, Gal and Ghahramani, 2016) and proposes a new data noising technique motivated by Bayesian RNNs. Specifically, the key contribution is to extend Gal and Ghahramani (2016) to word embedding noising, while drawing inspiration from trational data smoothing techniques. Some variants are discussed, including those motivated by linear interpolation and Kneser-Ney smoothing, just as in Xie et al., 2017. Empirical evaluation is performed with language modeling experiments, and the proposed methods outperforms comparable baselines. One can imagine such a technique could be useful in many other sequence tasks.

summaries: In this submission, the authors present a variational smoothing interpretation of the data noising approach presented in (Xie et al., 2017). Although the theoretical coverage of the problem gives interesting insights. However, a comparison to related work w.r.t. alternative regularization approaches is missing. Similarly, the perplexity values reported in the experimental results on Penn Treebank are far away from state-of-the-art results published by many competitors on this task, e.g. see the current state-of-the-art results on Penn Treebank by (Yang et al., 2017, https://arxiv.org/pdf/1703.02573.pdf and references therein). It is bad practice to ignore existing work completely like this. The interesting question here would be, inhowfar the presented smoothing/regularization methods are complementary to existing approaches, and if the presented methods do provide improvements on top of these.

abstractText: We present a new theoretical perspective of data noising in recurrent neural network language models (Xie et al., 2017). We show that each variant of data noising is an instance of Bayesian recurrent neural networks with a particular variational distribution (i.e., a mixture of Gaussians whose weights depend on statistics derived from the corpus such as the unigram distribution). We use this insight to propose a more principled method to apply at prediction time and propose natural extensions to data noising under the variational framework. In particular, we propose variational smoothing with tied input and output embedding matrices and an element-wise variational smoothing method. We empirically verify our analysis on two benchmark language modeling datasets and demonstrate performance improvements over existing data noising methods.

authors:
Stanley F. Chen
Joshua Goodman

title: An empirical study of smoothing techniques for language modeling

venue: In Proc. of ACL,

authors:
Kyunghyun Cho
Bart van Merrienboer
Caglar Gulcehre
Dzmitry Bahdanau
Fethi Bougares
Holger Schwenk
Yoshua Bengio

title: Learning phrase representations using rnn encoder–decoder for statistical machine translation

venue: In Proc. of EMNLP,

authors:
Andrew M. Dai
Quoc V. Le

title: Semi-supervised sequence learning

venue: In Proc. of NIPS,

authors:
Yarin Gal
Zoubin Ghahramani

title: Dropout as a bayesian approximation: Representing model uncertainty in deep learning

venue: In Proc. of ICML,

title: A theoretically grounded application of dropout in recurrent neural networks

authors:
Edouard Grave
Armand Joulin
Nicolas Usunier

title: Improving neural language models with a continuous cache

venue: In Proc. of ICLR,

title: Neural networks for machine

authors:
Sepp Hochreiter
Jurgen Schmidhuber

title: Long short-term memory

venue: Neural Computation,

authors:
Arthur E. Hoerl
Robert W. Kennard

title: Ridge regression: Biased estimation for nonorthogonal problems

authors:
Hakan Inan
Khashayar Khosravi
Richard Socher

title: Tying word vectors and word classifiers: A loss framework for language modeling

authors:
Mohit Iyyer
Varun Manjunatha
Jordan Boyd-Graber
Hal Daume III

title: Deep unordered composition rivals syntactic methods for text classification

authors:
Ankit Kumar
Ozan Irsoy
Jonathan Su
James Bradbury
Robert English
Brian Pierce
Peter Ondruska
Ishaan Gulrajani
Richard Socher

title: Ask me anything: Dynamic memory networks for natural language processing

authors:
Mitchell Marcus
Grace Kim
Mary Ann Marcinkiewicz
Robert MacIntyre
Ann Bies
Mark Ferguson
Karen Katz
Britta Schasberger

title: The penn treebank: Annotating predicate argument structure

venue: In Proc. of the Workshop on Human Language Technology,

authors:
Gabor Melis
Charles Blundell
Tomas Kocisky
Karl Moritz Hermann
Chris Dyer
Phil Blunsom

title: Pushing the bounds of dropout

venue: arXiv preprint,

authors:
Gabor Melis
Chris Dyer
Phil Blusom

title: On the state of the art of evaluation in neural language models

authors:
Stephen Merity
Caiming Xiong
James Bradbury
Richard Socher

title: Pointer sentinel mixture models

authors:
Stephen Merity
Nitish Shirish Keskar
Richard Socher

title: Regularizing and optimizing lstm language models

authors:
Ofir Press
Lior Wolf

title: Using the output embedding to improve language models. In Proceedings of the 15th Conference of the European Chapter of the Association for Computational Linguistics: Volume 2, Short Papers, pp. 157–163

venue: Association for Computational Linguistics,

authors:
Stanislau Semeniuta
Aliaksei Severyn
Erhardt Barth

title: Recurrent dropout without memory loss

venue: In Proc. of COLING,

authors:
Nitish Srivastava
Geoffrey Hinton
Alex Krizhevsky
Ilya Sutskever
Ruslan Salakhutdinov

title: Dropout: A simple way to prevent neural networks from overfitting

venue: Journal of Machine Learning Research,

authors:
Ke Tran
Arianna Bisazza
Christof Monz

title: Recurrent memory networks for language modeling

venue: In Proc. of NAACL-HLT,

authors:
Ziang Xie
Sida I. Wang
Jiwei Li
Daniel Levy
Aiming Nie
Dan Jurafsky
Andrew Y. Ng

title: Data noising as smoothing in neural network language models

authors:
Zhilin Yang
Zihang Dai
Ruslan Salakhutdinov
William W Cohen

title: Breaking the softmax bottleneck: A high-rank rnn language model

venue: arXiv preprint arXiv:1711.03953,

authors:
Dani Yogatama
Yishu Miao
Gabor Melis
Wang Ling
Adhiguna Kuncoro
Chris Dyer
Phil Blunsom

title: Memory architectures in recurrent neural network language models

heading: 1 INTRODUCTION

text:
Recurrent neural networks (RNNs) are state of the art models in various language processing tasks. However, their performance heavily depends on proper regularization at training time (Melis et al., 2018b; Merity et al., 2018). The two predominant approaches to regularize RNNs are dropout (randomly zeroing out neurons; Srivastava et al., 2014) and `2 regularization (applying `2 penalty to model parameters; Hoerl & Kennard, 1970). Recently, Xie et al. (2017) proposed data noising to regularize language models. Their method is formulated as a data augmentation method that randomly replaces words with other words drawn from a proposal distribution. For example, we can use the unigram distribution which models the number of word occurrences in the training corpus, or a more sophisticated proposal distribution that takes into account the number of bigram types in the corpus. Data noising has been shown to improve perplexity on small and large corpora, and that this improvement is complementary to other regularization techniques and translates to improvements on downstream models such as machine translation.
Xie et al. (2017) derived connections between data noising and smoothing in classical n-gram language models, which we review in §2. In smoothing (Chen & Goodman, 1996), since empirical counts for unseen sequences are zero, we smooth our estimates by a weighted average of higher order and lower order n-gram models. There are various ways to choose the weights and the lower order models leading to different smoothing techniques, with Kneser-Ney smoothing widely considered to be the most effective. Xie et al. (2017) showed that the pseudocounts of a noised data correspond to a mixture of different n-gram models.
In this paper, we provide a new theoretical foundation for data noising and show that it can be understood as a form of Bayesian recurrent neural network with a particular variational distribution (§4.1 and §4.2). Our derivation relates data noising to dropout and variational dropout (Gal & Ghahramani, 2016b), and naturally leads to a data dependent `2 regularization coefficient. We use this insight to arrive at a more principled way to do prediction with data noising—i.e., by taking the mean of the variational distribution, as opposed to the mode—and propose several extensions under the variational framework in §4.3. Specifically, we show how to use variational smoothing with tied input and output embeddings and propose element-wise smooothing. In §5, we validate our analysis in language modeling experiments on the Penn Treebank (Marcus et al., 1994) and Wikitext-2 (Merity et al., 2017) datasets.

heading: 2 RECURRENT NEURAL NETWORK LANGUAGE MODELS

text:
We consider a language modeling problem where the goal is to predict the next word xt given previously seen context words x<t = {x0, x1, . . . , xt−1}. Let W be parameters of a recurrent neural network, and p(xt | x<t;W) = RNN(x<t;W). Following previous work in language modeling (Melis et al., 2018b; Merity et al., 2018), we use LSTM (Hochreiter & Schmidhuber, 1997) as our RNN function, although other variants such as GRU (Cho et al., 2014) can be used as well.
Given a training corpus D = {x0, . . . , xT }, the likelihood we would like to maximize is:
L(W;D) = T∏ t=1 p(xt | x<t;W).
Directly optimizing the (log) likelihood above often leads to overfitting. We typically augment the objective function with a regularizer (e.g., `2 regularizer where R(W) = λ‖W‖22) or use dropout by randomly zeroing out neurons (Srivastava et al., 2014).
Data noising as smoothing. Xie et al. (2017) proposed a method to regularize recurrent neural network language models by noising the data. For each input word i in x<t (sometimes also the corresponding output word), we replace it with another word sampled from a proposal distribution T with probability γi. They introduced various methods on how to choose γi and T.
For example, if γi = γ for all i ∈ V and T is the unigram distribution of words in the corpus, ET[p(xt | noiseγ,T(x<t);W)] corresponds to a mixture of n-gram models with fixed weights (i.e., linear interpolation smoothing).
Another option is to set γi = γ distinct(i,•)
count(i) , where distinct(i, •) denotes the number of distinct continuations preceded by word i (i.e., the number of bigram types that has i as the first word), and count(i) denotes the number of times i appears in the corpus. For this choice of γi, when T is the unigram distribution, the expectation corresponds to absolute discounting; whereas if Ti =
distinct(•,i)∑ v∈V distinct(•,v)
and we replace both the input and output words, it corresponds to bigram Kneser-Ney smoothing. We summarize their proposed methods in Table 1.
At prediction (test) time, Xie et al. (2017) do not apply any noising and directly predict p(xt | x<t;W). They showed that a combination of smoothing and dropout achieves the best result on their language modeling and machine translation experiments.

heading: 3 BAYESIAN RECURRENT NEURAL NETWORKS

text:
In Bayesian RNNs, we define a prior over our parameters W and consider the following maximization problem:
L(D) = ∫ T∏ t=1 p(xt | x<t;W)p(W)dW.
A common prior is the standard normal distribution p(W) = N(0, I).
For recurrent neural networks, the posterior over W given the data D is intractable. We approximate the posterior with a variational distribution q(W), and minimize the KL divergence between the variational distribution and the true posterior by:
KL(q(W)‖p(W | D)) = KL(q(W)‖p(W))− ∫ log p(D |W)q(W)dW + constant
= KL(q(W)‖p(W))− T∑ t=1 ∫ log p(xt | x<t;W)q(W)dW + constant
The integral is often approximated with Monte Carlo integration with one sample W̃ ∼ q(W):1∫ log p(xt | x<t;W)q(W)dW ≈ log p(xt | x<t;W̃)q(W̃)dW̃. (1)
At test time, for a new sequence y0, . . . , yT , we can either set W to be the mean of q(W) or sample and average the results: p(yt | y<t) = 1S ∑S s=1 p(yt | y<t;W̃s), where S is the number of samples and W̃s ∼ q(W).

heading: 4 VARIATIONAL SMOOTHING

text: We now provide theoretical justifications for data smoothing under the variational framework. In a recurrent neural network language model, there are three types of parameters: an input (word) embedding matrix E, an LSTM parameter matrix R, and an output embedding matrix O that produces logits for the final softmax function. We have W = {E,R,O}.

heading: 4.1 LINEAR INTERPOLATION SMOOTHING.

text:
We first focus on the simplest data noising method—linear interpolation smoothing—and show how to extend it to other data noising methods subsequently.
Word embedding matrix E. Denote the word vector corresponding to word i ∈ V in the input embedding matrix E by ei. We obtain a similar effect to linear interpolation smoothing by using the following mixture of Gaussians variational distribution for ei:
ei ∼ q(ei) = (1− γ)N(ei, σI) + γ ∑ v∈V UvN(ev, σI)
= (1− γ + γUi)N(ei, σI) + γ ∑ v 6=i∈V UvN(ev, σI), (2)
where Uv is the unigram probability of word v and σ is small. In other words, with probability γUv , we replace the embedding ei with another embedding sampled from a normal distribution centered at ev .
Note that noising the input word is equivalent to choosing a different word embedding vector to be used in a standard recurrent neural network. Under the variational framework, we sample a different word embedding matrix for every sequence t at training time, since the integral is approximated with Monte Carlo by sampling from q(ei) (Eq. 1).
The above formulation is related to word embedding dropout (Dai & Le, 2015; Iyyer et al., 2015; Kumar et al., 2016), although in word embedding dropout q(ei) = (1− γ)N(ei, σI) + γN(0, σI).
KL term. For p(ei) = N(0, I), we use Proposition 1 in Gal & Ghahramani (2016a) and approximate the KL divergence between a mixture of Gaussians q(ei) and p(ei) as:
KL(q(ei)‖p(ei)) ≈ λ ∑ v∈V θiv 2 ( e>v ev + tr(σI)− V (1 + log 2π)− log |σI| ) = λ
∑ v∈V θiv 2 ‖ev‖22 + constant,
1In the followings, we use x̃ to denote a sample from a distribution q(x).
where V is the vocabulary size and θiv is the mixture proportion for word i: (1 − γ + γUi) for v = i and γ otherwise. In practice, the KL term directly translates to an `2 regularizer on each word embedding vector, but the regularization coefficient is data dependent. More specifically, the regularization coefficient for word vector ei taking into account contributions from KL(q(ev)‖p(ev)) for v ∈ V is:
λ (V − 1)γ + (1− γ + γUi)
2 = λ 1− (V − 2 + Ui)γ 2
(3)
In other words, the variational formulation of data smoothing results in a regularization coefficient that is a function of corpus statistics (the unigram distribution).
Other parameters. For other parameters R and O, we can use either simple variational distributions such as q(R) ∼ N(R, σI) and q(O) ∼ N(O, σI), which become standard `2 regularizers on these parameters; or incorporate dropout by setting q(ri) ∼ (1− α)N(ri, σI) + αN(0, σI), where α is the dropout probability and ri is the i-th row of R (Gal & Ghahramani, 2016b).
Training. In the original noising formulation (Xie et al., 2017), given a sequence {a, b, a, c}, it is possible to get a noised sequence {d, b, a, c}, since the decision to noise at every timestep is independent of others. In the variational framework, we use Monte Carlo integration W̃ ∼ q(W) for each sequence to compute p(xt | x<t;W̃). We can either sample one embedding matrix E per sequence similar to Gal & Ghahramani (2016b), or sample one embedding matrix per timestep (Melis et al., 2018a). While the first method is computationally more efficient (we use it in our experiments), if we decide to noise the first a to d, we will have a noised sequence where every a is replaced by d: {d, b, d, c}. At training time, we go through sequence x<t multiple times (once per epoch) to get different noised sequences.
Predictions. For predictions, Xie et al. (2017) do not noise any input. In the variational framework, this corresponds to taking the mode of the variational distribution, since 1 − γ is almost always greater than γ. They also reported that they did not observe any improvements by sampling. Using the mode is rather uncommon in Bayesian RNNs. A standard approach is to take the mean, so we can set:
ēi = (1− γ + γUi)ei + ∑ v 6=i∈V γUvev.
Extension to absolute discounting. It is straightforward to extend linear interpolation smoothing to absolute discounting by setting γi = γ distinct(i,•) count(i) in Eq. 2 and Eq. 3 above.

heading: 4.2 KNESER-NEY SMOOTHING

text:
We now consider the variational analog of Kneser-Ney data noising. Instead of smoothing towards lower order n-grams, Kneser-Ney uses models that take into account contextual diversity. For example, for bigram Kneser-Ney smoothing, we replace the unigram distribution with Ki =
distinct(•,i)∑ v∈V distinct(•,v)
, where distinct(•, i) denotes the number of distinct bigrams that end with word i. As a result, even if bigrams such as “San Francisco” and “Los Angeles” appear frequently in the corpus; “Francisco” and “Angeles” will not have high probabilities in K since they often follow “San” and “Los”.
Recall that similar to absolute discounting, Kneser-Ney also uses γi = γ distinct(i,•)
count(i) . However, if we replace an input word xt−1, Xie et al. (2017) proposed to also replace the corresponding output word xt. The intuition behind this is that since the probability of replacing an input word xt−1 is proportional to the number of distinct bigrams that start with xt−1, when we replace the input word, we also need to replace the output word (e.g., if we replace “San” by a word sampled from K, we should also replace “Francisco”).
Word embedding matrix E. To get (bigram) Kneser-Ney smoothing, we use the following mixture of Gaussians variational distribution for ei:
ei ∼ q(ei) = (1− γi)N(ei, σI) + γi ∑ v∈V KvN(ev, σI)
= (1− γi + γiKi)N(ei, σI) + γi ∑ v 6=i∈V KvN(ev, σI). (4)
Output embedding matrix O. For the output embedding matrix, we also use the same variational distribution as the word embedding matrix:
oi ∼ q(oi) = (1− γi)N(oi, σI) + γi ∑ v∈V KvN(ov, σI)
= (1− γi + γiKi)N(oi, σI) + γi ∑ v 6=i∈V KvN(ov, σI). (5)
KL term. Following similar derivations in the previous subsection, it is straightforward to show that the approximated KL term introduces the following regularization term to the overall objective:
λ ∑ v∈V 1− (V − 2 +Kv)γv 2 ‖ev‖22 + constant,
and similarly for the output embedding matrix O.
Training. Recall that we use Monte Carlo integration W̃ ∼ q(W) for each sequence to compute p(xt | x<t;W̃) at training time. Since we need to noise the output word when the input word is replaced, we sample γ̃i at t− 1 (for xt−1 = i) and decide whether we will keep the original (input and output) words or replace them. If we decide to replace, we need to sample two new words from KvN(ev, σI) and KvN(ov, σI) respectively. For other words k that are not the target at the current timestep (i.e., k 6= xt), we can either assume that õk = ok or alternatively sample (which can be expensive since we need to sample additional V − 1 times per timestep).
Predictions For predictions, Xie et al. (2017) also take the mode of the variational distribution (assuming 1− γi is almost always greater than γi). We use the mean instead, similar to what we do in variational linear interpolation smoothing.

heading: 4.3 EXTENSIONS

text:
Our derivations above provide insights into new methods under the variational smoothing framework. We describe three variants in the followings:
Tying input and output embeddings. Inan et al. (2017) and Press & Wolf (2017) showed that tying the input and output embedding matrices (i.e., using the same embedding matrix for both) improves language modeling performance and reduces the number of parameters significantly. We take inspirations from this work and sample both O and E from the same base matrix. As a result, we have fewer parameters to train due to this sharing mechanism, but we still have different samples of input and output embedding matrices per sequence.2 Similar to previous results in language modeling (Inan et al., 2017; Melis et al., 2018b), our experiments demonstrate that tying improves the performance considerably.
Combining smoothing and dropout. We can combine variational smoothing and variational dropout by modifying the variational distribution to incorporate a standard Gaussian component:
ei ∼ q(ei) = (1− α) { (1− γi)N(ei, σI) + γi
∑ v∈V KvN(ev, σI)
} + αN(0, σI),
2Alternatively, we could sample one matrix for both the input and output embeddings. We found that this approach is slightly worse in our preliminary experiments.
where α is the dropout probability. Note that in this formulation, we either drop an entire word by setting its embedding vector to zero (i.e., similar to word embedding dropout and blank noising) or choose an embedding vector from the set of words in the vocabulary.
However, it is more common to apply dropout to each dimension of the input embedding and output embedding matrix. In this formulation, we have:
ei,j ∼ q(ei,j) = (1− αi,j) { (1− γi)N(ei,j , σ) + γi
∑ v∈V KvN(ev,j , σ)
} + αi,jN(0, σ). (6)
At training time, while we sample α̃i,j multiple times (once per word embedding dimension for each word), we only sample γ̃i once per word to ensure that when the element is not noised, we still use the same base embedding vector. We use this variant of smoothing and dropout throughout our experiments for our models.
Element-wise smoothing. The variational formulation above allows us to derive a variant of data smoothing that samples each element of the embedding vector ei independently (and similarly for oi). Consider the variational distribution in Eq. 6. At training time, if we sample both α̃i,j and γ̃i,j multiple times (once per word embedding dimension for each word) we arrive at a new element-wise smoothing method. The main difference between this model and the previous model is that each dimension in the input (and output) embedding vector is sampled independently. As a result, the vector that is used is a combination of elements from various word vectors. Notice that the mean under this new scheme is still the same as sampling per vector, so we do not need to change anything at test time. One major drawback about this model is that it is computationally expensive since we need to sample each element of each embedding vector.

heading: 5 EXPERIMENTS

heading: 5.1 SETUP

text: We evaluate our approaches on two standard language modeling datasets: Penn Treebank (PTB) and Wikitext-2. We use a two-layer LSTM as our base language model. We perform non-episodic training with batch size 64 using RMSprop (Hinton, 2012) as our optimization method. We tune the RMSprop learning rate and `2 regularization hyperparameter λ for all models on a development set by a grid search on {0.002, 0.003, 0.004} and {10−4, 10−3} respectively, and use perplexity on the development set to choose the best model. We also tune γ from {0.1, 0.2, 0.3, 0.4}. We use recurrent dropout (Semeniuta et al., 2016) for R and set it to 0.2, and apply (element-wise) input and output embedding dropouts for E and O and set it to 0.5 when E,O ∈ RV×512 and 0.7 when E,O ∈ RV×1024 based on preliminary experiments. We tie the input and output embedding matrices in all our experiments (i.e., E = O), except for the vanilla LSTM model, where we report results for both tied and untied.3

heading: 5.2 MODELS

text:
We compare the following methods in our experiments:
• Baseline: a vanilla LSTM language model. We evaluate two variants: with tied input and output embeddings and without.
• Data noising (DN): an LSTM language model trained with data noising using linear interpolation smoothing or bigram Kneser-Ney smoothing (Xie et al., 2017).
• Variational smoothing (VS): an LSTM language model with variational smoothing using linear interpolation or Kneser-Ney. For both models, we use the mean of the variational distribution at test time.4
3Our preliminary experiments are consistent with previous work (Inan et al., 2017; Melis et al., 2018b; Merity et al., 2018) that show tying the input and output embedding matrices results in better models with fewer numbers of parameters.
4We discuss results using sampling at test time in §5.4.
• Variational element-wise smoothing: for the smaller PTB dataset, we evaluate an LSTM language model that uses elementwise Kneser-Ney variational smoothing and dropout. We also use the mean at test time.

heading: 5.3 RESULTS

text:
Our results are summarized in Table 2. Consistent with previous work on tying the input and output embedding matrices in language models (Inan et al., 2017; Melis et al., 2018b; Merity et al., 2018), we see a large reduction in perplexity (lower is better) when doing so. While our numbers are generally better, the results are also consistent with Xie et al. (2017) that show linear interpolation data noising is slightly better than vanilla LSTM with dropout, and that Kneser-Ney data noising outperforms these methods for both the medium (512) and large (1024) models.
Variational smoothing improves over data noising in all cases, both for linear interpolation and Kneser-Ney. Recall that the main differences between variational smoothing and data noising are: (1) using the mean at test time, (2) having a data dependent `2 regularization coefficient that comes from the KL term,5 and (3) how each method interacts with the input and output embedding tying mechanism.6 Our results suggest that the choice of the proposal distribution to sample from is less important for variational smoothing. In our experiments, Kneser-Ney outperforms linear interpolation on PTB but linear interpolation is slightly better on Wikitext-2.
Element-wise variational smoothing performs the best for both small and large LSTM models on PTB. We note that this improvement comes at a cost, since this method is computationally expensive.7 It
5The data dependent `2 coefficient penalizes vectors that are sampled more often higher. 6 In the variational framework, if we sample one matrix for the input and output embeddings, it effectively noises the output words even for linear interpolation. If we sample two matrices from the same base matrix, these matrices can be different at training time even if the parameters are tied. As described in §4.3, we use the latter in our experiments.
7Element-wise dropout can be implemented efficiently by sampling a mask (zero or one with some probability) and multiply the entire embedding matrix with this mask. In element-wise smoothing, we need to sample an index for each dimension and reconstruct the embedding matrix for each timestep.
took about one day to train the smaller model as opposed to a couple hours without element-wise smoothing. As a result, we were unable to train this on a much bigger dataset such as Wikitext-2. Nonetheless, our results show that applying smoothing in the embedding (latent) space results in better models.

heading: 5.4 DISCUSSIONS

text:
Sampling at test time. In order to better understand another prediction method for these models, we perform experiments where we sample at test time for both data noising and variational smoothing (instead of taking the mode or the mean). We use twenty samples and average the log likelihood. Our results agree with Xie et al. (2017) that mentioned sampling does not provide additional benefits. In our experiments, the perplexities of 512- and 1024-dimensional DN-Kneser-Ney models increase to 96.5 and 85.1 (from 75.0 and 71.4) on the PTB validation set. For VS-Kneser-Ney models, the perplexities increase to 89.8 and 78.7 (from 74.5 and 70.8). Both our results and Xie et al. (2017) suggest that introducing data dependent noise at test time is detrimental for recurrent neural network language models.
Sensitivity to γ. We evaluate the sensitivity of variational smoothing to hyperparameter γ. Figure 1 shows perplexities on the PTB validation set for a variant of our models. The model generally performs well within the range of 0.1 and 0.3, but it becomes progressively worse as we increase γ since there is too much noise.
Other applications. While we focus on language modeling in this paper, the proposed technique is applicable to other language processing tasks. For example, Xie et al. (2017) showed that data smoothing improves machine translation, so our techniques can be used in that setup as well. It is also interesting to consider how variational smoothing interacts with and/or can be applied to memory augmented language models (Tran et al., 2016; Grave et al., 2017; Yogatama et al., 2018) and state-of-the-art language models (Yang et al., 2017). We leave these for future work.

heading: 6 CONCLUSION

text: We showed that data noising in recurrent neural network language models can be understood as a Bayesian recurrent neural network with a variational distribution that consists of a mixture of Gaussians whose mixture weights are a function of the proposal distribution used in data noising (e.g., the unigram distribution, the Kneser-Ney continuation distribution). We proposed using the mean of the variational distribution at prediction time as a better alternative to using the mode. We combined it with variational dropout, presented two extensions (i.e., variational smoothing with tied input and output embedding matrices and element-wise smoothing), and demonstrated language modeling improvements on Penn Treebank and Wikitext-2.

title: VARIATIONAL SMOOTHING IN RECURRENT NEURAL NETWORK LANGUAGE MODELS