Optimal Decision Trees For Interpretable Clustering with Constraints (Extended Version)

Table 1 describes the set of Boolean variables.

We encode the construction of the decision tree (Eqs. (7)-(15)) following Shati et al. Shati et al. (2021), a state-of-the-art SAT-based decision-tree classifier that naturally supports numeric features that are prevalent in clustering problems. These clauses guarantee that exactly one feature is chosen at each branching node (Eqs. (7)-(8)), the points are directed to the left or right child of each branching node based on their value of the chosen feature (Eqs. (9)-(10)), the appearance of points at leaves correctly correspond to the path that they are directed through within the tree (Eqs. (11)-(13)), and thresholds are non-trivial (Eqs. (14)-(15)).³³3The set $O_{j}(X)$ contains all consecutive pairs of points when ordered according to feature $j$ and the set $A_{l}(t)$ ($A_{r}(t)$) contains all nodes that have $t$ as descendent of their left (right) child.⁴⁴4For a detailed description of the clauses in Eqs. (7)-(15), we refer the reader to Shati et al. Shati et al. (2021).

The following clauses extend the basic decision tree encoding to support $\epsilon$-optimal clustering trees that satisfy the $ML$ and $CL$ constraints. The clauses in Eq. (16) guarantee well-formed unary encoding of cluster labels in each leaf.⁵⁵5A unary encoding is well-formed if it does not include the sequence 01 at any point.

Eqs. (17)-(18) guarantee that the label assigned to each data point matches the label of leaf the data point reaches.

Eqs. (19)-(20) break the ties between the equivalent cluster assignments. We consider two cluster assignments $\Theta^{1}$ and $\Theta^{2}$ equivalent if there exists a relabelling, $\Gamma:K\rightarrow K$, such that $\Gamma\circ\Theta^{1}=\Theta^{2}$. To break ties, we force each $x_{i}$ in the ascending order to be assigned to the first empty cluster, if it needs a new one. This guarantees that there are no two feasible solutions that are relabellings of each other. Note that we do not eliminate viable solutions, since any clustering can be renamed into an equivalent one that respects this property.

The clause in Eq. (21) alongside the tie-breaking ones guarantee that all clusters are non-empty, assuming $|X|\geq K$. If we substitute $k$ with $k^{\prime}<k$ in Eq. (21), $k^{\prime}$ clusters are guaranteed to be non-empty. Thus, we can enforce minimum $k^{\prime}$ and maximum $k$ clusters. In our experiments we focus on the setting where all clusters are non-empty.

Eqs. (22)-(24), namely the unconditional separating clauses, guarantee that pairs in $CL$ are clustered separately.

Eqs. (25)-(26), namely the unconditional co-clustering clauses, guarantee that pairs in $ML$ are clustered together.

Eqs. (27)-(29), namely the conditional separating clauses, guarantee that a $b^{-}_{w}$ variable being set to false forces the pairs in distance class $w$ to be clustered separately.

Eqs. (30)-(31), namely the conditional co-clustering clauses, guarantee that a $b^{+}_{w}$ variable being set to true forces the pairs in distance class $w$ to be clustered together.

Eqs. (32)-(34) guarantee that the distance classes are partitioned according to valid $\lambda^{+}$ and $\lambda^{-}$ values, i.e., that $D_{1}..D_{\lambda^{+}}$ distance classes are clustered together, and $D_{\lambda^{-}{+}1}..D_{\mu}$ distance classes are clustered separately.

Since each pair of points has corresponding clauses that handle being clustered together or separately, a naive encoding will be quadratic in the number of points $|X|$. This is significant since all the other parts of the encoding are linear in $|X|$. To reduce the number of clauses, we see the set of points as nodes in a graph and exploit connections between pairs. Pairs that are forced to be clustered together are represented by positive edges and pairs that are force to be clustered separately by negative edges. The positive edges imply a set of connected components of points that will be clustered together while a negative edge between nodes in different components indicates that each of the components are mutually exclusive, i.e., each component will be in a different cluster. The order $<^{*}$, which sorts the pairs based on distance and breaks ties arbitrary, is used to greedily build the set of positive edges ($E^{+}$) and the set of negative edges ($E^{-}$). As we incrementally build the sets, each new edge can be classified as inner or crossing, based on the previous members in $E^{+}$ and $E^{-}$, to help us detect infeasibility or redundancy of clauses.

Since $E^{+}$ ($E^{-}$) represents the set of pairs that are forced to be clustered together (separately), an inner pair is forced to have the same label and a crossing pair to have different labels. We make use of this fact to avoid adding co-clustering clauses Eqs. (25, 26, 30, 31) or separating clauses Eqs. (22, 23, 24, 27, 28, 29) when it is redundant to do so. Furthermore, we can conclude infeasibility when an inner (crossing) pair is forced to be clustered separately (together).

Our treatment of the unconditional clauses Eqs. (22-26) and of the conditional ones Eqs. (27-31) differ in two ways.

1. 1.

   For the unconditional clauses the $E^{+}$ and $E^{-}$ sets are respectively constructed from $ML$ and $CL$ links. For conditional ones however, we also include the pairs that are implied to be co-clustered (separated), due to the order of distance imposed by $\lambda^{+}$ and $\lambda^{-}$ values. Note that implied pairs for the processing of conditional co-clustering clauses are not valid for the processing of conditional separating clauses.

2. 2.

   If the unconditional separation or co-clustering of a pair is infeasible, the problem is infeasible. However, for conditional clauses, we only fix the corresponding $b^{+}_{w}$ or $b^{-}_{w}$ variable to satisfy the conditional clause.

The detailed procedure is described in the appendix.

In Section 3.1 we established that in order to find an $\epsilon$-approximation of an MS-MD Pareto optimal solution, we need to find a solution that is Pareto optimal with regards to maximizing $\lambda^{+}$ and minimizing $\lambda^{-}$. Similarly, for an $\epsilon$-approximation of an MD optimal solution we need to find a solution that minimizes $\lambda^{-}$. Note that $\lambda^{+}$ and $\lambda^{-}$ are the number of $b^{+}_{w}$ and $b^{-}_{w}$ variables set to true, respectively.

To minimize $\lambda^{-}$, we simply introduce a soft clause with unit weight for each $\lnot b^{-}_{w}$ (Eq. (37)). To obtain a Pareto-optimal solution, we can optimize any function that is increasing w.r.t $\lambda^{-}$ and decreasing w.r.t $\lambda^{+}$ as objective. We opt for minimizing the simple combination of $\lambda^{-}-\lambda^{+}$ and model it using the soft clauses in Eq. (37) and Eq. (38).

For the MD objective, the $b^{+}_{w}$ variables and the clauses in Eqs. (30)-(31) become irrelevant and can be removed.

We run experiments on seven real datasets from the UCI repository Dua and Graff (2017) and four synthetic datasets from FCPS Ultsch and Lötsch (2020). The datasets vary in size, number of features, and number of clusters, as presented in Table 2. For all datasets, we normalize the values of each feature in the range $[0,100]$ so that features with larger values do not dominate the pairwise distances.

We evaluate the performance of our approach for different number of constraints, relative to the dataset size. Specifically, for a given $\kappa$ value with $0\leq\kappa\leq\frac{|X|-1}{2}$, we generate a set of $\kappa\cdot|X|$ random clustering constraints following the methodology in Wagstaff and Cardie Wagstaff and Cardie (2000): (1) We generate $\kappa\cdot|X|$ pairs of data points without repetition; (2) For each pair, we generate a ML constraint if both data point share the same ground-truth label and a CL constraint otherwise.

We evaluate the quality of the obtained clusterings based on ground-truth labels using two well-known external clustering evaluation metrics: the Adjusted Rand Index (ARI) Hubert and Arabie (1985) and the Normalized Mutual Information (NMI) Strehl and Ghosh (2002).

We compare our approach with optimal constrained clustering formulation that is not restricted to conform to a decision tree. To do so, we remove the tree-related components of the encoding, namely the variables $[a_{t,j}]$, $[s_{i,t}]$, $[z_{i,t}]$, $[g_{t,c}]$ and the clauses in Eqs. (7)-(18). Instead, we introduce the clauses in Eq. (39) that guarantee a well-formed unary encoding of labels.⁶⁶6Note that these clauses are redundant in the encoding of tree clustering as the condition is already enforced for the leaves.

When used with the MD objective and $\epsilon=0.0$, this baseline is equivalent to the maximum diameter constrained clustering formulations in previous works Dao et al. (2016, 2017), i.e., it has the same set of (feasible and) optimal solutions. However, this baseline also supports $\epsilon$-approximation and the [MD,MS] objective for a fair comparison with our approach.

The most closely related work to ours is Bertsimas et al. Bertsimas et al. (2021) that presents a MIO model for (unconstrained) clustering trees. While their model does not support clustering constraints, we can extend it to incorporate such constraints. However, Bertsimas et al. note that their model does not scale well and therefore do not present any experimental results and instead focus on a heuristic procedure that does not have any solution quality guarantees and cannot naturally support clustering constraints. Due to the relevance of the MIO model to our work, we have implemented the model using the Gurobi v10 solver and extended it to support clustering constraints (detailed description of the extended MIO model can be found in the appendix). To our knowledge, this is the only approach in the literature for constructing tree-based clustering using discrete optimization.

Our experiments with the MIO baseline found that it is unable to find a feasible solution for any of the datasets for the depths specified in Table 2 across multiple runs, both in constrained and unconstrained settings. This result is consistent with Bertsimas et al.’s Bertsimas et al. (2021) observation on the limited scalability of their MIO formulation and emphasizes the strong performance of our approach.

Finally, we study the impact of the smart pairs procedure (Section 3.2) and the $\epsilon$-approximation (Section 3.1) on the performance of our approach. Unlike the $\epsilon$-approximation, smart pairs is guaranteed to not change the set of feasible solutions and the set of optimal solutions. However, since there could be multiple optimal solutions, the different encoding may still lead to an optimal solution with a slightly different ARI score. The results presented in Table 4 show that the aforementioned methods can reduce the number of clauses and runtimes without meaningful decrease in score.