Project: a neural network

No library, no framework, no external data: a two-layer network trained by backpropagation, in a hundred and seventy lines of the language you have now learned all of. It solves XOR, which a single layer cannot, and then a harder problem, telling two rings of points apart, on data it has never seen.

topo/code/nn.nx
//! A neural network from nothing: a two-layer perceptron trained by
//! backpropagation, in plain Nexium. It learns two problems: XOR, the classic
//! that a single layer cannot solve, and telling two rings of points apart.

/// A dense matrix, row major.
struct Matrix { rows: usize, cols: usize, data: List(f64) }

impl Matrix {
    fn zeros(rows: usize, cols: usize) -> Matrix {
        var data = List(f64).with_capacity(rows * cols)
        for _ in 0..rows * cols { data.append(0.0) }
        return Matrix{ .rows = rows, .cols = cols, .data = data }
    }

    /// Small random weights, the usual starting point.
    fn random(rows: usize, cols: usize, scale: f64) -> Matrix {
        var m = Matrix.zeros(rows, cols)
        for i in 0..m.data.len { m.data[i] = (random.float() * 2.0 - 1.0) * scale }
        return m
    }

    fn at(self: *Self, r: usize, c: usize) -> f64 { return self.data[r * self.cols + c] }
    fn set(self: *mut Self, r: usize, c: usize, v: f64) { self.data[r * self.cols + c] = v }
}

fn sigmoid(x: f64) -> f64 { return 1.0 / (1.0 + math.exp(-x)) }

/// One layer: `outputs` neurons, each with a weight per input and a bias.
struct Layer {
    w: Matrix
    b: List(f64)
    // what the layer saw and produced in the last forward pass, for backprop
    input: List(f64)
    output: List(f64)
}

impl Layer {
    fn new(inputs: usize, outputs: usize) -> Layer {
        var b = List(f64).new()
        for _ in 0..outputs { b.append(0.0) }
        return Layer{ .w = Matrix.random(outputs, inputs, 1.0), .b = b, .input = List(f64).new(), .output = List(f64).new() }
    }

    fn forward(self: *mut Self, x: []f64) -> List(f64) {
        self.input = x.to_owned()
        var out = List(f64).with_capacity(self.w.rows)
        for r in 0..self.w.rows {
            var s = self.b[r]
            for c in 0..self.w.cols { s += self.w.at(r, c) * x[c] }
            out.append(sigmoid(s))
        }
        self.output = out.clone()
        return out
    }

    /// Given dLoss/dOutput, update the weights and return dLoss/dInput for the layer below.
    fn backward(self: *mut Self, grad_out: []f64, lr: f64) -> List(f64) {
        var grad_in = List(f64).new()
        for _ in 0..self.w.cols { grad_in.append(0.0) }
        for r in 0..self.w.rows {
            let y = self.output[r]
            let delta = grad_out[r] * y *(1.0 - y) // the sigmoid's derivative
            for c in 0..self.w.cols {
                grad_in[c] += delta * self.w.at(r, c)
                self.w.set(r, c, self.w.at(r, c) - lr * delta * self.input[c])
            }
            self.b[r] -= lr * delta
        }
        return grad_in
    }
}

struct Network { layers: List(Layer) }

impl Network {
    fn new(sizes: []usize) -> Network {
        var layers = List(Layer).new()
        for i in 1..sizes.len { layers.append(Layer.new(sizes[i - 1], sizes[i])) }
        return Network{ .layers = layers }
    }

    fn predict(self: *mut Self, x: []f64) -> List(f64) {
        var cur = x.to_owned()
        for i in 0..self.layers.len { cur = self.layers[i].forward(cur[..]) }
        return cur
    }

    /// One example: forward, squared-error loss, backward through every layer.
    fn train_one(self: *mut Self, x: []f64, target: []f64, lr: f64) -> f64 {
        let out = self.predict(x)
        var loss = 0.0
        var grad = List(f64).new()
        for o, i in out {
            let d = o - target[i]
            loss += d * d
            grad.append(2.0 * d)
        }
        var i = self.layers.len
        while i > 0 {
            i -= 1
            grad = self.layers[i].backward(grad[..], lr)
        }
        return loss
    }
}

/// A labelled data set: inputs and the expected outputs, row by row.
struct Data { xs: List(List(f64)), ys: List(List(f64)) }

fn xor_data() -> Data {
    var d = Data{ .xs = List(List(f64)).new(), .ys = List(List(f64)).new() }
    for a in 0..2 {
        for b in 0..2 {
            d.xs.append(List(f64).from([a as f64, b as f64]))
            d.ys.append(List(f64).from([if a != b { 1.0 } else { 0.0 }]))
        }
    }
    return d
}

/// Points on an inner ring (class 0) and an outer ring (class 1), with noise.
fn rings_data(n: usize) -> Data {
    var d = Data{ .xs = List(List(f64)).new(), .ys = List(List(f64)).new() }
    for i in 0..n {
        let outer = i % 2 == 1
        let radius = (if outer { 2.0 } else { 0.7 }) + (random.float() - 0.5) * 0.4
        let angle = random.float() * math.TAU
        d.xs.append(List(f64).from([radius * math.cos(angle) / 2.5, radius * math.sin(angle) / 2.5]))
        d.ys.append(List(f64).from([if outer { 1.0 } else { 0.0 }]))
    }
    return d
}

fn train(net: *mut Network, d: *Data, epochs: usize, lr: f64, report_every: usize) {
    for epoch in 0..epochs {
        var loss = 0.0
        for i in 0..d.xs.len { loss += net.train_one(d.xs[i][..], d.ys[i][..], lr) }
        if epoch % report_every == 0 or epoch + 1 == epochs {
            println("  epoch {>4}  loss {:.3}", .{epoch, loss / d.xs.len as f64})
        }
    }
}

fn accuracy(net: *mut Network, d: *Data) -> f64 {
    var right: usize = 0
    for i in 0..d.xs.len {
        let out = net.predict(d.xs[i][..])
        let guess = if out[0] >= 0.5 { 1.0 } else { 0.0 }
        if guess == d.ys[i][0] { right += 1 }
    }
    return 100.0 * right as f64 / d.xs.len as f64
}

fn main() {
    random.seed(7)

    println("XOR", .{})
    let xor = xor_data()
    var net = Network.new([2, 4, 1][..])
    train(&mut net, &xor, 3000, 0.5, 1000)
    for i in 0..xor.xs.len {
        let out = net.predict(xor.xs[i][..])
        println("  {} xor {} = {:.2}", .{xor.xs[i][0], xor.xs[i][1], out[0]})
    }
    println("  accuracy {:.0}%", .{accuracy(&mut net, &xor)})

    println("two rings", .{})
    let training = rings_data(200)
    let held_out = rings_data(100)
    var rings = Network.new([2, 8, 1][..])
    train(&mut rings, &training, 300, 0.3, 100)
    println("  training accuracy {:.0}%, held-out accuracy {:.0}%", .{accuracy(&mut rings, &training), accuracy(&mut rings, &held_out)})
}
topo/code/nn.expected
XOR
  epoch    0  loss 0.288
  epoch 1000  loss 0.002
  epoch 2000  loss 0.001
  epoch 2999  loss 0.000
  0.0 xor 0.0 = 0.02
  0.0 xor 1.0 = 0.98
  1.0 xor 0.0 = 0.98
  1.0 xor 1.0 = 0.02
  accuracy 100%
two rings
  epoch    0  loss 0.281
  epoch  100  loss 0.001
  epoch  200  loss 0.001
  epoch  299  loss 0.000
  training accuracy 100%, held-out accuracy 100%

It trains in well under a second in debug mode, and nx run topo/code/nn.nx --mode fast is faster still; random.seed(7) makes the run reproducible, so the book can show its numbers.

The pieces#

A matrix is a struct with a List(f64) and two methods, at and set, that turn row and column into an index. Matrix.random fills it with small values in [-scale, scale]; that is the whole of the linear algebra, because a network this size needs no more.

A layer holds its weights and biases and, after each forward pass, what it saw and what it produced. forward computes sigmoid(W x + b) one neuron at a time. backward receives the gradient of the loss with respect to its output, folds in the sigmoid's derivative (y * (1 - y)), updates every weight by lr * delta * input, and returns the gradient with respect to its input for the layer below. That function is backpropagation; there is nothing else to it.

A network is a list of layers. predict runs them forward; train_one runs one example forward, measures the squared error, and runs the layers backward in reverse order, handing each the gradient the one above computed. The loop in train does that for every example, some hundreds of times.

Where the language shows#

    fn forward(self: *mut Self, x: []f64) -> List(f64) {
        self.input = x.to_owned()

Ownership decides what gets copied. x is a slice, a view of the caller's data; the layer needs to keep the input for backward, so it takes a copy with to_owned, and that is the only copy in the forward pass. predict threads a List(f64) through the layers, moving it into each forward call (cur = self.layers[i].forward(cur[..])); nothing is retained, nothing is garbage collected, and nx leaks topo/code/nn.nx reports that all seven hundred thousand allocations the run makes were released.

        var i = self.layers.len
        while i > 0 {
            i -= 1
            grad = self.layers[i].backward(grad[..], lr)
        }

Indexing self.layers[i] on a *mut Self and calling a *mut Self method on the element mutates the layer in place. The checker knows i is in range here from the loop's shape; the arithmetic on f64 cannot panic; the compiled loop is what you would have written in C.

fn xor_data() -> Data {
    ...
            d.xs.append(List(f64).from([a as f64, b as f64]))

A Data holds two List(List(f64)), rows of inputs and rows of targets; List(f64).from([...]) builds a row from an array literal. Nested owning containers own their contents all the way down and are freed all the way down, with no code written for it.

What to try#

Next: the tools.