LU分解

数値解析と線型代数において、下三角行列と上三角行列の積として行列を下上( LU )分解または因数分解で因数分解します (行列分解を参照)。積には、置換行列も含まれる場合があります。LU 分解は、ガウス消去法の行列形式と見なすことができます。コンピューターは通常、LU 分解を使用して 2 次方程式の正方行列を解きます。また、LU 分解は、行列の逆行列を求めるときや行列の行列式を計算するときにも重要な手順です。 LU分解は、1938年にポーランドの天文学者タデウシュ・バナキエヴィチによって導入されました^。[1]引用すると、「ガウスとドゥーリットルは、この[消去法]を対称方程式にのみ適用したようです。より最近の著者、たとえばエイトキン、バナキエヴィチ、ドワイヤー、クラウトは、この方法またはそのバリエーションを非対称問題と関連させて使用することを強調しています…バナキエヴィチは、基本的な問題は実際には行列因数分解、または彼が言うところの「分解」であるという点を理解していました。」^{[2]これは、} LR分解（左三角行列と右三角行列に因数分解する）と呼ばれることもあります。

A を正方行列とします。LU分解とは、適切な行および/または列の順序または順列を使用して、 Aを下三角行列Lと上三角行列U の2 つの因子に分解することを指します。

${\displaystyle A=LU.}$

下三角行列では対角線より上の要素はすべてゼロで、上三角行列では対角線より下の要素はすべてゼロです。たとえば、3 × 3 行列Aの場合、LU 分解は次のようになります。

${\displaystyle {\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}}={\begin{bmatrix}\ell _{11}&0&0\\\ell _{21}&\ell _{22}&0\\\ell _{31}&\ell _{32}&\ell _{33}\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}}.}$

行列の順序や順列が適切でないと、因数分解が実現しないことがあります。たとえば、（行列の乗算を展開することによって）は簡単に確認できます。 の場合、との少なくとも 1 つはゼロでなければなりません。これは、 LまたはU のいずれかが特異であることを意味します。 Aが特異でない（可逆である）場合は、これは不可能です。 これは手続き上の問題です。Aの行を並べ替えて、順列化された行列の最初の要素がゼロ以外になるようにするだけで、この問題は解消できます。後続の因数分解ステップでの同じ問題も同じ方法で解消できます。以下の基本手順を参照してください。 ${\textstyle a_{11}=\ell _{11}u_{11}}$${\textstyle a_{11}=0}$${\textstyle \ell _{11}}$${\textstyle u_{11}}$

LU 分解には行 (または列) の適切な順列で十分であることがわかります。部分ピボットによる LU 分解(LUP) は、行の順列のみによる LU 分解を指すことがよくあります。

${\displaystyle PA=LU,}$

ここで、LとUはそれぞれ下三角行列と上三角行列、P は置換行列で、これをAに左掛けするとAの行が並べ替えられます。すべての正方行列はこの形式で因数分解できることが分かっており、^[3]因数分解は実際には数値的に安定しています。^[4]このため、LUP分解は実際に役立つ手法となっています。

完全なピボットによる LU 分解には、行と列の両方の順列が含まれます。

${\displaystyle PAQ=LU,}$

ここで、L、U、Pは前と同じように定義され、QはAの列を並べ替える順列行列である。^[5]

下対角上（LDU）分解は、次の形式の分解である。

${\displaystyle A=LDU,}$

ここで、Dは対角行列、LとUは一三角行列です。つまり、 LとUの対角要素はすべて1 になります。

上記ではA が正方行列であることを条件としましたが、これらの分解はすべて長方行列にも一般化できます。 ^[6]その場合、LとD はどちらもAと同じ行数を持つ正方行列であり、U はAとまったく同じ次元を持ちます。上三角は、左上隅から始まる主対角線の下にゼロの要素のみを持つと解釈する必要があります。同様に、 Uのより正確な用語は、行列Aの行階段形式です。

次の 2 行 2 列の行列を因数分解します。

${\displaystyle {\begin{bmatrix}4&3\\6&3\end{bmatrix}}={\begin{bmatrix}\ell _{11}&0\\\ell _{21}&\ell _{22}\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}\\0&u_{22}\end{bmatrix}}.}$

この単純な行列のLU分解を見つける一つの方法は、線型方程式を単純に調べて解くことです。行列の乗算を展開 すると、

${\displaystyle {\begin{aligned}\ell _{11}\cdot u_{11}+0\cdot 0&=4\\\ell _{11}\cdot u_{12}+0\cdot u_{22}&=3\\\ell _{21}\cdot u_{11}+\ell _{22}\cdot 0&=6\\\ell _{21}\cdot u_{12}+\ell _{22}\cdot u_{22}&=3.\end{aligned}}}$

この連立方程式は不確定です。この場合、L行列とU行列の任意の 2 つの非ゼロ要素は解のパラメータであり、任意の非ゼロ値に設定できます。したがって、一意の LU 分解を見つけるには、L 行列とU行列に何らかの制限を課す必要があります。たとえば、下三角行列L が単位三角行列であることを都合よく要求して、その主対角要素がすべて 1 に設定されるようにすることができます。すると、連立方程式の解は次のようになります。

${\displaystyle {\begin{aligned}\ell _{11}=\ell _{22}&=1\\\ell _{21}&=1.5\\u_{11}&=4\\u_{12}&=3\\u_{22}&=-1.5\end{aligned}}}$

これらの値を上記のLU分解に代入すると、次の式が得られます。

${\displaystyle {\begin{bmatrix}4&3\\6&3\end{bmatrix}}={\begin{bmatrix}1&0\\1.5&1\end{bmatrix}}{\begin{bmatrix}4&3\\0&-1.5\end{bmatrix}}.}$

任意の正方行列は、 LUPおよびPLU分解が可能です。^[3]が可逆である場合、すべての主要な主小行列^[7]が非ゼロである場合に限り、LU (またはLDU ) 分解が可能です^[8] (たとえば、はLUまたはLDU分解 できません)。 が階数 の特異行列である場合、最初の主要な主小行列が非ゼロであればLU分解が可能ですが、その逆は成り立ちません。^[9]${\textstyle A}$${\textstyle A}$${\displaystyle {\begin{bmatrix}0&1\\1&0\end{bmatrix}}}$${\textstyle A}$${\textstyle k}$${\textstyle k}$

正方可逆行列がLDU （ LとUの対角要素がすべて1になる因数分解）を持つ場合、その因数分解は一意です。^[8]その場合、（または）の対角要素が1で構成されることを要求すれば、 LU因数分解も一意です。 ${\textstyle L}$${\textstyle U}$

一般に、正方行列は次のいずれかになります。 ${\displaystyle A_{n\times n}}$

1. 一意の LU 分解（前述のとおり）
2. 最初の ( n −1) 列のうち 2 列以上が線形従属であるか、最初の ( n −1) 列のいずれかが 0 である場合、LU 分解は無限に存在します。
3. 最初の( n −1)列がゼロ以外かつ線形独立であり、少なくとも1つの主要な主マイナーがゼロである場合、LU分解は行われません。

ケース3では、対角要素を に変更して ゼロ主マイナーを回避することでLU因数分解を近似することができます。 [ ^10]${\displaystyle a_{jj}}$${\displaystyle a_{jj}\pm \varepsilon }$

A が対称（または複素数の場合はエルミート）正定値行列である場合、 U がLの共役転置となるように配置することができます。つまり、A は次のように 書くことができます。

${\displaystyle A=LL^{*}.\,}$

この分解はコレスキー分解と呼ばれます。 が正定値であれば、コレスキー分解は存在し、一意です。さらに、コレスキー分解の計算は、他の LU 分解を計算するよりも 効率的で、数値的に安定しています。${\displaystyle A}$

任意の体上の（必ずしも可逆ではない）行列について、LU分解が成立するための正確な必要十分条件は既知である。条件は特定の部分行列の階数で表現される。LU分解を得るためのガウス消去法も、この最も一般的なケースに拡張されている。^[11]

LDU分解が存在し、それが一意である場合、元の行列Aの特定の部分行列の行列式の比に関して、L、D、Uの要素に対する閉じた（明示的な）式が存在します。^[12]特に、、およびの場合、は、番目の主部分行列と番目の主部分行列の比です。行列式の計算は計算コストが高いため、この明示的な式は実際には使用されません。 ${\textstyle D_{1}=A_{1,1}}$${\textstyle i=2,\ldots ,n}$${\textstyle D_{i}}$${\textstyle i}$${\textstyle (i-1)}$

以下のアルゴリズムは本質的にはガウス消去法の修正版である。このアルゴリズムを使用してLU分解を計算するには、低次の項を無視した浮動小数点演算が必要である。部分ピボットでは2次項のみが追加されるが、完全ピボットではそうではない。^[13]${\displaystyle {\tfrac {2}{3}}n^{3}}$

N × N行列が与えられた場合、 を行列 の元の変更されていないバージョンとして定義します。行列の括弧内の上付き文字 (例: )は、行列のバージョンです。行列 は、第 1 列のガウス消去法によって主対角線の下の要素が0 に消去された行列です。 ${\displaystyle A=(a_{i,j})_{1\leq i,j\leq N}}$${\displaystyle A^{(0)}}$${\displaystyle A}$${\displaystyle (0)}$${\displaystyle A}$${\displaystyle A^{(n)}}$${\displaystyle A}$${\displaystyle n}$

以下は、表記法を覚えるのに役立つ行列です (各行列は、行列内の 任意の実数を表します)。${\displaystyle *}$

${\displaystyle A^{(n-1)}={\begin{pmatrix}*&&&\cdots &&&*\\0&\ddots &&&&\\&\ddots &*&&&\\\vdots &&0&a_{n,n}^{(n-1)}&&&\vdots \\&&\vdots &a_{i,n}^{(n-1)}&*\\&&&\vdots &\vdots &\ddots \\0&\cdots &0&a_{i,n}^{(n-1)}&*&\cdots &*\end{pmatrix}}}$

このプロセスでは、行演算を使用して行列を徐々に変更し、主対角線の下の要素がすべてゼロになる行列にします。この間、となる2 つの別々の行列と を同時に作成します。 ${\displaystyle A}$${\displaystyle U}$${\displaystyle P}$${\displaystyle L}$${\displaystyle PA=LU}$

最終的な順列 行列を、行列 が行列 に変換されるときに、すべての同じ行が同じ順序で交換された単位行列として定義します。行列 の場合、行を交換して n 番目の列に必要な条件を提供することから始めます。たとえば、行を交換して部分的なピボットを実行したり、ガウス消去法を完了できるように主対角線のピボット要素を 0 以外の数値に設定したりする場合があります。 ${\displaystyle P}$${\displaystyle A}$${\displaystyle U}$${\displaystyle A^{(n-1)}}$${\displaystyle a_{n,n}}$

行列 では、以下のすべての要素をゼロに設定します(ここで は主対角線の n 番目の列の要素です)。以下の各要素を(ここで)と表記します。 をゼロに設定するには、各行を に設定します。この操作では、 となります。最初の列に対して行操作を実行すると、で表される上三角行列が得られます。 ${\displaystyle A^{(n-1)}}$${\displaystyle a_{n,n}^{(n-1)}}$${\displaystyle a_{n,n}^{(n-1)}}$${\displaystyle a_{n,n}^{(n-1)}}$${\displaystyle a_{i,n}^{(n-1)}}$${\displaystyle i=n+1,\dotsc ,N}$${\displaystyle a_{i,n}^{(n-1)}}$${\displaystyle row_{i}=row_{i}-(\ell _{i,n})\cdot row_{n}}$${\displaystyle i}$${\textstyle \ell _{i,n}:={a_{i,n}^{(n-1)}}/{a_{n,n}^{(n-1)}}}$${\displaystyle N-1}$${\displaystyle A^{(N-1)}}$${\displaystyle U}$

以前に計算した の値を下の式で 直接入力して、で表される下三角行列を作成することもできます。${\textstyle L}$${\displaystyle \ell _{i,n}}$

${\displaystyle L={\begin{pmatrix}1&0&\cdots &0\\\ell _{2,1}&\ddots &\ddots &\vdots \\\vdots &\ddots &\ddots &0\\\ell _{N,1}&\cdots &\ell _{N,N-1}&1\end{pmatrix}}}$

行列 が与えられている場合は、部分ピボットを実装して 1 行目と 2 行目を入れ替えることを選択します。その結果、行列と行列の最初の反復はそれぞれ次のようになります。行を入れ替えたら、 を実行することで、最初の列の主対角線の下にある要素を削除できます。これらの行を減算すると、行列から が導出されます。部分ピボットを実装しているため、導出した行列の 2 行目と 3 行目、および現在のバージョンの行列をそれぞれ入れ替えて、次のようになります。次に、 を実行することで、2 列目の主対角線の下にある要素を削除します。この行の減算後、の現在の反復では主対角線の下にゼロ以外の要素が存在しないため、この行の減算によって、最終行列 ( と表記) と最終行列が導出されます。対応する行も入れ替えると、最終行列が得られます。これで、これらの行列は という関係になります。 $${\displaystyle A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},}$$${\displaystyle A}$${\displaystyle P}$$${\displaystyle A^{(0)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\2&7&9\\\end{pmatrix}},\quad P^{(0)}={\begin{pmatrix}0&1&0\\1&0&0\\0&0&1\\\end{pmatrix}}.}$$$${\displaystyle {\begin{alignedat}{0}row_{2}=row_{2}-(\ell _{2,1})\cdot row_{1}\\row_{3}=row_{3}-(\ell _{3,1})\cdot row_{1}\end{alignedat}}}$$$${\displaystyle {\begin{alignedat}{0}\ell _{2,1}={\frac {0}{4}}=0\\\ell _{3,1}={\frac {2}{4}}=0.5\end{alignedat}}}$$${\displaystyle A^{(1)}}$$${\displaystyle A^{(1)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.}$$${\displaystyle P}$$${\displaystyle A^{(1)}={\begin{pmatrix}4&2&1\\0&6&8.5\\0&5&{\frac {22}{3}}\\\end{pmatrix}},\quad P^{(1)}={\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\\\end{pmatrix}}.}$$${\displaystyle row_{3}=row_{3}-(\ell _{3,2})\cdot row_{2}}$${\textstyle \ell _{3,2}={\frac {5}{6}}}$${\displaystyle A}$${\displaystyle A}$${\displaystyle U}$${\displaystyle P}$$${\displaystyle A^{(2)}=A^{(N-1)}=U={\begin{pmatrix}4&2&1\\0&6&8.5\\0&0&0.25\\\end{pmatrix}},\quad P={\begin{pmatrix}0&1&0\\0&0&1\\1&0&0\\\end{pmatrix}}.}$$${\displaystyle L}$$${\displaystyle L={\begin{pmatrix}1&0&0\\\ell _{3,1}&1&0\\\ell _{2,1}&\ell _{3,2}&1\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\0.5&1&0\\0&{\frac {5}{6}}&1\\\end{pmatrix}}}$$${\displaystyle PA=LU}$

この処理中に行を入れ替えなかった場合、各列に対して行操作を同時に実行できます。ここで、はN × N単位行列のn番目の列を転置ベクトルで置き換えたものです。言い換えると、下三角行列 ${\displaystyle n}$${\displaystyle A^{(n)}:=L_{n}^{-1}A^{(n-1)},}$${\displaystyle L_{n}^{-1}}$${\displaystyle {\begin{pmatrix}0&\dotsm &0&1&-\ell _{n+1,n}&\dotsm &-\ell _{N,n}\end{pmatrix}}^{\textsf {T}}.}$

${\displaystyle L_{n}^{-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&-\ell _{N,n}&&&1\end{pmatrix}}.}$

式を使用して最初の列のすべての行操作を実行することは、分解を見つけることと同等です 。となることを表します。 ${\displaystyle N-1}$${\displaystyle A^{(n)}:=L_{n}^{-1}A^{(n-1)}}$$${\displaystyle A=L_{1}L_{1}^{-1}A^{(0)}=L_{1}A^{(1)}=L_{1}L_{2}L_{2}^{-1}A^{(1)}=L_{1}L_{2}A^{(2)}=\dotsm =L_{1}\dotsm L_{N-1}A^{(N-1)}.}$$${\textstyle L=L_{1}\dotsm L_{N-1}}$${\displaystyle A=LA^{(N-1)}=LU}$

それでは、 のシーケンスを計算してみましょう。 には次の式があることがわかっています。 ${\displaystyle L_{1}\dotsm L_{N-1}}$${\displaystyle L_{i}}$

${\displaystyle L_{n}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&\ell _{N,n}&&&1\end{pmatrix}}}$

主対角線に 1 が含まれる 2 つの下三角行列があり、どちらの行列も主対角線の下に、もう 1 つの行列と同じ列に非ゼロ項目がない場合は、2 つの行列の積に、すべての非ゼロ項目を同じ位置に含めることができます。例:

${\displaystyle \left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&0&1&0&0\\63&0&0&1&0\\7&0&0&0&1\end{array}}\right)\left({\begin{array}{ccccc}1&0&0&0&0\\0&1&0&0&0\\0&22&1&0&0\\0&33&0&1&0\\0&44&0&0&1\end{array}}\right)=\left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&22&1&0&0\\63&33&0&1&0\\7&44&0&0&1\end{array}}\right)}$

最後に、を掛け合わせて、（前述のように）で表される融合行列を生成します。行列 を使用すると、次の式が得られます。${\displaystyle L_{i}}$${\textstyle L}$${\textstyle L}$${\displaystyle A=LU.}$

このアルゴリズムが機能するためには、各ステップで が必要であることは明らかです( の定義を参照)。この仮定がどこかの時点で失敗した場合は、続行する前にn番目の行をその下の別の行と交換する必要があります。これが、LU 分解が一般に のように見える理由です。 ${\displaystyle a_{n,n}^{(n-1)}\neq 0}$${\displaystyle \ell _{i,n}}$${\displaystyle P^{-1}A=LU}$

この手順で得られる分解はドゥーリトル分解であることに注意してください。つまり、 Lの主対角線は1のみで構成されています。行の倍数を追加して対角線の下にある要素を削除するのではなく、列の倍数を追加して主対角線の上にある要素を削除すると、 Uの主対角線が1であるCrout 分解が得られます。

与えられた行列Aの Crout 分解を生成する別の (同等の) 方法は、 Aの転置の Doolittle 分解を取得することです。実際、がこのセクションで紹介したアルゴリズムによって取得された LU 分解である場合、 および を取ると、はCrout 分解になります。 ${\textstyle A^{\textsf {T}}=L_{0}U_{0}}$${\textstyle L=U_{0}^{\textsf {T}}}$${\textstyle U=L_{0}^{\textsf {T}}}$${\displaystyle A=LU}$

Cormenら^[14]はLUP分解のための再帰アルゴリズムを説明した。

行列Aが与えられたとき、P _1を次のような順列行列と する。

${\displaystyle P_{1}A=\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\v&&A'&\\&&&\end{array}}\right)}$、

ここで、 Aの最初の列に非ゼロの要素がある場合は、それ以外の場合はP _1を単位行列とします。次に、の場合は、それ以外の場合は とします。 ${\textstyle a\neq 0}$${\textstyle c=1/a}$${\textstyle a\neq 0}$${\textstyle c=0}$

${\displaystyle P_{1}A=\left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\cv&&I_{n-1}&\\&&&\end{array}}\right)\left({\begin{array}{c|c}a&w^{\textsf {T}}\\\hline &\\0&A'-cvw^{\textsf {T}}\\&\end{array}}\right).}$

これでLUP分解を再帰的に見つけることができます。 とします。したがって ${\textstyle P'\left(A'-cvw^{\textsf {T}}\right)=L'U'}$${\textstyle v'=P'v}$

${\displaystyle \left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\0&&P'&\\&&&\end{array}}\right)P_{1}A=\left({\begin{array}{c|ccc}1&&0&\\\hline &&&\\cv'&&L'&\\&&&\end{array}}\right)\left({\begin{array}{c|ccc}a&&w^{\textsf {T}}&\\\hline &&&\\0&&U'&\\&&&\end{array}}\right),}$

これはAの LUP 分解です。

ランダム化アルゴリズムを使用して、LU分解の低ランク近似値を求めることができます。入力行列と望ましい低ランクが与えられると、ランダム化LUは、それぞれサイズがとの下台形行列と上台形行列を返し、高い確率で となります。ここで、 はアルゴリズムのパラメータに依存する定数で、 は入力行列の 番目の特異値です。^[15]${\textstyle A}$${\textstyle k}$${\textstyle P,Q}$${\textstyle L,U}$${\textstyle m\times k}$${\textstyle k\times n}$${\textstyle \left\|PAQ-LU\right\|_{2}\leq C\sigma _{k+1}}$${\textstyle C}$${\textstyle \sigma _{k+1}}$${\textstyle (k+1)}$${\textstyle A}$

2つのn次行列をM ( n )の時間で乗算できる場合（ M ( n )≥n a ^、 a >2の場合）、LU分解はO( M ( n ))の時間で計算できます。^[16]これは、たとえば、Coppersmith–Winogradアルゴリズムに基づくO(n2.376)アルゴリズムが存在することを^意味します。

大規模な疎行列を因数分解するための特別なアルゴリズムが開発されています。これらのアルゴリズムは疎因数LとU を見つけようとします。理想的には、計算コストは​​行列のサイズではなく、非ゼロのエントリの数によって決まります。

これらのアルゴリズムは、行と列を自由に交換して、フィルイン (アルゴリズムの実行中に初期のゼロからゼロ以外の値に変化するエントリ) を最小限に抑えます。

埋め込みを最小化する順序付けの一般的な処理は、グラフ理論を使用して対処できます。

行列形式の線形方程式系が与えられる

${\displaystyle A\mathbf {x} =\mathbf {b} ,}$

Aとbが与えられたとき、方程式をxについて解きたいとします。となるAの LUP 分解がすでに得られていると仮定します。 ${\textstyle PA=LU}$${\textstyle LU\mathbf {x} =P\mathbf {b} }$

この場合、解決は次の 2 つの論理的なステップで実行されます。

1. まず、 yについて方程式を解きます。${\textstyle L\mathbf {y} =P\mathbf {b} }$
2. 次に、 xについて方程式を解きます。${\textstyle U\mathbf {x} =\mathbf {y} }$

どちらの場合も、三角行列 ( LとU )を扱っており、ガウス消去法を使用せずに前方置換と後方置換によって直接解くことができます (ただし、 LU分解自体を計算するには、このプロセスまたは同等のプロセスが必要です)。

上記の手順は、異なるbについて方程式を複数回解くために繰り返し適用できます。この場合、毎回ガウス消去法を使用するよりも、行列Aの LU 分解を 1 回実行してから、異なるbについて三角行列を解く方が高速 (かつ便利) です。行列LとU は、ガウス消去法のプロセスを「エンコード」したものと考えることができます。

行列のサイズが の場合、線形方程式系を解くコストは、浮動小数点演算 程度です。これにより、ハウスホルダー反射を使用する場合に浮動小数点演算程度のコストがかかるQR 分解に基づくアルゴリズムよりも 2 倍高速になります。このため、LU 分解が通常は好まれます。^[17]${\textstyle {\frac {2}{3}}n^{3}}$${\textstyle A}$${\textstyle n}$${\textstyle {\frac {4}{3}}n^{3}}$

連立方程式を解く場合、b は通常、行列Aの高さに等しい長さのベクトルとして扱われます。ただし、逆行列では、ベクトルbの代わりに行列Bを使用します。ここで、 Bはn行p列の行列であるため、行列X (これもn行p列の行列)を求めます。

${\displaystyle AX=LUX=B.}$

行列Xの各列について解くには、先に示したのと同じアルゴリズムを使用することができます。ここで、Bがサイズnの単位行列であると仮定します。結果XはAの逆行列でなければなりません。^[18]

正方行列AのLUP分解が与えられれば、 Aの行列式は次のように簡単に計算できる。 ${\textstyle A=P^{-1}LU}$

${\displaystyle \det(A)=\det \left(P^{-1}\right)\det(L)\det(U)=(-1)^{S}\left(\prod _{i=1}^{n}l_{ii}\right)\left(\prod _{i=1}^{n}u_{ii}\right).}$

^{2 番目の方程式は、三角行列の行列式が単純に対角要素の積であり、置換行列の行列式が (-1) S}に等しい( Sは分解における行交換の数) という事実から導かれます。

完全なピボットによる LU 分解の場合、行と列の交換の総数を Sとすると、上記の式の右側の項も等しくなります。${\textstyle \det(A)}$

同じ方法は、P を単位行列に設定することで LU 分解にも簡単に適用できます。

モジュールmlu
暗黙的 None整数、パラメータ:: SP = Kind ( 1 d0 ) ! I/O 実数精度を設定非公開      公開luban 、lusolve      サブルーチンluban ( a 、tol 、g 、h 、ip 、condinv 、detnth )を含む! LU 分解により、平方 A=LU=G^TH となるような上三角形 L=G^T、および U=H を計算します。 部分ピボット IP(:) は、最新の加算です。! 通常は平方 A に対して使用しますが、RHS a が既にわかっている場合は、 ! 入力 (A|a)^T は (L|y^T)^T となり、L^Tx=y の x は Ax=a の解となります。Real ( SP )、Intent ( In ) :: a (:, :) !入力行列 A(m,n), n<=m Real ( SP ), Intent ( In ) :: tol ! ゼロに近いピボットの許容値Real ( SP ), Intent ( Out ) :: g ( size ( a , dim = 1 ), size ( a , dim = 2 )), & ! L(m,n) h ( size ( a , dim = 2 ), size ( a , dim = 2 )), & ! U(n,n) condinv , & ! 1/cond(A), 特異な A の場合は 0 detnth ! sign*Abs(det(A))**(1/n) Integer , Intent ( Out ) :: ip ( size ( a , dim = 2 )) !列の順列整数:: k , n , j , l , isig実数( SP ) :: tol0 , pivmax , pivmin , piv ! n =サイズ      
            
      
 

 




                  
                     
                 
                                       
                                      
                                         
                 
         
               
               

            ( a , dim = 2 ) tol0 = max ( tol , 3._SP * epsilon ( tol0 )) ! tol=0 の場合はデフォルトを使用します。 ! 長方形 A と G は、次の条件で許可されます:   If ( n > size ( a , dim = 1 ) 。または。n < 1 ) Stop 90 Forall ( k = 1 : n ) ip ( k ) = k h = 0._SP g = 0._SP isig = 1 detnth = 0._SP ! k = 1 、n とします。 Banachiewicz (1938) 式 (2.3) h ( k 、ip ( k : n )) = a ( k 、ip ( k :)) - Matmul ( g ( k 、: k - 1 )、h (: k - 1 、ip ( k :))) !行ピボットを検索しますj = ( Maxloc ( Abs ( h ( k , ip ( k :))), dim = 1 ) + k - 1 ) If ( j /= k ) Then ! 列 j と k を入れ替えますisig = - isig ! 順列のため Det(A) の符号を変更しますl = ip ( k ) ip ( k ) = ip ( j ) ip ( j ) = l End If piv = Abs ( h ( k , ip ( k ))) If ( k 
               

                  
             
         
         
           
           

             

                  

                
                     
                       
                 
                 
                 
            
               
              == 1 ) Then ! 条件を初期化pivmax = piv pivmin = piv Else ! 条件を調整pivmax = Max ( piv 、pivmax ) pivmin = Min ( piv 、pivmin ) End If                            If ( piv < tol0 ) Then ! 特異行列isig = 0 pivmax = 1._SP Exit             Else ! ピボットの寄与を Det(A) の符号と値に考慮If ( h ( k 、ip ( k )) < 0._SP ) isig = - isig detnth = detnth + Log ( piv ) End If ! Banachiewicz (1938) 式。 (2.4) g ( k + 1 :, k ) = ( a ( k + 1 :, ip ( k )) - Matmul ( g ( k + 1 :, : k - 1 ), h (: k - 1 , ip ( k )))) & / h ( k , ip ( k )) g ( k , k ) = 1._SP End Do detnth = isig * exp ( detnth / n ) condinv = Abs ( isig ) * pivmin / pivmax ! 以下のコメントを解除して、平方 A(n,n) をテストします       。 Print *, '|AP-LU| ',Maxval (Abs(a(:,ip(:))-Matmul(g, h(:,ip(:))))) End Subroutine luban Subroutine lusolve ( l , u , ip , x ) !三角因数 LU=A を使用して Ax=a 系を解く             
                 
                 
                             
                  
                  
            
     
                 
                 
               
 
                      
                   
               

                  
                
              
         
               
                


      
      
      

         Real ( SP ), Intent ( In ) :: l (:, :) ! 下三角行列 L(n,n) Real ( SP ), Intent ( In ) :: u (:, :) ! 上三角行列 U(n,n) Integer , Intent ( In ) :: ip (:) ! 列の順列 IP(n) Real ( SP ), Intent ( InOut ) :: x (:, :) ! m セットの入力に対して X(n,m) !右辺は出力未知数に置き換えられますInteger :: n , m , i , j n = size ( ip ) m = size ( x , dim = 2 ) If ( n < 1.または. m < 1.または. Any ([ n , n ] /= shape ( l )).または. Any ( shape ( l ) /= shape ( u )).または. & n /= size ( x , dim = 1 )) Stop 91 Do i = 1 , m Do j = 1 , n x ( j , i ) = x ( j , i ) - dot_product ( x (: j - 1 , i ), l ( j ,: j - 1 )) End Do             Do j = n , 1 , - 1 x ( j , i ) = ( x ( j           
                    
                      
                 
            
         
              
        
           
            
           
             
            
               
                 
            
    
                 , i ) - dot_product ( x ( j + 1 :, i ), u ( j , ip ( j + 1 :)))) / & u ( j , ip ( j )) End Do          End Do       End サブルーチンlusolve End モジュールmlu