最小二乗法による線形回帰モデル

リッジ回帰を説明する上で最小二乗法による線形回帰モデルの推定方法の理解が必要になりますので、まずこちらから説明します。

説明変数が2以上の線形回帰モデル（重回帰モデル）は以下で表されます。

$\displaystyle{ \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{u} }$

ここで $\boldsymbol{y}$ は目的変数、 $\boldsymbol{X}$ は説明変数、 $\boldsymbol{\beta}$ は切片項を含む回帰係数（偏回帰係数）、 $\boldsymbol{u}$ は誤差項です。これの回帰係数 $\boldsymbol{\beta}$ を推定する方法として最小二乗法があります。これは観測値 $\boldsymbol{y}$ とモデルによる予測値 $\boldsymbol{\hat{y}}$ の差の二乗が $\boldsymbol{\beta}$ について最小になるような $\beta$ を見つける方法です。観測値と予測値の差の二乗を数式で書くと、

$\displaystyle{ h(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left \{ y_i - (\beta_0 + \beta_1x_{i1} + \ldots + \beta_p x_{ip}) \right \}^2 }$

この式をベクトルと行列を用いて書き直すと、

$\displaystyle{ \begin{aligned} h(\boldsymbol{\beta}) &= \left \{\boldsymbol{y} - (\beta_0 + \beta_1 \boldsymbol{x}_1 + \ldots + \beta_p \boldsymbol{x}_p) \right \}^2 \\ &= \left \{\boldsymbol{y} - (\beta_0 + \beta_1 \boldsymbol{x}_1 + \ldots + \beta_p \boldsymbol{x}_p) \right \}^\mathrm{T} \left \{\boldsymbol{y} - (\beta_0 + \beta_1 \boldsymbol{x}_1 + \ldots + \beta_p \boldsymbol{x}_p) \right \} \\ &= (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta})^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}) \\ &= \|\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \|^2 \end{aligned} }$

となります。なお、１行目から２行目の式変形はベクトルの2乗は自身の内積になるので $\boldsymbol{x}^2 = \boldsymbol{x} \cdot \boldsymbol{x} = \boldsymbol{x}^\mathrm{T} \boldsymbol{x}$ という関係を用いています。この観測値と予測値の差の二乗を残差平方和(RSS)と呼びます。ではこのRSSが $\boldsymbol{\beta}$ について最小となるようにするので、 $\boldsymbol{\beta}$ で偏微分して $\boldsymbol{0}$ と置きます。ベクトルで微分しているので、0ベクトルとなることに注意してください。

$\displaystyle{ \frac{\partial h(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2\boldsymbol{X}^\mathrm{T} (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}) = \boldsymbol{0} }$

式変形をして $\boldsymbol{\beta}$ について解きます。

$\displaystyle{ \begin{aligned} \boldsymbol{X}^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) &= \boldsymbol{0} \\ \boldsymbol{X}^\mathrm{T} \boldsymbol{y} - \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} &= \boldsymbol{0} \\ \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} &= \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \end{aligned} }$

$\beta$ イコールの形にするためには $\boldsymbol{X}^\mathrm{T} \boldsymbol{X}$ の逆行列を両辺の左からかける必要があります。逆行列はすべての行列において計算できるわけではなく、行列がフルランクである必要があります。 $\boldsymbol{X}$ がフルランクであるとは、言い換えると説明変数同士で完全相関（相関係数が1）している組が無いということです。なお、フルランクでない行列の場合、特異と呼びます。よって、ここではフルランクであることを仮定します。

$\displaystyle{ \begin{aligned} (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \therefore \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \end{aligned} }$

これが最小二乗法を用いた場合の回帰係数（偏回帰係数）の推定量になります。

最小二乗法とバイアスについて

上の最小二乗法で求めた回帰係数は不偏推定量といって、バイアス（偏り）がありません。これは回帰係数の期待値（平均）が母数（真の値）に一致することを意味しています。しかしながら、回帰係数の選び方は無数にあるため、バイアスはあるもののもっと誤差を減らせる（精度を上げられる）ような回帰係数の選び方もあるかもしれません。そこで、バイアスを増やすことによってこれを実現することを考えます。変数選択もその一つであり、今回取り上げるリッジ回帰やラッソ回帰などの縮小推定があります。縮小推定では推定する回帰係数について制限を設けて縮小または0にする手法です。いずれにしても、縮小推定を用いて推定した回帰係数はバイアスを持つ可能性があるため、解釈には注意が必要になると考えられます。

# 最小二乗法による線形回帰をリッジ回帰に拡張する

最小二乗法による線形回帰は下式で表される残差平方和(RSS)を最小とするような $\boldsymbol{\beta}$ を求めることでした。

$\displaystyle{ h(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left \{ y_i - (\beta_0 + \beta_1x_{i1} + \ldots + \beta_p x_{ip}) \right \}^2 }$

リッジ回帰を含む縮小推定では回帰係数 $\boldsymbol{\beta}$ に制限を制限を設けると上で書きました。縮小推定という名前が示しているように回帰係数 $\boldsymbol{\beta}$ が小さくなるようにします。残差平方和を最小とするというアプローチは変わりませんが、回帰係数 $\boldsymbol{\beta}$ を小さくしたいので回帰係数 $\boldsymbol{\beta}$ の大きさに応じてペナルティを課すようにします。具体的には、上の残差平方和の式に回帰係数 $\boldsymbol{\beta}$ の大きさを足し合わせます。単純に回帰係数 $\boldsymbol{\beta}$ を足し合わせないのは回帰係数 $\boldsymbol{\beta}$ はマイナスになることもあり得るからです。大きさという表現は数学的には様々な表現方法がありますが、回帰係数 $\boldsymbol{\beta}$ はベクトルであるためノルムになり、絶対値を取る方法(L1ノルム)、2乗する方法（L2ノルム）などがあります。前者を用いるとラッソ回帰に、後者だとリッジ回帰になります。

回帰係数 $\boldsymbol{\beta}$ の大きさを表す式を罰則項と呼ぶことにします。 $\beta_0$ は切片項に該当しますが、これは $\boldsymbol{y}$ の原点であり今回の関心から外れるものであり、罰則項には含めません。また、説明変数を中心化することで切片の推定量は単純に $\boldsymbol{\bar{y}}$ になります。これは切片の推定量 $\hat{\alpha}$ が

$\displaystyle{ \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x} }$

から求めることができるので、説明変数が中心化されていることは説明変数の平均値が0になり、第2項が0になることからわかります。

以降では説明変数が中心化されていることを前提とします。つまり説明変数 $\boldsymbol{X}$ の次元が $n \times p$ （ $n$ 個のデータと $p$ 個の説明変数）です。罰則項を以下のように表現します。

$\displaystyle{ \lambda \sum_{j=1}^{p} \beta^2_j }$

ここで $\lambda$ は罰則項の強さを調整するハイパーパラメータです。ではこの項を残差平方和に足し合わせます。

$\displaystyle{ h(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left \{ y_i - (\beta_1x_{i1} + \ldots + \beta_p x_{ip}) \right \}^2 + \lambda \sum_{j=1}^{p} \beta^2_j }$

最小二乗法のときと同様にベクトルと行列を用いて書き換えます。

$\displaystyle{ \begin{aligned} h(\boldsymbol{\beta}) &= (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta})^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}) +\lambda\boldsymbol{\beta}^\mathrm{T}\boldsymbol{\beta} \\ &= \|\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \|^2 + \lambda \|\boldsymbol{\beta}\|^2 \end{aligned} }$

これを最小化するので偏微分して $\boldsymbol{0}$ とおきます。

$\displaystyle{ \frac{\partial h(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2 \boldsymbol{X}^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) + 2\lambda \boldsymbol{\beta} = \boldsymbol{0} }$

$\boldsymbol{\beta}$ について解きます。

$\displaystyle{ \begin{aligned} -2 \boldsymbol{X}^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) + 2\lambda \boldsymbol{\beta} &= \boldsymbol{0} \\ -\boldsymbol{X}^\mathrm{T} \boldsymbol{y} + \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} + \lambda \boldsymbol{\beta} &= \boldsymbol{0} \\ \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} + \lambda \boldsymbol{\beta} &= \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I}) \boldsymbol{\beta} &= \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \end{aligned} }$

$\boldsymbol{I}$ は $p \times p$ の単位行列です。最小二乗法の時と同様に逆行列を用いて $\boldsymbol{\beta}$ イコールの形にしますが、今回のリッジ回帰では $\lambda\boldsymbol{I}$ つまり $\lambda$ が対角成分に足されています。これによって、必ずフルランクになりますので逆行列を持ちます。（個人的な感想ですが、これを始めて見たときにとても感動しました。）よって、

$\displaystyle{ \begin{aligned} (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}(\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I}) \boldsymbol{\beta} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \therefore \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \end{aligned} }$

これがリッジ回帰における回帰係数の推定量です。最小二乗法のものと比べてみます。

$\displaystyle{ \begin{aligned} 最小二乗法 : \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \cr リッジ回帰 : \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \end{aligned} }$

その違いは $\lambda\boldsymbol{I}$ が $\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}$ に足されていることだけです。ハイパーパラメータ $\lambda$ が $0$ であれば最小二乗法の式と等価に、 $\lambda \to \infty$ で回帰係数が $0$ に縮小されます。