Title: DIRECTIONAL ANALYSIS OF STOCHASTIC GRADIENT DESCENT VIA VON MISES-FISHER DISTRIBUTIONS IN DEEP LEARNING

Abstract: Although stochastic gradient descent (SGD) is a driving force behind the recent success of deep learning, our understanding of its dynamics in a high-dimensional parameter space is limited. In recent years, some researchers have used the stochasticity of minibatch gradients, or the signal-to-noise ratio, to better characterize the learning dynamics of SGD. Inspired by these work, we here analyze SGD from a geometrical perspective by inspecting the stochasticity of the norms and directions of minibatch gradients. We propose a model of the directional concentration for minibatch gradients through von Mises-Fisher distribution and show that the directional uniformity of minibatch gradients increases over the course of SGD. We empirically verify our result using deep convolutional networks and observe a higher correlation between the gradient stochasticity and the proposed directional uniformity than that against the gradient norm stochasticity, suggesting that the directional statistics of minibatch gradients is a major factor behind SGD.

1. INTRODUCTION:

Stochastic gradient descent (SGD) has been a driving force behind the recent success of deep learning. Despite a series of work on improving SGD by incorporating the second-order information of the objective function (Roux et al., 2008; Martens, 2010; Dauphin et al., 2014; Martens & Grosse, 2015; Desjardins et al., 2015) , SGD is still the most widely used optimization algorithm for training a deep neural network. The learning dynamics of SGD, however, has not been well characterized beyond that it converges to an extremal point (Bottou, 1998) due to the non-convexity and highdimensionality of a usual objective function used in deep learning.

Gradient stochasticity, or the signal-to-noise ratio (SNR) of the stochastic gradient, has been proposed as a tool for analyzing the learning dynamics of SGD. Shwartz-Ziv & Tishby (2017) identified two phases in SGD based on this. In the first phase, "drift phase", the gradient mean is much higher than its standard deviation, during which optimization progresses rapidly. This drift phase is followed by the "diffusion phase", where SGD behaves similarly to Gaussian noise with very small means. Similar observations were made by Li & Yuan (2017) and Chee & Toulis (2018) who have also divided the learning dynamics of SGD into two phases. Shwartz-Ziv & Tishby (2017) have proposed that such phase transition is related to information compression. Unlike them, we notice that there are two aspects to the gradient stochasticity. One is the L 2 norm of the minibatch gradient (the norm stochasticity), and the other is the directional balance of minibatch gradients (the directional stochasticity). SGD converges or terminates when either the norm of the minibatch gradient vanishes to zeros, or when the angles of the minibatch gradients are uniformly distributed and their non-zero norms are close to each other. That is, the gradient stochasticity, or the SNR of the stochastic gradient, is driven by both of these aspects, and it is necessary for us to investigate not only the holistic SNR but also the SNR of the minibatch gradient norm and that of the minibatch gradient angles.

In this paper, we use a von Mises-Fisher (vMF hereafter) distribution, which is often used in directional statistics (Mardia & Jupp, 2009) , and its concentration parameter κ to characterize the directional balance of minibatch gradients and understand the learning dynamics of SGD from the perspective of directional statistics of minibatch gradients. We prove that SGD increases the direc-tional balance of minibatch gradients. We empirically verify this with deep convolutional networks with various techniques, including batch normalization (Ioffe & Szegedy, 2015) and residual connections (He et al., 2015) , on MNIST and CIFAR-10 (Krizhevsky & Hinton, 2009) . Our empirical investigation further reveals that the proposed directional stochasticity is a major drive behind the gradient stochasticity compared to the norm stochasticity, suggesting the importance of understanding the directional statistics of the stochastic gradient.

Contribution We analyze directional stochasticity of the minibatch gradients via angles as well as the concentration parameter of the vMF distribution. Especially, we theoretically show that the directional uniformity of the minibatch gradients modeled by the vMF distribution increases as training progresses, and verify this by experiments. In doing so, we introduce gradient norm stochasticity as the ratio of the standard deviation of the minibatch gradients to their expectation and theoretically and empirically show that this gradient norm stochasticity decreases as the batch size increases.

Related work:

Most studies about SGD dynamics have been based on two-phase behavior (Shwartz-Ziv & Tishby, 2017; Li & Yuan, 2017; Chee & Toulis, 2018) . Li & Yuan (2017) investigated this behavior by considering a shallow neural network with residual connections and assuming the standard normal input distribution. They showed that SGD-based learning under these setups has two phases; search and convergence phases. Shwartz-Ziv & Tishby (2017) on the other hand investigated a deep neural network with tanh activation functions, and showed that SGD-based learning has drift and diffusion phases. They have also proposed that such SNR transition (drift + diffusion) is related to the information transition divided into empirical error minimization and representation compression phases. However, Saxe et al. (2018) have reported that the information transition is not generally associated with the SNR transition with ReLU (Nair & Hinton, 2010) activation functions. Chee & Toulis (2018) instead looked at the inner product between successive minibatch gradients and presented transient and stationary phases.

Unlike our work here, the experimental verification of the previous work conducted under limited settings -the shallow network (Li & Yuan, 2017) , the specific activation function (Shwartz-Ziv & Tishby, 2017) , and only MNIST dataset (Shwartz-Ziv & Tishby, 2017; Chee & Toulis, 2018) -that conform well with their theoretical assumptions. Moreover, their work does not offer empirical result about the effect of the latest techniques including both batch normalization (Ioffe & Szegedy, 2015) layers and residual connections (He et al., 2015) .

2. PRELIMINARIES:

Norms and Angles Unless explicitly stated, a norm refers to L 2 norm.

• and •, • thus correspond to L 2 norm and the Euclidean inner product on R d , respectively. We use x n ⇒ x to indicate that "a random variable x n converges to x in distribution." Similarly, x n P → x means convergence in probability. An angle θ between d-dimensional vectors u and v is defined by

θ = 180 π cos -1 u,v u v

.

Loss functions A loss function of a neural network is written as f (w) = 1 n n i=1 f i (w), where w ∈ R d is a trainable parameter. f i is "a per-example loss function" computed on the i-th data point. We use I and m to denote a minibatch index set and its batch size, respectively. Further, we call f I (w) = 1 m i∈I f i (w) "a minibatch loss function given I". In Section 3.1, we use g i (w) and ĝ(w) to denote -∇ w f i (w) and -∇ w f I (w), respectively. In Section 3.3, the index i is used for the corresponding minibatch index set I i . For example, the negative gradient of f Ii (w) is written as ĝi (w). During optimization, we denote a parameter w at the i-th iteration in the t-th epoch as w i t , and w 0 0 is an initial parameter. We use n b to refer to the number of minibatches in a single epoch. 2005)) The pdf of the vMF(µ, κ) is given by f d (x; µ, κ) = C d (κ) exp(κµ x) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q x y -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q κ 0 5 50 If κ is zero, then it is a uniform distribution on the unit hypersphere, and as κ → ∞, it becomes a point mass on the unit hypersphere (Figure 1 ). The maximum likelihood estimates for µ and κ are μ = . The formula for κ is approximate since the exact computation is intractable (Banerjee et al., 2005) .

3. THEORETICAL MOTIVATION:



3.1. ANALYSIS OF THE GRADIENT NORM STOCHASTICITY:

It is a usual practice for SGD to use a minibatch gradient ĝ(w) = -∇ w f I (w) instead of a full batch gradient g(w) = -∇ w f (w). The minibatch index set I is drawn from {1,. . . ,n} randomly. ĝ(w) satisfies E[ĝ(w)] = g(w) and Cov(ĝ(w), ĝ(w)) ≈ 1 mn n i=1 g i (w)g i (w) for n m where n is the number of full data points and g i (w) = -∇ w f i (w) (Hoffer et al., 2017) . As the batch size m increases, the randomness of ĝ(w) decreases. Hence E ĝ(w) tends to g(w) , and Var( ĝ(w) ), which is the variance of the norm of the minibatch gradient, vanishes. The convergence rate analysis is as the following: Theorem 1. Let ĝ(w) be a minibatch gradient induced from the minibatch index set I of batch size m from {1, . . . , n} and suppose γ = max i,j∈{1,...,n} | g i (w), g j (w) |. Then

0 ≤ E ĝ(w) -g(w) ≤ 2(n -m) m(n -1) × γ E ĝ(w) + g(w) ≤ (n -m)γ m(n -1) g(w)

and

Var( ĝ(w) ) ≤ 2(n -m) m(n -1) γ .

Hence,

Var( ĝ(w) ) E ĝ(w) ≤ 2(n -m) m(n -1) × γ g(w) 2 .

(1)

Proof. See Supplemental A.

According to Theorem 1, a large batch size m reduces the variance of ĝ(w) centered at E ĝ(w) with convergence rate O(1/m). We empirically verify this by estimating the gradient norm stochasticity at random points while varying the minibatch size, using a fully-connected neural network (FNN) with MNIST, as shown in Figure 2 (a) (see Supplemental E for more details.) q q q q q q q q q q q q q q q q Batch size Gradient norm stochasticity This theorem however only demonstrate that the gradient norm stochasticity is (l.h.s. of ( 1)) is low at random initial points. It may blow up after SGD updates, since the upper bound (r.h.s. of ( 1)) is inversely proportional to g(w) . This implies that the learning dynamics and convergence of SGD, measured in terms of the vanishing gradient, i.e., n b i=1 ĝi (w) ≈ 0, is not necessarily explained by the vanishing norms of minibatch gradients, but rather by the balance of the directions of ĝi (w)'s, which motivates our investigation of the directional statistics of minibatch gradients. See Figure 2 (b) as an illustration.

3.2. UNIFORMITY MEASUREMENT VIA ANALYSIS OF ANGLES:

In order to investigate the directions of minibatch gradients and how they balance, we start from an angle between two vectors. First, we analyze an asymptotic behavior of angles between uniformly random unit vectors in a high-dimensional space.

Theorem 2. Suppose that u and v are mutually independent d-dimensional uniformly random unit vectors. Then,

√ d × 180 π cos -1 u, v -90 ⇒ N 0, 180 π 2 as d → ∞.

Proof. See Supplemental B.

According to Theorem 2, the angle between two independent uniformly random unit vectors is normally distributed and becomes increasingly more concentrated as d grows (Figure 3(a) ). If SGD iterations indeed drive the directions of minibatch gradients to be uniform, then, at least, the distribution of angles between minibatch gradients and a given uniformly sampled unit vector follows asymptotically

N 90, 180 π √ d 2 .

(2)

Figures 3(b) and 3(c) show that the distribution of the angles between minibatch gradients and a given uniformly sampled unit vector converges to an asymptotic distribution (2) after SGD iterations. Although we could measure the uniformity of minibatch gradients how the angle distribution between minibatch gradients is close to (2), it is not as trivial to compare the distributions as to compare numerical values. This necessitates another way to measure the uniformity of minibatch gradients. We draw a density plot θ(u, ĝj (w)) ĝj (w)) ) for 3, 000 minibatch gradients (black) at w = w 0 0 (b) and w = w 0 final , with training accuracy of > 99.9%, (c) when u is given. After SGD iterations, the density of θ(u, ĝj (w)) converges to an asymptotic density (red). The dimension of FNN is 635,200.

3.3. UNIFORMITY MEASUREMENT VIA VMF DISTRIBUTION:

To model the uniformity of minibatch gradients, we propose to use the vMF distribution in Definition 1. The concentration parameter κ measures how uniformly the directions of unit vectors are distributed. By Theorem 1, with a large batch size, the norm of minibatch gradient is nearly deterministic, and μ is almost parallel to the direction of full batch gradient. In other words, κ measures the concentration of the minibatch gradients directions around the full batch gradient.

The following Lemma 1 introduces the relationship between the norm of averaged unit vectors and κ, the approximate estimator of κ. Lemma 1. The approximated estimator of κ induced from the d-dimensional unit vectors 

{x 1 , x 2 , • • • , x n b }, κ = r(d -r2 ) 1 -r2 ,

Consider:

κ(w) = h n b i=1 p i -w p i -w ,

which is measured from the directions from the current location w to the fixed points p i 's, where h(•) is a function defined in Lemma 1. Since h(•) is an increasing function, we may focus only on

n b i=1 pi-w pi-w

to see how κ behaves with respect to its argument. Lemma 2 implies that the estimated directional concentration κ decreases if we move away from w 0 0 to w = w 0 0 + i pi-w 0 0 pi-w 0 0 with a small (Figure 4 (a)). In other words, κ(w ) < κ(w 0 0 ). Lemma 2. Let p 1 , p 2 , • • • , p n b be d-dimensional vectors. If all p i 's are not on a single ray from the current location w, then there exists positive number η such that )/ p i -w 0 0 and x i = (p i -w )/ p i -w , then i x i < i x i which is equivalent to κ(w ) < κ(w 0 0 ); (b) If ĝi (w 0 0 )'s are sufficiently parallel to (p i -w 0 0 )'s for each i, then w ≈ w 0 0 + η i ĝi (w 0 0 ). When w and w 0 1 are sufficiently close to each other, we also have κ(w 0 1 ) < κ(w 0 0 ).

n b j=1 p j -w - n b i=1 pi-w pi-w p j -w - n b i=1 pi-w pi-w < n b i=1 p i -w p i -w

for all ∈ (0, η].

Proof. See Supplemental C.2.

We make the connection between the observation above and SGD by first viewing p i 's as local minibatch solutions.

Definition 2. For a minibatch index set I i , p i (w) = arg min w ∈N (w;ri) f Ii (w ) is a local minibatch solution of I i at w, where N (w; r i ) is a neighborhood of radius r i at w. Here, r i is determined by w and I i for p i (w) to exist uniquely.

Under this definition, p i (w) is local minimum of a minibatch loss function f Ii near w. Then we reasonably expect that the direction of ĝi (w) = -∇ w f Ii (w) is similar to that of p i (w) -w.

Each epoch of SGD with a learning rate η computes a series of w j t = w 0 t + η j i=1 ĝi (w i-1 t ) for all j ∈ {1, . . . , n b }. If ĝi (•)'s are Lipschitz continuous for all i ∈ {1, . . . , n b }, then we have ĝi (w i-1 t ) ≈ ĝi (w 0 t ) for a small η. Moreover, Theorem 1 implies ĝi (w 0 t ) ≈ τ for all i ∈ {1, . . . , n b } with a large batch size or at the early stage of SGD iterations. Combining these approximations, ĝi (w i-1 t ) ≈ τ for all i ∈ {1, . . . , n b }.

For example, suppose that t = 0, n b = 3 and τ = 1, and assume that p i (w 0 0 ) = p i (w 0 1 ) = p i for all i = 1, 2, 3. Then,

EQUATION

κ(w 0 1 ) ≈ h    3 j=1 p j -w 0 0 -η 3 i=1 ĝi(w i-1 0 ) ĝi(w i-1 0 ) p j -w 0 0 -η 3 i=1 ĝi(w i-1 0 ) ĝi(w i-1 0 )    . If η i pi-w 0 0 pi-w 0 0 ≈ η i ĝi(w i-1 0 ) ĝi(w i-1 0 ) , then 3 j=1 p j -w 0 0 -η 3 i=1 ĝi(w i-1 0 ) ĝi(w i-1 0 ) p j -w 0 0 -η 3 i=1 ĝi(w i-1 0 ) ĝi(w i-1 0 ) ≈ 3 j=1 p j -w 0 0 -η 3 i=1 pi-w 0 0 pi-w 0 0 p j -w 0 0 -η 3 i=1 pi-w 0 0 pi-w 0 0 .

Hence, we have κ(w 0 1 ) < κ(w 0 0 ) by Lemma 2. A trivial case satisfying this condition would be for each pair of p i -w 0 0 and ĝi (w i-1 0 ) to be approximately parallel, as illustrated in Figure 4 (b).

Theorem 3. Let p 1 (w 0 t ), p 2 (w 0 t ), • • • , p n b (w 0 t ) be d-dimensional vectors,

and all p i (w 0 t )'s are not on a single ray from the current location w 0 t . If

n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t - n b i=1 ĝi (w i-1 t ) ĝi (w i-1 t ) ≤ ξ (3)

for a sufficiently small ξ > 0, then there exists positive number η such that

n b j=1 p j (w 0 t ) -w 0 t - n b i=1 ĝi(w i-1 t ) ĝi(w i-1 t ) p j (w 0 t ) -w 0 t - n b i=1 ĝi(w i-1 t ) ĝi(w i-1 t ) < n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t (4)

for all ∈ (0, η].

Proof. See Supplemental C.3.

This Theorem 3 asserts that κ(•) decreases even with some perturbation along the averaged direction i pi(w)-w pi(w)-w . With additional assumptions on each minibatch loss functions, we have a sufficient condition for (3), summarized in Corollary 3.1. Corollary 3.1. Let p i be the local minibatch solution of each f Ii . Suppose a region R satisfying:

For all w, w ∈ R, p i (w) = p i (w ) = p i for all i = 1, • • • , n b .

Further, assume that Hessian matrices of f Ii 's are positive definite, wellconditioned, and bounded in the sense of matrix L 2 -norm on R. If SGD moves from w 0 t to w 0 t+1 on R with a large batch size and a small learning rate, then κ(w 0 t ) > κ(w 0 t+1 ). Moreover, we can estimate κ(w 0 t ) and κ(w 0 t+1 ) by minibatch gradients at w 0 t and w 0 t+1 , respectively.

Proof. See Supplemental D.

Without the corollary above, we need to solve p i (w 0 t ) = arg min w∈N (w 0 t ;r) f Ii (w) for all i ∈ {1, . . . , n s }, where n s is the number of samples to estimate κ, in order to compute κ(w 0 t ). Corollary 3.1 however implies that we can compute κ(w 0 t ) by using ĝi(w 0 t ) ĝi(w 0 t ) instead of pi(w 0 t )-w 0 t pi(w 0 t )-w 0 t , significantly reducing computational overhead.

In Practice Although the number of all possible minibatches in each epoch is n b = n m , it is often the case to use n b ≈ n/m minibatches at each epoch in practice to go from w 0 t to w 0 t+1 . Assuming that these n b minibatches were selected uniformly at random, the average of the n b normalized minibatch gradients is the maximum likelihood estimate of µ, just like the average of all n b normalized minibatch gradients. Thus, we expect with a large n b ,

( n m ) i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t - n b i=1 ĝi (w i-1 t ) ĝi (w i-1 t ) ≤ ξ,

and that SGD in practice also satisfies κ(w 0 t ) > κ(w 0 t+1 ).

4. EXPERIMENTS:



4.1. SETUP:

In order to empirically verify our theory on directional statistics of minibatch gradients, we train various types of deep neural networks using SGD and monitor the following metrics for analyzing the learning dynamics of SGD:

• Training loss Epoch κ ^1 5 20 100 1000 10 4 10 5 10 6 q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q (a) Batch size: 64

Epoch κ ^1 5 20 100 1000 10 4 10 5 10 6 10 7 q q q q q q q q q q q q q q qq qq q qq q qq qq qq q q q q qq q q q q q q q q q q q q q q q q (b) Batch size: 256

Epoch κ ^1 5 20 100 1000 10 4 10 5 10 6 10 7 q q q q q q q q q q q q q q q q qq q q q q qq qq qq q q q q qq q q q q q q q q q q q q q q q q (c) Batch size: 1,024

Epoch κ ^1 5 20 100 1000 10 5 10 6 10 7 10 8 q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q (d) Batch size: 4,096

Figure 5 : We show the average κ (black curve) ± std. (shaded area), as the function of the number of training epochs (in log-log scale) across various batch sizes in MNIST classifications using FNN with fixed learning rate 0.01 and 5 random initializations. Although κ with the large batch size decreases more smoothly rather than the small batch size, we observe that κ still decreases well with minibatches of size 64. We did not match the ranges of the y-axes across the plots to emphasize the trend of monotonic decrease.

• Validation loss

• Gradient stochasticity (GS)↑ E∇ w f Ii (w) / tr(Cov(∇ w f Ii (w), ∇ w f Ii (w)) ↓ • Gradient norm stochasticity (GNS) ↑ E ∇ w f Ii (w) / Var( ∇ w f Ii (w) ) ↓ • Directional Uniformity↑ κ ↓

The latter three quantities are statistically estimated using n s = 3, 000 minibatches. We use κ to denote the κ estimate.

We train the following types of deep neural networks (Supplemental E):

• FNN: a fully connected network with a single hidden layer • DFNN: a fully connected network with three hidden layers • CNN: a convolutional network with 14 layers (Krizhevsky et al., 2012) In the case of the CNN, we also evaluate its variant with skip connections (+Res) (He et al., 2015) .

As it was shown recently by Santurkar et al. (2018) that batch normalization (Ioffe & Szegedy, 2015) improves the smoothness of a loss function in terms of its Hessian, we also test adding batch normalization to each layer right before the ReLU (Nair & Hinton, 2010 ) nonlinearity (+BN). We use MNIST for the FNN, DFNN and their variants, while CIFAR-10 ( Krizhevsky & Hinton, 2009) for the CNN and its variants.

Our theory suggests a sufficiently large batch size for verification. We empirically analyze how large a batch size is needed in Figure 5 . From these plots, κ decreases monotonically regardless of the minibatch size, but the variance over multiple training runs is much smaller with a larger minibatch size. We thus decide to use a practical size of 64. With this fixed minibatch size, we use a fixed learning rate of 0.01, which allows us to achieve the training accuracy of > 99.9% for every training run in our experiments. We repeat each setup five times starting from different random initial parameters and report both the mean and standard deviation.

4.2. DIRECTIONAL UNIFORMITY INCREASES:

FNN and DFNN We first observe that κ decreases over training regardless of the network's depth in Figure 6 (a, b ). We however also notice that κ decrease monotonically with the FNN, but less so with its deeper variant (DFNN). We conjecture this is due to the less-smooth loss landscape of a deep neural network. This difference between FNN and DFNN however almost entirely vanishes when batch normalization (+BN) is applied (Figure 6 (e,f)). This was expected as batch normalization is known to make the loss function behave better, and our theory assumes a smooth objective function.

CNN The CNN is substantially deeper than either FNN or DFNN and is trained on a substantially more difficult problem of CIFAR-10. In other words, the assumptions underlying our theory may not hold as well. Nevertheless, as shown in Figure 6 (c), κ eventually drops below its initial point, although this trend is not monotonic and κ fluctuates significantly over training. The addition of batch normalization (+BN) helps with the fluctuation but κ does not monotonically decrease q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Epoch κ ^1

5 20 50 150 2 × 5 × 2 × 5 × (a) FNN

MNIST 98.18 ± 0.04% q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Epoch κ ^1 5 20 50 150 10 5 2 × 10 5 5 × 10 5 10 6 2 × 10 6 (b) DFNN MNIST 98.30 ± 0.05% q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Epoch κ ^1 (c) CNN CIFAR-10 53.78 ± 3.22% q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q Epoch κ ^1 CIFAR-10 63.39 ± 1.68% q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Epoch κ ^1

5 20 50 150 2 × 5 × 2 × 5 × (e) FNN+BN

MNIST 98.37 ± 0.08% q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Epoch κ ^1 (f) DFNN+BN MNIST 98.51 ± 0.11% q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q Epoch κ ^1 (g) CNN+BN CIFAR-10 73.00 ± 1.16% q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q Epoch κ ^1 Effect of +BN and +Res Based on our observations that the uniformity of minibatch gradients increases monotonically, when a deep neural network is equipped with residual connection (+Res) and trained with batch normalization (+BN), we conjecture that the loss function induced from these two techniques better satisfies the assumptions underlying our theoretical analysis, such as its well-behavedness. This conjecture is supported by for instance Santurkar et al. (2018) , who demonstrated batch normalization guarantees the boundedness of Hessian, and Orhan & Pitkow (2017), who showed residual connections eliminate some singularities of Hessian.

κ near the end of training The minimum average κ of DFNN+BN, which has 1, 920, 000 parameters, is 71, 009.20, that of FNN+BN, which has 636, 800 parameters, is 23, 059.16, and that of CNN+BN+Res, which has 207, 152 parameters, is 20, 320.43. These average κ are within a constant multiple of estimated κ using 3,000 samples from the vMF distribution with true κ = 0 (35, 075.99 with 1, 920, 000 dimensions, 11, 621.63 with 636, 800 dimensions, and 3, 781.04 with 207, 152 dimensions.) This implies that we cannot say that the underlying directional distribution of minibatch gradients in all these cases at the end of training is not close to uniform (Cutting et al., 2017) . For more detailed analysis, see Supplementary F.

4.3. DIRECTIONAL UNIFORMITY AND OTHER METRICS:

The gradient stochasticity (GS) was used by Shwartz-Ziv & Tishby (2017) as a main metric for identifying two phases of SGD learning in deep neural networks. This quantity includes both the gradient norm stochasticity (GNS) and the directional uniformity κ, implying that either or both of GNS and κ could drive the gradient stochasticity. We thus investigate the relationship among these three quantities as well as training and validation losses. We focus on CNN, CNN+BN and CNN+Res+BN trained on CIFAR-10. 1.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.1 (a) CNN q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (b) CNN+BN q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q qq q q q q q q q SNR κ ^0.05 These plots suggest that the SNR is largely driven by the directional uniformity.

From Figure 7 (First row), it is clear that the proposed metric of directional uniformity κ correlates better with the gradient stochasticity than the gradient norm stochasticity does. This was especially prominent during the early stage of learning, suggesting that the directional statistics of minibatch gradients is a major explanatory factor behind the learning dynamics of SGD. This difference in correlations is much more apparent from the scatter plots in Figure 7 (Second row). We show these plots created from other four training runs per setup in Supplemental G.

5. CONCLUSION:

Stochasticity of gradients is a key to understanding the learning dynamics of SGD (Shwartz-Ziv & Tishby, 2017) and has been pointed out as a factor behind the success of SGD (see, e.g., LeCun et al., 2012; Keskar et al., 2016) . In this paper, we provide a theoretical framework using von Mises-Fisher distribution, under which the directional stochasticity of minibatch gradients can be estimated and analyzed, and show that the directional uniformity increases over the course of SGD. Through the extensive empirical evaluation, we have observed that the directional uniformity indeed improves over the course of training a deep neural network, and that its trend is monotonic when batch normalization and skip connections were used. Furthermore, we demonstrated that the stochasticity of minibatch gradients is largely determined by the directional stochasticity rather than the gradient norm stochasticity.

Our work in this paper suggests two major research directions for the future. First, our analysis has focused on the aspect of optimization, and it is an open question how the directional uniformity relates to the generalization error although handling the stochasticity of gradients has improved SGD (Neelakantan et al., 2015; Hoffer et al., 2017; Smith et al., 2017; Jin et al., 2017) . Second, we have focused on passive analysis of SGD using the directional statistics of minibatch gradients, but it is not unreasonable to suspect that SGD could be improved by explicitly taking into account the directional statistics of minibatch gradients during optimization.

SUPPLEMENTARY MATERIAL A PROOFS FOR THEOREM 1:

In proving Theorem 1, we use Lemma A.1. Define selector random variables (Hoffer et al., 2017) as below:

s i = 1, if i ∈ I 0, if i / ∈ I .

Then we have

ĝ(w) = 1 m n i=1 g i (w)s i .

Lemma A.1. Let ĝ(w) be a minibatch gradient induced from the minibatch index set I with batch size m from {1, . . . , n}. Then

EQUATION

)

where γ = max i,j∈{1,...,n} | g i (w), g j (w) |.

Proof. By Jensen's inequality, 0 ≤ E ĝ(w) 2 -g(w) 2 . Note that

E ĝ(w) 2 = n i=1 n j=1 1 m 2 g i (w), g j (w) E[s i s j ]. Since E[s i s j ] = m n δ ij + m(m-1) n(n-1) (1 -δ ij ), E ĝ(w) 2 -g(w) 2 = 1 mn - m -1 mn(n -1) n i=1 g i (w), g i (w) + m -1 mn(n -1) - 1 n 2 n i=1 n j=1 g i (w), g j (w) = 1 mn - m -1 mn(n -1) n i=1 g i (w), g i (w) + m -n mn 2 (n -1) n i=1 n j=1 g i (w), g j (w) ≤ 1 mn - m -1 mn(n -1) nγ + n -m mn 2 (n -1) n 2 γ = 2(n -m) m(n -1) γ

where γ = max i,j∈{1,...,n} | g i (w), g j (w) |.

Theorem 1. Let ĝ(w) be a minibatch gradient induced from the minibatch index set I of batch size m from {1, . . . , n} and suppose γ = max i,j∈{1,...,n} | g i (w), g j (w) |. Then

0 ≤ E ĝ(w) -g(w) ≤ 2(n -m) m(n -1) × γ E ĝ(w) + g(w) ≤ (n -m)γ m(n -1) g(w)

and

Var( ĝ(w) ) ≤ 2(n -m) m(n -1) γ .

Hence, Var( ĝ(w) )

E ĝ(w) ≤ 2(n -m) m(n -1) × γ g(w) 2 .

Proof. By Jensen's inequality, we have g(w) = E[ĝ(w)] ≤ E ĝ(w) and (E ĝ(w) ) 2 ≤ E ĝ(w) 2 . From the second inequality and Lemma A.1,

(E ĝ(w) ) 2 ≤ E ĝ(w) 2 ≤ g(w) 2 + 2(n -m) m(n -1) γ or E ĝ(w) -g(w) E ĝ(w) + g(w) ≤ 2(n -m) m(n -1) γ . Hence E ĝ(w) ≤ g(w) + 2(n -m) m(n -1) × γ E ĝ(w) + g(w) ≤ (n -m)γ m(n -1) g(w)

.

Further,

Var( ĝ(w) ) = E ĝ(w) 2 -(E ĝ(w) ) 2 ≤ E ĝ(w) 2 -Eĝ(w) 2 ≤ g(w) 2 + 2(n -m) m(n -1) γ -g(w) 2 = 2(n -m) m(n -1) γ.

B PROOFS FOR THEOREM 2:

For proofs, Slutsky's theorem and delta method are key results to describe limiting behaviors of random variables in distributional sense. Theorem B.1. (Slutsky's theorem, Casella & Berger (2002) ) Let {x n }, {y n } be a sequence of random variables that satisfies x n ⇒ x and y n P → ρ when n goes to infinity and ρ is constant. Then x n y n ⇒ cx Theorem B.2. (Delta method, Casella & Berger (2002) ) Let y n be a sequence of random variables that satisfies √ n(y n -µ) ⇒ N (0, σ 2 ). For a given smooth function f : R → R, suppose that f (µ) exists and is not 0 where f is a derivative. Then √ n × (f (y n ) -f (µ)) ⇒ N (0, σ 2 (f (µ)) 2 ). Lemma B.1. Suppose that u and v are mutually independent d-dimensional uniformly random unit vectors. Then, √ d u, v ⇒ N (0, 1) as d → ∞.

Proof. Note that d-dimensional uniformly random unit vectors u can be generated by normalization of d-dimensional multivariate standard normal random vectors x ∼ N (0,

I d ). That is, u ∼ x x .

Suppose that two independent uniformly random unit vector u and v are generated by two independent d-dimensional standard normal vector

x = (x 1 , x 2 , • • • , x d ) and y = (y 1 , y 2 , • • • , y d ). Denote them u = x x and v = y y .

By SLLN, we have x i y i ⇒ N (0, 1).

x √ d → 1 a.s. (Use 1 d d i=1 x 2 i → Ex 2 1 = 1

Therefore, by Theorem B.1 (Slutsky's theorem), √ d u, v ⇒ N (0, 1).

Theorem 2. Suppose that u and v are mutually independent d-dimensional uniformly random unit vectors. Then,

√ d × 180 π cos -1 u, v -90 ⇒ N 0, 180 π 2 as d → ∞.

Proof. Suppose that µ = 0, σ = 1, and Recall that any continuous function on the compact interval [0, n b (1 -)] is bounded. Hence the derivative of κ with respect to u is bounded. This implies the Lipschitz continuity of h(•).

f (•) = 180 π cos -1 (•). Since 180 π d dx cos -1 (x) = -180 π √ 1-x 2 , we have f (µ) = -

h(•) is strictly increasing since r = u n b . Further,

h (u) = 1 n b 2 d 2 κ dr 2 = 2r 5 + (4 -8d)r 3 + (8d -6)r n b 2 (1 -r2 ) 4 > 0 due to r ∈ [0, 1]. Therefore h (•) is also increasing on [0, n b ). C.2 PROOF OF LEMMA 2 Lemma 2. Let p 1 , p 2 , • • • , p n b be d-dimensional vectors.

If all p i 's are not on a single ray from the current location w, then there exists positive number η such that

n b j=1 p j -w - n b i=1 pi-w pi-w p j -w - n b i=1 pi-w pi-w < n b i=1 p i -w p i -w for all ∈ (0, η].

Proof. Without loss of generality, we regard w as the origin. Let

f ( ) = n b j=1 pj - n b i=1 p i p i pj - n b i=1 p i p i 2 , then f (0) = n b i=1 pi pi 2

. Therefore, we only need to show

f (0) < 0. Now denote x j = pj pj , p j ( ) = p j - n b i=1 x i and u = - n b i=1 x i . That is, p j ( ) = p j + u. Since f ( ) = n b j=1 p j ( ) p j ( ) , n b j=1 p j ( ) p j ( ) ,

we have

f ( ) = 2 n b j=1 p j ( ) p j ( ) , d d n b j=1 p j ( ) p j ( ) and d d n b j=1 p j ( ) p j ( ) = n b j=1 p j ( ) u - u,pj ( ) pj ( ) p j ( ) p j ( ) 2 . Hence f ( ) = 2 n b j=1 p j ( ) p j ( ) , n b j=1 p j ( ) u - u,pj ( ) pj ( ) p j ( ) p j ( ) 2 .

Note that p j (0) = p j and x j = 1. We have

f (0) = 2 n b j=1 p j p j , n b j=1 p j u - u,pj pj p j p j 2 = 2 n b j=1 p j p j , n b j=1 1 p j u -u, p j p j p j p j = 2 n b j=1 x j , n b j=1 1 p j u -u, x j x j = 2 -u, n b j=1 1 p j u -u, x j x j = -2 n b j=1 u 2 -u, x j 2 p j ≤ -2 n b j=1 u 2 -u 2 x j 2 p j = 0

Since the equality holds when u, x j 2 = u 2 x j 2 for all j, we have strict inequality when all p i 's are not located on a single ray from the origin.

C.3 PROOF OF THEOREM 3:

The proof of Theorem 3 is very similar to that of Lemma 2. Theorem 3. Let p 1 (w 0 t ), p 2 (w 0 t ), • • • , p n b (w 0 t ) be d-dimensional vectors, and all p i (w 0 t )'s are not on a single ray from the current location w 0 t . If

n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t - n b i=1 ĝi (w i-1 t ) ĝi (w i-1 t ) ≤ ξ (6)

for sufficiently small ξ > 0, then there exists positive number η such that

n b j=1 p j (w 0 t ) -w 0 t - n b i=1 ĝi(w i-1 t ) ĝi(w i-1 t ) p j (w 0 t ) -w 0 t - n b i=1 ĝi(w i-1 t ) ĝi(w i-1 t ) < n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t (7)

for all ∈ (0, η].

Proof. We regard w 0 t as the origin 0. For simplicity, write p i (0) and ĝi (w i-1 t ) as p i and ĝi , respectively. Let f ( ) =

n b j=1 pj - n b i=1 p i p i pj - n b i=1 p i p i 2 and f ( ) = n b j=1 pj - n b i=1 ĝi ĝi | pj - n b i=1 ĝi ĝi 2 . Denote u = - n b j=1 pi pi , t = n b i=1 pi pi - n b i=1

ĝi ĝi and pj ( ) = p j + (u + t). Then

f ( ) = n b i=1 pj ( ) pj ( ) 2 .

Now we differentiate f ( ) with respect to , that is,

f ( ) = 2 n b j=1 pj ( ) pj ( ) , n b j=1 pj ( ) (u + t) - u+t, pj ( ) pj ( ) pj ( ) pj ( ) 2 .

Recall that pj (0) = p j . Rewrite pj pj = x j and use f (0) in the proof of Lemma 2

f (0) = 2 n b j=1 p j p j , n b j=1 p j (u + t) - u+t,pj pj p j p j 2 = 2 -u, n b j=1 u + t -u + t, pj pj pj pj p j = 2 -u, n b j=1 1 p j u + t -u + t, x j x j = 2 -u, n b j=1 1 p j u -u, x j x j + 2 -u, n b j=1 1 p j t -t, x j x j = f (0) -2 n b j=1 1 p j u, t -t, x j u, x j

Since f (0) < 0 by the proof of Lemma 2,

f (0) < 0 ⇐⇒ 2 n b j=1 1 p j t, x j u, x j -u, t < |f (0)|.

By using x j = 1 and applying the Cauchy inequality,

2 n b j=1 1 p j t, x j u, x j -u, t ≤ 2 n b j=1 t x j u x j + u t p j = 4 n b j=1 u t p j ≤ 4n b u t min j p j ≤ 4n 2 b ξ min j p j ∵ u ≤ n b i=1 x i = n b . Define r = min j p j . If ξ < |f (0)|r 4n 2 b , then f ( ) < 0.

D PROOFS FOR COROLLARY 3.1

Corollary 3.1. Let p i be local minibatch solutions for each f Ii . Suppose a region R satisfying:

For all w, w ∈ R, p i (w) = p i (w ) = p i for all i = 1, • • • , n b . Further, assume that Hessian matrices of f Ii 's are positive definite, wellconditioned, and bounded in the sense of matrix L 2 -norm on R. If SGD moves w 0 t to w 0 t+1 on R with a large batch size and a small learning rate, then κ(w 0 t ) > κ(w 0 t+1 ). Moreover, we can estimate κ(w 0 t ) and κ(w 0 t+1 ) by minibatch gradients at w 0 t and w 0 t+1 , respectively.

Proof. Recall that

w 0 t+1 = w 0 t + η n b i=1 ĝi (w i-1 t )

where η is a learning rate. To prove Corollary 3.1, we need to show κ(w 0 t+1 ) < κ(w 0 t ) which is equivalent to Bottou, 2010) . If the batch size is sufficiently large and the learning rate η is sufficiently small, ĝi (w i-1 t ) ≈ ĝi (w 0 t ) ≈ τ for all i by Theorem 1. Therefore, we have

EQUATION

) Since ∇ w 2 f Ii (•) 2 is bounded on R, ∇ w f Ii (•) is Lipschitz continuous on R(

η n b i=1 ĝi (w i-1 t ) ≈ τ η n b i=1 ĝi (w i-1 t ) ĝi (w i-1 t )

If we denote τ η as , we can convert (8) to (9).

EQUATION

)

Since both w 0 t+1 and w 0 t are in R for small learning rate, we have p i (w 0 t+1 ) = p i (w 0 t ) = p i by the assumption. That is, ( 9) is equivalent to (7). In (7), ĝi (w i-1 t )/ ĝi (w i-1 t ) cannot be replaced by (p i -w 0 t )/ p i -w 0 t in general. Hence we introduce Definition D.1 and Lemma D.1 to connect the direction of the minibatch gradient with the corresponding local minibatch solution.

Definition D.1. The condition number c(A) of a matrix A is defined as

c(A) = σ max (A) σ min (A)

where σ max (A) and σ min (A) are maximal and minimal singular values of A, respectively. If A is positive-definite matrix, then

c(A) = λ max (A) λ min (A) .

Here λ max (A) and λ min (A) are maximal and minimal eigenvalues of A, respectively.

Lemma D.1. If the condition number of the positive definite Hessian matrix of f Ii at a local minibatch solution p i , denoted by H i = ∇ w 2 f Ii (p i ), is close to 1 (well-conditioned), then the direction to p i from w is approximately parallel to its negative gradient at w. That is, for all w ∈ R,

p i -w p i -w - ĝi (w) ĝi (w) ≈ 0

where ĝi (w) = -∇ w f Ii (w).

Proof. By the second order Taylor expansion,

f Ii (w) ≈ f Ii (p i ) + 1 2 (w -p i ) H i (w -p i ). Hence, ĝi (w) = -∇ w f Ii (w) ≈ -H i (w -p i )

Denote p i -w as x. Then, we only need to show

x x - H i x H i x 2 ≈ 0

Since H i is positive definite, we can diagonalize it as H i = P i Λ i P i where P i is an orthonormal transition matrix for H i .

x x - H i x H i x 2 = 2 -2 x H i x x H i x = 2 -2 (P i x) Λ i P i x P i x P i Λ i P i x = 2 -2 (P i x) Λ i P i x P i x Λ i P x ≤ 2 -2 j λ j (P i x) 2 j P i x Λ i P i x ≤ 2 -2 λ min P i x 2 P i x λ max P i x = 2 -2 λ min λ max ≈ 0

Lemma D.1 suggests that a well-conditioned Hessian matrix of f Ii at p i allows ĝi (w)/ ĝi (w) to be replaced by (p i -w)/ p i -w for all w ∈ R. Using this, we prove Lemma D.2.

Lemma D.2. Let w be a parameter in R. If the condition number of Hessian matrix of f Ii is sufficiently close to 1 (well-conditioned) and w-w 0 t pi-w 0 t is sufficiently close to 0, then

p i -w 0 t p i -w 0 t - ĝi (w) ĝi (w) ≤ ξ n b

for all sufficiently small ξ.

Proof. We have

p i -w 0 t p i -w 0 t - ĝi (w) ĝi (w) ≤ p i -w 0 t p i -w 0 t - g(w) g(w) + g(w) g(w) - ĝi (w) ĝi (w) ≤ p i -w 0 t p i -w 0 t - g(w) g(w) + g(w) g(w) - p i -w 0 t p i -w 0 t + p i -w 0 t p i -w 0 t - p i -w p i -w + p i -w p i -w - ĝi (w) ĝi (w) ≤ + + 2 1 - p i -w 0 t , p i -w p i -w 0 t p i -w + = 3 + 2 1 - p i -w 0 t 2 -p i -w 0 t , w -w 0 t p i -w 0 t p i -w

for sufficently small (See Lemma D.1). Now we only need to show

EQUATION

Since w-w 0 t pi-w 0 t is sufficiently small, we have

p i -w 0 t 2 -p i -w 0 t , w -w 0 t = p i -w 0 t p i -w 0 t - p i -w 0 t p i -w 0 t , w -w 0 t ≥ p i -w 0 t ( p i -w 0 t -w -w 0 t ) = p i -w 0 t 2 1 - w -w 0 t p i -w 0 t

≥ 0. By using the above non-negativeness, we have the following inequality.

EQUATION

As 11) is monotonically increasing to 1. This implies that (10) holds for sufficiently small . With a small learning rate, w i-1 t 's are in R for all i ∈ {1, . . . , n b }. As a result, by Lemma D.2, we have

w-w 0 t pi-w 0 t → 0 + , (

EQUATION

for sufficiently small ξ. This implies (6) since

n b i=1 p i -w p i -w - n b i=1 ĝi (w i-1 t ) ĝi (w i-1 t ) ≤ n b i=1 p i -w p i -w - ĝi (w i-1 t ) ĝi (w i-1 t )

.

Then we can apply Theorem 3 and κ(w 0 t ) > κ(w 0 t+1 ) holds. For the last statement, "Moreover, we can estimate κ(w 0 t ) and κ(w 0 t+1 ) by minibatch gradients at w 0 t and w 0 t+1 , respectively.", recall that

κ(w 0 t ) = h n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t

where h(•) is increasing and Lipschitz continuous(Lemma 1). By Lemma D.1, we have

p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t - ĝi (w 0 t ) ĝi (w 0 t ) < ξ n b for sufficiently small ξ > 0. Therefore, n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t - n b i=1 ĝi (w 0 t ) ĝi (w 0 t ) ≤ n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t - n b i=1

ĝi (w 0 t ) ĝi (w 0 t ) where rhs is bounded by ξ. Hence, Lipschitz continuity of h(•) implies that

h n b i=1 p i (w 0 t ) -w 0 t p i (w 0 t ) -w 0 t -h n b i=1 ĝi (w 0 t ) ĝi (w 0 t ) → 0 as ξ → 0. That is, κ(w 0 t ) ≈ h n b i=1 ĝi (w 0 t ) ĝi (w 0 t )

.

Since t is arbitrary, we can apply this for all w ∈ R including w 0 t+1 .

E EXPERIMENTAL DETAILS E.1 MODEL ARCHITECTURE:

For all cases, their weighted layers do not have biases, and dropout (Srivastava et al., 2014) is not applied. We use Xavier initializations (Glorot & Bengio, 2010) and cross entropy loss functions for all experiments.

FNN The FNN is a fully connected network with a single hidden layer. It has 800 hidden units with ReLU (Nair & Hinton, 2010) activations and a softmax output layer.

DFNN The DFNN is a fully connected network with three hidden layers. It has 800 hidden units with ReLU activations in each hidden layers and a softmax output layer.

CNN The network architecture of CNN is similar to the network introduced in He et al. ( 2016) as a CIFAR-10 plain network. The first layer is 3 × 3 convolution layer and the number of output filters are 16. After that, we stack of {4, 4, 3, 1} layers with 3 × 3 convolutions on the feature maps of sizes {32, 16, 8, 4} and the numbers of filters {16, 32, 64, 128}, respectively. The subsampling is performed with a stride of 2. All convolution layers are activated by ReLU and the convolution part ends with a global average pooling (Lin et al., 2013) , a 10-way fully-conneted layers, and softmax. Note that there are 14 stacked weighted layers.

+BN We apply batch normalization right before the ReLU activations on all hidden layers.

+Res The identity skip connections are added after every two convolution layers before ReLU nonlinearity (After batch normalization, if it is applied on it.). We concatenate zero padding slices backwards when the number of filters increases.

E.2 DATA:

We use neither data augmentations nor preprocessings except scaling pixel values into [0, 1] both MNIST and CIFAR-10. In the case of CIFAR-10, for validation, we randomly choose 5000 images out of 50000 training images.

q q q q q q q q q q q q q κ κ ^0 100 500 4,000 200 500 1,000 2,000 4,000 q Sample size 1,000 2,000 3,000

Figure 8 : We show κ estimated from {1, 000 (black), 2, 000 (blue), 3, 000 (red)} random samples of the vMF distribution with underlying true κ in 10, 000-dimensional space, as the function of κ (in log-log scale except 0). For large κ, it is well-estimated by κ regardless of sample sizes. When the true κ approaches 0, we need a larger sample size to more accurately estimate this. F SOME NOTES ABOUT THE κ ESTIMATE We point out that, for a small κ, the absolute value of κ is not a precise indicator of the uniformity due to its dependence on the dimensionality, as was investigated earlier by Cutting et al. (2017) . In order to verify this claim, we run some simulations. First, we vary the number of samples and the true underlying κ with the fixed dimensionality (Unfortunately, we could not easily go over 10, 000 dimensions due to the difficulty in sampling from the vMF distribution with positive κ.). We draw {1, 000, 2, 000, 3, 000} random samples from the vMF distribution with the designated κ. We compute κ from these samples.

As can be seen from Figure 8 , the κ approaches the true κ from above as the number of samples increases. When the true κ is large, the estimation error rapidly becomes zero as the number of samples approaches 3, 000. When the true κ is low, however, the gap does not narrow completely even with 3, 000 samples.

While fixing the true κ to 0 and the number of samples to {1, 000, 2, 000, 3, 000}, we vary the dimensionality to empirically investigate the κ. We choose to use 3, 000 samples to be consistent with our experiments in this paper. We run five simulations each and report both mean and standard deviation (Table 1 ).

We clearly observe the trend of increasing κ's with respect to the dimensions. This suggests that we should not compare the absolute values of κ's across different network architectures due to the differences in the number of parameters. This agrees well with Cutting et al. (2017) which empirically showed that the threshold for rejecting the null hypothesis of κ = p by using κ where p is a fixed value grows with respect to the dimensions.

G OTHER FOUR TRAINING RUNS IN FIGURE 7:

We show plots from other four training runs in Figure 6 . For all runs, the curves of GS (inverse of SNR) and κ are strongly correlated while GNS (inverse of normSNR) is less correlated to GS. (a) CNN, initialization 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.1 (c) CNN+BN, initialization 1 q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (e) CNN+Res+BN, initialization 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (a) CNN, initialization 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q SNR κ ^0.1 (c) CNN+BN, initialization 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (e) CNN+Res+BN, initialization 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (a) CNN, initialization 3 q q q q qq qq q q q q q q q q q qq q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.1 (c) CNN+BN, initialization 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (e) CNN+Res+BN, initialization 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (a) CNN, initialization 4 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.1 (c) CNN+BN, initialization 4 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 (e) CNN+Res+BN, initialization 4 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SNR κ ^0.05 

Epoch:



:



Figures and Tables:

von Mises-Fisher Distribution We use the von Mises-Fisher (vMF) distribution to model the directions of vectors. The definition of the vMF distribution is as follows: Definition 1. (von Mises-Fisher Distribution, Banerjee et al. (

Figure 1: Characteristics of the vMF distribution in a 2-dimensional space. 100 random samples are drawn from vMF(µ, κ) where µ = (1, 0) and κ = {0, 5, 50}.

and κ ≈ r(d-r 2 ) 1-r 2 where x i 's are random samples from the vMF distribution and r = n i=1 xi n

Figure 2: The directions of minibatch gradients become more important than their lengths as the batch size increases. (a) The gradient norm stochasticity of ĝ(w) with respect to various batch sizes at 5 random points w with mean (black line) and mean±std.(shaded area) in a log-linear scale; (b) If the gradient norm stochasticity is sufficiently low, then the directions of ĝi (w)'s need to be balanced to satisfy 3 i=1 ĝi (w) ≈ 0.

Figure 3: (a) Asymptotic angle densities (2) of θ(u, v) = 180 π cos -1 u, v where u and v are independent uniformly random unit vectors for each large dimension d. As d → ∞, θ(u, v) tends to less scattered from 90 (in degree). (b-c) We apply SGD on FNN for MNIST classification with the batch size 64 and the fixed learning rate 0.01 starting from five randomly initialized parameters.We draw a density plot θ(u, ĝj (w)) ĝj (w)) ) for 3, 000 minibatch gradients (black) at w = w 0 0 (b) and w = w 0 final , with training accuracy of > 99.9%, (c) when u is given. After SGD iterations, the density of θ(u, ĝj (w)) converges to an asymptotic density (red). The dimension of FNN is 635,200.

If we consider κ = h(u) as a function of u = n b i=1 x i , then h(•) is Lipschitz continuous on [0, n b (1 -)] for any > 0. Moreover, h(•) and h (•) are strictly increasing and increasing on [0, n b ), respectively. Proof. See Supplemental C.1.

Figure 4: (a) If the point is slightly moved from w 0 0 to w by i x i wherex i = (p iw 0 0 )/ p i -w 0 0 and x i = (p i -w )/ p i -w , then i x i < i x i which is equivalent to κ(w ) < κ(w 0 0 ); (b) If ĝi (w00 )'s are sufficiently parallel to (p i -w 0 0 )'s for each i, then w ≈ w 0 0 + η i ĝi (w 0 0 ). When w and w 0 1 are sufficiently close to each other, we also have κ(w 0 1 ) < κ(w 0 0 ).

Figure6: We show the average κ (black curve) ± std. (shaded area), as the function of the number of training epochs(in log-log scale) for all considered setups. We report the name of architecture, dataset and maximum valid accuracy(mean±std.) of all epochs. Although κ decreases eventually over time in all the cases, batch normalization (+BN) significantly reduces the variance of the directional stochasticity ((a-d) vs. (e-h)). We also observe that the skip connections make κ decrease monotonically ((c,g) vs. (d,h)). Note the differences in the y-scales.

Figure7: (First row) We plot the evolution of the training loss (Train loss), validation loss (Valid loss), inverse of gradient stochasticity (SNR), inverse of gradient norm stochasticity (normSNR) and directional uniformity κ. We normalized each quantity by its maximum value over training for easier comparison on a single plot. In all the cases, SNR (orange) and κ (red) are almost entirely correlated with each other, while normSNR is less correlated. (Second row) We further verify this by illustrating SNR-κ scatter plots (red) and SNR-normSNR scatter plots (blue) in log-log scales. These plots suggest that the SNR is largely driven by the directional uniformity.

180 π . Hence, by Lemma B.1 and Theorem B.2 (Delta method), the desired convergence in distribution holds. C PROOFS FOR THEOREM 3 C.1 PROOF OF LEMMA 1Lemma 1. The approximated estimator of κ induced from the d-dimensional unit vectors{x 1 , x 2 , • • • , x n b }, increasing function on [0, 1]. If we consider κ = h(u) as function of u = n b i=1 x i . Then h(•) is Lipschitz continuous on [0, n b (1 -)]for any > 0. Moreover, h(•) and h (•) are strict increasing and increasing on [0, n b ), respectively.Proof. Note thatn b i=1 x i ≤ n b i=1 ||x i || = n b . Therefore, we have r ∈ [0, 1]. If d = 1, then κ = r and this increases on [0, 1]. For d > 1, dκ dr = d + r4 + (d -3r 2 )(1 -r2 ) 2 and its numerator is always positive for d > 2. When d = 2, r2 ) 2 > 0. So κ increases as r increases.The Lipschitz continuity of h(•) directly comes from the continuity of dκ dr since

Figure9: (a,c,e) We plot the evolution of the training loss (Train loss), validation loss (Valid loss), inverse of gradient stochasticity (SNR), inverse of gradient norm stochasticity (normSNR) and directional uniformity κ. We normalized each quantity by its maximum value over training for easier comparison on a single plot. In all the cases, SNR (orange) and κ (red) are almost entirely correlated with each other, while normSNR is less correlated. (b,d,f) We further verify this by illustrating SNRκ scatter plots (red) and SNR-normSNR scatter plots (blue) in log-log scales. These plots suggest that the SNR is largely driven by the directional uniformity.

Figure10: (a,c,e) We plot the evolution of the training loss (Train loss), validation loss (Valid loss), inverse of gradient stochasticity (SNR), inverse of gradient norm stochasticity (normSNR) and directional uniformity κ. We normalized each quantity by its maximum value over training for easier comparison on a single plot. In all the cases, SNR (orange) and κ (red) are almost entirely correlated with each other, while normSNR is less correlated. (b,d,f) We further verify this by illustrating SNRκ scatter plots (red) and SNR-normSNR scatter plots (blue) in log-log scales. These plots suggest that the SNR is largely driven by the directional uniformity.

Figure11: (a,c,e) We plot the evolution of the training loss (Train loss), validation loss (Valid loss), inverse of gradient stochasticity (SNR), inverse of gradient norm stochasticity (normSNR) and directional uniformity κ. We normalized each quantity by its maximum value over training for easier comparison on a single plot. In all the cases, SNR (orange) and κ (red) are almost entirely correlated with each other, while normSNR is less correlated. (b,d,f) We further verify this by illustrating SNRκ scatter plots (red) and SNR-normSNR scatter plots (blue) in log-log scales. These plots suggest that the SNR is largely driven by the directional uniformity.

Figure12: (a,c,e) We plot the evolution of the training loss (Train loss), validation loss (Valid loss), inverse of gradient stochasticity (SNR), inverse of gradient norm stochasticity (normSNR) and directional uniformity κ. We normalized each quantity by its maximum value over training for easier comparison on a single plot. In all the cases, SNR (orange) and κ (red) are almost entirely correlated with each other, while normSNR is less correlated. (b,d,f) We further verify this by illustrating SNRκ scatter plots (red) and SNR-normSNR scatter plots (blue) in log-log scales. These plots suggest that the SNR is largely driven by the directional uniformity.

The value of κ (mean±std.) estimated from {1, 000, 2, 000, 3, 000} random samples of the vMF distribution with κ = 0 across various dimensions.

Table content: <table><tr><td>Network</td><td colspan="2">Dimension κ (1,000 samples)</td><td>κ (2,000 samples)</td><td>κ (3,000 samples)</td></tr><tr><td>CNN</td><td>206, 128</td><td>6, 527.09 ± 5.98</td><td>4, 608.16 ± 6.18</td><td>3, 762.45 ± 7.23</td></tr><tr><td>CNN+Res</td><td>206, 128</td><td>6, 527.09 ± 5.98</td><td>4, 608.16 ± 6.18</td><td>3, 762.45 ± 7.23</td></tr><tr><td>CNN+BN</td><td>207, 152</td><td>6, 563.62 ± 8.68</td><td>4, 633.84 ± 5.60</td><td>3, 781.04 ± 3.12</td></tr><tr><td colspan="2">CNN+Res+BN 207, 152</td><td>6, 563.62 ± 8.68</td><td>4, 633.84 ± 5.60</td><td>3, 781.04 ± 3.12</td></tr><tr><td>FNN</td><td>635, 200</td><td colspan="3">20, 111.90 ± 13.04 14, 196.89 ± 14.91 11, 607.39 ± 9.27</td></tr><tr><td>FNN+BN</td><td>636, 800</td><td colspan="3">20, 157.57 ± 14.06 14, 259.83 ± 16.38 11, 621.63 ± 6.83</td></tr><tr><td>DFNN</td><td colspan="4">1, 915, 200 60, 619.02 ± 13.49 42, 849.86 ± 18.90 34, 983.31 ± 15.62</td></tr><tr><td>DFNN+BN</td><td colspan="4">1, 920, 000 60, 789.84 ± 17.93 42, 958.71 ± 25.61 35, 075.99 ± 12.39</td></tr></table>

Peer review (iclr2019):
 Summary: This work provides an analysis of the directional distribution of of stochastic gradients in SGD. The basic claim is that the distribution, when modeled as a von Mises-Fisher distribution, becomes more uniform as training progresses. There is experimental verification of this claim, and some results suggesting that the SNR is more correlated with their measure of uniformity than with the norm of the gradients. Quality: The proofs appear correct to me. Clarity: The paper is generally easy to read. Originality & Significance: I don*t know of this specific analysis existing in the literature, so in that sense it may be original. Nonetheless, I think there are serious issues with the significance. The idea that there are two phases of optimization is not particularly new (see for example Bertsekas 2015) and the paper*s claim that uniformity of direction increases as SGD convergence is easy to see in a simple example. Consider f_i(x) = |x-b_i|^2 quadratics with different centers. Clearly the minimum will be the centroid. Outside of a ball of certain radius from the centroid all of the gradients grad f_i point in the same direction, closer to the minimum they will point towards their respective centers. It is pretty clear, then that uniformity goes up as convergence proceeds, depending on the arrangement of the centers. The analysis in the paper is clearly more general and meaningful than the toy example, but I am not seeing what the take-home is other than the insight generated by the toy example. The paper would be improved by clarifying how this analysis provides additional insight, providing more analysis on the norm SNR vs uniformity experiment at the end. Pros: - SGD is a central algorithm and further analysis laying out its properties is important - Thorough experiments. Cons: - It is not entirely clear what the contribution is. Specific comments: - The comment at the top of page 4 about the convergence of the minibatch gradients is a bit strange. This could also be seen as the reason that analysis of the convergence of SGD rely on annealed step sizes. Without annealing step-sizes, it*s fairly clear that SGD will converge to a kind of stochastic process. - The paper would be stronger if the authors try to turn this insight into something actionable, either by providing a theoretical result that gives guidance or some practical algorithmic suggestions that exploit it. Dimitri P. Bertsekas. Incremental Gradient, Subgradient, and Proximal Methods for Convex Optimization: A Survey. ArXiv 2015.