Matrices in Python

We love matrices. I know what you’re thinking “ewwe no we don’t”. But yes, we do. Because when we use matrices, we dont have to do maths – the computer can do it all for us. Big win.

• Setting up a matrix
  • Making a square matrix
• Matrix functions
  • Transpose
  • Multiplication
• Gaussian Elimination

Setting up a Matrix

A list of integers or floats can be converted into an m x n matrix, where m is the number of rows and n the number of columns, using for loops or NumPy.

Without NumPy, using just basic loops and no library inputs, the code looks something like this:

#Defining an n x m matrix from a list
def Matrix(m, n, List):
    #defining an empty row:
    row = []
    matrix = [0 for height in range(0,m)]
    #setting rows:
    for p in range(0, m, 1):
        #setting columns:
        for i in range((p*n), ((p+1)*n), 1):
            #adding each element of a row to the row:
            row = row + [List[i]]
        matrix[p] = row
        row = []
    return matrix

This will output a list with multiple sub lists. They will not print on separate lines, so the matrix will look like [[a,b,c], [d,e,f], [g,h,i]]

However, if you use NumPy, it is a lot shorter:

#defining a matrix from a list using NumPy
import numpy as np
matrix = np.array(List).reshape(m,n)

NumPy arrays automatically format into rows and columns when printed, so appear more familiar and intuitive.

Making a Square Matrix

We can use the NumPy library to easily make a square matrix from a list of integers or floats:

def SquareMatrix(List):
    n = len(List)**0.5
    if int(len(List)**0.5)==n:
        matrix = np.array(List).reshape(int(n),int(n))
        return matrix
    else:
        return 'Cannot produce square matrix'

Matrix Functions

The following functions are pre-programmed in NumPy, where the matrix is assigned the name M. Always remember to import NumPy!

import numpy as np

However, if you want to see the code without using libraries, click on the ’empirical code’ tab below each snippet.

Transpose

np.transpose(M)

def Transpose(M):
    n = len(M)
    m = len(M[0])
    row = []
    transpose = [0 for height in range(0, m)]
    #B is the column number
    for B in range(0, m):
        #A is the row number
        for A in range(0, n):
            row = row + [M[A][B]]
        transpose[B] = row
        #reset row
        row = []
    return transpose

Multiplying Two Matrices

The two matrices being multiplied are A and B.

np.matmul(A,B)

def MatMat(A,B):
    #check if matricies can be multiplied
    if len(A[0]) == len(B):
        add = 0
        C = [[0 for width in range(len(B[0]))] for height in range(len(A))]
        for x in range(0, len(A)):
            for y in range(0, len(B[0])):
                for z in range(0, len(A[0])):
                    add = add + (A[x][z]*B[z][y])
                C[x][y] = add
                add = 0
        return C
    else:
    #if matrix cannot be multiplied
        return "Matricies cannot be multiplied."

Gaussian Elimination

The Gaussian Elimination is a method used to solve systems of linear equations in the form:

This system can be written in terms of three matrices:

The matrix on the left (a-i) is and that on the right (α-γ) is .

Assuming we know the values in these matrices, Gaussian elimination can be used to find , , and .

def Gaussian(M1,M2):
    import numpy as np
    
    #m is the height of the matrix
    m = len(M1)
    #n is the width of the matrix
    n = len(M1[0])
    
    if len(M1)!=len(M2) or len(M2[0])!=1:
        return 'The matrices are of incorrect dimensions'
    
    else:
        #making augmented matrix as numpy array
        A = M1
        for i in range(0,n):
            A[i] = A[i] + M2[i]
        A = np.asarray(A)
        
        #finding upper triangular form
        M = np.zeros((n,m+1))
        M[0] = A[0]
        for j in range(0,n):
            for i in range(j+1,m):
                M[i] = A[i]-A[i][j]*A[j]/A[j][j]
            A = M
        
        #finding lower triangular form
        M = np.zeros((n,m+1))
        M[-1] = A[-1]
        for j in range(-2,-n-1,-1):
            for i in range(j,-m-1,-1):
                M[i] = A[i]-A[i][j]*A[j+1]/A[j+1][j]
            A = M
        
        #cleaning up, removing storage errors
        for i in range(0,m):
            for j in range(0,n):
                if -1e-15 < A[i][j] < 1e-15:
                    A[i][j]=0
        
        #identifying
        for j in range(0,m):
            A[j][j] = A[j][-1]/A[j][j]
        
        #making column matrix
        S = np.zeros((m,1))
        for i in range(0,m):
            S[i] = A[i][i]
        return S

M1 = [[8,-2,1,3],[1,-5,2,1],[-1,2,7,2],[2,-1,3,8]]
M2 = [[9],[-7],[-1],[5]]
print(Gaussian(M1,M2))
#the output it prints is:
[[ 1.32323232]
 [ 1.56565657]
 [-0.60606061]
 [ 0.71717172]]

Using Gaussian Elimination to find an Inverse

Gaussian elimination can be used to find the inverse of a matrix by augmenting the matrix with the identity matrix:

The function in python is:

def GaussianInverse(M1):
    import numpy as np
    
    #m is the height of the matrix
    m = len(M1)
    #n is the width of the matrix
    n = len(M1[0])
    
    if len(M1)!=len (M1[0]):
        return 'The matrix is not square.'
    
    else:
        #making right zeros matrix
        I = [0 for i in range(m)]
        I = [I for i in range(m)]
        
        #making augmented matrix as numpy array
        A = M1
        for i in range(0,n):
            A[i] = A[i] + I[i]
        A = np.asarray(A)
        #making right identity matrix
        for i in range(m):
            A[i][i+m]=1
        #finding upper triangular form
        M = np.zeros((n,2*m))
        M[0] = A[0]
        for j in range(0,n):
            for i in range(j+1,m):
                M[i] = A[i]-A[i][j]*A[j]/A[j][j]
            A = M
        #finding lower triangular form
        for j in range(-1,-m-1,-1):
            for i in range(j-1,-m-1,-1):
                M[i] = A[i]-A[i][-m+j]*A[j]/A[j][-m+j]
            A = M
        #cleaning up, removing storage errors
        for i in range(0,m):
            for j in range(0,2*n):
                if -1e-15 < A[i][j] < 1e-15:
                    A[i][j]=0
        
        for i in range(m):
            A[i] = A[i]/A[i][i]
        
        return np.delete(A,range(m),1)