How do Newton-Raphson updates work in Gradient Boosting?

Question

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?

As reference for the Newton-Raphson Algorithm I've been reading the following links: https://www.stat.washington.edu/adobra/classes/536/Files/week1/newtonfull.pdf https://math.stackexchange.com/questions/1153655/newtons-method-vs-gradient-descent-with-exact-line-search

Answer 1

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.

Answer 2

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} $