Equivariance with Learned Canonicalization Functions

Symmetry-based neural networks often constrain the architecture in order to achieve invariance or equivariance to a group of transformations. In this paper, we propose an alternative that avoids this architectural constraint by learning to produce a canonical representation of the data. These canonicalization functions can readily be plugged into non-equivariant backbone architectures. We offer explicit ways to implement them for many groups of interest. We show that this approach enjoys universality while providing interpretable insights. Our main hypothesis is that learning a neural network to perform canonicalization is better than using predefined heuristics. Our results show that learning the canonicalization function indeed leads to better results and that the approach achieves excellent performance in practice.

Introduction
The problem of designing machine learning models that properly exploit the structure and symmetry of the data is becoming more important as the field is broadening its scope to more complex problems. In multiple applications, the transformations with respect to which we require a model to be invariant or equivariant are known and provide a strong inductive bias (Bronstein et al., 2021;Bogatskiy et al., 2022;van der Pol et al., 2020;Mondal et al., 2020;Celledoni et al., 2021).
As is often the case, taking a step back and drawing analogies with human cognition is fruitful here. Human pattern recognition handles some symmetries with relative ease. When data is transformed in a way that preserves its essential characteristics, we precisely know if and how we should adapt our response. One context in which this has been particularly well-studied in cognitive science is visual shape recognition. Experiments have shown that subjects can accurately make the difference between different orientations of an object and actual modifications to the structure of an object (Shepard and Metzler, 1971;Carpenter and Eisenberg, 1978).
There are multiple ways in which this could be achieved. According to Tarr and Pinker (1989), theories of invariant shape recognition broadly fall into three categories: viewpointindependent models, in which object representations depend only on invariants features, multiple-view models in which an object is represented as a set of representations corresponding to different orientations, and single-view-plus-transformation models in which an object is converted to a canonical orientation by a transformation process.
Correspondingly, similar ideas have been explored in deep learning to achieve equivariance. Models that impose equivariance through constraints in the architecture (Shawe-Taylor, 1989;Cohen and Welling, 2016; or that only use invariants as inputs (Villar et al., 2021) can be seen as belonging to the viewpoint-independent type. The multiple-view approach includes models that symmetrize the input by averaging over all the transformations or a subset of them (Manay et al., 2006;Benton et al., 2020;Yarotsky, 2022;Puny et al., 2022). By contrast, the transformation approach has seen less interest and, to our knowledge, has yet to be used to achieve exact equivariance. This is all the more surprising considering that evidence from cognitive science suggests that this approach is used in human visual cognition (Shepard and Metzler, 1971;Carpenter and Eisenberg, 1978;Hinton and Parsons, 1981). When presented with a rotated version of an original pattern, the time taken by humans to do the association is proportional to the rotation angle, which is more consistent with the hypothesis that we perform a mental rotation.
Some works (Jaderberg et al., 2015;Qi et al., 2017) have proposed to learn transformations of input to facilitate processing in a downstream task, but these approaches are closer to regularizers and provide no guarantees. Another method is to use heuristics to standardize inputs (Yüceer and Oflazer, 1993;Lowe, 2004;Aslan et al., 2022), but this approach requires significant hand-engineering and is hardly generalizable.
Present work We introduce a systematic and general method for equivariant machine learning based on learning mappings to canonical samples. We hypothesize that among all valid canonicalization functions, some will lead to better downstream performance than others. Rather than trying to hand-engineer these functions, we may as well let them be learned in an end-to-end fashion with a prediction neural network. Our method can readily be used as an independent module that can be plugged into existing architectures to make them equivariant to arbitrary transformation groups, discrete or continuous. Our approach enjoys similar expressivity advantages to methods like frame averaging by (Puny et al., 2022), but has several added benefits. It is simpler, more efficient, and replaces hand-engineered frames for each group by a systematic end-to-end learning approach.
Our contributions are as follows: • Novel Framework: We introduce a general framework for equivariance to arbitrary groups based on mappings to canonical samples. This framework can be plugged into any existing non-equivariant architecture.
• Theoretical Guarantees: We prove that in some settings, such models are universal approximators of equivariant functions.
• Efficient Implementations: We provide multiple variants of efficient implementations of this framework to specific domains.
• Practical Performance: We perform experiments that show that the proposed method achieves excellent results in the image domain. We also support our hypothesis that learning the canonicalization function is a better strategy than fixing it.

Problem Setting
We are interested in learning functions φ : X → Y with inputs x ∈ X and outputs y ∈ Y belonging to arbitrary finite-dimensional vector spaces. We will consider a set of linear symmetry transformations T ⊂ GL (X ), where GL (X ) is the set of invertible matrices over the vector space X . This is described by a group representation ρ : G → T , where G is an abstract group. Without loss of generality, we can assume that ρ is a group isomorphism. Therefore, the inverse ρ −1 : where the group action ρ on the input and the group action ρ on the output will be clear from the context. In particular, when ρ (g) = I, we say that φ is invariant. We use ρ (G) to denote the image of ρ. We call ρ (G) x = {ρ (g) x | ∀ g ∈ G} the orbit of the element x. It is the set of elements to which x can be mapped to by the group action. The set of orbits, denoted X /G forms a partition of the set X .

General Formulation
The invariance requirement on a function φ amounts to having all the members of a group orbit mapped to the same image by φ. It is thus possible to achieve invariance by appropriately mapping all elements to a canonical sample from their orbit before applying any function. For equivariance, elements can be mapped to a canonical sample and, after a function is applied, transformed back according to their original position in the orbit. This can be formalized by writing the equivariant function φ in canonicalized form as where the function f : X → Y is called the prediction function and the function h : X → ρ (G) is called the canonicalization function. Here h (x) −1 is the inverse of the representation matrix and h (x) = ρ ρ −1 (h (x)) is the counterpart of h (x) on the output. Equivariance in equation 2 is obtained for any prediction function if the canonicalization function is itself G-equivariant 1 , h (ρ (g) x) = ρ (g) h (x) ∀ g, x ∈ G × X .
It may seem like the problem of obtaining an equivariant function has merely been transferred in this formulation. This is however not the case: in equation 2, the equivariance and prediction components are effectively decoupled. The canonicalization function h can therefore be chosen as a simple and inexpressive equivariant function, while the heavy-lifting is done by the prediction function f .
A more general condition can be formulated, such that the decoupling is partial.
and the prediction function f is K-equivariant, then φ defined in equation 2 is G-equivariant.
If K is a normal subgroup such that G J K, this can be realized with a canonicalization function that has image ρ (J), and that is J-equivariant and K-invariant.
The proof follows in Appendix B. This is equivalent to saying that the canonicalization function should output a representation of a coset in G/K in an equivariant way, the applied transformation being chosen arbitrarily within the coset. When K = {e}, only the canonicalization function is constrained, which is the case described above. The other extreme, given by K = G, corresponds to transforming the input in an arbitrary way and constraining the prediction function as is usually done in equivariant architectures. These are respectively the single-view-plus-transformation and the viewpointindependent implementations described in the introduction. Subgroups {e} < K < G offer intermediary options; the lattice of subgroups of G therefore defines a family of models. Since equivariance to a smaller group is less constraining for the prediction function, set inclusion in the subgroup lattice corresponds to increased expressivity.

Universality Result
We can now introduce a more formal result on the expressivity of equivariant functions obtained with canonicalization functions. A parameterized function φ is a universal approximator of G-equivariant functions if for any G-equivariant continuous function ψ, any compact set K ⊆ X and any > 0, there exists a choice of parameters such that ψ (x) − φ (x) < ∀ x ∈ K. We make the additional assumption that the set K is closed under the group action.
Theorem 2 A G-equivariant parameterized function φ written in canonicalized form (equation 2) and satisfying equation 3 with K ≤ G, is a universal approximator of G-equivariant functions if the prediction function is a universal approximator of K-equivariant functions.
The proof follows in Appendix C. The following corollary is especially relevant.
Corollary 3 A G-equivariant parameterized function φ written in canonicalized form with a G-equivariant canonicalization function and a multilayer perceptron (MLP) as a prediction function is a universal approximator of G-equivariant functions.
This result can significantly simplify the design of universal approximators of equivariant functions since a non-universal equivariant architecture for the canonicalization function can be combined with an MLP. solves this issue. However, because the subset of symmetric inputs is of measure zero, using standard equivariance is not expected to be a problem in practice.

Model Design for Canonicalization Functions
In the next section, we elaborate on how suitable canonicalization functions can be obtained in different settings.
The canonicalization function can be chosen as any existing equivariant neural network architecture with the output being a group element; we call this the direct approach. For permutation groups and Lie groups, an equivariant multilayer perceptron (Shawe-Taylor, 1989;Finzi et al., 2021) can be used. We provide examples of implementations in the next section.
We also introduce an alternative method that can prove especially useful when dealing with continuous groups. The canonicalization function can be defined as where s : X → R is a neural network. Such a mapping is G-equivariant provided that s has a unique minimum in each orbit. We call this the optimization approach. Intuitively, the function s represents the unnormalized distance between the input and the canonical sample of the orbit and is therefore minimized when ρ (g) is the transformation that takes to the canonical sample. For continuous groups, a function with unique minimum can be learned with Input Convex Neural Networks (ICNN) (Amos et al., 2017), but satisfying this condition for each orbit is more challenging. The minimization in equation 4 can then be performed using any standard convex optimization method. For practical application, we are typically interested in low-dimensional Lie groups, which makes second-order optimization methods usable. The gradient of the neural network with respect to the parameters can be obtained efficiently (e.g. without having to differentiate through all the optimization steps) with implicit differentiation (Blondel et al., 2022). The optimization implementation presents a close analogy with the mental rotation phenomenon described in the introduction, as humans try to minimize the distance between their representation of an object and the canonical one. As such, it is expected that the optimization process will take more iterations when the input sample is farther away in orbit from the canonical sample. This is consistent with the experimental evidence for mental rotation (Shepard and Metzler, 1971;Carpenter and Eisenberg, 1978).

Euclidean Group
The Euclidean group E (d) is used to capture symmetry to rotations, translations, and reflections. Domains in which this type of symmetry is especially relevant include point cloud modelling, image processing and applications in physics. Below we give the design principles to get an equivariant model for image and point cloud inputs.
Image Input. In order to build networks equivariant to E(2), we need a canonicalization function that outputs an element of E(2) given an image. This can be achieved by using a G-CNN (Cohen and Welling, 2016) and taking an argmax of the output feature map. While G-CNN only considers four rotations, this can be easily extended to arbitrary rotations. This approach can be further simplified if we use a translation equivariant prediction network which is generally the case with most CNN-based architectures. As the translation group T (2) is a normal subgroup of the Euclidean group E(2), Theorem 1 indicates that we only require the canonicalization function to be equivariant to O(2). This means we can eliminate the spatial dimension in the output feature map of the canonicalization function and only need to take an argmax along the point group dimension to identify the correct orientation of the image. However, there are two problems with this approach. First, extending G-CNNs to higher-order rotations will lead to more computation and artifacts due to the finer rotation of filters. Second, we can not backpropagate through the canonicalization function as the argmax operation is not differentiable.
We can avoid the first problem by using a shallower network with a larger filter size. We empirically show why this is a sound choice for canonicalization function in Section 4. Furthermore, to solve the second problem, we can use the straight-through estimator trick (Bengio et al., 2013). Appendix E contains a PyTorch code snippet to perform the canonicalization function of images in a differentiable way using a G-CNN. Similarly, the output of the canonicalization function can be used to invert the feature maps from the prediction network back to their original orientation in a differentiable way.
Since CNNs are universal approximators of T (2)-equivariant functions (Yarotsky, 2022), it follows from Theorem 2 that a CNN augmented with an O(2) equivariant canonicalization function is a universal approximator of E(2)-equivariant functions.
Point Cloud Input. Elements of the Euclidean group can be written as (O, t), where O ∈ R n×n is an orthogonal matrix and t ∈ R n is an arbitrary translation vector. The vector representation on n + 1 dimensional vectors (defined by concatenating a constant 1 to the original vectors) is then defined in the following way This representation is bounded We seek to define an E (d)-equivariant canonicalization function. This can be done by defining it as h (x) = ρ h O (x) , h t (x) , where the functions h O : X → R n×n and h t : X → R n output the point transformation and the translation respectively. By the equivariance condition,  (Zaheer et al., 2017) as a backbone architecture would result in a universal approximator of E (d) and permutation equivariant functions, following Theorem 2 and Theorem 1 of (Segol and Lipman, 2020).

Symmetric Group
The symmetric group S n over a finite set of n elements contains all the permutations of that set. This group captures the inductive bias that input order should not matter. Domains for which S n -equivariance is desirable include object modelling and detection, graph representation learning, and applications in language modelling. S n -equivariant canonicalization functions can be obtained with a direct approach. In particular, the methods introduced by Mena  We propose an alternative of S n -equivariance for the canonicalization function that circumvents these issues. Instead of producing a permutation of the input, the canonicalization function can instead be designed to obtain a clustering of the input, with a function c : X → {1, . . . , k}, where k is the number of possible clusters. This can be made differentiable using the straight-through estimator described for the images above. Then, a corresponding permutation is obtained by sorting elements by cluster index, with elements belonging to the same cluster being sorted arbitrarily. The canonicalization function would therefore respect equation 3 with respect to the subgroup K = S l 1 × · · · × S l k , where l i is the number of elements assigned to cluster i, the product of symmetric groups within clusters. Heuristically, this enforces the inductive bias that set elements should be mapped to high-level classes, across which interactions will be evaluated in an expressive way. This is close in spirit to what is done in the ClusterFormer (Wang et al., 2022) and Reformer (Kitaev et al., 2020) transformer architectures.
The prediction function f would then have to be S l 1 × · · · × S l k equivariant to obtain permutation equivariance. Such a function can readily be obtained with an equivariant multilayer perceptron, as derived by Wang et al. (2020).

Method
Error % ↓ We perform a preliminary empirical analysis of the proposed framework in the image domain. Experiments on other domains will be the subject of more extensive work. We selected the Rotated MNIST dataset (Larochelle et al., 2007), which is used as a benchmark dataset that uses a classification task to test equivariant architectures in prior work (Cohen and Welling, 2016). In all our experiments, we only change the network architecture. We train the models by minimizing the cross entropy loss for 100 epochs using Adam (Kingma and Ba, 2014) with a learning rate of 0.001. We perform early stopping based on the classification performance of the validation dataset with patience of 20 epochs. We first compare our method with different CNN and G-CNN baselines. For CNN (base), we choose an architecture with 7 layers where layer 1 to 3 has 32, 4 to 6 has 64 and layer 7 has 128 channels, respectively. Instead of pooling, we use convolution filters of size 5 × 5 with a stride 2 at layers 4 and 7. The remaining convolutions have filters of size 3 × 3 and stride 1. All the layers are followed by batch-norm and ReLU activation with dropout(p=0.4) only at layers 4 and 7. For the canonicalization function, we choose a shallow G-CNN with three layers. We start with a lifting layer with a filter that is the same size as the input image. This is followed by ReLU and group equivariant layers with 1 × 1 filters.
We consider two variants: an untrained canonicalization function with frozen weights and a canonicalization function learned end-to-end with the prediction function. For a pure G-CNN-based baseline, we provide the value reported by Cohen and Welling (2016) and design a variant which has similar architecture to CNN (base) while matching the number of parameters of our model with canonicalization function. We call this G-CNN (p4 & =params). We consider two more variants of the CNN: one with the same number of parameters as our model and the other where the canonicalization function is done by finding the orientation of the digits using Principal Component Analysis (PCA). We see that using a fixed canonicalization function technique like PCA or canonicalization function with frozen parameters leads to an improvement in performance over the CNN baseline. However, learning the canonicalization function provides a significant performance improvement. Our approach outperforms all the CNN-based baselines and is comparable to G-CNNs. We further notice that adding a canonicalization function for reflections hurts performance. This is, in fact, not surprising as reflection augmentations were not used to build the rotated-MNIST dataset. Adding reflection invariance while the data does not possess it may introduce ambiguity between certain digits like 2 and 5. Next, we perform experiments to understand the role of different components in our model using the group of n discrete rotations (pn).

Ablation Study
First, we vary both the number of layers of the canonicalization network and the number of rotations it is equivariant to. For this, we extend the layers of GCNN to any arbitrary rotations. As we noticed that using a larger filter leads to better performance for higher order rotations we stick to architecture with a lifting layer that has image-sized filters followed by 1 × 1 filters. From Table 4, we notice that adding equivariance to higher order rotation in the canonicalization function leads to significant performance improvement in comparison to adding more layers. Figure 1 shows the canonical orientation resulting from the learnt canonicalization function with a single lifting layer on 110 randomly sampled images of class 7 from the test dataset. This suggests that having a shallow network is enough to learn the correct canonicalization function with a sufficiently high order of discrete rotations. For p64, we see that all the similar-looking samples are aligned in one particular orientation. In contrast, although techniques like PCA or freezing parameters of the canonicalization function finds the correct canonicalization function for simple digits like 1 (see Appendix 3) they struggle to find stable mappings for more complicated digits like 7. Next, we compare the inference time of our model with pure G-CNN-based architectures. For this experiment, we take the CNN architecture of our predictor network and replace the convolutions with group convolutions. As increasing the order of rotation in G-CNN requires more copies of rotated filters in the lifting layer and more parameters in the subsequent group convolution layers, we decreased the number of channels to keep the number of parameters the same as our model. Figure 2 shows that although G-CNN is slightly better for the p4 group increasing the order of discrete rotations improves our model's performance significantly compared to G-CNN. In addition to performance gain, our model's inference speed remains more or less constant while encoding invariance to higher-order rotations due to the shallow canonicalization network. This makes our approach more suitable for building equivariance to bigger groups.

Conclusion
In this work, we propose using learned canonicalization function to obtain equivariant machine learning models. These canonicalization functions can conveniently be plugged into existing architectures, resulting in highly expressive models. We have described general approaches to obtain canonicalization functions and specific implementations for the Euclidean group (for images and point clouds) and the symmetric group.
We performed preliminary experimental studies in the image domain to test our initial hypotheses regarding this framework. First, we show that our approach achieves comparable or better performance than baselines on invariant tasks. Importantly, learning the canonical network is a better approach than using a fixed mapping, either a frozen neural network or a heuristic approach like PCA. Our results also show that the canonicalization function can be realized with a shallow equivariant network without hindering performance. Finally, we show that this approach reduces inference time and is more suitable for bigger groups then G-CNNs.
Areas of future work include an experimental study of this framework on other domains, such as point clouds and graphs. In the image domain, we will also explore building canonicalization functions using steerable networks. The function would output an orientation fibre that transforms by the irreducible representation of the special orthogonal group. Understanding how design choices for canonicalization functions affect downstream performance would also be a potentially fruitful research direction.

Appendix A. Symmetric inputs and relaxed equivariance
We say that an input x ∈ X is symmetric if its stabilizer subgroup G x = {g ∈ G | ρ (g) x = x} is non-trivial. Given any g 1 , g 2 ∈ G, a necessary and sufficient condition for is that g 1 and g 2 are part of the same coset for the stabilizer, e.g. g 1 , g 2 ∈ gG x . This is follows from the well-known relation between orbits and stabilizers. Therefore, symmetric inputs are always fixed by multiple group elements. Symmetric inputs are problematic when using the standard definition of equivariance for the canonicalization function because for g 1 , g 2 ∈ gG x , we have For general group actions, there will not exist a h (x) ∈ ρ (G) such that the last equality is satisfied.
A relaxed version of equivariance can be defined to address this issue.
Definition 4 (Relaxed equivariance) Given group representations ρ : G → GL (X ) and ρ : G → GL (Y), a function h : X → Y satisfies the relaxed equivariance condition if ∀g 1 , x ∈ G × X there exists a g 2 ∈ g 1 G x such that h (ρ (g 1 ) x) = ρ (g 2 ) h (x) This is a generalization of the idea of multiset-equivariance introduced by  to arbitrary group representations. When G x = {e}, standard equivariance is recovered. Canonicalization functions satisfying this condition do not suffer from the aforementioned problem. In addition, this condition is sufficient to obtain relaxed equivariance for canonicalized functions (Eq. 2). This is because, for g 2 ∈ g 1 G x : Using the fact that g 2 −1 g 1 ∈ G x , The fact that we obtain relaxed equivariance for the canonicalized function is a feature rather than a bug. This captures the desideratum that a function should be able to output asymmetric outputs from symmetric inputs, which is not the case for standard equivariant functions.
Since any orbit member can conversely be written as x = ρ (k, j)x for some k ∈ K and j ∈ J, if the canonicalization function is J-equivariant and K-invariant, we have h (x) = ρ (e, j) h (x) h (x) = ρ (e, j) which completes the proof.

Appendix C. Proof of Theorem 2
Proof The proof is inspired by the symmetrization approach of Yarotsky (2022). Let ψ be an arbitrary G-equivariant function, and φ be defined by equation 2. By the equivariance of ψ, we have Since Y is finite-dimensional, we know that linear operators in GL (Y) are bounded. This means that for every representation matrix there exists a positive number c that bounds the induced operator norm, e.g. ∀g ∈ G, ∃c > 0, ρ (g) ≤ c.
We therefore obtain where c > 0. If f is a universal approximator of K-equivarant functions, then it is also a universal approximator of G-equivariant functions. We therefore have Using the fact that K is closed under the group action, we obtain the desired result


The computational complexity of some explainable clustering problems

We study the computational complexity of some explainable clustering problems in the framework proposed by [Dasgupta et al., ICML 2020], where explainability is achieved via axis-aligned decision trees. We consider the $k$-means, $k$-medians, $k$-centers and the spacing cost functions. We prove that the first three are hard to optimize while the latter can be optimized in polynomial time.

Introduction
Machine learning models and algorithms have been used in a number of systems that take decisions that affect our lives. Thus, explainable methods are desirable so that people are able to have a better understanding of their behavior, which allows for comfortable use of these systems or, eventually, the questioning of their applicability [1].
Recently, there has been some effort to devise explainable methods for unsupervised learning tasks, in particular, for clustering [2,3]. We investigate the framework discussed by [2], where an explainable clustering is given by a partition, induced by the leaves of an axis-aligned decision tree, that optimizes some predefined objective function. Figure 1 shows a decision tree that defines a clustering for the Iris dataset. The clustering has three groups, each of them corresponding to a leaf. The explanation of the group associated with the rightmost leaf is Sepal Length >0.4 AND Petal Width < 0.5.
Following [2], a series of papers [4,5,6,7,8] provided algorithms, with provable guarantees, to build decision trees that induce explainable clustering. Several cost functions to guide the clustering construction were investigated as the k-means, k-centers, k-medians and maximum-spacing. Despite this active research, the only work on the field that tackles the computational complexity of building explainable clustering is [9], where it was proved that optimizing the k-means and the k-medians cost functions is NP-Complete. Here, we improve

Problem definition
Let X be a finite set of points in R d . We say that a decision tree is standard if each internal node v is associated with a test (cut), specified by a coordinate i v ∈ [d] and a real value θ v , that partitions the points in X that reach v into two sets: those having the coordinate i v smaller than or equal to θ v and those having it larger than θ v . The leaves of a standard decision tree induce a partition of R d into axis-aligned boxes and, naturally, a partition of X into clusters.
Let k ≥ 2 be an integer. The clustering problems considered here consist of finding a partition of X into k groups, among those that can be induced by a standard decision tree with k leaves, that optimizes a given objective function. For the k-means, k-medians and k-centers cost functions, in addition to the partition, a representative µ(C) ∈ R d for each group C must also be output.
For the k-means problem the objective (cost function) to be minimized is the Sum of the Squared Euclidean Distances (SSED) between each point x ∈ X and the representative of the cluster where x lies. Mathematically, the cost (SSED) of a partition C = (C 1 , . . . , C k ) for X is given by For the k-medians problem the cost of a partition C = (C 1 , . . . , C k ) is given by The k-centers problem is also a minimization problem; its cost function for a partition C = (C 1 , . . . , C k ) is given by Let dist : X × X → R + be a distance function. The meximum-spacing problem consists of finding a partition with at least k groups that has maximum spacing, where the spacing sp(C) of a partition C is defined as sp(C) = min{dist(x, y) : x and y lie in distinct groups of C} In contrast to the other criteria, the spacing is an inter-clustering criterion.
We note that an optimal solution of the unrestricted version of any of these problems, in which the decision tree constraint is not enforced, might be a partition that is hard to explain in terms of the input features. Thus, the motivation for using standard decision trees.
For the sake of simplicity, throughout of this text, by explainable clustering we mean a clustering that is obtained via decision trees.

Our contributions
In Section 2, we first show that the problem of building a partition via decision trees that minimizes the k-means cost function does not admit an (1 + )approximation in polynomial time, for some > 0, unless P = N P . Then, we show that analogous results hold for both the k-median and k-centers cost functions. Our results for both the k-means and k-medians are stronger than the NP-Hardness result established recently by [9] and they formally help to justify the quest for approximation algorithms and/or heuristics for these cost functions.
In Section 3 we propose a polynomial time algorithm that produces an explainable clustering of maximum spacing. As far as we know, this is the first efficient method that produces optimal explainable clustering with respect to some well studied metric.

Related work
Our research is inspired by the recent work of [2], where the problem of building explainable clusterings, via standard decision trees, for both the k-means and the k-medians cost functions are studied. This paper proposes algorithms with provable approximation bounds for building explainable clusterings. In addition, it investigates the price of explainability for these cost functions, which is the unavoidable gap between the cost of the optimal explainable and the optimal unconstrained clustering. Among their results, they showed that the price of explainability for the k-means and k-median are respectively O(k 2 ) and O(k).
Their results were refined/improved by a series of recent papers [4,5,7,8,6]. Currently, the best upper bound for the k-medians is O(log k log log k) [5,7] while for the k-means is O(k log k) [7]. The study of bounds that depend on the dimension d was initiated in [4], where the authors present an O(d log k) upper bound for the k-medians and an O(dk log k) upper bound for the kmeans. These bounds were improved to O(d log 2 d) for the k-medians [7] and O(k 1−2/d poly(d, log k)) [6] for the k-means.
The price of explainability was also investigated for other cost functions. In [4], Laber and Murtinho considered the k-centers and maximum-spacing cost functions. In [5], Makarychev and Shan considered the k-medoids problem (k-median with 2 objective). Finally, in [8], Gamlath et. al addressed p p objectives.
The aforementioned papers, except [4] which also presents experiments, are mainly theoretical. However, there are also a number of papers that propose algorithms (without theoretical guarantees) for building explainable clustering, among them we cite [10,11,3].
The computational complexity of building explainable clustering via decision trees for both the k-means and the k-medians problems is studied in [9]. It is shown that both problems admit polynomial time algorithms when either k or d is constant and they are NP-Complete for arbitrary k and d. In addition, they show that an optimal explainable clustering cannot be found in f (k)·|X | o(k) time for any computable function f (), unless Exponential Time Hypothesis (ETH) fails.
When we turn to standard (non-explainable) clustering, the problems of optimizing the k-means, k-medians and k-centers cost functions are APX-Hard [12,13,14] and all of them admit polynomial time algorithms with constant approximation [15,16,17]. With regards to the spacing cost function, the single-link algorithm, a very popular algorithm to build hierarchical clustering, produces a partition with maximum spacing [18,Chapter 4].
2 Hardness of k-means, k-medians and k-centers cost function

Background
We start by recalling some basic definitions and facts that are useful for studying the hardness of optimization problems (see, e.g., [19, chapter 29]). Given a minimization problem A and a parameter > 0 we define the -Gap-A problem as the problem of deciding for an instance I of A and a parameter k whether: (i) I admits a solution of value at most k; or (ii) every solution of I have value at least (1 + )k. In such a gap decision problem it is tacitly assumed that the instances are either of type (i) or of type (ii). Fact 1. If for a minimization problem A there exists > 0 such that the -Gap-A problem is N P -hard, then no polynomial time (1+ )-approximation algorithm exists for A unless P = N P.
We will use the following definition of a gap-preserving reduction. Definition 1. Let A, B be minimization problems. A gap-preserving reduction from A to B is a polynomial time algorithm that, given an instance x of A and a value k, produces an instance y of B and a value κ such that there exist constants , η > 0 for which If there exists such that the -Gap-A problem is N P -hard and there exists a gap-preserving reduction from A to B then there exists η such that the η-Gap-B problem is N P -hard We will now specialize the above definitions for restricted variants of the problem of finding a minimum vertex cover in a graph and for our clustering problems.
Definition 2. For every > 0, the -Gap-MinVC-B-TF (gap) decision problem is defined as follows: given a triangle-free graph G = (V, E), with bounded degree, and an integer k, decide whether G has a vertex cover of size k or all vertex covers of G have size at least k(1 + ).
The -Gap-MinVC-3B-TF (gap) decision problem has a similar definition, the only differences is that, in addition of being triangle-free, the graphs are required to be 3-bounded, that is, all of its vertexes have degree at most 3.
The NP-Hardness of -Gap-MinVC-B-TF and -Gap-MinVC-3B-TF were established in [13] and [20], respectively. Definition 3. For every η > 0, the η-Gap-Explainable-kmeans (gap) decision problem is defined as follows: given a set of points X , an integer k, and a value κ, decide whether there exists an explainable k-clustering C = (C 1 , . . . C k ) of the points in X such that the k-means cost of C is at most κ or for each explainable k-clustering C of X it holds that the k-means cost of C is at least The η-Gap-Explainable-kmedians and η-Gap-Explainable-kcenters decision problems are analogously defined.
To prove the hardness for the k-means we use a gap-preserving reduction from the -Gap-MinVC-B-TF decision problem. To handle both the k-centers and k-medians, we use the -Gap-MinVC-3B-TF decision problem.
Our reductions have some common ingredients that we explain here. For all of them, given a graph G = (V, E), where V = {1, . . . , n}, we build an instance of the clustering problem under consideration by mapping every edge otherwise. This is exactly the mapping proposed in [13] to establish that the (standard) k-means problem is APX-Hard. We use X G := {v e |e ∈ E} to denote the input of the resulting clustering instance.
Let S = {i 1 , i 2 , . . . , i k } be a cover of size k for G, where each i j is an integer in [n] and i j < i j+1 . We define C S = (E 1 , . . . , E k ) as the k-clustering induced by S on the points in X G , where the group E j includes all points v that simultaneously satisfy: its component i j is 1 and its component i j , for j < j, is 0. Proposition 1. The clustering C S is explainable.
Proof. C S is the clustering induced by a decision tree with k − 1 internal nodes, with exactly one internal node per level. The internal node of level j is associated with cut (i j , 1/2).

Hardness of k-means cost function
We prove that the problem of finding an explainable clustering with minimum k-means cost function is hard to approximate. The reduction employed here is the one used by [13] to show that it is hard to find an (1 + )-approximation for the k-means cost function. The extra ingredient in our proof is the construction of an explainable clustering C S from a vertex cover S that was described in the previous section. Proof. Let G be a triangle-free graph that satisfies one of the following cases: (i) G has a vertex cover of size k or (ii) all vertex covers of G have size > k(1 + ).
First, we consider the case where G has a vertex cover S = {i 1 , i 2 , . . . , i k } of size k. We show that, in this case, the cost of C S = (E 1 , . . . , E k ) is at most |E| − k. Let us consider the mean of the points in E j as the representative of E j , that is, a point that has 1 at coordinate i j and 1/|E j | in the remaining |E j | coordinates with non-zero values. The squared distance of each point in E j to its representative is given by Thus, E j contributes to the total cost with |E j |−1. The cost of the clustering C S is, then, given by k j=1 |E j | − 1 = |E| − k Now, it remains to argue that if the minimum vertex cover for G has size at least (1 + )k then every explainable clustering for the corresponding instance has cost at least |E| − (1 − Ω( ))k. This follows from [13], as in this case every clustering (and, in particular, every explainable one) has cost at least |E| − (1 − Ω( ))k.

Hardness of k-medians cost function
We prove that the problem of finding an explainable clustering with minimum k-medians cost function is hard to approximate. We show a gap preserving reduction from the -Gap-MinVC-3B-TF problem to the η-Gap-Explainablekmedians problem.
The following well-known fact will be useful. The following lemma will be also useful.
Lemma 1. Let G be a 3-bounded triangle free graph and let let C ⊆ X G be a group of points corresponding to p edges of G. We have that: (i) if C is a star then its k-medians cost is p and (ii) if C is not a star then its k-medians cost is at least (4/3)p.
Proof. From the previous fact, the representative of C that yields to the minimum k-medians cost is a point in {0, 1} n , where the coordinate i has value 1 if and only if the number of edges that touch vertex i is larger than p/2. Thus, the cost of a cluster C is given by where d C (i) is the number of edges that touch vertex i in C.
If C is a star centred on vertex j then min{p − d C (j), d C (j)} = 0 and min{p − d C (i), d C (i)} = 1 for the other vertexes i in the star. Thus, If C is not a star then we have some cases: Note that the above case covers the case p ≥ 6 since the maximum degree in G is at most 3.
Case 2) p = 5 and d C (j) = 3 for exactly one j. We have Case 3) p = 5 and d C (j) = d C (j ) = 3 for exactly two values j and j . We have Note that we cannot have 3 vertexes with degree 3 and p = 5.
Case 4) p = 4 and d C (j) = 3 for some j. We must have exactly one j with d C (j) = 3, otherwise we would have more than 4 edges. Thus, Case 5) p = 3 and d C (j) = 2 for some j. We have two possible nonisomorphic graphs. One of them consists of a path with 2 edges and an additional edge while the other is a path with 3 edges. For both cases we have The problem of building an explainable clustering, via decision trees, that minimizes the k−medians cost function does not admit an (1 + )approximation, for some > 0, in polynomial time unless P = N P .
Proof. Let G be a triangle-free graph with maximum degree not larger than 3 that satisfies one of the following cases: (i) G has a vertex cover of size k or (ii) all vertex covers of G have size at least k(1 + ).
First, consider the case where G has a vertex cover S of size k. Since the clustering C S consists of stars, it follows from the previous lemma that its cost is |E|. Now, assume that all vertex covers for G have size at least k(1 + ). Let C be a clustering with k groups for the corresponding k-medians instance.
Let t be the number of groups in C that are stars and let p be the total number of edges in the remaining clusters. Since there is no vertex cover for G of size smaller than k(1 + ) we must have t + p ≥ k(1 + ), otherwise we could obtain a cover for G with size smaller than k(1 + ) by using one vertex per star and one additional vertex for each of the p edges. Since t ≤ k it follows that p ≥ k . Moreover, we must have k ≥ |E|/3 because the degree of every vertex in G is at most 3. Thus, from the previous lemma, the cost of clustering C is at least We have concluded a gap preserving reduction from -Gap-MinVC-3B-TF to η-Gap-Explainable-kmedians.

Hardness of k-centers cost function
In this section we discuss the computational complexity of minimizing the kcenters cost function. We show a gap preserving reduction from the -Gap-MinVC-3B-TF problem to the η-Gap-Explainable-kcenters problem.
Theorem 3. The problem of building an explainable clustering, via decision trees, that minimizes the k−centers cost function does not admit an (1 + )approximation, for some > 0, in polynomial time unless P = N P .
Proof. Let G = (V, E) be a triangle-free graph with maximum degree not larger than 3 that satisfies one of the following cases: (i) G has a vertex cover of size k or (ii) all vertex covers of G have size at least k(1 + ).
First, consider the case where G has a vertex cover S of size k. In this case, the clustering C S = (E 1 , . . . , E k ) consists of stars with at most 3 edges. For the representative of E j , as in the proof of Theorem 1, we use the mean of the points that lie in E j .
Thus, the distance of each point in E i to its representative is the square root of the rightmost term of (1), which is at most 1/3 since G has maximum degree 3. Now, we assume that G does not have a vertex cover with k vertex. Let C be a clustering with k groups for the edges of E. One of the groups, say A, does not have a vertex that touches all the edges in A. Pick the vertex, say v, that touches the largest number of edges in A. Consider an edge e = yz in A that does not touch v. We show that there is another edge in A, say e , that does not have intersection with e. In fact, pick an edge f = vw. If f does not intersect e (w is not an endpoint of e) we set e = f . Otherwise, we assume w.l.o.g. that f intersects e at point y, that is, w = y. We know that vz is not an edge for otherwise we would have a triangle vwz in G. Since v is the vertex that touches the largest number of edges in A then v must touch an edge f = vz , with z = y and z = z. We set e = f . We can argue that the distance of the representative µ(A) of A to either e or e is at least 1. For that, we consider the values of µ(A) at the components of the vertexes that define the edges e and e. Let µ 1 , µ 2 , µ 3 and µ 4 be these values. We have that ) ≥ 2 Thus, either e or e is at distance at least 1 from the representative of A 3 A polynomial time algorithm for the maximumspacing cost function We describe MaxSpacing, a simple greedy algorithm that finds an explainable partition of maximum spacing in polynomial time.
To simplify its description we introduce some notation. For a set of leaves L in a decision tree, we use sp(L) to refer to the spacing of the partition of the points in X induced by the leaves in L. Given a set of leaves L, a leaf ∈ L and an axis-aligned cut γ = (i, θ), we use L γ, to denote the set of leaves obtained when γ is applied to split the points that reach . More precisely, L γ, is obtained from L by removing and adding the two leaves that are created by using γ to split the points that reach .
A pseudo-code for MaxSpacing is presented in Algorithm 1. The algorithm adopts a natural greedy strategy that at each step chooses the cut that yields to the partition of maximum spacing. We note that it runs in polynomial time because in Step 1 we just need to test at most (|X | − 1)d axis-aligned cuts: for each ∈ L and each dimension i ∈ [d] we sort the | | points that reach according to their coordinate i and consider the cuts (i, θ j ), for j = 1, . . . , | |, where θ j is the midpoint between the values of the i-th coordinate of the jth and (j + 1)th points in the sorted list.
Algorithm 1 MaxSpacing(X : set of points; k: integer) Initialize a decision tree with only one leaf and associate it with X L ← { } Repeat k − 1 times: 1. Find a cut γ and a leaf ∈ L that simultaneously satisfy: (i) γ splits the points that reach into 2 non-empty groups (ii) sp(L γ, ) ≥ sp(L γ , ) for every ∈ L and every axis-aligned cut γ that splits the points that reach into two non-empty sets.
2. Split leaf using cut γ In what follows, we show that MaxSpacing produces an explainable partition with maximum possible (optimal) spacing. The following simple fact will be useful. Proof. Let D * be a decision tree that induces a partition with i + 2 groups that has spacing d * i+1 . Let D be a decision tree obtained by removing two leaves that are siblings in D * and turning their parent into a leaf. Let x and y be two closest points among those that reach different leaves in D. Since these points also reach distinct leaves in D * we have that the spacing of the leaves in D is not smaller than that of the leaves in D * . Thus, d * i ≥ sp(Leaves of D) ≥ d * i+1 . Theorem 4. For every 1 ≤ i ≤ k − 1, the partition induced by the leaves of MaxSpacing algorithm by the end of iteration i has the maximum spacing, among the explainable partitions with i + 1 groups for X .
Proof. Let C * i be an optimal explainable partition with i + 1 groups and let d * i be its spacing. Moreover, let L i , with i < k − 1, be the set of leaves by the end of iteration i of MaxSpacing algorithm. By the greedy choice sp(L 1 ) = d * 1 . We assume by induction that the spacing of L i is d * i and show that the spacing of L i+1 is d * i+1 . For a node ν in a decision tree, let P (ν) be the set of points that reach ν. Let ν * be a node in the decision tree D * for C * i+1 that satisfies the following: (i) for each ∈ L i either P ( ) ⊆ P (ν * ) or P ( ) ∩ P (ν * ) = ∅ and (ii) some child of ν * does not satisfy (i). We will use the cut associated with ν * in D * to argue that we can properly split L i .
To prove the existence of a node ν * with such properties, it suffices to show that the root r * of D * satisfies (i) and some leaf * from D * does not satisfy (i) since, in this case, we can set ν * as the last node in the path from r * to * that satisfies (i). Clearly, r * satisfies (i). It remains to argue that some leaf * ∈ D * does not satisfy (i). Since the number of leaves in L i is smaller than the number of leaves in D * , by the pigeonhole principle, there are two leaves, say * 1 and * 2 , in D * that contain points from the same leaf in L i . Thus, neither P ( ) ∩ P ( * 1 ) = ∅ nor P ( ) ⊆ P ( * 1 ). We set * to * 1 and ν * as the last node in the path from r * to * 1 that satisfies (i). Let ν * ch be a child of ν * in D * that does not satisfy (i). Moreover, let γ be the cut associated with ν * and let be a leaf in L i such that P ( ) ⊆ P (ν * ), P ( ) ∩ P (ν * ch ) = ∅ and P ( ) ⊂ P (ν * ch ). We show that the spacing of the set of leaves L i obtained from L i by applying cut γ to is at least d * i+1 . Let 1 and 2 be the two new leaves that are created by applying γ to and let x and y be the two closest points (according to dist) among those that reach different leaves in L i . If x reaches 1 (resp. 2 ) and y reaches 2 (resp. 1 ) then because, due to the application of cut γ on ν * , x and y lies in different groups in C * i+1 . If some of them, say x, does not reach then x and y reach different leaves in L i and, thus, , where the last inequality follows from Fact 4.
We have shown that there exists a leaf in L i and a cut γ such that the application of γ to yields to a partition of spacing d * i+1 . Thus, due to the greedy choice, MaxSpacing obtains a partition of spacing d * i+1 by the end of iteration i + 1.

Conclusions
We have showed that the problems of finding explainable clustering (via decision trees) that optimize the classical k-means, k-medians and k-centers cost functions do not admit polynomial time (1 + )-approximations. These results help to formally justify the quest for heuristics and/or approximation algorithms.
The algorithms recently proposed in the literature for building explainable clustering compare their costs with the costs of optimal unrestricted clustering [2,4,5,6,7,8]. A major open question in this line of research is whether better bounds can be obtained when the comparison is made against the optimal explainable clustering.
For the spacing cost function we provided a simple polynomial time algorithm that computes the explainable partition with maximum spacing. An interesting note is that we have not used the fact that the cuts are axis-aligned in the proof of Theorem 4 and, thus, our result holds for any family of cuts.