LU 分解に関して初歩的な内容から網羅的にまとめた.

ガウスの消去法

ガウスの消去法は, 前進消去による上三角行列の形成と後退代入の組み合わせのことをいい, その本質は行列の行基本変形, すなわち中学校で習う連立方程式の解法そのものである. 例えば, 何の芸もないが, 次の連立方程式

$$\begin{aligned} \begin{cases} x+2y&=&3 \\ 3x+4y&=&5 \end{cases}\leftrightarrow \left(\begin{array}{cc|c} 1 & 2 & 3 \\ 3 & 4 & 5 \end{array}\right)\tag{1} \end{aligned}$$

をガウスの消去法では次のようにして解くのであった.

$$\left(\begin{array}{cc|c} 1 & 2 & 3 \\ 3 & 4 & 5 \end{array}\right)\to \underbrace{ \left(\begin{array}{cc|c} \overbrace{1}^{\times (-3)} & \overbrace{2}^{\times (-3)} & \overbrace{3}^{\times (-3)} \\ \underline{3} & 4 & 5 \end{array}\right) }_{\rm 前進消去}\to \left(\begin{array}{cc|c} 1 & 2 & 3 \\ 0 & -2 & -4 \end{array}\right)\to \underbrace{ \left(\begin{array}{cc|c} 1 & \underline{2} & 3 \\ \overbrace{0}^{\times 1} & \overbrace{-2}^{\times 1} & \overbrace{-4}^{\times 1} \end{array}\right) }_{\rm 後退代入} \\ \to \left(\begin{array}{cc|c} 1 & 0 & -1 \\ 0 & -2 & -4 \end{array}\right) \\ \therefore (x,y)=(-1\div 1,-4\div (-2))=(-1,2)$$

この前進消去の操作は次のように行列の積で表現できる.

$$\left(\begin{array}{cc} 1 & 0 \\ -3 & 1 \end{array}\right) \left(\begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}\right)= \left(\begin{array}{cc} 1 & 2 \\ 0 & -2 \end{array}\right)$$

これを理解しておくと後述する LU 分解の理解に容易くなる. ここで, この時間計算量について, 一般の $n$ 次線形連立方程式 $X{\boldsymbol a}={\boldsymbol y}\ {\rm where}\ X\in\mathbb{R}^{n\times n}, {\boldsymbol a}\in\mathbb{R}^{n\times 1}, {\boldsymbol y}\in\mathbb{R}^{n\times 1}$ を用いて考えるする.

$$\left(\begin{array}{ccccc} x_{11} & x_{12} & x_{13} & \cdots & x_{1n} \\ x_{21} & x_{22} & x_{23} & \cdots & x_{2n} \\ x_{31} & x_{32} & x_{33} & \cdots & x_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & x_{n3} & \cdots & x_{nn} \end{array}\right) \left(\begin{array}{c}a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_n\end{array}\right)= \left(\begin{array}{c}y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n\end{array}\right)$$

まず前進消去の操作について考える. 前進消去では, $x_{11}$ を使って残りの $x_{21}, x_{31}, \cdots, x_{n1}$ を掃き出していくのであった. ただし $x_{21}, x_{31}, \cdots, x_{n1}$ は必ず $0$ となるのだから, これは実質の計算量として考える必要はないだろう. すると $x_{21}, x_{31}, \cdots, x_{n1}$ を掃きだすのに伴って, 実際の計算を要される部分というのは

$$\begin{array}{cccc} x_{22} & x_{23} & \cdots & x_{2n} \\ x_{32} & x_{33} & \cdots & x_{3n} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n2} & x_{n3} & \cdots & x_{nn} \end{array}$$

の箇所であることがいえる. この部分の成分数は $(n-1)^2$ なので, 一回目の前進消去では $(n-1)^2$ 回の計算を行うことがわかる. 次に, 二回目の前進消去では $x_{22}$ を使って, それと同列のそれよりも下の行の成分 $x_{32},\cdots,x_{n2}$ を先と同様に消していくわけだが, ここでも同様, $x_{32},\cdots,x_{n2}$ は必ず $0$ なので, 実際に計算を行う部分は

$$\begin{array}{ccc} x_{33} & \cdots & x_{3n} \\ \vdots & \ddots & \vdots \\ x_{n3} & \cdots & x_{nn} \end{array}$$

の箇所である. この部分の成分数は $(n-2)^2$ なので, 二回目の全身消去では $(n-2)^2$ 回の計算を行うことがわかる. これを繰り返すと, 計算を行う部分が残り一成分となるまで同様の計算回数がかかることがわかるから

$$(n-1)^2+(n-2)^2+\cdots+2^2+1^2 =\sum^{n}_{i=1}(i-1)^2$$

ここで, 高校数学の教科書にも載っている定理(証明略) $\sum^n_{j=1}j^2=\frac{n(n+1)(2n+1)}{6}$ を思い出せば, 先の式は $\sum^n_{i=1}(i-1)^2=\frac{n^3}{3}+n^2$ と表せることがわかる.

続いて, 後退代入について考える. 後退代入は, 上三角行列となっている係数行列に対して代入を繰り返し, 対角行列を形成していく操作であった.

$$\left(\begin{array}{ccccc} x_{11} & x_{12} & x_{13} & \cdots & x_{1n} \\ & x_{22} & x_{23} & \cdots & x_{2n} \\ & & x_{33} & \cdots & x_{3n} \\ & & & \ddots & \vdots \\ &&&& x_{nn} \end{array}\right) \left(\begin{array}{c}a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_n\end{array}\right)= \left(\begin{array}{c}y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n\end{array}\right)$$

この場合, まず $x_{nn}$ を使って, $x_{1n},x_{2n},x_{3n},\cdots,x_{(n-1)n}$ を消していくわけだが, この際, 係数行列に対する操作というのはただ $0$ にしていくということだけである. これを $x_{22}$ まで繰り返し行うわけだが, その間の係数行列に関する操作はただ $0$ にしていくということだけなので, これを計算量に含む必要はない. 実際に計算が発生するのは, 右辺ベクトルの部分である. まず $y_n$ に対する計算は, 後退代入の操作を考えれば, 当然なにもする必要はない. 次に $y_{n-1}$ に対する計算は $(y_{n-1})'=y_{n-1}-\frac{x_{(n-1)n}}{x_{nn}}y_n$ であり, 減算, 除算, 積算が各一回ずつ行われることがいえる. 次の $y_{n-2}$ に対する計算も同様 $y_{n-2}-\frac{x_{(n-2)n}}{x_{(n-1)n}}(y_{n-1})'$ であり, 計算の形式は先と全く同じであるが, その前の後退代入の結果 $(y_{n-1})'$ を用いている点で計算量は異なる. つまり, いま計算したい右辺ベクトルの成分を計算するのには, その一つ下の右辺ベクトルの成分に対する計算結果が必要となる(後退代入の操作そのもの)ことから, これが $y_1$ にまで及ぶことを考えると, 後退代入の総回数は

$$1+2+\cdots+(n-1)=\frac{n(n-1)}{2}$$

とかける(右辺への式変形の証明は高校数学の教科書で扱われているので略). ガウスの消去法は前進消去と後退代入の組み合わせであるので, その時間計算量は,

$$O(\frac{n^3}{3}+\frac{n(n-1)}{2})\tag{2}$$

ただし各項の中で最高次の係数だけを考えればよいので, ガウスの消去法の時間計算量は

$$\frac{1}{3}O(n^3)$$

である.

ガウス・ジョルダン法

ガウスの名がつく行列を用いた線形方程式の直接解法には, 今述べたガウスの消去法のほかにガウス・ジョルダン法というものがある. これらは明確に区別されたところであまり意味がないような気もするが, ガウスの消去法が上記のように係数行列となる部分を単位行列でない対角行列へと変形していったのに対し, ガウス・ジョルダン法はそれを直接単位行列となるように変形していく点で異なる. ガウス・ジョルダン法で同様にして計算していくと,

$$\left(\begin{array}{cc|c} 1 & 2 & 3 \\ 3 & 4 & 5 \end{array}\right)\to \left(\begin{array}{cc|c} \overbrace{1}^{\times (-3)} & \overbrace{2}^{\times (-3)} & \overbrace{3}^{\times (-3)} \\ \underline{3} & 4 & 5 \end{array}\right)\to \left(\begin{array}{cc|c} 1 & 2 & 3 \\ 0 & -2 & -4 \end{array}\right)\to \left(\begin{array}{cc|c} 1 & 2 & 3 \\ \overbrace{0}^{\times (-\frac{1}{2})} & \overbrace{-2}^{\times (-\frac{1}{2})} & \overbrace{-4}^{\times (-\frac{1}{2})} \end{array}\right) \\ \to \left(\begin{array}{cc|c} 1 & 2 & 3 \\ 0 & 1 & 2 \end{array}\right)\to \left(\begin{array}{cc|c} 1 & \underline{2} & 3 \\ \overbrace{0}^{\times (-2)} & \overbrace{1}^{\times (-2)} & \overbrace{2}^{\times (-2)} \end{array}\right)\to \left(\begin{array}{cc|c} 1 & 0 & -1 \\ 0 & 1 & 2 \end{array}\right) \\ \therefore (x,y)=(-1,2)$$

これは掃き出し法とも言われることがある. ガウス・ジョルダン法の計算量は $\frac{1}{2}O(n^3)$ なので(導出は省略), ガウス・ジョルダン法をわざわざ用いるシーンはあまりない.

ピボッティング

ところで, いま示したガウスの消去法, ガウスジョルダン法の手順は, このままでは困る場合がある. 例えば, 次の線形方程式をガウスの消去法で解いてみると

$$\left(\begin{array}{cccc|c} 1 & 2 & 7 & 6 & 6 \\ 2 & 4 & 4 & 2 & 2 \\ 1 & 8 & 5 & 2 & 12 \\ 2 & 4 & 3 & 3 & 5 \end{array}\right)\to \left(\begin{array}{cccc|c} \overbrace{1}^{\times (-2)} & \overbrace{2}^{\times (-2)} & \overbrace{7}^{\times (-2)} & \overbrace{6}^{\times (-2)} & \overbrace{6}^{\times (-2)} \\ \underline{2} & 4 & 4 & 2 & 2 \\ 1 & 8 & 5 & 2 & 12 \\ 2 & 4 & 3 & 3 & 5 \end{array}\right)\to \left(\begin{array}{cccc|c} 1 & 2 & 7 & 6 & 6 \\ 0 & 0 & -10 & -10 & -10 \\ 1 & 8 & 5 & 2 & 12 \\ 2 & 4 & 3 & 3 & 5 \end{array}\right) \\ \to \left(\begin{array}{cccc|c} \overbrace{1}^{\times (-1)} & \overbrace{2}^{\times (-1)} & \overbrace{7}^{\times (-1)} & \overbrace{6}^{\times (-1)} & \overbrace{6}^{\times (-1)} \\ 0 & 0 & -10 & -10 & -10 \\ \underline{1} & 8 & 5 & 2 & 12 \\ 2 & 4 & 3 & 3 & 5 \end{array}\right)\to \left(\begin{array}{cccc|c} 1 & 2 & 7 & 6 & 6 \\ 0 & 0 & -10 & -10 & -10 \\ 0 & 6 & -2 & -4 & 6 \\ 2 & 4 & 3 & 3 & 5 \end{array}\right)\to \left(\begin{array}{cccc|c} \overbrace{1}^{\times (-2)} & \overbrace{2}^{\times (-2)} & \overbrace{7}^{\times (-2)} & \overbrace{6}^{\times (-2)} & \overbrace{6}^{\times (-2)} \\ 0 & 0 & -10 & -10 & -10 \\ 0 & 6 & -2 & -4 & 6 \\ \underline{2} & 4 & 3 & 3 & 5 \end{array}\right) \\ \to \left(\begin{array}{cccc|c} 1 & 2 & 7 & 6 & 6 \\ 0 & 0 & -10 & -10 & -10 \\ 0 & 6 & -2 & -4 & 6 \\ 0 & 0 & -11 & -9 & -7 \end{array}\right)\to \left(\begin{array}{cccc|c} 1 & 2 & 7 & 6 & 6 \\ 0 & \overbrace{0}^{\times (-\frac{6}{0})} & \overbrace{-10}^{\times (-\frac{6}{0})} & \overbrace{-10}^{\times (-\frac{6}{0})} & \overbrace{-10}^{\times (-\frac{6}{0})} \\ 0 & \underline{6} & -2 & -4 & 6 \\ 0 & 0 & -11 & -9 & -7 \end{array}\right)$$

というように $0$ で割るということが起きてしまうのである. これを防ぐために考えられる方法としては, 掃きだすのに利用する値がその列の中で絶対値最大となるように行を入れ替える. これは, 部分ピボット選択付きガウスの消去法といわれる. ピボットとは, いま述べた掃きだすのに利用する値のことである. 部分ピボット選択といわれる理由は, 入れ替えの操作が行に対してのみ行われるからである. 列に対する入れ替え操作をも含んだ方法は完全ピボット選択といわれるが, 当然それは行基本変形の範疇でないので, そのまま計算を続行して単に右辺ベクトルを取り出せばよいという話ではなくなる. 完全ピボット選択は, 絶対値最大の値の選択の余地が部分ピボット選択よりも当然広がるので, 直感的に, より大きな絶対値の値をピボットとして選択できる確率が上がることは考えられるだろう. しかしながら, このアドバンテージは当然行列に依存したものであり, プログラムの複雑度が上がることに釣り合っていないため, 実際に用いられるようなことはあまりないと思われる.

LU 分解

漸く本題の LU 分解(LR 分解, 三角分解)について. 簡単のために式 $(1)$ を使って LU 分解の導出をする. 式 $(1)$ は次の式と同値である.

$$A^{(0)}\left(\begin{array}{c}x \\ y\end{array}\right)={\boldsymbol v} \ {\rm where}\ A^{(0)}=\left(\begin{array}{cc}1 & 2 \\ 3 & 4\end{array}\right),{\boldsymbol v}=\left(\begin{array}{c}3\\ 5\end{array}\right)$$

${A^{(0)}}^{-1}{\boldsymbol v}$ とすれば $(x,y)^T$ は求まるが, 逆行列の計算はガウスの消去法により $\frac{1}{3}O(n^3)$ の時間計算量がかかる. 一定の条件下でそれよりも高速に求める方法を考えることとする. いま $A^{(0)}$ を徐に上三角行列にすることを考えると, ガウスの消去法の前進消去より

$$A^{(1)}=L^{(1)}A^{(0)}=\left(\begin{array}{cc}1 & 2 \\ 0 & -2\end{array}\right) \leftrightarrow A^{(0)}={L^{(1)}}^{-1}A^{(1)} \ {\rm where}\ L^{(1)}=\left(\begin{array}{cc}1 & 0 \\ -3 & 1\end{array}\right)$$

従って

$${L^{(1)}}^{-1}A^{(1)}\left(\begin{array}{c}x \\ y\end{array}\right)={\boldsymbol v}$$

ここで, ${\boldsymbol b}=A^{(1)}(x,y)^T$ とおくと, 上の式は ${L^{(1)}}^{-1}{\boldsymbol b}={\boldsymbol v}$ と同値であり, この式を用いて ${\boldsymbol b}$ について解くことができる. まずこの時間計算量を考えるとする. $L^{(1)}$ は元々前進消去のための行列であり, それは必ず下三角行列である. 正則な下三角行列の逆行列は下三角行列であり(証明略), いま $L^{(1)}$ が正則であるとする(これが特異となるような場合には後述する PLU 分解が有効)と, その計算は前進代入(上記後退代入の下三角行列バージョンと考えればよい)を実行すればよいので, 時間計算量は $\frac{1}{2}O(n^2)\ \because$ となる.

その後に $A^{(1)}(x,y)^T={\boldsymbol b}$ を $(x,y)^T$ について解くわけであるが, $A^{(1)}$ は上三角行列であるので, その計算にはガウスの消去法の後退代入を実行すれば良く, 従ってその時間計算量は $\frac{1}{2}O(n^2)\ \because$ である.

よって, この一連の操作における時間計算量は $\frac{1}{3}O(n^3)$ であり, 部分ピボットつきガウスの消去法を実行した場合と変わらない.

しかし, $L U$ を流用できる(つまり, 共通の $A^{(0)}$ に対し異なる右辺ベクトル ${\boldsymbol v}$ から成る連立方程式を解く)とすればどうだろう. この場合, やらなければならない計算は前進代入および後退代入のみなので, 全体の時間計算量は $\frac{1}{2}O(n^2)$ となり, 先よりも高速に解を得ることができる.

いまの説明では, 式 \(\) において $A^{(0)}$ を ${L^{(1)}}^{-1}$ と $A^{(1)}$ に分解したが, これを LU 分解($L={L^{(1)}}^{-1},U=A^{(1)}$)といい, $L^{(i)}$ が正則ならば, 一般の場合においても同様にしていうことができる.

補足すると, $A^{(0)}$ に対し $n-1$ 回の前進消去をするというのは, $L^{(n-1)}\cdots L^{(2)}L^{(1)}A^{(0)}$ ということであり, $\lambda=L^{(n-1)}\cdots L^{(2)}L^{(1)}$ とおくと $\lambda A^{(0)}=U=A^{(n-1)}$ だから $A^{(0)}=\lambda^{-1}U$. ここで逆行列の公式 $(X Y)^{-1}=Y^{-1} X^{-1}$ より(証明略)上式となる. 一般論を得たところで, 実際に一つ LU 分解を実践してみることとする.

$$A^{(0)}= \left(\begin{array}{ccc} 3 & 1 & 0 \\ 6 & 1 & -2 \\ -3 & 0 & 3 \end{array}\right)\to \underbrace{\left(\begin{array}{ccc} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 1 & 0 & 1 \end{array}\right)}_{L^{(1)}} \left(\begin{array}{ccc} 3 & 1 & 0 \\ 6 & 1 & -2 \\ -3 & 0 & 3 \end{array}\right)= \underbrace{ \left(\begin{array}{ccc} 3 & 1 & 0 \\ 0 & -1 & -2 \\ 0 & 1 & 3 \end{array}\right) }_{A^{(1)}} \\ \to \underbrace{ \left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \end{array}\right) }_{L^{(2)}} A^{(1)}= \underbrace{ \left(\begin{array}{ccc} 3 & 1 & 0 \\ 0 & -1 & -2 \\ 0 & 0 & 1 \end{array}\right)}_{A^{(2)}}$$

より $$L=(L^{(2)}L^{(1)})^{-1}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{array}\right),U=A^{(2)}$$ 実際には, すべてを減算で考えることで, $(L^{(2)}L^{(1)})^{-1}$ の計算は楽に済む(つまり $A^{0}=({\boldsymbol a_1}^T,{\boldsymbol a_2}^T,{\boldsymbol a_3}^T)^T$ としたとき $\underline{2}{\boldsymbol a_1}-{\boldsymbol a_2}, \underline{-1}{\boldsymbol a_1}-{\boldsymbol a_3}, \underline{-1}{\boldsymbol a_2}-{\boldsymbol a_3}$ より $L$ が導けるということ).

この導出過程を見ればなんとなく LU 分解が一意となることは直感的にも納得できるが, 一応証明を与えておく.

次に, 次のような行列に対する LU 分解を考えてみる.

$$\begin{aligned} A=\left(\begin{array}{ccc} 0 & 1 & 0 \\ -8 & 8 & 1 \\ 2 & -2 & 0 \end{array}\right) \end{aligned}$$

見てわかるように, 前進消去の段階で $-8\div 0$ となってしまい計算できない. しかし, これも部分ピボット選択付きガウスの消去法と同様に, 絶対値最大の値がピボットとなるように行を予め入れ替えておけば, 計算が続行できる. そのような手続きのある LU 分解は PLU 分解といわれ, 置換行列 $P\in\mathbb{R}^{m\times n}$ をつかって, $A=P L U$ とする. 以下 $A$ をつかって導出してみることとする. ${\boldsymbol a_1}^T$ と ${\boldsymbol a_2}^T$ を入れ替えれば良いので,

$$\begin{aligned} \underbrace{\left(\begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right)}_{P^{(1)}}A= \underbrace{\left(\begin{array}{ccc} -8 & 8 & 1 \\ 0 & 1 & 0 \\ 2 & -2 & 0 \end{array}\right)}_{A'}\to \underbrace{\left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \frac{1}{4} & 0 & 1 \end{array}\right)}_{L^{(1)}}A'= \underbrace{\left(\begin{array}{ccc} -8 & 8 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & \frac{1}{4} \end{array}\right)}_{A'^{(1)}} \end{aligned}$$

より

$$L = {L^{(1)}}^{-1}=\left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ -\frac{1}{4} & 0 & 1\end{array}\right),U={A'}^{(1)}$$

とおくと $A'={L^{(1)}}^{-1}{A'}^{(1)}$ だから

$$P^{(1)}A=L U$$

$P^{(1)}$ は元々置換行列であるから正則であり, また直行行列でもある. すなわち ${P=P^{(1)}}^{-1}={P^{(1)}}^T$ とおけて

$$\begin{array}{rlr}{\displaystyle A}& {\displaystyle =}& {\displaystyle PLU}\\ {\displaystyle \left(\begin{array}{ccc}0& 1& 0\\ -8& 8& 1\\ 2& -2& 0\end{array}\right)}& {\displaystyle =}& {\displaystyle (\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& \end{array}}\end{array}$$