A very simple benchmark for binary classification gives:

PR:   768 ms ± 9.76 ms
main: 937 ms ± 128 ms

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
n_samples, n = 10, 10_000
y = np.tile(np.arange(n_samples) % 2, n)
x1 = np.minimum(y, n_samples / 2)
x2 = np.minimum(-y, -n_samples / 2)
X = np.c_[x1, x2]
%timeit GradientBoostingClassifier(n_estimators=100).fit(X, y)

Same for multiclass classification (10 classes) gives:

PR:   9.53 s ± 76.7
main: 28.8 s ± 324 ms

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
n_samples, n = 10, 10_000
y = np.tile(np.arange(n_samples), n)
x1 = np.minimum(y, n_samples / 2)
x2 = np.minimum(-y, -n_samples / 2)
X = np.c_[x1, x2]
%timeit GradientBoostingClassifier(n_estimators=100).fit(X, y)

Given the size of this PR ...

Therefore, I provided the different test_XXX_exact_backward_compat tests.

Dear future 2nd reviewer
This is much easier to review than it looks by just having a look at the tests.

It looks like one of the example broke after merging this PR, namely examples/ensemble/plot_gradient_boosting_regularization.py, because of NaNs.

From build log:

Unexpected failing examples:
/home/circleci/project/examples/ensemble/plot_gradient_boosting_regularization.py failed leaving traceback:
Traceback (most recent call last):
  File "/home/circleci/project/examples/ensemble/plot_gradient_boosting_regularization.py", line 77, in <module>
    test_deviance[i] = 2 * log_loss(y_test, y_proba[:, 1])
  File "/home/circleci/project/sklearn/utils/_param_validation.py", line 211, in wrapper
    return func(*args, **kwargs)
  File "/home/circleci/project/sklearn/metrics/_classification.py", line 2831, in log_loss
    y_pred = check_array(
  File "/home/circleci/project/sklearn/utils/validation.py", line 958, in check_array
    _assert_all_finite(
  File "/home/circleci/project/sklearn/utils/validation.py", line 123, in _assert_all_finite
    _assert_all_finite_element_wise(
  File "/home/circleci/project/sklearn/utils/validation.py", line 172, in _assert_all_finite_element_wise
    raise ValueError(msg_err)
ValueError: Input contains NaN.

@lesteve has been faster to me to report.
I will try to make a fix if I find what is the issue :).

OK so it is due to the _safe_divide where we expect numerator/denominator to raise an error using the errstate. However, it will not work when dividing two scalar np.float64:

In [5]: with np.errstate(divide="raise"):
   ...:     np.float64(0.0) / np.float64(0.0)
   ...: 
<ipython-input-5-941bf0f2d1fa>:2: RuntimeWarning: invalid value encountered in scalar divide
  np.float64(0.0) / np.float64(0.0)

In [6]: with np.errstate(divide="raise"):
   ...:     np.array([1.0, 2.0]) / np.float64(0.0)
   ...: 
---------------------------------------------------------------------------
FloatingPointError                        Traceback (most recent call last)
Cell In[6], line 2
      1 with np.errstate(divide="raise"):
----> 2     np.array([1.0, 2.0]) / np.float64(0.0)

FloatingPointError: divide by zero encountered in divide

Looking at the code, it seems that we expect to have the second case. If we are sure to always have some scalar, then we could convert them to a Python scalar that would always raise an error.

Actually, the only case that go side ways is when both the numerator and denominator goes to 0.0. I assume that it can happen when at the previous iteration we return 0.0 because of the FloatingPointError.