Selection
Cross-validation folds as masks, a cross-validated ridge regression with a penalty grid search, and simple column scores for feature selection.
This module covers model selection: splitting rows into folds, estimating out-of-sample error by k-fold cross-validation, choosing a regularisation penalty from a grid, and scoring columns to decide which features to keep.
Cross-validation here moves no data. A fold assignment is a vector of numbers, so a fold’s test set is the mask “fold equals f” and its training set is the complementary mask. Fitting on the training rows is an ordinary weighted fit in which the held-out rows carry a weight of zero, so no subset of the data is ever copied.
Conventions. The design matrix x has shape n×p, one row per observation; include a column of ones yourself if you want an intercept. The response y is a vector of length n. The number of folds k, a fold number f (counting from 0) and penalties l are reals. Folds are assigned by position, in contiguous blocks, never at random, so a cross-validated score is reproducible without a seed; if the rows should be shuffled, shuffle them before the call (a permutation can be drawn as in Prelude.Resampling).
Limitations. The language has no function-valued parameters, so cross-validation cannot take an arbitrary model as an argument. The fold masks and fold arithmetic are general, but the cross-validated estimator is ridge regression specifically; to cross-validate a different model, write your own versions of foldError, cvLoop and crossError with that model’s fit in place of heldRidge. The ridge penalty applies to every column, including an intercept column if you supply one. The linear solve is conjugate gradient, which requires the penalised normal equations to be symmetric positive-definite; a positive penalty guarantees this.
open Prelude.Selection
The folds
With \(n\) rows and \(k\) folds, row \(i\) (counting from 0) belongs to fold
\[f_i = \left\lfloor \frac{i\,k}{n} \right\rfloor\]which splits the rows into \(k\) contiguous blocks of nearly equal size. Stratified folds instead number the rows within each class and deal them out in turn, \(f_i = r_i \bmod k\) where \(r_i\) counts the earlier rows of the same class, so every fold receives its share of a rare class.
foldIds
foldIds(x, k) · x design matrix n×p · k number of folds scalar → fold numbers vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix, with rows at positions \(i = 0, 1, \dots, n-1\), and let \(k \in \{1, 2, \dots\}\) be the number of folds. Then
The block fold number of every row, from 0 to \(k - 1\). Only the number of rows of x matters.
Intuition Cut the rows into k consecutive slices of nearly equal size and label each row with its slice. Each slice will take one turn as the unseen test set while the others train, so every row is tested exactly once.
Use it when you set up k-fold cross-validation; 5 or 10 folds is the usual choice, and \(k = n\) gives leave-one-out. Because the split follows row order, data sorted by time, class or response give unrepresentative folds: shuffle the rows first unless you deliberately want contiguous blocks (as a rough guard for time-ordered data). For a rare class, use stratified.
cross-validationmodel-selection
def foldIds(x, k) ≔ floor(rowIds(x) · k / rows(x))
inFold
inFold(x, k, f) · x design matrix n×p · k number of folds scalar · f fold number scalar → 0/1 mask vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix, with rows at positions \(i = 0, 1, \dots, n-1\), and let \(k \in \{1, 2, \dots\}\) be the number of folds, and let \(f \in \{0, \dots, k-1\}\) be a fold number. Then
The test mask of fold f: 1 for the rows in the fold.
Intuition A test set that is a weight vector rather than a copy of rows: multiply a per-row error by it and only fold f counts. Selecting data becomes arithmetic, so no subset is ever built.
Use it when you score a model on one fold’s held-out rows, for example summing squared residuals times this mask as foldError does, or writing a fold loop for a model other than ridge. Divide by foldSize to turn the masked sum into a mean. For the training rows use outFold.
cross-validationheld-out-evaluation
def inFold(x, k, f) ≔ foldIds(x, k) = f
outFold
outFold(x, k, f) · x design matrix n×p · k number of folds scalar · f fold number scalar → 0/1 mask vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix, with rows at positions \(i = 0, 1, \dots, n-1\), and let \(k \in \{1, 2, \dots\}\) be the number of folds, and let \(f \in \{0, \dots, k-1\}\) be a fold number. Then
The training mask of fold f: 1 for every row outside the fold. Use it as the weight vector of a weighted fit.
Intuition The complement of inFold: every row the model may learn from in this round. Passed as weights, a zero makes a row vanish from the fit exactly, as if it had been deleted.
Use it when you train on everything except fold f, as the weights of heldRidge or weightedCoefficients in Multivariate. The model must never see the test rows during fitting, and that includes any preprocessing learned from data (scaling, feature scores): compute those with the same weights, or the cross-validated error will be optimistic.
cross-validationoverfitting-control
def outFold(x, k, f) ≔ 1 − inFold(x, k, f)
foldSize
foldSize(x, k, f) · x design matrix n×p · k number of folds scalar · f fold number scalar → row count scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix, with rows at positions \(i = 0, 1, \dots, n-1\), and let \(k \in \{1, 2, \dots\}\) be the number of folds, and let \(f \in \{0, \dots, k-1\}\) be a fold number. Then
The number of rows in fold f.
Intuition Summing a 0/1 mask counts its ones. It is the divisor that turns a masked total error into a mean error for that fold.
Use it when you average a fold’s errors in your own cross-validation loop, or check that no fold is empty or tiny. With more folds than rows some folds are empty; foldError guards that division and reports 0, which then pulls crossError down, so keep \(k \le n\).
cross-validationsanity-checking
def foldSize(x, k, f) ≔ einsum("i->", inFold(x, k, f))
stratified
stratified(x, g, k) · x design matrix n×p · g one-hot class membership matrix n×C · k number of folds scalar → fold numbers vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) have rows at positions \(i = 0, \dots, n-1\), let \(g \in \{0, 1\}^{n \times C}\) be a one-hot class membership matrix, and let \(k \in \{1, 2, \dots\}\) be the number of folds. Let \(r_i = \sum_{j=0}^{n-1} \bigl(\sum_{c=1}^{C} g_{ic}\, g_{jc}\bigr)\, [\, i > j \,]\) be the number of earlier rows in the class of row \(i\). Then
Stratified fold numbers: within each class, the rows are numbered in order and assigned to folds 0, 1, …, k − 1, 0, 1, … in turn. This keeps each class’s share roughly equal across folds, which matters when a class is rare. The result is a fold vector like foldIds; note that crossError and the grid search use the block folds, so stratified folds are for a fold loop of your own. It builds n×n comparison matrices.
Intuition Deal each class out like a deck of cards: the first member of a class goes to fold 0, the next to fold 1, and so on, separately within each class. Every fold then gets close to its fair share of every class, however rare.
Use it when you cross-validate a classifier with imbalanced classes, where a block split could leave a fold with no positives and a meaningless score. Rows keep their order within each class, so shuffle beforehand if that order carries structure. crossError and the grid search use block folds, so these fold numbers are for your own fold loop. Memory grows as \(n^2\).
cross-validationclass-imbalance
def stratified(x, g, k) ≔ withinClass(x, g) − floor(withinClass(x, g) / k) · k
Building blocks
rowIds(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rowIds}(X)_i = i\) for \(i = 0, 1, \dots, n-1\). The row positions of a matrix.sameClass(g)— Let \(g \in \{0, 1\}^{n \times C}\). Then \(\operatorname{sameClass}(g)_{ij} = \sum_{c=1}^{C} g_{ic}\, g_{jc}\), which for a one-hot \(g\) is 1 exactly when rows \(i\) and \(j\) share a class. The n×n 0/1 matrix marking pairs of rows in the same class.earlier(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{earlier}(X)_{ij} = [\, i > j \,]\) for positions \(i, j = 0, \dots, n-1\). The n×n 0/1 matrix marking pairs where the second row comes before the first.withinClass(x, g)— Let \(X \in \mathbb{R}^{n \times p}\) and \(g \in \{0, 1\}^{n \times C}\) one-hot. Then \(\operatorname{withinClass}(X, g)_i = \sum_{j} \operatorname{sameClass}(g)_{ij}\, [\, i > j \,]\). How many earlier rows share each row’s class: its position within its class.
A cross-validated ridge
For a penalty \(\lambda\) and a 0/1 training weight \(w\), the ridge coefficients solve the weighted, penalised normal equations
\[\left( X^\top W X + \lambda I \right) \beta = X^\top W y, \qquad W = \operatorname{diag}(w)\]The error of fold \(f\) is the mean squared residual on its own rows, using coefficients fitted on the others, and the cross-validated error is the average over the \(k\) folds:
\[\mathrm{CV}(\lambda) = \frac{1}{k} \sum_{f=0}^{k-1} \frac{1}{\lvert F_f \rvert} \sum_{i \in F_f} \left( y_i - x_i^\top \hat{\beta}^{(-f)} \right)^2\]heldRidge
heldRidge(l, x, y, w) · l penalty scalar · x design matrix n×p · y response vector n · w row weights vector n → coefficients vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), let \(y \in \mathbb{R}^n\) be the response, let \(w \in [0, \infty)^n\) be row weights, and let \(\lambda > 0\) (the argument l) be the penalty. Let \(W = \operatorname{diag}(w)\). Then
where \(X^{\top} W X = \sum_i w_i\, x_i x_i^{\top}\), \(X^{\top} W y = \sum_i w_i\, y_i\, x_i\), and \(\operatorname{solve}(A, b)\) is \(p\) turns of Jacobi-preconditioned conjugate gradient from Prelude.LinearAlgebra, which returns \(A^{-1} b\) in exact arithmetic for a symmetric positive-definite \(A\).
Weighted ridge regression. With a 0/1 weight vector it fits on the rows weighted 1 only; a row weighted 0 contributes nothing, as if it were absent.
Intuition Ridge regression and weighted least squares in one: the weights decide which rows the fit listens to, and the penalty pulls coefficients toward zero. With a 0/1 fold mask it is “ridge on the training rows only” without building the training set.
Use it when you fit ridge on a subset or with unequal row weights: inside a custom cross-validation loop, for down-weighting unreliable rows, or for bootstrap-style reweighting. Keep \(\lambda > 0\) so the system stays positive-definite. The penalty does not scale with the total weight, so the same \(\lambda\) shrinks harder when fewer rows carry weight. With all weights 1 it equals ridge in Multivariate.
regressioncross-validationoverfitting-control
def heldRidge(l, x, y, w) ≔ solve(weightedGram(w, x) + l · identity(gram(x)), weightedMoment(w, x, y))
foldError
foldError(l, x, y, k, f) · l penalty scalar · x design matrix n×p · y response vector n · k number of folds scalar · f fold number scalar → mean squared error scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), let \(y \in \mathbb{R}^n\) be the response, let \(\lambda > 0\) (the argument l) be the penalty, and let \(k \in \{1, 2, \dots\}\) be the number of folds, with rows at positions \(i = 0, \dots, n-1\) and fold sets \(F_f = \{ i : \lfloor i k / n \rfloor = f \}\). Let \(f \in \{0, \dots, k-1\}\), let \(\hat{\beta}^{(-f)} = \operatorname{heldRidge}(\lambda, X, y, \operatorname{outFold}(X, k, f))\) be the ridge fitted without fold \(f\), and let \(r_i = y_i - x_i^{\top} \hat{\beta}^{(-f)}\). Then
The test error of fold f: fit ridge on the other folds, then average the squared residuals over the rows of fold f.
Intuition An honest exam: the model is fitted without these rows and then graded on them, so it cannot score well by memorising them. The mean squared residual on unseen rows estimates how the model will do on new data.
Use it when you want one fold’s score, to see whether the folds disagree (a single bad fold points to unusual rows or an unshuffled order) or to build a fold loop by hand. An empty fold returns 0 instead of failing. For the average over all folds use crossError.
cross-validationheld-out-evaluation
def foldError(l, x, y, k, f) ≔ einsum("i,i,i->", inFold(x, k, f), foldResidual(l, x, y, k, f), foldResidual(l, x, y, k, f)) / safe(foldSize(x, k, f))
crossError
crossError(l, x, y, k) · l penalty scalar · x design matrix n×p · y response vector n · k number of folds scalar → cross-validated mean squared error scalar
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix with rows \(x_i^{\top}\), let \(y \in \mathbb{R}^n\) be the response, let \(\lambda > 0\) (the argument l) be the penalty, and let \(k \in \{1, 2, \dots\}\) be the number of folds. Then
The k-fold cross-validated mean squared error of ridge regression with penalty l: the average of the k fold errors. It performs k linear solves.
Intuition In-sample error rewards memorising noise, so it keeps falling as a model gets more flexible. Averaging the error over folds that each played “unseen data” once estimates prediction error on genuinely new rows, and it rises again when a model overfits.
Use it when you compare ridge settings or feature sets by how well they predict, rather than by rsquared. It averages the fold means equally, so unequal fold sizes weigh a row in a small fold a little more. The estimate is itself noisy on small data; differences smaller than the spread of foldError across folds are not meaningful. To search many penalties use gridErrors.
cross-validationmodel-selectionoverfitting-control
def crossError(l, x, y, k) ≔ cvLoop(k, l, x, y, k, 0) / k
Building blocks
foldResidual(l, x, y, k, f)— Let \(X, y, \lambda, k\) and \(F_f\) be as infoldError. Then \(\operatorname{foldResidual}(\lambda, X, y, k, f) = y - X \hat{\beta}^{(-f)}\), where \(\hat{\beta}^{(-f)} = \operatorname{heldRidge}(\lambda, X, y, \operatorname{outFold}(X, k, f))\). The residuals on all rows of the ridge fitted without fold f.cvLoop(f, l, x, y, k, run)— Let \(X, y, \lambda, k\) be as incrossError, \(f \in \mathbb{N}\) and \(\mathit{run} \in \mathbb{R}\). Then \(\operatorname{cvLoop}(0, \lambda, X, y, k, \mathit{run}) = \mathit{run}\) and \(\operatorname{cvLoop}(f, \lambda, X, y, k, \mathit{run}) = \operatorname{cvLoop}(f - 1, \lambda, X, y, k, \mathit{run} + \operatorname{foldError}(\lambda, X, y, k, f - 1))\), so \(\operatorname{cvLoop}(f, \dots, \mathit{run}) = \mathit{run} + \sum_{g=0}^{f-1} \operatorname{foldError}(\lambda, X, y, k, g)\). The loop that adds up the fold errors from fold f − 1 down to fold 0; used bycrossError.
The grid
A grid search evaluates the cross-validated error at each penalty in a vector ls and picks the one with the smallest error. The best penalty is found with a mask, not an index: the errors are compared with their own minimum and the penalty is read out of the matching positions. If several penalties tie for the minimum, the result is their mean.
gridErrors
gridErrors(ls, x, y, k) · ls candidate penalties vector m · x design matrix n×p · y response vector n · k number of folds scalar → cross-validated errors vector m
Let \(\ell = (\ell_0, \dots, \ell_{m-1}) \in (0, \infty)^m\) (the argument ls) be candidate penalties at positions \(0, \dots, m-1\), let \(X \in \mathbb{R}^{n \times p}\) (the argument x) be the design matrix, \(y \in \mathbb{R}^n\) the response, and \(k \in \{1, 2, \dots\}\) the number of folds. Then
The cross-validated error of every candidate penalty, in the order of ls. It performs \(m \cdot k\) linear solves.
Intuition The validation curve: held-out error as a function of the penalty. Too little penalty overfits and too much underfits, so the curve typically dips to a minimum somewhere in between.
Use it when you want to see that curve rather than only its minimum, to check the choice is not at the edge of the grid or that a wide range of penalties is nearly as good. Space candidates on a log scale (for example 0.001, 0.01, …, 1000), since the useful penalty can span orders of magnitude. The cost is \(m \cdot k\) solves; for just the winning value use bestPenalty.
hyperparameter-choicecross-validationmodel-selection
def gridErrors(ls, x, y, k) ≔ gridLoop(count(ls), ls, positions(ls), x, y, k, zeroed(ls))
bestPenalty
bestPenalty(ls, x, y, k) · ls candidate penalties vector m · x design matrix n×p · y response vector n · k number of folds scalar → penalty scalar
Let \(\ell \in (0, \infty)^m\), \(X\), \(y\) and \(k\) be as in gridErrors, with errors \(E = \operatorname{gridErrors}(\ell, X, y, k)\), and let \(b_q = [\, E_q = \min_{q'} E_{q'} \,]\) be the mask of the best candidates. Then
the mean of the penalties that attain the minimum error.
The candidate penalty with the lowest cross-validated error; the mean of the tied candidates if there is a tie.
Intuition Pick the penalty that predicted held-out rows best. It is found by comparing every error with the minimum, so exact ties average their penalties instead of favouring whichever came first.
Use it when you need a single penalty to report or reuse. If it lands on the smallest or largest candidate, the real optimum may lie outside the grid, so extend it. A tie averages the tied values, which on a log grid need not be a grid point. It recomputes the whole grid; when you also want the curve, call gridErrors once and pick the minimum yourself. To fit with the chosen value directly, use tunedRidge.
hyperparameter-choicemodel-selectionoverfitting-control
def bestPenalty(ls, x, y, k) ≔ einsum("i,i->", ls, gridBest(ls, x, y, k)) / safe(einsum("i->", gridBest(ls, x, y, k)))
tunedRidge
tunedRidge(ls, x, y, k) · ls candidate penalties vector m · x design matrix n×p · y response vector n · k number of folds scalar → coefficients vector p
Let \(\ell \in (0, \infty)^m\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^n\) and \(k\) be as in gridErrors, and let \(\lambda^* = \operatorname{bestPenalty}(\ell, X, y, k)\). Then
where \(\operatorname{solve}(A, b)\) is \(p\) turns of Jacobi-preconditioned conjugate gradient from Prelude.LinearAlgebra, which returns \(A^{-1} b\) in exact arithmetic for a symmetric positive-definite \(A\).
Ridge regression refitted on all the rows with the penalty chosen by bestPenalty.
Intuition Cross-validation is only for choosing the setting; once it is chosen, every row should help estimate the coefficients. So the final model is refitted on all the data with the winning penalty.
Use it when you want a regularised linear model with its penalty picked automatically from a candidate list. Its cross-validated error is a slightly optimistic estimate of its performance, since the same folds chose the penalty; for an unbiased figure, score it on data kept out of the whole procedure. Each fold fits on about \((k-1)/k\) of the rows, so the same penalty acts slightly less strongly in the refit. Standardise columns first.
regressionoverfitting-controlhyperparameter-choice
def tunedRidge(ls, x, y, k) ≔ ridge(bestPenalty(ls, x, y, k), x, y)
Building blocks
gridAt(ls, ids, j)— Let \(\ell \in \mathbb{R}^m\) (the argumentls), \(\mathit{ids} \in \mathbb{R}^m\) and \(j \in \mathbb{R}\). Then \(\operatorname{gridAt}(\ell, \mathit{ids}, j) = \sum_{q} \ell_q\, [\, \mathit{ids}_q = j \,]\), which is \(\ell_j\) when \(\mathit{ids} = (0, 1, \dots, m-1)\). The j-th candidate penalty, read out with a mask.gridLoop(j, ls, ids, x, y, k, acc)— Let \(\ell, \mathit{ids} \in \mathbb{R}^m\), \(X, y, k\) as ingridErrors, \(j \in \mathbb{N}\) and \(\mathit{acc} \in \mathbb{R}^m\). Then \(\operatorname{gridLoop}(0, \dots, \mathit{acc}) = \mathit{acc}\) and \(\operatorname{gridLoop}(j, \ell, \mathit{ids}, X, y, k, \mathit{acc}) = \operatorname{gridLoop}\bigl(j - 1, \ell, \mathit{ids}, X, y, k, \mathit{acc} + [\, \mathit{ids} = j - 1 \,] \operatorname{crossError}(\operatorname{gridAt}(\ell, \mathit{ids}, j - 1), X, y, k)\bigr)\), the comparison taken entrywise. The loop that fills the error vector one penalty at a time; used bygridErrors.gridBest(ls, x, y, k)— Let \(\ell, X, y, k\) be as ingridErrors, with errors \(E\). Then \(\operatorname{gridBest}(\ell, X, y, k)_q = [\, E_q = \min_{q'} E_{q'} \,]\). The 0/1 mask of the penalties whose error equals the minimum.
Choosing features
Feature selection here is a score per column followed by a mask over the scores. The variance says which columns vary at all; the alignment score says which vary with the response:
\[\text{align}_j = \frac{\left\lvert \sum_i (x_{ij} - \bar{x}_j)(y_i - \bar{y}) \right\rvert}{n\, s_j}\]where \(s_j\) is the sample standard deviation of column \(j\). This equals \(\tfrac{n-1}{n}\, \lvert r_j \rvert\, s_y\), so for a fixed response it ranks the columns exactly as the absolute correlation does, although its values are in the units of y. A constant column scores 0.
spreadScore
spreadScore(x) · x data matrix n×p → column variances vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with column means \(\bar{x}_j = \frac{1}{n} \sum_i X_{ij}\). Then
The sample variance of each column. Columns with (near) zero variance carry no information.
Intuition A column that barely changes cannot tell rows apart, so it cannot help predict anything. Variance measures how much a column moves, ignoring the response entirely.
Use it when you remove constant or nearly constant columns before fitting, a cheap cleanup that also avoids singular normal equations. Variance is in squared units, so a threshold only makes sense among columns on the same scale; do not rank mixed-unit columns by it. Remember that an intercept column of ones has variance 0: exclude it from the filter. To judge relevance to the response, use alignScore.
feature-selectiondata-cleaning
def spreadScore(x) ≔ variances(x)
alignScore
alignScore(x, y) · x data matrix n×p · y response vector n → column scores vector p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) with column means \(\bar{x}_j\) and sample standard deviations \(s_j\), and let \(y \in \mathbb{R}^n\) with mean \(\bar{y}\). Then
How strongly each column moves with the response, as defined above: proportional to the absolute Pearson correlation.
Intuition A one-column-at-a-time relevance check: does this column rise and fall with the response? The score is the column’s covariance with y divided by its own spread, so column units cancel out, and it ranks columns exactly as \(\lvert r \rvert\) does.
Use it when you want a fast filter over many candidate features before a model. It sees only linear, individual effects: it misses features that matter only in combination and keeps redundant copies of the same signal. Selecting features on all rows and then cross-validating leaks information and flatters the score, so compute it on training rows only. Pick columns with keepBest or keepAbove.
feature-selectioncorrelation-analysisdimensionality-reduction
def alignScore(x, y) ≔ abs(einsum("ij,i->j", centred(x), y − mean(y))) / safe(sqrt(max(variances(x), 0)) · count(y))
keepAbove
keepAbove(s, t) · s scores vector p · t threshold scalar → 0/1 mask vector p
Let \(s \in \mathbb{R}^p\) be scores and \(t \in \mathbb{R}\) a threshold. Then
\[\operatorname{keepAbove}(s, t)_j = [\, s_j \ge t \,], \qquad j = 1, \dots, p .\]Marks the columns whose score is at least t.
Intuition A pass mark: every column that scores at least t is kept, however many that turns out to be. The number of features adapts to the data.
Use it when the threshold has a meaning of its own, such as a minimum variance from spreadScore or a minimum relevance from alignScore. Scores must be on a comparable scale for one threshold to be fair. When you need a fixed number of features instead, use keepBest. Apply the mask with keptColumns.
feature-selectiondata-cleaning
def keepAbove(s, t) ≔ s ≥ t
keepBest
keepBest(s, k) · s scores vector p · k number of columns to keep scalar → 0/1 mask vector p
Let \(s \in \mathbb{R}^p\) be scores and \(k \in \{0, \dots, p\}\), and let \(\rho_j = \tfrac{1}{2}\bigl(\#\{l : s_l < s_j\} + \#\{l : s_l \le s_j\} - 1\bigr)\) be the zero-based midrank of \(s_j\). Then
\[\operatorname{keepBest}(s, k)_j = [\, \rho_j \ge p - k \,], \qquad j = 1, \dots, p .\]Marks the k columns with the highest scores. Scores are ranked with midranks, so ties at the boundary can mark slightly fewer or more than k.
Intuition Top-k by rank: keep the k highest-scoring columns, whatever their absolute scores. Only the order of the scores matters, so it works equally for any score scale.
Use it when you have a budget of features (for speed, interpretability, or \(p\) well below \(n\)), and choose k itself by cross-validation rather than by eye. Check the number of ones in the mask when scores can tie, since midranks can keep slightly more or fewer than k. For a score cut-off instead of a count, use keepAbove.
feature-selectiondimensionality-reduction
def keepBest(s, k) ≔ rank(s) ≥ count(s) − k
keptColumns
keptColumns(x, m) · x data matrix n×p · m column mask vector p → masked data matrix n×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) and let \(m \in \{0, 1\}^p\) be a column mask. Then
The data with the unselected columns set to zero. The shape is unchanged; the columns are zeroed, not removed.
Intuition Feature selection as a mask rather than a copy: dropped columns become zeros, so the matrix keeps its shape and the kept columns keep their positions.
Use it when you try a subset of features with an existing fit without rebuilding the design matrix. A zero column still counts as a column: with ridge or heldRidge its coefficient is simply 0, but with plain coefficients in Multivariate it makes \(X^\top X\) singular, and the degrees of freedom in adjusted and spread still count it. Remove the columns outside Hilbert when those matter.
feature-selectionregression
def keptColumns(x, m) ≔ x · einsum("i,j->ij", rowOnes(x), m)
Examples
Several of these functions change the rank of their input (a matrix in, a vector out). Call them with (Tensor) casts on the arguments, as shown, so the general Tensor overload is used.
Block and stratified folds
Six rows, three block folds; then two stratified folds for four rows of one class and two of another.
open Prelude.Selection
fn foldOf(x, k) ↦ foldIds(x, k)
fn testMask(x, k, f) ↦ inFold(x, k, f)
fn stratifiedFolds(x, g, k) ↦ stratified(x, g, k)
using MLambda.Hilbert.Runtime;
using Examples;
var x = new double[,] { {1,0},{1,1},{1,2},{1,3},{1,4},{1,5} };
var folds = Tuning.FoldOf((Tensor)x, 3.0);
var test1 = Tuning.TestMask((Tensor)x, 3.0, 1.0);
var g = new double[,] { {1,0},{1,0},{1,0},{1,0},{0,1},{0,1} };
var strat = Tuning.StratifiedFolds((Tensor)x, (Tensor)g, 2.0);
Output (verified):
folds = [0, 0, 1, 1, 2, 2]
test1 = [0, 0, 1, 1, 0, 0]
strat = [0, 1, 0, 1, 0, 1]
The first class’s rows (0 to 3) alternate between folds 0 and 1, and so do the second class’s (4 and 5), so each fold receives two rows of the first class and one of the second.
Cross-validated error of one ridge penalty
A straight-line fit with an intercept column, three folds, penalty 0.1.
open Prelude.Selection
fn heldOutError(l, x, y, k) ↦ crossError(l, x, y, k)
using Examples;
var x = new double[,] { {1,0},{1,1},{1,2},{1,3},{1,4},{1,5} };
var y = new[] { 1.1d, 2.9, 5.2, 6.8, 9.1, 11.0 };
var cv = Tuning.HeldOutError(0.1, x, y, 3.0);
Output (verified, rounded to 6 decimal places):
cv = 0.042582
This matches a direct computation: fit ridge on rows 2 to 5 and score rows 0 and 1, and likewise for the other two folds, then average the three mean squared errors.
Choose a penalty from a grid and refit
open Prelude.Selection
fn penaltyErrors(ls, x, y, k) ↦ gridErrors(ls, x, y, k)
fn chosenPenalty(ls, x, y, k) ↦ bestPenalty(ls, x, y, k)
fn tunedFit(ls, x, y, k) ↦ tunedRidge(ls, x, y, k)
using Examples;
var x = new double[,] { {1,0},{1,1},{1,2},{1,3},{1,4},{1,5} };
var y = new[] { 1.1d, 2.9, 5.2, 6.8, 9.1, 11.0 };
var ls = new[] { 0.01d, 1, 100 };
var grid = Tuning.PenaltyErrors(ls, x, y, 3.0);
var best = Tuning.ChosenPenalty(ls, x, y, 3.0);
var fit = Tuning.TunedFit(ls, x, y, 3.0);
Output (verified, rounded to 6 decimal places):
grid = [0.038085, 0.226962, 31.660899]
best = 0.01
fit = [1.035516, 1.99177]
The data are nearly a line, so the smallest penalty wins; a penalty of 100 shrinks both coefficients, the intercept included, towards zero. The refitted line is roughly \(y = 1.036 + 1.992\,x\).
Score columns and keep the best two
The first column follows the response, the second is constant, and the third moves only slightly with it.
open Prelude.Selection
fn movingColumns(x) ↦ spreadScore(x)
fn relevance(x, y) ↦ alignScore(x, y)
fn topColumns(s, k) ↦ keepBest(s, k)
fn masked(x, m) ↦ keptColumns(x, m)
using MLambda.Hilbert.Runtime;
using Examples;
var z = new double[,] { {1,5,2},{2,5,1},{3,5,2},{4,5,1} };
var t = new[] { 2d, 4, 6, 8 };
var spread = Tuning.MovingColumns((Tensor)z);
var align = Tuning.Relevance((Tensor)z, (Tensor)t);
var top = Tuning.TopColumns(new[] { 0.2d, 0.9, 0.5 }, 2.0);
var kept = Tuning.Masked((Tensor)z, (Tensor)new[] { 1d, 0, 1 });
Output (verified, rounded to 6 decimal places):
spread = [1.666667, 0, 0.333333]
align = [1.936492, 0, 0.866025]
top = [0, 1, 1]
masked = shape [4, 3] values [1, 0, 2, 2, 0, 1, 3, 0, 2, 4, 0, 1]
For the first column the centred cross-products sum to 10 and the standard deviation is 1.291, so the score is \(10 / (4 \cdot 1.291)\). keepBest with scores 0.2, 0.9 and 0.5 keeps the second and third. The masked matrix keeps its shape, with the second column zeroed.