Introduction
The data collected at the Large Hadron Collider (LHC) up to the now-completed third data-taking period (Run 3) have brought the field of particle physics into a regime where sub-percent accuracy is becoming the benchmark for both experimental measurements and theoretical predictions. These data have underpinned the most recent validations of the Standard Model (SM) through the measurement of its parameters [1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 ] , making precision itself a non-negotiable condition to probe physics beyond the SM (BSM) [11 ] .
These results, however, are rarely read off the raw data straight away. The quantities of interest have to be reconstructed by solving an inverse problem [12 ] , so that the reliability of the determination is bounded not only by the data but also by how well the methodology – the inference procedure – is kept under control. Many methodologies rest on machine learning (ML) tools [13 ] , and their reliability – uncertainty quantification, calibration, and validation – has itself drawn significant attention [14 , 15 , 16 , 17 , 18 ] . With the High-Luminosity LHC (HL-LHC) era [19 ] and promising future collider facilities [11 , 20 ] on the horizon, the requirements on the methodology will become even more stringent.
Meeting these requirements is, however, inherently difficult. Most of the inverse problems that challenge particle physics are indeed plagued by their ill-posed nature [12 ] ; the determination of parton distribution functions (PDFs) from global analyses of experimental data is no exception. PDFs are fundamental quantities that encode the inner structure of a hadron – the proton, in our case – and their extraction from data relies on an intricate interplay between three main ingredients: the data that constrain the model, the theory that predicts the data, and the methodology that dictates how the two combine. The problem is ill-posed: a finite amount of data cannot pin down the continuous functional form of the PDFs, leaving the methodology to fix the unconstrained directions. As in the case of unfolding [21 , 22 ] , any methodology therefore specifies a set of implicit and explicit assumptions that regularise the problem, the former being more often than not hidden within the methodology itself, and therefore difficult to spell out precisely.
That the choice of methodology has a significant impact on the result is now well established. Comparisons across different methodologies reveal sizeable differences both at the level of PDFs [23 , 24 ] and of the theoretical predictions [25 ] , even though the data entering the fits are largely the same. Understanding the origin of these deviations is crucial because PDFs are the primary input to essentially every theoretical prediction that involves a proton in the initial state. Moreover, they are used in the inference of SM parameters [26 , 27 ] , and are at the core of indirect searches for new physics, where an inaccurate determination of the PDFs can absorb a BSM signal [28 , 29 ] .
It is thus not surprising that, together with the remarkable achievements in the theoretical calculations now approaching N 3 L O N^{3}\mathrm{LO} N 3 LO [30 , 31 , 32 ] , the faithfulness of delivered PDF uncertainties is becoming the subject of active investigations [33 , 34 , 35 ] . Pinning down which features of the final results are dictated by the data, and which are inherited from the methodology, is key to distinguishing genuine physical signals from uncontrolled artefacts of the inference procedure.
Efforts to disentangle the role of the methodology in PDF determinations are not lacking. The Colibri framework [36 ] allows different methodologies and inference strategies to be compared within the same infrastructure. Bayesian approaches [37 , 38 , 39 , 40 ] have proven successful in laying out a framework where the assumptions are explicit – and so are their implications – while tools originally developed in the context of explainable AI have been applied to PDFs either to expose the inner features of the PDFs [41 , 42 ] or to make the role of the data more transparent [43 ] . Most of these studies, however, propose alternative ways to approach the problem, while a systematic study of the methodologies that set the current standard in PDF determination is still lacking.
We began to address this point in Ref. [44 ] where we focused on the NNPDF4.0 methodology [45 , 46 ] , which rests on the propagation of the experimental uncertainties into PDFs through Monte Carlo sampling of the data and on the parametrisation of PDFs by means of a neural network (NN). Using a single-flavour toy model, we studied the training dynamics of the NN by inspecting the evolution of the neural tangent kernel (NTK) [47 , 48 , 49 ] . We showed that, under the assumption of a lazy training regime [50 ] , characterised by a slowly varying NTK, the gradient flow can be integrated exactly, providing an analytical description of the training from the onset of the lazy regime onwards. Whether these results hold in a more realistic fit remained, however, an open question.
Our contribution
In this work we extend the previous analysis to a realistic fit of PDFs, based on a comprehensive dataset of deep inelastic scattering (DIS) measurements that allows for the simultaneous fit of eight independent flavours. We adopt the NNPDF4.0 methodology with plain gradient descent, with no sum rules and positivity constraints, and generate synthetic data from a known set of PDFs, so that the results can be validated against the underlying law. Within these settings, we aim at quantifying the impact of the implicit assumptions discussed above and at understanding the training dynamics of the NN. Our contribution can be summarised as follows.
We first generalise the analytical solution of the gradient flow to the multi-flavour case and validate it against the empirically trained network. The implicit regularisation of the inverse problem is thus made manifest as the projection of the training flow onto the sub-space of non-vanishing modes of the NTK. In this picture, training is viewed as an incremental inclusion of NTK modes, each learned at the timescale set by the corresponding eigenvalue. We show that only a few of the leading modes of the NTK spectrum are needed to reconstruct the underlying law, while the sub-leading ones contribute to the noise of the solution. The choice of when to stop the training – usually delegated to a stopping criterion – thus becomes a trade-off between the signal and noise of the data. The empirical analysis is supplemented with a derivation of the Fisher information matrix of the PDFs, from which we recover, for DIS observables, the Shapley value of Ref. [43 ] . We then outline a geometric interpretation of the NTK and of the analytical solution.
Using the analytical solution, we demonstrate that the central value and most of the error budget of the PDFs are driven by the data. The initial condition, although it does not bias the central value, contributes to the uncertainty by an amount that depends on the flavour and on the region in x x x
The NTK was first introduced in Ref. [47 ] , which showed that, in the infinite-width limit in which the network at initialisation is a Gaussian process [51 , 52 , 53 ] , the NTK is deterministic and constant during training, and gradient descent reduces to kernel learning. This was then recast in weight space [48 ] , showing that wide NNs train like their linearisation around the parameters at initialisation. The properties of the NTK during training have been studied in Refs. [49 , 54 ] , while examples of methods borrowed from statistical and quantum field theory used to describe wide networks and their finite-width corrections can be found in Refs. [55 , 56 , 57 ] . The connection with the Fisher information matrix and the geometric interpretation of § 3 were inspired by Refs. [58 , 59 , 60 ] . Early stopping as a regularisation has been studied in linear models in Refs. [61 , 62 ] . The study of the dependence on the initialisation was inspired by Refs. [63 , 64 , 65 ] .
The remainder of this paper is organised as follows. In § 2 we review the dataset and the methodology, generalise the analytical solution of the gradient flow to the multi-flavour case, and use it to interpret early stopping as a bias-variance trade-off. § 3 is devoted to the derivation of the Fisher information matrix of the PDFs and to the geometric interpretation of the NTK. In § 4 we present the NTK analysis of the DIS-only fit and compare the analytical solution with the trained network, while in § 5 we study how the training dynamics depends on the initialisation. We summarise our findings in § 6 , where we also propose future improvements and applications. A closed-form expression of the NTK for deep linear networks is derived in Section A, while in Section B we collect supplementary material to complement the discussion of § 5 .
Methodological and theoretical framework
In this section, we provide a brief overview of the methodological and experimental ingredients that enter this analysis. We then summarise in § 2.2 the main results of Ref. [44 ] and show how the analytical solution of the gradient flow generalises to the case of a multi-flavour PDF fit. In doing so, we introduce the NTK and the operator M M M § 3 . Finally, in § 2.6 we interpret early stopping through the lens of the NTK, thus exposing the role of each NTK mode in the bias-variance trade-off along the training dynamics.
Dataset and methodology
The experimental input comprises all the DIS measurements included in the NNPDF4.0 global analysis, both neutral- and charged-current, restricted to observables whose theoretical predictions depend linearly on the PDFs. After kinematic cuts, essential for ensuring accuracy in the theoretical predictions, the dataset amounts to N d a t a = 2968 N_{\mathrm{data}} = 2968 N data = 2968 C Y C_Y C Y [45 ] and references therein for a detailed description of each experimental input entering the dataset. In the following, we provide a concise account of the methodological aspects involved in this analysis, summarised in Table 1 .
Table 1. Settings of the baseline fit.
Theoretical predictions. These are computed as in the previous study, although the multi-flavour structure introduced here adds a layer of complexity to the NTK analysis. In fact, the FK tables are now rank-3 tensors that approximate the convolution integral of the factorisation theorem by way of a Riemann sum over a grid of points x α x_\alpha x α α = 1 , … , N g r i d \alpha= 1,\ldots,N_{\mathrm{grid}} α = 1 , … , N grid 1 f f f N f l a v = 14 N_{\mathrm{flav}} = 14 N flav = 14 x x x
T I = ∑ α = 1 N g r i d ∑ i = 1 N f l a v ( F K ) I , α i f i α , (1) T_I = \sum_{\alpha=1}^{N_{\mathrm{grid}}} \sum_{i=1}^{N_{\mathrm{flav}}} (\mathrm{FK})^i_{I,\alpha} f_{i\alpha}, \tag*{(1)} T I = α = 1 ∑ N grid i = 1 ∑ N flav ( FK ) I , α i f i α , (1) where I = 1 , … , N d a t a I = 1,\ldots,N_{\mathrm{data}} I = 1 , … , N data i = 1 , … , N f l a v i = 1,\ldots,N_{\mathrm{flav}} i = 1 , … , N flav N f l a v N_{\mathrm{flav}} N flav x f xf x f f f f f i α ≡ x α f i ( x α , Q 0 ) f_{i\alpha} \equiv x_\alpha f_i(x_\alpha,Q_0) f i α ≡ x α f i ( x α , Q 0 ) flavour-major ordered index a = 1 , … , N v e c a = 1,\ldots,N_{\mathrm{vec}} a = 1 , … , N vec N v e c = N f l a v × N g r i d N_{\mathrm{vec}} = N_{\mathrm{flav}} \times N_{\mathrm{grid}} N vec = N flav × N grid
T I = ∑ a = 1 N v e c ( F K ) I a f a , or T = ( F K ) f , (2) T_I = \sum_{a=1}^{N_{\mathrm{vec}}} (\mathrm{FK})_{Ia} f_a,\quad\text{or}\quad\mathbf{T} = (\mathrm{FK})\mathbf{f}, \tag*{(2)} T I = a = 1 ∑ N vec ( FK ) I a f a , or T = ( FK ) f , (2) where the FK table is reshaped as a matrix
( F K ) = [ F K Σ F K g … F K V 35 ] ∈ R N d a t a × N v e c , F K i ∈ R N d a t a × N g r i d ∀ i ∈ { Σ , g , … , V 35 } , (3) (\mathrm{FK}) = [\mathrm{FK}^{\Sigma}\ \mathrm{FK}^{g}\ \ldots\ \mathrm{FK}^{V_{35}}] \in\mathbb{R}^{N_{\mathrm{data}}\times N_{\mathrm{vec}}},\quad\mathrm{FK}^{i} \in\mathbb{R}^{N_{\mathrm{data}}\times N_{\mathrm{grid}}}\quad\forall i \in\{\Sigma,g,\ldots,V_{35}\}, \tag*{(3)} ( FK ) = [ FK Σ FK g … FK V 35 ] ∈ R N data × N vec , FK i ∈ R N data × N grid ∀ i ∈ { Σ , g , … , V 35 } , (3) and PDFs are cast into a single column vector
f = [ f Σ f g ⋮ f V 35 ] , f i ∈ R N g r i d ∀ i ∈ { Σ , g , … , V 35 } . (4) \mathbf{f} =
\begin{bmatrix}
f_{\Sigma}\\
f_g\\
\vdots\\
f_{V_{35}}
\end{bmatrix},\quad f_i \in\mathbb{R}^{N_{\mathrm{grid}}}\quad\forall i \in\{\Sigma,g,\ldots,V_{35}\}. \tag*{(4)} f = f Σ f g ⋮ f V 35 , f i ∈ R N grid ∀ i ∈ { Σ , g , … , V 35 } . (4) The NNPDF4.0 parametrisation. The PDFs are parametrised by means of a feed-forward NN with two hidden layers having 25 and 20 nodes, respectively. The activation function is the hyperbolic tangent in the hidden layers, while a linear function is used in the output layer. We parametrise N o u t = 8 N_{\mathrm{out}} = 8 N out = 8 Σ \Sigma Σ g g g V V V V 3 V_{3} V 3 V 8 V_{8} V 8 T 3 T_{3} T 3 T 8 T_{8} T 8 T 15 T_{15} T 15 Q 0 = 1.65 G e V Q_{0} = 1.65\ \mathrm{GeV} Q 0 = 1.65 GeV x 1 − δ ( 1 − x ) ϵ x^{1-\delta}(1-x)^{\epsilon} x 1 − δ ( 1 − x ) ϵ x x x x x x x x x ( x , log x ) (x,\log x) ( x , log x ) x x x N N α j = N N ( x α ) j \mathrm{NN}_{\alpha j} = \mathrm{NN}(x_{\alpha})_{j} NN α j = NN ( x α ) j j j j x α x_{\alpha} x α
f i α = A i ∑ j = 1 N o u t x α 1 − δ j ( 1 − x α ) ϵ j N N α j R j i ( F K ) , (5) f_{i\alpha} = A_i \sum_{j=1}^{N_{\mathrm{out}}} x_{\alpha}^{1-\delta_j}(1-x_{\alpha})^{\epsilon_j}\mathrm{NN}_{\alpha j}R_{ji}^{(\mathrm{FK})}, \tag*{(5)} f i α = A i j = 1 ∑ N out x α 1 − δ j ( 1 − x α ) ϵ j NN α j R j i ( FK ) , (5) where R ( F K ) R^{(\mathrm{FK})} R ( FK ) N o u t N_{\mathrm{out}} N out N f l a v = 14 N_{\mathrm{flav}} = 14 N flav = 14 A i A_i A i V V V V 3 V_{3} V 3 V 8 V_{8} V 8
A g = 1 − ∫ 0 1 d x x Σ ^ ( x ) ∫ 0 1 d x x g ^ ( x ) , A V = 3 ∫ 0 1 d x V ^ ( x ) , A V 3 = 1 ∫ 0 1 d x V ^ 3 ( x ) , A V 8 = 3 ∫ 0 1 d x V ^ 8 ( x ) . (6) A_g = \frac{1-\int_{0}^{1} \mathrm{d}x\,x\widehat{\Sigma}(x)}{\int_{0}^{1} \mathrm{d}x\,x\widehat{g}(x)}, \qquad
A_V = \frac{3}{\int_{0}^{1} \mathrm{d}x\,\widehat{V}(x)}, \qquad
A_{V_3} = \frac{1}{\int_{0}^{1} \mathrm{d}x\,\widehat{V}_3(x)}, \qquad
A_{V_8} = \frac{3}{\int_{0}^{1} \mathrm{d}x\,\widehat{V}_8(x)}. \tag*{(6)} A g = ∫ 0 1 d x x g ( x ) 1 − ∫ 0 1 d x x Σ ( x ) , A V = ∫ 0 1 d x V ( x ) 3 , A V 3 = ∫ 0 1 d x V 3 ( x ) 1 , A V 8 = ∫ 0 1 d x V 8 ( x ) 3 . (6) Here a hat denotes the un-normalised distribution, i.e. x f ^ i ( x ) x\widehat{f}_i(x) x f i ( x ) (5) with A i = 1 A_i = 1 A i = 1 x Σ x\Sigma x Σ x g xg xg x x x § 2.2 . For this reason, in the following we do not impose the sum rules and set A i = 1 A_i = 1 A i = 1
Monte Carlo replicas. Experimental uncertainties are propagated to the fitted PDFs by means of the Monte Carlo replica method [33 ] . We adopt the closure-test procedure and train an ensemble of neural networks on a set of synthetic data generated using a known underlying PDF set f ( i n ) f^{(\mathrm{in})} f ( in ) [45 ] . Working with synthetic data, for which we know the underlying true law, allows us to assess the performance of the model without introducing spurious effects of potentially inconsistent data [34 ] . We work throughout with level-2 (L2) pseudodata. These are generated by first smearing the predictions computed with the true law f ( i n ) f^{(\mathrm{in})} f ( in ) η \eta η ξ ( k ) \xi^{(k)} ξ ( k ) k k k
Y ( k ) = ( F K ) f ( i n ) + η + ξ ( k ) , η , ξ ( k ) ∼ N ( 0 , C Y ) , (7) Y^{(k)} = (\mathrm{FK})f^{(\mathrm{in})} + \eta+ \xi^{(k)}, \qquad\eta,\xi^{(k)} \sim\mathcal{N}(0,C_Y), \tag*{(7)} Y ( k ) = ( FK ) f ( in ) + η + ξ ( k ) , η , ξ ( k ) ∼ N ( 0 , C Y ) , (7) where k k k [44 ] and to Refs. [37 , 45 ] for a detailed treatment of the closure-test methodology.
Initialisation of the NN parameters. At initialisation, the exponents of the small- and large-x x x Table 2 . These ranges were chosen in the hyperoptimisation run of Ref. [45 ] , to which we refer for a detailed discussion. For the NN, the weights are drawn from the Glorot normal distribution [66 ] , while the biases are deliberately set to zero. The Glorot normal distribution is a zero-mean Gaussian distribution with a layer-dependent variance such that
Table 2. Ranges of the preprocessing exponents.
E [ W i j ( ℓ ) ] = 0 , E [ W i j ( ℓ ) W i ′ j ′ ( ℓ ′ ) ] = δ i i ′ δ j j ′ δ ℓ ℓ ′ σ ℓ 2 , σ ℓ 2 = 2 n ℓ − 1 + n ℓ , (8) \mathbb{E}\left[W_{ij}^{(\ell)}\right]=0,\qquad\mathbb{E}\left[W_{ij}^{(\ell)}W_{i'j'}^{(\ell')}\right]=\delta_{ii'}\delta_{jj'}\delta_{\ell\ell'}\sigma_{\ell}^{2},\qquad\sigma_{\ell}^{2}=\frac{2}{n_{\ell-1}+n_{\ell}}, \tag*{(8)} E [ W ij ( ℓ ) ] = 0 , E [ W ij ( ℓ ) W i ′ j ′ ( ℓ ′ ) ] = δ i i ′ δ j j ′ δ ℓ ℓ ′ σ ℓ 2 , σ ℓ 2 = n ℓ − 1 + n ℓ 2 , (8) where W i j ( ℓ ) W_{ij}^{(\ell)} W ij ( ℓ ) ℓ \ell ℓ i = 1 , … , n ℓ i=1,\ldots,n_{\ell} i = 1 , … , n ℓ j = 1 , … , n ℓ − 1 j=1,\ldots,n_{\ell-1} j = 1 , … , n ℓ − 1 n ℓ n_{\ell} n ℓ n ℓ − 1 n_{\ell-1} n ℓ − 1 ℓ \ell ℓ 2 § 5 .
Loss function and diagonal basis. Each replica of the model is trained to the corresponding replica of the data by minimising an appropriate loss function. Unlike the NNPDF4.0 methodology, here the fit is performed in the diagonal basis of the experimental covariance matrix C Y C_Y C Y C Y = R D R T C_Y=RDR^T C Y = R D R T D D D C Y C_Y C Y R R R C Y C_Y C Y
Y ~ = R T Y , T ~ = R T T . (9) \widetilde{Y}=R^T Y,\qquad\widetilde{T}=R^T T. \tag*{(9)} Y = R T Y , T = R T T . (9) For a given replica k k k D D D
L ( k ) = 1 N d a t ∥ Y ~ ( k ) − T ~ ( k ) ∥ D − 1 2 + L p o s + L i n t , (10) \mathcal{L}^{(k)}=\frac{1}{N_{\mathrm{dat}}}\left\|\widetilde{Y}^{(k)}-\widetilde{T}^{(k)}\right\|_{D^{-1}}^{2}+\mathcal{L}_{\mathrm{pos}}+\mathcal{L}_{\mathrm{int}}, \tag*{(10)} L ( k ) = N dat 1 Y ( k ) − T ( k ) D − 1 2 + L pos + L int , (10) where ∥ ⋅ ∥ D − 1 2 \|\cdot\|_{D^{-1}}^{2} ∥ ⋅ ∥ D − 1 2 D D D L p o s \mathcal{L}_{\mathrm{pos}} L pos L i n t \mathcal{L}_{\mathrm{int}} L int (10) are Lagrange multipliers that constrain the PDFs to satisfy the positivity and integrability conditions, respectively. However, because the analytical treatment of the training dynamics has not yet been extended to include these constraints, and because we also want to keep the fitting procedure as close as possible to the previous study, these terms are neglected in this analysis.
Minimisation and early stopping. The loss is minimised using vanilla gradient descent (GD), and no stochasticity is introduced in the training. This is done to keep the training steps as close as possible to the analytical solution of the gradient flow, which is derived without any stochasticity, momentum terms, or adaptive learning rates. Only the parameters of the neural network, namely weights and biases, are updated in the minimisation, while the preprocessing factors are kept fixed after initialisation. In order to prevent over-fitting, we split the dataset into a training and a validation set – the former is minimised during training while the latter is monitored. The splitting is performed at the level of the diagonalised dataset, where the data points are uncorrelated, so that the split does not spoil the correlation structure of the original data. The training is left to run for a fixed number of epochs, τ e n d \tau_{\mathrm{end}} τ end
The analytical solution of the gradient flow
The training of the network amounts to minimising the loss function in eq. (10) by iterative steps of GD. For the purpose of the analytical treatment, it is convenient to describe the discrete parameter updates as a continuous gradient flow , in which the parameters θ μ \theta_{\mu} θ μ t t t
d θ t , μ d t = − ∂ μ L t , (11) \frac{d\theta_{t,\mu}}{dt}=-\partial_{\mu}\mathcal{L}_{t}, \tag*{(11)} d t d θ t , μ = − ∂ μ L t , (11) where ∂ μ = ∂ / ∂ θ μ \partial_\mu= \partial/\partial\theta_\mu ∂ μ = ∂ / ∂ θ μ t t t θ μ \theta_\mu θ μ t t t η \eta η [67 ] , with the training time t t t τ \tau τ t = η τ t = \eta\tau t = η τ (9) .
The continuous gradient flow
Because the theory predictions are linear in the PDFs, the derivative of the predictions T T T ∂ T / ∂ f = ( F K ) \partial T/\partial f = (\mathrm{FK}) ∂ T / ∂ f = ( FK )
d f t d t = Θ t ( F K ) T C Y − 1 ( Y − ( F K ) f t ) , (12) \frac{d f_t}{d t} = \Theta_t (\mathrm{FK})^T C_Y^{-1} \left(Y - (\mathrm{FK}) f_t\right), \tag*{(12)} d t d f t = Θ t ( FK ) T C Y − 1 ( Y − ( FK ) f t ) , (12) where the NTK Θ \Theta Θ [44 , 47 ] is defined as
Θ t , i 1 α 1 i 2 α 2 = ∑ j 1 , j 2 = 1 N o u t R j 1 i 1 ( F K ) R j 2 i 2 ( F K ) x α 1 1 − δ j 1 ( 1 − x α 1 ) ϵ j 1 x α 2 1 − δ j 2 ( 1 − x α 2 ) ϵ j 2 ∑ μ ∂ μ N N t , α 1 j 1 ∂ μ N N t , α 2 j 2 . (13) \Theta_{t,i_1\alpha_1i_2\alpha_2} = \sum_{j_1,j_2=1}^{N_{\mathrm{out}}} R^{(\mathrm{FK})}_{j_1i_1} R^{(\mathrm{FK})}_{j_2i_2} x_{\alpha_1}^{1-\delta_{j_1}}(1-x_{\alpha_1})^{\epsilon_{j_1}} x_{\alpha_2}^{1-\delta_{j_2}}(1-x_{\alpha_2})^{\epsilon_{j_2}} \sum_\mu\partial_\mu\mathrm{NN}_{t,\alpha_1j_1}\,\partial_\mu\mathrm{NN}_{t,\alpha_2j_2}. \tag*{(13)} Θ t , i 1 α 1 i 2 α 2 = j 1 , j 2 = 1 ∑ N out R j 1 i 1 ( FK ) R j 2 i 2 ( FK ) x α 1 1 − δ j 1 ( 1 − x α 1 ) ϵ j 1 x α 2 1 − δ j 2 ( 1 − x α 2 ) ϵ j 2 μ ∑ ∂ μ NN t , α 1 j 1 ∂ μ NN t , α 2 j 2 . (13) The NTK carries two pairs of flavour and grid indices. With the same flavour-major vectorisation used for the FK tables and the PDFs, each pair is collapsed into a single composite index, so that Θ t \Theta_t Θ t
Θ t ∈ R N v e c × N v e c , (14) \Theta_t \in\mathbb{R}^{N_{\mathrm{vec}}\times N_{\mathrm{vec}}}, \tag*{(14)} Θ t ∈ R N vec × N vec , (14) acting on the same space as the vectorised PDFs f f f (12) in the compact form
d f t d t = − Θ t ( M f t − b ) , b = ( F K ) T C Y − 1 Y , (15) \frac{d f_t}{d t} = -\Theta_t\left(Mf_t-b\right), \qquad b = (\mathrm{FK})^T C_Y^{-1}Y, \tag*{(15)} d t d f t = − Θ t ( M f t − b ) , b = ( FK ) T C Y − 1 Y , (15) where b b b pull that collects the data-driven source term, and the metric M M M
M = ( F K ) T C Y − 1 ( F K ) ∈ R N v e c × N v e c . (16) M = (\mathrm{FK})^T C_Y^{-1}(\mathrm{FK}) \in\mathbb{R}^{N_{\mathrm{vec}}\times N_{\mathrm{vec}}}. \tag*{(16)} M = ( FK ) T C Y − 1 ( FK ) ∈ R N vec × N vec . (16) Since the parameters are updated by minimising the training loss alone, the data Y Y Y C Y C_Y C Y (15) are those of the training set.
The flow is thus driven by the interplay of two matrices: the NTK Θ t \Theta_t Θ t M M M Figure 1 we show the matrix M M M § 3.1 we give a statistical interpretation of the operator M M M [43 ] is discussed in § 3.2 .
Figure 1. The block matrices of the operator M M M x x x M M M
The frozen NTK and the training Hamiltonian
The NTK in eq. (13) depends on the training time t t t (15) . However, as discussed at length in Ref. [44 ] , we expect the training dynamics to exhibit two distinct regimes: an early rich phase, in which the NTK evolves rapidly and aligns its spectrum with the physically relevant directions, followed by a lazy phase, in which the NTK becomes approximately constant. In the lazy regime, the NTK can be frozen to its value at a reference time t N T K t_{\mathrm{NTK}} t NTK
Θ t = Θ t N T K ≡ Θ , ∀ t ≥ t N T K , (17) \Theta_t = \Theta_{t_{\mathrm{NTK}}} \equiv\Theta, \qquad\forall\, t \geq t_{\mathrm{NTK}}, \tag*{(17)} Θ t = Θ t NTK ≡ Θ , ∀ t ≥ t NTK , (17) so that eq. (15) reduces to a linear differential equation with constant coefficients. Note that in the infinite-width limit [47 , 48 ] the NTK is constant from initialisation, and the linearised dynamics holds along the
whole trajectory. This is not necessarily true in our context, as also previously shown in Ref. [44 ] , so that eq. (17) holds from t N T K t_{\mathrm{NTK}} t NTK
The solution to the flow equation is obtained by projecting the flow onto the eigenbasis of the frozen NTK. Being symmetric and positive semi-definite, the NTK admits the spectral decomposition
Θ = Z Λ Z T , Λ , Z ∈ R N v e c × N v e c , (18) \Theta= Z\Lambda Z^{T}, \qquad\Lambda, Z \in\mathbb{R}^{N_{\mathrm{vec}}\times N_{\mathrm{vec}}}, \tag*{(18)} Θ = Z Λ Z T , Λ , Z ∈ R N vec × N vec , (18) where Λ \Lambda Λ Z Z Z
Z a k ≡ z a ( k ) , a , k = 1 , … , N v e c . (19) Z_{ak} \equiv z_{a}^{(k)}, \qquad a,k=1,\ldots,N_{\mathrm{vec}}. \tag*{(19)} Z ak ≡ z a ( k ) , a , k = 1 , … , N vec . (19) Note that each eigenvector can have non-vanishing entries for multiple flavours at the same time, so that the label k k k Ker ( Θ ) \operatorname{Ker}(\Theta) Ker ( Θ ) d ⊥ d^{\perp} d ⊥ Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ )
Λ = [ Λ ⊥ 0 0 0 ] , Λ ⊥ ∈ R d ⊥ × d ⊥ , (20) \Lambda= \begin{bmatrix} \Lambda_{\perp} & 0 \\ 0 & 0 \end{bmatrix}, \qquad\Lambda_{\perp} \in\mathbb{R}^{d^{\perp}\times d^{\perp}}, \tag*{(20)} Λ = [ Λ ⊥ 0 0 0 ] , Λ ⊥ ∈ R d ⊥ × d ⊥ , (20) which induces a similar decomposition of the eigenvectors
Z = [ Z ⊥ Z ∥ ] , Z ⊥ ∈ R N v e c × d ⊥ , Z ∥ ∈ R N v e c × ( N v e c − d ⊥ ) , (21) Z = \begin{bmatrix} Z_{\perp} & Z_{\parallel} \end{bmatrix}, \qquad Z_{\perp} \in\mathbb{R}^{N_{\mathrm{vec}}\times d^{\perp}}, \qquad Z_{\parallel} \in\mathbb{R}^{N_{\mathrm{vec}}\times(N_{\mathrm{vec}}-d^{\perp})}, \tag*{(21)} Z = [ Z ⊥ Z ∥ ] , Z ⊥ ∈ R N vec × d ⊥ , Z ∥ ∈ R N vec × ( N vec − d ⊥ ) , (21) satisfying the following identities:
Z s T Z s = 1 s , P s = Z s Z s T ≠ 1 s , s ∈ { ∥ , ⊥ } (22) Z_{s}^{T}Z_{s} = \mathbb{1}_{s}, \qquad P_{s}=Z_{s}Z_{s}^{T}\ne\mathbb{1}_{s}, \qquad s\in\{\parallel,\perp\} \tag*{(22)} Z s T Z s = 1 s , P s = Z s Z s T = 1 s , s ∈ { ∥ , ⊥ } (22) Z ⊥ Z ⊥ T + Z ∥ Z ∥ T = P ⊥ + P ∥ = 1 N v e c , (23) Z_{\perp}Z_{\perp}^{T}+Z_{\parallel}Z_{\parallel}^{T}=P_{\perp}+P_{\parallel}=\mathbb{1}_{N_{\mathrm{vec}}}, \tag*{(23)} Z ⊥ Z ⊥ T + Z ∥ Z ∥ T = P ⊥ + P ∥ = 1 N vec , (23) where P ⊥ P_{\perp} P ⊥ P ∥ P_{\parallel} P ∥ Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) Ker ( Θ ) \operatorname{Ker}(\Theta) Ker ( Θ ) [44 ] , the vector of PDFs f f f
f ∥ = Z ∥ T f , f ⊥ = Λ ⊥ − 1 / 2 Z ⊥ T f . (24) f^{\parallel}=Z_{\parallel}^{T}f, \qquad f^{\perp}=\Lambda_{\perp}^{-1/2}Z_{\perp}^{T}f. \tag*{(24)} f ∥ = Z ∥ T f , f ⊥ = Λ ⊥ − 1/2 Z ⊥ T f . (24) By construction, the components that live in the kernel do not evolve during lazy training and remain at their initial condition, f t ∥ = f 0 ∥ f_{t}^{\parallel}=f_{0}^{\parallel} f t ∥ = f 0 ∥ training Hamiltonian
H ⊥ = Λ ⊥ 1 / 2 Z ⊥ T M Z ⊥ Λ ⊥ 1 / 2 , (25) H^{\perp}=\Lambda_{\perp}^{1/2}Z_{\perp}^{T}MZ_{\perp}\Lambda_{\perp}^{1/2}, \tag*{(25)} H ⊥ = Λ ⊥ 1/2 Z ⊥ T M Z ⊥ Λ ⊥ 1/2 , (25) which is symmetric and positive definite under the mild assumption Ker ( M ) ∩ Im ( Θ ) = { 0 } \operatorname{Ker}(M)\cap\operatorname{Im}(\Theta)=\{0\} Ker ( M ) ∩ Im ( Θ ) = { 0 }
We emphasise that eq. (25) is where the regularisation of the inverse problem is manifest. The matrix M M M M M M M M M
With this consideration in mind, we can write
H ⊥ w ( i ) = h ( i ) w ( i ) , i = 1 , … , d ⊥ , w ( i ) ∈ R d ⊥ , (26) H^{\perp}w^{(i)}=h^{(i)}w^{(i)}, \qquad i=1,\ldots,d^{\perp}, \qquad w^{(i)}\in\mathbb{R}^{d^{\perp}}, \tag*{(26)} H ⊥ w ( i ) = h ( i ) w ( i ) , i = 1 , … , d ⊥ , w ( i ) ∈ R d ⊥ , (26) and it is convenient to map the eigenvectors w ( i ) w^{(i)} w ( i )
Q = Z ⊥ Λ ⊥ 1 / 2 W , or Q a i ≡ q a ( i ) = ∑ k = 1 d ⊥ λ ( k ) z a ( k ) w k ( i ) ∈ R N v e c , (27) Q=Z_{\perp}\Lambda_{\perp}^{1/2}W, \qquad\text{or} \qquad Q_{ai}\equiv q_{a}^{(i)}=\sum_{k=1}^{d^{\perp}}\sqrt{\lambda^{(k)}}z_{a}^{(k)}w_{k}^{(i)}\in\mathbb{R}^{N_{\mathrm{vec}}}, \tag*{(27)} Q = Z ⊥ Λ ⊥ 1/2 W , or Q ai ≡ q a ( i ) = k = 1 ∑ d ⊥ λ ( k ) z a ( k ) w k ( i ) ∈ R N vec , (27) where k k k q ( i ) q^{(i)} q ( i ) 1 / h ( i ) 1/h^{(i)} 1/ h ( i )
q ( i ) T M q ( j ) = h ( i ) δ i j , ∑ i = 1 d ⊥ q ( i ) q ( i ) T = Z ⊥ Λ ⊥ Z ⊥ T = Θ . (28) q^{(i)T}Mq^{(j)}=h^{(i)}\delta_{ij}, \qquad\sum_{i=1}^{d^{\perp}}q^{(i)}q^{(i)T}=Z_{\perp}\Lambda_{\perp}Z_{\perp}^{T}=\Theta. \tag*{(28)} q ( i ) T M q ( j ) = h ( i ) δ ij , i = 1 ∑ d ⊥ q ( i ) q ( i ) T = Z ⊥ Λ ⊥ Z ⊥ T = Θ. (28) The first equality implies that the vectors q ( i ) q^{(i)} q ( i ) M M M
We finally observe that if we write q ( i ) = ∥ q ( i ) ∥ n ( i ) q^{(i)}=\lVert q^{(i)}\rVert n^{(i)} q ( i ) = ∥ q ( i ) ∥ n ( i ) n ( i ) n^{(i)} n ( i ) ∥ q ( i ) ∥ \lVert q^{(i)}\rVert ∥ q ( i ) ∥
h ( i ) = ∥ q ( i ) ∥ 2 n ( i ) T M n ( i ) , (29) h^{(i)}=\lVert q^{(i)}\rVert^{2}n^{(i)T}Mn^{(i)}, \tag*{(29)} h ( i ) = ∥ q ( i ) ∥ 2 n ( i ) T M n ( i ) , (29) which shows that each eigenvalue of the training Hamiltonian is the product of two terms. The first one is the weight in the mode expansion of the NTK, as in the second expression of eq. (28) , and can be small if the direction is poorly covered by the parametrisation. The second term is discussed more extensively in § 3.1 , and quantifies the resolution power of the data along a given direction n n n h ( i ) h^{(i)} h ( i )
The analytical solution
Integrating the flow in the two subspaces (see Ref. [44 ] for the explicit derivation) yields the central result
f t = U ( t ) f 0 + V ( t ) Y , (30) f_t=U(t)f_0+V(t)Y, \tag*{(30)} f t = U ( t ) f 0 + V ( t ) Y , (30) which expresses the trained solution at any time t t t U ( t ) U(t) U ( t ) f 0 f_0 f 0 V ( t ) V(t) V ( t ) [44 ] , the evolution operator U ( t ) U(t) U ( t )
U ( t ) = U ^ ⊥ ( t ) + U ˇ ⊥ ( t ) + U ∥ , (31) U(t)=\hat{U}^{\perp}(t)+\check{U}^{\perp}(t)+U^{\parallel}, \tag*{(31)} U ( t ) = U ^ ⊥ ( t ) + U ˇ ⊥ ( t ) + U ∥ , (31) where
U ^ ⊥ ( t ) = Z ⊥ Λ ⊥ 1 / 2 e − H ⊥ t Λ ⊥ − 1 / 2 Z ⊥ T , U ∥ = Z ∥ Z ∥ T . (32) \hat{U}^{\perp}(t)=Z_{\perp}\Lambda_{\perp}^{1/2}e^{-H^{\perp}t}\Lambda_{\perp}^{-1/2}Z_{\perp}^{T},\qquad U^{\parallel}=Z_{\parallel}Z_{\parallel}^{T}. \tag*{(32)} U ^ ⊥ ( t ) = Z ⊥ Λ ⊥ 1/2 e − H ⊥ t Λ ⊥ − 1/2 Z ⊥ T , U ∥ = Z ∥ Z ∥ T . (32) The contributions from U ˇ ⊥ ( t ) \check{U}^{\perp}(t) U ˇ ⊥ ( t ) V ( t ) V(t) V ( t )
M ( t ) = Z ⊥ Λ ⊥ 1 / 2 ( 1 − e − H ⊥ t ) ( H ⊥ ) − 1 Λ ⊥ 1 / 2 Z ⊥ T = ∑ i = 1 d ⊥ 1 − e − h ( i ) t h ( i ) q ( i ) q ( i ) T , (33) \begin{aligned}
\mathcal{M}(t)&=Z_{\perp}\Lambda_{\perp}^{1/2}\left(1-e^{-H^{\perp}t}\right)\left(H^{\perp}\right)^{-1}\Lambda_{\perp}^{1/2}Z_{\perp}^{T}\\
&=\sum_{i=1}^{d^{\perp}}\frac{1-e^{-h^{(i)}t}}{h^{(i)}}q^{(i)}q^{(i)T}, \tag*{(33)}
\end{aligned} M ( t ) = Z ⊥ Λ ⊥ 1/2 ( 1 − e − H ⊥ t ) ( H ⊥ ) − 1 Λ ⊥ 1/2 Z ⊥ T = i = 1 ∑ d ⊥ h ( i ) 1 − e − h ( i ) t q ( i ) q ( i ) T , (33) in terms of which we can write
U ˇ ⊥ ( t ) = − M ( t ) M U ∥ , V ( t ) = M ( t ) ( F K ) T C Y − 1 . (34) \check{U}^{\perp}(t)=-\mathcal{M}(t)MU^{\parallel},\qquad V(t)=\mathcal{M}(t)(FK)^{T}C_{Y}^{-1}. \tag*{(34)} U ˇ ⊥ ( t ) = − M ( t ) M U ∥ , V ( t ) = M ( t ) ( F K ) T C Y − 1 . (34) As observed in Ref. [44 ] , the operator M ( t ) \mathcal{M}(t) M ( t ) M M M t t t
M ∞ = lim t → ∞ M ( t ) = Z ⊥ ( Z ⊥ T M Z ⊥ ) − 1 Z ⊥ T , (35) \mathcal{M}_{\infty}=\lim_{t\to\infty}\mathcal{M}(t)=Z_{\perp}\left(Z_{\perp}^{T}MZ_{\perp}\right)^{-1}Z_{\perp}^{T}, \tag*{(35)} M ∞ = t → ∞ lim M ( t ) = Z ⊥ ( Z ⊥ T M Z ⊥ ) − 1 Z ⊥ T , (35) which satisfies M ∞ M M ∞ = M ∞ \mathcal{M}_{\infty}M\mathcal{M}_{\infty}=\mathcal{M}_{\infty} M ∞ M M ∞ = M ∞
The regularisation induced by the training time t t t ( 1 − e − h ( i ) t ) \left(1-e^{-h^{(i)}t}\right) ( 1 − e − h ( i ) t ) h ( i ) t ≫ 1 h^{(i)}t\gg1 h ( i ) t ≫ 1 h ( i ) t ≪ 1 h^{(i)}t\ll1 h ( i ) t ≪ 1 h ( i ) t h^{(i)}t h ( i ) t t t t [44 ] for an extensive discussion of these operators. Here we provide a second reading of the analytical expression, which ultimately facilitates our discussion in § 2.6 .
An alternative form. The three contributions to U ( t ) U(t) U ( t ) (25) , one can prove that M ( t ) M Z ⊥ Z ⊥ T = Z ⊥ Z ⊥ T − U ^ ⊥ ( t ) \mathcal{M}(t)MZ_{\perp}Z_{\perp}^{T}=Z_{\perp}Z_{\perp}^{T}-\hat{U}^{\perp}(t) M ( t ) M Z ⊥ Z ⊥ T = Z ⊥ Z ⊥ T − U ^ ⊥ ( t ) (23) yields
U ( t ) = 1 N v e c − M ( t ) M . (36) U(t)=\mathbb{1}_{N_{\mathrm{vec}}}-\mathcal{M}(t)M. \tag*{(36)} U ( t ) = 1 N vec − M ( t ) M . (36) The analytical solution in eq. (30) then takes the compact form
f t = f 0 − M ( t ) ( M f 0 − b ) , (37) f_t=f_0-\mathcal{M}(t)(Mf_0-b), \tag*{(37)} f t = f 0 − M ( t ) ( M f 0 − b ) , (37) where b b b (15) . This form offers a second reading of the solution. The combination M f 0 − b \mathcal{M}f_0-b M f 0 − b f f f f 0 f_0 f 0 M ( t ) \mathcal{M}(t) M ( t ) t t t M ( t ) ≃ t Θ \mathcal{M}(t)\simeq t\Theta M ( t ) ≃ t Θ (37) reduces to the first discrete step of the flow in eq. (15) .
Early stopping as a bias-variance trade-off
We conclude this section by showing how the analytical solution disentangles the role of early stopping in the training process. To do so, we use the L2 pseudodata of eq. (7) to recast eq. (37) into the teacher-student form, which reads
f t = f 0 + M ( t ) M ( f ( i n ) − f 0 ) + M ( t ) ( F K ) T C γ − 1 ( η + ξ ( k ) ) (38) f_t=f_0+\mathcal{M}(t)\mathcal{M}\left(f^{(\mathrm{in})}-f_0\right)+\mathcal{M}(t)(\mathrm{FK})^{T}C_{\gamma}^{-1}\left(\eta+\xi^{(k)}\right) \tag*{(38)} f t = f 0 + M ( t ) M ( f ( in ) − f 0 ) + M ( t ) ( FK ) T C γ − 1 ( η + ξ ( k ) ) (38) If we neglect the fluctuations of the NTK, and treat f 0 f_0 f 0 M \mathcal{M} M
The first term – the teacher in the language of teacher-student settings – depends on f ( i n ) f^{(\mathrm{in})} f ( in ) f 0 f_0 f 0 f ( i n ) − f 0 f^{(\mathrm{in})}-f_0 f ( in ) − f 0 (33) , we see that this is achieved by incrementally including new modes as they are resolved by the filter factor. For t → ∞ t\to\infty t → ∞ M ∞ \mathcal{M}_{\infty} M ∞ M \mathcal{M} M Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) ( F K ) f ( i n ) (\mathrm{FK})f^{(\mathrm{in})} ( FK ) f ( in ) f 0 + Im ( Θ ) f_0+\operatorname{Im}(\Theta) f 0 + Im ( Θ )
The second term, on the other hand, injects the noise of the data and carries no information about f ( i n ) f^{(\mathrm{in})} f ( in ) η \eta η ξ ( k ) \xi^{(k)} ξ ( k ) q ( i ) q^{(i)} q ( i ) (33) , the shift to the central value becomes
M ( t ) ( F K ) T C γ − 1 η = ∑ i = 1 d ⊥ 1 − e − h ( i ) t n ( i ) T M n ( i ) n ( i ) n ( i ) T ( F K ) T C γ − 1 η (39) \mathcal{M}(t)(\mathrm{FK})^{T}C_{\gamma}^{-1}\eta=\sum_{i=1}^{d_{\perp}}\frac{1-e^{-h^{(i)}t}}{n^{(i)T}\mathcal{M}n^{(i)}}n^{(i)}n^{(i)T}(\mathrm{FK})^{T}C_{\gamma}^{-1}\eta \tag*{(39)} M ( t ) ( FK ) T C γ − 1 η = i = 1 ∑ d ⊥ n ( i ) T M n ( i ) 1 − e − h ( i ) t n ( i ) n ( i ) T ( FK ) T C γ − 1 η (39) while the contribution to the variance is
Cov [ M ( t ) ( F K ) T C γ − 1 ξ ( k ) ] = M ( t ) M M ( t ) = ∑ i = 1 d ⊥ ( 1 − e − h ( i ) t ) 2 n ( i ) T M n ( i ) n ( i ) n ( i ) T (40) \operatorname{Cov}\left[\mathcal{M}(t)(\mathrm{FK})^{T}C_{\gamma}^{-1}\xi^{(k)}\right]=\mathcal{M}(t)\mathcal{M}\mathcal{M}(t)=\sum_{i=1}^{d_{\perp}}\frac{\left(1-e^{-h^{(i)}t}\right)^2}{n^{(i)T}\mathcal{M}n^{(i)}}n^{(i)}n^{(i)T} \tag*{(40)} Cov [ M ( t ) ( FK ) T C γ − 1 ξ ( k ) ] = M ( t ) MM ( t ) = i = 1 ∑ d ⊥ n ( i ) T M n ( i ) ( 1 − e − h ( i ) t ) 2 n ( i ) n ( i ) T (40) where in both expressions we used eq. (29) to express the weights in terms of the inverse of the resolution power of the data along the i i i n ( i ) T M n ( i ) n^{(i)T}\mathcal{M}n^{(i)} n ( i ) T M n ( i ) 3 t ≫ 1 / h ( i ) t\gg1/h^{(i)} t ≫ 1/ h ( i ) i i i (39) and (40) reduce to the same factor. Therefore, those directions that are poorly constrained by the data are also those with the largest impact on the central value and on the variance. By eq. (29) , we also see that these poorly constrained directions yield small h ( i ) h^{(i)} h ( i ) 4
We conclude that the inclusion of new modes, while improving the teacher term, can actually lead to a poor performance of the overall solution, as we shall explicitly see in § 4.5 (Figure 11 and Figure 13 ). Early stopping can thus be seen as a prescription to set the trade-off between the bias (the teacher plus the noise shift) and the variance (injected by ξ ( k ) \xi^{(k)} ξ ( k ) [62 ] for a simplified, solvable setting.
Figure 11. The analytical solution (blue) of the eight fitted flavours (rows) computed with the NTK frozen at τ N T K = 15000 \tau_{\mathrm{NTK}} = 15000 τ NTK = 15000 f 0 f_0 f 0 τ max = 20000 \tau_{\max} = 20000 τ m a x = 20000 f ( i n ) f^{(\mathrm{in})} f ( in )
Figure 13. Decomposition of the analytical solution of eq. (30) into the contribution from the initial condition, U ( t ) f 0 U(t)f_0 U ( t ) f 0 V ( t ) Y V(t)Y V ( t ) Y τ \tau τ Figure 11 (columns). The NTK is frozen at epoch 15000. The initial condition and the statistics over the ensemble of replicas are as in Figure 11 . The dashed black line is the underlying law f ( i n ) f^{(\mathrm{in})} f ( in )
Information content and geometry of a PDF fit
In the previous section, we introduced the two main objects of the present analysis, namely the NTK and the operator M M M § 3.1 we show that M M M § 3.2 we recover the Shapley value of Ref. [43 ] from the operator M M M § 3.3 we present a novel geometric interpretation of a PDF fit. This section is complementary to the numerical analysis presented in this work. The reader interested in the main results can jump to § 4 .
Resolving power of the data
We show that the matrix M M M (16) , is the FIM that measures the amount of information carried by the data Y Y Y f f f p ( Y ∣ f ) ≡ p f ( Y ) p(Y \mid f) \equiv p_f(Y) p ( Y ∣ f ) ≡ p f ( Y ) Y Y Y f f f
− 2 log p f ( Y ) ∝ ( Y − T [ f ] ) T C Y − 1 ( Y − T [ f ] ) (41) -2 \log p_f(Y) \propto(Y-T[f])^{T} C_Y^{-1}(Y-T[f]) \tag*{(41)} − 2 log p f ( Y ) ∝ ( Y − T [ f ] ) T C Y − 1 ( Y − T [ f ]) (41) which is the usual quadratic form that compares predictions against data in units of the experimental uncertainties. Note that eq. (41) is our deliberate choice in defining a notion of distance between two points in the space of the data, with the metric given by the inverse of the experimental covariance matrix C Y C_Y C Y (41) , in the space of PDFs.5
To do so, we use a variational approach and study how much a small perturbation δ f \delta f δ f p f ( Y ) p_f(Y) p f ( Y )
D K L ( p f ( Y ) ∥ p f + δ f ( Y ) ) ≡ ∫ Y p f ( Y ) log p f ( Y ) p f + δ f ( Y ) (42) D_{\mathrm{KL}}\left(p_f(Y)\mathbin\Vert p_{f+\delta f}(Y)\right) \equiv\int_Y p_f(Y)\log\frac{p_f(Y)}{p_{f+\delta f}(Y)} \tag*{(42)} D KL ( p f ( Y ) ∥ p f + δ f ( Y ) ) ≡ ∫ Y p f ( Y ) log p f + δ f ( Y ) p f ( Y ) (42) where p f + δ f ( Y ) p_{f+\delta f}(Y) p f + δ f ( Y ) δ f ∈ R N v e c \delta f \in\mathbb{R}^{N_{\mathrm{vec}}} δ f ∈ R N vec score function s ≡ s ( Y ; f ) ∈ R N v e c s \equiv s(Y;f) \in\mathbb{R}^{N_{\mathrm{vec}}} s ≡ s ( Y ; f ) ∈ R N vec
s i α = 1 p f ( Y ) ∂ p f ( Y ) ∂ f i α = ∂ log p f ( Y ) ∂ f i α (43) s_{i\alpha} = \frac{1}{p_f(Y)}\frac{\partial p_f(Y)}{\partial f_{i\alpha}} = \frac{\partial\log p_f(Y)}{\partial f_{i\alpha}} \tag*{(43)} s i α = p f ( Y ) 1 ∂ f i α ∂ p f ( Y ) = ∂ f i α ∂ log p f ( Y ) (43) which measures how sensitive the probability of the observed data is to a change in the PDFs, relative to the original distribution p f ( Y ) p_f(Y) p f ( Y ) (42) in powers of δ f \delta f δ f
D K L ( p f ( Y ) ∥ p f + δ f ( Y ) ) = δ f T E [ s ] + 1 2 δ f T F δ f + O ( δ f 3 ) (44) D_{\mathrm{KL}}\left(p_f(Y)\mathbin\Vert p_{f+\delta f}(Y)\right) = \delta f^{T}\mathbb{E}[s] + \frac{1}{2}\delta f^{T}\mathcal{F}\,\delta f + O\left(\delta f^{3}\right) \tag*{(44)} D KL ( p f ( Y ) ∥ p f + δ f ( Y ) ) = δ f T E [ s ] + 2 1 δ f T F δ f + O ( δ f 3 ) (44) where the expectation values are evaluated with respect to p f ( Y ) p_f(Y) p f ( Y ) F \mathcal{F} F
F = − E [ ∂ 2 log p f ( Y ) ∂ f ∂ f T ] (45) \mathcal{F}=-\mathbb{E}\left[\frac{\partial^{2}\log p_f(Y)}{\partial f\,\partial f^{T}}\right] \tag*{(45)} F = − E [ ∂ f ∂ f T ∂ 2 log p f ( Y ) ] (45) which is the expected local curvature of the log-likelihood. Assuming mild regularity conditions on p f ( Y ) p_f(Y) p f ( Y ) [68 ] , one can show that:
The score has zero mean, and therefore the first term in eq. (44) disappears;
The FIM admits an alternative form, which is F = E [ s s T ] \mathcal{F}=\mathbb{E}\left[ss^{T}\right] F = E [ s s T ]
From these two properties, it follows that the FIM is the variance of the score function: a large fluctuation of the score means that the data are very sensitive to f f f F = M \mathcal{F}=M F = M (44) can finally be written as
D K L ( p f ( Y ) ∥ p f + δ f ( Y ) ) = 1 2 δ f T M δ f + O ( δ f 3 ) (46) D_{\mathrm{KL}}\left(p_f(Y)\mathbin\Vert p_{f+\delta f}(Y)\right)=\frac{1}{2}\delta f^{T}M\,\delta f+O\left(\delta f^{3}\right) \tag*{(46)} D KL ( p f ( Y ) ∥ p f + δ f ( Y ) ) = 2 1 δ f T M δ f + O ( δ f 3 ) (46) This result shows that the matrix M M M (41) . A different prescription for comparing the data with the theory would imply a different notion of distance for the PDFs.
We emphasise that eq. (46) can be read in two complementary ways. Read forward, it is a statement on the distribution of the data given f f f M M M Ker ( M ) \operatorname{Ker}(M) Ker ( M ) M M M
The Shapley value in the DIS case
In game theory, the Shapley value is an estimator that quantifies the contribution to the total payoff of each player among a collaborating group. Adapted to the context of PDFs, it has been recently used in Ref. [43 ] to measure how much each individual flavour and kinematic region in x x x
Here we show that the quadratic form of eq. (46) is exactly the object that is computed, albeit in a different language, by the Shapley-value analysis: there the perturbations δ f \delta f δ f x x x M M M
In the referenced work, a localised Gaussian bump centred at a chosen x ^ \hat{x} x ^ S S S
f ( S ± ) = f + δ f ( S ± ) , (47) f^{(S^\pm)} = f + \delta f^{(S^\pm)}, \tag*{(47)} f ( S ± ) = f + δ f ( S ± ) , (47) where f f f (4) , and δ f ( S ± ) \delta f^{(S^\pm)} δ f ( S ± ) S S S j ∈ S j \in S j ∈ S x ^ \hat{x} x ^ [43 ]
δ f j α ( ± ) = ± Δ j ( ± ) ( x ^ ) exp [ − ( log 10 x α − log 10 x ^ ) 2 2 w 2 ] , (48) \delta f_{j\alpha}^{(\pm)} = \pm\Delta_j^{(\pm)}(\hat{x}) \exp\left[-\frac{(\log_{10} x_\alpha-\log_{10}\hat{x})^2}{2w^2}\right], \tag*{(48)} δ f j α ( ± ) = ± Δ j ( ± ) ( x ^ ) exp [ − 2 w 2 ( log 10 x α − log 10 x ^ ) 2 ] , (48) where Δ j ( ± ) ( x ^ ) \Delta_j^{(\pm)}(\hat{x}) Δ j ( ± ) ( x ^ ) x ^ \hat{x} x ^ S S S χ 2 \chi^2 χ 2 D \mathcal{D} D v ( S ) v(S) v ( S )
v ( S ) = 1 2 [ χ 2 ( f ( S + ) ; D ) + χ 2 ( f ( S − ) ; D ) ] , (49) v(S) = \frac{1}{2}\left[\chi^2\left(f^{(S^+)};\mathcal{D}\right) + \chi^2\left(f^{(S^-)};\mathcal{D}\right)\right], \tag*{(49)} v ( S ) = 2 1 [ χ 2 ( f ( S + ) ; D ) + χ 2 ( f ( S − ) ; D ) ] , (49) where χ 2 ( f ( S ± ) ; D ) \chi^2\left(f^{(S^\pm)};\mathcal{D}\right) χ 2 ( f ( S ± ) ; D ) S S S D \mathcal{D} D δ T ( S ± ) \delta T^{(S^\pm)} δ T ( S ± ) δ f ( S ± ) \delta f^{(S^\pm)} δ f ( S ± ) (49) can be decomposed into four components
v ( S ) = χ 2 ( f ; D ) + ( Y − T ) T C Y − 1 ( δ T ( S + ) + δ T ( S − ) ) + 1 2 [ ( δ T ( S + ) ) T C Y − 1 δ T ( S + ) + ( δ T ( S − ) ) T C Y − 1 δ T ( S − ) ] , (50) \begin{aligned}
v(S) ={}& \chi^2(f;\mathcal{D}) + (Y-T)^T C_Y^{-1}\left(\delta T^{(S^+)}+\delta T^{(S^-)}\right) \\
&+ \frac{1}{2}\left[\left(\delta T^{(S^+)}\right)^T C_Y^{-1}\delta T^{(S^+)} + \left(\delta T^{(S^-)}\right)^T C_Y^{-1}\delta T^{(S^-)}\right], \tag*{(50)}
\end{aligned} v ( S ) = χ 2 ( f ; D ) + ( Y − T ) T C Y − 1 ( δ T ( S + ) + δ T ( S − ) ) + 2 1 [ ( δ T ( S + ) ) T C Y − 1 δ T ( S + ) + ( δ T ( S − ) ) T C Y − 1 δ T ( S − ) ] , (50) where T T T (2) , this expression simplifies drastically if we further assume symmetric bounds in eq. (48) , which implies δ f ( S ) ≡ δ f ( S + ) = − δ f ( S − ) \delta f^{(S)} \equiv\delta f^{(S^+)} = -\delta f^{(S^-)} δ f ( S ) ≡ δ f ( S + ) = − δ f ( S − )
v ( S ) = χ 2 ( f ; D ) + ( δ T ( S ) ) T C Y − 1 δ T ( S ) = χ 2 ( f ; D ) + ( δ f ( S ) ) T M δ f ( S ) , (51) v(S) = \chi^2(f;\mathcal{D}) + \left(\delta T^{(S)}\right)^T C_Y^{-1}\delta T^{(S)} = \chi^2(f;\mathcal{D}) + \left(\delta f^{(S)}\right)^T M\,\delta f^{(S)}, \tag*{(51)} v ( S ) = χ 2 ( f ; D ) + ( δ T ( S ) ) T C Y − 1 δ T ( S ) = χ 2 ( f ; D ) + ( δ f ( S ) ) T M δ f ( S ) , (51) where in the second equality we used eq. (2) to express the second term in PDF space. For a given coalition S S S i i i
ϕ i ( S ) = v ( S ∪ { i } ) − v ( S ) , (52) \phi_i(S) = v(S \cup\{i\}) - v(S), \tag*{(52)} ϕ i ( S ) = v ( S ∪ { i }) − v ( S ) , (52) that is, the difference in the value functions with and without i i i Φ i ( x ^ ) \Phi_i(\hat{x}) Φ i ( x ^ ) 2 N f l a v − 1 2^{N_{\mathrm{flav}}-1} 2 N flav − 1 B ∖ { i } \mathcal{B} \setminus\{i\} B ∖ { i } N f l a v − 1 N_{\mathrm{flav}} - 1 N flav − 1 6
Φ i = ∑ S ⊆ B ∖ { i } w ( S ) ϕ i ( S ) , (53) \Phi_i = \sum_{S \subseteq\mathcal{B}\setminus\{i\}} w(S)\,\phi_i(S), \tag*{(53)} Φ i = S ⊆ B ∖ { i } ∑ w ( S ) ϕ i ( S ) , (53) where w ( S ) w(S) w ( S ) S S S
w ( S ) = ∣ S ∣ ! ( N f l a v − ∣ S ∣ − 1 ) ! N f l a v ! , (54) w(S) = \frac{|S|!\left(N_{\mathrm{flav}}-|S|-1\right)!}{N_{\mathrm{flav}}!}, \tag*{(54)} w ( S ) = N flav ! ∣ S ∣ ! ( N flav − ∣ S ∣ − 1 ) ! , (54) where ∣ S ∣ |S| ∣ S ∣
Recall that M M M M i j M_{ij} M ij Figure 1 . Since δ f ( S ) \delta f(S) δ f ( S ) S S S x x x x ^ \hat{x} x ^ M M M S S S (53) reduces to the following expression
Φ i = δ f i T M i i δ f i + ∑ j ≠ i δ f i T M i j δ f j , (55) \Phi_i = \delta f_i^T M_{ii}\delta f_i + \sum_{j\ne i}\delta f_i^T M_{ij}\delta f_j, \tag*{(55)} Φ i = δ f i T M ii δ f i + j = i ∑ δ f i T M ij δ f j , (55) which is the sum of quadratic forms that measure the actual contribution of the perturbation as seen by the FK tables and weighted by the experimental uncertainty. The diagonal entries M i i M_{ii} M ii M i j M_{ij} M ij M M M
A geometric interpretation of training
We now propose a geometric interpretation of the training dynamics and of the analytical solution of eq. (30) . Information Geometry (IG) [69 , 70 ] has long been proposed as a theoretical framework that aims at describing statistical models by using the language of differential geometry. For instance, in Ref. [60 ] IG is used to provide generalisation bounds in ReLU networks that are invariant under reparametrisations. Natural gradient descent [71 ] in machine learning follows from geometric considerations applied to statistical manifolds. Finally, in collider physics IG has been used to analyse the latent geometry of jet taggers [72 ] and to construct bases of correlator observables [73 ] . Yet, to the best of our knowledge, applications in the context of PDF determinations are still missing. Here we do not aim at a rigorous treatment of the geometric structure of a parametric fit to PDFs, which we leave to future work, but rather focus on a more intuitive visualisation of these ideas.
The role of the model
The functional model that parametrises the PDFs is a forward map from the parameter space P ⊆ R N p a r a m s \mathcal{P} \subseteq\mathbb{R}^{N_{\mathrm{params}}} P ⊆ R N params R N v e c \mathbb{R}^{N_{\mathrm{vec}}} R N vec M f \mathcal{M}_f M f R N v e c \mathbb{R}^{N_{\mathrm{vec}}} R N vec f f f M f \mathcal{M}_f M f x x x T f M T_f\mathcal{M} T f M i.e. by the tangent vectors
J μ = ∂ μ f ∈ R N v e c , μ = 1 , … , N p a r a m s , (56) J_\mu= \partial_\mu f \in\mathbb{R}^{N_{\mathrm{vec}}}, \qquad\mu= 1,\ldots,N_{\mathrm{params}}, \tag*{(56)} J μ = ∂ μ f ∈ R N vec , μ = 1 , … , N params , (56) each encoding the sensitivity of the function to a displacement of the corresponding parameter. Note that the tangent vectors are not necessarily orthogonal, nor linearly independent. We use the word frame to denote this overcomplete set of vectors that live in the tangent plane. In the NN parametrisation of eq. (5) , N p a r a m s = 763 N_{\mathrm{params}} = 763 N params = 763 N v e c = N g r i d ⋅ N f l a v = 700 N_{\mathrm{vec}} = N_{\mathrm{grid}} \cdot N_{\mathrm{flav}} = 700 N vec = N grid ⋅ N flav = 700 M f \mathcal{M}_f M f intrinsic geometric objects, and therefore invariant under reparametrisations.
The choice of the parametric model not only determines the manifold M f \mathcal{M}_f M f M f \mathcal{M}_f M f
Example. Consider the forward map f ( x ; θ 1 , θ 2 ) = θ 1 θ 2 x f(x;\theta_1,\theta_2) = \theta_1\theta_2 x f ( x ; θ 1 , θ 2 ) = θ 1 θ 2 x J 1 = θ 2 x J_1 = \theta_2 x J 1 = θ 2 x J 2 = θ 1 x J_2 = \theta_1 x J 2 = θ 1 x x x x span { x } \operatorname{span}\{x\} span { x }
δ f ( x ) = ∂ 1 f ( x ) δ θ 1 + ∂ 2 f ( x ) δ θ 2 = ( δ θ 1 θ 2 + θ 1 δ θ 2 ) x . (57) \delta f(x) = \partial_1 f(x)\,\delta\theta_1 + \partial_2 f(x)\,\delta\theta_2 = (\delta\theta_1\theta_2 + \theta_1\delta\theta_2)\,x . \tag*{(57)} δ f ( x ) = ∂ 1 f ( x ) δ θ 1 + ∂ 2 f ( x ) δ θ 2 = ( δ θ 1 θ 2 + θ 1 δ θ 2 ) x . (57) Any displacement along the direction ( δ θ 1 , δ θ 2 ) ∝ ( θ 1 , − θ 2 ) (\delta\theta_1,\delta\theta_2) \propto(\theta_1,-\theta_2) ( δ θ 1 , δ θ 2 ) ∝ ( θ 1 , − θ 2 ) ( θ 1 , θ 2 ) → ( c θ 1 , θ 2 / c ) (\theta_1,\theta_2) \to(c\theta_1,\theta_2/c) ( θ 1 , θ 2 ) → ( c θ 1 , θ 2 / c ) f ( x ; θ ) = θ x f(x;\theta) = \theta x f ( x ; θ ) = θ x x x x J = x J = x J = x
Θ = ( θ 1 2 + θ 2 2 ) x x T and Θ = x x T , (58) \Theta= (\theta_1^2 + \theta_2^2)\,x x^T \qquad\text{and} \qquad\Theta= x x^T , \tag*{(58)} Θ = ( θ 1 2 + θ 2 2 ) x x T and Θ = x x T , (58) respectively. These have the same eigenvector, but different eigenvalues. In fact, in the first case the strength with which the direction x x x
From this viewpoint, the NTK is nothing but the Gram matrix of the tangent vectors,
Θ = ∑ μ J μ J μ T , (59) \Theta= \sum_\mu J_\mu J_\mu^T , \tag*{(59)} Θ = μ ∑ J μ J μ T , (59) and as such it encodes the covering capability of the frame. We can understand how tangent vectors make up the eigenvectors of the NTK by taking the quadratic form (also known as the Rayleigh quotient) with any unit direction n n n
n T Θ n = ∑ μ ( J μ , n ) 2 , (60) n^T\Theta n = \sum_\mu(J_\mu,n)^2 , \tag*{(60)} n T Θ n = μ ∑ ( J μ , n ) 2 , (60) where the spectrum is given by the stationary points of eq. (60) . Therefore, if n n n Θ \Theta Θ J μ J_\mu J μ n n n n n n θ \theta θ
Training as a flow on the manifold
A gradient flow on the manifold of parameters M θ \mathcal{M}_\theta M θ
θ ˙ = − G − 1 ∇ θ L , (61) \dot{\theta} = -G^{-1}\nabla_\theta\mathcal{L}, \tag*{(61)} θ ˙ = − G − 1 ∇ θ L , (61) where G G G M θ \mathcal{M}_\theta M θ ∇ θ \nabla_\theta ∇ θ (15) , is thus the push-forward of eq. (61) through the Jacobian J J J
f ˙ = J θ ˙ = − J G − 1 J T ∇ f L , J ∈ R N v e c × N p a r a m s , (62) \dot{f} = J\dot{\theta} = -JG^{-1}J^T\nabla_f\mathcal{L}, \qquad J \in\mathbb{R}^{N_{\mathrm{vec}}\times N_{\mathrm{params}}} , \tag*{(62)} f ˙ = J θ ˙ = − J G − 1 J T ∇ f L , J ∈ R N vec × N params , (62) where ∇ θ L = J T ∇ f L \nabla_{\theta}\mathcal{L} = J^{T}\nabla_{f}\mathcal{L} ∇ θ L = J T ∇ f L ∇ f L = M f − b \nabla_{f}\mathcal{L} = Mf - b ∇ f L = M f − b (15) and (59) then follow from choosing the metric in eq. (11) to be the Euclidean one, G = 1 G = 1 G = 1 7 G G G (62) becomes
f ˙ = − ∑ μ J μ ( J μ , ∇ f L ) , (63) \dot{f} = -\sum_{\mu} J_{\mu}\left(J_{\mu}, \nabla_{f}\mathcal{L}\right), \tag*{(63)} f ˙ = − μ ∑ J μ ( J μ , ∇ f L ) , (63) from which we conclude that each vector J μ J_{\mu} J μ G G G G G G (62) is known as natural gradient flow [71 ] .
The velocity of the flow along any direction n n n f f f (62) , along n n n
n T f ˙ = n T Θ ∇ f L ≤ n T Θ n ∇ f L T Θ ∇ f L , (64) n^{T}\dot{f} = n^{T}\Theta\nabla_{f}\mathcal{L} \le\sqrt{n^{T}\Theta n}\sqrt{\nabla_{f}\mathcal{L}^{T}\Theta\nabla_{f}\mathcal{L}}, \tag*{(64)} n T f ˙ = n T Θ ∇ f L ≤ n T Θ n ∇ f L T Θ ∇ f L , (64) where the right-hand side factorises into two contributions. The first term is the coverage of eq. (60) , which measures the alignment between n n n (25) . If the latter is small, then the data constrain directions that are not accessible to the model, and therefore the rate of change of f f f n n n n n n n n n (60) , results in a larger bound indicating that this direction can be learned fast. Finally, if n ∈ Ker ( Θ ) n \in\operatorname{Ker}(\Theta) n ∈ Ker ( Θ )
Example. Consider the one-dimensional manifold of linear functions f ( x ) = a x f(x) = ax f ( x ) = a x a = θ a = \theta a = θ a = θ 3 a = \theta^{3} a = θ 3 a = θ 1 θ 2 a = \theta_{1}\theta_{2} a = θ 1 θ 2 span { x } \operatorname{span}\{x\} span { x } x x x Θ = x x T \Theta= xx^{T} Θ = x x T 3 θ 2 x 3\theta^{2}x 3 θ 2 x Θ = 9 θ 4 x x T \Theta= 9\theta^{4}xx^{T} Θ = 9 θ 4 x x T θ = 0 \theta= 0 θ = 0 { θ 2 x , θ 1 x } \{\theta_{2}x,\theta_{1}x\} { θ 2 x , θ 1 x } Θ = ( θ 1 2 + θ 2 2 ) x x T \Theta= (\theta_{1}^{2}+\theta_{2}^{2})xx^{T} Θ = ( θ 1 2 + θ 2 2 ) x x T
The frozen NTK dynamics as a tangent plane
Freezing the NTK at some epoch is thus equivalent to selecting the tangent plane – and its frame – at one point along the trajectory. From this point onwards, the analytical solution in eq. (37) then continues the training on this plane rather than on the manifold itself. If, over the region explored by the fit, the manifold can be approximated with its tangent plane, the approximation is faithful. Otherwise, the curvature of the manifold introduces an error that scales with the learning rate. This is the geometric visualisation of the core idea underlying the NTK dynamics first introduced in Refs. [47 , 48 ] .
Note that in the next section we freeze the NTK at some point along the trajectory, but then evolve the flow from an initial condition f 0 f_{0} f 0 (37) , can help us visualise this geometric picture. Since M ( t ) \mathcal{M}(t) M ( t ) Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) f t − f 0 f_{t}-f_{0} f t − f 0 Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) f 0 f_{0} f 0
f 0 + Im ( Θ ) , (65) f_{0} + \operatorname{Im}(\Theta), \tag*{(65)} f 0 + Im ( Θ ) , (65) Figure 2. Schematic representation of the training dynamics in the functional manifold M f \mathcal{M}_f M f M f \mathcal{M}_f M f t N T K t_{\mathrm{NTK}} t NTK T f t N T K M T_{f_{t_{\mathrm{NTK}}}\mathcal{M}} T f t NTK M f 0 f_0 f 0 f 0 + Im ( Θ ) f_0+\operatorname{Im}(\Theta) f 0 + Im ( Θ ) Ker ( Θ ) \operatorname{Ker}(\Theta) Ker ( Θ )
which is the plane through the initial condition f 0 f_0 f 0 T f t N T K M T_{f_{t_{\mathrm{NTK}}}\mathcal{M}} T f t NTK M f 0 f_0 f 0 ∥ P ∥ ( f 0 − f t N T K ) ∥ \lVert P_{\parallel}(f_0-f_{t_{\mathrm{NTK}}})\rVert ∥ P ∥ ( f 0 − f t NTK )∥ (37) , which is driven by the gradient of the loss and by the pull of the data at the frozen point, traces the trajectory on the affine space.
The physical directions
We emphasise that what really drives the analytic solution of a PDF fit is the training Hamiltonian H ⊥ H^\perp H ⊥ (25) . In fact, H ⊥ H^\perp H ⊥ (60) , but with the Euclidean metric replaced by M M M 8 vice versa . The training Hamiltonian is thus the combination of two interplaying components: on the one hand the spectrum of the NTK dictating what the parametrisation sees, on the other the matrix M M M
The NTK analysis of a DIS-only fit
The methodology outlined in § 2 forms our baseline for the NTK analysis of a DIS-only fit. We begin this section with a brief discussion of the distribution of PDFs at initialisation in § 4.1 , followed by a comprehensive study of the training dynamics in terms of the NTK, including the evolution of its spectrum and the corresponding eigenvectors, in § 4.2 . Finally, in § 4.5 we compare the analytical solution of the training dynamics with the empirical training of the network.
The model at initialisation
The probability distribution at initialisation of the model parameters – weights and biases of the NN and the preprocessing factors of the PDFs – induces a corresponding distribution in the functional space of the PDFs, which we show in Figure 3 and Figure 4 . In the eight panels, one for each fitted flavour, we show the mean (solid line), variance (shaded bands), and 68% c.l. (dashed lines) of the output at initialisation across all values in x x x N r e p = 120 N_{\rm rep} = 120 N rep = 120 x x x
Figure 3. The output of the ensemble of models at initialisation at the initial scale Q 0 Q_0 Q 0 Σ \Sigma Σ g g g V V V V 3 V_3 V 3
Figure 4. Continuation of Figure 3 with V 8 V_8 V 8 T 3 T_3 T 3 T 8 T_8 T 8 T 15 T_{15} T 15
First, we note that when the preprocessing factors are removed altogether, the distributions at initialisation appear very similar irrespective of the flavour. This is indeed expected, as the NN alone treats all flavours on the same footing at initialisation. Furthermore, the inclusion of log x \log x log x x x x [44 ] . On the other hand, the preprocessing factors act as a flavour-dependent deformation of the shapes of the PDFs and, with the exception of g g g Σ \Sigma Σ x x x
Note that in a Bayesian approach to the solution of the inverse problem, the prior distribution of the parameters at initialisation induces a corresponding prior distribution in the functional space of the PDFs. The prior distribution at initialisation is therefore a crucial ingredient in the Bayesian inference, as it encodes our prior knowledge and assumptions about the PDFs before any data are taken into account. In particular, the choice of the prior distribution can have a significant impact on the posterior distribution of the PDFs, especially in regions where data are scarce or absent. The analytical solution discussed in § 4.5 allows us to study the impact of the prior on the final result.
Diagnostic of the NTK during training
Our analysis of the training dynamics starts by inspecting the training (blue) and validation (orange) loss curves along the training trajectory, as shown in Figure 5 . Here we use the methodology outlined in § 2 to train the NNs, and we recall that momentum and valence sum rules, as well as positivity and integrability constraints, are not imposed in any of the fits here considered.
Figure 5. The evolution of the loss function during training for both the training (blue) and validation (orange) sets of data. The loss is the first term of eq. (10) , normalised to the number of data points of the corresponding set. The lines represent the mean value computed over an ensemble of 120 replicas, while the shaded bands represent the one-sigma uncertainty computed as the variance over the same ensemble. Note the logarithmic scale in the horizontal axis.
The behaviour of the two loss functions is significantly different from the one observed in the previous study [44 ] , where the loss function was characterised by a steady plateau followed by a sudden drop. In contrast, here the loss functions decay smoothly, with a rapid decrease in the first few epochs followed by a slower variation until reaching a plateau. In the following, we show how this behaviour is reflected in the internal structure of the NTK. We then devote § 5 to exploring how the training dynamics is affected by the initialisation prescription, recovering the behaviour originally observed in Ref. [44 ] .
Spectrum of the NTK
In Figure 6 we show the first 200 eigenvalues of the NTK as functions of the training time, distributed evenly across the four rows. The left and right columns display the eigenvalues of the NTK computed without and with the preprocessing factors, respectively.9 Figure 7 and Figure 8 , the first showing a close-up view of the leading ten eigenvalues where the 68% c.l. is also shown, the second tracking the dimension of the NTK, with and without preprocessing factors, over the entire training process.10
Figure 6. Eigenvalues of the NTK during training for both cases without (left column) and with (right column) preprocessing factors. The figure displays the first 200 ranks of the NTK spectrum, distributed evenly across the four rows. The colour of each curve encodes the rank i i i i.e. the numerical-zero relative threshold ε t o l λ 1 ( t ) \varepsilon_{\mathrm{tol}} \lambda_1(t) ε tol λ 1 ( t ) ε t o l = 10 − 12 \varepsilon_{\mathrm{tol}} = 10^{-12} ε tol = 1 0 − 12 Figure 8 is precisely the number of eigenvalues lying above that line at each training time.
Figure 7. Close-up view of the leading ten eigenvalues of the NTK during training for both cases without (left column) and with (right column) preprocessing factors. As in Figure 6 , the colour of each curve encodes the rank i i i
Figure 8. Dimension of the image of the NTK during training as a function of the training time. We show both the case of the full model with preprocessing factors (blue) and the case of the NN alone (orange). The solid line represents the mean value computed over an ensemble of 120 replicas, while the shaded bands represent the one-sigma uncertainty computed as the variance over the same ensemble. At each training time, this dimension is the number of eigenvalues lying above the image/kernel threshold ε t o l λ 1 \varepsilon_{\mathrm{tol}} \lambda_1 ε tol λ 1 Figure 6 .
From inspecting these figures, we observe that the main results of Ref. [44 ] remain qualitatively true in a DIS-only fit: the spectrum of the NTK is characterised by a hierarchical structure that spans several orders of magnitude and, as training progresses, the NN discovers new modes represented by eigenvalues that cross the kernel threshold and enter the image of the NTK. Unlike the single-flavour toy model, however, such behaviour is localised at the beginning of the training, which is reminiscent of the rapid decrease of the loss function observed in Figure 5 . After this discovery phase, which we refer to as the rich regime [44 , 74 ] , the spectrum of the NTK stabilises with the leading eigenvalues remaining approximately constant and the lower ones increasing slowly. We identify this latter phase as the lazy regime , where we assume the NTK to be approximately constant and the training dynamics to be well described by the analytical solution, as discussed in § 4.5 .
The comparison between the two columns of Figure 6 reveals that the inclusion of the preprocessing factors does not affect the qualitative behaviour of the spectrum, although it contributes to the magnitude of the eigenvalues. In fact, a closer inspection of Figure 7 reveals that when preprocessing factors are excluded, the first ten eigenvalues are all confined within one order of magnitude and the gap between consecutive curves is reduced. On the other hand, upon inclusion of the preprocessing factors, the eigenvalues span about three orders of magnitude, with increasing spacing between consecutive curves. Note also that the inclusion of the preprocessing factors increases the spread over replicas too, as shown by the 68% c.l. bands.
The same pattern is also reflected in the dimension of Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) Figure 8 . Here we see that the number of modes that populate Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ )
Modes of the training Hamiltonian
We now turn to discussing the modes q ( i ) q^{(i)} q ( i ) (27) , whose role in this context is twofold. First, they expose the functional forms that the parametrisation can represent when solving the inverse problem. Second, they shape the analytical solution through the mode expansion of eq. (37) . We show their components for each of the eight fitted flavours in Figure 9 for a representative subset of ranks, here chosen to be r = 1 , 2 , 3 , 5 , 10 r = 1, 2, 3, 5, 10 r = 1 , 2 , 3 , 5 , 10 τ = 0 , 1000 \tau= 0, 1000 τ = 0 , 1000 10000 10000 10000
Figure 9. The modes q ( r ) q^{(r)} q ( r ) (27) , as functions of x x x τ = 0 \tau= 0 τ = 0 τ = 1000 \tau= 1000 τ = 1000 τ = 10000 \tau= 10000 τ = 10000 r = 1 , 2 , 3 , 5 , 10 r = 1, 2, 3, 5, 10 r = 1 , 2 , 3 , 5 , 10 f ( i n ) f^{(\mathrm{in})} f ( in ) f ( i n ) f^{(\mathrm{in})} f ( in )
Starting from the upper block, we see that the modes exhibit a non-trivial structure already at initialisation, before any optimisation step is taken. An inspection by rank reveals some associations between ranks and flavours. For instance, the first mode makes up most of Σ \Sigma Σ g g g T 3 T_3 T 3 g g g V 3 V_3 V 3 f ( i n ) f^{(\mathrm{in})} f ( in )
We do not observe significant patterns emerging when looking across subsequent epochs, but rather a reorganisation of the flavours covered by each rank. However, we note that the first mode remains tied to Σ \Sigma Σ
In conclusion, the inspection of these modes is a useful tool to gain insight into the kind of primitive functions the parametrisation relies on. We remark that these modes depend on the parametrisation and,
as discussed extensively in § 5 , they change drastically if a different initialisation is adopted.
Analytical training
The analytical solution reviewed in § 2.2 can be seen as a linear combination of the Hamiltonian modes q ( i ) q^{(i)} q ( i )
τ ( i ) = 1 2 η h ( i ) , (66) \tau^{(i)} = \frac{1}{2\eta h^{(i)}}, \tag*{(66)} τ ( i ) = 2 η h ( i ) 1 , (66) where h ( i ) h^{(i)} h ( i ) H ⊥ H^{\perp} H ⊥ η \eta η 11 τ ( i ) \tau^{(i)} τ ( i ) η \eta η
Freezing the NTK at a specific epoch τ N T K \tau_{\mathrm{NTK}} τ NTK learning timescales . In the remainder of this paper we work with the reference epoch τ N T K = 15000 \tau_{\mathrm{NTK}} = 15000 τ NTK = 15000 Figure 10 as medians over the 120 replicas (leaving about 10 6 10^{6} 1 0 6 h ( i ) h^{(i)} h ( i ) τ = 20000 \tau= 20000 τ = 20000 training horizon . Everything above this line (grey shaded area) is excluded from the expansion.
Figure 10. Characteristic resolution times τ ( i ) = 1 / ( 2 η h ( i ) ) \tau^{(i)} = 1/(2\eta h^{(i)}) τ ( i ) = 1/ ( 2 η h ( i ) ) i i i Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) h ( i ) h^{(i)} h ( i ) ϵ t o l = 10 − 12 \epsilon_{\mathrm{tol}} = 10^{-12} ϵ tol = 1 0 − 12 Figure 8 , giving d ⊥ ≃ 10 6 d_{\perp} \simeq10^{6} d ⊥ ≃ 1 0 6 τ max = 20000 \tau_{\max} = 20000 τ m a x = 20000 τ \tau τ Figure 11 . Points are medians over the ensemble of replicas.
The hierarchical structure is not surprising, as it is a direct consequence of the spectrum of the NTK shown in Figure 6 . However, we observe that the points cover a significant range of timescales, spanning about 14 orders of magnitude. For instance, if the analytical solution is evolved up to the training horizon, we retain less than a third of the available modes. This is another point of contrast with the previous study [44 ] , where within the training horizon the analytical solution encompassed all the available modes. Such a difference is however expected, as it reflects the increased complexity of the current setting. We now see how these modes shape the PDFs in the analytical solution.
Convergence of the analytical solution
In order to quantify the amount of information encoded in the NTK at τ N T K \tau_{\mathrm{NTK}} τ NTK evolved – according to eq. (30) and compared against the empirically trained network. With the help of Figure 10 , we select three representative evolution times which are τ = 50 \tau= 50 τ = 50 Figure 10 , and on average they amount to 8, 17, and 30 resolved modes, respectively. The resulting PDFs are collected in Figure 11 , where the analytical solution is computed for each of the eight fitted flavours (rows) at the three representative evolution times (columns). In addition, we show two additional columns where the analytical evolution is computed at τ = 10 5 \tau= 10^{5} τ = 1 0 5 τ = 10 7 \tau= 10^{7} τ = 1 0 7 τ max = 20000 \tau_{\max} = 20000 τ m a x = 20000
The agreement between the analytical and trained solutions is generally good across all flavours and for the first three evolution times. Moreover, the two solutions converge onto the underlying law at the same pace and, at τ = 20000 \tau= 20000 τ = 20000
These results are consistent with the slowly varying trajectories of the leading NTK eigenvalues (Figure 6 ) and the monotonic, rapid decay of the loss function (Figure 5 ), the latter well approximated by the exponential behaviour of the analytical solution, eq. (30) . This is significantly different from our previous work [44 ] , where the analytical solution converged to the underlying law much faster than the trained network. This, however, was due to the inconsistency between the smooth behaviour of the loss function induced by the analytical solution and the abrupt learning phase that characterised the trained solution. In § 5 we explore what happens to the validity of the analytical solution when this sudden learning phase is induced in the training dynamics.
The only notable exception across the eight flavours is the small-x x x T 3 T_{3} T 3 x x x
We show in the last two columns of Figure 11 that this deviation does not improve upon inclusion of additional modes in the analytical solution. At τ = 10 5 \tau= 10^{5} τ = 1 0 5 x ∼ 10 − 2 x \sim10^{-2} x ∼ 1 0 − 2 τ = 10 7 \tau= 10^{7} τ = 1 0 7 § 2.6 , and demonstrates that the modes learned at later stages during training – those further down in the spectrum – are poorly constrained by the data. Their inclusion thus introduces noise in both central values and uncertainties with a noise that depends on the matrix M M M
Impact of the choice of τ N T K \tau_{\mathrm{NTK}} τ NTK
One might wonder how the choice of τ N T K \tau_{\mathrm{NTK}} τ NTK Figure 12 , where we show the analytical solution of the singlet Σ \Sigma Σ τ N T K \tau_{\mathrm{NTK}} τ NTK τ N T K = 0 \tau_{\mathrm{NTK}} = 0 τ NTK = 0 τ \tau τ Figure 11 , represented by the columns. The underlying solution is also shown.
Figure 12. Analytical (blue) versus empirically trained (orange) evolution of the singlet Σ \Sigma Σ τ N T K = 0 , 5000 , 10000 \tau_{\mathrm{NTK}} = 0, 5000, 10000 τ NTK = 0 , 5000 , 10000 τ \tau τ Figure 11 (columns). Statistics and the underlying law (dashed black) are as in Figure 11 .
Surprisingly, the analytical solution performs well even when the NTK is frozen at initialisation, as it is capable of tracking the trained network and reconstructing the underlying law. This was already anticipated when discussing the modes q ( i ) q^{(i)} q ( i ) Figure 9 , where we observed that these have non-trivial shapes even at initialisation. We remark that this result is different from the behaviour observed in our previous study [44 ] , and more similar to the notion of lazy training of Refs. [47 , 48 ] . The mild deviations that are observed are mostly ascribed to the lack of sub-leading eigenvalues that are not yet discovered at initialisation (see Figure 6 ). In fact, the agreement with the trained solution improves as τ N T K \tau_{\mathrm{NTK}} τ NTK Figure 10 ). Later in § 5 we show to what extent this behaviour is tied to the choice of initialisation.
The initial condition and the data
The analytical solution in eq. (30) splits into two contributions: the first one accounting for the initial condition, U ( t ) f 0 U(t)f_{0} U ( t ) f 0 V ( t ) Y V(t)Y V ( t ) Y f U f_{U} f U f V f_{V} f V
stochastic quantities, in that the evolution operators are built from the NTK of each replica, f 0 f_0 f 0 Y Y Y [44 ] ,
Cov [ f t , f t T ] = Cov [ f U , f U T ] ⏟ C t ( 00 ) + Cov [ f V , f V T ] ⏟ C t ( Y Y ) + Cov [ f U , f V T ] + Cov [ f V , f U T ] ⏟ C t ( 0 Y ) , (67) \operatorname{Cov}[f_t, f_t^T] = \underbrace{\operatorname{Cov}[f_U, f_U^T]}_{C_t^{(00)}} + \underbrace{\operatorname{Cov}[f_V, f_V^T]}_{C_t^{(YY)}} + \underbrace{\operatorname{Cov}[f_U, f_V^T] + \operatorname{Cov}[f_V, f_U^T]}_{C_t^{(0Y)}}, \tag*{(67)} Cov [ f t , f t T ] = C t ( 00 ) Cov [ f U , f U T ] + C t ( Y Y ) Cov [ f V , f V T ] + C t ( 0 Y ) Cov [ f U , f V T ] + Cov [ f V , f U T ] , (67) namely the uncertainty inherited from the initial condition, the one propagated from the data, and their correlation. In Figure 13 we show f U f_U f U f V f_V f V f ( i n ) f^{(\mathrm{in})} f ( in ) Figure 14 we show the diagonal components of the three terms of eq. (67) as fractions of the total variance. The layout, the reference τ N T K \tau_{\mathrm{NTK}} τ NTK τ \tau τ Figure 11 .
Figure 14. Fractions of the variance of the analytical solution carried by the initial condition (green), by the data (purple), and by their correlation (dotted grey). Each curve is the ratio between the diagonal entries of the corresponding terms in eq. (67) and the total covariance, so that all fractions sum to one at each point in x x x Figure 11 . The two horizontal grey lines mark zero and one of the y y y x = 1 x = 1 x = 1
A first inspection of Figure 13 reveals that f V f_V f V x x x f U f_U f U f V f_V f V Σ \Sigma Σ τ = 50 \tau= 50 τ = 50
The relative weight of the uncertainties of the two contributions – and their correlation – is better appreciated in Figure 14 , and it varies across panels. At τ = 50 \tau= 50 τ = 50 x x x V V V V 3 V_3 V 3 V 8 V_8 V 8 T 3 T_3 T 3 T 8 T_8 T 8 f U f_U f U f V f_V f V Σ \Sigma Σ T 15 T_{15} T 15 f U f_U f U x ≲ 10 − 4 x \lesssim10^{-4} x ≲ 1 0 − 4 x ≳ 0.5 x \gtrsim0.5 x ≳ 0.5 Σ \Sigma Σ x ≲ 10 − 5 x \lesssim10^{-5} x ≲ 1 0 − 5 g g g f U f_U f U x x x f V f_V f V
As training proceeds, the analytical solution suppresses the contribution of f U f_U f U
the data, and by τ = 20000 \tau= 20000 τ = 20000 x ∼ 10 − 1 x \sim10^{-1} x ∼ 1 0 − 1 Σ \Sigma Σ T 3 T_3 T 3 x x x Figure 13 . In the latter, the contribution of f U f_U f U x ≲ 10 − 1 x \lesssim10^{-1} x ≲ 1 0 − 1 Figure 13 .
We already observed in Figure 11 that pushing the evolution beyond the training horizon does not help reduce the uncertainties of T 3 T_3 T 3 Figure 13 and Figure 14 we now see that much of the worsening previously noted is not due to an amplification of the initial condition, but rather to an inflation of the term driven by the data. Indeed, we see from the last column of Figure 14 that while the relative size of f U f_U f U f V f_V f V x x x Figure 11 , this agrees with the interpretation given in § 2.6 , confirming that sub-leading modes can lead to undesirable effects.12
Taken together, Figure 13 and Figure 14 show that the initial condition, and therefore the dependence on the initialisation, does not yield a significant bias with the initialisation shown in Figure 3 and Figure 4 , but it does contribute to the total uncertainty with a relative weight that depends on the flavour and on the training time. We stress that the validity of this conclusion is not confined to the analytical solution, but extends to the trained network to the extent that the former is a good approximation of the training dynamics. Furthermore, small corrections due to the discrete updates and to the deviations from eq. (17) do not change the qualitative picture of this result. How the latter changes with the initialisation is discussed in the next section.
Dependence on the initialisation
The initialisation of the parameters of the NN defines the initial condition from which every replica starts its training course. Inspired by results found with Deep Linear Networks (DLNs) [63 , 64 , 65 ] and more complicated architectures such as those used in Large Language Models [75 ] , in this section we explore the dependence of the training dynamics upon variations of the initialisation distribution of the NN parameters.
To do so, we introduce to the setup of § 2 a power-law scaling of the variance of the Glorot normal distribution of the NN weights, defined as
σ ℓ 2 = ( 2 n ℓ + n ℓ − 1 ) γ (68) \sigma_\ell^2 = \left(\frac{2}{n_\ell+n_{\ell-1}}\right)^\gamma \tag*{(68)} σ ℓ 2 = ( n ℓ + n ℓ − 1 2 ) γ (68) where γ \gamma γ γ = 1 \gamma=1 γ = 1 Table 1 , eq. (68) allows us to explore two different regimes of initialisation widths.13 γ ≲ 1 \gamma\lesssim1 γ ≲ 1 γ ≳ 1 \gamma\gtrsim1 γ ≳ 1 γ → ∞ \gamma\to\infty γ → ∞
In § 5.1 we look at the PDFs at initialisation with different values of γ \gamma γ § 5.2 we inspect the impact on the training dynamics. We validate the analytical solution at different values of γ \gamma γ § 5.3 , and in § 5.6 we explore which figures of merit can be used to identify an optimal value of γ \gamma γ § 5.7 where we clarify the connection between the limit γ → ∞ \gamma\to\infty γ → ∞
A linear model at initialisation
We begin our discussion by studying the dependence on γ \gamma γ Figure 15 we show the one-sigma widths of the network at initialisation as functions of x x x γ \gamma γ γ = 0.001 \gamma=0.001 γ = 0.001 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 Figure 16 the ratio of the one-sigma curves for consecutive values of γ \gamma γ Figure 3 , and the dependence on γ \gamma γ
Figure 15. One-sigma width σ [ x f ( x ) ] \sigma[x f(x)] σ [ x f ( x )] x x x Q 0 Q_0 Q 0 N r e p = 120 N_{\mathrm{rep}}=120 N rep = 120 γ \gamma γ (68) . Both axes are logarithmic.
Figure 16. Ratio of the widths of Figure 15 for consecutive values γ < γ ′ \gamma< \gamma^\prime γ < γ ′ 1 / r = 65.21 1/r = 65.21 1/ r = 65.21 (69) . The ratio is close to this value for γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 γ \gamma γ x x x log x \log x log x
First, the width of the ensemble at initialisation decreases systematically with increasing values of γ \gamma γ x x x x = 1 x = 1 x = 1 x x x γ \gamma γ x Σ x\Sigma x Σ σ ∼ 10 0 ÷ 10 1 \sigma\sim10^0 \div10^1 σ ∼ 1 0 0 ÷ 1 0 1 γ = 0.001 \gamma= 0.001 γ = 0.001 γ = 1 \gamma= 1 γ = 1 σ ∼ 10 − 7 ÷ 10 − 6 \sigma\sim10^{-7} \div10^{-6} σ ∼ 1 0 − 7 ÷ 1 0 − 6 γ = 5 \gamma= 5 γ = 5
Second, the ratio of the widths for consecutive values of γ \gamma γ γ \gamma γ x x x γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 σ γ / σ γ + 1 ≈ 65.21 \sigma_\gamma/\sigma_{\gamma+1} \approx65.21 σ γ / σ γ + 1 ≈ 65.21 Figure 16 with a dotted horizontal line – whereas smaller values of γ \gamma γ x x x γ \gamma γ x x x Figure 15 appears as a rigid vertical translation.
Both facts can be understood by noting that, for large values of γ \gamma γ i.e. a product of weight matrices acting on the input vector X = ( x , log x ) X = (x, \log x) X = ( x , log x ) W ( ℓ ) = σ ℓ Ω ( ℓ ) W^{(\ell)} = \sigma_\ell\Omega^{(\ell)} W ( ℓ ) = σ ℓ Ω ( ℓ ) Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ ) γ \gamma γ γ \gamma γ f γ f_\gamma f γ
f γ ( x ) = r γ Ω ( L ) ⋯ Ω ( 1 ) X , r ≡ [ ∏ ℓ = 1 L 2 n ℓ + n ℓ − 1 ] 1 / 2 , (69) f_\gamma(x) = r^\gamma\Omega^{(L)} \cdots\Omega^{(1)} X,\qquad r \equiv\left[\prod_{\ell=1}^{L} \frac{2}{n_\ell+n_{\ell-1}}\right]^{1/2}, \tag*{(69)} f γ ( x ) = r γ Ω ( L ) ⋯ Ω ( 1 ) X , r ≡ [ ℓ = 1 ∏ L n ℓ + n ℓ − 1 2 ] 1/2 , (69) where r ≤ 1 r \le1 r ≤ 1 f γ = r γ − γ ′ f γ ′ f_\gamma= r^{\gamma-\gamma^\prime}f_{\gamma^\prime} f γ = r γ − γ ′ f γ ′ γ \gamma γ r γ r^\gamma r γ γ \gamma γ 1 / r = 65.21 1/r = 65.21 1/ r = 65.21 Figure 16 . For smaller values of γ \gamma γ x x x ∣ log x ∣ \lvert\log x\rvert ∣ log x ∣ (69) and the saturation of the internal layers are discussed in detail in Section A.
In summary, as long as the network operates in its linear regime, γ \gamma γ x x x
Training dynamics at different initialisations
We now turn to the implications of the initialisation on the training dynamics. We train six fits with plain GD and the same settings as the baseline fit, Table 1 , which differ only in the value of γ \gamma γ Figure 17 as functions of the training epoch for a single representative replica of each of the six values of γ \gamma γ
Figure 17. Training (solid) and validation (dashed) losses as functions of the epoch, for a single representative replica. Each curve represents one of the six values of γ \gamma γ γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 γ \gamma γ
Two regimes can be clearly distinguished from this comparison. For γ ≲ 2 \gamma\lesssim2 γ ≲ 2 Figure 5 . By contrast, the cases with γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 γ \gamma γ τ e s c \tau_{\mathrm{esc}} τ esc L = 3 \mathcal{L} = 3 L = 3 τ e s c ≈ 2500 \tau_{\mathrm{esc}} \approx2500 τ esc ≈ 2500 5400 5400 5400 ≳ 10 4 \gtrsim10^{4} ≳ 1 0 4 γ = 3 \gamma= 3 γ = 3 4 4 4 5 5 5 14 γ = 2 \gamma= 2 γ = 2 γ \gamma γ γ \gamma γ
The two regimes are mirrored by the spectrum of the NTK, shown in Figure 18 for the leading 50 eigenvalues. For γ ≲ 1 \gamma\lesssim1 γ ≲ 1 γ = 0.001 \gamma= 0.001 γ = 0.001 γ = 1 \gamma= 1 γ = 1 γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 τ e s c \tau_{\mathrm{esc}} τ esc γ = 2 \gamma= 2 γ = 2 10 3 10^{3} 1 0 3 Figure 31 , we show that during this first plateau the NN is operating in its linear regime, and therefore it is unable to reconstruct the underlying law. The sharp transient observed in both the loss function (Figure 17 ) and eigenvalues (Figure 18 ) signals the departure from this linear regime.
Figure 18. Training trajectories of the leading 50 eigenvalues λ ( r ) \lambda^{(r)} λ ( r ) γ \gamma γ (68) . The colour encodes the eigenvalue rank r r r Figure 6 . The curves are medians over the N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120 τ e s c \tau_{\mathrm{esc}} τ esc τ = 0 \tau= 0 τ = 0 10 − 10 10^{-10} 1 0 − 10 Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) Figure 8 for any of these fits.
Figure 31. The singlet (upper block) and gluon (lower block) components of six modes q ( r ) q^{(r)} q ( r ) (27) , as functions of x x x γ \gamma γ τ e s c \tau_{\mathrm{esc}} τ esc τ e n d = 2 × 10 4 \tau_{\mathrm{end}} = 2 \times10^{4} τ end = 2 × 1 0 4 γ \gamma γ N o u t N_{\mathrm{out}} N out ∣ q ( r ) ∣ |q^{(r)}| ∣ q ( r ) ∣ + 1 +1 + 1 log x \log x log x
The displacement of the parameters from their initial values, shown in Figure 19 , completes the picture. By the end of the training the parameters have moved by an amount of the same order for all values of γ \gamma γ ∥ θ 0 ∥ \lVert\theta_0\rVert ∥ θ 0 ∥ γ ≲ 1 \gamma\lesssim1 γ ≲ 1 ∥ θ 0 ∥ \lVert\theta_0\rVert ∥ θ 0 ∥ ∥ θ τ − θ 0 ∥ / ∥ θ 0 ∥ ≈ 0.2 \lVert\theta_{\tau} - \theta_0\rVert/\lVert\theta_0\rVert\approx0.2 ∥ θ τ − θ 0 ∥ / ∥ θ 0 ∥ ≈ 0.2 0.6 0.6 0.6 γ = 0.001 \gamma= 0.001 γ = 0.001 γ = 1 \gamma= 1 γ = 1
Figure 19. Displacements of the parameters in Euclidean space as a function of training time and for the different initialisation widths. (Left) The distance ∥ θ τ − θ 0 ∥ \lVert\theta_{\tau} - \theta_{0}\rVert ∥ θ τ − θ 0 ∥ ∥ θ τ − θ 0 ∥ / ∥ θ 0 ∥ \lVert\theta_{\tau} - \theta_{0}\rVert/\lVert\theta_{0}\rVert ∥ θ τ − θ 0 ∥ / ∥ θ 0 ∥ N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120 τ e s c \tau_{\mathrm{esc}} τ esc τ = 0 \tau= 0 τ = 0
training is well described by the NTK equations. This also explains why, in Figure 12 , the analytical solution with the NTK frozen at initialisation provides a good description of the training. For γ ≥ 2 \gamma\ge2 γ ≥ 2 ∥ θ 0 ∥ \lVert\theta_0\rVert ∥ θ 0 ∥ γ = 5 \gamma= 5 γ = 5 § 5.3 .
We remark that, for γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 [63 ] for deep linear networks with a small initialisation. We however clarify that, despite the qualitative similarity, the two phenomena differ in nature. In Ref. [63 ] , the network is linear by construction: the successive plateaus separate saddles corresponding to linear maps of increasing rank, and the transitions between them remain linear maps throughout. Here, instead, we observe a single plateau, along which the network operates in the linear regime discussed in Section A, and the escape coincides with the saturation of the hidden layers (see the right panel of Figure 30 ). What we observe is therefore the departure from the deep linear description, rather than a transition between two of its saddles.
Figure 30. (Left) The dimension of the image of the NTK, dim Im ( Θ ) \dim\operatorname{Im}(\Theta) dim Im ( Θ ) γ \gamma γ Figure 8 . Preprocessing factors are excluded from the NTK. The curve with γ = 1 \gamma= 1 γ = 1 Figure 8 without preprocessing factors. (Right) The saturation s α ( ℓ ) s_{\alpha}^{(\ell)} s α ( ℓ ) (82) , averaged over the input grid in x x x τ = 0 \tau= 0 τ = 0 (83) . In both panels, solid lines are medians over the N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120 τ e s c \tau_{\mathrm{esc}} τ esc τ = 0 \tau= 0 τ = 0
We therefore conclude that changing γ \gamma γ
Validity of the frozen-NTK description
An immediate question is whether and to what extent the analytical solution discussed in § 4.5 can be equally applied to the fits with different initialisations. We answer this question by investigating two complementary aspects: the first is whether the analytical solution can reproduce the numerical training , and the second is whether it can reproduce the underlying law f ( i n ) f^{(\mathrm{in})} f ( in ) γ \gamma γ τ N T K = 0 \tau_{\mathrm{NTK}} = 0 τ NTK = 0 τ N T K = τ e n d \tau_{\mathrm{NTK}} = \tau_{\mathrm{end}} τ NTK = τ end f 0 f_0 f 0 τ N T K \tau_{\mathrm{NTK}} τ NTK
Freezing at initialisation
We begin with the NTK frozen at initialisation. In Figure 20 we compare the loss function of the trained network with that obtained with the analytical solution along the whole training (τ e n d = 20000 \tau_{\mathrm{end}} = 20000 τ end = 20000 γ \gamma γ γ = 0.001 \gamma= 0.001 γ = 0.001 γ = 1 \gamma= 1 γ = 1 Figure 12 and supported by Figure 19 – that the fit is lazy from the beginning with γ ∼ 1 \gamma\sim1 γ ∼ 1 γ = 2 \gamma= 2 γ = 2
Figure 20. Training loss of the analytical solution of eq. (30) with the NTK frozen at initialisation (dashed) against the loss of the empirically trained network (solid), one panel per initialisation exponent γ \gamma γ N r e p = 120 N_{\rm rep}=120 N rep = 120 τ e s c \tau_{\rm esc} τ esc Figure 17 , up to the replacement of the single representative replica by the ensemble median.
The origin of this behaviour can be understood by inspecting the shapes of the resulting PDFs. In Figure 21 we compare the analytical and trained solutions of the singlet for γ = 1 \gamma= 1 γ = 1 3 3 3 5 5 5 γ = 1 \gamma= 1 γ = 1 Figure 12 for the baseline fit and explaining the agreement of the two loss functions in Figure 20 . On the other hand, for γ = 5 \gamma= 5 γ = 5 τ = 500 \tau= 500 τ = 500 τ = 5000 \tau= 5000 τ = 5000 x x x τ = 20000 \tau= 20000 τ = 20000 γ = 3 \gamma= 3 γ = 3 τ e s c ≃ 2500 \tau_{\mathrm{esc}} \simeq2500 τ esc ≃ 2500 τ = 500 \tau= 500 τ = 500 τ = 5000 \tau= 5000 τ = 5000
Figure 21. Analytical (blue) versus empirically trained (orange) solution of the singlet Σ \Sigma Σ τ N T K = 0 \tau_{\mathrm{NTK}} = 0 τ NTK = 0 γ = 1 \gamma= 1 γ = 1 γ = 3 \gamma= 3 γ = 3 γ = 5 \gamma= 5 γ = 5 τ \tau τ τ \tau τ N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120 f ( i n ) f^{(\mathrm{in})} f ( in )
The two examples given above show that, even for large values of γ \gamma γ Figure 18 and Figure 31 ), and the analytical solution provides an accurate description of the training dynamics. However, the dimension of the image of the NTK of the full model at initialisation, computed with the tolerance of Figure 8 ,
γ 0.001 1 2 3 4 5 dim Im ( Θ ) 134 77 51 33 25 15 , (70) \begin{aligned}
\gamma&\quad0.001\quad1\quad2\quad3\quad4\quad5\\
\dim\operatorname{Im}(\Theta) &\quad134\quad77\quad51\quad33\quad25\quad15, \tag*{(70)}
\end{aligned} γ dim Im ( Θ ) 0.001 1 2 3 4 5 134 77 51 33 25 15 , (70) decreases rapidly with γ \gamma γ γ \gamma γ f ( i n ) f^{(\mathrm{in})} f ( in )
V l i n ≡ span { 1 , x , log x } , (71) \mathcal{V}_{\rm lin}\equiv\operatorname{span}\{1,x,\log x\}, \tag*{(71)} V lin ≡ span { 1 , x , log x } , (71) which bounds the image of the network alone to 3 N o u t = 24 3N_{\rm out}=24 3 N out = 24 15 V l i n \mathcal{V}_{\rm lin} V lin
Freezing at the end of training
At the opposite end of the training trajectory the situation is reversed. Also in this case we show, in Figure 22 , the analytical loss function compared with the trained one. We clearly see that for γ = 0.001 \gamma=0.001 γ = 0.001 γ = 1 \gamma=1 γ = 1 γ \gamma γ
Figure 22. Same as Figure 20 , but with the NTK frozen at the last stored epoch, τ N T K = τ e n d = 2 × 10 4 \tau_{\mathrm{NTK}}=\tau_{\mathrm{end}}=2\times10^{4} τ NTK = τ end = 2 × 1 0 4
The resulting PDFs are shown in Figure 23 with the same layout as in Figure 21 , now with the NTK frozen at the last stored epoch. We observe that in the last column, where the dominant modes are indeed learned by the analytical solution, the two solutions agree to a good approximation for all three values of γ \gamma γ γ \gamma γ γ = 1 \gamma= 1 γ = 1 Figure 33 and k ≥ 16 k \ge16 k ≥ 16 τ e s c \tau_{\mathrm{esc}} τ esc γ \gamma γ (17) becomes an accurate approximation again. The analytical solution therefore converges much faster than the trained network, and the behaviour of Ref. [44 ] is recovered.
Figure 23. Same as Figure 21 , but with the NTK frozen at the last stored epoch, τ N T K = τ e n d = 2 × 10 4 \tau_{\mathrm{NTK}}=\tau_{\mathrm{end}}=2\times10^{4} τ NTK = τ end = 2 × 1 0 4
Figure 33. Subspace drift d k d_k d k (150) , measuring how far the k k k Θ τ \Theta_\tau Θ τ k k k k k k 24 = 3 × 8 24 = 3 \times8 24 = 3 × 8 k = 1 k = 1 k = 1 k = 2 k = 2 k = 2 k k k dim Im ( Θ 0 ) \dim\operatorname{Im}(\Theta_0) dim Im ( Θ 0 ) N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120 τ e s c \tau_{\mathrm{esc}} τ esc Figure 18 , with the neural network alone. The initialisation is omitted, d k d_k d k
Finally, the first two columns of Figure 23 reveal a complementary aspect of this story. Evolved for only 500 epochs, the analytical solution built on the late NTK has already reached the underlying law at γ = 3 \gamma= 3 γ = 3 γ = 5 \gamma= 5 γ = 5 V l i n \mathcal{V}_{\mathrm{lin}} V lin
Is there an optimal initialisation?
We provided evidence that the choice of γ \gamma γ Figure 17 ) and the expressiveness of the network in its early learning steps (Figure 21 and Section B.2). We also showed that these two aspects are related, so that changing one inevitably affects the other. We therefore ask whether an optimal value of γ \gamma γ
the lowest validation loss attained, the efficiency of training, and the analytical tractability of the training dynamics. We illustrate these properties through Figure 24 and Figure 25 , which we discuss in turn.
Figure 24. (Left) Minimum over the first 2 × 10 4 2 \times10^{4} 2 × 1 0 4 γ = 1 \gamma= 1 γ = 1 ± 0.05 \pm0.05 ± 0.05 γ = 1 \gamma= 1 γ = 1 γ = 5 \gamma= 5 γ = 5 2 × 10 4 2 \times10^{4} 2 × 1 0 4 η h max \eta h_{\max} η h m a x h max h_{\max} h m a x H ⊥ H^{\perp} H ⊥ (73) . We use black circles for the fixed learning rate, η = 10 − 6 \eta= 10^{-6} η = 1 0 − 6 η h max \eta h_{\max} η h m a x γ = 1 \gamma= 1 γ = 1 η = 10 − 7 \eta= 10^{-7} η = 1 0 − 7 γ = 0.001 , 0.1 \gamma= 0.001, 0.1 γ = 0.001 , 0.1 2 × 10 − 7 2 \times10^{-7} 2 × 1 0 − 7 γ = 0.5 \gamma= 0.5 γ = 0.5
Figure 25. Comparison between the trained solution with the contribution of the data removed, f − V ( T ) Y f - V(T)Y f − V ( T ) Y U ( T ) f 0 U(T)f_0 U ( T ) f 0 γ = 0.5 \gamma= 0.5 γ = 0.5 0.75 0.75 0.75 1 1 1 1.25 1.25 1.25 1.5 1.5 1.5 τ N T K = τ e n d = 2 × 10 4 \tau_{\mathrm{NTK}} = \tau_{\mathrm{end}} = 2 \times10^4 τ NTK = τ end = 2 × 1 0 4 T T T τ = τ e n d \tau= \tau_{\mathrm{end}} τ = τ end N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120
We start from the two left panels of Figure 24 , which show, for different values of γ \gamma γ γ = 1 \gamma= 1 γ = 1 L t r = 2.3 L_{\mathrm{tr}} = 2.3 L tr = 2.3 16 γ min ∼ 0.6 \gamma_{\min} \sim0.6 γ m i n ∼ 0.6 γ \gamma γ γ = 1 \gamma= 1 γ = 1 γ \gamma γ γ = 1 \gamma= 1 γ = 1 γ \gamma γ
Comparing different initialisations at a fixed learning rate requires, however, some care. In fact, the analytical solution shows that the size of a GD step along a given mode q ( i ) q^{(i)} q ( i ) (66) , with t = η τ t = \eta\tau t = η τ e τ ( i ) e_{\tau}^{(i)} e τ ( i ) i i i
e τ + 1 ( i ) = ( 1 − 2 η h ( i ) ) e τ ( i ) , (72) e_{\tau+1}^{(i)} = \left(1 - 2\eta h^{(i)}\right)e_{\tau}^{(i)}, \tag*{(72)} e τ + 1 ( i ) = ( 1 − 2 η h ( i ) ) e τ ( i ) , (72) so that every mode converges only if ∣ 1 − 2 η h ( i ) ∣ < 1 \lvert1 - 2\eta h^{(i)}\rvert< 1 ∣ 1 − 2 η h ( i ) ∣ < 1
η h max < 1 , (73) \eta h_{\max} < 1, \tag*{(73)} η h m a x < 1 , (73) where h max h_{\max} h m a x H ⊥ H^{\perp} H ⊥ 1 / 2 < η h < 1 1/2 < \eta h < 1 1/2 < η h < 1 η h > 1 \eta h > 1 η h > 1
Being derived for a frozen NTK, eq. (73) does not imply that a network starting above the bound diverges: the large steps also modify the NTK, and can lower h max h_{\max} h m a x η h max ≃ 1 \eta h_{\max} \simeq1 η h m a x ≃ 1 edge of stability [76 , 77 , 78 ] . In the right-most panel of Figure 24 we show η h max \eta h_{\max} η h m a x γ \gamma γ γ = 1 \gamma= 1 γ = 1 η h max = 0.95 \eta h_{\max} = 0.95 η h m a x = 0.95 γ < 1 \gamma< 1 γ < 1 η h max \eta h_{\max} η h m a x
to one by the end of the training. For γ > 1 \gamma> 1 γ > 1 γ \gamma γ γ ≥ 2 \gamma\ge2 γ ≥ 2 2 × 10 4 2 \times10^{4} 2 × 1 0 4 γ \gamma γ [76 , 77 ] , and that this condition is reached within the first few steps in those fits for which the NTK is almost constant.17 η h max \eta h_{\max} η h m a x
For a frozen NTK, the stability condition of eq. (73) thus guarantees that gradient descent converges to the same solution as the gradient flow of eq. (30) , the two differing only by a short transient in the directions of the larger modes. Any discrepancy between the analytical and the trained solutions must therefore be ascribed to the evolution of the NTK: if the latter is constant – or slowly varying – during training, then the analytical solution is guaranteed to follow the trained network. To check for which values of γ \gamma γ f − V ( t ) Y f - V(t)Y f − V ( t ) Y U ( t ) f 0 U(t)f_0 U ( t ) f 0 Figure 25 . We observe that the residual obtained by removing the contribution of the data from the trained network is closest to the analytical prediction for γ = 0.5 \gamma= 0.5 γ = 0.5 γ = 0.75 \gamma= 0.75 γ = 0.75 γ \gamma γ γ \gamma γ Figure 24 . The initialisations that yield the best and the fastest fits are therefore also those for which we gain the most analytical control over the training, and for which the contribution of the initial condition to the final solution can be predicted. The only exception is T 3 T_3 T 3
We finally note that the observations made here hold insofar as vanilla GD is adopted. Whether the same behaviour is expected once the optimisation algorithm is replaced by a more efficient one (e.g. Adam or Nadam as in the NNPDF case) is unclear from the present analysis. This aspect will also be investigated in future work.
A wide NN is not a substitute for a small initialisation
We devote the last part of this section to laying out the connection between the parameter γ \gamma γ (68) and the width of the NN, that is, the number of nodes in its hidden layers. In fact, one may be led to believe that increasing the number of internal nodes in eq. (68) is equivalent to increasing the value of γ \gamma γ
In Figure 26 we compare the spectrum of the NTK at initialisation for γ = 1 \gamma= 1 γ = 1 2 2 2 3 3 3 γ = 1 \gamma= 1 γ = 1 γ \gamma γ (68) , the variance of the parameter distribution of the wide network is σ w i d e 2 = 2 × 10 − 3 \sigma_{\mathrm{wide}}^{2} = 2 \times10^{-3} σ wide 2 = 2 × 1 0 − 3 γ ∼ 2 \gamma\sim2 γ ∼ 2 γ \gamma γ γ = 2 \gamma= 2 γ = 2 γ \gamma γ
Figure 26. Eigenvalues of the NTK of the neural network (preprocessing factors excluded) at initialisation, against their rank. Coloured curves are the baseline architecture [25 , 20 , 8 ] at γ = 1 , 2 , 3 \gamma= 1, 2, 3 γ = 1 , 2 , 3 γ = 1 \gamma= 1 γ = 1 N r e p = 120 N_{\rm rep} = 120 N rep = 120 Figure 8 .
The size of the spectrum of the wide network can be understood by considering the kernel of the inner layers, defined as [44 ]
K i α , i ′ α ′ ( ℓ ) = E [ ϕ i α ( ℓ ) ϕ i ′ α ′ ( ℓ ) ] , with ϕ i α ( ℓ ) = ∑ j = 1 n ℓ − 1 W i j ( ℓ − 1 ) ρ j α ( ℓ − 1 ) + b i ( ℓ ) (74) K^{(\ell)}_{i\alpha,i'\alpha'} = \mathbb{E}\left[\phi^{(\ell)}_{i\alpha}\phi^{(\ell)}_{i'\alpha'}\right], \qquad\text{with}\qquad\phi^{(\ell)}_{i\alpha} = \sum_{j=1}^{n_{\ell-1}} W^{(\ell-1)}_{ij}\rho^{(\ell-1)}_{j\alpha} + b^{(\ell)}_i \tag*{(74)} K i α , i ′ α ′ ( ℓ ) = E [ ϕ i α ( ℓ ) ϕ i ′ α ′ ( ℓ ) ] , with ϕ i α ( ℓ ) = j = 1 ∑ n ℓ − 1 W ij ( ℓ − 1 ) ρ j α ( ℓ − 1 ) + b i ( ℓ ) (74) where ρ i α ( ℓ ) ≡ ρ ( ϕ i α ( ℓ ) ) \rho^{(\ell)}_{i\alpha} \equiv\rho\left(\phi^{(\ell)}_{i\alpha}\right) ρ i α ( ℓ ) ≡ ρ ( ϕ i α ( ℓ ) ) ρ \rho ρ ϕ i α ( ℓ ) \phi^{(\ell)}_{i\alpha} ϕ i α ( ℓ ) ℓ \ell ℓ
K i α , i ′ α ′ ( ℓ ) = ∑ j = 1 n ℓ − 1 ∑ j ′ = 1 n ℓ − 1 E [ W i j ( ℓ ) W i ′ j ′ ( ℓ ) ] E [ ρ j α ( ℓ − 1 ) ρ j ′ α ′ ( ℓ − 1 ) ] = δ i i ′ σ ℓ 2 ∑ j = 1 n ℓ − 1 E [ ρ j α ( ℓ − 1 ) ρ j α ′ ( ℓ − 1 ) ] = δ i i ′ σ ℓ 2 n ℓ − 1 E [ ρ 1 α ( ℓ − 1 ) ρ 1 α ′ ( ℓ − 1 ) ] . (75) \begin{aligned}
K^{(\ell)}_{i\alpha,i'\alpha'} &= \sum_{j=1}^{n_{\ell-1}}\sum_{j'=1}^{n_{\ell-1}} \mathbb{E}\left[W^{(\ell)}_{ij}W^{(\ell)}_{i'j'}\right]\mathbb{E}\left[\rho^{(\ell-1)}_{j\alpha}\rho^{(\ell-1)}_{j'\alpha'}\right] \\
&= \delta_{ii'}\sigma_\ell^2\sum_{j=1}^{n_{\ell-1}}\mathbb{E}\left[\rho^{(\ell-1)}_{j\alpha}\rho^{(\ell-1)}_{j\alpha'}\right] \\
&= \delta_{ii'}\sigma_\ell^2 n_{\ell-1}\mathbb{E}\left[\rho^{(\ell-1)}_{1\alpha}\rho^{(\ell-1)}_{1\alpha'}\right]. \tag*{(75)}
\end{aligned} K i α , i ′ α ′ ( ℓ ) = j = 1 ∑ n ℓ − 1 j ′ = 1 ∑ n ℓ − 1 E [ W ij ( ℓ ) W i ′ j ′ ( ℓ ) ] E [ ρ j α ( ℓ − 1 ) ρ j ′ α ′ ( ℓ − 1 ) ] = δ i i ′ σ ℓ 2 j = 1 ∑ n ℓ − 1 E [ ρ j α ( ℓ − 1 ) ρ j α ′ ( ℓ − 1 ) ] = δ i i ′ σ ℓ 2 n ℓ − 1 E [ ρ 1 α ( ℓ − 1 ) ρ 1 α ′ ( ℓ − 1 ) ] . (75) The first equality holds because the weights are i.i.d. within and across layers and the biases are initialised to zero. The second one follows from the definition of σ ℓ \sigma_\ell σ ℓ ℓ − 1 \ell- 1 ℓ − 1 i.i.d. across neurons, all neurons in that layer are statistically equivalent, so the expectation value E [ ρ j α ( ℓ − 1 ) ρ j α ′ ( ℓ − 1 ) ] \mathbb{E}\left[\rho^{(\ell-1)}_{j\alpha}\rho^{(\ell-1)}_{j\alpha'}\right] E [ ρ j α ( ℓ − 1 ) ρ j α ′ ( ℓ − 1 ) ] j j j 18 j j j n ℓ − 1 n_{\ell-1} n ℓ − 1
We also mention that the NTK can be computed by iterating expressions similar to eq. (75) layer by layer, a recursion also known as the forward equation of the NTK (see e.g. Ref. [79 ] ). Being recursive, this construction is such that the scaling of each layer accumulates through the layers. We therefore conclude that the reduction of the variance induced by γ \gamma γ (75) . Specifically, we can show that
K γ ( ℓ ) ∼ O ( n 1 − γ ) , and K w i d e ( k ) ∼ O ( 1 ) , (76) K^{(\ell)}_{\gamma} \sim O\left(n^{1-\gamma}\right), \qquad\text{and}\qquad K^{(k)}_{\rm wide} \sim O(1), \tag*{(76)} K γ ( ℓ ) ∼ O ( n 1 − γ ) , and K wide ( k ) ∼ O ( 1 ) , (76) with n n n O ( 1 ) O(1) O ( 1 ) γ > 1 \gamma> 1 γ > 1 γ \gamma γ
The different scaling properties affect the training dynamics as well. In Figure 27 we show the training losses of the baseline fit and of the wide network just discussed (n 1 = n 2 = 500 n_1 = n_2 = 500 n 1 = n 2 = 500 γ = 1 \gamma= 1 γ = 1 t = η τ t = \eta\tau t = η τ η = 10 − 7 \eta= 10^{-7} η = 1 0 − 7 Figure 26 , in accordance with the stability condition of eq. (73) . We clearly see that the wide network has not reached convergence within the 2 × 10 4 2 \times10^{4} 2 × 1 0 4 γ \gamma γ Figure 17 , confirming that the two prescriptions lead to different training dynamics. We also note that the fair like-for-like comparison is actually given by the right panel, where the difference in the learning rates is accounted for.
Figure 27. Training loss at γ = 1 \gamma= 1 γ = 1 [25 , 20 , 8 ] (learning rate η = 10 − 6 \eta= 10^{-6} η = 1 0 − 6 η = 10 − 7 \eta= 10^{-7} η = 1 0 − 7 t = η τ t = \eta\tau t = η τ
Conclusions and outlook
In this paper we have investigated the training dynamics of a NN through a quantitative analysis of the NTK in a realistic PDF determination, extending a previous study based on a single-flavour toy model. We have generalised the analytical solution of the gradient flow of Ref. [44 ] to the multi-flavour case and validated it against the numerically trained network. Within the NNPDF4.0 methodology [45 , 46 ] , we have exposed the implicit regularisation induced by both gradient descent and parametrisation, and characterised the behaviour of the NN during the training process in terms of the NTK and a second matrix M M M M M M [43 ] , and outlined a geometric interpretation of training. Here we summarise our findings.
First, the analysis of Ref. [44 ] can be extended to a more realistic case, proving the NTK a valuable diagnostic tool to characterise the features learned by the NN during training. Consistent with the previous study, the spectrum of the NTK exhibits a hierarchical structure (Figure 6 ). Contrary to it, however, the rich phase is confined to the first few epochs, where sub-leading modes populate the image of the NTK incrementally; after this short initial transient, the NTK becomes lazy and the training flow can be integrated analytically (eq. (37) ). A similar behaviour is also reflected in the modes of the training Hamiltonian (eq. (27) ), which are able to capture much of the shape of the underlying law already at initialisation (Figure 9 ).
Second, the extension of the analytical solution to the multi-flavour case reproduces the numerically trained network to a good approximation for all flavours, except for the uncertainty of T 3 T_3 T 3 x x x Figure 11 ). Within the lazy training regime, we have used the analytical solution to quantify the impact of the initial condition and of the data on the final solution (Figure 13 ). We have found that the central value of the PDFs is essentially shaped by the contribution of the data, while the bias induced by the initialisation prescription, as defined in § 2 , is negligible. For the uncertainties, the contribution of the initial condition is suppressed in the region covered by the data – although not for T 3 T_3 T 3 Figure 14 ).
Third, the implicit regularisation of gradient descent amounts to the incremental inclusion of NTK modes, each learned at a timescale set by the size of the corresponding eigenvalue (eq. (37) ). We have demonstrated that only about 30 of the leading modes are needed to reconstruct the underlying law (Figure 11 ), which is in stark contrast with the 763 parameters of the network. This suggests that the total number of parameters largely overestimates the effective degrees of freedom used to reconstruct the underlying law, and provides an empirical evidence for the information-theoretic arguments of Ref. [58 ] . Beyond the leading modes, we have shown that the sub-leading ones are poorly constrained by the data (eq. (38) ) and their inclusion carries noise into the solution (Figure 11 ). This shows that early stopping sets the trade-off between the signal and the noise of the data.
Finally, we have shown that the width of the distribution of the parameters at initialisation has a considerable impact on the dynamics and efficiency of training. A distribution that is too narrow forces the network to operate in a regime linear in its input variables (eq. (69) ). The loss function then stalls on a plateau until the inner layers saturate (Figure 30 ), the loss function drops (Figure 17 ), and the NN departs from the linear regime. This is reflected in the spectrum of the NTK: the eigenvalues closely follow the behaviour of the loss function (Figure 18 ), while the corresponding eigenvectors are confined to the span { 1 , x , log x } \{1, x, \log x\} { 1 , x , log x } Figure 32 ). Conversely, an initialisation slightly wider than the standard one yields the best and fastest fits, which are also those over which we gain the most analytical control. We have noted that the sharp transition in both loss function and eigenvalues was already observed in Ref. [44 ] , and that, despite the qualitative similarity, it is different in nature from the saddle-to-saddle dynamics of Refs. [63 , 64 , 65 ] . We have then clarified that a small initialisation is not equivalent to a large width, as the NTK scales differently in the two cases.
Figure 32. Residual ρ ( r ) \rho^{(r)} ρ ( r ) (149) as a function of the rank r r r τ e s c \tau_{\rm esc} τ esc N r e p = 120 N_{\rm rep} = 120 N rep = 120 γ = 4 \gamma= 4 γ = 4 dim Im ( Θ ) \dim\operatorname{Im}(\Theta) dim Im ( Θ ) Figure 30 ). The dotted vertical lines mark the boundaries of the tiers as in Figure 29 . The wide band on the γ = 5 \gamma= 5 γ = 5 V l i n \mathcal{V}_{\rm lin} V lin
Our analysis can be extended in many respects. The analytical treatment can be improved by solving the gradient flow with more sophisticated optimisers or stochastic variants of gradient descent, and by restoring the sum rules and physical constraints removed here. Concerning possible applications, the operator M M M [34 ] , identifying the directions in PDF space affected by the inconsistencies. Moreover, the NTK can be computed for any parametrisation, such as the fixed functional forms of FPPDF [24 ] , as well as for NNs used in other tasks in particle physics.
We view these results as part of a broader change of perspective aiming at a quantitative interpretation of ML [18 ] , where the relationship with physics is reciprocal. In fact, physics can help ML with formal tools inherited from it [55 , 56 , 72 , 79 , 80 ] , and offers in turn realistic case studies – such as the PDF fit analysed here – against which these ideas can be tested, contributing to an emerging theory of deep learning [81 ] .
Code and data availability
The fits presented in this paper were performed with the open-source NNPDF code [46 , 82 ] , using the ntk_colibri branch (commit 16786e8e7e). The NTK analysis was carried out with ntkpdf [83 ] (commit 90bda63b35), a specialised interface built on top of Colibri [36 ] (ntk branch, commit fb82175704). The fit outputs used in this work are available from the authors upon request.
Acknowledgements
We are very grateful to M Ubiali, M N Costantini, V S Sánchez, and E Cole for helpful discussions during the early stages of this work and for support with the Colibri framework. LDD is supported by an STFC Consolidated Grant (ST/T000600/1, ST/X000494/1). AC used the AI assistant Claude (Anthropic) for manuscript editing and for ensuring consistency throughout the text, and for assistance with coding and the production of the figures. All scientific content, including the methodology, results and their interpretation, is the sole responsibility of the authors.
The deep linear network approximation
Throughout Section 5 we relied on the DLN approximation, in which the network is assumed to operate in the linear regime of its activations. Here we collect the derivations behind the statements made there. We first derive the scaling with γ \gamma γ V l i n V_{\mathrm{lin}} V lin (71) , compute the leading tier of the NTK spectrum exactly, and finally derive a compact form for the NTK of a DLN that reduces its eigenvalue problem to 3 N o u t 3N_{\mathrm{out}} 3 N out
Deep linear networks
In the linear regime, the pre-activations can be written as
ϕ i α ( 1 ) = ∑ j = 1 2 W i j ( 1 ) X j α + b i ( 1 ) , \phi_{i\alpha}^{(1)} = \sum_{j=1}^{2} W_{ij}^{(1)} X_{j\alpha} + b_i^{(1)}, ϕ i α ( 1 ) = j = 1 ∑ 2 W ij ( 1 ) X j α + b i ( 1 ) , ϕ i α ( ℓ ) = ∑ j = 1 n ℓ − 1 W i j ( ℓ ) ϕ j α ( ℓ − 1 ) + b i ( ℓ ) , with ϕ i α ( L ) ≡ N N i ( x α ) , (77) \phi_{i\alpha}^{(\ell)} = \sum_{j=1}^{n_{\ell-1}} W_{ij}^{(\ell)} \phi_{j\alpha}^{(\ell-1)} + b_i^{(\ell)}, \qquad\text{with}\qquad\phi_{i\alpha}^{(L)} \equiv\mathrm{NN}_i(x_\alpha), \tag*{(77)} ϕ i α ( ℓ ) = j = 1 ∑ n ℓ − 1 W ij ( ℓ ) ϕ j α ( ℓ − 1 ) + b i ( ℓ ) , with ϕ i α ( L ) ≡ NN i ( x α ) , (77) where X α = ( x α , log x α ) X_\alpha= (x_\alpha,\log x_\alpha) X α = ( x α , log x α ) W ( ℓ ) ∈ R n ℓ × n ℓ − 1 W^{(\ell)} \in\mathbb{R}^{n_\ell\times n_{\ell-1}} W ( ℓ ) ∈ R n ℓ × n ℓ − 1 b ( ℓ ) ∈ R n ℓ b^{(\ell)} \in\mathbb{R}^{n_\ell} b ( ℓ ) ∈ R n ℓ ℓ \ell ℓ (77) becomes matrix-vector multiplications and additions. With this convention, one can verify that the output of the network can be written as
ϕ ( L ) = W ( L ) ϕ ( L − 1 ) + b ( L ) = W ( L ) ( W ( L − 1 ) ϕ ( L − 2 ) + b ( L − 1 ) ) + b ( L ) = W ( L ) ( W ( L − 1 ) ( W ( L − 2 ) ( … ) + b ( L − 2 ) ) + b ( L − 1 ) ) + b ( L ) = ⋮ = b ( L ) + ∑ ℓ = 1 L − 1 ( ∏ m = L ℓ + 1 W ( m ) ) b ( ℓ ) + ( ∏ ℓ = 0 L − 1 W ( L − ℓ ) ) X , (78) \begin{aligned}
\phi^{(L)} &= W^{(L)}\phi^{(L-1)} + b^{(L)} \\
&= W^{(L)}\left(W^{(L-1)}\phi^{(L-2)} + b^{(L-1)}\right) + b^{(L)} \\
&= W^{(L)}\left(W^{(L-1)}\left(W^{(L-2)}(\ldots) + b^{(L-2)}\right) + b^{(L-1)}\right) + b^{(L)} \\
&\phantom{=} \vdots\\
&= b^{(L)} + \sum_{\ell=1}^{L-1}\left(\prod_{m=L}^{\ell+1} W^{(m)}\right)b^{(\ell)} + \left(\prod_{\ell=0}^{L-1} W^{(L-\ell)}\right)X, \tag*{(78)}
\end{aligned} ϕ ( L ) = W ( L ) ϕ ( L − 1 ) + b ( L ) = W ( L ) ( W ( L − 1 ) ϕ ( L − 2 ) + b ( L − 1 ) ) + b ( L ) = W ( L ) ( W ( L − 1 ) ( W ( L − 2 ) ( … ) + b ( L − 2 ) ) + b ( L − 1 ) ) + b ( L ) = ⋮ = b ( L ) + ℓ = 1 ∑ L − 1 ( m = L ∏ ℓ + 1 W ( m ) ) b ( ℓ ) + ( ℓ = 0 ∏ L − 1 W ( L − ℓ ) ) X , (78) which, once the biases vanish as they do at initialisation, reduces to
N N ( x ) = W ( L ) ⋯ W ( 1 ) X = A x + B log x , (79) \mathrm{NN}(x) = W^{(L)}\cdots W^{(1)}X = Ax + B\log x, \tag*{(79)} NN ( x ) = W ( L ) ⋯ W ( 1 ) X = A x + B log x , (79) with A , B ∈ R n L A,B \in\mathbb{R}^{n_L} A , B ∈ R n L
The output of the network in eq. (78) decomposes into three terms. The first one is x x x N o u t N_{\mathrm{out}} N out x x x ℓ \ell ℓ ℓ + 1 \ell+1 ℓ + 1 (79) , which carries the x x x (78) , only the first one is not affected by the scaling with γ \gamma γ
The scaling with γ \gamma γ (79) . The weights are independent Gaussian variables with zero mean and variance σ ℓ 2 \sigma_\ell^2 σ ℓ 2 (68) , so that we can write W ( ℓ ) = σ ℓ Ω ( ℓ ) W^{(\ell)} = \sigma_\ell\Omega^{(\ell)} W ( ℓ ) = σ ℓ Ω ( ℓ ) Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ ) γ \gamma γ f γ f_\gamma f γ γ \gamma γ γ \gamma γ
f γ ( x ) = ( ∏ ℓ = 1 L σ ℓ ) Ω ( L ) ⋯ Ω ( 1 ) X = r γ Ω ( L ) ⋯ Ω ( 1 ) X , (80) f_\gamma(x) = \left(\prod_{\ell=1}^{L}\sigma_\ell\right)\Omega^{(L)}\cdots\Omega^{(1)}X = r^\gamma\Omega^{(L)}\cdots\Omega^{(1)}X, \tag*{(80)} f γ ( x ) = ( ℓ = 1 ∏ L σ ℓ ) Ω ( L ) ⋯ Ω ( 1 ) X = r γ Ω ( L ) ⋯ Ω ( 1 ) X , (80) with r r r (69) . The product of the matrices Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ ) γ \gamma γ 19 Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ ) γ \gamma γ f γ = r γ − γ ′ f γ ′ f_{\gamma} = r^{\gamma-\gamma'} f_{\gamma'} f γ = r γ − γ ′ f γ ′ γ \gamma γ r ≤ 1 r \le1 r ≤ 1 r γ r^\gamma r γ γ \gamma γ 1 / r 1/r 1/ r Table 1 this evaluates to
r = 1.533 × 10 − 2 , 1 r = 65.21. (81) r = 1.533 \times10^{-2}, \qquad\frac{1}{r} = 65.21. \tag*{(81)} r = 1.533 × 1 0 − 2 , r 1 = 65.21. (81) This is the asymptotic value marked by the dotted line in Figure 16 , which the ratio between the γ = 4 \gamma= 4 γ = 4 γ = 5 \gamma= 5 γ = 5 Figure 15 are parallel on a logarithmic scale for γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 σ γ = r γ − γ ′ σ γ ′ \sigma_\gamma= r^{\gamma-\gamma'}\sigma_{\gamma'} σ γ = r γ − γ ′ σ γ ′ ( γ − γ ′ ) log 10 r (\gamma-\gamma')\log_{10} r ( γ − γ ′ ) log 10 r x x x
Saturation of the hidden layers
The DLN approximation holds as long as the hyperbolic tangents operate in their linear regime. Decreasing γ \gamma γ (69) is no longer valid.
The degree of saturation can be quantified by measuring in which regime the activation functions operate at different γ \gamma γ s ( ℓ ) ( x α ) ≡ s α ( ℓ ) s^{(\ell)}(x_\alpha) \equiv s^{(\ell)}_\alpha s ( ℓ ) ( x α ) ≡ s α ( ℓ ) ℓ \ell ℓ x α x_\alpha x α
s α ( ℓ ) ≡ ⟨ 1 − tanh ′ ϕ α ( ℓ ) ⟩ = ⟨ tanh 2 ϕ α ( ℓ ) ⟩ ≡ 1 n ℓ ∑ i = 1 n ℓ tanh 2 ϕ i α ( ℓ ) , (82) s^{(\ell)}_\alpha\equiv\left\langle1-\tanh' \phi^{(\ell)}_\alpha\right\rangle= \left\langle\tanh^2 \phi^{(\ell)}_\alpha\right\rangle\equiv\frac{1}{n_\ell}\sum_{i=1}^{n_\ell}\tanh^2 \phi^{(\ell)}_{i\alpha}, \tag*{(82)} s α ( ℓ ) ≡ ⟨ 1 − tanh ′ ϕ α ( ℓ ) ⟩ = ⟨ tanh 2 ϕ α ( ℓ ) ⟩ ≡ n ℓ 1 i = 1 ∑ n ℓ tanh 2 ϕ i α ( ℓ ) , (82) that is, one minus the average of the derivative of the hyperbolic tangent over the units of the layer. Here, ϕ i ( ℓ ) ( x α ) ≡ ϕ i α ( ℓ ) \phi^{(\ell)}_i(x_\alpha) \equiv\phi^{(\ell)}_{i\alpha} ϕ i ( ℓ ) ( x α ) ≡ ϕ i α ( ℓ ) i i i ℓ \ell ℓ x α x_\alpha x α Figure 28 we show the saturation of the first (left panel) and second (right panel) hidden layers as functions of x x x γ \gamma γ (82) averaged over x x x
Figure 28. Saturation of the first (left) and second (right) hidden layers as a function of x x x (82) . The six values of γ \gamma γ N r e p = 120 N_{\mathrm{rep}} = 120 N rep = 120 γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 x x x γ \gamma γ x x x log x \log x log x
γ 0.001 1 2 3 4 5 s ( 1 ) 0.732 0.453 0.162 0.022 0.002 0.000 s ( 2 ) 0.806 0.245 0.008 0.000 0.000 0.000 (83) \begin{array}{c|rrrrrr}
\gamma& 0.001 & 1 & 2 & 3 & 4 & 5 \\
s^{(1)} & 0.732 & 0.453 & 0.162 & 0.022 & 0.002 & 0.000 \\
s^{(2)} & 0.806 & 0.245 & 0.008 & 0.000 & 0.000 & 0.000 \tag*{(83)}
\end{array} γ s ( 1 ) s ( 2 ) 0.001 0.732 0.806 1 0.453 0.245 2 0.162 0.008 3 0.022 0.000 4 0.002 0.000 5 0.000 0.000 (83) which provide a more compact summary of the saturation at initialisation. Taken together, Figure 28 and eq. (83) support our previous observations: for γ ≳ 3 \gamma\gtrsim3 γ ≳ 3 x x x γ \gamma γ x x x log x \log x log x γ ≥ 1 \gamma\ge1 γ ≥ 1 γ = 0.001 \gamma= 0.001 γ = 0.001
The saturation of the internal layers is thus closely related to the regime in which the network operates. In fact, we can identify a non-linear regime for γ ≲ 2 \gamma\lesssim2 γ ≲ 2 γ ≳ 3 \gamma\gtrsim3 γ ≳ 3
The span of the tangent plane
It can be easily verified from eq. (78) that the derivatives of the network with respect to the parameters (the tangent vectors of § 3.3 ) can be written as
J i ′ α , μ ( n e t ) ≡ ∂ μ ϕ i ′ α ( L ) = u μ i ′ + v μ i ′ x α + w μ i ′ log x α , μ = 1 , … , N p a r a m s , (84) J^{(\mathrm{net})}_{i'\alpha,\mu} \equiv\partial_{\mu}\phi^{(L)}_{i'\alpha} = u_{\mu i'} + v_{\mu i'}x_{\alpha} + w_{\mu i'}\log x_{\alpha}, \qquad\mu= 1,\ldots,N_{\mathrm{params}}, \tag*{(84)} J i ′ α , μ ( net ) ≡ ∂ μ ϕ i ′ α ( L ) = u μ i ′ + v μ i ′ x α + w μ i ′ log x α , μ = 1 , … , N params , (84) where we have restored the indices for the sake of clarity. Here we use primed indices to label the N o u t = 8 N_{\mathrm{out}} = 8 N out = 8 (5) , whose N f l a v N_{\mathrm{flav}} N flav R ( F K ) ∈ R N o u t × N f l a v R^{(\mathrm{FK})} \in\mathbb{R}^{N_{\mathrm{out}}\times N_{\mathrm{flav}}} R ( FK ) ∈ R N out × N flav x x x
J i α , μ = ∑ i ′ ( R ( F K ) ) i i ′ x α 1 − δ i ′ ( 1 − x α ) ϵ i ′ J i ′ α , μ ( n e t ) , (85) J_{i\alpha,\mu} = \sum_{i'}\left(R^{(\mathrm{FK})}\right)_{ii'}x_{\alpha}^{1-\delta_{i'}}(1-x_{\alpha})^{\epsilon_{i'}}J^{(\mathrm{net})}_{i'\alpha,\mu}, \tag*{(85)} J i α , μ = i ′ ∑ ( R ( FK ) ) i i ′ x α 1 − δ i ′ ( 1 − x α ) ϵ i ′ J i ′ α , μ ( net ) , (85) where the un-primed indices label the N f l a v = 14 N_{\mathrm{flav}} = 14 N flav = 14 δ i ′ = 1 \delta_{i'} = 1 δ i ′ = 1 ϵ i ′ = 0 \epsilon_{i'} = 0 ϵ i ′ = 0 (85) is again a combination of constant and non-constant terms whose coefficients are linear in u u u v v v w w w R ( F K ) R^{(\mathrm{FK})} R ( FK ) N o u t N_{\mathrm{out}} N out
R ( F K ) ( R ( F K ) ) T = diag ( m i ′ ) , i ′ = 1 , … , N o u t . (86) R^{(\mathrm{FK})}\left(R^{(\mathrm{FK})}\right)^{T} = \operatorname{diag}(m_{i'}), \qquad i' = 1,\ldots,N_{\mathrm{out}}. \tag*{(86)} R ( FK ) ( R ( FK ) ) T = diag ( m i ′ ) , i ′ = 1 , … , N out . (86) Here m i ′ m_{i'} m i ′
V = V 15 = V 24 = V 35 , Σ = T 24 = T 35 , photon = 0 , (87) \begin{aligned}
V &= V_{15} = V_{24} = V_{35},\\
\Sigma&= T_{24} = T_{35},\\
\text{photon} &= 0, \tag*{(87)}
\end{aligned} V Σ photon = V 15 = V 24 = V 35 , = T 24 = T 35 , = 0 , (87) and all the other flavours unchanged. We then conclude that m V = 4 m_V = 4 m V = 4 m Σ = 3 m_{\Sigma} = 3 m Σ = 3 (85) when needed.
The tangent vectors in eq. (84) are generated by derivatives with respect to three groups of parameters:
G1) The N o u t N_{\mathrm{out}} N out u μ i u_{\mu i} u μ i
∂ ϕ i ′ α ( L ) ∂ b j ′ ( L ) = δ i ′ j ′ , (88) \frac{\partial\phi^{(L)}_{i'\alpha}}{\partial b^{(L)}_{j'}} = \delta_{i'j'}, \tag*{(88)} ∂ b j ′ ( L ) ∂ ϕ i ′ α ( L ) = δ i ′ j ′ , (88) that is, the identity 1 N o u t \mathbb{1}_{N_{\mathrm{out}}} 1 N out γ \gamma γ G2) The biases of the hidden layers b ( ℓ ) b^{(\ell)} b ( ℓ )
∂ ϕ i ′ α ( L ) ∂ b j ( ℓ ) = [ W ( L ) ⋯ W ( ℓ + 1 ) ] i ′ j = r ℓ + 1 γ [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] i ′ j , ℓ ≠ L , (89) \frac{\partial\phi^{(L)}_{i'\alpha}}{\partial b^{(\ell)}_j}=\left[W^{(L)}\cdots W^{(\ell+1)}\right]_{i'j}^{ }=r_{\ell+1}^{\gamma}\left[\Omega^{(L)}\cdots\Omega^{(\ell+1)}\right]_{i'j}^{ },\qquad\ell\ne L, \tag*{(89)} ∂ b j ( ℓ ) ∂ ϕ i ′ α ( L ) = [ W ( L ) ⋯ W ( ℓ + 1 ) ] i ′ j = r ℓ + 1 γ [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] i ′ j , ℓ = L , (89) where we made explicit the scaling with γ \gamma γ W ( ℓ ) = σ ℓ Ω ( ℓ ) W^{(\ell)}=\sigma_{\ell}\Omega^{(\ell)} W ( ℓ ) = σ ℓ Ω ( ℓ ) r ℓ 2 = ∏ m = ℓ L 2 / ( n m + n m − 1 ) r_{\ell}^{2}=\prod_{m=\ell}^{L}2/(n_m+n_{m-1}) r ℓ 2 = ∏ m = ℓ L 2/ ( n m + n m − 1 ) r L − 1 > r L − 2 > ⋯ > r 1 ≡ r r_{L-1}>r_{L-2}>\cdots>r_1\equiv r r L − 1 > r L − 2 > ⋯ > r 1 ≡ r (69) . These vectors are also non-zero for vanishing biases, but the chain of weights introduces a dependence on γ \gamma γ
G3) The weights of the model for all layers, W ( ℓ ) W^{(\ell)} W ( ℓ ) (84) . The tangent vectors are
∂ ϕ i ′ α ( L ) ∂ W j k ( ℓ ) = ∑ ℓ ′ = 1 ℓ − 1 ∑ p = 1 n ℓ ′ [ W ( L ) ⋯ W ( ℓ + 1 ) ] i ′ j [ W ( ℓ − 1 ) ⋯ W ( ℓ ′ + 1 ) ] k p b p ( ℓ ′ ) + ∑ p = 1 n i n [ W ( L ) ⋯ W ( ℓ + 1 ) ] i ′ j [ W ( ℓ − 1 ) ⋯ W ( 1 ) ] k p X p α = ∑ ℓ ′ = 1 ℓ − 1 r ℓ + 1 γ σ ℓ ∑ p = 1 n ℓ ′ [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] i ′ j [ Ω ( ℓ − 1 ) ⋯ Ω ( ℓ ′ + 1 ) ] k p b p ( ℓ ′ ) + r 1 γ σ ℓ ∑ p = 1 n i n [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] i ′ j [ Ω ( ℓ − 1 ) ⋯ Ω ( 1 ) ] k p X p α , (90) \begin{aligned}
\frac{\partial\phi^{(L)}_{i'\alpha}}{\partial W^{(\ell)}_{jk}}={}&\sum_{\ell'=1}^{\ell-1}\sum_{p=1}^{n_{\ell'}}\left[W^{(L)}\cdots W^{(\ell+1)}\right]_{i'j}\left[W^{(\ell-1)}\cdots W^{(\ell'+1)}\right]_{kp}b^{(\ell')}_p\\
&+\sum_{p=1}^{n_{\mathrm{in}}}\left[W^{(L)}\cdots W^{(\ell+1)}\right]_{i'j}\left[W^{(\ell-1)}\cdots W^{(1)}\right]_{kp}X_{p\alpha}\\
={}&\sum_{\ell'=1}^{\ell-1}\frac{r_{\ell+1}^{\gamma}}{\sigma_{\ell}}\sum_{p=1}^{n_{\ell'}}\left[\Omega^{(L)}\cdots\Omega^{(\ell+1)}\right]_{i'j}\left[\Omega^{(\ell-1)}\cdots\Omega^{(\ell'+1)}\right]_{kp}b^{(\ell')}_p\\
&+\frac{r_1^\gamma}{\sigma_{\ell}}\sum_{p=1}^{n_{\mathrm{in}}}\left[\Omega^{(L)}\cdots\Omega^{(\ell+1)}\right]_{i'j}\left[\Omega^{(\ell-1)}\cdots\Omega^{(1)}\right]_{kp}X_{p\alpha}, \tag*{(90)}
\end{aligned} ∂ W j k ( ℓ ) ∂ ϕ i ′ α ( L ) = = ℓ ′ = 1 ∑ ℓ − 1 p = 1 ∑ n ℓ ′ [ W ( L ) ⋯ W ( ℓ + 1 ) ] i ′ j [ W ( ℓ − 1 ) ⋯ W ( ℓ ′ + 1 ) ] k p b p ( ℓ ′ ) + p = 1 ∑ n in [ W ( L ) ⋯ W ( ℓ + 1 ) ] i ′ j [ W ( ℓ − 1 ) ⋯ W ( 1 ) ] k p X p α ℓ ′ = 1 ∑ ℓ − 1 σ ℓ r ℓ + 1 γ p = 1 ∑ n ℓ ′ [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] i ′ j [ Ω ( ℓ − 1 ) ⋯ Ω ( ℓ ′ + 1 ) ] k p b p ( ℓ ′ ) + σ ℓ r 1 γ p = 1 ∑ n in [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] i ′ j [ Ω ( ℓ − 1 ) ⋯ Ω ( 1 ) ] k p X p α , (90) where we made the dependence on γ \gamma γ (90) cancels out.
We conclude by observing that eq. (84) serves as an example to demonstrate the difference between the span of the tangent plane and the frame induced by the parametrisation. Here, the frame is made of the N p a r a m s = 763 N_{\mathrm{params}}=763 N params = 763 V l i n \mathcal{V}_{\mathrm{lin}} V lin (71) , which is the span. Since the eigenvectors of the NTK are linear combinations of the tangent vectors, they are an aggregation of these directions and are therefore contained in V l i n \mathcal{V}_{\mathrm{lin}} V lin Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) q ( r ) q^{(r)} q ( r ) (27) for the DLN.
The NTK of a deep linear network
In the DLN case the computation of the NTK is analytically tractable. For the sake of simplicity, in the derivation we neglect preprocessing factors. They can be added back in using eq. (85) .
We start by separating the x x x (84) from its coefficients u u u v v v w w w C ∈ R N p a r a m s × 3 N o u t C\in\mathbb{R}^{N_{\mathrm{params}}\times3N_{\mathrm{out}}} C ∈ R N params × 3 N out C C C μ = 1 , … , N p a r a m s \mu=1,\ldots,N_{\mathrm{params}} μ = 1 , … , N params ( r , i ′ ) (r,i') ( r , i ′ ) input-major ordering, the input index r r r
( r , i ′ ) ⟼ N o u t ( r − 1 ) + i ′ , r = 1 , 2 , 3 , i ′ = 1 , … , N o u t , (91) (r,i')\longmapsto N_{\mathrm{out}}(r-1)+i',\qquad r=1,2,3,\quad i'=1,\ldots,N_{\mathrm{out}}, \tag*{(91)} ( r , i ′ ) ⟼ N out ( r − 1 ) + i ′ , r = 1 , 2 , 3 , i ′ = 1 , … , N out , (91) where primed indices label the outputs of the network as in eq. (84) . Note that this is not the flavour-major convention defined in § 2.1 where the pair ( i , α ) (i,\alpha) ( i , α ) (91) the matrix C C C
C = [ u ∣ v ∣ w ] , u , v , w ∈ R N p a r a m s × N o u t , (92) C=\left[u\mid v\mid w\right],\qquad u,v,w\in\mathbb{R}^{N_{\mathrm{params}}\times N_{\mathrm{out}}}, \tag*{(92)} C = [ u ∣ v ∣ w ] , u , v , w ∈ R N params × N out , (92) where each block collects the coefficients for all N o u t N_{\mathrm{out}} N out V l i n \mathcal{V}_{\mathrm{lin}} V lin (71) .
In order to introduce the x x x (84) , we define the matrix of input variables Φ ∈ R N g r i d × 3 \Phi\in\mathbb{R}^{N_{\mathrm{grid}}\times3} Φ ∈ R N grid × 3 V l i n \mathcal{V}_{\mathrm{lin}} V lin
Φ = [ ϕ 1 ∣ ϕ 2 ∣ ϕ 3 ] , ϕ 1 , ϕ 2 , ϕ 3 ∈ R N g r i d , (93) \Phi=\left[\phi_1\mid\phi_2\mid\phi_3\right],\qquad\phi_1,\phi_2,\phi_3\in\mathbb{R}^{N_{\mathrm{grid}}}, \tag*{(93)} Φ = [ ϕ 1 ∣ ϕ 2 ∣ ϕ 3 ] , ϕ 1 , ϕ 2 , ϕ 3 ∈ R N grid , (93) where each ϕ r \phi_r ϕ r r r r V l i n \mathcal{V}_{\mathrm{lin}} V lin x x x Φ \Phi Φ N o u t N_{\mathrm{out}} N out
Ψ ( n e t ) = [ 1 N o u t ⊗ ϕ 1 | 1 N o u t ⊗ ϕ 2 | 1 N o u t ⊗ ϕ 3 ] ∈ R N o u t N g r i d × 3 N o u t , (94) \Psi^{(\mathrm{net})}=\left[\,\mathbf{1}_{N_{\mathrm{out}}}\otimes\phi_1\ \middle|\ \mathbf{1}_{N_{\mathrm{out}}}\otimes\phi_2\ \middle|\ \mathbf{1}_{N_{\mathrm{out}}}\otimes\phi_3\,\right]\in\mathbb{R}^{N_{\mathrm{out}}N_{\mathrm{grid}}\times3N_{\mathrm{out}}}, \tag*{(94)} Ψ ( net ) = [ 1 N out ⊗ ϕ 1 ∣ 1 N out ⊗ ϕ 2 ∣ 1 N out ⊗ ϕ 3 ] ∈ R N out N grid × 3 N out , (94) where each block of eq. (94) is a block-diagonal matrix with elements
1 N o u t ⊗ ϕ r = diag ( ϕ r , … , ϕ r ) N o u t , r = 1 , 2 , 3. (95) \mathbf{1}_{N_{\mathrm{out}}}\otimes\phi_r=\operatorname{diag}(\phi_r,\ldots,\phi_r)_{N_{\mathrm{out}}},\qquad r=1,2,3. \tag*{(95)} 1 N out ⊗ ϕ r = diag ( ϕ r , … , ϕ r ) N out , r = 1 , 2 , 3. (95) Equivalently, adopting the input-major notation of eq. (91) , eq. (94) can be expressed as
Ψ i ′ α , r j ′ ( n e t ) = δ i ′ j ′ Φ α r . (96) \Psi^{(\mathrm{net})}_{i' \alpha,rj'}=\delta_{i'j'}\Phi_{\alpha r}. \tag*{(96)} Ψ i ′ α , r j ′ ( net ) = δ i ′ j ′ Φ α r . (96) The rotation to the N f l a v N_{\mathrm{flav}} N flav R ( F K ) R^{(\mathrm{FK})} R ( FK ) (85) to all points in the grid
R = ( R ( F K ) ) T ⊗ 1 N g r i d ∈ R N v e c × N o u t N g r i d , (97) \mathcal{R}=\left(R^{(\mathrm{FK})}\right)^{T}\otimes\mathbf{1}_{N_{\mathrm{grid}}}\in\mathbb{R}^{N_{\mathrm{vec}}\times N_{\mathrm{out}}N_{\mathrm{grid}}}, \tag*{(97)} R = ( R ( FK ) ) T ⊗ 1 N grid ∈ R N vec × N out N grid , (97) and then by applying the resulting matrix to eq. (96) , which yields
Ψ = R Ψ ( n e t ) = [ ( R ( F K ) ) T ⊗ ϕ 1 | ( R ( F K ) ) T ⊗ ϕ 2 | ( R ( F K ) ) T ⊗ ϕ 3 ] ∈ R N v e c × 3 N o u t , (98) \Psi=\mathcal{R}\Psi^{(\mathrm{net})}=\left[\left(R^{(\mathrm{FK})}\right)^{T}\otimes\phi_1\ \middle|\ \left(R^{(\mathrm{FK})}\right)^{T}\otimes\phi_2\ \middle|\ \left(R^{(\mathrm{FK})}\right)^{T}\otimes\phi_3\right]\in\mathbb{R}^{N_{\mathrm{vec}}\times3N_{\mathrm{out}}}, \tag*{(98)} Ψ = R Ψ ( net ) = [ ( R ( FK ) ) T ⊗ ϕ 1 ( R ( FK ) ) T ⊗ ϕ 2 ( R ( FK ) ) T ⊗ ϕ 3 ] ∈ R N vec × 3 N out , (98) where each block follows from the mixed-product rule
( A ⊗ B ) ( C ⊗ D ) = ( A C ) ⊗ ( B D ) , (99) (A\otimes B)(C\otimes D)=(AC)\otimes(BD), \tag*{(99)} ( A ⊗ B ) ( C ⊗ D ) = ( A C ) ⊗ ( B D ) , (99) and can be seen as the application of the rotation matrix R ( F K ) R^{(\mathrm{FK})} R ( FK ) (95) .
With this notation, it can be easily verified that the tangent vectors factorise into the following product
J = Ψ C T = R Ψ ( n e t ) C T , (100) J=\Psi C^{T}=\mathcal{R}\Psi^{(\mathrm{net})}C^{T}, \tag*{(100)} J = Ψ C T = R Ψ ( net ) C T , (100) and the NTK is then expressed in the convenient form
Θ = J J T = Ψ C T C Ψ T = Ψ G Ψ T . (101) \Theta=JJ^{T}=\Psi C^{T}C\Psi^{T}=\Psi\mathcal{G}\Psi^{T}. \tag*{(101)} Θ = J J T = Ψ C T C Ψ T = Ψ G Ψ T . (101) Here G = C T C ∈ R 3 N o u t × 3 N o u t \mathcal{G}=C^{T}C\in\mathbb{R}^{3N_{\mathrm{out}}\times3N_{\mathrm{out}}} G = C T C ∈ R 3 N out × 3 N out (92) contracted over the parameter index. With the vectorisation convention adopted in eq. (92) , we can express this matrix as
G ≡ [ G 11 G 12 G 13 G 21 G 22 G 23 G 31 G 32 G 33 ] = [ u T v T w T ] [ u v w ] = [ u T u u T v u T w v T u v T v v T w w T u w T v w T w ] , G r s ∈ R N o u t × N o u t , (102) \mathcal{G}\equiv
\begin{bmatrix}
\mathcal{G}_{11} & \mathcal{G}_{12} & \mathcal{G}_{13}\\
\mathcal{G}_{21} & \mathcal{G}_{22} & \mathcal{G}_{23}\\
\mathcal{G}_{31} & \mathcal{G}_{32} & \mathcal{G}_{33}
\end{bmatrix}
=
\begin{bmatrix}
u^{T}\\
v^{T}\\
w^{T}
\end{bmatrix}
\begin{bmatrix}
u & v & w
\end{bmatrix}
=
\begin{bmatrix}
u^{T}u & u^{T}v & u^{T}w\\
v^{T}u & v^{T}v & v^{T}w\\
w^{T}u & w^{T}v & w^{T}w
\end{bmatrix},
\qquad\mathcal{G}_{rs}\in\mathbb{R}^{N_{\mathrm{out}}\times N_{\mathrm{out}}}, \tag*{(102)} G ≡ G 11 G 21 G 31 G 12 G 22 G 32 G 13 G 23 G 33 = u T v T w T [ u v w ] = u T u v T u w T u u T v v T v w T v u T w v T w w T w , G r s ∈ R N out × N out , (102) which ultimately motivates our choice of the input-major ordering: taking the input index as the slowest is what turns G \mathcal{G} G 3 × 3 3\times3 3 × 3 N o u t N_{\mathrm{out}} N out V l i n \mathcal{V}_{\mathrm{lin}} V lin
By introducing an auxiliary vector of coefficients ζ ∈ R 3 N o u t \zeta\in\mathbb{R}^{3N_{\mathrm{out}}} ζ ∈ R 3 N out (91) ,
ζ = ( ζ r i ′ ) , r = 1 , 2 , 3 , i ′ = 1 , … , N o u t , (103) \zeta=(\zeta_{ri'}),\qquad r=1,2,3,\qquad i'=1,\ldots,N_{\mathrm{out}}, \tag*{(103)} ζ = ( ζ r i ′ ) , r = 1 , 2 , 3 , i ′ = 1 , … , N out , (103) and writing the NTK eigenvectors as z = Ψ ζ z=\Psi\zeta z = Ψ ζ Θ z = λ z \Theta z=\lambda z Θ z = λ z
G S ζ = λ ζ , (104) \mathcal{G}S\zeta=\lambda\zeta, \tag*{(104)} G S ζ = λ ζ , (104) where we used eq. (99) and eq. (86) to define
S = Ψ T Ψ = Φ T Φ ⊗ R ( F K ) ( R ( F K ) ) T = S 3 ⊗ diag ( m i ′ ) N o u t ∈ R 3 N o u t × 3 N o u t , (105) S=\Psi^{T}\Psi=\Phi^{T}\Phi\otimes R^{(\mathrm{FK})}\left(R^{(\mathrm{FK})}\right)^{T}=S_3\otimes\operatorname{diag}(m_{i'})_{N_{\mathrm{out}}}\in\mathbb{R}^{3N_{\mathrm{out}}\times3N_{\mathrm{out}}}, \tag*{(105)} S = Ψ T Ψ = Φ T Φ ⊗ R ( FK ) ( R ( FK ) ) T = S 3 ⊗ diag ( m i ′ ) N out ∈ R 3 N out × 3 N out , (105) where m i ′ m_{i'} m i ′ (86) and S 3 S_3 S 3
S 3 = Φ T Φ = N g r i d ( 1 ⟨ x ⟩ ⟨ log x ⟩ ⟨ x ⟩ ⟨ x 2 ⟩ ⟨ x log x ⟩ ⟨ log x ⟩ ⟨ x log x ⟩ ⟨ log 2 x ⟩ ) , ⟨ g ( x ) ⟩ ≡ 1 N g r i d ∑ α = 1 N g r i d g ( x α ) . (106) S_3=\Phi^{T}\Phi=N_{\mathrm{grid}}
\begin{pmatrix}
1 & \langle x\rangle& \langle\log x\rangle\\
\langle x\rangle& \langle x^2\rangle& \langle x\log x\rangle\\
\langle\log x\rangle& \langle x\log x\rangle& \langle\log^2x\rangle
\end{pmatrix},
\qquad\langle g(x)\rangle\equiv\frac{1}{N_{\mathrm{grid}}}\sum_{\alpha=1}^{N_{\mathrm{grid}}}g(x_\alpha). \tag*{(106)} S 3 = Φ T Φ = N grid 1 ⟨ x ⟩ ⟨ log x ⟩ ⟨ x ⟩ ⟨ x 2 ⟩ ⟨ x log x ⟩ ⟨ log x ⟩ ⟨ x log x ⟩ ⟨ log 2 x ⟩ , ⟨ g ( x )⟩ ≡ N grid 1 α = 1 ∑ N grid g ( x α ) . (106) Note that, similarly to G \mathcal{G} G S S S 3 × 3 3\times3 3 × 3 N o u t × N o u t N_{\mathrm{out}}\times N_{\mathrm{out}} N out × N out ( r , s ) (r,s) ( r , s )
S r s = ( S 3 ) r s diag ( m i ′ ) N o u t . (107) S_{rs}=(S_3)_{rs}\operatorname{diag}(m_{i'})_{N_{\mathrm{out}}}. \tag*{(107)} S r s = ( S 3 ) r s diag ( m i ′ ) N out . (107) We pause here to make a few comments.
Factorisation. The NTK in eq. (101) factorises into a quadratic form where Ψ \Psi Ψ Φ \Phi Φ R ( F K ) R(\mathrm{FK}) R ( FK ) G \mathcal{G} G Ψ \Psi Ψ G \mathcal{G} G 3 N o u t 3N_{\mathrm{out}} 3 N out G \mathcal{G} G 3 N o u t ( 3 N o u t + 1 ) / 2 = 300 3N_{\mathrm{out}}(3N_{\mathrm{out}}+1)/2=300 3 N out ( 3 N out + 1 ) /2 = 300 N p a r a m s = 763 N_{\mathrm{params}}=763 N params = 763
Reduced NTK. The N v e c N_{\mathrm{vec}} N vec (104) . The non-zero eigenvalues – and the corresponding eigenvectors – are then completely specified by the 3 × N o u t = 24 3\times N_{\mathrm{out}}=24 3 × N out = 24 (104) . Note that although G S \mathcal{G}S G S S 1 / 2 G S 1 / 2 S^{1/2}\mathcal{G}S^{1/2} S 1/2 G S 1/2 ζ \zeta ζ S S S z = Ψ ζ z=\Psi\zeta z = Ψ ζ 24 × 24 24\times24 24 × 24 N v e c × N v e c N_{\mathrm{vec}}\times N_{\mathrm{vec}} N vec × N vec
Different set of input functions. The expressiveness of a DLN is uniquely determined by the input functions that make up Φ \Phi Φ x x x log x \log x log x 21 (93) to arbitrary n 1 n_1 n 1
Φ = [ ϕ 1 ∣ ϕ 2 ∣ ⋯ ∣ ϕ n 1 ] , ϕ r ∈ R N g r i d , r = 1 , ⋯ , n 1 , (108) \Phi=\left[\phi_1\mid\phi_2\mid\cdots\mid\phi_{n_1}\right],\qquad\phi_r\in\mathbb{R}^{N_{\mathrm{grid}}},\qquad r=1,\cdots,n_1, \tag*{(108)} Φ = [ ϕ 1 ∣ ϕ 2 ∣ ⋯ ∣ ϕ n 1 ] , ϕ r ∈ R N grid , r = 1 , ⋯ , n 1 , (108) where ϕ r \phi_r ϕ r r r r x x x e.g. , polynomials). The DLN – and the corresponding tangent vectors – then become equivalent to any fixed functional form that is a linear combination of basis functions weighted by free parameters, that is
f ( x ; θ ) = ∑ r = 1 n 1 θ r Φ ( r ) ( x ) , (109) f(x;\theta)=\sum_{r=1}^{n_1}\theta_r\Phi^{(r)}(x), \tag*{(109)} f ( x ; θ ) = r = 1 ∑ n 1 θ r Φ ( r ) ( x ) , (109) which is very similar to the expansion used in Ref. [39 ] . The only difference between the DLN case and eq. (109) is that, in the former, there are many equivalent combinations of weights and biases that realise the same weight θ r \theta_r θ r V l i n \mathcal{V}_{\mathrm{lin}} V lin
The reduced matrix G \mathcal{G} G
Turning G \mathcal{G} G S S S (104) in the following form
\begin{bmatrix}\mathcal{G}_{cc} & \mathcal{G}_{cn}\\ \mathcal{G}_{nc} & \mathcal{G}_{nn}\end{bmatrix}
\begin{bmatrix}S_{cc} & S_{cn}\\ S_{nc} & S_{nn}\end{bmatrix}
\begin{bmatrix}\zeta_c\\ \zeta_n\end{bmatrix}
=\lambda\begin{bmatrix}\zeta_c\\ \zeta_n\end{bmatrix}, \tag*{(110)}\end{bmatrix}
where the blocks of ζ \zeta ζ S S S G \mathcal{G} G
ζ c = ζ 1 , ζ n = [ ζ 2 ζ 3 ] , (111) \zeta_c=\zeta_1,\qquad\zeta_n=\begin{bmatrix}\zeta_2\\ \zeta_3\end{bmatrix}, \tag*{(111)} ζ c = ζ 1 , ζ n = [ ζ 2 ζ 3 ] , (111) S c c = S 11 , S n n = [ S 22 S 23 S 32 S 33 ] , S c n = [ S 12 S 13 ] ≡ S n c T , (112) S_{cc}=S_{11},\qquad S_{nn}=\begin{bmatrix}S_{22} & S_{23}\\ S_{32} & S_{33}\end{bmatrix},\qquad S_{cn}=\begin{bmatrix}S_{12} & S_{13}\end{bmatrix}\equiv S_{nc}^{T}, \tag*{(112)} S cc = S 11 , S nn = [ S 22 S 32 S 23 S 33 ] , S c n = [ S 12 S 13 ] ≡ S n c T , (112) G c c = G 11 , G n n = [ G 22 G 23 G 32 G 33 ] , G c n = [ G 12 G 13 ] ≡ G n c T . (113) \mathcal{G}_{cc}=\mathcal{G}_{11},\qquad\mathcal{G}_{nn}=\begin{bmatrix}\mathcal{G}_{22} & \mathcal{G}_{23}\\ \mathcal{G}_{32} & \mathcal{G}_{33}\end{bmatrix},\qquad\mathcal{G}_{cn}=\begin{bmatrix}\mathcal{G}_{12} & \mathcal{G}_{13}\end{bmatrix}\equiv\mathcal{G}_{nc}^{T}. \tag*{(113)} G cc = G 11 , G nn = [ G 22 G 32 G 23 G 33 ] , G c n = [ G 12 G 13 ] ≡ G n c T . (113) We now compute the entries of G \mathcal{G} G Off-diagonal terms. The parameter index runs over the three classes of parameters defined above – last- and hidden-layer biases and the weights – so that each of the three matrices in eq. (92) splits by rows as
u = [ u ( b L ) u ( b h ) u ( w ) ] , and likewise for v and w , (114) u=\begin{bmatrix}u^{(\mathrm{bL})}\\u^{(\mathrm{bh})}\\u^{(\mathrm{w})}\end{bmatrix},
\qquad\text{and likewise for }v\text{ and }w, \tag*{(114)} u = u ( bL ) u ( bh ) u ( w ) , and likewise for v and w , (114) the blocks having N o u t N_{\mathrm{out}} N out ∑ ℓ < L n ℓ \sum_{\ell<L}n_\ell ∑ ℓ < L n ℓ ∑ ℓ n ℓ n ℓ − 1 \sum_\ell n_\ell n_{\ell-1} ∑ ℓ n ℓ n ℓ − 1 G = C T C \mathcal{G}=C^{T}C G = C T C G \mathcal{G} G
G 21 = v T u = ( v ( b L ) ) T u ( b L ) + ( v ( b h ) ) T u ( b h ) + ( v ( w ) ) T u ( w ) . (115) \mathcal{G}_{21}=v^{T}u=\left(v^{(\mathrm{bL})}\right)^{T}u^{(\mathrm{bL})}+\left(v^{(\mathrm{bh})}\right)^{T}u^{(\mathrm{bh})}+\left(v^{(\mathrm{w})}\right)^{T}u^{(\mathrm{w})}. \tag*{(115)} G 21 = v T u = ( v ( bL ) ) T u ( bL ) + ( v ( bh ) ) T u ( bh ) + ( v ( w ) ) T u ( w ) . (115) From G1 and G2 , the last- and hidden-layer biases can only generate constant terms, implying v ( b L ) = v ( b h ) = 0 v^{(\mathrm{bL})}=v^{(\mathrm{bh})}=0 v ( bL ) = v ( bh ) = 0 (115) survives. However, since u ( w ) u^{(\mathrm{w})} u ( w ) (90) , when biases are set to zero at initialisation this term also disappears. We can then conclude that
G 12 = G 21 = 0 , if b ( ℓ ) = 0 ∀ ℓ . (116) \mathcal{G}_{12}=\mathcal{G}_{21}=0,\qquad\text{if }b^{(\ell)}=0\ \forall\ell. \tag*{(116)} G 12 = G 21 = 0 , if b ( ℓ ) = 0 ∀ ℓ . (116) Exactly the same argument applies to G 13 = w T u \mathcal{G}_{13}=w^{T}u G 13 = w T u w w w
G c n = G n c T = 0. (117) \mathcal{G}_{\mathrm{cn}}=\mathcal{G}_{\mathrm{nc}}^{T}=0. \tag*{(117)} G cn = G nc T = 0. (117) Constant term. The same decomposition applied to the diagonal block of the constant sector gives
G c c = G 11 = u T u = ( u ( b L ) ) T u ( b L ) + ( u ( b h ) ) T u ( b h ) + ( u ( w ) ) T u ( w ) , (118) \mathcal{G}_{\mathrm{cc}}=\mathcal{G}_{11}=u^{T}u=\left(u^{(\mathrm{bL})}\right)^{T}u^{(\mathrm{bL})}+\left(u^{(\mathrm{bh})}\right)^{T}u^{(\mathrm{bh})}+\left(u^{(\mathrm{w})}\right)^{T}u^{(\mathrm{w})}, \tag*{(118)} G cc = G 11 = u T u = ( u ( bL ) ) T u ( bL ) + ( u ( bh ) ) T u ( bh ) + ( u ( w ) ) T u ( w ) , (118) where the last term vanishes at initialisation for the reason given above, provided that biases are set to zero. The first one is fixed by G1 : eq. (88) gives u ( b L ) = 1 N o u t u^{(\mathrm{bL})}=\mathbf{1}_{N_{\mathrm{out}}} u ( bL ) = 1 N out N o u t N_{\mathrm{out}} N out γ \gamma γ G2 and eq. (89) , each hidden layer contributing the Gram matrix of its own chain. Altogether
G c c = 1 N o u t + ∑ ℓ = 1 L − 1 r ℓ + 1 2 γ [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] T . (119) \mathcal{G}_{\mathrm{cc}}=\mathbf{1}_{N_{\mathrm{out}}}+\sum_{\ell=1}^{L-1}r_{\ell+1}^{2\gamma}\left[\Omega^{(L)}\cdots\Omega^{(\ell+1)}\right]\left[\Omega^{(L)}\cdots\Omega^{(\ell+1)}\right]^{T}. \tag*{(119)} G cc = 1 N out + ℓ = 1 ∑ L − 1 r ℓ + 1 2 γ [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] [ Ω ( L ) ⋯ Ω ( ℓ + 1 ) ] T . (119) The correction is exponentially small in γ \gamma γ ℓ = L − 1 \ell=L-1 ℓ = L − 1 n L − 1 r L 2 γ n_{L-1}r_{L}^{2\gamma} n L − 1 r L 2 γ
Non-constant terms. The blocks with r , s ≠ 1 r,s\ne1 r , s = 1 v v v w w w (90) survives. Denoting by
Ω > ℓ = Ω ( L ) ⋯ Ω ( ℓ + 1 ) , Ω < ℓ = Ω ( ℓ − 1 ) ⋯ Ω ( 1 ) , (120) \Omega_{>\ell}=\Omega^{(L)}\cdots\Omega^{(\ell+1)},\qquad\Omega_{<\ell}=\Omega^{(\ell-1)}\cdots\Omega^{(1)}, \tag*{(120)} Ω > ℓ = Ω ( L ) ⋯ Ω ( ℓ + 1 ) , Ω < ℓ = Ω ( ℓ − 1 ) ⋯ Ω ( 1 ) , (120) the two chains of rescaled weights that the derivative leaves above and below the differentiated layer ℓ \ell ℓ
G r s = ∑ ℓ = 1 L r ℓ 2 γ σ ℓ 2 ( Ω < ℓ T Ω < ℓ ) r − 1 , s − 1 Ω > ℓ Ω > ℓ T , r , s = 2 , 3 , (121) \mathcal{G}_{rs}=\sum_{\ell=1}^{L}\frac{r_{\ell}^{2\gamma}}{\sigma_{\ell}^{2}}\left(\Omega_{<\ell}^{T}\Omega_{<\ell}\right)_{r-1,s-1}\Omega_{>\ell}\Omega_{>\ell}^{T},\qquad r,s=2,3, \tag*{(121)} G r s = ℓ = 1 ∑ L σ ℓ 2 r ℓ 2 γ ( Ω < ℓ T Ω < ℓ ) r − 1 , s − 1 Ω > ℓ Ω > ℓ T , r , s = 2 , 3 , (121) the shift by one accounting for the constant direction r = 1 r=1 r = 1 X p α X_{p\alpha} X p α
r ℓ 2 γ σ ℓ 2 = ∏ m ≠ ℓ ρ m γ , ρ m = 2 n m + n m − 1 , (122) \frac{r_{\ell}^{2\gamma}}{\sigma_{\ell}^{2}}=\prod_{m\ne\ell}\rho_m^{\gamma},\qquad\rho_m=\frac{2}{n_m+n_{m-1}}, \tag*{(122)} σ ℓ 2 r ℓ 2 γ = m = ℓ ∏ ρ m γ , ρ m = n m + n m − 1 2 , (122) that is, the product of all the layer variances but the one being differentiated. Every entry of G n n \mathcal{G}_{\mathrm{nn}} G nn L − 1 L-1 L − 1 ρ γ \rho^\gamma ρ γ (119) .
The matrix G \mathcal{G} G Collecting the three results, at initialisation G \mathcal{G} G
G c c = 1 N o u t + O ( r L 2 γ ) , G c n = 0 , G n n = O ( ∏ m ≠ ℓ ρ m γ ) . (123) \mathcal{G}_{cc} = \mathbb{1}_{N_{\mathrm{out}}} + \mathcal{O}\left(r_L^{2\gamma}\right), \qquad\mathcal{G}_{cn} = 0, \qquad\mathcal{G}_{nn} = \mathcal{O}\left(\prod_{m\ne\ell}\rho_m^\gamma\right). \tag*{(123)} G cc = 1 N out + O ( r L 2 γ ) , G c n = 0 , G nn = O m = ℓ ∏ ρ m γ . (123) The two sectors are thus separated by the initialisation width alone: as γ \gamma γ G c c \mathcal{G}_{cc} G cc G n n \mathcal{G}_{nn} G nn (110) then factorises as
G S = [ G c c 0 0 G n n ] [ S c c S c n S n c S n n ] = [ G c c S c c G c c S c n G n n S n c G n n S n n ] . (124) \mathcal{G}S =
\begin{bmatrix}
\mathcal{G}_{cc} & 0 \\
0 & \mathcal{G}_{nn}
\end{bmatrix}
\begin{bmatrix}
S_{cc} & S_{cn} \\
S_{nc} & S_{nn}
\end{bmatrix}
=
\begin{bmatrix}
\mathcal{G}_{cc}S_{cc} & \mathcal{G}_{cc}S_{cn} \\
\mathcal{G}_{nn}S_{nc} & \mathcal{G}_{nn}S_{nn}
\end{bmatrix}. \tag*{(124)} G S = [ G cc 0 0 G nn ] [ S cc S n c S c n S nn ] = [ G cc S cc G nn S n c G cc S c n G nn S nn ] . (124) The two diagonal blocks are of order unity and of order ∏ m ≠ ℓ ρ m γ \prod_{m\ne\ell}\rho_m^\gamma ∏ m = ℓ ρ m γ γ \gamma γ S c n ≠ 0 S_{cn}\ne0 S c n = 0
( G S ) n c = G n n S n c = O ( ∏ m ≠ ℓ ρ m γ ) , (125) (\mathcal{G}S)_{nc}=\mathcal{G}_{nn}S_{nc}=\mathcal{O}\left(\prod_{m\ne\ell}\rho_m^\gamma\right), \tag*{(125)} ( G S ) n c = G nn S n c = O m = ℓ ∏ ρ m γ , (125) which is negligible against the block of order unity sitting above it. To that accuracy G S \mathcal{G}S G S N o u t N_{\mathrm{out}} N out
spec ( G S ) ∣ O ( 1 ) = spec ( G c c S c c ) . (126) \left.\operatorname{spec}(\mathcal{G}S)\right|_{\mathcal{O}(1)}=\operatorname{spec}(\mathcal{G}_{cc}S_{cc}). \tag*{(126)} spec ( G S ) ∣ O ( 1 ) = spec ( G cc S cc ) . (126) The same argument does not apply to the remaining 2 N o u t 2N_{\mathrm{out}} 2 N out (125) carries one factor of G n n \mathcal{G}_{nn} G nn ζ c \zeta_c ζ c
The constant sector
By eq. (126) the N o u t N_{\mathrm{out}} N out G c c S c c \mathcal{G}_{cc}S_{cc} G cc S cc (107) for S c c = N g r i d diag ( m i ′ ) S_{cc}=N_{\mathrm{grid}}\operatorname{diag}(m_{i'}) S cc = N grid diag ( m i ′ ) (119) for G c c \mathcal{G}_{cc} G cc
λ i ′ ( c ) = m i ′ N g r i d [ 1 + O ( r L 2 γ ) ] , i ′ = 1 , … , N o u t , (127) \lambda_{i'}^{(c)}=m_{i'}N_{\mathrm{grid}}\left[1+\mathcal{O}\left(r_L^{2\gamma}\right)\right], \qquad i'=1,\ldots,N_{\mathrm{out}}, \tag*{(127)} λ i ′ ( c ) = m i ′ N grid [ 1 + O ( r L 2 γ ) ] , i ′ = 1 , … , N out , (127) that is, in the limit of vanishing initialisation width,
{ λ i ′ ( c ) } = { m i ′ N g r i d } i ′ = 1 , … , N o u t = { 200 , 150 , 50 , 50 , 50 , 50 , 50 , 50 } , (128) \left\{\lambda_{i'}^{(c)}\right\}=\left\{m_{i'}N_{\mathrm{grid}}\right\}_{i'=1,\ldots,N_{\mathrm{out}}}=\{200,150,50,50,50,50,50,50\}, \tag*{(128)} { λ i ′ ( c ) } = { m i ′ N grid } i ′ = 1 , … , N out = { 200 , 150 , 50 , 50 , 50 , 50 , 50 , 50 } , (128) which reproduces, once normalised to the largest eigenvalue, the first tier of the spectrum observed in Figure 29 . The multiplicities m i ′ m_{i'} m i ′ (86) are therefore the only structure in this tier: they come from the flavour map R ( F K ) R^{(\mathrm{FK})} R ( FK )
Figure 29. The spectrum of Θ \Theta Θ r r r γ \gamma γ λ r / λ r + 1 \lambda_r/\lambda_{r+1} λ r / λ r + 1 N r e p = 120 N_{\mathrm{rep}}=120 N rep = 120 N o u t N_{\mathrm{out}} N out 2 N o u t 2N_{\mathrm{out}} 2 N out 3 N o u t 3N_{\mathrm{out}} 3 N out (71) . Note that the curve for each γ \gamma γ dim Im ( Θ ) \dim\operatorname{Im}(\Theta) dim Im ( Θ ) Ker ( Θ ) \operatorname{Ker}(\Theta) Ker ( Θ )
The eigenvectors follow at once. To leading order G c c = 1 N o u t \mathcal{G}_{cc}=\mathbb{1}_{N_{\mathrm{out}}} G cc = 1 N out S c c S_{cc} S cc ζ c \zeta_c ζ c i ′ i' i ′ ζ n = 0 \zeta_n=0 ζ n = 0 z = Ψ ζ z=\Psi\zeta z = Ψ ζ (98) , the normalised modes are
z α ( i ′ ) = ( R ( F K ) ) i ′ i m i ′ N g r i d , (129) z_{\alpha}^{(i')}=\frac{\left(R^{(\mathrm{FK})}\right)_{i'i}}{\sqrt{m_{i'}N_{\mathrm{grid}}}}, \tag*{(129)} z α ( i ′ ) = m i ′ N grid ( R ( FK ) ) i ′ i , (129) independent of the grid point α \alpha α x x x ( z ( i ′ ) , z ( j ′ ) ) = ζ ( i ′ ) T S ζ ( j ′ ) \left(z^{(i')},z^{(j')}\right)=\zeta^{(i')T}S\zeta^{(j')} ( z ( i ′ ) , z ( j ′ ) ) = ζ ( i ′ ) T S ζ ( j ′ )
( z ( i ′ ) , z ( j ′ ) ) = ( S c c ) i ′ j ′ m i ′ N g r i d = δ i ′ j ′ . (130) \left(z^{(i')},z^{(j')}\right)=\frac{(S_{cc})_{i'j'}}{m_{i'}N_{\mathrm{grid}}}=\delta_{i'j'}. \tag*{(130)} ( z ( i ′ ) , z ( j ′ ) ) = m i ′ N grid ( S cc ) i ′ j ′ = δ i ′ j ′ . (130) The non-constant sector
The remaining 2 N o u t 2N_{\mathrm{out}} 2 N out (110) , G n n ( S n c ζ c + S n n ζ n ) = λ ζ n \mathcal{G}_{nn}(S_{nc}\zeta_c+S_{nn}\zeta_n)=\lambda\zeta_n G nn ( S n c ζ c + S nn ζ n ) = λ ζ n G n n \mathcal{G}_{nn} G nn λ = O ( G n n ) \lambda=\mathcal{O}(\mathcal{G}_{nn}) λ = O ( G nn ) (128) . The first block row can then be solved for ζ c \zeta_c ζ c G c c ( S c c ζ c + S c n ζ n ) = λ ζ c \mathcal{G}_{cc}(S_{cc}\zeta_c+S_{cn}\zeta_n)=\lambda\zeta_c G cc ( S cc ζ c + S c n ζ n ) = λ ζ c G c c S c c = O ( 1 ) \mathcal{G}_{cc}S_{cc}=\mathcal{O}(1) G cc S cc = O ( 1 ) G c c \mathcal{G}_{cc} G cc
ζ c = − S c c − 1 S c n ζ n . (131) \zeta_c=-S_{cc}^{-1}S_{cn}\zeta_n . \tag*{(131)} ζ c = − S cc − 1 S c n ζ n . (131) The two components are therefore of the same order, as anticipated. Substituting eq. (131) into the second block row leaves a problem of dimension 2 N o u t 2N_{\mathrm{out}} 2 N out
G n n S ^ ζ n = λ ζ n , (132) \mathcal{G}_{nn}\widehat{S}\zeta_n=\lambda\zeta_n , \tag*{(132)} G nn S ζ n = λ ζ n , (132) where the elimination has replaced S n n S_{nn} S nn (107) ,
S ^ = S n n − S n c S c c − 1 S c n = S ^ 2 ⊗ diag ( m i ′ ) N o u t , S ^ 2 = N g r i d ( ⟨ x 2 ⟩ − ⟨ x ⟩ 2 ⟨ x log x ⟩ − ⟨ x ⟩ ⟨ log x ⟩ ⟨ x log x ⟩ − ⟨ x ⟩ ⟨ log x ⟩ ⟨ log 2 x ⟩ − ⟨ log x ⟩ 2 ) . (133) \widehat{S}=S_{nn}-S_{nc}S_{cc}^{-1}S_{cn}=\widehat{S}_2\otimes\operatorname{diag}(m_i')_{N_{\mathrm{out}}},\qquad
\widehat{S}_2=N_{\mathrm{grid}}
\begin{pmatrix}
\langle x^2\rangle-\langle x\rangle^2 & \langle x\log x\rangle-\langle x\rangle\langle\log x\rangle\\
\langle x\log x\rangle-\langle x\rangle\langle\log x\rangle& \langle\log^2x\rangle-\langle\log x\rangle^2
\end{pmatrix}. \tag*{(133)} S = S nn − S n c S cc − 1 S c n = S 2 ⊗ diag ( m i ′ ) N out , S 2 = N grid ( ⟨ x 2 ⟩ − ⟨ x ⟩ 2 ⟨ x log x ⟩ − ⟨ x ⟩ ⟨ log x ⟩ ⟨ x log x ⟩ − ⟨ x ⟩ ⟨ log x ⟩ ⟨ log 2 x ⟩ − ⟨ log x ⟩ 2 ) . (133) The grid moments of eq. (106) are thus replaced by centred moments.
Recall that the blocks of G n n \mathcal{G}_{nn} G nn (121) . Collected into a single matrix, each layer contributes one Kronecker product of a 2 × 2 2\times2 2 × 2 N o u t × N o u t N_{\mathrm{out}}\times N_{\mathrm{out}} N out × N out
G n n = r 2 γ ∑ ℓ = 1 L 1 σ ℓ 2 A ℓ ⊗ B ℓ , A ℓ ≡ Ω < ℓ T Ω < ℓ ∈ R 2 × 2 , B ℓ ≡ Ω > ℓ Ω > ℓ T ∈ R N o u t × N o u t . (134) \mathcal{G}_{nn}=r^{2\gamma}\sum_{\ell=1}^{L}\frac{1}{\sigma_\ell^2}A_\ell\otimes B_\ell,\qquad
A_\ell\equiv\Omega_{<\ell}^{T}\Omega_{<\ell}\in\mathbb{R}^{2\times2},\qquad
B_\ell\equiv\Omega_{>\ell}\Omega_{>\ell}^{T}\in\mathbb{R}^{N_{\mathrm{out}}\times N_{\mathrm{out}}}. \tag*{(134)} G nn = r 2 γ ℓ = 1 ∑ L σ ℓ 2 1 A ℓ ⊗ B ℓ , A ℓ ≡ Ω < ℓ T Ω < ℓ ∈ R 2 × 2 , B ℓ ≡ Ω > ℓ Ω > ℓ T ∈ R N out × N out . (134) Note that the matrix S ^ \widehat{S} S (132) has the same tensor structure as each term of the sum in eq. (134) , but the sum over ℓ \ell ℓ A ℓ A_\ell A ℓ B ℓ B_\ell B ℓ
A ℓ + 1 = Ω < ℓ T ( Ω ( ℓ ) T Ω ( ℓ ) ) Ω < ℓ , B ℓ − 1 = Ω > ℓ ( Ω ( ℓ ) Ω ( ℓ ) T ) Ω > ℓ T . (135) A_{\ell+1}=\Omega_{<\ell}^{T}\left(\Omega^{(\ell)T}\Omega^{(\ell)}\right)\Omega_{<\ell},\qquad
B_{\ell-1}=\Omega_{>\ell}\left(\Omega^{(\ell)}\Omega^{(\ell)T}\right)\Omega_{>\ell}^{T}. \tag*{(135)} A ℓ + 1 = Ω < ℓ T ( Ω ( ℓ ) T Ω ( ℓ ) ) Ω < ℓ , B ℓ − 1 = Ω > ℓ ( Ω ( ℓ ) Ω ( ℓ ) T ) Ω > ℓ T . (135) Since Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ )
E [ A ℓ + 1 ∣ Ω < ℓ ] = Ω < ℓ T E [ Ω ( ℓ ) T Ω ( ℓ ) ] Ω < ℓ = n ℓ A ℓ . (136) \mathbb{E}\left[A_{\ell+1}\mid\Omega_{<\ell}\right]=\Omega_{<\ell}^{T}\mathbb{E}\left[\Omega^{(\ell)T}\Omega^{(\ell)}\right]\Omega_{<\ell}=n_\ell A_\ell. \tag*{(136)} E [ A ℓ + 1 ∣ Ω < ℓ ] = Ω < ℓ T E [ Ω ( ℓ ) T Ω ( ℓ ) ] Ω < ℓ = n ℓ A ℓ . (136) where in the second equality we used the fact that the entries of Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ ) Δ A ℓ + 1 ≡ A ℓ + 1 − n ℓ A ℓ \Delta A_{\ell+1}\equiv A_{\ell+1}-n_\ell A_\ell Δ A ℓ + 1 ≡ A ℓ + 1 − n ℓ A ℓ Ω ( ℓ ) \Omega^{(\ell)} Ω ( ℓ )
E [ ( Δ A ℓ + 1 ) a b 2 ∣ Ω < ℓ ] = n ℓ [ ( A ℓ ) a a ( A ℓ ) b b + ( A ℓ ) a b 2 ] . (137) \mathbb{E}\left[(\Delta A_{\ell+1})_{ab}^{2}\mid\Omega_{<\ell}\right]=n_\ell\left[(A_\ell)_{aa}(A_\ell)_{bb}+(A_\ell)_{ab}^{2}\right]. \tag*{(137)} E [ ( Δ A ℓ + 1 ) ab 2 ∣ Ω < ℓ ] = n ℓ [ ( A ℓ ) aa ( A ℓ ) bb + ( A ℓ ) ab 2 ] . (137) Compared with the average n ℓ A ℓ n_\ell A_\ell n ℓ A ℓ n ℓ − 1 / 2 n_\ell^{-1/2} n ℓ − 1/2 Ω > ℓ \Omega_{>\ell} Ω > ℓ E [ B ℓ − 1 ∣ Ω > ℓ ] = n ℓ − 1 B ℓ \mathbb{E}[B_{\ell-1}\mid\Omega_{>\ell}]=n_{\ell-1}B_\ell E [ B ℓ − 1 ∣ Ω > ℓ ] = n ℓ − 1 B ℓ Δ B ℓ − 1 ≡ B ℓ − 1 − n ℓ − 1 B ℓ \Delta B_{\ell-1}\equiv B_{\ell-1}-n_{\ell-1}B_\ell Δ B ℓ − 1 ≡ B ℓ − 1 − n ℓ − 1 B ℓ n ℓ − 1 − 1 / 2 n_{\ell-1}^{-1/2} n ℓ − 1 − 1/2 B ℓ B_\ell B ℓ N o u t 2 N_{\mathrm{out}}^2 N out 2 Δ B ℓ − 1 \Delta B_{\ell-1} Δ B ℓ − 1 ( N o u t / n ℓ − 1 ) 1 / 2 (N_{\mathrm{out}}/n_{\ell-1})^{1/2} ( N out / n ℓ − 1 ) 1/2
Iterating the recursions from A 1 = 1 2 A_1=\mathbb{1}_2 A 1 = 1 2 B L = 1 N o u t B_L=\mathbb{1}_{N_{\mathrm{out}}} B L = 1 N out
A ^ ℓ ≡ A ℓ a ℓ = 1 2 + δ A ℓ , δ A ℓ = ∑ ℓ ′ = 2 ℓ Δ A ℓ ′ a ℓ ′ , a ℓ = ∏ ℓ ′ = 1 ℓ n ℓ ′ , (138) \widehat{A}_\ell\equiv\frac{A_\ell}{a_\ell}=\mathbb{1}_2+\delta A_\ell,\qquad
\delta A_\ell=\sum_{\ell'=2}^{\ell}\frac{\Delta A_{\ell'}}{a_{\ell'}},\qquad
a_\ell=\prod_{\ell'=1}^{\ell}n_{\ell'}, \tag*{(138)} A ℓ ≡ a ℓ A ℓ = 1 2 + δ A ℓ , δ A ℓ = ℓ ′ = 2 ∑ ℓ a ℓ ′ Δ A ℓ ′ , a ℓ = ℓ ′ = 1 ∏ ℓ n ℓ ′ , (138) B ^ ℓ ≡ B ℓ b ℓ = 1 N o u t + δ B ℓ , δ B ℓ = ∑ ℓ ′ = ℓ L − 1 Δ B ℓ ′ b ℓ ′ , b ℓ = ∏ ℓ ′ = ℓ + 1 L n ℓ ′ − 1 . (139) \widehat{B}_\ell\equiv\frac{B_\ell}{b_\ell}=\mathbb{1}_{N_{\mathrm{out}}}+\delta B_\ell,\qquad
\delta B_\ell=\sum_{\ell'=\ell}^{L-1}\frac{\Delta B_{\ell'}}{b_{\ell'}},\qquad
b_\ell=\prod_{\ell'=\ell+1}^{L}n_{\ell'-1}. \tag*{(139)} B ℓ ≡ b ℓ B ℓ = 1 N out + δ B ℓ , δ B ℓ = ℓ ′ = ℓ ∑ L − 1 b ℓ ′ Δ B ℓ ′ , b ℓ = ℓ ′ = ℓ + 1 ∏ L n ℓ ′ − 1 . (139) Each term in these sums has zero mean given the layers already summed, so the terms are uncorrelated and their variances add: δ A ℓ \delta A_\ell δ A ℓ ( ∑ ℓ ′ < ℓ 1 / n ℓ ′ ) 1 / 2 \left(\sum_{\ell'<\ell}1/n_{\ell'}\right)^{1/2} ( ∑ ℓ ′ < ℓ 1/ n ℓ ′ ) 1/2 δ B ℓ \delta B_\ell δ B ℓ ( N o u t ∑ ℓ ′ ≥ ℓ 1 / n ℓ ′ ) 1 / 2 \left(N_{\rm out}\sum_{\ell'\geq\ell}1/n_{\ell'}\right)^{1/2} ( N out ∑ ℓ ′ ≥ ℓ 1/ n ℓ ′ ) 1/2
Substituting into eq. (134) gives G n n = r 2 γ ∑ ℓ w ℓ A ^ ℓ ⊗ B ^ ℓ G_{nn}=r^{2\gamma}\sum_\ell w_\ell\widehat{A}_\ell\otimes\widehat{B}_\ell G nn = r 2 γ ∑ ℓ w ℓ A ℓ ⊗ B ℓ w ℓ = a ℓ b ℓ / σ ℓ 2 w_\ell=a_\ell b_\ell/\sigma_\ell^2 w ℓ = a ℓ b ℓ / σ ℓ 2
G 2 = ∑ ℓ w ℓ A ^ ℓ ∑ ℓ w ℓ , F = r 2 γ ∑ ℓ w ℓ B ^ ℓ . (140) G_2=\frac{\sum_\ell w_\ell\widehat{A}_\ell}{\sum_\ell w_\ell},\qquad\mathcal{F}=r^{2\gamma}\sum_\ell w_\ell\widehat{B}_\ell. \tag*{(140)} G 2 = ∑ ℓ w ℓ ∑ ℓ w ℓ A ℓ , F = r 2 γ ℓ ∑ w ℓ B ℓ . (140) When G 2 ⊗ F G_2\otimes\mathcal{F} G 2 ⊗ F δ A ℓ \delta A_\ell δ A ℓ δ B ℓ \delta B_\ell δ B ℓ G n n G_{nn} G nn
G n n = G 2 ⊗ F + R , R = r 2 γ ∑ ℓ w ℓ ( δ A ℓ − δ A ‾ ) ⊗ ( δ B ℓ − δ B ‾ ) , (141) G_{nn}=G_2\otimes\mathcal{F}+R,\qquad R=r^{2\gamma}\sum_\ell w_\ell(\delta A_\ell-\overline{\delta A})\otimes(\delta B_\ell-\overline{\delta B}) , \tag*{(141)} G nn = G 2 ⊗ F + R , R = r 2 γ ℓ ∑ w ℓ ( δ A ℓ − δ A ) ⊗ ( δ B ℓ − δ B ) , (141) where the bars denote averages with weights w ℓ w_\ell w ℓ ∥ X ⊗ Y ∥ = ∥ X ∥ ∥ Y ∥ \|X\otimes Y\|=\|X\|\|Y\| ∥ X ⊗ Y ∥ = ∥ X ∥∥ Y ∥ δ A ℓ \delta A_\ell δ A ℓ δ B ℓ \delta B_\ell δ B ℓ
∥ R ∥ ∥ G n n ∥ = O ( N o u t ∑ ℓ = 1 L − 1 1 n ℓ ) , (142) \frac{\|R\|}{\|G_{nn}\|}=\mathcal{O}\left(\sqrt{N_{\rm out}\sum_{\ell=1}^{L-1}\frac{1}{n_\ell}}\right), \tag*{(142)} ∥ G nn ∥ ∥ R ∥ = O N out ℓ = 1 ∑ L − 1 n ℓ 1 , (142) i.e. a few percent for the widths used here. In other words, all terms in the sum of eq. (134) are proportional to one another up to corrections that shrink as the hidden layers get wider, and the leading corrections are absorbed into G 2 G_2 G 2 F \mathcal{F} F (132) factorises accordingly, so that we can write
G n n S ^ ≃ ( G ^ 2 S ^ 2 ) ⊗ ( F diag ( m i ′ ) N o u t ) . (143) G_{nn}\widehat{S}\simeq(\widehat{G}_2\widehat{S}_2)\otimes(\mathcal{F}\operatorname{diag}(m_{i'})_{N_{\rm out}}). \tag*{(143)} G nn S ≃ ( G 2 S 2 ) ⊗ ( F diag ( m i ′ ) N out ) . (143) By the definition of the spectrum of the Kronecker product
spec ( A ⊗ B ) = { a i b j } j = 1 , … , m i = 1 , … , n , (144) \operatorname{spec}(A\otimes B)=\{a_i b_j\}_{\substack{j=1,\ldots,m\\i=1,\ldots,n}}, \tag*{(144)} spec ( A ⊗ B ) = { a i b j } j = 1 , … , m i = 1 , … , n , (144) where a i a_i a i b j b_j b j A A A B B B 2 N o u t 2N_{\rm out} 2 N out
λ ± , k ( n ) = μ ± ν k , μ ± ∈ spec ( G ^ 2 S ^ 2 ) , ν k ∈ spec ( F diag ( m i ′ ) N o u t ) , (145) \lambda_{\pm,k}^{(n)}=\mu_\pm\nu_k,\qquad\mu_\pm\in\operatorname{spec}(\widehat{G}_2\widehat{S}_2),\qquad\nu_k\in\operatorname{spec}(\mathcal{F}\operatorname{diag}(m_{i'})_{N_{\rm out}}), \tag*{(145)} λ ± , k ( n ) = μ ± ν k , μ ± ∈ spec ( G 2 S 2 ) , ν k ∈ spec ( F diag ( m i ′ ) N out ) , (145) which reproduce the exact eigenvalues to a few percent in the regime of large γ \gamma γ (123) holds. The non-constant sector therefore splits into two tiers of N o u t N_{\rm out} N out μ + / μ − \mu_+/\mu_- μ + / μ −
The eigenvectors are obtained by projecting back, as in the constant sector. Using the three column blocks of eq. (98) , z = Ψ ζ z=\Psi\zeta z = Ψ ζ
z i α = ∑ i ′ ( R ( F K ) ) i ′ i [ ζ c , i ′ + ζ 2 , i ′ x α + ζ 3 , i ′ log x α ] , (146) z_{i\alpha}=\sum_{i'}\left(R^{(\mathrm{FK})}\right)_{i'i}\left[\zeta_{c,i'}+\zeta_{2,i'}x_\alpha+\zeta_{3,i'}\log x_\alpha\right], \tag*{(146)} z i α = i ′ ∑ ( R ( FK ) ) i ′ i [ ζ c , i ′ + ζ 2 , i ′ x α + ζ 3 , i ′ log x α ] , (146) and ζ c \zeta_c ζ c (131) . The multiplicities drop out of that step: both S c c S_{cc} S cc S c n S_{cn} S c n diag ( m i ′ ) \operatorname{diag}(m_{i'}) diag ( m i ′ ) (107) , which therefore cancels,
S c c − 1 S c n = ( ⟨ x ⟩ , ⟨ log x ⟩ ) ⊗ 1 N o u t , ζ c , i ′ = − ⟨ x ⟩ ζ 2 , i ′ − ⟨ log x ⟩ ζ 3 , i ′ , (147) S_{cc}^{-1}S_{cn}=(\langle x\rangle,\langle\log x\rangle)\otimes\mathbf{1}_{N_{\rm out}},\qquad\zeta_{c,i'}=-\langle x\rangle\zeta_{2,i'}-\langle\log x\rangle\zeta_{3,i'}, \tag*{(147)} S cc − 1 S c n = (⟨ x ⟩ , ⟨ log x ⟩) ⊗ 1 N out , ζ c , i ′ = − ⟨ x ⟩ ζ 2 , i ′ − ⟨ log x ⟩ ζ 3 , i ′ , (147) so that the constant term in eq. (146) does nothing but centre the other two. By eq. (143) the eigenvector ζ n \zeta_n ζ n ζ n = η ± ⊗ χ k \zeta_n=\eta_\pm\otimes\chi_k ζ n = η ± ⊗ χ k η ± \eta_\pm η ± χ k \chi_k χ k x x x
z i α ( ± , k ) ∝ [ ∑ i ′ χ k , i ′ ( R ( F K ) ) i ′ i ] [ cos φ ± ( x α − ⟨ x ⟩ ) + sin φ ± ( log x α − ⟨ log x ⟩ ) ] , (148) z_{i\alpha}^{(\pm,k)}\propto\left[\sum_{i'}\chi_{k,i'}\left(R^{(\mathrm{FK})}\right)_{i'i}\right]\left[\cos\varphi_\pm(x_\alpha-\langle x\rangle)+\sin\varphi_\pm(\log x_\alpha-\langle\log x\rangle)\right], \tag*{(148)} z i α ( ± , k ) ∝ [ i ′ ∑ χ k , i ′ ( R ( FK ) ) i ′ i ] [ cos φ ± ( x α − ⟨ x ⟩) + sin φ ± ( log x α − ⟨ log x ⟩) ] , (148) where ( cos φ ± , sin φ ± ) ≡ η ± (\cos\varphi_\pm,\sin\varphi_\pm)\equiv\eta_\pm ( cos φ ± , sin φ ± ) ≡ η ± μ + / μ − \mu_+/\mu_- μ + / μ − γ \gamma γ Figure 31 , which will be discussed later. The figure displays the PDF modes q ( r ) q^{(r)} q ( r ) (27) (first row of each block is at initialisation), which are combinations of the NTK eigenvectors z ( r ) z^{(r)} z ( r ) λ ( r ) \sqrt{\lambda^{(r)}} λ ( r ) z ( r ) z^{(r)} z ( r )
Note that S ^ 2 \widehat{S}_2 S 2 x x x log x \log x log x x x x r = N o u t r=N_{\rm out} r = N out log x \log x log x (148) , and look like straight lines in the logarithmic scale used in Figure 31 .
The tier structure of the NTK spectrum
The derivation above has two consequences for the NTK at initialisation. First, the dimension of its image cannot exceed three times the number of independent flavours, dim Im ( Θ ) ≤ 3 N o u t = 24 \dim\operatorname{Im}(\Theta) \le3N_{\mathrm{out}} = 24 dim Im ( Θ ) ≤ 3 N out = 24 V l i n \mathcal{V}_{\mathrm{lin}} V lin N o u t N_{\mathrm{out}} N out γ \gamma γ
The first one is readily verified. The median dimension of the image of the NTK at initialisation, computed for the bare network and at the tolerance of Figure 8 , is
The predicted bound is approached monotonically as the initialisation width shrinks, and is exact at γ = 5 \gamma=5 γ = 5 dim Im ( Θ ) = 24 = 3 × 8 \dim\operatorname{Im}(\Theta)=24=3\times8 dim Im ( Θ ) = 24 = 3 × 8 γ \gamma γ (83) reports a negligible saturation, so that the two diagnostics agree on where the linear regime sets in. For the full model, preprocessing factors included, the dimensions reduce to those of eq. (70) : the preprocessing factors shrink the image throughout, consistently with what we already observed in Figure 8 for the baseline fit.
The spectrum of the NTK at initialisation is shown in Figure 29 against the eigenvalue rank r r r γ \gamma γ γ \gamma γ γ = 5 \gamma=5 γ = 5 γ = 4 \gamma=4 γ = 4 γ = 3 \gamma=3 γ = 3 r = 16 r=16 r = 16 γ = 2 \gamma=2 γ = 2
Insofar as the size of the eigenvalues is concerned, we see that the first 8 eigenvalues predicted above are reproduced in Figure 29 , that is eq. (128) . This precisely the 1, 0.75, 0.25 × \times × γ = 5 \gamma=5 γ = 5 Figure 29 , and corresponds to the subspace of constant functions in V l i n \mathcal{V}_{\mathrm{lin}} V lin x x x log x \log x log x (148) .
NTK dynamics at different initialisation widths
In this appendix we complement the discussion of § 5.2 and § 5.3 with a closer look at the NTK during the training of the six fits with different initialisation widths. Here we focus on the escape from the linear regime in Section B.1: we relate the growth of the image of the NTK to the saturation of the hidden layers, and show how the escape changes the eigenvectors of the NTK, in their functional form and in their orientation. Finally, we return to the analytical solution with the frozen NTK, at initialisation (Section B.2) and along the whole training (Section B.3).
Escape from the linear regime
We now look at the escape from the linear regime from three complementary angles: the dimension of the image of the NTK, which we relate to the saturation of the hidden layers, the functional form of the modes, and the orientation of the eigenvectors.
Image of the NTK and saturation. The dimension of the image of the NTK is shown in the left panel of Figure 30 as a function of the training epoch. When γ = 1 \gamma= 1 γ = 1 Figure 8 ). The case γ = 0.001 \gamma= 0.001 γ = 0.001 Figure 18 . The sudden increase signalling the escape from the linear regime is instead clearly visible starting from γ = 2 \gamma= 2 γ = 2 γ ≥ 3 \gamma\ge3 γ ≥ 3 τ e s c \tau_{\mathrm{esc}} τ esc γ = 5 \gamma= 5 γ = 5 dim Im ( Θ ) \dim\operatorname{Im}(\Theta) dim Im ( Θ ) τ e s c \tau_{\mathrm{esc}} τ esc τ e n d \tau_{\mathrm{end}} τ end
Such an increase can be traced back to the saturation of the internal layers, which sets in at the same epochs (right panel of Fig. 30). When this happens, the deep linear network approximation of Section A breaks down, and the higher-order terms of the Taylor expansion of the activations become relevant. This supports the picture of § 5.2 : with a small initialisation width the network is confined to the linear regime for the first phase of the training, and it is the saturation of the hidden layers that allows it to explore a richer functional space.
Functional form of the modes. In Figure 31 we show the modes of the training Hamiltonian, eq. (27) , for the singlet and the gluon, at initialisation (first row of each block), at τ e s c \tau_{\rm esc} τ esc τ e n d \tau_{\rm end} τ end N o u t N_{\rm out} N out
At initialisation and for γ = 4 \gamma= 4 γ = 4 5 5 5 V l i n \mathcal{V}_{\rm lin} V lin x x x r = N o u t r = N_{\rm out} r = N out log x \log x log x x ≳ 0.1 x \gtrsim0.1 x ≳ 0.1 x x x 3 N o u t 3N_{\rm out} 3 N out Figure 29 , seen now in the functional forms rather than in the eigenvalues. At γ = 3 \gamma= 3 γ = 3 x x x γ ≲ 2 \gamma\lesssim2 γ ≲ 2 γ = 1 \gamma= 1 γ = 1 γ = 0.001 \gamma= 0.001 γ = 0.001
Note that the simple functional forms at initialisation with a large γ \gamma γ γ \gamma γ Figure 18 : a network confined to the linear regime at the beginning of the training can eventually leave it, provided that the fit is given enough time.
Distance from the linear span. We quantify how well a mode q ( r ) q^{(r)} q ( r )
ρ ( r ) = ∥ q ( r ) − q l i n e a r ( r ) ∥ ∥ q ( r ) ∥ , (149) \rho^{(r)} = \frac{\left\lVert q^{(r)} - q_{\rm linear}^{(r)} \right\rVert}{\left\lVert q^{(r)} \right\rVert}, \tag*{(149)} ρ ( r ) = q ( r ) q ( r ) − q linear ( r ) , (149) where q l i n e a r ( r ) q_{\rm linear}^{(r)} q linear ( r ) q ( r ) q^{(r)} q ( r ) V l i n \mathcal{V}_{\rm lin} V lin 22 V l i n \mathcal{V}_{\rm lin} V lin Figure 32 displays it as a function of the rank at initialisation (left) and at the escape epoch τ e s c \tau_{\rm esc} τ esc γ \gamma γ N o u t N_{\rm out} N out Figure 29 . As γ \gamma γ γ ≲ 2 \gamma\lesssim2 γ ≲ 2 V l i n \mathcal{V}_{\rm lin} V lin γ \gamma γ V l i n \mathcal{V}_{\rm lin} V lin Figure 30 .
Rotation of the eigenvectors. The stationary behaviour of the spectrum on the initial plateau for large γ \gamma γ k k k
d k ( τ ) = 1 − 1 k ∑ i , j = 1 k ( z τ ( i ) , z 0 ( j ) ) 2 . (150) d_k(\tau) = 1 - \frac{1}{k}\sum_{i,j=1}^{k}\left(z_\tau^{(i)}, z_0^{(j)}\right)^2. \tag*{(150)} d k ( τ ) = 1 − k 1 i , j = 1 ∑ k ( z τ ( i ) , z 0 ( j ) ) 2 . (150) Here z τ ( i ) z_\tau^{(i)} z τ ( i ) z 0 ( j ) z_0^{(j)} z 0 ( j ) i i i j j j τ \tau τ ( ⋅ , ⋅ ) (\cdot,\cdot) ( ⋅ , ⋅ ) k k k d k d_k d k k k k 23
Figure 33 shows d k d_k d k γ \gamma γ γ \gamma γ γ = 5 \gamma= 5 γ = 5 τ ≃ 8 × 10 3 \tau\simeq8 \times10^3 τ ≃ 8 × 1 0 3 d 1 = 0.002 d_1 = 0.002 d 1 = 0.002 d 8 = 0.000 d_8 = 0.000 d 8 = 0.000 d 1 = 0.26 d_1 = 0.26 d 1 = 0.26 d 8 = 0.11 d_8 = 0.11 d 8 = 0.11 τ e s c \tau_{\rm esc} τ esc γ = 1 \gamma= 1 γ = 1 d 1 d_1 d 1
the first few hundred epochs and grows smoothly from the first epochs onwards, while γ = 3 \gamma= 3 γ = 3
The end-of-training values of d k d_k d k
{#eq:152}
Apart from the γ = 0.001 \gamma= 0.001 γ = 0.001 k k k d 1 = 0.90 d_1 = 0.90 d 1 = 0.90 γ \gamma γ d 24 d_{24} d 24 γ \gamma γ
Taken together, these results complete the picture of the linear regime drawn in Figure 18 and Figure 29 , thus proving that the NTK is stationary when confined to this regime and not just slowly varying. Its eigenvalues are pinned at the values fixed at initialisation, its eigenvectors do not move, and its image does not grow beyond the dimension of eq. (149) – exactly 3 N o u t 3N_{\rm out} 3 N out γ = 5 \gamma= 5 γ = 5 (128) is also exact.
Timescales of the frozen NTK at initialisation
In Section 5.3 we showed that, for large values of γ \gamma γ Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) τ ( i ) \tau^{(i)} τ ( i ) (66) for every mode in Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) τ N T K = 0 \tau_{\rm NTK} = 0 τ NTK = 0 Figure 34 against the training horizon τ e n d \tau_{\rm end} τ end γ \gamma γ
Figure 34. Characteristic resolution times τ ( i ) = 1 / ( 2 η h ( i ) ) \tau^{(i)} = 1/(2\eta h^{(i)}) τ ( i ) = 1/ ( 2 η h ( i ) ) (66) , in epochs, of every mode in Im ( Θ ) \operatorname{Im}(\Theta) Im ( Θ ) i i i γ \gamma γ Figure 10 , so that each curve terminates at the corresponding entry of eq. (70) . Points are medians over the N r e p = 120 N_{\rm rep} = 120 N rep = 120 τ e n d = 2 × 10 4 \tau_{\rm end} = 2 \times10^4 τ end = 2 × 1 0 4 N o u t = 8 N_{\rm out} = 8 N out = 8 τ e s c \tau_{\rm esc} τ esc
Within the training horizon, the frozen kernel of the γ = 4 \gamma= 4 γ = 4 γ = 5 \gamma= 5 γ = 5 γ \gamma γ
eq. (128) . Every direction carrying a dependence on x x x log x \log x log x τ ( i ) ≳ 3 × 10 6 \tau^{(i)} \gtrsim3 \times10^6 τ ( i ) ≳ 3 × 1 0 6 V l i n \mathcal{V}_{\mathrm{lin}} V lin
The frozen NTK along the training
Here we supplement the discussion of § 5.3 by assessing the frozen-NTK description at intermediate freezing epochs. Specifically, we consider τ N T K = 0.5 \tau_{\mathrm{NTK}} = 0.5 τ NTK = 0.5 0.75 0.75 0.75 1 1 1 1.5 1.5 1.5 2 2 2 τ e s c \tau_{\mathrm{esc}} τ esc γ \gamma γ τ N T K = 0 \tau_{\mathrm{NTK}} = 0 τ NTK = 0 τ N T K = τ e n d \tau_{\mathrm{NTK}} = \tau_{\mathrm{end}} τ NTK = τ end § 5.3 .24 τ N T K \tau_{\mathrm{NTK}} τ NTK f ^ \hat{f} f ^ f 0 f_0 f 0 § 5.3 , up to the end of the training, where the two solutions are compared. The distance between the analytical solution f ^ \hat{f} f ^
D ( τ N T K ) = ∥ ⟨ f ^ τ e n d ⟩ − ⟨ f τ e n d ⟩ ∥ M ∥ ⟨ f τ e n d ⟩ ∥ M (151) D(\tau_{\mathrm{NTK}}) = \frac{\left\lVert\left\langle\hat{f}_{\tau_{\mathrm{end}}} \right\rangle- \left\langle f_{\tau_{\mathrm{end}}} \right\rangle\right\rVert_M}{\left\lVert\left\langle f_{\tau_{\mathrm{end}}} \right\rangle\right\rVert_M} \tag*{(151)} D ( τ NTK ) = ∥ ⟨ f τ end ⟩ ∥ M ⟨ f ^ τ end ⟩ − ⟨ f τ end ⟩ M (151) where ⟨ ⋅ ⟩ \langle\cdot\rangle ⟨ ⋅ ⟩ ∥ v ∥ M = ( v T M v ) 1 / 2 \lVert v\rVert_M = (v^{T} M v)^{1/2} ∥ v ∥ M = ( v T M v ) 1/2 M M M (16) . Here M M M
Figure 35 shows D D D τ N T K \tau_{\mathrm{NTK}} τ NTK γ ≲ 1 \gamma\lesssim1 γ ≲ 1 τ N T K \tau_{\rm NTK} τ NTK γ = 1 \gamma= 1 γ = 1 D ≃ 0.08 D \simeq0.08 D ≃ 0.08 τ N T K = 0 \tau_{\rm NTK} = 0 τ NTK = 0 D ≃ 0.01 D \simeq0.01 D ≃ 0.01 τ N T K = τ e n d \tau_{\rm NTK} = \tau_{\rm end} τ NTK = τ end Figure 11 and Figure 12 , later kernels give better approximations, but the kernel at initialisation is already a reasonable one. With γ = 0.001 \gamma= 0.001 γ = 0.001 D D D § 5.2 and in agreement with [63 ] .
Figure 35. The discrepancy D D D (151) between the analytical and the trained solution at the end of the training, as a function of the epoch τ N T K \tau_{\rm NTK} τ NTK γ \gamma γ τ N T K \tau_{\rm NTK} τ NTK τ e s c \tau_{\rm esc} τ esc τ N T K = 0 \tau_{\rm NTK} = 0 τ NTK = 0 τ N T K \tau_{\rm NTK} τ NTK τ N T K = τ e s c \tau_{\rm NTK} = \tau_{\rm esc} τ NTK = τ esc D D D N r e p = 120 N_{\rm rep} = 120 N rep = 120
The picture is different from γ = 2 \gamma= 2 γ = 2 D ≃ 0.4 D \simeq0.4 D ≃ 0.4 0.5 0.5 0.5 γ ≲ 1 \gamma\lesssim1 γ ≲ 1 § 5.3 . As τ N T K \tau_{\mathrm{NTK}} τ NTK D D D τ N T K = 0.5 τ e s c \tau_{\mathrm{NTK}} = 0.5\,\tau_{\mathrm{esc}} τ NTK = 0.5 τ esc τ N T K = 1.5 τ e s c \tau_{\mathrm{NTK}} = 1.5\,\tau_{\mathrm{esc}} τ NTK = 1.5 τ esc 1.5 1.5 1.5 2 τ e s c 2\,\tau_{\mathrm{esc}} 2 τ esc D ≃ 0.03 D \simeq0.03 D ≃ 0.03 0.07 0.07 0.07 τ N T K \tau_{\mathrm{NTK}} τ NTK τ e s c \tau_{\mathrm{esc}} τ esc γ \gamma γ γ \gamma γ D D D γ \gamma γ τ e s c \tau_{\mathrm{esc}} τ esc τ N T K = τ e s c \tau_{\mathrm{NTK}} = \tau_{\mathrm{esc}} τ NTK = τ esc D ≃ 0.08 D \simeq0.08 D ≃ 0.08 0.12 0.12 0.12 0.18 0.18 0.18 0.15 0.15 0.15 γ = 2 \gamma= 2 γ = 2 3 3 3 4 4 4 5 5 5
Together, these results confirm the picture drawn by the subspace drift, eq. (152): the initialisation controls when the frozen-NTK description becomes accurate, but not how accurate it eventually is.
