Asymmetric Scaling Laws from Sparse Features
arXiv:2605.23591v1 [stat.ML] 22 May 2026
John Sous∗ Department of Applied Physics, Yale University, New Haven, Connecticut 06511, USA & Energy Sciences Institute, Yale University, West Haven, Connecticut 06516, USA [email protected] Michael Winer† Institute for Advanced Study, Princeton, NJ 08540, USA [email protected] May 25, 2026 Abstract We introduce a model for neural scaling laws under sparse activations. In the model, test loss is often dominated by rare coordinates that are never observed in the training input. This mechanism induces a novel bottleneck absent from dense models. We derive the asymptotic population loss in both the underparameterized and overparameterized regimes, and show that the loss exhibits a double-descent peak near the interpolation threshold—where the number of parameters is just sufficient to fit the training data—resulting in a loss curve governed by two distinct scaling exponents—one for the overparameterized regime and one for the underparameterized regime—with a gap determined by the degree of sparsity. Additionally, we derive a compute-optimal frontier that favors increasing dataset size over model capacity under fixed compute budgets. We also analyze gradient-descent dynamics and identify a scaling law for the probability that fixed-step gradient descent becomes unstable. We further show that the sparsity-induced effect persists under nonlinear activations. Experiments validating the theory can be found at SparseScaling.
1
Introduction
Scaling laws describe predictable relationships between model performance and key resource variables: model size (number of parameters) N , dataset size (number of training examples or tokens) D, and compute budget C. These empirical laws provide crucial guidance for efficiently allocating resources when training large machine learning models, particularly large language models (LLMs). Foundational work by Kaplan et al. showed that test loss ℓ decreases as a power law with respect to N and D when scaling one while holding the other fixed [Kaplan et al., 2020]. Building on this, Hoffmann et al. showed that under the constraint of a fixed total compute budget C ∝ N D, both N and D should be scaled jointly in proportion to minimize ℓ [Hoffmann et al., 2022]. This leads to a compute-optimal scaling law that offers a principled prescription for balancing model and data scale in large-scale training. ∗ †
Author to whom correspondence should be addressed. Present affiliation: Alignment Research Center.
1
Despite intense research, a comprehensive theoretical understanding of scaling laws remains elusive. One approach to this challenge models the data as Gaussian with power-law covariance, applies a random embedding into a lower-dimensional representation, and uses linear regression to characterize the resulting scaling laws [Bordelon et al., 2020, Spigler et al., 2020, Bahri et al., 2021, Maloney et al., 2022]. This framework [Maloney et al., 2022] has been shown, both theoretically and empirically, to 1 α produce the scaling law ℓ ∝ N1 + D and the “Chinchilla" compute-optimal allocation in which, under a compute budget C ∝ N D, the optimal choice scales as N ∗ (C) ≍ D∗ (C) ≍ C 1/2 .1 An important aspect not determined by a framework that yields a compute-optimal exponent close to 12 is whether the power-law exponent governing loss decay with model size, αN , matches that governing decay with data size, αD , as empirical studies suggest need not be the case [Hoffmann et al., 2022].2 In this work, we propose a simple theoretical model that reproduces the observed scaling asymmetry and naturally explains the emergence of two distinct exponents in the population loss.3 Our model is a random embedding followed by a linear readout, with a sparse input activation structure such that only a subset of coordinates of x is active. This setup is motivated by the idea that the data-generating process excites only a sparse subset of high-dimensional feature directions. The resulting loss displays distinct scaling with N and D. We derive an asymptotic scaling law in which the difference between αN and αD arises continuously as the sparsity level is varied.4 Our main contributions are as follows. • We introduce the notion of scaling laws for sparse activations by introducing a model that captures the impact of sparsity on optimization and resulting scaling behavior. • We derive the population loss for the sparse random feature model, yielding a two-exponent scaling law with an intrinsic asymmetry between the underparameterized and overparameterized regimes and a double-descent peak near the interpolation threshold [Geiger et al., 2019, Belkin et al., 2019, 2020, Hastie et al., 2022, Bartlett et al., 2020, Mei and Montanari, 2022, Nakkiran et al., 2020]. • We derive the compute-optimal frontier under a fixed compute budget, and show that increasing sparsity improves compute-efficiency while shifting the optimal allocation toward larger datasets. • We analyze the training dynamics and convergence properties of the loss during optimization, deriving a scaling law for the failure of fixed-step gradient descent (GD). • We experimentally verify in a nonlinear two-layer network that scaling asymmetry continues to arise from sparsity rather than from nonlinearity. Related Work. Empirical scaling laws show power-law loss improvement with model size and data, and motivate compute-optimal training rules [Kaplan et al., 2020, Hoffmann et al., 2022]. On the theory side, solvable models based on random features, kernels, and high-dimensional regression derive closed-form or deterministic-equivalent scaling predictions [Bahri et al., 2021, Maloney et al., 2022, Bordelon et al., 2020, Spigler et al., 2020, Defilippis et al., 2024, Lin et al., 2024]. Complementary dynamical models study protocol- and time-dependent scaling (e.g., one-pass vs. multi-pass training) and can exhibit multiple scaling phases [Bordelon et al., 2024, Paquette et al., 2024, Bordelon et al., 1
This compute-allocation law arises in this framework both in the Bayes-optimal [Bahri et al., 2021, Maloney et al., 2022] and one-pass SGD [Bordelon et al., 2024, Paquette et al., 2024] settings. In the Bayes-optimal case, this scaling can be found by analytically extremizing the loss under the compute constraint. Results based on random matrix theory support the same scaling in one-pass stochastic gradient descent (SGD) dynamics. 2 Frameworks that predict αN = αD = α necessarily imply N ∗ (C) ≍ D∗ (C) ≍ C 1/2 under the fixed-compute constraint C ∝ N D. However, a small asymmetry between αN and αD can still yield compute-optimal exponents close to 21 . Thus, near-square-root compute-optimal scaling does not by itself imply symmetric scaling. 3 Bordelon et al. [2025] also report scaling asymmetry, but there the relevant regimes are controlled by task difficulty rather than by the sparse-feature mechanism studied here. 4 We note that in training scenarios such as few-pass SGD, the underlying sparsity may be partially or entirely obscured. This disparity suggests a novel phase transition as the number of epochs increases.
2
2025]. In contrast, we focus on sparse feature activation, and show it can intrinsically yield different exponents in the model-limited versus data-limited regimes even in Bayes-optimal learning. We thus isolate rare feature coverage, rather than optimizer dynamics or feature learning, as a distinct source of scaling asymmetry. See Appendix A for more discussion of related works.
2
Statement of Problem
In this section, we motivate the model and outline the main goals of our analysis.
2.1
Learning under Sparse High-Dimensional Data
We first specify the data-generation process that yields a sparse power-law structure in the inputs, then describe the random feature model used to learn from these high-dimensional inputs. 2.1.1
Data Generation Process: Bernoulli-Random Activations with Power-Law Covariance
We consider input data X = [ x1 , . . . , xD ] ∈ RM ×D , with each xd ∈ RM , where D is the number of training examples (data size) and M is the data dimension. We henceforth drop the d subscript and simply refer to a single training example as x. We consider a structured data-generation process in which each input coordinate of x ∈ RM is randomly activated according to a heavy-tailed sparsity pattern. Specifically, for j = 1, 2, . . . , M , the jth coordinate of x is drawn as P(xj = 0) = 1 − j −α1 −1 , P xj = ±j −(α2 +1)/2 = 12 j −α1 −1 . (1) Hence E[xj ] = 0 and Var(xj ) = j −α1 −α2 −2 . Let pj := P(xj ̸= 0) = j −α1 −1 . We require α1 ≥ −1 to ensure that pj ≤ 1. Additionally, to ensure that the target variance is finite, the sum over the P∞ −α1 −α 2 −2 must converge, which holds iff α + α + 1 > 0. coordinate-wise variances j=1 j 1 2 This formulation introduces sparsity through Bernoulli input activations: most coordinates of x are zero, while the few active ones follow a heavy-tailed scaling controlled by α1 and α2 . 2.1.2
Random Feature Model
We study a linear model built from N ≪ M features randomly embedded from these inputs, with M ≫ max{N, D} so that truncation does not affect the power-law tail asymptotics. Each input x ∈ RM is mapped to features ϕ via a fixed random embedding u ∈ RN ×M , with uij ∼ N (0, 1/N ) and its N × M elements picked i.i.d., ϕ = ux, and the (per-example) prediction is ŷ(x) = θ ⊤ ϕ = θ ⊤ ux (and the batched prediction is ŷ = θ ⊤ Φ = θ ⊤ uX), with trainable θ ∈ RN .5 The target function P M is y(x) = w⊤ x = j=1 wj xj , for random w ∈ RM with wj i.i.d. and Var(wj ) = 1. The model parameters and their specifications are summarized in Table 1. This generative model allows us to examine how sparsity and heavy-tailed structure in the data influence performance. As we show below, these structural properties induce a loss scaling law with two distinct exponents, αN = ̸ αD , a sparsity-dependent compute-optimal frontier with an exponent 5
Our model is, in effect, a sketched linear regression. It is random-feature-like in the sense of a frozen random map followed by a trained linear readout, but it differs from canonical random-features models, e.g. Hastie et al. [2022], Gerace et al. [2020], Goldt et al. [2020], which study nonlinear feature maps in high-dimensional regimes. Extending those frameworks to the sparse, strongly non-Gaussian/heavy-tailed designs studied here would require going beyond the Gaussian-equivalence tools used in the nonlinear random-features literature Hu and Lu [2023], Mei and Montanari [2022], Goldt et al. [2022].
3
Table 1: Model parameters: Shapes and elementwise scaling. Parameter Description
Shape
Scaling
X u ϕ Φ θ w y y
M ×D N ×M N ×1 N ×D N ×1 M ×1 scalar 1×D
Per-coordinate std dev. j −(α1 +α2 +2)/2 uij ∼ N (0, 1/N ), i.i.d.
Training data (columns are examples) Fixed random embedding Features for a single input Features for the full dataset Trainable readout Target linear functional Per-example target function value Per-batch target function
O(1) O(1) O(1)
2
αC that shifts the optimal allocation toward data, and a regime in which GD fails to converge. The resulting regimes are summarized in Figure 1. This scaling asymmetry also persists in experiments with nonlinear MLPs. Interpretation: High-Dimensional RepPhase diagram in ( 1, 2) 2 resentations with Rare but Informative Features. Sparse activations of x can be understood as rare but highly informative fea1 One-exponent Scaling Two-exponent scaling tures that appear only in a small fraction 1 1 0 1>0 of examples yet carry disproportionate pre0 dictive value, e.g., an extremely infrequent Undefined keyword such as “pheochromocytoma,”, rare 1< 1 financial shocks such as the 2008 crisis, or an 1 uncommon diagnostic biomarker. The random feature model can be motiGD may fail Infinite variance 2< 1 2 1+ 2+1<0 vated by interpreting the high-dimensional vector x as a nonlinear expansion of a much lower-dimensional latent variable x̃ [Rahimi 3 3 2 1 0 1 2 and Recht, 2007]. The components of x could 1 correspond to all polynomial or Fourier baFigure 1: Phase diagram in (α1 , α2 ). Solid black sis functions of x̃, while the features ϕ = ux lines mark the boundary beyond which the model represent the smaller subset of directions acis well defined. Within it we identify three regimes: tually captured by our model. Although the symmetric one-exponent scaling, asymmetric twomapping through u is entirely linear, introexponent scaling, and a GD-failure regime with its ducing randomness, it can be interpreted as own scaling law. the first layer of a two-layer network (with no nonlinearity) whose weights are frozen at random initialization [Cho and Saul, 2009]. In this view, the model explores how much predictive structure can be recovered when learning is restricted to a fixed, randomly oriented subspace of the high-dimensional representation. Thus, this random feature model with sparse activations reveals how rare, high-impact signals shape optimization and give rise to distinct scaling laws, highlighting how sparsity fundamentally alters achievable performance.
4
2.2
Main Goals
Our goal is to characterize the trained-to-completion loss for the random sparse feature model described above. The population loss of the trained estimator is (Appendix B) h 2 i ℓBayes = Ex (y − ŷ(x))2 = Ex w⊤ x − θ ⋆⊤ ux , (2) where θ ⋆ denotes the minimum-ℓ2 -norm empirical-risk minimizer, equivalently the solution reached by GD from zero initialization when the empirical loss is trained to completion. Note that ŷ(x) = θ ⋆⊤ ux corresponds to a linear function of x whose weight vector lies in the row space of u. Consequently, when w∈ / rowspan(u), this feature map cannot recover y = w⊤ x exactly, and a nonzero approximation error is unavoidable. We aim to derive the scaling laws and compute-optimal frontier for this Bayes-optimal loss, emphasizing the dependence on the data distribution and its effect on optimization. In particular, we show that a two-exponent loss scaling law and a sparsity-dependent compute-optimal exponent can arise purely from the structure of the input distribution, in which x contains sparse but informative coordinates. Our objective is to identify exponents αN , αD , and αC such that ℓBayes ≍ N −αN + D−αD , and the compute-optimal frontier satisfies ℓ⋆Bayes ≍ C −αC , where ℓ⋆Bayes is the loss achieved by compute-optimal allocation. We also study GD in this setting, examining how optimization dynamics interact with the random feature map and the sparse data distribution, and whether it attains the Bayes rates or instead exhibits additional computational constraints, which we dub a scaling law of failure of GD. All proofs are deferred to the appendices.
3
Scaling Asymptotics of Sparse Random Features
Our strategy is to analyze the scaling of the Bayes-optimal loss, Eq. (2), for our model of sparse random features (specified in Section 2.1) in the underparameterized regime, where N is much smaller than D and in the overparameterized regime, where N vastly exceeds D. We will show that these regimes are governed by two distinct power laws, with different exponents, demonstrating that the loss is asymmetric and thus follows a two-exponent scaling law. This symmetry is broken precisely due to the sparsity of the input activation pattern. See Appendix C for details.
3.1
Loss Scaling from Unmodeled Features
The random feature model reduces to linear regression in the effective weight vector ŵ = u⊤ θ. The ⊤ ⊤ 2 population mean-squared error can be expressed as ℓ(ŵ) = Ex (w x − ŵ x) = (w − ŵ)⊤ Σx (w − ŵ), where Σx = E[xx⊤ ] is diagonal since the input coordinates are independent. To expose the scaling mechanism, consider a predictor that matches the first k coordinates of w and sets all remaining coordinates to zero: ŵj = wj 1{j≤k} . Then w − ŵ has entries 0 for j ≤ k and wj for P j > k, and substituting into the expression for ℓ(ŵ) yields ℓ(ŵ) = j>k wj2 Var(xj ). Under our assumptions wj2P = O(1), the population loss is therefore governed, up to constants, by the unmodeled input variance j>k Var(xj ). This calculation suggests that the asymptotic behavior of the loss is governed by the tail of the variance spectrum of the input distribution: X ℓBayes ≍ Var(xj ). (3) j>k
The argument above is only an intuition-building asymptotic device to make the mechanism transparent. Appendix C, particularly Sections 3.1 and C.2.2, provide a rigorous proof. 5
3.2
Prior Results: Scaling Law for Uniformly Activated Data with Power-Law Covariance
Before analyzing our model with sparse activations, we first review the scaling laws established for random feature models with fully active (non-sparse) input coordinates [Maloney et al., 2022]. These results form the baseline against which the effects of sparsity will be contrasted. Uniformly Activated Data with Power-Law Covariance. In the original formulation of Maloney et al. [2022], every input coordinate is active and is modeled as an independent Gaussian variable with variance decaying as a power law: xj ∼ N (0, j −α−1 ), j = 1, . . . , M, α > 0. This matches our model’s fully activated covariance structure under the choice α1 = −1 and α2 = α. Symmetry between Underparameterized and Overparameterized Regimes. For this non-sparse data structure, the optimal-loss scaling law was analyzed in detail in [Maloney et al., 2022]. We produce a simple intuitive argument showing why, in this setting, a single exponent governs the asymptotic behavior of the loss in both the underand overparameterized regimes, 1 1 α yielding a symmetric one-exponent law ℓBayes ∝ N + D with αN = αD = α. As explained above, the random feature model reduces to linear regression in the effective weight vector. Thus, if a learned predictor resolves the leading k ≪ M coefficients of w, the residual (unmodeled) variance R∞ P is j>k j −α−1 ≍ k j −α−1 dj = α1 k −α . With finitely many parameters N and finitely many samples D, the model can identify at most O(min{N, D}) coefficients, since this minimum represents the true bottleneck. This yields the scaling ℓBayes ∼ min{N, D}−α , so that the exponents satisfy αN = αD = α. In the sparse model of Section 2.1, one still gets a symmetric one-exponent scaling law when −1 ≤ α1 ≤ 0. In this regime, with high probability the low-index coordinates are activated sufficiently often across D samples that sparsity is effectively subdominant. The scaling is therefore still controlled by min{N, D}, as in the dense case, yielding αN = αD = α1 + α2 + 1 (Appendix C.1).
3.3
Two-Exponent Scaling Law for Sparse Activations (α1 > 0)
We now analyze our model of sparse activations for α1 > 0 (Appendix C.2). We begin by analyzing the asymptotic behavior of the loss in the underparameterized regime, where we will show that sparsity is effectively hidden. We then turn to the overparameterized regime, where sparsity becomes consequential and leads to a distinct asymptotic loss behavior. We work in the joint limit N, D → ∞, and we consider asymptotic regimes corresponding to different relative scaling of N and D.6 3.3.1
Underparameterized Regime (N ≪ D).
In this regime, the model fits only the leading high-variance coordinates of w, and the resulting scaling law coincides with that of the non-sparse case. Proposition 3.1 (Underparameterized Scaling). Under the sparse activation model described in Section 2.1, the Bayes-optimal loss in the underparameterized regime (N ≪ D) satisfies ℓBayes,N ∼ N −(α1 +α2 +1) ,
so that αN = α1 + α2 + 1.
(4)
Interpretation. The model captures only the leading high-variance coordinates of w. Since sparse coordinates are too rarely observed to influence estimation at this scale, they do not affect the loss. Consequently, the exponent αN = α1 + α2 + 1 coincides with the non-sparse case, so sparsity is effectively masked in the underparameterized regime. 6
In the sparse activation model, the relevant comparison is ultimately between N and the number of input coordinates ever observed in the dataset, K(D) (quantified below), which grows sublinearly in D when α1 > 0.
6
3.3.2
Overparameterized Regime (D ≪ N ).
A key quantity in this regime is the number of input coordinates that are ever observed (i.e., activated at least once). We therefore begin by quantifying the effective number of learnable coordinates as a function of D, which will set the correct scale separating the two asymptotic regimes. Lemma 3.2 (Learnable Coordinates in the Sparse Model when α1 > 0). Under the sparse activation distribution pj = P(xj ̸= 0) = j −α1 −1 , the expected number of coordinates that are active in at least one of the D training samples is asymptotic to the scale 1 1 D α1 +1 , 7 K(D) = Γ 1 − (5) α1 + 1 where Γ(·) is the Gamma function. Moreover, K(D) sets the scale of the number of coefficients of w that can be estimated from data: a coordinate that is never activated carries no information about its associated weight, and hence cannot contribute to the learned predictor. Theorem 3.3 (Overparameterized Scaling). Under the sparse activation model with α1 > 0 and coordinate variances Var(xj ) = j −α1 −α2 −2 , the Bayes-optimal loss in the overparameterized regime satisfies α +α +1 α1 + α 2 + 1 − 1 2 ℓBayes,D ≍ D α1 +1 , so that αD = . (6) α1 + 1 A rigorous proof of Theorem 3.3 is given in Appendix C.2.2 via a continuous argument. Remark 3.4 (Intuition). A simple heuristic recovers the same sample-limited exponent as in Theorem 3.3. Coordinate j is active in a single sample with probability pj = j −α1 −1 , so across D samples the expected number of activations is Dpj . A coordinate is observed often enough to be estimable once Dpj ≳ 1, which yields the cutoff j ≲ jD ≡ D1/(α1 +1) . Thus the number of learnable coordinates scales as K(D) ≍ jD , in agreement with Lemma 3.2. Coordinates with j > jD are effectively unobserved. Their contribution to the test error is the unmodeled variance α1 +α2 +1 α +α2 +1 P − 1α +1 −α1 −α2 −2 ≍ j −(α1 +α2 +1) = D − α1 +1 . Consequently, ℓ 1 , with Bayes,D ≍ D j>jD j D α1 +α2 +1 αD = α1 +1 , matching the rigorous sample-limited exponent in Theorem 3.3. Corollary 3.5 (Asymmetry of Exponents). For all α1 > 0 and α1 + α2 + 1 > 0, the exponents in the underparameterized and overparameterized regimes satisfy αD =
α1 + α 2 + 1 < α1 + α2 + 1 = αN . α1 + 1
(7)
Thus, for α1 > 0, the sparse model exhibits an intrinsic asymmetry between the sample-limited and parameter-limited regimes, giving rise to the two-exponent scaling law. Together, Proposition 3.1 and Theorem 3.3 yield: ℓBayes (N, D) ≍ N −(α1 +α2 +1) + D−(α1 +α2 +1)/(α1 +1) .
(8)
This expression makes the asymmetry explicit: model-limited error decays sharply, whereas datalimited error decays slowly. Since the mechanism depends on sparsity and variance profiles, we expect it to extend beyond Bernoulli masks to broader sparse distributions. 7
When α1 > 0, K(D) ≍ D1/(α1 +1) = o(D), i.e. only a sublinear number of coordinates are ever observed. Interpolation is still possible since the data lie in a K(D)-dimensional subspace, which a linear predictor in RN can fit. Unobserved coordinates receive no signal and are typically suppressed by mild regularization.
7
4
Empirical Study of the Scaling Laws
We now empirically validate the theoretical predictions of the previous section by evaluating the population loss over a wide range of D and N in the sparse random-feature model of Section 2.1. Computational details are provided in Appendices B and F.
4.1
Scaling Collapse and Universality
Bayes N 1 + 2 + 1
In many models, including ours and oth105 ers such as the Chinchilla scaling law [Hoffmann et al., 2022], the test loss exhibits dis104 tinct scaling in two regimes: ℓ ∼ D−αD −α N in the data-limited regime and ℓ ∼ N in the parameter-limited regime. A natu103 ral question is how to interpolate between N these two limits. A common ansatz is that 10 102 25 the loss takes the form ℓ = L̄(D−αD , N −αN ) 63 for some scale-invariant function L̄ satisfying 158 101 398 L̄(cL1 , cL2 ) = cL̄(L1 , L2 ). Simple choices in1000 clude L̄(L1 , L2 ) = max(L1 , L2 ) or L1 +L2 , but 0 10 more refined analyses [Maloney et al., 2022, 10 1 101 10 2 crit 100 D 1/(1 + 1)/N Zhang, 2024] show that L̄ often exhibits a peak near L1 ≈ L2 , a phenomenon known Figure 2: Scaling collapse and double descent. as double descent [Belkin et al., 2019, Hastie We plot the rescaled loss ℓ αN as a function of Bayes · N 8 et al., 2022]. Using the scaleinvariance of (DαD /N αN )1/αN , across a wide range of values for D, D−αD α N L̄, we may write ℓ · N = L̄ N −αN , 1 , so N , for fixed (α1 = 1.0, α2 = 0.3). All curves collapse α N the rescaled loss ℓ · N becomes a universal onto a single universal function Sα1 ,α2 , verifying a D αD universal scale-invariant structure consistent with the function of the compute ratio Ξ := N αN . We therefore expect the rescaled loss curves from predicted scaling law. A sharp peak emerges near different (D, N ) to collapse onto a single curve ξcrit consistent with the double descent phenomenon Sα1 ,α2 (Ξ): a dimensionless scaling function of at the interpolation threshold. the rescaled compute ratio. Empirical Scaling Collapse. In Figure 2, we find that the full loss curve across a wide range of (D, N ) pairs collapses onto a single universal function after appropriate rescaling: ℓBayes · 1
−α
D D N α1 +α2 +1 = Sα1 ,α2 N1 D α1 +1 , where S(u) = L̄(u−αN , 1), so that L̄( N −αN , 1) = Sα1 ,α2 (ξ), with 1/(α +1) DαD 1/αN ξ := N = D N1 . This implies that knowing the loss curve for one setting of (D, N ) αN suffices to predict its shape at other scales via a simple rescaling transformation. This collapse supports the two-regime scaling theory and recovers the predicted asymptotic behaviors: in the underparameterized regime, ℓ ∼ N −αN with αN = α1 + α2 + 1, while in the overparameterized regime, ℓ ∼ D−αD with αD = (α1 + α2 + 1)/(α1 + 1).
Double Descent Peak. In the underparameterized limit, the loss scales as ℓ ∼ N −(α1 +α2 +1) , while in the overparameterized limit it scales as ℓ ∼ D−(α1 +α2 +1)/(α1 +1) . This suggests a crossover at N ∼ D1/(α1 +1) , which coincides with the point at which the number of model parameters N 8
The double-descent peak is a feature of benign interpolation under ridgeless fitting, rather than of the scaling-law picture itself. Regularization or early stopping smooths it out [Hastie et al., 2022].
8
matches the number of activated (and hence learnable) coordinates K(D) from Lemma 3.2. The shape of Sα1 ,α2 (ξ) in Figure 2 reveals a prominent double descent phenomenon near: 1
ξcrit = Γ
α1 1+α1
,
(9)
corresponding to the threshold where K(D) matches N , with location determined universally by the sparsity exponent α1 .
5
Compute-Optimal Scaling Laws
We now analyze how to allocate data and model size under a fixed compute budget so as to minimize the loss. We adopt the compute model C = N D · min{N, D} as a unified proxy for the cost of training to completion in our setting; see Appendix D for discussion. We then minimize ℓ(N, D) ∼ N −αN + D−αD subject to fixed compute C to derive the compute-optimal scaling laws. Proposition 5.1 (Compute-Optimal Frontier). Let the Bayes-optimal loss, for α1 > 0, scale as 2 +1 ℓ(N, D) ∼ N −αN + D−αD , with αN = α1 + α2 + 1 and αD = α1α+α . Under a fixed compute budget 1 +1 C = N D · min{N, D}, the unique compute-optimal allocation lies in the underparameterized regime N < D, and the optimal allocation and resulting loss scale as ℓ∗ (C) ∼ C −αC , N ∗ (C) ∼ C αD /(αN +2αD ) = C 1/(α1 +3) , D∗ (C) ∼ C 1−2αD /(αN +2αD ) = C (α1 +1)/(α1 +3) , αN αD α1 + α 2 + 1 = > 0. αC := αN + 2αD α1 + 3
(10)
Remark 5.2 (Absence of Overparameterized Optimum). The compute-optimal solution always lies in the underparameterized regime: N ∗ (C) < D∗ (C) for all C. An overparameterized optimum N > D, with C = N D2 , yields N ∗ (C)/D∗ (C) ∼ C (αD −αN )/(2αN +αD ) ≪ 1, contradicting the assumed regime N > D. Hence the only valid optimum satisfies N < D. Remark 5.3 (Impact of Sparsity on Optimal Scaling). As α1 increases, individual samples become more sparse, carrying less information per example. At the same time, increasing α1 also steepens the marginal variance spectrum Var(xj ) = j −α1 −α2 −2 , so the residual tail loss decays more rapidly once coordinates are resolved. Consequently, for fixed α2 < 2, the compute-optimal exponent 2 +1 αC = α1α+α increases monotonically with α1 , approaching 1 as α1 → ∞. As α1 increases, the 1 +3 optimal model size, N ∗ (C), grows more slowly with compute, while the optimal dataset size, D∗ (C), grows more rapidly. As a result, the allocation shifts toward spending a larger fraction of compute on data rather than parameters. Empirical Validation. To validate Proposition 5.1, Figure 3 shows test loss versus total compute C for different N . Each curve corresponds to a fixed N , while C is varied by sweeping D. The dashed line shows the predicted compute-optimal frontier, ℓ∗ (C) ∼ C −αC . As C increases, the empirical curves approach this frontier, confirming the predicted scaling.
6
Gradient Descent Training Dynamics
9
We next study the training dynamics of the readout weights θ under full-batch GD, a broadly meaningful computational model with some generality beyond the toy setting, on the empirical squared loss ℓ(θ) = 2 ⊤ 1 ⊤ 2D w X − θ uX 2 . With step size η, one η GD step is ∆θ = D uXX⊤ w − u⊤ θ . It is convenient to rewrite the dynamics in terms of the input-space residual rt := w − ŵt , where ŵt := u⊤ θ t ∈ RM . Then η rt+1 = IM − u⊤ u XX⊤ rt , (11) D
N
10 15 25 39 63 100 158 251 398 630 1000 Theory
102 101 100 10 1 10 2 10 3 10 4 10 5 103
104
105
106
C
107
108
109
1010
so convergence is controlled by the spectrum Figure 3: Empirical scaling of loss with compute. Each curve corresponds to test loss ℓ versus total of u⊤ u XX⊤9 ; see Appendix E. compute C for a fixed model size N , with C varied 6.1 Convergence with High Proba- by sweeping D, for fixed (α1 = 1.0, α2 = 0.3). The dashed line denotes the predicted compute-optimal bility scaling ℓ∗ (C) ∼ C −αC . The empirical loss converges A first question is how η must scale for (11) to toward this frontier at high compute, just past the converge. Since ℓ(θ) is quadratic, full-batch double descent spike for each curve demonstrating GD is a linear iteration, and a standard stabil- agreement between theory and experiment. ity criterion implies that (11) converges when2 ever 0 < η < λ , where λmax (·) denotes the largest eigenvalue [Nesterov, 2004]. 1 ⊤ ⊤ max( D u u XX ) 1 Equivalently, this is the usual condition 0 < η < 2/λmax (∇2θ ℓ), since ∇2θ ℓ(θ) = D uXX⊤ u⊤ .10 In the underparameterized regime, a standard concentration argument for the feature-space Hessian 1 ΦΦ⊤ yields λmax (∇2θ ℓ) = 1 + oP (1) and hence any fixed η ∈ (0, 2) is admissible with ∇2θ ℓ(θ) = D high probability. The overparameterized case, however, is more subtle and is controlled by whether the random embedding acts isometrically on the data span. Proposition 6.1 (Step-size Stability with High Probability). Assume the sparse activation model of Section 2.1 with α1 ≥ −1 and α1 + α2 + 1 > 0. Assume further that the random feature map is sufficiently wide that it is nearly isometric on the (low-dimensional) span of the data,11 so that 1 ⊤ ⊤ λmax (12) u u XX = 1 + oP (1). D Consequently, for any fixed η ∈ (0, 2) the iteration (11) is stable and converges with high probability. The data dependence enters through the empirical covariance XX⊤ . For dense inputs, rank(XX⊤ ) = rank(X) ≈ min{M, D}, which is typically ≈ D when D ≲ M , whereas under our sparse model, rank(XX⊤ ) ≤ K(D) ≍ D1/(α1 +1) , so it grows only sublinearly in D. 10 Here u⊤ u and XX⊤ are positive semidefinite. Moreover, u⊤ u XX⊤ has the same nonzero eigenvalues as the symmetric positive semidefinite matrix (XX⊤ )1/2 (u⊤ u)(XX⊤ )1/2 , so its spectrum on the data span is real and nonnegative. 11 Concretely, since rank(XX⊤ ) ≤ D, standard random projection heuristics suggest that for N ≫ D the operator ⊤ u u acts approximately like IM when sandwiched against XX⊤ . 9
10
6.2
A Scaling Law for Failure of Gradient Descent
Proposition 6.1 is a high-probability statement; rare datasets can still contain a single rare activation spike that produces an anomalously large top eigenvalue and causes divergence.12 GD becomes η ⊤ ⊤ unstable when λmax D u uXX > 2. Under the near-isometry approximation u⊤ u ≈ IM on the data span, this reduces to λmax XX⊤ > 2D η . This is only a concern when activated amplitudes grow with index, i.e. when α2 < −1 (since then x2j = j −(α2 +1) increases with j on activation). Theorem 6.2 (Failure Probability of GD). Fix η ∈ (0, 2) and suppose α2 < −1. Under the sparse 1 +α2 +1 activation model of Section 2.1, with α1 ≥ −1 and α1 + α2 + 1 > 0, let ν := α−α > 0. Then the 2 −1 probability that full-batch GD is unstable due to a rare spike on a random dataset of size D obeys Prare (diverge) ≍ D−ν ,
(13)
up to η-dependent constants. This theorem identifies instability under a fixed step size, not a fundamental optimization barrier. Reducing the learning rate or using clipped updates can restore stability [Nguyen et al., 2023].
Sparse 10 1 10 1 10 2
10 2 10 3 10 4
7 Experiments with Nonlinear Activations
10 3 101
N
101
102
102
D
103
Dense 10 1
10 2 A natural question is whether the twoexponent asymmetry we identify is specific 10 2 10 3 to the linear setting or reflects a more general property of sparse data distributions. In Fig10 4 10 3 ure 4, we test robustness to nonlinearity by 101 101 102 102 103 N D training a two-layer network with hidden layer ϕ = σ(ux) (ReLU activation, frozen random Figure 4: Asymmetry persists under nonlinfirst-layer weights) and a linear readout, using earity. Test loss scaling under the ReLU feature full-batch accelerated GD with Nesterov momap ϕ = σ(ux), trained with Nesterov + adaptive mentum and adaptive restart; see Appendix F restart; (α1 , α2 ) = (1.0, 0.3), 5 seeds. Top: sparse; for experimental details. The fitted exponents bottom: dense. Left: N -sweep at D = 50,000; right: shift downward in the model-limited regime, D-sweep at N = 8000. A single exponent α ≈ 1.5 consistent with the general expectation that (dashed lines, jointly fitted with separate intercepts) nonlinear feature maps smooth power-law indescribes the sparse N -, dense N -, and dense Dput spectra [Cho and Saul, 2009, Louart et al., sweeps. The sparse D-sweep requires its own shal2018, Bietti and Mairal, 2019]. Crucially, the lower fit, αD ≈ 1.2, breaking the dense-baseline symqualitative asymmetry persists: sparse N -, metry αN = αD predicted by linear theory. dense N -, dense D-sweeps are well-described by a single shared exponent, while the sparse D-sweep is better described by a distinct exponent. This supports our main result: the bottleneck mechanism—unobserved coordinates carrying no signal—is robust to nonlinearity and first-order optimization.
Rare activations can create spectral outliers in XX⊤ that are not controlled by the high-probability spectral concentration arguments used in recent approximation/bias/variance analyses of linear-regression scaling laws [Lin et al., 2024, 2025, Yan et al., 2026], so we instead adopt a probabilistic treatment. 12
11
8
Conclusion
We introduced a model that gives rise to two distinct scaling exponents under sparse input activations. By analyzing the population loss in both under- and overparameterized regimes, we showed how sparsity breaks the symmetry of scaling laws in the random feature setting. We derived a computeoptimal scaling law showing that sparsity enhances efficiency by accelerating the decay of test loss with compute, with a frontier that prioritizes data over model size. Analyzing GD dynamics, we proved high-probability stability for all step sizes below 2 and identified a suppressed sublinear scaling law for its probability of failure in the overparameterized regime, suggesting a role for gradient clipping in practical optimization. We validated these results experimentally in the theoretical setting and extended them to a nonlinear feature map, confirming robustness of the sparsity-induced asymmetry to nonlinearity. Our main takeaway is that sparsity as a structural property fundamentally shapes scaling behavior in high-dimensional learning. Limitations. Our results are built around data sparsity as the organizing principle, yet several extensions remain for future work. Incorporating exogenous label noise, which we expect to impose an error floor rather than remove the sparsity-driven bottleneck until it dominates the optimization scale, and analyzing ridge regularization, which would refine the statistical picture in the presence of explicit shrinkage, would extend our picture to noisy and explicitly regularized settings. More ambitiously, a theoretical treatment of nonlinear feature maps—complementing our empirical validation—and an extension to the dynamical regime of few- to multi-pass SGD both confront a common technical obstacle: the non-Gaussian, heavy-tailed nature of our features requires extending existing deterministic-equivalent tools beyond the Gaussian setting [Mei and Montanari, 2022, Hu and Lu, 2023]. Relatedly, extending our analysis to the feature-learning regime is an important direction, especially since feature learning may either concentrate capacity on rare informative coordinates [Bordelon et al., 2025] or induce sparse representations that hurt generalization in smooth tasks [Petrini et al., 2022]. A further natural direction is to test these effects at scale, in deeper networks and in transformers, particularly in light of sparse-feature observations from sparse autoencoders [Bricken et al., 2023, Templeton et al., 2024] and recent connections between superposition [Elhage et al., 2022] and scaling laws [Liu et al., 2025].
References Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprint arXiv:2001.08361, 2020. Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, Tom Hennigan, Eric Noland, Katie Millican, George van den Driessche, Bogdan Damoc, Aurelia Guy, Simon Osindero, Karen Simonyan, Erich Elsen, Jack W. Rae, Oriol Vinyals, and Laurent Sifre. Training compute-optimal large language models. In Advances in Neural Information Processing Systems (NeurIPS 2022), volume 35, pages 30016–30030, 2022. Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan. Spectrum dependent learning curves in kernel regression and wide neural networks. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 1024–1034, 2020.
12
Stefano Spigler, Mario Geiger, and Matthieu Wyart. Asymptotic learning curves of kernel methods: Empirical data versus teacher–student paradigm. Journal of Statistical Mechanics: Theory and Experiment, 2020(12):124001, 2020. doi: 10.1088/1742-5468/abc61d. Yasaman Bahri, Ethan Dyer, Jared Kaplan, Jaehoon Lee, and Utkarsh Sharma. Explaining neural scaling laws. arXiv preprint arXiv:2102.06701, 2021. Andrew Maloney, Daniel A. Roberts, and James Sully. A solvable model of neural scaling laws. arXiv preprint arXiv:2210.16859, 2022. Blake Bordelon, Alexander Atanasov, and Cengiz Pehlevan. A dynamical model of neural scaling laws. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 4345–4382, 2024. Elliot Paquette, Courtney Paquette, Lechao Xiao, and Jeffrey Pennington. 4+3 phases of computeoptimal neural scaling laws. In Advances in Neural Information Processing Systems (NeurIPS 2024), volume 37, pages 16459–16537, 2024. Blake Bordelon, Alexander Atanasov, and Cengiz Pehlevan. How feature learning can improve neural scaling laws. Journal of Statistical Mechanics: Theory and Experiment, 2025(8):084002, 2025. doi: 10.1088/1742-5468/adefb1. Mario Geiger, Stefano Spigler, Stéphane d’Ascoli, Levent Sagun, Marco Baity-Jesi, Giulio Biroli, and Matthieu Wyart. Jamming transition as a paradigm to understand the loss landscape of deep neural networks. Physical Review E, 100(1):012115, 2019. doi: 10.1103/PhysRevE.100.012115. Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019. doi: 10.1073/pnas.1903070116. Mikhail Belkin, Daniel Hsu, and Ji Xu. Two models of double descent for weak features. SIAM Journal on Mathematics of Data Science, 2(4):1167–1180, 2020. doi: 10.1137/20M1336072. Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani. Surprises in highdimensional ridgeless least squares interpolation. The Annals of Statistics, 50(2):949–986, 2022. doi: 10.1214/21-AOS2133. Peter L. Bartlett, Philip M. Long, Gábor Lugosi, and Alexander Tsigler. Benign overfitting in linear regression. Proceedings of the National Academy of Sciences, 117(48):30063–30070, 2020. doi: 10.1073/pnas.1907378117. Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75 (4):667–766, 2022. doi: 10.1002/cpa.22008. Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever. Deep double descent: Where bigger models and more data hurt. In International Conference on Learning Representations (ICLR), 2020. doi: 10.48550/arXiv.1912.02292. Leonardo Defilippis, Bruno Loureiro, and Theodor Misiakiewicz. Dimension-free deterministic equivalents and scaling laws for random feature regression. In Advances in Neural Information Processing Systems (NeurIPS 2024), volume 37, pages 104630–104693, 2024. 13
Licong Lin, Jingfeng Wu, Sham M. Kakade, Peter L. Bartlett, and Jason D. Lee. Scaling laws in linear regression: Compute, parameters, and data. In Advances in Neural Information Processing Systems (NeurIPS 2024), volume 37, pages 60556–60606, 2024. Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mezard, and Lenka Zdeborova. Generalisation error in learning with random features and the hidden manifold model. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 3452–3462, 2020. Sebastian Goldt, Marc Mezard, Florent Krzakala, and Lenka Zdeborova. Modelling the influence of data structure on learning in neural networks: the hidden manifold model. Physical Review X, 10 (4):041044, 2020. doi: 10.1103/PhysRevX.10.041044. Hong Hu and Yue M. Lu. Universality laws for high-dimensional learning with random features. IEEE Transactions on Information Theory, 69(3):1932–1964, 2023. doi: 10.1109/TIT.2022.3217698. Sebastian Goldt, Bruno Loureiro, Galen Reeves, Florent Krzakala, Marc Mézard, and Lenka Zdeborová. The Gaussian equivalence of generative models for learning with shallow neural networks. In Proceedings of The 33rd International Conference on Algorithmic Learning Theory, volume 167 of Proceedings of Machine Learning Research, pages 426–471, 2022. Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems (NeurIPS 2007), volume 20, pages 1177–1184, 2007. Youngmin Cho and Lawrence K. Saul. Kernel methods for deep learning. In Advances in Neural Information Processing Systems 22 (NeurIPS 2009), volume 22, pages 342–350, 2009. Zhengkang Zhang. Neural scaling laws from large-N field theory: Solvable model beyond the ridgeless limit. arXiv preprint arXiv:2405.19398, 2024. Yurii Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Kluwer Academic Publishers, 2004. Licong Lin, Jingfeng Wu, and Peter L. Bartlett. Improved scaling laws in linear regression via data reuse. In Advances in Neural Information Processing Systems (NeurIPS 2025), 2025. Tingkai Yan, Haodong Wen, Binghui Li, Kairong Luo, Wenguang Chen, and Kaifeng Lyu. Larger datasets can be repeated more: A theoretical analysis of multi-epoch scaling in linear regression. In The Fourteenth International Conference on Learning Representations, 2026. Ta Duy Nguyen, Thien H. Nguyen, Alina Ene, and Huy Nguyen. Improved convergence in high probability of clipped gradient methods with heavy tailed noise. In Advances in Neural Information Processing Systems (NeurIPS 2023), volume 36, pages 24191–24222, 2023. Cosme Louart, Zhenyu Liao, and Romain Couillet. A random matrix approach to neural networks. The Annals of Applied Probability, 28(2):1190–1248, 2018. doi: 10.1214/17-AAP1328. Alberto Bietti and Julien Mairal. On the inductive bias of neural tangent kernels. In Advances in Neural Information Processing Systems (NeurIPS 2019), volume 32, pages 12873–12884, 2019. Leonardo Petrini, Francesco Cagnetta, Eric Vanden-Eijnden, and Matthieu Wyart. Learning sparse features can lead to overfitting in neural networks. In Advances in Neural Information Processing Systems (NeurIPS 2022), volume 35, pages 9403–9416, 2022. 14
Trenton Bricken, Adly Templeton, Joshua Batson, Brian Chen, Adam Jermyn, Tom Conerly, Nick Turner, Cem Anil, Carson Denison, Amanda Askell, Robert Lasenby, Yifan Wu, Shauna Kravec, Nicholas Schiefer, Tim Maxwell, Nicholas Joseph, Zac Hatfield-Dodds, Alex Tamkin, Karina Nguyen, Brayden McLean, Josiah E. Burke, Tristan Hume, Shan Carter, Tom Henighan, and Christopher Olah. Towards monosemanticity: Decomposing language models with dictionary learning. Transformer Circuits Thread, 2023. Adly Templeton, Tom Conerly, Jonathan Marcus, Jack Lindsey, Trenton Bricken, Brian Chen, Adam Pearce, Craig Citro, Emmanuel Ameisen, Adam Jermyn, Cem Anil, Carson Denison, Amanda Askell, Robert Lasenby, Yifan Wu, Shauna Kravec, Nicholas Schiefer, Tim Maxwell, Nicholas Joseph, Alex Tamkin, Karina Nguyen, Brayden McLean, Josiah E. Burke, Tristan Hume, Shan Carter, Tom Henighan, and Christopher Olah. Scaling monosemanticity: Extracting interpretable features from claude 3 sonnet. Transformer Circuits Thread, 2024. Nelson Elhage, Tristan Hume, Catherine Olsson, Nicholas Schiefer, Tom Henighan, Shauna Kravec, Zac Hatfield-Dodds, Robert Lasenby, Dawn Drain, Carol Chen, Roger Grosse, Sam McCandlish, Jared Kaplan, Dario Amodei, Martin Wattenberg, and Christopher Olah. Toy models of superposition. Transformer Circuits Thread, 2022. Yizhou Liu, Ziming Liu, and Jeff Gore. Superposition yields robust neural scaling. In Advances in Neural Information Processing Systems (NeurIPS 2025), 2025. Tom Henighan, Jared Kaplan, Mor Katz, Mark Chen, Christopher Hesse, Jacob Jackson, Heewoo Jun, Tom B. Brown, Prafulla Dhariwal, Scott Gray, Chris Hallacy, Benjamin Mann, Alec Radford, Aditya Ramesh, Nick Ryder, Daniel M. Ziegler, John Schulman, Dario Amodei, and Sam McCandlish. Scaling laws for autoregressive generative modeling. arXiv preprint arXiv:2010.14701, 2020. Joel Hestness, Sharan Narang, Newsha Ardalani, Gregory Diamos, Heewoo Jun, Hassan Kianinejad, Md. Mostofa Ali Patwary, Yang Yang, and Yanqi Zhou. Deep learning scaling is predictable, empirically. arXiv preprint arXiv:1712.00409, 2017. Xiaohua Zhai, Alexander Kolesnikov, Neil Houlsby, and Lucas Beyer. Scaling vision transformers. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12104–12113, 2022. Tomer Porian, Mitchell Wortsman, Jenia Jitsev, Ludwig Schmidt, and Yair Carmon. Resolving discrepancies in compute-optimal scaling of language models. In Advances in Neural Information Processing Systems (NeurIPS 2024), volume 37, 2024. Shane Bergsma, Nolan Dey, Gurpreet Gosal, Gavia Gray, Daria Soboleva, and Joel Hestness. Power lines: Scaling laws for weight decay and batch size in llm pre-training. In Advances in Neural Information Processing Systems (NeurIPS 2025), 2025. Tanishq Kumar, Zachary Ankner, Benjamin Frederick Spector, Blake Bordelon, Niklas Muennighoff, Mansheej Paul, Cengiz Pehlevan, Christopher Re, and Aditi Raghunathan. Scaling laws for precision. In The Thirteenth International Conference on Learning Representations, 2025. Nikhil Sardana, Jacob Portes, Sasha Doubov, and Jonathan Frankle. Beyond Chinchilla-optimal: Accounting for inference in language model scaling laws. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 43445–43460, 2024. 15
Peihao Wang, Rameswar Panda, and Zhangyang Wang. Data efficient neural scaling law via model reusing. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 36193–36204, 2023. Niklas Muennighoff, Alexander M. Rush, Boaz Barak, Teven Le Scao, Nouamane Tazi, Aleksandra Piktus, Sampo Pyysalo, Thomas Wolf, and Colin A. Raffel. Scaling data-constrained language models. In Advances in Neural Information Processing Systems (NeurIPS 2023), volume 36, pages 50358–50376, 2023. Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová. Optimal errors and phase transitions in high-dimensional generalized linear models. Proceedings of the National Academy of Sciences, 116(12):5451–5460, 2019. doi: 10.1073/pnas.1802705116. Benjamin Aubin, Antoine Maillard, Jean Barbier, Florent Krzakala, Nicolas Macris, and Lenka Zdeborová. The committee machine: Computational to statistical gaps in learning a two-layers neural network. Journal of Statistical Mechanics: Theory and Experiment, 2019(12):124023, 2019. doi: 10.1088/1742-5468/ab43d2. Song Mei, Theodor Misiakiewicz, and Andrea Montanari. Learning with invariances in random features and kernel models. In Proceedings of the Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 3351–3418, 2021. Song Mei, Theodor Misiakiewicz, and Andrea Montanari. Generalization error of random feature and kernel methods: Hypercontractivity and kernel matrix concentration. Applied and Computational Harmonic Analysis, 59:3–84, 2022. doi: 10.1016/j.acha.2021.12.003. Hugo Latourelle-Vigeant and Elliot Paquette. Dyson equation for correlated linearizations and test error of random features regression. arXiv preprint arXiv:2312.09194, 2023. Alexander Atanasov, Jacob A. Zavatone-Veth, and Cengiz Pehlevan. Scaling and renormalization in high-dimensional regression. arXiv preprint arXiv:2405.00592, 2024. Arie Wortsman and Bruno Loureiro. Kernel ridge regression under power-law data: Spectrum and generalization. arXiv preprint arXiv:2510.04780, 2025. Elliot Paquette, Ke Liang Xiao, and Yizhe Zhu. Power-law spectrum of the random feature model. arXiv preprint arXiv:2603.14578, 2026. Jean Barbier, Francesco Camilli, Minh-Toan Nguyen, Mauro Pastore, and Rudy Skerk. Statistical physics of deep learning: Optimal learning of a multi-layer perceptron near interpolation. arXiv preprint arXiv:2510.24616, 2025. Leonardo Defilippis, Yizhou Xu, Julius Girardin, Emanuele Troiani, Vittorio Erba, Lenka Zdeborová, Bruno Loureiro, and Florent Krzakala. Scaling laws and spectra of shallow neural networks in the feature learning regime. In The Fourteenth International Conference on Learning Representations, 2026. Dechen Zhang, Xuan Tang, Yingyu Liang, and Difan Zou. Scaling laws for precision in highdimensional linear regression. arXiv preprint arXiv:2602.19241, 2026. Francesco Cagnetta, Allan Raventós, Surya Ganguli, and Matthieu Wyart. Deriving neural scaling laws from the statistics of natural language. arXiv preprint arXiv:2602.07488, 2026. 16
Maissam Barkeshli, Alberto Alfarano, and Andrey Gromov. On the origin of neural scaling laws: From random graphs to natural language. arXiv preprint arXiv:2601.10684, 2026. Binghui Li, Fengling Chen, Zixun Huang, Lean Wang, and Lei Wu. Functional scaling laws in kernel regression: Loss dynamics and learning rate schedules. In Advances in Neural Information Processing Systems (NeurIPS 2025), 2025. Blake Bordelon and Francesco Mori. Theory of optimal learning rate schedules and scaling laws for a random feature model. arXiv preprint arXiv:2602.04774, 2026. Jakob Kramp, Javed Lindner, and Moritz Helias. Dynamics of neural scaling laws in random feature regression with powerlaw-distributed kernel eigenvalues. arXiv preprint arXiv:2602.23039, 2026. Yunwei Ren, Eshaan Nichani, Denny Wu, and Jason D. Lee. Emergence and scaling laws in SGD learning of shallow neural networks. In Advances in Neural Information Processing Systems (NeurIPS 2025), 2025. Gérard Ben Arous, Murat A. Erdogdu, N. Mert Vural, and Denny Wu. Learning quadratic neural networks in high dimensions: SGD dynamics and scaling laws. In Advances in Neural Information Processing Systems (NeurIPS 2025), 2025. Roman Worschech and Bernd Rosenow. Analyzing neural scaling laws in two-layer networks with power-law data spectra. In The Thirteenth International Conference on Learning Representations, 2025. Xavier Glorot, Antoine Bordes, and Yoshua Bengio. Deep sparse rectifier neural networks. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics (AISTATS), volume 15 of Proceedings of Machine Learning Research, pages 315–323, 2011. Bruno A. Olshausen and David J. Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 381(6583):607–609, 1996. doi: 10.1038/381607a0. Noam Shazeer, Azalia Mirhoseini, Krzysztof Maziarz, Andy Davis, Quoc V. Le, Geoffrey E. Hinton, and Jeff Dean. Outrageously large neural networks: The sparsely-gated mixture-of-experts layer. In The Fifth International Conference on Learning Representations, 2017. Jan Ludziejewski, Jakub Krajewski, Kamil Adamczewski, Maciej Pióro, Michał Krutul, Szymon Antoniak, Kamil Ciebiera, Krystian Król, Tomasz Odrzygóźdź, Piotr Sankowski, Marek Cygan, and Sebastian Jaszczur. Scaling laws for fine-grained mixture of experts. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 33270–33288, 2024. Samira Abnar, Harshay Shah, Dan Busbridge, Alaaeldin El-Nouby, Joshua M. Susskind, and Vimal Thilak. Parameters vs FLOPs: Scaling laws for optimal sparsity for mixture-of-experts language models. In Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Proceedings of Machine Learning Research, pages 204–230, 2025. Yurii Nesterov. A method for solving the convex programming problem with convergence rate o(1/k 2 ). Doklady Akademii Nauk SSSR, 269:543–547, 1983. Brendan O’Donoghue and Emmanuel Candès. Adaptive restart for accelerated gradient schemes. Foundations of Computational Mathematics, 15(3):715–732, 2015. doi: 10.1007/s10208-013-9150-3. 17
A
Additional Related Work
Empirical scaling laws relating test loss to parameters, data, and compute were systematically documented in language modeling [Kaplan et al., 2020] and refined into compute-optimal prescriptions emphasizing data-limited training in LLMs [Hoffmann et al., 2022]. Related scaling behavior has also been observed across tasks and modalities, including autoregressive generative modeling [Henighan et al., 2020], multi-domain studies [Hestness et al., 2017] and vision transformer scaling analyses [Zhai et al., 2022]. More recent empirical work has expanded scaling-law analyses to address computeoptimal discrepancies across experimental protocols [Porian et al., 2024], practical scaling laws for weight decay and batch size in LLM pretraining [Bergsma et al., 2025], precision-aware scaling under quantization and training precision [Kumar et al., 2025], inference-aware compute-optimal scaling [Sardana et al., 2024], data-scarcity regimes via model reuse [Wang et al., 2023], and data-constrained scaling under token repetition [Muennighoff et al., 2023]. On the theory side, a growing literature develops tractable models that reproduce and explain scaling-law phenomenology. Statistical-physics teacher-student models characterize optimal errors and phase transitions in high-dimensional generalized linear models [Barbier et al., 2019] and computational-to-statistical gaps in committee machines [Aubin et al., 2019]. Solvable neural-scaling models connect power-law learning curves to spectral structure and compute-optimal frontiers [Bahri et al., 2021, Maloney et al., 2022]. Random-feature and kernel analyses address invariances [Mei et al., 2021], hypercontractivity and kernel matrix concentration [Mei et al., 2022], correlated linearizations and test error [Latourelle-Vigeant and Paquette, 2023], and dimension-free deterministic equivalents [Defilippis et al., 2024]. Scaling and renormalization approaches provide a complementary high-dimensional regression perspective [Atanasov et al., 2024]. Recent work also extends this powerlaw spectral perspective to kernel ridge regression under anisotropic power-law data [Wortsman and Loureiro, 2025] and spectral inheritance through random-feature maps [Paquette et al., 2026]. Related statistical-physics analyses study feature learning and specialization in finite-width teacherstudent neural networks near interpolation [Barbier et al., 2025], as well as scaling-law phase diagrams in shallow feature-learning networks [Defilippis et al., 2026]. A complementary line studies compute-optimal scaling in linear regression [Lin et al., 2024] and precision-aware variants in highdimensional linear regression [Zhang et al., 2026]. Recent work also studies the origin of scaling laws in simplified language and sequence-modeling settings, including theories based on natural-language statistics [Cagnetta et al., 2026] and controlled graph/random-walk generative processes [Barkeshli et al., 2026]. Beyond Bayes-optimal analyses, several works study training dynamics. Dynamical models of one-pass or finite-time training exhibit multiple scaling phases and shift compute-optimal allocations [Bordelon et al., 2024, Paquette et al., 2024]. Feature learning can also change scaling exponents and compute-optimal behavior [Bordelon et al., 2025]. Data reuse and multi-epoch training shift effective exponents and optimal allocations [Lin et al., 2025, Yan et al., 2026]. Functional loss dynamics and learning-rate schedules in kernel or random-feature regression are studied in [Li et al., 2025, Bordelon and Mori, 2026], while early stopping in random-feature regression with power-law spectra is analyzed in [Kramp et al., 2026]. Related analyses of SGD in shallow or two-layer networks with heterogeneous teacher components [Ren et al., 2025], quadratic feature-learning models [Ben Arous et al., 2025], or power-law data spectra [Worschech and Rosenow, 2025] derive sample-, time-, and parameter-dependent scaling exponents. Our work is motivated by a different axis: sparse or conditional feature activation [Glorot et al., 2011, Olshausen and Field, 1996], in which only a subset of features is observed or used for any example. Such sparsity is common in modern representation learning and conditional computation (e.g., mixture-of-experts) [Shazeer et al., 2017]. Related empirical work studies scaling laws for sparse 18
conditional computation in mixture-of-experts models [Ludziejewski et al., 2024, Abnar et al., 2025]. Rather than refining existing spectral-decay accounts, we isolate how sparsity in what the dataset reveals can itself produce an intrinsic asymmetry between model-limited and data-limited scaling exponents.
B
Population Objective of Sparse Model
We begin by reformulating the objective in Eq. (2) as the following loss function: ℓ =
D 2 γ 1 X ∥θ∥22 + y(xa ) − θ ⊤ uxa D D a=1
=
=
D 1 X
γ ∥θ∥22 + D D 1 γ ∥θ∥22 + D D
w⊤ xa − θ ⊤ uxa
2
(14)
a=1 D h X
i ⊤ ⊤ ⊤ ⊤ ⊤ w ⊤ xa x⊤ w − 2 θ ux x w + θ ux x u θ , a a a a a
a=1
which we wish to minimize with respect to θ. While our focus is on the ridge-less limit, γ = 0, we retain it for mathematical convenience and take the limit γ → 0 when appropriate. Using Φ = uX, we find the optimal solution θ ⋆ satisfies: θ ⋆⊤ = yΦ⊤ q = yQΦ⊤ ,
(15)
where q = (γ IN + ΦΦ⊤ )−1 , Q = (γ ID + Φ⊤ Φ)−1 . We are interested in the difference between the effective weight vector weff , defined by ⊤ weff = θ ⋆⊤ u,
and the true target weights w⊤ . In particular, the generalization error is given by ⊤ Egen = weff − w⊤ Σx weff − w . Table 2 lists the shapes and elementwise scalings used in our computations.
C
Scaling Asymptotics of Sparse Random Features
C.1
−1 ≤ α1 ≤ 0: One-Exponent Scaling Law
In our sparse model of Section 2.1, we obtain a symmetric one-exponent scaling law when −1 ≤ α1 ≤ 0, similar to dense random feature models [Bahri et al., 2021, Maloney et al., 2022]. The key point is that generalization can be bottlenecked by three resources: the number of parameters N , the number of samples D, and the set of input coordinates that are actually observed (activated at least once) in the D training examples. Recall that the activation probability of coordinate j is pj := P(xj ̸= 0) = j −(α1 +1) . 19
Table 2: Shapes and elementwise scaling. Parameter
Description
Shape
Scaling (elementwise)
X
Training data (columns are examples)
M ×D
u θ w y Φ Q q
Fixed random embedding Trainable readout Target linear functional Target function Feature matrix uX Feature Gram (γID + Φ⊤ Φ)−1 Dual Gram (γIN + ΦΦ⊤ )−1
N ×M N ×1 M ×1 1×D N ×D D×D N ×N
per-coordinate standard deviation j −(α1 +α2 +2)/2 uij ∼ N (0, 1/N ) i.i.d. O(1) O(1) O(1) N/A N/A N/A
Across D i.i.d. samples, let Nj denote the number of times coordinate j activates. Then P(coordinate j is never active in the training set) = P(Nj = 0) = (1 − pj )D ≈ e−Dpj . Thus, whenever Dpj ≫ 1 (equivalently, D j −(α1 +1) ≫ 1), we have P(Nj = 0) ≤ e−Ω(1) ≪ 1, so coordinate j is activated at least once with high probability. When −1 ≤ α1 < 0, this condition holds uniformly for all j up to order D. At the boundary α1 = 0, the coordinate j ∼ D has Dpj = O(1), so the never-observed probability is also O(1); this boundary case only affects constants and does not change the one-exponent scaling. Therefore the Bayes-optimal loss is governed by the unmodeled-variance tail beyond the first min{N, D} coordinates: X X ℓBayes (N, D) ≍ Var(xj ) ≍ j −(α1 +α2 +2) ≍ min{N, D}−(α1 +α2 +1) , j>min{N,D}
j>min{N,D}
and the scaling is symmetric in N and D in this regime, with a single exponent αN = αD = α1 +α2 +1.
C.2
α1 > 0: Two-Exponent Scaling Law
Here we briefly elaborate on supplementary details omitted from the main text regarding the scaling asymptotics of the two-exponent scaling law. 1 In a random set of D data points, only O(D 1+α1 ) features are likely to activate when α1 > 0. Thus, for positive α1 , there are naively three distinct regimes: 1
N < D 1+α1 ,
1
and N > D.
D 1+α1 < N < D,
The first is unambiguously underparameterized, the third is clearly overparameterized, and the intermediate regime exhibits characteristics of both. C.2.1
Underparameterized Regime 1
We first consider the underparameterized regime N ≪ D 1+α1 , where learning is primarily limited by model capacity.
20
Proof of Proposition 3.1 Proof. A resolvent analysis of random-feature regression [Maloney et al., 2022] reveals that the loss is governed by the covariance mass missed by an N -dimensional random feature subspace. In the ridgeless limit, the resolvent acts as a projection onto the latent directions inaccessible to the random features, so the leading-order loss is controlled by the unresolved tail of the covariance spectrum. Thus, at the level of scaling, the random-feature bottleneck resolves the leading N spectral directions and leaves the remaining covariance tail unresolved. In the present regime N ≪ D1/(1+α1 ) , sparsity is masked, so the relevant spectrum is Var(xj ) = j −α1 −α2 −2 . Hence the residual variance scales as Z ∞ X −α1 −α2 −2 j ≍ j −α1 −α2 −2 dj N
j>N
=
1 N −(α1 +α2 +1) . α1 + α 2 + 1
Combining this with the general form of the Bayes-optimal loss (Eq. (3)) yields (4). C.2.2
Overparameterized Regime 1
We now turn to the overparameterized regime N ≫ D 1+α1 . Heavily Overparameterized Regime We first focus on the maximally overparameterized case. We will show that in this regime, the performance of the random feature model is Bayes-optimal among all models, and we compute its generalization error. Optimal Predictor. An upper bound on the performance of any model would be for weff to equal w on indices that appear in the dataset, and be 0 on all other indices. We will show that this is indeed what happens in this regime. If we write out our expression for weff , we have −1 ⊤ ⊤ ⊤ weff = w⊤ X γID + X⊤ u⊤ u X X u u, (16) If we post-multiply by X b for any vector b, we obtain ⊤ weff X b = w⊤ X b.
This follows from the full rank of Q. Thus, in the overparameterized regime, for every possible input in the subspace spanned by the training data, the random feature model incurs zero loss. Coefficients Outside the Training Data Go to 0. While the analysis showing that weff perfectly interpolates the training data applies for any overparameterized model in the limit γ → 0, the argument in this subsection uses the stronger assumption that we are deep in the overparameterized regime. Since uij ∼ N (0, 1/N ), the diagonal entries of u⊤ u are 1 + OP (N −1/2 ), while the off-diagonal entries are OP (N −1/2 ). Thus, for any subspace S with dim(S) ≪ N , the restriction of u⊤ u to S is IS + oP (1). Since the row and column spaces of X are at most D-dimensional, when D ≪ N we can replace the occurrences of u⊤ u in Eq. (16) by I. Therefore, weff simplifies to ⊤ weff = w⊤ X (γID + X⊤ X)−1 X⊤ .
21
In the limit γ → 0, this is precisely the projection of w⊤ onto the space spanned by the data. We can see this using a singular value decomposition. Proof. Let the singular value decomposition (SVD) of the data matrix be X = U ΣV ⊤ , where Σ is the diagonal matrix of singular values σi . Substituting this into the term acting on w⊤ yields −1 X (γID + X⊤ X)−1 X⊤ = U ΣV ⊤ γID + V Σ⊤ ΣV ⊤ V Σ⊤ U ⊤ −1 = U ΣV ⊤ V (γID + Σ⊤ Σ)V ⊤ V Σ⊤ U ⊤ = U Σ(γID + Σ⊤ Σ)−1 Σ⊤ U ⊤ . {z } | Λγ
The inner matrix Λγ is diagonal with entries (Λγ )ii = σi2 /(γ + σi2 ). Taking the limit γ → 0 we find ( 1, if σi ̸= 0, σi2 lim = γ→0 γ + σ 2 0, if σi = 0. i Thus, limγ→0 Λγ acts as the identity on the subspace associated with non-zero singular values. ⊤ Consequently, the full expression becomes Urange Urange , which is the orthogonal projection operator ⊤ is precisely the projection of w⊤ onto the onto the column space of X. Therefore, in this limit, weff space spanned by the training data. Expected Test Loss. In the heavily overparameterized limit, the distribution of w − weff is straightforward. The j-th coefficient is zero if feature j appears at least once in the dataset, and is normally distributed with variance 1 if feature j does not appear. The contribution to the test loss is thus X ℓtest = Egen = Etest [x2j ] Etrain (w − weff )2j . j
Proof of Lemma 3.2 Here we derive the asymptotic expression for the expected number of coordinates that are active at least once across D i.i.d. samples under the sparse activation model. Proof (learnable-coordinate count K(D)). In the sparse activation model, coordinate j is active in a single sample with probability pj = j −α1 −1 and inactive with probability 1 − pj . The probability that coordinate j is never activated in D independent samples is Pnever (j) = (1 − pj )D . For pj ≪ 1, we use (1 − pj )D ≈ e−Dpj as D → ∞; this follows from (1 − Pnever (j) ≈ exp(−Dpj ) = exp −Dj −α1 −1 .
Dpj D −Dpj . Hence D ) →e
Thus the probability that coordinate j is activated at least once is 1 − exp(−Dj −α1 −1 ), and the expected number of such coordinates is K(D) ≈
∞ X
1 − e−Dj
−α1 −1
.
j=1
Passing to the continuum limit, we obtain Z ∞ K(D) ≍
1 − e−Dj
0
22
−α1 −1
dj.
We now perform the change of variables t = Dj
−(α1 +1)
=⇒
j=
D t
1 α1 +1
,
dj = −
α +2 1 1 − 1 D α1 +1 t α1 +1 dt, α1 + 1
which yields K(D) ≍
1 1 D α1 +1 α1 + 1
Z ∞
− α1 +2 1 − e−t t α1 +1 dt.
0
Let β ≡ α11+1 ∈ (0, 1). By integration by parts, Z ∞
(1 − e−t ) t−1−β dt =
0
1 β
Z ∞
e−t t−β dt =
0
Γ(1 − β) . β
Substituting β = α11+1 gives K(D) ≍ Γ 1 −
1 α1 + 1
1
D α1 +1 .
This matches Eq. (5) in Lemma 3.2. Proof of Theorem 3.3 Here we give two derivations of Theorem 3.3: a discrete argument via the cutoff K(D), which provides useful intuition, and a fully rigorous continuous argument. Informal discrete argument via K(D). By Lemma 3.2, when α1 > 0, at most K(D) coordinates of the target vector w can be reliably estimated from D samples. All coordinates with indices j > K(D) remain unobserved and therefore unmodeled. By the unmodeled-variance principle (Eq. (3)), the Bayes-optimal loss is X ℓBayes,D ≍ j −α1 −α2 −2 . j>K(D)
Approximating the sum by an integral, X
j −α1 −α2 −2 ≍
Z ∞
j −α1 −α2 −2 dj
K(D)
j>K(D)
=
1 K(D)−(α1 +α2 +1) . α1 + α2 + 1
Substituting the scaling K(D) ≍ D1/(α1 +1) from Eq. (5) yields ℓBayes,D ≍ D−(α1 +α2 +1)/(α1 +1) , which is Eq. (6). Rigorous proof (continuous argument). An equivalent expression for the Bayes-optimal loss can be written as a continuous expectation over coordinates, weighted by the probability that each coordinate is never activated. This leads to the integral representation Z ∞ ℓBayes,D = exp −Dj −α1 −1 j −(α1 +α2 +2) dj, (17) 0
23
We again apply the change of variables 1 D α1 +1 −(α1 +1) , t = Dj =⇒ j = t
dj = −
α +2 1 1 − 1 D α1 +1 t α1 +1 dt. α1 + 1
Substituting into (17) gives α +α +1 1 − 1 2 ℓBayes,D = D α1 +1 α1 + 1
Z ∞
e−t t
α1 +α2 +1 −1 α1 +1
dt.
0
The remaining integral is exactly the Gamma function, Z ∞ e−t tβ−1 dt = Γ(β), 0
with
α1 + α 2 + 1 . α1 + 1 α +α +1 α1 + α 2 + 1 1 − 1 2 Γ D α1 +1 , ℓBayes,D = α1 + 1 α1 + 1 β=
Hence
which matches the scaling law in Eq. (6) in the main text. When α1 > 0, sparsity reduces the effective number of observed coordinates to scale as D1/(α1 +1) , thereby becoming the dominant data bottleneck. This concludes the proof. Subtleties Concerning Theorem 3.3: Co-activation and Identifiability. A subtle complication arises when two rare coordinates co-activate, i.e., are both nonzero in the same training example. In such cases, it becomes impossible to disentangle their individual contributions to the label, creating ambiguity in estimating the corresponding weights. The coordinates that activate exactly once 1 have indices on the order of j ∼ D α1 +1 . Consequently, the total number of such coordinates also 1 scales as D α1 +1 . A nonzero number of pairwise co-activations is therefore likely whenever this quantity exceeds D1/2 (by the birthday paradox), which occurs when α1 < 1. However, the number of co-activating coordinates itself scales only as 2 1 1−α D α1 +1 o n 1 1 1 2 1 −1 α +1 α +1 α +1 α +1 α +1 min , D 1 = min D 1 , D 1 . , D 1 = min D 1 D −α1 −α2 −2
Each such coordinate contributes variance on the order of D α1 +1 , so the total variance contributed by co-activating coordinates scales as −2α −α −1 −α1 −α2 −1 1 2 α +1 α +1 1 1 min D , D , which vanishes as D → ∞. Thus, while co-activation may create isolated ambiguity, it does not affect the asymptotic loss scaling in Theorem 3.3.
D
Compute
This section provides an overview of the compute model we adopt for studying the compute-optimal frontier under training to completion: C = N D · min{N, D}. 24
(18)
D.1
Computational Efficiency
Here, “compute” refers to the leading-order training cost, measured in FLOPs up to hardwaredependent constants. To approach the population-loss scaling studied in the main text, we focus on the setting in which the readout is trained to convergence. This can be done either with gradientbased methods, which are closer to practical training, or by directly computing the min-norm solution in Eq. (15). D.1.1
Gradient-Based Methods
A single gradient step on the loss (14) costs O(N D). The iteration complexity is controlled by the condition number of the feature Gram matrix, defined on its nonzero spectrum, κ(ΦΦ⊤ ) = + ⊤ λmax (ΦΦ⊤ )/λ+ min (ΦΦ ), where λmin denotes the smallest nonzero eigenvalue. For plain GD, the number of steps scales linearly in κ, up to logarithmic accuracy factors. Accelerated first-order √ methods improve this dependence to κ. We therefore write the leading gradient-based training cost as Cgrad ∼ N D κp , (19) with p = 1 for GDand p = 1/2 for accelerated methods. Condition number, κ. We estimate the condition number from the effective variance profile of the active features. The number of active coordinates in a dataset of size D is K(D), so the readout effectively uses roughly min{N, K(D)} features. The variance explained by the smallest feature scales as min{N, K(D)}−(α1 +α2 +2) , while the largest is O(1). Thus the effective condition number scales as κ ∼ min{N, K(D)}α1 +α2 +2 . (20) Consequently, if the iteration complexity scales as κp , the number of optimization steps scales as min{N, K(D)}(α1 +α2 +2)p . D.1.2
Direct Least-Squares Solve
Eq. (15) gives two equivalent expressions for θ ⋆ , related by a standard matrix identity. Using the first expression, y Φ⊤ q, requires computing ΦΦ⊤ (cost N 2 D), inverting it (cost N 3 ) to form q, computing yΦ⊤ (cost N D), and then multiplying the resulting row vector by q (cost N 2 ). The total cost is therefore O(N 2 D + N 3 ). By the same logic, using the second expression, y Q Φ⊤ , costs O(D2 N + D3 ). The cheaper of these two options is Cdirect ∼ N D · min{N, D}. We summarize these estimates in Table 3, which compares the costs of gradient descent (GD), accelerated first-order methods, and direct solvers across the three regimes.
25
Table 3: Compute cost, normalized by N D, across the three parameter regimes. The GD and accelerated columns assume iteration complexity scaling as κp , with p = 1 for GD and p = 1/2 for accelerated first-order methods such as Nesterov acceleration. The direct-solve column corresponds to computing the min-norm least-squares solution in closed form. Regime
GD
Heavily overparameterized, N > D
D
D
α1 +α2 +2 2(α1 +1) α1 +α2 +2 2(α1 +1)
α1 +α2 +2 α1 +1
Intermediate, D > N > K(D) Underparameterized, K(D) > N
D.2
Accelerated
α1 +α2 +2 α1 +1
D N α1 +α2 +2
D α1 +α2 +2 2 N
Direct solve D N N
Proxy Compute Model
Table 3 implies different computational scalings depending on both the solver and the regime. Let a = α1 + α2 + 2 and b = α1 + 1, with b > 0. The comparison between direct solvers and first-order methods is solver-dependent. Since accelerated methods reduce the condition-number dependence from κ to κ1/2 , they are asymptotically cheaper than GD whenever a > 0, up to constants. The direct solve is cheaper than GD in the heavily overparameterized and intermediate regimes when a/b > 1, equivalently α2 + 1 > 0, while in the underparameterized regime the finite-variance condition a > 1 already guarantees that the direct solve is cheaper than GD. The comparison with accelerated methods is stricter. In the heavily overparameterized regime, the direct solve beats acceleration only when a/(2b) > 1, equivalently α2 > α1 . In the underparameterized regime, it beats acceleration only when a/2 > 1, equivalently α1 + α2 > 0. Thus direct solvers are uniformly cheapest under the stronger condition α2 > α1
and
α1 + α2 > 0.
(21)
There is also a mixed regime in which GD is the most expensive method, but the direct solve beats accelerated methods only in part of the phase diagram. For example, when α1 > 0 and max{−1, −α1 } < α2 < α1 ,
(22)
acceleration is cheaper than the direct solve in the heavily overparameterized regime, while the direct solve is cheaper in the underparameterized regime. In the intermediate regime, the crossover occurs at Nc (D) ∼ D(α1 +α2 +2)/(2(α1 +1)) : direct solution is cheaper for N < Nc (D), while acceleration is cheaper for N > Nc (D). For example, the representative choice (α1 , α2 ) = (1, 0.3) used in our experiments lies in this mixed regime: GD is asymptotically the most expensive method, while the cheapest method switches from acceleration in the heavily overparameterized regime to direct solution in the underparameterized regime, with a crossover between the two in the intermediate regime. This solver dependence motivates using a simple unified compute proxy, C = N D · min{N, D},
(23)
rather than tying the main scaling analysis to a particular optimization algorithm. The qualitative conclusion we emphasize below is that sparsity shifts the compute-optimal allocation toward data, and this conclusion does not depend on these solver-dependent distinctions within reasonable compute models.
26
D.2.1
Robustness to Alternative Compute Models
The solver-dependent comparisons above change the mapping from (N, D) to compute, but not the basic allocation principle. The loss has the form ℓ(N, D) ∼ N −αN + D−αD ,
(24)
so the compute-optimal allocation balances the two terms: N −αN ∼ D−αD
D ∼ N αN /αD = N α1 +1 .
=⇒
(25)
Thus sparsity fixes the relative allocation between data and model size independently of the particular solver. Different compute models only change how this balanced allocation scales with the total budget C. For example, suppose that in the underparameterized branch the compute model takes the more general form C ∼ N r D. (26) Combining this with D ∼ N α1 +1 gives N ∗ (C) ∼ C 1/(r+α1 +1) , and
D∗ (C) ∼ C (α1 +1)/(r+α1 +1) ,
ℓ∗ (C) ∼ C −αN /(r+α1 +1) .
(27) (28)
The direct-solve proxy used in Proposition 5.1 corresponds to r = 2, recovering Eq. (10). An accelerated first-order method in the same branch would instead have r = 1+ α1 +α2 2 +2 , which changes the numerical compute exponent but leaves the balanced allocation D∗ ∼ (N ∗ )α1 +1 unchanged. This distinction is especially useful in the mixed regime discussed above, where the cheapest solver can switch across the (N, D) plane. In that regime, acceleration is cheaper than the direct solve only in sufficiently model-heavy regions. The crossover in the intermediate regime occurs at α1 +α2 +2
Nc (D) ∼ D 2(α1 +1) . By contrast, the loss-balanced allocation satisfies Nbal (D) ∼ D1/(α1 +1) . In the 1 +α2 +2 mixed regime with α1 + α2 > 0, we have α11+1 < α2(α , so Nbal (D) ≪ Nc (D) asymptotically. 1 +1) Therefore the compute-optimal allocation lies on the direct-solve side of the crossover, even though acceleration is cheaper in more heavily overparameterized regions. For the representative experimental choice (α1 , α2 ) = (1, 0.3), this comparison gives Nbal (D) ∼ D1/2 and Nc (D) ∼ D0.825 . Thus the balanced allocation is well below the accelerated/direct crossover. Consequently, even under a piecewise cheapest-solver compute model, the asymptotic compute-optimal allocation is governed by the same direct-solve branch used in Proposition 5.1. More generally, alternative solver models can modify the compute exponent αC , but they do not alter the qualitative conclusion that sparsity shifts the compute-optimal allocation toward data.
D.3
Compute-Optimal Allocation
Proposition 5.1 adopts the proxy compute budget in Eq. (23), C = N D · min{N, D}, and shows that the unique compute-optimal allocation lies in the underparameterized regime, N < D. The optimal allocation and resulting loss scale as ℓ∗ (C) ∼ C −αC , N ∗ (C) ∼ C αD /(αN +2αD ) = C 1/(α1 +3) , D∗ (C) ∼ C αN /(αN +2αD ) = C (α1 +1)/(α1 +3) , αN αD α1 + α2 + 1 αC := = > 0. αN + 2αD α1 + 3 27
(29)
D.3.1
Sketch of Proposition 5.1
Sketch (compute-optimal scaling). We optimize the loss under the compute model C = N D · min{N, D}, Underparameterized Side (N < D). Here the compute constraint becomes C = N D · N = N 2 D. Assuming the population loss when α1 > 0 takes the form ℓ(N, D) ∼ N −αN + D−αD , we eliminate D using the compute constraint, D = NC2 : ℓ(C, N ) ∼ N
−αN
+
C N2
−αD
= N −αN + C −αD N 2αD .
For fixed C, the first term decreases with N while the second increases with N (since larger N forces smaller D = C/N 2 ). Thus the optimum balances the two contributions. Equating their scalings, N −αN ≍ C −αD N 2αD yields
αD
N ∗ (C) ≍ C αN +2αD ,
N αN +2αD ≍ C αD ,
⇐⇒
D∗ (C) =
C (N ∗ (C))2
≍ C
1− α
2αD N +2αD
.
At this optimum, both loss terms scale equally, so ℓ∗ (C) ∼ (N ∗ (C))−αN ≍ C −αC ,
αC =
αN αD , αN + 2αD
which gives the exponents reported in Eq. (29). Overparameterized Side (D < N ). In this case min{N, D} = D and the compute constraint becomes C = N D2 . Optimizing under this branch produces a candidate scaling that is not selfconsistent with D < N , and is therefore ruled out by Remark 5.2, so the unique compute-optimal allocation lies in the underparameterized regime.
E
Gradient Descent
Gradient descent (GD) provides a natural procedure for minimizing the training loss. In favorable settings, it converges to the trained-to-completion min-norm solution, though rare failures can occur due to sparsity as shown in the main text. Below, we provide a proof of Proposition 6.1 and a proof sketch for Theorem 6.2, which together summarize our main findings on GD, as reported in the main text.
28
E.1
Proof of Proposition 6.1
Proof (spectral bound for stability). Under the near-isometry approximation on span(X), we may 1 ⊤ ⊤ ⊤ ⊤ replace u u by IM inside u u XX , so it suffices to control λmax D XX . 1 If empirical second moments concentrate, then D XX⊤ ≈ Σx with Σx = E[xx⊤ ] diagonal and 1 Var(xj ) = j −(α1 +α2 +2) , whose top eigenvalue equals Var(x1 ) = 1. Although Σx is diagonal, D XX⊤ need not be. The off-diagonal entries have mean zero and concentrate entrywise; we assume that their aggregate contribution is subleading, so they do not change the leading O(1) scale of the top eigenvalue governed by the diagonal part. The only potential obstruction is an anomalously large diagonal entry caused by a single rare activation at a very large index when α2 < −1 (so activated amplitudes grow with j). Let jmax be the largest index that activates at least once in the dataset. For a given threshold J, the expected number of activations at indices larger than J is D
X
P(xj ̸= 0) = D
j>J
X
j −(α1 +1) ≍
j>J
1 D J −α1 . α1
Equivalently, writing N>J :=
D X X
(d)
1{xj
̸= 0},
d=1 j>J
a Poisson approximation yields D −α1 , P(N>J = 0) ≈ exp − J α1
D −α1 P(jmax ≤ J) = P(N>J = 0) ≈ exp − J , α1
and hence P(jmax = J) ≈ D J
−(α1 +1)
D −α1 exp − J . α1
In particular, jmax ≍ D1/α1 with high probability. 1 Conditional on an activation at jmax , the corresponding diagonal contribution to D XX⊤ is of order 1 2 1 −(α2 +1) xjmax = jmax ∼ D−(α1 +α2 +1)/α1 . D D But, given that α2 < −1, and since α1 + α2 + 1 > 0, we have D−(α1 +α2 +1)/α1 < 1, so such rare spikes cannot dominate (inparticular j = 1). Therefore the O(1) contribution from the leading coordinates 1 1 ⊤ ⊤ λmax D XX = 1 + oP (1), and the same holds for λmax D u u XX⊤ under near-isometry. The conclusion of Proposition 6.1 then follows from the usual spectral stability condition on the step size.
E.2
Sketch of Theorem 6.2
Sketch (rare-activation spike). We focus on the event that a single “late” coordinate activation produces a spiked top eigenvalue. Write the data matrix as X = [ Xreg xspec ], where Xreg ∈ RM ×(D−1) contains the regular samples and xspec ∈ RM is a special sample containing a single unusually large activation. Then ⊤ XX⊤ = Xreg X⊤ reg + xspec xspec .
29
In a spiked covariance scenario, the top eigenvector aligns with xspec and the corresponding eigenvalue is bounded below by the Rayleigh quotient λmax (XX⊤ ) ≥
⊤ ⊤ x⊤ x⊤ spec XX xspec spec Xreg Xreg xspec = + ∥xspec ∥22 . ∥xspec ∥22 ∥xspec ∥22
(30)
Under the finite-variance condition α1 + α2 + 1 > 0, the first term in (30) is o(D) at the rare-spike scale, so the condition for GD instability, λmax XX⊤ > 2D η , is to leading order in D equivalent to requiring 2D ∥xspec ∥22 ≳ . η In our model, a single-coordinate activation at index j has magnitude ∥xspec ∥22 ≈ x2j = j −(α2 +1) . Hence (E.2) requires 1 2D −α2 −1 ∗ j ≳ j := (α2 < −1). (31) η Finally, let N>j ∗ denote the (random) number of activations with index larger than j ∗ in the dataset. By definition of j ∗ , the rare-spike divergence event is the existence of at least one activation with index > j ∗ , so Prare (diverge) ≍ P(N>j ∗ ≥ 1). The expected number of activations at indices larger than j ∗ is E[N>j ∗ ] = D
X
P(xj ̸= 0) = D
j>j ∗
X
j −(α1 +1) ≍
j>j ∗
1 D (j ∗ )−α1 . α1
Equivalently (as in Subsection E.1), a Poisson approximation gives P(N>j ∗ ≥ 1) = 1 − exp(−E[N>j ∗ ]) ≍ E[N>j ∗ ] when E[N>j ∗ ] is small, hence P(N>j ∗ ≥ 1) ≍ α11 D (j ∗ )−α1 . Substituting (31) gives 1 Prare (diverge) ≍ D α1
2D η
−
α1 −α2 −1
α +α +1
≍D
1 2 − −α −1 2
,
which is Eq. (13) up to η-dependent constants. E.2.1
Empirical Considerations for Theorem 6.2
We discuss the practical difficulty of empirically validating Theorem 6.2’s rare-spike prediction Prare (diverge) ≍ D−ν , and report our finite-D measurements, which are consistent with the theorem’s qualitative form. A clean validation would estimate P λmax (XX⊤ /D) > ϵ for fixed threshold ϵ = 2/η across a range of D values spanning the asymptotic regime, then fit the slope on log-log axes and compare to the predicted exponent ν = (α1 + α2 + 1)/(−α2 − 1). However, this represents a statistical challenge: the theorem characterizes a tail probability, and at any given D, accurate estimation of Prare (diverge) = p requires Θ(1/p) independent samples to observe events at the relevant rate, with sample size scaling further with ν for accurate exponent estimation. For the theorem’s asymptotic regime to dominate observation, D must be large enough that finite-D corrections to the spectral concentration are negligible compared to the leading D−ν behavior; in our experiments, this asymptotic crossover occurs beyond D ≈ 104 . Combining these requirements, a quantitative validation at a predicted exponent ν = 9 (for chosen parameters α1 = 2, α2 = −1.2) would require roughly 109 to 1012 independent dataset samples per D value—many orders of magnitude beyond computational tractability with eigenvalue computation on D × D Gram matrices. 30
We perform experiments in which we sweep D ∈ {500, 1000, 2000, 5000, 10000} at fixed threshold ϵ = 1.0005 and feature dimension M = 5000, sampling between 500 and 4000 independent dataset replicates per D. We chose α1 = 2, α2 = −1.2 (giving ν = 9) as a representative point well inside the regime α2 < −1 and α1 + α2 + 1 > 0 where Theorem 6.2 applies (smaller values of ν are more easily detected in finite samples but tend to give less clean separation between tail-event probabilities and bulk concentration). The threshold ϵ was chosen empirically so that the smallest-D measurement gives a moderate (non-saturated) failure rate. Our results are shown in Table 4. Table 4: Empirical divergence probability as a function of dataset size D, with α1 = 2, α2 = −1.2, M = 5000, threshold ϵ = 1.0005. D
failures / seeds
P(diverge)
500 1000 2000 5000 10000
185/500 132/1000 40/2000 1/3000 0/4000
0.370 0.132 0.020 0.00033 < 0.00025
The decay spans more than three orders of magnitude in P(diverge) over an order of magnitude in D. Computing local exponents from successive D pairs gives νlocal ≈ 1.5, 2.7, 4.5, increasing monotonically—consistent with finite-D corrections that vanish as the asymptotic regime is approached. Two observations matter for the comparison to theory. First, the empirical P(diverge) exhibits a clean, monotone power-law decay with D—consistent in functional form and direction with Theorem 6.2. Second, the empirical exponent at the largest measurable D (νlocal ≈ 4.5) is approaching, but has not reached, the asymptotic prediction ν = 9. Both observations are consistent with the theorem: the asymptotic exponent is a statement about D → ∞, and the local empirical exponent’s monotonic increase suggests we are inside the finite-D regime where the asymptotic rate has not yet emerged. Distinguishing the asymptotic exponent ν = 9 from any other large-ν value (e.g., ν = 6 or ν = 12) requires sample sizes and D ranges substantially beyond what is computationally accessible with this experimental design. We conclude that our measurements are quantitatively non-contradictory with the theorem and qualitatively confirm its prediction of a sharp power-law decay in the divergence probability with dataset size.
F
Experiments
F.1
Min-norm least-squares solutions via gradient-based optimization
Our scaling-law theory predicts the test loss of the min-norm random-feature readout as N and D vary. We use two independent methods to recover this solution: 1 1. Closed-form pseudoinverse. For the linear regression problem minθ 2D ∥Φ⊤ θ − y⊤ ∥22 on (linear ⋆ or ReLU) random features, the min-norm least-squares solution is θ = (Φ⊤ )+ y⊤ (equivalently θ ⋆⊤ = yΦ+ ), computed via SVD-based pseudoinverse with relative tolerance 10−10 . 2. Nesterov-accelerated GD with adaptive restart [Nesterov, 1983, O’Donoghue and Candès, 2015]. The readout θ is initialized at zero and trained on the empirical mean-squared error with Nesterov momentum and gradient-based adaptive restart. From zero initialization, the converged solution is provably the same min-norm least-squares solution as the closed-form one.
31
Sparse Figures 2 and 3 are computed using the closed-form pseudoinverse, while Figure 4 10 2 is computed using the Nesterov-accelerated 10 2 solver. 3 10 The closed-form and Nesterov approaches 10 3 10 4 give the same qualitative two-regime struc10 5 10 4 ture, supporting the interpretation that the 101 101 102 102 103 N D observed exponent asymmetry is primarily Dense a property of the regression problem rather 1 10 than of a particular solver. In the linear 10 2 10 2 case (Figure 5), the joint exponent across the 10 3 10 3 sparse N -sweep, dense N -sweep, and dense 10 4 10 4 D-sweep recovers α ≈ 2.0, in close agree10 5 10 5 ment with the linear-theory prediction αN = 10 6 10 6 α1 + α2 + 1 = 2.30, while the sparse D-sweep 101 101 102 103 102 103 N D gives αD ≈ 1.11, matching the theoretical prediction αD = (α1 + α2 + 1)/(α1 + 1) = 1.15. Figure 5: Two-exponent scaling validates linear The asymmetry ratio αN /αD ≈ 1.80 is close theory. Test loss scaling under the linear feature to the predicted 2.30/1.15 = 2.00, with a residmap ϕ(x) = ux, computed via the closed-form minual gap plausibly attributable to finite-size cornorm least-squares solution; (α1 , α2 ) = (1.0, 0.3), 20 rections and possible crossover effects near the seeds. Top: sparse; bottom: dense. Left: N -sweep at double-descent peak. Under ReLU (Figures 6 D = 50,000; right: D-sweep at N = 16,000. A single and 4) both methods give the same qualiexponent α ≈ 2.0 (dashed lines, jointly fitted with tative structure: a joint exponent describing separate intercepts) describes the sparse N -, dense three of the four sweeps and a shallower sparse N -, and dense D-sweeps, in close agreement with the D-exponent breaking the dense-baseline symlinear-theory prediction αN = α1 + α2 + 1 = 2.30. metry, with the magnitudes αN /αD ≈ 1.13 The sparse D-sweep requires its own shallower fit, (closed-form) and ≈ 1.25 (Nesterov). These αD ≈ 1.11, matching the predicted sparse exponent finite-size corrections are roughly consistent αD = (α1 + α2 + 1)/(α1 + 1) = 1.15 and breaking across the linear and ReLU cases, indicating the dense-baseline symmetry αN = αD . that nonlinearity reduces the magnitude of the asymmetry but does not alter its qualitative structure.
Setup. All experiments use the data-generating model defined in Section 2.1, with sparse parameters (α1 , α2 ) = (1.0, 0.3) and ambient dimension M = 10,000. The teacher is y(x) = w⊤ x with wj ∼ N (0, 1) i.i.d. The first-layer weights u ∈ RN ×M are drawn i.i.d. from N (0, 1/N ) and held frozen throughout training. The dense baseline uses Gaussian inputs with per-coordinate variance j −(α+1) , α = α1 + α2 + 1, matching the underparameterized exponent of the sparse model so that any asymmetry in scaling is attributable purely to sparsity. Test loss is evaluated on ntest fresh samples drawn from the same distribution. For each (N, D) configuration we run multiple seeds (varying both data sampling and the random first layer) and report mean test loss with ±1 standard error. Power-law exponents are fitted on log-log axes by least squares, restricted to the asymptotic tail (last 4–6 values of N or D) to mitigate finite-size deviations. As a check, we also compare full-batch GD with minibatch SGD at small learning rate; the resulting trajectories and scaling exponents agree closely, suggesting that the sparse scaling laws reported here are not driven primarily by minibatch gradient noise, but by the underlying optimization/statistical structure.
32
Experiment 1: Linear random features. The forward map is ŷ = θ⊤ (ux). The min-norm readout is solved in float64 via torch.linalg.pinv with rtol = 10−10 to handle rank-deficient feature matrices robustly. We sweep N ∈ {4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256} at fixed D = 50,000, and D ∈ {10, 14, 20, 28, 40, 56, 80, 112, 160, 224, 320, 448, 640, 896, 1280} at fixed N = 16,000, with 20 seeds per configuration and ntest = 50,000. Fits use the last 6 points of each sweep.
Sparse 10 1
10 1 10 2
10 2
10 3 10 3
10 4 101
N
102
101
103
102
D
103
Dense 10 1 10 1
10 2
10 2 10 3 10 3
Experiment 2: ReLU random features. 10 4 10 4 ⊤ The forward map is ŷ = θ σ(ux) with σ 101 101 102 103 102 103 N D the ReLU activation. Two solver variants are used: Figure 6: Asymmetry persists under nonlinearity (closed-form). Test loss scaling under Closed-form (pinv). Identical to Ex- the ReLU feature map ϕ(x) = σ(ux), computed periment 1, with ReLU applied to the features via the closed-form min-norm least-squares solution; before the pseudoinverse. Same sweep grids, (α1 , α2 ) = (1.0, 0.3), 20 seeds. Top: sparse; bot20 seeds, and ntest = 50,000. tom: dense. Left: N -sweep at D = 50,000; right: D-sweep at N = 16,000. A single exponent α ≈ 1.25 Iterative (Nesterov). The ReLU fea- (dashed lines, jointly fitted with separate intercepts) tures Φ = σ(uX) ∈ RN ×D are precomputed describes the sparse N -, dense N -, and dense Donce in float64 on GPU. The step size is set sweeps. The sparse D-sweep requires its own shalto η = 1/λmax (ΦΦ⊤ /D), where λmax is esti- lower fit, αD ≈ 1.10, breaking the dense-baseline mated by 50 steps of power iteration. Nes- symmetry αN = αD predicted by linear theory; the terov updates with O’Donoghue–Candès adap- smaller magnitude of the asymmetry compared to tive restart [O’Donoghue and Candès, 2015] the linear case (Figure 5) reflects spectral smoothing are applied until the relative gradient norm under nonlinear feature maps. satisfies ∥∇∥ < 10−9 ∥∇0 ∥, with a maximum of 5 × 105 iterations. We sweep N ∈ {4, 8, 16, 32, 64, 128, 256} at fixed D = 50,000, and D ∈ {10, 20, 40, 80, 160, 320, 640, 1280} at fixed N = 8000, with 5 seeds per configuration and ntest = 20,000. Fits use the last 4 points of each sweep.
Hardware and reproducibility. All experiments were run on a single NVIDIA T4 or A100 GPU (Google Colab). Total wall-clock time was approximately 4 hours for the closed-form sweeps (Experiments 1 and 2 closed-form) and 30 minutes for the Nesterov sweeps. Float64 precision is used throughout the optimization and pseudoinverse computations; data sampling and feature evaluation use float32 with TF32 matrix multiplications enabled. The first layer u is regenerated per seed.
33