LU 分解

この記事は, 旧ブログから移植された記事です. よって, その内容として, 旧ブログに依存した文脈が含まれている可能性があります. 予めご了承下さい.

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

ガウスの消去法

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

{x+2y=33x+4y=5(123345)(1) \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}

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

(123345)(1×(3)2×(3)3×(3)345)前進消去(123024)(1230×12×14×1)後退代入(101024)(x,y)=(1÷1,4÷(2))=(1,2) \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)

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

(1031)(1234)=(1202) \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 分解の理解に容易くなる. ここで, この時間計算量について, 一般の nn 次線形連立方程式 Xa=y where XRn×n,aRn×1,yRn×1X{\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} を用いて考えるする.

(x11x12x13x1nx21x22x23x2nx31x32x33x3nxn1xn2xn3xnn)(a1a2a3an)=(y1y2y3yn)\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)

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

x22x23x2nx32x33x3nxn2xn3xnn \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}

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

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

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

(n1)2+(n2)2++22+12=i=1n(i1)2(n-1)^2+(n-2)^2+\cdots+2^2+1^2 =\sum^{n}_{i=1}(i-1)^2

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

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

(x11x12x13x1nx22x23x2nx33x3nxnn)(a1a2a3an)=(y1y2y3yn)\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)

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

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

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

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

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

13O(n3)\frac{1}{3}O(n^3)

である.

ガウス・ジョルダン法

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

(123345)(1×(3)2×(3)3×(3)345)(123024)(1230×(12)2×(12)4×(12))(123012)(1230×(2)1×(2)2×(2))(101012)(x,y)=(1,2) \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)

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

ピボッティング

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

(127662442218521224335)(1×(2)2×(2)7×(2)6×(2)6×(2)2442218521224335)(127660010101018521224335)(1×(1)2×(1)7×(1)6×(1)6×(1)0010101018521224335)(12766001010100624624335)(1×(2)2×(2)7×(2)6×(2)6×(2)001010100624624335)(127660010101006246001197)(1276600×(60)10×(60)10×(60)10×(60)06246001197) \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)

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

というわけで, 下記に部分ピボット選択付きガウスの消去法による計算過程を省略することなく書き下したが, とくに面白みもない上に長いので隠しておいた. クリックで展開.

LU 分解

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

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

A(1)=L(1)A(0)=(1202)A(0)=L(1)1A(1) where L(1)=(1031) 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)1A(1)(xy)=v {L^{(1)}}^{-1}A^{(1)}\left(\begin{array}{c}x \\ y\end{array}\right)={\boldsymbol v}

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

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

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

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

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

LU 分解 (外積形式ガウス法)

すべての前進消去の行列 L(i)L^{(i)} が正則ならば ARm×nA\in\mathbb{R}^{m\times n} に対する LU 分解は A=A(0)=LU where L=i=1n1L(i)1,U=A(n1)A=A^{(0)}=L U\ {\rm where}\ L=\prod_{i=1}^{n-1}{L^{(i)}}^{-1}, U=A^{(n-1)}

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

A(0)=(310612303)(100210101)L(1)(310612303)=(310012013)A(1)(100010011)L(2)A(1)=(310012001)A(2) 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=(100210111),U=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(L^{(2)}L^{(1)})^{-1} の計算は楽に済む(つまり A0=(a1T,a2T,a3T)TA^{0}=({\boldsymbol a_1}^T,{\boldsymbol a_2}^T,{\boldsymbol a_3}^T)^T としたとき 2a1a2,1a1a3,1a2a3\underline{2}{\boldsymbol a_1}-{\boldsymbol a_2}, \underline{-1}{\boldsymbol a_1}-{\boldsymbol a_3}, \underline{-1}{\boldsymbol a_2}-{\boldsymbol a_3} より LL が導けるということ).

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

定理 1

LU 分解された L,UL,U は一意に決まる

定理 1

正則行列 AA の LU 分解 A=LUA= L U に対して A=L1U1=L2U2L21L1=U2U11A= L_{1}U_{1} = L_{2}U_{2} \leftrightarrow L_{2}^{-1}L_{1} = U_{2}U_{1}^{-1} とおく. LL は元々下三角行列であり下三角行列の逆行列は下三角行列, また下三角行列の積は下三角行列だから L21L1L_{2}^{-1}L_{1} は下三角行列である. UU は元々上三角行列であり上三角行列の逆行列は上三角行列, また上三角行列の積は上三角行列だから U2U11U_{2}U_{1}^{-1} は上三角行列である. 従って, 両行列は上および下三角行列でなければならず, それを満たす唯一の行列は対角行列であり, LL の対角成分は元々 11 であるから L21L1=I=U2U11L_{2}^{-1}L_{1} = I = U_{2}U_{1}^{-1}. 故に L2=L1,U2=U1L_{2}=L_{1}, U_{2}=U_{1}.

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

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

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

(010100001)P(1)A=(881010220)A(1000101401)L(1)A=(8810100014)A(1) \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=(1000101401),U=A(1)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)1A(1)A'={L^{(1)}}^{-1}{A'}^{(1)} だから

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

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

A=PLU(010881220)=(010100001)(1000101401)(8810100014) \begin{aligned} A&=&P L U\\ \left(\begin{array}{ccc} 0 & 1 & 0 \\ -8 & 8 & 1 \\ 2 & -2 & 0 \end{array}\right)&=& \left(\begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right) \left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -\frac{1}{4} & 0 & 1 \end{array}\right) \left(\begin{array}{ccc} -8 & 8 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & \frac{1}{4} \end{array}\right) \end{aligned}

ここで aRn×1,vRn×1{\boldsymbol a}\in\mathbb{R}^{n\times 1},{\boldsymbol v}\in\mathbb{R}^{n\times 1} に対して PLUa=vP L U {\boldsymbol a}={\boldsymbol v} というように, PLU 分解を用いて連立方程式解くことを考えると, L(Ua)b=P1vL\underbrace{(U{\boldsymbol a})}_{{\boldsymbol b}}=P^{-1}{\boldsymbol v} だから先と同様にまず前進代入によって b{\boldsymbol b} を求め(このときの P1vP^{-1}{\boldsymbol v}P1=P(1)P^{-1}=P^{(1)} であり, また置換行列であるので, その計算は単に v{\boldsymbol v} を並び替えるだけである), Ua=bU{\boldsymbol a}={\boldsymbol b}a{\boldsymbol a} について後退代入によって求めればよい.

さて, 数学的な言葉では, 以上のように書き下すことで十分であるが, これをプログラムで組むことを考えると, 様々な工夫やアプローチが考えられる. まず L,UL, U はそれぞれ必ず上三角行列, 下三角行列であり 0011 の部分をそのまま持っているのは無駄である. 従って L,UL, U の値は次のように一つの行列として持っておけば十分.

(88101014014) \left(\begin{array}{ccc} -8 & 8 & 1 \\ 0 & 1 & 0 \\ -\frac{1}{4} & 0 & \frac{1}{4} \end{array}\right)

また, 置換行列 PP はただの並び替えでなので, これも各インデックスへの対応関係をテーブルにでもしておけば十分. サンプル実装, 実行例を下記に示す. Haskell で書いたわけだが, 普通の ST モナドによる実装とリストによる実装をそれぞれ行った. クリックで実装を開く.
λ> :m Data.Array Math.Matrix.LU
λ> lu [[0,1,0],[-8,8,1],[2,-2,0]] :: Maybe (PLU Array Int Double)
Just (array (0,2) [(0,1),(1,0),(2,2)],
{       -8.0    8.0     1.0     }
{       -0.0    1.0     0.0     }
{       -0.25   0.0     0.25    }
)
λ> luST' [[0,1,0],[-8,8,1],[2,-2,0]] :: Maybe (PLU Array Int Double)
Just (array (0,2) [(0,1),(1,0),(2,2)],
{       -8.0    8.0     1.0     }
{       -0.0    1.0     0.0     }
{       -0.25   0.0     0.25    }
)
連立方程式を解く. クリックで実装を開く.
λ> resolveLinearEq' [[1,2,7,6],[2,4,4,2],[1,8,5,2],[2,4,3,3]] [6,2,12,5] :: Maybe (Array Int Rational)
Just (array (0,3) [(0,(-3) % 1),(1,2 % 1),(2,(-1) % 1),(3,2 % 1)])
λ> plu = luST' [[1,2,7,6],[2,4,4,2],[1,8,5,2],[2,4,3,3]] :: Maybe (PLU Array Int Double)
λ> plu
Just (array (0,3) [(0,1),(1,2),(2,0),(3,3)],
{       2.0     4.0     4.0     2.0     }
{       0.5     6.0     3.0     1.0     }
{       0.5     0.0     5.0     5.0     }
{       1.0     0.0     -0.2    2.0     }
)
λ> fromJust plu `assign` (listArray (0, 3) [6,2,12,5]) :: Maybe (Array Int Double)
Just (array (0,3) [(0,-3.0),(1,2.0),(2,-1.0),(3,2.0)])
λ> fromJust plu `assign` (listArray (0, 3) [1,2,3,4]) :: Maybe (Array Int Double)
Just (array (0,3) [(0,0.6666666666666667),(1,0.6666666666666666),(2,-1.0),(3,1.0)])
λ> fromJust plu `assign` (listArray (0, 3) [5,6,7,8]) :: Maybe (Array Int Double)
Just (array (0,3) [(0,1.666666666666667),(1,0.8666666666666667),(2,-0.8),(3,1.2)])
ところで, LU 分解をしておくと逆行列も簡単に求めることがすぐに示せる. 逆行列とはそもそも AA1=IA A^{-1} =I であり, A1=(a11,a21,,a1m),I=(I1,I2,,Im)A^{-1}=({\boldsymbol a^{-1}_1},{\boldsymbol a^{-1}_2},\cdots,{\boldsymbol a^{-1}}_m), I=({\boldsymbol I_1},{\boldsymbol I_2},\cdots,{\boldsymbol I_m}) とすると行列の積の定義より Aai1=Ii where iZ+,1imA {\boldsymbol a^{-1}_i}={\boldsymbol I_i}\ {\rm where\ } i\in\mathbb{Z}^{+}, 1\leq i\leq m だから, A=LUA=L U と分解して 11 から mm までのすべての a1i{\boldsymbol a^{-1}}_i を得てそれらをそのまま 1 つの行列とすればよい. 当然ながら, 構成される方程式のうち変わる部分は a1i{\boldsymbol a^{-1}}_iIiI_i の部分だけなので, LU 分解は一度行うだけで済む. プログラムでの実行例. クリックで実装を開く.
λ> inverse' [[3,1,1],[5,1,3],[2,0,1]] :: Maybe (Matrix Array Int Rational)
Just
{       1 % 2   (-1) % 2        1 % 1   }
{       1 % 2   1 % 2   (-2) % 1        }
{       (-1) % 1        1 % 1   (-1) % 1        }

ただし, 逆行列の計算には今述べたようにすべての I1i{\boldsymbol I^{-1}}_i に関して代入操作を行わなければならないので, Ax=vA{\boldsymbol x}={\boldsymbol v} といった方程式を解く目的で逆行列 A1A^{-1} を求めることはただの愚行である.

また, LU 分解は行列式の計算も簡単にする. A=LUA= L U ならば積の行列式は行列式の積(証明略)なので A=LU=LU\left|A\right|=\left|L U\right|=\left|L\right|\left|U\right| であるが, 上および下三角行列の行列式は対角成分の積(証明略)であるので L=1\left|L\right|=1 である. よって A=i=1nuii\left|A\right|=\prod_{i=1}^{n}{\boldsymbol u}_{ii} である. 置換行列 PP を考慮すれば, いま SS を LU 分解の過程で行の入れ替えを行った回数としたとき, A=PLU=(1)Si=1nuii\left|A\right|=\left|P\right|\left|L\right|\left|U\right|=(-1)^S \prod_{i=1}^{n}{\boldsymbol u}_{ii} となる. また, 後述する Crout 法では UU のすべての対角成分を 11 とするので, その場合 A=(1)Si=1nlii\left|A\right| = (-1)^S \prod_{i=1}^{n}{\boldsymbol l}_{ii} となる. クリックで実装を開く.
λ> determinant $ toMat $ listArray ((0,0),(2,2)) [3,1,1,5,1,3,2,0,1]
2.0
λ> determinant $ toMat $ listArray ((0,0),(3,3)) [1,2,7,6,2,4,4,2,1,8,5,2,2,4,3,3]
120.0

なお, LU 分解は LDU 分解ともいわれることがある. その場合, 上記の UU に含まれる対角成分を対角行列 DD に分離して A=LDUA = L D U とする.

LDU=(100l2110lm1lm21)diag(d1,,dm)(1u12u1n01u2n001) L D U= \left(\begin{array}{cccc} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{m1} & l_{m2} & \cdots & 1 \end{array}\right) \mathrm{diag}(d_1,\cdots,d_m) \left(\begin{array}{cccc} 1 & u_{12} & \cdots & u_{1n} \\ 0 & 1 & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{array}\right)

ところで, ARm×nA\in\mathbb{R}^{m\times n} が対称行列ならば, この LDU 分解は A=LDLTA=L D L^{T} と計算することができる.

定理 2

A=AT where ARm×nA=A^T\ {\rm where}\ A\in\mathbb{R}^{m\times n} ならば AA の LDU 分解は A=LDLTA=L D L^{T}

定理 2

MT=D1UM^T = D^{-1}U とし AA の LDU 分解を A=LDMTA= L D M^T とおくと A=AT=(LDMT)T=MDLT=LUA = A^T = (L D M^T)^T = M D L^T = L U ここで, 第二辺から第三辺への変形は, 積の転置は積の左右を入れ替えた転置の積なる公式を用いた. このとき M(DLT)M(D L^T)LUL U は LU 分解の 2 つの表現であるが, 定理 1 より LU 分解は一意であるから M=LM=L でなければならない(後述する Crout 法の LU 分解ならば M=LDM=L D でなければならない. 導かれる結論は同じ). 従って LDMT=LDLTL D M^T = L D L^T.

ここまで述べてきた LU 分解の方法は, 外積形式ガウス法といわれるものであるが, LU 分解の他の方法として内積形式ガウス法(以下 Doolittle 法), クラウト法がある. 簡単のために n=3n=3 として, 行列 XX を次のように下三角行列 LR3×3L\in\mathbb{R}^{3\times 3} と上三角行列 UR3×3U\in\mathbb{R}^{3\times 3} に分解することを考える. このとき X=LDUX = L' D U' に対して L=LDL = L' D とおいて X=LUX = L U' というように分解できることを過程して行列 L,UL,U' を導出することを Doolittle 法, また U=DUU=D U' とおいて X=LUX = L' U というように分解できることを過程して行列 LUL' U を導出することを Crout 法という.

X=LU(x11x12x13x21x22x23x31x32x33)=(100l2110l31l321)(u11u12u130u22u2300u33):=Doolittle 法X=LU(x11x12x13x21x22x23x31x32x33)=(l1100l21l220l31l32l33)(1u12u1301u23001):=Crout 法 \begin{aligned} X= L U'\leftrightarrow \left(\begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{array}\right)=\left(\begin{array}{ccc} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{array}\right)\left(\begin{array}{ccc} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{array}\right)&:=&{\rm Doolittle\ 法} \\ X= L' U\leftrightarrow \left(\begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{array}\right)=\left(\begin{array}{ccc} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{array}\right)\left(\begin{array}{ccc} 1 & u_{12} & u_{13} \\ 0 & 1 & u_{23} \\ 0 & 0 & 1 \end{array}\right)&:=&{\rm Crout\ 法} \end{aligned}

いま行列 XX を Doolittle 法により LU 分解できたならば LUL U' を単に計算して X=LUX=L U' は次のようにかけるはずである.

(u11u12u13l21u11l21u12+u22l21u13+u23l31u11l31u12+l32u22l31u13+l32u23+u33)=(100l2110l31l321)(u11u12u130u22u2300u33) \left(\begin{array}{ccc} u_{11} & u_{12} & u_{13} \\ l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} \\ l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33} \end{array}\right)= \left(\begin{array}{ccc} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{array}\right)\left(\begin{array}{ccc} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{array}\right)

よって, 行列 XX の成分で行列 L,UL, U' を次のように書き換えることができる.

(x11x12x13x21x22x23x31x32x33)=(100x21u1110x31u11l321)(u11u12u130u22u2300u33)x21=l21u11l21=x21u11,l31 についても同様=(100x21x1110x31x11l321)(x11x12x130u22u2300u33) u11=x11,u12=x12,u13=x13=(100x21x1110x31x11l321)(x11x12x130x22x21x11x12x23x21x11x1300u33)x22=l21u12+u22u22=x22l21u12,u23 についても同様=(100x21x1110x31x11x32x31x11x12x22x21x11x121)(x11x12x130x22x21x11x12x23x21x11x1300x33x31x11x13x32x31x11x12x22x21x11x12(x23x21x11x13))x32=l31u12+l32u22l32=x32l31u12u22,u33 についても同様 \begin{aligned} \left(\begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{array}\right)&=& \left(\begin{array}{ccc} 1 & 0 & 0 \\ \frac{x_{21}}{u_{11}} & 1 & 0 \\ \frac{x_{31}}{u_{11}} & l_{32} & 1 \end{array}\right) \left(\begin{array}{ccc} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{array}\right)\because \begin{array}{l} x_{21}=l_{21}u_{11}\leftrightarrow l_{21}=\frac{x_{21}}{u_{11}}, \\ l_{31}\ {\rm についても同様} \end{array} \\ &=& \left(\begin{array}{ccc} 1 & 0 & 0 \\ \frac{x_{21}}{x_{11}} & 1 & 0 \\ \frac{x_{31}}{x_{11}} & l_{32} & 1 \end{array}\right)\left(\begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{array}\right)\because\ u_{11}=x_{11},u_{12}=x_{12},u_{13}=x_{13} \\ &=& \left(\begin{array}{ccc} 1 & 0 & 0 \\ \frac{x_{21}}{x_{11}} & 1 & 0 \\ \frac{x_{31}}{x_{11}} & l_{32} & 1 \end{array}\right)\left(\begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ 0 & x_{22}-\frac{x_{21}}{x_{11}}x_{12} & x_{23}-\frac{x_{21}}{x_{11}}x_{13} \\ 0 & 0 & u_{33} \end{array}\right)\because \begin{array}{l} x_{22}=l_{21}u_{12}+u_{22} \leftrightarrow u_{22}=x_{22}-l_{21}u_{12},\\ u_{23}\ {\rm についても同様} \end{array} \\ &=& \left(\begin{array}{ccc} 1 & 0 & 0 \\ \frac{x_{21}}{x_{11}} & 1 & 0 \\ \frac{x_{31}}{x_{11}} & \frac{x_{32}-\frac{x_{31}}{x_{11}}x_{12}}{x_{22}-\frac{x_{21}}{x_{11}}x_{12}} & 1 \end{array}\right)\left(\begin{array}{ccc} x_{11} & x_{12} & x_{13} \\ 0 & x_{22}-\frac{x_{21}}{x_{11}}x_{12} & x_{23}-\frac{x_{21}}{x_{11}}x_{13} \\ 0 & 0 & x_{33}-\frac{x_{31}}{x_{11}}x_{13}-\frac{x_{32}-\frac{x_{31}}{x_{11}}x_{12}}{x_{22}-\frac{x_{21}}{x_{11}}x_{12}}(x_{23}-\frac{x_{21}}{x_{11}}x_{13}) \end{array}\right)\\ &\because& \begin{array}{l} x_{32}=l_{31}u_{12}+l_{32}u_{22}\leftrightarrow l_{32}=\frac{x_{32}-l_{31}u_{12}}{u_{22}}, \\ u_{33}\ {\rm についても同様} \end{array} \end{aligned}

これをみると, 一行目, 一列目, 二行目, 二列目 \cdots と展開していくことで, 芋づる式に L,UL,U' が決まっていくことがわかる. この作業を一般化すると, uiju_{ij} の導出およびそれによって得られた値で lijl_{ij} を導出する部分に分けることができる. それぞれをいま漸化式で書くと

Doolittle:={{u1k=x1kuik=x1kj=1i1lijujk, (i=2,3,,k)lik=(xikj=1k1lijujk)ukk, (i=k+1,k+2,,n) \begin{aligned} {\rm Doolittle 法} &:=& \begin{cases} \begin{cases} u_{1k}&=&x_{1k} \\ u_{ik}&=&x_{1k}-\sum^{i-1}_{j=1}l_{ij}u_{jk},\ (i=2,3,\cdots,k) \end{cases} \\ l_{ik}=\frac{(x_{ik}-\sum^{k-1}_{j=1}l_{ij}u_{jk})}{u_{kk}},\ (i=k+1,k+2,\cdots,n) \end{cases} \end{aligned}

ただし ukk=0u_{kk}=0 の場合は計算できないので, 実際にはピボッティングを要することになるわけであるが, UUkk 番目の行が LL の対応する列の前に計算されるという Doolittle 法の性質上, このままではどの行が kk 番目に来るのかを処理以前に知ることができない. この問題は likl_{ik} の分子を次のように計算することで自明に克服できる.

si=xikj=1k1lijujk, (i=k,,n)s_i = x_{ik}-\sum^{k-1}_{j=1}l_{ij}u_{jk},\ (i=k,\cdots,n)

これにより sis_i の最大値を求め, 対応する行を入れ替えて最大の要素を kk 行目に入れることができる. 交換後は ukk=sku_{kk}=s_{k} となるが, UUkk 番目の行の他の要素はそれ以前と同様に計算することができ, 対応する LL の要素は lik=siukkl_{ik}=\frac{s_i}{u_{kk}} と得られる.

Crout 法も同様に, X=LUX=L' U を次のように書けるはずなので,

(l11l11u12l11u13l21l21u12+l22l21u13+l22u23l31l31u12+l32l31u13+l32u23+l33)=(l1100l21l220l31l32l33)(1u12u1301u23001) \begin{aligned} \left(\begin{array}{ccc} l_{11} & l_{11}u_{12} & l_{11}u_{13} \\ l_{21} & l_{21}u_{12}+l_{22} & l_{21}u_{13}+l_{22}u_{23} \\ l_{31} & l_{31}u_{12}+l_{32} & l_{31}u_{13}+l_{32}u_{23}+l_{33} \end{array}\right)= \left(\begin{array}{ccc} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{array}\right) \left(\begin{array}{ccc} 1 & u_{12} & u_{13} \\ 0 & 1 & u_{23} \\ 0 & 0 & 1 \end{array}\right) \end{aligned}

従ってこの作業の一般形は結果的に

Crout:={lik=xikj=1k1lijujk, (i=k,k+1,,n)ukj=(xkji=1k1lkiuij)lkk, (j=k,k+1,,n) \begin{aligned} {\rm Crout 法} &:=& \begin{cases} l_{ik}&=&x_{ik}-\sum^{k-1}_{j=1}l_{ij}u_{jk},\ (i=k,k+1,\cdots,n)\\ u_{kj}&=&\frac{(x_{kj}-\sum^{k-1}_{i=1}l_{ki}u_{ij})}{l_{kk}},\ (j=k,k+1,\cdots,n) \end{cases} \end{aligned}

LLkk 番目の列の要素を計算した後に最大値を求め, LL の要素を含む最初の k1k-1 列に対応する行列の行を交換できるため, Crout 法はピボットを簡単に選択できる.

参考文献

  • Richard Hamming (1987) “Numerical Methods for Scientists and Engineers (Dover Books on Mathematics)” Dover Publications, ISBN 9780486652412

活動継続のためのご支援を募集しています