diff --git a/sklearn/ensemble/_hist_gradient_boosting/grower.py b/sklearn/ensemble/_hist_gradient_boosting/grower.py
index 7eec680082e97..b18c2e5ab2ef0 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/grower.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/grower.py
@@ -266,8 +266,8 @@ def _intilialize_root(self, gradients, hessians, hessians_are_constant):
             self._finalize_leaf(self.root)
             return
 
-        self.root.histograms = self.histogram_builder.compute_histograms_brute(
-            self.root.sample_indices)
+        self.root.histograms, _ = self.histogram_builder.compute_histograms_brute(
+            self.root.sample_indices, None)
         self._compute_best_split_and_push(self.root)
 
     def _compute_best_split_and_push(self, node):
@@ -374,12 +374,10 @@ def split_next(self):
             # smallest number of samples, and the subtraction trick O(n_bins)
             # on the other one.
             tic = time()
-            smallest_child.histograms = \
+            smallest_child.histograms, largest_child.histograms = \
                 self.histogram_builder.compute_histograms_brute(
-                    smallest_child.sample_indices)
-            largest_child.histograms = \
-                self.histogram_builder.compute_histograms_subtraction(
-                    node.histograms, smallest_child.histograms)
+                    smallest_child.sample_indices, node.histograms)
+
             self.total_compute_hist_time += time() - tic
 
             tic = time()
diff --git a/sklearn/ensemble/_hist_gradient_boosting/histogram.pyx b/sklearn/ensemble/_hist_gradient_boosting/histogram.pyx
index c83fa0c79db71..307e81c997dd2 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/histogram.pyx
+++ b/sklearn/ensemble/_hist_gradient_boosting/histogram.pyx
@@ -103,7 +103,9 @@ cdef class HistogramBuilder:
 
     def compute_histograms_brute(
             HistogramBuilder self,
-            const unsigned int [::1] sample_indices):  # IN
+            const unsigned int [::1] sample_indices,
+            hist_struct [:, ::1] parent):  # IN
+
         """Compute the histograms of the node by scanning through all the data.
 
         For a given feature, the complexity is O(n_samples)
@@ -122,6 +124,7 @@ cdef class HistogramBuilder:
             int n_samples
             int feature_idx
             int i
+
             # need local views to avoid python interactions
             unsigned char hessians_are_constant = \
                 self.hessians_are_constant
@@ -130,10 +133,13 @@ cdef class HistogramBuilder:
             G_H_DTYPE_C [::1] gradients = self.gradients
             G_H_DTYPE_C [::1] ordered_hessians = self.ordered_hessians
             G_H_DTYPE_C [::1] hessians = self.hessians
-            hist_struct [:, ::1] histograms = np.zeros(
+            # hist_struct [:, ::1] histograms = np.zeros(
+            hist_struct [:, ::1] histograms = np.empty(
                 shape=(self.n_features, self.max_bins),
-                dtype=HISTOGRAM_DTYPE
-            )
+                dtype=HISTOGRAM_DTYPE)
+            hist_struct [:, ::1] sibling = np.empty(
+                shape=(self.n_features, self.max_bins),
+                dtype=HISTOGRAM_DTYPE)
 
         with nogil:
             n_samples = sample_indices.shape[0]
@@ -150,18 +156,21 @@ cdef class HistogramBuilder:
                         ordered_gradients[i] = gradients[sample_indices[i]]
                         ordered_hessians[i] = hessians[sample_indices[i]]
 
+        with nogil:
+            n_samples = sample_indices.shape[0]
             for feature_idx in prange(n_features, schedule='static'):
                 # Compute histogram of each feature
                 self._compute_histogram_brute_single_feature(
-                    feature_idx, sample_indices, histograms)
-
-        return histograms
+                    feature_idx, sample_indices, histograms, parent, sibling)
+        return histograms, sibling
 
     cdef void _compute_histogram_brute_single_feature(
             HistogramBuilder self,
             const int feature_idx,
             const unsigned int [::1] sample_indices,  # IN
-            hist_struct [:, ::1] histograms) nogil:  # OUT
+            hist_struct [:, ::1] histograms,
+            hist_struct [:, ::1] parent,
+            hist_struct [:, ::1] sibling) nogil:  # OUT
         """Compute the histogram for a given feature."""
 
         cdef:
@@ -176,6 +185,11 @@ cdef class HistogramBuilder:
             unsigned char hessians_are_constant = \
                 self.hessians_are_constant
 
+        for i in range(self.max_bins):
+            histograms[feature_idx, i].sum_gradients = 0
+            histograms[feature_idx, i].sum_hessians = 0
+            histograms[feature_idx, i].count = 0
+
         if root_node:
             if hessians_are_constant:
                 _build_histogram_root_no_hessian(feature_idx, X_binned,
@@ -194,6 +208,7 @@ cdef class HistogramBuilder:
                 _build_histogram(feature_idx, sample_indices,
                                  X_binned, ordered_gradients,
                                  ordered_hessians, histograms)
+            _subtract_histograms(feature_idx, self.max_bins, parent, histograms, sibling)
 
     def compute_histograms_subtraction(
             HistogramBuilder self,
@@ -225,12 +240,14 @@ cdef class HistogramBuilder:
         cdef:
             int feature_idx
             int n_features = self.n_features
-            hist_struct [:, ::1] histograms = np.zeros(
+            # hist_struct [:, ::1] histograms = np.zeros(
+            hist_struct [:, ::1] histograms = np.empty(
                 shape=(self.n_features, self.max_bins),
                 dtype=HISTOGRAM_DTYPE
             )
 
-        for feature_idx in prange(n_features, schedule='static', nogil=True):
+        # for feature_idx in prange(n_features, schedule='static', nogil=True):
+        for feature_idx in range(n_features):
             # Compute histogram of each feature
             _subtract_histograms(feature_idx,
                                  self.max_bins,