GridTensor

src/grid_tensor.rs and the src/grid_tensor/ module tree implement one fitted separable component — the rank-1 product in two-tensor form. This is the heart of the codebase.

The GridTensor struct

pub struct GridTensor {
    pub splits: Vec<Vec<f64>>,              // split thresholds per axis
    pub observation_counts: Vec<Vec<usize>>,// samples per interval per axis (for normalization)
    pub intervals: Vec<Vec<(f64, f64)>>,    // interval bounds per axis (±inf at the ends)
    pub scaling: f64,                       // legacy; ignored in two-tensor mode (= 1.0)
    pub backbone_values: Vec<Vec<f64>>,     // b_j ≥ 0  per interval per axis
    pub tilt_values: Vec<Vec<f64>>,         // d_j ∈ ℝ  per interval per axis
    pub lambda_plus: f64,                   // λ₊ ≥ 0
    pub lambda_minus: f64,                  // λ₋ ≥ 0
}

Prediction follows the model directly — the difference of the two non-negative branches, \(\hat{m}_+ - \hat{m}_- = \lambda_+\prod_j b_j e^{d_j} - \lambda_-\prod_j b_j e^{-d_j}\):

Method	Returns
new_two_tensor(...)	construct from the two-tensor fields (sets scaling = 1.0)
predict_single_unscaled(x) / predict_unscaled(X)	unscaled \(\tilde{m}_+ - \tilde{m}_-\) (one point / batch)
predict_single_backbone_and_tilt(x) / predict_backbone_and_tilt(X)	decompose into \((\prod_j b_j,\ \sum_j d_j)\)

Use the unscaled path

predict_unscaled returns \(\tilde{m}_+ - \tilde{m}_-\) without \(\lambda\) scaling beyond what is baked into the grid; the OLS stage scaling lives on StagePredictor. See the scaling invariant.

The fitting loop

grid_tensor::fit (src/grid_tensor/fit.rs) iterates the action/reducer pattern:

1. build a FittingState from (x, y);
2. initialize the RefinementStrategy (precomputes all candidate scores) and the SplitStrategy (builds allowed intervals);
3. optionally apply the histogram-binning mask (max_bins);
4. score eligible split, resplit, and merge actions under one fixed-scale cost-complexity objective, then apply the highest positive-gain action through the reducer;
5. at the n_iter fineness budget, exclude new splits while leaving resplits and merges eligible; a merge frees capacity for another split;
6. normalize with l2_identify and return the GridTensor.

With complexity_penalty=0.0, merge statistics are not maintained and fitting retains the split/resplit-only budget behavior.

FittingState (state.rs)

The mutable per-grid context threaded through the loop: the in-progress two-tensor grid (backbone_values, tilt_values, lambda_plus, lambda_minus); per-point caches (f_plus, f_minus, f, r_tilde = within-stage residuals); sorted row/rank mappings; and a PrecomputedStatistics cache of prefix sums and per-interval sufficient statistics enabling \(O(1)\) split-candidate scoring. Holds an optional LoggingState.

Actions (action.rs) and splitting (splitting.rs)

FittingAction enumerates the state mutations: ApplySplit, ApplyResplit, ApplyMerge, ApplyScaling, Terminate, Composite. The SplitStrategy enum proposes them:

• Random — sample up to split_try valid positions over colsample_bytree of the features (the default; usually best for speed);
• Best — greedily take the highest-gain split (better on small data);
• TopK — sample from the top-k candidates (must_fill_all_k controls strictness).

A SplitCandidate carries the column, position, parameter-regularized gain, and the left/right \((u_+, u_-)\) updates. A MergeCandidate identifies the removed boundary and its parameter-regularized gain; the absolute merged factors live in the aligned precomputed cache. SplitStrategyState tracks valid positions and the fixed boundary cost. After a split at position idx, the neighborhood is forbidden from further splitting.

The reducer (reducer.rs)

fitting_reducer matches on the action and applies it: it inserts the split into the grid, updates backbone/tilt on the two new sub-intervals, recomputes the affected per-point caches and the prefix-sum/interval statistics, and (under evo-logging) emits log events. For a merge, it removes one boundary, stores the candidate's absolute factors, and assigns the affected branch predictions directly from the other-axis partial products.

The two-tensor solver

solve_two_tensor (src/grid_tensor/two_tensor_solver.rs) is the load-bearing primitive. It builds a locally regularized \(2\times2\) proposal system from the five sufficient statistics \(S_{11}, S_{22}, S_{12}, t_1, t_2\),

\[ A = \begin{pmatrix}S_{11}&S_{12}\\S_{12}&S_{22}\end{pmatrix} +\frac14 \begin{pmatrix}\alpha+\tau&\alpha-\tau\\\alpha-\tau&\alpha+\tau\end{pmatrix}, \qquad t=\begin{pmatrix}t_1\\t_2\end{pmatrix}, \]

solves \(\hat{u}=A^{-1}t\) (or the three-case subgradient system when \(\rho>0\)), and clamps \(v_\pm = \mathrm{clamp}(1+\hat{u}_\pm, v_{\min}, v_{\max})\) to preserve positivity. The candidate score subtracts the exact log-coordinate cost \(\alpha\beta^2+\tau\delta^2+\rho|\delta|\), where \(\beta=\tfrac12(\log v_++\log v_-)\) and \(\delta=\tfrac12(\log v_+-\log v_-)\). The full derivation is in Fitting → the closed-form solver.

Refinement strategies (refinement.rs)

RefinementStrategy has two variants that differ only in the weights \(w_i\):

• L2Refinement — \(w_i = 1\);
• HuberRefinement — robust weights \(w_i = \min(1, c/|r_i|)\), \(c\approx1.345\).

Both expose alpha (log-backbone \(\ell_2\)), tilt_tau (tilt \(\ell_2\)), tilt_rho (tilt \(\ell_1\)), prior_sample_size, and update_clamp. initialize precomputes the penalized gain of every candidate position.

Identification

l2_identify (src/grid_tensor/identification.rs) applies the prediction-preserving normalization after a fit: it rescales each axis's backbone to unit weighted \(L^2\) (pushing the scale into \(\lambda_\pm\)) and centers each tilt by its weighted mean (pushing the offset into \(\lambda_+ \mathrel{*}= e^{\sum c_j}\), \(\lambda_- \mathrel{*}= e^{-\sum c_j}\)). It does not flip orientation (which would change predictions). The same routine runs before averaging bagged grids (Bagging & aggregation).

Histogram binning

apply_histogram_binning_mask (src/grid_tensor/histogram_binning.rs) is an optional prologue enabled by max_bins: Some(u16). It computes deterministic quantile bin edges per feature, forbids non-edge split positions, and replaces per-position prefix sums with per-bin cumulative sums. This cuts the per-candidate cost from \(O(n)\) to \(O(\text{bins})\) with little accuracy loss on larger datasets. With max_bins: None the exact path (every position a candidate) is used. Parity between the two paths is tested in tests/histogram_binning.rs.

Parameters (params.rs)

GridTensorParams bundles n_iter (the split budget), a SplitStrategyParams (split_try, colsample_bytree, min_interval_samples, min_split_loss, complexity_penalty, kind), a RefinementStrategyParams (alpha, tilt_tau, tilt_rho, …, kind), and max_bins. All have fluent …Builder types. The Python flat hyperparameters map onto these — see Hyperparameters.