はじめに
行列分解は、もとの行列Aを積の形式に分解する操作です。行列分解の種類に代表的なものをピックアップしてあります。今回はこの中からLU分解について見ていきますが、LU分解の操作はライブラリに任せることにしてアルゴリズムの詳細には触れません。
LU分解の用途
LU分解の主な用途は次の通りです。今回はこの中から「1. 連立1次方程式の求解」と「2. 行列式の計算」について見ていきます。
- 連立1次方程式の求解
- 行列式の計算
- 逆行列の計算
- 行列のランク判定
- 条件数の推定（数値計算の信頼性を評価する指標）
- 前処理付き共役勾配法
- ガウス分布の尤度計算
LU分解の定義
LU分解は、正方行列Aを下三角行列Lと上三角行列Uの積に分解します。
成分表示は次の通りです。
\tag{1.2}
\begin{pmatrix}
a_{11}&a_{12}&\cdots&a_{1n}\\
a_{21}&a_{22}&\cdots&a_{2n}\\
\vdots&\vdots&\ddots&\vdots\\
a_{n1}&a_{n2}&\cdots&a_{nn}
\end{pmatrix}=
\begin{pmatrix}
1&0&\cdots&0\\
l_{21}&1&\cdots&0\\
\vdots&\vdots&\ddots&\vdots\\
l_{n1}&l_{n2}&\cdots&1
\end{pmatrix}
\begin{pmatrix}
u_{11}&u_{12}&\cdots&u_{1n}\\
0&u_{22}&\cdots&u_{2n}\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\cdots&a_{nn}
\end{pmatrix}
LU分解のアルゴリズムには「Doolittle法」や「Crout法」などがありますが、冒頭に書いたように詳細には触れません。
また、現実的な問題になりますが、実際の数値計算では数値的な安定性を図るためにLU分解ではなくPLU分解（部分ピボット選択付きLU分解）が使われます。実際のところ、今回の実装に使っているLinearAlgebra.lu関数はデフォルトでPLU分解を行います。しかし、PLU分解だと狙い通りの出力形式にならないのでLU分解（LinearAlgebra.lu関数の第2引数に「Val(false)」を渡す）で話を進めます。
連立１次方程式の求解
計算方法
1. 方程式の行列化
次のような連立3元1次方程式を解きたいとします。
\tag{2.1}
\begin{cases}
a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}=b_{1}\\
a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}=b_{2}\\
a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}=b_{3}\\
\end{cases}
これを行列とベクトルを用いて表現します。
\tag{2.2}
\begin{pmatrix}
a_{11}&a_{12}&a_{13}\\
a_{21}&a_{22}&a_{23}\\
a_{31}&a_{32}&a_{33}
\end{pmatrix}
\begin{pmatrix}
x_{1}\\x_{2}\\x_{3}
\end{pmatrix}=
\begin{pmatrix}
b_{1}\\b_{2}\\b_{3}
\end{pmatrix}
これらの 3つの項について、左から順番に A,\bold{x},\bold{b} と置きます。
\tag{2.3}A\bold{x}=\bold{b}
-
A\in\R^{n×n}は係数行列
-
\bold{x}\in\R^{n}は解ベクトル
-
\bold{b}\in\R^{n}は定数項ベクトル
2. LU分解
係数行列AをLU分解します。式(1.1)を再掲しましょう。
これを式(2.3)の方程式に代入します。
\tag{2.5}LU\bold{x}=\bold{b}
3. 前進代入
ここで中間変数ベクトル\bold{y}を導入して U\bold{x}=\bold{y} と置きます。
\tag{2.6}L\bold{y}=\bold{b}
まずは、この方程式のベクトル\bold{y}を求めます。行列Lは下三角行列なので、ベクトル\bold{y}を求めるには順方向（ y_1からy_nの順）に解いていきます。
4. 後退代入
中間変数ベクトル\bold{y}を求めたら、これを定数項ベクトルに見立てて、最終的な解であるベクトル\bold{x}を求めます。
\tag{2.7}U\bold{x}=\bold{y}
行列Uは上三角行列なので、ベクトル\bold{x}を求めるには逆方向（ x_nからx_1の順）に解いていきます。
実例
1. 方程式の行列化
次の連立3元1次方程式を解きます。
\tag{3.1}
\begin{array}{l}
\begin{cases}
2x_{1}&+\,\,\,x_{2}&+\,\,\,x_{3}&=4\,\,\,\\
4x_{1}&+3x_{2}&+3x_{3}&=10\\
8x_{1}&+7x_{2}&+9x_{3}&=36\\
\end{cases}
\end{array}
A\bold{x}=\bold{b} の形にします。
\tag{3.2}
\def\arraystretch{1.2}
\begin{array}{cccc}
&\underline{\begin{pmatrix}2&1&1\\4&3&3\\8&7&9\end{pmatrix}}
&\underline{\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}}&=
&\underline{\begin{pmatrix}4\\10\\36\end{pmatrix}}\\
&A&\bold{x}&&\bold{b}
\end{array}
2. LU分解
係数行列AをLU分解します。
Julia
using LinearAlgebra
A = [ 2 1 1;
4 3 3;
8 7 9 ]
b = [4, 10, 36]
F = lu(A, Val(false))
実行結果
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
1.0 0.0 0.0
2.0 1.0 0.0
4.0 3.0 1.0
U factor:
3×3 Matrix{Float64}:
2.0 1.0 1.0
0.0 1.0 1.0
0.0 0.0 2.0
下三角行列Lと上三角行列Uが求められたので書き写しておきます。
\tag{3.3}
L=\begin{pmatrix}1&0&0\\2&1&0\\4&3&1\end{pmatrix},\enspace
U=\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix}
3. 前進代入
L\bold{y}=\bold{b} の形にします。
\tag{3.4}
\def\arraystretch{1.2}
\begin{array}{cccc}
&\underline{\begin{pmatrix}1&0&0\\2&1&0\\4&3&1\end{pmatrix}}
&\underline{\begin{pmatrix}y_{1}\\y_{2}\\y_{3}\end{pmatrix}}&=
&\underline{\begin{pmatrix}4\\10\\36\end{pmatrix}}\\
&{L}&{\bold{y}}&&{\bold{b}}
\end{array}
ベクトル\bold{y}をy_1（先頭）から解いていきます。
\tag{3.5}
\begin{array}{l}
y_{1}:\enspace y_{1}=4\\
y_{2}:\enspace 2y_{1}+y_{2}=10\enspace\Rightarrow\enspace y_{2}=10-2\cdot4=2\\
y_{3}:\enspace 4y_{1}+3y_{2}+y_{3}=36\enspace\Rightarrow\enspace y_{3}=36-4\cdot4-3\cdot2=14
\end{array}
式(3.5)をベクトル\bold{y}にまとめます。
\tag{3.6}\bold{y}=\begin{pmatrix}4\\2\\14\end{pmatrix}
前進代入の実装と実行結果は次のようになります。
実行結果
3-element Vector{Float64}:
4.0
2.0
14.0
4. 後退代入
U\bold{x}=\bold{y} の形にします。
\tag{3.7}
\def\arraystretch{1.2}
\begin{array}{cccc}
&\underline{\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix}}
&\underline{\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}}&=&
\underline{\begin{pmatrix}4\\2\\14\end{pmatrix}}\\
&{U}&{\bold{x}}&&{\bold{y}}
\end{array}
ベクトル\bold{x}をx_3（最後）から解いていきます。
\tag{3.8}
\begin{array}{l}
x_{3}:\enspace 2x_{3}=14\enspace\Rightarrow\enspace x_{3}=7\\
x_{2}:\enspace x_{2}+x_{3}=2\enspace\Rightarrow\enspace x_{2}=2-7=-5\\
x_{1}:\enspace 2x_{1}+x_{2}+x_{3}=4\enspace\Rightarrow\enspace 2x_{1}-5+7=4\enspace\Rightarrow\enspace 2x_{1}=2\enspace\Rightarrow\enspace x_{1}=1
\end{array}
連立方程式の最終的な解は、式(3.8)を解ベクトル\bold{x}にまとめて完成です。
\tag{3.9}\bold{x}=\begin{pmatrix}1\\-5\\7\end{pmatrix}
後退代入の実装と実行結果は次のようになります。
実行結果
3-element Vector{Float64}:
1.0
-5.0
7.0
5. 検算
検算は式(3.9)を式(3.2)に代入します。
\tag{3.10}
\begin{array}{l}
\begin{cases}
2\cdot1+1\cdot(-5)+1\cdot7=2-5+7=4 ✓\\
4\cdot1+3\cdot(-5)+3\cdot7=4-15+21=10 ✓\\
8\cdot1+7\cdot(-5)+9\cdot7=8-35+63=36 ✓\\
\end{cases}
\end{array}
6. バックスラッシュ（\）演算子
LU分解に関心がなく、単に連立方程式を解きたいときはバックスラッシュ（\）演算子を使うのが簡単です。この演算子を行列に適用すると A\bold{x}=\bold{b} を満たすベクトル\bold{x}を求めることができます。この演算では浮動小数点演算による丸め誤差が発生するので、これを避けたいときは演算に先立って行列を有理数（Rational型）に変換しておきます。
Julia
using LinearAlgebra
A = [ 2 1 1;
4 3 3;
8 7 9 ]
b = [4, 10, 36]
x = Rational.(A) \ Rational.(b)
実行結果
3-element Vector{Rational{Int64}}:
1
-5
7
結果は式(3.9)の解ベクトル\bold{x}と同じです。
行列式の計算
計算方法
LU分解によって求められた行列について、その行列式には次のような関係があります。
\tag{4.1}\det(A)=\det(L)\det(U)
下三角行列Lの対角成分はすべて1なのでこれらを掛けても1です。結局のところ、行列式を計算するには上三角行列Uの対角成分を掛けるだけで良いことになります。
\tag{4.2}\det(A)=\prod_{i=1}^{n}U_{ii}
実例
連立方程式の式を流用します。式(3.2)から係数行列Aを、式(3.3)から下三角行列Lと上三角行列Uを再掲します。
\tag{5.1}
A=\begin{pmatrix}2&1&1\\4&3&3\\8&7&9\end{pmatrix}, L=\begin{pmatrix}1&0&0\\2&1&0\\4&3&1\end{pmatrix}, U=\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix}
式(4.2)を用いて行列式を計算します。これは、上三角行列Uの対角成分を掛けるだけなので次のようになります。
\tag{5.2}\det(A)=2\cdot1\cdot2=4
サラスの方法を用いても同じ結果が得られることを確認します。
\tag{5.3}
\def\arraystretch{1.5}
\begin{array}{ccc:ccc}
2&1&1&2&1&1\\
4&3&3&4&3&3\\
8&7&9&8&7&9\\[2pt]
\hline \\[-18pt]
-42&-36&-24&+54&+24&+28&=4
\enspace ✓
\end{array}
Discussion