Attention from scratch: NumPy vs PyTorch
Most people stall on the same page of the Transformer paper: Attention(Q, K, V) = softmax(QK^T / sqrt(d_k))V. You can memorize the formula and still not be able to write it. Where do Q, K, V come from? Which step divides by d_k? Does the mask go before or after the softmax? Why is the reshape in multi-head attention ordered the way it is?
I tested this the slow way: I wrote the whole thing from scratch in NumPy, then ran it against PyTorch 2.14.1 scaled_dot_product_attention and nn.MultiheadAttention with identical weights and identical inputs. The largest absolute difference between my version and PyTorch's was 2.220446049250313e-16, which is floating-point noise. Every number below was produced on my own Windows CPU on October 12, 2026, and the code is copy-paste runnable.
Disclosure: there are no affiliate links in this article and nothing here is sponsored.
TL;DR
- One step in the formula is genuinely tricky: the division by sqrt(d_k) before the softmax. Without it, at d_k=64 I measured a max attention weight of 0.921946, which means the head looks at essentially one position.
- Multi-head is not "run attention a few more times": it splits d_model into h slices, runs each one independently, then concatenates and projects again.
- The mask goes into the softmax input, and its value is -inf, not 0. With 0, masked positions still receive weight.
- Matching the official implementation takes about 40 lines: single-head difference 2.22e-16, multi-head difference 4.44e-16.
The division is the part worth reading twice
Section 3.2.1 of "Attention Is All You Need" (arXiv:1706.03762, v1 submitted June 12, 2017; version v7 cited here revised August 2, 2023) states the formula as:
Attention(Q, K, V) = softmax(QK^T / sqrt(d_k)) V
The same section gives the reason, in the authors' words: "We suspect that for large values of d_k, the dot products grow large in magnitude, pushing the softmax function into regions where it has extremely small gradients."
The footnote does the math: if the components of q and k are independent with mean 0 and variance 1, then q·k = Σ q_i·k_i has mean 0 and variance d_k. Higher dimension, larger dot products, and the softmax gets pushed into a region with almost no gradient — the output is nearly one-hot and backpropagation has almost nothing to work with.
Section 3.2.2 gives the base configuration: d_model = 512, h = 8 heads, d_k = d_v = d_model / h = 64, N = 6 identical layers in each stack, dropout rate 0.1.
What happens without it: watching softmax collapse
Rather than restate the claim, I ran it. This generates 2,000 keys for d_k = 16 / 64 / 512 and compares the softmax shape with and without scaling:
import numpy as np
def softmax(x):
x = x - x.max()
e = np.exp(x)
return e / e.sum()
for dk in (16, 64, 512):
r = np.random.default_rng(7)
q = r.normal(size=(dk,))
k = r.normal(size=(2000, dk))
dots = k @ q
w_raw = softmax(dots)
w_scl = softmax(dots / np.sqrt(dk))
print(f"dk={dk:4d} dot std={dots.std():6.2f} "
f"max weight unscaled={w_raw.max():.6f} max weight scaled={w_scl.max():.6f} "
f"effective terms unscaled={1/np.sum(w_raw**2):7.2f} scaled={1/np.sum(w_scl**2):7.2f}")Output on my machine:
dk= 16 dot std= 2.47 max weight unscaled=0.227685 max weight scaled=0.004241 effective terms unscaled= 14.85 scaled=1341.59
dk= 64 dot std= 7.32 max weight unscaled=0.921946 max weight scaled=0.009399 effective terms unscaled= 1.17 scaled= 829.74
dk= 512 dot std= 20.99 max weight unscaled=0.933175 max weight scaled=0.005796 effective terms unscaled= 1.15 scaled= 917.82"Effective terms" is my own metric, defined as 1 / Σ w_i², read as "how many of those 2,000 keys are actually being attended to."
The d_k=64 row makes the point: unscaled, the max weight is 92.19% and effective terms is 1.17 — across a 2,000-token sequence, attention lands almost entirely on one position and the gradient for the rest is close to zero. Scaled, the max weight drops to 0.94% and effective terms rises to 829.74, so the head is actually aggregating. At d_k=16 the gap is smaller (14.85 vs 1341.59), which lines up with the paper's own "for small values of d_k the two mechanisms perform similarly."
The minimal runnable implementation: under 40 lines of NumPy
import numpy as np
def softmax(x, axis=-1):
x = x - np.max(x, axis=axis, keepdims=True) # keeps exp from overflowing
e = np.exp(x)
return e / np.sum(e, axis=axis, keepdims=True)
def scaled_dot_product_attention(Q, K, V, mask=None):
"""Q:(..., Tq, dk) K:(..., Tk, dk) V:(..., Tk, dv)
mask: bool, True means this position gets masked out (replaced by -inf)
"""
dk = Q.shape[-1]
scores = Q @ np.swapaxes(K, -1, -2) / np.sqrt(dk)
if mask is not None:
scores = np.where(mask, -np.inf, scores)
w = softmax(scores, axis=-1)
return w @ V, wThat is a one-to-one translation of the paper's equation, plus np.where(mask, -np.inf, scores) and the max subtraction inside softmax (numerical stability, not decoration).
With 5 tokens and d_k=4 (seed fixed at 20261012), the first row of attention weights comes out as:
[0.027843 0.167136 0.246698 0.211041 0.347282]Five positions, five weights, summing to 1. I asserted the row sums: the largest deviation from 1.0 was 1.1102230246251565e-16. Measure this yourself — if you ever see that deviation at 1e-3, suspect the softmax axis before you suspect NumPy.
Multi-head: split, attend separately, concatenate
Section 3.2.2 describes it as projecting Q, K, V h times with different learned projections down to d_k, d_k, d_v dimensions, attending in parallel, concatenating, and projecting once more. In code that is a reshape plus a transpose:
def multi_head_attention(X, WQ, WK, WV, WO, h, causal=False):
T, d_model = X.shape
dh = d_model // h
Q, K, V = X @ WQ, X @ WK, X @ WV
# (T, d_model) -> (h, T, dh)
Qh = Q.reshape(T, h, dh).transpose(1, 0, 2)
Kh = K.reshape(T, h, dh).transpose(1, 0, 2)
Vh = V.reshape(T, h, dh).transpose(1, 0, 2)
mask = np.triu(np.ones((T, T), dtype=bool), k=1) if causal else None
heads, w = scaled_dot_product_attention(Qh, Kh, Vh, mask=mask)
# (h, T, dh) -> (T, d_model)
return heads.transpose(1, 0, 2).reshape(T, d_model) @ WO, wWith 5 tokens, d_model=16, h=4: output shape (5, 16), weights shape (4, 5, 5) — one 5×5 attention table per head.
One thing worth knowing about reshape(T, h, dh): it slices along the last axis in order, so head 0 gets columns 0..dh-1. That is why WQ, WK, WV should stay full d_model × d_model matrices and let each head learn its own block, rather than pre-slicing the weights into h pieces and initializing them separately. The latter runs, but it is not what the paper describes.
Section 3.2.3 also lists three uses of attention: encoder self-attention (Q, K, V from the same place), decoder self-attention with the causal mask, and encoder-decoder attention (Q from the previous decoder layer, K and V from the encoder output stack). The implementation above passes X as all three, which covers the first two; for the third, hand it the encoder output as K and V and the function needs no changes.
Causal masking: the paper says -inf, not 0
Section 3.2.3 is blunt about it: "We need to prevent leftward information flow in the decoder to preserve the auto-regressive property. We implement this inside of scaled dot-product attention by masking out (setting to −∞) all values in the input of the softmax which correspond to illegal connections."
Two facts in that sentence: the mask applies to the input of the softmax, and the value is −∞.
I verified it on the same 5-token example with a causal mask: the largest weight at a masked position was 0.0, and row 2 had only 3 positions with weight above zero, meaning position 2 can see 0, 1, 2 — matching "up to and including that position."
If your masked positions come out as 1e-9 rather than exactly 0, you used a float mask with -1e9 instead of -inf. It works numerically, but they are different things: -inf makes exp produce exactly 0, -1e9 only makes it produce a very small number.
Checking my code against PyTorch
My NumPy version and F.scaled_dot_product_attention, same weights and inputs in float64: max absolute difference 2.220446049250313e-16. Loading WQ/WK/WV into nn.MultiheadAttention's in_proj_weight and WO into out_proj.weight and running the full multi-head path: max absolute difference 4.440892098500626e-16.
One honest performance number from the same CPU: at T=128, d=64, the per-row for-loop version took 5.669 ms and the matrix version 1.050 ms, 5.4x faster, with outputs differing by 4.00e-15. The loop version is not wrong, it just splits one matrix multiply into 128 pieces — which is what the paper means by dot-product attention being "implemented using highly optimized matrix multiplication code."
Here is the trap. The PyTorch 2.14 docs flag it in a Note on the scaled_dot_product_attention page: boolean mask semantics are the inverse of MultiheadAttention's key_padding_mask. In scaled_dot_product_attention, True means the element takes part in attention; in MultiheadAttention, True means masked out.
I ran both ways with the same causal mask:
correct convention (inverted before passing): 2.220446049250313e-16
reversed convention (passed as-is): 3.668259537831219A difference of 3.67, and no error is raised — the model quietly attends to future tokens. If you take one thing from this article, take that.
Troubleshooting: five real errors
All five were captured while I was writing this code, copied verbatim. They are not illustrative examples.
Error 1: forgot to transpose K (NumPy)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0,
with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 5 is different from 4)Q is (5, 4) and so is K, so Q @ K cannot work. Write Q @ K.T, or np.swapaxes(K, -1, -2) for batched input. In self-attention Q and K have identical shapes, so the error never says "transpose" — if you see "size 5 is different from 4", check the transpose first.
Error 2: forgot to transpose K (PyTorch)
RuntimeError: mat1 and mat2 shapes cannot be multiplied (5x4 and 5x4)Same mistake, more direct message. Build the habit of writing K.transpose(-2, -1) in full.
Error 3: is_causal and attn_mask at the same time
RuntimeError: _scaled_dot_product_attention: Explicit attn_mask should not be set when is_causal=TrueThe docs state "An error is thrown if both attn_mask and is_causal are set." Pick one: is_causal=True for a standard square causal mask, or pass attn_mask for non-square or custom masks and leave is_causal at its default.
Error 4: exp underflows in float32 and softmax returns nan
dot = -1802.824
exp(dot) = 0.0
e/e = nan
RuntimeWarning: invalid value encountered in scalar divideMeasured with 1024-dimensional float32 vectors scaled by 10. The dot product is -1802.824, exp underflows to 0.0, numerator and denominator are both zero, and the division gives nan. With a large positive dot product you get inf instead, and inf/inf is nan too. The fix is not float64, it is subtracting the max before exp: after x - x.max() the largest argument to exp is 0, so exp(0)=1 and the denominator is at least 1. The implementation above already does this.
Error 5: wrong softmax axis — silent, and wrong
axis=0 row sums: [0.838357 1.937897 0.662364 0.992172 0.569209]
axis=-1 row sums: [1. 1. 1. 1. 1.]
output difference: 1.2147038513899666np.exp(s) / np.exp(s).sum(axis=0) runs fine, warns about nothing, and has the right shape — it just normalizes over columns instead of rows. The check is one line: assert that w.sum(axis=-1) is all ones. This is the case where an assertion is not optional, because nothing else will tell you.
Where to go next
- To actually run models instead of simulating them in NumPy: Docker Model Runner architecture explained
- Why retrieval misses relevant content even when it is in the index: Why RAG retrieval is wrong
- For a reading list that covers the fundamentals: Technical books for programmers, ranked
FAQ
Q1: Can I train with the from-scratch version?
No, and you shouldn't want to. Its value is confirming that you understood the algorithm: matching PyTorch to 1e-16 means the formula, axis order, and mask placement are all right. For training, use F.scaled_dot_product_attention, which selects among FlashAttention, Memory-Efficient attention, and the C++ implementation (the three listed in the docs), and supports GQA.
Q2: Why is my difference 1e-7 instead of 1e-16?
Check the dtype first. Everything above was measured in float64; float32 lands around 1e-7, which is expected. The docs also note that fused floating-point kernels can produce different results from each other, and that the C++ implementation supports float64.
Q3: Can I use -1e9 instead of -inf in the mask?
It runs, but it is not exact. -inf makes exp return exactly 0; -1e9 leaves a tiny non-zero value in float32. The 0.0 I measured at masked positions came from using -inf. Also note that PyTorch's float attn_mask is added to the scores, so adding -inf is what masks.
Q4: Must d_k equal d_model / h?
That is how the paper sets it (d_k = d_v = d_model/h = 64, h=8, d_model=512), with the stated reason that with the reduced per-head dimension, "the total computational cost is similar to that of single-head attention with full dimensionality." It is a convention, not a mathematical requirement, and modern models do deviate.
Sources and how the numbers were produced
The formula, the justification for dividing by sqrt(d_k), h=8 / d_k=d_v=64 / d_model=512 / N=6 / dropout 0.1, and the -inf masking all come from sections 3.2.1–3.2.3 of arXiv:1706.03762, "Attention Is All You Need" (v1 June 12, 2017; v7 cited here, August 2, 2023).
Every measured number (2.22e-16, 4.44e-16, 0.921946, effective terms 1.17, the 5.4x timing, 3.668259537831219, and the five error messages) was produced on my Windows CPU on October 12, 2026 with Python 3.13, NumPy 2.5.4, and PyTorch 2.14.1+cpu. Seeds are hardcoded in the snippets, so the runs are reproducible.
The boolean mask semantics, the is_causal/attn_mask exclusivity, and the backend selection statements come from the PyTorch 2.14 documentation page for scaled_dot_product_attention, cross-checked against what I measured locally.
No affiliate links and no product recommendations here. The timing numbers depend on CPU model and load and will differ on your machine; they are here to show the order of magnitude, not to serve as a benchmark.
👉 Join Xiaomi MiMo Platform: Leading AI model platform with cost-effective inference
👉 Join Aliyun AI: Top AI products with exclusive coupons for business innovation
📌 This article was AI-assisted generated and human-reviewed | TechPassive — An AI-driven content testing site focused on real tool reviews
🔗 Recommended Tools
These are carefully selected tools. Using our affiliate links supports us to keep producing quality content: