Model, sample, and epoch-wise descents: exact solution of gradient flow in the random feature model
Year: 2021

Abstract:
Recent evidence has shown the existence of a so-called double-descent and even triple-descent behavior for the generalization error of deep-learning models. This important phenomenon commonly appears in implemented neural network architectures, and also seems to emerge in epoch-wise curves during the training process. A recent line of research has highlighted that random matrix tools can be used to obtain precise analytical asymptotics of the generalization (and training) errors of the random feature model. In this contribution, we analyze the whole temporal behavior of the generalization and training errors under gradient flow for the random feature model. We show that in the asymptotic limit of large system size the full time-evolution path of both errors can be calculated analytically. This allows us to observe how the double and triple descents develop over time, if and when early stopping is an option, and also observe time-wise descent structures. Our techniques are based on Cauchy complex integral representations of the errors together with recent random matrix methods based on linear pencils.

heading: 1 Introduction
text: Deep learning models have vastly increased in terms of number of parameters in the architecture and data sample sizes with recent applications using unprecedented numbers with as much as 175 billions parameters trained over billions of tokens [1]. Such massive amounts of data and growing training budgets have spurred research seeking empirical power laws to scale model sizes appropriately with available resources [2], and nowadays it is common wisdom among practitioners that "larger models are better". This ongoing trend has been challenging classical statistical modeling where it is thought that increasing the number of parameters past an interpolation threshold (at which the training error vanishes while the test error usually increases) is doomed to over-fit the data [3]. We refer to [4] for a recent extensive discussion on this contradictory state of affairs. Progress towards rationalizing this situation came from a series of papers [5, 6, 7, 8, 9, 10, 11] evidencing the existence of phases where increasing the number of parameters beyond the interpolation threshold can actually achieve good generalization, and the characteristic U curve of the bias-variance tradeoff is followed by a "descent" of the generalization error. This phenomenon has been called the double descent and was analytically corroborated in linear models [12, 13, 14, 15, 16] as well as random feature (RF) (or random feature regression) shallow network models [17, 18, 19, 20]. Many of these works provide rigorous proofs with precise asymptotic expressions of double descent curves. Further developments have brought forward rich phenomenology, for example, a triple-descent phenomenon [21] linked to the degree of non-linearity of the activation function. Further empirical evidence [22] has also shown that a similar effect occurs while training (ResNet18s on CIFAR10 trained using Adam) and has been
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
called epoch-wise double descent. Moreover the authors of [22] extensively test various CIFAR data sets, architectures (CNNs, ResNets, Transformers) and optimizers (SGD, Adam) and classify their observations into three types of double descents: (i) model-wise double descent when the number of network parameters is varied; (ii) sample-wise double descent when the data set size is varied; and (iii) epoch-wise double descent which occurs while training. We wish to note that sample-wise double descent was derived long ago in precursor work on single layer perceptron networks [23, 24]. An important theoretical challenge is to unravel all these structures in a unified analytical way and understand how generalization error evolves in time.
In this contribution we achieve a detailed analytical analysis of the gradient flow dynamics of the RF model (or regression) in the high-dimensional asymptotic limit. The model was initially introduced in [25] as an approximation of kernel machines; more recently it has been recognized as an important playground for theoretical analysis of the model-wise double descent phenomenon, using tools from random matrix theory [17, 18, 26]. Following [17] we view the RF model as a 2-layer neural network with fixed-random-first-layer-weights and dynamical second layer learned weights. The data is given by n training pairs constituted of d-dimensional input vectors and output given by a linear function with additive gaussian noise. The data is fed through N neurons with a non-linear activation function and followed by one linear neuron whose weights we learn by gradient descent over a quadratic loss function. The high-dimensional asymptotic limit is defined as the regime n, d,N → +∞ while the ratios tend to finite values Nd → ψ and n d → φ. As the training loss is convex one expects that the least-squares predictor (with Moore-Penrose inversion) gives the long time behavior of gradient descent. This has led to the calculation of highly non-trivial analytical algebraic expressions for training and generalization errors which describe (model-wise and sample-wise) double and triple descent curves [17, 21]. However, to the best of our knowledge, there is no complete analytical derivation of the whole time evolution of the two errors.
We analyze the gradient flow equations in the high-dimensional regime and deduce the whole time evolution of the training and generalization errors. Numerical simulations show that the gradient flow is an excellent approximation of gradient descent in the high-dimensional regime as long as the step size is small enough (see Fig. 2). Main contributions presented in detail in Sect. 3 comprise:
a. Results 3.1 and 3.2 give expressions of the time evolution of the errors in terms of one and two-dimensional integrals over spectral densities whose Stieltjes transforms are given by a closed set of purely algebraic equations. The expressions lend themselves to numerical computation as illustrated in Fig. 1 and more extensively in Sect. 3 and the supplementary material.
b. Model and sample-wise double descents develop after some definite time at the interpolation threshold and are preceded by a dip or minimum before the spike develops. This indicates that early stopping is beneficial for some parameter regimes. A similar behavior also occurs for the triple descent. (See Figs. 3, 4 and the 3D version Fig. 1).
c. We observe two kinds of epoch-wise "descent" structures. The first is a double plateau monotonously descending structure at widely different time scales in the largely overparameterized regime (see Fig. 3). The second is an epoch-wise double descent similar to the one found in [22]. In fact, as in [22], rather than a spike, this double descent appears to be an elongated bump over a wide time scale (see Fig. 5 and the 3D version Fig. 1).
Let us say a few words about the techniques used in this work. We first translate the gradient flow equations for the learned weights of the second layer into a set of integro-differential equations for generating functions, as in [27], involving the resolvent of a random matrix (constructed out of the fixed first layer weights, the data, and the non-linear activation). The solution of the integrodifferential equations and the time evolution of the errors can then be expressed in terms of Cauchy complex integral representation which has the advantage to decouple the time dependence and static contributions involving traces of algebraic combinations of standard random matrices (see [28] for related methods). This is the content of propositions 2.0.1 and 2.0.2. With a natural concentration hypothesis in the high-dimensional regime, it remains to carry out averages over the static traces involving random matrices. This is resolved using traces of sub-blocks from the inverse of a larger nontrivial block-matrix, a so-called linear pencil. To the best of our knowledge linear pencils have been introduced in the machine learning community only recently in [29]. This theory is developed in the context of random matrix theory in [30, 31] and [32] using operator valued free-probability. A side-contribution in the SM is also an independent (but non-rigorous) derivation of a basic result of this theory (a fixed point equation) using a replica symmetric calculation from statistical physics.
The non-linearity of the activation function is addressed using the gaussian equivalence principle [33, 34, 29]. Finally, our analysis is not entirely mathematically controlled mainly due to the concentration hypothesis in Sect. 2.3 but comparison with simulations (see Fig. 2 and SM) confirms that the analytical results are exact.
In the conclusion we briefly discuss possible extensions of the present analysis and open problems among which is the comparison with a dynamical mean-field theory approach.

heading: 2 Random feature model

heading: 2.1 Model description
text: Generative model and neural network: We consider the problem of learning a linear function fd(x) = d
− 12 βTx with x, β ∈ Rd column vectors. The vector x is interpreted as a random input and β as a random hidden vector; both with distribution N (0, Id), Id the d× d identity matrix. We assume having access to the hidden function through the noisy outputs y = fd(x) + with additive gaussian noise ∼ N (0, s2), s ∈ R+. We suppose that we have n data-points (xi, yi)1≤i≤n. This data can be represented as the n×dmatrixX ∈ Rn×d where xTi is the i-th row ofX , and the column vector vector Y ∈ Rn with i-th entry yi. Therefore, we have the matrix notation Y = d− 1 2Xβ + ξ where ξ ∼ N (0, s2In) and In the n× n identity matrix. We learn the data with a shallow 2-layer neural network. There are N hidden neurons with weight (column) vectors θi ∈ Rd, i = 1, · · · , N each connected to the d input neurons. Out of these we form the matrix (of the first layer connecting input and hidden neurons) Θ ∈ RN×d where θTi is the i-th row of Θ . Its entries are assumed independent and sampled through a standard gaussian distribution N (0, 1); they are not learned but fixed once for all. The data-points in X are applied linearly to the parameters Θ, and the output Z ∈ Rn×N of the first layer is the pointwise application of an activation function σ : R→ R, Z = σ(d− 12XΘT ). We use the notation zTi to express the i-th row of Z. The second layer consists in a weight (column) vector at ∈ RN to be learned, indexed by time t ≥ 0, with components initially sampled at t = 0 i.i.dN (0, r2), r ∈ R+. The prediction vector is expressed as Ŷt = N− 1 2Zat.
We assume that the activation function belongs to L2(e− x2
2 dx) with inner product denoted 〈 , 〉. It can be expanded on the basis of Hermite polynomials, so σ ∈ Span ((Hek)k≥0), where Hek(x) = (−1)ke x 2 2 dk
dxk e−
x2
2 (so He0(x) = 1, He1(x) = x, He2(x) = x 2 − 1, He3(x) = x3 − 3x, ...).
Furthermore we take σ centered with 〈σ,He0〉 = 0, and set µ = 〈σ,He1〉, ν2 = 〈σ, σ〉 − µ2. For instance, σ = id has (µ, ν) = (1, 0) while σ = Relu − 1√
2π has (µ, ν) = ( 12 , 1 2 (1 − 2 π ) 1/2) ' (0.5, 0.3). Finally, we recall that we are interested in the high dimensional regime where the parameters N,n, d tend to infinity with the ratios Nd → ψ and n d → φ.
Training and test errors: For a new input x0 ∈ Rd, we define the predictor ŷt(x0) = 1√ N z(x0)
Tat where z(x0) = σ( 1√dΘx0). We will further define the standard training and test errors with a penalization term λ > 0 and the quadratic loss:
Htrain(a) = 1 n ∥∥∥Y − Ŷ ∥∥∥2 + λ N ‖a‖2 , Htest(a) = Ex0∼N (0,1) [ (y(x0)− ŷ(x0))2 ] (1)
Note that because of the λ-penalization term, in this context, the training error can be above the test error in some configurations of parameters. Also, we will slightly abuse this notation throughout the paper by usingHtraint ,Htestt to designateHtrain(at),Htest(at).
Gradient flow: Minimizing the training error of this shallow-network is equivalent to a standard Tikhonov regularization problem with a design matrix Z for which the optimal weights are given by a∞ = ( ZTZ N + n N λIN ) −1 ZT√ N Y . The errors generated by the predictors with weights a∞ have been analytically calculated in the high-dimensional regime in [17] and further analyzed in [21]. Here we study the whole time evolution of the gradient flow and thus introduce an additional time dimension in our model. Of course as t→ +∞ one recovers the errors generated by the predictors with weights a∞. The output vector at is updated through the ordinary differential equation datdt = −η∇aH
train(at) with a fixed learning rate parameter η > 0. As η can be absorbed in the time parameter, from now on we consider without loss of generality that η = n2 . Setting δ = λ n N , we find that the gradient flow for at is a first order linear matrix differential equation,
dat dt
= − ( ZTZ
N + δIN
) at +
ZTY√ N . (2)
Recall the initial condition a0 is a vector with i.i.d N (0, r2) components.

heading: 2.2 Cauchy integral representations of the training and test errors
text: An important step of our analysis is the representation ofHtrain andHtest in terms of Cauchy contour integrals in the complex plane. To this end we decompose both errors in elementary contributions and derive contour integrals for each of them. Details are found in section 4 and the SM.
We begin with the test error which is more complicated. We have
Htestt = 1 + s2 − 2µg(t) + µ2h(t) + ν2l(t) + od(1) (3)
where limd→+∞ od(1) = 0 with high probability, and g(t) = β T
√ d ΘT√ d at√ N , h(t) = ‖Θ T √ d at√ N ‖2, and
l(t) = ‖ at√ N ‖2. To describe Cauchy’s integral representation of the elementary functions g, h, l we introduce the resolvent R(z) = (Z TZ N − zIN ) −1 for all z ∈ C \ Sp(Z TZ N ).
Proposition 2.0.1 (Test error) Let Rz be the functional acting on holomorphic functions f : C \ Sp(Z
TZ N ) → C as Rz{f(z)} = − ∮ Γ dz 2πif(z) over a contour Γ encircling the spectrum Sp( ZTZ N ) in the counterclockwise direction. Similarly, let Rx,y be the functional acting on two-variable holomorphic functions f : (C \ Sp(Z
TZ N )) 2 → C as Rx,y{f(x, y)} = ∮ Γ ∮ Γ dx 2πi dy 2πif(x, y). Let
Gt(z) = βT√ d ΘT√ d R(z) at√ N and K(z) = β T √ d ΘT√ d R(z)Z TY N . We have for all t ≥ 0
g(t) = Rz { e−t(z+δ)G0(z) + 1− e−t(z+δ)
z + δ K(z)
} . (4)
Let Lt(z) = aTt√ N R(z) at√ N and Ut(z) = Y TZ N R(z) at√ N and V (z) = Y TZ N R(z) ZTY N . For all t ≥ 0
l(t)=Rz { e−2t(z+δ)L0(z)+2e −t(z+δ) ( 1− e−t(z+δ)
δ + z
) U0(z)+ ( 1− e−t(δ+z)
δ + z
)2 V (z) } . (5)
Let Ht(x, y) = aTt√ N R(x)ΘΘ T d R(y) at√ N , Qt(x, y) = aTt√ N R(x)ΘΘ T d R(y) ZTY N and W (x, y) = Y TZ N R(x) ΘΘT d R(y) ZTY N . For all t ≤ 0
h(t) = Rx,y { e−t(2δ+x+y) ( et(δ+y) − 1 δ + y Q0(x, y) + et(δ+x) − 1 δ + x Q0(y, x) )} +Rx,y { 1− e−t(x+δ)
x+ δ
1− e−t(y+δ)
y + δ W (x, y)
} +Rx,y { e−t(x+y+2δ)H0(x, y) } . (6)
A similar but much simpler representation holds for the training error.
Proposition 2.0.2 (Training error) With the same definitions than in proposition 2.0.1 we have
Htraint = ‖Y ‖2
n +
1 c Rz { (z + δ)e−2t(z+δ)L0(z)−2e−2t(z+δ)U0(z)− 1− e−2t(δ+z) δ + z V (z) } . (7)

heading: 2.3 High-dimensional framework
text: The Cauchy integral representation involves a set of one-variable functions S1 = {G0,K, L0, U0, V } : C \ Sp(Z TZ N ) → C and a set of two-variable functions S2 = {H0,W,Q0} : C \ Sp(Z TZ N ))
2 → C so that g, h, l and thus alsoHtest andHtrain are actually functions of (t;S1,S2). Thus we can write for instance: Htest(at) = Htest(t;S1,S2). We simplify the problem by considering the high-dimensional regime where N,n, d→∞ with ratios Nd → ψ, n d → φ tending to fixed values of order one. In this regime we expect that the functions in S1 and S2 concentrate and can therefore be replaced by their averages over randomness. These averages can be carried out using recent progress in random matrix theory [30], [31], and we are able to compute pointwise asymptotic values of the functions in S1, S2, and eventually substitute them in the Cauchy integral representations for the training and test error. In general, rigorously showing concentration of the various functions involved is not easy and we will make the following assumptions:
Assumptions 2.1 In the high dimensional limit with d,N, n→∞ and Nd → ψ, n d → φ:
1. The random functions in S1, S2 are assumed to concentrate. We let S̄1 = {Ḡ0, K̄, L̄0, Ū0, V̄ } and S̄2 = {H̄0, W̄ , Q̄0} be the pointwise limit of the functions.
2. There exists a bounded subset C ⊂ R+ such that the functions in S̄1 and S̄2 are holomorphic on C \ C and (C \ C)2 respectively
3. The gaussian equivalence principle (see sect. 4.2) can be applied to the limiting quantities.
It is common that the closure of the spectrum of suitably normalized random matrices concentrates on a deterministic set. Thus the bounded set C can be understood as the limit of the finite interval [0, limd max Sp(Z TZ N )]. In the sequel we will distinguish the theoretical high-dimensional regime from the finite dimensional regime using the upper-bar notation.
Definition 2.1 (High-dimensional framework) Under the assumptions 2.1, we define the theoretical test error H̄testt = Htest(t; S̄1, S̄2) and the theoretical training error H̄traint = Htrain(t; S̄1, S̄2)
We conjecture that limdHtraint = H̄traint and limdHtestt = H̄testt at all times t ∈ R. We verify that this conjecture stands experimentally for sufficiently large d on different configurations (see additional figures in the SM). This also lends experimental support on the assumption 2.1. Furthermore we conjecture that the d→ +∞ and t→∞ limits commute, namely limd limtHtraint = limt H̄traint and limd limtHtestt = limt H̄testt .

heading: 3 Results and insights

heading: 3.1 Main results
text: In this section we provide the main results of this work: analytical formulas tracking the test and training errors during gradient flow of the random feature model for all times in the high-dimensional theoretical framework.
Result 3.1 Under the assumption 2.1, the theoretical test and training errors of definition 2.1 are given for all times t ≥ 0 by the formulas
H̄testt = 1 + s2 − 2µḡ(t) + µ2h̄(t) + ν2 l̄(t), (8)
H̄traint = 1 + s2 + 1
c ∫ R [ (δ + ω)e−2t(ω+δ)ρL̄0(ω)− 1− e−2t(δ+ω) δ + ω ρV̄ (ω) ] dω, (9)
with c = φψ , δ = cλ, and the functions ḡ, h̄, l̄ given by
ḡ(t) = ∫ R 1− e−t(ω+δ) ω + δ ρK̄(ω)dω, (10)
l̄(t) = ∫ R [ e−2t(ω+δ)ρL̄0(ω) + ( 1− e−t(ω+δ) ω + δ )2 ρV̄ (ω) ] dω, (11)
h̄(t) = ∫∫ R2 [ e−t(u+v+2δ)ρH̄0(u, v) + 1− e−t(u+δ) u+ δ 1− e−t(v+δ) v + δ ρW̄ (u, v) ] dudv, (12)
where the measures ρK̄ , ρL̄0 , ρV̄ , ρH̄0 , ρW̄ (are possibly signed) are characterized by their Stieltjes transforms given by K̄, L̄0, V̄ , H̄0, W̄{
K̄(x) = tx1 , L̄0(x) = r 2gx1 , V̄ (x) = s 2 (1 + xgx1 ) + (c− hx4) , H̄0(x, y) = r 2q1, W̄ (x, y) = s 2cq4 + q2
(13)
where for each x, y ∈ C+ (the upper half complex plane) gx1 , hx4 , tx1 , g y 1 , h y 4, t y 1 and q1, q2, q4, q5 (which depend symmetrically on (x, y), e.g., q1 = q x,y 1 = q y,x 1 ) are solutions of a purely algebraic system of equations (see SM for the criterion to select the relevant solution) 0 = µψgx1h x 4 − tx1 0 = µψgy1h y 4 − t y 1 0 = (c− 1− xgx1 ) ( c− µ2φgx1hx4 ) − chx4 0 = (c− 1− ygy1 ) ( c− µ2φgy1h y 4 ) − chy4 0 = 1− gx1 ( µ2hx4 + (c− 1− xgx1 )ν2 − x ) 0 = 1− gy1 ( µ2hy4 + (c− 1− yg y 1 )ν 2 − y ) 0 = −µ2gy1q2 + µ2hx4q1 + µg y 1 t x 1 + µg y 1 t y 1 − cg y 1q4ν 2 − gy1 − q1x+ q1ν2 (c− gx1x− 1) 0 = µ (φ− ψgx1x− ψ) (−µgx1 q2 + µh y 4q1 + g x 1 t y 1) + cq4(1− µt y 1)− q2 0 = −µ2φgx1 (1− µtx1)q4 + µ2q5 (c− g y 1y − 1)− ν2φgx1 q4 − φq4 + q1ν2 (φ− ψg y 1y − ψ)
0 = ψ(µ2φgx1g y 1q4 + ψg x 1g y 1 + q1)(1− µt y 1)− µ2ψgx1 q5 (c− gx1x− 1)− q5
We can also deduce the limiting training error and test errors in the infinite time limit:
Result 3.2 In the limit t→∞ we find:
lim t→∞
H̄testt =1 + s2 − 2µK̄(−δ) + µ2W̄ (−δ,−δ) + ν2 dV̄
dx (−δ), lim t→∞ H̄traint =1 + s2 −
1 c V̄ (−δ)
Interestingly, in the limit t → ∞, the expressions become simpler and completely algebraic in the sense that we do not need to compute integrals (or double-integrals) over the supports of the eigenvalue distributions. It is not obvious to see on the analytical expressions that the result is the same as the algebraic expressions obtained in [17] but Fig. 2 shows an excellent match with simulation experiments. We note here that checking that two sets of complicated algebraic equations are equivalent is in general a non-trivial problem of computational algebraic geometry [35].

heading: 3.2 Insights and illustrations of results
text: The set of analytical formulas allows to compute numerically the measures ρK , ρL0 , ρV , ρH0 , ρW and in turn the full time evolution of the test and training errors. The result matches the simulation of a large random feature model where d is taken large as can be seen on Figs. 2 for the infinite time limit (experimental check of result 3.2) and additional figures in the SM (experimental check of result 3.1). Below we illustrate numerical computations obtained with analytical formulas of result 3.1
for various sets of parameters (t, µ, ν, ψ, φ, r, s, λ). For instance, we can freely choose two of these parameters and plot the generalization error in 3D as in Fig. 1, or as a heat-map in the following. We describe three important phenomena which are observed with our analysis.
Double descent and early-stopping benefits: while [17] mostly analyze the minimum least-squares estimator of the random feature model which displays the double-descent at ψ = φ, we are predicting the whole time evolution of the gradient flow as in Fig. 3. We clearly observe the double-descent curve at t = 1010 for ψ = φ; but we now notice that if we stop the training earlier, say at times 1 < t < 10, the generalization error performs better than the minimum least squares estimator. Actually, in the time interval t ∈ (1, 10) for ψ ≈ φ the test error even has a dip or minimum just before the spike develops. We also notice a two-steps descent structure with the test error which is non-existent in the training error and materializes long after the training error has stabilized in the overparameterized regime ψ φ. This is also reminiscent but not entirely similar to the abrupt grokking phenomenon described in [36].
Triple descent: We can observe a triple descent phenomenon materialized by two spikes as seen in Fig. 2 at t =∞ (we also check that the theoretical result matches very well the empirical prediction of the minimum least squares estimator both for training and test errors). This triple descent phenomenon is already contained in the formulas of [17] (although not discussed in this reference) and has been analyzed in detail in [21]. The test error contains a so-called linear spike for φ = 1 (n = d) and a non-linear spike for ψ = φ (N = n). The two spikes are often not seen together as this requires certain conditions to be met, and they tend to materialize together for specific values µ, ν of the activation function where µ ν. Here we further observe the evolution through time of the triple descent and the two spikes and how they develop in Fig. 4. There, we notice that the linear-spike seems to appear earlier than the non-linear one.
Epoch-wise descent structures: Important phenomena that we uncover here are two time-wise "descent structures". (i) As can be seen in Fig. 3, the test error develops a double plateau structure at widely different time scales in the over-parameterized regime (ψ φ) while there seems to be only one time scale for the training error. This kind of double plateau descent is different from the "usual" double-descent. (ii) Moreover, on Fig. 5 for well chosen parameters (in particular for noises with s and r "larger" and ψ = 2φ), we can also observe an elongated bump (rather than a thin spike) for small λ’s. Notice the logarithmic time-scale which clearly shows that here we need to wait exponentially longer to attain the "second descent" after the bump. This is very reminiscent of the epoch-wise double descent described in [22] for deep networks (which happens on similar time scales).

heading: 4 Sketch of proofs and analytical derivations
text: The analysis is threefold. Firstly, we decompose the training and test errors in elementary terms and establish Cauchy’s integral representation for each of them, as provided in proposition 2.0.1. A crucial advantage of this form is that it dissociates a scalar time-wise component and static matrix terms. Secondly, we switch to the high-dimensional framework where the matrix terms are substituted by their limit using the gaussian equivalence principle. Thirdly, we can compute the expectations of
matrix terms thanks to a random matrix technique based on linear pencils. In this section we only sketch the main ideas for each step and provide details in the supplementary material.

heading: 4.1 Cauchy’s integral representation
text: We sketch the derivation for the test error and leave details to appendices. The derivation for the training error is entirely found in the SM. Expanding the square in Equ. (1) and carrying out averages we find Equ. (3) for Htestt with g(t) = βT√ d ΘT√ d at√ N , h(t) = ‖Θ T √ d at√ N ‖2, and l(t) = ‖ at√ N ‖2 (see SM for this derivation).
We show how to derive the Cauchy integral representation for g(t). For h(t), l(t) the steps are similar and found in SM. Let us consider the function (t, z) 7→ Gt(z) as in 2.0.1. Then we have the relation g(t) = −12iπ ∮ Γ dzGt(z) where Γ is a loop in C enclosing the spectrum of Z TZ N . This can easily be seen by decomposing the symmetric Z TZ N in an orthonormal basis v1, . . . , vN with the eigenvalues
λ1, . . . , λN : then we have Gt(z) = ∑N i=1 1 λi−z ( βT√ d ΘT√ d viv T i at√ N ) and because λi are all encircled
by Γ, we find − ∮
Γ dz 2πiGt(z) = ∑N i=1 βT√ d ΘT√ d viv T i at√ N = g(t). Now, the ODE derived for at in (2),
can be written slightly differently using the fact that Z TZ N = R(z) −1 +zI for any z outside Sp(Z TZ N ). Namely, datdt = ZTY√ N −R(z)−1at − (z + δ)at. Then, we can derive an integro-differential equation for Gt(z) involving g(t) and K(z):
∂tGt(z)
∂t = K(z)− g(t)− (z + δ)Gt(z) (14)
In the following, we let L be the Laplace transform operator (Lf)(p) = ∫ +∞
0 dt e−ptf(t), Re p large
enough. Note that the contour integral is performed over a compact set Γ so for Re p large enough, by Fubini’s theorem, the operations L andRz commute. Applying L to (14) and rearranging terms we find for Re(p+ z + δ) 6= 0:
LGp(z) = G0(z)
p+ z + δ +
K(z) p(p+ z + δ) − Lg(p) p+ z + δ
(15)
Now, we can always choose Γ such that−(p+δ) is outside of the contour if we assume Re(p+δ) > 0 (since mini λi ≥ 0). Thus, applyingRz to (15) nullifies the last term because the pole is outside Γ, and using commutativityRzLGp = LRzGp,
RzLGp = LRz { e−t(z+δ)G0(z) + 1− e−t(z+δ)
z + δ K(z)
} = Lg(p). (16)
Finally, using the inverse Laplace transform leads to (4).

heading: 4.2 Gaussian equivalence principle
text: The matrix terms must be estimated in the limit d → ∞ with {β, a0, ξ,Θ, X} all independently distributed. As per assumptions 2.1 all the matrix terms in S1,S2 are assumed to concentrate. So for instance we assume that the following limit exists K̄(z) ≡ limd→∞K(z) = limd→∞ Eβ,ξ,Θ,X [K(z)]. Using cyclicity of the trace we easily perform averages over β, ξ to find
K̄(z) = lim d
Eβ,Θ,XTr [
ΘT√ d R(z) ZTX N ββT d
] = lim
d
1 d Eβ,Θ,XTr [ ΘT√ d R(z) ZTX N ] . (17)
After these reductions, the expressions of all functions in S̄1, S̄2 essentially involve products of random matrices Θ, X and pointwise applications of the non-linear activation σ. This can be further reduced to simpler algebraic expressions using the gaussian equivalence principle. This principle states that: there exists a standard gaussian random matrix Ω ∈ Rn×N independent of {X,Θ} such that in the infinite dimensional limit we can make the substitution Z = σ ( d−1/2XΘT ) −→ µd−1/2XΘT + νΩ in the expressions of all functions in S̄1, S̄2. This approach is quite general and is well described in [29, 37] (and formerly in [33] and [34]). Thus it remains to compute expectations of traces containing only products, and inverses of products and sums, of gaussian matrices.

heading: 4.3 Expectations over random matrices using linear pencils
text: We explain how to compute the limit of (17) once the gaussian equivalence principle has been applied. A powerful approach is to design a so-called linear pencil. In the present context this is a suitable block-matrix containing gaussian random matrices and multiples of the identity matrix, for which full block-inversion gives back the products of terms in the traces that are being sought. This approach has been described in [30, 31, 38]. We have found a suitable linear pencil which contains fortuitously all the terms required in S̄1, S̄2. It is described by the 13× 13 block-matrix M , and pursuing with
our example, we get for instance with the block (7, 12) that limdK(y) = limd 1dTr[(M −1)(7,12)]
M=  −xI −µ Θ√ d −I 0 0 0 Θ√ d 0 0 0 0 0 0 0 I 0 X T √ N 0 0 0 0 0 0 0 0 0 0 0 I ν Ω T √ N 0 0 0 0 0 0 0 0 0 0 0 0 I X√ N ν Ω√ N 0 0 0 0 0 0 0 µΘ T √ d 0 0 0 I 0 0 0 0 0 0 0 0 I 0 0 0 0 I 0 0 0 0 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 Θ T √ d 0 0 0 0 0 0 0 I 0 ν Ω T √ N 0 0 0 0 0 0 0 0 0 0 0 I X T √ N 0 0 0 0 0 0 0 0 0 0 0 0 I ν Ω√ N X√ N 0 0 0 0 0 0 0 0 0 0 0 I 0 −I 0 0 0 0 0 0 0 0 0 0 0 I −µΘ
T √ d
0 0 0 0 0 0 0 I µ Θ√ d
0 0 0 −yI  Next, the great advantage of the linear pencil is that (as described in [30, 31, 38]) it allows to write a
fixed point equation F (G) = G for a "small" 13× 13 matrix G with scalar matrix elements. We also provide in the SM an independent derivation of the fixed point equations using the replica method (a technique from statistical physics [39]).The components of G are linked to the limiting traces of the blocks of M−1 as in [G]2,7 = K̄(z). The action of F can be completely described as an algebraic function leading to (a priori) 13× 13 = 169 equations over the matrix elements of G. The number of equations can be immediately reduced to 39 because many elements vanish, and with the help of a computer algebra system the number of equations can be further brought down to 10. We refer to the SM for all the details about the method.

heading: 5 Conclusion
text: We believe that our analysis could be extended to study the learning of non-linear functions, the effect of multilayered structures, and potentially different layers such as convolutions, as long as they are not learned. A challenging task is to extend the present methods to learned multilayers. A further question is the application of our analysis in teacher-student scenarios with realistic datasets (See [40, 37]).
Finally we wish to point out that a comparison of the approach of the present paper (and the similar but simpler one of [27]) with the dynamical mean-field theory (DMFT) approach of statistical physics remains to be investigated. DMFT has a long history originating in studies of complex systems (turbulent fluids, spin glasses) where one eventually derives a set of complicated integro-differential equations for suitable correlation and response functions capturing the whole dynamics of the system (we refer to the recent book [41] and references therein). This is a powerful formalism but the integral equations must usually be solved entirely numerically which itself is not a trivial task. For problems close to the present context (neural networks,generalized linear models, phase retrieval) DMFT has been developed in the recent works [42, 43, 44, 45, 46]. We think that comprehensively comparing this formalism with the present approach is an interesting open problem. It would be desirable to connect the DMFT equations to our closed form solutions for the training and generalization errors expressed in terms of a set of algebraic equations of suitable Stieltjes transforms.
Acknowledgments and Disclosure of Funding
The work of A. B has been supported by Swiss National Science Foundation grant no 200020 182517.

Review summaries:
In this paper, authors proposed an analytical model to understand the multiple-descent curves in generalization errors of over-parameterized models (random feature models optimized over gradient flow). More specifically, they provided exact solutions to interpret the model-wise, epoch-wise and sample-wise double descent phenomenon.   In Section 2., authors define the learning problem of over-parameterized models with random features as a regression problem of Ridge-alike, then define the optimization procedure as an ODE. In Section 2.2., authors represent  training and testing errors with  Cauchy integral representations with assumptions prespecified. They report the result as the solution of ODE under assumptions while providing some simulation results to backup their theoretical results. The sketched proofs are attached. 

The authors analytically study model, sample, and epochwise descents of the random feature model. Model descent meaning how the test error changes as the model complexity is varied; sample descent meaning how the test error changes as the number of samples are varied; and epochwise descent meaning how the test error changes as the model trains. They derive analytic expressions for the test and train errors and study what happens when the number of neurons $N$, the number of samples $n$, and the sample dimensionality $d$ all tend to infinity. Since this can be done in varying ways, they study when $\frac{N}{d} \rightarrow \psi$ and $\frac{n}{d} \rightarrow \phi$ for fixed values of $\psi$ and $\phi$. They observe double descent which is beneficial (lower test error when $\psi \gg \phi$) if training is performed for a long time. On the other hand, if training is cut short (early stopping) the test error is low and even has a dip where later there develops a peak. The authors also observe triple descent if the training time tends to $\infty$.

This article provides, modulo some important technical assumptions, computationally tractable formulas for the test and train error of a random features model trained by gradient flow on \ell_2-regularized MSE during the entire course of training. The features are those computed by the first layer of a randomly initialized neural network and the regime of interest for this article is a triple scaling limit in which the input dimension, dataset size, and hidden layer width tend to infinity but their ratios remain order 1. Unlike prior work this article provides formulas for the evolution in time of both the training and test error, instead only of only their behaviors are t = \infty. While the results are not completely rigorous (they assume concentration for spectral statistics of certain random matrices), I found them intriguing and am fairly convinced by plots in which the predictions fit empirics quite well. The main weakness for me is a lack of further formulation of takeaways from the theoretical results (see weakness #2). Overall, I thought this was a solid paper. 

This paper presents an analytical derivation of the train and test loss curves of the random feature model. The main novelty lies in the fact that the authors not only capture the behavior at final times, but also the dynamics. This leads to various interesting observations, since such models have been shown to exhibit epoch-wise double descent.