Matrix Multiplication

Composing two linear maps into one: the row-by-column rule that underlies nearly all numerical computing.

The rule

The product C = AB of an m-by-k matrix A and a k-by-n matrix B is an m-by-n matrix whose entry C[i,j] is the dot product of row i of A with column j of B. Explicitly, C[i,j] = sum over p of A[i,p] * B[p,j]. The inner dimension k must match; the outer dimensions m and n become the shape of the result.

Matrix multiplication is not commutative: in general AB is not equal to BA, and one product may be defined while the other is not. It is associative, (AB)C = A(BC), and distributive over addition, A(B+C) = AB + AC. These properties let you rearrange long products freely as long as you preserve order.

Three ways to see it

First, as row-by-column dot products, the definition above. Second, as a linear combination of columns: the j-th column of AB is A times the j-th column of B. Third, as composition of transformations: if A and B represent linear maps, then AB represents doing B first and then A. This last view is the deepest, because it explains why order matters and why associativity holds automatically.

Cost and performance

The naive algorithm for multiplying two n-by-n matrices costs on the order of n^3 scalar multiplications. Because this operation dominates so much scientific and machine-learning computation, highly tuned libraries (BLAS level 3) exploit cache blocking and vector hardware to reach near-peak throughput. Sub-cubic algorithms such as Strassen's exist but are used mainly for very large dense problems.

import numpy as np
A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 1], [1, 0]])
print(A @ B)   # [[2 1] [4 3]]
print(B @ A)   # [[3 4] [1 2]]  -- different

In large-scale plasma and structural simulations, each timestep applies operator matrices to state vectors; the aggregate cost of these products often sets the runtime of the whole model, so choosing algorithms that respect sparsity and structure is essential.