Movatterモバイル変換

LU decomposition

From Wikipedia, the free encyclopedia

Type of matrix factorization

Innumerical analysis andlinear algebra,lower–upper (LU)decomposition orfactorization factors amatrix as the product of a lowertriangular matrix and an upper triangular matrix (seematrix multiplication andmatrix decomposition). The product sometimes includes apermutation matrix as well. LU decomposition can be viewed as the matrix form ofGaussian elimination. Computers usually solve squaresystems of linear equations using LU decomposition, and it is also a key step when inverting a matrix or computing thedeterminant of a matrix. It is also sometimes referred to asLR decomposition (factors into left and right triangular matrices). The LU decomposition was introduced by the Polish astronomerTadeusz Banachiewicz in 1938,^[1] who first wrote product equation $LU=A=h^{T}g$ (The last form in his alternate yet equivalent matrix notation appears as $g\times h$ ). While most ellimination and LU algorithms employ equations of the form $element:=element-vector\times vector$ Banachiewicz wrote it as $vector:=vector-submatrix\times vector$ , in this way revealing additional problem symmetry.

Definitions

[edit]

LetA be a square matrix. AnLU factorization refers to expression ofA into product of two factors – a lower triangular matrixL and an upper triangular matrixU:A =LU. Sometimes factorization is impossible without prior reordering ofA to prevent division by zero or uncontrolled growth of rounding errors hence alternative expression becomesPAQ =LU, where in formal notationpermutation matrices factorsP andQ indicate permutation of rows (or columns) ofA. In theoryP (orQ) are obtained by permutations of rows (or columns) of theidentity matrix, in practice the corresponding permutations are applied directly to rows (or columns) ofA.

MatrixA of siden has $n^{2}$ coefficients while two triangle matrices combined containn(n + 1) coefficients, thereforen coefficients of matricesLU are not independent. Usual convention is to setL unitriangular, i.e. with alln main diagonal elements equal one. However, setting insteadU matrix unitriangular reduces to the same procedure aftertranspose of matrix product (cf. properties of matrix transposition): $B=A^{\textsf {T}}=(LU)^{\textsf {T}}=U^{\textsf {T}}L^{\textsf {T}}.$ After transposition,U^T is lower triangle whileL^T is upper unitriangular factor ofB. This demonstrates also, that operations on rows (e.g. pivoting) are equivalent to those on columns of a transposed matrix, and in general choice of row or column algorithm offers no advantage.

In the lower triangular matrix all elements above the main diagonal are zero, in the upper triangular matrix, all the elements below the diagonal are zero. For example, for a3 × 3 matrixA, its LU decomposition looks like this: ${\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}}.$

Without a proper ordering or permutations in the matrix, the factorization may fail to materialize. For example, it is easy to verify (by expanding thematrix multiplication) that ${\textstyle a_{11}=\ell _{11}u_{11}}$ . If ${\textstyle a_{11}=0}$ , then at least one of ${\textstyle \ell _{11}}$ and ${\textstyle u_{11}}$ has to be zero, which implies that eitherL orU issingular. This is impossible ifA is nonsingular (invertible). In terms of operations, zeroing/elimination of remaining elements of first column ofA involves division of ${\textstyle a_{21},a_{31}}$ with ${\textstyle a_{11}}$ , impossible if it is 0. This is a procedural problem. It can be removed by simply reordering the rows ofA so that the first element of the permuted matrix is nonzero. The same problem in subsequent factorization steps can be removed the same way. For numerical stability against rounding errors/division by small numbers it is important to select ${\textstyle a_{11}}$ of large absolute value (cf. pivoting).

LU Through recursion

[edit]

The above example of3 × 3 matrices demonstrates that matrix product of top row and leftmost columns of involved matrices plays special role forLU to succeed. Let us mark consecutive versions of matrices with $(0),\;(1),\dots$ and then let us write matrix product $A\equiv A^{(0)}=L^{(0)}U^{(0)}$ in such way that these rows and columns are separated from the rest. In doing so we shall useblock matrix notation, such that e.g. $a\equiv a_{11}$ is an ordinary number, ${\bf {w}}^{\textsf {T}}\equiv (a_{12},a_{13})^{\textsf {T}}$ is a row vector and ${\bf {v}}=(a_{21},a_{31})$ is a column vector and $A^{'} {\displaystyle A'}$ is sub-matrix of matrix $A^{(0)}$ without top row and leftmost column. Then we can replace $A^{(0)}=L^{(0)}U^{(0)}$ with a blockmatrix product. Namely it turns out that one can multiply matrix blocks in such way as if they were ordinary numbers, i.e. row times column, except that now their components are sub-matrices, sometimes reduced to scalars or vectors. Thus $u{\bf {l}}$ denotes a vector obtained from ${\bf {l}}$ after multiplication of each component by a number $u {\displaystyle u}$ , ${\bf {lu}}^{\textsf {T}}$ is anouter product of vectors ${\bf {l,u}}$ , i.e. a matrix which first column is $u_{12}{\bf {l}}$ , next is $u_{13}{\bf {l}}$ and so on for all components of ${\bf {u}}$ and $L^{(1)}U^{(1)}$ is a product of sub-matrices of $L^{(0)},\;U^{(0)}$ ${\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}}$

From equality of first and last matrices follow final $u=a$ , ${\bf {u=w}}$ , ${\bf {l}}={(1/a)}{\bf {v}}$ while matrix $A^{'} {\displaystyle A'}$ becomes updated/replaced with $A^{(1)}\equiv L^{(1)}U^{(1)}=$ $A'-{\bf {lu}}^{\textsf {T}}$ . Now comes the crucial observation: nothing prevents us to treat $A^{(1)}$ the same way as we did with $A^{(0)}$ , repeatedly. If dimension of $A {\displaystyle A}$ isn × n, aftern− 1 such steps all columns ${\bf {v}}$ form sub-diagonal part of triangle matrix $L {\displaystyle L}$ and all pivots $a {\displaystyle a}$ combined with rows ${\bf {w}}^{\textsf {T}}$ form upper triangle matrix $U {\displaystyle U}$ , as required. In the above examplen = 3 so only two steps suffice.

The above procedure demonstrates that at no step the top diagonal pivot element $a {\displaystyle a}$ of consecutive sub-matrices can be zero. To avoid it columns or rows may be swapped so that $a {\displaystyle a}$ becomes nonzero.Such procedure involving permutation is calledLUP, decomposition with pivoting.

Permutation of columns corresponds to matrix product $AQ^{(0)}$ where $Q^{(0)}$ is a permutation matrix, i.e. the identity matrix $I {\displaystyle I}$ after the same column permutation. After all steps such LUP decomposition applies to $AQ^{(0)}\cdots Q^{(n-1)}\equiv AQ=LU$ . Present computation scheme and similar in Cormen et al.^[2] are examples ofrecurrence algorithms. They demonstrate two general properties of LU factorization:

the need for pivoting at each step; and
final values ofL andU matrices are obtained gradually, one row or column per step.

Recurrence algorithms are not overly costly in terms of algebraic operations yet they suffer from practical disadvantage due to need to update and store most elements ofA at each step. It will be seen that by reordering calculations it is possible to dispose with storage of intermediate values.

LU factorization with partial pivoting

[edit]

It turns out that a proper permutation of rows (or columns) to select column (or row) absolute maximal pivota₁₁ is sufficient for numerically stable LU factorization, except for known pathological cases. It is calledLU factorization with partial pivoting (LUP): $PA=LU,\quad (AQ=LU),$ whereL andU are again lower and upper triangular matrices, andP andQ are correspondingpermutation matrices, which, when correspondingly left- and right-multiplied toA, reorder the rows and columns ofA. It turns out that all square matrices can be factorized in this form,^[3] and the factorization is numerically stable in practice.^[4] This makes LUP decomposition a useful technique in practice.

A variant calledrook pivoting at each step involves search of maximum element the way rook moves on a chessboard, along column, row, column again and so on till reaching a pivot maximal in both its row and column. It can be proven that for large matrices of random elements its cost of operations at each step is similarly to partial pivoting proportional to the length of matrix side unlike its square for full pivoting.

LU factorization with full pivoting

[edit]

AnLU factorization with full pivoting involves both row and column permutations to find absolute maximum element in the whole submatrix: $PAQ=LU,$ whereL,U, andP are defined as before, andQ is a permutation matrix that reorders the columns ofA.^[5]

Lower-diagonal-upper (LDU) decomposition

[edit]

Alower-diagonal-upper decomposition (LDU) is a decomposition of the form $A=LDU,$ whereD is adiagonal matrix, andL andU areunitriangular matrices, meaning that all the entries on the diagonals ofL andU are one.

Rectangular matrices

[edit]

Above we required thatA be a square matrix, but these decompositions can all be generalized to rectangular matrices as well.^[6] In that case,L andD are square matrices both of which have the same number of rows asA, andU has exactly the same dimensions asA. 'Upper triangular' should be interpreted as having only zero entries below the main diagonal, which starts at the upper left corner. Similarly, the more precise term forU is that it is therow echelon form of the matrixA.

Example

[edit]

We factor the following2 × 2 matrix: ${\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}}.$

One way to find the LU decomposition of this simple matrix would be to simply solve the linear equations by inspection. Expanding thematrix multiplication gives $\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.$

This system of equations isunderdetermined. In this case any two nonzero elements ofL andU matrices are parameters of the solution and can be set arbitrarily to any nonzero value. Therefore, to find the unique LU decomposition, it is necessary to put some restriction onL andU matrices. For example, we can conveniently require the lower triangular matrixL to be a unit triangular matrix, so that all the entries of its main diagonal are set to one. Then the system of equations has the following solution: ${\begin{aligned}\ell _{11}&=\ell _{22}=1\\\ell _{21}&=1.5\\u_{11}&=4\\u_{12}&=3\\u_{22}&=-1.5\end{aligned}}$

Substituting these values into the LU decomposition above yields ${\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

[edit]

Square matrices

[edit]

Any square matrixA admits LUP and PLU factorizations.^[3] IfA isinvertible, then it admits an LU (or LDU) factorization if and only if all its leading principalminors are nonzero^[7]^[8] (for example $\left[{\begin{smallmatrix}0&1\\1&0\end{smallmatrix}}\right]$ does not admit an LU or LDU factorization). IfA is asingular matrix of rankk, then it admits an LU factorization if the firstk leading principal minors are nonzero, although the converse is not true.^[9]

If a square, invertible matrix has an LDU factorization (with all diagonal entries ofL andU equal to1), then the factorization is unique.^[8] In that case, the LU factorization is also unique if we require that the diagonal of eitherL orU consists of ones.

In general, any square matrixA_n×n could have one of the following:

a unique LU factorization (as mentioned above);
infinitely many LU factorizations if any of the first(n − 1) columns are linearly dependent;
no LU factorization if the first(n − 1) columns are linearly independent and at least one leading principal minor is zero.

In Case 3, one can approximate an LU factorization by changing a diagonal entrya_ij toa_ij ±ε to avoid a zero leading principal minor.^[10]

Symmetric positive-definite matrices

[edit]

IfA is a symmetric (orHermitian, ifA is complex)positive-definite matrix, we can arrange matters so thatU is theconjugate transpose ofL. That is, we can writeA as $A=LL^{*}\,.$

This decomposition is called theCholesky decomposition. IfA is positive definite, then the Cholesky decomposition exists and is unique. Furthermore, computing the Cholesky decomposition is more efficient andnumerically more stable than computing some other LU decompositions.

General matrices

[edit]

For a (not necessarily invertible) matrix over any field, the exact necessary and sufficient conditions under which it has an LU factorization are known. The conditions are expressed in terms of the ranks of certain submatrices. The Gaussian elimination algorithm for obtaining LU decomposition has also been extended to this most general case.^[11]

Algorithms

[edit]

Closed formula

[edit]

When an LDU factorization exists and is unique, there is a closed (explicit) formula for the elements ofL,D, andU in terms of ratios of determinants of certain submatrices of the original matrixA.^[12] In particular,D₁ =A_1,1, and fori = 2, ... ,n,D_i is the ratio of thei-th principal submatrix to the(i − 1)-th principal submatrix. Computation of the determinants iscomputationally expensive, so this explicit formula is not used in practice.

Using Gaussian elimination

[edit]

The following algorithm is essentially a modified form ofGaussian elimination. Computing an LU decomposition using this algorithm requires⁠2/3⁠n³ floating-point operations, ignoring lower-order terms. Partialpivoting adds only a quadratic term; this is not the case for full pivoting.^[13]

Generalized explanation

[edit]

Notation

[edit]

Given anN × N matrix $A=(a_{i,j})_{1\leq i,j\leq N}$ , define $A^{(0)}$ as the original, unmodified version of the matrixA. The parenthetical superscript (e.g.,(0)) of the matrixA is the version of the matrix. The matrixA⁽ⁿ⁾ is theA matrix in which the elements below themain diagonal have already been eliminated to 0 through Gaussian elimination for the firstn columns.

Below is a matrix to observe to help us remember the notation (where each∗ represents anyreal number in the matrix):

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

Procedure

[edit]

During this process, we gradually modify the matrixA using row operations until it becomes the matrixU in which all the elements below the main diagonal are equal to zero. During this, we will simultaneously create two separate matricesP andL, such thatPA =LU.

We define the finalpermutation matrixP as the identity matrix which has all the same rows swapped in the same order as theA matrix while it transforms into the matrixU. For our matrixA⁽ⁿ⁻¹⁾, we may start by swapping rows to provide the desired conditions for then-th column. For example, we might swap rows to perform partial pivoting, or we might do it to set the pivot elementa_n,n on the main diagonal to a nonzero number so that we can complete the Gaussian elimination.

For our matrixA⁽ⁿ⁻¹⁾, we want to set every element below $a_{n,n}^{(n-1)}$ to zero (where $a_{n,n}^{(n-1)}$ is the element in then-th column of the main diagonal). We will denote each element below $a_{n,n}^{(n-1)}$ as $a_{i,n}^{(n-1)}$ (wherei =n+1, ... ,N). To set $a_{i,n}^{(n-1)}$ to zero, we setrow_i =row_i − (ℓ_i,n)⋅row_n for each rowi. For this operation, ${\textstyle \ell _{i,n}:={a_{i,n}^{(n-1)}}/{a_{n,n}^{(n-1)}}}$ . Once we have performed the row operations for the firstN − 1 columns, we have obtained anupper triangular matrixA^(N−1) which is denoted byU.

We can also create thelower triangular matrix denoted asL, by directly inputting the previously calculated values ofℓ_i,n via the formula below.

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

Example

[edit]

If we are given the matrix $A={\begin{pmatrix}0&5&{\frac {22}{3}}\\4&2&1\\2&7&9\\\end{pmatrix}},$ we will choose to implement partial pivoting and thus swap the first and second row so that our matrixA and the first iteration of ourP matrix respectively become $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}}.$ Once we have swapped the rows, we can eliminate the elements below the main diagonal on the first column by performing ${\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}}$ such that, ${\begin{alignedat}{0}\ell _{2,1}={\frac {0}{4}}=0\\\ell _{3,1}={\frac {2}{4}}=0.5\end{alignedat}}$ Once these rows have been subtracted, we have derived fromA⁽¹⁾ the matrix $A^{(1)}={\begin{pmatrix}4&2&1\\0&5&{\frac {22}{3}}\\0&6&8.5\\\end{pmatrix}}.$

Because we are implementing partial pivoting, we swap the second and third rows of our derived matrix and the current version of ourP matrix respectively to obtain $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}}.$ Now, we eliminate the elements below the main diagonal on the second column by performingrow₃ =row₃ − (ℓ_3,2)⋅row₂ such thatℓ_3,2 =⁠5/6⁠. Because no nonzero elements exist below the main diagonal in our current iteration ofA after this row subtraction, this row subtraction derives our finalA matrix (denoted asU) and finalP matrix: $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}}.$ After also switching the corresponding rows, we obtain our finalL matrix: $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}}$

Now these matrices have a relation such thatPA =LU.

Relations when no rows are swapped

[edit]

If we did not swap rows at all during this process, we can perform the row operations simultaneously for each columnn by setting $A^{(n)}:=L_{n}^{-1}A^{(n-1)},$ where $L_{n}^{-1}$ is theN × N identity matrix with itsn-th column replaced by thetransposed vector(0 ⋯ 0 1 −ℓ_n+1,n ⋯ −ℓ_N,n)^T.

In other words, the lower triangular matrix $L_{n}^{-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&-\ell _{N,n}&&&1\end{pmatrix}}.$

Performing all the row operations for the firstN − 1 columns using the $A^{(n)}:=L_{n}^{-1}A^{(n-1)}$ formula is equivalent to finding the decomposition $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)}.$ DenoteL =L₁ ⋯L_N−1 so thatA =LA^(N−1) =LU.

Now let's compute the sequence ofL₁ ⋯L_N−1. We know thatL_i has the following formula: $L_{n}={\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 nonzero item below the main diagonal in the same column as the other, then we can include all nonzero 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, multiplyL₁ together and generate the fused matrix denoted asL (as previously mentioned). Using the matrixL, we obtainA =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ℓ_i,n. If this assumption fails at some point, one needs to interchangen-th row with another row below it before continuing. This is why an LU decomposition in general looks likeP⁻¹A =LU.

LU Banachiewicz decomposition

[edit]

Illustration of operation of Banachiewicz LU algorithm to obtain 3-rd row and column of respectively matricesU andL. Involved matrices are named above squares marking their content. Matrix products and subtractions are applied only to elements in the thick frame boxes. Green filled thin frame boxes indicate values already known, from previous stages. Blue boxes indicate places inU andL matrices for storing of results.

Although Banachiewicz (1938) LU decomposition algorithm preceded the advent of programmed electronic computers, it wasready made for direct implementation into code as index swapping, transpose and column by column multiplication remain native built capabilities of the most programming languages and they are handled by compilers alone with little delay of the actual execution. The peculiar matrix notation used by Banachiewicz enabled him to multiply matrices column by column, a convenient feature for mechanical calculations as he could reveal consecutive factors by sliding a ruler to next rows of matrices. For human readers however, his equations are best transformed into standard matrix notation. To obtain from a full matrixA triangle matricesU andL calculations, start by copying top row and leftmost column ofA respectively into corresponding positions of matricesU andL. The known unit diagonal elements ofL are not stored nor used throughout the whole process.Next calculations continue for the subsequent rows and columns till the bottom right corner ofA.

The figure illustrates calculations for 3rd row and column, assuming previous stages were already completed. Involved matrices are named above squares marking their content. Matrix products and subtractions are applied only to elements in the thick frame boxes.Green filled thin frame boxes indicate values already known, from previous stages. Blue boxes indicate places inU andL matrices for storing of results. Note that at each stage the result elements ofL need to be divided by the corresponding pivot element on the main diagonal ofU. This applies to the leftmost column ofL too.

Note that after completion of 3rd stage the involved elements of matrixA are no longer used and neither those from the previous stages. This enables replacement of these elements with the result values ofU andL, i.e. execution of LU decompositionin place, so that the wholeA is replaced withU andL except for the unit diagonal ofL. Banachiewicz LU algorithm is well suited for partial pivoting by choosing the absolute maximum pivot from the newly calculated row ofU and subsequently swapping its columns so that it lands on the main diagonal. More details can be figured out from inspection of the enclosed Fortran90 code.

All partial pivoting LU algorithms cost roughly the same amount, of order ${\textstyle O\left({2 \over 3}n^{3}\right)}$ operations, wheren is the number of rows or columns ofA.

LU Crout decomposition

[edit]

Note that the decomposition obtained through this procedure is aDoolittle decomposition: the main diagonal ofL is composed solely of ones. If one would proceed by removing elementsabove the main diagonal by adding multiples of thecolumns (instead of removing elementsbelow the diagonal by adding multiples of therows), we would obtain aCrout decomposition, where the main diagonal ofU is of ones.

Another (equivalent) way of producing a Crout decomposition of a given matrixA is to obtain a Doolittle decomposition of the transpose ofA. Indeed, ifA^T =L₀U₀ is the LU-decomposition obtained through the algorithm presented in this section, then by takingL =U^T
₀ andU =L^T
₀, we have thatA =LU is a Crout decomposition.

Randomized algorithm

[edit]

It is possible to find a low rank approximation to an LU decomposition using arandomized algorithm. Given an input matrixA and a desired low rankk, the randomized LU returns permutation matricesP,Q and lower/upper trapezoidal matricesL,U of sizem × k andk × n respectively, such that with high probability‖PAQ −LU‖₂ ≤Cσ_k+1, whereC is a constant that depends on the parameters of the algorithm andσ_k+1 is the(k+1)-th singular value of the input matrixA.^[14]

Theoretical complexity

[edit]

If two matrices of ordern can be multiplied in timeM(n), whereM(n) ≥n^a for somea > 2, then an LU-decomposition can be computed in timeO(M(n)).^[15] This means, for example, that anO(n^2.376) algorithm exists based on theCoppersmith–Winograd algorithm.

Sparse-matrix decomposition

[edit]

Special algorithms have been developed for factorizing largesparse matrices. These algorithms attempt to find sparse factorsL andU. 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 nonzero value during the execution of an algorithm).

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

Applications

[edit]

Solving linear equations

[edit]

Given a system of linear equations in matrix form $A\mathbf {x} =\mathbf {b} ,$

we want to solve the equation forx, givenA andb. Suppose we have already obtained the LUP decomposition ofA such thatPA =LU, soLUx =Pb.

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

First, we solve the equationLy =Pb fory.
Second, we solve the equationUx =y forx.

In both cases we are dealing with triangular matrices (L andU), which can be solved directly byforward and backward substitution without using theGaussian 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 differentb. In this case it is faster (and more convenient) to do an LU decomposition of the matrixA once and then solve the triangular matrices for the differentb, rather than using Gaussian elimination each time. The matricesL andU could be thought to have "encoded" the Gaussian elimination process.

The cost of solving a system of linear equations is approximately⁠2/3⁠n³ floating-point operations if the matrixA has sizen. This makes it twice as fast as algorithms based onQR decomposition, which costs about⁠4/3⁠n³ floating-point operations whenHouseholder reflections are used. For this reason, LU decomposition is usually preferred.^[16]

Inverting a matrix

[edit]

When solving systems of equations,b is usually treated as a vector with a length equal to the height of matrixA. In matrix inversion however, instead of vectorb, we have matrixB, whereB is ann × p matrix, so that we are trying to find a matrixX (also an × p matrix): $AX=LUX=B.$

We can use the same algorithm presented earlier to solve for each column of matrixX. Now suppose thatB is the identity matrix of sizen,I_n. It follows that the resultX must be the inverse ofA.^[17]

Computing the determinant

[edit]

Given the LUP decompositionA =P⁻¹LU of a square matrixA, the determinant ofA can be computed straightforwardly as ${\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}}$

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 whereS is the number of row exchanges in the decomposition.

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

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

History

[edit]

LU decomposition: LU factors and their product in original Banachiewicz(1938) matrix notation

The LU decomposition is related to elimination of linear systems of equations, as e.g. described by Ralston.^[18] The solution ofN linear equations inN unknowns by elimination was already known to ancient Chinese.^[19] Before Gauss many mathematicians in Eurasia were performing and perfecting it yet as the method became relegated to school grade, few of them left any detailed descriptions. Thus the name Gaussian elimination is only a convenient abbreviation of a complex history.

The Polish astronomer Tadeusz Banachiewicz introduced the LU decomposition in 1938.^[1] About Banachiewicz, Paul Dwyer stated:^[20]^{[page needed]}

It appears that Gauss and Doolittle applied the method [of elimination] only to symmetric equations. More recent authors, for example, Aitken, Banachiewicz, Dwyer, and Crout ... have emphasized the use of the method, or variations of it, in connection with non-symmetric problems ... Banachiewicz ... saw the point ... that the basic problem is really one of matrix factorization, or "decomposition" as he called it.

— Paul Dwyer,Linear Computations (1951)

Banachiewicz^[1] was the first to consider elimination in terms of matrices and in this way formulated LU decomposition, as demonstrated by his graphic illustration. His calculations follow ordinary matrix ones, yet notation deviates in that he preferred to write one factor transposed, to be able to multiply them mechanically column by column, by sliding ruler down consecutive rows of both (usingarithmometer). Combined with swapped order of indices his formulae in modern notation read ${\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}}$

whereIA →A^T;x ≡ [x₁, ... ,x_n, −1];A′ refers toA extended with the last column; and the last component ofx is−1. Matrix formulae to calculate rows and columns of LU factors by recursion are given in the remaining part of Banachiewicz's paper as Eq. (2.3) and (2.4)(see§ Fortran90 code example below). This paper by Banachiewicz contains both derivation ofLU andR^TR factors of respectively non-symmetric and symmetric matrices. They are sometimes confused as later publications tend to tie his name solely with the rediscovery of Cholesky decomposition. Banachiewicz himself can be excused of inaction as already next year he suffered from persecution by occupiers, spending three month inthe Sachsenhausen Concentration Camp, on release from which he carried himself from a train his collaborator and co-prisoner Antoni Wilk, who died of exhaustion a week later.

Code examples

[edit]

Fortran90 code example

[edit]

ModulemluImplicit NoneInteger,Parameter::SP=Kind(1d0)! set I/O real precisionPrivate      Publicluban,lusolveContains      Subroutineluban(a,tol,g,h,ip,condinv,detnth)! By Banachiewicz (1938, hereafter B38) LU decomposition method calculates such! triangles L=G^T, and U=H  that square B=A^T=G^TH=LU. Partial pivoting! by column permutation IP(:) is modern addition.! Within the code a, g correspond to B38 A^T and G^T, so that a=gh holds.!! Normal use is for square A, however for RHS l already known! input of (A|l)^T yields (L|y^T)^T where x in L^Tx=y is solution of Ax=l.Real(SP),Intent(In)::a(:,:)! input matrix A(m,n), n<=mReal(SP),Intent(In)::tol! tolerance for near zero pivotReal(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)! note U columns are permutedReal(SP),Intent(Out)::condinv! 1/cond(A), 0 for singular AReal(SP),Intent(Out)::detnth! sign*Abs(det(A))**(1/n)Integer,Intent(Out)::ip(size(a,dim=2))! columns permutation!Integer::k,n,j,l,isigReal(SP)::tol0,pivmax,pivmin,piv!n=size(a,dim=2)tol0=Max(tol,3._SP*epsilon(tol0))! use default for tol=0!! Rectangular A and G are permitted under condition:If(n>size(a,dim=1).Or.n<1)Stop91Forall(k=1:n)ip(k)=kh=0._SPg=0._SPisig=1detnth=0._SPpivmax=Maxval(Abs(a(1,:)))pivmin=pivmax!Dok=1,n! Banachiewicz (1938) Eq. (2.3)h(k,ip(k:))=a(k,ip(k:))-Matmul(g(k,:k-1),h(:k-1,ip(k:)))!! Find row pivotj=(Maxloc(Abs(h(k,ip(k:))),dim=1)+k-1)If(j/=k)Then! Swap columns j and kisig=-isig! Change Det(A) sign because of permutationl=ip(k)ip(k)=ip(j)ip(j)=lEnd Ifpiv=Abs(h(k,ip(k)))pivmax=Max(piv,pivmax)! Adjust condinvpivmin=Min(piv,pivmin)If(piv<tol0)Then! singular matrixisig=0pivmax=1._SPExit            Else! Account for pivot contribution to Det(A) sign and valueIf(h(k,ip(k))<0._SP)isig=-isigdetnth=detnth+Log(piv)End If!! Transposed Banachiewicz (1938) Eq. (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._SPEnd Do!detnth=isig*Exp(detnth/n)condinv=Abs(isig)*pivmin/pivmax! Test for square A(n,n) by uncommenting below!         Print *, '|AQ-LU| ',Maxval (Abs(a(:,ip(:))-Matmul(g, h(:,ip(:)))))End SubroutinelubanSubroutinelusolve(l,u,ip,x)! Solves Ax=b system using triangle factors LU=AReal(SP),Intent(In)::l(:,:)! lower triangle matrix L(n,n)Real(SP),Intent(In)::u(:,:)! upper triangle matrix U(n,n)Integer,Intent(In)::ip(:)! columns permutation IP(n)Real(SP),Intent(InOut)::x(:,:)! Input: m sets of RHSs B(n,m),! Output: the corresponding sets of unknowns X(n,m)Integer::n,m,i,jn=size(ip)m=size(x,dim=2)If(n<1.Or.m<1.Or.Any([n,n]/=shape(l)).Or.Any(shape(l)/=shape(u)).Or.&n/=size(x,dim=1))Stop91Doi=1,mDoj=1,nx(j,i)=x(j,i)-dot_product(x(:j-1,i),l(j,:j-1))End Do            Doj=n,1,-1x(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 SubroutinelusolveEnd Modulemlu

C code example

[edit]

/* 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 */intLUPDecompose(double**A,intN,doubleTol,int*P){inti,j,k,imax;doublemaxA,*ptr,absA;for(i=0;i<=N;i++)P[i]=i;//Unit permutation matrix, P[N] initialized with Nfor(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)return0;//failure, matrix is degenerateif(imax!=i){//pivoting Pj=P[i];P[i]=P[imax];P[imax]=j;//pivoting rows of Aptr=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];}}return1;//decomposition done}/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension * OUTPUT: x - solution vector of A*x=b */voidLUPSolve(double**A,int*P,double*b,intN,double*x){for(inti=0;i<N;i++){x[i]=b[P[i]];for(intk=0;k<i;k++)x[i]-=A[i][k]*x[k];}for(inti=N-1;i>=0;i--){for(intk=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 */voidLUPInvert(double**A,int*P,intN,double**IA){for(intj=0;j<N;j++){for(inti=0;i<N;i++){IA[i][j]=P[i]==j?1.0:0.0;for(intk=0;k<i;k++)IA[i][j]-=A[i][k]*IA[k][j];}for(inti=N-1;i>=0;i--){for(intk=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 */doubleLUPDeterminant(double**A,int*P,intN){doubledet=A[0][0];for(inti=1;i<N;i++)det*=A[i][i];return(P[N]-N)%2==0?det:-det;}

C# code example

[edit]

publicclassSystemOfLinearEquations{publicdouble[]SolveUsingLU(double[,]matrix,double[]rightPart,intn){// decomposition of matrixdouble[,]lu=newdouble[n,n];doublesum=0;for(inti=0;i<n;i++){for(intj=i;j<n;j++){sum=0;for(intk=0;k<i;k++)sum+=lu[i,k]*lu[k,j];lu[i,j]=matrix[i,j]-sum;}for(intj=i+1;j<n;j++){sum=0;for(intk=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 = bdouble[]y=newdouble[n];for(inti=0;i<n;i++){sum=0;for(intk=0;k<i;k++)sum+=lu[i,k]*y[k];y[i]=rightPart[i]-sum;}// find solution of Ux = ydouble[]x=newdouble[n];for(inti=n-1;i>=0;i--){sum=0;for(intk=i+1;k<n;k++)sum+=lu[i,k]*x[k];x[i]=(1/lu[i,i])*(y[i]-sum);}returnx;}}

MATLAB code example

[edit]

functionLU=LUDecompDoolittle(A)n=length(A);LU=A;fork=2:nfori=1:k-1lamda=LU(k,i)/LU(i,i);LU(k,i)=lamda;LU(k,i+1:n)=LU(k,i+1:n)-LU(i,i+1:n)*lamda;endendendfunctionx=SolveLinearSystem(LU, B)n=length(LU);y=zeros(size(B));% find solution of Ly = Bfori=1:ny(i,:)=B(i,:)-LU(i,1:i)*y(1:i,:);end% find solution of Ux = yx=zeros(size(B));fori=n:(-1):1x(i,:)=(y(i,:)-LU(i,(i+1):n)*x((i+1):n,:))/LU(i,i);endendA=[433;633;343]LU=LUDecompDoolittle(A)B=[123;456;789;101112]'x=SolveLinearSystem(LU,B)A*x

Notes

[edit]

^^a ^b ^cBanachiewicz (1938).
^Cormen et al. (2009), p. 819, 28.1: Solving systems of linear equations.
^^a ^bOkunev & Johnson (1997), Corollary 3.
^Trefethen & Bau (1997), p. 166.
^Trefethen & Bau (1997), p. 161.
^Banachiewicz (1938);Lay, Lay & McDonald (2021), p. 133, 2.5: Matrix Factorizations.
^Rigotti (2001), Leading Principal Minor.
^^a ^bHorn & Johnson (1985), Corollary 3.5.5
^Horn & Johnson (1985), Theorem 3.5.2.
^Nhiayi, Ly; Phan-Yamada, Tuyetdong (2021). "Examining Possible LU Decomposition".North American GeoGebra Journal.9 (1).
^Okunev & Johnson (1997).
^Householder (1975).
^Golub & Van Loan (1996), pp. 112, 119.
^Shabat, Gil; Shmueli, Yaniv; Aizenbud, Yariv; Averbuch, Amir (2016). "Randomized LU Decomposition".Applied and Computational Harmonic Analysis.44 (2):246–272.arXiv:1310.7202.doi:10.1016/j.acha.2016.04.006.S2CID 1900701.
^Bunch & Hopcroft (1974).
^Trefethen & Bau (1997), p. 152.
^Golub & Van Loan (1996), p. 121.
^Ralston (1965).
^Hart (2011).
^Dwyer (1951).

References

[edit]

Banachiewicz, T. (1938),"Méthode de résolution numérique des équations linéaires …"(PDF),Bull. Intern. De l'Acad. Polonaise, Série A. Sc. Math.:393–404.
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,hdl:1813/6003,ISSN 0025-5718,JSTOR 2005828.
Cormen, Thomas H.;Leiserson, Charles E.;Rivest, Ronald L.;Stein, Clifford (2009),Introduction to Algorithms (3rd ed.), MIT Press and McGraw-Hill,ISBN 978-0-262-03293-3.
Dwyer, Paul S. (1951),Linear Computations, New York: Wiley.
Golub, Gene H.;Van Loan, Charles F. (1996),Matrix Computations (3rd ed.), Baltimore: Johns Hopkins,ISBN 978-0-8018-5414-9.
Hart, Roger (2011),The Chinese Roots of Linear Algebra, Baltimore: Johns Hopkins,ISBN 978-0801897559.
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.
Lay, David C.; Lay, Steven R.; McDonald, Judi J. (2021),Linear Algebra and its Applications (Sixth ed.), Pearson,ISBN 978-0-13-585125-8.
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.
Ralston, Anthony (1965),A First Course in Numerical Analysis, New York: McGraw-Hill, Inc,ISBN 978-0-070-51157-6.
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

[edit]

References

LU decomposition onMathWorld.
LU decomposition onMath-Linux.
LU decomposition atHolistic 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

v t e Numerical linear algebra
Key concepts	Floating point Numerical stability
Problems	System of linear equations Matrix decompositions Matrix multiplication (algorithms) Matrix splitting Sparse problems
Hardware	CPU cache TLB Cache-oblivious algorithm SIMD Multiprocessing
Software	ATLAS MATLAB Basic Linear Algebra Subprograms (BLAS) LAPACK Specialized libraries General purpose software