How did I get from towers on a coordinate plane to writing my own regression + optimization engine?

Towers, a house, and a cost nobody had defined yet

Before any math, there was a plain-language problem. I wrote it down first, in my own words, before I knew what tool I'd need to solve it.

The problem, in my own words

The first cell of regression_one_tower_problem.ipynb is not code — it's a docstring where I stated the problem before touching NumPy:

So the original question was about one unknown point — the house — and n known points — the towers. The "cost" of a candidate house location was the sum of squared distances from that single point to every tower. This is a facility-location problem: where do I place one point so that, on average (in a squared sense), it's as close as possible to a set of fixed points?

What information was available, what was unknown

• Known: the coordinates of every tower — an x and a y value each.
• Unknown: a location that minimizes total squared distance to all of them.
• Assumed: that "cost" could be captured entirely by geometric distance, squared.

Why I needed a mathematical model at all

I couldn't just "look" at five towers scattered on a plane and guess the best point — with more than two towers there's no single obvious answer, and "best" only means something once cost is a number I can compare across candidate locations. That's the whole reason a cost function enters the picture at Stage 0, before any regression vocabulary shows up: I needed a way to turn "is this a good location?" into "how small is this number?"

The coordinate interpretation

Every tower became a pair (x, y). That single move — turning a physical tower into a coordinate pair — is what makes the rest of the journey possible. Once towers are points in a plane, any relationship I hypothesize between them (a line, a curve, a location) becomes a question I can answer with algebra instead of intuition.

How the idea of a line arose

Once the towers are treated as (x, y) points that I want to relate to each other, the simplest possible relationship I can hypothesize is a straight line: y = mx + b. It has exactly two unknowns, it's easy to differentiate, and if it fits the towers reasonably well I can use it to reason about the whole configuration — including, eventually, a "best" point along it. This is the geometric idea that Stage 1 formalizes.

→ Next: once towers are coordinates, what's the simplest relationship I can propose between them?

Representing towers as data, and proposing y = mx + b

What I knew before this stage: towers have coordinates. What I didn't know yet: how to write "a line through them" in code.

Representing the tower coordinates

Across the notebooks the same small dataset keeps reappearing — five towers used again and again as the canonical test case, first in regression_multi_tower_problem.ipynb and then reused in optimization_using_newton_method.ipynb:

tx = [2,4,6,8,9]
ty = [3,6,8,2,5]

tx holds the x-coordinate of each tower, ty the y-coordinate. Tower i is the point (tx[i], ty[i]). This is the same 5-point dataset driving the hero widget at the top of this page.

Understanding x and y

It's worth being precise about what x and y even are here, because the notebooks never fully pin this down (a good example of a real gap between "problem statement" and "implementation" — see Stage 0). Whatever x and y represent physically, mathematically I only need two things from them: a value I treat as input (x) and a value I'm trying to relate it to (y).

The model: y = mx + b

This is a hypothesis, not a fact: I'm proposing that y can be approximated as a linear function of x. Two unknowns fully describe this hypothesis:

• m — the slope. Geometrically, how much ŷ changes per unit of x. In the tower context, it's the rate at which the y-coordinate trends upward or downward as x increases across towers.
• b — the intercept. Where the line crosses x = 0. It's the baseline value of ŷ before x has any effect.

I named the residual function fxy in regression_one_tower_problem.ipynb, even though it isn't really "a function of x and y" — it's a function of m, b, and one point:

# function for calculating a single value f(x,y)
# for two variable for a point in place (a1,b1) cost function c(m,b) = (m·a1 + b - b1)**2

def fxy(m,b,a1,b1):
    return m*a1 + b - b1

This naming is worth keeping rather than "fixing," because it captures exactly where my head was at this point: I was still thinking of this as "a function involving x and y," even though what it actually computes — m·a1 + b − b1 — is the residual: how far the line's prediction at a1 is from the true value b1. Naming it after its inputs rather than what it means is a very natural first-pass mistake, and it's the seed of the vocabulary I properly separate out in Stage 2 (error vs. residual vs. squared error vs. cost).

What happens when the line doesn't pass through the points

Almost always, it doesn't — and it can't, in general, pass through all five towers exactly, because five points impose five constraints on a line that only has two degrees of freedom (m and b). For most datasets a line has to miss some or all of the points. Drag m and b in the hero widget above and watch: for any choice, some towers sit above the line and some below. Those vertical gaps are exactly what Stage 2 needs to measure.

→ Next: the line almost never fits perfectly. How do I turn "how far off" into a single number I can minimize?

Error, residual, squared error, and total cost are four different things

What I knew before this stage: a line usually misses the towers. What I didn't know: how to turn "misses" into one number to minimize.

Why a line needs an objective at all

"This line looks pretty good" isn't something an algorithm can act on. I needed a single number that gets smaller as the line gets better and bigger as it gets worse, so that "improving the line" becomes "decreasing a number" — a problem calculus and search algorithms actually know how to solve.

Four words I was using loosely, made precise

• Geometric distance / error for a point: in the facility-location framing from Stage 0, this would be the full Euclidean distance from a candidate point to a tower, √((x−h)² + (y−k)²). This is not what my code ends up computing.
• Residual: what my code actually computes — the signed vertical gap between the line's prediction and the true value at one tower's x-position: m·tx[i] + b − ty[i]. This is exactly what fxy returns in Stage 1. It can be positive (line too high) or negative (line too low).
• Squared error: the residual, squared: (m·tx[i] + b − ty[i])². Always non-negative — one tower's individual contribution to "badness."
• Total cost: the sum of squared errors over every tower — one number describing how bad the whole line is, for the whole dataset, for one specific (m, b).

The cost function, exactly as I wrote it

# cost function

def cost(tx,ty,m,b):
    c = 0
    for i in range(0,np.size(tx)):
        c += (m*tx[i] + b - ty[i])**2

    return c

Notice the loop: this sums one squared error per tower into a running total c. That loop is going to matter a lot later — Stage 10 is entirely about replacing it with a vectorized NumPy expression.

Why the cost becomes a function of m and b, not of x and y

Once tx and ty are fixed (they're the towers — they don't change), the only things left free to vary are m and b. So C is a function of the parameters, evaluated using the fixed data. This is the same shift in perspective from Stage 1: I stop asking "what's the residual at this x" and start asking "what's the total cost for this choice of line."

Why square the error instead of using distance directly?

This genuinely confused me for a while, and it's worth spelling out because the notebooks never resolve it explicitly — it's implicit in every cost function I wrote, always squared, never square-rooted. A few real reasons squaring wins over using the raw (or absolute) residual or a full square-root distance:

• Sign cancellation: residuals can be positive or negative. Summing raw residuals lets errors on opposite sides of the line cancel out, making a badly-fit line look artificially good. Squaring removes the sign, so every error contributes positively.
• Differentiability: |residual| (absolute value) also removes the sign, but it has a sharp corner at residual = 0, where the derivative doesn't exist. Squaring produces a smooth parabola with a well-defined derivative everywhere — essential once Stage 4 needs to differentiate the cost.
• No square root needed: a full distance formula would require √(...), which is more expensive to differentiate and doesn't change where the minimum sits (minimizing squared distance and minimizing distance itself give the same optimal point, since square-root is monotonic for non-negative inputs) — so squaring is "free" simplification, not an approximation.

Watching cost change as m, b change

Go back to the hero widget and move the sliders — the "Total squared error (cost)" readout is C(m,b) evaluated live. Notice it never goes negative, and it has one clear region where it's smallest. That region is the subject of Stage 3.

→ Next: cost is a function of two numbers, m and b. What does that function actually look like?

C(m,b) isn't a curve — it's a surface, and "minimum" means something specific on it

What I knew before this stage: cost is one number per (m,b) pair. What I didn't know: what the shape of "all possible costs" looks like.

From one number to a whole surface

Stage 2 gave me a way to compute C(m,b) for one specific line. But nothing stops me from computing it for every possible (m,b) — and if I do, and plot the result with m along one axis, b along another, and cost as height, I get a 3D surface. This is the object Newton's method's opening markdown cell in optimization_using_newton_method.ipynb is actually about:

What "minimum" means on this surface

Because every term in C(m,b) = Σ(m·xᵢ + b − yᵢ)² is a square, and a sum of squares (each convex in m,b) is itself convex, the surface for enough well-spread towers is bowl-shaped: one lowest point, no other dips to get trapped in. "Minimum" means the (m,b) pair sitting at the bottom of that bowl — the line that makes total squared error as small as it can possibly be.

What it means for the gradient to be zero

At the very bottom of a bowl, the surface is momentarily flat in every direction — nudging m slightly, or b slightly, doesn't change the height to first order. That "flatness in every direction" is exactly what "gradient equals zero" captures, and it's why Stage 4's derivatives matter: the minimum is defined as the point where both partial derivatives vanish simultaneously.

A hidden case where the bowl isn't a bowl

Look closely at cell 9 of regression_one_tower_problem.ipynb: I tested gradient descent using only one tower, the point (2, 3):

x,y = symbols('x y', positive='Integer')
x = 2
y = 3
cost_exp = (m*x + b - y)**2

With only one point, C(m,b) = (2m + b − 3)² is zero along the entire line 2m + b = 3 — not at one point. That's a valley, not a bowl: infinitely many (m,b) pairs achieve the exact minimum (cost = 0), because one equation (one tower) can't pin down two unknowns (m and b). Sure enough, the notebook's gradient descent converged to (m, b) ≈ (0.4, 2.2) — check it: 2(0.4) + 2.2 = 3.0, right on that line, but not the same point I'd get starting from a different initial guess. The bowl only becomes a true single-point bowl once enough towers (with different x-values) are included to fully determine both m and b — which is exactly what the 5-tower dataset does.

→ Next: to find the bottom of the bowl mathematically, I need the tool that finds where a surface is flat — differentiation.

Differentiation tells me which way is downhill

What I knew before this stage: the minimum is where the surface is flat. What I didn't know: how to find that point without checking every (m,b) by hand.

Why differentiation is needed

Checking every possible (m,b) pair to find the smallest cost isn't feasible — there are infinitely many. Differentiation gives something far more useful than a brute-force scan: at any single (m,b), it tells me the direction in which cost is increasing fastest. If I know which way is "uphill," I automatically know which way is "downhill," and I can move that way to make cost smaller — without ever having evaluated the whole surface.

The derivative with respect to m

Starting from C(m,b) = Σᵢ(m·xᵢ + b − yᵢ)², treat b as fixed and differentiate with respect to m using the chain rule — each term is (something)², so its derivative is 2·(something)·(derivative of something):

The xᵢ at the end comes from ∂/∂m of (m·xᵢ), which is just xᵢ. So each tower contributes its own error, scaled by its own x-position.

The derivative with respect to b

Same chain rule, but ∂/∂b of (b) is just 1, so no extra factor appears — every tower contributes its error unscaled.

The gradient, as a pair of partials

The gradient is nothing more mysterious than "both partial derivatives, collected into one vector." Geometrically, it points in the direction of steepest ascent on the cost surface — the direction that increases C the fastest, from wherever you're currently standing. Moving in the opposite direction is therefore the direction of steepest descent — this single geometric fact is the entire justification for gradient descent in Stage 5.

My first implementation — where I actually made a mistake

# now calculate the gradient that use for gradient descent

def gradient(m,b):
    dm = np.zeros(len(tx))
    db = np.zeros(len(tx))

    for i in range(len(tx)):
        error = m*tx[i] + b - ty[i]
        dm[i] = 2*error*tx[i]
        db[i] = 2*error

    return dm,db


dcdm , dcdb = gradient(2,3)
print(dcdm,dcdb)

There's also a second, quieter issue in this same notebook: the line right after it, defining gradient_descent(dcdm, dcdb, m, b, alpha, iteration), was written to take dcdm, dcdb as arguments directly — but it also contains print("m: "+ m + " b: " + b), string-concatenating numbers, which throws a TypeError the moment it actually runs. There is no output cell underneath it in the notebook — it was written but never successfully executed. It gets quietly superseded a cell later by a version that works.

The version that actually works — a single point, summed correctly

def dfdm(a,c):
    return diff(cost_exp,m).evalf(subs={m:a,b:c})

def dfdb(a,c):
    return diff(cost_exp,b).evalf(subs={m:a,b:c})

Here I switched strategy entirely: instead of manually looping and manually differentiating by hand, I let SymPy's diff() compute the exact symbolic derivative of cost_exp, then substitute in numeric values with .evalf(subs=...). This sidesteps the "forgot to sum" mistake completely, because SymPy differentiates the already-summed cost expression — there's no per-point array to accidentally leave unsummed. This is the moment symbolic differentiation enters the journey (formalized in Stage 7).

→ Next: with a working gradient, I can finally build an algorithm that walks downhill step by step.

Walking downhill, one step at a time — and watching it go wrong in three different ways

What I knew before this stage: the negative gradient points downhill. What I didn't know: how big a step to take, or what happens when I get that wrong.

The mathematical idea

If the negative gradient points downhill from wherever I'm standing, then repeatedly taking a small step in that direction should walk me toward the minimum. That's the entire idea — everything else is bookkeeping around how big the step is and how many times I take it.

The update rule

• α (alpha), the learning rate: how big a step to take. Too small and progress crawls; too large and the step overshoots the minimum entirely.
• Initial parameters: where the walk starts. Since (for a true bowl) the surface is convex, the starting point only affects how long the walk takes, not where it ends up — unlike the single-tower valley case from Stage 3, where the starting point determines which minimum is found.
• Iterations: how many steps to take before stopping.

Why does this move toward the minimum at all? Because at every step, moving a small enough distance in the negative-gradient direction is guaranteed to decrease the cost (that's what "steepest descent direction" means, to first order). Repeating that enough times, with small enough steps, has to end up at a point where there's no more downhill left to go — the minimum.

Mistake — gradient ascent instead of descent

def gradient_descent(dcdm,dcdb,m,b,alpha=0.2,iteration=100):
    for i in range(iteration):
        m , b = m + alpha*dcdm , b + alpha*dcdb
        print("m: "+ m , " b: " + b)

    return m,b

The corrected direction shows up immediately in the very next cell, grad_desc(dfdm,dfdb,a,c,alpha,iteration), which subtracts: a, c = a - alpha*dfdm(a,c) , c - alpha*dfdb(a,c) — this is the version that actually converges.

Experiment — overshoot as a stable oscillation

Cell 10 of the same notebook tries several learning rates on the single-tower problem starting from (m,b) = (18, 23). One of the runs produces this, forever:

This is overshoot in its purest form: the step size is large enough that each update leaps clean over the minimum and lands on the opposite side, symmetric enough that the next update leaps right back. Instead of converging, the walk bounces between two points forever — neither approaching the minimum nor diverging to infinity, just perfectly overshooting every single time. Mathematically, this happens because for this single-point cost (a valley, from Stage 3), the update along the direction perpendicular to the valley behaves like x ← x − α·k·x for some constant k; when α·k = 2, that update becomes x ← −x — an exact flip every step, which is precisely a 2-cycle.

Experiment — real divergence

In linear_regression.ipynb, the same 5-tower dataset with a learning rate of 0.5 blows up completely within 10 iterations:

Here α·k is greater than 2 (not just equal to it), so each step doesn't just flip sign — it flips sign and grows in magnitude. Every update makes the next update larger, and the sequence runs away to infinity exponentially fast. This is the mathematical reason "too large a learning rate" doesn't just slow convergence — past a threshold, it guarantees divergence.

Experiment — painfully slow convergence

The opposite failure shows up a few cells later: 1000 randomly-scaled towers (x up to 500), learning rate shrunk all the way to 1e-8, run for up to 1,000,000 iterations — and it still hadn't fully converged (b was creeping toward 7 at a rate of about +0.00004 per 100 iterations by the end). A learning rate small enough to avoid divergence on the large-scale feature (m) turns out to be far too small to make timely progress on the small-scale one (b). This exact experiment is what Stage 13 comes back to explain properly, once feature scale is named as the real culprit.

The α slider is scaled so tiny and huge learning rates are both reachable. Try α around 0.03–0.06 (converges smoothly), around 0.09–0.1 (oscillates, similar to the two-value repeat above), and above ~0.13 (diverges, similar to the 1e23 blow-up above).

→ Next: watching the numbers scroll by isn't very satisfying. What does this optimization actually look like?

Turning a scrolling wall of numbers into something I could actually watch happen

What I knew before this stage: gradient descent produces a sequence of (m,b) values. What I didn't know: how to see that sequence as motion.

Storing parameter history

To animate anything, I needed every intermediate (m,b), not just the final answer — so gradient_descent in regression_multi_tower_problem.ipynb started appending each step to a list:

def gradient_descent(dcdm,dcdb,a,c,alpha=0.1,iteration=100):
    history = []
    cost_data = []
    for i in range(iteration):
        a , c = a - alpha*dcdm(a,c) , c - alpha*dcdb(a,c)
        cos = cost.evalf(subs={m:a,b:c})

        history.append((a,c))
        cost_data.append(cos)
        print("m:",a,"b:",c,"cost:",cos)

    return a,c,history,cost_data

This version is safe without any special care: a and c are plain Python floats, and a, c = a - ... rebinds them to brand-new float objects every iteration. Appending the tuple (a,c) stores a snapshot of whatever those floats were at that moment — there's nothing left to mutate later, because floats and tuples are immutable.

Why .copy() became necessary later

By the time I wrote the LinearRegression class in linear_regression.ipynb, the parameters had become a single mutable container, m_b (a NumPy array), updated by assigning into its slots rather than rebinding the whole thing:

def gradient_descent(self,m_b,alpha,iteration):
    history = []
    grad = self.gradient()

    for i in range(iteration):
        dm = grad[0].evalf(subs={self.m:m_b[0],self.b:m_b[1]})
        db = grad[1].evalf(subs={self.m:m_b[0],self.b:m_b[1]})
        m_b[0],m_b[1] = m_b[0] - alpha*dm , m_b[1] - alpha*db
        history.append(m_b.copy())
        print(f"m: {m_b[0]} b:{m_b[1]}")

    return history

There's real, visible evidence in the notebook of exactly this kind of aliasing having consequences, even though the failure mode above isn't literally shown mid-crash: in cell 4, a_c is passed into l2.gradient_descent(a_c, ...), which mutates it in place and returns it, converged, around (3.02, 0.0097). Two cells later, cell 6 calls l2.gradient_optimized(a_c, 0.00000001, 1000000, 100) without redefining a_c first — and its very first printed line is already m: 3.000147... b: 6.950..., picking up mid-flight from wherever the array was left, not starting fresh from (0,0). That's the same "the array remembers what happened to it" behavior that makes .copy() necessary for history — just showing up between cells instead of inside a single loop.

The animation itself

from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(figsize=(6, 6))
ax.set_xlim(0, 10); ax.set_ylim(0, 10)
ax.scatter(tx, ty)
line, = ax.plot([], [], lw=2)
x = [0, 10]

def init():
    line.set_data([], [])
    return line,

def animate(frame):
    m = history[frame][0]
    b = history[frame][1]
    y = [m * i + b for i in x]
    line.set_data(x, y)
    return line,

anim = FuncAnimation(fig, animate, init_func=init, frames=len(history), interval=50, blit=True)
anim.save("gradientdescentline.mp4", writer="ffmpeg", fps=30)

Conceptually, FuncAnimation just calls animate(frame) repeatedly, once per frame, handing it an increasing frame number. animate looks up history[frame] — the (m,b) recorded at that exact step of gradient descent — recomputes the two endpoints of the line at x=0 and x=10, and redraws it. So each frame is literally one iteration of gradient descent: frame 0 is the starting line, frame 1 is after one update, and so on, with the towers plotted once as a fixed scatter underneath. Before committing to the full animation, cell 7 sanity-checked this by hand, manually plotting the line at a handful of chosen iterations (0, 5, 50, 100, 1000, 99999) to see it visibly walking toward the towers before generating the full video.

The cost-vs-iteration plot

iteration = [i for i in range(10000)]
plt.xlim(0,500); plt.ylim(0,2000)
plt.xlabel('iteration number'); plt.ylabel('cost')
plt.plot(iteration,cost_data)
plt.scatter(iteration,cost_data)
plt.show()

This is a different, complementary view from the line animation: instead of watching the line move in data-space, it watches the single number cost_data[i] — the height on the Stage 3 surface — fall over iterations. The two plots are two projections of the same underlying walk down the bowl.

→ Next: all of this leaned on SymPy's symbolic .evalf() inside the loop. Why did that eventually need to change?

From "let SymPy differentiate it" to "I don't need SymPy at all anymore"

What I knew before this stage: SymPy's diff() gives an exact derivative. What I didn't know: that evaluating it repeatedly, symbolically, is slow — and eventually unnecessary.

Step 1 — pure symbolic, evaluated with .evalf() every iteration

The earliest working gradient descent (Stage 4/5) called .evalf(subs=...) inside the loop, every single iteration:

dm = diff(cost_exp,m).evalf(subs={m:a,b:c})

This works, but diff(cost_exp,m) is a symbolic SymPy expression, and .evalf(subs=...) re-walks that whole expression tree and substitutes numbers into it from scratch, every single call. For a handful of iterations on 5 towers this is invisible; for hundreds of thousands of iterations on 1000 towers, it's the main reason the slow experiments in Stage 5 and Stage 13 take so long to run.

Why a NumPy array can't just use .evalf()

SymPy's .evalf() is a method that exists on SymPy expression and number objects — it doesn't know anything about NumPy arrays, and NumPy arrays don't know how to call it element-by-element automatically. This is directly visible in optimization_using_newton_method.ipynb:

def gradient_c(a_c):
    return np.array([dm.evalf(subs={m:a_c[0],b:a_c[1]}),db.evalf(subs={m:a_c[0],b:a_c[1]})])

Step 2 — lambdify: compiling symbolic expressions into numeric functions

regression_multi_tower_problem.ipynb introduces the fix:

from sympy.utilities.lambdify import lambdify

dcdm_exp = diff(cost,m)
dcdb_exp = diff(cost,b)

dcdm = lambdify((m,b),dcdm_exp,'numpy')
dcdb = lambdify((m,b),dcdb_exp,'numpy')

diff() is still doing the actual calculus — nothing about the math changes. What changes is that lambdify takes the resulting symbolic expression once and compiles it into an ordinary Python function that uses plain NumPy/Python arithmetic internally. From that point on, dcdm(a,c) is just a normal function call — no expression tree walking, no substitution dictionary, no SymPy Float objects. The symbolic work happens exactly once, up front; every iteration afterward is pure numeric evaluation.

The LinearRegression class in linear_regression.ipynb makes this speed difference explicit by keeping both versions side by side: gradient_descent (still using .evalf(subs=...) every step) and gradient_optimized (using lambdify once, then calling the compiled function every step) — literally the "before" and "after" of this exact upgrade, kept in the same file.

One inconsistency worth preserving: even after discovering lambdify and using it in gradient_optimized, newton_optimization in that same class still calls .evalf(subs=...) inside its loop. The lesson from Step 1 hadn't yet been applied everywhere — a very normal, uneven kind of refactoring where a better pattern is adopted in the place it was just learned, but not yet propagated to every other function that could use it.

Step 3 — dropping symbolic computation entirely

The final version, MultiRegressionGeneral in multiple_regression.ipynb, doesn't import SymPy at all. Once the gradient of a sum-of-squares cost has a known closed form (Stage 10's 2·Xᵀe), there's no need to symbolically differentiate anything ever again — the formula can just be written directly in NumPy:

def gradient(self):
    prediction = (self.param[:-1, None] * self.data).sum(axis=0) + self.param[-1]
    error = prediction - self.output
    grad_w = 2 * (self.data @ error)
    grad_b = 2 * np.sum(error)
    return np.append(grad_w,grad_b)

This is the real end-point of Stage 7's arc, and it's a better sequence than "symbolic → lambdify → done": symbolic differentiation (SymPy) was a tool for discovering and checking the formula, not something a finished implementation needs to carry forward. Once the formula for the gradient of squared error is known and trusted, computing it is pure linear algebra — faster than even the lambdified version, because there's no compiled-expression call overhead at all, just array operations.

Symbolic vs. numeric vs. vectorized numeric, side by side

• Symbolic representation: diff(cost,m) — an expression tree, exact, general, but not evaluated.
• Numerical evaluation: .evalf(subs={...}) or a lambdify'd function called on scalars — one specific number, computed once.
• Vectorized numerical computation: 2 * (self.data @ error) — every tower's contribution computed simultaneously as array operations, with no Python-level loop or per-call overhead at all (fully explored in Stage 10).

→ Next: gradient descent only uses first-order (slope) information. What if I use curvature too?

Using curvature, not just slope — and reaching the minimum in one jump

What I knew before this stage: gradient descent needs many small steps. What I didn't know: a method that uses how the slope itself is changing, not just the slope.

Why gradient descent only uses first-order information

The gradient tells you the slope at exactly one point — it says nothing about whether that slope is about to get steeper or flatter as you move. Gradient descent has to find that out the hard way, by taking a step and re-measuring the slope again. Curvature — how fast the slope itself changes — is second-order information, and if it's known, it can be used to jump much more directly toward the minimum instead of repeatedly re-measuring and re-stepping.

The Hessian: curvature as a matrix

For a function of two parameters, there isn't just one curvature — there's how ∂C/∂m changes as m changes, how it changes as b changes, and the same again for ∂C/∂b. That's four numbers (two of which turn out equal), naturally arranged as a matrix:

It's a matrix (not a vector, like the gradient) precisely because curvature is inherently about pairs of directions — how the slope in one direction changes as you move in another — and with n parameters there are n×n such pairings.

My implementation, from the Newton notebook

def gradient(cost):
    dm = diff(cost,m)
    db = diff(cost,b)
    return dm,db

def hessian(dm,db):
    dmb = diff(dm,b)
    dmm = diff(dm,m)
    dbm = diff(db,m)
    dbb = diff(db,b)
    return np.array([[dmm,dmb],[dbm,dbb]],dtype=float)

Note dtype=float here — this only works because dmm, dmb, dbm, dbb come out as constants with no remaining m or b symbols in them. That's not a coincidence; it's the whole reason Newton behaves so well on this problem (see below).

The Newton update

def newton_optimization(a_c,grad,hess,iteration):
    for i in range(iteration):
        a_c = a_c - np.matmul(np.linalg.inv(hess),grad(a_c))
        print("m:",a_c[0],"b:",a_c[1])
    return a_c

Every remaining iteration after the first prints the exact same (m,b) — Newton's method reached the minimum on its very first step and then just kept confirming it was already there.

Gradient, Hessian, inverse Hessian, and the Newton step — what each one is

• Gradient (∇C): a vector, direction of steepest ascent. Shape: (n,) for n parameters.
• Hessian (H): a matrix of second derivatives, describing curvature in every pair of directions. Shape: (n,n).
• Inverse Hessian (H⁻¹): "undoes" the curvature — it rescales the gradient direction by how flat or steep the surface actually is in each direction, instead of taking a fixed-size step regardless of curvature.
• Newton step (H⁻¹∇C): the actual displacement applied to θ. Geometrically, it's "the gradient direction, corrected for how the surface curves," rather than "the gradient direction, scaled by a fixed learning rate α."

Why Newton reaches the optimum in essentially one step here

This only happens because of a specific chain of facts about squared-error regression, worth spelling out in full:

1. The objective is quadratic in the parameters: C(m,b) = Σ(mxᵢ+b−yᵢ)² expands into terms with m², b², and mb — a quadratic form, nothing higher-order.
2. Its gradient is therefore linear in the parameters: differentiating a quadratic once always produces something linear (compare: d/dx(x²) = 2x, linear in x).
3. Its Hessian is therefore constant: differentiating something linear once more produces a constant — no m or b left in it at all, which is exactly why dtype=float worked directly on hessian(dm,db) above.
4. Newton can jump directly to the optimum: Newton's update is exact whenever the true function is exactly quadratic, because a quadratic function is fully described by its value, gradient, and (constant) curvature at any single point — Newton's step effectively solves "where would this exact quadratic hit zero gradient" directly, in one linear-algebra step, rather than approximating that answer iteratively.

This generalizes cleanly, as confirmed independently in multiple_regression.ipynb's MultipleRegression.newton_method, which reaches the known parameters [2, 3, 4] starting from [0, 0, 0] immediately, then repeats that answer for every further iteration — the exact same one-step behavior, now with three parameters instead of two. My own comment in that notebook records the insight directly: for a quadratic objective, Newton's method reaches the target in one iteration no matter how many parameters are involved, because the Hessian stays constant regardless of how many dimensions θ has.

Extension — a numerical improvement I didn't implement

→ Next: everything so far used one feature, x. What changes when a tower — or any data point — has more than one number describing it?

y = mx + b becomes y = w₁x₁ + w₂x₂ + ... + wₙxₙ + b

What I knew before this stage: one input, one output, two parameters. What I didn't know: how the same idea scales when a tower is described by more than one number.

Why one feature stopped being enough

Every tower so far had exactly one x-value. But in multiple_regression.ipynb, I explicitly moved to a setting with two inputs per sample — described in my own comment as "two dependent variable distance D and price/unit P" — and asked whether an output could be predicted from both of them at once, not just one. That requires more than one slope: one weight per input feature, plus a single shared bias.

Vocabulary, precisely

• Feature: one measured input quantity (e.g. distance, or price-per-unit). Each sample has n of them.
• Sample: one complete data point — a full set of feature values, paired with one target value. In the tower context, one tower's full description.
• Target / output (y): the value being predicted for a sample.
• Parameter / weight (wⱼ): how much feature j contributes to the prediction, per unit of that feature — the direct generalization of "m".
• Bias (b): the same role as before — a constant baseline added regardless of feature values.

Matrix dimensions — the part I have to get exactly right

This is the part of the notebooks where I actually used two different data orientations in two different classes, and I'm not silently picking one — both are real, both appear in my code, and knowing which is which matters.

w1,w2,b = symbols('w1 w2 b')
param = [w1,w2,b]
data = [[1,2,0],[1,3,0],[3,7,0],[6,4,0],[4,9,0],[10,2,0],[12,15,0]]
# data[i]      -> one full sample: [x1, x2, y]
# data[i][j]   -> feature j of sample i, for j in 0..1

distance = np.array([1,3,5,10,8,11,16,21,34,2,0], dtype=float)
price    = np.array([12,3,6,8,9,0,5,6,14,1,56], dtype=float)
data = [distance, price]
# data[0]      -> the ENTIRE distance feature, across all samples
# data[1]      -> the ENTIRE price feature, across all samples

Neither orientation is "wrong" — plenty of real ML code uses samples × features (it's the more common convention, e.g. scikit-learn's X), and features × samples is also legitimate and, as Stage 10 shows, made my particular vectorized gradient formula (data @ error) fall out very naturally. What matters is being explicit about which one a given piece of code assumes, since the two are not interchangeable without transposing.

Why the same gradient-descent algorithm can work for any number of features

Nothing about the update rule θ ← θ − α·∇C(θ) cares how many entries θ has. Whether θ = [m, b] (2 entries) or θ = [w₁,...,wₙ,b] (n+1 entries), it's still "subtract a scaled gradient from a parameter vector." Stage 10 and 11 show exactly how the loop-based version generalizes into this shape.

→ Next: with n features, per-point Python loops start to look expensive. What does the same math look like as matrix operations?

Every Python for-loop over samples is secretly a matrix operation waiting to happen

What I knew before this stage: for-loops over towers compute the right numbers. What I didn't know: that the same numbers come out of a couple of matrix operations, much faster.

Starting point — the loop-based cost

def cost(self):
    c = 0
    for i in range(len(self.data)):
        gather = 0
        for j in range(len(self.param)-1):
            gather += self.param[j]*self.data[i][j]
        c += (gather + self.param[-1] - self.data[i][-1])**2
    return c

Two nested loops: the outer one walks samples, the inner one walks features within a sample, manually accumulating gather = Σⱼ wⱼxᵢⱼ before adding the bias and squaring the error. Every single multiplication and addition here has a direct matrix-operation counterpart.

The same calculation as matrix operations

def cost(self):
    prediction = (self.param[:-1, None] * self.data).sum(axis=0) + self.param[-1]
    error = prediction - self.output
    return np.sum(error**2)

Reading this line by line, against the loop it replaces:

• self.param[:-1, None] — reshapes the weight vector from shape (n,) to (n,1) so it can be broadcast against every sample's column at once. This is broadcasting: NumPy stretches the (n,1) array across the sample axis without actually copying it n times.
• self.param[:-1, None] * self.data — an element-wise multiplication (not matrix multiplication): for a features×samples-shaped self.data, this multiplies every feature's row by its own weight, for every sample simultaneously. This replaces the inner loop's self.param[j]*self.data[i][j], done for every (i,j) pair at once instead of one at a time.
• .sum(axis=0) — sums down each column (over the features axis, axis 0, in a features×samples array), collapsing n feature-contributions into one prediction per sample. This is exactly the inner loop's gather accumulator, computed for every sample in one call instead of looped.
• + self.param[-1] — adds the bias to every sample's prediction at once (broadcasting a single scalar across the whole array).
• error**2 then np.sum(...) — squares every residual element-wise, then adds them all up. This replaces the outer loop's running total c.

The gradient — where @ (matrix multiplication / dot product) enters

def gradient(self):
    prediction = (self.param[:-1, None] * self.data).sum(axis=0) + self.param[-1]
    error = prediction - self.output
    grad_w = 2 * (self.data @ error)
    grad_b = 2 * np.sum(error)
    return np.append(grad_w,grad_b)

self.data @ error is a matrix-vector product. Because self.data is features×samples (shape (n, m)) and error is (m,), the product has shape (n,) — one number per feature. What does that number mean? Row j of self.data is feature j's value across every sample; taking its dot product with error computes exactly Σᵢ xⱼᵢ · errorᵢ — which is precisely ∂C/∂wⱼ (up to the factor of 2), for every feature j, simultaneously. This single @ replaces an entire loop over both samples and features.

Vocabulary, tied directly to the code above

• Element-wise multiplication (*): multiplies matching positions independently; shapes must match or broadcast.
• Broadcasting: NumPy's rule for stretching smaller-shaped arrays across a larger one without copying data, e.g. turning a (n,1) weight column into an effective (n,m) grid to multiply against (n,m) data.
• Matrix multiplication / dot product (@): sums products across a shared axis — turns "multiply then sum" (two separate loop-conceptual steps) into a single operation.
• np.sum(..., axis=...): collapses one specific axis by summation; axis=0 sums down columns, axis=1 sums across rows — getting this backwards silently produces a result with the wrong shape or the wrong meaning.
• @ vs. plain np.sum: @ is "multiply corresponding entries, then sum" fused into one call; when a computation is exactly of that shape, @ replaces both a multiplication and a subsequent sum.

Why vectorization is faster

The loop version pays Python's per-iteration overhead (bytecode dispatch, bounds checks, object creation for each intermediate float) once for every single sample-feature pair. The vectorized version dispatches into NumPy's compiled C loops once for the whole array, regardless of how many samples or features there are — the actual arithmetic is comparable, but the overhead around each individual operation collapses from "once per number" to "once total." This is a big part of why the slow 1,000,000-iteration experiments in Stage 5 and Stage 13 were even feasible to run at all once the code moved to this form.

How the feature count generalizes for free

Nothing above hard-codes "2 features." self.param[:-1, None] * self.data works whether self.data has 2 rows or 20 — the shapes just have to agree. This is the concrete mechanism behind Stage 9's claim that the same algorithm handles any number of features: the loop bounds that used to be explicit (range(len(self.param)-1)) are replaced by array shapes that NumPy checks automatically.

→ Next: with vectorized operations and n features in hand, what does the fully general implementation look like?

From hard-coded w1, w2, b to a parameter vector of any length

What I knew before this stage: two named weights plus a bias, written by hand. What I didn't know: how to stop hard-coding the parameter count entirely.

Where I started — named parameters

w1,w2,b = symbols('w1 w2 b')
param = [w1,w2,b]

MultipleRegression works, but every one of its three parameters exists because I typed its name. Adding a third feature would mean editing the class by hand — adding a w3 symbol, extending the inner loop's range, and so on. That's fine for exploring the idea with 2 features, but it isn't a general implementation yet — it's a specific one that happens to have 2 features baked in structurally (even though the loop itself, for j in range(len(self.param)-1), was already written generically).

Where I ended up — an arbitrary-length parameter vector

class MultiRegressionGeneral:
    def __init__(self,data,output,param):
        self.data = np.asarray(data, dtype=float)
        self.output = np.asarray(output, dtype=float)
        self.param = np.asarray(param, dtype=float)

param is no longer a hand-written list of named SymPy symbols — it's a plain NumPy array, and nothing in the class cares how long it is. The transition looks like this:

This is visible directly in how differently the class gets used across two experiments in the same notebook. First, with 2 features:

data = [distance, price]
param = np.array([2,3,4],dtype=float)
m = MultiRegressionGeneral(data,output,param)

Then, without a single change to the class itself, with 5 features:

data = [distance, price, area, age, rooms]
output = 2*distance + 3*price + 0.5*area - 4*age + 10*rooms + 20
param = np.zeros(6)   # 5 weights + 1 bias — the "6" is the only thing that changed
m2 = MultiRegressionGeneral(data,output,param)

The only number that changed between a 2-feature and a 5-feature regression is the length of param (and, correspondingly, the number of rows in data). No method inside the class needed editing.

Why the same gradient-descent algorithm still works

Trace it through: self.param[:-1, None] * self.data broadcasts across however many feature-rows exist; .sum(axis=0) collapses however many there are; self.data @ error produces a gradient vector exactly as long as the number of feature-rows in self.data. Every shape in the computation is derived from the shapes of self.data and self.param, never hard-coded — which is precisely what makes "generalized" true in a structural sense, not just in the sense of "happens to work for the numbers I tried."

→ Next: a general implementation is only useful if I can trust it. How did I actually verify it was correct?

If I already know the answer, my optimizer had better find it

What I knew before this stage: the code runs without crashing. What I didn't know: whether it's actually computing the right thing.

The debugging trick — generate data from a known equation

A regression implementation can run cleanly, produce plausible-looking numbers, and still be subtly wrong — the one_tower gradient bug in Stage 4 proved that. The fix I used throughout multiple_regression.ipynb: instead of using real, noisy, unknown data, generate data from a known formula, and check whether the optimizer can find that exact formula back out.

def cost(x1,x2):
    return 2*x1 + 3*x2 + 4

data = [[1,2,0],[1,3,0],[3,7,0],[6,4,0],[4,9,0],[10,2,0],[12,15,0]]
for i in range(len(data)):
    data[i][2] = cost(data[i][0],data[i][1])   # fill in y using the KNOWN true formula

Every target value here comes from a formula I chose myself: y = 2x₁ + 3x₂ + 4. So the "correct" parameters are known in advance to be exactly [2, 3, 4] — I'm not hoping the optimizer finds a good fit, I'm checking whether it finds the fit I built the data to have.

Known parameters → generate y → initialize wrong → optimize → compare

opt = np.array([0,0,0],dtype=float)   # deliberately wrong: nowhere near [2,3,4]
m.gradient_descent(opt,0.001,10000)

Starting from [0,0,0], gradient descent visibly drifts toward the true [2,3,4] — the sequence of printed vectors is the tell: if the implementation were wrong, it would converge to some other triplet, or not converge at all, or (as in Stage 5's oscillation and divergence experiments) never settle down. Newton's method, tested the same way, jumps straight to [2. 3. 4.] on the very first iteration — independent confirmation from a completely different method, which matters, since if both a slow iterative method and an exact one-step method agree on the same answer, that's much stronger evidence than either alone.

This exact same test was repeated at larger scale with MultiRegressionGeneral: a known formula cost(x1,x2) = 2·x1 + 3·x2 + 4, applied to 11 synthetic distance/price samples, recovered as array([2.00000007, 3.00000004, 3.99999839]) after 200,000 iterations — accurate to about 6 decimal places, which is exactly what I'd expect from a converging-but-not-yet-perfectly-converged iterative method.

The stronger verification: correct parameters ⇒ cost = 0 and gradient = 0

Recovering approximately the right answer after optimizing is good evidence, but there's a sharper, more direct test that doesn't even require running the optimizer: if the parameters are exactly right, the cost function should evaluate to exactly zero (every prediction matches every target exactly), and the gradient should evaluate to exactly the zero vector (there's no downhill direction left, because you're already at the minimum). I checked both, directly:

print("Current:")
print(m2.param)
print("Cost:", m2.cost())
print("Gradient:", m2.gradient())

true_param2 = np.array([2., 3., 0.5, -4., 10., 20.])
m2.param = true_param2

print("\nTrue parameters:")
print(m2.param)
print("Cost:", m2.cost())
print("Gradient:", m2.gradient())

This is exactly the confirmation I wanted: plugging in parameters that hadn't yet converged gives a non-zero cost and a non-zero gradient — there's still room to improve, and the gradient correctly points toward that room. Plugging in the exact known true parameters gives cost exactly 0.0 and gradient exactly the zero vector, with no rounding residue at all. That's strong evidence the cost and gradient formulas are mathematically consistent with each other — the gradient really is the derivative of the cost I implemented, not a plausible-looking but subtly disconnected formula.

Why this kind of test matters more than it looks like it does

It's tempting to trust code because it runs and produces numbers in a reasonable range — the Stage 4 gradient bug did exactly that. A synthetic test with a known answer converts "does this look plausible?" into "does this match a number I can check exactly?" — turning debugging from inspection into verification.

→ Next: the implementation is verified correct. So why did it get so slow on realistic, differently-scaled features?

One learning rate, two very different curvatures

What I knew before this stage: gradient descent can oscillate, diverge, or crawl. What I didn't know: why the exact same algorithm could be both too aggressive and too timid on the exact same run.

Coming back to the experiment that didn't add up

Stage 5 ended with an experiment I described but didn't explain: 1000 towers, x-values scaled up to 500, a learning rate shrunk all the way to 1e-8 just to keep things stable, run for 1,000,000 iterations — and b was still creeping toward its true value at roughly +0.00004 per 100 iterations at the end. Nothing about that run was broken. Every single update was a legitimate step downhill. It was just almost unusably slow. That's a different failure mode from oscillation or divergence, and it deserves its own explanation.

Revisiting the per-parameter curvature from Stage 8

Stage 8 already derived the exact Hessian for this cost function, without naming why it mattered yet:

The diagonal entries are the closest thing to "curvature along m alone" and "curvature along b alone": roughly k_m ≈ 2Σxᵢ² and k_b ≈ 2n (this is an approximation — the true axes of the bowl are rotated whenever Σxᵢ ≠ 0, but it's close enough to see the point). Stage 5's own derivation already showed that gradient descent's stability condition is |1 − α·k| < 1, i.e. roughly 0 < α < 2/k, for whichever curvature k is in play.

Why one α can't serve both parameters at once

Plug in the 1000-tower experiment's shape: with x-values up to 500, Σxᵢ² is enormous — on the order of tens of millions — so k_m is huge, and the stable range for updating m alone is a razor-thin 0 < α < 2/k_m, forcing α down into the 1e-8 territory just to avoid the divergence pattern from Stage 5. But k_b = 2n = 2000 is comparatively tiny — the stable range for updating b alone would happily allow an α many orders of magnitude larger. Gradient descent only has one α to share between both updates. Set small enough to keep m stable, that same α makes every single b-update absurdly small too — even though b itself would have been perfectly happy converging fast. The whole optimization is throttled to the pace of its worst-conditioned direction.

Why Newton's method never had this problem

This is the payoff of Stage 8's H⁻¹∇C update, stated plainly: multiplying by the inverse Hessian rescales each direction by its own curvature before stepping. A huge k_m and a tiny k_b get corrected for automatically and separately — Newton's method doesn't share one global step size across mismatched directions the way gradient descent's scalar α does. That's precisely why Newton converged in a single step regardless of how the 5-tower or 1000-tower data happened to be scaled, while gradient descent's single α had to compromise between them.

The standard fix I didn't implement: feature scaling

The other way out: not iterating at all

There's a second, completely different escape from this whole problem, and I'd already built half of it without naming it: the hero widget's "Snap to best fit" button doesn't run gradient descent — it solves for the minimizing (m,b) directly, algebraically, using the same closed-form normal equations that make ordinary least-squares a solved problem: θ = (XᵀX)⁻¹Xᵀy. This is the same underlying idea as Newton's method — using the known, constant Hessian to jump straight to the answer — just phrased as "solve a linear system" instead of "take one optimizer step." It sidesteps learning rates and condition numbers entirely, at the cost of an explicit matrix inversion that gets expensive once the number of features grows large — which is exactly why iterative methods like gradient descent remain the practical choice once n is large, even though they're the ones sensitive to conditioning.

→ Next: this whole shape — model, cost, gradient, optimizer — didn't stop existing once the towers were solved. Where does it show up next?

The same four-part skeleton, wearing different clothes

What I knew before this stage: how to fit a line, and then a hyperplane, to towers. What I didn't know: that I'd just built the actual skeleton underneath most of modern machine learning.

Naming the skeleton explicitly

Strip away everything specific to towers and every stage of this journey reduces to the same four-part shape, repeated with different specifics each time:

1. A model — a hypothesis with adjustable parameters, mapping inputs to a prediction. Mine was ŷ = w·x + b.
2. A cost — a single number, computed from the parameters and the data, that gets smaller as predictions get better. Mine was Σ(ŷᵢ − yᵢ)².
3. A gradient — the cost's sensitivity to each parameter, telling me which way is downhill.
4. An optimizer — a rule for using the gradient (and sometimes curvature) to actually update the parameters, repeatedly, until the cost stops improving.

Nothing about this shape is specific to lines, towers, or even regression. It's the shape of "learning parameters from data" in general.

Swap the model → logistic regression

Wrap the same linear combination w·x + b in a sigmoid function to squash it into (0,1), and the model becomes a probability instead of a raw number — the beginning of classification instead of regression. The cost usually changes too (cross-entropy instead of squared error, because squared error stops being the "right" measure of badness once predictions are probabilities), but the gradient is still found by differentiating that cost with respect to the parameters, and the optimizer is still gradient descent (or Newton, or a variant) updating parameters by a scaled gradient. Nothing structurally new — a different model, a matching cost, the same machinery.

Swap the model → a neural network

Stack several of these linear-combination-plus-nonlinearity units in layers, feeding one layer's output into the next layer's input, and the "model" becomes a neural network. The cost is still a single number measuring prediction badness. The gradient is still needed for every parameter — except now there are potentially millions of them, spread across many layers, and computing ∂C/∂(each weight) by hand the way Stage 4 did for m and b would be hopeless. Backpropagation is the name for doing exactly that differentiation efficiently, by repeatedly applying the chain rule backward through the layers — which is the same chain rule Stage 4 used to differentiate a sum of squares, just applied through a much deeper composition of functions instead of one.

Swap the optimizer → momentum and Adam

Stage 5 and Stage 13 both surfaced the same weakness in plain gradient descent: a single global learning rate is either unstable in high-curvature directions or too slow in low-curvature ones. Rather than fixing this by rescaling features, a whole family of optimizers fixes it by making the update rule itself smarter — momentum accumulates a running average of past gradients so the walk builds speed through consistently-downhill directions and damps oscillation in noisy ones; Adam additionally tracks a per-parameter estimate of gradient magnitude and effectively gives each parameter its own adaptive learning rate — a softer, cheaper cousin of what Newton's method does exactly with the full Hessian. Both are still, at their core, "use the gradient to decide how to update the parameters" — just with a more sophisticated rule than "subtract α times the gradient."

Add a term to the cost → regularization

Nothing stops the cost function from being more than just squared error. Adding a penalty term like λΣwⱼ² (ridge regression) or λΣ|wⱼ| (lasso) discourages weights from growing arbitrarily large, which helps prevent a model from fitting noise in the training data too closely. The mechanical impact on everything built in this journey is small and direct: the new cost still needs a gradient, and since differentiation distributes over addition, the penalty term just contributes an extra 2λwⱼ (or similar) added onto the gradient Stage 4 already derived — the optimizer loop itself doesn't change at all.

Swap "all the data every step" → stochastic / mini-batch gradient descent

Every gradient in this journey was computed using the entire dataset at once — Stage 10's vectorized data @ error sums over every sample simultaneously. That's fine for 5 or even 1000 towers, but for datasets with millions of samples, computing one exact gradient before taking a single step becomes itself the bottleneck. Stochastic and mini-batch gradient descent compute an approximate gradient from a small random subset of the data each step instead — noisier, but far cheaper per step, and correct on average. This is a direct extension of exactly the vectorized sum from Stage 10, just applied to a subset of rows instead of all of them.

What this journey actually was

None of these are new ideas layered on top of what came before — they're the same four-part skeleton (model, cost, gradient, optimizer) with one part swapped at a time. Learning that skeleton on the smallest possible example — five towers and a straight line — is what makes it recognizable everywhere else it shows up.

→ Last: if I strip away every notebook, every bug, and every stage number, what's the one mental model I actually want to keep?

The whole journey, compressed back down to one sentence per stage

The central question from the top of this page, answered end to end.

The chain, compressed

1. A vague real-world problem only becomes solvable once it's translated into coordinates — and even then, my own problem statement and my own code quietly solved two different (both valid) problems.
2. The simplest relationship between coordinates worth proposing is a line, y = mx + b — two parameters, easy to reason about, almost never a perfect fit.
3. "Almost never a perfect fit" demands a number for badness — residual, then squared error, then total cost — precise enough to compare any two candidate lines.
4. A cost that depends on two parameters is really a surface, not a curve — and "best line" means the lowest point of that surface, which only has one answer once there's enough data to fully determine both parameters.
5. Finding the bottom of a surface without checking every point by hand requires differentiation — the gradient, pointing uphill, and its negative, pointing down.
6. Walking downhill repeatedly is gradient descent — and getting the step size wrong doesn't just slow it down, it can make it oscillate forever or blow up to infinity, depending on exactly how wrong.
7. Watching numbers scroll by isn't understanding — recording every step's parameters (carefully, respecting the difference between rebinding and mutating) turns a sequence of numbers into a video of a line visibly walking toward the data.
8. Differentiating symbolically every single iteration is correct but wasteful — compiling the derivative once (lambdify), and eventually writing the known closed-form gradient directly, removes that waste entirely.
9. Slope alone tells you which way, but curvature (the Hessian) tells you how far — and for a cost this well-behaved (exactly quadratic), curvature-aware Newton's method reaches the exact answer in a single step.
10. One input number stops being enough the moment towers (or anything) are described by more than one measurement — the model grows from y = mx + b to ŷ = w₁x₁ + ... + wₙxₙ + b, with no new ideas, just more parameters.
11. More parameters means more loops — until every one of those loops is recognized as a matrix operation in disguise: broadcasting, element-wise products, axis-sums, and dot products, all faster and all exactly equivalent to the loop they replace.
12. A class that "works for 2 features" isn't the same as a class that's actually general in the number of features — that only happens once every shape in the computation is derived from data, not hard-coded.
13. Code that runs isn't code that's correct — the sharpest test available is synthetic data with a known answer, checked two ways: does the optimizer converge to it, and does cost/gradient evaluate to exactly zero at the true parameters?
14. A working, verified optimizer can still be painfully slow — because one learning rate has to serve every parameter, and parameters with very different curvature (often from features on very different scales) can't all be served well by the same number.
15. None of this was really about towers — model, cost, gradient, optimizer is the same skeleton underneath logistic regression, neural networks, regularization, and modern optimizers; only the specific pieces change.

If I forget every notebook, this is what I want to still know

Back to the beginning: the hero widget at the top of this page is still just five towers and one line. Everything between here and there is the answer to "why did that ever need to become more complicated than dragging two sliders?"