Lasso Screening Rules via Dual Polytope Projection

In plain words

Lasso is a common way to fit a model that uses only a few input features: it pushes the weights of unhelpful features to exactly zero. Solving it is slow when there are very many features, and it has to be solved again for each penalty setting you try. A screening rule is a quick test that says, before solving, which features are certain to end up at zero, so they can be dropped from the problem.

The rules here come from the geometry of the problem's dual form. The right answer must lie inside a small region that can be computed from the solution at a nearby penalty setting. The enhanced version draws a tighter region, so it can drop more features. Because the test never drops a feature that would have been used, the smaller problem gives the same answer. It is a pre-step, so any Lasso solver can run after it. The authors report large speedups on real data sets, and they also extend the idea to group Lasso.

TL;DR

• Safe screening for the Lasso: DPP bounds the dual optimum in a ball around the dual solution at a known penalty, then discards features whose dual scores fall inside a threshold §2.2 [Thm. 3].
• EDPP (enhanced DPP) shrinks that ball by combining two facts about projections onto a convex polytope: a ray-projection bound and firm nonexpansiveness §2.3 [Thm. 7] [Thm. 13] [Thm. 15].
• The sequential form uses the solution at the previous penalty on a decreasing grid and screens the next one [Cor. 5] [Cor. 17]. The authors report that it works with any Lasso solver §1.
• A group Lasso extension screens whole groups with the same ideas §3 [Thm. 20] [Cor. 21]. All performance numbers are author-reported.

Problem

Solving the Lasso is hard when the feature count and the sample count are both large, possibly too large to load the data matrix into memory §1. Screening aims to find the features whose coefficients are zero at the solution and remove them first §1. Heuristic screens such as strong rules can wrongly discard active features and must check KKT conditions afterwards; safe screens such as SAFE and DOME guarantee that discarded features are inactive §1. Safe screens are only as strong as their estimate of the dual optimum: a more accurate estimate discards more features §1 §2.1.

Model selection needs the Lasso solved along a grid λ₁ > λ₂ > … > λ_K, which is slow without screening §1 §2.2. The authors build sequential rules that reuse the previous solution on that grid §2.2 §1.

Method

Setup §2.1. The Lasso is ½‖y − Xβ‖₂² + λ‖β‖₁ with X ∈ ℝ^{N×p} eq. 1. No unit-length assumption on y or the columns of X is made §2.1. Its dual is a projection problem: maximise ½‖y‖² − (λ²/2)‖θ − y/λ‖² subject to |x_iᵀθ| ≤ 1 for all i eq. 2. The KKT conditions are y = Xβ*(λ) + λθ*(λ) eq. 3 and x_iᵀθ*(λ) equal to sign(β_i*) when β_i* ≠ 0, otherwise in [−1, 1] eq. 4. So if |x_iᵀθ*(λ)| < 1 then β_i*(λ) = 0 §2.1, unnumbered rule R1. Since θ* is unknown, the rule is relaxed: if a region Θ that contains θ*(λ) has sup over θ ∈ Θ of |x_iᵀθ| < 1, then β_i*(λ) = 0 §2.1, unnumbered rule R1′. Smaller Θ means more features discarded §2.1.

Dual as projection §2.1. The dual optimum is the projection θ*(λ) = P_F(y/λ) onto the feasible polytope F, the intersection of 2p half-spaces eq. 6. The threshold λmax = max_i |x_iᵀy| eq. 7 is the smallest penalty with β* = 0 §2.1. For λ ≥ λmax, θ*(λ) = y/λ eq. 9.

DPP §2.2.

1. Nonexpansiveness of projections gives ‖P_C(w₂) − P_C(w₁)‖ ≤ ‖w₂ − w₁‖ [Thm. 1] eq. 10. Applied to θ*(λ) and θ*(λ₀), this gives ‖θ*(λ) − θ*(λ₀)‖ ≤ |1/λ − 1/λ₀| ‖y‖ [Thm. 2] eq. 11. So θ*(λ) lies in the ball B(θ*(λ₀), |1/λ − 1/λ₀| ‖y‖) eq. 12.
2. Maximising |x_iᵀθ| over that ball gives |x_iᵀθ*(λ₀)| + ρ‖x_i‖ with ρ the radius [Thm. 3 proof] eq. 14. Therefore β_i*(λ) = 0 if |x_iᵀθ*(λ₀)| < 1 − ‖x_i‖ ‖y‖ |1/λ − 1/λ₀| [Thm. 3] eq. 13.
3. Basic DPP takes λ₀ = λmax, where θ* = y/λmax [Cor. 4]. Sequential DPP takes the previous grid point as λ₀ and uses θ*(λ_k) = (y − Xβ*(λ_k))/λ_k [Cor. 5].
4. The paper contrasts DPP with SAFE, which discards when |x_iᵀy| < λ − ‖x_i‖‖y‖(λmax − λ)/λmax eq. 15. Both have the same radius, but DPP's centre is y/λmax, not y/λ §2.2 remark.

Improvement 1: ray projection §2.3.1.

1. Lemma 6: for w outside a closed convex set C, the projection of the whole ray P_C(w) + t(w − P_C(w)), t ≥ 0, is P_C(w) [Lemma 6].
2. Define v₁(λ₀) = y/λ₀ − θ*(λ₀) if λ₀ < λmax, and v₁ = sign(xᵀy) x at λ₀ = λmax, where x* = argmax_i |x_iᵀy| [Thm. 7]. Define v₂(λ, λ₀) = y/λ − θ*(λ₀), and v₂⊥ = v₂ − (⟨v₁, v₂⟩/‖v₁‖²) v₁ [Thm. 7].
3. Theorem 7: θ*(λ) ∈ B(θ*(λ₀), ‖v₂⊥‖), and this ball lies inside the DPP ball. The proof needs ⟨v₁, v₂⟩ ≥ 0 [Thm. 7 proof] and uses Lemma 10 when λ₀ = λmax [Lemma 10].
4. Rule (Improvement 1): β_i*(λ) = 0 if |x_iᵀθ*(λ₀)| < 1 − ‖v₂⊥(λ, λ₀)‖ ‖x_i‖ [Thm. 11]. The sequential form replaces θ*(λ₀) with (y − Xβ*(λ_k))/λ_k and λ₀ with λ_k §2.3.1, displayed rule.

Improvement 2: firm nonexpansiveness §2.3.2.

1. Projections are firmly nonexpansive [Thm. 12]. Applied to θ*(λ) and θ*(λ₀) this gives θ*(λ) ∈ B(θ*(λ₀) + ½(1/λ − 1/λ₀) y, ½|1/λ − 1/λ₀| ‖y‖) [Thm. 13]. This ball is inside the DPP ball and has half its radius [Thm. 13].
2. Rule (Improvement 2): β_i*(λ) = 0 if |x_iᵀ(θ*(λ₀) + ½(1/λ − 1/λ₀) y)| < 1 − ½(1/λ − 1/λ₀) ‖y‖ ‖x_i‖ [Thm. 14]. The sequential form is displayed in §2.3.2.

EDPP §2.3.3.

1. Combining the two ideas gives ‖θ*(λ) − (θ*(λ₀) + ½v₂⊥(λ, λ₀))‖ ≤ ½‖v₂⊥(λ, λ₀)‖ [Thm. 15]. This is the smallest of the three balls §2.3.3.
2. Rule: β_i*(λ) = 0 for λ ∈ (0, λ₀] if |x_iᵀ(θ*(λ₀) + ½v₂⊥(λ, λ₀))| < 1 − ½‖v₂⊥(λ, λ₀)‖ ‖x_i‖ [Thm. 16].
3. Sequential EDPP [Cor. 17] uses θ*(λ_k) = (y − Xβ*(λ_k))/λ_k and screens the next grid point λ_{k+1}. The pseudocode below implements this.

# Sequential EDPP screening along a decreasing lambda path  [§2.3.3] [Cor. 17]
# Input:  X (N, p), y (N,), lambdas (decreasing, all below lambda_max)
# Output: kept[k] = column indices to solve with at lambdas[k]

def edpp_path(X, y, lambdas):
    # lambda_max = max_i |x_i^T y|, attained at column best  [eq. 7]
    correlations = X.T @ y                                    # (p,)
    best = argmax(abs(correlations))
    lambda_max = abs(correlations[best])
    # Path starts at lambda_max, where beta* = 0 and theta* = y / lambda_max
    path = [lambda_max] + lambdas
    theta_prev = y / lambda_max                               # (N,)  [eq. 9]
    # Sign-adjusted best column gives the ray direction at lambda_max  [Thm. 7]
    v1_at_max = sign(correlations[best]) * X[:, best]         # (N,)
    kept = []
    for k in range(len(path) - 1):
        lam_prev = path[k]
        lam = path[k + 1]
        # Ray direction from theta*(lam_prev)  [Thm. 7]
        if k == 0:
            v1 = v1_at_max
        else:
            v1 = y / lam_prev - theta_prev
        # Vector from theta*(lam_prev) towards y / lam  [Thm. 7]
        v2 = y / lam - theta_prev
        # Keep only the part of v2 orthogonal to v1  [Thm. 7]
        v2_perp = v2 - (v1 @ v2) / (v1 @ v1) * v1
        # Ball centred at theta_prev + v2_perp / 2 with half the radius  [Thm. 15]
        centre = theta_prev + 0.5 * v2_perp
        radius = 0.5 * norm(v2_perp)
        # Discard column i if |x_i^T centre| < 1 - radius * ||x_i||  [Cor. 17]
        scores = abs(X.T @ centre)                            # (p,)
        threshold = 1 - radius * column_norms(X)              # (p,)
        keep = scores >= threshold                            # (p,) bool
        kept.append(nonzero(keep))
        # Solve on the kept columns, then recover theta*(lam)  [eq. 3]
        beta = lasso_solve(X[:, keep], y, lam)
        theta_prev = (y - X[:, keep] @ beta) / lam            # (N,)
    return kept

Group Lasso §3. The group problem is ½‖y − Σ_g X_g β_g‖² + λ Σ_g √n_g ‖β_g‖₂ §3.1. Its dual has constraints ‖X_gᵀθ‖₂ ≤ √n_g §3.1. The KKT conditions give ‖X_gᵀθ*(λ)‖₂ < √n_g ⇒ β_g*(λ) = 0 §3.1, rule R2. The group penalty value is λ̄max = max_g ‖X_gᵀy‖₂ / √n_g, and θ*(λ) = y/λ for λ ≥ λ̄max §3.1.

1. Lemma 18: at λ₀ < λ̄max, v̄₁ = y/λ₀ − θ*(λ₀); at λ₀ = λ̄max, v̄₁ = X_* X_ᵀ y, where X_ is the group with largest ‖X_gᵀy‖₂/√n_g [Lemma 18].
2. Theorem 19: θ*(λ) ∈ B(θ*(λ₀) + ½v̄₂⊥, ½‖v̄₂⊥‖), with v̄₂⊥ defined as for the Lasso but using v̄₁ [Thm. 19].
3. Rule: β_g*(λ) = 0 if ‖X_gᵀ(θ*(λ₀) + ½v̄₂⊥)‖₂ < √n_g − ½‖v̄₂⊥‖ ‖X_g‖₂ [Thm. 20]. Sequential form: [Cor. 21].

Hyperparameters

The paper's method has no tuned hyperparameters. The table lists the experimental settings that the results depend on.

Data sets. Sizes are the data-matrix dimensions given in the paper.

Results

Author-reported. Setup: Lasso solved along a 100-point path, solver SLEP, times in seconds for the whole path. Screening times are listed separately. The paper's text gives speedups, and the card's ratios below are computed from the table values.

Family of DPP rules, sequential path tab. 1 fig. 1. Prostate, PIE and MNIST:

The paper says EDPP gives about 150× (PIE) and 230× (MNIST). The table ratios are 152× and 231×. For Prostate the ratio is 33×; the paper says 30× for Prostate in §4.1.2. The paper also says Improvement 1 reaches about 20×, 60× and 70×, which matches the table ratios of 19×, 56× and 68×.

Basic versions, rejection ratio only fig. 2. Values are in the figure, not in a table. The paper reports that EDPP beats SAFE, DOME and the strong rule on five of six data sets, and that on PIE EDPP and DOME are similar. Both beat SAFE and the strong rule there.

Synthetic, sequential path tab. 2 fig. 3. Solver and screening times (seconds):

The paper reports that strong rule and EDPP have comparable rejection ratios, and that the strong rule's time is about twice EDPP's §4.1.2 fig. 3. The table's screen times give 1.8×.

Real data, sequential path tab. 3 fig. 4. Times in seconds:

The paper reports EDPP speedups of about 10× (Breast), 17× (Leukemia) and 30× (Prostate), and about two orders of magnitude on PIE, MNIST and SVHN. The table gives 10×, 17×, 33×, 152×, 231× and 185×. For SVHN the paper says the solver takes about 3 hours and EDPP's path takes under 1 minute §4.1.2.

EDPP with LARS tab. 4. Times in seconds:

The speedup curves are in fig. 5.

Group Lasso tab. 5 fig. 6. X is 250 × 200000. Times in seconds for the group path:

The paper says the solver speeds up by about 80×, 120× and 160× with EDPP. The table gives 84×, 120× and 163×. The paper reports that EDPP discards more groups than the strong rule and is more robust to n_g §4.2 fig. 6. The paper does not compare with any safe group rule, because it knows of none §4 §4.2.

Paper vs code

No official code.

Searched: the paper and its arXiv metadata have no code or project link. A GitHub repository search for "edpp", "lasso screening", "dual polytope projection" and "jieping lasso" found no repository owned by an author. Third-party repositories were not used.

Pitfalls

• Scaling. The rules are derived for the unnormalised objective ½‖y − Xβ‖² + λ‖β‖₁ eq. 1. The thresholds use ‖y‖, ‖x_i‖ and 1/λ. A solver with a different loss scaling or λ convention must be mapped to this λ first, or the screens are wrong.
• Exact previous solution assumed. The sequential rules use β*(λ_k) as the exact solution [Cor. 5] [Cor. 17]. The theorems take θ*(λ₀) as exact [Thm. 3] [Thm. 7]. The paper does not analyse approximate solutions, solver tolerance or duality gap. The tolerance it used is not stated.
• Degenerate directions. v₂⊥ divides by ‖v₁‖² [Thm. 7] [Thm. 19]. v₁ is nonzero for λ₀ < λmax, and at λmax it is nonzero whenever x* ≠ 0 [Lemma 10]. The paper does not discuss zero columns, ties in argmax, or sign(x*ᵀy) = 0.
• Group norm. The group rule uses ‖X_g‖ as an operator norm, because the proof bounds ‖X_gᵀu‖₂ ≤ ‖X_g‖₂‖u‖₂ §3.2, proof of Thm. 20. Use the spectral norm, the largest singular value. The LaTeX source also has a commented-out draft that used the Frobenius norm. That draft is not in the compiled paper. A Frobenius norm is also safe, since it is larger and so discards fewer groups.
• Typo in the text. §2.3.3 says the second speedup is "roughly 10 times" for "Improvement 1". It means Improvement 2. The table gives Improvement 2 about 7×, 11× and 11× and Improvement 1 about 19×, 56× and 68×.
• Leukemia size. The text says 52 samples, but the matrix is 55 × 11225 §4.1.2 fig. 4. The paper does not resolve this.
• Synthetic 2 generation. The text says the entries are i.i.d. Gaussian and also that the pairwise correlation is 0.5^{|i−j|}, which contradicts itself. The text also has "corr(x_i, x_i) = 0", which should be i ≠ j. The method for building correlated columns is not stated.
• Rejection ratio. It is defined against "the actual number of zero features in the ground truth" §2.3.3. The paper does not say how that zero set is found, for example what solver tolerance decides a zero.
• Time accounting. The +solver columns may or may not include screening time. The captions say the last columns are screening totals, but do not say whether the +solver columns include them. Whether the strong rule's time includes its KKT re-check loop is also not stated.
• Preprocessing. Feature normalisation is stated only for the basic-version comparison §4.1.1. It is not stated for the real-data path comparison §4.1.2 or the group experiments §4.2.
• Missing baseline formulas. The paper writes only the basic SAFE test eq. 15. It does not write the sequential SAFE test or the strong-rule test. Those must come from the works it cites.
• Hardware and solvers. The paper does not give hardware, implementation language, solver tolerances, warm starts or random seeds. Speedups depend on the machine and cannot be reproduced exactly from the text.
• Group comparison. The group experiments compare only with the strong rule, which is heuristic. The paper gives no safe group baseline §4.
• Figures. Rejection ratios appear only in fig. 1 to fig. 4 and fig. 6, not in tables. Read them from the figures; the card describes them in words.

Lineage

• Builds on: SAFE (El Ghaoui et al. 2012) and DOME (Xiang et al. 2011, 2012) for safe screening, and the strong rules (Tibshirani et al. 2012) for sequential screening §1. Their arXiv ids are not in this wiki, so none are linked.
• Compared against: basic and sequential SAFE, DOME, strong rule, and EDPP with LARS §4.1 §4.1.1 §4.1.2.