最小二乗法(least−squares fitting)

(x,f(x))なる点列にフィットする線(linear fitting of points(x,f(x))

例: $\{(x_{i},y_{i})\}_{i=1}^{m}$ で表される点列を $y=Ax+B$ でfittingさせることを考える。差の二乗を考えてそれが最小となる $A,B$ を求める。この「差の二乗」を考えるというのが、この方法の名前たるゆえんである。

$E(A,B)=\displaystyle \sum_{i=1}^{m}\left[\left(Ax_{i}+B\right)-y_{i}\right]^{2}$

最小値は、 $(0,0)=\nabla E=$ (合成関数の微分法により) $2\displaystyle \sum_{i=1}^{m}\left[\left(Ax_{i}+B\right)-y_{i}\right](x_{i},1)$

行列表現で書けば $\left[\begin{array}{l}\sum_{i=1}^{m}x_{i}^{2} \sum_{i=1}^{m}x_{i}\\\sum_{i=1}^{m}x_{i} 1\end{array}\right]\left[\begin{array}{l}A\\B\end{array}\right]=\left[\begin{array}{l}\sum_{i=1}^{m}x_{i}y_{i}\\\sum_{i=1}^{m}y_{i}\end{array}\right]$

あとは逆行列つかえば $A,B$ が求まる。

直交回帰を用いた点列にフィットする線(linear fitting of points using orthogonal regression)

以下の議論ではn次元での点列と線を考える。

直線: $\vec{L}(t)=\vec{A}+t\vec{D} , \vec{D}$ は単位ベクトル

点列: $\vec{X}_{i}=\vec{A}+d_{i}\vec{D}+p_{i}\vec{D}_{i}\perp$ , $d_{i}=\vec{D}\text{・}(\vec{X}_{i}-\vec{A})$ , $\vec{D}_{i}\perp$ は $\vec{D}$ と直交する単位ベクトル, $p_{i}$ は適当な係数

$\vec{A}$ を起点として考えると計算がすんなりいきそうなので、 $\vec{Y}_{i}=\vec{X}_{i}-\vec{A}$ とおく。

こいつを使って、 $\vec{X}_{i}$ からこの直線 $\vec{L}(t)$ を法線とする平面への投影は

$\vec{Y}_{i}-d_{i}\vec{D}=p_{i}\vec{D}_{i}\perp$

このベクトルの長さの二乗は $p_{i}^{2}=(\vec{Y}_{i}-d_{i}D)^{2}$

よって最小化すべきエネルギー関数は $E(\displaystyle \vec{A},\vec{D})=\sum_{i=1}^{m}p_{i}^{2}$

これは次の二通りの書き方が出来る。

１．

$E(\displaystyle \vec{A},\vec{D})=\sum_{i=1}^{m}(\vec{Y}_{i}^{T}\left[I-\vec{D}\vec{D}^{T}\right]\vec{Y}_{i})$

２．

$\vec{E}(\vec{A},\vec{D})=\vec{D}^{T}\left(\sum_{i=1}^{m}\left[\left(\vec{Y}_{i}\text{・}\vec{Y}_{i}\right)I - \vec{Y}_{i}\vec{Y}_{i}^{T}\right]\right)\vec{D}=\vec{D}^{T}M(A)\vec{D}$

ちょっと補足が必要だろう。

１．は、 $d_{i}=\vec{D}\text{・}(\vec{X}_{i}-\vec{A})=\vec{D}\text{・}\vec{Y}$ を $p_{i}^{2}=(\vec{Y}_{i}-d_{i}\vec{D})^{2}$ に代入し、 $\vec{D}$ , $\vec{Y}_{i}$ だけで表したもの。

この式変形がわかりにくい人のために少しベクトルの計算をおさらいしておこう。ベクトルは行列の一種なので行列で成り立つものはベクトルでも成り立つが、ベクトル特有の法則がある。

1. $\left|\vec{x}\right|^{2}=\vec{x}\text{・}\vec{x}=\vec{x}^{T}\vec{x}$ (ただしxの各要素は実数であるとする→複素数に拡張するなら、転置ではなく、複素共役転置ベクトル $\vec{x}^{*}$ を用いて $\left|\vec{x}\right|^{2}=\vec{x}^{*}\vec{x}$ と書く。複素共役転置ベクトルについても以下の3.の転置に公式とかが使える。しかし、当分関係ないのでこれ以上詳しくは立ち入らない。)

これを使えば、ベクトルの内積を行列の掛け算として書ける。行列の掛け算に変形すれば、行列の掛け算のさまざまな公式を使うことが出来るようになる。行列の掛け算に関する公式は以下のようなもの。

2. $ABC=A(BC)=(AB)C$ ただし $A,B,C$ は行列。

結合則。

3. $(AB)^{T}=B^{T}A^{T}$ , $(A+B)^{T}=A^{T}+B^{T}$ , $(A^{T})^{T}=A$

転置してあるものを展開するのに重要な公式である。

4. $A=IA=AI , II=I , I^{T}=I$ , ただし $I$ は単位行列

$I$ は数字の１の行列のやつだと思うとわかりやすい。

5. $cAB=(cA)B=A(cB)$ ただし cはスカラー。

行列の掛け算に対して交換則。たとえば、 $\vec{x}^{T}\vec{y}$ はスカラーなのでこの部分に関しては交換則が成り立つ。例) $\vec{x}^{T}\vec{y}\vec{z}=\vec{z}\vec{x}^{T}\vec{y}$

ただし、行列においては一般的には交換則は成り立たないことに注意。

6. $\vec{x}^{T}\vec{y}=\vec{y}^{T}\vec{x}$

$\vec{x}^{T}\vec{y}$ の形がベクトルの内積であることから、交換則が成り立つのは明らか。

7. $(A+B)(C+D)=AC+AD+BC+BD$

分配則。このへんは、実数の時のと同じなので説明はいらんでしょう。

よく使うのはそんだけ。これらを使うと、たとえば、

ベクトルの内積の2乗: $(\vec{x}\text{・}\vec{y})^{2}=(\vec{x}^{T}\vec{y})^{2}=(\vec{x}^{T}\vec{y})^{T}(\vec{x}^{T}\vec{y})=\vec{y}^{T}\vec{x}\vec{x}^{T}\vec{y}=\vec{y}^{T}\left(\vec{x}\vec{x}^{T}\right)\vec{y}$

などと自在に変形できたりする。行列の結合則から、この $\left(\vec{x}\vec{x}^{T}\right)$ の部分を先に計算することも許されることになる。この $\left(\vec{x}\vec{x}^{T}\right)$ はなんじゃい、というと、これは $\vec{x}$ が3次元のベクトルなら、 $\left(\vec{x}\vec{x}^{T}\right)$ は3×3の行列だ。テンソルなんて呼ばれることもある。

もうちょい練習。

ベクトルの和の2乗: $(\vec{x}+\vec{y})^{2}=(\vec{x}+\vec{y})^{T}(\vec{x}+\vec{y})=(\vec{x}^{T}+\vec{y}^{T})(\vec{x}+\vec{y})=\vec{x}^{T}\vec{x}+2\vec{x}^{T}\vec{y}+\vec{y}^{T}\vec{y}$

行列の積の2乗: $(AB)^{2}=(AB)^{T}(AB)=B^{T}A^{T}AB$

行列の和の2乗: $(A+B)^{2}=(A+B)^{T}(A+B)=(A^{T}+B^{T})(A+B)=(A^{T}A+B^{T}A+A^{T}B+B^{T}B)$

ちょっと前置きが長くなったが、問題の部分、計算してみよう。

$p_{i}^{2}=(\vec{Y}_{i}-d_{i}\vec{D})^{2}=(\vec{Y}_{i}-\vec{D}^{T}\vec{Y}_{i}\vec{D})^{2}=(\vec{Y}_{i}-\vec{D}\vec{D}^{T}\vec{Y}{}_{i})^{2}=\left(\left(I-\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}\right)^{2}$

行列の積の2乗の形になったので、↑の公式で展開してみる。

$=\vec{Y}_{i}^{T}\left(I-\vec{D}\vec{D}^{T}\right)^{T}\left(I-\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}$
ところで $\left(I-\vec{D}\vec{D}^{T}\right)^{T}=I^{T}-\left(\vec{D}\vec{D}^{T}\right)^{T}=I-\vec{D}\vec{D}^{T}$ なので

$=\vec{Y}_{i}^{T}\left(I-\vec{D}\vec{D}^{T}\right)\left(I-\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}$

行列の分配則で
$=\vec{Y}_{i}^{T}\left(I-2\vec{D}\vec{D}^{T}+\vec{D}\vec{D}^{T}\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}$
ところで、 $\vec{D}\vec{D}^{T}\vec{D}\vec{D}^{T}$ の部分、結合則から $\vec{D}\left(\vec{D}^{T}\vec{D}\right)\vec{D}^{T}$ と計算してもいいはずだ。それで、 $\vec{D}^{T}\vec{D}$ は実のところ内積なので $=\vec{D}^{2}=1$ ( $\vec{D}$ は単位ベクトルなので)だったりする。となれば、

$=\vec{Y}_{i}^{T}\left(I-2\vec{D}\vec{D}^{T}+\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}=\vec{Y}_{i}^{T}\left(I+\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}$
てなことになる。

２．のほうは、１．からの変形だ。 $\vec{D}$ を $\displaystyle \sum$ の外に追い出すように外側に来るように変形する。

$\vec{Y}_{i}^{T}\left(I+\vec{D}\vec{D}^{T}\right)\vec{Y}_{i}=\vec{Y}_{i}^{T}\vec{Y}_{i}+\vec{Y}_{i}^{T}\vec{D}\vec{D}^{T}\vec{Y}_{i}=\vec{Y}_{i}^{T}\vec{Y}_{i}+\vec{D}^{T}\vec{Y}_{i}\vec{Y}_{i}^{T}\vec{D}$

$\vec{D}^{T}$ (...) $\vec{D}$ の形でくくってしまいたいのだが、 $\vec{Y}_{i}^{T}\vec{Y}_{i}$ が邪魔だ。 $\vec{D}^{T}\vec{D}=1$ なので、これをかけても値は変わらない。 $\vec{Y}_{i}^{T}\vec{Y}_{i}=\vec{D}^{T}\vec{D}\vec{Y}_{i}^{T}\vec{Y}_{i}$ あと $\vec{Y}_{i}^{T}\vec{Y}_{i}$ はスカラーであることから、この部分で交換則が成り立って、 $\vec{D}^{T}\left(\vec{Y}_{i}^{T}\vec{Y}_{i}I\right)\vec{D}$ とか出来てしまう。よって、

$E(\displaystyle \vec{A},\vec{D})=\sum_{i=1}^{m}(\vec{Y}_{i}^{T}\left[I-\vec{D}\vec{D}^{T}\right]\vec{Y}_{i})=\vec{D}^{T}\left(\sum_{i=1}^{m}\left[\left(\vec{Y}_{i}^{T}\vec{Y}_{i}\right)I - \vec{Y}_{i}\vec{Y}_{i}^{T}\right]\right)\vec{D}$

あと、この $\displaystyle \sum$ の中身はAをパラメータとする行列なのでM(A)と置くと

$=\vec{D}^{T}M(A)\vec{D}$

さて、式変形が終わったところで計算していこう。

１．は $\vec{A}$ で微分すると

$\displaystyle \frac{\partial E}{\partial A}=\frac{\partial Y}{\partial A}\frac{\partial E}{\partial Y}=-2\left[I-\vec{D}\vec{D}^{T}\right]\sum_{i=1}^{m}\vec{Y}_{i}$

これは $\displaystyle \sum_{i=1}^{m}\vec{Y}_{i}=0$ のときは、常に偏微分が0になる。つまり、 $\displaystyle \vec{A}=\frac{1}{m}\sum_{i=1}^{m}\vec{X}_{i}$ , 早い話が点列の重心。

※　このベクトルの微分、少し説明がいるかも知れない。詳しいことはベクトル解析の教科書を見ればいいのだが、 $\displaystyle \frac{\partial}{\partial Y}\left(Y^{T}Y\right)=2Y$ $,\displaystyle \frac{\partial}{\partial Y}\left(aC^{T}Y\right)=aC$ , $\displaystyle \frac{\partial}{\partial Y}\left(Y^{T}(aI)Y\right)=2aY,\frac{\partial}{\partial Y}\left(Y^{T}(IC)Y\right)=2CY,\frac{\partial}{\partial Y}\left(Y^{T}AY\right)=2AY,\frac{\partial}{\partial Y}\left(C^{T}AY\right)=AC$ (ただし $a$ はスカラ, $C$ はベクトル, $A$ は対称行列)などと $\displaystyle \frac{d}{dx}(x^{2})=2x$ とか $\displaystyle \frac{d}{dx}(ax^{2})=2ax$ と同じ感覚で計算できることを言っておこう。あと、 $\vec{D}\vec{D}^{T}$ は対称行列であることにも注意か。

２．の式で $\vec{A}$ を与えると $M(A)$ が決まる。 $\vec{D}^{T}M(A)\vec{D}$ は二次形式なので、この最小値は、 $M(A)$ のもっとも小さい固有値である。

※　 $\displaystyle \vec{x}^{T}A\vec{x}=\sum_{i,j=1}^{n}a_{ij}x_{i}x_{j}$ これを行列Aに関する行列の二次形式(quadratic form)という。「主軸変換」でよく出てくる形だ。計算が簡略になるのでふつうは $A$ が対称行列と仮定する。これを仮定しても一般性を失わないし、これを仮定しないと行列からめて微分する時にややこしいからである。(さっきの $\displaystyle \frac{\partial}{\partial Y}\left(Y^{T}AY\right)=2AY$ みたいなの) またこれは $A$ の固有ベクトルの方向が主軸となるような回転変換(鏡像変換もこれに含む)である。詳しいことは、線形代数の教科書に書いてある。

$n=2,\vec{A}=(a,b)$ のとき、 $M(A)$ は

$M(A)=\left(\sum_{i=1}^{m}(x_{i}-a)^{2}+\sum_{i=1}^{n}(y_{i}-b)^{2}\right)\left[\begin{array}{l}1 0\\0 1\end{array}\right]-\left[\begin{array}{l}\sum_{i=1}^{m}(x_{i}-a)^{2} \sum_{i=1}^{m}(x_{i}-a)(y_{i}-b)\\\sum_{i=1}^{m}(x_{i}-a)(y_{i}-b) \sum_{i=1}^{m}(y_{i}-b)^{2}\end{array}\right]$

$n=3$ のときは省略。

直交回帰を用いた点列にフィットする超平面(hyperplanar fitting of points using orthogonal regression)

同様にn次元での平面でフィッティングさせてみる

平面: $\vec{N}\text{・}(\vec{X}-\vec{A})=0 , \vec{N}:\text{単位ベクトル}$

点列: $\vec{X}_{i}=\vec{A}+\lambda_{i}\vec{N}+p_{i}\vec{N}_{i}\perp$ , $\lambda_{i}=\vec{N}\text{・}(\vec{X}_{i}-\vec{A})$ , $\vec{N}_{i}\perp$ は $\vec{N}$ と直交する単位ベクトル, $p_{i}$ は適当な係数

よって $\lambda_{i}^{2}=(\vec{N}\text{・}\vec{Y}_{i})^{2}$

最小二乗法のためのエネルギー関数 $E(\displaystyle \vec{A},\vec{N})=\sum_{i=1}^{m}\lambda_{i}^{2}$

さっきの議論と同様、二通りに書けて

１． $E(\displaystyle \vec{A},\vec{N})=\sum_{i=1}^{m}\left(\vec{Y}_{i}^{T}\left[\vec{N}\vec{N}^{T}\right]\vec{Y}_{i}\right)$

２． $E(\vec{A},\vec{N})=\vec{N}^{T}\left(\sum_{i=1}^{m}\vec{Y}_{i}\vec{Y}_{i}^{T}\right)\vec{N}$ $= \vec{N}^{T}M(A)\vec{N}$

１．は $\vec{A}$ で微分すると

$\displaystyle \frac{\partial E}{\partial A}=-2\left[\vec{N}\vec{N}^{T}\right]\sum_{i=1}^{m}\vec{Y}_{i}$

以下、線にフィットさせる時と同様の議論なので割愛。

円の2次元の点列へのフィッティング(fitting a circle to 2D points)

点列: $\left\{(x_{i},y_{i})\right\}_{i=1}^{m},m\geq 3$ 　, ただしすべての点が同一直線上(collinear)にはないとする。

円: $(x-a)^{2}+(y-b)^{2}=r^{2}$

エネルギー関数は
$E(a,b,r)=\displaystyle \sum_{i=1}^{m}(L_{i}-r)^{2} , L_{i}=\sqrt{(x_{i}-a)^{2}+(y_{i}-b)^{2}}$

$\displaystyle \frac{\partial E}{\partial r}=-2\sum_{i=1}^{m}(L_{i}-r)$ これが0になるのは、 $r=\displaystyle \frac{1}{m}\sum_{i=1}^{m}L_{i}$ のとき。

$\displaystyle \frac{\partial E}{\partial a}=-2\sum_{i=1}^{m}(L_{i}-r)\frac{\partial L_{i}}{\partial a}=2\sum_{i=1}^{m}\left((x_{i}-a)+r\frac{\partial L_{i}}{\partial a}\right)$

$\displaystyle \frac{\partial E}{\partial b}=-2\sum_{i=1}^{m}(L_{i}-r)\frac{\partial L_{i}}{\partial b}=2\sum_{i=1}^{m}\left((y_{i}-b)+r\frac{\partial L_{i}}{\partial b}\right)$

これらが0になるのは、
$a = \displaystyle \frac{1}{m}\sum_{i=1}^{m}x_{i}+r\frac{1}{m}\sum_{i=1}^{m}\frac{\partial L_{i}}{\partial a}$
$b = \displaystyle \frac{1}{m}\sum_{i=1}^{m}y_{i}+r\frac{1}{m}\sum_{i=1}^{m}\frac{\partial L_{i}}{\partial b}$
のとき。