Analysis of one-hidden-layer Neural Networks via the Resolvent Method
Year: 2022

Abstract:
In this work, we investigate the asymptotic spectral density of the random feature matrix M = Y Y ∗ with Y = f(WX) generated by a single-hidden-layer neural network, where W and X are random rectangular matrices with i.i.d. centred entries and f is a non-linear smooth function which is applied entry-wise. We prove that the Stieltjes transform of the limiting spectral distribution approximately satisfies a quartic self-consistent equation, which is exactly the equation obtained by Pennington and Worah [22] and Benigni and Péché [6] with the moment method. We extend the previous results to the case of additive bias Y = f(WX +B) with B being an independent rank-one Gaussian random matrix, closer modelling the neural network infrastructures encountered in practice. Our key finding is that in the case of additive bias it is impossible to choose an activation function preserving the layer-to-layer singular value distribution, in sharp contrast to the bias-free case where a simple integral constraint is sufficient to achieve isospectrality. To obtain the asymptotics for the empirical spectral density we follow the resolvent method from random matrix theory via the cumulant expansion. We find that this approach is more robust and less combinatorial than the moment method and expect that it will apply also for models where the combinatorics of the former become intractable. The resolvent method has been widely employed, but compared to previous works, it is applied here to non-linear random matrices.

heading: 1 Introduction
text: Machine learning has seen many successful achievements in recent years. Applications in face identification, object and speech recognition, translation, email spam filtering, navigation, medical diagnosis, etc. have proved the enormous potential of machine learning for day-to-day live [16, 11]. Deep neural networks have turned out to be a particularly powerful machine learning method, and understanding the theoretical underpinning of their success has received tremendous attention in mathematics, physics and computer science.
A fully-connected, feed-forward neural network with L hidden layers of dimensions n1, . . . , nL can be modelled as follows:
fθ(x) = β ∗f(W (L) f(W (L−1) f(. . . f(W (1)x) . . . ))) ∈ Rd,
where x ∈ Rn0 denotes the input data vector and f : R → R is a non-linear activation function which is applied entry-wise. We denote the parameters of the network by θ := (W (1), . . . ,W (L), β), where W (l) ∈ Rnl×nl−1 for 1 ≤ l ≤ L and β ∈ RnL×d are the matrices of the weights. In the classical setting of supervised learning, we are given a training set of (say) m samples of input feature vectors xi ∈ Rn0 with associated target vectors zi ∈ Rd. For example, xi may encode the pixels of a photograph of an animal and the target zi may label the species of the
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
animal in the image. Roughly speaking, the goal of supervised learning is to learn the mapping between the feature and the target vectors based on a given training set in order to predict the output of new unlabelled data. Let X = (x1 . . . xm) ∈ Rn0×m be the matrix of the data and let Z = (z1 . . . zm) ∈ Rd×m be the target matrix. Then, the aim of the network is to find optimal parameters θ such that fθ(X) = (fθ(x1), . . . , fθ(xm)) ∈ Rd×m approximates the target Z optimally. During the training phase, weights are adjusted in order to minimize the empirical risk R(θ) = EL(fθ(X), Z), where L(·, ·) is a given loss function, usually involving some penalty for large weights θ in order to avoid over-fitting. Stochastic gradient descent (SGD) and its variants with back-propagation are the most commonly used algorithms for training multilayer networks by iteratively updating the parameters into the direction of the negative of the gradient of the empirical risk. For a much more complete survey, we refer the reader to [11].
In the present paper, we will focus on a single-hidden-layer neural network of the form fθ(X) = β∗Y with Y = f(WX). This model was first studied by Louart, Liao, and Couillet [19] for the case where the data matrix X is deterministic and W is a matrix of random weights (in particular, the weights are given by functions of standard Gaussian random variables), and by Pennington and Worah [22] for the case where X and W are independent random matrices with both centred Gaussian entries. In both papers, the matrix β ∈ Rn1×d is the only parameter to be learned and is chosen as the unique minimizer of the ridge-regularized least squares loss function
L(fθ(X), Z) = 1
2dm ∥Z − β∗Y ∥2F + γ∥β∥2F ,
where γ > 0 is the learning rate. The unique minimizing weight matrix β̂ is then equal to β̂ = Y G(−γ)Z∗, where
G(−γ) = ( 1
m Y ∗Y + γ )−1 is the resolvent of 1mY
∗Y . As proved in [19, 22], the expected training loss Etrain is related to −γ ∂∂γG(−γ), and thus also to the Stieltjes transform of the limiting spectral measure of 1 mY
∗Y . Here the Stieltjes transformmµ of a probability measure µ on R is defined asmµ(z) := ∫ R(x−z)
−1 dµ(x) for z ∈ C such that ℑz ≥ 0, and for µn1 being the empirical probability measure of the n1 eigenvalues of 1mY ∗Y is related to the resolvent via mµn1 (z) = 1 n1
TrG(z). The performance of one-hiddenlayer neural networks depends on the asymptotic spectral properties of the matrix 1mY
∗Y . Pennington and Worah [22] investigated the limiting spectral measure of the random matrix M = 1mY Y
∗ and derived the quartic self-consistent equation
1 + zg∞ = θ1(f)g∞
( 1− ϕ
ψ (1 + zg∞)
) − θ2(f)
ψ g∞(1 + zg∞)
( 1− ϕ
ψ (1 + zg∞) ) + θ2(f)(θ1(f)− θ2(f))
ψ g2∞
( 1− ϕ
ψ (1 + zg∞)
)2 ,
(1)
where g∞(z) := limn1→∞ g(z) and g(z) := 1 n1 TrG(z) is the Stieltjes transform, which is approximately satisfied, g∞(z) ≈ g(z), by g(z) in case of Gaussian W,X . It is notable that the asymptotic spectrum of Y ∗Y for large dimensions such that n0/m → ϕ ∈ (0,∞) and n0/n1 → ψ ∈ (0,∞) depends on the non-linear function f only through two integral parameters θ1(f) and θ2(f), where θ1(f) is the Gaussian mean of f2 and θ2(f) is the square of the Gaussian mean of f ′ (c.f. (5)). Benigni and Péché [6] then extended this model to random matrices W and X with general i.i.d. centred entries, and obtained the same self-consistent equation (1). We mention that (1) may be reduced, for some special cases, to the quadratic equation that is satisfied by the Stieltjes transform mµMP of the Marchenko-Pastur distribution µMP [20]. This means that for some activation functions, the non-linear random matrix model has the same limiting spectral distribution as that of sample covariance matrices XX∗ (in other cases, the equation can simplify to the cubic equation approximately satisfied by product Wishart matrices [7, 10]). This can be generalised to multilayer networks: [22] found experimentally that the singular value distribution is preserved through multiple layers by activation functions with θ2(f) = 0 and is given by the Marchenko-Pastur distribution in each layer. This conjecture was proved in [6] for the general case of bounded activation functions. Moreover, [19] performed a spectral analysis on the Gram matrix model with general training data and proved that, in the large dimensional regime, the resolvent of Y ∗Y has a similar behaviour as that
observed in sample covariance matrix models. This was extended in [17] by considering Gaussian mixture of data. We also refer to the recent paper [18]. In the context of multilayer feedforward neural networks, Fan and Wang [9] analysed the eigenvalue distribution of the Gram matrix model, where the weights are at random and the input vectors are assumed to be approximately pairwise orthogonal. In particular, they showed that the limiting spectral distribution converges to a deterministic limit and, at each intermediate layer, this limit corresponds to the Marchenko-Pastur map of a linear transformation of that of the previous layer.
In recent years, there has been some progress in the asymptotic analysis of the eigenvalue distribution of another Gram matrix, the so-called Neural Tangent Kernel (NTK). Consider a multilayer neural network and denote by J = ∇θfθ(X) the Jacobian matrix of the network outputs with respect to the weights θ. Then, the NTK is the Gram matrix of J , defined by KNTK = J∗J . It was shown in [14] that the NTK at random initialization converges, in the infinite-width limit, to a deterministic kernel and it remains constant during the whole training time of the network. Subsequently, [23] analysed the spectrum of the sample covariance matrix JJ∗ in a single-hidden-layer neural network, and provided an exact asymptotic characterization of the spectral distribution of JJ∗ with random Gaussian weights and data. Recently, [9] proved that the limiting spectral measure of the NTK converges to a deterministic measure, which may be described by recursive fixed-point equations that extend the Marchenko-Pastur distribution.
The present paper is structured as follows. In the first part we consider the non-linear random matrix model studied in [6] and we compute its asymptotic spectral density. We follow the resolvent method via the cumulant expansion which, together with the moment method, is a standard approach to obtain the asymptotics for the empirical spectral density. In particular, we compute the self-consistent equation that is approximately satisfied by the Stieltjes transform of the limiting spectral distribution. This is a quartic equation and is the same as that found in [6]. In [6, 22] the authors relied on the method of moments: they approximated general non-linear functions by polynomial ones and then computed the asymptotics of high moments ETr(Y k) with Y = f(WX) to obtain the limiting measure via its moments. Conversely, we approach matters in a more robust and less combinatorial fashion by applying the resolvent method: we consider Y as a random matrix with correlated entries and then we directly derive a self-consistent equation for its resolvent. In particular, we prove that the random matrix Y has cycle correlations, in the sense that the joint cumulant does not vanish when the random variables Yij’s are joined by a cycle graph. We find that the variance of Yij is given by the parameter θ1(f), whereas for k > 1 the cumulants κ(Yi1i2 , Y ∗ i2i3 , Yi3i4 , . . . , Y ∗ i2ki1
) are powers of θ2(f). We note that in the random matrix literature matrices with general decaying correlations have been studied previously, see e.g. [3, 8, 1]. However, the cycle correlations of Y considered in the present paper are much stronger compared to these previous results. The second part of this paper concerns the additive bias case which is a more realistic model for machine learning applications. More precisely, we consider the random feature matrix Y = f(WX +B), where B is a rectangular rank-one Gaussian random matrix, and derive a characterization of the Stieltjes transform of the limiting spectral density. We chose B to be rank-one since for the most commonly used neural network architectures the added bias is chosen equal for each sample. [2] studied the bias case for deterministic data and i.i.d. Gaussian random weights, and computed the exact training error of a ridge-regularized noisy autoenconder in the high-dimensional regime. Interestingly we find that in the case of additive bias it is impossible to choose an activation function f such that the eigenvalue distribution is preserved throughout multiple layers, unlike in the bias-free case where θ2(f) = 0 yields the Marchenko-Pastur distribution in each layer. Finally, we remark that in the bias-free case our proof via the resolvent method has no significant advantage compared to the moment method, beyond requiring less combinatorics. The main advantage of the resolvent approach is that it allows to include an additive bias without much additional effort.

heading: 2 Model and main results
text: We consider a random data matrix X ∈ Rn0×m with i.i.d. random variables Xij with distribution ν1 and a random weight matrix W ∈ Rn1×n0 with i.i.d. weights Wij with distribution ν2. We assume that both distributions are centred with variance EX2ij = σ 2 x and EW 2 ij = σ 2 w. Moreover, we assume
that the distributions ν1, ν2 have finite moments of all orders1. Since for 1 ≤ i ≤ n1 and 1 ≤ j ≤ m we have (
WX √ n0 ) ij = 1 √ n0 n0∑ k=1 WikXkj ,
we note that in light of the central limit theorem the entries of the matrix WX√n0 are approximately N (0, σ2wσ2x)-normally distributed random variables. Therefore, for any t > 0, we have the large deviation estimate
P ( max i,j ∣∣∣∣ (WX)ij√n0 ∣∣∣∣ > t) ≲ n20 e−t2/2σ2wσ2x ,
where we use the notation A ≲ B as shorthand for the inequality A ≤ cB for some constant c. Let f : R → R be a C∞ function with zero mean with respect to the Gaussian density of standard deviation σwσx, i.e. ∫
R f(σwσxx)
e−x 2/2
√ 2π dx = 0. (2)
We consider the random feature model generated by a single-hidden-layer neural network,
M = 1
m Y Y ∗ ∈ Rn1×n1 with Y = f ( WX √ n0 ) , (3)
where the activation function f is applied entry-wise. Let χ : R → R be a smooth cut-off function that is equal to one for |x| ≤ 1 and zero for |x| ≥ 2. We then replace f by f(·)χ(log−1(n0) ·). In particular, we now have that f is smooth with compact support. Moreover, for any l > 0 and n0 large enough, with probability 1− n−l0 , the singular values of Y remain the same. We are interested in the eigenvalue density of the random matrix M in the infinite size limit. So, we assume that the dimensions of both the columns and the rows of each matrix are large and grow at the same speed, i.e. we introduce some positive constants ϕ and ψ such that
n0 m −→ ϕ and n0 n1 −→ ψ as n0, n1,m→ ∞. (4)
We denote by (λ1, . . . , λn1) the eigenvalues of M and define its empirical spectral distribution by µn1 = 1 n1 ∑n1 i=1 δλi . Then, as n1 grows large, the empirical distribution of eigenvalues converges in distribution to some deterministic limiting density.
Theorem 2.1. There exists a deterministic measure µ = µϕ,ψ(θ1, θ2) such that almost surely weakly
µn1 −→ µ as n1 → ∞.
We notice that if m < n1, then rank(M) = min(n1,m) = m and M has n1 −m zero eigenvalues. In this case, since ϕ/ψ > 1, there exists an atom at 0 with mass µn1(0) = 1−ψ/ϕ > 0, and we have
µn1 = n1 −m n1 δ0 + 1 n1 n1∑ i=1 δλi .
Conversely, if n1 < m, the matrix M has full rank and it is invertible. Since the nonzero eigenvalues of Y Y ∗ and Y ∗Y are the same, the limiting measure µ of Theorem 2.1 turns out to be
µ = ( 1− ψ
ϕ ) + δ0 + µ̃,
where (·)+ = max(0, ·), and µ̃ is the limiting spectral measure of 1mY ∗Y .
We will prove that the deterministic measure µ of Theorem 2.1 is characterized through a quartic self-consistent equation for the Stieltjes transform g(z) = 1n1 TrG(z) of the empirical spectral measure µn1 , where
G(z) = (M − z)−1 ∈ Cn1×n1
1This assumption can be relaxed by a customary cut-off argument, but we refrain from doing so for simplicity.
is the resolvent of the random matrix M and the spectral parameter z lies in the upper half plane H = {z ∈ C | ℑz ≥ 0}. We set
θ1(f) := ∫ R f2(σwσxx) e−x 2/2 √ 2π dx and θ2(f) := ( σwσx ∫ R f ′(σwσxx) e−x 2/2 √ 2π dx )2 . (5)
Then, the following theorem characterizes g as the solution to a quartic equation which depends only on the two parameters θ1(f) and θ2(f).
Theorem 2.2. For some δ, ϵ > 0 and any z ∈ H with ℑz > n− 1 4+ϵ
1 , the measure µ is characterized through the following self-consistent equation∣∣∣∣∣1 + zg − ( θ1 − θ2 ψ (1 + zg) ) g ( 1− ϕ ψ (1 + zg) ) − θ2(θ1 − θ2) ψ g2 ( 1− ϕ ψ (1 + zg) )2∣∣∣∣∣ ≤ n−δ1 (6)
almost surely.
Remark 2.3.
(i) We obtain an analogous result for complex feature and weight matrices W,X , c.f. Remark E.1.
(ii) Note that the quartic self-consistent equation (6) may not have a unique solution such that ℑg(z) > 0. However, it has a unique solution which is analytic in the upper half-plane and satisfies g(z) ∼ −1/z for large |z|.
(iii) Since the resolvent itself satisfies TrG(z)/n1 ∼ −1/z for large |z| and is analytic in the upper half-plane, by continuity Theorem 2.2 implies that g(z) is approximately given by a properly chosen solution of (6). Then, the limiting spectral measure itself can be recovered via the Stieltjes inversion formula,
µ(λ) = lim ϵ→0+
1 π ℑg(λ+ iϵ),
and Theorem 2.1 follows from Theorem 2.2.
(iv) It follows from the self-consistent equation (6) that the limiting spectral measure µ̃ is absolutely continuous w.r.t. the Lebesgue measure, and therefore so is µ away from zero. Moreover, for large z, equation (6) has real solutions and thus via Stieltjes inversion the limiting measure µ is compactly supported.
Remark 2.4. It should be noted that Theorem 2.1 and Theorem 2.2 were proven in [22, 6] under different assumptions and with a different method. The result in [22] was obtained for i.i.d. Gaussian features and weights, whereas [6] extends the result to the case where both the inputs and the random weights have sub-Gaussian tails but are not necessarily Gaussian.
Observing equation (6), we note that if θ2(f) = 0, then the limiting measure µ is exactly the Marchenko-Pastur µMP distribution with parameter ϕ/ψ. Indeed, in this case, g(z) approximately satisfies the quadratic equation
1 + ( z + θ1(f) ( ϕ ψ − 1 )) g(z) + θ1(f) ϕ ψ zg(z)2 ≈ 0, (7)
which corresponds to the self-consistent equation satisfied by the Stieltjes transform of µMP [20]. As discussed in the introduction, this consideration is relevant when studying multilayer networks. Pennington and Worah [22] conjectured that the asymptotic spectral distribution is preserved through multiple layers only by activation functions with θ2(f) = 0 and is given by the Marchenko-Pastur distribution in each layer. Benigni and Péché [6] then proved this conjecture for bounded activation functions satisfying θ2(f) = 0. Moreover, if θ1(f) = θ2(f), then equation (6) becomes cubic. In particular, the equality θ1(f) = θ2(f) holds if and only if f is a linear function (for more details, we refer to the supplementary material in [22]). In this case, M = 1mY Y
∗ with Y =WX , and thus the limiting measure µ corresponds to the limiting spectral distribution of a product Wishart matrix. The spectral density for matrices of this type has been computed in [7, 10].

heading: 2.1 Additive bias case
text: The previous model can be generalised by adding random biases. In neural networks, the bias is an additional parameter that allows the model to better fit the given data. In this case, for each input data x ∈ Rn0 , a bias vector b ∈ Rn1 is added to the vector Wx ∈ Rn1 . We then apply a non-linear function f : R → R in an element-wise fashion to its vector arguments Wx+ b in order to obtain n1 random features f(Wx+ b) ∈ Rn1 . We consider a random bias matrix B ∈ Rn1×m of i.i.d. Gaussian random variables Bij = Bi with zero mean and variance EB2i = σ 2 b . Note that the random matrix B has rank 1. Let X ∈ Rn0×m and W ∈ Rn1×n0 be random matrices with i.i.d. entries, defined as before. Moreover, let f : R → R be a C∞ function satisfying
∫ R f (√ σ2wσ 2 x + σ 2 b x ) e−x 2/2 √ 2π dx = 0. (8)
Just as before, without loss of generality, upon replacing f by f(·)χ(log−1(n0) ·), we may assume that f is a smooth function with compact support. We then define the random matrix M by
M = 1
m Y Y ∗ ∈ Rn1×n1 with Y = f ( WX √ n0 +B ) , (9)
where f is applied entry-wise. We introduce the parameter
σ̃ :=
√ σ2wσ 2 x(σ 2 wσ 2 x + 2σ 2 b )
σ2wσ 2 x + σ 2 b
,
and we define the following integral parameters:
θ1(f) := ∫ R f2 (√ σ2wσ 2 x + σ 2 bx ) e−x 2/2 √ 2π dx,
θ1,b(f) := 1 2πσ̃ √ σ2wσ 2 x + σ 2 b ∫ R2 f(x1)f(x2) exp ( −x 2 1 + x 2 2 2σ̃2 + σ2bx1x2 σ̃2(σ2wσ 2 x + σ 2 b ) ) dx,
θ2(f) := σwσx 2πσ̃ √ σ2wσ 2 x + σ 2 b ∫ R2 f ′(x1)f ′(x2) exp ( −x 2 1 + x 2 2 2σ̃2 + σ2bx1x2 σ̃2(σ2wσ 2 x + σ 2 b ) ) dx.
(10)
We can now state the analogue of Theorem 2.2 in the additive bias case. In particular, the following theorem shows that the normalized trace of the resolvent of M approximately satisfies the selfconsistent equation (6) with parameters given by (10). Theorem 2.5. The Stieltjes transform g satisfies (6) with parameters given by (10), where θ1(f) is replaced by θ1(f)− θ1,b(f). Moreover, there exists a single outlier eigenvalue λmax = n1θ1,b(1 + O(n−1/21 )) of M that is separated from the support of the rest of the spectrum.
We remark that the parameters θ1,b(f), θ2(f) can be alternatively expressed as infinite series, directly demonstrating that for σb ̸= 0 and non-trivial f both coefficients are strictly positive, θ1,b(f), θ2(f) > 0. For notational implicitly, we introduce the Hermite inner product
⟨f, g⟩He := 1√ 2π ∫ R f(x)g(x)e−x 2/2 dx.
Remark 2.6. We have
θ1,b(f) = σ̃√
σ2wσ 2 x + σ 2 b ∑ k≥0 1 k! ( σ2b σ2wσ 2 x + σ 2 b )k ⟨xk, f(σ̃·)⟩2He
θ2(f) = σ2wσ 2 xσ̃√
σ2wσ 2 x + σ 2 b ∑ k≥0 1 k! ( σ2b σ2wσ 2 x + σ 2 b )k ⟨xk, f ′(σ̃·)⟩2He
(11)
and therefore θ1,b(f) = 0, σb ̸= 0 implies that f(σ̃·) is orthogonal to Hermite polynomials of any order, and consequently f ≡ 0. Similarly, θ2(f) = 0, σb ̸= 0 implies that f ≡ const.

heading: 2.2 Multiple layers
text: In [22] it was observed empirically that in the bias-free case activation functions with θ2(f) = 0 have the remarkable property that for multiple layers
Y (l+1) := f(W (l)Y (l)), Y (0) := X (12)
the singular value distributions of Y (1), Y (2), . . . all asymptotically agree (up to scaling) with the probability distribution µ(θ1, θ2) = µ(θ1, 0) from Theorem 2.1. This observation is very natural from our point of view since we find that Y (1) is approximately an i.i.d. random matrix if θ2(f) = 0, c.f. Proposition 3.2 below.
An interesting corollary of our Theorem 2.5 is that a similar isospectral property cannot be ensured for the case of additive bias
Y (l+1) := f(W (l)Y (l) +B(l)), Y (0) := X. (13)
Indeed, in light of Remark 2.6, for σb ̸= 0 we have θ1,b(f), θ2(f) > 0 for all activation functions f , and therefore already the random matrix Y (1) necessarily has leading order correlations, c.f. Proposition 3.3 below. Hence, convergence of the spectral density to the solution of (6) is not expected beyond the first layer. In Fig. 2 we test this result experimentally and choose the activation function f(x) = c1|x| − c2 with c1, c2 such that (2) is satisfied and θ1(f) = 1. We find that in the bias-free case (left), irrespective of the network depth, the eigenvalues of the covariance matrix Y (l)(Y (l))∗ converge to their theoretical limit from Theorem 2.1, exactly as in [22, Fig. 1]2. In the case of an additive bias (right), no such convergence is observed, and this provides empirical evidence of our result.
2In the notation of [22], f = f1.
The spectrum of the covariance matrix Y (l)(Y (l))∗ reflects the distortion of input data through the network and highly skewed distributions indicate poor conditioning which may impede learning performance [22]. Batch normalization seeks to remedy the distortion by normalising by the trace of the covariance matrix Y (l)(Y (l))t in each layer. In [22] it was suggested that choosing activation functions with θ2(f) = 0, i.e. functions which naturally preserve the singular value distribution, may serve as an alternative method of tuning networks for fast optimisation. Our result indicates that in the case of additive bias this alternative is not present. However, batch normalization seems to help stabilising the singular value distribution also in the additive bias case, c.f. Fig. 3.

heading: 3 Outline of proof of Theorems 2.2 and 2.5
text: The proof of both Theorem 2.2 and 2.5 can be broken into two distinct parts. The first step is to show that Y = f ( WX√ n0 ) ∈ Rn1×m can be viewed as a correlated random matrix with cycle correlations, c.f. Propositions 3.2 and 3.3 below. The second step is to prove the global law for the random matrix M = 1mY Y
∗ with the cycle correlations. In the following, we will sketch the derivation of the self-consistent equation. A more detailed proof is provided in the supplementary material.
The key idea is to use a multivariate cumulant expansion formula. Cumulants of a random vector X = (X1, . . . , Xn) can be defined in a combinatorial way by
κ(X1, . . . , Xn) = ∑ π (−1)|π|−1(|π| − 1)! ∏ B∈π E (∏ i∈B Xi ) , (14)
where the sum runs over all partitions π of the set [n] = {1, . . . , n}, the product runs over the blocks B of the partition π, and |π| is the number of blocks in the partition. The following expansion is commonly referred to as a cumulant expansion and generalises the Gaussian integration by parts. In the context of random matrix theory, the usefulness of this expansion was first observed in [15] and later revived in [12, 13]. A proof of the following lemma is provided in Appendix C for completeness. Lemma 3.1 (Cumulant expansion). If X = (X1, . . . , Xn) is a random vector with finite moments of all orders, then
EX1f(X) = ∑ l≥1 ∑ i1,...,il κ(X1, Xi1 , . . . , Xil) l! E∂i1 · · · ∂ilf(X),
where f : Rn → R is smooth.
We start with the defining identity of the resolvent, 1n1 + zG =MG, where 1n1 denotes the n1 ×n1 identity matrix, and we compute its average trace:
1 + zg = 1
n1 Tr
Y Y ∗G
m =
1
n1 n1∑ i=1 m∑ j=1 Yij ( Y ∗G m ) ji , (15)
where g(z) = 1n1 Tr(M − z1n1) −1 is the normalized trace of the resolvent of M . Since the random variable (Y ∗G)ji can be seen as a function of Yij , we can take the expectation on both sides of (15) and apply Lemma 3.1:
1+zEg = 1
n1 ∑ k≥1 ∑ i1,...,i2k κ(Yi1i2 , Yi3i4 , . . . , Yi2k−1i2k) (k − 1)! E∂Yi3i4 · · · ∂Yi2k−1i2k ( Y ∗G m ) i2i1 . (16)
The main goal now is to show that Y can be viewed as a random matrix with cycle correlations given as in the Propositions 3.2 and 3.3 below: Prop. 3.2 refers to the bias-free case and Prop. 3.3 to the additive bias case. We postpone the proof of both propositions to Subsections A.2 and B.2, resp.
Proposition 3.2 (Correlation structure without bias). The random matrix Y defined by (3) has joint cumulants given by
κ(Yi1i2) = O(n −1/2 0 ),
κ(Yi1i2 , Y ∗ i2i1) ≈ θ1(f),
κ(Yi1i2 , Y ∗ i2i3 , Yi3i4 , . . . , Y ∗ i2ki1 ) ≈ θ2(f)kn1−k0 , k > 1 (17)
where i1, . . . , i2k are all distinct, and we write X ≈ Y as a shorthand notation for X = Y (1 + O(n−1/20 )).
Proposition 3.3 (Correlation structure with bias). The random matrix Y defined by (9) has joint cumulants given by
κ(Yi1i2) = O(n −1/2 0 ),
κ(Yi1i2 , Y ∗ i2i1) ≈ θ1(f), κ(Yi1i2 , Y ∗ i3i1) ≈ θ1,b(f)
κ(Yi1i2 , Y ∗ i2i3 , Yi3i4 , . . . , Y ∗ i2ki1 ) ≈ θ2(f)kn1−k0 , k > 1
(18)
where i1, . . . , i2k are all distinct.
Applying Propositions 3.2 and 3.3 to (16), computing the partial derivatives and doing some bookkeeping, we get the desired equation (6) as n0, n1,m→ ∞. To complete the proofs of Theorems 2.2 and 2.5, one has to show the concentration of g around Eg, as stated in the following lemma. Lemma 3.4. For the random matrix M = 1mY Y
∗ and a complex number z ∈ H such that ℑz > n− 1 4+ϵ
1 , for some ϵ > 0, it holds that
EW |g(z)−EW g(z)|4 ≲ 1
n21(ℑz)4 (19)
with high probability in X , and analogously
EX |g(z)−EXg(z)|4 ≲ 1
n21(ℑz)4 (20)
with high probability in W , where EX (resp. EW ) is the expectation in the X-space (resp. W -space).
The proof of this lemma relies on a standard argument (e.g. see the proof of the concentration inequality in [5, Subsection 3.3.2]) and is given in Appendix D.

heading: 4 Conclusion
text: In this paper, we analysed the singular value distribution of fully random neural networks and found that in the case of additive biases it is impossible to achieve isospectrality by tuning the activation function. In addition, we showed that the resolvent method from random matrix theory also applies to the neural network analysis, despite the non-linearities and we expect that this robust method will prove to be useful in contexts where the conventionally used moment method becomes intractable.
Broader impact
Our result is a purely theoretical one for fully random features, weights and biases. Therefore, we do not expect our contribution to have ethical concerns or adverse future societal consequences.
Acknowledgments and Disclosure of Funding
D. Schröder would like to thank L. Benigni for illuminating discussions on the subject and both authors would like to thank him for his helpful comments on an early version of this preprint. Both authors thank the referees for their careful reading of our manuscript. This work was carried out when the first author was a research assistant at ETH Zurich in the group of W. Werner. The second author is supported by Dr. Max Rössler, the Walter Haefner Foundation and the ETH Zürich Foundation.

Review summaries:
This paper studies the limiting eigenvalue-distribution of the gram matrix $M$ associated with a random feature model, and extends current theory (Pennington and Worah 2017; Benigni and Péché 2019, etc.) to the case of multiple layers with additive layer-wise additive bias. Understanding the distribution of eigenvalues of the gram matrix useful for the understanding of the training dynamics and generalization in neural networks.  

This paper studies the limiting singular value distribution of non-linear random matrix of the form Y=f(WX), where W and X both consist of iid entries and f is applied entrywise, in the proportional growth asymptotics. The limiting law is characterized in terms of its Stieltjes transform as the solution to a self-consistency equation. The same result was obtained before by Pennington and Worah [21] in the Gaussian case and later extended by Benigni and Péché [6] to the general case. The difference is that the method in the current paper is resolvent-based. However, it is unclear what the advantage of this approach is. For example, is the result stronger than that of [6] in the sense of the required assumptions? The authors should clarify the main contribution of this work and what is new.  

This paper studies the spectral distribution of random matrix inspired by a one-hidden layer neural networks. In the first model, the matrix is obtained by applying a nonlinear function entrywise to the product of large rectangular matrices with zero mean IID entries. In the second model, a rank-one bias term is added before applying the non-linearity.  The main results  establish convergence of the emprical spectral distribution to a non-random limiting measure in the regime where the dimensions of the matrices scale linearly. This limiting measure is characterized in terms of the the second moments of the distributions on the matrices, the relative scaling of the dimensions, and certain properties of the non-linearity (described by Gaussian integrals).   While the result for the first model has been derived previously (using moment based methods), the result for the second model is novel. Furthermore, in comparison to the prior work, this paper uses the resolvent method, which requires an evaluation of the cumulants instead of the moments -- a potentially benefit of this approach is that is can be applied more broadly. There is some discussion about the implications of the bias term in the context of multilayer networks. Specially, a non-zero bias preclude the certain invariance properties seen in the setting without bias. 

The authors study the asymptotic eigenvalue distribution of the post-activation Gram matrix in a one-hidden-layer neural network, in a regime where all matrix dimensions are growing proportionally. For input data with i.i.d. centered entries and a weight matrix with i.i.d. centered weights, the authors provide a new proof (under weaker distributional assumptions) of a result first due to Pennington and Worah that the asymptotic eigenvalue distribution is characterized by a quartic equation in its Stieltjes transform. The authors then provide an extension of this result to a model with i.i.d. Gaussian biases in its hidden layer. A qualitative difference in the setting with biases is that---in contrast to the bias-free setting, where special choices of activation function can allow the eigenvalue distribution to remain the same as that of the input Gram matrix, and hence remain constant across multiple layers---the authors show that the eigenvalue distribution with biases must necessarily evolve over layers (cf. Remark 2.6 and the discussion on page 7).  The proof is based on deriving the fixed-point equation for the expected Stieltjes transform using the cumulant-expansion generalization of Gaussian integration-by-parts, and explicitly analyzing the asymptotic mixed cumulants of the entries of the post-activation matrix. A technical challenge is to translate the more easily-obtained understanding of mixed cumulants of the pre-activation entries to the post-activations, and the authors accomplish this using a Fourier-analytic technique.