Solving Ax=b with the LU decomposition

Solving $Ax=b$ problems is one of the main uses of the LU decomposition. To solve $Ax=b$ with the LU decomp, we substitute the factorization $A=LU$ into $Ax=b$ to get

\begin{equation*} LU x = b \end{equation*}

We solve this by first letting $y = Ux$ (where $y$ and $x$ are both unknown vectors) and solving

\begin{equation*} Ly = b \end{equation*}

for $y$. Now $y$ is known, and we can find $x$ by solving

\begin{equation*} Ux = y \end{equation*}

for $x$. It turns out that solving $Ly=b$ and $Ux=b$ is much easier than the original $Ax=b$ problem due to the special lower-triangular and upper-triangular structure of $L$ and $U$.

Forward substitution: the 3 x 3 case.

To solve $Ly=b$, where $L$ is a given lower-triangular matrix, $b$ is a known right-hand-side vector, and $y$ is unknown, we use an algorithm known as forward substitution.

We'll illustrate with a 3 x 3 system $Ly=b$, where $L$ is lower-triangular with 1's on the diagonal.

\begin{align*} \left[\begin{array}{ccc} 1 & 0 & 0 \\ \ell_{21} & 1 & 0 \\ \ell_{31} & \ell_{32} & 1 \end{array}\right] \left[\begin{array}{c} y_1 \\ y_2 \\ y_3 \end{array} \right] = \left[ \begin{array}{c} b_1 \\ b_2 \\ b_3 \end{array} \right] \end{align*}

We solve this system in three steps.

Step 1 The top row of the 3 x 3 system is the equation

\begin{align*} y_1 = b_1 \end{align*}

Thus we find the first unknown $y_1$ simply given by setting it to the known value $b_1$.

Step 2 The second row of the system is the equation

\begin{align*} \ell_{21} y_1 + y_2 = b_2 \end{align*}

Solving for $y_2$ gives

\begin{align*} y_2 = b_2 - \ell_{21} y_1 \end{align*}

Since we know the value of $y_1$ from the last step, we determine the value of $y_2$ by calculating the right-hand-side of this equation.

Step 3 The third row of the system is

\begin{align*} \ell_{31} y_1 + + \ell_{32} y_2 + y_3 = b_3 \end{align*}

Solving for $y_3$ gives

\begin{align*} y_3 = b_3 - \ell_{31} y_1 - \ell_{32} y_2 \end{align*}

At this stage we know $y_1$ and $y_2$, so we determine the value of $y_3$ by evaluating the expression on the right-hand side.

Forward substitution, general case

The general $m \times m$ of $Ly=b$ is a straightforward generalization of the $3 \times 3$ case shown above. The first three equations are solved exactly as shown above. The fourth equation is

\begin{align*} \ell_{41} y_1 + \ell_{42} y_2 + \ell_{43} y_3 + y_4 = b_4 \end{align*}

which when solved for $y_4$ becomes

\begin{align*} y_4 = b_4 - \ell_{41} y_1 - \ell_{42} y_2 - \ell_{43} y_3 \end{align*}

Generalizing, we get

\begin{align*} y_i = b_i - \sum_{j=1}^{i-1} \ell_{ij} y_j \end{align*}

Forward substitution, pseudo code

Given an $m \times m$ lower-triangular matrix $L$ and a right-hand-side vector $b$, the above formula for the solution $y$ to the system $Ly=b$ can be calculated via the following pseudo code

# for i = 1 to m # yᵢ = bᵢ # for j = 1 to i-1 # yᵢ = yᵢ - ℓᵢⱼ yⱼ # end # end

Your next homework assignment, HW3, will call for you to

  • implement a forward substitution algorithm in Julia
  • derive the formula for backsubstitution to solve $Ux=y$
  • implement a backward substitution algorithm in Julia
  • solve a number of $Ax=b$ problems using your LU decomp and forward & backward substitution algorithms
In [ ]: