Safe screening

Overview

A safe screening rule is a test that runs before the solver. It names features whose coefficients must be exactly zero at the solution, so they can be removed from the problem. "Safe" means the test never removes a feature that would be non-zero. The solver then sees a smaller matrix and returns the same answer 1211.3966 §1.

The test rests on the dual form of the Lasso. The dual optimum lies in a polytope. If a region that contains the optimum at penalty λ can be found, and every point of that region gives a feature score below the boundary, the feature is inactive 1211.3966 §2.1. The region here is a ball, centred at the dual solution at a nearby penalty λ₀ that is already known. A tighter ball drops more features, and the variants differ only in how the ball is drawn 1211.3966 §2.2 1211.3966 §2.3.

The main use is a decreasing grid of penalties, where each solve can start from the previous answer. Each screen costs little next to the solve it saves 1211.3966 tab. 1 1211.3966 tab. 3. Only one paper is in the wiki so far, so the evidence comes from a single set of authors.

Variants

The ball test is shared by all variants (see Implementation sketch). They differ in the centre and radius of the ball.

Basic DPP from λmax

The ball is centred on y/λmax, the dual solution where every coefficient is zero. The radius is set by the gap between 1/λ and 1/λmax times ‖y‖. The test needs only X and y, so no earlier solve is required 1211.3966 §2.2 1211.3966 Cor. 4. Card: Lasso Screening Rules via Dual Polytope Projection.

Sequential DPP

The same ball, but centred on the dual solution at the previous grid point λ₀, computed from that solve's residual 1211.3966 Cor. 5. The paper's sequential form is the base for the later variants.

Ray-projection bound (Improvement 1)

The ball is shrunk by using the direction of the ray from the projection at λ₀. Only the part of the move to y/λ that is orthogonal to that ray is kept. The result lies inside the DPP ball 1211.3966 §2.3.1 [arxiv:1211.3966 Thm. 7] 1211.3966 Thm. 11.

Firm-nonexpansiveness bound (Improvement 2)

The ball is moved halfway towards y/λ and its radius is halved. It lies inside the DPP ball with half the radius 1211.3966 §2.3.2 1211.3966 Thm. 13 1211.3966 Thm. 14.

Enhanced DPP (EDPP)

Combines the two improvements. The centre is the ray-bound centre moved halfway, and the radius is half the orthogonal move. This is the smallest of the three balls in the paper's analysis, and the paper reports it as the fastest rule in its experiments 1211.3966 §2.3.3 [arxiv:1211.3966 Thm. 15] 1211.3966 tab. 1.

Group Lasso extension

The same ideas apply to whole groups of features. A group is dropped when the norm of its scores falls below a group-sized threshold. The paper gives the group rule and its sequential form, but its experiments compare only against a heuristic group rule 1211.3966 §3 1211.3966 Thm. 20 1211.3966 §4.2.

When to use which

Caveats for implementers.

• The rules are derived for ½‖y − Xβ‖² + λ‖β‖₁. Other scalings of λ must be mapped to this one first 1211.3966 eq. 1.
• The ball uses ‖y‖ and the column norms ‖x_i‖. The rules do not need the columns to be unit length 1211.3966 §2.1. The paper's basic-version comparison normalises the features, and the other experiments do not state whether they do 1211.3966 §4.1.1.
• The improvement terms divide by ‖v₁‖². At λ = λmax, v₁ is taken from the column most correlated with y. The paper does not cover ties or zero columns 1211.3966 Thm. 7 1211.3966 Lemma 10.
• For groups, ‖X_g‖ is the spectral norm, the largest singular value of the group's columns 1211.3966 §3.2.

Implementation sketch

The shared core is a ball test. Given a centre c and radius r, feature i is dropped when its score at the furthest point of the ball is below 1. The variants only change how c and r are built.

# Ball test, common to all variants  [arxiv:1211.3966 Thm. 3]
# Input:  X (N, p), centre c (N,), radius r (scalar)
# Output: keep (p,) bool, True = feature stays in the problem

def ball_screen(X, c, r):
    # Column norms ||x_i||, one per feature  [Thm. 3]
    col_norms = norm(X, axis=0)                     # (p,)
    # Largest |x_i^T theta| over the ball is |x_i^T c| + r ||x_i||
    # [Thm. 3 proof]
    worst = abs(X.T @ c) + r * col_norms           # (p,)
    # Drop feature i when even the worst case stays below 1  [Thm. 3]
    return worst >= 1                               # (p,) bool

# Ball for EDPP at penalty lam, from the solve at lam0  [Thm. 7] [Thm. 15]
# Input:  y (N,), lam, lam0, theta0 = theta*(lam0) (N,), v1 (N,)
# Output: centre (N,), radius (scalar)

def edpp_ball(y, lam, theta0, v1):
    # Move from theta0 towards y / lam  [Thm. 7]
    v2 = y / lam - theta0                           # (N,)
    # Keep only the part orthogonal to the ray direction  [Thm. 7]
    v2_perp = v2 - (v1 @ v2) / (v1 @ v1) * v1       # (N,)
    # Half-way centre, half-radius  [Thm. 15]
    centre = theta0 + 0.5 * v2_perp                 # (N,)
    radius = 0.5 * norm(v2_perp)
    return centre, radius

# Sequential path loop  [Cor. 17]
# lambdas: decreasing, all below lambda_max. solve() is any Lasso solver.

def screened_path(X, y, lambdas, solve):
    # lambda_max and the column that attains it  [eq. 7]
    corr = X.T @ y
    best = argmax(abs(corr))
    lam_max = abs(corr[best])
    # At lambda_max: theta* = y / lam_max, ray from the best column  [eq. 9]
    theta0 = y / lam_max
    v1 = sign(corr[best]) * X[:, best]
    kept_per_lambda = []
    for lam in lambdas:
        c, r = edpp_ball(y, lam, theta0, v1)
        keep = ball_screen(X, c, r)
        kept_per_lambda.append(keep)
        # Solve on the kept features only; then update the dual point
        beta = solve(X[:, keep], y, lam)
        theta0 = (y - X[:, keep] @ beta) / lam
        # Next ray direction is the residual direction at lam  [Thm. 7]
        v1 = y / lam - theta0
    return kept_per_lambda

Notes on where variants differ:

• Basic DPP: centre y / lam_max, radius |1/lam − 1/lam_max| ‖y‖ [Thm. 3].
• Sequential DPP: centre theta0, radius |1/lam − 1/lam0| ‖y‖ [Cor. 5].
• Improvement 1: centre theta0, radius ‖v2⊥‖ [Thm. 11].
• Improvement 2: centre theta0 + ½(1/lam − 1/lam0) y, radius ½|1/lam − 1/lam0|‖y‖ [Thm. 13].
• EDPP: the function above [Thm. 15].
• At lam_max the EDPP ray v1 is sign(xᵀy) x. Later steps use y/lam − theta0 [Thm. 7].

Contradictions

None known. One paper so far.

Within the paper, some numbers do not match each other. These are noted for checking, not resolved:

• The text gives a 30× speedup for Prostate. The table ratio is about 33× 1211.3966 §4.1.2 1211.3966 tab. 1.
• The Leukemia text says 52 samples, but the matrix is 55 × 11225 1211.3966 §4.1.2 1211.3966 tab. 3.
• The synthetic-2 generation method is described inconsistently. The text gives both i.i.d. Gaussian columns and pairwise correlation 0.5^{|i−j|} 1211.3966 §4.1.2.
• Whether the "+solver" times in the tables include screening time is not stated 1211.3966 tab. 1 1211.3966 tab. 3.

Papers

• Lasso Screening Rules via Dual Polytope Projection: introduces the DPP family and EDPP for the Lasso, and extends the idea to group Lasso. Author-reported speedups; the solver is SLEP.