Ridge回帰#
モデル#
リッジ回帰（ridge regression）は誤差関数に正則化項を追加した線形回帰モデルである。
モデル自体は線形回帰と同様で目的変数\(y\)をパラメータ\(\boldsymbol{\beta}=(\beta_1, \beta_2, ..., \beta_d)^\top\)と特徴量\(\boldsymbol{x}=(x_1, x_2, ..., x_d)^\top\)の線形関数と誤差\(\varepsilon\)で表現するものになる。
サンプルサイズが\(n\)のデータセット\(\{\boldsymbol{x}_i, y_i\}^n_{i=1}\)があるとして、目的変数を\(\boldsymbol{y} = (y_1, y_2, ..., y_n)^\top\)、特徴量を\(\boldsymbol{X}=(\boldsymbol{x}_1, \boldsymbol{x}_2, ..., \boldsymbol{x}_n)^\top\)とおくと
と表記することもできる。
パラメータの推定#
制約付き最小化問題#
リッジ回帰が線形回帰と違う点は、パラメータの推定（誤差関数の最小化）に関して制約条件があること。
線形回帰は目的変数の実測値\(y_i\)と予測値\(\hat{y}_i=\sum^d_{j=1} x_{ij} \beta_j\)の誤差二乗和\(SSE=\sum^n_{i=1} (y_i - \hat{y}_i)^2\)を最小にするパラメータを求めるものであった。
リッジ回帰はこれに「パラメータ\(\beta_j\)の二乗和がある値\(R\)以下である」という制約条件が付いた下でパラメータを推定する。
この問題を解くために、ラグランジュの未定乗数法を利用する。
Show code cell source
# 二乗和が円になるイメージ
import numpy as np
import matplotlib.pyplot as plt
import itertools
R = 0.5
b1 = np.linspace(-1, 1, 100)
b2 = np.linspace(-1, 1, 100)
b = np.array(list(itertools.product(b1, b2)))
is_in_R = (b[:, 0]**2 + b[:, 1]**2) <= R
fig, ax = plt.subplots(figsize=[4, 4])
ax.scatter(b[~is_in_R, 0], b[~is_in_R, 1])
ax.scatter(b[is_in_R, 0], b[is_in_R, 1], color="orange")
Show code cell output
<matplotlib.collections.PathCollection at 0x7f1d0ccda9b0>
誤差関数の整理#
ラグランジュの未定乗数法を利用して、誤差関数を以下のように書くことができる。
この右辺第二項は正則化項、\(\lambda\)は正則化パラメータと呼ばれる。\(\lambda=0\)のときは最小二乗法と同じ誤差関数になり、得られる推定量も最小二乗推定量と等しくなる。
一般的なラグランジュ双対問題とは異なり、リッジ回帰では最適な\(\lambda\)を推定することはせず、あらかじめ\(\lambda\)を指定して\(x\)を推定する。（誤差を最小にする\(\lambda\)はゼロであり、正則化の意味がなくなるためだと思われる。）
また、最初の制約問題で登場した\(R\)は誤差関数に含めない（おそらく誤差関数を\(\beta\)について微分したときに定数項の\(R\)は消えて特に影響をもたらさないため）
ちなみに行列表記するとこんな感じに表される。
リッジ回帰の推定量\(\hat{\boldsymbol{\beta}}^{\text{Ridge}}\)は解析的に解くことができ、
なので、誤差関数の傾きがゼロ（誤差が最小値）の点を求めると
となり、
となる。
リッジ回帰の特徴#
リッジ回帰は制約をかけたことにより通常の線形回帰（最小二乗法）とは異なる特徴をもっている。
正則化#
線形回帰の最小二乗推定量\(\hat{\boldsymbol{\beta}}^{\text{LS}}\)は次のようなものだった。
このとき、
\(\boldsymbol{X}\)の列数が行数よりも多い（サンプルサイズより特徴量の次元数のほうが多い）
特徴量間の相関が非常に強い
といった状況においては、\(\boldsymbol{X}^\top \boldsymbol{X}\)が正則でなくなって逆行列が計算できなくなったり、あるいは推定が不安定になることがある。
そこでリッジ回帰の推定量\(\hat{\boldsymbol{\beta}}^{\text{Ridge}}\)では
と、\(\boldsymbol{X}^\top \boldsymbol{X}\)の対角成分に\(\lambda\)を足すことで尾根（ridge）を作り、正則にすることで逆行列が計算できるようにしている。
リッジ推定量が正則となる証明
\(X\)は\(n\)次の実行列とする。実対称行列\(X^\top X\)は非負値定符号行列であるため、
と分解可能。ここで\(P\)は直交行列であり、\(\Gamma = \mathrm{diag}(\gamma_1, \cdots, \gamma_n)\)は\(X^\top X\)の固有値（\(\gamma_1 \geq \cdots \geq \gamma_n \geq 0\)）を対角成分にもつ対角行列。
もし\(\gamma_n = 0\)なら\(X^\top X\)の逆行列は存在せず、\(\gamma_n > 0\)なら逆行列は存在し、
となる。\(X^\top X\)の最小固有値が\(\gamma_n \to 0\)の場合、\(1/\gamma_n \to \infty\)になり逆行列が計算できない。
一方、リッジ推定量のように\(X^\top X + \lambda I\)とする（\(\lambda \in\mathbb{R}\)）と、その逆行列は
となる。こちらは\(\gamma_n \to 0\)の場合であっても\(1 / (\gamma_p + \lambda I)\)は無限大に発散することがないため、\(X^\top X + \lambda I\)は正則となる。
過学習の抑制#
リッジ回帰の誤差関数
は、\(\lambda\)を大きくすると誤差に占める第二項の比重が大きくなり、パラメータ\(\beta_j\)の値をゼロに向けて縮小（shrink）させる。
パラメータが学習データに過学習して大きい値になっているような状況ではパラメータをうまく縮小させることによって過学習を抑制してモデルの汎化誤差を下げることができる。
推定量のバイアスとバリアンス#
正則化しない通常の最小二乗推定量はBLUE（最良線形不偏推定量：線形不偏推定量のなかでバリアンスが最小）だった。 リッジ回帰やlassoの推定量は不偏ではない（バイアスがある）推定量であるが、最小二乗推定量よりも小さなバリアンスとなる可能性がある。 そのため、最小二乗推定法を使用する線形回帰よりもリッジ回帰のほうが予測の2乗誤差を小さくする可能性がある。
Ridge Regression: Biased Estimation for Nonorthogonal Problems
モンテカルロ・シミュレーション#
多項式をつかったデータに、多項式回帰をfittingしてみる
import numpy as np
beta = np.array([3, 2, 1])
def gen_data(seed=0, n=500):
np.random.seed(seed)
x = np.random.uniform(-5, 3, size=n)
X = np.array([np.ones(n), x, x**2]).T
u = np.random.normal(scale=3, size=n)
y = X @ beta + u
return X, y
import matplotlib.pyplot as plt
X, y = gen_data(seed=0, n=30)
plt.scatter(X[:, 1], y)
<matplotlib.collections.PathCollection at 0x7f1d0cc38eb0>
Show code cell source
class OLS:
def fit(self, X, y):
self.beta_ = np.linalg.inv(X.T @ X) @ X.T @ y
return self
def predict(self, X):
return X @ self.beta_
class Ridge:
def __init__(self, lamb = 0.1):
self.lambda_ = lamb
def fit(self, X, y):
d = X.shape[1]
I = np.ones((d, d))
self.beta_ = np.linalg.inv((X.T @ X) + self.lambda_ * I) @ X.T @ y
return self
def predict(self, X):
return X @ self.beta_
from sklearn.metrics import mean_squared_error
results = []
for i in range(500):
result = {}
X_train, y_train = gen_data(seed=i, n=30)
X_test, y_test = gen_data(seed=i+1, n=100)
reg = OLS().fit(X_train, y_train)
y_pred = reg.predict(X_test)
result["ols_beta"] = reg.beta_[1]
result["ols_error"] = mean_squared_error(y_test, y_pred)
reg = Ridge(lamb=1).fit(X_train, y_train)
y_pred = reg.predict(X_test)
result["ridge_beta"] = reg.beta_[1]
result["ridge_error"] = mean_squared_error(y_test, y_pred)
results.append(result)
import pandas as pd
results = pd.DataFrame(results)
print(f"""
beta1
bias:
- ols: {beta[1] - results["ols_beta"].mean():.3f}
- ridge: {beta[1] - results["ridge_beta"].mean():.3f}
variance:
- ols: {results["ols_beta"].var():.3f}
- ridge: {results["ridge_beta"].var():.3f}
""")
beta1
bias:
- ols: -0.026
- ridge: 0.038
variance:
- ols: 0.115
- ridge: 0.110
ridgeのほうがvarianceが若干小さくbiasが大きくなっている
Show code cell source
fig, ax = plt.subplots()
ax.hist(beta[1] - results["ols_beta"], label="OLS", alpha=0.5, bins=20)
ax.hist(beta[1] - results["ridge_beta"], label="Ridge", alpha=0.5, bins=20)
ax.axvline(beta[1] - results["ols_beta"].mean(), label="Bias(OLS)", color="steelblue")
ax.axvline(beta[1] - results["ridge_beta"].mean(), label="Bias(Ridge)", color="darkorange")
ax.legend()
ax.set(xlabel=r"$\beta_1 - \hat\beta_1$", ylabel="Frequency")
fig.show()
fig, ax = plt.subplots()
ax.hist(results["ols_error"], label="OLS", alpha=0.5, bins=20)
ax.hist(results["ridge_error"], label="Ridge", alpha=0.5, bins=20)
ax.axvline(results["ols_error"].mean(), label="Mean(OLS)", color="steelblue")
ax.axvline(results["ridge_error"].mean(), label="Mean(Ridge)", color="darkorange")
ax.legend()
ax.set(xlabel="Mean Squared Error", ylabel="Frequency")
fig.show()
参考文献#
リッジ回帰#
杉山将『イラストで学ぶ機械学習』
スパース推定
ESLⅡ
佐和隆光. (1979). 回帰分析.
ラグランジュの未定乗数法#
『言語処理のための機械学習入門』pp.14-19
『これならわかる最適化数学』