Matrix Multiplication and Linear Transformations

What this page teaches:

import numpy as np

A = np.array([[1, 2, 3],
              [4, 5, 6]])   # shape (2, 3)

B = np.array([[7,  8],
              [9,  10],
              [11, 12]])    # shape (3, 2)

# Matrix multiplication — the @ operator is the modern way
result = A @ B
print(result)
# [[ 58,  64],
#  [139, 154]]

print(result.shape)  # (2, 2)

# np.matmul() is identical to @
result2 = np.matmul(A, B)
print(np.array_equal(result, result2))  # True

# np.dot() ALSO works for 2D arrays — but confusing for higher dims.
# Use @ everywhere. It is the clearest and most consistent.
result3 = np.dot(A, B)
print(np.array_equal(result, result3))  # True

# Verify manually: result[0,0] = dot(A[0], B[:,0])
manual = np.dot(A[0], B[:, 0])
print(manual)   # 58  ✓

# Valid multiplications
A = np.random.randn(2, 3)
B = np.random.randn(3, 4)
C = A @ B
print(C.shape)   # (2, 4)  ← outer dims: 2 and 4

# Real ML sizes
X = np.random.randn(1000, 4)    # 1000 training examples, 4 features
W = np.random.randn(4, 64)      # weight matrix for 64 neurons
out = X @ W
print(out.shape)   # (1000, 64) ← 1000 examples, each mapped to 64 values

# The most useful trick: check shapes BEFORE multiplying
def check_matmul(A, B):
    if A.shape[1] != B.shape[0]:
        print(f"Shape mismatch! A={A.shape}, B={B.shape}")
        print(f"A has {A.shape[1]} columns, B has {B.shape[0]} rows — must match.")
        print(f"Try: B = B.T  (current B.T shape would be {B.T.shape})")
    else:
        result = A @ B
        print(f"{A.shape} @ {B.shape} → {result.shape}  ✓")
    return None

check_matmul(np.random.randn(5, 3), np.random.randn(3, 7))   # (5,3) @ (3,7) → (5,7)  ✓
check_matmul(np.random.randn(5, 3), np.random.randn(7, 3))   # Shape mismatch!

The key insight: linear layers (matrix multiplications) give neural networks their capacity — the ability to transform data into useful representations. Non-linear activation functions give neural networks their power — the ability to learn patterns that no straight line could capture. You need both. One without the other is useless.

# Proving linearity: f(ax) = a*f(x)
W = np.array([[2, 1], [0, 3], [1, 2]])  # (3, 2) weight matrix

x = np.array([1.0, 2.0])   # input vector

# Transform the input
f_x = W @ x
print(f_x)   # [4., 6., 5.]

# Scale input by 3, then transform
f_3x = W @ (3 * x)
print(f_3x)  # [12., 18., 15.]

# 3 × transform = same result
print(3 * f_x)  # [12., 18., 15.]  ✓  same!

# Proving: f(x+y) = f(x) + f(y)
y = np.array([0.5, 1.0])
print(np.allclose(W @ (x + y),  (W @ x) + (W @ y)))  # True ✓

# Now prove non-linear activations break this linearity:
def relu(z): return np.maximum(0, z)

x1 = np.array([2.0, -1.0])
x2 = np.array([-3.0, 4.0])

# relu(x1 + x2) ≠ relu(x1) + relu(x2)
print(relu(x1 + x2))           # [0., 3.]
print(relu(x1) + relu(x2))     # [2., 3.]  ← DIFFERENT
# This non-linearity is exactly why neural networks can learn complex patterns

import numpy as np

np.random.seed(42)

# 32 DoorDash orders, 4 features each
X = np.random.randn(32, 4)

# ── Layer 1: 4 features → 64 neurons ──────────────────────────────────
# Weight matrix: random initialisation (will be updated during training)
# Shape: (input_features, output_neurons) = (4, 64)
W1 = np.random.randn(4, 64) * 0.01   # small init to avoid exploding activations
b1 = np.zeros(64)                     # bias initialised to zero

Z1 = X @ W1 + b1                      # (32,4) @ (4,64) + (64,) → (32,64)
A1 = np.maximum(0, Z1)                # ReLU activation: max(0, x)

print(f"Layer 1 output shape: {A1.shape}")   # (32, 64)

# ── Layer 2: 64 neurons → 32 neurons ─────────────────────────────────
W2 = np.random.randn(64, 32) * 0.01
b2 = np.zeros(32)

Z2 = A1 @ W2 + b2   # (32,64) @ (64,32) + (32,) → (32,32)
A2 = np.maximum(0, Z2)

print(f"Layer 2 output shape: {A2.shape}")   # (32, 32)

# ── Output layer: 32 neurons → 1 prediction ──────────────────────────
W3 = np.random.randn(32, 1) * 0.01
b3 = np.zeros(1)

Z3 = A2 @ W3 + b3   # (32,32) @ (32,1) + (1,) → (32,1)
# No activation on output for regression — the raw number IS the prediction

predictions = Z3.squeeze()   # (32,) — one predicted delivery time per order
print(f"Predictions shape: {predictions.shape}")  # (32,)
print(f"Sample prediction: {predictions[0]:.4f} minutes")

# The total parameter count in this tiny network:
total_params = (4*64 + 64) + (64*32 + 32) + (32*1 + 1)
print(f"Total learnable parameters: {total_params}")  # 2,369
# GPT-4 has ~1.8 trillion. Same structure, vastly larger.

# Architecture as shape decisions
# Narrow network: low capacity, fast, less likely to overfit
narrow_shapes = [(4, 16), (16, 8), (8, 1)]

# Wide network: high capacity, slow, more likely to overfit
wide_shapes = [(4, 512), (512, 256), (256, 1)]

# Deep network: many transformations
deep_shapes = [(4, 32), (32, 32), (32, 32), (32, 32), (32, 1)]

for shapes in [narrow_shapes, wide_shapes, deep_shapes]:
    params = sum(rows * cols + cols for rows, cols in shapes)
    print(f"Shapes: {shapes}")
    print(f"Total params: {params:,}
")

# Output:
# Narrow: 16*4+16 + 16*8+8 + 8*1+1 = 201 params
# Wide: 512*4+512 + 256*512+256 + 1*256+1 = 133,889 params
# Deep: each (32,32) layer: 32*32+32 = 1,056, total ~4,257 params
#
# Wide adds capacity via bigger matrices.
# Deep adds capacity via more transformations.
# Modern models (Transformers) do both.

# Batched matmul: apply matmul across a batch dimension

# Scenario: multi-head attention with 12 heads
# Each head works on queries and keys of shape (seq_len, d_k)
# We process all 12 heads at once with a 3D tensor

batch_size = 32          # 32 sentences in this batch
seq_len    = 128         # each sentence has 128 tokens
num_heads  = 12
d_k        = 64          # dimension per head

# Q and K are 3D: (batch, seq_len, d_k)
Q = np.random.randn(batch_size, seq_len, d_k)
K = np.random.randn(batch_size, seq_len, d_k)

# Attention scores: Q @ K^T for each item in the batch
# K needs to be transposed on the last two dims: (batch, d_k, seq_len)
K_transposed = K.transpose(0, 2, 1)   # (32, 64, 128)

# @ operator applies matmul to the last two dims, broadcast over batch
scores = Q @ K_transposed              # (32, 128, 64) @ (32, 64, 128) → (32, 128, 128)
print(scores.shape)   # (32, 128, 128)

# This is the attention score matrix:
# scores[i, t1, t2] = how much token t1 attends to token t2 in sentence i
# One (128, 128) attention map per sentence, 32 sentences at once

# Scale by sqrt(d_k) to prevent vanishing gradients in softmax
scores = scores / np.sqrt(d_k)

# Softmax to get attention weights (probabilities that sum to 1)
def softmax(x):
    e = np.exp(x - x.max(axis=-1, keepdims=True))  # stable softmax
    return e / e.sum(axis=-1, keepdims=True)

weights = softmax(scores)   # (32, 128, 128) — attention weights
print(f"Attention weights shape: {weights.shape}")
print(f"Row sums (should all be 1.0): {weights[0, 0].sum():.4f}")

# SVD in NumPy
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9],
              [10, 11, 12]], dtype=float)   # (4, 3)

U, sigma, Vt = np.linalg.svd(A, full_matrices=False)

print(f"U shape:     {U.shape}")       # (4, 3)
print(f"sigma shape: {sigma.shape}")   # (3,)  ← the singular values
print(f"Vt shape:    {Vt.shape}")      # (3, 3)
print(f"Singular values: {sigma.round(2)}")
# [2.54e+01, 1.72e+00, 1.41e-15]
# The third singular value is essentially 0 — this matrix has rank 2!
# It means all information lives in just 2 directions.

# Reconstruct A from SVD (verify it works)
A_reconstructed = U @ np.diag(sigma) @ Vt
print(f"Reconstruction error: {np.abs(A - A_reconstructed).max():.2e}")  # ~1e-14

# ── Low-rank approximation — the key application ─────────────────────
# Keep only the top k singular values, set the rest to 0
# This compresses the matrix while preserving the most important structure

k = 1   # keep only the largest singular value
A_compressed = (
    sigma[0] * np.outer(U[:, 0], Vt[0, :])
)
print(f"Compression error with k=1: {np.abs(A - A_compressed).max():.4f}")
# Very small — because A was already nearly rank-1

# ── Why this matters for ML ────────────────────────────────────────────
# PCA uses SVD to find the k most important directions in your data
# Recommender systems decompose the user-item matrix with SVD
# Weight matrices in neural nets are initialised using SVD principles
# Attention matrices can be compressed using low-rank approximation

# Check: sigma values tell you importance of each direction
total_variance = sigma.sum()
for i, s in enumerate(sigma):
    print(f"Direction {i+1}: {100 * s / total_variance:.1f}% of total variance")

# PCA IS SVD — shown directly
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

np.random.seed(42)
X = np.random.randn(100, 4)   # 100 orders, 4 features

# Standardise first (PCA is sensitive to scale)
X_scaled = StandardScaler().fit_transform(X)

# Method 1: sklearn PCA
pca = PCA(n_components=2)
pca.fit(X_scaled)
print("PCA components (sklearn):")
print(pca.components_.round(3))

# Method 2: manual SVD
X_mean = X_scaled - X_scaled.mean(axis=0)
U, sigma, Vt = np.linalg.svd(X_mean, full_matrices=False)
print("
PCA components (manual SVD — Vt):")
print(Vt[:2].round(3))

# The rows of Vt ARE the principal components.
# Small differences due to sign conventions — the magnitude is identical.
# sklearn PCA = SVD. No magic.

# The debugging template — use this every time shapes behave unexpectedly
def matmul_debug(A, B, name=""):
    print(f"{'─'*40}")
    if name: print(f"Operation: {name}")
    print(f"  A shape: {A.shape}  dtype: {A.dtype}")
    print(f"  B shape: {B.shape}  dtype: {B.dtype}")

    if A.ndim >= 2 and B.ndim >= 2:
        if A.shape[-1] != B.shape[-2]:
            print(f"  ✗ SHAPE MISMATCH: A cols ({A.shape[-1]}) ≠ B rows ({B.shape[-2]})")
            print(f"  Try: B = B.T  → B.T shape would be {B.T.shape}")
            return None

    result = A @ B
    print(f"  ✓ Result shape: {result.shape}")
    return result

# Usage
X = np.random.randn(32, 4)
W = np.random.randn(4, 64)
out = matmul_debug(X, W, "Input layer")   # ✓ Result shape: (32, 64)

W_wrong = np.random.randn(64, 4)
out = matmul_debug(X, W_wrong, "Wrong weight shape")  # ✗ SHAPE MISMATCH

import time

# Compare: Python loop vs NumPy vs (conceptually) GPU

n = 1000
A = np.random.randn(n, n)
B = np.random.randn(n, n)

# Method 1: Pure Python loop — NEVER DO THIS
def naive_matmul(A, B):
    m, k = A.shape
    _, n = B.shape
    result = np.zeros((m, n))
    for i in range(m):
        for j in range(n):
            for p in range(k):
                result[i, j] += A[i, p] * B[p, j]
    return result

# Method 2: NumPy — uses BLAS, highly optimised C under the hood
t0 = time.time()
result_numpy = A @ B
numpy_time = time.time() - t0
print(f"NumPy:  {numpy_time*1000:.1f}ms")

# For n=100 the naive loop takes ~10 seconds.
# NumPy takes ~0.1ms. That's a 100,000× speedup.
# A GPU would take ~0.001ms — another 100× faster.

# The rule: NEVER loop over matrix elements in Python.
# Always express operations as NumPy vectorised operations.
# The difference is not 2× or 10× — it is 100,000×.

# Vectorisation means: express the operation as array operations
# so the CPU/GPU can execute them in parallel using SIMD instructions.

# Wrong (loop)
manual_row_sums = np.array([A[i].sum() for i in range(len(A))])

# Right (vectorised)
vectorised_row_sums = A.sum(axis=1)

print(np.allclose(manual_row_sums, vectorised_row_sums))  # True
# Both give the same answer. The vectorised version is ~1000× faster.

The practical rule: if you find yourself writing a for loop over rows or columns of a NumPy array, stop. There is almost certainly a vectorised operation that does the same thing. The entire skill of "writing efficient ML code" is learning to express computations as matrix operations instead of loops. This track shows you how to do that at every step.

# The estimate an engineer runs before a design review -- not after

def layer_cost(batch, in_dim, out_dim):
    """FLOPs and parameter count for one X @ W + b layer."""
    flops = 2 * batch * in_dim * out_dim   # the matmul itself
    params = in_dim * out_dim + out_dim    # W and b
    return flops, params

# Current production layer: 1000-example batch, 256 -> 256
current_flops, current_params = layer_cost(1000, 256, 256)
print(f"Current: {current_flops:,} FLOPs, {current_params:,} params")

# Proposed change: double the hidden dimension to 512
proposed_flops, proposed_params = layer_cost(1000, 256, 512)
print(f"Proposed: {proposed_flops:,} FLOPs, {proposed_params:,} params")

print(f"Compute increase: {proposed_flops / current_flops:.1f}x")
print(f"Parameter increase: {proposed_params / current_params:.1f}x")

# Output:
# Current: 131,328,000 FLOPs, 65,792 params
# Proposed: 262,656,000 FLOPs, 131,584 params
# Compute increase: 2.0x
# Parameter increase: 2.0x
#
# This is the actual conversation in the design review:
# "Doubling the hidden dimension doubles both serving cost and
#  model size -- is the accuracy gain worth 2x the inference bill?"

Where SVD shows up in this same conversation: when a proposed layer is too expensive as a full (n × p) matrix, teams often reach for a low-rank approximation instead — this is the entire idea behind LoRA (Low-Rank Adaptation), the dominant technique for cheaply fine-tuning large language models. Instead of updating the full weight matrix, LoRA learns two much smaller matrices whose product approximates the update, exactly the U-Σ-Vᵀ decomposition from the SVD section above, applied to make fine-tuning a multi-billion parameter model affordable on a single GPU.

That operation exists and is useful, but it is called element-wise multiplication, and it is a completely different result from matrix multiplication. Matrix multiplication takes a row from the left matrix and a column from the right matrix and computes their dot product for each position in the result — as walked through step by step earlier on this page. The two operations even use different NumPy syntax on purpose: @ or np.matmul() for real matrix multiplication, plain * for element-wise. Mixing them up does not raise an error when shapes happen to allow both — it just silently produces the wrong numbers.

Two same-shaped square matrices are exactly the dangerous case, because both A @ B and A * B run without error and return a same-shaped result — just a numerically different one. There is no shape mismatch to warn you. This is precisely the bug in the 'Result has unexpected shape' entry in the errors section above: the code runs fine, the shapes look fine, and the output is simply the wrong computation, often not discovered until model accuracy quietly underperforms for reasons nobody can explain.

Matrix multiplication is not commutative, and the callout above states this directly: A @ B and B @ A are generally different, and one of the two often is not even a valid operation at all once the shapes stop lining up in reverse. The reason is structural, not a quirk of notation — the left matrix's rows drive the result's rows and the right matrix's columns drive the result's columns, so swapping the order changes which vectors get dotted with which. Order is not a stylistic choice in matrix multiplication; it changes the answer.

The 'Matrix shapes define the network architecture' section above makes the reverse point directly: when a paper describes 'a hidden dimension of 512,' that sentence is entirely describing the number of columns in a weight matrix. Narrow versus wide versus deep are not abstract design philosophies competing with the shapes — they are literally different choices of matrix shape, and the resulting parameter count and compute cost follow mechanically from that choice, as the FLOPs estimate above shows.

GPUs move the practical ceiling enormously higher, but they do not remove it — the efficiency section above shows the same operation still taking measurably longer as matrix sizes grow, and the FLOPs formula (2 × m × n × p) still applies whether the hardware is a laptop CPU or a data-centre GPU cluster. This is exactly why the industry spends serious engineering effort on techniques like low-rank approximation (SVD) and LoRA rather than simply 'throwing more GPU at it' — at the scale of models like GPT-4, the multiplication cost is large enough that shrinking it meaningfully changes what is financially feasible to train or serve.

To compute result[i, j], take row i of the left matrix and column j of the right matrix, multiply them element by element, and sum the products — a single dot product. Repeat that for every (i, j) position and you have the full result. The shape rule that makes this possible: for an (m × n) matrix times an (n × p) matrix, the inner dimensions — both n — must match, since a row and a column can only be dotted if they have the same length. The result takes the outer dimensions, (m × p).

It is not commutative because the operation is structurally asymmetric: the left matrix's rows and the right matrix's columns get dotted together, so swapping the order changes which vectors are being combined, and often changes whether the shapes are even compatible at all. It matters constantly in practice — a layer computes X @ W, not W @ X, because X's rows are examples and W's columns are output neurons; reversing the order would either throw a shape error or silently compute something with a completely different meaning.

The standard estimate is 2 × m × n × p FLOPs for an (m × n) times (n × p) multiplication — two operations, a multiply and an add, per term in each output's dot product. In practice this is exactly what gets computed before approving an architecture change: doubling a hidden dimension roughly doubles both the FLOPs and the parameter count of that layer, which translates directly into serving latency and GPU cost at inference time. Teams run this arithmetic before committing to a more expensive architecture, not after.

A fully-connected layer is Z = X @ W + b followed by an activation function, and nothing more is hidden inside it. If X is a batch of examples with shape (batch, in_features) and W is the layer's weight matrix with shape (in_features, out_features), the matmul produces Z with shape (batch, out_features) — every example transformed into a new representation simultaneously, in one operation. The activation function is then applied element-wise to Z without changing its shape at all.

Singular Value Decomposition breaks any matrix into three simpler matrices, U, Σ, and Vᵀ, equivalent to a rotation, a scaling, and another rotation — and the singular values in Σ rank how much information each direction of the data carries. A concrete production application is LoRA (Low-Rank Adaptation), the standard technique for cheaply fine-tuning large language models: instead of updating a full, enormous weight matrix, LoRA learns two much smaller matrices whose product approximates the needed update, using far fewer parameters — the same low-rank idea that underlies SVD-based compression and PCA.