Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 8: Top 50% of accepted NeurIPS papers, clear accept
review_content: summary: 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.


main_review: Strengths
---

- A stark contrast between the bias-free case and the case with additive bias is shown: in the former case, the eigenvalue-distribution of $M$ converges almost-surely to a compactly supported deterministic probability distribution, whose Stieltjes transform satisfies a certain fourth order equation. In the later case, there is no such convergence (if there is more than 1 hidden layer), irrespective the chose activation function.

- The main technique used to prove the theorems is the so-called "resolvent method". This is in contrast with the very combinatorial and cumbersome "moment method" used in prior works. 

- The paper is well-written and easy to follow. Related works are discussed in great detail, and the authors' contributions are well-situated.

Questions
----
L220 "However, in the case of an additive bias no such convergence can be observed."  Is this a conjecture, an empirical observation, or is this rigorously proved somewhere in the paper ?

Minor Issues
-----
- The super-script numbering of the references is non-standard (page 10).

limitations_and_societal_impact: Yes, this has been taken care of.
review_title: null

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: 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.


main_review: The paper is in general well-written, but I have some difficulty understanding the technical details and I felt a few steps are either missing or skipped. 

- The proof starts with the cumulant expansion in (16) and then applies Prop 3.2. But Prop 3.2 only deals with three special cases of cumulants. Why the contribution of other configurations are negligible? It seems it is assumed the only non-negligible contributions are those from even cycles. For example, it is claimed on p. 15 of the supplement "Finally, we compute the cumulants of the entries of the random matrix Y. Since the activation function f is applied component-wise, it follows from the previous results that the only contribution comes from κ(Yi1, Yi2i3^*, ..., Yi2ki1^*) for k ≥ 1 and i 1 ,...,i 2k distinct, thus proving that Y has cycle correlations. " This is not clear to me why this is the case and what the *previous results* being referred to are. For example, it is clear that odd cycles can be ignored due to the zero-mean assumption. But how about four cycles with double edges, something like \kappa(Y_12,Y_12,Y_23,Y_23,Y_34,Y_34,Y_41,Y_41)?

- Along the same line, in (23), I cannot understand the second line. Note that here the kappa in the first line has k arguments and kappa in the second line has 2k arguments. So this is not a matter of counting. I presume at this point we are removing all other contributions except for even cycles. So this needs to formally justified.

- Also, from the writing it seems in the second line of (23) the indices i_1,...,i_{2k} are required to be distinct. What about the first line and also those summations in Lemma 3.1 and in (16)? None of the these are specified.

- Another crucial step that is not justified is the last step of the proof of Prop 3.2, p16, which claims that the cumulant of the four tuple \approx expectation of their product. Is this true for all other cumulants? In fact, \approx is used throughout the paper without definition.

Other comments:
- line 106: The author should clarify \kappa(Y_{i_1,i_2},Y_{i_2,i_3}^*,...,Y_{i_2k,i_1}^*). Do you mean 
\kappa(Y_{i_1,i_2},Y_{i_2,i_3}^*,Y_{i_3,i_4}^*,Y_{i_4,i_5}^*,...,Y_{i_2k,i_1}^*) or \kappa(Y_{i_1,i_2},Y_{i_2,i_3}^*,Y_{i_3,i_4},Y_{i_4,i_5}^*,...,Y_{i_2k,i_1}^*)? I presume the latter.

- Existence and uniqueness of the solution to the self-consistency equation (1) is not discussed. Perhaps this is in the previous work [21,6]?
- In equation (1) the dummy variable m conflicts with the sample size m which is also used everywhere. One should consider changing the notations.
- It would be good to comment early the meaning of the integral parameters. For example, theta_1(f) is the variance of the entry of Y.
- It is not very clear what "invariance" means as mentioned in the abstract and also in the paper. Perhaps you meant the limiting spectral distributions of Y=f(WX) and X coincide?
- Lemma A.1. What is \mathbb{H}? Is it only used for complex valued W/X?
limitations_and_societal_impact: N/A
review_title: null

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: 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. 
main_review: The high level ideas of this paper are clearly explained. Also, the contribution of the paper is clear in the context of previous work. While the basic idea of the resolvent method itself it not new, this paper shows how it can be used to simplify and generalize the analysis of a problem of current interest. Hence, I think this paper will be of interest to researchers working in this area.

On the technical side, I have some concerns about the lack of higher order moment assumptions on the entries of the matrices. Unless I missed something, the author(s) require only the exists of finite second moments. But clearly some more assumptions are needed justify the CLT statement in lines 125-127. The nonexistence of moment conditions is also potentially problematic later (e.g., in the proof of approximate equality in Proposition 3.2) where higher order terms are said to be negligible. Overall, I would not be surprised if the all of the analysis goes through under some tail conditions (such as the ones given in citation [6]).

There is also room to improve the clarity of the technical results. I list some points below

	-What does the notation less-than-approximately-equal-to mean?

	- What is the definition of 1_m? (I assume identity matrix)
	
	- In line 146 - 147  there are some further statements about the limiting spectral measure (e.g., the nonzero part is continuous). These statements do not seem to follow from Theorem 1 which simply asserts the existence of a limit so I am not sure where they come from.
	
	- In line 148 , the function g(z) is the Stieltjes transform of the empirical measure, not the limit mu, and so the wording is a bit confusing. 
	
	- Remark 2.4 part (ii) says that Theorem 2.1 follows from Theorem 2.2. However, doesn't one also need to establish uniqueness of the solution (in the proper limiting sense)? 
	
	- It would be useful to provide full details of the batch normalization used in Figure 3.

limitations_and_societal_impact: I think the discussion is appropriate.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: 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.
main_review: I think this is a nice contribution to this body of literature: To my knowledge, it is the first work to use this method of proof to study this type of non-linear matrix model. The results of Theorems 2.1 and 2.2 are perhaps not so new in the context of existing work, but the authors illustrate how this method can be used to obtain a new result in the form of Theorem 2.5. There is a high-level component of the proof that I don't understand---see below---but if this is a misunderstanding on my part and can be easily clarified or addressed, my opinion is that this work is above the bar for NeurIPS.

----------------------------------------

Major comment: In Propositions 3.2, 3.3, and Lemma A.3, what exactly is the claim that is being made on which mixed cumulants are 0 and which are non-zero, and how is this used on lines 371 and 456 in the proofs of the main results?

My confusion lies in e.g. the following example: Consider the 4th-order cumulant 1/n_0^2 kappa( (WX)_{1,2}, (WX)_{1,2}, (WX)_{1,2}, (WX)_{1,2} ) where all four entries of WX coincide. Is this case treated in the discussion preceding Lemma A.3, and is Lemma A.3 claiming that this cumulant is 0? My understanding is that it's not 0, but since (WX)_{1,2}/sqrt(n_0) converges to a Gaussian, this is asymptotically 0. However now consider kappa( Y_{1,2}, Y_{1,2}, Y_{1,2}, Y_{1,2} ) for the post-activations: Is Proposition 3.2 claiming that this is 0? Since Y_{1,2} = f((WX)_{1,2}/sqrt(n_0)) converges to a non-Gaussian limit, I guess this is non-zero even in the asymptotic limit. But then perhaps I don't understand how this mixed cumulant (and analogous terms for larger k) is treated in the second equality of line 371. My confusion is related to the sentence "Since the activation function f is applied component-wise, it follows from the previous results... that Y has cycle correlations" on lines 420-423; I'm not sure what is meant by this sentence.

A higher-level comment is that I feel a bit uneasy about the extent of casework that has been done---not only for Lemma A.3 as discussed above, but also in the "Proof of (20)" starting on line 495, and the proofs of Lemmas A.4 and D.1 in appendix D.

------------------------

Other minor comments:

(1) I don't know if it's fair to characterize the moment method as "more commonly used" (in the abstract and conclusion). For example Refs. [19] and [9] also used the resolvent method, and I think an appendix of Ref. [9] derived also the quartic equation of Theorem 2.2 using a resolvent approach (although using a different method to analyze the resolvent).
(2) Line 51, is it common to call gamma here the learning rate?
(3) Can the authors state more precisely the moment conditions needed for nu_1 and nu_2 at the bottom of page 3, for their result to hold.
(4) In Eqs. (18) and (20) of Propositions 3.2 and 3.3 the pattern of indices and conjugations in "..." are not clear, and it would be helpful to write out a few more terms.
(5) Should the approximations of Proposition 3.2 and 3.3 be interpreted for each fixed k as n -> infinity? In the analysis of Eq. (16) one needs to consider the infinite series in k for each fixed n---can the authors elaborate on how the error terms are treated in this series?
limitations_and_societal_impact: N/A
review_title: null

Meta-review: This paper considers the limited singular value distribution of a random feature model, which consists of an entrywise nonlinearity and a product of two i.i.d. matrices. In particular, the generalization performance of a single hidden layer neural network where the input data and hidden layer weights are i.i.d. random can be determined through this analysis. The authors show that the resolvent method can be used to simplify and generalize the existing techniques for this problem. The reviews all agreed that the paper contains strong theoretical results, only expressed minor concerns and suggestions, and recommended acceptance. As Reviewer o1S8 noted, adding a discussion of the advantage of the resolvent approach and clarifying the novelty of the results w.r.t. existing work (e.g., dealing with additive bias) will improve the paper. Please take into account the updated reviews when preparing the final version to accommodate the requested changes. Thank you for your submission to NeurIPS.
Decision: Accept (Poster)