# Chapter 75. LU Decomposition # Chapter 75. LU Decomposition LU decomposition is a factorization of a matrix into a lower triangular part and an upper triangular part. It is one of the standard matrix decompositions used to solve linear systems, compute determinants, and organize Gaussian elimination in matrix form. For a square matrix \(A\), an LU decomposition has the form $$ A = LU, $$ where \(L\) is lower triangular and \(U\) is upper triangular. A lower triangular matrix has zero entries above the main diagonal. An upper triangular matrix has zero entries below the main diagonal. This is the standard lower-upper factorization described in numerical linear algebra references. ## 75.1 Triangular Matrices A lower triangular matrix has the form $$ L = \begin{bmatrix} \ell_{11} & 0 & 0 & \cdots & 0 \\ \ell_{21} & \ell_{22} & 0 & \cdots & 0 \\ \ell_{31} & \ell_{32} & \ell_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \ell_{n1} & \ell_{n2} & \ell_{n3} & \cdots & \ell_{nn} \end{bmatrix}. $$ An upper triangular matrix has the form $$ U = \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & u_{nn} \end{bmatrix}. $$ Triangular systems are easy to solve because each equation contains only a controlled set of unknowns. A lower triangular system is solved from top to bottom. An upper triangular system is solved from bottom to top. This is the reason LU decomposition is useful. It replaces one general linear system by two triangular systems. ## 75.2 The Basic Factorization Let $$ A = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}. $$ An LU decomposition writes \(A\) as $$ A = LU. $$ In many treatments, \(L\) is chosen to have ones on the diagonal: $$ L = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ \ell_{21} & 1 & 0 & \cdots & 0 \\ \ell_{31} & \ell_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \ell_{n1} & \ell_{n2} & \ell_{n3} & \cdots & 1 \end{bmatrix}. $$ Such a matrix is called unit lower triangular. With this convention, the diagonal scaling is placed in \(U\). This convention is often called the Doolittle form of LU decomposition. A related convention, called the Crout form, places the unit diagonal in \(U\) instead. ## 75.3 LU and Gaussian Elimination LU decomposition is Gaussian elimination recorded as a matrix factorization. Suppose Gaussian elimination transforms \(A\) into an upper triangular matrix \(U\). At each step, a multiple of one row is subtracted from a lower row to eliminate an entry below the pivot. The multipliers used in these eliminations become the entries of \(L\). For example, to eliminate \(a_{21}\) using \(a_{11}\), one uses the multiplier $$ \ell_{21} = \frac{a_{21}}{a_{11}}. $$ To eliminate \(a_{31}\), one uses $$ \ell_{31} = \frac{a_{31}}{a_{11}}. $$ These multipliers are stored below the diagonal in \(L\). The remaining matrix after elimination is \(U\). Thus \(U\) records the result of elimination, while \(L\) records how the elimination was performed. ## 75.4 A Two by Two Example Let $$ A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}. $$ We seek $$ A = LU, $$ with $$ L = \begin{bmatrix} 1 & 0 \\ \ell_{21} & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} u_{11} & u_{12} \\ 0 & u_{22} \end{bmatrix}. $$ Multiplying gives $$ LU = \begin{bmatrix} u_{11} & u_{12} \\ \ell_{21}u_{11} & \ell_{21}u_{12} + u_{22} \end{bmatrix}. $$ Equating this with \(A\), we obtain $$ u_{11} = 2, \qquad u_{12} = 3, $$ $$ \ell_{21}u_{11} = 4, $$ so $$ \ell_{21} = 2. $$ The final entry gives $$ \ell_{21}u_{12} + u_{22} = 7. $$ Therefore $$ 2 \cdot 3 + u_{22} = 7, $$ so $$ u_{22} = 1. $$ Hence $$ L = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix}. $$ Indeed, $$ \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}. $$ ## 75.5 Solving Linear Systems with LU Suppose $$ Ax = b $$ and $$ A = LU. $$ Then $$ LUx = b. $$ Introduce an intermediate vector $$ y = Ux. $$ Then the original system becomes two triangular systems: $$ Ly = b $$ and $$ Ux = y. $$ The first system is solved by forward substitution. The second system is solved by backward substitution. This procedure is especially useful when many systems have the same coefficient matrix \(A\) but different right-hand sides \(b\). The factorization \(A = LU\) is computed once. Each new right-hand side then requires only two triangular solves. ## 75.6 Forward Substitution Let $$ Ly = b, $$ where \(L\) is lower triangular. Written componentwise, $$ \ell_{11}y_1 = b_1, $$ $$ \ell_{21}y_1 + \ell_{22}y_2 = b_2, $$ $$ \ell_{31}y_1 + \ell_{32}y_2 + \ell_{33}y_3 = b_3, $$ and so on. The first equation gives \(y_1\). The second equation then gives \(y_2\). The third equation gives \(y_3\). Continuing in this order gives all entries of \(y\). For a general lower triangular matrix, $$ y_i = \frac{ b_i - \sum_{j=1}^{i-1} \ell_{ij}y_j }{ \ell_{ii} }. $$ If \(L\) is unit lower triangular, then \(\ell_{ii} = 1\), and the formula becomes $$ y_i = b_i - \sum_{j=1}^{i-1} \ell_{ij}y_j. $$ ## 75.7 Backward Substitution Let $$ Ux = y, $$ where \(U\) is upper triangular. Written componentwise, the last equation is $$ u_{nn}x_n = y_n. $$ This gives \(x_n\). The previous equation gives \(x_{n-1}\), since \(x_n\) is already known. The process continues upward. For a general upper triangular matrix, $$ x_i = \frac{ y_i - \sum_{j=i+1}^{n} u_{ij}x_j }{ u_{ii} }. $$ Backward substitution works only when the diagonal entries needed for division are nonzero. In the LU setting, this means the pivots must be nonzero. ## 75.8 Example: Solving a System Let $$ A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}, \qquad b = \begin{bmatrix} 8 \\ 18 \end{bmatrix}. $$ From the previous example, $$ A = LU, $$ where $$ L = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix}. $$ First solve $$ Ly = b. $$ That is, $$ \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} 8 \\ 18 \end{bmatrix}. $$ The equations are $$ y_1 = 8, $$ $$ 2y_1 + y_2 = 18. $$ Therefore $$ y_2 = 18 - 16 = 2. $$ So $$ y = \begin{bmatrix} 8 \\ 2 \end{bmatrix}. $$ Now solve $$ Ux = y. $$ That is, $$ \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 8 \\ 2 \end{bmatrix}. $$ The second equation gives $$ x_2 = 2. $$ The first equation gives $$ 2x_1 + 3x_2 = 8. $$ Thus $$ 2x_1 + 6 = 8, $$ so $$ x_1 = 1. $$ The solution is $$ x = \begin{bmatrix} 1 \\ 2 \end{bmatrix}. $$ ## 75.9 Existence of LU Decomposition An LU decomposition without row exchanges does not always exist. The elimination process requires nonzero pivots. If a pivot is zero, division by that pivot cannot be performed. For example, $$ A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} $$ cannot be decomposed as \(A = LU\) with \(L\) unit lower triangular and \(U\) upper triangular. The first pivot would be zero. However, row exchanges fix this problem. Instead of factoring \(A\) directly, one factors a row-permuted version of \(A\): $$ PA = LU. $$ Here \(P\) is a permutation matrix. This form is called LU decomposition with pivoting, or LUP decomposition. Standard numerical algorithms usually include pivoting because it avoids zero pivots and improves numerical stability. ## 75.10 Pivoting Pivoting means exchanging rows before or during elimination. The most common form is partial pivoting. At each step, the algorithm chooses a pivot row by looking down the current column and selecting an entry with large absolute value. That row is swapped into the pivot position. The factorization then has the form $$ PA = LU. $$ Equivalently, $$ A = P^{-1}LU. $$ Since a permutation matrix satisfies $$ P^{-1} = P^T, $$ one may also write $$ A = P^TLU. $$ Pivoting has two purposes. First, it prevents division by zero when a possible nonzero pivot is available. Second, it can reduce the growth of rounding errors in floating point arithmetic. ## 75.11 Determinants from LU LU decomposition also simplifies determinant computation. If $$ A = LU, $$ then $$ \det(A) = \det(L)\det(U). $$ The determinant of a triangular matrix is the product of its diagonal entries. Therefore, $$ \det(A) = \left(\prod_{i=1}^{n} \ell_{ii}\right) \left(\prod_{i=1}^{n} u_{ii}\right). $$ If \(L\) is unit lower triangular, then $$ \det(L) = 1. $$ Hence $$ \det(A) = \prod_{i=1}^{n} u_{ii}. $$ If pivoting is used and $$ PA = LU, $$ then $$ \det(P)\det(A) = \det(L)\det(U). $$ Thus $$ \det(A) = \det(P)\det(L)\det(U), $$ because \(\det(P) = \pm 1\). Each row swap changes the sign of the determinant. ## 75.12 Inverses from LU LU decomposition can be used to compute the inverse of a matrix, although direct inversion is often avoided in numerical work when only a solution of \(Ax=b\) is needed. To find \(A^{-1}\), solve $$ AX = I. $$ If $$ A = LU, $$ then $$ LUX = I. $$ This is equivalent to solving $$ LY = I $$ and then $$ UX = Y. $$ Each column of \(X\) is found by solving two triangular systems. The columns of \(X\) are the columns of \(A^{-1}\). This method is systematic, but it is usually better to solve \(Ax=b\) directly by LU rather than compute \(A^{-1}\) explicitly. ## 75.13 Algorithm The following algorithm computes an LU decomposition without pivoting, using the convention that \(L\) has ones on the diagonal. ``` lu_decomposition(A): n = number of rows of A L = identity matrix of size n U = zero matrix of size n by n for k = 1 to n: for j = k to n: U[k,j] = A[k,j] for s = 1 to k-1: U[k,j] = U[k,j] - L[k,s] * U[s,j] for i = k+1 to n: L[i,k] = A[i,k] for s = 1 to k-1: L[i,k] = L[i,k] - L[i,s] * U[s,k] L[i,k] = L[i,k] / U[k,k] return L, U ``` The division by \(U[k,k]\) requires a nonzero pivot. If the pivot is zero or very small, the algorithm should use pivoting. ## 75.14 Operation Count For a dense \(n \times n\) matrix, LU decomposition requires on the order of $$ \frac{2}{3}n^3 $$ arithmetic operations. Solving the two triangular systems requires on the order of $$ 2n^2 $$ operations for each right-hand side. Thus the expensive step is the factorization. Once the factorization is available, solving for additional right-hand sides is much cheaper. This explains why LU decomposition is useful in repeated-solve problems. ## 75.15 Interpretation LU decomposition separates elimination into two pieces. The matrix \(U\) is the upper triangular matrix obtained after elimination. The matrix \(L\) stores the multipliers that produced the eliminations. The equation $$ A = LU $$ says that the original matrix can be reconstructed from these two pieces. In computational terms, LU decomposition is a reusable form of Gaussian elimination. In algebraic terms, it is a factorization into simpler matrices. In geometric terms, it decomposes a linear transformation into transformations with triangular structure. ## 75.16 Common Variants | Variant | Form | Description | |---|---|---| | LU | \(A = LU\) | Factorization without pivoting | | LUP | \(PA = LU\) | LU with row pivoting | | PLU | \(A = PLU\) | Equivalent convention using a permutation factor | | LDU | \(A = LDU\) | Separates diagonal scaling into \(D\) | | Crout | \(A = LU\) | \(U\) has unit diagonal | | Doolittle | \(A = LU\) | \(L\) has unit diagonal | Different conventions place diagonal entries in different factors. The essential idea remains the same: a general matrix is represented by triangular factors. ## 75.17 Summary LU decomposition factors a matrix into lower and upper triangular matrices: $$ A = LU. $$ With pivoting, the more stable and general form is $$ PA = LU. $$ The decomposition is Gaussian elimination written as a matrix product. It is used to solve linear systems, compute determinants, and perform repeated solves efficiently. Its practical value comes from replacing one general system by two triangular systems: $$ Ly = b, \qquad Ux = y. $$ The triangular systems are solved by forward and backward substitution. This makes LU decomposition one of the central tools of computational linear algebra.