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}$を求める。
    • まず、 $ i=1,2,\cdots,j$について、次式に従い $ \beta_{ij}$を 計算する。

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

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

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

    • これで、 $ \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{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
平成16年11月9日


no counter