Implicit Sparse Regularization: The Impact of Depth and Early Stopping
Year: 2021

Abstract:
In this paper, we study the implicit bias of gradient descent for sparse regression. We extend results on regression with quadratic parametrization, which amounts to depth-2 diagonal linear networks, to more general depth-N networks, under more realistic settings of noise and correlated designs. We show that early stopping is crucial for gradient descent to converge to a sparse model, a phenomenon that we call implicit sparse regularization. This result is in sharp contrast to known results for noiseless and uncorrelated-design cases. We characterize the impact of depth and early stopping and show that for a general depth parameter N , gradient descent with early stopping achieves minimax optimal sparse recovery with sufficiently small initialization w0 and step size ⌘. In particular, we show that increasing depth enlarges the scale of working initialization and the early-stopping window so that this implicit sparse regularization effect is more likely to take place.

text: In this paper, we study the implicit bias of gradient descent for sparse regression. We extend results on regression with quadratic parametrization, which amounts to depth-2 diagonal linear networks, to more general depth-N networks, under more realistic settings of noise and correlated designs. We show that early stopping is crucial for gradient descent to converge to a sparse model, a phenomenon that we call implicit sparse regularization. This result is in sharp contrast to known results for noiseless and uncorrelated-design cases. We characterize the impact of depth and early stopping and show that for a general depth parameter N , gradient descent with early stopping achieves minimax optimal sparse recovery with sufficiently small initialization w0 and step size ⌘. In particular, we show that increasing depth enlarges the scale of working initialization and the early-stopping window so that this implicit sparse regularization effect is more likely to take place.

heading: 1 Introduction
text: Motivation. Central to recent research in learning theory is the insight that the choice of optimization algorithms plays an important role in model generalization [1, 2, 3]. A widely adopted view is that (stochastic) gradient descent — the most popular optimization algorithm in machine learning — exhibits some implicit form of regularization. Indeed for example, in the classical under-determined least squares setting, gradient descent (with small step size) starting from the origin converges to the model with minimum Euclidean norm. Similar implicit biases are also observed in deep neural network training in which the networks typically have many more parameters than the sample size. There, gradient descent without explicit regularization finds solutions that not only interpolate the training data points but also generalize well on test sets [2, 4, 5, 6, 7].
This insight, combined with the empirical success stories of deep learning, has sparked significant interest among theoretical researchers to rigorously understand implicit regularization. The majority of theoretical results focus on well-understood problems such as regression with linear models [8, 9, 10, 11, 12, 13] and matrix factorization [14, 15, 16, 17], and show that the parametrization (or architecture) of the model plays a crucial role. For the latter, Gunasekar et al. [15] conjectured that gradient descent on factorized matrix representations converges to the solution with minimum nuclear norm. The conjecture was partially proved by Li et al. [16] under the Restricted Isometry Property (RIP) and the absence of noise. Arora et al. [17] further show the same nuclear-norm implicit bias using depth-N linear networks (i.e., the matrix variable is factorized into N components).
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Parallel work on nonlinear models and classification [5, 9] has shown that gradient descent biases the solution towards the max-margin/minimum `2-norm solutions over separable data. The scale of initialization in gradient descent leads to two learning regimes (dubbed “kernel” and “rich”) in linear networks [18], shallow ReLU networks [19] and deep linear classifiers [20]. Li et al. [21] showed that depth-2 network requires an exponentially small initialization, whereas depth-N network (N 3) only requires a polynomial small initialization, to obtain low-rank solution in matrix factorisation. Woodworth et al. [18] obtained a similar result for high dimensional sparse regression.
The trend in the large majority of the above works has been to capture implicit regularization of gradient descent using some type of norm with respect to the working parametrization [18, 22, 23, 24]. On the other hand, progress on understanding the trajectory of gradient descent has been somewhat more modest. [11, 12] study the sparse regression problem using quadratic and Hadamard parametrization respectively and show that gradient descent with small initialization and careful early stopping achieves minimax optimal rates for sparse recovery. Unlike [16, 18] that study noiseless settings and require no early stopping, [11, 12] mathematically characterize the role of early stopping and empirically show that it may be necessary to prevent gradient descent from over-fitting to the noise. These works suggest that the inductive bias endowed by gradient descent may be influenced not only by the choice of parametrization, but also algorithmic choices such as initialization, learning rate, and the number of iterations. However, our understanding of such gradient dynamics is incomplete, particularly in the context of deep architectures; see Table 1 for some comparisons.
Contributions. Our focus in this paper is the implicit regularization of (standard) gradient descent for high dimensional sparse regression, namely implicit sparse regularization. Let us assume a groundtruth sparse linear model and suppose we observe n noisy samples (xi, yi), such that y = Xw⇤+⇠; a more formal setup is given in Section 2. Using the samples, we consider gradient descent on a squared loss kXw yk2 with no explicit sparsity regularization. Instead, we write the parameter vector w in the form w = uN vN with N 2. Now, the regression function f(x,u,v) = hx,uN vN i can be viewed as a depth-N diagonal linear network [18]. Minimizing the (now non-convex) loss over u and v with gradient descent is then analogous to training this depth-N network.
Our main contributions are the following. We characterize the impact of both the depth and early stopping for this non-convex optimization problem. Along the way, we also generalize the results of [11] for N > 2. We show that under a general depth parameter N and an incoherence assumption on the design matrix, gradient descent with early stopping achieves minimax optimal recovery with sufficiently small initialization w0 and step size ⌘. The choice of step size is of order O(1/N2). Moreover, the upper bound of the initialization, as well as the early-stopping window, increase with N , suggesting that depth leads to a more accessible generalizable solution on gradient trajectories.
Techniques. At a high level, our work continues the line of work on implicit bias initiated in [11, 12, 18] and extends it to the deep setting. Table 1 highlights key differences between our work and [11, 13, 18]. Specifically, Woodworth et al. [18] study the interpolation given by the gradient flow of the squared-error loss function. Vaskevicius et al. [11] analyze the finite gradient descent and characterize the implicit sparse regularization on the recovery of true parameters with N = 2. Lastly, Gissin et al. [13] discover the incremental learning dynamic of gradient flow for general N but in an idealistic model setting where u ⌫ 0,v = 0, ⇠ = 0, uncorrelated design and with infinitely many samples.
At first glance, one could attempt a straightforward extension of the proof techniques in [11] to general settings of N > 2. However, this turns out to be very challenging. Consider even the simplified case where the true model w? is non-negative, the design matrix is unitary (i.e., n 1XTX = I), and the noise is absent (⇠ = 0); this is the setting studied in [13]. For each entry wi of w, the tth iterate of gradient descent over the depth-N reparametrized model is given by:
wi,t+1 = wi,t ⇣ 1 + w 1 2N i,t (w ? i wi,t) ⌘N ,
which is no longer a simple multiplicative update. As pointed out in [13] (see their Appendix C), the recurrence relation is not analytically solvable due to the presence of the (pesky) term w1 2 N
i,t when
N > 2. Moreover, this extra term w1 2 N
i,t leads to widely divergent growth rates of weights with different magnitudes, which further complicates analytical bounds. To resolve this and rigorously analyze the dynamics for N > 2, we rely on a novel first order, continuous approximation to study growth rates without requiring additional assumptions on gradient flow, and carefully bound the approximation error due to finite step size; see Section 4.

heading: 2 Setup
text: Sparse regression/recovery. Let w? 2 Rp be a p-dimensional sparse vector with k non-zero entries. Assume that we observe n data points (xi, yi) 2 Rp⇥R such that yi = hxi,w?i+⇠i for i = 1, . . . , n, where ⇠ = (⇠1, . . . , ⇠n) is the noise vector. We do not assume any particular scaling between the number of observations n and the dimension p. Due to the sparsity of w?, however, we allow n ⌧ p.
The linear model can be expressed in the matrix-vector form:
y = Xw? + ⇠, (1)
with the n⇥p design matrix X = [x>1 , . . . ,x>n ]>, where xi denotes the ith row of X. We also denote X = [X1, . . . ,Xp], where Xi denotes the ith column of X.
The goal of sparse regression is to estimate the unknown, sparse vector w? from the observations. Over the past two decades, this problem has been a topic of active research in statistics and signal processing [25]. A common approach to sparse regression is penalized least squares with sparsityinduced regularization such as `0 or `1 penalties/constraints, leading to several well-known estimators [25, 26, 27] and algorithms [28, 29]. Multiple estimators enjoy optimal statistical and algorithmic recovery guarantees under some conditions of the design matrix X (e.g., RIP [30]) and the noise ⇠.
We deviate from the standard penalized least squares formulation and instead learn w⇤ via a polynomial parametrization: w = uN vN , u,v 2 Rp, where N 2 and zN = [zN1 , . . . , zNp ]> for any z = [z1, . . . , zN ]> 2 Rp. The regression function f(x,u,v) = hx,uN vN i induced by such a parametrization is equivalent to a N -layer diagonal linear network [18] with 2p hidden neurons and the diagonal weight matrix shared across all layers.
Given the data {X,y} observed in (1), we analyze gradient descent with respect to the new parameters u and v over the mean squared error loss without explicit regularization:
L(u,v) = 1
n
X(uN vN ) y 2 2 , u,v 2 Rp.
Even though the loss function yields the same value for the two parametrizations, L(u,v) is nonconvex in u and v. Unlike several recent studies in implicit regularization for matrix factorization and regression [16, 18, 13], we consider the noisy setting, which is more realistic and leads to more insights into the bias induced during the optimization. Because of noise, the loss evaluated at the ground truth (i.e., any u,v such that w? = uN vN ) is not necessarily zero or even minimal.
Gradient descent. The standard gradient descent update over L(u,v) reads as: u0 = v0 = ↵1,
(ut+1,vt+1) = (ut,vt) ⌘ @L(ut,vt)
@(ut,vt) , t = 0, 1, . . . . (2)
Here, ⌘ > 0 is the step size and ↵ > 0 is the initialization of u,v. In general, we analyze the algorithm presented in (2), and at each step t, we can estimate the signal of interest by simply calculating wt = uNt vNt . We consider constant initialization for simplicity sake. Our results apply for random initialization concentrating on a small positive region with a probabilistic statement.
Vaskevicius et al. [11] establish the implicit sparse regularization of gradient descent for N = 2 and show minimax optimal recovery, provided sufficiently small ↵ and early stopping. Our work aims to generalize that result to N > 2 and characterize the role of N in convergence.
Notation. We define S = {i 2 {1, . . . , p} : w?i 6= 0} and Sc = {1, . . . , p}\S. The largest and smallest absolute value on the support is denoted as w?max = maxi2S |w?i | and w?min = mini2S |w?i |. We use 1 to denote the vector of all ones and 1S denotes the vector whose elements on S are all one and 0 otherwise. Also, denotes coordinate-wise multiplication. We denote st = 1S wt and et = 1Sc wt meaning the signal part and error part at each time step t. We use ^ and _ to denote the pointwise maximum and minimum. The coordinate-wise inequalities are denoted as <. We denote inequalities up to multiplicative absolute constants by ., which means that they do not depend on any parameters of the problem. Definition 1. Let X 2 Rn⇥p be a matrix with `2-normalized columns X1, . . . ,Xp, i.e., kXik2 = 1 for all i. The coherence µ = µ(X) of the matrix X is defined as
µ := max 1i 6=jp |hXi,Xji|.
The matrix X is said to be satisfying µ-incoherence.
The coherence is a measure for the suitability of the measurement matrix in compressive sensing [31]. In general, the smaller the coherence, the better the recovery algorithms perform. There are multiple ways to construct a sensing matrix with low-incoherence. One of them is based on the fact that sub-Gaussian matrices satisfy low-incoherence property with high probability [32, 33]. In contrast to the coherence, the Restricted Isometry Property (RIP) is a powerful performance measure for guaranteeing sparse recovery and has been widely used in many contexts. However, verifying the RIP for deterministically constructed design matrices is NP-hard. On the other hand, coherence is a computationally tractable measure and its use in sparse regression is by now classical [33, 34]. Therefore, in contrast with previous results [11] (which assumes RIP), the assumptions made in our main theorems are verifiable in polynomial time.

heading: 3 Main Results
text: We now introduce several quantities that are relevant for our main results. First, the condition number r := w?max/w ? min plays an important role when we work on the incoherence property of the design matrix. Next, we require an upper bound on the initialization ↵, which depends on the following terms:
(w?max, w ? min, ✏, N) :=
✓ 1
8
◆2/(N 2) ^ (w?max) (N 2)/N
log w ? max ✏
!2/(N 2) ^ (w?min) (N 2)/N
log w?min
✏
!4/(N 2) ,
(w?min, N) := (2 2 N 2 N ) 1 N 2 (w?min) 1 N ^ 2 3 N (2 1 N 1) 1 N 2 (w?min) 1 N .
Finally, define
⇣ := 1
5 w
? min _
200
n
XT⇠ 1 _ 200✏.
We are now ready to state the main theorem: Theorem 1. Suppose that k 1 and X/ p n satisfies µ-incoherence with µ . 1/kr. Take any precision ✏ > 0, and let the initialization be such that
0 < ↵ 
✓ ✏
p + 1
◆4/N ^ (w?max, w ? min, ✏, N) ^ (w ? min, N). (3)
For any iteration t that satisfies Tl(w ? , ↵, N, ⌘, ⇣, ✏)  t  Tu(w ? , ↵, N, ⌘, ⇣, ✏), (4)
where Tl(·) and Tu(·) are given in (25) of the Appendix, the gradient descent algorithm (2) with step size ⌘  ↵ N
8N2(w?max) (3N 2)/N yields the iterate wt with the following property:
|wt,i w ? i | .
8 <
:
1 nX T⇠ 1 _ ✏ if i 2 S and w ? min . 1 nX T⇠ 1 _ ✏, 1 n (X T⇠)i _ kµ 1 nX T⇠ 1S 1 _ ✏ if i 2 S and w ? min & 1 nX T⇠ 1 _ ✏,
↵ N/4 if i /2 S.
(5)
In the special case w? = 0, if ↵  ⇣
✏ p+1
⌘4/N , ⌘  1
N(N 1)⇣↵(N 2)/2 and t  Tu(w ? , ↵, N, ⌘, ⇣, ✏),
then we have |wt,i w?i |  ↵N/4, 8i.
Theorem 1 states the convergence of the gradient descent algorithm (2) in `1-norm. The exact formula of Tl(·) and Tu(·) is omitted here due to the space limitation. We ensure that Tu(·) > Tl(·) so that there indeed exists some epochs to early stop at. The error bound on the signal is invariant to the choice of N 2, and the overall bound generalizes that of [11] for N = 2. We also establish the convergence result in `2-norm in the following corollary: Corollary 1. Suppose the noise vector ⇠ has independent 2-sub-Gaussian entries and ✏ = 2 q
2 log(2p) n . Under the assumptions of Theorem 1, the gradient descent algorithm (2) would
produce iterate wt satisfying kwt w?k 2 2 . (k 2 log p)/n with probability at least 1 1/(8p3).
Note that the error bound we obtain is minimax-optimal, which is the same as [11] in the N = 2 case. However, with some calculation, the sample complexity we obtain here is n & k2r2, while the sample complexity in [11] is n & k2 log2 r log p/k. Although neither our work nor [11] achieved the optimal sample complexity k log p/k, the goal of this work is to understand how the depth parameter affects implicit sparse regularization.
Let us now discuss the implications of Theorem 1 and the role of initialization and early stopping:
(a) Requirement on initialization. To roughly understand the role of initialization and the effect of N , we look at the non-negative case where w? < 0 and w = uN . This simplifies our discussion while still capturing the essential insight of the general setting. At each step, the “update” on w can be translated from the corresponding gradient update of u = w1/N as
w0 = ↵ N1,
wt+1 = wt
✓ 1 2N⌘
n
✓ XTX(wt w ⇤) XT⇠ ◆ w(N 2)/Nt ◆N .
(6)
In order to guarantee the convergence, we require the initialization ↵ to be sufficiently small so the error outside the support can be controlled. On the other hand, too small initialization slows down the convergence of the signal. Interestingly, the choice of N affects the allowable initialization ↵ that results in guarantees on the entries inside and outside the support.
Specifically, the role of N is played by the term w(N 2)/Nt in (6), which simply disappears as N = 2. Since this term only affects the update of wt+1 entry-wise, we only look at a particular entry wt of wt. Let wt represent an entry outside the support. For N > 2, the term w (N 2)/N t is increasingly small as N increases and wt < 1. Therefore, with a small initialization, it remains true that wt < 1 for the early iterations. Intuitively, this suggests that the requirement on the upper bound of the initialization would become looser when N gets larger. This indeed aligns with the behavior of the upper bound we derive in our theoretical results. Since ↵ = w1/N0 increases naturally with N , we fix w0 = ↵N instead of ↵ to mimic the same initialization in terms of w0, for the following comparison.
We formalize this insight in Theorem 3 in Appendix A and show the convergence of (6) under the special, non-negative case. Note that, in terms of initialization requirement, the only difference from Theorem 1 is that we no longer require the term (w?min, N) in (3). Remark 1. We investigate how the depth N influences the requirement on initialization due to the change on gradient dynamics. We rewrite (w?max, w?min, ✏, N) in terms of w0 = ↵
N , and therefore the upper bound for w0 under the simplified setting of non-negative signals (Theorem 3) is
w0 
✓ 1
8
◆2N/(N 2) ^ (w?max) (N 2)/N
log w ? max ✏
!2N/(N 2) ^ (w?min) (N 2)/N
log w?min
✏
!4N/(N 2) .
We start by analyzing each term in the upper bound. First, we notice that 1 8 2N/(N 2) is increasing with respect to N . For the second term,
(w?max) (N 2)/N
log w ? max ✏
!2N/(N 2) =
(w?max) 2
(log w ? max ✏ )
2N/(N 2) ,
the denominator gets smaller as N increases when we pick the error tolerance parameter ✏ small. Therefore, we get that the second term is getting larger as N increases. The last term
✓ (w?min) (N 2)/N
log w?min
✏
◆4N/(N 2) follows a similar argument. We see that it is possible to pick a larger
initialization w0 = ↵N for larger N . We will demonstrate that below in our experiments.
(b) Early stopping. Early stopping is shown to be crucial, if not necessary, for implicit sparse regularization [11, 12]. Interestingly, [13, 18] studied the similar depth-N polynomial parametrization but did not realize the need of early stopping due to an oversimplification in the model. We will discuss this in details in Section 4.1. We are able to explicitly characterize the window of the number of iterations that are sufficient to guarantee the optimal result. In particular, we get a lower bound of the window size for early stopping to get a sense of how it changes with different N . Theorem 2 (Informal). Define the early stopping window size as Tu(w?, ↵, N, ⌘, ⇣, ✏) Tl(w?, ↵, N, ⌘, ⇣, ✏), the difference between the upper bound and lower bound of the number of iterations in (4) of Theorem 1. Fixing ↵ and ⌘ for all N , the early stopping window size is increasing with N under mild conditions.
We defer the formal argument and proof of Theorem 2 to Appendix D.3. We note that the window we obtain in Theorem 1 is not necessarily the largest window that allows the guarantee, and hence the early stopping window size can be effectively regarded a lower bound of that derived from the largest window. We note that a precise characterization of the largest window is difficult. Although we only show that this lower bound increases with N , we see that the conclusion matches empirically with the largest window. Fix the same initialization ↵ = 0.005 and step size ⌘ = 0.01 for N = 2, 3, 4, we show the coordinate path in Figure 1. We can see that as N increases, the early stopping window increases and the error bound captures the time point that needs stopping quite accurately. The experimental details and more experiments about early stopping is presented in Section 5.
Remark 2. Similar to Theorem 2, we look at how initialization scale affects the early stopping window for any fixed N > 2. With ⌘ fixed, the early stopping window is increasing as the initialization ↵ decreases.
We defer the detailed calculation to Section D.4. This generalizes the finding that vanishing initialization increases the gap between the phase transition times in [14] from N = 2 to any N > 2.

heading: 4 Proof Ingredients
text: The goal of this paper is to understand how generalization and gradient dynamics change with different N > 2. For N = 2, gradient descent yields both statistically and computationally optimal recovery under the RIP assumption [11]. The matrix formulation of the same type of parametrization is considered in the setting of low-rank matrix recovery, and exact recovery can be achieved in the noiseless setting [15, 16]. The key proof ingredient is to reduce the convergence analysis to one-dimensional iterates and differentiate the convergence on the support from the error outside the support. Before we get into that, we conduct a simplified gradient flow analysis.

heading: 4.1 A Simplified Analysis
text: Consider a simplified problem where the target signal w? is non-negative, n 1XTX = I and the noise is absent. We omit the reparametrization of vN like before and the gradient descent updates on u will be independent for each coordinate. The gradient flow dynamics of w = uN is derived as
@wi
@t =
@wi
@ui
@ui
@t =
@wi
@ui
@L @ui = 2N2(w?i wi)w 2 2N i , (7)
for all i 2 {1, 2, . . . , p}. Notice that wi increases monotonically and converges to w?i if w?i is positive or otherwise keeps decreasing and converges to 0 if w?i = 0. As such, we can easily distinguish the support and non-support. In fact, gradient flow with dynamics as in (7) would exhibit a behavior of “incremental learning” — the entries are learned separately, one at a time [13]. However, with the presence of noise and perturbation arising from correlated designs, the gradient flow may end up over-fitting the noise. Therefore, early stopping as well as the choice of step size is crucial for obtaining the desired solution [11]. We use (7) to obtain a gradient descent update:
wi,t+1 = wi,t(1 + 2N 2 ⌘(w?i wi,t)w 1 2N i,t ). (8)
The gradient descent with N = 2 is analyzed in [11]. However, when N > 2, the presence of w1 2 N i,t imposes an asymmetrical effect on the gradient dynamics. The difficulty of analyzing such gradient descent (8) is pointed out in [13]. More specifically, the recurrence relation is not solvable. However, gradient descent updates still share similar dynamics with the idealized gradient flow in (7). Inspired by this effect, we are able to show that the entries inside the support and those outside the support are learned separately with a practical optimization algorithm shown in (2) and (12). As a result, we are able to explore how the depth N affects the choice of step size and early stopping criterion.

heading: 4.2 Proof Sketch
text: Growth rate of gradient descent. We adopt the same decomposition as illustrated in [11], and define the following error sequences:
bt = 1
n XTXet
1 n XT⇠, pt =
✓ 1
n XTX I
◆ (st w ?) , (9)
where et and st stand for error and signal accordingly, and the definitions can be found in (13) in Appendix. We can then write the updates on st and et as
st+1 = st (1 2N⌘(st w ? + pt + bt) s (N 2)/N t ) N ,
et+1 = et (1 2N⌘(pt + bt) e (N 2)/N t ) N .
(10)
To illustrate the idea, we think of the one-dimensional updates {st}t 0 and {et}t 0, ignore the error perturbations pt and bt in the signal updates {st}t 0, and treat kpt + btk1  B in the error updates {et}t 0.
st+1 = st(1 2N⌘(st w ?)s(N 2)/Nt ) N , et+1 = et(1 2N⌘Be (N 2)/N t ) N . (11)
We use the continuous approximation to study the discrete updates. Therefore, we can borrow many insights from the analysis about gradient flow to overcome the difficulties caused by w1 2 N
i,t as pointed out in equation (8). With a proper choice of step size ⌘, the number of iterations Tl for st converging to w? is derived as
Tl 
Tl 1X
t=0
st+1 st
2N2⌘(w? st)s (2N 2)/N t
 1
N2⌘w?
Z w?
↵N
1
s(2N 2)/N ds + O
✓ w ? ↵ N
↵2N 2
◆ .
The number of iterations Tu for et staying below some threshold ↵N/4 is derived as
Tu
Tu 1X
t=0
et+1 et
4N2⌘Be(2N 2)/Nt
1
4N2⌘B
Z ↵N/4
↵N
1
e(2N 2)/N de.
With our choice of coherence µ in Theorem 1, we are able to control B to be small so that Tl is smaller than Tu. This means the entries on the support converge to the true signal while the entries outside the support stay around 0, and we are able to distinguish signals and errors.
Dealing with negative targets. We now illustrate the idea about how to generalize the result about non-negative signals to general signals. The exact gradient descent updates on u and v are given by:
ut+1 = ut
✓ 1 2N⌘ ✓ 1
n XT(X(wt w ?) ⇠) uN 2t
◆◆ ,
vt+1 = vt
✓ 1 + 2N⌘ ✓ 1
n XT(X(wt w ?) ⇠) vN 2t
◆◆ .
(12)
The basic idea is to show that when w?i is positive, v?i remains small up to the early stopping criterion, and when w?i is negative, u?i remains small up to the early stopping criterion. We turn to studying the gradient flow of such dynamics. Write r(t) = 1nX
T(X(w(t) w?) ⇠. It is easy to verify that the gradient flow has a solution:
u(t) = ✓ ↵ 2 N1 + 2N(N 2)⌘ Z t
0 r( )d
◆ 1 2 N
,
v(t) = ✓ ↵ 2 N1 2N(N 2)⌘ Z t
0 r( )d
◆ 1 2 N
.
We may observe some symmetry here, when ui,t is large, vi,t must be small. For the case wi > 0, to ensure the increasing of ui,t and decreasing of vi,t as we desire, the initialization needs to be smaller than wi, which leads to the extra constraint on initialization (w?min, ✏) with order of O(w?min) as defined before. It remains to build the connection between gradient flow and gradient descent, where again we uses the continuous approximation as before. The detailed derivation is presented in Appendix B.3.

heading: 5 Simulation Study
text: We conduct a series of simulation experiments1 to further illuminate our theoretical findings. Our simulation setup is described as follows. The entries of X are sampled as i.i.d. Rademacher random variables and the entries of the noise vector ⇠ are i.i.d. N(0, 2) random variables. We let w? = 1S . The values for the simulation parameters are: n = 500, p = 3000, k = 5, = 1, = 0.5 unless otherwise specified. For `2-plots each simulation is repeated 30 times, and the median `2 error is depicted. The shaded area indicates the region between 25th and 75th percentiles pointwisely.
Convergence results. We start by showing that the general choice of N leads to the sparse recovery, similar to N = 2 in [11], as shown in our main theorem. We choose different values of N to illustrate the convergence of the algorithm. The result on simulated data is shown in Figure 2, and we defer the result on MNIST to Appendix E. Note that the ranges in the x-axes of these figures differ due to different choice of N and ⌘. We observe that as N increases, the number of iterations increases significantly. This is due to the term uN 2 and vN 2 in (12), and the step size ⌘ ⇡ 1N2 . With a very small initialization, it takes a large number of iterations to escape from the small region (close to 0).
Larger initialization. As discussed in Remark 1, the upper bound on initialization gets larger with larger N . We intentionally pick a relatively large ↵N = 2 ⇥ 10 3 where the algorithm fails to converge for N = 2. With the same initialization, the recovery manifests as N increases (Figure 3).
Early stopping window size. Apart from the coordinate path shown in Figure 1, we obtain multiple runs and plot the log-`2 error (the logarithm of the `2-error) of the recovered signals to further confirm the increase of early stopping window, as shown in Section 3. Note that for both Figures 1 and 4, we set n = 100 and p = 200. Since ↵N would decrease quickly with N , which would cause the algorithm takes a large number of iterations to escape from the small region. We fix ↵N = 10 5 instead of fixing ↵ for Figure 4.
Incremental learning dynamics. The dynamics of incremental learning for different N is discussed in [13]. The distinct phases of learning are also observed in sparse recovery (Figure 5), though we do
1The code is available on https://github.com/jiangyuan2li/Implicit-Sparse-Regularization.
not provide a theoretical justification. Larger values of N would lead to more distinct learning phases for entries with different magnitudes under the same initialization ↵N and step size ⌘.
Kernel regime. As pointed out in [18], the scale of initialization determines whether the gradient dynamics obey the “kernel” or “rich” regimes for diagonal linear networks. We have carefully analyzed and demonstrated the sparse recovery problem with small initialization, which corresponds to the “rich” regime. To explore the "kernel" regime in a more practical setting, we set n = 500, p = 100, and the entries of w? are i.i.d. N (0, 1) random variables. The noise level is = 25, and the initialization and step size is set as ↵N = 1000 and ⌘ = 10 7 for all N . Note that we are not working
in the case n ⌧ p as [18]. We still observe that the gradient dynamics with large initialization (Figure 6) can be connected to ridge regression if early stopping is deployed.

heading: 6 Conclusions and Future Work
text: In this paper, we extend the implicit regularization results in [11] from N = 2 to general N > 2, and further study how gradient dynamics and early stopping is affected by different choice N . We show that the error bound is invariant with different choice of N and yields the minimax optimal rate. The step size is of order O(1/N2). The initialization and early stopping window gets larger when increasing N due to the changes on gradient dynamics. The incremental learning dynamics and kernel regime of such parametrizations are empirically shown, however not theoretically justified, which is left for future work.
The convergence result can be further improved by relaxing the requirement on the incoherence of design matrix from µ . 1kw?max/w?min to µ . 1 k log(w?max/w ? min)
, similar to [11]. Overall, we believe that such an analysis and associated techniques could be applied for studying other, deeper nonlinear models in more practical settings.

heading: Acknowledgements
text: This work was supported in part by the National Science Foundation under grants CCF-1934904, CCF-1815101, CCF-2005804, and DMS-1711952.

Review summaries:
In this work, the authors study the implicit bias of gradient descent for sparse linear regression when parameterizing the ground truth signal as a depth-$N$ diagonal linear network. Previous works have shown that when utilizing a Hadamard parameterization $w = u \odot u - v \odot v$, gradient descent exhibits an implicit bias towards sparse solutions, without any explicit regularization encouraging sparsity. Follow-up work has also considered the depth-$N$ case in simplified settings (e.g., without noise, positive signals, and idealized measurement setups). The authors extend this analysis to the general depth-$N$ case $w = u^N - v^N$ in the presence of noise. They prove that when the measurement matrix is $\mu$-incoherent with sufficiently small $\mu$ and noise is present, there is an early stopping window depending on the problem parameters where any iterate during this time is close in $\ell_{\infty}$ norm to the underlying ground truth sparse signal. These results hold for any depth $N \geqslant 2$ and showcase that depth allows for larger initialization sizes and that the early stopping window grows with deeper parameterizations. Experimental evidence on toy examples is also provided corroborating several elements of the proposed theory.

The paper under review studies the implicit sparse regularization induced by gradient descent on a diagonal linear network-like of arbitrary depth. Doing so, they extend a past work on depth 2 showing the possible benefit of larger depth. More precisely the authors show the role of early stopping to recover a sparse signal.

In this paper the authors study implicit sparse regularization for diagonal linear networks. Specifically, the authors extend the previous result of Vaskevicius et al. [9] to the case of depth>2, showing that for noisy regression under a relaxed assumption on the inputs, increasing depth enlarges the scale of working initialization and the early-stopping window. 

This paper studies the implicit bias of gradient descent on the depth-N diagonal linear network for the problem of sparse recovery. It shows the convergence results for gradient descent and provides error bounds to the ground truth sparse signal. The result highlights the importance of early stopping in preventing overfitting and shows the role of network depth in extending the stopping window.