L1 (Lasso) Selection

September 2026 · 16 min read

Lasso fits a linear model while charging a penalty for coefficient magnitude. Some fitted coefficients become exactly zero. Their columns then contribute nothing to that model's prediction, giving us a way to select features during training.

A zero coefficient does not prove a feature is useless. It reflects a particular model, training sample and penalty. We will watch three coefficients shrink, then look at what a zero can and cannot tell us.

Balance prediction error against coefficient size

For row $i$, predict the target with an intercept $b$ plus a weighted sum of its inputs: $\hat y_i=b+\sum_{j=1}^{P}x_{ij}w_j$. Here $P$ is the number of features and $w_j$ is the coefficient of feature $j$. A coefficient's sign sets the direction of its contribution.

With $N$ training rows, lasso chooses $b$ and $w$ to minimize two costs: half the mean squared prediction error, plus the sum of absolute coefficient values multiplied by a penalty strength $\alpha\geq0$:

The intercept is not penalized. A larger $\alpha$ makes coefficient size more expensive relative to fit. The factor $1/2$ simplifies differentiation; this convention matches scikit-learn's Lasso. Standardize continuous inputs using training statistics so their units do not arbitrarily change the penalty tradeoff.

Watch features leave a small model

Take the eight possible combinations of three inputs A, B and C, each equal to −1 or +1. Across these rows, each column has mean 0 and variance 1, and distinct columns have zero correlation. Set the target to $y=2A-0.8B+0.2C$, with no added noise. This deliberately simple dataset lets us see the mechanism exactly.

Without a penalty, the coefficients are $(2,-0.8,0.2)$ and predictions match every target. Increase the penalty below. B's negative coefficient is a real contribution, not a negative importance score.

Why does lasso reach zero instead of only getting smaller?

For these orthogonal columns, each coefficient can be solved independently. Call its unpenalized best value $z_j$; here those values are 2, −0.8 and 0.2. The solution moves each value toward zero by $\alpha$ and stops at zero rather than crossing it. This operation is called soft thresholding:

At $\alpha=0.5$, A moves from 2 to 1.5, B moves from −0.8 to −0.3, and C reaches zero. The absolute-value penalty has a corner at zero, so a whole interval of small unpenalized values maps to that same zero. An L2 penalty instead gives smooth shrinkage; the regularization chapter compares the geometry.

This direct thresholding of the initial coefficients relies on orthogonality. With correlated inputs, fit the coefficients together. Coordinate descent repeatedly updates one coefficient using the residual left after accounting for the others. A feature can enter or leave as those other coefficients change.

The smaller model wins the penalized objective

At $\alpha=0.5$, compare three candidate coefficient vectors on the same eight rows. The objective uses half the MSE, not the MSE itself:

A zero may simply mean another column took the role

Consider a separate task: one centered, unit-variance input $x$, target $y=2x$, and two identical copies of $x$ offered to the model. At $\alpha=0.5$, all three allocations below are optimal:

With near-duplicates, small data changes can change which column survives. Evaluate prediction quality and selection stability separately. A selected set is not a causal explanation. If correlated groups make selection unstable, compare elastic net or a domain-chosen representative rather than treating one run's support as definitive.

Turn coefficients into a usable input subset

The retained-column mask is often called the support. Read coefficient magnitudes, not their signs. Use a small, explicit numerical threshold so floating-point residue near zero does not masquerade as a meaningful effect. For example, the theoretical zero at $\alpha=0.2$ can appear as a value around $10^{-16}$.

SelectFromModel retains magnitudes greater than or equal to its threshold. A threshold of zero therefore also retains exact zeros. The code uses 1e-6, separate from the lasso penalty, to identify the nonzero effects in this toy. On real data, document both the input scaling and the cutoff.

import numpy as np
from sklearn.linear_model import Lasso, LinearRegression, Ridge
from sklearn.feature_selection import SelectFromModel
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# All eight combinations: columns are centered, unit-variance and orthogonal.
X = np.array([[-1,-1,-1], [-1,-1,1], [-1,1,-1], [-1,1,1],
              [1,-1,-1],  [1,-1,1],  [1,1,-1],  [1,1,1]], dtype=float)
y = X @ np.array([2, -.8, .2])
z = X.T @ y / len(X)
for alpha in [0, .1, .2, .5, .8, 1, 2]:
    model = (Lasso(alpha=alpha, tol=1e-12, max_iter=10000)
             if alpha else LinearRegression()).fit(X, y)
    expected = np.sign(z) * np.maximum(np.abs(z) - alpha, 0)
    np.testing.assert_allclose(model.coef_, expected, atol=1e-10)
    print(alpha, np.round(model.coef_, 4))

# Select columns, then train a separate downstream predictor.
pipeline = Pipeline([
    ("scale", StandardScaler()),
    ("select", SelectFromModel(
        Lasso(alpha=.5, tol=1e-12), threshold=1e-6)),
    ("predict", Ridge(alpha=1)),
]).fit(X, y)
print(pipeline.named_steps["select"].get_support())  # [True, True, False]
# This eight-row toy checks the mechanics; it is not a performance evaluation.

The pipeline first fits scaling, then lasso selection, then a ridge predictor on the surviving columns. The downstream fit has its own coefficients; it does not reuse the lasso weights. If you keep lasso itself as the predictor, no second model is required. If a penalty removes every column, handle that case explicitly, for example by comparing an intercept-only baseline.

Choose the penalty inside validation

Choose $\alpha$ by validating the complete procedure against its intended metric and feature cost. A sparse model with poor predictions is not a success. Compare with an appropriate model using all eligible inputs, and inspect whether the chosen subset is stable across training folds.

Pass the entire preprocessing–selection–prediction pipeline to the tuning procedure so each validation fold gets a fresh fit. Reserve an untouched evaluation set, or use nested cross-validation, to assess the selected procedure. Respect group and time boundaries where needed. A scaler fitted before an internally cross-validating estimator has already seen that estimator's inner validation rows.

Optional implementation: update one coefficient at a time

The C++ function accepts a dense matrix and target already centered using training means. For each column, it temporarily adds back that column's current contribution to the residual, soft-thresholds its association with the residual, and updates the predictions. The denominator accounts for the column's mean squared magnitude, so the solver also handles columns with different scales.

#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <vector>
using Matrix = std::vector<std::vector<double>>;

double soft_threshold(double z, double alpha) {
    return std::copysign(std::max(std::abs(z) - alpha, 0.0), z);
}

// Dense training columns and target must already be centered.
// Standardize columns first when comparable penalties are intended.
std::vector<double> fit_lasso(const Matrix& X, const std::vector<double>& y,
                              double alpha) {
    if (X.empty() || X[0].empty() || y.size() != X.size() ||
        !std::isfinite(alpha) || alpha <= 0)
        throw std::invalid_argument("invalid shape or positive penalty");
    const auto n = X.size(), p = X[0].size();
    std::vector<double> scale(p, 0), w(p, 0), residual = y;
    for (const auto& row : X) {
        if (row.size() != p) throw std::invalid_argument("ragged matrix");
        for (std::size_t j = 0; j < p; ++j) {
            if (!std::isfinite(row[j])) throw std::invalid_argument("finite X required");
            scale[j] += row[j] * row[j] / n;
        }
    }
    for (double value : y) if (!std::isfinite(value))
        throw std::invalid_argument("finite y required");
    for (int sweep = 0; sweep < 10000; ++sweep) {
        double change = 0;
        for (std::size_t j = 0; j < p; ++j) {
            double z = 0;
            for (std::size_t i = 0; i < n; ++i)
                z += X[i][j] * (residual[i] + X[i][j] * w[j]) / n;
            double next = scale[j] == 0 ? 0 : soft_threshold(z, alpha) / scale[j];
            if (!std::isfinite(next)) throw std::runtime_error("nonfinite update");
            double delta = next - w[j];
            for (std::size_t i = 0; i < n; ++i) residual[i] -= X[i][j] * delta;
            change = std::max(change, std::abs(delta));
            w[j] = next;
        }
        if (change < 1e-10) return w;
    }
    throw std::runtime_error("coordinate updates did not converge");
}

References and further reading

• Tibshirani: Regression Shrinkage and Selection via the Lasso — the original formulation and its sparse linear-model interpretation.
• Scikit-learn: Lasso — objective convention, coefficient fitting and convergence controls.
• Scikit-learn: SelectFromModel — coefficient magnitudes, inclusive thresholds and support masks.
• Friedman, Hastie and Tibshirani: Regularization Paths for Generalized Linear Models via Coordinate Descent — coordinate updates and regularization paths.