능형회귀 보기¶

(1) 해석적으로. \(X = UDV^\top\)를 특이값분해라 하면 (연습문제 2에서 유도하듯)

이고 \(v_j\)가 정규직교이므로

이다. 각 항의 분모가 \(\lambda\)에 대해 증가하므로 모든 항이 따로따로 감소한다. 합도 감소한다. 여기에 U자가 생길 수 없으며, 이것은 자료와 무관한 항등식 수준의 사실이다.

참값과의 거리는 다르다. \(\lambda\)가 크면 \(\hat\beta \to 0\)이고 거리는 \(\lVert\beta\rVert = \sqrt{9+4+2.25} = 3.9051\)로 올라가므로, 거리가 끝까지 줄어들 수는 없다. (2)에서 보듯 중간에 바닥이 있다.

(2) 해석적으로. \(S = X^\top X\), \(A_\lambda = (S+\lambda I)^{-1}\)이라 두자. \(X\)를 고정하면 \(\hat\beta = A_\lambda X^\top y\)이고 \(y = X\beta + \varepsilon\), \(\operatorname{Var}(\varepsilon) = \sigma^2 I\)이므로

이다. 마지막 등식은 \(A_\lambda S - I = A_\lambda(S + \lambda I - \lambda I) - I = -\lambda A_\lambda\)에서 나온다. 분산 쪽은 \(\operatorname{Var}(\hat\beta) = \sigma^2 A_\lambda S A_\lambda\)이므로

다. \(\lambda = 0\)에서 값은 \(\sigma^2\operatorname{tr}(S^{-1})\)으로 최소제곱의 오차다.

이제 \(\lambda = 0\)에서의 도함수를 본다. 편향항은 \(\lambda^2\)에 비례하므로 미분하면 \(2\lambda(\cdots)\) 꼴이라 \(\lambda = 0\)에서 0이다. 분산항은 \(\frac{d}{d\lambda}A_\lambda = -A_\lambda^2\)이므로

이다. 편향은 이차로 자라고 분산은 일차로 줄어든다. 그러므로 \(\lambda = 0\) 바로 오른쪽에서 오차가 반드시 감소하며, 최소제곱보다 나은 \(\lambda > 0\)이 언제나 존재한다. 이것이 능형회귀의 존재정리이고, 이 보기의 표에서 \(\lambda = 1\)이 \(\lambda = 0\)을 이기는 까닭이다.

(2) 수치적으로.

import numpy as np
from sklearn.linear_model import LinearRegression, Ridge

rng = np.random.default_rng(42)

# 설명변수 20개 중 참으로 쓰이는 것은 앞의 셋뿐이다. 관측은 40개로 변수 수에
# 비해 적어, 최소제곱이 잡음까지 따라가기 좋은 상황이다.
n, p = 40, 20
X = rng.normal(size=(n, p))
beta_true = np.zeros(p)
beta_true[:3] = [3.0, -2.0, 1.5]
y = X @ beta_true + rng.normal(0, 1, n)

ols = LinearRegression().fit(X, y)

# lambda 를 키울수록 계수가 0 쪽으로 줄어든다. 라쏘와 달리 정확히 0 이 되지는
# 않고 작아지기만 한다. 그래서 능형은 변수를 고르지 못한다.
print(f"{'lambda':>8}  {'계수의 L2 크기':>14}  {'참값과의 거리':>14}")
print(f"{0.0:>8.1f}  {np.linalg.norm(ols.coef_):>14.3f}  "
      f"{np.linalg.norm(ols.coef_ - beta_true):>14.3f}")
for lam in [0.1, 1.0, 10.0, 100.0]:
    ridge = Ridge(alpha=lam).fit(X, y)
    print(f"{lam:>8.1f}  {np.linalg.norm(ridge.coef_):>14.3f}  "
          f"{np.linalg.norm(ridge.coef_ - beta_true):>14.3f}")

# 참값과의 거리가 lambda=0 일 때보다 중간 어딘가에서 작아진다. 편향을 조금
# 받아들이는 대가로 분산을 크게 줄인 결과이며, 이것이 편향-분산 절충이다.

# --- 유도한 기댓값 곡선과 맞춰 본다 (sigma = 1, X 는 고정) ---
S = X.T @ X

def expected_sq_error(lam):
    """E||beta_hat(lam) - beta||^2 = 편향^2 + 분산.  X 를 고정한 조건부 기댓값."""
    A = np.linalg.inv(S + lam * np.eye(p))
    bias2 = lam ** 2 * beta_true @ (A @ A) @ beta_true
    var = np.trace(A @ S @ A)
    return bias2, var

lams = np.logspace(-2, 2, 400)
tot = np.array([sum(expected_sq_error(l)) for l in lams])
k = tot.argmin()
b2, v = expected_sq_error(lams[k])
print(f"이론: lambda = 0 에서 E||.||^2 = tr(S^-1) = {np.trace(np.linalg.inv(S)):.4f}")
print(f"이론: 최솟값 {tot[k]:.4f} (편향^2 {b2:.4f} + 분산 {v:.4f}) at lambda = {lams[k]:.4f}")
print(f"0 에서의 도함수 = -2*tr(S^-2) = {-2 * np.trace(np.linalg.inv(S) @ np.linalg.inv(S)):.4f}  (< 0)")

grid = np.logspace(-2, 2, 2000)
dist = np.array([np.linalg.norm(Ridge(alpha=l).fit(X, y).coef_ - beta_true) for l in grid])
print(f"이 표본의 거리 곡선 최솟값 {dist.min():.4f} at lambda = {grid[dist.argmin()]:.4f}")

출력:

  lambda       계수의 L2 크기         참값과의 거리
     0.0           4.202           1.081
     0.1           4.169           1.048
     1.0           3.929           0.880
    10.0           2.891           1.404
   100.0           1.091           3.026
이론: lambda = 0 에서 E||.||^2 = tr(S^-1) = 1.3683
이론: 최솟값 1.0621 (편향^2 0.1914 + 분산 0.8707) at lambda = 1.6810
0 에서의 도함수 = -2*tr(S^-2) = -0.5116  (< 0)
이 표본의 거리 곡선 최솟값 0.8499 at lambda = 1.7584

유도와 코드가 맞는다. 기댓값 곡선은 \(\lambda = 1.6810\)에서 바닥을 치고, 이 한 표본의 거리 곡선은 \(\lambda = 1.7584\)에서 바닥을 친다. 둘이 정확히 같을 이유는 없다. 앞의 것은 \(y\)에 대한 기댓값이고 뒤의 것은 \(y\)를 한 번 뽑은 실현값이기 때문이다. 바닥의 위치가 가까운 것으로 충분하며, 바닥의 값은 \(1.0305\)(기댓값의 제곱근) 대 \(0.8499\)로 더 벌어진다. 이 표본이 운 좋게 평균보다 잘 맞은 것이다.

\(\lambda = 0\)에서의 도함수가 \(-0.5116\)으로 음수인 것이 (2)의 결론을 그대로 확인한다. 최적점에서 편향제곱 \(0.1914\)와 분산 \(0.8707\)의 합 \(1.0621\)이 최소제곱의 \(1.3683\)보다 \(22\%\) 작다. 편향 \(0.19\)를 사서 분산 \(0.50\)을 깎은 거래다.

위 표의 다섯 줄을 촘촘한 격자로 채워 그린 것이 위 그림이다. 파란 곡선은 계수벡터의 크기 \(\lVert\hat{\boldsymbol{\beta}}\rVert_2\), 주황 곡선은 참값과의 거리 \(\lVert\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}\rVert_2\)다. 두 곡선의 모양이 다르다는 것이 이 그림의 전부다.

파란 곡선은 끝까지 단조감소한다. \(4.202\)에서 출발해 \(\lambda = 100\)에서 \(1.091\)까지, 그리고 그 뒤로도 계속 내려간다. (1)에서 항별로 보인 그대로이며 여기에는 U자가 있을 수 없다.

주황 곡선은 다르다. \(\lambda = 0.01\)에서 \(1.081\)(OLS와 사실상 같다)로 시작해 \(\lambda = 1.76\)에서 \(0.850\)까지 내려갔다가 다시 올라간다. \(\lambda = 10\)에서 벌써 \(1.404\)로 OLS보다 나빠지고 \(\lambda = 100\)에서는 \(3.026\)이다. 회색 점선으로 그린 OLS 수준 \(1.081\)과 주황 곡선이 만나는 지점이 \(\lambda = 5.51\)인데, 그보다 큰 \(\lambda\)는 손해라는 뜻이다. 이득이 나는 구간은 생각보다 좁다.

여기서 실무적 함정 하나가 보인다. 계수가 작아지는 것은 눈에 잘 띄지만 정확도가 나빠지는 것은 눈에 띄지 않는다. \(\lambda = 100\)에서 계수의 크기는 \(1.091\)로 "아주 깔끔하게 정리된" 모형처럼 보이는데, 실제로는 참값에서 \(3.026\)만큼 떨어져 OLS보다 세 배 나쁘다. 계수가 작다는 것은 좋은 모형의 증거가 아니다. 어디가 바닥인지는 오직 자료로 추정해야 하며, 그 방법이 교차검증이다.

목적함수는

이다. 전개한 뒤 \(\beta\)로 미분하면

이고, 이를 0으로 놓으면

를 얻는다. \(\lambda > 0\)이면 \(X^\top X + \lambda I_p\)는 양정치이므로 가역이고, 따라서

\(X = U D V^\top\)이고 \(D = \text{diag}(d_1, \dots, d_p)\)라 하자. 그러면 \(X^\top X = V D^2 V^\top\), \(X^\top y = V D U^\top y\)이므로

이다. 성분으로 쓰면 \(V^\top \hat{\beta}^{\text{ridge}}\)의 \(j\)번째 원소가 \(\frac{d_j}{d_j^2 + \lambda} u_j^\top y\)이므로

를 얻는다. 인자 \(d_j^2 / (d_j^2 + \lambda) \in [0, 1)\)은 특이값이 작은 방향을 더 강하게 축소한다. \(\square\)

\(X^\top X = I_p\)이면 OLS 추정량은 \(\hat{\beta}^{\text{OLS}} = X^\top y\)이고, 능형 추정량은

이 된다. 즉 모든 계수가 \(1/(1 + \lambda)\)배로 균일하게 축소된다. 이는 능형회귀가 비례 축소를 수행함을 확인해 준다. \(\square\)

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

np.random.seed(42)
n, p = 200, 10

# 상관된 설계행렬
rho = 0.9
Sigma = rho * np.ones((p, p)) + (1 - rho) * np.eye(p)
L = np.linalg.cholesky(Sigma)
X = np.random.randn(n, p) @ L.T
beta_true = np.array([3, -2, 1.5, 0, 0, 0, 0, 0, 0, 0])
y = X @ beta_true + np.random.randn(n)

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

# OLS
ols_scores = cross_val_score(
    LinearRegression(), X_s, y, cv=5,
    scoring="neg_mean_squared_error"
)
ols_rmse = np.sqrt(-ols_scores.mean())

# 능형회귀의 격자탐색
best_rmse, best_lam = np.inf, None
for exp in range(-3, 4):
    lam = 10.0 ** exp
    scores = cross_val_score(
        Ridge(alpha=lam), X_s, y, cv=5,
        scoring="neg_mean_squared_error"
    )
    rmse = np.sqrt(-scores.mean())
    if rmse < best_rmse:
        best_rmse, best_lam = rmse, lam

print(f"OLS CV RMSE:  {ols_rmse:.4f}")
print(f"Best lambda:  {best_lam}")
print(f"Ridge CV RMSE: {best_rmse:.4f}")

출력:

OLS CV RMSE:  1.0171
Best lambda:  1.0
Ridge CV RMSE: 1.0156

실행하면 OLS의 교차검증 RMSE는 1.0171이고, 최적 \(\lambda = 1\)에서 능형회귀의 RMSE는 1.0156으로 조금 더 작다. \(\lambda\)를 더 키우면(10, 100, 1000) RMSE는 각각 1.0768, 1.4158, 1.7891로 오히려 나빠진다. 즉 정칙화는 도움이 되지만 그 이득의 크기와 최적 \(\lambda\)의 위치는 자료에 따라 다르며, 여기서는 \(n = 200\)이 \(p = 10\)에 비해 충분히 커서 OLS 자체가 이미 안정적이므로 개선폭이 작다. 개선폭은 \(p/n\)이 커질수록 뚜렷해진다. \(\square\)

KKT 조건에 의해 능형 문제는 어떤 \(t > 0\)에 대해

와 동치다. OLS 해는 아무 제약 없이 손실을 최소화하므로, \(\hat{\beta}^{\text{OLS}}\)는 제약영역 안에 있거나(이 경우 \(\hat{\beta}^{\text{ridge}} = \hat{\beta}^{\text{OLS}}\)이고 등호가 성립한다) 밖에 있다. 밖에 있으면 제약 최적해는 경계 \(\|\beta\|_2^2 = t < \|\hat{\beta}^{\text{OLS}}\|_2^2\) 위에 놓인다. 어느 경우든