13. Basic matrix algebra in numpy#

Saleh Rezaeiravesh, saleh.rezaeiravesh@manchester.ac.uk
Department of Mechanical and Aerospace Engineering, The University of Manchester, Manchester, UK

Overview: In this notebook, we explore tools provided by `linalg`, a submodule of `numpy`. Specifically, we cover functions for computing the norm, determinant, and inverse of arrays/matrices. Additionally, we demonstrate how to write scripts to perform these computations directly.

For further details and additional methods, please refer to the official numpy.linalg documentation:

https://numpy.org/devdocs/reference/routines.linalg.html

13.1. Intended Learning Outcomes#

After reading and running this notebook, the students should be able to:

  • Understand how to import the numpy subclass linalg and its methods.

  • Use numpy.linalg for computing norm, determinant and inverse of matrices.

  • Implement basic algorithms in Python for computing norm, determinant and inverse of matrices.

13.2. Import the library#

We import numpy in order to be able to define arrays and apply operations on them.

import numpy as np

To import linalg submodule in a Python script, we have two options:

  • Option 1: Import linalg as a whole. Then in the script, we can call specific methods (functions) from that.

  • Option 2: Import only the methods that needed in the script from numpy.linalg.

Both of these are used in practice, and you can decide which one to use. But within the present notebook, we adopt the first option.

# Option 1: import the whole linalg
import numpy.linalg as la

# Option 2: Import a specific method from linalg
from numpy.linalg import norm

13.3. Norm of arrays#

Mathematically speaking, a norm of an array is a non-negative single value that represents the “size” or “magnitude” of that array.

There are various ways and definitions for measuring the norm of an array. Below, we provide definitions of some standard norms for 1D and 2D arrays, and then explain how they can be avaulated in Python.

13.3.1. Mathematical definitions#

13.3.1.1. 1D Arrays (Vectors)#

Consider a 1D array \(\mathbf{a}\) of size \(n\):

\[ \mathbf{a} = [a_1,a_2,\cdots,a_n]\,,\quad \mathbf{a}=[a_i]\,,\quad \text{where}\,\, i=1,2,\cdots,n \]

Some most frequently-used norms for the vectors are given bellow:

  • \(L_1\) norm: sum of the absolute values of elements

\[ \|\mathbf{a}\|_1 = \sum_{i=1}^n |a_i| \]

  • \(L_2\) norm: square root of the sum of squares of elements

\[ \|\mathbf{a}\|_2 = \left (\sum_{i=1}^n |a_i|^2 \right)^{1/2} \]

  • \(L_p\) norm: for \(p\geq 1\), and \(p\) integer, this is a general form that includes the above two:

\[ \|\mathbf{a}\|_p = \left (\sum_{i=1}^n |a_i|^p \right)^{1/p} \]

  • \(L_\infty\) norm: maximum absolute value of the elements

\[ \|\mathbf{a}\|_\infty = \max_{1\leq i\leq n} |a_i| \]

13.3.1.2. 2D Arrays#

Consider a general \(m\times n\) matrix \(\mathbf{A}\):

\[\begin{split} \mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \ddots & \cdots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \\ \end{bmatrix} \,,\quad \mathbf{A}=[a_{ij}] \,\, \text{where} \,\, i=1,2,\cdots,m \,,\, j=1,2,\cdots,n \end{split}\]

For this matrix, the above norms are defined as:

  • \(L_1\) norm: maximum column-sum of the absolute values of elements:

\[ \|\mathbf{A}\|_1 = \max_{1\leq j\leq n}\sum_{i=1}^m |a_{ij}| \]

  • \(L_2\) norm: square root of the maximum eigenvalue of \(\mathbf{A^*} \mathbf{A}\) where \(A^*\) is the conjugate transpose of \(\mathbf{A}\):

\[ \|\mathbf{A}\|_2 = \sqrt{\lambda_{\max} (\mathbf{A}^*\mathbf{A})} \]

  • \(L_\infty\) norm: maximum row-sum of the absolute values of elements:

\[ \|\mathbf{A}\|_\infty = \max_{1\leq i\leq m}\sum_{j=1}^n |a_{ij}| \]

13.3.2. Computing Norms#

13.3.2.1. numpy built-in function#

To find the norm of an array in numpy, we use numpy.linalg.norm(a,s) where s can be

  • non-negative integer: to compute \(L_p\) norm,

  • numpy.inf: to compute \(L_\infty\) norm.

Note that the above function is the same for arrays of any size and dimension.

Example: Norm of a real 1D array:

a = np.array([-1.6,3.5,4.2,-5.8,9.7])

print('L1-norm of a:',la.norm(a,1))  
print('L2-norm of a:',la.norm(a,2))
print('L5-norm of a:',la.norm(a,5))
print('L infiniy-norm of a:',la.norm(a,np.inf))

L1-norm of a: 24.8
L2-norm of a: 12.656223765404908
L5-norm of a: 9.882880278755472
L infiniy-norm of a: 9.7

Example: Norm of a complex 1D array:

a = np.array([-1.6+3.j, 3.5+1.j,4.2-3.8j, -5.8+4.6j, 9.7+0.5j])

print('L1-norm of a:',la.norm(a,1))  
print('L2-norm of a:',la.norm(a,2))
print('L5-norm of a:',la.norm(a,5))
print('L infiniy-norm of a:',la.norm(a,np.inf))

L1-norm of a: 29.81955610664194
L2-norm of a: 14.353745155881791
L5-norm of a: 10.29413759208615
L infiniy-norm of a: 9.712878049270463

Example: Norms of a 2D array:

A = np.array([[-1.,3.2,4.2-3.7j],[2.0+5.6j,3.j,6.0+2.j]])

print(A)
print('L1-norm of A:',la.norm(A,1))
print('L2-norm of A:',la.norm(A,2))
print('Linf-norm of A:',la.norm(A,np.inf))

[[-1. +0.j   3.2+0.j   4.2-3.7j]
 [ 2. +5.6j  0. +3.j   6. +2.j ]]
L1-norm of A: 11.92187610799937
L2-norm of A: 10.337464314431093
Linf-norm of A: 15.270982819264162

13.3.2.2. Direct implementation#

We can also write Python function to compute different norms of input arrays.

Here we write a Python function that computes the \(L_p\) norm of 1D arrays:

def LpNorm(x,p):
    """
    Finds l-p norm of 1D array x
    
    Args:
       `x`: 1d numpy array
       `p`: int
       
    Returns:
       `f`: float, ||x||_p       
    """
    n = len(x)   #size of the array
    f = 0.0      #initialize the norm
    
    for i in range(n):
        f += abs(x[i])**p     #update norm by adding x_i^p
        
    f = f**(1./p)             #compute the norm by f^(1/p)
    
    return f

Let’s test this function and validate its results against the numpy built-in function:

a = np.array([-1.6+3.j, 3.5+1.j,4.2-3.8j, -5.8+4.6j, 9.7+0.5j])

print('numpy L1-norm of a:',la.norm(a,1))  
print('numpy L2-norm of a:',la.norm(a,2))
print('numpy L5-norm of a:',la.norm(a,5))

numpy L1-norm of a: 29.81955610664194
numpy L2-norm of a: 14.353745155881791
numpy L5-norm of a: 10.29413759208615

print('Our L1-norm of a:',LpNorm(a,1))  
print('Our L2-norm of a:',LpNorm(a,2))
print('Our L5-norm of a:',LpNorm(a,5))

Our L1-norm of a: 29.81955610664194
Our L2-norm of a: 14.353745155881791
Our L5-norm of a: 10.29413759208615

Next, we write a Python function that computes the \(L_\infty\) norm of 1D arrays:

def LinfNorm(x):
    """
    Finds l-infinity norm of 1D array x
    
    Args:
       `x`: 1d numpy array
       
    Returns:
       `f`: float, norm
    """
    n = len(x)   #size of the array
    f = 0.0      #initialize the norm
    
    for i in range(n):
        if abs(x[i]) > f:
           f= abs(x[i]) 
    
    return f

Let’s test this function and validate its results against the numpy built-in function:

print('numpy L infiniy-norm of a:',la.norm(a,np.inf))
print('Our L infiniy-norm of a:',LinfNorm(a))

numpy L infiniy-norm of a: 9.712878049270463
Our L infiniy-norm of a: 9.712878049270463

Exercise: Can you combine the above two function, so it returns \(L_p\) or \(L_\infty\) norm as chosen by the user?

13.4. Determinant of square matrices#

Determinant of a \(n\times n\) matrix \(\mathbf{A}\) is defined as,

\[ \det(\mathbf{A}) = \sum_{j=1}^n (-1)^{i+j} a_{ij} m_{ij} \]

where \(m_{ij}\) is called the minor of the element \(a_{ij}\) and is computed as the determinant of \(\mathbf{A}\) after removing its \(i\)-th row and \(j\)-th column.

Note that, although the above expansion has been written around the \(i\)-th row of \(\mathbf{A}\), we can alternatively find the determinant of \(\mathbf{A}\) by expanding around its \(j\)-th column.

Below, we show how to compute the determinant of an array using numpy, and we also write a Python function based on the above definition.

13.4.1. numpy built-in function#

To find the determinant of a square matrix, use numpy.linalg.det(A):

A = np.array([[-1.,3.,4.],[2.,3.,6.],[6.,-9.,0.]])

print(A)

detA = la.det(A)
print('det(A)=',detA)

[[-1.  3.  4.]
 [ 2.  3.  6.]
 [ 6. -9.  0.]]
det(A)= -90.0

We can also apply this function to matrices with complex elements:

B = np.array([[-1.j,3.-2.j,4.],[2.,3.,6.-3.j],[6.+7.j,-9.,0.]])

print(B)

detB = la.det(B)
print('det(B)=',detB)

[[-0.-1.j  3.-2.j  4.+0.j]
 [ 2.+0.j  3.+0.j  6.-3.j]
 [ 6.+7.j -9.+0.j  0.+0.j]]
det(B)= (48-180.0000000000001j)

13.4.2. Direct implementation#

We can implement the above definition to compute the determinant of a square matrix. First, we write a function to compute the minors \(m_{ij}\) of \(\mathbf{A}\) for given indices \(i\) and \(j\). Note that,

  • To remove the elements on row \(i\) of \(\mathbf{A}\), we use numpy.delete(A,i,0), where 0 refers to axis=0 (rows) of the array.

  • To remove the elements on column \(j\) of \(\mathbf{A}\), we use numpy.delete(A,j,1), where 1 refers to axis=1 (columns) of the array.

def minor(A,i,j):
    """
    Find the minor m_{ij} of array A
    """
    B = np.delete(A,i,0)   #remove the i-th row of A
    B = np.delete(B,j,1)   #remove the j-th column of A
    m = la.det(B)          #compute minor mij
    return m

Given the minors \(m_{ij}\), we implement the above expression in the following function to compute the determninant of \(\mathbf{A}\).

def determinant(A,i=0):
    """
    Find the determinant of square matrix A by expanding around its i-th row
    """
    n = A.shape[0] #size of the array
    if n != A.shape[1] or A.ndim < 2:
       raise ValueError('A should be a square matrix!')
        
    s = 0.   #initialize the determinant
    for j in range(n):
        mij = minor(A,i,j)
        s += (-1)**(i+j) * A[i,j] * mij
        
    return s

Now, call our function determinant to compute the determinant of \(\mathbf{A}\). In the following lines, the determinant is computed by expanding around different rows. If no row index is provided, the default value of 0 (first row) is considered.

print('1st row expand, det(A)=',determinant(A,0))
print('2nd row expand, det(A)=',determinant(A,1))
print('3rd row expand, det(A)=',determinant(A,2))

print('numpy det(A)=',la.det(A))

1st row expand, det(A)= -90.0
2nd row expand, det(A)= -89.99999999999999
3rd row expand, det(A)= -89.99999999999997
numpy det(A)= -90.0

Do you see any difference in the value of the \(\det(\mathbf{A})\) using different rows in the expansion? Why do you think this happens? Expansion around which row do you think is considered by default in the numpy.linalg.det() function?

Exercise: Modify the function determinant to compute the determinant of a matrix \(\mathbf{A}\) by expanding around its \(j\)-th column.

13.5. Inverse of square matrices#

Consider the \(n\times n\) matrix \(\mathbf{A}\). The inverse of this matrix is obtained from,

\[ \mathbf{A}^{-1} = \frac{adj(\mathbf{A})}{\det{\mathbf{A}}} \,, \]

where \(adj(\mathbf{A})\) is the adjoint of \(\mathbf{A}\), defined as,

\[\begin{split} adj(\mathbf{A}) = \begin{bmatrix} c_{11} & c_{12} & \cdots & c_{1n} \\ c_{21} & c_{22} & \cdots & c_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ c_{n1} & c_{n2} & \cdots & c_{nn} \\ \end{bmatrix}^T \,, \end{split}\]

with \(c_{ij} = (-1)^{i+j} m_{ij}\).

Note that if \(\det(\mathbf{A}) = 0\), then \(\mathbf{A}^{-1}\) cannot be computed (\(\mathbf{A}\) is not invertible).

13.5.1. numpy built-in function#

To find the inverse of a square matrix, we can use numpy.linalg.inv(A):

print('A =',A)

A_inv = la.inv(A)

print('inv(A)=',A_inv)

A = [[-1.  3.  4.]
 [ 2.  3.  6.]
 [ 6. -9.  0.]]
inv(A)= [[-0.6         0.4        -0.06666667]
 [-0.4         0.26666667 -0.15555556]
 [ 0.4        -0.1         0.1       ]]

We can check if \(\mathbf{A} \mathbf{A}^{-1} = \mathbf{I}\), where \(\mathbf{I}\) is the identity matrix of the same size as \(\mathbf{A}\):

print(A @ A_inv)

[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 2.22044605e-16  1.00000000e+00  5.55111512e-17]
 [ 5.55111512e-16 -5.55111512e-17  1.00000000e+00]]

The above identity holds, although noting some of the off-diagonal elements of \(\mathbf{I}\) are not exactly zero. Why?

13.5.2. Direct implementation#

We can write a Python code that computes \(\mathbf{A}^{-1}\) following the above definition, see function inverse below. To this end, we need the following to be computed:

  • \(\det(\mathbf{A})\): This is implemented in the function determinant, see above.

  • \(m_{ij}\): Minors of elements of \(\mathbf{A}\) can be computed from function minor, see above.

  • \(adj(\mathbf{A})\): See the function below.

def adjoint(A):
    """
    Compute adjoint(A)
    """
    n=A.shape[0]
    
    adj = np.zeros((n,n))   #initialize adj(A)
    for j in range(n):      #columns of A, rows of adj(A)
        for i in range(n):  #rows of A, columns of adj(A)
            adj[j,i] = (-1)**(i+j)*minor(A,i,j)   
            
    return adj

Now, we have all the components to compute \(\mathbf{A}^{-1}\) using the following function:

def inverse(A):
    """
    Compute the inverse of A
    """
    n= A.shape[0]
    
    #check if the matrix A is square
    if n != A.shape[1] or A.ndim < 2:
       raise ValueError('A should be a square matrix!')
    
    #compute the determinant of A
    detA = determinant(A,0)
    
    #compute the adjoint of A
    adjA = adjoint(A)
    
    #compute the inverse of A
    if detA == 0:
       raise ErrorValue("A is not invertible!")         
    else:
       invA = adjA/detA
    
    return invA

Let’s validate our implementation by comparing its output to that of numpy.linalg.inv():

print('A =',A)

A_inv_numpy = la.inv(A)   #numpy inverse
A_inv_ours = inverse(A)   #our implementation

print('numpy inv(A):')
print(A_inv_numpy)

print('Our inv(A):')
print(A_inv_ours)

A = [[-1.  3.  4.]
 [ 2.  3.  6.]
 [ 6. -9.  0.]]
numpy inv(A):
[[-0.6         0.4        -0.06666667]
 [-0.4         0.26666667 -0.15555556]
 [ 0.4        -0.1         0.1       ]]
Our inv(A):
[[-0.6         0.4        -0.06666667]
 [-0.4         0.26666667 -0.15555556]
 [ 0.4        -0.1         0.1       ]]