Ridge Regression: L2 Regularization for Robust Linear Modeling

Table of Contents

1. Introduction and Motivation
2. Mathematical Foundation
3. Understanding Regularization
4. Ridge Regression vs Linear Regression
5. Hyperparameter Tuning
6. Feature Scaling and Ridge
7. Python Implementation from Scratch
8. Advanced Ridge Techniques
9. Two Comprehensive Case Studies
10. Multicollinearity and Ridge
11. Model Validation and Diagnostics
12. Comparative Analysis with Other Methods

---

1. Introduction and Motivation

Ridge regression, also known as L2 regularization or Tikhonov regularization, addresses one of the fundamental challenges in linear regression: overfitting due to multicollinearity and high-dimensional data. By adding a penalty term to the loss function, Ridge regression shrinks coefficients toward zero, leading to more stable and generalizable models.

Why Ridge Regression?

Traditional linear regression can suffer from several issues:

• Multicollinearity: When predictors are highly correlated, coefficients become unstable
• Overfitting: Models with many features can memorize training data
• Poor Generalization: High variance in coefficient estimates
• Numerical Instability: Ill-conditioned matrices in high-dimensional settings

Key Advantages

• Stability: Reduces coefficient variance and improves model stability
• Multicollinearity Handling: Effective for correlated features
• Generalization: Better performance on unseen data
• Bias-Variance Trade-off: Balances bias and variance optimally
• Mathematical Elegance: Well-understood theoretical foundation

Limitations to Consider

• Coefficient Shrinkage: All coefficients are shrunk, even important ones
• No Feature Selection: Doesn't perform automatic feature selection
• Hyperparameter Tuning: Requires careful selection of regularization strength
• Interpretability: Coefficients are biased (shrunken) estimates

---

2. Mathematical Foundation

Ridge Regression Objective Function

Ridge regression modifies the ordinary least squares (OLS) objective by adding an L2 penalty term:

minimize: ||y - Xβ||² + α||β||²

Where:

• y = target variable vector (n × 1)
• X = feature matrix (n × p)
• β = coefficient vector (p × 1)
• α = regularization parameter (α ≥ 0)
• ||β||² = L2 norm of coefficients (sum of squared coefficients)

Closed-Form Solution

The Ridge regression solution has a closed form:

β̂_ridge = (X^T X + αI)^(-1) X^T y

Where I is the identity matrix of size p × p.

Geometric Interpretation

The L2 penalty can be interpreted as:

1. Constraint Space: Restricts coefficients to lie within a hypersphere
2. Shrinkage: Pulls coefficients toward zero
3. Stability: Reduces the impact of multicollinearity

Bias-Variance Decomposition

Ridge regression introduces bias but reduces variance:

E[β̂_ridge] = (X^T X + αI)^(-1) X^T X β_true

As α increases:

• Bias increases: Coefficients are more biased toward zero
• Variance decreases: Coefficient estimates become more stable

---

3. Understanding Regularization

The Regularization Effect

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

# Set style and random seed
plt.style.use('seaborn-v0_8')
np.random.seed(42)

def demonstrate_regularization_effect():
    """Demonstrate how Ridge regularization affects coefficients and predictions"""
    
    # Generate synthetic data with multicollinearity
    n_samples = 100
    n_features = 10
    
    # Create correlated features
    X = np.random.normal(0, 1, (n_samples, n_features))
    
    # Introduce multicollinearity
    X[:, 2] = 0.8 * X[:, 0] + 0.2 * X[:, 1] + 0.1 * np.random.normal(0, 1, n_samples)
    X[:, 3] = 0.7 * X[:, 1] + 0.3 * X[:, 2] + 0.1 * np.random.normal(0, 1, n_samples)
    
    # True coefficients (only first 5 are important)
    true_coefficients = np.array([2.0, -1.5, 1.0, -0.8, 0.5, 0, 0, 0, 0, 0])
    
    # Generate target
    y = X @ true_coefficients + np.random.normal(0, 0.5, n_samples)
    
    # Scale features
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # Test different alpha values
    alphas = [0, 0.1, 1, 10, 100, 1000]
    
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    axes = axes.ravel()
    
    coefficient_data = []
    
    for i, alpha in enumerate(alphas):
        # Fit Ridge regression
        if alpha == 0:
            model = LinearRegression()
        else:
            model = Ridge(alpha=alpha)
        
        model.fit(X_scaled, y)
        
        # Get coefficients
        coefs = model.coef_
        coefficient_data.append(coefs)
        
        # Plot coefficients
        axes[i].bar(range(n_features), coefs, alpha=0.7, color='skyblue')
        axes[i].axhline(y=0, color='black', linestyle='-', alpha=0.3)
        axes[i].set_title(f'α = {alpha}\nCoefficient Magnitude')
        axes[i].set_xlabel('Feature Index')
        axes[i].set_ylabel('Coefficient Value')
        axes[i].grid(True, alpha=0.3)
        
        # Add true coefficients for comparison
        if alpha == 0:
            axes[i].plot(range(n_features), true_coefficients, 'ro', 
                        markersize=8, label='True Coefficients')
            axes[i].legend()
    
    plt.tight_layout()
    plt.show()
    
    # Coefficient shrinkage analysis
    print("Coefficient Shrinkage Analysis:")
    print("=" * 60)
    print(f"{'Alpha':<8} {'L2 Norm':<12} {'Max Coef':<12} {'Min Coef':<12} {'R²':<8}")
    print("-" * 60)
    
    for i, alpha in enumerate(alphas):
        coefs = coefficient_data[i]
        l2_norm = np.sqrt(np.sum(coefs**2))
        max_coef = np.max(np.abs(coefs))
        min_coef = np.min(np.abs(coefs))
        
        # Calculate R²
        if alpha == 0:
            model = LinearRegression()
        else:
            model = Ridge(alpha=alpha)
        
        model.fit(X_scaled, y)
        y_pred = model.predict(X_scaled)
        r2 = r2_score(y, y_pred)
        
        print(f"{alpha:<8} {l2_norm:<12.4f} {max_coef:<12.4f} {min_coef:<12.4f} {r2:<8.4f}")
    
    # Coefficient paths
    plt.figure(figsize=(12, 8))
    alphas_log = np.logspace(-3, 3, 50)
    coef_paths = []
    
    for alpha in alphas_log:
        ridge = Ridge(alpha=alpha)
        ridge.fit(X_scaled, y)
        coef_paths.append(ridge.coef_)
    
    coef_paths = np.array(coef_paths)
    
    for i in range(n_features):
        plt.plot(alphas_log, coef_paths[:, i], label=f'Feature {i+1}', linewidth=2)
    
    plt.xscale('log')
    plt.xlabel('Alpha (Regularization Strength)')
    plt.ylabel('Coefficient Value')
    plt.title('Ridge Coefficient Paths')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()

demonstrate_regularization_effect()

Bias-Variance Trade-off in Ridge Regression

def demonstrate_bias_variance_tradeoff():
    """Demonstrate bias-variance tradeoff with different alpha values"""
    
    # Generate synthetic data
    np.random.seed(42)
    n_samples = 100
    n_features = 20
    
    # Create features with some correlation
    X = np.random.normal(0, 1, (n_samples, n_features))
    
    # Introduce multicollinearity
    for i in range(5, n_features, 2):
        X[:, i] = 0.8 * X[:, i-1] + 0.2 * np.random.normal(0, 1, n_samples)
    
    # True coefficients (sparse)
    true_coefficients = np.zeros(n_features)
    true_coefficients[:5] = [2.0, -1.5, 1.0, -0.8, 0.5]
    
    # Generate target
    y = X @ true_coefficients + np.random.normal(0, 0.5, n_samples)
    
    # Scale features
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # Split data
    X_train, X_test, y_train, y_test = train_test_split(
        X_scaled, y, test_size=0.3, random_state=42
    )
    
    # Test different alpha values
    alphas = [0, 0.01, 0.1, 1, 10, 100]
    
    results = []
    
    for alpha in alphas:
        # Fit model
        if alpha == 0:
            model = LinearRegression()
        else:
            model = Ridge(alpha=alpha)
        
        model.fit(X_train, y_train)
        
        # Predictions
        y_train_pred = model.predict(X_train)
        y_test_pred = model.predict(X_test)
        
        # Metrics
        train_mse = mean_squared_error(y_train, y_train_pred)
        test_mse = mean_squared_error(y_test, y_test_pred)
        train_r2 = r2_score(y_train, y_train_pred)
        test_r2 = r2_score(y_test, y_test_pred)
        
        # Coefficient statistics
        coefs = model.coef_
        l2_norm = np.sqrt(np.sum(coefs**2))
        max_coef = np.max(np.abs(coefs))
        
        results.append({
            'alpha': alpha,
            'train_mse': train_mse,
            'test_mse': test_mse,
            'train_r2': train_r2,
            'test_r2': test_r2,
            'l2_norm': l2_norm,
            'max_coef': max_coef
        })
    
    # Plot results
    fig, axes = plt.subplots(2, 2, figsize=(15, 12))
    
    # MSE comparison
    alphas_plot = [r['alpha'] for r in results]
    train_mse_plot = [r['train_mse'] for r in results]
    test_mse_plot = [r['test_mse'] for r in results]
    
    axes[0, 0].plot(alphas_plot, train_mse_plot, 'bo-', label='Train MSE', linewidth=2)
    axes[0, 0].plot(alphas_plot, test_mse_plot, 'ro-', label='Test MSE', linewidth=2)
    axes[0, 0].set_xscale('log')
    axes[0, 0].set_xlabel('Alpha')
    axes[0, 0].set_ylabel('Mean Squared Error')
    axes[0, 0].set_title('MSE vs Alpha')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    
    # R² comparison
    train_r2_plot = [r['train_r2'] for r in results]
    test_r2_plot = [r['test_r2'] for r in results]
    
    axes[0, 1].plot(alphas_plot, train_r2_plot, 'bo-', label='Train R²', linewidth=2)
    axes[0, 1].plot(alphas_plot, test_r2_plot, 'ro-', label='Test R²', linewidth=2)
    axes[0, 1].set_xscale('log')
    axes[0, 1].set_xlabel('Alpha')
    axes[0, 1].set_ylabel('R² Score')
    axes[0, 1].set_title('R² vs Alpha')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    
    # L2 norm of coefficients
    l2_norm_plot = [r['l2_norm'] for r in results]
    
    axes[1, 0].plot(alphas_plot, l2_norm_plot, 'go-', linewidth=2)
    axes[1, 0].set_xscale('log')
    axes[1, 0].set_xlabel('Alpha')
    axes[1, 0].set_ylabel('L2 Norm of Coefficients')
    axes[1, 0].set_title('Coefficient Shrinkage')
    axes[1, 0].grid(True, alpha=0.3)
    
    # Maximum coefficient magnitude
    max_coef_plot = [r['max_coef'] for r in results]
    
    axes[1, 1].plot(alphas_plot, max_coef_plot, 'mo-', linewidth=2)
    axes[1, 1].set_xscale('log')
    axes[1, 1].set_xlabel('Alpha')
    axes[1, 1].set_ylabel('Max |Coefficient|')
    axes[1, 1].set_title('Maximum Coefficient Magnitude')
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Summary table
    print("Ridge Regression Performance Summary:")
    print("=" * 80)
    print(f"{'Alpha':<8} {'Train MSE':<12} {'Test MSE':<12} {'Train R²':<10} {'Test R²':<10} {'L2 Norm':<10}")
    print("-" * 80)
    
    for result in results:
        print(f"{result['alpha']:<8} {result['train_mse']:<12.4f} {result['test_mse']:<12.4f} "
              f"{result['train_r2']:<10.4f} {result['test_r2']:<10.4f} {result['l2_norm']:<10.4f}")

demonstrate_bias_variance_tradeoff()

---

4. Ridge Regression vs Linear Regression

When to Use Ridge Regression

Ridge regression is particularly effective in these scenarios:

1. Multicollinearity: When features are highly correlated
2. High-dimensional data: When p > n (more features than samples)
3. Overfitting: When linear regression shows poor generalization
4. Numerical instability: When OLS fails due to singular matrices

Comparison with Linear Regression

def compare_ridge_vs_linear():
    """Compare Ridge regression with Linear regression on different scenarios"""
    
    np.random.seed(42)
    
    # Scenario 1: Multicollinearity
    print("Scenario 1: Multicollinearity")
    print("=" * 50)
    
    n_samples = 50
    n_features = 5
    
    # Create highly correlated features
    X = np.random.normal(0, 1, (n_samples, n_features))
    X[:, 1] = 0.95 * X[:, 0] + 0.05 * np.random.normal(0, 1, n_samples)
    X[:, 2] = 0.9 * X[:, 0] + 0.1 * np.random.normal(0, 1, n_samples)
    
    # True coefficients
    true_coefs = np.array([2.0, -1.5, 1.0, 0.5, -0.8])
    y = X @ true_coefs + np.random.normal(0, 0.5, n_samples)
    
    # Scale features
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # Fit models
    linear_reg = LinearRegression()
    ridge_reg = Ridge(alpha=1.0)
    
    linear_reg.fit(X_scaled, y)
    ridge_reg.fit(X_scaled, y)
    
    # Compare coefficients
    print("Coefficient Comparison:")
    print(f"{'Feature':<8} {'True':<8} {'Linear':<8} {'Ridge':<8}")
    print("-" * 40)
    
    for i in range(n_features):
        print(f"{i+1:<8} {true_coefs[i]:<8.3f} {linear_reg.coef_[i]:<8.3f} {ridge_reg.coef_[i]:<8.3f}")
    
    # Compare predictions
    y_linear_pred = linear_reg.predict(X_scaled)
    y_ridge_pred = ridge_reg.predict(X_scaled)
    
    linear_mse = mean_squared_error(y, y_linear_pred)
    ridge_mse = mean_squared_error(y, y_ridge_pred)
    
    print(f"\nMSE Comparison:")
    print(f"Linear Regression MSE: {linear_mse:.4f}")
    print(f"Ridge Regression MSE: {ridge_mse:.4f}")
    print(f"Improvement: {((linear_mse - ridge_mse) / linear_mse * 100):.2f}%")
    
    # Scenario 2: High-dimensional data
    print(f"\nScenario 2: High-dimensional data (p > n)")
    print("=" * 50)
    
    n_samples = 30
    n_features = 50
    
    # Generate high-dimensional data
    X_high = np.random.normal(0, 1, (n_samples, n_features))
    
    # True coefficients (sparse)
    true_coefs_high = np.zeros(n_features)
    true_coefs_high[:10] = np.random.normal(0, 2, 10)
    
    y_high = X_high @ true_coefs_high + np.random.normal(0, 0.5, n_samples)
    
    # Scale features
    scaler_high = StandardScaler()
    X_high_scaled = scaler_high.fit_transform(X_high)
    
    # Split data
    X_train, X_test, y_train, y_test = train_test_split(
        X_high_scaled, y_high, test_size=0.3, random_state=42
    )
    
    # Fit models
    linear_reg_high = LinearRegression()
    ridge_reg_high = Ridge(alpha=10.0)
    
    linear_reg_high.fit(X_train, y_train)
    ridge_reg_high.fit(X_train, y_train)
    
    # Compare performance
    y_train_linear = linear_reg_high.predict(X_train)
    y_test_linear = linear_reg_high.predict(X_test)
    y_train_ridge = ridge_reg_high.predict(X_train)
    y_test_ridge = ridge_reg_high.predict(X_test)
    
    print("Performance Comparison:")
    print(f"{'Model':<20} {'Train MSE':<12} {'Test MSE':<12} {'Train R²':<10} {'Test R²':<10}")
    print("-" * 70)
    
    linear_train_mse = mean_squared_error(y_train, y_train_linear)
    linear_test_mse = mean_squared_error(y_test, y_test_linear)
    linear_train_r2 = r2_score(y_train, y_train_linear)
    linear_test_r2 = r2_score(y_test, y_test_linear)
    
    ridge_train_mse = mean_squared_error(y_train, y_train_ridge)
    ridge_test_mse = mean_squared_error(y_test, y_test_ridge)
    ridge_train_r2 = r2_score(y_train, y_train_ridge)
    ridge_test_r2 = r2_score(y_test, y_test_ridge)
    
    print(f"{'Linear Regression':<20} {linear_train_mse:<12.4f} {linear_test_mse:<12.4f} "
          f"{linear_train_r2:<10.4f} {linear_test_r2:<10.4f}")
    print(f"{'Ridge Regression':<20} {ridge_train_mse:<12.4f} {ridge_test_mse:<12.4f} "
          f"{ridge_train_r2:<10.4f} {ridge_test_r2:<10.4f}")
    
    # Coefficient stability analysis
    print(f"\nCoefficient Stability Analysis:")
    print(f"Linear Regression - Max |coef|: {np.max(np.abs(linear_reg_high.coef_)):.4f}")
    print(f"Ridge Regression - Max |coef|: {np.max(np.abs(ridge_reg_high.coef_)):.4f}")
    print(f"Linear Regression - L2 norm: {np.sqrt(np.sum(linear_reg_high.coef_**2)):.4f}")
    print(f"Ridge Regression - L2 norm: {np.sqrt(np.sum(ridge_reg_high.coef_**2)):.4f}")

compare_ridge_vs_linear()

---

5. Hyperparameter Tuning

Selecting the Optimal Alpha

def hyperparameter_tuning_demo():
    """Demonstrate hyperparameter tuning for Ridge regression"""
    
    from sklearn.model_selection import GridSearchCV, cross_val_score
    from sklearn.pipeline import Pipeline
    
    np.random.seed(42)
    
    # Generate synthetic data
    n_samples = 200
    n_features = 15
    
    # Create features with some correlation
    X = np.random.normal(0, 1, (n_samples, n_features))
    
    # Introduce multicollinearity
    for i in range(5, n_features, 2):
        X[:, i] = 0.8 * X[:, i-1] + 0.2 * np.random.normal(0, 1, n_samples)
    
    # True coefficients (sparse)
    true_coefficients = np.zeros(n_features)
    true_coefficients[:8] = [2.0, -1.5, 1.0, -0.8, 0.5, 1.2, -0.9, 0.7]
    
    # Generate target
    y = X @ true_coefficients + np.random.normal(0, 0.5, n_samples)
    
    # Split data
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=42
    )
    
    # Method 1: Grid Search with Cross-Validation
    print("Method 1: Grid Search with Cross-Validation")
    print("=" * 60)
    
    # Define parameter grid
    param_grid = {
        'ridge__alpha': [0.001, 0.01, 0.1, 1, 10, 100, 1000]
    }
    
    # Create pipeline
    pipeline = Pipeline([
        ('scaler', StandardScaler()),
        ('ridge', Ridge())
    ])
    
    # Grid search
    grid_search = GridSearchCV(
        pipeline, param_grid, cv=5, scoring='neg_mean_squared_error',
        return_train_score=True
    )
    
    grid_search.fit(X_train, y_train)
    
    print(f"Best alpha: {grid_search.best_params_['ridge__alpha']}")
    print(f"Best CV score: {-grid_search.best_score_:.4f}")
    
    # Method 2: Validation Curve
    print(f"\nMethod 2: Validation Curve Analysis")
    print("=" * 60)
    
    from sklearn.model_selection import validation_curve
    
    # Create pipeline for validation curve
    pipe = Pipeline([
        ('scaler', StandardScaler()),
        ('ridge', Ridge())
    ])
    
    # Define alpha range
    alphas = np.logspace(-3, 3, 20)
    
    # Compute validation curve
    train_scores, val_scores = validation_curve(
        pipe, X_train, y_train, param_name='ridge__alpha', 
        param_range=alphas, cv=5, scoring='neg_mean_squared_error'
    )
    
    # Convert to positive MSE
    train_scores = -train_scores
    val_scores = -val_scores
    
    # Plot validation curve
    plt.figure(figsize=(12, 8))
    
    plt.subplot(2, 2, 1)
    plt.plot(alphas, train_scores.mean(axis=1), 'bo-', label='Train MSE', linewidth=2)
    plt.plot(alphas, val_scores.mean(axis=1), 'ro-', label='Validation MSE', linewidth=2)
    plt.fill_between(alphas, train_scores.mean(axis=1) - train_scores.std(axis=1),
                     train_scores.mean(axis=1) + train_scores.std(axis=1), alpha=0.3)
    plt.fill_between(alphas, val_scores.mean(axis=1) - val_scores.std(axis=1),
                     val_scores.mean(axis=1) + val_scores.std(axis=1), alpha=0.3)
    plt.xscale('log')
    plt.xlabel('Alpha')
    plt.ylabel('Mean Squared Error')
    plt.title('Validation Curve')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    # Method 3: Learning Curves
    print(f"\nMethod 3: Learning Curves")
    print("=" * 60)
    
    from sklearn.model_selection import learning_curve
    
    # Best model from grid search
    best_alpha = grid_search.best_params_['ridge__alpha']
    best_pipeline = Pipeline([
        ('scaler', StandardScaler()),
        ('ridge', Ridge(alpha=best_alpha))
    ])
    
    # Compute learning curves
    train_sizes = np.linspace(0.1, 1.0, 10)
    train_sizes_abs, train_scores_lc, val_scores_lc = learning_curve(
        best_pipeline, X_train, y_train, train_sizes=train_sizes, cv=5,
        scoring='neg_mean_squared_error'
    )
    
    # Convert to positive MSE
    train_scores_lc = -train_scores_lc
    val_scores_lc = -val_scores_lc
    
    plt.subplot(2, 2, 2)
    plt.plot(train_sizes_abs, train_scores_lc.mean(axis=1), 'bo-', label='Train MSE', linewidth=2)
    plt.plot(train_sizes_abs, val_scores_lc.mean(axis=1), 'ro-', label='Validation MSE', linewidth=2)
    plt.fill_between(train_sizes_abs, train_scores_lc.mean(axis=1) - train_scores_lc.std(axis=1),
                     train_scores_lc.mean(axis=1) + train_scores_lc.std(axis=1), alpha=0.3)
    plt.fill_between(train_sizes_abs, val_scores_lc.mean(axis=1) - val_scores_lc.std(axis=1),
                     val_scores_lc.mean(axis=1) + val_scores_lc.std(axis=1), alpha=0.3)
    plt.xlabel('Training Set Size')
    plt.ylabel('Mean Squared Error')
    plt.title('Learning Curves')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    # Method 4: Coefficient Paths
    plt.subplot(2, 2, 3)
    
    alphas_path = np.logspace(-3, 3, 50)
    coef_paths = []
    
    for alpha in alphas_path:
        ridge = Ridge(alpha=alpha)
        ridge.fit(StandardScaler().fit_transform(X_train), y_train)
        coef_paths.append(ridge.coef_)
    
    coef_paths = np.array(coef_paths)
    
    for i in range(min(10, n_features)):
        plt.plot(alphas_path, coef_paths[:, i], label=f'Feature {i+1}', linewidth=2)
    
    plt.xscale('log')
    plt.xlabel('Alpha')
    plt.ylabel('Coefficient Value')
    plt.title('Coefficient Paths')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, alpha=0.3)
    
    # Method 5: Final Model Performance
    plt.subplot(2, 2, 4)
    
    # Fit best model
    best_pipeline.fit(X_train, y_train)
    y_train_pred = best_pipeline.predict(X_train)
    y_test_pred = best_pipeline.predict(X_test)
    
    # Plot predictions vs actual
    plt.scatter(y_train, y_train_pred, alpha=0.6, label='Train', s=50)
    plt.scatter(y_test, y_test_pred, alpha=0.6, label='Test', s=50)
    plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=2)
    plt.xlabel('Actual Values')
    plt.ylabel('Predicted Values')
    plt.title('Predictions vs Actual')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Final performance summary
    print(f"\nFinal Model Performance:")
    print("=" * 60)
    print(f"Best Alpha: {best_alpha}")
    print(f"Train MSE: {mean_squared_error(y_train, y_train_pred):.4f}")
    print(f"Test MSE: {mean_squared_error(y_test, y_test_pred):.4f}")
    print(f"Train R²: {r2_score(y_train, y_train_pred):.4f}")
    print(f"Test R²: {r2_score(y_test, y_test_pred):.4f}")
    
    # Compare with linear regression
    linear_pipeline = Pipeline([
        ('scaler', StandardScaler()),
        ('linear', LinearRegression())
    ])
    
    linear_pipeline.fit(X_train, y_train)
    y_train_linear = linear_pipeline.predict(X_train)
    y_test_linear = linear_pipeline.predict(X_test)
    
    print(f"\nComparison with Linear Regression:")
    print(f"Linear Train MSE: {mean_squared_error(y_train, y_train_linear):.4f}")
    print(f"Linear Test MSE: {mean_squared_error(y_test, y_test_linear):.4f}")
    print(f"Linear Train R²: {r2_score(y_train, y_train_linear):.4f}")
    print(f"Linear Test R²: {r2_score(y_test, y_test_linear):.4f}")

hyperparameter_tuning_demo()

---

Summary and Best Practices

Key Takeaways

Ridge Regression Strengths:

• Effectively handles multicollinearity
• Reduces overfitting through coefficient shrinkage
• Maintains all features (no feature selection)
• Provides stable coefficient estimates
• Works well with high-dimensional data

Common Challenges:

• Hyperparameter tuning: Requires careful selection of alpha
• Feature scaling: Sensitive to feature scales
• No feature selection: All coefficients are shrunk
• Interpretation: Coefficients are biased estimates

Best Practices:

1. Data Preprocessing

   • Always scale features before Ridge regression
   • Handle missing values appropriately
   • Check for multicollinearity

2. Hyperparameter Tuning

   • Use cross-validation to select optimal alpha
   • Consider log-scale search for alpha
   • Monitor both train and validation performance

3. Model Evaluation

   • Use multiple metrics (MSE, R², MAE)
   • Check for overfitting using validation curves
   • Examine coefficient stability

4. Feature Engineering

   • Create meaningful interaction terms
   • Consider polynomial features when appropriate
   • Remove highly correlated features if interpretability is important

5. Alternative Approaches

   • Use Lasso for automatic feature selection
   • Consider Elastic Net for balanced regularization
   • Try other regularization methods (L1, L2, or both)

Decision Framework:

• Use Linear Regression when data is well-conditioned and no multicollinearity
• Use Ridge Regression when multicollinearity is present or for high-dimensional data
• Use Lasso when feature selection is important
• Use Elastic Net when you want both regularization benefits
• Use Other Methods (Random Forest, Neural Networks) for complex non-linear relationships

This comprehensive guide provides both theoretical understanding and practical implementation skills for Ridge regression, enabling you to apply these techniques effectively to real-world problems while avoiding common pitfalls.