最小二乗法による近似直線の係数を行列計算で求めてみた。証明もしてみた

スポンサーリンク

最小二乗法による近似直線の係数を行列計算で求める方法

最小二乗法を使って近似直線を引くには、行列計算を使うと考え方が簡単です。

まず準備運動として、平面上の2点\left(x_1, y_1\right), \left(x_2, y_2\right)を通る直線y=a x+bを引いてみましょう。このとき、

    $$ Y=\left[\begin{array}{c} y_1 \\ y_2 \end{array}\right], X=\left[\begin{array}{cc} x_1 & 1\\ x_2 & 1 \end{array}\right], A=\left[\begin{array}{c} a \\ b \end{array}\right] $$

とおけば、

    $$ \left[\begin{array}{c} y_1 \\ y_2 \end{array}\right] = \left[\begin{array}{cc} x_1 & 1\\ x_2 & 1 \end{array}\right] \left[\begin{array}{c} a \\ b \end{array}\right] $$

と書けます。よって、

    $$ Y=X A $$

と表すことができます。求めたいAベクトルは、左からXの逆行列をかけて、

    $$ X^{-1} Y=A $$

とすれば求まります。ここまでは問題ありませんよね?

最小二乗法を使う場合は、点の数が3点以上になる場合であって、

    $$ Y=\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right], X=\left[\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{array}\right], A=\left[\begin{array}{c} a \\ b \end{array}\right], \, \left(n \geq 3\right) $$

として、

    $$ Y \simeq XA $$

と書くことができます(直線が必ずしもすべての点を通らないのでニアリイコールを使っています)。でも、Xが正方行列ではないので、逆行列を求めることができません。
ところが、左から転置行列をかけてしまえば、正方行列になるではありませんか。
つまり、

    $$ X^\mathsf{T} Y=X^\mathsf{T} XA \eqno (1) $$

とすれば、X^\mathsf{T}2 \times nの行列で、Xn \times 2の行列なので、X^\mathsf{T} X2 \times 2の正方行列になります。
この正方行列が正則であれば、逆行列が存在します。それを左から掛けると、

    $$ \left( X^\mathsf{T} X \right)^{-1} X^\mathsf{T} Y=A \eqno (2) $$

となり、何かしらのAが解析的に求まるではありませんか!
実は、これが最小二乗法で求めた解と一致するのです。

何と素晴らしい!

ちなみに、ここでは、X-Yの2次元平面について述べていますが、3次元以上の多次元でも同様に考えることができます。

念のため、実データで検算してみる

行列を使わずに係数を求めた前回の記事で使ったデータ\left(-2,-6\right), \left(1,3\right), \left(4,0\right)で検算してみます。この答えは\left(a, b\right)=\left(1,-2\right)でした。

まず、データから配列を作ります。

    $$ Y=\left[\begin{array}{r} -6 \\ 3 \\ 0 \end{array}\right], X=\left[\begin{array}{rr} -2 & 1\\ 1 & 1 \\ 4 & 1 \end{array}\right] $$

よって、

    \begin{eqnarray*} \left[\begin{array}{r} a \\ b \end{array}\right] &=&\left( X^\mathsf{T} X \right)^{-1} X^\mathsf{T} Y \\ &=&\left( \left[\begin{array}{rrr} -2 & 1 & 4\\ 1 & 1 & 1 \end{array}\right] \left[\begin{array}{rrr} -2 & 1 \\ 1 & 1 \\ 4 & 1 \end{array}\right] \right)^{-1} \left[\begin{array}{rrr} -2 & 1 & 4\\ 1 & 1 & 1 \end{array}\right] \left[\begin{array}{r} -6 \\ 3 \\ 0 \end{array}\right] \\ &=& \left[\begin{array}{rr} 21 & 3\\ 3 & 3 \end{array}\right]^{-1} \left[\begin{array}{r} 15 \\ -3 \end{array}\right]\\ &=&\frac{1}{54} \left[\begin{array}{rrr} 3 & -3 \\ -3 & 21 \end{array}\right] \left[\begin{array}{r} 15 \\ -3 \end{array}\right] \\ &=&\frac{1}{54} \left[\begin{array}{r} 54 \\ -108\end{array}\right] \\ &=&\left[\begin{array}{r} 1 \\ -2\end{array}\right] \end{eqnarray*}

となり、あらかじめ示しておいた答えと一致しました。

証明

念のため、証明します。

まず、誤差の二乗和Eを考えます。

    \begin{eqnarray*} E &=& \sum_{i=1}^{n}\left\{y_i -\left(a x_i + b\right)\right\}^{2} \\ &=& \left(\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right]- \left[\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{array}\right] \left[\begin{array}{c} a \\ b \end{array}\right] \right)^\mathsf{T} \left(\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right]- \left[\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{array}\right] \left[\begin{array}{c} a \\ b \end{array}\right] \right)\end{eqnarray*}

ここで、

    $$ Y=\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right], X=\left[\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{array}\right], A=\left[\begin{array}{c} a \\ b \end{array}\right] $$

とおくと、

    \begin{eqnarray*} E &=& \left(Y - X A\right)^\mathsf{T}\left(Y - X A\right) \\ &=& \left(Y^\mathsf{T} - A^\mathsf{T} X^\mathsf{T}\right)\left(Y - X A\right) \\ &=& Y^\mathsf{T} Y-Y^\mathsf{T} X A - A^\mathsf{T} X^\mathsf{T} Y + A^\mathsf{T} X^\mathsf{T} X A \\ &=& Y^\mathsf{T} Y - 2 A^\mathsf{T} X^\mathsf{T} Y + A^\mathsf{T} X^\mathsf{T} X A \end{eqnarray*}

となります。ただし、最後の式変形では、

    $$ Y^\mathsf{T} X A = \left(\left( Y^\mathsf{T} X A \right)^\mathsf{T} \right)^\mathsf{T} = \left( A^\mathsf{T} X^\mathsf{T} Y \right)^\mathsf{T} =A^\mathsf{T} X^\mathsf{T} Y $$

であることを使いました。A^\mathsf{T} X^\mathsf{T} Yは、1 \times 22 \times nn \times 1の行列の積なのでスカラ値であり、転置をとっても値が変わりません。

ここで、

    $$ A^\mathsf{T} X^\mathsf{T} Y=\begin{bmatrix}a & b \end{bmatrix}\begin{bmatrix}x_1 & x_2 & \cdots & x_n \\1 & 1 & \cdots & 1 \end{bmatrix}\begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}=\sum_{i=1}^n\left(ax_i+b\right)y_i $$

および

    $$ A^\mathsf{T} X^\mathsf{T} XA=\begin{bmatrix}a & b \end{bmatrix}\begin{bmatrix}x_1 & x_2 & \cdots & x_n \\1 & 1 & \cdots & 1 \end{bmatrix}\begin{bmatrix}x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix}\begin{bmatrix}a \\b \end{bmatrix} =\sum_{i=1}^n\left(ax_i+b\right)^2 $$

であることに注意して、Ea,\ bそれぞれで偏微分して、勾配ベクトルを求めると、

    \begin{align*} \begin{bmatrix}\frac{\partial E}{\partial a} \\ \frac{\partial E}{\partial b} \end{bmatrix} &=\begin{bmatrix}\frac{\partial\left(Y^\mathsf{T} Y - 2 A^\mathsf{T} X^\mathsf{T} Y + A^\mathsf{T} X^\mathsf{T} X A\right)}{\partial a} \\\frac{\partial \left(Y^\mathsf{T} Y - 2 A^\mathsf{T} X^\mathsf{T} Y + A^\mathsf{T} X^\mathsf{T} X A\right)}{\partial b} \end{bmatrix} \\ &=\begin{bmatrix}\frac{\partial \left(0-2\sum_{i=1}^n\left(a x_i+b\right)y_i+\sum_{i=1}^n\left(a x_i+b\right)^2\right)}{\partial a} \\\frac{\partial \left(0-2\sum_{i=1}^n\left(a x_i+b\right)y_i+\sum_{i=1}^n\left(a x_i+b\right)^2\right)}{\partial b} \end{bmatrix} \\ &=\begin{bmatrix}-2\sum_{i=1}^n x_i y_i \\-2\sum_{i=1}^n y_i \end{bmatrix}+\begin{bmatrix}\sum_{i=1}^n 2x_i\left(a x_i+b\right) \\\sum_{i=1}^n 2\left(a x_i+b\right) \end{bmatrix}\\ &=-2\begin{bmatrix}x_1 & x_2 & \cdots & x_n \\1 & 1 & \cdots & 1 \end{bmatrix}\begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}+2\begin{bmatrix}x_1 & x_2 & \cdots & x_n \\1 & 1 & \cdots & 1 \end{bmatrix}\begin{bmatrix}x_1 & 1 \\x_2 & 1 \\ \vdots & \vdots \\ x_n & 1\end{bmatrix}\begin{bmatrix}a \\b \end{bmatrix} \\ &=-2 X^\mathsf{T} Y+2 X^\mathsf{T} X A \end{align*}

となります。

誤差の二乗和が最小値になる条件は、この勾配ベクトルが0になるときで、

    $$ X^\mathsf{T} Y = X^\mathsf{T} X A$$

です。これは、(1)に示した条件になるので、(2)で係数が計算できることが示されました。

コメント

タイトルとURLをコピーしました