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

HistGradientBoosting* interaction constraints 
#### Describe the workflow you want to enable
I'd like to use `HistGradientBoostingClassifier` and `HistGradientBoostingRegressor` with the possibility to set interaction constraints for certain features. As said in https://github.com/microsoft/LightGBM/issues/2884#issue-577348269, it is one way to make those *black boxes* more intuitive and interpretable. In addition, it makes it much more easy to marginalize over those features.

#### Additional context
LightGBM has `interaction_constraints`, see [their docs](https://lightgbm.readthedocs.io/en/latest/Parameters.html#interaction_constraints). XGBoost has them, see [their docs](https://xgboost.readthedocs.io/en/latest/parameter.html#parameters-for-tree-booster).
Have also a look at the [XGBoost tutorial on interaction constraints](https://xgboost.readthedocs.io/en/latest/tutorials/feature_interaction_constraint.html) for a nice visualization and for potential benefits:

> - Better predictive performance from focusing on interactions that work – whether through domain specific knowledge or algorithms that rank interactions
> - Less noise in predictions; better generalization
> - More control to the user on what the model can fit. For example, the user may want to exclude some interactions even if they perform well due to regulatory constraints



Functions to edit:
sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:BaseHistGradientBoosting
sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:BaseHistGradientBoosting.__init__
sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:BaseHistGradientBoosting.fit
sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:HistGradientBoostingRegressor.__init__
sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:HistGradientBoostingClassifier.__init__
sklearn/ensemble/_hist_gradient_boosting/grower.py:TreeNode.__init__
sklearn/ensemble/_hist_gradient_boosting/grower.py:TreeGrower.__init__
sklearn/ensemble/_hist_gradient_boosting/grower.py:TreeGrower._intilialize_root
sklearn/ensemble/_hist_gradient_boosting/grower.py:TreeGrower._compute_best_split_and_push
sklearn/ensemble/_hist_gradient_boosting/grower.py:TreeGrower.split_next
Functions added: sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:BaseHistGradientBoosting._check_interaction_cst
sklearn/ensemble/_hist_gradient_boosting/grower.py:TreeGrower._compute_interactions

Patch:
diff --git a/doc/modules/ensemble.rst b/doc/modules/ensemble.rst
index 3aad30d792e2a..0c4159165e181 100644
--- a/doc/modules/ensemble.rst
+++ b/doc/modules/ensemble.rst
@@ -317,9 +317,9 @@ to the prediction function.
 
 .. topic:: References
 
- .. [L2014] G. Louppe,
-         "Understanding Random Forests: From Theory to Practice",
-         PhD Thesis, U. of Liege, 2014.
+ .. [L2014] G. Louppe, :arxiv:`"Understanding Random Forests: From Theory to
+    Practice" <1407.7502>`,
+    PhD Thesis, U. of Liege, 2014.
 
 .. _random_trees_embedding:
 
@@ -711,7 +711,7 @@ space.
   accurate enough: the tree can only output integer values. As a result, the
   leaves values of the tree :math:`h_m` are modified once the tree is
   fitted, such that the leaves values minimize the loss :math:`L_m`. The
-  update is loss-dependent: for the absolute error loss, the value of 
+  update is loss-dependent: for the absolute error loss, the value of
   a leaf is updated to the median of the samples in that leaf.
 
 Classification
@@ -1174,6 +1174,44 @@ Also, monotonic constraints are not supported for multiclass classification.
 
   * :ref:`sphx_glr_auto_examples_ensemble_plot_monotonic_constraints.py`
 
+.. _interaction_cst_hgbt:
+
+Interaction constraints
+-----------------------
+
+A priori, the histogram gradient boosting trees are allowed to use any feature
+to split a node into child nodes. This creates so called interactions between
+features, i.e. usage of different features as split along a branch. Sometimes,
+one wants to restrict the possible interactions, see [Mayer2022]_. This can be
+done by the parameter ``interaction_cst``, where one can specify the indices
+of features that are allowed to interact.
+For instance, with 3 features in total, ``interaction_cst=[{0}, {1}, {2}]``
+forbids all interactions.
+The constraints ``[{0, 1}, {1, 2}]`` specifies two groups of possibly
+interacting features. Features 0 and 1 may interact with each other, as well
+as features 1 and 2. But note that features 0 and 2 are forbidden to interact.
+The following depicts a tree and the possible splits of the tree:
+
+.. code-block:: none
+
+      1      <- Both constraint groups could be applied from now on
+     / \
+    1   2    <- Left split still fulfills both constraint groups.
+   / \ / \      Right split at feature 2 has only group {1, 2} from now on.
+
+LightGBM uses the same logic for overlapping groups.
+
+Note that features not listed in ``interaction_cst`` are automatically
+assigned an interaction group for themselves. With again 3 features, this
+means that ``[{0}]`` is equivalent to ``[{0}, {1, 2}]``.
+
+.. topic:: References
+
+  .. [Mayer2022] M. Mayer, S.C. Bourassa, M. Hoesli, and D.F. Scognamiglio.
+     2022. :doi:`Machine Learning Applications to Land and Structure Valuation
+     <10.3390/jrfm15050193>`.
+     Journal of Risk and Financial Management 15, no. 5: 193
+
 Low-level parallelism
 ---------------------
 
diff --git a/doc/whats_new/v1.2.rst b/doc/whats_new/v1.2.rst
index 7734e771c0c1e..0de25ff7f45e6 100644
--- a/doc/whats_new/v1.2.rst
+++ b/doc/whats_new/v1.2.rst
@@ -242,6 +242,12 @@ Changelog
 :mod:`sklearn.ensemble`
 .......................
 
+- |Feature| :class:`~sklearn.ensemble.HistGradientBoostingClassifier` and
+  :class:`~sklearn.ensemble.HistGradientBoostingRegressor` now support
+  interaction constraints via the argument `interaction_cst` of their
+  constructors.
+  :pr:`21020` by :user:`Christian Lorentzen <lorentzenchr>`.
+
 - |Feature| Adds `class_weight` to :class:`ensemble.HistGradientBoostingClassifier`.
   :pr:`22014` by `Thomas Fan`_.
 
diff --git a/examples/inspection/plot_partial_dependence.py b/examples/inspection/plot_partial_dependence.py
index d7480a7898424..a7ef29edef183 100644
--- a/examples/inspection/plot_partial_dependence.py
+++ b/examples/inspection/plot_partial_dependence.py
@@ -255,6 +255,85 @@
 # house age, whereas for values less than two there is a strong dependence on
 # age.
 #
+# Interaction constraints
+# .......................
+#
+# The histogram gradient boosters have an interesting option to constrain
+# possible interactions among features. In the following, we do not allow any
+# interactions and thus render the model as a version of a tree-based boosted
+# generalized additive model (GAM). This makes the model more interpretable
+# as the effect of each feature can be investigated independently of all others.
+#
+# We train the :class:`~sklearn.ensemble.HistGradientBoostingRegressor` again,
+# now with `interaction_cst`, where we  pass for each feature a list containing
+# only its own index, e.g. `[[0], [1], [2], ..]`.
+
+print("Training interaction constraint HistGradientBoostingRegressor...")
+tic = time()
+est_no_interactions = HistGradientBoostingRegressor(
+    interaction_cst=[[i] for i in range(X_train.shape[1])]
+)
+est_no_interactions.fit(X_train, y_train)
+print(f"done in {time() - tic:.3f}s")
+
+# %%
+# The easiest way to show the effect of forbidden interactions is again the
+# ICE plots.
+
+print("Computing partial dependence plots...")
+tic = time()
+display = PartialDependenceDisplay.from_estimator(
+    est_no_interactions,
+    X_train,
+    ["MedInc", "AveOccup", "HouseAge", "AveRooms"],
+    kind="both",
+    subsample=50,
+    n_jobs=3,
+    grid_resolution=20,
+    random_state=0,
+    ice_lines_kw={"color": "tab:blue", "alpha": 0.2, "linewidth": 0.5},
+    pd_line_kw={"color": "tab:orange", "linestyle": "--"},
+)
+
+print(f"done in {time() - tic:.3f}s")
+display.figure_.suptitle(
+    "Partial dependence of house value with Gradient Boosting\n"
+    "and no interactions allowed"
+)
+display.figure_.subplots_adjust(wspace=0.4, hspace=0.3)
+
+# %%
+# All 4 plots have parallel ICE lines meaning there is no interaction in the
+# model.
+# Let us also have a look at the corresponding 2D-plot.
+
+print("Computing partial dependence plots...")
+tic = time()
+_, ax = plt.subplots(ncols=3, figsize=(9, 4))
+display = PartialDependenceDisplay.from_estimator(
+    est_no_interactions,
+    X_train,
+    ["AveOccup", "HouseAge", ("AveOccup", "HouseAge")],
+    kind="average",
+    n_jobs=3,
+    grid_resolution=20,
+    ax=ax,
+)
+print(f"done in {time() - tic:.3f}s")
+display.figure_.suptitle(
+    "Partial dependence of house value with Gradient Boosting\n"
+    "and no interactions allowed"
+)
+display.figure_.subplots_adjust(wspace=0.4, hspace=0.3)
+
+# %%
+# Although the 2D-plot shows much less interaction compared with the 2D-plot
+# from above, it is much harder to come to the conclusion that there is no
+# interaction at all. This might be a cause of the discrete predictions of
+# trees in combination with numerically precision of partial dependence.
+# We also observe that the univariate dependence plots have slightly changed
+# as the model tries to compensate for the forbidden interactions.
+#
 # 3D interaction plots
 # --------------------
 #
diff --git a/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py b/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py
index 68931e50b16ca..c5fe46496a577 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py
@@ -2,6 +2,7 @@
 # Author: Nicolas Hug
 
 from abc import ABC, abstractmethod
+from collections.abc import Iterable
 from functools import partial
 from numbers import Real, Integral
 import warnings
@@ -91,6 +92,7 @@ class BaseHistGradientBoosting(BaseEstimator, ABC):
         "min_samples_leaf": [Interval(Integral, 1, None, closed="left")],
         "l2_regularization": [Interval(Real, 0, None, closed="left")],
         "monotonic_cst": ["array-like", None],
+        "interaction_cst": [Iterable, None],
         "n_iter_no_change": [Interval(Integral, 1, None, closed="left")],
         "validation_fraction": [
             Interval(Real, 0, 1, closed="neither"),
@@ -121,6 +123,7 @@ def __init__(
         max_bins,
         categorical_features,
         monotonic_cst,
+        interaction_cst,
         warm_start,
         early_stopping,
         scoring,
@@ -139,6 +142,7 @@ def __init__(
         self.l2_regularization = l2_regularization
         self.max_bins = max_bins
         self.monotonic_cst = monotonic_cst
+        self.interaction_cst = interaction_cst
         self.categorical_features = categorical_features
         self.warm_start = warm_start
         self.early_stopping = early_stopping
@@ -252,6 +256,42 @@ def _check_categories(self, X):
 
         return is_categorical, known_categories
 
+    def _check_interaction_cst(self, n_features):
+        """Check and validation for interaction constraints."""
+        if self.interaction_cst is None:
+            return None
+
+        if not (
+            isinstance(self.interaction_cst, Iterable)
+            and all(isinstance(x, Iterable) for x in self.interaction_cst)
+        ):
+            raise ValueError(
+                "Interaction constraints must be None or an iterable of iterables, "
+                f"got: {self.interaction_cst!r}."
+            )
+
+        invalid_indices = [
+            x
+            for cst_set in self.interaction_cst
+            for x in cst_set
+            if not (isinstance(x, Integral) and 0 <= x < n_features)
+        ]
+        if invalid_indices:
+            raise ValueError(
+                "Interaction constraints must consist of integer indices in [0,"
+                f" n_features - 1] = [0, {n_features - 1}], specifying the position of"
+                f" features, got invalid indices: {invalid_indices!r}"
+            )
+
+        constraints = [set(group) for group in self.interaction_cst]
+
+        # Add all not listed features as own group by default.
+        rest = set(range(n_features)) - set().union(*constraints)
+        if len(rest) > 0:
+            constraints.append(rest)
+
+        return constraints
+
     def fit(self, X, y, sample_weight=None):
         """Fit the gradient boosting model.
 
@@ -308,6 +348,9 @@ def fit(self, X, y, sample_weight=None):
 
         self.is_categorical_, known_categories = self._check_categories(X)
 
+        # Encode constraints into a list of sets of features indices (integers).
+        interaction_cst = self._check_interaction_cst(self._n_features)
+
         # we need this stateful variable to tell raw_predict() that it was
         # called from fit() (this current method), and that the data it has
         # received is pre-binned.
@@ -595,6 +638,7 @@ def fit(self, X, y, sample_weight=None):
                     has_missing_values=has_missing_values,
                     is_categorical=self.is_categorical_,
                     monotonic_cst=self.monotonic_cst,
+                    interaction_cst=interaction_cst,
                     max_leaf_nodes=self.max_leaf_nodes,
                     max_depth=self.max_depth,
                     min_samples_leaf=self.min_samples_leaf,
@@ -1191,6 +1235,22 @@ class HistGradientBoostingRegressor(RegressorMixin, BaseHistGradientBoosting):
 
         .. versionadded:: 0.23
 
+    interaction_cst : iterable of iterables of int, default=None
+        Specify interaction constraints, i.e. sets of features which can
+        only interact with each other in child nodes splits.
+
+        Each iterable materializes a constraint by the set of indices of
+        the features that are allowed to interact with each other.
+        If there are more features than specified in these constraints,
+        they are treated as if they were specified as an additional set.
+
+        For instance, with 5 features in total, `interaction_cst=[{0, 1}]`
+        is equivalent to `interaction_cst=[{0, 1}, {2, 3, 4}]`,
+        and specifies that each branch of a tree will either only split
+        on features 0 and 1 or only split on features 2, 3 and 4.
+
+        .. versionadded:: 1.2
+
     warm_start : bool, 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
@@ -1315,6 +1375,7 @@ def __init__(
         max_bins=255,
         categorical_features=None,
         monotonic_cst=None,
+        interaction_cst=None,
         warm_start=False,
         early_stopping="auto",
         scoring="loss",
@@ -1334,6 +1395,7 @@ def __init__(
             l2_regularization=l2_regularization,
             max_bins=max_bins,
             monotonic_cst=monotonic_cst,
+            interaction_cst=interaction_cst,
             categorical_features=categorical_features,
             early_stopping=early_stopping,
             warm_start=warm_start,
@@ -1505,6 +1567,22 @@ class HistGradientBoostingClassifier(ClassifierMixin, BaseHistGradientBoosting):
 
         .. versionadded:: 0.23
 
+    interaction_cst : iterable of iterables of int, default=None
+        Specify interaction constraints, i.e. sets of features which can
+        only interact with each other in child nodes splits.
+
+        Each iterable materializes a constraint by the set of indices of
+        the features that are allowed to interact with each other.
+        If there are more features than specified in these constraints,
+        they are treated as if they were specified as an additional set.
+
+        For instance, with 5 features in total, `interaction_cst=[{0, 1}]`
+        is equivalent to `interaction_cst=[{0, 1}, {2, 3, 4}]`,
+        and specifies that each branch of a tree will either only split
+        on features 0 and 1 or only split on features 2, 3 and 4.
+
+        .. versionadded:: 1.2
+
     warm_start : bool, 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
@@ -1653,6 +1731,7 @@ def __init__(
         max_bins=255,
         categorical_features=None,
         monotonic_cst=None,
+        interaction_cst=None,
         warm_start=False,
         early_stopping="auto",
         scoring="loss",
@@ -1674,6 +1753,7 @@ def __init__(
             max_bins=max_bins,
             categorical_features=categorical_features,
             monotonic_cst=monotonic_cst,
+            interaction_cst=interaction_cst,
             warm_start=warm_start,
             early_stopping=early_stopping,
             scoring=scoring,
diff --git a/sklearn/ensemble/_hist_gradient_boosting/grower.py b/sklearn/ensemble/_hist_gradient_boosting/grower.py
index 4dc4a0ee3a7a4..5e3010fa4a509 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/grower.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/grower.py
@@ -55,6 +55,8 @@ class TreeNode:
         The sum of the hessians of the samples at the node.
     split_info : SplitInfo or None
         The result of the split evaluation.
+    is_leaf : bool
+        True if node is a leaf
     left_child : TreeNode or None
         The left child of the node. None for leaves.
     right_child : TreeNode or None
@@ -66,6 +68,14 @@ class TreeNode:
         start position of the node's sample_indices in splitter.partition.
     partition_stop : int
         stop position of the node's sample_indices in splitter.partition.
+    allowed_features : None or ndarray, dtype=int
+        Indices of features allowed to split for children.
+    interaction_cst_indices : None or list of ints
+        Indices of the interaction sets that have to be applied on splits of
+        child nodes. The fewer sets the stronger the constraint as fewer sets
+        contain fewer features.
+    children_lower_bound : float
+    children_upper_bound : float
     """
 
     split_info = None
@@ -92,6 +102,8 @@ def __init__(self, depth, sample_indices, sum_gradients, sum_hessians, value=Non
         self.sum_hessians = sum_hessians
         self.value = value
         self.is_leaf = False
+        self.allowed_features = None
+        self.interaction_cst_indices = None
         self.set_children_bounds(float("-inf"), float("+inf"))
 
     def set_children_bounds(self, lower, upper):
@@ -170,6 +182,8 @@ class TreeGrower:
           - -1: monotonic decrease
 
         Read more in the :ref:`User Guide <monotonic_cst_gbdt>`.
+    interaction_cst : list of sets of integers, default=None
+        List of interaction constraints.
     l2_regularization : float, default=0.
         The L2 regularization parameter.
     min_hessian_to_split : float, default=1e-3
@@ -222,6 +236,7 @@ def __init__(
         has_missing_values=False,
         is_categorical=None,
         monotonic_cst=None,
+        interaction_cst=None,
         l2_regularization=0.0,
         min_hessian_to_split=1e-3,
         shrinkage=1.0,
@@ -308,6 +323,7 @@ def __init__(
         self.max_leaf_nodes = max_leaf_nodes
         self.has_missing_values = has_missing_values
         self.monotonic_cst = monotonic_cst
+        self.interaction_cst = interaction_cst
         self.is_categorical = is_categorical
         self.l2_regularization = l2_regularization
         self.n_features = X_binned.shape[1]
@@ -402,6 +418,14 @@ def _intilialize_root(self, gradients, hessians, hessians_are_constant):
         self.root.histograms = self.histogram_builder.compute_histograms_brute(
             self.root.sample_indices
         )
+
+        if self.interaction_cst is not None:
+            self.root.interaction_cst_indices = range(len(self.interaction_cst))
+            allowed_features = set().union(*self.interaction_cst)
+            self.root.allowed_features = np.fromiter(
+                allowed_features, dtype=np.uint32, count=len(allowed_features)
+            )
+
         self._compute_best_split_and_push(self.root)
 
     def _compute_best_split_and_push(self, node):
@@ -414,13 +438,14 @@ def _compute_best_split_and_push(self, node):
         """
 
         node.split_info = self.splitter.find_node_split(
-            node.n_samples,
-            node.histograms,
-            node.sum_gradients,
-            node.sum_hessians,
-            node.value,
-            node.children_lower_bound,
-            node.children_upper_bound,
+            n_samples=node.n_samples,
+            histograms=node.histograms,
+            sum_gradients=node.sum_gradients,
+            sum_hessians=node.sum_hessians,
+            value=node.value,
+            lower_bound=node.children_lower_bound,
+            upper_bound=node.children_upper_bound,
+            allowed_features=node.allowed_features,
         )
 
         if node.split_info.gain <= 0:  # no valid split
@@ -477,6 +502,19 @@ def split_next(self):
         right_child_node.partition_start = left_child_node.partition_stop
         right_child_node.partition_stop = node.partition_stop
 
+        # set interaction constraints (the indices of the constraints sets)
+        if self.interaction_cst is not None:
+            # Calculate allowed_features and interaction_cst_indices only once. Child
+            # nodes inherit them before they get split.
+            (
+                left_child_node.allowed_features,
+                left_child_node.interaction_cst_indices,
+            ) = self._compute_interactions(node)
+            right_child_node.interaction_cst_indices = (
+                left_child_node.interaction_cst_indices
+            )
+            right_child_node.allowed_features = left_child_node.allowed_features
+
         if not self.has_missing_values[node.split_info.feature_idx]:
             # If no missing values are encountered at fit time, then samples
             # with missing values during predict() will go to whichever child
@@ -578,6 +616,48 @@ def split_next(self):
 
         return left_child_node, right_child_node
 
+    def _compute_interactions(self, node):
+        r"""Compute features allowed by interactions to be inherited by child nodes.
+
+        Example: Assume constraints [{0, 1}, {1, 2}].
+           1      <- Both constraint groups could be applied from now on
+          / \
+         1   2    <- Left split still fulfills both constraint groups.
+        / \ / \      Right split at feature 2 has only group {1, 2} from now on.
+
+        LightGBM uses the same logic for overlapping groups. See
+        https://github.com/microsoft/LightGBM/issues/4481 for details.
+
+        Parameters:
+        ----------
+        node : TreeNode
+            A node that might have children. Based on its feature_idx, the interaction
+            constraints for possible child nodes are computed.
+
+        Returns
+        -------
+        allowed_features : ndarray, dtype=uint32
+            Indices of features allowed to split for children.
+        interaction_cst_indices : list of ints
+            Indices of the interaction sets that have to be applied on splits of
+            child nodes. The fewer sets the stronger the constraint as fewer sets
+            contain fewer features.
+        """
+        # Note:
+        #  - Case of no interactions is already captured before function call.
+        #  - This is for nodes that are already split and have a
+        #    node.split_info.feature_idx.
+        allowed_features = set()
+        interaction_cst_indices = []
+        for i in node.interaction_cst_indices:
+            if node.split_info.feature_idx in self.interaction_cst[i]:
+                interaction_cst_indices.append(i)
+                allowed_features.update(self.interaction_cst[i])
+        return (
+            np.fromiter(allowed_features, dtype=np.uint32, count=len(allowed_features)),
+            interaction_cst_indices,
+        )
+
     def _finalize_leaf(self, node):
         """Make node a leaf of the tree being grown."""
 
diff --git a/sklearn/ensemble/_hist_gradient_boosting/splitting.pyx b/sklearn/ensemble/_hist_gradient_boosting/splitting.pyx
index 89f4e5baebc3e..f6630efd28a0f 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/splitting.pyx
+++ b/sklearn/ensemble/_hist_gradient_boosting/splitting.pyx
@@ -421,6 +421,7 @@ cdef class Splitter:
             const Y_DTYPE_C value,
             const Y_DTYPE_C lower_bound=-INFINITY,
             const Y_DTYPE_C upper_bound=INFINITY,
+            const unsigned int [:] allowed_features=None,
             ):
         """For each feature, find the best bin to split on at a given node.
 
@@ -453,6 +454,9 @@ cdef class Splitter:
         upper_bound : float
             Upper bound for the children values for respecting the monotonic
             constraints.
+        allowed_features : None or ndarray, dtype=np.uint32
+            Indices of the features that are allowed by interaction constraints to be
+            split.
 
         Returns
         -------
@@ -461,38 +465,53 @@ cdef class Splitter:
         """
         cdef:
             int feature_idx
-            int best_feature_idx
-            int n_features = self.n_features
+            int split_info_idx
+            int best_split_info_idx
+            int n_allowed_features
             split_info_struct split_info
             split_info_struct * split_infos
             const unsigned char [::1] has_missing_values = self.has_missing_values
             const unsigned char [::1] is_categorical = self.is_categorical
             const signed char [::1] monotonic_cst = self.monotonic_cst
             int n_threads = self.n_threads
+            bint has_interaction_cst = False
+
+        has_interaction_cst = allowed_features is not None
+        if has_interaction_cst:
+            n_allowed_features = allowed_features.shape[0]
+        else:
+            n_allowed_features = self.n_features
 
         with nogil:
 
             split_infos = <split_info_struct *> malloc(
-                self.n_features * sizeof(split_info_struct))
+                n_allowed_features * sizeof(split_info_struct))
+
+            # split_info_idx is index of split_infos of size n_features_allowed
+            # features_idx is the index of the feature column in X
+            for split_info_idx in prange(n_allowed_features, schedule='static',
+                                         num_threads=n_threads):
+                if has_interaction_cst:
+                    feature_idx = allowed_features[split_info_idx]
+                else:
+                    feature_idx = split_info_idx
 
-            for feature_idx in prange(n_features, schedule='static',
-                                      num_threads=n_threads):
-                split_infos[feature_idx].feature_idx = feature_idx
+                split_infos[split_info_idx].feature_idx = feature_idx
 
                 # For each feature, find best bin to split on
                 # Start with a gain of -1 (if no better split is found, that
                 # means one of the constraints isn't respected
                 # (min_samples_leaf, etc) and the grower will later turn the
                 # node into a leaf.
-                split_infos[feature_idx].gain = -1
-                split_infos[feature_idx].is_categorical = is_categorical[feature_idx]
+                split_infos[split_info_idx].gain = -1
+                split_infos[split_info_idx].is_categorical = is_categorical[feature_idx]
 
                 if is_categorical[feature_idx]:
                     self._find_best_bin_to_split_category(
                         feature_idx, has_missing_values[feature_idx],
                         histograms, n_samples, sum_gradients, sum_hessians,
                         value, monotonic_cst[feature_idx], lower_bound,
-                        upper_bound, &split_infos[feature_idx])
+                        upper_bound, &split_infos[split_info_idx])
                 else:
                     # We will scan bins from left to right (in all cases), and
                     # if there are any missing values, we will also scan bins
@@ -508,7 +527,7 @@ cdef class Splitter:
                         feature_idx, has_missing_values[feature_idx],
                         histograms, n_samples, sum_gradients, sum_hessians,
                         value, monotonic_cst[feature_idx],
-                        lower_bound, upper_bound, &split_infos[feature_idx])
+                        lower_bound, upper_bound, &split_infos[split_info_idx])
 
                     if has_missing_values[feature_idx]:
                         # We need to explore both directions to check whether
@@ -518,12 +537,14 @@ cdef class Splitter:
                             feature_idx, histograms, n_samples,
                             sum_gradients, sum_hessians,
                             value, monotonic_cst[feature_idx],
-                            lower_bound, upper_bound, &split_infos[feature_idx])
+                            lower_bound, upper_bound, &split_infos[split_info_idx])
 
             # then compute best possible split among all features
-            best_feature_idx = self._find_best_feature_to_split_helper(
-                split_infos)
-            split_info = split_infos[best_feature_idx]
+            # split_info is set to the best of split_infos
+            best_split_info_idx = self._find_best_feature_to_split_helper(
+                split_infos, n_allowed_features
+            )
+            split_info = split_infos[best_split_info_idx]
 
         out = SplitInfo(
             split_info.gain,
@@ -549,18 +570,19 @@ cdef class Splitter:
         return out
 
     cdef unsigned int _find_best_feature_to_split_helper(
-            self,
-            split_info_struct * split_infos) nogil:  # IN
-        """Returns the best feature among those in splits_infos."""
+        self,
+        split_info_struct * split_infos,  # IN
+        int n_allowed_features,
+    ) nogil:
+        """Return the index of split_infos with the best feature split."""
         cdef:
-            unsigned int feature_idx
-            unsigned int best_feature_idx = 0
-
-        for feature_idx in range(1, self.n_features):
-            if (split_infos[feature_idx].gain >
-                    split_infos[best_feature_idx].gain):
-                best_feature_idx = feature_idx
-        return best_feature_idx
+            unsigned int split_info_idx
+            unsigned int best_split_info_idx = 0
+
+        for split_info_idx in range(1, n_allowed_features):
+            if (split_infos[split_info_idx].gain > split_infos[best_split_info_idx].gain):
+                best_split_info_idx = split_info_idx
+        return best_split_info_idx
 
     cdef void _find_best_bin_to_split_left_to_right(
             Splitter self,

Test patch: diff --git a/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py b/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py
index 852b7386c7709..dcdd01c4f28ec 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/tests/test_gradient_boosting.py
@@ -47,6 +47,47 @@ def _make_dumb_dataset(n_samples):
     return X_dumb, y_dumb
 
 
+@pytest.mark.parametrize(
+    "GradientBoosting, X, y",
+    [
+        (HistGradientBoostingClassifier, X_classification, y_classification),
+        (HistGradientBoostingRegressor, X_regression, y_regression),
+    ],
+)
+@pytest.mark.parametrize(
+    "params, err_msg",
+    [
+        (
+            {"interaction_cst": "string"},
+            "",
+        ),
+        (
+            {"interaction_cst": [0, 1]},
+            "Interaction constraints must be None or an iterable of iterables",
+        ),
+        (
+            {"interaction_cst": [{0, 9999}]},
+            r"Interaction constraints must consist of integer indices in \[0,"
+            r" n_features - 1\] = \[.*\], specifying the position of features,",
+        ),
+        (
+            {"interaction_cst": [{-1, 0}]},
+            r"Interaction constraints must consist of integer indices in \[0,"
+            r" n_features - 1\] = \[.*\], specifying the position of features,",
+        ),
+        (
+            {"interaction_cst": [{0.5}]},
+            r"Interaction constraints must consist of integer indices in \[0,"
+            r" n_features - 1\] = \[.*\], specifying the position of features,",
+        ),
+    ],
+)
+def test_init_parameters_validation(GradientBoosting, X, y, params, err_msg):
+
+    with pytest.raises(ValueError, match=err_msg):
+        GradientBoosting(**params).fit(X, y)
+
+
 # TODO(1.3): remove
 @pytest.mark.filterwarnings("ignore::FutureWarning")
 def test_invalid_classification_loss():
@@ -1082,6 +1123,72 @@ def test_uint8_predict(Est):
     est.predict(X)
 
 
+@pytest.mark.parametrize(
+    "interaction_cst, n_features, result",
+    [
+        (None, 931, None),
+        ([{0, 1}], 2, [{0, 1}]),
+        ([(1, 0), [5, 1]], 6, [{0, 1}, {1, 5}, {2, 3, 4}]),
+    ],
+)
+def test_check_interaction_cst(interaction_cst, n_features, result):
+    """Check that _check_interaction_cst returns the expected list of sets"""
+    est = HistGradientBoostingRegressor()
+    est.set_params(interaction_cst=interaction_cst)
+    assert est._check_interaction_cst(n_features) == result
+
+
+def test_interaction_cst_numerically():
+    """Check that interaction constraints have no forbidden interactions."""
+    rng = np.random.RandomState(42)
+    n_samples = 1000
+    X = rng.uniform(size=(n_samples, 2))
+    # Construct y with a strong interaction term
+    # y = x0 + x1 + 5 * x0 * x1
+    y = np.hstack((X, 5 * X[:, [0]] * X[:, [1]])).sum(axis=1)
+
+    est = HistGradientBoostingRegressor(random_state=42)
+    est.fit(X, y)
+    est_no_interactions = HistGradientBoostingRegressor(
+        interaction_cst=[{0}, {1}], random_state=42
+    )
+    est_no_interactions.fit(X, y)
+
+    delta = 0.25
+    # Make sure we do not extrapolate out of the training set as tree-based estimators
+    # are very bad in doing so.
+    X_test = X[(X[:, 0] < 1 - delta) & (X[:, 1] < 1 - delta)]
+    X_delta_d_0 = X_test + [delta, 0]
+    X_delta_0_d = X_test + [0, delta]
+    X_delta_d_d = X_test + [delta, delta]
+
+    # Note: For the y from above as a function of x0 and x1, we have
+    # y(x0+d, x1+d) = y(x0, x1) + 5 * d * (2/5 + x0 + x1) + 5 * d**2
+    # y(x0+d, x1)   = y(x0, x1) + 5 * d * (1/5 + x1)
+    # y(x0,   x1+d) = y(x0, x1) + 5 * d * (1/5 + x0)
+    # Without interaction constraints, we would expect a result of 5 * d**2 for the
+    # following expression, but zero with constraints in place.
+    assert_allclose(
+        est_no_interactions.predict(X_delta_d_d)
+        + est_no_interactions.predict(X_test)
+        - est_no_interactions.predict(X_delta_d_0)
+        - est_no_interactions.predict(X_delta_0_d),
+        0,
+        atol=1e-12,
+    )
+
+    # Correct result of the expressions is 5 * delta**2. But this is hard to achieve by
+    # a fitted tree-based model. However, with 100 iterations the expression should
+    # at least be positive!
+    assert np.all(
+        est.predict(X_delta_d_d)
+        + est.predict(X_test)
+        - est.predict(X_delta_d_0)
+        - est.predict(X_delta_0_d)
+        > 0.01
+    )
+
+
 # TODO(1.3): Remove
 @pytest.mark.parametrize(
     "old_loss, new_loss, Estimator",
diff --git a/sklearn/ensemble/_hist_gradient_boosting/tests/test_grower.py b/sklearn/ensemble/_hist_gradient_boosting/tests/test_grower.py
index 5ae497d6e5845..c4ae90b7e7d96 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/tests/test_grower.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/tests/test_grower.py
@@ -567,3 +567,85 @@ def test_ohe_equivalence(min_samples_leaf, n_unique_categories, target):
         assert predictor.get_max_depth() < predictor_ohe.get_max_depth()
 
     np.testing.assert_allclose(preds, preds_ohe)
+
+
+def test_grower_interaction_constraints():
+    """Check that grower respects interaction constraints."""
+    n_features = 6
+    interaction_cst = [{0, 1}, {1, 2}, {3, 4, 5}]
+    n_samples = 10
+    n_bins = 6
+    root_feature_splits = []
+
+    def get_all_children(node):
+        res = []
+        if node.is_leaf:
+            return res
+        for n in [node.left_child, node.right_child]:
+            res.append(n)
+            res.extend(get_all_children(n))
+        return res
+
+    for seed in range(20):
+        rng = np.random.RandomState(seed)
+
+        X_binned = rng.randint(
+            0, n_bins - 1, size=(n_samples, n_features), dtype=X_BINNED_DTYPE
+        )
+        X_binned = np.asfortranarray(X_binned)
+        gradients = rng.normal(size=n_samples).astype(G_H_DTYPE)
+        hessians = np.ones(shape=1, dtype=G_H_DTYPE)
+
+        grower = TreeGrower(
+            X_binned,
+            gradients,
+            hessians,
+            n_bins=n_bins,
+            min_samples_leaf=1,
+            interaction_cst=interaction_cst,
+            n_threads=n_threads,
+        )
+        grower.grow()
+
+        root_feature_idx = grower.root.split_info.feature_idx
+        root_feature_splits.append(root_feature_idx)
+
+        feature_idx_to_constraint_set = {
+            0: {0, 1},
+            1: {0, 1, 2},
+            2: {1, 2},
+            3: {3, 4, 5},
+            4: {3, 4, 5},
+            5: {3, 4, 5},
+        }
+
+        root_constraint_set = feature_idx_to_constraint_set[root_feature_idx]
+        for node in (grower.root.left_child, grower.root.right_child):
+            # Root's children's allowed_features must be the root's constraints set.
+            assert_array_equal(node.allowed_features, list(root_constraint_set))
+        for node in get_all_children(grower.root):
+            if node.is_leaf:
+                continue
+            # Ensure that each node uses a subset of features of its parent node.
+            parent_interaction_cst_indices = set(node.interaction_cst_indices)
+            right_interactions_cst_indices = set(
+                node.right_child.interaction_cst_indices
+            )
+            left_interactions_cst_indices = set(node.left_child.interaction_cst_indices)
+
+            assert right_interactions_cst_indices.issubset(
+                parent_interaction_cst_indices
+            )
+            assert left_interactions_cst_indices.issubset(
+                parent_interaction_cst_indices
+            )
+            # The features used for split must have been present in the root's
+            # constraint set.
+            assert node.split_info.feature_idx in root_constraint_set
+
+    # Make sure that every feature is used at least once as split for the root node.
+    assert (
+        len(set(root_feature_splits))
+        == len(set().union(*interaction_cst))
+        == n_features
+    )
diff --git a/sklearn/ensemble/_hist_gradient_boosting/tests/test_splitting.py b/sklearn/ensemble/_hist_gradient_boosting/tests/test_splitting.py
index 0d19bdc6df72b..d1da34015a2a4 100644
--- a/sklearn/ensemble/_hist_gradient_boosting/tests/test_splitting.py
+++ b/sklearn/ensemble/_hist_gradient_boosting/tests/test_splitting.py
@@ -856,3 +856,100 @@ def test_splitting_categorical_sanity(
     left_mask = np.isin(X_binned.ravel(), expected_categories_left)
     assert_array_equal(sample_indices[left_mask], samples_left)
     assert_array_equal(sample_indices[~left_mask], samples_right)
+
+
+def test_split_interaction_constraints():
+    """Check that allowed_features are respected."""
+    n_features = 4
+    # features 1 and 2 are not allowed to be split on
+    allowed_features = np.array([0, 3], dtype=np.uint32)
+    n_bins = 5
+    n_samples = 10
+    l2_regularization = 0.0
+    min_hessian_to_split = 1e-3
+    min_samples_leaf = 1
+    min_gain_to_split = 0.0
+
+    sample_indices = np.arange(n_samples, dtype=np.uint32)
+    all_hessians = np.ones(1, dtype=G_H_DTYPE)
+    sum_hessians = n_samples
+    hessians_are_constant = True
+
+    split_features = []
+
+    # The loop is to ensure that we split at least once on each allowed feature (0, 3).
+    # This is tracked by split_features and checked at the end.
+    for i in range(10):
+        rng = np.random.RandomState(919 + i)
+        X_binned = np.asfortranarray(
+            rng.randint(0, n_bins - 1, size=(n_samples, n_features)),
+            dtype=X_BINNED_DTYPE,
+        )
+        X_binned = np.asfortranarray(X_binned, dtype=X_BINNED_DTYPE)
+
+        # Make feature 1 very important
+        all_gradients = (10 * X_binned[:, 1] + rng.randn(n_samples)).astype(G_H_DTYPE)
+        sum_gradients = all_gradients.sum()
+
+        builder = HistogramBuilder(
+            X_binned,
+            n_bins,
+            all_gradients,
+            all_hessians,
+            hessians_are_constant,
+            n_threads,
+        )
+        n_bins_non_missing = np.array([n_bins] * X_binned.shape[1], dtype=np.uint32)
+        has_missing_values = np.array([False] * X_binned.shape[1], dtype=np.uint8)
+        monotonic_cst = np.array(
+            [MonotonicConstraint.NO_CST] * X_binned.shape[1], dtype=np.int8
+        )
+        is_categorical = np.zeros_like(monotonic_cst, dtype=np.uint8)
+        missing_values_bin_idx = n_bins - 1
+        splitter = Splitter(
+            X_binned,
+            n_bins_non_missing,
+            missing_values_bin_idx,
+            has_missing_values,
+            is_categorical,
+            monotonic_cst,
+            l2_regularization,
+            min_hessian_to_split,
+            min_samples_leaf,
+            min_gain_to_split,
+            hessians_are_constant,
+        )
+
+        assert np.all(sample_indices == splitter.partition)
+
+        histograms = builder.compute_histograms_brute(sample_indices)
+        value = compute_node_value(
+            sum_gradients, sum_hessians, -np.inf, np.inf, l2_regularization
+        )
+
+        # with all features allowed, feature 1 should be split on as it is the most
+        # important one by construction of the gradients
+        si_root = splitter.find_node_split(
+            n_samples,
+            histograms,
+            sum_gradients,
+            sum_hessians,
+            value,
+            allowed_features=None,
+        )
+        assert si_root.feature_idx == 1
+
+        # only features 0 and 3 are allowed to be split on
+        si_root = splitter.find_node_split(
+            n_samples,
+            histograms,
+            sum_gradients,
+            sum_hessians,
+            value,
+            allowed_features=allowed_features,
+        )
+        split_features.append(si_root.feature_idx)
+        assert si_root.feature_idx in allowed_features
+
+    # make sure feature 0 and feature 3 are split on in the constraint setting
+    assert set(allowed_features) == set(split_features)