Gradient Boosting as Functional Gradient Descent

In ordinary parametric gradient descent, we optimize a loss $L(\theta)$ over a parameter vector $\theta \in \mathbb{R}^p$ by repeatedly moving in the direction $-\nabla_\theta L(\theta)$. Gradient boosting applies the identical logic, but replaces the parameter vector with the vector of function values $\mathbf{F} = (F(x_1), \dots, F(x_n))$ evaluated at the training points. The empirical risk

is viewed as a function of the $n$-dimensional vector $\mathbf{F}$, not of any finite parameter set. Its gradient at the current predictions is

The ideal update would move each $F(x_i)$ by $-g_i$, exactly like unconstrained gradient descent — but this "point cloud" of updates is not a function we can evaluate at new test inputs. Gradient boosting's key move is to project this ideal but non-generalizing gradient onto the space of outputs of a weak learner (e.g. shallow regression trees), by fitting the weak learner $h_m(x)$ to the pseudo-residuals $-g_i$ using ordinary least squares (or whatever loss the weak learner itself minimizes). The fitted $h_m$ is the best available finite-dimensional approximation to the true functional steepest-descent direction, and it generalizes to new $x$ because it is a genuine function, not a lookup table of per-example corrections. Having found a direction $h_m$, boosting performs a (possibly approximate) line search for a step size $\rho_m$ that minimizes $\sum_i \ell(y_i, F_{m-1}(x_i) + \rho\, h_m(x_i))$, and updates

The final model $F_M = F_0 + \sum_{m=1}^M \rho_m h_m$ is an additive ensemble, structurally identical to AdaBoost's weighted combination of weak classifiers, but now derived for any differentiable loss rather than being tied to the exponential loss. When $\ell$ is squared error, the pseudo-residual $-g_i$ reduces exactly to the ordinary residual $y_i - F_{m-1}(x_i)$, so "fit a tree to the residuals" — the folk description of boosting — is the squared-error special case of a much more general procedure.

The generic algorithm requires three ingredients: a loss $\ell(y, F)$ differentiable in its second argument, a family of weak learners $\mathcal{H}$ (typically shallow regression trees), and a fitting procedure that can approximately minimize squared error against arbitrary targets (this is what "fit a weak learner" means operationally, regardless of the outer loss).

Step 1 — Initialize. Set $F_0(x) = \arg\min_c \sum_i \ell(y_i, c)$, a constant model (e.g. the mean for squared error, the log-odds of the base rate for logistic loss).

Step 2 — For each stage $m = 1, \dots, M$: 1. Compute pseudo-residuals $r_{im} = -\left.\partial \ell(y_i, F)/\partial F\right|_{F=F_{m-1}(x_i)}$ for every training example. 2. Fit a weak learner $h_m$ to the training set $\{(x_i, r_{im})\}_{i=1}^n$ by least-squares regression — regardless of the outer loss, this inner fit is always squared-error regression of the weak learner onto the residuals. 3. Choose a step size $\rho_m = \arg\min_\rho \sum_i \ell(y_i, F_{m-1}(x_i) + \rho\, h_m(x_i))$ (exact line search, or a fixed shrinkage constant $\nu \in (0,1]$ used instead of a full search). 4. Update $F_m(x) = F_{m-1}(x) + \nu\, \rho_m\, h_m(x)$.

Step 3 — Output $F_M$, the sum of the initial constant and all weighted weak learners. ``` Algorithm: GradientBoosting Input: data {(x_i, y_i)}_{i=1..n}, loss L(y, F), weak learner class H, number of stages M, shrinkage rate nu in (0, 1] Output: additive model F_M(x) F_0(x) := argmin_c sum_i L(y_i, c) # constant baseline for m = 1 to M:

# 1. functional gradient (pseudo-residuals) for i = 1 to n: r_i := - dL(y_i, F)/dF evaluated at F = F_{m-1}(x_i)

# 2. fit weak learner to the negative gradient (least squares) h_m := FitLeastSquares( {(x_i, r_i)}_{i=1..n}, H )

# 3. line search for optimal step length along h_m rho_m := argmin_rho sum_i L( y_i, F_{m-1}(x_i) + rho * h_m(x_i) )

# 4. shrink and update F_m(x) := F_{m-1}(x) + nu rho_m h_m(x) return F_M ``` Two loss-specific instantiations illustrate the generality:

• Squared error, $\ell(y,F) = \tfrac12 (y-F)^2$: gradient is $F - y$, so $r_i = y_i - F_{m-1}(x_i)$, the ordinary residual. This recovers "L2Boosting."
• Logistic loss for binary classification with $y \in \{-1,+1\}$, $\ell(y,F) = \log(1+e^{-2yF})$: the negative gradient is $r_i = \dfrac{2y_i}{1+e^{2y_i F_{m-1}(x_i)}}$, a smoothly-weighted margin residual, and $F_M(x)$ is passed through a sigmoid to obtain class probabilities. This is the algorithm underlying gradient-boosted logistic classifiers (and is closely related to, but distinct from, AdaBoost's exponential-loss update).

The diagram below shows the computational structure: predictions flow forward, the loss gradient flows backward to produce residual targets, and a new weak learner is trained on those targets at every stage.

Time complexity. Each of the $M$ stages requires (i) computing $n$ pseudo-residuals in $O(n)$ time, (ii) fitting a weak learner, and (iii) a line search. If the weak learner is a regression tree of maximum depth $d$ built by the standard greedy split-finding procedure, fitting one tree costs $O(n \cdot p \cdot d)$ where $p$ is the number of features (each of the $d$ levels considers, for each of $p$ features, a sort/scan over $n$ points; with pre-sorted features this is $O(n\log n \cdot p)$ once, reused across levels via level-order data structures, giving $O(n \cdot p \cdot d)$ amortized in practice). Line search for many common losses (squared error, exponential) has a closed form and costs $O(n)$; for losses without closed form (e.g. absolute or Huber loss for arbitrary residual shapes), a small number of Newton or golden-section steps each costing $O(n)$ suffices. Total training time is

Space complexity. $O(n)$ for the residual vector and current predictions per stage, plus $O(M \cdot \text{size}(h_m))$ to store the ensemble; a depth-$d$ tree needs $O(2^d)$ nodes, so total model storage is $O(M \cdot 2^d)$, typically small since $d$ is kept shallow (3–8) by design.

Correctness — why each stage decreases the loss. The argument mirrors why parametric gradient descent decreases a convex differentiable objective, adapted to the function-space setting. Treat $L$ as a function of the vector $\mathbf{F} \in \mathbb{R}^n$ of predictions at the training points. A first-order Taylor expansion around $\mathbf{F}_{m-1}$ gives, for a candidate update $\mathbf{F}_{m-1} + \rho\,\mathbf{h}$,

where $\mathbf{g} = \nabla L(\mathbf{F}_{m-1})$ is the vector of true pointwise gradients ($g_i$ as defined above). Among all unit-length directions $\mathbf{h}$, the choice $\mathbf{h} = -\mathbf{g}/\|\mathbf{g}\|$ makes $\mathbf{g}^\top\mathbf{h}$ most negative — this is the standard steepest-descent argument, applied here coordinatewise over the $n$ training points rather than over parameter coordinates. Since we cannot realize the ideal $-\mathbf{g}$ exactly (it isn't a function of $x$ we can evaluate elsewhere), we instead fit $h_m \in \mathcal{H}$ by least squares to $-\mathbf{g}$:

This maximizes the correlation of the achievable direction $\mathbf{h}_m = (h_m(x_1), \dots, h_m(x_n))$ with the ideal descent direction $-\mathbf{g}$, subject to the constraint that $h_m$ lies in the weak-learner family. Provided $\mathbf{g}^\top \mathbf{h}_m < 0$ (i.e., the weak learner is even weakly correlated with the true negative gradient — the same "better than random" condition required in AdaBoost), the line search finds a $\rho_m > 0$ small enough that the first-order Taylor approximation is accurate and

Hence the sequence of training losses $L(F_0) > L(F_1) > \cdots$ is strictly decreasing as long as the weak learner can achieve nonzero correlation with the residual — exactly analogous to the weak learning assumption in boosting theory. If $L$ is convex in $\mathbf{F}$ (true for squared error, logistic loss, and exponential loss) and the weak learner family is rich enough that the achievable directions eventually span the relevant subspace, the greedy stagewise procedure converges to the global minimum of $L$ over the additive span of $\mathcal{H}$; more generally it is a form of coordinate-free steepest descent and inherits the standard non-increasing property of gradient descent under a sufficiently small step size (this is the role shrinkage $\nu$ plays: it enforces the small-step regime where the first-order approximation used above remains valid, at the cost of needing more stages $M$ to converge).

Honest tradeoffs. The procedure is greedy and stagewise: once $h_m$ is fixed it is never revisited, so the algorithm does not jointly optimize all weak learners together, unlike, say, fully retraining a neural network end to end. This greediness is precisely what makes it computationally tractable — $O(Mnpd)$ instead of a combinatorial search over ensembles — but it means gradient boosting can get stuck fitting noise late in training if $M$ is too large or $\nu$ too close to 1, which is why regularization (subsampling, shrinkage, tree-depth limits, early stopping on a validation set) is not optional in practice but integral to the method as actually deployed.

Consider a tiny 1-D regression problem with squared-error loss, $\ell(y,F) = \tfrac12(y-F)^2$, so residuals are just $y - F$. Weak learners are decision stumps (single-split regression trees). Data:

Stage 0. $F_0(x) = \bar y = (2.0+2.2+5.0+5.4)/4 = 3.65$.

Stage 1. Residuals $r_i = y_i - F_0$:

Fit a stump: the best split is between $x=2$ and $x=3$ (it separates the negative residuals from the positive ones). Left region $\{x\le 2\}$ predicts the mean residual $(-1.65-1.45)/2 = -1.55$; right region $\{x>2\}$ predicts $(1.35+1.75)/2 = 1.55$. So

Exact line search for squared error with a stump target reduces to $\rho_1=1$ (the stump already predicts the mean residual per region, which is the exact least-squares fit). Apply shrinkage $\nu = 0.5$:

For $x\le2$: $F_1 = 3.65 - 0.775 = 2.875$. For $x>2$: $F_1 = 3.65+0.775=4.425$.

Stage 2. New residuals against $F_1$:

The same split location is again optimal (residual sign still changes between $x=2$ and $x=3$). New stump: left mean $= (-0.875-0.675)/2 = -0.775$, right mean $=(0.575+0.975)/2=0.775$. Update with $\nu=0.5$:

Observe the pattern. Each stage halves the remaining residual magnitude (because $\nu=0.5$ and the stump exactly captures the group means), converging geometrically toward the true group means $2.1$ (average of $y_1,y_2$) and $5.2$ (average of $y_3,y_4$):

This numerically confirms strict monotone decrease of training loss, converging toward the piecewise-constant function that perfectly separates the two clusters — exactly the additive-model global optimum over stumps split at $x=2.5$.

Gradient boosting's core formulation — casting boosting as functional gradient descent for arbitrary differentiable losses — was developed by Jerome H. Friedman at Stanford in the late 1990s, unifying the AdaBoost framework (loss-specific, exponential loss) into a general-purpose machine that works with any smooth loss and any weak learner supporting least-squares fitting.

Gradient boosted trees are the default choice for structured/tabular prediction problems across industry: credit risk scoring, fraud detection, click-through-rate and conversion prediction in online advertising, search ranking (e.g. learning-to-rank systems built on gradient boosting), and demand forecasting in logistics. Libraries such as XGBoost, LightGBM, and CatBoost implement the functional gradient descent framework described here, augmented with second-order Newton corrections, regularized leaf weights, and engineering optimizations (histogram binning, GPU acceleration, categorical feature handling) — but the mathematical core remains: fit a weak learner to the negative gradient, take a controlled step, repeat. Kaggle-style competitions on tabular data are dominated almost exclusively by gradient boosting variants, and production ML pipelines routinely combine boosted trees with other model families, a practice formalized under Stacking and Model Blending.

Ask the tutor: How does the line-search step size interact with shrinkage, and why do modern libraries often replace exact line search with a fixed learning rate plus second-order (Newton) leaf weights?