Nearly any engineer learned some day to solve by hand linear systems of equations in its academic background. The Gauss method is usually the main algorithm taught. While it might seem basic, the process is an omnipresent piece of linear algebra in nearly every engineering of science field. I enjoyed diving through once again with new insights.

This article will go through the matrix formulation of Gauss elimination, the (P)LU factorization and how it is used every day to solve linear systems of equations.

Gauss elimination

The Gauss elimination algorithm aims at transforming a system of linear equations in an "easier to work with" version. This version is called row echelon form and is basically a setup where the answer can be found by straightforward substitutions.

Let me introduce an example of linear system of equation which has a unique solution:

$$ \begin{aligned} 1x + 1y + 1z + 0w &= b_0\cr 0x + 3y + 1z + 2w &= b_1 \cr 2x + 3y + 1z + 0w &= b_2 \cr 1x + 0y + 2z + 3w &= b_3 \cr \end{aligned} $$

If you remember correctly, the first approach one would take to solve such a system is to manipulate rows via multiplications, subtractions and permutations, ideally getting to the row echelon form (REF):

$$ \begin{aligned} 2x + 3y + 1z + 0w &= b'_0\cr 3y + 1z + 2w &= b'_1 \cr 2z + 4w &= b'_2 \cr -1w &= b'_3 \cr \end{aligned} $$

From here, the succession of substitutions starting from $w$ leads to the final system solution.

While it's perfectly fine to think about this problem as a system of independent equations, the matrix formulation represents an elegant way to frame the computations involved.

Gauss elimination with matrices

First, we'll place the $x$, $y$, $z$ and $w$ coefficients in a matrix $\mathbf{A}$.

$$ \mathbf{Ax}= \begin{bmatrix} 1 & 1 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 2 & 3 & 1 & 5 \cr 1 & 0 & 2 & 3 \cr \end{bmatrix}\begin{bmatrix} x \cr y \cr z \cr w \cr \end{bmatrix}=\mathbf{b} $$

Using the matrix versions of row permutations, subtraction and multiplication, one can obtain the REF matrix $\mathbf{U}$. I call it $\mathbf{U}$ here because in case of non-singular matrices, the REF form is an Upper triangular matrix. Let's review the basic matrix operations.

Elementary permutation

The term elementary induces that only one permutation is performed while applying the matrix $\mathbf{P}$. This constraint has for consequence the following property of $\mathbf{P}$: $$ \mathbf{P}^{-1} = \mathbf{P} $$ In order to build a permutation matrix that permutes the rows $i$ and $j$, swap the rows $i$ and $j$ of the identity matrix $\mathbf{I}$.

$$ \mathbf{I}= \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr \end{bmatrix} \rightarrow \begin{bmatrix} 0 & 0 & 1 & 0 \cr 0 & 1 & 0 & 0 \cr 1 & 0 & 0 & 0 \cr 0 & 0 & 0 & 1 \cr \end{bmatrix}=\mathbf{P}_{13} $$

A python code for the generation of such a matrix could be:

def permutation(i, j, N):
    P = np.eye(N)
    i_row, j_row = np.zeros(N), np.zeros(N)
    i_row[j-1], j_row[i-1] = 1, 1
    P[i-1], P[j-1] = i_row, j_row
    return P

Multiply & Substract

One way to view a matrix product is that each row of the left matrix represents a linear combination of the rows of the right matrix. Using this knowledge, building a matrix that replaces the second row by the sum of the row $1$ and $2$ can be written like:

$$ \mathbf{E}_{12}(1)= \begin{bmatrix} 1 & 0 & 0 & 0 \cr 1 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr \end{bmatrix} $$

If I want to subtract half of it, I simply change the coefficient to $-0.5$.

$$ \mathbf{E}_{12}(-0.5) \begin{bmatrix} 1 & 0 & 0 & 0 \cr -0.5 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr \end{bmatrix} $$

We will notate this matrix $\mathbf{E}_{ij}(c)$ which is like the identity matrix $\mathbf{I}$ but with the element at position $(i,j)$ equal to $c$. Because we are working top-down on the simplification of the matrix, only matrices $\mathbf{E}_{ij}$ with $j < i$ (under the diagonal) will be allowed.

Again, the code could be:

def elimination(i, j, c, N):
    assert j-1 < i-1, "Elimination only on lower triangular part."
    E = np.eye(N)
    E[i-1][j-1] = c
    return E

Again, a nice observation is that the inverse of such lower triangular can be obtain by changing the sign of the lower triangular part, this will be handy for factorization.

$$\mathbf{E}^{-1}_{ij}(c)=\mathbf{E}_{ij}(-c)$$

Factor $\mathbf{A}=\mathbf{LU}$

While it's possible to apply the Gauss elimination for each $\mathbf{b}$ you are trying to solve, this is computationally inefficient if $\mathbf{A}$ doesn't change (you have to compute the elimination each time). One way to record the history of manipulations of the Gauss elimination is to factor $\mathbf{A}$ in two matrices $\mathbf{L}$ and $\mathbf{U}$.

• $\mathbf{U}$ is the Gauss elimination result as an upper triangular matrix.
• $\mathbf{L}$ is the lower triangular matrix with coefficients of row multiplies and subtracts.

We will see in next section how $\mathbf{L}$ and $\mathbf{U}$ can be used to solve this equation multiple times without going through elimination: $$\mathbf{PAx}=\mathbf{LUx}=b$$

In a generic point of view, factorizing $\mathbf{A}=\mathbf{LU}$ (without thinking about permutations for now) consists in first applying the row elimination matrices $\mathbf{E}_{ij}(c_{ij})$ on $\mathbf{A}$ to get $\mathbf{U}$. $$\mathbf{E}_{43}\mathbf{E}_{33}\ldots\mathbf{E}_{41}\mathbf{E}_{31}\mathbf{E}_{21}\mathbf{A}=\mathbf{U}$$

Then, we can invert the left side: $$ \begin{aligned} (\mathbf{E}_{43}\mathbf{E}_{33}\ldots\mathbf{E}_{41}\mathbf{E}_{31}\mathbf{E}_{21})^{-1}\mathbf{A}&=\mathbf{U} \cr \mathbf{A}&=\mathbf{E}^{-1}_{21}\mathbf{E}^{-1}_{31}\mathbf{E}^{-1}_{41}\ldots\mathbf{E}^{-1}_{33}\mathbf{E}^{-1}_{43}\mathbf{U} \cr \mathbf{A}&=\mathbf{LU} \end{aligned} $$

With $\mathbf{L}=\mathbf{E}^{-1}_{21}\mathbf{E}^{-1}_{31}\mathbf{E}^{-1}_{41}\ldots\mathbf{E}^{-1}_{33}\mathbf{E}^{-1}_{43}$ and using the following property of $\mathbf{E}^{-1}_{ij}(c)=\mathbf{E}_{ij}(-c),$ the inverses of $\mathbf{E}_{ij}$ are easy to compute.

Solve $\mathbf{LUx}=\mathbf{b}$

Once $\mathbf{L}$ and $\mathbf{U}$ (and $\mathbf{P}$) are computed by one pass of Gauss elimination, the following system is yet to solve:

$$\mathbf{Ax}=\mathbf{LUx}=\mathbf{b}$$

This can actually be worked on as two separate systems to solve: $$ \left\{ \begin{aligned} \mathbf{Ly}&=\mathbf{b}\cr \mathbf{Ux}&=\mathbf{y} \end{aligned} \right. $$

In order to get the solution, first solve $y = b \ \backslash \ L$ with forward substitution, then $x = y \ \backslash \ U$ with backward substitution.

What about permutations?

The Gauss method may need row permutations in order to place pivots in the right places. Also, partial pivoting, the act of taking the greatest absolute value of a column as the pivot for numerical stability also heavily permutes the matrix under elimination.

Our objective would be here to factorize the permutations too, for example in this form:

$$\mathbf{PA}=\mathbf{LU}$$

Where we would get the global permutation matrix $\mathbf{P}$ from the first elimination, then solving for $\mathbf{PA}$ multiple times using the $\mathbf{LU}$ decomposition.

Under this knowledge, we can reframe the factorization with permutations before eliminations (if no permutation is needed for a given column under elimination, consider $\mathbf{P}_{ij}$ being $\mathbf{I}$).

$$\mathbf{E}_{43}\mathbf{E}_{33}\mathbf{P}_3\ldots\mathbf{E}_{41}\mathbf{E}_{31}\mathbf{E}_{21}\mathbf{P}_1\mathbf{A}=\mathbf{U}$$

This looks a bit messy, combining each $\mathbf{E}_{ij}$ for a given column $j$ in a matrix $\mathbf{L}_{j}$, we rewrite the same equation as:

$$\mathbf{L}_3\mathbf{P}_3\mathbf{L}_2\mathbf{P}_2\mathbf{L}_1\mathbf{P}_1\mathbf{A}=\mathbf{U}$$

Now a trick is needed in order to get $\mathbf{L}_i$ and $\mathbf{P}_i$ stacked. Starting from what we want, we would like to get at some point:

$$\mathbf{L'}_3\mathbf{L'}_2\mathbf{L'}_1\mathbf{P}_3\mathbf{P}_2\mathbf{P}_1\mathbf{A}=\mathbf{U}$$

Which is equal to previous equation by setting: $$ \left\{ \begin{aligned} \mathbf{L'}_1&=\mathbf{P}_3\mathbf{P}_2\mathbf{L}_1(\mathbf{P}_3\mathbf{P}_2)^{{-1}}\cr \mathbf{L'}_2&=\mathbf{P}_3\mathbf{L}_2\mathbf{P}_3^{{-1}}\cr \mathbf{L'}_3&=\mathbf{L}_3\cr \end{aligned} \right. $$

Or, in a generic manner: $$\mathbf{L'}_i=\mathbf{P}_{n-1}\ldots\mathbf{P}_{i+1}\mathbf{L}_i(\mathbf{P}_{n-1}\ldots\mathbf{P}_{i+1})^{-1}$$

From there, we are able to write the following:

$$ \begin{aligned} \mathbf{P}_3\mathbf{P}_2\mathbf{P}_1\mathbf{A}&=\mathbf{L'}^{-1}_1\mathbf{L'}^{-1}_2\mathbf{L'}^{-1}_3\mathbf{U} \cr \mathbf{P}\mathbf{A}&=\mathbf{L}\mathbf{U} \cr \end{aligned} $$

We now know how to factorize any non-singular matrix using the $\mathbf{PLU}$ decomposition. It is now time for the step-by-step example.

Step-by-step example

I will take the same matrix $\mathbf{A}$ that was introduced at the beginning of the post:

$$ \mathbf{A}= \begin{bmatrix} 1 & 1 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 2 & 3 & 1 & 5 \cr 1 & 0 & 2 & 3 \cr \end{bmatrix} $$

I will work out factorization column in order to see how the Gauss elimination looks like. Then I will compute matrices $\mathbf{P}$ and $\mathbf{L}$.

Let's start!

First column

Permuting rows $1$ and $3$ in order to get the $2$ as pivot, then eliminating rows $2$, $3$ and $4$ from the first column. $$ \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr \scriptstyle -1/2 & 0 & 1 & 0 \cr \scriptstyle -1/2 & 0 & 0 & 1 \cr \end{bmatrix}}_{\mathbf{L}_1} \underbrace{ \begin{bmatrix} 0 & 0 & 1 & 0 \cr 0 & 1 & 0 & 0 \cr 1 & 0 & 0 & 0 \cr 0 & 0 & 0 & 1 \cr \end{bmatrix}}_{\mathbf{P}_1} \underbrace{ \begin{bmatrix} 1 & 1 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 2 & 3 & 1 & 5 \cr 1 & 0 & 2 & 3 \cr \end{bmatrix}}_{\mathbf{A}} = \underbrace{ \begin{bmatrix} 2 & 3 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 0 & \scriptstyle -1/2 & \scriptstyle 1/2 & 0 \cr 0 & \scriptstyle -3/2 & \scriptstyle 3/2 & 3 \cr \end{bmatrix}}_{\mathbf{U}_1} $$

Second column

No permutation and eliminating rows $3$ and $4$ from the second column. $$ \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & \scriptstyle 1/6 & 1 & 0 \cr 0 & \scriptstyle 1/2 & 0 & 1 \cr \end{bmatrix}}_{\mathbf{L}_2} \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr \end{bmatrix}}_{\mathbf{P}_2} \underbrace{ \begin{bmatrix} 2 & 3 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 0 & -\scriptstyle 1/2 & \scriptstyle 1/2 & 0 \cr 0 & -\scriptstyle 3/2 & \scriptstyle 3/2 & 3 \cr \end{bmatrix}}_{\mathbf{U}_1} = \underbrace{ \begin{bmatrix} 2 & 3 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 0 & 0 & \scriptstyle 2/3 & \scriptstyle 1/3 \cr 0 & 0 & 2 & 4 \cr \end{bmatrix}}_{\mathbf{U}_2} $$

Third column

Finally permuting rows $3$ and $4$ and proceeding to elimination on third column in order to get the final $\mathbf{U}$. $$ \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & \scriptstyle -1/3 & 1 \cr \end{bmatrix}}_{\mathbf{L}_3} \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 0 & 1 \cr 0 & 0 & 1 & 0 \cr \end{bmatrix}}_{\mathbf{P}_3} \underbrace{ \begin{bmatrix} 2 & 3 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 0 & 0 & \scriptstyle 2/3 & \scriptstyle 1/3 \cr 0 & 0 & 2 & 4 \cr \end{bmatrix}}_{\mathbf{U}_2} = \underbrace{ \begin{bmatrix} 2 & 3 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 0 & 0 & 2 & 4 \cr 0 & 0 & 0 & -1 \cr \end{bmatrix}}_{\mathbf{U}} $$

Computing P

We previously set the expression for the global permutation matrix $\mathbf{P}$.

$$ \mathbf{P} = \mathbf{P}_3\mathbf{P}_2\mathbf{P}_1 = \begin{bmatrix} 0 & 0 & 1 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 0 & 1 \cr 1 & 0 & 0 & 0 \cr \end{bmatrix} $$

Computing L

Then, following the formula for $\mathbf{L'}_i$ derived before:

$$\mathbf{L'}_i=\mathbf{P}_{n-1}\ldots\mathbf{P}_{i+1}\mathbf{L}_i(\mathbf{P}_{n-1}\ldots\mathbf{P}_{i+1})^{-1},$$

with

$$ \mathbf{L}=\mathbf{L'}^{-1}_1\mathbf{L'}^{-1}_2\mathbf{L'}^{-1}_3 $$

We find the $\mathbf{L}$ matrix.

$$ \mathbf{L}= \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr \scriptstyle 1/2 & -\scriptstyle 1/2 & 1 & 0 \cr \scriptstyle 1/2 & -\scriptstyle 1/6 & \scriptstyle 1/3 & 1 \cr \end{bmatrix} $$

The final factorization

And that's it, we get the final factorization.

$$ \underbrace{ \begin{bmatrix} 0 & 0 & 1 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 0 & 1 \cr 1 & 0 & 0 & 0 \cr \end{bmatrix}}_{\mathbf{P}} \underbrace{ \begin{bmatrix} 1 & 1 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 2 & 3 & 1 & 5 \cr 1 & 0 & 2 & 3 \cr \end{bmatrix}}_{\mathbf{A}}= \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr \scriptstyle 1/2 & -\scriptstyle 1/2 & 1 & 0 \cr \scriptstyle 1/2 & -\scriptstyle 1/6 & \scriptstyle 1/3 & 1 \cr \end{bmatrix}}_{\mathbf{L}} \underbrace{ \begin{bmatrix} 2 & 3 & 1 & 0 \cr 0 & 3 & 1 & 2 \cr 0 & 0 & 2 & 4 \cr 0 & 0 & 0 & -1 \cr \end{bmatrix}}_{\mathbf{U}} $$

Using $\mathbf{P}$, $\mathbf{L}$ and $\mathbf{U}$ one can now solve linear systems of equation using forward and back substitution.

I hope this post gave you a clearer insight on how the $\mathbf{PLU}$ decomposition happens under the hood everyday :)

Python code

Below is the Python code I wrote for educational purposes (do not use, not optimized) that computes the $\mathbf{PLU}$ factorization with partial pivoting.

def plu_factor(A):
    # Gauss elimination
    N, U = A.shape[0], A.copy()
    Pi, Li = [], []
    for col in range(N-1):
        Ei = []
        # Partial pivoting
        max_abs_indice = np.argmax(np.abs(U[col:, col])) + col + 1 
        P = permutation(col + 1, max_abs_indice, N)
        U = P @ U
        Pi.append(P)
        for row in range(col+1, N):
            i, j, c = row+1, col+1, -U[row,col]/U[col,col]
            E = elimination(i, j, c, N)
            U = E @ U
            Ei.append(E)
        Li.append(np.sum(np.tril(np.stack(Ei), k=-1), axis = 0) + np.eye(N))
    
    # Permutation matrix
    P = np.linalg.multi_dot(list(reversed(Pi)))
    
    # Lower triangular matrix
    Li_prime = np.zeros((len(Li), N, N))
    for i in range(len(Li)):
        P_seq = [np.eye(N)] if i+1 >= len(Pi) else Pi[i+1:]
        Li_prime[i] = mprod(list(reversed(P_seq))) @ Li[i] @ mprod(P_seq)
    L = np.sum(-np.tril(Li_prime, k = -1), axis = 0) + np.eye(N)
    
    return P, L, U

With its utility functions:

def permutation(i, j, N):
    P = np.eye(N)
    i_row, j_row = np.zeros(N), np.zeros(N)
    i_row[j-1], j_row[i-1] = 1, 1
    P[i-1], P[j-1] = i_row, j_row
    return P

def elimination(i, j, c, N):
    E = np.eye(N)
    E[i-1][j-1] = c
    return E

def mprod(mat):
        return mat[0] if len(mat) == 1 else np.linalg.multi_dot(mat)

Get the results of this post back by calling the function on $\mathbf{A}$:

A = np.array(
    [[1, 1, 1, 0], 
     [0, 3, 1, 2], 
     [2, 3, 1, 0], 
     [1, 0, 2, 3]]
)

P, L, U = plu_factor(A)