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

2021年7月21日 （ 2023年12月11日）

この記事では、線形回帰モデルのフィッティングでよく利用する最小二乗法(least squares method)について、正規方程式を用いて最小二乗法の解を導出します。 私の理解のため表現、理解等に誤りがあれば、コメント等でご指摘いただけますと嬉しいです。

目次

線形回帰モデル

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

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

単回帰

単回帰は、説明変数が1変数のみであり、説明変数(1次元)から

xT=(1,x1)\mathbf{x}^\mathrm{T} = (1, x_1)

と考え、

y^(x,w)=w0+w1x1\hat{y}(\mathbf{x}, \mathbf{w}) = w_0 + w_1x_1

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

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

y^(x,w)=w0+w1x1=(w0,w1)(1x1)=wTϕ(x)\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}

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

重回帰

重回帰は、説明変数が多変数となり、説明変数(D次元)から

xT=(1,x1,x2,…,xD)\mathbf{x}^\mathrm{T} = (1, x_1, x_2, \ldots, x_D)

と考え、

y^(x,w)=w0+w1x1+w2x2+⋯+wDxD\hat{y}(\mathbf{x}, \mathbf{w})=w_0 + w_1 x_1 + w_2 x_2 + \cdots+w_Dx_D

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

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

y^(x,w)=w0+w1x1+w2x2+⋯+wDxD=(w0,w1,w2,…,wD)(1x1x2⋮xD)=wTϕ(x)\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次多項式や三角関数等)を基底関数として、

y^(x,w)=w0ϕ0(x)+w1ϕ1(x)+⋯+wHϕH(x)\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})

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

例として、多項式

y^(x,w)=w0+w1x1++w2x12+w3sin⁡(x1)\hat{y}(\mathbf{x}, \mathbf{w})= w_0 + w_1 x_1 + +w_2 {x_1}^2 +w_3 \sin(x_1)

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

この場合、基底関数はϕ0(x)=1,ϕ1(x)=x1,ϕ2(x)=x12,ϕ3(x)=sin⁡(x1)\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)、特徴ベクトルはϕ(x)=(ϕ0(x),ϕ1(x),ϕ2(x),ϕ3(x))T\boldsymbol{\phi}\mathbf{(x)}=(\phi_0(\mathbf{x}), \phi_1(\mathbf{x}), \phi_2(\mathbf{x}), \phi_3(\mathbf{x}))^\mathrm{T}となるため、以下のように表現できます。

y^(x,w)=w0+w1x1++w2x12+w3sin⁡(x1)=(w0,w1,w2,w3)(1x1x12sin⁡(x1))=wTϕ(x)\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}

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

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

最小二乗法

ここからが本題です。

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

E=∑n=1N(yn−y^n)2=∑n=1N(yn−wTϕ(xn))2\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}

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

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

事前準備

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

∂E∂w=∂∂w∑n=1N(yn−y^n)2=∂∂w∑n=1N(yn−wTϕ(xn))2=0\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}

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

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

計画行列

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

(y^1y^2⋮y^N)=(wTϕ(x1)wTϕ(x2)⋮wTϕ(xN))=(ϕ(x1)Tϕ(x2)T⋮ϕ(xN)T)w=(ϕ0(x1)ϕ1(x1)⋯ϕH(x1)ϕ0(x2)ϕ1(x2)⋯ϕH(x2)⋮⋮ϕ0(xN)ϕ1(xN)⋯ϕH(xN))w\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}

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

ここで、

y^=(y^1y^2⋮y^N),Φ=(ϕ0(x1)ϕ1(x1)⋯ϕH(x1)ϕ0(x2)ϕ1(x2)⋯ϕH(x2)⋮⋮ϕ0(xN)ϕ1(xN)⋯ϕH(xN))\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}

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

ベクトル要素の二乗和

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

x12+x22+⋯xn2=(x1,x2,…,xn)(x1x2⋮xn)=xTx\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}

内積の微分

w=(w1,w2,…,wn)T,x=(x1,x2,…,xn)T\mathbf{w}=(w_1,w_2,\ldots,w_n)^\mathrm{T}, \mathbf{x}=(x_1,x_2,\ldots,x_n)^\mathrm{T}とすると、

wTx=w1x1+w2x2⋯+wnxn\mathbf{w}^\mathrm{T} \mathbf{x} = w_1 x_1 + w_2 x_2 \cdots + w_n x_n

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

∂∂wwTx=(∂∂w1wTx∂∂w2wTx⋮∂∂wnwTx)=(x1x2⋮xn)=x\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}

となります。

二次形式の微分

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

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

wTXw=wT(∑j=1NX1jwj⋮∑j=1NXNjwj)=w1∑j=1NX1jwj+w2∑j=1NX2jwj+⋯+wN∑j=1NXNjwj=∑i=1N∑j=1NXijwiwj\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}

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

∂∂wkwTXw=∂∂wk(∑i=1N∑j=1NXijwiwj)=∑i=1N∑j=1NXij(∂∂wkwi)wj+∑i=1N∑j=1NXijwi∂∂wkwj=∑j=1NXkjwj+∑i=1NXikwi\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}

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

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

∂∂wwTXw=(∑j=1NX1jwj+∑i=1NXi1wi∑j=1NX2jwj+∑i=1NXi2wi⋮∑j=1NXNjwj+∑i=1NXiNwi)=Xw+XTw=(X+XT)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}

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

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

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

E=∑n=1N(yn−y^n)2=∑n=1N(yn−wTϕ(xn))2=(y−Φw)T(y−Φw)=yTy−yTΦw−(Φw)Ty+(Φw)TΦ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}

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

E=yTy−(yTΦw)T−(Φw)Ty+(Φw)TΦw=yTy−(Φw)Ty−(Φw)Ty+(wTΦT)Φw=yTy−2(Φw)Ty+wTΦTΦw=yTy−2(wTΦT)y+wTΦTΦw=yTy−2wT(ΦTy)+wTΦTΦw\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}

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

∂E∂w=∂∂w(yTy−2wT(ΦTy)+wTΦTΦw)=−2ΦTy+(ΦTΦ+(ΦTΦ)T)w=−2ΦTy+(ΦTΦ+ΦTΦ)w=−2ΦTy+2ΦTΦ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}

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

ΦTΦw=ΦTy\Phi^{\mathrm{T}} \Phi \mathbf{w} = \Phi^{\mathrm{T}} \mathbf{y}

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

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

w=(ΦTΦ)−1ΦTy\mathbf{w} = (\Phi^{\mathrm{T}} \Phi)^{-1} \Phi^{\mathrm{T}} \mathbf{y}

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

終わりに

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

Posted by mako
関連記事
関連記事はありません
コメント
...

Copyright © 2018-2023 makoのノート All Rights Reserved.