Title: EMPIRICAL RISK LANDSCAPE ANALYSIS FOR UNDER-STANDING DEEP NEURAL NETWORKS

Abstract: This work aims to provide comprehensive landscape analysis of empirical risk in deep neural networks (DNNs), including the convergence behavior of its gradient, its stationary points and the empirical risk itself to their corresponding population counterparts, which reveals how various network parameters determine the convergence performance. In particular, for an l-layer linear neural network consisting of d i neurons in the i-th layer, we prove the gradient of its empirical risk uniformly converges to the one of its population risk, at the rate of
Here d is the total weight dimension, s is the number of nonzero entries of all the weights and the magnitude of weights per layer is upper bounded by r. Moreover, we prove the one-to-one correspondence of the non-degenerate stationary points between the empirical and population risks and provide convergence guarantee for each pair. We also establish the uniform convergence of the empirical risk to its population counterpart and further derive the stability and generalization bounds for the empirical risk. In addition, we analyze these properties for deep nonlinear neural networks with sigmoid activation functions. We prove similar results for convergence behavior of their empirical risk gradients, non-degenerate stationary points as well as the empirical risk itself. To our best knowledge, this work is the first one theoretically characterizing the uniform convergence of the gradient and stationary points of the empirical risk of DNN models, which benefits the theoretical understanding on how the neural network depth l, the layer width d i , the network size d, the sparsity in weight and the parameter magnitude r determine the neural network landscape.
In this work, we consider both linear and nonlinear DNNs. Suppose both networks consist of l layers. We use u (j) and v (j) to respectively denote the input and output of the j-th layer, ∀j = 1, . . . , l.
Deep linear neural networks: The function of the j-th layer is formulated as u (j) W (j) v (j-1) ∈ R dj , v (j) u (j) ∈ R dj , ∀j = 1, • • • , l,

1. INTRODUCTION:

Deep learning has achieved remarkable success in many fields, such as computer vision (Hinton et al., 2006; Szegedy et al., 2015; He et al., 2016) , natural language processing (Collobert & Weston, 2008; Bakshi & Stephanopoulos, 1993) , and speech recognition (Hinton et al., 2012; Graves et al., 2013) . However, theoretical understanding on the properties of deep learning models still lags behind their practical achievements (Shalev-Shwartz et al., 2017; Kawaguchi, 2016) due to their high non-convexity and internal complexity. In practice, parameters of deep learning models are learned by minimizing the empirical risk via (stochastic-)gradient descent. Therefore, some recent works (Bartlett & Maass, 2003; Neyshabur et al., 2015) analyzed the convergence of the empirical risk to the population risk, which are however still far from fully understanding the landscape of the empirical risk in deep learning models. Beyond the convergence properties of the empirical risk itself, the convergence and distribution properties of its gradient and stationary points are also essential in landscape analysis. A comprehensive landscape analysis can reveal important information on the optimization behavior and practical performance of deep neural networks, and will be helpful to designing better network architectures. Thus, in this work we aim to provide comprehensive landscape analysis by looking into the gradients and stationary points of the empirical risk.

Formally, we consider a DNN model f (w; x, y) : R d0 × R d l → R parameterized by w ∈ R d consisting of l layers (l ≥ 2) that is trained by minimizing the commonly used squared loss function over sample pairs {(x, y)} ⊂ R d0 × R d l from an unknown distribution D, where y is the target output for the sample x. Ideally, the model can find its optimal parameter w * by minimizing the population risk through (stochastic-)gradient descent by backpropagation: min w J (w) E (x,y)∼D f (w; x, y), where f (w; x, y) = 1 2 v (l) -y 2 2 is the squared loss associated to the sample (x, y) ∼ D in which v (l) is the output of the l-th layer. In practice, as the sample distribution D is usually unknown and only finite training samples (x (i) , y (i) ) n i=1 i.i.d. drawn from D are provided, the network model is usually trained by minimizing the empirical risk:

min w Ĵn (w) 1 n n i=1 f (w; x (i) , y (i) ).

(1)

Understanding the convergence behavior of Ĵn (w) to J (w) is critical to statistical machine learning algorithms. In this work, we aim to go further and characterize the landscape of the empirical risk Ĵn (w) of deep learning models by analyzing the convergence behavior of its gradient and stationary points to their corresponding population counterparts. We provide analysis for both multi-layer linear and nonlinear neural networks. In particular, we obtain following new results.

• We establish the uniform convergence of empirical gradient ∇ w Ĵn (w) to its population counterpart ∇ w J (w). Specifically, when the sample size n is not less than O max(l 3 r 2 /(ε 2 s log(d/l)), s log(d/l)/l) , with probability at least 1 -ε the convergence rate is O(r 2l l √ max i d i s log(d/l)/n), where there are s nonzero entries in the parameter w, the output dimension of the i-th layer is d i and the magnitude of the weight parameter of each layer is upper bounded by r. This result implies that as long as the training sample size n is sufficiently large, any stationary point of Ĵn (w) is also a stationary point of J (w) and vise versa, although both Ĵn (w) and J (w) are very complex.

• We then prove the exact correspondence of non-degenerate stationary points between Ĵn (w)

and J (w). Indeed, the corresponding non-degenerate stationary points also uniformly converge to each other at the same convergence rate as the one revealed above with an extra factor 2/ζ. Here ζ > 0 accounts for the geometric topology of non-degenerate stationary points (see Definition 1).

Based on the above two new results, we also derive the uniform convergence of the empirical risk Ĵn (w) to its population risk J (w), which helps understand the generalization error of deep learning models and stability of their empirical risk. These analyses reveal the role of the depth l of a neural network model in determining its convergence behavior and performance. Also, the results tell that the width factor √ max i d i , the nonzero entry number s of weights, and the total network size d are also critical to the convergence and performance. In addition, controlling magnitudes of the parameters (weights) in DNNs are demonstrated to be important for performance. To our best knowledge, this work is the first one theoretically characterizing the uniform convergence of empirical gradient and stationary points in both deep linear and nonlinear neural networks.

2. RELATED WORK:

To date, only a few theories have been developed for understanding DNNs which can be roughly divided into following three categories. The first category aims to analyze training error of DNNs. Baum (1988) pointed out that zero training error can be obtained when the last layer of a neural network has more units than training samples. Later, Soudry & Carmon (2016) proved that for DNNs with leaky rectified linear units (ReLU) and a single output, the training error achieves zero at any of their local minima as long as the product of the number of units in the last two layers is larger than the training sample size.

The second category of analysis works (Dauphin et al., 2014; Choromanska et al., 2015a; Kawaguchi, 2016; Tian, 2017) focus on analyzing loss surfaces of DNNs, e.g., how the stationary points are distributed. Those results are helpful to understanding performance difference of large-and small-size networks (Choromanska et al., 2015b) . Among them, Dauphin et al. (2014) experimentally verified that a large number of saddle points indeed exist for DNNs. With strong assumptions, Choromanska et al. (2015a) connected the loss function of a deep ReLU network with the spherical spin-class model and described locations of the local minima. Later, Kawaguchi (2016) proved the existence of degenerate saddle points for deep linear neural networks with squared loss function. They also showed that any local minimum is also a global minimum. By utilizing techniques from dynamical system analysis, Tian (2017) gave guarantees that for two-layer bias-free networks with ReLUs, the gradient descent algorithm with certain symmetric weight initialization can converge to the ground-truth weights globally, if the inputs follow Gaussian distribution. Recently, Nguyen & Hein (2017) proved that for a fully connected network with squared loss and analytic activation functions, almost all the local minima are globally optimal if one hidden layer has more units than training samples and the network structure after this layer is pyramidal. Besides, some recent works, e.g., (Zhang et al., 2016; 2017) , tried to alleviate analysis difficulties by relaxing the involved highly nonconvex functions into ones easier.

In addition, some existing works (Bartlett & Maass, 2003; Neyshabur et al., 2015) analyze the generalization performance of a DNN model. Based on the Vapnik-Chervonenkis (VC) theory, Bartlett & Maass (2003) proved that for a feedforward neural network with one-dimensional output, the best convergence rate of the empirical risk to its population risk on the sample distribution can be bounded by its fat-shattering dimension. Recently, Neyshabur et al. (2015) adopted Rademacher complexity to analyze learning capacity of a fully-connected neural network model with ReLU activation functions and bounded inputs.

However, although gradient descent with backpropagation is the most common optimization technique for DNNs, none of existing works analyzes convergence properties of gradient and stationary points of the DNN empirical risk. For single-layer optimization problems, some previous works analyze their empirical risk but essentially differ from our analysis method. For example, Negahban et al. (2009) proved that for a regularized convex program, the minimum of the empirical risk uniformly converges to the true minimum of the population risk under certain conditions. Gonen & Shalev-Shwartz (2017) proved that for nonconvex problems without degenerated saddle points, the difference between empirical risk and population risk can be bounded. Unfortunately, the loss of DNNs is highly nonconvex and has degenerated saddle points (Fyodorov & Williams, 2007; Dauphin et al., 2014; Kawaguchi, 2016) , thus their analysis results are not applicable. Mei et al. (2017) analyzed the convergence behavior of the empirical risk for nonconvex problems, but they only considered the single-layer nonconvex problems and their analysis demands strong sub-Gaussian and subexponential assumptions on the gradient and Hessian of the empirical risk respectively. Their analysis also assumes a linearity property on gradient which is difficult to hold or verify. In contrast, our analysis requires much milder assumptions. Besides, we prove that for deep networks which are highly nonconvex, the non-degenerate stationary points of empirical risk can uniformly converge to their corresponding stationary points of population risk at the rate of O( s/n) which is faster than the rate O( d/n) for single-layer optimization problems in (Mei et al., 2017) . Also, Mei et al. (2017) did not analyze the convergence rate of the empirical risk, stability or generalization error of DNNs as this work.

3. PRELIMINARIES:

Throughout the paper, we denote matrices by boldface capital letters, e.g. A. Vectors are denoted by boldface lowercase letters, e.g. a, and scalars are denoted by lowercase letters, e.g. a. We define the r-radius ball as B d (r) {z ∈ R d | z 2 ≤ r}. To explain the results, we also need the vectorization operation vec(•). It is defined as vec(A) = (A(:, 1); • • • ; A(:, t)) ∈ R st that vectorizes A ∈ R s×t along its columns. We use d = l j=1 d j d j-1 to denote the total dimension of weight parameters, where d j denotes the output dimension of the j-th layer.

where v (0) = x is the input and W (j) ∈ R dj ×dj-1 is the weight matrix of the j-th layer.

Deep nonlinear neural networks: We adopt the sigmoid function as the non-linear activation function. The function within the j-th layer can be written as u (j) W (j) v (j-1) ∈ R dj , v (j) h j (u (j) ) = (σ(u

(j) 1 ); • • • ; σ(u (j) dj )) ∈ R dj , ∀j = 1, • • • , l,

where u (j) i denotes the i-th entry of u (j) and σ(•) is the sigmoid function, i.e., σ(a) = 1/(1 + e -a ). Following the common practice, both DNN models adopt the squared loss function defined as f (w; x, y) = 1 2 v (l) -y 2 2 , where w = (w (1) ; • • • ; w (l) ) ∈ R d contains all the weight parameters and w (j) = vec W (j) ∈ R dj dj-1 . Then the empirical risk Ĵn (w) is Ĵn (w) =

1 n n i=1 f (w; x (i) , y (i) ) = 1 2n n i=1 v (l) (i) -y (i) 2 2

, where v (l) (i) is the network's output of x (i) .

4. RESULTS FOR DEEP LINEAR NEURAL NETWORKS:

We first analyze linear neural network models and present following new results: (1) the uniform convergence of the empirical risk gradient to its population counterpart and (2) the convergence properties of non-degenerate stationary points of the empirical risk. As a corollary, we also derive the uniform convergence of the empirical risk to the population one, which further gives stability and generalization bounds. In the next section, we extend the analysis to non-linear neural network models.

We assume the input datum x is τ 2 -sub-Gaussian and has bounded magnitude, as formally stated in Assumption 1. Assumption 1. The input datum x ∈ R d0 has zero mean and is τ 2 -sub-Gaussian, i.e.,

E[exp ( λ, x )] ≤ exp 1 2 τ 2 λ 2 2 , ∀λ ∈ R d0 .

Besides, the magnitude x is bounded as x 2 ≤ r x , where r x is a positive universal constant.

Note that any random vector z consisting of independent entries with bounded magnitude is sub-Gaussian and satisfies Assumption 1 (Vershynin, 2012) . Moreover, for such a random z, we have τ = z ∞ ≤ z 2 ≤ r x . Such an assumption on bounded magnitude generally holds for natural data, e.g., images and speech signals. Besides, we assume the weight parameters w (j) of each layer are bounded as w ∈ Ω = {w | w (j) ∈ B dj dj-1 (r j ), ∀j = 1, • • • , l} where r j is a constant. For notational simplicity, we let r = max j r j . Such an assumption is common (Xu & Mannor, 2012) .

Here we assume the entry value of y falls in [0, 1]. For any bounded target output y, we can always scale it to satisfy such a requirement.

The results presented for linear neural networks here can be generalized to deep ReLU neural networks by applying the results from Choromanska et al. (2015a) and Kawaguchi (2016) , which transform deep ReLU neural networks into deep linear neural networks under proper assumptions.

4.1. UNIFORM CONVERGENCE OF EMPIRICAL RISK GRADIENT:

We first analyze the convergence of gradients for the DNN empirical and population risks. To our best knowledge, these results are the first ones giving guarantees on gradient convergence, which help better understand the landscape of DNNs and their optimization behavior. The results are stated blow. Theorem 1. Suppose Assumption 1 on the input datum x holds and the activation functions in a deep neural network are linear. Then the empirical gradient uniformly converges to the population gradient in Euclidean norm. Specifically, there exist two universal constants c g and c g such that if n ≥ c g max(l 3 r 2 r 4 x /(c q s log(d/l)ε 2 τ 4 log(1/ε)), s log(d/l)/(lτ 2 )) where c q = max 0≤i≤l d i , then

sup w∈Ω ∇ Ĵn (w) -∇J (w) 2 ≤ g c g τ ω g lc q

s log(dn/l) + log(12/ε) n holds with probability at least 1 -ε, where s denotes the number of nonzero entries of all weight parameters and ω g = max τ r 2l-1 , r 2l-1 , r l-1 .

From Theorem 1, one can observe that with an increasingly larger sample size n, the difference between empirical risk and population risk gradients decreases monotonically at the rate of O(1/ √ n) (up to a log factor). Theorem 1 also characterizes how the depth l contributes to obtaining small difference between the empirical and population risk gradients. Specifically, a deeper neural network needs more training samples to mitigate the difference. Also, due to the factor d, training a network of larger size using gradient descent also requires more training samples. We observe a factor of √ max i d i (i.e. c q ), which prefers a DNN architecture of balanced layer sizes (without extremely wide layers). This result also matches the trend and empirical performance in deep learning applications advocating deep but thin networks (He et al., 2016; Szegedy et al., 2015) .

By observing Theorem 1, imposing certain regularizations on the weight parameters is useful. For example, reducing the number of nonzero entries s encourages sparsity regularization like w 1 . The results also suggest not choosing large-magnitude weights w in order for a smaller factor r by adopting regularization like w 2 2 . Theorem 1 also reveals the point derived from optimizing that the empirical and population risks have similar properties when the sample size n is sufficiently large. For example, an /2-stationary point w of Ĵn (w) is also an -stationary point of J (w) with probability 1 -ε if n ≥ c (τ ω g / ) 2 lc q s log(d/l) with c being a constant. Here -stationary point for a function F means the point w satisfying ∇ w F 2 ≤ . Understanding such properties is useful, since in practice one usually computes an -stationary point of Ĵn (w). These results guarantee the computed point is at most a 2 -stationary point of J (w) and is thus close to the optimum.

4.2. UNIFORM CONVERGENCE OF STATIONARY POINTS:

We then proceed to analyze the distribution and convergence properties of stationary points of the DNN empirical risk. Here we consider non-degenerate stationary points which are geometrically isolated and thus unique in local regions. Since degenerate stationary points are not unique in a local region, we cannot expect to establish one-to-one corresponding relationship (see below) between them in empirical risk and population risk. Definition 1. (Non-degenerate stationary points) (Gromoll & Meyer, 1969) If a stationary point w is said to be a non-degenerate stationary point of J (w), then it satisfies

inf i λ i ∇ 2 J (w) ≥ ζ,

where λ i ∇ 2 J (w) denotes the i-th eigenvalue of the Hessian ∇ 2 J (w) and ζ is a positive constant.

Non-degenerate stationary points include local minima/maxima and non-degenerate saddle points, while degenerate stationary points refer to degenerate saddle points. Then we introduce the index of non-degenerate stationary points which can characterize their geometric properties. Definition 2. (Index of non-degenerate stationary points) (Dubrovin et al., 2012) The index of a symmetric non-degenerate matrix is the number of its negative eigenvalues, and the index of a non-degenerate stationary point w of a smooth function F is simply the index of its Hessian ∇ 2 F (w).

Suppose that J (w) has m non-degenerate stationary points that are denoted as {w (1) , w (2) , • • • , w (m) }. We prove following convergence behavior of these stationary points. Theorem 2. Suppose Assumption 1 on the input datum x holds and the activation functions in a deep neural network are linear. Then if n ≥ c h max(l 3 r 2 r 4

x /(c q s log(d/l)ε 2 τ 4 log(1/ε)), s log(d/l)/ζ 2 ) where c h is a constant, for k ∈ {1, • • • , m}, there exists a non-degenerate stationary point w (k) n of Ĵn (w) which corresponds to the non-degenerate stationary point w (k) of J (w) with probability at least 1 -ε. In addition, w (k) n and w (k) have the same non-degenerate index and they satisfy

w (k) n -w (k) 2 ≤ 2c g τ ω g ζ lc q s log(dn/l) + log(12/ε) n , (k = 1, • • • , m)

with probability at least 1 -ε, where the parameters c q , ω g , and c g are given in Theorem 1.

Theorem 2 guarantees the one-to-one correspondence between the non-degenerate stationary points of the empirical risk Ĵn (w) and the popular risk J (w). The distances of the corresponding pairs become smaller as n increases. In addition, the corresponding pairs have the same non-degenerate index. This implies that the corresponding stationary points have the same geometric properties, such as whether they are saddle points. Accordingly, we can develop more efficient algorithms, e.g. escaping saddle points (Ge et al., 2015) , since Dauphin et al. (2014) empirically proved that saddle points are usually surrounded by high error plateaus. Also when n is sufficiently large, the properties of stationary points of Ĵn (w) are similar to the points of the population risk J (w) in the sense that they have exactly matching local minima/maxima and non-degenerate saddle points. By comparing Theorems 1 and 2, we find that the requirement for sample number in Theorem 2 is more restrict, since establishing exact one-to-one correspondence between the non-degenerate stationary points of Ĵn (w) and J (w) and bounding their uniform convergence rate to each other are more challenging. From Theorems 1 and 2, we also notice that the uniform convergence rate of non-degenerate stationary points has an extra factor 1/ζ. This is because bounding stationary points needs to access not only the gradient itself but also the Hessian matrix. See more details in proof. Kawaguchi (2016) pointed out that degenerate stationary points indeed exist for DNNs. However, since degenerate stationary points are not isolated, such as forming flat regions, it is hard to establish the unique correspondence for them as for non-degenerate ones. Fortunately, by Theorem 1, the gradients at these points of Ĵn (w) and J (w) are close. This implies that a degenerate stationary point of J (w) will also give a near-zero gradient for Ĵn (w), i.e., it is also a stationary point for Ĵn (w).

In the proof, we consider the essential multi-layer architecture of the deep linear network, and do not transform it into a linear regression model and directly apply existing results (see Loh & Wainwright (2015) and Negahban et al. (2011)) . This is because we care more about deep ReLU networks which cannot be reduced in this way. Our proof technique is more suitable for analyzing the multi-layer neural networks which paves a way for analyzing deep ReLU networks. Also such an analysis technique can reveal the role of network parameters (dimension, norm, etc.) of each weight matrix in the results which may benefit the design of networks. Besides, the obtained results are more consistent with those for deep nonlinear networks (see Sec. 5).

4.3. UNIFORM CONVERGENCE, STABILITY AND GENERALIZATION OF EMPIRICAL RISK:

Based on the above results, we can derive the uniform convergence of empirical risk to population risk easily. In this subsection, we first give the uniform convergence rate of empirical risk for deep linear neural networks in Theorem 3, and then use this result to derive the stability and generalization bounds for DNNs in Corollary 1. Theorem 3. Suppose Assumption 1 on the input datum x holds and the activation functions in a deep neural network are linear. Then there exist two universal constants c f and

c f such that if n ≥ c f max(l 3 r 4 x /(d l s log(d/l)ε 2 τ 4 log(1/ε)), s log(d/l)/(τ 2 d l )), then sup w∈Ω Ĵn (w) -J (w) ≤ f c f τ max d l τ r 2l , r l s log(dn/l) + log(8/ε) n (2)

holds with probability at least 1 -ε. Here l is the number of layers in the neural network, n is the sample size and d l is the dimension of the final layer.

From Theorem 3, when n → +∞, we have | Ĵn (w) -J (w)| → 0. According to the definition of uniform convergence (Vapnik & Vapnik, 1998; Shalev-Shwartz et al., 2010) , under the distribution D, the empirical risk of a deep linear neural network converges to its population risk uniformly at the rate of O(1/ √ n). Theorem 3 also explains the roles of the depth l, the network size d, and the number of nonzero weight parameters s in a DNN model.

Based on VC-dimension techniques, Bartlett & Maass (2003) proved that for a feedforward neural network with polynomial activation functions and one-dimensional output, with probability at least

1 -ε the convergence bound satisfies | Ĵn (w) -inf f J (w)| ≤ O( (γ log 2 (n) + log(1/ε))/n).

Here γ is the shattered parameter and can be as large as the VC-dimension of the network model, i.e. at the order of O(ld log(d)+l 2 d) (Bartlett & Maass, 2003) . Note that Bartlett & Maass (2003) did not reveal the role of the magnitude of weight in their results. In contrast, our uniform convergence bound is sup w∈Ω | Ĵn (w)-J (w)| ≤ O( (s log(dn/l) + log(1/ε))/n). So our convergence rate is tighter. Neyshabur et al. (2015) proved that the Rademacher complexity of a fully-connected neural network model with ReLU activation functions and one-dimensional output is O r l / √ n (see Corollary 2 in (Neyshabur et al., 2015) ). Then by applying Rademacher complexity based argument (Shalev-Shwartz & Ben-David, 2014a), we have | sup f ( Ĵn (w) -J (w))| ≤ O((r l + log(1/ε))/ √ n) with probability at least 1 -ε where the loss function is the training error g = 1 (v (l) =y) in which v (l) is the output of the l-th layer in the network model f (w; x, y). The convergence rate in our theorem is O(r 2l (s log(d/l) + log(1/ε))/n) and has the same convergence speed O(1/ √ n) w.r.t. sample number n. Note that our convergence rate involves r 2l since we use squared loss instead of the training error in (Neyshabur et al., 2015) . The extra parameters s and d are involved since we consider the parameter space rather than the function hypothesis f in (Neyshabur et al., 2015) , which helps people more transparently understand the roles of the network parameters. Besides, the Rademacher complexity cannot be applied to analyzing convergence properties of the empirical risk gradient and stationary points as our techniques.

Based on Theorem 3, we proceed to analyze the stability property of the empirical risk and the convergence rate of the generalization error in expectation. Let S = {(x (1) , y (1) ), • • • , (x (n) , y (n) )} denote the sample set in which the samples are i.i.d. drawn from D. When the optimal solution w n to problem (1) is computed by deterministic algorithms, the generalization error is defined as g = Ĵn (w n ) -J (w n ). But one usually employs randomized algorithms, e.g. stochastic gradient descent (SGD), for computing w n . In this case, stability and generalization error in expectation defined in Definition 3 are more applicable. Definition 3. (Stability and generalization in expectation) (Vapnik & Vapnik, 1998; Shalev-Shwartz et al., 2010; Gonen & Shalev-Shwartz, 2017 ) Assume a randomized algorithm A is employed,

((x (1) , y (1) ), • • • , (x (n) , y (n) )) ∼ D and w n = argmin w Ĵn (w) is the empirical risk minimizer (ERM). For every j ∈ [n], suppose w j * = argmin w 1 n-1 i =j f i (w; x (i) , y (i)

). We say that the ERM is on average stable with stability rate k under distribution D if E S∼D,A,(x (j) ,y (j) )∼D 1 n n j=1 f j (w j * ; x (j) , y (j) ) -f j (w n ; x (j) , y (j) ) ≤ k . The ERM is said to have generalization error with convergence rate k under distribution D if we have

E S∼D,A J (w n ) -Ĵn (w n ) ≤ k .

Stability measures the sensibility of the empirical risk to the input and generalization error measures the effectiveness of ERM on new data. Generalization error in expectation is especially important for applying DNNs considering their internal randomness, e.g. from SGD optimization. Now we present the results on stability and generalization performance of deep linear neural networks. Corollary 1. Suppose Assumption 1 on the input datum x holds and the activation functions in a deep neural network are linear. Then with probability at least 1 -ε, both the stability rate and the generalization error rate of ERM of a deep linear neural network are at least f :

E S∼D,A,(x (j) ,y (j) )∼D 1 n n j=1 f * j -f j ≤ f and E S∼D,A J (w n ) -Ĵn (w n ) ≤ f ,

where f * j and f j respectively denote f j (w j * ; x (j) , y (j) ) and f j (w n ; x (j) , y (j) ), and f is defined in Eqn. (2).

According to Corollary 1, both the stability rate and the convergence rate of generalization error are O( f ). This result indicates that deep learning empirical risk is stable and its output is robust to small perturbation over the training data. When n is sufficiently large, small generalization error of DNNs is guaranteed.

5. RESULTS FOR DEEP NONLINEAR NEURAL NETWORKS:

In the above section, we analyze the empirical risk optimization landscape for deep linear neural network models. In this section, we extend our analysis to deep nonlinear neural networks which adopt the sigmoid activation function. Our analysis techniques are also applicable to other third-order differentiable activation functions, e.g., tanh function with different convergence rate. Here we assume the input data are i.i.d. Gaussian variables. Assumption 2. The input datum x is a vector of i.i.d. Gaussian variables from N (0, τ 2 ).

Since for any input, the sigmoid function always maps it to the range [0, 1]. Thus, we do not require the input x to have bounded magnitude. Such an assumption is common. For instance, Tian (2017) and Soudry & Hoffer (2017) both assumed that the entries in the input vector are from Gaussian distribution. We also assume w ∈ Ω as in (Xu & Mannor, 2012) . Here we also assume that the entry value of the target output y falls in [0, 1]. Similar to the analysis of deep linear neural networks, here we also aim to characterize the empirical risk gradient, stationary points and empirical risk for deep nonlinear neural networks.

5.1. UNIFORM CONVERGENCE OF GRADIENT AND STATIONARY POINTS:

Here we analyze convergence properties of gradients of the empirical risk for deep nonlinear neural networks. Theorem 4. Assume the input sample x obeys Assumption 2 and the activation functions in a deep neural network are sigmoid functions. Then the empirical gradient uniformly converges to the population gradient in Euclidean norm. Specifically, there are two universal constants c y and c y such that if n ≥ c y c d l 3 r 2 /(s log(d)τ 2 ε 2 log(1/ε)) where c d = max 0≤i≤l d i , then with probability at least 1 -ε

sup w∈Ω ∇ Ĵn (w) -∇J (w) 2 ≤ l τ 512 729 c y l(l + 2) (lc r + 1) c d c r s log(dn/l) + log(4/ε) n ,

where c r = max(r 2 /16, r 2 /16 l-1 ), and s denotes the nonzero entry number of all weights.

Similar to deep linear neural networks, the layer number l, width d i , number of nonzero parameter entries s, network size d and magnitude of weights are all critical to the convergence rate. Also, since there is a factor max i d i in the convergence rate, it is better to avoid choosing an extremely wide layer. Interestingly, when analyzing the representation ability of deep learning, Eldan & Shamir (2016) also suggested non-extreme-wide layers, though the conclusion was derived from a different perspective. By comparing Theorems 1 and 4, one can observe that there is a factor (1/16) n and w (k) have the same non-degenerate index and they obey

w (k) n -w (k) 2 ≤ 2τ ζ 512 729 c y l(l + 2) (lc r + 1) c d c r s log(dn/l) + log(4/ε) n , (k = 1, • • • , m)

with probability at least 1 -ε, where c y , c d and c r are the same parameters in Theorem 4.

According to Theorem 5, there is one-to-one correspondence between the non-degenerate stationary points of Ĵn (w) and J (w). Also the corresponding pair has the same non-degenerate index, implying they have exactly matching local minima/maxima and non-degenerate saddle points. When n is sufficiently large, the non-degenerate stationary point w

(k)

n in Ĵn (w) is very close to its corresponding non-degenerate stationary point w (k) in J (w). As for the degenerate stationary points, Theorem 4 guarantees the gradients at these points of J (w) and Ĵn (w) are very close to each other.

Here we first give the uniform convergence analysis of the empirical risk and then analyze its stability and generalization. Theorem 6. Assume the input sample x obeys Assumption 2 and the activation functions in a deep neural network are the sigmoid functions. If n ≥ 18l 2 r 2 /(s log(d)τ 2 ε 2 log(1/ε)), then

sup w∈Ω Ĵn (w) -J (w) ≤ n τ 9 8 c y c d (1 + c r (l -1)) s log(nd/l) + log(4/ε) n (3)

holds with probability at least 1-ε, where c y , c d and c r are given in Theorem 4.

From Theorem 6, we obtain that under the distribution D, the empirical risk of a deep nonlinear neural network converges at the rate of O(1/ √ n) (up to a log factor). Theorem 6 also gives similar results as Theorem 3, including the inclination of regularization penalty on weight and suggestion on non-extreme-wide layers. Similar to linear networks, our risk convergence rate is also tighter than the convergence rate on the networks with polynomial activation functions and one-dimensional output in (Bartlett & Maass, 2003) since ours is at the order of O( (l -1)(s log(dn/l) + log(1/ε))/n), while the later is O( (γ log 2 (n) + log(1/ε))/n) where γ is at the order of O(ld log(d) + l 2 d) (Bartlett & Maass, 2003) .

We then establish the stability property and the generalization error of the empirical risk for nonlinear neural networks. By Theorem 6, we can obtain the following results. Corollary 2. Assume the input sample x obeys Assumption 2 and the activation functions in a deep neural network are sigmoid functions. Then with probability at least 1 -ε, we have

E S∼D,A,(x (j) ,y (j) )∼D 1 n n j=1 f * j -f j ≤ n and E S∼D,A J (w n ) -Ĵn (w n ) ≤ n ,

where n is defined in Eqn.

(3). The notations f * j and f j here are the same in Corollary 1.

By Corollary 2, we know that both the stability convergence rate and the convergence rate of generalization error are O(1/ √ n). This result accords with Theorems 8 and 9 in (Shalev-Shwartz et al., 2010) which implies O(1/ √ n) is the bottleneck of the stability and generalization convergence rate for generic learning algorithms. From this result, we have that if n is sufficiently large, the empirical risk can be expected to be very stable. This also dispels misgivings of the random selection of training samples in practice. Such a result indicates that the deep nonlinear neural network can offer good performance on testing data if it achieves small training error.

6. PROOF ROADMAP:

Here we briefly introduce our proof roadmap. Due to space limitation, all the proofs of Theorems 1 ∼ 6 and Corollaries 1 and 2 as well as technical lemmas are deferred to the supplementary material.

The proofs of Theorems 1 and 4 are similar but essentially differ in some techniques for bounding probability due to their different assumptions.

For explanation simplicity, we define four events:

E = {sup w∈Ω ∇ Ĵn (w) -∇J (w) 2 > t}, E 1 = {sup w∈Ω 1 n n i=1 ∇f (w, x (i) )-∇f (w kw , x (i) ) 2 > t/3}, E 2 = {sup w i kw ∈Ni, i∈[l] 1 n n i=1 ∇f (w kw , x (i) ) -E∇f (w kw , x) 2 > t/3}, and E 3 = {sup w∈Ω E∇f (w kw , x) -E∇f (w, x) 2 > t/3}, where w kw = [w 1 kw ; w 2 kw ; • • • ; w l kw ] is constructed by selecting w i kw ∈ R didi-1 from d i d i-1 /d-net N i such that w -w kw 2 ≤ .

Note that in Theorems 1 and 4, t is respectively set to g and l . Then we have P(E) ≤ P(E 1 ) + P(E 2 ) + P(E 3 ). So we only need to separately bound P(E 1 ), P(E 2 ) and P(E 3 ). For P(E 1 ) and P(E 3 ), we use the gradient Lipschitz constant and the properties of -net to prove P(E 1 ) ≤ ε/2 and P(E 3 ) = 0, while bounding P(E 2 ) needs more efforts. Here based on the assumptions, we prove that P(E 2 ) has sub-exponential tail associated to the sample number n and the networks parameters, and it satisfies P(E 2 ) ≤ ε/2 with proper conditions. Finally, combining the bounds of the three terms, we obtain the desired results.

To prove Theorems 2 and 5, we first prove the uniform convergence of the empirical Hessian to its population Hessian. Then, we define such a set D = {w ∈ Ω : ∇J (w) 2 < and inf i λ i ∇ 2 J (w) ≥ ζ}. In this way, D can be decomposed into countably components, with each component containing either exactly one or zero non-degenerate stationary point. For each component, the uniform convergence of gradient and the results in differential topology guarantee that if J (w) has no stationary points, then Ĵn (w) also has no stationary points and vise versa. Similarly, for each component, the uniform convergence of Hessian and the results in differential topology guarantee that if J (w) has a unique non-degenerate stationary point, then Ĵn (w) also has a unique non-degenerate stationary point with the same index. After establishing exact correspondence between the non-degenerate stationary points of empirical risk and population risk, we use the uniform convergence of gradient and Hessian to bound the distance between the corresponding pairs. We adopt a similar strategy to prove Theorems 3 and 6. Specifically, we divide the event sup w∈Ω | Ĵn (w)-∇J (w)| > t into E 1 , E 2 and E 3 which have the same forms as their counterparts in the proofs of Theorem 1 with the gradient replaced by the loss function. To prove P(E 1 ) ≤ ε/2 and P(E 3 ) = 0, we can use the Lipschitz constant of the loss function and the -net properties. The remaining is to prove P(E 2 ). We also prove that it has sub-exponential tail associated to the sample number n and the networks parameters and it obeys P(E 2 ) ≤ ε/2 with proper conditions. Then we utilize the uniform convergence of Ĵn (w) to prove the stability and generalization bounds of Ĵn (w) (i.e. Corollaries 1 and 2).

7. CONCLUSION:

In this work, we provided theoretical analysis on the landscape of empirical risk optimization for deep linear/nonlinear neural networks with (stochastic-)gradient descent, including the properties of the gradient and stationary points of empirical risk as well as the uniform convergence, stability, and generalization of the empirical risk itself. To our best knowledge, most of the results are new to deep learning community. These results also reveal that the depth l, the nonzero entry number s of all weights, the network size d and the width of a network are critical to the convergence rates. We also prove that the weight parameter magnitude is important to the convergence rate. Indeed, small magnitude of the weights is suggested. All the results are consistent with the widely used network architectures in practice.

SUPPLEMENTARY MATERIAL OF EMPIRICAL RISK LANDSCAPE ANALYSIS FOR UNDERSTANDING DEEP NEURAL NETWORKS A STRUCTURE OF THIS DOCUMENT:

This document gives some other necessary notations and preliminaries for our analysis in Sec. B. Then we prove Theorems 1∼ 3 and Corollary 1 for deep linear neural networks in Sec. C. Then we present the proofs of Theorems 4 ∼ 6 and Corollary 2 for deep nonlinear neural networks in Sec. D.

In both Sec. C and D, we first present the technical lemmas for proving our final results and subsequently present the proofs of these lemmas. Then we utilize these technical lemmas to prove our desired results. Finally, we give the proofs of other auxiliary lemmas.

B NOTATIONS AND PRELIMINARY TOOLS:

Beyond the notations introduced in the manuscript, we need some other notations used in this document. Then we introduce several lemmas that will be used later.

B.1 NOTATIONS:

Throughout this document, we use •, • to denote the inner product. A ⊗ C denotes the Kronecker product between A and C. Note that A and C in A ⊗ C can be matrices or vectors. For a matrix

A ∈ R n1×n2 , we use A F = i,j A 2 ij to denote its Frobenius norm, where A ij is the (i, j)-th entry of A. We use A op = max i |λ i (A)| to denote the operation norm of a matrix A ∈ R n1×n1

, where λ i (A) denotes the i-th eigenvalue of the matrix A. For a 3-way tensor A ∈ R n1×n2×n3 , its operation norm is computed as

A op = sup λ 2≤1 λ ⊗ 3 , A = i,j,k A ijk λ i λ j λ k ,

where A ijk denotes the (i, j, k)-th entry of A. Also we denote the vectorization of W (j) (the weight matrix of the j-th layer) as

w (j) = vec W (j) ∈ R dj dj-1 .

We denote I k as the identity matrix of size k × k.

For notational simplicity, we further define e v (l) -y as the output error vector. Then the squared loss is defined as f (w; x, y) = 1 2 e 2 2 , where w = (w (1) ; • • • ; w (l) ) ∈ R d contains all the weight parameters.

B.2 TECHNICAL LEMMAS:

We first introduce Lemmas 1 and 2 which are respectively used for bounding the 2 -norm of a vector and the operation norm of a matrix. Then we introduce Lemmas 3 and 4 which discuss the topology of functions. In Lemma 5, we give the relationship between the stability and generalization of empirical risk. Lemma 1. (Vershynin, 2012) For any vector x ∈ R d , its 2 -norm can be bounded as

x 2 ≤ 1 1 - sup λ∈λ λ, x .

where λ = {λ 1 , . . . , λ kw } be an -covering net of B d (1). Lemma 2. (Vershynin, 2012) For any symmetric matrix X ∈ R d×d , its operator norm can be bounded as

X op ≤ 1 1 -2 sup λ∈λ | λ, Xλ | .

where λ = {λ 1 , . . . , λ kw } be an -covering net of B d (1).

Lemma 3. (Mei et al., 2017) Let D ⊆ R d be a compact set with a C 2 boundary ∂D, and f, g : A → R be C 2 functions defined on an open set A, with D ⊆ A. Assume that for all w ∈ ∂D and all t ∈ [0, 1], t∇f (w) + (1 -t)∇g(w) = 0. Finally, assume that the Hessian ∇ 2 f (w) is non-degenerate and has index equal to r for all w ∈ D. Then the following properties hold:

(1) If g has no critical point in D, then f has no critical point in D.

(2) If g has a unique critical point w in D that is non-degenerate with an index of r, then f also has a unique critical point w in D with the index equal to r. Lemma 4. (Mei et al., 2017) Suppose that F (w) : Θ → R is a C 2 function where w ∈ Θ. Assume that {w (1) , . . . , w (m) } is its non-degenerate critical points and let D = {w ∈ Θ : ∇F (w) 2 < and inf i λ i ∇ 2 F (w) ≥ ζ}. Then D can be decomposed into (at most) countably components, with each component containing either exactly one critical point, or no critical point. Concretely, there exist disjoint open sets

{D k } k∈N , with D k possibly empty for k ≥ m + 1, such that D = ∪ ∞ k=1 D k . Furthermore, w (k) ∈ D k for 1 ≤ k ≤ m and each D i , k ≥ m + 1 contains no stationary points. Lemma 5. (Shalev-Shwartz & Ben-David, 2014b; Gonen & Shalev-Shwartz, 2017) Assume that D is a sample distribution and randomized algorithm A is employed for optimization. Suppose that ((x (1) , y (1) ), • • • , (x (n) , y (n) )) ∼ D and w n = argmin w Ĵn (w). For every j ∈ {1, • • • , n}, suppose w j * = argmin w 1 n-1 i =j f i (w; x (i) , y (i) ). For arbitrary distribution D, we have E S∼D,A,(x (j) ,y (j) )∼D 1 n n j=1 f * j -f j = E S∼D, A J (w n ) -Ĵn (w n ) .

where f * j and f j respectively denote f j (w j * ; x (j) , y (j) ) and f j (w n ; x (j) , y (j) ).

C PROOFS FOR DEEP LINEAR NEURAL NETWORKS:

In this section, we first present the technical lemmas in Sec. C.1 and then we give the proofs of these lemmas in Sec. C.2. Next, we utilize these lemmas to prove the results in Theorems 1∼ 3 and Corollary 1 in Sec. C.3. Finally, we give the proofs of other lemmas in Sec. C.4.

C.1 TECHNICAL LEMMAS:

Here we present the technical lemmas for proving our desired results. For brevity, we also define B j:s as follows:

EQUATION

Lemma 6. Assume that the activation functions in the deep neural network f (w, x) are linear functions. Then the gradient of f (w, x) with respect to w (j) can be written as

∇ w (j) f (w, x) = (B j-1:1 x) ⊗ B T l:j+1 e, (j = 1, • • • , l)

, where ⊗ denotes the Kronecke product. Then we can compute the Hessian matrix as follows:

∇ 2 f (w, x) =      ∇ w (1) ∇ w (1) f (w, x) • • • ∇ w (1) ∇ w (l) f (w, x) ∇ w (2) ∇ w (1) f (w, x) • • • ∇ w (2) ∇ w (l) f (w, x) . . . . . . . . . ∇ w (l) ∇ w (1) f (w, x) • • • ∇ w (l) ∇ w (l) f (w, x)      , where Q st ∇ w (s) ∇ w (t) f (w, x) is defined as Q st =      B T t-1:s+1 ⊗ B s-1:1 xe T B T l:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 , if s < t, B s-1:1 xx T B s-1:1 ⊗ B l:s+1 T B l:s+1 , if s = t, B T l:s+1 ex T B T t-1:1 ⊗B s-1:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 , if s > t.

Lemma 7. Suppose Assumption 1 on the input data x holds and the activation functions in deep neural network are linear functions. Then for any t > 0, the objective f (w, x) obeys

P 1 n n i=1 f (w, x (i) )-E(f (w, x (i) )) > t ≤ 2 exp   -c f n min   t 2 ω 2 f max d l ω 2 f τ 4 , τ 2 , t ω 2 f τ 2     ,

where c f is a positive constant and ω f = r l . Lemma 8. Suppose Assumption 1 on the input data x holds and the activation functions in deep neural network are linear functions. Then for any t > 0 and arbitrary unit vector λ ∈ S d-1 , the gradient ∇f (w, x) obeys

P 1 n n i=1 λ, ∇ w f (w, x (i) ) -E∇ w f (w, x (i) ) > t ≤ 3 exp -c g n min t 2 l max (ω g τ 2 , ω g τ 4 , ω g τ 2 ) , t lω g max (τ, τ 2 ) ,

where c g is a constant; ω g = c q r 2(2l-1) and ω g = c q r 2(l-1) in which c q = max 0≤i≤l d i .

Lemma 9. Suppose Assumption 1 on the input data x holds and the activation functions in deep neural network are linear functions. Then for any t > 0 and arbitrary unit vector λ ∈ S d-1 , the Hessian ∇ 2 f (w, x) obeys

P 1 n n i=1 λ, (∇ 2 w f (w, x (i) ) -E∇ 2 w f (w, x (i) ))λ > t ≤ 5 exp -c h n min t 2 τ 2 l 2 max (ω g , ω g τ 2 , ω h ) , t √ ω g l max (τ, τ 2 )

,

where ω g = r 4(l-1) and ω h = r 2(l-2) . Lemma 10. Suppose the activation functions in deep neural network are linear functions. Then for any w ∈ B d (r) and x ∈ B d0 (r x ), we have

∇ w f (w, x) 2 ≤ √ α g

, where α g = c t lr 4 x r 4l-2 . in which c t is a constant. Further, for any w ∈ B d (r) and x ∈ B d0 (r x ), we also have

∇ 2 f (w, x) op ≤ ∇ 2 f (w, x) F ≤ l √ α l

, where α l c t r 4 x r 4l-2 . in which c t is a constant. With the same condition, we can bound the operation norm of ∇ 3 f (w, x). That is, there exists a universal constant α p such that ∇ 3 f (w, x) op ≤ α p .

Lemma 11. Suppose Assumption 1 on the input data x holds and the activation functions in deep neural network are linear functions. Then there exist two universal constant c g and c h such that the sample Hessian converges uniformly to the population Hessian in operator norm. Specifically, there exit two universal constants c h1 and c h2 such that if n ≥ c h2 max(

α 2 p r 2 τ 2 l 2 ω 2 h ε 2 s log(d/l) , s log(d/l)/(lτ 2 )), then sup w∈Ω ∇ 2 Ĵn (w)-∇ 2 J (w) op ≤ c h1 τ lω h d log(nl)+log(20/ε) n

holds with probability at least 1 -ε, where ω h = max τ r 2(l-1) , r 2(l-2) , r l-2 .

C.2 PROOFS OF TECHNICAL LEMMAS:

To prove the above lemmas, we first introduce some useful results. Lemma 12. (Rudelson & Vershynin, 2013) Assume that x = (x 1 ; x 2 ; • • • ; x k ) ∈ R k is a random vector with independent components x i which have zero mean and are independent τ 2 i -sub-Gaussian variables. Here max i τ 2 i ≤ τ 2 . Let A be an k × k matrix. Then we have

E exp   λ   i,j:i =j A ij x i x j -E( i,j:i =j A ij x i x j )     ≤ exp 2τ 2 λ 2 A 2 F , |λ| ≤ 1/(2τ A 2 ).

Lemma 13. Assume that x = (x 1 ; x 2 ; • • • ; x k ) ∈ R k is a random vector with independent components x i which have zero mean and are independent τ 2 i -sub-Gaussian variables. Here max i τ 2 i ≤ τ 2 . Let a be an n-dimensional vector. Then we have

E exp λ k i=1 a i x 2 i -E k i=1 a i x 2 i ≤ E exp 128λ 2 τ 4 k i=1 a 2 i , |λ| ≤ 1 τ 2 max i a i .

Lemma 14. For B j:t defined in Eqn. (4), we have the following properties:

B s:t op ≤ B s:t F ≤ ω r and B l:1 op ≤ B l:1 F ≤ ω f ,

where ω r = r s-t+1 ≤ max r, r l and ω f = r l .

Lemma 13 is useful for bounding probability. The two inequalities in Lemma 14 can be obtained by using

w (j) 2 ≤ r (∀j = 1, • • • , l).

We defer the proofs of Lemmas 13 and 14 to Sec. C.4.2.

C.2.1 PROOF OF LEMMA 6:

Proof. When the activation functions are linear functions, we can easily compute the gradient of f (w, x) with respect to w (j) :

∇ w (j) f (w, x) = (B j-1:1 x) ⊗ B T l:j+1 e, (j = 1, • • • , l)

, where ⊗ denotes the Kronecker product. Now we consider the computation of the Hessian matrix. For brevity, let Q s = (B s-1:1 x) ⊗ B T l:s+1 . Then we can compute ∇ 2 w (s) f (w, x) as follows:

∇ 2 w (s) f (w, x) = ∂ 2 f (w, x) ∂w T (s) ∂w (s) = ∂ 2 f (w, x) ∂w T (s) ∂w (s) = ∂(Q s e) ∂w T (s) = ∂vec (Q s e) ∂w T (s) = ∂vec Q s B l:s+1 W (t) B s-1:1 x ∂w T (s) = ∂ (B s-1:1 x) T ⊗ (Q s B l:s+1 ) vec W (s) ∂w T (s) =(B s-1:1 x) T ⊗ (B s-1:1 x) ⊗ B T l:s+1 B l:s+1 x =(B s-1:1 x) T ⊗ (B s-1:1 x) ⊗ B T l:s+1 B l:s+1 y = (B s-1:1 x) T ⊗ (B s-1:1 x) ⊗ B T l:s+1 B l:s+1 z = (B s-1:1 x)(B s-1:1 x) T ⊗ B T l:s+1 B l:s+1

, where x holds since B j-1:1 x is a vector and for any vector x, we have (x ⊗ A)B = x ⊗ (AB). y holds because for any four matrices Z 1 ∼ Z 3 of proper sizes, we have

(Z 1 ⊗ Z 2 ) ⊗ Z 3 = Z 1 ⊗ (Z 2 ⊗ Z 3 ). z holds because for any two matrices z 1 , z 2 of proper sizes, we have z 1 z T 2 = z 1 ⊗ z T 2 = z T 2 ⊗ z 1 .

Then, we consider the case s > t:

∇ w (t) ∇ w (s) f (w, x) = ∂ 2 f (w, x) ∂w T (t) ∂w (s) = ∂ 2 f (w, x) ∂w T (t) ∂w (s) = ∂(Q s e) ∂w T (t) = ∂vec (Q s e) ∂w T (t) = ∂vec Q s B l:t+1 W (t) B t-1:1 x ∂w T (t) + ∂vec (B s-1:1 x) ⊗ B T l:s+1 e ∂w T (t)

.

Notice, here we just think that Q s in the

∂vec(QsB l:t+1 W (t) Bt-1:1x) ∂w T (t)

is a constant matrix and is not related to W (t) . Similarly, we also take e in ∂vec(((Bs-1:1x)⊗B T l:s+1 )e)

∂w T (t)

as a constant vector. Since we have

∂vec Q s B l:t+1 W (t) B t-1:1 x ∂w T (t) = B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 ,

we only need to consider

∂vec (B s-1:1 x) ⊗ B T l:s+1 e ∂w T (t) = ∂vec (B s-1:1 x) ⊗ B T l:s+1 e ∂w T (t) = ∂vec (B s-1:1 x) B T l:s+1 e T ∂w T (t) = ∂vec B s-1:t+1 W (t) B t-1:1 xe T B l:s+1 ∂w T t = ∂ B t-1:1 xe T B l:s+1 T ⊗ B s-1:t+1 vec W (t) ∂w T t = B t-1:1 xe T B l:s+1 T ⊗ B s-1:t+1 .

Therefore, for s > t, by combining the above two terms, we can obtain

∇ w (t) ∇ w (s) f (w, x) = B T l:s+1 ex T B T t-1:1 ⊗B s-1:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 .

Then, by similar method, we can compute the Hessian for the case s < t as follows:

∇ w (t) ∇ w (s) f (w, x) = B T t-1:s+1 ⊗ B s-1:1 xe T B T l:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 .

The proof is completed.

C.2.2 PROOF OF LEMMA 7:

Proof. We first prove that v (l) , which is defined in Eqn. ( 5), is sub-Gaussian.

EQUATION

Then by the convexity in λ of exp(λt) and Lemma 14, we can obtain

EQUATION

where x uses the conclusion that B l:1 op ≤ B l:1 F ≤ ω f in Lemma 14. This means that v (l) is centered and is ω 2 f τ 2 -sub-Gaussian. Accordingly, we can obtain that the k-th entry of v (l) is also

z k τ 2 -sub-Gaussian, where z k is a universal positive constant. Note that max k z k ≤ ω 2 f . Let v (l) i

denotes the output of the i-th sample x (i) . By Lemma 13, we have that for s > 0,

P 1 n n i=1 v (l) i 2 2 -E v (l) i 2 2 > t 2 = P s n i=1 v (l) i 2 2 -E v (l) i 2 2 > nst 2 x ≤ exp - snt 2 E s n i=1 v (l) i 2 2 -E v (l) i 2 2 y ≤ exp - snt 2 n i=1 E s v (l) i 2 2 -E v (l) i 2 2 z ≤ exp - snt 2 n i=1 exp 128d l s 2 ω 4 f τ 4 |s| ≤ 1 ω 2 f τ 2 { ≤ exp -c n min t 2 d l ω 4 f τ 4 , t ω 2 f τ 2 .

Note that x holds because of Chebyshev's inequality. y holds since x (i) are independent. z is established by applying Lemma 13. We have { by optimizing s. Since v (l) is sub-Gaussian, we have

P 1 n n i=1 y T v (l) i -Ey T v (l) i > t 2 ≤P s n i=1 y T v (l) i -Ey T v (l) i > nst 2 ≤ exp - nst 2 E exp s n i=1 y T v (l) i -Ey T v (l) i ≤ exp - nst 2 n i=1 E exp s y T v (l) i -Ey T v (l) i x ≤ exp - nst 2 n i=1 exp ω 2 f τ 2 s 2 y 2 2 2 y ≤ exp - nt 2 8ω 2 f τ 2 y 2 2 ,

where x holds because of Eqn. ( 6) and we have y since we optimize s.

Since the loss function f (w, x) is defined as l) . Therefore, we have

f (w, x) = v (l) -y 2 2 , we have f (w, x) -E(f (w, x)) = v (l) -y 2 2 -E( v (l) -y 2 2 ) = v (l) 2 2 -E v (l) 2 2 + y T v (l) -Ey T v (

P 1 n n i=1 f (w, x (i) ) -E(f (w, x (i) )) > t ≤P 1 n n i=1 v (l) i 2 2 -E v (l) i 2 2 > t 2 + P 1 n n i=1 y T v (l) i -Ey T v (l) i > t 2 ≤2 exp -c f n min t 2 d l ω 4 f τ 4 , t 2 ω 2 f τ 2 , t ω 2 f τ 2 .

where c f is a constant. Note that y 2 2 is the label of x, then it can also be bounded. The proof is completed.

C.2.3 PROOF OF LEMMA 8:

Proof. For brevity, let Q j denote ∇ w (j) f (w, x). Then, by Lemma 6 we have

EQUATION

where x holds since B j-1:1 x is a vector, and y holds because for any four matrices Z 1 ∼ Z 4 of proper sizes, we have

(Z 1 Z 3 )⊗(Z 2 Z 4 ) = (Z 1 ⊗Z 2 )(Z 3 ⊗Z 4 ).

Note that e = v (l) -y = B l:1 x-y. Then we know that the i-th entry Q i j has the form Q i j = p,q z ij pq x p x q + p y ij p x p + r ij (Step 1 blow will give the detailed analysis) where x p denotes the p-th entry in x. Note that z ij pq , y ij p and r ij are constants and independent on x.

We divide λ ∈ R l j=1 dj dj-1 into λ = (λ 1 ; • • • ; λ l ) where λ j ∈ R dj dj-1 . Let λ i j denote the i-th entry in λ j . Accordingly, we have

E λ, ∇ w f (w, x) -E∇ w f (w, x) = l j=1 λ j , Q j -EQ j = E 1 + E 2 + E 3 ,

where E 1 , E 2 , and E 3 are defined as

EQUATION

Thus, we can further separate the event as:

P (E > t) ≤ P 1 n n k=1 E k 1 > t 3 +P 1 n n k=1 E k 2 > t 3 +P 1 n n k=1 E k 3 > t 3 .

Thus, to prove our conclusion, we can respectively establish the upper bounds of the three events.

To the end, for each input sample x (i) , we divide its corresponding Q j -EQ j into E 1 , E 2 and E 3 . Then we bound the three events separately. Before that, we first give several equalities. Since B j:s = W (j) W (j-1) • • • W (s) (j ≥ s), by Lemma 14 we have

EQUATION

These two inequalities can be obtained by using

W (i) 2 F = w (i) 2 2 ≤ r 2 .

Step 1. Divide Q j -EQ j : Note that e = v (l) -y = B l:1 xy. Let H j = B j-1:1 ⊗ B T l:j+1 . Then we can further write Eqn. (7) as

EQUATION

)

where I d0 ∈ R d0×d0 is the identity matrix. According to Eqn. (10), we can write the i-th entry of Q j as the form Q i j = p,q z ij pq x p x q + p y ij p x p + r ij where x p denotes the p-th entry in x. Let Z j = H j (I d0 ⊗ B l:1 ) ∈ R dj dj-1×d 2 0 . Then, we know that the i-th entry

Q i j = Z(i, :)x , where x = x ⊗ x = [x 1 x; x 2 x; • • • , x d0 x] ∈ R d 2 0 .

In this way, we have

EQUATION

where c q = max 0≤i≤l d i .

We divide the i-th row H j (i, :

EQUATION

)

Let λ i j denote the i-th entry of λ j . Then, by Eqn. ( 8), we can obtain

j λ j , (Q j -E(Q j )) = p,q:p =q a pq (x p x q -Ex p x q )+ p a pp x 2 p -Ex 2 p + p b p (x p -Ex p ) = E 1 + E 2 + E 3 ,

where a pq and b p are defined as

a pq = l j=1 dj dj-1 i=1 λ i j z ij pq and b p = l j=1 dj dj-1 i=1 λ i j y ij p .

Note that for any four matrices of proper sizes, we have

EQUATION

Note that Eqn. ( 13) uses the conclusion in Eqn. ( 9). Therefore, we can have the following bound:

EQUATION

where x uses Eqn. ( 11). Then we can utilize Eqn. ( 14) and l j=1 dj dj-1 i=1

(λ i j ) 2 = 1 to bound a pq as follows:

a 2 pq ≤ l    l j=1   dj dj-1 i=1 λ i j z ij pq   2    ≤ l l j=1   dj dj-1 i=1 (λ i j ) 2     dj dj-1 i=1 (z ij pq ) 2   ≤ lω. which further gives p,q a 2 pq ≤ l l j=1   dj dj-1 i=1 (λ i j ) 2     dj dj-1 i=1 p,q (z ij pq ) 2   x ≤ lω.

where x uses Eqn. ( 14).

Similarly, we can obtain

EQUATION

)

So we can bound b p as

b 2 p ≤ l l j=1   dj dj-1 i=1 λ i j y ij p   2 ≤ l l j=1   dj dj-1 i=1 (λ i j ) 2     dj dj-1 i=1 (y ij p ) 2   ≤ lω ,

where ω = c q y 2 2 r 2(l-1) . Accordingly, we can have

p b 2 p ≤ l l j=1   dj dj-1 i=1 (λ i j ) 2     dj dj-1 i=1 p (y ij p ) 2   x ≤ lω ,

where x uses (15).

Step 2. Bound P(E 1 > t/3), P(E 2 > t/3) and P(E 3 > t/3): Let E k h1 denotes the E h1 which corresponds to the k-th sample x (k) . Therefore, we can bound

P 1 n n k=1 E k 1 > t 3 =P   s n k=1   p,q:p =q a k pq x k p x k q -Ex k p x k q   > snt 3   x ≤ exp - nst 3 E exp   s n k=1   p,q:p =q a k pq x k p x k q -Ex k p x k q     y ≤ exp - nst 3 n k=1 E exp   s   p,q:p =q a k pq x k p x k q -Ex k p x k q     z ≤ exp - nst 3 n k=1 exp   2τ 2 s 2 p,q:p =q (a k pq ) 2   |s| ≤ 1 2τ √ lω ≤ exp - nst 3 n j=1 exp 2τ 2 s 2 lω { ≤ exp -c n min t 2 ωlτ 2 , t √ lωτ ,

where x holds because of Chebyshev's inequality. y holds since x (i) are independent. z is established by applying Lemma 12. We have { by optimizing s. Similarly, by Lemma 13 we can bound P 1 n n k=1 E k 2 > t 3 as follows:

P 1 n n k=1 E k 2 > t 3 ≤ exp - nst 3 n k=1 E exp s p a k pp (x k p ) 2 -E(x k p ) 2 ≤ exp - nst 3 n k=1 exp 128τ 4 s 2 lω |s| ≤ 1 τ 2 √ lω ≤ exp -c n min t 2 ωlτ 4 , t √ lωτ 2 .

Finally, since x (i) are independent sub-Gaussian, we can use Hoeffding inequality and obtain

P 1 n n k=1 E k 3 > t 3 ≤ P 1 n n k=1 p b k p x k p -Ex k p > t 3 exp - c nt 2 ω lτ 2 .

Step 3. Bound P(E > t): By comparing the values of ω and ω , we can obtain

P (E > t) ≤P 1 n n k=1 E j 1 > t 3 +P 1 n n k=1 E j 2 > t 3 +P 1 n n k=1 E j 3 > t 3 ≤3 exp -c g n min t 2 l max (ω g τ 2 , ω g τ 4 , ω g τ 2 ) , t lω g max (τ, τ 2 )

,

where ω g = c q r 2(2l-1) and ω g = c q r 2(l-1) in which c q = max 0≤i≤l d i . The proof is completed.

C.2.4 PROOFS OF LEMMA 9:

Proof. For brevity, let Q ts denote ∇ w (t) ∇ w (s) f (w, x) . Then, by Lemma 6 we have

Q ts =      B T l:s+1 ex T B T t-1:1 ⊗ B s-1:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 , if s > t, B s-1:1 xx T B s-1:1 ⊗ B l:s+1 T B l:s+1 , if s = t, B T t-1:s+1 ⊗ B s-1:1 xe T B T l:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 , if s < t.

Then we know that the (i, k)-th entry Q ik ts has the form Q ik ts = p,q z ik pq x p x q + p y ik p x p + r ik (explained in the following Step 1. I) where x p denotes the p-th entry in x. Note that z ik pq , y ik p and r ik are constant and independent on x. For convenience, we let Q ts = H ts + G ts , where G ts = B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 and H ts is defined as

H ts =    B T l:s+1 ex T B T t-1:1 ⊗ B s-1:t+1 , if s > t, 0, if s = t, B T t-1:s+1 ⊗ B s-1:1 xe T B T l:t+1 , if s < t. Let E = 1 n n j=1 λ, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) λ , E h = 1 n n j=1 t,s λ t ,(H ts -E(H ts )) λ s , E g = 1 n n j=1 t,s λ t ,(G ts -E(G ts )) λ s .

Then we divide the event as two events:

P (E > t) = P (E h + E g > t) ≤ P (E h > t/2) + P (E g > t/2) .

Now we look each event separately. Similar to Q ts , the (i, k)-th entry H ik ts has the form H ik ts = p,q z ik pq x p x q + p y ik p x p + r ik . We divide the unit vector

λ ∈ R d as λ = (λ 1 ; • • • ; λ l ) where λ j ∈ R dj dj-1 . For input vector x, let t,s λ t , (H ts -E(H ts )) λ s = E h1 + E h2 + E h3

, where

EQUATION

where x p denotes the p-th entry in x and λ i j denotes the i-th entry of λ j . Let E j h1 , E j h2 , and E j h3

denote the E h1 , E h2 , and E j h3 of the j-th sample. Thus, considering n samples, we can further separately divide the two events above as:

P E h > t 2 ≤ P   1 n n j=1 E j h1 > t 6   +P   1 n n j=1 E j h2 > t 6   +P   1 n n j=1 E j h3 > t 6   .

Similarly, we can define E g1 , E g2 and E g3 .

P E g > t 2 ≤ P   1 n n j=1 E j g1 > t 6   +P   1 n n j=1 E j g2 > t 6   +P   1 n n j=1 E j g3 > t 6   .

Thus, to prove our conclusion, we can respectively establish the upper bounds of P E h > t 2 and P E g > t 2 .

Step 1:

Bound P E h > t 2

To achieve our goal, for each input sample x (i) , we divide its corresponding t,s (H ts -EH ts ) as E h1 , E h2 and E h3 . Then we bound the three events separately. Before that, we first give two equalities. Since B j:s = W (j) W (j-1) • • • W (s) (j ≥ s), by Lemma 14 we have

EQUATION

These two inequalities can be obtained by using

W (i) 2 F = w (i)

I. Divide H ts -EH ts : For t = s, we can write the (i, k)-th entry H ik ts as the form H ik ts = p,q z ik pq x p x q + p y ik p x p + r ik . Now we try to bound z ik pq and y ik p . We first consider the case s < t. Note that e = v (l) -y = B l:1 xy. Specifically, we have

H ts = B T t-1:s+1 ⊗ B s-1:1 xx T B T l:1 B T l:t+1 -B s-1:1 xy T B T l:t+1 .

So the (i , k )-th entry in the matrix B s-1:

EQUATION

where x uses Eqn. ( 17). Similarly, we can bound

EQUATION

where x uses Eqn. ( 17). Note that for the case s ≥ t, Eqn. ( 18) and ( 19) also holds. Let λ i j denote the i-th entry of λ j . Then, by Eqn. ( 16), we can obtain t,s

( λ t , (H ts -E(H ts ))λ s ) = p,q:p =q a pq (x p x q -Ex p x q )+ p a pp x 2 p -Ex 2 p + p b p (x p -Ex p ) = E h1 + E h2 + E h3 ,

where a pq and b p are defined as

a pq = t,s i,k (λ i t λ k s )z ik pq and b p = t,s i,k (λ i t λ k s )y ik p .

Then according to Eqn. (18) and t,s i,k (λ i t λ k s ) 2 = 1, we can bound a pq as follows:

a 2 pq ≤ l 2 t,s   i,k (λ i t λ k s )z ik pq   2 ≤ l 2 t,s   i,k (λ i t λ k s ) 2     i,k (z ik pq ) 2   ≤ ωl 2 t,s   i,k (λ i t λ k s ) 2   ≤ ωl 2 . which further yields p,q a 2 pq ≤ l 2 t,s   i,k (λ i t λ k s ) 2     i,k p,q (z ik pq ) 2   ≤ ωl 2 t,s   i,k (λ i t λ k s ) 2   ≤ ωl 2 .

Similarly, by using Eqn. ( 19), we have

b 2 p ≤ l 2 t,s   i,k (λ i t λ k s )y ik p   2 ≤ l 2 t,s   i,k (λ i t λ k s ) 2     i,k (y ik p ) 2   x ≤ ω l 2 .

Accordingly, we can have

p b 2 p ≤ l 2 t,s   i,k (λ i t λ k s ) 2     i,k p (y ik p ) 2   ≤ ω l 2 .

II. Bound P(E h1 > t/6), P(E h2 > t/6) and P(E h3 > t/6): Let E j h1 denotes the E j h1 which corresponds to the j-th sample x (i) . Therefore, we can bound

P   1 n n j=1 E j h1 > t 6   ≤P   s n j=1   p,q:p =q a j pq x j p x j q -Ex j p x j q   > snt 6   x ≤ exp - nst 6 E exp   s n j=1   p,q:p =q a j pq x j p x j q -Ex j p x j q     y ≤ exp - nst 6 n j=1 E exp   s   p,q:p =q a j pq x j p x j q -Ex j p x j q     z ≤ exp - nst 6 n j=1 exp   2τ 2 s 2 p,q:p =q (a j pq ) 2   |s| ≤ 1 2τ l √ ω ≤ exp - nst 6 n j=1 exp 2τ 2 s 2 l 2 ω { ≤ exp -c n min t 2 ωl 2 τ 2 , t √ ωlτ ,

where x holds because of Chebyshev's inequality. y holds since x (i) are independent. z is established because of Lemma 12. We have { by optimizing s. Similarly, we can bound

P 1 n n j=1 E j h2 > t 6 as follows: P   1 n n j=1 E j h2 > t 6   ≤ exp - nst 6 n j=1 E exp s p a j pp (x j p ) 2 -E(x j p ) 2 ≤ exp - nst 6 n j=1 exp 128τ 4 s 2 l 2 ω |s| ≤ 1 τ 2 l √ ω ≤ exp -c n min t 2 ωl 2 τ 4 , t √ ωlτ 2 .

Finally, since x (i) are independent sub-Gaussian, we can use Hoeffding inequality and obtain

P   1 n n j=1 E j h3 > t 6   = P   1 n n j=1 p b j p x j p -Ex j p > t 6   ≤ exp - c nt 2 ω l 2 τ 2 .

Since for s = t, P 1

n n j=1 E j h1 > t 6 = P 1 n n j=1 E j h2 > t 6 = P 1 n n j=1 E j h3 > t 6

= 0, the above upper bounds also hold.

III: Bound:

P E h > t

2 By comparing the values of ω and ω , we can obtain

P E h > t 2 ≤P   1 n n j=1 E j h1 > t 6   +P   1 n n j=1 E j h2 > t 6   +P   1 n n j=1 E j h3 > t 6   ≤3 exp -c 2 n min t 2 l 2 max (ωτ 2 , ωτ 4 , ω q τ 2 ) , t √ ωl max (τ, τ 2 ) ,

where ω q = r 2(l-2) .

Step 2: Bound P E g > t 2 To achieve our goal, for each input sample x (i) , we also divide its corresponding t,s (G ts -EG ts ) as E h1 , E h2 and E h3 . Then we bound the three events separately. Before that, we first give several equalities.

I. Divide G ts -EG ts : Dividing G ts -EG ts is more easy than dividing H ts -EH ts since the later has more complex form. Since G ts = B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 . we also can write the (i, k)-th entry G ik ts as the form G ik ts = p,q z ik pq x p x q + p y ik p x p + r ik . But here y ik p = 0. Then similar to the step in dividing H ts -EH ts , we can bound

a 2 pq ≤ ω g l 2 and

p,q a 2 pq ≤ ω g l 2 where ω g = r 4(l-1) .

II. Bound P(E g1 > t/6), P(E g2 > t/6) and P(E g3 > t/6): Since y ik p = 0, P(E h3 > t/6) = 0. Similar to the above methods, we can bound

P   1 n n j=1 E j g1 > t 6   ≤ exp -c 1 n t 2 ω g l 2 τ 2 , t √ ω g lτ ,

and

P   1 n n j=1 E j g2 > t 6   ≤ exp -c 1 n t 2 ω g l 2 τ 4 , t √ ω g lτ 2

.

III: Bound P E h > t 2 We can obtain P E g > t 2 as follows:

P E g > t 2 ≤P   1 n n j=1 E j g1 > t 6   +P   1 n n j=1 E j g2 > t 6   +P   1 n n j=1 E j g3 > t 6   ≤2 exp -c 2 n min t 2 ω g l 2 max (τ 2 , τ 4 ) , t √ ω g l max (τ, τ 2 )

.

Step 3: Bound P(E > t) Finally, we combine the above results and obtain

P (E > t) ≤P E h > t 2 +P E g > t 2 ≤5 exp -c h n min t 2 τ 2 l 2 max (ω g , ω g τ 2 , ω h ) , t √ ω g l max (τ, τ 2 )

,

where ω g = r 4(l-1) and ω h = r 2(l-2) .

C.2.5 PROOF OF LEMMA 10:

Proof. Before proving our conclusion, we first give an inequality:

e 2 2 = B l:1 x -y 2 2 ≤ B l:1 x 2 2 + 2 y T B l:1 x + y 2 2 x ≤ r 2 x ω 2 f + 2r x ω f y 2 + y

where ω f = r l . Notice, x holds since by Lemma 14, we have B l:1 2 F ≤ r 2l . Then we consider ∇ w f (w, x). Firstly, by Lemma 6 we can bound ∇ w (j) f (w, x) 2 2 as follows:

∇ w (j) f (w, x) 2 2 = (B j-1:1 x) ⊗ B T l:j+1 e 2 2 ≤ B j-1:1 2 2 x 2 2 B l:j+1 2 2 e 2 2 x ≤r 2 x ω 2 f1 r 2 x ω 2 f + 2r x ω f y 2 + y 2 2 ,

where ω f1 = r (l-1) . x holds since we have B l:j+1

2 F B j-1:1 2 F ≤ r 2(l-1) by using W (i) 2 F = w (i) 2 2 ≤ r 2 . Therefore, we can further obtain ∇ w f (w, x) 2 2 = l i=1 ∇ w (i) f (w, x) 2 2 ≤ lr 2 x ω 2 f1 r 2 x ω 2 f + 2r x ω f y 2 + y 2 2 .

Notice, y is the label of sample and the weight magnitude r is usually lager than 1. Then we have y 2 ≤ r l . Also, the values in input data are usually smaller than r l . Thus, we have

∇ w f (w, x) 2 2 ≤ c t lr 4

x r 4l-2 α g , where c t is a constant. Then we use the inequality ∇ 2 f (w, x) op ≤ ∇ 2 f (w, x) F to bound ∇ 2 f (w, x) op . Next we only need to give the upper bound of ∇ 2 f (w, x) F . Let ω f2 = r l-2 . We first consider Q st ∇ w (s) ∇ w (t) f (w, x) . By Lemma 6, if s < t, we have

Q st 2 F = B T t-1:s+1 ⊗ B s-1:1 xe T B T l:t+1 + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 2 F ≤2 B T t-1:s+1 ⊗ B s-1:1 xe T B T l:t+1 2 F + B s-1:1 xx T B T t-1:1 ⊗ B T l:s+1 B l:t+1 2 F ≤2 B t-1:s+1 2 F B s-1:1 2 F x 2 2 e 2 2 B l:t+1 2 F + 2 B s-1:1 2 F x 2 2 x 2 2 B t-1:1 2 F B l:s+1 2 F B l:t+1 2 F x ≤2ω 2 f2 r 2 x r 2 x ω 2 f + r x ω f y 2 + y 2 2 + 2ω 4 f1 r 4 x ,

where x holds since we use B l:t+1

2

F B t-1:s+1 2 F B s-1:1 2 F ≤ ω 2 f2 and B s-1:1 2 F B l:s+1 2 F

≤ ω 2 f1 . Note that when s ≥ t, the above inequality also holds. Similarly, consider the values in input data and the values in label, we have

Q st 2 F ≤ c t r 4

x r 4l-2 α l , where c t is a constant. Therefore, we can bound

∇ 2 f (w, x) op ≤ ∇ 2 f (w, x) F ≤ l s=1 l t=1 Q st 2 F ≤ l √ α l .

On the other hand, if the activation functions are linear functions, f (w, x) is fourth order differentiable when l ≥ 2. This means that ∇ x ∇ 3 w f (w, x) exists. Also since for any input x ∈ B d0 (r x ) and w ∈ Ω, we can always find a universal constant α p such that

∇ 3 w f (w, x) op = sup λ 2≤1 λ ⊗ 3 , ∇ 3 w f (w, x) = i,j,k [∇ 3 w f (w, x)] ijk λ i λ j λ k ≤ α p < +∞.

We complete the proofs.

C.2.6 PROOF OF LEMMA 11:

Proof. Recall that the weight of each layer has magnitude bound separately, i.e. w (j) 2 ≤ r.

Assume that w (j) has s j non-zero entries. Then we have l j=1 s j = s. So here we separately assume w j = {w j 1 , • • • , w j n j } is the d j d j-1 /d-covering net of the ball B dj dj-1 (r) which corresponds to the weight w (j) of the j-th layer. Let n j be the /l-covering number. By -covering theory in (Vershynin, 2012) , we can have

n j ≤ d j d j-1 s j 3r d j d j-1 /d sj ≤ exp s j log 3rd j d j-1 d j d j-1 /d = exp s j log 3rd .

Let w ∈ Ω be an arbitrary vector. Since w = [w (1) , • • • , w (l) ] where w (j) is the weight of the j-th layer, we can always find a vector w j kj in w j such that w (j) -w j kj 2 ≤ d j d j-1 /d. For brevity, let j w ∈ [n j ] denote the index of w j kj in -net w j . Then let

w kw = [w j k1 ; • • • ; w j kj ; • • • ; w j k l ]

. This means that we can always find a vector w kw such that ww kw 2 ≤ . Now we use the decomposition strategy to bound our goal:

∇ 2 Ĵn (w) -∇ 2 J (w) op = 1 n n i=1 ∇ 2 f (w, x (i) ) -E(∇ 2 f (w, x)) op = 1 n n i=1 ∇ 2 f (w, x (i) ) -∇f (w kw , x (i) ) + 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) + E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) op ≤ 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) op + 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) op + E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) op .

Here we also define four events E 0 , E 1 , E 2 and E 3 as

E 0 = sup w∈Ω ∇ 2 Ĵn (w) -∇ 2 J (w) op ≥ t , E 1 =    sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) op ≥ t 3    , E 2 =    sup jw∈[n j ],j=[l] 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) op ≥ t 3    , E 3 = sup w∈Ω E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) op ≥ t 3 .

Accordingly, we have

P (E 0 ) ≤ P (E 1 ) + P (E 2 ) + P (E 3 ) .

So we can respectively bound P (E 1 ), P (E 2 ) and P (E 3 ) to bound P (E 0 ).

Published as a conference paper at ICLR 2018:

Step 1. Bound P (E 1 ): We first bound P (E 1 ) as follows:

P (E 1 ) =P sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) 2 ≥ t 3 x ≤ 3 t E sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) 2 ≤ 3 t E sup w∈Ω ∇ 2 f (w, x) -∇ 2 f (w kw , x) 2 ≤ 3 t E sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) w -w kw 2 sup w∈Ω w -w kw 2 y ≤ 3α p t ,

where x holds since by Markov inequality and y holds because of Lemma 10.

Therefore, we can set t ≥ 6α p ε .

Then we can bound P(E 1 ):

P(E 1 ) ≤ ε 2 .

Step 2. Bound P (E 2 ): By Lemma 2, we know that for any matrix X ∈ R d×d , its operator norm can be computed as

X op ≤ 1 1 -2 sup λ∈λ | λ, Xλ | .

where λ = {λ 1 , . . . , λ kw } be an -covering net of B d (1).

Let λ 1/4 be the 1 4 -covering net of B d (1) but it has only s nonzero entries. So the size of its -net is d s

3 1/4 s ≤ exp (s log (12d)) .

Recall that we use j w to denote the index of w j kj in -net w j and we have j w ∈ [n j ], (n j ≤ exp s j log 3rd ). Then we can bound P (E 2 ) as follows:

P (E 2 ) =P sup jw∈[n j ] j∈[l] 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) 2 ≥ t 3 ≤P sup jw∈[n j ] j∈[l],λ∈λ 1/4 2 λ, 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E ∇ 2 f (w kw , x) λ ≥ t 3 ≤ exp (s log (12d)) exp   l j=1 s j log 3rd   sup jw∈[n j ] j∈[l],λ∈λ 1/4 P 1 n n i=1 λ, ∇ 2 f (w kw , x (i) )-E ∇ 2 f (w kw , x) ≥ t 6 x ≤ exp s log 36rd 2 10 exp -c h n min t 2 36τ 2 l 2 max (ω g , ω g τ 2 , ω h ) , t 6 √ ω g l max (τ, τ 2 ) ,

where x holds since by Lemma 9, we have

P 1 n n i=1 λ, (∇ 2 w f (w, x) -E∇ 2 w f (w, x))λ > t ≤ 10 exp -c h n min t 2 τ 2 l 2 max (ω g , ω g τ 2 , ω h ) , t √ ω g l max (τ, τ 2 )

,

where ω g = r 4(l-1) and ω h = r 2(l-2) . Let d = s log(36d 2 r/ )+log(20/ε). Thus, if we set

t ≥ max   36τ 2 l 2 max (ω g , ω g τ 2 , ω h ) d c h n , 6 √ ω g l max τ, τ 2 d c h n   ,

then we have

P (E 2 ) ≤ ε 2 .

Step 3. Bound P (E 3 ): We first bound P (E 3 ) as follows:

P (E 3 ) =P sup w∈Ω E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) 2 ≥ t 3 ≤P E sup w∈Ω (∇ 2 f (w kw , x) -∇ 2 f (w, x) 2 ≥ t 3 ≤P sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) w -w kw 2 sup w∈Ω w -w kw 2 ≥ t 3 x ≤P α p ≥ t 3 ,

where x holds because of Lemma 10. We set enough small such that α p < t/3 always holds.

Then it yields P (E 3 ) = 0.

Step 4. Final result: For brevity, let ω 2 = 36τ 2 l 2 max ω g , ω g τ 2 , ω h and ω 3 = 6 √ ω g l max τ, τ 2 .

To ensure P(E 0 ) ≤ ε, we just set = 36rl/n and

t ≥ max   6α p ε , 3α p , ω 2 (s log(36d 2 r/ )+log(20/ε)) c h n , ω 3 (s log(36d 2 r/ )+log(20/ε)) c h n   = max   216α p r nε , ω 2 (s log(d 2 nl)+log(20/ε)) c h n , ω 3 (s log(36d 2 n/l)+log(20/ε)) c h n   .

Thus, there exit two universal constants c h1 and c h2 such that if n ≥ c h2 max(

α 2 p r 2 τ 2 l 2 ω 2 h ε 2 s log(d/l) , s log(d/l)/(lτ 2 )), then sup w∈Ω ∇ 2 Ĵn (w)-∇ 2 J (w) op ≤ c h1 τ lω h d log(nl)+log(20/ε) n

holds with probability at least 1 -ε, where ω h = max τ r 2(l-1) , r 2(l-2) , r l-2 . The proof is completed.

C.3 PROOFS OF MAIN THEOREMS C.3.1 PROOF OF THEOREM 1:

Proof. Recall that the weight of each layer has magnitude bound separately, i.e. w (j) 2 ≤ r.

Assume that w (j) has s j non-zero entries. Then we have l j=1 s j = s. So here we separately assume w j = {w j 1 , • • • , w j n j } is the d j d j-1 /d-covering net of the ball B dj dj-1 (r) which corresponds to the weight w (j) of the j-th layer. Let n j be the /l-covering number. By -covering theory in (Vershynin, 2012) , we can have

n j ≤ d j d j-1 s j 3r d j d j-1 /d sj ≤ exp s j log 3rd j d j-1 d j d j-1 /d = exp s j log 3rd .

Let w ∈ Ω be an arbitrary vector. Since w = [w (1) , • • • , w (l) ] where w (j) is the weight of the j-th layer, we can always find a vector w j kj in w j such that w (j) -w j kj 2 ≤ d j d j-1 /d. For brevity, let j w ∈ [n j ] denote the index of w j kj in -net w j . Then let

w kw = [w j k1 ; • • • ; w j kj ; • • • ; w j k l ]

. This means that we can always find a vector w kw such that ww kw 2 ≤ . Accordingly, we can decompose ∇ Ĵn (w) -∇J (w)

2 as ∇ Ĵn (w) -∇J (w) 2 = 1 n n i=1 ∇f (w, x (i) ) -E(∇f (w, x)) 2 = 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) + 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) + E(∇f (w kw , x)) -E(∇f (w, x)) 2 ≤ 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 + 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) 2 + E(∇f (w kw , x)) -E(∇f (w, x)) 2 .

Here we also define four events E 0 , E 1 , E 2 and E 3 as

E 0 = sup w∈Ω ∇ Ĵn (w) -∇J (w) 2 ≥ t , E 1 = sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 ≥ t 3 , E 2 = sup jw∈[n j ],j=[l] 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) 2 ≥ t 3 , E 3 = sup w∈Ω E(∇f (w kw , x)) -E(∇f (w, x)) 2 ≥ t 3 .

Accordingly, we have P (E 0 ) ≤ P (E 1 ) + P (E 2 ) + P (E 3 ) . So we can respectively bound P (E 1 ), P (E 2 ) and P (E 3 ) to bound P (E 0 ).

Step 1. Bound P (E 1 ): We first bound P (E 1 ) as follows:

P (E 1 ) =P sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 ≥ t 3 x ≤ 3 t E sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 ≤ 3 t E sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 w -w kw 2 sup w∈Ω w -w kw 2 ≤ 3 t E sup w∈Ω ∇ 2 Ĵn (w, x) 2 ,

where x holds since by Markov inequality, we have that for an arbitrary nonnegative random variable x, then P(x ≥ t) ≤ E(x) t .

Published as a conference paper at ICLR 2018 Now we only need to bound E sup w∈Ω ∇ 2 Ĵn (w, x)

2:

. Now we utilize Lemma 10 to achieve this goal:

E sup w∈Ω ∇ 2 Ĵn (w, x) 2 ≤= E sup w∈Ω ∇ 2 f (w, x) -∇ 2 f (w * , x) 2 ≤ l √ α l .

where α l = c t r 4 x r 4l-2 . Therefore, we have

P (E 1 ) ≤ 3l √ α l t .

We further let

t ≥ 6l √ α l ε .

Then we can bound P(E 1 ):

P(E 1 ) ≤ ε 2 .

Step 2. Bound P (E 2 ): By Lemma 1, we know that for any vector x ∈ R d , its 2 -norm can be computed as

x 2 ≤ 1 1 - sup λ∈λ λ, x .

where λ = {λ 1 , . . . , λ kw } be an -covering net of B d (1).

Let λ 1/2 be the 1 2 -covering net of B d (1) but it has only s nonzero entries. So the size of its -net is d s

3 1/2 s ≤ exp (s log (6d)) .

Recall that we use j w to denote the index of w j kj in -net w j and we have j w ∈ [n j ], (n j ≤ exp s j log 3rd ). Then we can bound P (E 2 ) as follows:

P (E 2 ) =P sup jw∈[n j ],j=[l] 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) 2 ≥ t 3 =P sup jw∈[n j ],j=[l],λ∈λ 1/2 2 λ, 1 n n i=1 ∇f (w kw , x (i) ) -E (∇f (w kw , x)) ≥ t 3 ≤ exp (s log (6d)) exp   l j=1 s j log 3rd   sup jw∈[n j ],j=[l],λ∈λ 1/2 P 1 n n i=1 λ, ∇f (w kw , x (i) ) -E (∇f (w kw , x)) ≥ t 6 x ≤ exp s log 18rd 6 exp -c g n min t 2 36l max (ω g τ 2 , ω g τ 4 , ω g τ 2 ) , t 6 lω g max (τ, τ 2 ) ,

where x holds since by Lemma 8, we have

P 1 n n i=1 λ, ∇ w f (w, x (i) ) -E∇ w f (w, x (i) ) > t ≤ 3 exp -c g n min t 2 l max (ω g τ 2 , ω g τ 4 , ω g τ 2 ) , t lω g max (τ, τ 2 ) ,

where c g is a constant; ω g = c q r 2(2l-1) and ω g = c q r 2(l-1) in which c q = max 0≤i≤l d i .

Let ω 2 = 36l max ω g τ 2 , ω g τ 4 , ω g τ 2 and ω 3 = 6 lω g max τ, τ 2 . Thus, if we set t ≥ max ω 2 (s log(18dr/ )+log(12/ε)) c g n , ω 3 (s log(18dr/ )+log(12/ε)) c g n , then we have

P (E 2 ) ≤ ε 2 .

Step 3. Bound P (E 3 ): We first bound P (E 3 ) as follows:

P (E 3 ) =P sup w∈Ω E(f (w kw , x)) -E(f (w, x)) 2 ≥ t 3 =P sup w∈Ω E (f (w kw , x) -f (w, x) 2 ) w -w kw 2 sup w∈Ω w -w kw 2 ≥ t 3 ≤P E sup w∈Ω ∇ 2 Ĵn (w, x) 2 ≥ t 3 ≤P l √ α l ≥ t 3 .

We set enough small such that l √ α l < t/3 always holds. Then it yields P (E 3 ) = 0.

Step 4. Final result: Finally, to ensure P(E 0 ) ≤ ε, we just set = 18lr/n and

t ≥ max 6l √ α l ε , 3l √ α l , ω 2 (s log(18dr/ )+log(12/ε)) c g n , ω 3 (s log(18dr/ )+log(12/ε)) c g n = max 108l 2 √ α l r nε , ω 2 (s log(dn/l)+log(12/ε)) c g n , ω 3 (s log(dn/l)+log(12/ε)) c g n .

Notice, we have α l = c t r 4 x r 4l-2 where c t is a constant. Therefore, there exists two universal constants c g and c g such that n ≥ c g max(

l 3 r 2 r 4 x cqs log(d/l)ε 2 τ 4 log(1/ε) , s log(d/l)/(lτ 2 )), then sup w∈Ω ∇ Ĵn (w)-∇J (w) 2 ≤ c g τ ω g lc q

s log(dn/l)+log(12/ε) n

holds with probability at least 1 -ε, where ω g = max τ r 2l-1 , r 2l-1 , r l-1 .

C.3.2 PROOF OF THEOREM 2:

Proof. Suppose that {w (1) , w (2) , • • • , w (m) } are the non-degenerate critical points of J (w). So for any w (k) , it obeys

inf i λ k i ∇ 2 J (w (k) ) ≥ ζ,

where

λ k i ∇ 2 J (w (k) ) denotes the i-th eigenvalue of the Hessian ∇ 2 J (w (k) ) and ζ is a constant. We further define a set D = {w ∈ R d | ∇J (w) 2 ≤ and inf i |λ i ∇ 2 J (w (k) ) | ≥ ζ}. According to Lemma 4, D = ∪ ∞ k=1 D k where each D k is a disjoint component with w (k) ∈ D k for k ≤ m

and D k does not contain any critical point of J (w) for k ≥ m + 1. On the other hand, by the continuity of ∇J (w), it yields ∇J (w) 2 = for w ∈ ∂D k . Notice, we set the value of blow which is actually a function related to n.

Then by utilizing Theorem 1, we let sample number n sufficient large such that

sup w∈Ω ∇ Ĵn (w) -∇J (w) 2 ≤ z g 2

holds with probability at least 1 -ε, where if n ≥ c g max(

l 3 r 2 r 4 x cqs log(d/l)ε 2 τ 4 log(1/ε) , s log(d/l) lτ 2

), z g = c g τ ω g lc q s log(dn/l)+log(12/ε) n

.

This further gives that for arbitrary w ∈ D k , we have

EQUATION

)

Similarly, by utilizing Lemma 11, let n be sufficient large such that

sup w∈Ω ∇ 2 Ĵn (w) -∇ 2 J (w) op ≤ z s ≤ ζ 2

holds with probability at least 1 -ε, where if n ≥ c h2 max(

α 2 p r 2 τ 2 l 2 ω 2 h ε 2 s log(d/l) , s log(d/l)/(lτ 2 )), z s = c h1 τ lω h s log(nl)+log(20/ε) n .

Assume that b ∈ R d is a vector and satisfies b T b = 1. In this case, we can bound

λ k i ∇ 2 Ĵn (w)

for arbitrary w ∈ D k as follows:

inf w∈D k λ k i ∇ 2 Ĵn (w) = inf w∈D k min b T b=1 b T ∇ 2 Ĵn (w)b = inf w∈D k min b T b=1 b T ∇ 2 Ĵn (w) -∇ 2 J (w) b + b T ∇ 2 J (w)b ≥ inf w∈D k min b T b=1 b T ∇ 2 J (w)b -min b T b=1 b T ∇ 2 Ĵn (w) -∇ 2 J (w) b ≥ inf w∈D k min b T b=1 b T ∇ 2 J (w)b -max b T b=1 b T ∇ 2 Ĵn (w) -∇ 2 J (w) b = inf w∈D k inf i |λ k i ∇ 2 f (w (k) , x) | -∇ 2 Ĵn (w) -∇ 2 J (w) op ≥ ζ 2 .

This means that in each set D k , ∇ 2 Ĵn (w) has no zero eigenvalues. Then, combine this and Eqn. ( 20), by Lemma 3 we know that if the population risk J (w) has no critical point in D k , then the empirical risk Ĵn (w) has also no critical point in D k ; otherwise it also holds. By Lemma 3, we can also obtain that in D k , if J (w) has a unique critical point w (k) with non-degenerate index s k , then Ĵn (w) also has a unique critical point w n (k) in D k with the same non-degenerate index s k . The first conclusion is proved. Now we bound the distance between the corresponding critical points of J (w) and Ĵn (w). Assume that in D k , J (w) has a unique critical point w (k) and Ĵn (w) also has a unique critical point w (k) n . Then, there exists t ∈ [0, 1] such that for any z ∈ ∂B d (1), we have

≥ ∇J (w (k) n ) 2 = max z T z=1 ∇J (w (k) n ), z = max z T z=1 ∇J (w (k) ), z + ∇ 2 J (w (k) + t(w (k) n -w (k) ))(w (k) n -w (k) ), z x ≥ ∇ 2 J (w (k) ) 2 (w (k) n -w (k) ), (w (k) n -w (k) ) 1/2 y ≥ζ w (k) n -w (k)

2 , where x holds since ∇J (w (k) ) = 0 and y holds since w (k) + t(w

(k) n -w (k) ) is in D k and for any w ∈ D k we have inf i |λ i ∇ 2 J (w) | ≥ ζ.

Consider the conditions in Lemma 11 and Theorem 1, we can obtain that if n ≥ c h max(

l 3 r 2 r 4 x cqs log(d/l)ε 2 τ 4 log(1/ε) , s log(d/l)/ζ 2 ) where c h is a constant, then w (k) n -w (k) 2 ≤ 2c g τ ω g ζ lc q

s log(dn/l)+log(12/ε) n

holds with probability at least 1 -ε.

C.3.3 PROOF OF THEOREM 3:

Proof. Recall that the weight of each layer has magnitude bound separately, i.e. w (j) 2 ≤ r.

Assume that w (j) has s j non-zero entries. Then we have l j=1 s j = s. So here we separately assume w j = {w j 1 , • • • , w j n j } is the d j d j-1 /d-covering net of the ball B dj dj-1 (r) which corresponds to the weight w (j) of the j-th layer. Let n j be the /l-covering number. By -covering theory in (Vershynin, 2012) , we can have

n j ≤ d j d j-1 s j 3r d j d j-1 /d sj ≤ exp s j log 3rd j d j-1 d j d j-1 /d = exp s j log 3rd .

Let w ∈ Ω be an arbitrary vector. Since w = [w (1) , • • • , w (l) ] where w (j) is the weight of the j-th layer, we can always find a vector w j kj in w j such that w (j) -w j kj 2 ≤ d j d j-1 /d. For brevity, let j w ∈ [n j ] denote the index of w j kj in -net w j . Then let

w kw = [w j k1 ; • • • ; w j kj ; • • • ; w j k l ]

. This means that we can always find a vector w kw such that ww kw 2 ≤ . Now we use the decomposition strategy to bound our goal:

Ĵn (w) -J (w) = 1 n n i=1 f (w, x (i) ) -E(f (w, x)) = 1 n n i=1 f (w, x (i) )-f (w kw , x (i) ) + 1 n n i=1 f (w kw , x (i) )-Ef (w kw , x)+Ef (w kw , x)-Ef (w, x) ≤ 1 n n i=1 f (w, x (i) )-f (w kw , x (i) ) + 1 n n i=1 f (w kw , x (i) )-Ef (w kw , x) + Ef (w kw , x)-Ef (w, x) .

Then, we define four events E 0 , E 1 , E 2 and E 3 as

E 0 = sup w∈Ω Ĵn (w) -J (w) ≥ t , E 1 = sup w∈Ω 1 n n i=1 f (w, x (i) ) -f (w kw , x (i) ) ≥ t 3 , E 2 = sup jw∈[n j ],j=[l] 1 n n i=1 f (w kw , x (i) )-E(f (w kw , x)) ≥ t 3 , E 3 = sup w∈Ω E(f (w kw , x))-E(f (w, x)) ≥ t 3 .

Accordingly, we have P (E 0 ) ≤ P (E 1 ) + P (E 2 ) + P (E 3 ) . So we can respectively bound P (E 1 ), P (E 2 ) and P (E 3 ) to bound P (E 0 ).

Step 1. Bound P (E 1 ): We first bound P (E 1 ) as follows:

P (E 1 ) =P sup w∈Ω 1 n n i=1 f (w, x (i) ) -f (w kw , x (i) ) ≥ t 3 x ≤ 3 t E sup w∈Ω 1 n n i=1 f (w, x (i) ) -f (w kw , x (i) ) ≤ 3 t E sup w∈Ω 1 n n i=1 f (w, x (i) ) -f (w kw , x (i) ) w -w kw 2 sup w∈Ω w -w kw 2 ≤ 3 t E sup w∈Ω ∇ Ĵn (w, x) 2 ,

where x holds since by Markov inequality, we have that for an arbitrary nonnegative random variable x, then

P(x ≥ t) ≤ E(x) t .

Now we only need to bound E sup w∈Ω ∇ Ĵn (w, x)

. Therefore, by Lemma 10, we have

E sup w∈Ω ∇ Ĵn (w, x) 2 = E sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) 2 = E sup w∈Ω ∇f (w, x) 2 ≤ √ α g .

where α g = c t lr 4 x r 4l-2 . Therefore, we have

P (E 1 ) ≤ 3 √ α g t .

We further let

t ≥ 6 √ α g ε .

Then we can bound P(E 1 ):

P(E 1 ) ≤ ε 2 .

Step 2. Bound P (E 2 ): Recall that we use j w to denote the index of w j kj in -net w j and we have j w ∈ [n j ], (n j ≤ exp s j log 3rd ). We can bound P (E 2 ) as follows:

P (E 2 ) =P sup jw∈[n j ] j∈[l] 1 n n i=1 f (w kw , x (i) ) -E(f (w kw , x)) ≥ t 3 ≤ exp   l j=1 s j log 3rd   sup jw∈[n j ] j∈[l] P 1 n n i=1 f (w kw , x (i) ) -E(f (w kw , x)) ≥ t 3 x ≤4 3dr s exp   -c f n min   t 2 9ω 2 f max d l ω 2 f τ 4 , τ 2 , t 3ω 2 f τ 2     ,

where x holds because in Lemma 7, we have

P 1 n n i=1 f (w, x (i) )-E(f (w, x (i) )) > t ≤ 2 exp   -c f n min   t 2 ω 2 f max d l ω 2 f τ 4 , τ 2 , t ω 2 f τ 2     ,

where c f is a positive constant and ω f = r l . Thus, if we set

t ≥ max    9ω 2 f (s log(3rd/ ) + log(8/ε)) max d l ω 2 f τ 4 , τ 2 c f n , 3ω 2 f τ 2 (s log(3rd/ ) + log(8/ε)) c f n    ,

then we have

P (E 2 ) ≤ ε 2 .

Step 3. Bound P (E 3 ): We first bound P (E 3 ) as follows:

P (E 3 ) =P sup w∈Ω E(f (w kw , x)) -E(f (w, x)) 2 ≥ t 3 =P sup w∈Ω E (f (w kw , x) -f (w, x) 2 ) w -w kw 2 sup w∈Ω w -w kw 2 ≥ t 3 ≤P E sup w∈Ω ∇J w (w, x) 2 ≥ t 3 x ≤P √ α g ≥ t 3 ,

where x holds since we utilize Lemma 10. We set enough small such that √ α g < t/3 always holds. Then it yields P (E 3 ) = 0.

Step 4. Final result: To ensure P(E 0 ) ≤ ε, we just set = 3rl/n. Note that

6 √ αg ε > 3 √ α g . Thus we can obtain t ≥ max    6 √ α g ε , 9ω 2 f (s log(3rd/ )+log(8/ε)) max d l ω 2 f τ 4 , τ 2 c f n , 3ω 2 f τ 2 (s log(3rd/ )+log(8/ε)) c f n    = max    18l √ α g r nε , 9ω 2 f (s log(dn/l)+log(8/ε)) max d l ω 2 f τ 4 , τ 2 c f n , 3ω 2 f τ 2 (s log(dn/l)+log(8/ε)) c f n    .

Note that we have α g = c t lr 4 x r 4l-2 where c t is a constant. Then Then there exist four universal constants c f and c f such that if n ≥ c f max

l 3 r 4 x d l s log(d)ε 2 τ 4 log(1/ε) , s log(d)/(τ 2 d l ) , then sup w∈Ω Ĵn (w) -J (w) 2 ≤ c f ω f τ max d l ω f τ, 1 s log(dn/l) + log(8/ε) n

holds with probability at least 1 -ε.

C.3.4 PROOF OF COROLLARY 1:

Proof. By Lemma 5, we know s = g . Thus, the remaining work is to bound s . Actually, we can have

E S∼D,A,(x (j) ,y (j) )∼D 1 n n j=1

f j (w j * ;x (j) ,y (j) )-f j (w n ;x (j) ,y (j) ) ≤E S∼D sup w∈Ω Ĵn (w) -J (w)

≤ sup w∈Ω Ĵn (w) -J (w) ≤ f .

Thus, we have g = s ≤ f . The proof is completed.

C.4 PROOF OF OTHER LEMMAS C.4.1 PROOF OF LEMMA 13

Lemma 15. (Rigollet, 2015) Suppose a random variable x is τ 2 -sub-Gaussian, then the random variable x 2 -Ex 2 is sub-exponential and obeys:

EQUATION

)

Proof. Here we utilize Lemma 15 to prove our conclusion. We have

E exp λ k i=1 a i x 2 i -E k i=1 a i x 2 i x = k i=1 E exp λa i x 2 i -Ex 2 i y ≤ k i=1 E exp 128λ 2 a 2 i τ 4 i , |λ| ≤ 1 max i a i τ 2 ≤E exp 128λ 2 τ 4 k i=1 a 2 i

, where x holds since x i are independent and y holds because of Lemma 15.

C.4.2 PROOF OF LEMMA 14:

Proof. Since the 2 -norm of each w (j) is bounded, i.e. w (j) 2 ≤ r (1 ≤ j ≤ l), we can obtain

B s:t 2 F ≤ W (s) 2 F W (s-1) 2 F • • • W (t) 2 F ≤ r 2(t-s+1) ω 2 r x ≤ max r 2 , r 2l ,

where x holds since the function r 2x obtains its maximum at two endpoints x = 1 and x = l for case r < 1 and r ≥ 1, respectively. On the other hand, we have B s:t op ≤ B s:t F ≤ ω r . Specifically,

we have B l:1 2 F ≤ r 2l ω 2 f .

D PROOFS FOR DEEP NONLINEAR NEURAL NETWORKS:

In this section, we first present the technical lemmas in Sec. D.1. Then in Sec. D.2 we give the proofs of these lemmas. Next, we utilize these technical lemmas to prove the results in Theorems 4 ∼ 6 and Corollary 2 in Sec. D.3. Finally, we give the proofs of other lemmas in Sec. D.4.

D.1 TECHNICAL LEMMAS:

Here we present the key lemmas and theorems for proving our desired results. For brevity, we define an operation G which maps an arbitrary vector z ∈ R k into a diagonal matrix G(z) ∈ R k×k with its i-th diagonal entry equal to σ(z i )(1 -σ(z i )) in which z i denotes the i-th entry of z. We further define A i ∈ R di-1×di as follows:

EQUATION

where W (i) is the weight matrix in the i-th layer and u (i) is the linear output of the i-th layer. In this section, we define

EQUATION

Lemma 16. Suppose that the activation function in deep neural network are sigmoid functions. Then the gradient of f (w, x) with respect to w (j) can be formulated as

∇ w (j) f (w, x) = vec G(u (j) )B j+1:l (v (l) -y) (v (j-1) ) T , (j = 1, • • • , l -1),

and

∇ w (l) f (w, x) = vec G(u (l) )(v (l) -y) (v (l-1) ) T .

Besides, the loss f (w, x) is α-Lipschitz,

∇ w f (w, x) 2 ≤ α, where α = 1 16 c y c d (1 + c r (l -1)

) in which c y , c d and c r are defined as

v (l) -y 2 2 ≤ c y < +∞, c d = max(d 0 , d 1 , • • • , d l ) and c r = max r 2 16 , r 2 16 l-1 .

Lemma 17. Suppose that the activation functions in deep neural network are sigmoid functions.

Then there exists two universal constants c s1 and c s2 such that

∇ 2 w f (w, x) op ≤ ∇ 2 w f (w, x) F ≤ ς, where ς = c s1 c r c 2 d l 4 in which c d = max i d i and c r = max r 2 16 , r 2 16 l-1 . Moreover, the gradient ∇ w f (w, x) is ς-Lipschitz, i.e. ∇ w f (w 1 , x) -∇ w f (w 2 , x) 2 ≤ς w 1 -w 2 2 .

Similarly, there also exist a universal constant ξ such that

∇ 3 w f (w, x) op ≤ ∇ 3 w f (w, x) F ≤ ξ.

Lemma 18. Suppose that the activation function in deep neural network are sigmoid functions. Then we have

∇ w ∇ x f (w, x) op ≤ ∇ w ∇ x f (w, x) F ≤β,

where β = 2 6 3 8 l(l + 2)c y c r c d (lc r + 1) in which c y , c d and c r are defined in Lemma 16. Lemma 19. Suppose that the input sample x obeys Assumption 2 and the activation functions in deep neural network are sigmoid functions. The gradient of the loss is 8β 2 τ 2 -sub-Gaussian. Specifically, for any λ ∈ R d , we have

E ( λ, ∇ w f (w, x) -E∇ w f (w, x) ) ≤ exp 8β 2 τ 2 λ 2 2 2 ,

where β = 2 6 3 8 l(l + 2)c y c r c d (lc r + 1) in which c y , c d and c r are defined in Lemma 16. Lemma 20. Suppose that the input sample x obeys Assumption 2 and the activation functions in deep neural network are sigmoid functions. The Hessian of the loss, evaluated on a unit vector, is sub-Gaussian. Specifically, for any unit λ ∈ S d-1 (i.e. λ 2 = 1), there exist universal constant γ such that

E t λ, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) λ ≤ exp 8t 2 γ 2 τ 2 2 .

Notice, γ obeys γ ≥ ∇ x ∇ 2 w f (w, x) op . Lemma 21. Assume that the input sample x obeys Assumption 2 and the activation functions in deep neural network are sigmoid functions. Then the sample Hessian uniformly converges to the population Hessian in operator norm. That is, there exists such two universal constants c m and c

m such that if n ≥ c m ξ 2 l 2 r 2 γ 2 τ 2 ε 2 s log(d) log(1/ε) , then sup w∈Ω ∇ 2

Ĵn (w)-∇ 2 J (w) op ≤ c m γτ s log(dn/l)+log(4/ε) n

holds with probability at least 1 -ε. Here γ is the same parameter in Lemma 20.

D.2 PROOFS OF TECHNICAL LEMMAS:

For brevity, we also define

D s:t = W (s) 2 F • • • W (t) 2 F (s ≤ t) and D s:t = 1, (s > t).

We define a matrix

P k ∈ R d 2 k ×d k whose ((s -1)d k + s, s) (s = 1, • • • , d k ) entry equal to σ(u (k) s )(1 -σ(u (k) s ))(1 -2σ(u (k)

s )) and rest entries are all 0. On the other hand, since the values in v (l) belong to the range [0, 1] and y is the label, v (l) -y 2 2 can be bounded:

v (l) -y 2 2 ≤ c y < +∞, where c y is a universal constant. We further define c d = max(d 0 , d 1 , • • • , d l ).

Then we give a lemma to summarize the properties of G(u (i) ) defined in Eqn. (22), B s:t defined in Eqn. (23), D s:t and P k .

Lemma 22. For G(u (i) ) defined in Eqn. (22), B s:t defined in Eqn. (23), D s:t and P k , we have the following properties:

(1) For arbitrary matrices M and N of proper sizes, we have

G(u (i) )M 2 F ≤ 1 16 M 2 F and N G(u (i) ) 2 F ≤ 1 16 N 2 F .

(2) For arbitrary matrices M and N of proper sizes, we have

P k M 2 F ≤ 2 6 3 8 M 2 F and N P k 2 F ≤ 2 6 3 8 N 2 F .

(3) For arbitrary matrices M and N of proper sizes, we have

B s:t 2 F ≤ 1 16 t-s+1 D s:t and 1 16 t-s+1 D s:t ≤ c st ≤ c r ,

where c st = r 4 2(t-s+1) and c r = max r 2 16 , r 2 16 l-1 .

(4) For arbitrary matrices M , N and I of proper sizes, let m = vec (M ). Then we have

(N ⊗ I)m 2 F ≤ M 2 F N 2 F and (I ⊗ N )m 2 F ≤ M 2 F N 2 F .

It should be pointed out that we defer the proof of Lemma 22 to Sec. D.4.

D.2.1 PROOF OF LEMMA 16:

Proof. We use chain rule to compute the gradient of f (w, x) with respect to w (j) . We first compute several basis gradient. According to the relationship between u (j) , v (j) , W (j) and f (w, x), we have

EQUATION

Then by chain rule, we can easily compute the gradient of f (w, x) with respect to w (j) which is formulated as

∇ w (j) f (w, x) = vec v (j-1) G(u (j) )A j+1 A j+2 • • • A l (v (l) -y) T , (j = 1, • • • , l -1),

and

∇ w (l) f (w, x) = vec v (l-1) G(u (l) )(v (l) -y) T .

Besides, since the values in v (l) belong to the range [0, 1]. Combine with Lemma 22, we can bound ∇ w f (w, x) 2 as follows:

∇ w f (w, x) 2 2 = l j=1 ∇ w (j) f (w, x) 2 2 = v (l-1) G(u (l) )(v (l) -y) T 2 F + l-1 j=1 v (j-1) G(u (j) )B j+1:l (v (l) -y) T 2 F ≤ 1 16 d l-1 v (l) -y 2 2 + 1 16 v (l) -y 2 2 l-1 j=1 d j-1 B j+1:l 2 F x ≤ 1 16 c y c d + 1 16 c y c d c r (l -1),

where c y , c d , c r are defined as

v (l) -y 2 2 ≤ c y , c d = max(d 0 , d 1 , • • • , d l ) and c r = max r 2 16 , r 2 16 l-1

.

Notice, x holds since in Lemma 22, we have

B s:t 2 F ≤ r 4 2(t-s+1) ≤ max r 2 16 , r 2 16 l-1

.

Thus, we can obtain

∇ w f (w, x) 2 ≤ 1 16 c y c d (1 + c r (l -1)) α.

The proof is completed.

D.2.2 PROOF OF LEMMA 17:

For convenience, we first give the computation of some gradients. Lemma 23. Assume the activation functions in deep neural network are sigmoid functions. Then the following properties hold:

(1) We can compute the gradients ∂f (w,x) ∂u (i) and ∂f (w,x) ∂v (i) as ∂f (w, x)

∂u (i) = G(u (i) )B i+1:l (v (l) -y) and ∂f (w, x) ∂v (i) = B i+1:l (v (l) -y).

(2) We can compute the gradient ∂u (i) ∂w (j) as

∂u (i) ∂w (j) = (v (j-1) ) T ⊗ G(u (j) )B j+1:i-1 (W (i) ) T T ∈ R di×dj dj-1 , (i > j). ∂u (i) ∂w (i) = (v (i-1) ) T ⊗ I di ∈ R di×didi-1 , (i = j).

(3) We can compute the gradient ∂v (i) ∂w (j) as

∂v (i) ∂w (j) = (v (j-1) ) T ⊗ G(u (j) )B j+1:i T ∈ R di×dj dj-1 , (i ≥ j).

It should be pointed out that the proof of Lemma 23 can be founded Sec. D.4.

Proof. To prove our conclusion, we have two steps: computing the Hessian and bounding its operation norm.

Step 1. Compute the Hessian: We first consider the computation of ∂ 2 f (w,x)

EQUATION

.

Recall that we define

B s:t = A s A s+1 • • • A t ∈ R ds-1×dt , (s ≤ t) and B s:t = I, (s > t).

Then we have

∂ 2 f (w, x) ∂w T (i) ∂w (j) = v (j-1) (v (l) -y) T B T j+1:l ⊗ I dj ∂vec G(u (j) ) ∂w T (i) ( Q ij 1 ) + l k=j+1 v (j-1) (v (l) -y) T B T k+1:l ⊗ G(u (j) )B j+1:k-1 W T k ∂vec G(u (k) ) ∂w T (i) ( Q ij 2 ) + v (j-1) (v (l) -y) T B T i+1:l G(u (i) ) ⊗ G(u (j) )B j+1:i-1 ∂vec W T i ∂w T (i) ( Q ij 3 ) + v (j-1) ⊗ G(u (j) )B j+1:l ∂(v (l) -y) ∂w T (i) ( Q ij 4 ) + I dj-1 ⊗ G(u (j) )B j+1:l (v (l) -y) ∂v (j-1) ∂w T (i) ( Q ij 5 )

Case I: i > j. We first consider the case that i > j. In this is case,

Q ij 1 = 0 since ∂vec(G(u (j) )) ∂w T (i) = 0.

Computing Q ij 2 needs more efforts. By utilizing the computation of ∂u (k) ∂w (i) in Lemma 23, we have

∂vec G(u (k) ) ∂w (i) = ∂vec G(u (k) ) ∂u (k) ∂u (k) ∂w (i) = P k v (i-1) ) T ⊗ G(u (i) )B i+1:k-1 (W (k) ) T T , (k > i)

where

P k is a matrix of size d 2 k × d k whose ((s -1)d k + s, s) (s = 1, • • • , d k ) entry equal to σ(u (k) s )(1 -σ(u (k) s ))(1 -2σ(u (k) s )) and rest entries are all 0. When k = i, ∂vec G(u (k) ) ∂w (k) = ∂vec G(u (k) ) ∂u (k) ∂u (k) ∂w (k) = P k (v (k-1) ) T ⊗ I d k ∈ R d 2 k ×d k d k-1 .

Note that for k < i, we have ∂G(u (k) ) ∂w (i) = 0. For brevity, let

EQUATION

Therefore, we have

Q ij 2 = D i P i (v (i-1) ) T ⊗ I di + l k=i+1 D k P k (v (i-1) ) T ⊗ G(u (i) )B i+1:k-1 (W (k) ) T T .

Then we consider Q ij 3 .

Q ij 3 = v (j-1) (v (l) -y) T B T i+1:l G(u (i) ) ⊗ G(u (j) )B j+1:i-1 .

Also we can use the computation of ∂v (l) ∂w (i) in Lemma 23 and compute Q ij 4 as follows:

Q ij 4 =v (j-1) ⊗ G(u (j) )B j+1:l ∂(v (l) -y) ∂w T (i) = v (j-1) ⊗ G(u (j) )B j+1:l (v (i-1) ) T ⊗ G(u (i) )B i+1:l T .

Finally, since i > j, we can compute Q ij 5 = 0.

Case II: i = j. We first consider ∂G(u (k) )

∂w (k) : ∂vec G(u (k) ) ∂w T (k) = ∂vec G(u (k) ) ∂u (k) ∂u (k) ∂w T (k) = P k (v (k-1) ) T ⊗ I d k ∈ R d 2 k ×d k d k-1 ,

where P k is a matrix of size d 2 k ×d k whose (s, (s-1)d k +s) entry equal to σ(u

(k) s )(1-σ(u (k) s ))(1- 2σ(u (k)

s )) and rest entries are all 0. Q jj 1 can be computed as

Q jj 1 = v (j-1) (v (l) -y) T B T j+1:l ⊗ I dj ∂vec G(u (j) ) ∂w T (j) = v (j-1) (v (l) -y) T B T j+1:l ⊗ I dj P j (v (j-1) ) T ⊗ I dj .

As for Q jj 2 , by Eqn. ( 25) we have

Q jj 2 = l k=j+1 D k P k v (j-1) ) T ⊗ G(u (j) )B j+1:k-1 (W (k) ) T T .

Since i = j, Q jj 3 does not exist. For convenience, we just set Q jj 3 = 0. Now we consider Q jj 4 which can be computed as follows:

Q jj 4 =v (j-1) ⊗ G(u (j) )B j+1:l ∂(v (l) -y) ∂w T (j) = v (j-1) ⊗ G(u (j) )B j+1:l (v (j-1) ) T ⊗ G(u (j) )B j+1:l T .

Finally, since i = j, we can compute Q jj 5 = 0.

Case III: i < j. Since ∂ 2 f (w,x) ∂w∂w T is symmetrical, we have

Q ij k = Q ji k (k = 1, • • • , 5).

Step 2. Bound the operation norm of Hessian: We mainly use Lemma 22 to achieve this goal. From Lemma 22, we have

(1) For arbitrary matrices M and N of proper size, we have

G(u (i) )M 2 F ≤ 1 16 M 2 F and N G(u (i) ) 2 F ≤ 1 16 N 2 F .

(2) For arbitrary matrices M and N of proper size, we have

P k M 2 F ≤ 2 6 3 8 M 2 F and N P k 2 F ≤ 2 6 3 8 N 2 F .

(3) For B s:t and D s:t , we have (4) For arbitrary matrices M , N and I of proper sizes, let m = vec (M ). Then we have

B s:t 2 F ≤ 1 16 t-s+1

(N ⊗ I)m 2 F ≤ M 2 F N 2 F and (I ⊗ N )m 2 F ≤ M 2 F N 2 F .

The values of entries in v (h) are bounded by 0 ≤ σ(u

(i) h ) ≤ 1 which leads to v (h) 2 F ≤ d h ≤ c d , where c d = max i d i .

On the other hand, since the values in v (l) belong to the range [0, 1] and y is the label, v (l) -y 2 2 can be bounded: v (l) -y 2 2 ≤ c y < +∞, where c y is a universal constant.

We first define

C ij k = D k P k v (i-1) ) T ⊗ G(u (i) )B i+1:k-1 (W (k) ) T T and C ij =D i P i (v (i-1) ) T ⊗ I di = v (i-1) ⊗ I di (D i P i ) T T x = v (i-1) ⊗ (D i P i ) T T =(v (i-1) ) T ⊗ (D i P i ) ,

where D k is defined in Eqn. (25). x holds since for an arbitrary vector u ∈ R k and an arbitrary matrix M ∈ R k×k , we have (u

⊗ I k ) M = u ⊗ M .

Case I: i > j. According to the definition of C ij and C ij k , we have

Q ij 2 = C ij + l k=i+1 C ij k . So we have ∂ 2 f (w, x) ∂w T (i) ∂w (j) 2 F = Q ij 1 + Q ij 2 + Q ij 3 + Q ij 4 + Q ij 5 2 F = C ij + l k=i+1 C ij k + Q ij 3 + Q ij 4 2 F =(l -i + 3) C ij 2 F + l k=i+1 C ij k 2 F + Q ij 3 2 F + Q ij 4 2 F .

Here we bound each term separately:

C ij 2 F ≤ v (j-1) 2 F v (l) -y 2 F B i+1:l 2 F 1 16 B j+1:i-1 W T i 2 F 2 6 3 8 v (i-1) 2 F ≤ 2 6 3 8 c y d j-1 d i-1 1 16 l-i D i+1:l 1 16 i-j D j+1:i ≤ 2 6 3 8 c y d j-1 d i-1 1 16 l-j D j+1:l ≤ 2 6 3 8 c y d j-1 d i-1 c r .

Similarly, we can bound C ij k 2 F as follows:

C ij k 2 F ≤ v (j-1) 2 F v (l) -y 2 F B k+1:l 2 F 1 16 B j+1:k-1 W T k 2 F 2 6 3 8 v (i-1) 2 F 1 16 B i+1:k-1 (W (k) ) T 2 F ≤ 2 6 3 8 c y d j-1 d i-1 1 16 l-k D k+1:l 1 16 k-j-1 D j+1:k 1 16 k-i-1 D i+1:k = 2 6 3 8 c y d j-1 d i-1 1 16 l-j-1 D j+1:l 1 16 k-i-1 D i+1:k ≤ 2 14 3 8 c y d j-1 d i-1 c 2 r .

We also bound

Q ij 3 2 F as Q ij 3 2 F ≤ v (j-1) 2 F v (l) -y 2 F 1 16 B i+1:l 2 F 1 16 B j+1:i-1 2 F ≤ 1 2 8 c y d j-1 c r .

Finally, we bound Q ij 4 2 F as follows:

Q ij 4 2 F ≤ v (j-1) 2 F 1 16 B j+1:l 2 F v (i-1) 2 F 1 16 B i+1:l 2 F ≤ 1 2 8 d j-1 d i-1 c 2 r .

Note that d i ≤ c d . Thus, we can bound ∂ 2 f (w,x)

∂w (j) ∂w T (i) 2 F as ∂ 2 f (w, x) ∂w T (i) ∂w (j) 2 F ≤(l -i + 3) 2 6 3 8 c y d j-1 d i-1 c r + l k=i+1 2 14 3 8 c y d j-1 d i-1 c 2 r + 1 2 8 c y d j-1 c r + 1 2 8 d j-1 d i-1 c 2 r ≤(l + 1) 64 6561 c y c 2 d c r + 4096 6561 c y (l -2)c 2 d c 2 r + 1 256 c y c d c r + 1 256 c d c 2 r .

Case II: i = j. According to the definition of C ij and C ij k , we have

Q jj 2 = l k=j+1 C jj k . Similarly, we have ∂ 2 f (w, x) ∂w T (i) ∂w (j) 2 F = Q jj 1 + Q jj 2 + Q jj 3 + Q jj 4 + Q jj 5 2 F = Q jj 1 + l k=j+1 C jj k + Q ij 4 2 F ≤(l -j + 2)   Q jj 1 2 F + l k=j+1 C jj k 2 F + Q jj 4 2 F   .

Thus, we can bound Q jj 1 2 F first:

Q jj 1 2 F ≤ v (j-1) 2 F v (l) -y 2 F B j+1:l 2 F 2 6 3 8 v (j-1) 2 F ≤ 2 6 3 8 c y d 2 j-1 c r .

As for Q jj 2 , we have

C ij k 2 F ≤ v (j-1) 2 F v (l) -y 2 F B k+1:l 2 F 1 16 B j+1:k-1 W T k 2 F 2 6 3 8 v (j-1) 2 F 1 16 B j+1:k-1 (W (k) ) T 2 F = 2 6 3 8 c y d 2 j-1 1 16 l-k D k+1:l 1 16 k-j-1 D j+1:k 1 16 k-j-1 D j+1:k ≤ 2 14 3 8 c y d 2 j-1 c 2 r .

Then we bound Q jj 4 2 F :

Q jj 4 2 F ≤ v (j-1) 2 F 1 16 B j+1:l 2 F v (j-1) 2 F 1 16 B j+1:l 2 F ≤ 1 2 8 d 2 j-1 c 2 r .

Note that for any input, we have

EQUATION

F can be bounded by a constant c y . Thus, we can bound

∂ 2 f (w,x) ∂w T (i) ∂w (j) 2 F as ∂ 2 f (w, x) ∂w T (i) ∂w (j) 2 F ≤(l -i + 3) 2 6 3 8 c y d 2 j-1 c r + l k=i+1 2 14 3 8 c y d 2 j-1 c 2 r + 1 2 8 d 2 j-1 c 2 r ≤(l + 2) 64 6561 c y c 2 d c r + 4096 6561 c y (l -1)c 2 d c 2 r + 1 256 c 2 d c 2 r .

Case III: i < j. Since ∂ 2 f (w,x) ∂w∂w T is symmetrical, we have

Q ij k = Q ji k (k = 1, • • • , 5). Thus, it yields ∂ 2 f (w, x) ∂w T (i) ∂w (j) 2 F ≤(l + 1) 64 6561 c y c 2 d c r + 4096 6561 c y (l -2)c 2 d c 2 r + 1 256 c y c d c r + 1 256 c d c 2 r .

Final result: Thus we can bound

∇ 2 w f (w, x) op ≤ ∇ 2 w f (w, x) F ≤ (l -1)l max i,j:i =j ∂ 2 f (w, x) ∂w (j) ∂w T (i) 2 F + l j=1 ∂ 2 f (w, x) ∂w (j) ∂w T (i) 2 F ≤ (l -1)l(l+1) 64 6561 c y c 2 d c r + 4096 6561 c y (l -2)c 2 d c 2 r + 1 256 c y c d c r + 1 256 c d c 2 r +(l + 2) 64 6561 c y c 2 d c r + 4096 6561 c y (l -1)lc 2 d c 2 r + 1 256 lc 2 d c 2 r 1 2 ≤ c s1 c r c 2 d l 4

, where c s1 and c s2 are two constants.

Since ∇ 2 w f (w, x) op ≤ ∇ 2 w f (w, x) F , we know that the gradient ∇ w f (w, x) is ς-Lipschitz, where ς = c s1 c r c 2 d l 4 .

On the other hand, since for any input x, σ(x) belongs to [0, 1], the values of the entries of ∇ 3 w f (w, x) can be bounded. Thus, we can bound

∇ 3 w f (w, x) op = sup λ 2≤1 λ ⊗ 3 , ∇ 3 w f (w, x) = [∇ 3 w f (w, x)] ijk λ i λ j λ k ≤ ξ < +∞.

We complete the proof.

D.2.3 PROOF OF LEMMA 18:

For convenience, we first give the computation of some gradients.

Lemma 24. Assume the activation functions in deep neural network are sigmoid functions. Then we can compute the gradients ∂u (j) ∂u (1) and ∂v (j) ∂u (1) as

∂u (j) ∂u (1) = G(u (1) )A 2 • • • A j-1 (W j ) T T ∈ R dj ×d1 , (j > 1)

.

∂v (j) ∂u (1) = G(u (1) )A 2 • • • A j T ∈ R dj ×d1 , (j > 1).

It should be pointed out that the proof of Lemma 24 can be founded Sec. D.4.

Proof. To prove our conclusion, we have two steps: computing ∇ x ∇ w f (w, x) and bounding its operation norm.

Step 1. Compute ∇ x ∇ w f (w, x):

We first consider the computation of ∂ 2 f (w,x) ∂x T ∂w (j) :

∂ 2 f (w, x) ∂x T ∂w (j) = ∂ vec G(u (j) )A j+1 A j+2 • • • A l (v (l) -y) (v (j-1) ) T ∂x T .

Recall that we define

A i = (W (i) ) T G(u (i) ) ∈ R di-1×di . B s:t = A s A s+1 • • • A t ∈ R ds-1×dt

, (s ≤ t) and B s:t = I, (s > t).

Then we have

∂ 2 f (w, x) ∂x T ∂w (j) = v (j-1) (v (l) -y) T B T j+1:l ⊗ I dj ∂vec G(u (j) ) ∂x T ( Q j 1 ) + l k=j+1 v (j-1) (v (l) -y) T B T k+1:l ⊗ G(u (j) )B j+1:k-1 W T k ∂vec G(u (k) ) ∂x T ( Q j 2 ) + v (j-1) ⊗ G(u (j) )B j+1:l ∂(v (l) -y) ∂x T ( Q j 3 ) + I dj-1 ⊗ G(u (j) )B j+1:l (v (l) -y) ∂v (j-1) ∂x T ( Q j 4 )

By using Lemma 24, we can compute

Q ij 1 as ∂vec G(u (k) ) ∂x T = ∂vec G(u (k) ) ∂u (k) ∂u (k) ∂x T = P k G(u (1) )B 2:k-1 (W k ) T T .

Thus, we have

Q j 1 = v (j-1) (v (l) -y) T B T j+1:l ⊗ I dj ∂vec G(u (j) ) ∂x T = v (j-1) (v (l) -y) T B T j+1:l ⊗ I dj P k G(u (1) )B 2:k-1 (W k ) T T .

As for Q j 2 , we also can utilize Lemma 24 to compute it:

Q j 2 = l k=j+1 v (j-1) (v (l) -y) T B T k+1:l ⊗ G(u (j) )B j+1:k-1 W T k ∂vec G(u (k) ) ∂x T = l k=i+1 v (j-1) (v (l) -y) T B T k+1:l ⊗ G(u (j) )B j+1:k-1 W T k P k G(u (1) )B 2:k-1 (W k ) T T .

Then we consider Q ij 3 .

Q j 3 = v (j-1) ⊗ G(u (j) )B j+1:l ∂(v (l) -y) ∂x T = v (j-1) ⊗ G(u (j) )B j+1:l G(u (1) )B 2:l T .

Q j 4 can be computed as follows:

Q j 4 = I dj-1 ⊗ G(u (j) )B j+1:l (v (l) -y) ∂v (j-1) ∂x T = I dj-1 ⊗ G(u (j) )B j+1:l (v (l) -y) G(u (1) )B 2:j T .

Step 2. Bound the operation norm of Hessian: We mainly use Lemma 22 to achieve this goal.

From Lemma 22, we have

(1) For arbitrary matrices M and N of proper size, we have

G(u (i) )M 2 F ≤ 1 16 M 2 F and N G(u (i) ) 2 F ≤ 1 16 N 2 F .

(2) For arbitrary matrices M and N of proper size, we have

P k M 2 F ≤ 2 6 3 8 M 2 F and N P k 2 F ≤ 2 6 3 8 N 2 F .

(3) For B s:t and D s:t , we have (4) For arbitrary matrices M , N and I of proper sizes, let m = vec (M ). Then we have

(N ⊗ I)m 2 F ≤ M 2 F N 2 F and (I ⊗ N )m 2 F ≤ M 2 F N 2 F .

The values of entries in v (h) are bounded by 0 ≤ σ(u

(i) h ) ≤ 1 which leads to v (h) 2 F ≤ d h ≤ c d , where c d = max i d i .

On the other hand, since the values in v (l) belong to the range [0, 1] and y is the label, v (l) -y 2 2 can be bounded: v (l) -y 2 2 ≤ c y < +∞, where c y is a universal constant.

We first define

C j k = v (j-1) (v (l) -y) T B T k+1:l ⊗ G(u (j) )B j+1:k-1 W T k P k G(u (1) )B 2:k-1 (W k ) T T .

Then we have

Q j 2 = l k=j+1 C j k . So we have ∂ 2 f (w, x) ∂x T ∂w (j) 2 F = Q j 1 + Q j 2 + Q j 3 + Q j 4 2 F = Q j 1 + l k=j+1 C j k + Q j 3 + Q j 4 2 F =(l -j + 3)   Q j 1 2 F + l k=j+1 C j k 2 F + Q j 3 2 F + Q j 4 2 F   .

Then we bound each term separately:

Q j 1 2 F ≤ v (j-1) 2 F v (l) -y 2 F B j+1:l 2 F 2 6 3 8 1 16 B 2:k-1 (W k ) T 2 F ≤ 2 6 3 8 c y d j-1 c 2 r .

Similarly, we bound

C j k 2 F : C j k 2 F = v (j-1) 2 F v (l) -y 2 F B k+1:l 2 F 1 16 B j+1:k-1 W T k 2 F 2 6 3 8 1 16 B 2:k-1 (W (k) ) T 2 F = 2 6 3 8 c y d j-1 1 16 l-k D k+1:l 1 16 k-j-1 D j+1:k 1 16 k-1 D 2:k ≤ 2 6 3 8 c y d j-1 c 2 r .

We also bound

Q ij 3 2 F as Q ij 3 2 F ≤ v (j-1) 2 2 1 16 B j+1:l 2 F 1 16 B 2:l 2 F ≤ 1 2 8 d j-1 c 2 r .

Finally, we bound Q j 4 2 F as follows:

Q j 4 2 F = 1 16 B j+1:l 2 F v (l) -y 2 F 1 16 B 2:j 2 F ≤ 1 2 8 c y c r . Since c d = max i d i , we can bound ∂ 2 f (w,x) ∂w (j) ∂x T 2 F as ∂ 2 f (w, x) ∂x T ∂w (j) 2 F ≤(l -j + 3)   2 6 3 8 c y d j-1 c 2 r + l k=j+1 2 6 3 8 c y d j-1 c 2 r + 1 2 8 c y d j-1 c r + 1 2 8 c y c r   ≤(l + 2)   2 6 3 8 c y d j-1 c 2 r + l k=j+1 2 6 3 8 c y d j-1 c 2 r + 1 2 8 c y d j-1 c r + 1 2 8 c y c r   .

Final result: Thus we can bound

∇ w ∇ x f (w, x) op ≤ ∇ w ∇ x f (w, x) F ≤ l j=1 ∂ 2 f (w, x) ∂w (j) ∂x T 2 F ≤ l j=1 (l + 2)   2 6 3 8 c y d j-1 c 2 r + l k=j+1 2 6 3 8 c y d j-1 c 2 r + 1 2 8 c y d j-1 c r + 1 2 8 c y c r   ≤ 2 6 3 8 l(l + 2)c y c r c d (lc r + 1)

, where c d = max j d j . The proof is completed.

D.2.4 PROOF OF LEMMAS 19 AND 20:

Lemma 25. (Alessandro, 2016; Rigollet, 2015) Let (x 1 , • • • , x k ) be a vector of i.i.d. Gaussian variables from N (0, τ 2 ) and let f : R d0 → R be L-Lipschitz. Then the variable f (x) -Ef (x) is sub-Gaussian. That is, we have

P (f (x) -Ef (x) > t) ≤ exp - t 2 2L 2 τ 2 , (∀t ≥ 0), or E (λ(f (x) -Ef (x))) ≤ exp 4λ 2 L 2 τ 2 , (∀λ ≥ 0).

Remarkably, this is a dimension free inequality.

Proof of Lemma 19. We first define a function g(x) = z T ∇ w f (w, x) where z ∈ R d is a constant vector. Then we have ∇

x g(x) = ∇ x z T ∇ w f (w, x) = ∇ x ∇ w f (w, x)z. Then by Lemma 18, we can obtain ∇ x g(x) 2 = ∇ x ∇ w f (w, x)z 2 ≤ β z 2 , where β = 2 6

3 8 l(l + 2)c y c r c d (lc r + 1) in which c y , c d and c r are defined in Lemma 18. This means g(x) is β z 2 -Lipschitz. Thus, by Lemma 25, we have

E (t z, ∇ w f (w, x) -E∇ w f (w, x) ) = E (t (g(x) -Eg(x))) ≤ exp 4t 2 β 2 z 2 2 τ 2 .

Let λ = tz. This further gives

E ( λ, ∇ w f (w, x) -E∇ w f (w, x) ) ≤ exp 4β 2 τ 2 λ 2 2 , which means λ, ∇ w f (w, x) -E∇ w f (w, x) is 8β 2 τ 2 -sub-Gaussian.

Proof of Lemma 20. We first define a function h

(x) = z T ∇ 2 w f (w, x)z where z ∈ S d-1 , i.e. z 2 = 1. Then h(w) is a γ-Lipschitz function, where γ = ∇ x ∇ 2

w f (w, x) op . Note that since the sigmoid function is infinitely differentiable function, ∇ x ∇ 2 w f (w, x) exists. Also since for any input x, σ(x) belongs to [0, 1]. Thus, the values of the entries in ∇ x ∇ 2 w f (w, x) can be bounded. So according to the definition of the operation norm of a 3-way tensor, the operation norm of ∇ x ∇ 2 w f (w, x) can be bounded by a constant. Without loss of generality, let ∇ x ∇ 2 w f (w, x) op ≤ γ < +∞. Thus, by Lemma 25, we have

E t z, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) z = E (t (h(x) -Eh(x))) ≤ exp 8t 2 γ 2 τ 2 2 .

This means that the hessian of the loss evaluated on a unit vector is 8γ 2 τ 2 -sub-Gaussian.

D.2.5 PROOF OF LEMMA 21:

Proof. Recall that the weight of each layer has magnitude bound separately, i.e. w (j) 2 ≤ r.

Assume that w (j) has s j non-zero entries. Then we have l j=1 s j = s. So here we separately assume

EQUATION

r) which corresponds to the weight w (j) of the j-th layer. Let n j be the /l-covering number. By -covering theory in (Vershynin, 2012) , we can have

n j ≤ d j d j-1 s j 3r d j d j-1 /(ld) sj ≤ exp s j log 3rd j d j-1 d j d j-1 /d = exp s j log 3rd .

Let w ∈ Ω be an arbitrary vector. Since w = [w (1) , • • • , w (l) ] where w (j) is the weight of the j-th layer, we can always find a vector w j kj in w j such that w (j) -w j kj 2 ≤ d j d j-1 /d. For brevity, let j w ∈ [n j ] denote the index of w j kj in -net w j . Then let

w kw = [w j k1 ; • • • ; w j kj ; • • • ; w j k l ]

. This means that we can always find a vector w kw such that ww kw 2 ≤ . Accordingly, we can decompose ∇ 2 Ĵn (w) -∇ 2 J (w) op as follows:

∇ 2 Ĵn (w) -∇ 2 J (w) op = 1 n n i=1 ∇ 2 f (w, x (i) ) -E(∇ 2 f (w, x)) op = 1 n n i=1 ∇ 2 f (w, x (i) ) -∇f (w kw , x (i) ) + 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) + E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) op ≤ 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) op + 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) op + E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) op .

Here we also define four events E 0 , E 1 , E 2 and E 3 as

E 0 = sup w∈Ω ∇ 2 Ĵn (w) -∇ 2 J (w) op ≥ t , E 1 =    sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) op ≥ t 3    , E 2 =    sup jw∈[n j ],j=[l] 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) op ≥ t 3    , E 3 = sup w∈Ω E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) op ≥ t 3 .

Accordingly, we have P (E 0 ) ≤ P (E 1 ) + P (E 2 ) + P (E 3 ) . So we can respectively bound P (E 1 ), P (E 2 ) and P (E 3 ) to bound P (E 0 ).

Step 1. Bound P (E 1 ): We first bound P (E 1 ) as follows:

P (E 1 ) =P sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) 2 ≥ t 3 x ≤ 3 t E sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) 2 ≤ 3 t E sup w∈Ω 1 n n i=1 ∇ 2 f (w, x (i) ) -∇ 2 f (w kw , x (i) ) 2 w -w kw 2 sup w∈Ω w -w kw 2 ≤ 3 t E   sup w∈Ω 1 n n i=1 ∇ 3 f (w, x (i) ) op   y ≤ 3ξ t ,

where x holds since by Markov inequality and y holds because of Lemma 17.

Therefore, we can set t ≥ 6ξ ε .

Then we can bound P(E 1 ):

P(E 1 ) ≤ ε 2 .

Step 2. Bound P (E 2 ): By Lemma 2, we know that for any matrix X ∈ R d×d , its operator norm can be computed as

X op ≤ 1 1 -2 sup λ∈λ | λ, Xλ | .

where λ = {λ 1 , . . . , λ kw } be an -covering net of B d (1).

Let λ 1/4 be the 1 4 -covering net of B d (1) but it has only s nonzero entries. So the size of its -net is

d s 3 1/4 s ≤ exp (s log (12d)) .

Recall that we use j w to denote the index of w j kj in -net w j and we have j w ∈ [n j ], (n j ≤ exp s j log 3rd ). Then we can bound P (E 2 ) as follows:

P(E 2 ) =P sup jw∈[n j ],j=[l] 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E(∇ 2 f (w kw , x)) 2 ≥ t 3 =P sup jw∈[n j ],j=[l],λ∈λ 1/4 2 λ, 1 n n i=1 ∇ 2 f (w kw , x (i) ) -E ∇ 2 f (w kw , x) λ ≥ t 3 ≤ exp (s log (12d)) exp   l j=1 s j log 3rd   sup jw∈[n j ],j=[l],λ∈λ 1/4 P 1 n n i=1 λ, ∇ 2 f (w kw , x (i) )-E ∇ 2 f (w kw , x) λ ≥ t 6 . Since by Lemma 20, λ, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) λ where λ ∈ B d (1) is 8γ 2 τ 2 -sub- Gaussian, i.e. E t λ, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) λ ≤ exp 8t 2 γ 2 τ 2 2 . Thus, 1 n n i=1 λ, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) λ is 8γ 2 τ 2 /n-sub-Gaussian random variable. So we can obtain P 1 n n i=1 y, ∇ 2 w f (w, x) -E∇ 2 w f (w, x) y ≥ t 6 ≤ 2 exp - nt 2 72γ 2 τ 2 . Note d = j d j d j-1 . Then the probability of E 2 is upper bounded as P (E 2 ) ≤ 2 exp - nt 2 72γ 2 τ 2 + s log 36d 2 r .

Thus, if we set t ≥ γτ 72 (s log(36d 2 r/ ) + log(4/ε)) n ,

then we have

P (E 2 ) ≤ ε 2 .

Step 3. Bound P (E 3 ): We first bound P (E 3 ) as follows:

P (E 3 ) =P sup w∈Ω E(∇ 2 f (w kw , x)) -E(∇ 2 f (w, x)) 2 ≥ t 3 ≤P E sup w∈Ω (∇ 2 f (w kw , x) -∇ 2 f (w, x) 2 ≥ t 3 =P E sup w∈Ω ∇ 2 f (w, x) -∇ 2 f (w kw , x) 2 w -w kw 2 sup w∈Ω w -w kw 2 ≥ t 3 ≤P E sup w∈Ω ∇ 3 f (w, x) op ≥ t 3 ≤P ξ ≥ t 3 .

We set enough small such that ξ < t/3 always holds. Then it yields P (E 3 ) = 0.

Step 4. Final result: To ensure P(E 0 ) ≤ ε, we just set = 36rl Proof. Recall that the weight of each layer has magnitude bound separately, i.e. w (j) 2 ≤ r.

Assume that w (j) has s j non-zero entries. Then we have l j=1 s j = s. So here we separately assume

EQUATION

r) which corresponds to the weight w (j) of the j-th layer. Let n j be the /l-covering number. By -covering theory in (Vershynin, 2012) , we can have

n j ≤ d j d j-1 s j 3r d j d j-1 /(ld) sj ≤ exp s j log 3rd j d j-1 d j d j-1 /d = exp s j log 3rd .

Let w ∈ Ω be an arbitrary vector. Since w = [w (1) , • • • , w (l) ] where w (j) is the weight of the j-th layer, we can always find a vector w j kj in w j such that w (j) -w j kj 2 ≤ d j d j-1 /d. For brevity, let j w ∈ [n j ] denote the index of w j kj in -net w j . Then let

w kw = [w j k1 ; • • • ; w j kj ; • • • ; w j k l ]

. This means that we can always find a vector w kw such that ww kw 2 ≤ . Accordingly, we can decompose ∇ Ĵn (w) -∇J (w) 2 as follows:

∇ Ĵn (w) -∇J (w) 2 = 1 n n i=1 ∇f (w, x (i) ) -E(∇f (w, x)) 2 = 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) + 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) + E(∇f (w kw , x)) -E(∇f (w, x)) 2 ≤ 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 + 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) 2 + E(∇f (w kw , x)) -E(∇f (w, x)) 2 .

Here we also define four events E 0 , E 1 , E 2 and E 3 as

E 0 = sup w∈Ω ∇ Ĵn (w) -∇J (w) 2 ≥ t , E 1 = sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 ≥ t 3 , E 2 = sup jw∈[n j ],j=[l] 1 n n i=1 ∇f (w kw , x (i) ) -E(∇f (w kw , x)) 2 ≥ t 3 , E 3 = sup w∈Ω E(∇f (w kw , x)) -E(∇f (w, x)) 2 ≥ t 3 .

Accordingly, we have P (E 0 ) ≤ P (E 1 ) + P (E 2 ) + P (E 3 ) . So we can respectively bound P (E 1 ), P (E 2 ) and P (E 3 ) to bound P (E 0 ).

Step 1. Bound P (E 1 ): We first bound P (E 1 ) as follows:

P (E 1 ) =P sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 ≥ t 3 x ≤ 3 t E sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 ≤ 3 t E sup w∈Ω 1 n n i=1 ∇f (w, x (i) ) -∇f (w kw , x (i) ) 2 w -w kw 2 sup w∈Ω w -w kw 2 ≤ 3 t E sup w∈Ω ∇ 2 Ĵn (w, x) 2 ,

where x holds because of Markov inequality. Then, we bound E sup w∈Ω ∇ 2 Ĵn (w, x) 2 as follows:

E sup w∈Ω ∇ 2 Ĵn (w, x) 2 ≤ E sup w∈Ω 1 n n i=1 ∇ 2 f (w,x) 2 = E sup w∈Ω ∇ 2 f (w,x) 2 x ≤ ς,

where x holds since by Lemma 17, we have

∇ 2 w f (w, x) op ≤ ∇ 2 w f (w, x) F ≤ς, where ς = c s1 c r c 2 d l 4 in which c d = max i d i and c r = max r 2 16 , r 2 16 l-1

. Therefore, we have

P (E 1 ) ≤ 3ς t .

We further let t ≥ 6ς ε .

Then we can bound P(E 1 ):

P(E 1 ) ≤ ε 2 .

Step 2. Bound P (E 2 ): By Lemma 1, we know that for any vector x ∈ R d , its 2 -norm can be computed as

x 2 ≤ 1 1 - sup λ∈λ λ, x .

where λ = {λ 1 , . . . , λ kw } be an -covering net of B d (1).

Let λ 1/2 be the 1 2 -covering net of B d (1) but it has only s nonzero entries. So the size of its -net is d s 3 1/2 s ≤ exp (s log (6d)) .

Recall that we use j w to denote the index of w j kj in -net w j and we have j w ∈ [n j ], (n j ≤ exp s j log 3rd . Then we can bound P (E 2 ) as follows: Notice, j d j d j-1 = d. In this case, the probability of E 2 is upper bounded as

P (E 2 ) =P sup jw∈[n j ],j=[l]

P (E 2 ) ≤ exp - nt 2 72β 2 τ 2 + d log 18r

.

Thus, if we set t ≥ βτ 72 (s log(18d 2 r/ ) + log(4/ε)) n , then we have P (E 2 ) ≤ ε 2 .

Step 3. Bound P (E 3 ): We first bound P (E 3 ) as follows:

P (E 3 ) =P sup where x holds since by Lemma 17. We set enough small such that ς < t/3 always holds. Then it yields P (E 3 ) = 0.

Step 4. Proof. Suppose that {w (1) , w (2) , • • • , w (m) } are the non-degenerate critical points of J (w). So for any w (k) , it obeys

inf i λ k i ∇ 2 J (w (k) ) ≥ ζ,

where λ k i ∇ 2 J (w (k) ) denotes the i-th eigenvalue of the Hessian ∇ 2 J (w (k) ) and ζ is a constant. We further define a set D = {w ∈ R d | ∇J (w) 2 ≤ and inf i |λ i ∇ 2 J (w to the weight w (j) of the j-th layer. Let n j be the /l-covering number. By -covering theory in (Vershynin, 2012) , we can have n j ≤ d j d j-1 s j 3r d j d j-1 /(ld) sj ≤ exp s j log 3rd j d j-1 d j d j-1 /d = exp s j log 3rd .

Let w ∈ Ω be an arbitrary vector. Since w = [w (1) , • • • , w (l) ] where w (j) is the weight of the j-th layer, we can always find a vector w j kj in w j such that w (j) -w j kj 2 ≤ d j d j-1 /d. For brevity, let j w ∈ [n j ] denote the index of w j kj in -net w j . Then let w kw = [w j k1 ; • • • ; w j kj ; • • • ; w j k l ]. This means that we can always find a vector w kw such that ww kw 2 ≤ . Accordingly, we can decompose Ĵn (w) -J (w) as f (w kw , x (i) )-Ef (w kw , x) + Ef (w kw , x)-Ef (w, x) .

Then, we define four events E 0 , E 1 , E 2 and E 3 as Accordingly, we have P (E 0 ) ≤ P (E 1 ) + P (E 2 ) + P (E 3 ) . So we can respectively bound P (E 1 ), P (E 2 ) and P (E 3 ) to bound P (E 0 ).

Step 1. Bound P (E 1 ): We first bound P (E 1 ) as follows: We further let t ≥ 6α ε .

P (E 1 ) =P sup w∈Ω 1 n n i=1 f (w, x (i) ) -f (w kw , x (i) ) ≥ t 3 x ≤ 3 t E

Then we can bound P(E 1 ):

P(E 1 ) ≤ ε 2 .

Step 2. Bound P (E 2 ): Recall that we use j w to denote the index of w j kj in -net w j and we have j w ∈ [n j ], (n j ≤ exp s j log 3rd . We can bound P (E 2 ) as follows:

P (E 2 ) =P sup f (w j , x (i) ) -E(f (w j , x)) ≥ t 3 .

Since when the activation functions are sigmoid functions, the loss f (w, x) is α-Lipschitz. Besides, we assume x to be a vector of i.i.d. Gaussian variables from N (0, τ 2 ). Then by Lemma 25, we know that the variable f (x) -Ef (x) is 8α 2 τ 2 -sub-Gaussian. Thus, we have f (w j , x (i) ) -E(f (w j , x)) ≥ t 3 ≤ 2 exp -nt 2 18α 2 τ 2 . Notice l j=1 s j = s. In this case, the probability of E 2 is upper bounded as

P (E 2 ) ≤ 2 exp - nt 2 18α 2 τ 2 + s log 3dr

.

Thus, if we set t ≥ ατ 18 (s log(3dr/ ) + log(4/ε)) n ,

then we have P (E 2 ) ≤ ε 2 .

Step 3. Bound P (E 3 ): We first bound P (E 3 ) as follows:

P (E 3 ) =P sup where x holds since by Lemma 16, for arbitrary x and w ∈ Ω, we have ∇ w f (w, x) 2 ≤ α. We set enough small such that α < t/3 always holds. Then it yields P (E 3 ) = 0.

Step 4. Final result: Notice, we have 6α ε ≥ 3α . To ensure P(E 0 ) ≤ ε, we just set = 3r/n and t ≥ max 6α ε , ατ 18 (s log(3dr/ ) + log(4/ε)) n = max 18αr nε , ατ 18 (s log(nd)+log(4/ε)) n .

Therefore, if n ≥ 18l 2 r 2 /(s log(d)τ 2 ε 2 log(1/ε)), then Proof. Since G(u (i) ) is a diagonal matrix and its diagonal values are upper bounded by σ(u

(i) h )(1 - σ(u (i) h )) ≤ 1/4 where u (i)

h denotes the h-th entry of u (i) , we can conclude

G(u (i) )M 2 F ≤ 1 16 M 2 F and N G(u (i) ) 2 F ≤ 1 16 N 2 F .

Note that P k is a matrix of size d 2 k × d k whose ((s -1)

d k + s, s) (s = 1, • • • , d k ) entry equal to σ(u (k) s )(1 -σ(u (k) s ))(1 -2σ(u (k)

s )) and rest entries are all 0. This gives

σ(u (k) s )(1 -σ(u (k) s ))(1 -2σ(u (k) s )) = 1 3 (3σ(u (k) s ))(1 -σ(u (k) s ))(1 -2σ(u (k) s )) ≤ 1 3 3σ(u (k) s ) + 1 -σ(u (k) s ) + 1 -2σ(u (k) s ) 3 3 ≤ 2 3 3 4 .

This means the maximal value in P k is at most 2 3 3 4 . Consider the structure in P k , we can obtain

P k M 2 F ≤ 2 6 3 8 M 2 F and N P k 2 F ≤ 2 6 3 8 N 2 F .

As for B s:t , we have

B s:t 2 F ≤ A s 2 F A s+1 2 F • • • A t 2 F = (W s ) T G(u (s) ) 2 F (W (s+1) ) T G(u (s+1) ) 2 F • • • (W (t) ) T G(u (t) ) 2 F ≤ 1 16 t-s+1 W (s) 2 F W (s+1) 2 F • • • W (t) 2 F = 1 16 t-s+1 D s:t .

Since the 2 -norm of each w (j) is bounded, i.e. w (j) 2 ≤ r, we can obtain 1 16 t-s+1 D s:t ≤ 1 16 t-s+1 r 2(t-s+1) = r 4 2(t-s+1) c st . Now we prove the final result. According to the property of Kronecker product that for any matrices A, B and X of proper sizes, vec (AXB) = (B T ⊗ A)vec (X), we have vec M N T = (N ⊗ I)vec (M ) = (N ⊗ I)m. This further yields

(N ⊗ I)m 2 F = vec M N T 2 F = M N T 2 F ≤ M 2 F N 2 F .

By similar way, we can obtain

(I ⊗ N )m 2 F ≤ M 2 F N 2 F .

The proof is completed.

D.4.2 PROOF OF LEMMA 23:

Proof. By utilizing the chain rule in Eqn. (24) in Sec. D.2.1, we can easily compute ∂f (w,x) ∂u (i) and ∂f (w,x) ∂v (i) as follows: ∂f (w, x) ∂u (i) = G(u (i) )A i+1 • • • A l (v (l) -y) = G(u (i) )B i+1:l (v (l) -y) and ∂f (w, x) l) y) = B i+1:l (v (l) -y). Therefore, we can further obtain ∂f (w, x)

∂v (i) = A i+1 • • • A l (v (

∂w (j) =vec G(u (j) )A j+1 A j+2 • • • A l (v (l) -y) (v (j-1) ) T =vec G(u (j) )A j+1 A j+2 • • • A i-1 (W (i) ) T G(u (i) )A i+1 • • • A l (v (l) -y) (v (j-1) ) T = v (j-1) ⊗ G(u (j) )A j+1 A j+2 • • • A i-1 (W (i) ) T vec G(u (i) )A i+1 • • • A l (v (l) -y) = v (j-1) ⊗ G(u (j) )A j+1 A j+2 • • • A i-1 (W (i) ) T ∂f (w, x)

∂u (i) .

Note that we have ∂f (w,x) ∂w (j) = ∂u (i) ∂w (j) ∂f (w,x) ∂u (i) . This gives ∂u (i) ∂w (j) = (v (j-1) ) T ⊗ G(u (j) )B j+1:i-1 (W (i) ) T T ∈ R di×dj dj-1 (i > j).

When i = j, we have ∂u (i) ∂w (i) = (v (i-1) ) T ⊗ I di ∈ R di×didi-1 .

Similarly, we can obtain ∂v (i) ∂w (j) = (v (j-1) ) T ⊗ G(u (j) )A j+1 A j+2 • • • A i T =(v (j-1) ) T ⊗ G(u (j) )B j+1:i T ∈R di×dj dj-1 (i ≥ j).

The proof is completed.

:

≤ r 2 .

Figures and Tables:

D s:t and1 16 t-s+1 D s:t ≤ c r ,where c r = max r 2 16 , r 2 16 l .

s+1 D s:t and1 16 t-s+1 D s:t ≤ c r ,where c r = max r 2 4 ,

kw , x (i) ) -E(∇f (w kw , x)) kw , x (i) ) -E (∇f (w kw , x)) kw , x (i) ) -E (∇f (w kw , x)) ≥ t 6 .Since by Lemma 19, y, ∇f (w, x) is 8β 2 τ 2 -sub-Gaussian, i.e.E ( λ, ∇ w f (w, x) -E∇ w f (w, x) ) (l + 2)c y c r c d (lc r + 1) in which c y , c d and c r are defined in Lemma 16. Thus, ∇f (w, x) is 8β 2 τ 2 /n-sub-Gaussian random variable. Thus, we can obtain (w kw , x (i) ) -E (∇f (w kw , x))

w∈Ω E(∇f (w kw , x)) -E(∇f (w, x)) (∇f (w kw , x) -∇f (w, x) 2 )w -

Final result: To ensure P(E 0 ) ≤ ε, we just set = 18r/n and t ≥ max 6ς ε , βτ 72 (s log(18d 2 r/ ) + log(4/ε)) n = max 108ςr nε , βτ 72 (s log(nl) + log(4/ε)) n .Note that ς = O( √ lc d β). Therefore, there exists a universal constant c y such that if n ≥ c y c d l 3 r 2 /(s log(d)τ 2 ε 2 log(1/ε)), (l + 2) (lc r + 1) c r c d s log(dn/l)+log(4/ε) n holds with probability at least 1 -ε, where c y , c d and c r are defined in Lemma 16. D.3.2 PROOF OF THEOREM 5

(k) ) | ≥ ζ}. According to Lemma 4, D = ∪ ∞ k=1 D k where each D k is a disjoint component with w (k) ∈ D k for k ≤ m and D k does not contain any critical point of J (w) for k ≥ m + 1. On the other hand, by the continuity of ∇J (w), it yields ∇J (w) 2 = for w ∈ ∂D k . Notice, we set the value of blow which is actually a function related n. Then by utilizing Theorem 4, we let sample number n sufficient large such that sup w∈Ω ∇ Ĵn (w) -∇J (w) at least 1 -ε, where β = τ 512 729 c y l(l + 2) (lc r + 1) c r c d s log(dn/l)+log(4/ε) n . This further gives that for arbitrary w ∈ D k , we have inf w∈D k t∇ Ĵn (w) + (1 -t)∇J (w) 2 = inf w∈D k t ∇ Ĵn (w) -∇J (w) + ∇J (w) 2 ≥ inf w∈D k ∇J (w) 2 -sup w∈D k t ∇ Ĵn (w) -∇J (w)

, x (i) )-f (w kw , x (i) ) + 1 n n i=1 f (w kw , x (i) )-Ef (w kw , x)+Ef (w kw , x)-Ef (w, x) , x (i) )-f (w kw , x (i) ) + 1 n n i=1

x (i) ) -f (w kw , x (i) kw , x (i) )-E(f (w kw , x)) ≥ t 3 , E 3 = sup w∈Ω E(f (w kw , x))-E(f (w, x)) ≥ t 3 .

c d (1 + c r (l -1)) in which c y , c d and c r are defined in Lemma 16. Therefore, we haveP (E 1 ) ≤ 3α t .

kw , x (i) ) -E(f (w kw , x))

(|f (x) -Ef (x)| > t) ≤ 2 exp -t 2 2α 2 τ 2 , (∀t ≥ 0),where α = 1 16 c y c d (1 + c r (l -1)) in which c y , c d and c r are defined in Lemma 16. Thus,1 n n i=1 f (w j , x (i) ) -E(f (w j , x)) is 8α 2 τ 2 /n-sub-Gaussian random variable. Thus, we can obtain

w∈Ω |E(f (w kw , x)) -E(f (w, x))| ≥ t 3 =P sup w∈Ω |E (f (w kw , x) -f (w, x))| ww kw 2 sup w∈Ω ww kw 2 ≥ t 3 ≤P E sup w∈Ω ∇J w (w, x) 2 ≥ t 3 x ≤P α ≥ t 3 ,

d (1 + c r (l -1)) s log(nd/l) + log(4/ε) nholds with probability at least 1 -ε, where c y , c d , and c r are defined asv (l) -y 2 2 ≤ c y < +∞, c d = max(d 0 , d 1 , • • • , d l ) and c r = maxThe proof is completed.D.3.4 PROOF OF COROLLARY 2Proof. By Lemma 5, we know s = g . Thus, the remaining work is to bound s . Actually, we can haveE S∼D,A,(x (j) ,y (j) w j * ;x (j) ,y (j) )-f j (w n ;x (j) ,y(j) ) ≤E S∼D sup w∈Ω Ĵn (w) -J (w) ≤ sup w∈Ω Ĵn (w) -J (w) ≤ n . Thus, we have g = s ≤ n . The proof is completed. D.4 PROOF OF OTHER LEMMAS D.4.1 PROOF OF LEMMA 22

sup Now we only need to bound E sup w∈Ω ∇ Ĵn (w, x)

Table content: <table><tr><td/><td/><td/><td/><td/><td/><td/><td>. Then by Lemma 16, we can bound it as</td></tr><tr><td>follows:</td><td/><td/><td/><td/><td/><td/><td>2</td></tr><tr><td colspan="4">E sup</td><td colspan="3">∇ Ĵn (w, x)</td></tr><tr><td colspan="4">w∈Ω</td><td/><td/><td/></tr><tr><td/><td/><td/><td/><td>w∈Ω</td><td>1 n</td><td>n i=1</td><td>f (w, x (i) ) -f (w kw , x (i) )</td></tr><tr><td>≤</td><td>3 t</td><td colspan="3">E sup w∈Ω</td><td>1 n</td><td colspan="2">n i=1 f (w, x (i) ) -f (w kw , x (i) ) w -w kw 2</td><td>sup w∈Ω</td><td>w -w kw 2</td></tr><tr><td>≤</td><td colspan="2">3 t</td><td colspan="2">E sup</td><td/><td/></tr></table>

Peer review (iclr2018):
 This paper provides the analysis of empirical risk landscape for GENERAL deep neural networks (DNNs). Assumptions are comparable to existing results for OVERSIMPLIFED shallow neural networks. The main results analyzed: 1) Correspondence of non-degenerate stationary points between empirical risk and the population counterparts. 2) Uniform convergence of the empirical risk to population risk. 3) Generalization bound based on stability. The theory is first developed for linear DNNs and then generalized to nonlinear DNNs with sigmoid activations. Here are two detailed comments: 1) For deep linear networks with squared loss, Kawaguchi 2016 has shown that the global optima are the only non-degerenate stationary points. Thus, the obtained non-degerenate stationary deep linear network should be equivalent to the linear regression model Y=XW. Should the risk bound only depends on the dimensions of the matrix W? 2) The comparison with Bartlett & Maass’s (BM) work is a bit unfair, because their result holds for polynomial activations while this paper handles linear activations. Thus, the authors need to refine BM*s result for comparison.