Ridge回帰

論文：Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55-67.

モデル

リッジ回帰（ridge regression）は誤差関数に正則化項を追加した線形回帰モデルである。

モデル自体は線形回帰と同様で目的変数\(y\)をパラメータ\(\boldsymbol{\beta}=(\beta_1, \beta_2, ..., \beta_d)^\top\)と特徴量\(\boldsymbol{x}=(x_1, x_2, ..., x_d)^\top\)の線形関数と誤差\(\varepsilon\)で表現するものになる。

\[ y = \beta_0 + \beta_1 x_1 + \cdots + \beta_d x_d + \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\)とおくと

\[ \boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{\varepsilon} \]

と表記することもできる。

パラメータの推定

制約付き最小化問題

リッジ回帰が線形回帰と違う点は、パラメータの推定（誤差関数の最小化）に関して制約条件があること。

線形回帰は目的変数の実測値\(y_i\)と予測値\(\hat{y}_i=\sum^d_{j=1} x_{ij} \beta_j\)の誤差二乗和\(SSE=\sum^n_{i=1} (y_i - \hat{y}_i)^2\)を最小にするパラメータを求めるものであった。

\[ \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \hat{\boldsymbol{\beta}}^{\text{LS}} = \argmin_{\boldsymbol{\beta}} \sum^n_{i=1} (y_i - \hat{y}_i)^2 \]

リッジ回帰はこれに「パラメータ\(\beta_j\)の二乗和がある値\(R\)以下である」という制約条件が付いた下でパラメータを推定する。

\[\begin{split} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \begin{align} \hat{\boldsymbol{\beta}}^{\text{Ridge}} = \text{ } & \argmin_{\boldsymbol{\beta}} \sum^n_{i=1} (y_i - \hat{y}_i)^2 \\ & \text{subject to } \sum^d_{i=1} \beta_j^2 \leq R \end{align} \end{split}\]

この制約条件は円の方程式と呼ばれるもので、2次元で描くと円の形になる。円の範囲内で誤差を最小化するパラメータを選ぶという条件付き最適化問題を解くことになる。

この問題を解くために、ラグランジュの未定乗数法を利用する。

Show code cell source

Hide code cell source

# 二乗和が円になるイメージ
import numpy as np
import matplotlib.pyplot as plt
import itertools

R = 0.5
b1 = np.linspace(-1, 1, 300)
b2 = np.linspace(-1, 1, 300)
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], color="lightblue", label=r"$\beta_1^2 + \beta_2^2 \leq 0.5$")
ax.axvline(color="black", linewidth=1)
ax.axhline(color="black", linewidth=1)
ax.set(xlabel=r"$\beta_1$", ylabel=r"$\beta_2$", xlim=(-1,1), ylim=(-1,1))
ax.legend()
fig.show()

Show code cell output

Hide code cell output

誤差関数の整理

ラグランジュの未定乗数法を利用して、誤差関数を以下のように書くことができる。

\[ J({\boldsymbol{\beta}}) = \sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^d \beta_j^2 \]

この右辺第二項は正則化項、\(\lambda\)は正則化パラメータと呼ばれる。\(\lambda=0\)のときは最小二乗法と同じ誤差関数になり、得られる推定量も最小二乗推定量と等しくなる。

一般的なラグランジュ双対問題とは異なり、リッジ回帰では最適な\(\lambda\)を推定することはせず、あらかじめ\(\lambda\)を指定して\(x\)を推定する。（誤差を最小にする\(\lambda\)はゼロであり、正則化の意味がなくなるためだと思われる。）

また、最初の制約問題で登場した\(R\)は誤差関数に含めない（おそらく誤差関数を\(\beta\)について微分したときに定数項の\(R\)は消えて特に影響をもたらさないため）

ちなみに行列表記するとこんな感じに表される。

\[\begin{split} \begin{align} J({\boldsymbol{\beta}}) &= \sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^d \beta_j^2 \\ &= \| \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\|^2_2 + \lambda\| \boldsymbol{\beta}\|^2_2 \\ &= (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})^\top (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}) + \lambda \boldsymbol{\beta}^\top \boldsymbol{\beta} \end{align} \end{split}\]

リッジ回帰の推定量\(\hat{\boldsymbol{\beta}}^{\text{Ridge}}\)は解析的に解くことができ、

\[\begin{split} \begin{align} J(\boldsymbol{\beta}) &= (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})^\top (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}) + \lambda \boldsymbol{\beta}^\top \boldsymbol{\beta} \\ &= \boldsymbol{y}^\top \boldsymbol{y} - \boldsymbol{y}^\top \boldsymbol{X}\boldsymbol{\beta} - (\boldsymbol{X}\boldsymbol{\beta})^\top \boldsymbol{y} + (\boldsymbol{X}\boldsymbol{\beta})^\top (\boldsymbol{X}\boldsymbol{\beta}) + \lambda \boldsymbol{\beta}^\top \boldsymbol{\beta} \\ &= \boldsymbol{y}^\top \boldsymbol{y} - 2 \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{y} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} + \lambda \boldsymbol{\beta}^\top \boldsymbol{\beta} \end{align} \end{split}\]

なので、誤差関数の傾きがゼロ（誤差が最小値）の点を求めると

\[\begin{split} \frac{\partial J(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2\boldsymbol{X}^\top\boldsymbol{y} + 2(\boldsymbol{X}^\top\boldsymbol{X})\boldsymbol{\beta} + 2\lambda\boldsymbol{\beta} =\boldsymbol{0} \\ \Longrightarrow (\boldsymbol{X}^\top\boldsymbol{X})\boldsymbol{\beta} + \lambda\boldsymbol{\beta} =\boldsymbol{X}^\top\boldsymbol{y} \\ \Longrightarrow (\boldsymbol{X}^\top\boldsymbol{X} + \lambda \boldsymbol{I})\boldsymbol{\beta} =\boldsymbol{X}^\top\boldsymbol{y} \end{split}\]

となり、

\[ \hat{\boldsymbol{\beta}}^{\text{Ridge}} = (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{y} \]

となる。

リッジ回帰の特徴

リッジ回帰は制約をかけたことにより通常の線形回帰（最小二乗法）とは異なる特徴をもっている。

正則化

線形回帰の最小二乗推定量\(\hat{\boldsymbol{\beta}}^{\text{LS}}\)は次のようなものだった。

\[ \hat{\boldsymbol{\beta}}^{\text{LS}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y} \]

このとき、

1. \(\boldsymbol{X}\)の列数が行数よりも多い（サンプルサイズより特徴量の次元数のほうが多い）

2. 特徴量間の相関が非常に強い

といった状況においては、\(\boldsymbol{X}^\top \boldsymbol{X}\)が正則でなくなって逆行列が計算できなくなったり、あるいは推定が不安定になることがある。

そこでリッジ回帰の推定量\(\hat{\boldsymbol{\beta}}^{\text{Ridge}}\)では

\[ \hat{\boldsymbol{\beta}}^{\text{Ridge}} = (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{y} \]

と、\(\boldsymbol{X}^\top \boldsymbol{X}\)の対角成分に\(\lambda\)を足すことで尾根（ridge）を作り、正則にすることで逆行列が計算できるようにしている。

リッジ推定量が正則となる証明

\(X\)は\(n\)次の実行列とする。実対称行列\(X^\top X\)は非負値定符号行列であるため、

\[ X^\top X = P \Gamma P^\top \]

と分解可能。ここで\(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\)なら逆行列は存在し、

\[\begin{split} \begin{aligned} (X^\top X)^{-1} &= (P \Gamma P^\top)^{-1}\\ &= P \Gamma^{-1} P^\top \quad (\because (ABC)^{-1} = (C^{-1} B^{-1} A^{-1}) であり、 Pは直交行列なのでP^\top = P^{-1})\\ &= P \operatorname{diag}(1 / \gamma_1, \ldots, 1 / \gamma_n) P^\top \end{aligned} \end{split}\]

となる。\(X^\top X\)の最小固有値が\(\gamma_n \to 0\)の場合、\(1/\gamma_n \to \infty\)になり逆行列が計算できない。

一方、リッジ推定量のように\(X^\top X + \lambda I\)とする（\(\lambda \in\mathbb{R}\)）と、その逆行列は

\[\begin{split} \begin{aligned} (X^\top X + \lambda I)^{-1} &= (P \Gamma P^\top + \lambda I)^{-1}\\ &= \{ P (\Gamma + \lambda I) P^\top \}^{-1}\\ &= P (\Gamma + \lambda I)^{-1} P^\top\\ &= P \operatorname{diag}[1 / (\gamma_1 + \lambda I), \ldots, 1 / (\gamma_n + \lambda I)] P^\top \end{aligned} \end{split}\]

となる。こちらは\(\gamma_n \to 0\)の場合であっても\(1 / (\gamma_p + \lambda I)\)は無限大に発散することがないため、\(X^\top X + \lambda I\)は正則となる。

過学習の抑制

リッジ回帰の誤差関数

\[ J(\beta)= \sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^d \beta_j^2 \]

は、\(\lambda\)を大きくすると誤差に占める第二項の比重が大きくなり、パラメータ\(\beta_j\)の値をゼロに向けて縮小（shrink）させる。

パラメータが学習データに過学習して大きい値になっているような状況ではパラメータをうまく縮小させることによって過学習を抑制してモデルの汎化誤差を下げることができる。

OLS推定量との関係

\(\boldsymbol{A}:=\boldsymbol{X}^\top \boldsymbol{X}\)とおくと

\[\begin{split} \begin{aligned} \hat{\boldsymbol{\beta}}^{\text{Ridge}} &= (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{y}\\ &= (\boldsymbol{A} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{y}\\ \end{aligned} \end{split}\]

となる。このとき

\[ \boldsymbol{A} + \lambda \boldsymbol{I} = \boldsymbol{A} (\boldsymbol{I} + \lambda \boldsymbol{A}^{-1}) \]

であり、この逆行列をとると\((AB)^{-1}=B^{-1}A^{-1}\)より

\[ (\boldsymbol{A} + \lambda \boldsymbol{I})^{-1} = [\boldsymbol{A} (\boldsymbol{I} + \lambda \boldsymbol{A}^{-1})]^{-1} = (\boldsymbol{I} + \lambda \boldsymbol{A}^{-1})^{-1} \boldsymbol{A}^{-1} \]

よって

\[\begin{split} \begin{aligned} \hat{\boldsymbol{\beta}}^{\text{Ridge}} &= (\boldsymbol{I} + \lambda \boldsymbol{A}^{-1})^{-1} \boldsymbol{A}^{-1} \boldsymbol{X}^\top \boldsymbol{y} \\ &= [\boldsymbol{I} + \lambda (\boldsymbol{X}^\top \boldsymbol{X})^{-1}]^{-1} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y} \\ &= [\boldsymbol{I} + \lambda (\boldsymbol{X}^\top \boldsymbol{X})^{-1}]^{-1} \hat{\boldsymbol{\beta}}^{\text{OLS}} \\ &= \boldsymbol{Z} \hat{\boldsymbol{\beta}}^{\text{OLS}} \end{aligned} \end{split}\]

と変形することができる。ここで

\[ \boldsymbol{Z} := [\boldsymbol{I} + \lambda (\boldsymbol{X}^\top \boldsymbol{X})^{-1}]^{-1} \]

とおいた。つまりRidge推定量はOLS推定量を\(Z\)で変換したもの。

なお\(\lambda=0\)のときにRidge推定量とOLS推定量が一致することは以下のようにわかる。

\[\begin{split} \begin{aligned} \hat{\boldsymbol{\beta}}^{\text{Ridge}} &= [\boldsymbol{I} + \underbrace{ \lambda (\boldsymbol{X}^\top \boldsymbol{X})^{-1} }_{=0} ]^{-1} \hat{\boldsymbol{\beta}}^{\text{OLS}} \\ &= \boldsymbol{I}^{-1} \hat{\boldsymbol{\beta}}^{\text{OLS}} \\ &= \hat{\boldsymbol{\beta}}^{\text{OLS}} \quad (\because I^{-1} = I)\\ \end{aligned} \end{split}\]

推定量のバイアスとバリアンス

正則化しない通常の最小二乗推定量はBLUE（最良線形不偏推定量：線形不偏推定量のなかでバリアンスが最小）だった。 リッジ回帰やlassoの推定量は不偏ではない（バイアスがある）推定量であるが、最小二乗推定量よりも小さなバリアンスとなる可能性がある。 そのため、最小二乗推定法を使用する線形回帰よりもリッジ回帰のほうが予測の2乗誤差を小さくする可能性がある。

バイアス・バリアンスのトレードオフ

Ridge推定量のバイアスとバリアンスは次のようになっている（鈴木 (2018)）

\[\begin{split} \begin{aligned} \mathrm{Bias} &= \lambda^2 \| ( X^\top X + \lambda I)^{-1} \boldsymbol{\beta} \|^2 \\ \mathrm{Variance} &= \sigma^2 \operatorname{Tr}\left[\left(X^{\top} X+\lambda I\right)^{-2} X^{\top} X\right] \end{aligned} \end{split}\]

これらにはトレードオフの関係がある。

(1) \(\lambda \to 0\)のとき（OLSに近づくとき）、\(\mathrm{Bias}\to 0, \mathrm{Variance} \to \sigma^2 \operatorname{Tr}[(X^{\top} X)^{-1}]\)となる。

\[\begin{split} \begin{aligned} \lim_{\lambda \to 0} \mathrm{Bias} &= 0 \times \| ( X^\top X )^{-1} \boldsymbol{\beta} \|^2 \to 0\\ \lim_{\lambda \to 0} \mathrm{Variance} &= \sigma^2 \operatorname{Tr}\left[(X^{\top} X)^{-2} X^{\top} X\right] = \sigma^2 \operatorname{Tr}[(X^{\top} X)^{-1}] \end{aligned} \end{split}\]

(2) \(\lambda \to \infty\)のとき、\(\mathrm{Bias} \to \|\beta\|^2, \mathrm{Variance} \to 0\)となる。

Biasについては、 \(\left(X^{\top} X+\lambda I\right)^{-1}=\lambda^{-1}\left(I+\frac{1}{\lambda} X^{\top} X\right)^{-1}\) となることとノルムの斉次性 \(\| c \boldsymbol{x} \|^2 = c^2 \| \boldsymbol{x} \|^2\)を用いると

\[ \begin{aligned} \lambda^2\left\|\left(X^{\top} X+\lambda I\right)^{-1} \beta\right\|^2 =\lambda^2\left\|\lambda^{-1}\left(I+\frac{1}{\lambda} X^{\top} X\right)^{-1} \beta\right\|^2 =\left\|\left(I+\frac{1}{\lambda} X^{\top} X\right)^{-1} \beta\right\|^2 \end{aligned} \]

となる。

\(\left(X^{\top} X+\lambda I\right)^{-1}=\lambda^{-1}\left(I+\frac{1}{\lambda} X^{\top} X\right)^{-1}\) について

\[ \lambda I + A = \lambda\left(I+\frac{1}{\lambda} A\right) \]

ゆえに

\[ (A+\lambda I)^{-1} =\left(\lambda\left(I+\frac{1}{\lambda} A\right)\right)^{-1} =\lambda^{-1}\left(I+\frac{1}{\lambda} A\right)^{-1} \]

となる。\(A=X^\top X\)とおけばよい

ここで\(\lambda \to \infty\)なので、 \(\left\|\left(I+\frac{1}{\lambda} X^{\top} X\right)^{-1} \beta\right\|^2\) は \(\frac{1}{\lambda} X^{\top} X \to 0\) となり、\(I^{-1} = I\)のため

\[\begin{split} \begin{aligned} \lim_{\lambda \to \infty} \mathrm{Bias} &= \left\|\left(I + \frac{1}{\lambda} X^{\top} X \right)^{-1} \beta\right\|^2 \\ &= \left\|I^{-1} \beta\right\|^2 \quad (\because \frac{1}{\lambda} X^{\top} X \to 0 ) \\ &= \left\| \beta \right\|^2 \quad (\because I^{-1} = I) \end{aligned} \end{split}\]

Variance \(\sigma^2 \operatorname{Tr}\left[\left(X^{\top} X+\lambda I\right)^{-2} X^{\top} X\right]\) については、

\[ (X^{\top} X +\lambda I)^{-2} =\lambda^{-2}\left(I+\frac{1}{\lambda} X^{\top} X \right)^{-2} \]

であり、 \(\left(I+\frac{1}{\lambda} X^{\top} X \right)^{-2}\) は \(\lambda \to \infty\) で\(I\)になるため

\[\begin{split} \begin{aligned} \lim_{\lambda \to \infty} \mathrm{Variance} &= \sigma^2 \operatorname{Tr}\left[\lambda ^{-2} X^{\top} X\right]\\ &= \lambda^{-2} \cdot \sigma^2 \operatorname{Tr}\left[X^{\top} X\right]\\ &= \frac{1}{\lambda^{2}} \cdot \sigma^2 \operatorname{Tr}\left[X^{\top} X\right]\\ &= 0 \end{aligned} \end{split}\]

バリアンス

\[\begin{split} \begin{aligned} \operatorname{Var}(\hat{\beta}^{\text{Ridge}}) &= \sigma^2 (X^\top X + \lambda I)^{-1} X^\top X(X^\top X + \lambda I)^{-1}\\ &\leq \sigma^2 (X^\top X)^{-1} = \operatorname{Var}(\hat{\beta}^{\text{LS}}) \end{aligned} \end{split}\]

https://hastie.su.domains/StatLearnSparsity_files/SLS.pdf

OLS推定量のバリアンス

線形回帰モデル\(y = X\beta + u\)をOLS推定量\(\hat{\beta}^{OLS} = (X^\top X)^{-1} X^\top y\)に代入して変形すると

\[\begin{split} \begin{align} \hat{\beta}^{OLS} &= (X^\top X)^{-1} X^\top y\\ &= (X^\top X)^{-1} X^\top (X\beta + u)\\ &= (X^\top X)^{-1} X^\top X\beta + (X^\top X)^{-1} X^\top u\\ &= \beta + (X^\top X)^{-1} X^\top u\\ \end{align} \end{split}\]

\(\hat{\beta}^{OLS} = \beta + (X^\top X)^{-1} X^\top u\) は \(u\sim N(0, \sigma^2 I)\)の仮定より、 \(\hat{\beta}^{OLS} \sim N(\beta, \sigma^2 (X^\top X)^{-1})\)となる。これは以下の定理を用いる。

定理

\(u\)を確率変数ベクトルとし、\(\mu\in\mathbb{R}^n\)、\(b\in\mathbb{R}^p\)、\(C\in\mathbb{R}^{p\times n}\)、\(\operatorname{rank} C = p\)とする。

\(u \sim N(\mu, \Sigma)\) のとき、\(Cu + b \sim N(C\mu + b, C\Sigma C^\top)\)

したがって \(\operatorname{Var}[\hat{\beta}^{OLS}] = \sigma^2 (X^\top X)^{-1}\) となる

証明

\(C := (X^\top X)^{-1} X^\top\)とおけば

\[\begin{split} \begin{aligned} C^\top = [(X^\top X)^{-1} X^\top]^\top &= X [(X^\top X)^{-1}]^\top \quad (\because (AB)^\top = B^\top A^\top) \\ &= X [(X^\top X)^\top]^{-1} \quad (\because (A^{-1})^\top = (A^\top)^{-1} ) \\ &= X (X^\top X)^{-1} \quad (\because (X^\top X)^\top=X^\top X ) \\ \end{aligned} \end{split}\]

であるため、\(u\sim N(0, \sigma^2 I)\) の仮定が満たされるとき、

\[\begin{split} \begin{aligned} &\beta + C u \sim N(\beta, \sigma^2 CC^\top)\\ &= \beta + (X^\top X)^{-1} X^\top u \sim N(\beta, \sigma^2 (X^\top X)^{-1} X^\top X (X^\top X)^{-1} )\\ &= \beta + (X^\top X)^{-1} X^\top u \sim N(\beta, \sigma^2 (X^\top X)^{-1}) \end{aligned} \end{split}\]

Ridge推定量のバリアンス

（導出部分は書籍を参考にしたのではなくChatGPTと会話しながら考えたので誤ってる可能性あり）

線形回帰モデル\(\boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{u}\)をRidge推定量に代入して変形すると

\[\begin{split} \begin{aligned} \hat{\boldsymbol{\beta}}^{\text{Ridge}} &= (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{y}\\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top (\boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{u})\\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{u}\\ \end{aligned} \end{split}\]

第1項は標本分布によるばらつきがないバイアス項なので第2項だけが関わる。つまり、

\[ \operatorname{Var}[\hat{\boldsymbol{\beta}}^{\text{Ridge}}] = \operatorname{Var}[ (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{u} ] \]

である。

分散には行列\(A\)と確率変数ベクトル\(u\)について、\(\operatorname{Var}[Au]= A\operatorname{Var}[u]A^\top\)となる定理があるので利用すると、\(\boldsymbol{A} :=(\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top\)と考えると

\[\begin{split} \begin{aligned} \operatorname{Var}[\hat{\boldsymbol{\beta}}^{\text{Ridge}}] &= \operatorname{Var}[ (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \boldsymbol{u} ]\\ &= \operatorname{Var}[ \boldsymbol{A} \boldsymbol{u} ]\\ &= \boldsymbol{A} \operatorname{Var}[\boldsymbol{u}] \boldsymbol{A}^\top \\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \operatorname{Var}[\boldsymbol{u}] [(\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top]^\top \\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \boldsymbol{X}^\top \operatorname{Var}[\boldsymbol{u}] \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X} + \lambda \boldsymbol{I})^{-1} \\ \end{aligned} \end{split}\]

モンテカルロ・シミュレーション

多項式をつかったデータに、多項式回帰をfittingしてみる

\[ y = \beta_0 + \beta_1 x + \beta_2 x^2 + u \]

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 0x7f9c881bbe20>

Show code cell source

Hide 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

Hide 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()

参考文献

リッジ回帰

• Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55-67.

• Ridge | L2 正則化を用いて多重共線性を回避してパラメーター推定を行う方法

• リッジ回帰 - 機械学習の「朱鷺の杜Wiki」

• 佐和隆光. (1979). 回帰分析.

ラグランジュの未定乗数法

• ラグランジュ緩和と双対問題

• 『言語処理のための機械学習入門』pp.14-19

• ラグランジュ双対問題について解説 - Qiita

• 『これならわかる最適化数学』