Matrix Multiplication — Silent Shape Mismatch in Production

Matrices are the fundamental data structure of numerical computing. Every machine learning model is a stack of matrix multiplications. Every 3D graphics transformation is a matrix product. Every linear system of equations is a matrix equation. NumPy's core operation (np.dot, @) and every GPU-accelerated deep learning framework (PyTorch, TensorFlow) are built around efficient matrix multiplication.

Understanding matrix operations from first principles — not just API calls — means understanding why matrix multiplication is O(n³), why it is not commutative (AB ≠ BA in general), and why the memory layout matters as much as the algorithm for performance.

Core Operations

Matrix operations are the building blocks of linear algebra. Here we cover the three fundamental operations: addition, multiplication, and transpose. Each has specific shape requirements and complexity characteristics.

Addition – Element-wise, requires same dimensions (m×n). Complexity O(mn). The result is simply C[i][j] = A[i][j] + B[i][j]. In NumPy: A + B.

Multiplication – The dot product of rows from the first matrix with columns from the second. If A is m×n and B is n×p, result C is m×p. Complexity O(mnp). In NumPy: A @ B or np.dot(A, B). This is NOT element-wise — the inner dimension (n) must match.

Transpose – Flips rows and columns. If A is m×n, A^T is n×m. Complexity O(mn), memory-bound. In NumPy: A.T or np.transpose(A).

import numpy as np

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

# Addition: element-wise, O(nm)
print('A + B =\n', A + B)

# Multiplication: (m×n) @ (n×p) = (m×p), O(mnp)
# AB[i][j] = dot product of row i of A and col j of B
print('A @ B =\n', A @ B)

# Non-commutativity: AB ≠ BA in general
print('BA =\n', B @ A)
print('AB == BA:', np.allclose(A@B, B@A))  # False

# Transpose: rows become columns
print('A^T =\n', A.T)

# Inverse: A @ A^-1 = I
if np.linalg.det(A) != 0:
    print('A^-1 =\n', np.linalg.inv(A))

Naive Matrix Multiplication Implementation

The textbook algorithm for matrix multiplication uses three nested loops. It's the reference implementation — correct but wildly impractical for large matrices. The complexity is O(n³) for square matrices of size n.

Why is it O(n³)? Because for each element in the result (there are n²), you compute a dot product of length n, giving n³ operations. Real-world implementations use block decomposition and parallelisation to reduce the constant factor drastically.

This implementation is useful for understanding the mechanics but should never be used in production.

def matmul(A, B):
    """O(n^3) naive matrix multiplication."""
    rows_A, cols_A = len(A), len(A[0])
    rows_B, cols_B = len(B), len(B[0])
    assert cols_A == rows_B, 'Incompatible dimensions'
    C = [[0]*cols_B for _ in range(rows_A)]
    for i in range(rows_A):
        for j in range(cols_B):
            for k in range(cols_A):
                C[i][j] += A[i][k] * B[k][j]
    return C

A = [[1,2],[3,4]]
B = [[5,6],[7,8]]
print(matmul(A, B))  # [[19,22],[43,50]]

Why Cache Layout Dominates Performance

Matrix multiplication in Python (pure loops) for n=500 takes ~10 seconds. NumPy for the same: ~1ms. The 10,000x difference is not the algorithm — it is memory access patterns and hardware optimisation.

NumPy uses BLAS (Basic Linear Algebra Subprograms) with: block decomposition for cache efficiency (fitting submatrices in L1/L2 cache), SIMD vectorisation (processing 8 floats per instruction with AVX2), and multi-threading. The algorithm is still O(n³) — but the constant factor is ~1000x smaller than naive Python loops.

Rule: never implement matrix multiplication yourself. Use NumPy (@) or torch.matmul for GPU. The gap between a correct naive implementation and a production one is orders of magnitude.

The Mathematics of Matrix Multiplication: Dot Product View

Every element of the product matrix C[i][j] is the dot product of the i-th row of A and the j-th column of B. This view helps understand why the inner dimensions must match: you need the same number of entries in the row and column to compute the sum of products.

Formally: C[i][j] = Σ_{k=1 to n} A[i][k] * B[k][j]

This is the definition used in the naive implementation. From this, two important properties follow:

Associativity: (AB)C = A(BC). This is why matrix chain multiplication can be optimised — you can parenthesise differently.

Non-commutativity: AB ≠ BA in general. Only when both matrices are square and their product commutes (rare).

Transpose and Symmetric Matrices

Transposing a matrix swaps its rows and columns. The transpose of an m×n matrix is n×m. Important properties:

• (A^T)^T = A
• (A + B)^T = A^T + B^T
• (AB)^T = B^T A^T (note the reversed order)
• For a symmetric matrix: A = A^T. This means A is square and a[i][j] = a[j][i].

Symmetric matrices are common in machine learning (covariance matrices, kernel matrices). They have nice numerical properties: real eigenvalues, orthogonal eigenvectors.

In NumPy, A.T returns a view (no copy). To force a copy: A.T.copy().

import numpy as np
A = np.array([[1,2],[2,1]])  # symmetric
print('A.T:\n', A.T)
print('Is symmetric:', np.allclose(A, A.T))

# Transpose of product: (AB)^T = B^T A^T
B = np.array([[3,4],[5,6]])
lhs = (A @ B).T
rhs = B.T @ A.T
print('(AB)^T == B^T A^T:', np.allclose(lhs, rhs))

Matrix Inverse and Solving Linear Systems

The inverse of a square matrix A, denoted A^{-1}, satisfies A A^{-1} = I. Not all matrices have an inverse — only those with non-zero determinant (full rank).

When to invert? Almost never. For solving Ax = b, you should use np.linalg.solve(A, b) instead of np.linalg.inv(A) @ b. Why? Because solve uses LU decomposition, which is faster and more numerically stable than computing the full inverse.

Inverse existence: Check the determinant: np.linalg.det(A) != 0. But float precision can make this tricky — prefer checking the condition number: np.linalg.cond(A) > 1e12 indicates near-singular.

import numpy as np
A = np.array([[1,2],[3,4]])
b = np.array([5,6])

# Correct way
x = np.linalg.solve(A, b)
print('solve:', x)

# Incorrect way (slower, less stable)
x_bad = np.linalg.inv(A) @ b
print('inv @ b:', x_bad)

# Check residual
print('Residual:', np.linalg.norm(A @ x - b))

Numerical Stability: When Floating Point Breaks Your Matrix Math

Matrices in computers are stored as floating point numbers (float32, float64). This means every operation introduces rounding errors. In large matrix multiplications, these errors accumulate and can produce wildly wrong results.

Common problems: - Catastrophic cancellation: subtracting two nearly equal numbers causes loss of significant digits. Happens when computing (1 - cos(x)) for small x. - Overflow/underflow: intermediate products can exceed float range. Use np.clip or higher precision (float64). - Ill-conditioned matrices: small changes in input cause large changes in output. Detected via condition number.

Best practices: - Use float64 by default (NumPy default). float32 can be 2x faster on GPU but loses accuracy. - Normalise matrices before multiplication to avoid overflow. - Use np.linalg.cond to check conditioning of matrices you invert.