the gradient boosting (GBM) algorithm computes the residuals (negative gradient) and then fit them by using a regression tree with mean square error (MSE) as the splitting criterion. How is that different from the XGBoost algorithm?
Both indeed fit a regression tree to minimize MSE w.r.t. a pseudo-response variable in every boosting iteration. This pseudo response is computed based on the original binary response, and predictions from the regression trees of previous iterations.
The two methods differ in how the pseudo response is computed. GBM uses a first-order derivative of the loss function at the current boosting iteration, while XGBoost uses both the first- and second-order derivatives. The latter is also known as Newton boosting.
gbm uses gradient boosting, by default. At each boosting iteration, the regression tree minimizes the least squares approximation to the negative gradient.
Is the only difference between GBM and XGBoost the regularization terms or does XGBoost use another split criterion to determine the regions of the regression tree?
I do not fully understand what you mean by the regularization term. The split criterion for the regression tree indeed differs, because XGBoost computes the pseudo-response variable differently. At each boosting iteration, it fits a regression tree to minimize a weighted least squares approximation. The pseudo response is calculated so that the more extreme the predictions from previous trees (i.e., linear predictor far from zero; predicted probability close to either 0 or 1), the less influence the observation will receive in the current iteration. Thus, the gradient is weighted by the uncertainty of predictions from previous trees.
Sigrist (2021) provides a good discussion of the differences, both in terms of predictive performance and the functions being optimized.
Sigrist, F. (2021). Gradient and Newton boosting for classification and regression. Expert Systems With Applications, 167, 114080. https://arxiv.org/abs/1808.03064