, ,

# Uncertainty quantification for iterative algorithms in linear models with
application to early stopping

Pierre C<EMAIL_ADDRESS>[ Kai
<EMAIL_ADDRESS>[ Department of Statistics, Rutgers
Universitypresep=, ]e1,e2

###### Abstract

This paper investigates the iterates
${\widehat{\bm{b}}}^{1},\dots,{\widehat{\bm{b}}}^{T}$ obtained from iterative
algorithms in high-dimensional linear regression problems, in the regime where
the feature dimension $p$ is comparable with the sample size $n$, i.e.,
$p\asymp n$. The analysis and proposed estimators are applicable to Gradient
Descent (GD), proximal GD and their accelerated variants such as Fast
Iterative Soft-Thresholding (FISTA). The paper proposes novel estimators for
the generalization error of the iterate ${\widehat{\bm{b}}}^{t}$ for any fixed
iteration $t$ along the trajectory. These estimators are proved to be
$\sqrt{n}$-consistent under Gaussian designs. Applications to early-stopping
are provided: when the generalization error of the iterates is a U-shape
function of the iteration $t$, the estimates allow to select from the data an
iteration $\hat{t}$ that achieves the smallest generalization error along the
trajectory. Additionally, we provide a technique for developing debiasing
corrections and valid confidence intervals for the components of the true
coefficient vector from the iterate ${\widehat{\bm{b}}}^{t}$ at any finite
iteration $t$. Extensive simulations on synthetic data illustrate the
theoretical results.

high-dimensional models,

risk estimate,

Lasso,

generalization error,

###### keywords:

###### Contents

  1. 1 Introduction
    1. 1.1 Iterative algorithms
    2. 1.2 Statistical properties of iterates for small $t$ or when convergence fails
    3. 1.3 Early stopping
    4. 1.4 Contributions
    5. 1.5 Related literature
    6. 1.6 Notation
  2. 2 Main results
    1. 2.1 Iterates and derivatives
    2. 2.2 Lipschitz assumption and the chain rule
    3. 2.3 Memory matrix
    4. 2.4 Probabilistic assumptions and proportional regime
    5. 2.5 Estimation of the generalization error of iterates
    6. 2.6 An oracle inequality for early stopping
    7. 2.7 Alternative weights and covariance-adaptive weights
    8. 2.8 Iterative debiased estimation
  3. 3 Concrete examples
    1. 3.1 Gradient descent
    2. 3.2 Nesterov’s accelerated gradient descent (AGD)
    3. 3.3 Iterative Shrinkage-Thresholding Algorithm (ISTA)
    4. 3.4 Fast Iterative Shrinkage-Thresholding Algorithm (FISTA)
    5. 3.5 Local quadratic approximation in non-convex penalized regression
  4. 4 Efficient computation of the memory matrix
    1. 4.1 Iteratively computing rows of $\smash{\widehat{\mathbf{A}}}$
    2. 4.2 Hutchinson’s trace approximation
  5. 5 Discussion
  6. A Additional simulation results
    1. A.1 Under parametrization: $(n,p)=(1200,500)$
    2. A.2 Equal parametrization: $(n,p)=(1200,1200)$
  7. B Preliminary
    1. B.1 Notation
    2. B.2 Review of Kronecker product
  8. C Derivative formula
  9. D Change of variables
  10. E Useful intermediate results
    1. E.1 Deterministic operator norm bound
    2. E.2 Moment bounds
    3. E.3 Frobenius norm bounds on Jacobians
  11. F Proof of Theorem 2.2
    1. F.1 Proofs of Proposition F.1
    2. F.2 Proofs of Proposition F.2
    3. F.3 Proofs of supporting lemmas
  12. G Proof of Theorem 2.1
  13. H Proof of Corollary 2.3
  14. I Proof of Theorem 2.4
  15. J Proof of Corollary 2.5
  16. K Asymptotic normality and state evolution in the separable case

##  1 Introduction

Consider the linear model

(1.1) $\bm{y}=\bm{X}\bm{b}^{*}+{\bm{\varepsilon}},$

where $\bm{y}\in{\mathbb{R}}^{n}$ is the response vector,
$\bm{X}\in{\mathbb{R}}^{n\times p}$ is the design matrix with i.i.d. rows with
covariance $\bm{\Sigma}\in{\mathbb{R}}^{p\times p}$, and
${\bm{\varepsilon}}\in{\mathbb{R}}^{n}$ is the noise vector with i.i.d.
entries from $\mathsf{N}(0,\sigma^{2})$. To estimate the coefficient vector in
the linear model with dimension $p$ larger or comparable to the sample size
$n$, a common strategy is to consider the penalized least-squares estimator

(1.2) ${\widehat{\bm{b}}}\in\mathop{\rm
arg\,min}_{\bm{b}\in{\mathbb{R}}^{p}}\frac{1}{2n}\|\bm{y}-\bm{X}\bm{b}\|^{2}+\rho(\bm{b}).$

This estimator could be ordinary least-squares estimate (if $\rho(\bm{b})=0$),
Ridge regression, Lasso (Tibshirani, 1996), Elastic Net (Zou and Hastie,
2005), group Lasso (Yuan and Lin, 2006), nuclear norm penalization
(Koltchinskii et al., 2011) or Slope (Bogdan et al., 2015), to name a few
convex examples. Concave penalty functions $\rho$ in (1.2) have also been
considered, including SCAD (Fan and Li, 2001) and MCP (Zhang, 2010). The
properties of the penalized estimator ${\widehat{\bm{b}}}$ in (1.2) are now
well understood in high-dimensional regimes where the dimension $p$ is larger
than $n$. For instance, optimal L1 and L2 bounds of Lasso estimator
${\widehat{\bm{b}}}$ have been studied extensively (Wainwright, 2009; Bickel
et al., 2009; Ye and Zhang, 2010; Bayati and Montanari, 2012; Dalalyan et al.,
2017; Bellec et al., 2018, among others). Beyond risk bounds, confidence
intervals for single entries of $\bm{b}^{*}$ have been developed in Zhang and
Zhang (2014b); Javanmard and Montanari (2014); van de Geer et al. (2014a), and
more recently in Cai and Guo (2017); Javanmard and Montanari (2018); Bellec
and Zhang (2022b). We refer readers to Bühlmann and Van De Geer (2011); van de
Geer (2016) for a comprehensive review of such results. In most results cited
above, $p\gg n$ and the L2 risk of (1.2) typically converges to 0, i.e.,
${\widehat{\bm{b}}}$ is consistent. Consistency, on the other hand, is not
granted in the regime where dimension and sample size are of the same order,
i.e.,

(1.3) $p\asymp n,\qquad\text{ or }\qquad\lim(p/n)\text{ exists and is
finite}.$

In this regime, Bayati and Montanari (2012); Miolane and Montanari (2021)
showed that the Lasso L2 risk, under suitable assumptions, converges to a
finite positive constant that can be characterized by a system of nonlinear
equations. Regime (1.3) will be in force in the present paper. Since the early
works (Bayati and Montanari, 2012; El Karoui et al., 2013; Donoho and
Montanari, 2016) on this proportional regime, significant progress was made to
develop a unified theory explaining the behaviors of M-estimators and
penalized estimators when (1.3) holds, see for instance Thrampoulidis et al.
(2018); Celentano et al. (2023); Loureiro et al. (2021). Compared with the
early debiasing literature (Zhang and Zhang, 2014a; van de Geer et al., 2014a;
Javanmard and Montanari, 2014), new debiasing techniques are needed to
construct confidence intervals for components of $\bm{b}^{*}$ in this regime,
in particular with the need to account for the degrees-of-freedom of the
estimator (1.2) for the debiasing correction (Bellec and Zhang, 2022a;
Celentano et al., 2023; Bellec and Zhang, 2023).

###  1.1 Iterative algorithms

As we have seen, the properties of the estimator ${\widehat{\bm{b}}}$ in (1.2)
are now well understood, both in the classical high-dimensional regime with
$p\gg n$ and in the proportional regime (1.3). In practice, an exact minimizer
${\widehat{\bm{b}}}$ of (1.2) is typically unavailable in closed-form, and
numerical approximations are used instead. To solve the minimization problem
(1.2) numerically, the practitioner resorts to iterative algorithms developed
in the optimization literature to obtain an approximate solution to (1.2).

For smooth objective functions, this includes first-order methods such as
Gradient Descent (GD), its accelerated version known as Accelerated Gradient
Descent (AGD) (Nesterov, 1983), or its randomized variants (Bottou, 2010). For
non-smooth objective functions, the practitioner typically resorts to
coordinate descent (Friedman et al., 2010) or iterative algorithms involving
the proximal operator of $\rho$ such as Iterative Soft-Thresholding Algorithm
(ISTA) (Daubechies et al., 2004) and its accelerated version Fast Iterative
Soft-Thresholding Algorithm (FISTA) (Beck and Teboulle, 2009). Such iterative
algorithm starts with an initializer ${\widehat{\bm{b}}}^{1}$ and then
iteratively produces iterates
${\widehat{\bm{b}}}^{2},{\widehat{\bm{b}}}^{3},...$ by applying a nonlinear
function to a weighted sum of the previous iterate and the gradient of the
first term in (1.2) (the smooth part of the objective function) at the
previous iterate. For instance, for a convex penalty function $\rho$ in (1.2)
the proximal gradient descent algorithm solves (1.2) with the iterations

(1.4)
${\widehat{\bm{b}}}^{t+1}=\operatorname{prox}[\tfrac{1}{L}\rho]\bigl{(}{\widehat{\bm{b}}}^{t}+\tfrac{1}{nL}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t})\bigr{)}{,}$

where the proximal operator is defined as
$\operatorname{prox}[\frac{1}{L}\rho](\bm{v})=\mathop{\rm
arg\,min}_{\bm{b}\in{\mathbb{R}}^{p}}\|\bm{b}-\bm{v}\|^{2}/2+\frac{1}{L}\rho(\bm{b})$
and the scalar $L>0$ is usually taken (Beck and Teboulle, 2009) as an upper
bound on the Lipschitz constant of the gradient of the first term in the
objective function (1.2), i.e., such that $L\geq\|\bm{X}^{\top}\bm{X}/n\|_{\rm
op}$. For the Lasso problem with $\rho(\bm{b})=\lambda\|\bm{b}\|_{1}$ in
(1.2), the iterations (1.4) become the ISTA iterations

(1.5)
${\widehat{\bm{b}}}^{t+1}=\operatorname{soft}_{\lambda/L}\bigl{(}{\widehat{\bm{b}}}^{t}+\tfrac{1}{nL}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t})\bigr{)}{,}$

where $\operatorname{soft}_{\lambda/L}(\cdot)$ applies the soft-thresholding
$u\mapsto{\rm sign}(u)(|u|-\lambda/L)_{+}$ componentwise. For AGD and FISTA,
the recursion formula defining ${\widehat{\bm{b}}}^{t+1}$ involves the
iterates $({\widehat{\bm{b}}}^{t},{\widehat{\bm{b}}}^{t-1})$ as we will detail
in Sections 3.2 and 3.4 below.

More generally, the focus of the present paper is on iterative algorithm of
the form

(1.6)
${\widehat{\bm{b}}}^{t+1}=\bm{g}_{t+1}({\widehat{\bm{b}}}^{t},{\widehat{\bm{b}}}^{t-1},\bm{v}^{t},\bm{v}^{t-1}),\quad\text{
where}\quad\bm{v}^{t}=\tfrac{1}{n}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t})$

for nonlinear Lipschitz functions $\bm{g}_{t+1}$, in order to include
accelerated algorithms from the optimization literature that are of the form
(1.6), i.e., using the previous two iterations, including AGD (Nesterov, 1983)
and FISTA (Beck and Teboulle, 2009). Our theory actually applies to recursions
of the form

(1.7)
${\widehat{\bm{b}}}^{t+1}=\bm{g}_{t+1}({\widehat{\bm{b}}}^{t},{\widehat{\bm{b}}}^{t-1},\dots,{\widehat{\bm{b}}}^{2},{\widehat{\bm{b}}}^{1},~{}\bm{v}^{t},\bm{v}^{t-1},\dots,\bm{v}^{2},\bm{v}^{1}),$

although the typical optimization algorithms mentioned in the previous
paragraph only use the previous two iterates, as in (1.6).

These iterative algorithms are typically stopped when a stopping criterion is
met. Existing theories in the optimization literature on iterative algorithms
examine their convergence properties regarding the training error,
specifically how fast (in $t$) does the objective function at the iterate
${\widehat{\bm{b}}}^{t}$ approach the minimal value in (1.2), or how fast does
the iterate ${\widehat{\bm{b}}}^{t}$ approach an ideal minimizer
${\widehat{\bm{b}}}$ of (1.2) as $t\to+\infty$. For instance, Beck and
Teboulle (2009) showed that FISTA enjoys a faster rate of convergence compared
to ISTA.

###  1.2 Statistical properties of iterates for small $t$ or when convergence
fails

If ${\widehat{\bm{b}}}^{t}$ converges to ${\widehat{\bm{b}}}$ as $t$
increases, it is reasonable to expect that the properties of
${\widehat{\bm{b}}}$ (e.g., risk bounds, or debiased confidence intervals for
components of $\bm{b}^{*}$) are also valid for ${\widehat{\bm{b}}}^{t}$
provided $t$ is sufficiently large. One motivation of the present paper is the
study of iterates $({\widehat{\bm{b}}}^{t})_{t\geq 1}$ for small values of
$t$, or when convergence of these iterates fail for the optimization problem
(1.2). In such cases, because ${\widehat{\bm{b}}}^{t}$ and an ideal solution
${\widehat{\bm{b}}}$ to (1.2) significantly differ from each other, we should
not expect that the statistical properties of ${\widehat{\bm{b}}}$ are
inherited by ${\widehat{\bm{b}}}^{t}$. This calls for a statistical analysis
of ${\widehat{\bm{b}}}^{t}$ itself, rather than ideal minimizers
${\widehat{\bm{b}}}$ of the optimization problem (1.2). Two concrete
situations for which convergence fails are the following.

  * •

The dataset $(\bm{X},\bm{y})$ is so large that a single iteration of the
recursion (1.7) is slow and/or computationally costly. Running the iterations
(1.7) until convergence ($t\to+\infty$) is then not realistic. Still, in such
situation, the practitioner may only run a few iterations and use
${\widehat{\bm{b}}}^{t}$ for a small $t$ (e.g., $t=5$ or $t=10$) for
statistical purposes, including making predictions at new test data points
$\bm{x}_{new}$.

  * •

If the optimization problem (1.2) is non-convex, as is the case with SCAD (Fan
and Li, 2001) or MCP (Zhang, 2010), convergence of the iterates to a local or
global solution to (1.2) is a subtle problem. Iterative algorithms may
converge to a local optima. Specific algorithms have been developed to ensure
good statistical properties under suitable assumptions, such as Restricted
Eigenvalue conditions (Zhang, 2010; Feng and Zhang, 2019). However, such
conditions are often not verifiable in practice (Bandeira et al., 2013).
Despite the absence of any convergence guarantees, the algorithm used to solve
(1.2) for SCAD or MCP is still implementable, yielding iterates in the form of
(1.7). Despite the lack of verifiable conditions to guarantee convergence, at
a given $t$ the iterate ${\widehat{\bm{b}}}^{t}$ can still be used for
downstream tasks, for instance predictions at test data points $\bm{x}_{new}$.

In both situations, convergence does not occur and statistical properties of
${\widehat{\bm{b}}}^{t}$ at a finite $t$ are expected to significantly departs
from those of ${\widehat{\bm{b}}}$. Using ${\widehat{\bm{b}}}^{t}$ for
statistical purposes, instead of waiting for convergence, is also appealing in
order to save computational resources. This calls for a statistical theory of
${\widehat{\bm{b}}}^{t}$ itself, rather than ideal minimizers
${\widehat{\bm{b}}}$ of (1.2). The present makes contributions in this
direction with two major applications on the properties of the iterate
${\widehat{\bm{b}}}^{t}$ at a fixed $t$: assessing the predictive performance
of ${\widehat{\bm{b}}}^{t}$, and constructing debiased confidence intervals
for components of $\bm{b}^{*}$ based on ${\widehat{\bm{b}}}^{t}$.

###  1.3 Early stopping

It is important to recognize that the predictive performance of
${\widehat{\bm{b}}}^{t}$ does not necessarily improve as $t$ grows. In fact,
early-stopped gradient descent has been found to be effective in several
machine learning applications, including non-parametric regression in a
reproducing kernel Hilbert space (Raskutti et al., 2014), boosting Bühlmann
and Yu (2003); Rosset et al. (2004) and ridge regression (Ali et al., 2019),
to name a few. The intuition is that early stopping provides a form of
regularization that prevents overfitting (see, e.g., Ali et al. (2019) and the
references therein). In these settings, the population risk curve versus the
iteration number is U-shaped (first decreasing and then increasing); we will
observe instances of this phenomenon in Figures 1-5. For practical purposes,
the major question in such situations is to infer from the data if early
stopping is beneficial, and when to stop the iterations. The practitioner
typically wishes to stop the iteration as soon as the population risk starts
to increase; this is only possible if an estimator of the population risk of
${\widehat{\bm{b}}}^{t}$ is available. Our results below provide such an
estimator that can be effectively computed from the data $(\bm{X},\bm{y})$.

###  1.4 Contributions

The aforementioned situations motivate studying the properties of the iterates
${\widehat{\bm{b}}}^{t}$ at each step $t$, in particular for smaller values of
$t$ where ${\widehat{\bm{b}}}^{t}$ is not yet close to the optimizer
${\widehat{\bm{b}}}$. Our focus centers on two key questions:

  * (Q1)

Can we quantify the predictive performance of ${\widehat{\bm{b}}}^{t}$ for
each $t$ along the trajectory?

  * (Q2)

Can we use the iterates ${\widehat{\bm{b}}}^{t}$ to construct confidence
intervals for the entries of $\bm{b}^{*}$, for instance by deriving asymptotic
normality results around ${\widehat{\bm{b}}}^{t}$?

This paper provides affirmative answers to these questions. Our contributions
are summarized as follows:

  * i)

We introduce a reliable method for estimating the generalization error of each
iterate ${\widehat{\bm{b}}}^{t}$ obtained from widely-used algorithms that
solve (1.2). This estimator takes the form $\frac{1}{n}\|\sum_{s\leq
t}\hat{w}_{s,t}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s})\|^{2}$, i.e., the mean
square error of the weighted average of the residual vector
$\\{\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s}\\}_{s\leq t}$ with carefully chosen
data-driven weights $\hat{w}_{s,t}$ defined in Theorem 2.1 below. These
weights are algorithm-specific and are computed from the data. The theory also
allows us to derive estimators for the generalization error of weighted
averages of the iterates. For instance, our approach enables the estimation of
the generalization error for the average
$\frac{1}{T}\sum_{t=1}^{T}{\widehat{\bm{b}}}^{t}$. The corresponding theory is
presented in Theorem 2.2.

  * ii)

The proposed estimator of the generalization error serves as a practical proxy
for minimizing the true generalization error. To identify the iteration index
$t$ such that ${\widehat{\bm{b}}}^{t}$ has the lowest generalization error,
one can choose the $t$ that minimizes the estimated generalization error. This
approach offers a viable method for pinpointing the optimal stopping point in
the algorithm’s trajectory. This application is presented in Corollary 2.3.

  * iii)

We develop a method for constructing valid confidence intervals for the
elements of $\bm{b}^{*}$ based on a given iteration $t$ and iterate
${\widehat{\bm{b}}}^{t}$. Consequently, one needs not wait for the convergence
of the algorithm to construct valid confidence intervals for $\bm{b}^{*}$. In
cases where early-stopping is beneficial for achieving smaller generalization
error, confidence intervals based on the early-stopped iterate
${\widehat{\bm{b}}}^{\hat{t}}$ (where $\hat{t}$ is chosen by minimizing the
estimated risk) are shorter than the confidence interval based on the fully
converged minimizer ${\widehat{\bm{b}}}$. The theory on asymptotic normality
and confidence intervals is provided in Theorem 2.4.

  * iv)

The proposed method is applicable to a broad spectrum of regression techniques
and algorithms. We provide detailed expressions of the weights $\hat{w}_{s,t}$
for several popular algorithms in Section 3. It offers a practical solution
for determining the optimal stopping iteration of the algorithm to enjoy the
smallest generalization error as a function of $t$ along the trajectory.
Extensive simulation studies on synthetic data corroborate our theoretical
findings. These studies demonstrate the efficiency of our proposed risk
estimator for ${\widehat{\bm{b}}}^{t}$, and confirm the validity of the
asymptotic normality results for the entries of $\bm{b}^{*}$ derived from
${\widehat{\bm{b}}}^{t}$. The simulations are presented and discussed in of
Section 3.

###  1.5 Related literature

For regression problems in the proportional regime (1.3), the most extensively
studied iterative algorithm is Approximate Message Passing (AMP) (Donoho et
al., 2009; Bayati and Montanari, 2011, 2012). In the current framework, AMP
can be formulated as
${\widehat{\bm{b}}}^{t+1}=\bm{\eta}_{t+1}({\widehat{\bm{b}}}^{t}+\sum_{s\leq
t}w_{t,s}^{\textsc{amp}}\bm{v}^{s})$, where $w_{t,s}^{\textsc{amp}}$ are
specific scalar weights. These weights ensure that at each step, the input of
$\bm{\eta}_{t+1}$ is approximately normally distributed. This construction
enables tracking the evolution of the error of ${\widehat{\bm{b}}}^{t}$ as a
function of $t$ using a simple recursion, referred to as state evolution. AMP
has found numerous applications, including for instance in high-dimensional
statistics (Bayati and Montanari, 2012; Celentano and Montanari, 2022), spiked
models/low-rank matrix denoising (Montanari and Venkataramanan, 2021, and
references therein) or sampling (El Alaoui et al., 2022). A proof of the state
evolution for non-separable $\bm{\eta}_{t}$ is given in Berthier et al.
(2020). The focus of the present is to study algorithms from the optimization
literature such as proximal gradient (1.4), AGD or ISTA/FISTA that do not use
the special AMP weights.

Celentano et al. (2020); Montanari and Wu (2022) exhibit fundamental limits
for the performance of iterations that encompass (1.7). These works show that
iterates can be viewed as a Lipschitz function $g_{*}$ of some AMP algorithm.
They characterize the fundamental lower bound obtained by taking the infimum
over all $g_{*}$ and corresponding AMP algorithms, and show that the lower
bound is attained by a specific AMP algorithm, termed Bayes-AMP. Our goal in
this paper is different: to develop a data-driven estimate of the error of the
iterates ${\widehat{\bm{b}}}^{t}$.

Still in the proportional regime, Chandrasekher et al. (2022, 2023); Lou et
al. (2024) consider iterations that, in order to compute the $(t+1)$-th
iterate, use a fresh random batch $(\bm{x}_{i}^{t+1},y_{i}^{t+1})_{i\in
N^{t+1}}$ independent of all the past (here, the past refers to
${\widehat{\bm{b}}}^{t},{\widehat{\bm{b}}}^{t-1},...{\widehat{\bm{b}}}^{1}$ as
well as the random samples used to compute these). These works characterize
the evolution of the iterates’ performance (measured by mean square error, or
inner product with the ground truth) by a sequence of deterministic scalars
defined by simple recursions, and they show that the iterates’ performance is
close to the deterministic scalars. If samples are reused across iterations,
as in the present paper, it is still possible to characterize deterministic
equivalents of the performance of the iterates in some cases (see Celentano et
al. (2021) for continuous-time iterates, Gerbelot et al. (2022) for discrete-
time ones). The deterministic equivalents of the iterates’ performance are
necessarily more complex than in the fresh random batch setting. This
complexity is captured by limiting Gaussian processes indexed by $t\in[T]$
(Celentano et al., 2021; Gerbelot et al., 2022).

For optimization of the least-squares problem, with $\rho(\cdot)=0$ in (1.2),
Ali et al. (2020, 2019); Sheng and Ali (2022) (among others) study
equivalences between early stopping and L2 regularization: stopping the
iterations of GD or AGD at a fixed $t$ implicitly performs shrinkage of a
similar nature to Ridge regression. For gradient descent iterations solving
this least-squares problem, Patil et al. (2024) develops estimates of the
generalization error of the iterates using leave-one-out cross-validation, and
show consistency of these leave-one-out estimates along the GD trajectory.
While leave-one-out estimates are close to the goal of the present paper,
i.e., to develop estimates of the generalization error along the trajectory of
the algorithm, they are not computationally practical without further
modifications since they require to running $n$ algorithms in parallel to
obtain the leave-one-out estimate.

Outside of the proportional regime (1.3), Luo et al. (2023) develop an
efficient approximation of leave-one-out cross validation for iterates
${\widehat{\bm{b}}}^{t}$ solving empirical risk minimization problems with
additive regularization. The normalized Hessian of the objective function is
assumed to be Lipschitz and non-singular at the iterates, which precludes
settings of interest for the present paper with $p>n$ and non-smooth penalty
in (1.2). For uncertainty quantification of the iterates, Hoppe et al. (2023)
extends the classical debiased lasso estimator (Zhang and Zhang, 2014a;
Javanmard and Montanari, 2014; van de Geer et al., 2014b) to the iterates of a
variant of ISTA. The asymptotic normality in Hoppe et al. (2023) requires a
bound on the size of the supports of ${\widehat{\bm{b}}}^{t}$ and
$\bm{b}^{*}$, as well as $\lVert{\widehat{\bm{b}}}^{t}-\bm{b}^{*}\rVert$ to
vanish as $n\to\infty$. This is not the case in the regime (1.3) of interest
here where $p$ and $n$ are of the same order. The present paper departs from
this analysis of a variant ISTA on several fronts. First, our analysis is
grounded in a more general iteration formula (see (1.6) and (1.7)), which
encompasses a broader spectrum of algorithms. Second, our results establish
asymptotic normality of ${\widehat{\bm{b}}}^{t}$ for any $t$, and do not
require $\lVert{\widehat{\bm{b}}}^{t}-\bm{b}^{*}\rVert$ to vanish as
$n\to\infty$. In contrast to these studies, the present paper introduces new
estimator of the generalization error of algorithm iterates that solve the
regression problem (1.2) for $p>n$ and with non-smooth penalty functions.
Additionally, we propose a method to construct confidence intervals for the
entries of $\bm{b}^{*}$ by debiasing the iterate ${\widehat{\bm{b}}}^{t}$.

###  1.6 Notation

Regular variables like $a,b,...$ refer to scalars, bold lowercase letters such
as $\bm{a},\bm{b},...$ represent vectors, and bold uppercase letters like
$\bm{A},\bm{B},...$ indicate matrices. For an integer $n\geq 1$, we use the
notation $[n]=\\{1,\ldots,n\\}$. The vectors
$\bm{e}_{i}\in{\mathbb{R}}^{n},\bm{e}_{j}\in{\mathbb{R}}^{p},\bm{e}_{t}\in{\mathbb{R}}^{T}$
denote the canonical basis vector of the corresponding index. For clarity, we
always use implicit index $i$ to loop or sum over $[n]$, index $j$ to loop
over $[p]$, and indices $t,t^{\prime},s$ to loop over $[T]$. For a real vector
$\bm{a}\in{\mathbb{R}}^{p}$, $\left\lVert\bm{a}\right\rVert$ denotes its
Euclidean norm. For a matrix $\bm{A}$,
$\left\lVert\bm{A}\right\rVert_{\rm{F}}$,
$\left\lVert\bm{A}\right\rVert_{\rm{op}}$, $\left\lVert\bm{A}\right\rVert_{*}$
denote its Frobenius, operator and nuclear norm, respectively. For a matrix
$\bm{A}$, we denote its maximum and minimum singular values by
$\sigma_{\max}(\bm{A})$ and $\sigma_{\min}(\bm{A})$, respectively. If $\bm{A}$
is symmetric, $\lambda_{\max}(\bm{A})$ and $\lambda_{\min}(\bm{A})$ are its
maximum and minimum eigenvalues. Let $\bm{A}\otimes\bm{B}$ be the Kronecker
product of matrices $\bm{A}$ and $\bm{B}$. Let ${\bf 1}_{n}$ denote the all-
ones vector in ${\mathbb{R}}^{n}$, and $\bm{I}_{n}$ denote the identity matrix
of size $n$. For an event $E$, ${\mathbb{I}}(E)$ denotes the indicator
function of $E$. It takes the value 1 if the event $E$ occurs and 0 otherwise.
Let $\mathsf{N}(\mu,\sigma^{2})$ denote the Gaussian distribution with mean
$\mu$ and variance $\sigma^{2}$, and $\mathsf{N}_{k}(\bm{\mu},\bm{\Sigma})$
denote the $k$-dimensional Gaussian distribution with mean $\bm{\mu}$ and
covariance matrix $\bm{\Sigma}$. For a random sequence $\xi_{n}$, we write
$\xi_{n}=O_{P}(a_{n})$ if $\xi_{n}/a_{n}$ is stochastically bounded, and
${\smash{\xrightarrow{~{}\mathrm{p}~{}}}}$ for convergence in probability and
${\smash{\xrightarrow{~{}\mathrm{d}~{}}}}$ for convergence in distribution. We
reserve the letters $c$ and $C$ to denote generic constants. Additionally, we
use $C(\zeta,T,\gamma,\kappa)$ to denote a positive constant that depends only
on $\zeta,T,\gamma,\kappa$. The exact values of these constants may vary from
place to place.

##  2 Main results

Throughout the paper, given iterates
${\widehat{\bm{b}}}^{1},...,{\widehat{\bm{b}}}^{t},...$, define
$\bm{v}^{t}\in{\mathbb{R}}^{p}$ by

(2.1)
$\bm{v}^{t}=n^{-1}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t})\qquad\text{
for }t\geq 1.$

###  2.1 Iterates and derivatives

Our first task is to establish some notation (namely, matrices
${\mathcal{J}},{\mathcal{D}}\in{\mathbb{R}}^{pT\times pT}$) that will be used
to construct the proposed estimators. As a warm-up, we will first introduce
this notation for proximal gradient descent (Section 2.1.1) and iterates
constructed from the two previous ones (Section 2.1.2). Notation for the
general form (1.7) will be given in Section 2.1.3.

####  2.1.1 Proximal gradient descent

Consider the proximal gradient descent iterates (1.4). With this notation and
starting with ${\widehat{\bm{b}}}^{1}=\bf 0$, the proximal gradient descent
(1.4) iterations can be rewritten as

(2.2)
$\displaystyle{\widehat{\bm{b}}}^{t}=\bm{g}_{t}({\widehat{\bm{b}}}^{t-1},\bm{v}^{t-1}),\quad\text{
where
}\bm{g}_{t}(\bm{b}^{t},\bm{v}^{t})=\operatorname{prox}[\tfrac{1}{L}\rho](\bm{b}^{t}+\tfrac{1}{L}\bm{v}^{t}).$

The Jacobian of the function $\bm{g}_{t}$ is a matrix of size $p\times 2p$ and
can be partitioned into the two $p\times p$ blocks

(2.3)
$\bm{J}_{t,t-1}=\pdv{\bm{g}_{t}}{\bm{b}^{t-1}}({\widehat{\bm{b}}}^{t-1},\bm{v}^{t-1}),\quad\bm{D}_{t,t-1}=\pdv{\bm{g}_{t}}{\bm{v}^{t-1}}({\widehat{\bm{b}}}^{t-1},\bm{v}^{t-1}).$

Next, define ${\mathcal{J}}\in{\mathbb{R}}^{pT\times pT}$ and
${\mathcal{D}}\in{\mathbb{R}}^{pT\times pT}$ with $T\times T$ blocks of size
$p\times p$ as follows

(2.4)
${\mathcal{J}}=\begin{bmatrix}{\mathbf{0}}&{\mathbf{0}}&\cdots&\cdots&\cdots&{\mathbf{0}}\\\
\bm{J}_{2,1}&{\mathbf{0}}&\ddots&\cdots&\cdots&\vdots\\\
{\mathbf{0}}&\bm{J}_{3,2}&{\mathbf{0}}&\ddots&\cdots&\vdots\\\
{\mathbf{0}}&{\mathbf{0}}&\bm{J}_{4,3}&{\mathbf{0}}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&{\mathbf{0}}\\\
{\mathbf{0}}&\cdots&{\mathbf{0}}&{\mathbf{0}}&\bm{J}_{T,T-1}&{\mathbf{0}}\\\
\end{bmatrix},~{}{\mathcal{D}}=\begin{bmatrix}{\mathbf{0}}&{\mathbf{0}}&\cdots&\cdots&\cdots&{\mathbf{0}}\\\
\bm{D}_{2,1}&{\mathbf{0}}&\ddots&\cdots&\cdots&\vdots\\\
{\mathbf{0}}&\bm{D}_{3,2}&{\mathbf{0}}&\ddots&\cdots&\vdots\\\
{\mathbf{0}}&{\mathbf{0}}&\bm{D}_{4,3}&{\mathbf{0}}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&{\mathbf{0}}\\\
{\mathbf{0}}&\cdots&{\mathbf{0}}&{\mathbf{0}}&\bm{D}_{T,T-1}&{\mathbf{0}}\end{bmatrix}.$

Each $\mathbf{0}$ block in the above ${\mathcal{J}},{\mathcal{D}}$ is the zero
matrix of size $p\times p$. For iterations depending not only on
${\widehat{\bm{b}}}^{t-1},\bm{v}^{t-1}$ but also on previous iterates, the
lower triangular blocks of ${\mathcal{D}},{\mathcal{J}}$ will be filled with
the corresponding non-zero blocks of the Jacobian of $\bm{g}_{t}$, as shown in
the next subsections.

####  2.1.2 Combining two previous iterates

To incorporate more general algorithms beyond proximal gradient descent, we
consider the iteration formula:

(2.5)
$\displaystyle{\widehat{\bm{b}}}^{t}=\bm{g}_{t}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\bm{v}^{t-1},\bm{v}^{t-2})\qquad\text{
for }t\geq 2.$

Here, ${\widehat{\bm{b}}}^{1}$ serves as an initial vector, and we set
${\widehat{\bm{b}}}^{0}=\bm{v}^{0}={\mathbf{0}}_{p}$. The iteration functions
$\bm{g}_{t}(\cdot):{\mathbb{R}}^{4p}\to{\mathbb{R}}^{p}$ are determined by
user-specified algorithms. We include ${\widehat{\bm{b}}}^{t-2}$ and
$\bm{v}^{t-2}$ into the argument of (2.5) to encompass optimization algorithms
that update iterates by considering information from the two preceding steps.
This general formulation accommodates popular algorithms such as Nesterov’s
accelerated gradient (AGD) method (Nesterov, 1983, 2003) and the Fast
Iterative Shrinkage-Thresholding Algorithm (FISTA) (Beck and Teboulle, 2009).
We provide detailed expressions of $\bm{g}_{t}$ for several important
algorithms in Section 3.

For algorithms with iteration function $\bm{g}_{t}$ in (2.5) depending on the
previous two iterates, the Jacobian of $\bm{g}_{t}$ at a point where
$\bm{g}_{t}$ is differentiable is a matrix of size ${\mathbb{R}}^{p\times
4p}$. This Jacobian matrix in ${\mathbb{R}}^{p\times 4p}$ is naturally
partitioned into the four $p\times p$ blocks:

(2.6) $\displaystyle\bm{J}_{t,t-1}$
$\displaystyle=\pdv{\bm{g}_{t}}{\bm{b}^{t-1}}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\bm{v}^{t-1},\bm{v}^{t-2}),\quad\bm{D}_{t,t-1}=\pdv{\bm{g}_{t}}{\bm{v}^{t-1}}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\bm{v}^{t-1},\bm{v}^{t-2}),$
$\displaystyle\bm{J}_{t,t-2}$
$\displaystyle=\pdv{\bm{g}_{t}}{\bm{b}^{t-2}}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\bm{v}^{t-1},\bm{v}^{t-2}),\quad\bm{D}_{t,t-2}=\pdv{\bm{g}_{t}}{\bm{v}^{t-2}}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\bm{v}^{t-1},\bm{v}^{t-2}).$

We further define larger matrices
${\mathcal{J}},{\mathcal{D}}\in{\mathbb{R}}^{pT\times pT}$ by setting the
$(t,t^{\prime})$ block of ${\mathcal{J}}$ as $\bm{J}_{t,t^{\prime}}$ and that
of ${\mathcal{D}}$ as $\bm{D}_{t,t^{\prime}}$ for all $t,t^{\prime}\in[T]$ for
all $t,t^{\prime}$ corresponding to matrices in (2.6), with zeros everywhere
else:

(2.7)
${\mathcal{J}}=\begin{bmatrix}{\mathbf{0}}&{\mathbf{0}}&\cdots&\cdots&\cdots&{\mathbf{0}}\\\
\bm{J}_{2,1}&{\mathbf{0}}&\ddots&\cdots&\cdots&\vdots\\\
\bm{J}_{3,1}&\bm{J}_{3,2}&{\mathbf{0}}&\ddots&\cdots&\vdots\\\
{\mathbf{0}}&\bm{J}_{4,2}&\bm{J}_{4,3}&{\mathbf{0}}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&{\mathbf{0}}\\\
{\mathbf{0}}&\cdots&{\mathbf{0}}&\bm{J}_{T,T-2}&\bm{J}_{T,T-1}&{\mathbf{0}}\\\
\end{bmatrix},~{}{\mathcal{D}}=\begin{bmatrix}{\mathbf{0}}&{\mathbf{0}}&\cdots&\cdots&\cdots&{\mathbf{0}}\\\
\bm{D}_{2,1}&{\mathbf{0}}&\ddots&\cdots&\cdots&\vdots\\\
\bm{D}_{3,1}&\bm{D}_{3,2}&{\mathbf{0}}&\ddots&\cdots&\vdots\\\
{\mathbf{0}}&\bm{D}_{4,2}&\bm{D}_{4,3}&{\mathbf{0}}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&{\mathbf{0}}\\\
{\mathbf{0}}&\cdots&{\mathbf{0}}&\bm{D}_{T,T-2}&\bm{D}_{T,T-1}&{\mathbf{0}}\end{bmatrix}.$

An equivalent definition of ${\mathcal{J}},{\mathcal{D}}$ using Kronecker
products is

${\mathcal{J}}=\sum_{t=2}^{T}\sum_{s=t-2}^{t-1}(\bm{e}_{t}\bm{e}_{s}^{\top})\otimes\bm{J}_{t,s},\qquad{\mathcal{D}}=\sum_{t=2}^{T}\sum_{s=t-2}^{t-1}(\bm{e}_{t}\bm{e}_{s}^{\top})\otimes\bm{D}_{t,s},$

where $\bm{e}_{t},\bm{e}_{s}$ are the $t$-th and $s$-th canonical basis
vectors in ${\mathbb{R}}^{T}$.

####  2.1.3 General form: combining all previous iterates

More generally, our theory applies to iterations that depend on all previous
iterates, i.e., iterates of the form

(2.8)
${\widehat{\bm{b}}}^{t}=\bm{g}_{t}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\dots,{\widehat{\bm{b}}}^{2},{\widehat{\bm{b}}}^{1},~{}\bm{v}^{t-1},,\dots,\bm{v}^{2},\bm{v}^{1})$

with $\bm{v}^{s}={\widehat{\bm{b}}}^{s}=\mathbf{0}_{p}$ for all $s\leq 0$ by
convention. With the more general iterations (2.8) depending on all previous
iterates, the Jacobian of $\bm{g}_{t}$ in (2.8) is now a matrix in
${\mathbb{R}}^{p\times 2p(t-1)}$ partitioned into blocks

(2.9)
$\bm{J}_{t,s}=\frac{\partial\bm{g}_{t}}{\partial\bm{b}^{s}}\Bigl{(}{\widehat{\bm{b}}}^{t-1},...,{\widehat{\bm{b}}}^{1},\bm{v}^{t-1},...,\bm{v}^{1}\Bigr{)},\quad\bm{D}_{t,s}=\frac{\partial\bm{g}_{t}}{\partial\bm{v}^{s}}\Bigl{(}{\widehat{\bm{b}}}^{t-1},...,{\widehat{\bm{b}}}^{1},\bm{v}^{t-1},...,\bm{v}^{1}\Bigr{)}$

for all $s=1,...,t-1$. The large matrices
${\mathcal{J}},{\mathcal{D}}\in{\mathbb{R}}^{pT\times pT}$ are then defined as

(2.10)
${\mathcal{J}}=\sum_{t=2}^{T}\sum_{s=1}^{t-1}(\bm{e}_{t}\bm{e}_{s}^{\top})\otimes\bm{J}_{t,s},\qquad{\mathcal{D}}=\sum_{t=2}^{T}\sum_{s=1}^{t-1}(\bm{e}_{t}\bm{e}_{s}^{\top})\otimes\bm{D}_{t,s}.$

This generalizes (2.7) when ${\widehat{\bm{b}}}^{t}$ depends on all previous
iterates. The matrices ${\mathcal{J}}$ and ${\mathcal{D}}$ are always lower
triangular block matrices, as in (2.7).

###  2.2 Lipschitz assumption and the chain rule

Throughout, we assume the following.

######  Assumption 2.1.

The algorithm starts with ${\widehat{\bm{b}}}^{1}={\mathbf{0}}$. The iteration
function $\bm{g}_{t}$ in (2.8) is $\zeta$-Lipschitz and satisfies
$\bm{g}_{t}({\mathbf{0}})={\mathbf{0}}$.

The notation in (2.6) and (2.9) naturally arises from the application of the
chain rule in the recursive computation of derivatives for the iterates (2.5)
and (2.8). While assuming that $\bm{g}_{t}$ is Lipschitz (Assumption 2.1)
implies that the iterates are differentiable almost everywhere with respect to
$(\bm{X},\bm{y})$ by Rademacher’s theorem (a locally Lipschitz function of
$(\bm{y},\bm{X})$ is differentiable almost everywhere), the chain rule may
fail for composition of multivariate Lipschitz functions. This technical issue
is clearly apparent with the example

(2.11)
$\phi^{*}(u,v)=(u,u),\quad\phi^{**}(u,v)=\max(u,{v}),\quad\phi^{**}\circ\phi^{*}(u,v)=u.$

Although $\frac{\partial}{\partial
u}\bigl{(}\phi^{**}\circ\phi^{*}\bigr{)}(u,v)=1$ clearly holds, the chain rule
is undefined because $\phi^{**}$ is not differentiable at $\phi^{*}(u,v)$ for
all $(u,v)\in{\mathbb{R}}^{2}$. A fix to this issue for the composition of
locally Lipschitz functions is given in Corollary 3.2 of Ambrosio and Dal Maso
(1990): the chain rule for for the composition of two Lipschitz functions
$g:{\mathbb{R}}^{k}\to{\mathbb{R}}^{q},u:{\mathbb{R}}^{q}\to{\mathbb{R}}^{p}$
holds almost surely in the modified form $\nabla(g\circ
u)(x)=\nabla(g|_{T_{x}^{u}})\nabla u(x)$ where $g|_{T_{x}^{u}}$ is the
restriction of $g$ to the affine space defined in (Ambrosio and Dal Maso,
1990, Corollary 3.2). Thus, in order to leverage the chain rule in our proofs,
we will assume that the matrices in (2.6) or (2.9) are all modified, if
necessary, in order to grant the chain rule as given in (Ambrosio and Dal
Maso, 1990, Corollary 3.2). In many practical examples, the Lipschitz (but
not-necessarily differentiable everywhere) nonlinear function used to
recursively obtain the iterates is separable, i.e.,
$\bm{g}^{sep}:{\mathbb{R}}^{p}\to{\mathbb{R}}^{p}$ of the form

$\bm{g}^{sep}(\bm{x})_{j}=g_{j}(x_{j})\quad\text{ for all
}\bm{x}\in{\mathbb{R}}^{p},j\in[p]\quad\text{for some functions
}g_{j}:{\mathbb{R}}\to{\mathbb{R}}.$

This is the case for the soft-thresholding in (1.5) as well as its variants,
FISTA and LQA, detailed in Sections 3.4 and 3.5. In this case by (Ziemer,
1989, Theorem 2.1.11), the chain rule holds almost everywhere for the
composition of Lipschitz functions of the form $g_{j}\circ h$ for
$h:{\mathbb{R}}^{k}\to{\mathbb{R}}$ and $g_{j}:{\mathbb{R}}\to{\mathbb{R}}$.
Consequently, the usual (unmodified) chain rule is granted almost everywhere
for compositions of Lipschitz functions that are either separable or
everywhere differentiable.

###  2.3 Memory matrix

Equipped with the above notation in (2.9) and (2.10), we now introduce the
memory matrix. The memory matrix is a key ingredient in the construction of
our estimator of the prediction error of the iterate ${\widehat{\bm{b}}}^{t}$
obtained from the iterations (2.5) or (2.8).

######  Definition 2.1.

For iterates
$\\{{\widehat{\bm{b}}}^{1},{\widehat{\bm{b}}}^{2},...,{\widehat{\bm{b}}}^{T}\\}$
obtained from the iteration (2.8), with ${\mathcal{D}},{\mathcal{J}}$ defined
in (2.10), define the memory matrix
$\smash{\widehat{\mathbf{A}}}\in{\mathbb{R}}^{T\times T}$ as

(2.12)
$\smash{\widehat{\mathbf{A}}}=\frac{1}{n}\sum_{i=1}^{n}\Bigl{(}\bm{I}_{T}\otimes(\bm{e}_{i}^{\top}\bm{X})\Bigr{)}\Bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\Bigr{)}^{-1}{\mathcal{D}}\Bigl{(}\bm{I}_{T}\otimes(\bm{X}^{\top}\bm{e}_{i})\Bigr{)}$

where $\bm{e}_{i}\in{\mathbb{R}}^{n}$ is the $i$-th canonical basis vector.

By definition, the memory matrix $\smash{\widehat{\mathbf{A}}}$ only depends
on the data $(\bm{X},\bm{y})$, the iterates ${\widehat{\bm{b}}}^{t}$, the
vectors $\bm{v}^{t}$ in (2.1), and the derivative matrices in (2.9) evaluated
at
$({\widehat{\bm{b}}}^{t-1},...,{\widehat{\bm{b}}}^{1},\bm{v}^{t-1},...,\bm{v}^{1})$.
The specific computation of $\smash{\widehat{\mathbf{A}}}$ thus depends on the
iteration function $\bm{g}_{t}$ through the derivative matrices in (2.9) that
appear as blocks of ${\mathcal{D}}$ and ${\mathcal{J}}$. The matrix
$\smash{\widehat{\mathbf{A}}}$ will differ significantly for different choices
of nonlinear functions $\bm{g}_{t}$. For the same nonlinear functions
$\bm{g}_{t}$, the matrix $\smash{\widehat{\mathbf{A}}}$ will also be different
for different realizations of the data $(\bm{X},\bm{y})$, although we observe
in practice that the entries of $\smash{\widehat{\mathbf{A}}}$ concentrate and
have small variance with respect to the randomness of $(\bm{X},\bm{y})$.

In Section 3, we provide detailed expressions of ${\mathcal{J}}$ and
${\mathcal{D}}$ for several important algorithms in the optimization
literature. While definition (2.12) involves the inverse of a matrix of size
${\mathbb{R}}^{pT\times pT}$, due to its specific structure, computation of
$\smash{\widehat{\mathbf{A}}}$ requires much fewer resources that performing a
matrix inversion or solving a linear system in ${\mathbb{R}}^{pT\times pT}$.
We provide efficient methods to compute $\smash{\widehat{\mathbf{A}}}$ row by
row in Section 4.

An equivalent definition of the matrix $\smash{\widehat{\mathbf{A}}}$ is to
consider the matrix

(2.13)
$\frac{1}{n}\Bigl{(}\bm{I}_{T}\otimes\bm{X}\Bigr{)}\Bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\Bigr{)}^{-1}{\mathcal{D}}\Bigl{(}\bm{I}_{T}\otimes(\bm{X}^{\top})\Bigr{)}$

in ${\mathbb{R}}^{nT\times nT}$ as a matrix by block, with $T\times T$ blocks
of size $n\times n$. The matrix
$\smash{\widehat{\mathbf{A}}}\in{\mathbb{R}}^{T\times T}$ is then obtained by
taking the trace of each $n\times n$ block.

Since ${\mathcal{D}}$ and ${\mathcal{J}}$ are lower triangular matrices with
diagonal blocks equal to $\mathbf{0}_{p\times p}$, both
$\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\bm{X}^{\top}\bm{X}/n)-{\mathcal{J}}$
and its inverse are lower triangular with diagonal blocks equal to
$\bm{I}_{p}$. Consequently, (2.13) is block lower triangular with diagonal
blocks equal to $\mathbf{0}_{n\times n}$. By taking the trace of each $n\times
n$ block in (2.13), we see that
$\smash{\widehat{\mathbf{A}}}\in{\mathbb{R}}^{T\times T}$ is lower triangular
with zero diagonal entries. It thus always holds that the matrix
$\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n$ is a lower triangular with all
diagonal entries equal to 1, invertible, with
$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$ also lower triangular with
all diagonal entries equal to 1. The matrix

$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$

plays a critical role in our estimator of the generalization error of each
iteration ${\widehat{\bm{b}}}^{t}$ in Theorems 2.1 and 2.2 below.

###  2.4 Probabilistic assumptions and proportional regime

Our main theorems below hold under the following assumptions.

######  Assumption 2.2.

The design matrix $\bm{X}$ has i.i.d. rows from $\mathsf{N}_{p}(\bf
0,\bm{\Sigma})$ for some positive definite matrix $\bm{\Sigma}$ satisfying
$0<\lambda_{\min}(\bm{\Sigma})\leq 1\leq\lambda_{\max}(\bm{\Sigma})$ and
$\lVert\bm{\Sigma}\rVert_{\rm{op}}\|\bm{\Sigma}^{-1}\|_{\rm op}\leq\kappa$.

######  Assumption 2.3.

The noise ${\bm{\varepsilon}}$ is independent of $\bm{X}$ and has i.i.d.
entries from $\mathsf{N}(0,\sigma^{2})$.

######  Assumption 2.4.

The sample size $n$ and predictor dimension $p$ satisfy $p/n\leq\gamma$ for a
constant $\gamma\in(0,\infty)$.

Recall that $\zeta$ is the Lipschitz constant in Assumption 2.1. In the
results below, $(\zeta,T,\gamma,\kappa)$ can be thought of as constant problem
parameters as $n,p\to+\infty$; our upper bounds involve constants
$C(\zeta,T,\gamma,\kappa)$ that only depend on these four quantities.

###  2.5 Estimation of the generalization error of iterates

We are interested in quantifying the performance of each iterate
${\widehat{\bm{b}}}^{t}$. We assess the performance of the iterate
${\widehat{\bm{b}}}^{t}$ using the out-of-sample prediction error (also
referred to as prediction risk or generalization error)

(2.14) $r_{t}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}{\mathbb{E}}\Bigl{[}\bigl{(}y_{new}-\bm{x}_{new}^{\top}{\widehat{\bm{b}}}^{t}\bigr{)}^{2}\mid(\bm{X},\bm{y})\Bigr{]}=\|\bm{\Sigma}^{1/2}({\widehat{\bm{b}}}^{t}-\bm{b}^{*})\|^{2}+\sigma^{2},$

where $(\bm{x}_{new},y_{new})$ is a new test sample that has the same
distribution as $(\bm{x}_{1},y_{1})$, the first row of $(\bm{X},\bm{y})$.
Since ${\widehat{\bm{b}}}^{t}$ is measurable with respect to
$(\bm{X},\bm{y})$, the last equality above follows from Assumptions 2.2 and
2.3. Note that $(\bm{x}_{new},y_{new})$ is only used as a mathematical device
to define $r_{t}$ in (2.14); we do not assume that independent test samples
are available. Our first result concerns the estimation of prediction risk
$r_{t}$ for each iterate ${\widehat{\bm{b}}}^{t}$.

######  Theorem 2.1 (Estimation of prediction risk).

Let Assumptions 2.2, 2.3, 2.4 and 2.1 be fulfilled. For each $t\in[T]$, define
the estimate $\hat{r}_{t}$ of $r_{t}$ by

(2.15)
$\hat{r}_{t}=\frac{1}{n}\Big{\|}\sum_{s=1}^{t}\hat{w}_{t,s}\bigl{(}\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s}\bigr{)}\Big{\|}^{2},$

where
$\hat{w}_{t,s}=\bm{e}_{t}^{\top}(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}\bm{e}_{s}$.
We have for any $t\in[T]$,

(2.16) ${\mathbb{E}}\big{|}\hat{r}_{t}-r_{t}\big{|}\leq n^{-\frac{1}{2}}{\rm
var}(y_{1})C(\zeta,T,\gamma,\kappa).$

Here ${\rm
var}(y_{1})=\lVert\bm{\Sigma}^{1/2}\bm{b}^{*}\rVert^{2}+\sigma^{2}$, and
$C(\zeta,T,\gamma,\kappa)$ is a constant depending only on
$\zeta,T,\gamma,\kappa$.

Theorem 2.1 establishes that $\hat{r}_{t}$ is $\sqrt{n}$-consistent when ${\rm
var}(y_{1}),\zeta,T,\gamma,\kappa$ are constant as $n,p\to+\infty$. The risk
estimate $\hat{r}_{t}$ for ${\widehat{\bm{b}}}^{t}$ is determined by both the
current residual vector $\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t}$ and all
preceding residual vectors
$\\{\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s}\\}_{s=1,...,t-1}$ through the weighted
sum inside the mean squared norm in (2.15). These weights
$(\hat{w}_{t,1},\ldots,\hat{w}_{t,t})$ are the first $t$ entries of the $t$-th
row of the matrix $(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$. The
proof of Theorem 2.1 is given in Appendix G.

By construction of the matrices $\smash{\widehat{\mathbf{A}}}$ and
$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$, the definition of these
weights are not influenced by the total number of iterations $T$ (e.g., the
weights $(\hat{w}_{s,t})_{s\leq t}$ for any $t\leq T$ are the same if the
total number of iterations is $T$ or any $T^{\prime}>T$). In other words,
opting for a larger $T$ simply enlarges the dimensions of the matrix
$\smash{\widehat{\mathbf{A}}}$ and the inverse matrix
$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$, yet the first $t$ entries
in the $t$-th row of $(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$ remain
unchanged. The necessity of a total number of iteration, $T$, is arguably
artificial; an equivalent presentation would increase the size of
$\smash{\widehat{\mathbf{A}}}$ and
$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$ by one row and one column
at every new iteration. We opted for a fixed total number of iterations $T$
with matrices $\smash{\widehat{\mathbf{A}}}$ and
$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$ having fixed size $T\times
T$ for clarity, to avoid dealing with matrices changing in size.

######  Remark 2.1 (Risk of initialization).

For $t=1$, $\hat{w}_{t,s}={\mathbb{I}}(s=1)$ since
$(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}$ is lower triangular with
diagonal entries all equal to 1. As a sanity check, we obtain
$\hat{r}_{1}=n^{-1}\|\bm{y}-\bm{X}{\widehat{\bm{b}}}^{1}\|^{2}$ from (2.15),
which is an unbiased estimate of the prediction error of the initialization
${\widehat{\bm{b}}}^{1}$ since ${\widehat{\bm{b}}}^{1}$ is independent of
$(\bm{X},\bm{y})$.

Besides the individual iterate ${\widehat{\bm{b}}}^{t}$, one might also
consider the estimate $\bar{\bm{b}}$ defined as the average of $m$ consecutive
iterates starting from ${\widehat{\bm{b}}}^{t_{0}}$, i.e.,
$\bar{\bm{b}}=\frac{1}{m}\sum_{t=t_{0}}^{t_{0}+m-1}{\widehat{\bm{b}}}^{t}$.
The prediction risk of the estimate $\bar{\bm{b}}$ is then

$\sigma^{2}+\|\bm{\Sigma}^{1/2}(\bar{\bm{b}}-\bm{b}^{*})\|^{2}=\frac{1}{m^{2}}\sum_{t=t_{0}}^{t_{0}+m-1}\sum_{t^{\prime}=t_{0}}^{t_{0}+m-1}\Bigl{[}\sigma^{2}+({\widehat{\bm{b}}}^{t}-\bm{b}^{*})^{\top}\bm{\Sigma}({\widehat{\bm{b}}}^{t^{\prime}}-\bm{b}^{*})\Bigr{]},$

simply by expanding the square. Estimation of the generalization error of
$\bar{\bm{b}}$ is thus possible provided that we can estimate the terms

(2.17) $r_{tt^{\prime}}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\sigma^{2}+({\widehat{\bm{b}}}^{t}-\bm{b}^{*})^{\top}\bm{\Sigma}({\widehat{\bm{b}}}^{t^{\prime}}-\bm{b}^{*})$

for each $t,t^{\prime}$ inside the double sum. The following result derives a
$\sqrt{n}$-consistent estimate of $r_{tt^{\prime}}$, using the same weights as
in Theorem 2.1.

######  Theorem 2.2 (Proof is given in Appendix F).

Let Assumptions 2.2, 2.3, 2.4 and 2.1 be fulfilled. For two integers
$t,t^{\prime}\leq T$ define the estimate of $r_{tt^{\prime}}$ as

(2.18)
$\hat{r}_{tt^{\prime}}=\frac{1}{n}\Bigl{(}\sum_{s=1}^{t}\hat{w}_{t,s}\bigl{(}\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s}\bigr{)}\Bigr{)}^{\top}\Bigl{(}\sum_{s^{\prime}=1}^{t^{\prime}}\hat{w}_{t^{\prime},s^{\prime}}\bigl{(}\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s^{\prime}}\bigr{)}\Bigr{)},$

where
$\hat{w}_{a,b}=\bm{e}_{a}^{\top}(\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}/n)^{-1}\bm{e}_{b}$
for all $a,b\in[T]$ as in Theorem 2.1. Then

$\displaystyle{\mathbb{E}}|\hat{r}_{tt^{\prime}}-r_{tt^{\prime}}|+{\rm
var}(r_{tt^{\prime}})^{1/2}\leq n^{-\frac{1}{2}}C(\zeta,T,\gamma,\kappa){\rm
var}(y_{1}).$

###  2.6 An oracle inequality for early stopping

Since Theorem 2.1 provides upper bounds of the form
${\mathbb{E}}|\hat{r}_{t}-r_{t}|\leq n^{-1/2}C(\zeta,T,\gamma,\kappa){\rm
var}(y_{1})$. A straightforward application of Markov’s inequality yields the
following.

######  Corollary 2.3 (Proof is given in Appendix H).

Let Assumptions 2.2, 2.3, 2.4 and 2.1 be fulfilled. Select an iteration
$\hat{t}\in[T]$ by minimizing (2.15), that is, $\hat{t}=\mathop{\rm
arg\,min}_{t\in[T]}\|\sum_{s\leq
t}\hat{w}_{t,s}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s})\|^{2}/n$. Then for any
constant $c\in(0,1/2)$,

${\mathbb{P}}\Bigl{(}\|\bm{\Sigma}^{1/2}({\widehat{\bm{b}}}^{{\hat{t}}}-\bm{b}^{*})\|^{2}\leq\min_{s\in[T]}\|\bm{\Sigma}^{1/2}({\widehat{\bm{b}}}^{s}-\bm{b}^{*})\|^{2}+\frac{{\rm
var}(y_{1})}{n^{1/2-c}}\Bigr{)}\geq
1-\frac{C(\zeta,\gamma,T,\kappa)}{n^{c}}\to 1.$

This means that after running $T$ iterations, we can pick the iteration
$\hat{t}$ that minimizes the criterion (2.15), and this choice achieves the
smallest prediction error among the first $T$ iterates up to a negligible
error term. This is appealing for the settings illustrated in Figures 1, 2, 3,
4 and 5 where the generalization error of the iterates
${\widehat{\bm{b}}}^{t}$ is first decreasing in $t$ up to some $t_{*}$, before
increasing for $t\geq t_{*}$. The above selection of $\hat{t}$ guarantees that
the risk of ${\widehat{\bm{b}}}^{\hat{t}}$ is close to that of
${\widehat{\bm{b}}}^{t_{*}}$.

Given $K_{n}$ candidate sequences of nonlinear functions, say
$(\bm{g}_{t}^{k})_{t\in[T]}$ for each $k\in[K_{n}]$, a similar argument using
Markov’s inequality yields

$\|\bm{\Sigma}^{1/2}({\widehat{\bm{b}}}^{{\hat{t}}}-\bm{b}^{*})\|^{2}\leq\min_{k\in[K_{n}]}\min_{s\in[T]}\|\bm{\Sigma}^{1/2}({\widehat{\bm{b}}}^{s}-\bm{b}^{*})\|^{2}+{\rm
var}(y_{1})o_{P}(1)$

if $(\zeta,\gamma,T,\kappa)$ are constants as $n,p\to+\infty$ and
$K_{n}=o(\sqrt{n})$.

###  2.7 Alternative weights and covariance-adaptive weights

The proofs of Theorems 2.1 and 2.2 reveal that if $\bm{\Sigma}=\bm{I}_{p}$,
the weights $\hat{w}_{t,s}$ in (2.15) and (2.18) can be replaced with the
alternative weights

$\begin{aligned}
\check{w}_{t,s}&=\mathbb{I}\\{t=s\\}+{\frac{1}{n}}\operatorname{Tr}\Bigl{[}(\bm{e}_{t}^{\top}\otimes\bm{I}_{p})\bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\bigr{)}^{-1}{\mathcal{D}}(\bm{e}_{s}\otimes\bm{I}_{p})\Bigr{]}\\\
&=\mathbb{I}\\{t=s\\}+{\frac{1}{n}}\sum_{j=1}^{p}(\bm{e}_{t}^{\top}\otimes\bm{e}_{j}^{\top})\bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\bigr{)}^{-1}{\mathcal{D}}(\bm{e}_{s}\otimes\bm{e}_{j})\end{aligned}\text{
(for $\bm{\Sigma}=\bm{I}_{p}$).}$

Indeed, these alternative weights $\check{w}_{t,s}$ are the $(t,s)$ entries of
the matrix $\bm{I}_{T}+{\widehat{\bm{C}}}/n$ in the proof of Theorem 2.2 with
${\widehat{\bm{C}}}\in{\mathbb{R}}^{T\times T}$ defined in (F.2). These
weights satisfy

(2.19)
${\mathbb{E}}\Bigl{[}\sum_{s=1}^{T}\frac{1}{n}\Big{\|}\sum_{s=1}^{t}(\hat{w}_{t,s}-\check{w}_{t,s})(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s})\Big{\|}^{2}\Bigr{]}\leq
C(\zeta,T,\gamma)n^{-1/2},\qquad\text{ when }\bm{\Sigma}=\bm{I}_{p},$

by (F.10) in the appendix. This means that $\hat{w}_{t,s}$ and
$\check{w}_{t,s}$ can be used interchangeably as weights for
$(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s})_{s\leq T}$ in the definition of
$\hat{r}_{t}$ when $\bm{\Sigma}=\bm{I}_{p}$.

This interchangeability is lost as soon as $\bm{\Sigma}\neq\bm{I}_{p}$. In
this case, the alternative weights $\check{w}_{t,s}$ have expression

(2.20)
$\check{w}_{t,s}=\mathbb{I}\\{t=s\\}+{\tfrac{1}{n}}\operatorname{Tr}\bigl{[}(\bm{e}_{t}^{\top}\otimes\bm{\Sigma}^{1/2})\bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\bigr{)}^{-1}{\mathcal{D}}(\bm{e}_{s}\otimes\bm{\Sigma}^{1/2})\bigr{]}$

using the correspondence (D.6). The alternative weights $\check{w}_{t,s}$ thus
require the knowledge of $\bm{\Sigma}$. This is the main reason we presented
the main results using the weights $\hat{w}_{t,s}$ and the memory matrix
$\smash{\widehat{\mathbf{A}}}$: The expressions of $\hat{w}_{t,s}$ and
$\smash{\widehat{\mathbf{A}}}$ do not require prior knowledge of $\bm{\Sigma}$
or estimating $\bm{\Sigma}$ from the data. The weights $\hat{w}_{t,s}$ are
thus preferable, and more broadly applicable than $\check{w}_{t,s}$. This
duality between two interchangeable scalar adjustments, one requiring the
knowledge of $\bm{\Sigma}$ and one not requiring it, was already observed in
regularized M-estimation (Bellec and Shen, 2022, Section 5).

###  2.8 Iterative debiased estimation

This subsection focuses on asymptotic normality and statistical inference for
the entries of $\bm{b}^{*}$ using iterate ${\widehat{\bm{b}}}^{t}$. The next
theorem shows that the following debiased estimate,

(2.21)
${\widehat{\bm{b}}}^{t,\rm{debias}}_{j}:=\underbrace{{\widehat{\bm{b}}}^{t}_{j}}_{\text{iterate}}~{}+~{}\underbrace{n^{-1}\sum_{s=1}^{t}\hat{w}_{t,s}\bigl{(}\bm{y}-\bm{X}{\widehat{\bm{b}}}^{s}\bigr{)}^{\top}\bm{X}\bm{\Sigma}^{-1}\bm{e}_{j}}_{\text{bias
correction}},$

is approximately normally distributed and centered at $\bm{b}_{j}^{*}$. Above,
the observable weights $\hat{w}_{t,s}$ are the same as in Theorem 2.1. The
approximate normality below is quantified using the 2-Wasserstein distance
between probability distributions in ${\mathbb{R}}^{T}$, defined as
$W_{2}(\mu,\nu)=\inf_{(\bm{u},\bm{w})}{\mathbb{E}}[\|\bm{u}-\bm{v}\|^{2}]^{1/2}$
where the infimum is taken over all couplings $(\bm{u},\bm{w})$ of the two
probability measures $(\mu,\nu)$. We refer to (Villani et al., 2009,
Definition 6.8) for several characterizations of convergence with respect to
$W_{2}$.

The following result is understood asymptotically as $n,p\to+\infty$.
Implicitly, we assume that a sequence of regression problems (1.1) is given,
together with nonlinear functions $\bm{g}_{t}$ in (2.8). It is implicit that
$n$ serves as the index of the sequence, and $p=p^{(n)}$,
$\bm{b}^{*}{}^{(n)}$, $\bm{g}_{t}{}^{(n)}$ all implicitly depend on $n$ and
may change values as $n$ increases, as long as Assumptions 2.2, 2.1, 2.3 and
2.4 are fulfilled for every $n$. The constants
$(T,\zeta,\gamma,\kappa,\sigma^{2})$ do not depend on $n$.

######  Theorem 2.4.

Let Assumptions 2.2, 2.1, 2.3 and 2.4 be fulfilled and assume that both $T$
and ${\rm var}(y_{1})$ are bounded from above by a fixed constant as
$n,p\to+\infty$. Then there exists a random vector $\bm{\zeta}_{j}$ in
${\mathbb{R}}^{T}$ with

(2.22)
$\max_{j\in[p]}W_{2}\Bigl{(}\mathsf{Law}(\bm{\zeta}_{j}),~{}~{}\mathsf{N}\Bigl{(}\mathbf{0}_{T},\bigl{(}{\mathbb{E}}[r_{tt^{\prime}}]\bigr{)}_{(t,t^{\prime})\in[T]\times[T]}\Bigr{)}\Bigr{)}\to
0$

such that for each $t\in[T]$,

(2.23)
$\sum_{j=1}^{p}{\mathbb{E}}\Bigl{[}\Bigl{(}\sqrt{\frac{n}{\|\bm{\Sigma}^{-1/2}\bm{e}_{j}\|^{2}}}\Bigl{(}{\widehat{\bm{b}}}^{t,\rm{debias}}_{j}-\bm{b}^{*}_{j}\Bigr{)}-\bm{\zeta}_{jt}\Bigr{)}^{2}\Bigr{]}\leq{C(\zeta,T,\kappa,\gamma){\rm
var}(y_{1})}.$

The proof of Theorem 2.4 is given in Appendix I. If ${\rm var}(y_{1})$ stays
bounded as $n,p\to+\infty$, the sum over $p$ entries in the left-hand side of
(2.23) stays bounded, and at most a constant number of entries may not
converge to 0. This is made precise in the following corollary.

######  Corollary 2.5.

Let Assumptions 2.2, 2.1, 2.3 and 2.4 be fulfilled and assume that $T$, ${\rm
var}(y_{1})$ are bounded from above by a fixed constant as $n,p\to+\infty$.
Let $a_{p}\to+\infty$ be a slowly increasing sequence (e.g., $a_{p}=\log p$).
There exists a subset $J_{n,p}$ of size at least $p-a_{p}$ such that

(2.24) $\max_{j\in
J_{n,p}}W_{2}\Bigl{(}\sqrt{\frac{n}{\|\bm{\Sigma}^{-1/2}\bm{e}_{j}\|^{2}}}\bigl{(}{\widehat{\bm{b}}}^{t,\rm{debias}}_{j}-\bm{b}^{*}_{j}\bigr{)},~{}~{}\mathsf{N}\bigl{(}0,{\mathbb{E}}[r_{t}]\bigr{)}\Bigr{)}\to
0.$

The proof of Corollary 2.5 is given in Appendix J. Convergence in
2-Wasserstein distance implies convergence in distributions (Villani et al.,
2009, Def. 6.8 and Theorem 6.9). If ${\rm var}(y_{1})/\sigma^{2}$ is bounded
as $n,p\to+\infty$, Theorem 2.2 further shows that $\hat{r}_{t}$ consistently
estimates ${\mathbb{E}}[r_{t}]$, so that by Slutsky’s theorem and (2.24), the
z-score

(2.25)
$\frac{\sqrt{n}({\widehat{\bm{b}}}^{t,\rm{debias}}_{j}-\bm{b}^{*}_{j})}{\|\bm{\Sigma}^{-1/2}\bm{e}_{j}\|\sqrt{\hat{r}_{t}}}$

converge to $\sf N(0,1)$ in distributions uniformly over $j\in J_{n,p}$. For
those overwhelming majority of components $j\in J_{n,p}$, Theorem 2.2 and
Corollary 2.5 thus provide the $1-\alpha$ confidence interval for
$\bm{b}_{j}^{*}:$

$\Big{[}{\widehat{\bm{b}}}_{j}^{t,\rm{debias}}-z_{\alpha/2}\sqrt{\frac{\hat{r}_{t}}{n}(\bm{\Sigma}^{-1})_{jj}},\quad{\widehat{\bm{b}}}^{t,\rm{debias}}+z_{\alpha/2}\sqrt{\frac{\hat{r}_{t}}{n}(\bm{\Sigma}^{-1})_{jj}}\Big{]},$

where $z_{\alpha/2}$ is the standard normal quantile defined by
${\mathbb{P}}(|\mathsf{N}(0,1)|\geq z_{\alpha/2})=\alpha$.

Asymptotic normality “on average” over coordinates, or only over a large
subsets of coordinates, is typical in asymptotic normality results in the
proportional regime, see for instance Bayati and Montanari (2012); Sur and
Candès (2019); Lei et al. (2018); Berthier et al. (2020); Celentano et al.
(2023); Bellec and Zhang (2023). Studying the remaining coordinates ($j\notin
J_{n,p}$ above) remains a challenge, and in some situations the remaining
coordinates exhibit a “variance spike” where the standard deviation estimate
in (2.24) is incorrect and the asymptotic variance of (2.25) is bounded away
from 1 (Bellec and Zhang, 2023, Section 3.7).

One application of Corollary 2.5 is early rejection in hypothesis testing. In
order to perform a hypothesis test of the form

(2.26) $H_{0j}:\bm{b}_{j}^{*}=0\qquad\text{ against }\qquad
H_{1j}:\bm{b}_{j}^{*}\neq 0,$

for the $j$-th component of the true regression vector $\bm{b}^{*}$ in (1.1).
Since the asymptotic normality of the z-score (2.25) is maintained throughout
the iterations, one may perform early tests, without waiting for the
convergence of ${\widehat{\bm{b}}}^{t}$: reject $H_{0j}$ at the $t$-th
iteration if the z-score under $H_{0j}$ is larger than the two-sided quantile
of the normal distribution. We may for instance perform this test at
$t=10,10^{2},10^{3}$ with an appropriate Bonferroni multiple testing
correction. While the estimates $\hat{r}_{t}$ in (2.15) and
$\hat{r}_{tt^{\prime}}$ in (2.18) do not require the knowledge of the
covariance matrix $\bm{\Sigma}$, the construction of the debiased estimate
(2.21) requires the knowledge of $\bm{\Sigma}^{-1}\bm{e}_{j}$ or an estimate
of it. If unlabeled samples $(\bm{x}_{i})_{i=n+1,...,N}$ are available for
instance, $\bm{\Sigma}^{-1}\bm{e}_{j}$ can be estimated by regressing the
$j$-th column onto the others. We emphasize that the knowledge of
$\bm{\Sigma}^{-1}\bm{e}_{j}$ is only used in the construction of debiased
estimate (2.21), and that the estimate $\hat{r}_{t}$ of the generalization
error is readily usable without any knowledge of $\bm{\Sigma}$.

##  3 Concrete examples

In this section, we present the analysis of a few popular algorithms aimed at
solving minimization problems of the form (1.2), including the least-squares
problem, Lasso and MCP. For the five algorithms present in the five
subsections below, we present explicit expressions for iteration functions
$\bm{g}_{t}$ as well as for the corresponding derivative matrices
${\mathcal{D}},{\mathcal{J}}$ and the memory matrix
$\smash{\widehat{\mathbf{A}}}$. For each algorithm, these expressions are then
used to compute the risk estimate $\hat{r}_{t}$ in (2.15) and the z-score in
(2.25). To illustrate our theoretical results, we conduct extensive simulation
studies comparing the proposed risk estimator $\hat{r}_{t}$ with the actual
risk $r_{t}$ at each iteration $t$ and constructing normal quantile-quantile
plots for z-score defined in (2.25) compared to the standard normal
distribution. We start this section with the simulation settings.

Simulation setup. We consider three distinct scenarios based on the
relationship between the sample size and feature dimension, denoted as
$(n,p)$, to generate datasets from the linear model (1.1). These scenarios
are:

  * •

Over-parametrized regime: $(n,p)=(1200,1500)$, in which the number of features
surpasses the number of samples.

  * •

Equal-parametrized regime: $(n,p)=(1200,1200)$, characterized by an equal
number of samples and features.

  * •

Under-parametrized regime: $(n,p)=(1200,500)$, where the number of samples
exceeds the number of features.

For each scenario, the design matrix $\bm{X}$ is generated from a multivariate
normal distribution with zero mean and a covariance matrix
$\bm{\Sigma}=(\bm{\Sigma}_{jk})_{j,k\in[p]}$, where
$\bm{\Sigma}_{jk}=0.5^{|j-k|}$. The noise vector ${\bm{\varepsilon}}$ follows
the standard normal distribution, namely, $\sigma^{2}=1$. The coefficient
vector $\bm{b}^{*}$ is chosen with $p/20$ nonzero entries, set to a constant
value such that the signal-to-noise ratio
$\|\bm{\Sigma}^{1/2}\bm{b}^{*}\|^{2}/\sigma^{2}$ equals 5. For each algorithm
in the following subsections, the initial vector is set to
${\widehat{\bm{b}}}^{1}={\mathbf{0}}$, and the algorithm runs for $T=500$
iterations, except for GD, where $T=3000$ is necessary to achieve convergence.
For each iterate $({\widehat{\bm{b}}}^{t})_{t\in[T]}$, we compute both the
actual risk $r_{t}$ following (2.14) and our proposed risk estimator
$\hat{r}_{t}$ as per (2.15), along with the proposed z-score defined in (2.25)
with $j=1$ (i.e., the first coordinate). Each simulation is replicated 100
times. We present the simulation results for $(n,p)=(1200,1500)$ for each
algorithm in the subsequent subsections, and leave the other $(n,p)$ pair
scenarios to Appendix A since their results are very similar to
$(n,p)=(1200,1500)$. Each simulation we tried confirmed the accuracy of the
estimated risk $\hat{r}_{t}$ for estimating $r_{t}$ as well as the closeness
of the z-scores (2.25) to the standard normal distribution.

###  3.1 Gradient descent

Consider estimating $\bm{b}^{*}$ by minimizing the squared loss function
$\bm{f}(\bm{b})=\frac{1}{2n}\|\bm{y}-\bm{X}\bm{b}\|^{2}.$ The gradient descent
(GD) method finds the minimizer of $\bm{f}(\bm{b})$ by iterations

${\widehat{\bm{b}}}^{t}={\widehat{\bm{b}}}^{t-1}-\eta\nabla\bm{f}({\widehat{\bm{b}}}^{t-1})={\widehat{\bm{b}}}^{t-1}+\frac{\eta}{n}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t-1}),$

with an initialization ${\widehat{\bm{b}}}^{1}\in{\mathbb{R}}^{p}$. Since the
function $\bm{f}$ is an $L$-smooth with
$L=\lVert\bm{X}\rVert_{\rm{op}}^{2}/n$, one can take fixed step size
$\eta=\frac{1}{L}$. Using the definition
$\bm{v}^{t}=\frac{1}{n}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t})$,
the iteration of GD can be written as

(3.1)
${\widehat{\bm{b}}}^{t}=\bm{g}_{t}({\widehat{\bm{b}}}^{t-1},\bm{v}^{t-1})={\widehat{\bm{b}}}^{t-1}+\eta\bm{v}^{t-1}\text{
for }t\geq 2.$

Therefore, the function $\bm{g}_{t}$ is Lipschitz continuous, and the matrices
in (2.3) are given by

(3.2) $\bm{J}_{t,t-1}=\bm{I}_{p},\quad\bm{D}_{t,t-1}=\eta\bm{I}_{p}.$

Hence, for GD with iteration (3.1), the expressions of ${\mathcal{J}}$ and
${\mathcal{D}}$ in (2.7) become

${\mathcal{J}}=\begin{bmatrix}{\mathbf{0}}&{\mathbf{0}}&\cdots&\cdots&{\mathbf{0}}\\\
\bm{I}_{p}&{\mathbf{0}}&\ddots&&\vdots\\\
{\mathbf{0}}&\bm{I}_{p}&{\mathbf{0}}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&{\mathbf{0}}\\\
{\mathbf{0}}&\cdots&{\mathbf{0}}&\bm{I}_{p}&{\mathbf{0}}\\\
\end{bmatrix}=\bm{L}\otimes\bm{I}_{p},\qquad{\mathcal{D}}=\eta{\mathcal{J}}=\eta(\bm{L}\otimes\bm{I}_{p}),$

where $\bm{L}\in{\mathbb{R}}^{T\times T}$ is the strictly lower triangular
matrix $\bm{L}=\sum_{t=2}^{T}\bm{e}_{t}\bm{e}_{t-1}^{\top}$. Now we proceed to
derive the expression for each entry of $\smash{\widehat{\mathbf{A}}}$ in
(2.12) for GD. First, note that ${\mathcal{D}}=\eta{\mathcal{J}}$, so that

$\displaystyle{\mathcal{J}}-{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})=\bm{L}\otimes(\bm{I}_{p}-\eta\bm{X}^{\top}\bm{X}/n).$

Since $\bm{L}$ is strictly lower triangular, the above display is also lower
triangular. With $\bm{\Gamma}=\bm{I}_{p}-\tfrac{\eta\bm{X}^{\top}\bm{X}}{n}$,
we have

$\displaystyle\Bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\Bigr{)}^{-1}=\Bigl{(}\bm{I}_{pT}-(\bm{L}\otimes\bm{\Gamma})\Bigr{)}^{-1}=\sum_{k=0}^{\infty}\bigl{(}\bm{L}\otimes\bm{\Gamma}\bigr{)}^{k}=\sum_{k=0}^{T-1}\bm{L}^{k}\otimes\bm{\Gamma}^{k}$

thanks to $(\bm{L}\otimes\bm{\Gamma})^{k}=\bm{L}^{k}\otimes\bm{\Gamma}^{k}$ by
the mixed-product property and since $\bm{L}^{k}={\mathbf{0}}$ for $k\geq T$.
In this case, by the definition (2.12) of $\smash{\widehat{\mathbf{A}}}$, we
get

$\displaystyle\smash{\widehat{\mathbf{A}}}$
$\displaystyle=\sum_{k=0}^{T-1}\bm{L}^{k+1}\operatorname{Tr}\Bigl{[}\frac{\eta}{n}\bm{X}\bm{\Gamma}^{k}\bm{X}^{\top}\Bigr{]}=\sum_{k=0}^{T-1}\bm{L}^{k+1}\operatorname{Tr}\Bigl{[}(\bm{I}_{p}-\bm{\Gamma})\bm{\Gamma}^{k}\Bigr{]}=\sum_{k=0}^{T-2}\bm{L}^{k+1}\operatorname{Tr}\Bigl{[}(\bm{I}_{p}-\bm{\Gamma})\bm{\Gamma}^{k}\Bigr{]}.$

More visually, since each incremental power of $\bm{L}$ moves the diagonal one
step towards the bottom left, we have

$\displaystyle\smash{\widehat{\mathbf{A}}}$
$\displaystyle=\begin{bmatrix}0&0&0&0&0&\\\
\operatorname{Tr}(\bm{I}_{p}-\bm{\Gamma})&0&0&0&0&\\\
\operatorname{Tr}((\bm{I}_{p}-\bm{\Gamma})\bm{\Gamma})&\operatorname{Tr}(\bm{I}_{p}-\bm{\Gamma})&0&0&0&\\\
\vdots&\ddots&\ddots&\ddots&0&\\\
\operatorname{Tr}((\bm{I}_{p}-\bm{\Gamma})\bm{\Gamma}^{T-2})&\ldots&\operatorname{Tr}((\bm{I}_{p}-\bm{\Gamma})\bm{\Gamma})&\operatorname{Tr}(\bm{I}_{p}-\bm{\Gamma})&0&\end{bmatrix}.$

The above simplifications are specific to GD and will typically not be
possible for examples involving nonlinear transformations such as soft-
thresholding. Computation of $\smash{\widehat{\mathbf{A}}}$ is straightforward
once the eigenvalues of $\bm{\Gamma}$ is computed.

######  Remark 3.1.

A closely related estimator to the GD iterates (3.1) is the minimum $\ell_{2}$
norm least-squares estimator, is defined as

(3.3) $\displaystyle{\widetilde{\bm{b}}}=\mathop{\rm
arg\,min}\big{\\{}\|\bm{b}\|:\bm{b}\text{ minimizes
}\|\bm{y}-\bm{X}\bm{b}\|^{2}\big{\\}}.$

For $n>p$,
${\widetilde{\bm{b}}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}$, and for
$n\leq p$,
${\widetilde{\bm{b}}}=\bm{X}^{\top}(\bm{X}\bm{X}^{\top})^{-1}\bm{y}$. It is
equivalent to write
${\widetilde{\bm{b}}}=(\bm{X}^{\top}\bm{X})^{\dagger}\bm{X}^{\top}\bm{y}$
where $\dagger$ denotes the Moore-Penrose inverse. The GD iterates
${\widehat{\bm{b}}}^{t}$ in (3.1) converge to the min-norm least-squares
solution ${\widetilde{\bm{b}}}$ as $t\to\infty$ (Hastie et al., 2022,
Proposition 1). We observe this phenomenon in the simulations below where the
risk of ${\widetilde{\bm{b}}}$ is represented by an horizontal green curve.

Simulation results. Define
$r_{\infty}:=\lVert\bm{\Sigma}^{1/2}({\widehat{\bm{b}}}^{\infty}-\bm{b}^{*})\rVert+\sigma^{2}$,
where ${\widehat{\bm{b}}}^{\infty}:=\lim_{t\to\infty}{\widehat{\bm{b}}}^{t}$
is the min-norm least-squares estimator (3.3), and has closed-form expression
${\widehat{\bm{b}}}^{\infty}={\widetilde{\bm{b}}}=(\bm{X}^{\top}\bm{X})^{\dagger}\bm{X}^{\top}\bm{y}$.
At each iteration $t\in[T]$, we calculate and present the average estimated
risk $\hat{r}_{t}$, the actual risk $r_{t}$, and the limiting risk
$r_{\infty}$ over 100 repetitions, including 2-standard error bars, in Figure
1(a). We also provide the quantile-quantile (Q-Q) plot of the z-score (2.25)
at different iterations in Figure 1(b).

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score at different iterations.

Figure 1: Risk curves and qq-plots of z-score of GD for $(n,p)=(1200,1500)$.

Figure 1(a) shows a strong alignment between our risk estimate $\hat{r}_{t}$
and the actual risk $r_{t}$. Notably, both the risk and the estimated risk
reach their minimum at iteration 8, suggesting that ${\widehat{\bm{b}}}^{8}$
yields the lowest out-of-sample prediction risk. This observation indicates
that it is beneficial to terminate the algorithm as early as iteration 8,
rather than continuing with additional iterations, because the risk will blow
up quickly after iteration 8. As $t$ increases, we observe from Figure 1(a)
that $\hat{r}_{t}$ converges to $r_{\infty}$, indicating the effectiveness of
our risk estimator for large $t$. Furthermore, Figure 1(b) reveals that the
quantiles closely match the 45-degree line (shown in red). This alignment
supports the conclusion that the z-score (2.25) closely approximates a
standard normal distribution.

###  3.2 Nesterov’s accelerated gradient descent (AGD)

The Nesterov’s accelerated gradient method (Nesterov, 1983) is a remarkable
extension of the gradient descent by utilizing the momentum from previous
iterates. It is well known that AGD enjoys quadratic convergence rates, which
is faster than the linear convergence rate of the gradient descent. To
describe the iteration of AGD, define the sequence of scalars

(3.4) $a_{0}=0,\quad a_{t}=\tfrac{1+\sqrt{1+4a_{t-1}^{2}}}{2},\quad
w_{t}=\tfrac{1-a_{t}}{a_{t+1}}\quad\text{ for all }t\geq 1.$

From some initialization ${\widehat{\bm{b}}}^{1}\in{\mathbb{R}}^{p}$, AGD
iterates are defined as the weighted sum

(3.5)
$\displaystyle{\widehat{\bm{b}}}^{t}=(1-w_{t-1})({\widehat{\bm{b}}}^{t-1}+\eta\bm{v}^{t-1})+w_{t-1}({\widehat{\bm{b}}}^{t-2}+\eta\bm{v}^{t-2}),\qquad\text{
for all }t\geq 2.$

Hence, the Jacobian matrices defined in (2.6) are given by

$\displaystyle\bm{J}_{t,t-1}$ $\displaystyle=(1-w_{t-1})\bm{I}_{p},\quad$
$\displaystyle\bm{D}_{t,t-1}$ $\displaystyle=\eta(1-w_{t-1})\bm{I}_{p},$
$\displaystyle\bm{J}_{t,t-2}$ $\displaystyle=w_{t-1}\bm{I}_{p},\quad$
$\displaystyle\bm{D}_{t,t-2}$ $\displaystyle=\eta w_{t-1}\bm{I}_{p}.$

It follows that the expressions of ${\mathcal{J}}$ and ${\mathcal{D}}$ in
(2.7) become

${\mathcal{J}}=\begin{bmatrix}{\mathbf{0}}&{\mathbf{0}}&\cdots&\cdots&\cdots&{\mathbf{0}}\\\
(1-w_{1})\bm{I}_{p}&{\mathbf{0}}&\ddots&\cdots&\cdots&\vdots\\\
w_{2}\bm{I}_{p}&(1-w_{2})\bm{I}_{p}&{\mathbf{0}}&\ddots&\cdots&\vdots\\\
{\mathbf{0}}&w_{3}\bm{I}_{p}&(1-w_{3})\bm{I}_{p}&{\mathbf{0}}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&{\mathbf{0}}\\\
{\mathbf{0}}&\cdots&{\mathbf{0}}&w_{T-1}\bm{I}_{p}&(1-w_{T-1})\bm{I}_{p}&{\mathbf{0}}\end{bmatrix},\quad{\mathcal{D}}=\eta{\mathcal{J}}.$

######  Remark 3.2.

By Proposition 1 in Hastie et al. (2022), the AGD iterate
${\widehat{\bm{b}}}^{t}$ in (3.5) also converges to the min-norm least-squares
estimator ${\widetilde{\bm{b}}}$ as $t\to\infty$.

Simulation results. Similar to the plots for gradient descent, we provide
results for AGD applied to the least-squares problem. The risk curves are
shown in Figure 2(a) and Q-Q plots of the z-score (2.25) in Figure 2(b). The
simulation results for AGD are similar to those for GD. Figure 2(a) clearly
shows that the risk estimate $\hat{r}_{t}$ closely matches the actual risk
$r_{t}$. In addition, the risk curve suggests stopping the algorithm early at
iteration 8 as the estimated risk increases quickly after iteration 8. Both
$r_{t}$ and $\hat{r}_{t}$ converge to $r_{\infty}$ as $t$ increases. Figure
2(b) again confirms that the z-scores (2.25) are close to a standard normal
distribution.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score at different iterations.

Figure 2: Risk curves and qq-plots of z-score of AGD for $(n,p)=(1200,1500)$.

###  3.3 Iterative Shrinkage-Thresholding Algorithm (ISTA)

For regression problems with high-dimensional features, penalized regression
are useful to achieve sparse solutions. Consider the Lasso regression

(3.6) ${\widehat{\bm{b}}}\in\mathop{\rm
arg\,min}_{\bm{b}\in{\mathbb{R}}^{p}}\frac{1}{2n}\|\bm{y}-\bm{X}\bm{b}\|^{2}+\lambda\|\bm{b}\|_{1}.$

This objective function of the above optimization has two parts, one is the
squared loss, the other is an $\ell_{1}$ penalty. ISTA (Daubechies et al.,
2004) is a simple algorithm to solve (3.6) by imposing a soft-thresholding
nonlinearity at each iteration. Concretely, let
$\operatorname{soft}_{\theta}:{\mathbb{R}}^{p}\to{\mathbb{R}}^{p}$ denote the
elementwise soft-thresholding operator, i.e.
$\operatorname{soft}_{\theta}(\bm{b})_{j}=(|\bm{b}_{j}-\theta|)_{+}\text{sgn}(\bm{b}_{j}).$
ISTA can be viewed as the proximal gradient descent in Section 2.1.1, and its
iteration function $\bm{g}_{t}$ is given by

(3.7)
${\widehat{\bm{b}}}^{t}=\bm{g}_{t}({\widehat{\bm{b}}}^{t-1},\bm{v}^{t-1})=\operatorname{soft}_{\lambda/L}({\widehat{\bm{b}}}^{t-1}+L^{-1}\bm{v}^{t-1}).$

The function $\bm{g}_{t}$ is Lipschitz continuous since the soft-thresholding
is 1-Lipschitz. Let
$S^{t}=\big{\\{}j\in[p]:|{\widehat{\bm{b}}}^{t}_{j}+L^{-1}\bm{v}^{t}_{j}|>\lambda/L\big{\\}}$.
For ISTA, the expressions of $\bm{J}_{t,t-1}$ and $\bm{D}_{t,t-1}$ in (2.3)
are the diagonal matrices

$\displaystyle\bm{J}_{t,t-1}=\sum_{j\in
S^{t-1}}\bm{e}_{j}\bm{e}_{j}^{\top}=\text{diag}\Bigl{(}\bigl{(}\mathbb{I}\\{j\in
S^{t-1}\\}\bigr{)}_{j\in[p]}\Bigr{)},\qquad\bm{D}_{t,t-1}=L^{-1}\bm{J}_{t,t-1}.$

Substituting the above into (2.4) gives the expressions of ${\mathcal{J}}$ and
${\mathcal{D}}$ for ISTA.

Simulation results. We apply ISTA to solve the Lasso regression (3.6) with two
regularization parameters $\lambda\in\\{0.01,0.1\\}$. The ISTA iterates
converge to the Lasso estimator (3.6) (see, e.g., Beck and Teboulle (2009)).
We compute the Lasso using the Python module sklearn.linear_model.Lasso and
use this Lasso estimator to evaluate the limiting risk $r_{\infty}$. Figure
3(a) showcases the risk estimator $\hat{r}_{t}$ and actual risk $r_{t}$ for
each iteration $t$ using ISTA. Furthermore, Figure 3(b) displays the Q-Q plots
of the z-score (2.25) for ISTA. Again, the estimated risk curve closely aligns
with the actual risk curve for both values of $\lambda$, and the corresponding
z-scores closely approximate the standard normal distribution.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 3: Risk curves and qq-plots of z-score of ISTA for $(n,p)=(1200,1500)$.

###  3.4 Fast Iterative Shrinkage-Thresholding Algorithm (FISTA)

Similar to the extension from GD to AGD, FISTA (Beck and Teboulle, 2009) is an
accelerated version of ISTA, incorporating momentum with the weights (3.4).
One advantage of FISTA is that it enjoys faster convergence rate than ISTA
(Beck and Teboulle, 2009). Using the same definitions of $a_{t},w_{t}$ in
(3.4), FISTA iterates the following steps with some initialization
${\widehat{\bm{b}}}^{1}\in{\mathbb{R}}^{p}$:

$\displaystyle{\widehat{\bm{b}}}^{t}$
$\displaystyle=\bm{g}_{t}({\widehat{\bm{b}}}^{t-1},{\widehat{\bm{b}}}^{t-2},\bm{v}^{t-1},\bm{v}^{t-2})$
$\displaystyle=(1-w_{t-1})\operatorname{soft}_{\lambda/L}({\widehat{\bm{b}}}^{t-1}+L^{-1}\bm{v}^{t-1})+w_{t-1}\operatorname{soft}_{\lambda/L}({\widehat{\bm{b}}}^{t-2}+L^{-1}\bm{v}^{t-2})\text{
for }t\geq 2.$

Using
$S^{t}=\big{\\{}j\in[p]:|{\widehat{\bm{b}}}^{t}_{j}+L^{-1}\bm{v}^{t}_{j}|>\lambda/L\big{\\}}$,
the matrices in (2.6) are the diagonal matrices

$\displaystyle\bm{J}_{t,t-1}$ $\displaystyle=(1-w_{t-1})\sum_{j\in
S^{t-1}}\bm{e}_{j}\bm{e}_{j}^{\top},\quad$ $\displaystyle\bm{D}_{t,t-1}$
$\displaystyle=\frac{1}{L}\bm{J}_{t,t-1},$ $\displaystyle\bm{J}_{t,t-2}$
$\displaystyle=w_{t-1}\sum_{j\in S^{t-2}}\bm{e}_{j}\bm{e}_{j}^{\top},\quad$
$\displaystyle\bm{D}_{t,t-2}$ $\displaystyle=\frac{1}{L}\bm{J}_{t,t-2}.$

We obtain the expressions of ${\mathcal{D}}$ and ${\mathcal{J}}$ for FISTA by
substituting the above into (2.7). Similarly to the simulation results for
ISTA, we present the risk curves and Q-Q plots of the z-score (2.25) for FISTA
in Figure 4.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 4: Risk curves and qq-plots of z-score of FISTA for
$(n,p)=(1200,1500)$.

###  3.5 Local quadratic approximation in non-convex penalized regression

In this subsection, we consider the folded-concave penalized least-squares,

(3.8)
$\min_{\bm{b}\in{\mathbb{R}}^{p}}\frac{1}{2n}\|\bm{y}-\bm{X}\bm{b}\|^{2}+\sum_{j=1}^{p}\rho(|\bm{b}_{j}|),$

where $\rho:{\mathbb{R}}_{+}\to{\mathbb{R}}$ is a concave penalty function.
This encompasses SCAD (Fan and Li, 2001) and MCP (Zhang, 2010). For
simulations, we set $\rho(\cdot)$ to be the MCP penalty with two positive
tuning parameters $(\lambda,\tau)$, defined as

(3.9) $\displaystyle\rho(x;\lambda,\tau)=\begin{cases}\lambda
x-\frac{1}{2\tau}x^{2}&\text{if~{}}x\leq\tau\lambda\\\
\frac{1}{2}\tau\lambda^{2}&\text{if~{}}x>\tau\lambda\end{cases}.$

For simplicity, we will omit the parameters $(\lambda,\tau)$ in the notation,
using $\rho(x)$ instead of $\rho(x;\lambda,\tau)$. Consequently, the
derivative of $\rho$ becomes
$\rho^{\prime}(x)=(\lambda-x/\tau){\mathbb{I}}(x\leq\tau\lambda)$. As
$\tau\to\infty$, $\rho(x)$ becomes $\lambda x$, and the optimization problem
(3.8) then coincides with the Lasso (3.6).

In order to solve the non-convex penalized regression (3.8), we consider the
local isotropic quadratic approximation (LQA) to the least-squares loss
$f(\bm{b}):=\frac{1}{2n}\|\bm{y}-\bm{X}\bm{b}\|^{2}$ at a vector
${\widehat{\bm{b}}}^{t-1}$, namely

(3.10) $f(\bm{b})\approx
f({\widehat{\bm{b}}}^{t-1})+(\bm{b}-{\widehat{\bm{b}}}^{t-1})^{\top}\nabla
f({\widehat{\bm{b}}}^{t-1})+(L/2)\|\bm{b}-{\widehat{\bm{b}}}^{t-1}\|^{2},$

where $\nabla f({\widehat{\bm{b}}}^{t-1})$ is the gradient of $f$ evaluated at
${\widehat{\bm{b}}}^{t-1}$, and $L=n^{-1}\lVert\bm{X}\rVert_{\rm{op}}^{2}$ as
in previous sections. Applying the above quadratic approximation (3.10) to the
least-squares in (3.8) and ignoring the constant term yields the optimization
problem

(3.11) ${\widehat{\bm{b}}}^{t}=\mathop{\rm
arg\,min}_{\bm{b}\in{\mathbb{R}}^{p}}\frac{1}{2}\|\bm{b}-({\widehat{\bm{b}}}^{t-1}+(nL)^{-1}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t-1}))\|^{2}+\frac{1}{L}\sum_{j=1}^{p}\rho(|\bm{b}_{j}|).$

Let ${\widehat{\bm{b}}}_{j}^{t}$ and $\bm{v}_{j}^{t}$ be the $j$-th entry of
${\widehat{\bm{b}}}^{t}$ and $\bm{v}^{t}$, respectively. For $\tau\geq 1/L$,
the optimization problem (3.11) admits the closed-form solution

(3.12) $\displaystyle{\widehat{\bm{b}}}^{t}_{j}$
$\displaystyle=\begin{cases}\operatorname{soft}_{\lambda/L}\Bigl{(}{\widehat{\bm{b}}}^{t-1}_{j}+\frac{1}{L}\bm{v}^{t-1}_{j}\Bigr{)}(1-\tfrac{1}{\tau
L})^{-1}&\text{if
}|{\widehat{\bm{b}}}^{t-1}_{j}+\tfrac{1}{L}\bm{v}^{t-1}_{j}|\leq\tau\lambda,\\\
{\widehat{\bm{b}}}^{t-1}_{j}+\frac{1}{L}\bm{v}^{t-1}_{j}&\text{otherwise}.\end{cases}.$

In this case, $\bm{J}_{t,t-1}$ and $\bm{D}_{t,t-1}$ in (2.3) are the diagonal
matrices

$\displaystyle\bm{J}_{t,t-1}$
$\displaystyle=\sum_{j=1}^{p}\bm{e}_{j}\bm{e}_{j}^{\top}\Bigl{[}(1-\tfrac{1}{\tau
L})^{-1}{\mathbb{I}}\Bigl{(}|{\widehat{\bm{b}}}^{t-1}_{j}+\tfrac{1}{L}\bm{v}^{t-1}_{j}|\in[\lambda/L,\tau\lambda]\Bigr{)}+{\mathbb{I}}\Bigl{(}|{\widehat{\bm{b}}}^{t-1}_{j}+\tfrac{1}{L}\bm{v}^{t-1}_{j}|>\tau\lambda\Bigr{)}\Bigr{]}$

and $\bm{D}_{t,t-1}=\tfrac{1}{L}\bm{J}_{t,t-1}$.

######  Remark 3.3.

If $\tau=\infty$, the MCP function reduces to
$\sum_{j=1}^{p}\rho(|\bm{b}_{j}|)=\lambda\lVert\bm{b}\rVert_{1}$, so MCP is
the same as the Lasso. In this case, the LQA iterations (3.12) become
${\widehat{\bm{b}}}^{t}=\operatorname{soft}_{\lambda/L}({\widehat{\bm{b}}}^{t-1}+\frac{1}{L}\bm{v}^{t-1})$,
which is the same as ISTA iterations (3.7) in Section 3.3.

Simulation results. For the MCP penalty function (3.9), we consider
$\lambda\in\\{0.1,0.2\\}$ and $\tau=3$. We display the curves for the risk
estimator $\hat{r}_{t}$ and the actual risk $r_{t}$ for each iteration $t$ in
Figure 5(a). Additionally, the Q-Q plots of the z-score (2.25) are given in
Figure 5(b). Figure 5(a) shows that the estimated risk accurately estimates
the true risk. The z-scores (2.25) closely approximate the standard normal
distribution. Similar to observations in the aforementioned algorithms, LQA
reaches its lowest risk level at around iteration 10. This suggests that early
stopping could be beneficial for LQA for certain tuning parameters to improve
generalization performance. When comparing the lowest points on the risk
curves for LQA with other algorithms (GD, AGD, ISTA, and FISTA), LQA achieves
the lower risk among the tested algorithms at the given tuning parameters.
Figure 5(b) provides empirical support for the established asymptotic
normality of the debiased LQA iterates in Theorem 2.4.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 5: Risk curves and qq-plots of z-score of LQA for $(n,p)=(1200,1500)$.

##  4 Efficient computation of the memory matrix

Recall that the memory matrix $\smash{\widehat{\mathbf{A}}}$ is crucial in the
formulae of proposed risk estimator $\hat{r}_{t}$ in (2.15) and debiased
estimator ${\widehat{\bm{b}}}_{j}^{t,\rm debias}$ in (2.21). We have provided
specific expressions of $\smash{\widehat{\mathbf{A}}}$ for various algorithms
in Section 3. In this section, we provide an efficient way to compute the
memory matrix $\smash{\widehat{\mathbf{A}}}$ in (2.12) for algorithms with the
general iterations (2.8).

###  4.1 Iteratively computing rows of $\smash{\widehat{\mathbf{A}}}$

Recall the definition of $\smash{\widehat{\mathbf{A}}}$ in (2.12) is

$\displaystyle\smash{\widehat{\mathbf{A}}}=\frac{1}{n}\sum_{i=1}^{n}\bigl{(}\bm{I}_{T}\otimes(\bm{e}_{i}^{\top}\bm{X})\bigr{)}\Bigl{(}\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}\Bigr{)}^{-1}{\mathcal{D}}\bigl{(}\bm{I}_{T}\otimes(\bm{X}^{\top}\bm{e}_{i})\bigr{)}.$

The first apparent computational hurdle lies in inverting the large matrix
$\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}$
of size ${pT\times pT}$. We now provide an efficient way to compute the memory
matrix $\smash{\widehat{\mathbf{A}}}$ without explicitly inverting this large
matrix. Recall that ${\mathcal{D}}$ and ${\mathcal{J}}$ are $T\times T$ block
lower triangular matrices given in (2.7) or (2.10), where each block is of
size $p\times p$. It follows that
${\mathcal{J}}-{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})$
is also a $T\times T$ block lower triangular matrix with zero diagonal blocks.
Let $\bm{M}_{t}\in{\mathbb{R}}^{p\times(pT)}$ be the $t$-th block row of
$\mathcal{D}$. Consider the linear system with unknowns
$\bm{R}_{1}\in{\mathbb{R}}^{p\times
pT},\dots,\bm{R}_{T}\in{\mathbb{R}}^{p\times pT}$ given by

$\Bigl{(}\bm{I}_{pT}-{\mathcal{J}}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})\Bigr{)}\begin{bmatrix}\bm{R}_{1}\\\
\bm{R}_{2}\\\ \vdots\\\ \bm{R}_{T-1}\\\
\bm{R}_{T}\end{bmatrix}=\begin{bmatrix}\bm{M}_{1}\\\ \bm{M}_{2}\\\ \vdots\\\
\bm{M}_{T-1}\\\ \bm{M}_{T}\end{bmatrix}={\mathcal{D}}$

which can be rewritten with
$\bm{P}_{t,s}=\bm{J}_{t,s}-\bm{D}_{t,s}\bm{X}^{\top}\bm{X}/n$ as

(4.1) $\begin{bmatrix}\bm{I}_{p}&{\mathbf{0}}_{p\times
p}&\dots&\dots&{\mathbf{0}}_{p\times p}\\\
-\bm{P}_{2,1}&\bm{I}_{p}&{\mathbf{0}}_{p\times p}&&\vdots\\\
-\bm{P}_{3,1}&-\bm{P}_{3,2}&\bm{I}_{p}&\ddots&\vdots\\\ \
\vdots&\ddots&\ddots&\ddots&{\mathbf{0}}_{p\times p}\\\
-\bm{P}_{T,1}&&\dots&-\bm{P}_{T,T-1}&\bm{I}_{p}\end{bmatrix}\begin{bmatrix}\bm{R}_{1}\\\
\bm{R}_{2}\\\ \vdots\\\ \bm{R}_{T-1}\\\
\bm{R}_{T}\end{bmatrix}=\begin{bmatrix}\bm{M}_{1}\\\ \bm{M}_{2}\\\ \vdots\\\
\bm{M}_{T-1}\\\ \bm{M}_{T}\end{bmatrix}={\mathcal{D}}.$

The left-most matrix being lower-triangular with identity diagonal blocks,
forward substitution provides the unique solution: starting with
$\bm{R}_{1}=\bm{M}_{1}$ we obtain directly
$\bm{R}_{2}=\bm{M}_{2}+\bm{P}_{2,1}\bm{R}_{1}$, and more generally for all
$t\geq 1$,

(4.2)
$\bm{R}_{t}=\bm{M}_{t}+\sum_{s=1}^{t-1}\bm{P}_{t,s}\bm{R}_{s}=\sum_{s=1}^{t-1}\bm{e}_{s}^{\top}\otimes\bm{D}_{t,s}+\sum_{s=1}^{t-1}\bm{P}_{t,s}\bm{R}_{s},$

where the second equality follows from the observation that ${\mathcal{D}}$ is
block lower triangular with $t$-th row block
$\bm{M}_{t}=\sum_{s=1}^{t-1}\bm{e}_{s}^{\top}\otimes\bm{D}_{t,s}$. Given
$(\bm{R}_{s})_{s\leq t-1}$, compute $\bm{R}_{t}$ according to (4.2) and set

(4.3)
$\smash{\widehat{\mathbf{A}}}_{t,t^{\prime}}=\tfrac{1}{n}\operatorname{Tr}\bigl{[}\bm{X}\bm{R}_{t}(\bm{e}_{t^{\prime}}\otimes\bm{X}^{\top})\bigr{]}\quad\text{
for
}t^{\prime}<t,\quad\smash{\widehat{\mathbf{A}}}_{t,t^{\prime}}=0\quad\text{
for }t^{\prime}\geq t.$

If only one or two blocks $\bm{P}_{t,s}$ per row are nonzero as in (2.7),
which is satisfied for all examples of Section 3, the recursion (4.2)
simplifies in this case to

(4.4)
$\bm{R}_{t}=\sum_{s=t-2}^{t-1}\bm{e}_{s}^{\top}\otimes\bm{D}_{t,s}+\sum_{s=t-2}^{t-1}\bigl{(}\bm{J}_{t,s}-\bm{D}_{t,s}\tfrac{\bm{X}^{\top}\bm{X}}{n}\bigr{)}\bm{R}_{s},$

followed by (4.3). Notably, we may compute each $\bm{R}_{t}$ recursively while
only keeping in memory $\bm{R}_{t-1}$ and $\bm{R}_{t-2}$, both of size
$p\times(pT)$, at each step.

###  4.2 Hutchinson’s trace approximation

The above computation of $\smash{\widehat{\mathbf{A}}}$ using (4.3) can still
be prohibitive, as it requires storing in memory matrices
$\bm{R}_{t},\bm{R}_{t-1}$ of size at most $p\times(pT)$ at each step, and
perform matrix-matrix products with dimensions $p\times p$ and $p\times(pT)$
in (4.4). We now describe an efficient way to approximate the entries of
$\smash{\widehat{\mathbf{A}}}$ while avoiding to store intermediate matrices
of size $p\times(pT)$. Let $m\geq 1$ be a small constant integer; in
simulations we have noticed that $m=1,2$ or $3$ already gives accurate
estimates. We propose to approximate the trace (4.3) for the $(t,s)$-th entry
of $\smash{\widehat{\mathbf{A}}}$ using Hutchinson’s trace approximation
(Hutchinson, 1990)

$\operatorname{Tr}(\bm{M})\approx\operatorname{Tr}(\bm{W}^{\top}\bm{M}\bm{W})$

for any matrix $\bm{M}\in{\mathbb{R}}^{n\times n}$, where
$\bm{W}\in{\mathbb{R}}^{n\times m}$ is a random matrix with i.i.d. entries
uniformly distributed in $\\{\frac{1}{\sqrt{m}},\frac{-1}{\sqrt{m}}\\}$. The
Hanson-Wright inequality ensures that the above approximation holds with high-
probability when $\|\bm{M}\|_{\rm F}$ is negligible compared to
$\operatorname{Tr}(\bm{M})$. Using this approximation, the computation of the
trace of a matrix of size $n\times n$ is reduced to that of a matrix of size
$m\times m$.

In our case, we need to approximate the trace in (4.3) using Hutchinson’s
approximation. To this end, let $\bm{W}\in{\mathbb{R}}^{n\times m}$ be a
random matrix with i.i.d. entries uniformly distributed in
$\\{\frac{1}{\sqrt{m}},\frac{-1}{\sqrt{m}}\\}$ independently of everything
else. The matrix $\bm{W}$ is only sampled once, and fixed throughout the
following computation. Instead of computing recursively $\bm{R}_{t}$ in (4.2),
we compute
$\bar{\bm{R}}_{t}=\bm{R}_{t}(\bm{I}_{T}\otimes(\bm{X}^{\top}\bm{W}))$
iteratively from (4.2) with the recursion

(4.5)
$\bar{\bm{R}}_{t}=\sum_{s=1}^{t-1}\bm{e}_{s}^{\top}\otimes(\bm{D}_{t,s}\bm{X}^{\top}\bm{W})+\sum_{s=1}^{t-1}\Bigl{(}\bm{J}_{t,s}-\bm{D}_{t,s}\frac{\bm{X}^{\top}\bm{X}}{n}\Bigr{)}\bar{\bm{R}}_{s}.$

Once $\bar{\bm{R}}_{t}$ is available, we compute the Hutchinson’s
approximation for each entry of the $t$-th row of
$\smash{\widehat{\mathbf{A}}}$ in (4.3) by

$\displaystyle\smash{\widehat{\mathbf{A}}}_{t,t^{\prime}}^{H}$
$\displaystyle=\frac{1}{n}\operatorname{Tr}\Bigl{[}\bm{W}^{\top}\bm{X}\bar{\bm{R}}_{t}(\bm{e}_{t^{\prime}}\otimes\bm{I}_{m})\Bigr{]}$

for $t^{\prime}<t$ and $0$ for $t^{\prime}\geq t$. For AGD and FISTA or any
iterative algorithm with ${\mathcal{D}},{\mathcal{J}}$ given by (2.7),
$\bm{D}_{t,s}$ and $\bm{J}_{t,s}$ are 0 except for $s=t-1$ and $s=t-2$. In
this case the sums in the recursion for $\bm{R}_{t}$ in (4.2) and
$\bar{\bm{R}}_{t}$ in (4.5) are reduced to $\sum_{s=t-2}^{t-1}$ with only two
terms. For ISTA and GD, $\bm{D}_{t,s}$ and $\bm{J}_{t,s}$ are 0 except for
$s=t-1$ so the sum is reduced to only one term at $s=t-1$. In this case, the
recursion (4.5) only uses $mT$ matrix-vector products with matrix dimensions
smaller than $\max\\{n,p,T\\}$. In particular, if $mT\ll\min\\{n,p\\}$, it is
never needed to perform a matrix-matrix multiplication with two matrices with
both dimensions of order $n$ or $p$. In terms of memory footprint, in the case
(2.7), only the last two $\bar{\bm{R}}_{t-1}$ and $\bar{\bm{R}}_{t-2}$ are
needed to compute $\bar{\bm{R}}_{t}$. This is the same cost as storing a
matrix of size $p\times(2mT)$, which is negligible compared to storing
$\bm{X}\in{\mathbb{R}}^{n\times p}$ as long as $mT\ll n$.

Finally, remark that $\bar{\bm{R}}_{t}$ is 0 except in its first $t-1$ column
blocks, so that only these first $t-1$ columns blocks need to be stored and
computed in the recursion (4.5). More precisely, let
$\bm{V}_{t,t^{\prime}}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\bm{R}_{t}(\bm{e}_{t^{\prime}}\otimes\bm{X}^{\top}\bm{W})\in{\mathbb{R}}^{p\times
m}$ for $t^{\prime}\leq t-1$ to be $t^{\prime}$ column block of
$\bar{\bm{R}}_{t}$. Then Hutchinson’s approximation equals

(4.6)
$\displaystyle\smash{\widehat{\mathbf{A}}}_{t,t^{\prime}}^{H}=\tfrac{1}{n}\operatorname{Tr}[\bm{W}^{\top}\bm{X}\bm{V}_{t,t^{\prime}}].$

At step $t$, in order to compute $\bm{V}_{t}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}[\bm{V}_{t,1},..,\bm{V}_{t,t-1}]\in{\mathbb{R}}^{p\times m(t-1)}$
in the case (2.7) satisfied by all examples of Section 3, we have the
recursion formula

$\displaystyle\bm{V}_{t}=\bm{P}_{t,t-2}\bigl{[}\bm{V}_{t-2},~{}{\mathbf{0}}_{p},{\mathbf{0}}_{p}\bigr{]}+\bm{P}_{t,t-1}\bigl{[}\bm{V}_{t-1},~{}{\mathbf{0}}_{p}\bigr{]}+\bigl{[}{\mathbf{0}}_{p\times
m(t-3)},~{}\bm{D}_{t,t-2}\bm{X}^{\top}\bm{W},~{}\bm{D}_{t,t-1}\bm{X}^{\top}\bm{W}\bigr{]},$

where $\bm{P}_{t,s}=\bm{J}_{t,s}-\bm{D}_{t,s}\bm{X}^{\top}\bm{X}/n$ as above.
In the case of ISTA and FISTA where the nonlinear functions use the soft-
thresholding operator, these approximations let us compute the iterations
(1.6) and the memory matrix (2.12) with problem dimensions
$n=25,000;p=40,000;T=30$ on a laptop with 32GB of RAM within four minutes.

##  5 Discussion

This paper introduces a novel procedure to estimate the out-of-sample
prediction error for iterates of gradient descent type algorithms, in the
context of high-dimensional regression. As illustrated in Section 3, this
risk-estimation procedure is applicable to a wide range of algorithms commonly
used in optimization, including gradient descent, Nesterov’s accelerated
gradient descent, iterative shrinkage-thresholding algorithm, fast iterative
shrinkage-thresholding algorithm (Beck and Teboulle, 2009), and local
quadratic approximation. The proposed procedure do not require the knowledge
of the design covariance $\bm{\Sigma}$ or of the noise level $\sigma^{2}$. The
estimates allow the statistician to leverage the benefits of early stopping,
by selecting an early iteration $\hat{t}\in[T]$ that minimizes the
generalization error among the first $T$ iterations up to negligible error
(Corollary 2.3). We have further established the asymptotic normality of the
entries of the iterates after a debiasing correction, which can be used to
construct confidence intervals for the $j$-th of the ground-truth $\bm{b}^{*}$
when $\bm{\Sigma}^{-1}\bm{e}_{j}$ is known or can be estimated (Theorem 2.4).
Extensive numerical simulations in Section 3 demonstrate that the proposed
estimate is accurate, and that the z-scores defining the confidence intervals
are approximately standard normal.

Here we highlight several directions for further exploration. A first avenue
for future research is to improve the dependence of the bounds on the number
$T$ of iterations. Currently, upper bounds in Theorem 2.2 and other main
results involve constants that worst than exponential in $T$, while we observe
in simulations that the risk estimate is still accurate over the whole
trajectory for $T\gg\log n$. Another direction of interest is to generalize
the estimates of the present paper beyond the square loss in (1.2), for
instance with the Huber or least-absolute deviation loss in robust regression,
or the logistic loss in classification problems. Another generalization
concerns randomized versions of GD such as stochastic gradient descent.
Finally, our proofs rely crucially on the Gaussianity of the design, and it
would be of interest to extend the validity of the estimates proposed here to
different design distributions. Recent progress has been made in this
direction (Montanari and Saeed, 2022; Hu and Lu, 2022; Han and Shen, 2023;
Pesce et al., 2023; Dudeja et al., 2023) regarding the universality of the
training and generalization error of minimizers such as (1.2). An extra
challenge presented here is that the validity of $\hat{r}_{t}$ for estimating
$r_{t}$ requires not only universality of the training and generalization
error, but also universality the weights $\hat{w}_{t,s}$.

Supplement A This supplement contains addtional numerical results and proofs.
Supplement B This supplement contains the code and instruction to produce the
numerical results.

## References

  * Ali et al. (2019) Alnur Ali, J Zico Kolter, and Ryan J Tibshirani.  A continuous-time view of early stopping for least squares regression.  In _The 22nd international conference on artificial intelligence and statistics_ , pages 1370–1378. PMLR, 2019.
  * Ali et al. (2020) Alnur Ali, Edgar Dobriban, and Ryan Tibshirani.  The implicit regularization of stochastic gradient flow for least squares.  In _International conference on machine learning_ , pages 233–244. PMLR, 2020.
  * Ambrosio and Dal Maso (1990) Luigi Ambrosio and Gianni Dal Maso.  A general chain rule for distributional derivatives.  _Proceedings of the American Mathematical Society_ , 108(3):691–702, 1990.
  * Bandeira et al. (2013) Afonso S. Bandeira, Edgar Dobriban, Dustin G. Mixon, and William F. Sawin.  Certifying the restricted isometry property is hard.  _IEEE Trans. Inform. Theory_ , 59(6):3448–3450, 2013.  ISSN 0018-9448.
  * Bayati and Montanari (2011) Mohsen Bayati and Andrea Montanari.  The dynamics of message passing on dense graphs, with applications to compressed sensing.  _IEEE Trans. Inform. Theory_ , 57(2):764–785, 2011.  ISSN 0018-9448.
  * Bayati and Montanari (2012) Mohsen Bayati and Andrea Montanari.  The LASSO risk for Gaussian matrices.  _IEEE Trans. Inform. Theory_ , 58(4):1997–2017, 2012.  ISSN 0018-9448.
  * Beck and Teboulle (2009) Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM J. Imaging Sci._ , 2(1):183–202, 2009.
  * Bellec and Shen (2022) Pierre C Bellec and Yiwei Shen.  Derivatives and residual distribution of regularized m-estimators with application to adaptive tuning.  In _Conference on Learning Theory_ , pages 1912–1947. PMLR, 2022\.
  * Bellec and Zhang (2021) Pierre C. Bellec and Cun-Hui Zhang.  Second-order Stein: SURE for SURE and other applications in high-dimensional inference.  _Ann. Statist._ , 49(4):1864–1903, 2021.  ISSN 0090-5364.
  * Bellec and Zhang (2022a) Pierre C. Bellec and Cun-Hui Zhang.  De-biasing the lasso with degrees-of-freedom adjustment.  _Bernoulli_ , 28(2):713–743, 2022a.  ISSN 1350-7265.
  * Bellec and Zhang (2022b) Pierre C. Bellec and Cun-Hui Zhang.  De-biasing the lasso with degrees-of-freedom adjustment.  _Bernoulli_ , 28(2):713–743, 2022b.  ISSN 1350-7265.
  * Bellec and Zhang (2023) Pierre C. Bellec and Cun-Hui Zhang.  Debiasing convex regularized estimators and interval estimation in linear models.  _Ann. Statist._ , 51(2):391–436, 2023.  ISSN 0090-5364.
  * Bellec et al. (2018) Pierre C. Bellec, Guillaume Lecué, and Alexandre B. Tsybakov.  Slope meets lasso: Improved oracle bounds and optimality.  _Ann. Statist._ , 46(6B):3603–3642, 2018.  ISSN 0090-5364.
  * Berthier et al. (2020) Raphaël Berthier, Andrea Montanari, and Phan-Minh Nguyen.  State evolution for approximate message passing with non-separable functions.  _Inf. Inference_ , 9(1):33–79, 2020.  ISSN 2049-8764.
  * Bickel et al. (2009) Peter J. Bickel, Ya’acov Ritov, and Alexandre B. Tsybakov.  Simultaneous analysis of lasso and dantzig selector.  _Ann. Statist._ , 37(4):1705–1732, 08 2009.  10.1214/08-AOS620.
  * Bogdan et al. (2015) Mał gorzata Bogdan, Ewout van den Berg, Chiara Sabatti, Weijie Su, and Emmanuel J. Candès.  SLOPE—adaptive variable selection via convex optimization.  _Ann. Appl. Stat._ , 9(3):1103–1140, 2015.  ISSN 1932-6157.
  * Bottou (2010) Léon Bottou.  Large-scale machine learning with stochastic gradient descent.  In _Proceedings of COMPSTAT’2010: 19th International Conference on Computational StatisticsParis France, August 22-27, 2010 Keynote, Invited and Contributed Papers_ , pages 177–186. Springer, 2010.
  * Boucheron et al. (2013) Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  _Concentration inequalities: A nonasymptotic theory of independence_.  Oxford University Press, 2013.
  * Bühlmann and Van De Geer (2011) Peter Bühlmann and Sara Van De Geer.  _Statistics for high-dimensional data: methods, theory and applications_.  Springer Science & Business Media, 2011.
  * Bühlmann and Yu (2003) Peter Bühlmann and Bin Yu.  Boosting with the $L_{2}$ loss: regression and classification.  _J. Amer. Statist. Assoc._ , 98(462):324–339, 2003.  ISSN 0162-1459.
  * Cai and Guo (2017) T. Tony Cai and Zijian Guo.  Confidence intervals for high-dimensional linear regression: minimax rates and adaptivity.  _Ann. Statist._ , 45(2):615–646, 2017.  ISSN 0090-5364.
  * Celentano and Montanari (2022) Michael Celentano and Andrea Montanari.  Fundamental barriers to high-dimensional regression with convex penalties.  _Ann. Statist._ , 50(1):170–196, 2022.  ISSN 0090-5364.
  * Celentano et al. (2020) 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\.
  * Celentano et al. (2021) Michael Celentano, Chen Cheng, and Andrea Montanari.  The high-dimensional asymptotics of first order methods with random data.  _arXiv preprint arXiv:2112.07572_ , 2021.
  * Celentano et al. (2023) Michael Celentano, Andrea Montanari, and Yuting Wei.  The Lasso with general Gaussian designs with applications to hypothesis testing.  _Ann. Statist._ , 51(5):2194–2220, 2023.  ISSN 0090-5364.
  * Chandrasekher et al. (2022) Kabir Aladin Chandrasekher, Mengqi Lou, and Ashwin Pananjady.  Alternating minimization for generalized rank one matrix sensing: Sharp predictions from a random initialization.  _arXiv preprint arXiv:2207.09660_ , 2022.
  * Chandrasekher et al. (2023) Kabir Aladin Chandrasekher, Ashwin Pananjady, and Christos Thrampoulidis.  Sharp global convergence guarantees for iterative nonconvex optimization with random data.  _Ann. Statist._ , 51(1):179–210, 2023.  ISSN 0090-5364.
  * Dalalyan et al. (2017) Arnak S Dalalyan, Mohamed Hebiri, and Johannes Lederer.  On the prediction performance of the lasso.  _Bernoulli_ , 23(1):552–581, 2017.
  * Daubechies et al. (2004) Ingrid Daubechies, Michel Defrise, and Christine De Mol.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  _Comm. Pure Appl. Math._ , 57(11):1413–1457, 2004\.  ISSN 0010-3640.
  * Davidson and Szarek (2001) Kenneth R Davidson and Stanislaw J Szarek.  Local operator theory, random matrices and banach spaces.  _Handbook of the geometry of Banach spaces_ , 1(317-366):131, 2001.
  * Donoho and Montanari (2016) David Donoho and Andrea Montanari.  High dimensional robust M-estimation: asymptotic variance via approximate message passing.  _Probab. Theory Related Fields_ , 166(3-4):935–969, 2016.  ISSN 0178-8051.
  * Donoho et al. (2009) David L Donoho, Arian Maleki, and Andrea Montanari.  Message-passing algorithms for compressed sensing.  _Proceedings of the National Academy of Sciences_ , 106(45):18914–18919, 2009.
  * Dudeja et al. (2023) Rishabh Dudeja, Yue M. Lu, and Subhabrata Sen.  Universality of approximate message passing with semirandom matrices.  _The Annals of Probability_ , 51(5):1616–1683, 2023.
  * El Alaoui et al. (2022) Ahmed El Alaoui, Andrea Montanari, and Mark Sellke.  Sampling from the sherrington-kirkpatrick gibbs measure via algorithmic stochastic localization.  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 323–334. IEEE, 2022.
  * El Karoui et al. (2013) Noureddine El Karoui, Derek Bean, Peter J Bickel, Chinghway Lim, and Bin Yu.  On robust regression with high-dimensional predictors.  _Proceedings of the National Academy of Sciences_ , 110(36):14557–14562, 2013.
  * Fan and Li (2001) Jianqing Fan and Runze Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _J. Amer. Statist. Assoc._ , 96(456):1348–1360, 2001.  ISSN 0162-1459.
  * Feng and Zhang (2019) Long Feng and Cun-Hui Zhang.  Sorted concave penalized regression.  _Ann. Statist._ , 47(6):3069–3098, 2019.  ISSN 0090-5364.
  * Friedman et al. (2010) Jerome H. Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _J. Stat. Softw._ , 33(1):1–22, 2010.
  * Gerbelot et al. (2022) Cedric Gerbelot, Emanuele Troiani, Francesca Mignacco, Florent Krzakala, and Lenka Zdeborova.  Rigorous dynamical mean field theory for stochastic gradient descent methods.  _arXiv preprint arXiv:2210.06591_ , 2022.
  * Han and Shen (2023) Qiyang Han and Yandi Shen.  Universality of regularized regression estimators in high dimensions.  _The Annals of Statistics_ , 51(4):1799–1823, 2023.
  * Hastie et al. (2022) Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  _Ann. Statist._ , 50(2):949–986, 2022.  ISSN 0090-5364.
  * Hoppe et al. (2023) Frederik Hoppe, Claudio Mayrink Verdun, Hannah Laus, Felix Krahmer, and Holger Rauhut.  Uncertainty quantification for learned ista.  In _2023 IEEE 33rd International Workshop on Machine Learning for Signal Processing (MLSP)_ , pages 1–6. IEEE, 2023.
  * Hu and Lu (2022) Hong Hu and Yue M Lu.  Universality laws for high-dimensional learning with random features.  _IEEE Transactions on Information Theory_ , 69(3):1932–1964, 2022.
  * Hutchinson (1990) M. F. Hutchinson.  A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines.  _Comm. Statist. Simulation Comput._ , 19(2):433–450, 1990.  ISSN 0361-0918.
  * Javanmard and Montanari (2014) Adel Javanmard and Andrea Montanari.  Confidence intervals and hypothesis testing for high-dimensional regression.  _J. Mach. Learn. Res._ , 15:2869–2909, 2014.  ISSN 1532-4435.
  * Javanmard and Montanari (2018) Adel Javanmard and Andrea Montanari.  Debiasing the Lasso: optimal sample size for Gaussian designs.  _Ann. Statist._ , 46(6A):2593–2622, 2018.  ISSN 0090-5364.
  * Koltchinskii et al. (2011) Vladimir Koltchinskii, Karim Lounici, and Alexandre B. Tsybakov.  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  _Ann. Statist._ , 39(5):2302–2329, 2011.  ISSN 0090-5364.
  * Lei et al. (2018) Lihua Lei, Peter J. Bickel, and Noureddine El Karoui.  Asymptotics for high dimensional regression $M$-estimates: fixed design results.  _Probab. Theory Related Fields_ , 172(3-4):983–1079, 2018.  ISSN 0178-8051.
  * Lou et al. (2024) Mengqi Lou, Kabir Aladin Verchand, and Ashwin Pananjady.  Hyperparameter tuning via trajectory predictions: Stochastic prox-linear methods in matrix sensing.  _arXiv preprint arXiv:2402.01599_ , 2024.
  * Loureiro et al. (2021) Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mezard, and Lenka Zdeborová.  Learning curves of generic features maps for realistic datasets with a teacher-student model.  _Advances in Neural Information Processing Systems_ , 34:18137–18151, 2021.
  * Luo et al. (2023) Yuetian Luo, Zhimei Ren, and Rina Barber.  Iterative approximate cross-validation.  In _International Conference on Machine Learning_ , pages 23083–23102. PMLR, 2023.
  * Miolane and Montanari (2021) Léo Miolane and Andrea Montanari.  The distribution of the Lasso: uniform control over sparse balls and adaptive parameter tuning.  _Ann. Statist._ , 49(4):2313–2335, 2021.  ISSN 0090-5364.
  * Montanari and Saeed (2022) Andrea Montanari and Basil N Saeed.  Universality of empirical risk minimization.  In _Conference on Learning Theory_ , pages 4310–4312. PMLR, 2022\.
  * Montanari and Venkataramanan (2021) Andrea Montanari and Ramji Venkataramanan.  Estimation of low-rank matrices via approximate message passing.  _Ann. Statist._ , 49(1):321–345, 2021.  ISSN 0090-5364.
  * Montanari and Wu (2022) Andrea Montanari and Yuchen Wu.  Statistically optimal first order algorithms: A proof via orthogonalization.  _arXiv preprint arXiv:2201.05101_ , 2022.
  * Nesterov (2003) Yurii Nesterov.  _Introductory lectures on convex optimization: A basic course_ , volume 87.  Springer Science & Business Media, 2003.
  * Nesterov (1983) Yurii Evgen’evich Nesterov.  A method of solving a convex programming problem with convergence rate o$\backslash$bigl(k^2$\backslash$bigr).  In _Doklady Akademii Nauk_ , volume 269, pages 543–547. Russian Academy of Sciences, 1983.
  * Patil et al. (2024) Pratik Patil, Yuchen Wu, and Ryan J Tibshirani.  Failures and successes of cross-validation for early-stopped gradient descent.  _arXiv preprint arXiv:2402.16793_ , 2024.
  * Pesce et al. (2023) Luca Pesce, Florent Krzakala, Bruno Loureiro, and Ludovic Stephan.  Are gaussian data all you need? the extents and limits of universality in high-dimensional generalized linear estimation.  In _International Conference on Machine Learning_ , pages 27680–27708. PMLR, 2023.
  * Raskutti et al. (2014) Garvesh Raskutti, Martin J. Wainwright, and Bin Yu.  Early stopping and non-parametric regression: an optimal data-dependent stopping rule.  _J. Mach. Learn. Res._ , 15:335–366, 2014.  ISSN 1532-4435.
  * Rosset et al. (2004) Saharon Rosset, Ji Zhu, and Trevor Hastie.  Boosting as a regularized path to a maximum margin classifier.  _J. Mach. Learn. Res._ , 5:941–973, 2004.  ISSN 1532-4435.
  * Sheng and Ali (2022) Yue Sheng and Alnur Ali.  Accelerated gradient flow: Risk, stability, and implicit regularization.  _arXiv preprint arXiv:2201.08311_ , 2022.
  * Sur and Candès (2019) Pragya Sur and Emmanuel J Candès.  A modern maximum-likelihood theory for high-dimensional logistic regression.  _Proceedings of the National Academy of Sciences_ , 116(29):14516–14525, 2019.
  * Tan and Bellec (2024) Kai Tan and Pierre C Bellec.  Multinomial logistic regression: Asymptotic normality on null covariates in high-dimensions.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Tan et al. (2022) Kai Tan, Gabriel Romon, and Pierre C Bellec.  Noise covariance estimation in multi-task high-dimensional linear models.  _arXiv preprint arXiv:2206.07256_ , 2022.
  * Thrampoulidis et al. (2018) Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi.  Precise error analysis of regularized $M$-estimators in high dimensions.  _IEEE Trans. Inform. Theory_ , 64(8):5592–5628, 2018.  ISSN 0018-9448.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 58(1):267–288, 1996.
  * van de Geer (2016) Sara van de Geer.  Estimation and testing under sparsity: École d’été de probabilités de saint-flour xlv–2015.  _Lecture Notes in Mathematics_ , 2159, 2016.
  * van de Geer et al. (2014a) Sara van de Geer, Peter Bühlmann, Ya’acov Ritov, and Ruben Dezeure.  On asymptotically optimal confidence regions and tests for high-dimensional models.  _Ann. Statist._ , 42(3):1166–1202, 2014a.  ISSN 0090-5364.
  * van de Geer et al. (2014b) Sara van de Geer, Peter Bühlmann, Ya’acov Ritov, and Ruben Dezeure.  On asymptotically optimal confidence regions and tests for high-dimensional models.  _Ann. Statist._ , 42(3):1166–1202, 2014b.  ISSN 0090-5364.
  * Villani et al. (2009) Cédric Villani et al.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Wainwright (2009) Martin J. Wainwright.  Sharp thresholds for high-dimensional and noisy sparsity recovery using $\ell_{1}$-constrained quadratic programming (Lasso).  _IEEE Trans. Inform. Theory_ , 55(5):2183–2202, 2009.  ISSN 0018-9448.
  * Ye and Zhang (2010) Fei Ye and Cun-Hui Zhang.  Rate minimaxity of the Lasso and Dantzig selector for the $\ell_{q}$ loss in $\ell_{r}$ balls.  _J. Mach. Learn. Res._ , 11:3519–3540, 2010.  ISSN 1532-4435.
  * Yuan and Lin (2006) Ming Yuan and Yi Lin.  Model selection and estimation in regression with grouped variables.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 68(1):49–67, 2006.  ISSN 1369-7412.
  * Zhang (2010) Cun-Hui Zhang.  Nearly unbiased variable selection under minimax concave penalty.  _Ann. Statist._ , 38(2):894–942, 2010.  ISSN 0090-5364.
  * Zhang and Zhang (2014a) Cun-Hui Zhang and Stephanie S. Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  _J. R. Stat. Soc. Ser. B. Stat. Methodol._ , 76(1):217–242, 2014a.  ISSN 1369-7412.
  * Zhang and Zhang (2014b) Cun-Hui Zhang and Stephanie S. Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  _J. R. Stat. Soc. Ser. B. Stat. Methodol._ , 76(1):217–242, 2014b.  ISSN 1369-7412.
  * Ziemer (1989) William P Ziemer.  _Weakly differentiable functions: Sobolev spaces and functions of bounded variation_ , volume 120.  Springer-Verlag New York, 1989.
  * Zou and Hastie (2005) Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 67(2):301–320, 2005.  ISSN 1369-7412.

##  Appendix A Additional simulation results

In this section, we provide more simulation results for the other two
scenarios $(n,p)=(1200,500)$ and $(n,p)=(1200,1200)$, which are similar to the
results for $(n,p)=(1200,1500)$ presented in Section 3.

###  A.1 Under parametrization: $(n,p)=(1200,500)$

####  A.1.1 GD and AGD

For GD and AGD applied to solve least-squares problem, when
$(n,p)=(1200,500)$, we know that the iterate ${\widehat{\bm{b}}}^{t}$
converges to the ordinary least-squares estimate
$(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}$ as $t\to\infty$. Thus we are
able to compute the limiting risk $r_{\infty}$ for GD and AGD. We present the
simulation results for GD and AGD in Figure 6 and Figure 7, respectively.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 6: Risk curves and qq-plots of z-score of GD for $(n,p)=(1200,500)$.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 7: Risk curves and qq-plots of z-score of AGD for $(n,p)=(1200,500)$.

####  A.1.2 ISTA and FISTA

We present the simulation results for ISTA and FISTA in Figure 8 and Figure 9,
respectively.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 8: Risk curves and qq-plots of z-score of ISTA for $(n,p)=(1200,500)$.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 9: Risk curves and qq-plots of z-score of FISTA for $(n,p)=(1200,500)$.

####  A.1.3 LQA

We present the simulation results for LQA in Figure 10.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 10: Risk curves and qq-plots of z-score of LQA for $(n,p)=(1200,500)$.

###  A.2 Equal parametrization: $(n,p)=(1200,1200)$

In this section, we present the simulation results for the equal
parametrization scenario $(n,p)=(1200,1200)$.

####  A.2.1 GD and AGD

For GD and AGD, we know that the iterates converge to the min-norm least-
squares estimator (3.3) as $t\to\infty$ (Hastie et al., 2022, Proposition 1).
Under $n=p=1200$, we know that the risk of the min-norm least-squares
estimator is infinite (Hastie et al., 2022), i.e., $r_{\infty}=+\infty$. So
here we do not plot the horizontal line for $r_{\infty}$.

We present the simulation results for GD and AGD in Figure 11 and Figure 12,
respectively.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 11: Risk curves and qq-plots of z-score of GD for $(n,p)=(1200,1200)$.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 12: Risk curves and qq-plots of z-score of AGD for $(n,p)=(1200,1200)$.

####  A.2.2 ISTA and FISTA

We present the simulation results for ISTA and FISTA in Figure 13 and Figure
14, respectively.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 13: Risk curves and qq-plots of z-score of ISTA for
$(n,p)=(1200,1200)$.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 14: Risk curves and qq-plots of z-score of FISTA for
$(n,p)=(1200,1200)$.

####  A.2.3 LQA

We present the simulation results for LQA in Figure 15.

(a) Risk curves versus iteration number.

(b) Q-Q plots of z-score for $\lambda=0.1$.

Figure 15: Risk curves and qq-plots of z-score of LQA for $(n,p)=(1200,1200)$.

##  Appendix B Preliminary

###  B.1 Notation

Notation and definitions that will be used in the rest of the paper are given
here. Regular variables like $a,b,...$ refer to scalars, bold lowercase
letters such as $\bm{a},\bm{b},...$ represent vectors, and bold uppercase
letters like $\bm{A},\bm{B},...$ indicate matrices. Let
$[n]=\\{1,2,\ldots,n\\}$ for all $n\in{\mathbb{N}}$. The vectors
$\bm{e}_{i}\in{\mathbb{R}}^{n},\bm{e}_{j}\in{\mathbb{R}}^{p},\bm{e}_{t}\in{\mathbb{R}}^{T}$
denote the canonical basis vector of the corresponding index. For a real
vector $\bm{a}\in{\mathbb{R}}^{p}$, $\left\lVert\bm{a}\right\rVert$ denotes
its Euclidean norm. For any matrix $\bm{A}$, $\bm{A}^{\dagger}$ is its
Moore–Penrose inverse; $\left\lVert\bm{A}\right\rVert_{\rm{F}}$
,$\left\lVert\bm{A}\right\rVert_{\rm{op}}$,
$\left\lVert\bm{A}\right\rVert_{*}$ denote its Frobenius, operator and nuclear
norm, respectively. Let $\bm{A}\otimes\bm{B}$ be the Kronecker product of
$\bm{A}$ and $\bm{B}$. For $\bm{A}$ symmetric, $\phi_{\min}(\bm{A})$ and
$\phi_{\max}(\bm{A})$ denote its smallest and largest eigenvalues,
respectively. Let ${\bf 1}_{n}$ denote the all-ones vector in
${\mathbb{R}}^{n}$, and $\bm{I}_{n}$ denote the identity matrix of size $n$.
For a mapping
${\mathbb{R}}^{p}\to{\mathbb{R}}^{n}:\bm{x}\mapsto\bm{f}(\bm{x})$, we denote
its Jacobian by $\pdv{\bm{f}(\bm{x})}{\bm{x}}\in{\mathbb{R}}^{n\times p}$,
i.e., $(\pdv{\bm{f}(\bm{x})}{\bm{x}})_{i,j}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\pdv{\bm{f}_{i}(\bm{x})}{\bm{x}_{j}}$ for all $i\in[n],j\in[p]$.
For a random sequence $\xi_{n}$, we write $\xi_{n}=O_{P}(a_{n})$ if
$\xi_{n}/a_{n}$ is stochastically bounded. For two scalars $a,b$, $a\vee b$
denotes the maximum of $a$ and $b$.

Let ${\mathbb{I}}(\Omega)$ denote the indicator function of event $\Omega$. It
takes the value 1 if the event $\Omega$ occurs and 0 otherwise. Let $C$
represent an absolute constant. Additionally, we use $C(\zeta,\gamma)$ to
denote a positive constant that only depends on $\tau$ and $\gamma$.
Similarly, we extend this notation to $C(\zeta,\gamma,\kappa,\ldots)$,
representing positive constants dependent only on $\tau$, $\gamma$, $\kappa$,
and so forth. The exact value of these constants may vary from place to place.
We write $a\lesssim b$ if $a\leq Cb$ for some absolute constant $C$.

Let $\mathsf{N}(\mu,\sigma^{2})$ denote the univariate Gaussian distribution
with mean $\mu$ and variance $\sigma^{2}$, and
$\mathsf{N}_{k}(\bm{\mu},\bm{\Sigma})$ denote the k-dimensional Gaussian
distribution with mean $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$.

In this paper, we adopt the convention use of the expectation and conditional
expectation. For a random variable $X$, the expressions ${\mathbb{E}}X$,
${\mathbb{E}}(X)$, and ${\mathbb{E}}[X]$ all refer to the expectation of $X$.
Similarly, ${\mathbb{E}}X^{2}$, ${\mathbb{E}}(X^{2})$, and
${\mathbb{E}}[X^{2}]$ mean that we first square the random variable $X$ and
then compute its expectation. In contrast, expressions like
${\mathbb{E}}(X)^{2}$ and ${\mathbb{E}}[X]^{2}$ denote a different operation;
here, we first calculate the expectation of $X$ and then square the result.
This convention extends consistently to other power operations and conditional
expectations throughout the paper.

###  B.2 Review of Kronecker product

In this section, we review the definitions and properties of the Kronecker
product. Let $\bm{A}$ be an $m\times n$ matrix, and $\bm{B}$ be a $p\times q$
matrix. The Kronecker product of $\bm{A}$ and $\bm{B}$, denoted by
$\bm{A}\otimes\bm{B}$, is defined as:

$\bm{A}\otimes\bm{B}=\begin{bmatrix}a_{11}\bm{B}&\cdots&a_{1n}\bm{B}\\\
\vdots&\ddots&\vdots\\\ a_{m1}\bm{B}&\cdots&a_{mn}\bm{B}\end{bmatrix},$

where $a_{ij}$ represents the entry in the $i$-th row and $j$-th column of
matrix $\bm{A}$. Below we list a few properties of the Kronecker product that
will be useful in our proofs.

  1. (i)

The mixed-product property: If $\bm{A},\bm{B},\bm{C}$ and $\bm{D}$ are
matrices of such size that the matrix products $\bm{A}\bm{C}$ and
$\bm{B}\bm{D}$ make sense, then

(B.1)
$(\bm{A}\otimes\bm{B})(\bm{C}\otimes\bm{D})=(\bm{A}\bm{C})\otimes(\bm{B}\bm{D})\in{\mathbb{R}}^{mp\times
nq}.$

  2. (ii)

Inverse property:

(B.2)
$\displaystyle(\bm{A}\otimes\bm{B})^{-1}=\bm{A}^{-1}\otimes\bm{B}^{-1}\text{
and }(\bm{A}\otimes\bm{B})^{\dagger}=\bm{A}^{\dagger}\otimes\bm{B}^{\dagger},$

where $\bm{A}^{\dagger}$ means the Moore-Penrose inverse of $\bm{A}$.

  3. (iii)

Trace property:

(B.3)
$\displaystyle\operatorname{Tr}(\bm{A}\otimes\bm{B})=\operatorname{Tr}(\bm{A})\operatorname{Tr}(\bm{B}).$

  4. (iv)

Norm property:

(B.4)
$\displaystyle\lVert\bm{A}\otimes\bm{B}\rVert_{\rm{op}}=\lVert\bm{A}\rVert_{\rm{op}}\lVert\bm{B}\rVert_{\rm{op}}\text{
and
}\lVert\bm{A}\otimes\bm{B}\rVert_{\rm{F}}=\lVert\bm{A}\rVert_{\rm{F}}\lVert\bm{B}\rVert_{\rm{F}}.$

  5. (v)

Relationship with vectorization operator: If the matrix product
$\bm{A}\bm{B}\bm{C}$ makes sense, then

(B.5)
$\displaystyle\operatorname{\mathbf{vec}}(\bm{A}\bm{B}\bm{C})=(\bm{C}^{\top}\otimes\bm{A})\operatorname{\mathbf{vec}}(\bm{B}),$

where $\operatorname{\mathbf{vec}}(\cdot)$ is the vectorization operator such
that $\operatorname{\mathbf{vec}}(\bm{B})$ is the vector obtained by stacking
vertically the columns of $\bm{B}$ on top of one another.

##  Appendix C Derivative formula

In this section, we provide the derivative formulas for the iteration function
$\bm{g}_{t}$ in (2.5). These derivative formulas will be important ingredient
in the proof of the main results in the next sections. Recall
${\widehat{\bm{b}}}^{1},{\widehat{\bm{b}}}^{2},\ldots,{\widehat{\bm{b}}}^{T}$
represent the first $T$ iterates of an algorithm using the iteration function
$\bm{g}_{t}$ in (2.8), and
$\bm{v}^{t}=\tfrac{1}{n}\bm{X}^{\top}(\bm{y}-\bm{X}{\widehat{\bm{b}}}^{t})$ as
in (2.1). The error matrix and residual matrix for the first $T$ iterates
${\widehat{\bm{b}}}^{1},{\widehat{\bm{b}}}^{2},\ldots,{\widehat{\bm{b}}}^{T}$
are defined as

(C.1)
$\bm{H}=[{\widehat{\bm{b}}}^{1}-\bm{b}^{*},\ldots,{\widehat{\bm{b}}}^{T}-\bm{b}^{*}]\in{\mathbb{R}}^{p\times
T},\quad\bm{F}=[\bm{y}-\bm{X}{\widehat{\bm{b}}}^{1},\ldots,\bm{y}-\bm{X}{\widehat{\bm{b}}}^{T}]\in{\mathbb{R}}^{n\times
T}.$

We also define
${\widehat{\bm{B}}}=[{\widehat{\bm{b}}}^{1},\ldots,{\widehat{\bm{b}}}^{T}]$ is
a matrix of size $p\times T$, $\bm{B}^{*}=[\bm{b}^{*},\ldots,\bm{b}^{*}]$,
$\bm{E}=[{\bm{\varepsilon}},...,{\bm{\varepsilon}}]$ and
$\bm{Y}=[\bm{y},\ldots,\bm{y}]$ are matrices formed by repeating the vectors
$\bm{b}^{*}$, ${\bm{\varepsilon}}$, and $\bm{y}$ column-wise $T$ times.
Therefore, we have $\bm{H}={\widehat{\bm{B}}}-\bm{B}^{*}$, and
$\bm{F}=\bm{Y}-\bm{X}{\widehat{\bm{B}}}$. Recall also the matrices
$\bm{D}_{t,s},\bm{J}_{t,s}$ in (2.9) and the large matrices
$\mathcal{D},\mathcal{J}\in{\mathbb{R}}^{pT\times pT}$ in (2.10).

######  Lemma C.1 (Derivative of iterates).

For $({\widehat{\bm{b}}}^{1},...,{\widehat{\bm{b}}}^{T})$ in (2.8), for almost
every $(\bm{X},{\bm{\varepsilon}})$,

(C.2) $\displaystyle\pdv{{\widehat{\bm{b}}}^{t}}{x_{ij}}$
$\displaystyle=n^{-1}(\bm{e}_{t}^{\top}\otimes\bm{I}_{p}){\mathcal{M}}^{-1}{\mathcal{D}}\Bigl{[}\Bigl{(}(\bm{F}^{\top}\bm{e}_{i})\otimes\bm{e}_{j}\Bigr{)}-\Bigl{(}(\bm{H}^{\top}\bm{e}_{j})\otimes(\bm{X}^{\top}\bm{e}_{i})\Bigr{)}\Bigr{]},$
(C.3) $\displaystyle\pdv{{\widehat{\bm{b}}}^{t}}{\epsilon_{l}}$
$\displaystyle=n^{-1}(\bm{e}_{t}^{\top}\otimes\bm{I}_{p}){\mathcal{M}}^{-1}{\mathcal{D}}({\mathbf{1}}_{T}\otimes\bm{X}^{\top}\bm{e}_{l}),$

where
$\bm{e}_{i},\bm{e}_{l}\in{\mathbb{R}}^{n},\bm{e}_{j}\in{\mathbb{R}}^{p},\bm{e}_{t}\in{\mathbb{R}}^{T}$
and
${\mathcal{M}}=\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}$.

The derivative results in Lemma C.1 directly imply the following by the
product rule.

######  Corollary C.2 (Derivative of residuals).

Let $F_{lt}$ be the $(l,t)$-th entry of $\bm{F}$, i.e.,
$F_{lt}=\bm{e}_{l}^{\top}\bm{F}\bm{e}_{t}$. Under the conditions of Lemma C.1,
for each $i,l\in[n]$, $j\in[p]$, $t\in[T]$, we have

(C.4) $\frac{\partial F_{lt}}{\partial
x_{ij}}=D_{ij}^{lt}+\Delta_{ij}^{lt}\quad{and}\quad\pdv{\bm{F}\bm{e}_{t}}{\epsilon_{i}}=\bm{e}_{i}-n^{-1}(\bm{e}_{t}^{\top}\otimes\bm{X}){\mathcal{M}}^{-1}{\mathcal{D}}({\mathbf{1}}_{T}\otimes\bm{X}^{\top}\bm{e}_{i}),$

where the expressions of $D_{ij}^{lt}$ and $\Delta_{ij}^{lt}$ are given by

(C.5) $\displaystyle D_{ij}^{lt}$
$\displaystyle=-(\bm{e}_{t}^{\top}\otimes\bm{e}_{l}^{\top})[\bm{I}_{nT}-n^{-1}(\bm{I}_{T}\otimes\bm{X}){\mathcal{M}}^{-1}{\mathcal{D}}(\bm{I}_{T}\otimes\bm{X}^{\top})]((\bm{H}^{\top}\bm{e}_{j})\otimes\bm{e}_{i}),$
(C.6) $\displaystyle\Delta_{ij}^{lt}$
$\displaystyle=-n^{-1}(\bm{e}_{t}^{\top}\otimes\bm{e}_{l}^{\top})(\bm{I}_{T}\otimes\bm{X}){\mathcal{M}}^{-1}{\mathcal{D}}\bigl{(}\bm{F}^{\top}\otimes\bm{I}_{p}\bigr{)}(\bm{e}_{i}\otimes\bm{e}_{j}).$

Here ${\mathcal{D}}$ is defined in (2.7) and
${\mathcal{M}}=\bm{I}_{pT}+{\mathcal{D}}(\bm{I}_{T}\otimes\tfrac{\bm{X}^{\top}\bm{X}}{n})-{\mathcal{J}}$
is defined in Lemma C.1. It immediately follows that

(C.7)
$\displaystyle\sum_{i=1}^{n}D_{ij}^{it}=-\bm{e}_{t}^{\top}(n\bm{I}_{T}-\smash{\widehat{\mathbf{A}}})\bm{H}^{\top}\bm{e}_{j}\quad\text{and}\quad\sum_{i=1}^{n}\sum_{t=1}^{T}D_{ij}^{it}\bm{e}_{t}^{\top}=-\bm{e}_{j}^{\top}\bm{H}(n\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}^{\top}),$
(C.8)
$\displaystyle\operatorname{Tr}\Bigl{(}\pdv{\bm{X}{\widehat{\bm{b}}}^{t}}{{\bm{\varepsilon}}}\Bigr{)}=\sum_{l=1}^{n}\pdv{\bm{e}_{l}^{\top}\bm{X}{\widehat{\bm{b}}}^{t}}{\epsilon_{l}}=\bm{e}_{t}^{\top}\smash{\widehat{\mathbf{A}}}{\bf
1}_{T}\quad\text{ and
}\quad\operatorname{Tr}\Bigl{(}\pdv{\bm{F}\bm{e}_{t}}{{\bm{\varepsilon}}}\Bigr{)}=\bm{e}_{t}^{\top}(n\bm{I}_{T}-\smash{\widehat{\mathbf{A}}}){\bf
1}_{T}.$

###### Proof of Lemma C.1.

By composition of locally Lipschitz functions, the map
$(\bm{X},{\bm{\varepsilon}})\mapsto{\widehat{\bm{b}}}^{t}$ is locally
Lipschitz and thus differentiable almost everywhere. Let
$(\dot{\bm{X}},\dot{\bm{\varepsilon}})$ be a perturbation direction. We now
compute the directional derivative of ${\widehat{\bm{b}}}^{t}$ with respect to
this direction by taking, for any function $F({\bm{\varepsilon}},\bm{X})$ the
limit
$u^{-1}(F({\bm{\varepsilon}}+u\dot{\bm{\varepsilon}},\bm{X}+u\dot{\bm{X}})-F({\bm{\varepsilon}},\bm{X}))$
as $u\to 0$ and call the limit $\dot{F}$. For ${\widehat{\bm{b}}}^{t}$, we
denote the corresponding directional derivative by $\dot{\bm{b}}^{t}$, and for
$\bm{v}^{t}$ by $\dot{\bm{v}}^{t}$. By definition of ${\widehat{\bm{b}}}^{t}$
in (2.8), we have almost surely by the chain rule (if necessary with
modification from (Ambrosio and Dal Maso, 1990, Corollary 3.2) as explained
after (2.11)),

(C.9)
$\displaystyle\dot{\bm{b}}^{t}=\sum_{s=1}^{t-1}\bm{J}_{t,s}\dot{\bm{b}}^{s}+\bm{D}_{t,s}\dot{\bm{v}}^{s}.$

Since
$\bm{v}^{t}=n^{-1}\bm{X}^{\top}\bm{F}\bm{e}_{t}=n^{-1}\bm{X}^{\top}({\bm{\varepsilon}}-\bm{X}\bm{H}\bm{e}_{t})$
and $\dot{\bm{H}}\bm{e}_{t}=\dot{\bm{b}}^{t}$, we have by the product rule

(C.10)
$\dot{\bm{v}}^{t}=-\tfrac{1}{n}\bm{X}^{\top}\bm{X}\dot{\bm{b}}^{t}+\underbrace{\tfrac{1}{n}\big{[}\dot{\bm{X}}^{\top}\bm{F}\bm{e}_{t}+\bm{X}^{\top}(\dot{\bm{\varepsilon}}-\dot{\bm{X}}\bm{H}\bm{e}_{t})\big{]}}_{\bm{a}^{t}}$

where $\bm{a}^{t}\in{\mathbb{R}}^{p}$. Substituting the expressions for
$\dot{\bm{v}}^{s}$ into (C.9), we have

(C.11)
$\displaystyle\dot{\bm{b}}^{t}=\sum_{s=1}^{t-1}(\bm{J}_{t,s}-\bm{D}_{t,s}\tfrac{\bm{X}^{\top}\bm{X}}{n})\dot{\bm{b}}^{s}+\bm{D}_{t,s}\bm{a}^{s}.$

For the initial condition, we have $\dot{\bm{b}}^{1}={\mathbf{0}}$ since
${\widehat{\bm{b}}}^{1}$ is a constant independent of $\bm{X}$. Similarly
$\bm{D}_{2,0}=\bm{J}_{2,0}={\mathbf{0}}$ since ${\widehat{\bm{b}}}^{0}$ and
$\bm{v}^{0}$ are set as constant vector ${\mathbf{0}}$. Therefore we can write
(C.11) as a linear system of size $pT$:

$\displaystyle\underbrace{\begin{bmatrix}\bm{I}_{p}&&&&\\\
\bm{D}_{2,1}\tfrac{\bm{X}^{\top}\bm{X}}{n}-\bm{J}_{2,1}&\bm{I}_{p}&&&\\\
\bm{D}_{3,1}\tfrac{\bm{X}^{\top}\bm{X}}{n}-\bm{J}_{3,1}&\bm{D}_{3,2}\tfrac{\bm{X}^{\top}\bm{X}}{n}-\bm{J}_{3,2}&\bm{I}_{p}&&\\\
\vdots&\ddots&\ddots&\ddots&\\\
\bm{D}_{T,1}\tfrac{\bm{X}^{\top}\bm{X}}{n}-\bm{J}_{T,1}&\cdots&\bm{D}_{T,T-2}\tfrac{\bm{X}^{\top}\bm{X}}{n}-\bm{J}_{T,T-2}&\bm{D}_{T,T-1}\tfrac{\bm{X}^{\top}\bm{X}}{n}-\bm{J}_{T,T-1}&\bm{I}_{p}\end{bmatrix}}_{{\mathcal{M}}}\underbrace{\begin{bmatrix}\dot{\bm{b}}^{1}\\\
\dot{\bm{b}}^{2}\\\ \dot{\bm{b}}^{3}\\\ \vdots\\\ \dot{\bm{b}}^{T}\\\
\end{bmatrix}}_{\operatorname{\mathbf{vec}}(\dot{\bm{B}})}$