LU分解 (LU decomposition)


数学的定義

n×n 行列を考える.

行列積 C = A B : (行列 C の i 行 j 列目の要素を c_ij と書けば)

       n
c_ij = Σ a_ik b_kj
       k=1

行列・ベクトル積 b = A x : (ベクトル bの i 行目の要素を b_i と書けば)

      n
b_i = Σ a_ik x_k
      k=1        

連立一次方程式の解: 与えられた行列 A と ベクトル b に対して、 A x = b を満たす ベクトル x を求める。

LU分解 とは 行列 A を、下三角行列 L、上三角行列 U を 用いて A = LU と分解すること。
ただし、 L については l_ij = 0 (1 ≦ i < j ≦ n), Uについては u_ij = 0 (1 ≦ j < i ≦ n)である。

LU 分解を用いた連立一次方程式の解法: A = LU と分解した後、 L y = byについて解き、さらに、 U x = yxについて解く。 L, U とも三角行列なので簡単に解ける。 ( x = U^{-1} L^{-1} b)


LU分解の分類

Doolittle型(l_ii=1), クラウト(Crout)型(u_ii=1)

クラウト型で考えると A=LU の A, L, Uの各成分について

        j
a_ij =  Σ l_ik u_kj                    (i ≧ j)
        k=1
        i
a_ij =  Σ l_ik u_kj                    (i <  j)
        k=1
であるから、総和の最後の項を分けると、
        j-1
a_ij =  Σ l_ik u_kj + l_ij             (i ≧ j)
        k=1
        i-1
a_ij =  Σ l_ik u_kj + l_ii u_ij        (i <  j)
        k=1
となり、これから、関係式(☆):
               j-1
l_ij =  a_ij - Σ  l_ik u_kj            (i ≧ j)
               k=1
               i-1
u_ij = (a_ij - Σ  l_ik u_kj) / l_ii    (i <  j)
               k=1
が成り立つことがわかる。

配列の添字と連続アクセス

Fortran言語とC言語で異なる。

LU分解の3種類のアルゴリズム

密行列に対するアルゴリズムのオーダはO(n^3)である。 通常,行列AをLU分解した結果のLとUは,どちらかの対角成分は 1なので,行列AをLとUで上書きする形でLU分解の結果を得ることが多い. 関係式☆に基づくアルゴリズムなら どんなアルゴリズムでも(例えば,ブロック化や再帰アルゴリズムでも) 正しい解は得られるが、単純な3重ループの場合, 最外ループについて Dongarra氏 流に分類すると次のようになる:

それぞれのアルゴリズムを理解するには, 2次元配列(行列)を図に書いてみて, 各計算ステップで読み出す要素と書き込む要素の対応を考えるとよい. (常人が理解するには、図を書くべきである。)

例: right-looking (down-looking?) アルゴリズムの場合 n=8のとき:
k=0, 1 の場合の処理後,k = 2のとき:

     for(j=k+1;j<n;j++)
       a[k][j] /= a[k][k];
は a[][] において
........
........
..*!!!!!
........
........
........
........
........
(*から読み出し.各!からの読み出しとそこへの書き込み)

     for(i=k+1;i<n;i++)
       for(j=k+1;j<n;j++)
         a[i][j] -= a[i][k] * a[k][j];
は a[][] において
........
........
...*****
..*!!!!!
..*!!!!!
..*!!!!!
..*!!!!!
..*!!!!!
(*から読み出し.各!からの読み出しとそこへの書き込み)

left-looking法などの用語は、 Fortran言語で Doolittle型(l_ii = 1) を用いる場合に使われてきたものであることに注意。


実用上の考慮

実用上は: を考慮する必要があるが、本実験演習では無視することにする。 無視するためには入力となる行列Aがある意味で「よい」性質を持つ必要がある。
並列プログラミング, 先頭ページへ
Masahiro Yasugi