How is gradient boosting like gradient descent? I am reading the useful Wikipedia entry on gradient boosting (https://en.wikipedia.org/wiki/Gradient_boosting), and try to understand how / why we can approximate the residuals by the steepest descent step (also called the pseudo-gradient). Can anyone give me the intuition on how the steepest descent is linked / similar to the residuals? Help much appreciated! 

 A: Suppose we are in the following situation.  We have some data $\{ x_i, y_i \}$, where each $x_i$ can be a number or vector, and we would like to determine a function $f$ that approximates the relationship $f(x_i) \approx y_i$, in the sense that the least squares error:
$$ \frac{1}{2} \sum_i (y_i - f(x_i))^2 $$
is small.
Now, the question enters of what we would like the domain of $f$ to be.  A degenerate choice for the domain is just the points in our training data.  In this case, we may just define $f(x_i) = y$, covering the entire desired domain, and be done with it.  A round about way to arrive at this answer is by doing gradient descent with this discrete space as the domain.  This takes a bit of a change in point of view.  Let's view the loss as a function of the point true $y$ and the prediction $f$ (for the moment, $f$ is not a function, but just the value of the prediction)
$$ L(f; y) = \frac{1}{2} (y - f)^2 $$
and then take the gradient with respect to the prediction
$$ \nabla_f L(f; y) = f - y $$
Then the gradient update, starting from an initial value of $y_0$ is
$$ y_1 = y_0 - \nabla_f (y_0, y) = y_0 - (y_0 - y) = y $$
So we recover our perfect prediction in a gradient step with this setup, which is nice!
The flaw here is, of course, that we want $f$ to be defined at much more than just our training data points.  To do this, we must make a few concessions, for we are not able to evaluate the loss function, or its gradient, at any points other than our training data set.  
The big idea is to weakly approximate $\nabla L$.  
Start with an initial guess at $f$, almost always a simple constant function $f(x) = f_0$, this is defined everywhere.  Now generate a new working dataset by evaluating the gradient of the loss function at the training data, using the initial guess for $f$:
$$ W = \{ x_i, f_0 - y \} $$
Now approximate $\nabla L$ by fitting weak learner to $W$.  Say we get the approximation $F \approx \nabla L$.  We have gained an extension of the data $W$ across the entire domain in the form of $F(X)$, though we have lost precision at the training points, since we fit a small learner.
Finally, use $F$ in place of $\nabla L$ in the gradient update of $f_0$ over the entire domain:
$$ f_1(x) = f_0(x) - F(x) $$
We get out $f_1$, a new approximation of $f$, a bit better than $f_0$.  Start over with $f_1$, and iterate until satisfied.
Hopefully, you see that what is really important is approximating the gradient of the loss.  In the case of least squares minimization this takes the form of raw residuals, but in more sophisticated cases it does not.  The machinery still applies though.  As long as one can construct an algorithm for computing the loss and gradient of loss at the training data, we can use this algorithm to approximate a function minimizing that loss.
