Paper deep dive
Learning with Synthetic Data via SGD in High-Dimensional Linear Regression
Jichu li, Difan Zou
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 85%
Last extracted: 9/10/2026, 2:57:01 AM
Summary
This paper analyzes the impact of synthetic data on the generalization of one-pass Stochastic Gradient Descent (SGD) in high-dimensional linear regression. It compares two training protocols: mixed training (randomly mixing real and synthetic data) and two-stage training (synthetic data first, then real data). The authors prove that mixed training induces 'strong model collapse,' characterized by a non-vanishing excess risk floor due to fluctuation and drift errors caused by source mismatch. Conversely, two-stage training avoids this collapse, allowing the model to achieve vanishing excess risk. The study derives finite-sample risk bounds and scaling laws, showing that larger models may amplify degradation in mixed training, while high-quality synthetic pretraining can reduce bias in two-stage training. An exact condition is provided for when two-stage training strictly outperforms real-only training.
Entities (10)
Relation Signals (7)
Two-Stage Training → avoids → Model Collapse
confidence 95% · two-stage training avoids the floor by using synthetic data only in the first stage
Mixed Training → induces → Model Collapse
confidence 95% · mixed training induces strong model collapse, while two-stage training avoids the floor
Source Mismatch → causes → Fluctuation Error
confidence 90% · source-mismatch effects, namely fluctuation and persistent drift under mixing
Source Mismatch → causes → Drift Error
confidence 90% · separating standard bias and variance from source-mismatch effects, namely fluctuation and persistent drift
Two-Stage Training → outperforms → Real-Only Training
confidence 90% · two-stage training to strictly outperform real-only training under the same real-data budget
Larger Models → amplify → Synthetic-Induced Degradation
confidence 88% · larger models may amplify synthetic-induced degradation under mixing
High-Quality Synthetic Pretraining → reduces → Bias
confidence 85% · high-quality synthetic pretraining may reduce bias in two-stage training
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Synthetic data has become a promising way to scale model training beyond limited human-generated data but it may also induce strong model collapse (Dohmatob et al., 2024), where any fixed fraction of synthetic data prevents model performance from improving under data scaling, leaving a non-vanishing excess risk floor. In this paper, we study how synthetic data affects the generalization of one-pass SGD in high-dimensional linear regression with model shift. We establish finite-sample risk bounds for mixed and two-stage training, separating standard bias and variance from source-mismatch effects, namely fluctuation and persistent drift under mixing and filtered initialization bias under two-stage. These bounds reveal a sharp contrast: mixed training induces strong model collapse, while two-stage training avoids the floor by using synthetic data only in the first stage, showing that collapse is not inevitable under a simple data curriculum. Under a random sketch model, we further obtain scaling laws for both protocols, with tight results for mixed training in the optimization-saturated regime. These laws show that larger models may amplify synthetic-induced degradation under mixing, and quantify how high-quality synthetic pretraining may reduce bias in two-stage training. Finally, we establish an exact finite-sample necessary-and-sufficient condition for two-stage training to strictly outperform real-only training under the same real-data budget and identical real-stage updates. Overall, our results highlight that synthetic data is neither inherently harmful nor beneficial; its effect depends critically on both its quality and the training protocol used to incorporate it.
Tags
Links
- Source: https://arxiv.org/abs/2609.09572v1
- Canonical: https://arxiv.org/abs/2609.09572v1
Trouble viewing inline? Open PDF directly →
Full Text
278,211 characters extracted from source content.
Expand or collapse full text
Learning with Synthetic Data via SGD in High-Dimensional Linear Regression Jichu Li Affiliation: Renmin University of China Email: lijichu52@gmail.com Difan Zou Affiliation: The University of Hong Kong Email: dzou@hku.hk Abstract Synthetic data has become a promising way to scale model training beyond limited human-generated data but it may also induce strong model collapse (Dohmatob et al., 2024b), where any fixed fraction of synthetic data prevents model performance from improving under data scaling, leaving a non-vanishing excess risk floor. In this paper, we study how synthetic data affects the generalization of one-pass SGD in high-dimensional linear regression with model shift. We establish finite-sample risk bounds for mixed and two-stage training, separating standard bias and variance from source-mismatch effects, namely fluctuation and persistent drift under mixing and filtered initialization bias under two-stage. These bounds reveal a sharp contrast: mixed training induces strong model collapse, while two-stage training avoids the floor by using synthetic data only in the first stage, showing that collapse is not inevitable under a simple data curriculum. Under a random sketch model, we further obtain scaling laws for both protocols, with tight results for mixed training in the optimization-saturated regime. These laws show that larger models may amplify synthetic-induced degradation under mixing, and quantify how high-quality synthetic pretraining may reduce bias in two-stage training. Finally, we establish an exact finite-sample necessary-and-sufficient condition for two-stage training to strictly outperform real-only training under the same real-data budget and identical real-stage updates. Overall, our results highlight that synthetic data is neither inherently harmful nor beneficial; its effect depends critically on both its quality and the training protocol used to incorporate it. 1 Introduction Large language models (LLMs) have demonstrated remarkable capabilities, primarily driven by large-scale, high-quality pre-training data (Radford et al., 2019; Brown et al., 2020). Scaling laws suggest that increasing both model size and data volume consistently leads to improved performance (Kaplan et al., 2020; Hoffmann et al., 2022; Muennighoff et al., 2023). However, high-quality pre-training data is limited, and its supply is rapidly diminishing, posing a significant challenge for continued model scaling (Muennighoff et al., 2023; Villalobos et al., 2024). Consequently, synthetic data, often generated by existing models, has emerged as a promising supplement to human-generated data for training, particularly in scenarios where labeled data is scarce (Muennighoff et al., 2023; Chen et al., 2024; Long et al., 2024; Havrilla et al., 2024). The training paradigm thus evolves to operate on a mixture of real data from 1P_1 and synthetic data from 2P_2, while our ultimate interest remains performance on the real data distribution 1P_1 alone. This raises a natural question: Does synthetic data truly benefit model training, and what limitations or adverse effects might arise? On the one hand, a large body of work has demonstrated the utility of synthetic data. By amplifying real datasets, synthetic samples effectively increase the total training volume, enabling models to learn richer patterns and improve generalization (Gunasekar et al., 2023; Li et al., 2023b; Liu et al., 2024; Maini et al., 2024; Long et al., 2024; Abdin et al., 2024; Qin et al., 2025). On the other hand, Shumailov et al. (2023) showed that training with synthetic data can induce model collapse, a phenomenon where the iterative use of generated data for training leads to a critical degradation in model performance. This phenomenon arises as the model gradually overfits to patterns found in synthetic data, which may not fully represent the richness or variability of real-world data, ultimately leading to a deterioration of model quality over time, as extensively validated by prior empirical works (Hataya et al., 2023; Bohacek and Farid, 2023; Alemohammad et al., 2023; Guo et al., 2024; Briesch et al., 2024; Shumailov et al., 2024). Understanding training with synthetic data has become central to large-scale model development. Theoretical studies of learning with synthetic data typically formulate the problem as a regression task under distribution mismatch between real and synthetic data. Some of them focus on the phenomenon of model collapse, showing how training with synthetic data can degrade model performance (Dohmatob et al., 2024c; Dohmatob et al., 2024a; Seddik et al., 2024; Jain et al., 2024) while other studies investigate strategies to mitigate model collapse or to effectively leverage synthetic data to improve learning (Seddik et al., 2024; Firdoussi et al., 2024; Barzilai and Shamir, 2025; Garg et al., 2025; Ferbach et al., 2024). Of particular relevance to our work, Dohmatob et al. (2024b) introduce the concept of strong model collapse: even a small proportion of synthetic data prevents model performance from improving under data scaling, leaving a non-vanishing excess risk floor, which serves as a key motivation for our analysis. However, most of these works focus on static estimators and asymptotic regimes, without capturing the dynamics of practical training procedures based on iterative optimization algorithms (e.g., SGD), nor providing a detailed finite-sample characterization. In this paper, we study the effect of synthetic data on model generalization through the lens of one-pass SGD in high-dimensional linear regression. We consider a two-source setting where real and synthetic data share the same covariate distribution but differ in their labeling functions, capturing common synthetic-supervision pipelines where real inputs are labeled by a biased teacher, and analyze SGD under both mixed and two-stage training protocols. We establish finite-sample upper and lower bounds showing that mixed training leads to strong model collapse under SGD, while a simple two-stage protocol avoids this issue and achieves vanishing excess risk. Further, we derive scaling laws under a random sketch model, characterizing how model size may affect synthetic-induced degradation and how synthetic data may benefit learning through bias reduction. Finally We give an exact finite-sample condition for two-stage training to strictly outperform real-only training with the same real-data budget and real-stage updates. Our contributions can be summarized as follows: • In Section 4.1, we derive finite-sample upper and lower bounds for last-iterate one-pass SGD under mixed training. The upper bound decomposes the synthetic-data-induced degradation into fluctuation and drift, characterizing how source mismatch affects optimization dynamics. Together, the upper and lower bounds establish a non-vanishing excess risk floor under data scaling, which is exactly strong model collapse. Moreover, we show that the same behavior persists for iterate-averaged SGD, indicating that strong model collapse is fundamentally induced by sample mixing in the training protocol, rather than an artifact of a specific optimization scheme. • In Section 4.2, we show that strong model collapse is not inevitable. We analyze a two-stage training protocol, where the model is first trained on synthetic data and then on real data, and show that this simple data curriculum avoids the non-vanishing excess risk floor. This highlights that the ordering and scheduling of data can be as important as the data itself. • In Section 4.3, we further establish scaling-law bounds for the two training protocols under a random sketch model, with tight results for mixed training in the optimization-saturated regime, showing that larger models may amplify synthetic-data-induced degradation in mixed training, and providing an explicit characterization of how high-quality synthetic data may reduce bias in two-stage training, thereby benefiting learning. • In Section 4.4, we establish an exact finite-sample condition under which two-stage training strictly outperforms real-data-only training with a fixed real-data budget and identical real-stage training; we also use a one-dimensional toy model to illustrate the roles of synthetic sample size, source mismatch, and label noise. Notations. For two positive-valued functions f(x)f(x) and g(x)g(x), we write f(x)≲g(x)f(x) g(x) (and f(x)=(g(x))f(x)=O(g(x))) or f(x)≳g(x)f(x) g(x) (and f(x)=Ω(g(x))f(x)= (g(x))) if f(x)≤cg(x)f(x)≤ cg(x) or f(x)≥cg(x)f(x)≥ cg(x) holds for some absolute (if not otherwise specified) constant c>0c>0 respectively. We write f(x)≂g(x)f(x) g(x) (and f(x)=Θ(g(x))f(x)= (g(x))) if f(x)≲g(x)≲f(x)f(x) g(x) f(x). For two vectors u and v in a Hilbert space, we denote their inner product by ⟨,⟩ ,v or ⊤u v. For two matrices A and B of appropriate dimensions, we define their inner product by⟨,⟩:=tr(⊤) ,B := tr(A B). We use ∥⋅∥\|·\| to denote the operator norm for matrices and ℓ2 _2-norm for vectors. For a positive semi-definite (PSD) matrix A and a vector v of appropriate dimension, we write ‖2:=⊤\|v\|_A^2:=v Av. Kronecker/tensor product is denoted by ⊗ . 2 Related Work Synthetic data and model collapse. Training with synthetic data has been shown to induce model collapse (Shumailov et al., 2023), a phenomenon in which model performance degrades significantly when consistently trained with poor-quality synthetic data, as has been empirically verified by a number of subsequent works (Hataya et al., 2023; Bohacek and Farid, 2023; Alemohammad et al., 2023; Guo et al., 2024; Briesch et al., 2024; Shumailov et al., 2024). Theoretical studies of model collapse typically frame it as a regression problem under distribution mismatch between real and synthetic data (Dohmatob et al., 2024c; Dohmatob et al., 2024a; Seddik et al., 2024), showing how synthetic data degrade performance, while others explore how to avoid collapse or leverage synthetic data to benefit learning (Seddik et al., 2024; Firdoussi et al., 2024; Barzilai and Shamir, 2025; Garg et al., 2025). In particular, Dohmatob et al. (2024b) introduce the concept of strong model collapse, where even a small proportion of synthetic data prevents model performance from improving under data scaling, leaving a non-vanishing excess risk floor. Moreover, such collapse cannot be mitigated by simple strategies such as data reweighting ((Jain et al., 2024; Ferbach et al., 2024)) unless the proportion of synthetic data vanishes asymptotically. However, most of them do not capture the dynamics of practical training, where the order and interplay of data sources matter. In this paper, we study finite-sample one-pass SGD and show that while mixed training induces strong model collapse, it can be avoided by a simple two-stage training protocol. Learning with multiple data sources. Learning from multiple data sources has been widely studied in transfer learning, domain adaptation, and continual learning under covariate or model shift (Blitzer et al., 2007; Ben-David et al., 2010; Wang and Schneider, 2015; Lei et al., 2021; Kpotufe and Martinet, 2018; Hanneke and Kpotufe, 2019; Pathak et al., 2022; Li et al., 2023a; Tahir et al., 2024; Li et al., 2025b; Kim et al., 2025). However, most of them focus on statistical estimators rather than optimization dynamics. More recently, several works have analyzed SGD under covariate shift. Wu et al. (2022b) study pretraining-finetuning under covariate shift via SGD and establish excess risk bounds. Ding et al. (2024) extend the analysis to continual learning settings and show how covariate shifts across tasks affect SGD dynamics and forgetting. Liu et al. (2025) characterize the minimax optimality of SGD-type methods for high-dimensional linear regression under covariate shift. Beyond standard covariate shift, Deng and Kpotufe (2025) propose an adaptive mixed-sample SGD that further accounts for model shift while achieving transfer-optimal target risk. In this paper, we consider different data sources with shared covariates but mismatched labeling functions, as commonly arises with synthetic data, and characterize how synthetic data affects finite-sample SGD dynamics. SGD in high-dimensional linear regression. The generalization of stochastic gradient descent (SGD) in linear regression has been extensively studied. In the classical underparameterized regime, a large number of works have studied its behavior (Polyak and Juditsky, 1992; Bach and Moulines, 2013; Défossez and Bach, 2015; Dieuleveut et al., 2017; Jain et al., 2017; Jain et al., 2018; Pillaud-Vivien et al., 2018). One-pass SGD under different iterate schemes (e.g., last iterate, averaging) has also been rigorously studied in the overparameterized setting (Dieuleveut and Bach, 2016; Ge et al., 2019; Berthier et al., 2020; Varre et al., 2021; Wu et al., 2022a; Zou et al., 2023; Zhang et al., 2024; Zhang et al., 2025), providing framework for analyzing the generalization of SGD in high-dimensional regime. In parallel, Several works have extended the analysis to multipass SGD (Lin and Rosasco, 2017; Pillaud-Vivien et al., 2018; Mücke et al., 2019; Lei et al., 2021; Zou et al., 2022). More recently, a line of research has focused on the theoretical scaling laws of SGD in high-dimensional linear regression. Lin et al. (2024) analyzed the last-iterate test error of one-pass SGD in a random sketch model, providing the first systematic derivation of a finite-sample joint scaling law. Lin et al. (2025) further extended their analysis to multi-epoch SGD, followed by an even finer-grained characterization from Yan et al. (2025). Li et al. (2025a) established functional scaling laws and analyzed optimal learning rate and batch size schedules (Li et al., 2026; Wang et al., 2026). Overall, this paper builds on high-dimensional SGD and scaling-law analyses, and extend them from single-source learning to a two-source setting where source mismatch induces new drift and fluctuation effects. 3 Preliminaries Problem Setup. We use ∈ℋx to denote a feature vector, where ℋH is a finite d-dimensional or countably infinite dimensional Hilbert space, and y∈ℝy to denote its label. Linear Regression concerns the following objective: minℛ(),whereℛ():=12(⟨,⟩−y)2 _wR(w),\ where\ R(w):= 12E ( ,w -y )^2 Here ∈ℋw is a weight vector to be learned and the expectation is over (,y)∼(x,y) for some distribution P on ℋ×ℝH×R. Following the prior literature (Bartlett et al., 2020; Wu et al., 2022a; Zou et al., 2023; Lin et al., 2024), we make the following assumptions on the data distribution. Assumption 3.1 (Data Covariance). Let :=[⊤]H:=E[xx ] be the data covariance matrix and assume that tr() tr(H) and all entries of H are finite. In addition, we assume that all parameter vectors involved in the analysis have finite H-norm, i.e., ‖<∞\|w\|_H<∞. Assumption 3.2 (Fourth moment conditions). (1) There exists constant α>0α>0 such that for every PSD matrix A it holds that [⊤⊤]⪯αtr()E[xx Axx ] α tr(HA)H; (2) There is a constant β>0β>0, such that for every PSD matrix A, we have [⊤⊤]−⪰βtr()E[xx Axx ]-HAH β tr(HA)H. In this paper, motivated by (Dohmatob et al., 2024b), we consider a two-source setting, where data are drawn from either a real distribution 1P_1 or a synthetic distribution 2P_2. We assume that both distributions share the same feature marginal ∼x , and differ only in their labeling functions, which allows us to isolate the effect of model mismatch and obtain a clean characterization of the learning dynamics. This formulation naturally captures the prevalent self-training and knowledge distillation paradigms, where a massive pool of unlabeled real-world features is annotated by a potentially biased teacher model to generate supervisory signals. Assumption 3.3 (Well-specified models). For ℓ∈1,2 ∈\1,2\, data from distribution ℓP_ follow a linear model y=⟨,ℓ∗⟩+ξℓ,y= ,w^*_ + _ , where ℓ∗∈ℋw^*_ is the underlying parameter, and the noise ξℓ∼N(0,σℓ2) _ N(0, _ ^2) is independent of x. Define the model shift between the two sources as :=2∗−1∗. δ:=w^*_2-w^*_1. While training may involve samples from both 1P_1 and 2P_2, our goal is to minimize the population risk on the real distribution ℛ():=12(,y)∼1(⟨,⟩−y)2.R(w):= 12E_(x,y) _1( ,w -y)^2. Training protocol. We consider two training protocols for combining real and synthetic data. (1) Mixed training. We are given a dataset consisting of M+NM+N samples, among which M are drawn from the synthetic distribution 2P_2 and N from the real distribution 1P_1. Let p:=M/(M+N)p:=M/(M+N) denote the fraction of synthetic data. The training samples are presented in a random order. An equivalent fixed-budget formulation key to the proof is detailed in Appendix B. (2) Two-stage training. We first train on M synthetic samples drawn from 2P_2, followed by training on N real samples drawn from 1P_1. SGD iteration. We consider training via stochastic gradient descent (SGD) with step sizes γt\ _t\: t+1=t−γt(⟨t,t⟩−yt)t≥1.w_t+1=w_t- _t ( _t,w_t -y_t )x_t t≥ 1. where (t,yt)t=1M+N\(x_t,y_t)\_t=1^M+N follow one of the two protocols above, and γtt=1M+N\ _t\_t=1^M+N are the step sizes. We adopt a well-used geometric decaying stepsize scheduler (Wu et al., 2022a; Wu et al., 2022b; Lin et al., 2024; Lin et al., 2025): for t=1,…,T,γt:=γ/2ℓt, and ℓt=⌊t/(T/log(T))⌋.for t=1,…,T, _t:=γ/2 _t, and _t= t/(T/ (T)) . (3.1) Here the initial stepsize γ is a hyperparameter for the algorithm and T denotes the total number of iterations. In the mixed-training setting, we take T=M+NT=M+N, while in the two-stage setting the same schedule is applied separately within each stage. The output of the algorithm is the last iterate M+Nw_M+N and we study the excess risk on the real distribution: :=ℰ1(M+N)=ℛ(M+N)−ℛ(1∗). Excess:=E_1(w_M+N)=R(w_M+N)-R(w^*_1). Additional Notations. let =∑iλiii⊤H= _i _iv_iv_i be the eigen-decomposition of H, where λii=1∞\ _i\^∞_i=1 are the eigenvalues of H sorted in non-increasing order and iv_i are the corresponding eigenvectors. Following (Zou et al., 2023; Wu et al., 2022b), we denote: 0:k:=∑i=1kλi⊤,k:∞:=∑i>kλi⊤,=∑iδi.H_0:k:= _i=1^k _iv_iv_i , _k:∞:= _i>k _iv_iv_i , δ= _i _iv_i. 4 Main Results 4.1 Mixed training In this subsection, we study SGD trained on a mixture of real and synthetic data. We show that, in this setting, the presence of synthetic data induces a non-vanishing excess risk floor, leading to the phenomenon of strong model collapse Dohmatob et al. (2024b). Detailed proof is deferred to Appendix B and C. Theorem 4.1 (An upper bound for mixed training). Given M synthetic and N real samples. Consider last iterate SGD with stepsize scheme (3.1). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Let N~:=(M+N)/log(M+N) N_ eff:=(M+N)/ (M+N) and σ¯2=(1−p)σ12+pσ22 σ^2=(1-p) _1^2+p _2^2. Suppose γ<1/(4αtr())γ<1/(4α tr(H)). We have [ℰ1(M+N)] [E_1(w_M+N)] ≲∥∏t=1M+N(−γt)(0−1∗)∥2+α∥0−1∗∥20:k∗γN~+k∗:∞D~N~⏟ \| _t=1^M+N(I- _tH)(w_0-w_1^*) \|^2_H+α\|w_0-w_1^*\|^2_ I_0:k^*γ N_ eff+H_k^*:∞ D_ eff N_ eff_ BiasError +σ¯2D~N~⏟+αp‖2D~N~⏟+p2‖(−∏t=1M+N(−γt))‖2⏟, + σ^2 D_ eff N_ eff_ VarError+ α\,p\| δ\|^2_H D_ eff N_ eff_ FlucError+ p^2 \| (I- _t=1^M+N(I- _tH) ) δ \|^2_H_ DriftError, where k∗=maxk:λk≥1γN~k^*= \k: _k≥ 1γ N_ eff\ and D~:=k∗+γ2N~2∑i>k∗λi2. D_ eff:=k^*+γ^2 N_ eff^2 _i>k^* _i^2. Interpretation. This upper bound decomposes the excess risk into four terms: a bias term, a variance term, and two additional terms induced by the presence of synthetic data, which we refer to as the fluctuation error and the drift error. The bias and variance terms match those in classical analyses of SGD with geometrically decaying step sizes (Wu et al., 2022b), capturing the optimization error and stochastic noise, respectively. The fluctuation error is upper bounded by terms involving p(1−p),p2p(1-p),p^2 and p, and is therefore summarized at the order level as (p‖2⋅D~N~).O(p\| δ\|_H^2· D_ eff N_ eff). It arises from the stochastic mixing of the two data sources. Particularly, when ‖2≂σ¯2\| δ\|^2_H σ^2, this term is of the same order as the standard variance term and behaves similarly to it. In contrast, the drift error p2‖(−∏t=1M+N(−γt))‖2p^2\|(I- _t=1^M+N(I- _tH)) δ\|^2_H captures a systematic bias in the optimization trajectory toward the synthetic model, which reflects the accumulation of updates aligned with 2∗w^*_2, leading the iterate to drift away from the target 1∗w^*_1. Together, these two terms reveal a fundamental trade-off: while incorporating synthetic data can increase sample size to reduce bias and variance, it also introduces a bias proportional to the model mismatch δ, which may dominate and degrade performance. Further, when the total sample size grows to infinity with a fixed synthetic proportion p, all terms except the drift term vanish under the condition D~=o(N~) D_ eff=o( N_ eff). As a result, the upper bound converges to a non-zero floor of order p2‖2p^2\| δ\|^2_H. Next, we complement this upper bound with a lower bound, showing that this excess risk floor cannot be avoided in general. Theorem 4.2 (A lower bound for mixed training). Given M synthetic and N real samples. Consider last iterate SGD with stepsize scheme (3.1). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Let N~:=(M+N)/log(M+N) N_ eff:=(M+N)/ (M+N) and σ¯2=(1−p)σ12+pσ22 σ^2=(1-p) _1^2+p _2^2. Suppose M+N≥500M+N≥ 500 and γ<1/λ1γ<1/ _1, then [ℰ1(M+N)] [E_1(w_M+N)] ≳‖∏t=1M+N(−γt)(0−1∗)+p(−∏t=1M+N(−γt))‖2⏟− \| _t=1^M+N(I- _tH)(w_0-w_1^*)+p (I- _t=1^M+N(I- _tH) ) δ \|^2_H_ Non-vanishing\ term +(β∥0−1∗−p∥k∗:∞2+σ¯2)D~N~+p(1−p)γ2N~∥3k∗:∞2 + (β\|w_0-w^*_1-p δ\|_H_k^*:∞^2+ σ^2 ) D_ eff N_ eff+p(1-p)γ^2 N_ eff\| δ\|_H^3_k^*:∞^2 where k∗=maxk:λk≥1γN~k^*= \k: _k≥ 1γ N_ eff\ and D~:=k∗+γ2N~2∑i>k∗λi2. D_ eff:=k^*+γ^2 N_ eff^2 _i>k^* _i^2. Interpretation. The lower bound complements Theorem 4.1 by showing that the excess risk cannot vanish in the large-sample regime. In particular, its first term captures the mean optimization error, jointly determined by the residual initialization error and the drift toward the synthetic model. At finite sample sizes, the spectral filter determines which directions remain dominated by the initialization and which have accumulated synthetic-induced bias. The remaining terms further quantify label-noise variance and mismatch-induced fluctuations, providing a more complete finite-sample characterization. Comparing with the upper bound, we observe that both bounds contain a term governed by the quantity p(−∏t=1M+N(−γt)),p (I- _t=1^M+N(I- _tH) ) δ, which characterizes the drift of the optimization trajectory toward the synthetic model. While the upper bound isolates this effect as the drift error, the lower bound shows that such a term is unavoidable as the initialization is progressively forgotten. Crucially, when the total sample size scales to infinity with a fixed synthetic proportion p, and under the condition D~=o(N~) D_ eff=o( N_ eff), all remaining terms vanish, and both bounds converge to the same non-zero limit of order p2‖2.p^2\| δ\|_H^2. This establishes that the excess risk admits an non-vanishing floor, which is entirely induced by the mismatch between real and synthetic data. This directly leads to the following corollary, which formalizes the phenomenon of strong model collapse(Dohmatob et al., 2024b). Corollary 4.3 (Strong model collapse). Consider last iterate SGD with stepsize scheme (3.1). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Suppose γ<1/4αtr()γ<1/4α tr(H) and D~=o(N~) D_ eff=o( N_ eff). As the total sample size M+NM+N scales to infinity with a fixed synthetic proportion p=M/(M+N)p=M/(M+N), we have limM+N→∞[ℰ1(M+N)]≂p2‖2. _M+N ∞E[E_1(w_M+N)] p^2\| δ\|_H^2. Corollary 4.3 shows that the strong model collapse phenomenon persists under SGD with decaying stepsizes. In particular, even when the total sample size M+NM+N tends to infinity, as long as the proportion of synthetic data p remains non-vanishing, the excess risk does not converge to zero, but instead admits a strictly positive limit of order p2‖2p^2\| δ\|_H^2. This demonstrates that strong model collapse is not merely a property of batch estimators or scaling-law limits, but also arises intrinsically from the optimization dynamics of SGD. Here the condition D~=o(N~) D_ eff=o( N_ eff) requires the effective spectral dimension to grow sublinearly with the effective sample size, and holds for many standard eigendecay regimes; see (Zou et al., 2023; Wu et al., 2022a). Our result should be compared with Dohmatob et al. (2024b), who established strong model collapse for regression models trained on mixtures of real and synthetic data. In their classical linear estimator setting, the error floor is of order p2[‖2∗−1∗‖2].p^2E[\|w_2^*-w_1^*\|_H^2]. Hence, their expression can be viewed as taking an additional expectation over the model mismatch, while ours characterizes the floor for a fixed δ. In this sense, our result is fully consistent with the estimator-level picture of Dohmatob et al. (2024b), but recast in stochastic optimization dynamics. However, a key distinction is that their analysis focuses on classical batch estimators (ridge/OLS) under a proportional high-dimensional asymptotic regime, and the resulting error floor is an estimator-level asymptotic characterization. In contrast, we analyze SGD under finite-sample training and establish non-asymptotic upper and lower bounds. Our corollary is also asymptotic, but it follows from these finite-sample bounds, showing that the same non-vanishing floor arises from the optimization dynamics of SGD. More importantly, we show that strong model collapse is fundamentally induced by sample mixing in the training protocol, rather than an artifact of the optimization scheme. In particular, the same phenomenon arises for constant stepsize SGD with iterate averaging, where γt=γ _t=γ and the output is ¯M+N=1M+N∑i=0M+N−1i w_M+N= 1M+N _i=0^M+N-1w_i. In this case, we have limM+N→∞[ℰ1(¯M+N)]≂p2‖2, _M+N ∞E[E_1( w_M+N)] p^2\| δ\|_H^2, demonstrating that the same non-vanishing error floor persists; see Appendix G for details. This naturally raises an important question: Q: Is strong model collapse inevitable with fixed synthetic data proportion? In the next subsection, we answer this question by showing that it can in fact be avoided. Specifically, we demonstrate that a two-stage training protocol eliminates the non-vanishing error floor, leading to vanishing excess risk, highlighting the importance of how synthetic and real data are scheduled. 4.2 Two-stage Training In this subsection, we study a two-stage training protocol, where the model is first trained on synthetic data and then further trained on real data. We show that, in contrast to the mixed-training setting, this training scheme avoids strong model collapse and leads to vanishing excess risk under mild conditions. The detailed proof is deferred to Appendix D. Theorem 4.4 (Two-stage training). Given M synthetic samples and N real samples. Consider last iterate SGD with stepsize scheme (3.1). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Let N:=N/logN_ eff:=N/ N and ℰ2():=12(,y)∼2[(⟨,⟩−y)2]−12(,y)∼2[(⟨,2∗⟩−y)2]E_2(w):= 12E_(x,y) _2[( ,w -y)^2]- 12E_(x,y) _2[( ,w_2^* -y)^2] denote the excess risk on the synthetic distribution. Suppose γ<1/(4αtr())γ<1/(4α tr(H)). Then we have [ℰ1(M+N)] [E_1(w_M+N)] ≲[ℰ2(M)]+‖∏t=1N(−γt)‖2⏟+(α[ℰ2(M)]+α‖2+σ12)DN⏟. _t=1^N D_ effN_ effE[E_2(w_M)]+ \| _t=1^N(I- _tH) δ \|_H^2_ EffectiveBias+ _t=1^N D_ effN_ eff (α\,E[E_2(w_M)]+α\| δ\|_H^2+ _1^2 ) D_ effN_ eff_ EffectiveVariance. where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. As a direct consequence of Theorem 4.4, by considering the asymptotic regime where the sample size goes to infinity, we can establish that strong model collapse does not occur: Corollary 4.5 (No strong model collapse). Under the assumptions of Theorem 4.4. As the total sample size M+NM+N scales to infinity with a fixed synthetic proportion p=M/(M+N)p=M/(M+N), suppose that D=o(N)D_ eff=o(N_ eff), then [ℰ2(M)]→M→∞0,E[E_2(w_M)] M ∞0, and we have limM+N→∞[ℰ1(M+N)]=0. _M+N→∞E[E_1(w_M+N)]=0. Corollary 4.5 shows that, in sharp contrast to the mixed-training setting, strong model collapse can be avoided under a two-stage training protocol. The key distinction is that synthetic data no longer appears throughout the entire optimization trajectory. Instead, the synthetic stage only determines the initialization of the second stage, after which all subsequent updates are performed on real data. As a result, the synthetic model does not induce a persistent drift in the trajectory. Broadly, this suggests that successfully leveraging synthetic data may hinge crucially on the training curriculum, as the specific ordering and scheduling of data may be just as important as the data itself. Further, Theorem 4.4 sheds light on the mechanism through which synthetic data can be beneficial. The main benefit does not come from reducing variance: the EffectiveVariance term still scales as D/ND_ eff/N_ eff, which is the same order as in the standard single-source SGD bound trained on real data only (Wu et al., 2022b). Instead, the main benefit of synthetic data lies in reducing the bias. First, the synthetic training phase error [ℰ2(M)]E[E_2(w_M)] can be driven to zero as M increases. Second, the residual mismatch term ‖∏t=1N(−γt)‖2\| _t=1^N(I- _tH) δ\|_H^2 replaces the real-only initialization term ‖∏t=1N(−γt)(0−1∗)‖2.\| _t=1^N(I- _tH)(w_0-w_1^*)\|_H^2. Thus, two-stage training effectively replaces random initialization error by the source mismatch =2∗−1∗ δ=w_2^*-w_1^*. When the synthetic data is of high quality so that ‖2\| δ\|_H^2 is small, first-stage places the model substantially closer to the target 1∗w_1^*, yielding a smaller bias after second-stage training. Consequently, synthetic data primarily acts as a bias-reduction mechanism: it improves the initialization for the real-data optimization stage, thereby accelerating convergence and reducing the optimization error when real data is limited. This highlights that the effectiveness of synthetic data critically depends on its quality (i.e., the magnitude of δ), as well as on the training protocol used to incorporate it. In the next section, we make this intuition more precise by analyzing a random sketch model under a power-law eigendecay assumption on the data covariance, which allows us to characterize more explicitly when and how synthetic data improves performance. Learning-rate restart. Throughout the main text, including the scaling-law analysis below, we restart the learning-rate schedule to treat the real-data phase as a separate fine-tuning stage, with a larger cumulative stepsize that accelerates forgetting of the synthetic initialization. Without restart, the mismatch filter becomes ∏t=M+1M+N(−γt), _t=M+1^M+N(I- _tH), which generally forgets δ more slowly. Nevertheless, restart does not affect the qualitative behavior: for fixed p∈(0,1)p∈(0,1), the cumulative real-stage stepsize diverges and the excess risk still has no non-vanishing floor; see detailed discussion in Appendix D.1. 4.3 Random Sketch Model & Scaling laws The results in the previous subsections establish a sharp qualitative distinction between mixed training and two-stage training. However, these results do not yet reveal how the effect of synthetic data interacts with model size. To study this question and establish scaling laws for the two training protocols, we adopt the random sketch model with Gaussian sketch matrix (Lin et al., 2024; Lin et al., 2025; Zhang et al., 2026). Detailed proof of this subsection is deferred to Appendix E. Concretely, following Lin et al. (2024), we introduce a Gaussian random sketch operator :ℋ→ℝDS:H ^D, where entries of S are independently sampled from N(0,1/D)N(0,1/D), and replace the original feature vector x in our problem setup by its D-dimensional sketch Sx. We then train a linear predictor with D trainable parameters in the sketched space, f()=⟨,⟩,∈ℝD,f_v(x)= ,Sx , ^D, starting from zero initialization 0=v_0=0 and using the same two-source data model and training protocols introduced in Section 3. Accordingly, the population risk on the real distribution becomes ℛD():=12(,y)∼1(⟨,⟩−y)2.R_D(v):= 12E_(x,y) _1 ( ,Sx -y )^2. where the expectation is conditioned on S. The corresponding optimal parameters are ℓ∗:=(⊤)−1ℓ∗,ℓ∈1,2.v_ ^*:=(SHS )^-1SHw_ ^*, ∈\1,2\. Thus, the sketch dimension D plays the role of model size, allowing us to study how model scaling interacts with source mismatch and synthetic data. As in Lin et al. (2024), we can decompose the risk into irreducible, approximation, and excess terms: ℛD(M+N)=minℛ(⋅)⏟+minℛD(⋅)−minℛ(⋅)⏟+ℛD(M+N)−minℛD(⋅)⏟.R_D(v_M+N)= (·)_ Irreducible+ _D(·)- (·)_ Approx+ R_D(v_M+N)- _D(·)_ Excess. Following the prior literature (Lin et al., 2024; Lin et al., 2025; Li et al., 2025a; Zhang et al., 2026), we further make the following additional assumptions to establish scaling law upper bounds for the two training protocols and to explicitly characterize the effects of synthetic data and model capacity. Assumption 4.6 (Distributional conditions). We assume ∼(0,)x (0,H), and the random vectors 1∗w_1^* and :=2∗−1∗ δ:=w_2^*-w_1^* satisfy [1∗(1∗)⊤]=E[w_1^*(w_1^*) ]=I, [⊤]=E[ δ δ ]= , and [1∗⊤]=E[w_1^* δ ]=0. Assumption 4.7 (Power-law spectrum and source condition). Let (λi,i)i≥1( _i,v_i)_i≥ 1 be the eigenvalue-eigenvector pairs of H with (λi)i≥1( _i)_i≥ 1 in non-increasing order. We assume λi≂i−a _i i^-a for some a>1a>1, and for the mismatch δ we assume [⟨i,⟩⟨j,⟩]=0E[ _i, δ _j, δ ]=0 for i≠ji≠ j and [λi⟨i,⟩2]≂i−bE[ _i _i, δ ^2] i^-b for some b>1b>1. The exponent b in Assumption 4.7 characterizes how hard the mismatch δ is relative to the data spectrum H. A larger b implies a simpler task, as the mismatch energy decays more rapidly along H’s eigen-directions. Theorem 4.8 (Scaling law for mixed training). Given M synthetic and N real samples. Consider last iterate SGD with stepsize scheme (3.1). Suppose that Assumptions 3.1, 3.2, 3.3, 4.6 and 4.7 hold. Let N~=(M+N)/log(M+N) N_ eff=(M+N)/ (M+N). Suppose γ<1/4αtr(⊤)γ<1/4α tr(SHS ), σ12,σ22≂1 _1^2, _2^2 1 and choose γ≂1γ 1, then with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S, we have ℛD(M+N)−ℛ(1∗) _D(v_M+N)-R(w_1^*) ≲1Da−1+1(N~)1−1/a⏟++(M+N)2(1−1Db−1)⏟ 1D^a-1+ 1( N_ eff)^1-1/a_ Approx+Bias+ ( MM+N )^2(1- 1D^b-1) 1D^a-1+ 1(N_ eff+M_ eff)^1-1/a_ Drift\ Floor Approx+Bias Moreover, if the sketched problem is optimization-saturated that γN~≳Daγ N_ eff D^a, then the preceding upper bound is sharp up to constants: ℛD(M+N)−ℛ(1∗)≂D1−a+N~ 1/a−1+p2(1−D1−b). _D(v_M+N)-R(w_1^*) D^1-a+ N_ eff^\,1/a-1+p^2(1-D^1-b). Interpretation. The first part of Theorem 4.8 provides a clean upper-bound scaling picture for mixed training under the random sketch model. When γN~≳Daγ N_ eff D^a, the cumulative stepsize is large enough to learn even the smallest sketched eigendirection, whose eigenvalue is of order D−aD^-a, and a matching lower bound holds in this optimization-saturated regime. The fluctuation and variance terms are absorbed into higher-order contributions and can be bounded by the same order as the + Approx+Bias term. Hence, they do not affect the leading-order behavior, which is dominated by + Approx+Bias and the drift floor. Second, synthetic data improves the reducible part of the risk at the pre-asymptotic level. Compared with the real-only scaling 1(N)1−1/a 1(N_ eff)^1-1/a in Lin et al. (2024), the term 1(N~)1−1/a 1( N_ eff)^1-1/a shows that synthetic and real samples contribute jointly, effectively increasing the sample size and reducing + Approx+Bias. However, this benefit comes with a fundamental trade-off. While increasing the sketch dimension D decreases the approximation term D1−aD^1-a and therefore improves + Approx+Bias, it simultaneously increases the drift-floor term p2(1−1Db−1).p^2(1- 1D^b-1). Indeed, as D grows, the factor D1−bD^1-b vanishes and the drift contribution approaches its limiting level of order p2p^2. Therefore, larger models fit the real signal better, but they may amplify synthetic-induced degradation by exposing the mismatch between the real and synthetic sources more fully, revealing a fundamental trade-off in mixed training. Theorem 4.9 (Scaling law upper bound for two-stage training). Given M synthetic and N real samples. Consider last iterate SGD with stepsize scheme (3.1). Suppose that Assumptions 3.1, 3.2, 3.3, 4.6 and 4.7 hold. Let N=N/logN_ eff=N/ N and M=M/logM_ eff=M/ M. Suppose γ<1/4αtr(⊤)γ<1/4α tr(SHS ) and σ12,σ22≂1 _1^2, _2^2 1. Choose γ≂1γ 1, and assume b<a+1b<a+1 for simplicity and clarity, then with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S, we have ℛD(M+N)−ℛ(1∗) _D(v_M+N)-R(w^*_1) ≲1Da−1⏟+minD,(N)1/aN⏟ 1D^a-1_ Approx+ \D,(N_ eff)^1/a\N_ eff_ Variance +e−Ω(NDa)⋅max1Da∧b−1,1(M)a∧b−1a+max1Db−1,1(N)b−1a⏟ + e^- ( N_ effD^a)· \ 1D^a b-1, 1(M_ eff) a b-1a \+ \ 1D^b-1, 1(N_ eff) b-1a \_ Bias Interpretation. The variance term is of the same order as in the real-only scaling law of Lin et al. (2024), namely minD,N1/aN, \D,N_ eff^1/a\N_ eff, and therefore synthetic data does not improve the variance, which is consistent with our earlier analysis. The bias term is split into two parts, reflecting two distinct effects of synthetic training phase. The first term, e−Ω(N/Da)maxD1−a∧b,M−(a∧b−1)/a,e^- (N_ eff/D^a) \D^1-a b,M_ eff^-(a b-1)/a\, captures the residual error from the first stage, and is exponentially forgotten during training on real data. The second term, maxD1−b,N−(b−1)/a, \D^1-b,N_ eff^-(b-1)/a\, captures the remaining source mismatch after the second stage real-data training. Compared with the real-only bias upper bound maxD1−a,N−(a−1)/a \D^1-a,N_ eff^-(a-1)/a\(Lin et al., 2024), this term is strictly smaller when b>ab>a, corresponding to high-quality synthetic data. In this regime, the mismatch is easier than the original task, so synthetic training phase places the model closer to the target and reduces the bias more effectively than random initialization. So both bias contributions can be improved at the upper-bound level than in the real-only case. Therefore, unlike mixed training, two-stage training can genuinely benefit from synthetic data: high-quality synthetic data can reduce the bias without introducing any non-vanishing drift floor. This subsection establishes only an upper bound for two-stage training, since a uniform termwise matching lower bound is generally impossible, because upper-bound analysis separately bounds the propagated first-stage error and the source-mismatch contribution, although their cross term may be negative and substantially cancel in the actual final iterate. Accordingly, Theorem 4.9 should be interpreted as an upper-bound scaling characterization of how high-quality synthetic pretraining may reduce bias rather than claiming that it uniformly outperforms real-only training. 4.4 Finite-Sample Conditions for Synthetic Benefit We now ask when the potential synthetic bias-reduction translates into a strict improvement in the actual risk. In this subsection, under a fixed real-data budget and identical real-stage training, we establish an exact finite-sample necessary-and-sufficient condition for two-stage training to outperform real-only training, and then use a one-dimensional toy model to characterize the roles of synthetic sample size, label noise, and source mismatch. Detailed proof is deferred to Appendix F. Theorem 4.10 (Exact condition for strict synthetic benefit). Let Nrealw_N^real and M+Nw_M+N denote the outputs of real-only and two-stage SGD, respectively. Both procedures start from 0w_0; the two-stage procedure additionally uses M synthetic samples before the two procedures undergo the same N real updates, with identical samples, order, and stepsizes. Define ΦN:=∏t=N1(−γttt⊤) _N:= _t=N^1(I- _tx_tx_t ) and N:=1,…,N[ΦN⊤ΦN].K_N:=E_x_1,…,x_N[ _N H _N]. Under Assumptions 3.1 and 3.3, 2([ℰ1(Nreal)]−[ℰ1(M+N)])=‖0−1∗‖N2−‖[M]−1∗‖N2−tr(NCov(M)). 2 (E[E_1(w_N^real)]-E[E_1(w_M+N)] )=\|w_0-w_1^*\|_K_N^2-\|E[w_M]-w_1^*\|_K_N^2- tr (K_NCov(w_M) ). Consequently, two-stage training strictly outperforms real-only training if and only if ‖[M]−1∗‖N2+tr(NCov(M))<‖0−1∗‖N2.\|E[w_M]-w_1^*\|_K_N^2+ tr (K_NCov(w_M) )<\|w_0-w_1^*\|_K_N^2. This condition compares three quantities. ‖0−1∗‖N2\|w_0-w_1^*\|_K_N^2 is the real-only initialization error that remains after N real updates. ‖[M]−1∗‖N2\|E[w_M]-w_1^*\|_K_N^2 is the remaining bias of the synthetic initialization; it includes both finite-M optimization bias and source mismatch. tr(NCov(M)) tr (K_NCov(w_M) ) is the remaining uncertainty generated by finite synthetic training, including synthetic sampling and label noise. The role of NK_N follows from ‖N2=‖ΦN‖2\|u\|_K_N^2=E\| _Nu\|_H^2, which measures the expected error remaining after applying the same N real-data updates to an initial error u. Therefore, synthetic pretraining is beneficial exactly when Mw_M provides a better starting point than 0w_0 after both are followed by the same real-data training and it must improve the initialization in directions that cannot be efficiently corrected by the available real data. Theorem 4.10 gives an exact criterion once the mean and covariance of Mw_M are characterized. In general, however, these quantities depend jointly on the synthetic sample size, label noise, source mismatch, data spectrum, and SGD dynamics, making their individual effects difficult to separate. To disentangle these factors and obtain an explicit sample-size threshold, we next specialize the condition to a one-dimensional toy model. One-dimensional illustration. Let x be Rademacher, independent of the label noises, and consider yreal=θx+ξ1,ysyn=(θ+δ)x+ξ2,[ξℓ]=0,Var(ξℓ)=σℓ2.y^real=θ x+ _1, y^syn=(θ+δ)x+ _2, [ _ ]=0, ( _ )= _ ^2. Starting from w0w_0, the synthetic stage performs M SGD-updates with stepsize 1/s1/s at iteration s. Afterward, both procedures use the same N real samples and stepsizes 0<γt<10< _t<1. Corollary 4.11 (One-dimensional strict-benefit condition). For every M,N≥1M,N≥ 1, [ℰ1(wM+N)]−[ℰ1(wNreal)]=12(∏t=1N(1−γt))2[δ2+σ22M−(w0−θ)2]. [E_1(w_M+N)]-E[E_1(w_N^real)]= 12 ( _t=1^N(1- _t) )^2 [δ^2+ _2^2M-(w_0-θ)^2 ]. Hence, synthetic pretraining strictly improves the expected risk if and only if δ2+σ22/M<(w0−θ)2δ^2+ _2^2/M<(w_0-θ)^2. Such an improvement is achievable for some M if and only if δ2<(w0−θ)2δ^2<(w_0-θ)^2; under this condition, the exact requirement on the synthetic sample size is M>σ22/((w0−θ)2−δ2)M> _2^2/ ((w_0-θ)^2-δ^2 ). This condition displays the bias-variance trade-off directly. The mismatch δ2δ^2 is the error that remains even with infinitely many synthetic samples, whereas σ22/M _2^2/M is the uncertainty caused by finite and noisy synthetic data. Thus, increasing M can eliminate the synthetic variance but cannot compensate for a mismatch larger than the original initialization error. The required synthetic sample size increases linearly with σ22 _2^2 and diverges as δ2δ^2 approaches (w0−θ)2(w_0-θ)^2. Figure 1: Simulation results (i). (a) Excess risk of two protocols in full high-dimensional linear model: mixed training exhibits a clear p-dependent error floor as T scales, while two-stage training continues to improve without saturation. (b) Excess risk at fixed T=50000T=50000 versus p2p^2: under mixing, the excess risk scales approximately linearly with p2p^2, consistent with the theoretical prediction p2‖2p^2\| δ\|_H^2, whereas two-stage training shows a much weaker dependence, matching our theory. (c) Risk increase over p=0p=0 under mixed training in the random sketch model at fixed T=50000T=50000: the degradation grows with model size D and saturates, consistent with the theoretical scaling p2(1−D1−b)p^2(1-D^1-b), indicating that larger models amplify the contamination from synthetic mixed training. 5 Experiments Simulation. We first validate our theory in a Gaussian high-dimensional linear regression setting with dimension d=5000d=5000. Features are drawn from a Gaussian distribution with covariance spectrum λi∝i−1.5 _i i^-1.5, and the synthetic teacher satisfies 2∗=1∗+ w_2^*= w_1^*+ δ with source condition λi[δi2]∝i−b _iE[ _i^2] i^-b. We train last-iterate SGD under the mixed and two-stage protocols in Section 3, and evaluate excess risk on the real distribution. This experiments contain two parts. The first part varies the synthetic proportion p∈(0,1)p∈(0,1), sample size T∈[50,50000]T∈[50,50000], and sketch width D to validate the predicted synthetic-induced error floor. The second part keeps the same high-dimensional setup but varies the synthetic quality b, synthetic sample size M, and synthetic label noise σ22 _2^2 to study when synthetic data helps. All experiments are repeated with 5 random seeds, , and the results of the two parts are reported in Figure 1 and Figure 2, respectively. Figure 2: Simulation results (i). (a) Excess risk of mixed training under different synthetic qualities: low-quality synthetic data induces a larger error floor, while increasing b substantially reduces the degradation. (b) Excess risk of two-stage training under different synthetic qualities: two-stage training continues to improve with T for all b, and higher-quality synthetic data leads to lower final risk. (c) Risk change from adding synthetic samples under a fixed real-data budget: synthetic data helps when b is large, but low-quality synthetic data can eventually hurt as M increases. (d) Benefit region over synthetic quality and label noise: additional synthetic data improves real risk only when the teacher mismatch is sufficiently benign and σ22 _2^2 is not too large. Real-world data. We further extend the experiments to validate our theory on CIFAR10 (Krizhevsky et al., 2009) using a 5-layer CNN and on − WikiText-103 via language model training. Synthetic data is generated by a pretrained teacher model. We observe the same qualitative behavior as predicted by our theory. Full experimental details are provided in Appendix H. 6 Conclusion We studied how synthetic data affects the generalization of one-pass SGD in high-dimensional linear regression with model shift. We establish finite-sample risk bounds for mixed and two-stage training that characterize how source mismatch enters the SGD dynamics under different training protocols. As a consequence, mixed training exhibits strong model collapse while two-stage training avoids this by using synthetic data only in the first stage, showing that collapse is not inevitable through a simple data curriculum. Our random sketch model analysis further establishes scaling-law bounds for both protocols, showing that larger models may amplify synthetic-induced degradation and giving an explicit characterization of how high-quality synthetic pre-training may reduce bias. We also establish an exact finite-sample condition under which two-stage training strictly outperforms real-data-only training with a fixed real-data budget and identical real-stage training. Overall, our results highlight that synthetic data is not inherently harmful or beneficial; its effect depends critically on both its quality and how it is incorporated during training. Limitations. Our analysis focuses on one-pass SGD with shared covariates. Extensions to joint covariate and model shift, as well as multi-epoch SGD, are natural next steps beyond our current setting. In addition, the mixed-training scaling law is tight only in the optimization-saturated regime, while the two-stage result remains an upper bound; this disparity between the two protocols leaves tight characterizations for data-limited and two-stage settings open. References Abdin et al. (2024) M. Abdin, J. Aneja, H. Behl, S. Bubeck, R. Eldan, S. Gunasekar, M. Harrison, R. J. Hewett, M. Javaheripi, P. Kauffmann, et al. Phi-4 technical report. arXiv preprint arXiv:2412.08905. Cited by: §1. Alemohammad et al. (2023) S. Alemohammad, J. Casco-Rodriguez, L. Luzi, A. I. Humayun, H. Babaei, D. LeJeune, A. Siahkoohi, and R. Baraniuk Self-consuming generative models go mad. In The Twelfth International Conference on Learning Representations, Cited by: §1, §2. Bach and Moulines (2013) F. Bach and E. Moulines Non-strongly-convex smooth stochastic approximation with convergence rate o (1/n). Advances in neural information processing systems 26. Cited by: §2. Bartlett et al. (2020) P. L. Bartlett, P. M. Long, G. Lugosi, and A. Tsigler Benign overfitting in linear regression. Proceedings of the National Academy of Sciences 117 (48), p. 30063–30070. Cited by: §3. Barzilai and Shamir (2025) D. Barzilai and O. Shamir When models don’t collapse: on the consistency of iterative mle. arXiv preprint arXiv:2505.19046. Cited by: §1, §2. Ben-David et al. (2010) S. Ben-David, J. Blitzer, K. Crammer, A. Kulesza, F. Pereira, and J. W. Vaughan A theory of learning from different domains. Machine learning 79 (1), p. 151–175. Cited by: §2. Berthier et al. (2020) R. Berthier, F. Bach, and P. Gaillard Tight nonparametric convergence rates for stochastic gradient descent under the noiseless linear model. Advances in Neural Information Processing Systems 33, p. 2576–2586. Cited by: §2. Blitzer et al. (2007) J. Blitzer, K. Crammer, A. Kulesza, F. Pereira, and J. Wortman Learning bounds for domain adaptation. Advances in neural information processing systems 20. Cited by: §2. Bohacek and Farid (2023) M. Bohacek and H. Farid Nepotistically trained generative-ai models collapse. arXiv e-prints, p. arXiv–2311. Cited by: §1, §2. Briesch et al. (2024) M. Briesch, D. Sobania, and F. Rothlauf Large language models suffer from their own output: an analysis of the self-consuming training loop. External Links: 2311.16822, Link Cited by: §1, §2. Brown et al. (2020) T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, et al. Language models are few-shot learners. Advances in neural information processing systems 33, p. 1877–1901. Cited by: §1. Chen et al. (2024) H. Chen, A. Waheed, X. Li, Y. Wang, J. Wang, B. Raj, and M. I. Abdin On the diversity of synthetic data and its impact on training large language models. arXiv preprint arXiv:2410.15226. Cited by: §1. Défossez and Bach (2015) A. Défossez and F. Bach Averaged least-mean-squares: bias-variance trade-offs and optimal sampling distributions. In Artificial Intelligence and Statistics, p. 205–213. Cited by: §2. Deng and Kpotufe (2025) Y. Deng and S. Kpotufe Mixed-sample sgd: an end-to-end analysis of supervised transfer learning. arXiv preprint arXiv:2507.04194. Cited by: §2. Dieuleveut and Bach (2016) A. Dieuleveut and F. Bach Non-parametric stochastic approximation with large step sizes. External Links: 1408.0361, Link Cited by: §2. Dieuleveut et al. (2017) A. Dieuleveut, N. Flammarion, and F. Bach Harder, better, faster, stronger convergence rates for least-squares regression. Journal of Machine Learning Research 18 (101), p. 1–51. Cited by: §2. Ding et al. (2024) M. Ding, K. Ji, D. Wang, and J. Xu Understanding forgetting in continual learning with linear regression. arXiv preprint arXiv:2405.17583. Cited by: §2. Dohmatob et al. (2024a) E. Dohmatob, Y. Feng, and J. Kempe Model collapse demystified: the case of regression. Advances in Neural Information Processing Systems 37, p. 46979–47013. Cited by: §1, §2. Dohmatob et al. (2024b) E. Dohmatob, Y. Feng, A. Subramonian, and J. Kempe Strong model collapse. arXiv preprint arXiv:2410.04840. Cited by: §1, §2, §3, §4.1, §4.1, §4.1, Abstract. Dohmatob et al. (2024c) E. Dohmatob, Y. Feng, P. Yang, F. Charton, and J. Kempe A tale of tails: model collapse as a change of scaling laws. arXiv preprint arXiv:2402.07043. Cited by: §1, §2. Ferbach et al. (2024) D. Ferbach, Q. Bertrand, A. J. Bose, and G. Gidel Self-consuming generative models with curated data provably optimize human preferences. arXiv preprint arXiv:2407.09499. Cited by: §1, §2. Firdoussi et al. (2024) A. E. Firdoussi, M. E. A. Seddik, S. Hayou, R. Alami, A. Alzubaidi, and H. Hacid Maximizing the potential of synthetic data: insights from random matrix theory. arXiv preprint arXiv:2410.08942. Cited by: §1, §2. Garg et al. (2025) A. Garg, S. Bhattacharya, and P. Sur Preventing model collapse under overparametrization: optimal mixing ratios for interpolation learning and ridge regression. arXiv preprint arXiv:2509.22341. Cited by: §1, §2. Ge et al. (2019) R. Ge, S. M. Kakade, R. Kidambi, and P. Netrapalli The step decay schedule: a near optimal, geometrically decaying learning rate procedure for least squares. Advances in neural information processing systems 32. Cited by: §2. Gunasekar et al. (2023) S. Gunasekar, Y. Zhang, J. Aneja, C. C. T. Mendes, A. Del Giorno, S. Gopi, M. Javaheripi, P. Kauffmann, G. de Rosa, O. Saarikivi, et al. Textbooks are all you need. arXiv preprint arXiv:2306.11644. Cited by: §1. Guo et al. (2024) Y. Guo, G. Shang, M. Vazirgiannis, and C. Clavel The curious decline of linguistic diversity: training language models on synthetic text. In Findings of the Association for Computational Linguistics: NAACL 2024, p. 3589–3604. Cited by: §1, §2. Hanneke and Kpotufe (2019) S. Hanneke and S. Kpotufe On the value of target data in transfer learning. Advances in Neural Information Processing Systems 32. Cited by: §2. Hataya et al. (2023) R. Hataya, H. Bao, and H. Arai Will large-scale generative models corrupt future datasets?. In Proceedings of the IEEE/CVF International Conference on Computer Vision, p. 20555–20565. Cited by: §1, §2. Havrilla et al. (2024) A. Havrilla, A. Dai, L. O’Mahony, K. Oostermeijer, V. Zisler, A. Albalak, F. Milo, S. C. Raparthy, K. Gandhi, B. Abbasi, et al. Surveying the effects of quality, diversity, and complexity in synthetic data from large language models. arXiv preprint arXiv:2412.02980. Cited by: §1. Hoffmann et al. (2022) J. Hoffmann, S. Borgeaud, A. Mensch, E. Buchatskaya, T. Cai, E. Rutherford, D. Casas, L. A. Hendricks, J. Welbl, A. Clark, et al. Training compute-optimal large language models. arXiv preprint arXiv:2203.15556 10. Cited by: §1. Jain et al. (2024) A. Jain, A. Montanari, and E. Sasoglu Scaling laws for learning with real and surrogate data. Advances in Neural Information Processing Systems 37, p. 110246–110289. Cited by: §1, §2. Jain et al. (2017) P. Jain, S. M. Kakade, R. Kidambi, P. Netrapalli, V. K. Pillutla, and A. Sidford A markov chain theory approach to characterizing the minimax optimality of stochastic gradient descent (for least squares). arXiv preprint arXiv:1710.09430. Cited by: §2. Jain et al. (2018) P. Jain, S. M. Kakade, R. Kidambi, P. Netrapalli, and A. Sidford Parallelizing stochastic gradient descent for least squares regression: mini-batching, averaging, and model misspecification. Journal of machine learning research 18 (223), p. 1–42. Cited by: §2. Kaplan et al. (2020) J. Kaplan, S. McCandlish, T. Henighan, T. B. Brown, B. Chess, R. Child, S. Gray, A. Radford, J. Wu, and D. Amodei Scaling laws for neural language models. arXiv preprint arXiv:2001.08361. Cited by: §1. Kim et al. (2025) Y. Kim, I. Kim, and S. Park Transfer learning for benign overfitting in high-dimensional linear regression. arXiv preprint arXiv:2510.15337. Cited by: §2. Kpotufe and Martinet (2018) S. Kpotufe and G. Martinet Marginal singularity, and the benefits of labels in covariate-shift. In Conference On Learning Theory, p. 1882–1886. Cited by: §2. Krizhevsky et al. (2009) A. Krizhevsky G. Hinton et al. Learning multiple layers of features from tiny images.(2009). Cited by: §H.1, §5. Lei et al. (2021) Y. Lei, T. Hu, and K. Tang Generalization performance of multi-pass stochastic gradient descent with convex loss functions. Journal of Machine Learning Research 22 (25), p. 1–41. Cited by: §2, §2. Li et al. (2025a) B. Li, F. Chen, Z. Huang, L. Wang, and L. Wu Functional scaling laws in kernel regression: loss dynamics and learning rate schedules. arXiv preprint arXiv:2509.19189. Cited by: §2, §4.3. Li et al. (2026) B. Li, Z. Wang, F. Chen, S. Zhao, R. Zheng, and L. Wu Optimal learning-rate schedules under functional scaling laws: power decay and warmup-stable-decay. arXiv preprint arXiv:2602.06797. Cited by: §2. Li et al. (2023a) H. Li, J. Wu, and V. Braverman Fixed design analysis of regularization-based continual learning. In Conference on lifelong learning agents, p. 513–533. Cited by: §2. Li et al. (2025b) H. Li, J. Wu, and V. Braverman Memory-statistics tradeoff in continual learning with structural regularization. arXiv preprint arXiv:2504.04039. Cited by: §2. Li et al. (2023b) Y. Li, S. Bubeck, R. Eldan, A. Del Giorno, S. Gunasekar, and Y. T. Lee Textbooks are all you need i: phi-1.5 technical report. arXiv preprint arXiv:2309.05463. Cited by: §1. Lin and Rosasco (2017) J. Lin and L. Rosasco Optimal rates for multi-pass stochastic gradient methods. Journal of Machine Learning Research 18 (97), p. 1–47. Cited by: §2. Lin et al. (2025) L. Lin, J. Wu, and P. L. Bartlett Improved scaling laws in linear regression via data reuse. arXiv preprint arXiv:2506.08415. Cited by: §2, §3, §4.3, §4.3. Lin et al. (2024) L. Lin, J. Wu, S. M. Kakade, P. L. Bartlett, and J. D. Lee Scaling laws in linear regression: compute, parameters, and data. Advances in Neural Information Processing Systems 37, p. 60556–60606. Cited by: §E.1, §E.1, §E.1, §E.2, §E.3, §E.3, §E.3, §E.3, §E.3, §E.3, §E.3, Appendix E, Appendix E, §2, §3, §3, §4.3, §4.3, §4.3, §4.3, §4.3, §4.3. Liu et al. (2024) R. Liu, J. Wei, F. Liu, C. Si, Y. Zhang, J. Rao, S. Zheng, D. Peng, D. Yang, D. Zhou, et al. Best practices and lessons learned on synthetic data. arXiv preprint arXiv:2404.07503. Cited by: §1. Liu et al. (2025) Y. Liu, H. Zhang, Q. Chen, and C. Fang Optimal algorithms in linear regression under covariate shift: on the importance of precondition. arXiv preprint arXiv:2502.09047. Cited by: §2. Long et al. (2024) L. Long, R. Wang, R. Xiao, J. Zhao, X. Ding, G. Chen, and H. Wang On llms-driven synthetic data generation, curation, and evaluation: a survey. In Findings of the Association for Computational Linguistics: ACL 2024, p. 11065–11082. Cited by: §1, §1. Maini et al. (2024) P. Maini, S. Seto, R. Bai, D. Grangier, Y. Zhang, and N. Jaitly Rephrasing the web: a recipe for compute and data-efficient language modeling. In Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), p. 14044–14072. Cited by: §1. Mücke et al. (2019) N. Mücke, G. Neu, and L. Rosasco Beating sgd saturation with tail-averaging and minibatching. Advances in Neural Information Processing Systems 32. Cited by: §2. Muennighoff et al. (2023) N. Muennighoff, A. Rush, B. Barak, T. Le Scao, N. Tazi, A. Piktus, S. Pyysalo, T. Wolf, and C. A. Raffel Scaling data-constrained language models. Advances in Neural Information Processing Systems 36, p. 50358–50376. Cited by: §1. Pathak et al. (2022) R. Pathak, C. Ma, and M. Wainwright A new similarity measure for covariate shift with applications to nonparametric regression. In International Conference on Machine Learning, p. 17517–17530. Cited by: §2. Pillaud-Vivien et al. (2018) L. Pillaud-Vivien, A. Rudi, and F. Bach Statistical optimality of stochastic gradient descent on hard learning problems through multiple passes. Advances in Neural Information Processing Systems 31. Cited by: §2. Polyak and Juditsky (1992) B. T. Polyak and A. B. Juditsky Acceleration of stochastic approximation by averaging. SIAM journal on control and optimization 30 (4), p. 838–855. Cited by: §2. Qin et al. (2025) Z. Qin, Q. Dong, X. Zhang, L. Dong, X. Huang, Z. Yang, M. Khademi, D. Zhang, H. H. Awadalla, Y. R. Fung, et al. Scaling laws of synthetic data for language models. arXiv preprint arXiv:2503.19551. Cited by: §1. Radford et al. (2019) A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, I. Sutskever, et al. Language models are unsupervised multitask learners. OpenAI blog 1 (8), p. 9. Cited by: §1. Seddik et al. (2024) M. E. A. Seddik, S. Chen, S. Hayou, P. Youssef, and M. Debbah How bad is training on synthetic data? a statistical analysis of language model collapse. arXiv preprint arXiv:2404.05090. Cited by: §1, §2. Shumailov et al. (2023) I. Shumailov, Z. Shumaylov, Y. Zhao, Y. Gal, N. Papernot, and R. Anderson The curse of recursion: training on generated data makes models forget. arXiv preprint arXiv:2305.17493. Cited by: §1, §2. Shumailov et al. (2024) I. Shumailov, Z. Shumaylov, Y. Zhao, N. Papernot, R. Anderson, and Y. Gal AI models collapse when trained on recursively generated data. Nature 631 (8022), p. 755–759. Cited by: §1, §2. Tahir et al. (2024) J. Tahir, S. Ganguli, and G. M. Rotskoff Features are fate: a theory of transfer learning in high-dimensional regression. arXiv preprint arXiv:2410.08194. Cited by: §2. Varre et al. (2021) A. V. Varre, L. Pillaud-Vivien, and N. Flammarion Last iterate convergence of sgd for least-squares in the interpolation regime.. Advances in Neural Information Processing Systems 34, p. 21581–21591. Cited by: §2. Villalobos et al. (2024) P. Villalobos, A. Ho, J. Sevilla, T. Besiroglu, L. Heim, and M. Hobbhahn Position: will we run out of data? limits of llm scaling based on human-generated data. In Forty-first International Conference on Machine Learning, Cited by: §1. Wang et al. (2026) J. Wang, B. Li, Z. Zhou, M. Wang, Y. Sun, J. Zhang, X. Cai, and L. Wu Fast catch-up, late switching: optimal batch size scheduling via functional scaling laws. arXiv preprint arXiv:2602.14208. Cited by: §2. Wang and Schneider (2015) X. Wang and J. G. Schneider Generalization bounds for transfer learning under model shift.. In UAI, p. 922–931. Cited by: §2. Wu et al. (2022a) J. Wu, D. Zou, V. Braverman, Q. Gu, and S. Kakade Last iterate risk bounds of sgd with decaying stepsize for overparameterized linear regression. In International conference on machine learning, p. 24280–24314. Cited by: §B.2, §B.3, §B.3, §C.2, §C.2, Appendix D, §2, §3, §3, §4.1. Wu et al. (2022b) J. Wu, D. Zou, V. Braverman, Q. Gu, and S. Kakade The power and limitation of pretraining-finetuning for linear regression under covariate shift. Advances in Neural Information Processing Systems 35, p. 33041–33053. Cited by: §B.2, §2, §3, §3, §4.1, §4.2. Yan et al. (2025) T. Yan, H. Wen, B. Li, K. Luo, W. Chen, and K. Lyu Larger datasets can be repeated more: a theoretical analysis of multi-epoch scaling in linear regression. arXiv preprint arXiv:2511.13421. Cited by: §2. Zhang et al. (2025) D. Zhang, J. Su, and D. Zou Learning under quantization for high-dimensional linear regression. arXiv preprint arXiv:2510.18259. Cited by: §2. Zhang et al. (2026) D. Zhang, X. Tang, Y. Liang, and D. Zou Scaling laws for precision in high-dimensional linear regression. arXiv preprint arXiv:2602.19241. Cited by: §4.3, §4.3. Zhang et al. (2024) H. Zhang, Y. Liu, Q. Chen, and C. Fang The optimality of (accelerated) sgd for high-dimensional quadratic optimization. arXiv preprint arXiv:2409.09745. Cited by: §2. Zou et al. (2023) D. Zou, J. Wu, V. Braverman, Q. Gu, and S. M. Kakade Benign overfitting of constant-stepsize sgd for linear regression. Journal of Machine Learning Research 24 (326), p. 1–58. Cited by: Appendix A, Appendix A, Appendix A, Lemma A.1, §G.1, §G.1, §G.1, §G.1, §2, §3, §3, §4.1. Zou et al. (2022) D. Zou, J. Wu, V. Braverman, Q. Gu, and S. Kakade Risk bounds of multi-pass sgd for least squares in the interpolation regime. Advances in Neural Information Processing Systems 35, p. 12909–12920. Cited by: §2. Appendix Contents 1 Introduction 2 Related Work 3 Preliminaries 4 Main Results 4.1 Mixed training 4.2 Two-stage Training 4.3 Random Sketch Model & Scaling laws 4.4 Finite-Sample Conditions for Synthetic Benefit 5 Experiments 6 Conclusion References A Technical Operator Preliminaries B Mixed Training Upper Bound Analysis B.1 Excess Risk Decomposition B.2 Bias upper bound B.3 Variance upper bound B.4 Drift upper bound B.5 Fluctuation upper bound B.6 Proof of Theorem C Mixed Training Lower Bound Analysis C.1 The mean recursion C.2 The centered fluctuation term C.3 Proof of Theorem C.4 Proof of Corollary D Two-stage Training Analysis D.1 Two-stage training without learning-rate restart E Random Sketch Model E.1 Mixed Training Upper Bound E.2 Mixed Training Lower Bound E.3 Two-stage Training Analysis F Exact Synthetic-Benefit Condition G constant stepsize SGD with iterate averaging G.1 Upper bound analysis G.2 Lower bound analysis G.3 Strong model collapse behavior H Additional Experiments H.1 CIFAR10 CNN Experiments H.2 Language Model Experiments Appendix Contents 1 Introduction 2 Related Work 3 Preliminaries 4 Main Results 4.1 Mixed training 4.2 Two-stage Training 4.3 Random Sketch Model & Scaling laws 4.4 Finite-Sample Conditions for Synthetic Benefit 5 Experiments 6 Conclusion References A Technical Operator Preliminaries B Mixed Training Upper Bound Analysis B.1 Excess Risk Decomposition B.2 Bias upper bound B.3 Variance upper bound B.4 Drift upper bound B.5 Fluctuation upper bound B.6 Proof of Theorem C Mixed Training Lower Bound Analysis C.1 The mean recursion C.2 The centered fluctuation term C.3 Proof of Theorem C.4 Proof of Corollary D Two-stage Training Analysis D.1 Two-stage training without learning-rate restart E Random Sketch Model E.1 Mixed Training Upper Bound E.2 Mixed Training Lower Bound E.3 Two-stage Training Analysis F Exact Synthetic-Benefit Condition G constant stepsize SGD with iterate averaging G.1 Upper bound analysis G.2 Lower bound analysis G.3 Strong model collapse behavior H Additional Experiments H.1 CIFAR10 CNN Experiments H.2 Language Model Experiments Appendix A Technical Operator Preliminaries Operators. We first summarize the linear operators (on symmetric matrices) to be used in the proof: ℐ=⊗,ℳ=[⊗],ℳ~=⊗, =I , =E[x ], M=H , t=⊗+⊗−γtℳ,~t=⊗+⊗−γt⊗. T_t=H +I - _tM, T_t=H +I - _tH . We use the notation ∘O to denotes the operator O acting on a symmetric matrix A. For example, with these definitions, we have that for a symmetric matrix A, [Zou et al., 2023]: ℐ∘=,ℳ∘=[(⊤)⊤],ℳ~∘=, =A,\ \ \ M =E[(x Ax)xx ],\ \ \ M =HAH, (ℐ−γtt)∘=[(−γt⊤)(−γt⊤)],(ℐ−γt~t)∘=(−γt)(−γt). (I- _tT_t) =E[(I- _txx )A(I- _txx )],\ \ (I- _t T_t) =(I- _tH)A(I- _tH). (A.1) For the linear operators we have the following technical lemma from Zou et al. [2023]. Lemma A.1 (Lemma B.1, Zou et al. [2023]). An operator O defined on symmetric matrices is called PSD mapping, if ⪰0A 0 implies ∘⪰0O 0. Then we have 1. ℳM and ℳ~ M are both PSD mappings. 2. ℐ−γtI-γT_t and ℐ−γ~tI-γ T_t are both PSD mappings. 3. ℳ−ℳ~M- M and ~t−t T_t-T_t are both PSD mappings. 4. If 0<γt<1/λ10< _t<1/ _1, then t~−1 T_t^-1 exists, and is a PSD mapping. 5. If 0<γt<1/(αtr())0< _t<1/(α tr(H)), then t−1∘T_t^-1 exists for PSD matrix A, and t−1T_t^-1 is a PSD mapping. Proof. See proof of Lemma B.1 in Zou et al. [2023]. ∎ Appendix B Mixed Training Upper Bound Analysis Local notation. For notational simplicity, throughout this section we depart slightly from the notation in the main text and use N to denote the total sample size. Thus, the numbers of synthetic and real samples are pNpN and (1−p)N(1-p)N, respectively, where p∈[0,1]p∈[0,1] is the synthetic-data proportion. We assume that pNpN and (1−p)N(1-p)N are integers; replacing them by the corresponding floor or ceiling does not affect any of our results. Fixed-budget sampling model. We formalize the random-order mixed-training protocol used in our analysis. Let N,p:==(z1,…,zN)∈0,1N:∑t=1Nzt=pN,Z_N,p:= \z=(z_1,…,z_N)∈\0,1\^N: _t=1^Nz_t=pN \, where zt=1z_t=1 indicates that the sample at iteration t is synthetic and zt=0z_t=0 indicates that it is real. We first draw the source schedule z uniformly from N,pZ_N,p. Conditional on the complete source schedule z, the training samples are generated according to t∼i.i.d.,yt=⟨t,1∗+zt⟩+ξt,t=1,…,N,x_t .i.d. D, y_t= _t,w_1^*+z_t δ + _t, t=1,…,N, where :=2∗−1∗ δ:=w_2^*-w_1^* and ξt|zt∼N(0,σ12),zt=0,N(0,σ22),zt=1. _t z_t casesN(0, _1^2),&z_t=0,\\ N(0, _2^2),&z_t=1. cases Conditional on z, the noises ξtt=1N\ _t\_t=1^N are independent across iterations and independent of all covariates. This construction is equivalent in distribution to independently drawing pNpN complete samples from 2P_2 and (1−p)N(1-p)N complete samples from 1P_1, and then uniformly shuffling their union. Indeed, a uniform shuffle induces a source schedule that is uniform over N,pZ_N,p. Conditional on any such schedule, the observations at positions with zt=1z_t=1 are independent samples from 2P_2, while those at positions with zt=0z_t=0 are independent samples from 1P_1. Since samples within each source are exchangeable, the two constructions induce the same joint distribution over the ordered training sequence. Our analysis does not assume independence of the source indicators. To identify the conditional independence used in the SGD recursions, define the filtration ℱt−1:=σ(,1,ξ1,…,t−1,ξt−1).F_t-1:=σ\! (z,x_1, _1,…,x_t-1, _t-1 ). Since the real and synthetic distributions share the same covariate marginal D, the current covariate tx_t is independent of ℱt−1F_t-1 and satisfies [tt⊤∣ℱt−1]=.E[x_tx_t _t-1]=H. Moreover, the sample-wise noise construction gives [ξt∣ℱt−1,t]=0.E[ _t _t-1,x_t]=0. Consequently, the noise cross terms appearing in the bias–variance recursions vanish. Averaging over the source schedule further yields [ξt2tt⊤]=σ¯2,σ¯2:=(1−p)σ12+pσ22.E[ _t^2x_tx_t ]= σ^2H, σ^2:=(1-p) _1^2+p _2^2. Therefore, the covariance and noise recursions used below remain exact under the fixed-budget sampling protocol. They rely on fresh covariates and sample-wise conditionally independent noises, rather than on independence of the source indicators. B.1 Excess Risk Decomposition Starting from 0∈ℝℋw_0 ^H, SGD on the squared loss performs t=t−1+γt(yt−t⊤t−1)t,t≥1.w_t=w_t-1+ _t(y_t-x_t w_t-1)x_t, t≥ 1. Substituting the label model yields t=t−1+γt(t⊤(1∗+zt−t−1)+ξt)t.w_t=w_t-1+ _t\! (x_t (w_1^*+z_t δ-w_t-1)+ _t )x_t. Let t:=t−1∗ η_t:=w_t-w_1^* denote the parameter error relative to the real-data optimum. Then equation (B.1) serves as the starting point of the analysis: t=(−γttt⊤)t−1+γtzttt⊤+γtξtt. η_t=(I- _tx_tx_t ) η_t-1+ _tz_tx_tx_t δ+ _t _tx_t. (B.1) With as light abuse of probability spaces, one can view the centered SGD iterates as the sum of four random processes: t=biast+vart+driftt+fluctt,t=1,2,…N. η_t= η^bias_t+ η^var_t+ η^drift_t+ η^fluct_t, t=1,2,...N. where tbias=(−γttt⊤)t−1bias;0bias=0−1∗, cases η^bias_t=(I- _tx_tx_t ) η^bias_t-1;\\ η^bias_0=w_0-w_1^*, cases tvar=(−γttt⊤)t−1var+γtξtt;0var=, cases η^var_t=(I- _tx_tx_t ) η^var_t-1+ _t _tx_t;\\ η^var_0= 0, cases tdrift=(−γttt⊤)t−1drift+γtp;0drift=, cases η^drift_t=(I- _tx_tx_t ) η^drift_t-1+ _tpH δ;\\ η^drift_0= 0, cases tfluct=(−γttt⊤)t−1fluct+γtzttt⊤−γtp;0fluct=, cases η^fluct_t=(I- _tx_tx_t ) η^fluct_t-1+ _tz_tx_tx_t δ- _tpH δ;\\ η^fluct_0= 0, cases Moreover, we denote t=[tbias⊗tbias],t=[tvar⊗tvar], _t=E[ η^bias_t η^bias_t], _t=E[ η^var_t η^var_t], t=[tdrift⊗tdrift],t=[tfluct⊗tfluct]. _t=E[ η^drift_t η^drift_t], _t=E[ η^fluct_t η^fluct_t]. Using the definition of t,t,t,tB_t,V_t,D_t,F_t, we can then decompose Excess Risk using Cauchy–Schwarz Inequality: [ℛ(N)−ℛ(1∗)] [R(w_N)-R(w_1^*)] =12⟨,[N⊗N]⟩ = 12 ,E[ η_N η_N] =12⟨,[(Nbias+Nvar+Ndrift+Nfluct)⊗(Nbias+Nvar+Ndrift+Nfluct)]⟩ = 12 ,E[( η_N^bias+ η_N^var+ η_N^drift+ η_N^fluct) ( η_N^bias+ η_N^var+ η_N^drift+ η_N^fluct)] ≤2⟨,[Nbias⊗Nbias]⟩+2⟨,[Nvar⊗Nvar]⟩ ≤ 2 ,E[ η_N^bias η_N^bias] +2 ,E[ η_N^var η_N^var] +2⟨,[Ndrift⊗Ndrift]⟩+2⟨,[Nfluct⊗Nfluct]⟩ +2 ,E[ η_N^drift η_N^drift] +2 ,E[ η_N^fluct η_N^fluct] =2⟨,N⟩+2⟨,N⟩+2⟨,N⟩+2⟨,N⟩. =2 ,B_N +2 ,V_N +2 ,D_N +2 ,F_N . B.2 Bias upper bound Using the defined operators, the update rule of the iterates imply the following recursive form of tB_t: t=(ℐ−γtt)∘t−1,0=0⊗0,B_t=(I- _tT_t) _t-1, _0= η_0 η_0, (B.2) Note that this bias term ⟨,N⟩ ,B_N is the same as that in Wu et al. [2022a], so we can apply Lemma in Wu et al. [2022a] to bound it directly: Lemma B.1 (A bias upper bound). Suppose Assumptions 3.1 and 3.2 hold. Let N=N/logN_ eff=N/ N. Consider (B.2). Suppose γ<1/(4αtr())γ<1/(4α tr(H)). We have ⟨,N⟩≲∥∏t=1N(−γt)(0−∗)∥2+α∥0−∗∥0:k∗γN+k∗:∞2⋅DN ,B_N \| _t=1^N(I- _tH)(w_0-w^*) \|^2_H+α\|w_0-w^*\|^2_ I_0:k^*γ N_ eff+H_k^*:∞· D_ effN_ eff where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. See proof of Theorem D.1 in Wu et al. [2022b]. ∎ B.3 Variance upper bound Using the defined operators, the update rule of the iterates imply the following recursive form of tV_t: t=(ℐ−γtt)∘t−1+γt2t,0=.V_t=(I- _tT_t) _t-1+ _t^2 _t, _0= 0. (B.3) where t=[ξt2tt⊤] _t=E[ _t^2x_tx_t ]. Note that t\x_t\ are independent of zt\z_t\. So we have: t=[[ξt2tt⊤]|zt]=(1−p)1+p2⪯((1−p)σ12+pσ22). _t=E[E[ _t^2x_tx_t ]|z_t]=(1-p) _1+p _2 ((1-p) _1^2+p _2^2)H. Then we can apply Lemma C.2 in Wu et al. [2022a] to bound variance term ⟨,N⟩ ,V_N by simply replace σ2σ^2 with (1−p)σ12+pσ22(1-p) _1^2+p _2^2. Lemma B.2 (A variance bound). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Consider (B.3). Let N=N/logN_ eff=N/ N. Suppose γ<1/(αtr())γ<1/(α tr(H)). We have ⟨,N⟩≤8σ21−γαtr()⋅DN, ,V_N ≤ 8σ^21-γα tr(H)· D_ effN_ eff, where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. See proof of Theorem C.2 in Wu et al. [2022a]. ∎ B.4 Drift upper bound We now analyze the drift term ⟨,N⟩=⟨,[Ndrift⊗Ndrift]⟩. ,D_N = ,E[ η^drift_N η^drift_N] . Lemma B.3 (A drift bound). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Consider (B.3). Let N=N/logN_ eff=N/ N. Suppose γ<1/(αtr())γ<1/(α tr(H)). We have ⟨,N⟩≤p2‖(−∏t=1N(−γt))‖2+8αp2‖21−γαtr()⋅DN, ,D_N ≤ p^2 \| (I- _t=1^N(I- _tH) ) δ \|_H^2+ 8α p^2\| δ\|_H^21-γα tr(H)· D_ effN_ eff, where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. Recall that tdrift η^drift_t satisfies tdrift=(−γttt⊤)t−1drift+γtp,0drift=. η^drift_t=(I- _tx_tx_t ) η^drift_t-1+ _tpH δ, η^drift_0= 0. We first center tdrift η^drift_t by separating its mean and fluctuation parts. Define ¯tdrift:=[tdrift]. η^drift_t:=E[ η^drift_t]. Since tx_t is independent of t−1drift η^drift_t-1 and [tt⊤]=E[x_tx_t ]=H, taking expectation on both sides yields ¯tdrift=(−γt)¯t−1drift+γtp,¯0drift=. η^drift_t=(I- _tH) η^drift_t-1+ _tpH δ, η^drift_0= 0. (B.4) Define the centered drift fluctuation ~tdrift:=tdrift−¯tdrift η^drift_t:= η^drift_t- η^drift_t and subtracting (B.4) from the recursion of tdrift η^drift_t, we obtain ~tdrift η^drift_t =(−γttt⊤)t−1drift−(−γt)¯t−1drift =(I- _tx_tx_t ) η^drift_t-1-(I- _tH) η^drift_t-1 =(−γttt⊤)~t−1drift+γt(−tt⊤)¯t−1drift,~0drift=. =(I- _tx_tx_t ) η^drift_t-1+ _t(H-x_tx_t ) η^drift_t-1, η^drift_0= 0. (B.5) By definition, [~tdrift]=,∀t≥1.E[ η^drift_t]= 0, ∀ t≥ 1. Using tdrift=¯tdrift+~tdrift, η^drift_t= η^drift_t+ η^drift_t, we can expand tD_t as t _t =[(¯tdrift+~tdrift)⊗(¯tdrift+~tdrift)] =E\! [( η^drift_t+ η^drift_t) ( η^drift_t+ η^drift_t) ] =¯tdrift⊗¯tdrift+[~tdrift⊗~tdrift]+¯tdrift⊗[~tdrift]+[~tdrift]⊗¯tdrift. = η^drift_t η^drift_t+E[ η^drift_t η^drift_t]+ η^drift_t [ η^drift_t]+E[ η^drift_t] η^drift_t. Since [~tdrift]=E[ η^drift_t]= 0, the two cross terms vanish. Hence t=¯tdrift⊗¯tdrift+~t,~t:=[~tdrift⊗~tdrift].D_t= η^drift_t η^drift_t+ D_t, D_t:=E[ η^drift_t η^drift_t]. Therefore the drift contribution to the excess risk can be decomposed into two parts: ⟨,N⟩=‖¯Ndrift‖2+⟨,~N⟩. ,D_N = \| η^drift_N \|_H^2+ , D_N . (B.6) In the following, we will bound the two terms in (B.6) separately. Mean part of the drift term. We now derive the explicit form of ‖¯Ndrift‖2 \| η^drift_N \|_H^2. Recall from (B.4) that ¯tdrift=(−γt)¯t−1drift+γtp,¯0drift=. η^drift_t=(I- _tH) η^drift_t-1+ _tpH δ, η^drift_0= 0. Unrolling the recursion gives ¯Ndrift η^drift_N =p∑s=1Nγs(∏t=s+1N(−γt)). =p _s=1^N _s ( _t=s+1^N(I- _tH) )H δ. (B.7) Since each factor (−γt)(I- _tH) is a polynomial in H, it commutes with H. Hence (∏t=s+1N(−γt))=(∏t=s+1N(−γt)), ( _t=s+1^N(I- _tH) )H=H ( _t=s+1^N(I- _tH) ), and therefore ∑s=1Nγs(∏t=s+1N(−γt)) _s=1^N _s ( _t=s+1^N(I- _tH) )H =−∏t=1N(−γt). =I- _t=1^N(I- _tH). Substituting this identity into (B.7), we obtain ¯Ndrift=p(−∏t=1N(−γt)). η^drift_N=p (I- _t=1^N(I- _tH) ) δ. As a consequence, ‖¯Ndrift‖2=p2‖(−∏t=1N(−γt))‖2. \| η^drift_N \|_H^2=p^2 \| (I- _t=1^N(I- _tH) ) δ \|_H^2. (B.8) This is exactly the first drift term appearing in the Lemma. Fluctuation part of the drift term. We now bound the second term in (B.6), namely ⟨,~N⟩,~t:=[~tdrift⊗~tdrift]. , D_N , D_t:=E[ η^drift_t η^drift_t]. Recall from (B.5) that ~tdrift=(−γttt⊤)~t−1drift+γt(−tt⊤)¯t−1drift,~0drift=. η^drift_t=(I- _tx_tx_t ) η^drift_t-1+ _t(H-x_tx_t ) η^drift_t-1, η^drift_0= 0. Define ~tdrift:=[((−tt⊤)¯t−1drift)⊗((−tt⊤)¯t−1drift)]. ^drift_t:=E\! [ ((H-x_tx_t ) η^drift_t-1 ) ((H-x_tx_t ) η^drift_t-1 ) ]. Since ¯t−1drift η^drift_t-1 is deterministic, the recursion for ~t D_t takes the form ~t=(ℐ−γtt)∘~t−1+γt2~tdrift,~0=. D_t=(I- _tT_t) D_t-1+ _t^2 ^drift_t, D_0= 0. (B.9) Indeed, the cross term vanishes because [(−tt⊤)¯t−1drift]=(−[tt⊤])¯t−1drift=.E [(H-x_tx_t ) η^drift_t-1 ]=(H-E[x_tx_t ]) η^drift_t-1= 0. Let t:=¯t−1drift⊗¯t−1drift,A_t:= η^drift_t-1 η^drift_t-1, then ~tdrift ^drift_t =[(−tt⊤)t(−tt⊤)] =E\! [(H-x_tx_t )A_t(H-x_tx_t ) ] =[tt⊤ttt⊤]−t =E[x_tx_t A_tx_tx_t ]-HA_tH =(ℳ−ℳ~)∘t. =(M- M) _t. Since ℳ−ℳ~M- M is a PSD mapping by Lemma A.1, we have ~tdrift⪰0 ^drift_t 0. On the other hand, by Assumption 3.2, [tt⊤ttt⊤]=ℳ∘t⪯αtr(t)=α‖¯t−1drift‖2.E[x_tx_t A_tx_tx_t ]=M _t α\, tr(HA_t)\,H=α\,\| η^drift_t-1\|_H^2\,H. Therefore, ~tdrift⪯α‖¯t−1drift‖2. ^drift_t α\,\| η^drift_t-1\|_H^2\,H. (B.10) Next, using (B.4), ¯t−1drift=p(−∏s=1t−1(−γs)), η^drift_t-1=p (I- _s=1^t-1(I- _sH) ) δ, hence ‖¯t−1drift‖2 \| η^drift_t-1\|_H^2 =p2∑iλi(1−∏s=1t−1(1−γsλi))2δi2 =p^2 _i _i (1- _s=1^t-1(1- _s _i) )^2 _i^2 ≤p2∑iλiδi2=p2‖2. ≤ p^2 _i _i _i^2=p^2\| δ\|_H^2. (B.11) Combining (B.10) and (B.11), we obtain ~tdrift⪯αp2‖2. ^drift_t α p^2\| δ\|_H^2H. (B.12) The recursion (B.9) has exactly the same form as the variance recursion (B.3). Therefore, applying Lemma B.2 with σ2←αp2‖2,σ^2\;←\;α p^2\| δ\|_H^2, we immediately get ⟨,~N⟩≤8αp2‖21−γαtr()⋅DN, , D_N ≤ 8α p^2\| δ\|_H^21-γα tr(H)· D_ effN_ eff, (B.13) where k∗=maxk:λk≥1γN,D:=k∗+γ2N2∑i>k∗λi2.k^*= \k: _k≥ 1γ N_ eff \, D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Combining (B.8) and (B.13), we conclude that ⟨,N⟩=‖¯Ndrift‖2+⟨,~N⟩≤p2‖(−∏t=1N(−γt))‖2+8αp2‖21−γαtr()⋅DN. ,D_N = \| η^drift_N \|_H^2+ , D_N ≤ p^2 \| (I- _t=1^N(I- _tH) ) δ \|_H^2+ 8α p^2\| δ\|_H^21-γα tr(H)· D_ effN_ eff. ∎ B.5 Fluctuation upper bound We now analyze the drift term ⟨,N⟩=⟨,[Nfluct⊗Nfluct]⟩. ,F_N = ,E[ η^fluct_N η^fluct_N] . Lemma B.4 (A fluctuation bound). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Consider (B.3). Let N=N/logN_ eff=N/ N. Suppose γ<1/(αtr())γ<1/(α tr(H)). We have ⟨,N⟩≤8Np(1−p)‖2(N−1)(1−αγtr())⋅DN+16αp‖21−γαtr()⋅DN, ,F_N ≤ 8Np(1-p)\| δ\|_H^2(N-1) (1-αγ tr(H) )· D_ effN_ eff+ 16α p\| δ\|_H^21-γα tr(H)· D_ effN_ eff, where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. Recall that tfluct η^fluct_t satisfies tfluct=(−γttt⊤)t−1fluct+γtzttt⊤−γtp,0fluct=. η^fluct_t=(I- _tx_tx_t ) η^fluct_t-1+ _tz_tx_tx_t δ- _tpH δ, η^fluct_0= 0. We decompose the driving term as zttt⊤−p=zt(tt⊤−)+(zt−p).z_tx_tx_t δ-pH δ=z_t(x_tx_t -H) δ+(z_t-p)H δ. Accordingly, define two auxiliary processes tfluct=tfluct,1+tfluct,2,t≥0, η^fluct_t= η^fluct,1_t+ η^fluct,2_t, t≥ 0, where tfluct,1=(−γttt⊤)t−1fluct,1+γtzt(tt⊤−),0fluct,1=,tfluct,2=(−γttt⊤)t−1fluct,2+γt(zt−p),0fluct,2=. cases η^fluct,1_t=(I- _tx_tx_t ) η^fluct,1_t-1+ _tz_t(x_tx_t -H) δ,\\[2.15277pt] η^fluct,1_0= 0, cases cases η^fluct,2_t=(I- _tx_tx_t ) η^fluct,2_t-1+ _t(z_t-p)H δ,\\[2.15277pt] η^fluct,2_0= 0. cases Then by Cauchy–Schwarz, ⟨,N⟩ ,F_N =⟨,[(Nfluct,1+Nfluct,2)⊗(Nfluct,1+Nfluct,2)]⟩ = ,E [( η^fluct,1_N+ η^fluct,2_N) ( η^fluct,1_N+ η^fluct,2_N) ] ≤2⟨,N(1)⟩+2⟨,N(2)⟩, ≤ 2 ,F_N^(1) +2 ,F_N^(2) , (B.14) where t(1):=[tfluct,1⊗tfluct,1],t(2):=[tfluct,2⊗tfluct,2].F_t^(1):=E[ η^fluct,1_t η^fluct,1_t], _t^(2):=E[ η^fluct,2_t η^fluct,2_t]. We first analyze the term t(1)F_t^(1). Expanding the outer product gives t(1) _t^(1) =[((−γt⊤)t−1fluct,1+γtzt(t⊤−)) =E [ ((I- _tx_tx_t ) η^fluct,1_t-1+ _tz_t(x_tx_t -H) δ ) ⊗((−γt⊤)t−1fluct,1+γtzt(t⊤−))] ((I- _tx_tx_t ) η^fluct,1_t-1+ _tz_t(x_tx_t -H) δ ) ] =[(−γttt⊤)t−1fluct,1⊗(−γttt⊤)t−1fluct,1] =E [(I- _tx_tx_t ) η^fluct,1_t-1 (I- _tx_tx_t ) η^fluct,1_t-1 ] +γt2[zt2((tt⊤−))⊗((tt⊤−))] + _t^2E [z_t^2 ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ] +γt[(−γttt⊤)t−1fluct,1⊗zt(tt⊤−)] + _tE [(I- _tx_tx_t ) η^fluct,1_t-1 z_t(x_tx_t -H) δ ] +γt[zt(tt⊤−)⊗(−γttt⊤)t−1fluct,1]. + _tE [z_t(x_tx_t -H) δ (I- _tx_tx_t ) η^fluct,1_t-1 ]. Define tfluct,1:=[zt2((tt⊤−))⊗((tt⊤−))]. _t^fluct,1:=E [z_t^2 ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ]. (B.15) Then t(1) _t^(1) =(ℐ−γtt)∘t−1(1)+γt2tfluct,1+t(1)+(t(1))⊤, =(I- _tT_t) _t-1^(1)+ _t^2 _t^fluct,1+C_t^(1)+(C_t^(1)) , (B.16) where t(1):=γt[(−γttt⊤)t−1fluct,1⊗zt(tt⊤−)].C_t^(1):= _tE [(I- _tx_tx_t ) η^fluct,1_t-1 z_t(x_tx_t -H) δ ]. We now prove that t(1)=C_t^(1)= 0. Let :=σ(z1,…,zN)Z:=σ(z_1,…,z_N) be the sigma-field generated by the entire source schedule. We first claim that [tfluct,1∣]=,∀t≥0.E [ η^fluct,1_t ]= 0, ∀ t≥ 0. (B.17) Indeed, the claim is trivial for t=0t=0. Suppose it holds at time t−1t-1. Using the recursion tfluct,1=(−γttt⊤)t−1fluct,1+γtzt(tt⊤−), η^fluct,1_t=(I- _tx_tx_t ) η^fluct,1_t-1+ _tz_t(x_tx_t -H) δ, and using the independence of tx_t from ∨σ(1,…,t−1)Z σ(x_1,…,x_t-1), we have [tfluct,1∣] [ η^fluct,1_t ] =[[(−γt⊤)t−1fluct,1+γtzt(t⊤−)∣,1,…,t−1]∣] =E [E [(I- _tx_tx_t ) η^fluct,1_t-1+ _tz_t(x_tx_t -H) δ ,x_1,…,x_t-1 ] ] =[(−γt)t−1fluct,1+γtzt[(tt⊤−)]∣] =E [(I- _tH) η^fluct,1_t-1+ _tz_tE [(x_tx_t -H) δ ] ] =(−γt)[t−1fluct,1∣] =(I- _tH)E [ η^fluct,1_t-1 ] =. = 0. This proves (B.17) by induction. We now show that the cross term vanishes. Define, for any deterministic vector u, the linear map ℒt():=t[(−γttt⊤)⊗(tt⊤−)].L_t(u):=E_x_t [(I- _tx_tx_t )u (x_tx_t -H) δ ]. The map ℒt(⋅)L_t(·) is linear in u. Since ztz_t is Z-measurable, and since tx_t is independent of ∨σ(1,…,t−1)Z σ(x_1,…,x_t-1), we obtain t(1) _t^(1) =γt[(−γttt⊤)t−1fluct,1⊗zt(tt⊤−)] = _tE [(I- _tx_tx_t ) η^fluct,1_t-1 z_t(x_tx_t -H) δ ] =γt[[(−γt⊤)t−1fluct,1⊗zt(t⊤−)∣,1,…,t−1]] = _tE [E [(I- _tx_tx_t ) η^fluct,1_t-1 z_t(x_tx_t -H) δ ,x_1,…,x_t-1 ] ] =γt[ztℒt(t−1fluct,1)] = _tE [z_t\,L_t ( η^fluct,1_t-1 ) ] =γt[[ztℒt(t−1fluct,1)∣]] = _tE [E [z_t\,L_t ( η^fluct,1_t-1 ) ] ] =γt[ztℒt([t−1fluct,1∣])] = _tE [z_t\,L_t (E[ η^fluct,1_t-1 ] ) ] =, = 0, where the penultimate equality uses the linearity of ℒtL_t, and the last equality follows from (B.17). Hence t(1)=.C_t^(1)= 0. Taking transpose also gives (t(1))⊤=.(C_t^(1)) = 0. Therefore (B.16) reduces to t(1)=(ℐ−γtt)∘t−1(1)+γt2tfluct,1,0(1)=.F_t^(1)=(I- _tT_t) _t-1^(1)+ _t^2 _t^fluct,1, _0^(1)= 0. (B.18) Since zt∈0,1z_t∈\0,1\, we have [zt2]=[zt]=pE[z_t^2]=E[z_t]=p. Moreover, tfluct,1 _t^fluct,1 =[zt2((tt⊤−))⊗((tt⊤−))] =E [z_t^2 ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ] ⪯p[((tt⊤−))⊗((tt⊤−))]. p\,E [ ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ]. Write :=⊗A:= δ δ. Then [((tt⊤−))⊗((tt⊤−))] [ ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ] =[(tt⊤−)(tt⊤−)] =E [(x_tx_t -H)A(x_tx_t -H) ] =[tt⊤tt⊤]−. =E[x_tx_t Ax_tx_t ]-HAH. Since ⪰0HAH 0, it follows that tfluct,1 _t^fluct,1 ⪯p[tt⊤tt⊤]. p\,E[x_tx_t Ax_tx_t ]. (B.19) Now apply Assumption 3.2 with the PSD matrix =⊗A= δ δ: [tt⊤tt⊤] [x_tx_t Ax_tx_t ] ⪯αtr()=α‖2. \,α\, tr(HA)\,H=\,α\,\| δ\|_H^2\,H. (B.20) Combining (B.19) and (B.20), we obtain tfluct,1⪯pα‖2. _t^fluct,1 p\,α\| δ\|_H^2\,H. (B.21) Comparing (B.18) with the variance recursion (B.3), we see that t(1)F_t^(1) satisfies exactly the same recursion form as the variance term, with noise covariance tfluct,1 _t^fluct,1 in place of t _t. By (B.21), we may apply Lemma B.2 with σ2←pα‖2.σ^2← pα\| δ\|_H^2. Therefore, ⟨,N(1)⟩≤8pα‖21−γαtr()⋅DN, ,F_N^(1) ≤ 8pα\| δ\|_H^21-γα tr(H)· D_ effN_ eff, (B.22) where k∗=maxk:λk≥1γN,D:=k∗+γ2N2∑i>k∗λi2.k^*= \k: _k≥ 1γ N_ eff \, D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Consequently, by (B.14), ⟨,N⟩≤2⟨,N(1)⟩+2⟨,N(2)⟩, ,F_N ≤ 2 ,F_N^(1) +2 ,F_N^(2) , and the first term is controlled by (B.22). The second term ⟨,N(2)⟩ ,F_N^(2) is be bounded by Lemma B.6. ∎ Lemma B.5 (Fixed-budget variance identity). Let m∈0,1,…,Nm∈\0,1,…,N\, and set p:=mN.p:= mN. Suppose that (z1,…,zN)(z_1,…,z_N) is uniformly distributed over all binary vectors in 0,1N\0,1\^N with exactly m ones. Let 1,…,N∈ℝℋ v_1,…, v_N ^H be deterministic vectors, and define ¯:=1N∑s=1Ns. v:= 1N _s=1^N v_s. Then ‖∑s=1N(zs−p)s‖2=Np(1−p)N−1∑s=1N‖s−¯‖2.E \| _s=1^N(z_s-p) v_s \|_H^2= Np(1-p)N-1 _s=1^N\| v_s- v\|_H^2. (B.23) In particular, ‖∑s=1N(zs−p)s‖2≤Np(1−p)N−1∑s=1N‖s‖2.E \| _s=1^N(z_s-p) v_s \|_H^2≤ Np(1-p)N-1 _s=1^N\| v_s\|_H^2. (B.24) Proof. Since the law of (z1,…,zN)(z_1,…,z_N) is exchangeable and exactly m coordinates are equal to 11, we have [zs]=mN=p,s=1,…,N.E[z_s]= mN=p, s=1,…,N. Moreover, because zs∈0,1z_s∈\0,1\, [(zs−p)2]=Var(zs)=p(1−p).E[(z_s-p)^2]=Var(z_s)=p(1-p). For s≠rs≠ r, using the fixed-budget constraint ∑t=1Nzt=m _t=1^Nz_t=m, we also have [zszr]=m(m−1)N(N−1).E[z_sz_r]= m(m-1)N(N-1). Therefore [(zs−p)(zr−p)] [(z_s-p)(z_r-p)] =[zszr]−p2 =E[z_sz_r]-p^2 =m(m−1)N(N−1)−m2N2 = m(m-1)N(N-1)- m^2N^2 =−m(N−m)N2(N−1) =- m(N-m)N^2(N-1) =−p(1−p)N−1. =- p(1-p)N-1. Now expand the quadratic form: ‖∑s=1N(zs−p)s‖2 \| _s=1^N(z_s-p) v_s \|_H^2 =[∑s=1N∑r=1N(zs−p)(zr−p)⟨s,r⟩] =E [ _s=1^N _r=1^N(z_s-p)(z_r-p) v_s, v_r _H ] =∑s=1N[(zs−p)2]‖s‖2+∑s≠r[(zs−p)(zr−p)]⟨s,r⟩ = _s=1^NE[(z_s-p)^2]\| v_s\|_H^2+ _s≠ rE[(z_s-p)(z_r-p)] v_s, v_r _H =p(1−p)∑s=1N‖s‖2−p(1−p)N−1∑s≠r⟨s,r⟩. =p(1-p) _s=1^N\| v_s\|_H^2- p(1-p)N-1 _s≠ r v_s, v_r _H. Thus ‖∑s=1N(zs−p)s‖2=p(1−p)N−1[(N−1)∑s=1N‖s‖2−∑s≠r⟨s,r⟩].E \| _s=1^N(z_s-p) v_s \|_H^2= p(1-p)N-1 [(N-1) _s=1^N\| v_s\|_H^2- _s≠ r v_s, v_r _H ]. (B.25) On the other hand, ∑s=1N‖s−¯‖2 _s=1^N\| v_s- v\|_H^2 =∑s=1N‖s‖2−2∑s=1N⟨s,¯⟩+∑s=1N‖¯‖2 = _s=1^N\| v_s\|_H^2-2 _s=1^N v_s, v _H+ _s=1^N\| v\|_H^2 =∑s=1N‖s‖2−N‖¯‖2. = _s=1^N\| v_s\|_H^2-N\| v\|_H^2. Since N2‖¯‖2=‖∑s=1Ns‖2=∑s=1N‖s‖2+∑s≠r⟨s,r⟩,N^2\| v\|_H^2= \| _s=1^N v_s \|_H^2= _s=1^N\| v_s\|_H^2+ _s≠ r v_s, v_r _H, it follows that N∑s=1N‖s−¯‖2 N _s=1^N\| v_s- v\|_H^2 =N∑s=1N‖s‖2−(∑s=1N‖s‖2+∑s≠r⟨s,r⟩) =N _s=1^N\| v_s\|_H^2- ( _s=1^N\| v_s\|_H^2+ _s≠ r v_s, v_r _H ) =(N−1)∑s=1N‖s‖2−∑s≠r⟨s,r⟩. =(N-1) _s=1^N\| v_s\|_H^2- _s≠ r v_s, v_r _H. Substituting this identity into (B.25) yields ‖∑s=1N(zs−p)s‖2=Np(1−p)N−1∑s=1N‖s−¯‖2,E \| _s=1^N(z_s-p) v_s \|_H^2= Np(1-p)N-1 _s=1^N\| v_s- v\|_H^2, which proves (B.23). Finally, since ∑s=1N‖s−¯‖2=∑s=1N‖s‖2−N‖¯‖2≤∑s=1N‖s‖2, _s=1^N\| v_s- v\|_H^2= _s=1^N\| v_s\|_H^2-N\| v\|_H^2≤ _s=1^N\| v_s\|_H^2, we obtain (B.24). ∎ Lemma B.6 (An upper bound for the second fluctuation term). Recall that tfluct,2=(−γttt⊤)t−1fluct,2+γt(zt−p),0fluct,2=, η^fluct,2_t=(I- _tx_tx_t ) η^fluct,2_t-1+ _t(z_t-p)H δ, η^fluct,2_0= 0, and define t(2):=[tfluct,2⊗tfluct,2].F_t^(2):=E[ η^fluct,2_t η^fluct,2_t]. Suppose Assumption 3.2 holds and γ<1αtr().γ< 1α tr(H). Let N:=N/logN,N_ eff:=N/ N, Then ⟨,N(2)⟩≤8Np(1−p)(N−1)(1−αγtr())‖2⋅DN, ,F_N^(2) ≤ 8Np(1-p)(N-1) (1-αγ tr(H) )\| δ\|_H^2· D_ effN_ eff, where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. For each t=1,…,Nt=1,…,N, define t:=−γttt⊤.B_t:=I- _tx_tx_t . Iterating the recursion for tfluct,2 η^fluct,2_t yields Nfluct,2=∑s=1N−1⋯s+1γs(zs−p), η^fluct,2_N= _s=1^NB_NB_N-1·sB_s+1\, _s(z_s-p)H δ, where the empty product is interpreted as I. Therefore Nfluct,2=∑s=1N(zs−p)s,s:=γsNN−1⋯s+1. η^fluct,2_N= _s=1^N(z_s-p) v_s, v_s:= _sB_NB_N-1·sB_s+1H δ. Conditional on the feature sequence (1,…,N)(x_1,…,x_N), the vectors 1,…,N v_1,…, v_N are deterministic. Since the source indicators are independent of the features and follow the fixed-budget law, Lemma B.5 gives [∥Nfluct,2∥2|1,…,N]≤Np(1−p)N−1∑s=1N∥s∥2.E\! [\| η^fluct,2_N\|_H^2\, |\,x_1,…,x_N ]≤ Np(1-p)N-1 _s=1^N\| v_s\|_H^2. Taking expectation over the features gives ⟨,N(2)⟩=[‖Nfluct,2‖2]≤Np(1−p)N−1∑s=1N[‖s‖2]. ,F_N^(2) =E\! [\| η^fluct,2_N\|_H^2 ]≤ Np(1-p)N-1 _s=1^NE\! [\| v_s\|_H^2 ]. (B.26) For each t=0,1,…,Nt=0,1,…,N and each s∈1,…,Ns∈\1,…,N\, define s,t:=γstt−1⋯s+1,1≤s≤t,,s>t. v_s,t:= cases _sB_tB_t-1·sB_s+1H δ,&1≤ s≤ t,\\[2.84526pt] 0,&s>t. cases In particular, s,N=s v_s,N= v_s. Now define t:=∑s=1ts,ts,t⊤,t=0,1,…,N,S_t:= _s=1^t v_s,t v_s,t , t=0,1,…,N, with the convention 0:=S_0:= 0. Then ⟨,N⟩=∑s=1N⟨,ss⊤⟩=∑s=1N‖s‖2. ,S_N = _s=1^N , v_s v_s = _s=1^N\| v_s\|_H^2. Moreover, for s<ts<t, s,t=ts,t−1,whilet,t=γt. v_s,t=B_t v_s,t-1, v_t,t= _tH δ. Hence t=tt−1t⊤+γt2,:=⊤.S_t=B_tS_t-1B_t + _t^2O, :=H δ δ H. Define ¯t:=[t]. S_t:=E[S_t]. Since t−1S_t-1 is measurable with respect to σ(1,…,t−1)σ(x_1,…,x_t-1) and tx_t is independent of the past, taking conditional expectation yields ¯t=(ℐ−γtt)∘¯t−1+γt2,¯0=. S_t=(I- _tT_t) S_t-1+ _t^2O, S_0= 0. We first bound the forcing matrix O by a multiple of H. For any vector u, ⊤=(⊤)2≤(⊤)(⊤)=‖2⊤,u Ou=( δ Hu)^2≤( δ H δ)(u Hu)=\| δ\|_H^2\,u Hu, where we used Cauchy–Schwarz in the inner product induced by H. Therefore ⪯‖2.O \| δ\|_H^2\,H. (B.27) Now define an auxiliary matrix sequence (tδ)t=0N(C_t^δ)_t=0^N by 0δ:=,tδ=(ℐ−γt)∘t−1δ+γt2∥2,t=1,…,N.C_0^δ:= 0, _t^δ=(I- _tT_t) _t-1^δ+ _t^2\| δ\|_H^2\,H, t=1,…,N. By Lemma A.1, the operator (ℐ−γtt)(I- _tT_t) is a PSD mapping. Since ¯0=0δ= S_0=C_0^δ= 0 and (B.27) holds, an induction on t shows that ¯t⪯tδ,t=0,1,…,N. S_t _t^δ, t=0,1,…,N. Consequently, ∑s=1N[‖s‖2]=⟨,¯N⟩≤⟨,Nδ⟩. _s=1^NE\! [\| v_s\|_H^2 ]= , S_N ≤ ,C_N^δ . (B.28) The recursion defining tδC_t^δ has exactly the same form as the variance recursion (B.3), with noise covariance δ:=‖2. ^δ:=\| δ\|_H^2\,H. Applying Lemma B.2 with σ2←‖2σ^2←\| δ\|_H^2 yields ⟨,Nδ⟩≤8‖21−αγtr()⋅DN. ,C_N^δ ≤ 8\| δ\|_H^21-αγ tr(H)· D_ effN_ eff. Combining this with (B.26) and (B.28), we obtain ⟨,N(2)⟩≤Np(1−p)N−1⋅8‖21−αγtr()⋅DN, ,F_N^(2) ≤ Np(1-p)N-1· 8\| δ\|_H^21-αγ tr(H)· D_ effN_ eff, ∎ B.6 Proof of Theorem 4.1 Proof. Combining Bias upper bound, Variance upper bound, Drift upper bound and Fluctuation upper bound yields the desired result. ∎ Appendix C Mixed Training Lower Bound Analysis In this section, we derive a lower bound for the mixed-training excess risk. The starting point is again the error recursion t=(−γttt⊤)t−1+γtzttt⊤+γtξtt,0=0−1∗. η_t=(I- _tx_tx_t ) η_t-1+ _tz_tx_tx_t δ+ _t _tx_t, η_0=w_0-w_1^*. (C.1) Recall that t=t−1∗. η_t=w_t-w_1^*. By definition of excess risk, [ℛ(N)−ℛ(1∗)]=12⟨,[N⊗N]⟩.E[R(w_N)-R(w_1^*)]= 12\, ,E[ η_N η_N] . To derive a lower bound, we directly center the total error process t η_t around its mean. Define ¯t:=[t],~t:=t−¯t. η_t:=E[ η_t], η_t:= η_t- η_t. Then t=¯t+~t,[~t]=. η_t= η_t+ η_t, [ η_t]= 0. Therefore, [N⊗N] [ η_N η_N] =[(¯N+~N)⊗(¯N+~N)] =E [( η_N+ η_N) ( η_N+ η_N) ] =¯N⊗¯N+[~N⊗~N]+¯N⊗[~N]+[~N]⊗¯N = η_N η_N+E[ η_N η_N]+ η_N [ η_N]+E[ η_N] η_N =¯N⊗¯N+[~N⊗~N]. = η_N η_N+E[ η_N η_N]. Hence [ℛ(N)−ℛ(1∗)]=12‖¯N‖2+12⟨,[~N⊗~N]⟩.E[R(w_N)-R(w_1^*)]= 12\| η_N\|_H^2+ 12 ,E[ η_N η_N] . (C.2) C.1 The mean recursion Taking expectation on both sides of (C.1), and using that [tt⊤]=,[zt]=p,[ξtt]=,E[x_tx_t ]=H, [z_t]=p, [ _tx_t]= 0, we obtain ¯t=(−γt)¯t−1+γtp,¯0=0−1∗. η_t=(I- _tH) η_t-1+ _tpH δ, η_0=w_0-w_1^*. Unrolling the recursion gives ¯N η_N =∏t=1N(−γt)(0−1∗)+p∑s=1Nγs(∏t=s+1N(−γt)). = _t=1^N(I- _tH)(w_0-w_1^*)+p _s=1^N _s ( _t=s+1^N(I- _tH) )H δ. (C.3) Since each factor (−γt)(I- _tH) is a polynomial in H, it commutes with H. Therefore, ∑s=1Nγs(∏t=s+1N(−γt))=−∏t=1N(−γt). _s=1^N _s ( _t=s+1^N(I- _tH) )H=I- _t=1^N(I- _tH). Substituting this identity into (C.3), we get the closed form ¯N=∏t=1N(−γt)(0−1∗)+p(−∏t=1N(−γt)). η_N= _t=1^N(I- _tH)(w_0-w_1^*)+p (I- _t=1^N(I- _tH) ) δ. Consequently, ‖¯N‖2=‖∏t=1N(−γt)(0−1∗)+p(−∏t=1N(−γt))‖2.\| η_N\|_H^2= \| _t=1^N(I- _tH)(w_0-w_1^*)+p (I- _t=1^N(I- _tH) ) δ \|^2_H. (C.4) The remaining task is to lower bound the centered fluctuation term 12⟨,[~N⊗~N]⟩, 12 ,E[ η_N η_N] , which will be handled next. C.2 The centered fluctuation term First, we introduce the lower bound Theorem of ⟨,N⟩ ,V_N in Wu et al. [2022a], which will be used in our analysis. Lemma C.1 (A variance lower bound). Suppose Assumptions 3.1 and 3.3 hold. Let N=N/logN_ eff=N/ N. Suppose N≥10N_ eff≥ 10 and γ<1/λ1γ<1/ _1. We have ⟨,N⟩≥σ12400⋅DN, ,V_N ≥ _1^2400· D_ effN_ eff, where k∗=maxk:λk≥1γN,D:=k∗+γ2N2∑i>k∗λi2.k^*= \k: _k≥ 1γ N_ eff \, D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. See proof of Theorem C.2 in Wu et al. [2022a]. ∎ We now continue the analysis of the centered fluctuation term 12⟨,[~N⊗~N]⟩. 12 ,E[ η_N η_N] . Starting from the total recursion t=(−γttt⊤)t−1+γtzttt⊤+γtξtt,0=0−1∗, η_t=(I- _tx_tx_t ) η_t-1+ _tz_tx_tx_t δ+ _t _tx_t, η_0=w_0-w_1^*, and subtracting the mean recursion ¯t=(−γt)¯t−1+γtp,¯0=0−1∗, η_t=(I- _tH) η_t-1+ _tpH δ, η_0=w_0-w_1^*, we obtain ~t η_t =(−γttt⊤)~t−1+γtξtt+γt(zttt⊤−p)+γt(−tt⊤)¯t−1. =(I- _tx_tx_t ) η_t-1+ _t _tx_t+ _t (z_tx_tx_t δ-pH δ )+ _t(H-x_tx_t ) η_t-1. (C.5) We further decompose ~t=~noiset+~rest,t=0,1,…,N, η_t= η^noise_t+ η^res_t, t=0,1,…,N, where ~tnoise=(−γttt⊤)~t−1noise+γtξtt,~0noise=,~tres=(−γttt⊤)~t−1res+γt(zttt⊤−p)+γt(−tt⊤)¯t−1,~0res=. cases η^noise_t=(I- _tx_tx_t ) η^noise_t-1+ _t _tx_t,\\[2.15277pt] η^noise_0= 0, cases cases η^res_t=(I- _tx_tx_t ) η^res_t-1+ _t (z_tx_tx_t δ-pH δ )+ _t(H-x_tx_t ) η_t-1,\\[2.15277pt] η^res_0= 0. cases Thus ~N=~Nnoise+~Nres. η_N= η^noise_N+ η^res_N. Expanding the second moment gives [~N⊗~N] [ η_N η_N] =[~Nnoise⊗~Nnoise]+[~Nres⊗~Nres] =E[ η^noise_N η^noise_N]+E[ η^res_N η^res_N] +[~Nnoise⊗~Nres]+[~Nres⊗~Nnoise]. +E[ η^noise_N η^res_N]+E[ η^res_N η^noise_N]. (C.6) We now prove that [~tnoise⊗~tres]=,∀t≥0.E[ η^noise_t η^res_t]= 0, ∀ t≥ 0. (C.7) Let t:=−γt⊤,:=σ((s,zs):1≤s≤N)B_t:=I- _tx_tx_t , :=σ ((x_s,z_s):1≤ s≤ N ) be the sigma-field generated by the entire feature sequence and source schedule. We first observe that ~tresis -measurable for every t. η^res_t G-measurable for every t. Indeed, ~0res= η^res_0= 0, and the recursion ~tres=t~t−1res+γt(zttt⊤−p)+γt(−tt⊤)¯t−1 η^res_t=B_t η^res_t-1+ _t (z_tx_tx_t δ-pH δ )+ _t(H-x_tx_t ) η_t-1 only involves (s,zs)s≤t(x_s,z_s)_s≤ t and the deterministic vector ¯t−1 η_t-1. Hence, by induction, ~tres η^res_t is G-measurable. Next, we show that [~tnoise∣]=,∀t≥0.E [ η^noise_t ]= 0, ∀ t≥ 0. (C.8) To see this, unroll the noise recursion ~tnoise=t~t−1noise+γtξtt,~0noise=. η^noise_t=B_t η^noise_t-1+ _t _tx_t, η^noise_0= 0. This gives ~tnoise=∑s=1t(t−1⋯s+1)γsξs, η^noise_t= _s=1^t (B_tB_t-1·sB_s+1 ) _s _sx_s, where the empty product is understood as the identity matrix. Conditional on G, all matrices jB_j and vectors sx_s are deterministic. Moreover, by the zero-mean noise assumption, [ξs∣]=0,s=1,…,N.E[ _s ]=0, s=1,…,N. Therefore, [~tnoise∣]=∑s=1t(t−1⋯s+1)γs[ξs∣]=,E [ η^noise_t ]= _s=1^t (B_tB_t-1·sB_s+1 ) _sx_sE[ _s ]= 0, which proves (C.8). Now, since ~tres η^res_t is G-measurable, we have [~tnoise⊗~tres] [ η^noise_t η^res_t ] =[[~tnoise⊗~tres∣]] =E [E [ η^noise_t η^res_t ] ] =[[~tnoise∣]⊗~tres] =E [E [ η^noise_t ] η^res_t ] =. = 0. This proves (C.7). Taking transpose also gives [~tres⊗~tnoise]=.E [ η^res_t η^noise_t ]= 0. Consequently, [~N⊗~N]=[~Nnoise⊗~Nnoise]+[~Nres⊗~Nres].E[ η_N η_N]=E[ η^noise_N η^noise_N]+E[ η^res_N η^res_N]. Hence ⟨,[~N⊗~N]⟩ ,E[ η_N η_N] =⟨,[~Nnoise⊗~Nnoise]⟩+⟨,[~Nres⊗~Nres]⟩. = ,E[ η^noise_N η^noise_N] + ,E[ η^res_N η^res_N] . (C.9) We now identify the first term with a variance-type recursion. Define tmix:=[~tnoise⊗~tnoise].V_t^mix:=E[ η^noise_t η^noise_t]. Then tmixV_t^mix satisfies tmix=(ℐ−γtt)∘t−1mix+γt2¯,0mix=,V_t^mix=(I- _tT_t) _t-1^mix+ _t^2 , _0^mix= 0, where ¯:=[ξt2tt⊤]=(1−p)1+p2. :=E[ _t^2x_tx_t ]=(1-p) _1+p _2. Under Assumption 3.3, we have 1=σ12,2=σ22, _1= _1^2H, _2= _2^2H, and therefore ¯=((1−p)σ12+pσ22). = ((1-p) _1^2+p _2^2 )H. Denote σ¯2:=(1−p)σ12+pσ22. σ^2:=(1-p) _1^2+p _2^2. Then tmixV_t^mix is exactly the same variance recursion as in Lemma C.1, with σ12 _1^2 replaced by σ¯2 σ^2. Hence ⟨,[~Nnoise⊗~Nnoise]⟩=⟨,Nmix⟩≥σ¯2400⋅DN. ,E[ η^noise_N η^noise_N] = ,V_N^mix ≥ σ^2400· D_ effN_ eff. (C.10) The remaining term ⟨,[~Nres⊗~Nres]⟩ ,E[ η^res_N η^res_N] will be lower bounded separately by the following Lemma. Lemma C.2 (A lower bound for the residual fluctuation term). Recall that ~tres=(−γttt⊤)~t−1res+γt(zttt⊤−p)+γt(−tt⊤)¯t−1,~0res=, η^res_t=(I- _tx_tx_t ) η^res_t-1+ _t (z_tx_tx_t δ-pH δ )+ _t(H-x_tx_t ) η_t-1, η^res_0= 0, Suppose Assumptions 3.1, 3.2 and 3.3 hold. Let N:=N/logN,N_ eff:=N/ N, Assume N≥500N≥ 500 and γ<1/λ1γ<1/ _1. Then ⟨,[~Nres⊗~Nres]⟩≥βe−8400∥0−1∗−p∥k∗:∞2⋅DN+(e−4−2−7)240p(1−p)γ2N∑i>k∗λi3δi2, ,E[ η^res_N η^res_N] ≥ β e^-8400\,\|w_0-w_1^*-p δ\|_H_k^*:∞^2· D_ effN_ eff+ (e^-4-2^-7)^240\,p(1-p)γ^2N_ eff _i>k^* _i^3 _i^2, where k∗=maxk:λk≥1γN,D:=k∗+γ2N2∑i>k∗λi2.k^*= \k: _k≥ 1γ N_ eff \, D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. Lett−1=p−¯t−1,c_t-1=p δ- η_t-1, we can rewrite the driving term in the recursion of ~tres η^res_t as zttt⊤−p+(−tt⊤)¯t−1=(zt−p)tt⊤+(tt⊤−)t−1.z_tx_tx_t δ-pH δ+(H-x_tx_t ) η_t-1=(z_t-p)x_tx_t δ+(x_tx_t -H)c_t-1. Hence ~tres=(−γttt⊤)~t−1res+γt(zt−p)tt⊤+γt(tt⊤−)t−1,~0res=. η^res_t=(I- _tx_tx_t ) η^res_t-1+ _t(z_t-p)x_tx_t δ+ _t(x_tx_t -H)c_t-1, η^res_0= 0. (C.11) On the other hand, from the closed form of the mean process, ¯t=(−γt)¯t−1+γtp, η_t=(I- _tH) η_t-1+ _tpH δ, we obtain t=p−¯t=(−γt)(p−¯t−1)=(−γt)t−1,c_t=p δ- η_t=(I- _tH)(p δ- η_t-1)=(I- _tH)c_t-1, In particular, 0=p−(0−1∗)=−(0−1∗−p).c_0=p δ-(w_0-w_1^*)=-(w_0-w_1^*-p δ). Let :=σ(z1,…,zN),:=σ(1,…,N),Z:=σ(z_1,…,z_N), :=σ(x_1,…,x_N), and define Z:=[~Nres∣],X:=[~Nres∣],∘:=~Nres−Z−X.S_Z:=E[ η^res_N ], _X:=E[ η^res_N ], := η^res_N-S_Z-S_X. (C.12) We first show that [~Nres]=.E[ η^res_N]= 0. (C.13) Indeed, taking expectation in (C.11) and using [zt−p]=0,[tt⊤−]=,[tt⊤]=,E[z_t-p]=0, [x_tx_t -H]= 0, [x_tx_t ]=H, gives [~tres]=(−γt)[~t−1res],E[ η^res_t]=(I- _tH)E[ η^res_t-1], and since ~0res= η^res_0= 0, (C.13) follows. Next, Z and X are independent, and [X]=[~Nres]=0E[S_X]=E[ η^res_N]=0, hence [X∣]=[X]=0.E[S_X ]=E[S_X]=0. Therefore, ⟨Z,X⟩=[⟨Z,[X∣]⟩]=0.E _Z,HS_X =E [ _Z,H\,E[S_X ] ]=0. Moreover, by (C.12), [∘∣]=[~Nres∣]−Z−[X∣]=0,E[R ]=E[ η^res_N ]-S_Z-E[S_X ]=0, and similarly [∘∣]=0E[R ]=0. Thus ⟨Z,∘⟩=0,⟨X,∘⟩=0.E _Z,HR =0, _X,HR =0. Consequently, ‖~Nres‖2=‖Z‖2+‖X‖2+‖∘‖2≥‖Z‖2+‖X‖2.E\| η^res_N\|_H^2=E\|S_Z\|_H^2+E\|S_X\|_H^2+E\|R \|_H^2 \|S_Z\|_H^2+E\|S_X\|_H^2. (C.14) Since ⟨,[~Nres⊗~Nres]⟩=‖~Nres‖2 ,E[ η^res_N η^res_N] =E\| η^res_N\|_H^2, it remains to lower bound the two terms on the right-hand side of (C.14). Define t:=[~tres∣].f_t:=E[ η^res_t ].Taking conditional expectation in (C.11) and using [zt−p∣]=0E[z_t-p ]=0, we obtain t=(−γttt⊤)t−1+γt(tt⊤−)t−1,0=.f_t=(I- _tx_tx_t )f_t-1+ _t(x_tx_t -H)c_t-1, _0= 0. (C.15) Thus X=NS_X=f_N. Further define u:=(uu⊤−)u−1,u,N:=γu∏j=u+1N(−γjjj⊤). ζ_u:=(x_ux_u -H)c_u-1, _u,N:= _u _j=u+1^N(I- _jx_jx_j ). Unrolling (C.15) gives X=∑u=1Nu,u:=u,Nu.S_X= _u=1^NU_u, _u:=K_u,N ζ_u. (C.16) We first prove the cross terms vanish. Let u<vu<v. Since u,NK_u,N and vU_v are measurable with respect to σ(u+1,…,N)σ(x_u+1,…,x_N), while u ζ_u depends only on ux_u, we have [u]=([uu⊤]−)u−1=0.E[ ζ_u]=(E[x_ux_u ]-H)c_u-1=0. Hence [⟨u,v⟩∣u+1,…,N] [ _u,HU_v _u+1,…,x_N ] =[⟨u,u,N⊤v⟩∣u+1,…,N]=0, =E [ ζ_u,K_u,N HU_v _u+1,…,x_N ]=0, and therefore ⟨u,v⟩=0,u≠v.E _u,HU_v =0, u≠ v. It follows that ‖X‖2=∑u=1N‖u‖2.E\|S_X\|_H^2= _u=1^NE\|U_u\|_H^2. (C.17) Now let u:=u−1u−1⊤⪰0.A_u:=c_u-1c_u-1 0. Using the independence of u,NK_u,N and u ζ_u, we have ‖u‖2 \|U_u\|_H^2 =[u⊤u,N⊤u,Nu]=⟨[u,N⊤u,N],u⟩, =E [ ζ_u K_u,N HK_u,N ζ_u ]= [K_u,N HK_u,N],\,G_u , (C.18) where u:=[uu⊤]=[(uu⊤−)u(uu⊤−)].G_u:=E[ ζ_u ζ_u ]=E[(x_ux_u -H)A_u(x_ux_u -H)]. By expansion, u=[uu⊤uuu⊤]−u.G_u=E[x_ux_u A_ux_ux_u ]-HA_uH. Applying Assumption 3.2 with A=uA=A_u, we obtain u⪰βtr(u)=β‖u−1‖2.G_u β\, tr(HA_u)\,H=β\,\|c_u-1\|_H^2\,H. Substituting this into (C.18) yields ‖u‖2≥β‖u−1‖2[tr(u,Nu,N⊤)].E\|U_u\|_H^2≥β\,\|c_u-1\|_H^2\,E\! [ tr(HK_u,NHK_u,N ) ]. (C.19) Set u:=1/2u,N1/2.M_u:=H^1/2K_u,NH^1/2. Then tr(u,Nu,N⊤)=‖u‖F2. tr(HK_u,NHK_u,N )=\|M_u\|_F^2. By Jensen’s inequality, ‖u‖F2≥‖[u]‖F2.E\|M_u\|_F^2≥\|E[M_u]\|_F^2. Moreover, [u,N]=γu∏j=u+1N(−γj).E[K_u,N]= _u _j=u+1^N(I- _jH). Hence, in the eigenbasis of H, ‖[u]‖F2=∑i≥1γu2λi2∏j=u+1N(1−γjλi)2.\|E[M_u]\|_F^2= _i≥ 1 _u^2 _i^2 _j=u+1^N(1- _j _i)^2. Define κu:=∑i≥1γu2λi2∏j=u+1N(1−γjλi)2. _u:= _i≥ 1 _u^2 _i^2 _j=u+1^N(1- _j _i)^2. (C.20) Then (C.19) gives ‖u‖2≥β‖u−1‖2κu.E\|U_u\|_H^2≥β\,\|c_u-1\|_H^2\, _u. Combining this with (C.17), ‖X‖2≥β∑u=1N‖u−1‖2κu.E\|S_X\|_H^2≥β _u=1^N\|c_u-1\|_H^2\, _u. (C.21) We now lower bound ‖u−1‖2\|c_u-1\|_H^2 uniformly on the tail space. In the eigenbasis of H, (u−1)i=−∏j=1u−1(1−γjλi)(0−1∗−p)i.(c_u-1)_i=- _j=1^u-1(1- _j _i)\,(w_0-w_1^*-p δ)_i. Hence ∥u−1∥k∗:∞2=∑i>k∗λi∏j=1u−1(1−γjλi)2(0−1∗−p)i2.\|c_u-1\|_H_k^*:∞^2= _i>k^* _i _j=1^u-1(1- _j _i)^2(w_0-w_1^*-p δ)_i^2. (C.22) Since i>k∗i>k^*, we have γλi<1/N≤1/10γ _i<1/N_ eff≤ 1/10 for N≥10N_ eff≥ 10. Also ∑t=1Nγt≤2γN. _t=1^N _t≤ 2γ N_ eff. Using 1−a≥e−2a1-a≥ e^-2a for 0≤a≤1/20≤ a≤ 1/2, we obtain ∏j=1u−1(1−γjλi)2≥exp(−4λi∑j=1u−1γj)≥exp(−8γNλi)≥e−8. _j=1^u-1(1- _j _i)^2≥ \! (-4 _i _j=1^u-1 _j )≥ (-8γ N_ eff _i)≥ e^-8. Substituting this into (C.22) gives ∥u−1∥2≥∥u−1∥k∗:∞2≥e−8∥0−1∗−p∥k∗:∞2,u=1,…,N.\|c_u-1\|_H^2≥\|c_u-1\|_H_k^*:∞^2≥ e^-8\|w_0-w_1^*-p δ\|_H_k^*:∞^2, u=1,…,N. (C.23) Therefore (C.21) yields ∥X∥2≥βe−8∥0−1∗−p∥k∗:∞2∑u=1Nκu.E\|S_X\|_H^2≥β e^-8\|w_0-w_1^*-p δ\|_H_k^*:∞^2 _u=1^N _u. (C.24) Finally, by applying Lemma C.1, ∑u=1Nκu=⟨,∑u=1Nγu2∏j=u+1N(−γj)2⟩≥1400(k∗N+γ2N∑i>k∗λi2). _u=1^N _u= ,\; _u=1^N _u^2 _j=u+1^N(I- _jH)^2H ≥ 1400 ( k^*N_ eff+γ^2N_ eff _i>k^* _i^2 ). Substituting this into (C.24), we conclude that ∥X∥2≥βe−8400∥0−1∗−p∥k∗:∞2(k∗N+γ2N∑i>k∗λi2).E\|S_X\|_H^2≥ β e^-8400\,\|w_0-w_1^*-p δ\|_H_k^*:∞^2 ( k^*N_ eff+γ^2N_ eff _i>k^* _i^2 ). (C.25) Then we lower bound schedule component ZS_Z. For each t, define t:=[~tres∣].a_t:=E[ η^res_t ]. Taking conditional expectation in (C.11) and using independence of tx_t from (,1,…,t−1)(Z,x_1,…,x_t-1), we obtain t=(−γt)t−1+γt(zt−p),0=0.a_t=(I- _tH)a_t-1+ _t(z_t-p)H δ, _0=0. Thus Z=NS_Z=a_N, and unrolling gives Z=∑u=1N(zu−p)u,u:=γu(∏j=u+1N(−γj)).S_Z= _u=1^N(z_u-p)v_u, _u:= _u ( _j=u+1^N(I- _jH) )H δ. (C.26) In the eigenbasis of H, (u)i=γuλi∏j=u+1N(1−γjλi)δi.(v_u)_i= _u _i _j=u+1^N(1- _j _i) _i. By the fixed-budget variance identity, ‖Z‖2=p(1−p)2(N−1)∑u=1N∑r=1N‖u−r‖2.E\|S_Z\|_H^2= p(1-p)2(N-1) _u=1^N _r=1^N\|v_u-v_r\|_H^2. (C.27) Let I−:=t:ℓt=0,I+:=t:ℓt≥7.I_-:=\t: _t=0\, I_+:=\t: _t≥ 7\. For N≥500N≥ 500, one has |I−|≥N2,|I+|≥N10.|I_-|≥ N_ eff2, |I_+|≥ N10. Restricting (C.27) to pairs (u,r)∈I−×I+(u,r)∈ I_-× I_+ yields ‖Z‖2≥p(1−p)2(N−1)∑u∈I−∑r∈I+‖u−r‖2.E\|S_Z\|_H^2≥ p(1-p)2(N-1) _u∈ I_- _r∈ I_+\|v_u-v_r\|_H^2. (C.28) Fix i>k∗i>k^*. Since γλi<1/N≤1/10<1/2γ _i<1/N_ eff≤ 1/10<1/2, for u∈I−u∈ I_- we have γu=γ _u=γ, and |(u)i|=γλi∏j=u+1N(1−γjλi)|δi|≥e−4γλi|δi|.|(v_u)_i|=γ _i _j=u+1^N(1- _j _i)| _i|≥ e^-4γ _i| _i|. For r∈I+r∈ I_+, we have γr≤2−7γ _r≤ 2^-7γ, so |(r)i|≤2−7γλi|δi|.|(v_r)_i|≤ 2^-7γ _i| _i|. Therefore |(u−r)i|≥(e−4−2−7)γλi|δi|,u∈I−,r∈I+,i>k∗.|(v_u-v_r)_i|≥(e^-4-2^-7)γ _i| _i|, u∈ I_-, r∈ I_+, i>k^*. Hence ‖u−r‖2≥(e−4−2−7)2γ2∑i>k∗λi3δi2.\|v_u-v_r\|_H^2≥(e^-4-2^-7)^2γ^2 _i>k^* _i^3 _i^2. (C.29) Substituting (C.29) into (C.28), and using |I−|≥N/2|I_-|≥ N_ eff/2, |I+|≥N/10|I_+|≥ N/10, we obtain ‖Z‖2≥(e−4−2−7)240p(1−p)γ2N∑i>k∗λi3δi2.E\|S_Z\|_H^2≥ (e^-4-2^-7)^240\,p(1-p)γ^2N_ eff _i>k^* _i^3 _i^2. That is, ‖Z‖2≥(e−4−2−7)240p(1−p)γ2N∑i>k∗λi3δi2.E\|S_Z\|_H^2≥ (e^-4-2^-7)^240\,p(1-p)γ^2N_ eff _i>k^* _i^3 _i^2. (C.30) Finally, combining (C.14), (C.25), and (C.30), we obtain ⟨,[~Nres⊗~Nres]⟩≥βe−8400∥0−1∗−p∥k∗:∞2⋅DN+(e−4−2−7)240p(1−p)γ2N∑i>k∗λi3δi2, ,E[ η^res_N η^res_N] ≥ β e^-8400\,\|w_0-w_1^*-p δ\|_H_k^*:∞^2· D_ effN_ eff+ (e^-4-2^-7)^240\,p(1-p)γ^2N_ eff _i>k^* _i^3 _i^2, This completes the proof. ∎ C.3 Proof of Theorem 4.2 Proof. Combining Lemma C.2 with Eq. C.10 yields the desired result. ∎ C.4 Proof of Corollary 4.3 Proof. Let T:=M+N,T:=∏t=1T(−γt).T:=M+N, _T:= _t=1^T(I- _tH). Since the synthetic proportion p=M/(M+N)p=M/(M+N) is fixed, it suffices to study the limit as T→∞T→∞. We first show that TP_T vanishes in H-norm. Writing the eigendecomposition =∑iλiii⊤,H= _i _iv_iv_i , we have Ti=(∏t=1T(1−γtλi))i.P_Tv_i= ( _t=1^T(1- _t _i) )v_i. Under the assumption on the stepsize, for every i with λi>0 _i>0, 0≤1−γtλi<1.0≤ 1- _t _i<1. Moreover, under the geometric tail-decay schedule, ∑t=1Tγt≍TlogT∑ℓ=0⌊logT⌋2−ℓ≍TlogT→T→∞. _t=1^T _t T T _ =0 T 2^- T T [T→∞]∞. Hence for every λi>0 _i>0, ∏t=1T(1−γtλi)⟶0. _t=1^T(1- _t _i) 0. Therefore, for any w with ‖<∞\|w\|_H<∞, ‖T‖2=∑iλi(∏t=1T(1−γtλi))2⟨,i⟩2⟶0\|P_Tw\|_H^2= _i _i ( _t=1^T(1- _t _i) )^2 ,v_i ^2 0 by dominated convergence. We now turn to the upper bound in Theorem 4.1 which gives [ℰ1(T)] [E_1(w_T)] ≲‖T(0−1∗)‖2⏟(I)+α∥0−1∗∥20:k∗γN~+k∗:∞D~N~⏟(II) \|P_T(w_0-w_1^*)\|_H^2_(I)+ α\|w_0-w_1^*\|^2_ I_0:k^*γ N_ eff+H_k^*:∞ D_ eff N_ eff_(I) +σ¯2D~N~⏟(III)+αp‖2D~N~⏟(IV)+p2‖(−T)‖2⏟(V). + σ^2 D_ eff N_ eff_(I)+ α p\| δ\|_H^2 D_ eff N_ eff_(IV)+ p^2\|(I-P_T) δ\|_H^2_(V). By the assumption D~=o(N~) D_ eff=o( N_ eff), terms (II)(I), (III)(I) and (IV)(IV) all vanish as T→∞T→∞. By the argument above, (I)→0(I)→ 0. Finally, ‖(−T)‖→‖,\|(I-P_T) δ\|_H→\| δ\|_H, because ‖T‖→0\|P_T δ\|_H→ 0. Therefore, lim supT→∞[ℰ1(T)]≲p2‖2. _T→∞E[E_1(w_T)] p^2\| δ\|_H^2. Next, apply the lower bound in Theorem 4.2: [ℰ1(T)] [E_1(w_T)] ≳‖T(0−1∗)+p(−T)‖2⏟(A) \|P_T(w_0-w_1^*)+p(I-P_T) δ \|_H^2_(A) +(β∥0−1∗−p∥k∗:∞2+σ¯2)D~N~⏟(B)+p(1−p)γ2N~∥3k∗:∞2⏟(C). + (β\|w_0-w_1^*-p δ\|_H_k^*:∞^2+ σ^2 ) D_ eff N_ eff_(B)+ p(1-p)γ^2 N_ eff\| δ\|_H^3_k^*:∞^2_(C). Since (B)≥0(B)≥ 0 and (C)≥0(C)≥ 0, it is enough to keep only (A)(A). Rewrite T(0−1∗)+p(−T)=T(0−1∗−p)+p.P_T(w_0-w_1^*)+p(I-P_T) δ=P_T(w_0-w_1^*-p δ)+p δ. Again using ‖T‖→0\|P_Tw\|_H→ 0 for every w with finite H-norm, we obtain ‖T(0−1∗−p)+p‖2⟶p2‖2. \|P_T(w_0-w_1^*-p δ)+p δ \|_H^2 p^2\| δ\|_H^2. Hence, lim infT→∞[ℰ1(T)]≳p2‖2. _T→∞E[E_1(w_T)] p^2\| δ\|_H^2. Combining the upper and lower bounds yields limM+N→∞[ℰ1(M+N)]≂p2‖2. _M+N→∞E[E_1(w_M+N)] p^2\| δ\|_H^2. This proves the corollary. ∎ Appendix D Two-stage Training Analysis First, consider the upper bound for training on pure one-source data from [Wu et al., 2022a], which can also be recovered by Theorem 4.1. Theorem D.1 (Learning with only real data). Consider last iterate SGD with stepsize scheme (3.1). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Let N:=N/logN_ eff:=N/ N. Suppose γ<1/(4αtr())γ<1/(4α tr(H)). Then we have [ℛ(N)]−ℛ(1∗) [R(w_N)]-R(w_1^*) ≲∥(∏t=1N−γt)(0∗−1∗)∥2+(α∥0−1∗∥0:k∗γN+k∗:∞2+σ12)DN \|( _t=1^NI- _tH)(w_0^*-w_1^*) \|^2_H+ (α\|w_0-w_1^*\|^2_ I_0:k^*γ N_ eff+H_k^*:∞+ _1^2 ) D_ effN_ eff with k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. We now turn to the two-stage training procedure. The key idea is to condition on the first-stage training iterate and then apply the pure real-data bound from Theorem 4.1 to the real-data training phase. Suppose that in the first stage we run SGD on M synthetic samples and obtain Mw_M. In the second stage, starting from Mw_M, we run SGD on N real samples and output M+Nw_M+N. Our goal is to upper bound [ℛ(M+N)−ℛ(1∗)].E [R(w_M+N)-R(w_1^*) ]. Let ℱM:=σ(1,y1,…,M,yM)F_M:=σ(x_1,y_1,…,x_M,y_M) be the sigma-field generated by the first stage. Conditional on ℱMF_M, the iterate Mw_M is deterministic, and the second-stage procedure is exactly SGD trained on N real samples, initialized at Mw_M. Therefore, applying Theorem 4.1 conditionally, we obtain [ℛ(M+N)−ℛ(1∗)∣ℱM] [R(w_M+N)-R(w_1^*) _M ] ≲‖(∏t=1N(−γt))(M−1∗)‖2 \| ( _t=1^N(I- _tH) )(w_M-w_1^*) \|_H^2 +(α∥M−1∗∥0:k∗γN+k∗:∞2+σ12)DN, + (α\|w_M-w_1^*\|^2_ I_0:k^*γ N_ eff+H_k^*:∞+ _1^2 ) D_ effN_ eff, (D.1) Taking expectation with respect to ℱMF_M on both sides of (D.1), we get [ℛ(M+N)−ℛ(1∗)] [R(w_M+N)-R(w_1^*) ] =[[ℛ(M+N)−ℛ(1∗)∣ℱM]] =E [E [R(w_M+N)-R(w_1^*) _M ] ] ≲[‖(∏t=1N(−γt))(M−1∗)‖2] [ \| ( _t=1^N(I- _tH) )(w_M-w_1^*) \|_H^2 ] +(α[∥M−1∗∥0:k∗γN+k∗:∞2]+σ12)DN. + ( [\|w_M-w_1^*\|^2_ I_0:k^*γ N_ eff+H_k^*:∞ ]+ _1^2 ) D_ effN_ eff. (D.2) We now relate the terms to the stage-one excess risk under the synthetic distribution. Recall that M−1∗=(M−2∗)+.w_M-w_1^*=(w_M-w_2^*)+ δ. Hence [‖(∏t=1N(−γt))(M−1∗)‖2] [ \| ( _t=1^N(I- _tH) )(w_M-w_1^*) \|_H^2 ] =[‖(∏t=1N(−γt))(M−2∗)+(∏t=1N(−γt))‖2] =E [ \| ( _t=1^N(I- _tH) )(w_M-w_2^*)+ ( _t=1^N(I- _tH) ) δ \|_H^2 ] ≤2[‖(∏t=1N(−γt))(M−2∗)‖2]+2‖(∏t=1N(−γt))‖2. ≤ 2E [ \| ( _t=1^N(I- _tH) )(w_M-w_2^*) \|_H^2 ]+2 \| ( _t=1^N(I- _tH) ) δ \|_H^2. (D.3) Now denote N:=∏t=1N(−γt).P_N:= _t=1^N(I- _tH). Since NP_N is a polynomial in H, it commutes with H. Therefore ‖N‖2=⟨,NN⟩=⟨,1/2N21/2⟩≤‖2,\|P_Nv\|_H^2= ,P_NHP_Nv = ,H^1/2P_N^2H^1/2v ≤\|v\|_H^2, because 0⪯N⪯0 _N under the stepsize assumption γt<1/λ1 _t<1/ _1. Hence [‖N(M−2∗)‖2] [\|P_N(w_M-w_2^*)\|_H^2 ] ≤[‖M−2∗‖2] [\|w_M-w_2^*\|_H^2 ] =2[L2(M)−L2(2∗)]. =2\,E[L_2(w_M)-L_2(w_2^*)]. (D.4) Next, let :=0:k∗γN+k∗:∞.B:= I_0:k^*γ N_ eff+H_k^*:∞. By the definition of k∗k^*, λi≥1γN,i≤k∗. _i≥ 1γ N_ eff, i≤ k^*. Therefore 0:k∗γN⪯0:k∗, I_0:k^*γ N_ eff _0:k^*, and hence =0:k∗γN+k∗:∞⪯0:k∗+k∗:∞=.B= I_0:k^*γ N_ eff+H_k^*:∞ _0:k^*+H_k^*:∞=H. (D.5) It follows that [‖M−2∗‖2] [\|w_M-w_2^*\|_B^2 ] ≤[‖M−2∗‖2] [\|w_M-w_2^*\|_H^2 ] =2[L2(M)−L2(2∗)]. =2\,E[L_2(w_M)-L_2(w_2^*)]. (D.6) On the other hand, [‖M−1∗‖2] [\|w_M-w_1^*\|_B^2 ] =[‖(M−2∗)+‖2] =E [\|(w_M-w_2^*)+ δ\|_B^2 ] ≤2[‖M−2∗‖2]+2‖2 ≤ 2E [\|w_M-w_2^*\|_B^2 ]+2\| δ\|_B^2 ≤4[L2(M)−L2(2∗)]+2‖2, ≤ 4\,E[L_2(w_M)-L_2(w_2^*)]+2\| δ\|_B^2, (D.7) where in the last step we used (D.6). Substituting (D.3), (D.4), and (D.7) into (D.2), we obtain [ℛ(M+N)−ℛ(1∗)] [R(w_M+N)-R(w_1^*) ] ≲[ℛ2(M)−ℛ2(2∗)]+‖N‖2 [R_2(w_M)-R_2(w_2^*)]+ \|P_N δ \|_H^2 +(α[ℛ2(M)−ℛ2(2∗)]+α‖2+σ12)DN. + (α\,E[R_2(w_M)-R_2(w_2^*)]+α\| δ\|_B^2+ _1^2 ) D_ effN_ eff. (D.8) Finally, by (D.5), ‖2≤‖2.\| δ\|_B^2≤\| δ\|_H^2. Therefore [ℛ(M+N)−ℛ(1∗)] [R(w_M+N)-R(w_1^*) ] ≲[ℛ2(M)−ℛ2(2∗)]+‖N‖2 [R_2(w_M)-R_2(w_2^*)]+ \|P_N δ \|_H^2 +(α[ℛ2(M)−ℛ2(2∗)]+α‖2+σ12)DN. + (α\,E[R_2(w_M)-R_2(w_2^*)]+α\| δ\|_H^2+ _1^2 ) D_ effN_ eff. (D.9) Equivalently, writing 2(M):=[ℛ2(M)−ℛ2(2∗)], Excess_2(w_M):=E[R_2(w_M)-R_2(w_2^*)], we have (M+N) Excess(w_M+N) ≲2(M)+‖N‖2+(α2(M)+α‖2+σ12)DN. Excess_2(w_M)+ \|P_N δ \|_H^2+ (α\, Excess_2(w_M)+α\| δ\|_H^2+ _1^2 ) D_ effN_ eff. (D.10) This is the desired two-stage upper bound. D.1 Two-stage training without learning-rate restart The two-stage analysis restarts the geometrically decaying stepsize schedule at the beginning of the real-data stage. We briefly discuss how the argument changes when a single uninterrupted schedule is used across both stages. Let T:=M+NT:=M+N, and denote the global T-step schedule by γt(T):=γ2⌊t/(T/logT)⌋,t=1,…,T. _t^(T):= γ2 t/(T/ T) , t=1,…,T. Under the restarted protocol analyzed above, the real-data stage instead uses a fresh N-step schedule γs(N):=γ2⌊s/(N/logN)⌋,s=1,…,N. _s^(N):= γ2 s/(N/ N) , s=1,…,N. Thus, without restarting, the contraction operator appearing in the two-stage bound changes from Nrst:=∏s=1N(−γs(N))P_N^rst:= _s=1^N(I- _s^(N)H) to the real-data tail of the global schedule M,Ntail:=∏t=M+1M+N(−γt(T)).P_M,N^tail:= _t=M+1^M+N(I- _t^(T)H). The proof follows the same conditional argument as above. Conditional on the first-stage iterate Mw_M, the second stage remains single-source SGD on real data. The decomposition M−1∗=(M−2∗)+w_M-w_1^*=(w_M-w_2^*)+ δ is unchanged, as are the subsequent comparisons between the stage-one excess risk and the corresponding H-weighted errors. The only modification is that the proof of the single-source bound must be applied to the shifted stepsize sequence γM+1(T),…,γM+N(T)\ _M+1^(T),…, _M+N^(T)\ rather than to a freshly restarted sequence γ1(N),…,γN(N)\ _1^(N),…, _N^(N)\. Consequently, every occurrence of NrstP_N^rst is replaced by M,NtailP_M,N^tail, and the spectral cutoff and effective-dimension factor are defined using the real-stage tail schedule. In particular, the resulting bound has the same qualitative decomposition: [ℰ1(M+N)]≲ [E_1(w_M+N)] [ℰ2(M)]+‖M,Ntail‖2 [E_2(w_M)]+ \|P_M,N^tail δ \|_H^2 +(α[ℰ2(M)]+α‖2+σ12)DefftailTeff, + ( [E_2(w_M)]+α\| δ\|_H^2+ _1^2 ) D_eff^tailT_eff, where Teff:=T/logT_eff:=T/ T, ktail∗:=maxk:λk≥1/(γM+1(T)Teff)k_tail^*:= \k: _k≥ 1/( _M+1^(T)T_eff)\, and Defftail:=ktail∗+(γM+1(T)Teff)2∑i>ktail∗λi2D_eff^tail:=k_tail^*+( _M+1^(T)T_eff)^2 _i>k_tail^* _i^2. Nevertheless, the qualitative asymptotic conclusion is unchanged. For every eigendirection with λi>0 _i>0, ∏t=M+1M+N(1−γt(T)λi) _t=M+1^M+N (1- _t^(T) _i ) ≤exp(−λi∑t=M+1M+Nγt(T)). ≤ (- _i _t=M+1^M+N _t^(T) ). When p=M/(M+N)∈(0,1)p=M/(M+N)∈(0,1) is fixed, the real-data tail of the geometric schedule satisfies ∑t=M+1M+Nγt(T)→∞. _t=M+1^M+N _t^(T)→∞. It follows that ‖M,Ntail‖2→0. \|P_M,N^tail δ \|_H^2→ 0. Moreover, under the corresponding effective-dimension consistency condition, Defftail/T→0.D_eff^tail/T→ 0. Therefore, using an uninterrupted schedule slows the finite-sample forgetting of the synthetic initialization, but it does not introduce a persistent synthetic drift during the real-data stage and does not create a non-vanishing excess-risk floor. Appendix E Random Sketch Model Following [Lin et al., 2024], we begin by decomposing the population risk into three components: irreducible risk, approximation error, and excess risk: ℛD(M+N)=minℛ(⋅)⏟+minℛD(⋅)−minℛ(⋅)⏟+ℛD(M+N)−minℛD(⋅)⏟.R_D(v_M+N)= (·)_ Irreducible+ _D(·)- (·)_ Approx+ R_D(v_M+N)- _D(·)_ Excess. (E.1) For the Irreducible risk, under the well-specified model Assumption 4.6, =ℛ(1∗)=12σ12. Irreducible=R(w_1^*)= 12 _1^2. For Approximation Error, as established in Lemma C.4 in [Lin et al., 2024], under Assumption 4.6, with probability at least 1−eΩ(D)1-e (D), 1∗≂D1−a.E_w_1^* Approx D^1-a. So in the following subsections, we focus on derive the bound for E Excess using risk upper bound under Assumption 4.6 and 4.7. Unless other illustration, expectations are conditioned on S. E.1 Mixed Training Upper Bound Recall the upper bound in Theorem 4.1, where the data covariance becomes ⊤SHS and the optimal parameters becomes 1∗,2∗v_1^*,v_2^*. Suppose Assumptions 3.2, 4.6 hold. Let N~:=(M+N)/log(M+N) N_ eff:=(M+N)/ (M+N), σ¯2=(1−p)σ12+pσ22 σ^2=(1-p) _1^2+p _2^2 and λ~j λ_j be the eigenvalue of ⊤SHS . Suppose γ<1/(4αtr(⊤))γ<1/(4α tr(SHS )). Then we have Excess ≲Bias+Var+αp‖2∗−1∗‖⊤2D~N~⏟+p2‖2∗−1∗‖⊤2⏟, +Var+ α\,p\|v_2^*-v_1^*\|^2_SHS D_ eff N_ eff_ Fluct+ p^2 \|v_2^*-v_1^* \|^2_SHS _ Drift, with ℓ∗:=(⊤)−1ℓ∗,andD:=#λ~j≥1/(N~γ)+(N~γ)2∑λ~j<1/(N~γ)λ~j2. ^*_ :=(SHS )^-1SHw^*_ , D_ eff:=\#\ λ_j≥ 1/( N_ effγ)\+( N_ effγ)^2 _ λ_j<1/( N_ effγ) λ_j^2. Note that when σ12,σ22≂1 _1^2, _2^2 1, Theorem A.4 in Lin et al. [2024] shows that Bias+Var≲‖∏t=1M+N(−γt⊤)(0−1∗)‖⊤2+DN~Bias+Var \| _t=1^M+N(I- _tSHS )(v_0-v_1^*) \|^2_SHS + D_ eff N_ eff Future assume that Assumption 4.7 holds and zero initialiazation where 0=v_0= 0, Lin et al. [2024](Appendix D and E) shows that with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S 1∗‖∏t=1M+N(−γt⊤)(1∗)‖⊤2≲maxD1−a,(N~γ)1/a−1, _w_1^* \| _t=1^M+N(I- _tSHS )(v_1^*) \|^2_SHS \D^1-a,\ ( N_ effγ)^1/a-1 \, 1∗‖∏t=1M+N(−γt⊤)(1∗)‖⊤2≳(N~γ)1/a−1 when (N~γ)1/a≤D/c for some constant c>0, _w_1^* \| _t=1^M+N(I- _tSHS )(v_1^*) \|^2_SHS ( N_ effγ)^1/a-1 when ( N_ effγ)^1/a≤ D/c for some constant c>0, and DN~≂minD,(N~γ)1/a/N~(N~γ)1/a−1. D_ eff N_ eff \D,\ ( N_ effγ)^1/a \/ N_ eff( N_ effγ)^1/a-1. For the Fluct term, we have Fluct≲‖2∗−1∗‖⊤2⋅DN~Fluct \|v_2^*-v_1^*\|_SHS ^2· D_ eff N_ eff We verify that ‖2∗−1∗‖⊤2 \|v_2^*-v_1^*\|_SHS ^2 =‖12⊤(2∗−1∗)‖2 =\|H 12S (v_2^*-v_1^*)\|^2 =‖12⊤(⊤)−1‖2 =\|H 12S (SHS )^-1SH δ\|^2 ≤‖12‖2=‖2, ≤\|H 12 δ\|^2=\| δ\|^2_H, which implies that ‖2∗−1∗‖⊤2≲‖2≂∑i−b≂1E_ δ\|v_2^*-v_1^*\|_SHS ^2 \| δ\|^2_H _i\,i^-b 1 (E.2) So that with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S Fluct≲DN~≂minD,(N~γ)1/a/N~≲D1−a+(N~γ)1/a−1. _ δFluct D_ eff N_ eff \D,\ ( N_ effγ)^1/a \/ N_ eff D^1-a+( N_ effγ)^1/a-1. For the Drift term p2‖2∗−1∗‖⊤2p^2\|v_2^*-v_1^*\|_SHS ^2, ‖2∗−1∗‖⊤2 \|v_2^*-v_1^*\|_SHS ^2 =‖(⊤)−1(2∗−1∗)‖⊤2 =\|(SHS )^-1SH(w_2^*-w_1^*)\|_SHS ^2 =⊤⊤(⊤)−1(⊤)(⊤)−1 = δ HS (SHS )^-1(SHS )(SHS )^-1SH δ =⊤⊤(⊤)−1+⊤−⊤ = δ HS (SHS )^-1SH δ+ δ H δ- δ δ =⊤−⊤(−⊤(⊤)−1) = δ H δ- δ (H-HS (SHS )^-1SH) δ =‖2−‖(−12⊤(⊤)−1⊤12)12‖2 =\| δ\|^2_H- \| (I-H 12S (SHS )^-1S H 12 )H 12 δ \|^2 For the second minus term, Lin et al. [2024](Lemma C.5) shows that for any random vector x, let (λi,i)i≥1( _i,v_i)_i≥ 1 be the eigenvalue-eigenvector pairs of H, if the Assumption 4.7 (power-law decay) satisfies and [⟨i,⟩⟨j,⟩]=0E[ _i,x _j,x ]=0 for i≠ji≠ j and [λi⟨i,⟩2]≂i−bE[ _i _i,x ^2] i^-b for some b>1b>1, Then with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S (,,):=‖(−12⊤(⊤)−1⊤12)12‖2≂D1−bE_x Approx(S,H,x):= \| (I-H 12S (SHS )^-1S H 12 )H 12x \|^2 D^1-b So that under Assumption 4.7, ‖(−12⊤(⊤)−1⊤12)12‖2≂D1−bE_ δ \| (I-H 12S (SHS )^-1S H 12 )H 12 δ \|^2 D^1-b (E.3) Taking expectation in the decomposition and using (E.3), we obtain ‖2∗−1∗‖⊤2 _ δ\|v_2^*-v_1^*\|_SHS ^2 =‖2−Θ(D1−b). =E_ δ\| δ\|_H^2- (D^1-b). (E.4) Since ‖2≂1E_ δ\| δ\|_H^2 1, we suppress this problem-dependent constant in the scaling-law statement and write, for sufficiently large D, ‖2∗−1∗‖⊤2≂1−D1−b.E_ δ\|v_2^*-v_1^*\|_SHS ^2 1-D^1-b. (E.5) Therefore, _ δ Drift =p2‖2∗−1∗‖⊤2 =p^2E_ δ\|v_2^*-v_1^*\|_SHS ^2 ≂p2(1−D1−b). p^2(1-D^1-b). Finally, put , Irreducible, Approx and Excess together and choose γ≂1γ 1, we have that with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S: ℛD(M+N) _D(v_M+N) ≲σ12+1Da−1+1(N~)1−1/a+(M+N)2(1−1Db−1) _1^2+ 1D^a-1+ 1( N_ eff)^1-1/a+ ( MM+N )^2(1- 1D^b-1) E.2 Mixed Training Lower Bound We now establish the matching lower bound in the optimization-saturated regime γN~≳Da.γ N_ eff D^a. Applying Theorem 4.2 to the sketched problem, with H replaced by ⊤SHS and using 0=v_0=0, gives ≳1∗,‖−(∏t=1M+N(−γt⊤))1∗+p(−∏t=1M+N(−γt⊤))(2∗−1∗)‖⊤2, Excess _w_1^*, δ \|- ( _t=1^M+N(I- _tSHS ) )v_1^*+p (I- _t=1^M+N(I- _tSHS ) )(v_2^*-v_1^*) \|_SHS ^2, where we have dropped the remaining nonnegative terms in Theorem 4.2. By Assumption 4.6, 1∗,[1∗(2∗−1∗)⊤] _w_1^*, δ [v_1^*(v_2^*-v_1^*) ] =(⊤)−1[1∗⊤]⊤(⊤)−1 =(SHS )^-1SHE[w_1^* δ ]HS (SHS )^-1 =. =0. Therefore, the cross term vanishes after taking expectation over (1∗,)(w_1^*, δ), and hence Excess ≳p2‖(−∏t=1M+N(−γt⊤))(2∗−1∗)‖⊤2. p^2E_ δ \| (I- _t=1^M+N(I- _tSHS ) )(v_2^*-v_1^*) \|_SHS ^2. (E.6) It remains to control the spectral filter in (E.6). By Lemma 6.2 of Lin et al. [2024], with probability at least 1−e−Ω(D)1-e^- (D), the eigenvalues of ⊤SHS satisfy λ~j≂j−a,j=1,…,D. λ_j j^-a, j=1,…,D. Moreover, the geometric stepsize schedule satisfies ∑t=1M+Nγt≂γN~. _t=1^M+N _t γ N_ eff. Thus, for every j=1,…,Dj=1,…,D, ∏t=1M+N(1−γtλ~j) _t=1^M+N(1- _t λ_j) ≤exp(−λ~j∑t=1M+Nγt) ≤ (- λ_j _t=1^M+N _t ) ≤exp(−cγN~D−a)≤e−c, ≤ (-cγ N_ effD^-a )≤ e^-c, where the last inequality follows from γN~≳Daγ N_ eff D^a. Let (λ~j,~j)j=1D( λ_j, v_j)_j=1^D be the eigenvalue-eigenvector pairs of ⊤SHS . Expanding the norm in this eigenbasis gives ‖(−∏t=1M+N(−γt⊤))(2∗−1∗)‖⊤2 \| (I- _t=1^M+N(I- _tSHS ) )(v_2^*-v_1^*) \|_SHS ^2 =∑j=1Dλ~j(1−∏t=1M+N(1−γtλ~j))2⟨~j,2∗−1∗⟩2 = _j=1^D λ_j (1- _t=1^M+N(1- _t λ_j) )^2 v_j,v_2^*-v_1^* ^2 ≥(1−e−c)2∑j=1Dλ~j⟨~j,2∗−1∗⟩2 ≥(1-e^-c)^2 _j=1^D λ_j v_j,v_2^*-v_1^* ^2 =(1−e−c)2‖2∗−1∗‖⊤2. =(1-e^-c)^2\|v_2^*-v_1^*\|_SHS ^2. Therefore, ‖(−∏t=1M+N(−γt⊤))(2∗−1∗)‖⊤2≳‖2∗−1∗‖⊤2. _ δ \| (I- _t=1^M+N(I- _tSHS ) )(v_2^*-v_1^*) \|_SHS ^2 _ δ\|v_2^*-v_1^*\|_SHS ^2. The upper-bound analysis above has already established that ‖2∗−1∗‖⊤2≂1−D1−b.E_ δ\|v_2^*-v_1^*\|_SHS ^2 1-D^1-b. Substituting this result into (E.6) gives ≳p2(1−D1−b).E Excess p^2(1-D^1-b). In addition, the approximation-error result stated above gives 1∗≂D1−a.E_w_1^* Approx D^1-a. Combining this with (E.1), we obtain ℛD(M+N)−ℛ(1∗)≳D1−a+p2(1−D1−b).ER_D(v_M+N)-R(w_1^*) D^1-a+p^2(1-D^1-b). Finally, since γ≂1γ 1 and γN~≳Daγ N_ eff D^a, we have N~1/a−1≲D1−a. N_ eff^1/a-1 D^1-a. Therefore, ℛD(M+N)−ℛ(1∗)≳D1−a+N~1/a−1+p2(1−D1−b),ER_D(v_M+N)-R(w_1^*) D^1-a+ N_ eff^1/a-1+p^2(1-D^1-b), which matches the upper bound in Theorem 4.8 up to constants. E.3 Two-stage Training Analysis First, Recall the upper bound in Theorem 4.4. For a tighter analysis, we need a intermediate result from its proof. The results are for specific 1∗,2∗w_1^*,w_2^* and the expectation is taken with respect to the one-stage M-step SGD. Suppose Assumptions 3.1, 3.2 and 4.6 hold. Let N:=N/logN_ eff:=N/ N and M=M/logM_ eff=M/ M. Suppose γ<1/(4αtr(⊤))γ<1/(4α tr(SHS )). Then we have (M+N) Excess(v_M+N) ≲[‖(∏t=1N(−γt⊤))(M−2∗)‖⊤2]+‖(∏t=1N−γt⊤)(2∗−1∗)‖⊤2 [ \| ( _t=1^N(I- _tSHS ) )(v_M-v_2^*) \|_SHS ^2 ]+ \|( _t=1^NI- _tSHS )(v_2^*-v_1^*) \|^2_SHS +α(2(M)+‖2∗−1∗‖⊤2)DN+ +α ( Excess_2(v_M)+\|v_2^*-v_1^*\|^2_SHS ) D_ effN_ eff+ Var where ℓ∗:=(⊤)−1ℓ∗,andD:=#λ~j≥1/(Nγ)+(Nγ)2∑λ~j<1/(Nγ)λ~j2. ^*_ :=(SHS )^-1SHw^*_ , D_ eff:=\#\ λ_j≥ 1/(N_ effγ)\+(N_ effγ)^2 _ λ_j<1/(N_ effγ) λ_j^2. Suppose σ12,σ22≂1 _1^2, _2^2 1. For the Var term, similar to the analysis in mixed training analysis and [Lin et al., 2024](Lemma E.1), with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S ≂DN≂minD,(Nγ)1/a/N. Var D_ effN_ eff \D,\ (N_ effγ)^1/a \/N_ eff. Analysis of 2(M) Excess_2(v_M). We now analyze the first-stage excess risk when the teacher is 2∗w_2^*. Conditioning on the sketch matrix S, define 2∗:=(⊤)−12∗,M:=M/logM.v_2^*:=(SHS )^-1SHw_2^*, M_ eff:=M/ M. Using Theorem 6.1 in Lin et al. [2024], for γ<1/(4αtr(⊤))γ<1/(4α\, tr(SHS )) and zero initialization, 2(M)≲2+D1M,E Excess_2(v_M) Bias_2+ D_ eff_1M_ eff, where D1:=#λ~j≥1/(Mγ)+(Mγ)2∑λ~j<1/(Mγ)λ~j2,D_ eff_1:=\#\ λ_j≥ 1/(M_ effγ)\+(M_ effγ)^2 _ λ_j<1/(M_ effγ) λ_j^2, and the expectation is over the randomness of 2∗w_2^* and the first-stage SGD. By Assumption 4.7, let (λi,i)i≥1( _i,v_i)_i≥ 1 be the eigenvalue-eigenvector pairs of H, and shares the eigenvectors (i)i≥1(v_i)_i≥ 1 with H, i.e., i=τii,λi≂i−a,λiτi≂i−b, v_i= _iv_i, _i i^-a, _i _i i^-b, for some a>1a>1 and b>1b>1. Then ⟨i,2∗⟩2=⟨i,1∗⟩2+⟨i,⟩2=1+τi,E _i,w_2^* ^2=E _i,w_1^* ^2+E _i, δ ^2=1+ _i, where we used [1∗⊤]=0E[w_1^* δ ]=0. For the bias term, applying Lemma D.1 in Lin et al. [2024] with a cutoff k≤D/3k≤ D/3 yields 2≲∑i≤k⟨i,2∗⟩2Mγ+∑i>kλi⟨i,2∗⟩2.E Bias_2 _i≤ kE _i,w_2^* ^2M_ effγ+ _i>k _iE _i,w_2^* ^2. Substituting ⟨i,2∗⟩2=1+τiE _i,w_2^* ^2=1+ _i gives 2≲∑i≤k(1+τi)Mγ+∑i>kλi(1+τi).E Bias_2 _i≤ k(1+ _i)M_ effγ+ _i>k _i(1+ _i). Using τi≂ia−b _i i^a-b, we obtain ∑i≤k(1+τi)≂k+k1+a−b,∑i>kλi(1+τi)≂k1−a+k1−b, _i≤ k(1+ _i) k+k^1+a-b, _i>k _i(1+ _i) k^1-a+k^1-b, and hence 2≲k+k1+a−bMγ+k1−a+k1−b.E Bias_2 k+k^1+a-bM_ effγ+k^1-a+k^1-b. Let c:=mina,b.c:= \a,b\. and distinguish two cases. If b≥ab≥ a, then k1+a−b≲k^1+a-b k and k1−b≲k1−ak^1-b k^1-a, so 2≲kMγ+k1−a.E Bias_2 kM_ effγ+k^1-a. If b<ab<a, then k≲k1+a−bk k^1+a-b and k1−a≲k1−bk^1-a k^1-b, so 2≲k1+a−bMγ+k1−b.E Bias_2 k^1+a-bM_ effγ+k^1-b. Hence, in both cases, 2≲k1+a−cMγ+k1−c,c=mina,b.E Bias_2 k^1+a-cM_ effγ+k^1-c, c= \a,b\. Choosing k≂minD,(Mγ)1/a,k \D,(M_ effγ)^1/a\, yields 2≲maxD1−c,(Mγ)(1−c)/a.E Bias_2 \D^1-c,(M_ effγ)^(1-c)/a\. For the lower bound, note that [2∗(2∗)⊤]=+,E[w_2^*(w_2^*) ]=I+ , under the assumption [1∗⊤]=0E[w_1^* δ ]=0. Since shares the eigenvectors of H and λiτi≂i−b _i _i i^-b, the eigenvalues of (+)H(I+ ) satisfy λi(1+τi)≂i−a+i−b≂i−c,c=mina,b. _i(1+ _i) i^-a+i^-b i^-c, c= \a,b\. Therefore, applying the same argument as in Lemma D.2 and Lemma D.4 of Lin et al. [2024], we obtain 2≳(Mγ)(1−c)/awhen(Mγ)1/a≤D/c0E Bias_2 (M_ effγ)^(1-c)/a (M_ effγ)^1/a≤ D/c_0 for some constant c0>0c_0>0. Consequently, with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S, the upper bound is 2≲maxD1−c,(Mγ)(1−c)/a,c=mina,b,E Bias_2 \D^1-c,(M_ effγ)^(1-c)/a\, c= \a,b\, Moreover, in the data-limited regime (Mγ)1/a≤D/c0(M_ effγ)^1/a≤ D/c_0, this upper bound is tight up to constants. Finally, by Lemma E.1 in Lin et al. [2024], D1M≂minD,(Mγ)1/aM D_ eff_1M_ eff \D,(M_ effγ)^1/a\M_ eff with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S. Combining the above estimates and choose γ≂1γ 1, we conclude that with probability at least 1−e−Ω(D)1-e^- (D) over S, 2(M) Excess_2(v_M) ≲maxD1−a∧b,(Mγ)(1−a∧b)/a+minD,(Mγ)1/aM≲1 \D^1-a b,(M_ effγ)^(1-a b)/a\+ \D,(M_ effγ)^1/a\M_ eff 1 For the term α(2(M)+‖2∗−1∗‖⊤2)DN,α ( Excess_2(v_M)+\|v_2^*-v_1^*\|^2_SHS ) D_ effN_ eff, we claim that it is of the same order as D/ND_ eff/N_ eff. Indeed, by the analysis in mixed-training-analysis, with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S, ‖2∗−1∗‖⊤2≲1. _ δ\|v_2^*-v_1^*\|^2_SHS 1. Therefore, on the intersection of the above high-probability events, α(2(M)+‖2∗−1∗‖⊤2)DN α (E Excess_2(v_M)+E_ δ\|v_2^*-v_1^*\|^2_SHS ) D_ effN_ eff ≲DN≂minD,(Nγ)1/aN. D_ effN_ eff \D,(N_ effγ)^1/a\N_ eff. For the term ‖(∏t=1N−γt⊤)(2∗−1∗)‖⊤2 \|( _t=1^NI- _tSHS )(v_2^*-v_1^*) \|^2_SHS , under Assumption 4.7, we can apply Lemma D.4 in Lin et al. [2024] directly and get that probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S: ‖(∏t=1N−γt⊤)(2∗−1∗)‖⊤2≲maxD1−b,(Nγ)(1−b)/a, _ δ \|( _t=1^NI- _tSHS )(v_2^*-v_1^*) \|^2_SHS \D^1-b,(N_ effγ)^(1-b)/a \, and ‖(∏t=1N−γt⊤)(2∗−1∗)‖⊤2≳(Nγ)(1−b)/a _ δ \|( _t=1^NI- _tSHS )(v_2^*-v_1^*) \|^2_SHS (N_ effγ)^(1-b)/a when (Nγ)1/a≤D/c(N_ effγ)^1/a≤ D/c for some constant c>0c>0. Moreover, when b≥a+1b≥ a+1 ‖(∏t=1N−γt⊤)(2∗−1∗)‖⊤2≲logN⋅max(Nγ)(1−b)/a,D1−b _ δ \|( _t=1^NI- _tSHS )(v_2^*-v_1^*) \|^2_SHS N_ eff· \(N_ effγ)^(1-b)/a,\ D^1-b \ For the last term [‖(∏t=1N(−γt⊤))(M−2∗)‖⊤2],E [ \| ( _t=1^N(I- _tSHS ) )(v_M-v_2^*) \|_SHS ^2 ], let N:=∏t=1N(−γt⊤).P_N:= _t=1^N(I- _tSHS ). Conditioning on the sketch matrix S, we diagonalize ⊤=diag(λ~1,…,λ~D)⊤SHS =U diag( λ_1,…, λ_D)U with λ~1≥⋯≥λ~D>0 λ_1≥·s≥ λ_D>0. Then N⊤⊤N=diag(λ~i∏t=1N(1−γtλ~i)2)⊤⪯(∏t=1N(1−γtλ~D))2⊤.P_N SHS P_N=U diag\! ( λ_i _t=1^N(1- _t λ_i)^2 )U ( _t=1^N(1- _t λ_D) )^2SHS . Hence [‖N(M−2∗)‖⊤2|]≤(∏t=1N(1−γtλ~D))2[‖M−2∗‖⊤2|].E\! [\|P_N(v_M-v_2^*)\|_SHS ^2\, |\,S ]≤ ( _t=1^N(1- _t λ_D) )^2E\! [\|v_M-v_2^*\|_SHS ^2\, |\,S ]. Using 1−u≤e−u1-u≤ e^-u and ∑t=1Nγt≂Nγ _t=1^N _t N_ effγ, we further obtain [‖N(M−2∗)‖⊤2|]≲exp(−cNγμD(⊤))[‖M−2∗‖⊤2|]E\! [\|P_N(v_M-v_2^*)\|_SHS ^2\, |\,S ] \! (-c\,N_ effγ\, _D(SHS ) )\,E\! [\|v_M-v_2^*\|_SHS ^2\, |\,S ] for some constant c>0c>0. Since [‖M−2∗‖⊤2|]=[2(M)|],E\! [\|v_M-v_2^*\|_SHS ^2\, |\,S ]=E\! [ Excess_2(v_M)\, |\,S ], it follows that [‖N(M−2∗)‖⊤2|]≲exp(−cNγμD(⊤))[2(M)|].E\! [\|P_N(v_M-v_2^*)\|_SHS ^2\, |\,S ] \! (-c\,N_ effγ\, _D(SHS ) )\,E\! [ Excess_2(v_M)\, |\,S ]. Now let ℰspecE_spec be the event that μD(⊤)≂D−a, _D(SHS ) D^-a, and let ℰexE_ex be the event that [2(M)|]≲maxD1−a∧b,(Mγ)(1−a∧b)/a.E\! [ Excess_2(v_M)\, |\,S ] \D^1-a b,(M_ effγ)^(1-a b)/a\. By Lemma 6.2 in Lin et al. [2024] and the above analysis of 2(M) Excess_2(v_M), both events hold with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S. Therefore, on the intersection event ℰspec∩ℰexE_spec _ex, we have [‖(∏t=1N(−γt⊤))(M−2∗)‖⊤2]≲exp(−cNγDa)maxD1−a∧b,(Mγ)(1−a∧b)/a.E [ \| ( _t=1^N(I- _tSHS ) )(v_M-v_2^*) \|_SHS ^2 ] \! (-c\, N_ effγD^a ) \D^1-a b,(M_ effγ)^(1-a b)/a\. Consequently, with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S, [‖(∏t=1N(−γt⊤))(M−2∗)‖⊤2]≲exp(−cNγDa)maxD1−a∧b,(Mγ)(1−a∧b)/a.E [ \| ( _t=1^N(I- _tSHS ) )(v_M-v_2^*) \|_SHS ^2 ] \! (-c\, N_ effγD^a ) \D^1-a b,(M_ effγ)^(1-a b)/a\. Finally, put , Irreducible, Approx and Excess together. Choose γ≂1γ 1 and assume b<a+1b<a+1 for simplicity and clarity, we have that with probability at least 1−e−Ω(D)1-e^- (D) over the randomness of the sketch matrix S: ℛD(M+N) _D(v_M+N) ≲σ12+1Da−1+minD,(N)1/aN⏟ _1^2+ 1D^a-1+ \D,(N_ eff)^1/a\N_ eff_ Var +e−Ω(NDa)⋅max1Da∧b−1,1(M)a∧b−1a+max1Db−1,1(N)b−1a⏟ + e^- ( N_ effD^a)· \ 1D^a b-1, 1(M_ eff) a b-1a \+ \ 1D^b-1, 1(N_ eff) b-1a \_ Bias Appendix F Exact Synthetic-Benefit Condition Proof of Theorem 4.10. Reindex the N real-stage samples as (t,yt)t=1N\(x_t,y_t)\_t=1^N, where yt=⟨t,1∗⟩+ξ1,t.y_t= _t,w_1^* + _1,t. For any real-stage iterate, the error recursion is t−1∗=(−γttt⊤)(t−1−1∗)+γtξ1,tt.w_t-w_1^*=(I- _tx_tx_t )(w_t-1-w_1^*)+ _t _1,tx_t. Iterating this recursion for real-only training gives Nreal−1∗ _N^real-w_1^* =ΦN(0−1∗) = _N(w_0-w_1^*) +∑t=1Nγt[(−γN⊤)⋯(−γt+1t+1t+1⊤)]ξ1,t, + _t=1^N _t [(I- _Nx_Nx_N )·s(I- _t+1x_t+1x_t+1 ) ] _1,tx_t, (F.1) where the product inside the summation is interpreted as I when t=Nt=N. The two-stage procedure uses exactly the same real samples, label noises, and stepsizes, but starts the real stage from Mw_M. Therefore, M+N−1∗ _M+N-w_1^* =ΦN(M−1∗) = _N(w_M-w_1^*) +∑t=1Nγt[(−γN⊤)⋯(−γt+1t+1t+1⊤)]ξ1,t. + _t=1^N _t [(I- _Nx_Nx_N )·s(I- _t+1x_t+1x_t+1 ) ] _1,tx_t. (F.2) In particular, the real-stage noise accumulation in (F.1) and (F.2) is identical. Let N ν_N denote the common real-stage noise-accumulation term in (F.1) and (F.2). These two recursions can be written as Nreal−1∗=ΦN(0−1∗)+N,M+N−1∗=ΦN(M−1∗)+N.w_N^real-w_1^*= _N(w_0-w_1^*)+ ν_N, _M+N-w_1^*= _N(w_M-w_1^*)+ ν_N. Conditional on (M,1,…,N)(w_M,x_1,…,x_N), the coefficients of all ξ1,t _1,t in N ν_N are fixed. Since the real label noises are independent of the covariates and the synthetic stage, and satisfy [ξ1,t]=0E[ _1,t]=0, we have [N∣M,1,…,N]=.E[ ν_N _M,x_1,…,x_N]=0. Both ΦN(0−1∗) _N(w_0-w_1^*) and ΦN(M−1∗) _N(w_M-w_1^*) are measurable under this conditioning. Therefore, by the tower property, [(ΦN(0−1∗))⊤N] \! [ ( _N(w_0-w_1^*) ) H ν_N ] =0, =0, [(ΦN(M−1∗))⊤N] \! [ ( _N(w_M-w_1^*) ) H ν_N ] =0. =0. Expanding the two squared errors consequently gives ‖Nreal−1∗‖2 \|w_N^real-w_1^*\|_H^2 =‖ΦN(0−1∗)‖2+‖N‖2, =E\| _N(w_0-w_1^*)\|_H^2+E\| ν_N\|_H^2, ‖M+N−1∗‖2 \|w_M+N-w_1^*\|_H^2 =‖ΦN(M−1∗)‖2+‖N‖2. =E\| _N(w_M-w_1^*)\|_H^2+E\| ν_N\|_H^2. The final terms are identical because the two procedures use the same real samples, label noises, and stepsizes. They therefore cancel in the risk difference. Under Assumption 3.3, for any w, ℰ1()=12‖−1∗‖2.E_1(w)= 12\|w-w_1^*\|_H^2. Hence, 2([ℰ1(Nreal)]−[ℰ1(M+N)]) 2 (E[E_1(w_N^real)]-E[E_1(w_M+N)] ) =‖ΦN(0−1∗)‖2−‖ΦN(M−1∗)‖2. =E\| _N(w_0-w_1^*)\|_H^2-E\| _N(w_M-w_1^*)\|_H^2. Since 0w_0 is deterministic, the definition of NK_N gives ‖ΦN(0−1∗)‖2 \| _N(w_0-w_1^*)\|_H^2 =(0−1∗)⊤1,…,N[ΦN⊤ΦN](0−1∗) =(w_0-w_1^*) E_x_1,…,x_N[ _N H _N](w_0-w_1^*) =‖0−1∗‖N2. =\|w_0-w_1^*\|_K_N^2. Similarly, because Mw_M is independent of the real-stage covariates, the tower property yields ‖ΦN(M−1∗)‖2 \| _N(w_M-w_1^*)\|_H^2 =M[(M−1∗)⊤1,…,N[ΦN⊤ΦN](M−1∗)] =E_w_M\! [(w_M-w_1^*) E_x_1,…,x_N[ _N H _N](w_M-w_1^*) ] =‖M−1∗‖N2. =E\|w_M-w_1^*\|_K_N^2. Finally, decompose M−1∗=(M−[M])+([M]−1∗)w_M-w_1^*= (w_M-E[w_M] )+ (E[w_M]-w_1^* ) and Expandthe corresponding quadratic form to get ‖M−1∗‖N2 \|w_M-w_1^*\|_K_N^2 =‖M−[M]‖N2+‖[M]−1∗‖N2 =E\|w_M-E[w_M]\|_K_N^2+\|E[w_M]-w_1^*\|_K_N^2 +2[(M−[M])⊤N([M]−1∗)]. +2E\! [ (w_M-E[w_M] ) K_N (E[w_M]-w_1^* ) ]. The last term is zero because [M−[M]]=E[w_M-E[w_M]]=0. Moreover, ‖M−[M]‖N2 \|w_M-E[w_M]\|_K_N^2 =tr(N(M−[M])(M−[M])⊤) =E tr\! (K_N(w_M-E[w_M])(w_M-E[w_M]) ) =tr(NCov(M)). = tr\! (K_NCov(w_M) ). Combining the preceding identities yields 2([ℰ1(Nreal)]−[ℰ1(M+N)]) 2 (E[E_1(w_N^real)]-E[E_1(w_M+N)] ) =‖0−1∗‖N2−‖[M]−1∗‖N2−tr(NCov(M)). =\|w_0-w_1^*\|_K_N^2-\|E[w_M]-w_1^*\|_K_N^2- tr\! (K_NCov(w_M) ). Therefore, [ℰ1(M+N)]<[ℰ1(Nreal)]E[E_1(w_M+N)]<E[E_1(w_N^real)] if and only if ‖[M]−1∗‖N2+tr(NCov(M))<‖0−1∗‖N2,\|E[w_M]-w_1^*\|_K_N^2+ tr\! (K_NCov(w_M) )<\|w_0-w_1^*\|_K_N^2, which proves Theorem 4.10. ∎ Proof of Corollary 4.11. Since xs2=1x_s^2=1, the s-th synthetic update with stepsize 1/s1/s satisfies ws w_s =ws−1−1s(ws−1xs−(θ+δ)xs−ξ2,s)xs =w_s-1- 1s (w_s-1x_s-(θ+δ)x_s- _2,s )x_s =(1−1s)ws−1+1s(θ+δ+ξ2,sxs). = (1- 1s )w_s-1+ 1s (θ+δ+ _2,sx_s ). Induction over s therefore gives wM=θ+δ+1M∑s=1Mξ2,sxs.w_M=θ+δ+ 1M _s=1^M _2,sx_s. Because the synthetic noises are independent, centered, and independent of the Rademacher covariates, [wM]=θ+δ,Var(wM)=σ22M.E[w_M]=θ+δ, (w_M)= _2^2M. (F.3) In one dimension, =[x2]=1H=E[x^2]=1. Moreover, every real-stage update has multiplicative operator 1−γtxt2=1−γt1- _tx_t^2=1- _t. Thus, ΦN=∏t=1N(1−γt),N=(∏t=1N(1−γt))2. _N= _t=1^N(1- _t), _N= ( _t=1^N(1- _t) )^2. Applying Theorem 4.10 and using (F.3), we obtain 2([ℰ1(wNreal)]−[ℰ1(wM+N)]) 2 (E[E_1(w_N^real)]-E[E_1(w_M+N)] ) =(∏t=1N(1−γt))2[(w0−θ)2−δ2−σ22M]. = ( _t=1^N(1- _t) )^2 [(w_0-θ)^2-δ^2- _2^2M ]. Since 0<γt<10< _t<1, the multiplicative factor is strictly positive. Therefore, two-stage training has strictly smaller expected risk if and only if δ2+σ22M<(w0−θ)2.δ^2+ _2^2M<(w_0-θ)^2. If δ2≥(w0−θ)2δ^2≥(w_0-θ)^2, this inequality cannot hold for any M. Otherwise, it is equivalent to M>σ22(w0−θ)2−δ2,M> _2^2(w_0-θ)^2-δ^2, which completes the proof. ∎ Appendix G constant stepsize SGD with iterate averaging In this section, we also denote total sample size as N, where synthetic sample size is pNpN and real sample size is (1−p)N(1-p)N. This light abuse of notation will achieve more simplicity and clarity. G.1 Upper bound analysis We use the same definition in Appendix B. t=[tbias⊗tbias],t=[tvar⊗tvar], _t=E[ η^bias_t η^bias_t], _t=E[ η^var_t η^var_t], t=[tdrift⊗tdrift],t=[tfluct⊗tfluct]. _t=E[ η^drift_t η^drift_t], _t=E[ η^fluct_t η^fluct_t]. Using the definition of t,t,t,tB_t,V_t,D_t,F_t, we can then decompose Excess Risk using Cauchy–Schwarz Inequality: [ℛ(¯N)−ℛ(1∗)] [R( w_N)-R(w_1^*)] =12⟨,[¯N⊗¯N]⟩ = 12 ,E[ η_N η_N] =12⟨,[(¯Nbias+¯Nvar+¯Ndrift+¯Nfluct)⊗(¯Nbias+¯Nvar+¯Ndrift+¯Nfluct)]⟩ = 12 ,E[( η_N^bias+ η_N^var+ η_N^drift+ η_N^fluct) ( η_N^bias+ η_N^var+ η_N^drift+ η_N^fluct)] ≤2⟨,[¯Nbias⊗¯Nbias]⟩+2⟨,[¯Nvar⊗¯Nvar]⟩ ≤ 2 ,E[ η_N^bias η_N^bias] +2 ,E[ η_N^var η_N^var] +2⟨,[¯Ndrift⊗¯Ndrift]⟩+2⟨,[¯Nfluct⊗¯Nfluct]⟩ +2 ,E[ η_N^drift η_N^drift] +2 ,E[ η_N^fluct η_N^fluct] =2⟨,¯N⟩+2⟨,¯N⟩+2⟨,¯N⟩+2⟨,¯N⟩. =2 , B_N +2 , V_N +2 , D_N +2 , F_N . Bias upper bound Using the defined operators, the update rule of the iterates imply the following recursive form of tB_t: t=(ℐ−γ)∘t−1,0=0⊗0,B_t=(I-γT) _t-1, _0= η_0 η_0, (G.1) Note that this bias term ⟨,¯N⟩ , B_N is the same as that in Zou et al. [2023], so we can apply Lemma B.11 in Zou et al. [2023] to bound it directly: Lemma G.1 (A bias upper bound for Avg-SGD). Suppose Assumptions 3.1 and 3.2 hold. Let N=N_ eff=N. Consider (G.1). Suppose γ<1/(αtr())γ<1/(α tr(H)). We have 12⟨,¯N⟩ 12 , B_N ≤1γ2N2⋅∥0−∗∥0:k∗−12+∥0−∗∥k∗:∞2 ≤ 1γ^2N_ eff^2·\|w_0-w^*\|_H_0:k^*^-1^2+\|w_0-w^*\|_H_k^*:∞^2 +2α(∥0−∗∥0:k∗2+Nγ∥0−∗∥k∗:∞2)Nγ(1−γαtr())⋅DN + 2α (\|w_0-w^*\|_I_0:k^*^2+N_ effγ\|w_0-w^*\|_H_k^*:∞^2 )N_ effγ(1-γα tr(H))· D_ effN_ eff where k∗=maxk:λk≥1γNk^*= \k: _k≥ 1γ N_ eff\ and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. See proof of Lemma B.11 in Zou et al. [2023]. ∎ Variance upper bound Using the defined operators, the update rule of the variance process implies t=(ℐ−γ)∘t−1+γ2t,0=,V_t=(I-γT) _t-1+γ^2 _t, _0= 0, (G.2) where t:=[ξt2tt⊤]. _t:=E[ _t^2x_tx_t ]. Since tx_t is independent of ztz_t, conditioning on the source indicator gives t=[[ξt2tt⊤∣zt]]=(1−p)1+p2⪯((1−p)σ12+pσ22). _t=E\! [E[ _t^2x_tx_t z_t] ]=(1-p) _1+p _2 ((1-p) _1^2+p _2^2 )H. Define σ¯2:=(1−p)σ12+pσ22. σ^2:=(1-p) _1^2+p _2^2. Then (G.2) has exactly the same form as the variance recursion for constant-stepsize SGD with iterate averaging analyzed in Zou et al. [2023], with noise level upper bounded by σmix2 _mix^2. Therefore, by applying the averaged-SGD variance bound in Zou et al. [2023], we obtain the following result. The corresponding sharp variance characterization for iterate averaging is stated in Theorem 2.1 and developed in Lemma B.6 of Zou et al. [2023]. Lemma G.2 (A variance upper bound for Avg-SGD). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Consider (G.2). Let N=N_ eff=N and σ¯2=(1−p)σ12+pσ22 σ^2=(1-p) _1^2+p _2^2. Suppose γ<1/(αtr())γ<1/(α tr(H)). Then 12⟨,¯N⟩≤σ¯21−γαtr()⋅DN, 12 , V_N ≤ σ^21-γα tr(H)· D_ effN_ eff, where k∗=maxk:λk≥1γN,k^*= \k: _k≥ 1γ N_ eff \, and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. The recursion (G.2) is identical to the variance recursion of constant-stepsize SGD with iterate averaging, except that the noise covariance is now t _t instead of a single-source covariance matrix. Since t⪯σ¯2for all t, _t σ^2H all t, the proof of the variance upper bound in Zou et al. [2023] applies verbatim with σ2σ^2 replaced by σmix2 _mix^2. Hence 12⟨,¯N⟩≤σ¯21−γαtr()(k∗N+Nγ2∑i>k∗λi2)=σ¯21−γαtr()⋅DN. 12 , V_N ≤ σ^21-γα tr(H) ( k^*N_ eff+N_ effγ^2 _i>k^* _i^2 )= σ^21-γα tr(H)· D_ effN_ eff. This completes the proof. ∎ Drift upper bound We now analyze the drift contribution 12⟨,¯N⟩=12⟨,[¯Ndrift⊗¯Ndrift]⟩,¯Ndrift:=1N∑t=0N−1tdrift. 12 , D_N = 12 ,E[ η^drift_N η^drift_N] , η^drift_N:= 1N _t=0^N-1 η^drift_t. Lemma G.3 (A drift upper bound for Avg-SGD). Suppose Assumptions 3.1, 3.2 and 3.3 hold. Let N=N_ eff=N. Suppose γ<1/(αtr())γ<1/(α tr(H)). Then 12⟨,¯N⟩ 12 , D_N ≤p22‖(−1γN−1(−(−γ)N))‖2 ≤ p^22 \| (I- 1γ N_ effH^-1 (I-(I- )^N_ eff ) ) δ \|_H^2 +αp2‖21−γαtr()⋅DN, + α p^2\| δ\|_H^21-γα tr(H)· D_ effN_ eff, where k∗=maxk:λk≥1γN,D:=k∗+γ2N2∑i>k∗λi2.k^*= \k: _k≥ 1γ N_ eff \, D_ eff:=k^*+γ^2N_ eff^2 _i>k^* _i^2. Proof. Recall that the drift process satisfies tdrift=(−γtt⊤)t−1drift+γp,0drift=. η^drift_t=(I- _tx_t ) η^drift_t-1+γ pH δ, η^drift_0= 0. We first separate the drift process into its mean part and centered fluctuation part. Define t:=[tdrift].m_t:=E[ η^drift_t]. Since tx_t is independent of t−1drift η^drift_t-1 and [tt⊤]=E[x_tx_t ]=H, taking expectation on both sides yields t=(−γ)t−1+γp,0=.m_t=(I- )m_t-1+γ pH δ, _0= 0. (G.3) Let ~tdrift:=tdrift−t. η^drift_t:= η^drift_t-m_t. Subtracting (G.3) from the recursion of tdrift η^drift_t gives ~tdrift=(−γtt⊤)~t−1drift+γ(−tt⊤)t−1,~0drift=. η^drift_t=(I- _tx_t ) η^drift_t-1+γ(H-x_tx_t )m_t-1, η^drift_0= 0. (G.4) Moreover, by construction, [~tdrift]=,∀t≥0.E[ η^drift_t]= 0, ∀ t≥ 0. Now define the averaged quantities ¯N:=1N∑t=0N−1t,~¯Ndrift:=1N∑t=0N−1~tdrift. m_N:= 1N _t=0^N-1m_t, η^drift_N:= 1N _t=0^N-1 η^drift_t. Then ¯Ndrift=¯N+~¯Ndrift. η^drift_N= m_N+ η^drift_N. Since [~¯Ndrift]=E[ η^drift_N]= 0, the cross term vanishes and ¯N=¯N⊗¯N+~¯N,~¯N:=[~¯Ndrift⊗~¯Ndrift]. D_N= m_N m_N+ D_N, D_N:=E\! [ η^drift_N η^drift_N ]. (G.5) Therefore, 12⟨,¯N⟩=12‖¯N‖2+12⟨,~¯N⟩. 12 , D_N = 12\| m_N\|_H^2+ 12 , D_N . (G.6) Unrolling (G.3), we obtain t=p(−(−γ)t),t≥0.m_t=p (I-(I- )^t ) δ, t≥ 0. Hence ¯N m_N =pN∑t=0N−1(−(−γ)t) = pN _t=0^N-1 (I-(I- )^t ) δ =p(−1N∑t=0N−1(−γ)t). =p (I- 1N _t=0^N-1(I- )^t ) δ. Using the matrix geometric-series identity ∑t=0N−1(−γ)t=(γ)−1(−(−γ)N), _t=0^N-1(I- )^t=( )^-1 (I-(I- )^N ), we get ¯N=p(−1γN−1(−(−γ)N)). m_N=p (I- 1γ NH^-1 (I-(I- )^N ) ) δ. (G.7) Therefore, 12‖¯N‖2=p22‖(−1γN−1(−(−γ)N))‖2. 12\| m_N\|_H^2= p^22 \| (I- 1γ NH^-1 (I-(I- )^N ) ) δ \|_H^2. (G.8) We next bound 12⟨,~¯N⟩. 12 , D_N . Define ~t:=[~tdrift⊗~tdrift]. D_t:=E[ η^drift_t η^drift_t]. From (G.4), we obtain the recursion ~t=(ℐ−γ)∘~t−1+γ2~tdrift,~0=, D_t=(I-γT) D_t-1+γ^2 ^drift_t, D_0= 0, (G.9) where ~tdrift:=[((−tt⊤)t−1)⊗((−tt⊤)t−1)]. ^drift_t:=E\! [ ((H-x_tx_t )m_t-1 ) ((H-x_tx_t )m_t-1 ) ]. Indeed, the cross term vanishes because [(−tt⊤)t−1]=(−[tt⊤])t−1=.E[(H-x_tx_t )m_t-1]=(H-E[x_tx_t ])m_t-1= 0. Let t:=t−1⊗t−1.A_t:=m_t-1 _t-1. Then ~tdrift ^drift_t =[(−tt⊤)t(−tt⊤)] =E\! [(H-x_tx_t )A_t(H-x_tx_t ) ] =[tt⊤ttt⊤]−t =E[x_tx_t A_tx_tx_t ]-HA_tH =(ℳ−ℳ~)∘t. =(M- M) _t. Since ℳ−ℳ~M- M is a PSD mapping, we have ~tdrift⪰ ^drift_t 0. On the other hand, by Assumption 3.23.2, [tt⊤ttt⊤]=ℳ∘t⪯αtr(t)=α‖t−1‖2.E[x_tx_t A_tx_tx_t ]=M _t α\, tr(HA_t)\,H=α\,\|m_t-1\|_H^2\,H. Therefore, ~tdrift⪯α‖t−1‖2. ^drift_t α\,\|m_t-1\|_H^2\,H. (G.10) Next, by the explicit formula for tm_t, t−1=p(−(−γ)t−1),m_t-1=p (I-(I- )^t-1 ) δ, and hence ‖t−1‖2 \|m_t-1\|_H^2 =p2∑iλi(1−(1−γλi)t−1)2δi2 =p^2 _i _i (1-(1-γ _i)^t-1 )^2 _i^2 ≤p2∑iλiδi2=p2‖2. ≤ p^2 _i _i _i^2=p^2\| δ\|_H^2. Substituting this into (G.10) gives the uniform bound ~tdrift⪯αp2‖2,∀t. ^drift_t α p^2\| δ\|_H^2\,H, ∀ t. (G.11) Now (G.9) has exactly the same form as the variance recursion for constant-stepsize SGD with iterate averaging, with noise level upper bounded by σdrift2:=αp2‖2. _drift^2:=α p^2\| δ\|_H^2. Therefore, applying Lemma G.2 with σ¯2 σ^2 replaced by σdrift2 _drift^2, we obtain 12⟨,~¯N⟩≤αp2‖21−γαtr()⋅DN. 12 , D_N ≤ α p^2\| δ\|_H^21-γα tr(H)· D_ effN. (G.12) Finally, combining (G.6), (G.8) and (G.12), we conclude that 12⟨,¯N⟩ 12 , D_N ≤p22‖(−1γN−1(−(−γ)N))‖2 ≤ p^22 \| (I- 1γ NH^-1 (I-(I- )^N ) ) δ \|_H^2 +αp2‖21−γαtr()⋅DN. + α p^2\| δ\|_H^21-γα tr(H)· D_ effN. This completes the proof. ∎ Fluctuation upper bound Recall that the fluctuation process satisfies tfluct=(−γtt⊤)t−1fluct+γzttt⊤−γp,0fluct=. η^fluct_t=(I- _tx_t ) η^fluct_t-1+γ z_tx_tx_t δ-γ pH δ, η^fluct_0= 0. We decompose the driving term into two parts: zttt⊤−p=zt(tt⊤−)+(zt−p).z_tx_tx_t δ-pH δ=z_t(x_tx_t -H) δ+(z_t-p)H δ. Accordingly, define tfluct=tcov+tsch, η^fluct_t= η^cov_t+ η^sch_t, where tcov=(−γtt⊤)t−1cov+γzt(tt⊤−),0cov=,tsch=(−γtt⊤)t−1sch+γ(zt−p),0sch=. cases η^cov_t=(I- _tx_t ) η^cov_t-1+γ z_t(x_tx_t -H) δ,\\[2.15277pt] η^cov_0= 0, cases cases η^sch_t=(I- _tx_t ) η^sch_t-1+γ(z_t-p)H δ,\\[2.15277pt] η^sch_0= 0. cases Thus, for the averaged iterates, ¯Nfluct=¯Ncov+¯Nsch. η^fluct_N= η^cov_N+ η^sch_N. Hence, by (a+b)2≤2a2+2b2(a+b)^2≤ 2a^2+2b^2, 12⟨,¯N⟩=12[‖¯Nfluct‖2]≤[‖¯Ncov‖2]+[‖¯Nsch‖2]. 12 , F_N = 12E [\| η^fluct_N\|_H^2 ] [\| η^cov_N\|_H^2 ]+E [\| η^sch_N\|_H^2 ]. (G.13) Define tcov:=[tcov⊗tcov].F_t^cov:=E\! [ η^cov_t η^cov_t ]. Then tcovF_t^cov satisfies the recursion tcov=(ℐ−γ)∘t−1cov+γ2tcov,0cov=,F_t^cov=(I-γT) _t-1^cov+γ^2 _t^cov, _0^cov= 0, (G.14) where tcov:=[(zt(tt⊤−))⊗(zt(tt⊤−))]. _t^cov:=E\! [ (z_t(x_tx_t -H) δ ) (z_t(x_tx_t -H) δ ) ]. Indeed, the cross term vanishes as we have proved in last-iterate analysis. Next, since zt∈0,1z_t∈\0,1\ and [zt]=pE[z_t]=p, we have tcov _t^cov =[zt2((tt⊤−))⊗((tt⊤−))] =E\! [z_t^2 ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ] =p[((tt⊤−))⊗((tt⊤−))]. =p\,E\! [ ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ]. Let :=⊗A:= δ δ. Then [((tt⊤−))⊗((tt⊤−))] \! [ ((x_tx_t -H) δ ) ((x_tx_t -H) δ ) ] =[(tt⊤−)(tt⊤−)] =E\! [(x_tx_t -H)A(x_tx_t -H) ] =[tt⊤tt⊤]− =E[x_tx_t Ax_tx_t ]-HAH =(ℳ−ℳ~)∘. =(M- M) . Since ℳ−ℳ~M- M is a PSD mapping, and by Assumption 3.2, ℳ∘=[tt⊤tt⊤]⪯αtr()=α‖2,M =E[x_tx_t Ax_tx_t ] α\, tr(HA)\,H=α\| δ\|_H^2\,H, we obtain the upper bound tcov⪯pα‖2. _t^cov pα\| δ\|_H^2\,H. (G.15) Therefore, the recursion (G.14) has exactly the same form as the variance recursion for constant-stepsize SGD with iterate averaging, with noise level upper bounded by σcov2:=pα‖2. _cov^2:=pα\| δ\|_H^2. Applying Lemma G.2 with σ¯2 σ^2 replaced by σcov2 _cov^2, we get 12[‖¯Ncov‖2]=12⟨,¯Ncov⟩≤pα‖21−γαtr()⋅DN. 12E [\| η^cov_N\|_H^2 ]= 12 , F_N^cov ≤ pα\| δ\|_H^21-γα tr(H)· D_ effN. (G.16) We now bound the second component tsch=(−γtt⊤)t−1sch+γ(zt−p),0sch=. η^sch_t=(I- _tx_t ) η^sch_t-1+γ(z_t-p)H δ, η^sch_0= 0. Its averaged iterate is ¯Nsch:=1N∑t=1Ntsch. η^sch_N:= 1N _t=1^N η^sch_t. Then we use the following lemma to bound ⟨,¯Nsch⊗¯Nsch⟩ , η^sch_N η^sch_N . Lemma G.4 (A schedule-fluctuation upper bound for Avg-SGD). Suppose Assumptions 3.1 and 3.2 hold. Assume the source indicators (z1,…,zN)(z_1,…,z_N) follow the fixed-budget model with exactly m ones, and let p:=mN.p:= mN. Consider the schedule fluctuation recursion tsch=(−γtt⊤)t−1sch+γ(zt−p),0sch=, η^sch_t=(I- _tx_t ) η^sch_t-1+γ(z_t-p)H δ, η^sch_0= 0, and define its averaged iterate by ¯Nsch:=1N∑t=1Ntsch. η^sch_N:= 1N _t=1^N η^sch_t. Suppose γ<1/(αtr())γ<1/(α tr(H)). Let k∗:=maxk:λk≥1γN,D:=k∗+γ2N2∑i>k∗λi2.k^*:= \k: _k≥ 1γ N \, D_ eff:=k^*+γ^2N^2 _i>k^* _i^2. Then 12[‖¯Nsch‖2]≤p(1−p)‖2N−1(1+2αγtr()1−γαtr()D). 12E\! [\| η^sch_N\|_H^2 ]≤ p(1-p)\| δ\|_H^2N-1 (1+ 2αγ tr(H)1-γα tr(H)D_ eff ). (G.17) Proof. Let t:=−γtt⊤.B_t:=I- _tx_t . By repeated substitution of the recursion tsch=tt−1sch+γ(zt−p),0sch=0, η^sch_t=B_t η^sch_t-1+γ(z_t-p)H δ, η^sch_0=0, we obtain, for every t∈1,…,Nt∈\1,…,N\, tsch=γ∑s=1t(∏j=s+1tj)(zs−p), η^sch_t=γ _s=1^t ( _j=s+1^tB_j )(z_s-p)H δ, (G.18) where the empty product is interpreted as the identity operator. Therefore, ¯Nsch η^sch_N =1N∑t=1Ntsch = 1N _t=1^N η^sch_t =γN∑t=1N∑s=1t(∏j=s+1tj)(zs−p) = γN _t=1^N _s=1^t ( _j=s+1^tB_j )(z_s-p)H δ =∑s=1N(zs−p)s, = _s=1^N(z_s-p)v_s, (G.19) where s:=γN∑t=sN(∏j=s+1tj).v_s:= γN _t=s^N ( _j=s+1^tB_j )H δ. (G.20) Condition on the sample sequence (1,…,N)(x_1,…,x_N). Then 1,…,Nv_1,…,v_N are deterministic vectors. Since the feature sequence is independent of the fixed-budget schedule, the indicators (z1,…,zN)(z_1,…,z_N) still follow the same fixed-budget law conditional on (1,…,N)(x_1,…,x_N). Under the fixed-budget model, for every s, [zs]=p,[(zs−p)2]=p(1−p),E[z_s]=p, [(z_s-p)^2]=p(1-p), and for every s≠rs≠ r, [(zs−p)(zr−p)]=−p(1−p)N−1.E[(z_s-p)(z_r-p)]=- p(1-p)N-1. Hence, conditioned on (1,…,N)(x_1,…,x_N), [‖∑s=1N(zs−p)s‖2|1,…,N] \! [ \| _s=1^N(z_s-p)v_s \|_H^2\, |\,x_1,…,x_N ] =∑s=1N∑r=1N[(zs−p)(zr−p)]⟨s,r⟩ = _s=1^N _r=1^NE[(z_s-p)(z_r-p)] _s,v_r _H =p(1−p)∑s=1N‖s‖2−p(1−p)N−1∑s≠r⟨s,r⟩. =p(1-p) _s=1^N\|v_s\|_H^2- p(1-p)N-1 _s≠ r _s,v_r _H. (G.21) Let ¯:=1N∑s=1Ns. v:= 1N _s=1^Nv_s. Using the identity ∑s=1N‖s−¯‖2=∑s=1N‖s‖2−1N‖∑s=1Ns‖2, _s=1^N\|v_s- v\|_H^2= _s=1^N\|v_s\|_H^2- 1N \| _s=1^Nv_s \|_H^2, one checks that the right-hand side of (G.21) equals Np(1−p)N−1∑s=1N‖s−¯‖2≤Np(1−p)N−1∑s=1N‖s‖2. Np(1-p)N-1 _s=1^N\|v_s- v\|_H^2≤ Np(1-p)N-1 _s=1^N\|v_s\|_H^2. Combining this with (G.19), and taking expectation over the samples, we obtain [‖¯Nsch‖2]≤Np(1−p)N−1∑s=1N[‖s‖2].E\! [\| η^sch_N\|_H^2 ]≤ Np(1-p)N-1 _s=1^NE\! [\|v_s\|_H^2 ]. (G.22) Fix s∈1,…,Ns∈\1,…,N\ and let ns:=N−s+1.n_s:=N-s+1. Define an auxiliary bias-only process initialized at γ δ: (s)0:=γ,(s)r:=s+r(s)r−1,r=1,…,ns−1.u^(s)_0:= δ, ^(s)_r:=B_s+ru^(s)_r-1, r=1,…,n_s-1. Then r(s)=(∏j=s+1s+rj)γ.u^(s)_r= ( _j=s+1^s+rB_j ) δ. Therefore, s _s =γN∑t=sN(∏j=s+1tj) = γN _t=s^N ( _j=s+1^tB_j )H δ =1N∑r=0ns−1r(s)=nsN¯ns(s), = 1N _r=0^n_s-1u^(s)_r= n_sN\, u^(s)_n_s, (G.23) where ¯ns(s):=1ns∑r=0ns−1r(s). u^(s)_n_s:= 1n_s _r=0^n_s-1u^(s)_r. Since the samples are i.i.d., the distribution of (0(s),…,ns−1(s))(u^(s)_0,…,u^(s)_n_s-1) coincides with the distribution of a length-nsn_s bias-only SGD trajectory with constant stepsize γ and initialization γ δ. Hence [‖s‖2]=(nsN)2[‖¯ns(γ)‖2],E[\|v_s\|_H^2]= ( n_sN )^2E\! [\| u_n_s( δ)\|_H^2 ], (G.24) where ¯n(γ) u_n( δ) denotes the averaged bias-only SGD iterate of horizon n initialized at γ δ. Substituting (G.24) into (G.22), and changing variables from s to n=nsn=n_s, yields [‖¯Nsch‖2]≤Np(1−p)N−1∑n=1N(nN)2[‖¯n(γ)‖2].E\! [\| η^sch_N\|_H^2 ]≤ Np(1-p)N-1 _n=1^N ( nN )^2E\! [\| u_n( δ)\|_H^2 ]. (G.25) For each n∈1,…,Nn∈\1,…,N\, let kn:=maxk:λk≥1γn,Dn:=kn+n2γ2∑i>knλi2.k_n:= \k: _k≥ 1γ n \, D_n:=k_n+n^2γ^2 _i>k_n _i^2. Applying Lemma G.1 with horizon n and initialization 0=γ,u_0= δ, we obtain 12[‖¯n(γ)‖2]≤ 12E\! [\| u_n( δ)\|_H^2 ]≤\; 1γ2n2∥γ∥−10:kn2+∥γ∥kn:∞2 1γ^2n^2\| δ\|_H^-1_0:k_n^2+\| δ\|_H_k_n:∞^2 +2α(∥γ∥0:kn2+nγ∥γ∥kn:∞2)nγ(1−γαtr())⋅Dnn. + 2α (\| δ\|_I_0:k_n^2+nγ\| δ\|_H_k_n:∞^2 )nγ(1-γα tr(H))· D_nn. (G.26) We now simplify each term in the eigenbasis of H. Writing =∑iδii,γ=γ∑iλiδii, δ= _i _iv_i, δ=γ _i _i _iv_i, we have 1γ2n2∥γ∥−10:kn2 1γ^2n^2\| δ\|_H^-1_0:k_n^2 =1n2∑i≤knλiδi2, = 1n^2 _i≤ k_n _i _i^2, (G.27) ∥γ∥kn:∞2 \| δ\|_H_k_n:∞^2 =γ2∑i>knλi3δi2, =γ^2 _i>k_n _i^3 _i^2, (G.28) ∥γ∥0:kn2 \| δ\|_I_0:k_n^2 =γ2∑i≤knλi2δi2. =γ^2 _i≤ k_n _i^2 _i^2. (G.29) For the first two terms, since λi<1/(γn) _i<1/(γ n) for every i>kni>k_n, γ2λi3≤λin2.γ^2 _i^3≤ _in^2. Therefore, by (G.27) and (G.28), 1γ2n2∥γ∥−10:kn2+∥γ∥kn:∞2 1γ^2n^2\| δ\|_H^-1_0:k_n^2+\| δ\|_H_k_n:∞^2 ≤1n2∑i≤knλiδi2+1n2∑i>knλiδi2 ≤ 1n^2 _i≤ k_n _i _i^2+ 1n^2 _i>k_n _i _i^2 =‖2n2. = \| δ\|_H^2n^2. (G.30) For the third term, using again λi<1/(γn) _i<1/(γ n) for i>kni>k_n, we have nγλi3≤λi2.nγ _i^3≤ _i^2. Hence ∥γ∥0:kn2+nγ∥γ∥kn:∞2 \| δ\|_I_0:k_n^2+nγ\| δ\|_H_k_n:∞^2 =γ2∑i≤knλi2δi2+nγ3∑i>knλi3δi2 =γ^2 _i≤ k_n _i^2 _i^2+nγ^3 _i>k_n _i^3 _i^2 ≤γ2∑iλi2δi2 ≤γ^2 _i _i^2 _i^2 ≤γ2tr()∑iλiδi2=γ2tr()‖2. ≤γ^2 tr(H) _i _i _i^2=γ^2 tr(H)\| δ\|_H^2. (G.31) Substituting (G.30) and (G.31) into (G.26), we obtain 12[‖¯n(γ)‖2]≤‖2n2(1+2αγtr()1−γαtr()Dn). 12E\! [\| u_n( δ)\|_H^2 ]≤ \| δ\|_H^2n^2 (1+ 2αγ tr(H)1-γα tr(H)\,D_n ). (G.32) Substituting (G.32) into (G.25) yields 12[‖¯Nsch‖2] 12E\! [\| η^sch_N\|_H^2 ] ≤Np(1−p)N−1∑n=1N(nN)2⋅‖2n2(1+2αγtr()1−γαtr()Dn) ≤ Np(1-p)N-1 _n=1^N ( nN )^2· \| δ\|_H^2n^2 (1+ 2αγ tr(H)1-γα tr(H)\,D_n ) =p(1−p)‖2N(N−1)∑n=1N(1+2αγtr()1−γαtr()Dn). = p(1-p)\| δ\|_H^2N(N-1) _n=1^N (1+ 2αγ tr(H)1-γα tr(H)\,D_n ). (G.33) It remains to compare DnD_n with the final effective dimension DD_ eff. Observe that Dn=kn+n2γ2∑i>knλi2=∑imin1,n2γ2λi2.D_n=k_n+n^2γ^2 _i>k_n _i^2= _i \1,n^2γ^2 _i^2\. Since n≤Nn≤ N, Dn≤∑imin1,N2γ2λi2.D_n≤ _i \1,N^2γ^2 _i^2\. Now let k∗:=maxk:λk≥1γN.k^*:= \k: _k≥ 1γ N \. Then ∑imin1,N2γ2λi2≤k∗+γ2N2∑i>k∗λi2=D. _i \1,N^2γ^2 _i^2\≤ k^*+γ^2N^2 _i>k^* _i^2=D_ eff. Hence Dn≤D,∀n≤N.D_n≤D_ eff, ∀ n≤ N. Applying this to (G.33) gives 12[‖¯Nsch‖2] 12E\! [\| η^sch_N\|_H^2 ] ≤p(1−p)‖2N(N−1)∑n=1N(1+2αγtr()1−γαtr()D) ≤ p(1-p)\| δ\|_H^2N(N-1) _n=1^N (1+ 2αγ tr(H)1-γα tr(H)\,D_ eff ) =p(1−p)‖2N−1(1+2αγtr()1−γαtr()D), = p(1-p)\| δ\|_H^2N-1 (1+ 2αγ tr(H)1-γα tr(H)\,D_ eff ), which proves (G.17). ∎ An upper bound for excess risk. Combining the bounds for the bias, variance, drift, and fluctuation terms derived above, we obtain the following upper bound on the excess risk. Theorem G.5 (Excess risk upper bound for constant-stepsize SGD with iterate averaging). Suppose Assumptions 3.1, 3.2, and 3.3 hold. Let σ¯2:=(1−p)σ12+pσ22. σ^2:=(1-p) _1^2+p _2^2. Assume γ<1/(αtr())γ<1/(α tr(H)). Then [ℛ(¯N)−ℛ(1∗)] [R( w_N)-R(w_1^*)] ≲1γ2N2∥0−1∗∥0:k∗−12+∥0−1∗∥k∗:∞2 1γ^2N^2\|w_0-w_1^*\|_H_0:k^*^-1^2+\|w_0-w_1^*\|_H_k^*:∞^2 +(α∥0−1∗∥2+σ¯2+αp∥2)⋅DN+p(1−p)‖2N, + (α\|w_0-w_1^*\|_H^2+ σ^2+α p\| δ\|_H^2 )· D_ effN+ p(1-p)\| δ\|_H^2N, +p2‖(−1γN−1(−(−γ)N))‖2 +p^2 \| (I- 1γ NH^-1 (I-(I- )^N ) ) δ \|_H^2 where k∗=maxk:λk≥1γN,k^*= \k: _k≥ 1γ N \, and D:=k∗+γ2N2∑i>k∗λi2.D_ eff:=k^*+γ^2N^2 _i>k^* _i^2. G.2 Lower bound analysis The starting point is again the error recursion t=(−γtt⊤)t−1+γzttt⊤+γξtt,0=0−1∗, η_t=(I- _tx_t ) η_t-1+γ z_tx_tx_t δ+γ _tx_t, η_0=w_0-w_1^*, (G.34) where t=t−1∗. η_t=w_t-w_1^*. Recall that the averaged iterate is ¯N:=1N∑t=0N−1t, w_N:= 1N _t=0^N-1w_t, and correspondingly ¯N:=1N∑t=0N−1t=¯N−1∗. η_N:= 1N _t=0^N-1 η_t= w_N-w_1^*. By definition of the excess risk, [ℛ(¯N)−ℛ(1∗)]=12⟨,[¯N⊗¯N]⟩.E[R( w_N)-R(w_1^*)]= 12\, ,E[ η_N η_N] . To derive a lower bound, we center the averaged error around its mean. Define ¯N:=[¯N],¯~N:=¯N−¯N. μ_N:=E[ η_N], η_N:= η_N- μ_N. Then ¯N=¯N+¯~N,[¯~N]=. η_N= μ_N+ η_N, [ η_N]= 0. Therefore, [¯N⊗¯N] [ η_N η_N] =[(¯N+¯~N)⊗(¯N+¯~N)] =E [( μ_N+ η_N) ( μ_N+ η_N) ] =¯N⊗¯N+[¯~N⊗¯~N]+¯N⊗[¯~N]+[¯~N]⊗¯N = μ_N μ_N+E[ η_N η_N]+ μ_N [ η_N]+E[ η_N] μ_N =¯N⊗¯N+[¯~N⊗¯~N]. = μ_N μ_N+E[ η_N η_N]. Hence [ℛ(¯N)−ℛ(1∗)]=12‖¯N‖2+12⟨,[¯~N⊗¯~N]⟩.E[R( w_N)-R(w_1^*)]= 12\| μ_N\|_H^2+ 12 ,E[ η_N η_N] . (G.35) Since the second term is nonnegative, we immediately obtain [ℛ(¯N)−ℛ(1∗)]≥12‖¯N‖2.E[R( w_N)-R(w_1^*)]≥ 12\| μ_N\|_H^2. (G.36) We now characterize the mean trajectory of the SGD error process. Let t:=[t]. μ_t:=E[ η_t]. Taking expectation on both sides of (G.34), and using [tt⊤]=,[zt]=p,[ξtt]=,E[x_tx_t ]=H, [z_t]=p, [ _tx_t]= 0, we obtain t=(−γ)t−1+γp,0=0−1∗. μ_t=(I- ) μ_t-1+γ pH δ, μ_0=w_0-w_1^*. Unrolling the recursion gives t μ_t =(−γ)t(0−1∗)+p∑s=1tγ(−γ)t−s. =(I- )^t(w_0-w_1^*)+p _s=1^tγ(I- )^t-sH δ. (G.37) Since (−γ)(I- ) is a polynomial in H, it commutes with H. Therefore, ∑s=1tγ(−γ)t−s=−(−γ)t. _s=1^tγ(I- )^t-sH=I-(I- )^t. Substituting this identity into (G.37), we obtain the closed form t=(−γ)t(0−1∗)+p(−(−γ)t). μ_t=(I- )^t(w_0-w_1^*)+p (I-(I- )^t ) δ. Now average over t=0,…,N−1t=0,…,N-1. Since ¯N=[¯N]=1N∑t=0N−1t, μ_N=E[ η_N]= 1N _t=0^N-1 μ_t, we have ¯N μ_N =1N∑t=0N−1(−γ)t(0−1∗)+pN∑t=0N−1(−(−γ)t). = 1N _t=0^N-1(I- )^t(w_0-w_1^*)+ pN _t=0^N-1 (I-(I- )^t ) δ. Using the matrix geometric-series identity ∑t=0N−1(−γ)t=(γ)−1(−(−γ)N), _t=0^N-1(I- )^t=( )^-1 (I-(I- )^N ), we obtain ¯N=1γN−1(−(−γ)N)(0−1∗)+p(−1γN−1(−(−γ)N)). μ_N= 1γ NH^-1 (I-(I- )^N )(w_0-w_1^*)+p (I- 1γ NH^-1 (I-(I- )^N ) ) δ. Consequently, ‖¯N‖2=‖1γN−1(−(−γ)N)(0−1∗)+p(−1γN−1(−(−γ)N))‖2.\| μ_N\|_H^2= \| 1γ NH^-1 (I-(I- )^N )(w_0-w_1^*)+p (I- 1γ NH^-1 (I-(I- )^N ) ) δ \|_H^2. (G.38) Combining (G.36) and (G.38), we arrive at the lower bound [ℛ(¯N)−ℛ(1∗)] [R( w_N)-R(w_1^*)] ≳‖1γN−1(−(−γ)N)(0−1∗)+p(−1γN−1(−(−γ)N))‖2. \| 1γ NH^-1 (I-(I- )^N )(w_0-w_1^*)+p (I- 1γ NH^-1 (I-(I- )^N ) ) δ \|_H^2. G.3 Strong model collapse behavior Corollary G.6 (Strong model collapse in Avg-SGD). Consider averaged SGD with constant stepsize γ. Suppose Assumptions 3.1, 3.2 and 3.3 hold. Suppose γ<1/αtr()γ<1/α tr(H) and D=o(M+N)D_ eff=o(M+N). Further assume that ‖0−1∗‖22\|w_0-w_1^*\|_2^2 is finite. As the total sample size M+NM+N scales to infinity with a fixed synthetic proportion p, limM+N→∞[ℰ1(¯M+N)]≂p2‖2. _M+N ∞E[E_1( w_M+N)] p^2\| δ\|_H^2. Proof. Let T:=M+N,T:=1γT−1(−(−γ)T).T:=M+N, _T:= 1γ TH^-1 (I-(I- )^T ). Since the synthetic proportion p=M/(M+N)p=M/(M+N) is fixed, it suffices to study the limit as T→∞T→∞. We first show that TA_T vanishes in H-norm. Writing the eigendecomposition =∑iλiii⊤,H= _i _iv_iv_i , we have Ti=1−(1−γλi)TγTλii.A_Tv_i= 1-(1-γ _i)^Tγ T _i\,v_i. Since γ<1/(αtr())γ<1/(α tr(H)), we have 0≤γλi≤γtr()<10≤γ _i≤γ tr(H)<1, hence 0≤1−γλi<1.0≤ 1-γ _i<1. For every i with λi>0 _i>0, it follows that 1−(1−γλi)TγTλi⟶0as T→∞. 1-(1-γ _i)^Tγ T _i 0 T→∞. Moreover, using the elementary inequality 1−(1−a)T≤Ta,0≤a≤1,1-(1-a)^T≤ Ta, 0≤ a≤ 1, we get 0≤1−(1−γλi)TγTλi≤1.0≤ 1-(1-γ _i)^Tγ T _i≤ 1. Therefore, for any w with ‖<∞\|w\|_H<∞, ‖T‖2=∑iλi(1−(1−γλi)TγTλi)2⟨,i⟩2⟶0\|A_Tw\|_H^2= _i _i ( 1-(1-γ _i)^Tγ T _i )^2 ,v_i ^2 0 by dominated convergence. We now turn to the upper bound. By the upper bound theorem for averaged SGD, [ℰ1(¯T)] [E_1( w_T)] ≲1γ2T2∥0−1∗∥0:k∗−12⏟(I)+∥0−1∗∥k∗:∞2⏟(II) 1γ^2T^2\|w_0-w_1^*\|_H_0:k^*^-1^2_(I)+ \|w_0-w_1^*\|_H_k^*:∞^2_(I) +(α‖0−1∗‖2+σ¯2+αp‖2)⋅DT⏟(III)+p(1−p)‖2T⏟(IV) + (α\|w_0-w_1^*\|_H^2+ σ^2+α p\| δ\|_H^2 )· D_ effT_(I)+ p(1-p)\| δ\|_H^2T_(IV) +p2‖(−T)‖2⏟(V), + p^2\|(I-A_T) δ\|_H^2_(V), where k∗=maxk:λk≥1γT,D:=k∗+γ2T2∑i>k∗λi2.k^*= \k: _k≥ 1γ T \, D_ eff:=k^*+γ^2T^2 _i>k^* _i^2. By the assumption D=o(T)D_ eff=o(T), terms (III)(I) and (IV)(IV) vanish as T→∞T→∞. For term (I)(I), note that 1γ2T2∥0−1∗∥0:k∗−12≤∥0−1∗∥22⋅1γT, 1γ^2T^2\|w_0-w_1^*\|_H_0:k^*^-1^2≤\|w_0-w_1^*\|_2^2· 1γ T, because for every i≤k∗i≤ k^* we have λi≥1/(γT) _i≥ 1/(γ T), hence λi−1≤γT _i^-1≤γ T. Therefore (I)→0(I)→ 0. For term (II)(I), since k∗→∞k^*→∞ as T→∞T→∞ and ‖0−1∗‖<∞\|w_0-w_1^*\|_H<∞, we have ∥0−1∗∥k∗:∞2=∑i>k∗λi⟨0−1∗,i⟩2⟶0.\|w_0-w_1^*\|_H_k^*:∞^2= _i>k^* _i _0-w_1^*,v_i ^2 0. Finally, ‖(−T)‖→‖,\|(I-A_T) δ\|_H→\| δ\|_H, because ‖(−T)−‖=‖T‖→0.\|(I-A_T) δ- δ\|_H=\|A_T δ\|_H→ 0. Therefore, lim supT→∞[ℰ1(¯T)]≲p2‖2. _T→∞E[E_1( w_T)] p^2\| δ\|_H^2. Next, we apply the lower bound obtained from the mean term: [ℰ1(¯T)]≳‖T(0−1∗)+p(−T)‖2.E[E_1( w_T)] \|A_T(w_0-w_1^*)+p(I-A_T) δ \|_H^2. Rewrite the term inside the norm as T(0−1∗)+p(−T)=T(0−1∗−p)+p.A_T(w_0-w_1^*)+p(I-A_T) δ=A_T(w_0-w_1^*-p δ)+p δ. Again using ‖T‖→0\|A_Tw\|_H→ 0 for every w with finite H-norm, we obtain ‖T(0−1∗−p)+p‖2⟶p2‖2. \|A_T(w_0-w_1^*-p δ)+p δ \|_H^2 p^2\| δ\|_H^2. Hence lim infT→∞[ℰ1(¯T)]≳p2‖2. _T→∞E[E_1( w_T)] p^2\| δ\|_H^2. Combining the upper and lower bounds yields limM+N→∞[ℰ1(¯M+N)]≂p2‖2. _M+N→∞E[E_1( w_M+N)] p^2\| δ\|_H^2. This proves the corollary. ∎ Appendix H Additional Experiments H.1 CIFAR10 CNN Experiments CIFAR10 experiment setup. We conduct experiments on the CIFAR10 [Krizhevsky et al., 2009] dataset, which consists of 50,000 training images and 10,000 test images across 10 classes. All images are of size 32×3232× 32. All experiments are run on a single NVIDIA RTX 3090 GPU. We use a CNN with five convolutional layers. Each block consists of a convolution layer followed by batch normalization and ReLU activation. Model size is controlled by a width parameter D, which scales the number of channels in all convolutional layers proportionally. Specifically, the channel dimensions follow the pattern (D,D,2D,2D,4D)(D,D,2D,2D,4D). We vary D in 8,16,32,64,128,256\8,16,32,64,128,256\ to study the effect of model capacity. Synthetic data shares the same input distribution as the real data, but labels are generated by a weak teacher model. The teacher is a smaller CNN (width 1616) trained on a random subset of 5,000 training samples for 5 epochs. For each input x, the teacher produces a soft label distribution qteacher(y∣x)q_teacher(y x) with temperature scaling. The synthetic label is then defined as y~=qteacher(y∣x), y=q_teacher(y x), , corresponding to fully synthetic labels. Similar to the simulation experiments, we consider two training protocols as described in Section 3. For a fixed total training budget T, we allocate pTpT synthetic samples and (1−p)T(1-p)T real samples, where p∈[0,1]p∈[0,1] is the synthetic data proportion. We vary p over 0.0,0.1,…,0.9\0.0,0.1,…,0.9\ and T over a grid ranging from 200200 to 20,00020,000. To examine how the synthetic-induced floor depends on model size, we fix p=0.9p=0.9 and vary the model width D (while also varying T across the same grid). Evaluation. We evaluate performance on the real test set using both classification error and cross-entropy loss. Each experiment is repeated with 3 random seeds, and we report the mean and standard deviation. To isolate the effect of synthetic data from baseline model performance, we report the metrics difference relative to p=0p=0 under the same T: Δerror=Error(p=0.9)−Error(p=0),Δloss=Loss(p=0.9)−Loss(p=0), _error=Error(p=0.9)-Error(p=0), _loss=Loss(p=0.9)-Loss(p=0), which removes the intrinsic advantage of larger models and highlights the synthetic-data-induced degradation. Figure 3: Real-data experiments of different training protocols. Test cross-entropy loss on the real distribution as a function of the number of training samples T, under different synthetic proportions p. Left: random mixing of real and synthetic data. Right: two-stage training (synthetic-data training followed by real-data training). From Figure 3, we observe the same qualitative behavior predicted by our theory. Under random mixing, the test loss quickly saturates as T increases, forming a clear floor that depends on the synthetic proportion p. In particular, larger values of p lead to systematically higher asymptotic error, indicating that synthetic data introduces an irreducible bias that cannot be removed by additional samples. In contrast, two-stage training mitigates this phenomenon: the test loss continues to decrease with T, and consistently achieves slightly better performance than mixing in the large-sample regime. These results confirm that the synthetic-data-induced floor is a consequence of mixing, rather than an inherent limitation of the data itself. (a) Test error difference (b) Test cross-entropy difference Figure 4: Effect of model size on synthetic-data-induced degradation. We plot the performance gap between high synthetic proportion (p=0.9p=0.9) and pure real-data training (p=0p=0), as a function of the number of samples T, for different model widths. Across both test error and cross-entropy loss, the gap consistently increases with model size, especially in the large-sample regime. This indicates that larger models are more sensitive to synthetic data, leading to a stronger irreducible bias under mixing. Figure 4 shows how the impact of synthetic data varies with model size on real data. We measure the performance gap between high synthetic proportion (p=0.9p=0.9) and pure real-data training (p=0p=0), thereby isolating the effect of synthetic data. We observe a clear and consistent trend: as model size increases, the gap becomes significantly larger, particularly at large sample sizes. While all models exhibit some degradation due to synthetic data, larger models suffer substantially more, indicating that the synthetic-data-induced bias is amplified with model capacity. This behavior aligns with our theoretical predictions, suggesting that the floor effect induced by mixing is not only persistent but also worsens with increasing model size. H.2 Language Model Experiments Experiment setup. We conduct language-modeling experiments on WikiText-103, tokenized with the GPT-2 tokenizer. All models are decoder-only Transformers trained with context length 256 using next-token cross-entropy loss, and evaluated by validation loss on the real validation set. Unless otherwise stated, we use batch size 64 and AdamW with learning rate 6×10−46× 10^-4, weight decay 0.10.1, training on a single NVIDIA RTX 3090 GPU. Synthetic text is generated from a DistilGPT-2 teacher. For the synthetic corpus, we sample continuations of length 256 from real WikiText-103 prompts of length 32, using temperature 1.51.5, top-p=0.95p=0.95, top-k=0k=0, and store only the generated continuations. We first compare the same two training protocols as before. In mixed training, each step uses real or synthetic data so that the cumulative synthetic fraction tracks p. In two-stage training, for each total budget T, the model is first trained for approximately pTpT synthetic steps and then for the remaining real steps. Our main experiment uses a 4-layer Transformer with 4 attention heads and embedding dimension 256, containing 16.1M parameters. We vary p∈0,0.3,0.6,0.9p∈\0,0.3,0.6,0.9\ and evaluate 12 logarithmically spaced budgets from 100 to 24,000 steps, corresponding to up to 393M trained tokens. Results are shown in Figure 5. Figure 5: Language-model protocol comparison. We compare mixed training and two-stage training on WikiText-103 for synthetic fractions p∈0.3,0.6,0.9p∈\0.3,0.6,0.9\, with the real-only baseline shown in gray. Under mixed training, validation loss remains consistently worse than the real-only curve, especially at large p, indicating a persistent degradation from training on a fixed fraction of synthetic data. Two-stage training reduces this degradation and approaches the real-only curve more closely at large budgets. This matches the theoretical prediction that mixed training induces a non-vanishing synthetic-data bias floor, whereas using synthetic data only as an initialization stage can mitigate this floor. To study the effect of model capacity, we fix a high synthetic fraction p=0.9p=0.9 under the mixed protocol and compare against the matched real-only baseline p=0p=0. We train three Transformer models with sizes 4/256/4 (16.1M params), 6/384/6 (30.0M params), and 8/512/8 (51.1M params), where the three numbers denote layers, embedding dimension, and attention heads. We evaluate budgets from 3,000 to 48,000 steps with batch size 32, corresponding to 24.6M to 393.2M trained tokens. We report the degradation Δloss=ValLoss(p=0.9)−ValLoss(p=0), =ValLoss(p=0.9)-ValLoss(p=0), which isolates the synthetic-data-induced effect from the baseline improvement due to model size. Results are shown in Figure 6. Figure 6: Effect of model size on synthetic-data-induced degradation. We plot the validation-loss gap Δloss=ValLoss(p=0.9)−ValLoss(p=0) =ValLoss(p=0.9)-ValLoss(p=0) for models with 16.1M, 30.0M, and 51.1M parameters. The larger models exhibit a larger degradation gap at small and intermediate budgets, showing that the effect of synthetic mismatch is amplified by model capacity. This supports the theoretical prediction that higher-capacity models can fit the biased synthetic signal more strongly, leading to a larger synthetic-induced degradation under mixed training. Finally, we test how teacher quality affects the usefulness of synthetic pretraining under a fixed real-data budget. We use the same 4-layer student and compare three synthetic teachers: a high-quality DistilGPT-2 corpus sampled with temperature 0.90.9, top-p=0.95p=0.95, top-k=0k=0, and prompt length 32; the medium-quality corpus from the main experiment; and a low-quality corpus sampled with temperature 4.04.0, top-p=1.0p=1.0, top-k=0k=0, and prompt length 4. We pretrain on synthetic budgets of approximately 25M, 50M, 100M, and 200M tokens, then reset the optimizer and train on real-data budgets up to 50M tokens. Performance is evaluated by real validation loss and by the gap relative to the real-only curve at the same real-token budget. Figure 7. (a) Validation loss after synthetic pretraining followed by fixed real-data training. (b) Validation-loss gap relative to the real-only curve at the same real-token budget. Figure 7: Effect of synthetic teacher quality under a fixed real-data budget. We compare high-, medium-, and low-quality synthetic corpora while keeping the downstream real-data budget fixed. High-quality synthetic pretraining consistently improves validation loss relative to real-only training, and the improvement grows with the synthetic pretraining budget. Medium-quality synthetic data provides weaker but still positive gains at larger synthetic budgets. In contrast, low-quality synthetic data leads to a positive loss gap, meaning that it hurts performance even after subsequent real-data training. These results validate the theoretical bias-quality trade-off: synthetic data can help when its induced bias is small enough to reduce optimization error, but poor-quality synthetic data introduces a mismatch that dominates and causes degradation.