LU decomposition

Original article: w:LU decomposition

수치 해석(numerical analysis)과 선형 대수(linear algebra)에서, 아래쪽–위쪽(lower–upper, LU) 분해 또는 인수분해는 행렬(matrix)을 아래쪽 삼각 행렬(triangular matrix)과 위쪽 삼각 행렬의 곱으로 인수분해합니다 (행렬 분해를 참조). 그 곱은 때때로 순열 행렬(permutation matrix)도 포함됩니다. LU 분해는 가우스 소거법(Gaussian elimination)의 행렬 형식으로 볼 수 있습니다. 컴퓨터는 보통 LU 분해를 사용하여 선형 방정식의 제곱 시스템을 해결하고, 그것은 행렬을 반전시키거나 행렬의 행렬식(determinant)을 계산할 때 중요한 단계이기도 합니다. LU 분해는 1938년 폴란드 천문학자 타데우시 바나히어비츠(Tadeusz Banachiewicz)에 의해 도입되었습니다. 인용하자면: "가우스와 둘리틀은 [소거의] 방법을 대칭 방정식에만 적용한 것으로 보입니다. 더 최근의 저자, 예를 들어, Aitken, Banachiewicz, Dwyer, 및 Crout는 ... 비-대칭 문제와 관련하여 이 방법이나 그 변형의 사용을 강조했습니다 ... Banachiewicz는 ... 기본 문제가 실제로 행렬 인수분해, 또는 그가 부르는 "분해" 중 하나라는 ... 요점을 보았습니다." 그것은 역시 LR 분해 (왼쪽과 오른쪽 삼각 행렬로의 인수화)라고 참조됩니다.

Definitions

LDU decomposition of a Walsh matrix

$A$를 정사각 행렬이라고 놓습니다. LU 분해는 적절한 행 및/또는 열 순서화 또는 순열을 갖는 $A$를 두 가지 인수, 즉 아래쪽 삼각 행렬 $L$과 위쪽 삼각 행렬 $U$로 분해하는 것을 참조합니다:

$A=LU.$

아래쪽 삼각 행렬에서는 대각선 위의 모든 원소가 영이고, 위쪽 삼각 행렬에서는 대각선 아래의 모든 원소가 영입니다. 예를 들어 3 × 3 행렬 $A$에 대해, LU 분해는 다음과 같습니다:

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

행렬에 적절한 순서화나 순열 없이, 인수분해가 실현되지 않을 수 있습니다. 예를 들어, ${\textstyle a_{11}=\ell _{11}u_{11}}$임을 (행렬 곱셈을 확장함으로써) 쉽게 확인할 수 있습니다. 만약 ${\textstyle a_{11}=0}$이면, ${\textstyle \ell _{11}}$과 ${\textstyle u_{11}}$ 중 적어도 하나는 영이어야 하며, 이는 $L$ 또는 $U$가 특이(singular)함을 의미합니다. 이것은 $A$가 비-특이 (역-가능)이면 불가능합니다. 이것은 절차상의 문제입니다. 그것은 순열된 행렬의 첫 번째 원소가 비영이 되도록 $A$의 행을 간단히 다시-순서화함으로써 제거될 수 있습니다. 후속 인수분해 단계에서 같은 문제는 동일한 방법으로 제거될 수 있습니다; 아래의 기본 절차를 참조하십시오.

LU factorization with partial pivoting

행 (또는 열)의 적절한 순열은 LU 분해에 대해 충분 조건인 것으로 밝혀졌습니다. 부분 피봇팅을 갖는 LU 인수분해 (LUP)는 종종 행 순열만 갖는 LU 분해를 참조합니다:

$PA=LU,$

여기서 $L$과 $U$는 다시 아래쪽 삼각 행렬과 위쪽 삼각 행렬이고, $P$는 $A$에 왼쪽-곱셈을 하면 $A$의 행을 다시-순서화하는 순열 행렬(permutation matrix)입니다. 모든 정사각 행렬은 이 형식으로 인수화될 수 있음이 밝혀졌고, 그 인수분해는 실제로 수치적으로 안정적입니다. 이것은 LUP 분해를 실제로 유용한 기술로 만듭니다.

LU factorization with full pivoting

전체 피봇팅을 갖는 LU 인수분해(LU factorization with full pivoting)는 행과 열 순열 둘 다를 포함합니다:

$PAQ=LU,$

여기서 $L$, $U$, 및 $P$는 이전처럼 정의되고, $Q$는 $A$의 열을 다시 순서화하는 순열 행렬입니다.

Lower-diagonal-upper (LDU) decomposition

아래쪽-대각-위쪽(Lower-diagonal-upper, LDU) 분해는 다음 형식의 분해입니다:

$A=LDU,$

여기서 $D$는 대각 행렬(diagonal matrix)이고, $L$과 $U$는 단위삼각 행렬(unitriangular matrices)이며, $L$과 $U$의 대각선 위의 모든 엔트리가 1임을 의미합니다.

Rectangular matrices

위에서는 $A$가 정사각 행렬이어야 했지만, 이들 분해는 모두 직사각 행렬로 일반화될 수 있습니다. 해당 경우에서, $L$과 $D$는 둘 다 $A$와 같은 행의 개수를 가지는 정사각 행렬이고, $U$는 $A$와 정확하게 같은 차원을 가집니다. 위쪽 삼각(Upper triangular)은 주요 대각선 아래에 영 엔트리만 가지는 것으로 해석되어야 하며, 이는 위쪽 왼쪽 모서리에서 시작합니다. 유사하게, $U$에 대한 더 정확한 용어는 그것이 행렬 $A$의 "행 사다리꼴 형식"이라는 것입니다.

Example

다음 2x2 행렬을 인수화합니다:

${\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 분해를 찾는 한 가지 방법은 단순히 검사를 통해 선형 방정식을 푸는 것입니다. 행렬 곱셈(matrix multiplication)을 전개하면 다음이 제공됩니다:

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

이 방정식 시스템은 미달-결정된(underdetermined) 것입니다. 이 경우에서, $L$과 $U$ 행렬의 임의의 두 개의 비-영 원소는 해의 매개변수이고 임의의 비-영 값으로 임의적으로 설정될 수 있습니다. 그러므로, 고유한 LU 분해를 찾기 위해, $L$과 $U$ 행렬에 일부 제한을 가할 필요가 있습니다. 예를 들어, 아래쪽 삼각 행렬 $L$이 단위-삼각 행렬이 되도록 편리하게 요구할 수 있습니다 (즉, 주요 대각선의 모든 엔트리를 1로 설정합니다). 그런-다음 방정식의 시스템은 다음 해를 가집니다:

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

이들 값을 위의 LU 분해에 대체하면 다음이 생성됩니다:

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

Existence and uniqueness

Square matrices

임의의 정사각 행렬 ${\textstyle A}$는 LUP 및 PLU 인수분해를 허용합니다. 만약 ${\textstyle A}$가 역가능(invertible)이면, 그것이 LU (또는 LDU) 인수분해를 허용하는 것과 모든 선행하는 주요 소행렬식(minors)이 비-영인 것은 필요충분 조건입니다 (예를 들어 ${\begin{bmatrix}0&1\\1&0\end{bmatrix}}$는 LU 또는 LDU 인수분해를 허용하지 않습니다). 만약 ${\textstyle A}$가 랭크 ${\textstyle k}$의 특이 행렬이면, 그것은 처음 ${\textstyle k}$개의 선행하는 주요 소행렬식이 비-영이면 LU 분해를 허용되지만, 그 전환은 참이 아닙니다.

만약 정사각, 역 행렬이 LDU (${\textstyle L}$과 ${\textstyle U}$의 모든 대각선 엔트리가 1과 같은 값을 갖는 인수분해)를 가지면, 그 인수분해는 고유합니다. 해당 경우에서, LU 인수분해는 ${\textstyle L}$ (또는 ${\textstyle U}$)의 대각선이 1로 구성되어야 함을 요구하면 역시 고유합니다.

일반적으로, 임의의 정사각 행렬 $A_{n\times n}$은 다음 중 하나를 가질 수 있습니다:

1. 고유한 LU 인수분해 (위에서 언급한 것);
2. 임의의 처음 (n−1) 열 중 두 개 이상이 선형적으로 종속이거나 처음 (n−1) 열 중 하나가 0이면 무한하게 많은 LU 인수분해;
3. 처음 (n−1) 열이 비-영이고 선형적으로 독립이고 적어도 하나의 선행하는 주요 소행렬식이 영이면 LU 분해가 없음.

경우 3에서, 선행하는 주요 소행렬식이 0이 되지 않도록 대각선 엔트리 $a_{jj}$를 $a_{jj}\pm \varepsilon$로 변경함으로써 LU 인수분해를 근사화할 수 있습니다.

Symmetric positive-definite matrices

만약 ${\textstyle A}$가 대칭 (또는 ${\textstyle A}$가 복소수이면 에르미트) 양수-한정 행렬(positive-definite matrix)이면, 우리는 U가 L의 켤레 전치(conjugate transpose)가 되도록 문제를 배열할 수 있습니다. 즉, A를 다음으로 쓸 수 있습니다:

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

이 분해는 숄레스키 분해(Cholesky decomposition)라고 불립니다. 숄레스키 분해는 행렬이 양수 한정이라는 조건으로 하여 항상 존재하고 고유합니다. 게다가, 숄레스키 분해를 계산하는 것은 일부 다른 LU 분해를 계산하는 것보다 더 효율적이고 수치적으로 더 안정적(numerically more stable)입니다.

General matrices

임의의 필드에 걸쳐 (반드시 역행렬이 아닌) 행렬에 대해, LU 분해를 가지는 정확한 필요 조건과 충분 조건이 알려져 있습니다. 조건은 특정 부분행렬의 랭크의 관점에서 표현됩니다. LU 분해를 얻기 위한 가우스 제거 알고리즘도 이 가장 일반적인 경우로 확장되었습니다.

Algorithms

Closed formula

LDU 분해가 존재하고 고유할 때, 원래 행렬 ${\textstyle A}$의 특정 부분행렬의 행렬식의 비율 관점에서 L, D, 및 U의 원소에 대한 닫힌 (명시적) 공식이 있습니다. 특히, ${\textstyle D_{1}=A_{1,1}}$, 그리고 ${\textstyle i=2,\ldots ,n}$에 대해, ${\textstyle D_{i}}$는 ${\textstyle i}$-번째 주요 부분행렬과 ${\textstyle (i-1)}$-번째 주요 부분행렬의 비율입니다. 행렬식의 계산은 계산적으로 비용이 많이 들기 때문에, 이 명시적 공식은 실제로 사용되지 않습니다.

Using Gaussian elimination

다음 알고리즘은 본질적으로 가우스 소거법(Gaussian elimination)의 수정된 형식입니다. 이 알고리즘을 사용하여 LU 분해를 계산하는 것은 낮은-차수 항을 무시하여, ${\tfrac {2}{3}}n^{3}$개의 부동-점 연산이 필요합니다. 부분 피봇팅(pivoting)은 2차 항만 더합니다; 이는 전체 피벗팅에 대한 경우가 아닙니다.

Procedure

N × N 행렬 $A=(a_{i,j})_{1\leq i,j\leq N}$가 주어졌을 때, $A^{(0)}$을 첫 번째 열에 대해 (부분 피벗팅과 같은) 원하는 조건을 충족하기 위해 필요한 행이 교체된 행렬 $A$로 정의합니다. 행렬 $A$의 괄호-처진 위첨자 (예를 들어, $(0)$)는 행렬의 버전입니다. 행렬 $A^{(n)}$은 처음 $n$개 열에 대해 가우스 소거법을 통해 주요 대각선(main diagonal) 아래의 원소가 이미 0으로 제거되었고, $(n+1)^{th}$ 열에 대한 원하는 조건을 충족하도록 필요한 행을 교체한 $A$ 행렬입니다.

$A^{(n-1)}$의 n-번째 열에서 주요 대각선 아래에 원소 ($a_{i,n}^{(n-1)}$로 표시됨, 여기서 $i=n+1,\dotsc ,N$)를 갖는 각 행 $i$에 대해 $row_{i}=row_{i}-(\ell _{i,n})\cdot row_{n}$ 연산을 수행합니다. 이 연산을 위해,

$\ell _{i,n}:={\frac {a_{i,n}^{(n-1)}}{a_{n,n}^{(n-1)}}}.$

원소 $a_{i,n}^{(n-1)}$을 영으로 제거하기 위해 이들 행 연산을 수행합니다. 일단 이들 행을 뺀 후에는, 행을 교환하여 $(n+1)^{th}$ 열에 원하는 조건을 제공할 수 있습니다. 부분 피벗팅을 수행하기 위해 또는 주요 대각선의 $a_{n+1,n+1}$ 원소가 영이기 때문에 (그리고 따라서 가우스 소거법을 구현하기 위해 사용될 수 없기 때문에) 여기에서 행을 교환할 수 있습니다.

최종 순열 행렬(permutation matrix) $P$를 $A$ 행렬과 같은 순서로 교환된 모든 같은 행을 가지는 항등 행렬로 정의합니다.

일단 처음 $N-1$ 열에 대한 행 연산을 수행했으면, $U$로 표시되는 위쪽 삼각(upper triangular) 행렬 $A^{(N-1)}$을 얻습니다. 참고로 $A^{(N-1)}$의 N-번째 열에는 행을 교체해야 하는 조건이 없기 때문에 $U$를 $A^{(N-1)}$로 표시할 수 있습니다. 역시 아래 수식을 통해 $\ell _{i,n}$의 값ㄷ르의 값을 직접 입력함으로써 $PA=LU$를 만족하는 $PA=LU$로 표시되는 아래쪽 삼각(lower triangular) 행렬을 계산할 수 있습니다:

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

Example

만약 다음 행렬이 주어지면$$A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},$$부분 피벗팅을 구현하고 따라서 행렬 $A$와 $P$ 행렬의 첫 번째 반복이 각각 다음이 되도록 첫 번째 행과 두 번째 행을 교환합니다:$$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}}.$$일단 행을 교환한 후에는, 다음을 수행함으로써 첫 번째 열의 주요 대각선 아래 원소를 제거할 수 있습니다: $${\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}}$$이때 다음임을 만족합니다,$${\begin{alignedat}{0}\ell _{2,1}={\frac {0}{4}}=0\\\ell _{3,1}={\frac {2}{4}}=0.5\end{alignedat}}$$일단 이들 행을 빼면, $A^{(0)}$에서 다음 행렬이 유도됩니다: $${\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.$$부분 피벗팅을 구현하고 있기 때문에, 유도된 행렬의 두 번째와 세 번째 행과 $P$ 행렬의 현재 버전을 각각 교환하여 다음을 얻습니다:$$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}}.$$이제, ${\textstyle \ell _{3,2}={\frac {5}{6}}}$임을 만족하는 $row_{3}=row_{3}-(\ell _{3,2})\cdot row_{2}$를 수행함으로써 두 번째 열의 주요 대각선 아래 원소를 제거합니다. 이 행 빼기 이후 $A$의 현재 반복에서 주요 대각선 아래에 비-영 원소가 존재하지 않기 때문에, 이 행 빼기는 최종 $A$ 행렬 ($U$로 표시됨)과 최종 $P$ 행렬을 유도합니다:$$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}}.$$이제 최종 $L$ 행렬을 얻을 수 있습니다:$$L={\begin{pmatrix}1&0&0\\\ell _{2,1}&1&0\\\ell _{3,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$임을 만족하는 관계를 가집니다.

Relations when no rows are swapped

If we did not swap rows at all during this process, we can perform the row operations simultaneously for each column $n$ by setting $A^{(n)}:=L_{n}A^{(n-1)},$ where $L_{n}$ is the N × N identity matrix with its n-th column replaced by the transposed vector ${\begin{pmatrix}0&\dotsm &0&1&-\ell _{n+1,n}&\dotsm &-\ell _{N,n}\end{pmatrix}}^{\textsf {T}}.$ In other words, the lower triangular matrix

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

Performing all the row operations for the first $N-1$ columns using the $A^{(n)}:=L_{n}A^{(n-1)}$ formula is equivalent to finding the decomposition $$A=L_{1}^{-1}L_{1}A^{(0)}=L_{1}^{-1}A^{(1)}=L_{1}^{-1}L_{2}^{-1}L_{2}A^{(1)}=L_{1}^{-1}L_{2}^{-1}A^{(2)}=\dotsm =L_{1}^{-1}\dotsm L_{N-1}^{-1}A^{(N-1)}.$$ Denote ${\textstyle L=L_{1}^{-1}\dotsm L_{N-1}^{-1}}$ so that $A=LA^{(N-1)}=LU$.

Now let's compute the sequence of $L_{1}^{-1}\dotsm L_{N-1}^{-1}$. We know that $L_{i}^{-1}$ has the following formula.

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

If there are two lower triangular matrices with 1s in the main diagonal, and neither have a non-zero item below the main diagonal in the same column as the other, then we can include all non-zero items at their same location in the product of the two matrices. For example:

$\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)$

Finally, multiply $L_{i}^{-1}$ together and generate the fused matrix denoted as ${\textstyle L}$ (as previously mentioned). Using the matrix ${\textstyle L}$, we obtain $A=LU.$

It is clear that in order for this algorithm to work, one needs to have $a_{n,n}^{(n-1)}\neq 0$ at each step (see the definition of $\ell _{i,n}$). If this assumption fails at some point, one needs to interchange n-th row with another row below it before continuing. This is why an LU decomposition in general looks like $P^{-1}A=LU$.

LU Crout decomposition

Note that the decomposition obtained through this procedure is a Doolittle decomposition: the main diagonal of L is composed solely of 1s. If one would proceed by removing elements above the main diagonal by adding multiples of the columns (instead of removing elements below the diagonal by adding multiples of the rows), we would obtain a Crout decomposition, where the main diagonal of U is of 1s.

Another (equivalent) way of producing a Crout decomposition of a given matrix A is to obtain a Doolittle decomposition of the transpose of A. Indeed, if ${\textstyle A^{\textsf {T}}=L_{0}U_{0}}$ is the LU-decomposition obtained through the algorithm presented in this section, then by taking ${\textstyle L=U_{0}^{\textsf {T}}}$ and ${\textstyle U=L_{0}^{\textsf {T}}}$, we have that $A=LU$ is a Crout decomposition.

Through recursion

Cormen et al. describe a recursive algorithm for LUP decomposition.

Given a matrix A, let P₁ be a permutation matrix such that

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

where ${\textstyle a\neq 0}$, if there is a nonzero entry in the first column of A; or take P₁ as the identity matrix otherwise. Now let ${\textstyle c=1/a}$, if ${\textstyle a\neq 0}$; or ${\textstyle c=0}$ otherwise. We have

$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).$

Now we can recursively find an LUP decomposition ${\textstyle P'\left(A'-cvw^{\textsf {T}}\right)=L'U'}$. Let ${\textstyle v'=P'v}$. Therefore

$\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),$

which is an LUP decomposition of A.

Randomized algorithm

It is possible to find a low rank approximation to an LU decomposition using a randomized algorithm. Given an input matrix ${\textstyle A}$ and a desired low rank ${\textstyle k}$, the randomized LU returns permutation matrices ${\textstyle P,Q}$ and lower/upper trapezoidal matrices ${\textstyle L,U}$ of size ${\textstyle m\times k}$ and ${\textstyle k\times n}$ respectively, such that with high probability ${\textstyle \left\|PAQ-LU\right\|_{2}\leq C\sigma _{k+1}}$, where ${\textstyle C}$ is a constant that depends on the parameters of the algorithm and ${\textstyle \sigma _{k+1}}$ is the ${\textstyle (k+1)}$-th singular value of the input matrix ${\textstyle A}$.

Theoretical complexity

If two matrices of order n can be multiplied in time M(n), where M(n) ≥ n^a for some a > 2, then an LU decomposition can be computed in time O(M(n)). This means, for example, that an O(n^2.376) algorithm exists based on the Coppersmith–Winograd algorithm.

Sparse-matrix decomposition

Special algorithms have been developed for factorizing large sparse matrices. These algorithms attempt to find sparse factors L and U. Ideally, the cost of computation is determined by the number of nonzero entries, rather than by the size of the matrix.

These algorithms use the freedom to exchange rows and columns to minimize fill-in (entries that change from an initial zero to a non-zero value during the execution of an algorithm).

General treatment of orderings that minimize fill-in can be addressed using graph theory.

Applications

Solving linear equations

Given a system of linear equations in matrix form

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

we want to solve the equation for x, given A and b. Suppose we have already obtained the LUP decomposition of A such that ${\textstyle PA=LU}$, so ${\textstyle LU\mathbf {x} =P\mathbf {b} }$.

In this case the solution is done in two logical steps:

1. First, we solve the equation ${\textstyle L\mathbf {y} =P\mathbf {b} }$ for y.
2. Second, we solve the equation ${\textstyle U\mathbf {x} =\mathbf {y} }$ for x.

In both cases we are dealing with triangular matrices (L and U), which can be solved directly by forward and backward substitution without using the Gaussian elimination process (however we do need this process or equivalent to compute the LU decomposition itself).

The above procedure can be repeatedly applied to solve the equation multiple times for different b. In this case it is faster (and more convenient) to do an LU decomposition of the matrix A once and then solve the triangular matrices for the different b, rather than using Gaussian elimination each time. The matrices L and U could be thought to have "encoded" the Gaussian elimination process.

The cost of solving a system of linear equations is approximately ${\textstyle {\frac {2}{3}}n^{3}}$ floating-point operations if the matrix ${\textstyle A}$ has size ${\textstyle n}$. This makes it twice as fast as algorithms based on QR decomposition, which costs about ${\textstyle {\frac {4}{3}}n^{3}}$ floating-point operations when Householder reflections are used. For this reason, LU decomposition is usually preferred.

Inverting a matrix

When solving systems of equations, b is usually treated as a vector with a length equal to the height of matrix A. In matrix inversion however, instead of vector b, we have matrix B, where B is an n-by-p matrix, so that we are trying to find a matrix X (also a n-by-p matrix):

$AX=LUX=B.$

We can use the same algorithm presented earlier to solve for each column of matrix X. Now suppose that B is the identity matrix of size n. It would follow that the result X must be the inverse of A.

Computing the determinant

Given the LUP decomposition ${\textstyle A=P^{-1}LU}$ of a square matrix A, the determinant of A can be computed straightforwardly as

$\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).$

The second equation follows from the fact that the determinant of a triangular matrix is simply the product of its diagonal entries, and that the determinant of a permutation matrix is equal to (−1)^S where S is the number of row exchanges in the decomposition.

In the case of LU decomposition with full pivoting, ${\textstyle \det(A)}$ also equals the right-hand side of the above equation, if we let S be the total number of row and column exchanges.

The same method readily applies to LU decomposition by setting P equal to the identity matrix.

Code examples

C code example

/* INPUT: A - array of pointers to rows of a square matrix having dimension N
 *        Tol - small tolerance number to detect failure when the matrix is near degenerate
 * OUTPUT: Matrix A is changed, it contains a copy of both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
 *        The permutation matrix is not stored as a matrix, but in an integer vector P of size N+1 
 *        containing column indexes where the permutation matrix has "1". The last element P[N]=S+N, 
 *        where S is the number of row exchanges needed for determinant computation, 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; //Unit permutation matrix, P[N] initialized with N

    for (i = 0; i < N; i++) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k++)
            if ((absA = fabs(A[k][i])) > maxA) { 
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
            ptr = A[i];
            A[i] = A[imax];
            A[imax] = ptr;

            //counting pivots starting from N (for determinant)
            P[N]++;
        }

        for (j = i + 1; j < N; j++) {
            A[j][i] /= A[i][i];

            for (k = i + 1; k < N; k++)
                A[j][k] -= A[j][i] * A[i][k];
        }
    }

    return 1;  //decomposition done 
}

/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double **A, int *P, double *b, int N, double *x) {

    for (int i = 0; i < N; i++) {
        x[i] = b[P[i]];

        for (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];
    }
}

/* INPUT: A,P filled in LUPDecompose; N - dimension
 * OUTPUT: IA is the inverse of the initial matrix
 */
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];
        }
    }
}

/* INPUT: A,P filled in LUPDecompose; N - dimension. 
 * OUTPUT: Function returns the determinant of the initial matrix
 */
double LUPDeterminant(double **A, int *P, int N) {

    double det = A[0][0];

    for (int i = 1; i < N; i++)
        det *= A[i][i];

    return (P[N] - N) % 2 == 0 ? det : -det;
}

C# code example

public class SystemOfLinearEquations
{
    public double[] SolveUsingLU(double[,] matrix, double[] rightPart, int n)
    {
        // decomposition of matrix
        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]) * (matrix[j, i] - sum);
            }
        }

        // lu = L+U-I
        // find solution of 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;
        }
        // find solution of 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 code example

function LU = LUDecompDoolittle(A)
    n = length(A);
    LU = A;
    % decomposition of matrix, Doolittle's Method
    for i = 1:1:n
        for j = 1:(i - 1)
            LU(i,j) = (LU(i,j) - LU(i,1:(j - 1))*LU(1:(j - 1),j)) / LU(j,j);
        end
        j = i:n;
        LU(i,j) = LU(i,j) - LU(i,1:(i - 1))*LU(1:(i - 1),j);
    end
    %LU = L+U-I
end

function x = SolveLinearSystem(LU, B)
    n = length(LU);
    y = zeros(size(B));
    % find solution of Ly = B
    for i = 1:n
        y(i,:) = B(i,:) - LU(i,1:i)*y(1:i,:);
    end
    % find solution of 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

See also

• Block LU decomposition
• Bruhat decomposition
• Cholesky decomposition
• Crout matrix decomposition
• Incomplete LU factorization
• LU Reduction
• Matrix decomposition
• QR decomposition

References

• Bunch, James R.; Hopcroft, John (1974), "Triangular factorization and inversion by fast matrix multiplication", Mathematics of Computation, 28 (125): 231–236, doi:10.2307/2005828, ISSN 0025-5718, JSTOR 2005828.
• Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford (2001), Introduction to Algorithms, MIT Press and McGraw-Hill, ISBN 978-0-262-03293-3.
• Golub, Gene H.; Van Loan, Charles F. (1996), Matrix Computations (3rd ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
• Horn, Roger A.; Johnson, Charles R. (1985), Matrix Analysis, Cambridge University Press, ISBN 978-0-521-38632-6. See Section 3.5. N − 1
• Householder, Alston S. (1975), The Theory of Matrices in Numerical Analysis, New York: Dover Publications, MR 0378371.
• Okunev, Pavel; Johnson, Charles R. (1997), Necessary And Sufficient Conditions For Existence of the LU Factorization of an Arbitrary Matrix, arXiv:math.NA/0506382.
• Poole, David (2006), Linear Algebra: A Modern Introduction (2nd ed.), Canada: Thomson Brooks/Cole, ISBN 978-0-534-99845-5.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Section 2.3", Numerical Recipes: The Art of Scientific Computing (3rd ed.), New York: Cambridge University Press, ISBN 978-0-521-88068-8.
• Trefethen, Lloyd N.; Bau, David (1997), Numerical linear algebra, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.
• Rigotti, Luca (2001), ECON 2001 - Introduction to Mathematical Methods, Lecture 8

External links

References

• LU decomposition on MathWorld.
• LU decomposition on Math-Linux.
• LU decomposition at Holistic Numerical Methods Institute
• LU matrix factorization. MATLAB reference.

Computer code

• LAPACK is a collection of FORTRAN subroutines for solving dense linear algebra problems
• ALGLIB includes a partial port of the LAPACK to C++, C#, Delphi, etc.
• C++ code, Prof. J. Loomis, University of Dayton
• C code, Mathematics Source Library
• Rust code
• LU in X10

Online resources

• WebApp descriptively solving systems of linear equations with LU Decomposition
• Matrix Calculator with steps, including LU decomposition,
• LU Decomposition Tool, uni-bonn.de
• LU Decomposition by Ed Pegg, Jr., The Wolfram Demonstrations Project, 2007.