3次元空間の回帰直線の求め方

2次元平面上に存在する複数の点からその回帰直線を最小二乗法で求める問題はよくあるが、3次元空間に存在する複数の点からその回帰直線を最小二乗法で求める問題はあまり見ないようなので、ここにその解法を書いておきます。

原点を\(O\)とする3次元座標系において、点\({ P }_{ 0 } ({ x }_{ 0 },{ y }_{ 0 },{ z }_{ 0 })\)から点\(P (x,y,z)\)に至るベクトル\(\vec {P_0P}\)(これが求める直線)とすると、ベクトル\(\vec { OP }\)は次式で表される。
\[
\vec { OP } ={ \vec { O{ P }_{ 0 } }  }+\vec {P_0P} \qquad (1)\\
={ \vec { O{ P }_{ 0 } }  }+t\vec {e} \qquad (2)
\]
ここで、\(\vec {P_0P}=t\vec {e}\)であり、\(t\)はベクトル\(\vec {P_0P}\)の長さ、\(\vec {e}\)は\(\vec {P_0P}\)の単位ベクトルである。
\(\vec {e}\)の成分を\((a, b, c)\)とすると、点Pの各座標は
\[
\begin{array}{c}
x = x_0 + at \qquad (3) \\
y = y_0 + bt \qquad (4) \\
z = z_0 + ct \qquad (5)
\end{array}
\]
となる。ただし、\(a^2+b^2+c^2=1\)。
方向単位ベクトル\(\vec {e}\)のz軸方向の成分が\(0\)でない、つまり\(c≠0\)ならば式(5)から
\[t = (z_0-z) / c \qquad (6)\]となる。
もし、z軸方向の成分が\(0\)、つまり\(c=0\)ならばベクトル\(\vec { t }\)が示す直線は\(x,y\)平面上にあるので、通常の2次元平面上の直線となり、普通の最小二乗法で解を求めることができる。
さて、式(6)を(3)、(4)に代入すると式(7)、(8)が得られる。
\[
x = Az + B \qquad (7)\\
y = Cz + D \qquad (8)\\
y = Ex + F \qquad (9)
\]
式(7)、(8)から\(z\)を消去して\(x,y\)の関係式(9)を算出した。ここでは、\(A=a/c、B=(x_0-Az_0)、C=b/c、D=y_0-Bz_0\)と置き換えている。
式(7)、(8)、(9)はそれぞれ\(zx\)平面、\(zy\)平面、\(xy\)平面における直線式を表している。
この各々の2次元平面において、以前書いた擬似逆行列で回帰直線を求める方法を用いて式(7)、(8)、(9)の係数6個を求める。
つまり、3次元空間に存在する全ての点\((x,y,z)\)の中から、\((x,z)\)の座標を式(7)に入れて係数\(A,B\)を求め、\((y,z)\)の座標を式(8)に入れて係数\(C,D\)を求め、\((y,x)\)の座標を式(9)に入れて係数\(E,F\)を求める。

これにより最初に設定した\(P_0\)の未知数3個と、方向ベクトル\(\vec {P_0P}\)の3つの未知数を全て算出できる。なぜならば、\(\vec {P_0P}\)の未知数は、\(t\)の長さとその単位ベクトルの未知数2個であるからだ(\(a^2+b^2+c^2=1\)より未知数は2個)。このようにして擬似逆行列で求めた直線は各点からの距離が近い回帰直線となる。

<補足>
式(3)、(4)、(5)では、\(t\)はベクトル\(\vec {P_0P}\)の長さとし、\(\vec {e}\)を\(\vec {P_0P}\)の単位ベクトルとした。そして、単位ベクトルの条件(\(a^2+b^2+c^2=1\)から、係数を一つ消去することで、未知数が6個であると説明した。
この考え方でよいのだが、より簡単に計算するには、式(5)の\(ct\)→t’と置き換えることで式(3)、(4)、(5)は以下のように書き換えられる。
\[
\begin{array}{c}
x = x_0 + a’t’ \qquad (3′) \\
y = y_0 + b’t’ \qquad (4′) \\
z = z_0 + t’ \qquad (5′)
\end{array}
\]
この式から、未知数は\(x_0, y_0, z_0, a’, b’, t’\)の6個であることが分かる。この後は上式(6)以下と同様に変形して、未知数を算出することができ、3次元空間の回帰直線が求まる。

擬似逆行列と回帰方程式(最小二乗法)

前の投稿では擬似逆行列を用いて、回帰直線を得る方法を解説しました。

今回は、擬似逆行列を用いて、高次の回帰曲線を得る方法を説明します。

一般に、n次元曲線の方程式は次のように表現されます。

\[y={ a }_{ 0 }+{ a }_{ 1 }{ x }+{ a }_{ 1 }{ x }^{2}+\cdots +{ a }_{ n }{ x }^{ n } \qquad (1)\]

この\(x,y\)が、点\((x_1,y_1)\)、\((x_2,y_2)\)、\((x_3,y_3)\)、、、やその近傍を通過して欲しいならば、式(1)に代入して、その係数を求めます。その際、点の数は係数\(a_i\) \((i=0\)~\(n)\)の数\(n+1\)以上が必要です。今その点数を\(m\)個とすると、以下のように書き下せます。

\[{ y }_{ 1 }={ a }_{ 0 }+{ a }_{ 1 }{ x }_{ 1 }+{ a }_{ 2 }{ x }_{ 1 }^{ 2 }\cdots +{ a }_{ n }{ x }_{ 1 }^{ n }\qquad (2)\]

\[{ y }_{2}={ a }_{ 0 }+{ a }_{ 1 }{ x }_{2}+{ a }_{ 2 }{ x }_{2}^{ 2 }\cdots +{ a }_{ n }{ x }_{2}^{ n }\qquad (3)\]

\[\cdot \cdot \cdot \cdot \cdot \cdot \cdot \]

\[{ y }_{m}={ a }_{ 0 }+{ a }_{ 1 }{ x }_{m}+{ a }_{ 2 }{ x }_{m}^{ 2 }\cdots +{ a }_{ n }{ x }_{m}^{ n }\qquad (4)\]

これを行列で表すと

\[\left[ \begin{matrix} { y }_{ 1 } \\ { y }_{ 2 } \\\cdot\\ { y }_{ m } \end{matrix} \right] =\left[ \begin{matrix} 1 & { x }_{1} & { x }_{1}^2 & \cdots & { x }_{1}^{n} \\ 1 & { x }_{ 2 } & { x }_{2}^2 & \cdots & { x }_{2}^{n} \\ \cdot & \cdot & \cdot & \cdots & \cdot \\ 1 & { x }_{ m } & { x }_{m}^2 & \cdot & { x }_{ m }^{n} \end{matrix} \right] \left[ \begin{matrix} { a }_{ 0 } \\ { a }_{ 1 } \\ { a }_{2}\\\cdot\\ { a }_{n} \end{matrix} \right] \qquad (5)\]
となります。これは前の投稿の式(5)と同じ式に簡略化できます。

\[Y = XA\qquad (6)\]

従って、前の投稿と同じ方法で係数\(A\)を求めることができるので、前の投稿の式(15)と同じ次の結果が得られます。

\[A = (X^TX)^{-1}X^TY\qquad (7)\]

このように、擬似逆行列を用いれば、高次の回帰曲線であろうが、単なる回帰直線であろうが、同じ計算手法で簡易に算出することができます。なんと有り難い計算手法でしょう!これは最小二乗法などを用いた場合の煩雑な計算に較べ、大きな特長といえます。

擬似逆行列と回帰直線(最小二乗法)

数学の美しさ、完璧さに幾度となく驚いた私でしたが、擬似逆行列もその一つでした。その例を一つご紹介します。
直線の方程式は次のようにあらわされます。
\[y= ax + b\qquad (1)\]
この直線が点(\(x_1,y_1\))と点(\(x_2,y_2\))を通過するならば、
\[y_1= ax_1 + b\qquad (2)\]

\[y_2= ax_2 + b\qquad (3)\]が成立します。これを行列で書くと

\[\left[ \begin{matrix} y_1 \\ y_2 \end{matrix} \right]=\left[ \begin{matrix} { x }_{ 1 } & 1 \\ { x }_{ 2 } & 1 \end{matrix} \right]\left[ \begin{matrix} a \\ b \end{matrix} \right] \qquad (4)\]

となります。各行列を簡略化して表すと

\[Y= XA\qquad (5)\]となります。

\(X\)が逆行列\({ X }^{ -1 }\)を持つならば(正則行列)、式(5)の両辺に左側から\({ X }^{ -1 }\)を乗じて

\[{ X }^{ -1 }Y= { X }^{ -1 }XA\qquad (6)\]

ここで、

\[{ X }^{ -1 }X= E\qquad (7)\]

\(E\)は単位行列。よって、式(6)は下記となります。

\[{ X }^{ -1 }Y = A \qquad (8)\]

結果として得られた式(8)は、式(5)の\(X, Y\)が単なる変数である場合に得られる

\[A=\frac { Y }{ X }\qquad (9)\]

と同様の形をしています。

まあ、ここまでは普通の逆行列を用いた計算で、そういう風に表現できるんだ、で済む話ですが、ここから擬似逆行列が本領を発揮します。

式(1)は2点を通る直線という条件において、式(8)により係数\(a,b\)を出すことができます。しかし、一直線上にない3点以上が与えられた場合、式(8)から1組の係数\(a,b\)を出すことは出来ません。

それは、3点以上が与えられた場合、式(5)の各項は下記のようになり、

\[Y=\left[ \begin{matrix} y_1 \\ y_2 \\ y_3 \end{matrix} \right] \qquad (10)\]

\[A=\left[ \begin{matrix} a \\ b \end{matrix} \right]\qquad (11)\]

\[X=\left[ \begin{matrix} { x }_{ 1 } & 1 \\ { x }_{ 2 } & 1 \\ { x }_{ 3 } & 1 \end{matrix} \right] \qquad (12)\]

\(X\)は正方行列ではなくなり、逆行列をもてないからです。

しかし、行列演算をうまく行うと、式(10)、式(11)、式(12)で構成された\(Y=XA\)から、一組の\(A\)を出すことができるのです。

その方法は\(A\)の前にある行列項を、行列演算して正方行列に変換することです。

式(5)の両辺に左から、\(X\)の転置行列\(X^T\)を乗じます。
\[X^TY= X^TXA\qquad (13)\]

\(X^TX\)は正方行列になるので、これが正則行列であるならば、逆行列\((X^TX)^{-1}\)が存在します。逆行列を式(13)の左から掛けると
\[(X^TX)^{-1}X^TY= (X^TX)^{-1}X^TXA\qquad (14)\]
右辺\(A\)の前の行列演算結果は単位行列\(E\)となるので、次式が得られます。
\[(X^TX)^{-1}X^TY = A \qquad (15)\]

このようにして一組の係数\(A\)が得られました。

このように逆行列を持てない行列に対し、転置行列を乗じて正方行列に転換し、逆行列を作成することを擬似逆行列と呼びます。

ううう~ん、擬似逆行列、恐ろしや、3点以上を直線式に代入しても、一つの直線式を算出することができたのです。

ところで、この直線は何を意味しているのでしょうか?

これこそが代入した全ての点の回帰直線、即ち、点から直線までの距離の二乗の和が最も小さい直線なのです。

図1.回帰直線

私はこの手法を知った時に、数学の美しさというか完璧さというか一貫性というか、そういうものに感動しました。

回帰直線は、最小二乗法を用いて算出することがよく行われていますが、その導出過程の考え方の煩雑さに対し、擬似逆行列では式(5)から\(A\)を算出するという非常にシンプルな考え方になっています。