Why is xgboost so much faster than sklearn GradientBoostingClassifier?
Since you mention "numeric" features, I guess your features are not categorical and have a high arity (they can take a lot of different values, and thus there are a lot of possible split points). In such a case, growing trees is difficult since there are [a lot of features $\times$ a lot of split points] to evaluate.

My guess is that the biggest effect comes from the fact that XGBoost uses an approximation on the split points. If you have a continuous feature with 10000 possible splits, XGBoost consider only "the best" 300 splits by default (this is a simplification). This behavior is controlled by the sketch_eps parameter, and you can read more about it in the doc. You can try lowering it and check the difference it makes. Since there is no mention of it in the scikit-learn documentation, I guess it is not available. You can learn what XGBoost method is in the their paper (arxiv).

XGBoost also uses an approximation on the evaluation of such split points. I do not know by which criterion scikit learn is evaluating the splits, but it could explain the rest of the time difference.



Adressing Comments

Regarding the evaluation of split points


  However, what did you mean by "XGBoost also uses an approximation on the evaluation of such split points"? as far as I understand, for the evaluation they are using the exact reduction in the optimal objective function, as it appears in eq (7) in the paper.


In order to evaluate the split point, you would have to compute $L(y,H_{i-1}+h_i)$ where $L$ is the cost function, $y$ the target, $H_{i-1}$ the model built until now, and $h_i$ the current addition. Notice that this is not what XGBoost is doing; they are simplifying the cost function $L$ by a Taylor Expansion, which leads to a very simple function to compute. They have to compute the Gradient and the Hessian of $L$ with respect to $H_{i-1}$, and they can reuse those number for all potential splits at stage $i$, making the overral computation fast. You can check Loss function Approximation With Taylor Expansion (CrossValidated Q/A) for more details, or the derivation in their paper.

The point is that they have found a way to approximate $L(y,H_{i-1} + h_i)$ efficiently. If you were to evaluate $L$ fully, without insider knowledge allowing optimisation or avoidance or redundant computation, it would take more time per split. It this regard, it is an approximation. However, other gradient boosting implementations also use a proxy cost functions to evaluate the splits, and I do not know whether XGBoost approximation is quicker in this regards than the others.