Building Matrix Multiplication in Python

Introduction

In our last lesson, we developed a powerful geometric intuition for matrices as operators that encode linear transformations. We saw how matrix-vector multiplication applies a transformation, and how matrix-matrix multiplication corresponds to composing successive transformations. The computational rules, while abstract, were shown to be a direct consequence of this geometric view.

Today, we transition from the "why" to the "how." Our goal is to implement matrix-vector and matrix-matrix multiplication from scratch in Python. This exercise is fundamental. By building these operations from first principles, without relying on the optimized routines of libraries like NumPy, you will solidify your understanding of the underlying mechanics. This is a crucial step for anyone interested in systems programming and performance, as it reveals precisely what high-performance libraries are abstracting away.

Representing Matrices and Vectors in Python

Before we can multiply, we need a way to store our mathematical objects. In Python, the most straightforward way to represent a vector is with a list of numbers, and a matrix with a list of lists, where each inner list represents a row.

For example, the vector $\mathbf{v}=\left[\begin{array}{c}2\\ 5\end{array}\right]$ and the matrix $\mathbf{M}=\left[\begin{array}{cc}1& 3\\ 4& 0\end{array}\right]$ would be represented as:

v = [2, 5]
M = [ [1, 3], 
      [4, 0] ]

The number of rows in M is len(M), and the number of columns is len(M[0]).

While this is our "from scratch" approach, it's important to know that the industry standard is to use NumPy arrays. A key resource for getting comfortable with Python for linear algebra demonstrates how to define vectors and matrices, first as lists, and then using NumPy.

The Logic of Multiplication: The Row-Column Rule

As we established, the entry in the $i$-th row and $j$-th column of a product matrix $\mathbf{C}=\mathbf{AB}$ is the dot product of the $i$-th row of $\mathbf{A}$ and the $j$-th column of $\mathbf{B}$.

Let's break this down algorithmically. Consider multiplying a matrix $\mathbf{A}$ of size $m\times n$ by a matrix $\mathbf{B}$ of size $n\times p$.

First, a question to set up our thinking: To calculate every element in the resulting $m\times p$ matrix $\mathbf{C}$, we need to iterate through its rows and columns. For each element $C_{ij}$, we then compute a dot product. How many nested loops do you think are required to implement this entire process? What is each loop responsible for?

...

The answer is three. We need:

1. A loop to iterate over the rows of the resulting matrix $\mathbf{C}$ (from $0$ to $m-1$). Let's use i.
2. A loop to iterate over the columns of the resulting matrix $\mathbf{C}$ (from $0$ to $p-1$). Let's use j.
3. An inner loop to compute the dot product for the element $C_{ij}$. This loop sums the products of elements from row i of $\mathbf{A}$ and column j of $\mathbf{B}$. It runs over the shared dimension, $n$. Let's use k.

The following video provides an excellent walkthrough of this three-loop logic.

From Logic to Code: A Python Implementation

Now that we have the algorithm, we can translate it into a Python function. Our function should take two matrices (as lists of lists), check if their dimensions are compatible for multiplication, and if so, compute and return the resulting matrix.

The process will be:

1. Extract the dimensions of the input matrices $\mathbf{A}$ and $\mathbf{B}$.
2. Assert that the number of columns in $\mathbf{A}$ equals the number of rows in $\mathbf{B}$.
3. Create a new matrix $\mathbf{C}$ (e.g., a list of lists) of the correct output dimensions, initialized with zeros.
4. Implement the three nested for loops (i, j, k) to populate $\mathbf{C}$. The core calculation will be C[i][j] += A[i][k] * B[k][j].

The same video resource walks through a complete, line-by-line implementation.

For a concise textual reference of the final algorithm, the following article provides a very clean implementation.

The Performance Reality: Why We Use Optimized Libraries

You have now built a working matrix multiplication function from scratch. This is a fantastic exercise for understanding the mechanics. However, in any practical application, especially in performance-sensitive fields like graphics or high-frequency trading, you would never use this implementation.

Why? Python for loops are notoriously slow. The Python interpreter has a high overhead for each iteration. The "from scratch" implementation becomes computationally infeasible for even moderately sized matrices.

Libraries like NumPy, PyTorch, and TensorFlow delegate these operations to underlying code written in C, C++, or Fortran. These low-level languages can execute the same loops orders of magnitude faster and can leverage hardware-specific optimizations (like SIMD instructions) that operate on multiple numbers at once.

The article you just read provides a stark comparison.

In NumPy, matrix multiplication is a single, highly optimized operation. You can use the function np.matmul() or, more commonly, the dedicated matrix multiplication operator @.

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])




# Using the @ operator for matrix multiplication
C = A @ B 

print(C)



# [[19 22]
#  [43 50]]

This single line replaces our entire matmul function and runs dramatically faster.

Conclusion

In this lesson, we bridged the gap between the theory and practice of matrix multiplication. You have seen how to translate the formal row-column rule into a working Python algorithm from scratch, giving you a first-principles understanding of the computation involved.

Our key takeaways are:

• Matrices can be represented as lists of lists in Python, allowing us to manipulate them with basic language features.
• Matrix-matrix multiplication is algorithmically a triply nested loop that iterates through the output matrix's rows and columns, with the innermost loop computing a dot product.
• Implementing algorithms from scratch is invaluable for understanding but is not practical for performance.
• High-performance libraries like NumPy are essential because they execute these fundamental operations in highly optimized, low-level code, offering massive speed advantages.

In our next lesson, we will step back from matrix operations to look more closely at the properties of the vectors themselves. We will introduce the crucial concepts of linear independence, span, and basis. These ideas form the vocabulary for describing vector spaces and are essential for understanding the rank of a matrix, the nature of solutions to linear systems, and the process of dimensionality reduction.