Output-Constrained Decision Trees

## Abstract

Abstract When there is a correlation between any pair of targets, one needs a prediction method that can handle vector-valued output. In this setting, multi-target learning is particularly important as it is widely used in various applications. This paper introduces new variants of decision trees that can handle not only multi-target output but also the constraints among the targets. We focus on the customization of conventional decision trees by adjusting the splitting criteria to handle the constraints and obtain feasible predictions. We present both an optimization-based exact approach and several heuristics, complete with a discussion on their respective advantages and disadvantages. To support our findings, we conduct a computational study to demonstrate and compare the results of the proposed approaches.

## 1 Introduction

Decision trees (DTs) play a fundamental role in machine learning and data science for several key reasons. First, DTs have proven to be accurate predictive models in a wide range of applications. Moreover, the heuristic approaches used for their construction are computationally efficient, resulting in extremely fast training times. Accuracy and speed make them ideal for large datasets or real-time applications. Second, DTs, particularly those that are not too deep, are inherently explainable. By following the branches of such a trained tree, one can easily understand the logic behind the model’s predictions. This feature is particularly valuable in sectors where interpretability is critical, such as in healthcare or finance. Lastly, DTs serve as the base learners in powerful ensemble methods. These methods, which include techniques like random forests and gradient boosting, utilize numerous DTs to create robust and accurate predictive models.

When there is a correlation between any pair of targets, we need a prediction method that can handle vector-valued output. In this setting, multi-target learning is particularly important as it is widely used in fields such as forecasting, monitoring, manufacturing, and recommendation systems. DTs can easily be customized for handling multi-target regression by modifying the splitting criteria to handle vector predictions at each node.

#### Motivation.

Despite the wealth of existing literature on DTs, to the best of our knowledge, none consider cases where constraints exist among multiple targets. This is notable, as in applications, decision-makers are aware of these constraints, and data analysts should impose them for feasible predictions. Here is an illustrative example from the field of operations management where machine learning approaches are becoming more prominent in the context of decision-making; see e.g., [].

Consider a dataset of products in a retail store, each with distinct attributes. This dataset includes the inventory amounts of these products across five different warehouses, represented as five-dimensional target vectors (training set). The goal is to predict the inventory amounts for new products across these warehouses (test set).

To achieve this, one can train a DT and obtain the partitions of the existing products represented by the leaf nodes of the trained tree. The prediction for a new product is then simply the average of the inventory amounts (mean vector) within the leaf node that it falls into. This leads to the prediction of five different inventory amounts, i.e., one vector for each warehouse.

However, it is crucial to consider an additional constraint in our scenario: Due to budget restrictions, a store can only transport a product in a maximum of two out of the five warehouses. This is also evident in the training data set where each target vector has at most two nonzero values. Consequently, if one continues to use the mean vector approach for prediction, it may not only yield inaccurate results but also infeasible predictions because it does not account for this constraint.

This erroneous estimation would also persist with many well-known DT training algorithms, since their estimation at a leaf node is also based on taking the mean or the median of the target vectors within the node. However, one may argue that growing a deeper DT may resolve the infeasibility issue as there will be fewer products at each leaf node. Even if we overlook the loss of interpretability due to the large size of such a tree, it is well-known that deep trees can easily overfit to the training data, potentially leading to poor generalization for new products. For instance, Figure 1 shows our results with one of our datasets from our computational study in Section 3. The bar plot demonstrates that feasible predictions can be achieved with larger tree depths, but this comes with a significant decrease in accuracy as shown in the scatter plot. The scatter plot, on the other hand, shows that while shallow trees offer much better accuracy at the expense of most of their predictions are being infeasible.

Figure: Figure 1: The trade-off between accuracy and feasibility. The red circle in the right plot marks the depth of the tree which returns the lowest mean squared error value for the test set. The corresponding red bar in the left plot shows the corresponding percentage of infeasible predictions. The sample size is 4000.
Refer to caption: extracted/5895402/motivating_ex.png

#### Output-Constraints.

To obtain accurate and feasible predictions, we propose an optimization-based exact approach. This approach incorporates target constraints into the splitting procedure of common binary tree construction methods. Our approach ensures that the estimated target vector at each node meets the required constraints. However, we observe that the exact approach could have a lengthy computation time, especially for deep trees. To address this, we offer two heuristic approaches. These are based on using existing samples as predictions or relaxing the constraints with a penalty function at each split operation.

Up until now, our focus has been on regression problems. It is reasonable to question why we have not considered multi-target classification problems as well. The reason is, if the number of classes for each target and the number of targets are relatively small, we can assign a label to each feasible combination of target values. This approach effectively reduces the problem to a standard single-target classification problem, although possibly with many labels.

#### Related Work.

As machine learning methods improve in accuracy, they may still produce predictions and outcomes that lack physical consistency or feasibility, due to biases or data sparsity []. To overcome this challenge, hybrid methods that integrate optimization techniques with machine learning models have gained significant attention in recent years due to their ability to leverage the strengths of both approaches. One key aspect of hybrid approaches is their ability to handle constraints from various sources to enhance the feasibility and accuracy of machine learning methods. These constraints can originate from the dataset, user-defined specifications, or domain knowledge []. Such machine learning approaches with imposed constraints are developed as practical applications across different domains, such as mechanics [], geomechanics [], and electronics [].

Traditional supervised learning methods primarily focus on predicting a single output, using data comprised of instances represented by features and a single response variable. Multi-target regression (also referred to as multi-output regression), is used to predict multiple continuous target variables simultaneously. One of the earliest methods for multi-target regression with decision trees was introduced by []. He extended the univariate recursive partitioning approach from Classification and Regression Trees (CART) [] to handle multiple responses, calling this method Multi-Response Trees (MRT). MRTs are built similarly to CART, but the node impurity is redefined as the sum of squared errors across all response variables. [] developed a system for building multi-objective regression trees (MORTs), which include constraints on tree size and accuracy, allowing users to balance interpretability and accuracy by limiting tree size or setting a minimum accuracy threshold. [] introduced a novel method combining multi-target regression with unsupervised learning. Specifically, they modified tree construction by incorporating a similarity-based clustering indicator in ensemble methods.

Decision trees (DTs) offer some flexibility in enforcing constraints, and many hybrid DT approaches incorporate *instance-level* constraints to bind a subset of instances in generating meaningfully combined results. For instance, the clustering tree algorithm by [] integrates domain knowledge to impose must-link and cannot-link constraints. DT methods have also been developed to address adversarial examples, which can be considered instance-level constraints [, ], or attribute-level constraints to assign several features to the predictors []. [] was one of the first to apply *attribute-level* constraints by imposing monotonicity among the model’s predictors, and [] extended this in an ordinal classification setting.

One of the earliest hybrid approaches of imposed constraints was implemented by [], who introduced constraints to deep learning architectures to approximate ordinary and partial differential equations. This approach laid the foundation for the general framework known as physics-informed neural networks, proposed by [], where constraints and boundary conditions are incorporated []. Similarly, [] integrated constraints and physical principles to generative diffusion processes. To the best of our knowledge, our proposed approach is the first method designed to handle *target-level* constraints in multi-target regression. This represents a novel contribution, building on the principles of constraint handling and extending them to more complex, interdependent target structures in predictive modeling.

#### Contributions.

We introduce novel approaches to handle constraints among target variables in DTs. The first method lends itself to an optimization-based exact method that ensures the feasibility of predictions by incorporating target constraints into the splitting procedure of common binary tree construction methods. We also propose two heuristic methods. These methods offer practical solutions to the computation time issue that might be associated with the exact approach, especially for deep trees. We provide the technical details of our proposed methods, exploring their strengths, potential limitations, and the conditions under which they would perform sensibly. To validate our approaches, we provide a reproducible computational study using our publicly available implementation.

## 2 Output-Constrained Decision Trees

Suppose that we have the training dataset $\mathcal{D}=\{(\mbox{\boldmath{$x$}}_{i},\mbox{\boldmath{$y$}}_{i}):i\in
\mathcal{I}_{\mathcal{D}}\}$ with $\mbox{\boldmath{$x$}}_{i}\in\mathbb{R}^{p}$ and $\mbox{\boldmath{$y$}}_{i}\in\mathbb{Y}\subseteq\mathbb{R}^{k}$ denoting the input vector and the target (output) vector for the data point $i$, respectively. The set $\mathbb{Y}$ shows the feasible set for the target vectors. It is important to observe that the target vector of each sample lies in $\mathbb{Y}$, and hence, they are also feasible.

A DT is grown by splitting the dataset into non-overlapping subsets recursively []. Each subset is represented by a node of the tree. The most common approach of tree construction uses binary splits; that is, at each iteration at most two child nodes, left ($\leftharpoonup$) and right ($\rightharpoonup$), are grown from a parent node. Let $\mathcal{N}\subseteq\mathcal{D}$ be the subset denoted by the parent node, and $\mathcal{I}_{\mathcal{N}}\subseteq\mathcal{I}_{\mathcal{D}}$ indicate the indices of the samples in that node. To obtain a pair of child nodes, a feature $j$ and a split value $v$ is selected to obtain the left and right index sets

$$ $\mathcal{I}^{\leftharpoonup}_{\mathcal{N}}(j,v)=\{i\in\mathcal{I}_{\mathcal{N} }:x_{ij}\leq v\}$ $$

and

$$ $\mathcal{I}^{\rightharpoonup}_{\mathcal{N}}(j,v)=\{i\in\mathcal{I}_{\mathcal{N }}:x_{ij}>v\},$ $$

respectively. Then, the gain by such a split is evaluated by

$$ $\alpha^{\leftharpoonup}_{\mathcal{N}}(j,v)\sum_{i\in\mathcal{I}^{ \leftharpoonup}_{\mathcal{N}}(j,v)}\ell(\mbox{\boldmath{$y$}}_{\mathcal{N}}, \mbox{\boldmath{$y$}}_{i})+\alpha^{\rightharpoonup}_{\mathcal{N}}(j,v)\sum_{i \in\mathcal{I}^{\rightharpoonup}_{\mathcal{N}}(j,v)}\ell(\mbox{\boldmath{$y$}} _{\mathcal{N}},\mbox{\boldmath{$y$}}_{i}),$ (1) $$

where $\mbox{\boldmath{$y$}}_{\mathcal{N}}\in\mathbb{Y}$ is the prediction of the parent node, $\ell:\mathbb{R}^{k}\times\mathbb{R}^{k}\mapsto\mathbb{R}$ is the loss function, $\alpha^{\leftharpoonup}_{\mathcal{N}}(j,v)=\tfrac{|\mathcal{I}^{\leftharpoonup
}_{\mathcal{N}}(j,v)|}{|\mathcal{I}_{\mathcal{N}}|}$ as well as $\alpha^{\rightharpoonup}_{\mathcal{N}}(j,v)=1-\alpha^{\leftharpoonup}_{
\mathcal{N}}(j,v)$ are the scaling coefficients representing the proportional weights of the respective subsets. Among all possible splits, the one that achieves the largest gain is used to construct the child nodes. The splitting continues until one of the stopping conditions, such as reaching the maximum depth or the minimum number of samples in a node, is satisfied. Then, the samples in the last parent constitute a leaf node. Considering the feasibility requirements on the predictions, the output vector $\mbox{\boldmath{$y$}}_{\mathcal{N}}$ of such a leaf node should be feasible; i.e., it should belong to the feasible set, $\mbox{\boldmath{$y$}}_{\mathcal{N}}\in\mathbb{Y}$.

### 2.1 Exact Approach

The common loss functions (e.g., mean squared error, mean Poisson deviance) used in constructing binary DT regressors are separable. Moreover, the loss in the multi-target case is obtained by summing up individual losses for each target. Consider, for instance, the well-known mean squared loss for which the prediction at a node with the subset $\mathcal{N}$ is given by

$$ $\mbox{\boldmath{$y$}}_{\mathcal{N}}=\arg\min_{\mbox{\boldmath{$y$}}\in\mathbb{ Y}}~{}\frac{1}{2|\mathcal{I}_{\mathcal{N}}|}\sum_{i\in\mathcal{I}_{\mathcal{N} }}\|\mbox{\boldmath{$y$}}-\mbox{\boldmath{$y$}}_{i}\|^{2}.$ $$

If there are no restrictions on the output vectors, i.e., $\mathbb{Y}=\mathbb{R}^{k}$, then it is not difficult to observe by simple differentiation that $\mbox{\boldmath{$y$}}_{\mathcal{N}}$ is the mean of the target vectors in $\mathcal{N}$ minimizing the subset variance. This is, in fact, the output of a node used by the well-known implementations, like scikit-learn []. The same holds for the mean Poisson deviance loss function. In fact, the following theorem shows that the mean vector is still the optimal choice when the constraints on the output vectors result in a convex feasible set. The proof of this theorem follows from observing that the mean operation is equivalent to taking the convex combination of the feasible target vectors in the node. The details of this proof are given formally in Section A of the appendix. There, we also present a counter example showing that the same argument does not apply when the loss function is the mean absolute deviation and its unconstrained minimum, i.e., the median of the target vectors, is selected as the node’s output (see Remark 1 in Section A of the appendix).

###### Theorem 1

Let the feasible set $\mathbb{Y}\subseteq\mathbb{R}^{k}$ be convex. If the loss function is the mean squared error or the mean Poisson deviance, then the mean of the target vectors in a node minimizes the loss.

This theorem implies that a restricted class of feasible sets can be handled with standard DT implementations using certain loss functions. For instance, simple bounds or, more generally, linear constraints on the target vectors lead to convex feasible sets. However, the same conclusion does not hold when the feasible set is not convex. Going back to our illustrative example of inventory amount prediction, suppose for simplicity that each warehouse has the same capacity of $C$ units. If we use only these bounds as constraints, then any standard implementation returning the mean vector for each node output would not violate this constraint. However, if we also impose the budget restriction “at most two among five warehouses can be used for transportation,” then we need auxiliary binary variables to represent this restriction as a constraint, and hence, the feasible set $\mathbb{Y}$ becomes non-convex. Formally,

$$ $\mathbb{Y}=\left\{\mbox{\boldmath{$y$}}\in\mathbb{R}^{5}:\sum_{t=1}^{5}z_{t} \leq 2,\ 0\leq y_{t}\leq Cz_{t},\text{ and }\ z_{t}\in\{0,1\}\text{ for }t=1, \dots,5\right\},$ (2) $$

where the binary variable $z_{t}$ becomes one if $y_{t}$ units of the product is transported to warehouse $t$. Note that when $z_{t}$ is zero, then it also forces $y_{t}$ to become zero implying that product is not transported to the corresponding warehouse.

We can now state the formal prediction model that needs to be solved at each node:

$$ $\mbox{\boldmath{$y$}}_{\mathcal{N}}=\arg\min_{\mbox{\boldmath{$y$}}\in\mathbb{ Y}}~{}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\ell(\mbox{\boldmath{$y$}},\mbox{ \boldmath{$y$}}_{i}),$ (3) $$

where $\mathcal{I}_{\mathcal{N}}\subseteq\mathcal{I}_{\mathcal{D}}$ with $\mathcal{N}\subseteq\mathcal{D}$. Dealing with nonconvex sets can pose a challenge when solving large-scale constrained models. Fortunately, the dimension of the problem ($k$) is quite small for many applications. In case of a mixed integer linear optimization problem, like our example on warehouse inventories, both free and commercial powerful solvers are available, e.g., COIN-OR [], GUROBI []. However, it is also important to remember that this optimization problem must be solved for each pair of feature-split values at every node. This increases the computation time compared to basic DT implementation that ignores the constraints, especially for deep trees. Next, we will also propose several heuristic approaches to develop faster –albeit potentially less accurate or infeasible– alternatives.

###### Theorem 1

### 2.2 Heuristic Approaches

A common method for training machine learning models with restrictions involves relaxing these constraints and adding a term to the loss function. This term serves to penalize predictions that are not feasible. In our setting, this idea leads to the following model:

$$ $\mbox{\boldmath{$y$}}^{R}_{\mathcal{N}}=\arg\min_{\mbox{\boldmath{$y$}}\in \mathbb{R}^{k}}~{}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\ell(\mbox{\boldmath{$y$ }},\mbox{\boldmath{$y$}}_{i})+\lambda v_{\mathbb{Y}}(\mbox{\boldmath{$y$}}),$ (4) $$

where $\lambda>0$ is the predefined penalty coefficient, and $v_{\mathbb{Y}}:\mathbb{R}^{k}\mapsto\mathbb{R}_{+}$ is the penalty function that assigns a positive value to the infeasible predictions, $\mbox{\boldmath{$y$}}\notin\mathbb{Y}$. This technique is a reminiscent of the approach known as Lagrangean relaxation. However, it is important to note that Lagrangean relaxation would introduce a separate penalty coefficient, known as Lagrangean multiplier, for each constraint []. The objective function of (4) shows the trade-off between accuracy and feasibility. For instance, using the feasible region defined in (2) leads to the following splitting problem for the mean-squared loss:

$$ $\mbox{\boldmath{$y$}}^{R}_{\mathcal{N}}=\arg\min_{\mbox{\boldmath{$y$}}\in \mathbb{R}^{5},\mbox{\boldmath{$z$}}\in\{0,1\}^{5}}~{}\frac{1}{2|\mathcal{I}_{ \mathcal{N}}|}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\|\mbox{\boldmath{$y$}}- \mbox{\boldmath{$y$}}_{i}\|^{2}+\lambda\sum_{t=1}^{5}\big{(}(z_{t}-2)+(y_{t}- Cz_{t})\big{)}.$ (5) $$

Here, $\lambda$ in (4) plays the role of balancing optimality versus feasibility. As $\lambda$ increases, the optimization algorithm tends to favor feasible solutions over accurate ones. Conversely, for small values of $\lambda$, accuracy is prioritized. Therefore, one is not guaranteed to obtain a feasible solution with the relaxation approach. Additionally, depending on the structure of the penalty function $v_{\mathbb{Y}}$, problem (4) can be difficult to optimize. Note that tuning $\lambda$ requires attention as the optimization problem changes at each node. Attempting to adjust a different penalty coefficient for each node could quickly become a computationally formidable task.

Our second heuristic is based on observing that all samples in a node belong to the feasible set. Thus, selecting one of them results in a feasible prediction. Considering the spread of the output values in the subset, a natural candidate could be the medoid as it would be a representative sample within the cluster designated by the node. The medoid of the subset is obtained by solving

$$ $\mbox{\boldmath{$y$}}^{M}_{\mathcal{N}}=\arg\min_{\mbox{\boldmath{$y$}}_{i^{ \prime}}\in\mathcal{N}}~{}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\ell(\mbox{ \boldmath{$y$}}_{i^{\prime}},\mbox{\boldmath{$y$}}_{i}).$ (6) $$

The solution to this problem involves evaluating the pairwise loss values and then summing these values for each sample. As a result, the medoid of a subset of samples can be found very quickly. In case there are multiple samples that could be the medoids, one of them can be selected arbitrarily. Clearly, the medoid at each node is a feasible prediction but not necessarily the optimal solution of (3). At this point, we note that a tree can be grown very deep so that each leaf node contains a single sample. While this allows for feasible predictions, such deep trees are susceptible to overfitting and can quickly lose interpretability (see also Figure 1).

Similar to the exact approach of the previous subsection, both heuristic approaches above are applied to every feature and split value at each node. An alternative heuristic would be growing a tree as in basic implementation without considering the constraints. Then, the feasibility of DT’s prediction can be guaranteed by solving the exact model (3), the medoid model (6), or the relaxation model (4) –with an appropriate choice of the penalty coefficient– at the *leaf nodes only*. This approach would improve the training time significantly.

## 3 Computational Study

This section is dedicated to testing the performance of our exact and heuristic approaches proposed for training output-constrained decision trees (OCDTs). We obtain different implementation variants by using one of the proposed approaches at a decision node for splitting (S) or at a leaf (L) node for prediction. We also use the shorthand notation O, R, M, A for the exact *optimization*, *relaxation*, *medoid*, and *average*(^1^11This refers to the standard implementation of taking average of the target vectors of the samples in a node.) approaches, respectively. For instance, the variant that uses relaxation approach for splitting and the medoid approach for prediction is denoted by S(R):L(M). Likewise, S(A):L(A) simply refers to standard DT implementation in software packages that ignores the constraints. In our study, we have used the scikit-learn implementation only for this variant [].

Our objective in this section is two-fold: First, we aim to assess the effectiveness of our proposed approaches on an existing dataset. Second, we aim to evaluate their performance in a controlled setting using an extended version of the same dataset. Our computational experiments can be reproduced by using the dedicated GitHub repo (^2^22[https://github.com/sibirbil/OCDT](https://github.com/sibirbil/OCDT)). For all methods compared, we experimented with varying maximum depth values of 3, 6, 9, 12, and 15, while ensuring a minimum of 10 data points to split a node and at least 5 data points in each leaf. The selection of penalty coefficients, $\lambda$, for each variant can be guided by evaluating several values over a reasonable range through cross-validation. Ideally, the chosen coefficient should minimize MSE while ensuring that no infeasible predictions are generated. The results were obtained using a computer equipped with an Intel(R) Core(TM) i7-9750H CPU 2.60GHz processor, featuring 12 CPUs, and 64GB of RAM.

### 3.1 Score Dataset

As a part of our computational study, we employ the scores dataset available online []. This dataset consists of 1000 records, simulating various personal and socio-economic factors of students at a public school. The dataset features three target variables, representing math, reading, and writing scores. To tailor the data to our research objectives, we have imposed two constraints: First, if the reading score falls below 50, the writing score is set to zero. Second, if the combined reading and writing scores amount to less than 110, the math score is also set to zero. The full description of our mathematical model is available in Section B of the appendix.

We conducted a thorough five-fold cross-validation experiment on the dataset using five different tree-depth parameter values, employing a total of nine distinct variants which are trained with the same parameters. These variants are classified into two primary components: split evaluation and leaf predictions. Split evaluation entails the assessment of each candidate for binary splitting based on the prediction generated by a specific method. Conversely, leaf prediction refers to the prediction method employed by each leaf node. For instance, the methodology involving the medoid split evaluation and optimal leaf prediction entails that each split during tree growth is determined by evaluating the total mean squared error (MSE) values of the resulting child nodes using their medoids as predictions. Subsequently, the resulting DT solves an optimization problem specific to the dataset for data falling to each leaf to assign the prediction of that particular node. The penalty coefficient $\lambda$ is selected by monitoring the trade-off between feasibility and accuracy through five-fold cross-validation.

We evaluated the performance of each variant based on the MSE value obtained during the cross-validation process. Our analysis focused on comparing the performance of an OCDT variant with different split evaluation and leaf prediction approaches against the baseline DT implementation (scikit-learn). The results are visualized in the left side of Figure 2 illustrating the key metrics concerning the performance of the variants. The bar plot reflects the (scaled) average MSE values with the left axis indicating the corresponding values. On top of each bar, the MSE values are displayed alongside the average number of infeasibilities resulting from the methods’ predictions (in parenthesis). The scatter plot shows the average training time of each method with the right axis denoting the corresponding values.

In general, it appears that it is important to assess the performance of each variant individually. We find that the variant utilizing relaxation in the optimization problem in both split evaluation and leaf prediction, S(R):L(R) exhibits the worst MSE score. However, the variant employing relaxation in split evaluation and optimization in leaf prediction, S(R):L(O) achieves the best MSE score without violating any constraints. Overall, we observe that the variants solving optimization problems at the leaf nodes for predictions demonstrate the top three performances for this dataset with no constraint violations. This affirms both the efficiency of the proposed variants and their adherence to constraints. We note that the baseline scikit-learn implementation, S(A):L(A) violated numerous constraints (averaging 16.4 infeasible predictions), aligning with our earlier findings.

Figure: Figure 2: Performance metrics plot (left) and the Friedman test for statistical differences (right) on the scores dataset. The values on each bar represent the average MSE value and the average number of infeasible predictions (in parenthesis) on the test sets. The scatter plot shows the training time for each variant. Note that MSE and computation times are scaled according to the maximum value for the corresponding performance measure obtained by one of the variants.
Refer to caption: extracted/5895402/scores_plot_v2.png

The comparison of the proposed variants for statistical significance is conducted using the Friedman test [], and if a significant difference is detected by the Friedman test, it is followed by the Nemenyi test []. The right side of Figure 2 shows the test results. We observe that the $p$-value is practically zero indicating that the null hypothesis of identical variants can be rejected at the 0.95 confidence level. This means that there is a statistically significant difference in the performances of the variants with S(R):L(O) having a notably better performance when compared against the others.

### 3.2 Extended Datasets

In this section, we extended the previous dataset [] to mimic the motivating inventory control example that we have given in Section 1. To test the proposed approaches in this controlled environment, we imposed other restrictions on the dataset, like each target vector having at most one nonzero value. The details of our data generation process along with the underlying mathematical model is given in Section C of the appendix. We call this dataset class.

In accordance with the experiments in Section 3.1, we apply the same methodology with five-fold cross-validation. This is done to compare our nine variants using the same parameters. To leverage the advantage of data generation, we evaluated our proposed variants across nine different datasets, varying in size ($|\mathcal{I}_{\mathcal{D}}|$) and the number of targets ($k$). These datasets comprised of combinations of 500, 1000, and 2000 samples with total number of target values of 3, 5 and 7.

The MSE score and computation time are depicted on the left side of Figure 3, where the average metrics across all nine datasets are given. On the right side of Figure 3, we see that $p$-value is close to 0, which indicates that the methods have statistically different performances at the 0.95 confidence level. We observe that none of the proposed algorithms produce infeasible predictions, unlike the baseline scikit-learn method, S(A):L(A) which, despite yielding the lowest MSE among all approaches, generated the only positive average number of infeasible predictions. This suggests that the scikit-learn method compromised constraints to improve the prediction performance, a trade-off that is clearly not suitable for applications.

Figure: Figure 3: Performance metrics plot (left) and the Friedman test for statistical differences (right) on class datasets. The values on each bar of the left plot represent the average MSE value and the average number of infeasible predictions (in paranthesis) on the test sets. The scatter plot shows the training time for each variant. Note that MSE and computation times are scaled according to the maximum value for the corresponding performance measure obtained by one of the variants.
Refer to caption: extracted/5895402/class_plot_v2.png

The box plot in Figure 4 illustrates the MSE values of all the variants across nine datasets with different sizes in terms of sample-sizes and targets. The proposed variants show comparable MSE scores, with most of the methods achieving scores within a close range, as indicated by the overlapping interquartile ranges. While there is some variation, the results suggest that none of the methods stand out as consistently superior in terms of accuracy.

Figure: Figure 4: MSE box plot on the class datasets.
Refer to caption: x1.png

Figure 5 depicts that the training times of the variants vary significantly with those involving optimization and relaxation, such as S(R):L(R) and S(R):L(M), taking notably longer. This observation aligns with the computational complexity expected from optimization-based approaches. Another noteworthy observation is that the methods with relatively short training times, likeS(A):L(M), tend to have higher MSE values. This trade-off between model complexity and training time versus prediction accuracy is evident across the board. Furthermore, variants like S(O):L(M) strike a balance, achieving moderate MSE values while maintaining competitive training times, suggesting a promising option for practical use in time-sensitive applications.

Figure: Figure 5: Training time box plot on the class datasets.
Refer to caption: x2.png

Figure 6 illustrates a schematic representation of a tree trained on a class dataset. This tree structure is obtained by solving the optimization problem for the split evaluation as well as for the leaf prediction; i.e., the selected variant is S(O):L(O). In this figure, the green nodes signify the leaf nodes, while the red nodes represent the decision nodes where the data is split. Additionally, the splits after each decision node are depicted with arrows indicating the left and the right child nodes, accompanied by text specifying the split feature and threshold. It is easy to observe on the tree plot that the OCDT made splits are both accurate and feasible. As we outline in Section C of the appendix, the target values are generated in a way that they are dependent on Grade and GradeAvgPrevElec. This is clearly depicted in the tree plot as the majority of splits employed by the OCDT do use those features. Additionally, the second-level splits in the tree plot are dependent on feature Class being one or higher, which aligns with the conditions imposing constraints to assign some targets to zero. This visualization example shows that OCDTs have the potential to be interpretable like their conventional counterparts available in standard DT implementations.

Figure: Figure 6: Visualization of an OCDT model on a class dataset.
Refer to caption: extracted/5895402/ocrt_tree_class_updated_2.png

Our computational study shows that the proposed approach is highly effective in meeting constraints while delivering strong performance on various datasets. By incorporating optimization in split evaluations and leaf predictions, our methodology offers more accurate and feasible predictions compared to traditional methods, even on complex datasets. These results emphasize the usability and versatility of our approach, particularly in situations where compliance with constraints is essential for success.

## 4 Conclusion

We have introduced new methods to handle constraints among target variables in decision trees, an area yet unexplored in the existing literature. We have proposed an optimization-based exact approach and several heuristic approaches to ensure feasible predictions and maintaining computational efficiency. Through a computational study with existing and generated datasets, we have demonstrated the effectiveness of the proposed methods in producing feasible predictions, showcasing their potential for applications in multi-target learning.

#### Limitations.

The work presented here has a few limitations. First, the optimization-based exact method can be computationally demanding for deep trees, which may impede its application for larger datasets. Second, the heuristic approaches proposed may not always produce accurate or even feasible predictions, limiting their reliability. Third, when the constraints among targets are highly non-convex, solving the prediction problem at each node can prove challenging, potentially affecting the overall performance of the decision tree. Finally, in its current form, it is unclear how these output-constrained decision trees can serve as the base learners for ensemble methods.

#### Future research.

We aim to investigate the integration of our approach into bagging methods like random forests, as trees often act as base learners in these methods. We are also interested in exploring a warm-start method in the exact approach, where an optimization problem needs solving at each node. Given that these nodes share a parent-child relationship, the subsets they create could aid in devising this method for subsequent problems at the child nodes once the parent node is resolved. We are also considering the inclusion of target constraints in optimization-based algorithms for training trees like, for example, optimal classification trees []. Although these algorithms can currently handle smaller problem sizes compared to standard implementations like CART [], this is a promising research direction. We are also looking into the potential advantages of considering constraints in clustering, since training a decision tree can be seen as a supervised clustering method.

## Appendix A Proof of Theorem 1

Theorem 1
Let the feasible set $\mathbb{Y}\subseteq\mathbb{R}^{k}$ be convex. If the loss function is the mean squared error or the mean Poisson deviance, then the mean of the target vectors in a node minimizes the loss.

###### Proof.

Suppose that the node corresponds to the subset $\mathcal{N}\subseteq\mathcal{D}$. Since $\mathbb{Y}\subseteq\mathbb{R}^{k}$, we have

$$ $\min_{\mbox{\boldmath{$y$}}\in\mathbb{R}^{k}}~{}\sum_{i\in\mathcal{I}_{ \mathcal{N}}}\ell(\mbox{\boldmath{$y$}},\mbox{\boldmath{$y$}}_{i})\leq\min_{ \mbox{\boldmath{$y$}}\in\mathbb{Y}}~{}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\ell (\mbox{\boldmath{$y$}},\mbox{\boldmath{$y$}}_{i})$ $$

for any loss function $\ell:\mathbb{R}^{k}\times\mathbb{R}^{k}\mapsto\mathbb{R}$. We call the optimization model on the left as the unconstrained problem and the one on the right as the constrained problem. Note that if the optimal solution of the unconstrained problem is in the feasible set $\mathbb{Y}$, then it should also be the optimal solution of the constrained problem. We use the optimal solution of the unconstrained problem as the output vector of the node and refer to it as $\mbox{\boldmath{$y$}}_{\mathcal{N}}$. Both loss functions that we consider are separable in terms of the target vector. Thus, we write $\ell(\mbox{\boldmath{$y$}},\mbox{\boldmath{$y$}}_{i})=\sum_{t\in\mathcal{T}}
\ell(y_{t},y_{it})$, where $\mathcal{T}=\{1,\dots,k\}$ is the set of target indices.

The unconstrained problems with the mean squared error and the mean Poisson deviance loss functions are given by

$$ $\min_{\mbox{\boldmath{$y$}}\in\mathbb{R}^{k}}~{}\frac{1}{2|\mathcal{I}_{ \mathcal{N}}|}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\sum_{t\in\mathcal{T}}(y_{t} -y_{it})^{2}$ $$

and

$$ $\min_{\mbox{\boldmath{$y$}}\in\mathbb{R}^{k}}~{}\frac{1}{|\mathcal{I}_{ \mathcal{N}}|}\sum_{i\in\mathcal{I}_{\mathcal{N}}}\sum_{t\in\mathcal{T}}(y_{it }\ln(\frac{y_{it}}{y_{t}})+y_{t}-y_{it}),$ $$

respectively. Clearly, both loss functions are differentiable and convex. Thus, the necessary-and-sufficient optimality conditions for both problems are given by

$$ $\sum_{i\in\mathcal{I}_{\mathcal{N}}}\frac{\partial\ell(y_{t},y_{it})}{\partial y _{t}}=0,\ \ t\in\mathcal{T}.$ (7) $$

Solving next the conditions in (7) for both unconstrained problems leads to

$$ $\sum_{i\in\mathcal{I}_{\mathcal{N}}}\frac{\partial\ell(y_{t},y_{it})}{\partial y _{t}}=0,\ t\in\mathcal{T}\implies y_{t}=\frac{1}{|\mathcal{I}_{\mathcal{N}}|} \sum_{i\in\mathcal{I}_{\mathcal{N}}}y_{it},\ t\in\mathcal{T}\implies\mbox{ \boldmath{$y$}}_{\mathcal{N}}=\frac{1}{|\mathcal{I}_{\mathcal{N}}|}\sum_{i\in \mathcal{I}_{\mathcal{N}}}\mbox{\boldmath{$y$}}_{i}.$ $$

Thus, the optimal solutions is the mean vector. Because $\mbox{\boldmath{$y$}}_{i}\in\mathbb{Y}$ for $i\in\mathcal{I}_{\mathcal{N}}$ and the feasible set $\mathbb{Y}$ is convex, the convex combination of $\mbox{\boldmath{$y$}}_{i}$ vectors should also be feasible; that is, $\mbox{\boldmath{$y$}}_{\mathcal{N}}\in\mathbb{Y}$. Therefore, the mean vector is also the optimal solution of the constrained problem. This shows the desired result.
∎

###### Remark 1

Standard multi-target regression tree implementations, like that of scikit-learn, use the median of the target vectors in the node as prediction when the loss function is selected as the mean absolute deviation. However, our result in Theorem 1 does not hold in this case, because the median vector may lie outside the feasible set even if it is convex.

Figure 7 shows such an example for $k=2$, where the *convex* feasible set is shown by the shaded area, and the node has five feasible target vectors denoted by the black circles. The median vector, shown as the red circle, is clearly outside the feasible set.

Figure: Figure 7: A counter example showing that the median vector lies outside the convex feasible set.
Refer to caption: extracted/5895402/medianexample.png

###### Proof.

###### Remark 1

## Appendix B Optimization Model of The scores Dataset

The optimization model imposed on the scores dataset is formulated as follows:

$$ minimize $\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\sum_{t=1}^{k}({y}_{t}-y_{it})^{2}$ subject to $\displaystyle{y}_{t}\leq 110b_{t},\quad\quad t=1,\dots,k,\quad$ $\displaystyle{y}_{2}\geq 50b_{3},\quad$ $\displaystyle{y}_{1}+{y}_{2}\geq 110b_{1},$ $$

where $n$ is the number of instances, $k$ is the number of targets, ${y}_{t}$ is the prediction for the $t$-th target, $y_{it}$ is the true value of the $t$-th target for the $i$-th instance, and $b_{t}$ is a binary variable indicating whether the prediction for the $t$-th target is used ($b_{t}=1$) or not ($b_{t}=0$).

## Appendix C Data Generation and Optimization Model of The class Dataset

Table 1 outlines the distributions and their parameters used to generate the extended data, illustrating how the features of the dataset are created. We should note that the feature GradePerm consists of the permuted values of the feature Grade, created with the purpose of imposing randomness on the data. The target values are computed by taking the equally weighted sum of the Grade and GradeAvgPrevElec, with the addition of a uniform noise component ranging between -10 and 10 if GradeAvgPrevElec is non-zero. If GradeAvgPrevElec is zero, the target value is determined by Grade plus a normal noise with a mean of 10.0 and a standard deviation of 2.0. These computed values are then clipped to ensure they fall within the range of 0 to 100. Additionally, certain target values are deliberately set to zero based on random selection according to the features. For instance, if the Class feature is equal to one, all but a randomly chosen subset of the first two targets’ grades are set to zero. To fit the data to our research objectives, we have imposed two constraints on the created dataset. One constraint is created using binary decision variables to ensure that one non-zero prediction is generated at most. Another set of constraints has been applied to each target variable to limit the prediction based on the binary variable.

**Table 1: Features of the class dataset**
| Features | Distribution | Min | Max | Mean | Std.Dev. |
| --- | --- | --- | --- | --- | --- |
| EnrolledElectiveBefore | Discrete Uniform | 0 | 1 | - | - |
| GradeAvgPrevElec | Normal | 0 | 100 | 60 | 15 |
| Grade | Normal | 0 | 100 | 70 | 10 |
| Major | Discrete Uniform | 1 | 3 | - | - |
| Class | Discrete Uniform | 1 | 4 | - | - |
| GradePerm | Normal | 0 | 100 | 70 | 10 |

The optimization model is as follows:

$$ minimize $\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\sum_{t=1}^{k}({y}_{t}-y_{it})^{2}$ subject to $\displaystyle\sum_{j=1}^{k}b_{j}\leq 1,\quad$ $\displaystyle{y}_{t}\leq 110b_{t},\quad$ $$

where $n$ is the number of instances, $k$ is the number of targets, ${y}_{t}$ is the prediction for the $t$-th target, $y_{it}$ is the true value of the $t$-th target for the $i$-th instance, and $b_{t}$ is a binary variable indicating whether the prediction for the $t$-th target is used ($b_{t}=1$) or not ($b_{t}=0$).

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: Besides comments in Weaknesses section, I also list several questions here. 
1. The experiments are based on artificial constraints applied to a single dataset. Could the authors suggest specific real-world datasets where such constraints naturally occur? Additionally, could the authors discuss how their approaches might be adapted or evaluated in those contexts?
2. Since the dataset was modified to include specific constraints, it is unclear how the insights derived from Figure 1, particularly regarding the trade-off, remain reasonable and reliable. Adjusting a dataset with imposed constraints can introduce biases or alter its original distribution, potentially influencing the model’s behaviours. It would be beneficial for the authors to discuss methods for mitigating these potential biases to ensure the robustness of their findings or provide a comparative analysis between the original and modified datasets, showing how the distribution of key variables changed.
3. Could you clarify which specific experiment led to the results in Figure 1? In particular, it is not immediately clear how the 4000 samples were generated. It would be better to provide a specific section reference where the experiment details should be added or include a brief description of the experimental setup directly in the figure caption.
soundness: 2
strengths: 1. The paper poses an insightful and valuable research question regarding constrained multi-target learning in decision trees and presents effective approaches to address it.
2. The paper explain the problem and the methods clearly and intuitively.
summary: This paper investigates multi-target learning problems where constraints are imposed on the target variables. To address this challenge, the authors propose an exact optimization-based approach that enforces these constraints during the decision tree-building process, ensuring feasibility but at the cost of increased computational complexity. Additionally, they introduce two heuristic methods that offer more efficient alternatives. Through experiments on datasets with imposed constraints, the authors demonstrate the effectiveness of these approaches in maintaining feasible predictions under constrained settings.
weaknesses: 1. The paper would benefit from a more robust differentiation of its contributions from existing research, particularly through a deeper exploration of recent work. The Related Work section includes only a few papers from the last five years, which limits the context for evaluating the novelty of this approach. Additionally, the claim that this is the first method to address target-constrained problems may be overstated. Expanding on prior work in this area would help to clarify the unique aspects and positioning of this contribution.
2. In the SCORE dataset, two specific constraints were imposed on the target variables, but the rationale behind selecting these particular constraints is unclear. Could the authors provide justification for the appropriateness of these constraints for this dataset? If the constraints are too tailored to the proposed method, this may obscure potential limitations. 
3. It would be clear if the paper could provide intuition behind the two heuristic methods and explicitly highlight the benefit of each heuristic from a design perspective.
4. The reliance on a self-modified dataset and a baseline comparison with only the traditional decision tree limits the strength of the paper’s claims regarding novelty and contributions. Have the authors considered testing the proposed methods against state-of-the-art decision tree models or advanced tree-based methods that incorporate constraints? Broader comparisons with such models would provide a stronger foundation for evaluating the effectiveness and competitive advantage of the proposed approaches.
5. The paper lacks a detailed analysis of the computational advantages of the two heuristics. It would be helpful to specify the order of magnitude of the improvement and to quantify how many seconds the heuristics take relative to a traditional decision tree. Providing these details would clarify the efficiency gains of the heuristics.

---

confidence: 5
contribution: 2
initial_rating: 5
presentation: 2
questions: N/A
soundness: 2
strengths: 1. The proposed algorithm addresses an interesting and practical problem of a regression with constraints on outputs.
2. The authors propose heuristics to alleviate the computational time of MIO.
summary: The authors propose an algorithm for regression problems with constraints on tree outputs. Their approach modifies the greedy recursive partitioning algorithm to grow a binary tree until a specified maximum depth or minimum number of samples in a node is reached. To split a node, the algorithm uses a mixed integer optimization (MIO) solver to find a feasible vector that minimizes the mean squared error loss. To reduce computation time, the authors introduce two heuristics: (1) Lagrangian relaxation and (2) selecting a feasible medoid as an approximation. The experiments are conducted on a single dataset, modified to evaluate the algorithm under various regimes (e.g., different number of constraints and dataset sizes). The method is compared against variations of itself, differing in the types of node optimization (e.g., with or without relaxation), and a baseline CART-like tree where outputs are replaced with feasible values.
weaknesses: 1. The use of an MIO solver within the tree induction process raises concerns about practicality in real-world applications, as the MIO problem is NP-hard and may not terminate in a reasonable amount of time.
2. The training times are reported in scaled units, which makes it difficult to see the true computational cost of the algorithm.
3. The paper uses a single dataset, making it difficult to assess performance relative to a straightforward baseline (e.g., training a CART regression tree and substituting outputs with feasible values).

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: No questions.
soundness: 2
strengths: The proposed heuristic approaches are interesting though not thoroughly explored or analyzed.
summary: The paper considers the problem of learning decision trees where the leaf nodes predict (constant) vector-valued output with some constraint on them, for example $\ell_0$ constraint. It adapts the traditional CART-type greedy top-down induction method and when evaluating each possible split (feature index and threshold pair), it takes into account the constraints. The split evaluation is either solved exaclty (e.g. using mixed-integer optimization solvers), or approximately using a penalty method or using a best training point as a solution. The proposed method has been evaluated on smaller-scale problems.
weaknesses: 1. Unprecise problem formulation. The paper discusses constraints on leaf predictions in general without specifying the exact type of constraints considered. Some constraints are much harder (or easier) to solve than others. 
2. Limited algorithmic novelty. The proposed method is still based on the traditional greedy top-down induction trees, and just adapts the split evaluation to take into account the constraints.
3. Lack of comparison. Methods that attempt to learn exactly (optimally) decision trees can straightforwardly incorporate the constraints in the objective function (and perhaps this makes the search space smaller and accelerates those type of algorithms). Given that the paper considers smaller-scale problems (1000s training points), those optimal decision tree algorithms can be trained for these problems, and thus be compared.
4. The literature review does not consider methods on fairness where different type of constraints are used for learning decision trees. It also does not consider other methods for learning trees such as dynamic programming, MIO and alternating optimization.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: Q/S1: Do you assume the constraints are imposed at the instance level, i.e., different instances can have different sets of feasible outputs, or at the population level, i.e., all the instances should have the same set of feasible outputs? Please refer to W1 for my further comments on this point.

Q/S2: It might be beneficial to show the (average) scores on the cases where the DTs provide feasible and infeasible solutions separately. 

Q/S3: While extensions to ensemble methods and adaptations to multi-label and multi-dimensional classification problems seem to be not obvious, adding more discussions on potential directions would be appreciated.
soundness: 3
strengths: S1: I think learning and predicting under output constraints are interesting problems. 

S2: The problem of finding a representative output vector, which is an essential step of the proposed framework, is elaborated both exactly and approximately/heuristically.

S3: Experiments are conducted to assess the potential advantages of the proposed framework.
summary: As far as I understand, this submission proposes a framework for multivariate regression trees under output constraints. Under this setting, the representative output vector of each node or leaf should belong to some given feasible set, and should be as accurate as possible with respect to the evaluation loss. I think this setting is interesting and can be useful in different applications, provided that the constraints, that shape the feasible sets, are well-defined. 

The submission elaborates on the problem of finding a representative output vector both exactly and approximately/heuristically. Experiments on two data sets are conducted, where the loss is the mean squared error (MSE), to assess the potential advantages of the proposed framework.
weaknesses: W1: The problem statement is not entirely clear. In particular, it might be beneficial to further clarify whether the constraints are imposed at the instance level, i.e., different instances can have different sets of feasible outputs, or at the population level, i.e., all the instances should have the same set of feasible outputs. I think neither of the cases calls for further elaboration/clarification on the prediction phase. However, the earlier case might call for further elaboration/clarification on the training phase.

W2: I think more insight discussions about the potential (dis)advantages of the proposed framework might be beneficial. The author(s) might consider adding the (average) scores on the cases DTs provide feasible and infeasible solutions separately. Ideally, one might wish to see frameworks/methods, which provide promising scores on both parts, but a good trade-off between the two parts might also be appreciated. 

W3: More discussions on extensions to ensemble methods and adaptations to multi-label and multi-dimensional classification problems might be beneficial. For example, while the submission mentions the label powerset algorithm as a solution for multi-label and multi-dimensional classification problems, handling a possibly large multi-class classification problem, which can be both data-hungry and contain output constraints at the instance level might not be obvious.

Meta-review: The paper concerns the problem of multi-target learning in which constraints are imposed on target variables. The Authors introduce a variant of the decision tree algorithm to solve this problem. 

The Reviewers highlighted several weaknesses in the submission, including an imprecise problem definition, a limited discussion of related work, insufficient evidence regarding computational costs, and a lack of extensive empirical studies.

All ratings are below the bar and the Authors did not prepare any rebuttal.

Decision: Reject