Repository: scikit-learn/scikit-learn
Issue category: Feature Request

Issue:
Feature request: warm starting for histogram-based GBM
#### Description
This is a feature request to add the warm start parameter, which exists for [gradient boosting](https://scikit-learn.org/dev/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html#sklearn.ensemble.GradientBoostingClassifier), to the new [histogram-based gradient boosting](https://scikit-learn.org/dev/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html#sklearn.ensemble.GradientBoostingClassifier).

Rationale: We're using gradient boosting in [Auto-sklearn](https://automl.github.io/auto-sklearn/master/), and fit each model until either the time (given by the user) is up, or the number of trees is reached. For sufficiently large datasets, it is possible that certain configurations time out. Therefore, if possible, we train models iteratively, and only add more trees if time allows. Since the new GBM implementation is really great (not only faster, but also better default accuracy for the problems I tried) we would like to use it within Auto-sklearn as a drop-in, ideally together with iterative training.


Hints:
This is on my TODO list, but I don't know yet when I'll start working on this.

If anyone wants to give it a try I'll be happy to provide review and/or guidance.
@mfeurer thanks for the input!
@NicolasHug I think this would be great to prioritize. Shouldn't be too hard, right?
Honestly I'm happy to work on it but I've been flooding with a lot of PRs recently (still awaiting reviews) and I feel it'd be more helpful to the project if I go back in reviewing mode for a little while.

It shouldn't be too hard for someone to pick it up.

If I'm not the one writing it, we only need one more reviewer instead of 2.
@NicolasHug I agree re PRs but I feel this one would be relatively small compared to the rest? I might even be able to review it ;)
Also: this ties in with successive halving nicely!
If not one picks it up within a week I'll get started then :)
I'll work on this one, though will probably need some guidance along the way @NicolasHug. 
Cool. A good starting point is to look at how this is implemented in `ensemble/gradient_boosting.py`, and start translating it to the new implementation.
@NicolasHug I just realised that this issues is still beyond my current skill level. Looking forward to seeing and learning from the PR though. 
Thanks for letting us know, Matts

I can try to work on this if you did not change your mind @NicolasHug .
Cool go ahead!

Patch:
diff --git a/doc/whats_new/v0.22.rst b/doc/whats_new/v0.22.rst
index d1459bcb8caf3..de0572b9ffd7d 100644
--- a/doc/whats_new/v0.22.rst
+++ b/doc/whats_new/v0.22.rst
@@ -47,6 +47,11 @@ Changelog
   validation data separately to avoid any data leak. :pr:`13933` by
   `NicolasHug`_.
 
+- |Feature| :class:`ensemble.HistGradientBoostingClassifier` and
+  :class:`ensemble.HistGradientBoostingRegressor` have an additional
+  parameter called `warm_start` that enables warm starting. :pr:`14012` by
+  :user:`Johann Faouzi <johannfaouzi>`.
+
 :mod:`sklearn.linear_model`
 ..................
 
diff --git a/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py b/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py
index 466181f445ee8..620e8afe5f276 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py
@@ -26,8 +26,8 @@ class BaseHistGradientBoosting(BaseEstimator, ABC):
     @abstractmethod
     def __init__(self, loss, learning_rate, max_iter, max_leaf_nodes,
                  max_depth, min_samples_leaf, l2_regularization, max_bins,
-                 scoring, validation_fraction, n_iter_no_change, tol, verbose,
-                 random_state):
+                 warm_start, scoring, validation_fraction, n_iter_no_change,
+                 tol, verbose, random_state):
         self.loss = loss
         self.learning_rate = learning_rate
         self.max_iter = max_iter
@@ -36,9 +36,10 @@ def __init__(self, loss, learning_rate, max_iter, max_leaf_nodes,
         self.min_samples_leaf = min_samples_leaf
         self.l2_regularization = l2_regularization
         self.max_bins = max_bins
-        self.n_iter_no_change = n_iter_no_change
-        self.validation_fraction = validation_fraction
+        self.warm_start = warm_start
         self.scoring = scoring
+        self.validation_fraction = validation_fraction
+        self.n_iter_no_change = n_iter_no_change
         self.tol = tol
         self.verbose = verbose
         self.random_state = random_state
@@ -88,7 +89,6 @@ def fit(self, X, y):
         -------
         self : object
         """
-
         fit_start_time = time()
         acc_find_split_time = 0.  # time spent finding the best splits
         acc_apply_split_time = 0.  # time spent splitting nodes
@@ -97,7 +97,13 @@ def fit(self, X, y):
         acc_prediction_time = 0.
         X, y = check_X_y(X, y, dtype=[X_DTYPE])
         y = self._encode_y(y)
-        rng = check_random_state(self.random_state)
+
+        # The rng state must be preserved if warm_start is True
+        if (self.warm_start and hasattr(self, '_rng')):
+            rng = self._rng
+        else:
+            rng = check_random_state(self.random_state)
+            self._rng = rng
 
         self._validate_parameters()
         self.n_features_ = X.shape[1]  # used for validation in predict()
@@ -112,7 +118,6 @@ def fit(self, X, y):
         # data.
         self._in_fit = True
 
-
         self.loss_ = self._get_loss()
 
         self.do_early_stopping_ = (self.n_iter_no_change is not None and
@@ -124,9 +129,15 @@ def fit(self, X, y):
             # stratify for classification
             stratify = y if hasattr(self.loss_, 'predict_proba') else None
 
+            # Save the state of the RNG for the training and validation split.
+            # This is needed in order to have the same split when using
+            # warm starting.
+            if not (self._is_fitted() and self.warm_start):
+                self._train_val_split_seed = rng.randint(1024)
+
             X_train, X_val, y_train, y_val = train_test_split(
                 X, y, test_size=self.validation_fraction, stratify=stratify,
-                random_state=rng)
+                random_state=self._train_val_split_seed)
         else:
             X_train, y_train = X, y
             X_val, y_val = None, None
@@ -142,86 +153,127 @@ def fit(self, X, y):
         if self.verbose:
             print("Fitting gradient boosted rounds:")
 
-        # initialize raw_predictions: those are the accumulated values
-        # predicted by the trees for the training data. raw_predictions has
-        # shape (n_trees_per_iteration, n_samples) where
-        # n_trees_per_iterations is n_classes in multiclass classification,
-        # else 1.
         n_samples = X_binned_train.shape[0]
-        self._baseline_prediction = self.loss_.get_baseline_prediction(
-            y_train, self.n_trees_per_iteration_
-        )
-        raw_predictions = np.zeros(
-            shape=(self.n_trees_per_iteration_, n_samples),
-            dtype=self._baseline_prediction.dtype
-        )
-        raw_predictions += self._baseline_prediction
 
-        # initialize gradients and hessians (empty arrays).
-        # shape = (n_trees_per_iteration, n_samples).
-        gradients, hessians = self.loss_.init_gradients_and_hessians(
-            n_samples=n_samples,
-            prediction_dim=self.n_trees_per_iteration_
-        )
+        # First time calling fit, or no warm start
+        if not (self._is_fitted() and self.warm_start):
+            # Clear random state and score attributes
+            self._clear_state()
+
+            # initialize raw_predictions: those are the accumulated values
+            # predicted by the trees for the training data. raw_predictions has
+            # shape (n_trees_per_iteration, n_samples) where
+            # n_trees_per_iterations is n_classes in multiclass classification,
+            # else 1.
+            self._baseline_prediction = self.loss_.get_baseline_prediction(
+                y_train, self.n_trees_per_iteration_
+            )
+            raw_predictions = np.zeros(
+                shape=(self.n_trees_per_iteration_, n_samples),
+                dtype=self._baseline_prediction.dtype
+            )
+            raw_predictions += self._baseline_prediction
 
-        # predictors is a matrix (list of lists) of TreePredictor objects
-        # with shape (n_iter_, n_trees_per_iteration)
-        self._predictors = predictors = []
+            # initialize gradients and hessians (empty arrays).
+            # shape = (n_trees_per_iteration, n_samples).
+            gradients, hessians = self.loss_.init_gradients_and_hessians(
+                n_samples=n_samples,
+                prediction_dim=self.n_trees_per_iteration_
+            )
 
-        # Initialize structures and attributes related to early stopping
-        self.scorer_ = None  # set if scoring != loss
-        raw_predictions_val = None  # set if scoring == loss and use val
-        self.train_score_ = []
-        self.validation_score_ = []
-        if self.do_early_stopping_:
-            # populate train_score and validation_score with the predictions
-            # of the initial model (before the first tree)
+            # predictors is a matrix (list of lists) of TreePredictor objects
+            # with shape (n_iter_, n_trees_per_iteration)
+            self._predictors = predictors = []
 
-            if self.scoring == 'loss':
-                # we're going to compute scoring w.r.t the loss. As losses
-                # take raw predictions as input (unlike the scorers), we can
-                # optimize a bit and avoid repeating computing the predictions
-                # of the previous trees. We'll re-use raw_predictions (as it's
-                # needed for training anyway) for evaluating the training
-                # loss, and create raw_predictions_val for storing the
-                # raw predictions of the validation data.
-
-                if self._use_validation_data:
-                    raw_predictions_val = np.zeros(
-                        shape=(self.n_trees_per_iteration_,
-                               X_binned_val.shape[0]),
-                        dtype=self._baseline_prediction.dtype
-                    )
+            # Initialize structures and attributes related to early stopping
+            self.scorer_ = None  # set if scoring != loss
+            raw_predictions_val = None  # set if scoring == loss and use val
+            self.train_score_ = []
+            self.validation_score_ = []
+
+            if self.do_early_stopping_:
+                # populate train_score and validation_score with the
+                # predictions of the initial model (before the first tree)
 
-                    raw_predictions_val += self._baseline_prediction
+                if self.scoring == 'loss':
+                    # we're going to compute scoring w.r.t the loss. As losses
+                    # take raw predictions as input (unlike the scorers), we
+                    # can optimize a bit and avoid repeating computing the
+                    # predictions of the previous trees. We'll re-use
+                    # raw_predictions (as it's needed for training anyway) for
+                    # evaluating the training loss, and create
+                    # raw_predictions_val for storing the raw predictions of
+                    # the validation data.
 
-                self._check_early_stopping_loss(raw_predictions, y_train,
-                                                raw_predictions_val, y_val)
-            else:
-                self.scorer_ = check_scoring(self, self.scoring)
-                # scorer_ is a callable with signature (est, X, y) and calls
-                # est.predict() or est.predict_proba() depending on its nature.
-                # Unfortunately, each call to scorer_() will compute
-                # the predictions of all the trees. So we use a subset of the
-                # training set to compute train scores.
-                subsample_size = 10000  # should we expose this parameter?
-                indices = np.arange(X_binned_train.shape[0])
-                if X_binned_train.shape[0] > subsample_size:
-                    # TODO: not critical but stratify using resample()
-                    indices = rng.choice(indices, subsample_size,
-                                         replace=False)
-                X_binned_small_train = X_binned_train[indices]
-                y_small_train = y_train[indices]
-                # Predicting is faster on C-contiguous arrays.
-                X_binned_small_train = np.ascontiguousarray(
-                    X_binned_small_train)
-
-                self._check_early_stopping_scorer(
-                    X_binned_small_train, y_small_train,
-                    X_binned_val, y_val,
+                    if self._use_validation_data:
+                        raw_predictions_val = np.zeros(
+                            shape=(self.n_trees_per_iteration_,
+                                   X_binned_val.shape[0]),
+                            dtype=self._baseline_prediction.dtype
+                        )
+
+                        raw_predictions_val += self._baseline_prediction
+
+                    self._check_early_stopping_loss(raw_predictions, y_train,
+                                                    raw_predictions_val, y_val)
+                else:
+                    self.scorer_ = check_scoring(self, self.scoring)
+                    # scorer_ is a callable with signature (est, X, y) and
+                    # calls est.predict() or est.predict_proba() depending on
+                    # its nature.
+                    # Unfortunately, each call to scorer_() will compute
+                    # the predictions of all the trees. So we use a subset of
+                    # the training set to compute train scores.
+
+                    # Save the seed for the small trainset generator
+                    self._small_trainset_seed = rng.randint(1024)
+
+                    # Compute the subsample set
+                    (X_binned_small_train,
+                     y_small_train) = self._get_small_trainset(
+                        X_binned_train, y_train, self._small_trainset_seed)
+
+                    self._check_early_stopping_scorer(
+                        X_binned_small_train, y_small_train,
+                        X_binned_val, y_val,
+                    )
+            begin_at_stage = 0
+
+        # warm start: this is not the first time fit was called
+        else:
+            # Check that the maximum number of iterations is not smaller
+            # than the number of iterations from the previous fit
+            if self.max_iter < self.n_iter_:
+                raise ValueError(
+                    'max_iter=%d must be larger than or equal to '
+                    'n_iter_=%d when warm_start==True'
+                    % (self.max_iter, self.n_iter_)
                 )
 
-        for iteration in range(self.max_iter):
+            # Convert array attributes to lists
+            self.train_score_ = self.train_score_.tolist()
+            self.validation_score_ = self.validation_score_.tolist()
+
+            # Compute raw predictions
+            raw_predictions = self._raw_predict(X_binned_train)
+
+            if self.do_early_stopping_ and self.scoring != 'loss':
+                # Compute the subsample set
+                X_binned_small_train, y_small_train = self._get_small_trainset(
+                    X_binned_train, y_train, self._small_trainset_seed)
+
+            # Initialize the gradients and hessians
+            gradients, hessians = self.loss_.init_gradients_and_hessians(
+                n_samples=n_samples,
+                prediction_dim=self.n_trees_per_iteration_
+            )
+
+            # Get the predictors from the previous fit
+            predictors = self._predictors
+
+            begin_at_stage = self.n_iter_
+
+        for iteration in range(begin_at_stage, self.max_iter):
 
             if self.verbose:
                 iteration_start_time = time()
@@ -318,13 +370,38 @@ def fit(self, X, y):
         del self._in_fit  # hard delete so we're sure it can't be used anymore
         return self
 
+    def _is_fitted(self):
+        return len(getattr(self, '_predictors', [])) > 0
+
+    def _clear_state(self):
+        """Clear the state of the gradient boosting model."""
+        for var in ('train_score_', 'validation_score_', '_rng'):
+            if hasattr(self, var):
+                delattr(self, var)
+
+    def _get_small_trainset(self, X_binned_train, y_train, seed):
+        """Compute the indices of the subsample set and return this set.
+
+        For efficiency, we need to subsample the training set to compute scores
+        with scorers.
+        """
+        subsample_size = 10000
+        rng = check_random_state(seed)
+        indices = np.arange(X_binned_train.shape[0])
+        if X_binned_train.shape[0] > subsample_size:
+            # TODO: not critical but stratify using resample()
+            indices = rng.choice(indices, subsample_size, replace=False)
+        X_binned_small_train = X_binned_train[indices]
+        y_small_train = y_train[indices]
+        X_binned_small_train = np.ascontiguousarray(X_binned_small_train)
+        return X_binned_small_train, y_small_train
+
     def _check_early_stopping_scorer(self, X_binned_small_train, y_small_train,
                                      X_binned_val, y_val):
         """Check if fitting should be early-stopped based on scorer.
 
         Scores are computed on validation data or on training data.
         """
-
         self.train_score_.append(
             self.scorer_(self, X_binned_small_train, y_small_train)
         )
@@ -555,6 +632,11 @@ class HistGradientBoostingRegressor(BaseHistGradientBoosting, RegressorMixin):
         allows for a much faster training stage. Features with a small
         number of unique values may use less than ``max_bins`` bins. Must be no
         larger than 256.
+    warm_start : bool, optional (default=False)
+        When set to ``True``, reuse the solution of the previous call to fit
+        and add more estimators to the ensemble. For results to be valid, the
+        estimator should be re-trained on the same data only.
+        See :term:`the Glossary <warm_start>`.
     scoring : str or callable or None, optional (default=None)
         Scoring parameter to use for early stopping. It can be a single
         string (see :ref:`scoring_parameter`) or a callable (see
@@ -568,7 +650,7 @@ class HistGradientBoostingRegressor(BaseHistGradientBoosting, RegressorMixin):
     n_iter_no_change : int or None, optional (default=None)
         Used to determine when to "early stop". The fitting process is
         stopped when none of the last ``n_iter_no_change`` scores are better
-        than the ``n_iter_no_change - 1``th-to-last one, up to some
+        than the ``n_iter_no_change - 1`` -th-to-last one, up to some
         tolerance. If None or 0, no early-stopping is done.
     tol : float or None, optional (default=1e-7)
         The absolute tolerance to use when comparing scores during early
@@ -592,13 +674,13 @@ class HistGradientBoostingRegressor(BaseHistGradientBoosting, RegressorMixin):
     n_trees_per_iteration_ : int
         The number of tree that are built at each iteration. For regressors,
         this is always 1.
-    train_score_ : ndarray, shape (max_iter + 1,)
+    train_score_ : ndarray, shape (n_iter_ + 1,)
         The scores at each iteration on the training data. The first entry
         is the score of the ensemble before the first iteration. Scores are
         computed according to the ``scoring`` parameter. If ``scoring`` is
         not 'loss', scores are computed on a subset of at most 10 000
         samples. Empty if no early stopping.
-    validation_score_ : ndarray, shape (max_iter + 1,)
+    validation_score_ : ndarray, shape (n_iter_ + 1,)
         The scores at each iteration on the held-out validation data. The
         first entry is the score of the ensemble before the first iteration.
         Scores are computed according to the ``scoring`` parameter. Empty if
@@ -621,14 +703,16 @@ class HistGradientBoostingRegressor(BaseHistGradientBoosting, RegressorMixin):
     def __init__(self, loss='least_squares', learning_rate=0.1,
                  max_iter=100, max_leaf_nodes=31, max_depth=None,
                  min_samples_leaf=20, l2_regularization=0., max_bins=256,
-                 scoring=None, validation_fraction=0.1, n_iter_no_change=None,
-                 tol=1e-7, verbose=0, random_state=None):
+                 warm_start=False, scoring=None, validation_fraction=0.1,
+                 n_iter_no_change=None, tol=1e-7, verbose=0,
+                 random_state=None):
         super(HistGradientBoostingRegressor, self).__init__(
             loss=loss, learning_rate=learning_rate, max_iter=max_iter,
             max_leaf_nodes=max_leaf_nodes, max_depth=max_depth,
             min_samples_leaf=min_samples_leaf,
             l2_regularization=l2_regularization, max_bins=max_bins,
-            scoring=scoring, validation_fraction=validation_fraction,
+            warm_start=warm_start, scoring=scoring,
+            validation_fraction=validation_fraction,
             n_iter_no_change=n_iter_no_change, tol=tol, verbose=verbose,
             random_state=random_state)
 
@@ -723,6 +807,11 @@ class HistGradientBoostingClassifier(BaseHistGradientBoosting,
         allows for a much faster training stage. Features with a small
         number of unique values may use less than ``max_bins`` bins. Must be no
         larger than 256.
+    warm_start : bool, optional (default=False)
+        When set to ``True``, reuse the solution of the previous call to fit
+        and add more estimators to the ensemble. For results to be valid, the
+        estimator should be re-trained on the same data only.
+        See :term:`the Glossary <warm_start>`.
     scoring : str or callable or None, optional (default=None)
         Scoring parameter to use for early stopping. It can be a single
         string (see :ref:`scoring_parameter`) or a callable (see
@@ -736,7 +825,7 @@ class HistGradientBoostingClassifier(BaseHistGradientBoosting,
     n_iter_no_change : int or None, optional (default=None)
         Used to determine when to "early stop". The fitting process is
         stopped when none of the last ``n_iter_no_change`` scores are better
-        than the ``n_iter_no_change - 1``th-to-last one, up to some
+        than the ``n_iter_no_change - 1`` -th-to-last one, up to some
         tolerance. If None or 0, no early-stopping is done.
     tol : float or None, optional (default=1e-7)
         The absolute tolerance to use when comparing scores. The higher the
@@ -761,13 +850,13 @@ class HistGradientBoostingClassifier(BaseHistGradientBoosting,
         The number of tree that are built at each iteration. This is equal to 1
         for binary classification, and to ``n_classes`` for multiclass
         classification.
-    train_score_ : ndarray, shape (max_iter + 1,)
+    train_score_ : ndarray, shape (n_iter_ + 1,)
         The scores at each iteration on the training data. The first entry
         is the score of the ensemble before the first iteration. Scores are
         computed according to the ``scoring`` parameter. If ``scoring`` is
         not 'loss', scores are computed on a subset of at most 10 000
         samples. Empty if no early stopping.
-    validation_score_ : ndarray, shape (max_iter + 1,)
+    validation_score_ : ndarray, shape (n_iter_ + 1,)
         The scores at each iteration on the held-out validation data. The
         first entry is the score of the ensemble before the first iteration.
         Scores are computed according to the ``scoring`` parameter. Empty if
@@ -790,15 +879,16 @@ class HistGradientBoostingClassifier(BaseHistGradientBoosting,
 
     def __init__(self, loss='auto', learning_rate=0.1, max_iter=100,
                  max_leaf_nodes=31, max_depth=None, min_samples_leaf=20,
-                 l2_regularization=0., max_bins=256, scoring=None,
-                 validation_fraction=0.1, n_iter_no_change=None, tol=1e-7,
-                 verbose=0, random_state=None):
+                 l2_regularization=0., max_bins=256, warm_start=False,
+                 scoring=None, validation_fraction=0.1, n_iter_no_change=None,
+                 tol=1e-7, verbose=0, random_state=None):
         super(HistGradientBoostingClassifier, self).__init__(
             loss=loss, learning_rate=learning_rate, max_iter=max_iter,
             max_leaf_nodes=max_leaf_nodes, max_depth=max_depth,
             min_samples_leaf=min_samples_leaf,
             l2_regularization=l2_regularization, max_bins=max_bins,
-            scoring=scoring, validation_fraction=validation_fraction,
+            warm_start=warm_start, scoring=scoring,
+            validation_fraction=validation_fraction,
             n_iter_no_change=n_iter_no_change, tol=tol, verbose=verbose,
             random_state=random_state)
 


Test patch:
diff --git a/sklearn/ensemble/_hist_gradient_boosting/tests/test_warm_start.py b/sklearn/ensemble/_hist_gradient_boosting/tests/test_warm_start.py
new file mode 100644
index 0000000000000..806ad94ccee98
--- /dev/null
+++ b/sklearn/ensemble/_hist_gradient_boosting/tests/test_warm_start.py
@@ -0,0 +1,190 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+from numpy.testing import assert_allclose
+
+import pytest
+
+from sklearn.base import clone
+from sklearn.datasets import make_classification, make_regression
+
+# To use this experimental feature, we need to explicitly ask for it:
+from sklearn.experimental import enable_hist_gradient_boosting  # noqa
+from sklearn.ensemble import HistGradientBoostingRegressor
+from sklearn.ensemble import HistGradientBoostingClassifier
+
+
+X_classification, y_classification = make_classification(random_state=0)
+X_regression, y_regression = make_regression(random_state=0)
+
+
+def _assert_predictor_equal(gb_1, gb_2, X):
+    """Assert that two HistGBM instances are identical."""
+    # Check identical nodes for each tree
+    for (pred_ith_1, pred_ith_2) in zip(gb_1._predictors, gb_2._predictors):
+        for (predictor_1, predictor_2) in zip(pred_ith_1, pred_ith_2):
+            assert_array_equal(predictor_1.nodes, predictor_2.nodes)
+
+    # Check identical predictions
+    assert_allclose(gb_1.predict(X), gb_2.predict(X))
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+def test_max_iter_with_warm_start_validation(GradientBoosting, X, y):
+    # Check that a ValueError is raised when the maximum number of iterations
+    # is smaller than the number of iterations from the previous fit when warm
+    # start is True.
+
+    estimator = GradientBoosting(max_iter=50, warm_start=True)
+    estimator.fit(X, y)
+    estimator.set_params(max_iter=25)
+    err_msg = ('max_iter=25 must be larger than or equal to n_iter_=50 '
+               'when warm_start==True')
+    with pytest.raises(ValueError, match=err_msg):
+        estimator.fit(X, y)
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+def test_warm_start_yields_identical_results(GradientBoosting, X, y):
+    # Make sure that fitting 50 iterations and then 25 with warm start is
+    # equivalent to fitting 75 iterations.
+
+    rng = 42
+    gb_warm_start = GradientBoosting(
+        n_iter_no_change=100, max_iter=50, random_state=rng, warm_start=True
+    )
+    gb_warm_start.fit(X, y).set_params(max_iter=75).fit(X, y)
+
+    gb_no_warm_start = GradientBoosting(
+        n_iter_no_change=100, max_iter=75, random_state=rng, warm_start=False
+    )
+    gb_no_warm_start.fit(X, y)
+
+    # Check that both predictors are equal
+    _assert_predictor_equal(gb_warm_start, gb_no_warm_start, X)
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+def test_warm_start_max_depth(GradientBoosting, X, y):
+    # Test if possible to fit trees of different depth in ensemble.
+    gb = GradientBoosting(max_iter=100, min_samples_leaf=1,
+                          warm_start=True, max_depth=2)
+    gb.fit(X, y)
+    gb.set_params(max_iter=110, max_depth=3)
+    gb.fit(X, y)
+
+    # First 100 trees have max_depth == 2
+    for i in range(100):
+        assert gb._predictors[i][0].get_max_depth() == 2
+    # Last 10 trees have max_depth == 3
+    for i in range(1, 11):
+        assert gb._predictors[-i][0].get_max_depth() == 3
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+def test_warm_start_early_stopping(GradientBoosting, X, y):
+    # Make sure that early stopping occurs after a small number of iterations
+    # when fitting a second time with warm starting.
+
+    n_iter_no_change = 5
+    gb = GradientBoosting(
+        n_iter_no_change=n_iter_no_change, max_iter=10000,
+        random_state=42, warm_start=True, tol=1e-3
+    )
+    gb.fit(X, y)
+    n_iter_first_fit = gb.n_iter_
+    gb.fit(X, y)
+    n_iter_second_fit = gb.n_iter_
+    assert n_iter_second_fit - n_iter_first_fit < n_iter_no_change
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+def test_warm_start_equal_n_estimators(GradientBoosting, X, y):
+    # Test if warm start with equal n_estimators does nothing
+    gb_1 = GradientBoosting(max_depth=2)
+    gb_1.fit(X, y)
+
+    gb_2 = clone(gb_1)
+    gb_2.set_params(max_iter=gb_1.max_iter, warm_start=True)
+    gb_2.fit(X, y)
+
+    # Check that both predictors are equal
+    _assert_predictor_equal(gb_1, gb_2, X)
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+def test_warm_start_clear(GradientBoosting, X, y):
+    # Test if fit clears state.
+    gb_1 = GradientBoosting(n_iter_no_change=5, random_state=42)
+    gb_1.fit(X, y)
+
+    gb_2 = GradientBoosting(n_iter_no_change=5, random_state=42,
+                            warm_start=True)
+    gb_2.fit(X, y)  # inits state
+    gb_2.set_params(warm_start=False)
+    gb_2.fit(X, y)  # clears old state and equals est
+
+    # Check that both predictors have the same train_score_ and
+    # validation_score_ attributes
+    assert_allclose(gb_1.train_score_, gb_2.train_score_)
+    assert_allclose(gb_1.validation_score_, gb_2.validation_score_)
+
+    # Check that both predictors are equal
+    _assert_predictor_equal(gb_1, gb_2, X)
+
+
+@pytest.mark.parametrize('GradientBoosting, X, y', [
+    (HistGradientBoostingClassifier, X_classification, y_classification),
+    (HistGradientBoostingRegressor, X_regression, y_regression)
+])
+@pytest.mark.parametrize('rng_type', ('int', 'instance'))
+def test_random_seeds_warm_start(GradientBoosting, X, y, rng_type):
+    # Make sure the seeds for train/val split and small trainset subsampling
+    # are correctly set in a warm start context.
+    def _get_rng(rng_type):
+        # Helper to avoid consuming rngs
+        if rng_type == 'int':
+            return 42
+        else:
+            return np.random.RandomState(0)
+
+    random_state = _get_rng(rng_type)
+    gb_1 = GradientBoosting(n_iter_no_change=5, max_iter=2,
+                            random_state=random_state)
+    gb_1.fit(X, y)
+    train_val_seed_1 = gb_1._train_val_split_seed
+    small_trainset_seed_1 = gb_1._small_trainset_seed
+
+    random_state = _get_rng(rng_type)
+    gb_2 = GradientBoosting(n_iter_no_change=5, max_iter=2,
+                            random_state=random_state, warm_start=True)
+    gb_2.fit(X, y)  # inits state
+    train_val_seed_2 = gb_2._train_val_split_seed
+    small_trainset_seed_2 = gb_2._small_trainset_seed
+    gb_2.fit(X, y)  # clears old state and equals est
+    train_val_seed_3 = gb_2._train_val_split_seed
+    small_trainset_seed_3 = gb_2._small_trainset_seed
+
+    # Check that all seeds are equal
+    assert train_val_seed_1 == train_val_seed_2
+    assert small_trainset_seed_1 == small_trainset_seed_2
+
+    assert train_val_seed_2 == train_val_seed_3
+    assert small_trainset_seed_2 == small_trainset_seed_3