Solving Permuted Systems From PA = LU

Why the Right-Hand Side Must Be Permuted

When elimination uses row exchanges, the factorization records those exchanges in a permutation matrix. The resulting identity is \( \mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}\), not simply \(\mathbf{A}=\mathbf{L}\mathbf{U}\). Consequently, solving \(\mathbf{A}\mathbf{x}=\mathbf{b}\) requires applying the same row permutation to \(\mathbf{b}\) before using the triangular factors. Forgetting this step, or applying the permutation in the wrong direction, can produce a vector that solves a different system.

Here \(\mathbf{A}\in\mathbb{R}^{n\times n}\), \(\mathbf{b},\mathbf{x}\in\mathbb{R}^n\), and \(\mathbf{P}\) is an \(n\times n\) permutation matrix: it is obtained by reordering the rows of the identity matrix. Multiplying a vector by \(\mathbf{P}\) reorders its entries in that same row order. The factors in a standard pivoted LU factorization are a unit lower-triangular matrix \(\mathbf{L}\), with diagonal entries all equal to \(1\), and an upper-triangular matrix \(\mathbf{U}\).

To see why the permutation acts on the right-hand side, multiply the original equation on the left by \(\mathbf{P}\). The result is \(\mathbf{P}\mathbf{A}\mathbf{x}=\mathbf{P}\mathbf{b}\). Replacing \(\mathbf{P}\mathbf{A}\) by \(\mathbf{L}\mathbf{U}\) gives \(\mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{P}\mathbf{b}\). Set \(\mathbf{y}=\mathbf{U}\mathbf{x}\). Then the equation separates into \(\mathbf{L}\mathbf{y}=\mathbf{P}\mathbf{b}\) and \(\mathbf{U}\mathbf{x}=\mathbf{y}\).

The Solve Procedure and Why It Works

For a lower-triangular matrix, forward substitution determines entries from top to bottom. If \(\mathbf{L}\) is unit lower-triangular, the equations in \(\mathbf{L}\mathbf{y}=\mathbf{P}\mathbf{b}\) give

For an upper-triangular matrix with nonzero diagonal entries, back substitution determines entries from bottom to top:

Proof. Assume the stated factorization and let \(\mathbf{b}\) be any vector in \(\mathbb{R}^n\). A permutation matrix reorders rows, so its rows are the standard basis vectors in some order. Their dot products are \(1\) with themselves and \(0\) with one another; therefore \(\mathbf{P}\mathbf{P}^{\top}=\mathbf{I}_n\), and \(\mathbf{P}^{-1}=\mathbf{P}^{\top}\). Thus multiplying an equation by \(\mathbf{P}\) preserves equivalence: \(\mathbf{A}\mathbf{x}=\mathbf{b}\) holds exactly when \(\mathbf{P}\mathbf{A}\mathbf{x}=\mathbf{P}\mathbf{b}\) holds. Using the factorization, this last equation is \(\mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{P}\mathbf{b}\).

Because \(\mathbf{L}\) is unit lower-triangular, its \(i\)th equation has coefficient \(1\) on \(y_i\), involves only \(y_1,\dots,y_i\), and determines \(y_i\) once \(y_1,\dots,y_{i-1}\) are known. Starting at \(i=1\), the forward-substitution formula therefore determines a unique \(\mathbf{y}\) satisfying \(\mathbf{L}\mathbf{y}=\mathbf{P}\mathbf{b}\). Because each \(u_{ii}\ne0\), the \(i\)th equation of \(\mathbf{U}\mathbf{x}=\mathbf{y}\) determines \(x_i\) from \(x_{i+1},\dots,x_n\). Starting at \(i=n\), back substitution determines a unique \(\mathbf{x}\). For that vector, \(\mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{L}\mathbf{y}=\mathbf{P}\mathbf{b}\), so \(\mathbf{P}\mathbf{A}\mathbf{x}=\mathbf{P}\mathbf{b}\). Multiplying by \(\mathbf{P}^{\top}\) gives \(\mathbf{A}\mathbf{x}=\mathbf{b}\). Any solution of the original system would also give a solution of both triangular systems, whose solutions are unique; hence the original solution is unique as well. \(\square\)

The nonzero diagonal condition on \(\mathbf{U}\) matters: it ensures that every back-substitution step can divide by \(u_{ii}\). If a diagonal entry is zero, this procedure does not guarantee a unique solution; the system may be inconsistent or may have more than one solution. In the examples below, all diagonal entries of \(\mathbf{U}\) are nonzero.

Worked Solves

Common Mistake: Using the Unpermuted Right-Hand Side

In the first worked example, applying forward substitution to \(\mathbf{L}\mathbf{y}=\mathbf{b}\) instead of \(\mathbf{L}\mathbf{y}=\mathbf{P}\mathbf{b}\) would give \(y_1=3\), \(y_2=2\), and \(2(3)-2+y_3=5\), so \(y_3=1\). Back substitution would then give \(x_3=1/2\), \(x_2=3/4\), and \(x_1=9/4\). This is a concrete but incorrect answer: the first entry of \(\mathbf{A}\mathbf{x}\) would be \(2(3/4)+1/2=2\), not the required \(3\).

The correction is to apply \(\mathbf{P}\) to the right-hand side before either triangular solve. The matrix \(\mathbf{P}\) in that example swaps the first two entries, so \(\mathbf{P}\mathbf{b}=(2,3,5)^{\top}\); using this vector gives the checked solution \((3/2,1/2,2)^{\top}\). A second common error is to apply \(\mathbf{P}^{\top}\) instead of \(\mathbf{P}\). For a single swap these matrices happen to be equal, but a longer permutation need not equal its transpose. The safe rule follows directly from the displayed identity: \(\mathbf{P}\mathbf{A}\mathbf{x}=\mathbf{P}\mathbf{b}\), so the right-hand side is \(\mathbf{P}\mathbf{b}\).

In practice, the factorization is reusable: for each new right-hand side, form its permuted vector, solve the two triangular systems, and check the resulting \(\mathbf{x}\) in the original equation when a direct verification is appropriate. The permutation changes the order of equations; it does not change the unknown vector or the meaning of the original system.