[TREE] add interaction constraints (#3466)

* add interaction constraints

* enable both interaction and monotonic constraints at the same time

* fix lint

* add R test, fix lint, update demo

* Use dmlc::JSONReader to express interaction constraints as nested lists; Use sparse arrays for bookkeeping

* Add Python test for interaction constraints

* make R interaction constraints parameter based on feature index instead of column names, fix R coding style

* Fix lint

* Add BlueTea88 to CONTRIBUTORS.md

* Short circuit when no constraint is specified; address review comments

* Add tutorial for feature interaction constraints

* allow interaction constraints to be passed as string, remove redundant column_names argument

* Fix typo

* Address review comments

* Add comments to Python test

R-package/R/utils.R

@@ -74,6 +74,19 @@ check.booster.params <- function(params, ...) {
         params[['monotone_constraints']] = vec2str
   }
   
+  # interaction constraints parser (convert from list of column indices to string)
+  if (!is.null(params[['interaction_constraints']]) && 
+      typeof(params[['interaction_constraints']]) != "character"){
+    # check input class
+    if (class(params[['interaction_constraints']]) != 'list') stop('interaction_constraints should be class list')
+    if (!all(unique(sapply(params[['interaction_constraints']], class)) %in% c('numeric','integer'))) {
+      stop('interaction_constraints should be a list of numeric/integer vectors')
+    }
+
+    # recast parameter as string
+    interaction_constraints <- sapply(params[['interaction_constraints']], function(x) paste0('[', paste(x, collapse=','), ']'))
+    params[['interaction_constraints']] <- paste0('[', paste(interaction_constraints, collapse=','), ']')
+  }
   return(params)
 }
 

R-package/R/xgb.train.R

@@ -26,6 +26,7 @@
 #'   \item \code{colsample_bytree} subsample ratio of columns when constructing each tree. Default: 1
 #'   \item \code{num_parallel_tree} Experimental parameter. number of trees to grow per round. Useful to test Random Forest through Xgboost (set \code{colsample_bytree < 1}, \code{subsample  < 1}  and \code{round = 1}) accordingly. Default: 1
 #'   \item \code{monotone_constraints} A numerical vector consists of \code{1}, \code{0} and \code{-1} with its length equals to the number of features in the training data. \code{1} is increasing, \code{-1} is decreasing and \code{0} is no constraint.
+#'   \item \code{interaction_constraints} A list of vectors specifying feature indices of permitted interactions. Each item of the list represents one permitted interaction where specified features are allowed to interact with each other. Feature index values should start from \code{0} (\code{0} references the first column).  Leave argument unspecified for no interaction constraints.
 #' }
 #' 
 #' 2.2. Parameter for Linear Booster

R-package/demo/interaction_constraints.R

@@ -0,0 +1,105 @@
+library(xgboost)
+library(data.table)
+
+set.seed(1024)
+
+# Function to obtain a list of interactions fitted in trees, requires input of maximum depth
+treeInteractions <- function(input_tree, input_max_depth){
+  trees <- copy(input_tree)  # copy tree input to prevent overwriting
+  if (input_max_depth < 2) return(list())  # no interactions if max depth < 2
+  if (nrow(input_tree) == 1) return(list())
+
+  # Attach parent nodes
+  for (i in 2:input_max_depth){
+    if (i == 2) trees[, ID_merge:=ID] else trees[, ID_merge:=get(paste0('parent_',i-2))]
+    parents_left <- trees[!is.na(Split), list(i.id=ID, i.feature=Feature, ID_merge=Yes)]
+    parents_right <- trees[!is.na(Split), list(i.id=ID, i.feature=Feature, ID_merge=No)]
+
+    setorderv(trees, 'ID_merge')
+    setorderv(parents_left, 'ID_merge')
+    setorderv(parents_right, 'ID_merge')
+
+    trees <- merge(trees, parents_left, by='ID_merge', all.x=T)
+    trees[!is.na(i.id), c(paste0('parent_', i-1), paste0('parent_feat_', i-1)):=list(i.id, i.feature)]
+    trees[, c('i.id','i.feature'):=NULL]
+
+    trees <- merge(trees, parents_right, by='ID_merge', all.x=T)
+    trees[!is.na(i.id), c(paste0('parent_', i-1), paste0('parent_feat_', i-1)):=list(i.id, i.feature)]
+    trees[, c('i.id','i.feature'):=NULL]
+  }
+
+  # Extract nodes with interactions
+  interaction_trees <- trees[!is.na(Split) & !is.na(parent_1), 
+                             c('Feature',paste0('parent_feat_',1:(input_max_depth-1))), with=F]
+  interaction_trees_split <- split(interaction_trees, 1:nrow(interaction_trees))
+  interaction_list <- lapply(interaction_trees_split, as.character)
+
+  # Remove NAs (no parent interaction)
+  interaction_list <- lapply(interaction_list, function(x) x[!is.na(x)])
+
+  # Remove non-interactions (same variable)
+  interaction_list <- lapply(interaction_list, unique)  # remove same variables
+  interaction_length <- sapply(interaction_list, length)
+  interaction_list <- interaction_list[interaction_length > 1]
+  interaction_list <- unique(lapply(interaction_list, sort))
+  return(interaction_list)
+}
+
+# Generate sample data
+x <- list()
+for (i in 1:10){
+  x[[i]] = i*rnorm(1000, 10)
+}
+x <- as.data.table(x)
+
+y = -1*x[, rowSums(.SD)] + x[['V1']]*x[['V2']] + x[['V3']]*x[['V4']]*x[['V5']] + rnorm(1000, 0.001) + 3*sin(x[['V7']])
+
+train = as.matrix(x)
+
+# Interaction constraint list (column names form)
+interaction_list <- list(c('V1','V2'),c('V3','V4','V5'))
+
+# Convert interaction constraint list into feature index form
+cols2ids <- function(object, col_names) {
+  LUT <- seq_along(col_names) - 1
+  names(LUT) <- col_names
+  rapply(object, function(x) LUT[x], classes="character", how="replace")
+}
+interaction_list_fid = cols2ids(interaction_list, colnames(train))
+
+# Fit model with interaction constraints
+bst = xgboost(data = train, label = y, max_depth = 4,
+              eta = 0.1, nthread = 2, nrounds = 1000,
+              interaction_constraints = interaction_list_fid)
+
+bst_tree <- xgb.model.dt.tree(colnames(train), bst)
+bst_interactions <- treeInteractions(bst_tree, 4)  # interactions constrained to combinations of V1*V2 and V3*V4*V5
+
+# Fit model without interaction constraints
+bst2 = xgboost(data = train, label = y, max_depth = 4,
+               eta = 0.1, nthread = 2, nrounds = 1000)
+
+bst2_tree <- xgb.model.dt.tree(colnames(train), bst2)
+bst2_interactions <- treeInteractions(bst2_tree, 4)  # much more interactions
+
+# Fit model with both interaction and monotonicity constraints
+bst3 = xgboost(data = train, label = y, max_depth = 4,
+               eta = 0.1, nthread = 2, nrounds = 1000,
+               interaction_constraints = interaction_list_fid,
+               monotone_constraints = c(-1,0,0,0,0,0,0,0,0,0))
+
+bst3_tree <- xgb.model.dt.tree(colnames(train), bst3)
+bst3_interactions <- treeInteractions(bst3_tree, 4)  # interactions still constrained to combinations of V1*V2 and V3*V4*V5
+
+# Show monotonic constraints still apply by checking scores after incrementing V1
+x1 <- sort(unique(x[['V1']]))
+for (i in 1:length(x1)){
+  testdata <- copy(x[, -c('V1')])
+  testdata[['V1']] <- x1[i]
+  testdata <- testdata[, paste0('V',1:10), with=F]
+  pred <- predict(bst3, as.matrix(testdata))
+  
+  # Should not print out anything due to monotonic constraints
+  if (i > 1) if (any(pred > prev_pred)) print(i)
+  prev_pred <- pred 
+}

R-package/tests/testthat/test_interaction_constraints.R

@@ -0,0 +1,38 @@
+require(xgboost)
+
+context("interaction constraints")
+
+set.seed(1024)
+x1 <- rnorm(1000, 1)
+x2 <- rnorm(1000, 1)
+x3 <- sample(c(1,2,3), size=1000, replace=TRUE)
+y <- x1 + x2 + x3 + x1*x2*x3 + rnorm(1000, 0.001) + 3*sin(x1)
+train <- matrix(c(x1,x2,x3), ncol = 3)
+
+test_that("interaction constraints for regression", {
+  # Fit a model that only allows interaction between x1 and x2
+  bst <- xgboost(data = train, label = y, max_depth = 3,
+                 eta = 0.1, nthread = 2, nrounds = 100, verbose = 0,
+                 interaction_constraints = list(c(0,1)))
+  
+  # Set all observations to have the same x3 values then increment
+  #  by the same amount
+	preds <- lapply(c(1,2,3), function(x){
+		tmat <- matrix(c(x1,x2,rep(x,1000)), ncol=3)
+		return(predict(bst, tmat))
+	})
+
+  # Check incrementing x3 has the same effect on all observations
+  #   since x3 is constrained to be independent of x1 and x2
+  #   and all observations start off from the same x3 value
+	diff1 <- preds[[2]] - preds[[1]]
+	test1 <- all(abs(diff1 - diff1[1]) < 1e-4)
+	
+	diff2 <- preds[[3]] - preds[[2]]
+	test2 <- all(abs(diff2 - diff2[1]) < 1e-4)
+	
+  expect_true({
+    test1 & test2
+  }, "Interaction Contraint Satisfied")
+  
+})