LU分解
2．3．前処理付きガウスの消去法
本節では、n次正方行列に前処理としてピボッティングやスケーリングを施すガウスの消去法について説明します。ガウスの消去法を適用する際、ピボッティングやスケーリングを施すことで、0割りや計算精度上の問題を回避できます （確実ではありませんが）。
そこで、行ピボッティングを\( Pr_k \)、列ピボッティングを\( Pc_k \)、スケーリングを\( M_k \)として、これらを行列\( A \)に作用させると、ガウスの消去法は次のようになります。
\[
( G_{n-1} M_{n-1} Pr_{n-1} \cdots G_2 M_2 Pr_2 G_1 M_1 Pr_1 )
A
( Pc_1 Pc_2 \cdots Pc_{n-1} )
=
U
\tag{2.3-1}
\]
\( G_k \)を作用させる前に、行ピボッティングを\( Pr_k \)、列ピボッティング\( Pc_k \)、スケーリング\( M_k \)を行うことから、これらの操作を前処理といいます。
行側の前処理を1つにまとめて\( Tr_k = M_k Pr_k \)で表すことにします。
以下に4次元正方行列での計算例を挙げます（列ピボッティングは省略します）。
1列ごとに前処理を施しつつガウスの消去法を行うので、前処理付きガウスの消去法は、次のような行列表現で表せます。
\[
( G_4 Tr_4 G_3 Tr_3 G_2 Tr_2 G_1 Tr_1 ) A = U
\tag{2.3-2}
\]
＜1列目の計算＞
1行目と3行目を交換して、交換後の1行目を\( c_1 \)倍します （これはあくまで例です）。\[
\begin{eqnarray}
& Tr_1 &
=
\begin{pmatrix}
c_1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & 1 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
0 & 0 & 1 & 0
\\
0 & 1 & 0 & 0
\\
1 & 0 & 0 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
=
\begin{pmatrix}
0 & 0 & c_1 & 0
\\
0 & 1 & 0 & 0
\\
1 & 0 & 0 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\\
\\
& Tr_1 A &
=
\begin{pmatrix}
0 & 0 & c_1 & 0
\\
0 & 1 & 0 & 0
\\
1 & 0 & 0 & 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}
c_1 a_{31} & c_1 a_{32} & c_1 a_{33} & c_1 a_{34}
\\
a_{21} & a_{22} & a_{23} & a_{24}
\\
a_{11} & a_{12} & a_{13} & a_{14}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\\
\\
& G_1 ( Tr_1 A ) \ &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
-a_{21}/c_1 a_{31} & 1 & 0 & 0
\\
-a_{11}/c_1 a_{31} & 0 & 1 & 0
\\
-a_{41}/c_1 a_{31} & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
c_1 a_{31} & c_1 a_{32} & c_1 a_{33} & c_1 a_{34}
\\
a_{21} & a_{22} & a_{23} & a_{24}
\\
a_{11} & a_{12} & a_{13} & a_{14}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
c_1 a_{31} & c_1 a_{32} & c_1 a_{33} & c_1 a_{34}
\\
0 & \txc{red}{a_{22}^{(1)}} & \txc{red}{a_{23}^{(1)}} & \txc{red}{a_{24}^{(1)}}
\\
0 & \txc{red}{a_{32}^{(1)}} & \txc{red}{a_{33}^{(1)}} & \txc{red}{a_{34}^{(1)}}
\\
0 & \txc{red}{a_{42}^{(1)}} & \txc{red}{a_{43}^{(1)}} & \txc{red}{a_{44}^{(1)}}
\end{pmatrix}
=
A^{(1)}
\quad
( a_{ij}^{(1)} = -a_{3j} a_{i1}/a_{31} +a_{ij} \ )
\end{eqnarray}
\]
実質の計算は、スケーリングを除けば赤文字範囲のみで済みます。
＜2列目の計算＞
2行目と4行目を交換して、交換後の2行目を\( c_2 \)倍します （これはあくまで例です）。\[
\begin{eqnarray}
& Tr_2 &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & c_2 & 0 & 0
\\
0 & 0 & 1 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 0 & 0 & 1
\\
0 & 0 & 1 & 0
\\
0 & 1 & 0 & 0
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 0 & 0 & c_2
\\
0 & 0 & 1 & 0
\\
0 & 1 & 0 & 0
\end{pmatrix}
\\
\\
& Tr_2 A^{(1)} &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 0 & 0 & c_2
\\
0 & 0 & 1 & 0
\\
0 & 1 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)}
\\
0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)}
\\
0 & a_{42}^{(1)} & a_{43}^{(1)} & a_{44}^{(1)}
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)}
\\
0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)}
\end{pmatrix}
\\
\\
& G_2 ( Tr_2 A^{(1)} ) &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & -a_{32}^{(1)}/c_2 a_{42}^{(1)} & 1 & 0
\\
0 & -a_{22}^{(1)}/c_2 a_{42}^{(1)} & 0 & 1
\end{pmatrix}
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)}
\\
0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)}
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & 0 & \txc{red}{a_{33}^{(2)}} & \txc{red}{a_{34}^{(2)}}
\\
0 & 0 & \txc{red}{a_{43}^{(2)}} & \txc{red}{a_{44}^{(2)}}
\end{pmatrix}
=
A^{(2)}
\quad
( a_{ij}^{(2)} = -a_{4j}^{(1)} a_{i2}^{(1)}/a_{42}^{(1)} +a_{ij}^{(1)} \ )
\end{eqnarray}
\]
実質の計算は、スケーリングを除けば赤文字範囲のみで済みます。
＜3列目の計算＞
3行目と4行目を交換して、交換後の3行目を\( c_3 \)倍します （これはあくまで例です）。\[
\begin{eqnarray}
& Tr_3 &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & c_3 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & 0 & 1
\\
0 & 0 & 1 & 0
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & 0 & c_3
\\
0 & 0 & 1 & 0
\end{pmatrix}
\\
\\
& Tr_3 A^{(2)} &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & 0 & c_3
\\
0 & 0 & 1 & 0
\end{pmatrix}
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & 0 & a_{33}^{(2)} & a_{34}^{(2)}
\\
0 & 0 & a_{43}^{(2)} & a_{44}^{(2)}
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & 0 & c_3 a_{43}^{(2)} & c_3 a_{44}^{(2)}
\\
0 & 0 & a_{33}^{(2)} & a_{34}^{(2)}
\end{pmatrix}
\\
\\
& G_3 ( Tr_3 A^{(2)} ) &
=
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & 1 & 0
\\
0 & 0 & -a_{33}^{(2)}/c_3 a_{43}^{(2)} & 1
\end{pmatrix}
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & 0 & c_3 a_{43}^{(2)} & c_3 a_{44}^{(2)}
\\
0 & 0 & a_{33}^{(2)} & a_{34}^{(2)}
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
a_{31} & a_{32} & a_{33} & a_{34}
\\
0 & c_2 a_{42}^{(1)} & c_2 a_{43}^{(1)} & c_2 a_{44}^{(1)}
\\
0 & 0 & c_3 a_{43}^{(2)} & c_3 a_{44}^{(2)}
\\
0 & 0 & 0 & \txc{red}{-a_{33}^{(2)} a_{44}^{(2)}/a_{43}^{(2)} + a_{34}^{(2)}} \
\end{pmatrix}
=
U
\end{eqnarray}
\]
実質の計算は、スケーリングを除けば赤文字範囲のみで済みます。
以上で、行列\( A \)は上三角行列\( U \)に変形できました。
次に、下三角行列\( L \)について見ていきます。 そこで、（2.3-2）式を次のように変形します。
\[
\begin{eqnarray}
A^{(0)}
& = &
A
\\
\\
A^{(1)}
& = &
G_1 Tr_1 A^{(0)}
\\
& = & L_1 ( Tr_1 A )
\\
\\
A^{(2)}
& = &
G_2 Tr_2 A^{(1)}
\\
& = &
G_2 Tr_2 ( L_1 Tr_1 A )
\\
& = &
G_2 ( Tr_2 L_1 Tr_2^{-1} ) Tr_2 Tr_1 A
\\
& = &
L_2 ( Tr_2 Tr_1 A )
\\
\\
A^{(3)}
& = &
G_3 Tr_3 A^{(2)}
\\
& = &
G_3 Tr_3 ( L_2 Tr_2 Tr_1 A )
\\
& = &
G_3 ( Tr_3 L_2 Tr_3^{-1} ) Tr_3 Tr_2 Tr_1 A
\\
& = &
L_3 ( Tr_3 Tr_2 Tr_1 A )
\\
& = &
U
\end{eqnarray}
\]
ここで\( Tr_2 L_1 Tr_2^{-1} \)を計算すると、
\[
\begin{eqnarray}
Tr_2 L_1 Tr_2^{-1}
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 0 & 0 & c_2
\\
0 & 0 & 1 & 0
\\
0 & 1 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
l_{21}^{(1)} & 1 & 0 & 0
\\
l_{31}^{(1)} & 0 & 1 & 0
\\
l_{41}^{(1)} & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 0 & 0 & 1
\\
0 & 0 & 1 & 0
\\
0 & 1/c_2 & 0 & 0
\end{pmatrix}
\\
\\
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
c_2 l_{41}^{(1)} & 0 & 0 & c_2
\\
l_{31}^{(1)} & 0 & 1 & 0
\\
l_{21}^{(1)} & 1 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 0 & 0 & 1
\\
0 & 0 & 1 & 0
\\
0 & 1/c_2 & 0 & 0
\end{pmatrix}
\\
\\
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
c_2 l_{41}^{(1)} & 1 & 0 & 0
\\
l_{31}^{(1)} & 0 & 1 & 0
\\
l_{21}^{(1)} & 0 & 0 & 1
\end{pmatrix}
\end{eqnarray}
\]
行列\( L_1 \)に対して、成分の0、非0構造は変わりません。
これは\( Tr_3 L_2 Tr_3^{-1} \)も同様です。
以上を考慮すると、\( L_2,L_3 \)は
\[
\begin{eqnarray}
L_2
=
G_2 ( Tr_2 L_1 Tr_2^{-1} )
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & l_{32}^{(2)} & 1 & 0
\\
0 & l_{42}^{(2)} & 0 & 1
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
c_2 l_{41}^{(1)} & 1 & 0 & 0
\\
l_{31}^{(1)} & 0 & 1 & 0
\\
l_{21}^{(1)} & 0 & 0 & 1
\end{pmatrix}
\\
\\
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
c_2 l_{41}^{(1)} & 1 & 0 & 0
\\
c_2 l_{41}^{(1)} l_{32}^{(2)} + l_{31}^{(1)} & l_{32}^{(2)} & 1 & 0
\\
c_2 l_{41}^{(1)} l_{42}^{(2)} + l_{21}^{(1)} & l_{42}^{(2)} & 0 & 1
\end{pmatrix}
\\
\\
L_3
=
G_3 ( Tr_3 L_2 Tr_3^{-1} )
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & 1 & 0 & 0
\\
0 & 0 & 1 & 0
\\
0 & 0 & l_{34}^{(3)} & 1
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0
\\
l_{21}^{(3)} & 1 & 0 & 0
\\
l_{31}^{(3)} & l_{32}^{(3)} & 1 & 0
\\
l_{41}^{(3)} & l_{42}^{(3)} & 0 & 1
\end{pmatrix}
\\
& = &
\begin{pmatrix}
1 & 0 & 0 & 0
\\
l_{21}^{(3)} & 1 & 0 & 0
\\
l_{31}^{(3)} & l_{32}^{(3)} & 1 & 0
\\
l_{31}^{(3)} l_{34}^{(3)} + l_{41}^{(3)} & l_{34}^{(2)} l_{32}^{(3)} + l_{42}^{(3)} & l_{34}^{(3)} & 1
\end{pmatrix}
\end{eqnarray}
\]
となって、\( L_3 \)は下三角行列になります。
下三角行列の逆行列もまた下三角行列になるため、
\[
L_3 A = U
\
\rightarrow
\
A = L_3^{-1}U
= LU
\quad ( L = L_3^{-1} )
\tag{2.3-3}
\]
行列\( A \)は、前処理を施してもLU分解可能であることがわかります。
2．4．非正方行列のガウスの消去法
これまでは正方行列に対するガウスの消去法について説明してきました。 ただ、ガウスの消去法は非正方行列に対しても適用できます。（1）\( A \in R^{m \times n} \ \)が横長行列の場合\( \ (m \lt n)\)
\[
A
=
\begin{pmatrix}
a_{11} & \cdots & a_{1m} & \cdots & a_{1n}
\\
\vdots & & \vdots & & \vdots
\\
a_{m1} & \cdots & a_{mm} & \cdots & a_{mn}
\end{pmatrix}
\]
横長行列にガウスの消去法を適用すると、次のようになります
（例として（3×4）行列を使います）。
\[
\begin{eqnarray}
G_1 A
& = &
\begin{pmatrix}
1 & 0 & 0
\\
-a_{21}/a_{11} & 1 & 0
\\
-a_{31}/a_{11} & 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}
\end{pmatrix}
\\
\\
& = &
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)}
\\
0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)}
\end{pmatrix}
\end{eqnarray}
\]
\[
\begin{eqnarray}
G_2 G_1 A
& = &
\begin{pmatrix}
1 & 0 & 0
\\
0 & 1 & 0
\\
0 & -a_{32}^{(1)}/a_{22}^{(1)} & 1
\end{pmatrix}
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)}
\\
0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)}
\end{pmatrix}
\\
\\
& = &
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)}
\\
0 & 0 & a_{33}^{(1)} & a_{34}^{(1)}
\end{pmatrix}
\end{eqnarray}
\]
このように、横長行列にガウスの消去法を適用すると、LU分解は次のような形になります。
\[
\begin{eqnarray}
A
& = &
LU
\\
& = &
\begin{pmatrix}
l_{11} & & O
\\
\vdots & \ddots &
\\
l_{m1} & \cdots & l_{mm}
\end{pmatrix}
\begin{pmatrix}
u_{11} & \cdots & u_{1m} & \cdots & u_{1n}
\\
& \ddots & \vdots & & \vdots
\\
O & & u_{mm} & \cdots & u_{mn}
\end{pmatrix}
\end{eqnarray}
\tag{2.4-1}
\]
なお、前処理行列について、行用はm次正方行列\( Tr \)、列用はn次正方行列\( Tc \)になります。
\[
T_r A T_c
=
\begin{pmatrix}
tr_{11} & \cdots & tr_{1m}
\\
\vdots & & \vdots
\\
tr_{m1} & \cdots & tr_{mm}
\end{pmatrix}
\begin{pmatrix}
a_{11} & \cdots & a_{1m} & \cdots & a_{1n}
\\
\vdots & & \vdots & & \vdots
\\
a_{m1} & \cdots & a_{mm} & \cdots & a_{mn}
\end{pmatrix}
\begin{pmatrix}
tc_{11} & \cdots & \cdots & tc_{1n}
\\
\vdots & & & \vdots
\\
\vdots & & & \vdots
\\
tc_{n1} & \cdots & \cdots & tc_{nn}
\end{pmatrix}
\]
（2）\( A \in R^{m \times n} \ \)が縦長行列の場合\( \ (m \gt n)\)
\[
A
=
\begin{pmatrix}
a_{11} & \cdots & a_{1n}
\\
\vdots & & \vdots
\\
a_{n1} & \cdots & a_{nn}
\\
\vdots & & \vdots
\\
a_{m1} & \cdots & a_{mn}
\end{pmatrix}
\]
縦長行列にガウスの消去法を適用する場合、計算にはひと手間必要です。始めに\( A \)を転置することで、縦長行列を横長行列にします。
\[
A^t
=
\begin{pmatrix}
a_{11} & \cdots & a_{1n} & \cdots & a_{1m}
\\
\vdots & & \vdots & & \vdots
\\
a_{n1} & \cdots & a_{nn} & \cdots & a_{nm}
\end{pmatrix}
\]
この転置行列\( A^t \)に対してガウスの消去法を適用します。
\[
A^t =
LU
=
\begin{pmatrix}
l_{11} & & O
\\
\vdots & \ddots &
\\
l_{n1} & \cdots & l_{nn}
\end{pmatrix}
\begin{pmatrix}
u_{11} & \cdots & u_{1n} & \cdots & u_{1m}
\\
& \ddots & \vdots & & \vdots
\\
O & & u_{nn} & \cdots & u_{nm}
\end{pmatrix}
\]
この結果に対してさらに転置をとると、\( (A^t)^t=A \)、\( (LU)^t=U^tL^t \)を得ます。
これを成分表示すると、
\[
A
=
U^t L^t
=
\begin{pmatrix}
u_{11} & & O
\\
\vdots & \ddots &
\\
u_{1n} & \cdots & u_{nn}
\\
\vdots & & \vdots
\\
u_{1m} & & u_{nm}
\end{pmatrix}
\begin{pmatrix}
l_{11} & \cdots & l_{n1}
\\
& \ddots & \vdots
\\
O & \cdots & l_{nn}
\end{pmatrix}
\tag{2.4-2}
\]
となり、結果\( U^t \)は下三角行列に、\( L^t \)は上三角行列になります。また、前処理付きガウスの消去法を行列表現すれば、
\[
G_n ( Tr_n L_{n-1} Tr_n^{-1} ) Tr_n \cdots Tr_1 A^t Tc_1 \cdots Tc_n
=
U
\tag{2.4-3}
\]
となります。このとき、\( L^{-1} = G_n ( Tr_n L_{n-1} Tr_n^{-1} ) \)であることを考慮して、上式の転置をとると
\[
( Tr_n \cdots Tr_1 A^t Tc_1 \cdots Tc_n )^t = ( LU )^t
\\
\rightarrow
\
(Tc_1 \cdots Tc_n )^t ( A^t )^t (Tr_n \cdots Tr_1)^t = U^t L^t
\\
\rightarrow
\
(Tc_n \cdots Tc_1 ) A (Tr_1 \cdots Tr_n) = U^t L^t
\]
と変形でき、\( U^t \)は下三角行列、\( L^t \)は上三角行列になります。
この結果、行基本変形行列は\( Tc_n \cdots Tc_1 \)、列基本変形行列は\( Tr_1 \cdots Tr_n \)となります。
