gijyutsu-keisan.com

1．LU分解

LU分解とは、行列\( A \in R^{m \times n} \ \)を上三角行列L、下三角行列Uの積に分解することです。
\[ A = LU \tag{1-1} \]
成分表示すると、
\[ \begin{pmatrix} & a_{11} & a_{12} & \cdots & a_{1n} & \\ & a_{21} & a_{22} & \cdots & a_{2n} & \\ & \vdots & & \ddots & & \\ & a_{n1} & a_{n2} & \cdots & a_{nn} & \end{pmatrix} = \begin{pmatrix} & l_{11} & 0 & \cdots & 0 & \\ & l_{21} & l_{22} & \ddots & \vdots & \\ & \vdots & & \ddots & 0 & \\ & l_{n1} & \cdots & l_{n,n-1} & l_{nn} & \end{pmatrix} \begin{pmatrix} & u_{11} & u_{12} & \cdots & u_{1n} & \\ & 0 & u_{22} & \cdots & u_{2n} & \\ & \vdots & \ddots & \ddots & \vdots & \\ & 0 & \cdots & 0 & u_{nn} & \end{pmatrix} \]
LU分解は、連立一次方程式を数値計算で解く際に有効な手法です。
LU分解する方法は様々ありますが、本サイトでは次の3つについて説明します。
  • ガウスの消去法
  • クラウト法
  • コレスキー分解

2．ガウスの消去法

ガウスの消去法は、行列\( A \in R^{m \times n} \ \)を上三角化する1つの手法です。 1列ごとに対角成分より下の成分を0にしていくことで上三角化を実現します。
\[ \begin{pmatrix} * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} * & * & * & * \\ 0 & * & * & * \\ 0 & * & * & * \\ 0 & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \\ 0 & 0 & * & * \end{pmatrix} \rightarrow \begin{pmatrix} * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \\ 0 & 0 & 0 & * \end{pmatrix} \]

2．1．前処理なしガウスの消去法

本節では説明を簡単にするため、“対角成分が0でないn次正方行列”に対するガウスの消去法について説明します （“前処理なし”の意味は2.3節で明らかになります。）。
まずは1列目を次のように変換します。
\[ A = \begin{pmatrix} & a_{11} & a_{12} & \cdots & a_{1n} & \\ & a_{21} & a_{22} & \cdots & a_{2n} & \\ & \vdots & & \ddots & & \\ & a_{n1} & a_{n2} & \cdots & a_{nn} & \end{pmatrix} \rightarrow \begin{pmatrix} & a_{11} & a_{12} & \cdots & a_{1n} & \\ & \txc{red}{0} & a^{(1)}_{22} & \cdots & a^{(1)}_{2n} & \\ & \txc{red}{\vdots} & & \ddots & & \\ & \txc{red}{0} & a^{(1)}_{n2} & \cdots & a^{(1)}_{nn} & \end{pmatrix} \]
この変換を実現する行列は、次のように設定します。
1列\( i \)行成分（\( 2 \leq i \leq n \)）を0にするため、1行目を\( -a_{i1} / a_{11} \)倍して、それを\( i \)行目に加えます （\( a_{i1} = 0\)の場合は後ほど説明します）。 この操作は行基本変形行列\( A(i,1: -a_{i1} / a_{11}) \)で実現できます。
\[ A(i,1: -a_{i1} / a_{11} ) = \begin{eqnarray}\left( \begin{array}{c|ccccc} 1 & & & \bf{o} & & & \\ \hline 0 & & & & & & \\ \vdots & & & & & & \\ -a_{i1}/a_{11} & & & E_{n-1} & & & \\ \vdots & & & & & & \\ 0 & & & & & & \end{array} \right)\end{eqnarray} \tag{2.1-1} \]
これを\( i=2 \)から\( i=m \)まで掛け合わせると、
\[ G_1 = \displaystyle \prod_{ i = 2 }^m A(i,1:-a_{i1} / a_{11}) = \begin{eqnarray}\left( \begin{array}{c|ccccc} 1 & & & \bf{o} & & & \\ \hline -a_{21}/a_{11} & & & & & & \\ \vdots & & & E_{n-1} & & & \\ -a_{n1}/a_{11} & & & & & & \end{array} \right)\end{eqnarray} \tag{2.1-2} \]
1列目に対する変換行列が得られます。
この行列\( G_1 \)を\( A \)に作用すると、次のように1列目の、対角成分より下の成分が0になります。
\[ \begin{eqnarray} & G_1 A & = \begin{pmatrix} & 1 & 0 & \cdots & 0 & \\ & -a_{21}/a_{11} & 1 & \cdots & & \\ & \vdots & & \ddots & & \\ & -a_{n1}/a_{11} & & \cdots & 1 & \end{pmatrix} \begin{pmatrix} & a_{11} & a_{12} & \cdots & a_{1n} & \\ & a_{21} & a_{22} & \cdots & a_{2n} & \\ & \vdots & & \ddots & & \\ & a_{n1} & a_{n2} & \cdots & a_{nn} & \end{pmatrix} \\ \\ & & = \begin{pmatrix} & a_{11} & a_{12} & \cdots & a_{1n} & \\ & -a_{11}a_{21}/a_{11} + a_{21} & -a_{12}a_{21}/a_{11} + a_{22} & \cdots & -a_{1n}a_{21}/a_{11} + a_{2n} & \\ & \vdots & & \ddots & & \\ & -a_{11}a_{n1}/a_{11} + a_{n1} & -a_{12}a_{n1}/a_{11} + a_{n2} & \cdots & -a_{1n}a_{n1}/a_{11} + a_{nn} & \end{pmatrix} \\ \\ & & = \begin{pmatrix} & a_{11} & a_{12} & \cdots & a_{1n} & \\ & \txc{blue}{0} & \txc{red}{a^{(1)}_{22}} & \txc{red}{\cdots} & \txc{red}{a^{(1)}_{2n}} & \\ & \txc{blue}{\vdots} & & \txc{red}{\ddots} & & \\ & \txc{blue}{0} & \txc{red}{a^{(1)}_{n2}} & \txc{red}{\cdots} & \txc{red}{a^{(1)}_{nn}} & \end{pmatrix} \end{eqnarray} \tag{2.1-3} \]
このとき1行目は計算する必要がありません。 また1列目は0になることがわかっているので、実質計算を行うのは赤字の範囲のみです。 つまり、本来はn×n個の計算が必要なところを(n-1)×(n-1)個の計算で事足ります。
同様にして2列目に対する変換行列を作ると、次のようになります。
\[ G_2 = \displaystyle \prod_{ i = 3 }^n A(i,2:-a_{i2} / a_{22}) = \begin{eqnarray}\left( \begin{array}{cc|cccc} 1 & 0 & & & O & & \\ 0 & 1 & & & & & \\ \hline 0 & -a_{32}/a_{22} & & & & & \\ \vdots & \vdots & & & E(n-2) & & \\ 0 & -a_{n2}/a_{22} & & & & & \end{array} \right)\end{eqnarray} \tag{2.1-4} \]
この行列\( G_2 \)を\( G_1 A \)に作用すると、2列目の対角成分より下の成分が0になります。
\[ G_2 G_1 A = \begin{pmatrix} & a_{11} & a_{12} & \cdots & & a_{1n} & \\ & 0 & a^{(1)}_{22} & \cdots & & a^{(1)}_{2n} & \\ & 0 & \txc{blue}{0} & \txc{red}{a^{(2)}_{33}} & \txc{red}{\cdots} & \txc{red}{a^{(2)}_{3n}} & \\ & \vdots & \txc{blue}{\vdots} & & \txc{red}{\ddots} & & \\ & 0 & \txc{blue}{0} & \txc{red}{a^{(2)}_{3n}} & \txc{red}{\cdots} & \txc{red}{a^{(2)}_{nn}} & \end{pmatrix} \tag{2.1-5} \]
このとき1、2行目は計算する必要がありません。 また2列目は0になることがわかっているので、実質計算を行うのは赤字の範囲のみです。 つまり、本来はn×n個の計算が必要なところを(n-2)×(n-2)個の計算で事足ります。
これをn-1列まで繰り返し行うことで、\( A \)は上三角行列\( U \)に変換されます。
\[ G_{n-1} \cdots G_2 G_1 A = \begin{pmatrix} & a_{11} & a_{12} & \cdots & & a_{1n} & \\ & & a^{(1)}_{22} & \cdots & & a^{(1)}_{2n} & \\ & & & a^{(2)}_{33} & & a^{(2)}_{3n} & \\ & & O & & \ddots & \vdots & \\ & & & & & a^{(n-1)}_{nn} & \end{pmatrix} =U \tag{2.1-6} \]
このとき、実質計算する個数は
  • 行列計算として行う場合、\( n^2 \times (n-1) \)
  • ガウスの消去法として行う場合、\( \displaystyle\prod_{i = 1}^{(n-1)} (n-i)^2 \)
であり、ガウスの消去法は計算回数を減らすメリットがあります。

また、変換行列の積\( G_{n-1} \cdots G_2 G_1 \)が下三角行列になることは成分計算することで確認できます。 また、この変換行列は基本変換行列の積で構成されているため、逆行列を持ちます。 従って、この変換行列を\( L^{-1} \)で表すことにします。
\[ G_{n-1} \cdots G_2 G_1 = \begin{pmatrix} & 1 & & & & \\ & * & 1 & & & \\ & \vdots & \ddots & \ddots & & \\ & * & \cdots & * & 1 & \end{pmatrix} = L^{-1} \tag{2.1-7} \]
（2.1-6）式と（2.1-7）式を合わせると
\[ L^{-1} A = U \quad \leftrightarrow \quad A = L U \tag{2.1-8} \]
で表せ、行列\( A \)は下三角行列\( L \)と上三角行列\( U \)に分解できることがわかります。 この分解をLU分解といいます。 つまり、ガウスの消去法は“行列をLU分解する”ことに他なりません。

2．2．ピボッティングとスケーリング

前節のガウスの消去法では、「行列Aの対角成分が0でない」ことを前提に話を進めました。 しかしこれは保証されず、対角成分=0のとき計算不能となります。 また、対角成分の絶対値が非常に小さい場合、割り算による計算精度上の問題が生じます。 これらの問題を回避するために、“対角成分≠0”となるよう、行または列の入れ替え操作を行います。 この入れ替え操作をピボッティングと呼び、以下に示す方法で行います。
  • 行ピボッティング：\( k \)列目\( k+1～n \)行目の中から絶対値最大となる成分を含む\( i \)行と\( k \)行を交換
  • 列ピボッティング：\( k \)行目\( k+1～n \)列目の中から絶対値最大となる成分を含む\( j \)列と\( k \)列を交換
  • 完全ピボッティング：行ピボッティングと列ピボッティングを両方行う
さらに、各行または列成分の絶対値に極端な違いがあるとき、ピボッティングを行っても計算精度上の問題を回避できないことがあります。 この問題を回避するため、「各行内の絶対値最大成分でその他の成分を割る」操作を行います。 これをスケーリングと呼びます。
ピボッティングとスケーリングは基本変形行列によって実行可能です。

（1）行ピボッティングの場合

\( i \)行と\( j \)行を入れ替えるために、 交換行列\( P(i,j) \)を左から掛けます。
\[ \begin{eqnarray} & P(i,j) A \quad & = \begin{pmatrix} 0 & 0 & \txc{red}{1} & 0 \\ 0 & 1 & 0 & 0 \\ \txc{red}{1} & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{pmatrix} \\ \\ & & = \begin{pmatrix} \txc{red}{a_{31}} & \txc{red}{a_{32}} & \txc{red}{a_{33}} & \txc{red}{a_{34}} \\ a_{12} & a_{22} & a_{23} & a_{24} \\ \txc{red}{a_{11}} & \txc{red}{a_{12}} & \txc{red}{a_{13}} & \txc{red}{a_{14}} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{pmatrix} \end{eqnarray} \]

（2）列ピボッティングの場合

\( i \)列と\( j \)列を入れ替えるために、 交換行列\( P(i,j) \)を右から掛けます。
\[ \begin{eqnarray} & P(i,j) A \quad & = \begin{pmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{pmatrix} \begin{pmatrix} 0 & 0 & \txc{red}{1} & 0 \\ 0 & 1 & 0 & 0 \\ \txc{red}{1} & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \\ \\ & & = \begin{pmatrix} \txc{red}{a_{13}} & a_{12} & \txc{red}{a_{11}} & a_{14} \\ \txc{red}{a_{23}} & a_{22} & \txc{red}{a_{21}} & a_{24} \\ \txc{red}{a_{33}} & a_{32} & \txc{red}{a_{31}} & a_{34} \\ \txc{red}{a_{43}} & a_{42} & \txc{red}{a_{41}} & a_{44} \end{pmatrix} \end{eqnarray} \]

（3）スケーリングの場合

\( i \)行目を\( 1/a_{ij} \)倍するために、 スカラー倍行列\( M(i:1/a_{ij}) \)を左から掛けます。
\[ \begin{eqnarray} M(i:1/a_{ij}) A & = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \txc{red}{1/a_{23}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{pmatrix} \\ \\ & = \begin{pmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ \txc{red}{a_{12}/a_{23}} & \txc{red}{a_{22}/a_{23}} & \txc{red}{a_{23}/a_{23}} & \txc{red}{a_{24}/a_{23}} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{pmatrix} \end{eqnarray} \]
スケーリングやピボッティングを実施する必要のない場合は単位行列\( E \)を適用します。

参考文献

▲