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 Swiggy 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.