Paper deep dive
A Coulomb Particle Model for Learning Kernel Attention in Transformers
Masoud Badiei Khuzani, Sharath Honnaiah, Atiq Islam, Alex Cozzi, Abraham Bagherjeiran
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Randomized features provide a scalable approximation to kernel machines, but their performance depends strongly on the choice of feature distribution. We propose a particle-based method that learns this distribution by optimizing kernel-target alignment while regularizing particles with a Riesz/Coulomb repulsive potential. The resulting Hamiltonian yields diverse, task-adaptive random features and admits a mean-field description through a McKean--Vlasov equation. We instantiate the method in linearized Transformer attention by learning positive random-feature maps in a first alignment phase, then freezing the kernel and training the remaining network parameters with cross-entropy. Experiments on synthetic classification and sentence-level benchmarks show that learned kernelized attention can improve accuracy, calibration, and robustness for several feature maps while preserving linear-attention inference complexity.
Tags
Links
- Source: https://arxiv.org/abs/2607.23869v1
- Canonical: https://arxiv.org/abs/2607.23869v1
Trouble viewing inline? Open PDF directly →
Full Text
147,056 characters extracted from source content.
Expand or collapse full text
pdftitle=A Coulomb Particle Model for Learning Kernel Attention in Transformers, pdfauthor=Masoud Badiei Khuzani, Atiq Islam, Alex Cozzi, Abraham Bagherjeiran, hypertexnames=false Badiei Khuzani @ebay.com Inc., San Jose, CA, USA and Honnaiah @ebay.com Inc., San Jose, CA, USA and Islam @ebay.com Inc., San Jose, CA, USA and Cozzi @ebay.com Inc., San Jose, CA, USA and Bagherjeiran @ebay.com Inc., San Jose, CA, USA A Coulomb Particle Model for Learning Kernel Attention in Transformers Abstract Randomized features provide a scalable approximation to kernel machines, but their performance depends strongly on the choice of feature distribution. We propose a particle-based method that learns this distribution by optimizing kernel-target alignment while regularizing particles with a Riesz/Coulomb repulsive potential. The resulting Hamiltonian yields diverse, task-adaptive random features and admits a mean-field description through a McKean–Vlasov equation. We instantiate the method in linearized Transformer attention by learning positive random-feature maps in a first alignment phase, then freezing the kernel and training the remaining network parameters with cross-entropy. Experiments on synthetic classification and sentence-level benchmarks show that learned kernelized attention can improve accuracy, calibration, and robustness for several feature maps while preserving linear-attention inference complexity. 1 Introduction Kernel methods are a principled way to encode nonlinear similarity, but classical training and prediction can scale poorly with the number of samples. Random Fourier features alleviate this problem by replacing an implicit kernel with an explicit randomized map (Rahimi and Recht, 2007). Nevertheless, the kernel and its associated feature distribution are typically fixed before any labels are observed—a choice that can dominate downstream performance when the appropriate similarity structure is unknown. We study a supervised kernel-learning procedure that keeps the computational advantages of random features while learning the feature distribution from labels. The key idea is to view the random features as interacting particles. A target-alignment objective attracts particles toward features that explain the labels, while a Riesz/Coulomb repulsive energy prevents collapse and encourages diversity. This gives a concrete optimization algorithm, a statistical-mechanics interpretation, and a direct route to kernelized attention in Transformers. Our contributions are: (i) a particle Hamiltonian for target-aligned random-feature learning; (i) a Langevin optimization procedure with a mean-field continuum limit; (i) an instantiation for positive random-feature attention that preserves linear-time inference; and (iv) empirical evidence that the learned kernels improve several accuracy and calibration metrics on synthetic and NLP benchmarks. Full related work, proofs, and additional ablations are in the supplementary material. 2 Particle kernel learning Let ϕ:×Ω→[−1,1]φ:X× →[-1,1] be a random feature map and let μ be a probability measure over feature parameters. The induced kernel is Kμ(,′)=∼μ[ϕ()ϕ(′)]. K_μ( x, x )=E_ ω μ[ _ ω( x) _ ω( x )]. (2.1) Given training data (i,yi)i=1n\( x_i,y_i)\_i=1^n, we learn μ by maximizing kernel-target alignment (Cristianini et al., 2002). With particles N=(1,…,N) _N=( ω_1,…, ω_N) and empirical measure μN=N−1∑k=1Nδk _N=N^-1 _k=1^N _ ω_k, this becomes the finite-dimensional Hamiltonian ℋN(N)=−1n(n−1)∑i≠jyiyj1N∑k=1Nϕk(i)ϕk(j)+λN,s(N), _N( _N)=- 1n(n-1) _i≠ jy_iy_j 1N _k=1^N _ ω_k( x_i) _ ω_k( x_j)+ _N,s( _N), (2.2) where N,s(N)=12N(N−1)∑k≠ℓgs(k−ℓ),gs()=‖2−s,s>0,−log‖2,s=0. _N,s( _N)= 12N(N-1) _k≠ g_s( ω_k- ω_ ), g_s( ω)= cases\| ω\|_2^-s,&s>0,\\ - \| ω\|_2,&s=0. cases (2.3) The first term rewards features whose empirical kernel aligns with labels; the second term spreads particles across Ω . For random Fourier features, this energy reduces to a trigonometric alignment objective involving cos(N⊤) ( X _N ) and sin(N⊤) ( X _N ), which makes the objective differentiable and easy to optimize. algorithm2e [!t] Projected Langevin feature learning k0k=1N\ ω_k^0\_k=1^N, step η, inverse temperature β, threshold δ empirical feature law μN _N and kernel K K Initialize μN0=N−1∑kδk0 _N^0=N^-1 _k _ ω_k^0 (μNm,μNm−1)<δD( _N^m, _N^m-1)<δ =1,…,Nk=1,…,N Draw km∼(0,I) ξ_k^m (0,I) and set km+1=Ω(km−ηN∇kℋN(Nm)+2η/βkm) ω_k^m+1=P_ ( ω_k^m-η N _ ω_kH_N( _N^m)+ 2η/β\, ξ_k^m) μNm+1=N−1∑kδkm+1 _N^m+1=N^-1 _k _ ω_k^m+1 Return K^(,′)=D−1∑k=1Dϕk()ϕk(′) K( x, x )=D^-1 _k=1^D _ ω_k( x) _ ω_k( x ) 3 Kernelized attention Self-attention can be written as a normalized kernel smoother (Vaswani et al., 2017; Nadaraya, 1964; Watson, 1964): i()=∑j=0ℓK(i,j)j∑j=0ℓK(i,j),K(,)=exp(⊤/d). a_i( X)= _j=0 K( q_i, k_j) v_j _j=0 K( q_i, k_j), K( q, k)= ( q k/ d). (3.1) Replacing K by a positive random-feature kernel KN(,)=ϕN()⊤ϕN()K_ _N( q, k)= φ_ _N( q) φ_ _N( k) gives the linearized attention estimator ^i()=ϕN(i)⊤(∑jϕN(j)j⊤)ϕN(i)⊤(∑jϕN(j)). a_i( X)= φ_ _N( q_i) ( _j φ_ _N( k_j) v_j ) φ_ _N( q_i) ( _j φ_ _N( k_j) ). (3.2) This preserves the normalized kernel-smoothing form but reduces the per-head sequence-length dependence from quadratic to linear in ℓ , up to the feature dimension (Katharopoulos et al., 2020a; Choromanski et al., 2021; Peng et al., 2021). In Phase A, we learn the feature particles N _N by alignment on sequence representations. In Phase B, we freeze N _N and train the remaining Transformer parameters using cross-entropy. 4 Theoretical results: mean-field limit and large deviation principle The particle view gives both a continuum training dynamics and an equilibrium concentration result. Define the alignment-induced potential V()=−1n(n−1)∑i≠jyiyjϕ(i)ϕ(j), V_D( ω)=- 1n(n-1) _i≠ jy_iy_j _ ω( x_i) _ ω( x_j), (4.1) and the continuum energy ℰs(μ)=∫V()μ()+λ2∬gs(−′)μ()μ(′). _s(μ)= V_D( ω)\,dμ( ω)+ λ2 g_s( ω- ω )\,dμ( ω)dμ( ω ). (4.2) The finite-particle optimizer has the following continuum training law. Theorem 4.1 (Projected-particle McKean–Vlasov mean-field limit). Let Ω⊂ℝd ^d be compact and convex with C2C^2 boundary and outward normal n. Assume V∈C2(Ω¯)V_D∈ C^2( ), a Lipschitz regularized interaction drift, exchangeable μ0N⇒ρ0d _0^N _0d ω, and a vanishing projected-Euler error. If km+1=Ω¯(km−ηNN∇kℋN(Nm)+2ηN/βkm),ηN↓0, ω_k^m+1=P_ \! ( ω_k^m- _NN _ ω_kH_N( _N^m)+ 2 _N/β\, ξ_k^m ), _N 0, (4.3) then μtN=N−1∑kδk⌊t/ηN⌋⇒ρtd _t^N=N^-1 _k _ ω_k t/ _N _td ω in probability, uniformly on compact time intervals. With Ut=V+λgs∗ρtU_t=V_D+λ g_s _t, ∂tρt=∇⋅(ρt∇Ut)+β−1Δρt, _t _t=∇\!·( _t∇ U_t)+β^-1 _t, (4.4) ρt|t=0=ρ0,(ρt∇Ut+β−1∇ρt)⋅=0on ∂Ω, _t|_t=0= _0, ( _t∇ U_t+β^-1∇ _t)\!· n=0 ∂ , (4.5) equivalently ∂ρt+βρt∂Ut=0 _ n _t+β _t _ nU_t=0. Moreover ρt≥0 _t≥ 0 and ∫Ωρt=1 _ _t=1. The proof is provided in the supplementary material. The two parts of UtU_t have distinct roles. The data potential V_D pulls mass toward features that reduce the empirical objective, while the interaction term spreads mass according to the regularized Riesz geometry and prevents all features from collapsing onto the same locations. The diffusion term β−1Δρtβ^-1 _t is the continuum trace of the Langevin noise: at finite temperature it encourages exploration and contributes an entropic regularization; in the zero-temperature limit the equation reduces to the deterministic transport law driven by −∇Ut-∇ U_t. Equivalently, ∂tρt=∇⋅(ρt∇(Ut+β−1logρt)), _t _t=∇\!· ( _t∇ (U_t+β^-1 _t ) ), so (4.5) is the Wasserstein gradient flow of ℱβ(ρ)=defℰs(ρ)+β−1Ent(ρ∣ℓ).F_β(ρ) =E_s(ρ)+β^-1Ent(ρ ). For smooth positive solutions the no-flux condition removes the boundary contribution and gives the dissipation identity dtℱβ(ρt)=−∫Ωρt|∇(Ut+β−1logρt)|2≤0. ddtF_β( _t)=- _ _t |∇ (U_t+β^-1 _t ) |^2d ω≤ 0. This identity is useful conceptually: training decreases the continuum free energy, and the only stationary points are self-consistent Gibbs densities of the form ρ∞()=1Z∞exp−β[V()+λ(gs∗ρ∞)()], _∞( ω)= 1Z_∞ \-β[V_D( ω)+λ(g_s _∞)( ω)]\, with the boundary condition inherited from projection. Theorem 4.2 (Large deviations for learned feature measures). Let Ω⊂ℝd ^d be bounded and let s>ds>d. If N∼Problem 11Problem 1Problem 1._N,β_N and =μNN-1∑=k1Nδωk, then μN satisfies an LDP on P(Ω). When →/βN1, the speed is =rNN and(4.6)4.6=Js(μ)-+Es(μ)Ent(∣μℓ)infν+Es(ν)Ent(∣νℓ),=Js(μ)-+Es(μ)Ent(∣μℓ)infν+Es(ν)Ent(∣νℓ),where =Ent(∣μℓ)∫log(/dμdℓ)dμ for ≪μℓ. When →/βN∞, the speed is =rNβN and(4.7)4.7=Js(μ)-Es(μ)infνEs(ν).=Js(μ)-Es(μ)infνEs(ν).For Borel ⊂AP(Ω), the standard LDP hold with rate Js and speed rN.ThesupplementprovesthisbycombiningSanov′stheoremwithVaradhan′slemmaaftertruncatingthesingularRieszinteraction.Consequentlythelearnedkernelconcentratesaroundthetask−adaptivevariationalkernel;iftherelevantratehasauniqueminimizerμs⋆,thenKμN(,′)→Kμs⋆(,′)exponentiallyatspeedNorβN. _N _N, _N$ and $ _N=N^-1 _k=1^N _ ω_k$, then $\ _N\$ satisfies an LDP on $P( )$. When $ _N/N→ 1$, the speed is $r_N=N$ and \@@amsalignJ_s(μ)=E_s(μ)+Ent(μ )- _ν\E_s(ν)+Ent(ν )\, where $Ent(μ )= (dμ/d )dμ$ for $μ $. When $ _N/N→∞$, the speed is $r_N= _N$ and \@@amsalignJ_s(μ)=E_s(μ)- _νE_s(ν). For Borel $A ( )$, the standard LDP hold with rate $J_s$ and speed $r_N$. theoremThesupplementprovesthisbycombiningSanov stheoremwithVaradhan slemmaaftertruncatingthesingularRieszinteraction.Consequentlythelearnedkernelconcentratesaroundthetask-adaptivevariationalkernel;iftherelevantratehasauniqueminimizer$ _s $,then$K_ _N( x, x )→ K_ _s ( x, x )$exponentiallyatspeed$N$or$ _N$. 5 Experiments We evaluate on synthetic classification and SST-2, QQP, and Rotten Tomatoes. Synthetic results show that optimized particles outperform fixed random Fourier features and the importance-sampling baseline under moderate noise. The accuracy gap widens with feature budget: at 64 features per head, optimized particles retain roughly 90% of the performance achieved at 256, while fixed random features degrade more sharply. NLP experiments use a two-layer Transformer encoder with hidden size 128, two heads, 256 random features per head, BERT tokenization, batch size 64, and the two-phase alignment-then-cross-entropy protocol. Training uses AdamW with linear warmup and cosine decay; the alignment phase runs for 20% of total steps before switching to cross-entropy fine-tuning. All results are averaged over three random seeds. We report the full detailed result tables for the three text-classification tasks in the main paper. For all rows in the lower block of each table, the attention projection is constrained to Q=K=reshape(x)Q=K=reshape(x), so vanilla and learned-kernel rows are matched feature-map comparisons. The reported metrics include accuracy, micro/macro/weighted F1, ROC-AUC, PR-AUC, MCC, balanced accuracy, LogLoss, Brier score, ECE, and wall-clock training time. Table 1: SST-2 detailed metrics across attention variants. Bold marks the best value per column; lower is better for LogLoss, Brier, ECE, and train time. Model Acc↑ F1↑ μ F1M↑ F1w↑ ROC-AUC↑ PR-AUC↑ MCC↑ BalAcc↑ LogLoss↓ Brier↓ ECE↓ Train(sec)↓ Learned Wq,WkW_q,W_k baselines Vaswani-softmaxattn 0.8050 0.8050 0.8050 0.8051 0.8744 0.8856 0.6100 0.8050 0.5083 0.1498 0.3695 86.8 Vanilla-favor-WqWk 0.7833 0.7833 0.7832 0.7832 0.8587 0.8576 0.5676 0.7837 0.5472 0.1592 0.3584 96.5 Performer 0.7798 0.7798 0.7793 0.7791 0.8618 0.8766 0.5656 0.7810 0.5697 0.1666 0.3804 89.0 Linformer 0.7844 0.7844 0.7843 0.7843 0.8569 0.8673 0.5706 0.7850 0.5922 0.1646 0.3782 82.6 No learned Wq,WkW_q,W_k (Q=K=reshape(x)); feature-map variants Vanilla-favor 0.8050 0.8050 0.8050 0.8049 0.8841 0.8928 0.6124 0.8057 0.4362 0.1388 0.3420 80.5 Kernel-favor 0.7901 0.7901 0.7900 0.7901 0.8734 0.8834 0.5801 0.7899 0.4477 0.1459 0.3086 166.0 Vanilla-elu 0.8028 0.8028 0.8019 0.8017 0.8731 0.8803 0.6150 0.8042 0.4938 0.1505 0.3724 78.8 Kernel-elu 0.8108 0.8108 0.8107 0.8108 0.8879 0.8913 0.6216 0.8108 0.4502 0.1392 0.3609 164.3 Vanilla-softplus 0.8108 0.8108 0.8108 0.8108 0.8760 0.8837 0.6225 0.8112 0.4582 0.1429 0.3378 88.5 Kernel-softplus 0.8062 0.8062 0.8058 0.8060 0.8861 0.8864 0.6126 0.8057 0.5122 0.1470 0.3920 173.5 Vanilla-sigmoid2 0.7959 0.7959 0.7957 0.7958 0.8786 0.8881 0.5916 0.7956 0.5025 0.1487 0.3744 82.6 Kernel-sigmoid2 0.7924 0.7924 0.7924 0.7924 0.8785 0.8880 0.5863 0.7929 0.4915 0.1516 0.3790 170.1 Vanilla-softmaxfeat 0.8005 0.8005 0.8003 0.8004 0.8683 0.8728 0.6008 0.8001 0.5163 0.1530 0.3621 89.5 Kernel-softmaxfeat 0.8085 0.8085 0.8085 0.8084 0.8858 0.8982 0.6181 0.8089 0.4381 0.1389 0.3526 171.4 Vanilla-cos2 0.8142 0.8142 0.8142 0.8141 0.8897 0.8945 0.6305 0.8148 0.4709 0.1395 0.3866 81.5 Kernel-cos2 0.8119 0.8119 0.8118 0.8119 0.8883 0.8976 0.6237 0.8118 0.4422 0.1366 0.3557 165.0 Vanilla-porf-softplus 0.8142 0.8142 0.8142 0.8142 0.8883 0.8983 0.6285 0.8143 0.4506 0.1380 0.3569 85.7 Kernel-porf-softplus 0.8188 0.8188 0.8188 0.8188 0.8897 0.8967 0.6376 0.8189 0.4561 0.1381 0.3734 162.8 Table 2: QQP detailed metrics across attention variants. Bold marks the best value per column; lower is better for LogLoss, Brier, ECE, and train time. Model Acc↑ F1↑ μ F1M↑ F1w↑ ROC-AUC↑ PR-AUC↑ MCC↑ BalAcc↑ LogLoss↓ Brier↓ ECE↓ Train(sec)↓ Learned Wq,WkW_q,W_k baselines Vaswani-softmaxattn 0.8005 0.8005 0.7949 0.8038 0.8959 0.8266 0.6058 0.8133 0.4112 0.1350 0.4314 605.9 Vanilla-favor-WqWk 0.7683 0.7683 0.7346 0.7595 0.8313 0.7619 0.4852 0.7248 0.4863 0.1591 0.4383 539.7 Performer 0.7751 0.7751 0.7472 0.7693 0.8402 0.7691 0.5030 0.7388 0.4732 0.1547 0.4360 554.1 Linformer 0.8077 0.8077 0.7944 0.8082 0.8828 0.8157 0.5890 0.7960 0.4153 0.1337 0.4541 557.1 No learned Wq,WkW_q,W_k (Q=K=reshape(x)); feature-map variants Vanilla-favor 0.7626 0.7626 0.7323 0.7560 0.8168 0.7443 0.4741 0.7242 0.5028 0.1637 0.4182 567.8 Kernel-favor 0.7608 0.7608 0.7319 0.7551 0.8158 0.7438 0.4712 0.7246 0.5016 0.1638 0.4106 1197.1 Vanilla-elu 0.7797 0.7797 0.7553 0.7757 0.8458 0.7772 0.5155 0.7485 0.4659 0.1518 0.4386 585.1 Kernel-elu 0.7831 0.7831 0.7648 0.7821 0.8533 0.7863 0.5299 0.7625 0.4598 0.1496 0.4474 1110.3 Vanilla-softplus 0.7786 0.7786 0.7522 0.7735 0.8442 0.7747 0.5116 0.7444 0.4656 0.1522 0.4297 633.3 Kernel-softplus 0.7860 0.7860 0.7651 0.7836 0.8564 0.7891 0.5323 0.7603 0.4578 0.1481 0.4582 1100.7 Vanilla-sigmoid2 0.7772 0.7772 0.7580 0.7760 0.8439 0.7728 0.5165 0.7554 0.4664 0.1529 0.4190 591.9 Kernel-sigmoid2 0.7852 0.7852 0.7666 0.7840 0.8549 0.7885 0.5338 0.7639 0.4537 0.1479 0.4415 1156.0 Vanilla-softmaxfeat 0.7801 0.7801 0.7527 0.7744 0.8471 0.7772 0.5143 0.7441 0.4625 0.1512 0.4365 579.9 Kernel-softmaxfeat 0.7909 0.7909 0.7753 0.7909 0.8679 0.7956 0.5506 0.7753 0.4354 0.1418 0.4371 1072.6 Vanilla-cos2 0.7746 0.7746 0.7533 0.7724 0.8422 0.7715 0.5082 0.7492 0.4688 0.1536 0.4249 650.2 Kernel-cos2 0.7842 0.7842 0.7590 0.7795 0.8525 0.7842 0.5245 0.7511 0.4621 0.1497 0.4526 1254.0 Vanilla-porf_softplus 0.7782 0.7782 0.7531 0.7739 0.8434 0.7751 0.5117 0.7461 0.4685 0.1528 0.4401 588.4 Kernel-porf_softplus 0.7841 0.7841 0.7632 0.7817 0.8536 0.7855 0.5284 0.7585 0.4547 0.1483 0.4419 1110.4 Table 3: Rotten Tomatoes detailed metrics across attention variants. Bold marks the best value per column; lower is better for LogLoss, Brier, ECE, and train time. Model Acc↑ F1↑ μ F1M↑ F1w↑ ROC-AUC↑ PR-AUC↑ MCC↑ BalAcc↑ LogLoss↓ Brier↓ ECE↓ Train(sec)↓ Learned Wq,WkW_q,W_k baselines Vaswani-softmaxattn 0.6801 0.6801 0.6793 0.6793 0.7492 0.7522 0.3620 0.6801 0.6396 0.2138 0.2751 12.6 Vanilla-favor-WqWk 0.6764 0.6764 0.6758 0.6758 0.7262 0.7162 0.3539 0.6764 0.7265 0.2327 0.3052 13.0 Performer 0.6642 0.6642 0.6616 0.6616 0.7262 0.7075 0.3333 0.6642 0.7477 0.2347 0.3009 14.4 Linformer 0.6614 0.6614 0.6612 0.6612 0.7133 0.7062 0.3230 0.6614 0.6614 0.2260 0.2509 13.7 No learned Wq,WkW_q,W_k (Q=K=reshape(x)); feature-map variants Vanilla-favor 0.7083 0.7083 0.7066 0.7066 0.7732 0.7708 0.4214 0.7083 0.5840 0.1977 0.2229 13.7 Kernel-favor 0.6848 0.6848 0.6848 0.6848 0.7597 0.7479 0.3696 0.6848 0.5823 0.1997 0.2193 26.3 Vanilla-elu 0.6895 0.6895 0.6835 0.6835 0.7844 0.7928 0.3943 0.6895 0.5809 0.1995 0.2435 12.7 Kernel-elu 0.7148 0.7148 0.7140 0.7140 0.7650 0.7586 0.4320 0.7148 0.5860 0.1991 0.2594 26.1 Vanilla-softplus 0.7036 0.7036 0.7023 0.7023 0.7695 0.7580 0.4107 0.7036 0.5838 0.1977 0.2176 12.6 Kernel-softplus 0.7017 0.7017 0.7017 0.7017 0.7743 0.7677 0.4034 0.7017 0.5766 0.1959 0.2517 23.8 Vanilla-sigmoid2 0.7026 0.7026 0.7022 0.7022 0.7714 0.7675 0.4065 0.7026 0.5758 0.1958 0.2408 12.7 Kernel-sigmoid2 0.6876 0.6876 0.6831 0.6831 0.7562 0.7538 0.3864 0.6876 0.6118 0.2091 0.2733 24.0 Vanilla-softmaxfeat 0.6782 0.6782 0.6733 0.6733 0.7664 0.7703 0.3678 0.6782 0.6065 0.2071 0.2610 12.5 Kernel-softmaxfeat 0.7167 0.7167 0.7163 0.7163 0.7779 0.7724 0.4345 0.7167 0.5697 0.1929 0.2527 24.5 Vanilla-cos2 0.6876 0.6876 0.6855 0.6855 0.7622 0.7586 0.3803 0.6876 0.5897 0.2015 0.2154 12.4 Kernel-cos2 0.7092 0.7092 0.7091 0.7091 0.7676 0.7680 0.4185 0.7092 0.5836 0.1982 0.2594 26.7 Vanilla-porf_softplus 0.6839 0.6839 0.6832 0.6832 0.7649 0.7637 0.3692 0.6839 0.5858 0.1996 0.2343 12.6 Kernel-porf_softplus 0.6989 0.6989 0.6986 0.6986 0.7702 0.7700 0.3984 0.6989 0.5750 0.1961 0.2378 24.1 Tables 1, 2, and 3 show that kernel learning is feature-map dependent but often improves both discrimination and proper scoring metrics. On SST-2, Kernel-PORF-softplus gives the highest accuracy and MCC, while Kernel-FAVOR gives the lowest ECE. On QQP, Kernel-softmaxfeat is the strongest Q=KQ=K model and reduces LogLoss and Brier score relative to its vanilla counterpart. On Rotten Tomatoes, Kernel-softmaxfeat gives the best accuracy, MCC, LogLoss, and Brier score. The training-time column reflects the one-time alignment phase; after particles are learned, the serving-time attention formula remains linear in sequence length. Additional experiment details. Inputs use the bert-base-uncased WordPiece tokenizer, are padded or truncated to 128 tokens, and use mean pooling over final token representations. The Transformer has two encoder layers, hidden size 128, two attention heads, feed-forward dimension 256, dropout 0.1, and 256 random features per head. Phase A optimizes only the feature particles N _N using the alignment objective with Langevin noise, log-repulsion, norm projection, and early stopping; Phase B freezes the learned particles and trains the remaining classifier with cross-entropy. Thus the extra cost appears during training, while inference uses the same linear-attention structure as the corresponding vanilla feature map. 6 Conclusion We introduced a Coulomb/Riesz particle model for supervised random-feature learning and applied it to kernelized Transformer attention. The framework links label alignment, repulsive regularization, Langevin optimization, and equilibrium concentration in a single formulation. Empirically, the method improves several random-feature attention variants and provides a practical way to learn task-adaptive linear-attention kernels. LLM usage statement. The authors used large language model (LLM) tools for writing assistance and code development. All LLM-assisted content was reviewed, verified, and edited by the authors, who take full responsibility for the correctness, originality, citations, proofs, experiments, figures, and final content of this paper. References Bach et al. (2004) Francis R Bach, Gert RG Lanckriet, and Michael I Jordan. Multiple kernel learning, conic duality, and the SMO algorithm. In Proceedings of the twenty-first international conference on Machine learning, page 6, 2004. Chapelle et al. (2002) Olivier Chapelle, Vladimir Vapnik, Olivier Bousquet, and Sayan Mukherjee. Choosing multiple parameters for support vector machines. Machine Learning, 46(1–3):131–159, 2002. 10.1023/A:1012450327387. Chen et al. (2021) Yifan Chen, Qi Zeng, Heng Ji, and Yun Yang. Skyformer: Remodel self-attention with gaussian kernel and nyström method. In Advances in Neural Information Processing Systems, 2021. Choromanski et al. (2017) Krzysztof Choromanski, Mark Rowland, and Adrian Weller. Structured orthogonal random features. In AISTATS, 2017. Choromanski et al. (2020) Krzysztof Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tamás Sarlós, Peter Hawkins, Jared Q. Davis, Afroz Mohiuddin, Łukasz Kaiser, David Belanger, Lucy J. Colwell, and Adrian Weller. Rethinking attention with performers. arXiv preprint arXiv:2009.14794, 2020. URL https://arxiv.org/abs/2009.14794. Choromanski et al. (2021) Krzysztof Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tam’as Sarl’os, Peter Hawkins, Jared Q. Davis, Afroz Mohiuddin, Łukasz Kaiser, David Belanger, Lucy J. Colwell, and Adrian Weller. Rethinking attention with performers. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=Ua6zuk0WRH. Clevert et al. (2015) Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter. Fast and accurate deep network learning by exponential linear units (elus). arXiv:1511.07289, 2015. Cortes et al. (2009) Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh. Learning non-linear combinations of kernels. Advances in neural information processing systems, 22, 2009. Cortes et al. (2010) Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh. Two-stage learning kernel algorithms. In Proceedings of the 27th International Conference on Machine Learning (ICML-10), pages 239–246, 2010. Cristianini et al. (2002) Nello Cristianini, John Shawe-Taylor, André Elisseeff, and Jaswinder Singh Kandola. On kernel-target alignment. In Advances in neural information processing systems, volume 14, pages 367–373, 2002. Dembo and Zeitouni (2009) Amir Dembo and Ofer Zeitouni. Large Deviations Techniques and Applications, volume 38 of Stochastic Modelling and Applied Probability. Springer, Berlin, Heidelberg, 2nd edition edition, 2009. ISBN 978-3-642-03310-0. 10.1007/978-3-642-03311-7. Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. Bert: Pre-training of deep bidirectional transformers for language understanding. In NAACL-HLT, 2019. Dugas et al. (2001) Charles Dugas, Yoshua Bengio, François Bélisle, Claude Nadeau, and René Garcia. Incorporating second-order functional knowledge for better option pricing. In NeurIPS, 2001. Ellis (2005) Richard S. Ellis. Entropy, Large Deviations, and Statistical Mechanics. Springer Monographs in Mathematics. Springer, June 2005. 10.1007/0-387-28537-0. Kandola et al. (2002) Jaswinder Singh Kandola, John Shawe-Taylor, and Nello Cristianini. Optimizing kernel alignment over combinations of kernels. Technical report, University of Southampton, 2002. Katharopoulos et al. (2020a) Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret. Transformers are rnns: Fast autoregressive transformers with linear attention. In Proceedings of the 37th International Conference on Machine Learning (PMLR 119), 2020a. Katharopoulos et al. (2020b) Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret. Transformers are RNNs: Fast autoregressive transformers with linear attention. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 5156–5165. PMLR, 2020b. URL https://proceedings.mlr.press/v119/katharopoulos20a.html. Lanckriet et al. (2004) Gert RG Lanckriet, Nello Cristianini, Peter Bartlett, Laurent El Ghaoui, and Michael I Jordan. Learning the kernel matrix with semidefinite programming. Journal of Machine Learning Research, 5:27–72, 2004. Lauriola and Aiolli (2020) Ivano Lauriola and Fabio Aiolli. MKLpy: a python-based framework for multiple kernel learning. arXiv preprint arXiv:2007.09982, 2020. Luo et al. (2021) Shengjie Luo, Shanda Li, Tianle Cai, Di He, Dinglan Peng, Shuxin Zheng, Guolin Ke, Liwei Wang, and Tie-Yan Liu. Stable, fast and accurate: Kernelized attention with relative positional encoding. arXiv:2106.12566, 2021. Micchelli and Pontil (2005) Charles A Micchelli and Massimiliano Pontil. Learning the kernel function via regularization. Journal of Machine Learning Research, 6:1099–1125, 2005. Nadaraya (1964) E. A. Nadaraya. On estimating regression. Theory of Probability and its Applications, 1964. Ober et al. (2021) Stefan Ober, Carl Edward Rasmussen, and Mark Van der Wilk. The promises and pitfalls of deep kernel learning. Journal of Machine Learning Research, 22(179):1–65, 2021. Pang and Lee (2005) Bo Pang and Lillian Lee. Seeing stars: Exploiting class relationships for sentiment categorization with respect to rating scales. In Proceedings of the 43rd Annual Meeting of the Association for Computational Linguistics (ACL), pages 115–124. Association for Computational Linguistics, 2005. Peng et al. (2021) Hao Peng, Nikolaos Pappas, Dani Yogatama, Roy Schwartz, Noah A. Smith, and Lingpeng Kong. Random feature attention. arXiv:2103.02143, 2021. Rahimi and Recht (2009) Ali Rahimi and Ben Recht. Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning. In NeurIPS, 2009. Rahimi and Recht (2007) Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, volume 20, pages 1177–1184, 2007. Rakotomamonjy et al. (2008) Alain Rakotomamonjy, Francis R Bach, St’ephane Canu, and Yves Grandvalet. SimpleMKL. Journal of Machine Learning Research, 9:2491–2521, 2008. Sinha and Duchi (2016) Samarth Sinha and John C Duchi. Learning kernels with random features. In Advances in Neural Information Processing Systems, pages 1298–1306, 2016. Socher et al. (2013) Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D. Manning, Andrew Y. Ng, and Christopher Potts. Recursive deep models for semantic compositionality over a sentiment treebank. In Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 1631–1642, 2013. Tsai et al. (2019) Yao-Hung Hubert Tsai, Shaojie Bai, Makoto Yamada, Louis-Philippe Morency, and Ruslan Salakhutdinov. Transformer dissection: A unified understanding of transformer’s attention via the lens of kernel. In EMNLP-IJCNLP, 2019. Varadhan (2016) S. R. S. Varadhan. Large Deviations, volume 27 of Courant Lecture Notes. American Mathematical Society, 2016. ISBN 978-1-4704-2580-1. Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Information Processing Systems, 2017. Wang et al. (2019) Alex Wang, Amanpreet Singh, Julian Michael, et al. GLUE: A multi-task benchmark and analysis platform for natural language understanding. In Proceedings of the 7th International Conference on Learning Representations (Workshop Track), 2019. URL https://openreview.net/forum?id=rJ4km2R5t7. Watson (1964) G. S. Watson. Smooth regression analysis. Sankhyā, Series A, 1964. Wilson et al. (2016) Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial intelligence and statistics, pages 370–378. PMLR, 2016. Xiong et al. (2021) Yunyang Xiong, Zhanpeng Zeng, Rudrasis Chakraborty, Mingxing Tan, Glenn Fung, Yin Li, and Vikas Singh. Nyströmformer: A nyström-based algorithm for approximating self-attention. In AAAI Conference on Artificial Intelligence, 2021. Yu et al. (2016) Felix X. Yu, Ananda Theertha Suresh, Krzysztof Choromanski, Daniel Holtmann-Rice, and Sanjiv Kumar. Orthogonal random features. In NeurIPS, 2016. Supplementary Material Appendix A Related work To discuss related work, we first describe the kernel selection problem in the context of supervised learning problem. Consider a set of n feature vectors and labels (i,yi)i=1n,i∈,yi∈\( x_i,y_i)\_i=1^n, x_i ,y_i . We have a loss function L:×→L: × → , where L(⋅,y)L(·,y) is convex for y∈y , and a reproducing kernel Hilbert space (RKHS) of functions ℱF with kernel K. The ℓ2 _2-regularized optimization problem that underlies the learning task of finding a function f∈ℱf is as follows Primal:minf∈ℱ1n∑i=1nL(f(i),yi)+λ2‖f‖ℋ2,Dual:max∈ℝn−1n∑i=1nL∗(αi,yi)−12λ⊤, : _f 1n _i=1^nL(f( x_i),y_i)+ λ2\|f\|_H^2, : _ α ^n- 1n _i=1^nL^*( _i,y_i)- 12λ α K α, (A.1) where ∥⋅∥ℋ\|·\|_H is the Hilbert space norm, ∈ℝn α ^n are dual variables, L∗(α,y)=supz∈ℝαz−L(z,y)L^*(α,y)= _z \α z-L(z,y) \ is the Fenchel conjugate of the loss function L, and =def[Kij]∈n×n K =[K_ij]∈ ^n× n with Kij=defK(i,j)K_ij =K( x_i, x_j) is the kernel matrix. While the kernel matrix K could, in principle, be optimized jointly with the dual variables, much of the literature instead focuses on approaches that decouple kernel learning from the estimation of f∈ℱf . A common strategy is to first construct or adapt the kernel—often by maximizing kernel–target alignment, which quantifies the similarity between the kernel and the target—before solving the regularized risk minimization problem (see, e.g., Cortes et al. (2010); Cristianini et al. (2002); Kandola et al. (2002); Lanckriet et al. (2004)). This alignment is formulated as the following optimization problem: maxK∈1n(n−1)∑1≤i≠j≤nyiyjK(i,j), _K 1n(n-1) _1≤ i =j≤ ny_iy_jK( x_i, x_j), (A.2) where K denotes a predefined class of kernel functions. Several approaches have been proposed in the literature for defining the class of kernels K in kernel-based learning. In Table 4, we provide a summary of common kernel class choices, along with their computational and memory complexities, as well as the corresponding references. Among the proposed approaches, Sinha and Duchi (2016) stands out as a framework that seamlessly integrates with the random feature model, specifically by employing importance sampling of random features within the random feature-based kernel class. In contrast, we propose an alternative approach that yields improved performance. This enhancement is achieved by directly optimizing the distribution of random features within a particle optimization framework, as opposed to relying on importance sampling of random features. Table 4: Summary of kernel class K choices in the literature. Notes. n: number of training samples; m: number of base kernels; D: number of random features; d: dimensionality of the input feature vectors; L: number of parameters or layers in the nonlinear transformation. Kernel Class Definition Comp. Complexity Memory Complexity References Convex combination of base kernels =K=∑i=1mwiKi|wi≥0,∑i=1mwi=1K= \K= _i=1^mw_iK_i\, |\,w_i≥ 0,\ _i=1^mw_i=1 \ =(w1,…,wm)∈ℝm w=(w_1,…,w_m) ^m lies on the probability simplex. (nm2)O(nm^2) (nm)O(nm) Lanckriet et al. (2004); Bach et al. (2004) Linear combination of base kernels =K=∑i=1mwiKi|∈ℝmK= \K= _i=1^mw_iK_i\, |\, w ^m \ Allows negative weights; PSD constraints may be required. (nm2)O(nm^2) (nm)O(nm) Cortes et al. (2009) Nonlinear kernel combinations =K(,′)=σ(∑i=1mwiKi(,′)+b)K= \K( x, x )=σ ( _i=1^mw_iK_i( x, x )+b ) \ σ is a nonlinearity (e.g., ReLU, sigmoid); wi,b∈ℝw_i,b . (nmL)O(nmL) (nm)O(nm) Wilson et al. (2016); Ober et al. (2021) Parameterized kernels =Kγ(,′)=exp(−γ∥−′∥2)|γ∈ΓK= \K_γ( x, x )= (-γ\| x- x \|^2)\, |\,γ∈ \ Γ⊂ℝ+ _+ is a bounded interval over which γ is optimized. (n2d)O(n^2d) (n2)O(n^2) Chapelle et al. (2002) Random feature-based kernels =Kμ(,′)=∼μ[ϕ()ϕ(′)]|μ∈ℳK= \K_μ( x, x )=E_ ω μ[ _ ω( x) _ ω( x )]\, |\,μ \ ϕ() _ ω( x) is a random feature map; ℳM is a set of distributions over ω. (nD)O(nD) (nD)O(nD) Sinha and Duchi (2016) SDP-based kernel learning =∈+n|tr()≤cK= \ K _+^n\, |\,tr( K)≤ c \ +nS_+^n denotes the set of n×n× n PSD matrices; c>0c>0 is a trace constraint. (n6)O(n^6) (n2)O(n^2) Lanckriet et al. (2004) Structured or hierarchical kernels =K=∑i=1mwiKi|∈structuredK= \K= _i=1^mw_iK_i\, |\,w _structured \ structuredW_structured encodes priors such as group sparsity or tree-structured dependencies. (nmlogm)O(nm m) (nm)O(nm) Rakotomamonjy et al. (2008); Micchelli and Pontil (2005) Mean-field kernel approach =Kμ(,′)=∼μ[ϕ()ϕ(′)]|μ∈ℳK= \K_μ( x, x )=E_ ω μ[ _ ω( x) _ ω( x )]\, |\,μ \ ϕ() _ ω( x) is a random feature map; ℳM is a set of distributions over ω. (nD2)O(nD^2) (nD)O(nD) This work Appendix B Proposed approach At a high level, we begin with a feature mapping to represent the kernel. Next, we learn a distribution that aligns this mapping with the labels using the kernel-target alignment (KTA) optimization formulated in Eq. (A.2). From this distribution, we sample random features, which are then used in a standard supervised learning framework. Specifically, let ϕ:×Ω→[−1,1]φ:X× →[-1,1], and μ denotes a probability measure on Ω . We define the kernel Kμ(i,j)=∼μ[ϕ(i)ϕ(j)], K_μ( x_i, x_j)=E_ ω μ[ _ ω( x_i) _ ω( x_j)], (B.1) where we used the shorthand notation ϕ(⋅)=defϕ(⋅;) _ ω(·) =φ(·; ω). We optimize kernel KμK_μ over all distributions μ in some (large, nonparametric) set ℳM of possible distributions on random features supμ∈ℳ1n(n−1)∑0≤i≠j≤nyiyj∼μ[ϕ(i)ϕ(j)]. _μ 1n(n-1) _0≤ i =j≤ ny_iy_jE_ ω μ[ _ ω( x_i) _ ω( x_j)]. (B.2) We consider independent, identically distributed (i.i.d.) samples or particles 1,…,N∼i.i.d.μ ω_1,…, ω_N .i.d. μ.111To distinguish between training samples and random feature samples, we refer to the latter as particles. Define the configuration of particles N=def(1,…,N)∈ΩN _N =( ω_1,…, ω_N)∈ ^N, and the associated empirical distribution μN(N)=1N∑k=1Nδk() _N( _N)= 1N _k=1^N _ ω_k( ω), where δk(⋅) _ ω_k(·) is Dirac’s delta function concentrated at k ω_k. We consider the following regularized optimization problem to estimate the expectation term in Eq. (B.2) by optimizing the samples of the distribution, where the expectation is substituted by the Monte Carlo sample average approximation. In particular, we consider the following Gibbs point process with the Hamiltonian infμ∈ℳNℋN(N)=defℰN(N)+λN,s(N), _μ _NH_N( _N) =E_N( _N)+ _N,s( _N), (B.3) where λ>0λ>0 controls the strength of the interaction, the population loss function is approximated as follows: ℰN(N)=def−1n(n−1)∑0≤i≠j≤nyiyj1N∑k=1Nϕk(i)ϕk(j). _N( _N) =- 1n(n-1) _0≤ i =j≤ ny_iy_j 1N _k=1^N _ ω_k( x_i) _ ω_k( x_j). (B.4) Specifically, in the derivation of the empirical loss in Eq. (B.4), we replace the expectation under μ by integration with respect to the empirical measure μN=def1N∑k=1Nδk _N = 1N _k=1^N _ ω_k: Kμ(i,j) K_μ( x_i, x_j) =∫Ωϕ(i)ϕ(j)μ(d) = _ _ ω( x_i)\, _ ω( x_j)\,μ(d ω) ≈∫Ωϕ(i)ϕ(j)μN(d)=KμN(i,j) ≈ _ _ ω( x_i)\, _ ω( x_j)\, _N(d ω)=K_ _N( x_i, x_j) =1N∑k=1Nϕk(i)ϕk(j). = 1N _k=1^N _ ω_k( x_i)\, _ ω_k( x_j). (B.5) The following lemma provides the explicit form of this energy function for the random Fourier feature models in Rahimi and Recht (2007): Lemma B.1 (Empirical KTA equals trigonometric energy with random bias). Let (i,yi)i=1n\( x_i,y_i)\_i=1^n with i∈ℝd x_i ^d, yi∈−1,+1y_i∈\-1,+1\. Let N=def(1,…,N)∈ΩN _N =( ω_1,…, ω_N)∈ ^N be the particle configuration and =def(b1,…,bN)⊤ b =(b_1,…,b_N) with bk∼i.i.d.Unif[0,2π]b_k .i.d. Unif[0,2π]. Define ϕ,b()=2cos(⊤+b) _ ω,b( x)= 2 ( ω x+b). Consider ℰN(N)=def−1n(n−1)∑0≤i≠j≤nyiyj1N∑k=1Nϕk,bk(i)ϕk,bk(j). _N( _N) =- 1n(n-1) _0≤ i≠ j≤ ny_iy_j\; 1N _k=1^N _ ω_k,b_k( x_i)\, _ ω_k,b_k( x_j). (B.6) Let ∈ℝn×d X ^n× d stack the i⊤ x_i , and write =N⊤+n⊤,=cos()∈ℝn×N,=sin()∈ℝn×N, A\;=\; X _N^\! + 1_n b , C= ( A) ^n× N, S= ( A) ^n× N, where n 1_n is the all-ones vector, and cos(⋅) (·) and sin(⋅) (·) are applied to each element of the matrix. Then, up to an additive constant independent of (N,)( _N, b), ℰN(N)≡−2n2‖⊤‖22, _N( _N)\;≡\;- 2n^2\, \| y C \|_2^2, (B.7) and, taking expectation over the random phases b, [ℰN(N)] _ b\! [E_N( _N) ] ≡−1Nn(n−1)(‖⊤cos(N⊤)‖22+‖⊤sin(N⊤)‖22) \;≡\;- 1Nn(n-1) ( \| y ( X _N^\! ) \|_2^2+ \| y ( X _N^\! ) \|_2^2 ) (B.8) =−1Nn(n−1)∑k=1N|∑i=1nyieik⊤i|2. =- 1N\,n(n-1) _k=1^N | _i=1^ny_i\,e^\,i\, ω_k^\! x_i |^2. (B.9) Furthermore, the regularization term in Eq. (B.3) captures the interaction energy of every sample k _k with all the other samples ℓ _ in an infinite configuration N,s(N)=12N(N−1)∑1≤k≠ℓ≤Ngs(k−ℓ), _N,s( _N)= 12N(N-1) _1≤ k = ≤ Ng_s( ω_k- ω_ ), (B.10) where gsg_s is the homogeneous potential of degree s, gs()=def‖2−s,for s∈(0,+∞],−log‖2,for s=0,−‖2−s,for s∈(−2,0).g_s( ω) = cases\| ω\|_2^-s,&for s∈(0,+∞],\\ - \| ω\|_2,&for s=0,\\ -\| ω\|_2^-s,&for s∈(-2,0). cases algorithm2e [t] Riesz/Coulomb Particles for Kernel Estimation in Kernel Methods of particles N; initial particle positions k0k=1N\ ω_k^0\_k=1^N; step size η; inverse temperature β; threshold δ; divergence (⋅,⋅)D(·,·); number of random feature samples D<ND<N kernel K^(,′) K( x, x ) μN0←N−1∑k=1Nδk0 _N^0← N^-1 _k=1^N _ ω_k^0 m←0m← 0 (μNm,μNm−1)<δD( _N^m, _N^m-1)<δ =1k=1 sample km∼(,d) ξ_k^m ( 0, I_d) ∇ℋN←N∇kmℋN(Nm) _N← N _ ω_k^mH_N( _N^m) km+1←Ω(km−η∇ℋN+2η/βkm) ω_k^m+1 _ \! ( ω_k^m-η _N+ 2η/β\, ξ_k^m ) μNm+1←N−1∑k=1Nδkm+1 _N^m+1← N^-1 _k=1^N _ ω_k^m+1 m←m+1m← m+1 Compute weights wk←exp[−βℋ(km)]/∑j=1Nexp[−βℋ(jm)]w_k← [- ( ω_k^m)]/ _j=1^N [- ( ω_j^m)] Sample D particles k∗k=1D\ ω_k \_k=1^D according to wkk=1N\w_k\_k=1^N K^(,′)←D−1∑k=1Dϕk∗()ϕk∗(′) K( x, x )← D^-1 _k=1^D _ ω_k ( x) _ ω_k ( x ) ^(,′) K( x, x ) Note that the Hamiltonian in Eq. (B.3) uses the standard mean-field normalization of the pairwise interaction. Under this normalization, the interaction energy is an average over particle pairs, and its large-N limit is the corresponding continuum Riesz energy. Therefore, we do not introduce an additional factor Ns/dN^s/d in either the loss term or the interaction term. The heuristic spacing N−1/dN^-1/d describes nearest-neighbor distances among N well-distributed particles in d dimensions, whereas the distance between a typical pair of independently sampled particles from a fixed limiting distribution remains O(1)O(1). Thus, multiplying the Hamiltonian by Ns/dN^s/d would correspond to a different, N-dependent choice of interaction strength rather than to the mean-field scaling considered here. We regard the coordinates in Ω as nondimensionalized; equivalently, any fixed length-scale normalization, such as replacing gs()g_s( ω) by gs(/ℓ0)g_s( ω/ _0), can be absorbed into the coupling parameter λ for s>0s>0, while in the logarithmic case it changes the energy only by an additive constant and hence leaves the induced dynamics unchanged. Viewing the particles as charged samples in Ω , the interaction in Eq. (B) corresponds to the Riesz potential for general s>0s>0, with the Coulomb/Newtonian case obtained at s=d−2s=d-2 for d≥3d≥ 3 and the logarithmic Coulomb kernel at s=0s=0 in d=2d=2. The sign convention is chosen so that the induced pairwise force is repulsive. These repulsive interactions act as a regularization mechanism, preventing the particles from collapsing into a single point mass (Dirac delta) and promoting their dispersion across the domain Ω . In particular, they encourage the support of the empirical distribution μ^N μ_N to approximate the support of the true underlying distribution Ω=supp(μ) =supp(μ), where supp(⋅)supp(·) denotes the support of a distribution. From a statistical learning theory perspective, such repulsive forces mitigate sample clustering, thereby enhancing the expressiveness of the kernel. By maintaining spatial diversity among the samples, the kernel can better capture the underlying structure of the data. Numerical simulations confirm that this regularization indeed improves both the stability of the particle system and the generalization performance of kernel-based models. B.1 Langevin dynamics for efficiently solving Eq. (B.3) We optimize the positions of the samples in Eq. (B.3) using Langevin dynamics. From an optimization perspective, we model the evolution of the particle system over time steps m=0,1,…,T−1m=0,1,…,T-1 according to the following Langevin update rule: km+1=Ω(km−ηN∇kmℋN(Nm)+2ηβN⋅km),k=1,2,⋯,N, ω_k^m+1=P_ ( ω_k^m-η\,N _ ω_k^mH_N( _N^m)+ 2η _N· ξ^m_k ), k=1,2,·s,N, (B.11) where η=defη(n,N)η =η(n,N) is the step size that scales with both the number of training samples n and particles N, βN _N is the inverse temperature that depends on the number of particles, km∼i.i.d.(,d×d) ξ^m_k .i.d. N( 0, I_d× d) denotes i.i.d. Gaussian noise with zero mean and identity covariance matrix. Moreover, Ω(⋅)P_ (·) is the Euclidean projection onto the set Ω . The particles are initially sampled independently from a probability distribution μ0 _0, i.e., 10,…,N0∼i.i.d.μ0 ω_1^0,…, ω_N^0 i.i.d. _0. This stochastic update rule blends deterministic gradient descent (on the energy landscape defined by the Hamiltonian ℋNH_N) with random perturbations, allowing the system to approximate samples from a Gibbs distribution under appropriate conditions. The projection step ensures the dynamics are constrained to a feasible domain, which may encode structural or regularization constraints critical to the optimization problem. In Algorithm B.1, we summarize these steps in a kernel learning algorithm. Using the random feature samples, we construct the randomized feature map ϕD()=(ϕ1∗(),…,ϕD∗()), φ_D( x)= ( _ ω _1( x),…, _ ω _D( x) ), (B.12) and define the corresponding RKHS-based function class as ℱ=f:f(;)=1D⟨,ϕD()⟩,∈Θ⊂ℝD.F= \f:f( x; θ)= 1 D θ, φ_D( x) ,\; θ∈ ^D \. Under this construction, the primal objective in Eq. (A.1) transforms into the following finite-dimensional optimization problem: max∈Θ∑i=1nL(1DTϕD(i),yi)+λ2D‖22. _ θ∈ _i=1^nL ( 1 D θ^T φ_D( x_i),y_i )+ λ2D\| θ\|_2^2. Appendix C Kernelized attention in transformer architecture Transformers hinge on the self-attention operation (Vaswani et al., 2017). Let the (embedded) input sequence be =def(0,1,…,ℓ−1)∈ℓ×dx X =( x_0, x_1,…, x_ -1)∈ × d_x, where 0=defCLS x_0 = x_CLS is the CLSCLS token. Queries, keys, and values are produced by linear maps ()=Q,()=K,()=V, Q( X)= X W_Q, K( X)= X W_K, V( X)= X W_V, with Q,K∈dx×d W_Q, W_K∈ ^d_x× d and V∈dx×dv W_V∈ ^d_x× d_v. Standard (scaled dot-product) attention is ()=softmax(1d()()⊤)(). A( X)=softmax\! ( 1 d\, Q( X) K( X) )\, V( X). Elementwise, ()=(i())i=0ℓ A( X)=( a_i( X))_i=0 , where the i-th output i∈dv a_i∈ ^d_v is a normalized kernel smoother (a.k.a. Nadaraya–Watson estimator) i()=∑j=0ℓK(i,j)j∑j=0ℓK(i,j),K(,)=exp(1d⊤), a_i( X)= _j=0 K( q_i, k_j)\, v_j _j=0 K( q_i, k_j), K( q, k)= \! ( 1 d q k ), (C.1) which connects attention to classical nonparametric regression (Nadaraya, 1964; Watson, 1964) and to kernel interpretations of Transformers (Tsai et al., 2019). Thus, softmax attention is attention with a specific positive kernel K; the denominator enforces a convex combination and stabilizes the estimator. Replacing softmax with any positive kernel K yields a family of attentions with controllable inductive bias (locality, smoothness, anisotropy). When K is approximated via a random feature map ϕN:d→N,↦ϕN=defN−1/2(ϕi())i=1N φ_ _N: ^d→ ^N,\ x φ_ _N =N^-1/2( _ ω_i( x))_i=1^N such that KN(,)=ϕN()⊤ϕN()K_ _N( q, k)= φ_ _N( q) φ_ _N( k), we obtain a linearized form (Katharopoulos et al., 2020a; Choromanski et al., 2021; Peng et al., 2021): ^i() a_i( X) =ϕN(i)⊤(∑j=1nϕN(j)j⊤)ϕN(i)⊤(∑j=1nϕN(j))=ϕN(i)⊤ϕN(i)⊤, = φ_ _N( q_i) \! ( _j=1^n φ_ _N( k_j)\, v_j ) φ_ _N( q_i) \! ( _j=1^n φ_ _N( k_j) )= φ_ _N( q_i) G φ_ _N( q_i) z, (C.2) G =def∑j=0ℓϕN(j)j⊤∈N×dv,=def∑j=0ℓϕN(j)∈N. = _j=0 φ_ _N( k_j)\, v_j ∈ ^N× d_v, z = _j=0 φ_ _N( k_j)∈ ^N. (C.3) This reduces complexity from (ℓ2d)O( ^2d) to (ℓN+Ndv)O( N+Nd_v) per head, with memory (Ndv)O(Nd_v), while preserving the normalized kernel-smoothing structure; causal/padding masks apply by omitting the masked terms in the sums. Alternative sub-quadratic routes include Nyström approximations and kernelized attention with relative positional encoding (Xiong et al., 2021; Chen et al., 2021; Luo et al., 2021). C.1 Alignment for attention kernel. Since labels depend on the entire input sequence through ((),(),()) ( Q( X), K( X), V( X) ), the alignment problem departs from the usual random feature setup in Eq. (B.2). For a one-layer Transformer, collect the attention outputs into ()=def(^0(),…,^ℓ())∈ℓ×dv. ( X) = ( a_0( X),…, a_ ( X) )∈ × d_v. A linear classifier can then be applied either • at token level for token classification problem (e.g., NER): y^i=⊤i()+b y_i= w a_i( X)+b, i=1,2,…,ℓi=1,2,…, , • at sequence level via a pooling operator P:ℓ×dv→dvP: × d_v→ ^d_v: y^=⊤P(())+b,P(())=∑i=0ℓπii(),∈Δℓ. y= w P\! ( ( X) )+b, P\! ( ( X) )= _i=0 _i\, a_i( X),\ \ π∈ . Here Δℓ is the probability simplex, ensuring a convex combination. Examples: [CLS][CLS]: P()=^0()P( )= a_0( X) (i.e., π0=1 _0=1, others πk=0,∀0<k≤ℓ _k=0,∀ 0<k≤ ) where =def(0,1,…,ℓ) X =( x_0, x_1,…, x_ ); mean pooling: πi=1ℓ _i= 1 for all i=1,⋯,ℓi=1,·s, and π0=0 _0=0. (With padding/masks mi∈0,1m_i∈\0,1\, use πi=mi∑j=1nmj _i= m_i _j=1^nm_j for all i=1,⋯,ℓi=1,·s, .) This mirros the classical random-feature setting—where a random map ϕ() φ( x) feeds a linear classifier—but here the attention features () ( X) (and their pooled variant) are functions of the entire sequence X. For sequence-level multi-class classification, given training sequences and labels (i,yi)i=1n\( X_i,y_i)\_i=1^n where the labels are yi∈=1,2,⋯,my_i =\1,2,·s,m\, and fixed embedding matrices (Q,K,V)( W_Q, W_K, W_V), we optimize the target-alignment objective maxμ∈ℳNNseq(N)=1n(n−1)∑1≤i<j≤n(yi⊤yj)P((i))⊤P((j)), _μ _N\ V^seq_N( _N)\;= 1n(n-1) _1≤ i<j≤ n ( e_y_i^\! e_y_j )\,P( ( X_i))^\! P( ( X_j)), (C.4) where P((i))∈ℝdvP( ( X_i)) ^d_v is the pooled representation of sequence i X_i, and kk=1m\ e_k\_k=1^m denotes the standard basis of ℝmR^m which hot-encode the label. Note that the dependence on the particle set N _N (equivalently, on μ) enters only through the sequence embeddings (⋅) (·) induced by the attention features and the pooling operator P. In particular, N=(1,⋯,N) _N=( ω_1,·s, ω_N) is implicit in the definition of attention features () ( X) through each coordinate ^i(),i=0,⋯,ℓ a_i( X),i=0,·s, . Similarly, for token level classification, consider the training dataset (i,1,yi,1),⋯,(i,ℓ,yi,ℓ))i=1n\( x_i,1,y_i,1),·s,( x_i, ,y_i, ))\_i=1^n. Then, the target alignment problem reads maxμ∈ℳNNtok(N)=1||(||−1)∑(i,t)<(j,s)(i,t),(j,s)∈(yi,t⊤yj,s)^i,t(i)⊤^j,s(j). _μ _N\ V^tok_N( _N)= 1|U|(|U|-1) _ subarrayc(i,t)<(j,s)\\ (i,t),(j,s) subarray ( e_y_i,t^\! e_y_j,s )\, a_i,t( X_i)^\! a_j,s( X_j). (C.5) We minimize the energy Nseq(N)V^seq_N( _N) and Ntok(N)V^tok_N( _N) in conjunction with the Coloumb/Riesz regularizer, for sequence level and token level classification, respectively. Remark. In standard Transformer blocks, the attention sublayer is followed by a position-wise feed-forward network (FFN), with each sublayer wrapped by residual connections and layer normalization. A canonical FFN acts independently at each position: FFN(i)=2σ(1i()+1)+2,1∈dv×dff,2∈dff×dv,FFN( a_i)= W_2\,σ( W_1 a_i( X)+ b_1)+ b_2, W_1\!∈\! ^d_v× d_f,\ W_2\!∈\! ^d_f× d_v, with nonlinearity σ (e.g., GELU). Using LNLN for LayerNorm, the residual/normalization update is, schematically, ~i=LayerNorm(i+FFN(i)) a_i=LayerNorm\! ( a_i+FFN( a_i) ). Because the FFN is parameter-shared across positions and applied pointwise, these operations (i) preserve sequence length and token indices; (i) leave the attention weights and the kernel-smoothing form that produced i() a_i( X) unchanged; and (i) implement a learned, per-token reparameterization of attention outputs that improves expressivity and optimization stability. We include FFN, residual, and normalization components in our numerical experiments; the linear classifier discussed above is used to isolate the representation induced by attention and to draw a precise parallel with classical random-feature models. Appendix D Positive random-feature maps for linearized attention We use linearized attention in normalized Nadaraya–Watson form (cf. Katharopoulos et al., 2020b; Vaswani et al., 2017), instantiated by positive feature maps ϕ:ℝd→ℝ+M φ:R^d\!→\!R^M_+ so that the kernel K(,)=defϕ()⊤ϕ()K( q, k) = φ( q)^\! φ( k) is positive semidefinite and the normalization is well-posed. Throughout, we use a particle parameterization with columns 1,…,M ω_1,…, ω_M of ∈ℝd×M \!∈\!R^d× M learned from data (as opposed to fixed i.i.d. draws Choromanski et al., 2020; Peng et al., 2021). Stabilizations used in practice mirror the code: a temperature τ>0τ>0, a small positive floor, and per-token ℓ1 _1 re-normalization when stated. Positive exponential random features (favor). The exponential map ϕfavor()=def1M(exp(i⊤−12‖22))i=1M φ favor_ ( x)\; =\; 1 M ( ( ω_i^\! x- 12\| x\|_2^2) )_i=1^M yields K(,)=1M∑i=1Mexp(i⊤+i⊤−12‖2−12‖2)≥0K( q, k)= 1M _i=1^M ( ω_i^\! q+ ω_i^\! k- 12\| q\|^2- 12\| k\|^2)≥ 0. With i∼(0,) ω_i\! \!N(0,I) and M→∞M\!→\!∞ this recovers a Monte Carlo approximation to the softmax kernel (Choromanski et al., 2020; Peng et al., 2021); we instead learn to obtain a task-adaptive kernel. Deterministic positive features (ELU+1+1 with floor + ℓ1 _1 renorm). Following Katharopoulos et al. (2020b), we use ϕelu+1()=def1Z()(1+ELU((i⊤)/τ+bi))i=1M,Z()∝∑imaxELU(⋅)+1,floor, φ^elu+1_ ( x)\; =\; 1Z( x) (1+ELU\! (( ω_i^\! x)/τ+b_i ) )_i=1^M, Z( x) _i \!\ELU(·)+1,\,floor\, where ELUELU is from Clevert et al. (2015). We clamp to a small floor and re-normalize per token so that ∑iϕi()=M _i _i( x)= M. Softplus features (floor + ℓ1 _1 renorm). ϕsoftplus()=def1Z()(softplus((i⊤)/τ)+floor)i=1M, φ^softplus_ ( x)\; =\; 1Z( x) (softplus\! (( ω_i^\! x)/τ )+floor )_i=1^M, again strictly positive, floor-stabilized, and re-normalized; see Dugas et al. (2001) for softplus. Squared-sigmoid features (ℓ1 _1 renorm). ϕsigmoid2()=def1Z()(σ((i⊤)/τ)2)i=1M,σ(t)=11+e−t, φ^sigmoid2_ ( x)\; =\; 1Z( x) (σ\! (( ω_i^\! x)/τ )^2 )_i=1^M, σ(t)= 11+e^-t, which are non-negative and re-normalized per token. Softmax-over-features (per-token). We also consider a feature-softmax map ϕsoftmaxfeat()=defMsoftmax((⊤)/τ), φ^softmaxfeat_ ( x)\; =\; M\,softmax\! (( ^\! x)/τ ), which is strictly positive and sums to M M by construction. Cosine-squared random features. Motivated by random Fourier features for shift-invariant kernels (Rahimi and Recht, 2007, 2009), we use ϕ,cos2()=def1Z()(cos(i⊤+bi)2)i=1M, φ^cos2_ , b( x)\; =\; 1Z( x) ( ( ω_i^\! x+b_i)^2 )_i=1^M, with fixed phases b; non-negativity is immediate, and we apply per-token ℓ1 _1 renormalization. PORF-Softplus (orthogonal initialization). To reduce variance and improve conditioning, we initialize with orthogonal random features blocks (Yu et al., 2016; Choromanski et al., 2017) and then apply the softplus map above: ϕporf-softplus()≡ϕsoftplus(),⊤≈. φ^porf -softplus_ ( x)\;≡\; φ^softplus_ ( x), . The parameters remain learnable after orthogonal initialization. Learning the kernel via alignment/KTA (Phase A). Beyond a vanilla cross-entropy training of the classifier head (Phase B), we first adapt with a representation-level objective: either a within-class alignment loss (maximizing same-class similarity) or Kernel Target Alignment (KTA; Cristianini et al., 2002) using centered Gram matrices K and label kernel Y. This yields task-adaptive positive kernels while preserving linear-time forward/backward passes. We evaluate on SST-2 from GLUE (Wang et al., 2019) with BERT tokenization (Devlin et al., 2019), matching the experimental setup in our code. Implementation details. All maps use a temperature τ, small positive floors where applicable, and (except favor and feature-softmax) per-token ℓ1 _1 re-normalization to keep ∑iϕi()=M _i _i( x)= M. Columns of are optionally constrained by column-wise ℓ2 _2 clipping during Phase A to stabilize learning. Appendix E Theoretical results Before presenting our theoretical results, we first outline the key assumptions that underpin our analysis: Assumption A.1 (Initial distribution of particles). The particles are initially sampled independently from a probability distribution μ0 _0 that admits a Lebesgue density ρ0()=μ0(d)d _0( ω)= _0(d )d ω. Assumption A.2 (Constant Temperature). The temperature parameter βN _N remains constant and finite throughout the dynamics for a fixed number of particles N. Assumption A.3 (Projection space). The feature domain Ω⊂d ⊂ ^d is compact and convex with C2C^2 boundary. The projected dynamics use non-absorbing, reflective boundary conditions on ∂Ω∂ . Assumption A.4 (Bounded Random Feature Embedding). Let Φ()∈L2(Ω,μ0) ( x)∈ L^2( , _0) denote the random feature embedding associated with the kernel K(,′)K( x, x ). We assume that supx∈‖Φ(x)‖L2(Ω,μ0)=supx∈K(x,x)≤LΦ<∞. _x \| (x)\|_L^2( , _0)= _x K(x,x)≤ L_ <∞. (E.1) E.1 Continuity equations and equilibrium distribution Theorem E.1 (Projected-particle mean-field continuity equation). Let Ω⊂d ⊂ ^d be compact and convex with C2C^2 boundary and outward unit normal n. Assume that V∈C2(Ω¯)V∈ C^2( ), that gsg_s is smoothly regularized or the particle system remains collision-free, and that μN0⇒ρ0()d _N^0 _0( ω)d ω. Consider the mean-field-scaled projected Langevin iteration km+1=Ω¯(km−ηN[∇V(km)+λN−1∑ℓ≠k∇gs(km−ℓm)]+2ηN/βkm), ω_k^m+1=P_ \! ( ω_k^m- _N [∇ V( ω_k^m)+ λN-1 _ ≠ k∇ g_s( ω_k^m- ω_ ^m) ]+ 2 _N/β\, ξ_k^m ), (E.2) where km∼i.i.d.(,Id) ξ_k^m .i.d. N( 0,I_d), ηN↓0 _N 0, and the projected-Euler consistency error vanishes as N→∞N→∞. Let μtN=N−1∑k=1Nδk⌊t/ηN⌋ _t^N=N^-1 _k=1^N _ ω_k t/ _N . Then μtN⇒μt=ρt()d _t^N _t= _t( ω)d ω in probability, uniformly for t in compact intervals. The limiting density is governed by the McKean–Vlasov equation ∂ρt()∂t=∇⋅[ρt()∇Ut()]+1βΔρt(),(t,)∈(0,T]×Ω, ∂ _t( ω)∂ t= _ ω· [ _t( ω) _ ωU_t( ω) ]+ 1β _ ω _t( ω), (t, ω)∈(0,T]× , (E.3) where Ut()=V()+λ∫Ωgs(−′)ρt(′)′. U_t( ω)=V( ω)+λ _ g_s( ω- ω ) _t( ω )\,d ω . (E.4) The PDE is supplemented by the initial datum, the Robin/no-flux boundary condition, and conservation of probability: ρt|t=0=ρ0,(ρt∇Ut+β−1∇ρt)⋅=0on (0,T]×∂Ω,ρt≥0,∫Ωρt()=1. _t|_t=0= _0, ( _t _ ωU_t+β^-1 _ ω _t )· n=0 (0,T]×∂ , _t≥ 0, _ _t( ω)d ω=1. (E.5) Equivalently, for smooth solutions, ∂ρt+βρt∂Ut=0 _ n _t+β _t _ nU_t=0 on ∂Ω∂ . Here V()=−[yy′ϕ()ϕ(′)]V( ω)=-E[y _ ω( x) _ ω( x )] is the external potential induced by the kernel-target alignment loss. The boundary condition in (E.5) is the multidimensional analogue of the Robin condition: it says that the probability flux through the boundary is zero. In the zero-temperature limit (β→∞)(β→∞), the evolution reduces to the reflected deterministic continuity equation ∂ρt()∂t+∇⋅(ρt()t())=0,(ρtt)⋅=0on ∂Ω, ∂ _t( ω)∂ t+ _ ω· ( _t( ω) v_t( ω) )=0, ( _t v_t)· n=0 ∂ , (E.6) where t()=−∇Ut() v_t( ω)=- _ ωU_t( ω). This PDE formalism captures the macroscopic evolution of particle densities driven by repulsion, external alignment forces, reflection, and thermal fluctuations. At equilibrium, for any fixed βN>0 _N>0 and N>0N>0, the particle configuration is distributed according to the canonical Gibbs measure: ℙN,βN(dN)=1ZN,βNexp(−βNℋN(N))ΩN(N)dN, _N, _N(d _N)= 1Z_N, _N (- _NH_N( _N) ) 1_ ^N( _N)\,d _N, (E.7) where ΩN(N) 1_ ^N( _N) is the indicator function of ΩN ^N, dNd _N is the Lebesgue measure on (ℝd)N(R^d)^N, and ZN,βZ_N,β is the partition function, ZN,βN=def∫ΩNexp(−βNℋN(N))dN, Z_N, _N = _ ^N (- _NH_N( _N) )\,d _N, (E.8) which ensures that the Gibbs measure is properly normalized, i.e., ∫ΩNℙN,β(N)dN=1 _ ^NP_N,β( _N)\,d _N=1. Moreover, the free energy is defined as Fs(βN,N,Ω)=def−1βNlogZN,βN. F_s( _N,N, ) =- 1 _N Z_N, _N. (E.9) As βN→+∞ _N→+∞ with N→∞N→∞, the Gibbs measure increasingly concentrates around the minimizer (ground state) of the Hamiltonian ℋN(N)H_N( _N), which corresponds to the solution of optimization problem in Eq. (B.3). E.2 Thermodynamic limit of the equilibrated state in the short-range case, s>ds>d Since the kernel function is approximated by the empirical measure μN _N through the Monte Carlo approximation in Eq. (B), it is essential to analyze the asymptotic behavior of this empirical measure in order to characterize the limiting properties of the kernel approximation itself. In particular, the statistical fluctuations and concentration properties of μN _N directly determine the accuracy and stability of the resulting kernel-based quantities. In the setting of interest in Algorithm B.1, the empirical measure μN _N arises from a system of interacting charged particles evolving under Langevin dynamics. Once the dynamics have reached equilibrium (i.e., after a sufficiently large number m of iterations), the distribution of the particle system converges to the canonical Gibbs measure given in Eq. (E.7). This measure describes the statistical equilibrium of the system, incorporating both the deterministic interaction potential and the stochastic perturbations induced by thermal noise. From the perspective of statistical mechanics, such an equilibrium corresponds to a thermodynamically stable macroscopic state, in which relevant observables become stationary in distribution. The principal objective of this section is to investigate the asymptotic behavior of the random empirical measure μN _N in this equilibrium regime, particularly as the number of particles N tends to infinity. This regime, known as the thermodynamic limit (N→∞N→∞), is of fundamental importance in both statistical mechanics and probability theory, as it establishes the connection between microscopic particle interactions and macroscopic statistical laws. Our main result is the derivation of a large deviation principle (LDP) for μN _N, which characterizes the exponential decay of probabilities of rare deviations from the typical equilibrium distribution. The LDP is governed by a rate function, which assigns to each admissible probability measure a nonnegative “cost” quantifying the likelihood of its occurrence in the large-N limit. This framework provides a precise quantitative description of the concentration of μN _N around its equilibrium value, as well as the nature of its fluctuations. For a comprehensive background on the theory of large deviations, we refer the reader to Dembo and Zeitouni (2009) and Varadhan (2016). For completeness, we recall the formal definition below, which introduces the notion of an LDP and the associated rate function. Definition E.2 (Large Deviation Principle (LDP)). Let (νN)N≥1( _N)_N≥ 1 be a sequence of probability measures on a Polish space Ω equipped with the Borel σ-algebra ℬ(Ω)B( ). We say that (νN)( _N) satisfies a Large Deviation Principle (LDP) at speed rNr_N with rate function ℐ:Ω→ℝ+I: _+ if, for every Borel set B⊂ℬ(Ω)B ( ), −infx∈B̊ℐ(x)≤lim infN→∞1rNlogνN(B)≤lim supN→∞1rNlogνN(B)≤−infx∈B¯ℐ(x),- _x∈ BI(x)≤ _N→∞ 1r_N _N(B)≤ _N→∞ 1r_N _N(B)≤- _x∈ BI(x), where B̊ B and B¯ B denote the interior and closure of B, respectively. The functional ℐI is called a good rate function if it is lower semi-continuous and has compact sub-level sets. Equipped with Definition E.2, we are ready to state the following theorem: Theorem E.3 (Large deviations for empirical measures). Let Ω⊂ℝd ^d be a bounded open set. For each N≥1N≥ 1, let 1,…,N⊂Ω\ ω_1,…, ω_N\⊂ be a random configuration with law ℙN,βNP_N, _N as in (E.7), and define the empirical measure μN=1N∑k=1Nδk∈(Ω). _N\;=\; 1N _k=1^N _ ω_k\;∈\;P( ). Then μNN≥1\ _N\_N≥ 1 satisfies a large deviation principle on (Ω)P( ) endowed with the weak topology, with speed rNr_N and good rate function given as follows. Define the energy functional ℰs(μ)=def∫Ωϕ(i)ϕ(j)μ(d)+λ2∫Ω∫Ωgs(−′)μ(d)μ(d′),E_s(μ) = _ _ ω( x_i)\, _ ω( x_j)\,μ(d ω)\;+\; λ2 _ _ g_s( ω- ω )\,μ(d ω)\,μ(d ω ), and the (relative) entropy with respect to Lebesgue measure ℓ on Ω , Ent(μ∣ℓ)=def∫Ωlog(dμdℓ())μ(d)=∫Ωdμdℓ()log(dμdℓ())ℓ(d),μ≪ℓ,+∞,otherwise.Ent(μ ) = cases _ \! ( dμd ( ω) )\,μ(d ω)= _ dμd ( ω) \! ( dμd ( ω) )\, (d ω),&μ ,\\[8.61108pt] +∞,&otherwise. cases • If βN/N→1 _N/N→ 1 as N→∞N→∞, then the LDP holds with speed rN=Nr_N=N and rate function s(μ)=(ℰs(μ)+Ent(μ∣ℓ))−infν∈(Ω)(ℰs(ν)+Ent(ν∣ℓ)).J_s(μ)= (E_s(μ)+Ent(μ ) )- _ν ( ) (E_s(ν)+Ent(ν ) ). • If βN/N→∞ _N/N→∞ as N→∞N→∞, then the LDP holds with speed rN=βNr_N= _N and rate function s(μ)=ℰs(μ)−infν∈(Ω)ℰs(ν).J_s(μ)=E_s(μ)- _ν ( )E_s(ν). The proof is provided in Appendix G.4. Theorem E.3 implies that for any Borel set A⊂(Ω)A ( ), −infμ∈A∘s(μ)≤ - _μ∈ A J_s(μ)\;≤\; lim infN→∞1rNlogℙN,βN(μN∈A) _N→∞ 1r_N _N, _N( _N∈ A) ≤lim supN→∞1rNlogℙN,βN(μN∈A)≤−infμ∈A¯s(μ), \;≤\; _N→∞ 1r_N _N, _N( _N∈ A)\;≤\;- _μ∈ AJ_s(μ), where A∘A and A¯ A denote the interior and closure of A, respectively. The theorem states that, under the Gibbs law, the empirical measure μN _N concentrates (in the weak topology) around the minimizers of the relevant variational functional, and the probability of observing a macroscopic deviation decays exponentially fast at speed rNr_N. In the regime βN∼N _N N, the rate function contains both the interaction energy ℰs(μ)E_s(μ) and the entropy term Ent(μ∣ℓ)Ent(μ ), capturing the competition between energetic preference for structured configurations and entropic preference for spreading mass. In the low–temperature regime βN/N→∞ _N/N→∞, the entropy contribution becomes negligible, so the large deviations are governed purely by ℰsE_s and μN _N concentrates on energy minimizers, with fluctuations suppressed on the faster exponential scale set by βN _N. Figure 1: Accuracy (left) and F1 score (right) versus noise standard deviation σ2σ^2 for kernel-learning approaches. Confidence intervals summarize variability across 10 independent trials. Table 5: Observed wall-clock runtime per job in seconds. Values correspond only to measured sample sizes. Model Number of samples n 100100 500500 1,0001,000 2,0002,000 3,0003,000 4,0004,000 10,00010,000 20,00020,000 Linear Logistic 0.142 0.100 0.0882 0.184 0.173 0.204 0.117 0.121 SVM (RBF kernel) 0.00231 0.0425 0.140 0.0910 0.317 0.468 1.98 9.72 RFF (linear SVM) 0.00332 0.0250 0.0270 0.0326 0.0416 0.0571 0.160 0.335 MKL-SVM 2.06 9.14 29.6 131 482 874 7.21e3 4.49e4 Optimized RFF (lin SVM) 0.0346 0.127 0.299 0.533 0.855 1.04 2.94 4.87 Langevin (Alg-1, lin SVM) 21.6 26.1 45.8 80.5 95.9 124 417 922 Table 6: Power-law fits t(n)=cnγt(n)=c\,n^γ summarizing runtime scaling behavior. Exponent γ measures the growth rate, prefactor c sets the baseline runtime, and σlog _ indicates variability in log-space residuals. Model γ c σlog _ nobsn_obs MKL-SVM 1.961 8.00e-05 0.671 8 SVM (RBF kernel) 1.477 3.00e-06 0.503 8 Optimized RFF (lin SVM) 0.967 3.55e-04 0.085 8 RFF (linear SVM) 0.799 9.50e-05 0.319 8 Langevin (Alg-1, lin SVM) 0.730 3.82e-01 0.430 8 Linear Logistic 0.021 1.16e-01 0.290 8 Figure 2: Ablation study of Coulomb (top) and Riesz (bottom) potentials for charged particles. Increasing λ strengthens the repulsive interaction between particles. The baseline uses vanilla RFF with no sampling optimizations. Confidence intervals summarize variability across 10 independent trials. Appendix F Numerical Experiments We present simulations on synthetic datasets as well as experiments on NLP datasets. We plan to release the code for all experiments publicly upon publication. F.1 Nonlinear synthetic classification: data, estimators, and evaluation protocol F.1.1 Problem setting. We consider binary classification with a fixed design matrix ∈ℝn×p X ^n× p drawn once as Xij∼i.i.d.(0,1)X_ij i.i.d. N(0,1) under a global seed. Unless stated otherwise we use n=400n=400 and a configurable feature dimension p≥5p≥ 5. Writing i x_i for the i-th row of X, the log-odds is ℓ(i)=1.5sin(πxi,1)+ 0.8xi,2 2− 1.0xi,3xi,4+ 0.5sin(3xi,5)+ 1p>5i,6:p⊤, ( x_i)=1.5\, (π\,x_i,1)\;+\;0.8\,x_i,2^\,2\;-\;1.0\,x_i,3\,x_i,4\;+\;0.5\, (3\,x_i,5)\;+\;1_\p>5\\; x_i,6:p^\! w, where the “extra–dimensions” weight vector ∈p−5 w∈ ^p-5 is drawn once from (,0.32) N( 0,0.3^2 I) (same w for all trials), and ⋅ 1_\·\ is the indicator function. Labels are generated by corrupting the logit with additive Gaussian noise and passing through a logistic link: ξi∼(0,σ2),yi|i∼Bernoulli(ς(ℓ(i)+ξi)),ς(t)=11+e−t. _i N(0,σ^2), y_i\, |\, x_i \! ( \! ( ( x_i)+ _i ) ), (t)= 11+e^-t. We sweep the logit noise standard deviation σ over logarithmically spaced values in [10−5,105][10^-5,10^5]. A single 70%70\%/30%30\% train/test split (random state 0) of X is reused across all trials and noise levels; for each σ and trial, only the noise (ξi)( _i) (and any method-specific randomness) is resampled. F.1.2 Baselines and learned feature maps. All methods use the same train/test partitions of X. Random-feature methods are aligned to a common feature budget D=200D=200. (1) Linear logistic regression on the raw inputs X (baseline). (2) SVM with Gaussian (RBF) kernel K(,′)=exp(−γ∥−′∥22)K( x, x )= \! (-γ x- x _2^2 ) with fixed γ=0.5γ=0.5 and default hinge loss, C=1C=1. (3) Random Fourier features (RFF) + linear SVM. Draw W:,j∼(0, 2γp)W_:,j\! \! N(0,\,2γ\, I_p), bj∼Unif[0,2π]b_j\! \!Unif[0,2π], form z()=2/Dcos(W⊤+b)z( x)= 2/D\, (W x+b) with optimized γ=0.5γ=0.5 and D=200D=200, then train a linear SVM (hinge loss, C=1C=1). (4) Multiple-kernel learning (MKL) + SVM (precomputed). Build Gaussian Gram matrices with γ∈0.1,0.3,0.7,1.0,1.3,1.6γ∈\0.1,0.3,0.7,1.0,1.3,1.6\, apply the MKLpy normalization, learn nonnegative mixture weights via MEMO, and feed the weighted train/test kernels to an SVM with a precomputed kernel. (5) Optimized random features (importance sampling) + linear SVM. Following Sinha and Duchi (2016), sample Nw=104N_w=10^4 candidate features, use divergence threshold ρ=Nw⋅0.005ρ=N_w· 0.005 and tolerance 10−1010^-10 to obtain a reweighted set, then align to D=200D=200 via stratified subsampling (padding if needed). Train a linear SVM on the resulting features (C=1C=1). (6) Langevin spectral estimators (two variants) + linear SVM. Evaluate (i) an “Alg-1” Coulomb-gas variant with λreg=0.5 _reg=0.5 and (i) a robust variant with λreg=0 _reg=0. Both use N=300N=300 particles, feature budget D=200D=200, step size η=30η=30, inverse temperature β=102β=10^2, horizon Tmax=2000T_ =2000, logarithmic repulsion, scaling by N, maximum frequency norm 5.05.0, and gradient clipping 1.01.0. After fitting, we form Fourier features from the learned frequencies/phases and train a linear SVM (C=1C=1). F.1.3 Implementation. We rely on MKLpy Lauriola and Aiolli (2020) for MKL (mixture-of-RBFs) and on scikit-learn for linear/logistic models and SVMs. The importance-sampling optimized RFF pipeline of Sinha and Duchi (2016) is reimplemented in Python. The Langevin estimators are also implemented in Python. We run the simulations on CPU. F.1.4 Training protocol, metrics, and plots. For each noise level σ and each of three independent trials, we fit every available method on the training split and evaluate on the test split. We report, for each method and σ, the mean and standard deviation (over trials) of: • Test accuracy: Acc=1ntest∑i∈testy^i=yiAcc= 1n_test _i _test1\ y_i=y_i\. • Test F1-score: (binary, with zero_division=0zero\_division=0 in implementation). • Wall-clock time per job: (seconds), measured with a timer from the start of a method’s pipeline (feature construction included for RFF variants) to prediction on the test set. We produce three summary plots versus σ (logarithmic x-axis): (A) accuracy, (B) F1-score, and (C) runtime; shaded bands depict ±1± 1 standard deviation. F.1.5 Reproducibility and computation. A global seed (4242) fixes X, the train/test split, and (when p>5p>5) the extra-dimension weights w; per-job seeds are derived deterministically from the noise and trial indices. We cap BLAS threads to one in parent and worker contexts to avoid oversubscription. Jobs indexed by (σ,trial)(σ,trial) are executed in parallel using a process pool with up to eight workers; progress is tracked asynchronously. Across methods using random features, the feature budget is held fixed at D=200D=200 to ensure comparable capacity. F.1.6 Results. In Figure 1, we report accuracy and F1 score as functions of the additive-noise standard deviation σ. Across noise levels, the SVM with a Gaussian (RBF) kernel—when its bandwidth γ is properly tuned—achieves the best performance. This is expected because the random-feature model only approximates the kernel; with a small number of features (D=200D=200), approximation error degrades classification performance. The same figure also shows that the RBF SVM can underperform substantially when the bandwidth γ is mis-specified, e.g., γ=0.01γ=0.01 and γ=100γ=100 corresponding to green and red dashed lines, respectively. Moreover, relative to a linear SVM trained on fixed random Fourier features (RFF), jointly optimizing the random features yields clear gains in both accuracy and F1 (purple dashed vs. grey solid curves), confirming the benefit of random-feature optimization. Moreover, compared to alternative importance sampling of random features, we obtain a clear gain by optimization of random feature samples directly (brown line vs purple line). F.1.7 Ablation: Coulomb/Riesz interaction strength. We vary the potential strength λ∈0.0,0.2,0.5,0.8,1.0,1.5λ∈\0.0,0.2,0.5,0.8,1.0,1.5\ and evaluate robustness as input noise increases (x–axis: noise standard deviation on a log scale). Across both metrics (Accuracy and F1), small noise regimes (σ≤10−3σ\!≤\!10^-3) show minimal separation between curves, indicating that the potential has little effect when the signal is clean. As noise grows to the moderate range (100−10210^0\!\!-\!10^2), nonzero potentials consistently outperform the vanilla RFF baseline and λ=0λ=0, with λ≈0.5−1.0λ≈ 0.5\!-\!1.0 yielding the most reliable gains (typically a few points) and smaller variance, suggesting improved stability. In the extreme-noise regime (σ≥103σ\!≥\!10^3), performance degrades for all settings and the gaps narrow; very large strengths (λ=1.5λ\!=\!1.5) can oversmooth and underperform, while small-to-moderate λ remains competitive but offers diminishing returns. Overall, a moderate Coulomb/Riesz potential (λ≈0.5−1.0λ≈ 0.5\!-\!1.0) provides the best robustness–accuracy trade-off, whereas too weak or too strong potentials are less effective. F.2 Sentence-Level Classification Benchmarks Table 7: Statistics of sentence-level and sentence-pair classification datasets used in our experiments. Counts denote the number of examples per split. For sentence-pair tasks, each example consists of a pair of sentences. Dataset Train Dev Test Classes Task SST-2 Socher et al. (2013) 67,349 872 1,821 2 Sentiment classification QQP Wang et al. (2019) 363,846 40,430 390,965 2 Duplicate question detection Rotten Tomatoes Pang and Lee (2005) 8,530 1,066 1,066 2 Sentiment classification Dataset selection. We evaluate sentence-level and sentence-pair classification performance on a set of widely used English benchmarks that are standard in prior work on attention mechanisms. The selected datasets span multiple semantic phenomena and dataset scales, including binary sentiment analysis on SST-2 and Rotten Tomatoes, as well as semantic equivalence detection on the QQP benchmark. Together, these tasks cover both single-sentence and sentence-pair settings, enabling evaluation across varying supervision regimes and levels of linguistic complexity. Datasets and splits. SST-2 Socher et al. (2013) is a binary sentiment classification task derived from the Stanford Sentiment Treebank. We use the GLUE version of the dataset, which removes neutral examples and provides sentence-level annotations, and follow the standard GLUE train/dev/test splits. QQP Wang et al. (2019) (Quora Question Pairs) is a sentence-pair classification task in which each example consists of two questions from Quora annotated for semantic equivalence. We adopt the official GLUE splits and process each question pair by concatenation with a special separator token. Rotten Tomatoes Pang and Lee (2005) is a sentence-level binary sentiment classification dataset constructed from movie reviews. We use the standard polarity version of the dataset and follow the commonly used train/dev/test splits. Preprocessing. We apply minimal preprocessing to preserve the original linguistic structure of each dataset. All text is tokenized using a fixed WordPiece vocabulary, with original casing and punctuation retained. Inputs are truncated or padded to a maximum length of 128 tokens. For sentence-pair inputs in QQP, the two questions are concatenated using a special separator token. No segment embeddings or task-specific features are used unless otherwise stated. Training protocol. For each dataset, we train models on the official training split and use the validation split exclusively for model selection and early stopping. Hyperparameters are fixed across all experiments and are not tuned on validation data. For each method, we select the checkpoint that achieves the highest validation accuracy and report performance on the held-out test split using this checkpoint. Test data are never used for hyperparameter selection or model selection. Evaluation and metrics. Across all datasets, we report classification accuracy as the primary performance metric. For binary classification tasks (SST-2, QQP, and Rotten Tomatoes), we additionally report micro-, macro-, and weighted F1F_1 scores to account for potential class imbalance. To evaluate ranking quality and probabilistic reliability, we also report ROC-AUC, precision–recall AUC, Matthews correlation coefficient (MCC), balanced accuracy, LogLoss, Brier score, and expected calibration error (ECE). All metrics are computed on the test split using the checkpoint selected based on validation performance. Experimental setup. Inputs are tokenized using the bert-base-uncased WordPiece tokenizer and truncated or padded to a maximum sequence length of 128 tokens. Models use a two-layer Transformer encoder with hidden size 128, two attention heads, feedforward dimension 256, and m=256m=256 random features per head. In our particle-based view, each column of N∈ℝH×dk×m _N ^H× d_k× m is a particle, so this corresponds to 256 particles per head (512 per layer; 1024 particles total across the two-layer encoder). Dropout is set to 0.1, and mean pooling is applied over token representations. We use a log repulsion force between particles. Training uses a batch size of 64. When kernel optimization is enabled, training proceeds in two phases. In Phase A, the attention feature map parameters N _N (and optionally layer normalization and value projection parameters) are optimized using an alignment loss objective for up to 10 epochs via particle optimization (SGLD/Langevin dynamics with step size η=2×10−3η=2× 10^-3 and inverse temperature β=50β=50, i.e., Gaussian noise scale 2η/β 2η/β, plus a repulsion term with λ=10−3λ=10^-3 and gradient clipping at norm 10), while the remaining model parameters are frozen. Column-wise L2L_2 norms of N _N (i.e., per-particle norms) are constrained to 1.5. Training in Phase A stops early if the maximum change in N _N falls below 10−610^-6. In Phase B, the kernel parameters are frozen and the model is trained end-to-end with cross-entropy loss for 10 epochs. For vanilla baselines, only Phase B is used. Optimization in Phase B is performed with Adam using a learning rate of 2×10−42× 10^-4. Effect of Kernel Learning. Tables 1, 2, and 3 quantify the effect of target-alignment kernel learning across SST-2, QQP, and Rotten Tomatoes. On SST-2 (Table 1), where overall performance is near saturation, kernelization yields feature-map-dependent gains. Kernel-porf-softplus attains the best accuracy and MCC (0.8188/0.6376). Several kernel variants improve both discriminative and probabilistic metrics: for example, Kernel-elu increases accuracy from 0.8028 to 0.8108 while simultaneously reducing LogLoss (0.4938→ 0.4502) and Brier score (0.1505→ 0.1392), and Kernel-softmaxfeat improves accuracy (0.8005→ 0.8085) while achieving near-optimal LogLoss (0.4381). Calibration-specific improvements are also evident: Kernel-cos2 yields the lowest Brier score (0.1366), while Kernel-favor achieves the lowest ECE (0.3086). However, not all kernels are uniformly beneficial—Kernel-softplus degrades both Brier and ECE—highlighting sensitivity to the choice of feature map. On the more challenging QQP task (Table 2), kernel learning provides more consistent benefits across feature maps. Kernel-softmaxfeat is the strongest Q=KQ=K model (Acc 0.7909, MCC 0.5506), with substantial reductions in LogLoss (0.4625→ 0.4354) and Brier score (0.1512→ 0.1418), while Kernel-favor achieves the best calibration as measured by ECE (0.4106). Although kernelization roughly doubles training time, it does not affect inference-time complexity, indicating that alignment primarily reshapes representation geometry and probabilistic reliability rather than merely sharpening decision boundaries. Results on Rotten Tomatoes (Table 3) exhibit similar but more modest trends. Kernel-softmaxfeat achieves the highest accuracy and MCC (0.7167/0.4345) while also attaining the lowest LogLoss (0.5697) and Brier score (0.1929), suggesting that kernel alignment is particularly effective in improving probabilistic calibration on smaller, noisier datasets. As on SST-2, gains vary across feature maps, reinforcing the importance of the kernel choice. Despite these improvements, Vaswani-style softmax attention remains the strongest performer on QQP. Its advantage stems from explicitly modeling dense all-pairs token interactions through an ℓ×ℓ × attention matrix, incurring (ℓ2)O( ^2) time and memory complexity in the sequence length ℓ . This expressivity is well suited to paraphrase detection, where fine-grained cross-token alignment is critical, but the quadratic scaling limits practicality for longer sequences. In contrast, linear-attention variants trade some modeling capacity for (ℓ)O( )-type scaling (up to the feature dimension), making them preferable when sequence length or throughput constraints dominate. Appendix G Appendix We organize the appendix as follows: • In Section G.3, we present the proof of Theorem E.1. • In Section G.4 we present the large deviation result of Theorem E.3. G.1 Notation We define the notation as follows: • ℳ(Ω)M( ): The space of measures on the measurable space (Ω,ℬ)( ,B), where ℬB is the Borel σ-algebra on Ω . • μ^∈ℳ(Ω) μ ( ): A random (counting) measure, typically representing the realization of the point process. It is a random element in the space of measures on Ω . • μ∈ℳ(Ω)μ ( ): The deterministic (limiting) measure, often representing the stationary or equilibrium distribution of the point process as N→∞N→∞. • :Ξ→ℳ(Ω)P: ( ): A point process mapping from the probability space (Ξ,ℱ,ℙ)( ,F,P) to the space of measures on Ω . • ν=ℙ∘−1ν=P ^-1: The push-forward measure of the point process, governing the randomness of the point process and describing the distribution of realizations of the point process. • Λ : The intensity measure, defined as Λ(B)=μ∼ν[μ(B)]forB∈ℬ, (B)=E_μ ν[μ(B)] B , where μ is the random measure associated with a realization of the point process. • Ξ : The underlying probability space, often taken as a sample space of configurations for the point process. • μ^(B) μ(B): The counting measure or number of points in a subset B∈ℬB , corresponding to the realization of the random measure μ μ. • [⋅]E[·]: Expectation with respect to the probability measure governing the point process. In this context, it is typically the expectation under ν, the push-forward measure. • μ(B)μ(B): The number of points in the set B as determined by the limiting measure μ. • X: A random variable (or random element) associated with the point process, representing a random realization or observation. • x: A realization of the random variable X, i.e., an outcome or observation from the random process. G.2 Proof of Lemma B.1 Insert ϕk,bk()=2cos(k⊤+bk) _ ω_k,b_k( x)= 2 ( ω_k x+b_k) into (B.4): 1N∑k=1Nϕk,bk(i)ϕk,bk(j)=2N∑k=1Ncos(αik)cos(αjk),αik=defk⊤i+bk. 1N _k=1^N _ ω_k,b_k( x_i) _ ω_k,b_k( x_j)= 2N _k=1^N ( _ik) ( _jk), _ik = ω_k x_i+b_k. Substituting and swapping finite sums gives ℰN(N,)=−2n(n−1)N∑k=1N∑i≠jyiyjcos(αik)cos(αjk). _N( _N, b)=- 2n(n-1)N _k=1^N _i≠ jy_iy_j ( _ik) ( _jk). (G.1) For any ,∈ℝn a, y ^n, ∑i≠jyiyjaiaj=(∑iyiai)2−∑iyi2ai2. _i≠ jy_iy_ja_ia_j= ( _iy_ia_i )^2- _iy_i^2a_i^2. (G.2) Applying (G.2) with ai=cos(αik)a_i= ( _ik) yields ∑i≠jyiyjcos(αik)cos(αjk)=(⊤k)2−∑iyi2cos2(αik),(k)i=defcos(αik). _i≠ jy_iy_j ( _ik) ( _jk)= ( y c_k )^2- _iy_i^2 ^2( _ik), ( c_k)_i = ( _ik). Take expectation over the phases b, assumed i.i.d. Unif[0,2π]Unif[0,2π] and independent of (N,i)( _N,\ x_i\): bk[cos2(k⊤i+bk)]=12.E_b_k [ ^2( ω_k x_i+b_k) ]\;=\; 12. Therefore, bk[∑iyi2cos2(αik)]=12∑iyi2.E_b_k\! [ _iy_i^2 ^2( _ik) ]\;=\; 12 _iy_i^2. In binary classification (yi2=1y_i^2=1), this equals n/2n/2 and is independent of k ω_k. Next, write uik=defk⊤iu_ik = ω_k x_i and expand cos(uik+bk)=cosuikcosbk−sinuiksinbk. (u_ik+b_k)= u_ik b_k- u_ik b_k. Define the vectors k=(cosu1k,…,cosunk)⊤ u_k=( u_1k,…, u_nk) , and k=(sinu1k,…,sinunk)⊤. v_k=( u_1k,…, u_nk) .. Then ⊤k=⊤(cosbkk−sinbkk) y c_k= y \! ( b_k\, u_k- b_k\, v_k ) (G.3) and (⊤k)2=(cosbk)2(⊤k)2+(sinbk)2(⊤k)2−2cosbksinbk(⊤k)(⊤k). ( y c_k)^2=( b_k)^2( y u_k)^2+( b_k)^2( y v_k)^2-2 b_k b_k\,( y u_k)( y v_k). (G.4) Taking expectation over bkb_k and using [cos2bk]=[sin2bk]=12E[ ^2b_k]=E[ ^2b_k]= 12 and [sinbkcosbk]=0E[ b_k b_k]=0, we obtain bk[(⊤k)2]=12((⊤k)2+(⊤k)2).E_b_k\! [( y c_k)^2 ]= 12 (( y u_k)^2+( y v_k)^2 ). Combining the two expectations, bk[∑i≠jyiyjcos(αik)cos(αjk)]=12((⊤k)2+(⊤k)2)−12∑iyi2,E_b_k\! [ _i≠ jy_iy_j ( _ik) ( _jk) ]= 12 (( y u_k)^2+( y v_k)^2 )- 12 _iy_i^2, where the second term is a constant (equal to n/2n/2 in the binary case) and thus does not depend on k ω_k. Summing over k and inserting into (G.1), [ℰN(N,)]=−1n(n−1)N∑k=1N((⊤k)2+(⊤k)2),E_ b\! [E_N( _N, b) ]=- 1n(n-1)N _k=1^N (( y u_k)^2+( y v_k)^2 ), where the equality is up to an additive constant independent of N _N. Finally, stack columns =cos(N⊤)=[1⋯N],=sin(N⊤)=[1⋯N] C= ( X _N^\! )=[ u_1\,·s\, u_N],\; S= ( X _N^\! )=[ v_1\,·s\, v_N] to write ∑k=1N((⊤k)2+(⊤k)2)=‖⊤‖22+‖⊤‖22. _k=1^N (( y u_k)^2+( y v_k)^2 )=\| y C\|_2^2+\| y S\|_2^2. This yield the stated result [ℰN(N,)]≡−1n2(‖⊤cos(N⊤)‖22+‖⊤sin(N⊤)‖22).E_ b\! [E_N( _N, b) ]\;≡\;- 1n^2 ( \| y ( X _N^\! ) \|_2^2+ \| y ( X _N^\! ) \|_2^2 ). G.3 Proof of Theorems 4.1 and E.1 Proof G.1. We prove Theorem E.1; Theorem 4.1 follows as the special case V=V=V_D. We present a proof sketch based on the projected-particle mean-field approach. Specifically, we construct a continuous-time interpolation of the projected Langevin particles and compare it with a reflected Itô diffusion. A mean-field (propagation-of-chaos) estimate is then used to replace the empirical drift with the law-dependent drift. Finally, we identify the corresponding adjoint Fokker–Planck equation together with its Robin boundary condition. We describe the comparison step in some detail, as the reflection term is precisely what gives rise to the boundary condition. Step 1: identification of the mean-field drift. For μ∈(Ω)μ ( ), define U[μ]()=defV()+λ∫Ωgs(−′)μ(d′),b[μ]()=def−∇U[μ](). U[μ]( ω) =V( ω)+λ _ g_s( ω- ω )\,μ(d ω ), b[μ]( ω) =-∇ U[μ]( ω). (G.5) By compactness of Ω and the assumed smooth regularization of gsg_s (or the collision-free restriction), there is a constant L<∞L<∞, independent of N, such that for all ,~∈Ω ω, ω∈ and μ,ν∈(Ω)μ,ν ( ), ‖b[μ]()−b[ν](~)‖≤L(‖−~‖+W1(μ,ν))≤L(‖−~‖+W2(μ,ν)). \|b[μ]( ω)-b[ν]( ω)\|≤ L (\| ω- ω\|+W_1(μ,ν) )≤ L (\| ω- ω\|+W_2(μ,ν) ). (G.6) Let μN(−k)=(N−1)−1∑ℓ≠kδℓ _N^(-k)=(N-1)^-1 _ ≠ k _ ω_ . Since gsg_s is even and the finite interaction is normalized as 1/(2N(N−1))∑k≠ℓgs(k−ℓ)1/(2N(N-1)) _k≠ g_s( ω_k- ω_ ), −N∇kℋN(N) -N _ ω_kH_N( _N) =−∇V(k)−λN−1∑ℓ≠k∇gs(k−ℓ)=b[μN(−k)](k). =-∇ V( ω_k)- λN-1 _ ≠ k∇ g_s( ω_k- ω_ )=b[ _N^(-k)]( ω_k). (G.7) Moreover ‖b[μN(−k)](k)−b[μN](k)‖≤C/N\|b[ _N^(-k)]( ω_k)-b[ _N]( ω_k)\|≤ C/N. Thus the algorithm is the projected Euler scheme for the empirical McKean drift, up to a uniformly vanishing O(N−1)O(N^-1) self-interaction error. Step 2: continuous-time embedding and the reflected diffusion comparison. Construct Brownian motions k B_k so that k((m+1)ηN)−k(mηN)=ηNkm B_k((m+1) _N)- B_k(m _N)= _N\, ξ_k^m. Let ¯kN,η(t)=km ω_k^N,η(t)= ω_k^m for t∈[mηN,(m+1)ηN)t∈[m _N,(m+1) _N). The recursion in (E.2) can be written as ¯kN,η(mηN)=Ω¯(¯kN,η((m−1)ηN)+ηNb[μ¯(m−1)ηN,(−k)](¯kN,η((m−1)ηN))+2/βΔkm), ω_k^N,η(m _N)=P_ \! ( ω_k^N,η((m-1) _N)+ _Nb[ μ_(m-1) _N^N,(-k)]( ω_k^N,η((m-1) _N))+ 2/β\, B_k^m ), (G.8) where μ¯tN,(−k)=(N−1)−1∑ℓ≠kδ¯ℓN,η(t) μ_t^N,(-k)=(N-1)^-1 _ ≠ k _ ω_ ^N,η(t) and Δkm=k(mηN)−k((m−1)ηN) B_k^m= B_k(m _N)- B_k((m-1) _N). This is the natural cadlag embedding of the projected particle chain. The limiting continuous reflected particle system associated with (G.8) is dkN(t) d X_k^N(t) =b[μtN,(−k)](kN(t))dt+2/βdk(t)−(kN(t))dLkN(t),μtN=1N∑j=1NδjN(t), =b[ _t^N,(-k)]( X_k^N(t))\,dt+ 2/β\,d B_k(t)- n( X_k^N(t))\,dL_k^N(t), _t^N= 1N _j=1^N _ X_j^N(t), (G.9) kN(t) X_k^N(t) ∈Ω¯,LkN(t) is nondecreasing,∫0TΩ(kN(t))LkN(t)=0. ∈ , L_k^N(t) is nondecreasing, _0^T1_ ( X_k^N(t))\,dL_k^N(t)=0. (G.10) Here n is the outward normal, so −dLkN- n\,dL_k^N is the inward reflection. The Skorokhod problem on a compact convex C2C^2 domain has a unique reflected solution, and the projection map in (G.8) is precisely the Euler approximation of this reflected equation. The projected-Euler consistency assumption in the theorem means that, for each fixed T, εN,η(T)=defsup0≤t≤T1N∑k=1N‖¯kN,η(t)−kN(t)‖2⟶0. _N,η(T) =E _0≤ t≤ T 1N _k=1^N\| ω_k^N,η(t)- X_k^N(t)\|^2 0. (G.11) For smooth bounded drifts this follows from the standard Euler–Skorokhod estimate; the condition on ηN _N makes the discrete projection error negligible on the mean-field scale. Step 3: nonlinear reflected process and propagation of chaos. Let k(t) X_k(t), k≥1k≥ 1, be i.i.d. copies of the nonlinear reflected McKean–Vlasov process dk(t) d X_k(t) =b[μt](k(t))dt+2/βdk(t)−(k(t))dLk(t),μt=Law(k(t)), =b[ _t]( X_k(t))\,dt+ 2/β\,d B_k(t)- n( X_k(t))\,dL_k(t), _t=Law( X_k(t)), (G.12) k(t) X_k(t) ∈Ω¯,∫0TΩ(k(t))Lk(t)=0. ∈ , _0^T1_ ( X_k(t))\,dL_k(t)=0. (G.13) The Lipschitz bound (G.6) gives existence and uniqueness by a fixed-point argument on measure-valued curves. Couple kN X_k^N and k X_k with the same Brownian motion and the same initial particle. For convex Ω , the reflection map is monotone: (kN(t)−k(t))⋅(−(kN(t))dLkN(t)+(k(t))dLk(t))≤0. ( X_k^N(t)- X_k(t) )· (- n( X_k^N(t))\,dL_k^N(t)+ n( X_k(t))\,dL_k(t) )≤ 0. (G.14) Applying Itô’s formula to ‖kN(t)−k(t)‖2\| X_k^N(t)- X_k(t)\|^2, using (G.6), and averaging over k yields eN(t) e_N(t) =defsup0≤r≤t1N∑k=1N‖kN(r)−k(r)‖2 =E _0≤ r≤ t 1N _k=1^N\| X_k^N(r)- X_k(r)\|^2 (G.15) ≤CT∫0teN(s)s+CT∫0tW22(1N∑j=1Nδj(s),μs)s+CTN2. ≤ C_T _0^te_N(s)\,ds+C_T _0^tEW_2^2\! ( 1N _j=1^N _ X_j(s), _s )ds+ C_TN^2. (G.16) The last term is the leave-one-out/self-interaction error. Since Ω is compact and j(s) X_j(s) are i.i.d. with law μs _s, aN(T)=defsup0≤s≤TW22(1N∑j=1Nδj(s),μs)⟶0. a_N(T) = _0≤ s≤ TEW_2^2\! ( 1N _j=1^N _ X_j(s), _s ) 0. (G.17) Gronwall’s inequality therefore gives sup0≤t≤TW22(μtN,μt)≤CT(aN(T)+N−1+εN,η(T))⟶0. _0≤ t≤ TEW_2^2( _t^N, _t)≤ C_T (a_N(T)+N^-1+ _N,η(T) ) 0. (G.18) This is the propagation-of-chaos statement: any fixed finite subcollection of particles becomes asymptotically independent, and each coordinate has law μt _t. Step 4: the Girsanov/change-of-measure ingredient. The preceding synchronous estimate can equivalently be written as a compact-domain change-of-measure estimate. To make this explicit, define k(t)=β/2(b[μtN](kN(t))−b[μt](kN(t))). u_k(t)= β/2 (b[ _t^N]( X_k^N(t))-b[ _t]( X_k^N(t)) ). (G.19) Novikov’s condition holds because Ω is compact and the drift is bounded. Hence the exponential martingale T=exp(−∑k=1N∫0Tk(t)⋅k(t)−12∑k=1N∫0T‖k(t)‖2t) _T= \! (- _k=1^N _0^T u_k(t)· d B_k(t)- 12 _k=1^N _0^T\| u_k(t)\|^2dt ) (G.20) defines a probability measure under which ~k(t)=k(t)+∫0tk(s)s B_k(t)= B_k(t)+ _0^t u_k(s)ds are Brownian motions. The relative entropy of this tilted law with respect to the original law is DKL(ℙ~N,T∥ℙN,T)=12∑k=1N~∫0T‖k(t)‖2t≤CN∫0T~W22(μtN,μt)t. D_KL( P_N,T\|P_N,T)= 12 _k=1^N E _0^T\| u_k(t)\|^2dt≤ C\,N _0^T EW_2^2( _t^N, _t)dt. (G.21) By Pinsker’s inequality and the bounded diameter of Ω , this entropy bound gives, at the empirical-measure level, W22(μtN,emp,μtN,nl)≤diam(Ω)2(2NDKL(ℙ~N,t∥ℙN,t))1/2, _2^2( _t^N,emp, _t^N,nl) ( )^2 ( 2ND_KL( P_N,t\|P_N,t) )^1/2, (G.22) where μtN,emp _t^N,emp and μtN,nl _t^N,nl denote the empirical measures of the empirical-drift and decoupled reflected systems. This quantifies the cost of replacing the particle drift by the law-dependent McKean drift. Combining (G.11), (G.18), (G.22), and the triangle inequality gives sup0≤t≤TW22(1N∑k=1Nδk⌊t/ηN⌋,μt)⟶0. _0≤ t≤ TEW_2^2\! ( 1N _k=1^N _ ω_k t/ _N , _t ) 0. (G.23) Since W2W_2 convergence on compact Ω implies weak convergence, this proves the empirical-measure convergence in both theorem statements. Step 5: identification of the McKean–Vlasov PDE and the boundary condition. Let ψ∈C2(Ω¯)ψ∈ C^2( ) belong to the generator domain of the reflected diffusion, i.e. ∂ψ=0 _ nψ=0 on ∂Ω∂ . Applying Itô’s formula to ψ(t)ψ( X_t) in (G.12) gives dt∫Ωψ()μt(d)=∫Ω[b[μt]()⋅∇ψ()+β−1Δψ()]μt(d). ddt _ ψ( ω)\, _t(d ω)= _ [b[ _t]( ω)·∇ψ( ω)+β^-1 ψ( ω) ] _t(d ω). (G.24) The local-time term is −∂ψ(t)dLt- _ nψ( X_t)dL_t, and therefore vanishes for this generator domain. If μt(d)=ρt()d _t(d ω)= _t( ω)d ω, b[μt]=−∇Utb[ _t]=-∇ U_t, and Ut=U[μt]U_t=U[ _t], (G.24) becomes dt∫Ωψρt=∫Ω[−∇Ut⋅∇ψ+β−1Δψ]ρt. ddt _ ψ _t\,d ω= _ [-∇ U_t·∇ψ+β^-1 ψ ] _t\,d ω. (G.25) The adjoint of the reflected generator is therefore ∂tρt=−∇⋅(ρtb[μt])+β−1Δρt=∇⋅(ρt∇Ut)+β−1Δρtin Ω. _t _t=-∇·( _tb[ _t])+β^-1 _t=∇·( _t∇ U_t)+β^-1 _t . (G.26) To identify the boundary condition, integrate the last display against arbitrary smooth test functions and use Green’s formula. The boundary contribution is ∫∂Ωψ()(ρt()b[μt]()−β−1∇ρt())⋅()S(). _∂ ψ( ω) ( _t( ω)b[ _t]( ω)-β^-1∇ _t( ω) )· n( ω)\,dS( ω). (G.27) Reflection means that the probability current through ∂Ω∂ vanishes. Since b[μt]=−∇Utb[ _t]=-∇ U_t, this current condition is (ρt∇Ut+β−1∇ρt)⋅=0on (0,T]×∂Ω. ( _t∇ U_t+β^-1∇ _t )· n=0 (0,T]×∂ . (G.28) For smooth ρt _t, (G.28) is exactly the Robin form ∂ρt+βρt∂Ut=0on (0,T]×∂Ω, _ n _t+β _t _ nU_t=0 (0,T]×∂ , (G.29) which is exactly the Robin/no-flux boundary condition generated by the reflected dynamics. The initial condition follows from μ0N⇒ρ0d _0^N _0d ω. Taking ψ≡1ψ≡ 1 in the weak formulation gives ∫Ωρt=1 _ _t=1, and nonnegativity follows because ρt _t is the density of the law of the reflected process. Finally, sending β→∞β→∞ removes the diffusion term and reduces the no-flux condition to (ρtt)⋅=0( _t v_t)· n=0, where t=−∇Ut v_t=-∇ U_t, yielding the deterministic reflected continuity equation. G.4 Proof of Theorem E.3 The proof of Theorem E.3 follows standard application of Varadhan’s lemma Ellis (2005). We provide the proof in multiple steps. G.4.1 Rewrite the Hamiltonian as a functional of μN _N Let ν denote Lebesgue measure on Ω and set ρ=defν(Ω).ρ =\; ν( ). We use ρ as reference probability measure, i.e., ρ∈(Ω)ρ ( ). We also recall the definition of empirical measure μN=1N∑k=1Nδk _N= 1N _k=1^N _ ω_k. Under the product measure ρ⊗Nρ N on ΩN ^N, the empirical measure μN _N satisfies Sanov’s theorem: it obeys an LDP on (Ω)P( ) with speed N and good rate function I(μ)=Ent(μ∣ρ).I(μ)=Ent(μ ρ). (G.30) Now, define the bounded measurable function V()=def−1n(n−1)∑0≤i≠j≤nyiyjϕ(i)ϕ(j).V( ω) =- 1n(n-1) _0≤ i≠ j≤ ny_iy_j _ ω( x_i) _ ω( x_j). (G.31) Then by the Monte–Carlo substitution in (B) and the definition (B.4), the empirical loss is exactly a linear functional of μN _N: ℰN(N)=−1n(n−1)∑i≠jyiyj∫Ωϕ(i)ϕ(j)μN(d)=∫ΩV()μN(d),E_N( _N)=- 1n(n-1) _i≠ jy_iy_j _ _ ω( x_i)\, _ ω( x_j)\, _N(d ω)= _ V( ω)\, _N(d ω), (G.32) Assume the interaction term is of mean-field form (as in the Gibbs law (E.7)): N,s(N)=def12N(N−1)∑1≤k≠ℓ≤Ngs(k−ℓ).W_N,s( _N) = 12N(N-1) _1≤ k≠ ≤ Ng_s( ω_k- ω_ ). (G.33) Define the corresponding continuum interaction functional on (Ω)P( ): s(μ)=def12∫Ω∫Ωgs(−′)μ(d)μ(d′).W_s(μ) = 12 _ _ g_s( ω- ω )\,μ(d ω)μ(d ω ). (G.34) Then the full empirical Hamiltonian in (B.3) can be written as ℋN,s(N)=ℰN(N)+λN,s(N).H_N,s( _N)=E_N( _N)+λ\,W_N,s( _N). (G.35) With this normalization, N,s(N)W_N,s( _N) converges to s(μN)W_s( _N); after truncation the difference is only a harmless finite-size correction. To handle possible singularities of the kernel gsg_s at 0, we introduce a truncation. G.4.2 Truncation of the interaction kernel and removal of the diagonal constant The Columb/Riesz kernel gsg_s is singular at the origin. In order to work with bounded continuous functionals on (Ω)P( ) (so that Varadhan’s lemma applies), we introduce a truncation. For ε≥0 ≥ 0 define the truncated kernel gsε()=defmin‖−s,ε−s g_s ( ω) = \\| ω\|^-s, ^-s\ (G.36) and adopt the convention that gsε()g_s ( 0) is the (finite) value of this truncation at 0. In particular, gsε()=ε−sg_s ( 0)= ^-s. Define the truncated empirical interaction by N,sε(ΩN)=def12N(N−1)∑1≤k≠ℓ≤Ngsε(k−ℓ),W_N,s ( _N) = 12N(N-1) _1≤ k≠ ≤ Ng_s ( ω_k- ω_ ), (G.37) and the corresponding truncated mean-field interaction functional on (Ω)P( ) by sε(μ)=def12∫Ω∫Ωgsε(−′)μ(d)μ(d′).W_s (μ) = 12 _ _ g_s ( ω- ω )\,μ(d ω)\,μ(d ω ). (G.38) For fixed ε , the map μ↦sε(μ)μ _s (μ) is continuous on (Ω)P( ) since gsεg_s is bounded and continuous on ΩN _N. Define the truncated energy functional on (Ω)P( ) by ℋsε(μ)=def∫ΩV()μ(d)+λsε(μ).H_s (μ) = _ V( ω)\,μ(d ω)+λ\,W_s (μ). (G.39) and define the truncated empirical Hamiltonian by ℋN,sε(N)=defℰN(N)+λN,sε(N).H_N,s ( _N) =E_N( _N)+λ\,W_N,s ( _N). (G.40) Lemma G.2 (Finite-size correction for empirical measures). For every configuration N=(1,…,N)∈ΩN _N=( ω_1,…, ω_N)∈ ^N, sε(μN)=N−1NN,sε(N)+gsε()2N.W_s ( _N)= N-1NW_N,s ( _N)+ g_s ( 0)2N. (G.41) Equivalently, ℋN,sε(N)=ℋsε(μN)+λN−1sε(μN)−λgsε()2(N−1).H_N,s ( _N)=H_s ( _N)+ λN-1W_s ( _N)-λ\, g_s ( 0)2(N-1). (G.42) In particular, for fixed ε>0 >0, supμ∈(Ω)|λN−1sε(μ)|≤|λ|ε−s2(N−1). _μ ( ) | λN-1W_s (μ) |≤ |λ|\, ^-s2(N-1). Thus replacing ℋN,sε(N)H_N,s ( _N) by ℋsε(μN)H_s ( _N) changes the normalized logarithmic Laplace limits below by o(1)o(1) at speed N when βN/N→1 _N/N→ 1 and by o(1)o(1) at speed βN _N when βN/N→∞ _N/N→∞. Proof G.3. Since μN=1N∑k=1Nδk _N= 1N _k=1^N _ ω_k, we have sε(μN)=121N2∑k,ℓ=1Ngsε(k−ℓ)=121N2∑k≠ℓgsε(k−ℓ)+gsε()2N.W_s ( _N)= 12 1N^2 _k, =1^Ng_s ( ω_k- ω_ )= 12 1N^2 _k≠ g_s ( ω_k- ω_ )+ g_s ( 0)2N. The off-diagonal term equals N−1N12N(N−1)∑k≠ℓgsε(k−ℓ)=N−1NN,sε(N), N-1N\, 12N(N-1) _k≠ g_s ( ω_k- ω_ )= N-1NW_N,s ( _N), which proves (G.41). Solving this identity for N,sεW_N,s and adding the linear term ℰN=∫VμNE_N= V\,d _N gives (G.42). The uniform bound follows from 0≤gsε≤ε−s0≤ g_s ≤ ^-s. The deterministic constant in (G.42) is absorbed into the partition function. The remaining tilt is uniformly O(1/N)O(1/N); after multiplication by βN _N and division by the relevant speed, its contribution is O(1/N)O(1/N) both in the thermal scale βN∼N _N N and in the zero-temperature scale βN _N. Lemma G.2 shows that the N(N−1)N(N-1)-normalized Hamiltonian and the mean-field functional ℋsε(μN)H_s ( _N) have the same logarithmic Laplace limits at the speeds used below. The diagonal constant is absorbed by the partition function, and the remaining finite-size tilt is uniformly negligible. In what follows we therefore use the asymptotically equivalent representation ℙN,βNε(dN)≍1Z~N,βNεexp(−βNℋsε(μN))ρ⊗N(dN),P_N, _N (d _N) 1 Z_N, _N \! (- _N\,H_s ( _N) )\,ρ N(d _N), where ≍ denotes equivalence of the normalized logarithmic Laplace limits at the relevant speed. G.4.3 Case βN/N→1 _N/N→ 1 Assume βN/N→1 _N/N→ 1 and write βN=NαN _N=N _N with αN=ON(1) _N=O_N(1). We can rewrite the truncated Gibbs law as a tilt by ℋsε(μN)H_s ( _N): ℙN,βNε(dΩN)∝exp(−NαNℋsε(μN))ρ⊗N(dN).P _N, _N(d _N) \! (-N _N\,H_s ( _N) )\,ρ N(d _N). (G.43) Let F:(Ω)→ℝF:P( ) be bounded and continuous. We leverage Varadhan’s lemma (Ellis, 2005, p. 51), a rigorous formulation of the Laplace principle (or the saddle point technique) applied to measures satisfying a large deviations property: Lemma G.4 (Varadhan’s Lemma Ellis (2005)). Suppose a sequence QNN=1∞\Q_N\_N=1^∞ of probability measures on X satisfies a large deviations property with rate function I(x)I(x). Let F:→ℝF:X be a continuous function that satisfies the tail condition limL→∞lim supN→∞1Nlog∫x:F(x)≥Lexp(NF(x))QN(dx)=−∞. _L→∞ _N→∞ 1N _x:F(x)≥ L (NF(x) )\,Q_N(dx)=-∞. (G.44) Then limN→∞1Nlog∫exp(NF(x))QN(dx)=supx∈F(x)−I(x). _N→∞ 1N _X (NF(x) )\,Q_N(dx)= _x \F(x)-I(x) \. (G.45) To apply Varadhan’s lemma, we must verify the tail condition. Let QN=defρ⊗N∘μN−1Q_N\; =\;ρ N _N^-1 denote the push-forward law of the empirical measure μN _N on (Ω)P( ) under ρ⊗Nρ N. Concretely, for any measurable A⊂(Ω)A ( ), QN(A)=ρ⊗N(N:μN(N)∈A).Q_N(A)=ρ N (\ _N: _N( _N)∈ A\ ). The tail condition in Varadhan’s lemma for a functional Φ:(Ω)→ℝ :P( ) reads limL→∞lim supN→∞1Nlog∫μ:Φ(μ)≥LeNΦ(μ)QN(dμ)=−∞. _L→∞ _N→∞ 1N _\μ:\, (μ)≥ L\e^N (μ)\,Q_N(dμ)\;=\;-∞. (G.46) Equivalently, since μN _N has law QNQ_N under ρ⊗Nρ N, limL→∞lim supN→∞1Nlogρ⊗N[eNΦ(μN)Φ(μN)≥L]=−∞. _L→∞ _N→∞ 1N _ρ N [e^N ( _N)1_\ ( _N)≥ L\ ]\;=\;-∞. (G.47) In our application, we use the continuous functional for numerator ΦNε(μ)=def−(αNℋsε(μ)+F(μ)),αN=βN/N→1, _N(μ)\; =\;- ( _N\,H_s (μ)+F(μ) ), _N= _N/N→ 1, (G.48) and for the denominator the functional ΨNε(μ)=def−αNℋsε(μ). _N(μ)\; =\;- _N\,H_s (μ). (G.49) We verify (G.46) for ΦNε _N. The same argument applies to ΨNε _N. To do so, we first show that ℋsεH_s is bounded for fixed ε>0 >0. Recall ℋsε(μ)=∫ΩV()μ(d)+λsε(μ),withsε(μ)=12∬gsε(−′)μ(d)μ(d′).H_s (μ)= _ V( ω)\,μ(d ω)\;+\;λ\,W_s (μ), \ W_s (μ)= 12 g_s ( ω- ω )\,μ(d ω)μ(d ω ). Since V is bounded, for all μ∈(Ω)μ ( ), |∫ΩVdμ| | _ V\,dμ | ≤‖V‖∞ ≤\|V\|_∞ (G.50) =1n(n−1)‖∑0≤i≠j≤nyiyjϕ(i)ϕ(j)‖∞ = 1n(n-1) \| _0≤ i≠ j≤ ny_iy_j _ ω( x_i) _ ω( x_j) \|_∞ (G.51) ≤|yiyj|‖ϕ(i)‖∞‖ϕ(j)‖∞ ≤|y_iy_j| \| _ ω( x_i)\|_∞\| _ ω( x_j) \|_∞ (G.52) ≤Lϕ2. ≤ L_φ^2. (G.53) Moreover, by definition of truncated kernel gsε()g_s ( ω) in Eq. (G.36), we have gsε≤ε−sg_s ≤ ^-s. Therefore, for any μ∈(Ω)μ ( ), |sε(μ)|=12|∬gsε(ω−ω′)μ(dω)μ(dω′)|≤12|∬ε−sμ(dω)μ(dω′)|=ε−s2.|W_s (μ)|= 12 | g_s (ω-ω )\,μ(dω)μ(dω ) |≤ 12 | ^-s\,μ(dω)μ(dω ) |= ^-s2. (G.54) Combining (G.51)–(G.54) gives, for all μ, |ℋsε(μ)|≤Lϕ2+|λ|ε−s2=defCλ,ε.|H_s (μ)|≤ L_φ^2+|λ|\, ^-s2\; =\;C_λ, . (G.55) Since F is bounded and continuous, define ‖F‖∞=defsupμ∈(Ω)|F(μ)|<∞.\|F\|_∞\; =\; _μ ( )|F(μ)|<∞. Moreover, since αN=N(1) _N=O_N(1), there exists N0N_0 and constant C>0C>0 such that for all N≥N0N≥ N_0, |αN|≤C.| _N|≤ C. (G.56) Then for all N≥N0N≥ N_0 and all μ∈(Ω)μ ( ), |ΦN(μ)| | _N(μ)| =|−αNℋsε(μ)−F(μ)| = |- _NH_s (μ)-F(μ) | (G.57) ≤|αN||ℋsε(μ)|+|F(μ)| ≤| _N|\,|H_s (μ)|+|F(μ)| (G.58) ≤CCλ,ε+‖F‖∞=defKε,F, ≤ C_λ, +\|F\|_∞\; =\;K_ ,F, (G.59) where we used (G.55) and (G.56). Thus, ΦN _N is bounded above by Kε,FK_ ,F uniformly in μ, for all N≥N0N≥ N_0. Similarly, for the denominator functional, ΨN(μ)=−αNℋsε(μ)≤|αN||ℋsε(μ)|≤CCε=defKε,N≥N0. _N(μ)=- _NH_s (μ)≤| _N|\,|H_s (μ)|≤ C_ \; =\;K_ , N≥ N_0. (G.60) Fix any L>KM,FL>K_M,F. Then by (G.59), for all N≥N0N≥ N_0, μ∈(Ω):ΦN(μ)≥L= .\μ ( ): _N(μ)≥ L\= . Hence the tail integral vanishes, i.e., ∫μ:ΦN(μ)≥LeNΦN(μ)QN(dμ)=0,N≥N0. _\μ: _N(μ)≥ L\e^N _N(μ)\,Q_N(dμ)=0, N≥ N_0. Therefore, 1Nlog∫μ:ΦN(μ)≥LeNΦN(μ)QN(dμ)=−∞,N≥N0. 1N _\μ: _N(μ)≥ L\e^N _N(μ)\,Q_N(dμ)=-∞, N≥ N_0. This establishes the tail condition (G.46). limL→∞lim supN→∞1Nlog∫μ:ΦN(μ)≥LeNΦN(μ)QN(dμ)=−∞, _L→∞ _N→∞ 1N _\μ: _N(μ)≥ L\e^N _N(μ)\,Q_N(dμ)=-∞, (G.61) The same reasoning, using (G.60), shows that the tail condition also holds for the denominator functional ΨN _N. Now, by Varadhan’s lemma applied to Sanov’s LDP (G.30), we obtain limN→∞1Nlogρ⊗N[exp(−N(αNℋsε+F)(μN))] _N→∞ 1N _ρ N\! [ (-N( _NH_s +F)( _N) ) ] =−infμ∈(Ω)(ℋsε(μ)+F(μ)+Ent(μ∣ρ)), =- _μ ( ) (H_s (μ)+F(μ)+Ent(μ ρ) ), (G.62) limN→∞1Nlogρ⊗N[exp(−NαNℋsε(μN))] _N→∞ 1N _ρ N\! [ (-N _NH_s ( _N) ) ] =−infμ∈(Ω)(ℋsε(μ)+Ent(μ∣ρ)). =- _μ ( ) (H_s (μ)+Ent(μ ρ) ). (G.63) Subtracting (G.63) from (G.62) yields the Laplace principle for μN _N under ℙN,βNεP_N, _N at speed N: limN→∞−1NlogℙN,βNε[e−NF(μN)]=infμ(F(μ)+sε(μ)), _N→∞- 1N _P_N, _N \! [e^-NF( _N) ]= _μ (F(μ)+J _s(μ) ), with sε(μ)=(ℰsε(μ)+Ent(μ∣ρ))−infν∈(Ω)(ℰsε(ν)+Ent(ν∣ρ)).J _s(μ)= (E_s (μ)+Ent(μ ρ) )- _ν ( ) (E_s (ν)+Ent(ν ρ) ). (G.64) Since (Ω)P( ) is a Polish space and the Laplace principle holds for all F∈Cb((Ω))F∈ C_b(P( )), Bryc’s inverse Varadhan lemma (see, e.g., (Dembo and Zeitouni, 2009, Theorem 4.4.13)) implies that (μN)N≥1( _N)_N≥ 1 satisfies a large deviation principle with speed N and good rate function sεJ_s . Finally let ε↓0 0. Because gsε↑gsg_s g_s, we have ℋsε(μ)→ℋs(μ)H_s (μ) _s(μ) pointwise (possibly +∞+∞), so ℋsH_s is lower semi-continuous as a supremum of continuous functions. Standard monotone truncation arguments for Laplace principles yield the same Laplace limit with ℋsH_s in place of ℋsεH_s , and hence the LDP with rate s(μ)=(ℋs(μ)+Ent(μ∣ρ))−infγ(ℋs(γ)+Ent(γ∣ρ)).J_s(μ)= (H_s(μ)+Ent(μ ρ) )- _γ (H_s(γ)+Ent(γ ρ) ). To convert to Ent(⋅∣ν)Ent(· ν), note that Ent(μ∣ρ)=Ent(μ∣ν)+logν(Ω)Ent(μ ρ)=Ent(μ ν)+ ν( ), and the additive constant cancels when subtracting the infimum. This gives the first bullet of the theorem. G.4.4 Case βN/N→∞ _N/N→∞ Fix ε≥0 ≥ 0 and recall that under the equivalent representation (Lemma G.2), ℙN,βNε(dN)=1Z~N,βNεexp(−βNℋsε(μN))ρ⊗N(dN),P_N, _N (d _N)= 1 Z_N, _N \! (- _N\,H_s ( _N) )\,ρ N(d _N), where the partition function is Z~N,βNε=defρ⊗N[exp(−βNℋsε(μN))]. Z_N, _N \; =\;E_ρ N\! [ \! (- _N\,H_s ( _N) ) ]. As before, let QN=defρ⊗N∘μN−1Q_N\; =\;ρ N _N^-1 denote the law of μN _N on (Ω)P( ) under the reference measure ρ⊗Nρ N. We notice that Sanov probabilities are negligible at speed βN/N→∞ _N/N→∞. The key point is that QNQ_N satisfies an LDP at speed N, so any QNQ_N-probability is at worst exp(−cN) (-cN), which is negligible on the scale βN≫N _N N. Lemma G.5 (Sanov scale is negligible at speed βN _N). Assume βN/N→∞ _N/N→∞. Then for every nonempty open set U⊂(Ω)U ( ), limN→∞1βNlogQN(U)=0. _N→∞ 1 _N Q_N(U)=0. Proof G.6. Since QN(U)≤1Q_N(U)≤ 1, we have lim supN→∞1βNlogQN(U)≤0 _N→∞ 1 _N Q_N(U)≤ 0. For the lower bound, note that ρ has full support on Ω and the set μ∈(Ω):Ent(μ∣ρ)<∞\μ ( ):Ent(μ ρ)<∞\ is dense in (Ω)P( ) for the weak topology, so any nonempty open U contains some μ with finite entropy. Hence infμ∈UEnt(μ∣ρ)<∞. _μ∈ UEnt(μ ρ)<∞. By the lower bound in Sanov’s theorem, for open U, lim infN→∞1NlogQN(U)≥−infμ∈UEnt(μ∣ρ)>−∞. _N→∞ 1N Q_N(U)\;≥\;- _μ∈ UEnt(μ ρ)\;>\;-∞. Therefore, lim infN→∞1βNlogQN(U)=lim infN→∞NβN⋅1NlogQN(U)≥ 0, _N→∞ 1 _N Q_N(U)= _N→∞ N _N· 1N Q_N(U)\;≥\;0, because N/βN→0N/ _N→ 0 and 1NlogQN(U) 1N Q_N(U) is bounded below along a subsequence by a finite constant. Combining with the lim sup≤0 ≤ 0 gives the desired limit 0. Now, define the (truncated) ground state energy eε⋆=definfμ∈(Ω)ℋsε(μ).e_ \; =\; _μ ( )H_s (μ). Since gsεg_s and V are bounded, ℋsεH_s is bounded and continuous on (Ω)P( ), hence eε⋆∈ℝe_ . Lemma G.7 (Partition function at scale βN _N). Under βN/N→∞ _N/N→∞, limN→∞1βNlogZ~N,βNε=−eε⋆. _N→∞ 1 _N Z_N, _N =-\,e_ . Proof G.8. Upper bound. Since ℋsε(μ)≥eε⋆H_s (μ)≥ e_ for all μ, Z~N,βNε=ρ⊗N[e−βNℋsε(μN)]≤e−βNeε⋆, Z_N, _N =E_ρ N\! [e^- _NH_s ( _N) ]≤ e^- _Ne_ , so lim supN→∞1βNlogZ~N,βNε≤−eε⋆ _N→∞ 1 _N Z_N, _N ≤-e_ . Lower bound. Fix ε>0 >0 and choose με∈(Ω) _ ( ) such that ℋsε(με)≤eε⋆+εH_s ( _ )≤ e_ + . By continuity of ℋsεH_s there exists a nonempty open neighborhood U of με _ such that supμ∈Uℋsε(μ)≤eε⋆+2ε. _μ∈ UH_s (μ)≤ e_ +2 . Then Z~N,βNε≥ρ⊗N[μN∈Ue−βNℋsε(μN)]≥e−βN(eε⋆+2ε)QN(U). Z_N, _N _ρ N\! [1_\ _N∈ U\\,e^- _NH_s ( _N) ]≥ e^- _N(e_ +2 )\,Q_N(U). Taking log and dividing by βN _N yields 1βNlogZ~N,βNε≥−(eε⋆+2ε)+1βNlogQN(U). 1 _N Z_N, _N ≥-(e_ +2 )+ 1 _N Q_N(U). By Lemma G.5, 1βNlogQN(U)→0 1 _N Q_N(U)→ 0, hence lim infN→∞1βNlogZ~N,βNε≥−(eε⋆+2ε). _N→∞ 1 _N Z_N, _N ≥-(e_ +2 ). Letting ε↓0 0 gives the lower bound ≥−eM⋆≥-e_M and completes the proof. Let F⊂(Ω)F ( ) be closed. Using the ratio representation, ℙN,βNε(μN∈F)=ρ⊗N[μN∈Fexp(−βNℋsε(μN))]Z~N,βNε.P_N, _N ( _N∈ F)= E_ρ N\! [1_\ _N∈ F\ (- _NH_s ( _N) ) ] Z_N, _N . On the event μN∈F\ _N∈ F\ we have exp(−βNℋsε(μN))≤exp(−βNinfμ∈Fℋsε(μ)) (- _NH_s ( _N) )≤ (- _N _μ∈ FH_s (μ) ), hence ρ⊗N[μN∈Fexp(−βNℋsε(μN))]≤exp(−βNinfμ∈Fℋsε(μ)).E_ρ N\! [1_\ _N∈ F\ (- _NH_s ( _N) ) ]≤ \! (- _N _μ∈ FH_s (μ) ). Therefore, lim supN→∞1βNlogℙN,βNε(μN∈F)≤−infμ∈Fℋsε(μ)−lim infN→∞1βNlogZ~N,βNε. _N→∞ 1 _N _N, _N ( _N∈ F)≤- _μ∈ FH_s (μ)- _N→∞ 1 _N Z_N, _N . By Lemma G.7, limN1βNlogZ~N,βNε=−eε⋆ _N 1 _N Z_N, _N =-e_ , so lim supN→∞1βNlogℙN,βNε(μN∈F)≤−(infμ∈Fℋsε(μ)−eε⋆). _N→∞ 1 _N _N, _N ( _N∈ F)≤- ( _μ∈ FH_s (μ)-e_ ). (G.65) Now, let G⊂(Ω)G ( ) be open and fix μ∈Gμ∈ G. By continuity of ℋsεH_s , there exists a nonempty open neighborhood U⊂GU⊂ G of μ such that supν∈Uℋsε(ν)≤ℋsε(μ)+ε. _ν∈ UH_s (ν) _s (μ)+ . Then ℙN,βNε(μN∈G) _N, _N ( _N∈ G) ≥ℙN,βNε(μN∈U) _N, _N ( _N∈ U) (G.66) =ρ⊗N[μN∈Ue−βNℋsε(μN)]Z~N,βNε = E_ρ N\! [1_\ _N∈ U\e^- _NH_s ( _N) ] Z_N, _N (G.67) ≥e−βN(ℋsε(μ)+ε)QN(U)Z~N,βNε. ≥ e^- _N(H_s (μ)+ )\,Q_N(U) Z_N, _N . (G.68) Using the simple bound Z~N,βNε≤e−βNeε⋆ Z_N, _N ≤ e^- _Ne_ we obtain ℙN,βNε(μN∈G)≥exp(−βN(ℋsε(μ)−eε⋆+ε))QN(U).P_N, _N ( _N∈ G)≥ \! (- _N(H_s (μ)-e_ + ) )\,Q_N(U). Taking log and dividing by βN _N gives lim infN→∞1βNlogℙN,βNε(μN∈G)≥−(ℋsε(μ)−eε⋆+ε)+lim infN→∞1βNlogQN(U). _N→∞ 1 _N _N, _N ( _N∈ G)≥-(H_s (μ)-e_ + )+ _N→∞ 1 _N Q_N(U). By Lemma G.5, the last term is 0, hence lim infN→∞1βNlogℙN,βNε(μN∈G)≥−(ℋsε(μ)−eε⋆+ε). _N→∞ 1 _N _N, _N ( _N∈ G)≥-(H_s (μ)-e_ + ). Letting ε↓0 0 and then taking the infimum over μ∈Gμ∈ G yields lim infN→∞1βNlogℙN,βNε(μN∈G)≥−(infμ∈Gℋs,M(μ)−eε⋆). _N→∞ 1 _N _N, _N ( _N∈ G)≥- ( _μ∈ GH_s,M(μ)-e_ ). (G.69) Combining Eqs. (G.65) and (G.69) and applying Definition (E.2) establishes that (μN)N≥1( _N)_N≥ 1 satisfies an LDP under ℙN,βNεP_N, _N with speed βN _N and good rate function sε(μ)=defℋsε(μ)−eε⋆=ℋsε(μ)−infν∈(Ω)ℋsε(ν).J _s(μ)\; =\;H_s (μ)-e_ =H_s (μ)- _ν ( )H_s (ν). Finally, let ε↓0 0. Since gsε↑gsg_s g_s, we have ℋsε(μ)↑ℋs(μ)H_s (μ) _s(μ) pointwise (possibly +∞+∞), and eε⋆=infμℋsε(μ)↑infμℋs(μ)=defe⋆.e_ = _μH_s (μ) _μH_s(μ) =e . Thus sε(μ)↑s(μ)J _s(μ) _s(μ) pointwise, where s(μ)=defℋs(μ)−e⋆=ℋs(μ)−infν∈(Ω)ℋs(ν),J_s(μ)\; =\;H_s(μ)-e =H_s(μ)- _ν ( )H_s(ν), which is the claimed rate function in the regime βN/N→∞ _N/N→∞. G.4.5 goodness of rate function. In the first regime, Ent(⋅∣ℓ)Ent(· ) has compact level sets (Sanov rate is good), and ℰsE_s is lower semicontinuous, so ℰs+Ent(⋅∣ℓ)E_s+Ent(· ) has compact sublevel sets, hence sJ_s is good. In the second regime, ℰsE_s is lower semicontinuous and the underlying space is bounded, so the sublevel sets of ℰsE_s are tight; together with lower semicontinuity this yields goodness.