Incorrect delta for huber objective

[thebigspin]

Description

huber objective implementation in lightgbm uses alpha directly as delta. However, alpha is meant to be interpreted as a percentile. This means the default alpha=0.9 is often meaningless and not targeted to the problem.

Instead, delta should be calculated from alpha, e.g. as (C++ equivalent of) delta=np.quantile(np.abs(labels - labels.mean()), 0.9)

LightGBM/src/objective/regression_objective.hpp

Lines 314 to 339 in 64ab66e

	void GetGradients(const double* score, score_t* gradients,
	score_t* hessians) const override {
	if (weights_ == nullptr) {
	#pragma omp parallel for num_threads(OMP_NUM_THREADS()) schedule(static)
	for (data_size_t i = 0; i < num_data_; ++i) {
	const double diff = score[i] - label_[i];
	if (std::abs(diff) <= alpha_) {
	gradients[i] = static_cast<score_t>(diff);
	} else {
	gradients[i] = static_cast<score_t>(Common::Sign(diff) * alpha_);
	}
	hessians[i] = 1.0f;
	}
	} else {
	#pragma omp parallel for num_threads(OMP_NUM_THREADS()) schedule(static)
	for (data_size_t i = 0; i < num_data_; ++i) {
	const double diff = score[i] - label_[i];
	if (std::abs(diff) <= alpha_) {
	gradients[i] = static_cast<score_t>(diff * weights_[i]);
	} else {
	gradients[i] = static_cast<score_t>(Common::Sign(diff) * static_cast<score_t>(weights_[i]) * alpha_);
	}
	hessians[i] = static_cast<score_t>(weights_[i]);
	}
	}
	}

[jameslamb]

Thanks for using LightGBM and for opening this issue. I investigated it tonight and I agree... this sure looks like a bug.

I've replaced your copied-in code with a direct link to the code as of a specific link (you'll see it still renders). That's preferable because it helps anyone finding this from search to understand which commit you were looking at. If you're unfamiliar with that, see https://docs.github.com/en/repositories/working-with-files/using-files/getting-permanent-links-to-files

---

uses alpha directly as delta

I investigated this today. Back when this was first introduced (#174), it was actually named huber_delta. @chivee suggested renaming it to alpha in #1043 (comment), to align with scikit-learn, and it's been that way here (and in scikit-learn) ever since.

alpha is meant to be interpreted as a percentile

scikit-learn appears to use this parameter in the same way LightGBM does. It must not be treated as a percentile, as the docs say:

alpha : float, default=0.0001
    Strength of the squared L2 regularization. Note that the penalty is
    equal to ``alpha * ||w||^2``.
    Must be in the range `[0, inf)`.

source code: https://github.com/scikit-learn/scikit-learn/blob/fe2edb3cdbd75ae4e662fda67dcb19277258792b/sklearn/linear_model/_huber.py#L159-L162

LightGBM's implementation (remembering that the code linked in the description is the first and second derivatives) looks consistent with other links I found.

From https://en.wikipedia.org/wiki/Huber_loss

And from https://docs.pytorch.org/docs/stable/generated/torch.nn.HuberLoss.html

_ This means the default alpha=0.9 is often meaningless and not targeted to the problem._

I do think this looks true :/

Looking at this:

LightGBM/src/objective/regression_objective.hpp

Lines 319 to 324 in 64ab66e

	const double diff = score[i] - label_[i];
	if (std::abs(diff) <= alpha_) {
	gradients[i] = static_cast<score_t>(diff);
	} else {
	gradients[i] = static_cast<score_t>(Common::Sign(diff) * alpha_);
	}

It does seem like for good performance, alpha as LightGBM uses it should be on the scale of the label.

Tested like this, using lightgbm built from latest master (eb27a2f)

from sklearn.datasets import make_regression
from sklearn.metrics import r2_score
import numpy as np
import lightgbm as lgb

X, y = make_regression(n_samples=10_000, n_informative=10, random_state=312)

params = {
    "alpha": 0.9,
    "deterministic": True,
    "force_row_wise": True,
    "objective": "huber",
    "num_iterations": 10,
    "num_threads": 1,
    "seed": 708,
}

# initial dataset
bst_scale0 = lgb.train(
    params=params,
    train_set=lgb.Dataset(X, label=y)
)

# label multiplied by 10
bst_scale10 = lgb.train(
    params=params,
    train_set=lgb.Dataset(X, label=y*10)
)

# initial dataset, alpha=90th percentile of deviations from the mean
bst_scale0_alpha_pct = lgb.train(
    params={
        **params,
        "alpha": np.quantile(np.abs(y - np.mean(y)), 0.9)
    },
    train_set=lgb.Dataset(X, label=y)
)

# label multiplied by 10, alpha=90th percentile of deviations from the mean
bst_scale10_alpha_pct = lgb.train(
    params={
        **params,
        "alpha": np.quantile(np.abs(10 * y - np.mean(10 * y)), 0.9)
    },
    train_set=lgb.Dataset(X, label=10 * y)
)

r2_score(y, bst_scale0.predict(X))
# 0.0049

r2_score(y * 10, bst_scale10.predict(X))
# 0.00049 <-- worse!

r2_score(y, bst_scale0_alpha_pct.predict(X))
# 0.6294

r2_score(y * 10, bst_scale10_alpha_pct.predict(X))
# 0.6294

Instead, delta should be calculated from alpha

At this point, changing the meaning of the parameter alpha for Huber loss would be a pretty significant breaking change. I don't think LightGBM should do that.

But we should do something to fix this situation.

Some backwards-compatible ideas I can think of:

• introduce a huber_method with values like "alpha_is_percentile" and "alpha_is_residual" (default "alpha_is_residual", same behavior as today)
• introduce a new huber_alpha which is bounded in (0, 1] and calculate huber_delta as that percentile of the label's deviations from its mean (default value -1 means to preserve the current use of alpha)
• introduce a warning that suggests a better value of alpha at the beginning of training
• introduce a totally new objective huber_scaled or huber_autoscaled and a configuration parameter huber_autoscaled_alpha (names TBD) which implements this behavior