You go to the market and buy 3 apples; with parking, the trip costs 6. Next week you buy 5 apples; the trip costs 8. How much will a 10-apple trip cost?
Write the story as an equation:
$$Expense = parking + price \cdot apples \qquad\text{i.e.}\qquad y = \beta_1 x + \beta_0$$
Trip 1 says $3\beta_1 + \beta_0 = 6$; trip 2 says $5\beta_1 + \beta_0 = 8$. Subtract them: $2\beta_1 = 2$, so the price is 1 and parking is 3. Ten apples: 13.
That is the oldest prediction machine there is: turn observations into equations, solve, extrapolate. Long before gradient descent, this is how curves were fit --- and in this chapter you rebuild the whole machine, then push it to its breaking point (data that no curve passes through) and invent least squares.
Prerequisite: the Loops and Arrays chapter --- you will reuse your solve_equations from its Part 5.
A fitted model is just a list of coefficients --- a vector. Evaluating it at a point is a dot product. Build that toolkit first.
dot_product(a, b) multiplies element-by-element and sums: dot_product([1, 2], [3, 4]) is 1*3 + 2*4 = 11.
Tasks:
dot_product(a, b).assert dot_product([1, 2], [3, 4]) == 11
assert dot_product([1, 0, 2], [5, 5, 5]) == 15
assert dot_product([2], [7]) == 14
assert dot_product([1, 2], [3, 4]) == 11
assert dot_product([1, 0, 2], [5, 5, 5]) == 15
assert dot_product([2], [7]) == 14
Use zip to pair up elements from both lists, multiply each pair, and sum the results.
return sum(ai * bi for ai, bi in zip(a, b))
def dot_product(a, b):
return sum(ai * bi for ai, bi in zip(a, b))
The length of a vector is the distance formula in disguise: $\|v\| = \sqrt{v \cdot v}$ --- reuse your dot_product.
A unit vector keeps the direction but has length 1: divide every component by the length. convert_to_unit([3, 4]) -> [0.6, 0.8].
Tasks:
vector_length(v) using dot_product and math.sqrt.convert_to_unit(v) that divides each component by the length.assert vector_length([3, 4]) == 5.0
assert convert_to_unit([3, 4]) == [0.6, 0.8]
assert abs(vector_length(convert_to_unit([2, 3, 6])) - 1.0) < 1e-9
assert vector_length([3, 4]) == 5.0
assert convert_to_unit([3, 4]) == [0.6, 0.8]
assert abs(vector_length(convert_to_unit([2, 3, 6])) - 1.0) < 1e-9
For vector_length: the dot product of a vector with itself gives the sum of squares. Take the square root of that.
vector_length(v) = math.sqrt(dot_product(v, v)). For convert_to_unit(v): compute L = vector_length(v), then return [x / L for x in v].
import math
def vector_length(v):
return math.sqrt(dot_product(v, v))
def convert_to_unit(v):
L = vector_length(v)
return [x / L for x in v]
How aligned are two vectors? Take the dot product of their unit versions. The result is the cosine of the angle between them:
1.00.0-1.0This one number is how modern systems compare documents, faces and embeddings --- you will meet it again the moment you touch LLMs.
Tasks:
cosine_similarity(a, b) using convert_to_unit and dot_product.assert abs(cosine_similarity([1, 2], [2, 4]) - 1.0) < 1e-9 # same direction
assert abs(cosine_similarity([1, 0], [0, 1])) < 1e-9 # perpendicular
assert abs(cosine_similarity([1, 2], [-1, -2]) + 1.0) < 1e-9 # opposite
assert abs(cosine_similarity([1, 2], [2, 4]) - 1.0) < 1e-9 # same direction
assert abs(cosine_similarity([1, 0], [0, 1])) < 1e-9 # perpendicular
assert abs(cosine_similarity([1, 2], [-1, -2]) + 1.0) < 1e-9 # opposite
Convert both vectors to unit vectors first, then take their dot product. The result is the cosine of the angle between them.
return dot_product(convert_to_unit(a), convert_to_unit(b))
def cosine_similarity(a, b):
return dot_product(convert_to_unit(a), convert_to_unit(b))
A matrix is a list of rows. In the product C = A @ B, entry C[i][j] is the dot product of row i of A with column j of B. An (n x k) matrix times a (k x m) matrix gives an (n x m) matrix.
Write matrix_mul(A, B) using your dot_product.
Hint: column j of B is [row[j] for row in B].
Verify:
A = [[1, 2],
[3, 4]]
B = [[5, 6],
[7, 8]]
assert matrix_mul(A, B) == [[19, 22], [43, 50]]
M = [[1, 2, 3],
[4, 5, 6]]
v = [[1], [0], [2]] # a column vector
assert matrix_mul(M, v) == [[7], [16]] # (2x3) @ (3x1) -> (2x1)
A = [[1, 2],
[3, 4]]
B = [[5, 6],
[7, 8]]
assert matrix_mul(A, B) == [[19, 22], [43, 50]]
M = [[1, 2, 3],
[4, 5, 6]]
v = [[1], [0], [2]] # a column vector
assert matrix_mul(M, v) == [[7], [16]] # (2x3) @ (3x1) -> (2x1)
First figure out the dimensions: A has len(A) rows and len(A[0]) columns; B has len(B[0]) columns. The result has len(A) rows and len(B[0]) columns.
For each row i in A and each column j in B, compute dot_product(A[i], [row[j] for row in B]). Collect these into the result matrix.
def matrix_mul(A, B):
n_cols_B = len(B[0])
result = []
for row in A:
new_row = []
for j in range(n_cols_B):
col = [r[j] for r in B]
new_row.append(dot_product(row, col))
result.append(new_row)
return result
Two vectors have a cosine similarity of 0. What does this tell you about their relationship?
Cosine similarity measures the cosine of the angle between two vectors. What angle has cosine = 0?
cos(90 degrees) = 0, so cosine similarity of 0 means the vectors are perpendicular --- completely unrelated in direction. Identical vectors have cosine similarity of 1 (angle = 0 degrees). Opposite vectors have cosine similarity of -1 (angle = 180 degrees). A zero vector would make the computation undefined (division by zero in the denominator), not produce a result of 0.
We store the equation $a_1 x_1 + a_2 x_2 + \dots + a_n x_n = b$ as the list [a1, a2, ..., an, b], and a system as a list of such lists --- exactly the representation you used in Loops and Arrays, Part 5.
Paste your solve_equations (with its helpers solve_for_first_variable, eliminate_variable_pair, eliminate_variable) into the code cell. The test must pass before you go on.
(Skipped that chapter? Go invent it there first --- this whole chapter stands on it.)
Think: your solver eliminates one variable at a time, and each elimination touches every remaining equation, and each equation has that many coefficients --- three nested loops. Roughly how does the running time grow with n equations? If solving 4 equations takes 1 second, what would 8 take? Keep that number in mind --- it returns at the end of the chapter.
r = solve_equations([[1, 1, 1, 35], [3, 2, 1, 75], [4, 3, 1, 105]])
assert all(abs(a - b) < 1e-9 for a, b in zip(r, [10.0, 20.0, 5.0])), r
print("solver ready:", r)
The elimination algorithm has three nested loops, each of size roughly n. So the running time grows like $n^3$. If 4 equations take 1 second, 8 equations take roughly $8^3 / 4^3 = 8$ seconds.
The running time grows like $n^3$. If 4 equations take 1 second, 8 equations take roughly $8^3 / 4^3 = 8$ seconds.
def solve_for_first_variable(equation, vars):
rhs = equation[-1]
for i, v in enumerate(vars):
rhs -= equation[i + 1] * v
return rhs / equation[0]
def eliminate_variable_pair(eq1, eq2, var_index):
c1 = eq1[var_index]
c2 = eq2[var_index]
new_eq = []
for i in range(len(eq1)):
if i == var_index:
continue
new_eq.append(c2 * eq1[i] - c1 * eq2[i])
return new_eq
def eliminate_variable(equations, var_index):
pivot = equations[0]
result = []
for eq in equations[1:]:
result.append(eliminate_variable_pair(pivot, eq, var_index))
return result
def solve_equations(equations):
if len(equations) == 1:
return [equations[0][-1] / equations[0][0]]
reduced = eliminate_variable(equations, 0)
rest = solve_equations(reduced)
x1 = solve_for_first_variable(equations[0], rest)
return [x1] + rest
Redo the opening story with your tools:
(beta1, beta0): trip 1 -> [3, 1, 6], trip 2 -> [5, 1, 8].solve_equations gives you beta = [price, parking].dot_product(beta, [10, 1]).# YOUR CODE — set up the equations, solve, predict the 10-apple trip
eqns = ...
beta = ...
prediction = ...
assert all(abs(a - b) < 1e-9 for a, b in zip(beta, [1.0, 3.0])), beta
assert abs(prediction - 13.0) < 1e-9
print("a 10-apple trip will cost:", prediction)
The equation for trip 1 (3 apples, cost 6) is [3, 1, 6] meaning 3*price + 1*parking = 6. Similarly for trip 2.
eqns = [[3, 1, 6], [5, 1, 8]]. Then beta = solve_equations(eqns) and prediction = dot_product(beta, [10, 1]).
eqns = [[3, 1, 6], [5, 1, 8]]
beta = solve_equations(eqns)
prediction = dot_product(beta, [10, 1])
Traffic through a junction is measured three times:
| t | traffic |
|---|---|
| 1 | 6 |
| 2 | 12 |
| 3 | 20 |
The growth is not linear (check the gaps: +6, +8). Assume a quadratic: $traffic = a t^2 + b t + c$. Each measurement is one equation in $(a, b, c)$ --- for $t=1$: $a + b + c = 6$, and so on.
Tasks:
(a, b, c).t = 4.# YOUR CODE — three equations in (a, b, c), solve, predict t = 4
eqns = ...
abc = ...
prediction = ...
assert all(abs(x - y) < 1e-9 for x, y in zip(abc, [1.0, 3.0, 2.0])), abc
assert abs(prediction - 30.0) < 1e-9
print("traffic at t=4:", prediction)
For $t = 1$: $a(1)^2 + b(1) + c = 6$ -> [1, 1, 1, 6]. For $t = 2$: $a(4) + b(2) + c = 12$ -> [4, 2, 1, 12]. For $t = 3$: [9, 3, 1, 20].
To predict at $t = 4$: compute dot_product(abc, [16, 4, 1]) since $a \cdot 16 + b \cdot 4 + c$.
eqns = [[1, 1, 1, 6], [4, 2, 1, 12], [9, 3, 1, 20]]
abc = solve_equations(eqns)
prediction = dot_product(abc, [16, 4, 1])
Automate it. Write solve_poly(x, y, xp):
x and y are n observed points; assume a polynomial of degree n - 1.[x_i**(n-1), ..., x_i**2, x_i, 1, y_i] and solve for the coefficients.xp --- a dot_product of the coefficients with the matching powers of xp.Think: why do n points pin down exactly a degree-(n-1) polynomial --- what goes wrong with fewer coefficients, and what freedom is left with more?
assert abs(solve_poly([1, 2, 3], [6, 12, 20], 4) - 30.0) < 1e-6
# a line is just a degree-1 polynomial — the apple trips again:
assert abs(solve_poly([3, 5], [6, 8], 10) - 13.0) < 1e-6
# and it must reproduce the points it was built from:
assert abs(solve_poly([1, 2, 3], [6, 12, 20], 2) - 12.0) < 1e-6
print("solve_poly: OK")
For n points you need a polynomial with n coefficients (degree n-1). Each point gives one equation: the row is the powers of x_i from highest to lowest, ending with 1, then the y_i value.
Build each row as: [x_i**(n-1-j) for j in range(n)] + [y_i]. Solve for coefficients coefs, then return dot_product(coefs, [xp**(n-1-j) for j in range(n)]).
def solve_poly(x, y, xp):
n = len(x)
eqns = []
for i in range(n):
row = [x[i] ** (n - 1 - j) for j in range(n)] + [y[i]]
eqns.append(row)
coefs = solve_equations(eqns)
powers = [xp ** (n - 1 - j) for j in range(n)]
return dot_product(coefs, powers)
A factory reports:
| input | output |
|---|---|
| 1 | 8 |
| 2 | 11 |
| 3 | 16 |
| 4 | ?? |
Try a line through the first two points: slope 3, intercept 5 --- but then input 3 predicts 14, not 16. Through points 1 and 3: slope 4, intercept 4 --- but input 2 predicts 12, not 11. No line passes through all three points. Real measurements are noisy; perfect solving just broke.
New plan: accept an imperfect line, but pick the one that is least wrong. But how do you measure how wrong a line is? You'll figure that out in the next exercise.
For a candidate line $y = mx + c$, each data point $(x_i, y_i)$ has an error: the prediction $mx_i + c$ minus the actual $y_i$. But some errors are positive (prediction too high) and some are negative (prediction too low) --- if you just add them up, they cancel out.
Think: How do you prevent cancellation? What operation makes all errors positive? Then how do you combine them into a single "total wrongness" score?
Write total_error(m, c, data) where data is a list of (x, y) pairs.
factory_data = [(1, 8), (2, 11), (3, 16)]
assert total_error(3, 5, factory_data) == 4.0 # hits two points, misses the third by 2
assert total_error(0, 0, factory_data) == 8*8 + 11*11 + 16*16
print("total_error: OK")
Square each error to make it positive (and to penalize large errors more). The total squared error is: $$E(m, c) = \sum_i (m x_i + c - y_i)^2$$ For each (x, y) pair, compute m * x + c - y, square it, and sum.
return sum((m * x + c - y) ** 2 for x, y in data)
def total_error(m, c, data):
return sum((m * x + c - y) ** 2 for x, y in data)
At the minimum of $E$, nudging $m$ or $c$ changes nothing --- both partial derivatives are zero.
Your task:
solve_equations from Part 2 to find $m$ and $c$.Use solve_equations on the two equations you derived. The test also nudges your answer in every direction to confirm no neighbor does better --- the definition of a minimum.
# YOUR CODE — solve the two normal equations for m and c
m, c = ...
assert abs(m - 4.0) < 1e-9, m
assert abs(c - 11/3) < 1e-9, c
best = total_error(m, c, factory_data)
for dm, dc in [(0.01, 0), (-0.01, 0), (0, 0.01), (0, -0.01)]:
assert total_error(m + dm, c + dc, factory_data) > best
print(f"best line: y = {m:.3f}x + {c:.3f} (error {best:.4f})")
After setting each partial derivative to zero, you get two linear equations in $m$ and $c$. Compute the sums ($\sum x_i^2$, $\sum x_i$, $\sum x_i y_i$, $\sum y_i$, $n$) from the factory data and plug them in, then pass the two equations to solve_equations.
The partial derivatives of $E = \sum_i (mx_i + c - y_i)^2$:
$\frac{\partial E}{\partial m} = 2\sum_i x_i(mx_i + c - y_i) = 0$ gives: $m\sum x_i^2 + c\sum x_i = \sum x_i y_i$
$\frac{\partial E}{\partial c} = 2\sum_i (mx_i + c - y_i) = 0$ gives: $m\sum x_i + cn = \sum y_i$
factory_data = [(1, 8), (2, 11), (3, 16)]
sum_x2 = sum(x ** 2 for x, y in factory_data) # 14
sum_x = sum(x for x, y in factory_data) # 6
sum_xy = sum(x * y for x, y in factory_data) # 78
sum_y = sum(y for x, y in factory_data) # 35
n = len(factory_data) # 3
m, c = solve_equations([
[sum_x2, sum_x, sum_xy],
[sum_x, n, sum_y]
])
Wrap it up. Write best_line(data):
(x, y) pairs compute $\sum x^2$, $\sum x$, $\sum xy$, $\sum y$ and $n$.solve_equations.(m, c).On data that is perfectly linear it must recover the exact line --- least squares contains exact fitting as a special case.
m, c = best_line([(0, 3), (1, 5), (2, 7)]) # perfectly linear: y = 2x + 3
assert abs(m - 2.0) < 1e-9 and abs(c - 3.0) < 1e-9
m, c = best_line(factory_data)
assert abs(m - 4.0) < 1e-9 and abs(c - 11/3) < 1e-9
print("factory prediction for input 4:", round(m * 4 + c, 3))
print("best_line: OK")
import matplotlib.pyplot as plt
xs = [p[0] for p in factory_data]
ys = [p[1] for p in factory_data]
m, c = best_line(factory_data)
line_x = [0.5, 4.2]
line_y = [m * x + c for x in line_x]
plt.scatter(xs, ys, s=70, zorder=3, label="measurements")
plt.plot(line_x, line_y, "r-", label=f"y = {m:.2f}x + {c:.2f}")
plt.scatter([4], [m * 4 + c], marker="*", s=200, color="green", zorder=3,
label="prediction for input 4")
plt.xlabel("input"); plt.ylabel("output")
plt.title("Least squares: the least-wrong line")
plt.legend(); plt.grid(alpha=0.3)
plt.show()
Compute the five sums: sum_x2 = sum(x**2 for x, y in data), sum_x = sum(x for x, y in data), sum_xy = sum(x*y for x, y in data), sum_y = sum(y for x, y in data), and n = len(data).
Build [[sum_x2, sum_x, sum_xy], [sum_x, n, sum_y]] and pass to solve_equations. The result is [m, c].
def best_line(data):
sum_x2 = sum(x ** 2 for x, y in data)
sum_x = sum(x for x, y in data)
sum_xy = sum(x * y for x, y in data)
sum_y = sum(y for x, y in data)
n = len(data)
m, c = solve_equations([
[sum_x2, sum_x, sum_xy],
[sum_x, n, sum_y]
])
return (m, c)
In the least squares method, you minimize the sum of squared errors. Why not just minimize the sum of raw errors (without squaring)?
If actual = 20 and predicted = 30, the error is -10. If actual = 40 and predicted = 30, the error is +10. What's their sum?
Raw errors cancel: an overestimate of +10 and an underestimate of -10 sum to 0, suggesting 'no error' when both predictions are clearly wrong. Squaring makes all errors positive ((-10)^2 = 100, (+10)^2 = 100), so they can't cancel. It also penalizes large errors more heavily than small ones. Absolute values would also prevent cancellation (that's MAE), but squaring has a mathematical advantage: the squared-error function is smooth and differentiable everywhere, making it easier to minimize with calculus.
So far matrices held data (trips, prices). But matrix_mul has a second life: a matrix can be a machine that moves points. Every video game, every 3D animation, every graphics chip is built on this.
Keep points as rows, [x, y]. Multiplying a stack of points by one 2x2 matrix moves all of them at once with your existing matrix_mul --- no new code needed.
The matrix that rotates every point counter-clockwise by angle $\theta$ (with points as rows) is:
$$R(\theta) = \begin{pmatrix} \cos\theta & \sin\theta \ -\sin\theta & \cos\theta \end{pmatrix}$$
Write rotation_matrix(angle_degrees) returning that 2x2 list-of-lists. Use math.radians to convert, then math.cos / math.sin.
Sanity checks worth predicting before you run: the point [1, 0] (due east) rotated 90 degrees should land at [0, 1] (due north). Rotated 180 degrees, at [-1, 0].
p = matrix_mul([[1, 0]], rotation_matrix(90))[0]
assert abs(p[0]) < 1e-9 and abs(p[1] - 1) < 1e-9 # east -> north
p = matrix_mul([[1, 0]], rotation_matrix(180))[0]
assert abs(p[0] + 1) < 1e-9 and abs(p[1]) < 1e-9 # east -> west
p = matrix_mul([[0, 1]], rotation_matrix(90))[0]
assert abs(p[0] + 1) < 1e-9 and abs(p[1]) < 1e-9 # north -> west
print("rotation_matrix: OK")
Convert degrees to radians with math.radians(angle_degrees). Then build the 2x2 list using math.cos(rad) and math.sin(rad).
rad = math.radians(angle_degrees), then return [[math.cos(rad), math.sin(rad)], [-math.sin(rad), math.cos(rad)]].
import math
def rotation_matrix(angle_degrees):
rad = math.radians(angle_degrees)
return [[math.cos(rad), math.sin(rad)],
[-math.sin(rad), math.cos(rad)]]
Here is a kite, five points as rows (the last repeats the first so it draws closed):
kite = [[0, 0], [1, 2], [0, 5], [-1, 2], [0, 0]]
Rotate the entire kite by 45 degrees with one call to matrix_mul, storing the result in turned. Then run the provided plotting cell to see both.
A deep property to verify while you're here: rotation moves points but never stretches them. The test checks every point stays the same distance from the origin.
# Tests: same number of points, all distances to origin preserved
assert len(turned) == len(kite)
for before, after in zip(kite, turned):
d_before = (before[0] ** 2 + before[1] ** 2) ** 0.5
d_after = (after[0] ** 2 + after[1] ** 2) ** 0.5
assert abs(d_before - d_after) < 1e-9
print("rotation preserves all lengths — the shape is rigid")
import matplotlib.pyplot as plt
plt.plot([p[0] for p in kite], [p[1] for p in kite], "b-o", label="original")
plt.plot([p[0] for p in turned], [p[1] for p in turned], "r-o", label="rotated 45 degrees")
plt.axhline(0, color="gray", lw=0.5); plt.axvline(0, color="gray", lw=0.5)
plt.gca().set_aspect("equal")
plt.legend(); plt.grid(alpha=0.3)
plt.title("One matrix multiplication moves every point")
plt.show()
The kite is already a matrix (5x2). Multiply it by the 2x2 rotation matrix: turned = matrix_mul(kite, rotation_matrix(45)).
kite = [[0, 0], [1, 2], [0, 5], [-1, 2], [0, 0]]
turned = matrix_mul(kite, rotation_matrix(45))
Rotation matrices are special. What happens with an arbitrary 2x2 matrix?
Build one from random numbers, apply it to the kite with matrix_mul, and plot the result (reuse the plotting code above). Run it several times.
Observe, then answer: straight edges stay straight, and parallel edges stay parallel --- but lengths and angles change: the kite gets stretched, squashed, sheared, maybe flipped. Every 2x2 matrix is some such linear transformation, and rotations are exactly the ones that preserve all lengths. (When a neural network multiplies data by a learned weight matrix, this stretching and squashing of space is literally what it's doing.)
import random
R = [[random.uniform(-1, 1), random.uniform(-1, 1)],
[random.uniform(-1, 1), random.uniform(-1, 1)]]
warped = ... # replace this — apply R to the kite
# Test: still 5 points, each with 2 coordinates
assert len(warped) == 5 and all(len(p) == 2 for p in warped)
plt.plot([p[0] for p in kite], [p[1] for p in kite], "b-o", label="original")
plt.plot([p[0] for p in warped], [p[1] for p in warped], "g-o", label="warped")
plt.gca().set_aspect("equal")
plt.legend(); plt.grid(alpha=0.3)
plt.title("A random matrix: straight lines stay straight, nothing else survives")
plt.show()
Same pattern as the rotation: warped = matrix_mul(kite, R). The random matrix R is already 2x2.
Every 2x2 matrix applies a linear transformation: straight lines stay straight and parallel lines stay parallel, but lengths and angles generally change. Rotation matrices are the special subset that preserve all lengths (and therefore all angles too).
warped = matrix_mul(kite, R)
| You built | The secret |
|---|---|
dot_product, vector_length, convert_to_unit, cosine_similarity, matrix_mul |
models are vectors; evaluation is a dot product |
| apple trips, traffic at the junction | prediction = set up equations + solve |
solve_poly |
n points <-> one degree-(n-1) polynomial |
total_error, normal equations, best_line |
least squares --- when nothing fits, minimize the squared error |
What you just derived is Ordinary Least Squares, the same closed-form answer sklearn.linear_model.LinearRegression computes.
So why does the next chapter exist? Two cracks in the ancient method: solving grows like $n^3$ (your Part 2 estimate), and the derive-set-to-zero trick needs a formula for the derivative --- which complicated models don't offer. Gradient Descent fixes both by walking toward the minimum instead of jumping to it.