最小二乗法

$x_1, \ldots, x_N$についての連立方程式

$$\begin{align*} y_1 &= a^1_1 x_1 + \cdots + a^N_1 x_N, \\ &\ \ \vdots \\ y_M &= a^1_M x_1 + \cdots + a^N_M x_N \end{align*}$$

i.e.

$$\begin{pmatrix} y_1 \\ \vdots \\ y_M \end{pmatrix} = \begin{pmatrix} a_1^1 & \cdots & a_1^N \\ \vdots & \ddots & \vdots \\ a_M^1 & \cdots & a_M^N \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_N \end{pmatrix}.$$

を解きたいとする。$y = (y_1 \cdots y_M)^t, \ A = (a_m^n)_{1\leq n \leq N, 1 \leq m \leq M}^t$とすれば上式は

$$y = A x$$

とかける。さて、この連立方程式は一般に解が一意に定まらない。そこで我々は、

$$\begin{align} J = \|y - Ax\|^2 \ \ \ \ \mathrm{and} \ \ \ \ \|x\|^2 \end{align}$$

が最小となるような$x$を求める。このような$x$を方程式$y = Ax$の最小二乗解という。

最小二乗解の幾何的解釈

第一の条件を満たす$x$は、幾何的には「$y \in \mathbf{R}^M$と$Ax \in \mathbf{R}^M$との距離が最小となる$x \in \mathbf{R}^N$」と解釈できる。さらに幾何的な解釈を推し進めれば、

「$y \in \mathbf{R}^M$を$\mathbf{R}^M$の部分空間$\mathrm{Im} A = \{Az;\ z \in \mathbf{R}^N\}$に直交射影したベクトルが$Ax$と一致するような$x \in \mathbf{R}^N$」

である。

なお、この解釈は次のように計算によって直接たしかめることができる; $A$の終域$\mathbf{R}^M$は、$\mathrm{Im} A$を用いて

$$\mathbf{R}^M = \mathrm{Im} A \oplus (\mathrm{Im} A)^\perp$$

と直和分解できる。$\mathrm{Im}A$への直交射影を$P$とすれば

$$\begin{align} J &= \|(1 - P)y + (Py - A x) \|^2 \\ &= \|(1 - P)y\|^2 + \|Py - A x \|^2 \\ \end{align}$$

であるから（$(1-P)y \perp Py - A x$に注意）、求める$x$は$\|Py - A x\|^2$を最小化する、すなわち

$$Py = A x$$

を満たす$x \in \mathbf{R}^{N}$である。$P$と$\mathrm{Im} A$の定義から、このような$x$は常に存在する。

$A$の特異値分解

一般に零でない行列$A \in \mathrm{Mat}(M,N)$に対して

$$Av = \sigma u,\ \ \ \ A^t u = \sigma v$$

となる$\sigma > 0, u \in \mathbf{R}^M, v \in \mathbf{R}^{N}$が存在するとき、$\sigma$を$u,v$に対する$A$の特異値といい、$u,v$をそれぞれ$\sigma$に対する$A$の左特異ベクトル、右特異ベクトルという。

さて、$\mathrm{rank}A=R(\leq\min\{M,N\})$とすると、$\mathrm{Im}A, \mathrm{Im}A^t$はそれぞれ$\mathbf{R}^M, \mathbf{R}^{N}$の$R$次元部分空間である。さらに、それぞれの正規直交基底$\{u_r\} \subset \mathbf{R}^M, \{v_r\} \subset \mathbf{R}^{N}$で、それぞれを$A$の左特異ベクトル、右特異ベクトルとし、非負実数$\{\sigma_r\}_{r = 1}^R$を特異値とするものが存在する。すなわち、

$$Av_r = \sigma_r u_r, \ \ \ \ A^t u_r = \sigma_r v_r$$

となるものが存在する。上式の第一式から行列$A$の特異値分解

$$\begin{align} A = \sum_{r=1}^R \sigma_r u_r v_r^t \end{align}$$

を得ることができる。なお、上式は$U = (u_1 \cdots u_R), V = (v_1 \cdots v_R)$とおけば

$$A = U \mathrm{diag}(\sigma_1, \ldots, \sigma_R) V^t$$

と書ける。

$\mathrm{Im}A$への直交射影$P$

特異値分解(6)を用いれば、非零行列に対して(Moore-Penrose型)一般逆行列$A^-$を

$$A^- = \sum_r \dfrac1\sigma_r v_r u_r^t$$

で定義できる。簡単な計算から

$$AA^- = \sum_r u_r u_r^t (= U U^t)$$

である。これは$\mathbf{R}^M$から$\mathrm{Im}A$への直交射影$P$に等しい。実際、任意の$\phi \in \mathbf{R}^M$は、$\mathrm{Im} A$の正規直交基底$\{u_r\}$を適当に拡張することで得られる$\mathbf{R}^M$の正規直交基底$\{u_m\}_{m = 1}^M$と適当な係数$a_m\in \mathbf{R}$を用いて$\phi = \sum_m a_m u_m$と書けるので、

$$A A^- \phi = \sum_r a_r u_r \in \mathrm{Im} A$$

であり、かつ定義から明らかに$(A A^-)^2 = (A A^-)^t = A A^-$が成立する。

以降

$$\begin{align} P = A A^- \end{align}$$

とする。

ふたたび、最小二乗解

我々は、与えられた$y \in \mathbf{R}^M$と$A \in \mathrm{Mat}(M,N)$に対して$\|Py - A x\|^2$を最小化する$x \in \mathbf{R}^{N}$を求めたいのであった。$P$が$\mathrm{Im} A$への直交射影であることと、$\mathrm{Im} A$の定義から明らかに$Py = A x$となる$x$が存在する。従って式(6), (7)より

$$Py = \sum_r \langle u_r, y \rangle u_r = \sum_r \sigma_r \langle v_r, x \rangle u_r = A x$$

となり、$r = 1, \ldots, R$に対して

$$\langle v_r, x \rangle = \dfrac{\langle u_r, y \rangle}{\sigma_r}$$

が成立する。$\mathrm{Im} A^t$の正規直交基底$\{v_r\}_{r=1}^R$を適当に拡張した$\{v_1, \ldots, v_R, v_{R+1}, \ldots, v_{N}\}$を$\mathbf{R}^{N}$の正規直交基底とすれば、$x$のFourier展開は

$$\begin{align} x &= \sum_n \langle v_n, x \rangle v_n \\ &= \sum_{r=1}^R \dfrac{\langle u_r, y \rangle}{\sigma_r} v_r + \sum_{i = R+1}^N \langle v_i, x \rangle v_i \end{align}$$

となる。いま、$\| x \|^2$が最小となるようにするので、上式第二項は$0$となるようにしておけばよい。この条件の下で上式は

$$\begin{align} x &= \sum_{r=1}^R \dfrac{\langle u_r, y \rangle}{\sigma_r} v_r \\ &= \left( \sum_r \dfrac1\sigma_r v_r u_r^t \right) y \\ &= A^- y \end{align}$$

となる。これが求める$x$である。

特殊な場合の最小二乗解

$M > N = \mathrm{rank}A$の場合

この場合、$N \times N$行列$A^t A$の逆行列$(A^tA)^{-1}$が存在する。そして$A$の特異値分解をもちいれば

$$\begin{align*} (A^tA)^{-1}A^t &= \left( \sum_r \sigma_r v_r u_r^t \sum_{r'} \sigma_{r'} u_{r'} v_{r'}^t \right)^{-1} A^t \\ &= \left( \sum_r \dfrac{1}{\sigma_r^2} v_r v_r^t \right) A^t \\ &= \sum_{r} \dfrac{1}{\sigma_{r}^2} v_{r} v_{r}^t \sum_{r'} \sigma_{r'} v_{r'} u_{r'}^t \\ &= \sum_r \dfrac{1}{\sigma_r} v_r u_r^t \\ &= A^- \end{align*}$$

であるから、方程式$y=Ax$の最小二乗解$x$は

$$\begin{align} x = (A^tA)^{-1}A^t y \end{align}$$

となる。

また、$\mathrm{Im} A$ への直交射影 $P$ は

$$P = A A^- = A (A^t A)^{-1} A^t$$

となる。

$N > M = \mathrm{rank} A$の場合

この場合、$M \times M$行列$A A^t$の逆行列$(A A^t)^{-1}$が存在する。そして$A$の特異値分解をもちいれば

$$\begin{align*} A^t (A A^t)^{-1} &= \sum_r \sigma_r v_r u_r^t \sum_{r'} \dfrac{1}{\sigma_{r'}^2} u_{r'} u_{r'}^t \\ &= \sum_r \dfrac{1}{\sigma_r} v_r u_r^t \\ &= A^- \end{align*}$$

であるから、方程式$y=Ax$の最小二乗解$x$は

$$x = A^t (A A^t)^{-1} y$$

となる。