How does XGBoost calculate base_score?

Since XGBoost 2.0 base_score is automatically calculated if it is not specified when initialising an estimator. I naively thought it would simply use the mean of the target, but this does not seem to be the case:

import json
import shap # only for the dataset
import xgboost as xgb

print('shap.__version__:',shap.__version__)
print('xgb.__version__:',xgb.__version__)
print()

X, y = shap.datasets.adult()

estimator = xgb.XGBClassifier(
    objective='binary:logistic',
    n_estimators=200)

estimator.fit(X,y)

print('y.mean():',y.mean())
print("float(json.loads(estimator.get_booster().save_config())['learner']['learner_model_param']['base_score']):",float(json.loads(estimator.get_booster().save_config())['learner']['learner_model_param']['base_score']))

Output:

shap.__version__: 0.46.0
xgb.__version__: 2.1.0

y.mean(): 0.2408095574460244
float(json.loads(estimator.get_booster().save_config())['learner']['learner_model_param']['base_score']): 0.26177529

The difference is too big for a rounding error. So how is base_score calculated? I think this is the relevant commit, but it's hard to tell exactly what it does.

TL;DR

Before XGBoost 2.0,

• For classification (binary:logistic), base_score was initialized to 0.5 (i.e., a probability; the raw score is logit(0.5) = 0).
• For regression (reg:squarederror), base_score was initialized to 0.

Since XGBoost 2.0, these initializations remain unchanged at first. However, before training the first tree, XGBoost internally fits a one-node tree (a stump) and sets the optimal output weight of that node as the updated base_score.

This logic can be found in the source code here.

Simple Derivation

Classification (binary:logistic)

Assume:

• n positive samples (label = 1)
• m negative samples (label = 0)

Initial prediction for every sample is:

• Predicted probability: pred_prob = base_score = 0.5

Then the gradients and Hessians for each sample are:

• Gradient:
  • Positive sample: g_1 = pred_prob - label = 0.5 - 1 = -0.5
  • Negative sample: g_0 = pred_prob - label = 0.5 - 0 = 0.5
• Hessian (same for all): h = pred_prob * (1 - pred_prob) = 0.5 * (1 - 0.5) = 0.25

Total gradient and Hessian:

• G = -0.5n + 0.5m = 0.5(m - n)
• H = 0.25(m + n)

Optimal constant prediction (i.e., the stump leaf value) is:

• w* = -G/H = - 2(m - n)(m + n)

This raw value is then passed through the sigmoid function:

• base_score_new = sigmoid(w*)

So the new base_score estimated by XGBoost is the probability corresponding to that optimal constant.

import numpy as np  
  
def cal_base_score(train_label):
    n = sum(train_label)  # num of 1s
    m = len(train_label) - n  # num of 0s
    w = -2 * (m - n) / (m + n)  
    base_score = 1 / (1 + np.exp(-w))  
    return base_score

Regression (reg:squarederror)

Initial prediction for all samples:

• pred = base_score = 0

Then for each sample:

• Gradient: g_i = pred_i - label_i = 0 - label_i = -label_i
• Hessian: h_i = 1

Total gradient and Hessian:

• G = sum(-label_i), where label_i is the label of each sample
• H = N, where N is the total number of samples

Optimal constant prediction:

• w* = -G/H = sum(label_i)/N = mean(label)

So the new base_score used by XGBoost is simply the mean of the labels: base_score_new = w* = mean of label