Improved Conic Reformulations for K-means Clustering

TL;DR

This paper introduces a new conic reformulation of K-means clustering that leads to a tighter SDP relaxation and a novel approximation algorithm, improving cluster identification and outperforming existing methods.

Contribution

It presents a polynomial-sized conic reformulation of K-means, a tighter SDP relaxation, and a new approximation algorithm with empirical performance gains.

Findings

Tighter SDP relaxation for K-means clustering.

New approximation algorithm outperforms existing methods.

Reformulation enables better cluster identification.

Abstract

In this paper, we show that the popular K-means clustering problem can equivalently be reformulated as a conic program of polynomial size. The arising convex optimization problem is NP-hard, but amenable to a tractable semidefinite programming (SDP) relaxation that is tighter than the current SDP relaxation schemes in the literature. In contrast to the existing schemes, our proposed SDP formulation gives rise to solutions that can be leveraged to identify the clusters. We devise a new approximation algorithm for K-means clustering that utilizes the improved formulation and empirically illustrate its superiority over the state-of-the-art solution schemes.

Tables

Table 1. Table 1: Improvement of the true K 𝐾 K -means objective value of the cluster assignment generated from the Algorithm  1 relative to the ones generated from the algorithm in [ 21 ] (left) and the Lloyd’s Algorithm (right).

		Statistic
		Mean		5th Percentile		95th Percentile
$D$	2	47.4%	26.6%	-4.4%	17.6%	186.7%	36.5%
	3	21.3%	18.3%	-2.3%	10.9%	168.9%	25.5%
	4	5.7%	14.5%	-1.5%	9.5%	10.8%	20.8%
	5	9.5%	11.1%	-2.1%	7.3%	125.8%	14.5%
	6	4.8%	10.9%	-0.7%	7.5%	8.4%	13.8%

Equations

$$\begin{array}[]{cll}\min&\displaystyle\sum_{i=1}^{K}\sum_{n\in\mathcal{P}_{i}}\|\bm{x}_{n}-\bm{c}_{i}\|^{2}\\ \textnormal{s.t.}&\displaystyle\mathcal{P}_{i}\subseteq\{1,\dots,N\},\;\bm{c}_{i}\in\mathbb{R}^{D}&\forall i\in\{1,\dots,K\}\\ &\displaystyle\bm{c}_{i}=\frac{1}{|\mathcal{P}_{i}|}\sum_{n\in\mathcal{P}_{i}}\bm{x}_{n}\\ &\mathcal{P}_{1}\cup\dots\cup\mathcal{P}_{K}=\{1,\dots,N\},\;\;\mathcal{P}_{i}\cap\mathcal{P}_{j}=\emptyset&\forall i,j\in\{1,\dots,K\}:i\neq j.\end{array}$$

$$\begin{array}[]{cll}\min&\bm{p}^{\top}\bm{C}_{0}\bm{p}+2\bm{c}_{0}^{\top}\bm{p}\\ \textnormal{s.t.}&\bm{p}\in\mathcal{K}\\ &\bm{A}\bm{p}=\bm{b}\\ &\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}=\phi_{j}&\forall j\in[J]\end{array}$$

$$\mathcal{F}=\left\{\mathbf{p}\in \mathcal{K}:\mathbf{A}\mathbf{p}=\mathbf{b},{\mathbf{p}}^{\top }{\mathbf{C}}_j\mathbf{p}+2{\mathbf{c}}_j^{\top }\mathbf{p}={\varphi }_j\forall j\in [J]\right\}$$

$$\mathcal{Q}=\left\{\left[\begin{array}{c}\mathbf{p}\\ 1\end{array}\right]{\left[\begin{array}{c}\mathbf{p}\\ 1\end{array}\right]}^{\top }:\mathbf{p}\in \mathcal{F}\right\}\text{ and }{\mathcal{Q}}^{\infty }=\left\{\left[\begin{array}{c}\mathbf{d}\\ 0\end{array}\right]{\left[\begin{array}{c}\mathbf{d}\\ 0\end{array}\right]}^{\top }:\mathbf{d}\in {\mathcal{F}}^{\infty }\right\}.$$

$$\begin{array}[]{cll}\min&\textup{tr}(\bm{C}_{0}\bm{Q})+2\bm{c}_{0}^{\top}\bm{p}\\ \textnormal{s.t.}&\begin{bmatrix}\bm{Q}&\bm{p}\\ \bm{p}^{\top}&1\end{bmatrix}\in\textup{cl}\textnormal{conv}\left(\mathcal{Q}\right).\end{array}$$

$${\mathcal{F}}_m=\left\{\mathbf{p}\in \mathcal{K}:\mathbf{A}\mathbf{p}=\mathbf{b},{\mathbf{p}}^{\top }{\mathbf{C}}_j\mathbf{p}+2{\mathbf{c}}_j^{\top }\mathbf{p}={\varphi }_j\forall j\in {\mathcal{T}}_m\right\}\forall m\in [M]\cup \{0\},$$

$${\varphi }_j=\underset{\mathbf{p}\in {\mathcal{F}}_{m-1}}{\min }{\mathbf{p}}^{\top }{\mathbf{C}}_j\mathbf{p}+2{\mathbf{c}}_j^{\top }\mathbf{p}\text{or}{\varphi }_j=\underset{\mathbf{p}\in {\mathcal{F}}_{m-1}}{\max }{\mathbf{p}}^{\top }{\mathbf{C}}_j\mathbf{p}+2{\mathbf{c}}_j^{\top }\mathbf{p}\forall j\in {\mathcal{T}}_m\setminus {\mathcal{T}}_{m-1},$$

$${\mathbf{d}}^{\top }{\mathbf{C}}_j\mathbf{d}+2{\mathbf{d}}^{\top }({\mathbf{C}}_j\overline{\mathbf{p}}+{\mathbf{c}}_j)=0\forall \mathbf{d}\in {\mathcal{F}}^{\infty }\forall j\in [J].$$

$$\mathcal{R}=\left\{\begin{bmatrix}\bm{Q}&\bm{p}\\ \bm{p}^{\top}&1\end{bmatrix}\in\mathcal{C}(\mathcal{K}\times\mathbb{R}_{+}):\begin{array}[]{l}\displaystyle\bm{A}\bm{p}=\bm{b},\;\operatorname*{diag}(\bm{A}\bm{Q}\bm{A}^{\top})=\bm{b}\circ\bm{b}\\ \displaystyle\textup{tr}(\bm{C}_{j}\bm{Q})+2\bm{c}_{j}^{\top}\bm{p}=\phi_{j}\quad\forall j\in[J]\end{array}\right\}.$$

$$\mathcal{Q}_{m}=\left\{\begin{bmatrix}\bm{p}\\ 1\end{bmatrix}\begin{bmatrix}\bm{p}\\ 1\end{bmatrix}^{\top}:\bm{p}\in\mathcal{F}_{m}\right\}\;\textup{ and }\;\mathcal{R}_{m}=\left\{\begin{bmatrix}\bm{Q}&\bm{p}\\ \bm{p}^{\top}&1\end{bmatrix}\in\mathcal{C}(\mathcal{K}\times\mathbb{R}_{+}):\begin{array}[]{l}\displaystyle\bm{A}\bm{p}=\bm{b}\\ \operatorname*{diag}(\bm{A}\bm{Q}\bm{A}^{\top})=\bm{b}\circ\bm{b}\\ \displaystyle\textup{tr}(\bm{C}_{j}\bm{Q})+2\bm{c}_{j}^{\top}\bm{p}=\phi_{j}\;\;\forall j\in\mathcal{T}_{m}\end{array}\right\}$$

$$\left[\begin{array}{cc}\mathbf{Q}& \mathbf{p}\\ {\mathbf{p}}^{\top }& 1\end{array}\right]=\sum _{s\in \mathcal{S}}\left[\begin{array}{c}{\mathbf{\zeta }}_s\\ {\eta }_s\end{array}\right]{\left[\begin{array}{c}{\mathbf{\zeta }}_s\\ {\eta }_s\end{array}\right]}^{\top }=\sum _{s\in {\mathcal{S}}_{+}}{\eta }_s^2\left[\begin{array}{c}{\mathbf{\zeta }}_s/{\eta }_s\\ 1\end{array}\right]{\left[\begin{array}{c}{\mathbf{\zeta }}_s/{\eta }_s\\ 1\end{array}\right]}^{\top }+\sum _{s\in {\mathcal{S}}_0}\left[\begin{array}{c}{\mathbf{\zeta }}_s\\ 0\end{array}\right]{\left[\begin{array}{c}{\mathbf{\zeta }}_s\\ 0\end{array}\right]}^{\top }.$$

$$\text{tr}({\mathbf{C}}_j\mathbf{Q})+2{\mathbf{c}}_j^{\top }\mathbf{p}={\varphi }_j\forall j\in {\mathcal{T}}_m\setminus {\mathcal{T}}_{m-1}$$

$$({\mathbf{\zeta }}_s/{\eta }_s{)}^{\top }{\mathbf{C}}_j({\mathbf{\zeta }}_s/{\eta }_s)+2{\mathbf{c}}_j^{\top }({\mathbf{\zeta }}_s/{\eta }_s)={\varphi }_j\forall j\in {\mathcal{T}}_m\setminus {\mathcal{T}}_{m-1}.$$

$${\varphi }_j$$

$$=\sum _{s\in {\mathcal{S}}_{+}}{\eta }_s^2\left[({\mathbf{\zeta }}_s/{\eta }_s{)}^{\top }{\mathbf{C}}_j({\mathbf{\zeta }}_s/{\eta }_s)+2{\mathbf{c}}_j^{\top }({\mathbf{\zeta }}_s/{\eta }_s)\right]+\sum _{s\in {\mathcal{S}}_0}{\mathbf{\zeta }}_s^{\top }{\mathbf{C}}_j{\mathbf{\zeta }}_s$$

$$=\sum _{s\in {\mathcal{S}}_{+}}{\eta }_s^2\left[({\mathbf{\zeta }}_s/{\eta }_s{)}^{\top }{\mathbf{C}}_j({\mathbf{\zeta }}_s/{\eta }_s)+2{\mathbf{c}}_j^{\top }({\mathbf{\zeta }}_s/{\eta }_s)\right].$$

$$\begin{array}[]{cl}\min&\|\bm{X}-\bm{H}\bm{U}^{\top}\|_{F}^{2}\\ \textnormal{s.t.}&\bm{H}\in\mathbb{R}_{+}^{D\times K},\;\bm{U}\in\mathbb{R}_{+}^{N\times K}\\ &\bm{U}^{\top}\bm{U}=\mathbb{I}.\end{array}$$

$$\begin{array}[]{cl}\min&\textup{tr}\left(\bm{X}^{\top}\bm{X}-2\bm{X}\bm{U}\bm{H}^{\top}+\bm{H}^{\top}\bm{H}\right)\\ \textnormal{s.t.}&\bm{H}\in\mathbb{R}_{+}^{D\times K},\;\bm{U}\in\mathbb{R}_{+}^{N\times K}\\ &\bm{U}^{\top}\bm{U}=\mathbb{I}.\end{array}$$

$$\mathcal{W}(\mathcal{B},K)=\left\{\left(({\mathbf{p}}_i{)}_{i\in [K]},({\mathbf{Q}}_{ij}{)}_{i,j\in [K]}\right):\left[\begin{array}{cccc}{\mathbf{Q}}_{11}& \cdots & {\mathbf{Q}}_{1K}& {\mathbf{p}}_1\\ ⋮& \ddots & ⋮& ⋮\\ {\mathbf{Q}}_{K1}& \cdots & {\mathbf{Q}}_{KK}& {\mathbf{p}}_K\\ {\mathbf{p}}_1^{\top }& \cdots & {\mathbf{p}}_K^{\top }& 1\end{array}\right]\in \mathcal{C}\left({\mathcal{B}}^K\times {\mathbb{R}}_{+}\right)\right\},$$

$$\begin{array}[]{cll}\min&\displaystyle\textup{tr}(\bm{X}^{\top}\bm{X})+\sum_{i\in[K]}\textup{tr}(-2\bm{X}\bm{W}_{ii}+\bm{G}_{ii})\\ \textnormal{s.t.}&\left((\bm{p}_{i})_{i\in[K]},(\bm{Q}_{ij})_{i,j\in[K]}\right)\in\mathcal{W}\left(\mathbb{S}\mathbb{O}\mathbb{C}_{+}^{N+1}\times\mathbb{R}_{+}^{D},K\right)\\ &\bm{u}_{i}\in\mathbb{R}_{+}^{N},\;\bm{V}_{ij}\in\mathbb{R}_{+}^{N\times N},\;\bm{h}_{i}\in\mathbb{R}_{+}^{D},\;\bm{G}_{ij}\in\mathbb{R}_{+}^{D\times D},\;\bm{W}_{ij}\in\mathbb{R}_{+}^{N\times D}&\forall i,j\in[K]\\ &\bm{p}_{i}=\begin{bmatrix}\bm{u}_{i}\\ 1\\ \bm{h}_{i}\end{bmatrix},\;\bm{Q}_{ij}=\begin{bmatrix}\bm{V}_{ij}&\bm{u}_{i}&\bm{W}_{ij}\\ \bm{u}_{j}^{\top}&1&\bm{h}_{j}^{\top}\\ \bm{W}_{ji}^{\top}&\bm{h}_{i}&\bm{G}_{ij}\end{bmatrix}&\forall i,j\in[K]\\ &\textup{tr}(\bm{V}_{ii})=1&\forall i\in[K]\\ &\textup{tr}(\bm{V}_{ij})=0&\forall i,j\in[K]:i\neq j.\end{array}$$

$$\begin{array}[]{cll}\min&\displaystyle\textup{tr}(\bm{X}^{\top}\bm{X})-2\sum_{i\in[K]}\textup{tr}(\bm{X}\bm{u}_{i}\bm{h}_{i}^{\top})+\sum_{i\in[K]}\textup{tr}(\bm{h}_{i}\bm{h}_{i}^{\top})\\ \textnormal{s.t.}&\bm{h}_{i}\in\mathbb{R}_{+}^{D},\;\bm{u}_{i}\in\mathbb{R}_{+}^{N}&\forall i\in[K]\\ &\bm{u}_{i}^{\top}\bm{u}_{i}=1&\forall i\in[K]\\ &\bm{u}_{i}^{\top}\bm{u}_{j}=0&\forall i,j\in[K]:i\neq j.\end{array}$$

$${\mathbf{p}}_i=\left[\begin{array}{c}{\mathbf{u}}_i\\ t_i\\ {\mathbf{h}}_i\end{array}\right]\in \mathbb{SO}{\mathbb{C}}_{+}^{N+1}\times {\mathbb{R}}_{+}^D\forall i\in [K].$$

$$\mathcal{F}_{0}=\left\{\bm{p}\in\mathcal{K}:t_{i}=1\quad\forall i\in[K]\right\}\quad\textup{and}\quad\mathcal{F}_{1}=\mathcal{F}=\left\{\bm{p}\in\mathcal{F}_{0}:\begin{array}[]{ll}\bm{u}_{i}^{\top}\bm{u}_{i}=1&\forall i\in[K]\\ \bm{u}_{i}^{\top}\bm{u}_{j}=0&\forall i,j\in[K]:i\neq j\end{array}\right\}.$$

$${\overline{\mathbf{p}}}_i=\left[\begin{array}{c}{\overline{\mathbf{u}}}_i\\ 1\\ 0\end{array}\right]\in \mathbb{SO}{\mathbb{C}}_{+}^{N+1}\times {\mathbb{R}}_{+}^D\forall i\in [K],$$

$$\underset{\mathbf{p}\in {\mathcal{F}}_0}{\max }\left\{{\mathbf{u}}_i^{\top }{\mathbf{u}}_i\right\}=1\forall i\in [K],$$

$$\underset{\mathbf{p}\in {\mathcal{F}}_0}{\min }\left\{{\mathbf{u}}_i^{\top }{\mathbf{u}}_j\right\}=0\forall i,j\in [K]:i\ne j,$$

$$u_{in}>0⟹u_{jn}=0\text{ and }{u}_{jn}>0⟹u_{in}=0\forall n\in [N].$$

$${\mathbf{V}}_{ij}={\mathbf{u}}_i{\mathbf{u}}_j^{\top },{\mathbf{W}}_{ij}={\mathbf{u}}_i{\mathbf{h}}_j^{\top },\text{ and }{\mathbf{G}}_{ij}={\mathbf{h}}_i{\mathbf{h}}_j^{\top }\forall i,j\in [K].$$

$${\mathbf{Q}}_{ij}={\mathbf{p}}_i{\mathbf{p}}_j^{\top }=\left[\begin{array}{ccc}{\mathbf{V}}_{ij}& {\mathbf{u}}_i& {\mathbf{W}}_{ij}\\ {\mathbf{u}}_j^{\top }& 1& {\mathbf{h}}_j^{\top }\\ {\mathbf{W}}_{ji}^{\top }& {\mathbf{h}}_i& {\mathbf{G}}_{ij}\end{array}\right]\forall i,j\in [K].$$

$$\mathcal{R}=\left\{\begin{bmatrix}\bm{Q}_{11}&\cdots&\bm{Q}_{1K}&\bm{p}_{1}\\ \vdots&\ddots&\vdots&\vdots\\ \bm{Q}_{K1}&\cdots&\bm{Q}_{KK}&\bm{p}_{K}\\ \bm{p}_{1}^{\top}&\cdots&\bm{p}_{K}^{\top}&1\\ \end{bmatrix}\in\mathcal{C}(\mathcal{K}\times\mathbb{R}_{+}):\begin{array}[]{ll}\bm{p}_{i}=\begin{bmatrix}\bm{u}_{i}\\ 1\\ \bm{h}_{i}\end{bmatrix},\;\bm{Q}_{ij}=\begin{bmatrix}\bm{V}_{ij}&\bm{u}_{i}&\bm{W}_{ij}\\ \bm{u}_{j}^{\top}&1&\bm{h}_{j}^{\top}\\ \bm{W}_{ji}^{\top}&\bm{h}_{i}&\bm{G}_{ij}\end{bmatrix}&\forall i,j\in[K]\\ \textup{tr}(\bm{V}_{ii})=1&\forall i\in[K]\\ \textup{tr}(\bm{V}_{ij})=0&\forall i,j\in[K]:i\neq j\end{array}\right\}.$$

$$\begin{array}[]{cl}\min&\textup{tr}(\bm{X}^{\top}\bm{X}-\bm{X}^{\top}\bm{X}\bm{U}\bm{U}^{\top})\\ \textnormal{s.t.}&\bm{U}\in\mathbb{R}_{+}^{N\times K}\\ &\bm{U}^{\top}\bm{U}=\mathbb{I}.\end{array}$$

Full text

1 Introduction

Given an input set of data points, cluster analysis endeavors to discover a fixed number of disjoint clusters so that the data points in the same cluster are closer to each other than to those in other clusters. Cluster analysis is fundamental to a wide array of applications in science, engineering, economics, psychology, marketing, etc. [13, 14]. One of the most popular approaches for cluster analysis is $K$-means clustering [13, 18, 20]. The goal of $K$-means clustering is to partition the data points into $K$ clusters so that the sum of squared distances to the respective cluster centroids is minimized. Formally, $K$-means clustering seeks for a solution to the mathematical optimization problem

[TABLE]

Here, $\bm{x}_{1},\dots,\bm{x}_{N}$ are the input data points, while $\mathcal{P}_{1},\dots,\mathcal{P}_{K}\subseteq\{1,\dots,N\}$ are the output clusters. The vectors $\bm{c}_{1},\dots,\bm{c}_{K}\in\mathbb{R}^{D}$ in (1) determine the cluster centroids, while the constraints on the last row of (1) ensure that the subsets $\mathcal{P}_{1},\dots,\mathcal{P}_{K}$ constitute a partition of the set $\{1,\dots,N\}$.

Due to its combinatorial nature, the $K$-means clustering problem (1) is generically NP-hard [2]. A popular solution scheme for this intractable problem is the heuristic algorithm developed by Lloyd [18]. The algorithm initializes by randomly selecting $K$ cluster centroids. It then proceeds by alternating between the assignment step and the update step. In the assignment step the algorithm designates each data point to the closest centroid, while in the update step the algorithm determines new cluster centroids according to current assignment.

Another popular solution approach arises in the form of convex relaxation schemes [23, 4, 25]. In this approach, tractable semidefinite programming (SDP) lower bounds for (1) are derived. Solutions of these optimization problems are then transformed into cluster assignments via well-constructed rounding procedures. Such convex relaxation schemes have a number of theoretically appealing properties. If the data points are supported on $K$ disjoint balls then exact recovery is possible with high probability whenever the distance between any two balls is sufficiently large [4, 12]. A stronger model-free result is achievable if the cardinalities of the clusters are prescribed to the problem [25].

A closely related problem is the non-negative matrix factorization with orthogonality constraints (ONMF). Given an input data matrix $\bm{X}$, the ONMF problem seeks for non-negative matrices $\bm{F}$ and $\bm{U}$ so that both the product $\bm{F}\bm{U}^{\top}$ is close to $\bm{X}$ in view of the Frobenius norm and the orthogonality constraint $\bm{U}^{\top}\bm{U}=\mathbb{I}$ is satisfied. Although ONMF is not precisely equivalent to $K$-means, solutions to this problem have clustering property [9, 17, 10, 15]. In [24], it is shown that the ONMF problem is in fact equivalent to a weighted variant of the $K$-means clustering problem.

In this paper, we attempt to obtain equivalent convex reformulations for the ONMF and $K$-means clustering problems. To derive these reformulations, we adapt the results by Burer and Dong [7] who show that any (non-convex) quadratically constrained quadratic program (QCQP) can be reformulated as a linear program over the convex cone of completely positive matrices. The resulting optimization problem is called a generalized completely positive program. Such a transformation does not immediately mitigate the intractability of the original problem, since solving a generic completely positive program is NP-hard. However, the complexity of the problem is now entirely absorbed in the cone of completely positive matrices which admits tractable semidefinite representable outer approximations [22, 8, 16]. Replacing the cone with these outer approximations gives rise to SDP relaxations of the original problem that in principle can be solved efficiently.

As byproducts of our derivations, we identify a new condition that makes the ONMF and the $K$-means clustering problems equivalent and we obtain new SDP relaxations for the $K$-means clustering problem that are tighter than the well-known relaxation proposed by Peng and Wei [23]. The contributions of this paper can be summarized as follows.

We disclose a new connection between ONMF and $K$-means clustering. We show that $K$-means clustering is equivalent to ONMF if an additional requirement on the binarity of solution to the latter problem is imposed. This amends the previous incorrect result by Ding et al. [9, Section 2] and Li and Ding [17, Theorem 1] who claimed that both problems are equivalent.111To our best understanding, they have shown only one of the implications that establish an equivalence. 2. 2.

We derive exact conic programming reformulations for the ONMF and $K$-means clustering problems that are principally amenable to numerical solutions. To our best knowledge, we are the first to obtain equivalent convex reformulations for these problems. 3. 3.

In view of the equivalent convex reformulation, we derive tighter SDP relaxations for the $K$-means clustering problem whose solutions can be used to construct high quality estimates of the cluster assignment. 4. 4.

We devise a new approximation algorithm for the $K$-means clustering problem that leverages the improved relaxation and numerically highlight its superiority over the state-of-the-art SDP approximation scheme by Mixon et al. [21] and the Lloyd’s algorithm.

The remainder of the paper is structured as follows. In Section 2, we present a theorem for reformulating the QCQPs studied in the paper as generalized completely positive programs. In Section 3, we derive a conic programming reformulation for the ONMF problem. We extend this result to the setting of $K$-means clustering in Section 4. In Section 5, we develop SDP relaxations and design a new approximation algorithm for $K$-means clustering. Finally, we empirically assess the performance of our proposed algorithm in Section 6.

Notation:

For any $K\in\mathbb{N}$, we define $[K]$ as the index set $\{1,\dots,K\}$. We denote by $\mathbb{I}$ the identity matrix and by $\mathbf{e}$ the vector of all ones. We also define $\mathbf{e}_{i}$ as the $i$-th canonical basis vector. Their dimensions will be clear from the context. The trace of a square matrix $\bm{M}$ is denoted as $\textup{tr}(\bm{M})$. We define $\operatorname*{diag}(\bm{v})$ as the diagonal matrix whose diagonal components comprise the entries of $\bm{v}$. For any non-negative vector $\bm{v}\in\mathbb{R}_{+}^{K}$, we define the cardinality of all positive components of $\bm{v}$ by $\#\bm{v}=|\{i\in[K]:v_{i}>0\}|$. For any matrix $\bm{M}\in\mathbb{R}^{M\times N}$, we denote by $\bm{m}_{i}\in\mathbb{R}^{M}$ the vector that corresponds to the $i$-th column of $\bm{M}$. The set of all symmetric matrices in $\mathbb{R}^{K\times K}$ is denoted as $\mathbb{S}^{K}$, while the cone of positive semidefinite matrices in $\mathbb{R}^{K\times K}$ is denoted as $\mathbb{S}_{+}^{K}$. The cone of completely positive matrices over a set $\mathcal{K}$ is denoted as $\mathcal{C}(\mathcal{K})=\textup{cl}\textnormal{conv}\{\bm{x}\bm{x}^{\top}:\bm{x}\in\mathcal{K}\}$. For any $\bm{Q},\bm{R}\in\mathbb{S}^{K}$ and any closed convex cone $\mathcal{C}$, the relations $\bm{Q}\succeq\bm{R}$ and $\bm{Q}\succeq_{\mathcal{C}}\bm{R}$ denote that $\bm{Q}-\bm{R}$ is an element of $\mathbb{S}_{+}^{K}$ and $\mathcal{C}$, respectively. The $(K+1)$-dimensional second-order cone is defined as $\mathbb{S}\mathbb{O}\mathbb{C}^{K+1}=\{(\bm{x},t)\in\mathbb{R}^{K+1}:\|\bm{x}\|\leq t\}$, where $\|\bm{x}\|$ denotes the 2-norm of the vector $\bm{x}$. We denote by $\mathbb{S}\mathbb{O}\mathbb{C}^{K+1}_{+}=\mathbb{S}\mathbb{O}\mathbb{C}^{K+1}\cap\mathbb{R}^{K+1}_{+}$ the intersection of the $K+1$-dimensional second-order cone and the non-negative orthant.

2 Completely Positive Programming Reformulations of QCQPs

To derive the equivalent completely positive programming reformulations in the subsequent sections, we first generalize the results in [7, Theorem 1] and [6, Theorem 3]. Consider the (nonconvex) quadratically constrained quadratic program (QCQP) given by

[TABLE]

Here, $\mathcal{K}\subseteq\mathbb{R}^{D}$ is a closed convex cone, while $\bm{A}\in\mathbb{R}^{I\times D}$, $\bm{b}\in\mathbb{R}^{I}$, $\bm{C}_{0},\bm{C}_{j}\in\mathbb{S}^{D}$, $\bm{c}_{0},\bm{c}_{j}\in\mathbb{R}^{D}$, $\phi_{j}\in\mathbb{R}$, $j\in[J]$, are the respective input problem parameters. We define the feasible set of problem (2) as

[TABLE]

and the recession cone of the linear constraint system as $\mathcal{F}^{\infty}:=\{\bm{d}\in\mathcal{K}:\bm{A}\bm{d}=\bm{0}\}$. We further define the following subsets of $\mathcal{C}(\mathcal{K}\times\mathbb{R}_{+})$:

[TABLE]

A standard result in convex optimization enables us to reformulate the QCQP (2) as the linear convex program

[TABLE]

Recently, Burer [6] showed that, in the absence of quadratic constraints in $\mathcal{F}$, the set $\textup{cl}\textnormal{conv}\left(\mathcal{Q}\right)$ is equal to the intersection of a polynomial size linear constraint system and a generalized completely positive cone. In [7], Burer and Dong showed that such a reformulation is achievable albeit more cumbersome in the presence of generic quadratic constraints in $\mathcal{F}$. Under some additional assumptions about the structure of the quadratic constraints, one can show that the set $\textup{cl}\textnormal{conv}\left(\mathcal{Q}\right)$ is amenable to a much simpler completely positive reformulation (see [7, Theorem 1] and [6, Theorem 3]). Unfortunately, these assumptions are too restrictive to reformulate the quadratic programming instances we study in this paper. To that end, the following theorem provides the required extension that will enable us to derive the equivalent completely positive programs.

Theorem 1**.**

Suppose there exists an increasing sequence of index sets $\mathcal{T}_{0}=\emptyset\subseteq\mathcal{T}_{1}\subseteq\mathcal{T}_{2}\subseteq\cdots\subseteq\mathcal{T}_{M}=[J]$ with the corresponding structured feasible sets

[TABLE]

such that for every $m\in[M]$ we have

[TABLE]

and there exists a vector $\overline{\bm{p}}\in\mathcal{F}$ such that

[TABLE]

Then, $\textup{cl}\textnormal{conv}\left(\mathcal{Q}\right)$ coincides with

[TABLE]

Theorem 1 constitutes a generalization of the combined results of [7, Theorem 1] and [6, Theorem 3], which we state in the following proposition.

Proposition 1**.**

*Let $\mathcal{L}=\{\bm{p}\in\mathcal{K}:\bm{A}\bm{p}=\bm{b}\}$. Suppose $\phi_{j}=\min_{\bm{p}\in\mathcal{L}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$, and both $\min_{\bm{p}\in\mathcal{L}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$ and $\max_{\bm{p}\in\mathcal{L}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$ are finite for all $j\in[J]$. If there exists $\overline{\bm{p}}\in\mathcal{F}$ such that $\bm{d}^{\top}(\bm{C}_{j}\overline{\bm{p}}+\bm{c}_{j})=0$ for all $\bm{d}\in\mathcal{F}^{\infty}$ and $j\in[J]$, then $\textup{cl}\textnormal{conv}\left(\mathcal{Q}\right)$ coincides with $\mathcal{R}$. *

To see this, assume that all conditions in Proposition 1 are satisfied. Then, setting $M=1$ and $\mathcal{T}_{1}=[J]$, we find that the condition (6) in Theorem 1 is satisfied. Next, for every $j\in[J]$, the finiteness of both $\min_{\bm{p}\in\mathcal{L}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$ and $\max_{\bm{p}\in\mathcal{L}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$ implies that $\bm{d}^{\top}\bm{C}_{j}\bm{d}=0$ for all $\bm{d}\in\mathcal{F}^{\infty}$. Combining this with the last condition in Proposition 1, we find that there there exists a vector $\overline{\bm{p}}\in\mathcal{F}$ such that $\bm{d}^{\top}\bm{C}_{j}\bm{d}+2\bm{d}^{\top}(\bm{C}_{j}\overline{\bm{p}}+\bm{c}_{j})=0$ for all $\bm{d}\in\mathcal{F}^{\infty}$ and $j\in[J]$. Thus, all conditions in Theorem 1 are indeed satisfied.

In the remainder of the section, we define the sets

[TABLE]

for $m\in[M]\cup\{0\}$. The proof of Theorem 1 relies on the following lemma, which is established in the first part of the proof of [6, Theorem 3].

Lemma 1**.**

Suppose there exists a vector $\overline{\bm{p}}\in\mathcal{F}$ such that $\bm{d}^{\top}\bm{C}_{j}\bm{d}+2\bm{d}^{\top}(\bm{C}_{j}\overline{\bm{p}}+\bm{c}_{j})=0$ for all $\bm{d}\in\mathcal{F}^{\infty}$ and $j\in[J]$, then we have $\textnormal{conv}(\mathcal{Q}_{m})+\textnormal{cone}(\mathcal{Q}^{\infty})\subseteq\textup{cl}\textnormal{conv}(\mathcal{Q}_{m})$ for all $m\in[M]$.

Using this lemma, we are now ready to prove Theorem 1.

Proof of Theorem 1.

The proof follows if $\textup{cl}\textnormal{conv}(\mathcal{Q}_{m})=\mathcal{R}_{m}$ for all $m\in[M]$. By construction, we have $\textup{cl}\textnormal{conv}(\mathcal{Q}_{m})\subseteq\mathcal{R}_{m}$, $m\in[M]$. It thus remains to prove the converse inclusions. By Lemma 1, it suffices to show that $\mathcal{R}_{m}\subseteq\textnormal{conv}(\mathcal{Q}_{m})+\textnormal{cone}(\mathcal{Q}^{\infty})$ for all $m\in[M]$. We proceed via induction. The base case for $m=0$ follows from [6, Theorem 1]. Assume now that $\mathcal{R}_{m-1}\subseteq\textnormal{conv}(\mathcal{Q}_{m-1})+\textnormal{cone}(\mathcal{Q}^{\infty})$ holds for a positive index $m-1<M$. We will show that this implies $\mathcal{R}_{m}\subseteq\textnormal{conv}(\mathcal{Q}_{m})+\textnormal{cone}(\mathcal{Q}^{\infty})$. To this end, consider the following completely positive decomposition of an element of $\mathcal{R}_{m}$:

[TABLE]

Here, $\mathcal{S}_{+}=\{s\in\mathcal{S}:\eta_{s}>0\}$ and $\mathcal{S}_{0}=\{s\in\mathcal{S}:\eta_{s}=0\}$, where $\mathcal{S}$ is a finite index set. By our induction hypothesis, we have $\bm{\zeta}_{s}/\eta_{s}\in\mathcal{F}_{m-1}$, $s\in\mathcal{S}_{+}$, and $\bm{\zeta}_{s}\in\mathcal{F}^{\infty}$, $s\in\mathcal{S}_{0}$. The proof thus follows if the constraints

[TABLE]

in $\mathcal{R}_{m}$ imply

[TABLE]

Indeed, for every $j\in\mathcal{T}_{m}\setminus\mathcal{T}_{m-1}$, the decomposition (9) yields

[TABLE]

Here, the last equality follows from our assumption that there exists a vector $\overline{\bm{p}}\in\mathcal{F}$ such that $\bm{d}^{\top}\bm{C}_{j}\bm{d}+2\bm{d}^{\top}(\bm{C}_{j}\overline{\bm{p}}+\bm{c}_{j})=0$ for all $\bm{d}\in\mathcal{F}^{\infty}$. Thus, $\bm{d}^{\top}\bm{C}_{j}\bm{d}=0$ for all $\bm{d}\in\mathcal{F}^{\infty}$. Next, since $\bm{\zeta}_{s}/\eta_{s}\in\mathcal{F}_{m-1}$, the $j$-th identity in (6) implies that $(\bm{\zeta}_{s}/\eta_{s})^{\top}\bm{C}_{j}(\bm{\zeta}_{s}/\eta_{s})+2\bm{c}_{j}^{\top}(\bm{\zeta}_{s}/\eta_{s})\geq\phi_{j}$ if $\phi_{j}=\displaystyle\min_{\bm{p}\in\mathcal{F}_{m-1}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$ or $(\bm{\zeta}_{s}/\eta_{s})^{\top}\bm{C}_{j}(\bm{\zeta}_{s}/\eta_{s})+2\bm{c}_{j}^{\top}(\bm{\zeta}_{s}/\eta_{s})\leq\phi_{j}$ if $\phi_{j}=\displaystyle\max_{\bm{p}\in\mathcal{F}_{m-1}}\bm{p}^{\top}\bm{C}_{j}\bm{p}+2\bm{c}_{j}^{\top}\bm{p}$. The proof thus follows since $\eta_{s}^{2}>0$ and $\sum_{s\in\mathcal{S}_{+}}\eta_{s}^{2}=1$. ∎

3 Orthogonal Non-Negative Matrix Factorization

In this section, we first consider the ONMF problem given by

[TABLE]

Here, $\bm{X}\in\mathbb{R}^{D\times N}$ is a matrix whose columns comprise $N$ data points $\{\bm{x}_{n}\}_{n\in[N]}$ in $\mathbb{R}^{D}$. We remark that problem (10) is generically intractable since we are minimizing a non-convex quadratic objective function over the Stiefel manifold [1, 3]. By expanding the Frobenius norm in the objective function and noting that $\bm{U}^{\top}\bm{U}=\mathbb{I}$, we find that problem (10) is equivalent to

[TABLE]

We now derive a convex reformulation for problem (11). We remark that this problem is still intractable due to non-convexity of the objective function and the constraint system. Thus, any resulting convex formulation will in general remain intractable. In the following, to reduce the clutter in our notation, we define the convex set

[TABLE]

where $\bm{p}_{i}\in\mathcal{B}$ and $\bm{Q}_{ij}\in\mathbb{R}^{(N+1+D)\times(N+1+D)}_{+}$, $i,j\in[K]$. Here, $\mathcal{B}$ is a given convex cone, $K$ is a positive integer, and $\mathcal{B}^{K}$ is the direct product of $K$ copies of $\mathcal{B}$.

Theorem 2**.**

Problem (11) is equivalent to the following generalized completely positive program:

[TABLE]

Proof.

By utilizing the notation for column vectors $\{\bm{u}_{i}\}_{i\in[K]}$ and $\{\bm{h}_{i}\}_{i\in[K]}$, we can reformulate problem (11) equivalently as the problem

[TABLE]

We now employ Theorem 1 to show the equivalence of problems (13) and (12). We first introduce an auxiliary decision variable $\bm{p}=(\bm{p}_{1},\ldots,\bm{p}_{K})$ that satisfies

[TABLE]

Let $M=1$ in Theorem 1 and set $\mathcal{K}=(\mathbb{S}\mathbb{O}\mathbb{C}_{+}^{N+1}\times\mathbb{R}^{D}_{+})^{K}$. We then define the structured feasible sets

[TABLE]

Note that for every $i\in[K]$, the constraints $\|\bm{u}_{i}\|_{2}\leq t_{i}$ and $t_{i}=1$ in $\mathcal{F}_{0}$ imply that the variables $\bm{u}_{i}$ and $t_{i}$ are bounded. Thus, the recession cone of $\mathcal{F}_{0}$ coincides with the set $\mathcal{F}^{\infty}=\left\{\bm{p}\in\mathcal{K}:\bm{u}_{i}=\bm{0},\;t_{i}=0\;\forall i\in[K]\right\}$. Next, we set the vector $\overline{\bm{p}}=(\overline{\bm{p}}_{1},\ldots,\overline{\bm{p}}_{K})\in\mathcal{F}$ in Theorem 1 to satisfy

[TABLE]

where the subvectors $\{\overline{\bm{u}}_{i}\}_{i\in[K]}$ are chosen to be feasible in (13). In view of the description of recession cone $\mathcal{F}^{\infty}$ and the structure of quadratic constraints in $\mathcal{F}$, one can readily verify that such a vector $\overline{\bm{p}}$ satisfies the condition (7) in Theorem 1. It remains to show that condition (6) is also satisfied. Indeed, we have

[TABLE]

since the constraints $\|\bm{u}_{i}\|_{2}\leq 1$, $i\in[K]$, are implied by $\mathcal{F}_{0}$, while equalities are attained whenever the 2-norm of each vector $\bm{u}_{i}$ is 1. Similarly, we find that

[TABLE]

since the constraints $\bm{u}_{i}\geq\bm{0}$, $i\in[K]$, are implied by $\mathcal{F}_{0}$, while equalities are attained whenever the solutions $\bm{u}_{i}$ and $\bm{u}_{j}$ satisfy the complementarity property:

[TABLE]

Thus, all conditions in Theorem 1 are satisfied.

Next, we introduce new matrix variables that represent a linearization of the quadratic variables, as follows:

[TABLE]

We also define an auxiliary decision variable $\bm{Q}=(\bm{Q}_{ij})_{i,j\in[K]}$ satisfying

[TABLE]

Using these new terms, we construct the set $\mathcal{R}$ in Theorem 1 as follows:

[TABLE]

By Theorem 1, this set coincides with $\textup{cl}\textnormal{conv}\left(\mathcal{Q}\right)$, where the set $\mathcal{Q}$ is defined as in (3). Thus, by linearizing the objective function using the matrix variables in (14), we find that the generalized completely positive program (12) is indeed equivalent to (11). This completes the proof. ∎

Let us now consider a special case of problem (10); if all components of $\bm{X}$ are non-negative, then we can reduce the problem into a simpler one involving only the decision matrix $\bm{U}$.

Lemma 2**.**

If $\bm{X}$ is a non-negative matrix then problem (10) is equivalent to the non-convex program

[TABLE]

Proof.

Solving the minimization over ${\bm{H}}\in\mathbb{R}_{+}^{D\times K}$ analytically in (11), we find that the solution ${\bm{H}}=\bm{X}\bm{U}$ is feasible and optimal. Substituting this solution into the objective function of (11), we arrive at the equivalent problem (15). This completes the proof. ∎

By employing the same reformulation techniques as in the proof of Theorem 2, we can show that problem (15) is amenable to an exact convex reformulation.

Proposition 2**.**

Problem (15) is equivalent to the following generalized completely positive program:

[TABLE]

4 $K$-means Clustering

Building upon the results from the previous sections, we now derive an exact generalized completely positive programming reformulation for the $K$-means clustering problem (1). To this end, we note that the problem can equivalently be solved via the following mixed-integer nonlinear program [11]:

[TABLE]

Here, $\bm{c}_{i}$ is the centroid of the $i$-th cluster, while $\bm{\pi}_{i}$ is the assignment vector for the $i$-th cluster, i.e., $\pi_{in}=1$ if and only if the data point $\bm{x}_{n}$ is assigned to the cluster $i$. The last constraint in (17) ensures that each data point is assigned to a cluster, while the constraint system in the penultimate row ensures that there are exactly $K$ clusters. We now show that we can solve the $K$-means clustering problem by solving a modified problem (15) with an additional constraint $\sum_{i\in[K]}\bm{u}_{i}\bm{u}_{i}^{\top}\mathbf{e}=\mathbf{e}$. To further simplify our notation we will employ the sets

[TABLE]

[TABLE]

in all reformulations in the remainder of this section.

Theorem 3**.**

The following non-convex program solves the $K$-means clustering problem:

[TABLE]

Proof.

We first observe that the centroids in (17) can be expressed as

[TABLE]

Substituting these terms into the objective function and expanding the squared norm yield

[TABLE]

Thus, (17) can be rewritten as

[TABLE]

For any feasible solution $(\bm{\pi}_{i})_{i\in[K]}$ to (18) we define the vectors $(\bm{u}_{i})_{i\in[K]}$ that satisfy

[TABLE]

We argue that the solution $(\bm{u}_{i})_{i\in[K]}$ is feasible to $\mathcal{Z}$ and yields the same objective value. Indeed, we have

[TABLE]

because $\bm{\pi}_{i}\in\{0,1\}^{N}$ and $\mathbf{e}^{\top}\bm{\pi}_{i}\geq 1$ for all $i\in[K]$. We also have

[TABLE]

and

[TABLE]

since the constraint $\sum_{i\in[K]}\bm{\pi}_{i}=\mathbf{e}$ in (18) ensures that each data point is assigned to at most $1$ cluster. Verifying the objective value of this solution, we obtain

[TABLE]

Thus, we conclude that problem $\mathcal{Z}$ constitutes a relaxation of (18).

To show that $\mathcal{Z}$ is indeed an exact reformulation, consider any feasible solution $(\bm{u}_{i})_{i\in[K]}$ to this problem. For any fixed $i,j\in[K]$, the complementary constraint $\bm{u}_{i}^{\top}\bm{u}_{j}=0$ in $\mathcal{Z}$ means that

[TABLE]

Thus, in view of the last constraint in $\mathcal{Z}$, we must have $\bm{u}_{i}\in\{0,{1}/{\bm{u}_{i}^{\top}\mathbf{e}}\}^{N}$ for every $i\in[K]$. Using this observation, we define the binary vectors $(\bm{\pi}_{i})_{i\in[K]}$ that satisfy

[TABLE]

For every $i\in[K]$, we find that $\mathbf{e}^{\top}\bm{\pi}_{i}\geq 1$ since $\bm{u}_{i}^{\top}\bm{u}_{i}=1$. Furthermore, we have

[TABLE]

Substituting the constructed solution $(\bm{\pi}_{i})_{i\in[K]}$ into the objective function of (18), we obtain

[TABLE]

Thus, any feasible solution to $\mathcal{Z}$ can be used to construct a feasible solution to (18) that yields the same objective value. Our previous argument that (18) is a relaxation of $\mathcal{Z}$ then implies that both problems are indeed equivalent. This completes the proof. ∎

Remark 1**.**

The constraint $\sum_{i\in[K]}\bm{u}_{i}\bm{u}_{i}^{\top}\mathbf{e}=\mathbf{e}$ in $\mathcal{Z}$ ensures that there are no fractional values in the resulting cluster assignment vectors $(\bm{\pi}_{i})_{i\in[K]}$. While the formulation (15) is only applicable for instances of ONMF problem with non-negative input data $\bm{X}$, the reformulation $\mathcal{Z}$ remains valid for any instances of $K$-means clustering problem, even if the input data matrix $\bm{X}$ contains negative components.

Remark 2**.**

In [9, Section 2] and [17, Theorem 1], it was claimed that the ONMF problem (15) is equivalent to the $K$-means clustering problem (1). Theorem 3 above amends this result by showing that both problems become equivalent if and only if the constraint $\sum_{i\in[K]}\bm{u}_{i}\bm{u}_{i}^{\top}\mathbf{e}=\mathbf{e}$ is added to (15).

Remark 3**.**

We can reformulate the objective function of $\mathcal{Z}$ as $\frac{1}{2}\textup{tr}\left(\bm{D}\sum_{i\in[K]}\bm{u}_{i}\bm{u}_{i}^{\top}\right)$, where $\bm{D}$ is the matrix with components $D_{pq}=\|\bm{x}_{p}-\bm{x}_{q}\|^{2}$, $p,q\in[N]$. To obtain this reformulation, define $\bm{Y}=\sum_{i\in[K]}\bm{u}_{i}\bm{u}_{i}^{\top}$. Then we have

[TABLE]

Here, the fourth equality holds because of the last constraint in $\mathcal{Z}$ which ensures that $\sum_{q\in[N]}Y_{pq}=1$ for all $p\in[N]$.

We are now well-positioned to derive an equivalent generalized completely positive program for the $K$-means clustering problem.

Theorem 4**.**

The following generalized completely positive program solves the $K$-means clustering problem:

[TABLE]

Proof.

We consider the following equivalent reformulation of $\mathcal{Z}$ with two additional strengthening constraint systems.

[TABLE]

Since $\bm{s}_{i}\geq\bm{0}$, the last constraint system in (19) implies that $\bm{u}_{i}\leq w_{i}\mathbf{e}$, while the penultimate constraint system ensures that $\bm{u}_{i}$ is a binary vector, i.e., $\bm{u}_{i}\in\{0,w_{i}\}^{N}$ for some $w_{i}\in\mathbb{R}_{+}$. Since any feasible solution to $\mathcal{Z}$ satisfies these conditions, we may thus conclude that the problems $\mathcal{Z}$ and (19) are indeed equivalent. As we will see below, the exactness of the generalized completely positive programming reformulation is reliant on these two redundant constraint systems.

We now repeat the same derivation steps as in the proof of Theorem 2. First, we introduce an auxiliary decision variable $\bm{p}=(\bm{p}_{i})_{i\in[K]}$, that satisfies

[TABLE]

We then set $\mathcal{K}=(\mathbb{S}\mathbb{O}\mathbb{C}_{+}^{N+1}\times\mathbb{R}_{+}^{N+1})^{K}$, and define the structured feasible sets

[TABLE]

[TABLE]

and $\mathcal{F}_{2}=\mathcal{F}=\left\{\bm{p}\in\mathcal{F}_{1}:\sum_{i\in[K]}\bm{u}_{i}\bm{u}_{i}^{\top}\mathbf{e}=\mathbf{e}\right\}$. Here, we find that the recession cone of $\mathcal{F}_{0}$ is given by

[TABLE]

Next, we set the vector $\overline{\bm{p}}=(\overline{\bm{p}}_{1},\ldots,\overline{\bm{p}}_{K})\in\mathcal{F}$ in Theorem 1 to satisfy

[TABLE]

where the subvectors $\{\overline{\bm{u}}_{i}\}_{i\in[K]}$, $\{\overline{\bm{s}}_{i}\}_{i\in[K]}$ and $\{\overline{w}_{i}\}_{i\in[K]}$ are chosen so that they are feasible in (19). In view of the description of the recession cone $\mathcal{F}^{\infty}$ and the structure of the quadratic constraints in $\mathcal{F}$, one can verify that such a vector $\overline{\bm{p}}$ satisfies the condition (7) in Theorem 1.

It remains to show that condition (6) is also satisfied. To this end, it is already verified in the proof of Theorem 2 that

[TABLE]

We now show that

[TABLE]

We first demonstrate that the constraint $\bm{u}_{i}+\bm{s}_{i}=w_{i}\mathbf{e}$ in (20) implies $\bm{u}_{i}\circ\bm{u}_{i}\leq w_{i}\bm{u}_{i}$. Indeed, since $\bm{s}_{i}\geq\bm{0}$, we have $w_{i}\mathbf{e}-\bm{u}_{i}\geq\bm{0}$. Applying a componentwise multiplication with the components of $\bm{u}_{i}\geq\bm{0}$ on the left-hand side, we arrive at the desired inequality. Thus, we find that each equation in (22) indeed holds, where equality is attained whenever $u_{in}=0$. Finally, we verify that

[TABLE]

Note that the constraint $\bm{u}_{i}\circ\bm{u}_{i}=w_{i}\bm{u}_{i}$ in (21) implies that $\bm{u}_{i}\in\{0,w_{i}\}^{N}$, while the constraint $\bm{u}_{i}^{\top}\bm{u}_{i}=1$ further implies that $\#\bm{u}_{i}w_{i}^{2}=1$. Moreover, the complementary constraint $\bm{u}_{i}^{\top}\bm{u}_{j}=0$ ensures that

[TABLE]

Thus, for any feasible vector $\bm{p}\in\mathcal{F}_{1}$, we have

[TABLE]

for some $k\in[K]$ such that $u_{kn}=w_{k}$. Thus, the equalities (23) indeed hold. In summary, we have shown that all conditions in Theorem 1 are satisfied.

We now introduce new variables, in addition to the ones described in (14), that linearize the quadratic terms, as follows:

[TABLE]

We further define an auxiliary decision variable $\bm{Q}_{ij}$, $i,j\in[K]$, that satisfy

[TABLE]

Using these new terms, we construct the set $\mathcal{R}$ in Theorem 1 as follows:

[TABLE]

Here, the last constraint system arises from squaring the left-hand sides of the equalities

[TABLE]

which correspond to the last constraint system in (19). Finally by linearizing the objective function using variables in (14) and (24), we arrive at the generalized completely positive program $\overline{\mathcal{Z}}$. This completes the proof. ∎

5 Approximation Algorithm for $K$-means Clustering

In this section, we develop a new approximation algorithm for $K$-means clustering. To this end, we observe that in the reformulation $\overline{\mathcal{Z}}$ the difficulty of the original problem is now entirely absorbed in the completely positive cone $\mathcal{C}(\cdot)$ which has been well studied in the literature [5, 6, 8]. Any such completely positive program admits the hierarchy of increasingly accurate SDP relaxations that are obtained by replacing the cone $\mathcal{C}(\cdot)$ with progressively tighter semidefinite-representable outer approximations [8, 16, 22]. For the generalized completely positive program $\overline{\mathcal{Z}}$, we employ the simplest outer approximation that is obtained by replacing the completely positive cone $\mathcal{C}\left((\mathbb{S}\mathbb{O}\mathbb{C}_{+}^{N+1}\times\mathbb{R}_{+}^{N+1})^{K}\times\mathbb{R}_{+}\right)$ in $\overline{\mathcal{Z}}$ with its coarsest outer approximation [26], given by the cone

[TABLE]

where

[TABLE]

If $\bm{M}$ has the structure of the large matrix in $\overline{\mathcal{Z}}$ then the constraint $\textup{tr}(\mathbb{J}_{i}\bm{M})\geq 0$ reduces to $\textup{tr}(\bm{V}_{ii})\leq 1$, which is redundant and can safely be omitted in view of the stronger equality constraint $\textup{tr}(\bm{V}_{ii})=1$ in $\overline{\mathcal{Z}}$. In this case, the outer approximation can be simplified to the cone of doubly non-negative matrices given by

[TABLE]

To further improve computational tractability, we relax the large semidefinite constraint into a simpler system of $K$ semidefinite constraints. We summarize our formulation in the following proposition.

Proposition 3**.**

The optimal value of the following SDP constitutes a lower bound on $Z^{\star}$.

[TABLE]

Proof.

Without loss of generality, we can assign the first data point $\bm{x}_{1}$ to the first cluster. The argument in the proof of Theorem 3 indicates that the assignment vector for the first cluster is given by

[TABLE]

Thus, the data point $\bm{x}_{1}$ is assigned to the first cluster if and only if the first element of $\bm{\pi}_{1}$ is equal to $1$, i.e., $1=\mathbf{e}_{1}^{\top}\bm{\pi}_{1}=\mathbf{e}_{1}^{\top}\bm{V}_{11}\mathbf{e}$. Henceforth, we shall add this constraint to $\overline{\mathcal{Z}}$. While the constraint is redundant for the completely positive program $\overline{\mathcal{Z}}$, it will cut-off any symmetric solution in the resulting SDP relaxation.

We now replace the generalized completely positive cone in $\overline{\mathcal{Z}}$ with the corresponding cone of doubly non-negative matrices, which yields the following SDP relaxation:

[TABLE]

Since all principal submatrices of the large matrix are also positive semidefinite, we can further relax the constraint to a more tractable system

[TABLE]

Next, we eliminate the constraints $\textup{tr}(\bm{V}_{ij})=0$, $i,j\in[K]:i\neq j$, from (25). As the other constraints and the objective function in the resulting formulation do not involve the decision variables $\bm{V}_{ij}$ and $\bm{Q}_{ij}$, for any $i,j\in[K]$ such that $i\neq j$, we can safely omit these decision variables. Finally, by renaming all double subscript variables, e.g., $\bm{Q}_{ii}$ to $\bm{Q}_{i}$, we arrive at the desired semidefinite program ${\mathcal{R}}_{0}$. This completes the proof.

∎

The symmetry breaking constraint $\mathbf{e}_{1}^{\top}\bm{V}_{1}\mathbf{e}=1$ in ${\mathcal{R}}_{0}$ ensures that the solution $\bm{V}_{1}$ will be different from any of the solutions $\bm{V}_{i}$, $i\geq 2$. Specifically, the constraint $\sum_{i\in[K]}\bm{V}_{i}\mathbf{e}=\mathbf{e}$ in ${\mathcal{R}}_{0}$ along with the aforementioned symmetry breaking constraint implies that $\mathbf{e}_{1}^{\top}\bm{V}_{i}\mathbf{e}=0$ for all $i\geq 2$. Thus, any rounding scheme that identifies the clusters using the solution $(\bm{V}_{i})_{i\in[K]}$ will always assign the data point $\bm{x}_{1}$ to the first cluster. It can be shown, however, that there exists a partially symmetric optimal solution to ${\mathcal{R}}_{0}$ with $\bm{V}_{2}=\dots=\bm{V}_{K}$. This enables us to derive a further simplification to ${\mathcal{R}}_{0}$.

Corollary 1**.**

Problem ${\mathcal{R}}_{0}$ is equivalent to the semidefinite program given by

[TABLE]

Proof.

Any feasible solution to ${\mathcal{R}}_{0}$ can be used to construct a feasible solution to $\overline{\mathcal{R}}_{0}$ with the same objective value, as follows:

[TABLE]

Conversely, any feasible solution to $\overline{\mathcal{R}}_{0}$ can also be used to construct a feasible solution to ${\mathcal{R}}_{0}$ with the same objective value:

[TABLE]

Thus, the claim follows.

∎

By eliminating the constraints $\operatorname*{diag}(\bm{W}_{i})=\bm{\psi}_{i},\;\bm{\gamma}_{i}+\bm{\eta}_{i}=\rho_{i}\mathbf{e},\;\operatorname*{diag}(\bm{W}_{i}+\bm{\Sigma}_{i}+2\bm{\Theta}_{i})+\beta_{i}\mathbf{e}-2\bm{\psi}_{i}-2\bm{\theta}_{i}=0$, $i=1,2$, from $\overline{\mathcal{R}}_{0}$ we obtain an even simpler SDP relaxation.

Corollary 2**.**

The optimal value of the following SDP constitutes a lower bound on $R_{0}^{\star}$:

[TABLE]

We remark that the formulation ${\mathcal{R}}_{1}$ is reminiscent of the well-known SDP relaxation for $K$-means clustering [23]:

[TABLE]

We now derive an ordering of the optimal values of problems $\mathcal{Z}$, ${\mathcal{R}}_{0}$, ${\mathcal{R}}_{1}$, and $\mathcal{R}_{2}$.

Theorem 5**.**

We have

[TABLE]

Proof.

The first and the second inequalities hold by construction. To prove the third inequality, consider any feasible solution $(\bm{W}_{1},\bm{W}_{2})$ to ${\mathcal{R}}_{1}$. Then, the solution $\bm{Y}=\bm{W}_{1}+\bm{W}_{2}$ is feasible to $\mathcal{R}_{2}$ and yields the same objective value, which completes the proof. ∎

Obtaining any estimations of the best cluster assignment using optimal solutions of problem $\mathcal{R}_{2}$ is a non-trivial endeavor. If we have exact recovery, i.e., $Z^{\star}=R_{2}^{\star}$, then an optimal solution of $\mathcal{R}_{2}$ assumes the form

[TABLE]

where $\bm{\pi}_{i}$ is the assigment vector for the $i$-th cluster. Such a solution $\bm{Y}$ allows for an easy identification of the clusters. If there is no exact recovery then a few additional steps need to be carried out. In [23], an approximate cluster assignment is obtained by solving exactly another $K$-means clustering problem on a lower dimensional data set whose computational complexity scales with $\mathcal{O}(N^{(K-1)^{2}})$. If the solution of the SDP relaxation $\mathcal{R}_{2}$ is close to the exact recovery solution (26), then the columns of the matrix $\bm{Y}\bm{X}$ will comprise denoised data points that are near to the respective optimal cluster centroids. In [21], this strengthened signal is leveraged to identify the clusters of the original data points.

The promising result portrayed in Theorem 5 implies that any well-constructed rounding scheme that utilizes the improved formulation ${\mathcal{R}}_{0}$ (or ${\mathcal{R}}_{1}$) will never generate inferior cluster assignments to the ones from schemes that employ the formulation $\mathcal{R}_{2}$. Our new SDP relaxation further inspires us to devise an improved approximation algorithm for the $K$-means clustering problem. The central idea of the algorithm is to construct high quality estimates of the cluster assignment vectors $(\bm{\pi}_{i})_{i\in[K]}$ using the solution $(\bm{V}_{i})_{i\in[K]}$, as follows:

[TABLE]

To eliminate any symmetric solutions, the algorithm gradually introduces symmetry breaking constraints $\mathbf{e}_{n_{i}}^{\top}\bm{V}_{i}\mathbf{e}=1$, $i\geq 2$, to ${\mathcal{R}}_{0}$, where the indices $n_{i}$, $i\geq 2$, are chosen judiciously. The main component of the algorithm runs in $K$ iterations and proceeds as follows. It first solves the problem ${\mathcal{R}}_{0}$ and records its optimal solution $(\bm{V}^{\star}_{i})_{i\in[K]}$. In each of the subsequent iterations $k=2,\dots,K$, the algorithm identifies the best unassigned data point $\bm{x}_{n}$ for the $k$-th cluster. Here, the best data point corresponds to the index $n$ that maximizes the quantity $\mathbf{e}_{n}^{\top}\bm{V}^{\star}_{k}\mathbf{e}$. For this index $n$, the algorithm then appends the constraint $\mathbf{e}_{n}^{\top}\bm{V}^{\star}_{k}\mathbf{e}=1$ to the problem ${\mathcal{R}}_{0}$, which breaks any symmetry in the solution $(\bm{V}_{i})_{i\geq k}$. The algorithm then solves the augmented problem and proceeds to the next iteration. At the end of the iterations, the algorithm assigns each data point $\bm{x}_{n}$ to the cluster $k$ that maximizes the quantity $\mathbf{e}_{n}^{\top}\bm{V}^{\star}_{k}\mathbf{e}$. The algorithm concludes with a single step of Lloyd’s algorithm. A summary of the overall procedure is given in Algorithm 1.

6 Numerical Results

In this section, we assess the performance of the algorithm described in Section 5. All optimization problems are solved with MOSEK v8 using the YALMIP interface [19] on a 16-core 3.4 GHz computer with 32 GB RAM.

We compare the performance of Algorithm 1 with the classical Lloyd’s algorithm and the approximation algorithm222MATLAB implementation of the algorithm is available at https://github.com/solevillar/kmeans_sdp. proposed in [21] on $50$ randomly generated instances of the $K$-means clustering problem. While our proposed algorithm employs the improved formulation ${\mathcal{R}}_{0}$ to identify the clusters, the algorithm in [21] utilizes the existing SDP relaxation $\mathcal{R}_{2}$.

We adopt the setting of [4] and consider $N$ data points in $\mathbb{R}^{D}$ supported on $K$ balls of the same radius $r$. We set $K=3$, $N=75$, and $r=2$, and run the experiment for $D=2,\dots,6$. All results are averaged over $50$ trials generated as follows. In each trial, we set the centers of the balls to $\bm{0}$, ${\mathbf{e}}/{\sqrt{D}}$, and ${c\mathbf{e}}/{\sqrt{D}}$, where the scalar $c$ is drawn uniformly at random from interval $[10,20]$. This setting ensures that the first two balls are always separated by unit distance irrespective of $D$, while the third ball is placed further with a distance $c$ from the origin. Next, we sample $N/K$ points uniformly at random from each ball. The resulting $N$ data points are then input to the three algorithms.

Table 1 reports the quality of cluster assignments generated from Algorithm 1 relative to the ones generated from the algorithm in [21] and the Lloyd’s algorithm. The mean in the table represents average percentage improvement of the true objective value from Algorithm 1 relative to other algorithms. The $p$th percentile is the value below which $p\%$ of these improvements may be found. We find that our proposed algorithm significantly outperforms both the other algorithms in view of the mean and the $95$th percentile statistics. We further observe that the improvements deteriorate as problem dimension $D$ increases. This should be expected as the clusters become more apparent in a higher dimension, which makes them easier to be identified by all the algorithms. The percentile statistics further indicate that while the other algorithms can generate extremely poor cluster assignments, our algorithm consistently produces high quality cluster assignments and rarely loses by more than $5\%$.

Acknowledgements.

The authors are thankful to the Associate Editor and two anonymous referees for their constructive comments and suggestions that led to substantial improvements of the paper. This research was supported by the National Science Foundation grant no. 1752125.