I want to use the Gradient Boosting algorithm with exponential loss function and I am struggling to understand how to use the Newton-Raphson update step for predictions. In python's sklearn GradientBoostingClassifier the update step is the following:

    numerator = np.sum(y_ * sample_weight * np.exp(-y_ * pred))
denominator = np.sum(sample_weight * np.exp(-y_ * pred))

# prevents overflow and division by zero
if abs(denominator) < 1e-150:
tree.value[leaf, 0, 0] = 0.0
else:
tree.value[leaf, 0, 0] = numerator / denominator


The numerator is the negative sum of the partial first derivatives of the exponential loss function $$\ L(pred) = \sum(e^{-y \cdot pred}), \text{ and } \\ \ numerator = \sum\left(-\frac{\partial}{\partial pred} L(pred)\right) = \sum(y \cdot e^{-y \cdot pred}).$$ The denominator is the sum of the partial second derivatives of the exponential loss function \begin{align} \ denominator = \sum \left( \frac{\partial^2}{\partial pred^2} L(pred)\right) = \sum(y^2 \cdot e^{-y \cdot pred}) = \sum(e^{-y \cdot pred}), \text{ since } y^2 = 1. \end{align}

But according to Newton-Raphson algorithm the update in pred should be:

\begin{align} \ pred = pred - Hessian(L(pred))^{-1} \cdot \nabla L(pred) \\ \end{align}

where Hessian is the matrix of second-order partial derivatives.

Why does Python sums over the gradient and over the Hessian and then take the ratio of the two as the update step in predictions?

Tree-boosting builds ensembles of trees in an iterative way. In every iteration, a new tree is added to the ensemble conditional on the current ensemble of trees in such a way that the empirical risk (aka the training loss) is minimized. Since the latter can usually not be done in closed form, an approximate minimization is done. Broadly speaking, a new tree is found using either functional gradient descent or a functional version of Newton's method, which is conceptually similar to Newton's method in finite dimensions (https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization). Specifically, for Newton boosting, the empirical risk is replaced by a (functional) second order Taylor approximation, and a tree is found by minimizing this approximation. The latter can be done using weighted least squares. Note that for standard regression with a squared error, gradient and Newton boosting are equivalent.

For more details, this preprint https://arxiv.org/abs/1808.03064 explains the difference between Newton and gradient booting.

Disclaimer: I am the author of the article.

• Welcome to the site. We are trying to build a permanent repository of high-quality statistical information in the form of questions & answers. Thus, we're wary of link-only answers. Can you post a full citation & a summary of the information at the link? Feb 20 '19 at 14:26

I have formulated the problem wrong. According to

Greedy Function Approximation: A Gradient Boosting Machine

the Gradient Boosting algorithm applied to the exponential loss function includes the following steps:

1. Initialize $$F_0(x) = argmin_{\rho} = \sum e^{-y\rho}$$ and repeat the following steps 2 to 5 for the number of the week learners $$m$$

2. Compute the negative gradient of the exponential loss function $$\ L(y,F) = \sum e^{-yF}$$

$$\ res = -\frac{\partial L(y,F)}{\partial F} = \ ye^{-yF}$$ where $$y \in (-1,1)$$

3. Fit a weak learner $$h$$ (regression tree) to the negative gradient $$res$$

4. In each terminal node of the tree compute the optimal step size $$\rho$$

$$\rho = argmin_\rho = L(y,F_{m-1} + \rho h)$$

5. Update $$F$$ as $$F_m = F_{m-1} + \rho h$$

By applying a single Newton-Raphson step in order to find $$\rho$$ we have:

$$G(\rho) = L(y,F + \rho h) = L(y,F) + \rho h^T \frac{\partial L}{\partial F} + \frac{1}{2} \rho^2 h^T \frac{\partial {^2} L}{\partial F^2}h,$$ by applying a Taylor second order approximation

Computing the derivative of $$G(\rho)$$ and setting it to zero we have: $$\frac{dG}{d\rho} = h^T \frac{\partial L}{\partial F} + \rho h^T \frac{\partial {^2} L}{\partial F^2}h = 0$$

Hence, $$\rho = -\frac{h^T \frac{\partial L}{\partial F}}{h^T \frac{\partial L^2}{\partial {^2} F}h} = -\frac{h \sum \frac {\partial L}{\partial F}}{h^2 \sum \frac {\partial {^2} L}{\partial F^2}} = -\frac{\sum \frac {\partial L}{\partial F}}{h \sum \frac {\partial {^2} L}{\partial F^2}}$$ since $$h$$ vector has the same value across all observations (that was the tricky part I have missed before)

Finally, $$F_m = F_{m-1} -\frac{\sum \frac {\partial L}{\partial F}}{\sum \frac {\partial {^2} L}{\partial F^2}}$$ which is the same as the python code

where $$pred = F, numerator = \sum -\frac{\partial L}{\partial F} = \sum ye^{-yF}, denominator = \sum \frac{\partial {^2} L}{\partial F^2} = \sum e^{-yF}$$