VMLS 12.3 最小二乗法 その三

線形代数

前回→VMLS 12.2
テキスト→Introduction to Applied Linear Algebra – Vectors, Matrices, and Least Squares

12.3 Solving least squares problems

　最小二乗法はQR分解を用いて解くことができる。{\displaystyle A} をQR分解における {\displaystyle A=QR} とする (これは12.3の仮定にもとづく)。疑似逆行列 {\displaystyle A^{\dagger} } は {\displaystyle A^{\dagger}=R^{-1}Q^T } と表されるため、次式を得る。

{\displaystyle \hat{x}=R^{-1}Q^T b. \quad(12.10) }
{\displaystyle \hat{x} } を計算するために、まず {\displaystyle Q^T } に {\displaystyle b } を掛ける。それから back substitution*1 を用いて {\displaystyle R^{-1}Q^T b } を計算する。{\displaystyle A } と {\displaystyle b } が与えられたときに最小二乗解 {\displaystyle \hat{x} } を求めるアルゴリズムを以下にまとめる。

アルゴリズム 12.1 QR分解を用いた最小二乗法の解法
各行が線形独立な {\displaystyle m\times n } 型行列 {\displaystyle A } と{\displaystyle m } 次元ベクトル {\displaystyle b } が与えられたとき、

  1. {\displaystyle A } のQR分解 {\displaystyle A=QR } を計算する。
  2. {\displaystyle Q^T b } を計算する。
  3. Back substitution により上三角行列の方程式 {\displaystyle R\hat{x}=Q^T b } を解く。


係数行列が正方行列の場合との比較

　係数行列 {\displaystyle A } が正方行列である連立一次方程式 {\displaystyle A } の解は {\displaystyle x=A^{-1}b } であった。{\displaystyle x } は {\displaystyle A } のQR分解を用いて {\displaystyle x=R^{-1}Q^T b } と表される ((11.4) 参照)。この方程式は (12.10) と形式的に同一である。(12.10) における唯一の相違点は、{\displaystyle A } と {\displaystyle Q } が正方行列である必要がなく、{\displaystyle R^{-1}Q^T b } が最小二乗解、すなわち、{\displaystyle Ax=b } を一般に満たさないということである。
　実際に、アルゴリズム 12.1 は連立方程式をQR分解で解くアルゴリズム 11.2 と同一である (11.2 の唯一の相違点は {\displaystyle A } と {\displaystyle Q } が縦長でもよいことである) 。
　{\displaystyle A } が正方行列であるとき、方程式 {\displaystyle Ax=b } を解くことと最小二乗法 {\displaystyle minimize\quad ||Ax-b||^2 } は同値であり、アルゴリズム 11.2 とアルゴリズム 12.1 も同値となる。したがってアルゴリズム 12.1 は {\displaystyle A } が正方行列であるとき方程式 {\displaystyle Ax=b } を解くアルゴリズム 11.2 の一般化であると考えることができ、{\displaystyle A } が縦長であるときの最小二乗解を計算できる。

バックスラッシュ記法

　行列を扱ういくつかのソフトウェアパッケージには優決定系の方程式の最小二乗法を計算するためのバックスラッシュ演算が提供されている (209ページ参照)。これらのパッケージでは {\displaystyle A\backslash b } は {\displaystyle A } が正則であるときの {\displaystyle Ax=b } の解 {\displaystyle x=A^{-1}b }、または {\displaystyle A } が縦長で線形独立な列をもつときの最小二乗解 {\displaystyle x=A^{\dagger}b } を意味する (ただしバックスラッシュ記法は数学的に正式な記法ではない)。

複雑性

　アルゴリズム 12.1 の最初のステップの複雑性は {\displaystyle 2mn^2 } フロップである。次のステップでは行列とベクトルの積を含むため、複雑性は {\displaystyle 2mn } である。三番目のステップでは {\displaystyle n^2 } フロップを要する。合計フロップ数は

{\displaystyle 2mn^2+2mn+n^2\approx 2mn^2 }
となる。第二項と第三項を無視すればそれは最初の {\displaystyle n } と {\displaystyle 2m } の要素よりも小さくなる。アルゴリズムのオーダーは {\displaystyle mn^2 } である。複雑性は {\displaystyle A } の行の次元に比例し変数の数の二乗に比例する。

スパース最小二乗法

　疎行列 {\displaystyle A } に対する最小二乗法はさまざまな場面に応用され、たとえば、12.1 の一般的なアルゴリズムで疎行列用に工夫したQR分解 (190ページ参照) を用いることでより効率的に解が求まる。
　{\displaystyle A } が疎行列であることを活用するもう一つの簡単な手法は正規方程式 {\displaystyle A^T Ax=A^T b } を巨大な疎行列を用いて解くことである。

{ \displaystyle \begin{bmatrix} 0&A^T\\ A&I \end{bmatrix} \begin{bmatrix} \hat{x}\\ \hat{y} \end{bmatrix} = \begin{bmatrix} 0\\b \end{bmatrix}. \quad (12.11) }
これは {\displaystyle m\times n } 型の方程式を組み合わせ、係数行列を正方行列にしたものである。{\displaystyle A } が疎行列ならばその係数行列も疎行列である。{\displaystyle (\hat{x},\hat{y}) } がこれらの方程式を満たすとき、(12.11) を満たす {\displaystyle \hat{x} } を求めることは容易であり、逆に {\displaystyle \hat{x} } が正規方程式を満たすとき、{\displaystyle (\hat{x},\hat{y}) } は {\displaystyle \hat{y}=b-A\hat{x} } とすると (12.11) を満たす。係数行列が疎行列である方程式を解くどのような手法でも (12.11) を解くことができる。

行列の最小二乗法

　最小二乗法のシンプルな拡張は方程式 {\displaystyle ||AX-B||^2 } を最小化する {\displaystyle n\times k } 型行列 {\displaystyle X } を求めることである。ここで {\displaystyle A } は {\displaystyle m\times n } 型行列、{\displaystyle B } は {\displaystyle m\times k } 型行列、ノルムは行列ノルムである。これは行列の最小二乗法 (matrix least squares problem) と呼ばれる。{\displaystyle k=1 } のとき {\displaystyle x } と {\displaystyle b } はベクトルであり、行列の最小二乗法はもとの最小二乗法となる。
　行列の最小二乗法は、実は {\displaystyle k } 個の最小二乗法の組に過ぎない。これを確かめるために、

{\displaystyle ||AX-B||^2=||Ax_1-b_1||^2+\cdots+||Ax_k-b_k||^2 }
と表す。{\displaystyle x_j } は {\displaystyle X } の {\displaystyle j } 番目の列、{\displaystyle b_j } は {\displaystyle B } の {\displaystyle j } 番目の列である (ここで、我々は行列のノルムの二乗は列のノルムの二乗和に等しいという性質を用いた)。よって、目的関数はそれぞれが {\displaystyle X } の列のみに依存する {\displaystyle k } 個の項の和である。したがって、我々は独立に列 {\displaystyle x_j } を選ぶことができ、それらはそれぞれに対応する項 {\displaystyle ||Ax_j-b_j||^2 } を最小化する。{\displaystyle A } が線形独立な列をもつとすると、解は {\displaystyle \hat{x}_j=A^{\dagger}b_j } である。したがって、行列の最小二乗法の解法は
{\displaystyle \begin{eqnarray} \hat{X} &=& \begin{bmatrix} \hat{x}_1 \cdots \hat{x}_k \end{bmatrix} \\ &=& \begin{bmatrix} A^{\dagger}b_1 \cdots A^{\dagger}b_k \end{bmatrix} \\ &=& A^{\dagger} \begin{bmatrix} b_1 \cdots b_k \end{bmatrix} \\ &=& A^{\dagger}B \quad\quad\quad\quad\quad\quad (12.12) \end{eqnarray} }
行列の最小二乗法の非常にシンプルな解 {\displaystyle \hat{X}=A^{\dagger}B } は、{\displaystyle k=1 } であるのもとの最小二乗法の場合にも成り立つ。多くの線形代数のソフトウェアパッケージはバックスラッシュ演算子 {\displaystyle A\backslash B } を {\displaystyle A^{\dagger}B } を表すために使っているが、これは数学的に正式な記法ではない。
　行列の最小二乗法はアルゴリズム 12.1 が factor-solve アルゴリズムのもう一つの例であるという事実を用いることで効率的に解くことができる。{\displaystyle \hat{X}=A^{\dagger}B } を解くために我々は {\displaystyle A } のQR分解を行い、アルゴリズム 12.1 のステップ2とステップ3を {\displaystyle k } 個の {\displaystyle B } の各列について実行する。合計の計算コストは {\displaystyle 2mn^2+k(2mn+n^2) } フロップである。 {\displaystyle k } が {\displaystyle n } と比較して小さいときは、およそ {\displaystyle 2mn^2} フロップであり、一つの最小二乗法 (すなわち、右側にベクトルがある場合) の計算コストと等しくなる。

*1:上三角行列を係数行列にもつ連立一次方程式を一番下の行から順に求める手法。