どのようにA行列をLU行列へと分解するのか – アルゴリズム –

\(LU\)分解は、方程式\(Ax = b\)の解を計算に工夫が必要な逆行列を求める方法でなく、機械的に簡単な計算で求められる方法です。

ではどのように\(A = LU\)へと分解するのか、具体的に計算を行いながら分解アルゴリズムを見てみましょう。Pivot付きの計算も行うので少し複雑ですが、実装してみれば簡単なことが分かります。

LU分解の基本的な方法

まずは基本的な行列で\(LU\)分解を行ってみます。

\(A = LU \)の行列を(1)式のように置いて考えます。

$$
\left(\begin{array}{cccc}
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{array}\right) =
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{cccc}
u_{11} & u_{12} & u_{13} & u_{14} \\
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{1}
$$

1段目

この(1)式のA行列の1行目の計算から、

$$
\left(\begin{array}{cccc}
a_{11} & a_{12} & a_{13} & a_{14} \\
\hline
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{array}\right) =
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
\hline
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{c|c|c|c}
u_{11} & u_{12} & u_{13} & u_{14} \\
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right)
$$

U行列の1行目は(2)式のように求まります。

$$
\begin{eqnarray}
a_{11} &=& u_{11} \\
a_{12} &=& u_{12} \\
a_{13} &=& u_{13} \\
a_{14} &=& u_{14}
\end{eqnarray} \tag{2}
$$

続いてA行列の1列目から

$$
\left(\begin{array}{c|ccc}
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{array}\right) =
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\ \hline
l_{21} & 1 & 0 & 0 \\ \hline
l_{31} & l_{32} & 1 & 0 \\ \hline
l_{41} & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{c|ccc}
u_{11} & u_{12} & u_{13} & u_{14} \\
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right)
$$

L行列の1列目は(3)式のように求まります。

$$
\begin{eqnarray}
a_{11} = l_{11} u_{11} &\Rightarrow& l_{11} = a_{11} / u_{11} = 1 \\
a_{21} = l_{21} u_{11} &\Rightarrow& l_{21} = a_{21} / u_{11} \\
a_{31} = l_{31} u_{11} &\Rightarrow& l_{31} = a_{31} / u_{11} \\
a_{41} = l_{41} u_{11} &\Rightarrow& l_{41} = a_{41} / u_{11}
\end{eqnarray} \tag{3}
$$

ここで\(l_1, u_1\)を(4)式のように定義すると

$$
l_1 = (1, l_{21}, l_{31}, l_{41})^T, u_1 = (u_{11}, u_{12}, u_{13}, u_{14}) \tag{4}
$$

\(l_1 u_1\)は(5)式となるので、

$$
\begin{eqnarray}
l_1 u_1 &=&
\left(\begin{array}{c}
1 \\ l_{21} \\ l_{31} \\ l_{41}
\end{array}\right)
\left(\begin{array}{cccc}
u_{11} & u_{12} & u_{13} & u_{14}
\end{array}\right) \\ &=&
\left(\begin{array}{cccc}
u_{11} & u_{12} & u_{13} & u_{14} \\
l_{21} u_{11} & l_{21} u_{12} & l_{21} u_{13} & l_{21} u_{14} \\
l_{31} u_{11} & l_{31} u_{12} & l_{31} u_{13} & l_{31} u_{14} \\
l_{41} u_{11} & l_{41} u_{12} & l_{41} u_{13} & l_{41} u_{14}
\end{array}\right)\\ &=&
\left(\begin{array}{cccc}
a_{11} & a_{12} & a_{13} & a_{14} \\
a_{21} & l_{21} u_{12} & l_{21} u_{13} & l_{21} u_{14} \\
a_{31} & l_{31} u_{12} & l_{31} u_{13} & l_{31} u_{14} \\
a_{41} & l_{41} u_{12} & l_{41} u_{13} & l_{41} u_{14}
\end{array}\right)
\end{eqnarray} \tag{5}
$$

\(A – l_1 u_1 = A_{(2)}\)として、\(A_{(2)}\)は(6)式となります。

$$
A_{(2)} =
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & a_{22(2)} & a_{23(2)} & a_{24(2)} \\
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right) \tag{6}
$$

以上の関係から\(A = A_{(1)}\)として、\(A_{(1)} = l_1 u_1 + A_{(2)}\)となります。

2段目

\(A_{(2)}\)は\(A_{(1)}\)から一つ縮小した関係なので、これをさらにLU分解を行います。\(A_{(2)}\)より、LUも(7)式のように一つ縮小した形を考えます。

$$
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & a_{22(2)} & a_{23(2)} & a_{24(2)} \\
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right)=
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 1 & 0 & 0 \\
0 & l_{32} & 1 & 0 \\
0 & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{7}
$$

この(7)式より、1段目と同じようにU行列の2行目は(8)式、L行列の2列目は(9)式のように求まります。

$$
\begin{eqnarray}
a_{22(2)} &=& u_{22} \\
a_{23(2)} &=& u_{23} \\
a_{24(2)} &=& u_{24}
\end{eqnarray} \tag{8}
$$

$$
\begin{eqnarray}
a_{22(2)} = l_{22} u_{22} &\Rightarrow& l_{22} = a_{22(2)} / u_{22} = 1 \\
a_{32(2)} = l_{32} u_{22} &\Rightarrow& l_{32} = a_{32(2)} / u_{22} \\
a_{42(2)} = l_{42} u_{22} &\Rightarrow& l_{42} = a_{42(2)} / u_{22}
\end{eqnarray} \tag{9}
$$

ここで、\(l_2, u_2\)を(10)式のように定義して、

$$
l_2 = (0, 1, l_{32}, l_{42})^T, u_1 = (0, u_{22}, u_{23}, u_{24}) \tag{10}
$$

\(A_{(2)} – l_2 u_2 = A_{(3)}\)とすると\(A_{(3)}\)は(11)式のよう一つ縮小した形になります。

$$
A_{(3)} =
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{33(3)} & a_{34(3)} \\
0 & 0 & a_{43(3)} & a_{44(3)}
\end{array}\right) \tag{11}
$$

\(A_{(2)} = l_2 u_2 + A_{(3)}\)より、\(A_{(1)} = l_1 u_1 + l_2 u_2 + A_{(3)}\)という関係になります。

3段目

\(A_{(2)}\)から一つ縮小した\(A_{(3)}\)をさらにLU分解を行います。\(A_{(3)}\)の関係から\(LU\)行列も(12)式のようになります。

$$
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{33(3)} & a_{34(3)} \\
0 & 0 & a_{43(3)} & a_{44(3)}
\end{array}\right)=
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 1 & 0 \\
0 & 0 & l_{43} & 1
\end{array}\right)
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{12}
$$

(12)式の関係から\(U\)行列の3行目は(13)式、\(L\)行列の3列目は(14)式のよう求まります。

$$
\begin{eqnarray}
a_{33(3)} &=& u_{33} \\
a_{34(3)} &=& u_{34}
\end{eqnarray} \tag{13}
$$

$$
\begin{eqnarray}
a_{33(3)} = l_{33} u_{33} &\Rightarrow& l_{33} = a_{33(3)} / u_{33} = 1 \\
a_{43(3)} = l_{43} u_{33} &\Rightarrow& l_{43} = a_{43(3)} / u_{33}
\end{eqnarray} \tag{14}
$$

ここで、\(l_3, u_3\)を(15)式のように定義して、

$$
l_3 = (0, 0 , 1 , l_{43}) ^T , u_3 = (0, 0 , u_{33}, u_{34}) \tag{15}
$$

\(A_{(3)} – l_3 u_3 = A_{(4)}\)とすると\(A_{(4)}\)は(16)式のように一つ縮小した形になります。

$$
A_{(4)} =
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & a_{44(4)}
\end{array}\right) \tag{16}
$$

\(A_{(3)} = l_3 u_3 + A_{(4)}\)より、\(A_{(1)} = l_1 u_1 + l_2 u_2 + l_3 u_3 + A_{(4)}\)という関係になります。

4段目

\(A_{(3)}\)から一つ縮小した\(A_{(4)}\)をさらにLU分解を行います。\(A_{(4)}\)の関係から\(LU\)行列も(17)式のようになります。

$$
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & a_{44(4)}
\end{array}\right)=
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & 1
\end{array}\right)
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{17}
$$

(17)式の関係から\(U\)行列の4行目は(18)式、\(L\)行列の4列目は(19)式のよう求まります。

$$
\begin{eqnarray}
a_{44(4)} &=& u_{44}
\end{eqnarray} \tag{18}
$$

$$
\begin{eqnarray}
a_{44(4)} = l_{44} u_{44} &\Rightarrow& l_{44} = a_{44(4)} / u_{44} = 1 \\
\end{eqnarray} \tag{19}
$$

ここで、\(l_4, u_4\)を(20)式のように定義して、

$$
l_4 = (0, 0 , 0 , 1) ^T , u_3 = (0, 0 , 0, u_{44}) \tag{20}
$$

\(A_{(4)} – l_4 u_4 = A_{(5)}\)とすると\(A_{(5)}\)は(21)式のようにゼロ行列になります。

$$
A_{(5)} =
\left(\begin{array}{cccc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{array}\right) \tag{21}
$$

\(A_{(5)}\)はゼロ行列なので分解できないので、ここでLUへの分解は終了します。

\(A_{(4)} = l_4 u_4\)より、\(A_{(1)} = l_1 u_1 + l_2 u_2 + l_3 u_3 + l_4 u_4\)という関係になります。

LU行列

これまでの計算で \(A_{(1)} = l_1 u_1 + l_2 u_2 + l_3 u_3 + l_4 u_4\)と分解が出来ました。それから\( A_{(1)} = A = LU\)という(1)式の関係がありましたので、(4)式、(10)式、(15)式、(20)式から、L行列とU行列は(22)式のように構成できます。

$$
\begin{eqnarray}
L &=& (l_1 , l_2 , l_3 , l_4 ) \\
U &=& (u_1 , u_2, u_3, u_4 ) ^T
\end{eqnarray} \tag{22}
$$

以上が基本的なA行列をLU行列へ分解するアルゴリズムです。1段目から4段目まで同じ手順で分解していきましたので、アルゴリズムとしてOKですね。

計算上の問題と対策

これまで代数で計算したので現れなかったのですが、\(L\)行列の要素を計算する(3)式、(9)式、(14)式、(19)式は\(U\)行列の対角要素、つまり\(A_{(i)}\)行列の左上の対角要素で割った計算を行っていました。

例えば、2段目の計算で(23)式のように\(A_{(2)}\)行列の左上の要素がゼロになってしまった場合、

$$
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{23(2)} & a_{24(2)} \\
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right)=
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 1 & 0 & 0 \\
0 & l_{32} & 1 & 0 \\
0 & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{23}
$$

\(u_{22} = 0 \)となり、L行列の要素計算で(24)式のように\(\infty\)となってしまい、計算エラーとなってしまいます。

$$
\begin{eqnarray}
a_{22(2)} = l_{22} u_{22} &\Rightarrow& l_{22} = a_{22(2)} / u_{22} = \infty \\
a_{32(2)} = l_{32} u_{22} &\Rightarrow& l_{32} = a_{32(2)} / u_{22} = \infty \\
a_{42(2)} = l_{42} u_{22} &\Rightarrow& l_{42} = a_{42(2)} / u_{22} = \infty
\end{eqnarray} \tag{24}
$$

この計算エラーを\(Pivot\)という方法で乗り越えます!

ここでようやく半分です。もう少しだけ頑張っていきましょう。

Pivot付きのLU分解

1段目

1段目は普通に計算し、\(l_1, u_1\)が求められたとします。

$$
\left(\begin{array}{cccc}
a_{11(1)} & a_{12 (1) } & a_{13 (1) } & a_{14 (1) } \\
a_{21(1)} & a_{22 (1) } & a_{23 (1) } & a_{24 (1) } \\
a_{31(1)} & a_{32 (1) } & a_{33 (1) } & a_{34 (1) } \\
a_{41(1)} & a_{42 (1) } & a_{43 (1) } & a_{44 (1) }
\end{array}\right) =
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{cccc}
u_{11} & u_{12} & u_{13} & u_{14} \\
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right)
$$

\(A_{(1)} – 1_l u_1 = A_{(2)}\)の計算で(25)式のように、左上の対角要素がゼロになったとします。

$$
A_{(2)} = \left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{23(2)} & a_{24(2)} \\
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right) \tag{25}
$$

とりあえずこの段階では、\(A_{(1)} = l_1u_1 + A_{(2)}\)です。

2段目

\(A_{(2)}\)の左上の対角要素がゼロとなっていますので、(26)式をそのまま計算すると、Lの要素が \(\infty\) となってしまいます。

$$
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{23(2)} & a_{24(2)} \\
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right) =
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 1 & 0 & 0 \\
0 & l_{32} & 1 & 0 \\
0 & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{26}
$$

そこでLの要素が\(\infty\)とならないように、値がゼロでない行と入替え(Pivot)を行います。例えば2行目と3行目をPivotするとします。

Pivot(行を入れ替える)方法は例えば2行目と3行目を入れ替える場合、(27)式の\(S_{2,3}\)を\(A_{(1)} = l_1 u_1 + A_{(2)}\)の左側から掛けることで、

$$
S_{2,3} =
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1
\end{array}\right) \tag{27}
$$

$$
\begin{eqnarray}
S_{2,3} A_{(1)} &=& S_{2,3} l_1 u_1 + S_{2,3} A_{(2)} \\
&=&
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1
\end{array}\right)
\left(\begin{array}{c}
1 \\ l_{21} \\ l_{31} \\ l_{41}
\end{array}\right)
\left(\begin{array}{cccc}
u_{11} & u_{12} & u_{13} & u_{14}
\end{array}\right) \\
&&+
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1
\end{array}\right)
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{23(2)} & a_{24(2)} \\
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right) \\
&=&
\left(\begin{array}{c}
1 \\ l_{31} \\ l_{21} \\ l_{41}
\end{array}\right)
\left(\begin{array}{cccc}
u_{11} & u_{12} & u_{13} & u_{14}
\end{array}\right)
+
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & 0 & a_{23(2)} & a_{24(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right) \tag{28}
\end{eqnarray}
$$

(28)式のように、\(A_{(2)}\)の左上の対角要素がゼロではなくなります。このとき、Pivotは関係する全体にかかるので、1段目で計算した\(l_1\)も(28)式のように入れ替わることになります。

このPivotを行った\(S_{2,3} A_{(2)}\)を以下の関係で分解し、\(l_2, u_2\)を計算します。

$$
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & a_{32(2)} & a_{33(2)} & a_{34(2)} \\
0 & 0 & a_{23(2)} & a_{24(2)} \\
0 & a_{42(2)} & a_{43(2)} & a_{44(2)}
\end{array}\right) =
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & 1 & 0 & 0 \\
0 & l_{32} & 1 & 0 \\
0 & l_{42} & l_{43} & 1
\end{array}\right)
\left(\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\ \hline
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right)
$$

すると\(S_{2,3}A_{(1)} = S_{2,3} l_1 u_1 + l_2 u_2 + A_{(3)}\)となります。\(l_2 u_2\)以降は\(S_{2,3}\)が掛かっていませんが、\(l_2, u_2 , A_{(3)}\)を計算するときに既にかかっているので内包していると考えます。

3段目

3段目の計算で\(A_{(3)}\)の左上の対角要素がゼロとなってしまったとします。

$$
A_{(3)} = \left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & a_{34(3)} \\
0 & 0 & a_{43(3)} & a_{44(3)}
\end{array}\right) \tag{28}
$$

ここで3行目と4行目を入れ替える\(S_{3,4}\)を左から掛けることで、

$$
\begin{eqnarray}
S_{3,4} S_{2,3} A &=& S_{3,4} S_{2,3}l_1 u_1 + S_{3,4} l_2 u_2 +
\left(\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 &1 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0
\end{array}\right)
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & a_{34(3)} \\
0 & 0 & a_{43(3)} & a_{44(3)}
\end{array}\right) \\
&=& S_{3,4} S_{2,3}l_1 u_1 + S_{3,4} l_2 u_2 +
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{43(3)} & a_{44(3)} \\
0 & 0 & 0 & a_{34(3)}
\end{array}\right) \tag{29}
\end{eqnarray}
$$

(29)式のように対角要素がゼロでない\(S_{3,4}A_{(3)}\)となります。また、このとき既に計算した\(l_1, l_2\)もPivotが行われます。

この\(S_{3,4}A_{(3)} \)で以下の関係で分解を行うことで、\(l_3, u_3\)を求めることが出来ます。

$$
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & a_{43(3)} & a_{44(3)} \\
0 & 0 & 0 & a_{34(3)}
\end{array}\right) =
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 1 & 0 \\
0 & 0 & l_{43} & 1
\end{array}\right)
\left(\begin{array}{cc|cc}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right)
$$

ここまでの関係は \(S_{3,4} S_{2,3} A = S_{3,4} S_{2,3}l_1 u_1 + S_{3,4} l_2 u_2 + l_3 u_3 + A_{(4)}\)です。

4段目

4段目は(30)式の関係で分解を行い、\(l_4, u_4\)が求まったとします。

$$
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & a_{44(4)}
\end{array}\right)=
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & 1
\end{array}\right)
\left(\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 0 & u_{44}
\end{array}\right) \tag{30}
$$

そうするとLUへの分解が終了し、Pivotを行った最終的な式は(31)式のようになります。

$$
S_{3,4}S_{2,4} A = S_{3,4} S_{2,3} l_1 u_1 + S_{3,4} l_2 u_2 + l_3 u_3 + l_4 u_4 \tag{31}
$$

LU分解後の処理

\(P = S_{3,4} S_{2,4}\)とし、\( L= (S_{3,4} S_{2,3} l_1 , S_{3,4} l_2 , l_3, l_4) \)、\(U = (u_1, u_2, u_3, u_4)^T \)とすると、Pivotを行ったLU分解は(32)式のようになります。

$$
PA = L U \tag{32}
$$

入替えの入替えは元に戻るという関係:\(P^2 = I\)より、\(P^{-1} = P\)となるので、(32)式を(33)式になります。

$$
A = P LU \tag{33}
$$

\(A \boldsymbol{x} = \boldsymbol{b}\)の関係があるので、これに(33)式を代入すると(34)式となり、

$$
P L U \boldsymbol{x} = \boldsymbol{b} \tag{34}
$$

(34)式は(35)式となるので、この関係から解\(\boldsymbol{x}\)を求めます。

$$
L U \boldsymbol{x} = P\boldsymbol{b} \tag{35}
$$

プログラムコード – Python

LU行列から方程式を求めるコードも再掲して、LU分解により方程式を求める全体のコードを掲載します。

ここに到達するまですごく長かったですが、プログラムコードは非常にスッキリしています。

Pivotの関数

import numpy as np

# Pivotの影響を受けるA, L P行列をk行を基準に入替え
def Pivot(A, L , k , P):
    n = A.shape[0]
    n_max = abs(A[k,k])
    pivot_n = k
    for i in range(k+1, n): # Pivotを行うため、最大の行を探索
        if n_max < A[i,k]:
            n_max = A[i,k]
            pivot_n = i
    
    # Pivot行列の入替え
    P_k = np.copy(P[k,:]) # '='だけでは参照渡しになるので copyを使用する
    P[k, :] = P[pivot_n, :]
    P[pivot_n, :] = P_k
    
    # A行列の入替え
    A_k = np.copy(A[k, :]) # '='だけでは参照渡しになるので copyを使用する
    A[k, :] = A[pivot_n, :]
    A[pivot_n, :] = A_k
    
    # L行列の入替え
    L_k = np.copy(L[k, :]) # '='だけでは参照渡しになるので copyを使用する
    L[k, :] = L[pivot_n, :]
    L[pivot_n, :] = L_k
    
    return A, L , P

LU分解の関数

import numpy as np

# A行列をLU分解
def LU_decomposition(A):
    # 配列の作成
    L = np.zeros(A.shape)
    U = np.zeros(A.shape)
    P = np.eye(A.shape[0], A.shape[1])

    for i in range(A.shape[0]):
        # Pivot実行
        A, L, P = Pivot(A, L , i, P)
        
        # L, U 行列の計算
        for j in range(i, A.shape[1]):
            U[i, j] = A[i, j]  # U行列 i 行目を算出
            L[j, i] = A[j ,i] / U[i, i] # L行列 i 列目を算出
    
        # 次のA行列を計算
        for j in range(A.shape[0]):
            for k in range(A.shape[1]):
                A[j,k] =A[j,k] - L[j, i] * U[i, k]

    return L, U, P

LU行列から方程式を解く関数

import numpy as np

# step 1 計算
def solve_z(L, b):
    z = np.zeros(b.shape)
    
    for i in range(z.shape[0]):
        r = 0
        for j in range(0, i):
            r += L[i,j] * z[j]
        z[i] = b[i] - r
    
    return z


# step 2 計算
def solve_x(U, z):
    x = np.zeros(z.shape)
    
    for i in reversed(range(x.shape[0])):
        r = 0
        for j in range(i + 1, x.shape[0]):
            r += U[i,j] * x[j]
        x[i] = (z[i] - r)/U[i,i]
        
    return x

実行

# A行列の定義
A = np.array([[1.0,2.0,3.0],
              [2.0,4.0,-9.0],
              [5.0,8.0,1.0]])

# A行列をLU行列へ分解、Pivot行列Pも出力
[L, U, P] = LU_decomposition(A)

# b行列の定義
b = np.array([1,2,3])

# b行列をP行列でPivot
pb = np.dot(P, b)
            
# 方程式の解を求める
z = solve_z(L, pb)
x = solve_x(U, z)

# 求めた解を出力
print(x) # 出力:[-1.  1.  0.]
世界最大のオンライン学習サイトで人気のPython講座で学んでみませんか?

コメント

タイトルとURLをコピーしました