LightGBM with Tweedie loss; I'm confused on the Gradient and Hessians used

Question

I'm trying to figure out custom objective functions in LightGBM, and I figured a good place to start would be replicating the built-in functions. The equation LightGBM uses to calculate the Tweedie metric (https://github.com/microsoft/LightGBM/blob/1c27a15e42f0076492fcc966b9dbcf9da6042823/src/metric/regression_metric.hpp#L300-L318) seems to match definitions of the Tweedie loss I've found online (https://towardsdatascience.com/tweedie-loss-function-for-right-skewed-data-2c5ca470678f), though they do a weird exp(ln(score)) process, I'm guessing for numerical stability. However, their equations for the gradient and Hessian seem to be done on the log of score directly (https://github.com/microsoft/LightGBM/blob/1c27a15e42f0076492fcc966b9dbcf9da6042823/src/objective/regression_objective.hpp#L702-L732).

It seems like they are using the equation:

gradients[i] = -label_[i] * e^((1 - rho_) * score[i]) + e^((2 - rho_) * score[i]);

where I would expect the gradient to be:

gradients[i] = -label_[i] * score[i]^(- rho_) + score[i]^(1 - rho_);

My guess is somewhere LightGBM is processing score as ln(score), like using parameter reg_sqrt, but I can't find where in the documentation this is described.

Anyway I've tried recreating both their formula and my own calculations as custom objective functions, and neither seem to work:

library(lightgbm)
library(data.table)

# Tweedie gradient with variance = 1.5, according to my own math
CustomObj_t1 <- function(preds, dtrain) {
  labels <- dtrain$getinfo('label')
  grad <- -labels * preds^(-3/2) + preds^(-1/2)
  hess <- 1/2 * (3*labels*preds^(-5/2) - preds^(-3/2))
  return(list(grad = grad, hess = hess))
}

# Tweedie gradient with variance = 1.5, recreating code from LightGBM github
CustomObj_t2 <- function(preds, dtrain) {
  labels <- dtrain$getinfo('label')
  grad <- -labels*exp(-1/2*preds) + exp(1/2*preds)
  hess <- -labels*(-1/2)*exp(-1/2*preds) + 1/2*exp(1/2*preds)
  return(list(grad = grad, hess = hess))
}

params = list(objective = "tweedie",
              seed = 1,
              metric = "rmse")

params2 = list(objective = CustomObj_t1,
               seed= 1,
               metric = "rmse")

params3 = list(objective = CustomObj_t2,
               seed= 1,
               metric = "rmse")

# Create data
set.seed(321)
db_Custom = data.table(a=runif(2000), b=runif(2000))
db_Custom[,X := (a*4+exp(b))]

# break into test and training sets
db_Test = db_Custom[1:10]
db_Custom=db_Custom[11:nrow(db_Custom),]

FeatureCols = c("a","b")

# Create dataset
ds_Custom <- lgb.Dataset(data.matrix(db_Custom[, FeatureCols, with = FALSE]), label = db_Custom[["X"]])  

# Train
fit = lgb.train(params, ds_Custom, verb=-1)
#print("  ")
fit2 = lgb.train(params2, ds_Custom, verb=-1)
#print("  ")
fit3 = lgb.train(params3, ds_Custom, verb=-1)

# Predict
pred = predict(fit, data.matrix(db_Test[, FeatureCols, with = FALSE]))
db_Test[, prediction := pmax(0, pred)]

pred2 = predict(fit2, data.matrix(db_Test[, FeatureCols, with = FALSE]))
db_Test[, prediction2 := pmax(0, pred2)]

pred3 = predict(fit3, data.matrix(db_Test[, FeatureCols, with = FALSE]))
db_Test[, prediction3 := pmax(0, pred3)]

print(db_Test[,.(X,prediction,prediction2,prediction3)])

I get the results (would expect prediction2 or prediction3 to be very similar to prediction):

"X"                "prediction"   "prediction2" "prediction3"
4.8931646234958     4.89996556839721    0   1.59154656425556
6.07328897031702    6.12313647937047    0   1.81022588429474
2.05728566704078    2.06824004875244    0   0.740577102751491
2.54732526765174    2.50329903656292    0   0.932517774958986
4.07044099941395    4.07047912554207    0   1.39922723582939
2.74639568121359    2.74408567443232    0   1.01628212910587
3.47720295158928    3.49241414141969    0   1.23049599462599
2.92043718858535    2.90464303454649    0   1.0680618051659
4.44415913080697    4.43091665909845    0   1.48607456777287
4.96566318066753    4.97898586895233    0   1.60163901781479

Is there something I'm missing? Am I just doing the math or coding wrong?

Answer 1

It appears, per the linked git page, and your prediction3 column, that if you exponentiate this column, it becomes very close to columns 0 and 1.