Estimating Generalization Performance Along the Trajectory of Proximal SGD in Robust Regression
Authors: Kai Tan
Pierre C Bellec
Venue: NeurIPS 2024 poster (2024)
Keywords: Robust regression
generalization error
stochastic gradient descent
early stopping
Stein's formula
Primary area: learning_theory
TL;DR: This paper provides generalization error estimates for each iteration of the SGD and proximal SGD algorithms in the context of robust regression.

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.

Estimating Generalization Performance Along the Trajectory of Proximal SGD in Robust Regression
Kai Tan Department of Statistics Rutgers University Piscataway, NJ 08854
Pierre C. Bellec
Department of Statistics Rutgers University
Piscataway, NJ 08854
kai, tangrutgers..... edu
www.zeebbllee@@uuterrr.ccc
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 compa- rable to the sample size and errors may be heavy-tailed. We introduce estimators that precisely track the generalization error the iterates along the trajectory the iterative algorithm. These estimators provably consistent under suitable conditions. The results are illustrated through several examples, including Huber regression, ikeudo-Huber regression, and their penalized variants with non-smooth regularizer. We provide explicit generalization error estimates for iterates gener- ated from GD and SGD, from proximal SGD in the presence of non-smooth regularizer. The proposed risk estimates serve as effective proxies for the actual generalization error, allowing us determine the optimal stopping iteration that minimizes the generalization error. Extensive simulations confirm the effectiveness of the proposed generalization error estimates.
Introduction
Consider the linear model:
y Xb" +
(1)
where y R is the response vector, X E Ruxp is the design matrix, b" RRP is the unknown regression vector, and e E IR" the noise vector that we assume independent of X. The entries 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 E arg min BER
Σ p(u x..b)+((b)) g(b),
YI
where P R is data-fitting loss and g TRP IR is regularization function. In the present robust regression setting, typical examples of include the Huber 14] loss p(r; =
82 f min(1, 2)) dz, the Pseudo-Huber loss p(r: 8) (v1++ (r/5)2 1) or other Lipschitz
loss functions to combat the possible heavy-tails of the additive noise. Typical examples of penalty
functions include the Li/Lasso 27 penalty g(b) A|| group-Lasso penalty (29 for grouped
variables, or their non-convex variants including for instance SCAD [2] or MCP 30
In order solve the optimization problem practitioners resort iterative algorithms, for instance gradient descent, accelerated gradient descent, stochastic gradient descent, and the corresponding proximal methods (20) in the presence of non-smooth regularizer. Let the algorithm starts with
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

some initializer b' E R' (typically 8' = 0) followed by consecutive iterates D² be.b..... where b' is typically obtained, for gradient descent and its variants as will be detailed below, from D 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'.
We assume throughout that the covariance [ E of the feature vectors is finite. We measure the predictive performance of b' using the out-of-sample error
E ( TT -i b") (zi, 3/ setml]
12 b²)]l²
where Fnew is new feature vector, independent of the data (æ:; 1/4 jie(n) and has the same distribution as z;. The above squared metric is used because the noise En (and thus yı) is allowed to have infinite
variance, and in this case the squared prediction error [[[ new b' irrespective of the value of B'.
Vnew) )² (z), 34 ve]] +oo
The paper proposes to estimate the out-of-sample error /2/2(b' b*)ll of the t-th iterate using the right-hand side of the approximation
/2 -
[ (I'xa') +
Σ We wr. 1b x60)|//m
(3)
where vb Rd + R is the derivative of acting component-wsse on each coordinate in IR" and S. is a diagonal matrix of the form S = Eier, e,e; where 1. C [n]] is the batch for the s-th stochastic gradient update and e; E _n is the i-th canonical basis vector. Here the Wt,s are quantities, introduced in Section 3.3 below, that can be computed from data and do not require the knowledge of E. 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 in probability) for first set of weights (Wa,))<<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*) are computed from the iterative algorithm, the approximation (3) lets us compare the oft-llesample error of iterates
b' at different time up to the additive term || /n (which does not depend on 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 oft-of-sample error.
1.1 Related literature
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 cross- validation or leave-one-out schemes. In the proportional regime of interest here, where dimension p and sample size 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-ord 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-spliting including Approximate on sample-spliting and refitting; see [I 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 avoids data-splitting and refitting; it only needs to fit the full data once and then adjust the training error by multiplicative factor larger than 1. Beyond ridge regression, the extension of GCV using degrees-of of-freedomo has been studied for Lasso regression 2 18 and alternatives were developed for robust M-estimators [3 4). While ALO 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
(ALO) schemes (4] that do not rely
r' obtained by widely-used iterative algorithms such gradient descent (GD), stochastic gradient

---

descent (SGD) or their proximal variants: ALO and GCV focus on estimating the final (t + +00)
iterate of the algorithm, when solution b in 2 is found. Our goal in the present paper is develop risk estimation methodologies along the trajectory of the algorithm.
Luo et al. 17 developed methods to estimate the vross-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.c., the smallest and largest cigenvalues This condition is not satisfied for the regression problems we consider in this paper, such as high- dimensional 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 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, 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 [I3] 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.
bounded) along all iterates.
Most recently, and most closely related to the present paper, [5] proposed risk estimate for iterates 8 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 p 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 [S) we will explain in Remark (3.8 that directly generalizing (5) to SGD in robust regression leads to 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. (2 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 [T0 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 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 random batch is used for each update.
2 Problem setup
The paper studies iterative algorithms aimed solving the optimization problem 2 We consider the algorithm that generates iterates b' for t = 1,2,..., ..., T according to the following iteration:
6 Φ (b.' 6' + M) x Sub(y -x6')
(4)
[1
where S, E R is the diagonal matrix S, = Eigr, e,e.T for Io C [n] the f-th batch (independent of (X, y)),, where Φb RP + Rp and r/6 12 + R are two functions and t/. is the step size. Typically, sb : Th - RD is the componentwise application of p' (where p is the data-fitting loss in (2 and the matrix S, E gax= is diagonal with elements in {0, 1} encoding the observations i e [n]
3

---

used in the t-th stochastic gradient update. The presence of S and possibly nonlinear function rbs such that the above iteration scheme is not covered by previous related works including [5] which only tackle S; = I, (full batch gradient updates) and eb : R - IR" the identity map (p in 2 restricted to be the square loss). The iterative scheme on the other hand, covers SGD with robust loss functions.
In the next section, we first provide few examples of algorithms encompassed in the general iteration 44. 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 with non-smooth regularizers.
problem 22
2.1 Robust regression without penalty
If there are no penalties in 2 i.e., g(b) O, then the minimization problem becomes
1 Σ p(y: x. b).
b E arg min DE
To solve this problem, provided p 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:
8 = 8 The
bb - Xi'),
(5)
where vb is the derivative of acting component-wisg wise on its argument, and Tu is the step size (also known as learning rate). For the least squares loss p(z) _ /2, we have (b(u) u.
Example 2.2 (SGD). Suppose at f-th iteration, we use the batch Ie C [ni] to compute the gradient,
-'
T/t L
Σ iEL
aublyn
x; bD) = B' + The X
Sab(y X6%)
4
where S: = EiEl, eje, and e, is the i-th canonical vector in R". If 1, = [n] for each t, then [4l = n and S = In. hence this SGD method reduces to the GD method in
2.2 Robust regression with Lasso penalty
Regularized regression useful for high-dimensional regression problems where is larger than n. We consider the Lasso penalty g( x||b|l 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* sparse). While GD and SGD are not directly applicable to solve the optimization problem (2) with Lasso penalty due to || la lacking differentiability at O, 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) = A|| in 2). the Proximal GD gives the following iterations:
7++1 = soft (b* R X 1 Xb*)),
where s-fta(-) applies the soft-thresholding u + sign (u) (Iu) - 6)+ component-wise.
Example 2.4 (Proximal SGD). Similar to the Proximal GD, the Proximal SGD consists of the following iterations:
1 soft (6' + 11
X Sab(y - X(')).
Let p' R Rd be the function applies the derivative of : R to each of its component, i.e., p'(u) = (p) p...,,'( p'(ux)) Then the above examples can be summarized in the following table with different definition of v/ di, and S,

---

Table 1: Specification of Φ. Or. S₁ for each algorithm
GD
SGD
Proximal GD
Proximal SGD
(w)
('(u)
p'(u)
p'(u)
p'(u) soft (v) S.
(v)
V
v
softxm, (v)
S
In
S
In
To define the proposed estimators of the generalization error, further define the following Jacobian matrices:
D. = (v(())
(0.(())
ro XO E R B
D. = dv v==
(y (x E RPP..
Ör
Th. [1.
x
Then, we have Dj I, for GD and SGD, and D, ==EEE, ege] for Proximal GD and Proximal SGD based soft-thresholding, where S. (j [y] e b #0}.
3 Main results
Assumption 3.1. The design matrix X has i.i.d. rows from (0, Σ) for some positive definite matrix
E satisfying < Amin (E) S
Amax((2) and | |E S k. We assume Var(ær b*] S 82,
that is, the signal of the model [ is bounded from above.
Assumption 3.2. The noise E is independent of X and has i.i.d. entries from a fixed distribution independent of n, p. with [[[e>|] S 6, that is, bounded first moment.
Assumption 3.3. The data fitting loss p R R is convex, continuously differentiable and its derivative sa I-Lipschitz and [+<[=]] S for all R. The function Φ; is I-Lipschitz and satisfies (0) 0. The matrices S. Eier, e,e,". and Le 2 con for some positive constant co (0. 1]. Let Mmax = max Th--
Huber loss and Dsuedo-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 and feature dimension p satisfy p/n S 7 for constant E (0, 20).
3.1 Intuition regarding the estimates of the generalization error
This subsection provides the intuition behind the definition of the estimates define below. For the sake of clarify, and in this subsection only, assume that
Σ == Ip:
=
7// 1/n.
With the above working assumptions, the validity of the estimates defined below relies on the probabilistic approximation
P
a
2
1²1l2
x ( b²)
dz.;
which was developed in 3 for risk estimation purposes, but outside the context of iterative algorithms. Above, ej E IP' 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 the residual y; x! 6' 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' with respect to the entries (Zi)))e(=),jEfr) of X. In (4). each iterate is relatively simple function of the previous ones, with the simplifications 7 this is b++1 (b' X S Xb')/n) For I 1. given that D a constant initialization, azu;
a 0_ = Die,e, Siv(y X6*) D S (b, - bb);/n.
We find in the proof, that when summing these quantities over j E [p). the second term involving
5

---

(6' -b); is negligible, and the same negligibility holds at later iterations with terms involving
b6' bb);
) (or any b") s S t). By performing a similar simple calculation at the next iteration,
and ignoring these terms, we find with ft = e, Suabyy Xb') and s? = e, S2bbyy x by the chain rule
ab dzi
D1 T2
P
ab Ozy
D2 71
TF
It D2 x S,D2x/ f.
22211
Maa,
This reveals the weights (wa,t) sect in 3 iteration & 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., D, and the stochastic gradient matrix S,, becomes increasingly complex. This recursive unrolling of the derivatives can be performed numerically at the same time 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.
3.2 Formal matrix notation to capture recursive derivatives
We now set up the matrix notation that captures this recursive unrolling of the derivatives by the chain rule. Throughout, T is the final number of iterations. Define three block diagonal matrices D E
RaTxnI D E RTTx and S paIxnt by D = Ei=1((exe ) D.). D Σ ((ere) )&D:). D.). and S = =I= (((ee)) @ St). Now we are ready to introduce the following matrices of size x T:
W = Σ (Ir 8 ( ( "17) ee),,
(8)
j=1
A = x
(9)
d=
K Σ
Tree - Σ (Ir e, DD @ X)T(Ir & X SDD(D(I 0e),
(10)
where r = M- L(A & 1,) e ReTTpT L 2(eel 1 = Et= e
I P
I
M
and
P. P. D.
1. X S.D.X)..
P I
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
3.3 Main results: estimating the generalization error consistently
For each iterate b', define the target ra (generalization error) and its estimate fz by
||e| n
I/2 - 2²)
(y- Xb
Σ Wt. Wr. S.y - - (1))
where Wt,s e, We, and W E RXX is the matrix defined in 8 The following shows that | rel 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 Ve
P P(1re rel >))
max
C(T, Tmaxx Co, 8, x)
E min {1, |Iel| [1)
(12)
holds then E1min{1, lell }] + O, so that, n, - too while (T 7, Thnxx; Co2 8, K., 6) are held fixed, the right-hand side converges and it Tt converges to in probability.
If additionally Assumption 3.2

---

This establishes that Pe is consistent at estimating rg. The statement E6 E[min ||=|l//]]]]] + 0 is equivalent ||=l² /n2 _ (convergence in probability), and is proved in [2 under the assumption that Eley| < +o0 with e i.i.d. from a fixed distribution; this allows Var[e:] = +oo as long as the first moment finite. The expression of W involves Σ, which is typically unknown 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 K where K and A are the T x T matrices defined in 9
10
We define another estimate Ft by replacing W in with W = K-IA:
1-1
(y X0) + Σ
wr S.b.y X6*]
where wit.. = e. Wes.
Theorem 3.7 (Proved in Appendix C.3
Under Assumptions 3.1 and 3.3 to 3.5
for any >> 0,
P(rr - - rel re)
2e-nn/ n/18 + max C(T.N.hao..t.oaa)
vn + E[min(1, |I] [
s 2e
If additionally Assumption 3.2 holds then E1min(1, |ell + O, so that, as n,P - too while (T, ", Mmaxx Co. 5, K. c) are held fixed, the right-hand side converges 0 and Pt Tt converges to 0 in probability.
This establishes the consistency of ît. The simulations presented next confirm that the two estimates i and ft both accurate estimates of r.. The estimate ig has the advantage of not relying on the knowledge of E and are recommended in practice.
Remark 3.8. We highlight that directly generalizing the approach in 5 would lead to the approx- imation A KW, where A and K are given (23) and (26). respectively. From A x 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, Dj] in (26). where S, E {0, jaxn is the diagonal matrix with 1 in position (i, f) 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 (Dt) is 0 or small. Let Pquh denote the estimate of the same form as Fg but using the weight matrix K- A instead. Simulation results in Figure 4 confirm that 17 is suboptimal compared to our proposed n. For SGD and proximal SGD, we solved this issue regarding the invertibility of K by using out-of- samples in the construction of K and Å, in order to avoid S, in the diagonal elements of K in equation BII. 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, Y, Maaxx Cu, 5, k)) in the above results is not explicit. Inspection of the proof reveals that the dependence of this constant in T is currently 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 255 24 feature for instance the same log(n)/ log log 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 x (n)y(n) while still controlling certain AMP related quantities, although for the risk [15] equations (16)-(17) the condition required on T still logarithmic in n. This suggests that advances on this front are possible, least for isotropic design and specific loss and regularizer such as those studied in [13: 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 a 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 to cover a general framework allowing iterations of the form 4 with little restriction on the nonlinear functions except being Lipschitz.
7

---

4 Simulation
In this section, we present numerical experiments 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 / SCD- BoD-generliee milltiiiion | Our goal is to compare the performance of the proposed risk estimates with the true risk r, for different regression methods and iterative algorithms.
com kaitan365/
We generate the dataset (X, y) from the linear model that is, y Xb + E. Here, the rows of X E Raxy are sampled from centered multivariate normal distribution with covariance matrix E = Ip. The noise vector ε consists of i.i.d. entries drawn from a f distribution with two degrees of freedom so that the noise variance is infinite. The true regression vector 6* Rº is chosen with p/20 nonzero entries, set to a constant value such that the signal strength ||6*|| equals 10.
We explore two scenarios of the 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 noftttressholdin). Both the gradient descent (GD) and stochastic gradient descent (SGD) algorithms are implemented for each type of regression.
(ii) (n, (10000, 12000): Here, we investigate Huber regression and Pseudo-Hubes regres- sion with an LI penalty, A|[b[| (A = 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 D = Op and proceed with a fixed step size = (1 + / p/m where n, for GD and proximal GD, and n, n/5 for SGD and proximal SGD. We run each algorithm for = 100 steps. For SGD and Proximal SGD, batches 4 c {1,2,... ... n are randomly sampled without replacement and independently of (X, y, (1, eett), each with cardinality [ = -.
A crucial component of the proposed risk estimates i; and Fg involve the weight matrices W and W. The matrix W is defined in Theorem 3.6 and W --1A is defined in Theorem 3.7 We employ Hutchinson's '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 for Pt |///266' b²)l
||=|l//m, one is é in
Theorem 3.6 which requires knowing Σ = E[ J. and the other is F, in Theorem (37 which does |e|l2/1 /n remains constant along the algorithm trajectory, we only
not need Σ. Since the quantity ||e||²/
focus on the estimation of (/ b²) 12 We repeat each numerical experiment 100 times and
present the aggregated results in Figures
and 3
In Figure b²) and its two estimates A; ||=| // and r. ||=| /1 along with the 2 standard error bar for GD and SGD algorithms applied to both Huber and Pseudo Huber regression. In Figure [ 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 Ll-peallized Huber regression and Pseudo-Huber regression.
focus on the scenario with (n, p) (10000, 5000), and plot the actual risk /2
I and Figure confirm the three curves are in close agreement, indicating that the proposed estimates Å; ||=|/// /n and F; ||=||/ /n are consistent estimates of the actual risk 1/2 (6' b²) The two estimates closely capture the risk || |/// (6' b²)l[2 over the entire trajectory of the algorithms. For GD) and Proximal GD, the risk curves exhibit U-shape, first decreasing and then increasing, and the estimates f; and F. 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.
Figure
Additional experiments: varying step sizes for different iterations. We also conduct simulations to investigate the accuracy of the proposed risk estimates setting with varying step size. We consider two types of step sizes: 1). 1). 1 if is odd, and nt 0 if is even; 2). nt- 1 if is odd,
8

---

1 b" ?/1²
11 b -²
Fe f- |e|l/1
Ft - |e|B//n
Fr - |ee|//m
Fr - |ee|//n
Value
Value
200
Iteration number
Iteration number
(a) Risk curves of GD for Huber regression (b) Risk curves of SGD for Huber regression
1| 'II
Pt tell
Fe-11e1991
Fc [ee]
Value
Value
Iteration number
lteration number
(c) Risk curves of GD for Pscudo-Hulee re- (d) Risk curves of SGD for Pscudo-Huber gression
regression
Figure Risk curves for Huber and Pseudo-Huber regression with GD and SGD algorithms for the scenario (n. p) = (10000, 5000). Upper row: Huber regression, Lower row: Pseudo-Huber regression. Left column: GD, Right column: SGD.
and 77. 0.5 if is even. While the above choices of step size not preferred in practice, here the goal is show that the proposed risk estimates is able to accurately capture the dynamics the risk 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 flat when is even. The results are presented in Figure 3 illustrating that the risk estimates accurately capture the flat segments of the true risk curve.
Additional experiments: the estimate pauh suboptimal. We compare the performance of poub with our proposed estimates in Huber regression with (n, 7) (4000, 1000, 20) and batch size 0.2 for all / E [7]. It is clear from Figure that F, is more accurate than the
|/| = 1/10 and suboptimal estimator iuub, especially when increases.

---

11 - b²
Fc -- 1e|f/m
Ft - |el
F 11
Fr - |ee|//n
Value
Value
150
200
Iteration number
Iteration number
(a) Risk curves of proximal GD for Huber (b) Risk curves of proximal SGD for Huber regression
regression
|| - MI?
Pe- [e|]2m
i - [e]]
Fe- 1el|
||e|_
Value
Value
number
(c) Proximal GD for Pseudo-Huber regres- (d) Proximal SGD for Pseudo-Huber regres- sion
sion
Figure 2: Risk curves for LI-penalized Huber and Pseudo-Hute Huber regression with Proximal GD and Proximal SGD algorithms for the scenario (n, (10000, 12000). Upper row: LI Ll-penalized Huber regression, Lower row: LI-penalized Pscudo-Huber regression. Left column: Proximal GD, Right column: Proximal SGD.
5
Discussion
This paper proposes 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 I XT in (8) (8- 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 Φ;, and it would be of interest to study the extension in which our consistency results are universal, allowing non- con-Gaussian feature distributions. Second, it of interest to extend the current estimates more general optimization problems of the form (2) with non-smooth data-fitting loss, for instance the Least Absolute Deviation this case 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]. of independent interest derive the risk estimates for iterates obtained by such primal-dual methods.
loss yy Xb|11-
10

---

Acknowledgments
C. Bellec acknowledges partial support from the NSF Grants The authors acknowledge the Office of Advanced Research Computing (OARC) at Rutgers, The State University New Jersey for providing access to the Amarel cluster and associated research computing resources that have contributed the results reported here. The authors thank the anonymous reviewers for their valuable comments and suggestions that helped improve the presentation of the paper.
and DMS-2413679.
References
[1] Arna Auddy, Haolin Zou, Kamiar Rahnamarad, and Arian Maleki. Approximate leave-one-out cross validation for regression with 11 regularizers. In International Conference on Artificial Intelligence and Statistics, pages 2377-2385. PMLR, 2024.
[2] Mohsen Bayati, Murat A Erdogdu, and Andrea Montanari. Estimating lasso risk and noise level. In Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
[3] Pierre Bellec. Out-of-sample error estimate for robust m-estimators with convex penalty, Information and Inference: A Journal of the IMA, 12(4))2782 28177 10 2023.
[4] Pierre Bellec and Yiwei Shen. Derivatives and residual distribution of regularized m-estimators with application to adaptive tuning. In Conference Learning Theory, pages 1912 1947. PMLR, 2022.
[5] Pierre C Bellec and Kai Tan. Uncertainty quantification for iterative algorithms in linear models with application to early stopping. arXiv preprint arXiv:2404.17856 2024.
[6] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed opti- mization and statistical learning via the alternating direction method of multipliers. Foundations and Trends@ in Machine learning. 3(1):1-12)) 2011.
[7] Michael Celentano, Andrea Montanari, and Yuchen Wu. The estimation error of general first order methods. In Conference on Learning Theory, pages 1078-1141. PMLR, 2020.
[8] Michael Celentano, Chen Cheng, and Andrea Montanari. The high-dimensional asymptotics of first order methods with random data. arXiv preprint arXiv:21 2021.
[9] Michael Celentano, Andrea Montanari, and Yuting Wei. The Lasso with general Gaussian designs with applications to hypothesis testing. Ann, Statist. 51(5):2194 2023. ISSN 0090-5364.
[10] Kabir Aladin Chandrasekher, Ashwin Pananjady, and Christos Thrampoulidis. Sharp global convergence guarantees for iterative nonconvex optimization: A gaussian process perspective. arXiv preprint 2rXiv:2109.0985) 2021.
[11] Kenneth R Davidson and Stanislaw Szarek. Local operator theory, random matrices and banach spaces. Handbook of the geometry of Banach spaces, 1317-3666)1311 2001.
[12] Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. J. Amer. Statist. Assoc., 1348:1334811366 2001. ISSN 0162-1459.
[13] Cedric Gerbelot, Emanuele Troiani. Francesca Mignacco, Florent Krzakala, and Lenka Zde- borova. Rigorous dynamical mean-field theory for stochastic gradient descent methods. SIAM J. Math. Data Sci., 6(2):400 427, 2024.
[14] Peter Huber. Robust statistics, volume 523. John Wiley & Sons, 2004.
[15] Gen Li and Yuting Wei. non-asymptotic distributional theory of approximate message passing for sparse and robust regression. arXiv preprint arXiv:2401.0392. 2024.

---

[16] Mengqi Lou, Kabir Aladin Verchand, and Ashwin Pananjady. Hyperparameter tuning via trajectory predictions: Stochastic prox-lince methods in matrix sensing. arXiv preprint arXiv:2402.01599 2024.
[17] Yuetian Luo, Zhimei Ren, and Rina Barber. Iterative approximate cross-validation. In Interna- tional Conference on Machine Learning, pages 23083-23102, PMLR, 2023.
[18] Léo Miolane and Andrea Montanari. The distribution the Lasso: uniform control over sparse balls and adaptive parameter tuning. Ann. Statist., -9334)2313-22 2335, 2021. ISSN 0090-5364.
[19] Courtney Paquette, Kiwon Lee, Fabian Pedregosa, and Elliot Paquette. Sgd in the large: Average-case analysis, asymptotics, and stepsize criticality. In Conference on Learning Theory, pages 3548-3626. PMLR, 2021.
[20] Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and trends@) Optimization, 1(3):127-239, 2014.
[21] Pratik Patil, Yuchen Wu, and Ryan Tibshirani. Failures and successes of cross vols-vvaliotion for carly-stopped gradient descent. In International Conference on Artificial Intelligence and Statistics, pages 2260-2268. PMLR, 2024.
[22] Iosif Pinelis. Large deviations: Growth of empirical average of iid non-negative random varialbes with infinite expectations?, 2021. URL htt.aapp:/mmmtorreff.oo..ee/e/// Accessed: 2024-05-21.
[23] Kamiar Rahnama Rad and Arian Maleki. A scalable estimate of the out-of-sample prediction error via approximate leave-one-out cross-validation. J. R. Stat. Soc. Ser. B. Stat. Methodol,, 82 (4):965-999) 996, 2020. ISSN 1369-7412.
[24] Cynthia Rush. An asymptotic rate for the lasso loss. In International Conference on Artificial Intelligence and Statistics, pages 3664-3673, PMLR, 2020.
[25] Cynthia Rush and Ramji Venkataramanan. Finite sample analysis of approximate message passing algorithms. IEEE Trans. Inform. Theory, (4(11):7264472) 2018.
[26] Kai Tan, Gabriel Romon, and Pierre C. Bellec. Noise covariance estimation in multi-task high- dimensional linear models. Bernoulli, 30(3): 1695 1722, 2024. doi: 2.3150/23-BFJ166
[27] Robert Tibshirani. Regression shrinkage and selection via the lasso. R. Stat. Soc. Ser. Stat. Methodol., 58(1):667-288 288, 1996.
[28] Grace Wahba. A comparison of gev and gml for choosing the smoothing parameter in the generalized spline smoothing problem. Ann. Statist., pages 1378-1402, 1985.
[29] Ming Yuan and Lin. Model selection and estimation in regression with grouped variables. J. R. Stat. Soc. Ser. B Stat. Methodol,, (8(1):49-66) 67, 2006. ISSN 1369-7412.
[30] Cun-Hui Zhang. Nearly unbiased variable selection under minimax concave penalty. Ann. Statist., 38(2):894-942 2010. ISSN 0090-5364.
12

---

Supplementary Material of "Estimating Generalization Performance Along the Trajectory of Proximal SGD in Robust Regression"
A Additional simulation results
The following figures illustrate the proposed risk estimates accurately estimate the trajectory of the risk even when the step size changes at every step.
| b"ソ Ft |E||//n
|| 111
1E|1n - 1E|
8
Ft - |E|12/n
Value
Value
10 11 12 13 Iteration number
14 15 18 17 19
10 11 19 14 15 18 17
Iteration number
(b) if odd, and 1 0.5 if & is Figure 3: Risk curves for SGD applied to Huber regression with (nz, p) (3000, 1000) using different choices of step sizes. Left panel: 7}: if is odd, and 7). 0 if is even. Right panel: T/2 1 if is odd, and = 0.5 if is even.
(a) if L is odd. and fj: 0 if is
Figure 4 compares the performance of the proposed estimators Åt. F and the estimator jaub general- ized directly from (5) It confirms that proposed estimators outperforms jqub.
1 b )
|| 16 b²)
i |Ell2/n
ft |E|l/n
Ge - |
i - |E|2/n
Anwe 1/
pque ||E|/2/n
Value
Value
Iteration number
Iteration number
(a) Huber regression
(b) Pseudo-Huber regression
Figure 4: Risk curves for SGD applied to Huber and pseudo- Huber regression with (n, p. T) (4000, 1000, 20), |4| n/10 and 1/2 = 0.2 for all f.
13

---

B Auxiliary Results
Throughout, we define
E..... =
=
E] R
= (i x Sr xb E RaxI (13)
H 2 b", 6" E R
R= Ry ..... Xb' R (14)
Note that in the above E, H are not observable since E and b" are unknown. However, F and are observable and can be easily computed once the iterates (8' (E() are calculated.
B.1 Change of variables
In this section, conduct the change of variable to simplify the proof. Specifically, view the linear model = Xb" + E as a model with design matrix and the regression vector 0", i.e.
y = =b + XE
1/2 1/2 /// +e.
9"
This way, the design matrix G' has i.i.d. entries from standard normal distribution. Using the same argument in [5] Appendix DJ, 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 0* In other words, matrix X has i.dd NN, 1), or equivalently that the independent rows of X are normally distributed with covariance Σ Tp- We prove the main results using Σ I,. and the results for general E follow by this change of variable with the constant C(T, y, 7max co, 5) appearing in the bounds depending additionally
assume without of loss of generality that the design
(the upper bound of the condition number of Σ from Assumption 3.1)
B.2 Derivative formulae
In this section, we present derivative formulae that will be useful later proofs. The following formulas differ from [5] due the 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
Lemma B.1 (Proved in Appendix D.Ip. Let (b)) Jeery be the iterates generated from the recursion and the initial value b² is independent of X. Then the derivative bb' with respect X is given by
ab
(e, @ I,) ((FT,))
ee) 0e,) - (Ir 2 SD((H e,) & s.))
(15)
where r M L(A & D,) L Et= (( 4p), A Σ
[t. e and
I
P
1p
M
where P. = D.
nh X TTS S.D,X). [h.
Pr-1 Ip
Recall F = [S (3 (0-x6))1- S X67) I. we have Fu - X6'). -.b)) The following two corollaries are a direct consequence of Lemma B.I
Lemma B.2 (Proved Appendix D.2)). Under the same conditions of Lemma Let Fu
e, Fe, = e, Sab(y Xb'), we have
OF Dz.j
D A
(16)
where
D - e, S,D.e,e He, + ((e] H 8 e DS & XJI X (T)see @ @e),
A - ((e) F) & (Ir & X SDSDee & et).
14

---

Lemma B.3 (Proved in Appendix D.3). Let F = (y (oyy-bb'),.... (oyy-byyn)) and Ft. e Fer
Under the same conditions of Lemma B.1 We have
OF.
D + A
(17)
where
bj = e Die,e He, + ((e) H) 8 DDS(L @ (x) & x @e))
А ((e, F) 8e (Ir x Dlee @ ei).
Definition B.4. Define the matrices Yí RoxT, I2 E RexT I3 E RKT, Y1 E RTXT.
Ys RTXT XT
by the identities
de
Vj E [p].
ej HK -e,
(18)
Ozg
de, H
Vi E [n],
e, FW - e, e,
(19)
Orij
OF e H -KH FW -
(20)
Ori
HH X e H
(nIr
TT,
(21)
Ozij
OF e,e, X = - F
FÂT Ys.
(22)
dzy
where the matrices K, A, W are defined as follows
n
A= Σ
I (e)e)) & x)r(I @ SD & e,).
(23)
A = Σ (1 D(L @ xjr(Ir x (17
(24)
P
W Σ (Ir e [[r 0e,),
(25)
K Tr(S.D.)eee) T - Eur @ SD @ X)r(I, & x XJD(D 8 e,).
(26)
The matrices Y1: Y2... 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 YI...., Y, in Appendix D.4
Lemma B.5 (Proved in Appendix D.6). Under the same conditions of Theorem 3.6 with Ip. we have
E[[ Y & 1lop
E] S C(T, Mmax co) + 1|b*113)
max ke (1.3,4)
E 1
[)] I -CC(T 7, maxx co) + | 113),
E 2
E] < n² n'" 7, Tmxxx (62 11b*112)
15

---

We further define few matrices of size T x T:
θ F XH + KH H F' FW
2 = n [F" xx KH TH'X' X PF F F +F F
0 H X XH (nlT -
(27)
4 F F -KH + X +)
9 nH H (WF' H WWF H'
(WF
, E (E' XH ET-EEFFY FW
The next lemma provides the moment bounds for the Frobenius of these matrices.
Lemma B.6 (Proved in Appendix D.7 Under the same conditions of Theorem 3.6 have
max kE(1,23)
O
E] nC(I, Mmax. co) + |1b*112,,
(28)
max [8 (E(4,5)
E| I n C/T 7, Mmax co) ||6 11²),
(29)
E 1 E) - C(T, y: Mmax, Co (62 |bb 1(²),
(30)
almost surely, where [E[- E) is the conditional expectation given
We are able to prove the main theorems using Lemma B.6
C Proof of main results
C.1 Proof of Theorem 3.6
It suffices to prove this theorem for the case Σ = Ip. When E Ip. the result can be derived using a change of variables argument, as outlined in Appendix B.I By basic algebra, we have
O + P + +
=nH H (XH FW F (TH E FW + (XH-FWT) =nH H E E (E XH FWT) XH FWT)
=n. H E E -R+ FW T (R+ FWW.
is the t-th diagonal entry of (H H + E E/n). and is the 1-th diagonal entry of (R+ + FW T)T TTR + t/n.. Since | = /T||=l| using the previous display that conditionally e, we have
Notice that r, | --
E l
$) < n 'R O6s F |e]
2 le]
Using the moment bounds of θ5 and θ₂ in Lemma B.6 we have
B[rr B[mh-rr) - rol
E
C(T, n. Muaxx Co. 5)
(1 |E||
Vn
Furthermore, if E[e, is finite, have by [22] that |Ell/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 0 in probability. By enlarging the constant if necessary, assume C(T, 7/ Mmax, Co. 6) 2 1. To obtain a quantitative bound, by the conditional version of Markov's inequality, for any > O. and almost surely with respect that
P P(A. rel C|E) min (
C(T, Y, Tmaxx Co, 5)
(1 |||||
Vn
max
C(T, 7, Mmax. co. 8)
minf
1
|||||
vn
16

---

Taking expectation with respect to we obtain
P P(Ifs rel > <) max{1 C(T, y, Mmax: co. 8)
} min 1,
1 ||el| }]
vn
max(1, C(T,7 Thnax, Co, 8)
1
{ I
VT + + (1, ||=|] }D
1
with E[min ||s| H + (equivalently, || Pp 0) if the entries of are i.i.d. with a fixed distribution independent of n with [[[|,]]] < too by [22]. This finishes the proof of Theorem 3.6
C.2 Operator norm bound on K
We first recall the definition of K from (10) in the main text:
K Tr(
Σ (Ir o 2e
D(Lr @ X)I @ X JSD(LI & e,). (31)
Define two events: Ω1 = {X E REXP; | /op/ Vn 2 2 MH 2n }. where M is large enough constant such that P( > M) S 1/6.
vp/n) and i22 = (H E [n| |ail s
Lemma C.1. Under the same conditions of Theorem 3.6 with E = Ip we have in the event
Ω. 2 0 92 that
n| n||K llp S CC(T, y Muaxx Co 8, ||bb II)
Furthermore, Qt hus probability least 1 - e 1/18 n/2
Proof of Lemma C.7 Under Assumptions 3.1 and 3.5 we know that P(2)) 1 en/ from In the event Q1 we have by Lemma D.2 that
|X 1i/n S C(T, Y, 7max (o +1bb*113= = C..
Markov's inequality further implies
Hi [n] || H 2 3C, S n/33
In other words, Iii E [n] : ||= I < 3C, 2
2n 3
Recall that M is such that P(le;| > M) S 1/6.
By Hoeffding's inequality, we have
P
Me M} > P(e:: > M) + +a) S -2nan
Taking a = 1/6, we have
(P Σ
kf > M) 2 P( M) + 116) S nn/18 Furthermore,
using H E [n] |e. > M)l = >m 1[[=, > M). we have
lfe: > M] 2 1/6 1/66
c
1 lfee > M] 2 P((e M) 1)
i=1
Therefore,
P (Ift [n]
|.
My 2 n/3 -n/18
Equivalently, have P(K [n] |æ,| MP > 3m) > 1 - n/118 That is, at least 2n of the 3
entries of e are bounded by M with probability at least - n/18
Recall that 22 {l{i E [n] le,l S MM 2 ), then P(2)) > 1 c 1118 Hence, P(L, 0 0 2
1 n/18 n/2, In the event 22 0 Ωg the set
= fi E [n] e,ll < Mz H 2 3C.)
17

---

has size at least
3' For any i E I and E [T], we have
yi x. b' = e He & |e +
|0,
Hel S M + V3C..
(32)
By the definition of D1, under Assumption
3.4 have
Tr(Di) = Σ p" (y: x.(²)
Σ p"(y: x. 69)
since p" 2 0 by convexity
ief
2 |
min
p"(u)
due to 32
M+| 330.
2 n/3
min
p (u) = can
since has size at least ni/3.
Here e, a constant depending on p, M,, C. only.
By the definition of in (3). K K/n n is a lower triangular matrix with diagonal entries equal to Tr(D.)/n. It invertible if and only if all its diagonal entries are non-zero. Therefore, in the event 21 n92 we have R/n is invertible.
Let 7 = Et Tr. Then diagonal, | -lop maxier] Tr(D.] S (c.n)- and
(33)
where (S matrix. Using the upper bound [|Tlop in Lemma D.4 we have |Ellop S (T, Tmax Co) in the event 22, Now we rewrite K as
) Here L is a strictly lower triangular
K
KA A --
KA -
A
iλ 1)-1
Notice that ÊA-1 E T is strictly lower triangular matrix, thus
T-1
(T- - ビーー Σ (
1EA-11
k=
K=
By the triangle inequality,
T-1 Σ
[I - EX-1)-111 1lop S
LA
-1 11 S C(T, Mmax Co. 8, [[b")]).
k=0
Therefore, in the event 22 n S2g which has probability P(2, 0 2 1- -m/18 en/2
K - I llop (I ΣA - S - Y, Miuxx Cu, 5, 1|bb II).
Lemma C.2. Under the same conditions of Theorem 3.7 with Ip. we have
K op < n11 + |X 1 [(lo))),
|| I 2 [ W lop I n T
Proof of Lemma C.2
By the definition of K in using D < 1 and ||S||op S 1, we have
Jop S || + || loop S n(1
X op [|)p))
Similarly, by the definition of À in 24 we have
| llop n
X 2 r
Last, by the definition of W in 25 we have W llop S P|oo
18

---

Lemma C.3. Under the same conditions of Theorem 3.6 with Σ = Ip, we have
F|| F(A
KW)
Ir T
E] $ n (TT 7/ Mmxx, co) + 1|b^112); F e) S ( + v CT,, 7, Mmax co) -11b611))
B||R F(À
Kw)
Proof of Lemma C.3 First, using the definitions of 1, 2, e in we have
n1 F F(Ā - KW) = 0 K O O4.
(34)
Hence,
E|
F F(À KW) || |
e K
+ n + +n044) | |
by 34
E
ll
< K 1op E| nE 22/lp ||0+|| e|]
by the Cauchy-Sehwarz inequality
< <n 213/2C(T Tmax, (6 + ([b')
Here the last line uses the upper bounds of O e] from Lemma B.6 and the bound of
||||K|lop E] follows from Lemma C.2 and the bound of [T lop from Lemma D.4 This proves the
first inequality.
For the second inequality, we define
θ R XH
-
E
θ R XX KH X F FÄ
n El /|p//++1) θ R F (KH R (K + +FXX)T
n BIB\/++1)
where K = nIr
Σ n
-== [Ir @ (e) ( & X JSD(LI @ (X (TT)))) Using similar argument
that proves Lemma B.6
we can obtain the following bound of θ, 22 ④..
max E x $ S nC(T, 7, Mumax (6 | 1²). (E
(35)
E O1Il
$] S n 1/ /222 maxx co) + |bb
(36)
By the definitions of ④, θ 6, we have
[|E]||/ 1) K + θ O,1 = R F(À KW)T.
Therefore, conditional E, have
El|l R| KW) || [ = (| V O) K nnO, + n [
<n E V ()C Mmaxx ()) + 1bbb)
This finishes the proof of Lemma C.3
C.3 Proof of Theorem 3.7
In the event Ω4 = Ω2 0 12, we know from Lemma C.I that K is invertible and ||-- llop S 1C. Define W A. Using R + FWT = R + FW F(W WTT, have
(R + FW T (R FW
(R + FW T (R+ FW
= =(R+ + FW T F(W - W) + (W - W)F (R+ FW + (W W)F F(W- -W)T
We have by the triangle inequality
[| + FWT T (R + FWT (R+ FWT) ") (R+ FW )
(R T F(W W) + IIW W)F F(W -Willr
< R F(W - W) || + |WF F(W
W + IW - F F(W -)
19

---

Recall that in the event Q., we have
K 1oopp S. Using A KW = K(W W) and
Lemma we have
I(2 FW ) R FW
(R FW (R+ ( FW
2 ((2.) F(W W) |F
+ E ((0.) F(W W) 1|E |
+ ((( W)FFF(W' W)
1|
E(( F(A KW) (K) ||F
+ + ((O))) F(Ā KW) (K
)
+ +
+ T(O.)) KW)F F(A Kw) T (K
According to Lemma C.3 and the bound of K llop in Lemma C.1 the first conditional expectation is bounded from above by
n Hell + 1 )C(T 7max co (0 + |16*1
Using the bound of 1lop in Lemma C.1 the bound of W lop in Lemma C.2
and the bound of
E|[ F(Ā- KW) E] in Lemma C.3 the second conditional expectation is bounded from above by
n 1/2C( 7, Tmax (6 Co) + [1b*
Similarly, the third conditional expectation is bounded from above by
n 1/22 7, Tmax (8 + 11b*
In summary, we have
n ( + FW ) FR+FW - (R + FW F (R+ FW ")l
5 |لعاا vn + I)C(T, Tmaxx Co) (6² +11 12).
Since F, and fr are the t-th diagonal entries of (R+ (R+ (R+FWT) (RRTT))(RRR)///) and (R+ FR+FWT)T (R- FWT)/m )/n, respectively, have
T
B(((22)) 1ic - feller
< E ((2.)) R FW W
(R FW (R+ ||F
< + C) T, Mmaxx (5² 12)
Using the same argument in the proof of Theorem 3.6 have for any > 0,
P(T(2.)) - Ael > cl
& min (1, C(T,
V + ||=||]
S mmax(1, C(Trommo eo min(1, vn + Lell
Taking expectation with respect to E, we have
((((2_)) - Fxl > e) < max(1, C(T-N/maa [[mmi[[1, + [ell)
Using the union bound and PP(22)) 2 1 - e- nn/18 - n/2 > 2e n/18 we obtain
P P - Prl <) s 2e n/18 + max(1, ((I.n..Www....d.
vn min
Using the triangle inequality and the tail probability of I rıl in Theorem 3.6 we have
P lir - rel > e) S 2e n/18 + max(1, C(T. 1/ + [miin[1, |e|| []]
20

---

D Proof of results in Appendix B.2
D.1 Proof of Lemma B.1
By assumption, we know b' is independent of X and b Φ: ( +
X
Sy xb')) from
I
(4) Recall that Dr
aup(())
w= Xb and D,
ad,(w)
B xb Let a dot denote
du
de
b=
[r.
the derivative with respect to zij- we have
By product rule and chain rule and using y X x (b)-
T b'+1 D b' +
X SAb(y
XD) - x S, x69)
4
D; b' +
N
T Fe X S + xb
[1
where the last line uses Fi Subbyy XbO and H, = - b² Arranging terms gives
DD p [1
X S. S + b+t -i: D: F. X S.),
Let P: = D the above recursion of b' as linear system:
nt. X S,D,X) and a; D.(X F; X S,D,HH), we can rewrite [1:
I
0
a1
P1
Ip
02
:
Pr-1 Ip
I
0
ar
Solving the above system, have b' (e] M La. Since X further simplified as
0x
e,e]. a can be
Jayy
dt Di(e,e F; X S,Die,e H.). [1
a
Using D Dim((ee ((eee) & D.). D ((ee()))))) D.), ((ere))) St), and
A = Eí [T| e,e, we have
a = (A 2 D e F) - (Ir 8 X SD vec H)]
(A @ I, D [ e,) 2 [e,) - (Ir @ x JJDD((H e,) @ e.)]
Plugging this expression for a into b' = (c] @ IpM La gives
& 0b (e) M,M L(A @ I,D e,) ee, - (IT 2 X x SD((( e;) & e:)]
This finishes the proof Lemma B.I
D.2 Proof of Lemma B.2
By definition, Fi: = e, Sib(y XoO). By the chain rule of differentiation, we have
0Fu = - 0S XiO) -e S,Di(ese; He; + X Orgj
O
dzi
dr.z'
Notice that (e, @ (D.S)) DS(e, e @ I.). The desired formula then follows by plugging in
expression of ab in Lemma B.1
Bag
D.3 Proof of Lemma B.3
The desired identity follows by
aF az,
Oxbly xかO)
-c Dilere] Her X Øb'
e
Ozg
Oz,g
and the expression of 26 0z,
in Lemma B.
21

---

D.4
Alternative expressions for the matrices defined in Definition B.4
This section derives alteraativ expressions for the matrices Ii.... Y, defined in Definition B.4
T, we have by Lemma B.2 We first study Y in 18 Using F = Ei= Fe,e,
de, F
T
FF dzi
et
th e
ΣΣA
(37)
Ozi
t=1
i=1
t=1
i1
Now we compute the two terms in the right-hand side of (37 of Di in Lemma B.2
For the first term, using the expression
Σ D
ej H
Tr),,Di)e:ee" T + De (x D.S,e)))
[=1
e H
ee +ej H 2 (Ir 8 e, DS x
(1r ( x DDS(Ir 8e,)
HK
by (26)
Next, we compute the second term in the right-hand side of (37). Using the expression of A in Lemma B.2
-AjeT= = E ((eeFF) (e @ (x D(S,e))
i,t
= = [(e,, I (I x DS 0e))
The identity (18) then follows by substituting the above two expressions into (37
To study Y2 in (99), we use similar procedure. Using the mixed property of Kronecker product and the fact that the transpose scalar remains the same, we have
de dzss
de T'
c. -e (c] (L)((( ()(((F))))))) 0.
dzi
Σej(et @ [,)((r @ x XJSDD(HHHHH e,) e,jer by
15
j.t.
Σ (I 8e (Lr
e;
[(e]H) 0 (er @ejee
3.t
e, F E
[((e) (( x)r (lr 0e,)
= =e FWT ee,T
where W E1 0 (e])r((rr 8e,) and Y3 E, H L.)SS & xir e,).
22

---

To study Y3 in
20
we use the product rule of differentiation and
18 19
7
OF e,e; H
OF
de H
c
+
c
dzis
dzi
dzy
1=1
(KH's + e H - Σ ee FW T - e
K FW (Y H 2
For Y, in (1 by the product rule,
OH X e,e, H
Ozij
0H X ee, H H E eej H H x
HH dzij
Σ
Σ
e,e;
Ozij
1.
ij
0H
X
H H r2)
by 19
Ozij
1.
We then compute the first term of the last line as follows
0H
XT e,e;
X
H
Ozy
0HT
Σ
e,e,
X e,e) H
Ozij
ob
T
Σ isd
e:
X e.ee H
dzij
Ob'
Σ (t1
TeeT x
e, H
Urij
Eee, X(e [, & X SDD(((H (A((H)))))))) ()) H r,
by (L5)
(I & (e. x 3e) H
AH H T
by 23
where
r Σ e Xee @ (FF 0e))e H E &e (.](T(((FFFF e,) H.
Combining the above pieces, we have
H X XT c.e, H
T (FW - -22)
dzij
= (nIr
A H
FW
H X Y 81)
T
This provides alternative expression for T, in
21
23

---

Last, we study Ys in
22
We have
OF e F
OF ei Cej F e,e, eje F F e.e,ee] x
OF dwis
Ozij
Ozij
OF dzis
(KK + +
Σ
F e,e Feeee X
remains compute the third in the last display, Σ Ei= e Fere, we have by Lemma B.3 that
Using the fact that F
OF
dej Fe,
Σ 1,t
T
EeiDD + e
El
1,t
Using
1 e D He, + ((c) H) @ e DS @ x)r (Ir x @
А ((e, F) 8 (e (Ir @ x Dle, & @ei).
find
OF EF F Fee, e,j xTex
del Fei et
Σ eee ij
x x
ei
dri
dzys
i.s.
Σ F eee x e
Σ FFeee
x TeD.fe.T
We now compute the two terms in the above display. For the first term, we have
F x e
=E F FTee x
(e( (Ir & X TDD(e @ eijle
F
(I (
= F FÄT
by 24
For the second term, we have
Σ eeee x T
tight
Σ FTee, e,e, x xeel eil D.e.e] He, + ((e
x)r (x" De
Light
F D,XHe, + Σ ee XH 2 F
x
DS
x D
T2
Thus, we have established
OF ee,, X
F (T (T F r
dxy
This provides alternative expression for Y; in
22

---

D.5 Preparation results for proving Lemmas B.5 and B.6
Lemma D.1 (Moment bounds for H, F. F). Under Assumptions 3.1 H. F be defined in (3 we have for any finite integer
3.3 and 3.5 with E = Ip Let
B|I|X / Vn llppp I C(1,k),
B F
E)
S C(T, 77 7,,00 Tmmx 18 + 11b01))2*,
B F
2k |e] S RI|| 1/ E) C(T, k.)2e.
24
Proof of Lemma D.1 For the first inequality, according to Assumption 3.1 and Σ = 1,. X has ii.d. standard normal entries. By [1] Theorem I1], there exists random variable N(, 1) such that | llop S v + vD + almost surely. Under Assumption 3.5 that p/n n 7, we have n hp S C(+, k)) for any finite integer k.
X|||
For the second inequality, since | [6' b²11 it suffices bound | ba 1|2 for each /
[T]. Define the sequence of scalars d; max (1 8). Since D' 1 + X Sub(yy
X
I)) where == 4t/. Note that (yr. x), S V by Assumption 3.3 we have
I 0|| = |b' (0)1
since Φ1-1 0
< ||6² 15'-1_1 + 7.
The
x S - X6--1)1 I|
since Φ is 1 I-Lipschitz
16 11-1 1 + de
X ab (yI. X 17.60-1)1
T.
- 16 1 + m
X l
by the triangle inequality
< 1 + vioi
X || 8
by [4] 2 coT
S a1-1 + Aumns
X / 7p
= (1 +
Amax Vco
X /ntl--p)at-1-
Since a, max 1'11, 01 and 0 < 0/- 1, we have
at S (1+
Tmaa Vco
X /1lpooolp)a-11111
Notice 01 8 since 6'_ Op. we obtain
at S (1
X X/
x//n|llop) "
Hence, using the inequality |
bo1l2 < 2||6 6 112 + 2|| 1|2
we have
T
(I|b* + 1|b)
1=1
[1
Tauax dCe
X //////llppppp 21 21 )2
|
S T(2
13) +
Yaua vi
X
1 2T
(38)
Taking conditional expectation on both sides given e, the desired moment bound for H follows from the moment bound for ||X //nllopp-
For the third inequality, since |ي(²) < from Assumption 3.3 we have [[b((u)] S n22 for any E RR. By the definitions of F and F, we have
Elv(y- 21 < Tnd?.
Since the above display holds for any e, it implies the desired conditional moment bounds for F, F.
25

---

Lemma D.2 (Frobenius norm bound for X Under Assumptions 3.7
3.3 and 3.5
we have
|X P S C(T, Tmax co) + |6
with probability at least 1 - exxx(-n/n2)
Proof of Lemma D.2 By 11 Theorem II. II.13], under Assumption 3.1
with E 1p. we have
P[] P(| /V///lloop 2 - v1) 1 exp(-n/22
Using
X
X 1p
H and the bound 38 we have
X p nC 7, /aax (6 + |bb1)²)
holds with probability at least - exp( 22/
Lemma D.3 (Operator norm bound for M). Under Assumptions 3.7
3.3 and 3.5 have
IIM-llppp S ((())(++)),,
where { = Duma
X 1 lop
cpn
Proof of Lemma D.3
By the definition of M in Lemma B.1 we have
I
P1
I
M =
where P. = Dr(Ip
Be T. X
T S,D.X).
Pr-1 Ip
Hence, write M = IST A, where A is the lower triangular matrix with off aff-diagonal blocks
Pi..., .... Pr-1. Using the matrix identity (I A)) we have
E k= A and noticing A 0 for k > T.
P
O
M-1
k=
Pr
Taking operator norm on both sides, we obtain
T-
|M 1 lop
Σ k=0
Pa l
(39)
Since D; @v((w) du -=y--bi X6 and v is 1-Lipschitz, we have | I 1. By the definition that
S. = Eier, e,e,, we know || S 1. Since [ 2 con and nt S 7max for any E [T], we have
P oollop 1
1/t [1
X 1lop S 1 7max
X lop
def 1
con
Plugging this incquality into 39 gives
IM opp
Lemma D.4. Under the same conditions as Lemma D.3 we have
[ Dilp I n Mhaxx (1) + +)
where &
Tmax X lop- cp72
26

---

Proof of Lemma D.4 By the definition of I in Lemma B.1 we have r = L( @ I,) Notice
that A - T. e we have ||A||op max [7]] n Vimas using |4g Z con and
co
Ni S 7muax- Since Φ is I-Lipschitz, we have || op I By definition of L in Lemma B.1 we have
oooo = 1. Using these upper bounds of L op. || lopp DPllop and the upper bound of || M 'lopp
Lemma D.3 we obtain
|P op S || looll | Iop A @ Ipllop PPlpp
n C(T, 7mmax, Co) (1 + E)².
Lemma D.5 (Moment bounds for derivative of H, F, F). Under Assumptions 3.1 3.3 and 3.5 and
Σ Ip. we have for any finite integer k,
HH
k F
E
S C(T,7, Tmax co) - 1|6"113)24
n
OF
2 k
n
P
S
E
s C(T,7, 7, 7max co) + +|16"11322.
Swij
F)
Uzig
1=
Proof of Lemma D.5
We first prove the first bound. By Lemmal B.1 we have
Ob'
(e) - (Ir @ x (D((H)))) (,) 0.)]. (e) @ ())
Ozij
Hence, using
0c, Hee
da.y
T' Oze
we
have
de; He,
(IT & x = (e, @ -
x e, &e))
Ozij
= (e,
@ (p)e (Ir @ x SSD(H @ (,) @ e.)].
Using the above equality, Enj @ 0 A(e @ e;) = | p for A E RRTxxp, and the triangle inequality, we have
P
OH
2
P
P
de He:
Dmig
F
Ör.j
+ [[ &
X
SD
FFI + n | 2
X
2
||
| H }
op
$ pl
FI + M T op
X
2
H
op
< C(T, 7/ Tmmxx (1
+
FILE/n + C(T, 7, Muaax (o) 2) 2 /n/l |
where the last inequality uses Lemma D.4 Taking the conditional expectation both sides given the desired moment bound follows from the moment bounds X, H. in Lemma
Now we prove the second bound. By definition, the f-th column of F is Fi = S Xb'), it can be written using the t-th column of F as F; = S.F, Since
| lop S 1, we have
OF
OF:
O dzy
1| 0F 1| Ozi,
1/2h
12 <
1
By Lemma B.3 we have
de Fe:
D
À where
Dars
D = ej Dieee He: + ((e] H) @ el DDS @ XXT (Ir & x @ei)
((e) H) e 4r DS( x)r (1r & X xT)peee 8 ei).
A = ((e, F) 8 r (IT 8 x T ((e, 8 8e).
27

---

Using have
Eisualle] @e, A(ee & ee]l² = |All for A E R" and the triangle inequality, we
OF
76
P
D + 2
|I 1 H (1 + | [l + FI op X llop
DI & X)T (Ir @ x JDlli + (I & I,)r" (Ir @ x )DI
where the last inequality uses PPlpp - and |S 1. Taking the conditional expectation on
both sides given E, the desired moment bound follows from the bound of I in Lemma D.4 and the moment bounds of X, H, F Lemma D.I
D.6 Proof of Lemma B.5
Bound of T,
By the expression for the matrix Y E RRXT obtained in Appendix D.4
r Σ (f.TF) 2 1,) (Ir & x () 86)
T
((e Fe,e,) & (Ir x (I ee))
(e @ 11, (I, & (S Fe,
Σ (ef @ 11r (Ir x x DD & (Fe:)).
By the triangle inequality and |Pl| |S < 1. Y X
llop
< T op
F Using the
bound of [Top in Lemma D.4, the moment bound
lop, | |F in Lemma
D.I gives
Y ll |] s C(T, 7s Mmax, (o) + 1b)))
Bound of Y2. By the expression for the matrix Y2 E Rex obtained in Appendix D.4 we have
p r.
Σ ((e H) L.DS( (Ir 8e,)
((e] He,,,) 8 & @ x)r (In 8e,)
T
Σ E (eT & L.)D((I & x)r (Ir 0e,)e He
Σ (e] & I,DDS(I @ x)r (I, (He,)).
By the triangle inequality and | v ||S]lop $ 1, T2llop S T r lop X op H. Similar the moment bound of Y [lop- we obtain
Y 1
E] S C(T, 7, Tmax )n 182 + 1|6")))).
28

---

Bound of Y3. By the expression for the matrix Y3 E Rx obtained in Appendix D.4 we have Y3 (T) H T2). It directly follows that
|3 S || H| |F F( M2llop-
Using the triangle inequality and the moment bounds of bounds of || illops
H [p. || Lemma D.1 and the moment
|Y lllop we just obtained, we have
Tallpp
C(T, 7. Tmax (8 11b"112).
Bound of T. By the expression for the matrix Y, E R TXT obtained in Appendix D.4 we have
T4 H X T T,
where E, e (F e,) @ H. We can rewrite Yf as
n
Y Σ
Ir 8 e x)r((F H))
(1 2 e x)r((ere) F I,)
LI (e F ee x))[(ee & I)
Σ Ir F x))[(e, H.
We have by the triangle inequality,
| Ilppop - 7 D
[lop
X
lop
F
||
H IF-
By the triangle inequality and the upper bound of |T
we have
|[Tallop - Y Ill X
lop
H
|F
+
| llop S T x lop H +
[lop] H [|F.
Squaring both sides and taking conditional expectation, we have
[
E) < C(T,"
Mmax, (52 + (bb1)
thanks to the upper bound of |Tpp in Lemma and the moment bounds of |X lop [| || || in Lemma D.I
Bound of Ys. By the expression for the matrix Y3 TxT obtained in Appendix D.4 we have
Y; = T. x F Y2, where
T
r2 = FF D,XHe, (ee
XH
F xyr (Ir @ X (D([ @e).
1=1
By the triangle inequality.
llpp S X
|lop
F
F|| H
+
||
ol Top
F||r
H
P.
Using the moment bounds of Y /rll
X
XXllpp opp |
H
we have
[ Y 11
E]
< n
C(, %, Tmaxx (82 + 11)...
D.7 Proof of Lemma B.6
We first state three useful lemmas.
29

---

Lemma D.6 (Adopted from Lemma E.10 of (6)). Let U, V R - R be two locally Lipschitz,function of Z with i.i.d. (0, 1) entries, then
E U T zV
(UT (eeeeV)
| RF.
dziz
dU
0v dzij
B + E
+ + 1F
Ozy
Lemma D.7 (Adopted from Lemma F.5 of (.). Let U, V TR - RaXT be two locally Lipschitz
functions of Z with i.dd N(O, 1) entries. Provided the following expectations are finite, we have
d.g Joe,. ej ezzzv) [|]
U e, U Ze,
S <11+2/D + + |/ + B|| 1/2 EII 11
where = 00 -, /0z; and | ||a Cc Σ || U 1/2
We will use the above two lemmas, conditionally on e, to bound the conditional moments of 1 22 2g 4 and O₃ given E.
Bound of O By the definition of 1 we have
OF eie] H
F
+ Ys
by (20)
=0+ Y3.
Applying Lemma
D.6 conditionally on E to (Z,U, V) = (x, F. H) gives
0F T H
1
|e]
E
XH
|6| Blll 1 E)
OF Fz.
HH dzis
E
11
e
+
[B[
Y l
e]
S (T 7, Tmax (82 + 11),
where the last line uses the moment bounds of F F|rr, aH dess
||r in Lemma D.I the moment bounds of OF [| in Lemma D.5 the moment bound of || in Lemma B.5
[r H
[F.
0wg
Bound of O2.
By the definition of O2 we have by 21
FF
e,e, T
FXX
F K F -AF-PPTI FF +Ts
Ozij
=nQ + Ts.
Applying Lemma D.6 conditionally on E to to, UZ,U, V) (X, F. F. F) gives
OFT e,e; X
|
E
X
E
+
El|
Y
||
[=]
Jwij
OF
ax F
X
+ E
X F
1
F
|| |
+
|e]
F.:z
Oziz
OF
SEI 18p +E+EE (1
X
2
2
+
TE
Trs 2 E)
< <n C(T, 7, MInxx co) + 11b
30

---

Here, the penultimate line uses ||
FIFE.
OF B2cs
1F
aF 0es
||,
and
1 F
||
llefe e; aF
oF Uzuy
or
|
lop
aF
Ir. The last line uses the moment bounds of
8mg
Bang
[X|lop
||
in Lemma
D.1 the moment bound of
OF Jany
||
Lemma D.5 the moment bound of
Mallop in Lemma
B.5
Bound of 33 By the definition of O we have
0H
X e H
H H'X
dz:j
4
= X - A)H
H H x FW Y4
by 21
= = 4.
Applying Lemma D.6
conditionally on e to (2, V) = (X, XH, H) gives
HH
X T e H 2
E
1 |
XH
E
a1 E]
am
0H
OX H
X
E] + E
11
+
|
E + T
alll
Dz:
izi
OH
op
H
[ + (1 +
X
H
|
+
+T Y
Tall
|E]
ij
SnC(T, 7ss Imax; Co) (52 + 11b
Here, the penultimate line uses
0XH ||| ||e,e] T 0H
+ X OH
|| S (1+ (1+||X|(p))
0H dugy
||r, and the
dag
last line uses the moment bounds of |X (lop, H | in Lemma D.I the moment bounds of 2H
||F
Bagg
in Lemma D.5 the moment bound of || Tallop in Lemma B.5
Bound of 4. By definition, we have
"
TF e; (KH +rT)e, and
n OF e;
Σ
(KH + +
Ozij
amij
1=1
Using Σ}=1 c = Ip. we find
F e; F xe,)
OF T ei FT Xe,)
dzg
Uı.
(KH
F X + r
KKH + X TT)
This further implies
R OF e:
OF e FT Xe,)
PF F
F xTxee)
dzy
dzg
1=1
= =nO +
(K
F X) + (KH + F X + r
T.
Applying Lemma D.7
conditionally on E to (Z,U, V) = (X, F,F.FF.
OF e
F Xe,) P
OF FT xTxe,) 1| | + +EIl| E|ll T [|]
|
dzsj
dzi
1=1
21++ (E F |2 + E|| F e2
E||
e2 +
22 + TE|
< C/, 7, Imax, (52 + ||b)
Here, the last inequality uses the moment bounds of |F FlFF in Lemma D.I and the moment bounds of 1F... | in Lemma D.5 and the moment bound of || olop in Lemma B.5
31

---

Bound of 5.
By the identity 19 we have
OH ej H
0H
X
e;)
H X ai)
dzij
j=1
=( H r (WF
F H x - - T.)T
By the definition of s, we have
nH H OH ej H T x e,))
H ej
H x ei)
dzy
dzy
=, + Y 2 - H (WF H 12 -2
f,
Here Y. =2 ((e H)& & DI, & X)T (I, (1 &
Applying Lemma D.7
conditionally on e to (Z,U, V) = (XT, HH (i.e., consider the mapping
from R2 to [ X H) gives
E
E)
SE
HH Ej dz.j
H H es)
H ej
H
e:)
|1 F
[e]
+ + le]
Fzij
(1+
|/ + E|| H l =)//2) + TE
[||
ir
5
llap
Sn Tmaxx (82 + 11b1)
Here, the last line use the moment bounds of H |F in Lemma D.I the moment bound of | 13 in
Lemma D.5 and the moment bound of | lop in Lemma B.5
Bound of Og.
Using have
p O e H
E. (FFT T2))
Dzgj
It follows that
12 dE' ete] H
XH
E XH ET (F T2)
Jzij
|B||: O5 + E Y2.
Thus, using E = E/| [r. we have
1
O e H ET T2. TT- P
=5=E XH
Jzij
i=
Applying Lemma D.6
conditionally on to (Z,U.V V) (X, E, H giies
O e H
E| 1
e)
E
E HH
e + E r211
Fzg
ÖH
E] + E
11
+ E
Y2ll
Dris
0H
E| +
+
B|Jl
E]
D..
< C(T, 7. 7max cn (0² + bb" 12).
Here, the last line uses the moment bound of [ || in Lemma D.I the moment bounds
8H
[F in
Basy
Lemma D.5 and the moment bound of Y in Lemma B.5 his finishes the proof of Lemma
B.6
32

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract and introduction clearly state the claims made in the paper.
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 perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
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 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 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 few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might 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 address problems of privacy and faimess.
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 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
33

---

Justification:
All the proof' are provided in the supplemental material.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: All the code needed reproduce the results is in the supplemental material.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
34

---

Answer:
[Yes]
Justification:
All the code are provided in the supplemental material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.cm
cesssnnnoiccyy for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ ccpppi.ccppulli//uuuii wwwiicccguie///eeeeeeeeis CoodeSubeiiie ssi wswonnolli onnolicy)))) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
The choices of step size for GD and SGD are provided in Section 4.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All the figures include the 2-standard error bars.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, at least for the experiments that support main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
35

---

It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of enrors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: is provided the README file of the supplemental material.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code Of Ethics
Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Question:
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
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 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.
36

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is 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 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
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, scraped datasets)?
Question:
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 such risks.
Released models that have 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 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 require this, but we encourage authors to take this into account and make 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 other code 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 dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.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 Datasets.papersstt paperswi thcode. cow/ datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of 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.
37

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS 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.
38


Reviews:
[{"review_id": "TF1fa7hOzn", "reviewer": "Reviewer_3EXM", "summary": "#### Summary\nThis 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.", "questions": "#### Questions\n1. **Computational Complexity**:\n    - Could you provide more details on the computational complexity of the proposed estimators? How do they scale with increasing dataset size and dimensionality?\n\n2. **Generality**:\n    - 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?\n\n3. **Comparison with Existing Methods**:\n    - How do the proposed estimators compare empirically with existing methods for estimating generalization error? Are there scenarios where your approach significantly outperforms others?\n\n4. **Practical Implementation**:\n    - 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": "no", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "#### Strengths\n1. **Innovative Methodology**: The introduction of estimators that track the generalization error of iterates is novel and provides a significant contribution to the field of robust regression.\n2. **Theoretical Rigor**: The paper provides a thorough theoretical foundation for the proposed estimators, including consistency proofs and detailed mathematical derivations.\n3. **Practical Relevance**: The approach is applicable to a variety of robust regression problems, including those with heavy-tailed errors and non-smooth regularizers.\n4. **Empirical Validation**: Extensive simulations demonstrate the effectiveness of the proposed estimators in tracking the generalization error and determining the optimal stopping iteration.", "weaknesses": "#### Weaknesses\n1. **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.\n2. **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.\n3. **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.\n4. **Practical Guidelines**: While the theoretical results are robust, practical guidelines for implementing the estimators in real-world scenarios are not sufficiently detailed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the insightful and encouraging feedback. Here, we respond to your comments point by point.\n\n> **Q1**: Could you provide more details on the computational complexity of the proposed estimators? How do they scale with increasing dataset size and dimensionality?\n\n**A1:**\nWe first note that our proposed risk estimates, $ \\hat{r}_t $ and $ \\tilde{r}_t $, only require computing the iterates $ \\hat{b}^t $ and the weight matrices $ W $ and $ \\hat{W} $. Since $ W $ is a $ T \\times T $ lower triangular matrix, and the computational complexity of each entry of $ W $ is $ O(npT) $ using the Hutchinson trace estimator, the total computational complexity of $ W $ is $ O(npT^3) $.\n\nSimilarly, for the computation of $ \\hat{W} = \\hat{K}^{-1} \\hat{A} $, the computation cost of both $ \\hat{A} \\in \\mathbb{R}^{T \\times T} $ and $ \\hat{K} \\in \\mathbb{R}^{T \\times T} $ is $ O(npT^3) $. Thus, the overall computational complexity of $ \\hat{W} $ is $ O(npT^6) $, with $T^3$ coming from inverting the triangular matrix $\\hat{K}$. (Note that in practice, for numerical stability we do not compute the inverse directly, but instead solve the corresponding linear system). \n\nOverall, the implementation (provided in the supplementary material) for separable loss/penalty to compute $ \\hat{b}^t $ and the weight matrix avoids any operation that would be $O(\\min(n,p)^3)$ such as multiplying two matrices of sizes larger $\\min(n,p)\\times \\min(n,p)$, for instance computing the full Gram matrix $X^TX$. The implementation only performs matrix-vector products with matrices of size smaller than $\\max(n,p) \\times \\max(n,p)$, or multiplication of a dense matrix by a diagonal matrix both of sizes smaller than $\\max(n,p) \\times \\max(n,p)$. It never incurs an operation with complexity larger than $\\max(n,p)^2$ (ignoring here the dependence on $T$). This makes it possible to run the implementation even on a laptop for $n,p$ both 10,000. \n\n\n> **Q2**: 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?\n\n**A2:**\nOur analysis focuses on the performance of the iterations that can be solved using SGD or proximal SGD methods.\nThis includes any gradient-Lipschitz loss function (for instance Huber) with a non-smooth penalty as illustrated in the simulations.\n\nIf the regression problem has both a non-smooth data-fitting loss and a non-smooth penalty, for instance\n\n$$\n\\hat{b} = \\arg\\min_{b \\in \\mathbb{R}^p} ||y - Xb||_1 + \\lambda ||b||_1,\n$$\n\nwe expect that one cannot use the analysis and algorithms studied here (proximal GD or proximal SGD) due to the non-differentiability of the data-fitting loss. For such optimization problems, other primal-dual algorithms are needed, such as the alternating direction method of multipliers (ADMM) and the Chambolle-Pock algorithm. However, these algorithms have a different structure, and we expect for these algorithms a significantly different result than our Lemma B.1 for instance. Because of the different structure, we expect these algorithms to require different analysis and risk estimates.\nWe leave this as an interesting direction for future work.\n\n> **Q3**:\nHow do the proposed estimators compare empirically with existing methods for estimating generalization error? Are there scenarios where your approach significantly outperforms others?\n\n**A3:**\nWe are not aware of any existing methods that can estimate the generalization error of the iterates of proximal SGD algorithms in our settings of the proportional regime ($ n\\asymp p $). In this regime, cross-validation with a finite number of folds is known to be inconsistent for instance [2, Figure 1], and diverging number of folds would be impractical computationally. Thus we did not provide a direct comparison (though we would be happy to provide empirical comparisons if reviewers suggest existing competing risk estimates that we missed). One related work is by Luo et al. (2023), which proposed an estimate for the cross-validation error of the iterates of SGD algorithms. However, their method has many restrictions on loss function and does not work for high-dimensional regression settings with $p > n$.\n\n\n> **Q4:**\nCan 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?\n\n**A4:**\nThank you for this question.\n\nWe provide the following practical guidelines for implementing the proposed estimators:\n\n1. If practitioners are solving a regression problem with a smooth loss function and are using SGD or proximal SGD algorithms, they can use the proposed estimator $\\hat{r}_t$ to estimate the generalization error of the iterates if the covariance matrix of the features is known (either because lots of additional unlabeled data are available, making estimating $\\Sigma$ possible, or because the practitioner samples the design $X$ themselves). Otherwise, if the covariance is unknown, use $\\tilde{r}_t$ which uses out-of-bag samples to maintain good performance (cf. answer to Reviewer 2 for a discussion on the use of out-of-bag samples for SGD). \n\n2. Once $\\tilde{r}_t$ is computed, plot the estimated generalization error of the iterates as a function of the number of iterations. Use this plot to choose the stopping time that achieves the smallest out-of-sample error, or as an additional tool to diagnose to study the convergence of the algorithm at the population level. \n\n---\n[1] Luo, Yuetian, Zhimei Ren, and Rina Barber. \"Iterative approximate cross-validation.\" International Conference on Machine Learning. PMLR, 2023.\n\n[2] Rad, Kamiar Rahnama, and Arian Maleki. \"A scalable estimate of the out-of-sample prediction error via approximate leave-one-out cross-validation.\" Journal of the Royal Statistical Society Series B: Statistical Methodology 82.4 (2020): 965-996.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FxfTqf28EA", "reviewer": "Reviewer_1ZGt", "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.", "questions": "There are no further questions.", "limitations": "This work is purely theoretical and has no negative influences.", "rating": 4, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "This paper is generally quite well-written. The proof and result appears to be reasonable. Empirical evidences confirm that the proposed estimator is indeed effective and accurate.", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the comments and the opportunity to clarify our contributions.\n\n> **Weakness:** \nDespite being more general, this paper seems to constitute limited progress towards the community. When comparing it with the work by Bellec and Tan [5], it is not hard to notice the estimator 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. (...)\n\n## **Response:**\n\nWhile some notation and techniques (in particular probabilistic lemmas based on Gaussian integration  by parts) are reused compared the previous result of [5] for the square loss (we did not try to hide this with new notation, and we did not reinvent the wheel unless necessary), let us point out some significant technical differences.\n\n- In terms of applicability, estimators applicable to SGD (or proximal SGD) and robust errors are a significant step forward compared to the full-batch and square loss setting of [5].\n- For the square loss, residuals and gradients are the same. For robust loss functions or the iterative algorithms studied here,\nresiduals and gradients are different; this defines a different structure which can be seen in the estimators (involving both the residuals\nand the gradients), as well as in the weight matrices $\\hat A, \\hat K, W$. In the square loss only two weight matrices appear, and it was a surprise to us to find the necessity to introduce three matrices to analysis the problem for proximal SGD or robust losses.\n- [5] studies Gaussian noise. We allow heavy-tailed noise with infinite variance. This requires different tools to control the noise,\nand the resulting rate is also different, with the rate explicitly depending on the noise through the quantity $\\mathbb{E}[\\min\\\\\\{1, ||\\varepsilon||/n\\\\\\}]$\nappearing in Theorem 3.6.\n\n**More importantly, directly generalizing the approach in [5] fails for SGD.**\nThis failure for SGD resides in the difference between the matrices $\\hat K$ and $\\tilde K$\nin Equations (26) and (27). Generalizing the approach of [5] leads our research to the\nmatrix $\\tilde K$ in (26). In order to build an estimate of the weight matrix $W$ of interest\n(as discussed in the submission, $W$ cannot be used directly as $\\Sigma$ is typically unknown),\none wishes to invert $\\tilde K$ in the approximation $\\tilde A \\approx \\tilde K W$.\nThis inversion fails for SGD for small (but still very realistic) batch sizes of order $0.1n$.\n\nThe matrix $\\tilde K$ is lower triangular, and the reason for the lack of invertibility \nof $\\tilde K$ can be seen in the diagonal terms equal to $\\text{Tr}[S_tD_t]$ in (26),\nwhere $S_t \\in \\\\\\{0,1\\\\\\}^{n \\times n}$ is the diagonal matrix with 1 in position $ii$\nif and only if the $i$-th observation is used in the $t$-th batch.\nThis diagonal element of $\\tilde K$ can easily be small (or even 0) for small batches,\nif the batch only contains observations such that $(D_t)_{ii}$ is 0 or small.\n\nFor SGD and proximal SGD, we solved this failure of the invertibility of $\\tilde K$ by using out-of-batch\nsamples in the construction of $\\hat K$ and $\\hat A$, in order to avoid $S_t$ in the diagonal elements\nof $\\hat K$ in equation (27). This is the key to making these estimators work for SGD and proximal SGD,\nand this use of out-of-batch samples could be anticipated by reading or\ngeneralizing [5] (which only tackles the square loss with full-batch\ngradients).\n\nThis phenomenon is seen in generic SGD simulations, for instance with the Huber loss and\n$n, p, T = 4000, 1000, 20$ and batch_size equal to $n/10$. From iteration 10 to 19, the diagonal elements of $\\tilde K$ are close to 0, while using out-of-batch samples in $\\hat K$ provides diagonal values bounded away from 0, and thus numerically stable invertibility of the triangular matrix $\\hat K$:\n\n| $t$  | $\\tilde{K}_{tt}/n$ |   $\\hat{K}_{tt}/n$ |\n|---:|----:|----:|\n| 10 | 0.04 | 0.4  |\n| 11 | 0.05 | 0.42 |\n| 12 | 0.04 |  0.43 |\n| 13 | 0.04 | 0.44 |\n| 14 | 0.04 | 0.46 |\n| 15 | 0.05 | 0.47 |\n| 16 | 0.05 | 0.49 |\n| 17 | 0.05 | 0.5  |\n| 18 | 0.05 | 0.51 |\n| 19 | 0.06 | 0.52 |\n\nWhile using $\\tilde K,\\tilde A$ (without using out-of-batch samples) is successful for estimating the risk for large batch sizes (above $0.3n$), it quickly deteriorates as the batch size decreases: in the same setting as above, with 100 repetitions, the true and risk estimates using $\\tilde K$ and $\\hat K$ give\n\n| $t$ | True risk | Estimate using $W$ | Using $\\hat A \\hat K^{-1}$ | Using $\\tilde A \\tilde K^{-1}$  | \n|:---|---:|----:|----:|----:|\n| 10 | 4.67 | 4.66 | 4.66 | 4.4  | \n| 11 | 4.29 | 4.29 | 4.29 | 3.98 | \n| 12 | 3.94 | 3.94 | 3.94 | 3.64 | \n| 13 | 3.63 | 3.63 | 3.63 | 3.33 | \n| 14 | 3.35 | 3.35 | 3.35 | 3.04 | \n| 15 | 3.09 | 3.1  | 3.1  | 2.77 | \n| 16 | 2.86 | 2.87 | 2.87 | 2.54 | \n| 17 | 2.65 | 2.66 | 2.66 | 2.34 | \n| 18 | 2.46 | 2.47 | 2.47 | 2.17 | \n| 19 | 2.29 | 2.3  | 2.3  | 2    | \n| 20 | 2.13 | 2.14 | 2.14 | 1.83 | \n\nThe direct generalization from [5], that does not leverage out-of-batch samples (right column), is inconsistent while our proposed estimate leveraging out-of-batch samples is consistent.\n\nNote that beyond simulations, it wasn't clear at first that using out-of-batch samples would provably work. After significant trial and error we eventually found the serendipitous combination (cf. line 440) of the probabilistic identities that grants the approximation $\\hat A\\approx\\hat K W$ in the row space of $F$ in eq (35).\n\n**Thanks for raising this point.** Highlighting and explaining this significant departure from [5] on the use of out-of-batch samples is something that was admittedly overlooked in the main text of the initial submission, and we will use the extra page available for the camera ready version to clearly explain this.  \n **We kindly suggest to reconsider the review rating in light of this.** We would be happy to provide further clarifications if needed.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "fe5uO3GZnX", "reviewer": "Reviewer_yB9W", "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.", "questions": "- What are the dependencies on $T,\\eta_{max}$ in the numerators of Theorems 2.6 and 2.7?\n- 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": "The authors have adequately addressed the limitations of their work.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is overall very well written, including the overall motivation, the related work, and the technical ideas behind Theorems 3.6 and 3.7.", "weaknesses": "- 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.\n\n- 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.\n\nMinor Points:\n- line 29: starts -> start\n- 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)\n- line 170: all -> both\n- 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}}$?\n- line 269: the three ... reveal", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Q1:** \n(...) dependencies on $T, \\eta_{max}$ in numerators of Theorems 2.6,2.7?\n\n**A1:**\nThe dependence on $T$ is $T^T$ and the dependence on $\\eta_{\\text{max}}$ is $\\eta_{\\text{max}}^{T}$, as can be seen in Lemma D.4. We do not expect this bound to be tight. Simulation results confirm that the proposed risk estimate is still accurate for all iterations, even when $T > n$, at least in the simulation settings that we tried.\n\n> **Q2:** \nIs the reason that $\\hat W$ must be computed recursively simply to invert $\\mathcal{M}$ ? (...) how to connect the unrolling in section 3.1 to the Kronecker expressions in 3.2.\n\n**A2:**\nYes, the recursive computation of $\\hat W$ is related to the inversion of the triangular matrix $\\mathcal{M}$.\nTo expand on this, as intuited in Section 3.1, unrolling the derivatives by the chain rule brings a matrix product of previous\nJacobian matrices from the current iteration to the first. The connection to $\\mathcal{M}^{-1}$ is easier to see\nfrom the identity\n\n$\n\\begin{pmatrix}\n    1 & 0 & 0 & \\cdots & 0 & 0 \\\\\\\\\n    -p_1 & 1 & 0 & \\cdots & 0 & 0 \\\\\\\\\n    0 & -p_2 & 1 & \\cdots & 0 & 0 \\\\\\\\\n    \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\\\\\n    0 & 0 & 0 & \\cdots & 1 & 0 \\\\\\\\\n    0 & 0 & 0 & \\cdots & -p_{T-1} & 1\n\\end{pmatrix}^{-1}\n$\n$=\n\\begin{pmatrix}\n    1 & 0 & 0 & \\cdots & 0 & 0 \\\\\\\\\n    p_1 & 1 & 0 & \\cdots & 0 & 0 \\\\\\\\\n    p_1 p_2 & p_2 & 1 & \\cdots & 0 & 0 \\\\\\\\\n    p_1 p_2 p_3 & p_2 p_3 & p_3 & \\ddots & \\vdots & 0 \\\\\\\\\n    \\vdots & \\vdots & \\vdots & \\ddots & 1 & 0 \\\\\\\\\n    p_1 p_2 \\cdots p_{T-1} & p_2 p_3 \\cdots p_{T-1} & p_3 p_4 \\cdots p_{T-1} & \\cdots & p_{T-1} & 1\n\\end{pmatrix}.\n$\n\nOn the left, the matrix has ones on the diagonal and $-p_1,...,-p_{t-1}$ just below the diagonal, like our $\\mathcal M$\nbut with blocks of size 1. On the right, the matrix is similar to the product of Jacobian matrices brought by unrolling the derivatives\nby the chain rule. With this in mind, the matrix $W$ in equation (8) is obtained by\ncontinuing to unroll the derivatives by the chain rule for $t=3,4,...$ as for $t=2,3$ in Section 3.1 (see also equation (19)).\n\nUsing $\\mathcal M^{-1}$ is useful to manipulate a compact notation with no product over $T$ matrices.\nHowever, for the implementation the matrices are computed recursively by leveraging the above structure\n(for instance, if we have computed the $t-1$-th row of the inverse, most of the next row is easily obtained by multiplying the $t-1$-th row by $p_{T-1}$).\n\n> **Weakness 1:** \nNumerators in Theorems 3.6-3.7 include unspecified $C(T,...)$. From the appendix, it appears that $C(T,...)$ scales like $T^T$, so that you can only run for \n$T \\approx \\log(n)/\\log\\log(n)$ steps before the bound is vacuous (...). 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.\n\n**Response:** \nYes, our current analysis provides a bound on the iteration number as $T^T$. We expect this dependence in $T$ to be an artefact of the proof. As illustrated in Figures 1-2 and further supported by additional simulations with $T > n$, the simulations suggest that the generalization error estimate is accurate for all iterations, even as $T$ diverges. \n\nNote that even though the theory of the submission incurs these constants of order $T^T$, this still allows us to pinpoint \nwhat the risk estimate should be for all $T$.\nWhile the theory only holds for $T \\le \\log(n)/\\log\\log(n)$ for which the risk estimate is consistent,\nonce the definition of the estimator is known,\nrevealing the form of the estimator is still useful:\nonce we have identified the definition of the estimate:\n\n- simulation studies such as those provided in the submission allow\nus to check the empirical validity of the estimate beyond $T \\le \\log(n)/\\log\\log(n)$, and \n- it will be easier for future work to improve this dependence in $T$, now that identifying what the risk estimate should be has been done.\nIt is common for new estimators to be proposed first, before further works refine the bounds on these estimators.\n\nImproving the dependence on $T$ appears difficult and possibly out of reach of current tools,\neven in the well-studied Approximate Message Passing (AMP) algorithms (which is also an algorithm typically studied\nin the proportional regime of interest here). The papers [1, 2] feature for instance the same dependence \n$T \\le \\log(n)/\\log\\log(n)$ for approximating the risk of AMP.\nThe 2024 preprint [3] offers the latest advances on the dependence on $T$ in the bounds satisfied by AMP.\nIt allows $T\\asymp \\text{poly}(n)$ to control certain AMP related quantities, although for the risk [3, equations (16)-(17)]\nthe dependence $t=O(\\log n)$ is required which is still logarithmic. This suggests that advances on this front are possible, at least\nfor isotropic design for separable loss and regularizer such as those studied in [3]: Lasso or Robust M-estimation with no regularizer.\n\nSince these latest advances in [3] are obtained for specific estimates (Lasso or Robust M-estimation with no regularizer),\nit may be possible to follow a similar strategy and improve our bounds for specific examples of iterative algorithms closer to AMP \nor featuring only separable losses and penalty. But we believe following such a strategy for specific examples would be \nout of scope for the current submission, which tackles a general framework allowing iterations of the form\n(4) with little restriction on the nonlinear functions\n$\\phi_t,\\psi$ except being Lipschitz (and bounded by $\\psi$).\n\n**Minor points**: We agree and will implement the suggestion regarding line 207 in the camera-ready version. Thanks!\n\n* * *\n\n[1]: \"Finite Sample Analysis of Approximate Message Passing Algorithms\".\nRush and Venkataramanan, 2016.\n\n[2]: \"An Asymptotic Rate for the LASSO Loss\". Rush, 2020.\n\n[3]: \"A non-asymptotic distributional theory of approximate message passing for sparse and robust regression\".\nLi and Wei, 2024.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "uLvQeDAB0k", "reviewer": "Reviewer_boVy", "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].", "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?\n\n2) 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?\n\n3) 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?\n\n4) 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?\n\n5) 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?\n\n****\nI have increased my score to 'Accept' post rebuttal.", "limitations": "The paper discusses limitation of the current approach as a way to motivate follow up research in the conclusions.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors have developed a procedure that reliably determines the best stopping time for solvers in contexts where no consistent estimators were previously available, particularly in robust regression with heavy-tailed noise, Huber-like loss functions, and non-smooth penalty terms. This methodology offers an alternative to cross-validation techniques, which are generally inconsistent (e.g., V-fold) or computationally expensive (e.g., leave-one-out). The paper is well-written and includes clear proofs.", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the comments, we provide our response below. \n\n> **Q1:** \nIn 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?\n\n**A1:**\nIf the squared loss and full-batch gradient descent is used, our proposed $\\tilde{r}_t$ has the same formula as $\\hat{r}_t$ in [5]. However, our work allows the noise vector can be heavy-tailed (e.g., infinite variance), which is significantly different from the Gaussian noise condition in [5].\nTo handle heavy-tailed noise, we consider a data-fitting loss whose derivative is Lipschitz continuous and bounded (Assumption 3.3). Hence, the squared loss is excluded by Assumption 3.3 of our current paper, so the result of [5] is not formally implied by the theorems of the submission.\n\n> **Q2:** \n...can the authors give examples of how their methodology can be use to tune other tuning parameters?\n\n**A2:**\nLet's take the L1 regularization parameter $\\lambda$ as an example. In proximal SGD, $\\lambda$ appears in the soft-thresholding parameter of Example 2.4. The algorithm's iterates, $\\hat{b}^t$, depend on the choice of $\\lambda$, so we can write it as $\\hat{b}^t(\\lambda)$. For each candidate $\\lambda \\in \\{\\lambda_1, \\ldots, \\lambda_K\\}$ (a finite grid of tuning parameters), we compute the estimated risk $\\tilde{r}_t(\\lambda)$. This estimate serves as a criterion to choose the $(t^*, k^*)$ that minimizes $\\tilde{r}_t(\\lambda)$. As long as the grid of $\\lambda$ is finite, the risk estimate for $\\hat{b}^t(\\lambda)$ is consistent across all parameters in the grid. The selected $(t^*, k^*)$ leads to the smallest generalization error over the grid, up to a vanishing error term.\n\n> **Q3:** \n... a specific comparison of their methodology and consistency results to other classical methods ... would be beneficial.\n\n**A3:**\nWe are not 100% sure to understand the question, but try to provide some useful pointers below. If the question is not answered by the pointers below, we would be happy to provide more.\n\n1. Comparison of Methodology.\nMethods for estimating the generalization performance of ridge and lasso regression focus on the minimizer (denoted as $\\hat{b}$ in eq. (2)) of the corresponding optimization problem. Such estimates, say for the Lasso, are studied in [1, 3, 9, 22] (citations from the submission PDF). These risk estimates are only valid for the final minimizer $\\hat b$ in eq. (2). \nIn particular, such estimators are not applicable to estimate the risk of intermediate iterates of algorithms.\n\nFor Lasso and Ridge, the risk estimates take the form of $\\frac{||y - X \\hat{b}||^2/n}{(1-df/n)^2}$, which can be viewed as an adjusted training error, where the degrees of freedom $df$ is, for the Lasso, the number of nonzero coefficients of $\\hat b$.\nIn contrast, our proposed risk estimator provides a consistent estimate of the risk of $\\hat{b}^t$ at each $t$.\nThe proposed risk estimate takes a very different form (weighted average of residual and previous gradients in eq. (11)) than the simpler ones available for Lasso/Ridge that can be simply described as an adjusted training error.\n\n2. Computational comparison. \nAs we mentioned in our response to Question 1 from Reviewer 3EXM, the computational complexity of calculating our estimate $\\tilde{r}_t$ for all $t\\in[T]$ is $O(npT^6)$. The computational complexity of the Lasso estimate is $O(np^2 + p^3)$ (see, for example, [1]). Thus, for a fixed $T$, our risk estimate has a lower complexity than that of Lasso. \n\n> (...) in terms of non-asymptotic statistical rates (related to consistency results)\n\nThe rate of convergence in Theorem 3.6 depends on the tails of the entries of the noise vector.\nIf the variance of the noise is finite then the right-hand side of (12) in Theorem 3.6 reduces to $C/(\\varepsilon\\sqrt n)$ (where $C$ has the same dependence as the constant in the numerator of (12)). Thus, if the variance is finite, we recover an error term of order $O_P(1/\\sqrt n)$, same as the best known rate for the difference between the risk and its estimate ([3] from the submission PDF).\n\nNow if the noise has infinite variance and is heavy-tailed, the rate of convergence of $\\min\\{1,||\\varepsilon||/n\\}$ in the right-hand side in (12) will depend on which moment of the noise exists (with a finite moment of order one, (12) converges to 0. If the moment of order $1+\\delta$ is finite, explicit rates can be obtained).\n\n> **Q4:** \nThe paper does not seem to make any assumptions on the penalty term, i.e. the function g in (2). Can the authors elaborate on this?\n\n**A4:**\nYes, we did not make explicit assumptions about the penalty term. Instead, we focus on the algorithm iterates that have the form of Eq. (4). The only requirement on the functions $\\phi_t,\\psi$ are given in Assumption 3.3. As long as the regression problem can be solved by algorithms similar to Eq. (4), our proposed risk estimate can be used to track the generalization error of these algorithms.\n\nFor instance, for any convex penalty function $g$, consider the corresponding proximal SGD algorithm in Example 2.4, with the soft-thresholding\nreplaced with the proximal of the penalty $g$. Since the proximal of $g$ is Lipschitz for any convex $g$, Assumption 3.3 is satisfied.\nNo further assumption is needed on the penalty $g$ except convexity.\n\n> **Q5:** \nCould 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?\n\n**A5:**\nThe consistency of the proposed estimate $\\hat r_t$ requires the regime where the ratio $p/n \\le \\gamma$, because in the proof (e.g., the proof of Lemma D.5), we need to bound the finite moments of $ ||X||\\_{op} /\\sqrt{n}$  by a constant. This is not possible if $p\\ggg n$\nas $||X||_{op}$ is known to be of order $\\sqrt n + \\sqrt p$  [11].\n\n---\n[1] Efron, Bradley, et al. \"Least angle regression.\" (2004): 407-499.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: 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.
Decision: Accept (poster)