Multivariate
Principal components, whitening, Mahalanobis distance, ordinary, weighted and ridge least squares, and linear discriminants over a data matrix.
This module collects the classical multivariate methods: principal component analysis and whitening, Mahalanobis distances, ordinary least squares with its fit statistics, weighted and ridge regression, and Fisher’s linear discriminants. Each is a short composition of the sample covariance matrix from Statistics and the decompositions of Linear algebra, so opening this module gives you those two as well.
The data matrix x is always n×p: one row per observation, one column per variable. A response
y is a vector n, weights w are a vector n, and group membership g is a one-hot
matrix n×k (a 1 in the column of each observation’s group). Regression functions do not add an
intercept: if you want one, include a column of ones in x. rsquared measures the fit against
the mean of y, which is the usual definition when that intercept column is present.
The normal equations are solved with solve, which is conjugate gradient and needs a symmetric
positive-definite matrix — in practice, the columns of x must be linearly independent (or, for
ridge, the penalty must be positive). Inverses (for standard errors, Mahalanobis distances and
discriminants) use a fixed Newton–Schulz iteration and eigenvectors a fixed number of turns of the
QR algorithm, so both are accurate for reasonably conditioned matrices and less so when the
covariance is nearly singular or has two nearly equal eigenvalues. The sign of an eigenvector is
arbitrary.
Throughout this page, \(X \in \mathbb{R}^{n \times p}\) is the data matrix with rows \(x_i^{\top}\),
\(\bar{x} = \frac1n \sum_{i=1}^{n} x_i \in \mathbb{R}^p\) its column means (means),
\(S = \operatorname{covariances}(X) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(x_i - \bar{x})^{\top}\)
its sample covariance matrix, and \(\operatorname{eigenvectors}\), \(\operatorname{eigenvalues}\),
\(\operatorname{inverse}\) and \(\operatorname{solve}\) are the fixed-turn QR algorithm, Newton–Schulz
inverse and Jacobi-preconditioned conjugate gradient of Linear algebra. In exact
arithmetic and at convergence, \(\operatorname{inverse}(A) = A^{-1}\) and
\(\operatorname{solve}(A, b) = A^{-1} b\) for symmetric positive-definite \(A\).
open Prelude.Multivariate
Principal components
Principal component analysis rotates the data onto the eigenvectors of its covariance matrix \(S = V \Lambda V^\top\), ordered by descending eigenvalue, so the first component carries the most variance. The eigenvalues are the variances along each component.
components
components(x) · x data matrix n×p → principal axes p×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the data matrix, \(n \ge 2\), with sample covariance \(S\). Then
where \(v_j\) is the \(j\)-th column of \(V\).
The eigenvectors of the sample covariance matrix, one per column, in descending order of explained variance, so the first column is the first principal component.
Intuition Picture the data as a cloud of points, usually shaped like a tilted ellipsoid. The columns of \(V\) are the axes of that ellipsoid, longest first: a rotation of the coordinate system so the new variables are uncorrelated and ordered by how much the data varies along them.
Use it when you want to know which combinations of variables carry the structure, before plotting, compressing or regressing on fewer inputs. The covariance is in the original units, so a variable measured in large numbers dominates: standardise columns first when their units differ. Axes come from 40 fixed QR turns and their signs are arbitrary. For the variances along the axes use explained, for the data in these axes use scores.
dimensionality-reductionfeature-decorrelationexploratory-analysis
def components(x) ≔ eigenvectors(covariances(x))
explained
explained(x) · x data matrix n×p → variance per component vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(n \ge 2\), with sample covariance \(S\). Then
where \(V = \operatorname{eigenvectors}(S)\).
The eigenvalues of the sample covariance matrix (divisor \(n-1\)), in the same order as
components: the variance of the data along each principal axis.
Intuition Each eigenvalue is the squared length of one axis of the data ellipsoid: how much the points spread out along that principal direction. The values add up to the total variance of all columns, so rotating onto the components only reshuffles variance, it creates or loses none.
Use it when you need the absolute variance per component, for a scree plot, for spotting a near-zero eigenvalue that reveals collinear columns, or to scale scores yourself (as whitened does). To decide how many components to keep, the proportions from explainedPart are easier to read.
dimensionality-reductionmulticollinearityexploratory-analysis
def explained(x) ≔ eigenvalues(covariances(x))
explainedPart
explainedPart(x) · x data matrix n×p → proportion of variance per component vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) and \(\lambda = \operatorname{explained}(X)\). Then, for \(j = 1, \dots, p\),
The explained variances divided by their total, so they sum to 1. Use it to decide how many components to keep.
Intuition The share of the total variance that each component accounts for. A steep drop after the first few values means the data really live in a few dimensions and the rest is mostly noise.
Use it when you choose how many components to keep: a common rule is to keep the leading ones until their running total reaches 80–95%, or to stop where the values level off. The shares depend on scaling, so standardise columns first if their units differ. For the variances themselves use explained.
dimensionality-reductionmodel-selection
def explainedPart(x) ≔ explained(x) / total(explained(x))
scores
scores(x) · x data matrix n×p → component scores n×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with column means \(\bar{x}\), and \(V = \operatorname{components}(X)\). Then
The centred data expressed in the principal axes. Column j is each observation’s coordinate along component j; keep the first few columns for a reduced representation.
Intuition The same points, described in the ellipsoid’s own coordinates instead of the original variables. The columns are uncorrelated, and the first ones hold most of the variation, so dropping the trailing columns loses as little as any linear projection can.
Use it when you want a low-dimensional version of the data: a 2-D scatter plot of the first two columns, compact inputs for a regression or clustering (principal component regression), or removing correlation between features. It returns all p columns; keep the first k yourself. For unit variance in every column use whitened.
dimensionality-reductionvisualisationfeature-decorrelation
def scores(x) ≔ mul(centred(x), components(x))
restored
restored(x) · x data matrix n×p → reconstructed data n×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with column means \(\bar{x}\), \(V = \operatorname{components}(X)\) and \(Z = \operatorname{scores}(X)\). Then
The scores rotated back and the column means added again. Because V is orthogonal, \(V V^{\top} = I_p\) and this reproduces x; it is mainly a check that the rotation lost nothing.
Intuition Rotating back undoes the rotation: a round trip data → scores → data. Because every component is kept, any difference from x is numerical error from the eigenvector iteration, not information lost.
Use it when you want to confirm that components converged on your data: a large gap between restored(x) and x signals nearly equal eigenvalues or poor conditioning. It does not truncate, so it is not a low-rank reconstruction; for compression or denoising, zero the trailing columns of scores and rotate back yourself.
numerical-stabilitysanity-checking
def restored(x) ≔ mul(scores(x), transpose(components(x))) + means(x)
whitened
whitened(x) · x data matrix n×p → whitened data n×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(Z = \operatorname{scores}(X)\) and \(\lambda = \operatorname{explained}(X)\). Then, for \(i = 1, \dots, n\) and \(j = 1, \dots, p\),
The scores with each column divided by its own standard deviation (the square root of its eigenvalue), so the result has the identity as its sample covariance matrix whatever the data’s was. A zero eigenvalue is clamped to a tiny positive number so the division cannot fail, but a direction with no variance cannot be meaningfully whitened.
Intuition Rotate the ellipsoid onto its axes, then squash or stretch each axis to length one, and the cloud becomes a round ball. After that, ordinary Euclidean distance between rows is Mahalanobis distance in the original data, and no variable matters more just because of its scale or its correlations.
Use it when a method assumes uncorrelated, unit-variance inputs: independent component analysis, k-means or nearest neighbours on correlated features, or conditioning inputs for gradient descent. Directions with near-zero eigenvalues get blown up, which amplifies noise, so drop those columns first (see explainedPart). If you only need decorrelation without rescaling, use scores.
feature-decorrelationfeature-scalingnumerical-stability
def whitened(x) ≔ scores(x) / sqrt(max(explained(x), 1e-300))
Mahalanobis distance
The Mahalanobis distance measures how far each observation lies from the column means in units of the data’s own covariance, so correlated or differently scaled variables do not dominate: \(d_i^2 = (x_i - \bar{x})^\top S^{-1} (x_i - \bar{x})\). It is computed for every row at once.
precision
precision(x) · x data matrix n×p → precision matrix p×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with invertible sample covariance \(S\). Then
The inverse of the sample covariance matrix. The covariance must be invertible, which needs more observations than variables and no exactly collinear columns.
Intuition Covariance says how variables move together. Its inverse says how each pair relates once all the other variables are held fixed: a zero off-diagonal entry means the two are conditionally uncorrelated (independent, for normal data), even if they are correlated overall.
Use it when you need Mahalanobis distances (mahalanobis uses it), a Gaussian log-density, or a map of direct dependencies between variables (Gaussian graphical models). It is a fixed 40-step Newton–Schulz inverse, so accuracy falls as the covariance approaches singularity. With about as many variables as observations, shrink the covariance first, as ridge does for regression, or reduce dimensions with scores.
conditional-dependenceanomaly-detectionnumerical-stability
def precision(x) ≔ inverse(covariances(x))
mahalanobis
mahalanobis(x) · x data matrix n×p → squared distance per row vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with rows \(x_i^{\top}\) and column means \(\bar{x}\), and let \(P = \operatorname{precision}(X)\). Then, for \(i = 1, \dots, n\),
The squared Mahalanobis distance of each row from the column means, using the sample covariance of the same data.
Intuition A squared z-score for many variables at once. It measures distance from the centre in “standard deviations along the data’s own ellipsoid”, so a point that is ordinary on each variable separately but breaks their usual correlation still comes out far away.
Use it when you screen for multivariate outliers or score how atypical each observation is. For roughly normal data the squared distances are approximately chi-square on p degrees of freedom, so they average about p. The mean and covariance come from the same rows, including any outliers, which pull the centre and inflate the covariance and can hide the very points you are looking for. For the unsquared distance use apart.
anomaly-detectiondata-cleaningsimilarity-measurement
def mahalanobis(x) ≔ einsum("ij,jk,ik->i", centred(x), precision(x), centred(x))
apart
apart(x) · x data matrix n×p → distance per row vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) and \(m = \operatorname{mahalanobis}(X) \in \mathbb{R}^n\). Then, for \(i = 1, \dots, n\),
The Mahalanobis distance itself, the square root of mahalanobis (tiny negative rounding errors are
clamped to zero first). Large values flag multivariate outliers.
Intuition The same distance as mahalanobis, back in “standard deviation” units, the way a z-score is. A value of 3 reads as “three standard deviations out along the data’s own shape”, which is easier to put a threshold on than its square.
Use it when you report or threshold how unusual each row is, or rank rows by atypicality. For a statistical cut-off, compare the squared distances from mahalanobis with a chi-square quantile on p degrees of freedom. The same caution applies: outliers in x shape the covariance used to judge them.
anomaly-detectiondata-cleaning
def apart(x) ≔ sqrt(max(mahalanobis(x), 0))
Ordinary least squares
Least squares finds the coefficients \(\hat\beta\) minimising \(\lVert y - X\beta \rVert^2\) by solving the normal equations \(X^\top X \hat\beta = X^\top y\). The functions below give the coefficients, the fitted values and residuals, the goodness of fit, and the standard errors and t-statistics of the coefficients.
coefficients
coefficients(x, y) · x design matrix n×p · y response vector n → coefficients vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix, with \(X^{\top} X\) positive-definite, and let \(y \in \mathbb{R}^n\) be the response. Then
The least-squares coefficients, one per column of x. Include a column of ones in x for an intercept. \(X^\top X\) must be positive-definite, so the columns must be linearly independent.
Intuition The fitted values \(X\hat\beta\) are the point in the span of the columns closest to y: the residual is perpendicular to every column, which is exactly what the normal equations say. Each coefficient is the change in y per unit change of its column with the other columns held fixed.
Use it when you fit a linear model: trend lines, calibration curves, polynomial fits (put powers of t in the columns), or effects of several predictors at once. It is solved by conjugate gradient in p steps without forming an inverse. With collinear or nearly collinear columns the answer is unstable; use ridge. When observations have unequal reliability use weightedCoefficients.
regressionparameter-estimationcurve-fitting
def coefficients(x, y) ≔ solve(gram(x), moment(x, y))
fitted
fitted(x, y) · x design matrix n×p · y response vector n → fitted values vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(y \in \mathbb{R}^n\) and \(\hat\beta = \operatorname{coefficients}(X, y)\). Then
The predictions of the fitted model at the observed rows.
Intuition The shadow of y on the space the columns of x can reach: the part of the response the model can explain. The rest is residuals, and the two add back to y.
Use it when you plot the fit against the data, compute your own error measures, or feed smoothed values onward. It predicts only at the rows used for fitting, so it flatters the model; for rows the model has not seen, multiply their design rows by coefficients, and judge accuracy on held-out data (see crossError in Selection).
regressionprediction
def fitted(x, y) ≔ apply(x, coefficients(x, y))
residuals
residuals(x, y) · x design matrix n×p · y response vector n → residuals vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(y \in \mathbb{R}^n\) and \(\hat\beta = \operatorname{coefficients}(X, y)\). Then
Observed minus fitted.
Intuition What the model failed to explain. If the model is right, residuals look like structureless noise; any pattern left in them is structure the model missed.
Use it when you check the model’s assumptions: plot them against fitted values or each predictor to spot curvature (a missing term), a funnel shape (unequal variance, try weightedCoefficients) or a trend over time (correlated errors). Single large values point to outliers or data errors. They sum to zero only when x includes an intercept column. For a single summary number use rss.
model-validationregressionanomaly-detection
def residuals(x, y) ≔ y − fitted(x, y)
rss
rss(x, y) · x design matrix n×p · y response vector n → residual sum of squares scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(y \in \mathbb{R}^n\) and \(e = \operatorname{residuals}(X, y)\). Then
The sum of squared residuals.
Intuition The quantity least squares makes as small as possible: the squared length of the residual vector. Squaring punishes a few big misses more than many small ones, which is why a single outlier can move a least-squares fit a long way.
Use it when you compare fits on the same data, for example nested models (the drop in RSS is what an F-test weighs), or build other statistics: rsquared, spread and information criteria all start from it. It grows with n and is in squared units of y, so it is not comparable across data sets; use rsquared or spread for that.
goodness-of-fitmodel-selection
def rss(x, y) ≔ dot(residuals(x, y), residuals(x, y))
tss
tss(y) · y response vector n → total sum of squares scalar
Let \(y \in \mathbb{R}^n\) and \(\bar{y} = \frac1n \sum_{i=1}^{n} y_i\). Then
\[\operatorname{tss}(y) = \sum_{i=1}^{n} (y_i - \bar{y})^2 .\]The sum of squared deviations of y from its mean.
Intuition The RSS of the laziest possible model, one that ignores x and predicts the mean every time. It is the baseline any real model must beat, and it is \((n-1)\) times the sample variance of y.
Use it when you need that baseline, mostly to form \(R^2\) as the fraction of it a model removes (rsquared already does this). It depends only on y, so it is the same for every model you fit to one response and can be computed once.
goodness-of-fitbaseline-comparison
def tss(y) ≔ dot(y − mean(y), y − mean(y))
rsquared
rsquared(x, y) · x design matrix n×p · y response vector n → coefficient of determination scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) and \(y \in \mathbb{R}^n\) not constant. Then
The proportion of the variance of y the fit explains. It is measured against the mean of y, so it is meaningful when x contains an intercept column.
Intuition How much of the variation in y the model removes compared with just predicting the mean: 0 means no better than the mean, 1 means a perfect fit. It is a ratio, so it is unit-free and comparable across responses.
Use it when you summarise how well a linear model describes its data. It never decreases when you add a column, even a column of pure noise, so do not use it to choose between models with different numbers of columns; use adjusted or cross-validation (crossError in Selection). Without an intercept column it can be negative or misleading. A high value says nothing about whether the coefficients are significant (see statistics).
goodness-of-fitregression
def rsquared(x, y) ≔ 1 − rss(x, y) / tss(y)
adjusted
adjusted(x, y) · x design matrix n×p · y response vector n → adjusted R² scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with \(n > p + 1\), \(y \in \mathbb{R}^n\) and \(R^2 = \operatorname{rsquared}(X, y)\). Then
\(R^2\) penalised for the number of columns, with p read as the number of columns of x. The textbook formula counts predictors excluding the intercept; since here p is every column of x, a design matrix that includes an intercept column is penalised by one more degree of freedom than that formula.
Intuition \(R^2\) with a charge for every column: the unexplained fraction is scaled up by \((n-1)/(n-p-1)\), so a new column must remove more variance than chance alone would before the score goes up.
Use it when you compare linear models with different numbers of columns on the same data, and want a quick guard against adding useless predictors. It needs \(n > p + 1\) and becomes unstable as p approaches n. It is a rough correction; for choosing a model that predicts well, prefer held-out error (crossError in Selection).
model-selectiongoodness-of-fitoverfitting-control
def adjusted(x, y) ≔ 1 − (1 − rsquared(x, y)) · (rows(x) − 1) / (rows(x) − columns(x) − 1)
spread
spread(x, y) · x design matrix n×p · y response vector n → residual variance scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with \(n > p\), and \(y \in \mathbb{R}^n\). Then
The unbiased estimate of the error variance: p degrees of freedom went into the fit, so the divisor is \(n-p\), not n.
Intuition The typical size of the noise around the fitted line, squared. Fitting p coefficients lets the model bend toward the data and makes the residuals look smaller than the true errors, and dividing by \(n-p\) instead of n corrects for exactly that.
Use it when you need the noise level of a regression: its square root is the residual standard error, in the units of y, which tells you how far predictions typically miss. It is the scale factor inside errors, and it is the variance to use when simulating from the fitted model. It needs \(n > p\).
uncertainty-quantificationregressionparameter-estimation
def spread(x, y) ≔ rss(x, y) / (rows(x) − columns(x))
errors
errors(x, y) · x design matrix n×p · y response vector n → standard errors vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(y \in \mathbb{R}^n\), \(\hat\sigma^2 = \operatorname{spread}(X, y)\) and \(C = \operatorname{inverse}(X^{\top} X) \approx (X^{\top} X)^{-1}\). Then, for \(j = 1, \dots, p\),
The standard error of each coefficient, the square roots of the diagonal of \(\hat\sigma^2 (X^\top X)^{-1}\). This is the one place the module forms an inverse of \(X^\top X\), so it inherits the inverse’s accuracy on ill-conditioned designs.
Intuition How much each coefficient would wobble if you collected the data again. It grows with the noise level (spread) and with how little independent information a column carries: a column that is nearly a combination of the others has a large diagonal entry in \((X^\top X)^{-1}\), so its coefficient is poorly pinned down.
Use it when you want confidence intervals for coefficients (coefficient plus or minus a t critical value times its error) or to diagnose collinearity. The formula assumes independent errors with equal variance; with a funnel in the residuals or correlated observations these errors are too small. For tests, statistics divides through for you.
uncertainty-quantificationconfidence-intervalsmulticollinearity
def errors(x, y) ≔ sqrt(max(spread(x, y) · diag(inverse(gram(x))), 0))
statistics
statistics(x, y) · x design matrix n×p · y response vector n → t-statistics vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x), \(y \in \mathbb{R}^n\), \(\hat\beta = \operatorname{coefficients}(X, y)\) and \(s = \operatorname{errors}(X, y)\). Then, for \(j = 1, \dots, p\),
Each coefficient divided by its standard error, for testing whether it differs from zero. Compare
against Student’s t with \(n-p\) degrees of freedom (studentCdf in
Distributions).
Intuition Signal over noise for each coefficient: how many standard errors the estimate sits away from zero. A value near zero means the data cannot tell the effect apart from no effect; beyond about 2 in size it usually can.
Use it when you ask whether each predictor contributes once the others are in the model. For p-values use twoSided from Inference with \(n-p\) degrees of freedom. Testing many coefficients at once inflates false positives, so correct with bonferroni or discoveries. A small t for two correlated predictors can mean they share one effect, not that neither matters.
hypothesis-testingfeature-selectionregression
def statistics(x, y) ≔ coefficients(x, y) / errors(x, y)
Building blocks
moment(x, y)— Let \(X \in \mathbb{R}^{n \times p}\) (the argumentx) and \(y \in \mathbb{R}^n\). Then \(\operatorname{moment}(X, y) = X^{\top} y\), with entries \((X^{\top} y)_j = \sum_{i=1}^{n} X_{ij}\, y_i\). The right-hand side of the normal equations, a vector of length p; used bycoefficientsandridge.
Weighted least squares
Weighted least squares gives each observation a non-negative weight \(w_i\) and minimises \(\sum_i w_i (y_i - x_i^\top \beta)^2\). The normal equations become \(X^\top W X \hat\beta = X^\top W y\) with W diagonal. A weight of zero removes an observation.
weightedCoefficients
weightedCoefficients(w, x, y) · w weights vector n · x design matrix n×p · y response vector n → coefficients vector p
Let \(w \in \mathbb{R}^n\) with \(w_i \ge 0\) be the weights and \(W = \operatorname{diag}(w)\), let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with \(X^{\top} W X\) positive-definite, and let \(y \in \mathbb{R}^n\). Then
The weighted least-squares coefficients. With all weights equal it gives the same answer as
coefficients.
Intuition Least squares where each observation gets a vote of size \(w_i\): a row with weight 2 counts like two copies of itself, and weight 0 removes it. Only the relative sizes of the weights matter, so multiplying all of them by a constant changes nothing.
Use it when observations have different reliability (set \(w_i = 1/\sigma_i^2\) when you know each measurement’s variance), rows stand for grouped counts, you fit locally with kernel weights, or you want to drop rows without copying data (Selection’s cross-validation holds rows out this way). It is also the inner step of iteratively reweighted methods such as robust regression. With equal weights use coefficients.
regressionunequal-variancerobust-statistics
def weightedCoefficients(w, x, y) ≔ solve(weightedGram(w, x), weightedMoment(w, x, y))
Building blocks
weightedGram(w, x)— Let \(w \in \mathbb{R}^n\), \(W = \operatorname{diag}(w)\) and \(X \in \mathbb{R}^{n \times p}\) (the argumentx) with rows \(x_i^{\top}\). Then \(\operatorname{weightedGram}(w, X) = X^{\top} W X = \sum_{i=1}^{n} w_i\, x_i x_i^{\top}\), with entries \(\sum_{i} w_i X_{ij} X_{il}\). The weighted cross-product matrix, p×p; used byweightedCoefficients.weightedMoment(w, x, y)— Let \(w \in \mathbb{R}^n\), \(W = \operatorname{diag}(w)\), \(X \in \mathbb{R}^{n \times p}\) (the argumentx) and \(y \in \mathbb{R}^n\). Then \(\operatorname{weightedMoment}(w, X, y) = X^{\top} W y\), with entries \(\sum_{i=1}^{n} w_i X_{ij}\, y_i\). The weighted right-hand side, length p; used byweightedCoefficients.
Ridge regression
Ridge regression adds \(\lambda I\) to the normal equations. The answer is biased towards zero but better conditioned, which is the trade the method exists to make; the system also stays positive-definite, so it can be solved even when the columns of x are collinear.
ridge
ridge(l, x, y) · l penalty \(\lambda \ge 0\) scalar · x design matrix n×p · y response vector n → coefficients vector p
Let \(\lambda \ge 0\) (the argument l) be the penalty, \(X \in \mathbb{R}^{n \times p}\) (the argument x) with \(X^{\top} X + \lambda I_p\) positive-definite, and \(y \in \mathbb{R}^n\). Then
The penalty applies to every column
of x, including an intercept column if you supplied one; centre the data first if the intercept
should not be shrunk. With \(\lambda = 0\) it equals coefficients.
Intuition When columns are nearly collinear, least squares can explain the data equally well with huge positive and negative coefficients that cancel, and small changes in the data swing them wildly. Adding \(\lambda I\) charges for coefficient size, so among nearly equal fits it picks the one with small coefficients: a little bias bought for much lower variance.
Use it when predictors are correlated, there are many of them relative to n, or plain coefficients overfits. Standardise columns first, since the same penalty applies to every coefficient whatever its units. Choose \(\lambda\) by held-out error, not by in-sample fit: bestPenalty and tunedRidge in Selection do that over a grid.
multicollinearityoverfitting-controlregression
def ridge(l, x, y) ≔ solve(gram(x) + l · identity(gram(x)), moment(x, y))
Linear discriminants
Fisher’s linear discriminant analysis finds directions that separate k known groups. It splits the total scatter of the data into a between-group part B, from the group means, and a within-group part W, and looks for the eigenvectors of \(W^{-1}B\). Group membership is a one-hot matrix g n×k, so every per-group quantity is computed without a loop over groups. Every group must be non-empty.
In the statements below, \(G \in \{0, 1\}^{n \times k}\) (the argument g) is the one-hot membership
matrix, \(G_{ic} = 1\) exactly when observation \(i\) belongs to group \(c\), and groups are indexed
\(c = 1, \dots, k\).
sizes
sizes(g) · g one-hot group membership n×k → observations per group vector k
Let \(G \in \{0, 1\}^{n \times k}\) (the argument g) be one-hot. Then, for \(c = 1, \dots, k\),
The number of observations in each group.
Intuition Summing a one-hot matrix down its rows counts the ones in each column, which is the number of members of each group. No loop, no labels, just a column total.
Use it when you check a grouping before using it: an empty group makes groupMeans divide by zero, and very unequal sizes (class imbalance) mean the big groups dominate between. It is also the weight vector for turning group means back into totals.
class-imbalancesanity-checking
def sizes(g) ≔ einsum("ig->g", g)
groupMeans
groupMeans(g, x) · g one-hot group membership n×k · x data matrix n×p → group means k×p
Let \(G \in \{0, 1\}^{n \times k}\) (the argument g) be one-hot with every \(n_c > 0\), and \(X \in \mathbb{R}^{n \times p}\) (the argument x). Then, for \(c = 1, \dots, k\) and \(j = 1, \dots, p\),
The mean of each variable within each group; row c is the centroid \(\bar{x}_c\) of group c.
Intuition The centre of gravity of each group’s points. It summarises a whole group by one representative point, and the distances between these centroids are what “the groups differ” means geometrically.
Use it when you profile groups side by side, build a nearest-centroid classifier, or run the update step of k-means (with g as the current cluster assignment). Every group must have at least one member, or the division by its size fails; check with sizes. Means are pulled by outliers, so inspect groups with extreme points.
comparing-groupsclusteringclassification
def groupMeans(g, x) ≔ einsum("ig,ij,g->gj", g, x, 1 / sizes(g))
between
between(g, x) · g one-hot group membership n×k · x data matrix n×p → between-group scatter p×p
Let \(G\) (the argument g) and \(X\) (the argument x) be as in groupMeans, with group sizes \(n_c\), group means \(\bar{x}_c\) and column means \(\bar{x}\). Then
The scatter of the group means about the overall mean, each weighted by its group size.
Intuition Replace every point by its group’s centroid and measure how spread out that simplified data is: the variation that group membership alone explains. It is the multivariate version of ANOVA’s between-group sum of squares, and its rank is at most \(k-1\) because k centroids span at most \(k-1\) directions.
Use it when you measure how far apart groups are across several variables, as the numerator of linear discriminant analysis (discriminants) or of MANOVA-style comparisons. On its own it ignores how spread out each group is, so compare it with within before calling a separation large.
comparing-groupsclassification
def between(g, x) ≔ einsum("g,gj,gk->jk", sizes(g), groupApart(g, x), groupApart(g, x))
within
within(g, x) · g one-hot group membership n×k · x data matrix n×p → within-group scatter p×p
Let \(G\) (the argument g) and \(X\) (the argument x) be as in groupMeans, \(T = \operatorname{scatter}(X)\) and \(B = \operatorname{between}(G, X)\). Then
The scatter of the observations about their own group means, obtained as total scatter less between-group scatter. It must be invertible for the discriminants.
Intuition The noise that group membership does not explain: how much points scatter around their own centroid. Total scatter splits exactly into this part plus between, just as in one-way ANOVA.
Use it when you need the pooled shape of the groups: divided by \(n-k\) it is the pooled within-group covariance that a linear discriminant classifier assumes all groups share. It is the yardstick against which between is judged. It must be invertible for discriminants, which needs more observations than groups plus variables and no collinear columns.
comparing-groupsclassification
def within(g, x) ≔ scatter(x) − between(g, x)
discriminants
discriminants(g, x) · g one-hot group membership n×k · x data matrix n×p → discriminant directions p×p
Let \(G\) (the argument g) and \(X\) (the argument x) be as in groupMeans, with \(W = \operatorname{within}(G, X)\) invertible and \(B = \operatorname{between}(G, X)\). Then
Directions of greatest group separation, one per column, ordered by separations. The first column is
the leading eigenvector of \(W^{-1}B\) — Fisher’s discriminant direction, and the one to use for
two groups. Note that \(W^{-1}B\) is not symmetric and the eigenvector routine is the symmetric QR
algorithm: the later columns are an orthonormal completion of the first, not further eigenvectors
of \(W^{-1}B\) in general. At most \(\min(k-1, p)\) directions carry any separation.
Intuition PCA looks for the direction of most spread and ignores the labels. LDA looks for the direction along which group centroids lie far apart relative to how wide each group is, so projecting onto it makes the groups overlap as little as a straight line allows.
Use it when you have labelled groups and want a supervised projection: to classify (project and threshold, for two groups), to plot labelled data in the few dimensions that separate it, or to see which variables drive the difference (the entries of the direction). Rely on the first column; as the description notes, later columns are not further discriminants in general. It assumes all groups share one covariance shape; groups with very different spreads need a quadratic rule, which this module does not provide.
classificationdimensionality-reductioncomparing-groups
def discriminants(g, x) ≔ eigenvectors(mul(inverse(within(g, x)), between(g, x)))
separations
separations(g, x) · g one-hot group membership n×k · x data matrix n×p → separation per direction vector p
Let \(G\), \(X\), \(W\) and \(B\) be as in discriminants, \(M = \operatorname{inverse}(W)\, B\) and \(V = \operatorname{eigenvectors}(M)\). Then, for \(j = 1, \dots, p\),
The eigenvalues of \(W^{-1}B\) in descending order: the ratio of between-group to within-group variance achieved along each discriminant direction. All but at most \(k-1\) of them are zero.
Intuition A signal-to-noise ratio for group separation: along a direction, how much the centroids spread apart compared with how much each group spreads internally. A value of 10 means the between-group variation is ten times the within-group noise along that line; a value near 0 means the groups overlap completely there.
Use it when you judge whether the groups are separable at all, or how many discriminant directions are worth keeping (count the clearly non-zero values, at most \(k-1\)). The leading value is the one that matches the leading column of discriminants. It is descriptive, not a test: for a p-value on a single variable use pAnova in Inference.
comparing-groupsclassificationeffect-size
def separations(g, x) ≔ eigenvalues(mul(inverse(within(g, x)), between(g, x)))
Building blocks
groupApart(g, x)— Let \(G\) (the argumentg) and \(X\) (the argumentx) be as ingroupMeans, with group means \(\bar{x}_c\) and column means \(\bar{x}\). Then \(\operatorname{groupApart}(G, X)_{cj} = \bar{x}_{cj} - \bar{x}_j\). Each group mean less the overall column means, k×p; used bybetween.scatter(x)— Let \(X \in \mathbb{R}^{n \times p}\) (the argumentx) with rows \(x_i^{\top}\) and column means \(\bar{x}\). Then \(\operatorname{scatter}(X) = T = \sum_{i=1}^{n} (x_i - \bar{x})(x_i - \bar{x})^{\top} = (n - 1)\, S\). The total scatter matrix of the centred data, p×p; used bywithin.
Examples
All examples below are compiled from one file, Fits.hb, so the C# class is Fits. Inputs are
passed through the Tensor overloads with a (Tensor) cast, because most of these functions change
rank (a matrix in, a vector or scalar out).
Fit a straight line
Four points close to \(y = 1 + 2x\). The design matrix has a column of ones for the intercept.
open Prelude.Multivariate
fn lineFit(x, y) ↦ coefficients(x, y)
fn fitQuality(x, y) ↦ rsquared(x, y)
fn leftOver(x, y) ↦ residuals(x, y)
using MLambda.Hilbert.Runtime;
using Examples;
Tensor x = (Tensor)new double[,] { { 1, 1 }, { 1, 2 }, { 1, 3 }, { 1, 4 } };
Tensor y = (Tensor)new[] { 3.1, 4.9, 7.2, 8.8 };
Tensor beta = Fits.LineFit(x, y);
Tensor r2 = Fits.FitQuality(x, y);
Tensor resid = Fits.LeftOver(x, y);
Output (verified, rounded to six significant digits):
coefficients = [1.15, 1.94]
rsquared = 0.995661
residuals = [0.01, -0.13, 0.23, -0.11]
By hand: the slope is \(S_{xy}/S_{xx} = 9.7/5 = 1.94\) and the intercept \(6 - 1.94 \times 2.5 = 1.15\).
Standard errors and t-statistics
The same data as above.
open Prelude.Multivariate
fn standardErrors(x, y) ↦ errors(x, y)
fn tStatistics(x, y) ↦ statistics(x, y)
Tensor se = Fits.StandardErrors(x, y);
Tensor ts = Fits.TStatistics(x, y);
Output (verified, rounded to six significant digits):
errors = [0.247992, 0.090554]
statistics = [4.637248, 21.423716]
The residual sum of squares is 0.082, so \(\hat\sigma^2 = 0.082/2 = 0.041\); the diagonal of \((X^\top X)^{-1}\) is [1.5, 0.2], giving \(\sqrt{0.0615} \approx 0.2480\) and \(\sqrt{0.0082} \approx 0.0906\).
Down-weight an outlier, or shrink with ridge
The last observation is an outlier; a zero weight removes it and the exact line \(y = 1 + 2x\) is recovered. Ridge with \(\lambda = 1\) on the first data set shrinks the coefficients.
open Prelude.Multivariate
fn weightedFit(w, x, y) ↦ weightedCoefficients(w, x, y)
fn shrunk(l, x, y) ↦ ridge(l, x, y)
Tensor w = (Tensor)new[] { 1d, 1d, 1d, 0d };
Tensor yo = (Tensor)new[] { 3d, 5d, 7d, 100d };
Tensor robust = Fits.WeightedFit(w, x, yo);
Tensor ridged = Fits.Shrunk(1.0, x, y);
Output (verified, rounded to six significant digits):
weightedCoefficients = [1, 2]
ridge = [0.854545, 1.972727]
Principal components and whitening
Four points on the axes, spread twice as far along the first variable as the second.
open Prelude.Multivariate
fn varianceAlong(x) ↦ explained(x)
fn varianceShares(x) ↦ explainedPart(x)
fn projected(x) ↦ scores(x)
fn whiteCovariance(x) ↦ covariances(whitened(x))
Tensor d = (Tensor)new double[,] { { 2, 0 }, { 0, 1 }, { -2, 0 }, { 0, -1 } };
Tensor variances = Fits.VarianceAlong(d);
Tensor shares = Fits.VarianceShares(d);
Tensor sc = Fits.Projected(d);
Tensor white = Fits.WhiteCovariance(d);
Output (verified, rounded to six significant digits):
explained = [2.666667, 0.666667]
explainedPart = [0.8, 0.2]
scores = shape [4, 2] values [2, 0, 0, 1, -2, 0, 0, -1]
whitened cov = shape [2, 2] values [1, 0, 0, 1]
The column variances are \(8/3\) and \(2/3\) (divisor \(n-1 = 3\)), so the first component explains 80% of the variance; the data already lie on the principal axes, so the scores equal the data, and whitening gives the identity covariance.
Mahalanobis distances
The same four points.
open Prelude.Multivariate
fn squaredDistances(x) ↦ mahalanobis(x)
fn distances(x) ↦ apart(x)
Tensor d2 = Fits.SquaredDistances(d);
Tensor dm = Fits.Distances(d);
Output (verified, rounded to six significant digits):
mahalanobis = [1.5, 1.5, 1.5, 1.5]
apart = [1.224745, 1.224745, 1.224745, 1.224745]
Every point is equally far in covariance units: \(2^2 / (8/3) = 1.5\) and \(1^2 / (2/3) = 1.5\).
Separate two groups
Six observations of two variables in two groups of three.
open Prelude.Multivariate
fn groupCentres(g, x) ↦ groupMeans(g, x)
fn withinScatter(g, x) ↦ within(g, x)
fn betweenScatter(g, x) ↦ between(g, x)
fn separationPower(g, x) ↦ separations(g, x)
Tensor g = (Tensor)new double[,] { { 1, 0 }, { 1, 0 }, { 1, 0 }, { 0, 1 }, { 0, 1 }, { 0, 1 } };
Tensor gx = (Tensor)new double[,] { { 1, 2 }, { 2, 1 }, { 3, 3 }, { 6, 5 }, { 7, 7 }, { 8, 6 } };
Tensor centres = Fits.GroupCentres(g, gx);
Tensor wScat = Fits.WithinScatter(g, gx);
Tensor bScat = Fits.BetweenScatter(g, gx);
Tensor sep = Fits.SeparationPower(g, gx);
Output (verified):
groupMeans = shape [2, 2] values [2, 2, 7, 6]
within = shape [2, 2] values [4, 2, 2, 4]
between = shape [2, 2] values [37.5, 30, 30, 24]
separations = [10.5, 0]
With two groups only one direction separates them, so the second separation is zero. The leading
column of discriminants for these data is \((0.894427, 0.447214)\) (verified), the eigenvector of
\(W^{-1}B\) for 10.5.