読者です 読者をやめる 読者になる 読者になる

作業中のメモ

よく「計算機」を使って作業をする.知らなかったことを中心にまとめるつもり.

LU 分解による連立一次方程式の解法

前回は,ガウスの消去法について,簡単に説明した.今回は,LU 分解について説明する.

LU 分解

行列  {A}正則行列逆行列を持つもの)であるとする.このとき,行列  {A} を部分ピボット選択をする(この行列を  {P} とする)ことにより,

 { \displaystyle
  PA = LU
}

のように LU 分解できる.ここで,L,U はそれぞれ,

 { \displaystyle
  \begin{equation}
    L = \left(
      \begin{array}{ccccc}
        1 & 0 & 0 & \cdots & 0\\
        l_{2}^{(1)} & 1 & 0 & \cdots & 0\\
        l_{3}^{(1)} & l_{3}^{(2)} & 1 & \cdots & 0\\
        \vdots & \vdots & \vdots & \ddots & \vdots\\
        l_{n}^{(1)} & l_{n}^{(2)} & l_{n}^{(3)} & \cdots & 1\\
      \end{array}
    \right),\quad U = \left(
      \begin{array}{ccccc}
        a_{11}^{(0)} & a_{12}^{(0)} & a_{13}^{(0)} & \cdots & a_{1n}^{(0)}\\
        0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & \cdots a_{2n}^{(1)}\\
        0 & 0 & a_{33}^{(2)} & \cdots & a_{3n}^{(2)}\\
        \vdots & \vdots & \vdots & \ddots & \vdots\\
        0 & 0 & 0 & \cdots & a_{nn}^{(n - 1)}\\
      \end{array}
    \right)
  \end{equation}
}

である.ここから,  {A\boldsymbol{x} = \boldsymbol{b}} を解くことを考える.まず, {L \boldsymbol{y} = P\boldsymbol{b}} を解く.すなわち,

 { \displaystyle
  \begin{equation}
    \left(
      \begin{array}{ccccc}
        1 & 0 & 0 & \cdots & 0\\
        l_{2}^{(1)} & 1 & 0 & \cdots & 0\\
        l_{3}^{(1)} & l_{3}^{(2)} & 1 & \cdots & 0\\
        \vdots & \vdots & \vdots & \ddots & \vdots\\
        l_{n}^{(1)} & l_{n}^{(2)} & l_{n}^{(3)} & \cdots & 1\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        y_{1}\\
        y_{2}\\
        y_{3}\\
        \vdots\\
        y_{n}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        b_{p_{1}}\\
        b_{p_{2}}\\
        b_{p_{3}}\\
        \vdots\\
        b_{p_{n}}\\
      \end{array}
    \right)
  \end{equation}
}

を解くことになる.ここで, {b_{p_{k}}} は,部分ピボット選択により,入れ替えたものを表す.入れ替えが行われなかった場合は, {b_{p_{k}} = b_{k}} である. この方程式は,前進代入で解くことが出来る.つまり,

 { \displaystyle
  \begin{equation}
    y_{k} = b_{p_{k}} - \sum_{j = 1}^{k - 1} l_{k}^{(j)} y_{j},\quad (k = 1, \cdots, n)
  \end{equation}
}

で求められる.次に, {U \boldsymbol{x} = \boldsymbol{y}} を解く.すなわち,

 { \displaystyle
  \begin{equation}
    \left(
      \begin{array}{ccccc}
        a_{11}^{(0)} & a_{12}^{(0)} & a_{13}^{(0)} & \cdots & a_{1n}^{(0)}\\
        0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & \cdots a_{2n}^{(1)}\\
        0 & 0 & a_{33}^{(2)} & \cdots & a_{3n}^{(2)}\\
        \vdots & \vdots & \vdots & \ddots & \vdots\\
        0 & 0 & 0 & \cdots & a_{nn}^{(n - 1)}\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
        \vdots\\
        x_{n}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        y_{1}\\
        y_{2}\\
        y_{3}\\
        \vdots\\
        y_{n}\\
      \end{array}
    \right)
  \end{equation}
}

を解く.これは,ガウスの消去法のときと同様に,後退代入で解くことが出来る.つまり,

 { \displaystyle
  \begin{equation}
    x_{k} = \dfrac{1}{a_{kk}^{(k - 1)}} \left( y_{k} - \sum_{j = k + 1}^{n} a_{kj}^{(k - 1)}x_{j} \right),\quad (k = n - 1, \cdots, 1)
  \end{equation}
}

で求められる.

アルゴリズム

これを,C 言語で書き下すと以下のようになる.

参考文献

ガウスの消去法,LU 分解 http://www.ms.u-tokyo.ac.jp/~narutaka/lec/lud.pdf

LU 分解法 http://www.kata-lab.itc.u-tokyo.ac.jp/OpenLecture/SP20121218.pdf