Size Matters: Cardinality-Constrained Clustering and Outlier Detection via Conic Optimization

TL;DR

This paper introduces a novel approach to clustering that incorporates outlier detection and cardinality constraints using conic optimization, improving robustness and balance in clustering results.

Contribution

It formulates a joint outlier detection and clustering problem as a MILP with relaxations and deterministic rounding, providing optimal solutions under certain conditions.

Findings

Proposes a MILP formulation for combined outlier detection and clustering.

Develops relaxations and rounding schemes with optimality guarantees.

Addresses cluster balance and outlier sensitivity issues in K-means.

Abstract

Plain vanilla K-means clustering has proven to be successful in practice, yet it suffers from outlier sensitivity and may produce highly unbalanced clusters. To mitigate both shortcomings, we formulate a joint outlier detection and clustering problem, which assigns a prescribed number of datapoints to an auxiliary outlier cluster and performs cardinality-constrained K-means clustering on the residual dataset, treating the cluster cardinalities as a given input. We cast this problem as a mixed-integer linear program (MILP) that admits tractable semidefinite and linear programming relaxations. We propose deterministic rounding schemes that transform the relaxed solutions to feasible solutions for the MILP. We also prove that these solutions are optimal in the MILP if a cluster separation condition holds.

Figures

Tables

Table 1. Table 1: Comparison of Recovery Guarantees for K 𝐾 K -means Clustering Relaxations.

	Awasthi et al.	Iguchi et al.	Ames	Vinayak and Hassibi	This Paper
data generating model	stochastic ball	stochastic ball	stochastic block	stochastic block	none/arbitrary
type of relaxation	SDP + LP	SDP	SDP	SDP	SDP + LP
type of guarantee	stochastic	stochastic	stochastic	stochastic	deterministic
guarantee depends on $N$	yes	yes	yes	yes	no
guarantee depends on $d$	yes	yes	no	no	no
requires balancedness	yes	yes	no	no	yes
proof technique	primal-dual	primal-dual	primal-dual	primal-dual	valid cuts
access to cardinalities	no	no	no	no	yes
outlier detection	no	no	yes	yes	yes

Table 2. Table 2: Overview of the main dataset characteristics.

		$N$	$d$	$K$	$n_k$
ID	Dataset Name	(# datapoints)	(# dimensions)	(# clusters)	(cardinalities)	balanced
1	Iris	150	4	3	50, 50, 50	yes
2	Seeds	210	7	3	70, 70, 70	yes
3	Planning Relax	182	12	2	130, 52	no
4	Connectionist Bench	208	60	2	111, 97	no
5	Urban Land Cover	168	147	9	23, 29, 14, 15, 17, 25, 16, 14, 15	no
6	Parkinsons	195	22	2	48, 147	no
7	Glass Identification	214	9	6	70, 76, 17, 13, 9, 29	no

Table 3. Table 3: Performance of ℛ LP subscript ℛ LP \mathcal{R}_{\text{LP}} , ℛ SDP subscript ℛ SDP \mathcal{R}_{\text{SDP}} , Bennett et al., and Peng and Wei. The “–” signs indicate that the problem instance could not be solved within a time limit of three hours.

	${ℛ}_{\text{LP}}/{ℛ}_{\text{LP}}^{\text{b}}$			${ℛ}_{\text{SDP}}/{ℛ}_{\text{SDP}}^{\text{b}}$			Bennett et al.			$𝒫{𝒲}_1/𝒫{𝒲}_1^b$		$𝒫{𝒲}_2$
ID	UB	LB	time [s]	UB	LB	time [s]	UB	CV [%]	time [s]	LB	time [s]	LB	time [s]
1	81.4	78.8	17	81.4	81.4	584	81.4	0.0	6	81.4	154	15.2	0.02
2	620.7	539.0	46	605.6	605.6	3,823	605.6	0.0	7	604.5	1,320	19.0	0.03
3	325.9	297.0	24	315.7	315.7	2,637	315.8	0.3	9	299.0	510	273.7	0.02
4	312.6	259.1	49	280.6	280.1	3,638	280.6	0.4	6	270.0	1,376	246.2	0.04
5	3.61$e$9	3.17$e$9	2,241	3.54$e$9	3.44$e$9	10,754	3.64$e$9	9.2	13	2.05$e$9	460	1.94$e$8	0.02
6	1.36$e$6	1.36$e$6	22	1.36$e$6	1.36$e$6	2,000	1.36$e$6	15.1	7	1.11$e$6	777	6.31$e$5	0.02
7	469.0	377.2	232	–	–	–	438.2	28.4	13	321.9	1,500	23.8	0.03

Equations

$$\begin{array}[]{l@{\quad}l}\displaystyle\text{minimize}&\textstyle\sum_{k=1}^{K}\sum_{i\in I_{k}}\lVert\bm{\xi}_{i}-\frac{1}{n_{k}}(\sum_{j\in I_{k}}\bm{\xi}_{j})\rVert^{2}\\[8.53581pt] \displaystyle\text{subject to}&\displaystyle(I_{1},\ldots,I_{K})\in\mathfrak{P}(n_{1},\ldots,n_{K}),\end{array}$$

$$\mathfrak{P}(n_{1},\ldots,n_{K})=\Big{\{}(I_{1},\ldots,I_{K})\,:\,|I_{k}|=n_{k}\;\;\forall k,\;\;\cup_{k=1}^{K}I_{k}=\{1,\ldots,N\},\;\;I_{k}\cap I_{\ell}=\emptyset\;\;\forall k\neq\ell\Big{\}}$$

$$\sum _{i=1}^n\Vert {\mathbf{\xi }}_i-\frac{1}{n}({\sum }_{j=1}^n{\mathbf{\xi }}_j){\Vert }^2=\frac{1}{2n}\sum _{i,j=1}^n\Vert {\mathbf{\xi }}_i-{\mathbf{\xi }}_j{\Vert }^2.$$

$$\underset{\sigma \in {\mathfrak{S}}^N}{\text{emminimize}}\frac{1}{2}⟨\mathbf{W},{\mathbf{P}}_{\sigma }\mathbf{D}{\mathbf{P}}_{\sigma }^{\top }⟩,$$

$$\sum _{k=1}^K\sum _{i\in I_k}\Vert {\mathbf{\xi }}_i-\frac{1}{n_k}({\sum }_{j\in I_k}{\mathbf{\xi }}_j){\Vert }^2$$

$$\frac{1}{2}{\sum }_{k=1}^K\frac{1}{n_k}{\sum }_{i,j=1}^N{d}_{ij}{\eta }_{ij}^k$$

$${\pi }_i^k\in \{0,1\},{\eta }_{ij}^k\in {\mathbb{R}}_{+}$$

$${\sum }_{i=1}^N{\pi }_i^k=n_k$$

$${\sum }_{k=1}^K{\pi }_i^k=1$$

$${\eta }_{ij}^k\ge {\pi }_i^k+{\pi }_j^k-1$$

$$\frac{1}{2}{\sum }_{k=1}^K\frac{1}{n_k}{\sum }_{i,j=1}^N{d}_{ij}\max \{{\pi }_i^k+{\pi }_j^k-1,0\}$$

$${\pi }_i^k\in \{0,1\}$$

$${\sum }_{i=1}^N{\pi }_i^k=n_k$$

$${\sum }_{k=1}^K{\pi }_i^k=1$$

$$\sum _{k=1}^K\sum _{i\in I_k}\Vert {\mathbf{\xi }}_i-\frac{1}{n_k}({\sum }_{j\in I_k}{\mathbf{\xi }}_j){\Vert }^2=\frac{1}{2}\sum _{k=1}^K\frac{1}{n_k}\sum _{i,j\in I_k}{d}_{ij}=\frac{1}{2}\sum _{k=1}^K\frac{1}{n_k}\sum _{i,j=1}^N{d}_{ij}\max \{{\pi }_i^k+{\pi }_j^k-1,0\}.$$

$${\mathbf{h}}^{\top }({\mathbf{\xi }}_{i_k}-{\mathbf{\xi }}_{i_{\ell }})<0$$

$$⟺{\mathbf{\xi }}_{i_k}^{\top }{\mathbf{\zeta }}_k+{\mathbf{\xi }}_{i_{\ell }}^{\top }{\mathbf{\zeta }}_{\ell }<{\mathbf{\xi }}_{i_k}^{\top }{\mathbf{\zeta }}_{\ell }+{\mathbf{\xi }}_{i_{\ell }}^{\top }{\mathbf{\zeta }}_k$$

$$⟺\Vert {\mathbf{\xi }}_{i_{\ell }}-{\mathbf{\zeta }}_k{\Vert }^2+\Vert {\mathbf{\xi }}_{i_k}-{\mathbf{\zeta }}_{\ell }{\Vert }^2<\Vert {\mathbf{\xi }}_{i_k}-{\mathbf{\zeta }}_k{\Vert }^2+\Vert {\mathbf{\xi }}_{i_{\ell }}-{\mathbf{\zeta }}_{\ell }{\Vert }^2,$$

$$\sum _{i\in {\tilde{I}}_k}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_k{\Vert }^2+\sum _{i\in {\tilde{I}}_{\ell }}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_{\ell }{\Vert }^2+\sum _{\begin{array}{c}m=1,\dots ,K\\ m\notin \{k,\ell \}\end{array}}\sum _{i\in I_m}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_m{\Vert }^2$$

$$<\sum _{i\in I_k}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_k{\Vert }^2+\sum _{i\in I_{\ell }}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_{\ell }{\Vert }^2+\sum _{\begin{array}{c}m=1,\dots ,K\\ m\notin \{k,\ell \}\end{array}}\sum _{i\in I_m}\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_m{\Vert }^2.$$

$$\frac{1}{2}\sum _{k=1}^K\frac{1}{n_k}\sum _{i,j=1}^N{d}_{ij}{\eta }_{ij}^k=\frac{1}{2}\sum _{k=1}^K\frac{1}{n_k}\sum _{i,j=1}^N{d}_{ij}\max \{{\pi }_i^k+{\pi }_j^k-1,0\}=\frac{1}{2}\sum _{k=1}^K\frac{1}{n_k}\sum _{i,j=1}^N{d}_{ij}{\pi }_i^k{\pi }_j^k,$$

$$\frac{1}{8}{\sum }_{k=1}^K\frac{1}{n_k}{\sum }_{i,j=1}^N{d}_{ij}(1+x_i^k)(1+x_j^k)$$

$$x_i^k\in \{-1,+1\}$$

$${\sum }_{i=1}^N{x}_i^k=2n_k-N$$

$${\sum }_{k=1}^K{x}_i^k=2-K$$

$$\frac{1}{8}⟨\mathbf{D},{\sum }_{k=1}^K\frac{1}{n_k}\left({\mathbf{M}}^k+11^{\top }+{\mathbf{x}}^k{1}^{\top }+1({\mathbf{x}}^k{)}^{\top }\right)⟩$$

$$({\mathbf{x}}^k,{\mathbf{M}}^k)\in {\mathcal{C}}_{\text{SDP}}(n_k)k=1,\dots ,K$$

$${\sum }_{k=1}^K{\mathbf{x}}^k=(2-K)1,$$

$$\mathcal{C}_{\text{SDP}}(n)=\left\{(\bm{x},\mathbf{M})\in\mathbb{R}^{N}\times\mathbb{S}^{N}:~{}\begin{array}[]{l}\bm{1}^{\top}\bm{x}=2n-N,\;\mathbf{M}\bm{1}=(2n-N)\bm{x}\\ \text{diag}(\mathbf{M})=\bm{1},\;\mathbf{M}\succeq\bm{xx}^{\top}\\ \mathbf{M}+\bm{1}\bm{1}^{\top}+\bm{x}\bm{1}^{\top}+\bm{1}\bm{x}^{\top}\geq\bm{0}\\ \mathbf{M}+\bm{1}\bm{1}^{\top}-\bm{x}\bm{1}^{\top}-\bm{1}\bm{x}^{\top}\geq\bm{0}\\ \mathbf{M}-\bm{1}\bm{1}^{\top}+\bm{x}\bm{1}^{\top}-\bm{1}\bm{x}^{\top}\leq\bm{0}\\ \mathbf{M}-\bm{1}\bm{1}^{\top}-\bm{x}\bm{1}^{\top}+\bm{1}\bm{x}^{\top}\leq\bm{0}\end{array}\right\}.$$

$${\mathbf{M}}^{k}1={\mathbf{x}}^k({\mathbf{x}}^k{)}^{\top }1=(2n_k-N){\mathbf{x}}^k,\text{diag}({\mathbf{M}}^k)=1,{\mathbf{M}}^k⪰{\mathbf{x}}^k({\mathbf{x}}^k{)}^{\top }$$

$${\mathbf{M}}^k+{11}^{\top }+{\mathbf{x}}^k{1}^{\top }+1({\mathbf{x}}^k{)}^{\top }=+(1+{\mathbf{x}}^k)(1+{\mathbf{x}}^k{)}^{\top }\ge 0$$

$${\mathbf{M}}^k+{11}^{\top }-{\mathbf{x}}^k{1}^{\top }-1({\mathbf{x}}^k{)}^{\top }=+(1-{\mathbf{x}}^k)(1-{\mathbf{x}}^k{)}^{\top }\ge 0$$

$${\mathbf{M}}^k-{11}^{\top }+{\mathbf{x}}^k{1}^{\top }-1({\mathbf{x}}^k{)}^{\top }=-(1-{\mathbf{x}}^k)(1+{\mathbf{x}}^k{)}^{\top }\le 0$$

$${\mathbf{M}}^k-{11}^{\top }-{\mathbf{x}}^k{1}^{\top }+1({\mathbf{x}}^k{)}^{\top }=-(1+{\mathbf{x}}^k)(1-{\mathbf{x}}^k{)}^{\top }\le 0,$$

$$\text{(a) }{\sum }_{k=1}^K{x}_i^k=2-K\forall i,\text{(b) }{\sum }_{i=1}^N{x}_i^k=2n_k-N\forall k,\text{(c) }{m}_{ij}^k-x_i^k-x_j^k+1\ge 0\forall i,j,k.$$

$$\text{(a’) }{\sum }_{k=1}^K{\pi }_i^k=1\forall i,\text{(b’) }{\sum }_{i=1}^N{\pi }_i^k=n_k\forall k,\text{(c’) }{\eta }_{ij}^k\ge {\pi }_i^k+{\pi }_j^k-1\forall i,j,k.$$

$${\mathbf{\Pi }}^{\prime }\in \underset{\mathbf{\Pi }}{\mathrm{argmax}}\left\{\sum _{i=1}^N\sum _{k=1}^K{\pi }_i^k{x}_i^k:{\pi }_i^k\in \{0,1\},\sum _{i=1}^N{\pi }_i^k=n_k\forall k,\sum _{k=1}^K{\pi }_i^k=1\forall i\right\}.$$

$${\mathbf{\Pi }}^{\star }\in \underset{\mathbf{\Pi }}{\mathrm{argmin}}\left\{\sum _{i=1}^N\sum _{k=1}^K{\pi }_i^k\Vert {\mathbf{\xi }}_i-{\mathbf{\zeta }}_k{\Vert }^2:{\pi }_i^k\in \{0,1\},\sum _{i=1}^N{\pi }_i^k=n_k\forall k,\sum _{k=1}^K{\pi }_i^k=1\forall i\right\}.$$

$$\frac{1}{8n}⟨\mathbf{D},{\mathbf{M}}^1+{11}^{\top }+{\mathbf{x}}^1{1}^{\top }+1({\mathbf{x}}^1{)}^{\top }+(K-1)\left(\mathbf{M}+{11}^{\top }+{\mathbf{x}1}^{\top }+{1\mathbf{x}}^{\top }\right)⟩$$

Full text

\OneAndAHalfSpacedXI\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough\MANUSCRIPTNO

MS-14-00887

\RUNTITLE

Size Matters: Cardinality-Constrained Clustering and Outlier Detection via Conic Optimization \RUNAUTHORRujeerapaiboon, Schindler, Kuhn, Wiesemann \TITLESize Matters: Cardinality-Constrained Clustering and Outlier Detection via Conic Optimization

\ARTICLEAUTHORS\AUTHOR

Napat Rujeerapaiboon \AFFDepartment of Industrial Systems Engineering and Management, National University of Singapore, Singapore, \EMAILisenapa@nus.edu.sg \AUTHORKilian Schindler, Daniel Kuhn \AFFRisk Analytics and Optimization Chair, École Polytechnique Fédérale de Lausanne, Switzerland,

\EMAILkilian.schindler@epfl.ch, \EMAILdaniel.kuhn@epfl.ch \AUTHORWolfram Wiesemann \AFFImperial College Business School, Imperial College London, United Kingdom, \EMAILww@imperial.ac.uk

\ABSTRACT

Plain vanilla $K$-means clustering has proven to be successful in practice, yet it suffers from outlier sensitivity and may produce highly unbalanced clusters. To mitigate both shortcomings, we formulate a joint outlier detection and clustering problem, which assigns a prescribed number of datapoints to an auxiliary outlier cluster and performs cardinality-constrained $K$-means clustering on the residual dataset, treating the cluster cardinalities as a given input. We cast this problem as a mixed-integer linear program (MILP) that admits tractable semidefinite and linear programming relaxations. We propose deterministic rounding schemes that transform the relaxed solutions to feasible solutions for the MILP. We also prove that these solutions are optimal in the MILP if a cluster separation condition holds.

\KEYWORDS

Semidefinite programming, $K$-means clustering, outlier detection, optimality guarantee

1 Introduction

Clustering aims to partition a set of datapoints into a set of clusters so that datapoints in the same cluster are more similar to another than to those in other clusters. Among the myriad of clustering approaches from the literature, $K$-means clustering stands out for its long history dating back to 1957 as well as its impressive performance in various application domains, ranging from market segmentation and recommender systems to image segmentation and feature learning Jain (2010).

This paper studies the cardinality-constrained $K$-means clustering problem, which we define as the task of partitioning $N$ datapoints $\bm{\xi}_{1},\ldots,\bm{\xi}_{N}\in\mathbb{R}^{d}$ into $K$ clusters $I_{1},\ldots,I_{K}$ of prescribed sizes $n_{1},\ldots,n_{K}$, with $n_{1}+\ldots+n_{K}=N$, so as to minimize the sum of squared intra-cluster distances. We can formalize the cardinality-constrained $K$-means clustering problem as follows,

[TABLE]

where

[TABLE]

denotes the ordered partitions of the set $\{1,\ldots,N\}$ into $K$ sets of sizes $n_{1},\ldots,n_{K}$, respectively.

Our motivation for studying problem (1) is threefold. Firstly, it has been shown by Bennett et al. (2000) and Chen et al. (2006) that the algorithms commonly employed for the unconstrained $K$-means clustering problem frequently produce suboptimal solutions where some of the clusters contain very few or even no datapoints. In this context, cardinality constraints can act as a regularizer that avoids local minima of poor quality. Secondly, many application domains require the clusters $I_{1},\ldots,I_{K}$ to be of comparable size. This is the case, among others, in distributed clustering (where different computer clusters should contain similar numbers of network nodes), market segmentation (where each customer segment will subsequently be addressed by a marketing campaign) and document clustering (where topic hierarchies should display a balanced view of the available documents), see Banerjee and Ghosh (2006) and Balcan et al. (2013). Finally, and perhaps most importantly, $K$-means clustering is highly sensitive to outliers. To illustrate this, consider the dataset in Figure 1, which accommodates three clusters as well as three individual outliers. The $K$-means clustering problem erroneously merges two of the three clusters in order to assign the three outliers to the third cluster (top left graph), whereas a clustering that disregards the three outliers would recover the true clusters and result in a significantly lower objective value (bottom left graph). The cardinality-constrained $K$-means clustering problem, where the cardinality of each cluster is set to be one third of all datapoints, shows a similar behavior on this dataset (graphs on the right). We will argue below, however, that the cardinality-constrained $K$-means clustering problem (1) offers an intuitive and mathematically rigorous framework to robustify $K$-means clustering against outliers. A comprehensive and principled treatment of outlier detection methods can be found in the book of Aggarwal (2013).

To our best knowledge, to date only two solution approaches have been proposed for problem (1). Bennett et al. (2000) combine a classical local search heuristic for the unconstrained $K$-means clustering problem due to Lloyd (1982) with the repeated solution of linear assignment problems to solve a variant of problem (1) that imposes lower bounds on the cluster sizes $n_{1},\ldots,n_{K}$. Banerjee and Ghosh (2006) solve the balanced version of problem (1), where $n_{1}=\cdots=n_{K}$, by sampling a subset of the datapoints, performing a clustering on this subset, and subsequently populating the resulting clusters with the remaining datapoints while adhering to the cardinality constraints. Balanced clustering is also considered by Malinen and Fränti (2014) and Costa et al. (2017). Malinen and Fränti (2014) proceed similarly to Bennett et al. (2000) but take explicit advantage of the Hungarian algorithm to speed up the cluster assignment step within the local search heuristic. Costa et al. (2017) propose a variable neighbourhood search heuristic that starts from a random partition of the datapoints into balanced clusters and subsequently searches for better solutions in the neighbourhood obtained by an increasing number of datapoint swaps between two clusters. Although all of these heuristics tend to quickly produce solutions of high quality, they are not known to be polynomial-time algorithms, they do not provide bounds on the suboptimality of the identified solutions, and their performance may be sensitive to the choice of the initial solution. Moreover, neither of these local search schemes accommodates for outliers.

In recent years, several conic optimization schemes have been proposed to alleviate the shortcomings of these local search methods for the unconstrained $K$-means clustering problem Peng and Wei (2007), Awasthi et al. (2015). Peng and Wei (2007) develop two semidefinite programming relaxations of the unconstrained $K$-means clustering problem. Their weaker relaxation admits optimal solutions that can be characterized by means of an eigenvalue decomposition. They further use this eigenvalue decomposition to set up a modified $K$-means clustering problem where the dimensionality of the datapoints is reduced to $K-1$ (provided that their original dimensionality was larger than that). To obtain an upper bound, they solve this $K$-means clustering problem of reduced dimensionality, which can be done either exactly by enumerating Voronoi partitions, as described in Inaba et al. (1994), or by approximation methods such as those in Hasegawa et al. (1993). Using either approach, the runtime grows polynomially in the number of datapoints $N$ but not in the number of desired clusters $K$. Hence, this method is primarily suitable for small $K$. Similar conic approximation schemes have been developed by Elhamifar et al. (2012) and Nellore and Ward (2015) in the context of unconstrained exemplar-based clustering.

Awasthi et al. (2015) and Iguchi et al. (2017) develop probabilistic recovery guarantees for the stronger semidefinite relaxation of Peng and Wei (2007) when the data is generated by a stochastic ball model (i.e., datapoints are drawn randomly from rotation symmetric distributions supported on unit balls). More specifically, they use primal-dual arguments to establish conditions on the cluster separation under which the semidefinite relaxation of Peng and Wei (2007) recovers the underlying clusters with high probability as the number of datapoints $N$ increases. The condition of Awasthi et al. (2015) requires less separation in low dimensions, while the condition of Iguchi et al. (2017) is less restrictive in high dimensions. In addition, Awasthi et al. (2015) consider a linear programming relaxation of the unconstrained $K$-means clustering problem, and they derive similar recovery guarantees for this relaxation as well.

Two more papers study the recovery guarantees of conic relaxations under a stochastic block model (i.e., the dataset is characterized by a similarity matrix where the expected pairwise similarities of points in the same cluster are higher than those of points in different clusters). Ames (2014) considers the densest $K$-disjoint-clique problem whose aim is to split a given complete graph into $K$ subgraphs such as to maximize the sum of the average similarities of the resulting subgraphs. $K$-means clustering can be considered as a specific instance of this broader class of problems. By means of primal-dual arguments, the author derives conditions on the means in the stochastic block model such that his semidefinite relaxation recovers the underlying clusters with high probability as the cardinality of the smallest cluster increases. Vinayak and Hassibi (2016) develop a semidefinite relaxation and regularize it with the trace of the cluster assignment matrix. Using primal-dual arguments they show that, for specific ranges of the regularization parameter, their regularized semidefinite relaxation recovers the true clusters with high probability as the cardinality of the smallest cluster increases. The probabilistic recovery guarantees of Ames (2014) and Vinayak and Hassibi (2016) can also be extended to datasets containing outliers.

In this paper, we propose the first conic optimization scheme for the cardinality-constrained $K$-means clustering problem (1). Our solution approach relies on an exact reformulation of problem (1) as an intractable mixed-integer linear program (MILP) to which we add a set of valid cuts before relaxing the resulting model to a tractable semidefinite program (SDP) or linear program (LP). The set of valid cuts is essential in strengthening these relaxations. Both relaxations provide lower bounds on the optimal value of problem (1), and they both recover the optimal value of (1) whenever a cluster separation condition is met. The latter requires all cluster diameters to be smaller than the distance between any two distinct clusters and, in case of outlier presence, also smaller than the distance between any outlier and any other point. The same condition (in the absence of outliers) was used in Elhamifar et al. (2012) and Awasthi et al. (2015). Our relaxations also give rise to deterministic rounding schemes which produce feasible solutions that are provably optimal in (1) whenever the cluster separation condition holds. Table 1 compares our recovery guarantees to the results available in the literature. We emphasize that our guarantees are deterministic, that they apply to arbitrary data generating models, that they are dimension-independent, and that they hold for both our SDP and LP relaxations. Finally, our algorithms extend to instances of (1) that are contaminated by outliers and whose cluster cardinalities $n_{1},\ldots,n_{K}$ are not known precisely. We summarize the paper’s contributions as follows.

We derive a novel MILP reformulation of problem (1) that only involves $NK$ binary variables, as opposed to the standard MILP reformulation that contains $N^{2}$ binary variables, and whose LP relaxation is at least as tight as the LP relaxation of the standard reformulation. 2. 2.

We develop lower bounds which exploit the cardinality information in problem (1). Our bounds are tight whenever a cluster separation condition is met. Unlike similar results for other classes of clustering problems, our separation condition is deterministic, model-free and dimension-independent. Furthermore, our proof technique does not rely on the primal-dual argument of SDPs and LPs. 3. 3.

We propose deterministic rounding schemes that transform the relaxed solutions to feasible solutions for problem (1). The solutions are optimal in (1) if the separation condition holds. To our best knowledge, we propose the first tractable solution scheme for problem (1) with optimality guarantees. 4. 4.

We illustrate that our lower bounds and rounding schemes extend to instances of problem (1) that are contaminated by outliers and whose cluster cardinalities are not known precisely.

The remainder of the paper is structured as follows. Section 2 analyzes the cardinality-constrained $K$-means clustering problem (1) and derives the MILP reformulation underlying our solution scheme. Sections 3 and 4 propose and analyze our conic rounding approaches for problem (1) in the absence and presence of outliers, respectively. Section 5 presents numerical experiments, and Section 6 gives concluding remarks. Finally, a detailed description of the heuristic proposed by Bennett et al. (2000) for cardinality-constrained $K$-means clustering is provided in the appendix.

Notation:

We denote by $\bm{1}$ the vector of all ones and by $\|\cdot\|$ the Euclidean norm. For symmetric square matrices $\mathbf{A},\mathbf{B}\in\mathbb{S}^{N}$, the relation $\mathbf{A}\succeq\mathbf{B}$ means that $\mathbf{A}-\mathbf{B}$ is positive semidefinite, while $\mathbf{A}\geq\mathbf{B}$ means that $\mathbf{A}-\mathbf{B}$ is elementwise non-negative. The notation $\left<\mathbf{A},\mathbf{B}\right>=\text{Tr}(\mathbf{A}\mathbf{B})$ represents the trace inner product of $\mathbf{A}$ and $\mathbf{B}$. Furthermore, we use $\text{diag}(\mathbf{A})$ to denote a vector in $\mathbb{R}^{N}$ whose entries coincide with those of $\mathbf{A}$’s main diagonal. Finally, for a set of $N$ datapoints $\bm{\xi}_{1},\ldots,\bm{\xi}_{N}$, we use $\mathbf{D}\in\mathbb{S}^{N}$ to denote the matrix of squared pairwise distances $d_{ij}=\|\bm{\xi}_{i}-\bm{\xi}_{j}\|^{2}$.

2 Problem Formulation and Analysis

We first prove that the clustering problem (1) is an instance of a quadratic assignment problem and transform (1) to an MILP with $NK$ binary variables. Then, we discuss the complexity of (1) and show that an optimal clustering always corresponds to some Voronoi partition of $\mathbb{R}^{d}$.

Our first result relies on the following auxiliary lemma, which we state without proof.

Lemma 2.1

For any vectors $\bm{\xi}_{1},\ldots,\bm{\xi}_{n}\in\mathbb{R}^{d}$, we have

[TABLE]

Proof 2.2

ProofSee Zha et al. (2002, p. 1060). $\square$

Using Lemma 2.1, Costa et al. (2017) notice that the $K$-means objective can be stated as a sum of quadratic terms. In the following proposition, we elaborate on this insight and prove that problem (1) is a specific instance of a quadratic assignment problem.

Proposition 2.3** **(Quadratic Assignment Reformulation)

The clustering problem (1) can be cast as the quadratic assignment problem

[TABLE]

where $\mathbf{W}\in\mathbb{S}^{N}$ is a block diagonal matrix with blocks $\frac{1}{n_{k}}\bm{11}^{\top}\in\mathbb{S}^{n_{k}}$, $k=1,\ldots,K$, $\mathfrak{S}^{N}$ is the set of permutations of $\{1,\ldots,N\}$, and $\mathbf{P}_{\sigma}$ is defined through $(\mathbf{P}_{\sigma})_{ij}=1$ if $\sigma(i)=j$; {\color[rgb]{0,0,0}(\mathbf{P}_{\sigma})_{ij}=0} otherwise.

Proof 2.4

ProofWe show that for any feasible solution of (1) there exists a feasible solution of (2) which attains the same objective value and vice versa. To this end, for any partition $(I_{1},\ldots,I_{K})$ feasible in (1), consider any permutation $\sigma\in\mathfrak{S}^{N}$ that satisfies $\sigma(\{1+\sum_{i=1}^{k-1}n_{i},\ldots,\sum_{i=1}^{k}n_{i}\})=I_{k}$ for all $k=1,\ldots,K$, and denote its inverse by $\sigma^{-1}$. This permutation is feasible in (2), and it achieves the same objective value as $(I_{1},\ldots,I_{K})$ in (1) because

[TABLE]

where the first equality is implied by Lemma 2.1, the second equality is a consequence of the definition of $\sigma$, and the third equality follows from the definition of $\mathbf{W}$.

Conversely, for any $\sigma\in\mathfrak{S}^{N}$ feasible in (2), consider any partition $(I_{1},\ldots,I_{K})$ satisfying $I_{k}=\sigma(\{1+\sum_{i=1}^{k-1}n_{i},\ldots,\sum_{i=1}^{k}n_{i}\})$ for all $k=1,\ldots,K$. This partition is feasible in (1), and a similar reasoning as before shows that the partition achieves the same objective value as $\sigma$ in (2). $\square$

Generic quadratic assignment problems with $N$ facilities and $N$ locations can be reformulated as MILPs with $\Omega(N^{2})$ binary variables via the Kaufmann and Broeckx linearization; see, e.g., Burkard (2013, p. 2741). The LP relaxations of these MILPs are, however, known to be weak, and give a trivial lower bound of zero; see, e.g., Zhang et al. (2013, Theorem 4.1). In Proposition 2.5 below we show that the intra-cluster permutation symmetry of the datapoints enables us to give an alternative MILP reformulation containing only $NK\ll\Omega(N^{2})$ binary variables. We also mention that the related, yet different, cardinality-constrained exemplar-based clustering problem can be formulated as an MILP containing $\Omega(N^{2})$ binary variables; see Mulvey and Beck (1984).

Proposition 2.5** **(MILP Reformulation)

The clustering problem (1) is equivalent to the MILP

[TABLE]

The binary variable $\pi_{i}^{k}$ in the MILP $\mathcal{P}$ satisfies $\pi_{i}^{k}=1$ if $i\in I_{k}$; {\color[rgb]{0,0,0}\pi_{i}^{k}=0} otherwise. At optimality, $\eta^{k}_{ij}=\max\{\pi^{k}_{i}+\pi^{k}_{j}-1,0\}$ is equal to 1 if $i,j\in I_{k}$ (i.e., $\pi^{k}_{i}=\pi^{k}_{j}=1$) and 0 otherwise.

Proof 2.6

Proof of Proposition 2.5 At optimality, the decision variables $\eta_{ij}^{k}$ in problem $\mathcal{P}$ take the values $\eta^{k}_{ij}=\max\{\pi_{i}^{k}+\pi_{j}^{k}-1,0\}$. Accordingly, problem $\mathcal{P}$ can equivalently be stated as

[TABLE]

In the following, we show that any feasible solution of (1) gives rise to a feasible solution of $\mathcal{P}$′ with the same objective value and vice versa. To this end, consider first a partition $(I_{1},\ldots,I_{K})$ that is feasible in (1). Choosing $\pi_{i}^{k}=1$ if $i\in I_{k}$ and $\pi_{i}^{k}=0$ otherwise, $k=1,\ldots,K$, is feasible in $\mathcal{P}$′ and attains the same objective value as $(I_{1},\ldots,I_{K})$ in (1) since

[TABLE]

Here, the first equality is implied by Lemma 2.1, and the second equality follows from the construction of $\pi_{i}^{k}$. By the same argument, every $\pi_{i}^{k}$ feasible in $\mathcal{P}$′ gives rise to a partition $(I_{1},\ldots,I_{K})$, $I_{k}=\{i:\pi_{i}^{k}=1\}$ for $k=1,\ldots,K$, that is feasible in $\mathcal{P}$′ and that attains the same objective value. $\square$

Remark 2.7

Note that zero is a (trivial) lower bound on the objective value of the LP relaxation of the MILP $\mathcal{P}$. As a consequence, this LP relaxation is at least as tight as the LP relaxation of the Kaufmann and Broeckx exact MILP formulation of problem (2), which always yields the lower bound of zero. It is also possible to construct instances where the LP relaxation of the MILP $\mathcal{P}$ is strictly tighter.

$K$-means clustering with cardinality constraints is known to be NP-hard as it is a special case of cardinality-constrained $p$-norm clustering, which was shown to be NP-hard (for any $p>1$) by Bertoni et al. (2012). The restriction to the Euclidean norm (i.e., $p=2$), however, allows for a more concise proof, which is given in the following proposition.

Proposition 2.8

$K$-means clustering with cardinality constraints is NP-hard even for $K=2$. Hence, unless P $=$ NP, there is no polynomial time algorithm for solving problem (1).

Proof 2.9

ProofIn analogy to Proposition 2.5, one can show that the unconstrained $K$-means clustering problem can be formulated as a variant of problem $\mathcal{P}$ that omits the first set of assignment constraints, which require that $\sum_{i=1}^{N}\pi^{k}_{i}=n_{k}$ for all k$=1,\ldots,K$, and replaces the (now unconstrained) cardinality $n_{k}$ in the objective function by the size of $I_{k}$, which can be expressed as $\sum_{i=1}^{N}\pi^{k}_{i}$. If $K=2$, we can thus solve the unconstrained $K$-means clustering problem by solving problem $\mathcal{P}$ for all cluster cardinality combinations $(n_{1},n_{2})\in\{(1,N-1),(2,N-2),\ldots,(\lfloor N/2\rfloor,\lceil N/2\rceil)\}$ and selecting the clustering with the lowest objective value. Thus, in this case, if problem $\mathcal{P}$ were polynomial-time solvable, then so would be the unconstrained $K$-means clustering problem. This, however, would contradict Theorem 1 in Aloise et al. (2009), which shows that the unconstrained $K$-means clustering problem is NP-hard even for $K=2$ clusters. $\square$

In the context of balanced clustering, similar hardness results have been established by Pyatkin et al. (2017). Specifically, they prove that the balanced $K$-means clustering problem is NP-complete for $K\geq 2$ and $\frac{N}{K}\geq 3$ (i.e., the shared cardinality of all clusters is greater than or equal to three). In contrast, if $K\geq 2$ and $\frac{N}{K}=2$ (i.e., each cluster should contain two points), balanced $K$-means clustering reduces to a minimum-weight perfect matching problem that can be solved in polynomial-time by different algorithms; see Cook and Rohe (1999, Table I) for a review.

In $K$-means clustering without cardinality constraints, the convex hulls of the optimal clusters do not overlap, and thus each cluster fits within a separate cell of a Voronoi partition of $\mathbb{R}^{d}$; see e.g., Hasegawa et al. (1993, Theorem 2.1). We demonstrate below that this property is preserved in the presence of cardinality constraints.

Theorem 2.10** **(Voronoi Partition)

For every optimal solution to problem (1), there exists a Voronoi partition of $\mathbb{R}^{d}$ such that each cluster is contained in exactly one Voronoi cell.

Proof 2.11

ProofWe show that for every optimal clustering $(I_{1},\ldots,I_{K})$ of (1) and every $k,\ell\in\{1,\ldots,K\}$, $k<\ell$, there exists a hyperplane separating the points in $I_{k}$ from those in $I_{\ell}$. This in turn implies the existence of the desired Voronoi partition. Given a cluster $I_{m}$ for any $m\in\{1,\ldots,K\}$, define its cluster center as $\bm{\zeta}_{m}=\frac{1}{n_{m}}\sum_{i\in I_{m}}\bm{\xi}_{i}$, and let $\bm{h}=\bm{\zeta}_{k}-\bm{\zeta}_{\ell}$ be the vector that connects the cluster centers of $I_{k}$ and $I_{\ell}$. The statement holds if $\bm{h}^{\top}(\bm{\xi}_{i_{k}}-\bm{\xi}_{i_{\ell}})\geq 0$ for all $i_{k}\in I_{k}$ and $i_{\ell}\in I_{\ell}$ as $\bm{h}$ itself determines a separating hyperplane for $I_{k}$ and $I_{\ell}$ in that case. We thus assume that $\bm{h}^{\top}(\bm{\xi}_{i_{k}}-\bm{\xi}_{i_{\ell}})<0$ for some $i_{k}\in I_{k}$ and $i_{\ell}\in I_{\ell}$. However, this contradicts the optimality of the clustering $(I_{1},\ldots,I_{K})$ because

[TABLE]

where the last equivalence follows from multiplying both sides of the second inequality with $2$ and then completing the squares by adding $\bm{\xi}_{i_{k}}^{\top}\bm{\xi}_{i_{k}}+\bm{\zeta}_{k}^{\top}\bm{\zeta}_{k}+\bm{\xi}_{i_{\ell}}^{\top}\bm{\xi}_{i_{\ell}}+\bm{\zeta}_{\ell}^{\top}\bm{\zeta}_{\ell}$ on both sides. Defining $\tilde{I}_{k}=I_{k}\cup\{i_{\ell}\}\setminus\{i_{k}\}$ and $\tilde{I}_{\ell}=I_{\ell}\cup\{i_{k}\}\setminus\{i_{\ell}\}$, the above would imply that

[TABLE]

The left-hand side of the above inequality represents an upper bound on the sum of squared intra-cluster distances attained by the clustering $(I_{1},\ldots,\tilde{I}_{k},\ldots,\tilde{I}_{\ell},\ldots,I_{K})$ since $\bm{\zeta}_{k}$ and $\bm{\zeta}_{\ell}$ may not coincide with the minimizers $\tfrac{1}{n_{k}}\sum_{i\in\tilde{I}_{k}}\bm{\xi}_{i}$ and $\tfrac{1}{n_{\ell}}\sum_{i\in\tilde{I}_{\ell}}\bm{\xi}_{i}$, respectively. Recall that the cluster centers are chosen so as to minimize the sum of the distances from the cluster center to each point in the cluster. We thus conclude that the clustering $(I_{1},\ldots,\tilde{I}_{k},\ldots,\tilde{I}_{\ell},\ldots,I_{K})$ attains a strictly lower objective value than $(I_{1},\ldots,I_{K})$ in problem (1), which is a contradiction. $\square$

3 Cardinality-Constrained Clustering without Outliers

We now relax the intractable MILP $\mathcal{P}$ to tractable conic programs that yield efficiently computable lower and upper bounds on $\mathcal{P}$.

3.1 Convex Relaxations and Rounding Algorithm

We first eliminate the $\eta^{k}_{ij}$ variables from $\mathcal{P}$ by re-expressing the problem’s objective function as

[TABLE]

where the last equality holds because the variables $\pi^{k}_{i}$ are binary. Next, we apply the variable transformation $x^{k}_{i}\leftarrow 2\pi^{k}_{i}-1$, whereby $\mathcal{P}$ simplifies to

[TABLE]

Here, $x^{k}_{i}$ takes the value $+1$ if the $i$-th datapoint is assigned to cluster $k$ and $-1$ otherwise. Note that the constraints in (LABEL:eq:kmean_obj2) are indeed equivalent to the first two constraints in $\mathcal{P}$, respectively. In Theorem 3.1 below we will show that the reformulation (LABEL:eq:kmean_obj2) of the MILP $\mathcal{P}$ admits the SDP relaxation

[TABLE]

where, for any $n\in\mathbb{N}$, the convex set $\mathcal{C}_{\text{SDP}}(n)\subset\mathbb{R}^{N}\times\mathbb{S}^{N}$ is defined as

[TABLE]

Note that $\mathcal{C}_{\text{SDP}}(n)$ is semidefinite representable because Schur’s complement allows us to express the constraint $\mathbf{M}\succeq\bm{xx}^{\top}$ as a linear matrix inequality; see, e.g., Boyd and Vandenberghe (2004). Furthermore, we point out that the last four constraints in $\mathcal{C}_{\text{SDP}}(n)$ are also used in the reformulation-linearization technique for nonconvex programs, as described by Anstreicher (2009).

We can further relax the above SDP to an LP, henceforth denoted by $\mathcal{R}_{\text{LP}}$, where the constraints $(\bm{x}^{k},\mathbf{M}^{k})\in\mathcal{C}_{\text{SDP}}(n_{k})$ are replaced with $(\bm{x}^{k},\mathbf{M}^{k})\in\mathcal{C}_{\text{LP}}(n_{k})$, and where, for any $n\in\mathbb{N}$, the polytope $\mathcal{C}_{\text{LP}}(n)$ is obtained by removing the non-linear constraint $\mathbf{M}\succeq\bm{xx}^{\top}$ from $\mathcal{C}_{\text{SDP}}(n)$.

Theorem 3.1** **(SDP and LP Relaxations)

We have $\min\mathcal{R}_{\text{LP}}$ $\leq\min$ $\mathcal{R}_{\text{SDP}}$ $\leq\min$ $\mathcal{P}$.

Proof 3.2

ProofThe inequality $\min\mathcal{R}_{\text{LP}}\leq\min\mathcal{R}_{\text{SDP}}$ is trivially satisfied because $\mathcal{C}_{\text{SDP}}(n)$ is constructed as a subset of $\mathcal{C}_{\text{LP}}(n)$ for every $n\in\mathbb{N}$. To prove the inequality $\min\mathcal{R}_{\text{SDP}}\leq\min\mathcal{P}$, consider any set of binary vectors $\{\bm{x}^{k}\}_{k=1}^{K}$ feasible in (LABEL:eq:kmean_obj2) and define $\mathbf{M}^{k}=\bm{x}^{k}(\bm{x}^{k})^{\top}$ for $k=1,\ldots,K$. By construction, the objective value of $\{\bm{x}^{k}\}_{k=1}^{K}$ in (LABEL:eq:kmean_obj2) coincides with that of $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ in $\mathcal{R}_{\text{SDP}}$. Moreover, the constraints in (LABEL:eq:kmean_obj2) imply that

[TABLE]

and

[TABLE]

which ensures that $(\bm{x}^{k},\mathbf{M}^{k})\in\mathcal{C}_{\text{SDP}}(n_{k})$ for every $k$. Finally, the constraint $\sum_{k=1}^{K}\bm{x}^{k}=(2-K)\bm{1}$ in $\mathcal{R}_{\text{SDP}}$ coincides with the last constraint in (LABEL:eq:kmean_obj2). Thus, $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ is feasible in $\mathcal{R}_{\text{SDP}}$. The desired inequality now follows because any feasible point in (LABEL:eq:kmean_obj2) corresponds to a feasible point in $\mathcal{R}_{\text{SDP}}$ with the same objective value. Note that the converse implication is generally false. $\square$

Remark 3.3

In the special case when $K=2$, we can half the number of variables in $\mathcal{R}_{\text{SDP}}$ and $\mathcal{R}_{\text{LP}}$ by setting $\bm{x}^{2}=-\bm{x}^{1}$ and $\mathbf{M}^{2}=\mathbf{M}^{1}$ without loss of generality.

It is possible to show that $\mathcal{R}_{\rm{LP}}$ is at least as tight as the naïve LP relaxation $\mathcal{L}$ of the MILP $\mathcal{P}$, where the integrality constraints are simply ignored. One can also construct instances where $\mathcal{R}_{\rm{LP}}$ is strictly tighter than $\mathcal{L}$. We also emphasize that both LP relaxations entail $\mathcal{O}(N^{2}K)$ variables and $\mathcal{O}(N^{2}K)$ constraints.

Proposition 3.4

We have $\min\mathcal{R}_{\rm{LP}}$ $\geq\min\mathcal{L}$.

Proof 3.5

ProofConsider a feasible solution $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ of $\mathcal{R}_{\text{LP}}$. Its feasibility implies that

[TABLE]

Next, set $\pi^{k}_{i}=(x^{k}_{i}+1)/2$ and $\eta^{k}_{ij}=\frac{1}{4}(m^{k}_{ij}+x^{k}_{i}+x^{k}_{j}+1)$ for all $i,j,k$. Then,

[TABLE]

Hence, this solution is feasible in $\mathcal{L}$. A direct calculation also reveals that both solutions attain the same objective value in their respective optimization problems. This confirms that $\mathcal{R}_{\text{LP}}$ is a relaxation that is at least as tight as $\mathcal{L}$. $\square$

Next, we develop a rounding algorithm that recovers a feasible clustering (and thus an upper bound on $\mathcal{P}$) from an optimal solution of the relaxed problem $\mathcal{R}_{\text{SDP}}$ or $\mathcal{R}_{\text{LP}}$; see Algorithm 1.

Recall that the continuous variables $\bm{x}^{k}=(x_{1}^{k},\ldots,x^{k}_{N})^{\top}$ in $\mathcal{R}_{\text{SDP}}$ and $\mathcal{R}_{\text{LP}}$ correspond to the binary variables in (LABEL:eq:kmean_obj2) with identical names. This correspondence motivates us to solve a linear assignment problem in Step 3 of Algorithm 1, which seeks a matrix $\bm{\Pi}\in\{0,1\}^{N\times K}$ with $\pi_{i}^{k}\approx\frac{1}{2}(x_{i}^{k}+1)$ for all $i$ and $k$ subject to the prescribed cardinality constraints. Note that even though this assignment problem constitutes an MILP, it can be solved in polynomial time because its constraint matrix is totally unimodular, implying that its LP relaxation is exact. Alternatively, one may solve the assignment problem using the Hungarian algorithm; see, e.g., Burkard et al. (2009).

Note that Steps 5–7 of Algorithm 1 are reminiscent of a single iteration of Lloyd’s algorithm for cardinality-constrained $K$-means clustering as described by Bennett et al. (2000). Specifically, Step 5 calculates the cluster centers $\bm{\zeta}_{k}$, while Steps 6 and 7 reassign each point to the nearest center while adhering to the cardinality constraints. Algorithm 1 thus follows just one step of Lloyd’s algorithm initialized with an optimizer of $\mathcal{R}_{\text{SDP}}$ or $\mathcal{R}_{\text{LP}}$. This refinement step ensures that the output clustering is compatible with a Voronoi partition of $\mathbb{R}^{d}$, which is desirable in view of Theorem 2.10.

3.2 Tighter Relaxations for Balanced Clustering

The computational burden of solving $\mathcal{R}_{\text{SDP}}$ and $\mathcal{R}_{\text{LP}}$ grows with $K$. We show in this section that if all clusters share the same size $n$ (i.e., $n_{k}=n$ for all $k$), then $\mathcal{R}_{\text{SDP}}$ can be replaced by

[TABLE]

whose size no longer scales with $K$. Similarly, $\mathcal{R}_{\text{LP}}$ simplifies to the LP $\mathcal{R}^{\text{b}}_{\text{LP}}$ obtained from $\mathcal{R}_{\text{SDP}}^{\text{b}}$ by replacing $\mathcal{C}_{\text{SDP}}(n)$ with $\mathcal{C}_{\text{LP}}(n)$. This is a manifestation of how symmetry can be exploited to simplify convex programs, a phenomenon which is studied in a more general setting by Gatermann and Parrilo (2004).

Corollary 3.6** **(Relaxations for Balanced Clustering)

We have $\min\mathcal{R}^{\text{b}}_{\text{LP}}$ $\leq\min$ $\mathcal{R}_{\text{SDP}}^{\text{b}}$ $\leq\min$ $\mathcal{P}$.

Proof 3.7

ProofThe inequality $\min\mathcal{R}^{\text{b}}_{\text{LP}}\leq\min$ $\mathcal{R}_{\text{SDP}}^{\text{b}}$ is trivially satisfied. To prove the inequality $\min$ $\mathcal{R}_{\text{SDP}}^{\text{b}}$ $\leq\min$ $\mathcal{P}$, we first add the symmetry breaking constraint $x_{1}^{1}=1$ to the MILP $\mathcal{P}$. Note that this constraint does not increase the optimal value of $\mathcal{P}$. It just requires that the cluster containing the datapoint $\bm{\xi}_{1}$ should be assigned the number $k=1$. This choice is unrestrictive because all clusters have the same size. By repeating the reasoning that led to Theorem 3.1, the MILP $\mathcal{P}$ can then be relaxed to a variant of the SDP $\mathcal{R}_{\text{SDP}}$ that includes the (linear) symmetry breaking constraint $x_{1}^{1}=1$. Note that the constraints and the objective function of the resulting SDP are invariant under permutations of the cluster indices $k=2,\ldots,K$ because $n_{k}=n$ for all $k$. Note also that the constraints are not invariant under permutations involving $k=1$ due to the symmetry breaking constraint. Next, consider any feasible solution $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ of this SDP, and define

[TABLE]

Moreover, construct a permutation-symmetric solution $\{(\bm{x}_{\rm s}^{k},\mathbf{M}_{\rm s}^{k})\}_{k=1}^{K}$ by setting

[TABLE]

By the convexity and permutation symmetry of the SDP, the symmetrized solution $\{(\bm{x}_{\rm s}^{k},\mathbf{M}_{\rm s}^{k})\}_{k=1}^{K}$ is also feasible in the SDP and attains the same objective value as $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$. Moreover, as the choice of $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ was arbitrary, we may indeed restrict attention to symmetrized solutions with $\bm{x}^{k}=\bm{x}^{\ell}$ and $\mathbf{M}^{k}=\mathbf{M}^{\ell}$ for all $k,\ell\in\{2,\ldots,K\}$ without increasing the objective value of the SDP. Therefore, the simplified SDP relaxation $\mathcal{R}_{\text{SDP}}^{\text{b}}$ provides a lower bound on $\mathcal{P}$. $\square$

If $n_{k}=n$ for all $k$, then the SDP and LP relaxations from Section 3.1 admit an optimal solution where both $\bm{x}^{k}$ and $\mathbf{M}^{k}$ are independent of $k$, in which case Algorithm 1 performs poorly. This motivates the improved relaxations $\mathcal{R}_{\text{SDP}}^{\text{b}}$ and $\mathcal{R}^{\text{b}}_{\text{LP}}$ involving the symmetry breaking constraint $x^{1}_{1}=1$, which ensures that—without loss of generality—the cluster harboring the first datapoint $\bm{\xi}_{1}$ is indexed by $k=1$. As the symmetry between clusters $2,\ldots,K$ persists and because any additional symmetry breaking constraint would be restrictive, the optimal solutions of $\mathcal{R}_{\text{SDP}}^{\text{b}}$ and $\mathcal{R}^{\text{b}}_{\text{LP}}$ only facilitate a reliable recovery of cluster 1. To recover all clusters, however, we can solve $\mathcal{R}_{\text{SDP}}^{\text{b}}$ or $\mathcal{R}^{\text{b}}_{\text{LP}}$ $K-1$ times over the yet unassigned datapoints, see Algorithm 2. The resulting clustering could be improved by appending one iteration of Lloyd’s algorithm (akin to Steps 5–7 in Algorithm 1).

In contrast, the naïve relaxation $\mathcal{L}$ of $\mathcal{P}$ becomes significantly weaker when all cardinalities are equal. To see this, we note that a solution $\pi^{k}_{i}=1/K$ and $\eta^{k}_{ij}=0$ for all $i,j=1,\ldots,N$ and for all $k=1,\ldots,K$ is feasible in $\mathcal{L}$ (i.e., it satisfies all constraints in problem $\mathcal{P}$ except the integrality constraints which are imposed on $\pi^{k}_{i}$) whenever $K\geq 2$. Hence, the optimal objective value of $\mathcal{L}$ is zero. This could be avoided by adding a symmetry breaking constraint $\pi^{1}_{1}=1$ to problem $\mathcal{L}$ to ensure that the cluster containing the first datapoint $\bm{\xi}_{1}$ is indexed by $k=1$. However, the improvement appears to be marginal.

3.3 Comparison to existing SDP Relaxations

We now compare $\mathcal{R}_{\text{SDP}}$ and $\mathcal{R}_{\text{SDP}}^{\rm{b}}$ with existing SDP relaxations from the literature. First, we report the various SDP relaxations proposed by Peng and Wei (2007) and Awasthi et al. (2015). Then, we establish that two of them are equivalent. Finally, we show that $\mathcal{R}_{\text{SDP}}$ and $\mathcal{R}_{\text{SDP}}^{\rm{b}}$ are relaxations that are at least as tight as their corresponding counterparts from the literature. The numerical experiments in Section 5 provide evidence that this relation can also be strict.

Peng and Wei (2007) suggest two different SDP relaxations for the unconstrained $K$-means clustering problem and an SDP relaxation for the balanced $K$-means clustering problem. All of them involve a Gram matrix $\mathbf{W}\in\mathbb{S}^{N}$ with entries $w_{ij}=\bm{\xi}_{i}^{\top}\bm{\xi}_{j}$. Their stronger relaxation for the unconstrained $K$-means clustering problem takes the form

[TABLE]

where $\mathbb{I}$ denotes the identity matrix of dimension $N$. Note that the constraints $\mathbf{Z}\geq\mathbf{0}$ and $\mathbf{Z}\bm{1}=\bm{1}$ ensure that $\mathbf{Z}$ is a stochastic matrix, and hence all of its eigenvalues lie between 0 and 1. Thus, further relaxing the non-negativity constraints leads to the following weaker relaxation,

[TABLE]

Peng and Wei (2007) also demonstrate that $\mathcal{PW}_{2}$ essentially reduces to an eigenvalue problem, which implies that one can solve $\mathcal{PW}_{2}$ in $\mathcal{O}(KN^{2})$ time; see Golub and Loan (1996). Their SDP relaxation for the balanced $K$-means clustering problem is similar to $\mathcal{PW}_{1}$ and takes the form

[TABLE]

Awasthi et al. (2015) suggest another SDP relaxation for the unconstrained $K$-means clustering problem, based on the same matrix of squared pairwise distances $\mathbf{D}$ considered in this paper,

[TABLE]

The following observation asserts that the stronger relaxation $\mathcal{PW}_{1}$ of Peng and Wei (2007) and the relaxation $\mathcal{A}$ of Awasthi et al. (2015) are actually equivalent.

{observation}

The problems $\mathcal{PW}_{1}$ and $\mathcal{A}$ are equivalent.

Proof 3.8

ProofBegin by expressing the objective of $\mathcal{PW}_{1}$ in terms of the pairwise distance matrix $\mathbf{D}$,

[TABLE]

Here, (a) follows from the observation that the $ij$-th element of the matrix $2\mathbf{W}+\mathbf{D}$ can be written as $2\,\bm{\xi}_{i}^{\top}\bm{\xi}_{j}+\|\bm{\xi}_{i}-\bm{\xi}_{j}\|^{2}=\|\bm{\xi}_{i}\|^{2}+\|\bm{\xi}_{j}\|^{2}$, and (b) uses the insights that $\mathbf{Z}\bm{1}=\bm{1}$ and $\bm{1}^{\top}\mathbf{Z}=\bm{1}^{\top}$. Comparing $\mathcal{PW}_{1}$ and $\mathcal{A}$, identity (4) shows that the two relaxations are equivalent because their objective functions are the same (up to a factor two) while they share the same feasible set. $\square$

Next, we establish that $\mathcal{R}_{\text{SDP}}$ is at least as tight a relaxation of the cardinality-constrained $K$-means clustering problem (1) as the stronger relaxation $\mathcal{PW}_{1}$ of Peng and Wei (2007).

Proposition 3.9

We have $\min\mathcal{R}_{\rm{SDP}}\geq\min$ $\mathcal{PW}_{1}$.

Note that, through Observation 3.3, Proposition 3.9 also implies that $\mathcal{R}_{\text{SDP}}$ is at least as tight as the relaxation $\mathcal{A}$ of Awasthi et al. (2015).

Proof 3.10

Proof of Proposition 3.9 To prove that $\mathcal{R}_{\text{SDP}}$ is at least as tight a relaxation as $\mathcal{PW}_{1}$, we will argue that for every feasible solution $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ of $\mathcal{R}_{\text{SDP}}$ one can construct a solution

[TABLE]

which is feasible in $\mathcal{PW}_{1}$ and achieves the same objective value. We first verify the feasibility of the proposed solution $\overline{\mathbf{Z}}$. Note that $\overline{\mathbf{Z}}$ is symmetric by construction. Next, we can directly verify that $\overline{\mathbf{Z}}$ is positive semidefinite since

[TABLE]

where the third implication is due to the definition of $\mathcal{C}_{\text{SDP}}(n_{k})$, which requires that $\mathbf{M}^{k}\succeq\bm{x}^{k}(\bm{x}^{k})^{\top}$. The last statement holds trivially because any quadratic form is non-negative. Next, we can ensure the element-wise non-negativity of $\overline{\mathbf{Z}}$, again through the definition of $\mathcal{C}_{\text{SDP}}(n_{k})$:

[TABLE]

Furthermore, combining the definition of $\mathcal{C}_{\text{SDP}}(n_{k})$ and the constraint $\sum_{k=1}^{K}\bm{x}^{k}=(2-K)\bm{1}$ of $\mathcal{R}_{\text{SDP}}$, we can see that each row of $\overline{\mathbf{Z}}$ indeed sums up to one:

[TABLE]

Finally, the trace of $\overline{\mathbf{Z}}$ is uniquely determined as follows:

[TABLE]

Thus, $\overline{\mathbf{Z}}$ is feasible in $\mathcal{PW}_{1}$, and it remains to prove that it achieves the same objective value as the original solution $\{(\bm{x}^{k},\mathbf{M}^{k})\}_{k=1}^{K}$ in $\mathcal{R}_{\text{SDP}}$. Invoking relation (4), it is easy to see that

[TABLE]

The proof thus concludes. $\square$

Finally, we assert that $\mathcal{R}_{\text{SDP}}^{\rm{b}}$ is at least as tight a relaxation of the balanced $K$-means clustering problem as the corresponding relaxation $\mathcal{PW}_{1}^{\rm{b}}$ of Peng and Wei (2007).

Proposition 3.11

We have $\min\mathcal{R}_{\rm{SDP}}^{\rm{b}}\geq\min$ $\mathcal{PW}_{1}^{\rm{b}}$.

Proof 3.12

Proof To show that $\mathcal{R}_{\text{SDP}}^{\rm{b}}$ is at least as tight a relaxation as $\mathcal{PW}_{1}^{\rm{b}}$, we will again argue that for every feasible solution $\{(\bm{x}^{1},\mathbf{M}^{1}),(\bm{x},\mathbf{M})\}$ of $\mathcal{R}_{\text{SDP}}^{\rm{b}}$ one can construct a solution

[TABLE]

that is feasible in $\mathcal{PW}_{1}^{\rm{b}}$ and achieves the same objective value. Following similar steps as in the proof of Proposition 3.9, one can verify that $\overline{\mathbf{Z}}$ indeed satisfies $\overline{\mathbf{Z}}\succeq\bm{0}$, $\overline{\mathbf{Z}}\geq\bm{0}$, $\mathbf{\overline{Z}}\bm{1}=\bm{1}$ and $\text{Tr}(\mathbf{\overline{Z}})=K$. In order to see that $\overline{\mathbf{Z}}\leq(K/N)\bm{11}^{\top}$, note from the definition of $\mathcal{C}_{\rm{SDP}}(n)$ (where $n=N/K$ denotes the shared cardinality of all clusters) that

[TABLE]

Using this insight and the constraint $\bm{x}^{1}+(K-1)\bm{x}=(2-K)\bm{1}$ of $\mathcal{R}_{\text{SDP}}^{\rm{b}}$, any arbitrary element $\overline{z}_{ij}$ of $\overline{\mathbf{Z}}$ can be bounded above as desired,

[TABLE]

Finally, a direct calculation reveals that the objective of $\mathcal{R}^{\text{b}}_{\text{SDP}}$ evaluated at $\{(\bm{x}^{1},\mathbf{M}^{1}),(\bm{x},\mathbf{M})\}$ coincides with the objective of $\mathcal{PW}_{1}^{\rm{b}}$ evaluated at $\overline{\mathbf{Z}}$, which from (4) is equal to $\frac{1}{2}\langle\mathbf{D},\overline{\mathbf{Z}}\rangle$. Hence, $\mathcal{R}^{\text{b}}_{\text{SDP}}$ is at least as tight a relaxation as $\mathcal{PW}_{1}^{\rm{b}}$, and the proof concludes. $\square$

Note that while Propositions 3.9 and 3.11 demonstrate that our SDP relaxations $\mathcal{R}_{\rm{SDP}}$ and $\mathcal{R}_{\rm{SDP}}^{\rm{b}}$ are at least as tight as their respective counterparts by Peng and Wei (2007), similar tightness results cannot be established for our LP relaxations. Indeed, our numerical experiments based on real-world datasets in Section 5 show that both $\mathcal{R}_{\rm{LP}}$ and $\mathcal{R}_{\rm{LP}}^{\rm{b}}$ can be strictly weaker than $\mathcal{PW}_{1}$ and $\mathcal{PW}_{1}^{\rm{b}}$, respectively. Furthermore, it is possible to construct artificial datasets on which even $\mathcal{PW}_{2}$ outperforms $\mathcal{R}_{\rm{LP}}$ and $\mathcal{R}_{\rm{LP}}^{\rm{b}}$.

3.4 Perfect Recovery Guarantees

We now demonstrate that the relaxations of Section 3.2 are tight and that Algorithm 2 finds the optimal clustering if the clusters are perfectly separated in the sense of the following assumption.

(S) Perfect Separation:

There exists a balanced partition $(J_{1},\ldots,J_{K})$ of $\{1,\ldots,N\}$ where each cluster $k=1,\ldots,K$ has the same cardinality $|J_{k}|=N/K\in\mathbb{N}$, and

[TABLE]

Assumption (S) implies that the dataset admits the natural balanced clustering $(J_{1},\ldots,J_{K})$, and that the largest cluster diameter (i.e., $\max_{1\leq k\leq K}\max_{i,j\in J_{k}}d_{ij}$) is smaller than the smallest distance between any two distinct clusters (i.e., $\min_{1\leq k_{1}<k_{2}\leq K}\min_{i\in J_{k_{1}},j\in J_{k_{2}}}d_{ij}$).

Theorem 3.13

If Assumption (S) holds, then the optimal values of $\mathcal{R}^{\text{b}}_{\text{LP}}$ and $\mathcal{P}$ coincide. Moreover, the clustering $(J_{1},\ldots,J_{K})$ is optimal in $\mathcal{P}$ and is recovered by Algorithm 2.

Put simply, Theorem 3.13 states that for datasets whose hidden classes are balanced and well separated, Algorithm 2 will succeed in recovering this hidden, provably optimal clustering.

Proof 3.14

Proof of Theorem 3.13Throughout the proof we assume without loss of generality that the clustering $(J_{1},\ldots,J_{K})$ from Assumption (S) satisfies $1\in J_{1}$, that is, the cluster containing the datapoint $\bm{\xi}_{1}$ is assigned the number $k=1$. The proof now proceeds in two steps. In the first step, we show that the optimal values of the LP $\mathcal{R}^{\text{b}}_{\text{LP}}$ and the MILP $\mathcal{P}$ are equal and that they both coincide with the sum of squared intra-cluster distances of the clustering $(J_{1},\ldots,J_{K})$, which amounts to

[TABLE]

In the second step we demonstrate that the output $(I_{1},\ldots,I_{K})$ of Algorithm 2 coincides with the optimal clustering $(J_{1},\ldots,J_{K})$ from Assumption (S). As the algorithm uses the same procedure $K$ times to recover the clusters one by one, it is actually sufficient to show that the first iteration of the algorithm correctly identifies the first cluster, that is, it suffices to prove that $I_{1}=J_{1}$.

Step 1:

For any feasible solution $(\bm{x}^{1},\bm{x},\mathbf{M}^{1},\mathbf{M})$ of $\mathcal{R}^{\text{b}}_{\text{LP}}$, we define $\mathbf{H},\mathbf{W}\in\mathbb{S}^{N}$ through

[TABLE]

From the definition of $\mathcal{C}_{\text{LP}}(n)$ it is clear that $\mathbf{H},\mathbf{W}\geq\bm{0}$. Moreover, we also have that

[TABLE]

A similar calculation for $\mathbf{W}$ reveals that $\sum_{i\neq j}w_{ij}=4n(n-1)$. Next, we consider the objective function of $\mathcal{R}^{\text{b}}_{\text{LP}}$, which can be rewritten in terms of $\mathbf{W}$ and $\mathbf{H}$ as

[TABLE]

The sum on the right-hand side can be viewed as a weighted average of the squared distances $d_{ij}$ with non-negative weights $h_{ij}+(K-1)w_{ij}$, where the total weight is given by

[TABLE]

Furthermore each weight $h_{ij}+(K-1)w_{ij}$ is bounded above by 4 because

[TABLE]

where the inequality holds because $\mathbf{M}^{1},\mathbf{M}\leq\bm{11}^{\top}$ (which we know from the proof of Proposition 3.11) and the last equality follows from the constraint $\bm{x}^{1}+(K-1)\bm{x}=(2-K)\bm{1}$ in $\mathcal{R}^{\text{b}}_{\text{LP}}$.

Hence, the sum on the right hand side of (6) assigns each squared distance $d_{ij}$ with $i\neq j$ a weight of at most $4$, while the total weight equals $4Kn(n-1)$. A lower bound on the sum is thus obtained by assigning a weight of 4 to the $Kn(n-1)$ smallest values $d_{ij}$ with $i\neq j$. Thus, we have

[TABLE]

where the last equality follows from Assumption (S). By Lemma 2.1, the right-hand side of (8) represents the objective value of the clustering $(J_{1},\ldots,J_{K})$ in the MILP $\mathcal{P}$. Thus, $\mathcal{R}^{\text{b}}_{\text{LP}}$ provides an upper bound on $\mathcal{P}$. By Corollary 3.6, $\mathcal{R}^{\text{b}}_{\text{LP}}$ also provides a lower bound on $\mathcal{P}$. We may thus conclude that the LP relaxation $\mathcal{R}^{\text{b}}_{\text{LP}}$ is tight and, as a consequence, that the clustering $(J_{1},\ldots,J_{K})$ is indeed optimal in $\mathcal{P}$.

Step 2:

As the inequality in (8) is tight, any optimal solution to $\mathcal{R}^{\text{b}}_{\text{LP}}$ satisfies $h_{ij}+(K-1)w_{ij}=4$ whenever $i\neq j$ and $i,j\in J_{k}$ for some $k=1,\ldots,K$ (i.e., whenever the datapoints $\bm{\xi}_{i}$ and $\bm{\xi}_{j}$ belong to the same cluster). We will use this insight to show that Algorithm 2 outputs $I_{1}=J_{1}$.

For any $i\in J_{1}$, the above reasoning and our convention that $1\in J_{1}$ imply that $h_{1i}+(K-1)w_{1i}=4$. This in turn implies via (7) that $m_{1i}^{1}=m_{1i}=1$ for all $i\in J_{1}$.

From the definition of $\mathcal{C}_{\text{LP}}(n)$, we know that

[TABLE]

This allows us to conclude that

[TABLE]

where the first equality holds because $\mathbf{M}^{1}\bm{1}=(2n-N)\bm{x}^{1}$, which is one of the constraints in $\mathcal{R}^{\text{b}}_{\text{LP}}$, and because of our convention that $x_{1}^{1}=1$. Hence, the above inequality must be satisfied as an equality, which in turn implies that $m_{1i}^{1}=-1$ for all $i\notin J_{1}$.

For any $i\notin J_{1}$, the $1i$-th entry of the matrix inequality $\mathbf{M}^{1}+\bm{1}\bm{1}^{\top}-\bm{x}^{1}\bm{1}^{\top}-\bm{1}(\bm{x}^{1})^{\top}\geq\bm{0}$ from the definition of $\mathcal{C}_{\text{LP}}(n)$ can be expressed as

[TABLE]

where the implication holds because $m^{1}_{1i}=-1$ for $i\notin J_{1}$ and because $x_{1}^{1}=1$ due to the symmetry breaking constraint in $\mathcal{R}^{\text{b}}_{\text{LP}}$. Similarly, for any $i\in J_{1}$, the $ii$-th entry of the matrix inequality $\mathbf{M}^{1}+\bm{1}\bm{1}^{\top}-\bm{x}^{1}\bm{1}^{\top}-\bm{1}(\bm{x}^{1})^{\top}\geq\bm{0}$ can be rewritten as

[TABLE]

where the implication follows from the constraint $\text{diag}(\mathbf{M}^{1})=\bm{1}$ in $\mathcal{R}^{\text{b}}_{\text{LP}}$.

*As $x_{i}^{1}\leq 1$ for all $i\in J_{1}$ and $x_{i}^{1}\leq-1$ for all $i\notin J_{1}$, the equality constraint $\bm{1}^{\top}\bm{x}^{1}=2n-N$ from the definition of $\mathcal{C}_{\text{LP}}(n)$ can only be satisfied if $x_{i}^{1}=1$ for all $i\in J_{1}$ and $x_{i}^{1}=-1$ for all $i\notin J_{1}$. Since Algorithm 2 constructs $I_{1}$ as the index set of the $n$ largest entries of the vector $\bm{x}^{1}$, we conclude that it must output $I_{1}=J_{1}$ and the proof completes. *$\square$

Theorem 3.13 implies via Corollary 3.6 that the optimal values of $\mathcal{R}_{\text{SDP}}^{\text{b}}$ and $\mathcal{P}$ are also equal. Thus, both the LP and the SDP relaxation lead to perfect recovery.

In the related literature, Assumption (S) has previously been used by Elhamifar et al. (2012) to show that the natural clustering can be recovered in the context of unconstrained exemplar-based clustering whenever a regularization parameter is chosen appropriately. In contrast, our formulation does not rely on regularization parameters. Likewise, Theorem 3.13 is reminiscent of Theorem 9 by Awasthi et al. (2015) which formalizes the recovery properties of their LP relaxation for the unconstrained $K$-means clustering problem. Awasthi et al. (2015) assume, however, that the datapoints are drawn independently from a mixture of $K$ isotropic distributions and provide a probabilistic recovery guarantee that improves with $N$ and deteriorates with $d$. In contrast, our recovery guarantee for constrained clustering is deterministic, model-free and dimension-independent. If Assumption (S) holds, simpler algorithms than Algorithm 1 and 2 can be designed to recover the true clusters. For instance, a simple threshold approach (i.e., assigning datapoints to the same cluster whenever the distance between them falls below a given threshold) would be able to recover the true clusters whenever Assumption (S) holds. It seems unlikely, however, that such approaches would perform well in a setting where Assumption (S) is not satisfied. In fact, Awasthi et al. (2015) show that their LP relaxation fails to recover the true clusters with high probability if Assumption (S) is violated. In contrast, the numerical experiments of Section 5 suggest that Algorithms 1 and 2 perform well even if Assumption (S) is violated.

Remark 3.15

To our best knowledge, there is no perfect recovery result for the cardinality-constrained $K$-means clustering algorithm by Bennett et al. (2000), see Appendix, whose performance depends critically on its initialization. To see that it can be trapped in a local optimum, consider the $N=4$ two-dimensional datapoints $\bm{\xi}_{1}=(0,0)$, $\bm{\xi}_{2}=(a,0)$, $\bm{\xi}_{3}=(a,b)$ and $\bm{\xi}_{4}=(0,b)$ with $0<a<b$, and assume that we seek two balanced clusters. If the algorithm is initialized with the clustering $\{\{1,4\},\{2,3\}\}$, then this clustering remains unchanged, and the algorithm terminates and reports a suboptimal solution with relative optimality gap $b^{2}/a^{2}-1$. In contrast, as Assumption (S) holds, Algorithm 2 recovers the optimal clustering $\{\{1,2\},\{3,4\}\}$ by Theorem 3.13.

4 Cardinality-Constrained Clustering with Outliers

If the dataset is corrupted by outliers, then the optimal value of (1) may be high, indicating that the dataset admits no natural clustering. Note that the bounds from Section 3 could still be tight, i.e., it is thinkable that the optimal clustering is far from ‘ideal’ even if it can be found with Algorithm 2. If we gradually remove datapoints that are expensive to assign to any cluster, however, we should eventually discover an ‘ideal’ low-cost clustering. In the extreme case, if we omit all but $K$ datapoints, then the optimal value of (1) drops to zero, and Algorithm 2 detects the optimal clustering due to Theorem 3.13.

We now show that the results of Section 3 (particularly Proposition 2.5 and Theorem 3.1) extend to situations where $n_{0}$ datapoints must be assigned to an auxiliary outlier cluster indexed by $k=0$ ($\sum_{k=0}^{K}n_{k}=N$), and where neither the distances between outliers and retained datapoints nor the distances between different outliers contribute to the objective function. In fact, we could equivalently postulate that each of the $n_{0}$ outliers forms a trivial singleton cluster. The use of cardinality constraints in integrated clustering and outlier detection has previously been considered by Chawla and Gionis (2013) in the context of local search heuristics. Inspired by this work, we henceforth minimize the sum of squared intra-cluster distances of the $N-n_{0}$ non-outlier datapoints. We first prove that the joint outlier detection and cardinality-constrained clustering problem admits an exact MILP reformulation.

Proposition 4.1** **(MILP Reformulation)

The joint outlier detection and cardinality-constrained clustering problem is equivalent to the MILP

[TABLE]

Proof 4.2

ProofThis is an immediate extension of Proposition 2.5 to account for the outlier cluster. $\square$

In analogy to Section 3.1, one can demonstrate that the MILP $\mathcal{P}^{\text{o}}$ admits the SDP relaxation

[TABLE]

Moreover, $\mathcal{R}^{\text{o}}_{\text{SDP}}$ can be further relaxed to an LP, henceforth denoted by $\mathcal{R}_{\text{LP}}^{\text{o}}$, by replacing the semidefinite representable set $\mathcal{C}_{\text{SDP}}(n_{k})$ in $\mathcal{R}^{\text{o}}_{\text{SDP}}$ with the polytope $\mathcal{C}_{\text{LP}}(n_{k})$ for all $k=0,\ldots,K$.

Theorem 4.3** **(SDP and LP Relaxations)

We have $\min\mathcal{R}_{\text{LP}}^{\text{o}}$ $\leq\min$ $\mathcal{R}^{\text{o}}_{\text{SDP}}$ $\leq\min$ $\mathcal{P}^{\text{o}}$.

Proof 4.4

ProofThis result generalizes Theorem 3.1 to account for the additional outlier cluster. As it requires no fundamentally new ideas, the proof is omitted for brevity. $\square$

The relaxations $\mathcal{R}^{\text{o}}_{\text{SDP}}$ and $\mathcal{R}_{\text{LP}}^{\text{o}}$ not only provide a lower bound on $\mathcal{P}^{\text{o}}$, but they also give rise to a rounding algorithm that recovers a feasible clustering and thus an upper bound on $\mathcal{P}^{\text{o}}$; see Algorithm 3. Note that this procedure calls the outlier-unaware Algorithm 1 as a subroutine.

If all normal clusters are equally sized, i.e., $n_{k}=n$ for $k=1,\ldots,K$, then $\mathcal{R}^{\text{o}}_{\text{SDP}}$ can be replaced by

[TABLE]

whose size no longer scales with $K$. Similarly, $\mathcal{R}_{\text{LP}}^{\text{o}}$ simplifies to the LP $\mathcal{R}_{\text{LP}}^{\text{ob}}$ obtained from $\mathcal{R}^{\text{ob}}_{\text{SDP}}$ by replacing $\mathcal{C}_{\text{SDP}}(n)$ and $\mathcal{C}_{\text{SDP}}(n_{0})$ with $\mathcal{C}_{\text{LP}}(n)$ and $\mathcal{C}_{\text{LP}}(n_{0})$, respectively. Note that the cardinality $n_{0}=N-Kn$ may differ from $n$.

Corollary 4.5** **(Relaxations for Balanced Clustering)

We have $\min\mathcal{R}_{\text{LP}}^{\text{ob}}$ $\leq\min$ $\mathcal{R}^{\text{ob}}_{\text{SDP}}$ $\leq\min$ $\mathcal{P}^{\text{o}}$.

Proof 4.6

ProofThis follows from a marginal modification of the argument that led to Corollary 3.6. $\square$

If the normal clusters are required to be balanced, then Algorithm 3 should be modified as follows. First, in Step 2 the relaxations $\mathcal{R}^{\text{ob}}_{\text{SDP}}$ or $\mathcal{R}_{\text{LP}}^{\text{ob}}$ can be solved instead of $\mathcal{R}^{\text{o}}_{\text{SDP}}$ or $\mathcal{R}_{\text{LP}}^{\text{o}}$, respectively. Moreover, in Step 5 Algorithm 2 must be called as a subroutine instead of Algorithm 1.

In the presence of outliers, the perfect recovery result from Theorem 3.13 remains valid if the following perfect separation condition is met, which can be viewed as a generalization of Assumption (S).

(S’) Perfect Separation:

There exists a partition $(J_{0},J_{1},\ldots,J_{K})$ of $\{1,\ldots,N\}$ where each normal cluster $k=1,\ldots,K$ has the same cardinality $|J_{k}|=(N-n_{0})/K\in\mathbb{N}$, while

[TABLE]

Assumption (S’) implies that the dataset admits the natural outlier cluster $J_{0}$ and the natural normal clusters $(J_{1},\ldots,J_{K})$. It also postulates that the diameter of each normal cluster is strictly smaller than (i) the distance between any two distinct normal clusters and (ii) the distance between any outlier and any other datapoint. Under this condition, Algorithm 3 correctly identifies the optimal clustering.

Theorem 4.7

If Assumption (S’) holds, then the optimal values of $\mathcal{R}_{\text{LP}}^{\text{ob}}$ and $\mathcal{P}^{\text{o}}$ coincide. Moreover, the clustering $(J_{0},\ldots,J_{K})$ is optimal in $\mathcal{P}^{\text{o}}$ and is recovered by Algorithm 3.

Proof 4.8

ProofThe proof parallels that of Theorem 3.13 and can be divided into two steps. In the first step we show that the LP relaxation $\mathcal{R}_{\text{LP}}^{\text{ob}}$ for balanced clustering and outlier detection is tight, and in the second step we demonstrate that Algorithm 3 correctly identifies the clusters $(J_{0},\ldots,J_{K})$. As for the second step, it suffices to prove that the algorithm correctly identifies the outlier cluster $J_{0}$. Indeed, once the outliers are removed, the residual dataset satisfies Assumption (S), and Theorem 3.13 guarantees that the normal clusters $(J_{1},\ldots,J_{K})$ are correctly identified with Algorithm 2.

As a preliminary, note that $(\bm{x},\mathbf{M})\in\mathcal{C}_{\text{LP}}(n)$ implies

[TABLE]

where the implications use $\text{diag}(\mathbf{M})=\bm{1}$. Similarly, $(\bm{x}^{0},\mathbf{M}^{0})\in\mathcal{C}_{\text{LP}}(n_{0})$ implies $-\bm{1}\leq\bm{x}^{0}\leq+\bm{1}$.

Step 1:

For any feasible solution $(\bm{x}^{0},\bm{x},\mathbf{M}^{0},\mathbf{M})$ of $\mathcal{R}_{\text{LP}}^{\text{ob}}$, introduce the auxiliary matrix $\mathbf{H}=\mathbf{M}+\bm{11}^{\top}+\bm{1x}^{\top}+\bm{x1}^{\top}$. Recall from the proof of Theorem 3.13 that $\mathbf{H}\geq\bm{0}$ and

[TABLE]

The constraint $K\bm{x}+\bm{x}^{0}=(1-K)\bm{1}$ from $\mathcal{R}_{\text{LP}}^{\text{ob}}$ ensures via the inequality $-\bm{1}\leq\bm{x}^{0}$ that $\bm{x}\leq(\tfrac{2}{K}-1)\bm{1}$. Recalling from the proof of Theorem 3.13 that $\mathbf{M}\leq\bm{11}^{\top}$, we then find

[TABLE]

Similar arguments as in the proof of Theorem 3.13 reveal that the objective function of the joint outlier detection and (balanced) clustering problem $\mathcal{R}_{\text{LP}}^{\text{ob}}$ can be expressed as

[TABLE]

where the equality follows from Assumption (S’). By Lemma 2.1, the right-hand side of (10) represents the objective value of the clustering $(J_{0},\ldots,J_{K})$ in the MILP $\mathcal{P}^{\text{o}}$. Thus, $\mathcal{R}_{\text{LP}}^{\text{ob}}$ provides an upper bound on $\mathcal{P}^{\text{o}}$. By Corollary 4.5, $\mathcal{R}_{\text{LP}}^{\text{ob}}$ also provides a lower bound on $\mathcal{P}^{\text{o}}$. We may thus conclude that the LP relaxation $\mathcal{R}_{\text{LP}}^{\text{ob}}$ is tight and, as a consequence, that the clustering $(J_{0},\ldots,J_{K})$ is indeed optimal in $\mathcal{P}^{\text{o}}$.

Step 2:

As the inequality in (10) is tight, any optimal solution to $\mathcal{R}_{\text{LP}}^{\text{ob}}$ satisfies $h_{ij}=\frac{4}{K}$ whenever $i\neq j$ and $i,j\in J_{k}$ for some $k=1,\ldots,K$ (i.e., whenever $\bm{\xi}_{i}$ and $\bm{\xi}_{j}$ are not outliers and belong to the same cluster). This in turn implies via (9) that $x_{i}=\frac{2}{K}-1$ for all $i\in\cup_{k=1}^{K}J_{k}$. Furthermore, the constraint $\bm{1}^{\top}\bm{x}=2n-N$ from $\mathcal{C}_{\text{LP}}(n)$ implies

[TABLE]

where the inequality holds because $-\bm{1}\leq\bm{x}$. Thus, the above inequality must in fact hold as an equality, which implies that $x_{i}=-1$ for all $i\in J_{0}$. The constraint $K\bm{x}+\bm{x}^{0}=(1-K)\bm{1}$ from $\mathcal{R}_{\text{LP}}^{\text{ob}}$ further implies that $x^{0}_{i}=-1$ for all $i\in\cup_{k=1}^{K}J_{k}$ and $x^{0}_{i}=+1$ for all $i\in J_{0}$.

*Since Algorithm 3 constructs $I_{0}$ as the index set of the $n_{0}=N-Kn$ largest entries of the vector $\bm{x}^{0}$, we conclude that it must output $I_{0}=J_{0}$, and the proof completes. *$\square$

Remark 4.9** **(Unknown Cluster Cardinalities)

The joint outlier detection and cardinality-constrained clustering problem $\mathcal{P}^{\text{o}}$ can also be used when the number of outliers is not precisely known and only an estimate of the relative size (as opposed to the exact cardinality) of the clusters is available. To this end, we solve $\mathcal{P}^{\text{o}}$ for different values of $n_{0}$, respectively assigning the remaining $N-n_{0}$ datapoints to clusters whose relative sizes respect the available estimates. The value $n_{0}^{\star}$ representing the most reasonable number of outliers to be removed from the dataset can then be determined using the elbow method; see, e.g., Gareth et al. (2017, Chapter 10).

As an illustration, consider again the dataset depicted in Figure 1 which showcases the crux of outlier detection in the context of cardinality-constrained clustering. In Section 1, we inadvertently assumed to have the knowledge that the dataset under consideration was contaminated by three outliers. To demonstrate the practical usefulness of our approach, we will now employ the elbow method to determine the number of outliers $n_{0}$ without making any assumptions about the dataset. As elucidated in Remark 4, the ideal value of $n_{0}$ can be determined by solving problem $\mathcal{P}^{\text{o}}$ repeatedly. However, as $\mathcal{P}^{\text{o}}$ constitutes an intractable optimization problem, we solve its convex relaxations $\mathcal{R}_{\text{LP}}^{\text{ob}}$ and $\mathcal{R}^{\text{ob}}_{\text{SDP}}$ instead and plot the resulting objective values in logarithmic scale in Figure 2. It becomes apparent that $n_{0}^{\star}=3$ is most appropriate as it marks the transition from the initially steep decline pattern of the objective value to a substantially flatter decline pattern. Note that $n_{0}$ needs to be a multiple of $K=3$ to allow for balanced clustering.

5 Numerical Experiments

We now investigate the performance of our algorithms on synthetic as well as real-world clustering problems with and without outliers. All LPs and SDPs are solved with CPLEX 12.7.1 and MOSEK 8.0, respectively, using the YALMIP interface on a 3.40GHz i7 computer with 16GB RAM.

5.1 Cardinality-Constrained $K$-Means Clustering (Real-World Data)

We compare the performance of our algorithms from Section 3 with the algorithm of Bennett et al. (2000), see Appendix, and with the two SDP relaxations proposed by Peng and Wei (2007) on the classification datasets of the UCI Machine Learning Repository (http://archive.ics.uci.edu/ml/) with 150–300 datapoints, up to 200 continuous attributes and no missing values. Table 2 reports the main characteristics of these datasets. In our experiments, we set the cluster cardinalities to the numbers of true class occurrences in each dataset. It should be emphasized that, in contrast to the other methods, and with exception of the two balanced datasets, the SDP relaxations of Peng and Wei (2007) do not have access to the cluster cardinalities. They should thus be seen as a baseline for the performance of the other methods. Furthermore, we remark that all datasets severely violate Assumption (S). Indeed, the ratios of largest cluster diameter to smallest distance between clusters (when the clusters are determined by the true labels) vary from 7 to 149, while they should be smaller than one in order to satisfy Assumption (S). Also, only two datasets actually entail balanced clusters.

Table 3 reports the lower bounds provided by $\mathcal{R}_{\text{LP}}/\mathcal{R}^{\text{b}}_{\text{LP}}$ and $\mathcal{R}_{\text{SDP}}/\mathcal{R}_{\text{SDP}}^{\text{b}}$ (LB), the upper bounds from Algorithms 1 and 2 (UB), the objective value of the best of 10 runs of the algorithm of Bennett et al. (UB), randomly initialized by the cluster centers produced by the K-means++ algorithm of Arthur and Vassilvitskii (2007), the coefficient of variation across these 10 runs (CV), the respective lower bounds (LB) obtained from the SDP relaxations $\mathcal{PW}_{1}$/$\mathcal{PW}_{1}^{\rm{b}}$ and $\mathcal{PW}_{2}$ of Peng and Wei (2007), and the solution times for each of these methods. The latter was limited to a maximum of three hours, and in one case (namely, “Glass Identification”), $\mathcal{R}_{\text{SDP}}$ did not terminate within this limit. The “–” signs in Table 3 indicate this occurrence.

The obtained lower bounds of $\mathcal{R}_{\text{SDP}}$/$\mathcal{R}_{\rm{SDP}}^{\rm{b}}$ allow us to certify that the algorithm of Bennett et al. (2000) provides nearly optimal solutions in almost all instances. Also, both Algorithms 1 and 2 are competitive in terms of solution quality with the algorithm of Bennett et al. (2000) while providing rigorous error bounds. Moreover, as expected in view of Propositions 3.9 and 3.11, for all datasets $\mathcal{R}_{\text{SDP}}$/$\mathcal{R}_{\rm{SDP}}^{\rm{b}}$ yield better lower bounds than the SDP relaxations $\mathcal{PW}_{1}$/$\mathcal{PW}_{1}^{\rm{b}}$ and $\mathcal{PW}_{2}$ of Peng and Wei (2007). The lower bounds obtained from $\mathcal{R}_{\text{LP}}$/$\mathcal{R}_{\rm{LP}}^{\rm{b}}$ are competitive with those provided by the relaxations $\mathcal{PW}_{1}$/$\mathcal{PW}_{1}^{\rm{b}}$, and they are always better than the lower bounds provided by their relaxation $\mathcal{PW}_{2}$. It should be mentioned, however, that one can construct instances where the situation is reversed, i.e., both $\mathcal{PW}_{1}$/$\mathcal{PW}_{1}^{\rm{b}}$ and $\mathcal{PW}_{2}$ are tighter than $\mathcal{R}_{\text{LP}}$/$\mathcal{R}_{\rm{LP}}^{\rm{b}}$. Peng and Wei (2007) also suggest a procedure to compute a feasible clustering (and thus upper bounds) for the unconstrained $K$-means clustering problem. However, this procedure relies on an enumeration of all possible Voronoi partitions, which is impractical for $K\geq 3$; see Inaba et al. (1994). Furthermore, it is not clear how to impose cardinality constraints in this setting.

Specifically, in the context of the two balanced datasets (i.e., “Iris” and “Seeds”), we can enrich the preceding comparison with the heuristics proposed by Costa et al. (2017) and Malinen and Fränti (2014). As for the variable neighbourhood search method of Costa et al. (2017), we were provided with the executables of the C++ implementation used in that paper. For the “Iris” dataset, the best objective value out of 10 independent runs of this method was 81.4 (which is provably optimal thanks to the lower bounds provided by $\mathcal{R}_{\text{SDP}}$ and $\mathcal{PW}_{1}^{\rm{b}}$) and the time to execute all runs was 0.12 seconds. For the “Seeds” dataset, the best objective value out of 10 independent runs was 605.6 (again, provably optimal in view of the lower bound provided by $\mathcal{R}_{\text{SDP}}$) and the overall runtime was 0.53 seconds. The algorithm of Malinen and Fränti (2014) follows the same steps as the one of Bennett et al. (2000) with the improvement that the cluster assignment step is solved by the Hungarian algorithm, which provides better runtime guarantees and typically solves faster than interior-point methods for LPs. For this reason, the upper bounds of Malinen and Fränti (2014) for the “Iris” and “Seeds” dataset coincide with those of Bennett et al. (2000) while their algorithm can be expected to terminate faster. A direct comparison of the time complexity of these two methods can be found in Malinen and Fränti (2014).

5.2 Cardinality-Constrained $K$-Means Clustering (Synthetic Data)

We now randomly generate partitions of $10$, $20$ and $70$ datapoints in $\mathbb{R}^{2}$ that are drawn from uniform distributions over $K=3$ unit balls centered at $\bm{\zeta}_{1},\bm{\zeta}_{2}$ and $\bm{\zeta}_{3}$, respectively, such that $\|\bm{\zeta}_{1}-\bm{\zeta}_{2}\|=\|\bm{\zeta}_{1}-\bm{\zeta}_{3}\|=\|\bm{\zeta}_{2}-\bm{\zeta}_{3}\|=\delta$. Theorem 3.13 shows that $\mathcal{R}^{\text{b}}_{\text{LP}}$ is tight and that Algorithm 2 can recover the true clusters whenever $n_{1}=n_{2}=n_{3}$ and $\delta\geq 4$. Figure 3 demonstrates that in practice, perfect recovery is often achieved by Algorithm 1 even if $\delta\ll 4$ and $n_{1}\neq n_{2}\neq n_{3}$. We also note that $\mathcal{R}_{\text{SDP}}$ outperforms $\mathcal{R}_{\text{LP}}$ when $\delta$ is small, and that the algorithm of Bennett et al. frequently fails to determine the optimal solution even if it is run 10 times. In line with the results from the real-world datasets, $\mathcal{R}_{\text{SDP}}$ and $\mathcal{R}_{\text{LP}}$ are tighter than the stronger SDP relaxation of Peng and Wei (2007). Furthermore, it can be shown that in this setting the weaker relaxation of Peng and Wei (2007) always yields the trivial lower bound of zero. The average runtimes are 7s ($\mathcal{R}_{\text{LP}}$), 106s ($\mathcal{R}_{\text{SDP}}$), 11s (Bennett et al.) and 15.6s (Peng and Wei).

5.3 Outlier Detection

We use $\mathcal{R}_{\text{LP}}^{\text{o}}$ and Algorithm 3 to classify the Breast Cancer Wisconsin (Diagnostic) dataset. The dataset has $d=30$ numerical features, which we standardize using a Z-score transformation, and it contains 357 benign and 212 malignant cases of breast cancer. We interpret the malignant cases as outliers and thus set $K=1$. Figure 4 reports the prediction accuracy as well as the false positives (benign cancers classified as malignant) and false negatives (malignant cancers classified as benign) as we increase the number of outliers $n_{0}$ from [math] to $400$. The figure shows that while setting $n_{0}\approx 212$, the true number of malignant cancers, maximizes the prediction accuracy, any choice $n_{0}\in[156,280]$ leads to a competitive prediction accuracy above 80%. Thus, even rough estimates of the number of malignant cancer datapoints can lead to cancer predictors of decent quality. The average runtime is 286s, and the optimality gap is consistently below 3.23% for all values of $n_{0}$.

6 Conclusion

Clustering is a hard combinatorial optimization problem. For decades, it has almost exclusively been addressed by heuristic approaches. Many of these heuristics have proven to be very successful in practice as they often provide solutions of high, or at least satisfactory, quality within attractive runtimes. The common drawback of these methods is that there is typically no way of certifying the optimality of the provided solutions nor to give guaranteed bounds on their suboptimality.

Maybe precisely because of this shortcoming, more recently, convex optimization approaches have been proposed for solving relaxed versions of the clustering problem. These conic programs are polynomial-time solvable and offer bounds on the suboptimality of a given solution. Furthermore, the solutions of these conic relaxations can be “rounded” to obtain actually feasible solutions to the original clustering problem, which results in a new class of heuristic methods.

The results presented in this paper follow precisely this recent paradigm. Combined, conic relaxations and (rounding) heuristics offer solutions to the clustering problem together with a-posteriori guarantees on their optimality. Naturally, one would also wish for attractive a-priori guarantees on the performance of these combined methods. The conditions required to derive such a-priori guarantees are still quite restrictive, but the strong performance of these methods on practical instances makes us confident that this is a promising avenue for future research.

Acknowledgements:

This research was supported by the Singapore Ministry of Education Academic Research Fund Tier 1 R-266-000-131-133, the SNSF grant BSCGI0_157733 and the EPSRC grant EP/N020030/1. We are grateful to Leandro Costa and Daniel Aloise for sharing with us the executables of the variable neighbourhood search heuristic proposed in Costa et al. (2017).

Appendix: Algorithm of Bennett et al. (2000)

The algorithm of Bennett et al. (2000) is designed for a variant of problem (1), where only lower bounds on the clusters’ cardinalities are imposed. This algorithm has a natural extension to our cardinality-constrained clustering problem (1) as follows.

Algorithm 4 adapts a classical local search heuristic for the unconstrained $K$-means clustering problem due to Lloyd (1982) to problem (1). At initialization, it generates random cluster centers $\bm{\zeta}_{k}$, $k=1,\ldots,K$. Each subsequent iteration of the algorithm consists of two steps. The first step assigns every datapoint $\bm{\xi}_{i}$ to the nearest cluster center while adhering to the prescribed cluster cardinalities, whereas the second step replaces each center $\bm{\zeta}_{k}$ with the mean of the datapoints that have been assigned to cluster $k$. The algorithm terminates when the cluster centers $\bm{\zeta}_{1},\ldots,\bm{\zeta}_{K}$ no longer change.