LU分解

寄稿:東條遼平

下三角行列と上三角行列との分解

次のようにある行列Aを下三角行列Lと上三角行列Uに分解することを考えます.

img173.png

これは例えば行列Aが3×3の行列であれば次のようになります.

img174.png

このように分解することのメリットは,連立1次方程式を解くときを考えると, 下三角行列と上三角行列に分解していれば後は簡単に解くことができるという事です.

img175.png

のような連立1次方程式を解くことを考えてみます.LU分解を行うと,

img176.png

とでき,ベクトルcを用いて,

img177.png

とおくと,

img178.png

となります.ここで例えば3次であれば,

img179.png

のようになり,

img181.png

を解くことと同じ事になります.するとベクトルcはすぐにわかりますので, 今度は,

img177.png

より,

img180.png

を解けばベクトルxが求められます.これは先程と同じで,

img182.png

のようになりますので,ベクトルxが簡単に求まります.

このように,行列を下三角行列と上三角行列に分解できれば連立方程式が簡単に解けます.

分解する方法

行列Aを下三角行列Lと上三角行列Uに分解すれば連立方程式が簡単に解けることはわかりましたが, どうやったらそのように分解できるのでしょうか.とりあえず,

img183.png

を計算してみる事にします.

img184.png

どうやら対角成分より下と,それ以外で場合分けができそうです. もちろん,3次のときにどうなったからといって,数学的にそれが成り立つ事が証明されるわけではありませんが, ある程度自明なところがあります.

一般的に行列のかけ算を考えると,

img185.png

となりますが,ここで対角成分よりも下側とそれ以外で分けると,

img186.png

であり,行列Lの対角成分が1であることを考えると,

img187.png

とできるため,行列L,行列Uの成分は,

img188.png

となります.

実装の流れ

以上で一般的にプログラムが出来そうですが,よく見ると下三角行列の上の部分は無駄ですし, 上三角行列の下の部分も無駄です. しかも行列Lの対角成分は1ですのでここも無くても 問題ありません.そのため行列Lと行列Uの2つの行列をまとめて,

img189.png

のように値を保持します.それでは3次のときの例で流れを見てみます.

lu1.png

今この行列は元の行列Aを表しています.

lu2.png

このような順で計算すると,

lu3.png

このようになり,この右下2×2の小行列を考えると先程の3×3のときと 同じ状態になっています.そこでここまでを1セットと考え,同じ様に

lu4.png

と計算すると,

lu5.png

となってLU分解が完成します.

問題点

上のような流れでLU分解が進んでいきますが,それぞれの1セットの最初に 行の交換をする必要がある場合があります.というのは上の流れを見てもわかるように, そのときの対角成分の左上の値でその列の残りの部分を割っています. そのため割る値がゼロであれば困りますし,ゼロでなくてもゼロに近ければ 誤差も大きくなります.そこで毎回1セットが始まるときに,例えば上の流れの 最初の場合であれば,3行全ての行の内,1列目の要素の絶対値が最大であるものと 交換し,2セット目であれば2行目以降の中から2列目の値の絶対値が最大であるものと 交換します.もちろん交換が行われれば,どことどこが交換されたかが わかるようにどこかにその情報を保持しておかねばなりません.

ソースコード

width = m->width;
height = m->height;

if((order = (int*)malloc(sizeof(int)*height)) == NULL){
  fprintf(stderr, "Allocation Error\n");
  exit(1);
}

for(i=0; i<height; i++) order[i] = i;

動的に整数の配列を作成し,orderにポインタを保持しておきます. その後0,1,2・・・と初期化していますが, これは行を交換したときにこの値も同じく交換しておくことで, 例えば2行目と3行目を交換すれば,0,2,1・・・となることになります.

d = i;
max = fabs(m->a[i*width + i]);

for(j=i+1; j<height; j++){
  if(max < fabs(m->a[j*width + i])){
    d = j;
    max = fabs(m->a[j*width + i]);
  }
}

if(d!=i){
  temp = order[i];
  order[i] = order[d];
  order[d] = temp;

  ExchangeRowvector(m, i, d);
}

if(!(m->a[i*width + i])) continue;

for文までは現在注目している列の中から絶対値が最大の行を探しています. そして最大の行が現在の行でなければ現在の行と最大の行を入れ換えます. ExchangeRowvectorはi行目とd行目を交換する関数です. 交換しても対角成分がゼロであればその回はLU分解が終わっているといえるので, continueします.

肝心のLU分解の心臓部ですが,以上のような部分がごちゃごちゃしているだけで, 上の肝じゃない部分を除くと次のようになります.

for(i=0; i<height-1; i++){
  for(j=i+1; j<height; j++){
    m->a[j*width + i] /= m->a[i*width + i];
    for(k=i+1; k<width; k++){
      m->a[j*width + k] -= m->a[j*width + i]*m->a[i*width + k];
    }
  }
}

長々と書いて来ましたが,実際のLU分解の心臓部はたったのこれだけです.

Valid XHTML 1.1! home > コンピュータ > プログラミング >
リロード   新規 編集 凍結 差分 添付 複製 改名   トップ 一覧 検索 最終更新 バックアップ   ヘルプ   最終更新のRSS
Modified by 物理のかぎプロジェクト PukiWiki 1.4.5_1 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-jiPowered by PHP 5.2.17HTML convert time to 0.279 sec.