High-Dimensional Linear Regression via Implicit Regularization
Peng Zhao, Yun Yang, Qiao-Chu He
TL;DR
This paper demonstrates that overparameterized gradient descent with early stopping implicitly regularizes high-dimensional linear regression, achieving near-optimal sparse solutions without explicit penalties.
Contribution
It introduces a new class of estimators defined through discretized gradient dynamics that implicitly regularize, outperforming explicit regularization methods under certain conditions.
Findings
Implicit regularization via gradient descent leads to sparse, near-optimal solutions.
Early stopping ensures convergence to low-bias estimators.
Simulations show advantages over explicit regularization methods.
Abstract
Many statistical estimators for high-dimensional linear regression are M-estimators, formed through minimizing a data-dependent square loss function plus a regularizer. This work considers a new class of estimators implicitly defined through a discretized gradient dynamic system under overparameterization. We show that under suitable restricted isometry conditions, overparameterization leads to implicit regularization: if we directly apply gradient descent to the residual sum of squares with sufficiently small initial values, then under some proper early stopping rule, the iterates converge to a nearly sparse rate-optimal solution that improves over explicitly regularized approaches. In particular, the resulting estimator does not suffer from extra bias due to explicit penalties, and can achieve the parametric root-n rate when the signal-to-noise ratio is sufficiently high. We alsoâŠ
| S1 | S2 | S3 | S4 | S5 | S6 | S7 | S8 | ||
| Strong Signals | GD | 0.520 | 0.448 | 0.510 | 0.568 | 0.385 | 0.290 | 0.465 | 0.460 |
| (0.0428) | (0.0530) | (0.0607) | (0.0850) | (0.0533) | (0.0465) | (0.0858) | (0.0863) | ||
| lasso | 3.11 | 3.10 | 3.42 | 4.51 | 4.62 | 4.04 | 4.40 | 6.98 | |
| (0.219) | (0.173) | (0.242) | (0.274) | (0.279) | (0.205) | (0.306) | (0.452) | ||
| SCAD | 0.613 | 0.533 | 0.650 | 0.691 | 0.519 | 0.401 | 0.595 | 0.646 | |
| (0.0464) | (0.0679) | (0.0702) | (0.103) | (0.0527) | (0.0574) | (0.0837) | (0.0776) | ||
| MCP | 0.628 | 0.552 | 0.594 | 0.733 | 0.484 | 0.405 | 0.595 | 0.708 | |
| (0.0392) | (0.0779) | (0.0809) | (0.0902) | (0.0706) | (0.0597) | (0.0741) | (0.0680) | ||
| Weak Signals | GD | 0.996 | 1.226 | 0.586 | 0.651 | 0.699 | 0.474 | 0.556 | 0.368 |
| (0.259) | (0.255) | (0.0938) | (0.122) | (0.117) | (0.0609) | (0.127) | (0.0351) | ||
| lasso | 2.530 | 3.063 | 2.292 | 2.620 | 2.048 | 1.952 | 1.330 | 1.178 | |
| (0.195) | (0.151) | (0.119) | (0.168) | (0.208) | (0.121) | (0.100) | (0.115) | ||
| SCAD | 3.942 | 3.781 | 3.331 | 3.204 | 3.241 | 2.369 | 3.881 | 3.397 | |
| (0.409) | (0.137) | (0.197) | (0.192) | (0.219) | (0.150) | (0.270) | (0.313) | ||
| MCP | 2.629 | 2.153 | 2.248 | 1.730 | 2.878 | 1.570 | 6.149 | 4.984 | |
| (0.279) | (0.180) | (0.261) | (0.135) | (0.346) | (0.174) | (0.312) | (0.447) |
| S1 | S2 | S3 | S4 | S5 | S6 | S7 | S8 | ||
| Strong Signals | GD | 0.638 | 0.636 | 0.640 | 0.634 | 0.646 | 0.651 | 0.641 | 0.642 |
| (0.0597) | (0.0753) | (0.0718) | (0.0709) | (0.0491) | (0.0510) | (0.0406) | (0.0498) | ||
| lasso | 0.676 | 0.672 | 0.671 | 0.685 | 0.693 | 0.699 | 0.696 | 0.708 | |
| (0.0899) | (0.0981) | (0.0932) | (0.0510) | (0.0568) | (0.0918) | (0.0615) | (0.0488) | ||
| SCAD | 0.638 | 0.638 | 0.637 | 0.637 | 0.650 | 0.654 | 0.643 | 0.647 | |
| (0.0530) | (0.0724) | (0.0716) | (0.0713) | (0.0470) | (0.0451) | (0.0402) | (0.0434) | ||
| MCP | 0.637 | 0.639 | 0.638 | 0.637 | 0.650 | 0.652 | 0.644 | 0.647 | |
| (0.0516) | (0.0756) | (0.0684) | (0.0753) | (0.0475) | (0.0441) | (0.0435) | (0.0408) | ||
| Weak Signals | GD | 1.023 | 1.039 | 1.007 | 1.038 | 1.013 | 1.030 | 1.000 | 1.013 |
| (0.114) | (0.128) | (0.141) | (0.111) | (0.132) | (0.089) | (0.151) | (0.069) | ||
| lasso | 1.061 | 1.093 | 1.050 | 1.097 | 1.055 | 1.070 | 1.034 | 1.054 | |
| (0.0908) | (0.0841) | (0.138) | (0.0924) | (0.0921) | (0.0873) | (0.132) | (0.106) | ||
| SCAD | 1.082 | 1.108 | 1.066 | 1.110 | 1.074 | 1.091 | 1.070 | 1.094 | |
| (0.107) | (0.104) | (0.197) | (0.0868) | (0.0866) | (0.0914) | (0.141) | (0.110) | ||
| MCP | 1.057 | 1.067 | 1.038 | 1.071 | 1.055 | 1.063 | 1.066 | 1.092 | |
| (0.0953) | (0.0810) | (0.162) | (0.112) | (0.115) | (0.0961) | (0.139) | (0.0640) |
| 1e-10 | 1e-5 | 1e-3 | 0.1 | 1 | 10 | |
|---|---|---|---|---|---|---|
| 7.433e-13 | 5.703e-7 | 1.289e-4 | 2.884e-2 | 2.987e-1 | 8.823e-1 | |
| 1.492e-13 | 1.141e-7 | 2.577e-5 | 5.769e-3 | 5.974e-2 | 1.765e-1 | |
| 1.492e-13 | 1.141e-7 | 2.577e-5 | 5.769e-3 | 5.974e-2 | 1.765e-1 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
MethodsEarly Stopping · Linear Regression
High-Dimensional Linear Regression via Implicit Regularization
Peng Zhao, Department of Statistics, Texas A&M University
ââ
Yun Yang, Department of Statistics, University of Illinois Urbana-Champaign
ââ
Qiao-Chu He, Southern University of Science and Technology
Abstract
Many statistical estimators for high-dimensional linear regression are -estimators, formed through minimizing a data-dependent square loss function plus a regularizer. This work considers a new class of estimators implicitly defined through a discretized gradient dynamic system under over-parametrization. We show that under suitable Restricted Isometry conditions, over-parametrization leads to implicit regularization: if we directly apply gradient descent to the residual sum of squares with sufficiently small initial values, then under some proper early stopping rule, the iterates converge to a nearly sparse rate-optimal solution that improves over explicitly regularized approaches. In particular, the resulting estimator does not suffer from extra bias due to explicit penalties, and can achieve the parametric root- rate when the signal-to-noise ratio is sufficiently high. We also perform simulations to compare our methods with high dimensional linear regression with explicit regularization. Our results illustrate the advantages of using implicit regularization via gradient descent after over-parametrization in sparse vector estimation.
1 Introduction
In high dimensional linear regression with noise , where each component of noise follows an independent sub-Gaussian distribution with scale bounded by , the goal is to parsimoniously predict the response as a linear combination of numerous covariates , and conduct statistical inference on the coefficients (Tibshirani,, 1996; Donoho,, 2006). By leveraging certain lower dimensional structure in the regression coefficient such as a sparsity constraint , where counts the number of nonzero elements in , the number of covariates may be substantially larger than the sample size . Because of the intrinsic computational hardness in dealing with the metric, people instead use different metrics as surrogates, and cast the estimation problem into various convex or non-convex optimization problems. Many approaches have been proposed for high dimensional regression by solving certain penalized optimization problem, including basis pursuit (Chen et al.,, 2001), the lasso (Tibshirani,, 1996), the Dantzig selector (Candes and Tao,, 2007), smoothly clipped absolute deviation (Fan and Li,, 2001), minimax concave penalty (Zhang,, 2010) and so on. In this work, we focus on the estimation of without explicitly specifying a penalty. Recent work (Hoff,, 2017) showed that through a change-of-variable via Hadamard product over-parametrization, the non-smooth convex optimization problem for the lasso:
[TABLE]
can be reformulated as a smoothed optimization problem at a cost of introducing non-convexity. Due to the smoothness feature, simple and low-cost first-order optimization methods such as gradient descent and coordinate descent can recover . Despite the non-convexity and exponentially many stationary points induced by the change-of-variable, these first-order algorithms exhibit encouraging empirical performance (Hoff,, 2017).
In this work, we consider the same Hadamard product over-parametrization as in Hoff, (2017), where and denotes the Hadamard product operator, which performs element-wise product between two vectors. Instead of solving the penalized optimization problem (1), we consider directly applying the gradient descent to the quadratic loss function
[TABLE]
In the special case of noiseless responses when , minimizing jointly over is a highly non-convex optimization problem with exponentially many saddle points. Interestingly, our theory indicates that by initializing and arbitrarily close to zero, a properly tuned gradient method converges to least -norm solution when the design matrix satisfies the prominent Restricted Isometry Property (Candes,, 2008), thereby inducing algorithmically implicit regularization. In the general case when , we show that by combining gradient method with early stopping (Zhang and Yu,, 2005; Raskutti et al.,, 2014), the resulting estimator overcomes the large bias (Vito et al.,, 2005; Yao et al.,, 2007) suffered by usual explicitly penalized approaches and leads to more accurate estimation. In particular, we show that by iteratively updating for a certain number steps, the resulting estimator can adapt to an optimal convergence rate of when all signals are relatively strong. When both strong signals and weak signals exist, our estimator attains the rate , with denoting the number of strong signals and weak signals respectively. Our numerical studies also illustrate encouraging improvements over existing explicitly regularization-based approaches.
Our work complements the recent surge of literature on implicit regularization in first-order iterative methods for solving non-convex optimization in machine learning (Gunasekar et al.,, 2017, 2018; Soudry et al.,, 2018; Li et al.,, 2018) that exclusively focused on noiseless or perfectly separable data. We show that in the context of sparse vector estimation, through a simple change-of-variable, a simple gradient descent initialized near-zero induces sparsity, and solves the non-convex optimization problem with provable guarantees. More broadly, this paper presents a new way of designing statistical estimatorsâunlike usual estimators such as -estimators and -estimators that are defined either as optimizers or roots of some criterion functions or system of equations, the estimator considered in this paper is implicitly defined through a discretized gradient dynamic system. Such an algorithmically induced estimator has several advantages over existing estimation procedures for high dimensional linear regression under sparsity constraints. First, the proposed implicit regularized estimator suffers less bias than some penalty-induced estimators, leading to more accurate estimation. Also, despite the non-convexity nature, our method has a natural initialization that provably leads to the optimal solution. In comparison, state-of-the-art -estimators based on non-convex penalties require stringent conditions on their initializations to avoid bad local minima with bad estimation errors.
2 Background and Our Method
2.1 Setup and notations
Recall that is the unknown -sparse signal in to be recovered. Let denote the index set that corresponds to the nonzero components of , and the size of is then . For two vectors , we call as their Hadamard product, whose components are for . For two vectors , we use the notation to indicate element-wise âgreat than or equal toâ. When there is no ambiguity, we use to denote the self-Hadamard product of . For a function , , we use and to denote its partial derivative relative to and , respectively. For any index set and vector , we use to denote the sub-vector of formed by concatenating the components indexed by . Let denote the vector with all entries as , and as the identity matrix in . Let be the diagonal matrix with one on the th diagonal for and [math] elsewhere. For a vector , we use to denote its vector--norm, and its -norm. Let denote the uniform distribution over interval . For a symmetric matrix , let denote its smallest eigenvalue. For two sequences and , we use the notation or to mean there exist some constants and independent of such that or for all , respectively, and to mean and .
2.2 Related literature
Li et al., (2018) studies the theory for implicit regularization in matrix sensing, which requires the data to be perfectly measured and has different geometric structures as linear regression. Hoff, (2017) considers the Hadamard product parametrization to optimize the parameters in high-dimensional linear regression. In particular, their objective function involves an penalty on that is equivalent to the penalty on , and the solution is precisely the lasso solution. During the revision of our paper, we learned that a recent work (Vaskevicius et al.,, 2019) also studied high-dimensional linear regression via implicit regularization via a slightly different parameterization. Our work is different from Vaskevicius et al., (2019) in many aspects. A detailed comparison between the two works is provided in Section F of the supplementary material.
2.3 Gradient descent with Hadamard product parametrization
As we mentioned in the introduction, we consider augmenting the -dimensional vector into two -dimensional vectors through . Instead of solving the lasso problem with replaced with , we consider directly applying gradient descent to the quadratic loss function . In particular, we apply the updating formula , , with random initial values , for (notice that is a saddle point of the objective function, so we need to apply a small perturbation on the initial values). This leads to the standard gradient descent algorithm, and the iterates tend to converge to a stationary point of that satisfies the first order optimality condition and . However, stationary points of can be local minimum, local maximum, or saddle points (when the Hessian matrix contains both positive and negative eigenvalues). The following result provides the optimization landscape of , showing that does not have local maximums, all its local minimums are global minimum, and all saddle points are strict. The strict saddle points are saddle points with the negative smallest eigenvalues for the Hessian matrix.
Lemma 2.1**.**
* does not have local maximum, and all its local minimums are global minimum. In particular, is a global minimum of if and only if*
[TABLE]
In addition, any saddle point of is a strict saddle, that is, \lambda_{\min}\big{(}\nabla^{2}_{g,l}f(g^{\dagger},l^{\dagger})\big{)}<0.
According to the first order condition associated with
[TABLE]
there could be exponentially many saddle points as a solution to this equation, for example, for those satisfying
[TABLE]
for any non-empty subset of . Consequently, the gradient descent algorithm may converge to any of these bad saddle points. To see this, if we initialize in a way such that the components in the index set are zero, or , then these components will remain zero forever in the gradient iterations. Fortunately, the following result implies that as long as we use a random initialization for with continuous probability density functions over as in the gradient descent algorithm, then the gradient descent almost surely converges to a global minimum.
Lemma 2.2**.**
Suppose the step size is sufficiently small. Then with probability one, the above gradient descent algorithm converges to a global minimum of .
In the low-dimensional regime where the design matrix has full column rank, the solution to the normal equation is unique, which is also the least-squares estimator. Under this scenario, Lemma 2.1 and Lemma 2.2 together certify that the above gradient descent algorithm will converge to this optimal least-squares estimator. However, in the high-dimensional regime, which is the main focus of the paper, the normal equation has infinitely many solutions, and it is not clear which solution the algorithm tends to converge to. For example, if we consider instead applying the gradient descent to the original parameter in the objective function with initialization , then the iterates will converge to the minimal -norm solution of the normal equation. Interestingly, as we will illustrate in the following, under the Hadamard parametrization, the gradient descent algorithm now tends to converge to the minimal -norm solution under certain conditions on initialization and the design matrix, thereby inducing sparsity implicitly.
2.4 Gradient descent converges to sparse solution
In this subsection, we provide two different perspectives for understanding the following informal statement on the behavior of simple gradient descent for the loss function defined in (2) under the Hadamard product parameterization . We assume that the noise variance throughout this subsection and will discuss the case later.
Informal Statement: If we initialize the algorithm to be arbitrarily close to , then under suitable conditions on the design , a simple gradient descent converges to a solution of basis pursuit problem:
[TABLE]
Consider the under-determined system , where has full row rank. Our first intuition comes from the fact that a zero-initialized gradient descent algorithm over for solving finds a minimal -norm solution to . Then under the Hadamard product parameterization , the fact that gradient descent tends to find the minimal -norm solution suggests (this is not rigorous) that the gradient descent algorithm for jointly minimizing over tends to converge to a solution to with a minimal -norm . However, a minimal -norm solution to must satisfy for each (otherwise we can always construct another solution with strictly smaller -norm), which implies . As a consequence, should be the minimal -norm solution to .
Another way to understand the difference in the evolutions of gradient descents for and is by noticing that the gradient \nabla_{g_{j}}f(g,l)=l_{j}\cdot\nabla_{\beta_{j}}h(\beta)\big{|}_{\beta=g\circ l} in the new parametrization, for each , has an extra multiplicative factor of than the gradient in the usual least squares of minimizing . It is precisely this extra multiplicative factor that helps select important signals, which are the nonzero regression coefficients and prevent unimportant signals, which are the zero regression coefficients, to grow too fast at the early stage of the evolution when both and are close to zero. Our second perspective comes from considering the limiting gradient dynamical system of the problem with an infinitesimally small step size , which is studied by Gunasekar et al., (2017) and Gunasekar et al., (2018). In particular, the behavior of this limiting dynamical system is captured by the ordinary differential equations
[TABLE]
where r(t)=n^{-1}\big{[}X\{g(t)\circ l(t)\}-y\big{]}\in\mathbb{R}^{p}, and for simplicity we fixed the initialization. As shown in the Section A.1 of the supplementary material, the limiting point condition coincides with the Karush-Kuhn-Tucker condition for the basis pursuit problem (3). Moreover, under the same solution continuity assumption, if there is additive noise on , then the limiting point condition matches the Karush-Kuhn-Tucker condition for the lasso (1). Again, the above two arguments are based on the hard-to-check solution continuity assumption. In the next section, we provide a formal proof of the result without making this assumption, albeit under a somewhat strong Restricted Isometry Property on .
Finally, as suggested by one reviewer, our current framework can be extended to allow plus type elastic net regularization by adding a penalty term to the current sum of squares objective function. Details about its formulation can be found in Subsection A.4 in the supplementary material.
2.5 Gradient descent with early stopping
In this subsection, we consider the general case where the response contains noise: . In particular, we combine early stopping, a widely used implicit regularization technique (Zhang and Yu,, 2005; Raskutti et al.,, 2014), with the gradient descent to prevent overfitting. Algorithm 1 below provides the detailed implementation. Here in the algorithm, we adopted the hold-out method where the algorithm terminates once the testing error on the validation data started to increase, as the testing error is expected to display a U-shaped curve similar to the estimation error shown in Fig. 1.
Due to the connection of our method with the basis pursuit problem (3), one may naturally think that our method in the noisy case should be equivalent to a basis pursuit denoising problem:
[TABLE]
with some error tolerance level depending on the stopping criterion, which is equivalent to the lasso. Surprisingly, a simulation example below shows that the iterate path of the gradient descent Algorithm 1 contains estimates with much smaller errors than the lasso. Precisely, we adopt the simulation setting S2 in Section 4 . As comparisons, we also report the lasso solution path as a function of the regularization parameter and norm level. For our gradient descent algorithm, we set in the random initialization. From Figure 1, when the iteration number is around , even though the prediction errors in panel (a) of our algorithm and the optimal lasso solution in panel (b) are almost the same, the estimation error in panel (a) of our method is significantly lower than that of the lasso, illustrating the occurrence of the large bias of the lasso. Moreover, the stabilized region of our method in panel (a) corresponding to the optimal stopping time is relatively wide. Therefore, the performance of implicit regularization with early stopping tends to be robust to the stopping criterion. In addition, from panel (a) in Fig. 1, we can see that as the number of iterations increases, the norm of the iterate also increases. When the logarithm of the iteration number is within , the norm of the estimated coefficients tends to be stabilized around the norm of the true , corresponding to the most accurate estimation region in panel (a) of Fig. 1. In contrast, as we can see from panels  (b) and (c) of Fig. 1, the estimation error is very sensitive in the regularization parameter/ norm domain â the region corresponds to the smallest estimation accuracy is very narrow. A small change in the norm or regularization parameters in the solution leads to a drastic deterioration in the estimation accuracy. The above simulation results illustrate that our implicit regularization method is unaffected, or at least less suffered from the large bias.
We may also employ implicit regularization as a refining technique to improve upon an estimator based on explicit regularization, such as the lasso, by transforming the data so that our algorithm is âinitializedâ near . Algorithm 2 below summarizes the procedure. Numerical comparisons with Algorithm 1 is reported in subsection E of the supplementary material. In situations where signal strengths in are highly-unbalanced, meaning that the ratio between the largest and the smallest âs is large, Algorithm 1 may need iterations to converge according to Theorem 3.1. In comparison, Algorithm 2 only needs iterations in the worst case if the initial estimator satisfies \|\beta^{\dagger}-\beta^{*}\|_{\infty}=O\big{(}\sqrt{s\log p/n}\big{)}, which can be achieved, for example, by using the lasso, where recall that is the sparsity level of the true signal.
2.6 Early Stopping criterion
This subsection discusses several commonly used early stopping rules by treating the iteration number as a tuning parameter.
Hold-out or cross-validation: The simplest method is to use hold-out data as validation. For example, we can first randomly split the entire data into training data with size and validation data with size ; then run gradient descent on the training data while evaluating the prediction risk on the validation data for all . Since our goal is to minimize the estimation error, whose pattern is consistently captured by , we recommend deciding the final iteration number by:
[TABLE]
To make use of the entire dataset, we may perform cross-validation: first split data into fold; then use one fold as the validation set for evaluating the prediction risk and apply the gradient update to the rest folds; finally, the overall prediction risk at time is set to be the average prediction error across all choices of the validation set. Based on , criterion (6) can be used to obtain the estimated iteration number . Finally, we apply this from cross-validation to early stop the gradient algorithm with the entire dataset.
Steinâs unbiased risk estimate (SURE): Stein, (1981) suggested the use of the effective degrees of freedom as a surrogate to the population level prediction risk for selecting a best tuning parameter, which is the stopping time in our case. Under our settings, ignoring some higher-order discretization error term of order with being the step size, the updating formula of the prediction residual \tau_{t}=\big{\{}X(g_{t}\circ l_{t})-y\big{\}}\in\mathbb{R}^{n} induced by the gradient algorithm can be approximated by:
[TABLE]
Here, diag denotes the diagonal matrix whose diagonals are components of a vector . Let denote the accumulated multiplicative operator up to time , so that we have the approximation with . Following Zou et al., (2007), we define the (approximated) effective degrees of freedom (df) at time as . According to Efron, (2004), an -type statistic for approximating the prediction risk can be related to the effective df via
[TABLE]
Finally, the stopping time can be selected by criterion (6) with being replaced by . In practice, we may use any consistent plug-in estimator , such as the averaged residual sum of squares, to replace the unknown variance in . According to our simulation studies in Section A.3 of the supplementary material, early stopping based on has a very similar performance as the hold-out or cross-validation method. As a default method, we will use cross validation to select the iteration number in the rest of the paper to avoid estimating .
3 Theoretical Analysis
3.1 Assumptions
Recall that the true data generating model is where is a vector of independent sub-Gaussian random variables with scales bounded by , and the true parameter is -sparse. Within the nonzero signal components of , we define the index set of strong signals as and weak signals as for some constants , where , . According to the information-theoretic limits from Wainwright, (2009), weak signals of order in sparse linear regression are generally impossible to be jointly recovered or selected, but they can be detected in terms of the type I/II error control in hypothesis testings (Jin and Ke,, 2016). Therefore, our primary focus would be the estimation of strong signals. We use the notation to denote the -th largest absolute component value of , and let , which reflects the minimal strength for strong signals. We also use to denote the strong signal-condition number as the ratio between the largest absolute signal value to the smallest strong signal. We will also make use of the notation of Restricted Isometry Property (Candes, (2008)), which is a commonly used assumption (Candes and Tao,, 2007) in the high dimensional linear regression literature.
Definition 1** (Restricted Isometry Property).**
A matrix is said to satisfy the -Restricted Isometry Property if for any -sparse vector in , we have:
[TABLE]
As an easy consequence, if matrix satisfies -Restricted Isometry Property, then Euclidean inner-product between some sparse vectors is also approximately preserved, that is, \big{|}n^{-1}\,\langle Xu,Xv\rangle-\langle u,\,v\rangle\big{|}\leq\delta\,\|u\|\cdot\|v\| holds for any two -sparse vectors .
With these preparations, we make the following assumptions on the true data generating model, noise , the true parameter , design matrix , initialization parameter and step size in Algorithm 1.
Assumption 1**.**
The true data generating model is where is a vector of independent sub-Gaussian random variables with scales bounded by . For simplicity, we also assume that all norm of column vectors of are normalized to .
Assumption 2**.**
The true parameter is -sparse, and , that is, each nonzero signal in is either weak or strong. In addition, for some constant .
Assumption 3**.**
The design matrix satisfies -Restricted Isometry Property with .
Assumption 4**.**
The initial values for gradient descent are , , where the initialization parameter satisfies , and the step size satisfies .
Although our proof heavily relies on the Restricted Isometry Property in Definition 1, the numerical results in Section D of the supplementary materials provide some strong evidence suggesting that the Restricted Isometry Property is not necessary for the method to work. We conjecture that our theoretical conclusion in the following subsection remains valid as long as the much weaker restricted eigenvalue condition (Wainwright,, 2009; Bickel et al.,, 2009) is satisfied. More discussions about the consequences due to Assumptions 3 and 4 are provided in the next subsection.
3.2 Finite sample error bound
Our main theoretical result is summarized in the following theorem.
Theorem 3.1** (Finite sample error under early stopping).**
Suppose that Assumptions hold. Then there exist positive constants such that for any , it holds with probability at least that for every time with , the time -iterate satisfies
[TABLE]
Several remarks are in order. First, the error bound in the theorem suggests that in high-dimensional linear regression, combining early stopping with implicit regularization can significantly improve the estimation accuracy compared to explicitly regularized estimators. For example, when all signals are strong ( and ), the estimate attains a parametric rate of convergence that is independent of the dimension . In the more general situation where weak signals exist, the overall rate becomes dependent on the number of weak and strong signals, which is still nearly minimax-optimal (Zhao et al.,, 2018) in this regime.
Second, according to Theorem 3.1, the conditions and the final estimation error bound are not sensitive to the tuning parameters and , as long as they are chosen to be small enough. As a default choice for initialization parameter , we recommend using . For the step size , we may first select an based on minimizing the hold-out estimate of the prediction risk described in Section 2.6, and then set the step size for the entire training data by rescaling as since the choice of according to Theorem 3.1 satisfies , where is the training data size and the validation data size in the hold-out. In addition, the range of iteration number to achieve the desired estimation accuracy is wide, making the empirical performance less sensitive to the choice of the early stopping criterion.
Third, the proof of Theorem 3.1 is based on a similar strategy as in Li et al., (2018) for analyzing matrix factorized gradient descent in noiseless matrix sensing by dividing the gradient dynamic into two stages, where in the first stage, all significant signals stand out by quickly moving away from zero, as shown in Proposition 1 in Section C.3 in the supplement. In the second stage, all insignificant signals remain bounded and small while significant ones exponentially converge to their true values modulo statistical errors, as proved in Proposition 2 in Section C.3 in the supplement. To facilitate the analysis, in our proofs we adopt another reparametrization of and so that to keep track of the change of signs in components of parameter . Interestingly, after the first version of our work was posted on arXiv, Vaskevicius et al., (2019) proposed a similar implicit regularization algorithm for linear regression directly based on this parametrization. They presented a significantly improved analysis by proving an estimation error bound under weaker conditions on restricted isometry constant, step size, and strength of signals, and condition number than ours. However, their error bound is slightly worse than ours with an extra factor in front of the term corresponding to strong signals.
As a final remark, Theorem 3.1 also reflects some algorithmic benefit of Algorithm 2 based on a consistent initial estimator over Algorithm 1. Specifically, consider the largest signal conditional number regime where , which corresponds to the worse case scenario where the minimal signal strength is of order , i.e. the ratio between the largest component in absolute value in and the smallest component is large. In this extreme setting, Assumption 4 requires . As a consequence, Theorem 3.1 suggests that at least iterations is needed to guarantee Algorithm 1 to achieve the desired error bound. In comparison, suppose the initial estimator in Algorithm 2 satisfies some error bound condition as , then after the pre-processing step 2, the new conditional number can always be upper bounded by , where the worst case occurs when , so that now only iterations is needed for Algorithm 1. For example, if , which is satisfied by estimators from the lasso, then the iteration number can be reduced from to . This iteration complexity leads to a total computational complexity of , which has the same complexity scaling in terms of as the improved algorithm proposed in Vaskevicius et al., (2019).
4 Simulations and Real Data Analysis
We further demonstrate the advantages of our algorithm by considering simulation settings S1-S8 with , rows and the sparsity level . For S1-S4, we set , , with (here we denote ), while S5-S6 corresponds to S1-S4 with but everything else the same. For the case of strong signals we set signals as and noise level with ; while for weak-signals case, we set and the four signals are all . We generate observations independently and split into even parts, then use the first part for training, the second part for validation and the final part for testing. The evaluation metric is standardized estimation error and mean prediction error for the test data set. We compare the median of the standardized estimation errors and prediction errors with the lasso, smoothly clipped absolute deviation and minimax concave penalty by repeating times. We implement the lasso using glmnet R package (Friedman et al.,, 2010) while for smoothly clipped absolute deviation and minimax concave penalty, we use the R package ncvreg (Breheny and Huang,, 2011). The standard error of medians are calculated by bootstrapping the calculated errors times. For our algorithm, we use the initialization . For all other methods, the hyperparameter is chosen based on the prediction performance on the validation data from evenly spaced grids within . The simulation results in Table 1 below and Table 2 in Section A.6 in the supplement indicate that our methods consistently have the best performance over all explicit penalization-based competitors across all settings.
Now we consider variable selection when there exist some weak signals. Suppose the simulation settings are similar to S3 above, with only the true signals change. Let , and the strength of first signals is , while the other are , where . The first signals are too weak to be selected by all methods. However, since all methods are based on minimizing the prediction error, the effect of these weak signals pertains, and may increase the false discovery rate. Under the above settings, we perform a model selection based on minimized prediction errors through -fold cross-validation. For our method, we use the same regularization parameter as the lasso to perform a hard threshold after estimation. We repeat the process times and compare variable selection errors. From panel (a) in Figure 2, we can see our method are robust to the enhancement of false detection due to failure on detecting weak signals: although the true negative error of our method is , which means all weak signals can not be detected, the false detections of our methods are closed to zero. Minimax concave penalty performs similarly with our method in true negative but slightly worse in false positive, while for the lasso and smoothly clipped absolute deviation, although sometimes weak signals can be detected, the risk of false detections is high. Overall, our methods perform a consistent variable selection for strong signals and achieve better estimation than the competitors.
We compare our method with others to analyze the Riboflavin data set (BĂŒhlmann et al.,, 2014), which is available in hdi R package. The dataset contains observations of log-transformed riboflavin production rate versus the logarithm of the expression level of 4088 genes. Before estimation, we first perform independence screening (Fan and Lv,, 2008) based on the rank of the correlation strength for each predictor versus response to decrease the dimension of feature space into . Then we normalize and add the intercept column into the design matrix. For evaluation, we split the observations into training samples and testing samples, with performing -fold cross validation to select iteration steps and regularization parameters in the training data. Still, for our algorithm, we use the initial value for all training processes. We record the prediction errors for the testing data set and repeat times. From panel (b) in Fig. 2, our method and smoothly clipped absolute deviation obtain the least prediction errors.
5 Discussion
We discuss several important open problems on implicit regularization as our future directions. First, our theory heavily relies on the RIP of the design matrix, which is relatively strong comparing to the restricted eigenvalue condition (Bickel et al.,, 2009) as the minimal possible assumption in the literature. It would be interesting to theoretically investigate whether our results remain valid without the RIP. Second, it is of practical importance to study whether any computationally-efficient early stopping rule based on certain data-driven model complexity measures rather than the cross validation method can be applied to reliably and robustly select the tuning parameters, such as the iteration number and the step size in our algorithm.
Acknowledgment
Dr. Zhaoâs work was partially supported by National Science Foundation grant CCF-1934904. Dr. Yangâs work was partially supported by National Science Foundation grant DMS-1810831. We would like to thank the editor, associate editor and two anonymous reviewers from Biometrika for their careful comments and helpful suggestions that significantly improved the quality of the paper.
Supplementary material
Supplementary material contains detailed discussions on dynamic system interpretation of the implicit regularization, extension to adaptive step size and elastic net, proof of the theorems, extra simulations results beyond design with RIP, some reproducible simulation codes and a detailed comparison with Vaskevicius et al., (2019).
Appendix A Some additional results
A.1 Dynamical system interpretation for section 2.4
To emphasize the dependence of the solution on , we instead write as . For illustration purposes, we assume that the limiting point of this system is continuous and bounded as the initialization value , that is, both limits and exist in and are finite.
Let , then simple calculation leads to the relation
[TABLE]
Under the aforementioned assumption on the existence of limits as and , the preceding display implies one of the following three situations for each :
[TABLE]
Denote as the limit . Recall , and the previous three cases can be unified into
[TABLE]
This identity, together with the limiting point condition coincides with the Karush-Kuhn-Tucker condition for the basis pursuit problem (3).
Under the same aforementioned assumption, if there is additive noise on , then based on the integral mean value theorem, the preceding display implies one of the following three situations for each :
[TABLE]
A.2 Adaptive step size and variable selection
A nature extension of gradient descent is to assign different weights (step sizes) to different coordinates of , which is related to the adaptive Lasso (Zou,, 2006). It can be seen from the differential equation interpretation: by inserting a constant weighting matrix into the equation (4) in the main manuscript, we obtain the limiting dynamical system as
[TABLE]
Based on similar heuristic analysis as in Section 2.4 for the noiseless case, the limiting point of the dynamic system satisfies:
[TABLE]
which is the KKT condition for the dual form of the adaptive Lasso
[TABLE]
In the limiting case when the step size of a particular component tends to [math], we are equivalently adding an when . In contrast, if we apply a larger step size to , then tends to move faster and more freely in the parameter space, which is equivalent to a smaller penalty on . The original paper in Zou, (2006) constructed the weights based on the ordinary least square solution, which requires . In practice, when , we can construct weights through a preprocessing step. For example, variable screening can be applied to introduce sparse weights.
To enable variable selection in our method, we can perform a component-wise hard thresholding operation to the final estimator . Based on our theoretical analysis, since our method tries to shrink both weak signals and errors into very small order , it is more robust to false detections than other explicit regularizations when the same tuning parameter for noise level is applied. Let us consider a simple example to illustrate the basic idea: we set , , for and , , , and all other components are zeros in the data generating model with . Since the strength of the first components of truth is weak, it is hard to be detected by all methods we have tried. However, the effect of the weak signals on still pertains. In particular, when applying cross-validation, traditional penalized methods tend to over-select the predictors, leading to many false discoveries. In comparison, due to the implicit regularization, our method tends to be more robust to the weak signalsâour estimate is typically non-sparse. The effect of the non-detected weak signals can be distributed to all components of the estimated vector, and no component is particularly spiked. Consequently, our method tends to be more robust to false discoveries after applying the hard thresholding. The variable selection results are shown in Figure 3. As we can see, the Lasso can not detect the weak signal, and two components, indexed by , appear to be falsely detected through cross-validation (note that in Lasso, soft thresholding has already been applied). In contrast, in our method, most unimportant signals remain small. Performing hard thresholding with the same regularization parameter selected by the Lasso can erase all false detections, leading to the selection of strong signals only.
A.3 Empirical Comparisons between different early stopping criteria
We adopt the simulation framework S1-S4 for the strong-signal setting (only change the standard derivation to ) in section 4 to compare different early stopping criteria. We record the mean estimation errors averaging over trials and report the errors in figure 4.
A.4 Extension to Elastic Net
When both and type of regularizations are of interest to be applied, we can use the same Hadamard reparametrization , and the following objective function
[TABLE]
where controls the trade-off between and norms. The gradient of the problem is:
[TABLE]
In addition, denote , then the Hessian matrix is
[TABLE]
Since ridge regression parameterized by is convex, by using the same proof technique with lemma 1, we can obtain that does not have local maximums, all its local minimums are the global minimum, and all saddle points are strict saddle. In particular, is a global minimum of if and only if
[TABLE]
which reflects the effect of the penalty.
In addition, since we have the structure and , where with as the pseudo inverse of , we can adopt the dynamic system formula as follows:
[TABLE]
Then we have
[TABLE]
where , which reflects the effect of the penalty.
Therefore, the above formulation for can be seen as applying implicit regularization for the elastic net.
A.5 Choosing the step size and initial value
According to the main theorem, the parameters and are not sensitive to the final rate of estimation errors, as long as we choose and to be small enough, the optimal error rate in estimation can be obtained with properly tuned stopping time.
For , we can simply set so that will be term of a smaller order in the estimation error according to the main theorem. For choosing suitable , we can first choose based on a small subset of the whole training data with observation . Note that we require , so we can obtain .
A.6 Simulation result for median of mean prediction error
A.7 Simulations on the effect of error variance on the optimal stopping region
Let the truth be , , and . Note that we have and , by fixing the initial value and , we want to see the effect of the variance on the optimal stabilized region of the stopping time:
The above figure perfectly matches our theory: the start point of the optimal stopping time is at the same rate when the initial values, conditional numbers, and the smallest signals are the same; while the endpoint of the optimal stopping time is also affected by the error variance . While is becoming smaller, the region of the optimal stopping time will be more expanded.
Appendix B Proof of the results in the paper
B.1 Notation
Recall that and denote the vector -norm and -norm, respectively. Moreover, is the identity matrix in , and for any subset of , is the diagonal matrix with on the th diagonals for and [math] elsewhere. We use bold letter to denote an all-one vector. denote the -largest component of vector in absolute value. We use notation and to denote and up to some positive multiplicative constant, respectively. For two vectors and of the same dimension, we use and to denote element-wise and . Denote and be the maximal and minimal eigenvalues of matrix . Through this document, letters , and denote some constants whose meaning may change from line to line.
B.2 Some Useful Results
In our proof, we will constantly deal with the Hadamard product and the operation for two vectors . Therefore, we collect some useful properties in this section, some of which are consequences of the RIP condition.
The first property regarding the Hadamard product is a direct consequence of the Hölder inequality.
Lemma B.1**.**
For any two vectors and in , we have:
[TABLE]
Proof.
This follows since . â
The second lemma shows that under the RIP, the product can be well-approximated by for all sparse vectors and any vector .
Lemma B.2**.**
Let be a matrix in that satisfies -restricted isometry property (see Definition 2.1 in the paper). Then for any -sparse vectors and any in , we have:
[TABLE]
Proof.
Let be the diagonal matrix in with diagonal elements the same as components of correspondingly. Then can be represented as:
[TABLE]
where is a -dimensional vector whose -th component is and [math] elsewhere. Using the fact that satisfies -RIP and is -sparse, we have (see the remark right after Definition 2.1 in the paper):
[TABLE]
â
Our third lemma considers the case when and are both arbitrary.
Lemma B.3**.**
Let be a matrix in that satisfies -restricted isometry property. Then for any vectors , we have:
[TABLE]
Proof.
Since we can decompose , we have
[TABLE]
â
Our fourth and fifth lemma consider the concentration behaviors about the noise terms.
Lemma B.4**.**
Let is a vector of independent sub-Gaussian random variables with scales bounded by , and , with , then there exists a constant , with probability for any , such that for any :
[TABLE]
Proof.
By directly denoting as the constant in Hanson-Wright inequality (Rudelson et al.,, 2013):
[TABLE]
â
Lemma B.5**.**
Let is a vector of independent sub-Gaussian random variables with scales bounded by , and all norm of column vectors of are normalized to , then with probability such that:
[TABLE]
Proof.
is the maximum over sub-Gaussian random variables with scales bounded by , so considering the union bound:
[TABLE]
the result is obtained by taking . â
Appendix C Proofs of Lemmas, Propositions and Theorems in the Paper
We prove the results in the paper one by one. We assume the largest signal strength is at the most constant level for simplicity, which implies . Note that the condition can be satisfied under our two-step post debiasing process.
C.1 Proof of Lemma 2.1
Let . Since the least square problem is a convex optimization, it does not have local maximum and any local minimum is a global minimum.
Now suppose has a local maximum , meaning that there exists an open set centering at such for any , . It is easy to verify that the set is also an open set in centering at . Consequently, is also a local maximum of , which is a contradiction.
Due to the same argument, any local minimum of corresponds to a local minimum of . Since all local minimums of are global and satisfies X^{T}\big{(}X\bar{\beta}-y\big{)}=0, must be a global minimum of and satisfies X^{T}\big{(}X(\bar{g}\circ\bar{l}\,)-y\big{)}=0.
Now let us prove the last part concerning saddle points. According to the previous paragraph and the discussion after lemma 2.1, a point is a saddle point (that is, neither local minimum nor local maximum) if and only if there exists some such that
[TABLE]
where R^{\dagger}=X^{T}\big{(}X(g^{\dagger}\circ l^{\dagger})-y\big{)}\in\mathbb{R}^{p} (since otherwise this saddle point would satisfy X^{T}\big{(}X(\bar{g}\circ\bar{l}\,)-y\big{)}=0, which is the sufficient and necessary condition for a global minimum). For any vector , we use Diag to denote the -by- diagonal matrix whose diagonals are components of . By direct calculations, we can express the -by- Hessian matrix of at as
[TABLE]
Consequently, it is easy to verify that
[TABLE]
where denotes the vector whose -th component is one and zero elsewhere for , and sgn is the sign function. Consequently, there exists some nonzero vector such that the evaluation of the quadratic form , implying \lambda_{\min}\big{(}\nabla^{2}_{g,l}f(g^{\dagger},l^{\dagger})\big{)}<0.
C.2 Proof of Lemma 2.2
According to lemma 2.1, all saddle points of are strict, and all local minimums are global. Therefore, we can apply Theorem 4 in Lee et al., (2016) to finish the proof, which states that for any twice continuously differentiable function, if it has only strict saddle points, then gradient descent with a random initialization and sufficiently small constant step size almost surely converges to a local minimizer.
C.3 Gradient descent error dynamics
In this subsection, we provide two propositions for characterizing the dynamics of the gradient descent algorithm for minimizing the Hadamard product parametrized quadratic loss defined in (2). Theorem 3.1 in the main content is a direct consequence of these two propositions. Due to the lack of knowledge on signs of nonzero components in signal , the proof of Theorem 3.1 has the extra component in showing that tends to shoot towards the right direction (positive versus negative) during the first âburn-inâ stage before entering the second geometric convergence region. To deal with this extra complication, we introduce another reparametrization:
[TABLE]
Notice that , implying that the sign of each component is determined by the relative magnitude of and . These new variables are especially convenient for understanding the dynamics of the sign of due to a more tractable updating formula:
[TABLE]
More precisely, in the first âburn-inâ stage, we show that each component of the strong signal part increases at an exponential rate in until hitting , while the weak signal and error part remains bounded by .
Proposition 1** (Stage one dynamics).**
Under the assumptions of Theorem 3.1, with probability at least , there is a positive constant , such that for each , we have:
[TABLE]
Let \Theta_{LR}^{G}=\{(a,\,b)\in\mathbb{R}^{2p}:\,\mbox{sign}(\beta^{\ast}_{S_{1}})\circ(a_{S_{1}}^{2}-b_{S_{1}}^{2})\geq(m/2)\,\mathbf{1},\,\max\big{\{}\|a_{S_{1}}\|_{\infty},\,\\ \|b_{S_{1}}\|_{\infty}\big{\}}\lesssim 1,\,\|a_{S_{1}}^{2}-b_{S_{1}}^{2}-\beta^{\ast}_{S_{1}}\|\lesssim\sqrt{s_{1}},\,\max\big{\{}\|a_{S_{1}^{c}}\|_{\infty},\,\|b_{S_{1}^{c}}\|_{\infty}\big{\}}\lesssim 1/p\} as a good parameter region for . Proposition 1 tells that for all , iterate satisfies all constraints of except , but will enter this good region in at most iterations.
The second stage starts when first enters . The next result summarizes the behavior of the gradient dynamics in the second stageâonce it enters , it will stay in for a long time, where converges toward at a linear rate up to the statistical error and then stay in an neighborhood of up to time \Theta\big{(}1/(\eta\sigma\sqrt{\log p/n})\big{)}.
Proposition 2** (Stage two dynamics).**
Under the assumptions of Proposition 1, if with , then , and satisfies
[TABLE]
A combination of these two propositions leads to a proof of Theorem 3.1. Therefore, the time interval \big{[}c_{3}\log(1/\alpha)/(\eta m),\,c_{4}/(\max\{\log(1/\alpha)/(\eta m),1/(\eta\sigma\sqrt{\log p/n})\}\big{]} would be the theoretical âbest solution regionâ corresponding to the stabilized region in Figure 1.
C.4 Proof of Proposition 1
Throughout the proof, we let and
for any . Recall we have the updating formula:
[TABLE]
with . We use induction to show that for each ,
[TABLE]
where the set is the union of weak signals and errors. When , we have , . Therefore, under the assumption , we have , , and similar bounds for . Now suppose for time , we have (16)-(18). We divide into 3 steps to show the same bounds hold for time . Note that the following analysis is still based on the conditional events:
[TABLE]
where the first concentration is based on lemma B.4 and the RIP of . Also recall the following conditions on the step size , and RIP constant : , and .
Proposition 3** (Step 1: Analyzing Error Dynamics:).**
Based on RIP and the assumption for constants in theorem 3.1, assume the induction hypothesis (16)-(18) hold for time with , then we have:
[TABLE]
When , let , then
[TABLE]
Proof.
By the definition of and , we have the updating rules:
[TABLE]
Denote . Note that . For , we bound the term in the updating formula:
[TABLE]
where in step we use the fact that ; step follows by lemma B.2:
[TABLE]
lemma B.3:
[TABLE]
and lemma B.5:
[TABLE]
In step , we applied the induction hypothesis (16). Finally, note that and by the minimal strength of , then an application of the triangle inequality leads to the claimed bound on and then similarly .
Similarly, for ,
[TABLE]
where step follows from the orthogonality between and . The last inequality follows from the induction hypothesis and the definition of weak signals.
Then, when :
[TABLE]
where the last inequality follows from and .
[TABLE]
which is based on , since and the signals in set are weak. A similar proof can be applied to the dynamics of . â
Step 2: Analyzing Signal Dynamics: The proof becomes more complicated for the signal dynamics. Our goal is to show:
- âą
if the true signal is positive and , then increases and decreases, both in an exponential rate. Overall, the sign of th component of tends to grow to positive in an exponential rate;
- âą
if the true signal is negative and , then decreases and increases, both in an exponential rate. Overall, the sign of th component of tends to fall to negative in an exponential rate.
Recall that we assumed that the RIP constant , step size , and .
Proposition 4**.**
Under assumptions in theorem 3.1 and the induction hypothesis (16)-(18) at , when , then for any :
[TABLE]
Proof.
First we approximate by based on the RIP condition via lemmas B.2 and B.3,
[TABLE]
implying that under the condition , we have
[TABLE]
where the last inequality uses and .
In order to analyze , let us focus on and , separately. According to the updating rule of , we have
[TABLE]
where recall that for a vector , denote the sub-vector of with indices in . Applying lemmas B.2 with , we obtain
[TABLE]
where in step we used sand . Since , , we have for . Therefore, we can obtain an element-wise bound for ,
[TABLE]
as . Equivalently, we can write
[TABLE]
Now let us divide into two cases depending on the sign of :
Case : When , since , we have by equation (21),
[TABLE]
where the last inequality follows since and . Similarly, we can analyze to get
[TABLE]
Therefore, increases at an exponential rate faster than the noise term while decreases to zero at an exponential rate, and when increases to , decreases to correspondingly. A combination of these two leads to the first claimed bound for .
Case : The analysis for the case is similar: when , we have:
[TABLE]
which leads to the second claimed bound for . â
As a consequence of the proof in this step, after at most iterations, we are guaranteed to have have with and .
Step 3: Prove Induction Hypothesis: Our last piece shows that the induction hypothesis is kept for . This proposition combined with propositions 3-4 leads to a proof for proposition 1 in the paper.
Proposition 5**.**
Under assumptions in theorem 3.1, the induction hypothesis (16)-(18) hold at , we have:
[TABLE]
Then for , the above induction still holds.
Proof.
Our induction hypothesis
[TABLE]
implies and . By updating rules:
[TABLE]
and our condition . For weak signals and error components , we have:
[TABLE]
where in step we apply . When , since , when :
[TABLE]
For strong signal component , , we have the component-wise bound , then
[TABLE]
where follows by proof of proposition 4, since the smaller one of will be its initial value decreasing at an exponential rate until the rate . If such that increases to large numbers while converges to zero, then based on , we can apply equation (21) to get element-wise upper bound for :
[TABLE]
Since , and the function is nonnegative when . Therefore, as long as , we always have . Similarly, . Then we have
[TABLE]
Now we have:
[TABLE]
These inequalities also imply and . Then similarly for , since the error updating still holds, equation (21) can still help us to provide the same proof, so we also have:
[TABLE]
â
C.5 Proof of Proposition 2
Proof.
First, we can decompose into 5 parts according to the order of :
[TABLE]
First note that we have equation (22), so we have the decomposition:
[TABLE]
which is based on the orthogonality between different sub vectors and the induction hypothesis, so we only need to calculate the errors induced by the strong signals.
For the gradient, we have:
[TABLE]
Note that is sparse, so we have the following core error term characterization:
[TABLE]
Similarly, our analysis for errors is based on equation (25) and (26).
For the order one term, recall the bound
[TABLE]
where is based on Cauchy-Schwarz inequality and we use the fact , lemma B.1, equation (26).
For the rest of terms, we can see they are all based on the bound:
[TABLE]
First consider the order two terms and . For the former one, since , we have
[TABLE]
where we use the condition . The analysis of the latter one is similar based on and :
[TABLE]
For the order three terms, using inequality , we obtain:
[TABLE]
For the order four terms, we have
[TABLE]
Putting all pieces together, and using our condition , we can obtain
[TABLE]
Since for , we have . Therefore, under our conditions on with , we have:
[TABLE]
â
C.6 Proof of Theorem 3.1
As mentioned after the theorem in the paper, we divide the convergence into two stages. In the first âburn-inâ stage, we show that the smallest component of the strong signal part increases at an exponential rate in until hitting ,
[TABLE]
In the second stage, after iterate enters , we have the geometric convergence up to some high-order error term
[TABLE]
When , the convergence is in the first stage; and proposition 1 implies that the iterate enters (corresponding to the second stage) in at most iterations. proposition 1 and proposition 2 together imply for any ,
[TABLE]
For any , since is still controlled, combined with bound (C.5), we have:
[TABLE]
Appendix D Simulations Beyond RIP Conditions
In this section, we conduct some numerical experiments to illustrate the performance of our method. In our first example, we study the performance as we change the initialization under settings where the RIP may or may not hold. In the second example, we study the behavior of the algorithm when the null space property (Cohen et al.,, 2009) is violated, under which the -norm minimizer in the basis pursuit problem differs from the sparsest solution (the -norm minimizer).
We consider two classes of random designs with rows , where population p.s.d covariance matrix is
[TABLE]
for , where denotes the indicator function of event . We choose sample size , sparsity level , and signal dimension . For the independent (correlated) design, the RIP is satisfied (fails) with high probability for some small . In both scenarios, we choose true signal , and set . When implementing gradient descent, we choose step size for the independent (correlated) design, , and stopping threshold . Figure 6 shows the estimation error versus in log-log plots. As we can see, they tend to have a linear trend under the log-scale, which is consistent with our theoretical error bound estimate in Section 4. In addition, in the correlated design scenario where the RIP does not hold, the algorithm can still recover as , albeit under a slower convergence (due to a smaller allowable step size and a larger condition number of ). This observation provides evidence of the correctness of our informal statement made at the beginning of Section 3, even without RIP condition. We leave the proof of this conjecture open.
In this example, we study the empirical performance of the algorithm when the least -norm in the basis pursuit problem (1) is not the sparsest solution of (the null space property is violated). In particular, we choose
[TABLE]
so that . It is easy to verify that for this example, the sparsest solution of is , while the least -norm solution is . We use the same setting as before for implementing the algorithm with . Table 3 reports final outputs of the algorithm. Due to our intuition in Section 3, as expected, the algorithm still converges to the least -norm solution instead of the least -norm solution . Again, the estimation error decreases as the initialization level decreases.
Appendix E Simulation results for Algorithm 2
We adopt the simulation settings in S1-S4 for strong signals, except that now . For Algorithm 2, we use MCP as the pilot estimator and then use gradient descent for post-debiasing. The below is the median of standardized RMSE of the repeated simulations. The results show that Algorithm 2 could achieve comparable results with Algorithm 1.
Appendix F Point-to point Comparisons with Vaskevicius et al., (2019)
Theoretical results. We both consider investigating implicit regularizations in high dimensional regressions under the RIP assumption. We adopt parametrization as the problem formulation since it has nice interpretations both from gradient dynamic systems and implicit regularizations in settings. In the proof for general signals, we introduced the equivalent parametrization , where and , as a proof technique to facilitate our proof (which is NOT inspired by their work, but appeared in our first version before their work appeared). However, their work directly applies the parametrization without providing detailed intuitions and connections.
Compared with previous work (Li et al.,, 2018) in the area, both our works make three new contributions in the theoretical analysis: 1. adopting the parametrization such that linear regression models can be properly analyzed; 2. dividing the signals into strong and weak ones and connecting signal strength with the minimal number of iterations in order to recover them; 3. showing that the final error rate can be adaptively improved via early stopping. For the detailed proofs, we follow (Li et al.,, 2018) by also analyzing the convergence under the RIP assumption, while Vaskevicius et al., (2019) developed an improved analysis for handling noises but still requires the RIP condition. Below we compare the assumptions made in the two works.
Our conditions on the RIP constant, step size and signal strengths are more stringent than theirs. There are two reasons to explain this difference. First, our target is to show that the estimator achieves the benchmark optimal rate of ( are number of strong and weak signals) â where the rate becomes parameter root rate if there is no weak signals (i.e. ), which is strictly better than the LASSO. To achieve this goal, we need to make some stronger assumptions. In comparison, their rate is , which is slightly worse in terms of the dependence on the sparsity level . Second, Vaskevicius et al., (2019) developed an improved analysis which makes the conditions less stringent than ours. However, our more stringent conditions do not affect the application of our method if we mainly apply implicit regularization as a way for refining a rough initial estimator through our newly added two-stage debiasing approach.
Another important theoretical distinction is the different stabilized optimal stopping regions : when the estimator achieves the optimal rate. Specifically, the theorem in Vaskevicius et al., (2019) implies that both and are at the same rate (c.f. Corollary 2,3 in Vaskevicius et al., (2019), where is the initial value and is the signal conditional number). The results also mean that is independent of when is at the constant level. In contrast, our theorem indicates that the optimal region is and ( is the smallest signal), where the smaller value in will result in a larger , therefore a more expanded optimal stopping region. Our theoretical results match the simulation studies in section A.6 in the supplementary material better. Moreover, similar results (where and have different rates) for the optimal stopping region can also be seen in the related literature (e.g., Theorem 3.6 in Fan et al., (2020)). 2. 2.
Computation.
For our vanilla algorithm, we did not elaborate much on the computational aspects, while their work designed an adaptive algorithm that performs better in computation, which is in the number of iteration for convergence. In the revised paper, we added a different algorithm to achieve the same computational efficiency based on pre-processing the data. Moreover, we want to point out that the improved algorithm in Vaskevicius et al., (2019) requires applying different step sizes (that also depend on time) for updating different components in . Their variation on the algorithm loses the implicit regularization interpretation: its accompanied (continuous-time) limiting gradient dynamical system will effectively multiply an extra diagonal matrix in front for both and as compared to our limiting dynamical system (4) in the paper. Consequently, it is not clear whether their improved algorithm still converges to the minimal solution, for example, in the noiseless setting without the RIP condition. This means that their improved algorithm may truly rely on the RIP condition, not like our method where the RIP condition is mainly made for technical purposes (c.f. the heuristic illustration in Section 2.4 and numerical results in Section D of the supplement for details). 3. 3.
Methodology. From the gradient dynamic systems point of view, our work has many extensions in the proposed methodology. First, assigning different weights (step sizes) to different coordinates will result in similar results as the adaptive lasso; moreover, with the help of one reviewer, we also bring the extension to elastic-net, where an additional penalty is incorporated into the objective function. Second, we also provide an empirical comparison of different stopping criteria. Third, we discuss how to conduct variable selection based on the output (via proper thresholding (Appendix A.2). Last but not least, we propose the use of our implicit regularization method as a refining technique for improving explicit regularization based estimators, which boosts their convergence rate from to when there is no weak signal. 4. 4.
Data analysis. Both our works show numerical evidence to support our main theorems. In addition to that, both of us provide simulations that implicit regularization still works when a much relaxed restricted eigenvalue condition on the design matrix instead of the RIP holds. In our work, we also show that when null space property does not hold, our algorithm can not converge to the sparest solution even under the noiseless case, which provides high-level intuitions that the sparsity of the solution is induced from the norm instead of type of method (like forward-backward regression). The high-level messages are closely connected with our dynamic system, and norm regularized interpretations. 5. 5.
Others. In our work, we provide an analysis of the optimization landscape in the noiseless regime (Lemma 1 & 2) by showing that the non-convex objective function after the Hadamard product parametrization does not have local maximums, and all its local minimums are global. Consequently, the gradient method for optimizing the objective function with random initialization almost surely converges to a global minimum. We also perform simulations and real data analysis to study variable selection and illustrate the success of some empirical rules (using cross-validated from the lasso) for achieving good variable selection performance. In comparison, Vaskevicius et al., (2019) didnât discuss/study variable selection.
Appendix G Discussion
We discuss several important open problems on implicit regularization as our future directions. First, our theory heavily relies on the RIP of the design matrix, which is relatively strong comparing to the restricted eigenvalue condition (Bickel et al.,, 2009) as the minimal possible assumption in the literature. It would be interesting to theoretically investigate whether our results remain valid without the RIP. Second, it is of practical importance to study whether any computationally-efficient early stopping rule based on certain data-driven model complexity measures rather than the cross validation method can be applied to reliably and robustly select the tuning parameters, such as the iteration number and the step size in our algorithm.
Acknowledgment
Dr. Zhaoâs work was partially supported by National Science Foundation grant CCF-1934904. Dr. Yangâs work was partially supported by National Science Foundation grant DMS-1810831. We would like to thank the editor, associate editor and two anonymous reviewers from Biometrika for their careful comments and helpful suggestions that significantly improved the quality of the paper.
Supplementary material
Supplementary material contains detailed discussions on dynamic system interpretation of the implicit regularization, extension to adaptive step size and elastic net, proof of the theorems, extra simulations results beyond Restricted Isometry Property, some reproducible simulation codes and a detailed comparison with Vaskevicius et al., (2019).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Bickel et al., (2009) Bickel, P. J., Ritov, Y., Tsybakov, A. B., et al. (2009). Simultaneous analysis of lasso and dantzig selector. Ann. Statist. , 37(4):1705â1732.
- 2Breheny and Huang, (2011) Breheny, P. and Huang, J. (2011). Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Ann. Appl. Statist. , 5(1):232.
- 3BĂŒhlmann et al., (2014) BĂŒhlmann, P., Kalisch, M., and Meier, L. (2014). High-dimensional statistics with a view toward applications in biology. Annu. Rev. Stat. Appl. , 1:255â278.
- 4Candes and Tao, (2007) Candes, E. and Tao, T. (2007). The dantzig selector: Statistical estimation when p is much larger than n. Ann. Statist. , 35(6):2313â2351.
- 5Candes, (2008) Candes, E. J. (2008). The restricted isometry property and its implications for compressed sensing. C. R. Math. , 346(9-10):589â592.
- 6Chen et al., (2001) Chen, S. S., Donoho, D. L., and Saunders, M. A. (2001). Atomic decomposition by basis pursuit. SIAM Rev. , 43(1):129â159.
- 7Cohen et al., (2009) Cohen, A., Dahmen, W., and De Vore, R. (2009). Compressed sensing and best k-term approximation. J. Am. Math. Soc. , 22:211â231.
- 8Donoho, (2006) Donoho, D. L. (2006). Compressed sensing. IEEE T. Inform. Theory , 52(4):1289â1306.