Master Regularization Techniques

Lasso regression takes regularization a step further by using the absolute value of coefficients instead of their squares. This seemingly small change has a profound effect: Lasso can drive coefficients to exactly zero, effectively performing automatic feature selection. When you need a sparse model that uses only the most important features, Lasso is your tool of choice.

The Lasso Cost Function

The Lasso cost function replaces the squared penalty (L2) with an absolute value penalty (L1). This seemingly minor change fundamentally alters the optimization landscape and enables the sparse solutions that make Lasso unique.

Why L1 Produces Zeros (Geometric Intuition)

The reason Lasso produces sparse solutions while Ridge does not comes down to geometry. Imagine the constraint regions as shapes on a 2D plane where we're optimizing two coefficients. Ridge creates a circular constraint (all points at equal distance from origin), while Lasso creates a diamond shape with corners on the axes.

The optimization process finds where the elliptical contours of the MSE loss function first touch the constraint region. For Ridge's circle, this contact point is almost never on an axis (zero coefficient). For Lasso's diamond, the contact point often hits a corner, which lies on an axis - meaning one coefficient is exactly zero!

Lasso Regression in Python

import numpy as np
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

# Generate data with many irrelevant features
np.random.seed(42)
n_samples, n_features = 100, 20
X = np.random.randn(n_samples, n_features)

# Only first 3 features actually matter
y = 5*X[:, 0] - 3*X[:, 1] + 2*X[:, 2] + np.random.randn(n_samples) * 0.5

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Train Lasso
lasso = Lasso(alpha=0.1)
lasso.fit(X_scaled, y)

# Count non-zero coefficients
n_nonzero = np.sum(lasso.coef_ != 0)
print(f"Non-zero coefficients: {n_nonzero} out of {n_features}")
print(f"Selected features: {np.where(lasso.coef_ != 0)[0]}")
print(f"Coefficients: {lasso.coef_[lasso.coef_ != 0].round(3)}")

Cross-Validation for Alpha Selection

Choosing the right alpha value is critical for Lasso performance. Scikit-learn provides LassoCV which automatically finds the optimal alpha using cross-validation.

from sklearn.linear_model import LassoCV
import matplotlib.pyplot as plt

# LassoCV automatically finds optimal alpha
lasso_cv = LassoCV(alphas=np.logspace(-4, 1, 100), cv=5)
lasso_cv.fit(X_scaled, y)

print(f"Optimal alpha: {lasso_cv.alpha_:.6f}")
print(f"R² score: {lasso_cv.score(X_scaled, y):.4f}")
print(f"Non-zero features: {np.sum(lasso_cv.coef_ != 0)}")

# Plot alpha vs MSE
plt.figure(figsize=(10, 5))
plt.semilogx(lasso_cv.alphas_, lasso_cv.mse_path_.mean(axis=1))
plt.axvline(lasso_cv.alpha_, color='r', linestyle='--', label=f'Optimal alpha={lasso_cv.alpha_:.4f}')
plt.xlabel('Alpha')
plt.ylabel('Mean Squared Error')
plt.title('Lasso Cross-Validation')
plt.legend()
plt.show()

import numpy as np

coefficients = np.array([2.5, -1.8, 0.0, 0.9, 0.0])
alpha = 0.3

# L1 penalty = alpha * sum(|coefficients|)
l1_penalty = alpha * np.sum(np.abs(coefficients))
# = 0.3 * (2.5 + 1.8 + 0 + 0.9 + 0)
# = 0.3 * 5.2
# = 1.56

print(f"L1 Penalty: {l1_penalty}")
# Output: L1 Penalty: 1.56

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

np.random.seed(42)
n_samples, n_features = 200, 50
X = np.random.randn(n_samples, n_features)

# True important features: indices 0, 10, 20, 30, 40
important_features = [0, 10, 20, 30, 40]
true_coefs = [3.0, -2.5, 4.0, -1.5, 2.0]

y = sum(c * X[:, i] for i, c in zip(important_features, true_coefs))
y += np.random.randn(n_samples) * 0.5  # Add noise

# Scale and fit Lasso
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

lasso = LassoCV(cv=5)
lasso.fit(X_scaled, y)

# Check which features were selected
selected = np.where(lasso.coef_ != 0)[0]
print(f"True important features: {important_features}")
print(f"Lasso selected features: {list(selected)}")
print(f"Match: {set(selected) == set(important_features)}")

import numpy as np

def soft_threshold(rho, alpha):
    """Soft thresholding operator for L1 regularization"""
    if rho < -alpha:
        return rho + alpha
    elif rho > alpha:
        return rho - alpha
    else:
        return 0.0

def lasso_coordinate_descent(X, y, alpha, n_iterations=1000):
    """Lasso using coordinate descent"""
    n_samples, n_features = X.shape
    weights = np.zeros(n_features)
    
    for iteration in range(n_iterations):
        for j in range(n_features):
            # Calculate residual excluding feature j
            residual = y - X @ weights + weights[j] * X[:, j]
            
            # Calculate rho (correlation)
            rho = X[:, j] @ residual / n_samples
            
            # Apply soft thresholding
            z = X[:, j] @ X[:, j] / n_samples
            weights[j] = soft_threshold(rho, alpha) / z
    
    return weights

# Test on data
np.random.seed(42)
X = np.random.randn(100, 10)
y = 3*X[:, 0] - 2*X[:, 1] + np.random.randn(100) * 0.1

weights = lasso_coordinate_descent(X, y, alpha=0.1)
print("Non-zero weights:", np.where(weights != 0)[0])
print("Weight values:", weights[weights != 0].round(4))

What if you want both the feature selection of Lasso and the stability of Ridge? Elastic Net combines L1 and L2 regularization into a single model, giving you fine-grained control over the balance between sparsity and coefficient shrinkage. It's particularly powerful when dealing with correlated features where Lasso alone might be unstable.

The Elastic Net Cost Function

In scikit-learn, this is parameterized slightly differently using alpha and l1_ratio:

Elastic Net in Python

from sklearn.linear_model import ElasticNet, ElasticNetCV
from sklearn.preprocessing import StandardScaler
import numpy as np

# Generate data with correlated features
np.random.seed(42)
n_samples = 200
base_feature = np.random.randn(n_samples)

# Create correlated features
X = np.column_stack([
    base_feature,                              # Feature 0
    base_feature + np.random.randn(n_samples) * 0.1,  # Correlated with 0
    base_feature + np.random.randn(n_samples) * 0.1,  # Correlated with 0
    np.random.randn(n_samples),                # Independent
    np.random.randn(n_samples),                # Independent
])

y = 3 * base_feature + np.random.randn(n_samples) * 0.5

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Compare Lasso vs Elastic Net on correlated data
from sklearn.linear_model import Lasso

lasso = Lasso(alpha=0.1)
lasso.fit(X_scaled, y)
print(f"Lasso coefficients: {lasso.coef_.round(3)}")

elastic = ElasticNet(alpha=0.1, l1_ratio=0.5)
elastic.fit(X_scaled, y)
print(f"Elastic Net coefficients: {elastic.coef_.round(3)}")

Cross-Validation for Hyperparameter Tuning

from sklearn.linear_model import ElasticNetCV
import matplotlib.pyplot as plt

# ElasticNetCV searches over both alpha and l1_ratio
elastic_cv = ElasticNetCV(
    l1_ratio=[0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 0.99],
    alphas=np.logspace(-4, 1, 50),
    cv=5
)
elastic_cv.fit(X_scaled, y)

print(f"Optimal alpha: {elastic_cv.alpha_:.6f}")
print(f"Optimal l1_ratio: {elastic_cv.l1_ratio_:.2f}")
print(f"R² score: {elastic_cv.score(X_scaled, y):.4f}")
print(f"Non-zero coefficients: {np.sum(elastic_cv.coef_ != 0)}")
print(f"Coefficients: {elastic_cv.coef_.round(4)}")

When to Use Each Method

import numpy as np

coefficients = np.array([2.0, -1.0, 0.5])
alpha = 0.4
l1_ratio = 0.7

# L1 component: alpha * l1_ratio * sum(|coef|)
l1_penalty = alpha * l1_ratio * np.sum(np.abs(coefficients))
# = 0.4 * 0.7 * (2.0 + 1.0 + 0.5) = 0.28 * 3.5 = 0.98

# L2 component: alpha * (1-l1_ratio) * sum(coef^2)
l2_penalty = alpha * (1 - l1_ratio) * np.sum(coefficients ** 2)
# = 0.4 * 0.3 * (4.0 + 1.0 + 0.25) = 0.12 * 5.25 = 0.63

total_penalty = l1_penalty + l2_penalty
print(f"L1 penalty: {l1_penalty:.2f}")
print(f"L2 penalty: {l2_penalty:.2f}")
print(f"Total penalty: {total_penalty:.2f}")
# Output: Total penalty: 1.61

import numpy as np
from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler

# Generate dataset
np.random.seed(42)
X = np.random.randn(200, 20)
y = 3*X[:, 0] - 2*X[:, 1] + 1.5*X[:, 2] + np.random.randn(200) * 0.5

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Compare models
models = {
    'Ridge': Ridge(alpha=1.0),
    'Lasso': Lasso(alpha=0.1),
    'Elastic Net': ElasticNet(alpha=0.1, l1_ratio=0.5)
}

print(f"{'Model':<15} {'Non-zero':<10} {'R² (CV)':<10}")
print("-" * 35)

for name, model in models.items():
    model.fit(X_scaled, y)
    n_nonzero = np.sum(model.coef_ != 0)
    cv_score = cross_val_score(model, X_scaled, y, cv=5).mean()
    print(f"{name:<15} {n_nonzero:<10} {cv_score:.4f}")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import GridSearchCV
from sklearn.preprocessing import StandardScaler

# Generate data
np.random.seed(42)
X = np.random.randn(200, 15)
y = 2*X[:, 0] - 3*X[:, 1] + X[:, 2] + np.random.randn(200) * 0.3

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Grid search
param_grid = {
    'alpha': [0.001, 0.01, 0.1, 1.0, 10.0],
    'l1_ratio': [0.1, 0.3, 0.5, 0.7, 0.9]
}

grid_search = GridSearchCV(
    ElasticNet(), param_grid, cv=5, scoring='r2', return_train_score=True
)
grid_search.fit(X_scaled, y)

print(f"Best params: {grid_search.best_params_}")
print(f"Best R²: {grid_search.best_score_:.4f}")

# Create heatmap
scores = grid_search.cv_results_['mean_test_score'].reshape(5, 5)
plt.figure(figsize=(8, 6))
plt.imshow(scores, cmap='viridis', aspect='auto')
plt.colorbar(label='R² Score')
plt.xticks(range(5), param_grid['l1_ratio'])
plt.yticks(range(5), param_grid['alpha'])
plt.xlabel('l1_ratio')
plt.ylabel('alpha')
plt.title('Elastic Net Grid Search Results')
plt.show()

Overfitting is the nemesis of machine learning models - it occurs when a model learns the training data too well, including its noise and peculiarities, failing to generalize to new data. Regularization is our primary weapon against overfitting, constraining model complexity to find the sweet spot between underfitting and overfitting.

The Bias-Variance Tradeoff

Understanding the bias-variance tradeoff is essential to grasping why regularization works. Every model's prediction error can be decomposed into three components: bias (systematic error), variance (sensitivity to training data), and irreducible noise.

Visualizing Overfitting

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.pipeline import make_pipeline

# Generate noisy data
np.random.seed(42)
X = np.linspace(0, 1, 30).reshape(-1, 1)
y = np.sin(2 * np.pi * X).ravel() + np.random.randn(30) * 0.3

X_test = np.linspace(0, 1, 100).reshape(-1, 1)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Underfitting: Degree 1
model1 = make_pipeline(PolynomialFeatures(1), LinearRegression())
model1.fit(X, y)
axes[0].scatter(X, y, alpha=0.7)
axes[0].plot(X_test, model1.predict(X_test), 'r-', linewidth=2)
axes[0].set_title('Underfitting (Degree 1)')

# Just right: Degree 4 with regularization
model2 = make_pipeline(PolynomialFeatures(4), Ridge(alpha=0.1))
model2.fit(X, y)
axes[1].scatter(X, y, alpha=0.7)
axes[1].plot(X_test, model2.predict(X_test), 'g-', linewidth=2)
axes[1].set_title('Good Fit (Degree 4 + Ridge)')

# Overfitting: Degree 15 no regularization
model3 = make_pipeline(PolynomialFeatures(15), LinearRegression())
model3.fit(X, y)
axes[2].scatter(X, y, alpha=0.7)
axes[2].plot(X_test, model3.predict(X_test), 'orange', linewidth=2)
axes[2].set_title('Overfitting (Degree 15)')

plt.tight_layout()
plt.show()

Learning Curves: Diagnosing Overfitting

Learning curves plot model performance against training set size, revealing whether your model suffers from bias or variance. They're invaluable for diagnosing overfitting and guiding model selection.

from sklearn.model_selection import learning_curve
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

def plot_learning_curve(model, X, y, title):
    train_sizes, train_scores, val_scores = learning_curve(
        model, X, y, cv=5, 
        train_sizes=np.linspace(0.1, 1.0, 10),
        scoring='neg_mean_squared_error'
    )

    train_mean = -train_scores.mean(axis=1)
    val_mean = -val_scores.mean(axis=1)

    plt.figure(figsize=(8, 5))
    plt.plot(train_sizes, train_mean, 'o-', label='Training Error')
    plt.plot(train_sizes, val_mean, 'o-', label='Validation Error')
    plt.xlabel('Training Set Size')
    plt.ylabel('Mean Squared Error')
    plt.title(title)
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

# Generate data
np.random.seed(42)
X = np.random.randn(200, 1)
y = 0.5 * X.ravel() ** 2 + np.random.randn(200) * 0.3

# High variance (overfitting) - weak regularization
plot_learning_curve(
    make_pipeline(PolynomialFeatures(10), Ridge(alpha=0.001)),
    X, y, 'Weak Regularization (Overfitting)'
)

# Good fit - proper regularization
plot_learning_curve(
    make_pipeline(PolynomialFeatures(10), Ridge(alpha=1.0)),
    X, y, 'Proper Regularization (Good Fit)'
)

Regularization Path Visualization

A regularization path shows how coefficients change as regularization strength (alpha) varies. This visualization helps you understand which features are most important and how aggressive the regularization needs to be.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso

# Generate data
np.random.seed(42)
X = np.random.randn(100, 5)
y = 3*X[:, 0] + 2*X[:, 1] - X[:, 2] + np.random.randn(100) * 0.5

# Test many alpha values
alphas = np.logspace(-3, 3, 100)

# Store coefficients
ridge_coefs = []
lasso_coefs = []

for alpha in alphas:
    ridge = Ridge(alpha=alpha)
    ridge.fit(X, y)
    ridge_coefs.append(ridge.coef_)
    
    lasso = Lasso(alpha=alpha, max_iter=10000)
    lasso.fit(X, y)
    lasso_coefs.append(lasso.coef_)

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

for i in range(5):
    axes[0].plot(np.log10(alphas), [c[i] for c in ridge_coefs], label=f'Feature {i}')
    axes[1].plot(np.log10(alphas), [c[i] for c in lasso_coefs], label=f'Feature {i}')

axes[0].set_xlabel('log10(Alpha)')
axes[0].set_ylabel('Coefficient')
axes[0].set_title('Ridge Regularization Path')
axes[0].legend()

axes[1].set_xlabel('log10(Alpha)')
axes[1].set_ylabel('Coefficient')
axes[1].set_title('Lasso Regularization Path')
axes[1].legend()

plt.tight_layout()
plt.show()

# Analysis
train_r2 = 0.99
test_r2 = 0.45
gap = train_r2 - test_r2  # 0.54

print(f"Training R²: {train_r2}")
print(f"Test R²: {test_r2}")
print(f"Gap: {gap}")

# Diagnosis
diagnosis = """
This is a clear case of OVERFITTING:
- Training score (0.99) is almost perfect
- Test score (0.45) is much lower
- Large gap (0.54) between train and test

Solutions:
1. Increase regularization (higher alpha)
2. Reduce model complexity (fewer features)
3. Get more training data
4. Apply cross-validation for tuning
5. Try Lasso for feature selection
"""
print(diagnosis)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import validation_curve
from sklearn.linear_model import Ridge

# Generate overfitting-prone data
np.random.seed(42)
X = np.random.randn(100, 20)
y = 2*X[:, 0] - X[:, 1] + np.random.randn(100) * 0.5

# Validation curve
alphas = np.logspace(-4, 4, 50)
train_scores, val_scores = validation_curve(
    Ridge(), X, y, param_name='alpha',
    param_range=alphas, cv=5, scoring='r2'
)

train_mean = train_scores.mean(axis=1)
val_mean = val_scores.mean(axis=1)

# Find optimal alpha
optimal_idx = np.argmax(val_mean)
optimal_alpha = alphas[optimal_idx]

plt.figure(figsize=(10, 6))
plt.semilogx(alphas, train_mean, 'b-', label='Training')
plt.semilogx(alphas, val_mean, 'r-', label='Validation')
plt.axvline(optimal_alpha, color='g', linestyle='--', 
            label=f'Optimal alpha={optimal_alpha:.4f}')
plt.xlabel('Alpha (Regularization Strength)')
plt.ylabel('R² Score')
plt.title('Validation Curve for Ridge Regression')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

print(f"Optimal alpha: {optimal_alpha:.6f}")
print(f"Best validation R²: {val_mean[optimal_idx]:.4f}")

import numpy as np
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import GridSearchCV, train_test_split
from sklearn.metrics import mean_squared_error, r2_score

# Generate non-linear data
np.random.seed(42)
X = np.random.randn(200, 3)
y = 2*X[:, 0]**2 - 3*X[:, 1] + X[:, 0]*X[:, 2] + np.random.randn(200) * 0.5

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# Build pipeline
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('poly', PolynomialFeatures()),
    ('regressor', ElasticNet(max_iter=10000))
])

# Hyperparameter grid
param_grid = {
    'poly__degree': [1, 2, 3],
    'regressor__alpha': [0.01, 0.1, 1.0],
    'regressor__l1_ratio': [0.2, 0.5, 0.8]
}

# Grid search with cross-validation
grid_search = GridSearchCV(pipeline, param_grid, cv=5, scoring='r2', n_jobs=-1)
grid_search.fit(X_train, y_train)

# Evaluate best model
best_model = grid_search.best_estimator_
y_pred = best_model.predict(X_test)

print("Best Parameters:", grid_search.best_params_)
print(f"Best CV R²: {grid_search.best_score_:.4f}")
print(f"Test R²: {r2_score(y_test, y_pred):.4f}")
print(f"Test RMSE: {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")

With three regularization techniques at your disposal, choosing the right one can be daunting. This section provides a decision framework based on your data characteristics, model requirements, and performance goals. The right choice depends on understanding your features and what you need from your model.

Decision Framework

Complete Comparison Table

Practical Comparison Example

import numpy as np
from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
import warnings
warnings.filterwarnings('ignore')

# Create dataset with different characteristics
np.random.seed(42)
n_samples = 200

# Scenario: Mix of important, irrelevant, and correlated features
X = np.random.randn(n_samples, 15)

# Add correlated features
X[:, 5] = X[:, 0] + np.random.randn(n_samples) * 0.1
X[:, 6] = X[:, 0] + np.random.randn(n_samples) * 0.1

# Target depends on only a few features
y = 3*X[:, 0] - 2*X[:, 1] + 1.5*X[:, 2] + np.random.randn(n_samples) * 0.5

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Compare models
models = {
    'Ridge (alpha=1.0)': Ridge(alpha=1.0),
    'Lasso (alpha=0.1)': Lasso(alpha=0.1),
    'Elastic Net (0.1, 0.5)': ElasticNet(alpha=0.1, l1_ratio=0.5)
}

print(f"{'Model':<25} {'CV R²':>10} {'Non-zero':>10}")
print("=" * 47)

for name, model in models.items():
    cv_score = cross_val_score(model, X_scaled, y, cv=5).mean()
    model.fit(X_scaled, y)
    n_nonzero = np.sum(model.coef_ != 0)
    print(f"{name:<25} {cv_score:>10.4f} {n_nonzero:>10}")

print("\nCoefficient Analysis:")
for name, model in models.items():
    model.fit(X_scaled, y)
    print(f"\n{name}:")
    print(f"  Top 3 coefficients: {np.argsort(np.abs(model.coef_))[-3:][::-1]}")
    print(f"  Coefficient values: {model.coef_[np.argsort(np.abs(model.coef_))[-3:][::-1]].round(3)}")