Foundations all 4 written
The Python, the vectors and the calculus you actually use. Nothing decorative.
4 The Python an AI engineer writes
You do not need all of Python. You need the shapes data arrives in, the two or three constructs that keep the code short, and one honest understanding of how floats behave. That last one is not decoration: every similarity score, probability and gradient in the rest of this course is a float, and comparing floats carelessly is the most common way a correct implementation looks broken.
The shapes that actually turn up
Almost every AI pipeline moves the same object around: a list of dicts. Documents, chunks, messages, tool calls, evaluation cases. Learn to map, filter and index that shape and you can read most Python in this field.
# A list of dicts, and the four things you do to it constantly.
DOCS = [
{"id": "d1", "title": "Refund policy", "words": 210, "tags": ["policy", "support"]},
{"id": "d2", "title": "Shipping zones", "words": 940, "tags": ["logistics"]},
{"id": "d3", "title": "Returns FAQ", "words": 130, "tags": ["policy", "faq"]},
]
titles = [d["title"] for d in DOCS] # map
short = [d["id"] for d in DOCS if d["words"] < 300] # filter
total = sum(d["words"] for d in DOCS) # reduce, and lazily: no list is built
by_id = {d["id"]: d for d in DOCS} # index once, look up many times
# Counting without a library. dict.get with a default is the whole technique.
tag_counts = {}
for d in DOCS:
for tag in d["tags"]:
tag_counts[tag] = tag_counts.get(tag, 0) + 1
# Sorting: key= takes a function. Negating the number is how you sort descending.
longest = sorted(DOCS, key=lambda d: -d["words"])[0]
print("titles :", titles)
print("short :", short)
print("words :", total)
print("tags :", tag_counts)
print("longest:", longest["title"])
test("three documents go in and three titles come out", len(titles), 3)
test("two of the three documents are under three hundred words", short, ["d1", "d3"])
test("the word counts add up to 1280", total, 1280)
test("policy is the only tag that appears twice", [t for t, n in tag_counts.items() if n == 2], ["policy"])
test("shipping zones is the longest document", longest["title"], "Shipping zones")
test("the index finds a document by its id", by_id["d3"]["title"], "Returns FAQ")What it printed, when the build ran it
titles : ['Refund policy', 'Shipping zones', 'Returns FAQ']
short : ['d1', 'd3']
words : 1280
tags : {'policy': 2, 'support': 1, 'logistics': 1, 'faq': 1}
longest: Shipping zones
ok three documents go in and three titles come out
ok two of the three documents are under three hundred words
ok the word counts add up to 1280
ok policy is the only tag that appears twice
ok shipping zones is the longest document
ok the index finds a document by its id
If you write JavaScript, you already know all four operations under different names. The panel below is the same work in the language you came from, so you can read the Python by translation rather than by memorisation. It is the only JavaScript panel in this lesson, and unlike the Python it runs in your browser.
const same = (a, b) => JSON.stringify(a) === JSON.stringify(b);
const test = (name, got, want) => console.log(
(same(got, want) ? 'ok ' : 'FAIL ') + name +
(same(got, want) ? '' : ` got ${JSON.stringify(got)}, want ${JSON.stringify(want)}`)
);
// Never use === on a computed float. See the Foundations lesson on floating point.
const close = (a, b, tol = 1e-9) => Math.abs(a - b) < tol;
// The same four operations, in the language you already know.
const DOCS = [
{ id: 'd1', title: 'Refund policy', words: 210 },
{ id: 'd2', title: 'Shipping zones', words: 940 },
{ id: 'd3', title: 'Returns FAQ', words: 130 },
];
const titles = DOCS.map((d) => d.title); // [d["title"] for d in DOCS]
const short = DOCS.filter((d) => d.words < 300).map((d) => d.id);
const total = DOCS.reduce((n, d) => n + d.words, 0); // sum(d["words"] for d in DOCS)
const byId = Object.fromEntries(DOCS.map((d) => [d.id, d])); // {d["id"]: d for d in DOCS}
console.log('titles :', titles);
console.log('short :', short);
console.log('words :', total);
test('map is a comprehension with no condition', titles.length, 3);
test('filter then map is a comprehension with an if', short, ['d1', 'd3']);
test('reduce is sum(), spelled out at length', total, 1280);
test('an object built from entries is the dict comprehension', byId.d3.title, 'Returns FAQ');Output
One real difference is worth holding on to. A Python comprehension inside square
brackets builds the whole list in memory; the same expression inside round brackets is a
generator and builds nothing until something iterates it. sum(d["words"]
for d in DOCS) never allocates a list. On four documents that is trivia. On four
million rows it is the difference between working and not.
Floating point: the reason every float check here uses close()
A float is a binary approximation of a decimal number. Most decimals have no exact binary form, so arithmetic accumulates a tiny error. This is not a Python quirk, it is how the hardware works, and it is the same in JavaScript.
# The one panel in this course that compares floats with ==, so you can watch it fail.
a = 0.1 + 0.2
print("0.1 + 0.2 =", a)
print("is it 0.3? =", a == 0.3)
print("close to 0.3? =", close(a, 0.3))
# The error compounds. Add a tenth ten times and you do not arrive at one.
running = 0.0
for _ in range(10):
running += 0.1
print()
print("0.1 added ten times =", repr(running))
print("short of 1.0 by =", abs(running - 1.0))
# Worse: the order of a sum changes the answer, so addition is not associative.
up = sum([1e16, 1.0, 1.0, -1e16])
down = sum([1e16, -1e16, 1.0, 1.0])
print()
print("1e16 + 1 + 1 - 1e16 =", up, " <- both ones were swallowed")
print("1e16 - 1e16 + 1 + 1 =", down)
test("a tenth plus two tenths is not three tenths", a == 0.3, False)
test("but it is close enough for any real purpose", close(a, 0.3), True)
test("adding a tenth ten times does not reach one either", running == 1.0, False)
test("close() forgives that too", close(running, 1.0), True)
test("the error is vanishingly small, which is exactly why it hides", abs(a - 0.3) < 1e-15, True)
test("changing the order of a sum changes its result", up == down, False)What it printed, when the build ran it
0.1 + 0.2 = 0.30000000000000004 is it 0.3? = False close to 0.3? = True 0.1 added ten times = 0.9999999999999999 short of 1.0 by = 1.1102230246251565e-16 1e16 + 1 + 1 - 1e16 = 0.0 <- both ones were swallowed 1e16 - 1e16 + 1 + 1 = 2.0 ok a tenth plus two tenths is not three tenths ok but it is close enough for any real purpose ok adding a tenth ten times does not reach one either ok close() forgives that too ok the error is vanishingly small, which is exactly why it hides ok changing the order of a sum changes its result
The rule for the rest of this course. Never assert equality on a
computed float. Use close(a, b), which is the math.isclose idea
with an absolute tolerance. Two follow-on habits matter as much: do not accumulate a long
sum in a loop when the numbers vary wildly in size, and never test a computed value for
being exactly zero — test that it is smaller than a threshold you chose on purpose. The
vectors lesson has a division that blows up precisely there.
Generators, because model output arrives one token at a time
Streaming is the normal case when you call a model, so the generator is not an advanced
topic here, it is the default shape. A function with yield in it returns
something that produces values on demand and remembers where it stopped.
# Generators: the shape you meet the moment you stream tokens out of a model.
import itertools
def token_stream(text):
"""Yields one token at a time. Nothing happens until somebody asks for one."""
for word in text.split():
yield word
def take(iterable, n):
"""The first n items, and no more work than that."""
out = []
for item in iterable:
out.append(item)
if len(out) == n:
break
return out
stream = token_stream("the quick brown fox jumps over the lazy dog")
first_three = take(stream, 3)
rest = list(stream) # the SAME generator, resumed where it left off
print("first three :", first_three)
print("what remained:", rest)
print("second pass :", list(stream), "<- a generator is consumed, not rewound")
# A generator expression is a lazy comprehension. This allocates nothing.
squares = (n * n for n in range(1000000))
print("first five squares:", take(squares, 5))
# itertools writes the awkward plumbing for you, correctly.
counted = list(itertools.islice(itertools.count(10), 4))
print("four items from an endless counter:", counted)
test("asking for three tokens gives three tokens", first_three, ["the", "quick", "brown"])
test("the generator resumes at the fourth token, not the first", rest[0], "fox")
test("nothing was lost between the two reads", len(first_three) + len(rest), 9)
test("reading a used generator again yields nothing at all", list(stream), [])
test("a million squares cost nothing until you take five", take((n * n for n in range(1000000)), 5), [0, 1, 4, 9, 16])
test("an endless counter is safe as long as you slice it", counted, [10, 11, 12, 13])What it printed, when the build ran it
first three : ['the', 'quick', 'brown'] what remained: ['fox', 'jumps', 'over', 'the', 'lazy', 'dog'] second pass : [] <- a generator is consumed, not rewound first five squares: [0, 1, 4, 9, 16] four items from an endless counter: [10, 11, 12, 13] ok asking for three tokens gives three tokens ok the generator resumes at the fourth token, not the first ok nothing was lost between the two reads ok reading a used generator again yields nothing at all ok a million squares cost nothing until you take five ok an endless counter is safe as long as you slice it
The trap is in the third print. A generator can be read once. Code that iterates a
stream to count it and then iterates again to use it gets an empty second pass, and in a
pipeline that failure looks like "the model returned nothing" rather than "I consumed it
twice". If you need it twice, materialise it once with list(...) and pay the
memory on purpose.
Shapes with names, and two traps in the way you declare them
Dicts are fine inside a function. At a boundary — something you index, cache, log or hand to another module — give the shape a name. A dataclass is the cheapest way, and freezing it stops anything downstream editing a record that has already been indexed.
# Dataclasses for the shapes that cross a boundary, and the two traps nearby.
from dataclasses import dataclass, field, replace
from typing import Dict, List, Optional, Union
@dataclass(frozen=True)
class Chunk:
"""A slice of a document on its way to being embedded.
Python here is 3.9, which has no int | str unions and no match statement.
A maybe-vector is spelled Optional[List[float]], and a value that is either a
string or a number is Union[str, int].
"""
doc_id: str
text: str
vector: Optional[List[float]] = None
meta: Dict[str, Union[str, int]] = field(default_factory=dict)
def is_embedded(self) -> bool:
return self.vector is not None
raw = Chunk(doc_id="d1", text="refunds are accepted within fourteen days")
embedded = replace(raw, vector=[0.2, -0.1, 0.4], meta={"source": "policy.md", "page": 3})
print("raw :", raw.doc_id, "embedded?", raw.is_embedded())
print("embedded :", embedded.doc_id, "embedded?", embedded.is_embedded(),
"dims", len(embedded.vector))
try:
raw.text = "something else"
frozen = False
except Exception as error:
frozen = True
print("editing a frozen chunk raises", type(error).__name__)
# Trap one. A mutable default is built once, when the function is defined, and then
# shared by every call for the lifetime of the process.
def collect(item, into=[]):
into.append(item)
return into
first = collect("a")
second = collect("b")
print()
print("collect('a') ->", first)
print("collect('b') ->", second, "<- not ['b']")
def collect_safely(item, into: Optional[List[str]] = None) -> List[str]:
"""The fix, and the only correct spelling of it."""
if into is None:
into = []
into.append(item)
return into
# Trap two. bool is a subclass of int, so a boolean passes a numeric check.
def looks_numeric(value: object) -> bool:
return isinstance(value, (int, float))
def is_numeric(value: object) -> bool:
return isinstance(value, (int, float)) and not isinstance(value, bool)
test("a fresh chunk has no vector yet", raw.is_embedded(), False)
test("replace() returns a new chunk and leaves the original alone",
(embedded.is_embedded(), raw.vector), (True, None))
test("a frozen dataclass refuses to be edited in place", frozen, True)
test("a mutable default argument is shared between calls", second, ["a", "b"])
test("both calls handed back the very same list object", first is second, True)
test("defaulting to None instead keeps every call independent",
[collect_safely("a"), collect_safely("b")], [["a"], ["b"]])
test("True passes a naive number check and fails a careful one",
[looks_numeric(True), is_numeric(True)], [True, False])What it printed, when the build ran it
raw : d1 embedded? False
embedded : d1 embedded? True dims 3
editing a frozen chunk raises FrozenInstanceError
collect('a') -> ['a', 'b']
collect('b') -> ['a', 'b'] <- not ['b']
ok a fresh chunk has no vector yet
ok replace() returns a new chunk and leaves the original alone
ok a frozen dataclass refuses to be edited in place
ok a mutable default argument is shared between calls
ok both calls handed back the very same list object
ok defaulting to None instead keeps every call independent
ok True passes a naive number check and fails a careful one
What a real project would use instead
Nothing above needed a library, and that is deliberate: this course has to verify itself with the standard library alone. In a real pipeline the numeric work goes to numpy, and the tabular work to pandas, for one reason worth stating plainly. A Python-level loop over a million floats runs in the interpreter; numpy pushes the same loop into compiled code over a contiguous block of memory, which is typically tens to hundreds of times faster and, just as importantly, shorter to read.
# The same operations with the libraries a real project would use.
# NOT RUN HERE - numpy and pandas are not installed.
import numpy as np
import pandas as pd
# One embedding: a fixed-width array of float32, not a Python list of float64.
vector = np.array([0.2, -0.1, 0.4], dtype=np.float32)
# A whole corpus: 10,000 embeddings of 768 dimensions, as one 2-D block.
matrix = np.zeros((10000, 768), dtype=np.float32)
# Every similarity in the next lesson, for all 10,000 documents, in one call.
scores = matrix @ vector[:768] # the loop happens in compiled code
# The list-of-dicts from the first panel, as a dataframe.
docs = pd.DataFrame([
{"id": "d1", "title": "Refund policy", "words": 210},
{"id": "d2", "title": "Shipping zones", "words": 940},
])
short = docs[docs["words"] < 300] # the filter comprehension
total = docs["words"].sum() # the sum
by_id = docs.set_index("id") # the dict comprehensionnumpy and pandas are not installed in this environment, so this panel was never executed. Every operation in it has a plain-Python equivalent above that did run. Read it as the shape to reach for in a real project, not as verified output.
The useful reading order, if you want to be efficient about this: comprehensions,
dicts, generators, dataclasses, typing. That is roughly ninety per cent of
the Python you will write as an AI engineer. Decorators, metaclasses, async and the
descriptor protocol can wait until something you are reading actually uses them.
The two traps in this lesson are not equally visible. A frozen
dataclass or a shared mutable default fails loudly and once. Floating point fails
quietly and forever: a cosine similarity of 0.9999999999999998 compared with
1.0 is a test that fails for no reason a reader can see, and it is the
single most common bug in hand-written vector code. Every float assertion from here on
uses close, and now you know why.
5 Vectors, distance and similarity
A vector is a list of numbers. An embedding is a list of numbers produced by a model, where nearby lists mean similar things. That is the whole idea, and the arithmetic that exploits it fits in three short functions: a dot product, a length, and a division. This lesson writes all three from scratch and then shows the two places they mislead you.
The dot product, which everything else is made of
Multiply the two vectors position by position and add up the results. One number comes out. It is large and positive when the vectors point the same way, near zero when they are unrelated, and negative when they oppose.
# The dot product. Every similarity score in this course is built on this.
def dot(a, b):
if len(a) != len(b):
raise ValueError("vectors must be the same length, got {} and {}".format(len(a), len(b)))
return sum(x * y for x, y in zip(a, b))
a = [1.0, 2.0, 3.0]
b = [4.0, -5.0, 6.0]
print("a . b =", dot(a, b), " (4 - 10 + 18)")
print("a . a =", dot(a, a), " <- a vector dotted with itself is its squared length")
print("[1,0] . [0,1] =", dot([1.0, 0.0], [0.0, 1.0]), " <- perpendicular")
print("[1,1] . [-1,-1] =", dot([1.0, 1.0], [-1.0, -1.0]), " <- opposed")
# The length check is not defensive clutter. zip() stops at the shorter input and
# says nothing, so without it a 768-dimension query silently scores against a
# 1536-dimension document and returns a plausible, wrong number.
print()
print("zip would have quietly truncated:", list(zip([1, 2, 3], [1, 2])))
try:
dot([1.0, 2.0, 3.0], [1.0, 2.0])
refused = False
except ValueError as error:
refused = True
print("so we refuse it instead:", error)
test("the dot product of a and b is twelve", close(dot(a, b), 12.0), True)
test("a vector dotted with itself gives its squared length", close(dot(a, a), 14.0), True)
test("perpendicular vectors dot to zero", close(dot([1.0, 0.0], [0.0, 1.0]), 0.0), True)
test("opposed vectors give a negative number", dot([1.0, 1.0], [-1.0, -1.0]) < 0.0, True)
test("the order of the two vectors makes no difference", close(dot(a, b), dot(b, a)), True)
test("mismatched lengths are refused rather than truncated", refused, True)What it printed, when the build ran it
a . b = 12.0 (4 - 10 + 18) a . a = 14.0 <- a vector dotted with itself is its squared length [1,0] . [0,1] = 0.0 <- perpendicular [1,1] . [-1,-1] = -2.0 <- opposed zip would have quietly truncated: [(1, 1), (2, 2)] so we refuse it instead: vectors must be the same length, got 3 and 2 ok the dot product of a and b is twelve ok a vector dotted with itself gives its squared length ok perpendicular vectors dot to zero ok opposed vectors give a negative number ok the order of the two vectors makes no difference ok mismatched lengths are refused rather than truncated
Length, and what normalising actually does
The length of a vector — its norm — is Pythagoras in as many dimensions as you have. Divide a vector by its own length and you get a unit vector: same direction, length one. That division is where hand-written vector code breaks.
# Length, unit vectors, and the division that has to be guarded.
import math
def dot(a, b):
return sum(x * y for x, y in zip(a, b))
def norm(v):
"""The length of v. Pythagoras, in as many dimensions as you like."""
return math.sqrt(sum(x * x for x in v))
def unit(v):
"""Same direction, length one. Refuses anything too short to have a direction."""
length = norm(v)
if length < 1e-12:
raise ValueError("a vector this short has no reliable direction")
return [x / length for x in v]
print("norm([3, 4]) =", norm([3.0, 4.0]), " <- the 3-4-5 triangle")
print("norm([1, 2, 3]) =", norm([1.0, 2.0, 3.0]))
print("unit([3, 4]) =", unit([3.0, 4.0]))
print("norm of that unit =", norm(unit([3.0, 4.0])))
# The guard matters because the failure is not a clean error. A vector of tiny
# numbers - an empty chunk, a padded row, a de-meaned constant feature - divides
# into nonsense long before it divides by zero.
tiny = [1e-200, 1e-200]
print()
print("norm of a tiny vector:", norm(tiny), "<- underflowed to exactly zero")
try:
unit(tiny)
refused = False
except ValueError as error:
refused = True
print("so we refuse it:", error)
test("the three four five triangle has length five", close(norm([3.0, 4.0]), 5.0), True)
test("a unit vector has length one", close(norm(unit([3.0, 4.0])), 1.0), True)
test("the length is the square root of the vector dotted with itself",
close(norm([1.0, 2.0, 3.0]), math.sqrt(dot([1.0, 2.0, 3.0], [1.0, 2.0, 3.0]))), True)
test("doubling a vector doubles its length", close(norm([6.0, 8.0]), 2.0 * norm([3.0, 4.0])), True)
test("normalising keeps the direction, so the unit vectors agree",
close(dot(unit([3.0, 4.0]), unit([30.0, 40.0])), 1.0), True)
test("a vector too short to have a direction is refused", refused, True)What it printed, when the build ran it
norm([3, 4]) = 5.0 <- the 3-4-5 triangle norm([1, 2, 3]) = 3.7416573867739413 unit([3, 4]) = [0.6, 0.8] norm of that unit = 1.0 norm of a tiny vector: 0.0 <- underflowed to exactly zero so we refuse it: a vector this short has no reliable direction ok the three four five triangle has length five ok a unit vector has length one ok the length is the square root of the vector dotted with itself ok doubling a vector doubles its length ok normalising keeps the direction, so the unit vectors agree ok a vector too short to have a direction is refused
Notice the guard is length < 1e-12 and not length == 0.0.
An exact zero is the easy case. The one that ruins an afternoon is a vector of very small
numbers whose norm underflows or nearly does: the division produces enormous values or
inf, the similarity score comes out as nan, and
nan compares false against everything, so the document quietly ranks last
instead of raising.
Cosine similarity: the dot product with the length divided out
Cosine similarity is the dot product of the two unit vectors. It answers one question — do these point the same way — and deliberately ignores how long either one is. The answer runs from 1 for identical directions, through 0 for unrelated, to -1 for opposite.
# Cosine similarity, built from the two previous functions.
import math
def dot(a, b):
return sum(x * y for x, y in zip(a, b))
def norm(v):
return math.sqrt(sum(x * x for x in v))
def cosine(a, b):
denominator = norm(a) * norm(b)
if denominator < 1e-12:
raise ValueError("cannot compare direction when one vector has no length")
return dot(a, b) / denominator
a = [1.0, 2.0, 3.0]
b = [4.0, -5.0, 6.0]
a_scaled = [10.0 * x for x in a]
print("cosine(a, a) =", cosine(a, a))
print("cosine(a, 10a) =", cosine(a, a_scaled), " <- magnitude thrown away")
print("cosine(a, b) = {:.6f}".format(cosine(a, b)))
print("cosine([1,0], [0,1]) =", cosine([1.0, 0.0], [0.0, 1.0]))
print("cosine([1,2], [-1,-2])=", cosine([1.0, 2.0], [-1.0, -2.0]))
# Now the sharp edge. Some vectors come out more similar to themselves than is legal.
awkward = [0.3, 0.5, 0.7, 0.11]
self_score = cosine(awkward, awkward)
print()
print("cosine(awkward, awkward) =", repr(self_score), "<- larger than 1.0")
try:
math.acos(self_score)
blew_up = False
except ValueError as error:
blew_up = True
print("math.acos on it raises : ValueError:", error)
print("clamping to [-1, 1] first:", math.acos(max(-1.0, min(1.0, self_score))), "radians")
try:
cosine([0.0, 0.0], [1.0, 2.0])
refused = False
except ValueError as error:
refused = True
print("a zero vector is refused :", error)
test("a vector is perfectly similar to itself", close(cosine(a, a), 1.0), True)
test("making a vector ten times longer does not change its cosine", close(cosine(a, a_scaled), 1.0), True)
test("the opposite direction scores minus one", close(cosine([1.0, 2.0], [-1.0, -2.0]), -1.0), True)
test("perpendicular vectors score zero", close(cosine([1.0, 0.0], [0.0, 1.0]), 0.0), True)
test("rounding can push a self-similarity just above one", self_score > 1.0, True)
test("so turning that score into an angle raises, unless you clamp it first", blew_up, True)
test("a zero-length vector is refused instead of returning nan", refused, True)What it printed, when the build ran it
cosine(a, a) = 1.0 cosine(a, 10a) = 1.0 <- magnitude thrown away cosine(a, b) = 0.365487 cosine([1,0], [0,1]) = 0.0 cosine([1,2], [-1,-2])= -0.9999999999999998 cosine(awkward, awkward) = 1.0000000000000002 <- larger than 1.0 math.acos on it raises : ValueError: math domain error clamping to [-1, 1] first: 0.0 radians a zero vector is refused : cannot compare direction when one vector has no length ok a vector is perfectly similar to itself ok making a vector ten times longer does not change its cosine ok the opposite direction scores minus one ok perpendicular vectors score zero ok rounding can push a self-similarity just above one ok so turning that score into an angle raises, unless you clamp it first ok a zero-length vector is refused instead of returning nan
The awkward vector is the practical warning. Cosine similarity is
mathematically bounded by 1, but the computed value is not: rounding in the norms can leave
it at 1.0000000000000002. Anything that treats the score as a genuine cosine —
math.acos to get an angle, an angular distance metric, a
sqrt(1 - s*s) anywhere — raises or returns nan. Clamp to
[-1, 1] before you do arithmetic on a similarity score.
Euclidean distance, and how it relates to cosine
Distance is the other obvious question: how far apart are these two points. It is the length of the difference. Distance and cosine are not rivals — for unit vectors they are the same information, related by one identity that is worth verifying rather than believing.
# Euclidean distance, and the identity that ties it to cosine.
import math
def dot(a, b):
return sum(x * y for x, y in zip(a, b))
def norm(v):
return math.sqrt(sum(x * x for x in v))
def unit(v):
length = norm(v)
return [x / length for x in v]
def cosine(a, b):
return dot(a, b) / (norm(a) * norm(b))
def euclidean(a, b):
"""The length of the difference between the two vectors."""
return math.sqrt(sum((x - y) ** 2 for x, y in zip(a, b)))
a = [1.0, 2.0, 3.0]
b = [4.0, -5.0, 6.0]
a_scaled = [10.0 * x for x in a]
print("euclidean(a, b) = {:.6f}".format(euclidean(a, b)))
print("euclidean(a, a) =", euclidean(a, a))
print("euclidean([0,0],[3,4])=", euclidean([0.0, 0.0], [3.0, 4.0]))
# On unit vectors the two measures carry the same information:
# squared distance = 2 - 2 * cosine
d = euclidean(unit(a), unit(b))
c = cosine(a, b)
print()
print("unit distance squared = {:.15f}".format(d * d))
print("2 - 2 * cosine = {:.15f}".format(2.0 - 2.0 * c))
print("equal? == says", d * d == 2.0 - 2.0 * c, " close() says", close(d * d, 2.0 - 2.0 * c))
# And here is where they part company. Ten times longer, same direction.
print()
print("cosine(a, 10a) =", cosine(a, a_scaled), " <- identical")
print("euclidean(a, 10a) = {:.4f}".format(euclidean(a, a_scaled)), " <- very far apart")
test("the distance from a point to itself is zero", close(euclidean(a, a), 0.0), True)
test("the three four five triangle again, as a distance",
close(euclidean([0.0, 0.0], [3.0, 4.0]), 5.0), True)
test("distance does not care which way round you ask", close(euclidean(a, b), euclidean(b, a)), True)
test("for unit vectors, squared distance is two minus twice the cosine",
close(d * d, 2.0 - 2.0 * c), True)
test("the identity is true but not bit for bit, which is the whole point of close()",
d * d == 2.0 - 2.0 * c, False)
test("scaling a vector leaves its cosine alone but moves it a long way",
[close(cosine(a, a_scaled), 1.0), euclidean(a, a_scaled) > 1.0], [True, True])What it printed, when the build ran it
euclidean(a, b) = 8.185353 euclidean(a, a) = 0.0 euclidean([0,0],[3,4])= 5.0 unit distance squared = 1.269026115352193 2 - 2 * cosine = 1.269026115352193 equal? == says False close() says True cosine(a, 10a) = 1.0 <- identical euclidean(a, 10a) = 33.6749 <- very far apart ok the distance from a point to itself is zero ok the three four five triangle again, as a distance ok distance does not care which way round you ask ok for unit vectors, squared distance is two minus twice the cosine ok the identity is true but not bit for bit, which is the whole point of close() ok scaling a vector leaves its cosine alone but moves it a long way
That identity is why production vector databases normalise every embedding once, at write time, and then rank by plain dot product. Once every vector has length one, the dot product is the cosine, and nearest by distance is the same ordering as most similar by cosine. One multiplication per dimension, no square roots, same answer. If you ever see a store described as using "inner product" similarity, this is the assumption it is making, and it is your job to have normalised.
A search, and the bug you get for forgetting to normalise
Here is retrieval in miniature. Five documents, each hand-placed on four made-up axes so you can read the numbers, and one query. Ranking by cosine gives the sensible answer. Ranking by raw dot product on un-normalised vectors gives a different one, and the reason is not subtle.
# Retrieval in miniature: rank five documents against one query, two ways.
import math
def dot(a, b):
return sum(x * y for x, y in zip(a, b))
def norm(v):
return math.sqrt(sum(x * x for x in v))
def cosine(a, b):
return dot(a, b) / (norm(a) * norm(b))
AXES = ("animal", "tech", "food", "money")
LIBRARY = {
"cat photos": [0.9, 0.10, 0.0, 0.00],
"python tutorial": [0.0, 0.95, 0.0, 0.10],
"dog food prices": [0.7, 0.00, 0.6, 0.30],
"stock market news": [0.0, 0.20, 0.0, 0.95],
# Same direction as "cat photos", ten times longer. Think of a page that
# repeats itself, or a chunk that was never length-normalised.
"cat photos, long edition": [9.0, 1.00, 0.0, 0.00],
}
QUERY = [0.85, 0.05, 0.0, 0.0] # roughly "kittens"
def rank(query, library, score):
return sorted(((name, score(query, v)) for name, v in library.items()),
key=lambda pair: -pair[1])
by_cosine = rank(QUERY, LIBRARY, cosine)
by_dot = rank(QUERY, LIBRARY, dot)
print("axes:", AXES)
print()
print("ranked by cosine")
for name, score in by_cosine:
print(" {:26} {:.4f}".format(name, score))
print()
print("ranked by raw dot product")
for name, score in by_dot:
print(" {:26} {:.4f}".format(name, score))
print()
print("Length alone changed the winner. Nothing about the meaning changed.")
short_cos = cosine(QUERY, LIBRARY["cat photos"])
long_cos = cosine(QUERY, LIBRARY["cat photos, long edition"])
test("cosine puts a cat document at the top", by_cosine[0][0].startswith("cat photos"), True)
test("the two cat documents point the same way, so cosine cannot separate them",
close(short_cos, long_cos), True)
test("the raw dot product crowns the longer document instead",
by_dot[0][0], "cat photos, long edition")
test("so the two rankings disagree, purely because of length",
by_dot[0][0] == by_cosine[0][0], False)
test("the least relevant document is last on cosine", by_cosine[-1][0], "stock market news")
test("every cosine score sits inside minus one to one",
all(-1.0 <= s <= 1.0 for _, s in by_cosine), True)What it printed, when the build ran it
axes: ('animal', 'tech', 'food', 'money')
ranked by cosine
cat photos 0.9987
cat photos, long edition 0.9987
dog food prices 0.7207
python tutorial 0.0584
stock market news 0.0121
ranked by raw dot product
cat photos, long edition 7.7000
cat photos 0.7700
dog food prices 0.5950
python tutorial 0.0475
stock market news 0.0100
Length alone changed the winner. Nothing about the meaning changed.
ok cosine puts a cat document at the top
ok the two cat documents point the same way, so cosine cannot separate them
ok the raw dot product crowns the longer document instead
ok so the two rankings disagree, purely because of length
ok the least relevant document is last on cosine
ok every cosine score sits inside minus one to one
This is a real bug, not a toy one. Most embedding APIs return vectors that are already unit length, so the mistake stays hidden until you mix sources, cache a vector from an older model version, average several embeddings into one, or build your own from word counts. Then the longest document wins every search and the ranking looks vaguely plausible, which is worse than looking broken. Normalise at write time and assert the norm is 1 in a test.
What a real project would use instead
Everything above is one document at a time, and that is the wrong shape at scale. Real retrieval scores a whole corpus at once: put every embedding in one matrix and the query becomes a single matrix-vector product. numpy does that in compiled code over contiguous memory, and the line count goes down as well as the runtime.
# The same three measures, vectorised.
# NOT RUN HERE - numpy is not installed in this environment.
import numpy as np
a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, -5.0, 6.0])
a @ b # dot(a, b)
np.linalg.norm(a) # norm(a)
a / np.linalg.norm(a) # unit(a)
np.linalg.norm(a - b) # euclidean(a, b)
(a @ b) / (np.linalg.norm(a) * np.linalg.norm(b)) # cosine(a, b)
# The version that matters: 100,000 documents scored against one query at once.
corpus = np.random.randn(100000, 768).astype(np.float32)
corpus /= np.linalg.norm(corpus, axis=1, keepdims=True) # normalise ONCE, at write time
query = corpus[0]
scores = corpus @ query # every cosine similarity, one call
top10 = np.argpartition(-scores, 10)[:10] # the ten best, without sorting all 100,000
# Because the rows are unit length, the dot product IS the cosine, which is the
# identity the lesson verified by hand.numpy is not installed here, so this panel was never executed. Every line has a plain-Python counterpart above that did run and did assert. Treat it as the shape a real project uses, not as verified output.
Three functions, and you can now read any retrieval code you meet: dot for the raw agreement, norm for the length, cosine for direction alone. The choice between cosine and distance almost never matters once vectors are normalised, because the identity above makes them the same ordering. What does matter is remembering that a similarity score is a number about geometry, not about truth. Two texts can sit next to each other in embedding space and mean opposite things, which is the subject of the retrieval lessons later on.
6 Probability, softmax and sampling
A model's raw output is a list of unbounded real numbers called logits. They are not probabilities: they can be negative, they can be enormous, and they do not sum to anything in particular. Softmax converts them into a distribution, temperature decides how sharp that distribution is, and sampling picks one. Every generation setting you will ever tune is one of those three steps.
What makes a list of numbers a probability distribution
Two rules, and no others: nothing is negative, and the whole thing sums to one. The second rule is where floating point walks in, so the check has to be written with a tolerance from the very start.
# A distribution is two rules: nothing negative, and the whole lot sums to one.
def normalise(counts):
"""Turn counts into probabilities. The only honest way to make a distribution."""
total = sum(counts.values())
if total <= 0:
raise ValueError("cannot normalise counts that sum to zero")
return {name: n / total for name, n in counts.items()}
def is_distribution(p, tol=1e-9):
return all(x >= 0.0 for x in p) and close(sum(p), 1.0, tol)
observed = {"cat": 12, "dog": 7, "bird": 1}
probs = normalise(observed)
print("counts :", observed)
print("probabilities :", probs)
print("summing to :", repr(sum(probs.values())))
# Seven equally likely outcomes. Same two rules, and the sum now misses.
seven = normalise({name: 1 for name in "abcdefg"})
print()
print("one seventh :", repr(seven["a"]))
print("seven sevenths sum to :", repr(sum(seven.values())))
print("exactly 1.0? :", sum(seven.values()) == 1.0)
print()
print("[0.6, 0.3, 0.1] a distribution?", is_distribution([0.6, 0.3, 0.1]))
print("[0.6, 0.3, 0.2] a distribution?", is_distribution([0.6, 0.3, 0.2]), " <- sums to 1.1")
print("[0.9, -0.2, 0.3] a distribution?", is_distribution([0.9, -0.2, 0.3]), " <- negative")
test("twelve cats out of twenty animals is a probability of 0.6", close(probs["cat"], 0.6), True)
test("the three animal probabilities form a distribution", is_distribution(list(probs.values())), True)
test("seven sevenths do not add up to one in floating point", sum(seven.values()) == 1.0, False)
test("yet they are a distribution by any sane tolerance", is_distribution(list(seven.values())), True)
test("numbers that sum to more than one are not probabilities", is_distribution([0.6, 0.3, 0.2]), False)
test("and nothing containing a negative is either", is_distribution([0.9, -0.2, 0.3]), False)What it printed, when the build ran it
counts : {'cat': 12, 'dog': 7, 'bird': 1}
probabilities : {'cat': 0.6, 'dog': 0.35, 'bird': 0.05}
summing to : 1.0
one seventh : 0.14285714285714285
seven sevenths sum to : 0.9999999999999998
exactly 1.0? : False
[0.6, 0.3, 0.1] a distribution? True
[0.6, 0.3, 0.2] a distribution? False <- sums to 1.1
[0.9, -0.2, 0.3] a distribution? False <- negative
ok twelve cats out of twenty animals is a probability of 0.6
ok the three animal probabilities form a distribution
ok seven sevenths do not add up to one in floating point
ok yet they are a distribution by any sane tolerance
ok numbers that sum to more than one are not probabilities
ok and nothing containing a negative is either
If you ever write assert sum(probs) == 1.0 you have written a test that
fails on real data for no reason a reader can see. Seven equal probabilities are enough to
break it. Assert close(sum(probs), 1.0) instead, and when a library complains
that your probabilities do not sum to 1, this is usually why — not a bug in your
arithmetic.
Softmax, and the one line that stops it exploding
Softmax exponentiates every logit and divides by the total. Exponentiating makes
everything positive and exaggerates gaps; dividing makes it sum to one. Written literally
from the formula it is three lines, and it is unusable, because exp of a
realistic logit does not fit in a float.
# Softmax: the formula as written, then the version you can actually ship.
import math
def softmax_naive(logits):
"""The textbook formula, taken literally."""
exps = [math.exp(z) for z in logits]
total = sum(exps)
return [e / total for e in exps]
def softmax(logits):
"""The same function, shifted. Subtracting the max changes nothing and fixes everything."""
shift = max(logits)
exps = [math.exp(z - shift) for z in logits]
total = sum(exps)
return [e / total for e in exps]
small = [2.0, 1.0, 0.0]
big = [1000.0, 999.0, 998.0] # the same gaps, moved up by 998
p_small = softmax(small)
print("logits :", small)
print("softmax_naive :", softmax_naive(small))
print("softmax :", p_small, " <- the same answer")
print()
print("where a float64 gives up:")
print(" math.exp(709) =", math.exp(709))
try:
math.exp(710)
print(" math.exp(710) =", math.exp(710))
except OverflowError as error:
print(" math.exp(710) -> OverflowError:", error)
print()
print("the same gaps at a scale a real model reaches:", big)
try:
softmax_naive(big)
overflowed = False
except OverflowError as error:
overflowed = True
print(" softmax_naive -> OverflowError:", error)
p_big = softmax(big)
print(" softmax ->", p_big, " <- identical, and it never built exp(1000)")
test("softmax output is a probability distribution", close(sum(p_small), 1.0), True)
test("every probability lands strictly between zero and one", all(0.0 < x < 1.0 for x in p_small), True)
test("the largest logit receives the largest probability", p_small.index(max(p_small)), 0)
test("the naive formula overflows on logits a real model produces", overflowed, True)
test("the shifted version copes with exactly the same input", close(sum(p_big), 1.0), True)
test("and returns the same answer, because only the gaps between logits matter",
all(close(x, y) for x, y in zip(p_small, p_big)), True)
test("adding any constant to every logit leaves softmax unchanged",
all(close(x, y) for x, y in zip(p_small, softmax([z + 5.0 for z in small]))), True)What it printed, when the build ran it
logits : [2.0, 1.0, 0.0] softmax_naive : [0.6652409557748219, 0.24472847105479764, 0.09003057317038046] softmax : [0.6652409557748218, 0.24472847105479764, 0.09003057317038046] <- the same answer where a float64 gives up: math.exp(709) = 8.218407461554972e+307 math.exp(710) -> OverflowError: math range error the same gaps at a scale a real model reaches: [1000.0, 999.0, 998.0] softmax_naive -> OverflowError: math range error softmax -> [0.6652409557748218, 0.24472847105479764, 0.09003057317038046] <- identical, and it never built exp(1000) ok softmax output is a probability distribution ok every probability lands strictly between zero and one ok the largest logit receives the largest probability ok the naive formula overflows on logits a real model produces ok the shifted version copes with exactly the same input ok and returns the same answer, because only the gaps between logits matter ok adding any constant to every logit leaves softmax unchanged
The last two assertions are the whole justification, so they are worth stating in words.
Softmax depends only on the differences between logits: multiply top and bottom of
the fraction by exp(-c) and the c cancels. So you may subtract any
constant you like from every logit before exponentiating. Subtracting the maximum is the
useful choice, because it makes the largest exponent exactly exp(0) = 1 and every
other one smaller. Nothing can overflow, and the denominator is at least 1, so nothing can
divide by zero either. That is why every real implementation does it, and why the naive
three-line version appears in tutorials but never in libraries.
The other half of the story is underflow, and shifting handles that too. A very
negative logit becomes exp(-800), which is not an error — it quietly rounds to
0.0. A zero probability is fine until something takes its logarithm, which is
exactly what a cross-entropy loss does. That is why training code computes
log-softmax directly, as z - shift - log(sum(exp(z - shift))),
rather than taking the log of a softmax that has already rounded to zero.
Temperature: one number that decides how bold the model is
Temperature divides the logits before the softmax. Dividing by a small number spreads the logits apart, which sharpens the distribution; dividing by a large number squashes them together, which flattens it. That is the entire mechanism.
# Temperature: divide the logits, then softmax. That is all it is.
import math
def softmax(logits, temperature=1.0):
if temperature <= 0.0:
raise ValueError("temperature must be positive; zero means argmax, not sampling")
scaled = [z / temperature for z in logits]
shift = max(scaled)
exps = [math.exp(z - shift) for z in scaled]
total = sum(exps)
return [e / total for e in exps]
def entropy(p):
"""How undecided the distribution is. Zero means certain."""
return -sum(x * math.log(x) for x in p if x > 0.0)
TOKENS = ("cat", "dog", "bird", "quark")
LOGITS = [3.0, 1.0, 0.2, -1.0]
TEMPERATURES = (0.25, 0.5, 1.0, 2.0, 20.0)
rows = [(t, softmax(LOGITS, t)) for t in TEMPERATURES]
print("logits:", dict(zip(TOKENS, LOGITS)))
print()
print(" temp " + "".join("{:>9}".format(name) for name in TOKENS) + " entropy spread")
for t, p in rows:
print(" {:5} ".format(t)
+ "".join("{:9.4f}".format(x) for x in p)
+ " {:7.4f} {:.4f}".format(entropy(p), max(p) - min(p)))
print()
print("At 0.25 it says 'cat' essentially every time. At 20 it is close to a coin")
print("toss between four options, including the one it thinks is least likely.")
entropies = [entropy(p) for _, p in rows]
coldest = rows[0][1]
hottest = rows[-1][1]
try:
softmax(LOGITS, 0.0)
refused = False
except ValueError as error:
refused = True
print()
print("softmax(logits, 0.0) ->", error)
test("every temperature still produces a distribution", all(close(sum(p), 1.0) for _, p in rows), True)
test("a cold temperature puts almost all the mass on one token", max(coldest) > 0.99, True)
test("a hot temperature flattens the distribution towards uniform", max(hottest) - min(hottest) < 0.1, True)
test("entropy rises with every step up in temperature", entropies, sorted(entropies))
test("temperature never changes which token is most likely",
all(p.index(max(p)) == 0 for _, p in rows), True)
test("a temperature of zero is refused, because that is argmax and not sampling", refused, True)What it printed, when the build ran it
logits: {'cat': 3.0, 'dog': 1.0, 'bird': 0.2, 'quark': -1.0}
temp cat dog bird quark entropy spread
0.25 0.9997 0.0003 0.0000 0.0000 0.0032 0.9997
0.5 0.9781 0.0179 0.0036 0.0003 0.1166 0.9778
1.0 0.8234 0.1114 0.0501 0.0151 0.6177 0.8083
2.0 0.5715 0.2102 0.1409 0.0773 1.1217 0.4941
20.0 0.2783 0.2518 0.2420 0.2279 1.3836 0.0505
At 0.25 it says 'cat' essentially every time. At 20 it is close to a coin
toss between four options, including the one it thinks is least likely.
softmax(logits, 0.0) -> temperature must be positive; zero means argmax, not sampling
ok every temperature still produces a distribution
ok a cold temperature puts almost all the mass on one token
ok a hot temperature flattens the distribution towards uniform
ok entropy rises with every step up in temperature
ok temperature never changes which token is most likely
ok a temperature of zero is refused, because that is argmax and not sampling
Two things in that table are worth carrying into an interview. First, temperature never
reorders the tokens — the most likely token is the most likely at every temperature, so
"raise the temperature so it picks a different word" is only true in the sense that it makes
the other words reachable. Second, the API parameter called
temperature=0 is not this function with zero passed in; that would divide by
zero. It is a separate code path that skips sampling and takes the argmax.
Sampling, and making it reproducible
Sampling from a distribution needs one random number. Walk along the probabilities, keeping a running total, and stop when the total passes your random draw. This is called inverse transform sampling and it is four lines.
# Sampling: turn a distribution into a choice, reproducibly.
import random
def sample(probs, rng):
"""Inverse transform sampling: walk the running total until it passes r."""
r = rng.random()
running = 0.0
for index, p in enumerate(probs):
running += p
if r < running:
return index
# Reachable only through floating point. The running total can finish a hair
# under 1.0 while r sits just above it, so the last index is the fallback.
# Code that returns None here fails once in millions of tokens, at 3am.
return len(probs) - 1
def greedy(probs):
"""No dice at all. This is what temperature zero actually does."""
return probs.index(max(probs))
TOKENS = ("cat", "dog", "bird", "quark")
PROBS = [0.6, 0.3, 0.08, 0.02]
def draw(seed, n):
rng = random.Random(seed) # our own generator, our own seed, nothing global
return [sample(PROBS, rng) for _ in range(n)]
run_a = draw(42, 12)
run_b = draw(42, 12)
run_c = draw(7, 12)
print("seed 42 :", [TOKENS[i] for i in run_a])
print("seed 42 :", [TOKENS[i] for i in run_b], "<- identical, because the seed is")
print("seed 7 :", [TOKENS[i] for i in run_c], "<- a different stream, a different sentence")
N = 20000
draws = draw(42, N)
counts = [draws.count(i) for i in range(len(PROBS))]
freqs = [c / N for c in counts]
print()
print(" token asked for got in {} draws off by".format(N))
for name, p, f in zip(TOKENS, PROBS, freqs):
print(" {:8} {:6.4f} {:6.4f} {:.4f}".format(name, p, f, abs(f - p)))
worst = max(abs(f - p) for f, p in zip(freqs, PROBS))
print()
print("greedy always returns '{}': no seed, no variety, no surprises.".format(TOKENS[greedy(PROBS)]))
test("the same seed replays exactly the same draws", run_a, run_b)
test("a different seed produces a different sequence", run_a == run_c, False)
test("over twenty thousand draws every frequency lands within a point of its probability",
worst < 0.01, True)
test("even the two per cent token turns up, which is the point of sampling", counts[3] > 0, True)
test("the observed frequencies are themselves a distribution", close(sum(freqs), 1.0), True)
test("greedy decoding ignores the dice and picks the same token every time",
sorted(set(greedy(PROBS) for _ in range(5))), [0])What it printed, when the build ran it
seed 42 : ['dog', 'cat', 'cat', 'cat', 'dog', 'dog', 'dog', 'cat', 'cat', 'cat', 'cat', 'cat'] seed 42 : ['dog', 'cat', 'cat', 'cat', 'dog', 'dog', 'dog', 'cat', 'cat', 'cat', 'cat', 'cat'] <- identical, because the seed is seed 7 : ['cat', 'cat', 'dog', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat', 'cat'] <- a different stream, a different sentence token asked for got in 20000 draws off by cat 0.6000 0.5962 0.0038 dog 0.3000 0.3036 0.0036 bird 0.0800 0.0802 0.0002 quark 0.0200 0.0199 0.0001 greedy always returns 'cat': no seed, no variety, no surprises. ok the same seed replays exactly the same draws ok a different seed produces a different sequence ok over twenty thousand draws every frequency lands within a point of its probability ok even the two per cent token turns up, which is the point of sampling ok the observed frequencies are themselves a distribution ok greedy decoding ignores the dice and picks the same token every time
The reproducibility trap. Note that this code creates
random.Random(42) rather than calling random.seed(42). The
module-level functions share one hidden generator with every library in the process, so any
dependency that draws a number — a retry jitter, a shuffled batch, a cache eviction — moves
your stream and your "reproducible" run stops reproducing. Own your generator, pass it in
as an argument, and a failing generation is something you can replay tomorrow. This is also
the honest answer to "why did the model give a different answer to the same prompt": at
temperature above zero it is sampling, and unless the seed is pinned it will not repeat.
What a real project would use instead
The arithmetic above is the arithmetic a library performs; the difference is that it happens on a whole vocabulary at once, in compiled code, in float32. A vocabulary of 50,000 tokens per generated token makes a Python loop the wrong tool, but nothing conceptual changes.
# The same three steps with the libraries a real project uses.
# NOT RUN HERE - numpy and torch are not installed.
import numpy as np
import torch
logits = np.random.randn(50257) # one score per token in the vocabulary
# Softmax with the same shift, over the whole vocabulary at once.
shifted = logits - logits.max()
probs = np.exp(shifted) / np.exp(shifted).sum()
# Temperature is still just a division before the softmax.
def softmax(logits, temperature=1.0):
shifted = logits / temperature
shifted = shifted - shifted.max()
e = np.exp(shifted)
return e / e.sum()
# Sampling, seeded so a bad generation can be replayed.
rng = np.random.default_rng(42)
token = rng.choice(len(probs), p=probs)
# In torch, and note log_softmax: the log is computed directly rather than by
# taking the log of a probability that has already underflowed to zero.
t = torch.from_numpy(logits)
probs = torch.softmax(t / 0.7, dim=-1)
token = torch.multinomial(probs, num_samples=1)
logp = torch.log_softmax(t, dim=-1) # what a cross-entropy loss usesnumpy and torch are not installed in this environment, so this panel was never executed. Every operation in it was written out and asserted in plain Python above. Read it for the shape of the API, not as verified output.
Worth knowing before an interview: top-k keeps the k highest probabilities, zeroes the rest and re-normalises. top-p, or nucleus sampling, keeps the smallest set of tokens whose probabilities add up to p, and so adapts its own size — wide where the model is unsure, narrow where it is confident. Both are filters applied to the distribution this lesson built, before the sampling step this lesson wrote. Lesson 19 implements both and shows where each one goes wrong.
7 Gradients, and how a model learns
A gradient is a list of slopes: for each parameter, how much the loss changes if you nudge that parameter up. Once you have it, learning is a two-line loop — step against the gradient, repeat. Everything difficult about training is elsewhere: getting the gradient right, and choosing how big a step to take.
A derivative with no calculus in it
The slope of a function at a point is the rise over the run for a very small run. So nudge the input, see what the output does, divide. Two versions are worth knowing: the forward difference, which is the obvious one, and the central difference, which is barely more work and much more accurate.
# A derivative from nothing but arithmetic: nudge the input, watch the output.
def f(x):
return x * x * x - 2.0 * x + 1.0
def f_prime(x):
"""The analytic derivative, by the power rule. 3x^2 - 2."""
return 3.0 * x * x - 2.0
def forward(fn, x, h):
"""The obvious version. One step forward, divided by the step."""
return (fn(x + h) - fn(x)) / h
def central(fn, x, h):
"""A step each way. Same idea, and the errors mostly cancel."""
return (fn(x + h) - fn(x - h)) / (2.0 * h)
x = 1.5
exact = f_prime(x)
print("f(x) = x^3 - 2x + 1")
print("f'(x) = 3x^2 - 2, so at x = 1.5 the true slope is", exact)
print()
print(" h forward error central error")
rows = []
for h in (1e-1, 1e-3, 1e-5, 1e-8, 1e-12):
fw = forward(f, x, h)
ce = central(f, x, h)
rows.append((h, abs(fw - exact), abs(ce - exact)))
print(" {:7.0e} {:13.9f} {:.3e} {:13.9f} {:.3e}".format(
h, fw, abs(fw - exact), ce, abs(ce - exact)))
best = min(rows, key=lambda row: row[2])
print()
print("The central error falls, bottoms out near h = {:.0e}, then gets worse again.".format(best[0]))
print("Too large a step and the curve is not straight across it. Too small and the")
print("subtraction throws away almost every significant digit it had.")
# What it costs. Two evaluations, per parameter, per gradient.
calls = {"n": 0}
def counted(value):
calls["n"] += 1
return f(value)
central(counted, x, 1e-5)
test("the true slope at 1.5 is 4.75", close(exact, 4.75), True)
test("a central difference with a sensible step agrees to eight decimal places",
close(central(f, x, 1e-5), exact, 1e-8), True)
test("the central difference beats the forward difference at the same step size",
rows[1][2] < rows[1][1], True)
test("shrinking the step past the sweet spot makes the answer worse, not better",
rows[-1][2] > best[2], True)
test("so the best step is neither the largest nor the smallest in the table",
best[0] not in (rows[0][0], rows[-1][0]), True)
test("one numerical derivative costs two evaluations of the function", calls["n"], 2)What it printed, when the build ran it
f(x) = x^3 - 2x + 1
f'(x) = 3x^2 - 2, so at x = 1.5 the true slope is 4.75
h forward error central error
1e-01 5.210000000 4.600e-01 4.760000000 1.000e-02
1e-03 4.754501000 4.501e-03 4.750001000 1.000e-06
1e-05 4.750045000 4.500e-05 4.750000000 1.310e-10
1e-08 4.749999993 6.664e-09 4.749999971 2.887e-08
1e-12 4.750422278 4.223e-04 4.750422278 4.223e-04
The central error falls, bottoms out near h = 1e-05, then gets worse again.
Too large a step and the curve is not straight across it. Too small and the
subtraction throws away almost every significant digit it had.
ok the true slope at 1.5 is 4.75
ok a central difference with a sensible step agrees to eight decimal places
ok the central difference beats the forward difference at the same step size
ok shrinking the step past the sweet spot makes the answer worse, not better
ok so the best step is neither the largest nor the smallest in the table
ok one numerical derivative costs two evaluations of the function
The instinct that a smaller h is a more accurate h is wrong,
and the table shows it. Two errors pull in opposite directions: shrinking the step reduces
the error from pretending the curve is straight, and increases the error from subtracting
two nearly identical floats, which is called catastrophic cancellation. At
h = 1e-12 the answer is worse than at h = 0.1. If you ever write
a gradient check, h around 1e-5 with a central difference is the
sensible default.
The analytic gradient, and checking it against the numerical one
Two evaluations per parameter is fine for one parameter and hopeless for a billion, so real training uses derivatives worked out by hand — or by autodiff, which is the same thing mechanised. Here is a model with two parameters, its loss, and the gradient derived with the chain rule, checked against the numerical version.
# A two-parameter model, its loss, and the gradient - derived, then verified.
XS = [0.0, 1.0, 2.0, 3.0]
YS = [1.0, 3.0, 5.0, 7.0] # exactly y = 2x + 1, so we know the right answer
def predict(w, b, x):
return w * x + b
def loss(w, b):
"""Mean squared error: one number for how wrong the whole model is."""
n = len(XS)
return sum((predict(w, b, x) - y) ** 2 for x, y in zip(XS, YS)) / n
def gradient(w, b):
"""The analytic gradient, by the chain rule.
L = mean((w*x + b - y) ** 2)
dL/dw = mean(2 * (w*x + b - y) * x)
dL/db = mean(2 * (w*x + b - y))
"""
n = len(XS)
errors = [predict(w, b, x) - y for x, y in zip(XS, YS)]
dw = sum(2.0 * e * x for e, x in zip(errors, XS)) / n
db = sum(2.0 * e for e in errors) / n
return dw, db
def numerical_gradient(w, b, h=1e-5):
"""The same gradient with no calculus. A check, not a method."""
dw = (loss(w + h, b) - loss(w - h, b)) / (2.0 * h)
db = (loss(w, b + h) - loss(w, b - h)) / (2.0 * h)
return dw, db
analytic = gradient(0.0, 0.0)
numeric = numerical_gradient(0.0, 0.0)
print("at w = 0, b = 0 the loss is", loss(0.0, 0.0))
print()
print("analytic : dL/dw = {:.10f} dL/db = {:.10f}".format(*analytic))
print("numerical : dL/dw = {:.10f} dL/db = {:.10f}".format(*numeric))
print("difference : {:.2e} {:.2e}".format(
abs(analytic[0] - numeric[0]), abs(analytic[1] - numeric[1])))
print()
print("Both slopes are negative, meaning the loss falls as w and b grow.")
print("So we subtract the gradient. That is the whole of 'step against it'.")
step = 0.05
after = loss(0.0 - step * analytic[0], 0.0 - step * analytic[1])
print()
print("loss before one step : {:.6f}".format(loss(0.0, 0.0)))
print("loss after one step : {:.6f}".format(after))
at_answer = gradient(2.0, 1.0)
print()
print("gradient at the exact answer w = 2, b = 1:", at_answer, "<- flat, nowhere left to go")
test("the starting loss is twenty one", close(loss(0.0, 0.0), 21.0), True)
test("the analytic slope for the weight matches the numerical one",
close(analytic[0], numeric[0], 1e-6), True)
test("the analytic slope for the bias matches the numerical one too",
close(analytic[1], numeric[1], 1e-6), True)
test("both slopes are negative, so both parameters need to grow",
[analytic[0] < 0.0, analytic[1] < 0.0], [True, True])
test("one step against the gradient does reduce the loss", after < loss(0.0, 0.0), True)
test("at the exact answer the gradient is zero and so is the loss",
[close(at_answer[0], 0.0), close(at_answer[1], 0.0), close(loss(2.0, 1.0), 0.0)],
[True, True, True])What it printed, when the build ran it
at w = 0, b = 0 the loss is 21.0 analytic : dL/dw = -17.0000000000 dL/db = -8.0000000000 numerical : dL/dw = -16.9999999997 dL/db = -7.9999999997 difference : 2.88e-10 3.03e-10 Both slopes are negative, meaning the loss falls as w and b grow. So we subtract the gradient. That is the whole of 'step against it'. loss before one step : 21.000000 loss after one step : 7.058750 gradient at the exact answer w = 2, b = 1: (0.0, 0.0) <- flat, nowhere left to go ok the starting loss is twenty one ok the analytic slope for the weight matches the numerical one ok the analytic slope for the bias matches the numerical one too ok both slopes are negative, so both parameters need to grow ok one step against the gradient does reduce the loss ok at the exact answer the gradient is zero and so is the loss
That comparison is the single most useful debugging technique in this part of the field, and it has a name: gradient checking. A wrong analytic gradient is vicious because the model still trains. The loss still falls, just more slowly or to a worse place, and nothing raises. Recomputing one gradient numerically and comparing takes four lines and catches a class of bug that otherwise costs days.
Gradient descent, printed step by step
Now the loop. Start with the parameters wrong, compute the gradient, subtract a small multiple of it, repeat. There is nothing else in it, and the loss below is printed at every step so you can watch it fall.
# Gradient descent in full. Three lines of loop, nothing hidden.
XS = [0.0, 1.0, 2.0, 3.0]
YS = [1.0, 3.0, 5.0, 7.0]
def loss(w, b):
n = len(XS)
return sum((w * x + b - y) ** 2 for x, y in zip(XS, YS)) / n
def gradient(w, b):
n = len(XS)
errors = [w * x + b - y for x, y in zip(XS, YS)]
return (sum(2.0 * e * x for e, x in zip(errors, XS)) / n,
sum(2.0 * e for e in errors) / n)
def descend(lr, steps, w=0.0, b=0.0):
"""The trajectory: (w, b, loss) at the start and after every step."""
trail = [(w, b, loss(w, b))]
for _ in range(steps):
dw, db = gradient(w, b)
w -= lr * dw # against the gradient, so downhill
b -= lr * db
trail.append((w, b, loss(w, b)))
return trail
trail = descend(lr=0.1, steps=200)
losses = [row[2] for row in trail]
print("fitting y = w*x + b to y = 2x + 1, from w = 0, b = 0, learning rate 0.1")
print()
print(" step w b loss")
for step in [0, 1, 2, 3, 4, 5, 6, 7, 20, 50, 100, 200]:
w, b, value = trail[step]
print(" {:4} {:10.6f} {:10.6f} {:16.12f}".format(step, w, b, value))
final_w, final_b, final_loss = trail[-1]
print()
print("finished at w = {:.9f}, b = {:.9f}".format(final_w, final_b))
print("final loss {:.3e}, down from {:.1f}".format(final_loss, losses[0]))
print()
print("The first step does most of the work. The rest is the bias creeping up to 1,")
print("which is slower because x = 0 is the only row that isolates it.")
fell_every_step = all(losses[i + 1] < losses[i] for i in range(len(losses) - 1))
test("the loss fell on every one of the two hundred steps", fell_every_step, True)
test("the first step alone removed most of the loss", losses[1] < losses[0] / 30.0, True)
test("the weight converged on two", close(final_w, 2.0, 1e-6), True)
test("the bias converged on one", close(final_b, 1.0, 1e-6), True)
test("the final loss is effectively zero", final_loss < 1e-12, True)
test("nothing was hard-coded: it began at zero and zero", (trail[0][0], trail[0][1]), (0.0, 0.0))What it printed, when the build ran it
fitting y = w*x + b to y = 2x + 1, from w = 0, b = 0, learning rate 0.1
step w b loss
0 0.000000 0.000000 21.000000000000
1 1.700000 0.800000 0.535000000000
2 1.970000 0.930000 0.014350000000
3 2.012000 0.953000 0.001021000000
4 2.017700 0.958800 0.000606235000
5 2.017670 0.961730 0.000528701350
6 2.016782 0.964083 0.000467477941
7 2.015810 0.966232 0.000413509146
20 2.007125 0.984787 0.000083938972
50 2.001132 0.997584 0.000002117569
100 2.000053 0.999887 0.000000004595
200 2.000000 1.000000 0.000000000000
finished at w = 2.000000114, b = 0.999999756
final loss 2.164e-14, down from 21.0
The first step does most of the work. The rest is the bias creeping up to 1,
which is slower because x = 0 is the only row that isolates it.
ok the loss fell on every one of the two hundred steps
ok the first step alone removed most of the loss
ok the weight converged on two
ok the bias converged on one
ok the final loss is effectively zero
ok nothing was hard-coded: it began at zero and zero
The learning rate, and the cliff it falls off
The step size is the one number you cannot reason your way past by staring at the code. Too small and training crawls. Too large and the parameters overshoot, land further away than they started, and the loss grows instead of shrinking. For a squared loss the exact threshold is computable, which makes this the rare case where you can watch theory and arithmetic agree.
# The learning rate: too small crawls, too large explodes, and the edge is exact.
import math
XS = [0.0, 1.0, 2.0, 3.0]
YS = [1.0, 3.0, 5.0, 7.0]
def loss(w, b):
n = len(XS)
return sum((w * x + b - y) ** 2 for x, y in zip(XS, YS)) / n
def gradient(w, b):
n = len(XS)
errors = [w * x + b - y for x, y in zip(XS, YS)]
return (sum(2.0 * e * x for e, x in zip(errors, XS)) / n,
sum(2.0 * e for e in errors) / n)
def descend(lr, steps, w=0.0, b=0.0):
history = [loss(w, b)]
for _ in range(steps):
dw, db = gradient(w, b)
w -= lr * dw
b -= lr * db
history.append(loss(w, b))
return w, history
RATES = (0.02, 0.05, 0.1, 0.2, 0.23, 0.237, 0.239, 0.24, 0.3)
STEPS = 80
results = {}
print("every run starts from w = 0, b = 0, where the loss is {:.1f}".format(loss(0.0, 0.0)))
print()
print(" rate loss after {} steps w ended at verdict".format(STEPS))
for lr in RATES:
w, history = descend(lr, STEPS)
results[lr] = (w, history)
if history[-1] >= history[0]:
verdict = "DIVERGED"
elif history[-1] > 1e-3:
verdict = "converging, slowly"
else:
verdict = "converged"
print(" {:<6} {:20.6g} {:15.6g} {}".format(lr, history[-1], w, verdict))
# The edge is not a mystery. For a squared loss the curvature is a constant matrix,
# H = (2/n) * [[sum(x*x), sum(x)], [sum(x), n]]
# and plain gradient descent is stable exactly while lr < 2 / largest eigenvalue.
n = len(XS)
sum_x2 = sum(x * x for x in XS)
sum_x = sum(XS)
trace = (2.0 / n) * (sum_x2 + n)
determinant = (2.0 / n) ** 2 * (sum_x2 * n - sum_x * sum_x)
lambda_max = (trace + math.sqrt(trace * trace - 4.0 * determinant)) / 2.0
threshold = 2.0 / lambda_max
converging = [lr for lr in RATES if results[lr][1][-1] < results[lr][1][0]]
diverging = [lr for lr in RATES if results[lr][1][-1] >= results[lr][1][0]]
print()
print("largest curvature : {:.6f}".format(lambda_max))
print("predicted safe limit : {:.6f}".format(threshold))
print("largest rate that worked : {}".format(max(converging)))
print("smallest rate that failed: {}".format(min(diverging)))
test("a very small rate helps but does not finish the job",
[results[0.02][1][-1] < 1.0, results[0.02][1][-1] > results[0.2][1][-1]], [True, True])
test("a well-chosen rate lands on the answer", close(results[0.2][0], 2.0, 1e-5), True)
test("the largest rate that converged is below the predicted limit", max(converging) < threshold, True)
test("the smallest rate that diverged is above it", min(diverging) > threshold, True)
test("so the prediction brackets the observed cliff, to three decimal places",
close(threshold, 0.238, 0.001), True)
test("well past the limit the loss does not merely fail, it explodes",
results[0.3][1][-1] > 1e20, True)What it printed, when the build ran it
every run starts from w = 0, b = 0, where the loss is 21.0 rate loss after 80 steps w ended at verdict 0.02 0.000143781 2.00932 converged 0.05 7.78417e-06 2.00217 converged 0.1 5.34209e-08 2.00018 converged 0.2 1.53991e-12 2 converged 0.23 0.000328527 1.99199 converged 0.237 5.82222 0.934095 converging, slowly 0.239 85.6489 -2.08823 DIVERGED 0.24 323.044 -5.93972 DIVERGED 0.3 3.07201e+30 -7.74258e+14 DIVERGED largest curvature : 8.405125 predicted safe limit : 0.237950 largest rate that worked : 0.237 smallest rate that failed: 0.239 ok a very small rate helps but does not finish the job ok a well-chosen rate lands on the answer ok the largest rate that converged is below the predicted limit ok the smallest rate that diverged is above it ok so the prediction brackets the observed cliff, to three decimal places ok well past the limit the loss does not merely fail, it explodes
Notice how the failure looks. At 0.239 the loss ends four times higher than
it started; at 0.3 it reaches 1e30 and the weights are in the
hundreds of trillions. In a real training run this is the nan that appears in
your loss after two hundred steps, and the first thing to try is always dividing the
learning rate by ten. The threshold depends on the data — it is set by the curvature, which
is set by the scale of your inputs — which is the real reason features get normalised before
training. Rescaling the inputs changes what a safe learning rate even is.
Why nobody computes gradients numerically
The numerical gradient in the second panel cost two evaluations of the loss for each of two parameters. That ratio does not improve. A model with a million parameters needs two million forward passes for one gradient, and a large language model has billions. Reverse-mode automatic differentiation — backpropagation — computes every partial derivative in a single backward pass costing roughly the same as one forward pass, regardless of how many parameters there are. That single fact is why deep learning is possible at all, and lesson 14 works through it by hand.
# The same descent, with the library that does the derivative for you.
# NOT RUN HERE - torch is not installed.
import torch
xs = torch.tensor([0.0, 1.0, 2.0, 3.0])
ys = torch.tensor([1.0, 3.0, 5.0, 7.0])
w = torch.zeros(1, requires_grad=True) # requires_grad: track operations on this
b = torch.zeros(1, requires_grad=True)
optimiser = torch.optim.SGD([w, b], lr=0.1) # the two lines of loop, packaged
for step in range(200):
loss = ((w * xs + b - ys) ** 2).mean() # the same mean squared error
optimiser.zero_grad() # forget last step's gradients, or they add up
loss.backward() # every partial derivative, one backward pass
optimiser.step() # w -= lr * w.grad, and the same for b
print(w.item(), b.item()) # 2.0, 1.0
# torch.autograd.gradcheck compares the analytic gradient with a numerical one,
# which is the technique from the second panel, packaged and made rigorous.torch is not installed in this environment, so this panel was never executed. The plain-Python loop above computes the identical gradient and did run, with its loss printed at every step. Read this for the shape of the API, not as verified output.
What to take from this lesson, because it is the last of the foundations. A loss is one number saying how wrong you are. A gradient is one slope per parameter. Training subtracts a small multiple of the gradient, over and over. The two things that go wrong are a gradient that is subtly incorrect, which you catch by checking it numerically, and a step size on the wrong side of a cliff, which you catch by watching the loss. Everything from here on — linear regression, logistic regression, a neural network, a transformer — is that same loop with a more complicated function in the middle.