Eric Campbell 6 min read Statistics

Linear Algebra: Matrix Multiplication Explained

Matrix multiplication isn't just academic exercise—it's the workhorse of modern computing. Every time you use a recommendation system, apply a filter to an image, or run a neural network, matrix...

Key Insights

  • Matrix multiplication is the fundamental operation behind neural networks, computer graphics, and data transformations—understanding it deeply separates competent developers from exceptional ones
  • The dimension compatibility rule (m×n × n×p = m×p) is non-negotiable: the inner dimensions must match, and violating this causes runtime errors that waste debugging time
  • NumPy’s optimized implementations are 100-1000x faster than pure Python loops, making the choice between @ operator and manual implementation a performance-critical decision

Introduction: Why Matrix Multiplication Matters

Matrix multiplication isn’t just academic exercise—it’s the workhorse of modern computing. Every time you use a recommendation system, apply a filter to an image, or run a neural network, matrix multiplication is happening under the hood. In deep learning, a single forward pass through a neural network involves hundreds or thousands of matrix multiplications.

The difference between understanding matrix multiplication conceptually versus mechanically is substantial. You need to know not just what it does, but how it works, why dimensions matter, and when to optimize. This article strips away the mathematical formalism and focuses on practical implementation.

Matrix Multiplication Basics

Matrix multiplication combines two matrices to produce a third. Given matrix A with dimensions m×n and matrix B with dimensions n×p, their product AB yields a matrix C with dimensions m×p.

The critical rule: the number of columns in the first matrix must equal the number of rows in the second matrix. This isn’t a suggestion—it’s a mathematical requirement. A 3×2 matrix can multiply with a 2×4 matrix (result: 3×4), but not with a 3×4 matrix.

Here’s the mechanical process: each element C[i][j] is computed by taking the dot product of row i from matrix A and column j from matrix B.

def matrix_multiply_basic(A, B):
    """
    Multiply two matrices using the fundamental definition.
    A: m x n matrix
    B: n x p matrix
    Returns: m x p matrix
    """
    m, n = len(A), len(A[0])
    n2, p = len(B), len(B[0])
    
    if n != n2:
        raise ValueError(f"Incompatible dimensions: {n} != {n2}")
    
    # Initialize result matrix with zeros
    C = [[0 for _ in range(p)] for _ in range(m)]
    
    # Compute each element
    for i in range(m):
        for j in range(p):
            for k in range(n):
                C[i][j] += A[i][k] * B[k][j]
    
    return C

# Example: 2x2 matrices
A = [[1, 2],
     [3, 4]]

B = [[5, 6],
     [7, 8]]

result = matrix_multiply_basic(A, B)
print(result)  # [[19, 22], [43, 50]]

Let’s verify one element manually: C[0][0] = (1×5) + (2×7) = 5 + 14 = 19. This is the dot product of A’s first row [1, 2] and B’s first column [5, 7].

The Mechanics: Row-by-Column Operation

The triple-nested loop structure reveals the operation’s true nature. The outer two loops (i and j) iterate through the result matrix positions. The inner loop (k) computes the dot product.

Think of it as a systematic scanning process:

  1. Fix a row from the first matrix
  2. Fix a column from the second matrix
  3. Multiply corresponding elements and sum them
  4. Store the result
  5. Move to the next column, then the next row
def matrix_multiply_annotated(A, B):
    """
    Matrix multiplication with detailed logging to show the process.
    """
    m, n = len(A), len(A[0])
    n2, p = len(B), len(B[0])
    
    if n != n2:
        raise ValueError(f"Cannot multiply {m}x{n} with {n2}x{p}")
    
    C = [[0] * p for _ in range(m)]
    
    for i in range(m):
        for j in range(p):
            # Compute C[i][j]
            dot_product = 0
            terms = []
            
            for k in range(n):
                product = A[i][k] * B[k][j]
                dot_product += product
                terms.append(f"{A[i][k]}*{B[k][j]}")
            
            C[i][j] = dot_product
            print(f"C[{i}][{j}] = {' + '.join(terms)} = {dot_product}")
    
    return C

# Example with 2x3 and 3x2 matrices
A = [[1, 2, 3],
     [4, 5, 6]]

B = [[7, 8],
     [9, 10],
     [11, 12]]

result = matrix_multiply_annotated(A, B)
# Output shows: C[0][0] = 1*7 + 2*9 + 3*11 = 58, etc.

Common mistakes: attempting to multiply incompatible dimensions, confusing row and column indexing, or implementing element-wise multiplication instead of matrix multiplication (that’s the Hadamard product, a different operation entirely).

Practical Implementation

Pure Python loops are pedagogically useful but practically slow. NumPy provides highly optimized implementations that leverage BLAS (Basic Linear Algebra Subprograms) libraries written in C and Fortran.

import numpy as np
import time

# Create larger test matrices
size = 500
A_np = np.random.rand(size, size)
B_np = np.random.rand(size, size)

# NumPy's @ operator (Python 3.5+)
start = time.time()
C1 = A_np @ B_np
numpy_time = time.time() - start

# np.matmul() - equivalent to @
start = time.time()
C2 = np.matmul(A_np, B_np)
matmul_time = time.time() - start

# np.dot() - works but @ is preferred for matrices
start = time.time()
C3 = np.dot(A_np, B_np)
dot_time = time.time() - start

print(f"@ operator: {numpy_time:.4f}s")
print(f"np.matmul: {matmul_time:.4f}s")
print(f"np.dot: {dot_time:.4f}s")
print(f"Results equal: {np.allclose(C1, C2) and np.allclose(C2, C3)}")

# Pure Python comparison (smaller size for sanity)
A_list = [[float(x) for x in row] for row in np.random.rand(100, 100)]
B_list = [[float(x) for x in row] for row in np.random.rand(100, 100)]

start = time.time()
C_python = matrix_multiply_basic(A_list, B_list)
python_time = time.time() - start

print(f"\nPure Python (100x100): {python_time:.4f}s")
print(f"NumPy speedup factor: ~{python_time / numpy_time * (size/100)**3:.0f}x")

On my machine, NumPy is approximately 500-1000x faster for 500×500 matrices. The @ operator is clearest and should be your default choice. Use np.matmul() when you need explicit function calls (e.g., in functional programming contexts). The np.dot() function has quirky behavior with higher-dimensional arrays, so avoid it for matrix multiplication.

Real-World Application: Linear Transformations

Matrix multiplication represents linear transformations. Want to rotate a point? Multiply by a rotation matrix. Scale coordinates? Multiplication. This is how graphics engines transform vertices.

import numpy as np
import matplotlib.pyplot as plt

def create_rotation_matrix(angle_degrees):
    """Create a 2D rotation matrix."""
    theta = np.radians(angle_degrees)
    return np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta), np.cos(theta)]
    ])

# Create a simple shape (a square)
square = np.array([
    [0, 1, 1, 0, 0],  # x coordinates
    [0, 0, 1, 1, 0]   # y coordinates
])

# Rotate by 45 degrees
rotation_matrix = create_rotation_matrix(45)
rotated_square = rotation_matrix @ square

# Visualize
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

ax1.plot(square[0], square[1], 'b-o', linewidth=2)
ax1.set_title('Original Square')
ax1.grid(True)
ax1.axis('equal')

ax2.plot(rotated_square[0], rotated_square[1], 'r-o', linewidth=2)
ax2.set_title('Rotated 45°')
ax2.grid(True)
ax2.axis('equal')

plt.tight_layout()
plt.savefig('rotation_example.png', dpi=150)
print("Transformation applied successfully")

This is exactly how game engines rotate objects, how computer vision algorithms transform images, and how robotics systems convert between coordinate frames. The matrix encodes the transformation; multiplication applies it.

Performance and Optimization

The naive triple-nested loop has O(n³) time complexity for n×n matrices. This scales poorly: doubling matrix size increases computation time by 8x. For a 1000×1000 matrix, that’s one billion multiplications.

Strassen’s algorithm reduces this to approximately O(n^2.807) through clever algebraic tricks, though it’s rarely implemented manually due to complexity and overhead for smaller matrices. Modern BLAS libraries use cache-optimized algorithms and SIMD instructions.

For sparse matrices (mostly zeros), specialized data structures save memory and computation:

import numpy as np
from scipy.sparse import csr_matrix
import time

# Create a sparse matrix (95% zeros)
size = 1000
density = 0.05

# Dense representation
A_dense = np.random.rand(size, size)
A_dense[A_dense > density] = 0
B_dense = np.random.rand(size, size)
B_dense[B_dense > density] = 0

# Sparse representation
A_sparse = csr_matrix(A_dense)
B_sparse = csr_matrix(B_dense)

# Compare memory usage
print(f"Dense memory: {A_dense.nbytes + B_dense.nbytes:,} bytes")
print(f"Sparse memory: {A_sparse.data.nbytes + B_sparse.data.nbytes:,} bytes")

# Compare computation time
start = time.time()
C_dense = A_dense @ B_dense
dense_time = time.time() - start

start = time.time()
C_sparse = A_sparse @ B_sparse
sparse_time = time.time() - start

print(f"\nDense computation: {dense_time:.4f}s")
print(f"Sparse computation: {sparse_time:.4f}s")
print(f"Speedup: {dense_time/sparse_time:.2f}x")

For matrices with >90% sparsity, sparse representations can be 10-100x faster and use 10-50x less memory. This matters enormously in applications like graph algorithms, natural language processing, and recommendation systems where sparse matrices are common.

Conclusion and Next Steps

Matrix multiplication is deceptively simple in concept but rich in implications. You’ve now seen the mechanical process, understood dimension compatibility, and learned when to optimize. The key takeaways: always use NumPy for production code, understand the O(n³) complexity for capacity planning, and recognize sparse matrices as optimization opportunities.

From here, explore tensor operations (the higher-dimensional generalization), GPU acceleration with libraries like CuPy or JAX, and specialized applications like the attention mechanism in transformers (which is fundamentally matrix multiplication).

The mathematical elegance of matrix multiplication belies its computational power. Master it, and you’ve mastered a fundamental building block of modern computing.

Liked this? There's more.

Every week: one practical technique, explained simply, with code you can use immediately.