LU Factorization

ＬＵ分解

Last Update:2009年06月03日

概要

行列を解くのに直接法の一種であるＬＵ分解を使うことが多い．ＬＵ分解とは係数行列を下三角行列と上三角行列の積で表す方法である．ＬＵ分解の後は三角行列であることを利用して，前進・後退代入を使って高速に連立一次方程式を解くことができる．ここではＬＵ分解の方法について簡単に述べる．ＬＵ分解を説明するにあたって，ブロックＬＵ分解をまず説明する．ＬＵ分解はブロックＬＵ分解の特殊な場合を何度も繰り返すことで得られる．

ブロックＬＵ分解

次のような行列を考えよう．

$\begin{array}{l}\b{A} & \b{B}\\ \b{C} & \b{E}\end{array}$

この行列は次のように分解できることが簡単にわかる．

$\begin{array}{l}\b{A} & \b{B}\\ \b{C} & \b{E}\end{array}\] = \[\begin{array}{l}\b{I} & 0\\ \b{C}\b{A}^{-1} & I\end{array}\]\[\begin{array}{l}\b{A} & \b{B}\\ 0 & \b{E} - \b{C}\b{A}^{-1}\b{B}\end{array}$

これはブロックLU分解の一種(下行列の対角ブロックが単位行列となるもの)で，呼ぶ．ブロック上三角行列の右下に現れる項

$\b{E} - \b{C}\b{A}^{-1}\b{B} = \b{S}$

をシュール補元(Shur Complement)と呼び， $\b{S}$ と書くこともある．

ブロック下三角行列の逆行列

$\begin{array}{l}\b{I} & 0\\ \b{C}\b{A}^{-1} & I\end{array}\]\[\begin{array}{l}\b{I} & 0\\ -\b{C}\b{A}^{-1} & \b{I}\end{array}\] = \[\begin{array}{l}\b{I} & 0\\ 0 & \b{I}\end{array}$

が成り立つことから，

$\begin{array}{l}\b{I} & 0\\ \b{C}\b{A}^{-1} & I\end{array}\]^{-1} = \[\begin{array}{l}\b{I} & 0\\ -\b{C}\b{A}^{-1} & \b{I}\end{array}$

であることが分かる．

連立一次方程式のブロックごとの計算による求解

次のような連立一次方程式を考えよう．

$\[\begin{array}{l}\b{A} & \b{B}\\ \b{C} & \b{E}\end{array}\]\{x_A\\ x_B\} = \{y_A\\ y_B\}$

係数行列についてブロックＬＵ分解をすると，次のようになる．但しShur Complementを $\b{S}$ と書いた．

$\[\begin{array}{l}\b{I} & 0\\ \b{C}\b{A}^{-1} & I\end{array}\]\[\begin{array}{l}\b{A} & \b{B}\\ 0 & \b{S}\end{array}\]\{x_A\\ x_B\} = \{y_A\\ y_B\}$

左から下三角行列の逆行列を掛けて，

$\[\begin{array}{l}\b{A} & \b{B}\\ 0 & \b{S}\end{array}\]\{x_A\\ x_B\} = \[\begin{array}{l}\b{I} & 0\\ -\b{C}\b{A}^{-1} & I\end{array}\]\{y_A\\ y_B\}$

よって，次のように解が求まる

$\b{x}_B = \b{S}^{-1}\{-\b{C}\b{A}\b{y}_A + \b{y}_B\}$

$\b{x}_A = \b{A}^{-1}\{\b{y}_A - \b{B} \b{x}_B\}$

逆をとっているのは， $\b{A}$ と $\b{S}$ である．の逆が簡単に取れるような状況（例えばが対角行列）ではこのようなブロック分解における解法は解く行列の次数を減らすことができて効率が良い．

ブロック三角行列の固有値

上のブロックＬＵ分解の結果である，次のブロック上三角行列の固有値を求めてみよう．

$\[\begin{array}{l}\b{A} & \b{B}\\ 0 & \b{S}\end{array}\]\phi = \lambda \phi$

特性方程式は次のとおり

$\det\[\begin{array}{l}\b{A}-\lambda\b{I} & \b{B}\\ 0 & \b{S}-\lambda\b{I}\end{array}\] = 0 \\ \Leftright \det(\b{A}-\lambda\b{I})\det(\b{S}-\lambda\b{I})=0$

よって，固有値 $\lambda$ は $\b{A}$ の固有値， $\b{S}$ の固有値に等しくなる．ブロック三角行列の固有値は対角に位置する行列の固有値に等しいことが分かる．また，対角行列が単位行列であるブロック三角行列の固有値は１であることがわかる．以上から，もとの行列の固有値は $\b{A}$ と， $\b{S}$ の固有値に等しくなる,Shur Complementは元の行列の性質を保持しているので非常に重要である．

ブロックLDU分解

LU分解のブロック上三角行列をさらに次のように分解してみよう．

$\begin{array}{l}\b{A} & \b{B}\\ 0 & \b{E} - \b{C}\b{A}^{-1}\b{B}\end{array}\] = \[\begin{array}{l}\b{A} & 0 \\ 0 & \b{E} - \b{C}\b{A}^{-1}\b{B}\end{array}\]\[\begin{array}{l}\b{I} & \b{A}^{-1}\b{B}\\ 0 & \b{I}\end{array}$

これを用いると次のように対称に書ける．

$\begin{array}{l}\b{A} & \b{B}\\ \b{C} & \b{E}\end{array}\] = \[\begin{array}{l}\b{I} & 0\\ \b{C}\b{A}^{-1} & I\end{array}\]\[\begin{array}{l}\b{A} & 0 \\ 0 & \b{S}\end{array}\]\[\begin{array}{l}\b{I} & \b{A}^{-1}\b{B}\\ 0 & \b{I}\end{array}$

LL^T分解のブロックでの計算

もとの行列が正定値対称行列である場合を考える．つまり，

$\b{A}^T = \b{A}$

$\b{C} = \b{B}^T$

$\b{E}^T = \b{E}$

もとの行列の固有値は $\b{A}$ と $\b{S}$ の固有値に等しいので，もとの行列が正定値であるということは $\b{A}$ も $\b{S}$ も正定値ということである．このとき次のような実行列 L_A , L_S による分解が可能である．

A=L_A L_A^T , 　　 S=L_S L_S^T

これを使うとＬＤＵ分解の対角ブロック行列は次のように分解できる．

$\begin{array}{l}\b{A} & 0 \\ 0 & \b{S}\end{array}\] = \[\begin{array}{l}\b{L}_A & 0 \\ 0 & \b{L}_S\end{array}\]\[\begin{array}{l}L_A^T & 0 \\ 0 & L_S^T\end{array}$

これを上に代入して，

$\[\begin{array}{l}\b{A} & \b{B}\\ \b{B}^T & \b{D}\end{array}\] = \[\begin{array}{l}\b{I} & 0\\ \b{B}^T\b{A}^{-1} & I\end{array}\]\[\begin{array}{l}\b{L}_A & 0 \\ 0 & \b{L}_S\end{array}\]\[\begin{array}{l}L_A^T & 0 \\ 0 & L_S^T\end{array}\]\[\begin{array}{l}\b{I} & \b{A}^{-1}\b{B}\\ 0 & \b{I}\end{array}\]\\ \;\;\;\;\;\;\; = \[\begin{array}{l}\b{L}_A & 0\\ \b{B}^T\b{L}_A^{-T} & \b{L}_S\end{array}\] \[\begin{array}{l}\b{L}_A^T & \b{L}_A^{-T}\b{B} \\ 0 & \b{L}_S^T\end{array}\]\\ \;\;\;\;\;\;\; = \[\begin{array}{l}\b{L}_A & 0\\ \b{B}^T\b{L}_A^{-T} & \b{L}_S\end{array}\] \[\begin{array}{l}\b{L}_A & 0\\ \b{B}^T\b{L}_A^{-T} & \b{L}_S\end{array}\]^T$

よって，もとの行列を下三角行列とその転置の形で表すことができた．

LU分解

前節ではブロックＬＵ分解について解説した．ここではそれを用いてＬＵ分解を説明する．ＬＵ分解を説明するためにまず，ＬＤＵ分解について解説する．ＬＤＵ分解は処理に対称性があるので，理解しやすい．ＬＤＵ分解した後にＤをＵにかけたり，ＤをＬにかけたりすることで，ＬＵ分解が得られる．ＬＤＵ分解はブロックＬＤＵ分解の特殊な場合を繰り返した物として考えることができるので，ここではそれを使って解説する．

LDU分解

ＬＤＵ分解は行列を下三角行列，対角行列，上三角行列の積に分けることである．ブロックＬＤＵ分解は,もとの行列を大きくブロックに分けて，ブロックごとの計算で全体の分解を行った．ブロックの分け方を代えて，サイズが１のブロックと，n-1のブロックに分けるようなブロックＬＤＵ分解を考える．これを繰り返すことでLDU分解が得られる．

$\b{S}_0=\[\begin{array}{l}a_0 & \b{b}^T_0\\ \b{c}_0 & \b{E}_0\end{array}\] = \[\begin{array}{l}1 & 0\\ \frac{1}{a_0}\b{c}_0 & \b{I}\end{array}\]\[\begin{array}{l}a_0 & 0\\ 0 & \b{E}_0 - \frac{1}{a}\b{c}_0\b{b}_0^T\end{array}\]\[\begin{array}{l}1 & \frac{1}{a_0}\b{b}_0^T\\ 0 & \b{I}\end{array}\] = \b{L}_1 \b{D}_1 \b{U}_1$

さて，

$\b{D}_1 = \[\begin{array}{l}a_0 & 0\\ 0 & \b{S}_1\end{array}\],\;\;\;\b{S}_1 = \b{E}_0 - \frac{1}{a_0}\b{c_0}\b{b_0}^T$

であった．この右下の(n-1)次の正方行列である $\b{S}_1 =\b{E}_0 - \frac{1}{a_0}\b{c_0}\b{b_0}^T$ が次のようにＬＤＵ分解されるとする．

$\b{S}_1=\[\begin{array}{l}a_1 & \b{b}^T_1\\ \b{c}_1 & \b{E}_1\end{array}\] = \[\begin{array}{l}1 & 0\\ \frac{1}{a_1}\b{c}_1 & \b{I}\end{array}\]\[\begin{array}{l}a_1 & 0\\ 0 & \b{E}_1 - \frac{1}{a_1}\b{c}_1\b{b}_1^T\end{array}\]\[\begin{array}{l}1 & \frac{1}{a_1}\b{b}_1^T\\ 0 & \b{I}\end{array}\] = \bar{\b{L}}_2 \bar{\b{D}}_2 \bar{\b{U}}_2$

これを使うと，もとの行列 $\b{D}_1$ は次のように書ける．

$\b{D}_1 = \[\begin{array}{l}a_0 & 0\\ 0 & \bar{\b{L}}_2\bar{\b{D}}_2\bar{\b{U}}_2\end{array}\] = \[\begin{array}{l}1 & 0\\ 0 & \bar{\b{L}}_2\end{array}\]\[\begin{array}{l}a_0 & 0\\ 0 & \bar{\b{D}}_2\end{array}\]\[\begin{array}{l}1 & 0\\ 0 & \bar{\b{U}}_2\end{array}\] = \tilde{\b{L}}_2\b{D}_2\tilde{\b{U}}_2$

これを代入すると，もとの行列 $\b{S}_0$ は次のように書ける

$\b{S}_0 = (\tilde{\b{L}}_1 \tilde{\b{L}}_2) \b{D}_2 (\tilde{\b{U}}_2 \tilde{\b{U}}_1)$

ここで， $\b{L}_2=\tilde{\b{L}}_1\tilde{\b{L}}_2$ について具体的に計算してみよう．

$\b{L}_2 = \tilde{\b{L}}_1\tilde{\b{L}}_2 = \[\begin{array}{l}1 & 0\\ \frac{1}{a_0}\b{c}_0 & \b{I}\end{array}\]\[\begin{array}{l}1 & 0\\ 0 & \bar{\b{L}}_2\end{array}\] = \[\begin{array}{l}1 & 0\\ \frac{1}{a_0}\b{c}_0 & \bar{\b{L}}_2\end{array}\] = \[\begin{array}{l}1 & 0\\ \frac{1}{a_0}\b{c}_0 & \[\begin{array}{l}1 & 0\\ \frac{1}{a_1}\b{c}_1 & \b{I}\end{array}\] \end{array}\]$

これは単位行列の前２列の対角より下が埋まった下三角行列になっていることがわかる．同様にして， $\b{D}_2$ $\b{U}_2$ について具体的に書くと次のとおり，

$\b{D}_2 = \[\begin{array}a_0 & 0 & 0\\ 0 & a_1 & 0 \\ 0 & 0 & \b{E}_1-\frac{1}{a_1}\b{c}_1\b{b}_1^T\end{array}\]$ 　　　　 $\b{U}_2 = \tilde{\b{U}}_2\tilde{\b{U}}_1 = \[\begin{array}{l}1 & \frac{1}{a_0}\b{b}^T_0\\ 0 & \[\begin{array}{l}1 & \frac{1}{a_1}\b{b}^T_1\\ 0 & \b{I}\end{array}\] \end{array}\]$

$\b{D}_2$ は対角ブロック行列で，始めの２つは対角行列で，以降はブロック行列となっている． $\b{U}_2$ については単位行列の前２行の対角より右がうまった上三角行列になっていることがわかる．

アルゴリズム

これらの作業をもとの行列の次数であるn回繰り返せば， $\b{D}_n$ が対角行列になる．n回繰り返した時にＬＤＵ分解が得られる．

もとの行列を S_0 とするとき，順番にこれらを求めるアルゴリズムは次のとおり

i=0とする
は(n-i)次の正方行列である．の１行１列目の成分をとする.の１行目で２列目以降の成分をベクトル表記して $\b{b}^T_i$ とする．同様にの１列で２行目以降の成分をベクトル表記して $\b{c}_i$ とする． $\b{b}_i,\b{c}_i$ は(n-i-1)次のベクトルである．またの２行目以降，２列目以降からなる(n-i-1)次の正方行列を $\b{E}_i$ とする． $S_{i+1} = \b{E}_i - \frac{1}{a_i}\b{c}_i\b{b}^T_i$ とする．
を $\b{D}$ のi番目の対角に代入． $\frac{1}{a_i}\b{c}_i$ を $\b{L}$ のi列目の対角の下に並べる． $\frac{1}{a_i}\b{b}^T_i$ を $\b{U}$ のi行目の対角から右に並べる
i<nならiに１を加えて２へ戻る．

これを図示すると下の図のようになる．元の行列 $\b{S}_0$ から始まり，右下へ向けて逐次的に計算を進めていく．

図１．ＬＤＵ分解のアルゴリズム

成分の計算

上のアルゴリズムを成分毎に書いてみよう．元の行列 $\b{A}$ として，ＬＤＵ分解した行列の成分は次のようになる．

但し， A,L,D,U のij成分をそれぞれ $a_{ij}$ のように書くとする．

$d_{ij} = \{\begin{array}{ll} a_{ii} - \sum_{k=0}^{i-1}l_{ik} u_{ki} d_{kk}\;\; & (i=j)\\ 0 & (i\ne j)\end{array}\;\;\;\; (i=0,1,\ldots,n-1)$

$l_{ij} = \{\begin{array}{ll}$a_{ij} - \sum_{k=0}^{j-1} l_{ik} u_{kj} d_{kk}$/d_{jj}\;\; & (i>j)\\ 1 &(i=j)\\ 0 & (i<j)\end{array}\;\;\;\; (i,j=0,1,\ldots,n-1)$

$u_{ij} = \{\begin{array}{ll}$a_{ij} - \sum_{k=0}^{i-1} l_{ik} u_{kj} d_{kk}$/d_{ii}\;\; & (i<j)\\ 1 & (i=j)\\ 0 & (i>j)\end{array}\;\;\;\; (i,j=0,1,\ldots,n-1)$

ＬＵ分解

と下三角行列と上三角行列の積になる．これをＬＵ分解という． $\b{A}$ が非対称な場合はクラウト分解，対称の場合はコレスキー分解と呼ばれる．

ＬＵ分解その１(Lの対角が１な分解)

前節でLDU分解を用いて分解した行列を $\b{A}=\b{L}\bar{D}\bar{U}$ とする．

$\b{U} = \bar{D}\bar{U}$

とすると， $\b{U}$ は上三角行列となるので，

$\b{A} = \b{L}\b{U}$

LDU分解を参考にして上のアルゴリズムを成分毎に書いてみよう．元の行列 $\b{A}=\b{L}\bar{\b{D}}\bar{\b{U}}$ として，ＬＤＵ分解した行列の成分は次のようになる． $\b{U} = \bar{\b{D}}\bar{\b{U}}$ から成分ごとに得られる次の関係式

$\bar{d}_{ii} = u_{ii}$ ，　　　 $u_{ij} = \bar{u}_{ij}\bar{d}_{ij}$

を代入して，次が得られる．

$l_{ij} = \{\begin{array}{ll}$a_{ij} - \sum_{k=0}^{j-1} l_{ik} u_{kj}$/u_{jj}\;\; & (i>j)\\ 1 &(i=j)\\ 0 & (i<j)\end{array}\;\;\;\; (i,j=0,1,\ldots,n-1)$

$u_{ij} = \{\begin{array}{ll}a_{ij} - \sum_{k=0}^{i-1} l_{ik} u_{kj}\;\; & (i\le j)\\ 0 & (i>j)\end{array}\;\;\;\; (i,j=0,1,\ldots,n-1)$

ＬＵ分解その２(Ｕの対角が１な分解)

前節でLDU分解を用いて分解した行列を $\b{A}=\bar{\b{L}}\bar{D}\b{U}$ とする．

$\b{L} = \bar{L}\bar{D}$

とすると， $\b{U}$ は上三角行列となるので，

$\b{A} = \b{L}\b{U}$

LDU分解を参考にして上のアルゴリズムを成分毎に書いてみよう．元の行列 $\b{A}=\bar{\b{L}}\bar{\b{D}}\b{U}$ として，ＬＤＵ分解した行列の成分は次のようになる． $\b{L} = \bar{\b{L}}\bar{\b{D}}$ から成分ごとに得られる次の関係式

$\bar{d}_{ii} = u_{ii}$ ，　　　 $u_{ij} = \bar{u}_{ij}\bar{d}_{ij}$

を代入して，次が得られる．

$l_{ij} = \{\begin{array}{ll}a_{ij} - \sum_{k=0}^{j-1} l_{ik} u_{kj}\;\; & (i\ge j)\\ 0 & (i<j)\end{array}\;\;\;\; (i,j=0,1,\ldots,n-1)$

$u_{ij} = \{\begin{array}{ll}$a_{ij} - \sum_{k=0}^{i-1} l_{ik} u_{kj}$/l_{ii}\;\; & (i< j)\\ 1 & (i=j) \\ 0 & (i>j)\end{array}\;\;\;\; (i,j=0,1,\ldots,n-1)$

これは行列のデータ構造として，ＣＲＳデータ形式を用いる場合に良く使われる．

連立一次方程式の求解

次のような連立一次方程式を解くとする．

$\b{A}\b{x} = \b{y}$

ここで係数が次のようにＬＵ分解されているとすると，

$\b{A} = \b{L}\b{U}$

これを解く際には，次のようにベクトル $\b{z}$ を導入して，まず下三角行列について解き，次に上三角行列について順番に解くことで得られる．下三角行列や，上三角行列について解くときには，後に述べるように前進代入，後退代入によって解が得られる．

$\b{z} = \b{L}^{-1}\b{y}$

$\b{x} = \b{U}^{-1}\b{z}$

Lz=yを解いてｚを求める方法

前進代入を用いてできる

$\b{L}\b{z} = \b{y}$

$l_{00}z_0 = y_0\;\;\; \Leftright\;\;\; z_0 = y_0/l_{00}$

$l_{10}z_0 + l_{11}z_1 = y_1\;\;\; \Leftright\;\;\; z_1 = \{y_1-l_{10}z_0\}/l_{11}$

$l_{20}z_0 + l_{21}z_1 + l_{22}z_2 = y_2\;\;\; \Leftright\;\;\; z_2 = \{y_2-l_{20}z_0-l_{21}z_1\}/l_{22}$

$\sum_{i=0}^k l_{ki}z_i = y_k\;\;\; \Leftright\;\;\; z_k = \{y_k-\sum_{i=0}^{k-1}l_{ki}z_i\}/l_{kk}\;\;\;\;(k=0,1,\ldots,n-1)$

Ux=zを解いてxを求める方法

後退消去を用いて求めることができる．

$\b{U}\b{x} = \b{z}$

xの最後の要素から，順番に最初の要素に向けて求めていく

$u_{(n-1,n-1)}x_{n-1} = z_{n-1}\;\;\; \Leftright\;\;\; x_{n-1} = z_{n-1}/u_{(n-1,n-1)}$

$u_{(n-2,n-2)}x_{n-2} + u_{(n-2,n-1)}x_{n-1} = z_{n-2}\;\;\; \Leftright\;\;\; x_{n-2} = \{z_{n-2}-u_{(n-2,n-1)}x_{n-1}\}/u_{(n-2,n-2)}$

$u_{(n-3,n-3)}x_{n-3} + u_{(n-3,n-2)}x_{n-2} + u_{(n-3,n-1)}x_{n-1} = z_{n-3}\\\;\;\;\;\;\;\;\; \Leftright\;\;\; x_{n-3} = \{z_{n-3}-u_{(n-3,n-1)}x_{n-1}-u_{(n-3,n-2)}x_{n-2}\}/u_{(n-3,n-3)}$

$\sum_{i=1}^k u_{(n-k,n-i)}x_{n-i} = z_{n-k}\;\;\; \Leftright\;\;\; x_{n-k} = \{\sum_{i=1}^k z_{n-k}-u_{(n-k,n-i)}x_{n-k}\}/u_{(n-k,n-k)}\;\;\;(k=1,2,\ldots,n)$

圧縮表示と計算の高速化

分解する前のもとの行列 $\b{A}$ に加え，LとＵを表現するためのメモリを別に確保することは効率が悪い．Ｌは下三角行列で，上半分は０であるし，Ｕは下三角行列で下半分は０である．もとの行列を変形して，上半分がＵになるようにして，下半分がＬになるようにすると，これらを１つの行列の中に表示することができる．

いまここで，Ｕの対角が１になるようにＬＵ分解を行う場合を考える．つまり $A=$\bar{L}\bar{D}$\b{U} = LU$ このとき $\b{A}$ を変形して次のような行列 $\b{A}'$ になるようにする．但し， $\b{A}',\;\b{L},\;\b{U}$ のij成分をそれぞれ $a'_{ij}$ ， $l_{ij}$ ， $u_{ij}$ とする．

$a'_{ij} = \{\begin{array}{ll}l_{ij} & (i>j)\\ 1/l_{ii} & (i=j) \\ u_{ij} & (i<j)\end{array}$

の対角成分が $\b{L}$ の対角成分の逆になっていることに注目されたい．これは $\b{A}'$ の対角成分はＬＵ分解や対角成分は前進代入の際に使われるが，常に逆の値として使われるので，逆の値を保存しておくことで高速化が期待できる．さて，Ｕの対角が１になるようなＬＵ分解は行ごとの分解に向いている．これはＣＲＳのような行ごとのデータ構造にとって，最も適した分解方法だと言える．このような行列 $\b{A}'$ を求めるアルゴリズムは次のとおりとなる．