メインコンテンツへスキップ

正規方程式による最小二乗法の解の導出

··5 分·
Normal-Equation Least-Squares-Methods
Makoto Morinaga
著者
Makoto Morinaga
技術メモ、コーディング、環境構築のための個人ノート。
目次

線形回帰モデルのフィッティングでよく利用する最小二乗法(least squares method)について、正規方程式を用いて最小二乗法の解を導出します。

私の理解のため表現、理解等に誤りがあれば、コメント等でご指摘いただけますと嬉しいです。

線形回帰モデル
#

線形回帰モデル(liner regression)は、1つ以上の説明変数(\(\mathbf{x}\))と被説明変数(\(y\))との関係性をモデル化しています。 線形回帰モデルの基底関数(basis function)の違いによって以下のように呼ぶことができます。

  • 単回帰
  • 重回帰
  • 多項式回帰

単回帰
#

単回帰は、説明変数が1変数のみであり、説明変数(1次元)から \[ \mathbf{x}^\mathrm{T} = (1, x_1) \] と考え、 \[ \hat{y}(\mathbf{x}, \mathbf{w}) = w_0 + w_1x_1 \] で、説明変数と被説明変数の関係性をモデル化する。

ここで、線形回帰モデルの基底関数を \(\phi_0(\mathbf{x})=1, \phi_1(\mathbf{x})=x_1\) として、特徴ベクトル(feature vector)を \(\boldsymbol{\phi}\mathbf{(\mathbf{x})}=(\phi_0(\mathbf{x}), \phi_1(\mathbf{x}))^\mathrm{T}\) とおくと、

\begin{equation} \begin{split} \hat{y}(\mathbf{x}, \mathbf{w}) &= w_0 + w_1 x_1 \\ &=(w_0, w_1) \left(\begin{array}{c} 1 \\ x_1 \end{array} \right) \\ &=\mathbf{w}^\mathrm{T} \boldsymbol{\phi} \mathbf{(x)} \end{split} \end{equation}

となり、ベクトルの内積で表現できます。なお、\(\mathbf{w}\) はモデルの係数ベクトルです。

重回帰
#

重回帰は、説明変数が多変数となり、説明変数(D次元)から \[ \mathbf{x}^\mathrm{T} = (1, x_1, x_2, \ldots, x_D) \] と考え、 \[ \hat{y}(\mathbf{x}, \mathbf{w})=w_0 + w_1 x_1 + w_2 x_2 + \cdots+w_Dx_D \]

で、説明変数と被説明変数の関係性をモデル化します。

ここで、線形回帰モデルの基底関数を \(\phi_0(\mathbf{x})=1, \phi_1(\mathbf{x})=x_1, \ldots, \phi_D(\mathbf{x})=x_D\) として、特徴ベクトルを \(\boldsymbol{\phi}\mathbf{(x)}=(\phi_0(\mathbf{x}), \phi_1(\mathbf{x}), \ldots, \phi_D(\mathbf{x}))^\mathrm{T}\) とおくと、

\begin{equation} \begin{split} \hat{y}(\mathbf{x}, \mathbf{w}) &= w_0 + w_1 x_1 + w_2 x_2 + \cdots+w_Dx_D \\ &=(w_0, w_1, w_2, \ldots, w_D) \left( \begin{array}{c} 1 \\ x_1 \\ x_2 \\ \vdots \\ x_D \end{array} \right) \\ &=\mathbf{w}^\mathrm{T} \boldsymbol{\phi} \mathbf{(x)} \end{split} \end{equation}

となり、単回帰と同様に内積で表現できます。異なっている部分は基底関数の定義のみです。 つまり、重回帰は単回帰を多変数に拡張しています。

多項式回帰
#

多項式回帰は、説明変数に対して非線形な関数\(\phi\)(n次多項式や三角関数等)を基底関数として、 \[ \hat{y}(\mathbf{x}, \mathbf{w})=w_0 \phi_0(\mathbf{x}) + w_1 \phi_1(\mathbf{x}) + \cdots+w_H \phi_H(\mathbf{x}) \]

で、説明変数と被説明変数の関係性をモデル化します。

例として、多項式 \[ \hat{y}(\mathbf{x}, \mathbf{w})= w_0 + w_1 x_1 + +w_2 {x_1}^2 +w_3 \sin(x_1) \]

に回帰する場合を考えます。

この場合、基底関数は \(\phi_0(\mathbf{x})=1, \phi_1(\mathbf{x})=x_1, \phi_2(\mathbf{x})={x_1}^2, \phi_3(\mathbf{x})=\sin(x_1)\) 、特徴ベクトルは \(\boldsymbol{\phi}\mathbf{(x)}=(\phi_0(\mathbf{x}), \phi_1(\mathbf{x}), \phi_2(\mathbf{x}), \phi_3(\mathbf{x}))^\mathrm{T}\) となるため、以下のように表現できます。

\begin{equation} \begin{split} \hat{y}(\mathbf{x}, \mathbf{w}) &= w_0 + w_1 x_1 + +w_2 {x_1}^2 +w_3 \sin(x_1) \\ &=(w_0, w_1, w_2, w_3) \left( \begin{array}{c} 1 \\ x_1 \\ {x_1}^2 \\ \sin(x_1) \end{array} \right) \\ &=\mathbf{w}^\mathrm{T} \boldsymbol{\phi} \mathbf{(x)} \end{split} \end{equation}

多項式回帰では、多項式が \(\mathbf{x}\) に対して非線形になっていますが、 \(\mathbf{w}\) に対しては線形となっているため、線形回帰と呼ばれています。

上記から、基底関数を変更することで、単回帰、重回帰、線形回帰が同一のベクトル表現で表現できることが分かります。

最小二乗法
#

ここからが本題です。

線形回帰モデルのフィッティングでは、標準的に使われる手法として、最小二乗法があります。 最小二乗法は、実際の測定値 \(y\) が \(N\) 個ある場合に、測定値 \(y\) と線形回帰モデルでの予測値 \(\hat{y}\) の二乗誤差の総和

\begin{equation} \begin{split} E &= \sum_{n=1}^{N}(y_n - \hat{y}_n)^2 \\ &= \sum_{n=1}^{N}(y_n - \mathbf{w}^\mathrm{T} \boldsymbol{\phi} (\mathbf{x}_n))^2 \\ \end{split} \end{equation}

を最小にする \(\mathbf{w}\) を求めることで、線形回帰モデルをフィッティングさせます。

ここからは、この \(\mathbf{w}\) を導出していきます。

事前準備
#

二乗誤差の総和(\(E\))を最小にするには、 \(E\) をベクトル \(\mathbf{w}\) で微分して0とすれば、\(E\) を最小とする \(\mathbf{w}\) が導出できます。

\begin{equation} \begin{split} \frac{\partial{E}}{\partial{\mathbf{w}}} &= \frac{\partial}{\partial{\mathbf{w}}}\sum_{n=1}^{N}(y_n - \hat{y}_n)^2 \\ &= \frac{\partial}{\partial{\mathbf{w}}} \sum_{n=1}^{N}(y_n - \mathbf{w}^\mathrm{T} \boldsymbol{\phi} (\mathbf{x}_n))^2 = 0 \\ \end{split} \end{equation}

上記を解くために以下の公式等を準備します。

  • 計画行列
  • ベクトル要素の二乗和
  • 内積の微分
  • 二次形式の微分

計画行列
#

線形回帰モデルでは、単回帰、重回帰、多項式回帰のどれも \(\hat{y}(\mathbf{x}, \mathbf{w}) = \mathbf{w}^\mathrm{T} \boldsymbol{\phi} \mathbf{(x)}\) で表現できました。 そこで、実際の測定値 \(y\) が \(N\) 個ある場合、\(y\) の予測値 \(\hat{y}\) は、基底関数を用いて以下の行列で表現できます。

\begin{equation} \begin{split} \left( \begin{array}{c} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{array} \right) &= \left( \begin{array}{c} \mathbf{w}^\mathrm{T} \boldsymbol{\phi}(\mathbf{x}_1) \\ \mathbf{w}^\mathrm{T} \boldsymbol{\phi}(\mathbf{x}_2) \\ \vdots \\ \mathbf{w}^\mathrm{T} \boldsymbol{\phi}(\mathbf{x}_N) \end{array} \right) \\ &= \left( \begin{array}{c} \boldsymbol{\phi}(\mathbf{x}_1)^\mathrm{T} \\ \boldsymbol{\phi}(\mathbf{x}_2)^\mathrm{T} \\ \vdots \\ \boldsymbol{\phi}(\mathbf{x}_N)^\mathrm{T} \end{array} \right) \mathbf{w} \\ &= \begin{pmatrix} \phi_0(\mathbf{x}_1) & \phi_1(\mathbf{x}_1) & \cdots & \phi_H(\mathbf{x}_1)\\ \phi_0(\mathbf{x}_2) & \phi_1(\mathbf{x}_2) & \cdots & \phi_H(\mathbf{x}_2)\\ \vdots & & & \vdots \\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_H(\mathbf{x}_N)\\ \end{pmatrix} \mathbf{w} \end{split} \end{equation}

なお、\(\phi_0(\mathbf{x})\equiv 1\) と定義しています。

ここで、

\begin{equation} \hat{\mathbf{y}} =\left( \begin{array}{c} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{array} \right), \Phi = \begin{pmatrix} \phi_0(\mathbf{x}_1) & \phi_1(\mathbf{x}_1) & \cdots & \phi_H(\mathbf{x}_1)\\ \phi_0(\mathbf{x}_2) & \phi_1(\mathbf{x}_2) & \cdots & \phi_H(\mathbf{x}_2)\\ \vdots & & & \vdots \\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_H(\mathbf{x}_N)\\ \end{pmatrix} \end{equation}

とおくと、\(\hat{\mathbf{y}}= \Phi \mathbf{w}\) と表現でき、行列 \(\Phi\) を計画行列と呼びます。 この計画行列 \(\Phi\) は、線形回帰モデルの基底関数を行列で表現しています。

ベクトル要素の二乗和
#

ベクトル \(\mathbf{x} = (x_1, x_2, \ldots, x_n)^{\mathrm{T}}\) の要素の二乗和は

\begin{equation} \begin{split} x_1^2 + x_2^2 + \cdots x_n^2 &= (x_1, x_2, \ldots, x_n) \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \\ \end{array} \right) \\ &= \mathbf{x}^\mathrm{T} \mathbf{x} \end{split} \end{equation}

内積の微分
#

\(\mathbf{w}=(w_1,w_2,\ldots,w_n)^\mathrm{T}, \mathbf{x}=(x_1,x_2,\ldots,x_n)^\mathrm{T}\) とすると、

\[ \mathbf{w}^\mathrm{T} \mathbf{x} = w_1 x_1 + w_2 x_2 \cdots + w_n x_n \]

スカラーである内積 \(\mathbf{w}^\mathrm{T} \mathbf{x}\) をベクトル \(\mathbf{w}\) で微分することを考えます。ベクトルの各要素で微分すれば良いので、

\begin{equation} \begin{split} \frac{\partial}{\partial{\mathbf{w}}}\mathbf{w}^\mathrm{T} \mathbf{x} &= \left( \begin{array}{c} \frac{\partial}{\partial{w_1}} \mathbf{w}^\mathrm{T} \mathbf{x} \\ \frac{\partial}{\partial{w_2}} \mathbf{w}^\mathrm{T} \mathbf{x} \\ \vdots \\ \frac{\partial}{\partial{w_n}} \mathbf{w}^\mathrm{T} \mathbf{x} \\ \end{array} \right) \\ &= \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \\ \end{array} \right) \\ &= \mathbf{x} \end{split} \end{equation}

となります。

二次形式の微分
#

\(\mathbf{w}=(w_1,w_2,\ldots,w_n)^\mathrm{T}, \mathbf{X}\) を \(N\) 次正方行列として、スカラーである二次形式 \(\mathbf{w}^\mathrm{T}\mathbf{X}\mathbf{w}\) をベクトル \(\mathbf{w}\) で微分することを考えます。

まず、二次形式を変換します。

\begin{equation} \begin{split} \mathbf{w}^\mathrm{T}\mathbf{X}\mathbf{w} &= \mathbf{w}^\mathrm{T} \left( \begin{array}{c} \sum_{j=1}^N \mathbf{X}_{1j} w_j \\ \vdots \\ \sum_{j=1}^N \mathbf{X}_{Nj} w_j \\ \end{array} \right) \\ &= w_1 \sum_{j=1}^N \mathbf{X}_{1j} w_j + w_2 \sum_{j=1}^N \mathbf{X}_{2j} w_j + \cdots + w_N \sum_{j=1}^N \mathbf{X}_{Nj} w_j \\ &= \sum_{i=1}^N \sum_{j=1}^N \mathbf{X}_{ij} w_i w_j \\ \end{split} \end{equation}

ここで、\(\mathbf{w}\) の \(k\) 番目の成分 \(w_k\) で二次形式の微分を考えると、積の微分公式から

\begin{equation} \begin{split} \frac{\partial}{\partial{w_k}} \mathbf{w}^\mathrm{T}\mathbf{X}\mathbf{w} &= \frac{\partial}{\partial{w_k}} (\sum_{i=1}^N \sum_{j=1}^N \mathbf{X}_{ij} w_i w_j) \\ &= \sum_{i=1}^N \sum_{j=1}^N \mathbf{X}_{ij} (\frac{\partial}{\partial{w_k}} w_i) w_j + \sum_{i=1}^N \sum_{j=1}^N \mathbf{X}_{ij} w_i \frac{\partial}{\partial{w_k}} w_j \\ &= \sum_{j=1}^N \mathbf{X}_{kj} w_j + \sum_{i=1}^N \mathbf{X}_{ik} w_i \\ \end{split} \end{equation}

上記から二次形式の \(w_k\) での微分は、前半の項が \(\mathbf{Xw}\) の \(k\) 番目の成分、後半の項が \(\mathbf{X}^{\mathrm{T}}\mathbf{w}\) の \(k\) 番目の成分になっていることが分かります。

これを \(\mathbf{w}\) に拡張すると、二次形式の微分は以下のようになります。

\begin{equation} \begin{split} \frac{\partial}{\partial{\mathbf{w}}} \mathbf{w}^\mathrm{T}\mathbf{X}\mathbf{w} &= \left( \begin{array}{c} \sum_{j=1}^N \mathbf{X}_{1j} w_j + \sum_{i=1}^N \mathbf{X}_{i1} w_i \\ \sum_{j=1}^N \mathbf{X}_{2j} w_j + \sum_{i=1}^N \mathbf{X}_{i2} w_i \\ \vdots \\ \sum_{j=1}^N \mathbf{X}_{Nj} w_j + \sum_{i=1}^N \mathbf{X}_{iN} w_i \\ \end{array} \right) \\ &=\mathbf{Xw} + \mathbf{X}^{\mathrm{T}}\mathbf{w} \\ &=(\mathbf{X} + \mathbf{X}^{\mathrm{T}})\mathbf{w} \end{split} \end{equation}

二乗誤差の総和を最小にする値の導出
#

ここまでで準備が整ったので、\(\mathbf{w}\) を導出していきます。

まず、\(E\) を \(\mathbf{w}\) で微分がしやすいように、計画行列とベクトル要素の二乗和を使って展開します。

\begin{equation} \begin{split} E &= \sum_{n=1}^{N}(y_n - \hat{y}_n)^2 \\ &= \sum_{n=1}^{N}(y_n - \mathbf{w}^\mathrm{T} \boldsymbol{\phi} (\mathbf{x}_n))^2 \\ &= (\mathbf{y} - \Phi \mathbf{w})^{\mathrm{T}}(\mathbf{y} - \Phi \mathbf{w}) \\ &= \mathbf{y}^{\mathrm{T}} \mathbf{y} - \mathbf{y}^{\mathrm{T}} \Phi \mathbf{w} - (\Phi \mathbf{w})^{\mathrm{T}} \mathbf{y} + (\Phi \mathbf{w})^{\mathrm{T}} \Phi \mathbf{w} \\ \end{split} \end{equation}

ここで、 \(\mathbf{y}^{\mathrm{T}} \Phi \mathbf{w}\) はスカラーなので、 \(\mathbf{y}^{\mathrm{T}} \Phi \mathbf{w}=(\mathbf{y}^{\mathrm{T}} \Phi \mathbf{w})^{\mathrm{T}}\) が成り立ち、

\begin{equation} \begin{split} E &= \mathbf{y}^{\mathrm{T}} \mathbf{y} - (\mathbf{y}^{\mathrm{T}} \Phi \mathbf{w})^{\mathrm{T}} - (\Phi \mathbf{w})^{\mathrm{T}} \mathbf{y} + (\Phi \mathbf{w})^{\mathrm{T}} \Phi \mathbf{w} \\ &= \mathbf{y}^{\mathrm{T}} \mathbf{y} - (\Phi \mathbf{w})^{\mathrm{T}}\mathbf{y} - (\Phi \mathbf{w})^{\mathrm{T}} \mathbf{y} + (\mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}}) \Phi \mathbf{w} \\ &= \mathbf{y}^{\mathrm{T}} \mathbf{y} - 2(\Phi \mathbf{w})^{\mathrm{T}}\mathbf{y} + \mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}} \Phi \mathbf{w} \\ &= \mathbf{y}^{\mathrm{T}} \mathbf{y} - 2(\mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}})\mathbf{y} + \mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}} \Phi \mathbf{w} \\ &= \mathbf{y}^{\mathrm{T}} \mathbf{y} - 2\mathbf{w}^{\mathrm{T}} (\Phi^{\mathrm{T}}\mathbf{y}) + \mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}} \Phi \mathbf{w} \\ \end{split} \end{equation}

次に \(E\) を最小とする \(\mathbf{w}\) を導出するために、 \(\frac{\partial{E}}{\partial{\mathbf{w}}}=0\) とおくと、内積と二次形式の微分から

\begin{equation} \begin{split} \frac{\partial{E}}{\partial{\mathbf{w}}} &= \frac{\partial}{\partial{\mathbf{w}}}(\mathbf{y}^{\mathrm{T}} \mathbf{y} - 2\mathbf{w}^{\mathrm{T}} (\Phi^{\mathrm{T}}\mathbf{y}) + \mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}} \Phi \mathbf{w}) \\ &= -2 \Phi^{\mathrm{T}} \mathbf{y} + (\Phi^{\mathrm{T}}\Phi + (\Phi^{\mathrm{T}}\Phi)^{\mathrm{T}}) \mathbf{w} \\ &= -2 \Phi^{\mathrm{T}} \mathbf{y} + (\Phi^{\mathrm{T}}\Phi + \Phi^{\mathrm{T}}\Phi) \mathbf{w} \\ &= -2 \Phi^{\mathrm{T}} \mathbf{y} + 2\Phi^{\mathrm{T}}\Phi \mathbf{w} = 0\\ \end{split} \end{equation}

となる。よって、方程式

\[ \Phi^{\mathrm{T}} \Phi \mathbf{w} = \Phi^{\mathrm{T}} \mathbf{y} \]

が得られます。この導出された方程式を正規方程式(normal equation)といいます。

したがって、正規方程式から \(\Phi^{\mathrm{T}} \Phi\) の逆行列が存在するならば、

\[ \mathbf{w} = (\Phi^{\mathrm{T}} \Phi)^{-1} \Phi^{\mathrm{T}} \mathbf{y} \]

が、二乗誤差の総和を最小にする \(\mathbf{w}\) となる。これで \(\mathbf{w}\) の導出が完了です。

関連記事

M1 MacへのNodejs 14.xのインストール
··1 分
Macos Nodejs
Gatsby v2からGatsby v3への移行
··2 分
Blog Gatsby
Gatsbyによるブログ構築
··4 分
Netlify Gatsby Blog