Year: 2021

Abstract:
Generalization is one of the fundamental issues in machine learning. However, traditional techniques like uniform convergence may be unable to explain generalization under overparameterization (Nagarajan & Kolter, 2019). As alternative approaches, techniques based on stability analyze the training dynamics and drive algorithm-dependent generalization bounds. Unfortunately, the stability-based bounds are still far from explaining the surprising generalization in deep learning since neural networks usually suffer from unsatisfactory stability. This paper proposes a novel decomposition framework to improve the stability-based bounds via a more fine-grained analysis of the signal and noise, inspired by the observation that neural networks converge relatively slowly when fitting noise (which indicates better stability). Concretely, we decompose the excess risk dynamics and apply stability-based bound only on the noise component. The decomposition framework performs well in both linear regimes (overparameterized linear regression) and non-linear regimes (diagonal matrix recovery). Experiments on neural networks verify the utility of the decomposition framework.

heading: 1 INTRODUCTION
text: Generalization is one of the essential mysteries uncovered in modern machine learning (Neyshabur et al., 2014; Zhang et al., 2016; Kawaguchi et al., 2017), measuring how the trained model performs on unseen data. One of the most popular approaches to generalization is uniform convergence (Mohri et al., 2018), which takes supremum over parameter space to decouple the dependency between the training set and the trained model. However, Nagarajan & Kolter (2019) pointed out that uniform convergence itself might not be powerful enough to explain generalization, because the uniform bound can still be vacuous under overparameterized linear regimes even if we reduce the parameter space to the minimum.
One alternative solution beyond uniform convergence is to analyze the generalization dynamics, which measures the generalization gap during the training dynamics. Stability-based bound is among the most popular techniques in generalization dynamics analysis (Lei & Ying, 2020), which is derived from algorithmic stability (Bousquet & Elisseeff, 2002). Fortunately, one can derive nonvacuous bounds under general convex regimes using stability frameworks (Hardt et al., 2016). However, stability is still far from explaining the remarkable generalization abilities of neural networks, mainly due to two obstructions. Firstly, stability-based bound depends heavily on the gradient norm in non-convex regimes (Li et al., 2019), which is typically large at the beginning phase in training neural networks. Secondly, stability-based bound usually does not work well under general nonconvex regimes (Hardt et al., 2016; Charles & Papailiopoulos, 2018; Zhou et al., 2018b) but neural networks are usually highly non-convex.
The aforementioned two obstructions mainly stem from the coarse-grained analysis of the signal and noise. As Zhang et al. (2016) argued, neural networks converge fast when fitting signal but converge relatively slowly when fitting noise1, indicating that the training dynamics over signal and noise are significantly different. Consequently, on the one hand, the fast convergence of signal-related training contributes to a large gradient norm at the beginning phase (see Figure 1a), resulting in
1In this paper, we refer the signal to the clean data without the output noise, and the noise to the output noise. See Section 2 for the formal definitions.
poor stability. On the other hand, the training on signal forces the trained parameter away from the initialization, making the whole training path highly non-convex (see Figure 1b). The above two phenomena inspire us to decompose the training dynamics into noise and signal component and only apply the stability-based analysis over the noise component. To demonstrate that such decomposition generally holds in practice, we conduct several experiments of neural networks on both synthetic dataset and real-world dataset (see Figure 2 for more details).
Based on the above discussion, we improve the stability-based analysis by proposing a decomposition framework on excess risk dynamics2, where we handle the noise and signal components separately via bias-variance decomposition. In detail, we decompose the excess risk into variance excess risk (VER) and bias excess risk (BER), where VER measures how the model fits noise and BER measures how the model fits signal. Under the decomposition, we apply the stability-based techniques to VER and apply uniform convergence to BER inspired by Negrea et al. (2020). The decomposition framework accords with the theoretical and experimental evidence surprisingly well, providing that it outperforms stability-based bounds in both linear (overparameterized linear regression) and non-linear (diagonal matrix recovery) regimes.
We summarize our contributions as follows:
• We propose a new framework aiming at improving the traditional stability-based bounds, which is a novel approach to generalization dynamics analysis focusing on decomposing the excess risk dynamics into variance component and bias component. Starting from the overparameterized linear regression regimes, we show how to deploy the decomposition framework in practice, and the proposed framework outperforms the stability-based bounds.
• We theoretically analyze the excess risk decomposition beyond linear regimes. As a case study, we derive a generalization bound under diagonal matrix recovery regimes. To our best knowledge, this is the first work to analyze the generalization performance of diagonal matrix recovery.
• We conduct several experiments on both synthetic datasets and real-world datasets (MINIST, CIFAR-10) to validate the utility of the decomposition framework, indicating that the framework provides interesting insights into the generalization community.
2We decompose the excess risk, which is closely related to generalization, purely due to technical reasons. The excess risk dynamics tracks the excess risk during the training process.

heading: 1.1 RELATED WORK
text: Stability-based Generalization. The stability-based researches can be roughly split into two branches. One branch is about how algorithmic stability leads to generalization (Feldman & Vondrak, 2018; 2019; Bousquet et al., 2020). Another branch focus on how to calculate the stability parameter for specific problems, e.g., Hardt et al. (2016) prove a generalization bound scales linearly with time in convex regimes. Furthermore, researchers try to apply stability techniques into more general settings, e.g., non-smooth loss (Lei & Ying, 2020; Bassily et al., 2020), noisy gradient descent in non-convex regimes (Mou et al., 2018; Li et al., 2019) and stochastic gradient descent in non-convex regimes (Zhou et al., 2018b; Charles & Papailiopoulos, 2018; Zhang et al., 2021). In this paper, we mainly focus on applying the decomposition framework to improve the stability-based bound.
Uniform Convergence is widely used in generalization analysis. For bounded losses, the generalization gap is tightly bounded by its Rademacher Complexity (Koltchinskii & Panchenko, 2000; Koltchinskii, 2001; Koltchinskii et al., 2006). Furthermore, we reach faster rates under realizable assumption (Srebro et al., 2010). A line of work focuses on uniform convergence under neural network regimes, which is usually related to the parameter norm (Bartlett et al., 2017; Wei & Ma, 2019). However, as Nagarajan & Kolter (2019) pointed out, uniform convergence may be unable to explain generalization. Therefore, more techniques are explored to go beyond uniform convergence.
Other Approaches to Generalization. There are some other approaches to generalization, including PAC-Bayes (Neyshabur et al., 2017a; Dziugaite & Roy, 2017; Neyshabur et al., 2017b; Dziugaite & Roy, 2018; Zhou et al., 2018a; Yang et al., 2019), information-based bound (Russo & Zou, 2016; Xu & Raginsky, 2017; Banerjee & Montúfar, 2021; Haghifam et al., 2020; Steinke & Zakynthinou, 2020), and compression-based bound (Arora et al., 2018; Allen-Zhu et al., 2018; Arora et al., 2019).
Bias-Variance Decomposition. Bias-variance decomposition plays an important role in statistical analysis (Lehmann & Casella, 2006; Casella & Berger, 2021; Geman et al., 1992). Generally, high bias indicates that the model has poor predicting ability on average and high variance indicates that the model performs unstably. Bias-variance decomposition is widely used in machine learning analysis, e.g., adversarial training (Yu et al., 2021), double descent (Adlam & Pennington, 2020), uncertainty (Hu et al., 2020). Oymak et al. (2019) applied bias-variance decomposition on Jacobian of neural networks to explain their different performances on clean and noisy data. This paper considers a slightly different bias-variance decomposition following the analysis of SGD (Dieuleveut et al., 2016; Jain et al., 2018; Zou et al., 2021), focusing on the decomposition of the noisy output.
Matrix Recovery. Earlier methods for solving matrix recovery problems rely on the convex relaxation techniques for minimum norm solutions (Recht et al., 2010; Chandrasekaran et al., 2011). Recently, a branch of works focuses on the matrix factorization techniques with simple local search methods. It has been shown that there is no spurious local minima in the exact regime (Ge et al., 2016; 2017; Zhang et al., 2019). In the overparameterized regimes, it was first conjectured by Gunasekar et al. (2018) and then answered by Li et al. (2018) that gradient descent methods converge to the low-rank solution efficiently. Later, Zhuo et al. (2021) extends the conclusion to the noisy settings.

heading: 2 PRELIMINARY
text: In this section, we introduce necessary definitions, assumptions, and formally introduce previous techniques in dealing with generalization.
Data distribution. Let x ∈ X ⊂ Rp be the input and y ∈ Y ⊂ R be the output, where x,y are generated from a joint distribution (x,y) ∼ P . Define Px,Py , and Py|x as the corresponding marginal and conditional distributions, respectively. Given n training samples D ≜ {xi,yi}i∈[n] generated from distribution P , we denote its empirical distribution by Pn. To simplify the notations, we define X ∈ Rn×p as the design matrix and Y ∈ Rn as the response vector. Excess Risk. Given the loss function ℓ(θ;x,y) with parameter θ and sample (x,y), we define the population loss as L(θ;P) ≜ E(x,y)∼P [ℓ(θ;x,y)] and its corresponding training loss as L(θ;Pn) ≜ 1n ∑ i ℓ(θ;xi,yi). Let At denote the optimization algorithm which takes dataset D as input and return the trained parameter θ̂(t) at time t, namely, At(D) = θ̂(t). During the analysis, we focus on the excess risk dynamics EL(θ̂(t);P), which measures how the trained parameter θ̂(t) performs on the population loss:
EL(θ̂(t);P) ≜ L(θ̂(t);P)−min θ L(θ;P).
Although the minimizer can be not unique, we define L(θ∗;Pn) ≜ minθ L(θ;P) where θ∗ denotes arbitrarily one of the minimizers. We additionally remark that bounding the generalization gap L(θ̂(t);P)−L(θ̂(t);Pn) suffices to bound the excess risk under Empirical Risk Minimization (ERM) framework by the following split (see Equation 1). Therefore, we mainly discuss the excess risk following Bartlett et al. (2020).
EL(θ̂(t);P) = [ L(θ̂(t);P)− L(θ̂(t);Pn) ] ︸ ︷︷ ︸
generalization gap
+ [ L(θ̂(t);Pn)− L(θ∗;Pn) ] ︸ ︷︷ ︸
≤0 under ERM
+ [ L(θ∗;Pn)− L(θ∗;P) ] ︸ ︷︷ ︸ ≈0 by concentration inequalities .
(1)
VER and BER. As discussed in Section 1, we aim to decompose the excess risk into the variance excess risk (VER) and the bias excess risk (BER) defined in Definition 1. The decomposition focuses on the noisy output regimes and we split the output y into the signal component E[y|x] and the noise component y − E[y|x]. Informally, VER measures how the model performs on pure noise, and BER measures how the model performs on clean data.
Definition 1 (VER and BER) Given (x,y) ∼ P , let (x,E[y|x]) ∼ Pb denote the signal distribution and (x,y − E[y|x]) ∼ Pv denote the noise distribution. The variance excess risk (VER) EvL(θ;P) and bias excess risk (BER) EbL(θ;P) are defined as:
EvL(θ;P) ≜ EL(θ;Pv), EbL(θ;P) ≜ EL(θ;Pb).
To better illustrate VER and BER, we consider three surrogate training dynamics: Standard Training, Variance Training, and Bias Training, corresponding to ER, VER, and BER, respectively. Standard Training: Training process over the noisy data (X,Y ) from the initialization θ(0). We denote the trained parameter at time t by θ̂(t). Variance Training: Training process over the pure noise (X,Y − E[Y |X]) from the initialization θ (0) v . We denote the trained parameter at time t by θ̂ (t) v .
Bias Training: Training process over the clean data (X,E[Y |X]) from the initialization θ(0)b . We denote the trained parameter at time t by θ̂(t)b . When the context is clear, although the trained parameters θ̂(t), θ̂(t)v , θ̂ (t) b are related to the corresponding initialization and the algorithms, we omit the dependency. Besides, we denote θ∗, θ∗v and θ∗b the optimal parameters which minimize the corresponding population loss L(θ;P), L(θ;Pv), and L(θ;Pb), respectively. Techniques in generalization. We next introduce two techniques in generalization analysis including stability-based techniques in Proposition 1 and uniform convergence in Proposition 2. These techniques will be revisited in Section 3.
Proposition 1 (Stability Bound from Feldman & Vondrak (2019)) Assume that algorithm At is ϵ-uniformly-stable at time t, meaning that for any two dataset D and D′ with only one different data point, we have
sup (x,y)
EA [ℓ(At(D);x,y)− ℓ(At(D′);x,y)] ≤ ϵ.
Then the following inequality holds3 with probability at least 1− δ:
|L (At(D);P)− L (At(D);Pn)| = O ( ϵ log(n) log(n/δ) + √ log(1/δ)√
n
) .
Proposition 2 (Uniform Convergence from Wainwright (2019)) Uniform convergence decouples the dependency between the trained parameter and the training set by taking supremum over a parameter space that is independent of training data, namely
L(At(D);P)− L(At(D);Pn) ≤ sup θ∈B [L(θ;P)− L(θ;Pn)] .
where B is independent of dataset D and At(D) ∈ B for any time t.

heading: 3 WARM-UP: DECOMPOSITION UNDER LINEAR REGIMES
text: In this section, we consider the overparameterized linear regression with linear assumptions, where the sample size is less than the data dimension, namely, n < p. We will show how the decomposition framework works and the improvements compared to traditional stability-based bounds. Before diving into the details, we emphasize the importance of studying overparameterized linear regression. As the arguably simplest models, we expect a framework can at least work under such regimes. Besides, recent works (e.g., Jacot et al. (2018)) show that neural networks converge to overparametrized linear models (with respect to the parameters) as the width goes to infinity under some regularity conditions. Therefore, studying overparameterized linear regression provides enormous intuition to neural networks.
Settings. When the context is clear, we reuse the notations in Section 2. Without loss of generality, we assume that E[x] = 0. Besides, we assume a linear ground truth y = x⊤θ∗+ϵwhere E[ϵ|x] = 0. Let Σx ≜ E[xx⊤] denote the covariance matrix. Given ℓ2 loss ℓ(θ;x,y) = (y − x⊤θ)2, we apply gradient descent (GD) with constant step size λ and initialization θ(0) = 0, namely, θ̂(t+1) = θ̂(t) + λn ∑ i xi(yi − x⊤i θ̂(t)). We provide the main results of overparameterized linear regression in Theorem 1 and then introduce how to apply the decomposition framework.
Theorem 1 (Linear Regression Regimes) Under the overparameterized linear regression settings, assume that w = θ
∗,⊤x√ θ∗,⊤Σxθ∗ is σ2w-subgaussian. If ∥x∥ ≤ 1, |ϵ| ≤ V , supt∈[T ] ∥θ̂ (t) v ∥ ≤ B are all
bounded, the following inequality holds with probability at least 1− δ:
EL(θ̂(T );P) = Õ ( max { 1, θ∗,⊤Σxθ ∗σ2w, [V +B] 2 }√ log(4/δ)
n +
∥θ∗∥2 λT + Tλ[V +B]2 n
) ,
3In the rest of the paper, we use O or ≲ to omit the constant dependency and use Õ to omit the log dependency of n. For example, we write 2/n = O(1/n) or 2/n ≲ 1/n, log(n)/n = Õ(1/n), and log2(n)/n = Õ(1/n).
where the probability is taken over the randomness of training data.
By setting T = Θ( √ n) in Theorem 1, the upper bound is approximately Õ( 1√
n ), indicating that
the estimator at time T is consistent4. Below we show the comparison with some existing bounds, including benign overfitting bound and traditional stability-based bound5.
Comparison with benign overfitting. The bound in Bartlett et al. (2020) mainly focus on the minnorm solution (as T → ∞, the parameter trained in GD converges to the min-norm solution). Instead, the bound in Theorem 1 does not require the min-norm solution assumption.
Comparison with stability-based bound. The stability-based bound (without applying the decomposition framework) is Õ ( max{1, [V +B′]2} √ log(2/δ)
2n + Tλ n [V +B ′]2 ) , where B′ = supt ∥θ̂(t)∥
denotes the bound of the trained parameter in standard training. As a comparison, the notation B in Theorem 1 denote the bound of the trained parameter in variance training. One major difference between B and B′ is that: B is independent of ∥θ∗∥ while B′ is closely related to ∥θ∗∥. Therefore, with a large time T , the bound in Theorem 1 outperforms the stability-based bound when ∥θ∗∥ is large compared to V (namely, large signal-to-noise ratio)6.
Proof Sketch. We defer the proof details to Appendix C due to space limitations. Firstly, we provide Lemma 1 to show that the excess risk can be decomposed into VER and BER.
Lemma 1 Under the overparameterized linear regression settings with initialization θ(0) = θ(0)b = θ (0) v = 0, for any time T , the following decomposition inequality holds:
EL(θ̂(T );P) ≤ 2EvL(θ̂(T )v ;P) + 2EbL(θ̂ (T ) b ;P).
We next bound VER and BER separately in Lemma 2 and Lemma 3. When applying the stabilitybased bound on VER, we first convert VER into a generalization gap (see Equation 1) and then bound three terms individually.
Lemma 2 Under the overparameterized linear regression settings, assume that ∥x∥ ≤ 1, |ϵ| ≤ V , supt ∥θ̂ (t) v ∥ ≤ B are all bounded. Consider the gradient descent algorithm with constant stepsize λ, with probability at least 1− δ:
EvL ( θ̂(T )v ;P ) = O
( Tλ
n [V +B]
2 log(n) log(2n/δ) + max { 1, [V +B] 2 }√ log(2/δ)
2n
) . (2)
In Lemma 3, we apply uniform convergence to bound BER. Although we can bound it by directly calculating the excess risk, we still apply uniform convergence, because exact calculating is limited to a specific problem while uniform convergence is much more general. To link the uniform convergence and excess risk, we split the excess risk into two pieces inspired by Negrea et al. (2020). For the first piece, which stands for the generalization gap, we bound it by uniform convergence. For the second piece, which stands for the (surrogate) training loss, we derive a bound decreasing with time T .
Lemma 3 Under the overparameterized linear regression settings, assume that w = θ ∗,⊤x√
θ∗,⊤Σxθ∗ is
σ2w-subgaussian, then with probability at least 1− δ:
EbL(θ̂ (T ) b ;P) = O
( θ∗,⊤Σxθ
∗ √ n σ2w
√ log ( 2
δ
) +
1
λ[2T + 1] ∥θ∗∥2
) . (3)
4Consistency means the upper bound converges to zero as the number of training samples goes to infinity. 5We remark that the assumptions in Theorem 1 can be relaxed. Firstly, assuming x is subGaussian suffices to show that w is subGaussian. Secondly, assuming ϵ is subGaussian suffices to bound |ϵ| ≤ V with a log(n) cost.
6A large time t can guarantee that the dominating term in Theorem 1 is Tλ n [V +B′] 2. We remark that with
a large signal-to-noise ratio, we have Tλ n [V +B′] 2 ≥ Tλ n [V +B]2 ).
Although we do not present it explicitly, the results in Lemma 3 are derived under uniform convergence frameworks. Combining Lemma 1, Lemma 2 and Lemma 3 leads to Theorem 1. □

heading: 4 DECOMPOSITION UNDER GENERAL REGIMES
text: This section mainly discusses how to decompose the excess risk into variance component (VER) and bias component (BER) beyond linear regimes. Due to a fine-grained analysis on the signal and noise, the decomposition framework improves the stability-based bounds under the general regimes. We emphasize that the methods derived in this section are nevertheless the only way to reaching the decomposition goal and we leave more explorations for future work. To simplify the discussion, we assume the minimizer θ∗ of excess risk exists and is unique. Before diving into the main theorem, it is necessary to introduce the definition of local sharpness for the excess risk function, measuring the smoothness of the excess risk function around the minimizer θ∗.
Definition 2 (Local Sharpness) For a given excess risk function EL(θ;P) with EL(θ∗;P) = 0, we define s > 0 as the local sharpness factor when the following equation holds7:
0 < lim θ:∥θ−θ∗∥→0 EL(θ;P) ∥θ − θ∗∥s <∞.
Assumption 1 (Uniform Bound) Assume that for excess risk function EL(θ;P) with local sharpness s, given Px, for any Py|x and θ, the uniform bounds 0 < mu(Px) ≤Mu(Px) exist:
mu(Px) ≤ EL(θ;P) ∥θ − θ∗∥s ≤Mu(Px).
When the context is clear, we omit the dependency of Px and denote the bounds by mu and Mu.
Roughly speaking, the uniform bounds act as the maximal and minimal eigenvalues of the excess risk function EL(θ;P). For example, under the overparameterized linear regression regimes (Section 3), the uniform bounds are the maximum and minimum eigenvalues of the covariance matrix Σx.
Assumption 2 (Dynamics Decomposition Condition, DDC) We assume the trained parameters θ̂(T ) are (a,C,C ′)-bounded at time T , namely:
∥θ̂(T ) − θ∗∥ ≤ a(∥θ̂(T )v − θ∗v∥+ ∥θ̂ (T ) b − θ ∗ b∥) + C√ T + C ′√ n , (4)
where n denotes the size of training samples, T denotes the training time, and θ∗, θ∗v and θ ∗ b denote the minimizers of the corresponding population loss L(θ;P), L(θ;Pv), and L(θ;Pb), respectively.
We emphasize that although we use C√ T + C ′ √ n
in Assumption 2, it can be replaced with any o(1) term which converges to zero when T, n go to infinity. Empirically, several experiments demonstrate the generality of DDC and its variety (see Figure 3). From the theoretical perspective, we derive that overparameterized linear regression is (1, 0, 0)-bounded in DDC, and diagonal matrix recovery regimes satisfy DDC (see Section 4.1). We next show how to transfer DDC to the Decomposition Theorem in Theorem 2.
Theorem 2 (Decomposition Theorem) Assume that the excess risk has a unique minimizer with local sharpness factor s. Under Assumption 1 and Assumption 2, the following decomposition inequality holds:
EL(θ̂(T );P) ≤ [4a]s Mu mu
[ EvL(θ̂(T )v ;P) + EbL(θ̂ (T ) b ;P) ] +Mu [ 4C√ T ]s +Mu [ 4C ′√ n ]s . (5)
We defer the proofs to Appendix D due to limited space. In Theorem 2, Assumption 1 and Assumption 2 focus on the population loss and the training dynamics individually. We remark that our proposed generalization bound is algorithm-dependent due to parameters (a,C,C ′) in DDC.
7Here ∥ · ∥ refers to the ℓ2-norm.
By Theorem 2, we successfully decompose the excess risk dynamics into VER and BER. The remaining question is how to bound VER and BER individually. As we will derive in Section 4.1 via a case study on diagonal matrix recovery, we bound VER by stability techniques and BER by uniform convergence.

heading: 4.1 DIAGONAL MATRIX RECOVERY
text: From the previous discussion, we now conclude the basic procedure when applying the decomposition framework under Theorem 2. Generally, there are three questions to be answered in a specific problem: (a) Does the problem satisfy DDC? (b) How to bound VER using stability? (c) How to bound BER using uniform convergence? To better illustrate how to solve the three questions individually, we provide a case study on diagonal matrix recovery, a special and simplified case of a general low-rank matrix recovery problem but still retains the key properties. The diagonal matrix recovery regime and its variants are studied as a prototypical problem in Gissin et al. (2019); HaoChen et al. (2021).
Settings. Let X⋆ = Diag{σ⋆1 , · · · , σ⋆r , 0, · · · , 0} ∈ Rd×d denote a rank-r (low rank) diagonal matrix, where σ⋆1 ≥ · · · ≥ σ⋆r > 0. We define A ∈ Rd×d as the measurement matrix which is diagonal and y ∈ Rd as the measurement vector. Assume the ground truth yi = Aiσ⋆i + εi where yi, Ai denotes their corresponding i-th element, ϵi denotes the noise, and σ⋆i = 0 if r < i ≤ d. Our goal is to recover X⋆ from dataset D = {A(i), y(i)}i∈[n] with n training samples.
In the training process, we use the diagonal matrix U = Diag{u1, · · · , ud} to predict X⋆ via X̂ = UU⊤. Given the loss function8 Ln(U) = 1n ∑n j=1 ∥∥y(j) −A(j) ⊙ UU⊤∥∥2, we train the model with Gradient Flow9 U̇t = −∇Ln(Ut). During the analysis, we assume that Ai is O(1)uniformly bounded with unit variance, and εi is ν2-subGaussian with uniform bounds |ε(j)i | ≤ V .
8We use ⊙ to represent the Hadamard product. 9We use gradient flow instead of gradient descent mainly due to technical simplicity.
We emphasize that the training procedure of diagonal matrix recovery is fundamentally different from linear regimes since Ln(U) is non-convex with respect to U .
Theorem 3 (Decomposition Theorem in Diagonal Matrix Recovery) In the diagonal matrix recovery settings, if we set the initialization U0 = αI where α > 0 is the initialization scale, then for any 0 < t <∞, with probability at least 1− δ, we have:
EL(Ut;P) ≲ ( ∥X⋆∥2F + dV 2 + dα4 )√ log (d/δ)
n + r∑ j=1 σ4j σ2j + α 4eΩ(σjt) + α2d t
+ dV 2(t+ 1)
n log(n) log
( 2dn
δ
) + dα2 + log2 ( 1
α2
) rν2 log (r/δ)
n .
(6)
By setting α = Θ ( (d2n)− 1 4 ) and t = Θ(log (dnσr) /σr) in Theorem 3, the upper bound is
approximately Õ ( dV 2+∥X⋆∥2F√
n
) , indicating that the estimator at time t is consistent.
Comparison with stability-based bound. The stability-based bound (without applying the decomposition framework) is Õ ( dV 2+∥X⋆∥2F√
n + (dV 2+∥X⋆∥⋆)(t+1) n
) . As a comparison, the bound in The-
orem 3 outperforms the stability-based bound in the sense that we deduce the coefficient of the blowing-up term t+1n from dV
2 + ∥X⋆∥⋆ to dV 2, which is a significant improvement with a large signal-to-noise ratio under large time t. Besides, the bound in Theorem 3 captures the bias-variance tradeoff, meaning that the excess risk curve falls first and then rises.
Proof Sketch. We defer the proof details to Appendix E due to space limitations. Similar to Section 3, the proof of Theorem 3 consists of the following three lemmas, corresponding to the DDC, BER, and VER terms separately. We first prove that diagonal matrix recovery regimes satisfy DDCin Lemma 4.
Lemma 4 Under the assumptions in Theorem 3, for any 0 ≤ t <∞, with probability at least 1− δ, we have∥∥UtU⊤t −X⋆∥∥F ≤ ∥∥∥U bt U b,⊤t −X⋆∥∥∥F + ∥∥∥Uvt Uv,⊤t ∥∥∥F +O ( log ( 1 α2 )√ rν2 log (r/δ) n ) . (7) Different from the overparameterized linear regression regimes, Lemma 4 suffers from a O (
1√ n ) factor due to the non-linearity. This phenomenon can be verified empirically in Figure 3b and Figure 3e. We next apply stability techniques and uniform convergence to tackle VER and BER in Lemma 5 and Lemma 6, respectively.
Lemma 5 For VER, under the assumptions in Theorem 3, with probability at least 1− δ, we have
EvL(Ut;P) ≲ dV 2(t+ 1)
n log(n) log
( 2dn
δ
) + dα2 + dV 2 √ log (2d/δ)
n . (8)
Lemma 6 For BER, under the assumptions in Theorem 3, with probability at least 1− δ, we have
EbL(Ut;P) ≲ ( ∥X⋆∥2F + dα 4 )√ log (d/δ)
n + r∑ j=1 σ4j σ2j + α 4eΩ(σjt) + α2d t . (9)
Combining the above lemmas leads to Theorem 3. □
Remark: Directly applying Theorem 2 under overparameterized linear regression regimes costs a loss on condition number. The cost comes from the DDC condition where we transform the function space (excess risk function) to the parameter space (the distance between the trained parameter and the minimizer). We emphasize that the field of decomposition framework is indeed much larger than
Theorem 2, and there are other different approaches beyond Theorem 2. For example, informally, by introducing a pre-conditioner matrix A, changing DDC from the form ∥θ̂− θ⋆∥ to ∥Aθ̂−Aθ⋆∥ and setting A = Σ1/2x , one can avoid the cost of the condition number. We leave the technical discussion for future work.

heading: 5 CONCLUSION
text: In this paper, inspired by the distinct performances of neural networks on signal and noise, we improve the traditional stability-based bounds by proposing a decomposition framework. Under the framework, we decompose the excess risk into variance excess risk (VER) and bias excess risk (BER), representing the model ability to fit noise and signal, respectively. We then bound VER using stability-based techniques and bound BER using uniform convergence. Our bound captures the tendency of the excess risk and blows up much slower than the traditional stability-based bounds. We apply the decomposition framework into a linear case (overparameterized linear regression with gradient descent) and a non-linear case (matrix recovery with gradient flow), showing that the framework accords better with the theoretical and empirical evidence. We finally conduct several experiments on both synthetic datasets and real-world datasets to validate the utility of the decomposition framework on neural networks.

heading: Supplementary Materials
text: In this section, we first provide additional numerical experiment in Section A and figure illustration in Section B. We then give all the omitted proofs of lemmas and theorems in Section C, Section D, and Section E. We finally provide some further discussions, including why we apply stability on VER, the motivation behind DDC and a relaxed version of DDC.

heading: A NUMERICAL EXPERIMENT
text: In this section we provide additional numerical experiments on neural network to verify that our theory is satisfied for general neural network architectures and optimization algorithms.

heading: A.1 SYNTHETIC EXPERIMENT
text: In this section, our goal is to recover a sparse linear function f⋆(x) = ⟨θ⋆,x⟩ (x ∈ R30) with fully connected ReLU network. For each setting, we run 5 independent trials and report their average performances (the solid line) and the corresponding standard deviations (the error bar).
Effects of depth: We first explore the effect of depth and width for the excess risk decomposition property. In this experiment, we fix the width of each layer to be 64, and use SGD as the optimizer with Gaussian initialization (N (0, σ2) where σ = 1 × 10−3) and stepsize η = 1 × 10−2. For the dash line, we plot a(BER + VER) to verify the performance of our risk decomposition framework. Here a can be regarded as a measure of the success of our theory. Ideally a should be a small
constant so that our bound is non-vacuous. Specific to this experiment, for Figure 4a, 4b, 4c, we set a = 1.3, 2.0, 2.3, respectively. Overall, these results are well-predicted by our theory. Moreover, it is remarkable to notice that this constant becomes large once increasing the depth. This mainly contributes to the fact that deeper neural network has higher non-linearity.
Effects of width: In this experiment, we explore the effect of width for 2-layer ReLU network with SGD optimizer (the initialization and the stepsize are the same as above). For Figure 4d, 4e, 4e, we set a = 1.3, 1.0, 1.0, respectively. As can be seen, the wider the neural network, the better our risk decomposition works. We comment that this behavior also comes to the level of non-linearity. It is now well-known that when the width tends to infinity, the neural network converges to a linear function (e.g. (Jacot et al., 2018)).
Effects of optimizer: In this experiment, we try three different optimizers: SGD (the same hyperparameters as before), Adam (stepsize η = 0.002, (β1, β2) = (0.9, 0.999), ϵ = 1e− 08, no weight decay), Rprop (learning rate η = 5 × 10−4, (η1, η2) = (0.5, 1.2), stepsizes (1 × 10−6, 50)). For Figure 4g, 4h, 4i, we set a = 1.3, 1.2, 1.1. It can be seen that our risk decomposition framework is robust to different optimization algorithms, demonstrating its generality and potential to extend to general optimization problems.
A.2 CIFAR-10
In this section we deliver additional experiment in CIFAR-10 dataset. In this experiment, we use ResNet-18 and choose Adam as the optimizer. We regard that optimizing over the original dataset as the bias training (in the sense that all the label is noiseless, referring to the bias part). As for the variance training, we set the ground truth to be the uniform logit v = [0.1, · · · , 0.1]⊤ (here we use one-hot coding and choose the categorical cross entropy loss). As for the labels in the training dataset, we generate it uniformly from {e1, · · · , e10}, where ei is the standard basis in R10. Under such data generating model, its Bayes risk is provided by −10 × 0.1 × log(0.1) ≈ 2.3026 for categorical cross entropy loss. For the standard training, we generate the training label as follows: y = (1 − p)y⋆ + pej . Here p ∈ [0, 1] is the corruption probability, y⋆ is the true label and j is uniformly sampled from {1, · · · , 10}. In this experiment, we set the corruption probability to be 0.4. The numerical result is provided in Figure 5 . Different from the performance on MINIST, here the variance training actually converges to the ground truth. We conjecture that this is due to the structure of ResNet, resulting to some kind of implicit regularization towards the training process.10

heading: B FIGURE ILLUSTRATION

heading: B.1 FIGURE 1A
text: Figure 1a conveys an important message that variance training has a much smaller gradient norm at the beginning compared to standard training. This mainly comes from the difference in their landscapes. Since there is no signal term in the variance learning, we can regard the landscape corresponding to the variance training as just a translation transformation of the original landscape (Figure 1a, the blue line is just a translation transformation of the green line). Based on this difference, they pursue different optima. We denote θ̂v, θ̂ the empirical loss minimizer, respectively. Due to concentration, we have θ̂v ≈ 0, and θ̂ ≈ θ⋆. Note that we use a near-zero initialization θ(0) = θ (0) v ≈ 0, we immediately have that θ(0)v ≈ θ̂v , indicating that the gradient norm is small. On the other hand, since θ(0) is far away from the optimum, its initial gradient norm can be very large.
To better illustrate this principle, we use linear regression to provide a simple example. For the standard training, the empirical loss is Ln(θ;P) = 1n ∑n i=1 (⟨θ,Xi⟩ − ⟨θ⋆, Xi⟩ − εi) 2. Hence, we
10We do observe that for certain learning rate regimes, the algorithm will overfit the training dataset for variance training, leading to a large VER (∼ 10).
have
∇Ln(θ(0);P) = 2
n n∑ i=1 (〈 θ(0), Xi 〉 − ⟨θ⋆, Xi⟩ − εi ) Xi
= − 2 n n∑ i=1 (⟨θ⋆, Xi⟩+ εi)Xi ≈ −2θ⋆, (10)
if we assume Xi ∼ N (0, Id×d). On the other hand, for the variance training, the empirical loss is Ln(θ;Pv) = 1n ∑n i=1 (⟨θ,Xi⟩ − εi) 2. Therefore, we have
∇Ln(θ(0)b ;Pv) = 2
n n∑ i=1 (〈 θ (0) b , Xi 〉 − εi ) Xi = − 2 n n∑ i=1 εiXi ≈ 0. (11)
In conclusion, we have∥∥∥∇Ln(θ(0);P)∥∥∥ ≈ 2 ∥θ⋆∥ , ∥∥∥∇Ln(θ(0)v ;Pv)∥∥∥ ≈ 0. (12) Obviously, there is a huge initial gradient norm gap between standard and variance training, showing that our bound is significantly better than the direct stability-based bound under such regimes.

heading: B.2 FIGURE 1B
text: Figure 1b provides an example showing that the trajectory experienced by the variance training is likely to stay a convex region, though the global landscape can be highly non-smooth and nonconvex. The reason is similar to that of Figure 1a. Since θ(0)v ≈ θ̂v , the gradient norm is small and the training dynamics may stay close to the optimum, which lies on a smooth and convex region. On the contrary, the standard training needs to cross potentially many bad local minimum or saddle points, suffering from large cumulative gradient norms.

heading: B.3 FIGURE 2
text: The setting of Figure 2 is similar to that of Section A. For MINIST dataset, we use three layer convoluted neural network with Adam optimizer.

heading: B.4 FIGURE 3
text: Figure 3 provides the numerical verification of Dynamics Decomposition Condition, and our upper bounds for both linear regression and matrix recovery problems 11.
Linear Regression We set the dimension d = 500, the sample sizes are n = 300, 800, respectively. Samples X1 · · · , Xn are i.i.d. from standard Gaussian distribution N (0, Id×d), and the output yi = ⟨Xi, θ⋆⟩+εi where εi ∼ N (0, σ2) and σ = 2. We set the initialization θ(0) = θ(0)b = θ (0) v = 0 and the learning rate λ = 0.01.
Matrix Recovery We set the dimension d = 20, rank r = 3, the ground truth X⋆ = V ΣV ⊤, where Σ = Diag{5, 3, 1, 0, · · · , 0} and V is an orthonormal matrix. The numbers of the measurement matrices are 200, 600, respectively. The measurement matrices A1, · · · , An are i.i.d. standard Gaussian matrices. The corresponding measurement yi = ⟨Ai, X⋆⟩ + εi, where εi ∼ N (0, 1). In this numerical experiment, we use the initialization U(0) = Ub(0) = Uv(0) = αId×d where α = 0.01. We solve this problem via gradient descent with constant stepsize η = 0.1.
For the upper bound, it is calculated via upper bound = ∥∥Xbt −X⋆∥∥F + ∥Xvt ∥F + 8.5√n , where the last term is corresponding to the additional term in Lemma 4.
Neural Network We the dimension d = 30, the ground truth function is f⋆(x) = ⟨θ⋆,x⟩ where θ⋆ is sparse and ∥θ⋆∥ = Θ(1). The covariate x is generated from Gaussian distribution. The label for bias training is exactly yi = f⋆(xi). For the variance trainng, its label is generated from an additional Gaussian distribution εi ∼ N (0, σ2) where σ = 1.5. For the standard training, its label is yi = f ⋆(xi)+εi. For the upper bound, it is calculated via upper bound = ∥∥∥fθbt − f⋆∥∥∥P +∥∥fθvt ∥∥P . Here ∥·∥P is defined to be the L2(P) norm, i.e., ∥g∥ 2 P = Ex P [ f(x)2 ] . For these settings, we choose SGD (with stepsize η = 1 × 10−3) as the optimizer where the initialization is randomly generated from a Gaussian distribution.

heading: C PROOF FOR OVERPARAMETERIZED LINEAR REGRESSION

heading: C.1 PROOF OF THEOREM 1
text: In this section, we prove Theorem 1 stated as follows:
Theorem 1 (Linear Regression Regimes) Under the overparameterized linear regression settings, assume that w = θ
∗,⊤x√ θ∗,⊤Σxθ∗ is σ2w-subgaussian. If ∥x∥ ≤ 1, |ϵ| ≤ V , supt∈[T ] ∥θ̂ (t) v ∥ ≤ B are all
bounded, the following inequality holds with probability at least 1− δ:
EL(θ̂(T );P) = Õ ( max { 1, θ∗,⊤Σxθ ∗σ2w, [V +B] 2 }√ log(4/δ)
n +
∥θ∗∥2 λT + Tλ[V +B]2 n
) ,
where the probability is taken over the randomness of training data.

heading: Proof:
text: 11Here, we extend the diagonal matrix recovery problem to the general matrix recovery problem.
We first calculate the population loss under linear regression regime, namely, for any time t: L(θ̂(t);P) = Ex,y [ y − x⊤θ̂(t) ]2 (a) = E [ x⊤θ∗ + ϵ− x⊤θ̂(t)
]2 = E [ x⊤θ∗ − x⊤θ̂(t) ]2 + E [ ϵx⊤θ̂(t) ] + Eϵ2
= [ θ∗ − θ̂(t) ]⊤ Σx [ θ∗ − θ̂(t) ] + ExEϵ [ ϵx⊤θ̂(t)|x ] + Eϵ2
= [ θ∗ − θ̂(t) ]⊤ Σx [ θ∗ − θ̂(t) ] + Exx⊤θ̂(t)Eϵ [ϵ|x] + Eϵ2
(b) = [ θ∗ − θ̂(t) ]⊤ Σx [ θ∗ − θ̂(t) ] + Eϵ2,
where (a) we use the condition that y = x⊤θ∗ + ϵ and (b) we use the assumption that E[ϵ|x] = 0. Therefore, the excess risk can be written as
E(θ̂(t)) = [ θ∗ − θ̂(t) ]⊤ Σx [ θ∗ − θ̂(t) ] .
Now plugging the Lemma 1, Lemma 2, Lemma 3, into Theorem 2, we derive that with probability at least 1− δ:
EL(θ̂(T );P) ≲ EvL(θ̂(T )v ;P) + EbL(θ̂ (T ) b ;P)
≲ θ ∗,⊤Σxθ ∗ √ n σ2w
√ log ( 4
δ
) +
1
λ [2T + 1] ∥θ∗∥2
+ Tλ
n [V +B]
2 log(n) log(2n/δ) + max { 1, [V +B] 2 }√ log(4/δ)
2n .
= Õ ( max { 1, θ∗,⊤Σxθ ∗σ2w, [V +B] 2 }√ log(4/δ)
n +
∥θ∗∥2
λT + Tλ [V +B]
2
n
) ,
which completes the proof. □

heading: C.2 PROOF OF LEMMA 1
text: We first recall Lemma 1 as follows:
Lemma 1 Under the overparameterized linear regression settings with initialization θ(0) = θ(0)b = θ (0) v = 0, for any time T , the following decomposition inequality holds:
EL(θ̂(T );P) ≤ 2EvL(θ̂(T )v ;P) + 2EbL(θ̂ (T ) b ;P).
Proof: Firstly, when we use gradient descent, we have that:
θ̂(t+1) = θ̂(t) + λ
n ∑ i xi[yi − x⊤i θ̂(t)],
θ̂(t+1)v = θ̂ (t) v +
λ
n ∑ i xi[ϵi − x⊤i θ̂(t)v )],
θ̂ (t+1) b = θ̂ (t) b +
λ
n ∑ i xi[x ⊤ i θ ∗ − x⊤i θ̂ (t) b ].
We next prove the conclusion by induction. Since we choose θ(0) = θ(0)v = θ (0) b = 0, the conclusion holds at the initialization θ(0) = θ(0)v + θ (0) b .
Now assume that θ̂(t) = θ̂(t)v + θ̂ (t) b holds at time t, we have
θ̂(t+1) = θ̂(t) + λ
n ∑ i xi[yi − x⊤i θ̂(t)]
= θ̂(t)v + θ̂ (t) b +
λ
n ∑ i xi[ϵi + x ⊤ i θ ∗ − x⊤i θ̂(t)v + θ̂ (t) b ]
= [ θ̂(t)v + λ
n ∑ i xi[ϵi − x⊤i θ̂(t)v ]
] + [ θ̂ (t) b + λ
n ∑ i xi[x ⊤ i θ ∗ − x⊤i θ̂ (t) b ] ] = θ̂(t+1)v + θ̂ (t+1) b .
In summary, we have: θ̂(t) = θ̂(t)v + θ̂ (t) b holds for any time t. Besides, we easily calculate that θ∗ = θ∗v + θ ∗ b since θ ∗ = θ∗b = θ ∗ and θ∗v = 0.
Therefore, we have θ̂(t) − θ⋆ = θ̂(t)v − θ⋆v + θ̂ (t) b − θ⋆b holds for any t. By CauchySchwarz inequality, we derive that for a specific time T :
EL(θ̂(T );P) = [ θ∗ − θ̂(T ) ]⊤ Σx [ θ∗ − θ̂(T ) ] = ∥∥∥θ∗ − θ̂(T )∥∥∥2
Σx = ∥∥∥θ∗v − θ̂(T )v + θ∗b − θ̂(T )b ∥∥∥2
Σx
(a) ≤2 ∥∥∥θ∗v − θ̂(T )v ∥∥∥2
Σx + 2 ∥∥∥θ∗b − θ̂(T )b ∥∥∥2 Σx
≤2EvL(θ̂(T )v ;P) + 2EbL(θ̂ (T ) b ;P).
where we denote the norm ∥u∥2A = u⊤Au. □

heading: C.3 PROOF OF LEMMA 2
text: Lemma 2 Under the overparameterized linear regression settings, assume that ∥x∥ ≤ 1, |ϵ| ≤ V , supt ∥θ̂ (t) v ∥ ≤ B are all bounded. Consider the gradient descent algorithm with constant stepsize λ, with probability at least 1− δ:
EvL ( θ̂(T )v ;P ) = O
( Tλ
n [V +B]
2 log(n) log(2n/δ) + max { 1, [V +B] 2 }√ log(2/δ)
2n
) . (2)
Proof: When the context is clear, we omit v in this section and denote the trained parameter by θ̂(t). For any time t, we first split the excess risk into three components following Equation 1.
EL(θ̂(t);Pv) = [ L(θ̂(t);Pv)− L(θ̂(t);Pnv ) ] ︸ ︷︷ ︸
(I)
+ [ L(θ̂(t);Pnv )− L(θ∗;Pnv ) ] ︸ ︷︷ ︸
(II)
+ [L(θ∗;Pnv )− L(θ∗;Pv)]︸ ︷︷ ︸ (III) .
(13)
Part (I). We bound the first component in Equation 13 using stability.
Fact A: the loss ℓ(θ;x,y) is 2 [V +B]-Lipschitz.
∥∇θℓ(θ;x,y)∥ = 2∥x(ϵ− x⊤θ)∥ ≤ 2 [V +B] .
Fact B: the stability of the parameter is 2λtn [V +B]. Denote the parameter trained on dataset S and the parameter trained on dataset S′ as θ̂(t)S and θ̂ (t) S′ , respectively. We will show that, when S and S ′
have only one different sample (e.g., WLOG S = {x1,x2, · · · ,xn} and S′ = {x′1,x2, · · · ,xn}), θ̂ (t) S is close to θ̂
(t) S′ .∥∥∥θ̂(t+1)S − θ̂(t+1)S′ ∥∥∥
≤ ∥∥∥∥[I − λnX⊤X ] θ̂ (t) S + λ n X⊤Y − [ I − λ n X ′,⊤X ′ ] θ̂ (t) S′ − λ n X ′,⊤Y ′ ∥∥∥∥ ≤ ∥∥∥∥[I − λnX⊤X ] [ θ̂ (t) S − θ̂ (t) S′ ] ∥+ ∥ [ λ n X ′,⊤X ′ − λ n X⊤X ] θ̂ (t) S′ ∥+ λ n ∥X⊤Y −X ′,⊤Y ′
∥∥∥∥ ≤ ∥∥∥[θ̂(t)S − θ̂(t)S′ ]∥∥∥+ λn ∥∥∥[x1x⊤1 − x′1x⊤,′1 ] θ̂(t)S′ ∥∥∥+ λn ∥∥∥x⊤1 ϵ1 − x′,⊤1 ϵ′1∥∥∥
≤ ∥∥∥[θ̂(t)S − θ̂(t)S′ ]∥∥∥+ 2λn [V +B] .
Summation over time t provides the following result:
∥θ̂(t)S − θ̂ (t) S′ ∥ ≤
2tλ
n [V +B] .
Therefore, the total stability with respect to the ℓ2-loss is
|ℓ(θ̂(t)S )− ℓ(θ̂ (t) S′ )| ≤
4tλ
n [V +B]
2 .
By Proposition 1, we have that
L(θ̂(t);Pv)− L(θ̂(t);Pnv ) ≲ tλ
n [V +B]
2 log(n) log(2n/δ) +
√ log(1/δ)
n .
Part II. We next show that the second component in Equation 13 is less than zero based on Lemma 7.
L(θ(t);Pnv )− L(θ∗;Pnv )
= 1
n
[ ∥ϵ−Xθ(t)∥2 − ∥ϵ∥2 ] = 1
n ϵ⊤
[ −2X [ I − [ I − λ
n X⊤X
]t] X† +X†,⊤ [ I − [ I − λ
n X⊤X
]t] X⊤X [ I − [ I − λ
n X⊤X
]t] X† ] ϵ.
The i-th eigenvalue of the matrix [−2X[I− [I− λnX ⊤X]t]X†+X†,⊤[I− [I− λnX ⊤X]t]X⊤X[I− [I − λnX ⊤X]t]X†] can be calculated as follows
−2 [ 1− (1− λσi)t ] + [ 1− (1− λσi)t ]2 = (1− λσi)2t − 1 ≤ 0,
where σi refers to the i-th eigenvalue of the matrix 1nX ⊤X . Therefore, we have
L(θ(t);Pnv )− L(θ∗;Pnv ) ≤ 0.
Part III. We bound the third component in Equation 13 using concentration.
Note that |ℓ2(θ;x,y)| = |(y − x⊤θ)2| ≤ [V +B]2, so we can use Hoeffding’s inequality, P [|L(θ∗;Pnv )− L(θ∗;Pv)| ≥ ϵ] ≤ 2 exp [ − 2nϵ 2
(B + V )4
] .
By setting ϵ = [V +B]2 √
log(2/δ) 2n , then with probability at least 1− δ, we have
|L(θ∗;Pnv )− L(θ∗;Pv)| ≤ [V +B] 2
√ log(2/δ)
2n .
Combining Part (I), (II), and (III) together, at time T , we have:
EL(θ̂(T );Pv) ≲ Tλ
n [V +B]
2 log(n) log(2n/δ) +
√ log(1/δ)
n + [V +B]
2
√ log(2/δ)
2n
≲ Tλ n [V +B] 2 log(n) log(2n/δ) + max{1, [V +B]2}
√ log(2/δ)
2n .
□

heading: C.4 PROOF OF LEMMA 3
text: Lemma 3 Under the overparameterized linear regression settings, assume that w = θ ∗,⊤x√
θ∗,⊤Σxθ∗ is
σ2w-subgaussian, then with probability at least 1− δ:
EbL(θ̂ (T ) b ;P) = O
( θ∗,⊤Σxθ
∗ √ n σ2w
√ log ( 2
δ
) +
1
λ[2T + 1] ∥θ∗∥2
) . (3)
Proof: Since we use the noiseless data (x,E(y|x)) and E(y|x) = x⊤θ∗, by Lemma 7 we have
θ (t) b =
[ I − [ I − λ
n X⊤X
]t] X†Xθ∗ = θ∗ [ I − [ I − λ
n X⊤X
]t] θ∗,
where we use the fact that [ I − [ I − λnX ⊤X ]t] X†X = I − [ I − λnX ⊤X ]t .
Therefore, BER can be expressed as:
EL(θ;Pnb ) = [ θ∗ − θ(t)b ]⊤ Σx [ θ∗ − θ(t)b ] = θ∗,⊤ [ I − λ
n X⊤X
]t Σx [ I − λ
n X⊤X
]t θ∗
= θ∗,⊤ [ I − λ
n X⊤X
]t [ Σx − 1
n X⊤X
] [ I − λ
n X⊤X
]t θ∗
+ θ∗,⊤ [ I − λ
n X⊤X
]t [ 1
n X⊤X
] [ I − λ
n X⊤X
]t θ∗.
We split the excess risk into two parts. Intuitively, it is just like the decomposition in Equation 1. The first part is similar to the generalization gap to use uniform convergence to bound it. Note that here we require that λ < supX0 1
∥ 1nX0⊤X0∥ . Hence, we have
∣∣∣∣∣θ∗,⊤ [ I − λ n X⊤X ]t [ Σx − 1 n X⊤X ] [ I − λ n X⊤X ]t θ∗ ∣∣∣∣∣ ≤
∣∣∣∣∣supX0 θ∗,⊤ [ I − λ n X0 ⊤X0 ]t [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗ ∣∣∣∣∣ =
∣∣∣∣∣supX0 tr ( θ∗,⊤ [ I − λ n X0 ⊤X0 ]t [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗ )∣∣∣∣∣ =
∣∣∣∣∣supX0 tr ([ I − λ n X0 ⊤X0 ]t [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗θ∗,⊤ )∣∣∣∣∣ =
∣∣∣∣∣supX0 ∥∥∥∥∥ [ I − λ n X0 ⊤X0 ]t [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗θ∗,⊤ ∥∥∥∥∥ ∣∣∣∣∣
≤ ∣∣∣∣∣supX0 ∥∥∥∥∥ [ I − λ n X0 ⊤X0 ]t∥∥∥∥∥ ∥∥∥∥∥ [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗θ∗,⊤ ∥∥∥∥∥ ∣∣∣∣∣
≤ ∣∣∣∣∣supX0 ∥∥∥∥∥ [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗θ∗,⊤ ∥∥∥∥∥ ∣∣∣∣∣
= ∣∣∣∣∣supX0 tr ( θ∗,⊤ [ Σx − 1 n X⊤X ] [ I − λ n X0 ⊤X0 ]t θ∗ )∣∣∣∣∣ ≤ ∣∣∣∣sup X0 tr ( θ∗,⊤ [ Σx − 1 n X⊤X ] θ∗ )∣∣∣∣
= ∣∣∣∣θ∗,⊤ [Σx − 1nX⊤X ] θ∗ ∣∣∣∣ ,
where we apply Lemma 8 repeatedly since rank [ θ∗,⊤θ∗A ] = 1 whereA stands for arbitrary matrix.
We highlight the difference between X0 and X . The notation X represents the contribution of the training set to the training error. The notation X0 represents the contribution of the training set to the training estimator θ̂. The uniform convergence is taken over the estimator θ. Therefore, we take sup operator on X0.
We next bound θ∗,⊤ [ Σx − 1nX ⊤X ] θ∗. Assume that w = θ∗,⊤x/ √ θ∗,⊤Σxθ∗ is SubGaussian
with subGaussian norm σw. Obviously, E [w] = 0,E [ w2 ] = 1. Therefore, we have that: w2 − 1 is
sub-Exponential distributions with sub-Exponential norm ∥∥w2 − 1∥∥
ψ1 ≤ σ2w. While ϵ/ [ θ⊤Σxθ ] < σ2w, we have the following tail bound
P ( |θ∗,⊤ [ Σx − 1
n X⊤X
] θ∗| ≥ ϵ ) = P ( | 1 n n∑ i=1 [ (x⊤i θ) 2 − θ⊤Σxθ ] | ≥ ϵ )
= P ( | 1 n n∑ i=1 [ w2i − 1 ] | ≥ ϵ/ [ θ⊤Σxθ ])
≤ 2 exp [ −c [ ϵ/ [ θ⊤Σxθ ]]2 nσ4w ] ,
where c is a universal constant. Therefore, by setting ϵ = [ θ⊤Σxθ ] σ2w log[2/δ]√ c √ n with probability at least 1− δ, we have
∣∣∣∣θ∗,⊤ [Σx − 1nX⊤X ] θ∗ ∣∣∣∣ ≤ [θ⊤Σxθ] σ2w √ log [2/δ]√ c √ n .
With abuse to use c as a constant, we have with probability at least ≥ 1− δ∣∣∣∣θ∗,⊤ [Σx − 1nX⊤X ] θ∗ ∣∣∣∣ ≤ c [θ⊤Σxθ] σ2w √ log(2/δ)√ n .
For the second part, which is closely related to the empirical loss of θ(t)b , We will show that: for a general case, the empirical loss converges with rate O(1/t). Denote the eigenvalues of 1nX
⊤X as σ1, · · · , σn, 0, · · · , 0. Then the corresponding eigenvalues of matrix[ I − λnX ⊤X ]t [ 1 nX ⊤X ] [ I − λnX ⊤X ]t are (1− λσi)2tσi. Therefore, we have:
θ∗,⊤ [ I − λ
n X⊤X
]t [ 1
n X⊤X
] [ I − λ
n X⊤X
]t θ∗
≤∥θ∗∥2 max i (1− λσi)2tσi
≤ 1 λ(2t+ 1) ∥θ∗∥2 .
The maximum is attained while σi = 1λ(2t+1) .
Combining these terms leads to the conclusion. □

heading: C.5 AUXILIARY LEMMAS
text: Lemma 7 The dynamics of θ̂(t) using GD is:
θ̂(t) = [ I − λ
n X⊤X
]t [ θ(0) −X†Y ] +X†Y.
where X† represents the Pseudo inverse of matrix X .
Proof: The proof is based on induction. We notice that when t = 0, the equation directly follows. Assume when t = t, the equation holds. Then at time t+ 1, we have
θ(t+1) = θ̂(t) + λ
n X⊤(Y −Xθ̂(t))
= [ 1− λ
n X⊤X
] [[ I − λ
n X⊤X
]t [ θ0 −X†Y ] +X†Y ] + λ
n X⊤Y
= [ I − λ
n X⊤X
]t+1 [ θ0 −X†Y ] + [ 1− λ
n X⊤X
] X†Y + λ
n X⊤Y
= [ I − λ
n X⊤X
]t+1 [ θ0 −X†Y ] +X†Y,
where we use the fact that X⊤XX† = X⊤.
□
Lemma 8 For rank-1 matrix A, we have that ∥A∥ = |tr [A]|.

heading: D PROOF OF THEOREM 2
text: In this section, we provide the proof of Theorem 2. We first recall the formal statement of Theorem 2.
Theorem 2 (Decomposition Theorem) Assume that the excess risk has a unique minimizer with local sharpness factor s. Under Assumption 1 and Assumption 2, the following decomposition inequality holds:
EL(θ̂(T );P) ≤ [4a]s Mu mu
[ EvL(θ̂(T )v ;P) + EbL(θ̂ (T ) b ;P) ] +Mu [ 4C√ T ]s +Mu [ 4C ′√ n ]s . (5)
Proof: Under Assumption 1 (upper bound), we have:
EL(θ̂(t);P) = EL(θ̂(t);P) ∥θ̂(t) − θ̂∗∥s ∥θ̂(t) − θ̂∗∥s
≤ [ M1/su ∥θ̂(t) − θ̂∗∥ ]s .
Combining Assumption 2 (DDC) and the lower bound in Assumption 1, we have
∥θ̂(t) − θ̂∗∥
≤a∥θ̂(t)b − θ̂ ∗ b∥+ a∥θ̂(t)v − θ̂∗v∥+ C√ t + C ′√ n
≤[ 1 mu ]1/saE1/sL (θ̂ (t) b ;Pb) + [ 1 mu ]1/saE1/sL (θ̂ (t) v ;Pv) + C√ t + C ′√ n .
Therefore, we have that
EL(θ̂(t);P) ≤ [ Mu mu 1/s aE1/sL (θ̂ (t) b ;Pb) + Mu mu 1/s aE1/sL (θ̂ (t) v ;Pv) + C√ t + C ′√ n ]s ≤ 4s
[ Mu mu asEL(θ̂(t)b ;Pb) + Mu mu asEL(θ̂(t)v ;Pv) + [ M1/su C√ t ]s + [ M1/su C ′√ n ]s] ≤ [4a]s Mu
mu
[ EL(θ̂(t)b ;Pb) + EL(θ̂ (t) v ;Pv) ] +Mu [ 4C√ t ]s +Mu [ 4C ′√ n ]s .
where we use the fact that (a+ b+ c+ d)p ≤ max{(4a)p, (4b)p, (4c)p, (4d)p} ≤ (4a)p + (4b)p + (4c)p + (4d)p when a, b, c, d, p > 0.
□

heading: E PROOF FOR DIAGONAL MATRIX RECOVERY
text: In this section, we prove Theorem 3 via proving the proceeding three lemmas separately.

heading: E.1 PROOF OF THEOREM 3
text: Theorem 3 (Decomposition Theorem in Diagonal Matrix Recovery) In the diagonal matrix recovery settings, if we set the initialization U0 = αI where α > 0 is the initialization scale, then for any 0 < t <∞, with probability at least 1− δ, we have:
EL(Ut;P) ≲ ( ∥X⋆∥2F + dV 2 + dα4 )√ log (d/δ)
n + r∑ j=1 σ4j σ2j + α 4eΩ(σjt) + α2d t
+ dV 2(t+ 1)
n log(n) log
( 2dn
δ
) + dα2 + log2 ( 1
α2
) rν2 log (r/δ)
n .
(6)
Note that although we do optimization based on U , we consider its square X = UU⊤ as the parameter since we require the solution is unique during the analysis. Note that we do not change the update rule (where GD is still taken on U ). One can directly check that under Diagonal Matrix Recovery regimes, s = 2 and Mu = mu = 1. Therefore, combining Lemma 4, Lemma 5 and Lemma 6 leads to Theorem 3.

heading: E.2 PROOF OF LEMMA 4
text: Lemma 4 Under the assumptions in Theorem 3, for any 0 ≤ t <∞, with probability at least 1− δ, we have∥∥UtU⊤t −X⋆∥∥F ≤ ∥∥∥U bt U b,⊤t −X⋆∥∥∥F + ∥∥∥Uvt Uv,⊤t ∥∥∥F +O ( log ( 1 α2 )√ rν2 log (r/δ) n ) . (7)
Proof: Since all the measurement matrices A(1), · · · , A(n) are diagonal, we can analyze the dynamic of each singular value separately.
For a given singular value σ, assume ai ∼ N (0, 1), εi ∼ subG(ν2). Then we define σ1 = 1 n ∑n i=1 a 2 iσ, and σ2 = 1 n ∑n i=1 aiεi, denote σ = σ1 + σ2, and ξ = 1 n ∑n i=1 a 2 i . Denote u
2(t), u2b(t), u 2 v(t) the dynamics of standard training, bias training, and variance training, respectively. Then our goal is to show that |σ − u2(t)| ≤ |σ − u2b(t)| + u2v(t) + small term holds with high probability.
Before diving into the details, we first give the dynamics of u2(t), u2b(t), u 2 v(t). If σ > 0, we have
u2(t) = σα2e2σt
σ − ξα2 + ξα2e2σt .
u2b(t) = σ1α
2e2σ1t
σ1 − ξα2 + ξα2e2σ1t .
Furthermore, if σ2 > 0, we have
u2v(t) = σ2α
2e2σ2t
σ2 − ξα2 + ξα2e2σ2t .
Otherwise, if σ2 < 0 we have
u2v(t) = |σ2|α2
(|σ2|+ ξα2) e2|σ2|t − ξα2 .
If σ⋆ = 0, we have
u2b(t) = α2
ξα2t+ 1 .
Now we turn to the proof of Lemma 4. If σ = 0, then σ = σ2, which implies u2(t) = u2v(t), indicating that the conclusion holds.
Hence, we only consider the case that σ ≥ σr > 0. By triangle inequality, it suffices to show∣∣u2(t)− u2b(t)∣∣ = O (|σ2|) . Without loss of generality, we assume σ > σ1 > 0. Note that∣∣u2(t)− u2b(t)∣∣ = ∣∣∣∣ σα2e2σtσ + ξα2e2σt − σ1α2e2σ1tσ1 + ξα2e2σ1t
∣∣∣∣+O(|σ2|) = α2σσ1 ( e2σt − e2σ1t ) (σ + ξα2e2σt) (σ1 + ξα2e2σ1t) +O(|σ2|).
Via integration-by-parts formula, we have e2σt − e2σ1t = 2σ2e2σt ∫ t 0 e−2σ2sds ≤ 2σ2te2σt.
Therefore, it suffices to show that
σσ1α 2te2σt
(σ + ξα2e2σt) (σ1 + ξα2e2σ1t) = Õ(1).
We prove that there exists a universal constant C, such that
σσ1α 2te2σt
(σ + ξα2e2σt) (σ1 + ξα2e2σ1t) ≤ C,
which is equivalent to show ϕ(t) = C ( σ + ξα2e2σt ) ( σ1 + ξα 2e2σ1t ) − σσ1α2te2σt ≥ 0, ∀0 ≤ t <∞.
Note that ϕ(0) = C ( σ + ξα2 ) ( σ1 + ξα 2 ) ≥ 0, and its gradient is
ϕ′(t) = 2Cξ2α4(σ + σ1)e 2(σ+σ1)t + 2Cσσ1ξα 2 ( e2σt + e2σ1t ) − σσ1α2e2σt − 2σ2σα2te2σt.
By concentration of Chi-square distribution, we have ξ ≥ 1− √
log(1/δ) n = Ω(1) with probability at least 1− δ, providing that n ≳ log ( 1 δ ) . It suffices to show that ϕ′(t) ≥ 0, which is equivalent to
ψ(t) = Cα2e2σ1t + Cσ1 − σσ1t ≥ 0. It attains its minimum σ2 + Cσ1 − σ 2 log ( 1 2Cα2 ) at the point t = log ( 1 2Cα2 ) /2σ1. Therefore, it
suffices to set C = log (
1 α2
) .
Overall, we have ∣∣u2(t)− u2b(t)∣∣ = O(log 1α2 |σ2| )
if σ⋆ > 0;∣∣u2(t)− u2b(t)∣∣ ≤ u2v(t) if σ⋆ = 0. Hence, with proper choice of probability, we have with probability at least 1−δ, the following holds
∥Xt −X⋆∥F ≤ ∥∥Xbt −X⋆∥∥F + ∥Xvt ∥F +O ( log ( 1
α2
)√ rσ2 log (r/δ)
n
) .
To summarize, the case is ( a = 1, C = 0, C ′ = O ( log (
1 α2
)√ rσ2 log (r/δ) )) -bounded (See
DDC).
□

heading: E.3 PROOF OF LEMMA 5
text: Lemma 5 For VER, under the assumptions in Theorem 3, with probability at least 1− δ, we have
EvL(Ut;P) ≲ dV 2(t+ 1)
n log(n) log
( 2dn
δ
) + dα2 + dV 2 √ log (2d/δ)
n . (8)

heading: Proof:
text: Similar to the proof in overparameterized linear regression regimes, we first split the excess risk in three componenets:
EL(Ut;Pv) = L ( u2(t) ) − Ln ( u2(t) )︸ ︷︷ ︸ (I) +Ln ( u2(t) ) − Ln (0)︸ ︷︷ ︸ (I) +Ln (0)− L (0)︸ ︷︷ ︸ (III) . (14)
Part (I). We bound the first component in Equation 13 using stability. We define σ2 = 1n ∑n i=1 aiεi, σ′2 = 1 n ∑n−1 i=1 aiεi + a ′ nε ′ n. Similarly, we define ξ = 1 n ∑n i=1 a 2 i , and ξ ′ = 1n ∑n−1 i=1 a 2 i + a ′2 n . In the first case, we assume both σ2, σ′2 > 0, so that their difference can be bounded
|u2v(t)− u′2v (t)| ≲ ∣∣∣α2σ2σ′2 (e2σ2t − e2σ′2t)∣∣∣+ α4 ∣∣∣ξ′σ2e2σ2t − ξσ′2e2σ′2t∣∣∣+ α4(σ2ξ′ − σ′2ξ)e2(σ2+σ′2)t (σ2 − ξα2 + ξα2e2σ2t) ( σ′2 − ξ′α2 + ξ′α2e2σ ′ 2t ) . WLOG, we assume σ2 > σ′2. Then, similar to the proof in Lemma 4, we have∣∣∣α2σ2σ′2 (e2σ2t − e2σ′2t)∣∣∣ (σ2 − ξα2 + ξα2e2σ2t) ( σ′2 − ξ′α2 + ξ′α2e2σ ′ 2t ) = O (|σ2 − σ′2|t) .
The second term is dominated by the first term. Hence, we only need to handle the third term. Note that
|σ2ξ′ − σ′2ξ| = |σ2 (ξ′ − ξ) + ξ(σ2 − σ′2)| ≤ |σ2| n ∣∣a2n − a,2n ∣∣+ ∣∣∣∣ ξn (anεn − a′nε′n) ∣∣∣∣ ,
we have α4(σ2ξ ′ − σ′2ξ)e2(σ2+σ ′ 2)t
(σ2 − ξα2 + ξα2e2σ2t) ( σ′2 − ξ′α2 + ξ′α2e2σ ′ 2t ) ≤ |σ2ξ′ − σ′2ξ| ξξ′
≲ |σ2| n ∣∣a2n − a,2n ∣∣+ ∣∣∣∣ ξn (anεn − a′nε′n) ∣∣∣∣ .
Combined the above three parts, we finally have |u2v(t)− u′2v (t)| = O ( V (t+ 1)
n
) .
In the second case, we assume both σ2, σ′2 < 0. Since both of them would decrease to 0, we have |u2b(t)− u′2b (t)| ≤ α2 where α = O( 1n1/4 ).
In the last case, without loss of generality, we assume σ2 < 0 < σ′2. Since u 2 v(t) decreases to 0, and u′ 2 v(t) increases to σ ′ 2, we have |u2b(t)− u′ 2 b(t)| ≤ σ′2 ≤ |σ′2 − σ2|.
Overall, we have
|u2v(t)− u′2v (t)| = O ( V (t+ 1)
n
) .
Besides, since the loss is O(V )-Lipschitz: ∥∇u2ℓ(u)∥ = 2 ( au2 − ε ) a ≲ V.
where u2 is bounded upper bounded by O(V ) during the training analysis due to a small initialization. Then the total stability is
|ℓ(u2v(t))− ℓ(u′2v (t))| ≤ V 2(t+ 1)
n .
Summation over the dimension d, the whole loss ℓd(u2v(t)) has stability
|ℓd(U2v (t))− ℓd(U ′2v (t))| ≤ dV 2(t+ 1)
n .
By Proposition 1, we have that
L(θ̂(t)v ;P)− L(θ̂(t)v ;Pn) ≲ dV 2(t+ 1)
n log(n) log(2dn/δ) +
√ log(1/δ)√
n .
Part (II). We next calculate the second term. Take the explicit formula into the equation, we have
Ln ( u2(t) ) − Ln (0) = 1
n n∑ i=1 ( a2iu 4 v(t)− 2aiεiu2v(t) ) .
If σ2 = 1n ∑n i=1 aiεi > 0, we have
1
n n∑ i=1 ( a2iu 4 v(t)− 2aiεiu2v(t) ) = − (σ2 − ξα 2)σ22α 2e2σ2t (σ2 − ξα2 + ξα2e2σ2t)2 ≤ 0.
If σ2 < 0, we have
1
n n∑ i=1 ( a2iu 4 v(t)− 2aiεiu2v(t) ) = ξ σ22α 4( (|σ2|+ ξα2)e2|σ2|t − ξα2 )2 + 2|σ2|2α2(|σ2|+ ξα2)e2|σ2|t − ξα2
≤ ξα4 + 2|σ2|α2 ≲ α2 √ ν2
n .
In conclusion, we have Ln ( u2(t) ) − Ln (0) ≲ α2 √ ν2
n .
Part (III). For the third term, we just need one single concentration inequality. Note that for the variance training, the loss function is of the order O ( V 2 ) . Hence, via the Hoeffding inequality, with probability at least 1− δ, we have
Ln (0)− L (0) ≲ V √ log(1/δ)
n .
Combining the three terms leads to Lemma 5.
□

heading: E.4 PROOF OF LEMMA 6
text: Lemma 6 For BER, under the assumptions in Theorem 3, with probability at least 1− δ, we have
EbL(Ut;P) ≲ ( ∥X⋆∥2F + dα 4 )√ log (d/δ)
n + r∑ j=1 σ4j σ2j + α 4eΩ(σjt) + α2d t . (9)

heading: Proof:
text: Similar to linear regression, we decompose the excess risk into two parts:
E [ a2 ( u2b(t)− σ )2] = E [ a2 − 1
n n∑ i=1 a2i ] ( u2b(t)− σ )2 + 1 n n∑ i=1 a2i ( u2b(t)− σ )2 .
We use the uniform convergence to bound the first term. Assume that σ > 0, we have
σ − u2b(t) = σ − σ1α
2e2σ1t
σ1 − ξα2 + ξα2e2σ1t
= σσ1 − ξα2σ + α2e2σ1t (ξσ − σ)
σ1 − ξα2 + ξα2e2σ1t
= σσ1 − ξα2σ
σ1 − ξα2 + ξα2e2σ1t
= σ2 − α2σ
σ − α2 + α2e2ξσt ≤ σ(1− α2) ≤ σ.
Now assume σ = 0, we have
u2b(t) = α2
ξα2t+ 1 ≤ α2.
The second term refers to the training loss, if σ > 0, we have
1
n n∑ i=1 a2i ( u2b(t)− σ )2 = ξ
( σ2 − α2σ )2 (σ − α2 + α2e2ξσt)2 ≲ σ 4 σ2 + α4eΩ(σt) .
If σ = 0, we have
1
n n∑ i=1 a2i ( u2b(t)− σ )2 = ξ
α4
(ξα2t+ 1) 2 ≲
α4
1 + α4t2 .
Finally, by assigning probability properly, we have
EL(Ut;Pb) ≲ ( ∥X⋆∥2F + dα 4 )√ log (d/δ)
n + r∑ j=1 σ4j σ2j + α 4eΩ(σjt) + α2d t .
□

heading: F ADDITIONAL DISCUSSION
text: In this section, we discuss more the details in the main text. We first validate the rationality of applying stability-based bound by showing that the generalization gap (under variance regime) increases with time. We next provide some motivation behind Dynamics Decomposition Condition by showing that it indeed holds as the beginning of the training phase. We finally introduce a related version of Dynamics Decomposition Condition which does not require the information about θ∗, θ∗v , and θ ∗.

heading: F.1 THE GENERALIZATION GAP UNDER VARIANCE INCREASES WITH TIME
text: In this part, we will prove that the (variance regime) generalization gap indeed increases with time, which validates that applying stability-based bound (which also increases with time) is rational.
We provide the following Lemma 9 which validates the above statement.
Lemma 9 Let ∆v(t) = L(θ̂(t)v ;P)− L(θ̂(t)v ;Pn) be the generalization gap under variance regime. Assume that Σx = I . Then under linear regression settings as in Section 3, we have
∂ ∂t ∆v(t) ≥ 0.
Proof: Note that the population loss is calculated as follows by Equation C.1 L(θ̂(t)v ;P) = [ θ̂(t)v ]⊤ Σx [ θ̂(t)v ] + Eϵ2.
And the training loss is calculated as follows:
L(θ̂(t)v ;Pn) = 1
n ∥E −Xθ̂(t)v ∥2
= 1
n ∥Xθ̂(t)v ∥2 −
2 n E⊤
[ Xθ̂(t)v ] + 1
n ∥E∥2
= [ θ̂(t)v ]⊤ [ 1 n X⊤X ] [ θ̂(t)v ] − 2 n E⊤Xθ̂(t)v + 1 n ∥E∥2 .
Note that by Lemma 7, we have
θ̂(t)v =
[ I − [ I − λ
n X⊤X
]t] [ X†E ] .
Therefore, we have: ∆v(t) = L(θ̂(t)v ;P)− L(θ̂(t)v ;Pn)
= [ θ̂(t)v ]⊤ Σx [ θ̂(t)v ] + Eϵ2 − [ θ̂(t)v ]⊤ [ 1 n X⊤X ] [ θ̂(t)v ] + 2 n E⊤Xθ̂(t)v − 1 n ∥E∥2
= [ θ̂(t)v ]⊤ [ Σx − 1
n X⊤X
] [ θ̂(t)v ] + 2
n E⊤Xθ̂(t)v + Eϵ2 −
1 n ∥E∥2
= E⊤ [ X† ]⊤ [ I − [ I − λ
n X⊤X
]t] [ Σx − 1
n X⊤X
] [ I − [ I − λ
n X⊤X
]t] [ X†E
] + 2
n E⊤X
[ I − [ I − λ
n X⊤X
]t] [ X†E ] + Eϵ2 − 1
n ∥E∥2.
By omitting the terms in ∆v(t) that is uncorrelated to t, we have (denote as ∆̄v(t)): ∆̄v(t) = E ⊤ [X†]⊤ [I − λ
n X⊤X
]t [ Σx − 1
n X⊤X
] [ I − λ
n X⊤X
]t [ X†E ] − 2E⊤ [ X† ]⊤ [ Σx − 1
n X⊤X
] [ I − λ
n X⊤X
]t [ X†E ] − 2 n E⊤X [ I − λ n X⊤X ]t [ X†E
] = E⊤ [ X† ]⊤ [
I − λ n X⊤X
]t [ Σx − 1
n X⊤X
] [ I − λ
n X⊤X
]t [ X†E ] − 2E⊤ [ X† ]⊤
[Σx]
[ I − λ
n X⊤X
]t [ X†E ] = E⊤ [ X† ]⊤ [[
I − λ n X⊤X
]t [ Σx − 1
n X⊤X
] − 2Σx ] [ I − λ
n X⊤X
]t [ X†E ] .
where we use the fact that [ X† ]⊤ X⊤X = X .
We next focus on the matrix Ma ≜ [ X† ]⊤ [[
I − λ n X⊤X
]t [ I − 1
n X⊤X
] − 2I ] [ I − λ
n X⊤X
]t X†,
where we plug in Σx = I by assumption. Denote the SVD of 1nX ⊤X as U⊤ΣU where the i-th eigenvalue is σi. Then when σi = 0, the i − th eigenvalue of Ma is also equal to zero (due to the persudo inverse X†). When σi > 0, the i− th eigenvalue of Ma can be derived as
σi(t;Ma) = [(1− λσi)t(1− σi)− 2] (1− λσi)t
σi .
Its derivation is then
∂ ∂t σi(t;Ma) = 2(1− λσi)t log(1− λσi) σi
[ (1− σi)(1− λσi)t − 1 ] ≥ 0, (15)
since log(1−λσi) ≤ 0 and (1−σi)(1−λσi)t−1 ≤ 0. Therefore, the derivation of each eigenvalue is larger than zero.
We rewrite ∆̄v(t) = z⊤ΣMaz⊤ where z⊤ = UE is a vector independent of time t and ΣMa is the diagonal matrix with diagonal σi(Ma). As a result, we have
∂ ∂t ∆̄v(t) = z ⊤( ∂ ∂t ΣMa)z ⊤ ≥ 0,
since ∂∂tΣMa ⪰ 0 according to Equation 15.
Note that ∆̄v(t) and ∆v(t) only differ with a term independent of time, we have
∂ ∂t ∆v(t) = ∂ ∂t ∆̄v(t) ≥ 0.
□

heading: F.2 MOTIVATION BEHIND DYNAMICS DECOMPOSITION CONDITION
text: For a dataset S = {(Xi, yi)}ni=1, where yi = fθ⋆(Xi) + εi, and θ⋆ ∈ Θ. We choose the ℓ2loss to solve this problem L(θ) = 12n ∑n i=1 (fθ(Xi)− yi)
2. Then we can show that the Dynamics Decomposition Condition holds at the beginning phase if we initialize θ(0)v = 0, and assume f0(X) = 0, ∀X ∈ supp{P}. For sufficient small t > 0, the dynamics are well approximated by
θ̂(t) = θ(0) − 1 n n∑ i=1 (fθ(0)(Xi)− yi)∇fθ(0)(Xi)t;
θ̂ (t) b = θ (0) b −
1
n n∑ i=1 ( f θ (0) b (Xi)− fθ⋆(Xi) ) ∇f θ (0) b (Xi)t;
θ̂(t)v = θ (0) v −
1
n n∑ i=1 ( f θ (0) v (Xi)− εi ) ∇f θ (0) v (Xi)t.
(16)
Suppose θ(0) = θ(0)b , then based on Taylor expansion, we have the approximation
θ̂(t) ≈ θ̂(t)b + θ̂ (t) v +
1
n n∑ i=1 εi∇2fθ(0)(Xi) ( θ(0)v − θ(0) ) t ≈ θ̂(t)b + θ̂ (t) v . (17)
Therefore, the Dynamics Decomposition Condition holds at least at the beginning phase. Moreover, we discuss in the main text that for a branch of learning tasks, such as linear regression and matrix recovery, these Dynamics Decomposition Condition uniformly holds for arbitrary 0 ≤ t <∞.

heading: F.3 RELAXED VERSION OF DDC
text: One may notice that there exist θ∗, θ∗v , θ ∗ b in DDC (See Assumption 2), which is the optimal solution to the corresponding population loss. They are sometimes unavailable in practice and even have no explicit solution. For more convenient analysis, we propose a relaxed version of DDC in Lemma 10.
Lemma 10 (Relaxed DDC) Let the loss be ℓ2 loss where ℓ(y, ŷ) = (y − ŷ)2. We assume the predicting model ŷ = fθ̂(x) is unique with respect to θ̂ (when θi ̸= θj , fθi ̸= fθj ). Besides, assume fθ=0(x) = 0 for all x. For the ground truth, we assume an additive noise, namely, y = fθ∗(x) + ϵ where ϵ is independent of x. Then the following Equation 18 suffices to prove (max{a, 1}, C, C ′)- boundedness for DDC (See Assumption 2).
∥θ̂(t) − θ̂(t)b ∥ ≤ a∥θ̂ (t) v ∥+ C√ t + C ′√ n . (18)
Based on Lemma 10, one can verify DDC without knowing the optimal solution parameter θ∗, θ∗v , θ ∗ b . Besides, since we usually have the iteration equation (the relationship between θ̂(t+1) and θ̂(t)) given an algorithm, it is possible to verify DDC based on the iteration forms. Note that we may not directly verify the condition in practice based on the relaxed DDC if we do not have the explicit split on signal E[y|x] and noise y − E[y|x].

Review summaries:
This work proposes a new bound for the excess risk based on a decomposition into a bias and variance term. The authors use a stability-based analysis to bound the bias excess risk, and uniform convergence to bound the variance excess risk. They illustrate their framework in two contexts: overparametrised linear regression with SGD and low-rank diagonal matrix recovery with gradient flow. 

This paper proposes a new approach to characterize the excess risk by decomposing the risk into two parts, the risk of learning the exact response and the risk of fitting the noise signal. Under the linear regression setting, the learning dynamic can be decomposed exactly into the proposed two parts. Stability and uniform convergence arguments are applied to two parts respectively. Similar results are shown under the more general setting, though additional assumptions are made to guarantee a reasonable dynamic decomposition. 

This paper studies the generalization performance of learning algorithms. The basic idea is to decompose the training dynamics into noise and signal components, and then consider separately the behavior of models trained with noise and signals. The authors then proposed to use the stability approach to tackle the noise components, and the uniform convergence approach to tackle the signal part. The authors consider two specific algorithms: overparameterized linear regression and diagonal matrix recovery. Experimental results are also reported.