5 LU分解

5.1 LU分解による解の計算方法

係数行列 $\boldsymbol{A}$ が下三角行列 $\boldsymbol{L}$ と上三角行列 $\boldsymbol{U}$ の積に展開でき たとする。

$$\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}$$ (14)

下三角行列と上三角行列の要素を書き出すと

$$\begin{bmatrix}\alpha_{11} & 0 & 0 & 0 \alpha_{21} & \alpha_{2... ...1} & a_{32} & a_{33} & a_{34} a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}$$ (15)

となる。

このようにLU分解できると、連立1次方程式は

$$\boldsymbol{Ax}=(\boldsymbol{LU})\boldsymbol{x}=\boldsymbol{L}(\boldsymbol{Ux})=\boldsymbol{b}$$ (16)

と書ける。これをさらに書き換えると、

$$\boldsymbol{Ly}=\boldsymbol{b}$$ (17)

$$\boldsymbol{Ux}=\boldsymbol{y}$$ (18)

となる。これらの連立方程式の解 $\boldsymbol{y}$ と $\boldsymbol{x}$ は、それぞれの係数 が三角行列なので容易に計算できる(ガウス消去法と後退代入の説明を見よ)。 $y$ の方は、係数が下三角行列なので$1$ 〜$N$ まで前進代入により解ける。 具体的には、

$$\begin{aligned}y_1&=\frac{1}{\alpha_{11}}b_1 y_i&=\frac{1}{\a... ...}^{i-1}\alpha_{ij} y_j\right) \qquad i=2,3,\cdots,N \end{aligned}$$

である。この$y$ が求まったならば、今度は係数が上三角行列の式()の $\boldsymbol{x}$ を求める。これは、$N$ 〜$1$ の順序で計算する後 退代入を使う。

$$\begin{aligned}x_N&=\frac{1}{\beta_{NN}}y_N x_i&=\frac{1}{\be... ...1}^N\beta_{ij} x_j\right) \qquad i=N-1,N-2,\cdots,1 \end{aligned}$$

これらの前進代入と後退代入は、コンピューターにとって非常に簡単に計算できる。これ は、係数行列$A$ をLU分解できれば、連立方程式は簡単に解けると言っている。次節でLU 分解の方法を詳しく説明する。

いったんLU分解が出来てしまえば、式(1)の右辺 $\boldsymbol{b}$ が変わっ ても、そのLU分解の形を変える必要がない。右辺が変わっても、LU分解は1回で済む。こ れが、ガウスの消去法と後退代入を組み合わせた方法やガウス・ジョルダン法に比べて、 際立って優れている点である。

5.2 LU分解(クラウトのアルゴリズム)

LU分解するということは、式(15)の $\alpha_{ij}$ と $\beta_{ij}$ を計算することにほかならない。この式の行列方程式は、$N^2+N$ の未知数と$N^2$ の式を 含む。未知数の数が方程式の数より多いので、$N$ 個の未知数を勝手に決めて残りを計算 することが可能である。従って、LU分解は一意に決まらない、計算しやすいように$N$ 個 の未知数を決めることができる。ここでは、LU分解のクラウト(Crout)のアルゴリズムに 従い、

$$\alpha_{ii}=1 \qquad i=1,2,3,\cdots,N$$ (21)

とする。

それでは、クラウトのアルゴリズムによるLU分解の手順を示すことにする。

1. $\alpha_{ii}=1\quad(i=1,\cdots,N)$ とします。この操作により、解 くべき行列方程式(15)は

   $$\begin{bmatrix}1 & 0 & 0 & 0 \alpha_{21} & 1 & 0 & 0 \alpha... ...1} & a_{32} & a_{33} & a_{34} a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}$$ (22)

   
   
   と変形できる。
2. この式を見ると、 $\beta_{ij}$ と $\alpha_{ij}$ が次に示す順序で簡単に求められ ることが分かる。まずは式を見て分かるように、 $\beta_{11}$ が直ちに計算できる。 次に $\beta_{11}$ を利用して、 $\alpha_{21},\alpha_{31},\alpha_{41}$ を求める ことができる。これで、$L$ と$U$ の第1列目が求められた。次に第2列目である。こ れも $\beta_{12}$ は直ちに計算できる。そうして、これまで分かっている $\beta_{ij}$ と $\alpha_{ij}$ を使うと、 $\beta_{22},\alpha_{32},\alpha_{42}$ を求めることができる。これで第2列目は終わりで、同じことを繰り返すと、全て の $\beta_{ij}$ と $\alpha_{ij}$ が計算できる。これをアルゴリズムにすると次の ようになる。
   $j=1,2,3,\cdots,N$ という順序で計算する。 $\boldsymbol{L}$ と $\boldsymbol{U}$ の$j$ 列目を計算することになる。具体的には、以下のようにして、$j$ 列目の $\beta_{ij}$ と $\alpha_{ij}$ を求める。
   • まず、 $i=1,2,\cdots,j$ について、次式に従い $\beta_{ij}$ を 計算する。

     $$\begin{aligned}\beta_{1j}&=a_{1j} \beta_{ij}&=a_{ij}-\sum_{k=1}^{i-1}\alpha_{ik}\beta_{kj} \end{aligned}$$

     
     

   • 次に、 $i=j+1,j+2,\cdots,N$ について、 $\alpha_{ij}$ を計算する。

     $$\begin{aligned}\alpha_{ij}&=\frac{1}{\beta_{jj}}\left( a_{ij}-\sum_{k=1}^{j-1}\alpha_{ik}\beta_{kj} \right) \end{aligned}$$

     
     

   • これで、 $\boldsymbol{L}$ と $\boldsymbol{U}$ の$j$ 列目が完成したので、同じ操作を $j+1$ 列目に行う。同じことを繰り返して、LU分解の行列を 完成させる。

この方法により、LU分解ができる。次に示すピボット選択をしなければ、アルゴリズムは 非常に単純である。

5.3 ピボット選択

ここでも、ピボット選択の問題が出てくる。式(24)の $\beta_{jj}$ で割る部分である。安定な解を求めるためには、ピボット選択は必要不可欠ということで ある。完全ピボット選択は複雑なので、部分ピボット選択で十分であろう。

ではどうするかですが、これも途中($j$ 列)まで分解した行列は崩したくない。そのため には、行列 $\boldsymbol{L}$ の行を交換し、それに対応した行列 $\boldsymbol{A}$ の行を交換すれば問題 がもっとも少なくなる。当然、行列 $\boldsymbol{U}$ は行も列も変化しない。最終的には行を交換 した行列 $\boldsymbol{A}^\prime$ のLU分解が出来る。連立1次方程式を解くときには、同様に $\boldsymbol{b}$ の行も交換しておく。ただし、行の交換であるため、解 $\boldsymbol{x}$ の要素の順序は 入れ替わらない。

つぎに、どのようにして交換する行を決めるかである。一般的には、 $\beta_{jj}$ が大き くなるように選択すれば良い結果が得られる。クラウト法のピボット選択は、次のように 進める。$j$ 列目のピボットを選択する場合についてである。

1. まずは、$j-1$ 列目までの行列 $\boldsymbol{L}$ の各行の要素の最大値を1に規 格化する。同時に、対応する行列 $\boldsymbol{A}$ の行も同じ係数を掛ける。
2. そうして、

   $$\begin{aligned}\gamma_{i}&=a_{ij}-\sum_{k=1}^{i-1}\alpha_{ik}\beta_{kj} \qquad (i=j,j+1,\cdots,N) \end{aligned}$$

   
   
   を計算する。これは、式(23)と同じ、 式(24)と $\beta_{jj}$ で割ること以外は同じであ ることに注意が必要である。
3. 最大の $\gamma_{i}$ となるものをピボットとして選択する。
4. 最大のピボットとなる行が分かったので、後は元(規格化前)の $\boldsymbol{L}$ と $\boldsymbol{A}$ を用いて、式(24)と $\beta_{jj}$ を計算する。

これがピボット探索のルーチンである。実際には、ピボット作成用の関数を作成 して、計算をすることになる。
ホームページ: Yamamoto's laboratory
著者: 山本昌志