Year: 2020

Abstract:
A central ingredient in the impressive predictive performance of deep neural networks is optimization via stochastic gradient descent (SGD). While some theoretical progress has been made, the effect of SGD in neural networks is still unclear, especially during the early phase of training. Here we generalize the theory of thermophoresis from statistical mechanics and show that there exists an effective force from SGD that pushes to reduce the gradient variance in certain parameter subspaces. We study this effect in detail in a simple two-layer model, where the thermophoretic force functions to decreases the weight norm and activation rate of the units. The strength of this effect is proportional to squared learning rate and inverse batch size, and is more effective during the early phase of training when the model’s predictions are poor. Lastly we test our quantitative predictions with experiments on various models and datasets.

text: A central ingredient in the impressive predictive performance of deep neural networks is optimization via stochastic gradient descent (SGD). While some theoretical progress has been made, the effect of SGD in neural networks is still unclear, especially during the early phase of training. Here we generalize the theory of thermophoresis from statistical mechanics and show that there exists an effective force from SGD that pushes to reduce the gradient variance in certain parameter subspaces. We study this effect in detail in a simple two-layer model, where the thermophoretic force functions to decreases the weight norm and activation rate of the units. The strength of this effect is proportional to squared learning rate and inverse batch size, and is more effective during the early phase of training when the model’s predictions are poor. Lastly we test our quantitative predictions with experiments on various models and datasets.

heading: 1 INTRODUCTION
text: Deep neural networks have achieved remarkable success in the past decade on tasks that were out of reach prior to the era of deep learning. Yet fundamental questions remain regarding the strong performance of over-parameterized models and optimization schemes that typically involve only first-order information, such as stochastic gradient descent (SGD) and its variants.
In particular, optimization via SGD is known in many cases to result in models that generalize better than those trained with full-batch optimization. To explain this, much work has focused on how SGD navigates towards so-called flat minima, which tend to generalize better than sharp minima (Hochreiter & Schmidhuber, 1997; Keskar et al., 2017). This has been argued by nonvacuous PACBayes bounds (Dziugaite & Roy, 2017) and Bayesian evidence (Smith & Le, 2018). More recently, Wei & Schwab (2019) discuss how optimization via SGD pushes models to flatter regions within a minimal valley by decreasing the trace of the Hessian.
However, these perspectives apply to models towards the end of training, whereas it is known that proper treatment of hyperparameters during the early phase is vital. In particular, when training a deep network one typically starts with a large learning rate and small batch size if possible. After training has progressed, the learning rate is annealed and decreased so that the model can be further trained to better fit the training set (Krizhevsky et al., 2012; Simonyan & Zisserman, 2015; He et al., 2016b;a; You et al., 2017; Vaswani et al., 2017). Crucially, using a small learning rate during the first phase of training usually leads to poor generalization and also result in large gradient variance practically (Jastrzebski et al., 2020; Faghri et al., 2020).
However, limited theoretical work has been done to understand the effect of SGD on the early phase of training. Jastrzebski et al. (2020) argue for the existence of a “break-even” point on an SGD trajectory. This point depends strongly on the hyperparameter settings. They argue that the breakeven point with large learning rate and small batch size tends to have a smaller leading eigenvalue of the Hessian spectrum, and this eigenvalue sets an upper bound for the leading eigenvalue beyond this point. They also present experiments showing that large learning rate SGD will reduce the variance of the gradient. However their analysis focuses only on the leading eigenvalue of the Hessian spectrum and requires the strong assumption that the loss function in the leading eigensubspace is quadratic.
Meanwhile Li et al. (2020) studied the simple setting of two-layer neural networks. They demonstrate that in this model, training with large learning rate in the early phase tends to result in better generalization than training with small learning rate. To explain this, they hypothesize a separation of features in the data: easy-to-generalize yet hard-to-fit features, and hard-to-generalize, easierto-fit features. They argue that a model trained with small learning rate will memorize easy-togeneralize, hard-to-fit patterns during phase one, and then generalize worse on hard-to-generalize, easier-to-fit patterns, while the opposite scenario occurs when training with large learning rate. However, this work relies heavily on the existence of these two distinct types of features in the data and the specific network architecture. Moreover, their analysis focuses mainly on learning rate instead of the effect of SGD.
In this paper, we study the dynamics of model parameter motion during SGD training by borrowing and generalizing the theory of thermophoresis from physics. With this framework, we show that during SGD optimization, especially during the early phase of training, the activation rate of hidden nodes is reduced as is the growth of parameter weight norm. This effect is proportional to squared learning rate and inverse batch size. Thus, thermophoresis in deep learning acts as an implicit regularization that may improve the model’s ability to generalize.
We first give a brief overview of the theory of thermophoresis in physics in the next section. Then we generalize this theory to models beyond physics and derive particle mass flow dynamics microscopically, demonstrating the existence of thermophoresis and its relation to relevant hyperparameters. Then we focus on a simple two-layer model to study the effect of thermophoresis in detail. Notably, we find the thermophoretic force is strongest during the early phase of training. Finally, we test our theoretical predictions with a number of experiments, finding strong agreement with the theory.

heading: 2 THERMOPHORESIS IN PHYSICS
text: Thermophoresis, also known as the Soret effect, describes particle mass flow in response to both diffusion and a temperature gradient. The effect was first discovered in electrolyte solutions (Ludwig, 1859; Soret, 1897; Chipman, 1926). However it was discovered in other systems such as gases, colloids, and biological fluids and solid (Janek et al., 2002; Köhler & Morozov, 2016).
Thermophoresis typically refers to particle diffusion in a continuum with a temperature gradient. In one method of analysis, the non-uniform steady-state density ρ is given by the ”Soret Equilibrium” (Eastman, 1926; Tyrell & Colledge, 1954; Wurger, 2014),
∇ρ+ ρST∇T = 0 , (1)
where T is temperature and ST is called the Soret coefficient.
In other work by de Groot & Mazur (1962), mass flow was calculated by non-equilibrium theory. They considered two types of processes for entropy balance: a reversible process stands for the entropy transfer and an irreversible process corresponds to the entropy production, or dissipation. The resulting mass flow induced by diffusion and temperature gradient was found to be
J = −D∇ρ− ρDT∇T , (2)
where D is the Einstein diffusion coefficient and DT is defined as thermal diffusion coefficient. Comparing the steady state in 1 and setting the flow to be zero, the Soret coefficient is simply
ST = DT D . (3)
The Soret coefficient can be calculated from molecular interaction potentials based on specific molecular models (Wurger, 2014).

heading: 3 THERMOPHORESIS IN GENERAL
text: In this section, we first study a kind of generalized random walk that has evolution equations for a particle state with coordinate q = {qi}i=1,...,n as
qt+1 = qt − ηγf(qt, ξ) , (4)
where f is a vector function, γ and ξ are random variables, and η is a small number controlling the step size. Notice that this is a generalized inhomogeneous random walk for the particle. Before further analysis, it is noted that the evolution equations 4 is similar to SGD updates in machine learning and we will show this in the next section.
To isolate the effect of thermophoresis, we assume the random walk is unbiased, in which case
P (γf(q, ξ) = a) = P (γf(q, ξ) = −a), (5)
for an arbitrary vector a. Thus there is no explicit force exerted on the particle. This simplification was used to demonstrate a residual thermophoretic force in the absence of a gradient. Including gradients is straightforward and corresponds to an external field that creates a bias term. We also denote the probability density, which we also call the mass density, as ρ(q) and
gi(q) := √∫ γ2f2i (q, ξ)dµ(γ, ξ), (6)
so that ηgi(q) is the standard deviation of the random walk in the ith direction.
From a position q, we consider a subset of coordinate indices, U ⊆ {1, . . . , n}, wherein
sign(fi(q, x)) = sign(fj(q, x)) and ∂igj(q) ≥ 0 (7)
for all i, j ∈ U . We note here that indices will correspond to parameters when we study learning dynamics. The first property is necessary for our derivation. The second condition will be used at the end to conclude that each gi decreases.
In order to study the dynamics of the particle and its density function, we focus on the probability mass flow induced by the inhomogeneous random walk. We will show that there is always a flow from regions with larger gi(q) to those with smaller gi(q) for i ∈ U , which is a generalization of thermophoresis in physics.
Since η 1, the movement of the particle will have a mean free path of gi(q) in ith direction. Therefore the random walk equation 4 becomes
qi = qi − ηgi(q)ζi, (8)
where i = 1, . . . , n and ζi is a binary random variable with P (ζi = −1) = P (ζi = 1) = 0.5. Moreover, from Eq. 7, we also have that ζi = ζj for all i and j ∈ U . Next we will show that the flow projecting on the subspace U is always toward smaller gi(q). Notice that although U can be multi-dimensional, the degree of freedom of the particle dynamics is 1 within U due to the sharing of the ζs, and therefore the mass flow projecting on it is also 1-dimensional. For each i ∈ U , we define the average flow in this dimension to be the mass that enters qi from q−i minus the mass from the opposite direction q+i . From Eq. 8 and the assumption that η 1, only mass close to qi will move across qi at each step. We let the farthest mass that will flow across qi in step i be qi + ∆+i and qi − ∆ − i , where ∆ + i and ∆ + i are positive. ∆ + i and ∆ − i are thus defined implicitly by the equations: ∆+i = ηgi(q+∆ +) and ∆−i = ηgi(q−∆−), respectively. Notice that if the random walk were homogeneous, we would have ∆+i = ∆ − i . In our inhomogeneous case, we have ∆+i ∼ ∆ − i ∼ ηgi(q) up to leading order of η, and the next to leading order will be calculated in order to compute the difference between ∆+i and ∆ − i .
Now we are ready to calculate the mass flow through q. The mass flow projecting onto the subspace U is calculated by the mass through q from q + ∆+ minus the mass from q−∆− where ∆+i and ∆−i are as above for i ∈ U . It is straightforward to show that1
∆+i −∆ − i = 2η 2 ∑ j∈U gj(q)∂jgi(q) +O(η 3). (9)
With this, we can compute the flow density, J , through q, finding J = −η2 √∑ i∈U g2i (q) ∑ i∈U gi(q)∂iρ(q)− η2 ∑ i,j∈U gi(q)gj(q)∂jgi(q)√∑
i∈U g 2 i (q)
ρ(q) +O(η3), (10)
where the derivation can be found in Appendix A.3. This can be understood as described in Diagram 1. Notice that this probability mass flow consists of two terms at order η2. The first represents diffu-
sion and the second corresponds to our goal in this section, namely the flow due to thermophoresis. By the second property of the gi in Eq. 7, we find that the coefficient of thermophoresis (Soret coefficient), which is defined as
c := −η2 ∑ i,j∈U gi(q)gj(q)∂jgi(q)
2 √∑
i∈U g 2 i (q)
(11)
≤ 0, (12)
is negative. This means that there is an effective force exerted on a particle at position q towards the smaller variance regime (by analogy, the colder area). The coefficient is proportional to η2.

heading: 4 MATHEMATICAL MODEL

heading: 4.1 TWO-LAYER MODEL
text: To study the physics behind SGD optimization in detail, we consider the simple setting of onehidden layer neural networks. The network is a function f : RM → R parameterized as follows:
f(x; V,W,b) = Vσ(Wx + b),
= N∑ i=1 Viσ( M∑ j=1 Wijxj + bi).
We also write f(x) for simplicity. The network has a scalar output, which is widely used in regression and binary classification. x is the network input with dimension M , W and b are the weights and biases in the first layer with dimensionN×M andN respectively, whereN is the number of hidden nodes in the hidden layer, and σ is the ReLU activation function defined as σ(a) = max(0, a).
The dataset is drawn i.i.d. from the data distribution, {(x, y)|(x, y) ∼ D(x, y)}. In this paper we consider two cases, where either xi ≥ 02 or xi ∼ N (0, 1)3. Here y ∈ Y and we denote the marginal distribution of y as DY . Finally, we have the loss function L : R× Y→ R+.

heading: 4.2 TRAINING
text: We consider optimization via SGD, where the gradient of the loss on a batch of size |B| is given by 1A brief derivation can be found in Appendix A.2. 2Usually in convolutional neural networks or intermediate layers. 3Often found when the data are normalized.
∇LB(V,W,b) = 1
|B| |B|∑ i=1 ∇fL(f(xi), yi)∇f(xi). (13)
In our two-layer model, we have
∇Vif(x) = σ( M∑ j=1 Wijxj + bi), ∇Wijf(x) = Vixjσ′( M∑ k=1 Wikxk + bi),
∇bif(x) = Viσ′( M∑ k=1 Wikxk + bi).
For an input vector x, we call the hidden node i activated when σ′( ∑M k=1Wikxk + bi) = 1, or equivalently Wikxk + bi > 0. We thus define the activation rate of the network to be
σ′ = 1
N N∑ i=1 Exσ′( M∑ k=1 Wikxk + bi) . (14)
This is an important concept to which we will return.
Henceforth, we drop the index i, since the dynamical equations are invariant with respect to node index, and write V := Vi, Wj := Wij and b := bi by abuse of notation. We also denote
hv(V,W, b) := Ex[∇V f(x)]2, hw(V,W, b) := Ex[∇Wf(x)]2, hb(V,W, b) := Ex[∇bf(x)]2,
where Ex denotes average over input x. We have the following property for the functions h: Property 4.1. Given W, if V 21 ≤ V 22 and b1 ≤ b2, we have
hv(V1,W, b1) ≤ hv(V2,W, b2), hw(V1,W, b1) ≤ hw(V2,W, b2), hb(V1,W, b1) ≤ hb(V2,W, b2), σ′(V1,W, b1) ≤ σ′(V2,W, b2).
Here we define a ≤ b as min(b− a) ≥ 0.
It is straightforward to see the following: Property 4.2. When the case of xi ≥ 0 is considered, if V 21 ≤ V 22 , W1 ≤ W2 and b1 ≤ b2, we have
hv(V1,W1, b1) ≤ hv(V2,W2, b2), hw(V1,W2, b1) ≤ hw(V2,W, b2). hb(V1,W2, b1) ≤ hb(V2,W, b2), σ′(V1,W, b1) ≤ σ′(V2,W, b2).
In our analysis, we focus for simplicity on binary classification tasks, where the loss is typically binary cross-entropy: L(f, y) = y ln p(f) + (1− y) ln(1− p(f)) and p(f) = 1/(1 + exp(f)). We thus have ∇fL(f, y) = p(f)− y. (15) Substituting into Eq. 13, the mini-batch gradient becomes
∇LB(V,W,b) = 1
|B| |B|∑ i=1 (pi − yi)∇f(xi). (16)
Our results also hold straightforwardly for squared error.

heading: 5 THERMOPHORESIS IN DEEP LEARNING
text: In this section, we will show that the parameters in the one hidden layer model and their dynamics approximately satisfy the criteria of the previous section and that the biases are pushed negative and V 2 is suppressed during training, the effects of both of which are proportional to squared learning rate η2 and inverse batch size 1/|B|. The gradient that dominates model training is defined in 16. Because training samples are i.i.d., the variance of the gradient is
var [ ∇LB(V,W,b) ] = var
[ 1
|B| |B|∑ i=1 (pi − yi)∇f(xi) ] , (17)
= 1 |B| var [ (p− y)∇f(x) ] (18)
The gradient has two components: p−y corresponding to γ in equation 4 and∇f(x) corresponding to f(q, ξ). We assume that the dataset is unbiased, in which case P (y = 0) = P (y = 1) = 0.5 and P (p − y = a) = P (p − y = −a), and that p − y and ∇f(x) are independent in the first period of training given that the dataset is complex and can’t be learned by linear model. It is straightforward to see that it satisfies Eq. 5.
Next we will show that V and b are always in the set of U , i.e. they satisfy the conditions of Eq. 7. First, if Vi ≥ 0, we have
∇Vif(x) = σ( M∑ j=1 Wijxj + bi), (19)
≥ 0. (20) and
∇bif(x) = Viσ′( M∑ k=1 Wikxk + bi), (21)
≥ 0. (22) Since we also have Property 4.1, the conditions in Eq. 7 are satisfied. If Vi < 0, we consider a coordinate transform that maps Vi to V̄i = −Vi. It is easy to show that Eq. 7 is again satisfied after this transform.
Next we consider W. The gradient of f with respect to Wij is the product of ∇bif and xi. If xi for i = 1, . . . ,M are always ≥ 0, which is usually the case in convolutional neural networks, it is easy to show that Wij is also in set U and smaller Wij corresponds to smaller variance according to Property 4.2. If xi ∼ N (0, 1), on the other hand, W is excluded from U . For the following, we only consider the case where xi ∼ N (0, 1), and
gV (Vi,Wi, bi) = 1√ |B| √√√√∫ [(p− y)σ( M∑ j=1 Wijxj + bi) ]2 dµ(x, y) , (23)
:= 1√ |B| φ1(Wi, bi) , (24)
gb(Vi,Wi, bi) = 1√ |B| √√√√∫ [(p− y)Viσ′( M∑ j=1 Wijxj + bi) ]2 dµ(x, y) , (25)
:= Vi√ |B| φ2(Wi, bi) , (26)
where g is defined as in Eq. 6. Inserting these into Eq. 10, we find the thermophoresis flow density to be
Jt = η2
|B| ψ , (27)
whereψ = Viφ1φ 2 2+Viφ1φ2∂bφ1+V 3φ22∂bφ2
2 √ φ21+V 2 i φ 2 2
ρ. This flow biases the model toward smaller bi and smaller
Vi 4 with a strength proportional to squared learning rate η2 and inverse batch size. We also note that ψ can be bounded by a function multiplying with a scalar ∫
(p−y)2µ(x, y). It is clear that this scalar measures the L-2 distance between model predictions and sample labels and decreases on average during training as prediction getting better. Thus thermophoresis is more effective during the early phase of training.
Therefore there exists an effective force that pushes to decrease the model’s activation rate, defined in equation 14, and reduces the weight norm of the second layer. The strength of this force scales as
F ∝ η 2
|B| . (28)
In Li & Liang (2018), Theorem 4.1 presents a linear relation between learning rate and training iterations for a target training error and small learning rate. This implies that if one uses a learning rate k times larger, the model will require k times fewer optimization steps for the same training performance. Together with our results, this implies the following: for the same model and initialization, comparing two optimization schemes with η1 ≤ η2 each achieving a given training error, the activation rate for scheme 1 will be at least as large as that for scheme 2, i.e. σ1 ≥ σ2. Similarly, denoting the weight norm for scheme 1(2) by v1(v2), we have that v1 ≥ v2. Model sparsity can mean two different things: sparsity of the weights, and frequency with which units are activated, called the activation rate. Intuitively, a sparser model has a smaller capacity Bizopoulos & Koutsouris (2020); Kurtz et al. (2020); Aghasi et al. (2017); Lin et al. (2017). Furthermore, certain forms of model pruning have been shown to improve generalization Frankle & Carbin (2019); Frankle et al. (2020). Therefore one might surmise that a smaller activation rate in general correlates with generalization. Moreover, in Appendix A.5, we construct an upper bound of the Hessian norm and this upper bound monotonously depends on activation rate and weight norm. This also sheds light on the connection between sparsity, weight norm, and generalization.
Our theory can also be generalized beyond two-layer models. We have shown that there exists an effective force in deep neural networks from SGD that reduces the gradient variance and have derived quantitative properties of it.

heading: 6 EXPERIMENTS
text: The essential result from the previous section is that there exists an effective force from SGD, analogous to thermophoresis, that pushes to decrease the gradient variance, and in one-hidden-layer neural networks decreases the model’s activation rate and reduces the weight norm of the second layer. The strength of the force is proportional to squared learning rate and inverse batch size. In this section, we present experiments to test these results. Further experiments can be found in the appendix.
First we consider a one hidden layer model with input dimension 100 and 100 hidden units. The input data, x, is distributed as N (0, I) where I is identity matrix, and the label is randomly chosen from {0, 1}. Batch size is set to 1 and the learning rate is varied from 0.025 to 0.1. We calculate the activation rate and L2 norm of the vector V after each training iteration. The result for activation rate is shown in the first row of Fig. 2. The leftmost plot shows activation rate as a function of true iteration on the x-axis, and we see that activation rate decreases during training, and the decreasing is more rapid with larger learning rate. In the middle plot we rescale the x-axis by a factor proportional to learning rate η5. This rescaling factor is to offset the movement difference due to learning rate difference. It is clear that even after this rescaling, we still observe that larger learning rates decrease activation rate faster. Finally, on the rightmost plot we rescale the x-axis with a factor proportional to squared learning rate η2. We see that all trajectories now overlap, which matches our prediction in the previous section that decreasing rate is proportional to η2.
4larger Vi if Vi < 0. 5For example, if raw iteration number for η = 0.05 is 1000 and rescaled iteration number is also 1000, the
rescaled iteration number for η = 0.1 is 1000 then its true iteration number is 500.
We next test our results for deep neural networks beyond the two-layer model. Instead of activation rate and weight norm, we plot the gradient variance as predicted by our theory. Networks architectures are 6-layer fully-connected with hidden layer sizes of 100 and 6-layer convolutional with 10 channels with kernel size of 5*5 and stride 1 except the last fully-connected layer output. The results are shown in the second row of Fig. 2 and the third row of Fig. 2 respectively.
Further experiments can be found in Appendix A.6, where we show that the evolution of the weight norm, scaling with batch size, and other results are consistent with our theoretical predictions. We also study other models and other datasets including CIFAR10.

heading: 7 CONCLUSION
text: In this paper we generalized the theory of thermophoresis from statistical mechanics and showed that there exists an effective thermophoretic force from SGD that pushes to reduce the gradient variance. We studied this effect in detail for a simple two-layer model, where the thermophoretic force serves to decrease the weight norm and the activation rate of the units. We found that the strength of this effect is proportional to square of the learning rate, inversely proportional to batch size, and is more effective during the early phase of training when the model’s predictions are poor. We found good agreement between our predictions and experiments on various models and datasets.

heading: A APPENDIX

heading: A.1 PROOF OF PROPERTY 4.1
text: Proof. By definition, we have
σ( M∑ j=1 Wjxj + b1) ≤ σ( M∑ j=1 Wjxj + b2) , (29)
σ′( M∑ j=1 Wjxj + b1) ≤ σ′( M∑ j=1 Wjxj + b2) . (30)
For hv ,
hv(V1,W, b1) = Exg2v(x, V1,W, b1) ,
= Exσ2( M∑ j=1 Wjxj + b1) , ≤ Exσ2( M∑ j=1 Wjxj + b2) ,
= hv(V2,W, b2) .
Similarly, we have
hwi(V1,W, b1) = ExV 21 x2iσ′( M∑ k=1 Wkxk + b1) ,
≤ ExV 22 x2iσ′( M∑ k=1 Wkxk + b1) ,
≤ ExV 22 x2iσ′( M∑ k=1 Wkxk + b2) ,
= hwi(V2,W, b2) .
Clearly the inequality also holds for hb.
A.2 DERIVATION OF EQ. 9
∆+i −∆ − i = ηgi(q + ∆ +)− ηgi(q−∆−) , (31)
= η n∑ j∈U (∆+j + ∆ − j )∂jgi(q) +O(η∆ 2) , (32)
= 2η2 ∑ j∈U gj(q)∂jgi(q) +O(η 3) . (33)

heading: A.3 DERIVATION OF EQ. 10
text: J =− 1 2 |∆+|ρ(q + ∆+) + 1 2 |∆−|ρ(q−∆−) ,
= 1 2 |∆+|[ρ(q−∆−)− ρ(q + ∆+)] + 1 2 (|∆−| − |∆+|)ρ(q−∆−) ,
=− 1 2 |∆+|(|∆+ + ∆−|)ρ(q + ∆ +)− ρ(q−∆−) |∆+ + ∆−|
− 1 2 (|∆+|2 − |∆−|2) |∆+|+ |∆−| ρ(q−∆−) ,
≈− 1 2 |∆+|(∆+ + ∆−)∇ρ(q)− 1 2
∑ i∈U (∆ + i + ∆ − i )(∆ + i −∆ − i )
|∆+|+ |∆−| ρ(q−∆−) ,
=− η2 √∑ i∈U g2i (q) ∑ i∈U gi(q)∂iρ(q)− η2 ∑ i,j∈U gi(q)gj(q)∂jgi(q)√∑
i∈U g 2 i (q)
ρ(q) +O(η3) .

heading: A.4 SANITY CHECK OF GENERALIZED THEORY
text: If |U | = 1 and gi(q) = gi(qi), the model will reduce to aforementioned physics model and the Soret coefficient reduces to
c = η2
2 g(q)g′(q) , (34)
= [( ηg(q)
2 )2]′ , (35)
≈ ∇T , (36)
where T is the effective temperature in the model. This result is consistent with thermophoresis model in physics.

heading: A.5 SPARSITY, WEIGHT NORM AND THEIR RELATION TO GENERALIZATION
text: In this section, we demonstrate how sparsity is related to the Hessian norm. We first denote the model’s probabilistic prediction on a C-class classification as
pµk = exp zµk∑C l=1 exp z µ l , (37)
where k is the probability for label k, µ is the data index, z is model output and C is the total number of categories. We consider cross entropy loss of the form
L(w) = − 1 B B∑ µ=1 C∑ k=1 yµk log p µ k , (38)
where y is sample labels and p stands for model probability prediction, similar to the previous definition. We denote the loss for individual sample to be
Lµ = − C∑ k=1 yµk log p µ k . (39)
The gradient with respect to the model output is
(∇zLµ)k = −yµk + p µ k . (40)
And it is easy to show that the Hessian with respect to output is
(∇2zLµ)kl = δklp µ k − p µ kp µ l . (41)
Therefore the Hessian with respect to model parameters is
Hµ = ∇2wL(z(w)) , (42) = ∇w(∇zL ∗ ∇wz) , (43) = (∇wz)(∇2zL)(∇wz) +∇zL∇2wz , (44) ≈ (∇wzµ)ij(∇2zLµ)jk(∇wzµ)kl . (45)
To study the spectrum of the Hessian, we calculate the trace and have
Tr(Hµ) ≈ Tr((∇wzµ)(∇2zLµ)(∇wzµ)T ) , (46) = Tr((∇2zLµ)(∇wzµ)T (∇wzµ)) , (47) = Tr(P ∗K) , (48) ≤ Tr(P ) ∗ Tr(K) , (49)
where
P = ∇2zLµ , (50) Kµν = ∑ l ∑ ij ( ∂zµ ∂wlij )( ∂zν ∂wlij ) . (51)
The trace of K therefore can be calculated by chain rule, Tr(K) = ∑ µ ∑ l ∑ ij ( ∂zµ ∂wlij )2 , (52)
= ∑ l ( ∑ iu (δli[µ]) 2 ∑ j (hl−1j ) 2) , (53)
where δ and h carry backward and forward information respectively. They are defined as
δli[µ] = δµnLW L nLnL−1σ ′...W l+1nl+2iσ ′ , (54) hl−1j = σW l−1 jnl σ...Wn1n0xn0 . (55)
It can further be shown that∑ iu (δli[µ]) 2 = ∑ iu δµnLW L nLnL−1σ ′ . . .W l+1nl+2iσ ′ ∗ σ′W l+1inl+2 . . . σ′WLnL−1nLδnLµ , (56)
= Tr(WLσ′ . . .W l+1σ′σ′W l+1 . . . σ′W̄L) , (57)
≤ Tr(σ′WLWLσ′)Tr(σ′WL−1WL−1σ′) . . . T r(σ′W l+1W l+1σ′) , (58)
as well as ∑ j (hl−1j ) 2 ≤ ‖X‖2Tr(W 1σσW 1)Tr(W 2σσW 2)...T r(W l−1σσW l−1) . (59)
Together with the previous calculations and the definition of K, we have Tr(K) ≤ ‖X‖2 ∑ l Πl−1n=1Tr(W nσσWn)ΠLn=l+1Tr(σ ′WnWnσ′) , (60)
= ‖X‖2 ∑ l Πl−1n=1‖σWn‖2FΠLn=l+1‖Wnσ′‖2F . (61)
Finally, we derive an upper bound for the trace of the Hessian, Tr(Hµ) ≤ Tr(P )‖X‖2 ∑ l Πl−1n=1‖σWn‖2FΠLn=l+1‖Wnσ′‖2F . (62)
Notice that activation rate and weight norm control the magnitude of ‖σWn‖2F and ‖Wnσ′‖2F . Therefore smaller activation rate and weight norm lead to tiger upper bound of the Hessian trace and thus indicate smaller matrix norm. This analysis connects sparsity with Hessian norm, Hessian trace specifically.

heading: A.6 MORE EXPERIMENTS
text: In this section, we design extensive experiments to test the obtained results. First of all, we consider binary classification with BCE loss. The model has the form
f(x) = Vσ(Wx + b) . (63)
Input dimension is 10 and the number of hidden nodes is also 10. x is distributed as N (0, I) where I is identity matrix and the label is randomly chosen from {0, 1}. For the first hyperparameter setting, batch size is 1 and learning rate varies from 0.025 to 0.1. We calculate the activation rate and L2 norm of the vector V after each training iteration. The result for activation rate and weight norm are shown in Fig. 3 and Fig. 4 respectively. Both figures contain three plots. The leftmost plots correspond to plot with raw iteration as x-axis. It is shown that both activation rate and weight norm are decreasing for all cases. Additionally, the decreasing rate is larger with larger learning rate. The middle plots correspond to plot with rescaled iteration as x-axis, where the rescaled factor is proportional to learning rate η6. This rescalation factor is to offset the movement difference due to learning rate difference. It is clear that even after this rescalation, we still observe activation rate and weight norm difference for different learning rates. Lastly, the rightmost plots correspond to plot with rescaled iteration as x-axis, where the rescaled factor is proportional to squared learning rate η2. The result that all trajectories overlap with each other matches our prediction in the previous section, that decreasing rate is proportional to η2.
For the second hyperparameter setting, learning rate is fixed to be 0.05 and batch size varies from 1 to 3. Again we calculate the activation rate and L2 norm of the vector V after each training
6For example, if raw iteration number for η = 0.05 is 1000 and rescaled iteration number is also 1000, the rescaled iteration number for η = 0.1 is 1000 then its true iteration number is 500.
iteration. The result for activation rate and weight norm are shown in Fig. 5 and Fig. 6 respectively. We observe similar tendency as we discussed in the previous hyperparameter setting. The result shows that both activation rate and weight norm are decreasing for all cases. While there exist decreasing rate difference in the left plots due to batch size discrepancy. This difference can be offset be rescaling x-axis according to proportional factor 1/|B|, which is the result in the right plots and it is consistent with our theory prediction.
Subsequently we consider that second case discussed in the previous section. This is also a consider binary classification with BCE loss. The model, however, has the form
f(x) = Vσ(Wx) . (64)
Input dimension is 10 and the number of hidden nodes is also 10. x is distributed uniformly between 0 and 1 as U(0, 1) and the label is randomly chosen from {0, 1}. The first setting is the same as the first setting in the aforementioned experiment. The result for activation rate and weight norm are shown in Fig. 7 and Fig. 8 respectively. Similar to previous experiment, both figures contain three plots. The leftmost, middle, right plots correspond to plot with raw iteration, rescaled iteration with factor η and rescaled iteration with factor η2 respectively as xaxis. The result that all trajectories overlap with each other matches our prediction in the previous section, that decreasing rate is proportional to η2.
The second setting is also similar to the second setting in the previous experiment where we fix learning rate to be 0.05 and vary batch size from 1 to 3. The result for activation rate and weight norm are shown in Fig. 9 and Fig. 10 respectively. The result shows that both activation rate and weight norm are decreasing for all cases. While there exist decreasing rate difference in the left plots due to batch size discrepancy. This difference can be offset be rescaling x-axis according to
proportional factor 1/|B|, which is the result in the right plots and it is consistent with our theory prediction.
Furthermore, we use real image dataset CIFAR10 for the next experiment instead of artificial data. The model in this experiment has one hidden layer with 300 hidden nodes. We first fix batch size to be 1000 and vary learning rate η. The result of activation rate and weight norm decreasing are shown in Fig. 11 and Fig. 12 respectively. We then fix learning rate to be 0.02 and vary batch size. The result of activation rate and weight norm decreasing are shown in Fig. 13 and Fig. 14 respectively.
It is clear that the result matches our theoretical predictions. We skip detailed analysis here as it is similar to our discussion in previous experiments.

Review summaries:
This paper proposes that SGD has the implicitly bias of reducing gradient variance via the phenomenon of thermophoresis that masses tend to flow from regions with higher temperature / variance of random walk to regions with lower temperature / variance of random walk. In the setup of two-layer neural networks trained by SGD for binary classification, the authors show the analogous phenomenon that the model is biased towards smaller activation rate and the norm of its second layer weight. The dependence of the rate at which this phenomenon happens on the learning rate and batch size are verified in the experiments.