Trees
Regression and binary classification trees (CART), weighted trees, random forests from bootstrap counts, gradient boosting and penalised leaves.
This module grows decision trees without storing one. A tree is described by a membership matrix — which training row sits in which node — and every quantity the growth needs (node sizes, sums, the gain of every candidate split) is a contraction against that matrix. The candidate thresholds are the data themselves: for feature f, the cut points tried are the values that column already takes, and row i goes left under candidate t when \(x_{if} \le x_{tf}\). Splits minimise the squared error of the children; for 0/1 labels this is the Gini criterion up to a factor, so a classification tree is the same regression tree read at one half.
Conventions. The data matrix x is n×p — one row per sample, and a matrix even for a single feature. Targets y are a vector of length n; for classification they must be 0 or 1. Queries z are m×p. The depth d is the number of levels of splitting. A node table never has more useful slots than the data has rows, so nodes are numbered from zero within each level (the children of slot j are slots 2j and 2j+1 of the next level) and the usable depth is bounded by \(2^d \le n\). Among equally good splits the first is kept, ordered by row of the threshold and then by feature. A node whose best gain is no more than \(10^{-12}\) does not split and passes all its rows to its left child, so a leaf simply repeats down the remaining levels.
Limitations. Growth costs about \(n^3 p\) per level, so trees suit thousands of rows, not millions; subsample beyond that. Because no tree is stored, predicting a new row re-grows the tree from the training data, which also guarantees a prediction can never disagree with the fit. Multi-class trees, the isolation forest and cost-complexity pruning are not provided; for outlier detection see the local outlier factor in Neighbours and the Mahalanobis distance in Multivariate.
Notation. Throughout this page every index is counted from zero: data rows and threshold rows \(i, t \in \{0, 1, \dots, n-1\}\), node slots \(j, c \in \{0, 1, \dots, n-1\}\), features \(f \in \{0, 1, \dots, p-1\}\) and query rows \(i \in \{0, 1, \dots, m-1\}\). \(X \in \mathbb{R}^{n \times p}\) is the data (the argument x), \(y \in \mathbb{R}^{n}\) the targets, \(Z \in \mathbb{R}^{m \times p}\) the queries, and a membership matrix \(M \in \mathbb{R}^{n \times n}\) (the argument m) holds in \(M_{ij}\) the weight with which row \(i\) sits in slot \(j\). \(e_0\) is the basis vector of slot \(0\), \([\,\cdot\,]\) is the Iverson bracket, and \(\operatorname{safe}(N) = N + [N = 0]\) entrywise.
open Prelude.Trees
The candidate grid
Every row of the data is a candidate threshold for every feature, so “does row i go left under candidate (t, f)” is one n×n×p mask. Per-node and per-candidate quantities are placed on their axes by explicit outer products rather than by broadcasting.
Building blocks
sends(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{sends}(X)_{itf} = [\,X_{if} \le X_{tf}\,]\). The n×n×p mask, one where row i goes left under threshold row t on feature f.overRows(x, v)— Let \(X \in \mathbb{R}^{n \times p}\) and \(v \in \mathbb{R}^{k}\) (usually \(k = n\), one value per slot). Then \(\operatorname{overRows}(X, v) = \mathbf{1}_n v^{\top}\), that is \(\operatorname{overRows}(X, v)_{ij} = v_j\). Repeats a per-node vector across rows, giving an n×n matrix.overNodes(x, c)— Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{n \times p}\) (the argumentc). Then \(\operatorname{overNodes}(X, C)_{jtf} = C_{tf}\) for every slot \(j\). Repeats a per-candidate n×p matrix across node slots, giving n×n×p.perNode(x, v)— Let \(X \in \mathbb{R}^{n \times p}\) and \(v \in \mathbb{R}^{n}\). Then \(\operatorname{perNode}(X, v)_{jtf} = v_j\) for every candidate \((t, f)\). Repeats a per-node vector across all candidates, giving n×n×p.
The gain of every candidate at every node
For a node with membership weights \(m_i\), a side holding mass \(N\) and target sum \(S\) explains \(S^2/N\) of the sum of squares. The gain of a split is what the two children explain less what the parent does:
\[\operatorname{gain} = \frac{S_L^2}{N_L} + \frac{S_R^2}{N_R} - \frac{S^2}{N}\]No variance is ever formed. A candidate that would leave either child empty has gain zero.
splitGain
splitGain(m, x, y) · m membership, row i in node slot j n×n · x data n×p · y targets vector n → gain of every candidate at every node n×n×p
Let \(M \in \mathbb{R}^{n \times n}\) (the argument m) be a membership matrix with non-negative weights, \(X \in \mathbb{R}^{n \times p}\) (the argument x) the data and \(y \in \mathbb{R}^{n}\) (the argument y) the targets. Then
where
\[N_j = \sum_{i} M_{ij}, \quad S_j = \sum_{i} M_{ij} y_i, \quad N^L_{jtf} = \sum_{i} M_{ij} [\,X_{if} \le X_{tf}\,], \quad S^L_{jtf} = \sum_{i} M_{ij} [\,X_{if} \le X_{tf}\,]\, y_i ,\]\(N^R_{jtf} = N_j - N^L_{jtf}\) and \(S^R_{jtf} = S_j - S^L_{jtf}\).
The squared-error reduction above for node slot j, threshold row t and feature f, masked to zero where a child would be empty.
Intuition A side holding mass \(N\) with target sum \(S\) can be summarised by its mean, and doing so explains \(S^2/N\) of the sum of squares. A good split creates two sides whose means differ, so together they explain more than the parent’s single mean. The gain is that improvement, scored for every node, threshold and feature at once.
Use it when you want to see why a tree chose a split, rank features by their best gain at a node, or build your own stopping rule. Pass membership(k, x, y) for the nodes at level k, where membership(0, x, y) is the root. The result is n×n×p, so memory grows as \(n^2 p\).
split-selectionfeature-importanceinterpretable-models
def splitGain(m, x, y) ≔ (childScore(m, x, y) − perNode(x, sideScore(nodeSum(m, y), nodeMass(m)))) · viable(m, x)
Building blocks
nodeMass(m)— Let \(M \in \mathbb{R}^{n \times n}\). Then \(\operatorname{nodeMass}(M)_j = N_j = \sum_{i} M_{ij}\). The number (or total weight) of rows in each node slot.nodeSum(m, y)— Let \(M \in \mathbb{R}^{n \times n}\) and \(y \in \mathbb{R}^{n}\). Then \(\operatorname{nodeSum}(M, y)_j = S_j = \sum_{i} M_{ij} y_i\). The sum of the targets in each node slot.leftMass(m, x)— Let \(M \in \mathbb{R}^{n \times n}\) and \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{leftMass}(M, X)_{jtf} = N^L_{jtf} = \sum_{i} M_{ij} [\,X_{if} \le X_{tf}\,]\). The mass each candidate sends left, per node.leftSum(m, x, y)— Let \(M \in \mathbb{R}^{n \times n}\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\). Then \(\operatorname{leftSum}(M, X, y)_{jtf} = S^L_{jtf} = \sum_{i} M_{ij} [\,X_{if} \le X_{tf}\,]\, y_i\). The target sum each candidate sends left, per node.rightMass(m, x)— Let \(M \in \mathbb{R}^{n \times n}\) and \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rightMass}(M, X)_{jtf} = N_j - N^L_{jtf}\). The mass each candidate sends right: the node’s mass less the left mass.rightSum(m, x, y)— Let \(M \in \mathbb{R}^{n \times n}\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\). Then \(\operatorname{rightSum}(M, X, y)_{jtf} = S_j - S^L_{jtf}\). The target sum each candidate sends right.sideScore(s, n)— Let \(s\) and \(\nu\) (the argumentssandn) be tensors of the same shape. Then \(\operatorname{sideScore}(s, \nu) = \frac{s \odot s}{\nu + [\nu = 0]}\), entrywise. The explained sum of squares of one side, s squared over n, with a zero mass treated as one.viable(m, x)— Let \(M \in \mathbb{R}^{n \times n}\) and \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{viable}(M, X)_{jtf} = [\,N^L_{jtf} > 0\,]\,[\,N^R_{jtf} > 0\,]\). One where both children of a candidate would be non-empty.childScore(m, x, y)— Let \(M, X, y\) be as inleftSum. Then \(\operatorname{childScore}(M, X, y) = \operatorname{sideScore}(S^L, N^L) + \operatorname{sideScore}(S^R, N^R)\). The explained sum of squares of both children together.
Choosing one candidate
The best candidate of each node is found with two reductions and a comparison. Ties are broken by position, not averaged: averaging two different splits would send a fraction of a row down each branch, which a tree cannot hold. Each candidate is numbered \(t \cdot p + f\) and the lowest-numbered of the best is chosen.
In the items below, \(G = \operatorname{splitGain}(M, X, y)\) and \(T = \operatorname{topMask}(M, X, y)\).
Building blocks
bestGain(g)— Let \(G \in \mathbb{R}^{n \times n \times p}\) (the argumentg). Then \(\operatorname{bestGain}(G)_j = \max_{t} \max_{f} G_{jtf}\). The largest gain at each node slot.topMask(m, x, y)— Let \(M \in \mathbb{R}^{n \times n}\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\). Then \(\operatorname{topMask}(M, X, y)_{jtf} = [\,G_{jtf} = \operatorname{bestGain}(G)_j\,]\, \operatorname{viable}(M, X)_{jtf}\). One at every viable candidate that attains its node’s best gain.candidateIds(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{candidateIds}(X)_{tf} = t\,p + f\) for \(t = 0, \dots, n-1\) and \(f = 0, \dots, p-1\). The number t times p plus f of every candidate.firstOfTop(m, x, y)— Let \(M, X, y\) be as intopMaskand \(c = \operatorname{candidateIds}(X)\). Then \(\operatorname{firstOfTop}(M, X, y)_j = \min_{t} \min_{f} \bigl( T_{jtf}\, c_{tf} + (1 - T_{jtf})\, 10^{300} \bigr)\), which is \(10^{300}\) for a node with no top candidate. The lowest candidate number among each node’s best.chosen(m, x, y)— Let \(M, X, y\) and \(c\) be as infirstOfTop. Then \(\operatorname{chosen}(M, X, y)_{jtf} = T_{jtf}\, [\,c_{tf} = \operatorname{firstOfTop}(M, X, y)_j\,]\). A one-hot n×n×p mask of the single chosen candidate per node (all zero for a node with no viable candidate).splitsHere(m, x, y)— Let \(M, X, y\) be as intopMask. Then \(\operatorname{splitsHere}(M, X, y)_j = [\,\operatorname{bestGain}(G)_j > 10^{-12}\,]\). One for each node whose best gain exceeds 1e-12.
The split a node keeps
A node’s split is two tensors: a one-hot feature row and a threshold. A node that does not split has an all-zero feature row and a threshold of zero.
In the items below, \(C = \operatorname{chosen}(M, X, y)\) and \(h = \operatorname{splitsHere}(M, X, y)\).
Building blocks
featureAt(m, x, y)— Let \(M \in \mathbb{R}^{n \times n}\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\). Then \(\operatorname{featureAt}(M, X, y)_{jf} = h_j \sum_{t} C_{jtf}\). An n×p matrix, one-hot per node slot on the chosen feature.thresholdAt(m, x, y)— Let \(M, X, y\) be as infeatureAt. Then \(\operatorname{thresholdAt}(M, X, y)_j = h_j \sum_{t} \sum_{f} C_{jtf} X_{tf}\). A vector of the chosen threshold per node slot.
Moving a level down
Moving a row to the next level is a contraction against a link matrix built by comparing slot numbers: slot j links to slot 2j on the left and 2j+1 on the right. The membership matrix stays n×n on every level, so the loop’s shapes never change.
membership
membership(d, x, y) · d depth scalar · x data n×p · y targets vector n → membership n×n
Let \(d \in \mathbb{N}\) (the argument d) be the depth, \(X \in \mathbb{R}^{n \times p}\) (the argument x) the data and \(y \in \mathbb{R}^{n}\) (the argument y) the targets. Let \(M^{(0)} = \mathbf{1}_n e_0^{\top}\) and \(M^{(k+1)} = \operatorname{dropTo}(M^{(k)}, X, y)\) for \(k = 0, 1, \dots, d-1\). Then
After d levels of growth, row i has a one in the column of the level-d node slot that holds it. Slot numbers are local to the level.
Intuition A tree is a way of sorting rows into groups, and this matrix is that sorting. Nothing else is stored. Each level sends every row of a node to child \(2j\) or \(2j+1\). A node that cannot improve sends all of its rows left, so a finished leaf keeps its rows as the levels go on.
Use it when you want the segments a tree finds: rows that share a column share a leaf. That is useful for grouping customers by a target, one-hot leaf features for another model, or per-leaf statistics. A slot number means nothing without its depth. Keep \(2^d \le n\), since that is the usable depth.
data-segmentationinterpretable-modelsfeature-engineering
def membership(d, x, y) ≔ growLoop(d, x, y, atRoot(x))
Building blocks
slotIds(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{slotIds}(X) = (0, 1, \dots, n-1) \in \mathbb{R}^{n}\). The slot numbers.heldSlots(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{heldSlots}(X)_{jc} = j\). An n×n matrix holding the parent slot number in every row.laidSlots(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{laidSlots}(X)_{jc} = c\). An n×n matrix holding the child slot number in every column.leftLink(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{leftLink}(X)_{jc} = [\,c = 2j\,]\). One where child slot c equals 2j for parent slot j.rightLink(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rightLink}(X)_{jc} = [\,c = 2j + 1\,]\). One where child slot c equals 2j + 1.wentLeft(m, x, y)— Let \(M \in \mathbb{R}^{n \times n}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(C = \operatorname{chosen}(M, X, y)\). Then \(\operatorname{wentLeft}(M, X, y)_{ij} = \sum_{t} \sum_{f} [\,X_{if} \le X_{tf}\,]\, C_{jtf}\), which is \([\,X_{i f^*} \le X_{t^* f^*}\,]\) when node \(j\) chose candidate \((t^*, f^*)\). One where row i goes left under node j’s chosen split.leftward(m, x, y)— Let \(M, X, y\) be as inwentLeftand \(h = \operatorname{splitsHere}(M, X, y)\). Then \(\operatorname{leftward}(M, X, y)_{ij} = \operatorname{wentLeft}(M, X, y)_{ij}\, h_j + (1 - h_j)\).wentLeft, with every row of a node that does not split sent left.dropTo(m, x, y)— Let \(M, X, y\) be as inwentLeftand \(\Lambda = \operatorname{leftward}(M, X, y)\). Then \(\operatorname{dropTo}(M, X, y)_{ic} = \sum_{j} M_{ij} \bigl( \Lambda_{ij} [\,c = 2j\,] + (1 - \Lambda_{ij}) [\,c = 2j + 1\,] \bigr)\). The membership one level down.rootOf(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rootOf}(X)_j = [\,j = 0\,]\), that is \(\operatorname{rootOf}(X) = e_0 \in \mathbb{R}^{n}\). The one-hot vector of slot 0.atRoot(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{atRoot}(X) = \mathbf{1}_n e_0^{\top}\), that is \(\operatorname{atRoot}(X)_{ij} = [\,j = 0\,]\). The starting membership: every row in slot 0.growLoop(d, x, y, m)— Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(M \in \mathbb{R}^{n \times n}\). Then \(\operatorname{growLoop}(0, X, y, M) = M\) and \(\operatorname{growLoop}(d, X, y, M) = \operatorname{growLoop}\bigl(d - 1, X, y, \operatorname{dropTo}(M, X, y)\bigr)\). AppliesdropTod times.
What the tree answers
Each leaf predicts the mean target of the rows it holds.
leafValue
leafValue(d, x, y) · d depth scalar · x data n×p · y targets vector n → leaf means vector n
Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\) (the arguments d, x, y), and let \(M = \operatorname{membership}(d, X, y)\). Then
The mean of y in each level-d node slot; empty slots read zero.
Intuition Under squared error the best constant for a group is its mean, so that is what every leaf predicts. For 0/1 labels the mean is the fraction of positives in the leaf, which is the tree’s class probability there.
Use it when you want to inspect or report what each segment predicts, for example the average spend or default rate per leaf, together with membership for who is in it. An empty slot and a leaf whose mean is 0 both read 0. Check the column totals of membership to tell them apart. For per-row predictions use treeFit or treePredict.
regressioninterpretable-modelsdata-segmentation
def leafValue(d, x, y) ≔ einsum("ij,i->j", membership(d, x, y), y) / safe(einsum("ij->j", membership(d, x, y)))
treeFit
treeFit(d, x, y) · d depth scalar · x data n×p · y targets vector n → fitted values vector n
Let \(d, X, y\) and \(M = \operatorname{membership}(d, X, y)\) be as in leafValue. Then
The prediction of a depth-d regression tree at every training row.
Intuition A tree’s prediction is a step function: the feature space is cut into axis-aligned boxes, and each box answers its own mean. It captures thresholds and interactions without any feature scaling, but it is flat inside every box.
Use it when you want in-sample fitted values, residuals to inspect, or a residual target for the next boosting round (as boosted uses it). These are training predictions, so they look better than the model is. A tree deep enough for \(2^d\) to approach \(n\) can fit every row. Choose d by the error of treePredict on held-out rows.
regressionnonlinear-modellingresidual-diagnostics
def treeFit(d, x, y) ≔ einsum("ij,j->i", membership(d, x, y), leafValue(d, x, y))
Predicting a new row
A prediction re-walks the growth. One loop carries the training membership and the query membership together, so at every level each query is tested against the split the training rows have just made. A node that never split sends a query left, as it does the training rows.
treePredict
treePredict(d, x, y, z) · d depth scalar · x data n×p · y targets vector n · z queries m×p → predictions vector m
Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(Z \in \mathbb{R}^{m \times p}\) (the arguments d, x, y, z). Let \(M^{(k)}\) be the training memberships of membership, let \(V^{(0)} = \mathbf{1}_m e_0^{\top}\) and \(V^{(k+1)} = \operatorname{sendQuery}(Z, X, y, M^{(k)}, V^{(k)})\) for \(k = 0, 1, \dots, d-1\). Then
The depth-d regression tree grown on x and y, evaluated at each query row.
Intuition A new row goes down the same splits the training rows made and takes the mean of the leaf it lands in. The tree is grown again during the same walk, so the prediction cannot disagree with the fit.
Use it when you need a regression tree’s answer on new or held-out rows, including validation to choose d. Every call grows the tree again at about \(n^3 p\) per level, so put all queries in one z, not one call per row. Predictions never leave the range of the training leaf means: a query far beyond the data gets the edge leaf’s value, so trees do not extrapolate trends.
regressionnonlinear-modellingtabular-prediction
def treePredict(d, x, y, z) ≔ einsum("ij,j->i", queryAt(d, x, y, z), leafValue(d, x, y))
treeClassify
treeClassify(d, x, y, z) · d depth scalar · x data n×p · y 0/1 labels vector n · z queries m×p → predicted labels, 0 or 1 vector m
Let \(d, X, Z\) be as in treePredict and let \(y \in \{0, 1\}^{n}\) (the argument y) be the labels. Then
A binary classification tree: for 0/1 labels a leaf’s mean is the class probability, so the prediction is one where it is at least 0.5.
Intuition On 0/1 labels, cutting squared error is the same as cutting Gini impurity (up to a factor), so the regression tree already is a classification tree. Each leaf holds the fraction of positives, and the label is a vote at one half.
Use it when you need a readable yes/no classifier, such as a churn, fraud or approval rule someone can follow by hand. When errors cost unequal amounts or the classes are imbalanced, take treePredict (the leaf probability) and apply your own cutoff. Labels must be 0 or 1, since multi-class trees are not provided. For more stable decisions, use forestClassify.
classificationinterpretable-models
def treeClassify(d, x, y, z) ≔ treePredict(d, x, y, z) ≥ 0.5
Building blocks
readAt(z, f)— Let \(Z \in \mathbb{R}^{m \times p}\) and \(F \in \mathbb{R}^{n \times p}\) (the argumentf). Then \(\operatorname{readAt}(Z, F) = Z F^{\top}\), that is \(\operatorname{readAt}(Z, F)_{ij} = \sum_{f} Z_{if} F_{jf}\). The value of each query on each node’s chosen feature, m×n.queryLeft(z, x, y, m)— Let \(Z \in \mathbb{R}^{m \times p}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(M \in \mathbb{R}^{n \times n}\). Then \(\operatorname{queryLeft}(Z, X, y, M)_{ij} = [\,\operatorname{readAt}(Z, \operatorname{featureAt}(M, X, y))_{ij} \le \operatorname{thresholdAt}(M, X, y)_j\,]\). One where query i goes left at node j.sendQuery(z, x, y, m, w)— Let \(Z, X, y, M\) be as inqueryLeft, let \(V \in \mathbb{R}^{m \times n}\) (the argumentw) be the query membership and \(\Gamma = \operatorname{queryLeft}(Z, X, y, M)\). Then \(\operatorname{sendQuery}(Z, X, y, M, V)_{ic} = \sum_{j} V_{ij} \bigl( \Gamma_{ij} [\,c = 2j\,] + (1 - \Gamma_{ij}) [\,c = 2j + 1\,] \bigr)\). The query membership one level down.bothLoop(d, x, y, z, m, w)— Let \(d \in \mathbb{N}\), \(X, y, Z\) be as inqueryLeft, \(M \in \mathbb{R}^{n \times n}\) and \(V \in \mathbb{R}^{m \times n}\) (the argumentw). Then \(\operatorname{bothLoop}(0, X, y, Z, M, V) = V\) and \(\operatorname{bothLoop}(d, X, y, Z, M, V) = \operatorname{bothLoop}\bigl(d - 1, X, y, Z, \operatorname{dropTo}(M, X, y), \operatorname{sendQuery}(Z, X, y, M, V)\bigr)\). Moves the training and query memberships down together d times and answers the query membership.queryAt(d, x, y, z)— Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(Z \in \mathbb{R}^{m \times p}\). Then \(\operatorname{queryAt}(d, X, y, Z) = \operatorname{bothLoop}(d, X, y, Z, \mathbf{1}_n e_0^{\top}, \mathbf{1}_m e_0^{\top})\). The level-d query membership, m×n.
The root’s split
Because slot numbers are local to their level, there is no whole-tree table of splits to read. The root is the exception — level zero holds one node — and it is the split most often inspected.
rootFeature
rootFeature(x, y) · x data n×p · y targets vector n → one-hot feature n×p
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) and \(y \in \mathbb{R}^{n}\) (the argument y). Then
Row 0 is one-hot on the feature the root splits on (all zeros if the root does not split); the other rows are zero.
Intuition The root is the single best question you could ask about the data, so the feature it uses is the one that explains the most target variance on its own. Slot numbers below the root depend on their level, so the root is the one split that can be read unambiguously.
Use it when you want a quick screen of which feature matters most, a one-rule baseline, or a leakage check. A root split on an ID or timestamp column is a warning sign. Only row 0 carries information. Pair it with rootCut for the threshold, and use splitGain to see how close the runner-up features came.
feature-importanceinterpretable-modelsdata-leakage
def rootFeature(x, y) ≔ featureAt(atRoot(x), x, y)
rootCut
rootCut(x, y) · x data n×p · y targets vector n → threshold vector n
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) and \(y \in \mathbb{R}^{n}\) (the argument y). Then
Entry 0 is the root’s threshold: rows whose chosen feature is at most this value go left. The other entries are zero.
Intuition The threshold is an actual value from the data, the largest value on the left side, not a midpoint between two values. So a new value that falls between it and the next observed value goes right.
Use it when you want the cut point of the best single rule, such as “age ≤ 42” or “dose ≤ 3”. It is useful for turning a continuous measure into a two-level flag, or for reporting the root rule next to rootFeature. When the root does not split, entry 0 is 0, so check rootFeature before trusting a zero cut.
threshold-selectioninterpretable-modelsdiscretisation
def rootCut(x, y) ≔ thresholdAt(atRoot(x), x, y)
A tree grown on weighted rows
Every mass, sum and gain is a contraction against the membership matrix, so placing row i in the root with weight \(w_i\) instead of one weights all of them by exactly \(w_i\). A bootstrap sample is just a vector of counts, so no second implementation is needed for bagging.
treeAt
treeAt(d, x, y, w, z) · d depth scalar · x data n×p · y targets vector n · w row weights or bootstrap counts vector n · z queries m×p → predictions vector m
Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(Z \in \mathbb{R}^{m \times p}\) (the arguments d, x, y, z), and let \(w \in \mathbb{R}^{n}\) (the argument w) be non-negative row weights. Let \(M^{(0)} = w\, e_0^{\top}\), \(V^{(0)} = \mathbf{1}_m e_0^{\top}\) and, for \(k = 0, 1, \dots, d-1\), \(M^{(k+1)} = \operatorname{dropTo}(M^{(k)}, X, y)\) and \(V^{(k+1)} = \operatorname{sendQuery}(Z, X, y, M^{(k)}, V^{(k)})\). Then
A depth-d tree grown on rows weighted by w, evaluated at the queries. A row with weight zero takes no part.
Intuition Weighting a row by \(w_i\) acts as if it were copied \(w_i\) times: every node mass, target sum and gain counts it that much. A bootstrap sample is just such a vector of copy counts, and a 0/1 vector is a subset.
Use it when some rows should count more, for example to upweight a rare class, correct a biased sample, or weight by exposure. Use it for one bagged tree (inside forest), or to train on a cross-validation fold with a 0/1 mask instead of copying data. Weights must be non-negative. With all weights 1 it is treePredict.
sample-weightingclass-imbalancecross-validation
def treeAt(d, x, y, w, z) ≔ einsum("ij,j->i", queryFrom(d, x, y, z, w), leafOf(d, x, y, w))
Building blocks
seededAt(x, w)— Let \(X \in \mathbb{R}^{n \times p}\) and \(w \in \mathbb{R}^{n}\). Then \(\operatorname{seededAt}(X, w) = w\, e_0^{\top}\), that is \(\operatorname{seededAt}(X, w)_{ij} = w_i [\,j = 0\,]\). The starting membership with row i in slot 0 at weight w.grownFrom(d, x, y, w)— Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(w \in \mathbb{R}^{n}\). Then \(\operatorname{grownFrom}(d, X, y, w) = \operatorname{growLoop}(d, X, y, w\, e_0^{\top})\). The weighted membership after d levels.leafOf(d, x, y, w)— Let \(d, X, y, w\) be as ingrownFromand \(G = \operatorname{grownFrom}(d, X, y, w)\). Then \(\operatorname{leafOf}(d, X, y, w)_j = \frac{\sum_{i} G_{ij} y_i}{\operatorname{safe}(\sum_{i} G_{ij})}\). The weighted mean target of each leaf slot.queryFrom(d, x, y, z, w)— Let \(d, X, y, w\) be as ingrownFromand \(Z \in \mathbb{R}^{m \times p}\). Then \(\operatorname{queryFrom}(d, X, y, Z, w) = \operatorname{bothLoop}(d, X, y, Z, w\, e_0^{\top}, \mathbf{1}_m e_0^{\top})\). The query membership after d levels of the weighted tree.
The random forest
A random forest here is bagging: w is the B×n matrix of bootstrap counts, one row per tree, so the forest has as many trees as rows the caller drew and is a deterministic function of those draws. Every tree considers every feature; there is no per-split feature subsampling. The trees are visited by a loop rather than lifted by one more index, which holds one tree’s worth of memory at a time.
forest
forest(d, x, y, z, w) · d depth scalar · x data n×p · y targets vector n · z queries m×p · w bootstrap counts B×n → predictions vector m
Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and \(Z \in \mathbb{R}^{m \times p}\) (the arguments d, x, y, z), and let \(W \in \mathbb{R}^{B \times n}\) (the argument w) hold bootstrap counts, with row \(w_b\) for tree \(b = 0, 1, \dots, B-1\). Then
The average of the B weighted trees’ predictions at the queries.
Intuition One deep tree is jumpy: small changes in the data move its splits. Trees grown on different bootstrap resamples make different mistakes, and averaging them cancels much of that noise while keeping each tree’s flexibility. Every tree here sees every feature, so this is bagging, and the trees are more alike than in a random forest that subsamples features.
Use it when a single tree overfits and you want a smoother, more accurate regressor without tuning much. Deeper trees gain most from averaging. The cost is B trees, but memory stays at one tree. Draw w as bootstrap counts (one row per tree) from a seeded stream, and the forest is reproducible.
variance-reductionoverfitting-controlregression
def forest(d, x, y, z, w) ≔ forestLoop(rows(w), d, x, y, z, w, treeIds(w), zeroed(rowOnes(z))) / rows(w)
forestClassify
forestClassify(d, x, y, z, w) · arguments as forest, with 0/1 labels y → predicted labels vector m
Let \(d, X, Z, W\) be as in forest and let \(y \in \{0, 1\}^{n}\) (the argument y) be the labels. Then
One where the forest’s average is at least 0.5.
Intuition Each tree gives the fraction of positives in the leaf a query lands in, and the forest averages those fractions before voting at one half. This is a soft vote: a tree that is 90% sure counts more than one that is 55% sure, unlike a majority count of hard labels.
Use it when you need a binary classifier that is more stable than treeClassify and less sensitive to the exact training sample. For a different cutoff, or a probability score to rank by, use forest on the 0/1 labels directly. Labels must be 0 or 1.
classificationvariance-reductionoverfitting-control
def forestClassify(d, x, y, z, w) ≔ forest(d, x, y, z, w) ≥ 0.5
Building blocks
treeIds(w)— Let \(W \in \mathbb{R}^{B \times n}\). Then \(\operatorname{treeIds}(W) = (0, 1, \dots, B-1) \in \mathbb{R}^{B}\). The tree numbers.drawRow(w, ix, k)— Let \(W \in \mathbb{R}^{B \times n}\), let \(\iota \in \mathbb{R}^{B}\) (the argumentix) be tree numbers and \(k \in \{0, \dots, B-1\}\). Then \(\operatorname{drawRow}(W, \iota, k)_i = \sum_{b} W_{bi} [\,\iota_b = k\,]\), which is \(W_{ki}\) when \(\iota = (0, 1, \dots, B-1)\). Row k of the bootstrap matrix: the counts for tree k.forestLoop(b, d, x, y, z, w, ix, a)— Let \(b \in \mathbb{N}\), \(d, X, y, Z, W\) be as inforest, \(\iota \in \mathbb{R}^{B}\) and \(a \in \mathbb{R}^{m}\). Then \(\operatorname{forestLoop}(0, d, X, y, Z, W, \iota, a) = a\) and \(\operatorname{forestLoop}(b, d, X, y, Z, W, \iota, a) = \operatorname{forestLoop}\bigl(b - 1, d, X, y, Z, W, \iota, a + \operatorname{treeAt}(d, X, y, \operatorname{drawRow}(W, \iota, b - 1), Z)\bigr)\). Adds the predictions of the remaining b trees to the running total a.
Gradient boosting
Boosting fits trees in sequence, each to what the ones before it got wrong. Under squared loss the negative gradient is the residual, so each round fits a tree to \(y - p\) and adds a fraction r (the shrinkage) of it:
\[p_0 = \bar y, \qquad p_{k+1} = p_k + r\, T_k(y - p_k)\]Predicting a new row re-runs all the rounds, carrying the training and query sums together so they stay in step.
boosted
boosted(t, d, r, x, y) · t number of rounds scalar · d depth of each tree scalar · r shrinkage scalar · x data n×p · y targets vector n → fitted values vector n
Let \(t \in \mathbb{N}\) (the argument t) be the number of rounds, \(d \in \mathbb{N}\) (the argument d) the depth of each tree, \(r > 0\) (the argument r) the shrinkage, \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\) (the arguments x, y). Let \(s_0 = \bar y\, \mathbf{1}_n\) with \(\bar y = \frac{1}{n}\sum_{i} y_i\), and \(s_{k+1} = s_k + r\, \operatorname{treeFit}(d, X, y - s_k)\) for \(k = 0, 1, \dots, t-1\). Then
The boosted model’s prediction at every training row after t rounds, starting from the mean of y.
Intuition Start from the mean, then repeatedly fit a small tree to whatever is still unexplained and add a fraction r of it. Each weak tree fixes part of the remaining error. The shrinkage stops any one tree from being trusted fully, which is what makes many small steps generalise better than a few large ones.
Use it when you want accuracy on tabular regression beyond a single tree or a bagged forest. Common starting points are shallow trees (d of 1 to 3) and small r with more rounds. These are fitted values on the training rows only, and there is no early stopping. Choose t by held-out error from boostPredict.
regressionbias-correctiontabular-prediction
def boosted(t, d, r, x, y) ≔ boostRun(t, d, r, x, y, fromMean(y))
boostPredict
boostPredict(t, d, r, x, y, z) · arguments as boosted · z queries m×p → predictions vector m
Let \(t, d, r, X, y\) and the training sums \(s_k\) be as in boosted, and let \(Z \in \mathbb{R}^{m \times p}\) (the argument z) be the queries. Let \(q_0 = \bar y\, \mathbf{1}_m\) and \(q_{k+1} = q_k + r\, \operatorname{treePredict}(d, X, y - s_k, Z)\) for \(k = 0, 1, \dots, t-1\). Then
The same boosted model evaluated at the query rows.
Intuition Round k’s tree is fitted to residuals that depend on every earlier round, so a stored model does not exist to query. Instead the rounds are replayed, with the training sum and the query sum updated together so each query gets exactly the trees the fit built.
Use it when you need the boosted model’s answers for new or held-out rows, including the validation curve used to choose t, d and r. Every round grows its tree both for the training sum and for the queries, so it costs more than boosted. Put all queries in one call.
regressiontabular-predictionmodel-selection
def boostPredict(t, d, r, x, y, z) ≔ boostBoth(t, d, r, x, y, z, fromMean(y), zeroed(rowOnes(z)) + mean(y))
Building blocks
boostBoth(t, d, r, x, y, z, p, q)— Let \(t, d \in \mathbb{N}\), \(r > 0\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\), the training sum \(s \in \mathbb{R}^{n}\) (the argumentp), \(Z \in \mathbb{R}^{m \times p}\) and \(q \in \mathbb{R}^{m}\). Then \(\operatorname{boostBoth}(0, d, r, X, y, Z, s, q) = q\) and \(\operatorname{boostBoth}(t, d, r, X, y, Z, s, q) = \operatorname{boostBoth}\bigl(t - 1, d, r, X, y, Z, s + r \operatorname{treeFit}(d, X, y - s), q + r \operatorname{treePredict}(d, X, y - s, Z)\bigr)\). Runs t rounds updating the training sum p and query sum q together; answers q.boostRun(t, d, r, x, y, p)— Let \(t, d \in \mathbb{N}\), \(r > 0\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and the training sum \(s \in \mathbb{R}^{n}\) (the argumentp). Then \(\operatorname{boostRun}(0, d, r, X, y, s) = s\) and \(\operatorname{boostRun}(t, d, r, X, y, s) = \operatorname{boostRun}\bigl(t - 1, d, r, X, y, s + r \operatorname{treeFit}(d, X, y - s)\bigr)\). Runs t rounds on the training sum p and answers it.fromMean(y)— Let \(y \in \mathbb{R}^{n}\). Then \(\operatorname{fromMean}(y) = \bar y\, \mathbf{1}_n\) with \(\bar y = \frac{1}{n}\sum_{i} y_i\). A vector the shape of y filled with its mean.
Penalised leaves
The objective XGBoost adds to gradient boosting is second order with an L2 penalty \(\lambda\) on the leaf weights. Under squared loss the Hessian is one per row, so the optimal leaf weight becomes \(G/(H + \lambda)\) — the residual sum over the leaf’s count plus \(\lambda\) — which shrinks small leaves towards zero. In this module the penalty is applied to the leaf values only: the splits themselves are still chosen by the unpenalised criterion. The engineering that makes XGBoost fast (histogram binning, sparsity-aware split finding, out-of-core training) is not reproduced.
treePenalised
treePenalised(d, x, y, l) · d depth scalar · x data n×p · y targets vector n · l leaf penalty λ scalar → fitted values vector n
Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\) (the arguments d, x, y), let \(\lambda > 0\) (the argument l) be the leaf penalty, and let \(M = \operatorname{membership}(d, X, y)\). Then
A depth-d tree whose leaves predict \(\sum y / (\text{count} + \lambda)\), evaluated at the training rows. Unlike leafValue, the divisor has no zero guard, so \(\lambda > 0\) is assumed: with \(\lambda = 0\) every empty slot computes \(0/0\).
Intuition Adding \(\lambda\) to the count is like adding \(\lambda\) imaginary rows with target 0 to every leaf. A leaf of \(c\) rows keeps a fraction \(c/(c+\lambda)\) of its mean, so big leaves barely move and tiny leaves, whose means are least trustworthy, are pulled strongly toward zero.
Use it when the targets are residuals centred on zero, as inside gradientBoost, and you want small leaves to count less. On raw targets it biases every prediction toward 0, so use treeFit there. Only the leaf values are penalised, and the splits are the same as the unpenalised tree’s.
regularisationoverfitting-controlsmall-samples
def treePenalised(d, x, y, l) ≔ einsum("ij,j->i", membership(d, x, y), leafPenalised(d, x, y, l))
gradientBoost
gradientBoost(t, d, r, l, x, y) · t number of rounds scalar · d depth scalar · r shrinkage scalar · l leaf penalty λ scalar · x data n×p · y targets vector n → fitted values vector n
Let \(t \in \mathbb{N}\) (the argument t), \(d \in \mathbb{N}\) (the argument d), \(r > 0\) (the argument r), \(\lambda > 0\) (the argument l), \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^{n}\) (the arguments x, y). Let \(s_0 = \bar y\, \mathbf{1}_n\) and \(s_{k+1} = s_k + r\, \operatorname{treePenalised}(d, X, y - s_k, \lambda)\) for \(k = 0, 1, \dots, t-1\). Then
Gradient boosting with penalised leaves, starting from the mean of y, answered at the training rows.
Intuition The same residual-fitting loop as boosted, with a second brake added. Shrinkage r scales every tree equally, while \(\lambda\) scales each leaf by how many rows back it, so corrections based on a few rows are damped the most. This is XGBoost’s leaf rule under squared loss.
Use it when boosting overfits small groups or outliers, and a lower r alone would need too many rounds. A leaf of \(c\) rows keeps half its mean when \(\lambda = c\), which helps you pick the scale. It answers the training rows only, and the module has no query version, so for predictions on new rows use boostPredict (unpenalised).
regressionregularisationoverfitting-control
def gradientBoost(t, d, r, l, x, y) ≔ gradientRun(t, d, r, l, x, y, fromMean(y))
Building blocks
leafPenalised(d, x, y, l)— Let \(d \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\), \(\lambda > 0\) (the argumentl) and \(M = \operatorname{membership}(d, X, y)\). Then \(\operatorname{leafPenalised}(d, X, y, \lambda)_j = \frac{\sum_{i} M_{ij} y_i}{\sum_{i} M_{ij} + \lambda}\). Each leaf’s target sum divided by its count plus l.gradientRun(t, d, r, l, x, y, p)— Let \(t, d \in \mathbb{N}\), \(r > 0\), \(\lambda > 0\), \(X \in \mathbb{R}^{n \times p}\), \(y \in \mathbb{R}^{n}\) and the training sum \(s \in \mathbb{R}^{n}\) (the argumentp). Then \(\operatorname{gradientRun}(0, d, r, \lambda, X, y, s) = s\) and \(\operatorname{gradientRun}(t, d, r, \lambda, X, y, s) = \operatorname{gradientRun}\bigl(t - 1, d, r, \lambda, X, y, s + r \operatorname{treePenalised}(d, X, y - s, \lambda)\bigr)\). Runs t rounds of penalised boosting on the training sum p.
Examples
The examples below share one file, Grove.hb:
open Prelude.Trees
fn fittedTree(d, x, y) ↦ treeFit(d, x, y)
fn predictedTree(d, x, y, z) ↦ treePredict(d, x, y, z)
fn classifiedTree(d, x, y, z) ↦ treeClassify(d, x, y, z)
fn firstFeature(x, y) ↦ rootFeature(x, y)
fn firstCut(x, y) ↦ rootCut(x, y)
fn forestOf(d, x, y, z, w) ↦ forest(d, x, y, z, w)
fn weightedTree(d, x, y, w, z) ↦ treeAt(d, x, y, w, z)
fn boostedFit(t, d, r, x, y) ↦ boosted(t, d, r, x, y)
fn boostedQuery(t, d, r, x, y, z) ↦ boostPredict(t, d, r, x, y, z)
fn penalisedBoost(t, d, r, l, x, y) ↦ gradientBoost(t, d, r, l, x, y)
Fit a step and predict new points
Six points with one feature and a step in the targets. A single split (depth 1) finds the step at 3 exactly; the queries 2.5, 3.5 and 10 fall on either side.
using MLambda.Hilbert.Runtime;
using Examples;
var x = (Tensor)new double[,] { { 1 }, { 2 }, { 3 }, { 4 }, { 5 }, { 6 } };
var y = (Tensor)new[] { 1d, 1, 1, 5, 5, 5 };
var fit = Grove.FittedTree(1.0, x, y);
var pred = Grove.PredictedTree(1.0, x, y, (Tensor)new double[,] { { 2.5 }, { 3.5 }, { 10 } });
Output (verified):
fit = [1, 1, 1, 5, 5, 5]
pred = [1, 5, 5]
Grow two levels, and see a tie broken by position
With targets [1, 2, 3, 10, 11, 15] the root splits at 3. The left node can split off {1} or {3} equally well; the first candidate wins, giving leaves 1 and 2.5. The right node separates 15 from 10 and 11.
var y2 = (Tensor)new[] { 1d, 2, 3, 10, 11, 15 };
var fit2 = Grove.FittedTree(2.0, x, y2);
Output (verified):
fit2 = [1, 2.5, 2.5, 10.5, 10.5, 15]
Inspect the root split and classify
Two features; the labels are determined by the second one (at most 3 means class 0). The root’s feature row is [0, 1] and its cut is 3, both in slot 0. The query [1, 2.5] is class 0 and [9, 5] is class 1, whatever the first feature says.
var xs = (Tensor)new double[,] { { 5, 1 }, { 1, 2 }, { 6, 3 }, { 2, 4 }, { 4, 5 }, { 3, 6 } };
var labels = (Tensor)new[] { 0d, 0, 0, 1, 1, 1 };
var feature = Grove.FirstFeature(xs, labels);
var cut = Grove.FirstCut(xs, labels);
var classes = Grove.ClassifiedTree(1.0, xs, labels, (Tensor)new double[,] { { 1, 2.5 }, { 9, 5 } });
Output (verified; feature is 6×2, shown row by row):
feature = [[0, 1], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]]
cut = [3, 0, 0, 0, 0, 0]
classes = [0, 1]
Average two bootstrap trees
Two trees: the first sees every row once, the second the counts [2, 0, 1, 0, 1, 2]. The first tree predicts 2 and 12 for the queries 2 and 5. The weighted tree’s leaves are \((2 \cdot 1 + 3)/3 = 1.6667\) and \((11 + 2 \cdot 15)/3 = 13.6667\). The forest averages the two.
var counts = (Tensor)new double[,] { { 1, 1, 1, 1, 1, 1 }, { 2, 0, 1, 0, 1, 2 } };
var queries = (Tensor)new double[,] { { 2 }, { 5 } };
var second = Grove.WeightedTree(1.0, x, y2, (Tensor)new[] { 2d, 0, 1, 0, 1, 2 }, queries);
var bagged = Grove.ForestOf(1.0, x, y2, queries, counts);
Output (verified, six decimals):
weighted tree = [1.666667, 13.666667]
forest = [1.833333, 12.833333]
Boost stumps and watch the residual halve
Boosting starts at the mean, 3. Each depth-1 tree fits the residual exactly and shrinkage 0.5 closes half of it, so after three rounds the remaining error is \(2 \cdot 0.5^3 = 0.25\).
var start = Grove.BoostedFit(0.0, 1.0, 0.5, x, y);
var rounds3 = Grove.BoostedFit(3.0, 1.0, 0.5, x, y);
var ahead = Grove.BoostedQuery(3.0, 1.0, 0.5, x, y, (Tensor)new double[,] { { 0 }, { 9 } });
Output (verified):
t = 0 = [3, 3, 3, 3, 3, 3]
t = 3 = [1.25, 1.25, 1.25, 4.75, 4.75, 4.75]
queries = [1.25, 4.75]
Penalise the leaves
One round with shrinkage 1 and \(\lambda = 3\). Each leaf holds three rows with residual −2 or +2, so its value is \(\pm 6/(3 + 3) = \pm 1\) instead of \(\pm 2\), and the fit moves only halfway from the mean.
var penalised = Grove.PenalisedBoost(1.0, 1.0, 1.0, 3.0, x, y);
Output (verified):
penalised = [2, 2, 2, 4, 4, 4]