LU分解

Type of matrix factorization

数値解析と線型代数において、下三角行列と上三角行列の積として行列を因数分解する(行列の乗算と行列分解を参照)。 積には、置換行列が含まれることもある。 LU分解は、ガウス消去法の行列形式とみなすことができる。 コンピュータは通常、LU 分解を使用して 2 次線形方程式を解き、また、行列の逆行列を求めるときや行列の行列式を計算するときの重要なステップでもある。 また、 LR分解 (左三角行列と右三角行列に因数分解する)と呼ばれることもある。 LU 分解は、 1938 年にポーランドの天文学者タデウシュ バナチエヴィチによって導入された^{[1] 。}彼は、積方程式を初めて書き記した(彼の代替だが同等の行列表記の最後の形式は として現れる)。 ${\displaystyle LU=A=h^{T}g}$ ${\displaystyle g\times h.}$

定義

ウォルシュ行列のLDU分解

A を正方行列とします。LU分解とは、 Aを下三角行列Lと上三角行列Uの2つの因子の積として表現することです：A = LU 。ゼロ除算や丸め誤差の制御不能な増加を防ぐためにAを事前に並べ替えないと分解できない場合があり、そのため代替表現はPAQ = LUとなります。ここで、正式な表記法では、置換行列因子PとQ はAの行（または列）の置換を示します。理論上はP（またはQ ）は単位行列の行（または列）の置換によって得られますが、実際には対応する置換がAの行（または列）に直接適用されます。

n辺の行列Aには係数がありますが、結合された 2 つの三角行列にはn ( n + 1)個の係数が含まれるため、行列LUのn個の係数は独立ではありません。通常の慣例により、 L を単三角行列、つまりn 個の主対角要素がすべて 1 になるように設定します。ただし、代わりにU行列を単三角行列に設定すると、行列積の 転置後に同じ手順になります (行列転置の特性を参照)。 転置後、U ^TはBの 下三角因子になり、L ^Tは上単三角因子になります。これはまた、行の操作 (ピボットなど) は転置行列の列の操作と同等であり、一般に行アルゴリズムと列アルゴリズムのどちらを選択しても利点がないことも示しています。 ${\displaystyle n^{2}}$ $${\displaystyle B=A^{\textsf {T}}=(LU)^{\textsf {T}}=U^{\textsf {T}}L^{\textsf {T}}.}$$

下三角行列では主対角線より上側の要素はすべてゼロであり、上三角行列では主対角線より下側の要素はすべてゼロです。例えば、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 つは0 でなければなりません。これは、 LまたはUのいずれかが特異 であることを意味します。 Aが非特異（可逆）である場合、これは不可能です。 演算に関しては、Aの最初の列の残りの要素を 0 にすること/除去することはでのの除算を必要としますが、 が 0 の場合は不可能です。 これは手続き上の問題です。これは、並べ替えられた行列の最初の要素が 0 以外になるようにAの行を並べ替えるだけで解消できます。後続の因数分解ステップで発生する同じ問題も、同じ方法で解消できます。丸め誤差/小さい数による除算に対する数値的安定性を確保するには、絶対値の大きい を選択することが重要です（ピボットを参照）。 ${\textstyle a_{11}=\ell _{11}u_{11}}$ ${\textstyle a_{11}=0}$ ${\textstyle \ell _{11}}$ ${\textstyle u_{11}}$ ${\textstyle a_{21},a_{31}}$ ${\textstyle a_{11}}$ ${\textstyle a_{11}}$

LU 再帰を通じて

上記の3×3行列の例では、関係する行列の最上行と最左列の行列積が、LU法の成功に特別な役割を果たすことが示されています。連続する行列を でマークし、これらの行と列を他の行と区別するように行列積を記述します。その際、ブロック行列記法を使用します。例えば、は通常の数、は行ベクトル、 は列ベクトル、は最上行と最左列のない行列の部分行列です。次に、 を ブロック行列積で置き換えます。つまり、行列ブロックを通常の数、つまり行×列のように乗算できますが、その成分は部分行列であり、スカラーまたはベクトルに簡約される場合もあります。したがって、 は各成分に数 を乗算した後に得られるベクトルを表します。はベクトルの外積、つまり最初の列が、次の列がなどである行列であり、のすべての成分について以下同様になり 、は の部分行列の積です。 ${\displaystyle (0),\;(1),\dots }$ ${\displaystyle A\equiv A^{(0)}=L^{(0)}U^{(0)}}$ ${\displaystyle a\equiv a_{11}}$ ${\displaystyle {\bf {w}}^{\textsf {T}}\equiv (a_{12},a_{13})^{\textsf {T}}}$ ${\displaystyle {\bf {v}}=(a_{21},a_{31})}$ ${\displaystyle A'}$ ${\displaystyle A^{(0)}}$ ${\displaystyle A^{(0)}=L^{(0)}U^{(0)}}$ ${\displaystyle u{\bf {l}}}$ ${\displaystyle {\bf {l}}}$ ${\displaystyle u}$ ${\displaystyle {\bf {lu}}^{\textsf {T}}}$ ${\displaystyle {\bf {l,u}}}$ ${\displaystyle u_{12}{\bf {l}}}$ ${\displaystyle u_{13}{\bf {l}}}$ ${\displaystyle {\bf {u}}}$ ${\displaystyle L^{(1)}U^{(1)}}$ ${\displaystyle L^{(0)},\;U^{(0)}}$ $${\displaystyle {\begin{aligned}\left({\begin{array}{c|c}a&{\bf {w}}^{\textsf {T}}\\\hline \\[-0.5em]{\bf {v}}&\quad A'\quad \\[-0.5em]\\\end{array}}\right)&=\left({\begin{array}{c|c}{\rm {1}}&{\bf {0}}^{\textsf {T}}\\\hline \\[-0.5em]{\bf {l}}&\quad L^{(1)}\quad \\[-0.5em]\\\end{array}}\right)\;\left({\begin{array}{c|c}u&{\bf {u}}^{\textsf {T}}\\\hline \\[-0.5em]{\bf {0}}&\quad U^{(1)}\\[-0.5em]\\\end{array}}\right)\\&=\left({\begin{array}{c|c}u&{\bf {u}}^{\textsf {T}}\\\hline \\[-0.5em]u{\bf {l}}&\quad {\bf {lu}}^{\textsf {T}}+L^{(1)}U^{(1)}\\[-0.5em]\\\end{array}}\right)\end{aligned}}}$$

最初と最後の行列が等しいことから、最終的に が導かれ、行列は に更新/置換されます 。ここで重要な点に留意してください。を と同じように繰り返し処理しても問題ありません。 の次元がn  ×  nの場合、このようなn − 1 回のステップの後、すべての列が三角行列の対角部分を形成し、すべてのピボットと行が組み合わされて上三角行列 が形成されます。上記の例ではn = 3なので、2 つのステップだけで十分です。 ${\displaystyle u=a}$ ${\displaystyle {\bf {u=w}}}$ ${\displaystyle {\bf {l}}={(1/a)}{\bf {v}}}$ ${\displaystyle A'}$ ${\displaystyle A^{(1)}\equiv L^{(1)}U^{(1)}=}$ ${\displaystyle A'-{\bf {lu}}^{\textsf {T}}}$ ${\displaystyle A^{(1)}}$ ${\displaystyle A^{(0)}}$ ${\displaystyle A}$ ${\displaystyle {\bf {v}}}$ ${\displaystyle L}$ ${\displaystyle a}$ ${\displaystyle {\bf {w}}^{\textsf {T}}}$ ${\displaystyle U}$

上記の手順は、連続する部分行列の対角ピボット要素がどのステップでもゼロになることはないことを示しています。これを避けるには、列または行を入れ替えて、 が非ゼロになるようにします。このような順列置換を伴う手順は、 LUP（ピボット分解）と呼ばれます。 ${\displaystyle a}$ ${\displaystyle a}$

列の置換は行列積に対応し、ここでは置換行列、すなわち同じ列の置換後の単位行列である。すべてのステップの後、このようなLUP分解は に適用される。本計算スキームとCormenら^{[2]の類似のスキームは、}再帰アルゴリズムの例である。これらはLU分解の2つの一般的な性質を示している。 ${\displaystyle AQ^{(0)}}$ ${\displaystyle Q^{(0)}}$ ${\displaystyle I}$ ${\displaystyle AQ^{(0)}\cdots Q^{(n-1)}\equiv AQ=LU}$

1. 各ステップで方向転換する必要性
2. L 行列とU行列の最終値は、ステップごとに 1 行または 1 列ずつ段階的に取得されます。

再帰アルゴリズムは代数演算の観点からはそれほどコストがかからないものの、各ステップでAのほとんどの要素を更新して保存する必要があるため、実用上の欠点があります。計算順序を変更することで、中間値の保存を省略できることがわかります。

部分ピボットによるLU分解

既知の異常なケースを除いて、数値的に安定した LU 分解には、行 (または列) の絶対最大ピボットa ₁₁を選択するための行 (または列) の適切な順列化で十分であることが判明しています。これは部分ピボットによる LU 分解(LUP) と呼ばれます。ここで、LとUは再び下三角行列と上三角行列、PとQは対応する順列行列で、 Aにそれぞれ左乗算と右乗算すると、 Aの行と列の順序が変わります。すべての正方行列はこの形式で因数分解でき、^[3]因数分解は実際には数値的に安定していることが判明しています。^[4]このため、LUP 分解は実際に役立つ手法となっています。 $${\displaystyle PA=LU,\quad (AQ=LU),}$$

各ステップにおけるルークピボットと呼ばれる変種では、ルークがチェス盤上を列、行、列と移動していくのと同じように、最大​​要素を探索し、行と列の両方でピボット最大値に到達するまで続けます。ランダム要素からなる大規模な行列の場合、各ステップにおける演算コストは​​、フルピボットの場合のコストの2乗とは異なり、行列の辺の長さに比例する部分ピボットと同様に高いことが証明されています。

フルピボットによるLU分解

フルピボットによるLU分解では、行と列の両方の順列置換を行って、部分行列全体の絶対最大要素を見つけます。ここで、 L、U、Pは前と同じように定義され、QはAの列を並べ替える順列行列です。^[5] $${\displaystyle PAQ=LU,}$$

下対角上（LDU）分解

下対角上分解(LDU) は、形式の分解です。ここで、Dは対角行列、LとUは一三角行列です。つまり、LとUの対角線上のすべての要素は 1 つです。 $${\displaystyle A=LDU,}$$

長方形行列

上記ではA が正方行列であることを前提としましたが、これらの分解はすべて長方行列にも一般化できます。 ^[6]その場合、LとD はどちらもAと同じ行数を持つ正方行列となり、U はAと全く同じ次元を持ちます。「上三角行列」とは、左上隅から始まる主対角線の下に要素が0個しかない行列と解釈されます。同様に、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 \left\{{\begin{alignedat}{4}\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{alignedat}}\right.}$$

この連立方程式は不確定である。この場合、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}}.}$$

存在と唯一性

正方行列

任意の正方行列A は、LUP 分解と PLU 分解が可能です。^[3] Aが可逆行列である場合、そのすべての主小行列が非ゼロのときに限り、LU 分解（または LDU 分解）が可能です^{[7] [8]}（例えば、 はLU 分解も LDU 分解もできません）。A が k 階の特異行列である場合、最初のk個の主小行列が非ゼロのときに LU 分解が可能ですが、その逆は成り立ちません。^[9] $\left[{\begin{smallmatrix}0&1\\1&0\end{smallmatrix}}\right]$

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

一般に、任意の正方行列A _{n × n}は次のいずれかになります。

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

ケース3では、対角要素a _ijをa _ij ± εに 変更することでLU分解を近似し、ゼロ主マイナーを回避することができます。^[10]

対称正定値行列

A が対称正定値行列（複素行列の場合はエルミート行列）である場合、 UはLの共役転置行列となるように整理できる。つまり、Aは次のように 書ける。 $${\displaystyle A=LL^{*}\,.}$$

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

一般的な行列

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

アルゴリズム

閉じた式

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

ガウス消去法を使う

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

一般化された説明

表記

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

以下は、表記法を覚えるのに役立つ行列です（各∗ は行列内の任意の実数を表します）。

${\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}}}$

手順

この過程では、行演算を用いて行列Aを徐々に修正し、主対角線より下の要素がすべてゼロとなる行列Uを作成します。この過程で、 PA = LUとなる2つの別々の行列PとLを同時に作成します。

最終的な順列行列 Pは、行列 Aを行列Uに変換する際、すべての同じ行を行列Aと同じ順序で入れ替えた単位行列として定義します。行列A ^{( n −1)}の場合、まずn列目の望ましい条件を満たすために行を入れ替えます。例えば、部分ピボットを実行するために行を入れ替えたり、ガウス消去法を完了できるように主対角線上のピボット要素a _n,nを非ゼロの値に設定するために行を入れ替えたりします。

行列A ^{( n −1)}について、以下のすべての要素をゼロに設定します( は主対角線のn番目の列の要素です)。以下の各要素を( i = n +1, ... , N )と表記します。ゼロに設定するには、各行iについて、行_i =行_i − ( ℓ _i,n )⋅行_nと設定します。この操作では、 です。最初のN − 1列に対して行操作を実行すると、 Uで表される上三角行列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 a_{i,n}^{(n-1)}}$ ${\textstyle \ell _{i,n}:={a_{i,n}^{(n-1)}}/{a_{n,n}^{(n-1)}}}$

以前に計算したℓ _i,nの値を以下の式で 直接入力して、 Lで表される下三角行列を作成することもできます。

$${\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行目を入れ替え、行列AとP行列の最初の反復がそれぞれ次のようになるようにします。行を入れ替えたら、1列目の主対角線の下の要素を次のよう に削除します。これらの行を減算すると、 A ⁽¹⁾から次の行列が導かれます。 $${\displaystyle A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},}$$ $${\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)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.}$$

部分ピボットを実装しているので、導出した行列の2行目と3行目、および現在のP行列の行を交換して、次の式を得ます 。ここで、行₃ =行₃ − ( ℓ _3,2 )⋅行₂を実行して、2列目の主対角線の下の要素を削除します。ℓ _3,2 = ⁠ $${\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}}.}$$5/6⁠ 。この行の減算後、 Aの現在の反復では主対角線の下に非ゼロ要素が存在しないため、この行の減算によって最終的なA行列（ Uと表記）と最終的なP行列が導出されます。対応する行を入れ替えると、最終的なL行列が得られます。 $${\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={\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}}}$$

これらの行列には、 PA = LUという関係があります。

行が交換されない場合の関係

このプロセス中に行をまったく交換しなかった場合は、 を設定することで各列nに対して同時に行操作を実行できます。ここで、 はN  ×  N単位行列で、n番目の列は転置ベクトル(0 ⋯ 0 1 − ℓ _{n +1, n}   ⋯  − ℓ _N、n ) ^Tに置き換えられます。 ${\displaystyle A^{(n)}:=L_{n}^{-1}A^{(n-1)},}$ ${\displaystyle L_{n}^{-1}}$

言い換えれば、下三角行列 $${\displaystyle L_{n}^{-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&-\ell _{N,n}&&&1\end{pmatrix}}.}$$

式を使用して最初のN − 1列のすべての行操作を実行することは_、分解を見つけることと同じです。L = L 1 ⋯ L _{N −1}と表記し、A = LA ^{( N −1)} = LUとします。 ${\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)}.}$$

さて、L ₁ ⋯ L _{N −1}のシーケンスを計算してみましょう。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)}$

最後に、L ₁を掛け合わせて、 Lで表される融合行列を生成します（前述の通り）。行列Lを用いて、 A = LUを得ます。

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

LUバナチエヴィッチ分解

行列UとLのそれぞれ3行目と3列目を求める Banachiewicz LU アルゴリズムの動作を示す図。関係する行列は、その内容を示す四角形の上に名前が付けられています。行列の積と減算は、太枠で囲まれた要素に対してのみ適用されます。緑色で塗りつぶされた細枠は、前の段階で既にわかっている値を示します。青色の四角は、行列UとLにおける結果の格納場所を示します。

Banachiewicz (1938) の LU 分解アルゴリズムは、プログラムされた電子計算機の出現に先行していましたが、インデックス交換、転置、列ごとの乗算はほとんどのプログラミング言語に組み込まれたネイティブ機能であり、実際の実行をほとんど遅延させることなくコンパイラだけで処理できるため、コードに直接実装できるようになりました。 Banachiewicz が使用した独特の行列表記法によって、彼は行列を列ごとに乗算することができました。これは、定規を行列の次の行にスライドさせることで連続する因子を明らかにすることができるため、機械計算には便利な機能でした。 ただし、人間の読者にとっては、彼の方程式を標準の行列表記に変換するのが最適です。 完全な行列Aから三角行列UおよびLの計算を得るには、まず、 A の一番上の行と左端の列をそれぞれ行列UとLの対応する位置にコピーします。 Lの既知の単位対角要素は、処理全体を通して保存も使用もされません。

この図は、前の段階が既に完了していると仮定して、3行目と3列目の計算を示しています。関係する行列は、その内容を示す四角形の上に名前が付けられています。行列の積と減算は、太枠で囲まれた要素にのみ適用されます。緑色で塗りつぶされた細枠は、前の段階から既にわかっている値を示します。青色の四角は、U行列とL行列における結果の格納場所を示します。各段階で、 Lの結果要素をUの主対角線上の対応するピボット要素で割る必要があることに注意してください。これはLの左端の列にも適用されます。

第3段階の完了後、行列Aの関連する要素は使用されなくなり、以前の段階の要素も使用されなくなることに注意してください。これにより、これらの要素をUとLの結果値に置き換えることが可能になります。つまり、LU分解をその場で実行し、 Lの単位対角線を除くA全体をUとLに置き換えます。Banachiewicz LUアルゴリズムは、 Uの新しく計算された行から絶対最大ピボットを選択し、その後、その列を主対角線上に配置するように入れ替える部分ピボットに適しています。詳細は、添付のFortran90コードを調べることで理解できます。

すべての部分ピボット LU アルゴリズムは、ほぼ同じ量の順序操作を必要とします。ここで、nはAの行数または列数です。 ${\textstyle O\left({2 \over 3}n^{3}\right)}$

LU Crout分解

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

与えられた行列AのCrout分解を生成する別の（同等の）方法は、 Aの転置行列のDoolittle分解を得ることである。実際、A ^T = L ₀ U ₀がこの節で示したアルゴリズムによって得られるLU分解であるとすると、L = U^T
₀そしてU = L^T
₀すると、 A = LUは Crout 分解であることがわかります。

ランダム化アルゴリズム

ランダム化アルゴリズムを用いてLU分解の低階数近似を求めることが可能である。入力行列Aと所望の低階数kが与えられると、ランダム化LUは、それぞれサイズm×kとk×nの順列行列P、Qと下/上台形行列L、Uを返す。この場合、高い確率で‖PAQ − LU‖2≤Cσk  +  1が成立 する 。ここで、 Cはアルゴリズムのパラメータに依存する_定数であり_、 σk + 1は入力行列_Aの(k+1)番目の特異値である。 [ 14 ^]

理論的な複雑さ

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

スパース行列分解

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

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

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

アプリケーション

線形方程式を解く

行列形式の線形方程式系が与えられる $${\displaystyle A\mathbf {x} =\mathbf {b} ,}$$

Aとbが与えられたとき、xについて方程式を解きたいとします。AのLUP 分解がPA = LUとなるように既に得られていると仮定します。つまり、LU x = P bとなります。

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

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

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

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

線形方程式を解くコストはおよそ⁠2/3⁠行列A のサイズがnの場合、 n^{回の浮動小数点演算で済みます。これは、約}⁠のコストがかかるQR 分解に基づくアルゴリズムの 2 倍の速度です4/3⁠ nハウスホルダー反射^{を用いる場合、 3回の}浮動小数点演算が必要となる。このため、LU分解が通常は好まれる。^[16]

行列の反転

連立方程式を解く場合、bは通常、行列Aの高さに等しい長さのベクトルとして扱われます。しかし、逆行列法ではベクトルbの代わりに行列B（Bはn  ×  p行列）を用います。つまり、行列X（これもn  ×  p行列）を求めようとするのです。 $${\displaystyle AX=LUX=B.}$$

行列Xの各列について、先ほど示したのと同じアルゴリズムを用いることができる。ここで、Bがn次元の単位行列（I n ）であるとする_。したがって、結果のX はAの逆行列でなければならない。^[17]

行列式の計算

正方行列AのLUP分解A = P ⁻¹ LUが与えられれば、 Aの行列式は次のように簡単に計算できる。 $${\displaystyle {\begin{aligned}\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).\end{aligned}}}$$

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

フルピボットによるLU分解の場合、行と列の交換の総数をSとすると、 det( A )は上記の式の右側の項と等しくなります。

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

歴史

LU分解：LU因子とその積をBanachiewicz(1938)の行列表記法で表す

LU分解は、例えばラルストンによって記述された線形方程式の消去法と関連している。^[18] N個の未知数を持つN個の線形方程式を消去法で解くことは、古代中国で既に知られていた。 [ ^19]ガウス以前にも、ユーラシア大陸では多くの数学者がこの手法を実践し、完成させていたが、この手法が小学校レベルでしか扱われなくなったため、詳細な記述を残した数学者はほとんどいなかった。したがって、「ガウスの消去法」という名称は、複雑な歴史を便宜的に略したに過ぎない。

ポーランドの天文学者タデウシュ・バナキェヴィチは1938年にLU分解を導入した。^[1]バナキェヴィチについて、ポール・ドワイヤーは次のように述べている。^{[20] [ページが必要]}

ガウスとドゥーリトルは、消去法を対称方程式にのみ適用したようです。より最近の著者、例えばエイトキン、バナチェヴィッツ、ドワイヤー、クラウトなどは、非対称問題との関連でこの方法、あるいはその変種を用いることを強調しています。バナチェヴィッツは、基本的な問題は実際には行列分解、あるいは彼が「分解」と呼んだものの問題であるという点を理解していました。

— ポール・ドワイヤー『線形計算』（1951年）

バナチエヴィチ^[1]は、行列を用いて消去法を考察した最初の人物であり、そのようにしてLU分解を定式化した。これは彼の図解で示されている。彼の計算は一般的な行列の計算に従っているが、記法は異なっており、彼は1つの因子を転置して書くことを好んだ。これは、定規を両方の行の連続する行にスライドさせることで（算術計を用いて）、機械的に列ごとに掛け算できるようにするためである。添え字の順序を入れ替えたことにより、彼の式は現代の記法では以下のように表される。 $${\displaystyle {\begin{aligned}{\mathbf {x} }\cdot IA'={\mathbf {0} }&\rightarrow A'{\mathbf {x} =0}\equiv (A|{\mathbf {l} }){\mathbf {x} },\\A=G\cdot H&\rightarrow A^{T}=G^{T}H,\end{aligned}}}$$

ここでIA → A ^T ; x ≡ [ x ₁ , ... , x _n , −1] ; A ′はA を最後の列で拡張したものを参照し、xの最後の要素は−1です。 LU因子の行と列を再帰的に計算する行列式は、Banachiewiczの論文の残りの部分で式(2.3)と(2.4)として示されています（以下の§ Fortran90コード例を参照）。 Banachiewiczのこの論文には、それぞれ非対称行列と対称行列のLU因子とR ^T R因子の導出が含まれています。後の出版物では彼の名前がコレスキー分解の再発見とのみ結び付けられる傾向があるため、これらは混同されることがあります。バナキェヴィチ自身は、翌年には占領軍による迫害に苦しみ、ザクセンハウゼン強制収容所で3か月間過ごし、そこから釈放された後、協力者で共に囚人だったアントニ・ヴィルクを列車から運び出し、ヴィルクが1週間後に衰弱して亡くなったことから、何も行動を起こさなかったことは許されるだろう。

コード例

Fortran90コード例

モジュールmlu暗黙的 None整数、パラメータ:: SP = Kind ( 1 d0 ) ! I/O を実数精度に設定Private  Public luban 、lusolve サブルーチンluban ( a 、tol 、g 、h 、ip 、condinv 、detnth )を含む! Banachiewicz (1938、以降 B38) による LU 分解法は、! 三角形 L=G^T、U=H で正方形 B=A^T=G^TH=LU となるものを計算し! 列の置換 IP(:) による部分ピボットは現代の加算です。! コード内では、a、g は B38 A^T および G^T に対応するので、a=gh が成り立ちます。! ! 通常は正方形 A に使用しますが、RHS l については既に知られています! (A|l)^T を入力すると (L|y^T)^T が生成されます。ここで、L^Tx=y の x は Ax=l の解です。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) Real ( SP ), Intent ( Out ) :: h ( size ( a , dim = 2 ), size ( a , dim = 2 )) ! U(n,n) ! U 列は並べ替えられることに注意してくださいReal ( SP ), Intent ( Out ) :: condinv ! 1/cond(A), 特異値Aの場合は0 Real ( SP ), Intent ( Out ) :: detnth ! sign*Abs(det(A))**(1/n) Integer , Intent ( Out ) :: ip                                                                       ( size ( a , dim = 2 )) ! 列の順列! Integer :: k , n , j , l , isig Real ( SP ) :: tol0 , pivmax , pivmin , piv ! n = size ( a , dim = 2 ) tol0 = Max ( tol , 3._SP * epsilon ( tol0 )) ! tol=0 の場合はデフォルトを使用!長方形 A および G は、次の条件で許可されます。If ( n > size ( a 、dim = 1 ) 。または。n < 1 ) Stop 91 Forall ( k = 1 : n ) ip ( k ) = k h = 0._SP g = 0._SP isig = 1 detnth = 0._SP pivmax = Maxval ( Abs ( a ( 1 、:))) pivmin = pivmax ! Do k = 1 、n ! Banachiewicz (1938) 式 (2.3) h ( k 、ip ( k :)) = 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 ))) pivmax = Max ( piv , pivmax ) ! 条件を調整しますpivmin = Min ( piv , pivmin ) 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 *, '|AQ-LU| ',Maxval (Abs(a(:,ip(:))-Matmul(g, h(:,ip(:))))) End Subroutine luban Subroutine lusolve ( l , u , ip , x ) ! 三角因子を使用して Ax=a システムを解きます。 LU=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 (:, :) ! X(n,m) 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             

Cコードの例

/* 入力: A - N次元の正方行列の行へのポインタ配列* Tol - 行列が縮退に近い場合にエラーを検出するための小さな許容値* 出力: 行列Aは変更され、LEとUの両方の行列のコピーがA=(LE)+Uとして格納され、P*A=L*Uとなる。* 順列行列は行列としてではなく、順列行列が1である列のインデックスを含む、サイズN+1の整数ベクトルPに格納される。最後の要素P[N]=S+N、* ここでSは行列式の計算に必要な行交換の数、det(P)=(-1)^S */ int LUPDecompose ( double ** A , int N , double Tol , int * P ) {          int i 、j 、k 、imax ; double maxA 、* ptr 、absA ;          for ( i = 0 ; i <= N ; i ++ ) P [ i ] = i ; //単位順列行列、P[N]はNで初期化される            ( i = 0 ; i < N ; i ++ )の場合{ maxA = 0.0 ; imax = i ;               k = i ; k < N ; k ++ )の場合、(( absA = fabs ( A [ k ][ i ] )) > maxA ) { maxA = absA ; imax = k ; }                       if ( maxA < Tol ) return 0 ; //失敗、行列は退化している       if ( imax != i ) { //ピボット P j = P [ i ]; P [ i ] = P [ imax ]; P [ imax ] = j ;               //A の行をピボットptr = A [ i ]; A [ i ] = A [ imax ]; A [ imax ] = ptr ;          //Nから始まるピボットを数える（行列式用）P [ N ] ++ ; }   ( j = i + 1 ; j < N ; j ++ )の場合{ A [ j ][ i ] /= A [ i ][ i ] ;              ( k = i + 1 ; k < N ; k ++ )の場合、 A [ j ][ k ] -= A [ j ][ i ] * A [ i ][ k ] となります。} }                 return 1 ; //分解完了}  /* 入力: LUPDecompose に代入された A,P; b - 右辺ベクトル; N - 次元* 出力: x - A*x=b の解ベクトル*/ void LUPSolve ( double ** A , int * P , double * b , int N , double * x ) {            ( int i = 0 ; i < N ; i ++ )の場合{ x [ i ] = b [ P [ i ]];             ( int k = 0 ; k < i ; k ++ ) x [ i ] -= A [ i ][ k ] * x [ k ] ; }               for ( int i = N - 1 ; i >= 0 ; i -- ) { for ( int k = i + 1 ; k < N ; k ++ ) x [ i ] -= A [ i ][ k ] * x [ k ] ;                            x [ i ] /= A [ i ][ i ]; } }   /* 入力: A、P は LUPDecompose で埋められます。N 次元* 出力: IA は初期行列の逆行列です*/ void LUPInvert ( double ** A 、int * P 、int N 、double ** IA ) { for ( int j = 0 ; j < N ; j ++ ) { for ( int i = 0 ; i < N ; i ++ ) { IA [ i ][ j ] = P [ i ] == j ? 1.0 : 0.0 ;                                        for ( int k = 0 ; k < i ; k ++ ) IA [ i ][ j ] -= A [ i ][ k ] * IA [ k ][ j ]; }               for ( int i = N - 1 ; i >= 0 ; i -- ) { for ( int k = i + 1 ; k < N ; k ++ ) IA [ i ][ j ] -= A [ i ][ k ] * IA [ k ][ j ];                            IA [ i ][ j ] /= A [ i ][ i ]; } } }    /* 入力: A,P は LUPDecompose に入力されます。N 次元。* 出力: 関数は初期行列の行列式を返します*/ double LUPDeterminant ( double ** A , int * P , int N ) {        ダブルdet = A [ 0 ][ 0 ];    ( int i = 1 ; i < N ; i ++ )でdet *= A [ i ][ i ] を実行します。            return ( P [ N ] - N ) % 2 == 0 ?デット: -デット; }           

C#コード例

public class SystemOfLinearEquations { public double [] SolveUsingLU ( double [,] matrix , double [] rightPart , int n ) { // 行列の分解double [,] lu = new double [ n , n ]; double sum = 0 ; for ( int i = 0 ; i < n ; i ++ ) { for ( int j = i ; j < n ; j ++ ) { sum = 0 ; for ( int k = 0 ; k < i ; k ++ ) sum += lu [ i , k ] * lu [ k , j ]; lu [ i , j ] = matrix [ i , j ] - sum ; } for ( int j = i + 1 ; j < n ; j ++ ) { sum = 0 ; for ( int k = 0 ; k < i ; k ++ ) sum += lu [ j , k ] * lu [ k , i ]; lu [ j , i ] = ( 1 / lu [ i , i ]) * (行列[ j , i ] - sum ); } }                                                                                                                   // lu = L+UI // Ly = b の解を求めるdouble [] y = new double [ n ]; for ( int i = 0 ; i < n ; i ++ ) { sum = 0 ; for ( int k = 0 ; k < i ; k ++ ) sum += lu [ i , k ] * y [ k ]; y [ i ] = rightPart [ i ] - sum ; } // Ux = y の解を求めるdouble [] x = new double [ n ]; for ( int i = n - 1 ; i >= 0 ; i -- ) { sum = 0 ; for ( int k = i + 1 ; k < n ; k ++ ) sum += lu [ i , k ] * x [ k ]; x [ i ] = ( 1 / lu [ i , i ]) * ( y [ i ] - sum ); } return x ; } }                                                                                            

MATLABコード例

関数 LU = LUDecompDoolittle ( A ) n =長さ( A ); LU = A ; k = 2 の場合: n i = 1 の場合: k - 1ラムダ= LU ( k , i ) / LU ( i , i ); LU ( k , i ) =ラムダ; LU ( k , i + 1 : n ) = LU ( k , i + 1 : n ) - LU ( i , i + 1 : n ) *ラムダ;終わり、終わり、終わり                                     function  x = SolveLinearSystem ( LU, B ) n = length ( LU ); y = zeros ( size ( B )); % Ly = B の解を求める。for i = 1 : n y ( i ,:) = B ( i ,:) - LU ( i , 1 : i ) * y ( 1 : i ,:); end % Ux = y の解を求めるx = zeros ( size ( B )); for i = n :( - 1 ): 1 x ( i ,:) = ( y ( i ,:) - LU ( i ,( i + 1 ): n ) * x (( i + 1 ): n ,:)) / LU ( i , i ); end end                                       A = [ 4 3 3 ; 6 3 3 ; 3 4 3 ] LU = LUDecompDoolittle ( A ) B = [ 1 2 3 ; 4 5 6 ; 7 8 9 ; 10 11 12 ] ' x = SolveLinearSystem ( LU , B ) A * x                                  

参照

• ブロックLU分解
• ブルハット分解
• コレスキー分解
• クラウト行列分解
• 不完全LU分解
• LU削減
• 行列分解
• QR分解

注記

1. バナチエヴィッチ (1938).
2. Cormen et al. (2009)、p. 819、28.1: 線形方程式系の解法。
3. Okunev & Johnson (1997)、系3。
4. Trefethen & Bau (1997)、p. 166.
5. Trefethen & Bau (1997)、p. 161.
6. Banachiewicz (1938); Lay, Lay & McDonald (2021)、p. 133、2.5: 行列因数分解。
7. Rigotti (2001)、首席副主任。
8. Horn & Johnson (1985)、系3.5.5
9. Horn & Johnson (1985)、定理3.5.2。
10. Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). 「LU分解の可能性の検討」. North American GeoGebra Journal . 9 (1).
11. オクネフ＆ジョンソン（1997年）。
12. ハウスホルダー（1975年）。
13. ゴラブとヴァン・ローン (1996)、112、119 ページ。
14. シャバット、ギル;シュムエリ、ヤニフ。アイゼンバド、ヤリブ。アーフェルブッフ、アミール (2016)。 「ランダム化された LU 分解」。応用および計算高調波解析。44 (2 ) : 246–272。arXiv : 1310.7202 。土井：10.1016/j.acha.2016.04.006。S2CID  1900701。
15. バンチ＆ホップクロフト（1974年）。
16. Trefethen & Bau (1997)、p. 152.
17. ゴラブとヴァン・ローン (1996)、p. 121.
18. ラルストン（1965年）。
19. ハート（2011年）。
20. ドワイヤー（1951年）。

参考文献

• Banachiewicz、T. (1938)、「解決策の数値計算方法…」(PDF)、Bull.インターン。デ・ラカド。ポロネーズ、セリエ A. Sc.数学。 : 393～ 404。
• バンチ, ジェームズ・R.;ホップクロフト, ジョン(1974)、「高速行列乗算による三角因数分解と逆行列計算」、Mathematics of Computation、28 (125): 231– 236、doi : 10.2307/2005828、hdl : 1813/ 6003、ISSN  0025-5718、JSTOR  2005828。
• トーマス・H・コーメン;チャールズ・E・ライザーソン;ロナルド・L・リベスト; Stein、Clifford (2009)、Introduction to Algorithms (第 3 版)、MIT Press および McGraw-Hill、ISBN 978-0-262-03293-3。
• ドワイヤー、ポール・S.（1951）『線形計算』ニューヨーク：ワイリー。
• ゴルブ、ジーン・H. ;ヴァン・ローン、チャールズ・F. (1996) 『行列計算』（第3版）、ボルチモア：ジョンズ・ホプキンス、ISBN 978-0-8018-5414-9。
• ハート、ロジャー（2011）『線形代数の中国的ルーツ』ボルチモア：ジョンズ・ホプキンス、ISBN 978-0801897559。
• ホーン、ロジャー A.; ジョンソン、チャールズ R. (1985)、『マトリックス分析』、ケンブリッジ大学出版局、ISBN 978-0-521-38632-63.5節を参照。N −  1
• ハウスホルダー、アルストン S. (1975)、「数値解析における行列理論」、ニューヨーク：ドーバー出版、MR  0378371。
• レイ、デイビッド・C.; レイ、スティーブン・R.; マクドナルド、ジュディ・J. (2021) 『線形代数とその応用（第6版）』ピアソン、ISBN 978-0-13-585125-8。
• Okunev, Pavel; Johnson, Charles R. (1997),任意行列のLU分解の存在のための必要十分条件, arXiv : math.NA/0506382。
• プール、デイビッド（2006年）、線形代数：現代入門（第2版）、カナダ：トムソンブルックス/コール、ISBN 978-0-534-99845-5。
• ラルストン、アンソニー（1965年）、数値解析入門、ニューヨーク：マグロウヒル社、ISBN 978-0-070-51157-6。
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007)「セクション2.3」、Numerical Recipes: The Art of Scientific Computing (第3版)、ニューヨーク: Cambridge University Press、ISBN 978-0-521-88068-8。
• トレフェテン、ロイド・N.；バウ、デイビッド（1997）、数値線形代数、フィラデルフィア：産業応用数学協会、ISBN 978-0-89871-361-9。
• リゴッティ、ルカ（2001）、ECON 2001 - 数学的手法入門、講義8

外部リンク

参考文献

• MathWorldでの LU 分解。
• Math-Linuxでの LU 分解。
• Holistic Numerical Methods InstituteにおけるLU分解
• LU行列分解。MATLABリファレンス。

コンピュータコード

• LAPACKは、密な線形代数問題を解くためのFORTRANサブルーチンのコレクションです。
• ALGLIB には、LAPACK の C++、C#、Delphi などへの部分的な移植が含まれています。
• C++ コード、デイトン大学のJ. Loomis 教授
• Cコード、数学ソースライブラリ
• Rustコード
• X10のLU

オンラインリソース

• LU分解を用いた線形方程式のシステム記述的解法を提供するWebアプリ
• LU分解を含む手順付き行列計算機
• LU分解ツール、uni-bonn.de
• LU 分解、Ed Pegg, Jr.著、The Wolfram Demonstrations Project、2007 年。

Retrieved from "https://en.wikipedia.org/w/index.php?title=LU_decomposition&oldid=1320197026"