Learning
Activations, losses with their gradients, and the optimiser steps (momentum, RMSProp, Adam, soft thresholding) that training loops are built from.
This module holds the pieces a trained model is assembled from: activation functions and their derivatives, the common losses together with their gradients, and the update rules of the usual optimisers — plain gradient descent, momentum, RMSProp, Adam and the soft-thresholding step of ISTA.
Each optimiser is written as a step, a function of the current parameters, the gradient and whatever state the optimiser carries. The loop that repeats the step is written once per model, as a definition that names itself (see the last example below, and the fits in Prelude.Regression). There is no shared “train” loop, because a loop in Hilbert cannot take a function as a parameter; the cost is a three-line loop per model.
Conventions. Losses take a vector of predictions or margins and a vector of targets of the same length, and are averaged over the rows rather than summed, so the same learning rate works whatever the size of the data. Gradients returned by the d… functions include that same division by the row count. Every logarithm is taken of a probability floored at \(10^{-12}\), so a confident wrong prediction yields a large finite loss rather than NaN. The activation functions are written so that neither tail overflows: sigmoid(800) is 1 and softplus(800) is 800.
open Prelude.Learning
Safe logarithms
A probability of exactly zero is the one value whose logarithm must never be taken; a probability near zero is meaningful and is kept.
floorOf
floorOf(p) · p probabilities any shape → floored probabilities same shape
Let \(p\) be a tensor of any shape with entries \(p_\iota\). Then
\[\operatorname{floorOf}(p)_\iota = \max(p_\iota, 10^{-12}) .\]Raises every value below \(10^{-12}\) to \(10^{-12}\) and leaves the rest unchanged. Use it before ln whenever a probability may be zero.
Intuition A probability of exactly 0 has logarithm \(-\infty\), and \(-\infty\) times a zero label is NaN, which poisons a whole sum. Flooring at \(10^{-12}\) caps the cost of a certain-but-wrong prediction at about 27.6 (\(\ln 10^{12}\)) and leaves every realistic probability untouched.
Use it when you take ln of a probability you computed yourself: a custom likelihood, a log-score, the entropy of a histogram. The losses on this page already apply it. It is a floor, not a normaliser: it does not rescale values to sum to one, and it also lifts negative inputs to \(10^{-12}\), which can hide a bug upstream.
numerical-stabilitymodel-evaluation
def floorOf(p) ≔ max(p, 1e-12)
Activations
Element-wise non-linearities and their derivatives. Every function here applies to a real or to a tensor of any shape.
sigmoid
sigmoid(z) · z logits any shape → values in (0, 1) same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\). Then
\[\operatorname{sigmoid}(z)_\iota = \sigma(z_\iota) = \frac{1}{1 + e^{-z_\iota}} .\]The logistic function, written as \(1/(1 + e^{-z})\) so that large positive inputs give exactly 1 instead of \(\infty/\infty\).
\[\sigma(z) = \frac{1}{1 + e^{-z}}\]Intuition Maps any real score to a value in \((0, 1)\): 0 goes to 0.5, the curve is symmetric (\(\sigma(-z) = 1 - \sigma(z)\)), and it flattens at both ends. It is the inverse of the log-odds, which is why it turns a linear score into a probability.
Use it when you turn a logit into a probability (logistic regression, a binary output unit) or need a gate between 0 and 1. For several mutually exclusive classes use a softmax (softmaxRows in Prelude.Kernels), or give the logits straight to crossEntropy. As a hidden-layer activation it saturates and its slope is at most 0.25, which slows deep networks; prefer relu.
probability-calibrationclassificationneural-network-training
def sigmoid(z) ≔ 1 / (1 + exp(0 − z))
relu
relu(z) · z inputs any shape → outputs same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\). Then
\[\operatorname{relu}(z)_\iota = \max(z_\iota, 0) .\]The rectified linear unit, \(\max(z, 0)\).
Intuition Passes positive signals unchanged and blocks negative ones. A network of such units is piecewise linear, and where a unit is active its slope is exactly 1, so gradients do not shrink as they flow back through many layers.
Use it when you build hidden layers of a neural network, where it is the default. It is also a hinge-shaped feature, “the amount above a threshold”: \(\max(x - c, 0)\) terms give piecewise-linear regression. A unit whose input stays negative gets no gradient and stops learning; leaky or softplus avoid that.
neural-network-trainingnonlinear-regressionfeature-engineering
def relu(z) ≔ max(z, 0)
leaky
leaky(z, a) · z inputs any shape · a slope for negative inputs scalar → outputs same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\), and let \(a \in \mathbb{R}\). Then
\[\operatorname{leaky}(z, a)_\iota = \max(z_\iota, 0) + a \min(z_\iota, 0) .\]The leaky rectifier: \(z\) for positive inputs and \(a z\) for negative ones.
Intuition A relu with a small slope a on the negative side. An inactive unit still passes a fraction of the gradient, so it can recover instead of staying switched off for good.
Use it when many relu units output 0 for every input (“dead” units), which often follows a too-large learning rate. a = 0.01 is the common choice; a = 0 gives back relu and a = 1 is the identity. Backpropagate through it with dLeaky and the same a.
neural-network-trainingvanishing-gradients
def leaky(z, a) ≔ max(z, 0) + a · min(z, 0)
softplus
softplus(z) · z inputs any shape → outputs same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\). Then
\[\operatorname{softplus}(z)_\iota = \max(z_\iota, 0) + \ln\!\left(1 + e^{-\lvert z_\iota \rvert}\right) = \ln\!\left(1 + e^{z_\iota}\right) .\]The smooth rectifier \(\ln(1 + e^z)\), computed through an identity with a bounded exponent so it never overflows:
\[\operatorname{softplus}(z) = \max(z, 0) + \ln\!\left(1 + e^{-|z|}\right)\]Intuition A smoothed relu: close to 0 for very negative inputs, close to \(z\) for large ones, rounded near the origin (\(\operatorname{softplus}(0) = \ln 2\)). Its output is always positive and its derivative is exactly the sigmoid.
Use it when you need a strictly positive quantity from an unconstrained parameter, such as a variance, rate or scale (\(\sigma = \operatorname{softplus}(\theta)\)), or a smooth activation where the kink of relu gets in the way. It also writes log-likelihoods stably: \(-\ln \sigma(z) = \operatorname{softplus}(-z)\). It costs a little more than relu.
positivity-constraintsnumerical-stabilityneural-network-training
def softplus(z) ≔ max(z, 0) + ln(1 + exp(0 − abs(z)))
dSigmoid
dSigmoid(z) · z logits any shape → derivative same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\), and \(\sigma(s) = 1/(1 + e^{-s})\). Then
\[\operatorname{dSigmoid}(z)_\iota = \sigma(z_\iota)\,\big(1 - \sigma(z_\iota)\big) = \sigma'(z_\iota) .\]The derivative of the sigmoid, \(\sigma(z)(1 - \sigma(z))\).
Intuition How far the sigmoid’s output moves per unit of input: largest (0.25) at 0 and vanishing in both tails. A saturated unit barely learns, and stacking several of these factors shrinks gradients quickly.
Use it when you backpropagate through a sigmoid you applied yourself; pass the input z, not the output. When the sigmoid feeds a log-loss, skip it: dLogistic already gives the combined gradient \((p - y)/n\), which stays large when a prediction is confidently wrong, whereas multiplying by dSigmoid would nearly erase it.
gradient-based-trainingvanishing-gradientsneural-network-training
def dSigmoid(z) ≔ sigmoid(z) · (1 − sigmoid(z))
dRelu
dRelu(z) · z inputs any shape → 0/1 mask same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\). Then
\[\operatorname{dRelu}(z)_\iota = [\,z_\iota > 0\,] .\]The derivative of relu: 1 where \(z > 0\) and 0 elsewhere (including at zero).
Intuition A gate: the gradient flows back unchanged through units that were active and is blocked at the rest. relu has no derivative at exactly 0, and this picks 0 there.
Use it when you backpropagate through relu: multiply the incoming gradient by this mask entrywise. If a unit’s mask is 0 for every row of the data, the unit is dead and will never change; switch to leaky (with dLeaky) or lower the learning rate.
gradient-based-trainingneural-network-training
def dRelu(z) ≔ z > 0
dLeaky
dLeaky(z, a) · z inputs any shape · a slope for negative inputs scalar → derivative same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\), and let \(a \in \mathbb{R}\). Then
\[\operatorname{dLeaky}(z, a)_\iota = [\,z_\iota > 0\,] + a\,[\,z_\iota \le 0\,] .\]The derivative of leaky: 1 where \(z > 0\) and \(a\) where \(z \le 0\).
Intuition The gate of dRelu with a leak: active units pass the gradient in full and inactive ones pass the fraction a, so no unit is ever completely cut off from learning.
Use it when you backpropagate through leaky, with the same a as the forward pass. An input of exactly 0 counts as inactive and gets slope a, consistent with dRelu giving 0 there.
gradient-based-trainingvanishing-gradientsneural-network-training
def dLeaky(z, a) ≔ (z > 0) + a · (z ≤ 0)
dTanh
dTanh(z) · z inputs any shape → derivative same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\). Then
\[\operatorname{dTanh}(z)_\iota = 1 - \tanh^2 z_\iota .\]The derivative of the hyperbolic tangent, \(1 - \tanh^2 z\).
Intuition tanh is a zero-centred sigmoid with range \((-1, 1)\). Its slope is 1 at the origin and falls to 0 as the unit saturates, and the derivative can be read off the output \(h = \tanh z\) as \(1 - h^2\).
Use it when you backpropagate through tanh: recurrent cells, small networks, or outputs that must stay in \((-1, 1)\). Pass the pre-activation z; if you kept the output instead, \(1 - h^2\) is the same number. Its peak slope of 1 (against 0.25 for the sigmoid) and centred outputs train better than the sigmoid in hidden layers, but deep stacks still favour relu.
gradient-based-trainingvanishing-gradientsneural-network-training
def dTanh(z) ≔ 1 − tanh(z) · tanh(z)
Losses, each a number
Each loss reduces a batch to one number, averaged over the rows. The matching d… function is the gradient of the loss with respect to the predictions (or margins), already divided by the row count, ready to be carried back through a model.
squaredLoss
squaredLoss(p, y) · p predictions n · y targets n → loss scalar
Let \(p \in \mathbb{R}^n\) be the predictions and \(y \in \mathbb{R}^n\) the targets. Then
\[\operatorname{squaredLoss}(p, y) = \frac{1}{2n} \sum_{i=1}^{n} (p_i - y_i)^2 .\]Half the mean squared error:
\[L = \frac{1}{2n} \sum_i (p_i - y_i)^2\]Intuition The average squared miss, halved so that its gradient is simply the residual. Squaring makes a miss of 2 cost four times a miss of 1, so the best constant prediction is the mean and a few large errors dominate the total.
Use it when you fit or evaluate a regression with roughly Gaussian noise. As a reported metric, twice this value is the mean squared error and \(\sqrt{2L}\) the RMSE. With outliers or heavy-tailed noise use huberLoss; to predict a median or another quantile use pinballLoss.
regressionmodel-evaluationparameter-estimation
def squaredLoss(p, y) ≔ einsum("i,i->", p − y, p − y) / (2 · count(p))
dSquared
dSquared(p, y) · p predictions n · y targets n → gradient n
Let \(p, y \in \mathbb{R}^n\). Then
\[\operatorname{dSquared}(p, y) = \frac{1}{n}(p - y) = \nabla_p \operatorname{squaredLoss}(p, y) .\]The gradient of squaredLoss with respect to the predictions, \((p - y)/n\).
Intuition The residual divided by \(n\): each prediction is pushed toward its target in proportion to how far off it is, and a perfect prediction is not pushed at all.
Use it when you train by gradient descent under squaredLoss. Carry it back to the parameters: for a linear model \(p = Xw\) the parameter gradient is \(X^\top\) times this vector, as in the slope example below. The division by \(n\) means the learning rate need not change when the batch size does.
gradient-based-trainingregression
def dSquared(p, y) ≔ (p − y) / count(p)
logisticLoss
logisticLoss(p, y) · p predicted probabilities n · y 0/1 labels n → loss scalar
Let \(p \in \mathbb{R}^n\) be predicted probabilities, \(y \in \{0, 1\}^n\) the labels. Then
\[\operatorname{logisticLoss}(p, y) = -\frac{1}{n} \sum_{i=1}^{n} \Big[ y_i \ln \max(p_i, 10^{-12}) + (1 - y_i) \ln \max(1 - p_i, 10^{-12}) \Big] .\]The mean binary cross-entropy, with both logarithms floored (shown here without the floor):
\[L = -\frac{1}{n} \sum_i \big[ y_i \ln p_i + (1 - y_i) \ln(1 - p_i) \big]\]Intuition The average negative log-probability given to the true label. A coin-flip guess of 0.5 costs \(\ln 2 \approx 0.69\), a confident correct prediction costs almost nothing, and a confident wrong one is expensive: predicting 0.01 for a true 1 costs 4.6. It is the Bernoulli log-likelihood, so minimising it rewards honest probabilities.
Use it when you train or evaluate a binary probabilistic classifier (the “log-loss” metric), or compare how well models are calibrated. It takes probabilities; if you hold logits, apply sigmoid first, or use crossEntropy with two columns, which stays accurate for extreme logits. A certain-but-wrong row is capped by the floor at about 27.6.
classificationprobability-calibrationmodel-evaluation
def logisticLoss(p, y) ≔ 0 − einsum("i->", y · ln(floorOf(p)) + (1 − y) · ln(floorOf(1 − p))) / count(p)
dLogistic
dLogistic(p, y) · p predicted probabilities n · y 0/1 labels n → gradient n
Let \(p, y \in \mathbb{R}^n\). Then
\[\operatorname{dLogistic}(p, y) = \frac{1}{n}(p - y) .\]This is the gradient of logisticLoss with respect to the logits when \(p = \sigma(z)\) (away from the floor).
Intuition The sigmoid’s slope and the log-loss’s slope cancel, leaving just prediction minus label: the push is proportional to how wrong the probability is, and it does not vanish for a confidently wrong one.
Use it when you train logistic regression or a sigmoid output unit by gradient descent. It is the gradient with respect to the logit \(z\), so carry it straight into the weights (\(X^\top\) times it) without multiplying by dSigmoid. It is not the gradient with respect to p. The formula matches dSquared; what differs is that p here must be \(\sigma(z)\).
gradient-based-trainingclassification
def dLogistic(p, y) ≔ (p − y) / count(p)
hingeLoss
hingeLoss(m, y) · m margins (raw scores) n · y labels −1 or +1 n → loss scalar
Let \(m \in \mathbb{R}^n\) be the margins and \(y \in \{-1, +1\}^n\) the labels. Then
\[\operatorname{hingeLoss}(m, y) = \frac{1}{n} \sum_{i=1}^{n} \max(1 - y_i m_i, 0) .\]The mean hinge loss of a support-vector machine, \(\frac{1}{n}\sum_i \max(1 - y_i m_i, 0)\). Labels are −1 and +1, not 0 and 1.
Intuition Zero for rows on the correct side with a margin of at least 1, rising linearly as a row approaches or crosses the boundary. Unlike log-loss it stops caring once a row is safely classified, so only the hard cases near the boundary count.
Use it when you train or evaluate a max-margin linear classifier (svmFit in Prelude.Regression minimises this plus an L2 penalty) and need decisions rather than probabilities. Labels must be −1/+1 and m raw scores, not probabilities. For probabilistic outputs use logisticLoss.
classificationmodel-evaluation
def hingeLoss(m, y) ≔ einsum("i->", max(1 − y · m, 0)) / count(m)
dHinge
dHinge(m, y) · m margins n · y labels −1 or +1 n → subgradient n
Let \(m \in \mathbb{R}^n\) and \(y \in \{-1, +1\}^n\). Then, for \(i = 1, \dots, n\),
\[\operatorname{dHinge}(m, y)_i = -\frac{1}{n}\, y_i\, [\,y_i m_i < 1\,] .\]A subgradient of hingeLoss with respect to the margins: \(-y_i/n\) where \(y_i m_i < 1\) and 0 elsewhere.
Intuition Only rows inside the margin push: each moves its own score toward its label by \(1/n\), and every other row contributes 0. The hinge has a corner at \(y_i m_i = 1\) with no derivative, and this picks 0 there.
Use it when you train a hinge-loss model by subgradient descent: carry it back through the design (\(X^\top\) times it), as svmFit in Prelude.Regression does. Its size does not shrink as the fit improves, so the iterates keep jittering near the optimum; use a small or decaying learning rate.
gradient-based-trainingclassification
def dHinge(m, y) ≔ 0 − y · (y · m < 1) / count(m)
huberLoss
huberLoss(r, d) · r residuals n · d threshold scalar → loss scalar
Let \(r \in \mathbb{R}^n\) be the residuals and \(d > 0\) the threshold. Then
\[\operatorname{huberLoss}(r, d) = \frac{1}{n} \sum_{i=1}^{n} \Big( [\,\lvert r_i \rvert \le d\,]\, \frac{r_i^2}{2} + [\,\lvert r_i \rvert > d\,]\, d \big(\lvert r_i \rvert - \tfrac{d}{2}\big) \Big) .\]The mean Huber loss: quadratic for residuals up to d in size and linear beyond, which makes it less sensitive to outliers than the squared loss.
Intuition Squared error for small residuals and absolute error beyond d, joined without a kink at \(\lvert r \rvert = d\). An outlier’s pull on the fit is capped at d instead of growing with its size, while small errors are still treated as in least squares.
Use it when regression data has occasional gross errors (sensor glitches, data-entry mistakes) that should not dominate. d is in residual units: about 1.35 noise standard deviations is the classic choice. A large d behaves like squaredLoss, a small one like a scaled absolute loss (median regression). No gradient function is provided; with respect to the residuals it is \(\operatorname{clip}(r, -d, d)/n\).
robust-regressionregression
def huberLoss(r, d) ≔ einsum("i->", (abs(r) ≤ d) · r · r / 2 + (abs(r) > d) · d · (abs(r) − d / 2)) / count(r)
pinballLoss
pinballLoss(r, q) · r residuals n · q quantile level in (0, 1) scalar → loss scalar
Let \(r \in \mathbb{R}^n\) be the residuals and \(0 < q < 1\). Then
\[\operatorname{pinballLoss}(r, q) = \frac{1}{n} \sum_{i=1}^{n} \max\big(q\, r_i,\, (q - 1)\, r_i\big) .\]The mean pinball (quantile) loss, \(\frac{1}{n}\sum_i \max(q r_i, (q - 1) r_i)\). Minimising it gives a quantile regression at level q; with residuals taken as target minus prediction, under-prediction is charged \(q\) and over-prediction \(1 - q\) per unit.
Intuition A tilted absolute loss: under-prediction costs q per unit and over-prediction \(1 - q\). The prediction that minimises it has a fraction q of the targets below it, which is the q-quantile; at q = 0.5 it is half the absolute error and gives the median.
Use it when you need a quantile rather than a mean: prediction intervals (fit q = 0.05 and q = 0.95), stock levels that meet 90 % of demand, latency percentiles, or any cost that is asymmetric. Residuals must be target minus prediction for this reading. It is also the right score for evaluating quantile forecasts.
quantile-forecastingprediction-intervalsuncertainty-quantification
def pinballLoss(r, q) ≔ einsum("i->", max(q · r, (q − 1) · r)) / count(r)
crossEntropy
crossEntropy(z, g) · z logits n×K · g one-hot labels n×K → loss scalar
Let \(Z \in \mathbb{R}^{n \times K}\) (the argument z) be the logits and \(G \in \mathbb{R}^{n \times K}\) (the argument g) the labels, and let \(\mu_i = \max_k Z_{ik}\). Then
where \(\mu_i + \ln \sum_{k'} e^{Z_{ik'} - \mu_i} = \ln \sum_{k'} e^{Z_{ik'}}\). When every row of \(G\) sums to 1 (for example one-hot labels) this is the formula shown below; otherwise each row’s log-normaliser is weighted by that row’s sum.
The mean multiclass cross-entropy, computed from logits rather than probabilities. The log-normaliser is a stable log-sum-exp, so the loss never takes the logarithm of a softmax that has underflowed.
\[L = \frac{1}{n} \sum_i \Big( \log \sum_k e^{z_{ik}} - \sum_k g_{ik} z_{ik} \Big)\]Intuition Multiclass log-loss: the average of \(-\ln\) of the softmax probability given to the true class, computed from logits with a log-sum-exp so it never takes the log of an underflowed zero. Guessing uniformly over \(K\) classes costs \(\ln K\).
Use it when you train or evaluate any \(K\)-class probabilistic classifier (softmax regression, a network’s output layer). Pass raw logits, not probabilities. Soft labels whose rows sum to 1 (label smoothing, distillation) work too. The gradient with respect to the logits is \((\operatorname{softmax}(Z) - G)/n\), which this module does not define; softmaxRows in Prelude.Kernels gives the softmax. For two classes logisticLoss is the same objective.
classificationmodel-evaluationtext-classification
def crossEntropy(z, g) ≔ (einsum("ik->", g · lifted(z, logSumRows(z))) − einsum("ik,ik->", g, z)) / rows(z)
The steps
Each step turns the current parameters w and a gradient g into new parameters. Optimisers that carry state (a velocity, a running square, Adam’s two moments) have a separate function to update that state; the calling loop carries it from turn to turn.
descend
descend(w, g, r) · w parameters any shape · g gradient same shape · r learning rate scalar → parameters same shape
Let \(w\) and \(g\) be tensors of the same shape, and let \(r > 0\). Then
\[\operatorname{descend}(w, g, r) = w - r\,g .\]One plain gradient-descent step, \(w - r g\).
Intuition Walk downhill: move every parameter against its gradient by a step proportional to r. Too small a rate crawls; too large overshoots and diverges (on a quadratic, beyond 2 divided by the largest curvature).
Use it when the objective is smooth and well conditioned, or you want the simplest baseline for checking a gradient. Standardising features helps it more than anything else. When progress zigzags down a long narrow valley, use momentum with heavyStep; when gradients differ greatly in scale between parameters, use adamStep or rmsStep.
gradient-based-trainingparameter-estimation
def descend(w, g, r) ≔ w − r · g
momentum
momentum(v, g, b) · v velocity any shape · g gradient same shape · b momentum coefficient scalar → velocity same shape
Let \(v\) and \(g\) be tensors of the same shape, and let \(b \in \mathbb{R}\). Then
\[\operatorname{momentum}(v, g, b) = b\,v + g .\]Updates the heavy-ball velocity, \(v \leftarrow b v + g\). Start v at zero.
Intuition The velocity is a decaying sum of past gradients: directions that stay consistent add up and speed travel along a valley, while components that flip sign each turn cancel out. b acts like low friction; b = 0.9 remembers roughly the last 10 gradients.
Use it when plain descend zigzags on an ill-conditioned problem or crawls over a plateau. Update v every turn, then take heavyStep. b = 0.9 is typical. This form has no \((1 - b)\) factor, so on a steady slope the step is \(1/(1-b)\) times larger than descend would take: divide r by about that factor (10 at b = 0.9) when switching.
gradient-based-trainingslow-convergence
def momentum(v, g, b) ≔ b · v + g
heavyStep
heavyStep(w, v, r) · w parameters any shape · v velocity from momentum same shape · r learning rate scalar → parameters same shape
Let \(w\) and \(v\) be tensors of the same shape, and let \(r > 0\). Then
\[\operatorname{heavyStep}(w, v, r) = w - r\,v .\]The momentum update of the parameters, \(w - r v\).
Intuition The parameter move of the heavy-ball method: step along the accumulated velocity rather than the latest gradient, like a ball that keeps rolling through small bumps.
Use it when paired with momentum, once per turn: first \(v \leftarrow \operatorname{momentum}(v, g, b)\), then \(w \leftarrow \operatorname{heavyStep}(w, v, r)\). With b = 0 the pair is exactly descend. If the loss oscillates or climbs, lower r or b.
gradient-based-trainingslow-convergence
def heavyStep(w, v, r) ≔ w − r · v
rmsCache
rmsCache(s, g, b) · s running mean square any shape · g gradient same shape · b decay scalar → running mean square same shape
Let \(s\) and \(g\) be tensors of the same shape, and let \(0 \le b < 1\). Then
\[\operatorname{rmsCache}(s, g, b) = b\,s + (1 - b)\, g \odot g .\]The RMSProp cache, \(s \leftarrow b s + (1 - b) g^2\). Start s at zero; a typical decay is 0.9.
Intuition An exponentially weighted average of each parameter’s squared gradient: a running estimate of how large that coordinate’s gradients usually are. With decay b it averages roughly the last \(1/(1 - b)\) turns.
Use it when you run RMSProp: update it every turn, before rmsStep. b = 0.9 is typical; a larger b gives a smoother but slower-reacting scale. Starting from zero makes early values too small and there is no bias correction here (unlike adamStep), so the first steps are several times r; a smaller rate for the first turns helps.
gradient-based-trainingpoorly-scaled-gradients
def rmsCache(s, g, b) ≔ b · s + (1 − b) · g · g
rmsStep
rmsStep(w, g, s, r) · w parameters any shape · g gradient same shape · s cache from rmsCache same shape · r learning rate scalar → parameters same shape
Let \(w\), \(g\) and \(s\) be tensors of the same shape with entries indexed by \(\iota\), and let \(r > 0\). Then
\[\operatorname{rmsStep}(w, g, s, r)_\iota = w_\iota - \frac{r\, g_\iota}{\sqrt{\max(s_\iota, 0)} + 10^{-8}} .\]The RMSProp update, each coordinate scaled by the root of its running mean square:
\[w \leftarrow w - \frac{r\, g}{\sqrt{s} + 10^{-8}}\]Intuition Divides each coordinate’s gradient by its typical size, so every parameter moves by roughly r per turn whatever its gradient scale: parameters with small gradients speed up and those with large, noisy ones slow down.
Use it when gradient scales differ widely between parameters (mixed feature scales, recurrent networks, objectives that drift during training). A rate around \(10^{-3}\) is a common start. The \(10^{-8}\) keeps a zero cache from dividing by zero. adamStep adds momentum and bias correction and is usually the safer default.
gradient-based-trainingpoorly-scaled-gradientsneural-network-training
def rmsStep(w, g, s, r) ≔ w − r · g / (sqrt(max(s, 0)) + 1e-8)
adamMean
adamMean(m, g, b) · m first moment any shape · g gradient same shape · b decay β₁ scalar → first moment same shape
Let \(m\) and \(g\) be tensors of the same shape, and let \(0 \le b < 1\). Then
\[\operatorname{adamMean}(m, g, b) = b\,m + (1 - b)\,g .\]Adam’s running mean of the gradient, \(m \leftarrow b m + (1 - b) g\). Start at zero; β₁ is typically 0.9.
Intuition Momentum written as an average rather than a sum: a smoothed estimate of the gradient’s direction that stays on the gradient’s own scale. \(\beta_1 = 0.9\) averages roughly the last 10 turns, which filters out mini-batch noise.
Use it when you run Adam: update it every turn together with adamVar, and pass the raw value to adamStep, which corrects the bias from starting at zero. That correction assumes the start really is zero, so do not initialise it with a first gradient.
gradient-based-trainingnoisy-gradients
def adamMean(m, g, b) ≔ b · m + (1 − b) · g
adamVar
adamVar(v, g, b) · v second moment any shape · g gradient same shape · b decay β₂ scalar → second moment same shape
Let \(v\) and \(g\) be tensors of the same shape, and let \(0 \le b < 1\). Then
\[\operatorname{adamVar}(v, g, b) = b\,v + (1 - b)\, g \odot g .\]Adam’s running mean of the squared gradient, \(v \leftarrow b v + (1 - b) g^2\). Start at zero; β₂ is typically 0.999.
Intuition A slow running average of each coordinate’s squared gradient: the scale Adam divides by. \(\beta_2 = 0.999\) averages roughly the last 1000 turns, so the scale stays stable even when single gradients are noisy.
Use it when you run Adam alongside adamMean: start at zero and pass the raw value to adamStep. A lower \(\beta_2\) (for example 0.99) lets the scale adapt faster when gradient magnitudes change during training, at the cost of noisier steps.
gradient-based-trainingpoorly-scaled-gradients
def adamVar(v, g, b) ≔ b · v + (1 − b) · g · g
adamStep
adamStep(w, m, v, b1, b2, t, r) · w parameters any shape · m first moment same shape · v second moment same shape · b1 β₁ scalar · b2 β₂ scalar · t turn number, starting at 1 scalar · r learning rate scalar → parameters same shape
Let \(w\), \(m\) and \(v\) be tensors of the same shape with entries indexed by \(\iota\), let \(0 \le \beta_1, \beta_2 < 1\) (the arguments b1, b2), let \(t \in \{1, 2, \dots\}\) and \(r > 0\). Then
The Adam update. Because both moments start at zero they are biased towards zero early on; the step divides them by \(1 - \beta^t\) to undo that, so pass the raw moments and the current turn number:
\[w \leftarrow w - \frac{r\, \hat m}{\sqrt{\hat v} + 10^{-8}}, \qquad \hat m = \frac{m}{1 - \beta_1^t}, \quad \hat v = \frac{v}{1 - \beta_2^t}\]Intuition Momentum supplies the direction and RMSProp the scale, so each coordinate moves by about r along its smoothed gradient whatever the raw gradient size. The bias correction keeps the first turns from being distorted by moments that started at zero.
Use it when the problem is noisy, stochastic or badly scaled, neural networks included; it is the usual default. \(r = 10^{-3}\), \(\beta_1 = 0.9\), \(\beta_2 = 0.999\) are the standard values. t counts turns from 1: at t = 0 the correction divides by zero. On a small, smooth, deterministic problem a Newton fit such as logisticFit in Prelude.Regression converges in far fewer turns.
gradient-based-trainingnoisy-gradientsneural-network-training
def adamStep(w, m, v, b1, b2, t, r) ≔ w − r · adamBias(m, b1, t) / (sqrt(max(adamBias(v, b2, t), 0)) + 1e-8)
soften
soften(z, t) · z values any shape · t threshold scalar → values same shape
Let \(z\) be a tensor of any shape with entries \(z_\iota\), and let \(t \ge 0\). Then
\[\operatorname{soften}(z, t)_\iota = \operatorname{sign}(z_\iota)\, \max(\lvert z_\iota \rvert - t, 0) .\]The soft-thresholding operator, \(\operatorname{sign}(z)\max(\lvert z \rvert - t, 0)\). Values smaller than t in size become exactly zero, which is what makes a lasso fit sparse.
Intuition Moves every value toward zero by t and stops at zero: values smaller than t in size vanish and larger ones shrink by exactly t. It is the closest point to z once a price of t per unit of absolute size is charged, which is why the L1 penalty yields exact zeros.
Use it when you want sparsity: inside L1-penalised solvers (istaStep), denoising by thresholding small wavelet or Fourier coefficients, or pruning small weights. If survivors should keep their full size, use a hard threshold (a mask) instead, since this shrinks them all by t.
feature-selectiondenoisingmodel-compression
def soften(z, t) ≔ sign(z) · max(abs(z) − t, 0)
istaStep
istaStep(w, g, r, l) · w parameters any shape · g gradient of the smooth part same shape · r learning rate scalar · l L1 penalty scalar → parameters same shape
Let \(w\) and \(g\) be tensors of the same shape, and let \(r > 0\) and \(l \ge 0\). Then
\[\operatorname{istaStep}(w, g, r, l) = \operatorname{soften}(w - r\,g,\; r\,l), \qquad \operatorname{istaStep}(w, g, r, l)_\iota = \operatorname{sign}(w_\iota - r g_\iota)\max(\lvert w_\iota - r g_\iota \rvert - r l, 0) .\]One step of ISTA (iterative soft thresholding): a gradient step on the smooth part of the objective followed by soft thresholding at \(r l\). Used by the lasso and elastic-net fits in Prelude.Regression.
Intuition Splits the objective in two: the smooth part is handled by an ordinary gradient step, and the L1 penalty is handled exactly by soft thresholding at \(rl\). Repeated, it converges to the penalised minimum and produces exact zeros along the way.
Use it when you write your own L1-penalised fit for a loss other than squared error, such as sparse logistic regression: pass that loss’s gradient as g. The rate must stay below \(2/L\), where \(L\) bounds the curvature of the smooth part; \(1/L\) is the classic safe choice. For squared error, lassoFit and elasticFit in Prelude.Regression already run this loop.
feature-selectionregularisationgradient-based-training
def istaStep(w, g, r, l) ≔ soften(w − r · g, r · l)
Building blocks
adamBias(m, b, t)— Let \(m\) be a tensor of any shape, \(0 \le b < 1\) and \(t \in \{1, 2, \dots\}\). Then \(\operatorname{adamBias}(m, b, t) = m / (1 - b^t)\). The bias correction used insideadamStep.
Examples
Apply activations to a vector
open Prelude.Learning
fn squash(z) ↦ sigmoid(z)
fn smoothRamp(z) ↦ softplus(z)
fn ramp(z) ↦ relu(z)
using Examples;
var s = Activations.Squash(new[] { -2d, 0d, 2d });
var p = Activations.SmoothRamp(new[] { -800d, 0d, 800d }); // no overflow at either tail
var r = Activations.Ramp(new[] { -2d, 0d, 2d });
Output (verified, six decimals):
sigmoid = [0.119203, 0.5, 0.880797]
softplus = [0, 0.693147, 800]
relu = [0, 0, 2]
Score predictions with four losses
For predictions [0.9, 0.2, 0.6] against labels [1, 0, 1] the squared loss is \((0.01 + 0.04 + 0.16)/6 = 0.035\) and the logistic loss is \(-(\ln 0.9 + \ln 0.8 + \ln 0.6)/3 = 0.279777\). The hinge loss of margins [2, −0.5, 0.3] against labels [1, −1, 1] is \((0 + 0.5 + 0.7)/3 = 0.4\). The cross-entropy averages \(\ln(e^2 + 2) - 2\) for the first row and \(\ln 3\) for the second.
open Prelude.Learning
fn meanSquared(p, y) ↦ squaredLoss(p, y)
fn logLoss(p, y) ↦ logisticLoss(p, y)
fn marginLoss(m, y) ↦ hingeLoss(m, y)
fn softmaxLoss(z, g) ↦ crossEntropy(z, g)
using Examples;
var p = new[] { 0.9, 0.2, 0.6 };
var y = new[] { 1d, 0d, 1d };
var squared = Losses.MeanSquared(p, y);
var logistic = Losses.LogLoss(p, y);
var hinge = Losses.MarginLoss(new[] { 2d, -0.5, 0.3 }, new[] { 1d, -1d, 1d });
var z = new double[,] { { 2, 0, 0 }, { 0, 0, 0 } }; // logits
var g = new double[,] { { 1, 0, 0 }, { 0, 1, 0 } }; // one-hot labels
var cross = Losses.SoftmaxLoss(z, g);
Output (verified, six decimals):
squared = 0.035
logistic = 0.279777
hinge = 0.4
cross = 0.669079
Soft-threshold coefficients and take one Adam step
Soft thresholding at 1 moves every value one unit towards zero and stops at zero. On the first Adam turn the bias-corrected moments are exactly \(g\) and \(g^2\), so each coordinate moves by the learning rate against the sign of its gradient.
open Prelude.Learning
fn shrink(z, t) ↦ soften(z, t)
fn adamTurn(w, g, b1, b2, t, r) ↦ adamStep(w, adamMean(0 · w, g, b1), adamVar(0 · w, g, b2), b1, b2, t, r)
using Examples;
var shrunk = Steps.Shrink(new[] { 3d, 0.5, 1.2, -2d }, 1.0);
var next = Steps.AdamTurn(new[] { 1d, 1d }, new[] { 0.5, -4d }, 0.9, 0.999, 1.0, 0.1);
Output (verified):
soften = [2, 0, 0.2, -1]
adam = [0.9, 1.1]
Write a gradient-descent loop
A training loop is a definition that names itself. This one fits the slope of \(y = w x\) by descending the squared loss; the gradient with respect to \(w\) is dSquared contracted with x. For x = [1, 2, 3] and y = [2, 4, 6] the slope is 2.
open Prelude.Learning
def slopeLoop(t, r, x, y, w) ≔ match t { when 0 ↦ w other ↦ slopeLoop(t − 1, r, x, y, descend(w, einsum("i,i->", x, dSquared(w · x, y)), r)) }
fn fitSlope(t, r, x, y) ↦ slopeLoop(t, r, x, y, 0)
using MLambda.Hilbert.Runtime;
using Examples;
var slope = Descent.FitSlope(20.0, 0.2, (Tensor)new[] { 1d, 2d, 3d }, (Tensor)new[] { 2d, 4d, 6d });
Output (verified):
2