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.
Background:
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.
