Stationary Methods - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import matplotlib.pyplot as plt

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

The Intuition¶

Stationary methods are a class of iterative solvers that repeatedly apply the exact same mathematical operation at every iteration.

What does it mean for a vector $\mathbf{x}$ to satisfy $\mathbf{A}\mathbf{x}=\mathbf{b}$ ? It means that all the elements in $\mathbf{x}$ are perfectly consistent with each other. If we guess a solution where one element is inconsistent, we can algebraically solve for it based on the other elements!

The Algorithm:

Start with a random guess for $\mathbf{x}$ .
For each iterative step, identify a single element in $\mathbf{x}$ and force it to be consistent with the others.
Check the magnitude of the residual (how well the new guess solves the system).
Repeat until the residual falls below a specified tolerance.

Jacobi Iteration¶

Let’s formalize this by writing out the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ in index form:

\begin{aligned} \sum_{j} A_{ij} x_j &= b_i \\ \end{aligned}

(1)

If we extract the $i$ -th row, we can isolate and solve for the diagonal term $x_i$ :

\begin{aligned} A_{ii}x_i^{(k+1)} + \sum_{j \ne i} A_{ij} x_j^{(k)} &= b_i \\ x_i^{(k+1)} &= \frac{1}{A_{ii}}\left[b_i-\sum_{j \ne i} A_{ij} x_j^{(k)}\right] \end{aligned}

(2)

This equation explicitly updates $x_i^{(k+1)}$ to be consistent with all the other elements from the previous iteration $\mathbf{x}^{(k)}$ .

Matrix Formulation¶

Since we proceed row-by-row, we can vectorize this operation by pulling apart the coefficient matrix $\mathbf{A}$ into its diagonal $\mathbf{D}$ , strictly lower triangular $\mathbf{L}$ , and strictly upper triangular $\mathbf{U}$ components:

\mathbf{A} = \mathbf{L} + \mathbf{D} + \mathbf{U}

(3)

NB: This is not a matrix decomposition; we are literally just pulling it apart (blank elements are zeros):

\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} = \begin{bmatrix} & & \\ a_{21} & & \\ a_{22} & a_{23} & \end{bmatrix} + \begin{bmatrix} a_{11} & & \\ & a_{22} & \\ & & a_{33} \end{bmatrix} + \begin{bmatrix} & a_{12} & a_{13} \\ & & a_{23} \\ & & \end{bmatrix}

(4)

The algorithm becomes:

\begin{align} D x^{k+1} &= b-[L+U]x^k \\ x^{k+1} &= D^{-1} \big[b-[L+U]x^k\big] \end{align}

(5)

The Jacobi algorithm elegantly becomes:

\begin{aligned} \mathbf{D} \mathbf{x}^{(k+1)} &= \mathbf{b} - [\mathbf{L}+\mathbf{U}]\mathbf{x}^{(k)} \\ \mathbf{x}^{(k+1)} &= \mathbf{D}^{-1} \big[\mathbf{b} - [\mathbf{L}+\mathbf{U}]\mathbf{x}^{(k)}\big] \end{aligned}

(6)

(Note: Because the new vector $\mathbf{x}^{(k+1)}$ requires the entire previous vector $\mathbf{x}^{(k)}$ , we must store both in memory. However, because each element update is independent of the others in the same step, Jacobi iteration is incredibly easy to parallelize!)

Jacobi iterations are guaranteed to converge if the matrix $\mathbf{A}$ is diagonally dominant.

A = np.array([[4, 1], [1, 3]])
b = np.array([5, 6])
x = np.array([0., 0.])
iterations = [x.copy()]
tol = 1e-4

print('The true answer is:', np.linalg.solve(A, b))
print(f'Initial Guess: {x} | Residual: {np.linalg.norm(A@x-b):.4f}\n')

iteration_count = 0
while np.linalg.norm(A@x-b) > tol and iteration_count < 100:
    x = exe.jacobi_step(A, b, x)
    iterations.append(x.copy())
    iteration_count += 1
    print(f'Iteration {iteration_count:2d}: {x} | Residual: {np.linalg.norm(A@x-b):.4f}')

fig = exe.visualize_convergence_2d(A, b, iterations, surface_type='residual')
fig.show()

Gauss-Seidel¶

Gauss-Seidel operates on the exact same premise as Jacobi, but with one critical optimization: as soon as we calculate a new element $x_i^{(k+1)}$ , we immediately use it to calculate the next elements in the same iteration!

If we move downwards through the rows, the $i$ -th element is updated using the already-updated rows above it, and the old rows below it:

\begin{aligned} x_i^{(k+1)} &= \frac{1}{A_{ii}}\left[b_i-\sum_{j=1}^{i-1} A_{ij} x_j^{(k+1)} -\sum_{j=i+1}^n A_{ij} x_j^{(k)}\right] \end{aligned}

(7)

Matrix Formulation¶

Written in matrix form, we simply move the lower triangular matrix $\mathbf{L}$ to the left side since those elements are already updated:

\begin{aligned} [\mathbf{D} + \mathbf{L}] \mathbf{x}^{(k+1)} &= \mathbf{b} - \mathbf{U}\mathbf{x}^{(k)} \end{aligned}

(8)

Because we are modifying $\mathbf{x}$ in-place as we go, we only need to store a single vector in memory!

Gauss-Seidel is guaranteed to converge if the matrix is diagonally dominant or symmetric positive definite.

(Note: While Gauss-Seidel is vastly more memory-efficient and generally converges faster than Jacobi, the in-place updates mean each calculation depends on the previous one, making it extremely difficult to parallelize.)

A = np.array([[4, 1], [1, 3]])
b = np.array([5, 6])
tol = 1e-4

print('The true answer is:', np.linalg.solve(A, b))

# We can now use our generalized wrapper function
x_final, iterations = exe.iter_solve(A, b, exe.gauss_seidel_step, tol=tol, track_history=True)

fig = exe.visualize_convergence_2d(A, b, iterations, surface_type='residual')
fig.show()

Successive Over Relaxation (SOR)¶

The convergence of Gauss-Seidel can be accelerated through relaxation. Instead of simply accepting the new $x_i^{(k+1)}$ value, we take a weighted average of the newly calculated value and the old value using a relaxation parameter $\omega$ :

\begin{aligned} x_i &\leftarrow \frac{\omega}{A_{ii}}\left[b_i-\sum_{j \ne i} A_{ij} x_j \right] + (1-\omega) x_i \end{aligned}

(9)

If $\omega = 1$ : The method is exactly Gauss-Seidel.
If $\omega < 1$ : The method is under-relaxed (converges more smoothly but slower).
If $\omega > 1$ : The method is over-relaxed (convergence is artificially accelerated).

Choosing the Relaxation Parameter¶

SOR with a well-chosen $\omega$ will drastically outperform Gauss-Seidel. However, finding the optimal $\omega$ dynamically is non-trivial.

One heuristic approach involves monitoring the difference between successive iterations $\Delta \mathbf{x}^{(k)} = ||\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}||$ . After $p$ iterations, you can estimate the optimal relaxation factor as:

\omega_{opt} \approx \frac{2}{1+\sqrt{1+\left[ \frac{\Delta \mathbf{x}^{(k+p)}}{\Delta \mathbf{x}^{(k)}} \right]^p }}

(10)

A = np.array([[4, 1], [1, 3]])
b = np.array([5, 6])
tol = 1e-4

print('The true answer is:', np.linalg.solve(A, b))

# We can now use our generalized wrapper function
x_final, iterations = exe.iter_solve(A, b, exe.sor_step, tol=tol, track_history=True)

fig = exe.visualize_convergence_2d(A, b, iterations, surface_type='residual')
fig.show()

Performance Benchmarking¶

How do these stationary methods stack up against a highly optimized direct solver like scipy.linalg.solve?

Let’s generate large diagonally-dominant matrices and time the execution of our iterative loops versus the direct solver. Try changing the ‘n’ field in the code below!

n = 4  # Size of the matrix
A = exe.create_diagonally_dominant_matrix(n)
b = np.random.rand(n)  # Create a random right-hand side vector

print('LU')
%timeit -n1 -r1 np.linalg.solve(A, b)

print('\nJacobi')
%timeit -n1 -r1 exe.iter_solve(A, b, exe.jacobi_step)
print('\nGauss-Seidel')
%timeit -n1 -r1 exe.iter_solve(A, b, exe.gauss_seidel_step)
print('\nSOR')
%timeit -n1 -r1 exe.iter_solve(A, b, exe.sor_step)