Paper deep dive
Probability-Signature Dynamics: Unpacking Modular Addition Learning Within Two-Layer Networks
Yunji Wang, Junjie Yao, Linyu Liu, Pinyan Lu, Zhi-Qin John Xu
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 10/10/2026, 2:33:54 AM
Summary
This paper introduces a probability signature framework to explain the emergence of Fourier-structured representations in two-layer neural networks trained on modular addition tasks. By analyzing gradient interactions through conditional statistics of the training distribution, the authors demonstrate that these signatures correspond to cyclic shift operators diagonalized by the discrete Fourier transform. This leads to decoupled Fourier-mode dynamics, explaining phenomena such as Fourier sparsity, frequency matching, and phase alignment. The framework also resolves the puzzle of faster early loss decrease under label noise by linking it to increased conditional label collisions. The method is shown to be applicable to other operators like XOR.
Entities (8)
Relation Signals (6)
Probability Signature → appliesto → Modular Addition
confidence 98% · For modular addition, these signatures are cyclic shift operators and are diagonalized by the discrete Fourier transform
Probability Signature → explains → Fourier Circuit
confidence 95% · We utilize probability signatures to directly explain the emergence of dominant Fourier modes during the learning dynamics.
Fourier Circuit → supports → Exact Generalization
confidence 94% · Neural networks trained on modular addition tasks often develop Fourier-structured representations that support exact generalization.
Weight decay → determines → Fourier Structure Formation
confidence 92% · weight decay determines the formation of the Fourier structure as a threshold
Label Noise → causes → Faster Early Loss Decrease
confidence 90% · corrupted examples can show faster early loss decrease than clean examples... noise increases conditional label collisions, strengthening early shared-coordinate reinforcement.
Probability Signature → appliesto → XOR
confidence 85% · Taking XOR as an example, we observed the predicted frequency in experiments.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Neural networks trained on modular addition tasks often develop Fourier-structured representations that support exact generalization. While prior work has identified these Fourier circuits, the mechanism by which gradient-based training selects them from the data distribution remains unclear. We address this question using probability signatures, which express leading gradient interactions through conditional statistics of the training distribution. For modular addition, these signatures are cyclic shift operators and are diagonalized by the discrete Fourier transform, yielding approximately decoupled Fourier-mode dynamics. This explains the emergence of Fourier sparsity, frequency matching, and phase alignment. The same framework resolves a puzzle under label noise: corrupted examples can show faster early loss decrease than clean examples, despite lacking a coherent generalization rule. We show that noise increases conditional label collisions, strengthening early shared-coordinate reinforcement. Finally, this method can be applied to other operators. Taking XOR as an example, we observed the predicted frequency in experiments.
Tags
Links
- Source: https://arxiv.org/abs/2610.11833v1
- Canonical: https://arxiv.org/abs/2610.11833v1
Trouble viewing inline? Open PDF directly →
Full Text
86,666 characters extracted from source content.
Expand or collapse full text
Probability-Signature Dynamics: Unpacking Modular Addition Learning Within Two-Layer Networks Yunji Wang, Junjie Yao, Linyu Liu, Pinyan Lu, Zhi-Qin John Xu School of Mathematical Sciences, Shanghai Jiao Tong University Institute of Natural Sciences, MOE-LSC, Shanghai Jiao Tong University Huawei Technologies Ltd Shanghai University of Finance and Economics †These authors contributed equally to this work. *Corresponding author: xuzhiqin@sjtu.edu.cn Abstract Neural networks trained on modular addition tasks often develop Fourier-structured representations that support exact generalization. While prior work has identified these Fourier circuits, the mechanism by which gradient-based training selects them from the data distribution remains unclear. We address this question using probability signatures, which express leading gradient interactions through conditional statistics of the training distribution. For modular addition, these signatures are cyclic shift operators and are diagonalized by the discrete Fourier transform, yielding approximately decoupled Fourier-mode dynamics. This explains the emergence of Fourier sparsity, frequency matching, and phase alignment. The same framework resolves a puzzle under label noise: corrupted examples can show faster early loss decrease than clean examples, despite lacking a coherent generalization rule. We show that noise increases conditional label collisions, strengthening early shared-coordinate reinforcement. Finally, this method can be applied to other operators. Taking XOR as an example, we observed the predicted frequency in experiments. Keywords Two Layer Network, Modular Addition, Training Dynamics 1 Introduction The modular addition task is a canonical setting in which to study grokking. During the initial phase of training, an over-parameterized model can fit individual training examples yet fails to uncover an interpretable rule. After substantially prolonged training, however, its test performance may improve abruptly even though the training loss has already become very small. This delayed transition is commonly referred to as grokking. In modular addition, the onset of grokking is accompanied by the emergence of a distinctive structure in the weights. Prior studies have shown that trained networks develop Fourier circuits (Nanda et al., 2023; Gromov, 2023). Specifically, the input and output weights of many hidden units become concentrated at matching Fourier frequencies, and appropriate combinations of these Fourier components can compute (a+b)modP(a+b) P. Mechanistic interventions, including Fourier projection and ablation experiments, provide strong evidence that these components are functionally important for generalization. Their formation also reflects the transition from a memorizing solution to a generalizing one. These results answer an important reverse question: given a trained network, which circuit does it implement, and which components of that circuit support generalization? They have not fully answered the complementary forward question: why should gradient-based training produce Fourier-structured parameters from the modular-addition data distribution in the first place? The fact that Fourier features can represent modular addition does not imply that the training dynamics must discover them. Meanwhile, current work can already account for observed phenomena in neural networks, such as frequency competition and phase alignment (He et al., 2026). However, its theoretical analysis is based on a Fourier decomposition of the weights, which need to presuppose that the Fourier basis is related to modular addition. Our focus is therefore not to show again that Fourier-structured neurons can implement modular addition, nor to re-establish that these neurons support generalization. Instead, we study the missing forward direction training distribution⟹leading gradient terms⟹Fourier-mode selection.training distribution gradient terms -mode selection. More specifically, we ask the following questions. Which statistics of the training data enter the leading-order parameter dynamics? Why are Fourier modes approximately decoupled in the modular-addition problem? What determines whether a neuron is suppressed by weight decay or develops a dominant Fourier component? And, how to derive the weight structure in advance from the training distribution? We study these questions in a two-layer embedding MLP trained on modular addition dataset with some randomly corrupted labels. Our main analytical tool is the probability signature framework introduced by Yao and Xu (2025), which shows the correlation between data statistics and embedding without studying either modular addition learning or the key structure of data statistics. Our contributions are as follows. • We utilize probability signatures to directly explain the emergence of dominant Fourier modes during the learning dynamics. • We have explained the puzzle in the early memorization regime: uniformly random label noise consistently exhibits a faster initial loss decrease than clean data. Some statistics from training distribution can measure this difference. • The analysis framework can be extended to more tasks, such as XOR, enabling us to predict the structures that network parameters may form before training. 2 Related Work Grokking. Grokking, the delayed onset of generalization long after interpolation, was first reported on small algorithmic datasets, where weight decay proved essential (Power et al., 2022). Subsequent work explains the delay through properties of the optimization trajectory: effective theories of representation learning and norm growth (Liu et al., 2022), loss landscapes shaped by weight decay and initialization scale (Liu et al., 2023), instabilities of adaptive optimizers (Thilak et al., 2022), and the transition from lazy to rich training dynamics (Kumar et al., 2024). These characterize when and why generalization occurs, from the emergent Fourier circuits as a result of training. Based on this, we point out that weight decay determines the formation of the Fourier structure as a threshold, which establishes the connection between weight decay and Grokking. Modular arithmetic. Prior work has established that trained networks can implement modular addition through sparse Fourier circuits, and has introduced diagnostics such as frequency matching and phase alignment to characterize these circuits (Nanda et al., 2023). Analytic and scaling-oriented studies have further related the norms and spectral profiles of Fourier modes to the learned solution (Gromov, 2023; Michaud et al., 2023). These works provide important mechanistic and functional descriptions of Fourier structure in trained networks. Our analysis is based on these research results, but we present a forward perspective: starting from training distribution rather than from the Fourier circuits obtained through training. More recent work (He et al., 2026) provides a complete explanation for the frequency competition and phase alignment by decomposing neurons onto the Fourier basis. Compared with theirs, our contribution stands out in the fact that we do not require the prior knowledge of the Fourier circuit found in the modular-addition task. By studying the algebraic structure of the probability signature, we derive the Fourier basis as the common eigenvectors. We not only explain the above phenomena, but also we hope to provide a unified explanation for the puzzle of early noise and the extension to other training tasks. Learning with corrupted labels. Neural networks can both generalize from structured examples and memorize corrupted labels. This interaction has been studied broadly through early learning and memorization dynamics (Arpit et al., 2017; Liu et al., 2020). For algorithmic and modular-arithmetic tasks, prior work has documented robust generalization in the presence of label corruption and the coexistence of generalizing and memorizing subnetworks (Doshi et al., 2024; Liu and Lu, 2026). Our contribution is not a general theory of learning under label noise. Instead, we isolate one early-time phenomenon: randomly corrupted examples lead to faster initial loss decrease, and we relate it to probability signature, which is computable from the training data. 3 Preliminaries 3.1 Settings Let Ω=0,1,…,P−1 =\0,1,…,P-1\ be a finite alphabet, where P is a prescribed prime. We consider the prediction of a binary operation F:Ω×Ω→ΩF: × → from examples X=(a,b)X=(a,b) to the label F(a,b)F(a,b). All pairs in Ω×Ω × are divided into three disjoint subsets clean,noise,testD_clean,D_noise,D_test of cardinalities NcleanN_clean, NnoiseN_noise, and NtestN_test, respectively. The labels of clean and test are set to F(a,b)=(a+b)modPF(a,b)=(a+b) P, while the labels of noise are replaced by uniform random numbers. The training dataset consists of cleanD_clean and noiseD_noise. We define Ntrain=Nclean+NnoiseN_train=N_clean+N_noise, where Ntrain:Ntest=8:2N_train:N_test=8:2 and Nnoise=αNtrainN_noise=α N_train. α represents the noise ratio. Model. Let a,b∈ℝPe_a,e_b ^P be one-hot encodings. We study the two-layer MLP f(a,b)=Wσ(Ua+Vb)∈ℝPf(a,b)=W\,σ(Ue_a+Ve_b) ^P, where U,V∈ℝD×PU,V ^D× P, W∈ℝP×DW ^P× D, a,be_a,e_b means the one-hot encodings of input a,ba,b, and D≫PD P. The output distribution is softmax(f(a,b))softmax(f(a,b)). We use the cross-entropy loss function and employ L2L^2-weight decay with coefficient λ. Parameters are initialized independently according to Uij,Vij,Wij∼(0,D−2γ),U_ij,V_ij,W_ij (0,D^-2γ), where γ controls the initialization scale. Consider a continuous gradient flow model as an approximation, it is dθdt=−η(∇θL+λθ) dθdt=-η( _θL+λθ). For more details about initialization, learning rate or optimizer, please refer to the appendix A. 3.2 Probability signatures The Probability Signature proposed in the related work (Yao and Xu, 2025) is a statistical quantity derived from sampling and not dependent on training. It can characterize the main terms of the gradient at different training stages. In order to create the application environment for Probability Signature, we first extended the fully connected network with the same embedding to our model. Proposition 1 (Gradient flow). Assume the neural network f(x,y)=Wσ(Ux+Vy)f(x,y)=Wσ(Ue_x+Ve_y), using full batch GD with cross-entropy loss and ignoring weight decay, with the activation function σ(h)=0+c1h+c22h2+o(h2)σ(h)=0+c_1h+ c_22h^2+o(h^2) which is C3C^3 near zero. If we sample training data from (x,y,z)∼π(x,y,z) π, then we denote the probability that the first position is x by rxr_x and the random variable π[⋅|x]E_π[·|x] by the conditional expectation under the fixed first position. Under the assumption of small initialization that ∥U∥∞,∥V∥∞,∥W∥∞≤ϵ U _∞, V _∞, W _∞≤ε, we derived that dUxdt=rxη dU^xdt=r_xη \ c1W⊤(I−1PPP⊤)[π[z∣x]−c1PW(Ux+π[Vy∣x])] c_1W (I- 1P1_P1_P ) [E_π[e_z x]- c_1PW (U^x+E_π[V^y x] ) ] +c2π[(Ux+Vy)⊙W⊤(z−1P)|x]+δ. +c_2\,E_π [ (U^x+V^y ) W (e_z- 1P1_P)\, |\,x ] \+δ. where Ux∈ℝdU^x ^d represent the x-th column of U and ∥δ∥∞=O(ϵ3) δ _∞=O(ε^3) means the higher order terms. Proof. See Appendix E.1. From the gradient flow equation, we can further decompose the conditional expectation. Assume sampling from the training distribution (X,Y,Z)∼π(X,Y,Z) π, we define the probability signatures: φ(a)(x,z)=ℙ(Z=z∣X=x),φ(b)(x,z)=ℙ(Z=z∣Y=x), ^(a)(x,z)=P(Z=z X=x), ^(b)(x,z)=P(Z=z Y=x), the 3-order tensor joint probability signatures: Φ(a)(x,y,z)=ℙ(Z=z,Y=y∣X=x),Φ(b)(x,y,z)=ℙ(Z=z,X=y∣Y=x), ^(a)(x,y,z)=P(Z=z,Y=y X=x), ^(b)(x,y,z)=P(Z=z,X=y Y=x), and the trace-back probability signatures ψ(a)(x,y)=ℙ(Y=y∣X=x),ψ(b)(x,y)=ℙ(X=y∣Y=x).ψ^(a)(x,y)=P(Y=y X=x), ψ^(b)(x,y)=P(X=y Y=x). We also denote φx(a):=∑z∈Ωℙ(Z=z∣X=x)ez ^(a)_x:= _z∈ P(Z=z X=x)e_z, a column vector, and Φx(a):=∑y,z∈Ωℙ(Z=z,Y=y∣X=x)ey⊗ez ^(a)_x:= _y,z∈ P(Z=z,Y=y X=x)e_y e_z, a matrix. The same notation applies to the others. Then we can rewrite the gradient flow equation in Proposition 1 as dUxdt=rxη dU^xdt=r_xη \ c1W⊤[(φx(a)−1PP)−c1P(I−1PPP⊤)W(Ux+Vψx(a))] c_1W [( _x^(a)- 1P1_P)- c_1P(I- 1P1_P1_P )W(U^x+V _x^(a)) ] +c2[Ux⊙W⊤φx(a)−1P(Ux+Vψx(a))⊙(W⊤P)+diag(VΦ(a)xW)]+δ. +c_2 [U^x W _x^(a)- 1P(U^x+V _x^(a)) (W 1_P)+diag(V ^(a)_xW) ] \+δ. The above statistics φ,Φ,ψ , ,ψ are referred to as Probability Signatures. 4 Theoretical Analysis In this chapter, we hope to derive this frequency structure from a forward perspective. Reviewing the gradient flow equation 1, we divide it into φ -part, Φ -part and ψ-part based on the association with different probability signatures. Additionally, due to the weight decay, there is also a −λU-λ U term. The proof of the following theorems will be found in Appendix E.2,E.3, and E.4. Theorem 1 characterizes the dominant part of the gradient flow equation. Theorem 1 (Dominance of the gradient). Under the condition of 1 and non-zero weight decay, assume that η=1η=1 and φx(a)=φx(b)=ψx(a)=ψx(b)=1PP,∀x∈Ω, ^(a)_x= ^(b)_x=ψ^(a)_x=ψ^(b)_x= 1P1_P, ∀ x∈ , Under small initialization that ∥U∥∞,∥V∥∞,∥W∥∞≤ϵ U _∞, V _∞, W _∞≤ε we have dUxdt dU^xdt =−λUx+c[diag(VΦx(a)W)−V¯⊙W¯]+δU, =-λ U^x+c [diag(V ^(a)_xW)- V W ]+ _U, dVxdt dV^xdt =−λVx+c[diag(UΦx(b)W)−U¯⊙W¯]+δV, =-λ V^x+c [diag(U ^(b)_xW)- U W ]+ _V, dWzdt dW_zdt =−λWz+c[diag(UΨzV⊤)−U¯⊙V¯]+δW. =-λ W_z+c [diag(U _zV )- U V ]+ _W. where c=c2rx,∥δ∥∞≲ϵ3c=c_2r_x,\ δ _∞ ε^3. Here we define U¯:=P−1∑xUx,V¯:=P−1∑xVx,W¯:=P−1∑x(W⊤)x U:=P^-1 _xU^x,\ V:=P^-1 _xV^x,\ W:=P^-1 _x(W )^x, and Ψz:=∑x,y∈Ωℙπ[X=x,Y=y∣Z=z]ex⊗ey _z:= _x,y∈ P_π[X=x,Y=y Z=z]e_x e_y. 4.1 The Fourier basis derived from the signatures Let’s assume there is no noise and the training is conducted on the entire clean dataset. We define Φx+:=P−1∑j∈ℤPej⊗ej+x _x^+:=P^-1 _j _Pe_j e_j+x as the joint signature in this situation. Write Sx=PΦx+S_x=P _x^+, which is a shift-x permutation matrix. They satisfy SxSy=Sx+y,Sx∗=S−x,Sx∗Sx=I.S_xS_y=S_x+y, S_x^*=S_-x, S_x^*S_x=I. Here ∗ denotes the conjugate transpose and I is the P×P× P identity matrix. They form a commuting family of normal operators and therefore possess a common orthonormal eigenvectors over ℂC. The basis can be obtained directly from the simplest nontrivial shift: S1q=ωqS_1q=ω q implies F(y+1)=ωF(y)F(y+1)=ω F(y) for some unit root ωP=1ω^P=1, which yields Fr(s)=P−1/2e−2πirs/P,Φx+Fr=P−1e2πirx/PFr.F_r(s)=P^-1/2e^-2 /P, _x^+F_r=P^-1e^2 /PF_r. (1) These vectors are the normalized Fourier basis, and F:=[F0,F1,⋯,FP−1]F:=[F_0,F_1,·s,F_P-1] is the standard Fourier matrix. We use ^:=F u:=Fu to represent the Discrete Fourier Transform (DFT) of vector u. Estimate by the main term of Theorem 1, we find that this system is separable for each hidden unit d∈[D]d∈[D]. Thus we can consider this equation on one hidden state ,,∈ℝPu,v,w ^P, with system ˙x=c[diag(⊤Φx+)⏟νx−(P−1∑k=0P−1k)(P−1∑k=0P−1k)⏟νmean]−λx, u_x=c [ diag(v ^+_xw)_ _x- (P^-1 _k=0^P-1v_k)(P^-1 _k=0^P-1w_k)_ _mean ]-λu_x, We can expand the middle term νx _x by νx=⊤(F−1ΛxF)=(F∗)Λx(F)=(1,ωx,⋯,ωx(P−1))(F¯⊙F), _x=v (F^-1 ^xF)w=(F^*v) ^x(Fw)=(1,ω^x,·s,ω^x(P-1))( Fv Fw), At the same time, the Fourier transform of νmean _mean will be absorbed into the zero-frequency component, because it is independent of the index x. Multiplying on the left by F, we have: (F)r=C(F¯)r(F)r−λr,(r≠0),(Fu)_r=C( Fv)_r(Fw)_r-λu_r, (r≠ 0), Similarly, we find that for each frequency r≠0r≠ 0, ^˙r≈C^r^r¯−λ^r,^˙r≈C^r^r¯−λ^r,^˙r≈C^r^r−λ^r. u_r≈ C w_r v_r-λ u_r, v_r≈ C w_r u_r-λ v_r, w_r≈ C u_r v_r-λ w_r. where v means DFT and ¯ v means conjugate. That is, the system is decoupled at each nonzero frequency. Each 3-variable system ^r,^r,^r u_r, v_r, w_r corresponding to frequency r will evolve independently. The next theorem indicates that the conditions for neurons to grow depend on the initialization, and weight decay plays a filtering role. Moreover, growing neurons always develop a dominant frequency and adhere to specific phase patterns. Proposition 2 (Dominant frequency and phase alignment). Fix a hidden unit d. For each nonzero frequency r, let Ur(t)=U^d,r(t)U_r(t)= U_d,r(t), Vr(t)=V^d,r(t)V_r(t)= V_d,r(t), and Wr(t)=Wd^r(t)W_r(t)= W^d_r(t). Assume C>0C>0 and that the three mode amplitudes remain nonzero up to the first singularity. Let sr⋆s_r be the blow-up time of the corresponding undamped system. If there exists some r, sr⋆<C/λs_r <C/λ and sr⋆=mini≠0si,s_r = _i≠ 0\s_i\, then r is the first frequency to blow up. Hence it makes the dominant frequency. Moreover, it is going to hold argUr(t)+argVr(t)−argWr(t)→0 U_r(t)+ V_r(t)- W_r(t)→ 0, when t approaches the blow-up time. Through the proposition, we can summarize the mechanism of weight decay filtering. The argument sr⋆s _r depends on the initial weights of a group of neurons d,d,du_d,v_d,w_d, and weight decay provides a unified survival threshold, and the competition of patterns and the initial random asymmetry determine which frequencies cross this threshold. 4.2 Estimation of noise The actual training dataset is 80%80\% of the complete data with α noise inside. The following theorem will provide an approximation to bridge the gap caused by sampling. Theorem 2 (Estimation of noise). Under Proposition 1, sample D as above experiments with training data ratio ρ and uniform-label noise ratio α. Take U⋆x:=diag(VΦx+W)−V¯⊙W¯U^x_ :=diag(V ^+_xW)- V W as the main term under the full-addition signature. Then we estimate U˙x−ηc2(1−α)PU⋆x+λUx=ℰUx+δ, U^x- η c_2(1-α)PU^x_ +λ U^x=E_U^x+δ, where ∥δ∥∞≲ϵ3 δ _∞ ε^3. The sampling error is controlled by ℰ∙=0,Var(ℰ∙)≲[(1−ρ)+α]η2ρP3(c12ϵ2+c22ϵ4).E_DE_ =0, _D(E_ ) [(1-ρ)+α]η^2ρ P^3 (c_1^2ε^2+c_2^2ε^4 ). Similar results also hold for V and W. 5 Validation on Fourier frequency Figure 1: Accuracy-Epoch curves of λ=0,10−4,5×10−4,10−3λ=0,10^-4,5× 10^-4,10^-3. For λ=0λ=0, model reaches 100% on clean and noise, but 0 on test, overfitting. For λ=10−4λ=10^-4, model reaches almost 100% on clean, noise and test, generalization. For λ=5×10−4λ=5× 10^-4, model forgets the noise. And for λ=10−3λ=10^-3, model cannot converge. We regard all the column vectors of W or the row vectors of U and V as neurons, since the output of the neural network is f(a,b)=∑d=1DWd(Udea+Vdeb)f(a,b)= _d=1^DW^d(U_de_a+V_de_b). In analysis, we assume that they evolve approximately independently. We denote the concatenation of the d-th row of U,VU,V, the d-th column of W by θd=(Ud,Vd,Wd) _d=(U_d,V_d,W^d) and the collection of θd _d by V. We consider Fourier transform of some neuron u∈ℝPu ^P, it is u^k=∑x=0P−1uxe−i2πkx/P u_k= _x=0^P-1u_x\,e^-i2π kx/P for k=0,…,P−1.k=0,…,P-1. We define the dominant frequency and sparsity by τ(u):=argmaxk≠0|u^k|,ρ(u)=2|u^τ(u)|2∑k=1P−1|u^k|2∈[0,1].τ(u):=argmax_k≠ 0| u_k|, ρ(u)= 2| u_τ(u)|^2 _k=1^P-1| u_k|^2∈[0,1]. It is noted that frequency r and P−rP-r are regarded as the same under Fourier transform. Experiments. We conducted experiments under different weight decay settings, with the same σ(x)=GELU(x)σ(x)=GELU(x), noise ratio α=0.1α=0.1 and initialization scale γ=0.5γ=0.5, in Figure 1. Additional experiments with other activation functions are reported in Appendix B,D. It is clearly evident from Figure 1 that weight decay has a significant impact on the direction of the training process. We observe that λ=0.0005λ=0.0005 acts as a regularization. That is, the model generalizes well while ignoring all the noise. We naturally derived from weight decay to consider observing the norms of each group of neurons in the network. The growth of the norms can be described as: 12dt‖θd‖2=−η⟨θd,∇dL⟩⏟CE drive−λ‖θd‖2⏟decay. 12 ddt\| _d\|^2= -η _d, _dL _CE drive- λ\| _d\|^2_decay. Figure 2: (A) Scatter between Fourier sparsity and the norm of neurons at epoch 10, 500, 1000, 2000. (B) The highest and lowest sparsity 100 neurons are selected, denoted by tV_t and bV_b. We observe ℓdrivetop=∑θ∈tℓdrive¯ _drive^top= _θ _t _drive and similarly for ℓdecaytop,ℓdrivebot,ℓdecaybot _decay^top, _drive^bot, _decay^bot (C) The clean/noise dataset contributions to the drive. It is ℓdrivetop,clean=∑θ∈tℓdrive¯∣clean _drive^top,clean= _θ _t _drive _D_clean and similarly for others. In Figure 2A, we present the evolution of the norm and frequency sparsity of the neurons. It can be observed that due to weight decay, the norm of most neurons is continuously reduced to 0, while a few neurons can grow. Meanwhile, all the new-grown neurons have relatively high sparsity. In Figure 2B, we report that norm growth is concentrated on high-sparsity neurons, and the growth of the norm and the improvement of the accuracy occur simultaneously in terms of time, according to Figure 1. We also compare separately on clean and noisy examples in Figure 2C. It shows that CE–drive term of high-sparsity neurons is contributed primarily by clean data. And these driving forces are almost completely unrelated to the noise in the data. Cosine Similarity Check. Theorem 1 and 2 provide two estimates for the real gradient flow. We conducted numerical experiments in Appendix C.1, comparing the cosine similarity between the estimated gradient and the actual gradient. The results show that it can achieve cosine-similarity ≥0.9≥ 0.9 with the actual gradient. We find that the dominant frequencies of UdU_d, VdV_d, and WdW^d tend to agree if the neuron reaches high sparsity. If the common frequency is τ, we extract the components at the dominant frequency by (d,d,d)=((Ud^)τ,(Vd^)τ,(Wd^)τ),(u_d,v_d,w_d)=(( U_d)_τ,( V_d)_τ,( W^d)_τ), then their phases will tend to satisfy arg(d)+arg(d)=arg(d)(mod2π). (u_d)+ (v_d)= (w_d) 2π. Figure 3 shows that frequency matching and phase alignment emerge together during training, which is consistent with the result of Prop 2. Full-batch GD check. Using the full addition table, a GELU-activated network, and shared initialization, we compare exact CE training with the truncated model. We simulate both systems over 50 independently evolved neurons and report the dominant frequencies to which they ultimately converge, as well as the frequency sparsity. For the valid neurons, the dominant frequencies predicted by the truncated model exhibit 93% agreement with the converged dominant frequencies of the actual GELU-network. Experimental details are provided in the Appendix C.2. Figure 3: Frequency matching and phase alignment. In the odd columns, the horizontal coordinate is the dominant nonzero Fourier frequency of w and the vertical coordinate is its frequency sparsity; in the even columns, the coordinates are arg()+arg()mod2π (u)+ (v) 2π and arg()mod2π (w) 2π, respectively. Red points indicate exact agreement of the dominant frequencies of du_d, dv_d, and dw_d, while blue points indicate a mismatch; dashed lines show phase equality. 6 Application of Probability Signature 6.1 The puzzle of early noise Figure 4: Early(<100<100 epochs) clean and noise loss dynamics. α means the noise ratio. Noise loss always decreases more rapidly. We identify a puzzle in the early training: uniformly random label noise consistently exhibits a faster initial loss decrease than clean data. It was also recorded by (Liu and Lu, 2026) but has not been resolved. From Figure 4, this phenomenon holds under almost any ratio of noise. Although this may appear inconsistent with the frequency principle (Xu et al., 2020), which states that high-frequency data components are fitted later and low-frequency components earlier. In the discrete setting, the similar redefinition of frequency will be the probability that the predicted label will change while changing one of the input tokens. In this section, we employ probability signatures to quantify this speed difference. In the small-initialization regime, the loss derivative is approximately L˙∝−∑i=1Nσ(Uxi+Vyi)⊤(∑j=1Nzj=ziσ(Uxj+Vyj)), L - _i=1^Nσ(U^x_i+V^y_i) ( _j=1^N1_\z_j=z_i\σ(U^x_j+V^y_j) ), with sampling (xi,yi,zi)(x_i,y_i,z_i) from the training distribution. We ignore weight decay because we are comparing the loss of different parts within the same training. From this, we consider counting the numbers of sample pairs that share common coordinates. For pairs with distinct coordinates define κ(0)=[σ(α+β)⊤σ(γ+δ)]κ(0)=E[σ(α+β) σ(γ+δ)], while for pairs sharing one coordinate define κ(1)=[σ(α+β)⊤σ(α+γ)]κ(1)=E[σ(α+β) σ(α+γ)], where α,β,γ,δ∼(0,1/D)α,β,γ,δ (0,1/D) are i.i.d. vectors. Generally, we have κ(1)>κ(0)κ(1)>κ(0) (see appendix F.1). Let nzn_z be the number of examples labeled z and define K0=∑znz2K_0= _zn_z^2. Then the following theorem can compare the early loss for full training dataset. Note that φ(a),φ(b) ^(a), ^(b) here are understood as matrices. Theorem 3 (Early loss decrease). Suppose the training dataset contains every distinct pairs (a,b)∈Ω×Ω(a,b)∈ × , and the model structure is the same as that mentioned before. For different label distribution, the early-time loss-slope is positively correlated with L˙CE∝−κ(0)K0−(κ(1)−κ(0))K1,K1=P2(∥φ(a)∥22+∥φ(b)∥22−2).E L_CE -κ(0)K_0-(κ(1)-κ(0))K_1, K_1=P^2( ^(a) _2^2+ ^(b) _2^2-2). Proof. See Appendix E.5. For modular addition, conditioning on any fixed value of a (or b), the label is uniform on Ω . Hence, φ(a)=φ(b)=1PPP⊤,(K0,K1)=(P3,0), ^(a)= ^(b)= 1P1_P1_P , (K_0,K_1)=(P^3,0), For a fully random labeling, let Zx,y=ℙπ(Z=y∣X=x)Z_x,y=P_π(Z=y X=x). The P pairs with first coordinate x receive independent labels, so Zx,y∼P−1Binomial(P,P−1).Z_x,y P^-1Binomial(P,P^-1). According to the binomial distribution, we get [Zx,y2]=2P−2−P−3.E[Z_x,y^2]=2P^-2-P^-3. It further gives [K0]=P3+P2−P,[K1]=2(P2−P)E [K_0 ]=P^3+P^2-P,E [K_1 ]=2(P^2-P), which are both greater than those in modular-addition case. Thus, uniform random labels can produce a faster early loss decrease than modular-addition labels. Next we consider the mixture of clean and noise. Without loss of generality, we assume sampling in complete Ω×Ω × and dividing it into clean samples (modular addition) and noise (uniform random label). The loss for one example (xi,yi,zi)(x_i,y_i,z_i) can be regarded as L˙i≈−κ(0)nzi2−(κ(1)−κ(0))Δi;Δi:=∑j=1,j≠iNtrainzi=zj(xi=xj+yi=yj). L_i≈-κ(0)n_z_i^2-(κ(1)-κ(0)) _i;\ \ _i:= _j=1,\ j≠ i^N_train1_z_i=z_j(1_x_i=x_j+1_y_i=y_j). Here the coefficient Δi _i is regarded as Δclean _clean for a clean sample, or Δnoise _noise for a noise. Taking expectations, the counts satisfy Δclean¯=2α(1−P−1), _clean=2α(1-P^-1), and Δnoise¯=2(1−P−1). _noise=2(1-P^-1). Consequently, Δnoise¯−Δclean¯=2(1−α)(1−P−1)>0, _noise- _clean=2(1-α)(1-P^-1)>0, so the loss of noise decreases more rapidly. This prediction is consistent with the experimental observations and can also be applied to other forms of noise. In Appendix F.1 and F.2, we conducted additional experiments, thereby verifying that the theoretical predictions presented here are reliable. 6.2 Extension of the Method: XOR as an example Figure 5: Hadamard frequency structure in the XOR experiment. For sampled neurons at epochs 00, 5050, 200200, the horizontal axis is ‖θd‖2\| _d\|_2 and the vertical axis is the sparsity of wdw_d. Red points denote neurons whose three dominant frequencies agree. The right panel reports test accuracy and the number of high-sparsity, frequency-matched neurons; high-sparsity means ρ>0.7ρ>0.7. The Fourier structure observed in modular addition has previously been identified and functionally validated from trained networks. However, these results primarily provide a mechanistic description of the circuit after it has formed, rather than a general forward description to predict the representation basis associated with a new algebraic operation. Our probability-signature analysis suggests such a prescription. In the small-initialization regime, the leading parameter dynamic is governed by linear operators constructed from conditional statistics of the training distribution. When these probability-signature operators admit a common eigen basis, transforming the parameters into this basis approximately decouples the leading-order dynamics into independent modes. Thus, the relevant representation basis need not be assumed to be the Fourier basis in advance, instead it can be predicted by finding the common eigenvectors of the joint signature. This observation provides a natural route for transferring the analysis beyond modular addition. Starting from the Fourier transform, we notice that relevant probability signatures can be regarded as the operators on the Abelian group, whose common eigenvectors are given by the corresponding group characters whenever the operators are simultaneously diagonalizable. For example, consider the task on Ω=0,…,63 =\0,…,63\, with F(a,b)=axorbF(a,b)=axorb. Its joint signature is Φx(a)=∑y,z∈Ωℙ(Z=z,Y=y∣X=x)ey⊗ez=164∑y∈Ω(y,xxory). ^(a)_x= _y,z∈ P(Z=z,Y=y X=x)e_y e_z= 164 _y∈ E_(y,xxory). Hadamard basis simultaneously diagonalize these operators: there is an orthogonal basis H=[H0,H1,…,H63]H=[H_0,H_1,…,H_63] such that H−1Φx(a)H^-1 ^(a)_xH is diagonal for every x∈Ωx∈ . The properties of H can be referred to in Appendix E.6. Therefore, we claim that after training on XOR, the weights will evolve into a shape similar to the Hadamard basis. Validation. We write u=∑k=063Eu,kHku= _k=0^63E_u,kH_k, and define its sparsity by ρu=maxk|Eu,k|2∑k|Eu,k|2. _u= _k|E_u,k|^2 _k|E_u,k|^2. The existence of this orthogonal basis allows the above analysis to be applied to XOR. Based on the redefined frequency sparsity, in Figure 5 we present the observations regarding the frequency matching, as well as the high correlation between frequency sparsity and norm. The observed birth and growth of high-sparsity neurons supports our claim. 7 Conclusion and Discussion Using the probability-signature framework, we provide a unified data-side account of three phenomena. First, for modular addition, the conditional-distribution operators induced by the training data recover the Fourier basis as their common eigenbasis and yield approximately decoupled leading-order dynamics, thereby accounting for the emergence of dominant frequencies, frequency matching, and phase alignment. Second, under uniformly random label noise, conditional label collisions among examples sharing an input coordinate explain why corrupted examples can exhibit a faster initial loss decrease than clean examples. Third, the same operator viewpoint suggests a route for predicting the special structures that neurons may develop on datasets generated by other binary operations: for XOR, the joint signature identifies the Hadamard basis as the relevant candidate coordinate system, in agreement with the weight patterns observed after training. Limitations. We approximate the empirical gradient flow using idealized dynamics derived from discrete, full-distribution probability signatures, together with finite-sample and label-noise corrections. Moreover, our dynamical analysis concerns the dominant gradient terms within a particular training regime; we do not claim that this approximation remains valid throughout the entire training process. Further discussion of the scope and empirical accuracy of this approximation is provided in Appendix C.1. References Arpit et al. (2017) D. Arpit, S. Jastrzkebski, N. Ballas, D. Krueger, E. Bengio, M. S. Kanwal, T. Maharaj, A. Fischer, A. Courville, Y. Bengio, and S. Lacoste-Julien A closer look at memorization in deep networks. In International Conference on Machine Learning, Cited by: §2. Doshi et al. (2024) D. Doshi, A. Das, T. He, and A. Gromov To grok or not to grok: disentangling generalization and memorization on corrupted algorithmic datasets. In International Conference on Learning Representations, Cited by: §2. Gromov (2023) A. Gromov Grokking modular arithmetic. External Links: 2301.02679 Cited by: §1, §2. He et al. (2026) J. He, L. Wang, S. Chen, and Z. Yang On the mechanism and dynamics of modular addition: fourier features, lottery ticket, and grokking. External Links: 2602.16849 Cited by: §1, §2. Jacot et al. (2020) A. Jacot, F. Gabriel, and C. Hongler Neural tangent kernel: convergence and generalization in neural networks. External Links: 1806.07572, Link Cited by: §D.3. Kumar et al. (2024) T. Kumar, B. Bordelon, S. J. Gershman, and C. Pehlevan Grokking as the transition from lazy to rich training dynamics. In International Conference on Learning Representations, Cited by: §2. Liu and Lu (2026) L. Liu and P. Lu Unveiling memorization-generalization coexistence: a case study on arithmetic tasks with label noise. External Links: 2605.18022, Link Cited by: §2, §6.1. Liu et al. (2020) S. Liu, J. Niles-Weed, N. Razavian, and C. Fernandez-Granda Early-learning regularization prevents memorization of noisy labels. In Advances in Neural Information Processing Systems, Cited by: §2. Liu et al. (2022) Z. Liu, O. Kitouni, N. Nolte, E. J. Michaud, M. Tegmark, and M. Williams Towards understanding grokking: an effective theory of representation learning. External Links: 2205.10343 Cited by: §2. Liu et al. (2023) Z. Liu, E. J. Michaud, and M. Tegmark Omnigrok: grokking beyond algorithmic data. In International Conference on Learning Representations, Cited by: §2. Michaud et al. (2023) E. J. Michaud, Z. Liu, U. Girit, and M. Tegmark The quantization model of neural scaling. In Advances in Neural Information Processing Systems, Cited by: §2. Nanda et al. (2023) N. Nanda, L. Chan, T. Lieberum, J. Smith, and J. Steinhardt Progress measures for grokking via mechanistic interpretability. In International Conference on Learning Representations, Cited by: §1, §2. Power et al. (2022) A. Power, Y. Burda, H. Edwards, I. Babuschkin, and V. Misra Grokking: generalization beyond overfitting on small algorithmic datasets. External Links: 2201.02177 Cited by: §2. Thilak et al. (2022) V. Thilak, E. Littwin, S. Zhai, O. Saremi, R. Paiss, and J. Susskind The slingshot mechanism: an empirical study of adaptive optimizers and the grokking phenomenon. External Links: 2208.07339 Cited by: §2. Xu et al. (2020) Z. J. Xu, Y. Zhang, T. Luo, Y. Xiao, and Z. Ma Frequency principle: fourier analysis sheds light on deep neural networks. Communications in Computational Physics 28 (5), p. 1746–1767. Cited by: §6.1. Yao and Xu (2025) J. Yao and Z. J. Xu Probability signature: bridging data semantics and embedding structure in language models. External Links: 2509.20124, Link Cited by: §1, §3.2. Zhang et al. (2017) C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals Understanding deep learning requires rethinking generalization. External Links: 1611.03530, Link Cited by: §D.3. Appendix A Model and training details Figure 6: Accuracy curves using Adam, with weight decay λ=0,10−5,5⋅10−4,10−4,10−3λ=0,10^-5,5· 10^-4,10^-4,10^-3 and initialization scale γ=0.1,0.5,1.0γ=0.1,0.5,1.0. Figure 7: Loss curves using Adam. The neural network structure of this article can be written as f(a,b)=Wσ(Uea+Veb)=∑d=0D−1Wd(Ud,a+Vd,b),f(a,b)=Wσ(Ue_a+Ve_b)= _d=0^D-1W^d(U_d,a+V_d,b), Here, U and V are input embedding matrices in ℝD×PR^D× P, and W is the output matrix in ℝP×DR^P× D, with D representing the width of the hidden layer (D=2000D=2000, P=113P=113). The output probability can be expressed as p(Xi)zi=exp(f(a,b)zi)∑k=0P−1exp(f(a,b)k)p(X_i)_z_i= (f(a,b)_z_i) _k=0^P-1 (f(a,b)_k) The loss of the training data is given by the cross-entropy: LCE=1N∑i=1N−logp(Xi)ziL_CE= 1N _i=1^N- p(X_i)_z_i The total loss after adding weight decay is: L=LCE+λ2(‖U‖F2+‖V‖F2+‖W‖F2)L=L_CE+ λ2 (\|U\|_F^2+\|V\|_F^2+\|W\|_F^2 ) The learning rate strategy we employ is warmup + cosine decay: ηt=η⋅tTwarmup,t<Twarmupη⋅12(1+η0+(1−η0)cos(t−Twarmup)πToptim−Twarmup),t≥Twarmup _t= casesη· tT_warmup,&t<T_warmup\\ η· 12 (1+ _0+(1- _0) (t-T_warmup)πT_optim-T_warmup ),&t≥ T_warmup cases Figure 8: Accuracy curves using AdamW. Figure 9: Loss curves using AdamW. The minimum learning rate η0=10−5 _0=10^-5, and the initial learning rate η=10−3η=10^-3 if using the Adam/AdamW optimizer. Experiments using Adam can generally be reproduced under the AdamW optimizer with about 1000x weight decay. Figure 6,7 compares the specific accuracy or loss curves with Adam under different weight decay and initialization. If the weight decay is too small or absent, there will be no generalization. If the weight decay is too large, it may cause the training to fail to converge, and at this point, all parameters will be compressed to zero. With an appropriate weight decay, using a smaller initial value can shorten the training convergence time. Moreover, check Figure 8,9 for AdamW. Compared with Adam, greater weight decay needs to be used to achieve a similar effect, and we find that in this problem the Adam optimizer can provide better regularization, that is, a reduction of noise. Appendix B Discussion on Activation Functions Figure 10: (a)–(d): Final accuracy using activation functions σ(x)=xσ(x)=x (linear), σ(x)=x2σ(x)=x^2 (quadratic), σ(x)=tanh(x)σ(x)= (x), or σ(x)=max(0,x)σ(x)= (0,x), under different weight decay. We use σ(x)=GELU(x)σ(x)=GELU(x) as the activation function in the main text, which is smooth near zero. The analytical method in our text is based on the expansion of the activation function, hence we conducted the same experiments with other activation functions, such as σ(x)=xσ(x)=x (linear), σ(x)=x2σ(x)=x^2, σ(x)=tanh(x)σ(x)= (x), and σ(x)=ReLU(x)σ(x)=ReLU(x). The experimental results are presented in Figure 10. From the figure, it can be seen that using the linear activation is almost impossible to learn well on the training set. Using tanh(x) (x) as the activation function results in very low test generalization. However, using the quadratic or ReLU activation function can achieve great test generalization under an appropriate weight decay. Let’s review the mathematical derivation of the probability signature presented in the main text, we use the Taylor expansion of the activation function, σ(h)=0+σ′(0)h+12σ′(0)h2+o(h2),σ(h)=0+σ (0)h+ 12σ (0)h^2+o(h^2), to derive the probability signature. We can categorize all the expanded terms into 1-order or 2-order terms, based on whether they come from the 1st-order expansion or the 2nd-order expansion. In Theorem 1, the 1-order terms are eliminated because of the assumption of uniform marginal distribution. As a result, the second-order derivative term is the key to this problem. The theoretical explanation matches our observation of poor generalization with linear/Tanh. Because tanh(x)=x−x3/3+⋯ (x)=x-x^3/3+·s does not contain the second-order derivative term. Meanwhile, the quadratic activation function contains the second-order derivative term, which does well in the test case. Besides, the success of RELU is not in contradiction to the analysis, because it is not differentiable near 0. We will discuss the RELU activation function in detail in Appendix D. Appendix C Experimental Details C.1 Gradient Flow Approximation In this section, we present the experimental results with the aim of verifying the rationality of the approximation methods employed in Theorem 1 and Theorem 2. Figure 11: Cosine similarity between the approximated gradients and the actual gradients. The blue curve uses the main term from Theorem 1 and the red one uses the main term from Theorem 2. By Theorem 1, we can approximate the gradient flow of UxU^x as Ux˙≈−λUx+c⋅diag(VΦx(a)W) U^x≈-λ U^x+c·diag(V ^(a)_xW), where Φx(a) ^(a)_x is the probability signature of the actual training dataset. We extract its main term, ∇^UxL1:=diag(VΦx(a)W), ∇_U^xL_1:=diag(V ^(a)_xW), as an approximation for the cross-entropy gradient. Flatten this vector, we get ∇^θL1 ∇_θL_1, which is an approximation for the cross-entropy gradient ∇θL _θL. We compare the two gradients by cosine similarity cos⟨∇^θL1,∇θL⟩‖∇^θL1‖2⋅‖∇θL‖2, ∇_θL_1, _θL \| ∇_θL_1\|_2·\| _θL\|_2, as Approx 1 in Figure 11. Similarly, by Theorem 2, we can approximate the gradient flow of UxU^x as Ux˙≈−λUx+c⋅diag(VΦx+W) U^x≈-λ U^x+c·diag(V ^+_xW), where Φx+ ^+_x is the probability signature of the clean modular addition, it is, Φx+=1P(11⋱11)x ^+_x= 1P pmatrix&1&&&\\ &&1&&\\ &&& &\\ &&&&1\\ 1&&&& pmatrix^x We also extract and flatten its main term as ∇^θL2 ∇_θL_2, comparing which with the cross-entropy gradient ∇θL _θL by cosine similarity, as Approx 2 in Figure 11. The experimental results show that both the approximated gradients and the actual gradients have a high similarity of 90%. This verifies our theoretical analysis. Meanwhile, the control experiments using σ(x)=x2σ(x)=x^2 can serve as an ablation study, which demonstrates that using only the 2-order main terms is sufficient for achieving test generalization in the modular-addition task. Discussion on the small initialization assumption. In the main text, the approximation of the gradient is required to be expanded within a small parameter region. At the same time, we claim that the neurons we are studying can always grow. This note will clarify that this is not in contradiction. First of all, we acknowledge that there are higher-order terms in the gradient flow equation such as δ=O(h3)δ=O(h^3). Our theoretical derivation does not assume that these higher-order terms should be omitted. We focus on the dominant part. There exists a scale h0h_0. When the parameter is less than h0h_0, the second-order term becomes dominant, while when the parameter is greater than h0h_0, the higher-order terms become dominant. It is noted that when weight decay exists, there is a non-high-order term of −λθ-λθ in the gradient flow equation. Therefore, an appropriate λ will position the decay term between the above, and when the scale of the parameters expands beyond the region dominated by the second-order, the weight decay becomes the main term, thereby preventing the scale of the parameters from continuing to grow. C.2 Frequency Selectivity Figure 12: Neuron-wise frequency concentration in the truncated and original GELU systems. Heatmaps show the evolution of matched-frequency sparsity for individual neurons during training. Rows correspond to neurons, and columns represent training epochs. Numbers indicate the final winning Fourier mode, whereas ’x’ marks neurons whose winner had not stabilized by the specified threshold. Higher values indicate stronger concentration onto a matched nonzero frequency. We used complete addition tables with P=7,D=50P=7,D=50, GELU activation, Gaussian scale D−1/2D^-1/2. Both models started from identical initialized parameters. The exact update was θk+1=θk+0.2(−∇L(θk)−10−4θk), _k+1= _k+0.2(-∇ L( _k)-10^-4 _k), and the truncated update replaced −∇L-∇ L by the field in Theorem 1 with c=σ′(0)/7c=σ (0)/7 and the same weight decay λ=0.0001λ=0.0001. Recording every 25 epochs, the experiment tracked the evolution of the dominant frequency of each neuron, in Figure 12. We found that approximately 50% of the neurons developed a dominant frequency, and among this subset, a proportion of 0.93 converged to the same dominant frequency. Because shared initialization parameters were used, these neurons gave rise to the same frequency selection, which is consistent with our argument in Proposition 2. Appendix D Discussion on the ReLU activation In this section, we will discuss the ReLU activation function in detail. The ReLU activation function is not differentiable at 0, which makes it difficult to analyze using Taylor expansion. However, we can still analyze the behavior of the network with ReLU activation by considering the piecewise linear nature of ReLU. D.1 Fourier-mode is able to do modular addition under RELU We give a note on how 2L-MLP with ReLU activation do well in modular addition. Proposition 3 (A sufficient Fourier condition for modular addition). Let P be an odd prime, and let wd∈ℤP∖0, 0≤d<Dw_d _P \0\,\,0≤ d<D be prescribed frequencies. For each hidden unit d, let ϕ1,d,ϕ2,d _1,d, _2,d be mutually independent random variables distributed uniformly on [0,2π)[0,2π). If the parameters are set to Ud,x=cos(2πwdxP+ϕ1,d),Vd,x=cos(2πwdxP+ϕ2,d),Wc,d=cos(2πwdcP+ϕ1,d+ϕ2,d).U_d,x= \! ( 2π w_dxP+ _1,d ),\,V_d,x= \! ( 2π w_dxP+ _2,d ),\,W_c,d= \! ( 2π w_dcP+ _1,d+ _2,d ). then the prediction of the neural network is (a+b)modP(a+b) P with probability tending to 1 as D→∞D→∞. Proof. Fix a,b∈ℤPa,b _P, and write c⋆=(a+b)modP.c_ =(a+b) P. For each hidden unit d and output coordinate c, define Xd,c:=Wc,dmax0,Ud,a+Vd,b.X_d,c:=W_c,d \0,U_d,a+V_d,b\. Then fc(D)(a,b)=1D∑d=1DXd,c.f_c^(D)(a,b)= 1D _d=1^DX_d,c. Since |Wc,d|≤1|W_c,d|≤ 1 and 0≤max0,Ud,a+Vd,b≤2,0≤ \0,U_d,a+V_d,b\≤ 2, we have −2≤Xd,c≤2.-2≤ X_d,c≤ 2. We first compute the expectation of a single-unit contribution. Set αd=2πwdaP+ϕ1,d,βd=2πwdbP+ϕ2,d, _d= 2π w_daP+ _1,d, _d= 2π w_dbP+ _2,d, and let δd,c=2πwd(c−a−b)P. _d,c= 2π w_d(c-a-b)P. Because ϕ1,d+ϕ2,d=αd+βd−2πwd(a+b)P, _1,d+ _2,d= _d+ _d- 2π w_d(a+b)P, we obtain Wc,d=cos(αd+βd+δd,c).W_c,d= ( _d+ _d+ _d,c). Thus, [Xd,c]=[cos(αd+βd+δd,c)max0,cosαd+cosβd].E[X_d,c]=E [ ( _d+ _d+ _d,c) \0, _d+ _d\ ]. Since αd _d and βd _d are independent and uniformly distributed on [0,2π)[0,2π), we may set s=αd+βd2,t=αd−βd2.s= _d+ _d2, t= _d- _d2. Using cosαd+cosβd=2cosscost, _d+ _d=2 s t, we have [Xd,c]=1(2π)2∫02π∫02πcos(2s+δd,c)[2cosscost]+ts.E[X_d,c]= 1(2π)^2 _0^2π _0^2π (2s+ _d,c)\,[2 s t]_+\,dt\,ds. For fixed s, ∫02π[2cosscost]+t=4|coss|. _0^2π[2 s t]_+\,dt=4| s|. Consequently, [Xd,c] [X_d,c] =1π2∫02π|coss|cos(2s+δd,c)s = 1π^2 _0^2π| s| (2s+ _d,c)\,ds =cos(δd,c)π2∫02π|coss|cos(2s)s = ( _d,c)π^2 _0^2π| s| (2s)\,ds =43π2cos(δd,c). = 43π^2 ( _d,c). In the second equality, the sine term vanishes by symmetry, and in the last equality we used ∫02π|coss|cos(2s)s=43. _0^2π| s| (2s)\,ds= 43. Therefore, [Xd,c]=43π2cos(2πwd(c−a−b)P).E[X_d,c]= 43π^2 ( 2π w_d(c-a-b)P ). For the target coordinate c⋆c_ , we have c⋆−a−b≡0(modP),c_ -a-b≡ 0 P, and hence [Xd,c⋆]=43π2.E[X_d,c_ ]= 43π^2. Now consider any c≠c⋆c≠ c_ . Since P is prime, wd≢0(modP),c−a−b≢0(modP)w_d ≡ 0 P, c-a-b ≡ 0 P imply wd(c−a−b)≢0(modP).w_d(c-a-b) ≡ 0 P. Define ρP:=maxr∈ℤP∖0cos(2πrP). _P:= _r _P \0\ ( 2π rP ). Because ℤP∖0Z_P \0\ is finite and none of the corresponding angles is congruent to 00 modulo 2π2π, ρP<1. _P<1. It follows that, for every d and every c≠c⋆c≠ c_ , [Xd,c]≤43π2ρP.E[X_d,c]≤ 43π^2 _P. Hence, with ΔP:=43π2(1−ρP)>0, _P:= 43π^2(1- _P)>0, we have [Xd,c⋆]−[Xd,c]≥ΔPfor all d and all c≠c⋆.E[X_d,c_ ]-E[X_d,c]≥ _P all d and all c≠ c_ . Averaging over d gives [fc⋆(D)(a,b)]−[fc(D)(a,b)]≥ΔPfor every c≠c⋆.E[f_c_ ^(D)(a,b)]-E[f_c^(D)(a,b)]≥ _P every c≠ c_ . It remains to control the fluctuations around these expectations. Since the random phases are mutually independent across hidden units, the variables Xd,cd=1D\X_d,c\_d=1^D are independent for each fixed c. Moreover, Xd,c∈[−2,2]X_d,c∈[-2,2]. Therefore, Hoeffding’s inequality yields, for every ε>0 >0, ℙ(|fc(D)(a,b)−[fc(D)(a,b)]|≥ε)≤2exp(−Dε28).P ( |f_c^(D)(a,b)-E[f_c^(D)(a,b)] |≥ )≤ 2 (- D ^28 ). Taking ε=ΔP/3 = _P/3 and applying the union bound over the P output coordinates, we obtain ℙ(maxc∈ℤP|fc(D)(a,b)−[fc(D)(a,b)]|≥ΔP3)≤2Pexp(−DΔP272).P ( _c _P |f_c^(D)(a,b)-E[f_c^(D)(a,b)] |≥ _P3 )≤ 2P (- D _P^272 ). On the complementary event, for every c≠c⋆c≠ c_ , fc⋆(D)(a,b)−fc(D)(a,b) f_c_ ^(D)(a,b)-f_c^(D)(a,b) ≥([fc⋆(D)(a,b)]−[fc(D)(a,b)])−2ΔP3 ≥ (E[f_c_ ^(D)(a,b)]-E[f_c^(D)(a,b)] )- 2 _P3 ≥ΔP3>0. ≥ _P3>0. Thus, on this event, c⋆c_ is the unique maximizing output coordinate. Consequently, ℙ(argmaxc∈ℤPfc(D)(a,b)=c⋆)≥1−2Pexp(−DΔP272).P ( _c _Pf_c^(D)(a,b)=c_ )≥ 1-2P (- D _P^272 ). Since the right-hand side converges to 11 as D→∞D→∞, we conclude that ℙ(argmaxc∈ℤPfc(D)(a,b)=(a+b)modP)⟶1.P\! ( _c _Pf_c^(D)(a,b)=(a+b) P ) 1. ∎ D.2 Experiment Details of ReLU activation We conducted experiments using the ReLU activation function, and the results are shown in Figures 13 and 14. The experimental setup is consistent with that of the GELU activation function. The results indicate that the ReLU activation function exhibits similar behavior to GELU in terms of weight decay selectivity or frequency-phase scatter evolution. This suggests that the observed phenomena are not specific to those smooth activation function but rather reflect more general properties of neural network training dynamics. Figure 13: Experiment of Figure 2 with ReLU. Figure 14: Experiment of Figure 3 with ReLU. D.3 Early Memorization Figure 15: Accuracy curve under RELU activation and (λ,γ)=(0.001,0.1)(λ,γ)=(0.001,0.1), with 30% noise. Current research Zhang et al. (2017); Jacot et al. (2020) have already provided a thorough discussion on the early memory behavior of neural networks. When use σ(x)=max(0,x)σ(x)= (0,x), we observe clear traces of memory at (λ,γ)=(0.001,0.1)(λ,γ)=(0.001,0.1). The training curve from Figure 15 shows it reaches high accuracy on both clean and noise dataset first and then decreases due to the weight decay. Here, we merely use the gradient flow to offer a more intuitive understanding. Consider a single example (a,b,y)(a,b,y) and let h=Uea+Vebh=Ue_a+Ve_b and p=softmax(f(a,b))p=softmax(f(a,b)). By the chain rule, ∂L∂Wν,d=(p−ey)νσ(h)d+λWν,d,∂L∂Ud,x=x=a∧hd>0[W⊤(p−ey)]d+λUd,x. ∂ L∂ W_ν,d=(p-e_y)_ν\,σ(h)_d+λ W_ν,d, ∂ L∂ U_d,x=1_\x=a h_d>0\[W (p-e_y)]_d+λ U_d,x. Thus the gradient flow for the output weights is dWν,dt=(ey−p)νσ(h)d−λWν,d. dW_ν,ddt=(e_y-p)_νσ(h)_d-λ W_ν,d. Near initialization, the softmax expansion gives softmax(f)=1PP+1P(f−f¯)+O(∥f∥2/P),softmax(f)= 1P1_P+ 1P(f- f)+O( f ^2/P), so the correct output coordinate is amplified while the others are weakly suppressed. In particular, the active coordinates satisfy dWy⊤dt=σ(h)−λWy⊤, dW_y dt=σ(h)-λ W_y , hence the weights will increase in those dimensions where σ(h)>0σ(h)>0. For the example in Figure 16(a), σ(x)=max(0,x)σ(x)= (0,x), only d=3,5d=3,5 are active in the case, namely those satisfying hd=Ud,a+Vd,b>0h_d=U_d,a+V_d,b>0. This is only a toy model of 1 training sample. It behaves like an interpolation, storing the data onto these random dimensions. Figure 16: Example-wise interpolation and the memory score. (a). For the illustrative input 3+1=43+1=4 with (D,P)=(7,5),σ(x)=max(0,x)(D,P)=(7,5),σ(x)= (0,x), the activation feature K=e3+e5K=e_3+e_5 represents the initially active hidden coordinates; the displayed ±1± 1 values denote random initialization signs rather than fixed parameter values. (b). Evolution of the memory score on clean and noisy training examples. In the full-batch setting, the exact flows are dWc,dt=−1Ntrain∑i[pi,c−yi=c]max(0,Ud,ai+Vd,bi)−λWc,d, dW_c,ddt=- 1N_ train _i[p_i,c-1_\y_i=c\] (0,U_d,a_i+V_d,b_i)-λ W_c,d, dUd,adt=−1Ntrain∑i:ai=aUd,ai+Vd,bi>0∑r[pi,r−yi=r]Wr,d−λUd,a, dU_d,adt=- 1N_ train _i:a_i=a1_\U_d,a_i+V_d,b_i>0\ _r[p_i,r-1_\y_i=r\]W_r,d-λ U_d,a, dVd,bdt=−1Ntrain∑i:bi=bUd,ai+Vd,bi>0∑r[pi,r−yi=r]Wr,d−λVd,b. dV_d,bdt=- 1N_ train _i:b_i=b1_\U_d,a_i+V_d,b_i>0\ _r[p_i,r-1_\y_i=r\]W_r,d-λ V_d,b. These equations show that each example (X=(a,b),y)(X=(a,b),y) strengthens its own output coordinate (y)(y) on its active hidden dimensions and suppresses the other output coordinates (c≠y)(c≠ y), with the latter effect smaller by approximately a factor 1/P1/P. This example-wise interpolation motivates the following memory score: mem:=∑i=1N∑d=1DUd,ai+Vd,bi>0∑c≠yi(Wyi,d−Wc,d).mem := _i=1^N _d=1^D1_\U_d,a_i+V_d,b_i>0\ _c≠ y_i(W_y_i,d-W_c,d). From figure 16(b), we can observe that the memory score undergoes a rapid accumulation process in the early stages of training. Simultaneously, due to the weight decay’s effect of reducing the norm, it gradually decreases. Appendix E Proofs and Additional Details E.1 Calculation of the gradient flow Proof. Let hx,y:=Uex+Vey=Ux+Vy,px,y:=softmax(Wσ(hx,y)).h_x,y:=Ue_x+Ve_y=U^x+V^y, p_x,y:=softmax (Wσ(h_x,y) ). Without loss of generality, assume that the learning rate η=1η=1. For the softmax cross-entropy loss associated with one sample (x,y,z)(x,y,z), ℓ(x,y,z)=−log(px,y)z, (x,y,z)=- (p_x,y)_z, the derivative with respect to the logits f(x,y)=Wσ(hx,y)f(x,y)=Wσ(h_x,y) is ∇fℓ=px,y−ez. _f =p_x,y-e_z. Therefore, by the chain rule, ∇hx,yℓ=Diag(σ′(hx,y))W⊤(px,y−ez). _h_x,y =Diag (σ (h_x,y) )W (p_x,y-e_z). Since UxU^x contributes to hx,yh_x,y only when the first position of the sample is x, the full-batch gradient flow gives dUxdt=rxπ[σ′(hx,y)⊙W⊤(ez−px,y)|x]. dU^xdt=r_x\,E_π [σ (h_x,y) W (e_z-p_x,y)\, |\,x ]. Because σ is C3C^3 near zero and σ(h)=c1h+c22h2+o(h2),σ(h)=c_1h+ c_22h^2+o(h^2), where the powers are understood componentwise, we have σ′(h)=c1+c2h+o(‖h‖).σ (h)=c_11+c_2h+o(\|h\|). Moreover, σ(h)=c1h+O(‖h‖2).σ(h)=c_1h+O(\|h\|^2). Let ΠP:=I−1PPP⊤. _P:=I- 1P1_P1_P . The softmax function satisfies softmax(0)=1PP,softmax(0)= 1P1_P, and its Jacobian at zero is Dsoftmax(0) Dsoftmax(0) =Diag(1PP)−(1PP)(1PP)⊤ =Diag ( 1P1_P )- ( 1P1_P ) ( 1P1_P ) =1PI−1P2PP⊤ = 1PI- 1P^21_P1_P =1PΠP. = 1P _P. Consequently, for a logit vector g near zero, softmax(g)=1PP+1PΠPg+O(‖g‖2).softmax(g)= 1P1_P+ 1P _Pg+O(\|g\|^2). Taking g=Wσ(hx,y)g=Wσ(h_x,y) and using (3), we obtain g=c1Whx,y+O(‖hx,y‖2).g=c_1Wh_x,y+O(\|h_x,y\|^2). Substituting this expression into (4) gives px,y=1PP+c1PΠPWhx,y+O(‖hx,y‖2).p_x,y= 1P1_P+ c_1P _PWh_x,y+O(\|h_x,y\|^2). Thus, W⊤(ez−px,y)= W (e_z-p_x,y)= W⊤ez−1PW⊤P W e_z- 1PW 1_P −c1PW⊤ΠPWhx,y+O(‖hx,y‖2). - c_1PW _PWh_x,y+O(\|h_x,y\|^2). Combining (2) and (6), and retaining all constant and first-order terms in hx,yh_x,y, yields σ′(hx,y)⊙W⊤(ez−px,y) σ (h_x,y) W (e_z-p_x,y) = = (c1+c2hx,y+o(‖hx,y‖)) (c_11+c_2h_x,y+o(\|h_x,y\|) ) ⊙(W⊤ez−1PW⊤P−c1PW⊤ΠPWhx,y+O(‖hx,y‖2)) (W e_z- 1PW 1_P- c_1PW _PWh_x,y+O(\|h_x,y\|^2) ) = = c1W⊤ez−c1PW⊤P c_1W e_z- c_1PW 1_P +c2hx,y⊙(W⊤ez)−c2Phx,y⊙(W⊤P) +c_2h_x,y (W e_z)- c_2Ph_x,y (W 1_P) −c12PW⊤ΠPWhx,y+δx,y, - c_1^2PW _PWh_x,y+ _x,y, where δx,y _x,y contains all higher-order terms. In particular, the product −c1c2Phx,y⊙(W⊤ΠPWhx,y)- c_1c_2Ph_x,y (W _PWh_x,y ) is quadratic in hx,yh_x,y and is therefore included in δx,y _x,y. We now take the conditional expectation of each term in (7). First, π[c1W⊤ez−c1PW⊤P|x] _π [c_1W e_z- c_1PW 1_P\, |\,x ] =c1W⊤(π[ez∣x]−1PP). =c_1W (E_π[e_z x]- 1P1_P ). Since hx,y=Ux+Vyh_x,y=U^x+V^y, hx,y⊙(W⊤ez)= h_x,y (W e_z)= Ux⊙(W⊤ez)+Vy⊙(W⊤ez). U^x (W e_z)+V^y (W e_z). Hence, π[hx,y⊙(W⊤ez)|x] _π [h_x,y (W e_z)\, |\,x ] =Ux⊙(W⊤π[ez∣x])+π[Vy⊙(W⊤ez)∣x]. =U^x (W E_π[e_z x] )+E_π [V^y (W e_z) x ]. Furthermore, π[hx,y∣x]=Ux+π[Vey∣x].E_π[h_x,y x]=U^x+E_π[Ve_y x]. It follows from (10) that π[hx,y⊙(W⊤P)|x] _π [h_x,y (W 1_P)\, |\,x ] =(Ux+π[Vey∣x])⊙(W⊤P), = (U^x+E_π[Ve_y x] ) (W 1_P), and π[W⊤ΠPWhx,y|x] _π [W _PWh_x,y\, |\,x ] =W⊤ΠPW(Ux+π[Vey∣x]). =W _PW (U^x+E_π[Ve_y x] ). Substituting (8)–(12) into the gradient-flow equation (1), we obtain dUxdt=rx dU^xdt=r_x \ c1W⊤(π[ez∣x]−1PP) c_1W (E_π[e_z x]- 1P1_P ) +c2(Ux⊙(W⊤π[ez∣x])+π[Vy⊙(W⊤ez)∣x]) +c_2 (U^x (W E_π[e_z x] )+E_π [V^y (W e_z) x ] ) −c2P(Ux+π[Vey∣x])⊙(W⊤P) - c_2P (U^x+E_π[Ve_y x] ) (W 1_P) −c12PW⊤(I−1P⊤)W(Ux+π[Vey∣x])+δ, - c_1^2PW (I- 1P1_P1_P )W (U^x+E_π[Ve_y x] ) \+δ, where δ=O(ϵ3)δ=O(ε^3). The same calculation would be applied to V. Then we obtain the gradient flow for W. Denote rzr_z by the probability of samples labeled z. According to the definition, there is dWz⊤dt=(ez−px,y)σ(hx,y), dW_z dt=(e_z-p_x,y)\,σ(h_x,y), Then, the gradient flow of the z-th row of W is expressed in the expected form as dWz⊤dt= dW_z dt= c1(rzπ[hx,y∣Z=z]−1Pπ[hx,y]) c_1 (r_z\,E_π[h_x,y Z=z]- 1PE_π[h_x,y] ) +c22(rzπ[(hx,y⊙hx,y)∣Z=z]−1Pπ[hx,y⊙hx,y]) + c_22 (r_z\,E_π[(h_x,y h_x,y) Z=z]- 1PE_π[h_x,y h_x,y] ) −c12Pπ[hx,y⊤hx,y]W⊤ΠPez+δ. - c_1^2PE_π[h_x,y h_x,y]W _Pe_z+δ. Substituting hx,y=Ux+Vyh_x,y=U^x+V^y yields a similar result. ∎ E.2 Proof of Dominance of the gradient Proof. We give the proof for U and the same calculation would be able to applied to V and W. For convenience, write qx:=π[ez∣x],V¯x:=π[Vey∣x],q_x:=E_π[e_z x], V_x:=E_π[Ve_y x], and define ΠP:=IP−1PPP⊤. _P:=I_P- 1P1_P1_P . From proposition 1, we have dUxdt=rx dU^xdt=r_x \ c1W⊤(qx−1PP) c_1W (q_x- 1P1_P ) +c2Ux⊙(W⊤qx) +c_2U^x (W q_x) +c2π[Vy⊙(W⊤ez)∣x] +c_2E_π [V^y (W e_z) x ] −c2P(Ux+V¯x)⊙(W⊤P) - c_2P (U^x+ V_x ) (W 1_P) −c12PW⊤ΠPW(Ux+V¯x)+δ. - c_1^2PW _PW (U^x+ V_x ) \+δ. Since z|xz x is uniform on Ω , qx=π[ez∣x]=1PP.q_x=E_π[e_z x]= 1P1_P. Therefore, c1W⊤(qx−1PP)=0.c_1W (q_x- 1P1_P )=0. and c2Ux⊙W⊤(qx−1PP)=0.c_2U^x W (q_x- 1P1_P )=0. The joint product term is W⊤ΠPW(Ux+V¯x)=O(ϵ3).W _PW (U^x+ V_x )=O(ε^3). Substitute all these, we have dUxdt=rxc2[π[Vy⊙(W⊤ez)∣x]−V¯⋅W¯]+δ, dU^xdt=r_xc_2 [E_π [V^y (W e_z) x ]- V· W ]+δ, It remains to identify the conditional expectation with the stated matrix expression. By definition, diag(VΦx(a)W) (V _x^(a)W) =∑y,zℙ(Z=z,Y=y∣X=x)diag(Vy(ez⊤W)) = _y,zP(Z=z,Y=y X=x)\,diag (V^y(e_z W) ) =∑y,zℙ(Z=z,Y=y∣X=x)(Vy⊙W⊤ez) = _y,zP(Z=z,Y=y X=x)\, (V^y W e_z ) =π[Vy⊙(W⊤ez)∣x]. =E_π [V^y (W e_z) x ]. adding the weight decay term, we conclude that dUxdt=−λUx+rxc2[diag(VΦx(a)W)−V¯⋅W¯]+δ. dU^xdt=-λ U^x+r_xc_2 [diag (V _x^(a)W )- V· W ]+δ. ∎ E.3 Proof of Estimation of noise Proof. Let G_D be the empirical full-batch vector field. To avoid random conditional denominators, for every (x,y)∈Ω2(x,y)∈ ^2 let IxyI_xy be the indicator that the pair is included in the training set, and let Z~xy Z_xy be its observed label when Ixy=1I_xy=1. Uniform sampling and uniform label replacement give the exact identity [IxymZ~xy=z]=1P2((1−α)z=x+y+αP).E_D\! [ I_xym1_\ Z_xy=z\ ]= 1P^2 ((1-α)1_\z=x+y\+ αP ). The expectation in (1) is over the training pairs, the corrupted subset, and the replacement labels, with the parameters fixed. We prove the assertion for U first. Applying the expansion in Proposition 1 to the empirical measure yields GU,x+λUx=ηm∑y∈ΩIxy[ G^U_D,x+λ U^x= ηm _y∈ I_xy [ c1W⊤(eZ~xy−1PP) c_1W (e_ Z_xy- 1P1_P ) +c2(Ux+Vy)⊙W⊤(eZ~xy−1P)]+ΔU,x, +c_2(U^x+V^y) W (e_ Z_xy- 1P1_P ) ]+ ^U_D,x, For equation (1), summed over y, shows that [1m∑yIxy(eZ~xy−1PP)]=0.E_D\! [ 1m _yI_xy (e_ Z_xy- 1P1_P ) ]=0. Consequently, both the first-order term and the part of the second-order term containing UxU^x have zero sampling mean. For the remaining part, (1) gives [ηc2m∑yIxyVy⊙W⊤(eZ~xy−1PP)] _D\! [ η c_2m _yI_xyV^y W (e_ Z_xy- 1P1_P ) ] =ηc2(1−α)P[P−1∑yVy⊙W⊤ex+y−V¯⊙W¯]. = η c_2(1-α)P [P^-1 _yV^y W e_x+y- V W ]. This has already yielded the main components of the theorem. Exchanging U and V for another proof. For the x-th output row, set Ψx+:=∑jejex−j⊤ _x^+:= _je_je_x-j . The same Taylor expansion gives G,xW+λWx=ηm∑a,bIab(Z~ab=x−1P)[c1hab+c22(hab⊙hab)]+Δ,xW,G^W_D,x+λ W_x= ηm _a,bI_ab (1_\ Z_ab=x\- 1P ) [c_1h_ab+ c_22(h_ab h_ab) ]+ ^W_D,x, where hab=Ua+Vbh_ab=U^a+V^b. The linear term in (5) has zero expectation: on every addition diagonal a+b=xa+b=x, both ∑aUa _aU^a and ∑bVb _bV^b agree with their uniform-background sums. The same cancellation removes the Ua⊙UaU^a U^a and Vb⊙VbV^b V^b terms. The two equal cross terms leave ηc2(1−α)P[P−1∑aUa⊙Vx−a−U¯⊙V¯]. η c_2(1-α)P [P^-1 _aU^a V^x-a- U V ]. It remains to control the finite-sample fluctuation. Define ℰ∙E_ as the centered fluctuation of the displayed first- and second-order terms, and include the high order remainder. For simple random sampling without replacement, the finite-population variance formula gives a factor (1−ρ)/(ρP3)(1-ρ)/(ρ P^3) for a fixed input or output coordinate. Randomly corrupting an α fraction of the selected labels contributes at most an additional α/(ρP3)α/(ρ P^3). In (2) and (5), each coordinate of a centered first-order summand is bounded by 2|c1|ϵ2|c_1|ε, while each second-order summand is bounded by C|c2|ϵ2C|c_2|ε^2. The variance formula gives Var(ℰ∙)≤Cη2ρP3[(1−ρ)+α](c12ϵ2+c22ϵ4).Var_D(E_ )≤ Cη^2ρ P^3 [(1-ρ)+α ] (c_1^2ε^2+c_2^2ε^4 ). which is the claimed bound. ∎ It is noted that the error is caused by the ρ and α during the sampling. If the training set contains too few clean samples, the error between actual gradient and the predicted term will be very large, which is consistent with our intuition. E.4 Proof of dominant frequency and phase locking Proof. The Fourier modes obey U˙r=CVr¯Wr−λUr,V˙r=CUr¯Wr−λVr,W˙r=CUrVr−λWr. U_r=C V_rW_r-λ U_r, V_r=C U_rW_r-λ V_r, W_r=CU_rV_r-λ W_r. Set Ur=e−λtXr(s)U_r=e^-λ tX_r(s), Vr=e−λtYr(s)V_r=e^-λ tY_r(s), and Wr=e−λtZr(s)W_r=e^-λ tZ_r(s), with s=Cλ(1−e−λt).s= Cλ(1-e^-λ t). Then dXrds=Yr¯Zr,dYrds=Xr¯Zr,dZrds=XrYr. dX_rds= Y_rZ_r, dY_rds= X_rZ_r, dZ_rds=X_rY_r. If sr⋆<C/λs_r <C/λ is the undamped blow-up time, its damped time is tr⋆=−1λlog(1−λsr⋆C).t_r =- 1λ (1- λ s_r C ). This map is increasing, so sw⋆<si⋆s_w <s_i implies tw⋆<ti⋆t_w <t_i . As t↑tw⋆t t_w , the energy of mode w diverges while every competing mode remains finite, and its energy fraction tends to one. For phase locking, write U=AeiαU=Ae^iα, V=BeiβV=Be^iβ, W=GeiγW=Ge^iγ, and δ=α+β−γδ=α+β-γ. Direct separation of real and imaginary parts gives δ˙=−Csinδ(BGA+AGB+ABG). δ=-C δ ( BGA+ AGB+ ABG ). While the amplitudes remain nonzero, the factor in parentheses is positive. Thus δ=0δ=0 is stable and δ=πδ=π is unstable. A continuous random initialization avoids the unstable equilibrium with probability one, yielding δ(t)→0δ(t)→ 0 modulo 2π2π on the regular continuation interval. ∎ E.5 Proof of the early-loss theorem Lemma 1. Let a,b,c,da,b,c,d be i.i.d. random variables and let f:ℝ→ℝf:R be such that the expectations below exist, then [f(a+b)f(a+c)]≥[f(a+b)f(c+d)].E[f(a+b)f(a+c)] [f(a+b)f(c+d)]. Proof. Define h(a)=[f(a+b)∣a],μ=[f(a+b)].h(a)=E[f(a+b) a], μ=E[f(a+b)]. Since a+ba+b and c+dc+d are independent, [f(a+b)f(c+d)]=[f(a+b)][f(c+d)]=μ2.E[f(a+b)f(c+d)]=E[f(a+b)]\,E[f(c+d)]=μ^2. Given a, the variables b and c are independent, so [f(a+b)f(a+c)∣a]=[f(a+b)∣a][f(a+c)∣a]=h(a)2.E[f(a+b)f(a+c) a]=E[f(a+b) a]\,E[f(a+c) a]=h(a)^2. Hence [f(a+b)f(a+c)]=[h(a)2].E[f(a+b)f(a+c)]=E[h(a)^2]. Therefore, [f(a+b)f(a+c)]−[f(a+b)f(c+d)]=[h(a)2]−μ2.E[f(a+b)f(a+c)]-E[f(a+b)f(c+d)]=E[h(a)^2]-μ^2. But μ=[h(a)]μ=E[h(a)], so [h(a)2]−μ2=Var(h(a))≥0.E[h(a)^2]-μ^2=Var(h(a))≥ 0. Thus [f(a+b)f(a+c)]≥[f(a+b)f(c+d)]E[f(a+b)f(a+c)] [f(a+b)f(c+d)]. ∎ For activation functions such as GELU and RELU, the above inequality will hold strictly. It is κ(1)>κ(0)κ(1)>κ(0) in the main text. Theorem 3. Proof. The total number of all pairs with the same label is K0=∑znz2,K_0= _zn_z^2, Let nx,y(a)n^(a)_x,y count examples with first coordinate X=xX=x and label Z=yZ=y, and define nz,y(b)n^(b)_z,y analogously. In the complete table, nx,y(a)=Pφx,y(a)n^(a)_x,y=P ^(a)_x,y and nz,y(b)=Pφz,y(b)n^(b)_z,y=P ^(b)_z,y. The ordered pairs sharing the first coordinate and label number Ca=∑x,ynx,y(a)(nx,y(a)−1)=P2‖φ(a)‖22−P2,C_a= _x,yn^(a)_x,y(n^(a)_x,y-1)=P^2\| ^(a)\|_2^2-P^2, and similarly Cb=P2‖φ(b)‖22−P2C_b=P^2\| ^(b)\|_2^2-P^2. We define κ(l)=[σ(hi)σ(hj)]⋅D,l=0,1,2,κ(l)=E[σ(h_i)σ(h_j)]· D, l=0,1,2, Here (hi,hj)(h_i,h_j) are scalar jointly Gaussian, each of variance 2s22s^2 and covariance ls2ls^2. s denotes the initialization scale. The definition is equivalent to that in the main text. Thus, We can count it by using the principle of inclusion-exclusion as: −L˙/η=1N2(κ(0)K0+(κ(1)−κ(0))(Ca+Cb)+(κ(2)−κ(0))N),-E L/η= 1N^2 (κ(0)K_0+(κ(1)-κ(0))(C_a+C_b)+(κ(2)-κ(0))N ), The last item is independent of the training distribution. ∎ E.6 Propositions of Hadamard basis Let G=(ℤ/2ℤ)6G=(Z/2Z)^6 and identify the group operation on G with bitwise XOR. For each a∈Ga∈ G, define the matrix Φa∈ℝG×G _a ^G× G by Φa(y,z)=164z=a+y,y,z∈G. _a(y,z)= 1641_\z=a+y\, y,z∈ G. For u∈Gu∈ G, define the Walsh character χu(v)=(−1)u⋅v, _u(v)=(-1)^u· v, where u⋅v=∑j=16ujvj(mod2).u· v= _j=1^6u_jv_j 2. Let Hu∈ℝGH_u ^G be the normalized Walsh vector Hu(v)=18χu(v).H_u(v)= 18 _u(v). Then the family Φa:a∈G\ _a:a∈ G\ is simultaneously diagonalized by the orthonormal Walsh–Hadamard basis Hu:u∈G\H_u:u∈ G\. More precisely, ΦaHu=χu(a)64Hu=(−1)u⋅a64Hu. _aH_u= _u(a)64H_u= (-1)^u· a64H_u. Consequently, if H is the orthogonal matrix whose columns are the vectors HuH_u, then H⊤ΦaH=diag(χu(a)64)u∈GH _aH=diag ( _u(a)64 )_u∈ G for every a∈Ga∈ G. Indeed, for u,u′∈Gu,u ∈ G, ⟨Hu,Hu′⟩ H_u,H_u =164∑v∈Gχu(v)χu′(v) = 164 _v∈ G _u(v) _u (v) =164∑v∈G(−1)(u+u′)⋅v. = 164 _v∈ G(-1)^(u+u )· v. If u=u′u=u , then the summand is identically equal to 11, and hence ⟨Hu,Hu⟩=1. H_u,H_u =1. If u≠u′u≠ u , then u+u′≠0u+u ≠ 0. Choose an element v0∈Gv_0∈ G such that (u+u′)⋅v0=1.(u+u )· v_0=1. The map v↦v+v0v v+v_0 is a bijection of G, and therefore ∑v∈G(−1)(u+u′)⋅v _v∈ G(-1)^(u+u )· v =∑v∈G(−1)(u+u′)⋅(v+v0) = _v∈ G(-1)^(u+u )·(v+v_0) =−∑v∈G(−1)(u+u′)⋅v. =- _v∈ G(-1)^(u+u )· v. Thus the sum is zero. Hence ⟨Hu,Hu′⟩=u=u′, H_u,H_u =1_\u=u \, so Hu:u∈G\H_u:u∈ G\ is an orthonormal basis of ℝGR^G. Now fix a,u∈Ga,u∈ G. For every y∈Gy∈ G, we have (ΦaHu)(y) ( _aH_u)(y) =∑z∈GΦa(y,z)Hu(z) = _z∈ G _a(y,z)H_u(z) =164Hu(a+y). = 164H_u(a+y). Using the character identity χu(a+y)=χu(a)χu(y), _u(a+y)= _u(a) _u(y), we obtain (ΦaHu)(y) ( _aH_u)(y) =164⋅18χu(a+y) = 164· 18 _u(a+y) =164⋅18χu(a)χu(y) = 164· 18 _u(a) _u(y) =χu(a)64Hu(y). = _u(a)64H_u(y). Since this holds for every y∈Gy∈ G, ΦaHu=χu(a)64Hu. _aH_u= _u(a)64H_u. Thus every Walsh vector HuH_u is a common eigenvector of every matrix Φa _a, with eigenvalue λu(a)=χu(a)64=(−1)u⋅a64. _u(a)= _u(a)64= (-1)^u· a64. Let H denote the orthogonal matrix whose u-th column is HuH_u. Since the columns of H form a common orthonormal eigenbasis, we have H⊤ΦaH=diag(χu(a)64)u∈G.H _aH=diag ( _u(a)64 )_u∈ G. Therefore, the matrices Φa:a∈G\ _a:a∈ G\ are simultaneously orthogonally diagonalized by the normalized Walsh–Hadamard matrix. Appendix F Experiment of Various Label Noise F.1 Experiment on Modular Multiplication Task Figure 17: Early-training loss curves of the modular multiplication by 59,60,6159,60,61. In order to verify our theory, we designed the following experiment on the modular multiplication task. Consider N=61,Ω=1,2,…,NN=61, =\1,2,…,N\ and there are three modular-period that M1=59,M2=60,M3=61M_1=59,M_2=60,M_3=61. The total pairs of Ω×Ω × are divided into three parts S1,S2,S3S_1,S_2,S_3. The labels of data from Si=(xi,yi),i=1,2,3,S_i=\(x_i,y_i)\, i=1,2,3, are set to zi~=xi⋅yimodMi,zi~∈Ω∪0, z_i=x_i· y_i M_i, z_i∈ ∪\0\, We trained it using the same neural network as before. See Figure 17 for the early-training loss curves of the three tasks. We can see that the task with M2=60M_2=60 showed a significantly faster decline than the others. Our theory can easily explain this point. For M=60M=60, there exist a large number of pairs (x,y)(x,y) such that x⋅y≡x′⋅y(modM)x· y≡ x · y M, such as 5×15≡13×15(mod60),5× 15≡ 13× 15 60, which means that the label z is more concentrated. This leads to larger K0=∑znz2K_0= _zn_z^2 and larger ‖φ(a)‖22+‖φ(b)‖22\| ^(a)\|_2^2+\| ^(b)\|_2^2 , thus a faster early training loss decline according to Theorem 3. In contrast, for prime numbers M1=59M_1=59 and M3=61M_3=61, it satisfies x⋅y≡x′⋅y(modM)x· y≡ x · y M if and only if x≡x′(modM)x≡ x M for non-zero y, which means the signature of M=59,61M=59,61 are more uniformly distributed. This leads to a smaller K0K_0 and K1K_1 (we think that the values of M1=59M_1=59 and M3=61M_3=61 are very close to each other). As shown in Figure 17, the theoretical prediction is consistent with the actual observation. The early loss slope of M2=60M_2=60 is much greater than the others. F.2 Experiment on the Bijective Noise In the main text, we point out that in the scenario of uniform noise on the labels, during the early training phase, it is always observed that the loss of the noise decreases first, as in figure 4. Figure 18: Early-training loss curves of the bijective noise experiment. Here we present a specific noise construction that ensures it does not undergo such collisions. Just as our theory predicted, the phenomenon observed after replacing this kind of noise was completely different from what was observed before, as shown in Figure 18. Specifically, we construct a binary function g(x,y)g(x,y) such that: • g(x,⋅)g(x,·) is bijective for any fixed x and g(⋅,y)g(·,y) is bijective for any fixed y. • g is not isomorphic to addition, that is, g(x,y)≠σ3(σ1(x)+σ2(y)),g(x,y)≠ _3( _1(x)+ _2(y)), where σi _i is from the permutation group of Ω . • Generally, g cannot be accurately fitted by our 2L-MLP. Our final construction is g(x,y)=x,x=y,(x+(4+(y−x113))⋅(y−x))mod113,x≠y,g(x,y)= casesx,&x=y,\\ (x+(4+ ( y-x113 ))·(y-x)) 113,&x≠ y, cases where (⋅113) ( ·113 ) is the Quadratic Residue Legendre symbol (of modular 113). It satisfies the above three properties. Next, we replace the original uniform random noise with the function g(x,y)g(x,y). In this case, the loss of g does not always decrease first during the early training phase. By calculating the coefficient preceding κ(1)κ(1), We can give Δclean−Δg=2α−1, _clean- _g=2α-1, Inferring like that in chapter 6.1, a larger Δ gives a more negative early loss derivative. This indicates that when the noise ratio α<0.5α<0.5, the noise data (g(x,y)g(x,y)) decrease first. In contrast, when α>0.5α>0.5, the opposite occurs. As shown in Figure 18, the theoretical prediction is consistent with the actual observation.