# Chapter 76. PLU Decomposition # Chapter 76. PLU Decomposition PLU decomposition is LU decomposition with row permutation. It factors a row-reordered matrix into a lower triangular matrix and an upper triangular matrix. The standard form is $$ PA = LU. $$ Here \(A\) is the original matrix, \(P\) is a permutation matrix, \(L\) is lower triangular, and \(U\) is upper triangular. The factor \(P\) records row exchanges. The factors \(L\) and \(U\) record Gaussian elimination after those row exchanges. PLU is used because ordinary LU may fail when a pivot is zero, while row exchanges can move a usable pivot into position. ## 76.1 Motivation LU decomposition without pivoting assumes that the required pivots are nonzero. This assumption is too strong. Consider $$ A = \begin{bmatrix} 0 & 1 \\ 1 & 2 \end{bmatrix}. $$ The first pivot is \(0\). Gaussian elimination cannot divide by this pivot. However, the matrix itself is not defective for solving linear systems. We can exchange the two rows: $$ PA = \begin{bmatrix} 1 & 2 \\ 0 & 1 \end{bmatrix}. $$ This matrix is already upper triangular. The problem was the row order, not the matrix. PLU decomposition separates these two issues. The permutation matrix \(P\) chooses a workable row order. The triangular factors \(L\) and \(U\) then describe the elimination. ## 76.2 Permutation Matrices A permutation matrix is obtained by rearranging the rows of the identity matrix. Each row contains exactly one \(1\), and each column contains exactly one \(1\). All other entries are \(0\). For example, $$ P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ exchanges the first two rows of a matrix. If $$ A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32} \end{bmatrix}, $$ then $$ PA = \begin{bmatrix} a_{21} & a_{22} \\ a_{11} & a_{12} \\ a_{31} & a_{32} \end{bmatrix}. $$ Left multiplication by a permutation matrix permutes rows. Right multiplication by a permutation matrix permutes columns. In PLU decomposition, the usual purpose of \(P\) is row permutation. ## 76.3 The Factorization A PLU decomposition has the form $$ PA = LU. $$ The factors have the following roles. | Factor | Meaning | |---|---| | \(P\) | Records row exchanges | | \(L\) | Records elimination multipliers | | \(U\) | Records the upper triangular result | The matrix \(P\) is orthogonal in the real case: $$ P^{-1} = P^T. $$ Therefore the equation $$ PA = LU $$ can also be written as $$ A = P^TLU. $$ Some authors call this a PLU decomposition. Others reserve PLU for $$ A = PLU. $$ The difference is only convention. The essential point is that a permutation factor is included with the triangular factors. ## 76.4 Why Pivoting Is Needed Gaussian elimination proceeds by choosing a pivot, dividing by it, and using multiples of the pivot row to eliminate entries below it. At step \(k\), the pivot is the entry in position \((k,k)\) of the current matrix. If this entry is zero, elimination cannot proceed. If it is very small, elimination can proceed algebraically, but the computation may be unstable in floating point arithmetic. Pivoting exchanges rows so that a better pivot appears in the pivot position. In partial pivoting, the pivot row is chosen from the current column. The usual rule is to choose the row whose candidate pivot has largest absolute value. This is the most common version in practical LU factorization. ## 76.5 Partial Pivoting Suppose the elimination is at column \(k\). We examine the entries $$ a_{kk}, a_{k+1,k}, \ldots, a_{nk}. $$ Partial pivoting chooses an index \(r\) such that $$ |a_{rk}| = \max_{k \le i \le n} |a_{ik}|. $$ Then rows \(k\) and \(r\) are exchanged. After the row exchange, the pivot is as large as possible within the remaining part of the current column. This does not guarantee perfect numerical behavior in every artificial example, but it is reliable in ordinary dense computations and widely used in numerical software. ## 76.6 A Simple Example Let $$ A = \begin{bmatrix} 0 & 2 \\ 3 & 4 \end{bmatrix}. $$ Ordinary LU without pivoting fails because the first pivot is \(0\). Exchange the two rows. The permutation matrix is $$ P = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}. $$ Then $$ PA = \begin{bmatrix} 3 & 4 \\ 0 & 2 \end{bmatrix}. $$ This matrix is already upper triangular, so we may take $$ L = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 3 & 4 \\ 0 & 2 \end{bmatrix}. $$ Thus $$ PA = LU. $$ This example shows the purpose of \(P\). The permutation allows the triangular factorization to proceed. ## 76.7 A Three by Three Example Let $$ A = \begin{bmatrix} 0 & 2 & 1 \\ 4 & 8 & 3 \\ 2 & 6 & 2 \end{bmatrix}. $$ The first pivot is \(0\), so we exchange the first row with the second row. This gives $$ P_1A = \begin{bmatrix} 4 & 8 & 3 \\ 0 & 2 & 1 \\ 2 & 6 & 2 \end{bmatrix}. $$ Now eliminate the entry below the first pivot in row \(3\). The multiplier is $$ \ell_{31} = \frac{2}{4} = \frac{1}{2}. $$ Subtract \(\frac{1}{2}\) times row \(1\) from row \(3\): $$ \begin{bmatrix} 4 & 8 & 3 \\ 0 & 2 & 1 \\ 0 & 2 & \frac{1}{2} \end{bmatrix}. $$ At the second pivot position, the pivot is \(2\). The entry below it is also \(2\). No row exchange is necessary. The multiplier is $$ \ell_{32} = \frac{2}{2} = 1. $$ Subtract row \(2\) from row \(3\): $$ U = \begin{bmatrix} 4 & 8 & 3 \\ 0 & 2 & 1 \\ 0 & 0 & -\frac{1}{2} \end{bmatrix}. $$ The lower triangular factor is $$ L = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \frac{1}{2} & 1 & 1 \end{bmatrix}. $$ The permutation matrix that exchanges the first two rows is $$ P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}. $$ Therefore $$ PA = LU. $$ Indeed, $$ LU = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \frac{1}{2} & 1 & 1 \end{bmatrix} \begin{bmatrix} 4 & 8 & 3 \\ 0 & 2 & 1 \\ 0 & 0 & -\frac{1}{2} \end{bmatrix} = \begin{bmatrix} 4 & 8 & 3 \\ 0 & 2 & 1 \\ 2 & 6 & 2 \end{bmatrix} = PA. $$ ## 76.8 Solving Systems with PLU Suppose we want to solve $$ Ax = b. $$ If $$ PA = LU, $$ then multiply the equation \(Ax=b\) by \(P\): $$ PAx = Pb. $$ Since \(PA=LU\), this becomes $$ LUx = Pb. $$ Introduce $$ y = Ux. $$ Then solve the two triangular systems $$ Ly = Pb $$ and $$ Ux = y. $$ Thus PLU solving has three stages: | Stage | Operation | |---|---| | 1 | Permute the right-hand side: \(Pb\) | | 2 | Solve \(Ly = Pb\) by forward substitution | | 3 | Solve \(Ux = y\) by backward substitution | The row permutation must be applied to the right-hand side as well as the coefficient matrix. ## 76.9 Example: Solving a System Let $$ A = \begin{bmatrix} 0 & 2 \\ 3 & 4 \end{bmatrix}, \qquad b = \begin{bmatrix} 6 \\ 11 \end{bmatrix}. $$ We have $$ P = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \qquad L = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 3 & 4 \\ 0 & 2 \end{bmatrix}. $$ First compute $$ Pb = \begin{bmatrix} 11 \\ 6 \end{bmatrix}. $$ Since \(L\) is the identity matrix, $$ y = Pb = \begin{bmatrix} 11 \\ 6 \end{bmatrix}. $$ Now solve $$ Ux = y. $$ That is, $$ \begin{bmatrix} 3 & 4 \\ 0 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 11 \\ 6 \end{bmatrix}. $$ The second equation gives $$ 2x_2 = 6, $$ so $$ x_2 = 3. $$ The first equation gives $$ 3x_1 + 4x_2 = 11. $$ Substituting \(x_2=3\), $$ 3x_1 + 12 = 11, $$ so $$ x_1 = -\frac{1}{3}. $$ Therefore $$ x = \begin{bmatrix} -\frac{1}{3} \\ 3 \end{bmatrix}. $$ ## 76.10 Determinants from PLU PLU decomposition also gives an efficient determinant formula. If $$ PA = LU, $$ then $$ \det(P)\det(A) = \det(L)\det(U). $$ If \(L\) is unit lower triangular, then $$ \det(L)=1. $$ Since \(U\) is triangular, $$ \det(U)=\prod_{i=1}^{n} u_{ii}. $$ Therefore $$ \det(A) = \det(P)\prod_{i=1}^{n} u_{ii}. $$ The determinant of a permutation matrix is either \(1\) or \(-1\). Each row exchange changes the sign. If the factorization uses \(s\) row exchanges, then $$ \det(P)=(-1)^s. $$ Thus $$ \det(A) = (-1)^s \prod_{i=1}^{n} u_{ii}. $$ This formula is commonly used in numerical computation, often through the logarithm of the absolute determinant when the matrix is large or ill-scaled. ## 76.11 Invertibility A square matrix \(A\) is invertible exactly when all pivots in a complete PLU factorization are nonzero. If $$ PA = LU $$ and \(P\), \(L\), and \(U\) are all invertible, then \(A\) is invertible. The permutation matrix \(P\) is always invertible. A unit lower triangular matrix \(L\) is invertible. Therefore the decisive condition is whether \(U\) has nonzero diagonal entries. Since \(U\) is triangular, $$ U \text{ is invertible} $$ if and only if $$ u_{11}u_{22}\cdots u_{nn} \ne 0. $$ Equivalently, $$ u_{ii} \ne 0 $$ for every \(i\). ## 76.12 Algorithm The following algorithm computes a PLU decomposition using partial pivoting. ```text plu_decomposition(A): n = number of rows of A P = identity matrix of size n L = zero matrix of size n by n U = copy of A for k = 1 to n: pivot = k for i = k to n: if abs(U[i,k]) > abs(U[pivot,k]): pivot = i if U[pivot,k] == 0: stop: matrix is singular if pivot != k: swap rows k and pivot in U swap rows k and pivot in P swap entries L[k,1:k-1] and L[pivot,1:k-1] L[k,k] = 1 for i = k+1 to n: L[i,k] = U[i,k] / U[k,k] U[i,k:n] = U[i,k:n] - L[i,k] * U[k,k:n] return P, L, U ``` The row swap in \(L\) affects only columns already computed. This detail is important. Those entries store previous elimination multipliers, and their row order must remain consistent with the current row order. ## 76.13 Storage in Practice Numerical libraries usually do not store \(P\), \(L\), and \(U\) as three full dense matrices. Instead, they commonly store: | Object | Practical storage | |---|---| | \(P\) | A pivot vector | | \(L\) | Strict lower part of an array | | \(U\) | Upper part of the same array | | Diagonal of \(L\) | Implicit ones | The pivot vector records which rows were exchanged. The array that originally held \(A\) is overwritten by the entries of \(L\) and \(U\). This storage method avoids unnecessary memory use. It also reflects the fact that \(L\) and \(U\) occupy complementary triangular parts of the matrix, except for the diagonal of \(L\), which is known to be all ones. ## 76.14 Operation Count For a dense \(n \times n\) matrix, PLU decomposition with partial pivoting has the same leading arithmetic cost as ordinary LU: $$ \frac{2}{3}n^3 $$ operations, ignoring lower-order terms. Pivot search and row exchanges add lower-order work. The search for pivots costs on the order of \(n^2\) comparisons. Row exchanges also cost on the order of \(n^2\) data movement. Thus pivoting improves robustness without changing the leading cubic cost of dense LU factorization. ## 76.15 Numerical Stability Pivoting is not only a device for avoiding zero pivots. It also helps control numerical error. If a pivot is very small, the multipliers $$ \ell_{ik} = \frac{a_{ik}}{a_{kk}} $$ can become very large. Large multipliers can magnify rounding errors and produce unreliable results. Partial pivoting chooses a large available pivot in the current column. This keeps the multipliers bounded in magnitude by \(1\), since $$ |\ell_{ik}| \le 1 $$ after pivoting within the column. This bound is one reason partial pivoting is effective in practice. It does not eliminate all possible instability, but it prevents the most immediate source of large elimination multipliers. ## 76.16 PLU and Row Echelon Form PLU decomposition is closely related to row echelon form. Gaussian elimination transforms a matrix into an upper triangular or row echelon matrix. PLU records the same process in factored form: $$ PA = LU. $$ The matrix \(U\) is the echelon-like result. The matrix \(L\) stores the elimination multipliers. The matrix \(P\) stores the row exchanges. Thus PLU is not a different method from Gaussian elimination. It is Gaussian elimination with memory. ## 76.17 Comparison with LU | Feature | LU | PLU | |---|---|---| | Form | \(A=LU\) | \(PA=LU\) | | Row exchanges | No | Yes | | Handles zero pivots | Not generally | Usually, if a nonzero pivot exists | | Numerical stability | Weaker | Stronger | | Stored permutation | None | Pivot vector or \(P\) | | Common in software | Less common alone | Standard dense solver form | Ordinary LU is useful when the matrix is known to allow elimination without row exchanges. PLU is the standard default for general dense matrices. ## 76.18 Common Conventions There are several equivalent ways to write the same idea. | Convention | Meaning | |---|---| | \(PA = LU\) | Permute rows of \(A\), then factor | | \(A = P^{-1}LU\) | Move the permutation to the other side | | \(A = P^TLU\) | Use \(P^{-1}=P^T\) | | \(A = PLU\) | Alternative convention with \(P\) defined differently | When reading formulas, check which side of \(A\) the permutation matrix appears on. The mathematical content is the same, but the stored permutation may represent either the performed row swaps or their inverse. ## 76.19 Summary PLU decomposition extends LU decomposition by allowing row exchanges: $$ PA = LU. $$ The permutation matrix \(P\) records pivoting. The lower triangular matrix \(L\) records elimination multipliers. The upper triangular matrix \(U\) records the result of elimination. PLU is the standard practical form of LU decomposition for general square matrices. It avoids zero pivots when possible, improves numerical stability, and provides an efficient way to solve systems: $$ Ly = Pb, \qquad Ux = y. $$ It also gives direct access to determinant computation, invertibility tests, and reusable factorizations for multiple right-hand sides.