This looks a bit lax. Any idea why we can be stricter in the case of the squared error? I wouldn't require to make it stricter as a pre-requisite to merge this PR as the Poisson loss is already in main but lightgbm equivalence was previously just not tested.

But we could at least add a TODO comment to inform the reader that we are not 100% satisfied with the state of things :)

Maybe this is because, in the context of the squared error, we only test for the equivalence with the larger dataset with shallow trees (max_leaf_nodes < 10 and n_samples >= 1000): this condition might allow the averaging effect in leafs to mitigate the impact of rounding errors and maybe reduce the likelihood of ties when deciding on which features and thresholds to split on?

Maybe, LightGBM's algo deviates more from ours for general losses (=not squared error). For instance in the Poisson case, LightGBM has a poisson_max_delta_step which I set to 1e-12, but I'm not 100% sure of the effective difference this causes.

I'll add a comment as suggested.

And BTW, our losses have a better test suite than the ones in LightGBM:smirk:

Currently, this test fails because the HGBT with poisson loss out-performs all other models on the test set.
I observe the exact same behaviour with (often) similar number with LightGBM:

from lightgbm import LGBMRegressor

lgbm_gamma = LGBMRegressor(objective="gamma")
lgbm_pois = LGBMRegressor(objective="poisson")

I don't know how to construct a better test - for now.

Is this true only for the small sample sizes that we typically use in our tests our would be this still true for a very large n_samples?

Is it similar to the fact that the median can be a better estimate of the expected value than the sample mean for a finite sample of a long tailed distribution (robustness to outliers)? This was observed empirically when computing metrics on a left out validation set here: https://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_quantile.html

I played with n_sample and n_features without success. This problem that Poisson HGBT always beats Gamma HGBT out-of-sample prevails. I think one reason is that the log-link is not canonical to Gamma deciance, but it is canonical to the Poisson deviance. Therefore, the Poisson HGBT is better calibrated, which here also translates to a better out-of-sample performance. Mabye, the features are not informative enough in this example.

I can try to generate a decision tree and use it to generate a Gamma distributed sample per tree node. That should be easier to fit for tree based models than this GLM based data generation.

I can try to generate a decision tree and use it to generate a Gamma distributed sample per tree node. That should be easier to fit for tree based models than this GLM based data generation.

That would be interesting to know indeed.

I tried to create gamma distributed data with a decision tree structure for the mean value, see details.
This gives the same result: poisson loss beats gamma loss when comparing out-of-sample gamma deviance.

Proposition for solutions:

1. Remove this test and add (and rely on) a test comparing with LGBT.
2. Do not compare with poisson loss, just with squared error (which depends a lot on the clipping constant!)

I'm in favor of the first point.

import numpy as np
from sklearn._loss import HalfGammaLoss
from sklearn.dummy import DummyRegressor
from sklearn.datasets import make_low_rank_matrix
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_gamma_deviance
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeRegressor

from lightgbm import LGBMRegressor


### Data Generation
rng = np.random.RandomState(42)
n_train, n_test, n_features = 500, 500, 20
X = make_low_rank_matrix(
    n_samples=n_train + n_test, n_features=n_features, random_state=rng, 
)
# First, we create a normal distributed data
coef = rng.uniform(low=-10, high=20, size=n_features)
y_normal = np.abs(rng.normal(loc=np.exp(X @ coef)))

# Second, we fit a decision tree
dtree = DecisionTreeRegressor().fit(X, y_normal)

# Finally, we generate gamma distributed data according to the tree
# mean(y) = dtree.predict(X)
# var(y) = dispersion * mean(y)**2
dispersion = 0.5
y = rng.gamma(shape=1 / dispersion, scale=dispersion * dtree.predict(X))

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=n_test, random_state=rng
)


### Model fitting
params = {"early_stopping": False, "random_state": 123}
p_lgbt = {"min_child_weight": 0, "random_state": 123}

gbdt_gamma = HistGradientBoostingRegressor(loss="gamma", **params)
gbdt_pois = HistGradientBoostingRegressor(loss="poisson", **params)
gbdt_ls = HistGradientBoostingRegressor(loss="squared_error", **params)
lgbm_gamma = LGBMRegressor(objective="gamma", **p_lgbt)
lgbm_pois = LGBMRegressor(objective="poisson", **p_lgbt)
for model in (gbdt_gamma, gbdt_pois, gbdt_ls, lgbm_gamma, lgbm_pois):
    model.fit(X_train, y_train)
dummy = DummyRegressor(strategy="mean").fit(X_train, y_train)

for m in (dummy, gbdt_gamma, gbdt_ls, gbdt_pois, lgbm_gamma, lgbm_pois):
    message = f"training gamma deviance {m.__class__.__name__}"
    if hasattr(m, "loss"):
        message += " " + m.loss
    elif hasattr(m, "objective"):
        message += " " + m.objective
    print(f"{message: <68}: {mean_gamma_deviance(y_train, np.clip(m.predict(X_train), 1e-12, None))}")

for m in (dummy, gbdt_gamma, gbdt_ls, gbdt_pois, lgbm_gamma, lgbm_pois):
    message = f"test gamma deviance {m.__class__.__name__}"
    if hasattr(m, "loss"):
        message += " " + m.loss
    elif hasattr(m, "objective"):
        message += " " + m.objective
    print(f"{message: <68}: {mean_gamma_deviance(y_test, np.clip(m.predict(X_test), 1e-12, None))}")

I replaced Poisson with squared error as model to compare against, bb234ee. Now it works. Comments on comparison to Poisson are preserved.