Implicit Regularization for Optimal Sparse Recovery

1 Introduction

• Sparse regression seeks a k-sparse d-dimensional signal from n observations, often with d greater than n, where explicit ℓ1 methods provide a standard approach.In noisy settings, the introduction identifies the minimax squared-error rate as a central statistical target.Source: §1 · p. 1
• Implicit regularization for sparse regression had lacked results despite extensive study of solver-induced regularization for ℓ2 problems.The introduction highlights cross-validation for explicit regularization parameters as an added computational burden.Source: §1 · p. 2

1.1 Our Contributions

• The paper proposes non-convex gradient descent on w = u⊙u−v⊙v, whose multiplicative coordinate updates support sparse recovery under the RIP.The authors present this as the first statistically and computationally optimal implicit-regularization algorithm for sparse linear regression under the RIP.Source: §1.1 · p. 2
• Increasing, coordinate-wise step sizes use refined solution estimates to reach minimax rates with data-read-scale computation up to poly-logarithmic factors.For sub-Gaussian noise, the paper gives an iteration bound for this scheme and states that each iteration costs O(nd).Source: §1.1 · p. 3
• Simulations show gradient descent fitting coordinates one at a time rather than all at once as lasso tends to, with performance eventually matching oracle least squares.Across hard-to-easy instances, the paper reports a phase transition not observed for lasso.Source: §1.1 · p. 3

1.2 Related Work

• Related work covers explicit ℓ1 sparse-recovery methods and implicit-bias studies, while concurrent work [56] is distinguished by this paper’s computational optimality and treatment of general noise.The paper also relates its empirical findings to ℓ0 regularization.Source: §1.2 · pp. 3–4

2 Model and Algorithms

• The RIP requires squared norms of the normalized design’s action on sparse vectors to stay within factors 1−δ and 1+δ of their original squared norms.The text notes that verifying RIP for a given matrix is NP-hard, while random matrices satisfy it with high probability under stated conditions.Source: §2 · pp. 4–5
• Updating the factors u and v rather than w makes the optimization path contain sparse solutions, unlike direct gradient descent on w.The model uses w = u⊙u−v⊙v and gradient descent on the unregularized, underdetermined least-squares objective.Source: §2 · p. 5
• The section considers constant-step gradient descent and an adaptive scheme that increases step sizes for sufficiently small coordinates over time.Algorithm 1 uses a constant step size; Algorithm 2 doubles selected step sizes after intervals.Source: §2 · p. 5
• The effective condition number treats coefficients below the target precision or maximum noise term as zero, motivating an initialization-size restriction.The setup imposes restrictions on the RIP constant and initialization size.Source: §2 · p. 5

Theorem 1

• Theorem 1 establishes rates in ℓ∞ norm under a restricted-isometry condition, rather than only the ℓ2 rates commonly derived for lasso.The theorem's bounds depend on initialization, step size, and stopping time; α controls off-support coordinates at stopping.Source: §3 · pp. 5–6

Parameter selection

• A gradient step can estimate ||w*||∞ from the data within a factor of two, enabling data-based choices of initialization and learning rate.The proposed learning rate is η ≍ 1/||w*||∞; the iterates are used as statistical estimators, and Algorithm 2's iterates would explode if run indefinitely.Source: §3 · pp. 6–7

Noiseless recovery

• In the noiseless setting, Algorithm 1 achieves squared ℓ2 error at most approximately kε² after O(κ_eff(w*) log α^-1) iterations.This guarantee holds under Theorem 1's assumptions with the step size specified by Theorem 2.Source: §3 · p. 6

Noisy-setting rates

• With independent sub-Gaussian noise and sub-linear sparsity, Algorithm 1 achieves minimax-optimal ℓ2 rates, although exact reconstruction is not possible in the noisy setting.The stated sub-linear sparsity regime has k^γ ≤ d for some γ > 1.Source: §3 · pp. 6–7

Instance adaptivity

• When w*_min ≳ ||Xᵀξ||∞/n, the bounds become dimension-independent; with sub-Gaussian noise, the transition occurs at n ≳ (σ² log d)/(w*_min)^2.The paper contrasts this behavior with lasso performance, which it says deteriorates with d regardless of problem difficulty.Source: §3 · p. 7

Algorithm 2

• Algorithm 2 attains Algorithm 1's statistical bounds and is computationally optimal up to poly-logarithmic factors.Its stopping time is O(log κ_eff log α^-1) iterations.Source: §3 · p. 7

Simplified setting

• The simplified proof setting assumes XᵀX/n = I; the paper notes that noise and restricted isometry instead of isometry require a more involved analysis.The simplified setting permits tighter bounds than Theorem 1.Source: §4 · pp. 7–8

Reparameterization

• The reparameterized updates preserve u_t ⊙ v_t ≤ u_0 ⊙ v_0, so at most one of each coordinate pair grows beyond initialization.The paper interprets u_t as serving positive true coordinates and v_t as serving negative ones.Source: §4 · p. 7

Error decomposition

• The proof separates w_t into support signal s_t and off-support error e_t, then writes their updates separately.The displayed updates express both sequences using the interaction term p_t and error term b_t.Source: §4 · p. 8

Noise error

• The noise-driven error b_t is irreducible at scale ||Xᵀξ||∞/n, while signal coordinates above that scale can be fitted to obtain minimax-optimal rates.If w*_min ≳ ||Xᵀξ||∞/n, the resulting error depends on the noise restricted to the support, yielding dimension-independent bounds.Source: §4 · p. 8

Interaction error

• The interaction error p_t is controlled using RIP, and the proof relies on it decreasing before off-support errors grow too large or signal components become too small.The paper describes the signal as converging exponentially toward the true parameters in this argument.Source: §4 · p. 8

Tuning learning rates

• Coordinate-wise increasing step sizes wait for p_t to become small before accelerating coordinates, avoiding the instability that can occur with larger steps too early.The analysis contrasts this scheme with oracle coordinate-wise step sizes in the identity-design setting.Source: §4 · pp. 8–9

5 Simulations

• Small initialization lets exponentially growing coordinates on the true support dominate those outside it, enabling the algorithm to fit coordinates one at a time.The simulations illustrate this effect for k = 5, n = 100, η = 0.05, and σ = 0.5; the two right-hand plots use different x-axis scales because α differs.Source: §5 · p. 9, Figure 1 · p. 9
• At κ = 64, Algorithm 2 fits successive coordinates in approximately equal iterations, unlike Algorithm 1, which takes about twice as long for each coordinate; its total is O(log κ log α−1) iterations.The comparison uses n = 250, k = 7, and a step size of 1/(20 · 26) for both algorithms.Source: §5 · p. 9, Figure 2 · p. 10
• When w*_min ≳ ||X^Tξ||∞/n, the method matches oracle least squares, whose performance is independent of d, while lasso deteriorates as d increases.In simulations varying γ, σ, and n, the method approaches oracle-like performance once the stated threshold is met.Source: Figure 3 · p. 10, §5 · p. 10
• A log d-dependent bias illustrated in Figure 4 is absent at the algorithm’s optimal stopping time.The supplied passage does not identify the estimate associated with that bias.Source: §5 · p. 10

6 Further Improvements

• In a high-SNR setting, gradient descent recovers coordinates on S without visible bias, unlike the lasso’s dimension-dependent bias.Figure 4 illustrates the contrast.Source: Figure 4 · p. 11
• The authors suggest that both the sample-complexity constraint and reliance on RIP may be analysis artifacts rather than inherent algorithmic limitations.They identify restricted eigenvalue conditions as a potential alternative to RIP.Source: §6 · p. 11
• For sub-Gaussian isotropic designs, the stated RIP scaling requires n ≳ k^2 log(ed/k) samples, yielding sub-optimal sample complexity.The section identifies the RIP-parameter constraint as the source of this sample-complexity issue.Source: §6 · p. 11
• Figure 5 suggests linear sample complexity in k suffices for the method to match and eventually exceed the lasso in ℓ2 error, and for the main theorems to hold.The left plot supports the comparison with the lasso; the right plot suggests linear sample complexity suffices for the theorems.Source: §6 · p. 11
• At µ = 0.5, where RIP fails but RE holds with high probability [50, Chapter 7], experiments show fast rates and eventual outperformance of the lasso.Figure 6 uses Gaussian designs with covariance (1 − µ)I + µ11T.Source: Figure 6 · p. 12, §6 · p. 12

A.1 Set Up and Intuition

• Stopping before the error sequence exceeds an initialization-dependent threshold recovers support coordinates sufficiently above noise while keeping off-support coordinates near zero.This is the stated role of the stopping rule in the proof intuition.Source: §A.1 · p. 17
• The analysis identifies a hard early-iteration case: when γ||w*||∞ exceeds a support coefficient, errors can grow rapidly while that signal coordinate shrinks.This follows from the worst-case control on the convergence-distance error sequence.Source: §A.1 · p. 17
• The iterates start with u0 = v0 = α for small positive α, since setting both to zero would place the optimization at a saddle point.Source: §A.1 · p. 18

A.2 The Key Propositions

• Proposition 1 allows the signal sequence to fit w⋆ to precision proportional to ζ while keeping the error sequence arbitrarily small, when the RIP constant is sufficiently small.It applies when the perturbation sequence is bounded in ℓ∞ by ζ.Source: §A.2 · p. 18
• In Proposition 2, once the error term is below α, signal convergence no longer depends on α, allowing the signal iterates to approach w⋆ while the error iterates remain arbitrarily small.The proposition fixes ε > 0 and assumes a sufficiently small step size.Source: §A.2 · pp. 18–19

A.3 Technical Lemmas

• The technical lemmas bound the error sequence and control perturbation terms under RIP, including the maximum noise term.They provide bounds for the ℓ∞ norms of the perturbations and of (1/n)X^T X e_t, and an upper bound on the maximum noise term.Source: §A.3 · pp. 19–20

A.4 Proof of Theorem 1

• Theorem 1’s proof separates hard and easy regimes; in the easy regime, it uses the lemmas to maintain a small error sequence while applying Proposition 1.The proof identifies an error bound to maintain over iterations in the easy regime.Source: §A.4 · p. 21

A.5 Proofs of Corollaries

• The corollary proofs derive the zero-noise case directly from Theorem 1 and use a high-probability noise bound for the other corollaries.The cited proof passages give probability at least 1 − 1/(8d^3) for the bound used in Corollary 2.Source: §A.5 · pp. 21–22

B Understanding Multiplicative Update Sequences

• The appendix analyzes multiplicative-update sequences and uses their technical lemmas to prove Propositions 1 and 2.Source: B Understanding Multiplicative Update… · p. 22

B.1 Basic Lemmas

• When the signal target exceeds the noise target, a sufficiently small initialization can make the signal hit its target tolerance before the noise sequence crosses its threshold.The illustrated effect of decreasing initialization is that both hitting times move later, but the noise hitting time moves later faster.Source: §B.1 · pp. 22–23
• For small enough step sizes, multiplicative iterates are monotone and bounded, and the lemmas establish exponential-time bounds for approaching the target.The bounds distinguish iterates starting below the target from those starting above it.Source: §B.1 · pp. 23–25
• Combining the earlier iteration bound with the time needed to reach x⋆/2 gives a sufficient total iteration count.The exact total bound is garbled in the supplied passages.Source: §B.1 · pp. 25–26
• While x_t < x⋆/2, the iterates satisfy x_(t+1) ≥ x_t(1 + 2ηx⋆)^2, yielding x_t ≥ x⋆/2 after a sufficient number of iterations.The supplied bound is t ≥ 5/(16ηx⋆) log(x⋆/(2x_0)).Source: §B.1 · pp. 25–26
• Under η ≤ 1/(8x⋆) and x⋆ ≥ 12y⋆ > 0, the cited lemmas establish conditions for the target inequality involving T and α.The passage says these conditions answer the section’s opening question, but the inequality is truncated in the supplied text.Source: §B.1 · p. 26

B.2 Dealing With Bounded Errors

• Bounded-error convergence can only be expected to precision B; the stated lemmas require B ≤ x⋆/5.The analysis also gives exponential approach to within B + ε under specified initialization and step-size conditions.Source: §B.2 · pp. 27, 30
• With bounded errors, noisy iterates are sandwiched between comparison sequences targeting x⋆−B and x⋆+B.Monotonicity holds outside the B-tube, and iterates stay inside once they enter it.Source: §B.2 · pp. 27–28

B.3 Dealing With Errors Proportional to Convergence Distance

• Errors proportional to convergence distance are handled through their expected exponential decay, while the analysis separately controls initial signal shrinkage.The proof uses halving error bounds over doubling time intervals and relates error growth control to lower bounds on signal terms.Source: §B.3 · pp. 30, 32

B.4 Dealing With Negative Targets

• For negative targets, the analysis uses the squared-difference parametrization by treating positive and negative coordinate sequences together.At least one of each coordinate pair remains close to its initialization, allowing the general parametrization to be handled almost as the positive-target case.Source: §B.4 · pp. 32–33

B.5 Proof of Proposition 1

• When the RIP constant is sufficiently small, the error sequence proportional to convergence distance can essentially be ignored: the fitted signal reaches precision ζ while the other error remains arbitrarily small.The proof reduces the maximum convergence distance by half at successive induction steps while shrinking the corresponding error bound.Source: §B.5 · p. 34

B.6 Proof of Proposition 2

• Under sub-Gaussian noise, the resulting ℓ2 parameter-estimation error can depend on log k instead of log d, yielding oracle-like performance.The proof first controls convergence distance and then applies coordinate-wise convergence bounds on the true support.Source: §B.6 · p. 38

C Missing Proofs from Section A.3

• Section C supplies proofs for the technical lemmas stated in Section A.3.For the RIP-based bounds, the proof uses the (k + 1, δ)-RIP on the joint support of a k-sparse vector and a coordinate.Source: C Missing Proofs from Section A.3 · p. 39, §C.3 · p. 40

C.4 Proof of Lemma A.4

• The proof of Lemma A.4 obtains its stated result directly from the (k + 1, δ)-RIP.The passage also invokes Cauchy–Schwarz for the corresponding pairwise bound.Source: §C.4 · pp. 40–41

• The design columns satisfy ∥Xi∥2/√n ≤ C, and the noise components are independent σ²-subGaussian random variables.Source: §C.5 · p. 41
• The proof invokes a standard result for any ε > 0.Source: §C.5 · p. 41
• The stated event holds with probability at least 1 − 1/(8d³).Source: §C.5 · p. 41

• The proof starts from w0 = 0 and refers to the parameter updates in equations (3) and (4).Source: D Proof of Theorem 2 · p. 41
• One gradient descent iteration recovers w⋆max up to multiplicative constants.The proof identifies this computation as the way to recover the factor needed in the step-size argument.Source: D Proof of Theorem 2 · p. 42
• The new step size satisfies Theorems 1 and 3 while being at most two times smaller than required.Source: D Proof of Theorem 2 · p. 42

E Proof of Theorem 3

• Propositions 3 and 4 adapt the earlier propositions to different step sizes along each dimension.Proposition 3 considers Algorithm 2 with τ = 640; Proposition 4 states conditions on the error sequences.Source: E Proof of Theorem 3 · p. 42
• Theorem 3’s proof uses Propositions 3 and 4 in place of Propositions 1 and 2, respectively.The proof notes that the preconditions for Proposition 4 are met by Proposition 3.Source: E Proof of Theorem 3 · p. 43

E.1 Proof of Proposition 3

• A constant step size can make iterates explode if made larger, while doubling the iteration count across induction steps creates inefficiency for smaller signal coordinates.The proof proposes a different step size for each coordinate to address that inefficiency.Source: §E.1 · p. 43
• The increasing-step-size argument extends to the general parametrization w_t := u_t ⊙ u_t − v_t ⊙ v_t.The appendix says this extension yields a proposition that allows the increasing-step-size scheme.Source: §E.1 · p. 46

F Gradient Descent Updates

• The appendix derives gradient-descent updates for w = u ⊙ u − v ⊙ v.Source: F Gradient Descent Updates · p. 46

• In both papers, α controls noiseless ℓ2 or reconstruction error, and the iteration count is affected by O(log 1/α).Source: G Comparing Assumptions to [30] · p. 47
• The conditions on α and η in [30] are stronger than those required here, and the results do not follow from [30] by considering diag(w⋆).Source: G Comparing Assumptions to [30] · p. 47
• The theorem requires δ to scale as 1/log κ and permits data-based step-size selection without knowing κ or k.Source: G Comparing Assumptions to [30] · p. 47
• The paper states that [30]’s conditions imply δ = O(n^-3/2), generally impossible to satisfy, for example with i.i.d. Gaussian design entries.It identifies moving κ-dependence into a logarithmic factor as key to handling the general noisy setting.Source: G Comparing Assumptions to [30] · p. 47

H Comparing Our Results to [56]

• At w⋆ min ≍ σ√log d/√n, [56]'s stronger RIP requirement cannot be met by random design matrices such as Gaussian designs, so it cannot recover signals at that scale.The paper contrasts this with its own sufficient condition, which moves κ into a logarithmic factor.Source: H Comparing Our Results to [56] · p. 47
• The increasing-step-size scheme is computationally optimal up to poly-logarithmic factors, unlike the constant-step-size approach in [56].The paper reports that [56] does not yield a computationally optimal algorithm.Source: H Comparing Our Results to [56] · p. 47
• Theorem 2 tunes the step size using an estimate of w⋆ max computed from data, whereas [56] also depends on w⋆ min, an additional parameter that is unclear how to choose in practice.The paper describes its data-based tuning as avoiding an extra hyperparameter.Source: H Comparing Our Results to [56] · p. 48
• Both papers achieve minimax-optimal rates and dimension-independent rates when w⋆ min ≳ ∥Xᵀξ∥∞/n, but this paper's dimension-independent rate has an extra log k.The paper attributes the difference to stronger RIP assumptions in [56].Source: H Comparing Our Results to [56] · p. 48

I Table of Notation

• The notation defines w⋆ as the k-sparse true parameter with unknown support S, and uses ⊙ for coordinate-wise products and ≲ for inequalities up to absolute constants.It also uses w² for coordinate-wise squares.Source: I Table of Notation · p. 49