Exponential convergence rates for Batch Normalization: The power of length-direction decoupling in non-convex optimization

INTRODUCTION

One of the most important recent innovations for optimizing deep neural networks is Batch Normalization (Bn) (Ioffe and Szegedy,, 2015). This technique has been proven to successfully stabilize and accelerate training of deep neural networks and is thus by now standard in many state-of-the art architectures such as ResNets (He et al.,, 2016) and the latest Inception Nets (Szegedy et al.,, 2017). The success of Batch Normalization has promoted its key idea that normalizing the inner layers of a neural network stabilizes training which recently led to the development of many such normalization methods such as (Arpit et al.,, 2016; Klambauer et al.,, 2017; Salimans and Kingma,, 2016) and (Ba et al.,, 2016) to name just a few.

Yet, despite the ever more important role of Batch Normalization for training deep neural networks, the Machine Learning community is mostly relying on empirical evidence and thus lacking a thorough theoretical understanding that can explain such success. Indeed – to the best of our knowledge – there exists no theoretical result which provably shows faster convergence rates for this technique on any problem instance. So far, there only exists competing hypotheses that we briefly summarize below.

The most widespread idea is that Batch Normalization accelerates training by reducing the so-called internal covariate shift, defined as the change in the distribution of layer inputs while the conditional distribution of outputs is unchanged. This change can be significant especially for deep neural networks where the successive composition of layers drives the activation distribution away from the initial input distribution. Ioffe and Szegedy, (2015) argue that Batch Normalization reduces the internal covariate shift by employing a normalization technique that enforces the input distribution of each activation layer to be whitened - i.e. enforced to have zero means and unit variances - and decorrelated . Yet, as pointed out by Lipton and Steinhardt, (2018), the covariate shift phenomenon itself is not rigorously shown to be the reason behind the performance of Batch Normalization. Furthermore, a recent empirical study published by (Santurkar et al.,, 2018) provides strong evidence supporting the hypothesis that the performance gain of Batch Normilization is not explained by the reduction of internal covariate shift.

Recently, Santurkar et al., (2018) argue that under certain assumptions a normalization layer simplifies optimization by smoothing the loss landscape of the optimization problem of the preceding layer. Yet, we note that this effect may - at best - only improve the constant factor of the convergence rate of Gradient Descent and not the rate itself (e.g. from sub-linear to linear). Furthermore, the analysis treats only the largest eigenvalue and thus one direction in the landscape (at any given point) and keeps the (usually trainable) BN parameters fixed to zero-mean and unit variance. For a thorough conclusion about the overall landscape, a look at the entire eigenspectrum (including negative and zero eigenvalues) would be needed. Yet, this is particularly hard to do as soon as one allows for learnable mean and variance parameters since the effect of their interplay on the distribution of eigenvalues is highly non-trivial.

Finally, a different perspective was brought up by another normalization technique termed Weight Normalization (Wn) (Salimans and Kingma,, 2016). This technique performs a very simple normalization that is independent of any data statistics with the goal of decoupling the length of the weight vector from its direction. The optimization of the training objective is then performed by training the two parts separately. As discussed in Section 2, Bn and Wn differ in how the weights are normalized but share the above mentioned decoupling effect. Interestingly, weight normalization has been shown empirically to benefit from similar acceleration properties as Batch Normalization (Gitman and Ginsburg,, 2017; Salimans and Kingma,, 2016). This raises the obvious question whether the empirical success of training with Batch Normalization can (at least partially) be attributed to its length-direction decoupling aspect.

2 Contribution and organization

We contribute to a better theoretical understanding of Batch Normalization by analyzing it from an optimization perspective. In this regard, we particularly address the following question:

Can we find a setting in which Batch Normalization provably accelerates optimization with Gradient Descent and does the length-direction decoupling play a role in this phenomenon?

We answer both questions affirmatively. In particular, we show that the specific variance transformation of Bn decouples the length and directional components of the weight vectors in such a way that allows local search methods to exploit certain global properties of the optimization landscape (present in the directional component of the optimal weight vector). Using this fact and endowing the optimization method with an adaptive stepsize scheme, we obtain an exponential (or as more commonly termed linear) convergence rate for Batch Norm Gradient Descent on the (possibly) non-convex problem of Learning Halfspaces with Gaussian inputs (Section 4), which is a prominent problem in machine learning Erdogdu et al., (2016). We thereby turn Bn from an effective practical heuristic into a provably converging algorithm. Additionally we show that the length-direction decoupling can be considered as a non-linear reparametrization of the weight space, which may be beneficial for even simple convex optimization tasks such as logistic regressions. Interestingly, non-linear weightspace transformations have received little to no attention within the optimization community (see (Mikhalevich et al.,, 1988) for an exception).

Finally, in Section 5 we analyze the effect of Bn for training a multilayer neural network (MLP) and prove – again under a similar Gaussianity assumption – that Bn acts in such a way that the cross dependencies between layers are reduced and thus the curvature structure of the network is simplified. Again, this is due to a certain global property in the directional part of the optimization landscape, which BN can exploit via the length-direction decoupling. As a result, gradient-based optimization in reparametrized coordinates (and with an adaptive stepsize policy) can enjoy a linear convergence rate on each individual unit. We substantiate both findings with experimental results on real world datasets that confirm the validity of our analysis outside the setting of our theoretical assumptions that cannot be certified to always hold in practice.

BACKGROUND

and make the following (weak) assumption.

[Weak assumption on data distribution] We assume that ${\mathbf{E}}\left[{\bf x}\right]=0$. We further assume that the spectrum of the matrix ${\bf S}$ is bounded as

As a result, ${\bf S}$ is the symmetric positive definite covariance matrix of ${\bf x}$.

The part of our analysis presented in Section 4 and 5 relies on a stronger assumption on the data distribution. In this regard we consider the combined random variable

whose mean vector and covariance matrix are ${\bf u}$ and ${\bf S}$ as defined above in Eq. (1).

[Normality assumption on data distribution] We assume that ${\bf z}$ is a multivariate normal random variable distributed with mean ${\mathbf{E}}\left[{\bf z}\right]={\mathbf{E}}\left[-y{\bf x}\right]={\bf u}$ and second-moment ${\mathbf{E}}\left[{\bf z}{\bf z}^{\top}\right]-{\mathbf{E}}[{\bf z}]{\mathbf{E}}[{\bf z}]^{\top}={\mathbf{E}}\left[{\bf x}{\bf x}^{\top}\right]-{\bf u}{\bf u}^{\top}={\bf S}-{\bf u}{\bf u}^{\top}$.

In the absence of further knowledge, assuming Gaussian data is plausible from an information-theoretic point of view since the Gaussian distribution maximizes the entropy over the set of all absolutely continuous distributions with fixed first and second moment (Dowson and Wragg,, 1973). Thus, many recent studies on neural networks make this assumption on ${\bf x}$ (see e.g. (Brutzkus and Globerson,, 2017; Du and Lee,, 2018)). Here we assume Gaussianity on $y{\bf x}$ instead which is even less restrictive in some casesFor example, suppose that conditional distribution $P({\bf x}|y=1)$ is gaussian with mean $\mu$ for positive labels and $-\mu$ for negative labels (mixture of gaussians). If the covariance matrix of these marginal distributions are the same, $z=y{\bf x}$ is Gaussian while ${\bf x}$ is not..

2 Batch normalization as a reparameterization of the weight space

is normalized on the pre-activation level. That is, the input-output mapping of this unit becomes

As stated (in finite-sum terms) in Algorithm 1 of (Ioffe and Szegedy,, 2015) the normalization operation amounts to computing

and replacing this expression into the batch normalized output of Eq.(5) yields

In order to keep concise notations, we will often use the induced norm of the positive definite matrix ${\bf S}$ defined as $\|{\bf w}\|_{\bf S}:=\left({\bf w}^{\top}{\bf S}{\bf w}\right)^{1/2}$. Comparing Eq. (4) and (8) it becomes apparent that Bn can be considered as a reparameterization of the weight space. We thus define

Note that Weight Normalization (Wn) is another instance of the above reparametrization, where the covariance matrix ${\bf S}$ is replaced by the identity matrix $\mathbf{I}$ (Salimans and Kingma,, 2016). In both cases, the objective becomes invariant to linear scaling of ${\bf w}$. From a geometry perspective, the directional part of Wn can be understood as performing optimization on the unit sphere while Bn operates on the ${\bf S}$-sphere (ellipsoid) (Cho and Lee,, 2017). Note that one can compute the variance term (7) in a matrix-free manner, i.e. ${\bf S}$ never needs to be computed explicitly for Bn.

Of course, this type of reparametrization is not exclusive to applications in neural networks. In the following two sections we first show how reparametrizing the weight space of linear models can be advantageous from a classical optimization point of view. In Section 5 we extend this analysis to training Batch Normalized neural networks with adaptive-stepsize Gradient Descent and show that the length-direction split induces an interesting decoupling effect of the individual network layers which simplifies the curvature structure.

ORDINARY LEAST SQUARES

As a preparation for subsequent analyses, we start with the simple convex quadratic objective encountered when minimizing an ordinary least squares problem

Despite the non-convexity of this new objective, we will prove that Gradient Descent (Gd) enjoys a linear convergence rate. Interestingly, our analysis establishes a link between $f_{\text{{\sc ols}}}$ in reparametrized coordinates (Eq. (12)) and the well-studied task of minimizing (generalized) Rayleigh quotients as it is commonly encountered in eigenvalue problems (Argentati et al.,, 2017).

To simplify the analysis, note that, for a given ${\bf w}$, the objective of Eq. (12) is convex w.r.t. the scalar $g$ and thus the optimal value $g_{\bf w}^{*}$ can be found by setting $\frac{\partial f_{\text{{\sc ols}}}}{\partial g}=0$, which gives $g_{\bf w}^{*}:=-\left({\bf u}^{\top}{\bf w}\right)/\|{\bf w}\|_{\bf S}$. Replacing this closed-form solution into Eq. (12) yields the following optimization problem

which – as discussed in Appendix A.2 – is a special case of minimizing the generalized Rayleigh quotient for which an extensive literature exists (Knyazev,, 1998; D’yakonov and McCormick,, 1995). Here, we particularly consider solving (13) with Gd, which applies the following iterative updates to the parameters

Based upon existing results, the next theorem establishes a linear convergence rate for the above iterates to the minimizer ${\bf w}^{*}$ in the normalized coordinates.

where $\Delta\rho_{t}:=\rho({\bf w}_{t})-\rho({\bf w}^{*})$. Furthermore, the ${\bf S}^{-1}$-norm of the gradient $\nabla\rho({\bf w}_{t})$ relates to the suboptimality as

This convergence rate is of the same order as the rate of standard Gd on the original objective $f_{\text{{\sc ols}}}$ of Eq. (LABEL:eq:least_squares_objective) (Nesterov,, 2013). Yet, it is interesting to see that the non-convexity of the normalized objective does not slow gradient-based optimization down. In the following, we will repeatedly invoke this result to analyze more complex objectives for which Gd only achieves a sublinear convergence rate in the original coordinate space but is provably accelerated after using Batch Normalization.

LEARNING HALFSPACES

We now turn our attention to the problem of Learning Halfspaces, which encompasses training the simplest possible neural network: the Perceptron. This optimization problem can be written as

Furthermore, we need $f_{\text{{\sc lh}}}$ to be sufficiently smooth.

Since globally optimizing (17) is in general NP-hard (Guruswami and Raghavendra,, 2009), we instead focus on understanding the effect of the normalized parameterization when searching for a stationary point. Towards this end we now assume that ${\bf z}$ is a multivariate normal random variable (see Assumption 2 and discussion there).

The learning halfspaces objective $f_{\text{{\sc lh}}}$ – on Gaussian inputs – has a remarkable property: all critical points lie on the same line, independent of the choice of the loss function $\varphi$. We formalize this claim in the next lemma.

Interestingly, the optimal direction of these critical points spans the same line as the solution of a corresponding least squares regression problem (see Eq. (38) in Appendix A). In the context of convex optimization of generalized linear models, this fact was first pointed out in (Brillinger,, 2012). Although the global optima of the two objectives are aligned, classical optimization methods - which perform updates based on local information - are generally blind to such global properties of the objective function. This is unfortunate since Gradient Descent converges linearly in the quadratic least-squares setting but only sublinearly on general Learning Halfspace problems (Zhang et al.,, 2015).

To accelerate the convergence of Gradient Descent, Erdogdu et al., (2016) thus proposed a two-step global optimization procedure for solving generalized linear models, which first involves finding the optimal direction by optimizing a least squares regression as a surrogate objective and secondly searching for a proper scaling factor of that minimizer. Here, we show that running Gd in coordinates reparameterized as in Eq. (9) makes this two-step procedure redundant. More specifically, splitting the optimization problem into searching for the optimal direction and scaling separately, allows even local optimization methods to exploit the property of global minima alignment. Thus - without having to solve a least squares problem in the first place - the directional updates on the Learning Halfspace problem can mimic the least squares dynamics and thereby inherit the linear convergence rate. Combined with a fast (one dimensional) search for the optimal scaling in each step the overall convergence stays linear.

As an illustration, Figure 1 shows the level sets as well as the optimal direction of a least squares-, a logistic- and a sigmoidal regression problem on the same Gaussian dataset. Furthermore, it shows iterates of Gd in original coordinates and a sequential version of Gd in normalized coordinates that first optimizes the direction and then the scaling of its parameters (Gdnpseq). Both methods start at the same point and run with an infinitesimally small stepsize. It can be seen that, while Gd takes completely different paths towards the optimal points of each problem instance, the dynamics of Gdnpseq are exactly the same until the optimal directionDepicted by the dotted red line. Note that – as a result of Lemma 1 – this line is identical in all problems. is found and differ only in the final scaling.

2 Local optimization in normalized parameterization

3 Convergence result

We now show that Algorithm 1 can achieve a linear convergence rate to a critical point on the possibly non-convex objective $f_{\text{{\sc lh}}}$ with Gaussian inputs. Note that all information for computing the adaptive stepsize $s_{t}$ is readily available and can be computed efficiently.

Finally, we note that the proof of this result relies specifically on the ${\bf S}$-reparametrization done by Batch Normalization. In Appendix B.3.6 we detail out why our proof strategy is not suitable for the $\mathbf{I}$-reparametrization of Weight Normalization and thus leave it as an interesting open question if other settings (or proof strategies) can be found where linear rates for Wn are provable.

4 Experiments I

In order to substantiate the above analysis we compare the convergence behavior of Gd and Agd to three versions of Gradient Descent in normalized coordinates. Namely, we benchmark (i) Gdnp (Algorithm 1) with multiple gradient steps on $g$ instead of Bisection, (ii) a simpler version (Bn) which updates ${\bf w}$ and ${\bf g}$ with just one fixed step-size gradient stepThus Bn is conceptually very close to the classical Batch Norm Gradient Descent presented in (Ioffe and Szegedy,, 2015) and (iii) Weight Normalization (Wn) as presented in (Salimans and Kingma,, 2016). All methods use full batch sizes and – except for Gdnp on ${\bf w}$ – each method is run with a problem specific, constant stepsize.

We consider empirical risk minimization as a surrogate for $f_{\text{{\sc lh}}}$ (17) on the common real-world dataset a9a as well as on synthetic data drawn from a multivariate Gaussian distribution. We center the datasets and use two different functions $\varphi(\cdot)$. First, we choose the softplus which resembles the classical logistic regression (convex). Secondly, we use the sigmoid which is a commonly used (non-convex) continuous approximation of the 0-1 loss (Zhang et al.,, 2015). Further details can be found in Appendix D.

The Gaussian design experiments clearly confirm Theorem 2 in the sense that the loss in the convex-, as well as the gradient norm in the non-convex case decrease at a linear rate. The results on a9a show that Gdnp can accelerate optimization even when the normality assumption does not hold and in a setting where no covariate shift is present, which motivates future research of normalization techniques in optimization. Interestingly, the performance of simple Bn and Wn is similar to that of Gd, which suggests that the length-direction decoupling on its own does not capture the entire potential of these methods. Gdnp on the other hand takes full advantage of the parameter splitting, both in terms of multiple steps on $g$ and – more importantly – adaptive stepsizes in ${\bf w}$.

NEURAL NETWORKS

See Appendix C.2 for a discussion of possible implications for deep neural networks.

2 Convergence result

Optimizing each unit independently formally results in minimizing the function $f^{(i)}({\bf w}^{(i)},g^{(i)})$ as defined in Algorithm 2. In the next theorem, we prove that this version of Gdnp achieves a linear rate of convergence to optimize each $f^{(i)}$.

[Convergence of Gdnp on MLP] Suppose Assumptions 1– 4 hold. We consider optimizing the weights $({\bf w}^{(i)},g^{(i)})$ of unit $i$, assuming that all directions $\{{\bf w}^{(j)}\}_{j<i}$ are critical points of $f_{\text{{\sc nn}}}$ and ${\bf w}^{k}={\bf 0}$ for $k>i$. Then, Gdnp with step-size policy $s^{(i)}$ as in (98) and stopping criterion $h^{(i)}$ as in (99) yields a linear convergence rate on $f^{(i)}$ in the sense that

where the constant $C>0$ is defined in Eq. (102).

The result of Theorem 3 relies on the fact that each ${\bf w}^{(j)},j\neq i$ is either zero or has zero gradient. If we assume that an exact critical point is reached after optimizing each individual unit, then the result directly implies that the alternating minimization presented in Algorithm 2 reaches a critical point of the overall objective. Since the established convergence rate for each individual unit is linear, this assumption sounds realistic. We leave a more precise convergence analysis, that takes into account that optimizing each individual unit for a finite number of steps may yield numerical suboptimalities, for future work.

3 Experiments II

In the proof of Theorem 3 we show that Gdnp can leverage the length-direction decoupling in a way that lowers cross-dependencies between hidden layers and yields faster convergence. A central part of the proof is Lemma 2 which says that – given Gaussian inputs – the optimal direction of a given layer is independent of all downstream layers. Since this assumption is rather strong and since Algorithm 2 is intended for analysis purposes only, we test the validity of the above hypothesis outside the Gaussian setting by training a Batch Normalized multilayer feedforward network (Bn) on a real-world image classification task with plain Gradient Descent. For comparison, a second unnormalized network is trained by Gd. To validate Lemma 2 we measure the interdependency between the central and all other hidden layers in terms of the Frobenius norm of their second partial cross derivatives (in the directional component). Further details can be found in Appendix D.

Figure 3 confirms that the directional gradients of the central layer are affected far more by the upstream than by the downstream layers to a surprisingly large extent. Interestingly, this holds even before reaching a critical point. The downstream cross-dependencies are generally decaying for the Batch Normalized network (Bn) (especially in the first 1000 iterations where most progress is made) while they remain elevated in the un-normalized network (Gd), which suggest that using Batch Normalization layers indeed simplifies the networks curvature structure in ${\bf w}$ such that the length-direction decoupling allows Gradient Descent to exploit simpler trajectories in these normalized coordinates for faster convergence.

Of course, we cannot untangle this effect fully from other possible positive aspects of training with Bn (see introduction). Yet, the fact that the (de-)coupling increases in the distance to the middle layer (note how earlier (later) layers are more (less) important for ${\bf W}_{4}$) emphasizes the relevance of this analysis particularly for deep neural network structures, where downstream dependencies might vanish completely with depth. This does not only make gradient based training easier but also suggests the possibility of using partial second order information, such as diagonal Hessian approximations (e.g. proposed in (Martens et al.,, 2012)).

CONCLUSION

We took a theoretical approach to study the acceleration provided by Batch Normalization. In a somewhat simplified setting, we have shown that the reparametrization performed by Batch Normalization leads to a provable acceleration of gradient-based optimization by splitting it into subtasks that are easier to solve. In order to evaluate the impact of the assumptions required for our analysis, we also performed experiments on real-world datasets that agree with the results of the theoretical analysis to a surprisingly large extent.

We consider this work as a first step for two particular directions of future research. First, it raises the question of how to optimally train Batch Normalized neural networks. Particularly, our results suggest that different and adaptive stepsize schemes for the two parameters - length and direction - can lead to significant accelerations. Second, the analysis of Section 3 and 4 reveals that a better understanding of non-linear coordinate transformations is a promising direction for the continuous optimization community.

References

Appendix

Appendix A LEAST SQUARES ANALYSIS

Recall that, after normalizing according to (9) and using the closed form solution for the optimal scaling factor $g^{*}:=-\left({\bf u}^{\top}{\bf w}\right)/\|{\bf w}\|_{\bf S}$, optimizing the ordinary least squares objective can be written as the following minimization problem

We consider optimizing the above objective by Gd which takes iterative steps of the form

The proof of this result crucially relies on the insight that the minimization problem given in (12) resembles the problem of maximizing the generalized Rayleigh quotient which is commonly encountered in generalized eigenproblems. We will thus first review this area, where convergence rates are usually provided in terms of the angle of the current iterate with the maximizer, which is the principal eigenvector. Interestingly, this angle can be related to both, the current function value as well as the the norm of the current gradient. We will make use of these connections to prove the above Theorem in Section A.5. Although not necessarily needed for convex function, we introduce the gradient norm relation as we will later go on to prove a similar result for possibly non-convex functions in the learning halfspace setting (Theorem 2).

A.2 Background on eigenvalue problems

and $q({\bf w})$ is called the Rayleigh quotient. Notably, this quotient satisfies the so-called Rayleigh inequality

where $\lambda_{\min}({\bf B})$ and $\lambda_{1}({\bf B})$ are the smallest and largest eigenvalue of ${\bf B}$ respectively.

Maximizing $q({\bf w})$ is a non-convex (strict-saddle) optimization problem, where the $i$-th critical point ${\bf v}_{i}$ constitutes the $i$-th eigenvector with corresponding eigenvalue $\lambda_{i}=q({\bf w}_{i})$ (see (Absil et al.,, 2009), Section 4.6.2 for details). It is known that optimizing $q({\bf w})$ with Gd - using an iteration-dependent stepsize - converges linearly to the principal eigenvector ${\bf v}_{1}$. The convergence analysis is based on the ”minidimensional” method and yields the following result

under weak assumptions on ${\bf w}_{0}$. Details as well as the proof of this result can be found in (Knyazev and Shorokhodov,, 1991).

Generalized rayleigh quotient

The reparametrized least squares objective (13), however, is not exactly equivalent to (27) because of the covariance matrix that appears in the denominator. As a matter of fact, our objective is a special instance of the generalized Rayleigh quotient

Maximizing (29) is a generalized eigenproblem in the sense that it solves the task of finding eigenvalues $\lambda$ of the matrix pencil $({\bf B},{\bf A})$ for which $\det({\bf B}-\lambda{\bf A})=0$, i.e. finding a vector ${\bf v}$ that obeys ${\bf B}{\bf v}=\lambda{\bf A}{\bf v}$. Again we have

Among the rich literature on solving generalized symmetric eigenproblems, a Gd convergence rate similar to (28) has been established in Theorem 6 of (Knyazev and Neymeyr,, 2003), which yields

again under weak assumptions on ${\bf w}_{0}$.

Our contribution

More importantly, we here have a special case where the nominator of $\rho({\bf w})$ has a particular low rank structure. In fact, ${\bf B}:={\bf u}{\bf u}^{\top}$ is a rank one matrix. Instead of directly invoking the convergence rate in (Knyazev and Neymeyr,, 2003), this allows for a much simpler analysis of the convergence rate of Gd on $\rho({\bf w})$ since the rank one property yields a simpler representation of the relevant vectors. Furthermore, we establish a connection between suboptimality on function value and the ${\bf S}^{-1}$-norm of the gradient. As mentioned earlier, we need such a guarantee in our future analysis on learning halfspaces which is an instance of a (possibly) non-convex optimization problem.

A.3 Preliminaries

Notations Let ${\bf A}$ be a symmetric positive definite matrix. We introduce the following compact notations that will be used throughout the analysis.

These notations allow us to make the analysis similar to the simple Rayleigh quotient case. For example, the denominator in (29) can now be written as $\|{\bf w}\|_{\bf A}^{2}$.

Properties We will use the following elementary properties of the induced terms defined above.

$\sin^{2}\angle_{{\bf A}}({\bf w},{\bf v})=1-\cos^{2}\angle_{{\bf A}}({\bf w},{\bf v})$

If $\hat{{\bf w}}$ is the ${\bf A}$-orthogonal projection of ${\bf w}$ to span$\{{\bf v}\}$, then it holds that

The ${\bf A}$-spectral norm of a matrix can be written in the alternative form

A.4 Characterization of the LS minimizer

By setting the gradient of (LABEL:eq:least_squares_objective) to zero and recalling the convexity of $f_{\text{{\sc ols}}}$ we immediately see that the minimizer of this objective is

where $\lambda_{1}:=\|{\bf u}\|_{{\bf S}^{-1}}^{2}$ is the corresponding generalized eigenvalue. The associated eigenvector with $\lambda_{1}$ is

holds for all $i,j$. Let ${\bf V}_{2}:=[{\bf v}_{2},{\bf v}_{3},\dots,{\bf v}_{d}]$ be the matrix whose $(i-1)$-th column is ${\bf v}_{i},i\in\{2,\dots,d\}$. The matrix ${\bf B}$ is orthogonal to the matrix ${\bf V}_{2}$ since

and since $\alpha_{1}{\bf v}_{1}$ is the ${\bf S}$-orthogonal projection of ${\bf w}$ to span$\{{\bf v}_{1}\}$, the result of (P.2) implies

Clearly this metric is zero for the optimal solution ${\bf v}_{1}$ and else bounded by one from above. To justify it is a proper choice, the next proposition proves that suboptimality on $\rho$, i.e. $\rho({\bf w})-\rho({\bf v}_{1})$, relates directly to this angle.

The suboptimality of ${\bf w}$ on $\rho({\bf w})$ relates to $\sin^{2}\angle_{{\bf S}}({\bf w},{\bf v}_{1})$ as

where $\rho({\bf v}_{1})=\lambda_{1}$. This is equivalent to

We use the proposed eigenexpansion of Eq. (43) to rewrite

and replace the above result into $\rho({\bf w})$. Then

which proves the second part of the proposition. The first follows directly from property (P.1). ∎

Gradient-suboptimality connection Fermat’s first-order optimality condition implies that the gradient is zero at the minimizer of $\rho({\bf w})$. Considering the structure of $\rho({\bf w})$, we propose a precise connection between the norm of gradient and suboptimality. Our analysis relies on the representation of the gradient $\nabla\rho({\bf w})$ in the ${\bf S}$-orthonormal basis $\{{\bf v}_{1},\dots,{\bf v}_{d}\}$ which is described in the next proposition.

Using the ${\bf S}$-orthogonal basis as given in Eq. (41), the gradient vector can be expanded as

The above derivation is based on two results: (i) ${\bf v}_{1}$ is an eigenvector of $({\bf B},{\bf S})$ and (ii) the representation of $\rho({\bf w})$ in Proposition 1. We recall the definition of $\nabla\rho({\bf w})$ in (25) and write

Exploiting the gradient representation of the last proposition, the next proposition establishes the connection between suboptimality and the ${\bf S}^{-1}$-norm of gradient $\nabla_{\bf w}\rho({\bf w})$.

Suppose that $\rho({\bf w})\neq 0$, then the ${\bf S}^{-1}$-norm of the gradient $\nabla\rho({\bf w})$ relates to the suboptimality as

Multiplying the gradient representation in Proposition 2 by ${\bf S}^{-1}$ yields

By combining the above result with the ${\bf S}$-orthogonality of the basis $({\bf v}_{1},{\bf V}_{2})$, we derive the (squared) ${\bf S}^{-1}$-norm of the gradient as

It remains to simplify the terms $T_{1}$ and $T_{2}$. For $T_{1}$,

Replacing the simplified expression of $T_{1}$ and $T_{2}$ into Eq. (LABEL:eq:T1_T2_expansion) yields

A rearrangement of terms in the above equation concludes the proof. ∎

A.5 Convergence proof

We have seen: suboptimality in $\rho({\bf w})$ directly relates to $\sin^{2}\angle_{{\bf S}}({\bf w},{\bf v}_{1})$ for all ${\bf w}\setminus\{\mathbf{0}\}$. In the next lemma we prove that this quantity is strictly decreased by repeated Gd updates at a linear rate.

To prove the above statement, we relate the sine of the angle of a given iterate ${\bf w}_{t+1}$ with ${\bf v}_{1}$ in terms of the previous angle $\angle({\bf w}_{t},{\bf v}_{1})$. Towards this end, we assume for the moment that $\rho({\bf w}_{t})\not=0$ but note that this naturally always holds whenever $\rho({\bf w}_{0})\not=0$,as we will show later by induction. such that the angle relation can be recursively applied through all $t\geq 0$ to yield Eq. (52).

(i) We start by deriving an expression for $\sin\angle_{\bf S}({\bf v}_{1},{\bf w}_{t})$. By (31) and the definition $\rho({\bf w}_{t})$, we have that $-\rho({\bf w}_{t}){\bf w}_{t}$ is the ${\bf S}$-orthogonal projection of ${\bf S}^{-1}{\bf B}{\bf w}_{t}$ to $\text{span}\{{\bf w}_{t}\}$. Indeed,

Note that ${\bf S}^{-1}{\bf B}{\bf w}_{t}=\left({\bf S}^{-1}{\bf u}\right)\left({\bf u}^{\top}{\bf w}_{t}\right)$ is a nonzero multiple of ${\bf v}_{1}$ and thus $\sin\angle_{\bf S}({\bf S}^{-1}{\bf B}{\bf w}_{t},{\bf w}_{t})=\sin\angle_{\bf S}({\bf v}_{1},{\bf w}_{t})$.

By the Pythagorean theorem and (54), we get

(iii) To see how the two quantities on the right hand side of (53) and (56) relate, let us rewrite the Gd updates from Eq. (24) as follows

By taking the ${\bf S}$-norm we can conclude

where the first inequality is due to property (P.4) of the ${\bf S}$-spectral norm and the second is due to Assumption (1) and (P.3) , which allows us to bound the latter in term of the usual spectral norm as follows

(iv) Combining the above results yields the desired bound

Proof of Theorem (1): By combining the results of Lemma 3 as well as Proposition 1 and 3, we can complete the proof of the Theorem 1 as follows

Appendix B LEARNING HALFSPACES ANALYSIS

In this section, we provide a convergence analysis for Algorithm 1 on the problem of learning halfspaces

This result, which can be derived using a simple application of integration by parts, is called Stein’s lemma (Landsman and Nevslehová,, 2008). In the next lemma, we show that this allows us to simplify the expression of the gradient of Eq. 17.

Under the normality assumption on the data distribution (Assumption 2), the gradient of $f_{LH}$ (Eq. 17) can be expressed as

The gradient of $f_{LH}$ can be written as follows

A straight forward application of Stein’s lemma (Eq. (61)) yields

which –after rearrangement – proves the result. See detailed derivation in (Erdogdu et al.,, 2016). ∎

In addition to the assumption on the data distribution, the proposed analysis also requires a rather weak assumption on $f_{\text{{\sc lh}}}$ and loss function $\varphi$. See 3 See 4 Recall that $\zeta$-smoothness of $f_{\text{{\sc lh}}}$, which is mentioned in the last assumption, implies that the gradient of $f_{\text{{\sc lh}}}$ is $\zeta$-Lipschitz, i.e.

B.2 Global characterization

Here, we prove a result about a global property of the solution of the problem of learning halfspaces. See 1

Setting the gradient of the objective $f_{\text{{\sc lh}}}$ as given in Eq. (62) to zero directly gives the result. ∎

B.3 Established Convergence Rate

Based on this assumption, we derive a linear convergence rate for Gdnp presented in Algorithm 1. We first restate the convergence guarantee before providing a detailed proof.

As mentioned earlier, the objective $f_{\text{{\sc lh}}}$ on Gaussian inputs has a particular global property. Namely, all its critical points are aligned along the same direction. The key idea is that ${\bf S}$-reparameterization provides this global information to a local optimization method through an elegant length-direction decoupling. This allows Gdnp to mimic the behaviour of Gradient Descent on the above mentioned Rayleigh quotient for the directional updates and thereby inherit the linear convergence rate. At the same time, the scaling factor can easily be brought to a critical point by a fast, one dimensional search algorithm. We formalize and combine these intuitions in a detailed proof below.

Gradient in the normalized parameterization

Since Gdnp relies on the normalized parameterization, we first need to derive the gradient of the objective in this parameterization

Note that the vector ${\bf S}{\bf w}$ is orthogonal to the column space of ${\bf A}_{\bf w}$ since

We will repeatedly use the above property in our future analysis. In the next lemma, we establish a connection between the norm of gradients in different parameterizations.

Under the reparameterization (9), the following holds:

We introduce the vector ${\bf q}_{1}={\bf S}{\bf w}/\|{\bf w}\|_{\bf S}$ that has unit ${\bf S}^{-1}$-norm, i.e. $\|{\bf q}_{1}\|_{S^{-1}}=1$.

Let ${\bf Q}_{2}$ be a matrix whose columns are $\{{\bf q}_{2},\dots,{\bf q}_{d}\}$. The choice of ${\bf q}_{1}$ together with ${\bf S}^{-1}$-orthogonality of the basis imply that ${\bf w}$ is orthogonal to ${\bf Q}_{2}$:

Consider the gradient expansion in the new basis, i.e.

Plugging the above expansion into Eq. (67) yields

hence the ${\bf S}^{-1}$-norm of the directional gradient in the new parameterization is

Therefore, one can establish the following connection between the ${\bf S}^{-1}$-norm of gradient in the two different parameterizations:

For the particular case of learning halfspaces with Gaussian input, the result of Lemma 4 allows us to write the gradient $\nabla f_{\text{{\sc lh}}}$ as

where the constants $c_{1}$ and $c_{2}$ are determined by the choice of the loss. Replacing this expression in Eq. (67) yields the following formulation for the gradient in normalized coordinates

Convergence of the scalar g𝑔g

According to Algorithm 1, the length of the search space for $g$ is cut in half by each bisection step and thus reduces to

The next step is to relate the above distance to the partial derivative of $f_{\text{{\sc lh}}}({\bf w},g)$ w.r.t $g$. Consider the compact notation ${\bf w}^{\prime}_{t}={\bf w}_{t}/\|{\bf w}_{t}\|_{\bf S}$. Using this notation and the gradient expression in Eq. (67), the difference of partial derivatives can be written as

Using the smoothness assumption on $f_{\text{{\sc lh}}}$ we bound the above difference as follows

where the last inequality is due to Assumption 1.

Directional convergence

The norm of the gradient w.r.t. ${\bf w}$ of $f_{\text{{\sc lh}}}$ as in Eq. (66) converges at the following linear rate

The key insight for this proof is a rather subtle connection between the gradient of the reparametrized least squares objective (Eq. (13)) and the directional gradient of the learning halfspace problem (Eq. (66)):

The stepsize choice of Eq.(80) guarantees that

which exactly matches the Gd iterate sequence of Eq. (14) on $\rho({\bf w})$. At this point, we can invoke the result of Theorem 1 to establish the following convergence rate:

Combined convergence guarantee

Using Proposition 4 and combining the results obtained for optimizing the directional and scalar components, we finally obtain the following convergence guarantee:

A word on Weight Normalization

The improved convergence rate for Batch Normalization (Theorem 2) relies heavily on the fact that normalizing and backpropagating through the variance term resembles splitting the optimization task into a length- and directional component. As mentioned in the introduction, this feature is also present in Weight Normalization and it is thus an obvious question, whether Wn can achieve a similar convergence rate. From a theoretical perspective, we were not able to prove this which is essentially due to the subtle difference in how the normalization is done: While Bn normalizes the parameters to live on the ${\bf S}$-sphere, Wn brings all parameters to the unit sphere.

That said, we observe similar empirical convergence behaviour in terms of suboptimality for Bn and Wn (without any adaptive stepsizes, see Section 4.4) but as can be seen on the right of Figure 7 the path that the two methods take can be very different. We thus leave it as an interesting open question if other settings and proof strategies can be found where fast rates for Wn are provable.

Appendix C NEURAL NETWORKS

Recall the training objective of the one layer MLP presented in Section 5:

Figure 4 illustrates the considered architecture in this paper.

Since the activation function is assumed to be an odd function ($\tanh$), this choice allows us to equivalently rewrite the training objective as

where the scalars $\alpha^{(i)}$, $\beta^{(i)}$ and $\gamma^{(i,j)}$ are defined as

Recall the gradient of $f_{\text{{\sc nn}}}$ as given in Eq. (86). Computing a first order critical point requires setting the derivatives of all units to zero which amounts to solving the following system of non-linear equations:

where each row (i) represents a system of $d$ equations.

Matrix formulation of system of equations Let us rewrite (91) in matrix form. Towards this end, we define

Note that $\mathbf{\Gamma}=\mathbf{\Gamma^{\top}}$ since $\gamma^{(i,j)}=\gamma^{(j,i)},\forall i,j$.

Solving the system of equations Using the notation introduced above, we can write (91) as follows

where $({\bf B}+\mathbf{\Gamma})^{\dagger}$ is the pseudo-inverse of $({\bf B}+\mathbf{\Gamma})$.

and hence the critical points of the objective are of the following type

C.2 Possible implications for deep neural networks

From Eq. (22) in the Lemma 2 we can conclude that the optimal direction of any $\hat{{\bf w}}^{(i)}$ is independent of the corresponding output weight $\theta^{(i)}$, which only affects $\hat{{\bf w}}^{(i)}$ through the scaling parameter $\hat{c}^{(i)}$. This is a very appealing property: Take a multilayer network and assume (for the moment) that all layer inputs are Gaussian. Then, Lemma 2 still holds for any given hidden layer and gives rise to a decoupling of the optimal direction of this layer with all downstream weights, which in turn simplifies the curvature structure of the network since many Hessian blocks become zero.

However, classical local optimizers such as Gd optimize both, direction and scaling, at the same time and are therefore blind to the above global property. It is thus very natural that performing optimization in the reparametrized weight space can in fact benefit from splitting the subtasks of optimizing scaling and direction in two parts, since updates in the latter are no longer sensitive to changes in the downstream part of the network. In the next section, we theoretically prove that such a decoupling accelerates optimization of weights of each individual unit in the presence of Gaussian inputs. Of course, the normality assumption is very strong but remarkably the experimental results of Section 5.3 suggest the validity of this result beyond the Gaussian design setting and thus motivate future research in this direction.

C.3 Convergence analysis

Here, we prove the convergence result restated below. See 3

According to the result of Lemma 2, all critical points of $f_{\text{{\sc nn}}}$ are aligned along the same direction as the solution of normalized least-squares. This property is similar to the objective of learning halfspaces (with Gaussian inputs) and the proof technique below therefore follows similar steps to the convergence proof of Theorem 2.

Gradient in the original parameterization Recall the gradient of $f_{\text{{\sc nn}}}$ is defined as

Gradient in the normalized parameterization: Let us now consider the gradient of $f_{\text{{\sc nn}}}$ w.r.t the normalized weights, which relates to the gradient in the original parameterization in the following way

Replacing the expression given in Eq. (86) into the above formula yields

Note that the constants $\alpha^{(i)}$, $\beta^{(i)}$ and $\gamma^{(i,j)}$ all depend on the parameters ${\bf w}^{(i)}$ and $\theta^{(j)}$ of the respective units $i$ and $j$. The orthogonality of ${\bf S}{\bf w}^{(i)}$ to ${\bf A}_{{\bf w}^{(i)}}$ (see Eq. (69)) allows us to simplify things further:

We now use the initialization of weights $\{{\bf w}^{(k)}=c_{k}{\bf S}^{-1}{\bf u}\}_{k<i}$ and $\{{\bf w}^{(j)}={\bf 0}\}_{j>i}$ into the above expression to get

where $\nabla\rho({\bf w})$ is the gradient of the normalized ordinary least squares problem (Eq. (13)), i.e.

We conclude that the global characterization property described in Eq. (22) transfers to the gradient since the above gradient aligns with the gradient of $\rho({\bf w})$.

Choice of stepsize and stopping criterion We follow the same approach used in the proof for learning halfspaces and choose a stepsize to ensure that the gradient steps on $f_{\text{{\sc nn}}}^{(i)}$ match the gradient iterates on $\rho$, i.e.

which leads to the following choice of stepsize

If $\xi_{t}=0$, then the gradient is zero. Therefore, we choose the stopping criterion as follows

Gradient norm decomposition Proposition 4 relates the ${\bf S}^{-1}$-norm of the gradient in the original space to the normalized space as follows

In the following, we will establish convergence individually in terms of $g$ and ${\bf w}$ and then use the above result to get a global result.

Convergence in scalar $g^{(i)}$ Since the smoothness property defined in Assumption 4 also holds for $f_{\text{{\sc nn}}}^{(i)}$, we can directly invoke the result of Lemma 5 to establish a convergence rate for $g$:

Directional convergence By the choice of stepsize in Eq. (98), the gradient trajectory on $f_{\text{{\sc nn}}}$ reduces to the gradient trajectory on $\rho({\bf w})$. Hence, we can establish a linear convergence in ${\bf w}^{(i)}$ by a simple modification of Eq. (84):

The assumption 3 on loss with the choice of activation function as $\tanh$ allows us to bound the scalar $\xi_{t}^{2}$:

Combined convergence bound Combining the above results concludes the proof in the following way

Appendix D EXPERIMENTAL DETAILS

Setting We consider empirical risk minimization (ERM) as a surrogate for (17) in the binary classification setting and make two different choices for $\varphi(\cdot)$:

The first resembles classical convex logistic regression when ${\bf y}_{i}\in\{-1,1\}$. The second is a commonly used non-convex, continuous approximation of the zero-one loss in learning halfspaces (Zhang et al.,, 2015)

As datasets we use the common realworld dataset a9a ($n=32^{\prime}561,d=123$) as well a synthetic data set drawn from a multivariate gaussian distribution such that ${\bf z}\sim\mathcal{N}({\bf u},{\bf S})$ ($n=1^{\prime}000,d=50$).

Methods We compare the convergence behavior of Gd and Accelerated Gradient Descent (Agd) (Nesterov,, 2013) to Batch Normalization plus two versions of Gd as well as Weight Normalization. Namely, we assess

Gdnp as stated in Algorithm 1 but with the Bisection search replaced by multiple Gradient Descent steps on $g$ (10 per outer iteration)

Batch Norm plus standard Gd which simultaneously updates ${\bf w}$ and ${\bf g}$ with one gradient step on each parameter.

Weight Normalization plus standard Gd as above. (Salimans and Kingma,, 2016)

Results The Gaussian design experiments clearly confirm Theorem 2 in the sense that the loss in the convex-, as well as the gradient norm in the non-convex case decrease at a linear rate. The results on a9a show that Gdnp can accelerate optimization even when the normality assumption does not hold and in a setting where no covariate shift is present. This motivates future research of non-linear reparametrizations even in convex optimization.

Regarding Bn and Wn we found a clear trade-off between making fast progress in the beginning and optimizing the last couple of digits. In the above results of Figure 2 and 5 we report runs with stepsizes that were optimized for the latter case but we here note that early progress can easily be achieved in normalized parametrizations (which the linear a9a softplus plot actually confirms) e.g. by putting a higher learning rate on ${\bf g}$. In the long run similar performance to that of Gd sets in, which suggests that the length-direction decoupling does not fully do the trick. The superior performance of Gdnp points out that either an increased number of steps in the scaling factor $g$ or an adaptive stepsize scheme such as the one given in Eq. (26) (or both) may significantly increase the performance of Batch Normalized Gradient Descent Bn.

It is thus an exciting open question whether such simple modifications to Gd can also speed up the training of Batch Normalized neural networks. Finally, since Gdnp performs similar to Agd in the non-gaussian setting, it is a logical next step to study how accelerated gradient methods like Agd or Heavy Ball perform in normalized coordinates.

As a side note, Figure 7 shows how surprisingly different the paths that Gradient Descent takes before and after normalization can be.

D.2 Neural networks

Setting and methods We test the validity of Theorem 3 and Lemma 2 outside the Gaussian setting and a normalized and an unnormalized feedforward networks on the CIFAR10 image classification task. This dataset consists of 60000 32x32 images in 10 classes, with 6000 images per class (Krizhevsky and Hinton,, 2009). The networks have six hidden layers with 50 hidden units in each of them. Each hidden unit has a tanh activation function, except for the very last layer which is linear. These scores are fed into a cross entropy loss layer which combines softmax and negative log likelihood loss. The experiments are implemented using the PyTorch framework (Paszke et al.,, 2017).

The first network is trained by standard Gd and the second by Gd in normalized coordinates (i.e. Bn) with the same fixed stepsize on and ${\bf w}$, but we increase the learning rate on ${\bf g}$ by a factor of 10 which accelerates training significantly. The second network thus resembles performing standard Gd in a network where all hidden layers are Batch Normalized. We measure the cross-dependency of the central with all other layers in terms of the Frobenius norm of the second partial derivatives $\frac{\partial^{2}f_{\text{{\sc nn}}}}{\partial{\bf W}_{4}\partial{\bf W}_{i}}$. This quantity signals how the gradients of layer 4 change when we alter the direction of any other layer. From an optimization perspective, this is a sound measure for the cross-dependencies: If it is close to zero (high), that means that a change in layer $i$ induces no (a large) change in layer 4. Compared to gradient calculations, computing second derivatives is rather expensive $O(nd^{2})$ (where $d=66700$), which is why we evaluate this measure every only 250 iterations.

Results Figure 3 and 8 confirm that the directional gradients of the central layer are affected far more by the upstream than by the downstream layers to a surprisingly large extent. Interestingly, this holds even before reaching a critical point. The cross-dependencies are generally decaying for the Batch Normalized network (Bn) while they remain elevated in the un-normalized network (Gd), which suggest that using Batch Normalization layers indeed simplifies the networks curvature structure in ${\bf w}$ such that the length-direction decoupling allows Gradient Descent to exploit simpler trajectories in these normalized coordinates for faster convergence. Of course, we cannot untangle this effect fully from the covariate shift reduction that was mentioned in the introduction. Yet, the fact that the (de-)coupling increases in the distance to the middle layer (note how earlier (later) layers are more (less) important for the ${\bf W}_{4}$) emphasizes the relevance of this analysis particularly for deep neural network structures, where downstream dependencies might vanish completely with depth. This does not only make gradient based training easier but also suggests the possibility of using partial second order information, such as diagonal Hessian approximations (e.g. proposed in (Martens et al.,, 2012)).