前書き

記事の概要

今回は連立方程式を解くときに利用する LU分解 について解説していきたいと思います。
この記事には

数式によるアルゴリズムの解説
LU分解に関係する数学的な話
プログラム実装例

が含まれています。それぞれ自分に必要な部分を参照してもらえればいいなと思っています。

対象読者

大学1年生で勉強する「線形代数学」の基本的な事項は既知として書いていきます。
キーワードは 上・下三角行列 / ガウスの消去法 です。
もし記事中でわからない言葉が出てきた場合には、線形代数学の教科書を読んだり他のサイトを見たりして調べてみてください。

導入

LU分解って何？

$n$ 次正方行列 $A$ を$n$ 次下三角行列 $L$ と $n$ 次上三角行列 $R$ の積に分解することを LU分解 と言います。
すなわち

$$
A = L U
$$

を満たすような下三角行列 $L$ 、上三角行列 $U$ を求めるのがLU分解です。
$A$ が「すべての首座小行列式が0でない」という条件を満たすときにはこのような分解が存在することが知られています。この条件については後で詳しく見ていきますが、「正則である」という条件よりも強い条件を満たす必要があります。

LU分解できると何が嬉しいの？

行列をLU分解できると嬉しいことはたくさんあります。

例えば「連立方程式 $A {\bf x} = {\bf b}$ が解きやすくなる」ということが挙げられます。係数行列 $A$ をLU分解して $ A = L U $ と表現できたと仮定しましょう。このとき、連立方程式は次のように書き直すことができます。

\begin{eqnarray*}
&& A {\bf x} = {\bf b} \\
&\Leftrightarrow& L U {\bf x} = {\bf b} \\
&\Leftrightarrow& L {\bf y} = {\bf b}, U {\bf x} = {\bf y}
\end{eqnarray*}

最終行のような形式に変換できたら、まず第1式を解くことで ${\bf y}$ を求め、次にその ${\bf y}$ を用いて第2式を解いて ${\bf x}$ を得ることができます。(1ステップ目を前進消去、2ステップ目を後退代入と呼びます。)三角行列の逆行列は $O(n^2)$ で簡単に計算できること、一方で一般の連立方程式を解くのに $O(n^3)$ かかることを考えると、LU分解さえできれば連立方程式は簡単に解けることがわかると思います。¹

他にも「行列式や逆行列を簡単に求められる」という利点もあります。LU分解によって得られた $L, U$ はともに三角行列なので

$$
\det(A) = \det(L U) = \det(L) \det(U) = \left( \prod_{i=1}^n l_{ii} \right) \left( \prod_{i=1}^n u_{ii} \right)
$$

として $A$ の行列式を求めることができます。逆行列についても

$$
A^{-1} = (L U)^{-1} = U^{-1} L^{-1}
$$

とすれば簡単に求めることができます。

理論

アルゴリズム

LU分解の方法は様々ありますが、ここでは 外積形式ガウス法 / 内積形式ガウス法 / クラウト法 について紹介していきます。

外積形式ガウス法

外積形式ガウス法 はガウスの消去法の実行過程を行列積として表現したものになります。

ガウスの消去法は、与えられた行列 $A$ に対して行基本変形を繰り返し適用することで階段行列に変換するアルゴリズムです。特に、 $A$ が正則行列である場合には上三角行列にすることができます。ここで、行基本変形というのは以下の3つの操作です。

操作1 : $A$ の第 $i$ 行と第 $j$ 行を入れ替える
操作2 : $A$ の第 $i$ 行を $c$ 倍する
操作3 : $A$ の第 $i$ 行に第 $j$ 行の $c$ 倍を加える

各操作に対応する変換行列は基本行列と呼ばれています。

いま、 $A$ は(行基本変形の操作1を使わない)ガウスの消去法によって上三角行列に変換することができると仮定して、 $A$ のLU分解を求めるアルゴリズムを見ていきます。

まず $A_{11}$ をピボットとして第1列第2行以下に対して掃き出しを実行します。これは行基本変形の操作3を $n-1$ 回適用することで実現できます。各行基本変形に対応する基本行列は下三角行列になることから、変換全体に対応する表現行列は下三角行列になることがわかります。

次に $A_{22}$ をピボットとして第2列第3行以下に対して掃き出しを実行します。このとき、1ステップ目の変換で $A_{21} = 0$ となっているため、第1列は変化しないことに注意してください。この変換も1ステップ目と同様に下三角行列によって表現することができます。

これ以降も同様の手順を踏むことで、各ステップ $k$ の変換行列を下三角行列 $L_k$ とすることができます。最終的に得られた上三角行列を $U$ と書けば

$$
L_n \dots L_2 L_1 A = U
$$

となります。下三角行列の積および(正則な)下三角行列の逆行列が下三角行列になることから、 $L \equiv (L_k \dots L_2 L_1)^{-1}$ も下三角行列になります。したがって

$$
A = L U
$$

という式を得ることができます。これは $A$ のLU分解になっています。このように、ガウスの消去法を実行することでLU分解する方法を外積形式ガウス法と言います。²

内積形式ガウス法

内積形式ガウス法では、 $L$ の対角成分を1に固定した上で、 $A = LU$ の各成分に関する等式から $L, U$の各成分を求めることを考えます。

$A = L U$ の $(i, j)$ 成分について考えると、

$$
a_{ij} = \sum_{k=1}^n l_{ik} u_{kj}
$$

となります。 $L$ は下三角行列、 $U$ は上三角行列になることを仮定しているので

$$
l_{ij} = 0 \text{ if } i < j \\
u_{ij} = 0 \text{ if } i > j
$$

です。この制約から

\begin{eqnarray*}
a_{ij}
&=& \sum_{k: k \le i, k \le j} l_{ik} u_{kj}
= u_{ij} + \sum_{k: k < i, k \le j} l_{ik} u_{kj} \\
a_{ij}
&=& \sum_{k: k \le i, k \le j} l_{ik} u_{kj}
= l_{ij} u_{jj} + \sum_{k: k \le i, k < j} l_{ik} u_{kj} \\
\end{eqnarray*}

となります。ここで $\sum$ は $k$ を(条件を満たす範囲で)動かしたときの総和を表しています。
よって

\begin{eqnarray*}
u_{ij} &=& a_{ij} - \sum_{k: k < i} l_{ik} u_{kj} \\
l_{ij} &=& \frac{1}{u_{jj}} \left( a_{ij} - \sum_{k: k < j} l_{ik} u_{kj} \right) \\
\end{eqnarray*}

という式を導くことができます。左辺の $l, u$ のインデックスに着目すると、第1式では $i \le j$、第2式では $i \ge j$ のみを考慮すれば良いため、右辺の総和の条件が少し変わっていることに注意してください。
特に

\begin{eqnarray*}
u_{1j} &=& a_{1j} \\
l_{i1} &=& \frac{a_{i1}}{u_{11}} \\
\end{eqnarray*}

が成立します。このように、内積形式ガウス法では $L$ の第 $1 \dots k$ 列、 $U$ の第 $1 \dots k$ 行から $L$ の第 $k+1$ 列、 $U$ の第 $k+1$ 行を順に求めていきます。

クラウト法

内積形式ガウス法では、 $U$ の対角成分を1に固定した上で、 $A = LU$ の各成分に関する等式から $L, U$の各成分を求めることを考えます。

$L, U$ の各成分を求める式は内積形式ガウス法のときと同様に導出することができて

\begin{eqnarray*}
u_{ij} &=& \frac{1}{l_{ii}} \left( a_{ij} - \sum_{k: k < i} l_{ik} u_{kj} \right) \\
l_{ij} &=& a_{ij} - \sum_{k: k < j} l_{ik} u_{kj} \\
\end{eqnarray*}

となります。

数学的な話

LU分解可能であるための条件(命題)

与えられた行列 $A$ がLU分解可能であるための必要十分条件は「 $A$ のすべての首座小行列式が0でない」という条件です。

行列 $A$ の小行列式は、 $A$ から1つ以上の行/列を削除した結果得られる正方行列の行列式です。小行列式のうち、第1行から第 $k$ 行、第1列から第 $k$ 列だけを残した行列の行列式が ( $k$ 次)首座小行列式 です。 $n$ 次正方行列 $A$ に対して首座小行列式は $n$ 個存在します。

$A$ が $n$ 次正方行列ならば $A$ の $n$ 次首座小行列式は $A$ 自身になります。そのため、 $A$ がLU分解可能であるためには少なくとも $A$ は正則である必要があります。この意味で、「すべての首座小行列式が0でない」という条件は「正則である」という条件よりも真に強い制約となっています。³

LU分解可能であるための条件(証明)

$n$ 次正方行列 $A$ がLU分解可能であるための必要十分条件は、 $A$ のすべての首座小行列式が0でないことである。

この節では上の命題の証明をしていますが、長くなってしまうので折りたたんでいます。必要があればクリックして展開してください。

$\Rightarrow )$ LU分解可能ならばすべての首座小行列式が0でない

行列のサイズ $n$ に関する帰納法によって証明します。

まず、$n = 1$ の場合を考えます。 $A$ の首座小行列式は $\det(A)$ に限られますが、その値は0になることはありません。なぜならば、外積形式ガウス法のアルゴリズムから $A$ は正則である必要があるからです。

次に、 $n \ge 1$ の場合を考えます。 $n$ 次正方行列に対して条件 $\Rightarrow )$ が成立すると仮定します。
$n+1$ 次正方行列 $A$ が $A = L U$ とLU分解できたと仮定します。この $A$ をブロック分割すると

\begin{eqnarray*}
\left(
\begin{array}{c|c}
A_{1:n, 1:n} & A_{1:n, n+1} \\
\hline
A_{n+1, 1:n} & A_{n+1, n+1}
\end{array}
\right)
=
\left(
\begin{array}{c|c}
L_{1:n, 1:n} & L_{1:n, n+1} \\
\hline
L_{n+1, 1:n} & L_{n+1, n+1}
\end{array}
\right)

\left(
\begin{array}{c|c}
U_{1:n, 1:n} & U_{1:n, n+1} \\
\hline
U_{n+1, 1:n} & U_{n+1, n+1}
\end{array}
\right)
\\
\end{eqnarray*}

と書くことができて

$$
A_{1:n, 1:n} = L_{1:n, 1:n} U_{1:n, 1:n} + L_{1:n, n+1} U_{n+1, 1:n}
$$

となります。$L$ は下三角、 $U$ は上三角なので第2項は0となります。したがって

$$
A_{1:n, 1:n} = L_{1:n, 1:n} U_{1:n, 1:n}
$$

という式が得られます。この式から、 $n$ 次正方行列 $A_{1:n, 1:n}$ がLU分解可能であることがわかります。
すると仮定より $A_{1:n, 1:n}$ のすべての首座小行列式が非零であることが従います。 $A$ 自身の行列式が非零であることと合わせると、 $n+1$ 正方行列でも条件 $\Rightarrow )$ が成立します。

以上より、すべての $n$ に関して条件 $\Rightarrow )$ が成立しています。

$\Leftarrow )$ すべての首座小行列式が0でないならばLU分解可能である

外積形式ガウス法のアルゴリズムを実際に回すことで証明します。

$n$ 次正方行列 $A$ のすべての首座小行列式が0でないと仮定します。
まず $k = 1$ 列目の掃き出しを考えます。
$A_{11}$ は0でないため、 $(1, 1)$ 成分をピボットとして1列目を掃き出すことができます。

次に $k \ge 1$ 列目まで掃き出しが完了したと仮定して、 $k+1$ 列目の掃き出しを行うことを考えます。この時点で $k+1$ 次首座行列 $A_{1:k+1, 1:k+1}$ は対角成分が非零の上三角行列になっていることに注意します。
$k$ 列目までの掃き出し操作によって $A_{1:k+1, 1:k+1}$ の行列式が変化していないことを考えると、 $A_{k+1,k+1}$ は非零になることがわかります。⁴ よって $k+1$ 列目についてもそれまでと同様に対角成分 $(k+1, k+1)$ をピボットとして掃き出しを実行することができます。

以上より、外積形式ガウス法を正しく実行することができるため、条件 $\Leftarrow )$ が成立しています。

実装

LU分解(外積形式ガウス法・内積形式ガウス法・クラウト法)のPythonによる実装はGitHubにあります。

内積形式ガウス法・クラウト法の実装では通常、メモリを削減するために $L, U$ を1つの行列でまとめて管理することが多いですが、ここでは $L, U$ をそれぞれ別の行列で管理しています。また、外積形式ガウス法の実装では本記事の説明に寄せるために(基本行列による)行列演算を使っています。実行効率を犠牲にしてわかりやすさを優先しているため、これよりも実行効率の良い実装があるということを心に留めて見てもらえればと思います。

内積形式ガウス法・クラウト法の実装を見てもらうと、 for ループが3重になっていることからLU分解の時間計算量は $O(n^3)$ となることがわかると思います。

後書き

Qiita初投稿記事でしたが、頑張って書き上げることができました。
これからも「チョットワカル」記事を増やしていけたらいいなと思います。

参考文献・リンク

https://www.cspp.cc.u-tokyo.ac.jp/hanawa/class/spc2016s/sp20160614-2.pdf

LU分解それ自体に $O(n^3)$ の時間を要するため、1つの連立方程式を解くという場合にはそれほど嬉しくはありません。ただ、係数行列 $A$ が共通で $b$ が異なるような複数の連立方程式があったときにはLU分解の恩恵を受けることができます。 ↩
ガウスの消去法において行基本変形の操作1を許すときには PLU分解 となります。 ↩
内積形式ガウス法・クラウト法のアルゴリズムをよく見ると、 $L, U$ ともに対角成分が0でない三角行列になっており、ゆえに正則であることが従います。LU分解が正しく動作すると仮定すると、 $A$ も正則である必要があるということがわかると思います。 ↩
対角成分をスカラー倍する操作も含めて考えると、厳密には行列式は不変ではありません。しかしこの場合でも「行列式が非零かどうか」は変わることはありません。 ↩