Linear Systems of Equations

Linear Systems

Tridiagonal solver (Thomas algorithm); input has to be a tridiagonal square matrix and a solution vector. If the matrix is singular (i.e. determinant = 0), the algorithm will break (division by zero).

The algorithm is not generally stable (i.e. does not magnify small errors in the original matrix via, for example, rounding errors), but it is stable in special cases such as diagonally dominant matrices or symmetric positive definite matrices. In practice, one of the two is often the case.

First import the tools needed

import numpy as np
from math import factorial
import matplotlib.pyplot as plt
import copy

Thomas Algorithm

Create a function for solving tridiagonal matrix problems

def tridiag(B, z):
    s = np.shape(B)
    x = np.zeros(len(z))
    A = copy.deepcopy(B)
    y = copy.deepcopy(z)
    if s[0] == s[1]:
        n = s[0] - 1
        for i in range(1, n+1):
            w = A[i, i-1] / A[i-1, i-1]
            A[i, i] = A[i, i] - w * A[i-1, i]
            y[i] = y[i] - w * y[i-1]
        x[n] = y[n] / A[n, n]
        for i in range(n, 0, -1):
            j = i - 1
            x[j] = (y[j] - A[j, j+1] * x[j+1]) / A[j, j]
    else:
        print('Not a square matrix')
    return x

We first check if potential conditions for termination of the algorithm hold. This step might be as costly as or even more costly than the algorithm itself, which is why one wouldn’t do this in general. In many practical cases, there might be indications prior to using the algorithm that the conditions hold. Sometimes, one just simply applies the algorithm without checking at all and just checks if it fails. Below, we first apply the checks and then apply the algorithm.

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

Check for positive definiteness and symmetry

print('Matrix is pos. def. ', np.all(np.linalg.eigvals(A) > 0))
At = np.matrix.transpose(A)
print('Matrix is symmetric ', np.all(A==At))

Check for diagonal dominance: first get the absolute value of diagonal coefficients and then the sum of every row except the diagonal element

diag = np.diag(np.abs(A))
rowsum = np.sum(np.abs(A), axis=1) - diag
print('Matrix is diagonally dominant ', np.all(diag > rowsum))

Matrix is pos. def.  True
Matrix is symmetric  True
Matrix is diagonally dominant  True

Solve

x = tridiag(A, b)
print('Solution vector of the problem is ', x)

Solution vector of the problem is  [ 0.75 -0.5   2.75]

Cholesky decomposition

Cholesky decomposition of an input matrix $A$, which must be positive definite and symmetric. It returns the lower triangular matrix $L$. The upper one, $U$, is given by the transpose of $L$, i.e. $U=L^T$.

def cholesky(A):
    """
    Performs a Cholesky decomposition of A, which must be a symmetric and positive
    definite matrix. The function returns the lower triangular matrix, L.
    """

    # size
    n = len(A)

    # Create zero matrix for L
    L = np.zeros([n, n])

    # Perform the Cholesky decomposition
    for i in range(0, n):
        for k in range(0, i+1):

            # temporary variable
            tmp = sum(L[i, :] * L[k, :])

            # if statement for diagonal and nondiagonal elements
            if (i == k):
                L[i, k] = np.sqrt(A[i, i] - tmp)
            else:
                L[i, k] = 1.0 / L[k, k] * (A[i, k] - tmp)
    return L

We apply the Cholesky decomposition to the same matrix $A$ as before for the Thomas algorithm.

For solving an actual system of ${Ax=b}$ one now only needs to do substitution with the matrix $L$ and its transpose to find $x$.

L = cholesky(A)
print('Lower triangular matrix of Cholesky decomposition ')
print(L)

Lower triangular matrix of Cholesky decomposition
[[1.41421356 0.0 0.0 ]
[0.70710678 1.58113883 0.0]
[0.0 0.63245553 1.26491106]]

Check that $A = LL^T$

K = L.dot(L.transpose())
print('Reconstruction of original matrix by LL^T')
print(K)

Reconstruction of original matrix by LL^T
[[2. 1. 0.]
 [1. 3. 1.]
 [0. 1. 2.]]