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

2021年7月21日（ 2022年6月24日）

この記事では、線形回帰モデルのフィッティングでよく利用する最小二乗法(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{aligned} \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{aligned}

となり、ベクトルの内積で表現できます。なお、 $\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{aligned} \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{aligned}

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

多項式回帰

多項式回帰は、説明変数に対して非線形な関数 $\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{aligned} \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{aligned}

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

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

最小二乗法

ここからが本題です。

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

\begin{aligned} 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{aligned}

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

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

事前準備

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

\begin{aligned} \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{aligned}

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

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

計画行列

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

\begin{aligned} \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{aligned}

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

ここで、

\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}

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

ベクトル要素の二乗和

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

\begin{aligned} 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{aligned}

内積の微分

$\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{aligned} \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{aligned}

となります。

二次形式の微分

$\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{aligned} \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{aligned}

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

\begin{aligned} \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{aligned}

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

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

\begin{aligned} \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{aligned}

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

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

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

\begin{aligned} 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{aligned}

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

\begin{aligned} 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{aligned}

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

\begin{aligned} \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{aligned}

となります。よって、方程式

\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}$ の導出が完了です。

終わりに

今回は、正規方程式を用いて最小二乗法の解を導出しました。ただ、これを手計算で行うことは現実的に大変なので、次はPythonで正規方程式を用いた最小二乗法を実装してみます。

数学

最小二乗法

線形代数

正規方程式

Posted by mako

1	Pythonの_(アンダースコア)の扱いまとめ 1524件のビュー
2	M1 MacへのNodejs 14.xのインストール 236件のビュー
3	EmacsでのlspをベースとしたPython開発環境 169件のビュー
4	正規方程式による最小二乗法の解の導出 155件のビュー
5	EmacsでのElpyをベースとしたPython開発環境 127件のビュー

1	Native compilationを有効化したEmacs28のインストール 2022年4月27日（ 2022年6月24日）
2	textlintとFlycheckによるEmacs上でのテキスト校正 2021年12月13日（ 2022年6月24日）
3	EmacsでのlspをベースとしたPython開発環境 2021年12月2日（ 2022年6月24日）
4	Docker Compose上のNeo4jにcsvデータをインポートする方法 2021年7月31日（ 2022年6月24日）
5	正規方程式による最小二乗法の解の導出 2021年7月21日（ 2022年6月24日）