Year: 2019

Abstract:
One of the distinguishing characteristics of modern deep learning systems is that they typically employ neural network architectures that utilize enormous numbers of parameters, often in the millions and sometimes even in the billions. While this paradigm has inspired significant research on the properties of large networks, relatively little work has been devoted to the fact that these networks are often used to model large complex datasets, which may themselves contain millions or even billions of constraints. In this work, we focus on this high-dimensional regime in which both the dataset size and the number of features tend to infinity. We analyze the performance of a simple regression model trained on the random features F = f(WX + B) for a random weight matrix W and random bias vector B, obtaining an exact formula for the asymptotic training error on a noisy autoencoding task. The role of the bias can be understood as parameterizing a distribution over activation functions, and our analysis directly generalizes to such distributions, even those not expressible with a traditional additive bias. Intriguingly, we find that a mixture of nonlinearities can outperform the best single nonlinearity on the noisy autoecndoing task, suggesting that mixtures of nonlinearities might be useful for approximate kernel methods or neural network architecture design.

heading: 1 INTRODUCTION
text: It is undeniable that in recent years deep learning systems have found widespread success in their applications to a diverse and ever-expanding set of domains. The foundational results on many tasks such as image recognition (Krizhevsky et al., 2012), speech recognition (Hinton et al., 2012), and machine translation (Wu et al., 2016), have begun to make their way into higher-level products that people interact with and rely upon in their daily lives. Whether these products generate a medical diagnosis, a navigation decision, or some other important output, it is crucial to understand the inner-workings of the deep learning algorithms that generate them.
Unfortunately, our theoretical understanding of these deep learning algorithms continues to lag behind their impressive practical successes. One main challenge in building a fuller understanding stems from the fact that deep neural networks are complex nonlinear functions that employ millions or even billions (Shazeer et al., 2017) of parameters. Traditional wisdom would suggest that to this parameter complexity corresponds to an optimization difficulty. Recent work, however, suggests that as the width of a network’s hidden layers becomes large, the loss function simplifies and a theoretical analysis becomes tractable (Jacot et al., 2018; Mei et al., 2018; Chizat & Bach, 2018a; Mei et al., 2019; Rotskoff et al., 2019; Rotskoff & Vanden-Eijnden, 2018). In some scenarios, the simplification is such that throughout training the parameters of the model stay within an infinitesimal radius of their initial (random) values, implying that much about neural network training can be understood by studying the random initialization (Jacot et al., 2018; Chizat & Bach, 2018b; Lee et al., 2019).
Another main challenge in building a rich understanding of deep learning systems stems from the fact that they are often trained on very large, complex datasets: even if the models themselves are very large, they may not be large in comparison to the number of constraints they are designed to satisfy. Indeed, many important phenomena may become apparent only by examining the high-dimensional regime in which the dataset size and width are both large and of the same order.
In this work, we focus on the high-dimensional regime and analyze the performance of a regression model trained on the random features F = f(WX +B) for a random weight matrix W and random
bias vector B. We obtain an exact formula for the training error on a noisy autoencoding task in the limit that the width and dataset size both go to infinity. The result is determined by the resolvent of the kernel matrix F tF , whose properties we analyze via the resolvent method from random matrix theory. Our analysis also provides an exact formula for the eigenvalue density of the kernel matrix, which may be of independent interest since it provides a characterization for how spectral properties of the data covariance matrix propagate through neural network layers at initialization.

heading: 1.1 OUR CONTRIBUTIONS
text: The main contribution of our work is an exact characterization of the training error of a ridgeregularized random feature regression model on a noisy autoencoder task in the high-dimensional regime. This is one of the first non-trivial models to be solved exactly in the joint limit of large data and large width and provides an interesting testing ground in which to analyze this regime. Some of our additional contributions include,
• An exact characterization of the spectral density of the randm feature matrix F = f(WX + B), extending prior results of Pennington & Worah (2017) to non-Gaussian data distributions and to non-zero bias distributions. • One interpretation of the random additive bias is that it induces a distribution of activation
functions parameterized by B, i.e. f(Z;B) := f(Z +B). Our analysis trivially extends to any distribution of activation functions f(· ;B) parameterized by B. • We show that there exists a non-trivial distribution over activation functions that outperforms
the best possible single activation function on a noisy autoencoding task. • Our method of proof introduces a surrogate “linearization" of F , F lin, that possesses the
same spectral information as F . F lin and its properties are likely to be of further interest and utility in analyzing neural networks in high dimensions.

heading: 1.2 RELATED WORK
text: Neural networks have been studied from the perspective of high-dimensional statistics in a number of recent works. Most prior work has focused on the bias-free case. Pennington & Worah (2017) studied the spectrum of the activation matrix f(WX) for iid Gaussian data and derived an analytic expression for the training error of a ridge-regularized random feature model trained on pure noise. It is natural to consider incorporating biases by appending a constant feature 1 to X . Unfortunately, this leads to biases that are the same order as the weights, and so the effect disappears in the large dataset limit. Moreover, this modification on the data violates the assumptions of Pennington & Worah (2017). Hastie et al. (2019) study ridgeless interpolation in high-dimensional interpolation for linear features as well as nonlinear random features of iid Gaussian data. Louart et al. (2018) derived a deterministic equivalent for the resolvent of the kernel matrix F tF which allowed for a characterization of the asymptotic training and test performance of linear ridge regression of random feature models.
Other work has investigated learning dynamics and generalization in the high-dimensional regime (Liao & Couillet, 2018a; Lampinen & Ganguli, 2018; Advani & Ganguli, 2016; Advani & Saxe, 2017) as well as the spectra of more complicated objects such as the Hessian (Pennington & Bahri, 2017) and Fisher information matrix (Pennington & Worah, 2018). From the mathematical perspective, random matrix theory provides natural tools (Silverstein & Bai, 1995a) for analyzing the behavior of neural networks in the high-dimensional regime. Liao & Couillet (2018b) examined spectra for data drawn from Gaussian mixture models; see also (El Karoui, 2010) on the spectra of random kernel matrices.

heading: 2 PRELIMINARIES
text: Consider a dataset X ∈ Rn0×m and the random feature matrix, F = f(WX;B),
generated by a single hidden-layer network with iid Gaussian weights W ∈ Rn1×n0 (Wak ∼ N (0, σ2W /n0)), activation function f , and biases B = b1mT ∈ Rn1×m (for b ∈ Rn1). We regard
the second argument of f as parametrizing (continuously or discretely) an ensemble of activation functions. We refer to B (or b) as the bias, in reference to the important special case f(WX + B) (additive bias). In general, however, we only assume the parameters ba are such that the measure 1n1 ∑ a δba converges in distribution to some limiting distribution µB , effecting an arbitrary distribution over activation functions.
We assume that E |f(N ; b)|k for N ∼ N (0, σ) is finite for all 1 ≤ k ≤ 3, σ > 0, and b ∈ support(µB). When µB is a single Dirac mass at location b0, the activation function can be written as f(WX;B) = f(WX; b0) = g(WX) for some single-argument function g (for an additive bias, g(WX) = f(WX + b0)). When this is the case, we say the model has a single activation function, as opposed to a mixture or distribution of activation functions.
The quantity of interest for the investigations below is the kernel matrix 1n1F tF , and in particular its resolvent,
G(z) = (
1 n1 F tF − zI
)−1 . (1)
As we review in Sec. 4, the optimal regression coefficients of a linear model on the random features F is a simple function of this resolvent.
The high-dimensional regime that we study is the one in which the dataset size m, feature dimensionality n0, and hidden layer width n1 all go to infinity at the same rate. In particular, as is standard in the random matrix literature, we assume that we can parameterize the limit in terms of the dataset size m in such a way that there exist two positive constants,
φ := lim m→∞
n0(m)
m and ψ := lim m→∞
n0(m) n1(m) . (2)
Note that the resolvent is a random matrix, but as m grows large, its normalized trace becomes a deterministic quantity. In the limit that m→∞, this quantity is known as the Stieltjes transform1,
m(z) := lim m→∞
1 m trG(z) . (3)
Together with an auxiliary transform m̃(z), defined below, these deterministic quantities completely characterize the asymptotic training error of kernel ridge regression on a noisy autoencoder task in this high-dimensional regime.
The Stieltjes transform frequently arises in random matrix methods as a way to encode the spectra of matrices. In particular, if λi are the eigenvalues of 1n1F
tF and the empirical distribution of eigenvalues converges in distribution to some deterministic limiting density as m→∞,
1
m m∑ i=1 δλi → µ(λ)dλ , (4)
then (with appropriate technical assumptions), the limiting spectral density itself can be recovered from the Stieltjes transform m(z) via the inversion formula,
µ(λ) = lim →0+ m(λ− i )−m(λ+ i ) 2πi . (5)
The Stieltjes transform then substitutes convergence in distribution for pointwise convergence for all z such that =z > 0.

heading: 2.1 METHODS FOR COMPUTING THE STIELTJES TRANSFORM
text: We briefly review two standard methods for computing the Stieltjes transform m(z), the resolvent method and the moments method.
The resolvent method is an approach for computing the Stieltes transform based on the application of the Schur complement formula to the resolvent itself (or to a closely-related block matrix). Intuitively, as the matrix size becomes large, the minors of the matrix are similar in distribution to the larger matrix, and, moreover, the Cauchy interlacing theorem guarantees that their Stieltjes transforms are
1The distinction between the Stieltjes transform m(z) and the dataset size m should be clear from context.
close as well. This allows for the derivation of a self-consistent equation (SCE) in which the Stieltjes transform appears on the left-hand side as the trace of the resolvent, and on the right-hand side as the trace of one of its minors.
The moments method is more combinatorial in nature and involves expanding the resolvent for large z and computing the traces of each term,
m(z) = lim m→∞ 1 m trG(z) = − limm→∞ ∞∑ k=0 1 nk1 tr(F tF )k zk+1 . (6)
The traces themselves are expended out as tr(F tF )k = ∑
Fa1α1Fa1α2 · · ·Fakα1 , (7)
where the sum runs over matrix indices a1, . . . , ak, α1, . . . αk. The essence of the moment method involves analyzing the asymptotic contribution of each term in the sum based on its combinatorial type and the details of F , and resumming the results to obtain m(z).
We refer the reader to (Erdos & Yau, 2017; Tao, 2012) for more details about these methods and additional background on random matrix theory.

heading: 3 RESULT FOR STIELTJES TRANSFORM

heading: 3.1 MAIN THEOREM
text: We make the following assumptions on the data matrix X and bias vector b: 1. 1n0 ∑ aXaαXaβ = δαβσ 2 X +O ( 1/ √ n0 ) for all α and β
2. the empirical eigenvalue distribution of 1n0X tX converges in distribution to a measure µX 3. 1n1 ∑n1 a=1 δba → µB in distribution.
Theorem 1. Define σZ = σWσX and resolvent G(z) = (
1 n1 F tF − zI
)−1 . Then under the above
assumptions and for all z such that =z > 0, the transforms 1 m trG(z) and 1 m tr ( 1 n0 XtXtrG(z) ) , (8)
converge in probability to the unique solution, m(z) and m̃(z), of the Eqn. (9) that map C+ to C+: m(z) = ES∼µX [ 1
C0(z) + SC1(z)
] and m̃(z) = ES∼µX [ S
C0(z) + SC1(z)
] (9)
where,
C0(z) := −z + EB∼µB [ η(B)− ζ(B)
D(B)
] , C1(z) := EB∼µB [ ζ(B)
D(B)
] , (10)
D(B) := 1 + ψ
φ (ζ(B)m̃(z) + (η(B)− ζ(B))m(z)) , (11)
and where η(B) and ζ(B) are the Gaussian expectations
η(B) := EN∼N (0,σZ) [ f(N ;B)2 ] and ζ(B) := EN∼N (0,σZ) [σZf ′(N ;B)] 2 . (12)
The proof is quite involved and is presented in the supplementary material. The basic idea is to derive a multivariate Gaussian random matrix model with the same correlation structure as F , then derive a self-consistent equation (SCE) using the resolvent method for this linearized version of F , F lin.
Remark 1. The self-consistent equations consist of two coupled equations involving the Stieltjes transform m(z) and an auxiliary object m̃(z), (cf. (Paul & Silverstein, 2009, Eq. (2))), which we will see in Cor. 1 essentially measures the autoencoding capacity of the network.
Remark 2. Note that the self-consistent equations contain an expectation over the limiting spectral density of the input data. While the assumptions on the data matrix X in Thm. 1 are quite general, they may not be optimal. See Sec. 3.3, where we show strong agreement with empirical data from MNIST and CIFAR-10 and for a range of synthetic distributions. This suggests that the theorem may hold for even more general data distributions.

heading: 3.2 ALTERNATE REPRESENTATION AND LIMITING RESULTS
text: When the data distribution is iid Gaussian, the expectations in Eqn. (9) can be expressed in closed form, though one must be careful to choose the correct branch of the resulting function. For simplicity and future reference, we focus on the setting where 0 < φ ≤ ψ ≤ 1, in which case we have the
coupled algebraic equations,
m(z) = C1 − (C0 + C1)φ+
√ C21 + 2(C0 − C1)C1φ+ (C0 + C1)2φ2
2C0C1 (13)
m̃(z) = 1− C0m(z)
C1 . (14)
When µB in Eqn. (10) is trivial, i.e. a single Dirac mass, then the result should reduce to the single activation function case with F = f(WX), which was studied in (Pennington & Worah, 2017). Indeed, writing η = EB [η(B)] and ζ = EB [ζ(b)] for such a distribution, and eliminating m̃(z) from Eqns. (13) and (14), we find that m(z) satisfies the following quartic polynomial:
0 = (z2ζ2ψ2)m(z)4 + (2zζ2ψ(ψ − φ)m(z)3 + (ζ2(ψ − φ)2 + zζφψ + zηφ2ψ)m(z)2
+ (ζφ(ψ − φ) + φ2(zφ+ η(ψ − φ))m(z) + φ3 (15)
which agrees with the result in (Pennington & Worah, 2017) upon identifying m(z) = −(1 − φ/ψ)/z − φ/ψG(z).

heading: 3.3 SPECTRAL DENSITY ESTIMATES
text: The self-consistent equations in Thm. 1 can be solved numerically by iterating Eqn. (9) until convergence, using numerical integration. By utilizing the Stieltjes inversion formula, Eqn. (5), we can extract predictions for the limiting eigenvalue density of 1n1F
tF . The results show close agreement with empirical spectral simulations from several interesting practical datasets and synthetic distributions, see Figs. 1 and 2.

heading: 4 RIDGE-REGULARIZED NOISY AUTOENCODER
text: We consider the problem of kernel ridge regression with random features given by F = f(WX;B) and noisy regression targets given by Y = AX + for some A ∈ Rn2×n0 and independent Gaussian noise ∈ Rn2×m such that ij ∼ N (0, σ2 ). As is common in the literature of high-dimensional statistics, we assume an isotropic prior on A such that AtA→ n2n0σ 2 AI as m,n0, n2 →∞. Note that when σA = 0, we have the pure memorization setting studied by Pennington & Worah (2017). Corollary 1. Let W ∗2 be the minimizer for regularized training loss
L(W2) = 1n2m ||Y −W2F || 2 + n1n2mγ||W2|| 2 , (16)
with random features F = f(WX;B). Then the asymptotic training error converges in probability as
Etrain = 1 n2m ||Y −W ∗2 F ||2 → γ2
d
dγ
( σ2Am̃(−γ) + σ2 γ2m(−γ) ) . (17)
In particular, we see that the derivative of the Stieltjes transform m′(−γ) measures the capacity to learn noisy labels, whereas m̃′(−γ) measures pure autoencoding capacity. See Fig. 3 for a comparison between these theoretical predictions and simulation.
Proof. The optimal weights for the regularized loss are given by
W ∗2 = 1 n1 Y G(γ)F t for G(γ) = ( 1n1F tF + γI)−1,
resulting in training error
Etrain = 1 n2m ||Y −W ∗2 F ||2 = γ2 1n2m tr
( Y tY G(γ)2 ) = γ2σ2A 1 n0m tr ( XtXG(γ)2 ) + γ2σ2 1 m tr ( G(γ)2
) = γ2σ2A d dγ (m̃(−γ)) + γ 2σ2 d dγ (m(−γ))
= −γ2 ( σ2Am̃ ′(−γ) + σ2 m′(−γ) ) .
Remark 3. As in (Pennington & Worah, 2017), there is a scaling homogeneity in the Etrain: an increase in the regularization constant γ can be compensated by a decrease in scale of W2, which, in turn, can be compensated by increasing the scale of F , which is equivalent to increasing η(b) and ζ(b). Owing to this homogeneity, we are free to choose a normalization of the activation function for which Eb∼µB [η(b)] = 1.

heading: 4.1 NONLINEAR MIXTURES CAN OUTPERFORM SINGLE NONLINEARITIES
text: The bias term in our random feature model can be viewed as one way of defining a distribution over activation functions. The choice of distribution, in general, affects the performance of the model on a given task – as quantified by the expectations in Eqn. (10). A key benefit of our model and analytical approach is that it permits nontrivial distributions of nonlinear activation functions.
In this section, we build on our results for the training loss on noisy autoencoder tasks to examine the benefits of utilizing nonlinear mixtures. The goal here is not to identify “good" mixtures, since this will clearly be a dataset- and architecture-dependent question, nor is it to identify a large performance gap. Instead, we merely seek to demonstrate a proof-of-principle, namely that there exist non-trivial distributions over nonlinearities that can provably outperform the best possible single nonlinearity.
For this analysis, we consider the simplest possible nontrivial distribution over activation functions: a Bernoulli mixture of two different functions. To each of these functions we associate two constants, η and ζ, which derive from Eqn. (12) but have no B-dependence since each function is a single nonlinearity. Concretely, let
η = EN∼N (0,σZ) [ f(N)2 ] and ζ = EN∼N (0,σZ) [σZf ′(N)] 2 . (18)
For the two functions themselves, we utilize (i) a "pure linear" activation function with η = 1, i.e. the identity function and (ii) a "pure nonlinear" (Hastie et al., 2019) activation function with η = 1 and ζ = 0. (The particular purely nonlinear function in (ii) is irrelevant, as our theory predicts and our experiments confirm; see Fig. 4(a)). To be precise, we define for p ∼ Bernoulli(p),
fp(x) := { x if p = 0 gζ=0(x) if p = 1 , (19)
where gζ=0 is any function with η = 1 and ζ = 0 (see, e.g., the functions in Fig. (3) of (Pennington & Worah, 2017)).
The task of computing Etrain for fp is cumbersome but purely algebraic. To see how to proceed, notice that the expectations in Eqn. (10) are simple for fp:
C0(z) = −z + p
1 + ψ/φm(z) , C1(z) = 1− p 1 + ψ/φ m̃(z) . (20)
Plugging these equations into Eqns. (13) and (14), collecting terms and simplifying yields a set of coupled polynomial equations for m(z) and m̃(z). Taking the total derivative of these equations with respect to z yields two additional equations which can be solved to express m′(z) and m̃′(z) in terms of m(z) and m̃(z). Combining these results produces a polynomial system whose solution2 encodes Etrain through Eqn. (17). Fig. 4(a) shows the result of this calculation in solid lines for various values of γ, while the 1σ error bars show empirical simulations with finite networks. The red stars in the figure show that for many values of γ, the optimal mixture percentage is intermediate, i.e. 0 < p < 1.
2Special care must be taken in selecting the correct root of this equation, in accordance with the condition that m(z) ∼ −1/z for large |z|.
The question is, does a non-trivial mixture actually outperform a single nonlinearity? First, we must understand the performance of the optimal single nonlinearity. We note that owing to the homogeneity of the the training loss in η, ζ , and γ, we can assume without loss of generality that η = 1. Therefore the entire effect of the nonlinearity should be encoded in the single constant ζ. In Fig. 4(a), we plot our theoretical prediction for Etrain in solid lines and empirical simulations for finite networks as 1σ error bars. The activation function used for each simulation is chosen randomly, conditional on the value of ζ. So the good agreement in the left panel of Fig. 4(a) demonstrates not just the correctness of our theoretical result but also the fact that Etrain depends on the activation function solely through the constant ζ. The blue stars in this figure indicate that the optimal single nonlinearity is neither purely linear (ζ = 1), nor purely nonlinear (ζ = 0), but rather something in between.
For this particular problem setup, the performance of the optimal single nonlinearity and the optimal Bernoulli mixture are rather close, as indicated by the top panel of Fig. 4(b). However, owing to our precise analytical formulation, we can evaluate the training loss to high precision and observe that there is indeed a difference in performance between the two models, as shown in the bottom panel of Fig. 4(b). This result establishes that there are some problems for which even the best single nonlinearity is outperformed by a mixture of nonlinearities.

heading: 5 CONCLUSIONS
text: In this work we studied the feature matrix F = f(WX;B) where W is a random matrix with iid Gaussian entries. Under mild assumptions on X and B, we obtained an exact analytic formula, Eqn. (9), that characterizes the Stieltjes transform of the spectral density of F . The result allowed us to describe the exact training loss of a ridge-regularized noisy autoencoder in the high-dimensional, large-dataset limit, providing one of the first closed-form solutions to a non-trivial model in this limit. We found excellent agreement between the asymptotic predictions of Eqn. (9) and a variety of finite-dimensional empirical simulations.
We also advanced the interpretation of the bias B as one particular way of parameterizing a distribution of activation functions. Indeed, our derivations proceed completely unchanged whether this distribution is of the traditional additive form f(·+B) or the more general f(·;B). By examining the latter, we showed that there are configurations in which a non-trivial distribution over activation functions provably outperforms the best possible single activation function. This opens the door to future investigations regarding optimal methods for parameterizing distributions over activation functions for approximate kernel methods, and suggests the possibility that mixtures of nonlinearities could be a useful design consideration when constructing neural network architectures.

heading: A Random Matrix Perspective on Mixtures of Nonlinearities in High Dimensions
text: The main goal of the supplementary material is to prove Thm. 1. This is achieved in three steps. The first step is to note that the limiting Stieltjes for the kernel matrix F tF is equivalent to a multivariate a multivariate Gaussian random matrix model with an identical correlation structure. This was proven in great generality in Banna et al. (2015) (see Theorem 5.). Similar results for covariance matrices have been previously studied Bai & Zhou (2008), but since F does not have identically distributed rows (or columns) they are inadequate for our case. Specifically, Theorem 5 from Banna et al. (2015), implies that the two Stieltjes transforms have the same limit for all z such that =z > 0 almost surely. The second step is to derive the leading order correlation structure (or kernel) for F (Lem. 1), and specify a multivariate Gaussian model with the same structure (Subsec. A.1). We refer to this Gaussian model as a linearized model, since it removes the nonlinearity f from the random part of our random matrix ensemble. Without loss of generality, we may assume that EFaα = 0 for all a and α. This is because EFaα only depends on a (as we demonstrate in Sec. A), so we can center F using a rank-1 matrix—a transform that has no effect on the limiting Stieltjes transform as it represents an O(1/m) change in the spectrum (see for more detail Bai & Zhou (2008)). The third a final step of the proof is a derivation of a self-consistent equation specifying the limiting Stieltjes transform of the Gaussian model. We do this using the resolvent method (see Erdos & Yau (2017) for an introduction). While the application of this method to multivariate Gaussian covariance matrices is standard, we note the form of the SCE is new and that we include the derivation here for completeness.
At the end of the supplementary material, we include some additional figures.

heading: A CORRELATION STRUCTURE OF F
text: To begin we derive an asymptotic form for correlation structure of F . For simplicity of presentation, we derive the results for f(·, B) ≡ f(·+B), but the argument generalizes directly. We recall that W is an n1 × n0 random matrix with mean 0 and variance σ2w/n0 iid Gaussian entries. Define cαβ := ∑ kXkαXkβ/n0. The data matrix X and bias vector B are deterministic, so all expectations are over W (and can be thought of as conditional on X and B). Moreover, X and B satisfy the following assumptions:
• cαβ = δαβσ2W +O ( 1/ √ n0 ) for all α and β;
• The empirical spectral distribution of 1n0X tX converges in distribution to some measure
µX ; • 1n1 ∑n1 a=1 δba → µB in distribution, or equivalently 1 n1 ∑n1 a=1 g(ba) → EB∼µBg(B) for
bounded, continuous functions g.
Observe that
Zaα := n0∑ k=1 WakXkα and Zbβ := n0∑ k=1 WbkXkβ (21)
are jointly Gaussian. Moreover,
EZaα = 0 and EZaαZbβ = σ2W cαβ1 (a = b) . (22) Note in particular, that Zaα and Zbβ are independent if a 6= b. For convenience, we normalize cαβ and define c̃αβ := cαβ/σ2X and εα := c̃αα − 1. In order to specify the correlations, we employ transforms of the activation function f . The transforms are Gaussian integrals of f at different locations and scales:
ξk(x) := EσkXσkW f (k) (σXσWN + x) , (23)
ηi,j(x) := Eσi+jX σ i+j W f (i)(σXσWN + x)f (j)(σXσWN + x), (24)
and γi,j(x) := Eσi+jX σ i+j W f (i)(σXσWN + x)f (j)(σXσWN ′ + x) = ξi(x)ξj(x). (25)
Remark 4. For simplicity, while writing f (i) we assume differentiability of f up to ith order. However, this assumption can be relaxed, as since the derivatives are only used to take a Gaussian expectation, they can be expressed using Stein’s lemma. For example, ξ1(x) = ENf(σXσWN+x). Thus, similar results should hold for rectified linear units for example, regardless of their lack of differentiability at 0. Lemma 1. To leading order the correlation structure of F is
EFaαFbβ =  0 if a 6= b η00(ba) if α = β and a = b γ11(ba)cαβ if α 6= β and a = b . (26)
Proof of Lem. 1. Before we prove Lem. 1, we state a simple lemma that is easy to show using Taylor expansion.
Lemma 2. Suppose ε is a small constant and that N ∼ N (0, 1), then
ENf ( σXσW √ 1 + εN + b ) = ξ0(b) + ξ2(b)
2 ε+O(ε2), (27)
ENf ( σXσW √ 1 + εN + b )2 = η00(b) + (η11(b) + η02(b)) ε+O(ε2), (28)
and EN ( σXσW √ 1 + ε ) f ′ ( σXσW √ 1 + εN + b ) = ξ1(b) + ξ1(b) + ξ3(b)
2 ε+O(ε2). (29)
Recall ξ0(x) = 0 by assumption. Using the observation in Eqn. (21), we repeatedly rewrite the expectations in F tF as we can reduce the expectation over W to an expectation over at most two correlated Gaussian random variables. Let N and N ′ denote independent, standard Gaussian random variables. Using Lem. 2, we see
Ef(Zaα) = Ef(σW √ cααN + ba) = ξ2(ba)
2 εα +O(ε2α). (30)
Again, using Lem. 2, we get
Ef(Zaα)2 = η00(ba) + (η11(ba) + η02(ba)) εα +O(ε2α). (31) For the covariance calculation, we define ραβ := cαβ/ √ cαα and recall that ραβ is order 1/ √ m. Thus,
Ef(Zaα)f(Zaβ) = Ef(σW √ cααN ′ + ba)f ( σW ραβN ′ + σW √ cββ − ρ2αβN + ba ) (32)
= (EN ′f(σW √ cααN ′ + ba)) ( ENf ( σW √ cββ − ρ2αβN + ba )) (33)
+ (EN ′σW ραβN ′f(σW √ cααN ′ + ba)) ( ENf ′ ( σW √ cββ − ρ2αβN + ba )) .
(34)
For the first term, we get
ξ0(ba) 2 + ξ0(ba)ξ2(ba) εα + εβ 2 +O(m−1) = O(m−1). (35)
For the second term, we have
c̃αβ
( 1− εα + εβ
2 +O(m−1)
)( ξ1(ba) + ξ1(ba) + ξ3(ba)
2 εα +O(m−1)
) (36)
× ( ξ1(ba) + ξ1(ba) + ξ3(ba)
2 εβ +O(m−1)
) (37)
= γ1,1(ba)c̃αβ +O(m−1). (38) Thus, to second order we have
Ef(Zaα)f(Zaβ) = γ1,1(ba)c̃αβ +O(n−10 ). (39)
A.1 LINEARIZED MODEL
Define F lin := CΘ1Σ + ( V2 − C2 )1/2 Θ2, (40)
where: Θ1 and Θ2 are n1 ×m matrices that have iid Gaussian entries with mean 0 and variance 1/n1; 2) V and C are n1 × n1 diagonal matrices with entries
V2a = η00(ba) ≡ η(ba) and C2a = γ11(ba) ≡ ζ(ba), (41) for a bias vector b; 3) Σ is a matrix square root of 1n0X tX .
A simple calculation shows that the entries of the matrix F lin match the first and second mixed moments of F given in Lem. 1. However, unlike F , F lin is linear in the random matrices Θ1 and Θ2; in this sense, it is a linearization of F .

heading: B DERIVATION OF SELF-CONSISTENT EQUATION FOR F LIN
text: Theorem 2. Let mm(z) be the Stieltjes transform of 1n1F lintF lin, i.e. 1m trG(z). Then with probability 1, as n1 → ∞, for all z such that =z > 0, mn1(z) → m(z), where m(z) is the solution of the coupled equations
m(z) = ES∼MP (φ) [
1
C0(z) + SC1(z)
] and m̃(z) = ES∼MP (φ) [ S
C0(z) + SC1(z)
] (42)
with
C0(z) := −z + EB∼µB [ η(B)− ζ(B)
D(B)
] , C1(z) := EB∼µB [ ζ(B)
D(B)
] , (43)
D(B) := 1 + ψ
φ (ζ(B)m̃(z) + (η(B)− ζ(B))m(z)) . (44)
We want to study the eigenvalues of 1n1F lintF lin. Note without loss of generality it is sufficient to consdier diagonal Σ, since we can diagonalize Σ using some orthogonal matrices O and O′ as it is the square root of a positive definite matrix. Moreover, these orthogonal matrices, when applied to either Θ1 or Θ2 do not change their distributions. Thus, diagonalizing Σ, so that Σαα = √ λXα , where λXα are the eigenvalues of X tX/n0, results in an equivalent matrix ensemble in distribution (see Silverstein & Bai (1995b) for more detail). With this simplification, F lin has independent entries given by
F linaα = CaΘ1aα √ λXα + √ V2a − C2aΘ2aα. (45)
To make the derivation easier, we can partly linearize the problem by studying the matrix
H :=
[ −zI F lint
F lin −I
] . (46)
By the Schur complement formula, one easily finds
mm(z) := 1
m Tr(F lin t F lin/n1 − zI) =
1
m m∑ α=1 Gαα(z) (47)
and zm̃m(z) := 1
n1 Tr(F linF lin t /n1 − zI) =
1
n1 n1+m∑ a=n1+1 Gaa(z), (48)
where G is the inverse of H . Again by the Schur complement formula
1
Gαα = −z − n1∑ a,b=1 F linaαF lin bαG (α) m+a,m+b (49)
and 1
Gaa = −1− m∑ α,β=1 F linaαF lin aβG (a) αβ (50)
for α ∈ {1, . . . ,m} and a ∈ {m+ 1, . . . ,m+ n1} and G(a) is the inverse of the minor H(a). Since G(α) is independent of Θ11α, . . . ,Θ 1 n1α and Θ 2 1α, . . . ,Θ 2 n1α, we see by taking the expectation over these variables that
E n1∑
a,b=1
F linaαF lin bαG (α) m+a,m+b =
1
n1 n1∑ a=1 ( λXα ζ(ba) + η(ba)− ζ(ba) ) G (α) m+a,m+a. (51)
Moreover, standard concentration inequalities and the Ward identity (see Erdos & Yau (2017)) show∣∣∣∣∣∣ n1∑
a,b=1
F linaαF lin bαG (α) m+a,m+b − E n1∑ a,b=1 F linaαF lin bαG (α) m+a,m+b ∣∣∣∣∣∣ (52) ≤ C  1 m2 ∑ a,b ∣∣∣G(α)m+a,m+b∣∣∣2 1/2 ≤ C ( 1 m2 ∑ a =G(α)m+a,m+a =z )1/2 ≤ O(1/ √ m) (53)
with high probability. Similar bounds are easily obtained for ∑m α,β=1 F lin aαF lin aβG (a) αβ .
We may also replace G(α) with G at the expense of another small error that can be bounded using the Cauchy interlacing theorem: ∣∣∣G(α)m+a,m+a −Gm+a,m+a∣∣∣ ≤ O(1/m). Using this control over these sums, we see
1
Gαα = −z − n1∑ a,b=1 F linaαF lin bαG (α) ab (54)
= −z − 1 n1 n1∑ a=1 ( λXα ζ(ba) + η(ba)− ζ(ba) ) Gm+a,m+a +O(1/ √ m) (55)
and
1
Gaa = −1− m∑ α,β=1 F linaαF lin aβG (a) αβ (56)
= −1− 1 n1 m∑ α=1 ( λXα ζ(ba) + η(ba)− ζ(ba) ) Gαα +O(1/ √ m). (57)
Finally, we invert Eqn. (56), multiply by λXα ζ(ba) + η(ba)− ζ(ba), and average over a to find
1
n1 n1∑ a=1 ( λXα ζ(ba) + η(ba)− ζ(ba) ) G (α) m+a,m+a (58)
= − 1 n1 n1∑ a=1 λXα ζ(ba) + η(ba)− ζ(ba) 1 + ψφ (ζ(ba)m̃(z) + (η(ba)− ζ(ba))mm(z)) +O(1/ √ m)
(59)
= − 1 n1 n1∑ a=1 λXα ζ(ba) + η(ba)− ζ(ba) 1 + ψφ (ζ(ba)m̃(z) + (η(ba)− ζ(ba)mm(z)) +O(1/ √ m) (60)
= −EB∼N (0,σ2b )
[ λXα ζ(B) + η(B)− ζ(B)
1 + ψφ (ζ(B)m̃(z) + (η(B)− ζ(B))mm(z))
] + o(1), (61)
where m̃m(z) = 1m ∑ α λ X α Gαα, where we Taylor expanded in the second step, and where we used our assumption on B.
We can now invert Eqn. (54) and average over α to find
mm(z)→ ES∼µX  1 −z + EB∼N (0,σ2b ) [ Sζ(B)+η(B)−ζ(B)
1+ψφ (ζ(B)mσ(z)+(η(B)−ζ(B))m(z))
]  . (62)
Similarly,
m̃m(z)→ ES∼µX  S −z + EB∼N (0,σ2b ) [ Sζ(B)+η(B)−ζ(B)
1+ψφ (ζ(B)m̃(z)+(η(B)−ζ(B)m(z))
]  . (63)
Note that the integral over S here is an integral over the limiting distribution of the data µX , which we assume
1
m ∑ α δλXα → µX (64)
in distribution. In the case of iid Gaussian data, this is exactly given by the Marchenko-Pastur distribution.
C ADDITIONAL FIGURES

Review summaries:
This paper investigates the asymptotic spectral density of a random feature model F(Wx + B).  This is an extension of existing result that analyzed a model without the bias term, i.e., F(WX). This extension requires a modification of the proof technique. In addition to that, it analyzed a mixture of linear and non-linear activation functions, and show that mixture is better than single nonlinearity in terms of expected training error for ridge regression estimators.