26. Learning From Data: Gradient Descent in Mountain Goat¶

Mountain goats on the mountain at sunset

What you will understand: how a machine-learning model is trained, from the ground up, and how to do it with a language that has matrices but no loops. You will fit a straight line (a linear regression) to data by gradient descent: the repeated, small, downhill steps that train almost every neural network. The mathematics is written in Mountain Goat and compiled; a short C++ program supplies the loop, the data and the checks. The result is checked against a closed-form answer and against an independent C++ implementation, and the tests are tested by breaking the program seven ways.

What you need to know first: the language tour (matrices, @, transpose, broadcasting), Chapter 21 (the arithmetic operations), Chapter 22 (reductions such as col_sum) and Chapter 20 (calling compiled Mountain Goat from C++). No machine-learning background is assumed.

Compile and run

cd docs/part24/code && ./build.sh               # once: the newest compiler (Chapter 24's mg-opt)
cd ../../part26/code
cpp/run.sh                                      # compile gd.mg, build the C++ trainer against it, run it (-> cpp/run_out.txt)
./show.sh                                       # the MLIR the front end makes of gd.mg (-> show/)
./gd_mutation.sh > gd_mutation_out.txt          # breaks gd.mg on purpose, seven ways
cd ../../part15/code && ./run_lit.sh            # the whole test suite

Every listing and output on this page comes from these commands.

Why this page shows FAIL lines, and why that is good

The mutation output near the end comes from a negative control: gd.mg deliberately broken and run against the test. There a "caught" line (the test failed) is the expected, wanted result. The baseline line must say the test passes, and it does. The first part of the output (the loss falling, the parameters recovered) is the real, correct program.

Primer: learning as going downhill¶

The problem. You have data: for each of 200 samples, three measurements (the features) and one number to predict (the target). You believe the target is roughly a weighted sum of the features plus a constant: target ≈ w1·f1 + w2·f2 + w3·f3 + b. The weights w1, w2, w3 and the bias b are the parameters. Finding the parameters that fit the data best is linear regression.

A trick to keep it simple. Add a fourth column to the data that is the number 1 for every sample. Then the bias b is just a fourth weight (the weight of a feature that is always 1), and the whole model is one matrix product: predictions = x @ p, where x has one row per sample and four columns, and p is a column of four parameters. That is why the code below has tensor[?x4] for the data and tensor[4x1] for the parameters.

Measuring how wrong a model is. The residual of a sample is its prediction minus its target. Square each residual (so errors of either sign count, and big errors count a lot), add them up and divide by the number of samples: the mean squared error, or loss. The loss is a single number, and it depends on the parameters: it is high for bad parameters and zero for perfect ones.

Going downhill. Picture the loss as the height of a landscape whose coordinates are the four parameters. The best parameters are at the bottom of the valley. From wherever you stand, the gradient is the direction in which the loss rises fastest, so stepping the opposite way goes downhill. Gradient descent repeats one rule: new parameters = old parameters − (learning rate) × gradient. The learning rate is the size of each step. For this loss the gradient has a clean form: gradient = (2 / N) · xᵀ (x p − y), where xᵀ is the transpose of x, x p − y is the column of residuals and N is the number of samples. Read it as: for each parameter, add up (residual × that feature) over all samples.

Too small, too large. With a small learning rate the steps are safe but slow. With one that is too large, each step overshoots the valley, lands higher on the far side, and the loss grows; gradient descent diverges. Part C of the output shows exactly that.

What the language cannot do. Gradient descent repeats a step hundreds of times, and Mountain Goat has no loops. Rather than add one, this chapter splits the work along a clear line: the math of one step (gradient, update, loss) is Mountain Goat, compiled; the loop that repeats it, the data and the checks are a small C++ program that calls the compiled functions (as in Chapter 20).

The Mountain Goat part¶

# Linear regression by gradient descent: the math, in Mountain Goat. The loop that repeats it is in the C++ driver (the language has none).
#
# The model is  prediction = x @ p,  where x has one row per sample and four columns (three features and a constant 1 for the bias),
# and p is the 4x1 column of parameters. The number of samples is dynamic (?); the number of columns is fixed at 4.
#
# grad  = x^T (x p - y)                          (half the gradient of the summed squared error)
# step  = p - rate * grad                        (one gradient-descent update; rate is passed as a 1x1 matrix: lr * 2 / number_of_samples)
# loss_sum = the sum of squared residuals        (1x1)
def grad(x: tensor[?x4], y: tensor[?x1], p: tensor[4x1]) = transpose(x) @ (x @ p - y)
def step(x: tensor[?x4], y: tensor[?x1], p: tensor[4x1], rate: tensor[1x1]) = p - rate * grad(x, y, p)
def loss_sum(x: tensor[?x4], y: tensor[?x1], p: tensor[4x1]) = col_sum((x @ p - y) * (x @ p - y))

Three functions. grad is xᵀ (x p − y) written almost as it is on paper: transpose(x) @ (x @ p - y). step is the update rule, p - rate * grad(x, y, p); the factor 2 / N and the learning rate are folded into one number, rate, which the caller passes as a 1×1 matrix (so the function works for any number of samples: N appears nowhere in the code). loss_sum squares the residuals elementwise and adds them with col_sum; the caller divides by N.

The typed parameters carry the shapes: x: tensor[?x4] is "any number of rows, exactly four columns", so one compiled function serves any dataset size; p: tensor[4x1] and rate: tensor[1x1] are fixed. Note that rate * grad(...) multiplies a 1×1 matrix by a 4×1 one: that is broadcasting (Chapter 22), and it requires static shapes, which is why p and rate have them. The front end turns the file into this MLIR (Chapter 20); everything from Chapters 21 to 24 now lowers it:

func.func @grad(%x: tensor<?x4xf64>, %y: tensor<?x1xf64>, %p: tensor<4x1xf64>) -> tensor<4x1xf64> {
  %t1 = mg.transpose %x : tensor<?x4xf64> to tensor<4x?xf64>
  %t2 = mg.matmul %x, %p : tensor<?x4xf64>, tensor<4x1xf64> -> tensor<?x1xf64>
  %t3 = mg.sub %t2, %y : tensor<?x1xf64>, tensor<?x1xf64> -> tensor<?x1xf64>
  %t4 = mg.matmul %t1, %t3 : tensor<4x?xf64>, tensor<?x1xf64> -> tensor<4x1xf64>
  func.return %t4 : tensor<4x1xf64>
}
func.func @step(%x: tensor<?x4xf64>, %y: tensor<?x1xf64>, %p: tensor<4x1xf64>, %rate: tensor<1x1xf64>) -> tensor<4x1xf64> {
  %t1 = func.call @grad(%x, %y, %p) : (tensor<?x4xf64>, tensor<?x1xf64>, tensor<4x1xf64>) -> tensor<4x1xf64>
  %t2 = mg.broadcast %rate : tensor<1x1xf64> -> tensor<4x1xf64>
  %t3 = mg.mul %t2, %t1 : tensor<4x1xf64>, tensor<4x1xf64> -> tensor<4x1xf64>
  %t4 = mg.sub %p, %t3 : tensor<4x1xf64>, tensor<4x1xf64> -> tensor<4x1xf64>
  func.return %t4 : tensor<4x1xf64>
}
func.func @loss_sum(%x: tensor<?x4xf64>, %y: tensor<?x1xf64>, %p: tensor<4x1xf64>) -> tensor<1x1xf64> {
  %t1 = mg.matmul %x, %p : tensor<?x4xf64>, tensor<4x1xf64> -> tensor<?x1xf64>
  %t2 = mg.sub %t1, %y : tensor<?x1xf64>, tensor<?x1xf64> -> tensor<?x1xf64>
  %t3 = mg.matmul %x, %p : tensor<?x4xf64>, tensor<4x1xf64> -> tensor<?x1xf64>
  %t4 = mg.sub %t3, %y : tensor<?x1xf64>, tensor<?x1xf64> -> tensor<?x1xf64>
  %t5 = mg.mul %t2, %t4 : tensor<?x1xf64>, tensor<?x1xf64> -> tensor<?x1xf64>
  %t6 = mg.reduce %t5 {axis = 0 : i64, kind = "sum"} : tensor<?x1xf64> -> tensor<1x1xf64>
  func.return %t6 : tensor<1x1xf64>
}

Look at @step: the mg.broadcast is inserted automatically for rate, and the result of grad flows into the elementwise mg.mul and mg.sub. The transpose, two matrix products, a subtraction, a broadcast, a multiplication and a subtraction: that is one gradient-descent step.

The C++ part¶

The driver generates data, runs the loop, and checks the answers:

// Gradient descent for linear regression. The MATH (gradient, update, loss) runs in compiled Mountain Goat (gd.mg);
// the LOOP, the data, and every check live here in C++, because the language has no loops.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <vector>
#include "gd.h"

// A tiny deterministic random number generator, so the data (and the output) are the same on every run and machine.
static uint64_t state = 88172645463325252ULL;
static double uniform(double lo, double hi) {
  state ^= state << 13; state ^= state >> 7; state ^= state << 17;
  return lo + (hi - lo) * ((state >> 11) * (1.0 / 9007199254740992.0));
}

static const int N = 200, D = 4;                         // 200 samples; 3 features + a constant-1 column for the bias
static const double TRUE_P[D] = {2.0, -3.0, 0.5, 1.0};    // the parameters the data is generated from (last one is the bias)

struct Data { mg::Matrix x, y; };
static Data make_data(double noise) {
  std::vector<double> x(N * D), y(N);
  for (int i = 0; i < N; i++) {
    for (int j = 0; j < D - 1; j++) x[i * D + j] = uniform(-1.0, 1.0);
    x[i * D + (D - 1)] = 1.0;
    double t = 0; for (int j = 0; j < D; j++) t += x[i * D + j] * TRUE_P[j];
    y[i] = t + (noise > 0 ? uniform(-noise, noise) : 0.0);
  }
  return {mg::Matrix(N, D, x), mg::Matrix(N, 1, y)};
}

// The same gradient descent written as plain C++ loops, as an independent check of the compiled version.
static std::vector<double> cpp_descent(const Data &d, double lr, int steps) {
  std::vector<double> p(D, 0.0);
  for (int s = 0; s < steps; s++) {
    std::vector<double> g(D, 0.0);
    for (int i = 0; i < N; i++) {
      double r = -d.y.data[i]; for (int j = 0; j < D; j++) r += d.x.data[i * D + j] * p[j];
      for (int j = 0; j < D; j++) g[j] += d.x.data[i * D + j] * r;
    }
    for (int j = 0; j < D; j++) p[j] -= (2.0 * lr / N) * g[j];
  }
  return p;
}

// The exact least-squares answer from the normal equations (X^T X) p = X^T y, by Gaussian elimination.
static std::vector<double> closed_form(const Data &d) {
  double a[D][D + 1] = {};
  for (int i = 0; i < N; i++) for (int j = 0; j < D; j++) {
    for (int k = 0; k < D; k++) a[j][k] += d.x.data[i * D + j] * d.x.data[i * D + k];
    a[j][D] += d.x.data[i * D + j] * d.y.data[i];
  }
  for (int c = 0; c < D; c++) {
    int piv = c; for (int r = c + 1; r < D; r++) if (std::fabs(a[r][c]) > std::fabs(a[piv][c])) piv = r;
    for (int k = 0; k <= D; k++) std::swap(a[c][k], a[piv][k]);
    for (int r = 0; r < D; r++) if (r != c) { double f = a[r][c] / a[c][c]; for (int k = c; k <= D; k++) a[r][k] -= f * a[c][k]; }
  }
  std::vector<double> p(D); for (int j = 0; j < D; j++) p[j] = a[j][D] / a[j][j];
  return p;
}

// The mean squared error computed with plain loops, to check the compiled loss_sum independently.
static double cpp_loss(const Data &d, const std::vector<double> &p) {
  double sum = 0;
  for (int i = 0; i < N; i++) { double r = -d.y.data[i]; for (int j = 0; j < D; j++) r += d.x.data[i * D + j] * p[j]; sum += r * r; }
  return sum / N;
}
static double max_abs_diff(const std::vector<double> &a, const std::vector<double> &b) {
  double m = 0; for (size_t i = 0; i < a.size(); i++) m = std::max(m, std::fabs(a[i] - b[i])); return m;
}
static double mean_loss(const Data &d, const mg::Matrix &p) { return mg::loss_sum(d.x, d.y, p).data[0] / N; }
static void show(const char *label, const std::vector<double> &p) {
  std::printf("  %-34s", label); for (double v : p) std::printf(" %9.5f", v); std::printf("\n");
}

// Run `steps` updates in compiled Mountain Goat, returning the final parameters; losses[s] is the mean squared error before step s.
static mg::Matrix descend(const Data &d, double lr, int steps, std::vector<double> &losses,
                          std::vector<std::vector<double>> *trajectory = nullptr) {   // trajectory[s] = parameters after s steps
  mg::Matrix p(D, 1, std::vector<double>(D, 0.0));
  mg::Matrix rate(1, 1, {2.0 * lr / N});                 // the update is  p - (2 lr / N) * x^T (x p - y)
  for (int s = 0; s <= steps; s++) {
    losses.push_back(mean_loss(d, p));
    if (trajectory) trajectory->push_back(p.data);
    if (s < steps) p = mg::step(d.x, d.y, p, rate);
  }
  return p;
}

int main() {
  const double lr = 0.2; const int steps = 300;
  std::printf("true parameters (w1 w2 w3 bias):  ");  for (double v : TRUE_P) std::printf(" %9.5f", v); std::printf("\n\n");

  // ---- A: no noise. Gradient descent must recover the true parameters, and the loss must fall at every step. ----
  std::printf("A. noise-free data, learning rate %.2f, %d steps\n", lr, steps);
  Data a = make_data(0.0);
  std::vector<double> la; std::vector<std::vector<double>> traj; mg::Matrix pa = descend(a, lr, steps, la, &traj);
  std::printf("  mean squared error at steps:");
  for (int s : {0, 1, 2, 5, 10, 20, 50, 100, 200, 300}) std::printf("  [%d] %.3e", s, la[s]);
  std::printf("\n");
  show("learned (Mountain Goat)", pa.data);
  bool monotone = true; for (size_t s = 1; s < la.size(); s++) if (!(la[s] < la[s - 1]) && la[s - 1] > 1e-25) monotone = false;
  std::printf("  loss fell at every step: %s\n", monotone ? "yes" : "NO");
  std::printf("  parameters recovered to within 1e-6 of the truth: %s\n", max_abs_diff(pa.data, std::vector<double>(TRUE_P, TRUE_P + D)) < 1e-6 ? "yes" : "NO");
  std::printf("  matches the plain C++ gradient descent to within 1e-9: %s\n", max_abs_diff(pa.data, cpp_descent(a, lr, steps)) < 1e-9 ? "yes" : "NO");
  // The end state alone cannot tell a correct gradient from one that is, say, twice too big (both converge to the same answer),
  // so also compare the first steps one at a time: the size of every early step depends on getting the gradient exactly right.
  bool early = true; for (int k = 1; k <= 5; k++) if (max_abs_diff(traj[k], cpp_descent(a, lr, k)) > 1e-12) early = false;
  std::printf("  the first 5 steps match plain C++ step for step to within 1e-12: %s\n", early ? "yes" : "NO");
  std::printf("  the compiled loss equals a plain C++ loss (at the start and the end): %s\n",
              std::fabs(la[0] - cpp_loss(a, std::vector<double>(D, 0.0))) < 1e-9 * la[0] && std::fabs(la.back() - cpp_loss(a, pa.data)) < 1e-12 ? "yes" : "NO");

  // ---- B: noisy data. The truth is no longer the best fit; the best fit is the least-squares (closed-form) answer. ----
  std::printf("\nB. noisy data (noise up to +-0.2), learning rate %.2f, %d steps\n", lr, steps);
  Data b = make_data(0.2);
  std::vector<double> lb; mg::Matrix pb = descend(b, lr, steps, lb);
  std::vector<double> best = closed_form(b);
  show("learned (Mountain Goat)", pb.data);
  show("exact least squares (closed form)", best);
  show("true parameters", std::vector<double>(TRUE_P, TRUE_P + D));
  std::printf("  final mean squared error: %.5f (the noise alone would give about %.5f = 0.2^2/3)\n", lb.back(), 0.2 * 0.2 / 3);
  std::printf("  gradient descent matches the exact least-squares answer to within 1e-6: %s\n", max_abs_diff(pb.data, best) < 1e-6 ? "yes" : "NO");
  std::printf("  the compiled loss equals a plain C++ loss on the noisy data: %s\n", std::fabs(lb.back() - cpp_loss(b, pb.data)) < 1e-12 ? "yes" : "NO");
  std::printf("  closer to the exact answer than to the truth: %s\n", max_abs_diff(pb.data, best) < max_abs_diff(pb.data, std::vector<double>(TRUE_P, TRUE_P + D)) ? "yes" : "NO");

  // ---- C: a learning rate that is too large makes gradient descent diverge. ----
  const double bad = 1.2;
  std::printf("\nC. noise-free data, learning rate %.2f (too large), 12 steps\n", bad);
  std::vector<double> lc; descend(a, bad, 12, lc);
  std::printf("  mean squared error at steps:"); for (int s : {0, 1, 2, 4, 8, 12}) std::printf("  [%d] %.3e", s, lc[s]); std::printf("\n");
  std::printf("  the loss grew instead of shrinking: %s\n", lc[12] > 100 * lc[0] ? "yes" : "NO");
  return 0;
}

The data is generated from known parameters, 2.0, −3.0, 0.5 and a bias of 1.0, with a small deterministic random number generator so the output is the same every run. Three experiments:

A, noise-free data. Gradient descent should recover the true parameters exactly.
B, noisy data. Now the data is not exactly a line, so the truth is no longer the best fit: the best fit is the least-squares solution, which has an exact formula (the normal equations, (xᵀx) p = xᵀy, solved here by Gaussian elimination in C++). Gradient descent should land on that, not on the truth.
C, learning rate too large. The loss should grow.

Everything else is independent checking: the driver also implements gradient descent and the loss in plain C++ loops and compares the compiled Mountain Goat results against them.

The run¶

Built with mgc lib gd.mg -O2 --matmul-order ikj (Chapters 23 and 24's optimizations), compiled against the C++ driver, and run:

#!/bin/sh
# Compile gd.mg (optimized, with the faster matmul loop order from Chapter 24), build the C++ trainer against it, run it. Output: run_out.txt
HERE=$(cd "$(dirname "$0")" && pwd); W=$HERE/../work/train; mkdir -p $W
$HERE/../../../part24/code/mgc lib $HERE/gd.mg -o $W -O2 --matmul-order ikj > /dev/null
cp $W/gd.h $HERE/gd_generated.h
clang++-18 -std=c++17 -O2 -Wall -I$W $HERE/train.cpp $W/gd.o -o $W/train
$W/train

true parameters (w1 w2 w3 bias):     2.00000  -3.00000   0.50000   1.00000

A. noise-free data, learning rate 0.20, 300 steps
  mean squared error at steps:  [0] 5.431e+00  [1] 3.607e+00  [2] 2.556e+00  [5] 1.061e+00  [10] 2.708e-01  [20] 1.843e-02  [50] 6.827e-06  [100] 1.616e-11  [200] 1.040e-22  [300] 9.548e-31
  learned (Mountain Goat)              2.00000  -3.00000   0.50000   1.00000
  loss fell at every step: yes
  parameters recovered to within 1e-6 of the truth: yes
  matches the plain C++ gradient descent to within 1e-9: yes
  the first 5 steps match plain C++ step for step to within 1e-12: yes
  the compiled loss equals a plain C++ loss (at the start and the end): yes

B. noisy data (noise up to +-0.2), learning rate 0.20, 300 steps
  learned (Mountain Goat)              1.99714  -2.98902   0.49820   0.99592
  exact least squares (closed form)    1.99714  -2.98902   0.49820   0.99592
  true parameters                      2.00000  -3.00000   0.50000   1.00000
  final mean squared error: 0.01343 (the noise alone would give about 0.01333 = 0.2^2/3)
  gradient descent matches the exact least-squares answer to within 1e-6: yes
  the compiled loss equals a plain C++ loss on the noisy data: yes
  closer to the exact answer than to the truth: yes

C. noise-free data, learning rate 1.20 (too large), 12 steps
  mean squared error at steps:  [0] 5.431e+00  [1] 2.815e+00  [2] 5.208e+00  [4] 2.098e+01  [8] 3.424e+02  [12] 5.587e+03
  the loss grew instead of shrinking: yes

Reading it.

A. The loss starts at 5.43 and falls every step: 0.27 after 10 steps, 1.8×10⁻² after 20, 7×10⁻⁶ after 50, about 10⁻³⁰ after 300 (essentially the limit of the arithmetic). The learned parameters are 2.00000, −3.00000, 0.50000, 1.00000: the truth. The compiled Mountain Goat result matches the plain C++ version to within 10⁻⁹.
B. The learned parameters, 1.99714, −2.98902, 0.49820, 0.99592, are identical to five decimals to the exact least-squares answer, and differ from the true parameters by the amount the noise pushed them. The final loss, 0.01343, is almost exactly what the noise alone would give: the noise is uniform up to ±0.2, whose mean square is 0.2²/3 = 0.01333. The model has learned everything there is to learn; what is left is noise. The program checks that gradient descent is closer to the exact answer than to the truth, which is the correct behavior for noisy data.
C. With a learning rate of 1.2 the loss goes 5.4, 2.8, 5.2, 21, 342, 5587: each step overshoots further than the last. The limit can be computed: for a loss like this one, gradient descent is stable only when the learning rate is below 2 divided by the largest curvature of the loss, and for this data that curvature is 2.0148, so the limit is 0.9927 (stability_limit.py computes it, in plain Python, by reproducing the driver's data):

largest curvature of the loss on this data: 2.0148
gradient descent is stable only for learning rates below 2 / 2.0148 = 0.9927
  learning rate 0.2: stable
  learning rate 1.2: diverges

So 0.2 is safely below the limit and 1.2 is above it. The checks print "the loss grew instead of shrinking: yes", which here is the correct behavior.

Tests¶

One new lit test, test/learning/gradient-descent, builds the trainer twice (once optimized with the ikj loop order, once with the plain defaults) and runs both. Every property the driver checks must say yes, and the two builds' complete outputs must be identical: the compiler settings change the speed, not the answer. The suite now has 103 tests.

The driver's checks are the substance of the test, and the second half of the chapter is about whether they are good enough.

Are the checks good enough? Break the program seven ways.¶

gd_mutation.sh edits gd.mg to be wrong, one way at a time, and re-runs the test:

#!/bin/sh
# READ THIS FIRST: this script edits cpp/gd.mg to be WRONG, one way at a time, and re-runs the Chapter 26 test (test/learning).
#   A "caught" line (the test failed) is the EXPECTED, wanted result: it shows the test can notice that mistake.
#   "NOT CAUGHT" would be a test gap. gd.mg is restored afterwards. Output: gd_mutation_out.txt
HERE=$(cd "$(dirname "$0")" && pwd); D=$HERE/../..; LIT=$D/part15/code/run_lit.sh; F=$HERE/cpp/gd.mg; W=$HERE/work/mut; mkdir -p $W
echo "NOTE: this script deliberately breaks gd.mg. 'caught' lines are EXPECTED: they show the test can detect the mistake."
cp $F $W/original
run() { out=$(MG_TEST_TMP=$W/out sh $LIT --filter ':: learning/' 2>&1)
  if echo "$out" | grep -q '^FAIL'; then echo "$1: caught"; else echo "$1: NOT CAUGHT (a test gap)"; fi; }
mutate() {  # label sed-expr
  cp $W/original $F; sed -i "$2" $F
  if cmp -s $F $W/original; then echo "$1: MUTATION DID NOT APPLY"; else run "$1"; fi
  cp $W/original $F
}
echo "baseline (nothing broken): $(MG_TEST_TMP=$W/out sh $LIT --filter ':: learning/' 2>&1 | grep -q '^FAIL' && echo 'UNEXPECTED FAILURE' || echo 'test passes (expected)')"
mutate "update adds the gradient instead of subtracting it"        's/^def step\(.*\) = p - rate/def step\1 = p + rate/'
mutate "gradient has the wrong sign (y - x p)"                    's/^def grad\(.*\) = transpose(x) @ (x @ p - y)/def grad\1 = transpose(x) @ (y - x @ p)/'
mutate "the transpose is missing from the gradient"                's/^def grad\(.*\) = transpose(x) @/def grad\1 = x @/'
mutate "the learning rate is ignored"                              's/p - rate \* grad(x, y, p)/p - grad(x, y, p)/'
mutate "the gradient is twice too large (converges to the same answer)" 's/^def grad\(.*\) = transpose(x) @ (x @ p - y)/def grad\1 = 2 * (transpose(x) @ (x @ p - y))/'
mutate "the loss forgets to square the residuals"                  's/col_sum((x @ p - y) \* (x @ p - y))/col_sum(x @ p - y)/'
mutate "the loss sums the wrong way (row_sum)"                     's/col_sum((x @ p - y)/row_sum((x @ p - y)/'
echo "restored: $(cmp -s $F $W/original && echo 'gd.mg is back to its original' || echo 'RESTORE FAILED')"

NOTE: this script deliberately breaks gd.mg. 'caught' lines are EXPECTED: they show the test can detect the mistake.
baseline (nothing broken): test passes (expected)
update adds the gradient instead of subtracting it: caught
gradient has the wrong sign (y - x p): caught
the transpose is missing from the gradient: caught
the learning rate is ignored: caught
the gradient is twice too large (converges to the same answer): caught
the loss forgets to square the residuals: caught
the loss sums the wrong way (row_sum): caught
restored: gd.mg is back to its original

All seven are caught. One of them deserves a closer look, because the first version of the driver would not have caught it. The mutation "the gradient is twice too large" multiplies the gradient by 2. That is the same as doubling the learning rate, and gradient descent with a larger (but still stable) learning rate converges to the same answer, just faster. So a driver that only looked at the end state (the parameters after 300 steps) would pass: the parameters are still the truth, still match plain C++ at 300 steps (both have converged to the limit of the arithmetic), and the loss still falls at every step. I checked this by running the mutated program:

  loss fell at every step: yes
  parameters recovered to within 1e-6 of the truth: yes
  matches the plain C++ gradient descent to within 1e-9: yes
  the first 5 steps match plain C++ step for step to within 1e-12: NO        <- the only check that notices
  the compiled loss equals a plain C++ loss (at the start and the end): yes
  gradient descent matches the exact least-squares answer to within 1e-6: yes
  ...

The only check that notices is the one added for exactly this reason: compare the first few steps with plain C++, one at a time. The size of every early step depends on the gradient being exactly right, while the final answer does not depend on it at all. This is a general lesson about testing iterative algorithms: checking where a process ends up does not check how it got there, and a bug can change the path while leaving the destination alone. The driver also checks the compiled loss against a plain C++ loss for the same reason: a loss function with a bug (the mutation that forgets to square the residuals) could otherwise go unnoticed, because gradient descent never looks at the loss, only at its gradient.

The full suite with the newest compiler:

Testing Time: 7.31s

Total Discovered Tests: 103
  Passed: 103 (100.00%)

Limits and what is not established¶

A toy problem. One linear model, 200 samples, four parameters, synthetic data generated from known parameters. Nothing here says Mountain Goat can train anything realistic.
The loop is in C++. The language has no loops, so every step crosses from C++ into the compiled code and back, and the wrapper copies the data in and out on each call. That is correct and simple, and it is the opposite of fast; nothing was timed.
No other models. Logistic regression and neural networks need functions the language lacks (an exponential, a way to compute the derivative of relu, which is a comparison). relu itself exists (Chapter 22) but its gradient does not.
Fixed number of features. The feature count (4) is part of the types; only the sample count is dynamic.
Full-batch gradient descent only. Each step uses all samples. Mini-batches, momentum and adaptive learning rates were not tried.
The learning-rate limit is computed once, by a separate script, not tested. Only the rates 0.2 and 1.2 were run through the real program; that rates between (say 0.9 and 1.0) behave as the computed limit predicts was not checked.
Exact equalities are checked to tolerances (1e-6 to 1e-12 as printed), not bit for bit; floating-point arithmetic in a different order gives slightly different last digits.
Platform: x86-64 Linux, LLVM 18.1.3.

Reproducing¶

cd docs/part24/code && ./build.sh
cd ../../part26/code
cpp/run.sh > cpp/run_out.txt
./show.sh
./gd_mutation.sh > gd_mutation_out.txt
cd ../../part15/code && ./run_lit.sh          # 103 tests

Chapter summary¶

Linear regression is predictions = x @ p; the bias becomes a fourth weight with a column of ones; the loss is the mean squared residual; the gradient is (2/N) xᵀ (x p − y); gradient descent repeats p ← p − rate · gradient.
The gradient, the update and the loss are three one-line Mountain Goat functions; the loop, data and checks are C++ (the language has no loops).
On noise-free data it recovers the true parameters; on noisy data it reaches the exact least-squares answer to five decimals, and the final loss is what the noise alone predicts; with a learning rate that is too large it diverges.
The results match independent plain-C++ implementations (the descent, the loss and the closed-form answer).
Seven deliberate bugs are all caught; one of them (a gradient twice too large) is invisible to end-state checks and needed a step-for-step trajectory check, because testing where an iteration ends does not test how it got there.

Self-check questions¶

Each answer is collapsed; try the question first.

Why does the data get a column of ones, and what does that do to the model?

Answer

It turns the bias (the constant term) into an ordinary weight: the weight of a feature that is always 1. Then the whole model is a single matrix product, x @ p, with no separate "add the bias" step, and the gradient formula treats all four parameters the same way.
What is the gradient, and why does gradient descent step in the opposite direction?

Answer

The gradient is the direction in which the loss rises fastest, with a length that says how steeply. To make the loss smaller you go the opposite way, so each update subtracts a multiple of the gradient from the parameters: p ← p − rate · gradient.
In gd.mg, why is the learning rate passed in as a 1×1 matrix, and why must p and rate have fixed shapes?

Answer

The language has no scalar parameters (a scalar is a compile-time number, not a value a function can receive), so a one-number value travels as a 1×1 matrix. Multiplying it by the 4×1 gradient needs broadcasting, and broadcasting works only for static shapes (Chapter 22): the compiler must know the 1×1 is a 1 that can be repeated. So p and rate are fixed-shape, while x and y keep a dynamic row count.
Why do the noisy-data results match the exact least-squares answer rather than the true parameters?

Answer

With noise, the data is no longer exactly produced by the true parameters, so the true parameters are not the ones that minimize the loss on this data. Gradient descent minimizes the loss on the data it is given, so it finds the least-squares solution, which the normal equations compute exactly. The difference between that solution and the truth is the effect of the noise on this sample.
What happens with a learning rate that is too large, and why?

Answer

Each step overshoots the bottom of the valley and lands higher on the other side, where the gradient is bigger still, so the next step overshoots by more. The loss grows from step to step (5.4, 2.8, 5.2, 21, 342, 5587 at a rate of 1.2). Below a limit determined by the shape of the loss (0.9927 for this data, computed in stability_limit.py) the steps shrink toward the bottom instead.
The mutation "the gradient is twice too large" passed every check about the final answer. Why, and which check caught it?

Answer

A gradient twice too large is the same as a learning rate twice as large; as long as that is still stable, gradient descent converges to the same minimum, so the final parameters (and a comparison with plain C++ at 300 steps, when both have converged to the limit of the arithmetic) look right. Only the early steps differ in size. The check that compares the first five steps one at a time with plain C++ caught it.
What general lesson about testing iterative algorithms does that case teach?

Answer

Checking where a process ends up does not check how it got there: a bug can change the path and leave the destination unchanged. Tests for iterations should also look at intermediate states (here, the first few steps) against an independent implementation.
The driver checks the compiled loss against a plain C++ loss. Gradient descent never uses the loss, so why bother?

Answer

Because nothing else would notice a wrong loss: the parameters come from the gradient alone. A loss function that forgot to square the residuals, or summed along the wrong axis, would leave every check about the parameters passing while the printed losses (the thing a user watches to see whether training works) were wrong. Two of the seven mutations are exactly that.
Why does the test build the trainer twice and compare the complete outputs?

Answer

The two builds use different compiler settings (optimized with the ikj loop order, and the unoptimized defaults), which change only how fast the code runs. If the outputs differ, a compiler setting changed the answer, which would be a bug. Identical output at the printed precision is evidence that the optimizations preserve the results for this program.
What would you need to add to the language to train a logistic regression or a small neural network?

Answer

An exponential (for the sigmoid), a way to take the derivative of relu (a comparison that gives 1 where the input is positive and 0 elsewhere), probably a log for the usual loss, and ideally a loop so the whole training run is in the language rather than in C++. Each new operation needs a verifier, a lowering, front-end syntax and tests, as in Chapters 21 and 22.