Classical models

In short: every model has the same shape. Fit takes an f32 [n, d] tensor (and f32 targets or i32 labels), copies it into a double linalg::Matrix through one shared helper, and returns a Status; Predict (and PredictProba where it makes sense) converts the answer back to f32 or i32 tensors. Options are plain structs with scikit-learn's defaults where there is one. The solvers are the textbook ones, written out: Cholesky normal equations or gradient descent, Newton with a line search, Lloyd's algorithm after k-means++ seeding, an SVD, counting, a CART sweep, bootstrap forests, LightGBM-style histogram boosting and LIBSVM-style SMO. Randomness comes from seeded Philox generators and parallel loops from the thread pool, so results do not depend on the thread count. Nothing here trains through autograd, and the inference stack does not use it.

Each header links its concept page on /ml-models: linear regression and regularization, logistic regression, nearest neighbors, clustering, PCA, Bayes' rule, decision trees, random forests, gradient boosting and support vector machines.

The shared shape

The helpers in internal are the only shared code (common.h L12–31). Their checks are RML_CHECKs, so a wrong dtype or shape is a programmer error that aborts; data problems a caller cannot rule out in advance (singular systems, too few classes, k > n, non-convergence) come back as a Status.

The models at a glance

Linear models

LinearRegression

The objective is ‖y − Xw − b‖² + l2‖w‖², ordinary least squares at l2 = 0 (linear_regression.h L8–19). Fit centers X and y when fit_intercept is on, which removes the intercept from the problem and keeps it unpenalized, solves for w, then recovers b = mean(y) − mean(x)·w (linear_regression.cc L67–102).

• Normal equations. Forms XcᵀXc + l2·I and factors it with Cholesky; a failed factorization is the singular case and returns InvalidArgument with "set l2 > 0 for ridge regression" (L19–31).
• Gradient descent. Full-batch steps on the objective divided by n, stopping once the gradient's largest entry is below tolerance; a non-finite gradient returns "gradient descent diverged; lower the learning rate" (L33–63). The header states the convergence condition, lr < 2 / λ_max, and says to standardize first.

LogisticRegression

Binary problems fit one weight row through a sigmoid, K > 2 classes fit K rows through a softmax. The objective is the summed negative log-likelihood plus ‖W‖²/(2C), scikit-learn's, so the optimum is the same and the intercept is unpenalized (logistic_regression.h L7–18). Parameters are packed one output at a time as [w_0 … w_{d−1}, b] against an augmented input [x, 1] (logistic_regression.cc L20–34).

1. Each iteration builds the full gradient and Hessian: residual p − y and curvature p(1 − p) (binary) or diag(p) − ppᵀ (softmax), plus 1/C on the weights and a 10⁻⁸ jitter on the whole diagonal, since a softmax's biases are only fixed up to a constant (L93–143).
2. The Newton step solves H·d = g by Cholesky; then Armijo backtracking halves t until the objective drops by 10⁻⁴·t·gᵀd, up to 50 halvings (L171–189).
3. It stops when the largest gradient entry is under tolerance, or when no step lowers the objective at all, which it treats as converged at double precision (L163–199). The loss uses max(z, 0) + log1p(e^−|z|) − yz, so separable data never overflows (L72–75).

coef() is [1, d] for two classes, the class-1 weights as in scikit-learn, and [K, d] otherwise (L49–54). The Hessian has (outputs·(d + 1))² entries, so the header scopes the class to small d.

Neighbors and clusters

k-nearest neighbors

Fit stores the training matrix and fails when k > n (knn.cc L63–78). Each query computes every squared distance, partial_sorts the k smallest with ties broken by training-row order, and takes square roots only of those (L21–39). Weights are uniform or 1/r; when the nearest neighbor is at distance 0, only the exact matches count, since their weight would be infinite (L41–55). Queries run under ParallelFor with a grain of 16 (L61); each writes only its own output row.

KMeans

• Seeding. KMeansPlusPlus picks the first center uniformly, then each next one with probability proportional to its squared distance from the nearest chosen center; with fewer distinct points than k it falls back to a uniform pick (kmeans.cc L104–140).
• Lloyd. Assign every point to its nearest center (lowest index on ties), move every center to its points' mean, stop when the summed squared center shift is at most the tolerance (L60–98). An empty cluster restarts at the point farthest from its center (L78–84).
• Tolerance and restarts. tolerance is relative: it is multiplied by the data's mean per-feature variance (L151–160). Restart r draws from Generator(seed, r), and the run with the least inertia wins (L144–174).

PCA and naive Bayes

Pca

Fit centers X and takes its SVD rather than eigendecomposing the covariance, which would square the condition number (pca.h L7–15). The components are the first c columns of V, the explained variance is s²/(n − 1), and the ratio divides by the total over all singular values (pca.cc L19–62). The sign convention that makes components match scikit-learn, largest-magnitude entry positive, lives in SvdFactor itself (svd.h L24). Transform projects and, with whiten, divides by √variance; InverseTransform undoes both (L64–104).

GaussianNb and MultinomialNb

Both fit by counting and score in log space, normalizing each row with log-sum-exp (naive_bayes.cc L30–45). Every class in 0..max(y) needs at least one sample, or Fit returns InvalidArgument (L20–28).

• Gaussian. Per-class means and population variances, plus var_smoothing times the largest variance of any feature over the whole data, so a feature that is constant within a class does not divide by zero (L57–104).
• Multinomial. Smoothed log-probabilities log(count + α) − log(total + dα); negative counts return InvalidArgument, and a zero count skips its term so α = 0 cannot give 0·log 0 (L134–184).

Trees

The tree code is split into a public pair of estimators and an internal grower that the forest reuses (decision_tree.h L48–89). A Tree is a flat vector of TreeNodes with the root at index 0; a node with feature = −1 is a leaf, and value holds class frequencies or the one mean target.

Growing a CART tree

1. Leaf test. A node is a leaf at max_depth, below min_samples_split, below twice min_samples_leaf, or with impurity ≤ 10⁻¹² (decision_tree.cc L191–193).
2. Sweep. For each candidate feature, sort the node's (value, row) pairs, then move rows one at a time from a right Stats to a left one, which keeps counts (or sum and sum of squares) up to date in O(1). Each boundary between distinct values that leaves both sides ≥ min_samples_leaf is scored by n_L·I_L + n_R·I_R, and the threshold is the midpoint (L128–158).
3. Zero-gain splits. The best split is taken even when it lowers impurity by nothing, which XOR needs before its second level pays off (decision_tree.h L15–18).
4. Recurse. std::partition the rows in place, add n·I − (n_LI_L + n_RI_R) to the feature's importance, and recompute a regression child's variance in two passes, so a constant child is exactly pure (L180–220).

Worked example. The first tree test fits x = 1, 2, 3, 10, 11, 12 with labels 0, 0, 0, 1, 1, 1 (decision_tree_test.cc L23–33). The root's Gini impurity is 1 − (0.5² + 0.5²) = 0.5. The sweep scores five boundaries:

Limits and feature subsets

Random forests

GrowForest grows tree t inside a ParallelFor with grain 1, drawing the bootstrap rows and the feature subsets from Generator(seed, t), so the forest depends on the seed and not on which thread grew which tree (random_forest.cc L22–43). Bootstrap rows may repeat, which GrowTree accepts. The classifier's default limits are {-1, 2, 1, 0, FeatureSubset::kSqrt} (random_forest.h L32); the regressor uses all features. Predictions average the trees' leaf values (class frequencies or means), and importances average each tree's normalized importances (L45–55).

Gradient boosting

Boosting fits F(x) = F₀ + Σ lr·f_t(x), each tree on the per-sample gradient g and Hessian h of the loss at the current scores. A leaf with sums G and H outputs the Newton step −G/(H + l2), and a split's gain is G_L²/(H_L + l2) + G_R²/(H_R + l2) − G²/(H + l2) (gradient_boosting.h L7–24). The shared loop is internal::Booster; the regressor uses squared error, the classifier log loss with one tree per round for two classes and one per class otherwise (L76–102).

Worked example. One round with lr 1, l2 2 and leaf size 1 on x = 0, 1, 2, 3, y = 0, 0, 10, 10 (gradient_boosting_test.cc L52–70): the baseline is the mean, 5, so g = F − y = 5, 5, −5, −5 and h = 1. The split at 1.5 has G_L = 10, H_L = 2 on the left, so that leaf outputs −10/(2 + 2) = −2.5 and the right one +2.5, and the predictions are 2.5 and 7.5. Its gain is 10²/4 + 10²/4 − 0 = 50; the code drops the factor ½ the header mentions, which does not change which split wins.

Support vector machine

SupportVectorClassifier is binary, labels 0 and 1 mapped to −1 and +1, and solves the soft-margin dual, min ½ Σ α_iα_jy_iy_jK_ij − Σ α_i subject to 0 ≤ α ≤ C and Σ α_iy_i = 0, by sequential minimal optimization (svm.h L8–27). The full n × n kernel matrix is computed once, which the header scopes to a few thousand points.

1. Kernel and gamma. Linear x·z, polynomial (γx·z + coef0)^degree or RBF exp(−γ‖x − z‖²); gamma = 0 resolves to 1/(d·Var(x)) over every entry, scikit-learn's "scale" (svm.cc L24–49).
2. Working set. The gradient Qα − 1 starts at −1 with α = 0. Each step picks i, the largest violator −y_tG_t among indices that can move up, and j, the index in the low set with the largest guaranteed decrease b²/a, LIBSVM's second-order rule, flooring the curvature a at 10⁻¹² (L81–128).
3. Stop. When the maximal violation is under tolerance (1e-3, LIBSVM's default); reaching max_iterations first returns Internal (L129–134).
4. Step. Move α_i and α_j along the line that keeps Σαy fixed, to the parabola's minimum clipped to the box, snap values within 10⁻¹²·C of a bound onto it, and update the gradient in O(n) (L136–153).
5. Bias and support. b averages over the free support vectors, or takes the midpoint of the interval the bounded ones allow when none is free; the support vectors (α > 0) and their α·y are kept for prediction (L156–194).

Checked against scikit-learn

tools/make_sklearn_fixtures.py writes rml/classical/testdata/sklearn.json: small datasets from a default_rng(0) stream, rounded through float32 so both sides see identical inputs, and what scikit-learn fits on them (make_sklearn_fixtures.py L27–29). The committed file records scikit-learn 1.9.1. sklearn_fixture_test reads it with nlohmann/json, so bazel test needs no Python (sklearn_fixture_test.cc L29–36, BUILD.bazel L275–292).

Random forests, gradient boosting and the tree classifier have no fixture: their results depend on random draws and tie-breaking that scikit-learn makes differently, so their own tests check properties instead.

Who uses it

Two callers outside the module. Retrieval's IvfIndex::Train clusters a sample into its inverted lists with KMeans at one restart and 50 iterations (ivf.cc L39–47), and the time-series example fits its seasonal autoregression with LinearRegression (time_series.cc L115–130). The inference stack leaves classical ML out of scope.

Tests