Ordered Piecewise Additive Learner: A Hybrid Regression Algorithm Combining Linear Basis Functions and Regression Trees

1. Introduction

Algorithms based on gradient-boosted trees are a strong default for tabular regression [5,11,14]. A separate line of hybrid rule-based modeling has explored how tree-derived rule terms can be combined with simpler additive terms; RuleFit, for example, combines tree-derived rules with linear terms in a single regularized predictor [9]. One potential benefit of such hybrids is that they can pair the predictive strength of tree-based components with additional structure from other term families, such as additive or rule terms that are easier to inspect. OPAL (Ordered Piecewise Additive Learner) is motivated by a related question: can additive hinge and categorical terms be combined with residual regression-tree terms to retain competitive accuracy while exposing part of the fitted function through an inspectable structure?

Conceptually, OPAL seeks to improve post-fit transparency by separating Stage 1 additive terms on prepared columns from Stage 2 tree terms. Accordingly, the OPAL architecture defines two sequential stages of model fitting: Stage 1 adds prepared-feature hinge and categorical indicator terms to squared loss residuals, and Stage 2 adds regression-tree terms to the remaining residual structure. When enabled, a single-tree initializer is fit before these stages and its in-sample predictions are subtracted, so the two-stage model is trained on a residualized target. To guard against overfitting, OPAL then uses inner cross-validation to decide whether final predictions should be formed from the full hybrid predictor, the two-stage model alone, the optional single-tree initializer, the training-set mean, or a weighted blend of available predictions. The selected prediction scoring rule is stored with the fitted model and reused at prediction time.

The remainder of the paper evaluates whether this configured OPAL pipeline is both technically coherent and empirically competitive. First, it formalizes OPAL as a hybrid regression model whose Stage 1 terms include left and right hinge basis functions and categorical indicators, and whose Stage 2 terms add regression-tree corrections. Second, it describes the ordered residual-fitting procedure, including deterministic sample- and feature-aware budget shrinkage and the inner-CV prediction selector used for regularization. Finally, it evaluates OPAL on the 35-task CTR23 benchmark against XGBoost, LightGBM, and MARS, followed by a brief example of the interpretability afforded by the distilled additive term representation. The paper does not make broad interpretability claims for the full predictor; rather, it argues that OPAL’s Stage 1 additive component provides a more transparent, albeit partial, view of selected prepared-feature effects than traditional black-box tabular regression algorithms.

2. Related Work

OPAL sits at the intersection of additive spline and hinge models, rule-based hybrids, and modern tree ensembles. Generalized additive models (GAMs) [12,13,22] represent a response through sums of smooth feature functions, while multivariate adaptive regression splines (MARS) [8] builds adaptive spline regressions from piecewise-linear hinge basis functions. MARS therefore provides the hinge-basis lineage most directly connected to OPAL’s Stage 1 terms. RuleFit [9] combines tree-derived rules and linear terms in a single regularized predictor. OPAL differs from these approaches by using an ordered residual-fitting design: hinge and categorical indicator terms are fitted first, and regression-tree terms are then added to the residuals that remain.

GA2M models [16] and explainable boosting machines (EBMs) [18] extend the additive-model line by retaining inspectable main-effect functions while allowing selected pairwise interactions or boosted tree-shaped functions. Neural additive models (NAMs) [1] and GAMI-Net [23] replace some of the function-estimation machinery with neural subnetworks while preserving additive main effects and, in GAMI-Net, structured pairwise interactions. OPAL instead uses hinge and categorical terms in Stage 1, then uses regression-tree rules as a downstream residual term family.

The hinge and categorical terms retained by OPAL are auditable because each exposes its prepared feature, threshold or level, direction where applicable, and fitted weight. The position taken herein is a limited transparency claim rather than a claim of full-model interpretability, as generated pairwise columns can enter Stage 1, and the deployed predictor can also include regression-tree terms, an optional initial tree component, and an internally selected component blend.

3. The OPAL Framework

3.1. Feature Preparation and Deterministic Augmentation

OPAL first converts raw tabular inputs into a deterministic prepared matrix. Numeric predictors are coerced to numeric form and median-imputed. In the benchmark presented in a subsequent section, categorical predictors are converted to strings after replacing missing values with a reserved __MISSING__ token. If a training-split categorical column has more than max_cat_levels=1000 distinct levels, only the 999 most frequent training levels are retained and all other levels are mapped to __OTHER__; columns with at most 1000 levels keep all observed training levels. The mapped training column is then encoded with scikit-learn’s OrdinalEncoder using its fitted categories_ order. Evaluation levels not present after training-split mapping are assigned one code above the largest training code for that column. The numeric and categorical blocks are then merged into contiguous floating-point matrices. The implementation carries the categorical-column indices forward: Stage 1 treats those columns through equality indicators, and Stage 2 uses one-versus-level equality splits rather than ordered threshold splits on the ordinal codes.

The prepared matrix may be augmented with deterministic transformed numeric features; the enabled transform set is shown in Table 1. Only numeric columns participate in this construction. Remaining non-finite numeric values are mapped to zero, and each numeric column receives a deterministic continuous-target mutual information (MI) relevance score [4,15,19]. If MI scoring is unavailable, the implementation uses training variance as a fallback. For pairwise features, OPAL enumerates numeric-column pairs in increasing prepared-column order $(i,j)$ with $i<j$, ranks pairs by the sum of their two relevance scores, and relies on the stable sort to preserve this enumeration order for ties. The configured pair budget is therefore a budget on retained numeric pairs, not on final generated columns. In the benchmark this pair budget is 30. Each retained pair is expanded in the fixed order product, ratio, difference when those transforms are enabled; the ratio and difference are one-directional, $x_i/\left(\right|x_j|+{ϵ}_j)$ and $x_i-x_j$, using the same $i<j$ orientation. Unary transforms are added in a separate pass: OPAL ranks individual numeric columns by their own relevance scores, retains $min(d_{\mathrm{num}},max(5,\lfloor K/4\rfloor ))$ columns for pair budget K, and appends the enabled signed square-root, signed logarithm of one plus absolute value, and square transforms. With the benchmark value $K=30$, this unary pass considers seven numeric columns before the later variance and MI filters.

For numerical stability, generated columns are processed using training-split statistics before they enter the model. Product, difference, signed square-root, signed logarithm of one plus absolute value, and square columns are winsorized to the training 1st–99th percentiles. Ratio columns use the stabilized denominator ${ϵ}_j=max(0.1median|x_j|,{10}^{-6})$ and are then clipped to the same training percentile range.

After generation, transformed columns are standardized using training-set moments, clipped to $[-10,10]$, and stripped of near-constant columns. If generated columns remain, a final deterministic MI filter keeps columns whose generated-feature score is at least 10% of the maximum generated-feature score or belongs to the top half of generated columns. The surviving generated columns are appended to the prepared matrix. This is supervised, training-split feature preparation rather than a Stage 1 or Stage 2 fitted term family: the MI scores, retained generated columns, winsorization limits, and standardization moments are learned only inside the relevant training fold and then applied to that fold’s evaluation rows.

3.2. Prediction Function and Components

The stagewise main component of OPAL is an additive residual model over the prepared matrix. Before the final ridge calibration described below, it uses two high-level fitted term families, with categorical indicators included in the first stage:

$$\begin{array}{cc}\hfill {\tilde{f}}_{\mathrm{main}}\left(x\right)& ={\beta }_0+\underset{\mathrm{hinge}\mathrm{terms}}{\underbrace{{\sum }_{k\in \mathcal{H}}{w}_k{h}_{j_k,s_k,t_k}\left(x\right)}}+\underset{\mathrm{categorical}\mathrm{indicators}}{\underbrace{{\sum }_{m\in \mathcal{C}}{u}_m\mathbb{I}(x_{c_m}=v_m)}}+\underset{\mathrm{tree}\mathrm{terms}}{\underbrace{{\sum }_{t\in \mathcal{T}}{\gamma }_t{T}_t\left(x\right)}}.\hfill \end{array}$$

(1)

Here $\mathcal{H}$ indexes left and right hinge terms, $\mathcal{C}$ indexes categorical indicator terms selected during Stage 1, and $\mathcal{T}$ indexes regression-tree rules selected during Stage 2. The hinge basis $h_{j_k,s_k,t_k}\left(x\right)$ uses prepared-feature index $j_k$, threshold $t_k$, and direction $s_k\in \{\mathrm{left},\mathrm{right}\}$. The tree term $T_t\left(x\right)$ maps an input row to its assigned leaf value, and ${\gamma }_t$ is the corresponding ridge-scaled tree weight. All stagewise updates in the reported benchmark use squared loss on the standardized training target, so the negative gradients used for fitting are ordinary residuals. The weights in Equation (1) are the internal stagewise term weights before the final ridge calibration.

When the optional initial tree component is enabled, OPAL first standardizes the target, fits a single regression tree, and subtracts that tree’s predictions on the original training rows. The initializer uses the non-matrix fit path: for $n\le 2000$ it trains on a seeded permutation containing every training row exactly once, while for $n>2000$ it trains on a seeded with-replacement draw of size n. The fitted initializer is then evaluated on the original training and evaluation rows, and its original-training predictions are subtracted before fitting the stagewise main component in Equation (1). In the code this option is controlled by the pretrain_* configuration fields; the paper refers to it as the initial tree component or single-tree residual initializer.

After the stagewise main component has been fit, the canonical benchmark path applies a one-column ridge refit to calibrate the main-component prediction. Let $y_i^{★}$ denote the standardized target used by the main component: the standardized target itself when no initializer is used, and the residualized standardized target when the initializer is used. Let $p_i={\tilde{f}}_{\mathrm{main}}\left(x_i\right)$ be the uncalibrated main-component prediction on training row i. With $\overline{p}=n^{-1}{\sum }_i{p}_i$ and ${\overline{y}}^{★}=n^{-1}{\sum }_i{y}_i^{★}$, OPAL solves

$$\widehat{a}={\displaystyle \frac{{\sum }_i(p_i-\overline{p})(y_i^{★}-{\overline{y}}^{★})}{{\sum }_i{(p_i-\overline{p})}^2+{\alpha }_{\mathrm{refit}}}},\widehat{b}={\overline{y}}^{★}-\widehat{a}\overline{p},$$

(2)

using the fixed benchmark value ${\alpha }_{\mathrm{refit}}=1.0$. The calibrated main component is $g_{\mathrm{main}}\left(x\right)=\widehat{b}+\widehat{a}{\tilde{f}}_{\mathrm{main}}\left(x\right)$. This ridge refit does not introduce new basis functions and does not refit individual hinge, categorical, or tree terms; it only applies a global slope and intercept to the stagewise main-component prediction. Centered contribution curves for the main component can therefore be shown on the calibrated scale by multiplying the corresponding stagewise contributions by $\widehat{a}$.

What is referred to here as the full component includes the calibrated main component plus the configured scaled initial-tree contribution, while the main component excludes that initial tree. The benchmark subtracts the unscaled initializer during residualization but adds back 97% of it in the full component. This is a fixed implementation choice intended to damp the initializer when it is recombined with the main residual model; it is not tuned inside the benchmark and is not presented as an optimal value. The current experiments do not include an ablation or sensitivity analysis over this constant, so the empirical results presented herein should be read as results for the configured OPAL system with this scale fixed. The inner-CV selector can still reject the full component in favor of the main component, the initializer alone, the training-set mean, or their simplex blend.

3.3. Stage 1: Hinge and Categorical Terms

For a prepared numerical feature indexed by j and threshold t, OPAL considers both hinge directions:

$$\begin{array}{cc}\hfill h_{j,t}^{\mathrm{right}}\left(x\right)& =max(0,x_j-t),\hfill \\ \hfill h_{j,t}^{\mathrm{left}}\left(x\right)& =max(0,t-x_j).\hfill \end{array}$$

(3)

The candidate threshold grid for each prepared numerical feature is built once from the full training column before row subsampling is applied. Let $u_{j,1}<\dots <u_{j,L_j}$ denote the sorted unique finite values in that column and let $Q_{\mathrm{eff}}$ denote the effective value of n_quantiles after the internal capacity adjustment in Section 3.6. If $L_j\le Q_{\mathrm{eff}}$, all unique values are used as candidate thresholds. If $L_j>Q_{\mathrm{eff}}$, OPAL uses the linearly interpolated empirical percentiles of the unique-value vector at probabilities $\ell /(Q_{\mathrm{eff}}-1)$ for $\ell =0,\dots ,Q_{\mathrm{eff}}-1$, then removes duplicate threshold values. Categorical columns do not use hinge thresholds; they contribute equality-indicator candidates.

At Stage 1 iteration m, OPAL samples a row set ${\mathcal{R}}_m$ without replacement using the tuned subsample value. With n training rows and subsampling rate s, the sampled size is $\lfloor ns\rfloor$ when $0<\lfloor ns\rfloor <n$ and all n rows otherwise. It also samples a candidate feature set ${\mathcal{F}}_m$ using colsample_bytree; with p prepared features and column-sampling rate c, the sampled feature count is $max(1,\lfloor pc\rfloor )$ when $c<1$ and all p features otherwise. Candidate scoring scans only features in ${\mathcal{F}}_m$ and computes the score over rows in ${\mathcal{R}}_m$. Given current residual vector r, ridge parameter $\lambda$, and a candidate basis vector h, the unshrunk one-coordinate ridge coefficient and the associated regularized gain score are

$$a^{★}\left(h\right)={\displaystyle \frac{{\sum }_{i\in {\mathcal{R}}_m}{r}_i{h}_i}{{\sum }_{i\in {\mathcal{R}}_m}{h}_i^2+\lambda }},gain\left(h\right)={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\left({\sum }_{i\in {\mathcal{R}}_m}{r}_i{h}_i\right)}^2}{{\sum }_{i\in {\mathcal{R}}_m}{h}_i^2+\lambda }}.$$

(4)

This score is the one-coordinate ridge objective improvement used for candidate selection, not the unpenalized squared-loss reduction after shrinkage. The selected term receives the post-shrinkage weight $\eta a^{★}\left(h\right)$, where $\eta$ is the global learning rate for the stagewise update, and its prediction contribution is then applied to all training rows when residuals are updated. For categorical predictors, the candidate basis is $h_{j,v}\left(x\right)=\mathbb{I}(x_j=v)$; the same coefficient and gain formulas apply with $h_i\in \{0,1\}$. Stage 1 repeats this greedy residual-fitting step for the retained main-term budget, adding the best hinge or categorical indicator candidate at each successful iteration. If the best gain is non-positive, the implementation advances without adding a term. Algorithm 1 summarizes the Stage 1 update.

Algorithm 1:Stage 1: hinge and categorical residual fitting

3.4. Stage 2: Regression-Tree Rule Terms

After Stage 1, OPAL fits regression-tree rule terms to the remaining residuals. At tree iteration t, OPAL draws a row subset ${\mathcal{R}}_t^{\prime }$ and a tree feature subset ${\mathcal{F}}_t^{\prime }$ using the same sample-size rules as Stage 1. The tree is fit on ${\mathcal{R}}_t^{\prime }$ and is allowed to split only on features in ${\mathcal{F}}_t^{\prime }$. In the benchmarked histogram tree path, numeric prepared features reuse the same global quantile grids described above, with $Q_{\mathrm{eff}}$ thresholds per feature before duplicate removal; at each node OPAL evaluates the available grid thresholds and accepts only splits satisfying the minimum-leaf constraint on the sampled rows in that node. Categorical prepared features use equality splits of the form $x_j=v$ versus $x_j\ne v$, so their ordinal codes are not treated as ordered quantities inside Stage 2. Each tree greedily chooses splits that reduce within-node residual variance, subject to the configured depth and minimum-leaf constraints. Leaf values are the mean residuals among sampled rows in ${\mathcal{R}}_t^{\prime }$ that reach the leaf, not the mean over all training rows assigned by the final split structure. Once a tree output vector T is grown, OPAL applies a scalar ridge rescaling computed on the sampled rows,

$${\alpha }^{★}\left(T\right)={\displaystyle \frac{{\sum }_{i\in {\mathcal{R}}_t^{\prime }}{r}_i{T}_i}{{\sum }_{i\in {\mathcal{R}}_t^{\prime }}{T}_i^2+\lambda }},\gamma =\eta {\alpha }^{★}\left(T\right),$$

(5)

then applies the scaled tree to all training rows and updates the residuals. The tree stage stops growing individual trees at the configured maximum depth, when too few samples remain in a node, or when no valid split is available. Algorithm 2 summarizes the Stage 2 update.

Algorithm 2:Stage 2: regression-tree residual fitting

3.5. Internal Prediction Regularization

As a final regularization layer, OPAL uses inner cross-validation to choose among the full model, the initial tree component, the main component, the target-mean baseline, or a simplex-weighted blend of these predictions; the selected rule is stored as fitted model state. For each candidate configuration, OPAL fits the model on each inner-training fold, predicts the available components on the corresponding inner-validation fold, and stacks those out-of-fold component predictions against the inner-validation targets. Candidate components and the simplex stack are scored by RMSE on the pooled out-of-fold validation rows. The component selector or stack is learned and scored on the same pooled out-of-fold prediction matrix, so this value is a model-selection criterion for OPAL rather than an unbiased estimate of the selected rule’s generalization error. The selected rule is carried with the selected hyperparameter candidate and reused when the model is refit on the full outer-training split. It is therefore internal model selection/regularization rather than post-hoc tuning on the outer test data.

For auditing, OPAL materializes a full distilled additive prediction equation after the prediction rule has been selected. This equation applies the stored rule’s component weights to the fitted intercept, initializer, Stage 1 terms, and Stage 2 terms, so every displayed term contribution respects the deployed prediction rule. If the selected rule is the training-set mean or the initializer alone, the Stage 1 and Stage 2 terms receive zero weight in the distilled equation; if the selected rule is a simplex blend, each term is scaled by the corresponding component weight. Thus effect curves are summaries of terms that survive in the rule-weighted distilled equation, not explanations of components that the selector has excluded.

When the stack is selected, the code fits nonnegative component weights that sum to one and uses no free intercept, so the stacked prediction remains a convex combination of the available component predictions. Otherwise, the stored rule selects a single component by name. The target-mean component is included as a conservative baseline option, and the main component specifically denotes the Stage 1 plus Stage 2 fit after any residual initialization.

3.6. Budget Shrinkage and Stopping

OPAL applies deterministic capacity regularization to the provisional Stage 1 and Stage 2 iteration budgets. The rule is GCV-inspired but is not stated here as a claim of strict GCV optimality. For a provisional stage budget $n_s$, training size n, prepared feature count p, prune quantile $q_s$, and prune scale $c_s$, define $n_{\mathrm{eff}}=max(2,n)$ and $p_{\mathrm{eff}}=max(1,p)$ and set

$${\rho }_s={\displaystyle \frac{n_s}{max(1,n_{\mathrm{eff}}-1)}},\pi ={\displaystyle \frac{p_{\mathrm{eff}}}{n_{\mathrm{eff}}}},{\zeta }_s={\displaystyle \frac{1}{1+c_s({\rho }_s^2+0.25\pi )}}.$$

(6)

The retained budget is

$${\tilde{n}}_s=\left\{\begin{array}{cc}0,\hfill & n_s=0,\hfill \\ min\{n_s,max(1,round\left(n_s{q}_s{\zeta }_s\right))\},\hfill & n_s>0.\hfill \end{array}\right.$$

(7)

The feature count p is the number of columns in the prepared matrix after deterministic augmentation, near-constant generated-column removal, and the generated-feature MI filter. In the benchmark code, the stage1_prune_* controls apply to the Stage 1 hinge/indicator budget n_main_estimators, while the stage2_prune_* controls apply to the Stage 2 tree budget n_interact_estimators. The benchmarked configuration leaves the enhanced early-stopping mode disabled, so reported fits use fixed retained budgets after deterministic shrinkage.

The OPAL algorithm also applies a fixed small-sample capacity prior after the above budget shrinkage. This prior is not tuned by Optuna. With training size n, prepared feature count p, effective rows-per-feature ratio $n/p$, reference size 2000, and reference rows-per-feature value 10, define

$$a={\left(min\{1,n/2000\}\right)}^{1/4},g=min\{1,10p/n\},\delta =(1-a)g,\varphi =0.25+0.75(1-\delta ).$$

(8)

For the effective quantile count $Q_{\mathrm{eff}}$, OPAL sets

$$Q_{\mathrm{eff}}=min\{n,Q,max(10,round\left(Q\varphi \right))\},$$

(9)

with a final lower bound of one when n is very small. For the minimum leaf size $m_{\ell }$, it computes $L=min\{n,max(2,\lfloor n/2\rfloor \left)\right\}$ and raises the leaf size to $max\{m_{\ell },min(L,round\left(m_{\ell }(1+\delta )\right))\}$. If $\lambda <10$, the ridge parameter is moved toward 10 by replacing it with $\lambda +\delta (10-\lambda )$; otherwise it is unchanged. Finally, the retained Stage 2 budget is replaced by $min\{{\tilde{n}}_2,max(0,round\left({\tilde{n}}_2\varphi \right))\}$. The Stage 1 budget is not further reduced by this prior.

3.7. Workflow Summary

OPAL training proceeds in a fixed order. The estimator first prepares the numeric and categorical predictors and appends any enabled deterministic transformed features. If the initial tree component is enabled, a single tree is fit to the standardized target and the main component is trained on the residualized standardized target. The retained Stage 1 and Stage 2 budgets are then set by deterministic capacity shrinkage and the fixed small-sample prior. Stage 1 adds hinge and categorical indicator terms to residuals; Stage 2 adds regression-tree rule terms to the remaining residuals. The fitted main-component prediction is then ridge-calibrated with a single slope and intercept. Finally, OPAL forms the available component predictions, applies the stored inner-CV prediction rule, and de-standardizes the selected or blended prediction. Algorithm 3 lists the sequential workflow, and Figure 1 summarizes the same flow visually.

Algorithm 3:Sequential OPAL training and prediction workflow

4. Experimental Results

4.1. Benchmark Setup

The benchmark uses the 35 OpenML-CTR23 tabular regression tasks of Fischer et al. [7], accessed through OpenML [21] and the task-provided outer splits. Appendix Table A1 lists the data-set-level RMSE results. The present empirical comparison includes exactly four models: OPAL, XGBoost, LightGBM, and MARS. Five tasks expose 10 repeats, 10 folds, and one sample, giving 100 outer splits per task. The remaining 30 tasks expose one repeat, 10 fold indices, and one sample, giving 10 outer splits per task. The size-regime labels in Figure 2 are the nominal CTR23 labels; even tasks in the large-data/holdout regime are evaluated here using the OpenML task-provided split indices rather than a newly generated 33% holdout split. Figure 2 summarises the suite by OpenML row count and raw predictor count.

Within each outer training split, the benchmark driver uses deterministic inner tuning splits chosen by sample size: 10-fold inner cross-validation when $n_{\mathrm{train}}<1000$, 3-fold inner cross-validation when $1000\le n_{\mathrm{train}}\le 10000$, and a 33% holdout when $n_{\mathrm{train}}>10000$. Each model is tuned independently on each outer training split by an Optuna RandomSampler seeded from the split seed, with n_trials$=100$ and n_jobs$=1$ in the canonical benchmark invocation. XGBoost, LightGBM, and MARS minimize the unweighted mean RMSE across the deterministic inner splits. OPAL minimizes the pooled out-of-fold RMSE of its internally selected component rule for the candidate configuration, because the candidate includes the inner-CV component selector described in Section 3.5. These objectives are therefore not algebraically identical; the difference is part of OPAL’s internal model-selection procedure, and the OPAL inner score is not used as a benchmark performance estimate. The comparative results use only the common outer-split RMSE values. After optimization, the best completed trial is refit on the full outer training split and evaluated once on the corresponding outer test rows. For each data set and model, the reported data-set-level RMSE is the unweighted arithmetic mean of the resulting outer-split RMSE values; ranks and cross-data-set statistical tests use these data-set-level means rather than pooled out-of-fold RMSE or test-size-weighted RMSE. Table 2 lists the nine OPAL hyperparameters tuned in the benchmark.

For OPAL, only the nine user-exposed tuning parameters in Table 2 are optimized. The displayed n_main and n_interact ranges are sampled provisional budgets after the fixed benchmark estimator-scale factor is applied; deterministic capacity shrinkage is applied after sampling. The initial tree component, its 0.97 full-component scale, the one-column ridge-refit penalty ${\alpha }_{\mathrm{refit}}=1.0$, transformed-feature construction, and the transformed-feature pair budget $K=30$ are fixed. The enabled transforms are product, ratio, sqrt, log, square, and diff; as described in Section 3.1, $K=30$ means at most 30 ranked numeric pairs before expansion, with seven unary columns considered before filtering. The final number of appended generated columns can be smaller after near-constant removal and the generated-feature MI filter. The prune quantiles are both 1.0, and the prune scale parameters are both 12.0. The sampled raw Stage 1 and Stage 2 estimator budgets are first shrunk deterministically by the rule in Section 3.6; the canonical implementation path then applies the fixed small-sample capacity prior defined in the same section, which can further change the effective n_quantiles, min_leaf, $\lambda$, and Stage 2 budget. The optional internal early-stopping path is disabled. XGBoost and LightGBM are likewise fit without validation-set early stopping, and MARS has no learner-specific early-stop mechanism in the harness. Table 3 gives the corresponding baseline search spaces.

Preprocessing is deterministic within each split, as described in Section 3.1. For XGBoost and LightGBM, numeric columns are median imputed and categorical columns use the same missing-token and 1000-level rare-category map before one-hot encoding with unseen levels ignored. MARS receives a dense numeric and ordinal representation built from the same split-local imputation and rare-level handling as OPAL. The MARS baseline uses the in-repository py-earth implementation [20]. In addition to the tuned knobs in Table 3, the wrapper enables use_fast=True when that constructor argument is available and, unless overridden, sets fast_K=25 and fast_h=1; no benchmark override changes those fast-mode settings.

Implementation-wise, OPAL is exposed through a scikit-learn compatible Python estimator. Aside from that wrapper and benchmark orchestration, the core fit/predict path and much of the term-generation logic are implemented in Mojo, a programming language developed by Modular for high-performance systems and numerical code [17]. The reported run uses base seed 1, jobs=4, and optuna_jobs=1. For OPAL, the split seed is the stable hash of the base seed, task ID, dataset ID, OPAL’s split-seed stream label, the outer split index, and the string split. For XGBoost, LightGBM, and MARS, the split seed is the stable hash of the base seed, task ID, dataset ID, model name, outer split index, and split. Inner split seeds are stable hashes of the split seed and inner, and Optuna sampler seeds are stable hashes of the split seed and optuna. OPAL also uses a fixed additional stable-seed stream for the initial tree component.

The timing columns report wall-clock measurements from this local invocation, not portable speed claims. The run was executed on an Apple M3 Max system with 16 logical cores and 128 GiB RAM under macOS arm64. The Python environment used Python 3.9.6, NumPy 2.0.2, pandas 2.3.3, SciPy 1.13.1, scikit-learn 1.6.1, Optuna 4.8.0, XGBoost 2.1.4, LightGBM 4.6.0, and OpenML 0.15.1. Because OPAL uses a Mojo core while the baselines rely on external library implementations, the reported tuning times should be interpreted as run-specific implementation evidence rather than portable algorithmic speed comparisons.

Performance is measured by RMSE, and ranks are assigned within each data set so that rank 1 denotes the lowest RMSE. All four models produced results on all 35 CTR23 tasks, so the aggregate summaries and complete-case significance tests use the same 35 data sets.

4.2. Statistical Analysis

Average RMSE rank is the primary omnibus summary. We use a Friedman omnibus test over the four current models and then a Nemenyi critical-difference view in the standard multi-data-set comparison tradition [6,10]. Because mean-rank post-hoc procedures have known limitations [2,3], and because the OPAL–boosting gaps are small, we also report paired Wilcoxon signed-rank tests on per-data-set log RMSE ratios, $log({\mathrm{RMSE}}_{\mathrm{OPAL}}/{\mathrm{RMSE}}_{\mathrm{comparator}})$. This normalized test is invariant to target-scale changes across data sets. The Wilcoxon values reported in Table 5 are two-sided unadjusted p values; Holm-adjusted values are reported in the prose where they affect interpretation.

4.3. Aggregate Accuracy Results

Table 4 summarises the aggregate results. OPAL has the best average RMSE rank among the four evaluated methods: 1.8571 for OPAL, 2.0286 for XGBoost, 2.4286 for LightGBM, and 3.6857 for MARS. For each data set, the RMSE ratio divides a model’s RMSE by the lowest RMSE achieved by any of the four evaluated models on that data set; the table reports the geometric mean of those ratios across data sets. OPAL has the lowest geometric mean RMSE ratio, 1.0650, and the most first-place finishes, 17 of 35 tasks. XGBoost remains close, with the most top-two finishes (26 versus 24 for OPAL). In this run, MARS has the lowest median tuning time because its search space is lightweight, but it has the weakest average rank. OPAL’s recorded median tuning time is lower than the two boosted-tree baselines (440.9 s versus 595.8 s for XGBoost and 612.5 s for LightGBM), with the environment-specific caveat described above.

The OPAL/XGBoost comparison is close: OPAL has lower RMSE on 19 of 35 tasks, while XGBoost has lower RMSE on 16 of 35 tasks. The two-sided paired Wilcoxon test on log RMSE ratios gives $p=0.377$, so the two methods are not statistically distinguishable in this benchmark. OPAL has lower RMSE than LightGBM on 23 of 35 tasks; the log-ratio Wilcoxon value is $p=0.018$ and remains below 0.05 after Holm adjustment ($p=0.036$), although the Nemenyi rank comparison below does not separate the two methods. OPAL has lower RMSE than MARS on 33 of 35 tasks, with a much smaller log-ratio Wilcoxon value ($p=5.82\times {10}^{-10}$).

Average rank compresses meaningful variation in per-task placement. Figure 3 reports the full within-task rank distribution. OPAL and XGBoost occupy most of the upper ranks, LightGBM remains competitive but more often falls behind the top pair, and MARS is concentrated in fourth place.

The four-model Friedman test rejects equality of ranks (${\chi }_F^2=42.977$, $p=2.49\times {10}^{-9}$), and the Nemenyi critical distance is 0.793 at $\alpha =0.05$. OPAL and XGBoost differ by about 0.171 average-rank units and are not significantly different. OPAL and LightGBM differ by about 0.571 average-rank units, also below the critical distance, while OPAL and MARS differ by about 1.829 average-rank units, exceeding the critical distance. OPAL is effectively tied with XGBoost, has a scale-normalized paired advantage over LightGBM, and is stronger than the configured MARS baseline.

4.4. Stage 1 Interpretability Example

Figure 4 illustrates the narrow interpretability claim supported by OPAL’s term representation.

The image depicted in Figure 4 was generated from split 9 of OPAL’s fit on the cps88wages data set, where the stored prediction rule selected the full component and therefore gives the main Stage 1 terms nonzero weight in the distilled equation. Education was the top original Stage 1 hinge feature for this split, with 39 retained raw-feature hinge terms and thresholds at 13, 14, and 15 years of education. The curve therefore sweeps education across its observed range and uses only original-feature Stage 1 education terms; generated pairwise features, the initial tree component, and Stage 2 tree-rule corrections are excluded. The vertical axis is the centered Stage 1 contribution. Although this is not a full explanation of every prediction, it provides an auditable raw-feature slice that ordinary boosted-tree ensembles do not expose natively.

5. Conclusions

The present work establishes OPAL as a competitive two-stage hinge-and-tree learner for tabular regression. Deterministic preprocessing and optional transformed-feature construction prepare the design matrix; Stage 1 fits hinge and categorical indicator terms to residuals; Stage 2 fits regression-tree rules to the remaining residual structure. The current implementation also includes an optional initial tree component and an inner-CV prediction selector that chooses among model components or a simplex-weighted blend as fitted model state.

On the present 35-task CTR23 benchmark, OPAL attains the lowest average RMSE rank among OPAL, XGBoost, LightGBM, and MARS. The margin over XGBoost is small and not statistically significant. OPAL does not perform substantively better than LightGBM, as evidenced by the Nemenyi rank comparison, but the paired log-ratio Wilcoxon test favors OPAL on normalized RMSE differences. OPAL also improves over MARS under both views, but that conclusion is conditional on the MARS tuning range, dense numeric/ordinal representation, and fixed py-earth fast-mode settings used here [20]. The comparison should also not be read as a component ablation of OPAL; sensitivity to fixed implementation choices such as the 0.97 initializer scale remains a separate question that should be addressed through further research. The hinge terms remain directly auditable in the rule-weighted distilled equation, but the main contribution of this work is evidence for the configured OPAL pipeline rather than a broad claim about full-model interpretability. Overall, the results position OPAL as a strong alternative that sits between fully interpretable models and high-performance blackbox models, as it offers competitive performance with additional options for explainability using the Stage 1 main effect terms.