Regularisation — Ridge, Lasso, and Elastic Net¶

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

fig_scale = 1.6

def plot_loss_interpretation():
    line = np.linspace(-1.5, 1.5, 1001)
    xx, yy = np.meshgrid(line, line)

    l2          = xx**2 + yy**2
    l1          = np.abs(xx) + np.abs(yy)
    rho         = 0.7
    elastic_net = rho * l1 + (1 - rho) * l2

    plt.figure(figsize=(5 * fig_scale, 4 * fig_scale))
    ax = plt.gca()

    en_c = plt.contour(xx, yy, elastic_net, levels=[1], linewidths=2 * fig_scale, colors='darkorange')
    l2_c = plt.contour(xx, yy, l2,          levels=[1], linewidths=2 * fig_scale, colors='c')
    l1_c = plt.contour(xx, yy, l1,          levels=[1], linewidths=2 * fig_scale, colors='navy')
    ax.set_aspect('equal')
    ax.spines['left'].set_position('center')
    ax.spines['right'].set_color('none')
    ax.spines['bottom'].set_position('center')
    ax.spines['top'].set_color('none')

    plt.clabel(en_c, inline=1, fontsize=12 * fig_scale, fmt={1.0: 'Elastic Net'}, manual=[(-0.6, -0.6)])
    plt.clabel(l2_c, inline=1, fontsize=12 * fig_scale, fmt={1.0: 'L2 (Ridge)'},  manual=[(-0.5, -0.5)])
    plt.clabel(l1_c, inline=1, fontsize=12 * fig_scale, fmt={1.0: 'L1 (Lasso)'},  manual=[(-0.5, -0.5)])

    x1 = np.linspace(0.5, 1.5, 100)
    x2 = np.linspace(-1.0, 1.5, 100)
    X1, X2 = np.meshgrid(x1, x2)
    Y = np.sqrt(np.square(X1 / 2 - 0.7) + np.square(X2 / 4 - 0.28))
    cp = plt.contour(X1, X2, Y)
    plt.clabel(cp, inline=1, fontsize=3)

    ax.tick_params(axis='both', pad=0)
    ax.scatter(1,    0,    c='navy', s=50 * fig_scale, zorder=5, label='Lasso optimum (on axis = zero)')
    ax.scatter(0.89, 0.42, c='c',   s=50 * fig_scale, zorder=5, label='Ridge optimum (off axis)')
    ax.legend(loc='lower right', fontsize=9 * fig_scale)
    ax.set_title('Constraint geometry: L1 corners produce sparsity', fontsize=10 * fig_scale)

    plt.tight_layout()
    plt.show()

plot_loss_interpretation()

%matplotlib inline
import numpy as np
from scipy.linalg import solve
import matplotlib.pyplot as plt

np.random.seed(42)
m = 100
X_raw = 6 * np.random.rand(m, 1) - 3
y_raw = 0.5 * X_raw**2 + X_raw + 2 + np.random.randn(m, 1)
X_new = np.linspace(-3, 3, 100).reshape(-1, 1)

degree = 30
def create_polynomial_features(X, degree):
    X_poly = np.ones((X.shape[0], 1))
    for d in range(1, degree + 1):
        X_poly = np.hstack((X_poly, X**d))
    return X_poly

X_poly     = create_polynomial_features(X_raw, degree)
X_new_poly = create_polynomial_features(X_new, degree)
mean_ = np.mean(X_poly[:, 1:], axis=0)
std_  = np.std(X_poly[:, 1:], axis=0); std_[std_ < 1e-8] = 1.0
X_ps     = X_poly.copy();     X_ps[:, 1:]     = (X_poly[:, 1:]     - mean_) / std_
X_new_ps = X_new_poly.copy(); X_new_ps[:, 1:] = (X_new_poly[:, 1:] - mean_) / std_

def linear_normal_equation(X, y):
    try:
        return solve(X.T @ X, X.T @ y, assume_a='pos')
    except Exception:
        return np.linalg.pinv(X.T @ X) @ (X.T @ y)

def linear_gradient_descent(X, y, learning_rate=0.01, n_iterations=1000):
    n_s, n_f = X.shape
    theta = np.zeros((n_f, 1))
    for _ in range(n_iterations):
        theta -= learning_rate * (1/n_s) * X.T @ (X @ theta - y)
    return theta

def ridge_normal_equation(X, y, lambda_reg):
    n_f = X.shape[1]
    I = np.eye(n_f); I[0, 0] = 0
    return solve(X.T @ X + lambda_reg * I, X.T @ y, assume_a='pos')

def ridge_gradient_descent(X, y, lambda_reg, learning_rate=0.01, n_iterations=1000):
    n_s, n_f = X.shape
    theta = np.zeros((n_f, 1))
    for _ in range(n_iterations):
        penalty = lambda_reg * np.vstack([0, theta[1:]])
        theta  -= learning_rate * ((1/n_s) * X.T @ (X @ theta - y) + penalty)
    return theta

def lasso_coordinate_descent(X, y, lambda_reg, max_iter=1000, tol=1e-4):
    n_s, n_f = X.shape
    theta = np.zeros((n_f, 1))
    y = y.reshape(-1, 1)
    for _ in range(max_iter):
        old = theta.copy()
        for j in range(n_f):
            Xj  = X[:, j:j+1]
            rho = (Xj * (y - X @ theta + Xj * theta[j])).sum() / n_s
            if j == 0:
                theta[j] = rho
            else:
                denom = float(Xj.T @ Xj) / n_s
                if   rho < -lambda_reg: theta[j] = (rho + lambda_reg) / denom
                elif rho >  lambda_reg: theta[j] = (rho - lambda_reg) / denom
                else:                   theta[j] = 0.0
        if np.linalg.norm(theta - old) < tol:
            break
    return theta

lam = 0.1
lr  = 0.005
iters = 2000
th_ols_ne  = linear_normal_equation(X_ps, y_raw)
th_ols_gd  = linear_gradient_descent(X_ps, y_raw, lr, iters)
th_ridge_ne = ridge_normal_equation(X_ps, y_raw, lam)
th_ridge_gd = ridge_gradient_descent(X_ps, y_raw, lam, lr, iters)
th_lasso    = lasso_coordinate_descent(X_ps, y_raw, lam)

fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(X_raw, y_raw, color='black', alpha=0.4, s=15, label='Data', zorder=3)
for th, color, ls, label in [
    (th_ols_ne,   'orange',    '-',  'OLS (Normal Eq)'),
    (th_ols_gd,   'black',     '--', 'OLS (Grad Desc)'),
    (th_ridge_ne, 'steelblue', '-',  f'Ridge NE (λ={lam})'),
    (th_ridge_gd, 'green',     '--', f'Ridge GD (λ={lam})'),
    (th_lasso,    'tomato',    '-',  f'Lasso CD (λ={lam})'),
]:
    ax.plot(X_new, np.clip(X_new_ps @ th, -5, 20), color=color, linestyle=ls, linewidth=1.8, label=label)
ax.set_ylim(-3, 18)
ax.legend(fontsize=9)
ax.set_title(f'Degree-{degree} polynomial: OLS vs Ridge vs Lasso')
ax.set_xlabel('x'); ax.set_ylabel('y')
plt.tight_layout(); plt.show()

nz = int(np.sum(np.abs(th_lasso) > 1e-6))
print(f'Lasso: {nz} non-zero coefficients out of {th_lasso.shape[0]}')

%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
from sklearn.linear_model import Ridge, lars_path

X_d, y_d = load_diabetes(return_X_y=True)

alphas = np.logspace(-5, 2, 60)
ridge_coefs = [Ridge(alpha=a, fit_intercept=False).fit(X_d, y_d).coef_ for a in alphas]

_, _, lasso_coefs = lars_path(X_d, y_d, method='lasso')
xx = np.sum(np.abs(lasso_coefs.T), axis=1)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(alphas, ridge_coefs)
axes[0].set_xscale('log')
axes[0].set_xlabel(r'Regularisation strength $\lambda$')
axes[0].set_ylabel('Coefficient magnitude')
axes[0].set_title(r'Ridge: all coefficients shrink smoothly toward 0')
axes[0].axis('tight')

axes[1].plot(3500 - xx, lasso_coefs.T)
axes[1].set_xlabel(r'Regularisation strength $\lambda$')
axes[1].set_ylabel('Coefficient magnitude')
axes[1].set_title('Lasso: coefficients hit zero and stay there (sparsity)')
axes[1].axis('tight')

plt.suptitle('Coefficient paths — Ridge vs Lasso on the diabetes dataset', y=1.01)
plt.tight_layout()
plt.show()

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import make_pipeline

np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = (0.5 * X**2 + X + 2 + np.random.randn(m, 1)).ravel()
X_new = np.linspace(-3, 3, 200).reshape(-1, 1)

degree = 15
configs = [
    ('OLS (no reg)',    make_pipeline(PolynomialFeatures(degree), StandardScaler(), Ridge(alpha=0)),      'grey',      '--'),
    ('Ridge α=1',      make_pipeline(PolynomialFeatures(degree), StandardScaler(), Ridge(alpha=1)),      'steelblue', '-'),
    ('Ridge α=10',     make_pipeline(PolynomialFeatures(degree), StandardScaler(), Ridge(alpha=10)),     'blue',      '-'),
    ('Lasso α=0.1',    make_pipeline(PolynomialFeatures(degree), StandardScaler(), Lasso(alpha=0.1, max_iter=10000)),   'tomato', '-'),
    ('ElasticNet',     make_pipeline(PolynomialFeatures(degree), StandardScaler(), ElasticNet(alpha=0.5, l1_ratio=0.5, max_iter=10000)), 'purple', '-'),
]

fig, ax = plt.subplots(figsize=(11, 5))
ax.scatter(X, y, color='black', alpha=0.3, s=15, label='Data', zorder=3)
for label, pipe, color, ls in configs:
    pipe.fit(X, y)
    ax.plot(X_new, np.clip(pipe.predict(X_new), -5, 20), color=color, linewidth=2, linestyle=ls, label=label)

ax.set_ylim(-3, 18)
ax.legend(fontsize=9)
ax.set_title(f'Degree-{degree} polynomial with different regularisers')
ax.set_xlabel('x'); ax.set_ylabel('y')
plt.tight_layout(); plt.show()

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso
from sklearn.datasets import make_regression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler

X_cv, y_cv = make_regression(n_samples=50, n_features=100,
                             n_informative=5, noise=20, random_state=42)
scaler = StandardScaler()
X_cv_s = scaler.fit_transform(X_cv)

alphas = np.logspace(-2, 4, 100)
ridge_mse, lasso_mse = [], []

for alpha in alphas:
    ridge_mse.append(np.mean(-cross_val_score(
        Ridge(alpha=alpha, max_iter=10000), X_cv_s, y_cv,
        cv=5, scoring='neg_mean_squared_error')))
    lasso_mse.append(np.mean(-cross_val_score(
        Lasso(alpha=alpha, max_iter=10000), X_cv_s, y_cv,
        cv=5, scoring='neg_mean_squared_error')))

best_ridge = alphas[np.argmin(ridge_mse)]
best_lasso = alphas[np.argmin(lasso_mse)]

fig, ax = plt.subplots(figsize=(10, 5))
ax.semilogx(alphas, ridge_mse, color='steelblue', linewidth=2, label='Ridge')
ax.semilogx(alphas, lasso_mse, color='tomato',    linewidth=2, label='Lasso')
ax.scatter(best_ridge, min(ridge_mse), color='steelblue', s=100, zorder=5,
           label=f'Best Ridge (alpha={best_ridge:.1f})')
ax.scatter(best_lasso, min(lasso_mse), color='tomato',    s=100, zorder=5,
           label=f'Best Lasso (alpha={best_lasso:.1f})')
ax.axhline(min(ridge_mse), color='steelblue', linestyle=':', alpha=0.4)
ax.axhline(min(lasso_mse), color='tomato',    linestyle=':', alpha=0.4)
ax.set_xlabel(r'Regularisation strength $\alpha$')
ax.set_ylabel('5-fold CV MSE')
ax.set_title('Cross-validation alpha search\n(50 samples, 100 features, 5 informative)')
ax.legend(); ax.grid(True, which='both', alpha=0.3)
plt.tight_layout(); plt.show()

print(f"OLS (alpha~0) CV MSE : {ridge_mse[0]:.1f}")
print(f"Best Ridge CV MSE    : {min(ridge_mse):.1f}  (improvement: {100*(ridge_mse[0]-min(ridge_mse))/ridge_mse[0]:.1f}%)")
print(f"Best Lasso CV MSE    : {min(lasso_mse):.1f}  (improvement: {100*(lasso_mse[0]-min(lasso_mse))/lasso_mse[0]:.1f}%)")

%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

rng = np.random.default_rng(0)
n = 80
X_e = rng.uniform(-3, 3, (n, 1))
y_e = (0.5*X_e.ravel()**2 + X_e.ravel() + 2 + rng.normal(0, 1, n)).reshape(-1, 1)

pf = PolynomialFeatures(10, include_bias=True)
Xp = pf.fit_transform(X_e)
sc = StandardScaler(with_mean=False); Xp = sc.fit_transform(Xp)

# TODO: implement Ridge gradient descent
# gradient = (2/n)*X^T(X@theta - y) + 2*lambda*theta  (skip bias at index 0)
# Run 2000 steps with lr=0.01, lambda=0.5
# Plot fitted curve vs data

import numpy as np
from sklearn.linear_model import LassoCV, Ridge
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(7)
n, p = 100, 50
X2 = rng.normal(0, 1, (n, p))
true_coef = np.zeros(p)
true_coef[[3, 17, 42]] = [2.5, -1.8, 3.1]   # only 3 features matter
y2 = X2 @ true_coef + rng.normal(0, 1, n)

X2s = StandardScaler().fit_transform(X2)

# TODO: fit LassoCV, print which features have non-zero coefficients
# TODO: compare with Ridge — does Ridge zero out the noise features?

from sklearn.linear_model import LassoCV, RidgeCV, ElasticNetCV
import numpy as np

# Use X2s and y2 from Exercise 2
# TODO: fit LassoCV, RidgeCV, ElasticNetCV
# Print best alpha and 5-fold CV MSE for each