💭

最小二乗法

2021/09/29に公開約5,000字

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}^MAx \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}が存在するとき、\sigmau,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}^MA \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

となる。

Discussion

ログインするとコメントできます