Gradient descent step size and momentum: roll a ball down a narrow valley (Python)

Play: roll a ball down a narrow valley

A narrow valley runs diagonally across the map below, its walls 100 times more steeply curved than its floor, and a ball rolls down it by gradient descent. Each hop goes straight downhill, and its length is the local slope times a fixed number η, the step size: a steep spot gives a long hop, a flat one a short hop. The loss on screen is the height f of the ball's spot, 0 at the bottom. Tap the map to place the ball or keep it where it is, set the step size η and the momentum β with the sliders or their − and + buttons, and press Go. Step makes a single hop, Reset puts the ball back at its start, and Line fit swaps in a second landscape. Watch the trail and the step counter. At η = 0.019 the ball reaches the minimum after 385 steps. At η = 0.021, a step size only about 10 % larger, it flies off the map after 23.

Download the game

Show code

// game.js: roll a ball down a narrow valley by gradient descent with a step size η and a
// momentum β. The notebook writes both landscapes (Hessian, minimum, contour levels, the line
// fit's data) and a cividis color table into the #data element.
"use strict";

const DATA = JSON.parse(document.getElementById("data").textContent);
const TOL = 1e-3;                  // "at the minimum": closer than this
const MAX_STEPS = 20000;
const DOTS = 60;                   // the first steps get a dot each
const DECADES = 5;                 // the relief's colors span five decades of height
const RELIEF = 0.1;                // height of the shaded surface √(f − f*), for the hillshade
const LIGHT = [-1, 1, 1].map((c) => c / Math.sqrt(3));   // light from the upper left
const MAP_INK = "#1f2a44";        // ink and accent on the relief, which is the same in both themes
const MAP_ACCENT = "#c8553d";
const CIVIDIS = DATA.cividis.map((h) => [1, 3, 5].map((i) => parseInt(h.slice(i, i + 2), 16)));

// ---- the state: everything that changes, in one object
const state = {land: 0, eta: DATA.lands[0].eta.default, beta: 0, p: [0, 0], v: [0, 0],
  k: 0, status: "ready", start: [0, 0], trail: [], due: 0, hint: true};

function place(p) {                                     // the ball to p, at rest, trail cleared
  Object.assign(state, {p: [...p], v: [0, 0], k: 0, status: "ready", start: [...p], trail: [[...p]], due: 0});
}

function switchLand(i) {
  Object.assign(state, {land: i, eta: DATA.lands[i].eta.default});
  place(DATA.lands[i].start);
}

// ---- one gradient step with momentum; the notebook's descend() does the same arithmetic
function step(s) {
  const {H, pmin, half} = DATA.lands[s.land];
  const dx = s.p[0] - pmin[0], dy = s.p[1] - pmin[1];
  const g = [H[0][0] * dx + H[0][1] * dy, H[1][0] * dx + H[1][1] * dy];   // ∇f = H (p − p*)
  s.v = [s.beta * s.v[0] - s.eta * g[0], s.beta * s.v[1] - s.eta * g[1]];  // v ← βv − η∇f
  s.p = [s.p[0] + s.v[0], s.p[1] + s.v[1]];                                // p ← p + v
  s.k += 1;
  s.trail.push(s.p);
  const ex = s.p[0] - pmin[0], ey = s.p[1] - pmin[1];
  if (ex * ex + ey * ey < TOL * TOL) s.status = "arrived";
  else if (Math.max(Math.abs(ex), Math.abs(ey)) > half) s.status = "left";
  else if (s.k >= MAX_STEPS) s.status = "stopped";
}

const height = (lnd, dx, dy) => 0.5 * (lnd.H[0][0] * dx * dx + 2 * lnd.H[0][1] * dx * dy + lnd.H[1][1] * dy * dy);

// ---- the loop: slow at first, so that the zigzag can be followed, then faster
const rate = (k) => (k < 20 ? 8 : k < 400 ? 60 : 600);   // steps per second
function tick(dt) {
  if (state.status !== "running") return;
  state.due += Math.min(dt, 0.05) * rate(state.k);      // steps owed; the cap stops jumps
  while (state.due >= 1 && state.status === "running") {
    step(state);
    state.due -= 1;
  }
}

// Three significant digits: 4.9·10⁻⁷ when small, 0.00406 or 12.3 otherwise.
const SUP = {"-": "⁻", 0: "⁰", 1: "¹", 2: "²", 3: "³", 4: "⁴", 5: "⁵", 6: "⁶", 7: "⁷", 8: "⁸", 9: "⁹"};
function sig3(x) {
  if (x !== 0 && (Math.abs(x) >= 1e6 || Math.abs(x) < 1e-3)) {
    const [m, e] = x.toExponential(2).split("e");
    return `${Number(m)}·10${[...String(Number(e))].map((c) => SUP[c]).join("")}`;
  }
  return Number(x.toPrecision(3)).toLocaleString("en-US", {maximumFractionDigits: 6});
}

// ---- drawing
const canvas = document.getElementById("map");
const ctx = canvas.getContext("2d");
const $ = (id) => document.getElementById(id);
let view = {side: 100, dpr: 1};                        // the map's size in CSS px, set by layout()
let cache = {key: null, image: null};

function layout() {
  const W = document.documentElement.clientWidth, H = document.documentElement.clientHeight;
  const narrow = W < 600, dpr = window.devicePixelRatio || 1;
  let side = narrow ? W : Math.min(H, Math.round(0.6 * W));
  canvas.style.width = canvas.style.height = `${side}px`;
  if (narrow) {                                         // a short frame: a smaller map, so the panel fits
    const over = $("panel").scrollHeight - $("panel").clientHeight;
    if (over > 0) side = Math.max(160, side - over);
    canvas.style.width = canvas.style.height = `${side}px`;
  }
  canvas.width = canvas.height = Math.round(side * dpr);
  view = {side, dpr, narrow};
  render();
}

// Map units to CSS px: the square of half-width `half` around the minimum fills the canvas.
function toPixel(lnd, x, y) {
  const s = view.side / (2 * lnd.half);
  return [(x - lnd.pmin[0] + lnd.half) * s, (lnd.pmin[1] + lnd.half - y) * s];
}

// The relief, once per landscape and size: cividis of log height, lightened toward white,
// times a hillshade as if the surface √(f − f*) were lit from the upper left.
function relief(lnd, n) {
  const img = ctx.createImageData(n, n);
  const top = Math.max(height(lnd, lnd.half, lnd.half), height(lnd, lnd.half, -lnd.half));   // the highest corner
  const a = lnd.H[0][0], b = lnd.H[0][1], c = lnd.H[1][1];
  for (let j = 0; j < n; j++) for (let i = 0; i < n; i++) {
    const dx = ((2 * (i + 0.5)) / n - 1) * lnd.half, dy = (1 - (2 * (j + 0.5)) / n) * lnd.half;
    const h = height(lnd, dx, dy);
    const t = Math.min(1, Math.max(0, 1 + Math.log10(h / top) / DECADES));
    const r = RELIEF / (2 * Math.sqrt(h) || 1);          // slope of √h is ∇f / (2√h)
    const nx = -r * (a * dx + b * dy), ny = -r * (b * dx + c * dy);
    const shade = 1 + 0.3 * ((nx * LIGHT[0] + ny * LIGHT[1] + LIGHT[2]) / Math.hypot(nx, ny, 1) - LIGHT[2]);   // flat ground: 1
    const rgb = CIVIDIS[Math.round(63 * t)], o = 4 * (j * n + i);
    for (let q = 0; q < 3; q++) img.data[o + q] = (rgb[q] + 0.4 * (255 - rgb[q])) * shade;
    img.data[o + 3] = 255;
  }
  return img;
}

// The level sets of a quadratic are ellipses: ½ dᵀHd = c has semi-axes √(2c/λ) along H's main directions.
function drawContours(lnd, color) {
  const [a, b, c] = [lnd.H[0][0], lnd.H[0][1], lnd.H[1][1]];
  const angle = 0.5 * Math.atan2(2 * b, a - c), r = Math.hypot((a - c) / 2, b);
  const big = (a + c) / 2 + r, small = (a + c) / 2 - r;   // the two curvatures, big along angle
  const s = view.dpr * view.side / (2 * lnd.half);       // device pixels per map unit
  const [ox, oy] = toPixel(lnd, 0, 0);                    // the map's origin, in CSS px
  ctx.strokeStyle = color;
  ctx.lineWidth = 1;
  for (const level of lnd.levels) {
    ctx.setTransform(s, 0, 0, -s, ox * view.dpr, oy * view.dpr);   // draw in map units, y up
    ctx.beginPath();
    ctx.ellipse(lnd.pmin[0], lnd.pmin[1], Math.sqrt(2 * level / big), Math.sqrt(2 * level / small), angle, 0, 2 * Math.PI);
    ctx.setTransform(view.dpr, 0, 0, view.dpr, 0, 0);   // back to CSS px, so the line width stays even
    ctx.stroke();
  }
}

function drawInset(lnd, css) {                           // the data and the ball's line y = a + b x
  const w = 0.38 * view.side, h = 0.75 * w, x0 = view.side - w - 8, y0 = 8;
  const xs = [0.9, 2.1], ys = [Math.min(...lnd.y) - 1, Math.max(...lnd.y) + 1];
  const px = (x) => x0 + ((x - xs[0]) / (xs[1] - xs[0])) * w;
  const py = (y) => y0 + h - ((y - ys[0]) / (ys[1] - ys[0])) * h;
  ctx.save();
  ctx.fillStyle = css("--bg");
  ctx.globalAlpha = 0.92;
  ctx.fillRect(x0, y0, w, h);
  ctx.globalAlpha = 1;
  ctx.beginPath();
  ctx.rect(x0, y0, w, h);
  ctx.clip();
  const xbar = lnd.x.reduce((t, v) => t + v, 0) / lnd.x.length;
  ctx.strokeStyle = css("--muted");
  ctx.setLineDash([4, 3]);
  ctx.beginPath(); ctx.moveTo(px(xbar), y0); ctx.lineTo(px(xbar), y0 + h); ctx.stroke();
  ctx.setLineDash([]);
  ctx.fillStyle = css("--ink");
  lnd.x.forEach((x, i) => { ctx.beginPath(); ctx.arc(px(x), py(lnd.y[i]), 2.5, 0, 2 * Math.PI); ctx.fill(); });
  const [a, b] = state.p;
  ctx.strokeStyle = css("--accent");
  ctx.lineWidth = 2;
  ctx.beginPath(); ctx.moveTo(px(xs[0]), py(a + b * xs[0])); ctx.lineTo(px(xs[1]), py(a + b * xs[1])); ctx.stroke();
  const label = "y = a + b x";                          // in the lower right corner, below the data
  ctx.font = "14px system-ui, sans-serif";
  const tw = ctx.measureText(label).width;
  ctx.fillStyle = css("--bg");                          // a patch, so that the x̄ line does not cross it
  ctx.fillRect(x0 + w - tw - 9, y0 + h - 22, tw + 6, 19);
  ctx.fillStyle = css("--accent");
  ctx.textAlign = "right";
  ctx.fillText(label, x0 + w - 6, y0 + h - 7);
  ctx.restore();
}

function statusText(s) {
  const k = s.k.toLocaleString("en-US");
  if (s.status === "arrived") return `at the minimum after ${k} steps`;
  if (s.status === "left") return `left the map after ${k} steps`;
  if (s.status === "stopped") {
    const m = DATA.lands[s.land].pmin;
    return `stopped after ${k} steps, ${sig3(Math.hypot(s.p[0] - m[0], s.p[1] - m[1]))} from the minimum`;
  }
  return s.status === "running" ? "rolling" : s.status === "paused" ? "paused" : "ready";
}

function render() {
  const lnd = DATA.lands[state.land], root = getComputedStyle(document.documentElement);
  const css = (name) => root.getPropertyValue(name).trim();
  const n = canvas.width, key = `${state.land} ${n}`;
  if (cache.key !== key) cache = {key, image: relief(lnd, n)};
  ctx.setTransform(1, 0, 0, 1, 0, 0);
  ctx.putImageData(cache.image, 0, 0);
  drawContours(lnd, css("--muted"));

  // the minimum, the trail, dots at the first steps, and the ball
  const [mx, my] = toPixel(lnd, ...lnd.pmin);
  ctx.strokeStyle = MAP_INK;
  ctx.lineWidth = 2;
  ctx.beginPath(); ctx.moveTo(mx - 6, my - 6); ctx.lineTo(mx + 6, my + 6); ctx.moveTo(mx - 6, my + 6); ctx.lineTo(mx + 6, my - 6); ctx.stroke();
  const pts = state.trail.map((p) => toPixel(lnd, ...p));
  ctx.strokeStyle = ctx.fillStyle = MAP_ACCENT;
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  pts.forEach(([x, y], i) => (i ? ctx.lineTo(x, y) : ctx.moveTo(x, y)));
  ctx.stroke();
  for (const [x, y] of pts.slice(0, DOTS + 1)) { ctx.beginPath(); ctx.arc(x, y, 2.5, 0, 2 * Math.PI); ctx.fill(); }
  const [bx, by] = toPixel(lnd, ...state.p);
  ctx.beginPath(); ctx.arc(bx, by, 7, 0, 2 * Math.PI); ctx.fill();
  ctx.strokeStyle = "#ffffff";
  ctx.lineWidth = 2;
  ctx.stroke();
  if (lnd.x) drawInset(lnd, css);

  // axis names at the lower and left edges, and the hint until the first Go or tap
  ctx.font = "14px system-ui, sans-serif";
  ctx.fillStyle = css("--muted");
  ctx.textAlign = "center";
  ctx.fillText(lnd.axes[0], view.side / 2, view.side - 6);
  ctx.save(); ctx.translate(14, view.side / 2);
  if (lnd.axes[1].length > 1) ctx.rotate(-Math.PI / 2);   // a single letter stays upright
  ctx.fillText(lnd.axes[1], 0, 0); ctx.restore();
  if (state.hint) {
    ctx.fillStyle = MAP_INK;
    ctx.fillText("tap the map to place the ball", view.side / 2, view.side - 26);
  }
  ctx.textAlign = "left";

  // the panel
  const f = height(lnd, state.p[0] - lnd.pmin[0], state.p[1] - lnd.pmin[1]) + lnd.fmin;
  $("count").textContent = `step ${state.k.toLocaleString("en-US")}`;
  $("loss").textContent = `loss ${Number.isFinite(f) ? sig3(f) : "∞"}`;
  $("status").textContent = statusText(state);
  const parts = [`curvatures ${sig3(lnd.lmin)} and ${sig3(lnd.L)}`, `κ = ${Math.round(lnd.kappa)}`, `2/L = ${(2 / lnd.L).toFixed(3)}`];
  $("info").textContent = parts.map((t) => t.replaceAll(" ", "\u00a0")).join(" · ");   // breaks only between parts
  const e = lnd.eta, digits = Math.round(-Math.log10(e.step));
  Object.assign($("eta"), {min: e.min, max: e.max, step: e.step});
  $("eta").value = state.eta;
  $("eta-val").textContent = state.eta.toFixed(digits);
  // the thumb's center runs from thumb/2 to 100% − thumb/2; the 2 px tick is centered on 2/L
  $("tick").style.left = `calc(${(2 / lnd.L - e.min) / (e.max - e.min)} * (100% - var(--thumb)) + var(--thumb) / 2 - 1px)`;
  $("beta").value = state.beta;
  $("beta-val").textContent = state.beta.toFixed(2);
  $("go-label").textContent = state.status === "running" ? "Pause" : "Go";
  document.querySelectorAll("#lands button").forEach((b, i) => b.classList.toggle("on", i === state.land));
}

// ---- controls, keys, and the loop
function set(partial) {                                 // run: n replays n steps from the start at once
  const {land, run, ...rest} = partial;
  if (land !== undefined && land !== state.land) switchLand(land);
  Object.assign(state, rest);
  if (run) {
    place(state.start);
    state.status = "running";
    while (state.status === "running" && state.k < run) step(state);
    if (state.status === "running") state.status = "paused";
    state.hint = false;
  }
  render();
}

const clamp = (x, lo, hi) => Math.min(hi, Math.max(lo, x));
function setEta(eta) {
  const e = DATA.lands[state.land].eta;
  set({eta: Number(clamp(eta, e.min, e.max).toFixed(Math.round(-Math.log10(e.step))))});
}
const setBeta = (beta) => set({beta: Number(clamp(beta, 0, 0.95).toFixed(2))});
const nudgeEta = (sign) => setEta(state.eta + sign * DATA.lands[state.land].eta.step);

function go() {
  if (state.status === "running") return set({status: "paused"});
  if (["arrived", "left", "stopped"].includes(state.status)) place(state.start);   // Go again replays
  set({status: "running", hint: false});
}
function oneStep() {
  if (["arrived", "left", "stopped"].includes(state.status)) place(state.start);
  step(state);
  if (state.status === "ready" || state.status === "running") state.status = "paused";
  set({hint: false});
}
const reset = () => { place(state.start); render(); };
const toggleLand = (i = 1 - state.land) => { switchLand(i); render(); };

$("eta").addEventListener("input", (e) => setEta(Number(e.target.value)));
$("beta").addEventListener("input", (e) => setBeta(Number(e.target.value)));
$("eta-down").addEventListener("click", () => nudgeEta(-1));
$("eta-up").addEventListener("click", () => nudgeEta(1));
$("beta-down").addEventListener("click", () => setBeta(state.beta - 0.05));
$("beta-up").addEventListener("click", () => setBeta(state.beta + 0.05));
$("go").addEventListener("click", go);
$("step").addEventListener("click", oneStep);
$("reset").addEventListener("click", reset);
document.querySelectorAll("#lands button").forEach((b, i) => b.addEventListener("click", () => toggleLand(i)));

// A tap places the ball; a finger that moves 10 px or more is scrolling the page instead.
let down = null;
canvas.addEventListener("pointerdown", (e) => { down = [e.clientX, e.clientY]; });
canvas.addEventListener("pointercancel", () => { down = null; });
canvas.addEventListener("pointerup", (e) => {
  if (!down || Math.hypot(e.clientX - down[0], e.clientY - down[1]) >= 10) return;
  down = null;
  const lnd = DATA.lands[state.land], r = canvas.getBoundingClientRect(), s = (2 * lnd.half) / view.side;
  place([lnd.pmin[0] - lnd.half + (e.clientX - r.left) * s, lnd.pmin[1] + lnd.half - (e.clientY - r.top) * s]);
  set({hint: false});
});

document.addEventListener("keydown", (e) => {
  const actions = {" ": go, s: oneStep, r: reset, l: () => toggleLand(),
    ArrowUp: () => nudgeEta(1), ArrowDown: () => nudgeEta(-1),
    ArrowRight: () => setBeta(state.beta + 0.05), ArrowLeft: () => setBeta(state.beta - 0.05)};
  const act = actions[e.key.length === 1 ? e.key.toLowerCase() : e.key];
  const native = (e.target instanceof HTMLInputElement && e.key.startsWith("Arrow"))
    || (e.target instanceof HTMLButtonElement && e.key === " ");   // the focused control handles these
  if (act && !native) {
    e.preventDefault();
    act();
  }
});

let last = null;
function frame(now) {                                   // called by the browser once per screen refresh
  const dt = last === null ? 0 : (now - last) / 1000;
  last = now;
  if (state.status === "running") {
    tick(dt);
    render();
  }
  requestAnimationFrame(frame);
}

place(DATA.lands[0].start);
window.addEventListener("resize", layout);
layout();
layout();                                               // again, now that the panel holds its text
requestAnimationFrame(frame);
window.game = {state, set, render};

A step size above 2/L throws the ball out of the valley

At η = 0.019 each hop across the floor lands a little lower on the opposite wall, at 0.021 a little higher.

Show code

import json
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"

TOL, MAX_STEPS, HALF = 1e-3, 20_000, 2.5            # "at the minimum", the cap, the map's half-width
ANGLE = np.radians(30)                              # the valley floor's direction on the map

# ---- landscape 1: f = ½(u² + 100 w²), u along the floor and w across it, turned by 30°
R = np.array([[np.cos(ANGLE), -np.sin(ANGLE)],      # columns: along the floor, across it
              [np.sin(ANGLE), np.cos(ANGLE)]])
H_valley = R @ np.diag([1.0, 100.0]) @ R.T
start_valley = R @ np.array([-1.6, 0.3])            # (u, w) = (-1.6, 0.3)

# ---- landscape 2: a straight line y = a + b x through 20 calibration points
rng = np.random.default_rng(359)
x = np.linspace(1.0, 2.0, 20)
y = 0.5 + 2.0 * x + 0.1 * rng.standard_normal(x.size)
X = np.column_stack([np.ones_like(x), x])
H_fit = X.T @ X / x.size                            # Hessian of the cost (1/2n) Σ (a + b x - y)²
pmin_fit = np.linalg.lstsq(X, y, rcond=None)[0]
fmin_fit = 0.5 * np.mean((X @ pmin_fit - y) ** 2)

LANDS = [
    dict(name="Valley", H=H_valley, pmin=np.zeros(2), fmin=0.0, start=start_valley,
         axes=["x", "y"], eta=dict(min=0.001, max=0.025, step=0.001, default=0.019)),
    dict(name="Line fit", H=H_fit, pmin=pmin_fit, fmin=fmin_fit, start=pmin_fit + np.array([-2.0, 1.0]),
         axes=["intercept a", "slope b"], eta=dict(min=0.01, max=0.75, step=0.01, default=0.58)),
]
for land in LANDS:                                  # 14 contour levels c·k², the outermost touching the square
    top = HALF**2 / (2 * np.max(np.diag(np.linalg.inv(land["H"]))))
    land["levels"] = top * (np.arange(1, 15) / 14) ** 2
valley, fit = LANDS

def descend(land, eta, beta=0.0, p0=None):
    """The game's step() repeated, same arithmetic in the same order: the trail and how it ended."""
    (h00, h01), (h10, h11) = land["H"].tolist()
    m0, m1 = land["pmin"].tolist()
    px, py = land["start"].tolist() if p0 is None else p0
    vx = vy = 0.0
    trail, status = [(px, py)], "stopped"
    for k in range(MAX_STEPS):
        dx, dy = px - m0, py - m1
        gx, gy = h00 * dx + h01 * dy, h10 * dx + h11 * dy        # ∇f = H (p - p*)
        vx, vy = beta * vx - eta * gx, beta * vy - eta * gy      # v ← βv - η∇f
        px, py = px + vx, py + vy                                # p ← p + v
        trail.append((px, py))
        ex, ey = px - m0, py - m1
        if ex * ex + ey * ey < TOL * TOL:
            status = "arrived"
            break
        if max(abs(ex), abs(ey)) > HALF:
            status = "left"
            break
    return np.array(trail), status

def steps(land, eta, beta=0.0):
    trail, status = descend(land, eta, beta)
    return len(trail) - 1, status

print(f"valley: start (x, y) = ({start_valley[0]:.3f}, {start_valley[1]:.3f}), "
      f"{np.linalg.norm(start_valley):.3f} from the minimum; line fit: a = {pmin_fit[0]:.3f}, b = {pmin_fit[1]:.3f}")

valley: start (x, y) = (-1.536, -0.540), 1.628 from the minimum; line fit: a = 0.378, b = 2.105

Show code

# ---- the map, by the game's rule: cividis of log height, lightened, times a hillshade of √(f - f*)
CIVIDIS = [colors.to_hex(c) for c in plt.get_cmap("cividis")(np.linspace(0, 1, 64))]
CIVIDIS_RGB = np.array([[int(h[i:i + 2], 16) for i in (1, 3, 5)] for h in CIVIDIS], dtype=float)
DECADES, RELIEF, DOTS = 5, 0.1, 60
LIGHT = np.array([-1.0, 1.0, 1.0]) / np.sqrt(3)     # from the upper left

def height(land, px, py):
    (a, b), (_, c) = land["H"]
    dx, dy = px - land["pmin"][0], py - land["pmin"][1]
    return 0.5 * (a * dx * dx + 2 * b * dx * dy + c * dy * dy)

def relief(land, xs, ys):
    """RGB image, rows ys (top first) by columns xs, in map units."""
    (a, b), (_, c) = land["H"]
    dx, dy = xs[None, :] - land["pmin"][0], ys[:, None] - land["pmin"][1]
    h = height(land, xs[None, :], ys[:, None])
    top = max(height(land, *(land["pmin"] + HALF * np.array(corner))) for corner in [(1, 1), (1, -1)])
    with np.errstate(divide="ignore"):
        t = np.clip(1 + np.log10(h / top) / DECADES, 0, 1)
    r = RELIEF / np.where(h > 0, 2 * np.sqrt(h), 1)
    nx, ny = -r * (a * dx + b * dy), -r * (b * dx + c * dy)
    shade = 1 + 0.3 * ((nx * LIGHT[0] + ny * LIGHT[1] + LIGHT[2]) / np.sqrt(nx**2 + ny**2 + 1) - LIGHT[2])
    rgb = CIVIDIS_RGB[np.round(63 * t).astype(int)]
    return np.clip((rgb + 0.4 * (255 - rgb)) * shade[..., None], 0, 255).astype(np.uint8)

def draw_map(ax, land, trails, n=360, window=None):
    """The game's map on ax; trails is a list of (trail, color, line width, number of dots).
    window = (x0, x1, y0, y1) shows a part of the map, in map units; the default is the game's square."""
    m0, m1 = land["pmin"]
    ext = window or (m0 - HALF, m0 + HALF, m1 - HALF, m1 + HALF)
    scale = n / max(ext[1] - ext[0], ext[3] - ext[2])  # pixels per map unit
    xs = np.linspace(ext[0], ext[1], round(scale * (ext[1] - ext[0])))
    ys = np.linspace(ext[3], ext[2], round(scale * (ext[3] - ext[2])))
    ax.imshow(relief(land, xs, ys), extent=ext, interpolation="bilinear")
    gx, gy = np.meshgrid(xs, ys)
    ax.contour(gx, gy, height(land, gx, gy), levels=land["levels"], colors=MUTED, linewidths=0.7)
    ax.plot(m0, m1, "x", color=INK, ms=7, mew=2)
    for trail, color, lw, dots in trails:
        ax.plot(*trail.T, color=color, lw=lw)
        ax.plot(*trail[:dots + 1].T, "o", ms=2.5, color=color)
        ax.plot(*trail[-1], "o", ms=8, color=color, mec="white", mew=1.5)
    ax.set(xlim=ext[:2], ylim=ext[2:], aspect="equal", xticks=[], yticks=[])
    ax.set_xlabel(land["axes"][0], color=MUTED)
    ax.set_ylabel(land["axes"][1], color=MUTED)
    ax.spines[:].set_visible(False)
    ax.grid(False)
    return ax

def across(trail):                                   # w, the offset across the valley floor
    return trail @ R[:, 1]

def along(trail):                                    # u, the distance along it
    return trail @ R[:, 0]

Show code

fig, axes = plt.subplots(2, 2, figsize=(8, 7), gridspec_kw=dict(height_ratios=[1.35, 1], hspace=0.25))
for col, (eta, label) in enumerate([(0.019, "factor −0.9"), (0.021, "factor −1.1, left after 23 steps")]):
    trail, status = descend(valley, eta)
    first = trail[:26]
    draw_map(axes[0, col], valley, [(first, ACCENT, 1.2, 25)])
    axes[0, col].text(0.03, 0.97, f"η = {eta}: {label}", transform=axes[0, col].transAxes, va="top",
                      color=ACCENT, fontsize=10, bbox=dict(fc="white", ec="none", alpha=0.85, pad=2))
    ax = axes[1, col]
    k, w = np.arange(len(first)), across(first)
    ax.vlines(k, 0, w, color=INK, lw=1.2)
    ax.plot(k, w, "o", ms=4, color=INK)
    envelope = 0.3 * abs(1 - eta * 100) ** np.linspace(0, 25, 200)
    for sign in (1, -1):
        ax.plot(np.linspace(0, 25, 200), sign * envelope, ls="--", lw=1, color=MUTED)
    ax.set(xlabel="step k", xlim=(-0.5, 25.5), ylim=(-3, 3))
axes[1, 0].set_ylabel("offset across the floor w")
plt.show()

Top: the valley map with the first 25 hops at step size 0.019, a zigzag across the floor that dies down, and at 0.021, a zigzag that grows until the ball leaves after 23 steps. Bottom: the offset across the floor w against step k, shrinking by a factor 0.9 per step and growing by 1.1.

Cut the valley straight across and you get a parabola, ½ L w², with w the offset across the floor and L = 100 its curvature. A hop subtracts η times the slope L w, which multiplies the offset by a fixed factor:

\[w_{k+1} = (1 - \eta L)\, w_k.\]

At η = 0.019 the factor is −0.9: the ball overshoots the floor but lands 10 % closer each time. At 0.021 it is −1.1, and every overshoot outgrows the last. The factor stays between −1 and 1 only for η < 2/L = 0.020.

In two dimensions the gradient ∇f is the pair of slopes (∂f/∂x, ∂f/∂y), an arrow pointing straight uphill, as long as the slope is steep. Gradient descent steps against it:

\[p_{k+1} = p_k - \eta\, \nabla f(p_k).\]

In the valley's own coordinates, u along the floor and w across it, the height is

\[f = \tfrac12\,(u^2 + 100\, w^2),\]

so the gradient is (u, 100 w), and the rule above applies to each direction alone, with its own curvature: \(\lambda_{\min}\) = 1 along the floor, L = 100 across.

Every tilted quadratic bowl has two such perpendicular directions, and the eigenvalues of its Hessian H, the 2 × 2 table of second derivatives of f, are the curvatures along them. The largest sets the limit, because a step that is stable along the floor can still be unstable across it.

The cell below gets both curvatures from H with eigvalsh, NumPy's eigenvalues of a symmetric matrix. Then steps replays the game's own hop, shown in How the game computes it, until the ball is within 10⁻³ of the minimum or off the map.

curvatures = np.linalg.eigvalsh(H_valley)           # ascending: along the floor, across it
lam_min, L = curvatures
print(f"curvatures {lam_min:.0f} and {L:.0f}, limit 2/L = {2 / L:.3f}")
for eta in [0.01, 0.019, 0.021]:
    print(f"η = {eta:5.3f}: factor across 1 - ηL = {1 - eta * L:+.1f}")
for eta in [0.019, 0.021]:
    k, status = steps(valley, eta)
    print(f"η = {eta}: {'at the minimum' if status == 'arrived' else 'left the map'} after {k} steps")

curvatures 1 and 100, limit 2/L = 0.020
η = 0.010: factor across 1 - ηL = +0.0
η = 0.019: factor across 1 - ηL = -0.9
η = 0.021: factor across 1 - ηL = -1.1
η = 0.019: at the minimum after 385 steps
η = 0.021: left the map after 23 steps

Play: roll a ball down a narrow valley

A step size above 2/L throws the ball out of the valley

The second half of this tutorial needs a free account