2020.12.29 00:00

正則化（Ridge回帰とLasso回帰）について

はじめに

今回は、回帰分析において過学習を防いだり、多重共線性に対応したりするために使われるRidge回帰や、変数の数が標本数より多いような時に変数選択の方法として使われるLasso回帰について、理論を整理しようと思います。

これら2つの手法ともには正則化という手法で説明されるものです。まずは正則化の観点からRidge回帰とLasso回帰を理解し、さらにこれらの手法が確率モデルでも説明可能であることを示そうと思います。

確率モデルで説明可能な手法であるということは、データセットとモデルから尤度を計算しベイズ推定可能であるということなので、ゆくゆくはRstanでこれら2つの手法を実装してみたいと思っています。

本記事の構成は以下の通りです。

参考にした本や記事は以下のとおりです。

線形回帰の導入

以下の線形回帰を考えます。

$$ \hat{y} = w_0 + w_1\phi_1(\boldsymbol{x}) + \ldots, w_H\phi_H(\boldsymbol{x}) \tag{1} $$

ここで$\phi(\boldsymbol{x})$は$\boldsymbol{x}$の関数です。

これを行列で表現すると、

$$ \left( \begin{array}{ccc} \hat{y}_1 \\\\ \hat{y}_2 \\\\ \vdots \\\\ \hat{y}_N \end{array} \right) = \left( \begin{array}{ccc} \phi_0(\boldsymbol{x}_1) & \phi_1(\boldsymbol{x}_1) & \cdots & \phi_H(\boldsymbol{x}_1) \\\\ \phi_0(\boldsymbol{x}_2) & \phi_1(\boldsymbol{x}_2) & \cdots & \phi_H(\boldsymbol{x}_2) \\\\ \cdots & \cdots & \ddots & \cdots \\\\ \phi_0(\boldsymbol{x}_N) & \phi_1(\boldsymbol{x}_N) & \cdots & \phi_H(\boldsymbol{x}_N) \end{array} \right) \left( \begin{array}{ccc} w_0 \\\\ w_1 \\\\ \vdots \\\\ w_H \end{array} \right) \tag{2} $$

ただし、$\phi_0({\boldsymbol{x}_n}) = 1$です。

これを行列を使って以下のように表記します。

$$ \hat{\boldsymbol{y}} = \boldsymbol{\Phi}\boldsymbol{w} \tag{3} $$

一般に$\boldsymbol{\Phi}$を計画行列と呼びます。また$\phi_0(\boldsymbol{x}),\ldots,\phi_H(\boldsymbol{x})$を基底関数、$\boldsymbol{\phi}(\boldsymbol{x}) = \left(\phi_0(\boldsymbol{x}),\phi_1(\boldsymbol{x}),\ldots, \phi_H(\boldsymbol{x}) \right)^T$を$\boldsymbol{x}$の特徴ベクトルと呼びます。

memo

基底関数を、ベクトル$\boldsymbol{x}$をスカラーに変換する解析的な関数であると想定していますが、 通常の重回帰モデルも$\boldsymbol{\phi}(\boldsymbol{x}) = \boldsymbol{x}$とすれば$(3)$式で説明できます！その場合基底関数は各変数それぞれを出力する関数になります。

このとき、以下の公式が得られます。

◆線形回帰モデルの解 線形回帰モデル$\hat{\boldsymbol{y}} = \boldsymbol{\Phi}\boldsymbol{w}$における$\boldsymbol{w}$の最小二乗解は $$ \boldsymbol{w} = \left(\boldsymbol{\Phi}^T\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^T\boldsymbol{y} \tag{4} $$


◆証明 データ全体の誤差(損失関数)は、 $$ E = \sum_{n=1}^{N} \left( y_n - \boldsymbol{\phi}(\boldsymbol{x}_n)\boldsymbol{w} \right)^2 = (\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{w})^T(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}) = \boldsymbol{y}^T\boldsymbol{y} - 2\boldsymbol{w}^T(\boldsymbol{\Phi}^T\boldsymbol{y}) + \boldsymbol{w}^T\boldsymbol{\Phi}\boldsymbol{\Phi}\boldsymbol{w} \tag{5} $$ これを$\boldsymbol{w}$で微分すると、 $$ \cfrac{\partial E}{\partial\boldsymbol{w}} = -2\boldsymbol{\Phi}^T\boldsymbol{y} + 2\boldsymbol{\Phi}^T\boldsymbol{\Phi}\boldsymbol{w} \tag{6} $$ よって、$E$を最小化する$\boldsymbol{w}$は、$\cfrac{\partial E}{\partial\boldsymbol{w}}=0$より下記の通り得られる。 $$ \boldsymbol{w} = \left(\boldsymbol{\Phi}^T\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^T\boldsymbol{y} $$

正則化

線形回帰モデルの係数パラメータは$(4)$式で計算できることが分かりましたが、右辺には$\boldsymbol{\Phi}^T\boldsymbol{\Phi}$の逆行列が存在するため、$\boldsymbol{\Phi}^T\boldsymbol{\Phi}$が正則でない場合、係数パラメータの解を得ることが出来ません。

memo

$n$次正方行列$\boldsymbol{A}$について、 $$ \boldsymbol{A}\boldsymbol{B} = \boldsymbol{B}\boldsymbol{A} = \boldsymbol{I} $$ となる$n$次正方行列$\boldsymbol{B}$が存在するとき、$\boldsymbol{A}$は正則であるという。

では、$\boldsymbol{\Phi}^T\boldsymbol{\Phi}$が正則でないのは具体的にどんな場合かとうと、$\boldsymbol{\phi}_ {i}(\boldsymbol{x}) = \boldsymbol{\phi}_ {j}(\boldsymbol{x})$、もしくは$\boldsymbol{\phi}_ {i}(\boldsymbol{x}) = \alpha\boldsymbol{\phi}_ {j}(\boldsymbol{x})~~(\alpha \in \mathbb{R})$となるような$\boldsymbol{\phi}(\boldsymbol{x})$の組が存在する場合が典型的です。このような場合、$\boldsymbol{\phi}_ {i}(\boldsymbol{x})$と$\boldsymbol{\phi}_ {j}(\boldsymbol{x})$の相関は1または-1になっています。また$\boldsymbol{\phi}(\boldsymbol{x}_ {i}) \simeq \alpha\boldsymbol{\phi}(\boldsymbol{x}_ {j})$だったりする場合では、$\boldsymbol{\Phi}^T\boldsymbol{\Phi}^{-1}$の要素が大きくなり、結果的に$\boldsymbol{w}$の絶対値が大きくなってしまい、過学習が生じるようです。このような場合では、$\boldsymbol{\phi}_ {i}(\boldsymbol{x})$と$\boldsymbol{\phi}_ {j}(\boldsymbol{x})$はかなり強い相関関係になっています1。

memo

$\boldsymbol{\phi}_ {i}(\boldsymbol{x}) = \alpha\boldsymbol{\phi}_ {j}(\boldsymbol{x})~~(\alpha \in \mathbb{R})$となるような$\boldsymbol{\phi}(\boldsymbol{x})$の組が存在する場合、 $$ \mathrm{rank}\left(\boldsymbol{\Phi}^T\boldsymbol{\Phi}\right) < H $$ と行列$\boldsymbol{\Phi}^T\boldsymbol{\Phi}$はランク落ちしており、ランク落ちした行列は逆行列が存在しない。

係数パラメータの絶対値が大きくなるのを避けるための工夫が正則化です。正則化では、$(5)$式で得られたデータ全体の誤差に、$\boldsymbol{w}$の絶対値が大きくなることによるペナルティを課すように設定します。 このペナルティの設定方法が、Ridge回帰とLasso回帰の唯一の違いです。

◆Ridge回帰とLasso回帰の損失関数 Ridge回帰における損失関数は、 $$ E = (\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{w})^T(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}) + \alpha||\boldsymbol{w}||_2^2 \tag{7} $$ Lasso回帰における損失関数は、 $$ E = (\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{w})^T(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}) + \alpha||\boldsymbol{w}||_1^1 \tag{8} $$ ただし、$\alpha \geq 0$


memo

$||\boldsymbol{x}||_p$はベクトルの$p$ノルムといい $$ ||\boldsymbol{x}||_p = \sqrt[p]{|x_1|^p+\ldots,|x_n|^p} $$ を意味する。定義より$||\boldsymbol{x}||_1^1$は$\boldsymbol{x}$のマンハッタン距離を、$||\boldsymbol{x}||_2^2$は$\boldsymbol{x}$のユークリッド距離を意味する。

損失関数を上記のように設定することが$\boldsymbol{\Phi}^T\boldsymbol{\Phi}$が正則でない場合に有効である仕組みを把握するために、Ridge回帰の場合についてみていきます。

$(7)$式を$\boldsymbol{w}$で微分すると、

$$ \cfrac{\partial E}{\partial\boldsymbol{w}} = -2\boldsymbol{\Phi}^T\boldsymbol{y} + 2\boldsymbol{\Phi}^T\boldsymbol{\Phi}\boldsymbol{w} + 2\alpha\boldsymbol{w} \tag{9} $$

よって、$E$を最小化する$\boldsymbol{w}$は、$\cfrac{\partial E}{\partial\boldsymbol{w}}=0$より下記のとおり得られます。

$$ \boldsymbol{w} = \left(\boldsymbol{\Phi}^T\boldsymbol{\Phi} + \alpha\boldsymbol{I}\right)^{-1}\boldsymbol{\Phi}^T\boldsymbol{y} \tag{10} $$

リッジ回帰は、もともとは$(7)$式のように損失関数に$\boldsymbol{w}$の絶対値を小さくする為の項を設定したものでしたが、結果的に$\boldsymbol{\Phi}^T\boldsymbol{\Phi}$の対角要素に微小量$\alpha$を足すことで、逆行列計算の対象が正則行列であることを確実にし、計算を安定化させていることが分かります。

memo

任意の行列$\boldsymbol{A} ∈ \mathbb{R}^{m×n}$、ベクトル$\boldsymbol{x} ∈ \mathbb{R}^m \neq \boldsymbol{0}$について、 $$ \boldsymbol{x}^T\left(\boldsymbol{A}\boldsymbol{A}^T\right)\boldsymbol{x} = \left(\boldsymbol{A}^T\boldsymbol{x}\right)^T\left(\boldsymbol{A}^T\boldsymbol{x}\right) = ||\boldsymbol{A}^T\boldsymbol{x}||_2^2 \geq 0 $$ だから、$\boldsymbol{A}\boldsymbol{A}^T$は半正定値行列。 また $$ \boldsymbol{x}^T\left(\alpha\boldsymbol{I}\right)\boldsymbol{x} = \alpha||\boldsymbol{x}||_2^2 > 0 $$ だから、$\alpha\boldsymbol{I}$は正定値行列。このとき $$ \boldsymbol{x}^T\left(\boldsymbol{A}\boldsymbol{A}^T + \alpha\boldsymbol{I}\right)\boldsymbol{x} = ||\boldsymbol{A}^T\boldsymbol{x}||_2^2 + \alpha||\boldsymbol{x}||_2^2 > 0 $$ よって$\boldsymbol{A}^T\boldsymbol{A} + \alpha\boldsymbol{I}$は正定値行列であり、正則行列。

Lasso回帰については同様の方法では正則化との関係が見えてきませんでした。ただ係数パラメータの多くが0になる仕組みについてはこちらの記事に詳しかったです。

確率モデルでRidge回帰を捉える

前節で説明したRidge回帰は確率分布を用いた以下のモデルでも説明できます。

◆Ridge回帰の確率モデル $(7)$のRidge回帰と等価の確率モデルは、 $$ \begin{cases} p(y | \boldsymbol{w}, \boldsymbol{x}) = \mathrm{Normal}(y | \boldsymbol{w}^T\boldsymbol{\phi}(\boldsymbol{x}), \sigma^2) \\
p(\boldsymbol{w}) = \mathrm{Normal}(\boldsymbol{w}|0, \rho^{-1}\boldsymbol{I}_H) \end{cases} \tag{11} $$ ただし、$\alpha = \rho\sigma^2$


◆証明 $y$と$\boldsymbol{w}$の同時確率密度について以下の通り仮定する。 $$ p(y,\boldsymbol{w}|\boldsymbol{x}) = p(y | \boldsymbol{w},\boldsymbol{x})p(\boldsymbol{w}) $$ これを$\boldsymbol{w}$の関数$F(\boldsymbol{w})$とみたとき、独立同分布のデータセット$\boldsymbol{y}$、$\boldsymbol{X}=(\boldsymbol{x}_1,\ldots,\boldsymbol{x}_N)^T$がすべて観測された時の$F(\boldsymbol{w})$を最大化する$\boldsymbol{w}$が最尤推定解である。 $$ \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \begin{split} \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}} \log F(\boldsymbol{w}) &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}} \log \prod_{n=1}^{N}p(y_n|\boldsymbol{w}, \boldsymbol{x}_ n)p(\boldsymbol{w})\\\\ &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[ \sum_{n=1}^{N}\log\mathrm{Normal}(y_n| \boldsymbol{w}^T\boldsymbol{\phi}(\boldsymbol{x}),\sigma^2) + \log\mathrm{Normal}(\boldsymbol{w}|\boldsymbol{0}, \rho^{-1}\boldsymbol{I}_H)\right] \\\\ &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[\log\mathrm{MultiNormal}(\boldsymbol{y}| \boldsymbol{\Phi}\boldsymbol{w},\sigma^2\boldsymbol{I}_N) + \log\mathrm{Normal}(\boldsymbol{w}|\boldsymbol{0}, \rho^{-1}\boldsymbol{I}_H)\right] \\\\ &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[ -\cfrac{1}{2\sigma^2}\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right)^T\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right) -\cfrac{\rho}{2}\boldsymbol{w}^T\boldsymbol{w} + C \right] \\\\ &=\argmin_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right)^T\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right) + \rho\sigma^2||\boldsymbol{w}||_2^2\right] \\\\ &= \argmin_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right)^T\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right) + \alpha||\boldsymbol{w}||_2^2\right] \end{split} $$ これは$(7)$式のRidge回帰における損失関数$E$の$\boldsymbol{w}$についての最小化である。

確率モデルでLasso回帰を捉える

Lasso回帰は確率分布を用いた以下のモデルでも説明できます。

◆Lasso回帰の確率モデル $(8)$のLasso回帰と等価の確率モデルは、 $$ \begin{cases} p(y | \boldsymbol{w}, \boldsymbol{x}) = \mathrm{Normal}(y | \boldsymbol{w}^T\boldsymbol{\phi}(\boldsymbol{x}), \sigma^2) \\
p(\boldsymbol{w}) = \prod_{h=1}^{H}\mathrm{Laplace}(w_h | 0, b) \end{cases} \tag{12} $$ ただし、$\alpha = \cfrac{2\sigma^2}{b}$


◆証明 $y$と$\boldsymbol{w}$の同時確率密度について以下の通り仮定する。 $$ p(y,\boldsymbol{w}|\boldsymbol{x}) = p(y | \boldsymbol{w},\boldsymbol{x})p(\boldsymbol{w}) $$ これを$\boldsymbol{w}$の関数$F(\boldsymbol{w})$とみたとき、独立同分布のデータセット$\boldsymbol{y}$、$\boldsymbol{X}=(\boldsymbol{x}_1,\ldots,\boldsymbol{x}_N)^T$がすべて観測された時の$F(\boldsymbol{w})$を最大化する$\boldsymbol{w}$が最尤推定解である。 $$ \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \begin{split} \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}} \log F(\boldsymbol{w}) &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}} \log \prod_{n=1}^{N}p(y_n|\boldsymbol{w}, \boldsymbol{x}_ n)p(\boldsymbol{w})\\\\ &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[ \sum_{n=1}^{N}\log\mathrm{Normal}(y_n| \boldsymbol{w}^T\boldsymbol{\phi}(\boldsymbol{x}),\sigma^2) + \sum_{h=1}^{H}\log\mathrm{Laplace}(w_h|0, b)\right] \\\\ &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[\log\mathrm{MultiNormal}(\boldsymbol{y}| \boldsymbol{\Phi}\boldsymbol{w},\sigma^2\boldsymbol{I}_ N) + \sum_{h=1}^{H}\log\mathrm{Laplace}(w_h|0, b)\right] \\\\ &= \argmax_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[ -\cfrac{1}{2\sigma^2}\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right)^T\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right) - \cfrac{1}{b}\sum_{h=1}^{H}|w_h| + C \right] \\\\ &=\argmin_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right)^T\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right) + \cfrac{2\sigma^2}{b}||\boldsymbol{w}||_1^1\right] \\\\ &= \argmin_{\boldsymbol{w}∈\mathbb{R}^{H+1}}\left[\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right)^T\left(\boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w}\right) + \alpha||\boldsymbol{w}||_1^1\right] \end{split} $$ これは$(8)$式のLasso回帰における損失関数$E$の$\boldsymbol{w}$についての最小化である。

memo

確率密度$p(x)$が $$ p(x) = \cfrac{1}{2b}\exp\left[-\cfrac{|x - \mu|}{b}\right] $$ である確率分布をラプラス分布といい、確率変数$X$がラプラス分布に従うことを $$ X \sim \mathrm{Laplace}(\mu, b) $$ と表す。

蛇足

今回の記事はHugoのshortcode機能を駆使して文章に枠をもたせてみました。今までの記事よりかなり見やすくなっていませんか？？しばらくはこのスタイルで行こうと思います。

最後に今回の記事のイメージソングをあげておきます。ブログ名の由来になった曲です。

  1. 回帰分析で生じるこのような問題は多重共線性と呼ばれています。 ↩︎

回帰

コメントを書く


※ コメントは承認されると表示されます

承認されたコメント一覧