BENEFIT OF DEEP LEARNING WITH NON-CONVEX NOISY GRADIENT DESCENT: PROVABLE EXCESS RISK BOUND AND SUPERIORITY TO KERNEL METHODS
Year: 2021

Abstract:
Establishing a theoretical analysis that explains why deep learning can outperform shallow learning such as kernel methods is one of the biggest issues in the deep learning literature. Towards answering this question, we evaluate excess risk of a deep learning estimator trained by a noisy gradient descent with ridge regularization on a mildly overparameterized neural network, and discuss its superiority to a class of linear estimators that includes neural tangent kernel approach, random feature model, other kernel methods, k-NN estimator and so on. We consider a teacher-student regression model, and eventually show that any linear estimator can be outperformed by deep learning in a sense of the minimax optimal rate especially for a high dimension setting. The obtained excess bounds are so-called fast learning rate which is faster than O(1/ √ n) that is obtained by usual Rademacher complexity analysis. This discrepancy is induced by the non-convex geometry of the model and the noisy gradient descent used for neural network training provably reaches a near global optimal solution even though the loss landscape is highly non-convex. Although the noisy gradient descent does not employ any explicit or implicit sparsity inducing regularization, it shows a preferable generalization performance that dominates linear estimators.

text: √ n) that is obtained by usual Rademacher
complexity analysis. This discrepancy is induced by the non-convex geometry of the model and the noisy gradient descent used for neural network training provably reaches a near global optimal solution even though the loss landscape is highly non-convex. Although the noisy gradient descent does not employ any explicit or implicit sparsity inducing regularization, it shows a preferable generalization performance that dominates linear estimators.

heading: 1 INTRODUCTION
text: In the deep learning theory literature, clarifying the mechanism by which deep learning can outperform shallow approaches has been gathering most attention for a long time. In particular, it is quite important to show that a tractable algorithm for deep learning can provably achieve a better generalization performance than shallow methods. Towards that goal, we study the rate of convergence of excess risk of both deep and shallow methods in a setting of a nonparametric regression problem. One of the difficulties to show generalization ability of deep learning with certain optimization methods is that the solution is likely to be stacked in a bad local minimum, which prevents us to show its preferable performances. Recent studies tackled this problem by considering optimization on overparameterized networks as in neural tangent kernel (NTK) (Jacot et al., 2018; Du et al., 2019a) and mean field analysis (Nitanda & Suzuki, 2017; Chizat & Bach, 2018; Rotskoff & Vanden-Eijnden, 2018; 2019; Mei et al., 2018; 2019), or analyzing the noisy gradient descent such as stochastic gradient Langevin dynamics (SGLD) (Welling & Teh, 2011; Raginsky et al., 2017; Erdogdu et al., 2018).
The NTK analysis deals with a relatively large scale initialization so that the model is well approximated by the tangent space at the initial solution, and eventually, all analyses can be reduced to those of kernel methods (Jacot et al., 2018; Du et al., 2019b; Allen-Zhu et al., 2019; Du et al., 2019a; Arora et al., 2019; Cao & Gu, 2019; Zou et al., 2020). Although this regime is useful to show its global
convergence, the obtained estimator looses large advantage of deep learning approaches because the estimation ability is reduced to the corresponding kernel methods. To overcome this issue, there are several “beyond-kernel” type analyses. For example, Allen-Zhu & Li (2019; 2020) showed benefit of depth by analyzing ResNet type networks. Li et al. (2020) showed global optimality of gradient descent by reducing the optimization problem to a tensor decomposition problem for a specific regression problem, and showed the “ideal” estimator on a linear model has worse dependency on the input dimensionality. Bai & Lee (2020) considered a second order Taylor expansion and showed that the sample complexity of deep approaches has better dependency on the input dimensionality than kernel methods. Chen et al. (2020) also derived a similar conclusion by considering a hierarchical representation. The analyses mentioned above actually show some superiority of deep learning, but all of these bounds are essentially Ω(1/ √ n) where n is the sample size, which is not optimal for regression problems with squared loss (Caponnetto & de Vito, 2007). The reason why only such a sub-optimal rate is considered is that the target of their analyses is mostly the Rademacher complexity of the set in which estimators exist for bounding the generalization gap. However, to derive a tight excess risk bound instead of the generalization gap, we need to evaluate so called local Rademacher complexity (Mendelson, 2002; Bartlett et al., 2005; Koltchinskii, 2006) (see Eq. (2) for the definition of excess risk). Moreover, some of the existing analyses should change the target function class as the sample size n increases, for example, the input dimensionality is increased against the sample size, which makes it difficult to see how the rate of convergence is affected by the choice of estimators.
Another promising approach is the mean field analysis. There are also some work that showed superiority of deep learning against kernel methods. Ghorbani et al. (2019) showed that, when the dimensionality d of input is polynomially increasing with respect to n, the kernel methods is outperformed by neural network approaches. Although the situation of increasing d explains well the modern high dimensional situations, this setting blurs the rate of convergence. Actually, we can show the superiority of deep learning even in a fixed dimension setting.
There are several studies about approximation abilities of deep and shallow models. Ghorbani et al. (2020) showed adaptivity of kernel methods to the intrinsic dimensionality in terms of approximation error and discuss difference between deep and kernel methods. Yehudai & Shamir (2019) showed that the random feature method requires exponentially large number of nodes against the input dimension to obtain a good approximation for a single neuron target function. These are only for approximation errors and estimation errors are not compared.
Recently, the superiority of deep learning against kernel methods has been discussed also in the nonparametric statistics literature where the minimax optimality of deep learning in terms of excess risk is shown. Especially, it is shown that deep learning achieves better rate of convergence than linear estimators in several settings (Schmidt-Hieber, 2020; Suzuki, 2019; Imaizumi & Fukumizu, 2019; Suzuki & Nitanda, 2019; Hayakawa & Suzuki, 2020). Here, the linear estimators are a general class of estimators that includes kernel ridge regression, k-NN regression and Nadaraya-Watson estimator. Although these analyses give clear statistical characterization on estimation ability of deep learning, they are not compatible with tractable optimization algorithms.
In this paper, we give a theoretical analysis that unifies these analyses and shows the superiority of a deep learning method trained by a tractable noisy gradient descent algorithm. We evaluate the excess risks of the deep learning approach and linear estimators in a nonparametric regression setting, and show that the minimax optimal convergence rate of the linear estimators can be dominated by the noisy gradient descent on neural networks. In our analysis, the model is fixed and no explicit sparse regularization is employed. Our contributions can be summarized as follows:
• A refined analysis of excess risks for a fixed model with a fixed input dimension is given to compare deep and shallow estimators. Although several studies pointed out the curse of dimensionality is a key factor that separates shallow and deep approaches, we point out that such a separation appears in a rather low dimensional setting, and more importantly, the non-convexity of the model essentially makes the two regimes different.
• A lower bound of the excess risk which is valid for any linear estimator is derived. The analysis is considerably general because the class of linear estimators includes kernel ridge regression with any kernel and thus it also includes estimators in the NTK regime. • All derived convergence rate is a fast learning rate that is faster than O(1/ √ n). We show that
simple noisy gradient descent on a sufficiently wide two-layer neural network achieves a fast
learning rate by using a fact that the solution converges to a Bayes estimator with a Gaussian process prior, and the derived convergence rate can be faster than that of linear estimators. This is much different from such existing work that compared only coefficients with the same rate of convergence with respect to the sample size n.
Other related work Bach (2017) analyzed the model capacity of neural networks and its corresponding reproducing kernel Hilbert space (RKHS), and showed that the RKHS is much larger than the neural network model. However, separation of the estimation abilities between shallow and deep is not proven. Moreover, the analyzed algorithm is basically the Frank-Wolfe type method which is not typically used in practical deep learning. The same technique is also employed by Barron (1993). The Frank-Wolfe algorithm is a kind of sparsity inducing algorithm that is effective for estimating a function in a model with an L1-norm constraint. It has been shown that explicit or implicit sparse regularization such as L1-regularization is beneficial to obtain better performances of deep learning under certain situations (Chizat & Bach, 2020; Chizat, 2019; Gunasekar et al., 2018; Woodworth et al., 2020; Klusowski & Barron, 2016). For example, E et al. (2019b;a) showed that the approximation error of a linear model suffers from the curse of dimensionality in a setting where the target function is in the Barron class (Barron, 1993), and showed anL1-type regularization avoids the curse of dimensionality. However, our analysis goes in a different direction where a sparse regularization is not required.

heading: 2 PROBLEM SETTING AND MODEL
text: In this section, we give the problem setting and notations that will be used in the theoretical analysis. We consider the standard nonparametric regression problem where data are generated from the following model for an unknown true function fo : Rd → R:
yi = f o(xi) + i (i = 1, . . . , n), (1)
where xi is independently identically distributed from PX whose support is included in Ω = [0, 1]d, and i is an observation noise that is independent of xi and satisfies E[ i] = 0 and i ∈ [−U,U ] almost surely. The n i.i.d. observations are denoted by Dn = (xi, yi)ni=1. We want to estimate the true function fo through the training data Dn. To achieve this purpose, we employ the squared loss `(y, f(x)) = (y − f(x))2 and accordingly we define the expected and empirical risks as L(f) := EY,X [`(Y, f(X))] and L̂(f) := 1n ∑n i=1 `(yi, f(xi)) respectively. Throughout this paper, we are interested in the excess (expected) risk of an estimator f̂ defined by
(Excess risk) L(f̂)− inf f :measurable L(f). (2)
Since the loss function ` is the squared loss, the infimum of inff :measurable L(f) is achieved by fo: inff :measurable L(f) = L(fo). The population L2(PX)-norm is denoted by ‖f‖L2(PX) :=√
EX∼PX [f(X) 2] and the sup-norm on the support of PX is denoted by ‖f‖∞ := supx∈supp(PX) |f(x)|. We can easily check that for an estimator f̂ , the L2-distance ‖f̂ − f o‖2L2(PX) between the estimator f̂ and the true function fo is identical to the excess risk: L(f̂) − L(fo) = ‖f̂ − fo‖2L2(PX). Note that the excess risk is different from the generalization gap L(f̂) − L̂(f̂). Indeed, the generalization gap typically converges with the rate of O(1/ √ n) which is optimal in a typical setting (Mohri et al., 2012). On the other hand, the excess risk can be faster than O(1/ √ n), which is known as a fast learning rate (Mendelson, 2002; Bartlett et al., 2005; Koltchinskii, 2006; Giné & Koltchinskii, 2006).

heading: 2.1 MODEL OF TRUE FUNCTIONS
text: To analyze the excess risk, we need to specify a function class (in other words, model) in which the true function fo is included. In this paper, we only consider a two layer neural network model, whereas the techniques adapted in this paper can be directly extended to deeper neural network models. We consider a teacher-student setting, that is, the true function fo can be represented by a neural network defined as follows. For w ∈ R, let w̄ be a “clipping” of w defined as w̄ :=
R × tanh(w/R) where R ≥ 1 is a fixed constant, and let [x; 1] := [x>, 1]> for x ∈ Rd. Then, the teacher network is given by
fW (x) = ∑∞ m=1 amw̄2,mσm(w > 1,m[x; 1]),
where w1,m ∈ Rd+1 and w2,m ∈ R (m ∈ N) are the trainable parameters (where W = (w1,m, w2,m) ∞ m=1), am ∈ R (m ∈ N) is a fixed scaling parameter, and σm : R→ R is an activation function for the m-th node. The reason why we applied the clipping operation to the parameter of the second layer is just for a technical reason to ensure convergence of Langevin dynamics. The dynamics is bounded in high probability in practical situations and the boundedness condition would be removed if further theoretical development of infinite dimensional Langevin dynamics would be achieved.
Let H be a set of parameters W such that its squared norm is bounded: H := {W = (w1,m, w2,m) ∞ m=1 | ∑∞ m=1(‖w1,m‖2 + w22,m) < ∞}. Define ‖W‖H := [ ∑∞ m=1(‖w1,m‖2 + w22,m)] 1/2 for W ∈ H. Let (µm)∞m=1 be a regularization parameter such that µm ↘ 0. Accordingly
we defineHγ := {W ∈ H | ‖W‖Hγ <∞} where ‖W‖Hγ := [ ∑∞ m=1 µ −γ m (‖w1,m‖2 + w22,m)]1/2 for a given 0 < γ. Throughout this paper, we analyze an estimation problem in which the true function is included in the following model:
Fγ = {fW |W ∈ Hγ , ‖W‖Hγ ≤ 1}.
This is basically two layer neural network with infinite width. As assumed later, am is assumed to decrease as m → ∞. Its decreasing rate controls the capacity of the model. If the first layer parameters (w1,m)m are fixed, this model can be regarded as a variant of the unit ball of some reproducing kernel Hilbert space (RKHS) with basis functions amσm(w>1,m[x; 1]). However, since the first layer (w1,m) is also trainable, there appears significant difference between deep and kernel approaches. The Barron class (Barron, 1993; E et al., 2019b) is relevant to this function class. Indeed, it is defined as the convex hull of w2σ(w>1 [x; 1]) with norm constraints on (w1, w2) where σ is an activation function. On the other hand, we will put an explicit decay rate on am and the parameter W has an L2-norm constraint, which makes the model Fγ smaller than the Barron class.

heading: 3 ESTIMATORS
text: We consider two classes of estimators and discuss their differences: linear estimators and deep learning estimator with noisy gradient descent (NGD).
Linear estimator A class of linear estimators, which we consider as a representative of “shallow” learning approach, consists of all estimators that have the following form:
f̂(x) = ∑n i=1 yiϕi(x1, . . . , xn, x).
Here, (ϕi)ni=1 can be any measurable function (and L2(PX)-integrable so that the excess risk can be defined). Thus, they could be selected as the “optimal” one so that the corresponding linear estimator minimizes the worst case excess risk. Even if we chose such an optimal one, the worst case excess risk should be lower bounded by our lower bound given in Theorem 1. It should be noted that the linear estimator does not necessarily imply “linear model.” The most relevant linear estimator in the machine learning literature is the kernel ridge regression: f̂(x) = Y >(KX + λI)−1k(x) where KX = (k(xi, xj))ni,j=1 ∈ Rn×n, k(x) = [k(x, x1), . . . , k(x, xn)]> ∈ Rn and Y = [y1, . . . , yn]
> ∈ Rn for a kernel function k : Rd × Rd → R. Therefore, the ridge regression estimator in the NTK regime or the random feature model is also included in the class of linear estimators. The solution obtained in the early stopping criteria instead of regularization in the NTK regime under the squared loss is also included in the linear estimators. Other examples include the k-NN estimator and the Nadaraya-Watson estimator. All of them do not train the basis function in a nonlinear way, which makes difference from the deep learning approach. In the nonparametric statistics literature, linear estimators have been studied for estimating a wavelet series model. Donoho et al. (1990; 1996) have shown that a wavelet shrinkage estimator can outperform any linear estimator by showing suboptimality of linear estimators. Suzuki (2019) utilized such an argument to show superiority of deep learning but did not present any tractable optimization algorithm.
Noisy Gradient Descent with regularization As for the neural network approach, we consider a noisy gradient descent algorithm. Basically, we minimize the following regularized empirical risk:
L̂(fW ) + λ2 ‖W‖ 2 H1 .
Here, we employ H1-norm as the regularizer. We note that the constant γ controls the relative complexity of the true function fo compared to the typical solution obtained under the regularization. Here, we define a linear operator A as λ‖W‖H1 = W>AW , that is, AW = (λµ−1m w1,m, λµ −1 m w2,m) ∞ m=1. The regularized empirical risk can be minimized
by noisy gradient descent as Wk+1 = Wk − η∇(L̂(fWk) + λ2 ‖Wk‖ 2 H1) + √ 2η β ξk, where η > 0 is a step size and ξk = (ξk,(1,m), ξk,(2,m))∞m=1 is an infinite-dimensional Gaussian noise, i.e., ξk,(1,m) and ξk,(2,m) are independently identically distributed from the standard normal distribution (Da Prato & Zabczyk, 1996). Here, ∇L̂(fW ) = 1n ∑n i=1 2(fW (xi) − yi)(w̄2,mam[xi; 1]σ ′ m(w > 1,m[xi; 1]), am tanh ′(w2,m/R)σm(w > 1,m[xi; 1])) ∞ m=1. However, since ∇‖Wk−1‖2H1 is unbounded which makes it difficult to show convergence, we employ the semiimplicit Euler scheme defined by
Wk+1 =Wk−η∇L̂(fWk)−ηAWk+1+ √ 2η β ξk ⇔Wk+1 =Sη ( Wk−η∇L̂(fWk)+ √ 2η β ξk ) , (3)
where Sη := (I+ηA)−1. It is easy to check that this is equivalent to the following update rule: Wk = Wk−1 − η ( ∇L̂(fWk−1) + SηAWk−1 + √ 2η β ξk−1 ) . Therefore, the implicit Euler scheme can be seen as a naive noisy gradient descent for minimizing the empirical risk with a slightly modified ridge regularization. This can be interpreted as a discrete time approximation of the following infinite dimensional Langevin dynamics:
dWt = −∇(L̂(fWt) + λ2 ‖Wt‖ 2 H1)dt+ √ 2/βdξt, (4)
where (ξt)t≥0 is the so-called cylindrical Brownian motion (see Da Prato & Zabczyk (1996) for the details). Its application and analysis for machine learning problems with non-convex objectives have been recently studied by, for example, Muzellec et al. (2020); Suzuki (2020).
The above mentioned algorithm is executed on an infinite dimensional parameter space. In practice, we should deal with a finite width network. To do so, we approximate the solution by a finite dimensional one: W (M) = (w1,m, w2,m)Mm=1 where M corresponds to the width of the network. We identify W (M) to the “zero-padded” infinite dimensional one, W = (w1,m, w2,m)∞m=1 with w1,m = 0 and w2,m = 0 for all m > M . Accordingly, we use the same notation fW (M) to indicate fW with zero padded vector W . Then, the finite dimensional version of the update rule is given by W (M) k+1 = S (M) η ( W
(M) k − η∇L̂(fW (M)k ) + √ 2η β ξ (M) k ) , where ξ(M)k is the Gaussian noise vector
obtained by projecting ξk to the first M components and S (M) η is also obtained in a similar way.

heading: 4 CONVERGENCE RATE OF ESTIMATORS
text: In this section, we present the excess risk bounds for linear estimators and the deep learning estimator. As for the linear estimators, we give its lower bound while we give an upper bound for the deep learning approach. To obtain the result, we setup some assumptions on the model. Assumption 1.
(i) There exists a constant cµ such that µm ≤ cµm−2 (m ∈ N). (ii) There exists α1 > 1/2 such that am ≤ µα1m (m ∈ N).
(iii) The activation functions (σm)m is bounded as ‖σm‖∞ ≤ 1. Moreover, they are three times differentiable and their derivatives upto third order differentiation are uniformly bounded: ∃Cσ such that ‖σm‖1,3 := max{‖σ′m‖∞, ‖σ′′m‖∞, ‖σ′′′m‖∞} ≤ Cσ (∀m ∈ N).
The first assumption (i) controls the strength of the regularization, and combined with the second assumption (ii) and definition of the model Fγ , complexity of the model is controlled. If α1 and γ are large, the model is less complicated. Indeed, the convergence rate of the excess risk becomes
faster if these parameters are large as seen later. The decay rate µm ≤ cµm−2 can be generalized as m−p with p > 1 but we employ this setting just for a technical simplicity for ensuring convergence of the Langevin dynamics. The third assumption (iii) can be satisfied by several activation functions such as the sigmoid function and the hyperbolic tangent. The assumption ‖σm‖∞ ≤ 1 could be replaced by another one like ‖σm‖∞ ≤ C, but we fix this scaling for simple presentation.

heading: 4.1 MINIMAX LOWER BOUND FOR LINEAR ESTIMATORS
text: Here, we analyze a lower bound of excess risk of linear estimators, and eventually we show that any linear estimator suffers from curse of dimensionality. To rigorously show that, we consider the following minimax excess risk over the class of linear estimators:
Rlin(Fγ) := inf f̂ :linear sup fo∈Fγ EDn [‖f̂ − fo‖2L2(PX)],
where inf is taken over all linear estimators and EDn [·] is taken with respect to the training data Dn. This expresses the best achievable worst case error over the class of linear estimators to estimate a function in Fγ . To evaluate it, we additionally assume the following condition. Assumption 2. We assume that µm = m−2 and am = µα1m (m ∈ N) (and hence cµ = 1). There exists a monotonically decreasing sequence (bm)∞m=1 and s ≥ 3 such that bm = µα2m (∀m) with α2 > γ/2 and σm(u) = bsmσ(b −1 m u) (u ∈ R) where σ is the sigmoid function: σ(u) = 1/(1+e−u).
Intuitively, the parameter s controls the “resolution” of each basis function σm, and the relation between parameter α1 and α2 controls the magnitude of coefficient for each basis σm. Note that the condition s ≥ 3 ensures ‖σm‖1,3 is uniformly bounded and 0 < bm ≤ 1 ensures ‖σm‖∞ ≤ 1. Our main strategy to obtain the lower bound is to make use of the so-called convex-hull argument. That is, it is known that, for a function class F , the minimax risk R(F) over a class of linear estimators is identical to that for the convex hull of F (Hayakawa & Suzuki, 2020; Donoho et al., 1990):
Rlin(F) = Rlin(conv(F)), where conv(F) = { ∑N i=1 λifi | fi ∈ F , ∑N i=1 λi = 1, λi ≥ 0, N ∈ N} and conv(·) is the closure of conv(·) with respect to L2(PX)-norm. Intuitively, since the linear estimator is linear to the observations (yi)ni=1 of outputs, a simple application of Jensen’s inequality yields that its worst case error on the convex hull of the function class F does not increase compared with that on the original one F (see Hayakawa & Suzuki (2020) for the details). This indicates that the linear estimators cannot distinguish the original hypothesis class F and its convex hull. Therefore, if the class F is highly non-convex, then the linear estimators suffer from much slower convergence rate because its convex hull conv(F) becomes much “fatter” than the original one F . To make use of this argument, for each sample size n, we pick up appropriate mn and consider a subset generated by the basis function σmn , i.e., F (n) γ := {amnw̄2,mnσm(w>1,mn [x; 1]) ∈ Fγ}. By applying the convex hull argument to this set, we obtain the relation Rlin(Fγ) ≥ Rlin(F (n)γ ) = Rlin(conv(F (n)γ )). Since F (n)γ is highly non-convex, its convex hull conv(F (n)γ ) is much larger than the original set F (n)γ and thus the minimax risk over the linear estimators would be much larger than that over all estimators including deep learning. More intuitively, linear estimators do not adaptively select the basis functions and thus they should prepare redundantly large class of basis functions to approximate functions in the target function class. The following theorem gives the lower bound of the minimax optimal excess risk over the class of linear estimators. Theorem 1. Suppose that Var( ) > 0, PX is the uniform distribution on [0, 1]d, and Assumption 2 is satisfied. Let β̃ = α1+(s+1)α2α2−γ/2 . Then for arbitrary small κ ′ > 0, we have that
Rlin(Fγ) & n − 2β̃+d 2β̃+2dn−κ ′ . (5)
The proof is in Appendix A. We utilized the Irie-Miyake integral representation (Irie & Miyake, 1988; Hornik et al., 1990) to show there exists a “complicated” function in the convex hull, and then we adopted the technique of Zhang et al. (2002) to show the lower bound. The lower bound is characterized by the decaying rate (α1) of am relative to that (α2) of the scaling factor bm. Indeed, the faster am decays with increasing m, the faster the rate of the minimax lower bound becomes.
We can see that the minimax rate of linear estimators is quite sensitive to the dimension d. Actually, for relatively high dimensional settings, this lower bound becomes close to a slow rate Ω(1/ √ n), which corresponds to the curse of dimensionality.
It has been pointed out that the sample complexity of kernel methods suffers from the curse of dimensionality while deep learning can avoid that by a tractable algorithms (e.g., Ghorbani et al. (2019); Bach (2017)). Among them, Ghorbani et al. (2019) showed that if the dimensionality d is polynomial against n, then the excess risk of kernel methods is bounded away from 0 for all n. On the other hand, our analysis can be applied to any linear estimator including kernel methods, and it shows that even if the dimensionality d is fixed, the convergence rate of their excess risk suffers from the curse of dimensionality. This can be accomplished thanks to a careful analysis of the rate of convergence. Bach (2017) derived an upper bound of the Rademacher complexity of the unit ball of the RKHS corresponding to a neural network model. However, it is just an upper bound and there is still a large gap from excess risk estimates. Allen-Zhu & Li (2019; 2020); Bai & Lee (2020); Chen et al. (2020) also analyzed a lower bound of sample complexity of kernel methods. However, their lower bound is not for the excess risk of the squared loss. Eventually, the sample complexities of all methods including deep learning take a form of O(C/ √ n) and dependency of coefficient C to the dimensionality or other factors such as magnitude of residual components is compared. On the other hand, our lower bound properly involves the properties of squared loss such as strong convexity and smoothness and the bound shows the curse of dimensionality occurs even in the rate of convergence instead of just the coefficient. Finally, we would like to point out that several existing work (e.g., Ghorbani et al. (2019); Allen-Zhu & Li (2019)) considered a situation where the target function class changes as the sample size n increases. However, our analysis reveals that separation between deep and shallow occurs even if the target function class Fγ is fixed.

heading: 4.2 UPPER BOUND FOR DEEP LEARNING
text: Here, we analyze the excess risk of deep learning trained by NGD and its algorithmic convergence rate. Our analysis heavily relies on the weak convergence of the discrete time gradient Langevin dynamics to the stationary distribution of the continuous time one (Eq. (4)). Under some assumptions, the continuous time dynamics has a stationary distribution (Da Prato & Zabczyk, 1992; Maslowski, 1989; Sowers, 1992; Jacquot & Royer, 1995; Shardlow, 1999; Hairer, 2002). If we denote the probability measure onH corresponding to the stationary distribution by π∞, then it is given by
dπ∞ dνβ (W ) ∝ exp(−βL̂(fW )),
where νβ is the Gaussian measure in H with mean 0 and covariance (βA)−1 (see Da Prato & Zabczyk (1996) for the rigorous definition of the Gaussian measure on a Hilbert space). Remarkably, this can be seen as the Bayes posterior for a prior distribution νβ and a “log-likelihood” function exp(−βL̂(W )). Through this view point, we can obtain an excess risk bound of the solution Wk. The proofs of all theorems in this section are in Appendix B.
Under Assumption 1, the distribution ofWk derived by the discrete time gradient Langevin synamics satisfies the following weak convergence property to the stationary distribution π∞. This convergence rate analysis depends on the techniques by Bréhier & Kopec (2016); Muzellec et al. (2020). Proposition 1. Assume Assumption 1 holds and β > η. Then, there exist spectral gaps Λ∗η and Λ∗0 (defined in Eq. (10) of Appendix B.1) and a constant C0 such that, for any 0 < a < 1/4, the following convergence bound holds for almost sure observation Dn:
|EWk [L(fWk)|Dn]− EW∼π∞ [L(fW )|Dn]| ≤ C0 exp(−Λ∗ηηk) + C1 √ β
Λ∗0 η1/2−a =: Ξk, (6)
where C1 is a constant depending only on cµ, R, α1, Cσ, U, a (independent of η, k, β, λ, n).
This proposition indicates that the expected risk of Wk can be almost identical to that of the “Bayes posterior solution” obeying π∞ after sufficiently large iterations k with sufficiently small step size η even though L̂(fW ) is not convex. The definition of Λ∗η can be found in Eq. (10). We should note that its dependency on β is exponential. Thus, if we take β = Ω(n), then the computational cost until a sufficiently small error could be exponential with respect to the sample size n. The same convergence holds also for finite dimensional oneW (M)k with a modified stationary distribution. The
constants appearing in the bound are independent of the model size M (see the proof of Proposition 1 in Appendix B). In particular, the convergence can be guaranteed even ifW is infinite dimensional. This is quite different from usual finite dimensional analyses (Raginsky et al., 2017; Erdogdu et al., 2018; Xu et al., 2018) which requires exponential dependency on the dimension, but thanks to the regularization term, we can obtain the model size independent convergence rate. Xu et al. (2018) also analyzed a finite dimensional gradient Langevin dynamics and obtained a similar bound where O(η) appears in place of the second term η1/2−a which corresponds to time discretization error. In our setting the regularization term is ‖W‖2H1 = ∑ m(‖w1,m‖2 + w22,m)/µm with µm . m−2, but
if we employ ‖W‖2Hp/2 = ∑ m(‖w1,m‖2 + w22,m)/µ p/2 m for p > 1, then the time discretization error term would be modified to η(p−1)/p−a (Andersson et al., 2016). We can interpret the finite dimensional setting as the limit of p → ∞ which leads to η(p−1)/p → η that recovers the finite dimensional result (O(η)) as shown by Xu et al. (2018).
In addition to the above algorithmic convergence, we also have the following convergence rate for the excess risk bound of the finite dimensional solution W (M)k .
Theorem 2. Assume Assumption 1 holds, assume η < β ≤ min{n/(2U2), n}, and 0 < γ < 1/2 + α1. Then, if the width satisfies M ≥ min { λ1/4γ(α1+1)β1/2γ , λ−1/2(α1+1), n1/2γ } , the expected excess risk of Wk is bounded as
EDn [ E W
(M) k
[‖f W
(M) k
−fo‖2L2(PX)|Dn] ] ≤C max { (λβ) 1/γ 1+1/2γ n− 1 1+1/2γ, λ − 1 2(α1+1) β−1, λ γ 1+α1 } +Ξk,
where C is a constant independent of n, β, λ, η, k. In particular, if we set β = min{n/(2U2), n} and λ = β−1, then for M ≥ n1/2(α1+1), we obtain
EDn [ E W
(M) k
[‖f W
(M) k
− fo‖2L2(PX)|Dn] ] . n− γ α1+1 + Ξk.
In addition to this theorem, if we further assume Assumption 2, we obtain a refined bound as follows. Corollary 1. Assume Assumptions 1 and 2 hold and η < β, and let β = min{n/(2U2), n} and λ = β−1. Suppose that there exists 0 ≤ q ≤ s − 3 such that 0 < γ < 1/2 + α1 + qα2. Then, the excess risk bound of W (M)k for M ≥ n1/2(α1+qα2+1) can be refined as
EDn [ E W
(M) k
[‖f W
(M) k
− fo‖2L2(PX)|Dn] ] . n− γ α1+qα2+1 + Ξk. (7)
These theorem and corollary shows that the tractable NGD algorithm achieves a fast convergence rate of the excess risk bound. Indeed, if q is chosen so that γ > (α1+qα2+1)/2, then the excess risk bound converges faster thanO(1/ √ n). Remarkably, the convergence rate is not affected by the input dimension d, which makes discrepancy from linear estimators. The bound of Theorem 2 is tightest when γ is close to 1/2 + α1 (γ ≈ 1/2 + α1 + 3α2 for Corollary 1), and a smaller γ yields looser bound. This relation between γ and α1 reflects misspecification of the “prior” distribution. When γ is small, the regularization λ‖W‖2H1 is not strong enough so that the variance of the posterior distribution becomes unnecessarily large for estimating the true function fo ∈ Fγ . Therefore, the best achievable bound can be obtained when the regularization is correctly specified. The analysis of fast rate is in contrast to some existing work (Allen-Zhu & Li, 2019; 2020; Li et al., 2020; Bai & Lee, 2020) that basically evaluated the Rademacher complexity. This is because we essentially evaluated a local Rademacher complexity instead.

heading: 4.3 COMPARISON BETWEEN LINEAR ESTIMATORS AND DEEP LEARNING
text: Here, we compare the convergence rate of excess risks between the linear estimators and the neural network method trained by NGD using the bounds obtained in Theorem 1 and Corollary 1 respectively. We write the lower bound (5) of the minimax excess risk of linear estimators as R∗lin and the excess risk of the neural network approach (7) asR∗NN. To make the discussion concise, we consider a specific situation where s = 3, α1 = γ = 14α2. In this case, β̃ = 17/3 ≈ 5.667, which gives
R∗lin & n − ( 1+ d 2β̃+d )−1 n−κ ′ ≈ n − ( 1+ d 11.3+d )−1 n−κ ′ .
On the other hand, by setting q = 0, we have
R∗NN . n − α1α1+1 = n
− (
1+ 1 α1 )−1 .
Thus, as long as α1 > 11.3/d+ 1 ≈ 2β̃/d+ 1, we have that
R∗lin & R ∗ NN, and limn→∞ R∗NN R∗lin = 0.
In particular, as d gets larger, R∗lin approaches to Ω(n −1/2) while R∗NN is not affected by d and it gets close to O(n−1) as α1 gets larger. Moreover, the inequality α1 > 11.3/d + 1 can be satisfied by a relatively low dimensional setting; for example, d = 10 is sufficient when α1 = 3. As α1 becomes large, the model becomes “simpler” because (am)∞m=1 decays faster. However, the linear estimators cannot take advantage of this information whereas deep learning can. From the convex hull argument, this discrepancy stems from the non-convexity of the model. We also note that the superiority of deep learning is shown without sparse regularization while several existing work showed favorable estimation property of deep learning though sparsity inducing regularization (Bach, 2017; Chizat, 2019; Hayakawa & Suzuki, 2020). However, our analysis indicates that sparse regularization is not necessarily as long as the model has non-convex geometry, i.e., sparsity is just one sufficient condition for non-convexity but not a necessarily condition. The parameter setting above is just a sufficient condition and the lower bound R∗lin would not be tight. The superiority of deep learning would hold in much wider situations.

heading: 5 CONCLUSION
text: In this paper, we studied excess risks of linear estimators, as a representative of shallow methods, and a neural network estimator trained by a noisy gradient descent where the model is fixed and no sparsity inducing regularization is imposed. Our analysis revealed that deep learning can outperform any linear estimator even for a relatively low dimensional setting. Essentially, non-convexity of the model induces this difference and the curse of dimensionality for linear estimators is a consequence of a fact that the geometry of the model becomes more “non-convex” as the dimension of input gets higher. All derived bounds are fast rate because the analyses are about the excess risk with the squared loss, which made it possible to compare the rate of convergence. The fast learning rate of the deep learning approach is derived through the fact that the noisy gradient descent behaves like a Bayes estimator with model size independent convergence rate.

heading: ACKNOWLEDGMENTS
text: TS was partially supported by JSPS Kakenhi (18K19793, 18H03201, and 20H00576), Japan Digital Design and JST-CREST.

heading: A PROOF OF THEOREM 1
text: We basically combine the “convex hull argument” and the minimax optimal rate analysis for linear estimators developed by Zhang et al. (2002).
Zhang et al. (2002) essentially showed the following statement in their Theorem 1. Proposition 2 (Theorem 1 of Zhang et al. (2002)). Let µ be the Lebesgue measure. Suppose that the space Ω has even partition A such that |A| = 2K for an integer K ∈ N, each A has equivalent measure µ(A) = 2−K for all A ∈ A, andA is indeed a partition of Ω, i.e., ∪A∈A = Ω, A∩A′ = ∅ for A,A′ ∈ Ω and A 6= A′. Then, if K is chosen as n−γ1 ≤ 2−K ≤ n−γ2 for constants γ1, γ2 > 0 that are independent of n, then there exists an event E such that, for a constant C ′ > 0,
P (E) ≥ 1 + o(1) and |{xi | xi ∈ A (i ∈ {1, . . . , n})}| ≤ C ′n/2K (∀A ∈ A). Moreover, suppose that, for a class F◦ of functions on Ω, there exists ∆ > 0 that satisfies the following conditions:
1. There exists F > 0 such that, for any A ∈ A, there exists g ∈ F◦ that satisfies g(x) ≥ 1 2∆F for all x ∈ A,
2. There exists K ′ and C ′′ > 0 such that 1n ∑n i=1 g(xi)
2 ≤ C ′′∆22−K′ for any g ∈ F◦ on the event E .
Then, there exists a constant F1 such that at least one of the following inequalities holds:
F 2 4F1C ′′ 2K ′ n ≤ Rlin(F◦), (8a) F 3 32 ∆22−K ≤ Rlin(F◦), (8b)
for sufficiently large n.
Before we show the main assertion, we prepare some additional lemmas. For a sigmoid function σ, let F̃ (σ)C,τ := {x ∈ Rd 7→ aσ(τ(w>x+ b))) | |a| ≤ 2C, ‖w‖ ≤ 1, |b| ≤ 2 (a, b ∈ R, w ∈ Rd)} for C > 0, τ > 0. Lemma 1. Let ψ(x) = 12 (σ(x + 1) − σ(x − 1)) and ψ̂ be its Fourier transform: ψ̂(ω) := (2π)−1 ∫ e−iωxψ(x)dx. Let h > 0 and Dw > 0. Then, by setting τ = h−1(2 √ d + 1)Dw and C = (2 √ d+1)Dw
πh|ψ̂(1)| , the Gaussian RBF kernel can be approximated by
inf ǧ∈conv(F̃(σ)C,τ ) sup x∈[0,1]d
∣∣∣∣ǧ(x)− exp(−‖x− c‖22h2 )∣∣∣∣
≤ 4 |2πψ̂(1)|
[ CdD 2(d−2) w exp(−D2w/2) + exp(−Dw) ] for any c ∈ [0, 1]d, where Cd is a constant depending only on d. In particular, the right hand side is O(exp(−nκ)) if Dw = nκ.
Proof of Lemma 1. Let ψh(x) = ψ(h−1x). Suppose that, for f ∈ L1(Rd), its Fourier transform f̂(ω) = (2π)−d ∫ e−iω
>xf(x)dx (ω ∈ Rd) gives∫ Rd exp(iw>x)f̂(w)dw = f(x),
for every x ∈ Rd1. Then the Irie-Miyake itegral representation (Irie & Miyake (1988); see also the proof of Theorem 3.1 in Hornik et al. (1990)) gives
f(x) = ∫ a∈Rd ∫ b∈R ψ(a>x+ b)dν(a, b) (a.e.),
1If f̂ is integrable, this inversion formula holds for almost every x ∈ Rd (Rudin, 1987). However, we assume a stronger condition that it holds for every x ∈ Rd.
where dν(a, b) is given by
dν(a, b) = Re
( |ω|de−iwb
2πψ̂(ω)
) f̂(wa)dadb
for any ω 6= 0. Since the characteristic function of the multivariate normal distribution gives that∫ Rd exp(iw>(x− c)) √ h2d (2π)d exp ( −h 2‖w‖2 2 ) ︸ ︷︷ ︸
=f̂(w)
dw = exp ( −‖x− c‖ 2
2h2
) =: f(x) (∀x ∈ Rd),
we have that
exp ( −‖x− c‖ 2
2h2
) =
∫ a∈Rd ∫ b∈R ψh(a >(x− c) + b)Re ( e−iwb 2πψ̂h(ω) )√ |ωh|2d (2π)d exp ( − (ωh) 2‖a‖2 2 ) dadb,
for all x ∈ Rd. Since ψh(·) = ψ(h−1·) and ψ̂h(·) = hψ̂(h·) by its definition, the right hand side is equivalent to∫
a∈Rd ∫ b∈R ψ(h−1[a>(x− c) + b])Re ( e−iwb 2πhψ̂(hω) )√ |ωh|2d (2π)d exp ( − (ωh) 2‖a‖2 2 ) dadb.
Here, we set ω = h−1. Let Nσ2 be the probability measure corresponding to the multivariate normal with mean 0 and covariance σ2I, and let AD := {w ∈ Rd | ‖w‖ ≤ D}. Let Da > 0 and Db = Da( √ 2d+ 1), and define
fDa(x) := 1
2DbN1(ADa) ∫ ‖a‖≤Da,|b|≤Db ψ(h−1[a>(x− c) + b])Re ( e−ib/h 2πhψ̂(1) ) ×√
1
(2π)d exp
( −‖a‖ 2
2
) dadb.
Then, we can see that, for any x ∈ [0, 1]d, it holds that∣∣∣∣ 12DbN1(ADa)f(x)− fDa(x) ∣∣∣∣
≤ 1 2DbN1(ADa)|2πhψ̂(1)|
[ N1(A c Da) ∫ 2 exp(−h−1|x|)dx+ ∫ |b|>Db 2 exp(−[h−1(|b| − 2 √ dDa)])db ]
≤ 1 2DbN1(ADa)|2πhψ̂(1)|
[ 4hN1(A c Da) + 4h exp(−Da) ] = 4h
2DbN1(ADa)|2πhψ̂(1)|
[ CdD 2(d−2) a exp(−D2a/2) + exp(−Da) ] ,
where Cd > 0 is a constant depending on only d, and we used |a>(x − c) + b| ≥ |b| − |a>(x − c)| ≥ |b| − 2 √ dDa and ψ(x) ≤ 2 exp(−|x|). Note that if Da = nκ, then the right hand side is O(h exp(−nκ)). Therefore, since N1(ADa) ≤ 1, by setting τ = h−1Db, C = Dbπh|ψ̂(1)| , we have that
inf ǧ∈conv(F̃(σ)C,τ ) sup x∈[0,1]d
∣∣∣∣ǧ(x)− exp(−‖x− c‖22h2 )∣∣∣∣
≤ 4 |2πψ̂(1)|
[ CdD 2(d−2) a exp(−D2a/2) + exp(−Da) ] .
Hence, by rewriting Dw ← Da, we obtain the assertion. As noted above, the right hand is O(exp(−nκ)) if Da = nκ.
Proof of Theorem 1. For a sample size n, we fix mn which will be determined later and use Proposition 2 with F◦ = F (n)γ . If w2,mn = b √ µγmn/2 with |b| ≤ 1 and w1,m = µ γ/2 mn [u;−u>c]/( √ 2(d+ 1)) for u ∈ Rd such that ‖u‖ ≤ 1 and c ∈ [0, 1]d, then ‖(w1,mn , w2,mn)‖2 ≤ µγmn(1/2 + (1 + |u >c|2)/2(d + 1)) ≤ µγmn . Therefore, ϕ̃u,c(x) = amnw̄2,mnσmn(w > 1,mn [x; 1]) = µ α1 mn(bµ γ/2 mn / √ 2)µsα2mn σ ( µ −α2+γ/2 mn u >(x− c)/ √ 2(d+ 1) ) ∈
F (n)γ ⊂ Fγ for all b ∈ R with |b| ≤ 1, u ∈ Rd with ‖u‖ ≤ 1, and c ∈ [0, 1]d. In other words, µ α1+γ/2+sα2 mn (2C) −1F (σ)C,τ ⊂ F (n) γ for any C > 0 and τ = 1√ 2(d+1) µ −α2+γ/2 mn . Therefore, by setting C = ( √
2d + 1)Dw/(πh|ψ̂(1)|) for Dw > 0, Lemma 1 yields that for any c ∈ [0, 1]d and given h > 0, there exists g ∈ conv(F (n)γ ) such that∥∥∥∥∥∥µα1+γ/2+sα2mn ( 2( √ 2d+ 1)Dw πh|ψ̂(1)| )−1 exp ( −‖ · −c‖ 2 2h2 ) − g ∥∥∥∥∥∥ ∞
≤ µα1+γ/2+sα2mn
( 2( √
2d+ 1)Dw
πh|ψ̂(1)|
)−1 4
|2πψ̂(1)|
[ CdD 2(d−2) w exp(−D2w/2) + exp(−Dw) ] = µα1+γ/2+sα2mn h
( √
2d+ 1)Dw
[ CdD 2(d−2) a exp(−D2w/2) + exp(−Dw) ] .
We let Dw = nκ for any κ > 0 and choose µmn as τ ' µ −α2+γ/2 mn = Dwh −1 = h−1nκ. We write ∆ := µ α1+γ/2+sα2 mn (2C) −1 ' h α1+sα2+γ/2 α2−γ/2 +1 n −κ(α1+sα2+γ/2
α2−γ/2 +1). Then, it holds that∥∥∥∥∆ exp(−‖ · −c‖22h2 ) − g ∥∥∥∥ ∞ . ∆ exp(−nκ). (9)
Here, we set h as h = 2−k with a positive integer k. Accordingly, we define a partition A of Ω so that any A ∈ A can be represented as A = [2−kj1, 2−k(j1 + 1)] × · · · × [2−kjd, 2−k(jd + 1)] by non-negative integers 0 ≤ ji ≤ 2k − 1 (i = 1, . . . , d). Note that |A| = 2dk = h−d. For each A ∈ A, we define cA as cA = (2−k(j1 + 1/2), . . . , 2−k(jd + 1/2))> where (j1, . . . , jd) is a set of indexes that satisfies A = [2−kj1, 2−k(j1 + 1)] × · · · × [2−kjd, 2−k(jd + 1)]. For each A ∈ A, we define gA ∈ conv(F (n)γ ) as a function that satisfies Eq. (9) for c = cA.
Now, we apply Proposition 2 with F◦ = conv(F (n)γ ) and K = K ′ = dk. Let R∗ := Rlin(conv(F (n)γ )). First, we can see that there exits a constant F > 0 such that
gA(x) ≥ F∆ (∀x ∈ A),
where we used exp(−nκ) 1. Second, in the event E introduced in the statement of Proposition 2, there exists C such that |{i ∈ {1, . . . , n} | xi ∈ A′}| ≤ Cn/2−dk for all A′ ∈ A. In this case, we can check that
1
n n∑ i=1 [ ∆ exp ( −‖xi − cA‖ 2 2h2 )]2 . ∆2hd = ∆22−kd,
by the uniform continuity of the Gaussian RBF. Therefore, we also have
1
n n∑ i=1 gA(xi) 2 ≤ 2 n n∑ i=1 [ ∆ exp ( −‖xi − cA‖ 2 2h2 )]2 + c∆2 exp(−2nκ)
. ∆2(hd + exp(−2nκ)),
where c > 0 is a constant. Thus, as long as h is polynomial to n like h = Θ(n−a), the right hand side is O(∆2hd).
Now, if we write
β̃ = α1 + sα2 + γ/2
α2 − γ/2 + 1 = α1 + (s+ 1)α2 α2 − γ/2 ,
then we have ∆ ' hβ̃n−κβ̃ by its definition.
Here, we choose k as a maximum integer that satisfies F 3 32 ∆ 22−dk > R∗. In this situation, it holds that h2β̃+dn−2κβ̃ ' R∗.
Since Eq. (8b) is not satisfied, Eq. (8a) must hold, and hence we have
n−1h−d . R∗ ' h2β̃+dn−2κβ̃
⇒ h ' n− 1−2κβ̃ 2β̃+2d .
Therefore, we obtain that
R∗ & n − 2β̃+d 2β̃+2dn − 2κdβ̃ 2β̃+2d
≥ n− 2β̃+d 2β̃+2dn−κ ′ ,
by setting κ′ = κ 2dβ̃ 2β̃+2d . This gives the assertion.

heading: B PROOFS OF PROPOSITION 1, THEOREM 2 AND COROLLARY 1
text: Proposition 1, Theorem 2 and Corollary 1 can be shown by using Propositions 3 and 4 given in Appendix B.1 shown below.
Let TαW = (µαmw1,m, µ α mw2,m) ∞ m=1 for W = (w1,m, w2,m) ∞ m=1 for α > 0, and let us consider a model hW := fT−α/2W . Then, the training error can be rewritten as
L̂(fW ) = L̂(hTα/2W ).
For notational simplicity, we let L̂(W ) := L̂(fW ).
Let H(M) be {W (M) = (w1,m, w2,m)Mm=1 | w1,m ∈ Rd+1, w2,m ∈ R, 1 ≤ m ≤ M} and ι : H(M) → H be the zero padding of W (M), that is, ι(W (M)) = (w′1,m, w′2,m)∞m=1 ∈ H satisfies w′1,m = w1,m, w ′ 2,m = w2,m (m ≤ M) and w′1,m = 0, w′2,m = 0 (m > M). Moreover, we define ι∗ : H → H(M) as the map that extracts first M components. By abuse of notation, we write fW (M) for W (M) ∈ H(M) to indicate fι(W (M)). Finally, let A(M) : H(M) → H(M) be a linear operator such that A(M)W (M) = ι∗(Aι(W (M))), which is just a truncation of A. Similarly, let T aMW
(M) for W (M) ∈ H(M) be the operator corresponding to T aW for W ∈ H, i.e., T aMW (M) = ι∗(T aι(W (M))).
B.1 AUXILIARY LEMMAS
First, we show some key propositions to show the main results. To do so, we utilize the result by Muzellec et al. (2020) and Suzuki (2020).
Assumption 3.
(i) There exists a constant cµ such that µm ≤ cµm−2.
(ii) There exist B,U > 0 such that the following two inequalities hold for some a ∈ (1/4, 1) almost surely:
‖∇L̂(W )‖H ≤ B (∀W ∈ H), ‖∇L̂(W )−∇L̂(W ′)‖H ≤ L‖W −W ′‖H−a (∀W,W ′ ∈ H).
(iii) For any dataDn, L̂ is three times differentiable. Let∇3L̂(W ) be the third-order derivative of L̂(W ). This can be identified with a third-order linear form and∇3L̂(W )·(h, k) denotes the Riesz representor of l ∈ H 7→ ∇3L̂(W ) · (h, k, l). There exists α′ ∈ [0, 1), Cα′ ∈ (0,∞) such that ∀W,h, k ∈ H, ‖∇3L̂(W ) · (h, k)‖H−α′ ≤ Cα′‖h‖H‖k‖H, ‖∇
3L̂(W ) · (h, k)‖H ≤ Cα′‖h‖Hα′‖k‖H (a.s.).
Remark 1. In the analysis of Bréhier & Kopec (2016); Muzellec et al. (2020); Suzuki (2020), Assumption 3-(iii) is imposed for any finite dimensional projectionL(W (M)) as a function onH(M)) for all M ≥ 1 instead of L(W ) as a function of H. However, the condition on L(W ) gives a sufficient condition for any finite dimensional projection in our setting. Thus, we employed the current version.
Assumption 4. For the loss function `(y, f(x)) = (y − f(x))2, the following conditions holds:
(i) There exists C > 0 such that for any fW (W ∈ H), it holds that
EX,Y [(`(Y, fW (X))− `(Y, f∗(X)))2] ≤ C(L(fW )− L(f∗)).
(ii) β > 0 is chosen so that, for any h : Rd → R and x ∈ supp(PX), it holds that EY |X=x [ exp ( − βn (`(Y, h(x))− `(Y, f ∗(x))) )] ≤ 1.
(iii) There exists Lh > 0 such that ‖∇W `(Y, hW (X)) − ∇W `(Y, hW ′(X))‖H ≤ Lh‖W − W ′‖H (∀W,W ′ ∈ H) almost surely.
(iv) There exists Ch such that ‖hW − hW ′‖∞ ≤ Ch‖W −W ′‖H (W,W ′ ∈ H). Proposition 3. Assume Assumption 3 holds and β > η. Suppose that ∃R̄ > 0, 0 ≤ `(Y, fW (X)) ≤ R̄ for any W ∈ H (a.s.). Let ρ = 11+λη/µ1 and b = µ1 λ B + cµ βλ . Accordingly, let b̄ = max{b, 1}, κ = b̄+ 1 and V̄ = 4b̄/( √ (1+ρ1/η)/2−ρ1/η). Then, the spectral gap of the dynamics is given by
Λ∗η = min
( λ
2µ1 , 12 ) 4 log(κ(V̄ + 1)/(1− δ)) δ (10)
where 0 < δ < 1 is a real number satisfying δ = Ω(exp(−Θ(poly(λ−1)β))). We define Λ∗0 = limη→0 Λ ∗ η (i.e., V̄ is replaced by 4b̄/( √ (1+exp(− λµ1 ))/2−exp(− λ µ1 ))). We also define CW0 = κ[V̄ + 1] + √
2(R̄+b)√ δ . Then, for any 0 < a < 1/4, the following convergence bound holds for almost sure
observation Dn: for either L = L or L = L̂,
|EWk [L(Wk)|Dn]− EW∼π∞ [L(W )|Dn]| (11) ≤ C1 [ CW0 exp(−Λ∗ηηk) + √ β
Λ∗0 η1/2−a
] = Ξ′k, (12)
where C1 is a constant depending only on cµ, B, L,Cα′ , a, R̄ (independent of η, k, β, λ).
Proposition 4. Assume that Assumptions 3 and 4 hold. Let α̃ := 1/{2(α + 1)} for a given α > 0 and θ be an arbitrary real number satisfying 0 < θ < 1 − α̃. Assume that the true function fo can be represented by hW∗ = fo for W ∗ ∈ Hθ(α+1). Then, if M ≥ min { λα̃/2[θ(α+1)]β1/2[θ(α+1)], λ−1/2(α+1), n1/2[θ(α+1)] } , the expected excess risk is bounded by
EDn [ E W
(M) k [L(h T α/2 M W (M) k
)|Dn]− L(fo) ]
≤ C max { (λβ) 2α̃/θ 1+α̃/θ n− 1 1+α̃/θ , λ−α̃β−1, λθ, 1/n }
+ Ξ′k, (13)
where C is a constant independent of n, β, λ, η, k.
Proof. Repeating the same argument in Proposition 1 and using the same notation, Proposition 3 gives
|E W
(M) k [L(W (M)k )|Dn]− EW∼π(M)∞ [L(W )|Dn]| ≤ Ξ ′ k,
for any 1 ≤ M ≤ ∞. Therefore, we just need to bound the following quantity:∣∣∣EDn [EW (M)∼π(M)∞ [L(hTα/2M W (M))|Dn]]− L(fo)∣∣∣. We define ‖W (M)‖H(M) := ‖ι∗(W (M))‖H for W (M) ∈ H(M). For a > 0, we define H (M) a be the projection of Ha to the first M components, H(M)a = {ι(W ) | W ∈ Ha}, and we define ‖W (M)‖H(M)a := ‖ι ∗(W (M))‖Ha (note that since H (M) a is a finite dimensional linear space, it is same as H as a set). Let ν(M)β be the Gaussian measure on H(M) with mean 0 and covariance (βA(M))−1, and ν̃(M)β be the Gaussian measure corresponding to the random variable T α/2 M W (M) with W (M) ∼ ν(M)β . Let the concentration function be
φ (M) β,λ ( ) := inf
W∈H(M)α+1: L(hW )−L(fo)≤ 2
βλ‖W‖2 H(M)α+1 − log ν̃(M)β ({W ∈ H (M) : ‖W‖H(M) ≤ }) + log(2),
where, if there does not exist W ∈ H(M)α+1 that satisfies the condition inf , then we define φ (M) β,λ ( ) = ∞, then Let ∗ > 0 be
∗ := max{inf{ > 0 | φβ,λ( ) ≤ β 2}, 1/n}.
Then, Suzuki (2020) showed the following bound:∣∣∣∣EDn [EW (M)∼π(M)∞ [L(hTα/2(M)W (M))|Dn]− L(fo) ]∣∣∣∣
≤ C max { ∗2, ( β n ∗2 + n− 1 1+α̃/θ (λβ) 2α̃/θ 1+α̃/θ ) , 1
n
} . (14)
They also showed that, for M =∞, it holds that ∗2 . max { (λβ)−α̃β−(1−α̃), λθ, n−1 } = max { λ−α̃β−1, λθ, n−1 } .
Substituting this bound of ∗ to Eq. (14), we obtain Eq. (13) for M =∞. Moreover, in their proof, if M ≥ ( ∗)−1/[θ(α+1)], then
inf W∈H(M)α+1:
L(hW )−L(fo)≤ 2
βλ‖W‖2 H(M)α+1 . β( ∗)2.
Finally, since ν̃(M)β is a marginal distribution of ν̃ (∞) β , it holds that
− log ν̃(M)β ({W ∈ H (M) : ‖W‖H(M) ≤ }) ≤ − log ν̃ (∞) β ({W ∈ H : ‖W‖H ≤ }).
Therefore, as long as M ≥ ( ∗)−1/[θ(α+1)], the rate of ∗ is not deteriorated from M = ∞. In other words, if M ≥ min { λα̃/2[θ(α+1)]β1/2[θ(α+1)], λ−θ/2[θ(α+1)], n1/2[θ(α+1)] } , the bound (13) holds.
Remark 2. Suzuki (2020) showed Proposition 4 under a condition α > 1/2. However, this is used only to ensure Assumption 3. In our setting, we can show Assumption 3 holds directly and thus we may omit the condition α > 1/2.
B.2 PROOFS OF PROPOSITION 1, THEOREM 2 AND COROLLARY 1
Here, we give the proofs of Proposition 1 and Theorem 2 simultaneously. Proof of Proposition 1 and Theorem 2. Let R̄ = (2 ∑∞ m=1 amR + U)
2. Then, we can easily check that (yi − fW (xi))2 ≤ R̄. As stated above, we use Propositions 3 and 4 to show the statements.
First, we show Proposition 1 for the dynamics of W (M)k for any 1 ≤ M ≤ ∞. However, it suffices to show the statement only for M = ∞ because the finite dimensional version can be seen as a
specific case of the infinite dimensional one. Actually, the dynamics of W (M)k is same as that of ι(W̃k) where W̃k ∈ H obeys the following dynamics:
W̃k+1 = Sη ( W̃k − η∇L̂(fι(W̃k)) + √ 2η
β ξk
) .
This is because fι(W̃k) is determined by only the first M components ι(W̃k), ι(∇L̂(fι(W̃k))) = ∇W (M)L̂(fW (M))|W (M)=ι(W̃k) and Sη is a diagonal operator. Since the components of W̃k with indexes higher than M does not affect the objective, smoothness of the objective is not lost. The stationary distribution π(M)∞ of the continuous dynamics corresponding to W (M) is a probability measure onH(M) that satisfies
dπ (M) ∞
dν (M) β
(W (M)) ∝ exp(−βL̂(fW (M))),
where ν(M)β is the Gaussian measure on RM×(d+2) with mean 0 and covariance (βA(M))−1. We can notice that this is the marginal distribution of the stationary distribution of the continuous time counterpart of W̃k: dπ̃∞(W̃ ) ∝ exp(−βL̂(fι(W̃ )))dνβ . Therefore, we just need to consider an infinite dimensional one. For this reasoning, we show the convergence for the original infinite dimensional dynamics (Wk)∞k=1. The convergence of the finite dimensional one (W (M) k ) ∞ k=1 can be shown by the same manner using the argument above.
To show Proposition 1, we use Propositions 3. To do so, we need to check validity of Assumptions 3. First, we check Assumption 3. Assumption 3-(i) is ensured by Assumption 1. Next, we check Assumption 3-(ii). The boundedness of the gradient can be shown as follows:
‖∇L̂(fW )‖2H
= ∞∑ m=1 (∥∥∥ 1 n n∑ i=1 2(fW (xi)− yi)w̄2,mam[xi; 1]σ′m(w>1,m[xi; 1]) ∥∥∥2
+ ∣∣∣ 1 n n∑ i=1 2(fW (xi)− yi)am tanh′(w2,m/R)σm(w>1,m[xi; 1]))∞m=1 ∣∣∣2)
≤ ∞∑ m=1 4R̄R2a2m(d+ 1)C 2 σ + 4R̄a 2 m
(∵ |fW (xi)− yi| ≤ R̄, ‖σ′m‖∞ ≤ Cσ, ‖ tanh ′ ‖∞ ≤ 1)
≤4R̄[R2C2σ(d+ 1) + 1] ∞∑ m=1 a2m <∞.
Similarly, we can show the Lipschitz continuity of the gradient as
‖∇L̂(fW )−∇L̂(fW ′)‖2H
≤ ∞∑ m=1 µ−2α1m µ 2α1 m { 4R̄a2m(d+ 1)C 2 σ[(w2,m − w′2,m)2 +R2‖w1,m − w′1,m‖2]
+ 4R̄a2m[(w2,m − w′2,m)2/R2 + C2σ(d+ 1)‖w1,m − w′1,m‖2] } (∵ ‖ tanh′′ ‖∞ ≤ 1)
≤ 4R̄[(d+ 1)C2σ(1 +R2) + 1/R2 + C2σ(d+ 1)] max m∈N {µ−2α1m a2m}
× ∞∑ m=1 µ2α1m [(w2,m − w′2,m)2 + ‖w1,m − w′1,m‖2]
. ‖W −W ′‖2H−α1 .
We can also verify Assumption 3-(iii) in a similar way. Then, we have verified Assumption 3. Therefore, we may apply Proposition 3, and then we obtain Proposition 1.
Next, we show Theorem 2 by using Proposition 4. For that purpose, we need to we verify Assumption 4. The first condition can be verified as
EX,Y [((Y − fW (X))2 − (Y − fo(X))2)2] = EX, [((f
o(X) + − fW (X))2 − 2)2] = EX [((f
o(X)− fW (X))2 + 2 (fo(X)− fW (X)))2] = EX [(f
o(X)− fW (X))4 + 2 (fo(X)− fW (X))(fo(X)− fW (X))2 + 2(fo(X)− fW (X))2] = ‖fo − fW ‖2∞EX [(fo(X)− fW (X))2] + U2EX [(fo(X)− fW (X))2] ≤ R̄EX [(fo(X)− fW (X))2] = R̄(L(fW )− L(fo)).
The second condition can be checked as follows. Note that
EY |X=x
( exp { −β n [(Y − fW (x))2 − (Y − fo(x))2] })
= E ( exp [ −β n (fo(x)− fW (x))2 − 2 (fW (x)− fo(x))] })
= exp [ −β n (fo(x)− fW (x))2 ] E { exp [ 2β n (fW (x)− fo(x)) ]} ≤ exp [ −β n (fo(x)− fW (x))2 ] exp [ 1 8 4β2 n2 4U2(fW (x)− fo(x))2 ] .
Thus, under the condition β ≤ n/(2U2), the right hand side can be upper bounded by
exp [ −β n ( 1− 2U 2β n ) (fW (x)− fo(x))2 ] ≤ 1.
Next, we check the third and fourth conditions. Noting that
∇WhW (X) = ( am(µ −α/2 m w2,m)µ −α/2 m [xi; 1]σ ′ m(µ −α/2 m w > 1,m[xi; 1]),
amµ −α/2 m tanh ′(µ−α/2m w2,m/R)σm(µ −α/2 m w > 1,m[xi; 1])) ∞ m=1 )∞ m=1 ,
we have that
‖∇WhW (X)‖2H
≤ ∞∑ m=1 a2mµ −α m [(d+ 1)R 2C2σ + 1]
≤ [(d+ 1)R2C2σ + 1] ∞∑ m=1 µ−α+2α1m
≤ [(d+ 1)R2C2σ + 1]c−α+2α1µ ∞∑ m=1 m−2(−α+2α1) =: C1 <∞
(∵ −α+ 2α1 = α1 > 1/2), and
‖∇WhW (X)−∇WhW ′(X)‖2H
≤ ∞∑ m=1 a2mµ −α m (d+ 1)[µ −α m (w2,m − w′2,m)2 +R2µ−αm ‖w1,m − w′1,m‖2]
+ a2mµ −α m [µ −α m (w2,m − w′2,m)2/R2 + C2σ(d+ 1)µ−αm ‖w1,m − w′1,m‖2]
≤ ∞∑ m=1 a2mµ −2α m [(d+ 1)(1 +R 2) + 1/R2 + C2σ(d+ 1)][‖w1,m − w′1,m‖2 + (w2,m − w′2,m)2]
≤ c2α1µ max m {µ2(α1−α)m }[(d+ 1)(1 +R2) + 1/R2 + C2σ(d+ 1)]‖W −W ′‖2H =: C2‖W −W ′‖2H,
for a constant 0 < C2 <∞. Therefore, it holds that
|hW (X)− hW ′(X)|2 ≤ C1‖W −W ′‖2H,
which yields the forth condition, and we also have
‖∇W `(Y, hW (X))−∇W `(Y, hW ′(X))‖2H =‖2(hW (X)− Y )∇WhW (X)− 2(hW ′(X)− Y )∇WhW ′(X)‖2H ≤2‖2(hW (X)− Y )(∇WhW (X)−∇WhW (X))‖2H
+ 2‖2(hW (X)− hW ′(X))∇WhW ′(X)‖2H ≤ 8R̄C2‖W −W ′‖2H + 8C21‖W −W ′‖2H . ‖W −W ′‖2H,
which yields the third condition.
Since fo ∈ Fγ , there exists W ∗ ∈ Hγ such that fo = fW∗ . Therefore, applying Proposition 4 with α = α1 (α̃ = 1/[2(α1 + 1)]) and θ = γ/(1 + α1) (since γ < 1/2 + α1, the condition θ < 1 − α̃ is satisfied), we obtain that for M ≥ min { λ1/4γ(α1+1)β1/2γ , λ−1/2(α1+1), n1/2γ } , the following excess risk bound holds:
EDn [ E W
(M) k
[L(W (M)k )|Dn]− L(f ∗) ] . max { (λβ) 2α̃/θ 1+α̃/θ n− 1 1+α̃/θ , λ−α̃β−1, λθ, 1/n } + Ξk.
Finally, by noting L(W (M)k )− L(f∗) = ‖fW (M)k − f ∗‖2L2(PX), we obtain the assertion.
Finally, we give the proof of Corollary 1.
Proof of Corollary 1. Note that
fW (x)
= ∞∑ m=1 amw̄2,mσm(w > 1,m[x; 1])
= ∞∑ m=1 µα1m w̄2,mµ qα2 m µ −qα2 m µ sα2 m σ(µ −α2 m w > 1,m[x; 1]) (∵ am = µ α1 m , bm = µ α2 m )
= ∞∑ m=1 µα1+qα2m w̄2,mµ −(s−q)α2 m σ(µ −α2 m w > 1,m[x; 1]).
Therefore, we may redefine α′1 ← α1 +qα2 and s′ ← s−q so that we obtain another representation of the model Fγ :
Fγ = { fW (x) =
∞∑ m=1 µ α′1 m w̄2,mσ̌m(w > 1,m[x; 1]) ∣∣∣W ∈ Hγ , ‖W‖Hγ ≤ 1 } ,
where σ̌m(·) = µ−s ′α2 m σ(µ −α2 m ·). Note that the condition 0 ≤ q ≤ s−3 gives s− q ≥ 3. Therefore, Assumptions 3 and 4 are valid even for the redefined parameters α′1, s ′ and σ̌m instead of α1, s and σm. Therefore, we can apply Theorem 2 by simply replacing α1 by α′1 = α1 + qα2.

Review summaries:
The paper shows that a two-layer neural network (although an extension to deeper models seem unproblematic) may outperform a class of linear functions in terms of the excess risk learning rate, and in a minimax optimality analysis, and when approximating a target function from the neural network class. The paper essentially shows that linear functions have a problem with the non-convexity of the neural network class, and approximate the slow rate of 1/(n)^(1/2) for increasing dimension. A neural network trained with noisy stochastic gradient descent on the other hand has a faster rate, depending on several parameters.