Please find the code for this post on GitHub. As always, I hope you’ll clone it and make it your own. The main module in the repo that holds all the modules that we’ll cover is named LinearAlgebraPurePython.py. There’s a simple python file named BasicToolsPractice.py that imports that main module and illustrates the modules functions. It’d be great if you could clone or download that first to have handy as we go through this post.

Why This Post?

REMINDER: Our goal is to better understand principles of machine learning tools by exploring how to code them ourselves …

Meaning, we are seeking to code these tools without using the AWESOME python modules available for machine learning.

These efforts will provide insights and better understanding, but those insights won’t likely fly out at us every post. Rather, we are building a foundation that will support those insights in the future.

To streamline some upcoming posts, I wanted to cover some basic functions that will make those future posts easier. Obviously, if we are avoiding using numpy and scipy, we’ll have to create our own convenience functions / tools. This post covers those convenience tools. 

At one end of the spectrum, if you are new to linear algebra or python or both, I believe that you will find this post helpful among, I hope, a good group of saved links.

At the other end of the spectrum, if you have background with python and linear algebra, your reason to read this post would be to compare how I did it to how you’d do it. The review may give you some new ideas, or it may confirm that you still like your way better. 

The Tools

You’ll find documentation and comments in all of these functions. When more description is warranted, I will give it or provide directions to other resource to describe it in more detail. If there is a specific part you don’t understand, I am eager for you to understand it better. What’s the best way to do that? Rebuild these functions from the inner most operations yourself and experiment with them at that level until you understand them, and then add the next layer of looping, or code that repeats that inner most operation, and understand that, etc.

First up is zeros_matrix. When we just need a new matrix, let’s make one and fill it with zeros.

def zeros_matrix(rows, cols):
    """
    Creates a matrix filled with zeros.
        :param rows: the number of rows the matrix should have
        :param cols: the number of columns the matrix should have

        :return: list of lists that form the matrix
    """
    M = []
    while len(M) < rows:
        M.append([])
        while len(M[-1]) < cols:
            M[-1].append(0.0)

    return M

As you’ve seen from the previous posts, matrices and vectors are both being handled in Python as two dimensional arrays. Thus, the array of rows contains an array of the column values, and each column value is initialized to 0. Notice the -1 index to the matrix row in the second while loop. This is a simple way to reference the last element of an array, and in this case, it’s the last array (row) that’s been appended to the array.

Our Second helper function is identity_matrix used to create an identity matrix. And, as a good constructively lazy programmer should do, I have leveraged heavily on an initial call to zeros_matrix. All that’s left once we have an identity matrix is to replace the diagonal elements with 1.

def identity_matrix(n):
    """
    Creates and returns an identity matrix.
        :param n: the square size of the matrix

        :return: a square identity matrix
    """
    IdM = zeros_matrix(n, n)
    for i in range(n):
        IdM[i][i] = 1.0

    return IdM

Third is copy_matrix also relying heavily on zeros_matrix. We want this for those times where we need to work on a copy and preserve the original matrix. Here, we are simply getting the dimensions of the original matrix and using those dimensions to create a zeros matrix and then copying the elements of the original matrix to the new matrix element by element.

def copy_matrix(M):
    """
    Creates and returns a copy of a matrix.
        :param M: The matrix to be copied

        :return: A copy of the given matrix
    """
    # Section 1: Get matrix dimensions
    rows = len(M)
    cols = len(M[0])

    # Section 2: Create a new matrix of zeros
    MC = zeros_matrix(rows, cols)

    # Section 3: Copy values of M into the copy
    for i in range(rows):
        for j in range(cols):
            MC[i][j] = M[i][j]

    return MC

Fourth is print_matrix so that we can see if we’ve messed up or not in our linear algebra operations! Here, we are just printing the matrix, or vector, one row at a time. The “+0” in the list comprehension was mentioned in a previous post. Try the list comprehension with and without that “+0” and see what happens. 

def print_matrix(M, decimals=3):
    """
    Print a matrix one row at a time
        :param M: The matrix to be printed
    """
    for row in M:
        print([round(x,decimals)+0 for x in row])

In case you don’t yet know python list comprehension techniques, they are worth learning. There are tons of good blogs and sites that teach it. As I always, I recommend that you refer to at least three sources when picking up any new skill but especially when learning a new Python skill. Some brief examples would be …

some_new_list = [<method(x)> for x in list if <condition> else <other>]
# or
another_list = [s.method() for s in string_list if <condition>]
# or
one_more_list = [<method(x)> for x in list]

The point of showing one_more_list is to make it abundantly clear that you don’t actually need to have any conditionals in the list comprehension, and the method you apply can be one that you write.

Fifth is transpose. Transposing a matrix is simply the act of moving the elements from a given original row and column to a  row = original column and a column = original row. That is, if a given element of M is m_{i,j}, it will move to m_{j,i} in the transposed matrix, which is shown as

MT[j][i] = M[i][j]

in the code. In relation to this principle, notice that the zeros matrix is created with the original matrix’s number of columns for the transposed matrix’s number of rows and the original matrix’s number of rows for the transposed matrix’s number of columns. What a mouthful!

Notice that in section 1 below, we first make sure that M is a two dimensional Python array. Then we store the dimensions of M in section 2. Next, in section 3, we use those dimensions to create a zeros matrix that has the transposed matrix’s dimensions and call it MT. Finally, in section 4, we transfer the values from M to MT in a transposed manner as described previously.

def transpose(M):
    """
    Returns a transpose of a matrix.
        :param M: The matrix to be transposed

        :return: The transpose of the given matrix
    """
    # Section 1: if a 1D array, convert to a 2D array = matrix
    if not isinstance(M[0],list):
        M = [M]

    # Section 2: Get dimensions
    rows = len(M)
    cols = len(M[0])

    # Section 3: MT is zeros matrix with transposed dimensions
    MT = zeros_matrix(cols, rows)

    # Section 4: Copy values from M to it's transpose MT
    for i in range(rows):
        for j in range(cols):
            MT[j][i] = M[i][j]

    return MT

Sixth and Seventh are matrix_addition and matrix_subtraction. I am explaining them at the same time, because they are essentially identical with the exception of the single line of code where the element by element additions or subtractions take place. In section 1 of each function, you see that we check that each matrix has identical dimensions, otherwise, we cannot add them. Section 2 of each function creates a zeros matrix to hold the resulting matrix. Section 3 of each function performs the element by element operation of addition or subtraction, respectively. 

def matrix_addition(A, B):
    """
    Adds two matrices and returns the sum
        :param A: The first matrix
        :param B: The second matrix

        :return: Matrix sum
    """
    # Section 1: Ensure dimensions are valid for matrix addition
    rowsA = len(A)
    colsA = len(A[0])
    rowsB = len(B)
    colsB = len(B[0])
    if rowsA != rowsB or colsA != colsB:
        raise ArithmeticError('Matrices are NOT the same size.')

    # Section 2: Create a new matrix for the matrix sum
    C = zeros_matrix(rowsA, colsB)

    # Section 3: Perform element by element sum
    for i in range(rowsA):
        for j in range(colsB):
            C[i][j] = A[i][j] + B[i][j]

    return C

def matrix_subtraction(A, B):
    """
    Subtracts matrix B from matrix A and returns difference
        :param A: The first matrix
        :param B: The second matrix

        :return: Matrix difference
    """
    # Section 1: Ensure dimensions are valid for matrix subtraction
    rowsA = len(A)
    colsA = len(A[0])
    rowsB = len(B)
    colsB = len(B[0])
    if rowsA != rowsB or colsA != colsB:
        raise ArithmeticError('Matrices are NOT the same size.')

    # Section 2: Create a new matrix for the matrix difference
    C = zeros_matrix(rowsA, colsB)

    # Section 3: Perform element by element subtraction
    for i in range(rowsA):
        for j in range(colsB):
            C[i][j] = A[i][j] - B[i][j]

    return C

Eighth is matrix_multiply. The first rule in matrix multiplication is that if you want to multiply matrix A times matrix B, the number of columns of A MUST equal the number of rows of B. Thus, if A has dimensions of m rows and n columns (m\,x\,n for short) B must have n rows and it can have 1 or more columns. Let’s say it has k columns. Thus, the resulting product of the two matrices will be an m\,x\,k matrix, or the resulting matrix has the number of rows of A and the number of columns of B.  Hence, we create a zeros matrix to hold the resulting product of the two matrices that has dimensions of rows_A \, x \, cols_B in the code. Also, IF A and B have the same dimensions of n rows and n columns, that is they are square matrices, A \cdot B does NOT equal B \cdot A. Remember that the order of multiplication matters when multiplying matrices. Finally, the result for each new element c_{i,j} in C, which will be the result of A \cdot B, is found as follows using a 3\,x\,3 matrix as an example:

c_{i,j} = a_{i,0} \cdot b_{0,j} + a_{i,1} \cdot b_{1,j} + a_{i,2} \cdot b_{2,j}

That is, to get c_{i,j} we are multiplying each column element in each row i of A times each row element in each column j of B and adding up those products. Phew!

def matrix_multiply(A, B):
    """
    Returns the product of the matrix A * B
        :param A: The first matrix - ORDER MATTERS!
        :param B: The second matrix

        :return: The product of the two matrices
    """
    # Section 1: Ensure A & B dimensions are correct for multiplication
    rowsA = len(A)
    colsA = len(A[0])
    rowsB = len(B)
    colsB = len(B[0])
    if colsA != rowsB:
        raise ArithmeticError(
            'Number of A columns must equal number of B rows.')

    # Section 2: Store matrix multiplication in a new matrix
    C = zeros_matrix(rowsA, colsB)
    for i in range(rowsA):
        for j in range(colsB):
            total = 0
            for ii in range(colsA):
                total += A[i][ii] * B[ii][j]
            C[i][j] = total

    return C

Ninth is a function, multiply_matrices, to multiply out a list of matrices using matrix_multiply. Note that we simply establish the running product as the first matrix in the list, and then the for loop starts at the second element (of the list of matrices) to loop through the matrices and create the running product, matrix_product, times the next matrix in the list. 

def multiply_matrices(list):
    """
    Find the product of a list of matrices from first to last
        :param list: The list of matrices IN ORDER

        :return: The product of the matrices
    """
    # Section 1: Start matrix product using 1st matrix in list
    matrix_product = list[0]

    # Section 2: Loop thru list to create product
    for matrix in list[1:]:
        matrix_product = matrix_multiply(matrix_product, matrix)

    return matrix_product

Tenth, and I confess I wasn’t sure when it was best to present this one, is check_matrix_equality. There will be times where checking the equality between two matrices is the best way to verify our results. However, those operations will have some amount of round off error to where the matrices won’t be exactly equal, but they will be essentially equal. Thus, note that there is a tol (tolerance parameter), that can be set. If the default is used, the two matrices are expected to be exactly equal. If a tolerance is set, the value of tol is the number of decimal places the element values are rounded off to to check for an essentially equal state.

def check_matrix_equality(A, B, tol=None):
    """
    Checks the equality of two matrices.
        :param A: The first matrix
        :param B: The second matrix
        :param tol: The decimal place tolerance of the check

        :return: The boolean result of the equality check
    """
    # Section 1: First ensure matrices have same dimensions
    if len(A) != len(B) or len(A[0]) != len(B[0]):
        return False

    # Section 2: Check element by element equality
    #            use tolerance if given
    for i in range(len(A)):
        for j in range(len(A[0])):
            if tol == None:
                if A[i][j] != B[i][j]:
                    return False
            else:
                if round(A[i][j],tol) != round(B[i][j],tol):
                    return False

    return True

The dot product between two vectors or matrices is essentially matrix multiplication and must follow the same rules. It’s important to note that our matrix multiplication routine could be used to multiply two vectors that could result in a single value matrix. In such cases, that result is considered to not be a vector or matrix, but it is single value, or scaler. However, using our routines, it would still be an array with a one valued array inside of it. To read another reference, check HERE, and I would save that link as a bookmark – it’s a great resource.

The Eleventh function is the unitize_vector function. Let’s step through its sections. Section 1 ensures that a vector was input meaning that one of the dimensions should be 1. Also, it makes sure that the array is 2 dimensional. This tool kit wants all matrices and vectors to be 2 dimensional for consistency. Section 2 uses the Pythagorean theorem to find the magnitude of the vector.  Section 3 makes a copy of the original vector (the copy_matrix function works fine, because it still works on 2D arrays), and Section 4 divides each element by the determined magnitude of the vector to create a unit vector. 

def unitize_vector(vector):
    """
    Find the unit vector for a vector
        :param vector: The vector to find a unit vector for

        :return: A unit-vector of vector
    """
    # Section 1: Ensure that a vector was given
    if len(vector) > 1 and len(vector[0]) > 1:
        raise ArithmeticError(
            'Vector must be a row or column vector.')

    # Section 2: Determine vector magnitude
    rows = len(vector); cols = len(vector[0])
    mag = 0
    for row in vector:
        for value in row:
            mag += value ** 2
    mag = mag ** 0.5

    # Section 3: Make a copy of vector
    new = copy_matrix(vector)

    # Section 4: Unitize the copied vector
    for i in range(rows):
        for j in range(cols):
            new[i][j] = new[i][j] / mag

    return new

Using Numpy For The Above Operations

How would we do all of these actions with numpy? It’s pretty simple and elegant. The code below follows the same order of functions we just covered above but shows how to do each one in numpy. The code below is in the file NumpyToolsPractice.py in the repo. Copy the code below or get it from the repo, but I strongly encourage you to run it and play with it.

import numpy as np


# Zeros Matrix
print(np.zeros((3, 3)), '\n')

# Identity Matrix
print(np.identity(3), '\n')

A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
C = np.array([[2, 4, 6], [8, 10, 12], [14, 16, 18]])

# Matrix Copy
AC = A.copy()
print(AC, '\n')

# Transpose a matrix
AT = np.transpose(A)
print(AT, '\n')

# Add and Subtract
SumAC = A + C
print(SumAC, '\n')

DifCA = C - A
print(DifCA, '\n')

# Matrix Multiply
ProdAC = np.matmul(A, C)
print(ProdAC, '\n')

# Multiply a List of Matrices
arr = [A, C, A, C, A, C]
Prod = np.matmul(A, C)
num = len(arr)
for i in range(2, num):
    Prod = np.matmul(Prod, arr[i])
print(Prod, '\n')

ChkP = np.matmul(
            np.matmul(
                np.matmul(
                    np.matmul(
                        np.matmul(arr[0], arr[1]),
                        arr[2]), arr[3]), arr[4]), arr[5])
print(ChkP, '\n')

# Check Equality of Matrices
print(Prod == ChkP, '\n')
# https://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html
print(np.allclose(Prod, ChkP), '\n')

# Dot Product (follows the same rules as matrix multiplication)
# https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html
v1 = np.array([[2, 4, 6]])
v2 = np.array([[1], [2], [3]])
ans1 = np.dot(v1, v2)
ans2 = np.dot(v1, v2)[0, 0]
print(f'ans1 = {ans1}, ans2 = {ans2}\n')

# Unitize an array
mag1 = (1*1 + 2*2 + 3*3) ** 0.5
mag2 = np.linalg.norm(v2)
norm1 = v2 / mag1
norm2 = v2 / mag2
print(f'mag1 = {mag1}, mag2 = {mag2}, they are equal: {mag1 == mag2}\n')
print(norm1, '\n')
print(norm2, '\n')
print(norm1 == norm2)

Closing

That’s it for now. This library will grow of course with each new post. I’ll introduce new helper functions if and when they are needed in future posts, and have separate posts for those additions that require more explanation. But these functions are the most basic ones. Some of these also support the work for the inverse matrix post and for the solving a system of equations post.