Sum of hessians calculated incorrectly under hist method

[david-cortes]

When running with the hist tree method which is now default, if one manually splits the data according to the information in the trees, one would see that the sum of hessians would not always match with what one would get by passing observations down the trees.

In the case of squared error, the hessians are always equal to one, so their sum is equal to the number of rows, which makes it easy to make the calculations.

Example where there would be some slight mismatches:

import xgboost as xgb
import treelite
from sklearn.datasets import fetch_california_housing
import numpy as np
import polars as pl
import json

tmp = fetch_california_housing()
X = pl.DataFrame(tmp.data, schema=tmp.feature_names)
y = tmp.target

model = xgb.train(
    dtrain=xgb.DMatrix(X, y),
    num_boost_round=1,
    params={
        "objective": "reg:squarederror",
        "seed": 123,
        "tree_method": "hist",
    }
)
tl_model = treelite.frontend.from_xgboost(model)
tl_json = json.loads(tl_model.dump_as_json())

def assign_nrows_to_json():
    curr_node = 0
    curr_X = tmp.data
    assign_nrows_to_json_recursive(curr_node, curr_X, tl_json["trees"][0]["nodes"])


def assign_nrows_to_json_recursive(
    curr_node: int,
    curr_X: np.ndarray,
    tl_json_nodes: dict,
):
    node_info = tl_json_nodes[curr_node]
    node_info["n_rows"] = curr_X.shape[0]
    if "leaf_value" in node_info:
        return
    
    val = curr_X[:, node_info["split_feature_id"]]
    thr = node_info["threshold"]
    if node_info["comparison_op"] == "<":
        mask_left = val < thr
    elif node_info["comparison_op"] == "<=":
        mask_left = val <= thr
    elif node_info["comparison_op"] == ">":
        mask_left = val > thr
    elif node_info["comparison_op"] == ">=":
        mask_left = val >= thr
    
    left_X = curr_X[mask_left]
    right_X = curr_X[~mask_left]
    assign_nrows_to_json_recursive(
        node_info["left_child"],
        left_X,
        tl_json_nodes,
    )
    assign_nrows_to_json_recursive(
        node_info["right_child"],
        right_X,
        tl_json_nodes,
    )

assign_nrows_to_json()

nodes_df = pl.DataFrame(tl_json["trees"][0]["nodes"])

display(
    nodes_df
    .filter(
        pl.col("sum_hess") != pl.col("n_rows")
    )
)

shape: (44, 12)

node_id	split_feature_id	default_left	node_type	comparison_op	threshold	left_child	right_child	sum_hess	gain	n_rows	leaf_value
i64	i64	bool	str	str	f64	i64	i64	f64	f64	i64	f64
5	5	false	"numerical_test_node"	"<"	2.746993	11	12	3056.0	506.611816	3057	null
6	0	false	"numerical_test_node"	"<"	7.8239	13	14	1335.0	222.381836	1334	null
7	0	false	"numerical_test_node"	"<"	2.2174	15	16	3436.0	245.249878	3441	null
8	0	false	"numerical_test_node"	"<"	2.4167	17	18	4617.0	191.924561	4612	null
11	1	false	"numerical_test_node"	"<"	20.0	23	24	1270.0	119.098389	1271	null
…	…	…	…	…	…	…	…	…	…	…	…
104	null	null	null	null	null	null	null	205.0	null	206	0.046349
109	null	null	null	null	null	null	null	388.0	null	387	0.161952
123	null	null	null	null	null	null	null	3.0	null	4	-0.159801
124	null	null	null	null	null	null	null	15.0	null	14	0.511588
125	null	null	null	null	null	null	null	518.0	null	517	0.812895

If one were to change tree_method to exact in the code above, it would return an empty table as expected.

[trivialfis]

Hi, XGBoost converts data to float32 internally. You can obtain the same result by converting the data to float32 first: curr_X = tmp.data.astype(np.float32).

Also, XGBoost supports saving the model in JSON in save_model and save_raw.

[david-cortes]

Hi, XGBoost converts data to float32 internally. You can obtain the same result by converting the data to float32 first: curr_X = tmp.data.astype(np.float32).

Also, XGBoost supports saving the model in JSON in save_model and save_raw.

Thanks. But then curious as to how this would happen.

Is the histogram algorithm using exact values in the data as borders instead of midpoints between them?

[trivialfis]

Is the histogram algorithm using exact values in the data as borders instead of midpoints between them?

yes.

[david-cortes]

Wouldn't it be more logical to use midvalues like in the exact tree method? I imagine it should in theory generalize better in real-valued distributions, and would make it harder to run into cases like this when using the trees outside xgboost.

[trivialfis]

Wouldn't it be more logical to use midvalues

One could add quantile interpolation, but I'm not sure how much it would benefit decision trees. Using the exact data point so far works, and it is easy to "debug" the tree. I don't think there's a strong opinion on this.

[RAMitchell]

Actually @david-cortes might be right. Last week I debugged a tiny accuracy gap in cuml random forest between sklearn. I could resolve the gap by taking midpoints. The key in the histogram algorithm is that many bins can be completely empty in a node. It may generalise better to take the midpoint of the empty bin range for continuous data. I need to experiment more.

[trivialfis]

Taking the midpoint feels like shifting the bin boundaries. If your goal is to shift the boundary, you don't need the midpoint; it could be an arbitrary point between the samples as long as it's consistent.

And the midpoint is not necessarily consistent, given that the bin sizes are different.

[RAMitchell]

My idea is to take the middle empty bin because this represents the median of the empty range of training data. It can be an arbitrary point but this seems reasonable to me.

[david-cortes]

And what about cases where there's no empty bin? It would make sense to take either the midpoint or an interpolation based on number of observations on each bin.

[RAMitchell]

If there is no empty bin you cannot change the split without potentially moving training data from one side to the other.

[david-cortes]

I don't follow. The bins each have boundaries. The boundary between two bins doesn't need to correspond with values taken as-is from the data, can also be interpolated (by midpoint or some other method) between the maximum and minimum of two adjacent bins.

[RAMitchell]

I see what you mean - yes you dont need to use data points as bin boundaries. We could in theory change the quantile cuts to the midpoint of the two boundary values. Maybe there is a benefit. It could be easily tested on a synthetic normal/uniform distribution dataset and see if the test error improves.

[RAMitchell]

I've just realised the quantile sketch algorithm doesn't store these neighbouring data points. It stores a very small number of cut points and throws the rest away. So we would need a modified algorithm.

[trivialfis]

I don't follow the motivation of doing this. It changes the interpretation of the bin boundaries.

[david-cortes]

I don't follow the motivation of doing this. It changes the interpretation of the bin boundaries.

Yes, that's the idea.

At a high level, you want xgboost to create decision rules that generalize better to new data, and if data is real-valued and the reweighted target (gradient / hessian) is symmetric, it's likely that interpolation between points would make up for a better rule.

However, I'm now having second thoughts about using midpoints in the case of non-empty bins, since the reweighted 'target' to predict at a given iteration (gradients / hessians) might not be expected to follow a symmetric distribution for objectives like reg:gamma, and it might be more expectable to find right-skewed distributions than left-skewed ones in practice.

In the most extreme case, to illustrate the mathematics behind it, assume x=y, and only two bins (= one boundary) and leafs determined from two observations. If decision rules are made up from random samples like this, using exact values of points for boundaries will yield worse predictions than using midpoints.

Example with a normal distribution:

import numpy as np
rng = np.random.default_rng(seed=123)

x = rng.standard_normal(size=int(1e5))

num_simulations = int(1e3)
sim_err_exact_lower = np.empty(num_simulations) # rule: x_i <= x_t
sim_err_exact_higher = np.empty(num_simulations) # rule: x_i >= x_t
sim_err_midpoint = np.empty(num_simulations) # rule: x_i <= (x_{t} - x_{t_1}) / 2

for sim in range(num_simulations):
    x_lower_bin = rng.standard_normal()
    x_higher_bin = rng.standard_normal()
    x_midpoint = 0.5*x_lower_bin + 0.5*x_higher_bin

    sim_err_exact_lower[sim] = np.where(x <= x_lower_bin, x_lower_bin, x_higher_bin).mean()
    sim_err_exact_higher[sim] = np.where(x >= x_higher_bin, x_higher_bin, x_lower_bin).mean()
    sim_err_midpoint[sim] = np.where(x <= x_midpoint, x_lower_bin, x_higher_bin).mean()

def mean_sq_resid(x):
    return np.dot(x, x) / x.shape[0]

print(
    f"Error with lower threshold: {mean_sq_resid(sim_err_exact_lower)}\n"
    f"Error with higher threshold: {mean_sq_resid(sim_err_exact_higher)}\n"
    f"Error with midpoint threshold: {mean_sq_resid(sim_err_midpoint)}"
)

Error with lower threshold: 0.72171539353638
Error with higher threshold: 0.7505254313629268
Error with midpoint threshold: 0.5858239917142031

But this depends on the distribution of the data - for example, using a Gamma(1) distribution which is skewed towards one side would yield instead:

Error with lower threshold: 2.3695836915722843
Error with higher threshold: 0.9622682648816496
Error with midpoint threshold: 1.591325959002651

And if the distribution were skewed towards the other side, then lower threshold would be better.

Perhaps midpoints could be a theoretically better choice for reg:squarederror and binary:logsitic, but since xgboost is expected to operate on more objectives, maybe the current logic is better.

That being said, in the case of empty bins in a histogram, it would indeed make more sense to take the decision boundary from the middle bin.