On Linear Stability of SGD and Input-Smoothness of Neural Networks
Year: 2022

Abstract:
The multiplicative structure of parameters and input data in the first layer of neural networks is explored to build connection between the landscape of the loss function with respect to parameters and the landscape of the model function with respect to input data. By this connection, it is shown that flat minima regularize the gradient of the model function, which explains the good generalization performance of flat minima. Then, we go beyond the flatness and consider high-order moments of the gradient noise, and show that Stochastic Gradient Descent (SGD) tends to impose constraints on these moments by a linear stability analysis of SGD around global minima. Together with the multiplicative structure, we identify the Sobolev regularization effect of SGD, i.e. SGD regularizes the Sobolev seminorms of the model function with respect to the input data. Finally, bounds for generalization error and adversarial robustness are provided for solutions found by SGD under assumptions of the data distribution.

text: The multiplicative structure of parameters and input data in the first layer of neural networks is explored to build connection between the landscape of the loss function with respect to parameters and the landscape of the model function with respect to input data. By this connection, it is shown that flat minima regularize the gradient of the model function, which explains the good generalization performance of flat minima. Then, we go beyond the flatness and consider high-order moments of the gradient noise, and show that Stochastic Gradient Descent (SGD) tends to impose constraints on these moments by a linear stability analysis of SGD around global minima. Together with the multiplicative structure, we identify the Sobolev regularization effect of SGD, i.e. SGD regularizes the Sobolev seminorms of the model function with respect to the input data. Finally, bounds for generalization error and adversarial robustness are provided for solutions found by SGD under assumptions of the data distribution.

heading: 1 Introduction
text: Stochastic gradient descent (SGD) is the most widely used optimization algorithm to train neural networks [25, 4]. By taking mini-batches of training data instead of all the data in each iteration, it was firstly designed as a substitute of the gradient descent (GD) algorithm to reduce its computational cost. Extensive researches are conducted on the convergence of SGD, both on convex and non-convex objective functions [21, 22, 5, 3]. In these studies, convergence is usually proven in the cases where the learning rate is sufficiently small, hence the gradient noise is small. In practice, however, SGD is preferred over GD not only for the low computational cost, but also for the implicit regularization effect that produces solutions with good generalization performance [7, 16]. Since a trajectory of SGD tends to that of GD when the learning rate goes to 0, this implicit regularization effect must come from the gradient noise induced by mini-batch and a moderate learning rate.
When studying the gradient noise of SGD, a majority of work treat SGD as an SDE, and studies the noise of the SDE [15, 20, 17, 18, 32]. However, SGD is close to SDE only when the learning rate is small, and it is unclear that in practical setting whether the Gaussian noise of SDE can fully characterize the gradient noise of SGD. Some work resort to heavy-tailed noise like the Levy process [27, 37]. Another perspective to study the behavior of SGD is by its linear stability [30, 12]. This is relevant when the learning rate is not very small. The linear stability theory can explain the fast escape of SGD from sharp minima. The escape time derived from this theory depends on the logarithm of the barrier height, while the escape time derived from the diffusion theory based on SDE depends exponentially on the barrier height [32]. The former is more consistent with empirical observations [30].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
An important observation that connects the generalization performance of the solution with the landscape of the loss function is that flat minima tend to generalize better [13, 31]. SGD is shown to pick flat minima, especially when the learning rate is big and the batch size is small [16, 15, 30]. Algorithms that prefer flat minima are designed to improve generalization [6, 14, 19, 35]. On the other hand, though, the reason why flat minima generalize better is still unclear. Intuitive explanations from description length or Bayesian perspective are provided in [13] and [23]. In [10], the authors show sharp minimum can also generalize by rescaling the parameters at a flat minimum. Hence, flatness is not a necessary condition of good generalization performance, but it is still possible to be a sufficient condition. In the study of linear stability in [30], besides the sharpness (a quantity inversely proportional to flatness), another quantity named non-uniformity is proposed which roughly characterizes the second order moment of the gradient noise. It is shown that SGD selects solutions with both low sharpness and low non-uniformity.
In this paper, we build a complete theoretical pipeline to analyze the implicit regularization effect and generalization performance of the solution found by SGD. Our starting points are the following two questions: (1) Why SGD finds flat minima? (2) Why flat minima generalize better? Our answers to these two questions go beyond the flatness and cover the non-uniformity and higher-order moments of the gradient noise. This distinguishes SGD from GD, and is out of the scope that can be explained by SDE. For the first question, we extend the linear stability theory of SGD from the second-order moments of the iterator of the linearized dynamics to the high-order moments. At the interpolation solutions found by SGD, by the linear stability theory, we derive a set of accurate upper bounds of the gradients’ moment. For the second question, using the multiplicative structure of the input layer of neural networks, we show that the upper bounds obtained in the first step regularize the Sobolev seminorms of the model function. Finally, bridging the two components, our main result is a bound of generalization error under some assumptions of the data distribution. The bound works well when the distribution is supported on a low-dimensional manifold (or a union of low-dimensional manifolds). An informal statement of our main result is
(Main result) Around an interpolation solution of the neural network model, assume that (1)the k-th
order moment of SGD’s iterator of the linearized dynamics is stable, (2)with probability at least
1 " the testing data is close to a training data with distance smaller than , and (3) both the model function and the target function are upper bounded by a constant M . Then, at this interpolation
solution we have generalization error . n 1k 2 + "M2, where n is the number of data. The bound depends on the learning rate and the batch size of SGD as
constant factors.
The formal description is stated in Theorem 6 of Section 5. Our analysis also provide bounds for adversarial robustness. As a byproduct, we theoretically show that flatness of the minimum controls the gradient norm of the model function at the training data. Therefore, searching for flat minima has the effect of Lipschitz regularization, which is shown to be able to improve generalization [24].
Lying at the center of our analysis is the multiplicative structure of neural networks, i.e. in each layer the output features from the previous layer is multiplied with a parameter matrix. This structure is a rich source of implicit regularization. In this paper, we focus on the input layer, and build connection between the model function’s derivative with respect to the parameters and its derivatives with respect to the data. Concretely, let f(x,W ) be a neural network model, with x being the input data and W being the parameters. We split the parameters by W = (W1,W2), where W1 is the parameter matrix of the first layer, and W2 denotes all other parameters. Then, the neural network can be represented by the form f(x,W ) = f̃(W1x,W2) due to the multiplicative structure of W1 and x. Accordingly rW1 f̃(W1x,W2) and rxf̃(W1x,W2) take the form
rW1 f̃(W1x,W2) = @f(W1x,W2)
@(W1x) xT , rxf̃(W1x,W2) = WT1
@f(W1x,W2)
@(W1x) , (1)
where @f(W1x,W2)@W1x is usually a long expression produced by back propagation from the output to the first layer. By (1) we have
krxf̃(W1x,W2)k2  kW1k2 kxk2 krW1 f̃(W1x,W2)k2. (2)
Hence, rxf̃(W1x,W2) is upper bounded by rW1 f̃(W1x,W2), given that W1 is not too big and x is not too small. By the bound (2) we can derive the regularization effect of flatness at interpolation
solutions. To see this, let {(xi, yi)}ni=1 be the training data, and L̂(W ) := 12n nP
i=1 (f(xi,W ) yi)2
be the empirical risk given by the square loss. Let W ⇤ = (W ⇤1 ,W ⇤2 ) be an interpolation solution, i.e. f(xi,W ⇤) = yi for any i = 1, 2, ..., n, and L̂(W ⇤) = 0. Define the flatness of the minimum as the sum of the eigenvalues of r2L̂(W ⇤), i.e. flatness(W ⇤) = Tr(r2L̂(W ⇤)). W ⇤ being the interpolation solution implies
flatness(W ⇤) = 1
n
nX
i=1
krW f(xi,W ⇤)k2. (3)
Hence, from (2) we can obtain the following equation that gives bounds for the gradients of the model function with respect to the input data at the training data,
1
n
nX
i=1
krxf(xi,W ⇤)k2  kW ⇤1 k22 mini kxik22 1 n nX
i=1
krW1f(xi,W ⇤)k22  kW ⇤1 k22
mini kxik2 flatness(W ⇤).
(4)
The left hand side of (4) is usually used to regularize the Lipschitz constant of the model function, and such regularization can improve the generalization performance and adversarial robustness of the model. Hence, (4) reveals the regularization effect of flat minima. Later in the paper, we extend the analysis to higher-order moments of the gradient and combine the results with the linear stability theory of SGD to explain the implicit regularization effect of SGD. We also extend the bound on rxf(x,W ⇤) from training data to all x in a neighborhood of the training data, which implies the regularization of flatness is actually stronger than the left-hand-side of (4).
The strong correlation between rW f(x,W ) and rxf(x,W ) is justified by numerical experiments in practical settings. Specifically, in experiments we compare the following two quantities:
gW :=
1
n
nX
i=1
krW f(xi,W )k22
! 1 2
, gx :=
1
n
nX
i=1
krxf(xi,W )k22
! 1 2
. (5)
Figure 1 shows the results for a fully-connected network trained on FashionMNIST dataset and a VGG-11 network trained on CIFAR10 dataset. In each plot, gW and gx of different solutions found by SGD are shown by scatter plots. The plots show strong correlations between the two quantities. The colors of the points show that SGD with big learning rate and small batch size tends to converge to solutions with small gW and gx, which is consistent with our theoretical results on the Sobolev regularization effect of SGD. To summarize, the main contributions of this paper are:
1. We extend the linear stability analysis of SGD to high-order moments of the iterators. At the solutions selected by SGD, we find a class of conditions satisfied by the gradients of different training data. These conditions cover the flatness and non-uniformity, and also include higher order moments. They characterize the regularization effect of SGD beyond GD and SDE. 2. By exploring the multiplicative structure of the neural networks’ input layer, we build relations between the model function’s derivatives with respect to the parameters and with respect to the inputs. By these relations we turn the conditions obtained for SGD into bounds of different Sobolev (semi)norms of the model function. In particular, we show that flatness of the minimum regularizes the L2 norm of the gradient of the model function. This explains how the flatness (as well as other stability conditions for SGD) benefits generalization and adversarial robustness. 3. Still using the multiplicative structure, the bounds for Sobolev seminorms can be extended from the training data to a neighborhood around the training data, based on certain smoothness assumption of the model function (with respect to parameters). Then, bounds for generalization error and adversarial robustness are provided under reasonable assumptions of the data distribution. The bounds work well when the data are distributed effectively in a set that consists of low dimensional manifolds.

heading: 2 Preliminaries
text: Basic notations For any vector x = (x1, ..., xn)T 2 Rn and any p 1, kxkp is the conventional p-norm of x. If g(x) : Rd ! Rn is a function with vector output, kg(x)kp means the p-norm of the vector g(x), instead of a function norm. For any matrix A 2 Rm⇥n, kAkp is the operator norm
induced by the p-norm on Rn and Rm. kAkF denotes the Frobenius norm. Let g : Rn ! R be a function, ⌦ ⇢ Rn be a set in which g is defined, q 2 N⇤, and p 1. The Sobolev seminorm |g|q,p,⌦ is defined as |g|q,p,⌦ = ( P |↵|=q R ⌦ |D ↵ g(x)|pdx) 1 p , where ↵ = (↵1, ...,↵n) is a multi-index with n positive integers, |↵| = Pn
i=1 ↵i, and D ↵ g is defined as D↵g = @ |↵|g @x
↵1 1 ···@x ↵n n . The index q can be extended to fractions, but in this paper we always consider integers. When ⌦ = {x1, ...,xn} is a finite set, we define the Sobolev seminorm as |g|q,p,⌦ = ( P |↵|=q 1 n Pn i=1 |D↵g(xi)|p) 1 p .
For a symmetric matrix A 2 Rn⇥n, we let ⇤(A) be the set of eigenvalues of A and max(A) be the maximum eigenvalue of A. For two matrices A and B, A ⌦ B denotes the Kronecker product of A and B. For k 2 N⇤, let A⌦k = A ⌦ A ⌦ · · · ⌦ A where A is multiplied for k times. Therefore, if A 2 Rm⇥n, then A⌦k 2 Rmk⇥nk . However, for a vector x 2 Rn, with an abuse of notation sometimes we also use x⌦k to denote the rank-one tensor product. In this case, x⌦k 2 Rn⇥n⇥···⇥n. Problem settings In this paper, we consider the learning of the parameterized model f(x,W ), with x being the input data and W being the parameters. Let d be the dimension of the input data and w be the number of parameters in W , i.e. x 2 Rd and W 2 Rw. We assume that the model has scalar output. We consider the models with multiplicative structure on input data and part of the parameters. With an abuse of notation let W = (W1,W2), where W1 2 Rm⇥d is the part of the parameters multiplied with x and W2 denotes other parameters. Then, the model can be written as
f(x,W ) = f̃(W1x,W2). (6)
We remark that most neural networks have form (6), including convolutional networks, residual networks, recurrent networks, and even transformers. Note that since the x in (6) can also be understood as fixed features calculated using input data, (6) also includes random feature models.
We consider a supervised learning setting. Let µ be a data distribution supported within Rd, and f⇤ : Rd ! R be a target function. A set of n training data {(xi, yi)}ni=1 is obtained by sampling x1, ...,xn i.i.d. from µ, and letting yi = f⇤(xi). As mentioned in the introduction, the model is learned by minimizing the empirical loss function L̂(W ) = 12n Pn i=1(f(xi,W ) yi)2 in the space of parameters.
The population loss, or the generalization error, is defined as L(W ) = R (f(x,W ) f⇤(x))2dµ(x).
For any i = 1, 2, ..., n, let Li(W ) = 12 (f(xi,W ) yi) 2 be the loss at xi. Then, the iteration scheme of the SGD with learning rate ⌘ and batch size B is
Wt+1 = Wt ⌘
B
BX
i=1
rL⇠ti (Wt), (7)
where ⇠t = (⇠t1, ⇠t2, ..., ⇠tB) is a B-dimensional random variable uniformly distributed on the B-tuples in {1, 2, ..., n} and independent with Wt. In the paper, we study interpolation solutions W ⇤ found by SGD, which satisfies f(xi,W ⇤) = yi for any i = 1, 2, ..., n. Obtaining interpolation solutions is possible in the over-parameterized setting [34], and is widely studied in existing work [36, 33].
Some assumptions will be made when deriving the generalization error bounds. Firstly, we assume the model function with respect to the parameter W to be smooth around W ⇤. Definition 1. Let , C be positive numbers, and k be a positive integer. We say the model f(x,W ) satisfies the (C, , k)-local smoothness condition at data x and parameter W , if for any W 0 such
that kW 0 Wk2  , there is krW f(x,W 0)k2k  C(krW f(x,W )k2k + 1). (8)
Remark 1. In the definition above, we consider the 2k-norm of the gradients for the convenience of later analysis. When k = 1, the condition (8) becomes krW f(x,W 0)k2  C(krW f(x,W )k2+1). This is weaker than local approximation by Taylor expansion, which usually yields results like krW f(xi,W ) rW f(xi,W ⇤)k2  C . Here, we only require that the gradient with respect to the parameters does not get exceedingly large when W is close to W ⇤ .
Next, we need the data distribution µ to support roughly on a low dimensional manifold (or a union of low dimensional manifolds), thus the neighborhoods of training data can well cover all test data. Definition 2. Let µ be a probability distribution supported in Rd, and {xi} = {x1, ...,xn} be a set of points in Rd. For positive constants and ", we say {xi} ( , ")-covers µ if
Px⇠µ ✓
min 1in
kx xik2 > ◆ < ", (9)
i.e. with high probability a point sampled from µ lies close to a point in {xi}. Definition 3. Let {xi} = {x1, ...,xn} be a set of n points in Rd. For positive constants and K, we say {xi} is ( ,K)-scattered if the function (x) := Pn i=1 1B(xi, )(x) satisfies (x)  K. Here, B(xi, ) is the ball in Rd centered at xi with radius , and 1A(x) is the indicator function of set A. Remark 2. The uniform upper bound (x)  K in the above definition can be weakened as an integral upper bound
1 VX R X (x)dx  K, where X = Sn i=1 B(xi, ) and VX is the Lebesgue
volume of X . Definition 4. Let µ be a probability distribution supported in Rd. For positive integer n and positive constants , "1, "2, we say µ satisfies the (n, , "1, "2)-covered condition, if with probability at least 1 "1 over the choice of n i.i.d. sampled data from µ, {x1, ...,xn}, we have
Px⇠µ ✓
min 1in
kx xik2 > ◆ < "2. (10)
On the other hand, for positive integer n and positive constants , "1 and K, we say µ satisfies the (n, , "1,K)-scattered condition, if with probability at least 1 "1 over the choice of n i.i.d. sampled data {xi} from µ, {xi} is ( ,K)-scattered.
Later in the analysis of the generalization error (Section 5), we will assume the model and the data distribution satisfy the conditions in Definition 1 and 4 with appropriately chosen constants.

heading: 3 Linear stability theory of SGD
text: Compared with full-batch GD, SGD adds in each iteration a random “noise” to the gradient of the loss function. Hence, it is harder for SGD to be stable around minima than GD. This is shown by the linear stability theory of SGD.
Recall that the iteration scheme of SGD with learning rate ⌘ and batch size B is given by (7). Let W ⇤ be an interpolation solution for the learning problem. When {Wt} is close to W ⇤, the behavior of (7), including its stability around the minimum, can be characterized by the linearized dynamics at W ⇤:
W̃t+1 = W̃t ⌘
B
BX
i=1
rW f(x⇠ti ,W ⇤)rW f(x⇠ti ,W ⇤)T (W̃t W ⇤). (11)
The linearization is made by considering the quadratic approximation of Li near W ⇤. For ease of notation, let ai = rW f(xi,W ⇤), Hi = aiaTi . In the linearized dynamics (11), without loss of generality we can set W ⇤ = 0 by replacing W̃t by W̃t W ⇤. Then (11) becomes
W̃t+1 = W̃t ⌘
B
BX
i=1
a⇠tia T ⇠ti W̃t =
I ⌘
B
BX
i=1
H⇠ti
! W̃t. (12)
Next, we define the stability of the above dynamics.
Definition 5. For any k 2 N⇤, we say W ⇤ is k-th order linearly stable for SGD with learning rate ⌘ and batch size B, if there exists a constant C (which may depend on k) that satisfies
EW̃⌦kt F  C EW̃⌦k0 F ,
for any W̃t, t 2 N⇤, given by the dynamics (12) from any initialization distribution of W̃0.
{W̃t} is the trajectory of the linearized dynamics. It characterizes the performance of SGD around the minimum W ⇤, but does not equal to the original trajectory of SGD {Wt}. In practice, quadratic approximation of the loss function works in a neighborhood around the minimum. Hence, linear stability characterizes the local behavior of the iterators around the minimum. The W̃0 in Definition 5 is the starting point of the linearized dynamics. It can be a point in the trajectory of the real dynamics, but it should not be understood as the initialization of the real dynamics (W0). Since W̃t follows a linear dynamics, the W̃0 here does not need to be close to W ⇤. Without linear stability of {W̃t}, {Wt} can never converge to the minimum. It is possible that {Wt} oscillates around the minimum at a distance out of the reach of local quadratic approximation. However, empirical observations in [11] show that SGD usually converges to and oscillates in a region with good quadratic approximation. Hence, linear stability plays an important role in practice.
In [30], the following condition on the linear stability for k = 2 is provided. Proposition 1. (Theorem 1 in [30]) The global minimum W ⇤ is 2nd-order linearly stable for SGD
with learning rate ⌘ and batch size B if
max
⇢ (I ⌘H)2 + ⌘
2(n B) B(n 1) ⌃  1, (13)
where H = 1n nP
i=1 Hi, and ⌃ =
1 n nP i=1 H 2 i H2.
Shown by the proposition, the spectra of H and ⌃ influence the linear stability of SGD. Thus, in [30] the biggest eigenvalue of H and ⌃ are named sharpness and non-uniformity, respectively. The condition (13) is then be relaxed into a condition of sharpness and non-uniformity. Similar analysis was also conducted in [9] for linear least squares problems.
Remark 3. The definition of stability in [30] concerns EkW̃tk2, which is slightly different from the definition above when k = 2. However, since EkW̃tk2 and kEW̃⌦2t kF are equivalent, i.e. EkW̃tk2/ p w  kEW̃⌦2t kF  EkW̃tk2, the two ways to define linear stability are equivalent.
In the current work, we analyze the dynamics of EW̃⌦kt , which is a closed linear dynamics. As a comparison, the dynamics of EkW̃tk2 is not closed, i.e. EkW̃t+1k2 is not totally determined by EkW̃tk2. By (12) we have
W̃ ⌦k t+1 = I ⌘
B
BX
i=1
H⇠ti
!⌦k W̃
⌦k t . (14)
Let IB be the set of all subsets of {1, 2, ..., n} with cardinality B, and I = {i1, ..., iB} 2 IB be a batch. Note that ⇠t is independent with W̃t. Taking expectation for (14) gives
EW̃⌦kt+1 = E I ⌘
B
BX
i=1
H⇠ti
!⌦k EW̃⌦kt =
1 n B
X
I2I
0
@I ⌘ B
BX
j=1
Hij
1 A ⌦k
EW̃⌦kt (15)
Denote
T ⌘,B I,k :=
0
@I ⌘ B
BX
j=1
Hij
1 A ⌦k
(16)
for each batch I, and T ⌘,Bk = 1/ n B P I T ⌘,B I,k be the expectation of T ⌘,B I over the choice of batches. Then we have EW̃⌦kt+1 = T ⌘,B k EW̃ ⌦k t .
On the other hand, if EW̃⌦kt is understood as a tensor in Rw⇥w⇥···⇥w, then it is always a symmetric tensor (any permutation of the index does not change the entry value). Hence, it has the decomposition EW̃⌦kt = Pr i=1 iv ⌦k i , with i 2 R and vi 2 Rw [8]. Let Mk be the set of symmetric tensors in
Rw⇥w⇥···⇥w, and M+k = {A 2 M : A = Pr i=1 iv ⌦k i and i 0 for any i = 1, ..., r} be the set of “positive semidefinite” symmetric tensors. We provide the following conditions for the linear stability of SGD. Theorem 1. For any k 2 N⇤, the global minimum W ⇤ is kth order linearly stable for SGD with learning rate ⌘ and batch size B, if and only if
kT ⌘,Bk AkF  kAkF (17)
holds for any A 2 M+k if k is an even number, or for any A 2 Mk is k is an odd number.
All proofs are provided in the appendix. As a corollary, when k = 2, we can write down another sufficient condition for linear stability than Proposition 1. See Corollary 3 in the appendix.
Recall that we let ai = rW f(xi,W ⇤) and Hi = aiaTi . The linear stability conditions in Theorem 1 imply the following bound on the moments of ai. Specifically, we have Theorem 2. If a global minimum W ⇤ is k th order linearly stable for SGD with learning rate ⌘ and batch size B, then, for any j 2 {1, 2, ..., p}, we have
1
n
nX
i=1
a2ki,j  2kBk 1
⌘k , (18)
where ai,j means the jth entry of ai.
By Theorem 2, if a global minimum is stable for SGD with some order, then the gradient moment of this order is controlled. Summing j from 1 to w, we have 1n Pn i=1 krW f(xi,W )k 2k 2k  2w(2B)k 1
⌘k . For k = 1 and 2, this gives control to the flatness and non-uniformity of the minimum, respectively. For general k, applying the Hölder inequality on (18), we have the following corollary which bounds the mean 2k-norm of the gradients at the training data. Corollary 1. If a global minimum W ⇤ is k th order linearly stable for SGD with learning rate ⌘ and
batch size B, then
1
n
nX
i=1
krW f(xi,W ⇤)k2k = 1
n
nX
i=1
kaik2k  ⇣ w
B
⌘ 1 2k s 2B
⌘ . (19)

heading: 4 The Sobolev regularization effect
text: In this section, we build connection between rW f(x,W ) and rxf(x,W ) using the multiplication of W1 and x. For general parameterized model, it is hard to build connection between the two gradients. However, this is possible for neural network models, in which the input variable is multiplied with a set of parameters (the first-layer parameter) before any non-linear operation. By this connection, at an interpolation solution that is stable for SGD, we can turn the moments bounds derived in the previous section into the bounds on rxf(x,W ). For different k, the moment bound on rW f(x,W ) controls the Sobolev seminorms of f(·,W ) at the training data with different index p. When k = 1, this becomes a bound for the sum of gradient square at the training data, which is used in the literature as a regularization term to regularize the Lipschitz constant of the model function. This explains the favorable generalization performance of flat minima.
Recall f(x,W ) = f̃(W1x,W2), and we can express rxf̃ and rW1 f̃ as (1). Then, we easily have the following proposition.
Proposition 2. Consider the model f̃(W1x,W2). For any k 2 N⇤, any W = (W1,W2) and x 6= 0,
krxf̃(W1x,W2)k2k2k  kWT1 k2k2k kxk2k2k krW f̃(W1x,W2)k2k2k, (20)
As a corollary, we have the following bound for the Sobolev seminorm using rW f̃ .
Corollary 2. Let f(x,W ) = f̃(W1x,W2), and {xi} = {x1,x2, ...,xn}. Let |f(·,W )|1,2k,{xi} be the Sobolev seminorm of f(·,W ) on {xi} as defined at the beginning of Section 2. Then,
|f(·,W )|1,2k,{xi}  kWT1 k2k
min1in kxik2k
1
n
nX
i=1
krW f(xi,W )k2k2k
! 1 2k
. (21)
Combining Corollary 2 with the stability condition in Theorem 2, we have the following control for the Sobolev seminorm of the model function at interpolation solutions that are stable for SGD. Theorem 3. If a global minimum W ⇤ which interpolates the training data is k th order linearly stable for SGD with learning rate ⌘ and batch size B, and W ⇤ = (W ⇤1 ,W ⇤ 2 ). Then,
|f(·,W ⇤)|1,2k,{xi}  k(W ⇤1 )T k2k
min1in kxik2k
⇣ w
B
⌘ 1 2k s 2B
⌘ . (22)
By Corollary 2 and Theorem 3, as long as the data are not very small (in practical problems, the input data are usually normalized), and the input-layer parameters are not very large, the Sobolev seminorms of the model function (evaluated at the training data) is regularized by the linear stability of SGD. The regularization effect on Sobolev seminorm gets stronger for bigger learning rate and smaller batch size. When k is big, the dependence of the bound with p is negligible.
If the solution W ⇤ satisfies the local smoothness condition defined in Definition 1, the control on rxf̃(W ⇤1 x,W ⇤2 ) can be extended to a neighborhood around the training data. Then, we will be able to control the “population” Sobolev seminorms. First, around a certain training data, we have the following estimate. Proposition 3. Assume the model f(x,W ) satisfies (C, approx, k)-local smoothness condition for some C, approx > 0 and k 2 N⇤, at W ⇤ and x⇤. Then, for any x that satisfies kx x⇤k2  approxkx⇤k2/kW ⇤1 k2, we have
krxf(x,W ⇤)k2k  Ck(W ⇤1 )T k2k
kx⇤k2k (krW f(x⇤,W ⇤)k2k + 1) . (23)
Proposition 3 turns the local smoothness condition in the parameter space into a local smoothness condition in the data space. This is made possible by the multiplicative structure of the network’s input layer. Specifically, since in the first layer W1 and x are multiplied together, a perturbation of x can be turned to a perturbation of W1 without changing the value of their product. See the proof of the proposition in Appendix B.2 for more detail.
By the results above, we can obtain the following theorem which estimates the Sobolev seminorm of the model function on a neighborhood of the training data. Theorem 4. Let W
⇤ = (W ⇤1 ,W ⇤ 2 ) be an interpolation solution that is k th order linearly stable for
SGD with learning rate ⌘ and batch size B. Let {x1, ...,xn} be the training data. Suppose the model f(x,W ) satisfies (C, approx, k)-local smoothness condition at W ⇤ and xi for any i = 1, 2, ..., n. Consider the set X := Sn i=1 {x : kx xik  } with  approxkxik2/kW ⇤1 k2. Assume {x1, ...,xn} satisfies ( ,K)-scattered condition as defined in Definition 3. Then, we have
|f(x,W ⇤)|1,2k,X  (KVX )
1 2k 2Ck(W ⇤1 )T k2k
mini kxik2k
⇣ w
B
⌘ 1 2k s 2B
⌘ + 1
! , (24)
where VX is the Lebesgue volume of X .
By Theorem 4, the model function found by SGD is smooth in a neighborhood of the training data, given that the landscape of the model in the parameter space is smooth. When k = 1, the results implies that flatness regularizes |rxf(x,W ⇤)|1,2,X , which is stronger than the gradient square at training data.
Discussion on the tightness of the bounds The bounds in this section depend on the norm of W ⇤1 . This norm may be big, especially when homogeneous activation functions, such as ReLU, is used. In
the ReLU case, for example, we can simply scale the first layer parameters to be arbitrarily large and second layer small while not changing the model function.
However, in practice SGD does not find these solutions with unbalanced layers because the homogeneity of ReLU induces an invariant in parameters. Specifically, let Wl and Wl 1 be the parameter matrices in the lth and (l 1)th layers, then WTl 1Wl 1 WlWTl does not change significantly throughout the training (it is actually fixed in the zero learning rate limit). This invariant keeps layers balanced and prevents the first layer parameters from getting big alone.
On the other hand, we conduct numerical experiments to check the norm of W1 during training. Table 1 shows the average l2-norm of W1 for the same experiments in part (a) of Figure1. The model is a fully-connected neural network, and the dataset is FashionMNIST. Results show that the norm of W1 increases at the beginning, but becomes stable afterwards. It will not keep increasing during the training.

heading: 5 Generalization error and adversarial robustness
text: Regularizing the smoothness of the model function can improve generalization performance and adversarial robustness. This is confirmed in many practical studies [24, 29, 26]. Intuitively, a function with small gradient changes slowly as the input data changes, hence when the test data is close to one of the training data, the prediction error on the test data will be small. In this section, going along this intuition, we provide theoretical analysis of the generalization performance and adversarial robustness based on the Sobolev regularization effect derived in the previous sections. Still consider the model f(x,W ) = f̃(W1x,W2), and an interpolation solution W ⇤ = (W ⇤1 ,W ⇤2 ) found by SGD. First, based on Proposition 3 and Theorem 2, we show the following theorem, which is similar to Theorem 4 but estimates the maximum value of krxf(x,W ⇤)k in X . Theorem 5. Let W⇤ = (W ⇤1 ,W ⇤2 ) be an interpolation solution that is kth order linearly stable for SGD with learning rate ⌘ and batch size B. Let {x1, ...,xn} be the training data. Suppose the model f(x,W ) satisfies (C, approx, k)-local smoothness condition at W ⇤ and xi for any 1 = 1, 2, ..., n. Recall the definition of X in Theorem 4. Then, for any x 2 X , we have
krxf(x,W ⇤)k2k  Ck(W ⇤1 )T k2k mini kxik2k
✓ 2nw
B
◆ 1 2k s 2B
⌘ + 1
! . (25)
With the result above, we can bound the generalization error if most test data are close to the training data, i.e. the data distribution µ satisfies a covered condition. This happens for machine learning problems with sufficient training data, especially for those where the training data lie approximately on a union of some low-dimensional surfaces. This is common in practical problems [1].
Theorem 6. Suppose parameterized model f(x,W ) = f̃(W1x,W2) is used to solve a supervised learning problem with data distribution µ. Let f ⇤ be the target function. Assume µ sat-
isfies (n, , "1, "2)-covered condition. Assume with probability no less than 1 "1 SGD with learning rate ⌘ and batch size B can find an interpolation solution W ⇤ which is k th order linearly stable, and (C, approx, k)-local smoothness condition is satisfied at W ⇤ and x1, ...,xn with approx kW ⇤1 k2/mini{kxik2}. Further assume that both |f(·,W ⇤)| and |f⇤(·)| are upper bounded by M1, and krxf⇤(·)k2 is upper bounded by M2, for some constants M1 and M2. Then, with probability 1 2"1 over the sampling of training data {xi}ni=1, we have
Ex⇠µkf(x,W ⇤) f⇤(x)k22  2dC2k(W ⇤1 )T k22k
mini kxik22k
✓ 2nw
B
◆ 1 2k s 2B
⌘ + 1
!2 2+2M22 2+4M21 "2.
(26)
The bound (26) tends to 0 as n ! 1 as long as decays faster than n 12k . From a geometric perspective, this happens when the dimension of the support of µ is less than k. And the lower the dimension, the faster the decay. It is possible to get rid of the n dependency in the bound when k is sufficiently large. We may use the estimate of Sobolev functions with scattered zeros, e.g. Theorem 4.1 in [2]. We leave the analysis to future work. Adversarial robustness Neural network models suffer from adversarial examples [28], because the model function changes very fast in some directions so the function value becomes very different after a small perturbation. However, the results in Theorem 5 directly imply the adversarial robustness of the model at the training data. Hence, flatness, as well as high-order linear stability conditions, also imply adversarial robustness. Specifically, we have the following theorem. Theorem 7. Let W⇤ = (W ⇤1 ,W ⇤2 ) be an interpolation solution that is kth order linearly stable for SGD with learning rate ⌘ and batch size B. Let {x1, ...,xn} be the training data. Suppose the model f(x,W ) satisfies (C, approx, k)-local smoothness condition at W ⇤ and xi for any 1 = 1, 2, ..., n. Then, for any x that satisfies kx xik2  for some i 2 {1, ..., n} and  approx mini kxik2/kW ⇤1 k2, we have
|f(x,W ⇤) f⇤(xi,W ⇤)|  C p dk(W ⇤1 )T k2k
mini kxik2k
✓ 2nw
B
◆ 1 2k s 2B
⌘ + 1
! . (27)

heading: 6 Conclusion
text: In this paper, we connect the linear stability theory of SGD with the generalization performance and adversarial robustness of the neural network model. As a corollary, we provide theoretical insights of why flat minimum generalizes better. To achieve the goal, we explore the multiplicative structure of the neural network’s input layer, and build connection between the model’s gradient with respect to the parameters and the gradient with respect to the input data. We show that as long as the landscape on the parameter space is mild, the landscape of the model function with respect to the input data is also mild, hence the flatness (as well as higher order linear stability conditions) has the effect of Sobolev regularization. Our study reveals the significance of the multiplication structure between data (or features in intermediate layers) and parameters. It is an important source of implicit regularization of neural networks and deserves further exploration.

Review summaries:
The paper extends previous stability analysis ([30]) to high-order moments of the iterator. They characterize the implicit regularization of SGD by showing that a stable global minimum must satisfy conditions regarding its flatness and non-uniformity. In addition, the authors build relations between the model function’s derivative w.r.t. the parameters and w.r.t. to the inputs. Using these relations, they provide upper bounds for the Sobolev seminorm of the model function at stable SGD interpolating solutions. In addition, under additional smoothness conditions, the authors were able to extend that last result to a neighborhood of the training data. Finally, the authors built upon the Sobolev regularization and use it to provide bounds on the generalization performance and adversarial robustness.

Prior work on minima stability for SGD considered linear stability in expectation using second order moment. This paper considers linear stability in higher moments, while providing necessary and sufficient conditions for it. Using these conditions and the multiplicative structure of neural networks, the authors show that the Sobolev semi-norm of SGD solutions (w.r.t. the input) is bounded from above. Then, under certain smoothness assumptions on the learned function and the target function, the authors prove generalization bounds and adversarial robustness.

This paper shows that due to the multiplicative structure of input data and first layer parameters in a MLP neural network the flat minima regularizes the norm of the gradient of the model and then explains its generalizability. It also shows that SGD actually regularizes higher-order moments of the gradient of the noise when it is close to an interpolating minima. Utilizing this property of SGD and multiplicative structure shows that SGD minimizes the Sobolev seminorms of the model. The presented analysis is based on quadratic estimation of the objective around the minima and uses linear stability analysis.