Matrix multiplication

The Vectorization chapter treated arrays as bags of independent numbers. This one breaks that: matrix multiplication mixes elements — every output is a dot product woven from a whole row and a whole column. It is the engine behind every neural-network layer, every linear-regression fit, every PCA decomposition. NumPy gives you three ways to spell it — np.dot, np.matmul, and the @ operator — and they do not quite agree at the edges, which is exactly where the bugs live.

The three ways

import numpy as np

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

# Three spellings, same result for 2D arrays
print("np.dot:\n",    np.dot(A, B))
print("np.matmul:\n", np.matmul(A, B))
print("A @ B:\n",     A @ B)

np.dot:
 [[19 22]
 [43 50]]
np.matmul:
 [[19 22]
 [43 50]]
A @ B:
 [[19 22]
 [43 50]]

For two 2-D matrices they are identical. The differences surface at the edges — 1-D vectors and 3-D+ tensors.

The shape rule

For A @ B to work, the last axis of A must equal the second-to-last axis of B. The result drops those two and keeps the rest.

Why? Each output element is a dot product (sum of element-wise products) between a row of A and a column of B. Rows and columns must be the same length — that is why the inner k must agree.

When you see ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, that is NumPy telling you the inner dimensions don’t line up. Print the shapes and read right-to-left.

Where dot, matmul, and @ diverge

np.dot has legacy behaviour — for 1-D arrays it does an inner product, for higher dims it sums over a single axis. np.matmul (and @) treat the last two axes as the matrix and broadcast everything else. That matters the moment you have a batch dimension.

import numpy as np

# Batch of 4 matrices, each 2x3
A = np.arange(24).reshape(4, 2, 3)
# Batch of 4 matrices, each 3x5
B = np.arange(60).reshape(4, 3, 5)

# matmul / @ broadcast the leading "batch" axis
C = A @ B
print("A @ B shape:", C.shape)   # (4, 2, 5)

# np.dot does NOT broadcast batches the same way
D = np.dot(A, B)
print("np.dot shape:", D.shape)  # (4, 2, 4, 5) — surprising!

A @ B shape: (4, 2, 5)
np.dot shape: (4, 2, 4, 5)

That second result — a 4-D (4, 2, 4, 5) — is almost never what you want: np.dot paired every batch matrix of A with every batch matrix of B. Use @ (or np.matmul). PEP 465 added @ in Python 3.5 precisely because the ambiguity around dot was costing people hours of debugging.

A neural net layer in one line

Every dense layer is output = inputs @ W + b. Here a batch of 64 feature vectors of size 784 flows through a layer that produces 128 hidden units:

import numpy as np

rng = np.random.default_rng(42)

# A batch of 64 flattened 28x28 images
X = rng.standard_normal((64, 784))

# A weight matrix mapping 784 -> 128
W = rng.standard_normal((784, 128)) * 0.01
b = np.zeros(128)

# Forward pass through one dense layer
H = X @ W + b
print("Input:", X.shape, "Weights:", W.shape, "Output:", H.shape)
# (64, 784) @ (784, 128) = (64, 128) — the 784s cancel

Input: (64, 784) Weights: (784, 128) Output: (64, 128)

The 784s on either side of @ are the “inner” dimensions, and they cancel. You are left with (64, 128) — 64 samples, each now represented by 128 hidden activations. That is it. That is every dense layer.

Other “products” you’ll see

NumPy has three more product flavours that come up regularly:

import numpy as np

u = np.array([1, 2, 3])
v = np.array([4, 5, 6])

# Inner (dot) product — scalar
print("inner:", np.inner(u, v))        # 1*4 + 2*5 + 3*6 = 32

# Outer product — matrix shape (3, 3)
print("outer:\n", np.outer(u, v))

# Cross product — only defined for 3D vectors
print("cross:", np.cross(u, v))        # vector perpendicular to both

inner: 32
outer:
 [[ 4  5  6]
 [ 8 10 12]
 [12 15 18]]
cross: [-3  6 -3]

np.outer is the move when you need a rank-1 matrix from two vectors — common in low-rank updates, word-embedding tricks, and deriving backprop by hand. np.cross shows up in graphics and physics.

Reading shape errors

Almost every matmul error is a mismatched inner dimension. The fix is one of:

1. Transpose one operand (A.T @ B instead of A @ B).
2. Reshape to add or move an axis (v[:, None] to make a column).
3. Check the layout — did your data loader produce (features, samples) when you expected (samples, features)?

When in doubt, print(A.shape, B.shape) before the @, and the mismatch is almost always obvious.

In one breath

Matrix multiplication mixes elements: each output is the dot product of a row of A and a column of B, so the inner axes must match — (m, k) @ (k, n) → (m, n), the k cancels. Three spellings agree on 2-D (A @ B == np.dot == np.matmul), but on 3-D+ they diverge: @/matmul treat the last two axes as the matrix and broadcast the batch ((4,2,3) @ (4,3,5) → (4,2,5)), while np.dot pairs every batch with every batch ((4,2,4,5)) — so prefer @. Every dense layer is X @ W + b ((64,784) @ (784,128) → (64,128)); attention is Q @ K.T. Also handy: np.inner (scalar), np.outer (rank-1 matrix), np.cross (3-D). Most matmul errors are a mismatched inner dim — print the shapes.

Practice

A question to carry forward

Matrix multiplication builds a linear transformation — it takes A and applies it to data. But half of linear algebra is the inverse question: given the transformation A and a result b, what input x produced it? Solve Ax = b. Or: what are A’s hidden axes — its eigenvalues, its SVD — the decompositions that PCA and a hundred other methods quietly run on?

Those answers don’t come from @; they come from np.linalg. The next lesson opens that toolbox. What do np.linalg.solve, inv, eig, and svd actually give you — and why, even when the textbook writes x = A⁻¹b, should you almost never literally invert a matrix in code?