This is a very long question, but perhaps people who are trying to deeply understand the Gradient Boosting Machine algorithm will think it's interesting.

I've been working on understanding the thought process behind the gradient boosting machine (GBM) algorithm for binary classification. There is a step which is not perfectly clear to me, and my question (which I'll ask below) will be to clarify that particular step. But first, I'll explain my current understanding of how the GBM algorithm for binary classification works and how the algorithm might be discovered. Feel free to skip directly to the question at the end.


Let $\sigma: \mathbb R \to \mathbb R$ be the logistic function, defined by $$ \sigma(u) = \frac{ e^u}{1 + e^u}. $$ The logistic function is useful because it converts a real number into a probability. Note that $$ \sigma'(u) = \sigma(u)(1 - \sigma(u)). $$

Let $\ell$ be the cross-entropy loss function defined by $$ \ell(p,q) = -p \log(q) - (1 - p) \log(1 - q). $$ The inputs $p, q \in \mathbb R$ are required to be between $0$ and $1$. We can think of $p$ and $q$ as being probabilities, and $\ell(p,q)$ measures how closely $p$ and $q$ agree. Let's also define $$ L(p,u) = \ell(p, \sigma(u)). $$ Here the input $p$ is required to be between $0$ and $1$, but $u$ can be any real number. Note that $$ \partial_u L(p,u) = \sigma(u) - p, $$ which is a nicely simple formula. (Here $\partial_u L$ denotes the partial derivative of $L$ with respect to $u$.)

GBM algorithm for binary classification:

We are given a list of feature vectors $x_1, \ldots, x_N \in \mathbb R^d$ and corresponding labels $y_1, \ldots, y_N \in \{0, 1 \}$. As is typical in binary classification algorithms, our goal is to find a function $f: \mathbb R^d \to \mathbb R$ such that for each training example the label $y_i$ agrees closely with $\sigma(f(x_i))$. More precisely, we hope that $$ \sum_{i=1}^N \ell(y_i, \sigma(f(x_i)) $$ is small.

In the gradient boosting machine approach, we assume that $f$ is a sum of regression trees: $$ f(x) = \sum_{m=1}^M T_m(x), $$ where each function $T_m: \mathbb R^d \to \mathbb R$ is a regression tree. We grow the regression trees $T_m$ in a forward stagewise manner. Specifically, we define $f_0 = 0$, and then for $m = 1, \ldots, M$ we define $f_m = f_{m-1} + T_m$, where $T_m$ is a regression tree which is chosen so that $$ \tag{1} \sum_{i=1}^N L(y_i, f_{m-1}(x_i) + T_m(x_i)) $$ is small. To simplify the computation of $T_m$, we make the first-order approximation \begin{align} L(y_i, f_{m-1}(x_i) + T_m(x_i)) &\approx L(y_i, f_{m-1}(x_i)) + \partial_u L(y_i, f_{m-1}(x_i)) T_m(x_i) \\ &= L(y_i, f_{m-1}(x_i)) + (\sigma(f_{m-1}(x_i)) - y_i) T_m(x_i). \end{align} This approximation is only accurate when $T_m(x)$ is small. For that reason, we will also modify the objective function (1) by including penalty terms $(1/2t)T_m(x_i)^2$ which encourage $T_m(x_i)$ to be small. Thus, we replace the objective function (1) with the new objective function $$ \tag{2} \sum_{i=1}^N L(y_i, f_{m-1}(x_i)) + (\sigma(f_{m-1}(x_i)) - y_i) T_m(x_i) + \frac{1}{2t} T_m(x_i)^2. $$ A typical choice for the parameter $t$ is $t = 1$, but smaller values are more conservative and might work better.

Completing the square, we see that the objective function in (2) can be written as $$ \tag{2'} \sum_{i=1}^N \frac{1}{2t} \left[t\big(y_i - \sigma(f_{m-1}(x_i))\big) - T_m(x_i)\right]^2 + C_i $$ where the terms $C_i$ do not depend on $T_m(x_i)$. Growing a regression tree $T_m$ to minimize (2') is equivalent to growing $T_m$ to minimize $$ \sum_{i=1}^N \left[t \big(y_i - \sigma(f_{m-1}(x_i)) \big) - T_m(x_i)\right]^2. $$ In other words, to find $T_m$ we simply fit a regression tree to the residuals $$ r_i = t \big(y_i - \sigma(f_{m-1}(x_i)) \big). $$ The resulting regression tree partitions $\mathbb R^d$ into regions $R_1, \ldots, R_{J_m}$. If $x \in \mathbb R_j$, then $T_m(x) = \gamma_j$. The regions $R_j$ and scalars $\gamma_j$ are now known.

There is one final step in the construction of $T_m$ in the GBM algorithm. After the regions $R_j$ are found, the scalars $\gamma_j$ are then recomputed as follows: $$ \tag{3} \gamma_j = \arg \min_\gamma \sum_{x_i \in R_j} L(y_i, f_{m-1}(x_i) + \gamma). $$

Question: How exactly do we solve for $\gamma_j$ in (3)? Setting the derivative with respect to $\gamma$ equal to $0$, we find that \begin{align} & \sum_{x_i \in R_j} \partial_u L(y_i, f_{m-1}(x_i) + \gamma) = 0 \\ \tag{4} \implies & \sum_{x_i \in R_j} \sigma(f_{m-1}(x_i) + \gamma) - y_i = 0. \end{align} But this is a nonlinear equation for $\gamma$, and I think there is no closed-form solution. So what do we do?

One option would be to compute $\gamma$ using Newton's method. Another less expensive option would be to use the approximation $$ \tag{5} \sigma(u_i + \gamma) \approx \sigma(u_i) + \sigma'(u_i) \gamma = \sigma(u_i) + \sigma(u_i)(1 - \sigma(u_i)) \gamma. $$ Setting $u_i = f_{m-1}(x_i)$ and making the approximation (5), equation (4) reduces to \begin{align} &\sum_{x_i \in R_j} \sigma(u_i) + \sigma(u_i)(1 - \sigma(u_i)) \gamma - y_i = 0 \\ \implies & \gamma = \gamma_j = \frac{\sum_{x_i \in R_j} y_i - \sigma(u_i)}{\sum_{x_i \in \mathbb R_j} \sigma(u_i)(1 - \sigma(u_i))}. \end{align} Is this the correct formula that we should use for $\gamma_j$ in the GBM algorithm for binary classification?

For reference, here is the GBM algorithm as presented in The Elements of Statistical Learning by Hastie et al. Step 2 (c) is the step that I am asking about in this question.

enter image description here

  • 3
    $\begingroup$ A second-order Taylor series is a quadratic which can be optimized directly. This is the approach taken in xgboost, for example. xgboost.readthedocs.io/en/latest/tutorials/model.html Scroll down to "Structure Score" $\endgroup$ – Sycorax Jan 11 '20 at 5:17
  • $\begingroup$ @SycoraxsaysReinstateMonica Thanks, I just spent some time reading that page. Whoever wrote that XGBoost documentation is an exceptionally clear writer / clear thinker. $\endgroup$ – littleO Jan 11 '20 at 5:58
  • 1
    $\begingroup$ Believe me, it’s seen many drafts. $\endgroup$ – Sycorax Jan 11 '20 at 6:10
  • 3
    $\begingroup$ I can't answer this question, but I appreciate the clear presentation of your understanding of GBM. $\endgroup$ – Noah Jan 12 '20 at 2:39

I think your reasoning is basically correct: we could use an iterative solver like Newton's method to get a really good approximation to the optimal $\gamma$, but that's a lot of work.

The software uses a second-order Taylor expansion of $(3)$ because quadratics are easy to optimize. (Described in https://xgboost.readthedocs.io/en/latest/tutorials/model.html; scroll down to the section "Structure Score.") Because of its simplicity, I wouldn't be surprised if this is the approach taken in most boosting libraries.

An extensive answer about xgboost can be found here: XGBoost Loss function Approximation With Taylor Expansion

  • $\begingroup$ Thanks. Very helpful. In Hastie's presentation of the GBM algorithm (Algorithm 10.3 in Hastie), once the tree $T_m$ is grown (in step 2 (b) above), there is a "polishing" step 2 (c) where the constants $\gamma_{jm}$ are recomputed with the regions $R_{jm}$ held fixed. In step 2(c), as Hastie describes it in algorithm 10.3, we do not make the same approximation that led to step 2 (b). From the XGBoost documentation, I get the impression that the polishing step 2 (c) does not appear at all in the XGBoost algorithm. Is that correct? $\endgroup$ – littleO Jan 15 '20 at 6:56

I would say that expand equation 3 using taylor series until seconnd order polynomial and then set the derivative of this modified equation wrt gamma equal to zero to get the gamma optimal..this is equivalent to performing only one iteration of newtons method for optimization. you should get the gamma as sum(residuals)/sum(p*(1-p)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.