Gradient Boosting for Linear Regression - why does it not work? While learning about Gradient Boosting, I haven't heard about any constraints regarding the properties of a "weak classifier" that the method uses to build and ensemble model. However, I could not imagine an application of a GB that uses linear regression, and in fact when I've performed some tests - it doesn't work. I was testing the most standard approach with a gradient of sum of squared residuals and adding the subsequent models together.
The obvious problem is that the residuals from the first model are populated in such manner that there is really no regression line to fit anymore. My another observation is that a sum of subsequent linear regression models can be represented as a single regression model as well (adding all intercepts and corresponding coefficients) so I cannot imagine how that could ever improve the model. The last observation is that a linear regression (the most typical approach) is using sum of squared residuals as a loss function - the same one that GB is using.
I also thought about lowering the learning rate or using only a subset of predictors for each iteration, but that could still be summed up to a single model representation eventually, so I guess it would bring no improvement.
What am I missing here? Is linear regression somehow inappropriate to use with Gradient Boosting? Is it because the linear regression uses the sum of squared residuals as a loss function? Are there any particular constraints on the weak predictors so they can be applied to Gradient Boosting?
 A: 
What am I missing here?

I don't think you're really missing anything!

Another observation is that a sum of subsequent linear regression models can be represented as a single regression model as well (adding all intercepts and corresponding coefficients) so I cannot imagine how that could ever improve the model. The last observation is that a linear regression (the most typical approach) is using sum of squared residuals as a loss function - the same one that GB is using.

Seems to me that you nailed it right there, and gave a short sketch of a proof that linear regression just beats boosting linear regressions in this setting.
To be pedantic, both methods are attempting to solve the following optimization problem
$$ \hat \beta = \text{argmin}_\beta (y - X \beta)^t (y - X \beta) $$
Linear regression just observes that you can solve it directly, by finding the solution to the linear equation
$$ X^t X \beta = X^t y $$
This automatically gives you the best possible value of $\beta$ out of all possibilities.
Boosting, whether your weak classifier is a one variable or multi variable regression, gives you a sequence of coefficient vectors $\beta_1, \beta_2, \ldots$.  The final model prediction is, as you observe, a sum, and has the same functional form as the full linear regressor
$$ X \beta_1 + X \beta_2 + \cdots + X \beta_n = X (\beta_1 + \beta_2 + \cdots + \beta_n) $$
Each of these steps is chosen to further decrease the sum of squared errors.  But we could have found the minimum possible sum of square errors within this functional form by just performing a full linear regression to begin with.
A possible defense of boosting in this situation could be the implicit regularization it provides.  Possibly (I haven't played with this) you could use the early stopping feature of a gradient booster, along with a cross validation, to stop short of the full linear regression. This would provide a regularization to your regression, and possibly help with overfitting.  This is not particularly practical, as one has very efficient and well understood options like ridge regression and the elastic net in this setting.
Boosting shines when there is no terse functional form around.  Boosting decision trees lets the functional form of the regressor/classifier evolve slowly to fit the data, often resulting in complex shapes one could not have dreamed up by hand and eye.  When a simple functional form is desired, boosting is not going to help you find it (or at least is probably a rather inefficient way to find it).
A: OP is absolutely right. However I came across an algorithm called RegBoost by Li et al. (2020) which attempts to adapt Linear Regression (LR) to be used as the weak learner in Gradient Boosting by combining it with a non-linearity. This is done by constructing a decision tree of LR models based on the sign (+ve or -ve) of the error at each stage. They report results on par with Gradient Boosted Decision Trees (LightGBM) on 3 datasets, which seems promising. I haven't used/implemented it myself yet but it seems interesting so I will give a brief rundown of the core algorithm below.
Algorithm

*

*Train an LR model on training data

*Split the regression outputs into 2 categories: those above the target and those below the target

*Train an LR model on each of these splits with the new target being the residual error from the previous model

*Repeat steps 2 and 3 until a predetermined maximum depth is reached

Additional information

*

*The two key hyperparameters for tuning are the maximum depth of the tree and the learning rate (excluded here for simplicity)


*In order to determine the split when the target is not known (e.g. at inference time), the authors use a KNN model and assign a positive or negative class depending on the $k$ nearest neighbours from the training set at that particular split point


*In order to reduce dimensionality, the authors run a stepwise regression as a feature selection step prior to fitting each LR model. I personally think this step should not be part of the core algorithm. Feature selection may not be necessary depending on your number of features and you should be able to use any feature selection technique you want without affecting the algorithm as depicted above
A: The least squares projection matrix is given by
$X(X^{T}X)^{-1}X^{T}$
We can use this to directly obtain our predicted values $\hat{y}$, e.g.
$\hat{y} = X(X^{T}X)^{-1}X^{T}y $
Let's say you fit a regression and subsequently you calculate your residuals
$e = y - \hat{y} = y - X(X^{T}X)^{-1}X^{T}y $
And then you use this residual vector e as your new dependent variable in the next regression. Use the projection matrix again to directly calculate the predictions of this second regression and call these new predictions $\hat{y}_{2}$ :
$\hat{y}_{2} = X(X^{T}X)^{-1}X^{T}e \\
\quad = X(X^{T}X)^{-1}X^{T} (y - X(X^{T}X)^{-1}X^{T}y) \\
\quad =  X(X^{T}X)^{-1}X^{T}y - X(X^{T}X)^{-1}X^{T}X(X^{T}X)^{-1}X^{T}y \\
\quad = X(X^{T}X)^{-1}X^{T}y - X(X^{T}X)^{-1}X^{T}y \\
\quad = 0 $
A reason for this is that by construction the residual vector e from the initial regression is orthogonal to the X Space i. e. $\hat{y}$ is a orthogonal projection from y onto the X space (you'll find nice pictures visualizing this in the literature).
This means the simple approach of fitting a regression and then fitting a new regression on the residuals from the first regression will not result in anything senseful because X is entirely uncorrelated with e.
I write this because you said there is not really a new line to fit which corresponds to the derivations above.  
