この記事は,
旧ブログ から移植された記事です. よって, その内容として,
旧ブログ に依存した文脈が含まれている可能性があります. 予めご了承下さい.
LU 分解に関して初歩的な内容から網羅的にまとめた.
ガウスの消去法
ガウスの消去法は, 前進消去による上三角行列の形成と後退代入の組み合わせのことをいい, その本質は行列の行基本変形, すなわち中学校で習う連立方程式の解法そのものである. 例えば, 何の芸もないが, 次の連立方程式
{ x + 2 y = 3 3 x + 4 y = 5 ↔ ( 1 2 3 3 4 5 ) (1)
\begin{aligned}
\begin{cases}
x+2y&=&3 \\
3x+4y&=&5
\end{cases}\leftrightarrow
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
3 & 4 & 5
\end{array}\right)\tag{1}
\end{aligned}
{ x + 2 y 3 x + 4 y = = 3 5 ↔ ( 1 3 2 4 3 5 ) ( 1 )
をガウスの消去法では次のようにして解くのであった.
( 1 2 3 3 4 5 ) → ( 1 ⏞ × ( − 3 ) 2 ⏞ × ( − 3 ) 3 ⏞ × ( − 3 ) 3 ‾ 4 5 ) ⏟ 前進消去 → ( 1 2 3 0 − 2 − 4 ) → ( 1 2 ‾ 3 0 ⏞ × 1 − 2 ⏞ × 1 − 4 ⏞ × 1 ) ⏟ 後退代入 → ( 1 0 − 1 0 − 2 − 4 ) ∴ ( x , y ) = ( − 1 ÷ 1 , − 4 ÷ ( − 2 ) ) = ( − 1 , 2 )
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
3 & 4 & 5
\end{array}\right)\to
\underbrace{
\left(\begin{array}{cc|c}
\overbrace{1}^{\times (-3)} & \overbrace{2}^{\times (-3)} & \overbrace{3}^{\times (-3)} \\
\underline{3} & 4 & 5
\end{array}\right)
}_{\rm 前進消去}\to
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
0 & -2 & -4
\end{array}\right)\to
\underbrace{
\left(\begin{array}{cc|c}
1 & \underline{2} & 3 \\
\overbrace{0}^{\times 1} & \overbrace{-2}^{\times 1} & \overbrace{-4}^{\times 1}
\end{array}\right)
}_{\rm 後退代入} \\ \to
\left(\begin{array}{cc|c}
1 & 0 & -1 \\
0 & -2 & -4
\end{array}\right) \\
\therefore (x,y)=(-1\div 1,-4\div (-2))=(-1,2)
( 1 3 2 4 3 5 ) → 前進消去 1 × ( − 3 ) 3 2 × ( − 3 ) 4 3 × ( − 3 ) 5 → ( 1 0 2 − 2 3 − 4 ) → 後退代入 1 0 × 1 2 − 2 × 1 3 − 4 × 1 → ( 1 0 0 − 2 − 1 − 4 ) ∴ ( x , y ) = ( − 1 ÷ 1 , − 4 ÷ ( − 2 )) = ( − 1 , 2 )
この前進消去の操作は次のように行列の積で表現できる.
( 1 0 − 3 1 ) ( 1 2 3 4 ) = ( 1 2 0 − 2 )
\left(\begin{array}{cc}
1 & 0 \\
-3 & 1
\end{array}\right)
\left(\begin{array}{cc}
1 & 2 \\
3 & 4
\end{array}\right)=
\left(\begin{array}{cc}
1 & 2 \\
0 & -2
\end{array}\right)
( 1 − 3 0 1 ) ( 1 3 2 4 ) = ( 1 0 2 − 2 )
これを理解しておくと後述する LU 分解の理解に容易くなる. ここで, この時間計算量について, 一般の n n n X a = y w h e r e X ∈ R n × n , a ∈ R n × 1 , y ∈ R n × 1 X{\boldsymbol a}={\boldsymbol y}\ {\rm where}\ X\in\mathbb{R}^{n\times n}, {\boldsymbol a}\in\mathbb{R}^{n\times 1}, {\boldsymbol y}\in\mathbb{R}^{n\times 1} X a = y where X ∈ R n × n , a ∈ R n × 1 , y ∈ R n × 1
( x 11 x 12 x 13 ⋯ x 1 n x 21 x 22 x 23 ⋯ x 2 n x 31 x 32 x 33 ⋯ x 3 n ⋮ ⋮ ⋮ ⋱ ⋮ x n 1 x n 2 x n 3 ⋯ x n n ) ( a 1 a 2 a 3 ⋮ a n ) = ( y 1 y 2 y 3 ⋮ y n ) \left(\begin{array}{ccccc}
x_{11} & x_{12} & x_{13} & \cdots & x_{1n} \\
x_{21} & x_{22} & x_{23} & \cdots & x_{2n} \\
x_{31} & x_{32} & x_{33} & \cdots & x_{3n} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
x_{n1} & x_{n2} & x_{n3} & \cdots & x_{nn}
\end{array}\right)
\left(\begin{array}{c}a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_n\end{array}\right)=
\left(\begin{array}{c}y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n\end{array}\right)
x 11 x 21 x 31 ⋮ x n 1 x 12 x 22 x 32 ⋮ x n 2 x 13 x 23 x 33 ⋮ x n 3 ⋯ ⋯ ⋯ ⋱ ⋯ x 1 n x 2 n x 3 n ⋮ x nn a 1 a 2 a 3 ⋮ a n = y 1 y 2 y 3 ⋮ y n
まず前進消去の操作について考える. 前進消去では, x 11 x_{11} x 11 x 21 , x 31 , ⋯ , x n 1 x_{21}, x_{31}, \cdots, x_{n1} x 21 , x 31 , ⋯ , x n 1 x 21 , x 31 , ⋯ , x n 1 x_{21}, x_{31}, \cdots, x_{n1} x 21 , x 31 , ⋯ , x n 1 0 0 0 x 21 , x 31 , ⋯ , x n 1 x_{21}, x_{31}, \cdots, x_{n1} x 21 , x 31 , ⋯ , x n 1
x 22 x 23 ⋯ x 2 n x 32 x 33 ⋯ x 3 n ⋮ ⋮ ⋱ ⋮ x n 2 x n 3 ⋯ x n n
\begin{array}{cccc}
x_{22} & x_{23} & \cdots & x_{2n} \\
x_{32} & x_{33} & \cdots & x_{3n} \\
\vdots & \vdots & \ddots & \vdots \\
x_{n2} & x_{n3} & \cdots & x_{nn}
\end{array}
x 22 x 32 ⋮ x n 2 x 23 x 33 ⋮ x n 3 ⋯ ⋯ ⋱ ⋯ x 2 n x 3 n ⋮ x nn
の箇所であることがいえる. この部分の成分数は ( n − 1 ) 2 (n-1)^2 ( n − 1 ) 2 ( n − 1 ) 2 (n-1)^2 ( n − 1 ) 2 x 22 x_{22} x 22 x 32 , ⋯ , x n 2 x_{32},\cdots,x_{n2} x 32 , ⋯ , x n 2 x 32 , ⋯ , x n 2 x_{32},\cdots,x_{n2} x 32 , ⋯ , x n 2 0 0 0
x 33 ⋯ x 3 n ⋮ ⋱ ⋮ x n 3 ⋯ x n n
\begin{array}{ccc}
x_{33} & \cdots & x_{3n} \\
\vdots & \ddots & \vdots \\
x_{n3} & \cdots & x_{nn}
\end{array}
x 33 ⋮ x n 3 ⋯ ⋱ ⋯ x 3 n ⋮ x nn
の箇所である. この部分の成分数は ( n − 2 ) 2 (n-2)^2 ( n − 2 ) 2 ( n − 2 ) 2 (n-2)^2 ( n − 2 ) 2
( n − 1 ) 2 + ( n − 2 ) 2 + ⋯ + 2 2 + 1 2 = ∑ i = 1 n ( i − 1 ) 2 (n-1)^2+(n-2)^2+\cdots+2^2+1^2
=\sum^{n}_{i=1}(i-1)^2 ( n − 1 ) 2 + ( n − 2 ) 2 + ⋯ + 2 2 + 1 2 = i = 1 ∑ n ( i − 1 ) 2
ここで, 高校数学の教科書にも載っている定理(証明略) ∑ j = 1 n j 2 = n ( n + 1 ) ( 2 n + 1 ) 6 \sum^n_{j=1}j^2=\frac{n(n+1)(2n+1)}{6} ∑ j = 1 n j 2 = 6 n ( n + 1 ) ( 2 n + 1 ) ∑ i = 1 n ( i − 1 ) 2 = n 3 3 + n 2 \sum^n_{i=1}(i-1)^2=\frac{n^3}{3}+n^2 ∑ i = 1 n ( i − 1 ) 2 = 3 n 3 + n 2
続いて, 後退代入について考える. 後退代入は, 上三角行列となっている係数行列に対して代入を繰り返し, 対角行列を形成していく操作であった.
( x 11 x 12 x 13 ⋯ x 1 n x 22 x 23 ⋯ x 2 n x 33 ⋯ x 3 n ⋱ ⋮ x n n ) ( a 1 a 2 a 3 ⋮ a n ) = ( y 1 y 2 y 3 ⋮ y n ) \left(\begin{array}{ccccc}
x_{11} & x_{12} & x_{13} & \cdots & x_{1n} \\
& x_{22} & x_{23} & \cdots & x_{2n} \\
& & x_{33} & \cdots & x_{3n} \\
& & & \ddots & \vdots \\
&&&& x_{nn}
\end{array}\right)
\left(\begin{array}{c}a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_n\end{array}\right)=
\left(\begin{array}{c}y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n\end{array}\right)
x 11 x 12 x 22 x 13 x 23 x 33 ⋯ ⋯ ⋯ ⋱ x 1 n x 2 n x 3 n ⋮ x nn a 1 a 2 a 3 ⋮ a n = y 1 y 2 y 3 ⋮ y n
この場合, まず x n n x_{nn} x nn x 1 n , x 2 n , x 3 n , ⋯ , x ( n − 1 ) n x_{1n},x_{2n},x_{3n},\cdots,x_{(n-1)n} x 1 n , x 2 n , x 3 n , ⋯ , x ( n − 1 ) n 0 0 0 x 22 x_{22} x 22 0 0 0 y n y_n y n y n − 1 y_{n-1} y n − 1 ( y n − 1 ) ′ = y n − 1 − x ( n − 1 ) n x n n y n (y_{n-1})'=y_{n-1}-\frac{x_{(n-1)n}}{x_{nn}}y_n ( y n − 1 ) ′ = y n − 1 − x nn x ( n − 1 ) n y n y n − 2 y_{n-2} y n − 2 y n − 2 − x ( n − 2 ) n x ( n − 1 ) n ( y n − 1 ) ′ y_{n-2}-\frac{x_{(n-2)n}}{x_{(n-1)n}}(y_{n-1})' y n − 2 − x ( n − 1 ) n x ( n − 2 ) n ( y n − 1 ) ′ ( y n − 1 ) ′ (y_{n-1})' ( y n − 1 ) ′ y 1 y_1 y 1
1 + 2 + ⋯ + ( n − 1 ) = n ( n − 1 ) 2 1+2+\cdots+(n-1)=\frac{n(n-1)}{2} 1 + 2 + ⋯ + ( n − 1 ) = 2 n ( n − 1 )
とかける(右辺への式変形の証明は高校数学の教科書で扱われているので略). ガウスの消去法は前進消去と後退代入の組み合わせであるので, その時間計算量は,
O ( n 3 3 + n ( n − 1 ) 2 ) (2) O(\frac{n^3}{3}+\frac{n(n-1)}{2})\tag{2} O ( 3 n 3 + 2 n ( n − 1 ) ) ( 2 )
ただし各項の中で最高次の係数だけを考えればよいので, ガウスの消去法の時間計算量は
1 3 O ( n 3 ) \frac{1}{3}O(n^3) 3 1 O ( n 3 )
である.
ガウス・ジョルダン法
ガウスの名がつく行列を用いた線形方程式の直接解法には, 今述べたガウスの消去法のほかにガウス・ジョルダン法というものがある. これらは明確に区別されたところであまり意味がないような気もするが, ガウスの消去法が上記のように係数行列となる部分を単位行列でない対角行列へと変形していったのに対し, ガウス・ジョルダン法はそれを直接単位行列となるように変形していく点で異なる. ガウス・ジョルダン法で同様にして計算していくと,
( 1 2 3 3 4 5 ) → ( 1 ⏞ × ( − 3 ) 2 ⏞ × ( − 3 ) 3 ⏞ × ( − 3 ) 3 ‾ 4 5 ) → ( 1 2 3 0 − 2 − 4 ) → ( 1 2 3 0 ⏞ × ( − 1 2 ) − 2 ⏞ × ( − 1 2 ) − 4 ⏞ × ( − 1 2 ) ) → ( 1 2 3 0 1 2 ) → ( 1 2 ‾ 3 0 ⏞ × ( − 2 ) 1 ⏞ × ( − 2 ) 2 ⏞ × ( − 2 ) ) → ( 1 0 − 1 0 1 2 ) ∴ ( x , y ) = ( − 1 , 2 )
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
3 & 4 & 5
\end{array}\right)\to
\left(\begin{array}{cc|c}
\overbrace{1}^{\times (-3)} & \overbrace{2}^{\times (-3)} & \overbrace{3}^{\times (-3)} \\
\underline{3} & 4 & 5
\end{array}\right)\to
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
0 & -2 & -4
\end{array}\right)\to
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
\overbrace{0}^{\times (-\frac{1}{2})} & \overbrace{-2}^{\times (-\frac{1}{2})} & \overbrace{-4}^{\times (-\frac{1}{2})}
\end{array}\right) \\ \to
\left(\begin{array}{cc|c}
1 & 2 & 3 \\
0 & 1 & 2
\end{array}\right)\to
\left(\begin{array}{cc|c}
1 & \underline{2} & 3 \\
\overbrace{0}^{\times (-2)} & \overbrace{1}^{\times (-2)} & \overbrace{2}^{\times (-2)}
\end{array}\right)\to
\left(\begin{array}{cc|c}
1 & 0 & -1 \\
0 & 1 & 2
\end{array}\right) \\
\therefore (x,y)=(-1,2)
( 1 3 2 4 3 5 ) → 1 × ( − 3 ) 3 2 × ( − 3 ) 4 3 × ( − 3 ) 5 → ( 1 0 2 − 2 3 − 4 ) → 1 0 × ( − 2 1 ) 2 − 2 × ( − 2 1 ) 3 − 4 × ( − 2 1 ) → ( 1 0 2 1 3 2 ) → 1 0 × ( − 2 ) 2 1 × ( − 2 ) 3 2 × ( − 2 ) → ( 1 0 0 1 − 1 2 ) ∴ ( x , y ) = ( − 1 , 2 )
これは掃き出し法とも言われることがある. ガウス・ジョルダン法の計算量は 1 2 O ( n 3 ) \frac{1}{2}O(n^3) 2 1 O ( n 3 )
ピボッティング
ところで, いま示したガウスの消去法, ガウスジョルダン法の手順は, このままでは困る場合がある. 例えば, 次の線形方程式をガウスの消去法で解いてみると
( 1 2 7 6 6 2 4 4 2 2 1 8 5 2 12 2 4 3 3 5 ) → ( 1 ⏞ × ( − 2 ) 2 ⏞ × ( − 2 ) 7 ⏞ × ( − 2 ) 6 ⏞ × ( − 2 ) 6 ⏞ × ( − 2 ) 2 ‾ 4 4 2 2 1 8 5 2 12 2 4 3 3 5 ) → ( 1 2 7 6 6 0 0 − 10 − 10 − 10 1 8 5 2 12 2 4 3 3 5 ) → ( 1 ⏞ × ( − 1 ) 2 ⏞ × ( − 1 ) 7 ⏞ × ( − 1 ) 6 ⏞ × ( − 1 ) 6 ⏞ × ( − 1 ) 0 0 − 10 − 10 − 10 1 ‾ 8 5 2 12 2 4 3 3 5 ) → ( 1 2 7 6 6 0 0 − 10 − 10 − 10 0 6 − 2 − 4 6 2 4 3 3 5 ) → ( 1 ⏞ × ( − 2 ) 2 ⏞ × ( − 2 ) 7 ⏞ × ( − 2 ) 6 ⏞ × ( − 2 ) 6 ⏞ × ( − 2 ) 0 0 − 10 − 10 − 10 0 6 − 2 − 4 6 2 ‾ 4 3 3 5 ) → ( 1 2 7 6 6 0 0 − 10 − 10 − 10 0 6 − 2 − 4 6 0 0 − 11 − 9 − 7 ) → ( 1 2 7 6 6 0 0 ⏞ × ( − 6 0 ) − 10 ⏞ × ( − 6 0 ) − 10 ⏞ × ( − 6 0 ) − 10 ⏞ × ( − 6 0 ) 0 6 ‾ − 2 − 4 6 0 0 − 11 − 9 − 7 )
\left(\begin{array}{cccc|c}
1 & 2 & 7 & 6 & 6 \\
2 & 4 & 4 & 2 & 2 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
\overbrace{1}^{\times (-2)} & \overbrace{2}^{\times (-2)} & \overbrace{7}^{\times (-2)} & \overbrace{6}^{\times (-2)} & \overbrace{6}^{\times (-2)} \\
\underline{2} & 4 & 4 & 2 & 2 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
1 & 2 & 7 & 6 & 6 \\
0 & 0 & -10 & -10 & -10 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right) \\ \to
\left(\begin{array}{cccc|c}
\overbrace{1}^{\times (-1)} & \overbrace{2}^{\times (-1)} & \overbrace{7}^{\times (-1)} & \overbrace{6}^{\times (-1)} & \overbrace{6}^{\times (-1)} \\
0 & 0 & -10 & -10 & -10 \\
\underline{1} & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
1 & 2 & 7 & 6 & 6 \\
0 & 0 & -10 & -10 & -10 \\
0 & 6 & -2 & -4 & 6 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
\overbrace{1}^{\times (-2)} & \overbrace{2}^{\times (-2)} & \overbrace{7}^{\times (-2)} & \overbrace{6}^{\times (-2)} & \overbrace{6}^{\times (-2)} \\
0 & 0 & -10 & -10 & -10 \\
0 & 6 & -2 & -4 & 6 \\
\underline{2} & 4 & 3 & 3 & 5
\end{array}\right) \\ \to
\left(\begin{array}{cccc|c}
1 & 2 & 7 & 6 & 6 \\
0 & 0 & -10 & -10 & -10 \\
0 & 6 & -2 & -4 & 6 \\
0 & 0 & -11 & -9 & -7
\end{array}\right)\to
\left(\begin{array}{cccc|c}
1 & 2 & 7 & 6 & 6 \\
0 & \overbrace{0}^{\times (-\frac{6}{0})} & \overbrace{-10}^{\times (-\frac{6}{0})} & \overbrace{-10}^{\times (-\frac{6}{0})} & \overbrace{-10}^{\times (-\frac{6}{0})} \\
0 & \underline{6} & -2 & -4 & 6 \\
0 & 0 & -11 & -9 & -7
\end{array}\right)
1 2 1 2 2 4 8 4 7 4 5 3 6 2 2 3 6 2 12 5 → 1 × ( − 2 ) 2 1 2 2 × ( − 2 ) 4 8 4 7 × ( − 2 ) 4 5 3 6 × ( − 2 ) 2 2 3 6 × ( − 2 ) 2 12 5 → 1 0 1 2 2 0 8 4 7 − 10 5 3 6 − 10 2 3 6 − 10 12 5 → 1 × ( − 1 ) 0 1 2 2 × ( − 1 ) 0 8 4 7 × ( − 1 ) − 10 5 3 6 × ( − 1 ) − 10 2 3 6 × ( − 1 ) − 10 12 5 → 1 0 0 2 2 0 6 4 7 − 10 − 2 3 6 − 10 − 4 3 6 − 10 6 5 → 1 × ( − 2 ) 0 0 2 2 × ( − 2 ) 0 6 4 7 × ( − 2 ) − 10 − 2 3 6 × ( − 2 ) − 10 − 4 3 6 × ( − 2 ) − 10 6 5 → 1 0 0 0 2 0 6 0 7 − 10 − 2 − 11 6 − 10 − 4 − 9 6 − 10 6 − 7 → 1 0 0 0 2 0 × ( − 0 6 ) 6 0 7 − 10 × ( − 0 6 ) − 2 − 11 6 − 10 × ( − 0 6 ) − 4 − 9 6 − 10 × ( − 0 6 ) 6 − 7
というように 0 0 0
というわけで, 下記に部分ピボット選択付きガウスの消去法による計算過程を省略することなく書き下したが, とくに面白みもない上に長いので隠しておいた. クリックで展開.
( 1 2 7 6 6 2 4 4 2 2 1 8 5 2 12 2 4 3 3 5 ) → ( 2 4 4 2 2 1 2 7 6 6 1 8 5 2 12 2 4 3 3 5 ) → ( 2 ⏞ × ( − 1 2 ) 4 ⏞ × ( − 1 2 ) 4 ⏞ × ( − 1 2 ) 2 ⏞ × ( − 1 2 ) 2 ⏞ × ( − 1 2 ) 1 ‾ 2 7 6 6 1 8 5 2 12 2 4 3 3 5 ) → ( 2 4 4 2 2 0 0 5 5 5 1 8 5 2 12 2 4 3 3 5 ) → ( 2 ⏞ × ( − 1 2 ) 4 ⏞ × ( − 1 2 ) 4 ⏞ × ( − 1 2 ) 2 ⏞ × ( − 1 2 ) 2 ⏞ × ( − 1 2 ) 0 0 5 5 5 1 ‾ 8 5 2 12 2 4 3 3 5 ) → ( 2 4 4 2 2 0 0 5 5 5 0 6 3 1 11 2 4 3 3 5 ) → ( 2 ⏞ × ( − 1 ) 4 ⏞ × ( − 1 ) 4 ⏞ × ( − 1 ) 2 ⏞ × ( − 1 ) 2 ⏞ × ( − 1 ) 0 0 5 5 5 0 6 3 1 11 2 ‾ 4 3 3 5 ) → ( 2 4 4 2 2 0 0 5 5 5 0 6 3 1 11 0 0 − 1 1 3 ) → ( 2 4 4 2 2 0 6 3 1 11 0 0 5 5 5 0 0 − 1 1 3 ) → ( 2 4 4 2 2 0 6 3 1 11 0 0 5 5 5 0 0 − 1 1 3 ) → ( 2 4 4 2 2 0 6 3 1 11 0 0 5 ⏞ × 1 5 5 ⏞ × 1 5 5 ⏞ × 1 5 0 0 − 1 ‾ 1 3 ) → ( 2 4 4 2 2 0 6 3 1 11 0 0 5 5 5 0 0 0 2 4 ) ⏟ 上三角行列 → ( 2 4 4 2 ‾ 2 0 6 3 1 11 0 0 5 5 5 0 0 0 2 ⏞ × ( − 1 ) 4 ⏞ × ( − 1 ) ) → ( 2 4 4 0 − 2 0 6 3 1 11 0 0 5 5 5 0 0 0 2 4 ) → ( 2 4 4 0 − 2 0 6 3 1 ‾ 11 0 0 5 5 5 0 0 0 2 ⏞ × ( − 1 2 ) 4 ⏞ × ( − 1 2 ) ) → ( 2 4 4 0 − 2 0 6 3 0 9 0 0 5 5 5 0 0 0 2 4 ) → ( 2 4 4 0 − 2 0 6 3 0 9 0 0 5 5 ‾ 5 0 0 0 2 ⏞ × ( − 5 2 ) 4 ⏞ × ( − 5 2 ) ) → ( 2 4 4 0 − 2 0 6 3 0 9 0 0 5 0 − 5 0 0 0 2 4 ) → ( 2 4 4 ‾ 0 − 2 0 6 3 0 9 0 0 5 ⏞ × ( − 4 5 ) 0 ⏞ × ( − 4 5 ) − 5 ⏞ × ( − 4 5 ) 0 0 0 2 4 ) → ( 2 4 0 0 2 0 6 3 0 9 0 0 5 0 − 5 0 0 0 2 4 ) → ( 2 4 0 0 2 0 6 3 ‾ 0 9 0 0 5 ⏞ × ( − 3 5 ) 0 ⏞ × ( − 3 5 ) − 5 ⏞ × ( − 3 5 ) 0 0 0 2 4 ) → ( 2 4 0 0 2 0 6 0 0 12 0 0 5 0 − 5 0 0 0 2 4 ) → ( 2 4 ‾ 0 0 2 0 6 ⏞ × ( − 2 3 ) 0 ⏞ × ( − 2 3 ) 0 ⏞ × ( − 2 3 ) 12 ⏞ × ( − 2 3 ) 0 0 5 0 − 5 0 0 0 2 4 ) → ( 2 0 0 0 − 6 0 6 0 0 12 0 0 5 0 − 5 0 0 0 2 4 ) ∴ ( − 6 ÷ 2 , 12 ÷ 6 , − 5 ÷ 0 , 4 ÷ 2 ) T = ( − 3 , 2 , − 1 , 2 ) T
\left(\begin{array}{cccc|c}
1 & 2 & 7 & 6 & 6 \\
2 & 4 & 4 & 2 & 2 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
1 & 2 & 7 & 6 & 6 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
\overbrace{2}^{\times(-\frac{1}{2})} & \overbrace{4}^{\times(-\frac{1}{2})} & \overbrace{4}^{\times(-\frac{1}{2})} & \overbrace{2}^{\times(-\frac{1}{2})} & \overbrace{2}^{\times(-\frac{1}{2})} \\
\underline{1} & 2 & 7 & 6 & 6 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\\ \to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 0 & 5 & 5 & 5 \\
1 & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
\overbrace{2}^{\times(-\frac{1}{2})} & \overbrace{4}^{\times(-\frac{1}{2})} & \overbrace{4}^{\times(-\frac{1}{2})} & \overbrace{2}^{\times(-\frac{1}{2})} & \overbrace{2}^{\times(-\frac{1}{2})} \\
0 & 0 & 5 & 5 & 5 \\
\underline{1} & 8 & 5 & 2 & 12 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 0 & 5 & 5 & 5 \\
0 & 6 & 3 & 1 & 11 \\
2 & 4 & 3 & 3 & 5
\end{array}\right)\\ \to
\left(\begin{array}{cccc|c}
\overbrace{2}^{\times (-1)} & \overbrace{4}^{\times (-1)} & \overbrace{4}^{\times (-1)} & \overbrace{2}^{\times (-1)} & \overbrace{2}^{\times (-1)} \\
0 & 0 & 5 & 5 & 5 \\
0 & 6 & 3 & 1 & 11 \\
\underline{2} & 4 & 3 & 3 & 5
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 0 & 5 & 5 & 5 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & -1 & 1 & 3
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & -1 & 1 & 3
\end{array}\right) \\ \to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & -1 & 1 & 3
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & \overbrace{5}^{\times\frac{1}{5}} & \overbrace{5}^{\times\frac{1}{5}} & \overbrace{5}^{\times\frac{1}{5}} \\
0 & 0 & \underline{-1} & 1 & 3
\end{array}\right)\to
\underbrace{\left(\begin{array}{cccc|c}
2 & 4 & 4 & 2 & 2 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right)}_{上三角行列} \\ \to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & \underline{2} & 2 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & 0 & \overbrace{2}^{\times(-1)} & \overbrace{4}^{\times(-1)}
\end{array}\right) \to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 0 & -2 \\
0 & 6 & 3 & 1 & 11 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 0 & -2 \\
0 & 6 & 3 & \underline{1} & 11 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & 0 & \overbrace{2}^{\times (-\frac{1}{2})} & \overbrace{4}^{\times (-\frac{1}{2})}
\end{array}\right) \\ \to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 0 & -2 \\
0 & 6 & 3 & 0 & 9 \\
0 & 0 & 5 & 5 & 5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 0 & -2 \\
0 & 6 & 3 & 0 & 9 \\
0 & 0 & 5 & \underline{5} & 5 \\
0 & 0 & 0 & \overbrace{2}^{\times (-\frac{5}{2})} & \overbrace{4}^{\times (-\frac{5}{2})}
\end{array}\right)\to
\left(\begin{array}{cccc|c}
2 & 4 & 4 & 0 & -2 \\
0 & 6 & 3 & 0 & 9 \\
0 & 0 & 5 & 0 & -5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \\ \to
\left(\begin{array}{cccc|c}
2 & 4 & \underline{4} & 0 & -2 \\
0 & 6 & 3 & 0 & 9 \\
0 & 0 & \overbrace{5}^{\times (-\frac{4}{5})} & \overbrace{0}^{\times (-\frac{4}{5})} & \overbrace{-5}^{\times (-\frac{4}{5})} \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \to
\left(\begin{array}{cccc|c}
2 & 4 & 0 & 0 & 2 \\
0 & 6 & 3 & 0 & 9 \\
0 & 0 & 5 & 0 & -5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \to
\left(\begin{array}{cccc|c}
2 & 4 & 0 & 0 & 2 \\
0 & 6 & \underline{3} & 0 & 9 \\
0 & 0 & \overbrace{5}^{\times (-\frac{3}{5})} & \overbrace{0}^{\times (-\frac{3}{5})} & \overbrace{-5}^{\times (-\frac{3}{5})} \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \\ \to
\left(\begin{array}{cccc|c}
2 & 4 & 0 & 0 & 2 \\
0 & 6 & 0 & 0 & 12 \\
0 & 0 & 5 & 0 & -5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \to
\left(\begin{array}{cccc|c}
2 & \underline{4} & 0 & 0 & 2 \\
0 & \overbrace{6}^{\times(-\frac{2}{3})} & \overbrace{0}^{\times(-\frac{2}{3})} & \overbrace{0}^{\times(-\frac{2}{3})} & \overbrace{12}^{\times(-\frac{2}{3})} \\
0 & 0 & 5 & 0 & -5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \to
\left(\begin{array}{cccc|c}
2 & 0 & 0 & 0 & -6 \\
0 & 6 & 0 & 0 & 12 \\
0 & 0 & 5 & 0 & -5 \\
0 & 0 & 0 & 2 & 4
\end{array}\right) \\
\therefore (-6\div 2,12\div 6, -5\div 0, 4\div 2)^T =(-3,2,-1,2)^T
1 2 1 2 2 4 8 4 7 4 5 3 6 2 2 3 6 2 12 5 → 2 1 1 2 4 2 8 4 4 7 5 3 2 6 2 3 2 6 12 5 → 2 × ( − 2 1 ) 1 1 2 4 × ( − 2 1 ) 2 8 4 4 × ( − 2 1 ) 7 5 3 2 × ( − 2 1 ) 6 2 3 2 × ( − 2 1 ) 6 12 5 → 2 0 1 2 4 0 8 4 4 5 5 3 2 5 2 3 2 5 12 5 → 2 × ( − 2 1 ) 0 1 2 4 × ( − 2 1 ) 0 8 4 4 × ( − 2 1 ) 5 5 3 2 × ( − 2 1 ) 5 2 3 2 × ( − 2 1 ) 5 12 5 → 2 0 0 2 4 0 6 4 4 5 3 3 2 5 1 3 2 5 11 5 → 2 × ( − 1 ) 0 0 2 4 × ( − 1 ) 0 6 4 4 × ( − 1 ) 5 3 3 2 × ( − 1 ) 5 1 3 2 × ( − 1 ) 5 11 5 → 2 0 0 0 4 0 6 0 4 5 3 − 1 2 5 1 1 2 5 11 3 → 2 0 0 0 4 6 0 0 4 3 5 − 1 2 1 5 1 2 11 5 3 → 2 0 0 0 4 6 0 0 4 3 5 − 1 2 1 5 1 2 11 5 3 → 2 0 0 0 4 6 0 0 4 3 5 × 5 1 − 1 2 1 5 × 5 1 1 2 11 5 × 5 1 3 → 上三角行列 2 0 0 0 4 6 0 0 4 3 5 0 2 1 5 2 2 11 5 4 → 2 0 0 0 4 6 0 0 4 3 5 0 2 1 5 2 × ( − 1 ) 2 11 5 4 × ( − 1 ) → 2 0 0 0 4 6 0 0 4 3 5 0 0 1 5 2 − 2 11 5 4 → 2 0 0 0 4 6 0 0 4 3 5 0 0 1 5 2 × ( − 2 1 ) − 2 11 5 4 × ( − 2 1 ) → 2 0 0 0 4 6 0 0 4 3 5 0 0 0 5 2 − 2 9 5 4 → 2 0 0 0 4 6 0 0 4 3 5 0 0 0 5 2 × ( − 2 5 ) − 2 9 5 4 × ( − 2 5 ) → 2 0 0 0 4 6 0 0 4 3 5 0 0 0 0 2 − 2 9 − 5 4 → 2 0 0 0 4 6 0 0 4 3 5 × ( − 5 4 ) 0 0 0 0 × ( − 5 4 ) 2 − 2 9 − 5 × ( − 5 4 ) 4 → 2 0 0 0 4 6 0 0 0 3 5 0 0 0 0 2 2 9 − 5 4 → 2 0 0 0 4 6 0 0 0 3 5 × ( − 5 3 ) 0 0 0 0 × ( − 5 3 ) 2 2 9 − 5 × ( − 5 3 ) 4 → 2 0 0 0 4 6 0 0 0 0 5 0 0 0 0 2 2 12 − 5 4 → 2 0 0 0 4 6 × ( − 3 2 ) 0 0 0 0 × ( − 3 2 ) 5 0 0 0 × ( − 3 2 ) 0 2 2 12 × ( − 3 2 ) − 5 4 → 2 0 0 0 0 6 0 0 0 0 5 0 0 0 0 2 − 6 12 − 5 4 ∴ ( − 6 ÷ 2 , 12 ÷ 6 , − 5 ÷ 0 , 4 ÷ 2 ) T = ( − 3 , 2 , − 1 , 2 ) T
LU 分解
漸く本題の LU 分解(LR 分解, 三角分解)について. 簡単のために式 ( 1 ) (1) ( 1 ) ( 1 ) (1) ( 1 )
A ( 0 ) ( x y ) = v w h e r e A ( 0 ) = ( 1 2 3 4 ) , v = ( 3 5 )
A^{(0)}\left(\begin{array}{c}x \\ y\end{array}\right)={\boldsymbol v}
\ {\rm where}\ A^{(0)}=\left(\begin{array}{cc}1 & 2 \\ 3 & 4\end{array}\right),{\boldsymbol v}=\left(\begin{array}{c}3\\ 5\end{array}\right)
A ( 0 ) ( x y ) = v where A ( 0 ) = ( 1 3 2 4 ) , v = ( 3 5 )
A ( 0 ) − 1 v {A^{(0)}}^{-1}{\boldsymbol v} A ( 0 ) − 1 v ( x , y ) T (x,y)^T ( x , y ) T 1 3 O ( n 3 ) \frac{1}{3}O(n^3) 3 1 O ( n 3 ) A ( 0 ) A^{(0)} A ( 0 )
A ( 1 ) = L ( 1 ) A ( 0 ) = ( 1 2 0 − 2 ) ↔ A ( 0 ) = L ( 1 ) − 1 A ( 1 ) w h e r e L ( 1 ) = ( 1 0 − 3 1 )
A^{(1)}=L^{(1)}A^{(0)}=\left(\begin{array}{cc}1 & 2 \\ 0 & -2\end{array}\right)
\leftrightarrow
A^{(0)}={L^{(1)}}^{-1}A^{(1)}
\ {\rm where}\ L^{(1)}=\left(\begin{array}{cc}1 & 0 \\ -3 & 1\end{array}\right)
A ( 1 ) = L ( 1 ) A ( 0 ) = ( 1 0 2 − 2 ) ↔ A ( 0 ) = L ( 1 ) − 1 A ( 1 ) where L ( 1 ) = ( 1 − 3 0 1 )
従って
L ( 1 ) − 1 A ( 1 ) ( x y ) = v
{L^{(1)}}^{-1}A^{(1)}\left(\begin{array}{c}x \\ y\end{array}\right)={\boldsymbol v}
L ( 1 ) − 1 A ( 1 ) ( x y ) = v
ここで, b = A ( 1 ) ( x , y ) T {\boldsymbol b}=A^{(1)}(x,y)^T b = A ( 1 ) ( x , y ) T L ( 1 ) − 1 b = v {L^{(1)}}^{-1}{\boldsymbol b}={\boldsymbol v} L ( 1 ) − 1 b = v b {\boldsymbol b} b L ( 1 ) L^{(1)} L ( 1 ) L ( 1 ) L^{(1)} L ( 1 ) 1 2 O ( n 2 ) ∵ \frac{1}{2}O(n^2)\ \because 2 1 O ( n 2 ) ∵
その後に A ( 1 ) ( x , y ) T = b A^{(1)}(x,y)^T={\boldsymbol b} A ( 1 ) ( x , y ) T = b ( x , y ) T (x,y)^T ( x , y ) T A ( 1 ) A^{(1)} A ( 1 ) 1 2 O ( n 2 ) ∵ \frac{1}{2}O(n^2)\ \because 2 1 O ( n 2 ) ∵
よって, この一連の操作における時間計算量は 1 3 O ( n 3 ) \frac{1}{3}O(n^3) 3 1 O ( n 3 )
しかし, L U L U LU A ( 0 ) A^{(0)} A ( 0 ) v {\boldsymbol v} v 1 2 O ( n 2 ) \frac{1}{2}O(n^2) 2 1 O ( n 2 )
いまの説明では, 式 \(\) において A ( 0 ) A^{(0)} A ( 0 ) L ( 1 ) − 1 {L^{(1)}}^{-1} L ( 1 ) − 1 A ( 1 ) A^{(1)} A ( 1 ) L = L ( 1 ) − 1 , U = A ( 1 ) L={L^{(1)}}^{-1},U=A^{(1)} L = L ( 1 ) − 1 , U = A ( 1 ) L ( i ) L^{(i)} L ( i )
すべての前進消去の行列 L ( i ) L^{(i)} L ( i ) A ∈ R m × n A\in\mathbb{R}^{m\times n} A ∈ R m × n A = A ( 0 ) = L U w h e r e L = ∏ i = 1 n − 1 L ( i ) − 1 , U = A ( n − 1 ) A=A^{(0)}=L U\ {\rm where}\ L=\prod_{i=1}^{n-1}{L^{(i)}}^{-1}, U=A^{(n-1)} A = A ( 0 ) = LU where L = i = 1 ∏ n − 1 L ( i ) − 1 , U = A ( n − 1 )
補足すると, A ( 0 ) A^{(0)} A ( 0 ) n − 1 n-1 n − 1 L ( n − 1 ) ⋯ L ( 2 ) L ( 1 ) A ( 0 ) L^{(n-1)}\cdots L^{(2)}L^{(1)}A^{(0)} L ( n − 1 ) ⋯ L ( 2 ) L ( 1 ) A ( 0 ) λ = L ( n − 1 ) ⋯ L ( 2 ) L ( 1 ) \lambda=L^{(n-1)}\cdots L^{(2)}L^{(1)} λ = L ( n − 1 ) ⋯ L ( 2 ) L ( 1 ) λ A ( 0 ) = U = A ( n − 1 ) \lambda A^{(0)}=U=A^{(n-1)} λ A ( 0 ) = U = A ( n − 1 ) A ( 0 ) = λ − 1 U A^{(0)}=\lambda^{-1}U A ( 0 ) = λ − 1 U ( X Y ) − 1 = Y − 1 X − 1 (X Y)^{-1}=Y^{-1} X^{-1} ( X Y ) − 1 = Y − 1 X − 1
A ( 0 ) = ( 3 1 0 6 1 − 2 − 3 0 3 ) → ( 1 0 0 − 2 1 0 1 0 1 ) ⏟ L ( 1 ) ( 3 1 0 6 1 − 2 − 3 0 3 ) = ( 3 1 0 0 − 1 − 2 0 1 3 ) ⏟ A ( 1 ) → ( 1 0 0 0 1 0 0 1 1 ) ⏟ L ( 2 ) A ( 1 ) = ( 3 1 0 0 − 1 − 2 0 0 1 ) ⏟ A ( 2 )
A^{(0)}=
\left(\begin{array}{ccc}
3 & 1 & 0 \\
6 & 1 & -2 \\
-3 & 0 & 3
\end{array}\right)\to
\underbrace{\left(\begin{array}{ccc}
1 & 0 & 0 \\
-2 & 1 & 0 \\
1 & 0 & 1
\end{array}\right)}_{L^{(1)}}
\left(\begin{array}{ccc}
3 & 1 & 0 \\
6 & 1 & -2 \\
-3 & 0 & 3
\end{array}\right)=
\underbrace{
\left(\begin{array}{ccc}
3 & 1 & 0 \\
0 & -1 & -2 \\
0 & 1 & 3
\end{array}\right)
}_{A^{(1)}} \\
\to
\underbrace{
\left(\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 1 & 1
\end{array}\right)
}_{L^{(2)}}
A^{(1)}=
\underbrace{
\left(\begin{array}{ccc}
3 & 1 & 0 \\
0 & -1 & -2 \\
0 & 0 & 1
\end{array}\right)}_{A^{(2)}}
A ( 0 ) = 3 6 − 3 1 1 0 0 − 2 3 → L ( 1 ) 1 − 2 1 0 1 0 0 0 1 3 6 − 3 1 1 0 0 − 2 3 = A ( 1 ) 3 0 0 1 − 1 1 0 − 2 3 → L ( 2 ) 1 0 0 0 1 1 0 0 1 A ( 1 ) = A ( 2 ) 3 0 0 1 − 1 0 0 − 2 1
より L = ( L ( 2 ) L ( 1 ) ) − 1 = ( 1 0 0 2 1 0 − 1 − 1 1 ) , U = A ( 2 ) L=(L^{(2)}L^{(1)})^{-1}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{array}\right),U=A^{(2)} L = ( L ( 2 ) L ( 1 ) ) − 1 = 1 2 − 1 0 1 − 1 0 0 1 , U = A ( 2 ) ( L ( 2 ) L ( 1 ) ) − 1 (L^{(2)}L^{(1)})^{-1} ( L ( 2 ) L ( 1 ) ) − 1 A 0 = ( a 1 T , a 2 T , a 3 T ) T A^{0}=({\boldsymbol a_1}^T,{\boldsymbol a_2}^T,{\boldsymbol a_3}^T)^T A 0 = ( a 1 T , a 2 T , a 3 T ) T 2 ‾ a 1 − a 2 , − 1 ‾ a 1 − a 3 , − 1 ‾ a 2 − a 3 \underline{2}{\boldsymbol a_1}-{\boldsymbol a_2}, \underline{-1}{\boldsymbol a_1}-{\boldsymbol a_3}, \underline{-1}{\boldsymbol a_2}-{\boldsymbol a_3} 2 a 1 − a 2 , − 1 a 1 − a 3 , − 1 a 2 − a 3 L L L
この導出過程を見ればなんとなく LU 分解が一意となることは直感的にも納得できるが, 一応証明を与えておく.
LU 分解された L , U L,U L , U
正則行列 A A A A = L U A= L U A = LU A = L 1 U 1 = L 2 U 2 ↔ L 2 − 1 L 1 = U 2 U 1 − 1 A= L_{1}U_{1} = L_{2}U_{2} \leftrightarrow L_{2}^{-1}L_{1} = U_{2}U_{1}^{-1} A = L 1 U 1 = L 2 U 2 ↔ L 2 − 1 L 1 = U 2 U 1 − 1 L L L L 2 − 1 L 1 L_{2}^{-1}L_{1} L 2 − 1 L 1 U U U U 2 U 1 − 1 U_{2}U_{1}^{-1} U 2 U 1 − 1 L L L 1 1 1 L 2 − 1 L 1 = I = U 2 U 1 − 1 L_{2}^{-1}L_{1} = I = U_{2}U_{1}^{-1} L 2 − 1 L 1 = I = U 2 U 1 − 1 L 2 = L 1 , U 2 = U 1 L_{2}=L_{1}, U_{2}=U_{1} L 2 = L 1 , U 2 = U 1
次に, 次のような行列に対する LU 分解を考えてみる.
A = ( 0 1 0 − 8 8 1 2 − 2 0 )
\begin{aligned}
A=\left(\begin{array}{ccc}
0 & 1 & 0 \\
-8 & 8 & 1 \\
2 & -2 & 0
\end{array}\right)
\end{aligned}
A = 0 − 8 2 1 8 − 2 0 1 0
見てわかるように, 前進消去の段階で − 8 ÷ 0 -8\div 0 − 8 ÷ 0 P ∈ R m × n P\in\mathbb{R}^{m\times n} P ∈ R m × n A = P L U A=P L U A = P LU A A A a 1 T {\boldsymbol a_1}^T a 1 T a 2 T {\boldsymbol a_2}^T a 2 T
( 0 1 0 1 0 0 0 0 1 ) ⏟ P ( 1 ) A = ( − 8 8 1 0 1 0 2 − 2 0 ) ⏟ A ′ → ( 1 0 0 0 1 0 1 4 0 1 ) ⏟ L ( 1 ) A ′ = ( − 8 8 1 0 1 0 0 0 1 4 ) ⏟ A ′ ( 1 )
\begin{aligned}
\underbrace{\left(\begin{array}{ccc}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1
\end{array}\right)}_{P^{(1)}}A=
\underbrace{\left(\begin{array}{ccc}
-8 & 8 & 1 \\
0 & 1 & 0 \\
2 & -2 & 0
\end{array}\right)}_{A'}\to
\underbrace{\left(\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
\frac{1}{4} & 0 & 1
\end{array}\right)}_{L^{(1)}}A'=
\underbrace{\left(\begin{array}{ccc}
-8 & 8 & 1 \\
0 & 1 & 0 \\
0 & 0 & \frac{1}{4}
\end{array}\right)}_{A'^{(1)}}
\end{aligned}
P ( 1 ) 0 1 0 1 0 0 0 0 1 A = A ′ − 8 0 2 8 1 − 2 1 0 0 → L ( 1 ) 1 0 4 1 0 1 0 0 0 1 A ′ = A ′ ( 1 ) − 8 0 0 8 1 0 1 0 4 1
より
L = L ( 1 ) − 1 = ( 1 0 0 0 1 0 − 1 4 0 1 ) , U = A ′ ( 1 ) L = {L^{(1)}}^{-1}=\left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ -\frac{1}{4} & 0 & 1\end{array}\right),U={A'}^{(1)} L = L ( 1 ) − 1 = 1 0 − 4 1 0 1 0 0 0 1 , U = A ′ ( 1 )
とおくと A ′ = L ( 1 ) − 1 A ′ ( 1 ) A'={L^{(1)}}^{-1}{A'}^{(1)} A ′ = L ( 1 ) − 1 A ′ ( 1 )
P ( 1 ) A = L U P^{(1)}A=L U P ( 1 ) A = LU
P ( 1 ) P^{(1)} P ( 1 ) P = P ( 1 ) − 1 = P ( 1 ) T {P=P^{(1)}}^{-1}={P^{(1)}}^T P = P ( 1 ) − 1 = P ( 1 ) T
A = P L U ( 0 1 0 − 8 8 1 2 − 2 0 ) = ( 0 1 0 1 0 0 0 0 1 ) ( 1 0 0 0 1 0 − 1 4 0 1 ) ( − 8 8 1 0 1 0 0 0 1 4 )
\begin{aligned}
A&=&P L U\\
\left(\begin{array}{ccc}
0 & 1 & 0 \\
-8 & 8 & 1 \\
2 & -2 & 0
\end{array}\right)&=&
\left(\begin{array}{ccc}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1
\end{array}\right)
\left(\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
-\frac{1}{4} & 0 & 1
\end{array}\right)
\left(\begin{array}{ccc}
-8 & 8 & 1 \\
0 & 1 & 0 \\
0 & 0 & \frac{1}{4}
\end{array}\right)
\end{aligned}
A 0 − 8 2 1 8 − 2 0 1 0 = = P LU 0 1 0 1 0 0 0 0 1 1 0 − 4 1 0 1 0 0 0 1 − 8 0 0 8 1 0 1 0 4 1
ここで a ∈ R n × 1 , v ∈ R n × 1 {\boldsymbol a}\in\mathbb{R}^{n\times 1},{\boldsymbol v}\in\mathbb{R}^{n\times 1} a ∈ R n × 1 , v ∈ R n × 1 P L U a = v P L U {\boldsymbol a}={\boldsymbol v} P LU a = v L ( U a ) ⏟ b = P − 1 v L\underbrace{(U{\boldsymbol a})}_{{\boldsymbol b}}=P^{-1}{\boldsymbol v} L b ( U a ) = P − 1 v b {\boldsymbol b} b P − 1 v P^{-1}{\boldsymbol v} P − 1 v P − 1 = P ( 1 ) P^{-1}=P^{(1)} P − 1 = P ( 1 ) v {\boldsymbol v} v U a = b U{\boldsymbol a}={\boldsymbol b} U a = b a {\boldsymbol a} a
さて, 数学的な言葉では, 以上のように書き下すことで十分であるが, これをプログラムで組むことを考えると, 様々な工夫やアプローチが考えられる. まず L , U L, U L , U 0 0 0 1 1 1 L , U L, U L , U
( − 8 8 1 0 1 0 − 1 4 0 1 4 )
\left(\begin{array}{ccc}
-8 & 8 & 1 \\
0 & 1 & 0 \\
-\frac{1}{4} & 0 & \frac{1}{4}
\end{array}\right)
− 8 0 − 4 1 8 1 0 1 0 4 1
また, 置換行列 P P P クリックで実装を開く.
> : m Data.Array Math.Matrix.LU > lu [[0 ,1 ,0 ],[- 8 ,8 ,1 ],[2 ,- 2 ,0 ]] :: Maybe (PLU Array Int Double )Just (array (0 ,2 ) [(0 ,1 ),(1 ,0 ),(2 ,2 )],- 8.0 8.0 1.0 }- 0.0 1.0 0.0 }- 0.25 0.0 0.25 }> luST' [[0 ,1 ,0 ],[- 8 ,8 ,1 ],[2 ,- 2 ,0 ]] :: Maybe (PLU Array Int Double )Just (array (0 ,2 ) [(0 ,1 ),(1 ,0 ),(2 ,2 )],- 8.0 8.0 1.0 }- 0.0 1.0 0.0 }- 0.25 0.0 0.25 } クリックで実装を開く.
> resolveLinearEq' [[1 ,2 ,7 ,6 ],[2 ,4 ,4 ,2 ],[1 ,8 ,5 ,2 ],[2 ,4 ,3 ,3 ]] [6 ,2 ,12 ,5 ] :: Maybe (Array Int Rational )Just (array (0 ,3 ) [(0 ,(- 3 ) % 1 ),(1 ,2 % 1 ),(2 ,(- 1 ) % 1 ),(3 ,2 % 1 )])> plu = luST' [[1 ,2 ,7 ,6 ],[2 ,4 ,4 ,2 ],[1 ,8 ,5 ,2 ],[2 ,4 ,3 ,3 ]] :: Maybe (PLU Array Int Double )> pluJust (array (0 ,3 ) [(0 ,1 ),(1 ,2 ),(2 ,0 ),(3 ,3 )],2.0 4.0 4.0 2.0 }0.5 6.0 3.0 1.0 }0.5 0.0 5.0 5.0 }1.0 0.0 - 0.2 2.0 }> fromJust plu `assign` (listArray (0 , 3 ) [6 ,2 ,12 ,5 ]) :: Maybe (Array Int Double )Just (array (0 ,3 ) [(0 ,- 3.0 ),(1 ,2.0 ),(2 ,- 1.0 ),(3 ,2.0 )])> fromJust plu `assign` (listArray (0 , 3 ) [1 ,2 ,3 ,4 ]) :: Maybe (Array Int Double )Just (array (0 ,3 ) [(0 ,0.6666666666666667 ),(1 ,0.6666666666666666 ),(2 ,- 1.0 ),(3 ,1.0 )])> fromJust plu `assign` (listArray (0 , 3 ) [5 ,6 ,7 ,8 ]) :: Maybe (Array Int Double )Just (array (0 ,3 ) [(0 ,1.666666666666667 ),(1 ,0.8666666666666667 ),(2 ,- 0.8 ),(3 ,1.2 )])A A − 1 = I A A^{-1} =I A A − 1 = I A − 1 = ( a 1 − 1 , a 2 − 1 , ⋯ , a − 1 m ) , I = ( I 1 , I 2 , ⋯ , I m ) A^{-1}=({\boldsymbol a^{-1}_1},{\boldsymbol a^{-1}_2},\cdots,{\boldsymbol a^{-1}}_m), I=({\boldsymbol I_1},{\boldsymbol I_2},\cdots,{\boldsymbol I_m}) A − 1 = ( a 1 − 1 , a 2 − 1 , ⋯ , a − 1 m ) , I = ( I 1 , I 2 , ⋯ , I m ) A a i − 1 = I i w h e r e i ∈ Z + , 1 ≤ i ≤ m A {\boldsymbol a^{-1}_i}={\boldsymbol I_i}\ {\rm where\ } i\in\mathbb{Z}^{+}, 1\leq i\leq m A a i − 1 = I i where i ∈ Z + , 1 ≤ i ≤ m A = L U A=L U A = LU 1 1 1 m m m a − 1 i {\boldsymbol a^{-1}}_i a − 1 i a − 1 i {\boldsymbol a^{-1}}_i a − 1 i I i I_i I i クリックで実装を開く.
> inverse' [[3 ,1 ,1 ],[5 ,1 ,3 ],[2 ,0 ,1 ]] :: Maybe (Matrix Array Int Rational )Just 1 % 2 (- 1 ) % 2 1 % 1 }1 % 2 1 % 2 (- 2 ) % 1 }- 1 ) % 1 1 % 1 (- 1 ) % 1 }ただし, 逆行列の計算には今述べたようにすべての I − 1 i {\boldsymbol I^{-1}}_i I − 1 i A x = v A{\boldsymbol x}={\boldsymbol v} A x = v A − 1 A^{-1} A − 1
また, LU 分解は行列式の計算も簡単にする. A = L U A= L U A = LU ∣ A ∣ = ∣ L U ∣ = ∣ L ∣ ∣ U ∣ \left|A\right|=\left|L U\right|=\left|L\right|\left|U\right| ∣ A ∣ = ∣ LU ∣ = ∣ L ∣ ∣ U ∣ ∣ L ∣ = 1 \left|L\right|=1 ∣ L ∣ = 1 ∣ A ∣ = ∏ i = 1 n u i i \left|A\right|=\prod_{i=1}^{n}{\boldsymbol u}_{ii} ∣ A ∣ = ∏ i = 1 n u ii P P P S S S ∣ A ∣ = ∣ P ∣ ∣ L ∣ ∣ U ∣ = ( − 1 ) S ∏ i = 1 n u i i \left|A\right|=\left|P\right|\left|L\right|\left|U\right|=(-1)^S \prod_{i=1}^{n}{\boldsymbol u}_{ii} ∣ A ∣ = ∣ P ∣ ∣ L ∣ ∣ U ∣ = ( − 1 ) S ∏ i = 1 n u ii U U U 1 1 1 ∣ A ∣ = ( − 1 ) S ∏ i = 1 n l i i \left|A\right| = (-1)^S \prod_{i=1}^{n}{\boldsymbol l}_{ii} ∣ A ∣ = ( − 1 ) S ∏ i = 1 n l ii クリックで実装を開く.
> determinant $ toMat $ listArray ((0 ,0 ),(2 ,2 )) [3 ,1 ,1 ,5 ,1 ,3 ,2 ,0 ,1 ]2.0 > determinant $ toMat $ listArray ((0 ,0 ),(3 ,3 )) [1 ,2 ,7 ,6 ,2 ,4 ,4 ,2 ,1 ,8 ,5 ,2 ,2 ,4 ,3 ,3 ]120.0 なお, LU 分解は LDU 分解ともいわれることがある. その場合, 上記の U U U D D D A = L D U A = L D U A = L D U
L D U = ( 1 0 ⋯ 0 l 21 1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ l m 1 l m 2 ⋯ 1 ) d i a g ( d 1 , ⋯ , d m ) ( 1 u 12 ⋯ u 1 n 0 1 ⋯ u 2 n ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 )
L D U=
\left(\begin{array}{cccc}
1 & 0 & \cdots & 0 \\
l_{21} & 1 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
l_{m1} & l_{m2} & \cdots & 1
\end{array}\right)
\mathrm{diag}(d_1,\cdots,d_m)
\left(\begin{array}{cccc}
1 & u_{12} & \cdots & u_{1n} \\
0 & 1 & \cdots & u_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 1
\end{array}\right)
L D U = 1 l 21 ⋮ l m 1 0 1 ⋮ l m 2 ⋯ ⋯ ⋱ ⋯ 0 0 ⋮ 1 diag ( d 1 , ⋯ , d m ) 1 0 ⋮ 0 u 12 1 ⋮ 0 ⋯ ⋯ ⋱ ⋯ u 1 n u 2 n ⋮ 1
ところで, A ∈ R m × n A\in\mathbb{R}^{m\times n} A ∈ R m × n A = L D L T A=L D L^{T} A = L D L T
A = A T w h e r e A ∈ R m × n A=A^T\ {\rm where}\ A\in\mathbb{R}^{m\times n} A = A T where A ∈ R m × n A A A A = L D L T A=L D L^{T} A = L D L T
M T = D − 1 U M^T = D^{-1}U M T = D − 1 U A A A A = L D M T A= L D M^T A = L D M T A = A T = ( L D M T ) T = M D L T = L U A = A^T = (L D M^T)^T = M D L^T = L U A = A T = ( L D M T ) T = M D L T = LU M ( D L T ) M(D L^T) M ( D L T ) L U L U LU 定理 1 より LU 分解は一意であるから M = L M=L M = L M = L D M=L D M = L D L D M T = L D L T L D M^T = L D L^T L D M T = L D L T
ここまで述べてきた LU 分解の方法は, 外積形式ガウス法といわれるものであるが, LU 分解の他の方法として内積形式ガウス法(以下 Doolittle 法), クラウト法がある. 簡単のために n = 3 n=3 n = 3 X X X L ∈ R 3 × 3 L\in\mathbb{R}^{3\times 3} L ∈ R 3 × 3 U ∈ R 3 × 3 U\in\mathbb{R}^{3\times 3} U ∈ R 3 × 3 X = L ′ D U ′ X = L' D U' X = L ′ D U ′ L = L ′ D L = L' D L = L ′ D X = L U ′ X = L U' X = L U ′ L , U ′ L,U' L , U ′ U = D U ′ U=D U' U = D U ′ X = L ′ U X = L' U X = L ′ U L ′ U L' U L ′ U
X = L U ′ ↔ ( x 11 x 12 x 13 x 21 x 22 x 23 x 31 x 32 x 33 ) = ( 1 0 0 l 21 1 0 l 31 l 32 1 ) ( u 11 u 12 u 13 0 u 22 u 23 0 0 u 33 ) : = D o o l i t t l e 法 X = L ′ U ↔ ( x 11 x 12 x 13 x 21 x 22 x 23 x 31 x 32 x 33 ) = ( l 11 0 0 l 21 l 22 0 l 31 l 32 l 33 ) ( 1 u 12 u 13 0 1 u 23 0 0 1 ) : = C r o u t 法
\begin{aligned}
X= L U'\leftrightarrow
\left(\begin{array}{ccc}
x_{11} & x_{12} & x_{13} \\
x_{21} & x_{22} & x_{23} \\
x_{31} & x_{32} & x_{33}
\end{array}\right)=\left(\begin{array}{ccc}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{array}\right)\left(\begin{array}{ccc}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{array}\right)&:=&{\rm Doolittle\ 法} \\
X= L' U\leftrightarrow
\left(\begin{array}{ccc}
x_{11} & x_{12} & x_{13} \\
x_{21} & x_{22} & x_{23} \\
x_{31} & x_{32} & x_{33}
\end{array}\right)=\left(\begin{array}{ccc}
l_{11} & 0 & 0 \\
l_{21} & l_{22} & 0 \\
l_{31} & l_{32} & l_{33}
\end{array}\right)\left(\begin{array}{ccc}
1 & u_{12} & u_{13} \\
0 & 1 & u_{23} \\
0 & 0 & 1
\end{array}\right)&:=&{\rm Crout\ 法}
\end{aligned}
X = L U ′ ↔ x 11 x 21 x 31 x 12 x 22 x 32 x 13 x 23 x 33 = 1 l 21 l 31 0 1 l 32 0 0 1 u 11 0 0 u 12 u 22 0 u 13 u 23 u 33 X = L ′ U ↔ x 11 x 21 x 31 x 12 x 22 x 32 x 13 x 23 x 33 = l 11 l 21 l 31 0 l 22 l 32 0 0 l 33 1 0 0 u 12 1 0 u 13 u 23 1 := := Doolittle 法 Crout 法
いま行列 X X X L U ′ L U' L U ′ X = L U ′ X=L U' X = L U ′
( u 11 u 12 u 13 l 21 u 11 l 21 u 12 + u 22 l 21 u 13 + u 23 l 31 u 11 l 31 u 12 + l 32 u 22 l 31 u 13 + l 32 u 23 + u 33 ) = ( 1 0 0 l 21 1 0 l 31 l 32 1 ) ( u 11 u 12 u 13 0 u 22 u 23 0 0 u 33 )
\left(\begin{array}{ccc}
u_{11} & u_{12} & u_{13} \\
l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} \\
l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33}
\end{array}\right)=
\left(\begin{array}{ccc}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{array}\right)\left(\begin{array}{ccc}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{array}\right)
u 11 l 21 u 11 l 31 u 11 u 12 l 21 u 12 + u 22 l 31 u 12 + l 32 u 22 u 13 l 21 u 13 + u 23 l 31 u 13 + l 32 u 23 + u 33 = 1 l 21 l 31 0 1 l 32 0 0 1 u 11 0 0 u 12 u 22 0 u 13 u 23 u 33
よって, 行列 X X X L , U ′ L, U' L , U ′
( x 11 x 12 x 13 x 21 x 22 x 23 x 31 x 32 x 33 ) = ( 1 0 0 x 21 u 11 1 0 x 31 u 11 l 32 1 ) ( u 11 u 12 u 13 0 u 22 u 23 0 0 u 33 ) ∵ x 21 = l 21 u 11 ↔ l 21 = x 21 u 11 , l 31 についても同様 = ( 1 0 0 x 21 x 11 1 0 x 31 x 11 l 32 1 ) ( x 11 x 12 x 13 0 u 22 u 23 0 0 u 33 ) ∵ u 11 = x 11 , u 12 = x 12 , u 13 = x 13 = ( 1 0 0 x 21 x 11 1 0 x 31 x 11 l 32 1 ) ( x 11 x 12 x 13 0 x 22 − x 21 x 11 x 12 x 23 − x 21 x 11 x 13 0 0 u 33 ) ∵ x 22 = l 21 u 12 + u 22 ↔ u 22 = x 22 − l 21 u 12 , u 23 についても同様 = ( 1 0 0 x 21 x 11 1 0 x 31 x 11 x 32 − x 31 x 11 x 12 x 22 − x 21 x 11 x 12 1 ) ( x 11 x 12 x 13 0 x 22 − x 21 x 11 x 12 x 23 − x 21 x 11 x 13 0 0 x 33 − x 31 x 11 x 13 − x 32 − x 31 x 11 x 12 x 22 − x 21 x 11 x 12 ( x 23 − x 21 x 11 x 13 ) ) ∵ x 32 = l 31 u 12 + l 32 u 22 ↔ l 32 = x 32 − l 31 u 12 u 22 , u 33 についても同様
\begin{aligned}
\left(\begin{array}{ccc}
x_{11} & x_{12} & x_{13} \\
x_{21} & x_{22} & x_{23} \\
x_{31} & x_{32} & x_{33}
\end{array}\right)&=&
\left(\begin{array}{ccc}
1 & 0 & 0 \\
\frac{x_{21}}{u_{11}} & 1 & 0 \\
\frac{x_{31}}{u_{11}} & l_{32} & 1
\end{array}\right)
\left(\begin{array}{ccc}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{array}\right)\because
\begin{array}{l}
x_{21}=l_{21}u_{11}\leftrightarrow l_{21}=\frac{x_{21}}{u_{11}}, \\
l_{31}\ {\rm についても同様}
\end{array} \\
&=&
\left(\begin{array}{ccc}
1 & 0 & 0 \\
\frac{x_{21}}{x_{11}} & 1 & 0 \\
\frac{x_{31}}{x_{11}} & l_{32} & 1
\end{array}\right)\left(\begin{array}{ccc}
x_{11} & x_{12} & x_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{array}\right)\because\ u_{11}=x_{11},u_{12}=x_{12},u_{13}=x_{13} \\
&=&
\left(\begin{array}{ccc}
1 & 0 & 0 \\
\frac{x_{21}}{x_{11}} & 1 & 0 \\
\frac{x_{31}}{x_{11}} & l_{32} & 1
\end{array}\right)\left(\begin{array}{ccc}
x_{11} & x_{12} & x_{13} \\
0 & x_{22}-\frac{x_{21}}{x_{11}}x_{12} & x_{23}-\frac{x_{21}}{x_{11}}x_{13} \\
0 & 0 & u_{33}
\end{array}\right)\because
\begin{array}{l}
x_{22}=l_{21}u_{12}+u_{22} \leftrightarrow u_{22}=x_{22}-l_{21}u_{12},\\
u_{23}\ {\rm についても同様}
\end{array} \\
&=&
\left(\begin{array}{ccc}
1 & 0 & 0 \\
\frac{x_{21}}{x_{11}} & 1 & 0 \\
\frac{x_{31}}{x_{11}} & \frac{x_{32}-\frac{x_{31}}{x_{11}}x_{12}}{x_{22}-\frac{x_{21}}{x_{11}}x_{12}} & 1
\end{array}\right)\left(\begin{array}{ccc}
x_{11} & x_{12} & x_{13} \\
0 & x_{22}-\frac{x_{21}}{x_{11}}x_{12} & x_{23}-\frac{x_{21}}{x_{11}}x_{13} \\
0 & 0 & x_{33}-\frac{x_{31}}{x_{11}}x_{13}-\frac{x_{32}-\frac{x_{31}}{x_{11}}x_{12}}{x_{22}-\frac{x_{21}}{x_{11}}x_{12}}(x_{23}-\frac{x_{21}}{x_{11}}x_{13})
\end{array}\right)\\
&\because&
\begin{array}{l}
x_{32}=l_{31}u_{12}+l_{32}u_{22}\leftrightarrow l_{32}=\frac{x_{32}-l_{31}u_{12}}{u_{22}}, \\
u_{33}\ {\rm についても同様}
\end{array}
\end{aligned}
x 11 x 21 x 31 x 12 x 22 x 32 x 13 x 23 x 33 = = = = ∵ 1 u 11 x 21 u 11 x 31 0 1 l 32 0 0 1 u 11 0 0 u 12 u 22 0 u 13 u 23 u 33 ∵ x 21 = l 21 u 11 ↔ l 21 = u 11 x 21 , l 31 についても同様 1 x 11 x 21 x 11 x 31 0 1 l 32 0 0 1 x 11 0 0 x 12 u 22 0 x 13 u 23 u 33 ∵ u 11 = x 11 , u 12 = x 12 , u 13 = x 13 1 x 11 x 21 x 11 x 31 0 1 l 32 0 0 1 x 11 0 0 x 12 x 22 − x 11 x 21 x 12 0 x 13 x 23 − x 11 x 21 x 13 u 33 ∵ x 22 = l 21 u 12 + u 22 ↔ u 22 = x 22 − l 21 u 12 , u 23 についても同様 1 x 11 x 21 x 11 x 31 0 1 x 22 − x 11 x 21 x 12 x 32 − x 11 x 31 x 12 0 0 1 x 11 0 0 x 12 x 22 − x 11 x 21 x 12 0 x 13 x 23 − x 11 x 21 x 13 x 33 − x 11 x 31 x 13 − x 22 − x 11 x 21 x 12 x 32 − x 11 x 31 x 12 ( x 23 − x 11 x 21 x 13 ) x 32 = l 31 u 12 + l 32 u 22 ↔ l 32 = u 22 x 32 − l 31 u 12 , u 33 についても同様
これをみると, 一行目, 一列目, 二行目, 二列目 ⋯ \cdots ⋯ L , U ′ L,U' L , U ′ u i j u_{ij} u ij l i j l_{ij} l ij
D o o l i t t l e 法 : = { { u 1 k = x 1 k u i k = x 1 k − ∑ j = 1 i − 1 l i j u j k , ( i = 2 , 3 , ⋯ , k ) l i k = ( x i k − ∑ j = 1 k − 1 l i j u j k ) u k k , ( i = k + 1 , k + 2 , ⋯ , n )
\begin{aligned}
{\rm Doolittle 法} &:=&
\begin{cases}
\begin{cases}
u_{1k}&=&x_{1k} \\
u_{ik}&=&x_{1k}-\sum^{i-1}_{j=1}l_{ij}u_{jk},\ (i=2,3,\cdots,k)
\end{cases} \\
l_{ik}=\frac{(x_{ik}-\sum^{k-1}_{j=1}l_{ij}u_{jk})}{u_{kk}},\ (i=k+1,k+2,\cdots,n)
\end{cases}
\end{aligned}
Doolittle 法 := ⎩ ⎨ ⎧ { u 1 k u ik = = x 1 k x 1 k − ∑ j = 1 i − 1 l ij u jk , ( i = 2 , 3 , ⋯ , k ) l ik = u kk ( x ik − ∑ j = 1 k − 1 l ij u jk ) , ( i = k + 1 , k + 2 , ⋯ , n )
ただし u k k = 0 u_{kk}=0 u kk = 0 U U U k k k L L L k k k l i k l_{ik} l ik
s i = x i k − ∑ j = 1 k − 1 l i j u j k , ( i = k , ⋯ , n ) s_i = x_{ik}-\sum^{k-1}_{j=1}l_{ij}u_{jk},\ (i=k,\cdots,n) s i = x ik − j = 1 ∑ k − 1 l ij u jk , ( i = k , ⋯ , n )
これにより s i s_i s i k k k u k k = s k u_{kk}=s_{k} u kk = s k U U U k k k L L L l i k = s i u k k l_{ik}=\frac{s_i}{u_{kk}} l ik = u kk s i
Crout 法も同様に, X = L ′ U X=L' U X = L ′ U
( l 11 l 11 u 12 l 11 u 13 l 21 l 21 u 12 + l 22 l 21 u 13 + l 22 u 23 l 31 l 31 u 12 + l 32 l 31 u 13 + l 32 u 23 + l 33 ) = ( l 11 0 0 l 21 l 22 0 l 31 l 32 l 33 ) ( 1 u 12 u 13 0 1 u 23 0 0 1 )
\begin{aligned}
\left(\begin{array}{ccc}
l_{11} & l_{11}u_{12} & l_{11}u_{13} \\
l_{21} & l_{21}u_{12}+l_{22} & l_{21}u_{13}+l_{22}u_{23} \\
l_{31} & l_{31}u_{12}+l_{32} & l_{31}u_{13}+l_{32}u_{23}+l_{33}
\end{array}\right)=
\left(\begin{array}{ccc}
l_{11} & 0 & 0 \\
l_{21} & l_{22} & 0 \\
l_{31} & l_{32} & l_{33}
\end{array}\right)
\left(\begin{array}{ccc}
1 & u_{12} & u_{13} \\
0 & 1 & u_{23} \\
0 & 0 & 1
\end{array}\right)
\end{aligned}
l 11 l 21 l 31 l 11 u 12 l 21 u 12 + l 22 l 31 u 12 + l 32 l 11 u 13 l 21 u 13 + l 22 u 23 l 31 u 13 + l 32 u 23 + l 33 = l 11 l 21 l 31 0 l 22 l 32 0 0 l 33 1 0 0 u 12 1 0 u 13 u 23 1
従ってこの作業の一般形は結果的に
C r o u t 法 : = { l i k = x i k − ∑ j = 1 k − 1 l i j u j k , ( i = k , k + 1 , ⋯ , n ) u k j = ( x k j − ∑ i = 1 k − 1 l k i u i j ) l k k , ( j = k , k + 1 , ⋯ , n )
\begin{aligned}
{\rm Crout 法} &:=&
\begin{cases}
l_{ik}&=&x_{ik}-\sum^{k-1}_{j=1}l_{ij}u_{jk},\ (i=k,k+1,\cdots,n)\\
u_{kj}&=&\frac{(x_{kj}-\sum^{k-1}_{i=1}l_{ki}u_{ij})}{l_{kk}},\ (j=k,k+1,\cdots,n)
\end{cases}
\end{aligned}
Crout 法 := { l ik u kj = = x ik − ∑ j = 1 k − 1 l ij u jk , ( i = k , k + 1 , ⋯ , n ) l kk ( x kj − ∑ i = 1 k − 1 l ki u ij ) , ( j = k , k + 1 , ⋯ , n )
L L L k k k L L L k − 1 k-1 k − 1
参考文献
Richard Hamming (1987) “Numerical Methods for Scientists and Engineers (Dover Books on Mathematics)” Dover Publications, ISBN 9780486652412