4 LU分解

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

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

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

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

$\displaystyle \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次方程式は

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

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

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

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

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

\begin{equation*}\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}\end{equation*}

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

\begin{equation*}\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}\end{equation*}

これらの前進代入と後退代入は、コンピューターにとって非常に簡単に計算で きます。これは、係数行列$ 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)のアルゴリズムに従い、

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

とします。

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

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


    と変形できます。
  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}$を求めます。

この方法により、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{equation*}\begin{aligned}\gamma_{i}&=a_{ij}-\sum_{k=1}^{i-1}\alpha_{ik}\beta_{kj} \qquad (i=j,j+1,\cdots,N) \end{aligned}\end{equation*}

    を計算します。これは、式(27)と同じ、 式(28)と $ \beta_{jj}$で割ること以外は同じであ ることに注意してください。
  3. 最大の $ \gamma_{i}$となるものをピボットと選択します。
  4. 最大のピボットとなる行が分かったので、後は元(規格化前)の $ \boldsymbol{L}$ $ \boldsymbol{A}$を用いて、式(28)と $ \beta_{jj}$を計算します。
これがピボット探索のルーチンです。実際には、ピボット作成用の関数を作成 して、計算をすることになります。
ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成19年8月21日


no counter