'hist': Montonic Constraints (#3085)

* Extended monotonic constraints support to 'hist' tree method.

* Added monotonic constraints tests.

* Fix the signature of NoConstraint::CalcSplitGain()

* Document monotonic constraint support in 'hist'

* Update signature of Update to account for latest refactor

doc/tutorials/monotonic.md

@@ -76,3 +76,15 @@ Some other examples:
 
 - ```(1,0)```: An increasing constraint on the first predictor and no constraint on the second.
 - ```(0,-1)```: No constraint on the first predictor and a decreasing constraint on the second.
+
+**Choise of tree construction algorithm**. To use monotonic constraints, be
+sure to set the `tree_method` parameter to one of `'exact'`, `'hist'`, and
+`'gpu_hist'`.
+
+**Note for the `'hist'` tree construction algorithm**.
+If `tree_method` is set to either `'hist'` or `'gpu_hist'`, enabling monotonic
+constraints may produce unnecessarily shallow trees. This is because the
+`'hist'` method reduces the number of candidate splits to be considered at each
+split. Monotonic constraints may wipe out all available split candidates, in
+which case no split is made. To reduce the effect, you may want to increase
+the `max_bin` parameter to consider more split candidates.

src/tree/param.h

@@ -376,7 +376,7 @@ struct NoConstraint {
   inline static void Init(TrainParam *param, unsigned num_feature) {
     param->monotone_constraints.resize(num_feature, 0);
   }
-  inline double CalcSplitGain(const TrainParam &param, bst_uint split_index,
+  inline double CalcSplitGain(const TrainParam &param, int constraint,
                               GradStats left, GradStats right) const {
     return left.CalcGain(param) + right.CalcGain(param);
   }
@@ -421,6 +421,7 @@ template <typename param_t>
 template <typename param_t>
   XGBOOST_DEVICE inline double CalcSplitGain(const param_t &param, int constraint,
                               GradStats left, GradStats right) const {
+    const double negative_infinity = -std::numeric_limits<double>::infinity();
     double wleft = CalcWeight(param, left);
     double wright = CalcWeight(param, right);
     double gain =
@@ -429,9 +430,9 @@ template <typename param_t>
     if (constraint == 0) {
       return gain;
     } else if (constraint > 0) {
-      return wleft < wright ? gain : 0.0;
+      return wleft <= wright ? gain : negative_infinity;
     } else {
-      return wleft > wright ? gain : 0.0;
+      return wleft >= wright ? gain : negative_infinity;
     }
   }
 

src/tree/updater_fast_hist.cc

@@ -870,13 +870,13 @@ class FastHistMaker: public TreeUpdater {
             if (d_step > 0) {
               // forward enumeration: split at right bound of each bin
               loss_chg = static_cast<bst_float>(
-                  constraint.CalcSplitGain(param, fid, e, c) -
+                  constraint.CalcSplitGain(param, param.monotone_constraints[fid], e, c) -
                   snode.root_gain);
               split_pt = cut_val[i];
             } else {
               // backward enumeration: split at left bound of each bin
               loss_chg = static_cast<bst_float>(
-                  constraint.CalcSplitGain(param, fid, c, e) -
+                  constraint.CalcSplitGain(param, param.monotone_constraints[fid], c, e) -
                   snode.root_gain);
               if (i == imin) {
                 // for leftmost bin, left bound is the smallest feature value
@@ -961,10 +961,45 @@ class FastHistMaker: public TreeUpdater {
   std::unique_ptr<TreeUpdater> pruner_;
 };
 
+// simple switch to defer implementation.
+class FastHistTreeUpdaterSwitch : public TreeUpdater {
+ public:
+  FastHistTreeUpdaterSwitch() : monotone_(false) {}
+  void Init(const std::vector<std::pair<std::string, std::string> >& args) override {
+    for (auto &kv : args) {
+      if (kv.first == "monotone_constraints" && kv.second.length() != 0) {
+        monotone_ = true;
+      }
+    }
+    if (inner_.get() == nullptr) {
+      if (monotone_) {
+        inner_.reset(new FastHistMaker<GradStats, ValueConstraint>());
+      } else {
+        inner_.reset(new FastHistMaker<GradStats, NoConstraint>());
+      }
+    }
+
+    inner_->Init(args);
+  }
+
+  void Update(HostDeviceVector<bst_gpair>* gpair,
+              DMatrix* data,
+              const std::vector<RegTree*>& trees) override {
+    CHECK(inner_ != nullptr);
+    inner_->Update(gpair, data, trees);
+  }
+
+ private:
+  //  monotone constraints
+  bool monotone_;
+  // internal implementation
+  std::unique_ptr<TreeUpdater> inner_;
+};
+
 XGBOOST_REGISTER_TREE_UPDATER(FastHistMaker, "grow_fast_histmaker")
 .describe("Grow tree using quantized histogram.")
 .set_body([]() {
-    return new FastHistMaker<GradStats, NoConstraint>();
+    return new FastHistTreeUpdaterSwitch();
   });
 
 }  // namespace tree

tests/python/test_monotone_constraints.py

@@ -0,0 +1,85 @@
+import numpy as np
+import xgboost as xgb
+import unittest
+
+
+def is_increasing(y):
+    return np.count_nonzero(np.diff(y) < 0.0) == 0
+
+
+def is_decreasing(y):
+    return np.count_nonzero(np.diff(y) > 0.0) == 0
+
+
+def is_correctly_constrained(learner):
+    n = 100
+    variable_x = np.linspace(0, 1, n).reshape((n, 1))
+    fixed_xs_values = np.linspace(0, 1, n)
+
+    for i in range(n):
+        fixed_x = fixed_xs_values[i] * np.ones((n, 1))
+        monotonically_increasing_x = np.column_stack((variable_x, fixed_x))
+        monotonically_increasing_dset = xgb.DMatrix(monotonically_increasing_x)
+        monotonically_increasing_y = learner.predict(
+            monotonically_increasing_dset
+        )
+
+        monotonically_decreasing_x = np.column_stack((fixed_x, variable_x))
+        monotonically_decreasing_dset = xgb.DMatrix(monotonically_decreasing_x)
+        monotonically_decreasing_y = learner.predict(
+            monotonically_decreasing_dset
+        )
+
+        if not (
+            is_increasing(monotonically_increasing_y) and
+            is_decreasing(monotonically_decreasing_y)
+        ):
+            return False
+
+    return True
+
+
+number_of_dpoints = 1000
+x1_positively_correlated_with_y = np.random.random(size=number_of_dpoints)
+x2_negatively_correlated_with_y = np.random.random(size=number_of_dpoints)
+
+x = np.column_stack((
+    x1_positively_correlated_with_y, x2_negatively_correlated_with_y
+))
+zs = np.random.normal(loc=0.0, scale=0.01, size=number_of_dpoints)
+y = (
+    5 * x1_positively_correlated_with_y +
+    np.sin(10 * np.pi * x1_positively_correlated_with_y) -
+    5 * x2_negatively_correlated_with_y -
+    np.cos(10 * np.pi * x2_negatively_correlated_with_y) +
+    zs
+)
+training_dset = xgb.DMatrix(x, label=y)
+
+
+class TestMonotoneConstraints(unittest.TestCase):
+
+    def test_monotone_constraints_for_exact_tree_method(self):
+
+        # first check monotonicity for the 'exact' tree method
+        params_for_constrained_exact_method = {
+            'tree_method': 'exact', 'silent': 1,
+            'monotone_constraints': '(1, -1)'
+        }
+        constrained_exact_method = xgb.train(
+            params_for_constrained_exact_method, training_dset
+        )
+        assert is_correctly_constrained(constrained_exact_method)
+
+    def test_monotone_constraints_for_hist_tree_method(self):
+
+        # next check monotonicity for the 'hist' tree method
+        params_for_constrained_hist_method = {
+            'tree_method': 'hist', 'silent': 1,
+            'monotone_constraints': '(1, -1)'
+        }
+        constrained_hist_method = xgb.train(
+            params_for_constrained_hist_method, training_dset
+        )
+
+        assert is_correctly_constrained(constrained_hist_method)