Forward Substitution

Goal

Understand how to solve a lower triangular system by forward substitution. Master the $O(n^2)$ complexity, the simplification for unit lower triangular matrices, and the role in LU and Cholesky decompositions.

Prerequisites

  • Basic matrix and vector operations
  • Definition of lower and upper triangular matrices
Table of Contents

1. Definition and Principle

Forward substitution is an algorithm for solving a system of linear equations with a lower triangular coefficient matrix $L$:

$$L\mathbf{y} = \mathbf{b}$$

Since $L$ is lower triangular ($l_{ij} = 0$ for $i < j$), the first row gives a single-variable equation $l_{11} y_1 = b_1$, from which $y_1$ is immediately determined. Substituting the known $y_1$ into the second row yields $y_2$, and repeating this process through the last row determines all variables.

Written out explicitly:

$$\begin{pmatrix} l_{11} & 0 & \cdots & 0 \\ l_{21} & l_{22} & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}$$

A unique solution exists provided $l_{ii} \neq 0$ (i.e., $L$ is nonsingular).

Forward substitution: start from y₁ in row 1 and solve downward 2 0 0 1 3 0 −1 2 4 y₁ y₂ y₃ = 4 7 10 1 2 3 top to bottom 2·y₁ = 4 → y₁ = 2 1·y₁ + 3·y₂ = 7 −y₁ + 2·y₂ + 4·y₃ = 10 Every entry above the diagonal is 0, so each row adds one new unknown. ① Row 1: l₁₁y₁ = b₁ gives y₁ at once (green = solved). ② Substitute the known y₁ into the next row to find y₂. ③ In general yᵢ = (bᵢ − known terms) / lᵢᵢ, repeating downward. The exact mirror image of back substitution (solved upward).
Figure 1. Forward substitution. In the lower triangular system $L\mathbf{y}=\mathbf{b}$ the entries above the diagonal are $0$, so we start from $y_1$ in row 1 and solve from top to bottom, substituting each value into the rows below. Each $y_i = (b_i - \sum_{j<i} l_{ij}y_j)/l_{ii}$, the mirror image of back substitution.

2. Formula

$$y_1 = \dfrac{b_1}{l_{11}}$$
$$y_i = \dfrac{1}{l_{ii}} \left(b_i - \displaystyle\sum_{j=1}^{i-1} l_{ij} y_j\right), \quad i = 2, 3, \ldots, n$$

Each $y_i$ is obtained by subtracting the contributions of the already-known $y_1, \ldots, y_{i-1}$ from the right-hand side $b_i$, then dividing by the diagonal element $l_{ii}$. This proceeds in the opposite direction to back substitution, forming a symmetric pair.

3. Algorithm

function forwardSubstitution(L, b, n):
    y[1] = b[1] / L[1][1]
    for i = 2 to n:
        sum = 0
        for j = 1 to i-1:
            sum = sum + L[i][j] * y[j]
        y[i] = (b[i] - sum) / L[i][i]
    return y

4. Unit Lower Triangular Case

The $L$ factor obtained from LU decomposition is often a unit lower triangular matrix, meaning $l_{ii} = 1$. In this case, division is unnecessary and

$$y_i = b_i - \displaystyle\sum_{j=1}^{i-1} l_{ij} y_j$$

directly gives $y_i$. This saves $n$ division operations and avoids the rounding errors from division, leaving only those from the multiplications and additions/subtractions.

5. Computational Complexity

OperationGeneral $L$Unit Lower Triangular $L$
Divisions$n$$0$
Multiplications$\dfrac{n(n-1)}{2}$$\dfrac{n(n-1)}{2}$
Additions/Subtractions$\dfrac{n(n-1)}{2}$$\dfrac{n(n-1)}{2}$
Total$n^2$ ($O(n^2)$)$n^2 - n$ ($O(n^2)$)

6. Role in LU and Cholesky Decompositions

LU Decomposition

Given the LU decomposition $PA = LU$, the system $A\mathbf{x} = \mathbf{b}$ is solved in two steps:

  1. Solve $L\mathbf{y} = P\mathbf{b}$ by forward substitution.
  2. Solve $U\mathbf{x} = \mathbf{y}$ by back substitution.

Cholesky Decomposition

For a symmetric positive definite matrix $A = LL^T$ (Cholesky decomposition), the procedure is analogous:

  1. Solve $L\mathbf{y} = \mathbf{b}$ by forward substitution.
  2. Solve $L^T \mathbf{x} = \mathbf{y}$ by back substitution.

In the Cholesky decomposition, $L$ is a general lower triangular matrix ($l_{ii} \neq 1$), so divisions are required during forward substitution.

7. Numerical Stability

Forward substitution, like back substitution, is backward stable, and the growth of rounding errors is mild. For the unit lower triangular matrix $L$ obtained from LU decomposition with partial pivoting, $|l_{ij}| \leq 1$ is guaranteed, leading to particularly stable computation.

8. Worked Examples

Example 1: 3x3 Lower Triangular System

$\begin{pmatrix} 2 & 0 & 0 \\ 1 & 3 & 0 \\ -1 & 2 & 4 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 7 \\ 10 \end{pmatrix}$

StepComputationResult
$y_1$$4 / 2$$y_1 = 2$
$y_2$$(7 - 1 \times 2) / 3$$y_2 = 5/3$
$y_3$$(10 - (-1) \times 2 - 2 \times 5/3) / 4$$y_3 = 23/12$

The solution is $\mathbf{y} = (2,\; 5/3,\; 23/12)^T$.

Example 2: Unit Lower Triangular Matrix

$\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ -0.25 & 0.75 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 6 \\ 1 \\ 5 \end{pmatrix}$

StepComputationResult
$y_1$$6$$y_1 = 6$
$y_2$$1 - 0.5 \times 6$$y_2 = -2$
$y_3$$5 - (-0.25) \times 6 - 0.75 \times (-2)$$y_3 = 8$

Note that no divisions are needed.

9. Frequently Asked Questions

Q1. What is forward substitution?

It is an algorithm that solves a lower triangular system $L\mathbf{y} = \mathbf{b}$ by computing the variables sequentially from $y_1$ forward. It is a fundamental algorithm used in LU and Cholesky decompositions.

Q2. What is the computational complexity?

For an $n \times n$ lower triangular matrix, it is $O(n^2)$, the same order as back substitution.

Q3. What about unit lower triangular matrices?

Since $l_{ii} = 1$, division is unnecessary, and $y_i = b_i - \displaystyle\sum_{j=1}^{i-1} l_{ij} y_j$ gives the result directly.

10. References