Which data are used at each step of Stochastic Gradient Boosting? Subsample of the original training set or gradient of the loss function? The bag.fraction parameter in SGB controls the size of the random subsample of the original training set on which each successive weak learner is fitted:

At each iteration a subsample of the training data is drawn at random
  (without replacement) from the full training data set. This randomly
  selected subsample is then used, instead of the full sample, to fit
  the base learner and compute the model update for the current
  iteration.

(Quote from the original SGB paper by Friedman).
What I don't get is how this is compatible with the fact that in Gradient Boosting, each successive base model is fitted directly on the residuals (or gradient of the loss function) of the current model.
In other words, how can we fit each new base model on the residuals of the current model and at the same time on a random subsample of the original training set? What is the link between these two key steps?
If someone with enough Python literacy could find the information from the code, that'd be great. Also, the R code is not useful as it is only a wrapper for a C implementation.
 A: I've implemented this here in a way that performs comparably to other implementations at least in my testing.  
At each step you have a residual R calculated across the entire data set (before any base learners are fit its the residual from predicting the average).
You subsample without replacement and fit the learner to the values of R for the points in that subsample.
Then you apply that learner to get a prediction across the entire data set and update R using that prediction.
So you're still updating the residual across the entire data set, you're just randomizing the process of fitting each base learner a bit.
A: I tried to decipher the steps in the original Friedman 2002 paper. If I am not wrong, it appears that we actually both deal with the current residuals and the original training set. Especially, the algorithm is as follows:

Here is what I have understood, step by step. Please tell me if there is something amiss. See (*) below for some notation.


*

*line 1: $F_{0}(x)$ is just an initial guess (and that's not what I'm interested in here)

*line 4: we compute the current "pseudo-residuals" for a random subsample of the training set $(y_{\pi(i)},x_{\pi(i)})$ as the negative gradient of the loss function. By current, I mean that we want to add the $m^{th}$ base learner, and that the model so far is a sequence of $m-1$ base learners.

*line 5: the $m^{th}$ base learner is a regression tree that partitions the space of the $x_{\pi(i)}$'s into $L$ disjoint regions $R_{lm}$ so that the within-region variance in terms of pseudo-residuals (computed at line 4) is minimized (this is how a tree works). The predictions $\bar{y}_{lm}$ made by the tree are the means of the current pseudo-residuals in each region. This is what I was referring to in my question by: how can we fit each new base model on the residuals of the current model...

*line 6: within each region $R_{lm}$ found at line 5, we pick the expansion coefficient $\gamma_{lm}$ that minimizes the loss function (not its gradient, so we're not dealing with residuals here) evaluated on the random subsample of the target $y_{\pi(i)}$ and the corresponding predictions of the current model $F_{m-1}(x_{\pi(i)})$. This is what I was referring to in my question by: ...and at the same time on a random subsample of the original training set?

*line 7: we add the (contribution of) the new base learner to the sequence as its expansion coefficients for each region multiplied by the learning rate. This final step is still a bit weird to me, as we are not directly adding the new base learner itself (unlike what the widespread belief about Boosting might suggest). But I guess this is equivalent, since the $\gamma_{lm}$'s inherently encode the information about the new base learner (they were estimated within each leaf of the tree).


*Some notation:


*

*$(y_{i},x_{i})$ is the training set (target, predictor)

*the goal of Boosting (like any other ML algo) is to approximate the
function $F(x)$ that maps $x$ to $y$

*The approximation of $F$ proposed by Boosting takes the form
$F(x)=\sum_{m=0}^{M} \beta_{m}h(x,a_{m})$, where the $\beta$'s are
expansion coefficients and the $a$'s are parameters of the base
learners $h$'s (which are regression trees in Tree Boosting). So the $a$'s are the splits of the tree (optimizing the $a$'s=building the tree). $m$ is the number of base learners in the sequence. For trees, there is one $\beta_{m}$ per leaf, they are referred to as the $\gamma_{lm}$'s

*$\Psi$ is any differentiable loss function. We want the $F$ that
minimizes the expected value of $\Psi(y,F(x))$.

*At each step, $F_{m}(x)=F_{m-1}(x)+\nu\beta_{m}h(x,a_{m})$

A: The concept of residuals makes sense on any set of data for which:


*

*You can use your model to make predictions and

*You have a true response value on each record


So, if you, for example, split your full training data in two halves and then train on the first half, you can still predict on the second half and compute the residuals there.
I want to note though, that each base learner is fit to residuals in only the gaussian case of a gradient booster.  In the general case the weak learners are fit to the gradient of the loss function, which in the gaussian case turns out to be equivalent, but in all others are not.  You make this note as well, but I thought it was worth repeating.
A: I think you are confusing the features with the dependent variable. In gradient boosting models you fit the tree on the current gradient (which coincides with residuals for Gaussian models) using already existing features. Features do not change, only the gradient is changing at each step, because of the updating. Gradient is computed for each observation in the training set, so when you randomly select half of the training set it is easy to track which gradients to use. See Figure 2 in R package gbm vignette for better explanation, or the gradient tree boosting algorithm in section 10 in Elements of Statistical Learning.
