LASSO Regression

Ridge regression tames unstable coefficients by gently pulling all of them towards zero. It works, but it never lets go of any variable: a model built on 1000 wavelengths still has 1000 non-zero coefficients. Chemists often want to ask a sharper question:

Which wavelengths actually carry the information about my analyte?

The lasso answers it with a surprisingly small change to ridge regression. Replace the sum of squared coefficients in the penalty by the sum of their absolute values, and some coefficients stop shrinking gradually: they become exactly zero. The model selects variables while it fits.

In 1996, Robert Tibshirani, then at the University of Toronto, published Regression Shrinkage and Selection via the Lasso in the Journal of the Royal Statistical Society [1]. The name is an acronym, Least Absolute Shrinkage and Selection Operator, and it describes the method well: it shrinks coefficients with an absolute-value penalty and, in doing so, selects variables.

Tibshirani had trained at Stanford under Bradley Efron, and his work belongs to a period in which statisticians were increasingly facing datasets with more variables than observations. The ingredients were not entirely new: ridge regression had shown in 1970 that penalizing coefficients stabilizes ill-conditioned problems [3], and geophysicists had already used L1 penalties to obtain sparse solutions in seismic deconvolution [2]. Frank and Friedman had proposed bridge regression, a family of penalties on ${\sum }_j|{\beta }_j{|}^q$ that includes both ridge ($q=2$) and the absolute-value penalty ($q=1$) [14]. Most directly, Breiman's non-negative garrote (1995), which Tibshirani took as his starting point, had already shown that a constrained least squares fit could shrink some coefficients to exactly zero [13]. Tibshirani's contribution was to formulate the idea cleanly for regression, explain why it produces exact zeros, and show how well it predicts.

The lasso became truly practical a few years later, with the LARS algorithm, which computes the whole solution path at the cost of roughly one least squares fit in favourable cases [4], and with fast coordinate descent methods, which are what most software uses today [5].

Consider a near-infrared calibration: 200 tablets, spectra recorded at 1000 wavelengths, and the content of the active ingredient as the response. The data matrix $\mathbf{X}$ has 200 rows and 1000 columns.

Least squares cannot even produce a unique answer: with more unknowns than equations, infinitely many coefficient vectors fit the calibration data perfectly. Ridge regression, PCR and PLS all solve that problem, but in ways that keep every wavelength in the model: ridge with small coefficients, PCR and PLS through components that are combinations of all wavelengths. None of them hands you a short list such as "these 15 wavelengths are the ones that matter".

A short list is valuable. It can be checked against known absorption bands, it is easier to explain, and it tells an instrument designer which few wavelengths a cheap dedicated sensor would need. This is where the lasso enters.

The lasso minimizes the usual squared error plus a penalty on the absolute values of the coefficients:

• The first term, the data fidelity, is the mean squared residual (with a convenient factor of one half): it wants the model to fit.
• The second term, the L1 penalty, adds up the absolute values of the coefficients: it wants them small, and as we will see, it is happy to make them exactly zero.

The tuning parameter $\lambda \ge 0$ sets the balance. With $\lambda =0$ you are back to least squares. As $\lambda$ grows, coefficients shrink and drop out one after another, until at a large enough $\lambda$ every coefficient is zero and the model predicts the mean of $\mathbf{y}$. As in ridge regression, the intercept is not penalized, and the predictors are standardized first so that the penalty treats them fairly.

The smallest $\lambda$ for which every coefficient is zero has a simple formula. With centred $\mathbf{y}$ and standardized columns ${\mathbf{x}}_j$, it is

the largest covariance between a predictor and the response: below it, the data pull on that predictor harder than the penalty holds it at zero. Software starts the path there and works downwards.

This is the heart of the method, and there are two complementary ways to see it.

Just as for ridge regression, the penalized problem is equivalent to a constrained one:

With two coefficients, the ridge budget ${\beta }_1^2+{\beta }_2^2\le t^2$ is a disc, but the lasso budget $|{\beta }_1|+|{\beta }_2|\le t$ is a diamond whose corners lie on the axes. The contours of the squared error are ellipses around the least squares estimate, and the solution is the first point where an ellipse touches the allowed region.

Drag the least squares point around (or move it with the sliders). Over a large part of the plane, the first contact is a corner of the diamond, where one coefficient is exactly zero. A disc has no corners, so with ridge the contact point essentially never sits on an axis: tick "Show both shapes" to see the two solutions for the same budget side by side. In higher dimensions the L1 ball has corners, edges and faces lying in coordinate planes, so many coefficients can be zero at once.

The second explanation looks at the penalty one coefficient at a time. The ridge penalty ${\beta }_j^2$ has slope $2{\beta }_j$: its pull towards zero gets weaker and weaker as the coefficient gets smaller, so it never quite reaches zero. The lasso penalty $|{\beta }_j|$ has slope $\pm 1$ everywhere except at zero: it pulls with the same force no matter how small the coefficient already is. If the data push a coefficient away from zero less strongly than $\lambda$ pulls it back, the coefficient stays exactly at zero.

The cleanest case is when the predictors are uncorrelated and standardized, so that ${\mathbf{X}}^T\mathbf{X}=n\mathbf{I}$. Then each coefficient can be computed separately from its least squares value ${\hat{\beta }}_j$:

(for ridge written with the same $1/(2n)$ scaling, $\frac{1}{2n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }^2+\frac{\lambda }{2}\Vert \mathbf{\beta }{\Vert }^2$). Ridge divides every coefficient by the same factor, so nothing becomes zero. The lasso subtracts $\lambda$ from the size of every coefficient and sets to zero those that would cross it. This operation is called soft thresholding: small coefficients are removed, large ones survive but are shifted towards zero. (Hard thresholding, which keeps the large coefficients untouched, corresponds, in this orthogonal case, to best subset selection.)

Because of the absolute values, the lasso has no closed-form solution in general. The most common algorithm, coordinate descent, simply cycles through the coefficients, updating each one by a soft-thresholding step while the others are held fixed, until nothing changes [5]. For standardized predictors the update of coefficient $j$ is:

where ${\mathbf{r}}_{(j)}$ is the residual computed without the contribution of predictor $j$. Like gradient descent, it is an iterative method: start somewhere, improve, repeat. The demos on this page use exactly this algorithm.

As $\lambda$ decreases from a large value to zero, the lasso solution traces a regularization path:

Between the points where variables enter or leave, each coefficient changes linearly with $\lambda$: the lasso path is piecewise linear, a property that the LARS algorithm exploits [4]. On a logarithmic $\lambda$ axis, which is how paths are usually drawn, those straight pieces look curved; the demo lets you switch to a linear axis to see them. Explore the path for a simulated NIR calibration with ten wavelengths, and switch to ridge to compare.

A few things to notice. The analyte band (1730 nm, then 1715 nm) enters first. Some wavelengths with real but redundant information enter later and are pushed out again when a correlated neighbour takes over. At the cross-validation optimum the model uses only five of the ten wavelengths, with a prediction error (about 0.91) close to that of ridge (about 0.84), which keeps all ten. Note that the two methods use different scales for $\lambda$ (the lasso objective has the $1/(2n)$ factor, the ridge objective does not), so compare shapes and errors, not $\lambda$ values.

As for ridge regression, $\lambda$ is chosen by cross-validation:

Software usually reports two choices:

• ${\lambda }_{\min }$, the value that minimizes the cross-validated error;
• ${\lambda }_{1\text{SE}}$, the largest $\lambda$ whose error is within one standard error of the minimum. It gives a sparser model that predicts almost as well, and is often preferred.

The last row deserves a sentence. Ridge regression is the most probable solution under a Gaussian prior on the coefficients; the lasso is the most probable solution (the posterior mode) under a Laplace prior, which has a sharp peak at zero and heavier tails [1],[12]. The sharp peak is the probabilistic version of the diamond's corners.

Ridge tends to win when many variables carry small pieces of real information, which is often the case for broad spectral bands, and when groups of variables are strongly correlated. Faced with a group of nearly identical wavelengths, the lasso picks one more or less arbitrarily, and a slightly different dataset may make it pick another.

The lasso tends to win when only a few variables truly matter, when a compact and interpretable model is the goal, or when you need to know which wavelengths to measure.

When you want sparsity but your variables come in correlated groups, the elastic net mixes the two penalties [6]:

The mixing parameter $\alpha \in [0,1]$ moves between ridge ($\alpha =0$) and lasso ($\alpha =1$). Intermediate values keep the ability to produce zeros while encouraging correlated variables to enter or leave the model together, the so-called grouping effect.

The examples simulate a dataset with 500 correlated "wavelengths" of which only 10 are truly relevant, fit the lasso with cross-validation, and compare the selected variables with the true ones.

import numpy as np
from sklearn.linear_model import LassoCV, lasso_path
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
import matplotlib.pyplot as plt

# Simulate data: 200 samples, 500 correlated variables, only 10 truly relevant
rng = np.random.default_rng(42)
n_samples, n_wavelengths = 200, 500

X = rng.standard_normal((n_samples, n_wavelengths))
for j in range(1, n_wavelengths):                 # correlate neighbouring variables
    X[:, j] = 0.7 * X[:, j - 1] + 0.3 * X[:, j]

true_beta = np.zeros(n_wavelengths)
relevant_idx = [45, 102, 150, 203, 255, 310, 348, 400, 425, 470]
true_beta[relevant_idx] = rng.standard_normal(len(relevant_idx)) * 2
y = X @ true_beta + rng.standard_normal(n_samples) * 0.5

# --- Lasso with 10-fold cross-validation ---
# scikit-learn's Lasso minimizes 1/(2n) ||y - Xb||^2 + alpha ||b||_1, so alpha is our lambda.
# The scaler sits inside the pipeline, so each fold is standardized with its own training data.
model = make_pipeline(StandardScaler(),
                      LassoCV(cv=10, max_iter=10000, random_state=42)).fit(X, y)
lasso_cv = model[-1]

# lambda_1SE: the largest lambda whose CV error is within one standard error of the minimum
mse_mean = lasso_cv.mse_path_.mean(axis=1)
mse_se = lasso_cv.mse_path_.std(axis=1) / np.sqrt(lasso_cv.mse_path_.shape[1])
i_min = mse_mean.argmin()
lambda_1se = lasso_cv.alphas_[mse_mean <= mse_mean[i_min] + mse_se[i_min]].max()

selected = np.flatnonzero(lasso_cv.coef_)
print(f"lambda_min: {lasso_cv.alpha_:.4f} ({selected.size} non-zero of {n_wavelengths})")
print(f"Selected variables: {selected}")
print(f"True relevant variables: {relevant_idx}")

# --- Regularization path (on the whole standardized data set) ---
X_scaled = StandardScaler().fit_transform(X)
alphas, coefs, _ = lasso_path(X_scaled, y, alphas=np.logspace(-4, 1, 100))
coef_1se = lasso_path(X_scaled, y, alphas=[lambda_1se])[1][:, 0]
print(f"lambda_1SE: {lambda_1se:.4f} ({np.count_nonzero(coef_1se)} non-zero)")

plt.figure(figsize=(10, 5))
plt.semilogx(alphas, coefs.T, alpha=0.5)
plt.axvline(lasso_cv.alpha_, color="black", linestyle="--", label="lambda_min")
plt.axvline(lambda_1se, color="gray", linestyle=":", label="lambda_1SE")
plt.xlabel("lambda")
plt.ylabel("Coefficient")
plt.title("Lasso regularization path")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# --- Cross-validation curve ---
mse_std = lasso_cv.mse_path_.std(axis=1)
plt.figure(figsize=(8, 5))
plt.semilogx(lasso_cv.alphas_, mse_mean, "b-", linewidth=2, label="Mean CV error")
plt.fill_between(lasso_cv.alphas_, mse_mean - mse_std, mse_mean + mse_std, alpha=0.2)
plt.axvline(lasso_cv.alpha_, color="red", linestyle="--", label="lambda_min")
plt.axvline(lambda_1se, color="gray", linestyle=":", label="lambda_1SE")
plt.xlabel("lambda")
plt.ylabel("Mean squared error")
plt.title("Lasso cross-validation curve")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# --- True versus estimated coefficients ---
fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
axes[0].stem(true_beta, markerfmt="ro", linefmt="r-", basefmt="gray")
axes[0].set_title("True coefficients (10 non-zero)")
axes[1].stem(lasso_cv.coef_, markerfmt="bo", linefmt="b-", basefmt="gray")
axes[1].set_title(f"Lasso estimates at lambda_min ({selected.size} non-zero)")
axes[1].set_xlabel("Variable index")
for ax in axes:
    ax.set_ylabel("Coefficient")
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

% Lasso regression with cross-validation
% Requires the Statistics and Machine Learning Toolbox
rng(42);

% Simulate data: 200 samples, 500 correlated variables, only 10 truly relevant
n_samples = 200;
n_wavelengths = 500;
X = randn(n_samples, n_wavelengths);
for j = 2:n_wavelengths
    X(:, j) = 0.7 * X(:, j - 1) + 0.3 * X(:, j);
end

true_beta = zeros(n_wavelengths, 1);
relevant_idx = [46, 103, 151, 204, 256, 311, 349, 401, 426, 471];
true_beta(relevant_idx) = 2 * randn(numel(relevant_idx), 1);
y = X * true_beta + 0.5 * randn(n_samples, 1);

% lasso() standardizes the predictors by default and returns the whole path
[B, FitInfo] = lasso(X, y, 'CV', 10, 'NumLambda', 100);

beta_min = B(:, FitInfo.IndexMinMSE);
beta_1se = B(:, FitInfo.Index1SE);
fprintf('Lambda_min: %.4f (%d non-zero)\n', FitInfo.LambdaMinMSE, nnz(beta_min));
fprintf('Lambda_1SE: %.4f (%d non-zero)\n', FitInfo.Lambda1SE, nnz(beta_1se));

% Cross-validation curve and regularization path
lassoPlot(B, FitInfo, 'PlotType', 'CV');
lassoPlot(B, FitInfo, 'PlotType', 'Lambda', 'XScale', 'log');

% True versus estimated coefficients
figure;
subplot(2, 1, 1);
stem(true_beta, 'r', 'MarkerSize', 4); ylabel('Coefficient');
title('True coefficients'); grid on;
subplot(2, 1, 2);
stem(beta_min, 'b', 'MarkerSize', 4);
xlabel('Variable index'); ylabel('Coefficient');
title(sprintf('Lasso estimates (%d non-zero)', nnz(beta_min))); grid on;

fprintf('Selected: %s\n', mat2str(find(beta_min ~= 0)'));
fprintf('True relevant: %s\n', mat2str(relevant_idx));

library(glmnet)
set.seed(42)

# Simulate data: 200 samples, 500 correlated variables, only 10 truly relevant
n_samples <- 200
n_wavelengths <- 500
X <- matrix(rnorm(n_samples * n_wavelengths), n_samples, n_wavelengths)
for (j in 2:n_wavelengths) {
  X[, j] <- 0.7 * X[, j - 1] + 0.3 * X[, j]
}

true_beta <- rep(0, n_wavelengths)
relevant_idx <- c(46, 103, 151, 204, 256, 311, 349, 401, 426, 471)
true_beta[relevant_idx] <- rnorm(length(relevant_idx)) * 2
y <- as.vector(X %*% true_beta + rnorm(n_samples) * 0.5)

# --- Lasso with cross-validation (alpha = 1 selects the pure L1 penalty) ---
# glmnet standardizes the predictors internally by default
cv_fit <- cv.glmnet(X, y, alpha = 1, nfolds = 10)
beta_min <- as.numeric(coef(cv_fit, s = "lambda.min"))[-1]
beta_1se <- as.numeric(coef(cv_fit, s = "lambda.1se"))[-1]
cat(sprintf("lambda.min: %.4f (%d non-zero)\n", cv_fit$lambda.min, sum(beta_min != 0)))
cat(sprintf("lambda.1se: %.4f (%d non-zero)\n", cv_fit$lambda.1se, sum(beta_1se != 0)))

# Cross-validation curve
plot(cv_fit)

# Regularization path
fit <- glmnet(X, y, alpha = 1)
plot(fit, xvar = "lambda")
abline(v = log(c(cv_fit$lambda.min, cv_fit$lambda.1se)), lty = 2, col = c("red", "blue"))

# True versus estimated coefficients
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
plot(true_beta, type = "h", col = "red", lwd = 2,
     main = "True coefficients", xlab = "", ylab = "Coefficient")
plot(beta_min, type = "h", col = "blue", lwd = 2,
     main = sprintf("Lasso estimates (%d non-zero)", sum(beta_min != 0)),
     xlab = "Variable index", ylab = "Coefficient")

cat("Selected:", which(beta_min != 0), "\n")
cat("True relevant:", relevant_idx, "\n")

Do not expect the lasso to return exactly the ten true variables. At ${\lambda }_{\min }$ it selects several times more: 57 in the Python run and 65 in the MATLAB run. The true variables are among them (nine of ten in Python, all ten in MATLAB), together with some of their correlated neighbours and many small coefficients on irrelevant variables. Cross-validation chooses $\lambda$ for prediction, and a few extra small coefficients cost little in prediction error, so it does not deliver a clean selection. The ${\lambda }_{1\text{SE}}$ model is much sparser (26 variables in Python, 24 in MATLAB). In the Python run, variable 400 is never found, and that is correct behaviour: its simulated true coefficient happens to be 0.02, far too small to detect. If you need a reliable list of variables, use stability selection (see the tips below).

Good applications:

• Wavelength selection for NIR, MIR or Raman calibrations, to find the spectral regions that carry the analyte information.
• Compact sensor design: choosing a few wavelengths for a dedicated filter photometer.
• Biomarker discovery: finding which of hundreds of metabolites (from LC-MS or NMR) predict an outcome.
• Any problem where you have good reason to believe that only a few variables matter.

Consider alternatives when:

• the spectra are broad and strongly correlated and the information is spread over whole bands: the lasso's selections become unstable (refit it on replicate calibration sets and the chosen wavelength often jumps from one point of a band to its neighbour, just as the code example above picks neighbours of the true variables), and elastic net, ridge or PLS are usually more robust;
• the dataset is very small, so that cross-validation cannot reliably choose $\lambda$;
• prediction is all you need: in spectroscopy, where sparsity is rarely how the data were generated, PLS or ridge often predict as well or better.

Preprocessing still matters. The lasso selects variables; it does not correct baselines, scatter or noise. Apply your usual spectral preprocessing first.

Look at the path, not just the final model. A wavelength that enters early and stays is more trustworthy than one that appears only at very small $\lambda$.

Check stability. Refit the lasso on many bootstrap samples or subsamples and count how often each variable is selected. Stability selection turns this into a formal procedure [7]: variables chosen in a large fraction of the runs are robust; those selected only occasionally are probably noise.

Beware of the grouping effect. Neighbouring wavelengths are correlated, so the lasso tends to keep one representative of each band and drop the others. The selected wavelength is a representative of a region, not necessarily the important wavelength.

Compare against PLS. PLS is the reference method for spectroscopic calibration. The lasso is most valuable when it predicts about as well as PLS while also offering an interpretable selection.

The lasso is least squares with an absolute-value penalty. Its diamond-shaped constraint region has corners on the axes, and its penalty pushes with a force that does not fade near zero; both are ways of saying that some coefficients become exactly zero. That turns a regression method into a variable selection method, controlled by a single parameter chosen by cross-validation.

Its weakness is the flip side of its strength: with groups of highly correlated variables, which is exactly what spectra are, the choice of which variable to keep can be fragile. Several extensions address this:

• Elastic net mixes L1 and L2 penalties to select correlated variables together [6].
• Group lasso selects or discards predefined groups, such as whole spectral regions [8].
• Adaptive lasso uses data-dependent weights to make the selection more reliable [9].
• Sparse PLS combines the latent variables of PLS with L1 penalties, directly targeting the chemometric setting [10].

For a thorough modern treatment of sparse methods, see Hastie, Tibshirani and Wainwright [11].