Tutorial 03 - Linear algebra¶

Basic linear algebra operations, direct methods for solving linear equation systems, forward and backward substitution, Gaussian elimination, LU decomposition, Thomas algorithm.

import numpy as np
import scipy.linalg as la

Basic operations¶

Scalar product¶

The scalar product of two vectors \( x \in \mathbb{R}^{n} \) and \( y \in \mathbb{R}^{n} \) is defined as

\[ x \cdot y = \sum_{i=0}^{n-1} x_i y_i. \]

Exercise 03.1: Write a function that calculates the scalar product of two vectors (do not use the numpy.dot function nor the @ operator).

def scalar_product(x, y):
    """
    Calculates scalar product of two vectors.
    Args:
        x (array_like): Vector of size n
        y (array_like): Vector of size n
    Returns:
        numpy.float: scalar product of x and y
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the numpy.dot function.

x = np.random.rand(3)
y = np.random.rand(3)
try:
    np.testing.assert_array_almost_equal(scalar_product(x, y), np.dot(x, y), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Matrix-vector product¶

The product of matrix \( \mathbb{A} \in \mathbb{R}^{m \ \times \ n} \) and vector \( x \in \mathbb{R}^{n} \) is defined as

\[ (\mathbb{A} \cdot x)_i = \sum_{j = 0}^{n-1} a_{ij} x_j, \quad i = 0, 1, \dots, m-1. \]

Exercise 03.2: Write a function that calculates the product of a matrix and a vector (do not use the numpy.dot function nor the @ operator).

def matrix_vector_product(A, x):
    """
    Calculates matrix-vector product.
    Args:
        A (array_like): A m-by-n matrix
        x (array_like): Vector of size n
    Returns:
        numpy.ndarray: Matrix-vector product
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the numpy.dot function.

A = np.random.rand(3, 3)
x = np.random.rand(3)
try:
    np.testing.assert_array_almost_equal(matrix_vector_product(A, x), np.dot(A, x), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Matrix-matrix product¶

The product of two matrices \( \mathbb{A} \in \mathbb{R}^{m \ \times \ n} \) and \( \mathbb{B} \in \mathbb{R}^{n \ \times \ p} \) is defined as

\[ (\mathbb{A} \cdot \mathbb{B})_{ij} = \sum_{k = 0}^{n-1} a_{ik} b_{kj}, \qquad i = 0, 1, \dots, m-1, \quad j = 0, 1, \dots, p-1. \]

Exercise 03.3: Write a function that calculates the product of two matrices (do not use the numpy.dot function nor the @ operator).

def matrix_matrix_product(A, B):
    """
    Calculates matrix-matrix product.
    Args:
        A (array_like): A m-by-n matrix
        B (array_like): A n-by-p matrix
    Returns:
        numpy.ndarray: Matrix-matrix product
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the numpy.dot function.

A = np.random.rand(3, 3)
B = np.random.rand(3, 3)
try:
    np.testing.assert_array_almost_equal(matrix_matrix_product(A, B), np.dot(A, B), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Direct methods for linear algebraic equations¶

Forward substitution¶

A system of linear algebraic equations \( \mathbb{A} x = b \) with lower triangular coefficient matrix \( \mathbb{A} \in \mathbb{R}^{n \ \times \ n} \) can be solved by the forward substitution algorithm,

\[ x_i := \frac{1}{a_{ii}} \left( b_i - \sum_{j = 0}^{i - 1} a_{ij} x_j \right), \qquad i = 0, 1, \dots, n-1. \]

Exercise 03.4: Implement the forward substitution algorithm.

def forward_substitution(A, b):
    """
    Solves a system of linear equations with lower triangular matrix.
    Args:
        A (array_like): A n-by-n lower triangular matrix
        b (array_like): RHS vector of size n
    Returns:
        numpy.ndarray: Vector of solution
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the scipy.linalg.solve function.

A = np.tril(np.random.rand(3, 3)) # create a lower triangular matrix from a random matrix
b = np.random.rand(3)
try:
    np.testing.assert_array_almost_equal(forward_substitution(A, b), la.solve(A, b), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Backward substitution¶

A system of linear algebraic equations \( \mathbb{A} x = b \) with upper triangular coefficient matrix \( \mathbb{A} \in \mathbb{R}^{n \ \times \ n} \) can be solved by the backward substitution algorithm,

\[ x_i := \frac{1}{a_{ii}} \left( b_i - \sum_{j = i + 1}^{n-1} a_{ij} x_j \right), \qquad i = n-1, n-2, \dots, 0. \]

Exercise 03.5: Implement the backward substitution algorithm.

def backward_substitution(A, b):
    """
    Solves a system of linear equation with upper triangular matrix.
    Args:
        A (array_like): A n-by-n upper triangular matrix
        b (array_like): RHS vector of size n
    Returns:
        numpy.ndarray: Vector of solution
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the scipy.linalg.solve function.

A = np.triu(np.random.rand(3, 3)) # create an upper triangular matrix from a random matrix
b = np.random.rand(3)
try:
    np.testing.assert_array_almost_equal(backward_substitution(A, b), la.solve(A, b), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Gaussian elimination¶

A matrix \( \mathbb{A} \in \mathbb{R}^{n \ \times \ n} \) can be transformed into an upper triangular form using the Gaussian elimination algorithm,

\[ a_{jk} := a_{jk} - \frac{a_{ji}}{a_{ii}} a_{ik}, \qquad i = 0, 1, \dots, n-1, \quad j = i + 1, i + 2, \dots, n-1, \quad k = 0, 1, \dots, n-1. \]

Exercise 03.6: Implement the Gaussian elimination algorithm.

def gaussian_elimination(A):
    """
    Transforms given matrix into an upper triangular form using the Gaussian elimination algorithm.
    Args:
        A (array_like): A n-by-n matrix
    Returns:
        numpy.ndarray: Upper triangular matrix
    """

    # add your code here

Use the function above to transform the following matrix \( \mathbb{A} \) into an upper triangular form,

\[\begin{split} \mathbb{A} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}. \end{split}\]

# add your code here

Warning: The Gaussian elimination algorithm itself is unable to transform the matrix above into an upper triangular form because it leads to division by zero.

Pivoting: The pivot or pivot element is the element of a matrix, or an array, which is selected by an algorithm to do certain calculations. In the case of Gaussian elimination, the algorithm requires that pivot elements are distinct from zero. Therefore, interchanging rows or columns in the case of a zero pivot element is necessary. Furthermore, it is generally desirable to choose a pivot element with large absolute value, which dramatically reduces the round-off error.

In partial pivoting, the algorithm selects the entry with largest absolute value from the column of the matrix that is currently being considered as the pivot element. In complete pivoting, the algorithm interchanges both rows and columns in order to use the largest (by absolute value) element in the matrix as the pivot. For Gaussian elination, partial pivoting is usually sufficient.

Exercise 03.7: Implement pivoting into the Gaussian elimination algorithm.

def gaussian_elimination_with_pivoting(A):
    """
    Transforms given matrix into an upper triangular form using the Gaussian elimination algorithm with pivoting.
    Args:
        A (array_like): A n-by-n matrix
    Returns:
        numpy.ndarray: Upper triangular matrix
    """
    
    # add your code here

Use the function above to transform the following matrix \( \mathbb{A} \) into an upper triangular form,

\[\begin{split} \mathbb{A} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}. \end{split}\]

# add your code here

Exercise 03.8: Pass the right-hand side vector into the Gaussian elimination algorithm and perform the identical operations on it.

def gaussian_elimination_with_pivoting_vector(A, b):
    """
    Transforms given matrix into an upper triangular form using the Gaussian elimination algorithm with pivoting, 
    performs identical operations on RHS vector.
    Args:
        A (array_like): A n-by-n regular matrix
        b (array_like): RHS vector of size n
    Returns:
        numpy.ndarray: Upper triangular matrix
        numpy.ndarray: RHS vector corresponding to upper triangular matrix
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the scipy.linalg.solve function.

A = np.random.rand(3, 3)
b = np.random.rand(3)
A_upper, bb = gaussian_elimination_with_pivoting_vector(A, b)
try:
    np.testing.assert_array_almost_equal(backward_substitution(A_upper, bb), la.solve(A, b), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Note: The code above uses Gaussian elimination to transform a given random matrix into an upper triangular form, solves a system of linear algebraic equations with the transformed matrix using the backward substitution, and compare the result with the one obtained by scipy.linalg.solve function.

LU decomposition¶

LU decomposition factors a matrix as the product of a lower triangular matrix and an upper triangular matrix. LU decomposition can be achieved by Crout algorithm,

\[ u_{ij} = a_{ij} - \sum_{k = 1}^{i - 1} l_{ik} u_{kj}, \]

\[ l_{ji} = \frac{1}{u_{ii}} \left( a_{ji} - \sum_{k = 1}^{i - 1} l_{jk} u_{ki} \right), \]

\[ i = 1, 2, \dots, n, \quad j = i, i + 1, \dots, n. \]

A system of linear algebraic equations \( \mathbb{A} x = b \) can be solved using LU decomposition as follows,

\[ \mathbb{A} x = (\mathbb{L} \mathbb{U}) x = \mathbb{L} ( \mathbb{U} x) = b, \]

\[ \mathbb{L} y = b, \]

\[ \mathbb{U} x = y. \]

Exercise 03.9: Implement the LU decomposition algorithm.

def lu_decomposition(A):
    """
    Factors given matrix as the product of a lower and an upper triangular matrix using LU decomposition.
    Args:
        A (array_like): A n-by-n matrix
    Returns:
        numpy.ndarray: Lower triangular matrix
        numpy.ndarray: Upper triangular matrix
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the scipy.linalg.solve function.

A = np.random.rand(3, 3)
b = np.random.rand(3)
L, U = lu_decomposition(A)
try:
    np.testing.assert_array_almost_equal(backward_substitution(U, forward_substitution(L, b)), la.solve(A, b), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Note: The code above decomposes a given random matrix using the LU decomposition, solves the systems of algebraic equations using the forward and backward substitutions, and compare the result with the one obtained by scipy.linalg.solve function.

Thomas algorithm¶

A system of linear algebraic equations \( \mathbb{A} x = b \) with tridiagonal coefficient matrix \( \mathbb{A} \) can be efficiently solved by the Thomas algorithm. Tridiagonal matrix \( \mathbb{A} \) can be represented by three vectors \( p \), \( q \), and \( r \), as follows,

\[\begin{split} \mathbb{A} = \begin{pmatrix} q_0 & p_0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ r_1 & q_1 & p_1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & r_2 & q_2 & p_2 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & r_{n-2} & q_{n-2} & p_{n-2} \\ 0 & 0 & 0 & 0 & \cdots & 0 & r_{n-1} & q_{n-1} \end{pmatrix}, \end{split}\]

where \( p_{n-1} = 0 \) and \( r_0 = 0 \). The solution of \( \mathbb{A} x = b \) is then obtained iteratively as

\[ x_i = \mu_i x_{i+1} + \rho_i, \qquad i = n-1, n-2, \dots, 0, \]

where the coefficients \( \mu_i \) and \( \rho_i \) are

\[ \mu_i = -\frac{p_i}{r_i \mu_{i-1} + q_i}, \qquad \rho_i = \frac{b_i - r_i \rho_{i-1}}{r_i \mu_{i-1} + q_i}, \qquad i = 1, 2, \dots, n-1 \]

with \( \mu_0 = -p_0 \, / \, q_0 \) and \( \rho_0 = b_0 \, / \, q_0 \).

Exercise 03.10: Implement the Thomas algorithm.

def thomas_algorithm(A, b):
    """
    Solves system of linear equations with a tridiagonal matrix using Thomas algorithm.
    Args:
        A (array_like): A n-by-n regular matrix
        b (array_like): RHS vector of size n
    Returns:
        numpy.ndarray: Vector of solution
    """
    
    # add your code here

You may evaluate the correctness of your implementation using the scipy.linalg.solve function.

A = np.triu(np.tril(np.random.rand(5, 5), 1), -1) # create a tridiagonal matrix from random matrix
b = np.random.rand(5)
try:
    np.testing.assert_array_almost_equal(thomas_algorithm(A, b), la.solve(A, b), decimal=7)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

12NME1 - Numerical methods

Tutorial 03 - Linear algebra

Contents

Tutorial 03 - Linear algebra¶

Basic operations¶

Scalar product¶

Matrix-vector product¶

Matrix-matrix product¶

Direct methods for linear algebraic equations¶

Forward substitution¶

Backward substitution¶

Gaussian elimination¶

LU decomposition¶

Thomas algorithm¶