Sparsity-Promoting Regularization

Hide code cell source

 
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns

Sparsity-Promoting Regularization#

\(L^1\)-regularization is a method to prevent overfitting in linear models. It is also known as Lasso regression. It adds a penalty term to the loss function, which is the sum of the absolute values of the coefficients. This forces the coefficients to be small, and some of them to be exactly zero. This is useful for feature selection, as it allows us to discard some of the features.

Suppose we have \(N\) input–response pairs \((x_i,y_i)\) and \(M\) candidate features collected in \(\boldsymbol{\phi}(x)\). We fit the linear basis-function model:

\[ y = w_0 + \mathbf{w}^T\boldsymbol{\phi}(x) + \epsilon. \]

Here \(\mathbf{w}\) contains the feature coefficients, \(w_0\) is the bias or intercept, and \(\epsilon\) is observation noise.

We can add the \(L^1\)-regularization term to the usual sum of square errors loss function:

\[ L = \frac{1}{N}\sum_{i=1}^N(y_i - w_0 - \mathbf{w}^T\boldsymbol{\phi}(x_i))^2 + \lambda\|\mathbf{w}\|_1, \]

where \(\lambda\geq0\) controls the penalty strength. The \(L^1\)-norm of the vector \(\mathbf{w}\) is defined as:

\[ \|\mathbf{w}\|_1 = \sum_{j=1}^M|w_j|. \]

Notice that the bias term \(w_0\) is not regularized. We don’t want to push it to zero, as it is a constant term, and it is not related to the features.

Let’s illustrate \(L^1\)-regularization with a simple example. We generate noisy data from a cubic polynomial and fit a larger library of polynomial and Fourier features. The penalty helps us select a smaller subset of terms.

Here are the data:

import numpy as np

rng = np.random.default_rng(697)
num_samples = 100
f = lambda x: 0.2 * x + 0.3 * x**2 + 0.4 * x**3
sigma = 0.01
x = rng.uniform(-1, 1, num_samples)
y = f(x) + rng.normal(0, sigma, num_samples)

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(x, y, '.', color='k')
ax.set(xlabel='x', ylabel='y')
finalize_axes(keep_box=False)

array([<Axes: xlabel='x', ylabel='y'>], dtype=object)
Noisy samples from a cubic polynomial.

We will try to fit a ridiculously complicated model. We will have a degree-ten polynomial, and then we will also add a bunch of Fourier terms:

from jax import config
config.update("jax_enable_x64", True)
import jax.numpy as jnp
from jax import vmap

@vmap
def features(x):
    tmp1 = jnp.array([x ** i for i in range(1, 11)])
    tmp2 = jnp.array([jnp.sin(2 * np.pi * x * i) for i in range(1, 11)])
    return jnp.concatenate([tmp1, tmp2])

Visualize the features:

xs = jnp.linspace(-1, 1, 200)
Phi = features(xs)
fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(xs, Phi, lw=0.2)
ax.set(xlabel='x', ylabel='feature value')
finalize_axes(keep_box=False)

array([<Axes: xlabel='x', ylabel='feature value'>], dtype=object)
Twenty polynomial feature functions over the input interval.

Clearly an overkill. But Lasso regression will still be able to throw away the unnecessary terms.

In what follows we do Lasso regression using the sklearn library. We use LassoCV, which also does cross-validation to find the best value of the regularization parameter \(\lambda\).

from sklearn import linear_model

Phi_train = features(x)
clf = linear_model.LassoCV(fit_intercept=True, max_iter=100000, tol=1e-8) 
clf.fit(Phi_train, y)
lr = linear_model.RidgeCV(fit_intercept=True)
lr.fit(Phi_train, y)

The library contains ten nonconstant monomials and ten sine terms. Both fits estimate a separate, unpenalized intercept, which is excluded from the coefficient panels. Scikit-learn uses a factor \(1/(2N)\) for the Lasso squared-error term, so its parameter \(\alpha\) corresponds to \(\lambda/2\) in our convention.

For comparison, we do Ridge regression (which adds an L2-regularization term) also with cross-validation.

fig, axes = plt.subplots(2, 1, figsize=(7.35, 2.2))
axes[0].imshow(clf.coef_.reshape(1, 20), cmap='gray', vmin=-0.5, vmax=0.5)
axes[0].set(xlabel='feature index', title='Lasso',
            xticks=np.arange(0, 20, 2), xticklabels=np.arange(0, 20, 2),
            yticks=[])
axes[1].imshow(lr.coef_.reshape(1, 20), cmap='gray', vmin=-0.5, vmax=0.5)
axes[1].set(xlabel='feature index', title='Ridge',
            xticks=np.arange(0, 20, 2), xticklabels=np.arange(0, 20, 2),
            yticks=[])
finalize_axes(keep_box=True)

array([<Axes: title={'center': 'Lasso'}, xlabel='feature index'>,
       <Axes: title={'center': 'Ridge'}, xlabel='feature index'>],
      dtype=object)
Lasso and ridge coefficient values across the twenty polynomial and Fourier features.

We see that Ridge regression has a lot more non-zero coefficients than Lasso regression.

Let’s visualize the predictions as well:

Phi_test = features(xs)
lass_pred = clf.predict(Phi_test)
ridge_pred = lr.predict(Phi_test)

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(x, y, '.', color='k')
ax.plot(xs, f(xs), lw=2, label='true')
ax.plot(xs, lass_pred, '--', lw=2, label='Lasso')
ax.plot(xs, ridge_pred, '-.', lw=2, label='Ridge')
ax.set(xlabel='x', ylabel='y')
ax.legend(loc='best')
finalize_axes(keep_box=False)

array([<Axes: xlabel='x', ylabel='y'>], dtype=object)
Noisy data and the true, Lasso, and ridge fitted curves.

Geometry of the \(L^1\) penalty#

To understand how Lasso regression works, we need to look at the \(L^1\)-norm. Let’s work with just two parameters: \(\mathbf{w} = (w_1,w_2)\). The \(L^1\)-norm is:

\[ \|\mathbf{w}\|_1 = |w_1| + |w_2|. \]

This is also called the Manhattan norm. The name reflects the distance traveled along the grid of city blocks.

Let’s consider a simple quadratic objective function:

\[ f(\mathbf{w}) = (\mathbf{w}-\mathbf{w}^*)^T\mathbf{A}(\mathbf{w}-\mathbf{w}^*), \]

where \(\mathbf{A}\) is symmetric positive definite, so the minimum is at \(\mathbf{w}^*\). We seek the minimum of the \(L^1\)-regularized objective function

\[ g(\mathbf{w}) = f(\mathbf{w}) + \lambda\|\mathbf{w}\|_1. \]

To understand this minimum, first we need to mention that Lasso regression is equivalent to the following constrained optimization problem:

\[ \min_{\mathbf{w}}f(\mathbf{w})\quad\text{subject to}\quad \|\mathbf{w}\|_1\leq t. \]

To see that this is equivalent to the \(L^1\)-regularized problem, we can use the Lagrange multiplier method. For a penalized minimizer \(\widehat{\mathbf{w}}_\lambda\), setting \(t=\|\widehat{\mathbf{w}}_\lambda\|_1\) gives the same constrained minimizer. Conversely, a Lagrange multiplier supplies a corresponding penalty for a positive constraint radius. This correspondence need not be one-to-one: several penalties can give the zero solution, and sufficiently large radii leave the unpenalized solution unchanged.

Now, we can understand what is going on by looking at the level sets of \(f(\mathbf{w})\) and \(\|\mathbf{w}\|_1\).

Hide code cell source

 
w_star = jnp.array([0.5, 0.2])
A = jnp.array([[1, -0.4], [-0.4, 1]])
t = 0.25
w_constrained = jnp.array([t, 0.0])

from jax import jit

f = jit(vmap(lambda w: (w - w_star) @ A @ (w - w_star)))
minimum = float(f(w_constrained[None, :])[0])
# The corner satisfies the convex optimality condition:
# grad f = (-0.34, -0.20), with |grad f_2| <= -grad f_1.
grad_at_corner = 2 * A @ (w_constrained - w_star)
assert np.allclose(grad_at_corner, [-0.34, -0.20])
assert abs(float(grad_at_corner[1])) <= -float(grad_at_corner[0])

w1 = jnp.linspace(-0.45, 0.85, 500)
w2 = jnp.linspace(-0.45, 0.65, 500)
W1, W2 = jnp.meshgrid(w1, w2)
W_flat = jnp.stack([W1.ravel(), W2.ravel()], axis=1)
F = f(W_flat).reshape(W1.shape)

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_tall"], constrained_layout=True)
ax.contour(W1, W2, F, levels=[0.01, 0.03, minimum, 0.10],
           colors="0.45", linewidths=0.8)
ax.contour(W1, W2, F, levels=[minimum], colors="black", linewidths=1.25)
ax.plot([t, 0, -t, 0, t], [0, t, 0, -t, 0], color="black", linewidth=1.4)
ax.plot(*w_star, marker="x", color="black", label="Unconstrained minimum")
ax.plot(*w_constrained, marker="o", color="black", markerfacecolor="white",
        markeredgecolor="black", label="Constrained minimum")
ax.set(xlabel="$w_1$", ylabel="$w_2$", aspect="equal")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.18), fontsize=7)
finalize_axes(keep_box=True)

array([<Axes: xlabel='$w_1$', ylabel='$w_2$'>], dtype=object)
Quadratic-loss contours meeting the diamond-shaped L1 constraint at a sparse solution.

For the plotted radius \(t=0.25\), the constrained minimum is \((0.25,0)\), where the quadratic contour touches a corner of the diamond. You are trying to get to the minimum of \(f(\mathbf{w})\) while remaining inside the diamond \(\|\mathbf{w}\|_1 \leq t\). The level sets of \(\|\mathbf{w}\|_1\) are diamonds. The constrained minimum of \(f(\mathbf{w})\) will be where a level set of \(f(\mathbf{w})\) first touches the diamond. You will often find yourself at a corner of the diamond. Then one of the components of \(\mathbf{w}\) will be zero.

Bayesian Interpretation of Lasso Regression#

Lasso regression can be interpreted as a Bayesian method. To see this, put a Laplace prior on each one of the weights:

\[ w_j \sim \text{Laplace}(0, b). \]

The density of the Laplace distribution is:

\[ p(w_j) = \frac{1}{2b}\exp\left(-\frac{|w_j|}{b}\right). \]

Taking the logarithm of the Laplace prior, we get:

\[ \log p(w_j) = -\frac{|w_j|}{b} + \text{const}. \]

Therefore, we see that (assuming a Gaussian likelihood), the \(L^1\)-regularized objective function is (up to an additive constant) equal to the negative log posterior.

Let’s visualize the Laplace and Gaussian densities:

import scipy.stats as st

xs = np.linspace(-1, 1, 300)
fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(xs, st.norm.pdf(xs, 0, 0.1), color='black', linestyle='-',
        label='Normal')
ax.plot(xs, st.laplace.pdf(xs, 0, 0.1), color='black', linestyle='--',
        label='Laplace')
ax.legend()
finalize_axes(keep_box=False)

array([<Axes: >], dtype=object)
Normal and Laplace densities with solid and dashed black curves.