一般逆行列を使った最小二乗法 の履歴(No.1)

一般逆行列を使った最小二乗法†

$\def\bm#1{{\boldsymbol #1}}$

最小二乗法によるデータフィッティング†

例えば多数のデータ点 $(x_1, y_1), (x_2,y_2),\dots$ が与えられて、これを二次関数 $y=Ax^2+Bx+C$ でフィッティングしたい、という状況を考える。

未知数は $A,B,C$ の３つなので、３つの独立な点が与えられればそれらを通るようにパラメータを決められるのだけれど、ここではデータ点にはノイズが含まれているので、３つよりも多くの点が合が得られた上で、それらのデータ点をよりよく再現するパラメータ $A,B,C$ を求めたい。

具体的には誤差の二乗 $$ \sum_{n=1}^N (Ax_n^2+Bx_n+C-y_n)^2 $$ を最小化するよう $A,B,C$ を決める。

$A,B,C$ で偏微分して、、、というのも１つの手なのだけれど、ここではこの問題を線形代数で解く。

行列方程式に直す†

誤差が全くない場合、すべての $n$ に対して $$ Ax_n^2+Bx_n+C=y_n $$ が成り立つ。これを $$ \begin{pmatrix} x_1^2&x_1&1\\ x_2^2&x_2&1\\ \vdots&\vdots&\vdots\\ x_N^2&x_N&1\\ \end{pmatrix}\begin{pmatrix} A\\B\\C \end{pmatrix} =\begin{pmatrix} y_1\\y_2\\\vdots\\y_N \end{pmatrix} $$ と書こう。ここでは $A,B,C$ が未知数であることを考えると、この式は $$ A\bm x=\bm b $$ の形をしているが、 $N=3$ の場合を除いて $A$ は正則でないばかりか正方行列ですらないため、逆行列を使って $$ \bm x=A^{-1}\bm b $$ というように値を決めることはできないし、そもそもそのような $\bm x$ は必ずしも存在しない。最小二乗法は $$ A\bm x=\bm b $$ を目指すのではなく、 $$ \|A\bm x-\bm b\|^2 $$ を最小化することを目指すのである。実際、上で示した二乗誤差の式はこれらの変数を使えば $$ \sum_{n=1}^N (Ax_n^2+Bx_n+C-y_n)^2=\|A\bm x-\bm b\|^2 $$ と表される。

問題の一般化と適用限界†

ここまでの話を一般化すれば、パラメータは $x$ ひとつだけでなく複数次元のパラメータ空間を考えてもいいし、フィッティングもパラメータの単なる多項式で表すのではなく指数関数等複雑な関数を含むようにしてもかまわない。またデータ点も $y$ という１つの数値だけでなく複数次元空間中の点であってもよい。

最終的にパラメータから求めた何らかの値を並べた $A$ と係数ベクトル $\bm x$ との積でデータ点を表す形でありさえすれば以下の話をそのまま適用可能である。

一方、フィッティング関数の中に各項の「係数」以外に未知数があってはならないというのは重要な適用限界と言えそう。

例えば指数関数の減衰定数をこの方法でフィッティングすることはできない気がする。

一般逆行列

上述をまとめると、$A$ が正則でないとき $A\bm x=\bm b$ を解いて $\bm x=A^{-1}\bm b$ とはできないので、代わりに $\|A\bm x-\bm b\|$ を最小化する $\bm x$ を $\bm x=A^- \bm b$ として求められるような $A^-$ を求めたい。

このとき $A$ は正則でないばかりか、そもそも正方行列でなくてもいい。

このような $A^-$ を一般逆行列と呼び、その定義は $$ AA^-A=A $$ で与えられるそうだ？

ここから、$A$ が $m\times n$ 行列なら $A^-$ は $n\times m$ 行列になることが分かる。

$A$ が正則なら $A^-=A^{-1}$ であり、上記の定義も、 $\|A\bm x-\bm b\|=0$ もどちらも満たすことに注意せよ。

一般逆行列にはいろいろな種類がある：

• 反射型
• 最小二乗型
• ノルム最小型
• ムーア・ペンローズ

以下ではムーア・ペンローズ一般逆行列について考える。

とうのも、ムーア・ペンローズ一般逆行列を使うと最小二乗解の中でノルムが最小となる解が得られるそうだ(?)

ムーア・ペンローズ一般逆行列の定義†

以下ではムーア・ペンローズ一般逆行列を $A^\dagger$ で表す。

このサイトの他のページではエルミート共役（複素共役+転置）を $A^\dagger$ で表しているが、このページでは複素共役+転置は $A^*$ で表す。

ムーア・ペンローズ一般逆行列 $A^\dagger$ は以下の４つの性質を持つ行列として定義される。

1. $AA^\dagger A=A$ 2. $A^\dagger AA^\dagger=A^\dagger$ 3. $(AA^\dagger)^*=AA^\dagger$ → ここでは $^*$ が複素共役＋転置を表す 4. $(A^\dagger A)^*=A^\dagger A$

つまり、

• $(A^\dagger)^\dagger=A$ であること
• $AA^\dagger$ も $A^\dagger A$ も自己共役($X^*=X$ を満たす)であること が条件となる。

任意の行列 $A$ に対してムーア・ペンローズ一般逆行列は一意に定まる†

なぜなら $\mathrm{rank}\,A=r$ として、$A$ の特異値標準化を非ゼロ対角行列 $D$ とユニタリ行列 $U, V$ により

$$ A=U\begin{pmatrix} D_{r\times r}&O_{r\times(N-r)}\\ O_{(M-r)\times r}&O_{(M-r)\times(N-r)} \end{pmatrix}V^* $$ と表すとき、

$$ A^\dagger=V\begin{pmatrix} D_{r\times r}^{-1}&O_{r\times(M-r)}\\ O_{(N-r)\times r}&O_{(N-r)\times(M-r)} \end{pmatrix}U^* $$ を作ればこれがムーア・ペンローズ一般逆行列となる。

最小二乗解が得られる†

$\bm x=A^\dagger \bm b$ とすることで最小二乗解が得られる。

証明は２段階で行う:
まず、正規方程式の解が最小二乗解を与えることを示し、 $\hat{\bm x}=A^\dagger\bm b$ が正規方程式の解を与えることを示す。

(1) 正規方程式の解が最小二乗解を与えること: $$ A\bm x=\bm b $$ に対して、左から $A$ の複素共役転置 $A^*$ を掛けた $$ A^*A\bm x=A^*\bm b $$ を正規方程式と呼ぶ。

この解を $\bm x'$ と書くと、 $$ A^*(A\bm x'-\bm b)=0 $$ は $A$ のすべての列ベクトルが $A\bm x-\bm b$ と直交することを示している。

すると、任意の数ベクトルを $A$ にかけたベクトルも直交する。

そこで、 $$ \begin{aligned} \|A\bm z-\bm b\|^2 &=\|A(\bm z-\bm x')+A\bm x'-\bm b\|^2\\ &=\|A(\bm z-\bm x')\|^2+\cancel{2\big(A(\bm z-\bm x'),A\bm x'-\bm b\big)}+\|A\bm x'-\bm b\|^2\\ &=\|A(\bm z-\bm x')\|^2+\|A\bm x'-\bm b\|^2\\ &\ge \|A\bm x'-\bm b\|^2\\ \end{aligned} $$ すなわち、任意のベクトル $\bm z$ に対して $\|A\bm z-\bm b\|^2$ は $\|A\bm x'-\bm b\|^2$ を下回ることはなく、 $\bm x'$ が最小二乗解を与えることが確認できた。

(2) $\hat{\bm x}=A^\dagger\bm b$ が正規方程式の解になること

確かめてみよう： $$ A^*A\hat{\bm x}=A^*\bm b $$ の右辺に代入して左辺が導かれることを示す。 $$ A^*AA^\dagger\bm b $$ $AA^\dagger=(A^\dagger A)^*$ より、 $$ A^*(AA^\dagger)^*\bm b=A^*A^{\dagger *}A^* \bm b=(AA^\dagger A)^*\bm b=A^*\bm b $$ 一番右の等号は $A^\dagger$ が $A$ のムーア・ペンローズ一般逆行列であることを示している。

実は、最小二乗解というのは一般に一意には定まらないのだけれど、ムーア・ペンローズ一般逆行列を使って求めた $\hat{\bm x}=A^\dagger\bm b$ は、そのなかでノルム $\|\bm x\|$ が最小のものを与えるのだそうだ。

すごい。

数値解法の手順†

$A$ を特異値分解 (svd: Singular Value Decomposition)　する。

$$ A=U\begin{pmatrix} D_{r\times r}&O_{r\times(N-r)}\\ O_{(M-r)\times r}&O_{(M-r)\times(N-r)} \end{pmatrix}V^* $$

例えば numpy ならこれかな？
https://numpy.org/doc/2.2/reference/generated/numpy.linalg.svd.html

特異値行列は対角行列なので、その逆行列（対角成分の逆数を取る）を求め、

$$ A^\dagger=V\begin{pmatrix} D_{r\times r}^{-1}&O_{r\times(M-r)}\\ O_{(N-r)\times r}&O_{(N-r)\times(M-r)} \end{pmatrix}U^* $$ によって $A^\dagger$ を得る。

$\hat{\bm x}=A^\dagger\bm b$ によって $\hat{\bm x}$ を求める。

コメント・質問†