Estimating Generalization Performance Along the Trajectory of Proximal SGD in Robust Regression

Authors: kai  tan
pierre c bellec
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 00:03:20

Abstract: This paper studies the generalization performance of iterates obtained by Gradient Descent (GD), Stochastic Gradient Descent (SGD) and their proximal variants in high-dimensional robust regression problems. The number of features is comparable to the sample size and errors may be heavy-tailed. We introduce estimators that precisely track the generalization error of the iterates along the trajectory of the iterative algorithm. These estimators are provably consistent under suitable conditions. The results are illustrated through several examples, including Huber regression, pseudo-Huber regression, and their penalized variants with non-smooth regularizer. We provide explicit generalization error estimates for iterates generated from GD and SGD, or from proximal SGD in the presence of a non-smooth regularizer. The proposed risk estimates serve as effective proxies for the actual generalization error, allowing us to determine the optimal stopping iteration that minimizes the generalization error. Extensive simulations confirm the effectiveness of the proposed generalization error estimates.

Main paper: 
 introduction: 
 Consider the linear model:EQUATIONwhere y ∈ R n is the response vector, X ∈ R n×p is the design matrix, b * ∈ R p is the unknown regression vector, and ε ∈ R n is the noise vector that we assume independent of X. The entries of ε may be heavy-tailed, for instance our working assumptions allow for infinite second moment.For the estimation of b * , we consider the following regularized optimization problem b ∈ arg minEQUATIONwhere ρ : R → R is a data-fitting loss and g : R p → R is a regularization function. In the present robust regression setting, typical examples of ρ include the Huber [14] loss ρ(r; δ) = δ 2 |r/δ| some initializer b 1 ∈ R p (typically b 1 = 0) followed by consecutive iterates b 2 , b 3 , . . . , where b t is typically obtained, for gradient descent and its variants as will be detailed below, from b t-1 and by applying an additive correction involving the gradient of the objective function. Our goal of this paper is to quantify the predictive performance of each iterate b t .We assume throughout that the covariance E[x i x ⊤ i ] = Σ of the feature vectors is finite. We measure the predictive performance of b t using the out-of-sample errorE x ⊤ new b t -x ⊤ new b * 2 | (x i , y i ) i∈[n] = ∥Σ 1/2 ( b t -b * )∥ 2where x new is a new feature vector, independent of the data (x i , y i ) i∈ [n] and has the same distribution as x i . The above squared metric is used because the noise ε i (and thus y i ) is allowed to have infinite variance, and in this case the squared prediction errorE[(x ⊤ new b t -y new ) 2 | (x i , y i ) i∈[n] ] = +∞ irrespective of the value of b t .The paper proposes to estimate the out-of-sample error ∥Σ 1/2 ( b tb * )∥ 2 of the t-th iterate using the right-hand side of the approximationEQUATIONwhere ψ : R n → R is the derivative of ρ acting component-wise on each coordinate in R n and S s is a diagonal matrix of the form S s = i∈Is e i e ⊤ i where I s ⊂ [n] is the batch for the s-th stochastic gradient update and e i ∈ R n is the i-th canonical basis vector. Here the w t,s are quantities, introduced in Section 3.3 below, that can be computed from data and do not require the knowledge of Σ. The approximation (3) is made rigorous in Theorem 3.6, where the right-hand side is proved to be consistent (i.e., the difference between the two sides of the inequality converges to 0 in probability) for a first set of weights (w s,t ) s<t , and in Theorem 3.7 where a second set of weights are proposed.Because the right-hand side of ( 3) is observable from the data and the iterates ( b s ) s≤t are computed from the iterative algorithm, the approximation (3) lets us compare the out-of-sample error of iterates b t at different time t up to the additive term ∥ε∥ 2 /n (which does not depend on t nor on the choice of the iterative scheme or the choice of loss and penalty). It also lets us compare different tuning parameters, for instance learning rate, multiplicative parameter of the penalty function, batch size in Stochastic Gradient Descent (SGD). The right-hand side of (3) can serve as the criteria to choose the iteration number or tuning parameters that achieves the smallest out-of-sample error.Estimation of prediction risk of regression estimates has received significant attention in the last few decades. One natural avenue to estimate the generalization performance is to use V -fold crossvalidation or leave-one-out schemes. In the proportional regime of interest here, where dimension p and sample size n are of the same order, V -fold cross-validation with finite V , e.g., V = 5, 10 is known to fail at consistently estimate the risk of the estimator trained on the full dataset [23, Figure 1] ; this is simply explained because training with the biased sample size n(V -1)/V may behave differently than training with the full dataset. Leave-one-out schemes, or drastically increasing V , requires numerous refitting and is thus computationally expensive. This motivates computationally efficient estimates of the risk of an estimator trained on the full dataset without sample-splitting, including Approximate Leave-One-out (ALO) schemes [23] that do not rely on sample-splitting and refitting; see [1] and the references therein for recent developments. For ridge regression and other estimators constructed from the square loss, the Generalized Cross-Validation (GCV) [28] has been shown be to be effective, and it avoids data-splitting and refitting; it only needs to fit the full data once and then adjust the training error by a multiplicative factor larger than 1. Beyond ridge regression, the extension of GCV using degrees-of-freedom has been studied for Lasso regression [2, 3, 18, 9] , and alternatives were developed for robust M-estimators [3, 4] . While ALO or GCV and its extensions are good estimators of the predictive risk of a solution b to the optimization problem (2) , they are not readily applicable to quantify the prediction risk of iterates b t obtained by widely-used iterative algorithms such as gradient descent (GD), stochastic gradient descent (SGD) or their proximal variants: ALO and GCV focus on estimating the final (t → +∞) iterate of the algorithm, when a solution b in (2) is found. Our goal in the present paper is to develop risk estimation methodologies along the trajectory of the algorithm.Luo et al. [17] developed methods to estimate the cross-validation error of iterates that solves an empirical risk minimization problem. Their approach requires the Hessian of the objective function to be well-conditioned (i.e., the smallest and largest eigenvalues are bounded) along all iterates. This condition is not satisfied for the regression problems we consider in this paper, such as highdimensional robust regression with a Lasso penalty. In the context of least squares problems with both p and n being large, [7] studied the fundamental limits on the performance of first order methods, showing that these are dominated by a specific Approximate Message Passing algorithm. Paquette et al. [19] demonstrated that the dynamics of Stochastic Gradient Descent (SGD) become deterministic in the large sample and dimensional limit, providing explicit expressions for these dynamics when the design matrix is isotropic. Our work differs from [19] in two key ways: First, we address a more general regression problem incorporating a non-smooth regularizer, thereby considering both SGD and proximal SGD; second, we offer explicit risk estimates for each iteration, rather than focusing solely on the theoretical dynamics of the iterates. Celentano et al. [8] and [13] characterize the dynamical mean-field dynamics of iterative schemes, and identify that the limiting process involves a "memory" kernel, describing how the dynamics of early iterates affect later ones.Most recently, and most closely related to the present paper, [5] proposed risk estimate for iterates b t obtained by running gradient descent and proximal gradient descent methods for solving penalized least squares optimizations. However, [5] focuses exclusively on the square loss for ρ in (3), which is not readily applicable to robust regression with heavy tailed noise for which the Huber or other robust losses must be used. Bellec and Tan [5] is further restricted to gradient updates using the full dataset, which does not cover stochastic gradient descent. While several proof techniques used in the present paper are inspired by [5] , we will explain in Remark 3.8 that directly generalizing [5] to SGD in robust regression leads to a poor risk estimate for small batch sizes. The proposal of the present paper leverages out of batch samples to overcome this issue.For gradient descent for the square loss and without penalty, Patil et al. [21] demonstrates both the failure of GCV along the trajectory and the success of computationally expensive leave-one out schemes, and develops a proposal to reduce the computational cost. Finally, let us mention the works [10, 16] that characterize the dynamics of the iterates in phase retrieval and matrix sensing problems, assuming that a fresh batch of observations (independent of all previous updates) is used at each iteration. This is different from the usual SGD setting studied in the present paper where the observations used during a stochastic gradient update may be reused in future stochastic gradient updates, creating intricate probabilistic dependence between gradient updates at different iterations.Robust regression is highly valuable in real data analysis due to its ability to handle heavy-tailed noise effectively, and we will see below that the use of stochastic gradient updates and data-fitting loss functions different from the square loss require estimates that have a drastically different structure that in the square loss case. The present paper develops generalization error estimates in situations where no consistent estimate have been proposed: (1) we develop generalization error estimates along the trajectory of iterative algorithms aimed at solving (2) for robust loss functions including the pseudo-Huber loss; (2) the estimates are applicable not only to gradient updates involving the full dataset (gradient descent and its variants), but also to SGD and proximal SGD where a random batch is used for each update. 
problem setup: 
 The paper studies iterative algorithms aimed at solving the optimization problem (2) . We consider the algorithm that generates iterates b t for t = 1, 2, ..., T according to the following iteration:EQUATIONwhere S t ∈ R n×n is the diagonal matrix S t = i∈It e i e ⊤ i for I t ⊂ [n] the t-th batch (independent of (X, y)), where ϕ t : R p → R p and ψ : R n → R n are two functions and η t is the step size. Typically, ψ : R n → R n is the componentwise application of ρ ′ (where ρ is the data-fitting loss in (2)), and the matrix S t ∈ R n×n is diagonal with elements in {0, 1} encoding the observations i ∈ [n] used in the t-th stochastic gradient update. The presence of S t and possibly nonlinear function ψ is such that the above iteration scheme is not covered by previous related works including [5, 21] , which only tackle S t = I n (full batch gradient updates) and ψ : R n → R n the identity map (ρ in (2) restricted to be the square loss). The iterative scheme (4), on the other hand, covers SGD with robust loss functions.In the next section, we first provide a few examples of algorithms encompassed in the general iteration (4). This includes Gradient Descent (GD), Stochastic Gradient Descent (SGD), and their corresponding proximal methods [20] , Proximal GD and Proximal SGD. GD and SGD are widely used in practice, while the proximal methods are particularly useful for solving the optimization problem (2) with non-smooth regularizers.If there are no penalties in (2), i.e., g(b) = 0, then the minimization problem becomes b ∈ arg minb∈R p 1 n n i=1 ρ(y i -x ⊤ i b).To solve this problem, provided ρ is differentiable, one may use gradient descent (SGD) and stochastic gradient descent (SGD).Example 2.1 (GD). The GD method consists of the following iteration:EQUATION)where ψ is the derivative of ρ acting component-wise on its argument, and η t is the step size (also known as learning rate). For the least squares loss ρ(x) = x 2 /2, we have ψ(u) = u.Example 2.2 (SGD). Suppose at t-th iteration, we use the batchEQUATIONwhere S t = i∈It e i e ⊤ i and e i is the i-th canonical vector in R n . If I t = [n] for each t, then |I t | = n and S t = I n , hence this SGD method reduces to the GD method in (5).Regularized regression is useful for high-dimensional regression problems where p is larger than n. We consider the Lasso penalty g(b) = λ∥b∥ 1 to fight for the curse of dimensionality and obtain sparse estimates (our working assumptions, on the other hand, do not assume that the ground truth b * is sparse). While GD and SGD are not directly applicable to solve the optimization problem (2) with Lasso penalty due to ∥ • ∥ 1 lacking differentiability at 0, Proximal Gradient Descent (Proximal GD) [20] and Stochastic Proximal Gradient Descent (Proximal SGD) can be used to solve this optimization with Lasso penalty.Example 2.3 (Proximal GD). For g(b) = λ∥b∥ 1 in (2), the Proximal GD gives the following iterations:b t+1 = soft ληt b t + ηt n X ⊤ ψ(y -X b t ) , where soft θ (•) applies the soft-thresholding u → sign(u)(|u| -θ) + component-wise.Example 2.4 (Proximal SGD). Similar to the Proximal GD, the Proximal SGD consists of the following iterations:b t+1 = soft ληt b t + ηt |It| X ⊤ S t ψ(y -X b t ) .Let ρ ′ : R n → R n be the function applies the derivative of ρ : R → R to each of its component, i.e., ρ ′ (u) = (ρ ′ (u 1 ), ..., ρ ′ (u n )) ψ(u) ρ ′ (u) ρ ′ (u) ρ ′ (u) ρ ′ (u) ϕ t (v) v v soft ληt (v) soft ληt (v) S t I n S t I n S tTo define the proposed estimators of the generalization error, we further define the following Jacobian matrices:D t = ∂ψ(r) ∂r r=y-X b t ∈ R n×n , D t = ∂ϕ t (v) ∂v v= b t + ηt |It| X ⊤ Stψ(y-X b t ) ∈ R p×p .Then, we have D t = I p for GD and SGD, and D t = j∈ Ŝt e j e ⊤ j for Proximal GD and Proximal SGD based on soft-thresholding, where Ŝt = {j ∈ [p] : e ⊤ j b t+1 ̸ = 0}.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Our results only holds under several assumptions discussed in the paper.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
main results: 
 Assumption 3.1. The design matrix X has i.i.d. rows from N p (0, Σ) for some positive definite matrixΣ satisfying 0 < λ min (Σ) ≤ 1 ≤ λ max (Σ) and ∥Σ∥ op ∥Σ -1 ∥ op ≤ κ. We assume Var[x ⊤ i b * ] ≤ δ 2, that is, the signal of the model ( 1) is bounded from above. Assumption 3.2. The noise ε is independent of X and has i.i.d. entries from a fixed distribution independent of n, p, with E[|ε i |] ≤ δ, that is, bounded first moment. Assumption 3.3. The data fitting loss ρ : R → R is convex, continuously differentiable and its derivative ψ is 1-Lipschitz and |ψ(x)| ≤ δ for all x ∈ R. The function ϕ t is 1-Lipschitz and satisfies ϕ t (0) = 0. The matrices S t = i∈It e i e ⊤ i , and|I t | ≥ c 0 n for some positive constant c 0 ∈ (0, 1]. Let η max = max t∈[T ] η t .Huber loss and Psuedo-Huber loss all satisfy Assumption 3.3. Assumption 3.4. The data fitting loss ρ is twice continuously differentiable with positive second derivative. Assumption 3.5. The sample size n and feature dimension p satisfy p/n ≤ γ for a constant γ ∈ (0, ∞).This subsection provides the intuition behind the definition of the estimates define below. For the sake of clarify, and in this subsection only, assume thatEQUATIONWith the above working assumptions, the validity of the estimates defined below relies on the probabilistic approximation∥ b t -b * ∥ 2 ≈ 1 n n i=1 -x ⊤ i ( b t -b * ) + p j=1 e ⊤ j ∂ b t ∂x ij 2 ,which was developed in [3] for risk estimation purposes, but outside the context of iterative algorithms.Above, e j ∈ R p is the j-th canonical basis vector. In the present noiseless case with ε = 0, the first term inside the squared norm in the right-hand side is equal to the residual y i -x ⊤ i b t , so that the above display resembles (3) . Taking this probabilistic approximation for granted, to study the second term in the right-hand side, we must understand the derivatives of b t with respect to the entries (x ij ) i∈[n],j∈[p] of X. In (4), each iterate is a relatively simple function of the previous ones, with the simplifications (7) this is b t+1 = ϕ t ( b t + X ⊤ S t ψ(y -X b t )/n). For t = 1, given that b 1 is a constant initialization, ∂ ∂xij b 2 = D 1 e j e ⊤ i S 1 ψ(y -X b 1 )/n -D 1 X ⊤ S 1 D 1 e i ( b 1 -b * ) j /n.We find in the proof, that when summing these quantities over j ∈ [p], the second term involving ( b 1 -b * ) j is negligible, and the same negligibility holds at later iterations with terms involving ( b tb * ) j (or any ( b sb * ) j , s ≤ t). By performing a similar simple calculation at the next iteration, and ignoring these terms, we find withf 1 i = e ⊤ i S 1 ψ(y -X b 1 ) and f 2 i = e ⊤ i S 2 ψ(y -X b 2 ) by the chain rule p j=1 e ⊤ j ∂ b 2 ∂x ij ≈ Tr D 1 n w2,1 f 1 i , p j=1 e ⊤ j ∂ b 3 ∂x ij ≈ Tr D 2 n w3,2 f 2 i +Tr D 2 (X ⊤ S 2 D 2 X/n -I p ) D 1 /n w3,1 f 1 i .This reveals the weights (w s,t ) s<t in (3) at iteration t = 2 and t = 3. We could continue this further by successive applications of the chain rule, although for later iterations this unrolling of the derivatives, capturing the interplay between the Jacobians D t , D t and the stochastic gradient matrix S t , becomes increasingly complex. This recursive unrolling of the derivatives can be performed numerically at the same time as the computation of the iterates. On the other hand, for the mathematical proof, for the formal definition of the weights in (3) and for the proposed estimates of the generalization error, the matrix notation defined in the next subsection exactly captures this unrolling of the derivatives.We now set up the matrix notation that captures this recursive unrolling of the derivatives by the chain rule. EQUATIONEQUATIONK = T t=1Tr(D t )e t e ⊤ t -EQUATIONwhereΓ = M -1 L(Λ ⊗ I p ) D ∈ R pT ×pT , L = T t=2 (e t e ⊤ t-1 ) ⊗ I p , Λ = T t=1 ηt |It| e t e ⊤ t , M =     I p -P 1 I p . . . . . . -P T -1 I p     and P t = D t I p -ηt |It| X ⊤ S t D t X .Although notationally involved, the purpose of these matrices is just to formalize the recursive computation of the derivatives by the chain rule mentioned in Section 3.1.For each iterate b t , define the target r t (generalization error) and its estimate rt byEQUATIONwhere w t,s := e ⊤ t W e s and W ∈ R T ×T is the matrix defined in (8) . The following shows that |r t -r t | → P 0 (convergence to 0 in probability) under suitable assumptions. Theorem 3.6 (Proved in Appendix C.1). Let Assumptions 3.1, 3.3 and 3.5 be fulfilled. Then ∀ϵ > 0,EQUATION)If additionally Assumption 3.2 holds then E[min{1, ∥ε∥ n }] → 0, so that, as n, p → +∞ while (T, γ, η max , c 0 , δ, κ, ϵ) are held fixed, the right-hand side converges to 0 and rt -r t converges to 0 in probability.This establishes that rt is consistent at estimating r t . The statement E[min{1, ∥ε∥/n}] → 0 is equivalent to ∥ε∥ 2 /n 2 → P 0 (convergence in probability), and is proved in [22] under the assumption that E|ε i | < +∞ with ε i i.i.d. from a fixed distribution; this allows Var[ε i ] = +∞ as long as the first moment is finite. The expression of W involves Σ, which is typically unknown in practice. Our next result provides a consistent estimate of W using quantities that do not require the knowledge of Σ.We propose to estimate W by W def = K -1 A where K and A are the T × T matrices defined in ( 9)- (10) . We define another estimate rt by replacing W in (11) with W = K -1 A:rt = 1 n (y -X b t ) + t-1 s=1 wt,s S s ψ(y -X b s ) 2 ,where wt,s = e ⊤ t W e s . Theorem 3.7 (Proved in Appendix C.3). Under Assumptions 3.1 and 3.3 to 3.5, for any ϵ > 0,P |r t -r t | > ϵ ≤ 2e -n/18 + max 1, C(T,γ,ηmax,c0,δ,κ) ϵ 1 √ n + E[min(1, ∥ε∥ n )] .If additionally Assumption 3.2 holds then E[min{1, ∥ε∥ n }] → 0, so that, as n, p → +∞ while (T, γ, η max , c 0 , δ, κ, ϵ) are held fixed, the right-hand side converges to 0 and rt -r t converges to 0 in probability.This establishes the consistency of rt . The simulations presented next confirm that the two estimates rt and rt both are accurate estimates of r t . The estimate rt has the advantage of not relying on the knowledge of Σ and are recommended in practice. Remark 3.8. We highlight that directly generalizing the approach in [5] would lead to the approximation A ≈ KW , where A and K are given in (23) and (26) , respectively. From A ≈ KW , obtaining an estimate of W requires inverting K. However, this inversion fails for SGD for small (but still very realistic) batch sizes of order 0.1n in simulations (see Figure 4 ). The matrix K is lower triangular, and the reason for the lack of invertibility of K can be seen in the diagonal terms equal to Tr[S t D t ] in (26) , where S t ∈ {0, 1} n×n is the diagonal matrix with 1 in position (i, i) if and only if the i-th observation is used in the t-th batch. This diagonal element of K can easily be small (or even 0) for small batches, if the batch only contains observations such that (D t ) ii is 0 or small. Let rsub t denote the estimate of the same form as rt but using the weight matrix K -1 A instead. Simulation results in Figure 4 confirm that rsub t is suboptimal compared to our proposed rt . For SGD and proximal SGD, we solved this issue regarding the invertibility of K by using out-of-batch samples in the construction of K and A, in order to avoid S t in the diagonal elements of K in equation (31). This is the key to making these estimators work for SGD and proximal SGD, and this use of out-of-batch samples is new compared to [5] (which only tackles the square loss with full-batch gradients). Remark 3.9. The constant C(T, γ, η max , c 0 , δ, κ) in the above results is not explicit. Inspection of the proof reveals that the dependence of this constant in T is currently T T , allowing T of order log(n)/ log log n before the bound becomes vacuous. Improving this dependence in T appears challenging and possibly out of reach of current tools, even for the well-studied Approximate Message Passing (AMP) algorithms. The papers [25, 24] feature for instance the same log(n)/ log log n dependence for approximating the risk of AMP. The preprint [15] offers the latest advances on the dependence on T in the bounds satisfied by AMP. It allows T ≍ poly(n) while still controlling certain AMP related quantities, although for the risk [15, equations ( 16)-( 17)] the condition required on T is still logarithmic in n. This suggests that advances on this front are possible, at least for isotropic design and specific loss and regularizer such as those studied in [15] : Lasso or Robust M-estimation with no regularizer. Since these latest advances in [15] are obtained for specific estimates (Lasso or Robust M-estimation with no regularizer), it may be possible to follow a similar strategy and improve our bounds for specific examples of iterative algorithms closer to AMP, or algorithms featuring only separable losses and penalty. We leave such improvements for specific examples for future work, as the goal of the current paper is to cover a general framework allowing iterations of the form (4) with little restriction on the nonlinear functions except being Lipschitz.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification:Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
simulation: 
 In this section, we present numerical experiments to assess the performance of the proposed risk estimates. All necessary code for reproducing these experiments is provided in the supplementary material and is publicly available in the GitHub repository https://github.com/kaitan365/ SGD-generlization-errors. Our goal is to compare the performance of the proposed risk estimates with the true risk r t for different regression methods and iterative algorithms.We generate the dataset (X, y) from the linear model (1) , that is, y = Xb * + ε. Here, the rows of X ∈ R n×p are sampled from a centered multivariate normal distribution with covariance matrix Σ = I p . The noise vector ε consists of i.i.d. entries drawn from a t distribution with two degrees of freedom so that the noise variance is infinite. The true regression vector b * ∈ R p is chosen with p/20 nonzero entries, set to a constant value such that the signal strength ∥b * ∥ 2 equals 10.We explore two scenarios of the (n, p) pairs and corresponding iterative algorithms:(i) (n, p) = (10000, 5000): In this configuration, with n much larger than p, we examine Huber regression and Pseudo-Huber regression (without penalty or soft-thresholding). Both the gradient descent (GD) and stochastic gradient descent (SGD) algorithms are implemented for each type of regression. (ii) (n, p) = (10000, 12000): Here, we investigate Huber regression and Pseudo-Huber regression with an L1 penalty, λ∥b∥ 1 (λ = 0.002) and corresponding soft-thresholding step. For each penalized regression, we employ the Proximal Gradient Descent (Proximal GD) and Stochastic Proximal Gradient Descent (Proximal SGD) algorithms.In all algorithms, we start with the initial vector b 1 = 0 p and proceed with a fixed step size η = (1 + p/n * ) -2 where n * = n for GD and proximal GD, and n * = n/5 for SGD and proximal SGD. We run each algorithm for T = 100 steps. For SGD and Proximal SGD, batches I t ⊂ {1, 2, . . . , n} are randomly sampled without replacement and independently of (X, y, (I s ) s̸ =t ), each with cardinality |I t | = n 5 . A crucial component of the proposed risk estimates rt and rt involve the weight matrices W and W . The matrix W is defined in Theorem 3.6, and W = K -1 A is defined in Theorem 3.7. We employ Hutchinson's trace approximation to compute W , A, and K. This implementation is computationally efficient. We refer readers to [5, Section 4] for more details.Recall that we have proposed two estimates forr t = ∥Σ 1/2 ( b t -b * )∥ 2 + ∥ε∥ 2 /n, one is rt in Theorem 3.6 which requires knowing Σ = E[x i x ⊤ i ],and the other is rt in Theorem 3.7 which does not need Σ. Since the quantity ∥ε∥ 2 /n remains constant along the algorithm trajectory, we only focus on the estimation of ∥Σ 1/2 ( b tb * )∥ 2 . We repeat each numerical experiment 100 times and present the aggregated results in Figures 1, 2 and 3 .In Figure 1 , we focus on the scenario with (n, p) = (10000, 5000), and plot the actual risk ∥Σ 1/2 ( b tb * )∥ 2 , and its two estimates rt -∥ε∥ 2 /n and rt -∥ε∥ 2 /n along with the 2 standard error bar for GD and SGD algorithms applied to both Huber and Pseudo Huber regression. In Figure 2 , we focus on the scenario with (n, p) = (10000, 12000), and present the risk curves for the Proximal GD and Proximal SGD algorithms applied to both L1-penalized Huber regression and Pseudo-Huber regression.Figure 1 and Figure 2 confirm the three curves are in close agreement, indicating that the proposed estimates rt -∥ε∥ 2 /n and rt -∥ε∥ 2 /n are consistent estimates of the actual risk∥Σ 1/2 ( b t -b * )∥ 2 .The two estimates closely capture the risk ∥Σ 1/2 ( b t -b * )∥ 2 over the entire trajectory of the algorithms. For GD and Proximal GD, the risk curves exhibit a U-shape, first decreasing and then increasing, and the estimates rt and rt closely capture this pattern. This suggests that the proposed estimates are reliable and can be used to monitor the risk of the iterates and find the optimal iteration (the iteration minimizing the generalization error) along the trajectory of the algorithm.Additional experiments: varying step sizes for different iterations. We also conduct simulations to investigate the accuracy of the proposed risk estimates in a setting with varying step size. We consider two types of step sizes: 1). η t = 1 if t is odd, and and η t = 0.5 if t is even. While the above choices of step size are not preferred in practice, here the goal is show that the proposed risk estimates is able to accurately capture the dynamics of the risk even when the step size changes along the trajectory of the algorithm. For instance, the first choice of step size should produce a risk curve that is flat when t is even. The results are presented in Figure 3 , illustrating that the risk estimates accurately capture the flat segments of the true risk curve.η t = 0 if t is even; 2). η t = 1 if t is odd,Additional  
discussion: 
 This paper proposes a novel risk estimate for the generalization error of iterates generated by the proximal GD and proximal SGD algorithms in robust regression. The proposed risk estimates accurately capture the predictive risk of the iterates along the trajectory of the algorithms, and are provably consistent (Theorems 3.6 and 3.7). Three matrices in R T ×T in ( 8)- (10) reveal the interplay between the squared risk, the residuals and the gradients, so that the approximation (3) holds. This structure is different from the square loss case studied in [5] where only two matrices (inverse of each other) are sufficient.Let us mention some open questions along with potential future research directions. The first question regards the probabilistic model: we currently assume Gaussian features x i , and it would be of interest to study the extension in which our consistency results are universal, allowing non-Gaussian feature distributions. Second, it is of interest to extend the current estimates to more general optimization problems of the form (2) with non-smooth data-fitting loss, for instance the Least Absolute Deviation loss ∥y -Xb∥ 1 . In this case the gradient does not exist at the origin, which calls for different algorithms than the GD and SGD variants presented here, for instance the Alternating Direction Method of Multipliers (ADMM) [6] . It is of independent interest to derive the risk estimates for iterates obtained by such primal-dual methods.A Additional simulation resultsThe following figures illustrate the proposed risk estimates accurately estimate the trajectory of the risk even when the step size changes at every step. Throughout, we defineEQUATIONEQUATION)Note that in the above E, H are not observable since ε and b * are unknown. However, F and R are observable and can be easily computed once the iterates ( b t ) t∈[T ] are calculated.In this section, we conduct the change of variable to simplify the proof. Specifically, we view the linear model y = Xb * + ε as a model with design matrix G and the regression vector θ * , i.e.y = Xb * + ε = XΣ -1/2 G Σ 1/2 b * θ * +ε.This way, the design matrix G has i.i.d. entries from standard normal distribution. Using the same argument in [5, Appendix D], we can show that the matrices H, F , A, K remains the same under the change of variable. Therefore, we can prove the main results using the model with design matrix G and the regression vector θ * . In other words, we assume without of loss of generality that the design matrix X has i.i.d. N (0, 1), or equivalently that the independent rows of X are normally distributed with covariance Σ = I p . We prove the main results using Σ = I p , and the results for general Σ follow by this change of variable with the constant C(T, γ, η max , c 0 , δ) appearing in the bounds depending additionally on κ (the upper bound of the condition number of Σ from Assumption 3.1).In this section, we present derivative formulae that will be useful in later proofs. The following formulas differ from [5] due to the use of robust loss functions and the application of SGD with random batches at each iteration. The formulae are also significantly more complex than in the case of regularized M-estimators [3, 4] .Lemma B.1 (Proved in Appendix D.1). Let ( b t ) t∈[T ] be the iterates generated from the recursion (4) and the initial value b 1 is independent of X. Then the derivative of b t with respect to X is given byEQUATIONwhereΓ = M -1 L(Λ ⊗ I p ) D, L = T t=2 (e t e ⊤ t-1 ) ⊗ I p , Λ = T t=1ηt |It| e t e ⊤ t , and l F e t = e ⊤ l S t ψ(y -X b t ), we haveM =     I p -P 1 I p . . . . . . -P T -1 I p     where P t = D t (I p -ηt |It| X ⊤ S t D t X). Recall F = [S 1 ψ(y-X b 1 ), ..., S T ψ(y-X b T )], we have F lt = e ⊤ l S t ψ(y-X b t ).EQUATIONwhereD lt ij = -e ⊤ l S t D t e i e ⊤ j He t + ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )SD(e t ⊗ e l ), ∆ lt ij = -((e ⊤ i F ) ⊗ e ⊤ j )Γ ⊤ (I T ⊗ X ⊤ )SD(e t ⊗ e l ).Lemma B.3 (Proved in Appendix D.3). Let F = [ψ(y-X b 1 ), ..., ψ(y-X b T )] and F l,t = e ⊤ l F e t . Under the same conditions of Lemma B.1. We haveEQUATIONwhereDlt ij = -e ⊤ l D t e i e ⊤ j He t + ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e l ), ∆lt ij = -((e ⊤ i F ) ⊗ e ⊤ j )Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e l ).Definition B.4. Define the matricesEQUATION∀i ∈ [n],EQUATIONEQUATIONEQUATIONEQUATIONwhere the matrices K, A, A, W are defined as followsEQUATIONEQUATIONEQUATIONEQUATIONThe matrices Υ 1 , Υ 2 , . . . are negligible in the sense that their Frobenius norms are of smaller orders compared to their preceding terms in ( 18)- (22) . We provide the bounds in next lemma, which is obtained by deriving alternative expressions for Υ 1 , ..., Υ 5 in Appendix D.4.Lemma B.5 (Proved in Appendix D.6). Under the same conditions of Theorem 3.6 with Σ = I p , we havemax k∈{1,3,4} E[∥Υ k ∥ 2 op | ε] ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ), E[∥Υ 2 ∥ 2 op | ε] ≤ n -1 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ), E[∥Υ 5 ∥ 2 op | ε] ≤ n 2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).We further define a few matrices of size T × T :EQUATIONThe next lemma provides the moment bounds for the Frobenius norm of these matrices. Lemma B.6 (Proved in Appendix D.7). Under the same conditions of Theorem 3.6, we haveEQUATIONEQUATIONEQUATIONalmost surely, where E[• | ε] is the conditional expectation given ε.We are able to prove the main theorems using Lemma B.6.C.1 Proof of Theorem 3.6It suffices to prove this theorem for the case Σ = I p . When Σ ̸ = I p , the result can be derived using a change of variables argument, as outlined in Appendix B.1. By basic algebra, we haveΘ 5 + ∥E∥ F (Θ 6 + Θ ⊤ 6 ) = nH ⊤ H -(XH -F W ⊤ ) ⊤ (XH -F W ⊤ ) + E ⊤ (XH -F W ⊤ ) + (XH -F W ⊤ ) ⊤ E = nH ⊤ H + E ⊤ E -(E -XH + F W ⊤ ) ⊤ (E -XH + F W ⊤ ) = nH ⊤ H + E ⊤ E -(R + F W ⊤ ) ⊤ (R + F W ⊤ ). Notice that r t = ∥ b t -b * ∥ 2 + ∥ε∥ 2 /n is the t-th diagonal entry of (H ⊤ H + E ⊤ E/n), and rt is the t-th diagonal entry of (R + F W ⊤ ) ⊤ (R + F W ⊤ )/n. Since ∥E∥ F = √T ∥ε∥, using the previous display that conditionally on ε, we haveE rt -r t | ε ≤ n -1 E ∥Θ 5 ∥ F + 2∥E∥ F ∥Θ 6 ∥ F | ε = n -1 E ∥Θ 5 ∥ F + 2 √ T ∥ε∥∥Θ 6 ∥ F | ε .Using the moment bounds of Θ 5 and Θ 6 in Lemma B.6, we haveE rt -r t | ε ≤ C(T, γ, η max , c 0 , δ) √ n 1 + ∥ε∥ √ n . Furthermore, if E[|ε i |]is finite, we have by [22] that ∥ε∥/n → P 0 (convergence in probability) if ε has i.i.d. entries with a fixed distribution independent of n. Under this assumption, the right-hand side of the previous display converges to 0 in probability. By enlarging the constant if necessary, assume C(T, γ, η max , c 0 , δ) ≥ 1. To obtain a quantitative bound, by the conditional version of Markov's inequality, for any ϵ > 0, and almost surely with respect to ε thatP |r t -r t | > ϵ | ε ≤ min 1, C(T, γ, η max , c 0 , δ) ϵ 1 √ n + ∥ε∥ n ≤ max 1, C(T, γ, η max , c 0 , δ) ϵ min 1, 1 √ n + ∥ε∥ n .Taking expectation with respect to ε, we obtainP |r t -r t | > ϵ ≤ max 1, C(T, γ, η max , c 0 , δ) ϵ E min 1, 1 √ n + ∥ε∥ n ≤ max 1, C(T, γ, η max , c 0 , δ) ϵ 1 √ n + E min 1, ∥ε∥ nwith E[min{1, ∥ε∥ n }] → 0 (equivalently, ∥ε∥/n → P 0) if the entries of ε are i.i.d. with a fixed distribution independent of n with E[|ε i |] < +∞ by [22] . This finishes the proof of Theorem 3.6.We first recall the definition of K from (10) in the main text:EQUATIONDefine two events:Ω 1 = {X ∈ R n×p : ∥X∥ op / √ n ≤ 2 + p/n} and Ω 2 = {|{i ∈ [n] : |ε i | ≤ M }| ≥ 2n 3 }, where M is a large enough constant such that P(|ε i | > M ) ≤ 1/6. Lemma C.1. Under the same conditions of Theorem 3.6 with Σ = I p , we have in the eventΩ * = Ω 1 ∩ Ω 2 that n∥ K -1 ∥ op ≤ C(T, γ, η max , c 0 , δ, ∥b * ∥).Furthermore, Ω * has probability at least 1 -e -n/18 -e -n/2 . Proof of Lemma C.1. Under Assumptions 3.1 and 3.5, we know that P(Ω 1 ) ≥ 1 -e -n/2 from [11] . In the event Ω 1 , we have by Lemma D.2 that∥XH∥ 2 F /n ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ) := C * . Markov's inequality further implies |{i ∈ [n] : ∥x ⊤ i H∥ 2 > 3C * }| ≤ n/3. In other words, |{i ∈ [n] : ∥x ⊤ i H∥ 2 ≤ 3C * }| ≥ 2n 3 .Recall that M is such that P(|ε i | > M ) ≤ 1/6. By Hoeffding's inequality, we haveP 1 n n i=1 1{|ε i | > M } ≥ P(|ε i | > M ) + a ≤ e -2na 2.Taking a = 1/6, we haveP 1 n n i=1 1{|ε i | > M } ≥ P(|ε i | > M ) + 1/6 ≤ e -n/18 . Furthermore, using |{i ∈ [n] : |ε i | > M }| = n i=1 1{|ε i | > M }, we have |{i ∈ [n] : |ε i | > M }| ≥ n/3| = 1 n n i=1 1{|ε i | > M } ≥ 1/6 + 1/6 ⊆ 1 n n i=1 1{|ε i | > M } ≥ P(|ε i | > M ) + 1/6 .Therefore, 18 .P |{i ∈ [n] : |ε i | > M }| ≥ n/3 ≤ e -n/Equivalently, we have 18 . That is, at least 2n 3 of the entries of ε are bounded by M with probability at least 1 -e -n/18 . 18 . Hence, P(Ω 1 ∩ Ω 2 ) ≥ 1 -e -n/18 -e -n/2 . In the eventP(|{i ∈ [n] : |ε i | ≤ M }| ≥ 2n 3 ) ≥ 1 -e -n/Ω 2 = {|{i ∈ [n] : |ε i | ≤ M }| ≥ 2n 3 }, then P(Ω 2 ) ≥ 1 -e -n/Ω 1 ∩ Ω 2 , the set Î = {i ∈ [n] : |ε i | ≤ M, ∥x ⊤ i H∥ 2 ≤ 3C * }has size at least n 3 . For any i ∈ Î and t ∈ [T ], we haveEQUATIONBy the definition of D t , under Assumption 3.4, we haveTr(D t ) = n i=1 ρ ′′ (y i -x ⊤ i b t ) > i∈ Î ρ ′′ (y i -x ⊤ i b t ) since ρ ′′ ≥ 0 by convexity ≥ | Î| min u:|u|≤M + √ 3C * ρ ′′ (u) due to (32) ≥ n/3 min u:|u|≤M + √ 3C * ρ ′′ (u) := c * n since Î has size at least n/3.Here c * is a constant depending on ρ, M, C * only.By the definition of K in (31), K/n is a lower triangular matrix with diagonal entries equal to Tr(D t )/n. It is invertible if and only if all its diagonal entries are non-zero. Therefore, in the eventEQUATIONwhereL = n i=1 (I T ⊗e ⊤ i )D(I T ⊗X)Γ(I T ⊗X ⊤ )SD(I T ⊗e i ).Here L is a strictly lower triangular matrix. Using the upper bound of ∥Γ∥ op in Lemma D.4, we have ∥ L∥ op ≤ nC(T, γ, η max , c 0 ) in the event Ω 1 . Now we rewrite K -1 asK -1 = ( K Λ -1 Λ) -1 = Λ -1 ( K Λ -1 ) -1 = Λ -1 (I T -L Λ -1 ) -1 .Notice that L Λ -1 ∈ R T ×T is a strictly lower triangular matrix, thus(I T -L Λ -1 ) -1 = ∞ k=0 ( L Λ -1 ) k = T -1 k=0 ( L Λ -1 ) k .By the triangle inequality,∥(I T -L Λ -1 ) -1 ∥ op ≤ T -1 k=0 ∥ L Λ -1 ∥ k op ≤ C(T, γ, η max , c 0 , δ, ∥b * ∥).Therefore, in the event Ω 1 ∩ Ω 2 which has probability P(Ω 1 ∩ Ω 2 ) ≥ 1 -e -n/18 -e -n/2 , ∥ K -1 ∥ op ≤ ∥ Λ -1 ∥ op ∥(I T -L Λ -1 ) -1 ∥ op ≤ n -1 C(T, γ, η max , c 0 , δ, ∥b * ∥). Lemma C.2.E[∥F ⊤ F ( A -KW ) ⊤ ∥ F | ε] ≤ n 3/2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ), E[∥R ⊤ F ( A -KW ) ⊤ ∥ F | ε] ≤ n 2 ( ∥ε∥ n + 1 √ n )C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Proof of Lemma C.3. First, using the definitions of Θ 1 , Θ 2 , Θ 4 in ( 27), we haven -1 F ⊤ F ( A -KW ) ⊤ = n -1 Θ 1 K ⊤ + Θ 2 + Θ 4 . (34) Hence, E[∥F ⊤ F ( A -KW ) ⊤ ∥ F | ε] = E[∥(Θ 1 K ⊤ + nΘ 2 + nΘ 4 )∥ F | ε]by ( 34) For the second inequality, we define≤ E[∥Θ 1 ∥ 2 F | ε] 1/2 E[∥ K∥ 2 op | ε] 1/2 + nE[∥Θ 2 ∥ F + ∥Θ 4 ∥ F | ε] by the Cauchy-Schwarz inequality ≤ n 3/2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Θ1 = R ⊤ XH + ǨH ⊤ H -R ⊤ F W ⊤ ∥E∥ F / √ n + 1 , Θ2 = R ⊤ XX ⊤ F + ǨH ⊤ X ⊤ F -pR ⊤ F + R ⊤ F A ⊤ n(∥E∥ F / √ n + 1) , Θ4 = pR ⊤ F -( ǨH ⊤ + R ⊤ X)( KH ⊤ + F ⊤ X) ⊤ n(∥E∥ F / √ n + 1) , where Ǩ = nI T - n i=1 (I T ⊗ (e ⊤ i X))Γ(I T ⊗ X ⊤ )SD(I T ⊗ (X ⊤ e i )). Using similar argument that proves Lemma B.6, we can obtain the following bound of Θ1 , Θ2 , Θ4 .EQUATIONEQUATIONBy the definitions of Θ1 , Θ2 , Θ4 , we have(∥E∥ F / √ n + 1)[n -1 Θ1 K ⊤ + Θ2 + Θ4 ] = n -1 R ⊤ F ( A -KW ) ⊤ .Therefore, conditional on ε, we haveE[∥R ⊤ F ( A -KW ) ⊤ ∥ F | ε] = (∥E∥ F / √ n + 1)E[∥ Θ1 K ⊤ + n Θ2 + n Θ4 ∥ F | ε] ≤ n 3/2 (∥E∥ F / √ n + 1)C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).This finishes the proof of Lemma C.3.In the eventΩ * = Ω 1 ∩ Ω 2 , we know from Lemma C.1 that K is invertible and ∥ K -1 ∥ op ≤ n -1 C. Define W = K -1 A. Using R + F W ⊤ = R + F W ⊤ + F ( W -W ) ⊤ , we have (R + F W ⊤ ) ⊤ (R + F W ⊤ ) -(R + F W ⊤ ) ⊤ (R + F W ⊤ ) = (R + F W ⊤ ) ⊤ F ( W -W ) ⊤ + ( W -W )F ⊤ (R + F W ⊤ ) + ( W -W )F ⊤ F ( W -W ) ⊤ .We have by the triangle inequality∥(R + F W ⊤ ) ⊤ (R + F W ⊤ ) -(R + F W ⊤ ) ⊤ (R + F W ⊤ )∥ F ≤ 2∥(R + F W ⊤ ) ⊤ F ( W -W ) ⊤ ∥ F + ∥( W -W )F ⊤ F ( W -W )∥ F ≲ ∥R ⊤ F ( W -W ) ⊤ ∥ F + ∥W F ⊤ F ( W -W ) ⊤ ∥ F + ∥( W -W )F ⊤ F ( W -W )∥ F .Recall that in the event Ω * , we have ∥ K-1 ∥ op ≤ n -1 C. Using A -KW = K( W -W ) and Lemma C.3, we have E I(Ω * )∥(R + F W ⊤ ) ⊤ (R + F W ⊤ ) -(R + F W ⊤ ) ⊤ (R + F W ⊤ )∥ F ε ≲ E I(Ω * )∥R ⊤ F ( W -W ) ⊤ ∥ F ε + E I(Ω * )∥W F ⊤ F ( W -W ) ⊤ ∥ F ε + E I(Ω * )∥( W -W )F ⊤ F ( W -W ) ⊤ ∥ F ε = E I(Ω * )∥R ⊤ F ( A -KW ) ⊤ ( K ⊤ ) -1 ∥ F ε + E I(Ω * )∥W F ⊤ F ( A -KW ) ⊤ ( K ⊤ ) -1 ∥ F ε + E I(Ω * )∥ K -1 ( A -KW )F ⊤ F ( A -KW ) ⊤ ( K ⊤ ) -1 ∥ F ε .According to Lemma C.3 and the bound of ∥ K -1 ∥ op in Lemma C.1, the first conditional expectation is bounded from above byn( ∥ε∥ n + 1 √ n )C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Using the bound of ∥K -1 ∥ op in Lemma C.1, the bound of ∥W ∥ op in Lemma C.2, and the bound ofE[∥F ⊤ F ( A -KW ) ⊤ ∥ F | ε] in Lemma C.3, the second conditional expectation is bounded from above by n 1/2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Similarly, the third conditional expectation is bounded from above byn 1/2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).In summary, we haven -1 E[I(Ω * )∥(R + F W ⊤ ) ⊤ (R + F W ⊤ ) -(R + F W ⊤ ) ⊤ (R + F W ⊤ )∥ F | ε] ≤ 1 √ n ∥ε∥ √ n + 1 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Since rt and rt are the t-th diagonal entries of (R+F W ⊤ ) ⊤ (R+F W ⊤ )/n and (R+F W ⊤ ) ⊤ (R+ F W ⊤ )/n, respectively, we haveE[I(Ω * )|r t -rt | | ε] ≤ E I(Ω * )∥(R + F W ⊤ ) ⊤ (R + F W ⊤ ) -(R + F W ⊤ ) ⊤ (R + F W ⊤ )∥ F ε ≤ 1 √ n ∥ε∥ √ n + 1 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2). Using the same argument in the proof of Theorem 3.6, we have for any ϵ > 0,P I(Ω * )|r t -rt | > ϵ | ε ≤ min 1, C(T,γ,ηmax,c0,δ) ϵ ( 1 √ n + ∥ε∥ n ) ≤ max{1, C(T,γ,ηmax,c0,δ) ϵ } min(1, 1 √ n + ∥ε∥ n ).Taking expectation with respect to ε, we haveP I(Ω * )|r t -rt | > ϵ ≤ max{1, C(T,γ,ηmax,c0,δ) ϵ }E[min(1, 1 √ n + ∥ε∥ n )].Using the union bound and P(Ω * ) ≥ 1 -e -n/18 -e -n/2 ≥ 1 -2e -n/18 , we obtainP |r t -rt | > ϵ ≤ 2e -n/18 + max{1, C(T,γ,ηmax,c0,δ) ϵ } 1 √ n + E[min(1, ∥ε∥ n )] .Using the triangle inequality and the tail probability of |r t -r t | in Theorem 3.6, we haveP |r t -r t | > ϵ ≤ 2e -n/18 + max{1, C(T,γ,ηmax,c0,δ) ϵ } 1 √ n + E[min(1, ∥ε∥ n )] . D Proof of results in Appendix B.2 D.1 Proof of Lemma B.1By assumption, we know b 1 is independent of X and b t+1 = ϕ t b t + ηt |It| X ⊤ S t ψ(y -X b t ) from (4). Recall that D t = ∂ψ(u) ∂u u=y-X b t and D t = ∂ϕt(v) ∂v v= b t + ηt |It| X ⊤ Stψ(y-X b t ) . Let a dot denote the derivative with respect to x ij . By product rule and chain rule and using y-X b t = ε-X( b t -b * ), we have ḃt+1 = D t ḃt + η t |I t | Ẋ⊤ S t ψ(y -X b t ) -X ⊤ S t D t ( Ẋ( b t -b * ) + X ḃt ) = D t ḃt + η t |I t | Ẋ⊤ F t -X ⊤ S t D t ( ẊH t + X ḃt ) ,where the last line uses F t = S t ψ(y -X b t ) and H t = b tb * . Arranging terms gives-D t (I p -ηt |It| X ⊤ S t D t X) ḃt + ḃt+1 = ηt |It| D t ( Ẋ⊤ F t -X ⊤ S t D t ẊH t ). Let P t = D t (I p -ηt |It| X ⊤ S t D t X) and a t = ηt |It| D t ( Ẋ⊤ F t -X ⊤ S t D t ẊH t ), we can rewrite the above recursion of ḃt as a linear system:    I p -P 1 I p . . . . . . -P T -1 I p     M      ḃ1 ḃ2 . . . ḃT      =     0 I p 0 . . . . . . I p 0     L     a 1 a 2 . . . a T     a .Solving the above system, we have ḃt = (e ⊤ t ⊗ I p )M -1 La. Since Ẋ = ∂X ∂xij = e i e ⊤ j ,= (Λ ⊗ I p ) D vec(e j e ⊤ i F ) -(I T ⊗ X ⊤ )SD vec(e i e ⊤ j H) = (Λ ⊗ I p ) D ((F ⊤ e i ) ⊗ e j ) -(I T ⊗ X ⊤ )SD((H ⊤ e j ) ⊗ e i ) .a into ḃt = (e ⊤ t ⊗ I p )M -1 La gives ∂ b t ∂x ij = (e ⊤ t ⊗ I p )M -1 L(Λ ⊗ I p ) D[((F ⊤ e i ) ⊗ e j ) -(I T ⊗ X ⊤ )SD((H ⊤ e j ) ⊗ e i )].This finishes the proof of Lemma B.1.By definition, F lt = e ⊤ l S t ψ(y -X b t ). By the chain rule of differentiation, we have∂F lt ∂x ij = e ⊤ l ∂S t ψ(y -X b t ) ∂x ij = -e ⊤ l S t D t (e i e ⊤ j He t + X ∂ b t ∂x ij ).Notice that (e t ⊗ (D t S t )) = DS(e t ⊗ I n ). The desired formula then follows by plugging in the expression of ∂ b t ∂xij in Lemma B.1.The desired identity follows by∂ Flt ∂x ij = e ⊤ l ∂ψ(y -X b t ) ∂x ij = -e ⊤ l D t (e i e ⊤ j He t + X ∂ b t ∂x ij )and the expression of ∂ b t ∂xij in Lemma B.1.This section derives alternative expressions for the matrices Υ 1 , ..., Υ 5 defined in Definition B.4.We first study Υ 1 in (18) . Using F = T t=1 F e t e ⊤ t , we have by Lemma B.2EQUATIONNow we compute the two terms in the right-hand side of (37). For the first term, using the expression ofD lt ij in Lemma B.2, i,t D it ij e ⊤ t = -e ⊤ j H T t=1 Tr(S t D t )e t e ⊤ t + i,t ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (e t ⊗ (X ⊤ D t S t e i ))e ⊤ t = -e ⊤ j H T t=1 Tr(S t D t )e t e ⊤ t + e ⊤ j H n i=1 (I T ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )DS(I T ⊗ e i )= -e ⊤ j H K ⊤ by (26) .Next, we compute the second term in the right-hand side of (37). Using the expression of∆ lt ij in Lemma B.2, i,t ∆ it ij e ⊤ t = - i,t ((e ⊤ i F ) ⊗ e ⊤ j )Γ ⊤ (e t ⊗ (X ⊤ D t S t e i ))e ⊤ t = -e ⊤ j i ((e ⊤ i F ) ⊗ I p )Γ ⊤ (I T ⊗ X ⊤ )DS(I T ⊗ e i ) Υ1.The identity (18) then follows by substituting the above two expressions into (37).To study Υ 2 in ( 19), we use a similar procedure. Using the mixed property of Kronecker product and the fact that the transpose of a scalar remains the same, we havep j=1 ∂e ⊤ j H ∂x ij = p j=1 T t=1 ∂e ⊤ j b t ∂x ij e ⊤ t = j,t e ⊤ j (e ⊤ t ⊗ I p )Γ((F ⊤ e i ) ⊗ e j )e ⊤ t - j,t e ⊤ j (e ⊤ t ⊗ I p )Γ(I T ⊗ X ⊤ )SD((H ⊤ e j ) ⊗ e i )e ⊤ t by (15) = e ⊤ i F j (I T ⊗ e ⊤ j )Γ ⊤ (I T ⊗ e j ) - j,t ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (e t ⊗ e j )e ⊤ t = e ⊤ i F j (I T ⊗ e ⊤ j )Γ ⊤ (I T ⊗ e j ) -e ⊤ i j ((e ⊤ j H) ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ e j ) = e ⊤ i F W ⊤ -e ⊤ i Υ 2 ,where W = j (I T ⊗ e ⊤ j )Γ(I T ⊗ e j ) andΥ 2 = j ((e ⊤ j H) ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ e j ).To study Υ 3 in (20), we use the product rule of differentiation and ( 18)-( 19):n i=1 p j=1 ∂F ⊤ e i e ⊤ j H ∂x ij = n i=1 p j=1 ∂F ⊤ ∂x ij e i e ⊤ j H + F ⊤ n i=1 p j=1 e i ∂e ⊤ j H ∂x ij = - j ( KH ⊤ e j + Υ ⊤ 1 e j )e ⊤ j H + F ⊤ i e i (e ⊤ i F W ⊤ -e ⊤ i Υ 2 ) = -KH ⊤ H + F ⊤ F W ⊤ -(Υ ⊤ 1 H + F ⊤ Υ 2 ) Υ3 .For Υ 4 in ( 21), by the product rule,n i=1 p j=1 ∂H ⊤ X ⊤ e i e ⊤ j H ∂x ij = i,j ∂H ⊤ ∂x ij X ⊤ e i e ⊤ j H + H ⊤ i,j e j e ⊤ i e i e ⊤ j H + H ⊤ X ⊤ i,j e i e ⊤ j ∂H ∂x ij = i,j ∂H ⊤ ∂x ij X ⊤ e i e ⊤ j H + nH ⊤ H + H ⊤ X ⊤ (F W ⊤ -Υ 2 )by (19) .We then compute the first term of the last line as follows i,j∂H ⊤ ∂x ij X ⊤ e i e ⊤ j H = i,j,t e t e ⊤ t ∂H ⊤ ∂x ij X ⊤ e i e ⊤ j H = i,j,t e t ∂ b t ∂x ij ⊤ X ⊤ e i e ⊤ j H = i,j,t e t e ⊤ i X ∂ b t ∂x ij e ⊤ j H = - i,j,t e t e ⊤ i X(e ⊤ t ⊗ I p )Γ(I T ⊗ X ⊤ )SD((H ⊤ e j ) ⊗ e i )e ⊤ j H + Υ 1 by (15) = - i (I T ⊗ (e ⊤ i X))Γ(I T ⊗ X ⊤ )SD(I T ⊗ e i )H ⊤ H + Υ 1 = -AH ⊤ H + Υ 1 by (23),whereΥ 1 = i,j,t e t e ⊤ i X(e ⊤ t ⊗ I p )Γ((F ⊤ e i ) ⊗ e j )e ⊤ j H = i (I T ⊗ e ⊤ i X)Γ((F ⊤ e i ) ⊗ I p )H.Combining the above pieces, we haven i=1 p j=1 ∂H ⊤ X ⊤ e i e ⊤ j H ∂x ij = -AH ⊤ H + Υ 1 + nH ⊤ H + H ⊤ X ⊤ (F W ⊤ -Υ 2 ) = (nI T -A)H ⊤ H + H ⊤ X ⊤ F W ⊤ -(H ⊤ X ⊤ Υ 2 -Υ 1 )Υ4.This provides an alternative expression for Υ 4 in (21) .Last, we study Υ 5 in (22) . We haven i=1 p j=1 ∂F ⊤ e i e ⊤ j X ⊤ F ∂x ij = i,j ∂F ⊤ e i ∂x ij e ⊤ j X ⊤ F + F ⊤ e i e ⊤ j e j e ⊤ i F + F ⊤ e i e ⊤ j X ⊤ ∂ F ∂x ij = -( KH ⊤ + Υ ⊤ 1 )X ⊤ F + pF ⊤ F + i,j F ⊤ e i e ⊤ j X ⊤ ∂ F ∂x ij .It remains to compute the third term in the last display.Using the fact that F = n l=1T t=1 e l e ⊤ l F e t e ⊤ t , we have by Lemma B.3 that∂ F ∂x ij = l,t e l ∂e ⊤ l F e t ∂x ij e ⊤ t = l,t e l ( Dlt ij + ∆lt ij )e ⊤ t .UsingDlt ij = -e ⊤ l D t e i e ⊤ j He t + ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e l ), ∆lt ij = -((e ⊤ i F ) ⊗ e ⊤ j )Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e ),we find i,jF ⊤ e i e ⊤ j X ⊤ ∂ F ∂x ij = i,j,l,t F ⊤ e i e ⊤ j X ⊤ e l ∂e ⊤ l F e t ∂x ij e ⊤ t = i,j,l,t F ⊤ e i e ⊤ j X ⊤ e l ∆lt ij e ⊤ t + i,j,l,t F ⊤ e i e ⊤ j X ⊤ e l Dlt ij e ⊤ t .We now compute the two terms in the above display. For the first term, we have i,j,l,t (24) .F ⊤ e i e ⊤ j X ⊤ e l ∆lt ij e ⊤ t = i,j,l,t F ⊤ e i e ⊤ j X ⊤ e l -((e ⊤ i F ) ⊗ e ⊤ j )Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e l ) e ⊤ t = -F ⊤ F n l=1 (I T ⊗ (e ⊤ l X))Γ ⊤ (I T ⊗ X ⊤ )D(I T ⊗ e l ) = -F ⊤ F A ⊤ byFor the second term, we have i,j,l,tF ⊤ e i e ⊤ j X ⊤ e l Dlt ij e ⊤ t = i,j,l,t F ⊤ e i e ⊤ j X ⊤ e l [-e ⊤ l D t e i e ⊤ j He t + ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (e t ⊗ (X ⊤ D t e l ))]e ⊤ t = - T t=1 F ⊤ D t XHe t + n l=1 (e ⊤ l XH ⊗ F ⊤ )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )D(I T ⊗ e l ) Υ2.Thus, we have establishedi=1 p j=1 ∂F ⊤ e i e ⊤ j X ⊤ F ∂x ij = -KH ⊤ X ⊤ F + pF ⊤ F -F ⊤ F A ⊤ -(Υ ⊤ 1 X ⊤ F -Υ 2 )Υ5.This provides an alternative expression for Υ 5 in (22) .Lemma D.1 (Moment bounds for H, F , F ). Under Assumptions 3.1, 3.3 and 3.5 with Σ = I p . Let H, F be defined in (13), we have for any finite integer k,E[∥X/ √ n∥ 2k op ] ≤ C(γ, k), E[∥H∥ 2k F | ε] ≤ C(T, γ, c 0 , η max , k)(δ 2 + ∥b * ∥) 2k , E[∥F / √ n∥ 2k F | ε] ≤ E[∥ F / √ n∥ 2k F | ε] ≤ C(T, k)δ 2k .Proof of Lemma D.1. For the first inequality, according to Assumption 3.1 and Σ = I p , X has i.i.d. standard normal entries. By [11, Theorem II.13] , there exists a random variable z ∼ N (0, 1) such that ∥X∥ op ≤ √ n + √ p + z almost surely. Under Assumption 3.5 that p/n ≤ γ, we haveE[∥X/ √ n∥ k op ] ≤ C(γ, k) for any finite integer k.For the second inequality, since∥H∥ 2 F = T t=1 ∥ b t -b * ∥ 2 , it suffices to bound ∥ b t -b * ∥ 2 for each t ∈ [T ]. Define the sequence of scalars a t def = max{∥ b t ∥, δ}. Since b t = ϕ t-1 ( b t-1 + ηt nt X ⊤ S t ψ(y - X b t-1 )) where n t := |I t |. Note that ∥ψ(y It -X It b)∥ ≤ |I t |δ by Assumption 3.3, we have ∥ b t -0∥ = ∥ b t -ϕ t-1 (0)∥ since ϕ t-1 = 0 ≤ ∥ b t-1 + ηt |It| X ⊤ S t ψ(y -X b t-1 )∥ since ϕ t is 1-Lipschitz = ∥ b t-1 + ηt |It| X ⊤ It ψ(y It -X It b t-1 )∥ ≤ ∥ b t-1 ∥ + ηt √ |It| ∥X It ∥ op δ by the triangle inequality ≤ ∥ b t-1 ∥ + ηmax √ c0n ∥X∥ op δ by |I t | ≥ c 0 n ≤ a t-1 + ηmax √ c0 ∥X/ √ n∥ op a t-1 = (1 + ηmax √ c0 ∥X/ √ n∥ op )a t-1 .Since a t = max{∥ b t ∥, δ} and δ ≤ a t-1 , we havea t ≤ (1 + ηmax √ c0 ∥X/ √ n∥ op )a t-1 .Notice a 1 = δ since b 1 = 0 p , we obtaina t ≤ (1 + ηmax √ c0 ∥X/ √ n∥ op ) t-1 δ.Hence, using the inequalityEQUATIONTaking conditional expectation on both sides given ε, the desired moment bound for H follows from the moment bound for ∥X/ √ n∥ op .For the third inequality, since |ψ(x)| ≤ δ from Assumption 3.3, we have ∥ψ(u)∥ ≤ nδ 2 for any u ∈ R n . By the definitions of F and F , we have∥F ∥ 2 F ≤ ∥ F ∥ 2 F = T t=1 ∥ψ(y -X b t )∥ 2 ≤ T nδ 2 .Since the above display holds for any ε, it implies the desired conditional moment bounds for F , F .Lemma D.2 (Frobenius norm bound for XH). Under Assumptions 3.1, 3.3 and 3.5, we have∥XH∥ 2 F ≤ C(T, γ, η max , c 0 )n(δ 2 + ∥b * ∥ 2 ) with probability at least 1 -exp(-n/2).Proof of Lemma D.2. By [11, Theorem II.13], under Assumption 3.1 with Σ = I p , we haveP(∥X/ √ n∥ op ≤ 2 + √ γ) ≥ 1 -exp(-n/2).Using ∥XH∥ 2 F ≤ ∥X∥ 2 op ∥H∥ 2 F and the bound (38), we have∥XH∥ 2 F ≤ nC(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ) holds with probability at least 1 -exp(-n/2).Lemma D.3 (Operator norm bound for M). Under Assumptions 3.1, 3.3 and 3.5, we have∥M -1 ∥ op ≤ C(T )(1 + ξ) T , where ξ = ηmax c0n ∥X∥ 2 op .Proof of Lemma D.3. By the definition of M in Lemma B.1, we haveM =     I p -P 1 I p . . . . . . -P T -1 I p     where P t = D t (I p -ηt |It| X ⊤ S t D t X).Hence, we can write M = I pT -A, where A is the lower triangular matrix with off-diagonal blocks P 1 , ..., P T -1 . Using the matrix identity(I -A) -1 = ∞ k=0 A k and noticing A k = 0 for k ≥ T , we have M -1 = T -1 k=0     0 P 1 0 . . . . . . P T -1 0     k .Taking operator norm on both sides, we obtainEQUATION) Since D t = ∂ψ(u) ∂u u=y-X b t and ψ is 1-Lipschitz, we have ∥D t ∥ op ≤ 1. By the definition that S t = i∈It e i e ⊤ i , we know ∥S t ∥ op ≤ 1. Since |I t | ≥ c 0 n and η t ≤ η max for any t ∈ [T ],we have∥P t ∥ op ≤ 1 + η t |I t | ∥X It ∥ 2 op ≤ 1 + η max c 0 n ∥X∥ 2 op def = 1 + ξ.Plugging this inequality into (39) gives∥M -1 ∥ op ≤ T -1 k=0 (T (1 + ξ)) k ≤ C(T )(1 + ξ) T .Lemma D.4. Under the same conditions as Lemma D.3, we have ∥Γ∥ op ≤ n -1 C(T, η max , c 0 )(1 + ξ)∥Γ∥ op ≤ ∥M -1 ∥ op ∥L∥ op ∥Λ ⊗ I p ∥ op ∥ D∥ op ≤ n -1 C(T, η max , c 0 )(1 + ξ) T .Lemma D.5 (Moment bounds for derivative of H, F , F ). Under Assumptions 3.1, 3.3 and 3.5 and Σ = I p , we have for any finite integer k,E n i=1 p j=1 ∂H ∂x ij 2 F k ε ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ) 2k , E n i=1 p j=1 ∂F / √ n ∂x ij 2 F k ε ≤ E n i=1 p j=1 ∂ F / √ n ∂x ij 2 F k ε ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ) 2k .Proof of Lemma D.5. We first prove the first bound. By Lemma B.1, we have∂ b t ∂x ij = (e ⊤ t ⊗ I p )Γ[((F ⊤ e i ) ⊗ e j ) -(I T ⊗ X ⊤ )SD((H ⊤ e j ) ⊗ e i )].Hence, using∂e ⊤ k Het ∂xij = ∂e ⊤ k b t ∂xij , we have ∂e ⊤ k He t ∂x ij = (e ⊤ t ⊗ e ⊤ k )Γ[((F ⊤ e i ) ⊗ e j ) -(I T ⊗ X ⊤ )SD((H ⊤ e j ) ⊗ e i )] = (e ⊤ t ⊗ e ⊤ k )Γ[(F ⊤ ⊗ I p )(e i ⊗ e j ) -(I T ⊗ X ⊤ )SD(H ⊤ ⊗ I n )(e j ⊗ e i )]. Using the above equality, i,j,t,k [(e ⊤ t ⊗ e ⊤ k )A(e j ⊗ e i )] 2 = ∥A∥ 2 F for A ∈ R pT ×np , and the triangle inequality, we haven i=1 p j=1 ∂H ∂x ij 2 F = n i=1 p j=1 p k=1 T t=1 ∂e ⊤ k He t ∂x ij 2 ≲ ∥Γ(F ⊤ ⊗ I p )∥ 2 F + ∥Γ(I T ⊗ X ⊤ )SD(H ⊤ ⊗ I n )∥ 2 F ≤ p∥Γ∥ 2 op ∥F ∥ 2 F + n∥Γ∥ 2 op ∥X∥ 2 op ∥S∥ 2 op ∥D∥ 2 op ∥H∥ 2 F ≤ p∥Γ∥ 2 op ∥F ∥ 2 F + n∥Γ∥ 2 op ∥X∥ 2 op ∥H∥ 2 F ≤ C(T, γ, η max , c 0 )(1 + ξ) 2T ∥F ∥ 2 F /n + C(T, γ, η max , c 0 )(1 + ξ) 2T ∥X∥ 2 op /n∥H∥ 2 F, where the last inequality uses Lemma D.4. Taking the conditional expectation on both sides given ε, the desired moment bound follows from the moment bounds of X, H, F in Lemma D.1. Now we prove the second bound. By definition, the t-th column of F is F t = S t ψ(y -X b t ), it can be written using the t-th column of F as F t = S t Ft . Since ∥S t ∥ op ≤ 1, we have∥ ∂F ∂x ij ∥ 2 F = t ∥ ∂F t ∂x ij ∥ 2 ≤ t ∥ ∂ Ft ∂x ij ∥ 2 = ∥ ∂ F ∂x ij ∥ 2 F .By Lemma B.3, we have∂e ⊤ l F et ∂xij = Dlt ij + ∆lt ij where Dlt ij = -e ⊤ l D t e i e ⊤ j He t + ((e ⊤ j H) ⊗ e ⊤ i )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e l ), = ((e ⊤ j H) ⊗ e ⊤ i )[-I pT + DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )]D(e t ⊗ e l ), ∆lt ij = -((e ⊤ i F ) ⊗ e ⊤ j )Γ ⊤ (I T ⊗ X ⊤ )D(e t ⊗ e l ).Using i,j,t,k [(e ⊤ j ⊗ e ⊤ i )A(e t ⊗ e l )] 2 = ∥A∥ 2 F for A ∈ R np×nT and the triangle inequality, we haven i=1 p j=1 ∥ ∂ F ∂x ij ∥ 2 F = n i=1 p j=1 n l=1 T t=1 ( Dlt ij + ∆lt ij ) 2 ≲ ∥(H ⊗ I n )[-I pT + DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )]D∥ 2 F + ∥(F ⊗ I p )Γ ⊤ (I T ⊗ X ⊤ )D∥ 2 F ≲ n∥H∥ 2 F (1 + ∥X∥ 4 op ∥Γ∥ 2 op ) + ∥F ∥ 2 F ∥Γ∥ 2 op ∥X∥ 2 op ,where the last inequality uses ∥D∥ op ≤ 1 and ∥S∥ op ≤ 1. Taking the conditional expectation on both sides given ε, the desired moment bound follows from the bound of Γ in Lemma D.4 and the moment bounds of X, H, F in Lemma D.1.Bound of Υ 1 . By the expression for the matrix Υ 1 ∈ R p×T obtained in Appendix D.4, Υ 1 = n i=1 ((e ⊤ i F ) ⊗ I p )Γ ⊤ (I T ⊗ X ⊤ )DS(I T ⊗ e i ) = n i=1 T t=1 ((e ⊤ i F e t e ⊤ t ) ⊗ I p )Γ ⊤ (I T ⊗ X ⊤ )DS(I T ⊗ e i ) = n i=1 T t=1 (e ⊤ t ⊗ I p )Γ ⊤ (I T ⊗ X ⊤ )DS(I T ⊗ e i )e ⊤ i FE[∥Υ 1 ∥ 2 op | ε] ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Bound of Υ 2 . By the expression for the matrix Υ 2 ∈ R n×T obtained in Appendix D.4, we haveΥ 2 = p j=1 ((e ⊤ j H) ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ e j ) = p j=1 T t=1 ((e ⊤ j He t e ⊤ t ) ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ e j ) = p j=1 T t=1 (e ⊤ t ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ e j )e ⊤ j He t = T t=1 (e ⊤ t ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ (He t )).By the triangle inequality and ∥D∥ op ∨ ∥S∥ op ≤ 1, ∥Υ 2 ∥ op ≤ T ∥Γ∥ op ∥X∥ op ∥H∥ F . Similar to the moment bound of ∥Υ 1 ∥ op , we obtainE[∥Υ 2 ∥ 2 op | ε] ≤ C(T, γ, η max , c 0 )n -1 (δ 2 + ∥b * ∥ 2 ).Bound of Υ 3 . By the expression for the matrix Υ 3 ∈ R T ×T obtained in Appendix D.4, we haveΥ 3 = (Υ ⊤ 1 H + F ⊤ Υ 2 ). It directly follows that ∥Υ 3 ∥ op ≤ ∥Υ 1 ∥ op ∥H∥ F + ∥F ∥ F ∥Υ 2 ∥ op .Using the triangle inequality and the moment bounds of ∥H∥ F , ∥F ∥ F in Lemma D.1 and the moment bounds of ∥Υ 1 ∥ op , ∥Υ 2 ∥ op we just obtained, we haveE[∥Υ 3 ∥ 2 op | ε] ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Bound of Υ 4 . By the expression for the matrix Υ 4 ∈ R T ×T obtained in Appendix D.4, we haveΥ 4 = H ⊤ X ⊤ Υ 2 -Υ 1 ,whereΥ 1 = i (I T ⊗ e ⊤ i X)Γ((F ⊤ e i ) ⊗ I p )H.We can rewrite Υ 1 asΥ 1 = n i=1 (I T ⊗ e ⊤ i X)Γ((F ⊤ e i ) ⊗ I p )H = n i=1 T t=1 (I T ⊗ e ⊤ i X)Γ((e t e ⊤ t F ⊤ e i ) ⊗ I p )H = n i=1 T t=1 (I T ⊗ (e ⊤ t F ⊤ e i e ⊤ i X))Γ(e t ⊗ I p )H = T t=1 (I T ⊗ (e ⊤ t F ⊤ X))Γ(e t ⊗ I p )H.We have by the triangle inequality,∥ Υ 1 ∥ op ≤ T ∥Γ∥ op ∥X∥ op ∥F ∥ F ∥H∥ F .By the triangle inequality and the upper bound of ∥Υ 2 ∥ op , we have∥Υ 4 ∥ op ≤ ∥Υ 2 ∥ op ∥X∥ op ∥H∥ F + ∥ Υ 1 ∥ op ≤ T ∥Γ∥ op ∥X∥ 2 op ∥H∥ 2 F + T ∥Γ∥ op ∥X∥ op ∥F ∥ F ∥H∥ F .Squaring both sides and taking conditional expectation, we haveE[∥Υ 4 ∥ 2 op | ε] ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ),thanks to the upper bound of ∥Γ∥ op in Lemma D.4 and the moment bounds of ∥X∥ op , ∥F ∥ F , ∥H∥ F in Lemma D.1.Bound of Υ 5 . By the expression for the matrix Υ 5 ∈ R T ×T obtained in Appendix D.4, we haveΥ 5 = Υ ⊤ 1 X ⊤ F -Υ 2 ,whereΥ 2 = - T t=1 F ⊤ D t XHe t + n l=1 (e ⊤ l XH ⊗ F ⊤ )DS(I T ⊗ X)Γ ⊤ (I T ⊗ X ⊤ )D(I T ⊗ e l ).By the triangle inequality,∥ Υ 2 ∥ op ≤ T ∥X∥ op ∥F ∥ F ∥H∥ F + ∥X∥ 3 op ∥Γ∥ op ∥F ∥ F ∥H∥ F .Using the moment bounds of ∥Υ 1 ∥ op , ∥X∥ op , ∥F ∥ F , ∥ F ∥ F , ∥H∥ F , we haveE[∥Υ 5 ∥ 2 op | ε] ≤ n 2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).We first state three useful lemmas.Lemma D.6 (Adopted from Lemma E.10 of [26] ). Let U , V : R n×p → R n×T be two locally Lipschitz functions of Z with i.i.d. N (0, 1) entries, thenE U ⊤ ZV - p j=1 n i=1 ∂ ∂z ij U ⊤ e i e ⊤ j V 2 F ≤ E∥U ∥ 2 F ∥V ∥ 2 F + E ij 2∥V ∥ 2 F ∂U ∂z ij 2 F + 2∥U ∥ 2 F ∂V ∂z ij 2 F .Lemma D.7 (Adopted from Lemma F.5 of [5] ). Let U , V : R n×p → R n×T be two locally Lipschitz functions of Z with i.i.d. N(0, 1) entries. Provided the following expectations are finite, we haveE pU ⊤ V - p j=1 n i=1 ∂ ij U ⊤ e i -U ⊤ Ze j n i=1 ∂ ij e ⊤ i V -e ⊤ j Z ⊤ V F ≤ (1 + 2 √ p) E[∥U ∥ 4 F ] 1/2 + E[∥V ∥ 4 F ] 1/2 + E[∥U ∥ 4 ∂ ] 1/2 + E[∥V ∥ 4 ∂ ] 1/2 , where ∂ ij U = ∂U /∂z ij and ∥U ∥ ∂ = ( n i=1 p j=1 ∥∂ ij U ∥ 2 F ) 1/2 .We will use the above two lemmas, conditionally on ε, to bound the conditional moments of Θ 1 , Θ 2 , Θ 3 , Θ 4 and Θ 5 given ε.Bound of Θ 1 . By the definition of Θ 1 , we haveEQUATION= Θ 1 + Υ 3 .Applying Lemma D.6 conditionally on ε to (Z, U , V ) = (X, F , H) givesE[∥Θ 1 ∥ 2 F | ε] ≲ E F ⊤ XH - i,j ∂F ⊤ e i e ⊤ j H ∂x ij 2 | ε + E[∥Υ 3 ∥ 2 F | ε] ≲ E[∥F ∥ 2 F ∥H∥ 2 F | ε] + E ij ∥H∥ 2 F ∥ ∂F ∂x ij ∥ 2 F + ∥F ∥ 2 F ∥ ∂H ∂x ij ∥ 2 F | ε + E[∥Υ 3 ∥ 2 F | ε] ≤ nC(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ),where the last line uses the moment bounds of ∥F ∥ F , ∥H∥ F in Lemma D.1, the moment bounds of ∥ ∂H ∂xij ∥ F , ∥ ∂F ∂xij ∥ F in Lemma D.5, the moment bound of ∥Υ 3 ∥ op in Lemma B.5.Bound of Θ 2 . By the definition of Θ 2 , we have by ( 21)F ⊤ XX ⊤ F - i,j ∂F ⊤ e i e ⊤ j X ⊤ F ∂x ij = F ⊤ XX ⊤ F + KH ⊤ X ⊤ F -pF ⊤ F + F ⊤ F A ⊤ + Υ 5 = nΘ 2 + Υ 5 .Applying Lemma D.6 conditionally on ε to (Z, U , V ) = (X, F , X ⊤ F ) givesn 2 E[∥Θ 2 ∥ 2 F | ε] ≲ E F ⊤ XX ⊤ F - i,j ∂F ⊤ e i e ⊤ j X ⊤ F ∂x ij 2 | ε + E[∥Υ 5 ∥ 2 F | ε] ≲ E[∥F ∥ 2 F ∥X ⊤ F ∥ 2 F | ε] + E ij ∥X ⊤ F ∥ 2 F ∥ ∂F ∂x ij ∥ 2 F + ∥F ∥ 2 F ∥ ∂X ⊤ F ∂x ij ∥ 2 F | ε + E[∥Υ 5 ∥ 2 F | ε] ≲ E[∥ F ∥ 4 F ∥X∥ 2 op | ε] + E (1 + ∥X∥ 2 op )∥ F ∥ 2 F ij ∥ ∂ F ∂x ij ∥ 2 F | ε + T E[∥Υ 5 ∥ 2 op | ε] ≤ n 3 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Here, the penultimate line uses∥F ∥ F ≤ ∥ F ∥ F , ∥ ∂F ∂xij ∥ F ≤ ∥ ∂ F ∂xij ∥ F , and ∥ ∂X ⊤ F ∂xij ∥ F = ∥e ⊤ j e ⊤ i ∂ F ∂xij + X ⊤ ∂ F ∂xij ∥ F ≤ ∥ ∂ F ∂xij ∥ F + ∥X∥ op ∥ ∂ F ∂xij ∥ F .The last line uses the moment bounds of ∥X∥ op , ∥ F ∥ F in Lemma D.1, the moment bound of ∥ ∂ F ∂xij ∥ F in Lemma D.5, the moment bound of ∥Υ 5 ∥ op in Lemma B.5.Bound of Θ 3 . By the definition of Θ 3 , we haveH ⊤ X ⊤ XH - i,j ∂H ⊤ X ⊤ e i e ⊤ j H ∂x ij = H ⊤ X ⊤ XH -(nI T -A)H ⊤ H -H ⊤ X ⊤ F W ⊤ + Υ 4 by (21) = Θ 3 + Υ 4 .Applying Lemma D.6 conditionally on ε to (Z, U , V ) = (X, XH, H) givesE[∥Θ 3 ∥ 2 F | ε] ≲ E H ⊤ X ⊤ XH - i,j ∂H ⊤ X ⊤ e i e ⊤ j H ∂x ij 2 | ε + E[∥Υ 4 ∥ 2 F | ε] ≲ E[∥XH∥ 2 F ∥H∥ 2 F | ε] + E ij ∥XH∥ 2 F ∥ ∂H ∂x ij ∥ 2 F + ∥H∥ 2 F ∥ ∂XH ∂x ij ∥ 2 F | ε + E[∥Υ 4 ∥ 2 F | ε] ≲ E[∥X∥ 2 op ∥H∥ 4 F | ε] + E (1 + ∥X∥ 2 op )∥H∥ 2 F ij ∥ ∂H ∂x ij ∥ 2 F | ε + T E[∥Υ 4 ∥ 2 op | ε] ≤nC(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Here, the penultimate line uses ∥ ∂XH ∂xij ∥ F = ∥e i e ⊤ j ∂H ∂xij +X ∂H ∂xij ∥ F ≤ (1+∥X∥ op )∥ ∂H ∂xij ∥ F ,∂F ⊤ e i ∂x ij -F ⊤ Xe j n i=1 ∂ F ⊤ e i ∂x ij -F ⊤ Xe j ⊤ = ( KH ⊤ + F ⊤ X + Υ ⊤ 1 )( KH ⊤ + F ⊤ X + Υ ⊤ 1 ). This further implies pF ⊤ F - p j=1 n i=1 ∂F ⊤ e i ∂x ij -F ⊤ Xe j n i=1 ∂ F ⊤ e i ∂x ij -F ⊤ Xe j ⊤ = nΘ 4 + Υ ⊤ 1 ( KH ⊤ + F ⊤ X) + ( KH ⊤ + F ⊤ X) Υ ⊤ 1 + Υ ⊤ 1 Υ ⊤ 1 Υ4 . Applying Lemma D.7 conditionally on ε to (Z, U , V ) = (X, F , F ), nE[∥Θ 4 ∥ F | ε] ≲ E pF ⊤ F - p j=1 n i=1 ∂F ⊤ e i ∂x ij -F ⊤ Xe j n i=1 ∂ F ⊤ e i ∂x ij -F ⊤ Xe j ⊤ F | ε + E[∥ Υ 4 ∥ F | ε] ≲ (1 + 2 √ p) E[∥F ∥ 4 F | ε] 1/2 + E[∥ F ∥ 4 F | ε] 1/2 + E[∥F ∥ 4 ∂ | ε] 1/2 + E[∥ F ∥ 4 ∂ | ε] 1/2 + T E[∥ Υ 4 ∥ 2 op | ε] ≤ n 3/2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Here, the last inequality uses the moment bounds of ∥F ∥ F , ∥ F ∥ F in Lemma D. ∂H ⊤ e j ∂x ij -H ⊤ X ⊤ e i p j=1 ∂H ⊤ e j ∂x ij -H ⊤ X ⊤ e i ⊤ = (W F ⊤ -H ⊤ X ⊤ -Υ ⊤ 2 )(W F ⊤ -H ⊤ X ⊤ -Υ ⊤ 2 ) ⊤ . By the definition of Θ 5 , we have nH ⊤ H - n i=1 p j=1 ∂H ⊤ e j ∂x ij -H ⊤ X ⊤ e i p j=1 ∂H ⊤ e j ∂x ij -H ⊤ X ⊤ e i ⊤ = Θ 5 + Υ ⊤ 2 (W F ⊤ -H ⊤ X ⊤ ) ⊤ + (W F ⊤ -H ⊤ X ⊤ )Υ 2 -Υ ⊤ 2 Υ 2 Υ5.HereΥ 2 = j ((e ⊤ j H) ⊗ I n )DS(I T ⊗ X)Γ ⊤ (I T ⊗ e j ).Applying Lemma D.7 conditionally on ε to (Z, U , V ) = (X ⊤ , H, H) (i.e., consider the mapping from R p×n to R p×T :X ⊤ → H) gives E[∥Θ 5 ∥ F | ε] ≤ E nH ⊤ H - n i=1 p j=1 ∂H ⊤ e j ∂x ij -H ⊤ X ⊤ e i p j=1 ∂H ⊤ e j ∂x ij -H ⊤ X ⊤ e i ⊤ F | ε + E[∥ Υ 5 ∥ F | ε] ≲ (1 + 2 √ p) E[∥H∥ 4 F | ε] 1/2 + E[∥H∥ 4 ∂ | ε] 1/2 + T E[∥ Υ 5 ∥ 2 op | ε] ≤ n 1/2 C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Here, the last line use the moment bounds of ∥H∥ F in Lemma D.1, the moment bound of ∥H∥ 2 ∂ in Lemma D.5, and the moment bound of ∥ Υ 5 ∥ op in Lemma B.5.Bound of Θ 6 . Using (19) , we haven i=1 p j=1 ∂E ⊤ e i e ⊤ j H ∂x ij = E ⊤ (F W ⊤ -Υ 2 ).It follows thatE ⊤ XH - n i=1 p j=1 ∂E ⊤ e i e ⊤ j H ∂x ij = E ⊤ XH -E ⊤ (F W ⊤ -Υ 2 ) = ∥E∥ F Θ 6 + E ⊤ Υ 2 .Thus, using E = E/∥E∥ F , we haveΘ 6 = E ⊤ XH - n i=1 p j=1 ∂ E ⊤ e i e ⊤ j H ∂x ij -E ⊤ Υ 2 . Applying Lemma D.6 conditionally on ε to (Z, U , V ) = (X, E, H)gives E[∥Θ 6 ∥ 2 F | ε] ≲ E E ⊤ XH - n i=1 p j=1 ∂ E ⊤ e i e ⊤ j H ∂x ij 2 | ε + E[∥ E ⊤ Υ 2 ∥ 2 F | ε] ≲ E[∥ E∥ 2 F ∥H∥ 2 F | ε] + E ∥ E∥ 2 F ∥ ∂H ∂x ij ∥ 2 F | ε + E[∥ E ⊤ Υ 2 ∥ 2 F | ε] ≤ E[∥H∥ 2 F | ε] + E ∥ ∂H ∂x ij ∥ 2 F | ε + E[∥Υ 2 ∥ 2 op | ε] ≤ C(T, γ, η max , c 0 )(δ 2 + ∥b * ∥ 2 ).Here, the last line uses the moment bound of ∥H∥ F in Lemma D.1, the moment bounds of ∥ ∂H ∂xij ∥ F in Lemma D.5, and the moment bound of ∥Υ 2 ∥ op in Lemma B.5. This finishes the proof of Lemma B.6. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: Given the theoretical nature of our work, it does not have any negative societal impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: The paper is of the theoretical nature, it does not invent any new models or datasets.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [NA] Justification: We write all the code by ourselves, and did not use other code or data.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.• If this information is not available online, the authors are encouraged to reach out to the asset's creators. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: The paper introduces estimators that accurately track the generation error of iterates produced by gradient descent and SGD (and their proximal variants) on high-dimensional (proportional regime) robust regression problems.  The proposed estimators are very similar to the previous work of Bellec and Tan ’24 -- a distinction between this work and that of Bellec and Tan should be clarified.

Reviews: review: paper_summary: paper_summary: #### Summary
This paper examines the generalization performance of iterates produced by Gradient Descent (GD), Stochastic Gradient Descent (SGD), and their proximal variants in high-dimensional robust regression problems. The paper introduces estimators that accurately track the generalization error of the iterates along the trajectory of the iterative algorithms. These estimators are proven to be consistent under certain conditions and are illustrated through several examples, including Huber regression, pseudo-Huber regression, and their penalized variants with non-smooth regularizers.
strength_weakness: strength_weakness: #### Weaknesses
1. **Computational Complexity**: The proposed estimators involve complex calculations that may be computationally intensive, especially for large datasets. More discussion on computational efficiency and scalability is needed.
2. **Generality**: The paper focuses on specific types of robust regression problems. Extending the methodology to a broader range of regression problems would increase its impact.
3. **Comparison with Existing Methods**: The paper provides limited empirical comparisons with existing state-of-the-art methods for estimating generalization error. More comparative analysis would strengthen the validity of the proposed approach.
4. **Practical Guidelines**: While the theoretical results are robust, practical guidelines for implementing the estimators in real-world scenarios are not sufficiently detailed.
questions: questions: #### Questions
1. **Computational Complexity**:
    - Could you provide more details on the computational complexity of the proposed estimators? How do they scale with increasing dataset size and dimensionality?

2. **Generality**:
    - The paper focuses on specific types of robust regression problems. Are there any challenges in extending the methodology to other types of regression problems, such as those with different types of regularizers or loss functions?

3. **Comparison with Existing Methods**:
    - How do the proposed estimators compare empirically with existing methods for estimating generalization error? Are there scenarios where your approach significantly outperforms others?

4. **Practical Implementation**:
    - Can you provide practical guidelines or heuristics for implementing the proposed estimators in real-world scenarios? What are the key considerations practitioners should keep in mind?
limitations: limitations: no
paper_summary: paper_summary: This manuscript aims at finding computational efficient measure of generalization performance for high-dimensional robust regression with regularization. In this scenario, when loss function is not quadratic, estimating the out-of-sample error $\| \Sigma^{1/2}(b_t - b^* ) \|^2$ (where $\Sigma$ is the covariance matrix, $b_t$ is the $t$-iteration variable and $b^*$ is the ground truth) for iterative methods like GD, SGD or proximal GD can be very difficult and is not well-explored in previous literature. To overcome this issue, the authors propose a new estimator for the out-of-sample error for a general class of gradient methods and show that the difference converges to zero in probability when sample size goes to infinity.
strength_weakness: strength_weakness: Despite being more general, this paper seems to constitute *limited* progress towards the community. When comparing it with the work by Bellec and Tan, it is not hard to notice the estimator $\hat{r}_t$ in Theorem 3.6 is a trivial extension of Bellec and Tan's estimator from square loss to convex loss. The exact form of the proposed estimator and its guarantee in Theorem 3.6 are very similar to the result of the mentioned paper. Also, the proof looks quite straighforward and extending from square loss to convex loss does not rely on technical innovations. In this regard, I believe this paper falls short in getting accepted to NeurIPS, unless the authors could provide more convincing arguments how this result is *fundamentally* different from previous works.
questions: questions: There are no further questions.
limitations: limitations: This work is purely theoretical and has no negative influences.
paper_summary: paper_summary: This paper focuses on deriving consistent estimators of the generalization error for robust regression with Gaussian design and heavy tailed noise along the SGD trajectory. They propose two estimators, $\hat r$ (which requires knowledge of the covariance) and $\tilde r$ (which does not) and prove consistency of both estimators when the number of samples is at least linear in the dimension. They support this with experiments which demonstrate that both $\hat r$ and $\tilde r$ accurately capture the generalization error in various settings.
strength_weakness: strength_weakness: - The numerators in Theorems 3.6 and 3.7 include unspecified constants $C(T,\ldots)$. From the appendix, it appears that $C(T,\ldots)$ scales like $T^T$, so that you can only run for $T \approx \log(n)/\log\log(n)$ steps before the bound is vacuous, which is severely restricting. I'm concerned it may not be possible for SGD to significantly decrease the loss in this number of steps, rendering the generalization estimate vacuous.

- While some intuition for $\boldsymbol{W}$ (eq 8) is given in Section 3.1, none is given for $\widehat{\mathbf{A}}$ (eq 9) or $\widehat{\mathbf{K}}$ except that they are used to construct the weights for Theorem 3.7. This makes the notation in sections 3.2 and 3.3 somewhat difficult to follow.

Minor Points:
- line 29: starts -> start
- line 118: maybe provide some intuition here for $\phi_t, \psi$ rather than deferring it to sections 2.1/2.2 (e.g. just $\phi_t$ is a proximal operator)
- line 170: all -> both
- line 207: it is strange to define $\hat r$ in terms of $\mathbf{W}$ and $\tilde r$ in terms of $\widehat{\mathbf{W}}$. Perhaps use $\widetilde{\mathbf{W}}$ instead of $\widehat{\mathbf{W}}$?
- line 269: the three ... reveal
questions: questions: - What are the dependencies on $T,\eta_{max}$ in the numerators of Theorems 2.6 and 2.7?
- Is the reason that $W$ must be computed recursively simply to invert $\mathcal{M}$ block by block? I don't see how to connect the unrolling in section 3.1 to the Kronecker expressions in 3.2.
limitations: limitations: The authors have adequately addressed the limitations of their work.
paper_summary: paper_summary: The paper explores linear regression with Gaussian design corrupted by noise that has bounded first moments, focusing on the 'high-dimensional' regime where the feature dimensionality $p$ scales proportionally with the number of training data $n$. The authors propose two estimators that, aside from an additive noise term, approximate the out-of-sample error (measured with the square loss) achieved by a general iterative scheme. This scheme includes commonly used solvers like GD, SGD, proximal GD, and proximal SGD, which are used to compute the regularized empirical risk minimizer with a convex, differentiable, and smooth loss function (with Lipschitz-continuous gradients) and a potentially non-smooth regularizer. This framework encompasses robust linear regression with heavy-tailed noise using the Huber loss function. The authors provide non-asymptotic guarantees in probability for finite $n$ (and $p$), demonstrating the asymptotic consistency of their estimators (Theorems 3.6 and 3.7). They also present examples and numerical experiments. Their primary tool is a generalized probabilistic approximation of the generalization error previously discussed in specific settings in the references [3] and [5].
strength_weakness: strength_weakness: The main results of the paper can be seen as an extension of the findings in [5], applied to a broader context that encompasses the stochastic case and extends beyond the square loss. While the generality of the weaker conditions is appreciated and introduces a different structure, the work is somewhat marginal in its contribution. Although the new methodology can track the generalization performance over training time for various iterative solver schemes of interest, it remains uncertain how this methodology can be used to (provably) "optimally" tune parameters other than the stopping time. This limitation seems notable, as the increased flexibility of the algorithmic schemes—allowing for additional tuning parameters such as mini-batch size and loss-specific parameters—is not accompanied by an enhanced analysis of these tuning parameters.
questions: questions: 1) In the case of the square loss and full-batch gradient descent, do the authors recover exactly the same theoretical results in [5]? If not, what are the differences in this case?

2) The authors claim that their procedure can be used to choose the "tuning parameters that achieves the smallest out-of-sample error" (Page 2). While this is clear for stopping time, can the authors give examples of how their methodology can be use to tune other tuning parameters?

3) While the authors clearly highlight the relevance of their contributions to robust statistics and their goal of tracking solver performance over iteration time, a specific comparison of their methodology and consistency results to other classical methods in some fundamental models, such as ridge regression and lasso regression, would be beneficial. Although Section 1.1 provides a discussion on this topic, it lacks quantitative details, making it unclear how the proposed methods compare in terms of non-asymptotic statistical rates (related to consistency results) and computational efficiency. For example, the computational complexity of the data-dependent estimator the authors consider ($\tilde{r}_t$ on page 6) is not clearly compared to other estimators in basic settings. Can the authors address this?

4) The paper does not seem to make any assumptions on the penalty term, i.e. the function g in (2). This seems surprising. Can the authors elaborate on this?

5) Could the authors comment on the challenges of extending their analysis beyond the proportional regime they consider, where the data dimension is of the same order as the sample size?

****
I have increased my score to 'Accept' post rebuttal.
limitations: limitations: The paper discusses limitation of the current approach as a way to motivate follow up research in the conclusions.

Comments: comment: Thanks for reading the rebuttal and taking the time to re-read the proof. One challenge that arose when attempting to leverage out-of-batch samples, is that we are now manipulating quantities that are not involved in the iterative algorithm or its derivatives. This leads to some difficulty in obtaining bounds such as (35) around line 444 where the approximation between $\hat A$ and $\hat KW$ (involving out-of-bag samples for invertibility of $\hat K$) holds in the row space of $F$ (the stochastic gradient matrix, that does not involve out of bag samples). What we would like to point out as a final remark with this, is that it is subjective to argue about the challenges/difficulty of a proof once the final product is finished, since the final product does not showcase the challenges and difficulty necessary to produce it.

In any case, many thanks for your comments and work on this manuscript. We believe that case of (proximal) SGD with robust losses is important for the NeurIPS readership and community. We will emphasize the role and necessity to use out-of-bag samples, and how this departs from [5], in the camera-ready version (should the paper be accepted).
Thanks for your final comments and your work throughout the process.
Thanks for reading the rebuttal and the additional comment. We agree to include, in the camera-ready version, a discussion of C(T), the disconnect between the theory and simulations somewhere below Theorem 3.6, as well as the pointers to the literature on recent advances to lower dependence on T (including Li and Wei, 2024).
I thank the reviewers for addressing my concerns. The detailed explanations regarding the technical difficulties and the empirical evidence provided are very helpful. I would like to slightly increase my score, however, I am still very dubious about the contribution and novelty of this work. I took some time to re-read the proof and response of the authors. While I acknowledge that, directly generalizing [5] to SGD and non-square loss is not possible, I still question if this technique is truly innovative. The invertibility of $\tilde K$ seems to be solved by easy tricks and the rest proof still looks very similar to [5]. As a result, I change my score to 4.
I thank the authors for addressing my questions and concerns, but I prefer to keep my current score. I am skeptical of the “revealing the form of the estimator” claim as an estimator that works for small T at might require additional correction terms to work for large T. I understanding fixing this issue may be challenging and can be reasonably left to future work, but it would be helpful to include a discussion of C(T) and the disconnect between the theory and the experiments somewhere below Theorem 3.6.
I appreciate the authors' response and maintain a positive outlook about this work. As a result of the authors' reply, I am raising my evaluation to Accept.