Short answer to you question:
when the algorithm it fits the residual (or the negative gradient) is it using one feature at each step (i.e. univariate model) or all features (multivariate model)?
The algorithm is using one features or all features depends on your set up. In my long answer listed below, in both decision stump and linear learner examples, they uses all features, but if you want, you can also fit a subset of features. Sampling columns (features) are viewed as reducing the variance of the model or increasing the "robustness" of the model, especially if you have large number of features.
In xgboost, for tree base learner, you can set colsample_bytree to sample features to fit in each iteration. For linear base learner, there are not such options, so, it should be fitting all features. In addition, not too many people use linear learner in xgboost or gradient boosting in general.
Long answer for linear as weak learner for boosting:
In most cases, we may not use linear learner as a base learner. The reason is simple: adding multiple linear models together will still be a linear model.
In boosting our model is a sum of base learners:
$$
f(x)=\sum_{m=1}^M b_m(x)
$$
where $M$ is number of iterations in boosting, $b_m$ is the model for $m^{th}$ iteration.
If the base learner is linear, for example, suppose we just run $2$ iterations, and $b_1=\beta_0+ \beta_1x$ and $b_2=\theta_0+ \theta_1x$, then
$$
f(x)=\sum_{m=1}^2 b_m(x)=\beta_0+ \beta_1x+\theta_0+ \theta_1x=(\beta_0+\theta_0)+ (\beta_1+ \theta_1)x
$$
which is a simple linear model! In other words, the ensemble model have the "same power" with the base learner!
More importantly, if we use linear model as base learner, we can just do it one step by solving the linear system $ X^T X \beta = X^T y $ instead of go though multiple iterations in boosting.
Therefore, people would like to use other models than linear model as base learner. Tree is a good option, since adding two trees is not equal to one tree. I will demo it with a simple case: decision stump, which is a tree with 1 split only.
I am doing a function fitting, where the data is generated by a simple quadratic function, $f(x,y)=x^2+y^2$. Here is the filled contour ground truth (left) and final decision stump boosting fitting (right).
Now, check the first four iterations.
Note, different from linear learner, the model in 4th iteration cannot be achieved by one iteration (one single decision stump) with other parameters.
So far, I explained, why people are not using linear learner as base learner. However, nothing prevent people to do that. If we use linear model as base learner, and restrict number of iterations, it is equal to solving a linear system, but limit the number of the iterations during the solving process.
The same example, but in 3d plot, the red curve are the data, and the green plane is the final fit. You can easily see, final model is a linear model, and it is z=mean(data$label) which is parallel to x,y plane. (You can think why? this is because our data is "symmetric", so any tilt of the plane will increase the loss). Now, check out what happened in first 4 iterations: the fitted model is slowly going up to the optimal value (mean).
Final conclusion, linear learner is not widely used, but nothing prevent people to use it or implement it in a R library. In addition, you can use it and limit number of iterations to regularize the model.
Related post:
Gradient Boosting for Linear Regression - why does it not work?
Is a decision stump a linear model?
