最小二乗法とその周辺

はじめに

　今回は、最小二乗法を入り口にして、最尤推定、MAP推定、リッジ回帰について説明したい。

最小二乗法

　入力データ \{{\bf x}_1,\cdots,{\bf x}_N\},\;{\bf x}_n \in \mathbb{R}^Mと出力データ \{y_1,\cdots,y_N\},\;y\in\mathbb{R}が与えられているとする。ここで、{\bf x}_nはM次元列ベクトルである。これらを用いて次の量を定義する。

    \begin{equation*} X= \begin{pmatrix} {\bf x}_1^T \\ \vdots \\ {\bf x}_N^T \end{pmatrix}\in\mathbb{R}^{N\times M} \end{equation*}

    \begin{equation*} {\bf y}= \begin{pmatrix} y_1 \\ \vdots \\ y_N \end{pmatrix}\in\mathbb{R}^{N} \end{equation*}

Tは転置を表す。XはN\times M行列、{\bf y}はN次元列ベクトルである。後の計算を簡単にするため、次式が成り立つとする。

    \begin{align*} \sum_{n=1}^N{\bf x}_n&={\bf 0}\\ \sum_{n=1}^Ny_n&=0 \end{align*}

すなわち、データは前処理により「中心化」されている。このとき次の線形回帰モデルを考える。

(1)   \begin{align*} {\bf y}=X{\boldsymbol \omega}+{\boldsymbol \epsilon} \end{align*}

ここで、{\boldsymbol \omega}\in \mathbb{R}^M、{\boldsymbol \epsilon}\in \mathbb{R}^Nである。{\boldsymbol \epsilon}は誤差ベクトル（ノイズ）を表す。式(1)は

    \begin{align*} y_n={\boldsymbol \omega}^T{\bf x}_n+\epsilon_n,\;n=1,\cdots,N \end{align*}

をまとめて書いたものである。{\bf x}_n,y_nは中心化されているので、式(1)に定数項（バイアス項）は存在しないことに注意する。さて、{\boldsymbol \epsilon}={\bf y}-X{\boldsymbol \omega}が最小になるように{\boldsymbol \omega}を決めることを考える。その際、次の二乗誤差を用いる。

(2)   \begin{align*} L(\boldsymbol \omega)&=\left({\bf y}-X{\boldsymbol \omega}\right)^T\left({\bf y}-X{\boldsymbol \omega}\right)\\ &={\bf y}^T{\bf y}-2{\bf y}^TX{\boldsymbol \omega}+{\boldsymbol \omega}^TX^TX{\boldsymbol \omega} \end{align*}

上式を{\boldsymbol \omega}で偏微分し、その値を0と置くと

(3)   \begin{align*} \dfrac{\partial L({\boldsymbol \omega})}{\partial {\boldsymbol \omega}}=2X^TX{\boldsymbol \omega}-2X^T{\bf y}=0 \end{align*}

を得る。これより

    \begin{align*} X^TX{\boldsymbol \omega}=X^T{\bf y} \end{align*}

を得る。X^TXの逆行列が存在すれば

    \begin{align*} {\boldsymbol \omega}=\left(X^TX\right)^{-1}X^T{\bf y} \end{align*}

となる。ここまでの手順が最小二乗法である。すなわち、二乗誤差を最小にする{\boldsymbol\omega}を決定する手法である。

最尤推定

　最尤推定（Maximum Likelihood Estimation）の入り口は確率である。最初に{\bf y},X,{\boldsymbol \omega}が同時に実現する確率p\left({\bf y},X,{\boldsymbol \omega}\right)を考える。このとき、ベイズの定理より

    \begin{align*} p\left({\boldsymbol \omega}|{\bf y},X\right)=\dfrac{p({\bf y}|X,{\boldsymbol \omega})p({\boldsymbol \omega})}{p({\bf y}|X)} \end{align*}

が成り立つ。左辺は、観測値{\bf y},Xが与えられたとき{\boldsymbol \omega}が実現する確率を表す。この確率を最大にするような{\boldsymbol \omega}を求めることが目的である。
　さて、右辺の分母にある確率は{\boldsymbol \omega}に依存しないので定数として扱うことができる。

(4)   \begin{align*} p\left({\boldsymbol \omega}|{\bf y},X\right)\propto p({\bf y}|X,{\boldsymbol \omega})p({\boldsymbol \omega}) \end{align*}

さらに、p({\boldsymbol \omega})を定数とみなす近似を行う。

(5)   \begin{align*} p\left({\boldsymbol \omega}|{\bf y},X\right)\propto p({\bf y}|X,{\boldsymbol \omega}) \end{align*}

右辺は、X,{\boldsymbol\omega}が実現しているとき{\bf y}が実現する確率であるが、この確率を{\boldsymbol\omega}の関数と見るとき尤度関数と呼ぶ。尤度関数を最大にするように{\boldsymbol\omegaを決める手法を最尤推定と言う。
　各観測値は独立に同一の分布から得られると仮定すると

(6)   \begin{align*} p({\bf y}|X,{\boldsymbol \omega})=\prod_{n=1}^N p(y_n|{\bf x}_n,{\boldsymbol \omega}) \end{align*}

と書くことができる。各データの尤度を正規分布で近似する。

(7)   \begin{align*} p(y_n|{\bf x}_n,{\boldsymbol \omega})=\dfrac{1}{\sqrt{2\pi\sigma^2}} \exp{\left(-\dfrac{\left(y_n-{\boldsymbol \omega}^T{\bf x}_n\right)^2}{2\sigma^2}\right)} \end{align*}

式(7)を式(6)に代入し、その値を最大にするような{\boldsymbol\omega}を求めることが目的であるが、そのまま扱うより対数をとり最大値を探した方が楽である。そこで次の対数尤度を定義する。

    \begin{align*} L({\boldsymbol\omega})\equiv\ln{p({\bf y}|X,{\boldsymbol \omega})}=-\dfrac{1}{2\sigma^2}\sum_{n=1}^N\left(y_n-{\boldsymbol\omega}^T{\bf x}_n\right)^2-\dfrac{N}{2}\ln{2\pi\sigma^2} \end{align*}

右辺を最大化することは

    \begin{align*} \sum_{n=1}^N\left(y_n-{\boldsymbol\omega}^T{\bf x}_n\right)^2 \end{align*}

を最小化することと等価である。これは先に考えた最小二乗法の式(2)である。すなわち、最小二乗法とは、尤度関数として正規分布を採用した最尤推定と同じであることが分かる。

MAP推定とリッジ回帰

　ここまでの議論で最小二乗法と最尤推定が等価であることを示した。そして求めたい{\boldsymbol\omegaの値は次式で与えられた。

    \begin{align*} {\boldsymbol\omega}=\left(X^TX\right)^{-1}X^T{\bf y} \end{align*}

しかし、X^TXが逆行列を持たない場合、{\boldsymbol\omega}を求めることができない。X^TXが逆行列を持たない代表的な例として、Xを構成する2つのデータ{\bf x}_n,{\bf x}_mが線形従属となる場合を挙げることができる。このような場合を多重共線性があると言い、このときXの行列式とX^TXの行列式は0になる。
　多重共線性を回避する1つの方法としてMAP推定（Maximum A Posteriori Estimation）がある。MAP推定では、式(4)で定数とみなしたp({\boldsymbol\omega})に対しても正規分布を与える。

    \begin{align*} p({\boldsymbol \omega})=\dfrac{1}{\sqrt{2\pi\mu^2}} \exp{\left(-\dfrac{{\boldsymbol \omega}^T{\boldsymbol \omega}}{2\mu^2}\right)} \end{align*}

そして、先と同じように対数尤度を最大化することを考える。

    \begin{align*} L({\boldsymbol\omega})\equiv\ln{p({\bf y}|X,{\boldsymbol \omega})}=-\dfrac{1}{2\sigma^2}\sum_{n=1}^N\left(y_n-{\boldsymbol\omega}^T{\bf x}_n\right)^2&-\dfrac{N}{2}\ln{2\pi\sigma^2}\\ &-\dfrac{1}{2\mu^2}{\boldsymbol\omega}^T{\boldsymbol\omega}-\dfrac{1}{2}\ln{2\pi\mu^2} \end{align*}

{\boldsymbol\omega}に依存しない項を無視して整理すると

    \begin{align*} L({\boldsymbol\omega})\equiv\ln{p({\bf y}|X,{\boldsymbol \omega})}&=-\dfrac{1}{2\sigma^2}\sum_{n=1}^N\left(y_n-{\boldsymbol\omega}^T{\bf x}_n\right)^2-\dfrac{1}{2\mu^2}{\boldsymbol\omega}^T{\boldsymbol\omega}+\mbox{Const.}\\ &=-\dfrac{1}{2\sigma^2}\left\{ \sum_{n=1}^N\left(y_n-{\boldsymbol\omega}^T{\bf x}_n\right)^2+\dfrac{\sigma^2}{\mu^2}{\boldsymbol\omega}^T{\boldsymbol\omega}\right\}+\mbox{Const.} \end{align*}

を得る。改めて\eta=\sigma^2/\mu^2と置くと

(8)   \begin{align*} \sum_{n=1}^N\left(y_n-{\boldsymbol\omega}^T{\bf x}_n\right)^2+\eta\;{\boldsymbol\omega}^T{\boldsymbol\omega} \end{align*}

を最小化すれば良いことになる。この解も解析的に求めることができて

    \begin{align*} {\boldsymbol\omega}=\left(X^TX+\eta I\right)^{-1}X^T{\bf y} \end{align*}

を得る。ここで、Iは単位行列である。行列X^TX+\eta Iに対しては、Xに多重共線性の問題があっても常に逆行列を計算することができる。すなわち、MAP推定は多重共線性の問題を回避できる手法である。ところで、式(8)は最小二乗法の式に正則化項（\eta{\boldsymbol\omega^T\boldsymbol\omega}）を追加した式とみなすこともできる。この立場からの定式化をリッジ回帰と呼ぶ。つまり、この例においてMAP推定とリッジ回帰は等価である。
　パラメータ\etaは、あらかじめ決めておくパラメータ（ハイパーパラメータ）である。このパラメータを決めるにはまた別の議論が必要であるが、今回は割愛する。

まとめ

　今回は、最小二乗法から始めて最尤推定、MAP推定、リッジ回帰について説明した。お互いに関連した手法であること、入り口は異なるが同じ結果を得る手法であることを示した。
　今回の説明では、X^TXが逆行列を持たない例として多重共線性だけを考えたが、実はもうひとつ重要な例がある。N<Mとなる場合である。すなわち、サンプル数Nの方が{\bf x}_nの次元数Mより少なくなる場合である。このとき、多重共線性の有無とは無関係にX^TXの逆行列を計算できない。この場合の対処法の説明は機会を改めたい。

Kumada Seiya

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

最近の記事

  • 関連記事
  • おすすめ記事
  • 特集記事

アーカイブ

カテゴリー

PAGE TOP