4 LU分解

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

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

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

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

$$\begin{pmatrix}\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{pmatrix}$$ (19)

となります。

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

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

となります。これをさらに書き換えると、

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

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

となります。これらの連立方程式の解 $\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回で済みます。これが、ガウスの消去法と後退代入を組 み合わせた方法やガウス・ジョルダン法に比べて、際立って優れている点です。

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

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

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

とします。

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

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

   $$\begin{pmatrix}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{pmatrix}$$ (26)

   
   
   と変形できます。
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分解ができます。次に示すピボット選択をしなければ、ア ルゴリズムは非常に単純です。

4.3 ピボット選択

ここでも、ピボット選択の問題が出てきます。式(28)の $\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}$$

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

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