0x00 前言
本文是《GBDT源码分析》系列的第三篇，主要关注和GBDT本身以及Ensemble算法在scikit-learn中的实现。
0x01 整体说明
scikit-learn的ensemble模块里包含许多各式各样的集成模型，所有源码均在sklearn/ensemble文件夹里，代码的文件结构可以参考该系列的第一篇文章。其中 GradientBoostingRegressor 和 GradientBoostingClassifier分别是基于grandient_boosting的回归器和分类器。
0x02 源码结构分析
BaseEnsemble
ensemble中的所有模型均基于基类BaseEnsemble，该基类在sklearn/ensemble/base.py里。BaseEnsemble继承了两个父类，分别是BaseEstimator和MetaEstimatorMixin。BaseEnsemble里有如下几个方法，基本都是私有方法：
-
__init__: 初始化方法，共三个参数base_estimator, n_estimators, estimator_params。 -
_validate_estimator: 对n_estimators和base_estimator做检查，其中base_enstimator指集成模型的基模型。在GBDT中，base_estimator（元算法/基模型）是决策树。 -
_make_estimator: 从base_estimator中复制参数。 -
__len__: 返回ensemble中estimator的个数。 -
__getitem__: 返回ensemble中第i个estimator。 -
__iter__: ensemble中所有estimator的迭代器。
gradient_boosting.py
GradientBoostingClassifier和GradientBoostingRegressor这两个模型的实现均在gradient_boosting.py里。事实上，该脚本主要实现了一个集成回归树Gradient Boosted Regression Tree，而分类器和回归器都是基于该集成回归树的。
注意： 这里要先说明一个问题，GBDT本质上比较适合回归和二分类问题，而并不特别适用于多分类问题。在scikit-learn中，处理具有K类的多分类问题的GBDT算法实际上在每一次迭代里都构建了K个决策树，分别对应于这K个类别。我们在理论学习时并没有怎么接触过这个点，这里也并不对多分类问题做阐述。
实现Gradient Boosted Regression Tree的类是BaseGradientBoosting，它是GradientBoostingClassifier和GradientBoostingRegressor的基类。它实现了一个通用的fit方法，而回归问题和分类问题的区别主要在于损失函数LossFunction上。
下面我们梳理一下gradient_boosting.py的内容。
基本的estimator类
一些简单基本的estimator类，主要用于LossFunction中init_estimator的计算，即初始预测值的计算。举例来说：QuantileEstimator(alpha=0.5) 代表用中位数作为模型最初的预测值。
-
QuantileEstimator: 预测训练集target的alpha-百分位的estimator。 -
MeanEstimator: 预测训练集target的平均值的estimator。 -
LogOddsEstimator: 预测训练集target的对数几率的estimator（适合二分类问题）。 -
ScaledLogOddsEstimator: 缩放后的对数几率（适用于指数损失函数）。 -
PriorProbabilityEstimator: 预测训练集中每个类别的概率。 -
ZeroEstimator： 预测结果都是0的estimator。
LossFunction
LossFunction: 损失函数的基类，有以下一些主要方法
-
__init__: 输入为n_classes。当问题是回归或二分类问题时，n_classes为1；K类多分类为题时n_classes为K。 -
init_estimator: LossFunction的初始estimator，对应上述那些基本的estimator类，用来计算模型初始的预测值。在基类中不实现，并抛出NotImplementedError。 -
negative_gradient： 根据预测值和目标值计算负梯度。 -
update_terminal_regions: 更新树的叶子节点，更新模型的当前预测值。 -
_update_terminal_regions: 更新树的叶子节点的方法模板。
RegressionLossFunction
RegressionLossFunction: 继承LossFunction类，是所有回归损失函数的基类。
-
LeastSquaresError: init_estimator是MeanEstimator；负梯度是目标值y和预测值pred的差；唯一一个在update_terminal_regions中不需要更新叶子节点value的LossFunction。 -
LeastAbsoluteError: init_estimator是QuantileEstimator(alpha=0.5)；负梯度是目标值y和预测值pred的差的符号，适用于稳健回归。 -
HuberLossFunction: 一种适用于稳健回归Robust Regression的损失函数，init_estimator是QuantileEstimator(alpha=0.5)。 -
QuantileLossFunction: 分位数回归的损失函数，分位数回归允许估计目标值条件分布的百分位值。
ClassificationLossFunction
ClassificationLossFunction: 继承LossFunction，是所有分类损失函数的基类。有_score_to_proba（将分数转化为概率的方法模板）以及_score_to_decision（将分数转化为决定的方法模板）两个抽象方法。
-
BinomialDeviance: 二分类问题的损失函数，init_estimator是LogOddsEstimator。 -
MultinomialDeviance: 多分类问题的损失函数，init_estimator是PriorProbabilityEstimator。 -
ExponentialLoss: 二分类问题的指数损失，等同于AdaBoost的损失。init_estimator是ScaledLogOddsEstimator。
其它
gradient_boosting.py文件里还有以下几个内容，再下一章里我们将对这些内容做深入分析。
-
VerboseReporter: 输出设置。 -
BaseGradientBoosting: Gradient Boosted Regression Tree的实现基类。 -
GradientBoostingClassifier: Gradient Boosting分类器。 -
GradientBoostingRegressor: Gradient Boosting回归器。
0x03 GDBT主要源码分析
BaseGradientBoosting
下面我们具体看一下BaseGradientBoosting这个基类，该类有几个主要方法：
-
__init__: 除了来自决策树的参数外，还有n_estimators, learning_rate, loss, init, alpha, verbose, warm_start几个新的参数。其中loss是损失函数LossFunction的选择，learning_rate为学习率，n_estimators是boosting的次数（迭代次数或者Stage的个数）。通常learning_rate和n_estimators中需要做一个trade_off上的选择。init参数指的是BaseEstimator，即默认为每个LossFunction里面的init_estimator，用来计算初始预测值。warm_start决定是否重用之前的结果并加入更多estimators，或者直接抹除之前的结果。 -
_check_params: 检查参数是否合法，以及初始化模型参数（包括所用的loss等）。 -
_init_state: 初始化init_estimator以及model中的状态（包括init_estimator、estimators_、train_score_、oob_improvement_，后三个都是数组，分别存储每一个Stage对应的estimator、训练集得分和outofbag评估的进步）。 -
_clear_state: 清除model的状态。 -
_resize_state: 调整estimators的数量。 -
fit: 训练gradient boosting模型，会调用_fit_stages方法。 -
_fit_stages: 迭代boosting训练的过程，会调用_fit_stage方法。 -
_fit_stage: 单次stage训练过程。 -
_decision_function：略。 -
_staged_decision_function：略。 -
_apply：返回样本在每个estimator落入的叶子节点编号。
fit方法
# 如果不是warm_start，清除之前的model状态
if not self.warm_start:
____self._clear_state()
# ......检查输入、参数是否合法
# ......如果模型没有被初始化，则初始化模型，训练出初始模型以及预测值；
# ......如果模型已被初始化，判断n_estimators的大小并重新设置模型状态。
# boosting训练过程，调用_fit_stages方法
n_stages = self._fit_stages(X, y, y_pred, sample_weight, random_state, begin_at_stage, monitor, X_idx_sorted)
# ......当boosting训练次数与初始化的estimators_长度不一致时，修正相关变量/状态，包括estimators_、train_score_、oob_improvement_
# fit方法返回self
_fit_stages方法
# 获取样本数（为什么每次都要获取样本数而不作为self.n_samples呢）
n_samples = X.shape[0]
# 判断是否做oob（交叉验证），仅当有抽样时才做oob
do_oob = self.subsample < 1.0
# 初始化sample_mask，即标注某一轮迭代每个样本是否要被抽样的数组
sample_mask = np.ones((n_samples, ), dtype=np.bool)
# 计算inbag（用来训练的样本个数）
n_inbag = max(1, int(self.subsample * n_samples))
# 获取loss对象
loss_ = self.loss_
# ......设置min_weight_leaf、verbose、sparsity等相关参数
# 开始boosting迭代
# begin_at_stage是迭代初始次数，一般来说是0，如果是warm_start则是之前模型结束的地方
i = begin_at_stage
# 开始迭代
for i in range(begin_at_stage, self.n_estimators):
# 如果subsample < 1，do_oob为真，做下采样
if do_oob:
# _random_sample_mask是在_gradient_boosting.pyx里用cpython实现的一个方法，用来做随机采样，生成inbag/outofbag样本（inbag样本为True）
sample_mask = _random_sample_mask(n_samples, n_inbag, random_state)
# 获得之前的oob得分
old_oob_score = loss_(y[~sample_mask], y_pred[~sample_mask], sample_weight[~sample_mask])
# 调用_fit_stage来训练下一阶段的数
y_pred = self._fit_stage(i, X, y, y_pred, sample_weight, sample_mask, random_state, X_idx_sorted, X_csc, X_csr)
# 跟踪偏差/loss
# 当do_oob时，计算训练样本的loss和oob_score的提升值
if do_oob:
# inbag训练样本的loss
self.train_score_[i] = loss_(y[sample_mask], y_pred[sample_mask], sample_weight[sample_mask])
# outofbag样本的loss提升
self.oob_improvement_[i] = (old_oob_score - loss_(y[~sample_mask], y_pred[~sample_mask], sample_weight[~sample_mask]))
# subsample为1时
else:
self.train_score_[i] = loss_(y, y_pred, sample_weight)
# 若verbose大于0，更新标准输出
if self.verbose > 0:
verbose_reporter.update(i, self)
# ......若有monitor，检查是否需要early_stopping
# _fit_stages方法返回i+1，即迭代总次数（包括warm_start以前的迭代）
_fit_stage方法
# 判断sample_mask的数据类型
assert sample_mask.dtype == np.bool
# 获取损失函数
loss = self.loss_
# 获取目标值
original_y = y
# 这里K针对的是多分类问题，回归和二分类时K为1
for k in range(loss.K):
# 当问题是多分类问题时，获取针对该分类的y值
if loss.is_multi_class:
y = np.array(original_y == k, dtype=np.float64)
# 计算当前负梯度
residual = loss.negative_gradient(y, y_pred, k=k, sample_weight=sample_weight)
# 构造决策回归树（事实上是对负梯度做决策树模型）
tree = DecisionTreeRegressor(
criterion=self.criterion,
splitter='best',
max_depth=self.max_depth,
min_samples_split=self.min_samples_split,
min_samples_leaf=self.min_samples_leaf,
min_weight_fraction_leaf=self.min_weight_fraction_leaf,
min_impurity_split=self.min_impurity_split,
max_features=self.max_features,
max_leaf_nodes=self.max_leaf_nodes,
random_state=random_state,
presort=self.presort)
# 如果做sabsample，重新计算sample_weight
if self.subsample < 1.0:
sample_weight = sample_weight * sample_mask.astype(np.float64)
# 根据输入X是否稀疏，采用不同的fit方法，针对负梯度训练决策树
if X_csc is not None:
tree.fit(X_csc, residual, sample_weight=sample_weight, check_input=False, X_idx_sorted=X_idx_sorted)
else:
tree.fit(X, residual, sample_weight=sample_weight, check_input=False, X_idx_sorted=X_idx_sorted)
# 根据输入X是否稀疏，使用update_terminal_regions方法更新叶子节点（注意这是LossFunction里的一个方法）
if X_csr is not None:
loss.update_terminal_regions(tree.tree_, X_csr, y, residual, y_pred, sample_weight, sample_mask, self.learning_rate, k=k)
else:
loss.update_terminal_regions(tree.tree_, X, y, residual, y_pred, sample_weight, sample_mask, self.learning_rate, k=k)
# 将新的树加入到ensemble模型中
self.estimators_[i, k] = tree
# _fit_stage方法返回新的预测值y_pred，注意这里y_pred是在loss.update_terminal_regions计算的
LossFunction中的update_terminal_regions方法
为了加深理解，我么我们再看一下update_terminal_regions都做了什么。
# 计算每个样本对应到树的哪一个叶子节点
terminal_regions = tree.apply(X)
# 将outofbag的样本的结果都置为-1（不参与训练过程）
masked_terminal_regions = terminal_regions.copy()
masked_terminal_regions[~sample_mask] = -1
# 更新每个叶子节点上的value，tree.children_left == TREE_LEAF是判断叶子节点的方法。一个很关键的点是这里只更新了叶子节点，而只有LossFunction是LeastSquaresError时训练时生成的决策树上的value和我们实际上想要的某个节点的预测值是一致的。
for leaf in np.where(tree.children_left == TREE_LEAF)[0]:
# _update_terminal_region由每个具体的损失函数具体实现，在LossFunction基类中只提供模板
self._update_terminal_region(tree, masked_terminal_regions, leaf, X, y, residual, y_pred[:, k], sample_weight)
# 更新预测值，tree预测的是负梯度值，预测值通过加上学习率 * 负梯度来更新，这里更新所有inbag和outofbag的预测值
y_pred[:, k] += (learning_rate * tree.value[:, 0, 0].take(terminal_regions, axis=0))
笔者之前使用GBDT做回归模型时观察每颗树的可视化结果，发现对于损失函数是ls（LeastSquaresError）的情况，每棵树的任意一个节点上的value都是当前点的target预估值差（即residual，所有树叶子节点预测的都是residual，它们的和是最终的预测结果）；但使用lad损失函数时，只有叶子节点的结果是收入预估值差。原因应该就在这里：
- ls对应的LossFunction类是LeastSquaresError，每个节点的value就是当前点的target预估值差，叶子节点也不需要更新。这是因为ls的负梯度计算方法是预测值和目标值的差，这本身就是residual的概念，所以所有节点的value都是我们想要的值。
- lad对应的LossFunction类是LeastAbsoluteError，每个节点的value并不是当前点的target预估值差，而最后代码里也只更新了叶子节点，所以可视化时会有一些问题，也不能直接获得每个节点的value作为target预估值差。事实上，lad在训练的时候有点像一个“二分类”问题，它的负梯度只有两种取值-1和1，即预测值比目标值大还是小，然后根据这个标准进行分裂。
所以如果没有改源码并重新训练模型的话，若不是ls，其它已有的GBDT模型没有办法直接获取每个非叶子结点的target预估值差，这个在分析模型时会有一些不方便的地方。
feature_importances_的计算
# 初始化n_feautures_长度的数组
total_sum = np.zeros((self.n_features_, ), dtype=np.float64)
# 对于boosting模型中的每一个estimator（实际上就是一棵树，多分类是多棵树的数组）
for stage in self.estimators_:
# 当前stage每个feature在各个树内的所有的importance平均（多分类时一个stage有多棵树）
stage_sum = sum(tree.feature_importances_ for tree in stage) / len(stage)
# 累加各个stage的importance
total_sum += stage_sum
# 做归一化
importances = total_sum / len(self.estimators_)
GradientBoostingClassifier
GBDT分类器的loss可取deviance或exponential，分别对应MultinomialDeviance和ExponentialLoss这两个损失函数。分类器在predict时需要多加一步，把不同类别对应的树的打分综合起来才能输出结果。因此GBDT实际上不太适合做多分类问题。
GradientBoostingRegressor
GBDT回归器的loss可取ls, lad, huber, quantile，分别对应LeastSquaresError, LeastAbsoluteError, HuberLossFunction, QuantileLossFunction这几个损失函数。
0xFF 总结
至此，该系列三篇文章已结。
个人主页：http://cathyxlyl.github.io/
文章可以转载, 但必须以超链接形式标明文章原始出处和作者信息
