简单易学的机器学习算法——梯度提升决策树GBDT

梯度提升决策树（Gradient Boosting Decision Tree，GBDT）算法是近年来被提及比较多的一个算法，这主要得益于其算法的性能，以及该算法在各类数据挖掘以及机器学习比赛中的卓越表现，有很多人对GBDT算法进行了开源代码的开发，比较火的是陈天奇的XGBoost和微软的LightGBM。

1. 监督学习

1.1. 监督学习的主要任务

监督学习是机器学习算法中重要的一种，对于监督学习，假设有$m$个训练样本：

$$\left\{\left(X^{\left(1\right)},y^{\left(1\right)}\right),\left(X^{\left(2\right)},y^{\left(2\right)}\right),\cdots ,\left(X^{\left(m\right)},y^{\left(m\right)}\right)\right\}$$

其中，$X^{\left(i\right)}=\left\{x_1^{\left(i\right)},x_2^{\left(i\right)},\cdots ,x_m^{\left(i\right)}\right\}$称为第$i$个样本的特征，$y^{\left(i\right)}$称为第$i$个样本的标签，样本标签可以为离散值，如分类问题；也可以为连续值，如回归问题。在监督学习中，利用训练样本训练出模型，该模型能够实现从样本特征$X^{\left(i\right)}$到样本标签$y^{\left(i\right)}$的映射，即：

$$X^{\left(i\right)}\stackrel{F}{\to }{y}^{\left(i\right)}$$

为了能够对映射$F\left(X\right)$进行求解，通常对模型设置损失函数$L\left(y,F\left(X\right)\right)$，并求得在损失函数最小的情况下的映射为最好的映射：

$$F^{\ast }=\underset{F\left(X\right)}{\mathrm{argmin}}L\left(y,F\left(X\right)\right)$$

对于一个具体的问题，如线性回归问题，其映射函数的形式为：

$$F\left(X;W\right)=WX=w_0+w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n$$

此时对于最优映射函数$F\left(X;W\right)$的求解，实质是对映射函数中的参数$W$的求解。对于参数的求解方法有很多，如梯度下降法。

1.2. 梯度下降法

梯度下降法（Gradient Descent，GD）算法是求解最优化问题最简单、最直接的方法。梯度下降法是一种迭代的优化算法，对于优化问题：

$$minf\left(w\right)$$

其基本步骤为：

• 随机选择一个初始点$w_0$
• 重复以下过程：
  • 决定下降的方向：$d_i=-\frac{\partial }{\partial w}f\left(w\right){\mid }_{w_i}$
  • 选择步长$\rho$
  • 更新：$w_{i+1}=w_i+\rho \cdot d_i$
• 直到满足终止条件

梯度下降法的具体过程如下图所示：

由以上的过程，我们可以看出，对于最终的最优解$w^{\ast }$，是由初始值$w_0$经过$M$代的迭代之后得到的，在这里，设$w_0=d_0$，则$w^{\ast }$为：

$$w^{\ast }=\sum _{i=0}^M{\rho }_i\cdot d_i$$

1.3. 在函数空间的优化

以上是在指定的函数空间中对最优函数进行搜索，那么，能否直接在函数空间（function space）中查找到最优的函数呢？根据上述的梯度下降法的思路，对于模型的损失函数$L\left(y,F\left(X\right)\right)$，为了能够求解出最优的函数$F^{\ast }\left(X\right)$，首先，设置初始值为：

$$F_0\left(X\right)=f_0\left(X\right)$$

以函数$F\left(X\right)$作为一个整体，对于每一个样本$X^{\left(i\right)}$，都存在对应的函数值$F\left(X^{\left(i\right)}\right)$。与梯度下降法的更新过程一致，假设经过$M$代，得到最有的函数$F^{\ast }\left(X\right)$为：

$$F^{\ast }\left(X\right)=\sum _{i=0}^M{f}_i\left(X\right)$$

其中，$f_i\left(X\right)$为：

$$f_i\left(X\right)=-{\rho }_i{g}_m\left(X\right)$$

其中，$g_m\left(X\right)={\left[\frac{\partial L\left(y,F\left(X\right)\right)}{\partial F\left(X\right)}\right]}_{F\left(X\right)=F_{m-1}\left(X\right)}$。

由上述的过程可以得到函数$F\left(X\right)$的更新过程：

$$F_m\left(X\right)=\sum _{i=0}^m{f}_i\left(X\right)$$

与上面类似，函数$f\left(X\right)$是由参数$\text{a}$决定的，即：

$$f\left(X\right)=-\rho \cdot h\left(X;\text{a}\right)$$

2. Boosting

2.1. 集成方法之Boosting

Boosting方法是集成学习中重要的一种方法，在集成学习方法中最主要的两种方法为Bagging和Boosting，在Bagging中，通过对训练样本重新采样的方法得到不同的训练样本集，在这些新的训练样本集上分别训练学习器，最终合并每一个学习器的结果，作为最终的学习结果，Bagging方法的具体过程如下图所示：

在Bagging方法中，最重要的算法为随机森林Random Forest算法。由以上的图中可以看出，在Bagging方法中，$b$个学习器之间彼此是相互独立的，这样的特点使得Bagging方法更容易并行。与Bagging方法不同，在Boosting算法中，学习器之间是存在先后顺序的，同时，每一个样本是有权重的，初始时，每一个样本的权重是相等的。首先，第$1$个学习器对训练样本进行学习，当学习完成后，增大错误样本的权重，同时减小正确样本的权重，再利用第$2$个学习器对其进行学习，依次进行下去，最终得到$b$个学习器，最终，合并这$b$个学习器的结果，同时，与Bagging中不同的是，每一个学习器的权重也是不一样的。Boosting方法的具体过程如下图所示：

在Boosting方法中，最重要的方法包括：AdaBoost和GBDT。

2.2. Gradient Boosting

由上图所示的Boosting方法中，最终的预测结果为$b$个学习器结果的合并：

$$f\left(X\right)=\sum _{j=1}^b{\theta }_j{\phi }_j\left(X\right)$$

这与上述的在函数空间中的优化类似：

$$F_m\left(X\right)=\sum _{i=0}^m-{\rho }_i\cdot h\left(X;{\text{a}}_i\right)$$

根据如上的函数空间中的优化可知，每次对每一个样本的训练的值为：

$${\bar{y}}_i={\left[\frac{\partial L\left(y_i,F\left(X^{\left(i\right)}\right)\right)}{\partial F\left(X^{\left(i\right)}\right)}\right]}_{F\left(X\right)=F_{m-1}\left(X\right)}$$

上建立模型，由于上述是一个求解梯度的过程，因此也称为基于梯度的Boost方法，其具体过程如下所示：

3. Gradient Boosting Decision Tree

在上面简单介绍了Gradient Boost框架，梯度提升决策树Gradient Boosting Decision Tree是Gradient Boost框架下使用较多的一种模型，在梯度提升决策树中，其基学习器是分类回归树CART，使用的是CART树中的回归树。

3.1. 分类回归树CART

分类回归树CART算法是一种基于二叉树的机器学习算法，其既能处理回归问题，又能处理分类为题，在梯度提升决策树GBDT算法中，使用到的是CART回归树算法，对于CART树算法的更多信息，可以参考简单易学的机器学习算法——分类回归树CART。

对于一个包含了$m$个训练样本的回归问题，其训练样本为：

$$\left\{\left(X^{\left(1\right)},y^{\left(1\right)}\right),\left(X^{\left(2\right)},y^{\left(2\right)}\right),\cdots ,\left(X^{\left(m\right)},y^{\left(m\right)}\right)\right\}$$

其中，$X^{\left(i\right)}$为$n$维向量，表示的是第$i$个样本的特征，$y^{\left(i\right)}$为样本的标签，在回归问题中，标签$y^{\left(i\right)}$为一系列连续的值。此时，利用训练样本训练一棵CART回归树：

• 开始时，CART树中只包含了根结点，所有样本都被划分在根结点上：

此时，计算该节点上的样本的方差（此处要乘以$m$），方差表示的是数据的波动程度。那么，根节点的方差的$m$倍为：

$$s^2\cdot m={\left(y^{\left(1\right)}-\bar{y}\right)}^2+{\left(y^{\left(2\right)}-\bar{y}\right)}^2+\cdots +{\left(y^{\left(m\right)}-\bar{y}\right)}^2$$

其中，$\bar{y}$为标签的均值。此时，从$n$维特征中选择第$j$维特征，从$m$个样本中选择一个样本的值：$x_j$作为划分的标准，当样本$i$的第$j$维特征小于等于$x_j$时，将样本划分到左子树中，否则，划分到右子树中，通过以上的操作，划分到左子树中的样本个数为$m_1$，划分到右子树的样本的个数为$m_2=m-m_1$，其划分的结果如下图所示：

那么，什么样本的划分才是当前的最好划分呢？此时计算左右子树的方差之和：$s_1^2\cdot m_1+s_2^2\cdot m_2$：

$$s_1^2\cdot m_1+s_2^2\cdot m_2=\sum _{X^{\left(i\right)}\in left}{\left(y^{\left(i\right)}-{\bar{y}}_1\right)}^2+\sum _{X^{\left(j\right)}\in right}{\left(y^{\left(j\right)}-{\bar{y}}_2\right)}^2$$

其中，${\bar{y}}_1$为左子树中节点标签的均值，同理，${\bar{y}}_2$为右子树中节点标签的均值。选择其中$s_1^2\cdot m_1+s_2^2\cdot m_2$最小的划分作为最终的划分，依次这样划分下去，直到得到最终的划分，划分的结果为：

注意：对于上述最优划分标准的选择，以上的计算过程可以进一步优化。

首先，对于$s^2\cdot m$：

$$\begin{array}{rl}{s}^2\cdot m& =\sum _{X^{\left(i\right)}}{\left(y^{\left(i\right)}-\bar{y}\right)}^2\\ & =\sum _{X^{\left(i\right)}}\left({\left(y^{\left(i\right)}\right)}^2-2y^{\left(i\right)}\cdot \bar{y}+{\left(\bar{y}\right)}^2\right)\\ & =\sum _{X^{\left(i\right)}}{\left(y^{\left(i\right)}\right)}^2-\frac{2}{m}{\left(\sum _{X^{\left(i\right)}}{y}^{\left(i\right)}\right)}^2+\frac{1}{m}{\left(\sum _{X^{\left(i\right)}}{y}^{\left(i\right)}\right)}^2\\ & =\sum _{X^{\left(i\right)}}{\left(y^{\left(i\right)}\right)}^2-\frac{1}{m}{\left(\sum _{X^{\left(i\right)}}{y}^{\left(i\right)}\right)}^2\end{array}$$

而对于$s_1^2\cdot m_1+s_2^2\cdot m_2$：

$$\begin{array}{rl}{s}_1^2\cdot m_1+s_2^2\cdot m_2& =\sum _{X^{\left(i\right)}\in left}{\left(y^{\left(i\right)}-{\bar{y}}_1\right)}^2+\sum _{X^{\left(j\right)}\in right}{\left(y^{\left(j\right)}-{\bar{y}}_2\right)}^2\\ & =\sum _{X^{\left(i\right)}\in left}\left({\left(y^{\left(i\right)}\right)}^2-2y^{\left(i\right)}\cdot {\bar{y}}_1+{\left({\bar{y}}_1\right)}^2\right)+\sum _{X^{\left(j\right)}\in right}\left({\left(y^{\left(j\right)}\right)}^2-2y^{\left(j\right)}\cdot {\bar{y}}_2+{\left({\bar{y}}_2\right)}^2\right)\end{array}$$

$$\begin{array}{rl}& =\sum _{X^{\left(i\right)}}{\left(y^{\left(i\right)}\right)}^2-\frac{2}{m_1}{\left(\sum _{X^{\left(i\right)}\in left}{y}^{\left(i\right)}\right)}^2+\frac{1}{m_1}{\left(\sum _{X^{\left(i\right)}\in left}{y}^{\left(i\right)}\right)}^2\\ & -\frac{2}{m_2}{\left(\sum _{X^{\left(j\right)}\in right}{y}^{\left(j\right)}\right)}^2+\frac{1}{m_2}{\left(\sum _{X^{\left(j\right)}\in right}{y}^{\left(j\right)}\right)}^2\\ & =\sum _{X^{\left(i\right)}}{\left(y^{\left(i\right)}\right)}^2-\frac{1}{m_1}{\left(\sum _{X^{\left(i\right)}\in left}{y}^{\left(i\right)}\right)}^2-\frac{1}{m_2}{\left(\sum _{X^{\left(j\right)}\in right}{y}^{\left(j\right)}\right)}^2\end{array}$$

通过以上的过程，我们发现，划分前，记录节点的值为：

$$\frac{1}{m}{\left(\sum _{X^{\left(i\right)}}{y}^{\left(i\right)}\right)}^2$$

当划分后，两个节点的值的和为：

$$\frac{1}{m_1}{\left(\sum _{X^{\left(i\right)}\in left}{y}^{\left(i\right)}\right)}^2+\frac{1}{m_2}{\left(\sum _{X^{\left(j\right)}\in right}{y}^{\left(j\right)}\right)}^2$$

最好的划分，对应着两个节点的值的和的最大值。

3.2. GBDT——二分类

在梯度提升决策树GBDT中，通过定义不同的损失函数，可以完成不同的学习任务，二分类是机器学习中一类比较重要的分类算法，在二分类中，其损失函数为：

L ( y , F ) = l o g ( 1 + e x p ( − 2 y F ) ) ,    y ∈ { − 1 , 1 } L\left ( y,F \right )=log\left ( 1+exp\left ( -2yF \right ) \right ),\; y\in\left \{ -1,1 \right \} L(y,F)=log(1+exp(−2yF)),y∈{−1,1}

套用上面介绍的GB框架，得到下述的二分类GBDT的算法：

在构建每一棵CART回归树的过程中，对一个样本的预测值应与$\tilde{y}$尽可能一致，对于$\tilde{y}$，其计算过程为：

$$\begin{array}{rl}{\tilde{y}}^{\left(i\right)}& =-{\left[\frac{\partial L\left(y^{\left(i\right)},F\left(X^{\left(i\right)}\right)\right)}{\partial F\left(X^{\left(i\right)}\right)}\right]}_{F\left(X\right)=F_{m-1}\left(X\right)}\\ & =-{\left[\frac{\partial log\left(1+exp\left(-2y^{\left(i\right)}F\left(X^{\left(i\right)}\right)\right)\right)}{\partial F\left(X^{\left(i\right)}\right)}\right]}_{F\left(X\right)=F_{m-1}\left(X\right)}\\ & =-{\left[\frac{1}{1+exp\left(-2y^{\left(i\right)}F\left(X^{\left(i\right)}\right)\right)}\cdot exp\left(-2y^{\left(i\right)}F\left(X^{\left(i\right)}\right)\right)\cdot \left(-2y^{\left(i\right)}\right)\right]}_{F\left(X\right)=F_{m-1}\left(X\right)}\end{array}$$

$$\begin{array}{rl}& ={\frac{2y^{\left(i\right)}\cdot exp\left(-2y^{\left(i\right)}F\left(X^{\left(i\right)}\right)\right)}{1+exp\left(-2y^{\left(i\right)}F\left(X^{\left(i\right)}\right)\right)}}_{F\left(X\right)=F_{m-1}\left(X\right)}\\ & =\frac{2y^{\left(i\right)}}{1+exp\left(2y^{\left(i\right)}{F}_{m-1}\left(X^{\left(i\right)}\right)\right)}\end{array}$$

在$\tilde{y}$（通常有的地方称为残差，在这里，更准确的讲是梯度下降的方向）上构建CART回归树。最终将每一个训练样本划分到对应的叶子节点中，计算此时该叶子节点的预测值：

$${\gamma }_{jm}=\underset{\gamma }{\mathrm{argmin}}\sum _{X^{\left(i\right)}\in R_{jm}}log\left(1+exp\left(-2y^{\left(i\right)}\left(F_{m-1}\left(X^{\left(i\right)}\right)+\gamma \right)\right)\right)$$

由Newton-Raphson迭代公式可得：

$${\gamma }_{jm}=\frac{{\sum }_{X^{\left(i\right)}\in R_{jm}}{\tilde{y}}^{\left(i\right)}}{{\sum }_{X^{\left(i\right)}\in R_{jm}}\left|{\tilde{y}}^{\left(i\right)}\right|\left(2-\left|{\tilde{y}}^{\left(i\right)}\right|\right)}$$

以参考文献3 Idiots’ Approach for Display Advertising Challenge中提供的代码为例：

• GBDT训练的主要代码为：

void GBDT::fit(Problem const &Tr, Problem const &Va)
{
        bias = calc_bias(Tr.Y); //用于初始化的F

        std::vector<float> F_Tr(Tr.nr_instance, bias), F_Va(Va.nr_instance, bias);

        Timer timer;
        printf("iter     time    tr_loss    va_loss\n");
        // 开始训练每一棵CART树
        for(uint32_t t = 0; t < trees.size(); ++t)
        {
                timer.tic();

                std::vector<float> const &Y = Tr.Y;
                std::vector<float> R(Tr.nr_instance), F1(Tr.nr_instance); // 记录残差和F

                #pragma omp parallel for schedule(static)
                for(uint32_t i = 0; i < Tr.nr_instance; ++i)
                        R[i] = static_cast<float>(Y[i]/(1+exp(Y[i]*F_Tr[i]))); //计算残差，或者称为梯度下降的方向

                // 利用上面的残差值，在此函数中构造一棵树
                trees[t].fit(Tr, R, F1); // 分类树的生成

                double Tr_loss = 0;
                // 用上面训练的结果更新F_Tr，并计算log_loss
                #pragma omp parallel for schedule(static) reduction(+: Tr_loss)
                for(uint32_t i = 0; i < Tr.nr_instance; ++i)
                {
                        F_Tr[i] += F1[i];
                        Tr_loss += log(1+exp(-Y[i]*F_Tr[i]));
                }
                Tr_loss /= static_cast<double>(Tr.nr_instance);

                // 用上面训练的结果预测测试集，打印log_loss
                #pragma omp parallel for schedule(static)
                for(uint32_t i = 0; i < Va.nr_instance; ++i)
                {
                        std::vector<float> x = construct_instance(Va, i);
                        F_Va[i] += trees[t].predict(x.data()).second;
                }

                double Va_loss = 0;
                #pragma omp parallel for schedule(static) reduction(+: Va_loss)
                for(uint32_t i = 0; i < Va.nr_instance; ++i)
                        Va_loss += log(1+exp(-Va.Y[i]*F_Va[i]));
                Va_loss /= static_cast<double>(Va.nr_instance);

                printf("%4d %8.1f %10.5f %10.5f\n", t, timer.toc(), Tr_loss, Va_loss);
                fflush(stdout);
        }
}

• CART回归树的训练代码为：

void CART::fit(Problem const &prob, std::vector<float> const &R, std::vector<float> &F1){
	uint32_t const nr_field = prob.nr_field; // 特征的个数
	uint32_t const nr_sparse_field = prob.nr_sparse_field;
	uint32_t const nr_instance = prob.nr_instance; // 样本的个数

	std::vector<Location> locations(nr_instance); // 样本信息

	#pragma omp parallel for schedule(static)
	for(uint32_t i = 0; i < nr_instance; ++i)
		locations[i].r = R[i]; // 记录每一个样本的残差

	for(uint32_t d = 0, offset = 1; d < max_depth; ++d, offset *= 2){// d:深度

		uint32_t const nr_leaf = static_cast<uint32_t>(pow(2, d)); // 叶子节点的个数


		std::vector<Meta> metas0(nr_leaf); // 叶子节点的信息

		for(uint32_t i = 0; i < nr_instance; ++i){

			Location &location = locations[i]; //第i个样本的信息
 
			if(location.shrinked)
				continue;

			Meta &meta = metas0[location.tnode_idx-offset]; //找到对应的叶子节点

			meta.s += location.r; //残差之和
			++meta.n;
		}

		std::vector<Defender> defenders(nr_leaf*nr_field); //记录每一个叶节点的每一维特征
		std::vector<Defender> defenders_sparse(nr_leaf*nr_sparse_field);
		// 针对每一个叶节点

		for(uint32_t f = 0; f < nr_leaf; ++f){

			Meta const &meta = metas0[f]; // 叶子节点

			double const ese = meta.s*meta.s/static_cast<double>(meta.n); //该叶子节点的ese

			for(uint32_t j = 0; j < nr_field; ++j)
				defenders[f*nr_field+j].ese = ese;

			for(uint32_t j = 0; j < nr_sparse_field; ++j)
				defenders_sparse[f*nr_sparse_field+j].ese = ese;
		}
		
		std::vector<Defender> defenders_inv = defenders;

		std::thread thread_f(scan, std::ref(prob), std::ref(locations),
				std::ref(metas0), std::ref(defenders), offset, true);
		std::thread thread_b(scan, std::ref(prob), std::ref(locations),
				std::ref(metas0), std::ref(defenders_inv), offset, false);
		scan_sparse(prob, locations, metas0, defenders_sparse, offset, true);
		thread_f.join();
		thread_b.join();

		// 找出最佳的ese，scan里是每个字段的最佳ese，这里是所有字段的最佳ese，赋值给相应的tnode
		for(uint32_t f = 0; f < nr_leaf; ++f){
			// 对于每一个叶节点都找到最好的划分
			Meta const &meta = metas0[f];
			double best_ese = meta.s*meta.s/static_cast<double>(meta.n);

			TreeNode &tnode = tnodes[f+offset];
			for(uint32_t j = 0; j < nr_field; ++j){

				Defender defender = defenders[f*nr_field+j];//每一个叶节点都对应着所有的特征

				if(defender.ese > best_ese)
				{
					best_ese = defender.ese;
					tnode.feature = j;
					tnode.threshold = defender.threshold;
				}

				defender = defenders_inv[f*nr_field+j];
				if(defender.ese > best_ese)
				{
					best_ese = defender.ese;
					tnode.feature = j;
					tnode.threshold = defender.threshold;
				}
			}
			for(uint32_t j = 0; j < nr_sparse_field; ++j)
			{
				Defender defender = defenders_sparse[f*nr_sparse_field+j];
				if(defender.ese > best_ese)
				{
					best_ese = defender.ese;
					tnode.feature = nr_field + j;
					tnode.threshold = defender.threshold;
				}
			}
		}

		// 把每个instance都分配给树里的一个叶节点下
		#pragma omp parallel for schedule(static)
		for(uint32_t i = 0; i < nr_instance; ++i){

			Location &location = locations[i];
			if(location.shrinked)
				continue;

			uint32_t &tnode_idx = location.tnode_idx;
			TreeNode &tnode = tnodes[tnode_idx];
			if(tnode.feature == -1){
				location.shrinked = true;
			}else if(static_cast<uint32_t>(tnode.feature) < nr_field){

				if(prob.Z[tnode.feature][i].v < tnode.threshold)
					tnode_idx = 2*tnode_idx; 
				else
					tnode_idx = 2*tnode_idx+1; 
			}else{
				uint32_t const target_feature = static_cast<uint32_t>(tnode.feature-nr_field);
				bool is_one = false;
				for(uint64_t p = prob.SJP[i]; p < prob.SJP[i+1]; ++p) 
				{
					if(prob.SJ[p] == target_feature)
					{
						is_one = true;
						break;
					}
				}
				if(!is_one)
					tnode_idx = 2*tnode_idx; 
				else
					tnode_idx = 2*tnode_idx+1; 
			}
		}
	}
	
	// 用于计算gamma
	std::vector<std::pair<double, double>> 
		tmp(max_tnodes, std::make_pair(0, 0));
	for(uint32_t i = 0; i < nr_instance; ++i)
	{
		float const r = locations[i].r;
		uint32_t const tnode_idx = locations[i].tnode_idx;
		tmp[tnode_idx].first += r;
		tmp[tnode_idx].second += fabs(r)*(1-fabs(r));
	}

	for(uint32_t tnode_idx = 1; tnode_idx <= max_tnodes; ++tnode_idx)
	{
		double a, b;
		std::tie(a, b) = tmp[tnode_idx];
		tnodes[tnode_idx].gamma = (b <= 1e-12)? 0 : static_cast<float>(a/b);
	}

#pragma omp parallel for schedule(static)
	for(uint32_t i = 0; i < nr_instance; ++i)
		F1[i] = tnodes[locations[i].tnode_idx].gamma;// 重新更新F1的值
}

在参考文献A simple GBDT in Python中提供了Python实现的GBDT的版本。

参考文献

[1] Greedy Function Approximation: A Gradient Boosting Machine

[2] Gradient boosting machines, a tutorial

[3] 3 Idiots’ Approach for Display Advertising Challenge

[4] 《统计机器学习》

[5] GBDT：梯度提升决策树

[6] 随机森林&GBDT算法以及在MLlib中的实现

[7] A simple GBDT in Python