Kadomin

KADOMIN

Forward and Backward Substitution

To make practical use of the LU decomposition we’ve implemented, we need a way to solve the system for a given right-hand side vector \(\mathbf{b}\). Fortunately, this is straightforward. Consider the following:

\(
\begin{align}
\mathbf{Ax} &= \mathbf{b} \\
\mathbf{LUx} &= \mathbf{b} \\
\mathbf{Lc} &= \mathbf{b} \\
\end{align}
\)

with \(\mathbf{Ux} = \mathbf{c}\).

The solution we’re after is the vector \(\mathbf{x}\). To find it, we first solve \(\mathbf{Lc} = \mathbf{b}\). Once we have \(\mathbf{c}\), we then solve \(\mathbf{Ux} = \mathbf{c}\).

But what does it mean to “solve” these systems? Since both \(\mathbf{L}\) and \(\mathbf{U}\) are triangular matrices, these are triangular systems, which lets us use the earlier values to solve for the next ones. Let’s start with \(\mathbf{L}\) and some arbitrary \(\mathbf{b}\):

\(
\begin{align}
\mathbf{Lc} &= \mathbf{b} \\\\
\begin{bmatrix}
1 & 0 & 0\\
2 & 1 & 0\\
3 & 4 & 1\\
\end{bmatrix}
\begin{bmatrix}
c_1\\
c_2\\
c_3\\
\end{bmatrix}
&= \begin{bmatrix}
1\\
3\\
8\\
\end{bmatrix}
\end{align} \)

Notice that the first row gives us the equation \(c_1=1\). This variable is already isolated—it’s immediately solved. We can then substitute this value into the row below to solve for the next variable, since each one only depends on the variables above it. We continue this process until we’ve solved for every variable, ultimately arriving at the solution \(\mathbf{c}\).

\(
\begin{align}
2c_1 + c_2 &= 3 \\
c_2 &= 3 \,-\, 2c_1 \\
c_2 &= 1
\end{align}
\)

\(
\begin{align}
3c_1 + 4c_2 + c_3 &= 8 \\
c_3 &= 8 \,-\, 3c_1 \,-\, 4c_2\\
c_3 &= 1
\end{align}
\)

This gives us the solution for \(\mathbf{c}\) as:

\(
\mathbf{c} =
\begin{bmatrix}
1\\
1\\
1\\
\end{bmatrix}
\)

In general, in an upper triangular system, each variable depends only on the ones above it. For lower triangular systems, it’s the opposite—each variable depends only on the ones below. Once we’ve solved for \(\mathbf{c}\), we plug it into the equation \(\mathbf{Ux} = \mathbf{c}\) and solve for \(\mathbf{x}\).

\(
\begin{align}
\mathbf{Ux} &= \mathbf{c} \\\\
\begin{bmatrix}
2 & 4 & 9\\
0 & 3 & 4\\
0 & 0 & 1\\
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2\\
x_3\\
\end{bmatrix}
&= \begin{bmatrix}
1\\
1\\
1\\
\end{bmatrix}
\end{align} \)

Now it’s the last variable that is isolated:

\(x_3 = 1\)

\(
\begin{align}
3x_2 + 4x_3 &= 1 \\\\
x_2 &= \frac{1 \,-\, 4x_3}{3} \\\\
x_2 &= -1
\end{align}
\)

\(
\begin{align}
2x_1 + 4x_2 + 9x_3 &= 1 \\\\
x_1 &= \frac{1 \,-\, 4x_2 \,-\ 9x_3}{2} \\\\
x_1 &= -2
\end{align}
\)

This gives us the solution for \(\mathbf{x}\) as:

\(
\mathbf{x} =
\begin{bmatrix}
-2\\
-1\\
1\\
\end{bmatrix}
\)

Notice that, unlike in \(\mathbf{L}\), the diagonal entries in \(\mathbf{U}\) are no longer all 1. This means we need to divide the right-hand side of the equation by these diagonal values during substitution.

Solving \(\mathbf{Lc} = \mathbf{b}\) requires forward substitution (top to bottom), while solving \(\mathbf{Ux} = \mathbf{c}\) requires backward substitution (bottom to top). Now it’s time to implement both.

Implementation

The implementation is straightforward, but we need to clarify one important point. We will design the code to handle both forward and backward substitution in a general case, which means the diagonal elements can be any value (except zero of course), not just 1’s as seen in the special case of \(\mathbf{L}\) from LU decomposition.

We’ll begin with forward substitution (note that backward substitution follows the same principles but in reverse order, from bottom to top). The function will accept three parameters:

  1. A lower triangular matrix \(\mathbf{L}\) (which is not limited to the LU decomposition case, so it may contain non-1 values on the diagonal).
  2. The right-hand side vector \(\mathbf{b}\).
  3. A small tolerance value \(\epsilon\), which will serve as a threshold for numerical precision — any value smaller than epsilon will be treated as zero.
Python
def forward_substitution(L, b, epsilon: float = 1e-10):

Next, we need to validate that both L and b meet the required conditions. These are:

  1. \(\mathbf{L}\) must be a lower triangular matrix.
  2. \(\mathbf{L}\) must be a square matrix of size \(n \times n\).
  3. \(\mathbf{L}\) must be a 2D array.
  4. \(\mathbf{b}\) must be a 1D array.
  5. The length of \(\mathbf{b}\) must match the dimension of \(\mathbf{L}\) (i.e., \(\mathbf{L}\) must have size \(n\)).
  6. Neither \(\mathbf{L}\) nor \(\mathbf{b}\) can be empty.
  7. Both \(\mathbf{L}\) and \(\mathbf{b}\) must be convertible to NumPy arrays with numeric data types.

If any of these conditions are violated, the function will raise the appropriate error. It’s crucial to ensure the input is correct before proceeding with any further computations. These validation checks are implemented at the end of the full code provided.

Now, let’s dive into the logic. We’ll start with an empty vector c, which we will populate during the algorithm. First, we’ll solve for the first variable in the first row. This is easy since the only operation required is to divide the right-hand side by the factor corresponding to \(\mathbf{c}_1\), which is the (0,0) entry of \(\mathbf{L}\). Make sure to check that this entry is not zero, as that would cause problems.

Python
def forward_substitution(L, b, epsilon: float = 1e-10): 

    c = np.empty_like(b, dtype=np.result_type(L, b, float))

    # Solve first row
    d_0 = L[0,0]
    if abs(d_0) < epsilon:
        raise np.linalg.LinAlgError("Zero (or near-zero) pivot at row 0")
    c[0] = b[0] / d_0

Now that we have our starting point, we can loop through the remaining rows of the matrix. In each iteration, we grab the current row, trim off the part that contains zeros for efficiency, and then compute the dot product with the corresponding part of \(\mathbf{c}\). This is essentially a more efficient way of multiplying the previous variables by their respective coefficients in \(\mathbf{L}\). The dot product will give us a single value, which we subtract from the corresponding entry in the right-hand side vector \(\mathbf{b}\). Then, we divide the result by the coefficient in front of the variable we want to solve for.

Again, make sure that this coefficient is not zero, or close enough to zero that we consider it zero. Finally, we return the result. Now we have the complete code (including the checks for the properties of \(\mathbf{L}\) and \(\mathbf{b}\))

Python
def forward_substitution(L, b, epsilon: float = 1e-10): 

    try:
        L = np.asarray(L, dtype=float)
        b = np.asarray(b, dtype=float)
    except ValueError as e:
        raise TypeError(f"Input could not be converted to a numeric array. Got error: {e}")

    if L.ndim != 2:
        raise ValueError(f"L must be a 2D matrix, got {L.ndim}D input instead.")

    if L.size < 1:
        raise ValueError("Matrix L cannot be empty. Please provide a non-empty matrix.")

    if b.size < 1:
        raise ValueError("Vector b cannot be empty. Please provide a non-empty vector.")

    m, n = L.shape
    if m != n:
        raise ValueError(f"L must be square, got shape {L.shape!r}")

    if not np.allclose(L, np.tril(L), atol=0):
        raise ValueError("L needs to be lower triangular!")

    if b.ndim != 1 or b.size != n:
        raise ValueError(f"b must be a 1-D vector of length {n}, got shape {b.shape!r}")

    c = np.empty_like(b, dtype=np.result_type(L, b, float))

    # Solve first row
    d_0 = L[0,0]
    if abs(d_0) < epsilon:
        raise np.linalg.LinAlgError("Zero (or near-zero) pivot at row 0")
    c[0] = b[0] / d_0
    
    # Solve rows 1...n‑1
    for r in range(1, n):
        temp = L[r, :r] @ c[:r]
        div = L[r, r]
        if abs(div) < epsilon:
            raise np.linalg.LinAlgError(f"Zero (or near-zero) pivot at row {r}")
        c[r] = (b[r] - temp) / div

    return c

The code for backward substitution is almost identical. It starts at the bottom right corner of the matrix and works its way up. I’ll leave this as an exercise for the reader, as it’s mostly just a matter of adjusting the indices.

If you’re building something like a linear algebra module, as we are here, I’d recommend keeping forward and backward substitution in separate methods for clarity and modularity.

Example Run

We’ll run the example from the beginning of this page and check if we get the same results. The code for this looks something like this:

Python
L = np.array([[1,0,0], [1,2,0], [3,4,1]])
b = np.array([1,3,8])
c = forward_substitution(L, b)

print("solution c:\n", c)

The results are the same (thankfully):

Terminal
solution c:
 [1. 1. 1.]