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. 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.
references.authors: Sanjeev Arora Nadav Cohen Elad Hazan
references.title: On the optimization of deep networks: Implicit acceleration by overparameterization
references.venue: Proceedings of the 35th International Conference on Machine Learning,
references.authors: J. Bjorck C. Gomes B. Selman
references.title: Understanding Batch Normalization
references.venue: ArXiv e-prints,
references.authors: Léon Bottou Frank E Curtis Jorge Nocedal
references.title: Optimization methods for large-scale machine learning
references.venue: SIAM Review,
references.authors: Tim Cooijmans Nicolas Ballas César Laurent Aaron C. Courville
references.title: Recurrent batch normalization
references.venue: CoRR, abs/1603.09025,
references.authors: Xavier Glorot Yoshua Bengio
references.title: Understanding the difficulty of training deep feedforward neural networks
references.venue: In Proceedings of the thirteenth international conference on artificial intelligence and statistics,
references.authors: Kaiming He Xiangyu Zhang Shaoqing Ren Jian Sun
references.title: Deep residual learning for image recognition
references.venue: In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR),
references.authors: Sergey Ioffe
references.title: Batch renormalization: Towards reducing minibatch dependence in batch-normalized models
references.venue: CoRR, abs/1702.03275,
references.authors: Sergey Ioffe Christian Szegedy
references.title: Batch normalization: Accelerating deep network training by reducing internal covariate shift
references.venue: In International Conference on Machine Learning,
references.authors: J. Kohler H. Daneshmand A. Lucchi M. Zhou K. Neymeyr T. Hofmann
references.title: Towards a Theoretical Understanding of Batch Normalization
references.authors: Steven G Krantz Harold R Parks
references.title: A primer of real analytic functions
references.venue: Springer Science & Business Media,
references.authors: Alex Krizhevsky Geoffrey Hinton
references.title: Learning multiple layers of features from tiny images
references.venue: Technical report, Citeseer,
references.authors: Yann LeCun Léon Bottou Yoshua Bengio Patrick Haffner
references.title: Gradient-based learning applied to document recognition
references.venue: Proceedings of the IEEE,
references.authors: J.D. Lee M. Simchowitz M.I. Jordan B. Recht
references.title: Gradient Descent Converges to Minimizers
references.venue: ArXiv e-prints,
references.authors: Yintai Ma Diego Klabjan
references.title: Convergence analysis of batch normalization for deep neural nets
references.venue: CoRR, 1705.08011,
references.authors: Ioannis Panageas Georgios Piliouras
references.title: Gradient Descent Only Converges to Minimizers: NonIsolated Critical Points and Invariant Regions
references.venue: 8th Innovations in Theoretical Computer Science Conference (ITCS 2017),
references.authors: S.P. Ponomarev
references.title: Submersions and preimages of sets of measure zero
references.venue: Siberian Mathematical Journal,
references.authors: Yousef Saad
references.title: Iterative methods for sparse linear systems, volume 82
references.venue: siam,
references.authors: S. Santurkar D. Tsipras A. Ilyas A. Madry
references.title: How Does Batch Normalization Help Optimization? (No, It Is Not About Internal Covariate Shift)
references.authors: Michael Shub
references.title: Global stability of dynamical systems
references.venue: Springer Science & Business Media,
references.authors: Han Xiao Kashif Rasul Roland Vollgraf
references.title: Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017
sections.heading: 1 INTRODUCTION
sections.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.
sections.heading: 1.1 RELATED WORK
sections.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).
sections.heading: 1.2 ORGANIZATION
sections.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.
sections.heading: 2 BACKGROUND
sections.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.
sections.heading: 2.1 ORDINARY LEAST SQUARES
sections.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)
sections.heading: 2.2 BATCH NORMALIZATION
sections.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.
sections.heading: 3 MATHEMATICAL ANALYSIS OF BNGD ON OLS
sections.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.
sections.heading: 3.1 SCALING PROPERTY
sections.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).
sections.heading: 3.2 BATCH NORMALIZATION CONVERGES FOR ARBITRARY STEP SIZE
sections.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).
sections.heading: Sketch of Proof.
sections.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.
sections.heading: 3.3 CONVERGENCE RATE, ACCELERATION AND ASYMPTOTIC SENSITIVITY
sections.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.
sections.heading: 4 EXPERIMENTS
sections.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.
sections.heading: 4.1 EXPERIMENTS ON OLS
sections.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.
sections.heading: 4.2 EXPERIMENTS ON PRACTICAL DEEP LEARNING PROBLEMS
sections.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.
sections.heading: 5 CONCLUSION AND OUTLOOK
sections.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.
sections.heading: A PROOF OF THEOREMS
sections.heading: A.1 GRADIENTS AND HESSIAN MATRIX
sections.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.
sections.heading: A.2 SCALING PROPERTY
sections.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.
sections.heading: A.3 PROOF OF THEOREM 3.3
sections.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.
sections.heading: A.5 CONVERGENCE RATE
sections.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.
sections.heading: A.6 ESTIMATING THE EFFECTIVE STEP SIZE
sections.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)
sections.heading: B MODIFIED BNGD
sections.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.
sections.heading: B.1 PROOF
sections.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.
sections.heading: B.2 EXPERIMENT
sections.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).