Regression
Logistic and Poisson regression by Newton's method, the lasso and elastic net, a linear SVM, kernel ridge and Gaussian-process regression.
This module fits the linear models that ordinary least squares does not cover: logistic regression, Poisson regression, the lasso and the elastic net, a primal linear support-vector machine, kernel ridge regression and its probabilistic twin, Gaussian-process regression.
Logistic and Poisson regression are fitted by iteratively reweighted least squares (IRLS), which is Newton’s method on the log-likelihood: each turn solves a weighted least-squares problem with solve from Prelude.LinearAlgebra. There is no learning rate to tune, and about six turns reach convergence on ordinary data. The lasso, the elastic net and the SVM have a corner in their objective, so Newton’s method does not apply; they are fitted by (sub)gradient methods and take a learning rate r that you choose.
Conventions. The design matrix x is n×p, one row per observation; no intercept is added for you — include a column of ones if you want one (in the lasso, elastic net and SVM that column is penalised like any other). Responses y are a vector n: 0/1 labels for logistic regression, non-negative counts for Poisson regression, −1/+1 labels for the SVM. Every iterative fit takes the number of turns t as its first argument and starts from zero coefficients. Coefficients come back as a vector p.
Honest limitations. On perfectly separable data the logistic likelihood has no maximum; a ridge of \(10^{-8}\) on the Hessian turns the divergence into large but finite coefficients. The lasso and elastic net diverge if the rate exceeds \(2/\lambda_{\max}(X^\top X/n)\). The Gaussian-process variance assumes an RBF kernel (whose diagonal is one) and forms a matrix inverse by Newton–Schulz iteration, so it is less accurate than the mean when the kernel matrix is badly conditioned.
open Prelude.Regression
Logistic regression, by iteratively reweighted least squares
The model is \(P(y_i = 1) = \sigma(x_i^\top b)\). A Newton step on the log-likelihood is a weighted least-squares solve with weights \(w_i = p_i(1 - p_i)\) and working response \(z_i = x_i^\top b + (y_i - p_i)/w_i\):
\[b \leftarrow (X^\top W X + 10^{-8} I)^{-1} X^\top W z\]Weights are floored at \(10^{-10}\) so that a row the fit is already certain about cannot remove a column from the Hessian.
Throughout this page, \(\operatorname{solve}(A, v)\) is the Jacobi-preconditioned conjugate-gradient solution of \(A\beta = v\) from Prelude.LinearAlgebra (equal to \(A^{-1}v\) in exact arithmetic for symmetric positive-definite \(A\)), and \(\sigma(s) = 1/(1 + e^{-s})\) is the logistic function.
logisticOdds
logisticOdds(b, x) · b coefficients p · x design matrix n×p → log-odds n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), and let \(b \in \mathbb{R}^p\) be the coefficients. Then
The linear predictor \(Xb\), which is the log-odds of each row.
Intuition The raw, unbounded score of a logistic model. Each unit of feature \(j\) adds \(b_j\) to the log-odds \(\ln\frac{p}{1-p}\), so a coefficient is an additive effect on a log scale and \(e^{b_j}\) is its odds ratio. A score of 0 means even odds.
Use it when you need the decision value rather than a probability: ranking rows, thresholding at 0 (the same as probability 0.5), or reporting odds ratios. It keeps its resolution in the tails, where logisticProbability rounds to 0 or 1. It computes the same thing as svmMargin; only the coefficients you pass differ.
classificationrisk-scoringmodel-interpretation
def logisticOdds(b, x) ≔ apply(x, b)
logisticProbability
logisticProbability(b, x) · b coefficients p · x design matrix n×p → probabilities n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), and let \(b \in \mathbb{R}^p\). Then, for \(i = 1, \dots, n\),
The fitted probability that each row’s label is 1, \(\sigma(Xb)\). Pass the coefficients from logisticFit and the rows to score.
Intuition The sigmoid squeezes the log-odds onto \((0, 1)\): a score of 0 is 50 %, a score of \(\pm 2.2\) is about 90 % or 10 %, and each further unit multiplies the odds by \(e\). Maximum likelihood fits exactly this quantity, so when the model is well specified the outputs are calibrated probabilities.
Use it when you need a probability rather than a label: churn or default risk, expected-cost decisions with a threshold other than 0.5, or input to logisticLoss. Score rows with the same column layout, intercept column included, that logisticFit was given. For a hard yes/no alone, comparing logisticOdds with 0 is enough.
probability-calibrationclassificationrisk-scoring
def logisticProbability(b, x) ≔ sigmoid(logisticOdds(b, x))
logisticFit
logisticFit(t, x, y) · t number of Newton turns scalar · x design matrix n×p · y 0/1 labels n → coefficients p
Let \(t \in \mathbb{N}\), let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), and let \(y \in \{0, 1\}^n\) be the labels. For \(b \in \mathbb{R}^p\) write \(p_i(b) = \sigma(x_i^{\top} b)\), \(w_i(b) = \max(p_i(b)(1 - p_i(b)), 10^{-10})\), \(W(b) = \operatorname{diag}(w(b))\) and \(z_i(b) = x_i^{\top} b + (y_i - p_i(b))/w_i(b)\). Let \(b^{(0)} = 0 \in \mathbb{R}^p\) and
Then \(\operatorname{logisticFit}(t, X, y) = b^{(t)}\).
Maximum-likelihood logistic regression by t IRLS turns from zero. The first turn is an ordinary least-squares fit to \(4(y - \tfrac12)\); about six turns give a converged fit, and extra turns leave it unchanged.
Intuition Each Newton turn replaces the log-likelihood by a parabola and jumps to its peak, and that jump is a weighted least-squares fit. Rows the model is unsure about (\(p \approx 0.5\)) get the most weight; rows it is already sure about get almost none. Quadratic convergence is why a handful of turns suffices and no learning rate exists.
Use it when the outcome is binary (click, default, diagnosis) and you want interpretable coefficients plus probabilities. Each turn solves a p×p system, so it suits tens to hundreds of features. About six turns converge. If the coefficients grow huge, the classes are perfectly separable and the numbers mean little; the \(10^{-8}\) ridge only keeps them finite, so consider svmFit, whose L2 penalty gives a well-defined boundary.
classificationparameter-estimationprobability-calibration
def logisticFit(t, x, y) ≔ irlsLoop(t, x, y, zeroed(moment(x, y)))
Building blocks
irlsWeights(b, x)— Let \(b \in \mathbb{R}^p\) and \(X \in \mathbb{R}^{n \times p}\) with rows \(x_i^{\top}\), and \(p_i = \sigma(x_i^{\top} b)\). Then \(\operatorname{irlsWeights}(b, X)_i = w_i = \max(p_i(1 - p_i), 10^{-10})\). The IRLS weights p(1 − p), floored at 1e-10.irlsResponse(b, x, y)— Let \(b\), \(X\), \(p_i\) and \(w_i\) be as inirlsWeights, and \(y \in \{0,1\}^n\). Then \(\operatorname{irlsResponse}(b, X, y)_i = z_i = x_i^{\top} b + (y_i - p_i)/w_i\). The working response: log-odds plus residual divided by weight.irlsGram(b, x)— Let \(b\), \(X\) and \(w\) be as inirlsWeights, and \(W = \operatorname{diag}(w)\). Then \(\operatorname{irlsGram}(b, X) = X^{\top} W X + 10^{-8} I_p\), with entries \(\sum_{i} w_i X_{ij} X_{ik} + 10^{-8}[j = k]\). The weighted Gram matrix plus a 1e-8 ridge.irlsMoment(b, x, y)— Let \(b\), \(X\), \(w\), \(W\) and \(z\) be as above. Then \(\operatorname{irlsMoment}(b, X, y) = X^{\top} W z = \sum_{i=1}^{n} w_i z_i x_i\). The weighted moment.irlsStep(b, x, y)— Let \(b\), \(X\), \(y\) be as above. Then \(\operatorname{irlsStep}(b, X, y) = \operatorname{solve}(\operatorname{irlsGram}(b, X), \operatorname{irlsMoment}(b, X, y))\). One Newton step:solveof the weighted normal equations.irlsLoop(t, x, y, b)— Let \(t \in \mathbb{N}\), \(X\), \(y\) as above and \(b \in \mathbb{R}^p\). Then \(\operatorname{irlsLoop}(0, X, y, b) = b\) and \(\operatorname{irlsLoop}(t, X, y, b) = \operatorname{irlsLoop}(t - 1, X, y, \operatorname{irlsStep}(b, X, y))\). RepeatsirlsStepttimes.
The lasso, by iterative soft thresholding
The lasso minimises squared error plus an L1 penalty, and the elastic net adds an L2 penalty as well:
\[\min_b \; \frac{1}{2n}\lVert y - Xb \rVert^2 + l\,\lVert b \rVert_1 + \frac{a}{2}\lVert b \rVert^2\]with \(a = 0\) for the lasso. Each turn takes a gradient step on the smooth part and then soft-thresholds at \(r l\) (see istaStep in Prelude.Learning). A coefficient that the penalty outweighs becomes exactly zero, so the lasso selects variables where ridge regression only shrinks them. The L1 term is not differentiable at zero, which is why these fits need a learning rate r; keep it below \(2/\lambda_{\max}(X^\top X/n)\) (with the elastic net, below \(2/(\lambda_{\max} + a)\)).
lassoFit
lassoFit(t, r, l, x, y) · t number of turns scalar · r learning rate scalar · l L1 penalty scalar · x design matrix n×p · y responses n → coefficients p
Let \(t \in \mathbb{N}\), let \(r > 0\) be the learning rate, \(l \ge 0\) the L1 penalty, \(X \in \mathbb{R}^{n \times p}\) (the argument x) the design matrix and \(y \in \mathbb{R}^n\) the responses. Let \(S_\tau(v)_j = \operatorname{sign}(v_j)\max(\lvert v_j \rvert - \tau, 0)\) be soft thresholding. Let \(b^{(0)} = 0 \in \mathbb{R}^p\) and
Then \(\operatorname{lassoFit}(t, r, l, X, y) = b^{(t)}\).
The lasso coefficients after t ISTA turns from zero.
Intuition The L1 penalty charges every coefficient a flat price per unit of size, so a feature whose contribution to the fit is worth less than that price is set to exactly zero instead of being kept small. Each turn is a gradient step on the squared error, then every coefficient is pulled toward zero by \(rl\) and stopped at zero.
Use it when you have many candidate features and want a sparse, readable model, including when p is close to or above n. Standardise the columns first, since the penalty is in each coefficient’s units, and raise l to drop more variables. With standardised columns \(r = 1/p\) is always below the divergence bound. Budget far more turns than the Newton fits and check the result has stopped changing. With correlated features the lasso picks among them arbitrarily; use elasticFit.
feature-selectionregularisationoverfitting-controlregression
def lassoFit(t, r, l, x, y) ≔ istaLoop(t, r, l, x, y, zeroed(moment(x, y)))
elasticFit
elasticFit(t, r, l, a, x, y) · t number of turns scalar · r learning rate scalar · l L1 penalty scalar · a L2 penalty scalar · x design matrix n×p · y responses n → coefficients p
Let \(t\), \(r > 0\), \(l \ge 0\), \(X\), \(y\) and \(S_\tau\) be as in lassoFit, and let \(a \ge 0\) be the L2 penalty. Let \(b^{(0)} = 0 \in \mathbb{R}^p\) and
Then \(\operatorname{elasticFit}(t, r, l, a, X, y) = b^{(t)}\).
The elastic-net coefficients after t turns from zero: the L2 penalty enters the gradient, the L1 penalty enters the threshold.
Intuition Two penalties with different jobs. The L1 part, applied as the threshold, still zeroes weak coefficients. The L2 part, added to the gradient, shrinks all of them smoothly and makes the objective strictly convex, so a group of correlated features tends to be kept or dropped together instead of one member being picked at random.
Use it when features come in correlated groups (genes, neighbouring sensors, overlapping survey questions) or the lasso’s selection changes from sample to sample. With a = 0 it is exactly lassoFit; with l = 0 it is ridge regression solved by gradient descent. The standardisation advice of lassoFit applies, with the rate bound \(2/(\lambda_{\max} + a)\).
feature-selectioncorrelated-featuresregularisationoverfitting-control
def elasticFit(t, r, l, a, x, y) ≔ elasticLoop(t, r, l, a, x, y, zeroed(moment(x, y)))
Building blocks
leastGradient(b, x, y)— Let \(b \in \mathbb{R}^p\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^n\). Then \(\operatorname{leastGradient}(b, X, y) = \frac{1}{n}(X^{\top} X b - X^{\top} y)\). The gradient of the averaged squared error \(\frac{1}{2n}\lVert y - Xb \rVert^2\).istaTurn(r, l, x, y, b)— Let \(r > 0\), \(l \ge 0\), \(X\), \(y\), \(b\) as above and \(S_\tau\) as inlassoFit. Then \(\operatorname{istaTurn}(r, l, X, y, b) = S_{rl}(b - r\,\operatorname{leastGradient}(b, X, y))\). One lasso turn:istaSteponleastGradient.istaLoop(t, r, l, x, y, b)— Let \(t \in \mathbb{N}\) and the rest as above. Then \(\operatorname{istaLoop}(0, r, l, X, y, b) = b\) and \(\operatorname{istaLoop}(t, r, l, X, y, b) = \operatorname{istaLoop}(t - 1, r, l, X, y, \operatorname{istaTurn}(r, l, X, y, b))\). RepeatsistaTurnttimes.elasticGradient(a, b, x, y)— Let \(a \ge 0\) and \(b\), \(X\), \(y\) as above. Then \(\operatorname{elasticGradient}(a, b, X, y) = \frac{1}{n}(X^{\top} X b - X^{\top} y) + a\,b\).leastGradientplus the ridge term.elasticTurn(r, l, a, x, y, b)— Let \(r > 0\), \(l \ge 0\), \(a \ge 0\) and \(X\), \(y\), \(b\) as above. Then \(\operatorname{elasticTurn}(r, l, a, X, y, b) = S_{rl}(b - r\,\operatorname{elasticGradient}(a, b, X, y))\). One elastic-net turn:istaSteponelasticGradient.elasticLoop(t, r, l, a, x, y, b)— Let \(t \in \mathbb{N}\) and the rest as above. Then \(\operatorname{elasticLoop}(0, r, l, a, X, y, b) = b\) and \(\operatorname{elasticLoop}(t, r, l, a, X, y, b) = \operatorname{elasticLoop}(t - 1, r, l, a, X, y, \operatorname{elasticTurn}(r, l, a, X, y, b))\). RepeatselasticTurnttimes.
A Poisson GLM at the log link
Counts are modelled as \(y_i \sim \operatorname{Poisson}(\mu_i)\) with \(\ln \mu_i = x_i^\top b\). The fit is the same IRLS as logistic regression with a different weight and working response: at the canonical log link the weight is the rate itself, \(w_i = \max(\mu_i, 10^{-10})\), and \(z_i = x_i^\top b + (y_i - \mu_i)/w_i\). The linear predictor is clamped to \([-30, 30]\) before exponentiating, which covers rates from about \(10^{-13}\) to \(10^{13}\) and keeps a wild turn from overflowing (the working response uses the unclamped \(x_i^\top b\)).
poissonRate
poissonRate(b, x) · b coefficients p · x design matrix n×p → rates n
Let \(b \in \mathbb{R}^p\) be the coefficients and \(X \in \mathbb{R}^{n \times p}\) (the argument x) the design matrix with rows \(x_i^{\top}\). Then, for \(i = 1, \dots, n\),
The fitted mean count of each row, \(\exp(Xb)\) with the exponent clamped to \([-30, 30]\).
Intuition The log link makes effects multiplicative: raising feature \(j\) by one unit multiplies the expected count by \(e^{b_j}\), and the exponential keeps every rate positive whatever the linear predictor is. The clamp at \(\pm 30\) only bites on absurd predictors.
Use it when you want expected counts from fitted poissonFit coefficients: calls per hour, defects per batch, claims per policy. The output is a mean, not a whole number; the chance of an exact count follows from the Poisson distribution with that mean. When the outcome is yes/no rather than how many, use logisticProbability.
count-datademand-forecastingregression
def poissonRate(b, x) ≔ exp(min(max(apply(x, b), 0 − 30), 30))
poissonFit
poissonFit(t, x, y) · t number of Newton turns scalar · x design matrix n×p · y counts n → coefficients p
Let \(t \in \mathbb{N}\), let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), and let \(y \in \mathbb{R}^n\), \(y_i \ge 0\), be the counts. For \(b \in \mathbb{R}^p\) write \(\mu_i(b) = \operatorname{poissonRate}(b, X)_i\), \(w_i(b) = \max(\mu_i(b), 10^{-10})\), \(W(b) = \operatorname{diag}(w(b))\) and \(z_i(b) = x_i^{\top} b + (y_i - \mu_i(b))/w_i(b)\). Let \(b^{(0)} = 0 \in \mathbb{R}^p\) and
Then \(\operatorname{poissonFit}(t, X, y) = b^{(t)}\).
Maximum-likelihood Poisson regression by t IRLS turns from zero.
Intuition The same Newton machinery as logisticFit, tuned for counts. A Poisson count’s variance equals its mean, so each row is weighted by its predicted rate: a row expecting 100 events pins the coefficients down far more than one expecting 0.1.
Use it when the response is a non-negative count (events, visits, insurance claims) and effects plausibly multiply. It assumes variance equals mean; with much more spread (overdispersion) the coefficients are still sensible but look more certain than they are. Counts observed over unequal exposure (different durations or populations) need an offset, which this fit does not take: adding \(\ln(\text{exposure})\) as a column estimates its coefficient instead of fixing it at 1.
count-datademand-forecastingparameter-estimation
def poissonFit(t, x, y) ≔ poissonLoop(t, x, y, zeroed(moment(x, y)))
Building blocks
poissonWeights(b, x)— Let \(b \in \mathbb{R}^p\), \(X \in \mathbb{R}^{n \times p}\) and \(\mu_i = \operatorname{poissonRate}(b, X)_i\). Then \(\operatorname{poissonWeights}(b, X)_i = w_i = \max(\mu_i, 10^{-10})\). The IRLS weights: the rate, floored at 1e-10.poissonResponse(b, x, y)— Let \(b\), \(X\), \(\mu_i\) and \(w_i\) be as inpoissonWeights, and \(y \in \mathbb{R}^n\). Then \(\operatorname{poissonResponse}(b, X, y)_i = z_i = x_i^{\top} b + (y_i - \mu_i)/w_i\). The working response: linear predictor plus residual divided by the rate.poissonGram(b, x)— Let \(b\), \(X\) and \(w\) be as above, and \(W = \operatorname{diag}(w)\). Then \(\operatorname{poissonGram}(b, X) = X^{\top} W X + 10^{-8} I_p\). The weighted Gram matrix plus a 1e-8 ridge.poissonStep(b, x, y)— Let \(b\), \(X\), \(y\), \(W\) and \(z\) be as above. Then \(\operatorname{poissonStep}(b, X, y) = \operatorname{solve}(X^{\top} W X + 10^{-8} I_p,\; X^{\top} W z)\). One Newton step:solveof the weighted normal equations.poissonLoop(t, x, y, b)— Let \(t \in \mathbb{N}\), \(X\), \(y\) as above and \(b \in \mathbb{R}^p\). Then \(\operatorname{poissonLoop}(0, X, y, b) = b\) and \(\operatorname{poissonLoop}(t, X, y, b) = \operatorname{poissonLoop}(t - 1, X, y, \operatorname{poissonStep}(b, X, y))\). RepeatspoissonStepttimes.
The primal support-vector machine
A linear SVM minimises the mean hinge loss plus an L2 penalty:
\[\min_b \; \frac{1}{n}\sum_i \max(1 - y_i\, x_i^\top b,\, 0) + \frac{l}{2}\lVert b \rVert^2\]The hinge has a corner, so the fit is a subgradient descent with a fixed learning rate rather than a Newton method. Labels must be −1 and +1, not 0 and 1.
svmMargin
svmMargin(b, x) · b coefficients p · x design matrix n×p → margins n
Let \(b \in \mathbb{R}^p\) be the coefficients and \(X \in \mathbb{R}^{n \times p}\) (the argument x) the design matrix with rows \(x_i^{\top}\). Then
The signed score \(Xb\) of each row; its sign is the predicted label.
Intuition A signed score: its sign is the predicted side, its size how far past the boundary the row sits, in units where the margin edges are at \(\pm 1\) (divide by \(\lVert b \rVert\) for a geometric distance). Rows with \(y_i m_i \ge 1\) are safely classified and cost nothing in the hinge loss.
Use it when you predict with svmFit coefficients (take the sign), rank rows by confidence, or find the support vectors: the rows with \(y_i m_i < 1\), which are the only ones shaping the fit. It is not a probability; for that fit logisticFit and use logisticProbability.
classificationmodel-interpretation
def svmMargin(b, x) ≔ apply(x, b)
svmFit
svmFit(t, r, l, x, y) · t number of turns scalar · r learning rate scalar · l L2 penalty scalar · x design matrix n×p · y labels −1 or +1 n → coefficients p
Let \(t \in \mathbb{N}\), let \(r > 0\) be the learning rate, \(l \ge 0\) the L2 penalty, \(X \in \mathbb{R}^{n \times p}\) (the argument x) the design matrix with rows \(x_i^{\top}\), and \(y \in \{-1, +1\}^n\) the labels. Let \(b^{(0)} = 0 \in \mathbb{R}^p\) and
Then \(\operatorname{svmFit}(t, r, l, X, y) = b^{(t)}\). A row exactly on the margin, \(y_i x_i^{\top} b = 1\), contributes nothing to the subgradient.
The SVM coefficients after t subgradient steps of size r from zero.
Intuition The hinge loss ignores rows already beyond the margin and charges linearly for the rest, so only the hard cases near the boundary shape the answer. The L2 penalty prefers a small \(\lVert b \rVert\), which is the same as a wide margin. Each turn shrinks b by the factor \(1 - rl\) and pushes it toward the rows that are misclassified or too close.
Use it when you want a robust linear classifier and do not need probabilities, especially with many features such as text. Labels must be −1/+1, and features should be standardised. A fixed-step subgradient method hovers around the optimum rather than settling on it, so use a small rate (\(rl\) well below 1) and many turns. For probabilities use logisticFit.
classificationtext-classificationoverfitting-control
def svmFit(t, r, l, x, y) ≔ svmLoop(t, r, l, x, y, zeroed(moment(x, y)))
Building blocks
svmGradient(l, b, x, y)— Let \(l \ge 0\), \(b \in \mathbb{R}^p\), \(X \in \mathbb{R}^{n \times p}\) with rows \(x_i^{\top}\), and \(y \in \{-1, +1\}^n\). Then \(\operatorname{svmGradient}(l, b, X, y) = -\frac{1}{n}\sum_{i=1}^{n} y_i\,[\,y_i\, x_i^{\top} b < 1\,]\, x_i + l\,b\). The subgradient: \(X^{\top}\) applied todHingeof the margins, plus the penalty term.svmTurn(r, l, x, y, b)— Let \(r > 0\) and \(l\), \(X\), \(y\), \(b\) as above. Then \(\operatorname{svmTurn}(r, l, X, y, b) = b - r\,\operatorname{svmGradient}(l, b, X, y)\). Onedescendstep alongsvmGradient.svmLoop(t, r, l, x, y, b)— Let \(t \in \mathbb{N}\) and the rest as above. Then \(\operatorname{svmLoop}(0, r, l, X, y, b) = b\) and \(\operatorname{svmLoop}(t, r, l, X, y, b) = \operatorname{svmLoop}(t - 1, r, l, X, y, \operatorname{svmTurn}(r, l, X, y, b))\). RepeatssvmTurnttimes.
Kernel ridge, and the Gaussian process it turns into
Kernel ridge regression is ridge regression on a kernel matrix. With \(K\) the n×n kernel between training points, the dual weights and predictions are
\[\alpha = (K + \lambda I)^{-1} y, \qquad \hat f(z) = K_{z x}\, \alpha\]\(K + \lambda I\) is symmetric positive-definite for any \(\lambda > 0\), so solve (conjugate gradient) applies. A Gaussian process with an RBF kernel \(k(u, v) = \exp(-g\lVert u - v \rVert^2)\), zero prior mean and noise variance \(\lambda\) has exactly this prediction as its posterior mean, and it also reports a posterior variance:
The leading 1 is the prior variance of an RBF kernel. The variance is that of the underlying function (the noise \(\lambda\) is not added back) and is floored at zero, since rounding can make it slightly negative at a training point.
kernelWeights
kernelWeights(k, l, y) · k kernel matrix n×n · l ridge penalty, positive scalar · y responses n → dual weights n
Let \(K \in \mathbb{R}^{n \times n}\) (the argument k) be a symmetric positive semi-definite kernel matrix, let \(\lambda > 0\) (the argument l) be the ridge penalty, and let \(y \in \mathbb{R}^n\). Then
Solves \((K + lI)\alpha = y\). k is any symmetric positive semi-definite kernel matrix, for example rbf(x, x, g) or linear(x, x) from Prelude.Kernels.
Intuition Instead of one coefficient per feature there is one weight per training point, and a prediction is a similarity-weighted sum of those weights. The ridge l trades exact interpolation (small l, which chases noise) for smoothness (large l), and keeps the solve well posed. In the Gaussian-process reading, l is the noise variance.
Use it when the relationship is non-linear and n is moderate: it solves an n×n system, so thousands of points, not millions. Choose the kernel, its width and l by validation, and scale features first because the kernel works on distances. It gives no error bars; with an RBF kernel, gpMean and gpVariance are the same fit with uncertainty.
nonlinear-regressionregularisationregression
def kernelWeights(k, l, y) ≔ solve(k + l · identity(k), y)
kernelPredict
kernelPredict(kx, a) · kx kernel between queries and training points m×n · a dual weights n → predictions m
Let \(K_{zx} \in \mathbb{R}^{m \times n}\) (the argument kx) be the kernel between \(m\) queries and \(n\) training points, and let \(\alpha \in \mathbb{R}^n\) (the argument a) be the dual weights. Then
The kernel-ridge prediction \(K_{zx}\alpha\). Build kx with the same kernel as the training matrix, queries first: for example rbf(z, x, g).
Intuition Every training point contributes its dual weight, scaled by how similar the query is to it. With an RBF kernel distant points contribute nothing, so far from all the data the prediction falls back to 0.
Use it when you score new rows after kernelWeights. Build kx with exactly the kernel and width used for the training matrix, queries as rows. Every prediction touches every training point (m×n kernel values), so there is no compact model to ship. Because predictions revert to 0 away from the data, subtract the mean of y before fitting and add it back to the predictions.
nonlinear-regressioninterpolation
def kernelPredict(kx, a) ≔ einsum("ij,j->i", kx, a)
gpMean
gpMean(g, l, x, y, z) · g RBF inverse width scalar · l noise variance scalar · x training inputs n×p · y training targets n · z query inputs m×p → posterior means m
Let \(g > 0\) be the inverse width, \(\lambda > 0\) (the argument l) the noise variance, \(X \in \mathbb{R}^{n \times p}\) (the argument x) the training inputs with rows \(x_j^{\top}\), \(y \in \mathbb{R}^n\) the targets, and \(Z \in \mathbb{R}^{m \times p}\) (the argument z) the queries with rows \(z_i^{\top}\). Let \(K_{jk} = \exp(-g\lVert x_j - x_k \rVert^2)\) and \((K_{zx})_{ij} = \exp(-g\lVert z_i - x_j \rVert^2)\). Then
The posterior mean of a zero-mean Gaussian process with kernel \(\exp(-g\lVert u - v\rVert^2)\) at each query row.
Intuition A Gaussian process is a prior over smooth functions: g sets how fast the function may change (larger g, wigglier) and l how much each observation is trusted. The posterior mean is the best guess given the data and reverts to the prior mean 0 where there is none. It is kernel ridge regression with an RBF kernel.
Use it when you model a smooth unknown function from few, expensive observations: calibration curves, spatial fields from sparse sensors, surrogate models of slow simulations. Subtract the mean of y first and add it back, since the prior mean is 0. It solves an n×n system. Pair it with gpVariance for error bars.
nonlinear-regressionuncertainty-quantificationsurrogate-modelling
def gpMean(g, l, x, y, z) ≔ kernelPredict(gpCross(g, x, z), kernelWeights(gpPrior(g, x), l, y))
gpVariance
gpVariance(g, l, x, z) · g RBF inverse width scalar · l noise variance scalar · x training inputs n×p · z query inputs m×p → posterior variances m
Let \(g\), \(\lambda\), \(X\), \(Z\), \(K\) and \(K_{zx}\) be as in gpMean, and let \(k_i \in \mathbb{R}^n\) be the \(i\)-th row of \(K_{zx}\). Then, for \(i = 1, \dots, m\),
where \(\operatorname{inverse}\) is the 40-turn Newton–Schulz inverse of Prelude.LinearAlgebra.
The posterior variance at each query row: near 0 where the data pins the function down, rising to 1 far from every training point. It does not depend on the targets. The inverse of \(K + lI\) is formed explicitly (Newton–Schulz), so keep l large enough for the matrix to be reasonably conditioned.
Intuition How much the data has not pinned the function down: 1 (the prior variance) far from every training input, near 0 at a well-observed one. It depends only on where you have measured, not on the values seen, so it can be computed before taking a new measurement.
Use it when you need error bars (mean \(\pm 2\sqrt{\text{variance}}\) is a rough 95 % band for the function) or must decide where to measure next: sample where the variance is highest, as in active learning and Bayesian optimisation. It is the variance of the function; add l for a new noisy observation. The prior variance is fixed at 1, so standardise y to match. Avoid a tiny l: the explicit inverse loses accuracy on an ill-conditioned kernel matrix.
uncertainty-quantificationactive-learningexperimental-design
def gpVariance(g, l, x, z) ≔ max(1 − gpResolved(g, l, x, z), 0)
Building blocks
gpPrior(g, x)— Let \(g > 0\) and \(X \in \mathbb{R}^{n \times p}\) with rows \(x_j^{\top}\). Then \(\operatorname{gpPrior}(g, X) = K\) with \(K_{jk} = \exp(-g\lVert x_j - x_k \rVert^2)\). The n×n RBF kernel matrix of the training inputs.gpCross(g, x, z)— Let \(g > 0\), \(X \in \mathbb{R}^{n \times p}\) with rows \(x_j^{\top}\) and \(Z \in \mathbb{R}^{m \times p}\) with rows \(z_i^{\top}\). Then \(\operatorname{gpCross}(g, X, Z) = K_{zx}\) with \((K_{zx})_{ij} = \exp(-g\lVert z_i - x_j \rVert^2)\). The m×n RBF kernel between query and training inputs.gpGuard(g, l, x)— Let \(g > 0\), \(\lambda > 0\) (the argumentl) and \(X\) as above. Then \(\operatorname{gpGuard}(g, \lambda, X) = K + \lambda I_n\). The regularised kernel matrix.gpResolved(g, l, x, z)— Let \(g\), \(\lambda\), \(X\), \(Z\) be as above, \(M = \operatorname{inverse}(K + \lambda I_n)\) and \(k_i^{\top}\) the rows of \(K_{zx}\). Then \(\operatorname{gpResolved}(g, \lambda, X, Z)_i = k_i^{\top} M k_i = \sum_{j=1}^{n}\sum_{j'=1}^{n} (K_{zx})_{ij} M_{jj'} (K_{zx})_{ij'}\). The explained variance for each query row.
Examples
The C# snippets below pass matrices and vectors as Tensor (for example (Tensor)new double[,] { … }), because every fit returns a result of a different rank from its inputs.
Fit a logistic regression
Six observations with an intercept column and one feature; the labels overlap, so the maximum likelihood exists. The fit is already converged after 10 turns (25 turns give the same digits). A logistic fit with an intercept reproduces the number of positive labels: the fitted probabilities sum to 3.
open Prelude.Regression
fn fitLogistic(t, x, y) ↦ logisticFit(t, x, y)
fn fittedChances(t, x, y) ↦ logisticProbability(logisticFit(t, x, y), x)
fn fitCounts(t, x, y) ↦ poissonFit(t, x, y)
using MLambda.Hilbert.Runtime;
using Examples;
var x = (Tensor)new double[,] { { 1, 0.5 }, { 1, 1.5 }, { 1, 2 }, { 1, 3 }, { 1, 3.5 }, { 1, 4 } };
var y = (Tensor)new[] { 0d, 0d, 1d, 0d, 1d, 1d };
var b = Glm.FitLogistic(10.0, x, y);
var chances = Glm.FittedChances(10.0, x, y);
Console.WriteLine(chances.Values.Sum());
Output (verified, six decimals):
b = [-3.455893, 1.405061]
chances = [0.059897, 0.206148, 0.343943, 0.681199, 0.811811, 0.897002]
sum = 3
Fit a Poisson regression to counts
Two groups with an intercept and a group indicator: counts 1 and 3 in the first group (mean 2) and 4 and 8 in the second (mean 6). The fitted coefficients are \(\ln 2 = 0.693147\) and \(\ln(6/2) = 1.098612\). This uses fitCounts from the Glm.hb file above.
var x = (Tensor)new double[,] { { 1, 0 }, { 1, 0 }, { 1, 1 }, { 1, 1 } };
var y = (Tensor)new[] { 1d, 3d, 4d, 8d };
var b = Glm.FitCounts(10.0, x, y);
Output (verified):
[0.693147, 1.098612]
Select variables with the lasso
With an orthogonal design, \(X^\top X/n = \tfrac12 I\) and \(X^\top y/n = (1.5, 0.5)\), so the lasso solution is the soft threshold of \((1.5, 0.5)\) at \(l = 0.6\), divided by \(\tfrac12\): \((1.8, 0)\). The second variable is dropped exactly. Adding an L2 penalty \(a = 0.5\) doubles the curvature and halves the surviving coefficient to 0.9.
open Prelude.Regression
fn fitLasso(t, r, l, x, y) ↦ lassoFit(t, r, l, x, y)
fn fitElastic(t, r, l, a, x, y) ↦ elasticFit(t, r, l, a, x, y)
using MLambda.Hilbert.Runtime;
using Examples;
var x = (Tensor)new double[,] { { 1, 0 }, { 0, 1 }, { -1, 0 }, { 0, -1 } };
var y = (Tensor)new[] { 3d, 1d, -3d, -1d };
var lasso = Sparse.FitLasso(60.0, 1.0, 0.6, x, y);
var elastic = Sparse.FitElastic(60.0, 1.0, 0.6, 0.5, x, y);
Output (verified):
lasso = [1.8, 0]
elastic = [0.9, 0]
Separate two classes with a linear SVM
Four points on a line, labelled −1 on the left and +1 on the right, with an intercept column. By symmetry the intercept is zero; the signs of the margins recover every label.
open Prelude.Regression
fn fitMargin(t, r, l, x, y) ↦ svmFit(t, r, l, x, y)
fn marginSides(t, r, l, x, y) ↦ sign(apply(x, svmFit(t, r, l, x, y)))
using MLambda.Hilbert.Runtime;
using Examples;
var x = (Tensor)new double[,] { { 1, -2 }, { 1, -1 }, { 1, 1 }, { 1, 2 } };
var y = (Tensor)new[] { -1d, -1d, 1d, 1d };
var b = Margin.FitMargin(200.0, 0.5, 0.01, x, y);
var sides = Margin.MarginSides(200.0, 0.5, 0.01, x, y);
Output (verified, six decimals):
b = [0, 1.136044]
sides = [-1, -1, 1, 1]
Predict with uncertainty using a Gaussian process
Three training points \((0, 0), (1, 1), (2, 0)\), kernel width \(g = 1\) and noise \(l = 0.01\). At the training input 1 the mean is close to the target 1 and the variance is small; halfway to the next point the variance grows; at 10, far from all data, the mean returns to the prior mean 0 and the variance to the prior variance 1.
open Prelude.Regression
fn posteriorMean(g, l, x, y, z) ↦ gpMean(g, l, x, y, z)
fn posteriorVariance(g, l, x, z) ↦ gpVariance(g, l, x, z)
using MLambda.Hilbert.Runtime;
using Examples;
var x = (Tensor)new double[,] { { 0 }, { 1 }, { 2 } };
var y = (Tensor)new[] { 0d, 1d, 0d };
var z = (Tensor)new double[,] { { 1 }, { 1.5 }, { 10 } };
var mean = Process.PosteriorMean(1.0, 0.01, x, y, z);
var variance = Process.PosteriorVariance(1.0, 0.01, x, z);
Output (verified, rounded to six decimals):
mean = [0.986609, 0.619295, 0]
variance = [0.009866, 0.104812, 1]