6.5 Lab: 模型選擇與正則化（Lab: Model Selection & Regularization）

Theory

Formulas

Full Pipeline Code

Best Subset + Ridge CV + Lasso CV + 5-Fold CV + p>n Scenario

Complete end-to-end implementation: data generation, best subset search, Ridge/Lasso with cross-validated alpha, model comparison via 5-fold CV, and high-dimensional (p>n) stress test.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, Ridge, Lasso, RidgeCV, LassoCV
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
from itertools import combinations
import warnings
warnings.filterwarnings('ignore')

# 1. Generate synthetic data
np.random.seed(42)
n, p = 200, 30

# Only 5 features matter (sparse true model)
beta_true = np.zeros(p)
beta_true[[0, 1, 4, 9, 15]] = [3.0, -2.0, 1.5, 0.8, -1.2]

X = np.random.randn(n, p)
X[:, 5:] += 0.3 * X[:, :5].sum(axis=1, keepdims=True) / 5
y = X @ beta_true + np.random.randn(n) * 1.5

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=1)

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

print("Data:", X_train.shape, X_test.shape)
baseline = mean_squared_error(y_test, np.full_like(y_test, y_train.mean()))
print("Baseline MSE:", baseline)

# 2. Best Subset Selection (k=5)
print("\n--- Best Subset (k=5) ---")
best_score = float('inf')
best_combo = None
for combo in combinations(range(p), 5):
    lr = LinearRegression().fit(X_train_scaled[:, combo], y_train)
    mse = mean_squared_error(y_test, lr.predict(X_test_scaled[:, combo]))
    if mse < best_score:
        best_score = mse
        best_combo = combo

true_important = {0, 1, 4, 9, 15}
found = set(best_combo[:5])
print(f"Best 5-feature subset: {sorted(best_combo)}")
print(f"Test MSE: {best_score:.4f}")
print(f"True signal recall: {len(true_important & found)}/5")

# 3. Ridge Regression
print("\n--- Ridge ---")
ridge_cv = RidgeCV(alphas=np.logspace(-3, 3, 50))
ridge_cv.fit(X_train_scaled, y_train)
y_pred_ridge = ridge_cv.predict(X_test_scaled)
print(f"Best alpha: {ridge_cv.alpha_:.4f}")
print(f"Test MSE: {mean_squared_error(y_test, y_pred_ridge):.4f}")

# 4. Lasso Regression
print("\n--- Lasso ---")
lasso_cv = LassoCV(alphas=np.logspace(-3, 1, 50), cv=5,
                   max_iter=10000, random_state=42)
lasso_cv.fit(X_train_scaled, y_train)
y_pred_lasso = lasso_cv.predict(X_test_scaled)
nz = np.sum(lasso_cv.coef_ != 0)
selected = np.where(lasso_cv.coef_ != 0)[0]
print(f"Best alpha: {lasso_cv.alpha_:.4f}, Non-zero coefs: {nz}")
print(f"Selected: {sorted(selected)}")
recall = len(true_important & set(selected)) / len(true_important)
print(f"Recall: {recall:.0%}")
print(f"Test MSE: {mean_squared_error(y_test, y_pred_lasso):.4f}")

# 5. 5-Fold CV Comparison
print("\n--- 5-Fold CV ---")
models = {
    'OLS': LinearRegression(),
    'Ridge': Ridge(alpha=ridge_cv.alpha_),
    'Lasso': Lasso(alpha=lasso_cv.alpha_, max_iter=10000)
}
kf = KFold(n_splits=5, shuffle=True, random_state=42)
for name, model in models.items():
    scores = cross_val_score(model, X_train_scaled, y_train,
                             cv=kf, scoring='neg_mean_squared_error')
    print(f"{name}: CV MSE = {-scores.mean():.4f} (+/- {scores.std():.4f})")

# 6. High-Dimensional: p > n
print("\n--- High-Dim: p=100, n=50 ---")
np.random.seed(0)
beta_big = np.zeros(100)
beta_big[[0, 1, 2]] = [5.0, -3.0, 2.0]
Xb = np.random.randn(50, 100)
yb = Xb @ beta_big + np.random.randn(50) * 0.5
Xb_tr, Xb_te, yb_tr, yb_te = train_test_split(
    Xb, yb, test_size=15, random_state=1)
scb = StandardScaler()
Xb_tr_s = scb.fit_transform(Xb_tr)
Xb_te_s = scb.transform(Xb_te)

ridge_big = RidgeCV(alphas=np.logspace(-2, 2, 30)).fit(Xb_tr_s, yb_tr)
lasso_big = LassoCV(alphas=np.logspace(-2, 1, 30), cv=5,
                    max_iter=10000, random_state=42).fit(Xb_tr_s, yb_tr)

print(f"Ridge (p>n): MSE={mean_squared_error(yb_te, ridge_big.predict(Xb_te_s)):.4f}")
print(f"Lasso (p>n): MSE={mean_squared_error(yb_te, lasso_big.predict(Xb_te_s)):.4f}")
print(f"Lasso selected: {np.sum(lasso_big.coef_ != 0)}/100")
print("\nALL DONE")

Bonus: Lasso Coefficient Path Visualization

Show how true signal features (thick lines) survive regularization while noise features (thin gray) die at low alpha.

import matplotlib.pyplot as plt

alphas = np.logspace(-3, 1, 100)
coefs = []
for a in alphas:
    lasso = Lasso(alpha=a, max_iter=10000)
    lasso.fit(X_train_scaled, y_train)
    coefs.append(lasso.coef_)
coefs = np.array(coefs)

plt.figure(figsize=(12, 6))
for i in [0, 1, 4, 9, 15]:
    plt.plot(alphas, coefs[:, i], linewidth=2, label=f'Feature {i} (true)')
for i in [2, 3, 5, 6, 7, 8, 10, 11, 12]:
    plt.plot(alphas, coefs[:, i], alpha=0.3, linewidth=0.8, color='gray')
plt.xscale('log')
plt.axvline(lasso_cv.alpha_, color='red', linestyle='--',
            label=f'CV best alpha={lasso_cv.alpha_:.3f}')
plt.xlabel('Alpha (regularization strength)')
plt.ylabel('Coefficient value')
plt.title('Lasso Coefficient Paths: True Signals Survive Longer')
plt.legend(fontsize=8)
plt.tight_layout()
plt.show()

Method Comparison

Application Scenarios

Pros and Cons

Key Quote

Self-Reflection

For Hermes: This lesson mirrors our own agent architecture debates. Lasso-style "zero-out irrelevant features" maps to skill pruning (which skills actually matter?). Ridge-style "keep everything shrunk" maps to RAG retrieval (all memories contribute, but importance-weighted). The key insight from Ch6: there is no universally best method - the right choice depends entirely on your data structure and what you value (prediction vs interpretability, dense vs sparse truth). Same for agent memory systems: no single architecture dominates across all scenarios (as arXiv 2606.24775 also found this week).