LU 分解とは｜求め方と連立方程式への使い方

正方行列 $A$ を、下三角行列 $L$ と上三角行列 $U$ の積に書き直すのが LU 分解です。

$$A=LU$$

新しい手法というより、掃き出し法で消去した記録をそのまま行列にしたものです。消去のたびに使った倍率が $L$ に、変形後の行列が $U$ に残ります。

いちど分解しておけば、同じ係数行列に対する連立方程式が代入だけで解けます。行列式も対角成分を掛けるだけで求まります。

<svg viewBox="0 0 460 210" width="100%" style="max-width:460px;display:block;margin:0 auto">
	<rect x="30" y="45" width="100" height="100" fill="none" stroke="#c8c8c8" stroke-width="1"></rect>
	<text x="80" y="100" font-size="15" fill="#1f1f1f" text-anchor="middle">A</text>
	<text x="80" y="34" font-size="11" fill="#9a9a9a" text-anchor="middle">もとの行列</text>
	<text x="150" y="102" font-size="15" fill="#9a9a9a" text-anchor="middle">=</text>
	<polygon points="175,45 175,145 275,145" fill="#1b81e0" fill-opacity="0.12"></polygon>
	<rect x="175" y="45" width="100" height="100" fill="none" stroke="#c8c8c8" stroke-width="1"></rect>
	<text x="205" y="128" font-size="15" fill="#1b81e0" text-anchor="middle">L</text>
	<text x="225" y="34" font-size="11" fill="#9a9a9a" text-anchor="middle">下三角（対角は 1）</text>
	<polygon points="320,45 420,45 420,145" fill="#1b81e0" fill-opacity="0.12"></polygon>
	<rect x="320" y="45" width="100" height="100" fill="none" stroke="#c8c8c8" stroke-width="1"></rect>
	<text x="390" y="70" font-size="15" fill="#1b81e0" text-anchor="middle">U</text>
	<text x="370" y="34" font-size="11" fill="#9a9a9a" text-anchor="middle">上三角</text>
	<line x1="175" y1="45" x2="275" y2="145" stroke="#9a9a9a" stroke-width="1" stroke-dasharray="3 3"></line>
	<line x1="320" y1="45" x2="420" y2="145" stroke="#9a9a9a" stroke-width="1" stroke-dasharray="3 3"></line>
	<text x="230" y="180" font-size="11" fill="#9a9a9a" text-anchor="middle">対角より上と下に、それぞれ成分が集まった 2 つの行列に分けます</text>
	<text x="230" y="198" font-size="11" fill="#9a9a9a" text-anchor="middle">白い側はすべて 0 です</text>
</svg>

LU 分解とは

$n$ 次正方行列 $A$ に対して、対角成分がすべて $1$ の下三角行列 $L$ と、上三角行列 $U$ を見つけて $A=LU$ と書くことを LU 分解といいます。

$L$ の対角を $1$ に決めておくのは、分解を一通りに定めるためです。この決め方をドゥーリトル法といいます。逆に $U$ の対角を $1$ にする流儀もあり、そちらはクラウト法と呼ばれます。

対角を別に取り出した $A=LDU$ という書き方もあります。$L$ と $U$ の対角をともに $1$ にし、$D$ に大きさを集める形です。

掃き出しの前半がそのまま分解になる

上三角に変形する操作は、ある行の定数倍を下の行から引くことの繰り返しです。この倍率を覚えておくだけで $L$ が作れます。

第 $j$ 列で第 $i$ 行を消すときの倍率を ${\ell }_{ij}$ と書くと、$L$ の $(i,j)$ 成分がその ${\ell }_{ij}$ になります。変形し終えた行列が $U$ です。

消去は「引く」操作、復元は「足す」操作です。$L$ は引いた分を足し戻す行列で、だから対角より下に倍率がそのまま並びます。

例：3 次の行列を LU 分解する

次の行列を分解します。

$$A=\left(\begin{array}{ccc}2& 1& 1\\ 4& -6& 0\\ -2& 7& 2\end{array}\right)$$

第 1 列を消します。第 2 行からは第 1 行の $2$ 倍を引き、第 3 行からは $-1$ 倍を引きます。倍率は ${\ell }_{21}=2$、${\ell }_{31}=-1$ です。

第 2 行は $(0,-8,-2)$、第 3 行は $(0,8,3)$ になります。続けて第 2 列を消すと、倍率は ${\ell }_{32}=-1$ で、第 3 行が $(0,0,1)$ になります。

倍率を並べたものが $L$、変形後が $U$ です。

$$L=\left(\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ -1& -1& 1\end{array}\right),U=\left(\begin{array}{ccc}2& 1& 1\\ 0& -8& -2\\ 0& 0& 1\end{array}\right)$$

掛け戻して確かめる

$LU$ を計算すると、$L$ の第 $i$ 行が $U$ の各行をどれだけ混ぜるかを表していると分かります。

第 2 行なら $2\times (2,1,1)+1\times (0,-8,-2)=(4,-6,0)$ です。もとの第 2 行に戻ります。

第 3 行は $-1\times (2,1,1)-1\times (0,-8,-2)+1\times (0,0,1)=(-2,7,2)$ です。こちらも一致します。

消去で引いたものを足し戻しているだけなので、当然そうなります。

連立方程式は 2 段階に分かれる

$Ax=b$ に $A=LU$ を入れると $L(Ux)=b$ です。$Ux=y$ と置けば、三角の方程式が 2 つ並びます。

まず $Ly=b$ を上から順に解きます。$L$ の対角が $1$ なので、割り算すら要りません。これを前進代入といいます。

次に $Ux=y$ を下から順に解きます。こちらは後退代入です。どちらも $n^2$ 程度の手数で終わります。

<svg viewBox="0 0 460 190" width="100%" style="max-width:460px;display:block;margin:0 auto">
	<rect x="20" y="60" width="110" height="50" fill="none" stroke="#c8c8c8" stroke-width="1"></rect>
	<text x="75" y="90" font-size="13" fill="#1f1f1f" text-anchor="middle">A x = b</text>
	<line x1="140" y1="85" x2="185" y2="85" stroke="#9a9a9a" stroke-width="1"></line>
	<polygon points="185,85 177,81 177,89" fill="#9a9a9a"></polygon>
	<text x="162" y="76" font-size="11" fill="#9a9a9a" text-anchor="middle">分解</text>
	<rect x="195" y="25" width="110" height="50" fill="#1b81e0" fill-opacity="0.12" stroke="#1b81e0" stroke-width="1"></rect>
	<text x="250" y="55" font-size="13" fill="#1f1f1f" text-anchor="middle">L y = b</text>
	<rect x="195" y="105" width="110" height="50" fill="#1b81e0" fill-opacity="0.12" stroke="#1b81e0" stroke-width="1"></rect>
	<text x="250" y="135" font-size="13" fill="#1f1f1f" text-anchor="middle">U x = y</text>
	<line x1="250" y1="75" x2="250" y2="100" stroke="#9a9a9a" stroke-width="1"></line>
	<polygon points="250,105 246,97 254,97" fill="#9a9a9a"></polygon>
	<text x="320" y="55" font-size="11" fill="#9a9a9a">前進代入（上から順に）</text>
	<text x="320" y="135" font-size="11" fill="#9a9a9a">後退代入（下から順に）</text>
	<text x="230" y="180" font-size="11" fill="#9a9a9a" text-anchor="middle">三角の方程式は代入だけで解けます</text>
</svg>

例：LU 分解で連立方程式を解く

さきほどの $A$ に対して、右辺を $b=(5,-2,9)$ とします。

$Ly=b$ を上から解きます。第 1 式から $y_1=5$、第 2 式の $2y_1+y_2=-2$ から $y_2=-12$ です。

第 3 式は $-y_1-y_2+y_3=9$ なので、$y_3=9+5-12=2$ になります。割り算が 1 度も出てきません。

次に $Ux=y$ を下から解きます。$x_3=2$、$-8x_2-2\cdot 2=-12$ から $x_2=1$、$2x_1+1+2=5$ から $x_1=1$ です。

解は $x=(1,1,2)$ になります。$Ax$ を計算すると $(5,-2,9)$ に戻ります。

右辺が変わっても分解を使い回せる

同じ $A$ で右辺だけが違う方程式を何本も解く場面はよくあります。掃き出し法をやり直すと毎回 $2n^3/3$ 程度かかります。

LU 分解なら、分解は 1 回で済みます。右辺が増えても代入の $n^2$ が追加されるだけです。

逆行列を作って掛ける方法より速く、しかも誤差の面でも有利です。逆行列を陽に求める必要は、実はほとんどありません。

行列式は U の対角の積になる

三角行列の行列式は対角成分の積です。$L$ の対角はすべて $1$ なので $\det L=1$ になります。

$$\det A=\det L\cdot \det U=u_{11}{u}_{22}\cdots u_{nn}$$

余因子展開が $n!$ 個の項を持つのに対して、こちらは分解の $2n^3/3$ 回で済みます。大きな行列式を数値で求めるときの標準的なやり方です。

例：行列式を LU 分解から求める

さきほどの $U$ の対角は $2$、$-8$、$1$ でした。

$$\det A=2\times (-8)\times 1=-16$$

余因子展開で確かめても同じ値になります。分解ができていれば、行列式は掛け算 2 回で出ます。

分解が存在する条件

行の入れ替えなしで LU 分解ができるのは、左上から取った小行列式がすべて $0$ でないときです。これを首座小行列式といいます。

$$\det A_k\ne 0(k=1,2,\dots ,n-1)$$

このとき $L$ の対角を $1$ に決めれば、分解は一通りに決まります。$LDU$ の形なら、正則行列に対して一意性がそのまま成り立ちます。

さきほどの例では首座小行列式が $2$、$-16$、$-16$ で、どれも $0$ ではありません。だから入れ替えなしで分解できました。

例：ピボットが 0 で分解できない行列

次の行列は、行を入れ替えないと分解できません。

$$A=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$$

$A=LU$ と書けたとすると、$(1,1)$ 成分の比較から $u_{11}=0$ が出ます。

すると第 2 行の $(1,1)$ 成分は ${\ell }_{21}{u}_{11}=0$ になり、$1$ にはなりません。矛盾するので分解は存在しません。

正則であっても分解できない例がある、というのがこの行列の意味です。首座小行列式の $1$ 次が $0$ になっているのが原因です。

行の入れ替えを入れた PA = LU

行を先に入れ替えておけば、どんな正則行列でも分解できます。入れ替えを表す置換行列を $P$ と書きます。

$$PA=LU$$

$P$ は単位行列の行を並べ替えたもので、$\det P$ は入れ替えの回数が偶数なら $1$、奇数なら $-1$ です。行列式は符号のぶんだけ変わります。

$$\det A=\det P^{-1}\cdot u_{11}\cdots u_{nn}$$

方程式を解くときは、右辺にも同じ入れ替えを施して $Ly=Pb$ から始めます。

例：置換を入れて分解する

左上が $0$ の行列で試します。

$$A=\left(\begin{array}{ccc}0& 2& 1\\ 2& 1& 0\\ 1& 1& 1\end{array}\right)$$

第 1 行と第 2 行を入れ替えると、左上が $2$ になって消去を始められます。倍率は ${\ell }_{31}=1/2$、続いて ${\ell }_{32}=1/4$ です。

$$L=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ \frac{1}{2}& \frac{1}{4}& 1\end{array}\right),U=\left(\begin{array}{ccc}2& 1& 0\\ 0& 2& 1\\ 0& 0& \frac{3}{4}\end{array}\right)$$

$U$ の対角の積は $2\times 2\times \frac{3}{4}=3$ です。入れ替えが 1 回なので符号が変わり、$\det A=-3$ になります。

部分ピボット選択と数値の安定性

入れ替えは、ピボットが $0$ のときだけの応急処置ではありません。$0$ に近い値で割ると誤差が拡大するので、実際の計算では毎回いちばん絶対値の大きい成分を選んで上へ持ってきます。

この選び方を部分ピボット選択といいます。倍率がすべて $1$ 以下に収まるので、消去のあいだに値が暴れにくくなります。

入れ替えをしない素朴な消去は、理論上は正しくても数値計算では不安定です。数値解析の教科書がピボット選択を必ず扱うのは、このためです。

LDU 分解と対称行列

対角を分けた $A=LDU$ で $A$ が対称なら、$U$ は $L$ の転置になります。$A=LD{L}^{\top }$ と書ける形です。

さらに $A$ が正定値なら $D$ の成分がすべて正なので、平方根を $L$ に配って $A=L^{\prime }{L}^{\prime \top }$ にできます。これがコレスキー分解です。

コレスキー分解ではピボット選択が要らず、計算量も LU のおよそ半分の $n^3/3$ で済みます。対称正定値と分かっているなら、こちらを使います。

計算量

分解そのものは $2n^3/3$ 回程度の掛け算と足し算で終わります。前進代入と後退代入は、それぞれ $n^2$ 程度です。

つまり最初の 1 本だけが重く、2 本目からは軽くなります。右辺が $m$ 本あるときの総手数は $2n^3/3+m{n}^2$ の形になります。

$L$ と $U$ は $A$ と同じ配列に上書きできます。$L$ の対角は $1$ と決まっているので、記録する必要がないからです。

ブロックで見るとシューア補元

行を 2 つの組に分けて、ブロックのまま同じ消去をすると、右下にシューア補元が現れます。

$$\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)=\left(\begin{array}{cc}I& O\\ A_{21}{A}_{11}^{-1}& I\end{array}\right)\left(\begin{array}{cc}{A}_{11}& A_{12}\\ O& S\end{array}\right)$$

ここで $S=A_{22}-A_{21}{A}_{11}^{-1}{A}_{12}$ です。左の行列が $L$ にあたり、右の行列を同じ手続きで分解し続ければ LU 分解が組み上がります。

成分ごとの消去は、$1$ 次のブロックに区切った場合にあたります。詳しくはブロック行列とシューア補元の記事にあります。

他の話題とのつながり

消去の手順そのものは掃き出し法の記事で、変形を行列の積に書き直す話は基本行列の記事で扱っています。LU 分解は、その 2 つを「分解」という形にまとめたものです。

行列式を求める道具としては、クラメルの公式と対照的です。あちらは理論的な表示を与え、こちらは数値を求める手順を与えます。

直交行列を使う QR 分解は、計算量が LU の 2 倍ほどかかる代わりに数値的な安定性で勝ります。最小二乗法ではそちらが選ばれます。

正方でない行列や階数が落ちた行列では、分解の形が変わります。解の存在と自由度の話は、連立一次方程式の解の構造の記事にあります。

よくある誤り

参考文献