1. Introduction. Nowadays, deep neural networks have become one of the leading models in machine learning, and many theoretical results have been established to understand the training and generalization of neural networks. Among them, two kernel matrices are prominent in deep learning theory: conjugate kernel (CK) [24,27,44,54,58,69,71,74,82] and neural tangent kernel (NTK) [4,28,40]. The CK matrix defined in (5), which has been exploited to study the generalization of random feature regression, is the Gram matrix of the output of the last hidden layer on the training dataset. The NTK matrix, defined in (7), is the Gram matrix of the Jacobian of the neural network with respect to training parameters, characterizing the performance of a wide neural network through gradient flows. Both are related to the kernel machine and help us explore the generalization and training process of the neural network.

We are interested in the behaviors of CK and NTK matrices at random initialization. A recent line of work has proved that these two random kernel matrices will converge to their expectations when the width of the network becomes infinitely wide [7,40]. Although CK and NTK are usually referred to as these expected kernels in literature, we will always call CK and NTK the empirical kernel matrices in this paper, with a slight abuse of terminology.

In this paper, we study the random CK and NTK matrices of a two-layer fully connected neural network with input data X ∈ R d 0 ×n , given by f : R d 0 ×n → R n such that (1) f (X)

where W ∈ R d 1 ×d 0 is the weight matrix for the first layer, a ∈ R d 1 are the second layer weights, and σ is a nonlinear activation function applied to the matrix W X elementwisely. We assume that all entries of a and W are independently identically distributed by the standard Gaussian N (0, 1). We will always view the input data X as a deterministic matrix (independent of the random weights in a and W ) with certain assumptions. In terms of random matrix theory, we study the difference between these two kernel matrices (CK and NTK) and their expectations with respect to random weights, showing both asymptotic and nonasymptotic behaviors of these differences as the width of the first hidden layer d 1 is growing faster than the number of samples n. As an extension of [29], we prove that when n/d 1 → 0, the centered CK and NTK with appropriate normalization have the limiting eigenvalue distribution given by a deformed semicircle law, determined by the training data spectrum and the nonlinear activation function. To prove this global law, we further set up a limiting law theorem for centered sample covariance matrices with dependent structures and a nonlinear version of the Hanson-Wright inequality. These two results are very general, which makes them potentially applicable to different scenarios beyond our neural network model. For the nonasymptotic analysis, we establish concentration inequalities between the random kernel matrices and their expectations. As a byproduct, we provide lower bounds of the smallest eigenvalues of CK and NTK, which are essential for the global convergence of gradient-based optimization methods when training a wide neural network [59,60,63]. Because of the nonasymptotic results for kernel matrices, we can also describe how close the performances of the random feature regression and the limiting kernel regression are with a general dataset, which allows us to compute the limiting training error and generalization error for the random feature regression via its corresponding kernel regression in the ultra-wide regime.

1.1. Nonlinear random matrix theory in neural networks. Recently, the limiting spectra of CK and NTK at random initialization have received increasing attention from a random matrix theory perspective. Most of the papers focus on the linear-width regime d 1 ∝ n, using both the moment method and Stieltjes transforms. Based on moment methods, [67] first computed the limiting law of the CK for two-layer neural networks with centered nonlinear activation functions, which is further described as a deformed Marchenko-Pastur law in [64]. This result has been extended to sub-Gaussian weights and input data with real analytic activation functions by [19], even for multiple layers with some special activation functions. Later, [2] generalized their results by adding a random bias vector in pre-activation and a more general input data matrix. Similar results for the two-layer model with a random bias vector and random input data were analyzed in [68] by cumulant expansion. In parallel, by Stieltjes transform, [52] investigated the CK of a one-hidden-layer network with general deterministic input data and Lipschitz activation functions via some deterministic equivalent. [49] further developed a deterministic equivalent for the Fourier feature map. With the help of the Gaussian equivalent technique and operator-valued free probability theory, the limiting spectrum of NTK with one hidden layer has been analyzed in [3]. Then the limiting spectra of CK and NTK of a multi-layer neural network with general deterministic input data have been fully characterized in [29], where the limiting spectrum of CK is given by the propagation of the Marchenko-Pastur map through the network, while the NTK is approximated by the linear combination of CK's of each hidden layer. [29] illustrated that the pairwise approximate orthogonality assumption on the input data is preserved in all hidden layers. Such a property is useful to approximate the expected CK and NTK. We refer to [32] as a summary of the recent development in nonlinear random matrix theory.

Most of the results in nonlinear random matrix theory focus on the case when d 1 is proportional to n as n → ∞. We build a random matrix result for both CK and NTK under the ultra-wide regime, where d 1 /n → ∞ and n → ∞. As an intrinsic interest of this regime, this exhibits the connection between wide (or overparameterized) neural networks and kernel learning induced by limiting kernels of CK and NTK. In this article, we will follow general assumptions on the input data and activation function in [29] and study the limiting spectra of the centered and normalized CK matrix

where Y := σ (W X). Similar results for the NTK can be obtained as well. To complete the proofs, we establish a nonlinear version of the Hanson-Wright inequality, which has previously appeared in [49,52]. This nonlinear version is a generalization of the original Hanson-Wright inequality [1,36,72], and may have various applications in statistics, machine learning, and other areas. In addition, we also derive a deformed semicircle law for normalized sample covariance matrices without independence in columns. This result is of independent interest in random matrix theory as well.

1.2. General sample covariance matrices. We observe that the random matrix Y ∈ R d 1 ×n defined above has independent and identically distributed rows. Hence, Y Y is a generalized sample covariance matrix. We first inspect a more general sample covariance matrix Y whose rows are independent copies of some random vector y ∈ R n . Assuming n and d 1 both go to infinity but n/d 1 → 0, we aim to study the limiting empirical eigenvalue distribution of centered Wishart matrices in the form of (2) with certain conditions on y. This regime is also related to the ultra-high-dimensional setting in statistics [70].

This regime has been studied for decades starting in [14], where Y has i.i.d. entries and E[Y Y ] = d 1 Id. In this setting, by the moment method, one can obtain the semicircle law.

This normalized model also arises in quantum theory with respect to random induced states (see [8,9,26]). The largest eigenvalue of such a normalized sample covariance matrix has been considered in [22]. Subsequently, [21,45,70,87] analyzed the fluctuations for the linear spectral statistics of this model and applied this result to hypothesis testing for the covariance matrix. A spiked model for sample covariance matrices in this regime was recently studied in [30]. This kind of semicircle law also appears in many other random matrix models. For instance, [42] showed this limiting law for normalized sample correlation matrices. Also, the semicircle law for centered sample covariance matrices has already been applied in machine learning: [31] controlled the generalization error of shallow neural networks with quadratic activation functions by the moments of this limiting semicircle law; [35] derived a semicircle law of the fluctuation matrix between stochastic batch Hessian and the deterministic empirical Hessian of deep neural networks.

For general sample covariance, [80] considered the form Y = BXAfoot_0/2 with deterministic A and B, where X consists of i.i.d. entries with mean zero and variance one. The same result has been proved in [16] by generalized Stein's method. Unlike previous results, [85] tackled the general case, only assuming Y has independent rows with some deterministic covariance n . Though this is similar to our model in Section 4, we will consider more general assumptions on each row of Y , which can be directly verified in our neural network models.

1.3. Infinite-width kernels and the smallest eigenvalues of empirical kernels. Besides the above asymptotic spectral fluctuation of (2), we provide nonasymptotic concentrations of (2) in spectral norm and a corresponding result for the NTK. In the infinite-width networks, where d 1 → ∞ and n are fixed, both CK and NTK will converge to their expected kernels. This has been investigated in [27,44,54,74] for the CK and [4,7,28,40,47] for the NTK. Such kernels are also called infinite-width kernels in literature. In this current work, we present the precise probability bounds for concentrations of CK and NTK around their infinite-width kernels, where the difference is of order √ n/d 1 . Our results permit more general activation functions and input data X only with pairwise approximate orthogonality, albeit similar concentrations have been applied in [3,10,39,57,76].

A corollary of our concentration is the explicit lower bounds of the smallest eigenvalues of the CK and the NTK. Such extreme eigenvalues of the NTK have been utilized to prove the global convergence of gradient descent algorithms of wide neural networks since the NTK governs the gradient flow in the training process, see, for example, [6,23,28,59,60,63,76,83]. The smallest eigenvalue of NTK is also crucial for proving generalization bounds and memorization capacity in [6,57]. Analogous to Theorem 3.1 in [57], our lower bounds are given by the Hermite coefficients of the activation function σ . Besides, the lower bound of NTK for multi-layer ReLU networks is analyzed in [61].

1.4. Random feature regression and limiting kernel regression. Another byproduct of our concentration results is to measure the difference of performance between random feature regression with respect to 1

Y and corresponding kernel regression when d 1 /n → ∞. Random feature regression can be viewed as the linear regression of the last hidden layer, and its performance has been studied in, for instance, [33,38,49,50,52,53,55,56,67] under the linear-width regime. 1 In this regime, the CK matrix 1 d 1 Y Y is not concentrated around its expectation

under the spectral norm, where w is the standard normal random vector in R d 0 . But the limiting spectrum of CK is exploited to characterize the asymptotic performance and double descent phenomenon of random feature regression when n, d 0 , d 1 → ∞ proportionally. Several works have also utilized this regime to depict the performance of the ultra-wide random network by letting d 1 /n → ψ ∈ (0, ∞) first, getting the asymptotic performance and then taking ψ → ∞ (see [56,86]). However, there is still a difference between this sequential limit and the ultra-wide regime. Before these results, random feature regression has already attracted significant attention in that it is a random approximation of the Reproducing Kernel Hilbert Space (RKHS) defined by population kernel function K : R

when width d 1 is sufficiently large [11,12,71,73]. We point out that Theorem 9 of [10] has the same order √ n/d 1 of the approximation as ours, despite only for random Fourier features. In our work, the concentration between empirical kernel induced by 1 d 1 Y Y and the population kernel matrix K defined in (4) for X leads to the control of the differences of training/test errors between random feature regression and kernel regression, which were previously concerned by [10,41,55,57] in different cases. Specifically, [41] obtained the same kind of estimation but considered random features sampled from Gaussian processes. Our results explicitly show how large width d 1 should be so that the random feature regression gets the same asymptotic performance as kernel regression [55]. With these estimations, we can take the limiting test error of the kernel regression to predict the limiting test error of random feature regression as n/d 1 → 0 and d 0 , n → ∞. We refer [46,47,51,55], [18], Section 4.3, and references therein for more details in high-dimensional kernel ridge/ridgeless regressions. We emphasize that the optimal prediction error of random feature regression in linear-width regime is actually achieved in the ultra-wide regime, which boils down to the limiting kernel regression, see [53,55,56,86]. This is one of the motivations for studying the ultra-wide regime and the limiting kernel ridge regression.

In the end, we would like to mention the idea of spectral-norm approximation for the expected kernel , which helps us describe the asymptotic behavior of limiting kernel regression. For specific activation σ , kernel has an explicit formula, see [48,49,52], whereas generally, it can be expanded in terms of the Hermite expansion of σ [29,56,67]. Thanks to pairwise approximate orthogonality introduced in [29], Definition 3.1, we can approximate in the spectral norm for general deterministic data X. This pairwise approximate orthogonality defines how orthogonal is within different input vectors of X. With certain i.i.d. assumption on X, [47] and [18], Section 4.3, where the scaling d 0 ∝ n α , for α ∈ (0, 1], determined which degree of the polynomial kernel is sufficient to approximate . Instead, our theory leverages the approximate orthogonality among general datasets X to obtain a similar approximation. Our analysis presumably indicates that the weaker orthogonality X has, the higher degree of the polynomial kernel we need to approximate the kernel .

Notation. We use tr(A) = 1 n i A ii as the normalized trace of a matrix A ∈ R n×n and Tr(A) = i A ii . Denote vectors by lowercase boldface. A is the spectral norm for matrix A, A F denotes the Frobenius norm, and x is the 2 -norm of any vector x. A B is the Hadamard product of two matrices, that is,

and Var w [•] be the expectation and variance only with respect to random vector w. Given any vector v, diag(v) is a diagonal matrix where the main diagonal elements are given by v. λ min (A) is the smallest eigenvalue of any Hermitian matrix A.

Before stating our main results, we describe our model with assumptions. We first consider the output of the first hidden layer and empirical conjugate kernel (CK):

(5)

Observe that the rows of matrix Y are independent and identically distributed since only W is random and X is deterministic. Let the ith row of Y be y i , for 1 ≤ i ≤ d 1 . Then, we obtain a sample covariance matrix,

which is the sum of d 1 independent rank-one random matrices in R n×n . Let the second moment of any row y i be (3). Later on, we will approximate based on the assumptions of input data X.

Next, we define the empirical neural tangent kernel (NTK) for (1), denoted by H ∈ R n×n . From Section 3.3 in [29], the (i, j )th entry of H can be explicitly written as

where w r is the rth row of weight matrix W , x i is the ith column of matrix X, and a r is rth entry of the output layer a. In the matrix form, H can be written by

where the αth column of S is given by diag

We introduce the following assumptions for the random weights, nonlinear activation function σ , and input data. These assumptions are basically carried on from [29].

ASSUMPTION 1.1. The entries of W and a are i.i.d. and distributed by N (0, 1). ASSUMPTION 1.2. Activation function σ (x) is a Lipschitz function with the Lipschitz constant λ σ ∈ (0, ∞). Assume that σ is centered and normalized with respect to ξ ∼ N (0, 1) such that

Furthermore, σ satisfies either of the following:

for some constants a, b, c ∈ R such that (9) holds.

Analogously to [39], our Assumption 1.2 permits σ to be the commonly used activation functions, including ReLU, Sigmoid, and Tanh, although we have to center and normalize the activation functions to guarantee (9). Such normalized activation functions exclude some trivial spike in the limiting spectra of CK and NTK [19,29]. The foregoing assumptions ensure our nonlinear Hanson-Wright inequality in the proof. As a future direction, going beyond Gaussian weights and Lipschitz activation functions may involve different types of concentration inequalities.

Next, we present the conditions of the deterministic input data X and the asymptotic regime for our main results. Define the following (ε, B)-orthonormal property for our data matrix X. DEFINITION 1.3. For given any ε, B > 0, matrix X is (ε, B)-orthonormal if for any distinct columns x α , x ³ in X, we have

and also

The empirical spectral distribution μ0 of X X converges weakly to a fixed and nondegenerate probability distribution μ 0 = µ 0 on [0, ∞).

In above (b), the (ε n , B)-orthonormal property with nε 4 n = o(1) is a quantitative version of pairwise approximate orthogonality for the column vectors of the data matrix X ∈ R d 0 ×n . When d 0 n, it holds, with high probability, for many random X with independent columns x α , including the anisotropic Gaussian vectors x α ∼ N (0, ) with tr( ) = 1 and 1/n, vectors generated by Gaussian mixture models, and vectors satisfying the log-Sobolev inequality or convex Lipschitz concentration property. See [29], Section 3.1, for more details. Specifically, when x α 's are independently sampled from the unit sphere S d 0 -1 , X is (ε n , B)orthonormal with high probability where ε n = O( log(n) n ) and B = O(1) as n d 0 . In this case, for any > 2, we have nε n → 0. In our theory, we always treat X as a deterministic matrix. However, our results also work for random input X independent of weights W and a by conditioning on the high probability event that X satisfies (ε n , B)-orthonormal property. Unlike data vectors with independent entries, our assumption is promising to analyze realworld datasets [53] and establish some n-dependent deterministic equivalents like [49].

The following Hermite polynomials are crucial to the approximation of in our analysis.

DEFINITION 1.5 (Normalized Hermite polynomials). The rth normalized Hermite polynomial is given by

Here {h r } ∞ r=0 form an orthonormal basis of L 2 (R, ), where denotes the standard Gaussian distribution. For σ 1 , σ 2 ∈ L 2 (R, ), the inner product is defined by

Every function σ ∈ L 2 (R, ) can be expanded as a Hermite polynomial expansion

where ζ r (σ ) is the rth Hermite coefficient defined by

In the following statements and proofs, we denote ξ ∼ N (0, 1). Then for any k ∈ N, we have

We define the innerproduct kernel random matrix f k (X X) ∈ R n×n by applying f k entrywise to X X. Define a deterministic matrix (11) 0

where the αth

We will employ 0 as an approximation of the population covariance in (3) in the spectral norm when nε 4 n → 0. For any n × n Hermitian matrix A n with eigenvalues λ 1 , . . . , λ n , the empirical spectral distribution of A is defined by

We write lim spec(A n ) = μ if μ A n → μ weakly as n → ∞. The main tool we use to study the limiting spectral distribution of a matrix sequence is the Stieltjes transform defined as follows.

DEFINITION 1.6 (Stieltjes transform). Let μ be a probability measure on R. The Stieltjes transform of μ is a function s(z) defined on the upper half plane C + by

For any n × n Hermitian matrix A n , the Stieltjes transform of the empirical spectral distribution of A n can be written as tr(A nz Id) -1 . We call (A nz Id) -1 the resolvent of A n .

2.1. Spectra of the centered CK and NTK. Our first result is a deformed semicircle law for the CK matrix. Denote by μ0 = (1b σ ) 2 + b 2 σ μ 0 the distribution of (1b 2 σ ) + b 2 σ λ with λ sampled from the distribution μ 0 . The limiting law of our centered and normalized CK matrix is depicted by μ s μ0 , where μ s is the standard semicircle law and the notation is the free multiplicative convolution in free harmonic analysis. For full descriptions of free independence and free multiplicative convolution, see [62], Lecture 18,and [5], Section 5.3.3. The free multiplicative convolution was first introduced in [79], which later has many applications for products of asymptotic free random matrices. The main tool for computing free multiplicative convolution is the S-transform, invented by [79]. S-transform was recently utilized to study the dynamical isometry of deep neural networks [25,37,65,66,84]. Some basic properties and intriguing examples for free multiplicative convolution with μ s can also be found in [15], Theorems 1.2, 1.3. THEOREM 2.1 (Limiting spectral distribution for the conjugate kernel). Suppose Assumptions 1.1, 1.2, and 1.4 of the input matrix X hold, the empirical eigenvalue distribution of

also converges weakly to the probability measure μ almost surely, whose Stieltjes transform m(z) is defined by

for each z ∈ C + , where ³(z) ∈ C + is the unique solution to

Suppose that we additionally have b σ = 0, that is, E[σ (ξ )] = 0. In this case, our Theorem 2.1 shows that the limiting spectral distribution of (2) is the semicircle law, and from (13), the deterministic data matrix X does not have an effect on the limiting spectrum. See Figure 1 for a cosine-type σ with b σ = 0. The only effect of the nonlinearity in μ is the coefficient b σ in the deformation μ0 . Appendix B with different activation functions. The red curves are implemented by the selfconsistent equations ( 15) and ( 16) in Theorem 2.1. In Section 4, we present general random matrix models with similar limiting eigenvalue distribution as μ whose Stieltjes transform is also determined by ( 15) and ( 16).

Theorem 2.1 can be extended to the NTK model as well. Denote by

As an approximation of in the spectral norm, we define

where f k 's are defined in (11), a σ is defined in (10), and the kth Hermite coefficient of σ is

Then, a similar deformed semicircle law can be obtained for the empirical NTK matrix H . THEOREM 2.2 (Limiting spectral distribution of the NTK). Under Assumptions 1.1, 1.2, and 1.4 of the input matrix X, the empirical eigenvalue distribution of (20)

weakly converges to μ almost surely, where 0 and 0 are defined in (11) and (18), respectively.

With our nonlinear Hanson-Wright inequality (Corollary 3.5), we attain the following concentration bound on the CK matrix in the spectral norm.

Then with probability at least

where C > 0 is a universal constant. REMARK 2.4. Theorem 2.3 ensures that the empirical spectral measure μ n of the centered random matrix d 1 n ( 1 d 1 Y Y -) has a bounded support for all sufficiently large n. Together with the global law in Theorem 2.1, our concentration inequality (22) is tight up to a constant factor. Additionally, by the weak convergence of μ n to μ proved in Theorem 2.1, we can take a test function x 2 to obtain that

Based on the foregoing estimation, we have the following lower bound on the smallest eigenvalue of the empirical conjugate kernel, denoted by λ min ( 1

THEOREM 2.5. Suppose Assumptions 1.1 and 1.2 hold and σ is not a linear function, X is (ε n , B)-orthonormal. Then with probability at least 1 -4e -2n ,

where C B is a constant depending on B. In particular, if

REMARK 2.6. A related result in [63], Theorem 5, assumed

. We relax the assumption on the column vectors of X, and extend the range of d 1 down to

) is lower bounded by an absolute constant, with an extra assumption that E[σ (ξ )] = 0. This assumption can always be satisfied by shifting the activation function with a proper constant. Our bound for λ min ( ) is derived via Hermite polynomial expansion, similar to [63], Lemma 15. However, we apply an ε-net argument for concentration bound for 1 d 1 Y Y around , while a matrix Chernoff concentration bound with truncation was used in [63], Theorem 5.

Additionally, the concentration for the NTK matrix H can be obtained in the next theorem. Recall that H is defined by (7) and the columns of S are defined by ( 8) with Assumption 1.1. THEOREM 2.7. Suppose d 1 ≥ log n, and σ is λ σ -Lipschitz. Then with probability at least

Moreover, if the assumptions in Theorem 2.3 hold, then with probability at least

REMARK 2.8. Compared to Proposition D.3 in [39], we assume a is a Gaussian vector instead of a Rademacher random vector and attain a better bound. If a i ∈ {+1, -1}, then one can apply matrix Bernstein inequality for the sum of bounded random matrices. In our case, the boundedness condition is not satisfied. Section S1.1 in [3] applied matrix Bernstein inequality for the sum of bounded random matrices when a is a Gaussian vector, but the boundedness condition does not hold in equation (S7) of [3].

Based on Theorem 2.7, we get a lower bound for the smallest eigenvalue of the NTK. THEOREM 2.9. Under Assumptions 1.1 and 1.2, suppose that X is (ε n , B)-orthonormal, σ is not a linear function, and d 1 ≥ log n. Then with probability at least 1n -7/3 ,

where C B is a constant depending only on B, and η k (σ ) is defined in (19). In particular, if

, and d 1 = ω(log n), then with high probability,

REMARK 2.10. We relax the assumption in [61] to d 1 = ω(log n) for the 2-layer case and our result is applicable beyond the ReLU activation function and to more general assumptions on X. Our proof strategy is different from [61]. In [61], the authors used the inequality λ min ((S S) (X X)) ≥ min i S i 2 2 • λ min (X X) where S i is the ith column of S. Then, getting the lower bound is reduced to show the concentration of the 2-norm of the column vectors of S. Here we apply a matrix concentration inequality to (S S) (X X) and gain a weaker assumption on d 1 to ensure the lower bound on λ min (H ). REMARK 2.11. In Theorems 2.5 and 2.9, we exclude the linear activation function. When σ (x) = x, it is easy to check both 1 d 1 λ min (Y Y ) and λ min (H ) will trivially determined by λ min (X X), which can be vanishing. In this case, the lower bounds of the smallest eigenvalues of CK and NTK rely on the assumption of μ 0 or the distribution of X. For instance, when the entries of X are i.i.d. Gaussian random variables, λ min (X X) has been analyzed in [75].

Training and test errors for random feature regression. We apply the results of the preceding sections to a two-layer neural network at random initialization defined in (1), to estimate the training errors and test errors with mean-square losses for random feature regression under the ultra-wide regime where d 1 /n → ∞ and n → ∞. In this model, we take the random feature 1

σ (W X) and consider the regression with respect to θ ∈ R d 1 based on

with training data X ∈ R d 0 ×n and training labels y ∈ R n . Considering the ridge regression with ridge parameter λ ≥ 0 and squared loss defined by

we can conclude that the minimization θ := arg min θ L(θ ) has an explicit solution

where Y = σ (W X) is defined in (5). When σ is nonlinear, by Theorem 2.5, it is feasible to take inverse in (26) for any λ ≥ 0. Hence, in the following results, we will focus on nonlinear activation functions. 2 In general, the optimal predictor for this random feature with respect to (25) is

where we define an empirical kernel

The n-dimension row vector is given by

and the (i, j ) entry of K n (X, X) is defined by

The optimal kernel predictor with a ridge parameter λ ≥ 0 for the kernel ridge regression is given by (see [10,18,41,46,51,71] for more details)

where K(X, X) is an n × n matrix such that its (i, j ) entry is K(x i , x j ), and K(x, X) is a row vector in R n similarly with (29). We compare the characteristics of the two different predictors f (RF ) λ (x) and f (K) λ (x) when the kernel function K is defined in (4). Denote the optimal predictors for random features and kernel K on training data X by (5). We aim to compare the training and test errors for these two predictors in ultra-wide random neural networks, respectively. Let training errors of these two predictors be

In the following theorem, we show that, with high probability, the training error of the random feature regression model can be approximated by the corresponding kernel regression model with the same ridge parameter λ ≥ 0 for ultra-wide neural networks. THEOREM 2.12 (Training error approximation). Suppose Assumptions 1.1, 1.2, and 1.4 hold, and σ is not a linear function. Then, for all large n, with probability at least 1 -4e -2n ,

where constants C 1 and C 2 only depend on λ, B and σ .

Next, to investigate the test errors (or generalization errors), we introduce further assumptions on the data and the target function that we want to learn from training data. Denote the true regression function by f * : R d 0 → R. Then, the training labels are defined by

where " ∈ R n is the training label noise. For simplicity, we further impose the following assumptions, analogously to [50]. ASSUMPTION 2.13. Assume that the target function is a linear function f * (x) = ³ ³ ³ * , x , where random vector satisfies ³ ³ ³ * ∼ N (0, σ 2 ³ ³ ³ Id). Then, in this case, the training label vector is given by y

, and a test data x ∈ R d 0 is independent with X and y such that X := [x 1 , . . . , x n , x] ∈ R d 0 ×(n+1) is also (ε n , B)-orthonormal. For convenience, we further assume the population covariance of the test data is

REMARK 2.15. Our Assumption 2.14 of the test data x ensures the same statistical behavior as training data in X, but we do not have any explicit assumption of the distribution of x. It is promising to adopt such assumptions to handle statistical models with real-world data [48,49]. Besides, it is possible to extend our analysis to general population covariance for E x [xx ].

For any predictor f , define the test error (generalization error) by

We first present the following approximation of the test error of a random feature predictor via its corresponding kernel predictor. THEOREM 2.16 (Test error approximation). Suppose that Assumptions 1.1, 1.2, 2.13, and 2.14 hold, and σ is not a linear function. Then, for any ε ∈ (0, 1/2), the difference of test errors satisfies

Theorems 2.12 and 2.16 verify that the random feature regression achieves the same asymptotic errors as the kernel regression, as long as n/d 1 → ∞. This is closely related to [55], Theorem 1, with different settings. Based on that, we can compute the asymptotic training and test errors for the random feature model by calculating the corresponding quantities for the kernel regression in the ultra-wide regime where n/d 1 → 0. THEOREM 2.17 (Asymptotic training and test errors). Suppose Assumptions 1.1 and 1.2 hold, and σ is not a linear function. Suppose the target function f * , training data X, and test data x ∈ R d 0 satisfy Assumptions 2.13 and 2.14. For any λ ≥ 0, let the effective ridge parameter be (36) λ eff (λ, σ )

If the training data has some limiting eigenvalue distribution μ 0 = lim spec X X as n → ∞ and n/d 0 → ´∈ (0, ∞), then when n/d 1 → 0 and n → ∞, the training error satisfies

and the test error satisfies

where the bias and variance functions are defined by

We emphasize that in the proof of Theorem 2.17, we also get n-dependent deterministic equivalents for training/test errors of the kernel regression to approximate the performance of random feature regression. This is akin to [49], Theorem 3, and [18], Theorem 4.13, but in different regimes. In the following Figure 3, we present implementations of test errors for random feature regressions on standard Gaussian random data and their limits (38). For simplicity, we fix n, d 0 , only let d 1 → ∞, and use empirical spectral distribution of X X to approximate μ 0 in B K (λ eff (λ, σ )) and V K (λ eff (λ, σ )), which is actually the n-dependent deterministic equivalent. However, for Gaussian random matrix X, μ 0 is actually a Marchenko-Pastur law with ratio ´, so B K (λ eff (λ, σ )) and V K (λ eff (λ, σ )) can be computed explicitly according to [50], Definition 1. REMARK 2.18 (Implicit regularization). For nonlinear σ , the effective ridge parameter (36) can be viewed as an inflated ridge parameter since b 2 σ ∈ [0, 1) and λ eff > λ ≥ 0. This λ eff leads to implicit regularization for our random feature and kernel ridge regressions even for the ridgeless regression with λ = 0 [18,41,46,57]. This effective ridge parameter λ eff also shows the effect of the nonlinearity in the random feature and kernel regressions induced by ultra-wide neural networks. REMARK 2.19. For convenience, we only consider the linear target function f * , but in general, the above theorems can also be obtained for nonlinear target functions, for instance, Test errors in solid lines with error bars are computed using an independent test set of size 5000. We average our results over 50 repetitions. Limiting test errors in black dash lines are computed by (38), and we take μ 0 to be the corresponding Marchenko-Pastur distributions.

f * is a nonlinear single-index model. Under (ε n , B)-orthonormal assumption with nε 4 n → 0, our expected kernel K(X, X) ≡ is approximated in terms of (39) lim spec

Id , whence, this kernel regression can only learn linear functions. So if f * is nonlinear, the limiting test error should be decomposed into the linear part as (38) and the nonlinear component as a noise [18], Theorem 4.13. For more conclusions of this kernel machine, we refer to [46,47,51,55]. REMARK 2.20 (Neural tangent regression). In parallel to the above results, we can obtain a similar analysis of the limiting training and test errors for random feature regression in (27) with empirical NTK given by either K n (X, X) = 1 d 1 (S S) (X X) or K n (X, X) = H . This random feature regression also refers to neural tangent regression [57]. With the help of our concentration results in Theorem 2.7 and the lower bound of the smallest eigenvalues in Theorem 2.9, we can directly extend the above Theorems 2.12, 2.16, and 2.17 to this neural tangent regression. We omit the proofs in these cases and only state the results as follows.

If K n (X, X) = 1 d 1 (S S) (X X) with expected kernel K(X, X) = defined by ( 17), the limiting training and test errors of this neural tangent regression can be approximated by the kernel regression with respect to , as long as d 1 = ω(log n). Analogously to (39), we have an additional approximation (40) lim spec = lim spec

Under the same assumptions of Theorem 2.17 and replacing n/d 1 → 0 with d 1 = ω(log n), we can conclude that the test error of this neural tangent regression has the same limit as (38) but changing the effective ridge parameter (36)

. This result is akin to [57], Corollary 3.2, but permits more general assumptions on X. The limiting training error of this neural tangent regression can be obtained by slightly modifying the coefficient in (37).

Similarly, if K n (X, X) = H defined by ( 7) possesses an expected kernel K(X, X) = + , this neural tangent regression in ( 27) is close to kernel regression (30) with kernel

under the ultra-wide regime, n/d 1 → 0. Combining ( 39) and ( 40), Theorem 2.17 can directly be extended to this neural tangent regression but replacing (36) with

. Section 6.1 of [3] also calculated this limiting test error when data X is isotropic Gaussian.

Organization of the paper. The remaining parts of the paper are structured as follows. In Section 3, we first provide a nonlinear Hanson-Wright inequality as a concentration tool for our spectral analysis. Section 4 gives a general theorem for the limiting spectral distributions of generalized centered sample covariance matrices. We prove the limiting spectral distributions for the empirical CK and NTK matrices (Theorem 2.1 and Theorem 2.2) in Section 5. Nonasymptotic estimates in Section 2.2 are proved in Section 6. In Section 7, we justify the asymptotic results of the training and test errors for the random feature model (Theorem 2.12 and Theorem 2.16). Auxiliary lemmas and additional simulations are included in Appendices A and B.

We give an improved version of Lemma 1 in [52] with a simple proof based on a Hanson-Wright inequality for random vectors with dependence [1]. This serves as the concentration tool for us to prove the deformed semicircle law in Section 5 and provide bounds on extreme eigenvalues in Section 6. We first define some concentration properties for random vectors. DEFINITION 3.1 (Concentration property). Let X be a random vector in R n . We say X has the K-concentration property with constant K if for any 1-Lipschitz function f : R n → R, we have E|f (X)| < ∞ and for any t > 0,

There are many distributions of random vectors satisfying K-concentration property, including uniform random vectors on the sphere, unit ball, hamming or continuous cube, uniform random permutation, etc. See [78], Chapter 5, for more details. DEFINITION 3.2 (Convex concentration property). Let X be a random vector in R n . We say X has the K-convex concentration property with the constant K if for any 1-Lipschitz convex function f : R n → R, we have E|f (X)| < ∞ and for any t > 0,

We will apply the following result from [1] to the nonlinear setting. LEMMA 3.3 (Theorem 2.5 in [1]). Let X be a mean zero random vector in R n . If X has the K-convex concentration property, then for any n × n matrix A and any t > 0,

for some universal constant C > 1. THEOREM 3.4. Let w ∈ R d 0 be a random vector with K-concentration property, X = (x 1 , . . . , x n ) ∈ R d 0 ×n be a deterministic matrix. Define y = σ (w X) , where σ is λ σ -Lipschitz, and = Eyy . Let A be an n × n deterministic matrix.

where C > 0 is an absolute constant.

2. If E[y] = 0, for any t > 0,

PROOF. Let f be any 1-Lipschitz convex function. Since y = σ (w X) , f (y) = f (σ (w X) ) is a λ σ X -Lipschitz function of w. Then by the Lipschitz concentration property of w in (41), we obtain

Therefore, y satisfies the Kλ σ X -convex concentration property. Define f (x) = f (x -Ey), then f is also a convex 1-Lipschitz function and f (y) = f (y -Ey). Hence ỹ := y -Ey also satisfies the Kλ σ X -convex concentration property. Applying Theorem 3.3 to ỹ, we have for any t > 0,

.

Since Eỹ = 0, the inequality above implies (42). Note that

Hence,

Since y (A + A )Ey is a (2 A Ey X λ σ )-Lipschitz function of w, by the Lipschitz concentration property of w, we have

Then combining (43), (44), and (45), we have

This finishes the proof.

Since the Gaussian random vector w ∼ N (0, I d 0 ) satisfies the K-concentration inequality with K = √ 2 (see, e.g., [20]), we have the following corollary.

COROLLARY 3.5. Let w ∼ N (0, I d 0 ), X = (x 1 , . . . , x n ) ∈ R d 0 ×n be a deterministic matrix. Define y = σ (w X) , where σ is λ σ -Lipschitz, and = Eyy . Let A be an n × n deterministic matrix.

for some absolute constant C > 0.

2. If E[y] = 0, for any t > 0,

where

REMARK 3.6. Compared to [52], Lemma 1, we identify the dependence on A F and Ey in the probability estimate. By using the inequality A F ≤ √ n A , we obtain a similar inequality to the one in [52] with a better dependence on n. Moreover, our bound in t 0 is independent of d 0 , while the corresponding term t 0 in [52], Lemma 1, depends on X and d 0 . In particular, when Eσ (ξ ) = 0 and X is (ε n , B)-orthonormal, t 0 is of order 1. Hence, (47) with the special choice of t 0 is the key ingredient in the proof of Theorem 2.3 to get a concentration of the spectral norm for CK. PROOF OF COROLLARY 3.5. We only need to prove (47), since other statements follow immediately by taking K = √ 2. Let x i be the ith column of X. Then

Let ξ ∼ N (0, 1). We have

Therefore

and ( 47) holds.

We include the following corollary about the variance of y Ay, which will be used in Section 5 to study the spectrum of the CK and NTK. COROLLARY 3.7. Under the same assumptions of Corollary 3.5, we further assume that t 0 ≤ C 1 n, and A , X ≤ C 2 . Then as n → ∞,

Thanks to Theorem 3.5 (2), we have that for any t > 0,

where constant C > 0 only relies on C 1 , C 2 , λ σ , and K. Therefore, we can compute the variance in the following way:

as n → ∞. Here, we use the dominant convergence theorem for the first integral in the last line.

Limiting law for general centered sample covariance matrices. Independent of the subsequent sections, this section focuses on the generalized sample covariance matrix where the dimension of the feature is much smaller than the sample size. We will later interpret such sample covariance matrix specifically for our neural network applications. Under certain weak assumptions, we prove the limiting eigenvalue distribution of the normalized sample covariance matrix satisfies two self-consistent equations, which are subsumed into a deformed semicircle law. Our findings in this section demonstrate some degree of universality, indicating that they hold across various random matrix models and may have implications for other related fields. -→ 0, and

Then the empirical eigenvalue distribution of matrix A n weakly converges to μ almost surely, whose Stieltjes transform m(z) is defined by

for each z ∈ C + , where ³(z) ∈ C + is the unique solution to

In particular, μ = μ s μ . REMARK 4.2. In [85], it was assumed that d n 3 E|y D n y -Tr D n n | 2 → 0, where n 3 /d → ∞ and n/d → 0 as n → ∞. By martingale difference, this condition implies (49). However, we are not able to verify a certain step in the proof of [85]. Hence, we will not directly adopt the result of [85] but consider a more general situation without assuming n 3 /d → ∞. The weakest conditions we found are conditions ( 49) and (50), which can be verified in our nonlinear random model.

The self-consistent equations we derived are consistent with the results in [16,85], where they studied the empirical spectral distribution of separable sample covariance matrices in the regime n/d → 0 under different assumptions. When n → ∞ and n/d → 0, our goal is to prove that the Stieltjes transform m n (z) of the empirical eigenvalue distribution of A n and ³ n (z) := tr[R(z) n ] pointwisely converges to m(z) and ³(z), respectively.

For the rest of this section, we first prove a series of lemmas to get n-dependent deterministic equivalents related to (51) and ( 52) and then deduce the proof of Theorem 4.1 at the end of this section. Recall

, and y is a random vector independent of A n with the same distribution of y i .

, where y j 's are independent copies of y defined in Theorem 4.1. Notice that, for any deterministic matrix D ∈ R n×n ,

Without loss of generality, we assume D ≤ 1. Taking normalized trace, we have

For each 1

where the leave-one-out resolvent R (i) is defined as

.

Hence, by (55), we obtain

Combining equations ( 54) and ( 56), and applying expectation at both sides implies

because of the assumption that all y i 's have the same distribution as vector y for all i ∈ [d +1].

With (57), to prove (53), we will first show that when n, d → ∞,

and spectral norms R , R(z) ≤ 1/v due to Proposition C.2 in [29]. Notice that n ≤ C. Hence, we can deduce that

as n → ∞. The same argument can be applied to the error of

Therefore ( 58) and ( 59) hold. For (60), we denote ỹ := y/(nd) 1/4 and observe that

ỹ 2 µ λ i . Then, we can control the real part of ỹ R(z)ỹ by Lemma A.4: (61) Re

We now separately consider two cases in the following:

• If the right-hand side of the above inequality (61) is at most 1/2, then

where we exploit the fact that (see also equation (A.1.11) in [13])

Finally, combining ( 62) and (63) in the above two cases, we can conclude the asymptotic result ( 60) because E y 2 = Tr n ≤ Cn in terms of the assumptions of Theorem 4.1.

Then with ( 58), (59), and (60), we get

as n → ∞. We utilize the notion E y to clarify the expectation only with respect to random vector y, conditioning on other independent random variables. So the conditional expectation is E y [ 1 n y DR(z)y] = tr DR(z) n and

Therefore, based on (64), the conclusion (53) holds.

In the next lemma, we apply the quadratic concentration condition (50) to simplify (53).

for each z ∈ C + and any deterministic matrix D with D ≤ C.

In other words, µ n can be expressed by

.

Observe that

Thus, µ n has the same expectation as the last term

, since we can first take the expectation for y conditioning on the resolvent R(z) and then take the expectation for R(z). Besides, notice that | Q1 |, | Q2 | ≤ C v uniformly. Hence, n d Q2 converges to zero uniformly and there exists some constant C > 0 such that (66) 1

for all large d and n. In addition, observe that

where ỹ is defined in the proof of Lemma 4.3. In terms of ( 62) and ( 63), we verify that (67)

where C > 0 is some constant depending on v. Next, recall that condition (50) exposes that ( 68)

as n → ∞. The first convergence is derived by viewing D n = R(z) and taking expectation conditional on R(z). To sum up, we can bound | n | based on ( 66) and ( 67) in the subsequent way:

Here, | tr n | ≤ n and n is uniformly bounded by some constant. Then, by Hölder's inequality, (68) implies that E[| n |] → 0, as n approaching to infinity. This concludes Likewise, in Lemma 4.5, we consider D = ( 3n (z

Because n is uniformly bounded, D ≤ C U . In terms of Lemma A.6, we only need to provide a lower bound for the imaginary part of U . Observe that Im U = Im 3n (z) n + v Id v Id since λ min ( n ) ≥ 0 and Im 3n (z) > 0. Thus, D ≤ Cv -1 for all n. Meanwhile, we have the equation 3n (z) n D = n -zD and hence,

So applying Lemma 4.5 again, we have another limiting equation tr D + 3n (z) → 0. In other words, (71) lim

Thanks to the identity

we can modify ( 70) and (71) to get (72) lim n,d→∞ " mn (z) + tr " 3n (z) n + z Id -1 = 0.

Since 3n (z) and mn (z) are uniformly bounded, for any subsequence in n, there is a further convergent sub-subsequence. We denote the limit of such sub-subsequence by ³(z) and m(z) ∈ C + respectively. Hence, by ( 71) and ( 72), one can conclude

Because of the convergence of the empirical eigenvalue distribution of n , we obtain the fixed point equation (52) for ³(z). Analogously, we can also obtain (51) for m(z) and ³(z). The existence and the uniqueness of the solutions to (51) and ( 52) are proved in [15], Theorem 2.1, and [80], Section 3.4, which implies the convergence of mn (z) and 3n (z) to m(z) and ³(z) governed by the self-consistent equations ( 51) and ( 52) as n → ∞, respectively. Then, by virtue of condition (49) in Theorem 4.1, we know m n (z) -mn (z) a.s.

-→ 0 and ³ n (z) -3n (z) a.s.

-→ 0. Therefore, the empirical Stieltjes transform m n (z) converges to m(z) almost surely for each z ∈ C + . Recall that the Stieltjes transform of μ is m(z). By the standard Stieltjes continuity theorem (see, e.g., [13], Theorem B.9), this finally concludes the weak convergence of empirical eigenvalue distribution of A n to μ. Now we show μ = μ s μ . The fixed point equations ( 51) and ( 52) induce (73) ³ 2 (z) + 1 + zm(z) = 0, since ³(z) ∈ C + for any z ∈ C + . Together with (51), we attain the same self-consistent equations for the convergence of the empirical spectral distribution of the Wigner-type matrix studied in [15], Theorem 1.1. Define W n , the n-by-n Wigner matrix, as a Hermitian matrix with independent entries

The Wigner-type matrix studied in [15], Definition 1.2, is indeed 1

n . Hence, such Wigner-type matrix 1

n has the same limiting spectral distribution as A n defined in Theorem 4.1. Both limits are determined by self-consistent equations ( 51) and (73).

On the other hand, based on [5], Theorem 5.4.5, 1

√ n W n and n are almost surely asymptotically free, that is, the empirical distribution of { 1 √ n W n , n } converges almost surely to the law of {s, d}, where s and d are two free noncommutative random variables (s is a semicircle element and d has the law μ ). Thus, the limiting spectral distribution μ of 1

is the free multiplicative convolution between μ s and μ . This implies μ = μ s μ in our setting.

To prove Theorem 2.1, we first establish the following proposition to analyze the difference between Stieltjes transform of ( 12) and its expectation. This will assist us to verify condition (49) in Theorem 4.1. The proof is based on [29], Lemma E.6. PROPOSITION 5.1. Let D ∈ R n×n be any deterministic symmetric matrix with a uniformly bounded spectral norm. Following the notions in Theorem 2.1, assume X ≤ C for some constant C and Assumption 1.2 holds. Let R(z) be the resolvent

, for any fixed z ∈ C + . Then, there exist some constants s, n 0 > 0 such that for all n > n 0 and any t > 0,

PROOF. Define function F : R d 1 ×d 0 → R by F (W ) := tr R(z)D. Fix any W, ∈ R d 1 ×d 0 where F = 1, and let W t = W + t . We want to verify F (W ) is a Lipschitz function in W with respect to the Frobenius norm. First, recall

where the last two terms are deterministic with respect to W . Hence, vec( )

where is the Hadamard product, and σ is applied entrywise. Here we utilize the formula

Therefore, based on the assumption of D, we have vec( )

for some constant C > 0. For the first term in the product on the right-hand side,

For the second term,

This holds for every such that F = 1, so F (W ) is C/ √ n-Lipschitz in W with respect to the Frobenius norm. Then the result follows from the Gaussian concentration inequality for Lipschitz functions.

Next, we investigate the approximation of = E w [σ (w X) σ (w X)] via the Hermite polynomials {h k } k≥0 . The orthogonality of Hermite polynomials allows us to write as a series of kernel matrices. Then we only need to estimate each kernel matrix in this series. The proof is directly based on [34], Lemma 2. The only difference is that we consider the deterministic input data X with the (ε n , B)-orthonormal property, while in Lemma 2 of [34], the matrix X is formed by independent Gaussian vectors. LEMMA 5.2. Recall the definition of 0 in (11). If X is (ε n , B)-orthonormal and Assumption 1.2 holds, then we have the spectral norm bound

where

PROOF. By Assumption 1.2, we know that

For any fixed t, σ (tx) ∈ L 2 (R, ). This is because σ (x) ∈ L 2 (R, ) is a Lipschitz function and by triangle inequality |σ (tx)σ (x)| ≤ λ σ |tx -x|, we have, for ξ ∼ N (0, 1),

For 1 ≤ α ≤ n, let σ α (x) := σ ( x α x) and the Hermite expansion of σ a can be written as

where the coefficient

where (ξ α , ξ ³ ) = (w u α , w u ³ ) is a Gaussian random vector with mean zero and covariance

.

By the orthogonality of Hermite polynomials with respect to and Lemma A.5, we can obtain

which leads to

For any k ∈ N, let T k be an n-by-n matrix with (α, ³)th entry

Specifically, for any k ∈ N, we have

. At first, we consider twice differentiable σ in Assumption 1.2. Similar with [34], equation (26), for any ε > 0 and |t -1| ≤ ε, we take the Taylor approximation of σ (tx) at point x, then there exists η between tx and x such that

Replacing x by ξ and taking expectation, since σ is uniformly bounded, we can get

where constant C does not depend on k. As for piecewise linear σ , it is not hard to see (79) E σ (tξ )σ (ξ ) = E σ (ξ )ξ (t -1). Now, we begin to approximate T k separately based on ( 77), (78), and (79). Denote diag(A) the diagonal submatrix of a matrix A.

(1) Approximation for k≥4 (T kdiag(T k )). At first, we estimate the L 2 norm with respect to of the function σ α . Recall that

Hence, σ α L 2 is uniformly bounded with some constant for all large n. Next, we estimate the off-diagonal entries of T k when k ≥ 4. From (76), we obtain that

when n is sufficiently large.

(2) Approximation for T 0 . Recall E[σ (ξ )] = 0 and by Gaussian integration by part,

Then, we have

If σ is twice differentiable, then E[σ (ξ )] = √ 2ζ 2 (σ ) as well. Thus, taking t = x α in ( 77) and ( 79) implies that for any 1 ≤ α ≤ n, (11). Then, (83)

Hence the difference between T 0 and μμ is controlled by

(3) Approximation for T k for k = 1, 2, 3. For 0 ≤ k ≤ 3, Assumption 1.4 and (78) show that

when n is sufficiently large. Notice that T k = D k f k (X X)D k , where D k is the diagonal matrix. Hence, by (86),

And for k = 1, 2, 3, by the triangle inequality,

From Lemma A.1 in Appendix A, we have that

So the left-hand side of ( 87) is bounded. Analogously, we can verify f 3 (X X) is also bounded. Therefore, we have

for some constant C and k = 1, 2, 3 when n is sufficiently large. (4) Approximation for k≥4 diag(T k ). Since u α u α = 1, we know

.

First, by ( 80) and ( 81), we can claim that

Second, in terms of (86), we obtain

for k = 1, 2 and 3. Combining these together, we conclude that

In terms of approximations ( 82), ( 85), (88), and (89), we can finally manifest

for some constant C > 0 as √ nε 2 n → 0. The spectral norm bound of is directly deduced by the spectral norm bound of 0 based on ( 84) and ( 87), together with (90). REMARK 5.3 (Optimality of ε n ). For general deterministic data X, our pairwise orthogonality assumption with rate nε 4 n = o(1) is optimal for the approximation of by 0 in the spectral norm. If we relax the decay rate of ε n in Assumption 1.4, the above approximation may require including terms of higher-degree f k (X X) for k ≥ 4 in 0 , which will lead to the invalidation of some of our following results and simplifications. Subsequent to the initial completion of our paper, this weaker regime has been considered in our follow-up work [81].

Next, we continue to provide an additional estimate for , but in the Frobenius norm to further simplify the limiting spectral distribution of . LEMMA 5.4. If Assumptions 1.2 and 1.4 hold, then has the same limiting spectrum as

PROOF. By the definition of b σ , we know that b σ = ζ 1 (σ ). As a direct deduction of Lemma 5.2, the limiting spectrum of is identical to the limiting spectrum of 0 . To prove this lemma, it suffices to check the Frobenius norm of the difference between 0 and

By the definition of vector μ and the assumption of X, we have

For k = 2, 3, the Frobenius norm can be controlled by

""

Hence, as n → ∞, we have

Hence, lim spec is the same as lim spec

Moreover, the proof of Lemma 5.4 can be modified to prove (40), so we omit its proof. Now, based on Corollary 3.7, Proposition 5.1, Lemma 5.2, and Lemma 5.4, applying Theorem 4.1 for general sample covariance matrices, we can finish the proof of Theorem 2.1. PROOF OF THEOREM 2.1. Based on Corollary 3.7 and Proposition 5.1, we can verify the conditions ( 49) and ( 50) in Theorem 4.1. By Lemma 5.2 and Lemma 5.4, we know that the limiting eigenvalue distributions of and (1b 2 σ ) Id +b 2 σ X X are identical and is uniformly bounded. So the limiting eigenvalue distribution of denoted by μ is just

Hence, the first conclusion of Theorem 2.1 follows from Theorem 4.1. For the second part of this theorem, we consider the difference

where we employ Lemma 5.2 and the assumption d 1 ε 4 n = o(1). Thus, because of Lemma A.7,

) has the same limiting eigenvalue distribution as (12), μ s ((1b 2 σ ) + b 2 σ μ 0 ). This finishes the proof of Theorem 2.1.

Next, we move to study the empirical NTK and its corresponding limiting eigenvalue distribution. Similarly, we first verify that such NTK concentrates around its expectation and then simplify this expectation by some deterministic matrix only depending on the input data matrix X and nonlinear activation σ . The following lemma can be obtained from ( 23) in Theorem 2.7. LEMMA 5.5. Suppose that Assumption 1.1 holds, sup x∈R |σ (x)| ≤ λ σ and X ≤ B.

n n, where and 0 are defined in (17) and (18), respectively, and C B is a constant depending on B.

PROOF. We can directly apply methods in the proof of Lemma 5.2. Notice that ( 6) and ( 8) imply

for any standard Gaussian random vector w ∼ N (0, Id). Recall that (19) defines the kth coefficient of Hermite expansion of σ (x) by η k (σ ) for any k ∈ N. Then, Assumption 1.2 indicates b σ = η 0 (σ ) and a σ = ∞ k=0 η 2 k (σ ). For 1 ≤ α ≤ n, we introduce φ α (x) := σ ( x α x) and the Hermite expansion of this function as

where (ξ α , ξ ³ ) = (w u α , w u ³ ) is a Gaussian random vector with mean zero and covariance (74). Following the derivation of formula (75), we obtain

For any k ∈ N, let T k ∈ R n×n be an n-by-n matrix with (α, ³) entry

We can write

Then, adopting the proof of (88), we can similarly conclude that

for some constant C and k = 0, 1, 2, when n is sufficiently large. Likewise, (82) indicates

and a similar proof of (89) implies that

Based on these approximations, we can conclude the result of this lemma.

PROOF OF THEOREM 2.2. The first part of the statement is a straight consequence of (91) and Theorem 2.1. Denote by

Hence, (91) indicates B -A → 0 as n → ∞. This convergence implies that limiting laws of A and B are identical because of Lemma A.3. The second part is because of Lemma 5.2 and Lemma 5.6. From ( 7) and ( 17), E[H ] = + . Then almost surely,

as ε 4 n d 1 → 0 by the assumption of Theorem 2.2. Therefore, the limiting eigenvalue distribution of ( 21) is the same as (20).

6. Proof of the concentration for extreme eigenvalues. In this section, we obtain the estimates of the extreme eigenvalues for the CK and NTK we studied in Section 5. The limiting spectral distribution of 1

) tells us the bulk behavior of the spectrum. An estimation of the extreme eigenvalues will show that the eigenvalues are confined in a finite interval with high probability. We first provide a nonasymptotic bound on the concentration of 1 d 1 Y Y under the spectral norm. The proof is based on the Hanson-Wright inequality we proved in Section 3 and an ε-net argument.

PROOF OF THEOREM 2.3. Recall notation in Section 1. Define

where y i = σ (w i X).

For any fixed z ∈ S n-1 , we have

where

and column vector (y 1 , . . . , y d 1 ) ∈ R nd 1 is the concatenation of column vectors y 1 , . . . , y d 1 . Then (y 1 , . . . ,

Notice that

Denote ỹ = (y 1 , . . . , y d 1 ). With (48), we obtain

where the last line is from the assumptions on X and σ . When B = 0, applying (47) to (92) implies P " (y 1 , . . . , y d 1 ) A z (y 1 , . . . ,

Let N be a 1/2-net on S n-1 with |N | ≤ 5 n (see, e.g., [78], Corollary 4.2.13), then

Taking a union bound over N yields

to conclude

the upper bound in ( 22) is then verified. When B = 0, we can apply (46) and follow the same steps to get the desired bound.

By the concentration inequality in Theorem 2.3, we can get a lower bound on the smallest eigenvalue of the conjugate kernel 1 d 1 Y Y as follows.

LEMMA 6.1. Assume X satisfies n i=1 ( x i 2 -1) 2 ≤ B 2 for a constant B > 0, and σ is λ σ -Lipschitz with Eσ (ξ ) = 0. Then with probability at least 1 -4e -2n ,

PROOF. By Weyl's inequality [5], Corollary A.6, we have

Then (93) follows from (22).

The lower bound in (93) relies on λ min ( ). Under certain assumptions on X and σ , we can guarantee that λ min ( ) is bounded below by an absolute constant. LEMMA 6.2. Assume σ is not a linear function and σ (x) is Lipschitz. Then

PROOF. Suppose that sup{k ∈ N : ζ k (σ ) 2 > 0} is finite. Then σ is a polynomial of degree at least 2 from our assumption, which is a contradiction to the fact that σ is Lipschitz. Hence, (94) holds. LEMMA 6.3. Assume Assumption 1.2 holds, σ is not a linear function, and X satisfies (ε n , B)-orthonormal property. Then,

REMARK 6.4. This bound will not hold when σ is a linear function. Suppose σ is a linear function, under Assumption 1.2, we must have σ (x) = x and = X X. Then we will not have a lower bound on λ min ( ) based on the Hermite coefficients of σ .

PROOF OF LEMMA 6.3. From Lemma 5.2, under our assumptions, we know that

where 0 is given by (11). Thus, λ min ( ) ≥ λ min ( 0 ) -C B ε 2 n √ n, and, from Weyl's inequality [5], Theorem A.5, we have

0 ×n , and each column of K k is given by the kth Kronecker product

Since σ is nonlinear and Lipschitz, (94) holds for σ . Therefore,

and (95) holds. Theorem 2.5 then follows directly from Lemma 6.1 and Lemma 6.3. Next, we move on to nonasymptotic estimations for NTKs. Recall that the empirical NTK matrix H is given by (7) and the αth column of S is defined by diag(σ (W x α ))a, for 1 ≤ α ≤ n, in (8).

The ith row of S is given by z i := σ (w i X)a i , and E[z i ] = 0, where a i is the ith entry of a. Define D α = diag(σ (w α X)a α ), for 1 ≤ α ≤ d 1 . We can rewrite (S S) (X X) as

Let us define L and further expand it as follows:

Here Z i is a centered random matrix, and we can apply matrix Bernstein's inequality to show the concentration of L. Since Z i does not have an almost sure bound on the spectral norm, we will use the following sub-exponential version of the matrix Bernstein inequality from [77]. LEMMA 6.5 ([77], Theorem 6.2). Let Z k be independent Hermitian matrices of size n × n. Assume

for any integer p ≥ 2. Then for all t ≥ 0,

PROOF OF THEOREM 2.7. From (97), EZ i = 0, and

where C 1 = λ 2 σ X 2 and where a i ∼ N (0, 1) is the ith entry of the second layer weight a. Then

So we can take R = 2C 2 1 , a 2 = 8C 4 1 in (98) and obtain

.

Hence, L defined in (96) has a probability bound:

.

Under the assumption that d 1 ≥ log n, we conclude that, with high probability at least 1n -7/3 ,

Thus, as a corollary, the two statements in Lemma 5.5 follow from (99). Meanwhile, since

the bound in (24) follows from Theorem 2.3 and (99).

We now proceed to provide a lower bound of λ min (H ) from Theorem 2.7.

PROOF OF THEOREM 2.9. Note that from ( 7), (17), and (96), we have

Then with Lemma 5.6, we can get

Therefore, from Theorem 2.7, with probability at least 1n -7/3 ,

Since σ is Lipschitz and nonlinear, we know σ (x) is not a linear function (including the constant function) and |σ (x)| is bounded. Suppose that σ (x) has finite many nonzero Hermite coefficients, σ (x) is a polynomial, then we get a contradiction. Hence, the Hermite coefficients of σ satisfy

This finishes the proof.

7. Proofs of Theorem 2.12 and Theorem 2.17. By definitions, the random matrix K n (X, X) is 1 d 1 Y Y and the kernel matrix K(X, X) = is defined in (3). These two matrices have already been analyzed in Theorem 2.3 and Theorem 2.5, so we will apply these results to estimate how great the difference between training errors of random feature regression and its corresponding kernel regression. PROOF OF THEOREM 2.12. Denote K λ := (K + λ Id). From the definitions of training errors in (31) and (32), we have

(100)

Here, in (100), we employ the identity (101)

for A = (K(X, X) + λ Id) -2 and B = (K n (X, X) + λ Id) -2 , and the fact that (K(X, X) + λ Id) -1 ≤ λ -1 min (K(X, X)) and (K n (X, X) + λ Id) -1 ≤ λ -1 min (K n (X, X)). Next, before providing uniform upper bounds for λ -2 min (K(X, X)) and λ -2 min (K n (X, X)) in (100), we can first get a bound for the last term of (100) as follows:

for some constant C > 0, with probability at least 1 -4e -2n , where the last bound in (102) is due to Theorem 2.3 and Lemma A.9 in Appendix A. Additionally, combining Theorem 2.3 and Theorem 2.5, we can easily get (103)

for all large n and some universal constant C, under the same event that (102) holds. Theorem 6.3 also shows λ -1 min (K(X, X)) ≤ C for all large n. Hence, with the upper bounds for λ -2 min (K(X, X)) and λ -2 min (K n (X, X)), (33) follows from the bounds of (100) and (102).

For ease of notation, we denote K := K(X, X) and K n := K n (X, X). Hence, from (34), we can further decompose the test errors for K and K n into

As we can see, to compare the test errors between random feature and kernel regression models, we need to control |E 1 -Ē1 | and |E 2 -Ē2 |. First, it is necessary to study the concentrations of

LEMMA 7.1. Under Assumption 1.2 for σ and Assumption 2.14 for x and X, with probability at least 1 -4e -2n , we have

where C > 0 is a universal constant. Here, we only consider the randomness of the weight matrix in K n (x, X) defined by (28) and (29).

PROOF. We consider X = [x 1 , . . . , x n , x], its corresponding kernels K n ( X, X), and K( X, X) ∈ R (n+1)×(n+1) . Under Assumption 2.14, we can directly apply Theorem 2.3 to get the concentration of K n ( X, X) around K( X, X), namely,

with probability at least 1 -4e -2n . Meanwhile, we can write K n ( X, X) and K( X, X) as block matrices:

Since the 2 -norm of any row is bounded above by the spectral norm of its entire matrix, we complete the proof of (106).

LEMMA 7.2. Assume that training labels satisfy Assumption 2.13 and X ≤ B, then for any deterministic A ∈ R n×n , we have

F , where constant c only depends on σ ³ ³ ³ , σ ", and B. Moreover,

PROOF. We follow the idea in Lemma C.8 of [56] to investigate the variance of the quadratic form for the Gaussian random vector by (108) Var

F , for any deterministic square matrix A and standard normal random vector g. Notice that the quadratic form

where g is a standard Gaussian random vector in R d 0 +n . Similarly, the second quadratic form can be written as

and similarly Ã2 F ≤ c A 2 F for a constant c, we can complete the proof.

As a remark, in Lemma 7.2, for simplicity, we only provide a variance control for the quadratic forms to obtain convergence in probability in the following proofs of Theorems 2.16 and 2.17. However, we can apply Hanson-Wright inequalities in Section 3 to get more precise probability bounds and consider non-Gaussian distributions for ³ ³ ³ * and ".

PROOF OF THEOREM 2.16. Based on the preceding expansions of L( f (RF ) λ (x)) and L( f (K) λ (x)) in ( 104) and (105), we need to control the right-hand side of

In the subsequent procedure, we first take the concentrations of E 1 and E 2 with respect to normal random vectors ³ ³ ³ * and ", respectively. Then, we apply Theorem 2.3 and Lemma 7.1 to complete the proof of (35). For simplicity, we start with the second term

where I 1 and I 2 are quadratic forms defined below

and their expectations with respect to random vectors ³ ³ ³ * and " are denoted by

We first consider the randomness of the weight matrix in K n and define the event E where both (103) and (107) hold. Then, Theorem 2.5 and the proof of Lemma 7.1 indicate that event E occurs with probability at least 1 -4e -2n for all large n. Notice that E does not rely on the randomness of test data x.

We now consider A = E x [x(K n (x, X) -K(x, X))](K n + λ Id) -1 in Lemma 7.2. Conditioning on event E, we have

for some constant C, where we utilize the assumption E[ x 2 ] = 1. Hence, based on Lemma 7.2, we know Var ",³ ³ ³ * (I 1 ) ≤ cn/d 1 , for some constant c. By Chebyshev's inequality and event E , (111)

), we can apply ( 101) and Moreover, conditioning on the event E ,

for some constant C. In the same way, with (112

Analogously, the first term | Ē1 -E 1 | is controlled by the following four quadratic forms

where we define by J i := y A i y for 1 ≤ i ≤ 4 and

Similarly with (110) and (113), it is not hard to verify

d 1 conditioning on the event E . Then, like (111), we can invoke Lemma 7.2 for each A i to apply Chebyshev's inequality and conclude | Ē1 -E 1 | = o((n/d 1 ) 1/2-ε ) with probability 1o(1) when d 1 /n → ∞, for any ε ∈ (0, 1/2). LEMMA 7.3. With Assumptions 1.2 and 2.14, for (ε n , B)-orthonormal X, we have that

for some constant C > 0. PROOF. By Lemma A.8, we have an entrywise approximation

Then, we can verify (114) based on the following approximation

for some universal constant C. The same argument can also be employed to prove (115), so details will be omitted here.

PROOF OF THEOREM 2.17. From ( 33) and ( 35), we can easily conclude that train under the ultra-wide regime.

Recall that K λ = (K + λ Id) and the test error is given by

where

The spectral norm of K λ is bounded from above and the smallest eigenvalue is bounded from below by some positive constants.

We first focus on the last two terms L 1 and L 2 in the test error. Let us define

To obtain the last approximation, we define K(X, X) := b 2 σ X X + (1b 2 σ ) Id and (119) Kλ := b 2 σ X X + " 1 + λb 2 σ Id . We aim to replace K λ by Kλ in L1 and L2 . Recalling the identity (101), we have

K(X, X) -K(X, X) K -1 λ . Since σ is not a linear function, 1b 2 σ > 0. Then, with (103), the proof of Lemma 5.4 indicates

where we apply the fact that λ min ( K(X, X)) ≥ 1b 2 σ > 0. Let us denote Then, with the help of (120) and uniform bounds of the spectral norms of X X, K -1 λ and K-1 λ , we obtain that

as n → ∞, n/d 0 → ´and nε 4 n → 0. Combining all the approximations, we conclude that L i and L 0 i have identical limits in probability for i = 1, 2. On the other hand, based on the assumption of X and definitions in (119), (121), and (122), it is not hard to check that

Therefore, L 1 and L 2 converge in probability to the above limits, respectively, as n → ∞. In the end, we apply the concentration of the quadratic form ³ ³ ³ * ³ ³ ³ * in (118) to get 1 d 0 ³ ³ ³ * 2 P -→ σ 2 ³ ³ ³ . Then, by (117), we can get the limit in (38) for the test error L( f (RF ) λ

). As a byproduct, we can even use L 0 1 and L 0 2 to form an n-dependent deterministic equivalent of L( f (RF ) λ ) as well.

Thanks to Lemma 7.2, the training error, E (K,λ) train = λ 2 n y K -2 λ y, analogously, concentrates around its expectation with respect to ³ ³ ³ * and ", which is σ 2 ³ ³ ³ λ 2 tr K -2 λ X X + σ 2 " λ 2 tr K -2 λ . Moreover, because of (120), we can further substitute K -2 λ by K-2 λ defined in (119). Hence, we know that, asymptotically, E (K,λ) trainσ 2 ³ ³ ³ λ 2 tr K-2 λ X Xσ 2 " λ 2 tr K-2 λ P -→ 0, where as n, d 0 → ∞,

σ ) 2 dμ 0 (x). The last two limits are due to μ 0 = lim spec X X as n, d 0 → ∞. Therefore, by (116), we obtain our final result (37) in Theorem 2.17. APPENDIX A: AUXILIARY LEMMAS LEMMA A.1 (Equation (3.7.9) in [43]). Let A, B be two n × n matrices, A be positive semidefinite, and A B be the Hadamard product between A and B. Then,

LEMMA A.2 (Sherman-Morrison formula, [17]). Suppose A ∈ R n×n is an invertible square matrix and u, v ∈ R n are column vectors. Then

LEMMA A.3 (Theorem A.45 in [13]). Let A, B be two n × n Hermitian matrices. Then A and B have the same limiting spectral distribution if A -B → 0 as n → ∞.

LEMMA A.4 (Theorem B.11 in [13]). Let z = x + iv ∈ C, v > 0 and s(z) be the Stieltjes transform of a probability measure. Then | Re s(z)| ≤ v -1/2 √ Im s(z).

LEMMA A.5 (Lemma D.2 in [60]). Let x, y ∈ R d such that x = y = 1 and w ∼ N (0, I d ). Let h j be the j th normalized Hermite polynomial given in (1.5). Then PROOF. When σ is twice differentiable in Assumption 1.2, this result follows from Lemma D.3 in [29]. When σ is a piecewise linear function defined in case 2 of Assumption 1.2, the second inequality follows from (79) with t = x α . For the first inequality, the Hermite expansion of α³ is given by (75) "

for ε ∈ (0, 1) and (ε, B)-orthonormal X. This completes the proof of this lemma.

With the above lemma, the proof of Lemma D.4 in [29] yields the following lemma.

LEMMA A.9. Under the same assumptions as Lemma A.8, there exists a constant C such that K(X, X) ≤ CB 2 . Additionally, with Assumption 2.14, we have K( X, X) ≤ CB 2 . APPENDIX B: ADDITIONAL SIMULATIONS Figures 4 and 5 provide additional simulations for the eigenvalue distribution described in Theorem 2.1 with different activation functions and scaling. Here, we compute the empirical eigenvalue distributions of centered CK matrices in histograms and the limiting spectra in terms of self-consistent equations. All the input data X's are standard random Gaussian matrices. Interestingly, in Figure 5, we observe an outlier that emerges outside the bulk distribution for the piecewise linear activation function defined in case 2 of Assumption 1.2. The analysis of the emergence of the outlier, in this case, would be interesting for future work. Acknowledgments. Z.W. would like to thank Denny Wu for his valuable suggestions and comments. Both authors would like to thank Lucas Benigni, Ioana Dumitriu, and Kameron Decker Harris for their helpful discussion. Funding. Z.W. is partially supported by NSF DMS-2055340 and NSF DMS-2154099. This material is based upon work supported by the National Science Foundation under Grant No. DMS-1928930 while Y.Z. was in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the Fall 2021 semester for the program "Universality