忘れた時に振り返ることができるように最小二乗法によるパラメータ推定についてまとめた。

問題設定

$n$ このデータ点 $(x_{1}, y_{1}, z_{1}), (x_{2}, y_{2}, z_{2}), \cdots , (x_{n}, y_{n}, z_{n})$ が得られたとき、このデータ点を $z=a+bx+cy$ という関係式で表したい（すなわち $a, b, c$ の値を推定したい）。これを行列表記すると

$\begin{eqnarray} \left( \begin{array}{c} z_{1} \\ z_{2} \\ \vdots \\ \\ z_{n} \end{array} \right) &=& \left( \begin{array}{ccc} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ \vdots \\ \\ 1 & x_{n} & y_{n} \\ \end{array} \right) \left( \begin{array}{c} a \\ b \\ c \end{array} \right) + \vec{\varepsilon} \\ \vec{z} &=& A\vec{\theta} + \vec{\varepsilon} \end{eqnarray}$

となり、 $\vec{\varepsilon}$ が最小となるようなパラメータ $\vec{\theta}$ を求めたい。

パラメータ導出

$|\vec{\varepsilon}|^{2}$ を最小にする $\vec{\theta}$ は

$\begin{eqnarray} \vec{0} &=& \frac{\partial}{\partial \vec{\theta}} |\vec{\varepsilon}|^{2} \\ &=& \frac{\partial}{\partial \vec{\theta}} | \vec{z}-A\vec{\theta} |^{2} \\ &=& \frac{\partial}{\partial \vec{\theta}} (\vec{\theta}^\mathrm{T}A^\mathrm{T}A\vec{\theta} - 2\vec{z}^\mathrm{T}A\vec{\theta}) \\ &=& 2\vec{\theta}^\mathrm{T}A^\mathrm{T}A -2\vec{z}^\mathrm{T}A \\ \end{eqnarray}$

別の書き方をしてみると

$\begin{eqnarray} |\vec{\varepsilon}|^{2} &=& \sum_{i} e_{i}^{2} \\ \vec{0} &=& \frac{\partial |\vec{\varepsilon}|^{2}}{\partial \vec{\theta}} \\ &=& 2\sum_{i} e_{i} \frac{\partial e_{i}}{\partial \vec{\theta}} \end{eqnarray}$

途中で数学公式まとめの「ベクトルでの微分」を参考にしている。転置をとると

$\begin{eqnarray} A^\mathrm{T}A\vec{\theta} &=& A^\mathrm{T}\vec{z} \\ \vec{\theta} &=& (A^\mathrm{T}A)^{-1}A^\mathrm{T}\vec{z} \end{eqnarray}$

と求まる。これをnormal equationと呼ぶ。

reduced QR分解を使った解法

QR分解については数学公式まとめを参照。
$m\times n (m \ge n)$ 行列 $A$ を $A=QR$ とreduced QR分解（ $Q$ は $m \times n$ 行列、 $R$ は $n \times n$ 行列）できたとする。この場合、normal equationは

$\begin{eqnarray} (QR)^\mathrm{T}(QR)\vec{\theta} &=& (QR)^\mathrm{T} \vec{z} \\ R^\mathrm{T}Q^\mathrm{T} QR \vec{\theta} &=& R^\mathrm{T}Q^\mathrm{T} \vec{z} \\ \end{eqnarray}$

$Q^\mathrm{T}Q$ は単位行列になる（ただし $Q$ は正方行列ではないので直交行列とは呼ばない）。したがって

$\begin{eqnarray} R^\mathrm{T}R \vec{\theta} &=& R^\mathrm{T}Q^\mathrm{T} \vec{z} \end{eqnarray}$

Rは正方行列であり、 $A$ がフルランク行列（ $\mathrm{rank}(A)=n$ ）であることは前提なので逆行列 $R^{-1}$ が存在する。両辺に左から $(R^{-1})^\mathrm{T}$ をかけると

$\begin{eqnarray} (R^{-1})^\mathrm{T}R^\mathrm{T}R \vec{\theta} &=& (R^{-1})^\mathrm{T}R^\mathrm{T}Q^\mathrm{T} \vec{z} \end{eqnarray}$

数学公式まとめより、 $(R^{-1})^\mathrm{T}=(R^\mathrm{T})^{-1}$ なので、

$\begin{eqnarray} (R^\mathrm{T})^\mathrm{-1}R^\mathrm{T}R \vec{\theta} &=& (R^{T})^\mathrm{-1}R^\mathrm{T}Q^\mathrm{T} \vec{z} \\ R \vec{\theta} &=& Q^\mathrm{T} \vec{z} \\ QR \vec{\theta} &=& QQ^\mathrm{T} \vec{z} \\ \vec{\theta} &=& (QR)^{-1}\vec{z} \end{eqnarray}$

このことは、最小二乗法の解はnormal equationではなく $A\vec{\theta} = \vec{z}$ の方程式をQR分解で解けば良いことを示している。

参考

mathtrain.jp

ロバスト推定

上記の最小二乗法では誤差の重みを全て同等に扱っているが、これだと回帰値から大きく外れた点が $|\vec{\varepsilon}|^{2}$ に与える寄与が大きくなりアウトライヤー点の影響を受けやすくなる。そこで、そのような点の寄与を小さくするために重みを下げる方法がある。各誤差の寄与 $\rho$ を誤差の絶対値の関数として決めることにすると

$\begin{eqnarray} |\vec{\varepsilon}^{'}|^{2} &=& \sum_{i} \rho(e_{i}) \end{eqnarray}$

ここで、通常の最小二乗法では $\rho(e_{i})=e_{i}^{2}$ であることに注意せよ。
これを微分していくと

$\begin{eqnarray} \vec{0} &=& \sum_{i} \frac{\partial \rho}{\partial e_{i}}\frac{\partial e_{i}}{\partial \vec{\theta}} \\ &=& \sum_{i} \left[ \frac{1}{e_{i}} \frac{\partial \rho}{\partial e_{i}} \right] e_{i} \frac{\partial e_{i}}{\partial \vec{\theta}} \\ \end{eqnarray}$

この表式から、 $\frac{1}{e_{i}}\frac{\partial \rho}{\partial e_{i}}$ を各誤差に対する重みと考えることができる。
具体例として、Huber関数を考える。Huber関数とは、誤差関数がある正の定数 $k$ を用いて

$\begin{eqnarray} \rho(e) &=& \frac{1}{2} e^{2} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (|e| \lt k) \\ &=& k|e|-\frac{1}{2}k^{2} \ \ \ \ \ (|e| \ge k) \end{eqnarray}$

と定義される関数である。この重みは

$\begin{eqnarray} \frac{1}{e}\frac{\partial \rho(e)}{\partial e} &=& 1 \ \ \ \ \ \ \ \ \ (|e| \lt k) \\ &=& \frac{k}{|e|} \ \ \ \ \ (|e| \ge k) \\ \end{eqnarray}$