Towards Understanding the Spectral Bias of Deep Learning

1 Introduction

Over-parameterized neural networks have achieved great success in many applications. However, the success of deep learning has not been well understood. In order to understand neural network training, a recent work [Rahaman et al., 2019] pointed out an intriguing phenomenon called spectral bias, which says that during training, neural networks tend to learn the components of lower complexity faster. Similar observation has also been pointed out in [Xu et al., 2019; Xu et al., 2020]. The concept of spectral bias is appealing because this may intuitively explain why over-parameterized neural networks can achieve a good generalization performance without overfitting. During training, the networks fit the low complexity components first and thus lie in the concept class of low complexity. Arguments like this may lead to rigorous guarantee for generalization.

�Equal contribution

Great efforts have been made in search of explanations about the spectral bias. [Rahaman et al., 2019] evaluated the Fourier spectrum of ReLU networks and empirically showed that the lower frequencies are learned first; also lower frequencies are more robust to random perturbation. [Andoni et al., 2014] showed that for a sufficiently wide two-layer network, gradient descent with respect to the second layer can learn any low degree bounded polynomial. [Xu, 2018] provided Fourier analysis to two-layer networks and showed similar empirical results on one-dimensional functions and real data. [Nakkiran et al., 2019] used information theoretical approach to show that networks obtained by stochastic gradient descent can be explained by a linear classifier during early training. These studies provide certain explanations about why neural networks exhibit spectral bias in real tasks. But explanations in the theoretical aspect, if any, are to some extent limited. For example, Fourier analysis is usually done in the one-dimensional setting, and thus lacks generality.

Meanwhile, a recent line of work has taken a new approach to analyze neural networks based on the neural tangent kernel (NTK) [Jacot et al., 2018]. In particular, they show that under certain over-parameterization condition, the neural network trained by gradient descent behaves similarly to the kernel regression predictor using the neural tangent kernel. For training a neural network with hidden layer width m and sample size n, recent optimization results on the training loss in the so-called "neural tangent kernel regime" can be roughly categorized into the following two families: (i) Without any assumption on the target function (the function used to generate the true labels based on the data input), if the network width m ¥ polypn, λ�1 min q, where λmin is the smallest eigenvalue of the NTK Gram matrix, then square loss/crossentropy loss can be optimized to zero [Du et al., 2019b; Allen-Zhu et al., 2019; Du et al., 2019a; Zou et al., 2019; Zou and Gu, 2019; Lee et al., 2019]; and (ii) If the target function has bounded norm in the NTK-induced reproducing kernel Hilbert space (RKHS), then global convergence can be achieved with milder requirements on m [Arora et al., 2019; Cao and Gu, 2019].

Inspired by these works mentioned above in the neural tangent kernel regime, in this paper we study the spectral bias of over-parameterized two-layer neural networks and its connection to the neural tangent kernel. We show that, given a training data set that is generated based on a target function, a fairly narrow network, although cannot fit the training data well due to its limited width, can still learn certain low-complexity components of the target function in the eigenspace corresponding to large eigenvalues of neural tangent kernel. As the width of the network increases, more high-frequency components of the target function can be learned with a slower convergence rate. As a special case, our result implies that when the input data follows uniform distribution on the unit sphere, polynomials of lower degrees can be learned by a narrower neural network at a faster rate. We also conduct experiments to corroborate the theory we establish.

Our contributions are as follows:

1. We prove a generic theorem for arbitrary data distributions, which states that under certain sample complexity and over-parameterization conditions, the convergence of the training error along different eigendirections of NTK relies on the corresponding eigenvalues. This theorem gives a more precise control on the regression residual than [Su and Yang, 2019], where the authors focused on the case when the labeling function is close to the subspace spanned by the first few eigenfunctions.
2. We establish a rigorous explanation for spectral bias based on the aforementioned theoretical results without any specific assumptions on the target function. We show that the error terms from different frequencies are controlled by the eigenvalues of the NTK, and the lower-frequency components can be learned with less training examples and narrower networks at a faster convergence rate.
3. For uniformly distributed input data, we further give convergence results of two-layer ReLU networks by characterizing the spectra of the neural tangent kernel. We show that in this setting the eigenvalues of neural tangent kernel are µk � Ωpmaxtk�d�1 , d�k�1 uq, k ¥ 0, with corresponding eigenfunctions being the k-th order spherical harmonics. Our result is better than the bound Ωpk�d�1 q derived in [Bietti and Mairal, 2019] when d " k, which is in a more practical setting.

1.1 Additional RelatedWork

A few theoretical results have been established towards understanding the spectra of neural tangent kernels. To name a few, [Bach, 2017] studied two-layer ReLU networks by relating it to kernel methods, and proposed a harmonic decomposition for the functions in the reproducing kernel Hilbert space which we utilize in our proof. Based on the technique in [Bach, 2017], [Bietti and Mairal, 2019] studied the eigenvalue decay of integrating operator defined by the neural tangent kernel on unit sphere by using spherical harmonics. [Vempala and Wilmes, 2019] calculated the eigenvalues of neural tangent kernel corresponding to two-layer neural networks with sigmoid activation function. [Basri et al., 2019] established similar results as [Bietti and Mairal, 2019], but considered the case of training the first layer parameters of a two-layer networks with bias terms. [Yang and Salman, 2019] studied the the eigenvalues of integral operator with respect to the NTK on Boolean cube by Fourier analysis. Very recently, [Bordelon et al., 2020] gave a spectral analysis on the generalization error of NTK-based kernel ridge regression. [Basri et al., 2020] studied the convergence of full training residual with a focus on one-dimensional, non-uniformly distributed data. [Advani et al., 2020] gave an average case analysis of the learning dynamics of large neural networks trained by gradient descent and studied the generalization of over-parameterized neural networks.

2 Preliminaries

In this section we introduce the basic problem setup including the neural network structure and the training algorithm, as well as some background on the neural tangent kernel proposed recently in [Jacot et al., 2018].

Notation. We use lower case, lower case bold face, and upper case bold face letters to denote scalars, vectors and matrices respectively. For a vector v � pv1, . . . , vdqT P Rd and a number 1 / p � 8, we denote its p�norm by }v}p � p�di�1 |vi|p q1{p . We also define infinity norm by }v}8 � maxi |vi|. For a matrix A � pAi,jqm�n, we use }A}0 to denote the number of non-zero entries of A, and use }A}F � p�di,j�1 A2i,j q1{2 to denote its Frobenius norm. Let }A}p � max}v}p/1 }Av}p for p ¥ 1, and }A}max � maxi,j |Ai,j|. For two matrices A,B P Rm�n , we define xA, By � TrpAJBq. We use A ' B if A � B is positive semi-definite. In addition, we define the asymptotic notations Op�q, r Op�q, Ωp�q and rΩp�q as follows. Suppose that an and bn be two sequences. We write an � Opbnq if limsupn�8 |an{bn| � 8, and an � Ωpbnq if liminfn�8 |an{bn| ¡ 0. We use r Op�q and rΩp�q to hide the logarithmic factors in Op�q and Ωp�q.

Problem Setup. Here we introduce the basic problem setup. We consider two-layer fully connected neural networks of the form

whereW1 P Rm�pd�1q , W2 P R1�m are 1 the first and second layer weight matrices respectively, and σp�q � maxt0, �u is the entry-wise ReLU activation function. The network is trained according to the square loss on n training examples S � tpxi, yiq : i P rnsu:

where θ is a small coefficient to control the effect of initialization, and the data inputs txiuni�1 are assumed to follow some unknown distribution τ on the unit sphere Sd � Rd�1 . Without loss of generality, we also assume that |yi| / 1.

1 Here the input dimension is d�1 since in this paper we assume that all training data lie in the d-dimensional unit sphere Sd � Rd�1 .

Algorithm 1 GD for DNNs starting at Gaussian initialization

Input: Number of iterations T , step size η.

Generate each entry of Wp0q 1 and Wp0q 2 from Np0, 2{mq and Np0, 1{mq respectively.

end for

Output: WpT q .

We first randomly initialize the parameters of the network, and run gradient descent for both layers. We present our detailed neural network training algorithm in Algorithm 1. The initialization scheme for Wp0q given in Algorithm 1 is known as He initialization [He et al., 2015]. It is consistent with the initialization scheme used in [Cao and Gu, 2019].

Neural Tangent Kernel. Many attempts have beenmade to study the convergence of gradient descent assuming the width of the network is extremely large [Du et al., 2019b; Li and Liang, 2018]. When the width of the network goes to infinity, with certain initialization on the model weights, the limit of inner product of network gradients defines a kernel function, namely the neural tangent kernel [Jacot et al., 2018]. The neural tangent kernel is derived by linearizing the network around its random initialization, and is defined as

For two-layer networks, by the definition of Gaussian initialization, we can obtain the following calculations of the neural tangent kernel by the law of large numbers:

where

Since we apply gradient descent to both layers, the neural tangent kernel is the sum of the two different kernel functions and clearly it can be reduced to one layer training setting. These two kernels are arc-cosine kernels of degree 0 and 1 [Cho and Saul, 2009], which are given as κ1px, x1 q � pκ1pxx, x1 yp}x}2 }x1 }2 qq, κ2px, x1 q � pκ2pxx, x1 y{p}x}2 }x1 }2 qq, where

Integral Operator. The theory of integral operator with respect to kernel function has been well studied in literature [Smale and Zhou, 2007; Rosasco et al., 2010] thus we only give a brief introduction here. Let L2τ pXq be the Hilbert space of square-integrable functions with respect to a Borel measure τ from X � R. For any continuous kernel function κ : X �X � R and τ we can define an integral operator Lκ on L2τ pXq by

»

It has been pointed out in [Cho and Saul, 2009] that arccosine kernels are positive semi-definite. Thus the kernel function κ defined by (2.1) is positive semi-definite being a product and a sum of positive semi-definite kernels. Clearly this kernel is also continuous and symmetric, which implies that the neural tangent kernel κ is a Mercer kernel.

3 Main Results

In this section we present our main results. We first give a general result on the convergence rate of gradient descent along different eigendirections of neural tangent kernel, and then apply this result to give an explicit convergence rate for uniformly distributed data on the unit sphere.

3.1 Convergence Analysis of Gradient Descent

In this section we study the convergence of Algorithm 1. Instead of studying the standard convergence of loss function value, we provide a refined analysis on the speed of convergence along different directions defined by the eigenfunctions of Lκ. For simplicity, we focus on the setting where the data inputs have normalized lengths, i.e., x P Sd .

Let tλiui¥1 with λ1 ¥ λ2 ¥ � � � be the strictly positive eigenvalues of Lκ, and φ1p�q, φ2p�q, . . . be the corresponding orthonormal eigenfunctions (i.e., Lκpφiqpxq � λiφipxq for all x P Sd ). Set vi � n�1{2 pφipx1q, . . . , φipxnqqJ, i � 1, 2, . . .. Note that Lκ may have eigenvalues with multiplicities larger than 1 and λi, i ¥ 1 are not distinct. Therefore for any integer k, we define rk as the sum of the multiplicities of the first k distinct eigenvalues of Lκ. Define Vrk � pv1, . . . , vrk q. By definition, vi, i P rrks are rescaled restrictions of orthonormal functions in L2τ pSd q on the training examples. Therefore we can expect them to form a set of almost orthonormal bases in the vector space Rn . The following lemma follows by standard concentration inequality.

Lemma 3.1. Suppose that |φipxq| / M for all x P Sd and i P rrks. For any δ ¡ 0, with probability at least 1 � δ,

where C is an absolute constant.

Denote y � py1, . . . , ynqJ and pyptq � θ � pfWptqpx1q, . . . , fWptqpxnqqJ for t � 0, . . . , T. Then Lemma 3.1 shows that the convergence rate of }VJ rk py � pyptq q}2 roughly represents the speed gradient descent learns the components of the target function corresponding to the first rk eigenvalues. The following theorem gives the convergence guarantee of }VJ rk py � pyptq q}2.

Theorem 3.2. Suppose |φjpxq| / M for j P rrks and x P Sd . For any ϵ, δ ¡ 0 and integer k, if n ¥ rΩpϵ�2 �maxtpλrk � λrk�1q�2 , M 4 r2 k uq, m ¥ rΩppolypT, λ�1 rk , ϵ�1 qq, then with probability at least 1�δ, Algorithm 1 with η � r Opm�1 θ�2 q, θ � r Opϵq satisfies

Theorem 3.2 shows that the convergence rate of }VJ rk py � pyptq q}2 (or its normalized version, }VJ rk py � pyptq q}2{}VJ rk y}2) is determined by the rk-th eigenvalue λrk . This reveals the spectral bias of neural network training under the NTK regime. Specifically, when the network is wide enough and the sample size is large enough, gradient descent first learns the target function along the eigendirections of neural tangent kernel with larger eigenvalues, and then learns the rest components corresponding to smaller eigenvalues 2 . Moreover, by showing that learning the components corresponding to larger eigenvalues can be done with smaller sample size and narrower networks, our theory pushes the study of neural networks towards more practical settings. Therefore, Theorem 3.2 provides an explanation of the observations given in [Rahaman et al., 2019] in the NTK regime.

It is also worth noting that Theorem 3.2 is not based on the common NTK-type assumption that the target function belongs to the NTK-induced RKHS [Arora et al., 2019; Cao and Gu, 2019]. Instead, Theorem 3.2 works for arbitrary labeling, and is therefore rather general.

Comparison with existing results. The most relevant results to ours are [Arora et al., 2019; Su and Yang, 2019; Bordelon et al., 2020; Basri et al., 2020]. Compared with [Arora et al., 2019] which focuses on the setting where the network is wide enough to guarantee global convergence, our result works for narrower networks for which global convergence may not even be possible. Compared with [Su and Yang, 2019], the key difference is that while [Su and Yang, 2019] studied the full residual }y � pypTq }2 and required that the target function lies approximately in the eigenspace of large eigenvalues of the neural tangent kernel, our result in Theorem 3.2 works for arbitrary target function, and shows that even if the target function has very high frequency components, its components in the eigenspace of large eigenvalues can still be learned very efficiently by neural networks. More recently, [Bordelon et al., 2020] studied the solution of NTK-based kernel ridge regression. Compared with their result, our analysis is directly on the practical neural network training procedure, and specifies the width requirement for a network to successfully learn a certain component of the target function. Another recent work by [Basri et al., 2020] provides theoretical analysis for one-dimensional nonuniformly distributed data, and only studies the convergence of the full residual vector. In comparison, our results cover high-dimensional data, and provide a more detailed analysis on the convergence of different projections of the residual.

3.2 Spectral Analysis of NTK for UniformData

We now study the case when the data inputs are uniformly distributed over the unit sphere as an example where the eigendecompositon of NTK can be calculated. We present our results on the spectral analysis of neural tangent kernel in the form of a Mercer decomposition, which explicitly gives the eigenvalues and eigenfunctions of NTK.

2 Note that by the definition we have λrk � 1 for all k. More details are given in Section 3.2 and th proofs in the appendix

Theorem 3.3. For any x, x1 P Sd , the Mercer decomposition of the neural tangent kernel κ : Sd � Sd � R is

where Yk,j for j � 1, � � � ,Npd, kq are linearly independent spherical harmonics of degree k in d � 1 variables with Npd, kq � 2k�d�1 k �k�d�2 d�1 � and orders of µk are given by

where j P N�. Specifically, we have µk � Ω �k�d�1 � when k " d and µk � Ω �d�k�1 � when d " k, k � 2, 4, 6, . . ..

Remark 3.4. The µk's in Theorem 3.3 are the distinct eigenvalues of the integral operator Lκ on L2τd pSd q defined by

»

where τd is the uniform probability measure on unit sphere Sd . Therefore the eigenvalue λrk in Theorem 3.2 is just µk�1 given in Theorem 3.3 when τd is uniform distribution.

Remark 3.5. [Vempala and Wilmes, 2019] studied twolayer neural networks with sigmoid activation function, and proved that in order to achieve ϵ0 � ϵ error, it requires T � pd � 1qOpkq log p}f � }2{ϵq iterations and m � pd � 1qOpkqpolyp}f � }2{ϵq wide neural networks, where f � is the target function, and ϵ0 is certain function approximation error. Another highly related work is [Bietti and Mairal, 2019], which gives µk � Ωpk�d�1 q. The order of eigenvalues we present appears as µk � Ωpmaxpk�d�1 , d�k�1 qq. This is better when d " k, which is closer to the practical setting.

3.3 Convergence for Uniformly Distributed Data

In this subsection, we apply the result in the previous subsection and give explicit convergence rate for uniformly distributed data on the unit sphere. The following two corollaries are based on the calculation of the spectral decomposition of the NTK for uniform data in Theorem 3.1.

Corollary 3.6. Suppose that k " d, and the sample txiuni�1 follows the uniform distribution τd on the unit sphere Sd . For any ϵ, δ ¡ 0 and integer k, if n ¥ rΩpϵ�2 � maxtk2d�2 , k2d�2 r2 k uq, m ¥ rΩppolypT, kd�1 , ϵ�1 qq, then with probability at least 1 � δ, Algorithm 1 with η � r Oppmθ2 q�1 q, θ � r Opϵq satisfies

where rk � �k�1 k1�0 Npd, k1 q and Vrk � pn�1{2φjpxiqqn�rk with φ1, . . . , φrk being a set of orthonormal spherical harmonics of degrees up to k � 1.

Corollary 3.7. Suppose that d " k, and the sample txiuni�1 follows the uniform distribution τd on the unit sphere Sd . For any ϵ, δ ¡ 0 and integer k, if n ¥ rΩpϵ�2d2k�2 r2 k q, m ¥ rΩppolypT, dk�2 , ϵ�1 qq, then with probability at least 1 � δ, Algorithm 1 with η � r Oppmθ2 q�1 q, θ � r Opϵq satisfies

where rk � �k�1 k1�0 Npd, k1 q and Vrk � pn�1{2φjpxiqqn�rk with φ1, . . . , φrk being a set of orthonormal spherical harmonics of degrees up to k � 1.

Corollaries 3.6 and 3.7 further illustrate the spectral bias of neural networks by providing exact calculations of λrk , Vrk and M in Theorem 3.2. They show that if the input distribution is uniform over unit sphere, then spherical harmonics with lower degrees are learned first by wide neural networks.

4 Experiments

In this section we present experimental results to verify our theory. Across all tasks, we train a two-layer neural networks with 4096 hidden neurons and initialize it exactly as defined in the problem setup. The optimization method is vanilla gradient descent, and the training sample size is 1000.

We first verify our theoretical results by learning linear combinations of spherical harmonics with data uniformly distributed over unit sphere. Consider

where the Pkptq is the Gegenbauer polynomial, and ζk, k � 1, 2, 4 are fixed vectors that are generated uniformly over S9 . Note that by the addition formula �Npd,kq j�1 Yk,jpxqYk,jpyq � Npd, kqPkpxx, yyq, every normalized Gegenbauer polynomial is a spherical harmonic, so f � pxq is a linear combination of spherical harmonics of order 1,2 and 4. The higher odd-order Gegenbauer polynomials are omitted because the spectral analysis showed that µk � 0 for k � 3, 5, 7 . . . .

Following our theoretical analysis, we denote vk � n�1{2 pPkpx1q, . . . , Pkpxnqq, and define the (approximate) projection length of residual rptq onto vk at step t as pak � |vJ k rptq |, where rptq � pf � px1q�θfWptq px1q, . . . , f � pxnq� θfWptqpxnqq and fWptqpxq is the neural network function.

The results are shown in Figure 1. Following the settings in [Rahaman et al., 2019], we test target functions with different magnitudes in different components. We can see in (a),(b) that the residual at the lowest frequency (k � 1) converges to zero first and then the second lowest (k � 2). The highest frequency component is the last one to converge. The convergence curves in logarithmic scales are given in (d),(e). It is clear that the convergence of different projection length is close to linear convergence, which is wellaligned with our theoretical analysis. Note that we performed a moving average of range 20 on these curves to avoid the heavy oscillation especially at late stage. In (c)(f) we further consider more general target functions, which are given by f � pxq � �i cospaixζ, xyq and f � pxq � �i xζ, xypi for some fixed unit vector ζ. The results demonstrate that Theorem 3.2 still holds even when the target function has highfrequency components.

We also do experiments for non-uniformly distributed data. Note that the eigendecomposition of NTK for general multidimensional input distributions do not necessarily have good analytic forms. Therefore, here we treat the non-uniform distribution of the input data as model misspecification, and test whether residual projections onto spherical harmonics of different degrees can still exhibit different convergence rates. We consider three types of non-uniform distributions:

• Piece-wise uniform distribution (PUD). We divide the unit sphere into two semi-spheres along a randomly drawn direction ζ0. A data input is then generated as follows: with probability 1{4, draw the input uniformly over the first unit sphere; with probability 3{4, draw the input uniformly over the second unit sphere.
• Normalized anisotropic Gaussian (NAG). We generate the data by first generating non-zero mean, anisotropic Gaussian vectors, and then normalize them to unit length. The Gaussian mean vector is generated elementwise from Unifpr�1, 1sq; the covariance matrix is generated as Σ � AJA, where A P Rd�D(d � 10,D � 20) has i.i.d. standard Gaussian entries.
• Normalized Gaussian mixture (NGM). Here the inputs are drawn from a mixture of 3 anisotropic Gaussian distributions described in the second example.

Figure 2 shows that the residual components corresponding to lower order polynomials are still learned relatively faster, even for the non-uniform distributions described above, where spherical harmonics are no longer exactly the eigenfunctions of NTK. This suggests that our theoretical results can tolerate certain level of model misspecification, and therefore the spectral bias characterized in Theorem 3.2 and Corollaries 3.6, 3.7 holds in a variety of problem settings.

5 Conclusion and FutureWork

This paper gives theoretical justification for spectral bias through analyzing the convergence behavior of two-layer ReLU networks. We show that the convergence of gradient descent in different directions depends on the corresponding eigenvalues of the neural tangent kernel, and give the exact convergence rate on different directions. We also conduct experiments on synthetic data to support our theoretical result.

Our current analysis of spectral bias is mostly in the NTK regime. Extending our results to the general setting beyond NTK regime is an important future work direction. Studying spectral bias for the training of deep neural networks with regularization [Hu et al., 2021] is another future direction.

Acknowledgements

We would like to thank the anonymous reviewers for their helpful comments. YC and QG are partially supported by the National Science Foundation IIS-2008981 and Salesforce Deep Learning Research Award. DZ is supported partially by the Research Grants Council of Hong Kong [Project # CityU

(d) P1pxq �P2pxq �P4pxq, Log-scale (e) P1pxq�3P2pxq�5P4pxq, Log-scale

(d) PUD, P1pxq � 3P2pxq �5P4pxq (e) NAG, P1pxq �3P2pxq �5P4pxq (f) NGM, P1pxq � 3P2pxq �5P4pxq

Figure 2: Convergence curves of projection lengths for non-uniform data distributions.

11307319], Hong Kong Institute for Data Science, National Science Foundation of China [Project No. 12061160462], and an InnoHK project of the Hong Kong Innovation and Technology Commission.