Networks

This module states the mathematics of neural networks in the same language as the rest of the Prelude, so a network can be read, differentiated and reasoned about beside the statistics it rests on. A dense layer is one contraction, backpropagation is the same contraction with its subscripts read the other way, attention is two contractions and a softmax, and a recurrent network is a definition that names itself — a loop, not an unrolled graph.

It is not a training framework. Control flow runs on the host, so a fit of a thousand steps is a thousand host turns, each handing its contractions to the executor. That is correct and useful for small models, teaching and verification, but it will not match a framework that fuses a whole training step into one kernel launch. For a production network, load a trained model instead of fitting it here.

Conventions. A batch is a matrix with one row per sample: inputs x are n×p, targets y are n×q. A weight matrix maps columns to columns, so w1 is p×h and w2 is h×q, and the widths of the network are read from the shapes of the weights you pass. No layer carries a bias: to have one, append a column of ones to x, as least squares does elsewhere in the Prelude. Initial weights and random draws are always arguments — nothing here draws a random number — so every result is a deterministic function of its inputs. The step functions stepFirst and stepSecond are the ones a model … : Learner declaration calls from its train method, one gradient step per call; the mlp…, auto… and gan… definitions run a whole fit of t steps inside one call.

Throughout, \(\operatorname{relu}(z) = \max(z, 0)\), \(\sigma(z) = 1/(1 + e^{-z})\) and \(\tanh\) are applied entrywise, \(\odot\) is the entrywise product, and the row-wise softmax of a matrix \(M \in \mathbb{R}^{n \times q}\) is

\[\operatorname{softmaxRows}(M)_{ij} = \frac{e^{M_{ij}}}{\sum_{j'=1}^{q} e^{M_{ij'}}} ,\]

computed after subtracting each row’s maximum from that row, which leaves the value unchanged and avoids overflow.

open Prelude.Networks

A dense layer, forward

A dense layer multiplies a batch by a weight matrix, \(A W\). The two-layer network places a rectifier between two such layers:

\[\hat Y = \operatorname{relu}(X W_1)\, W_2\]

dense

dense(a, w) · a activations n×p · w weights p×q → activations n×q

Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) be a batch of activations, one row per sample, and let \(W \in \mathbb{R}^{p \times q}\) (the argument w) be a weight matrix. Then

\[\operatorname{dense}(A, W) = A W, \qquad (A W)_{ij} = \sum_{k=1}^{p} A_{ik} W_{kj} .\]

One linear layer without a bias, one output row per input row.

Intuition Each output column is a weighted sum of the input columns, and row \(k\) of \(W\) says how input feature \(k\) feeds each output. Stacking two of these with nothing in between gives just another matrix, \(A W_1 W_2\), which is why hiddenOf puts a rectifier between the layers.

Use it when you need a linear map on a batch: a linear regression head, a change of width, or the query, key and value projections that selfAttention makes. There is no bias, so append a column of ones to a when you need an offset. The cost grows as \(n \cdot p \cdot q\).

regressiondimensionality-reduction

def dense(a, w)      ≔ einsum("np,pq->nq", a, w)

hiddenOf

hiddenOf(x, w1) · x inputs n×p · w1 first-layer weights p×h → hidden activations n×h

Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the inputs and \(W_1 \in \mathbb{R}^{p \times h}\) (the argument w1) the first-layer weights. Then

\[\operatorname{hiddenOf}(X, W_1) = \operatorname{relu}(X W_1) .\]

The rectified hidden layer of the two-layer network.

Intuition Each hidden unit is a hinge: zero on one side of a hyperplane in input space and rising linearly on the other. The output layer adds those hinges up, and with enough of them the sum can bend to follow almost any curve. That is where the network’s non-linearity comes from.

Use it when you want to inspect the features a trained network learned, feed them to another model, or find dead units: a column that is zero for every row gets no gradient and never recovers. Starting w1 at zero makes every unit dead, so draw small random starting weights.

nonlinear-modellingfeature-extractiontraining-diagnostics

def hiddenOf(x, w1)  ≔ relu(preHidden(x, w1))

outputs

outputs(x, w1, w2) · x inputs n×p · w1 first-layer weights p×h · w2 second-layer weights h×q → outputs n×q

Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the inputs, \(W_1 \in \mathbb{R}^{p \times h}\) (the argument w1) and \(W_2 \in \mathbb{R}^{h \times q}\) (the argument w2) the weights. Then

\[\operatorname{outputs}(X, W_1, W_2) = \operatorname{relu}(X W_1)\, W_2 .\]

The forward pass of the two-layer network. The output is raw (linear); apply softmaxRows to it for class shares.

Intuition The whole network in one line: hinge features, then a linear read-out. The output is left linear so it can take any sign and scale, which is what a regression target needs. The result is piecewise linear in the input, so beyond the data it continues in a straight line.

Use it when you already have weights, from mlpFirst and mlpSecond or from elsewhere, and want predictions without retraining. That is much cheaper than calling mlpPredict repeatedly. For class shares apply softmaxRows. Nothing here trains: it is the forward pass only.

regressionnonlinear-modellingmodel-inference

def outputs(x, w1, w2) ≔ dense(hiddenOf(x, w1), w2)

Building blocks

preHidden(x, w1) — Let \(X \in \mathbb{R}^{n \times p}\) and \(W_1 \in \mathbb{R}^{p \times h}\). Then \(\operatorname{preHidden}(X, W_1) = X W_1\). The hidden layer before the rectifier.

Backward, and one gradient step

The network is trained on the squared error \(L = \frac{1}{2n}\lVert \hat Y - Y \rVert^2\). Its gradients are the forward contractions with the subscripts rearranged:

\[\delta_{\text{out}} = \frac{\hat Y - Y}{n}, \qquad \delta_{\text{hid}} = (\delta_{\text{out}} W_2^{\top}) \odot [X W_1 > 0]\] \[\nabla_{W_2} L = H^{\top} \delta_{\text{out}}, \qquad \nabla_{W_1} L = X^{\top} \delta_{\text{hid}}\]

where \(H\) is the hidden activation. A descent step is \(W \leftarrow W - r \nabla_W L\). The rectifier’s derivative is taken as \([z > 0]\), so it is \(0\) at \(z = 0\).

stepFirst

stepFirst(r, x, y, w1, w2) · r learning rate scalar · x inputs n×p · y targets n×q · w1 weights p×h · w2 weights h×q → updated first-layer weights p×h

Let \(r > 0\) (the argument r) be the learning rate, \(X \in \mathbb{R}^{n \times p}\) (the argument x) the inputs, \(Y \in \mathbb{R}^{n \times q}\) (the argument y) the targets, and \(W_1 \in \mathbb{R}^{p \times h}\), \(W_2 \in \mathbb{R}^{h \times q}\) (the arguments w1, w2) the current weights. Then

\[\operatorname{stepFirst}(r, X, Y, W_1, W_2) = W_1 - r\, X^{\top} \delta_{\text{hid}} ,\]

where \(\delta_{\text{out}} = \frac{1}{n}\bigl(\operatorname{relu}(X W_1) W_2 - Y\bigr)\) and \(\delta_{\text{hid}} = (\delta_{\text{out}} W_2^{\top}) \odot [X W_1 > 0]\), the bracket applied entrywise.

One gradient-descent step on w1, with the gradient evaluated at the current w1 and w2.

Intuition The output error travels back through w2 to each hidden unit, but only through units that were active: the rectifier’s derivative \([XW_1 > 0]\) closes the path for the others. Each first-layer weight then moves in proportion to its input times the error that reached its unit.

Use it when you want control over training that mlpFirst does not give: watching the loss, stopping early, changing the rate, or mini-batches (pass a different x, y batch on each call). It is the step a model … : Learner calls. Compute it and stepSecond from the same old weights before replacing either. A rate that is too large makes the loss grow.

gradient-based-trainingtraining-control

def stepFirst(r, x, y, w1, w2)  ≔ descend(w1, gradFirst(x, y, w1, w2), r)

stepSecond

stepSecond(r, x, y, w1, w2) · r learning rate scalar · x inputs n×p · y targets n×q · w1 weights p×h · w2 weights h×q → updated second-layer weights h×q

Let \(r > 0\), \(X \in \mathbb{R}^{n \times p}\), \(Y \in \mathbb{R}^{n \times q}\), \(W_1 \in \mathbb{R}^{p \times h}\) and \(W_2 \in \mathbb{R}^{h \times q}\) be as in stepFirst. Then

\[\operatorname{stepSecond}(r, X, Y, W_1, W_2) = W_2 - r\, \operatorname{relu}(X W_1)^{\top} \delta_{\text{out}} ,\]

where \(\delta_{\text{out}} = \frac{1}{n}\bigl(\operatorname{relu}(X W_1) W_2 - Y\bigr)\).

One gradient-descent step on w2, evaluated at the same current weights as stepFirst. Calling both with the weights a step was entered with makes one simultaneous step on both layers; this is how a model … : Learner declaration’s train method uses them.

Intuition Seen from the output layer, the hidden activations are just input features, so this is one gradient step of linear least squares on them. A weight grows when its hidden unit was active on rows where the output was too low.

Use it when you drive training yourself alongside stepFirst. Evaluate both at the same (w1, w2) and only then replace them, which is a simultaneous step. Passing the already-updated w1 into this call gives a different, sequential method. For a fixed number of full-batch steps, mlpFirst and mlpSecond run the loop for you.

gradient-based-trainingtraining-control

def stepSecond(r, x, y, w1, w2) ≔ descend(w2, gradSecond(x, y, w1, w2), r)

Building blocks

outDelta(x, y, w1, w2) — Let \(X \in \mathbb{R}^{n \times p}\), \(Y \in \mathbb{R}^{n \times q}\), \(W_1 \in \mathbb{R}^{p \times h}\) and \(W_2 \in \mathbb{R}^{h \times q}\). Then \(\operatorname{outDelta}(X, Y, W_1, W_2) = \frac{1}{n}\bigl(\operatorname{relu}(X W_1) W_2 - Y\bigr)\). The output error divided by the number of rows.
hidDelta(x, y, w1, w2) — Let \(X, Y, W_1, W_2\) be as in outDelta and \(\delta_{\text{out}} = \operatorname{outDelta}(X, Y, W_1, W_2)\). Then \(\operatorname{hidDelta}(X, Y, W_1, W_2) = (\delta_{\text{out}} W_2^{\top}) \odot [X W_1 > 0]\). The output error carried back through w2 and gated by the rectifier’s derivative.
gradSecond(x, y, w1, w2) — Let \(X, Y, W_1, W_2\) be as in outDelta. Then \(\operatorname{gradSecond}(X, Y, W_1, W_2) = \operatorname{relu}(X W_1)^{\top} \operatorname{outDelta}(X, Y, W_1, W_2) \in \mathbb{R}^{h \times q}\). The gradient of the loss with respect to w2.
gradFirst(x, y, w1, w2) — Let \(X, Y, W_1, W_2\) be as in outDelta. Then \(\operatorname{gradFirst}(X, Y, W_1, W_2) = X^{\top} \operatorname{hidDelta}(X, Y, W_1, W_2) \in \mathbb{R}^{p \times h}\). The gradient of the loss with respect to w1.

The fit

A fit runs t steps of full-batch gradient descent, updating both weight matrices together at every step. A loop answers one tensor, so the fit is walked twice — once to answer each weight matrix — and both walks step both matrices, so the two answers come from the same run. Each call therefore retrains from the starting weights; to watch or continue a fit, call stepFirst and stepSecond yourself.

mlpFirst

mlpFirst(t, r, x, y, w1, w2) · t number of steps scalar · r learning rate scalar · x inputs n×p · y targets n×q · w1 starting weights p×h · w2 starting weights h×q → trained first-layer weights p×h

Let \(t \in \mathbb{N}\) (the argument t) be the number of steps, \(r > 0\) (the argument r) the learning rate, \(X \in \mathbb{R}^{n \times p}\) and \(Y \in \mathbb{R}^{n \times q}\) (the arguments x, y) the data, and \(W_1 \in \mathbb{R}^{p \times h}\), \(W_2 \in \mathbb{R}^{h \times q}\) (the arguments w1, w2) the starting weights. Let \(W_1^{(0)} = W_1\), \(W_2^{(0)} = W_2\) and, for \(s = 0, 1, \dots, t-1\),

\[W_1^{(s+1)} = \operatorname{stepFirst}\bigl(r, X, Y, W_1^{(s)}, W_2^{(s)}\bigr), \qquad W_2^{(s+1)} = \operatorname{stepSecond}\bigl(r, X, Y, W_1^{(s)}, W_2^{(s)}\bigr).\]

Then

\[\operatorname{mlpFirst}(t, r, X, Y, W_1, W_2) = W_1^{(t)} .\]

The first-layer weights after t gradient steps.

Intuition After training, each column of w1 is a feature detector: the direction in input space that switches one hidden unit on. The loop steps both layers together but can answer only one tensor, so this call returns the first layer of that joint run.

Use it when you want to keep a trained network and reuse it through outputs or hiddenOf rather than retraining for every query. Every call retrains from the starting weights, so it and mlpSecond run two identical fits. They agree because the fits are deterministic, provided the arguments are identical. Draw small random starting weights, because an all-zero start never moves.

gradient-based-trainingnonlinear-modellingfeature-extraction

def mlpFirst(t, r, x, y, w1, w2)  ≔ firstLoop(t, r, x, y, w1, w2)

mlpSecond

mlpSecond(t, r, x, y, w1, w2) · arguments as mlpFirst → trained second-layer weights h×q

Let \(t, r, X, Y, W_1, W_2\) be as in mlpFirst, and let \(W_2^{(t)}\) be the iterate defined there. Then

\[\operatorname{mlpSecond}(t, r, X, Y, W_1, W_2) = W_2^{(t)} .\]

The second-layer weights after t gradient steps.

Intuition The trained read-out: how strongly each hidden feature pushes each output up or down. Together with mlpFirst it is the complete model.

Use it when you need the second half of the trained weights to predict with outputs or to inspect which hidden units drive which target. Pass exactly the arguments you gave mlpFirst, because a different t, rate or starting point is a different run. It performs its own full fit of t steps.

gradient-based-trainingnonlinear-modelling

def mlpSecond(t, r, x, y, w1, w2) ≔ secondLoop(t, r, x, y, w1, w2)

mlpPredict

mlpPredict(t, r, x, y, w1, w2, z) · arguments as mlpFirst · z queries m×p → outputs m×q

Let \(t, r, X, Y, W_1, W_2\) and the iterates \(W_1^{(t)}, W_2^{(t)}\) be as in mlpFirst, and let \(Z \in \mathbb{R}^{m \times p}\) (the argument z) be the queries. Then

\[\operatorname{mlpPredict}(t, r, X, Y, W_1, W_2, Z) = \operatorname{relu}\bigl(Z W_1^{(t)}\bigr)\, W_2^{(t)} .\]

Trains the network for t steps on x and y, then runs the forward pass on the queries z.

Intuition Fit, then ask, in one call. A single hidden rectifier layer can approximate any continuous function if it is wide enough, so this is a general-purpose non-linear regressor. Between training points it interpolates piecewise linearly, and outside them it extends in straight lines.

Use it when you have a small regression problem, are teaching or verifying, or need a one-shot answer. Each call runs the full fit twice (once per weight matrix), so for repeated queries get the weights once with mlpFirst and mlpSecond and call outputs. Put all queries in one z.

regressionnonlinear-modelling

def mlpPredict(t, r, x, y, w1, w2, z) ≔ outputs(z, mlpFirst(t, r, x, y, w1, w2), mlpSecond(t, r, x, y, w1, w2))

mlpFit

mlpFit(t, r, x, y, w1, w2) · arguments as mlpFirst → fitted outputs n×q

Let \(t, r, X, Y, W_1, W_2\) and the iterates \(W_1^{(t)}, W_2^{(t)}\) be as in mlpFirst. Then

\[\operatorname{mlpFit}(t, r, X, Y, W_1, W_2) = \operatorname{relu}\bigl(X W_1^{(t)}\bigr)\, W_2^{(t)} .\]

The trained network’s outputs on its own training inputs.

Intuition The model’s answers on the data it learned from. Comparing them with y shows how well the optimiser did, not how well the model generalises.

Use it when you are checking that training converged. Try a few values of t and rate, and look for a falling training error. A residual that stays large means the network is too narrow or the rate too small. A tiny training error says nothing about new data, so judge that with mlpPredict on held-out rows.

training-diagnosticsregression

def mlpFit(t, r, x, y, w1, w2) ≔ mlpPredict(t, r, x, y, w1, w2, x)

mlpShares

mlpShares(t, r, x, y, w1, w2, z) · arguments as mlpPredict → class shares m×q, each row summing to one

Let \(t, r, X, Y, W_1, W_2\) and \(Z \in \mathbb{R}^{m \times p}\) be as in mlpPredict. Then

\[\operatorname{mlpShares}(t, r, X, Y, W_1, W_2, Z) = \operatorname{softmaxRows}\bigl(\operatorname{mlpPredict}(t, r, X, Y, W_1, W_2, Z)\bigr) .\]

The raw outputs are treated as logits; the training itself is the same squared-error descent as mlpPredict.

Intuition Softmax turns a row of scores into positive shares that sum to one, with larger scores getting exponentially more. The network is still trained on squared error, so with one-hot targets its outputs sit near 0 and 1, and softmax of \([1, 0]\) is about \([0.73, 0.27]\), not \([1, 0]\).

Use it when you want a quick multi-class classifier from one-hot targets: the largest share is the predicted class. Do not read the shares as calibrated probabilities, since they are systematically soft. If confidence values matter, calibrate them on held-out data. It retrains on every call, like mlpPredict.

classificationprobability-calibration

def mlpShares(t, r, x, y, w1, w2, z) ≔ softmaxRows(mlpPredict(t, r, x, y, w1, w2, z))

Building blocks

firstLoop(t, r, x, y, w1, w2) — Let \(t \in \mathbb{N}\), \(r > 0\), \(X \in \mathbb{R}^{n \times p}\), \(Y \in \mathbb{R}^{n \times q}\), \(W_1 \in \mathbb{R}^{p \times h}\) and \(W_2 \in \mathbb{R}^{h \times q}\). Then \(\operatorname{firstLoop}(0, r, X, Y, W_1, W_2) = W_1\) and \(\operatorname{firstLoop}(t, r, X, Y, W_1, W_2) = \operatorname{firstLoop}\bigl(t - 1, r, X, Y, \operatorname{stepFirst}(r, X, Y, W_1, W_2), \operatorname{stepSecond}(r, X, Y, W_1, W_2)\bigr)\). Steps both weight matrices t times and answers the first.
secondLoop(t, r, x, y, w1, w2) — Let \(t, r, X, Y, W_1, W_2\) be as in firstLoop. Then \(\operatorname{secondLoop}(0, r, X, Y, W_1, W_2) = W_2\) and \(\operatorname{secondLoop}(t, r, X, Y, W_1, W_2) = \operatorname{secondLoop}\bigl(t - 1, r, X, Y, \operatorname{stepFirst}(r, X, Y, W_1, W_2), \operatorname{stepSecond}(r, X, Y, W_1, W_2)\bigr)\). Steps both weight matrices t times and answers the second.

The autoencoder

An autoencoder is a network trained to reproduce its own input through a narrow middle layer. This one is linear: \(\hat X = X W_1 W_2\) with w1 of shape p×k and w2 of shape k×p, trained on \(\frac{1}{2n}\lVert \hat X - X \rVert^2\). A linear autoencoder with a bottleneck of k units spans the same subspace as the first k principal components — so it reaches the same reconstruction error as PCA — although its units are neither orthogonal nor ordered.

encoded

encoded(x, w1) · x data n×p · w1 encoder weights p×k → codes n×k

Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the data and \(W_1 \in \mathbb{R}^{p \times k}\) (the argument w1) the encoder weights. Then

\[\operatorname{encoded}(X, W_1) = X W_1 .\]

The bottleneck representation.

Intuition Each row is squeezed to k numbers, its coordinates in a k-dimensional subspace. Training forces the subspace to keep whatever is needed to rebuild the row, so the codes are a compact summary of the data.

Use it when you need dimensionality reduction: plotting in 2 or 3 dimensions, compressing features before another model, or comparing rows in a smaller space. Use trained encoder weights from autoCode. The units are not orthogonal or ordered by importance, so if you need ranked, uncorrelated components use PCA (see Multivariate).

dimensionality-reductionfeature-extractiondata-compression

def encoded(x, w1)          ≔ dense(x, w1)

rebuilt

rebuilt(x, w1, w2) · x data n×p · w1 encoder weights p×k · w2 decoder weights k×p → reconstruction n×p

Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(W_1 \in \mathbb{R}^{p \times k}\) (the argument w1) and \(W_2 \in \mathbb{R}^{k \times p}\) (the argument w2). Then

\[\operatorname{rebuilt}(X, W_1, W_2) = X W_1 W_2 .\]

The reconstruction of the data through the bottleneck.

Intuition The row as the k-number summary remembers it. Anything the bottleneck cannot express, such as noise or unusual structure, is dropped, so the gap between a row and its rebuild measures how typical the row is.

Use it when you want a per-row reconstruction error as an anomaly score (rows unlike the training bulk rebuild poorly), a denoised copy of the data, or a check of whether k is large enough. Use the trained weights from autoCode and autoBuild.

anomaly-detectiondenoisingdimensionality-reduction

def rebuilt(x, w1, w2)      ≔ dense(encoded(x, w1), w2)

autoCode

autoCode(t, r, x, w1, w2) · t number of steps scalar · r learning rate scalar · x data n×p · w1 starting encoder p×k · w2 starting decoder k×p → trained encoder p×k

Let \(t \in \mathbb{N}\) (the argument t), \(r > 0\) (the argument r), \(X \in \mathbb{R}^{n \times p}\) (the argument x), and let \(W_1 \in \mathbb{R}^{p \times k}\), \(W_2 \in \mathbb{R}^{k \times p}\) (the arguments w1, w2) be the starting weights. Let \(W_1^{(0)} = W_1\), \(W_2^{(0)} = W_2\) and, for \(s = 0, 1, \dots, t-1\),

\[W_1^{(s+1)} = W_1^{(s)} - r\, X^{\top} \Delta_s W_2^{(s)\top}, \qquad W_2^{(s+1)} = W_2^{(s)} - r\, \bigl(X W_1^{(s)}\bigr)^{\top} \Delta_s ,\]

where \(\Delta_s = \frac{1}{n}\bigl(X W_1^{(s)} W_2^{(s)} - X\bigr)\). Then

\[\operatorname{autoCode}(t, r, X, W_1, W_2) = W_1^{(t)} .\]

The encoder weights after t simultaneous gradient steps on both matrices.

Intuition Gradient descent on reconstruction error, with no eigendecomposition, settles on the same subspace as the top k principal components. The encoder it returns is one basis for that subspace, not necessarily PCA’s basis.

Use it when you want an encoder to apply to new rows with encoded. There is no bias, so centre x first (subtract column means) if you want PCA’s subspace, because uncentred data pulls the first direction toward the mean. Draw small non-zero starting weights. Each call runs a full fit of t steps from the start.

dimensionality-reductiongradient-based-training

def autoCode(t, r, x, w1, w2)   ≔ codeLoop(t, r, x, w1, w2)

autoBuild

autoBuild(t, r, x, w1, w2) · arguments as autoCode → trained decoder k×p

Let \(t, r, X, W_1, W_2\) and the iterate \(W_2^{(t)}\) be as in autoCode. Then

\[\operatorname{autoBuild}(t, r, X, W_1, W_2) = W_2^{(t)} .\]

The decoder weights after the same t steps.

Intuition The decoder maps k codes back to p columns. Each of its k rows is a pattern over the original features, and every rebuilt row is a weighted mix of those patterns.

Use it when you want to turn codes back into data, for example after averaging or editing codes, or to look at the patterns the bottleneck kept. Pass the same arguments as autoCode so that both halves come from one run. It performs its own full fit.

dimensionality-reductiondata-compression

def autoBuild(t, r, x, w1, w2)  ≔ buildLoop(t, r, x, w1, w2)

autoFit

autoFit(t, r, x, w1, w2) · arguments as autoCode → reconstruction n×p

Let \(t, r, X, W_1, W_2\) and the iterates \(W_1^{(t)}, W_2^{(t)}\) be as in autoCode. Then

\[\operatorname{autoFit}(t, r, X, W_1, W_2) = X W_1^{(t)} W_2^{(t)} .\]

Trains the autoencoder and returns its reconstruction of x.

Intuition Train, then push the training data through the bottleneck and back. The result is the best rank-k linear approximation of x that the descent reached.

Use it when you want to choose k: fit several values and watch reconstruction error level off. It also gives anomaly scores or a denoised copy of the training data in one call. It trains twice (once for the encoder, once for the decoder), so for new rows get the weights once and use rebuilt.

dimensionality-reductionanomaly-detectiondenoising

def autoFit(t, r, x, w1, w2)    ≔ rebuilt(x, autoCode(t, r, x, w1, w2), autoBuild(t, r, x, w1, w2))

Building blocks

rebuildDelta(x, w1, w2) — Let \(X \in \mathbb{R}^{n \times p}\), \(W_1 \in \mathbb{R}^{p \times k}\) and \(W_2 \in \mathbb{R}^{k \times p}\). Then \(\operatorname{rebuildDelta}(X, W_1, W_2) = \frac{1}{n}(X W_1 W_2 - X)\). The reconstruction error divided by the number of rows.
rebuildSecond(x, w1, w2) — Let \(X, W_1, W_2\) be as in rebuildDelta. Then \(\operatorname{rebuildSecond}(X, W_1, W_2) = (X W_1)^{\top} \operatorname{rebuildDelta}(X, W_1, W_2) \in \mathbb{R}^{k \times p}\). The gradient of the reconstruction loss with respect to the decoder.
rebuildFirst(x, w1, w2) — Let \(X, W_1, W_2\) be as in rebuildDelta. Then \(\operatorname{rebuildFirst}(X, W_1, W_2) = X^{\top} \operatorname{rebuildDelta}(X, W_1, W_2)\, W_2^{\top} \in \mathbb{R}^{p \times k}\). The gradient of the reconstruction loss with respect to the encoder.
codeLoop(t, r, x, w1, w2) — Let \(t \in \mathbb{N}\), \(r > 0\) and \(X, W_1, W_2\) be as in rebuildDelta. Then \(\operatorname{codeLoop}(0, r, X, W_1, W_2) = W_1\) and \(\operatorname{codeLoop}(t, r, X, W_1, W_2) = \operatorname{codeLoop}\bigl(t - 1, r, X, W_1 - r \operatorname{rebuildFirst}(X, W_1, W_2), W_2 - r \operatorname{rebuildSecond}(X, W_1, W_2)\bigr)\). Steps both matrices t times and answers the encoder.
buildLoop(t, r, x, w1, w2) — Let \(t, r, X, W_1, W_2\) be as in codeLoop. Then \(\operatorname{buildLoop}(0, r, X, W_1, W_2) = W_2\) and \(\operatorname{buildLoop}(t, r, X, W_1, W_2) = \operatorname{buildLoop}\bigl(t - 1, r, X, W_1 - r \operatorname{rebuildFirst}(X, W_1, W_2), W_2 - r \operatorname{rebuildSecond}(X, W_1, W_2)\bigr)\). Steps both matrices t times and answers the decoder.

Layer normalisation

Layer normalisation standardises each row across its features, not each feature across the batch — which is why it can be applied to one sequence at a time. For row i with p columns:

\[\operatorname{LN}(a)_{ij} = \frac{a_{ij} - \mu_i}{\sqrt{\sigma_i^2 + 10^{-5}}}, \qquad \mu_i = \frac{1}{p}\sum_j a_{ij}, \quad \sigma_i^2 = \frac{1}{p}\sum_j (a_{ij} - \mu_i)^2\]

layerNorm

layerNorm(a) · a activations n×p → normalised activations n×p

Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) be a batch of activations. Then

\[\operatorname{layerNorm}(A)_{ij} = \frac{A_{ij} - \mu_i}{\sqrt{\sigma_i^2 + 10^{-5}}}, \qquad i = 1, \dots, n, \quad j = 1, \dots, p ,\]

where \(\mu_i = \frac{1}{p}\sum_{j=1}^{p} A_{ij}\) and \(\sigma_i^2 = \frac{1}{p}\sum_{j=1}^{p} (A_{ij} - \mu_i)^2\).

Each row centred to mean zero and scaled to unit (population) variance, with \(10^{-5}\) added to the variance for stability. There is no learned gain or shift; multiply and add your own if you want them.

Intuition Whatever scale or offset a row’s activations drift to, the next layer sees them centred with unit spread. That keeps layers deep in a stack from being overwhelmed or starved as training changes the weights. Each row is normalised on its own, so the result does not depend on the batch.

Use it when you are stacking attention or dense layers and training becomes unstable, when batches are tiny or have size one, or when sequence lengths vary. Add your own learned gain and shift if needed. Rows need more than one column, since a single column normalises to 0. To standardise each feature across samples, this is the wrong axis.

training-stabilitynumerical-stability

def layerNorm(a) ≔ rowCentre(a) / lifted(a, sqrt(rowVar(a) + 1e-5))

Building blocks

rowMean(a) — Let \(A \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rowMean}(A)_i = \frac{1}{p}\sum_{j=1}^{p} A_{ij}\) for \(i = 1, \dots, n\). The mean of each row.
rowCentre(a) — Let \(A \in \mathbb{R}^{n \times p}\) and \(\mu = \operatorname{rowMean}(A)\). Then \(\operatorname{rowCentre}(A) = A - \mu \mathbf{1}_p^{\top}\), that is \(\operatorname{rowCentre}(A)_{ij} = A_{ij} - \mu_i\). Each row minus its own mean.
rowVar(a) — Let \(A \in \mathbb{R}^{n \times p}\) and \(C = \operatorname{rowCentre}(A)\). Then \(\operatorname{rowVar}(A)_i = \frac{1}{p}\sum_{j=1}^{p} C_{ij}^2\). The population variance of each row.

Attention

Scaled dot-product attention weighs each value by how well its key matches the query:

\[\operatorname{Attention}(Q, K, V) = \operatorname{softmax}_{\text{rows}}\!\left(\frac{Q K^{\top}}{\sqrt{d}}\right) V\]

The division by \(\sqrt{d}\) matters: without it the dot products grow with the embedding width, the softmax saturates and its gradient vanishes. A causal mask lets position i attend only to positions \(j \le i\); the forbidden scores are lowered by \(10^{300}\) rather than removed, so each row is still a proper softmax.

attentionOf

attentionOf(q, k, v, d) · q queries m×e · k keys n×e · v values n×c · d scale, usually the key width scalar → attended values m×c

Let \(Q \in \mathbb{R}^{m \times e}\) (the argument q) be the queries, \(K \in \mathbb{R}^{n \times e}\) (the argument k) the keys, \(V \in \mathbb{R}^{n \times c}\) (the argument v) the values, and let \(d > 0\) (the argument d). Then

\[\operatorname{attentionOf}(Q, K, V, d) = \operatorname{softmaxRows}\!\left(\frac{Q K^{\top}}{\sqrt{d}}\right) V .\]

Each output row is a weighted average of the rows of v, with weights \(\operatorname{softmax}(q_i K^{\top} / \sqrt{d})\).

Intuition A soft dictionary lookup. Each query is compared with every key, softmax turns the similarities into weights, and the answer is the weighted average of the values. Dividing by \(\sqrt{d}\) keeps the scores from growing with width, so the softmax does not collapse onto one key and stop passing gradient.

Use it when outputs should draw on whichever inputs are relevant: cross-attention between two sequences, retrieval over a memory, or pooling a variable-size set into fixed outputs. Set d to the key width. It ignores order, so add position information to the inputs yourself. The score matrix is m×n.

attentionsequence-modellinginformation-retrieval

def attentionOf(q, k, v, d)  ≔ einsum("qk,ke->qe", softmaxRows(attentionScores(q, k, d)), v)

causalAttention

causalAttention(q, k, v, d) · q queries n×e · k keys n×e · v values n×c · d scale scalar → attended values n×c

Let \(Q, K \in \mathbb{R}^{n \times e}\) (the arguments q, k) be the queries and keys, \(V \in \mathbb{R}^{n \times c}\) (the argument v) the values, and let \(d > 0\) (the argument d). Then

\[\operatorname{causalAttention}(Q, K, V, d) = \operatorname{softmaxRows}\bigl(S - 10^{300} M\bigr)\, V ,\]

where \(S = Q K^{\top} / \sqrt{d}\) and \(M_{ij} = [\,j > i\,]\) marks the entries strictly above the diagonal.

Row i may only read positions \(0 \ldots i\) — a token never sees its own future. Queries and keys must have the same number of rows.

Intuition The same lookup with the future hidden. Because position i can see only positions up to i, a whole sequence can be trained to predict each next element in one pass without any position seeing its answer.

Use it when you are building autoregressive models: next-token language modelling, or forecasting where later values must not leak into earlier predictions. Rows of q and k must be the same positions in the same order. There is no causal version of selfAttention, so project the sequence with dense and pass the three projections here.

language-modellingattentiondata-leakage

def causalAttention(q, k, v, d) ≔ einsum("qk,ke->qe", softmaxRows(causalScores(q, k, d)), v)

selfAttention

selfAttention(a, wq, wk, wv, d) · a sequence n×p · wq query projection p×e · wk key projection p×e · wv value projection p×c · d scale scalar → attended sequence n×c

Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) be the sequence, \(W_q, W_k \in \mathbb{R}^{p \times e}\) and \(W_v \in \mathbb{R}^{p \times c}\) (the arguments wq, wk, wv) the projections, and let \(d > 0\) (the argument d). Then

\[\operatorname{selfAttention}(A, W_q, W_k, W_v, d) = \operatorname{attentionOf}(A W_q, A W_k, A W_v, d) = \operatorname{softmaxRows}\!\left(\frac{A W_q W_k^{\top} A^{\top}}{\sqrt{d}}\right) A W_v .\]

Attention whose queries, keys and values are three projections of the same sequence.

Intuition Every position gathers information from every other position of the same sequence in one step. The three projections separate roles: what a position is looking for (query), what it can be found by (key), and what it passes on when found (value).

Use it when each element’s meaning depends on context, such as words in a sentence or events in a log, and distant elements need to interact without a recurrence. This is one head with no output projection, residual or normalisation, so compose those from dense and layerNorm. It sees in both directions. For generation use causalAttention. Cost grows as \(n^2\).

attentionsequence-modellinglong-range-dependencies

def selfAttention(a, wq, wk, wv, d) ≔ attentionOf(dense(a, wq), dense(a, wk), dense(a, wv), d)

Building blocks

attentionScores(q, k, d) — Let \(Q \in \mathbb{R}^{m \times e}\), \(K \in \mathbb{R}^{n \times e}\) and \(d > 0\). Then \(\operatorname{attentionScores}(Q, K, d) = Q K^{\top} / \sqrt{d} \in \mathbb{R}^{m \times n}\). The score matrix, query-key dot products divided by the square root of d.
causalOf(s) — Let \(S \in \mathbb{R}^{n \times n}\). Then \(\operatorname{causalOf}(S)_{ij} = 1 - [\,j > i\,] = [\,j \le i\,]\); only the shape of \(S\) is read. A 0/1 mask the shape of s holding ones on and below the diagonal.
causalScores(q, k, d) — Let \(Q, K \in \mathbb{R}^{n \times e}\), \(d > 0\) and \(S = \operatorname{attentionScores}(Q, K, d)\). Then \(\operatorname{causalScores}(Q, K, d) = S - 10^{300}\bigl(1 - \operatorname{causalOf}(S)\bigr)\), that is \(S_{ij} - 10^{300} [\,j > i\,]\). The attention scores with every entry above the diagonal lowered by 1e300.

Dropout

Dropout zeroes each unit with probability p and scales the survivors by \(1/(1-p)\), so the layer’s expected output is unchanged and inference needs no separate code path. The random draws are an argument.

dropped

dropped(a, u, p) · a activations tensor · u uniform draws in [0, 1), same shape as a · p drop probability scalar → activations, same shape as a

Let \(A\) (the argument a) be a tensor of activations, let \(U\) (the argument u) be independent uniform draws in \([0, 1)\) of the same shape, and let \(0 \le p \le 1\) (the argument p). Then

\[\operatorname{dropped}(A, U, p) = \frac{A \odot [\,U \ge p\,]}{\max(1 - p,\ 10^{-12})} ,\]

with the bracket applied entrywise.

The entries where \(u \ge p\) are kept and scaled up; the rest are zeroed. Because the divisor is floored at \(10^{-12}\), p = 1 returns zeros rather than failing.

Intuition Silencing a random subset of units on each step stops the network from leaning on any single unit or fragile combination. It behaves like training many thinned networks that share weights. Scaling survivors by \(1/(1-p)\) keeps the expected output unchanged, so inference uses the activations as they are.

Use it when a network fits the training data much better than held-out data. Apply it during training only, and pass fresh draws u on every step, because reusing the same draws is just a fixed mask. Values of p between 0.1 and 0.5 are common. At inference, skip it or use p = 0.

regularisationoverfitting-control

def dropped(a, u, p) ≔ a · (u ≥ p) / max(1 − p, 1e-12)

Convolution and pooling

A one-dimensional convolution slides a kernel along each signal. Here it is the unflipped form (cross-correlation), as in most neural-network libraries:

\[y_{it} = \sum_{j=0}^{K-1} x_{i,\,t+j}\, k_j\]

The shift is built from a 0/1 window tensor that holds a one at \((t, j, s)\) exactly when \(s = t + j\), so the whole convolution is a single contraction. The output has the same length as the input: past the end of the signal no position matches, which pads the tail with zeros.

In this section positions along a signal are counted from zero, \(0, 1, \dots, L-1\), and kernel taps likewise, \(0, 1, \dots, K-1\); a signal is read as zero past its end, \(X_{i,s} = 0\) for \(s \ge L\).

convolve

convolve(x, k) · x signals b×L (one signal per row) · k kernel vector K → filtered signals b×L

Let \(X \in \mathbb{R}^{b \times L}\) (the argument x) hold one signal per row and let \(k \in \mathbb{R}^{K}\) (the argument k) be the kernel. Then, for \(i = 1, \dots, b\) and \(t = 0, 1, \dots, L-1\),

\[\operatorname{convolve}(X, k)_{it} = \sum_{s=0}^{L-1} \sum_{j=0}^{K-1} X_{is}\, [\,s = t + j\,]\, k_j = \sum_{j=0}^{K-1} X_{i,\,t+j}\, k_j .\]

The sliding weighted sum above, zero-padded at the tail.

Intuition One small pattern, the kernel, is slid along the signal, and the output at each position says how strongly the signal there matches it. The same weights are reused everywhere, so a feature is detected wherever it occurs. This is the unflipped (cross-correlation) form.

Use it when you need smoothing (equal weights \(1/K\) give a moving average), differencing or edge detection ([1, -1]), or matched filtering for a known shape. The window runs from \(t\) to \(t+K-1\), so on a time series it reads future values. The last \(K-1\) outputs see zero padding. The kernel is an argument, not learned, and memory grows as \(L^2 K\).

signal-processingsmoothingfeature-extraction

def convolve(x, k)   ≔ einsum("is,tjs,j->it", x, windowOf(x, k), k)

maxPool

maxPool(x, k) · x signals b×L · k any vector whose length K is the window width → pooled signals b×L

Let \(X \in \mathbb{R}^{b \times L}\) (the argument x) and let \(k \in \mathbb{R}^{K}\) (the argument k) be any vector of length \(K\). Then, for \(i = 1, \dots, b\) and \(t = 0, 1, \dots, L-1\),

\[\operatorname{maxPool}(X, k)_{it} = \max_{0 \le j \le K-1} X_{i,\,t+j} .\]

Only the length of k is used. Because the tail is padded with zeros, on a signal that goes negative the last windows read at least zero.

Intuition Keep the strongest response in each window and forget exactly where in the window it happened. A feature detected one or two steps earlier or later still gives the same pooled value, which makes later stages tolerant of small shifts.

Use it when you follow convolve and care whether a pattern appeared nearby, not its exact position, or when you need a running peak or envelope of a signal. The window moves one step at a time and the output keeps length L. To shrink the signal, keep every K-th column yourself. On signals that go negative, the last windows read at least 0 because of the padding.

signal-processingfeature-extractionshift-tolerance

def maxPool(x, k)    ≔ maximum(gathered(x, k), 2)

meanPool

meanPool(x, k) · x signals b×L · k window vector K → pooled signals b×L

Let \(X \in \mathbb{R}^{b \times L}\) (the argument x) and \(k \in \mathbb{R}^{K}\) (the argument k). Then, for \(i = 1, \dots, b\) and \(t = 0, 1, \dots, L-1\),

\[\operatorname{meanPool}(X, k)_{it} = \frac{1}{c} \sum_{j=0}^{K-1} X_{i,\,t+j} , \qquad c = \operatorname{columns}(k^0) ,\]

where \(k^0\) is the vector of ones shaped like \(k\).

Intended as the mean of each window, with \(c = K\). But a vector has no second axis, so columns(k^0) is undefined for it, and in this release a call with a vector window fails at run time with an axis error. Use convolve(x, k) with a kernel of equal weights \(1/K\) for a moving average instead.

Intuition Averaging each window keeps the local level of the signal and cancels noise that changes sign from sample to sample. It is the smooth counterpart of maxPool, which keeps only the peak.

Use it when you would want a moving average or window mean. In this release, do not call it: with a vector window it fails at run time. Call convolve(x, k) with a kernel of K equal weights \(1/K\), which computes the intended window mean (with the same tail padding).

signal-processingsmoothing

def meanPool(x, k)   ≔ einsum("itj->it", gathered(x, k)) / columns(k^0)

Building blocks

siteIds(x) — Let \(X \in \mathbb{R}^{b \times L}\). Then \(\operatorname{siteIds}(X) = (0, 1, \dots, L-1) \in \mathbb{R}^{L}\). The positions along a signal.
tapIds(k) — Let \(k \in \mathbb{R}^{K}\). Then \(\operatorname{tapIds}(k) = (0, 1, \dots, K-1) \in \mathbb{R}^{K}\). The positions along a kernel.
windowFrom(x, k) — Let \(X \in \mathbb{R}^{b \times L}\) and \(k \in \mathbb{R}^{K}\). Then \(\operatorname{windowFrom}(X, k) \in \mathbb{R}^{L \times K \times L}\) with \(\operatorname{windowFrom}(X, k)_{tjs} = t + j\). The tensor holding t + j at every entry.
windowTo(x, k) — Let \(X \in \mathbb{R}^{b \times L}\) and \(k \in \mathbb{R}^{K}\). Then \(\operatorname{windowTo}(X, k) \in \mathbb{R}^{L \times K \times L}\) with \(\operatorname{windowTo}(X, k)_{tjs} = s\). The tensor holding s at every entry.
windowOf(x, k) — Let \(X \in \mathbb{R}^{b \times L}\) and \(k \in \mathbb{R}^{K}\). Then \(\operatorname{windowOf}(X, k)_{tjs} = [\,s = t + j\,]\) for \(t, s = 0, \dots, L-1\) and \(j = 0, \dots, K-1\). The 0/1 window.
gathered(x, k) — Let \(X \in \mathbb{R}^{b \times L}\) and \(k \in \mathbb{R}^{K}\). Then \(\operatorname{gathered}(X, k)_{itj} = \sum_{s=0}^{L-1} X_{is} [\,s = t + j\,] = X_{i,\,t+j}\), which is \(0\) when \(t + j \ge L\). Every window of every signal, shape b×L×K.

Recurrence

A recurrent network carries a state through a sequence:

\[h_{t+1} = \tanh\!\left(h_t W_h + x_t\, w_x\right)\]

The definition names itself, so it compiles to a loop over time; there is no unrolling. Each sequence supplies one scalar reading per time step.

stateStep

stateStep(h, xt, wh, wx) · h state b×w · xt readings at this step vector b · wh state-to-state weights w×w · wx reading weights vector w → next state b×w

Let \(H \in \mathbb{R}^{b \times w}\) (the argument h) be the state, \(x \in \mathbb{R}^{b}\) (the argument xt) the readings at this step, \(W_h \in \mathbb{R}^{w \times w}\) (the argument wh) and \(w_x \in \mathbb{R}^{w}\) (the argument wx) the weights. Then

\[\operatorname{stateStep}(H, x, W_h, w_x) = \tanh\bigl(H W_h + x\, w_x^{\top}\bigr) .\]

One step of the recurrence above.

Intuition The state is a running summary of everything read so far. Each step mixes the old summary through \(W_h\), adds the new reading, and squashes the result with tanh so it stays bounded no matter how long the sequence is.

Use it when you write your own loop over time and need the state at every step, for per-step outputs, early exits or inspection. recurrent returns only the final state. Each sequence contributes one scalar per step, and there is no bias. For long-range memory use the gated cell (gateOf, cellMemory).

sequence-modellingtime-series

def stateStep(h, xt, wh, wx) ≔ tanh(einsum("ng,gh->nh", h, wh) + einsum("n,h->nh", xt, wx))

recurrent

recurrent(x, wh, wx, h0) · x sequences b×T (one sequence per row) · wh state weights w×w · wx reading weights vector w · h0 initial state b×w → final state b×w

Let \(X \in \mathbb{R}^{b \times T}\) (the argument x) hold one sequence per row, with time steps counted from zero, and write \(x^{(s)} \in \mathbb{R}^{b}\) for its column \(s = 0, 1, \dots, T-1\). Let \(W_h \in \mathbb{R}^{w \times w}\) and \(w_x \in \mathbb{R}^{w}\) (the arguments wh, wx) be the weights, and \(H_0 \in \mathbb{R}^{b \times w}\) (the argument h0) the initial state. Let

\[H_{s+1} = \operatorname{stateStep}\bigl(H_s, x^{(s)}, W_h, w_x\bigr) = \tanh\bigl(H_s W_h + x^{(s)} w_x^{\top}\bigr), \qquad s = 0, 1, \dots, T-1 .\]

Then

\[\operatorname{recurrent}(X, W_h, w_x, H_0) = H_T .\]

The state after the last time step, the sequence having been walked forwards.

Intuition A whole sequence is folded, step by step, into one fixed-size state vector, whatever its length. Recent readings dominate: each step multiplies older information by \(W_h\) and a tanh slope of at most 1, so its influence usually fades.

Use it when you want a fixed-length encoding of a sequence to feed a dense head for classification or regression, with weights you supply, since the module has no fit for recurrences. For dependencies over many steps, the plain recurrence forgets and its gradients vanish, so compose the gated cell instead.

sequence-modellingtime-seriesvanishing-gradients

def recurrent(x, wh, wx, h0) ≔ stateLoop(columns(x), x, siteIds(x), wh, wx, h0)

Building blocks

atTime(x, ids, t) — Let \(X \in \mathbb{R}^{b \times T}\), let \(\iota \in \mathbb{R}^{T}\) (the argument ids) be a vector of positions and \(t \in \{0, \dots, T-1\}\). Then \(\operatorname{atTime}(X, \iota, t)_i = \sum_{s} X_{is} [\,\iota_s = t\,]\), which is \(X_{i,t}\) when \(\iota = (0, 1, \dots, T-1)\). Column t of the sequences, read by contraction against a one-hot.
stateLoop(t, x, ids, wh, wx, h) — Let \(t \in \mathbb{N}\), \(X \in \mathbb{R}^{b \times T}\) with \(T = \operatorname{columns}(X)\), \(\iota \in \mathbb{R}^{T}\), \(W_h \in \mathbb{R}^{w \times w}\), \(w_x \in \mathbb{R}^{w}\) and \(H \in \mathbb{R}^{b \times w}\). Then \(\operatorname{stateLoop}(0, X, \iota, W_h, w_x, H) = H\) and \(\operatorname{stateLoop}(t, X, \iota, W_h, w_x, H) = \operatorname{stateLoop}\bigl(t - 1, X, \iota, W_h, w_x, \operatorname{stateStep}(H, \operatorname{atTime}(X, \iota, T - t), W_h, w_x)\bigr)\). Applies stateStep for the remaining t steps, walking the sequence forwards.

The gated cell

A long short-term memory (LSTM) cell is four copies of the same affine map followed by two elementwise rules. The map is stated once, with a sigmoid for gates and a tanh for the candidate, and a caller composes a step from them with separate weights for the forget gate f, input gate i, output gate o and candidate g:

\[c_{t+1} = f \odot c_t + i \odot g, \qquad h_{t+1} = o \odot \tanh(c_{t+1})\]

The memory update is additive, so a gradient flowing back through it is multiplied by f rather than by a saturating derivative — which is what lets the cell remember over long spans.

gateOf

gateOf(h, xt, wh, wx) · h state b×w · xt readings vector b · wh weights w×w · wx weights vector w → gate values in (0, 1) b×w

Let \(H \in \mathbb{R}^{b \times w}\) (the argument h) be the state, \(x \in \mathbb{R}^{b}\) (the argument xt) the readings, and \(W_h \in \mathbb{R}^{w \times w}\), \(w_x \in \mathbb{R}^{w}\) (the arguments wh, wx) the gate’s weights. Then

\[\operatorname{gateOf}(H, x, W_h, w_x) = \sigma\bigl(H W_h + x\, w_x^{\top}\bigr), \qquad \sigma(z) = \frac{1}{1 + e^{-z}} \text{ entrywise} .\]

Call it once each, with its own weights, for the forget, input and output gates.

Intuition A gate is a learned dimmer switch per memory unit. The sigmoid gives a value between 0 (closed) and 1 (open), computed from the current reading and the previous state, and it multiplies whatever flows through. The forget, input and output gates are the same map with different weights.

Use it when you assemble an LSTM step: call it three times with separate weights for \(f\), \(i\) and \(o\), then combine with candidateOf, cellMemory and cellHidden inside your own recursive definition, since no LSTM loop is provided. There is no bias, so with zero weights every gate starts half open.

sequence-modellinglong-range-dependenciesvanishing-gradients

def gateOf(h, xt, wh, wx)      ≔ sigmoid(einsum("ng,gh->nh", h, wh) + einsum("n,h->nh", xt, wx))

candidateOf

candidateOf(h, xt, wh, wx) · arguments as gateOf → candidate values in (−1, 1) b×w

Let \(H, x, W_h, w_x\) be as in gateOf. Then

\[\operatorname{candidateOf}(H, x, W_h, w_x) = \tanh\bigl(H W_h + x\, w_x^{\top}\bigr) .\]

The candidate memory.

Intuition The new content proposed for memory. It is the same map as stateStep, but in an LSTM it does not overwrite the state: the input gate decides how much of the proposal is written. Tanh’s range of \((-1, 1)\) lets it push each memory unit up or down.

Use it when you compute the \(g\) argument of cellMemory, with its own weights, separate from the three gates.

sequence-modellinglong-range-dependencies

def candidateOf(h, xt, wh, wx) ≔ tanh(einsum("ng,gh->nh", h, wh) + einsum("n,h->nh", xt, wx))

cellMemory

cellMemory(c, f, i, g) · c memory b×w · f forget gate b×w · i input gate b×w · g candidate b×w → next memory b×w

Let \(c, f, i, g \in \mathbb{R}^{b \times w}\) (the arguments of the same names) be the memory, the forget gate, the input gate and the candidate. Then

\[\operatorname{cellMemory}(c, f, i, g) = f \odot c + i \odot g .\]

The additive memory update.

Intuition Memory as a conveyor belt: keep a fraction \(f\) of what was there and add the gated proposal. Because the update adds rather than squashes, information and gradients pass back through a step multiplied only by \(f\). With \(f\) near 1 they survive for many steps, which is how an LSTM avoids vanishing gradients.

Use it when you need the core update of an LSTM step. With constant gates it is also a leaky integrator: \(i = 1 - f\) turns it into an exponential moving average of \(g\) with smoothing factor \(f\).

vanishing-gradientslong-range-dependenciestime-series

def cellMemory(c, f, i, g)     ≔ f · c + i · g

cellHidden

cellHidden(c, o) · c memory b×w · o output gate b×w → hidden state b×w

Let \(c, o \in \mathbb{R}^{b \times w}\) (the arguments c, o) be the memory and the output gate. Then

\[\operatorname{cellHidden}(c, o) = o \odot \tanh(c) .\]

The exposed state.

Intuition The memory can grow past \(\pm 1\) and holds things not needed right now. Tanh brings it back into range, and the output gate chooses which parts to reveal at this step, so the cell can remember something without acting on it yet.

Use it when you produce the hidden state of an LSTM step. That state feeds the next step’s gates and candidate, and a dense head reads it for predictions. Carry the memory \(c\) and this hidden state together through your loop.

sequence-modellinglong-range-dependencies

def cellHidden(c, o)           ≔ o · tanh(c)

The adversarial pair

A generative adversarial network pairs a critic D, which scores a sample as real, with a maker G, which produces samples from noise. Both are linear here, \(D(a) = \sigma(a\, w_d)\) and \(G(z) = z\, W_g\), which is enough to state the objectives exactly:

\[L_D = -\operatorname{mean}\ln D(x_{\text{real}}) - \operatorname{mean}\ln\bigl(1 - D(G(z))\bigr), \qquad L_G = -\operatorname{mean}\ln D(G(z))\]

The maker’s gradient passes through the critic: \(\partial L_G / \partial G(z) = -(1 - D)\, w_d\), so the critic’s weights appear in the maker’s update. Both are trained by simultaneous gradient descent. Nothing here guarantees that the pair converges; schedules, gradient penalties and similar fixes are training decisions outside this module.

Below, \(D(A) = \operatorname{critic}(A, w_d)\) for whichever critic weights \(w_d\) are in play.

critic

critic(a, wd) · a samples n×p · wd critic weights vector p → probability each sample is real vector n

Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) hold one sample per row and let \(w_d \in \mathbb{R}^{p}\) (the argument wd) be the critic weights. Then

\[\operatorname{critic}(A, w_d)_i = \sigma\Bigl(\sum_{j=1}^{p} A_{ij} (w_d)_j\Bigr) = \sigma(A w_d)_i, \qquad i = 1, \dots, n .\]

The logistic score for every row. There is no bias term.

Intuition Logistic regression that asks “is this sample real?”. It projects each sample onto one direction \(w_d\) and squashes the result to a probability. Being linear with no bias, it can tell real from fake only along a single direction through the origin.

Use it when you need the discriminator of the adversarial pair, or a quick two-sample check: if even a trained critic stays near 0.5, generated and real data are hard to separate along any direction. Append a column of ones to the samples if you need an offset.

generative-modellingclassificationdistribution-comparison

def critic(a, wd)      ≔ sigmoid(einsum("np,p->n", a, wd))

madeFrom

madeFrom(z, wg) · z noise m×q · wg maker weights q×p → generated samples m×p

Let \(Z \in \mathbb{R}^{m \times q}\) (the argument z) be the noise and \(W_g \in \mathbb{R}^{q \times p}\) (the argument wg) the maker weights. Then

\[\operatorname{madeFrom}(Z, W_g) = Z W_g .\]

The maker’s samples, one per row of noise.

Intuition The maker reshapes random noise into samples. Being linear, it can only stretch, rotate and mix the noise: Gaussian noise becomes Gaussian samples whose covariance is shaped by \(W_g\), centred where the noise is centred.

Use it when you generate samples, either during training or afterwards from new noise with weights from ganMaker. To learn a non-zero mean, append a column of ones to z. Distributions with several clusters are beyond this linear maker.

generative-modellingsynthetic-data

def madeFrom(z, wg)    ≔ einsum("mq,qp->mp", z, wg)

criticLoss

criticLoss(real, fake, wd) · real real samples n×p · fake generated samples m×p · wd critic weights vector p → loss scalar

Let \(X \in \mathbb{R}^{n \times p}\) (the argument real) be the real samples, \(F \in \mathbb{R}^{m \times p}\) (the argument fake) the generated samples, and \(w_d \in \mathbb{R}^{p}\) (the argument wd) the critic weights. Then

\[\operatorname{criticLoss}(X, F, w_d) = -\frac{1}{n}\sum_{i=1}^{n} \ln \max\bigl(D(X)_i,\ 10^{-12}\bigr) - \frac{1}{m}\sum_{i=1}^{m} \ln \max\bigl(1 - D(F)_i,\ 10^{-12}\bigr) .\]

The critic’s binary cross-entropy \(L_D\), with probabilities floored at \(10^{-12}\) before the logarithm.

Intuition The cross-entropy of a classifier that labels real samples 1 and fakes 0. It is low when the critic separates them confidently, and \(2 \ln 2 \approx 1.386\) when it says 0.5 to everything, which is where a perfect maker would leave it.

Use it when you monitor adversarial training alongside makerLoss. A critic loss that drops toward zero means the critic is winning and the maker’s gradient is weak. A loss stuck near 1.386 from the start can mean the critic is too weak to teach anything, not that the fakes are good.

generative-modellingtraining-diagnostics

def criticLoss(real, fake, wd) ≔ 0 − mean(ln(floorOf(critic(real, wd)))) − mean(ln(floorOf(1 − critic(fake, wd))))

makerLoss

makerLoss(fake, wd) · fake generated samples m×p · wd critic weights vector p → loss scalar

Let \(F \in \mathbb{R}^{m \times p}\) (the argument fake) be the generated samples and \(w_d \in \mathbb{R}^{p}\) (the argument wd) the critic weights. Then

\[\operatorname{makerLoss}(F, w_d) = -\frac{1}{m}\sum_{i=1}^{m} \ln \max\bigl(D(F)_i,\ 10^{-12}\bigr) .\]

The maker’s loss \(L_G\) (the non-saturating form), floored the same way.

Intuition The maker wants the critic to call its samples real. Minimising \(\ln(1 - D)\) directly gives almost no gradient early on, when the critic confidently rejects fakes and \(D\) is near 0. The non-saturating form \(-\ln D\) is steepest exactly there, so the maker learns fastest when it is worst.

Use it when you track the maker’s progress next to criticLoss. If one loss collapses while the other climbs steadily, the pair is out of balance: adjust the learning rate or the number of critic steps using criticStep and makerStep.

generative-modellingtraining-diagnosticsvanishing-gradients

def makerLoss(fake, wd)        ≔ 0 − mean(ln(floorOf(critic(fake, wd))))

criticStep

criticStep(r, real, z, wg, wd) · r learning rate scalar · real real samples n×p · z noise m×q · wg maker weights q×p · wd critic weights vector p → updated critic weights vector p

Let \(r > 0\) (the argument r), \(X \in \mathbb{R}^{n \times p}\) (the argument real), \(Z \in \mathbb{R}^{m \times q}\) (the argument z), \(W_g \in \mathbb{R}^{q \times p}\) (the argument wg) and \(w_d \in \mathbb{R}^{p}\) (the argument wd). Then

\[\operatorname{criticStep}(r, X, Z, W_g, w_d) = w_d - r\, \operatorname{criticGrad}(X, Z W_g, w_d) .\]

One gradient-descent step on \(L_D\), with the fake samples made from z by the current maker.

Intuition One nudge of the critic toward scoring real samples higher and the current maker’s fakes lower. The fakes are made from z by the maker as it is now, so the critic always trains against its present opponent.

Use it when you want a training schedule that ganCritic does not provide: several critic steps per maker step, separate learning rates for the two players, or fresh noise on every step. Loop over it and makerStep in your own recursive definition.

generative-modellinggradient-based-trainingtraining-control

def criticStep(r, real, z, wg, wd) ≔ descend(wd, criticGrad(real, madeFrom(z, wg), wd), r)

makerStep

makerStep(r, z, wg, wd) · r learning rate scalar · z noise m×q · wg maker weights q×p · wd critic weights vector p → updated maker weights q×p

Let \(r > 0\) (the argument r), \(Z \in \mathbb{R}^{m \times q}\) (the argument z), \(W_g \in \mathbb{R}^{q \times p}\) (the argument wg) and \(w_d \in \mathbb{R}^{p}\) (the argument wd). Then

\[\operatorname{makerStep}(r, Z, W_g, w_d) = W_g - r\, \operatorname{makerGrad}(Z, W_g, w_d) .\]

One gradient-descent step on \(L_G\) against the current critic.

Intuition One nudge of the maker toward samples the critic scores as more real. The gradient contains the critic’s weights \(w_d\), so the maker moves its samples in the very direction the critic uses to tell them apart.

Use it when you drive adversarial training yourself with criticStep. Passing new noise on each call keeps the maker from fitting one fixed batch. Compute both steps from the same current weights for a simultaneous update, which is what ganCritic and ganMaker do.

generative-modellinggradient-based-trainingtraining-control

def makerStep(r, z, wg, wd)        ≔ descend(wg, makerGrad(z, wg, wd), r)

ganCritic

ganCritic(t, r, real, z, wg, wd) · t number of steps scalar · r learning rate scalar · real real samples n×p · z noise m×q · wg starting maker weights q×p · wd starting critic weights vector p → critic weights vector p

Let \(t \in \mathbb{N}\) (the argument t), \(r > 0\) (the argument r), \(X \in \mathbb{R}^{n \times p}\) (the argument real), \(Z \in \mathbb{R}^{m \times q}\) (the argument z), and let \(W_g \in \mathbb{R}^{q \times p}\), \(w_d \in \mathbb{R}^{p}\) (the arguments wg, wd) be the starting weights. Let \(W_g^{(0)} = W_g\), \(w_d^{(0)} = w_d\) and, for \(s = 0, 1, \dots, t-1\),

\[W_g^{(s+1)} = \operatorname{makerStep}\bigl(r, Z, W_g^{(s)}, w_d^{(s)}\bigr), \qquad w_d^{(s+1)} = \operatorname{criticStep}\bigl(r, X, Z, W_g^{(s)}, w_d^{(s)}\bigr).\]

Then

\[\operatorname{ganCritic}(t, r, X, Z, W_g, w_d) = w_d^{(t)} .\]

The critic’s weights after t simultaneous steps of both players.

Intuition The critic at the end of a two-player game run for t simultaneous steps. Its weight vector is the direction that, at the end, still best separates real from generated samples.

Use it when you want to inspect what distinguishes the fakes, or to score samples with critic. The run reuses the same noise z and one learning rate throughout, and nothing guarantees it settles: losses can oscillate, so check criticLoss at the result. Calling it and ganMaker runs the game twice with identical results.

generative-modellingtraining-diagnostics

def ganCritic(t, r, real, z, wg, wd) ≔ criticLoop(t, r, real, z, wg, wd)

ganMaker

ganMaker(t, r, real, z, wg, wd) · arguments as ganCritic → maker weights q×p

Let \(t, r, X, Z, W_g, w_d\) and the iterate \(W_g^{(t)}\) be as in ganCritic. Then

\[\operatorname{ganMaker}(t, r, X, Z, W_g, w_d) = W_g^{(t)} .\]

The maker’s weights after the same t steps.

Intuition The trained maker is the game’s useful product: a map from noise to samples meant to resemble the real data.

Use it when you want to generate new samples later with madeFrom(z2, W) from fresh noise z2, without retraining. Pass the same arguments as ganCritic if you need the matching critic. Training used a single fixed noise batch, so check new samples against real data statistics before relying on them.

generative-modellingsynthetic-data

def ganMaker(t, r, real, z, wg, wd)  ≔ makerLoop(t, r, real, z, wg, wd)

ganSamples

ganSamples(t, r, real, z, wg, wd) · arguments as ganCritic → generated samples m×p

Let \(t, r, X, Z, W_g, w_d\) and the iterate \(W_g^{(t)}\) be as in ganCritic. Then

\[\operatorname{ganSamples}(t, r, X, Z, W_g, w_d) = Z\, W_g^{(t)} .\]

Trains the pair for t steps and returns the trained maker’s samples from z.

Intuition Train the pair, then show what the maker makes from the noise it trained on. It is the quickest way to see whether the game went anywhere.

Use it when you want a one-call look at generated data: compare its means and covariances with real. These samples come from the training noise z itself, so for fresh samples use madeFrom with new noise and the weights from ganMaker. Every call retrains from the starting weights.

synthetic-datagenerative-modelling

def ganSamples(t, r, real, z, wg, wd) ≔ madeFrom(z, ganMaker(t, r, real, z, wg, wd))

Building blocks

criticGrad(real, fake, wd) — Let \(X \in \mathbb{R}^{n \times p}\), \(F \in \mathbb{R}^{m \times p}\) and \(w_d \in \mathbb{R}^{p}\). Then \(\operatorname{criticGrad}(X, F, w_d) = \frac{1}{m} F^{\top} D(F) - \frac{1}{n} X^{\top} \bigl(\mathbf{1}_n - D(X)\bigr)\). The gradient of the critic’s loss with respect to its weights; it is the gradient of the unfloored loss.
makerGrad(z, wg, wd) — Let \(Z \in \mathbb{R}^{m \times q}\), \(W_g \in \mathbb{R}^{q \times p}\) and \(w_d \in \mathbb{R}^{p}\). Then \(\operatorname{makerGrad}(Z, W_g, w_d) = -\frac{1}{m} Z^{\top} \bigl(\mathbf{1}_m - D(Z W_g)\bigr)\, w_d^{\top}\). The gradient of the maker’s loss with respect to its weights, passing through the critic; like criticGrad it ignores the floor.
criticLoop(t, r, real, z, wg, wd) — Let \(t \in \mathbb{N}\), \(r > 0\), \(X \in \mathbb{R}^{n \times p}\), \(Z \in \mathbb{R}^{m \times q}\), \(W_g \in \mathbb{R}^{q \times p}\) and \(w_d \in \mathbb{R}^{p}\). Then \(\operatorname{criticLoop}(0, r, X, Z, W_g, w_d) = w_d\) and \(\operatorname{criticLoop}(t, r, X, Z, W_g, w_d) = \operatorname{criticLoop}\bigl(t - 1, r, X, Z, \operatorname{makerStep}(r, Z, W_g, w_d), \operatorname{criticStep}(r, X, Z, W_g, w_d)\bigr)\). Steps both players t times and answers the critic’s weights.
makerLoop(t, r, real, z, wg, wd) — Let \(t, r, X, Z, W_g, w_d\) be as in criticLoop. Then \(\operatorname{makerLoop}(0, r, X, Z, W_g, w_d) = W_g\) and \(\operatorname{makerLoop}(t, r, X, Z, W_g, w_d) = \operatorname{makerLoop}\bigl(t - 1, r, X, Z, \operatorname{makerStep}(r, Z, W_g, w_d), \operatorname{criticStep}(r, X, Z, W_g, w_d)\bigr)\). Steps both players t times and answers the maker’s weights.

Examples

The examples below use two files, Layers.hb and Fits.hb:

open Prelude.Networks

fn layerThrough(a, w) ↦ dense(a, w)
fn normalisedRows(a) ↦ layerNorm(a)
fn attend(q, k, v, d) ↦ attentionOf(q, k, v, d)
fn attendCausally(q, k, v, d) ↦ causalAttention(q, k, v, d)
fn kept(a, u, p) ↦ dropped(a, u, p)
fn filtered(x, k) ↦ convolve(x, k)
fn pooledMax(x, k) ↦ maxPool(x, k)

open Prelude.Networks

fn trainedOutputs(t, r, x, y, w1, w2) ↦ mlpFit(t, r, x, y, w1, w2)
fn trainedQuery(t, r, x, y, w1, w2, z) ↦ mlpPredict(t, r, x, y, w1, w2, z)
fn unrolled(x, wh, wx, h0) ↦ recurrent(x, wh, wx, h0)
fn squeezed(t, r, x, w1, w2) ↦ autoFit(t, r, x, w1, w2)
fn criticLossOf(real, fake, wd) ↦ criticLoss(real, fake, wd)

Matrix results are shown row by row.

Push a batch through a layer and normalise it

Two input rows pass through a 2×3 weight matrix. Layer normalisation then standardises each row on its own: [1, 2, 3] has mean 2 and variance 2/3, so it becomes about \([-1.2247, 0, 1.2247]\).

using MLambda.Hilbert.Runtime;
using Examples;

var a = (Tensor)new double[,] { { 1, 2 }, { 3, 4 } };
var w = (Tensor)new double[,] { { 1, 0, 1 }, { 0, 1, 1 } };
var layer = Layers.LayerThrough(a, w);

var rows = (Tensor)new double[,] { { 1, 2, 3 }, { 10, 10, 40 } };
var normal = Layers.NormalisedRows(rows);

Output (verified, six decimals):

dense     = [[1, 2, 3], [3, 4, 7]]
layerNorm = [[-1.224736, 0, 1.224736], [-0.707107, -0.707107, 1.414214]]

Attend over three tokens, with and without a causal mask

Queries and keys are the 3×3 identity and the values are 1, 2 and 3, with d = 1. Without a mask, token 0 weighs itself by e and the others by 1: \((e + 2 + 3)/(e + 2) = 1.6358\). With the causal mask, token 0 can only read itself (1), and token 1 reads tokens 0 and 1 with weights 1 and e: \((1 + 2e)/(1 + e) = 1.7311\). The last token sees everything either way.

var q = (Tensor)new double[,] { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
var v = (Tensor)new double[,] { { 1 }, { 2 }, { 3 } };

var open   = Layers.Attend(q, q, v, 1.0);
var causal = Layers.AttendCausally(q, q, v, 1.0);

Output (verified, six decimals):

attention = [[1.635825], [2], [2.364175]]
causal    = [[1], [1.731059], [2.364175]]

Apply dropout with explicit draws

With drop probability 0.5, the entries whose draw is at least 0.5 survive and are doubled.

var acts  = (Tensor)new double[,] { { 1, 2 }, { 3, 4 } };
var draws = (Tensor)new double[,] { { 0.1, 0.7 }, { 0.9, 0.3 } };

var kept = Layers.Kept(acts, draws, 0.5);

Output (verified):

dropped = [[0, 4], [6, 0]]

Filter and pool a signal

The kernel [1, -1] takes the difference \(x_t - x_{t+1}\); at the last position the padding supplies zero, so it reads 4. Max pooling with a window of two reads the larger of each neighbouring pair, and the last window is [4, 0].

var signal = (Tensor)new double[,] { { 1, 2, 3, 4 } };
var diff = Layers.Filtered(signal, (Tensor)new[] { 1d, -1d });

var bumpy = (Tensor)new double[,] { { 1, 5, 2, 4 } };
var pooled = Layers.PooledMax(bumpy, (Tensor)new[] { 1d, 1d });

Output (verified):

convolve = [[-1, -1, -1, 4]]
maxPool  = [[5, 5, 4, 4]]

Fit a small network to a line

Four points on \(y = 2x + 1\), with a column of ones as the bias. Before training (t = 0) the network outputs its starting guess; after 500 steps at rate 0.05 it reproduces the targets, and it extrapolates to \(x = 4\) as 9. Each call retrains from the same starting weights, so the result is deterministic.

var x  = (Tensor)new double[,] { { 0, 1 }, { 1, 1 }, { 2, 1 }, { 3, 1 } };
var y  = (Tensor)new double[,] { { 1 }, { 3 }, { 5 }, { 7 } };
var w1 = (Tensor)new double[,] { { 0.5, 0.3 }, { 0.2, 0.4 } };
var w2 = (Tensor)new double[,] { { 0.5 }, { 0.5 } };

var before = Fits.TrainedOutputs(0.0, 0.05, x, y, w1, w2);
var after  = Fits.TrainedOutputs(500.0, 0.05, x, y, w1, w2);
var ahead  = Fits.TrainedQuery(500.0, 0.05, x, y, w1, w2, (Tensor)new double[,] { { 4, 1 } });

Output (verified, six decimals):

t = 0     = [[0.3], [0.7], [1.1], [1.5]]
t = 500   = [[1], [3], [5], [7]]
query x=4 = [[9]]

Carry a state through a sequence

One sequence [1, 0, 0], a single state unit with \(W_h = 0.5\), \(w_x = 1\) and \(h_0 = 0\). The states are \(\tanh(1) = 0.7616\), then \(\tanh(0.3808) = 0.3634\), then \(\tanh(0.1817) = 0.1797\).

var seq = (Tensor)new double[,] { { 1, 0, 0 } };
var h = Fits.Unrolled(seq, (Tensor)new double[,] { { 0.5 } }, (Tensor)new[] { 1d }, (Tensor)new double[,] { { 0 } });

Output (verified, six decimals):

recurrent = [[0.179726]]

Compress points on a line through one unit

Three points lie on the line \(x_1 = x_2\), so a one-unit linear autoencoder can reconstruct them exactly. The starting weights reconstruct them poorly; 300 steps at rate 0.1 recover them.

var pts = (Tensor)new double[,] { { 1, 1 }, { 2, 2 }, { -1, -1 } };
var enc = (Tensor)new double[,] { { 0.5 }, { 0.1 } };
var dec = (Tensor)new double[,] { { 0.3, 0.4 } };

var start  = Fits.Squeezed(0.0, 0.1, pts, enc, dec);
var fitted = Fits.Squeezed(300.0, 0.1, pts, enc, dec);

Output (verified, six decimals):

t = 0   = [[0.18, 0.24], [0.36, 0.48], [-0.18, -0.24]]
t = 300 = [[1, 1], [2, 2], [-1, -1]]

Score a critic

With zero weights the critic says 0.5 for every sample, so its loss is \(2 \ln 2 = 1.3863\). With weights [2, -2], the real sample [1, 0] scores \(\sigma(2)\) and the fake [0, 1] scores \(\sigma(-2)\), and the loss falls to \(-2 \ln \sigma(2) = 0.2539\).

var real = (Tensor)new double[,] { { 1, 2 }, { 3, 4 } };
var fake = (Tensor)new double[,] { { 0, 1 } };
var undecided = Fits.CriticLossOf(real, fake, (Tensor)new[] { 0d, 0d });

var sharp = Fits.CriticLossOf((Tensor)new double[,] { { 1, 0 } }, fake, (Tensor)new[] { 2d, -2d });

Output (verified, six decimals):

zero weights = 1.386294
[2, -2]      = 0.253856