Suppose that $A$ is an upper triangular matrix $$ A = \begin{bmatrix} a_{1,1} & \cdots & a_{1, n} \\ \vdots & \ddots & \vdots \\ 0 & \cdots & a_{n, n} \end{bmatrix}, $$ where $a_{i, j} = 0$ for all $i > j$. We want to solve $$Ax = b,$$ where $b$ is a vector in $\mathbb{R}^n$. The back substitution algorithm gives the solution, according to the formula $$x_i = \frac{b_i - \sum_{j = i+1}^n a_{i, j}\, x_j}{a_{i,i}}.$$

In?[1]:

import numpy as np

In?[2]:

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

In?[3]:

print(A)

[[1 2 3]
 [0 3 4]
 [0 0 5]]

In?[4]:

print(b)

[1, 1, 1]

In?[5]:

def back_substitution(A, b):

    n = np.shape(A)[0]
    x = np.zeros(n)

    for i in range(n - 1, -1, -1):
        temp_sum = 0
        for j in range(i + 1, n):
            temp_sum += A[i, j] * x[j]
        x[i] = (b[i] - temp_sum) / A[i, i]

    return x

In?[6]:

x = back_substitution(A, b)
print(x)

[0.26666667 0.06666667 0.2       ]

In?[8]:

x_test = np.linalg.solve(A, b)
print(x_test)

[0.26666667 0.06666667 0.2       ]

In python, it is always smart to look for vectorized operations. This reduces the runtime of algorithms.

In?[9]:

def back_substitution(A, b):
    n = np.shape(A)[0]
    x = np.zeros(n)

    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

    return x

Some issues which can arise:

Singular matrix
Numerical instability

In?[28]:

N = 50
# Matrix of uniformly distributed numbers
A = np.random.rand(N, N) + 1

for j in range(0, N):
    for i in range(j+1, N):
        A[i, j] = 0

b = np.random.rand(N) + 1

The condition number hints to what happens for large matrices here.

In?[29]:

print(np.linalg.cond(A))

769.9036691798275

In?[30]:

x = back_substitution(A, b)

In?[31]:

x_test = np.linalg.solve(A, b)

In?[32]:

print(np.linalg.norm(x - x_test, 2))

4.334651545230797e-15

In?[?]: