# A Novel Convex Clustering Method for High-Dimensional Data Using Semiproximal ADMM

Authors: Huangyue Chen; Lingchen Kong; Yan Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Clustering is an important ingredient of unsupervised learning; classical clustering methods include K-means clustering and hierarchical clustering. These methods may suffer from instability because of their tendency prone to sink into the local optimal solutions of the nonconvex optimization model. In this paper, we propose a new convex clustering method for high-dimensional data based on the sparse group lasso penalty, which can simultaneously group observations and eliminate noninformative features. In this method, the number of clusters can be learned from the data instead of being given in advance as a parameter. We theoretically prove that the proposed method has desirable statistical properties, including a finite sample error bound and feature screening consistency. Furthermore, the semiproximal alternating direction method of multipliers is designed to solve the sparse group lasso convex clustering model, and its convergence analysis is established without any conditions. Finally, the effectiveness of the proposed method is thoroughly demonstrated through simulated experiments and real applications.


---

## Body

## 1. Introduction

Clustering is an important ingredient of unsupervised learning. It assigns samples into different clusters by minimizing the differences in the same cluster and maximizing the differences between different clusters. As an exploratory data analysis technique, it has been applied in many fields, such as image processing, energy engineering, and social networks [1–3]. To date, a wide variety of clustering methods have been introduced, including classical clustering methods such as K-means clustering [4, 5] and hierarchical clustering [6, 7]. However, the convexity of the corresponding optimization models cannot be guaranteed in general, so their global optimal solutions are hard to find. In addition, the performances of these methods are strongly dependent on their initial settings. To overcome the aforementioned disadvantages, convex clustering (CC), which provides a convex optimization perspective toward the task of clustering via shrinkage or regularization techniques, has been considered by many researchers, e.g., [8–11]. Its corresponding optimization model is given by(1)minX∈ℝn×pfX≔12X−AF2+γ1∑ℓ∈ΘωℓXℓ1⋅−Xℓ2⋅q,where A is a given data matrix with n observations and p features; ⋅F denotes the Frobenius norm of the matrix; γ1≥0 is a tuning parameter that controls the balance between the model fit and the number of clusters; Θ=ℓ=ℓ1,ℓ21≤ℓ1<ℓ2≤n is the index set; ωℓ≥0,ℓ∈Θ are given weights that are generally chosen based on the given data matrix A; Xℓ1⋅X⋅ℓ1 is the ℓ1th row (column) of matrix X; and the most common choices of q are 1,2, and ∞, which ensure the convexity of model (1). Let X^ be an optimal solution of (1); then, X^i⋅=X^j⋅ indicates that the observations Ai⋅ and Aj⋅ are assigned to the same cluster. Note from (1) that if we set γ1=0, this will lead to the trivial solution X^=A, i.e., the partition of the observations into singletons. As γ1 increases, some rows of X^ become identical. If we set γ1 to be sufficiently large, all rows of X^ become identical, which indicates that all observations are lumped into a single cluster. Hence, the clustering path can be obtained by solving (1) with a range of tuning parameters, and then the number of clusters is determined by learning from the data rather than being given in advance. This fascinating characteristic of convex clustering has attracted increasing attention, see, e.g., [12–21].The theoretical results on cluster recovery of (1) with ωℓ=1 had been established in [12, 13]. Later, Sun et al. [14] extended these results to the general weighted convex clustering model. Furthermore, in [15], the authors analyzed the statistical properties of (1) with ωℓ=1, such as an unbiased estimator of the degrees of freedom and a finite sample bound for the prediction error. The large sample behavior of (1) with l1 fusion penalization is presented in [16]. In [17], the authors present conditions that guarantee that the convex clustering solution path recovers a tree and explicitly describes how the affinity parameters modulate the structure of the recovered tree. Meanwhile, some researchers concentrated on the solution algorithms of (1). For example, Chi and Lange [10] adopted the Alternating Direction Method of Multipliers (ADMM) and Alternating Minimization Algorithm (AMA) to solve (1). More recently, in [14, 18], the semismooth Newton-based augmented Lagrangian method was applied to (1) and outperformed ADMM and AMA in efficiency, scalability, and robustness. Moreover, the idea of convex clustering was applied to process missing data [19], binary data [20], and outlier detection [21].Although convex clustering (1) has desirable theoretical results and performs well computationally, Wang et al. [22] pointed out that the clustering performance may be destroyed in high-dimensional scenarios where the number of features is much more than observations. The noninformative features included in the clustering are the key reason behind why convex clustering (1) is disheartening. Hence, eliminating some noninformative features is indispensable for high-dimensional data clustering. This motivates us to propose a more reasonable convex clustering method that can perform cluster analysis and feature screening simultaneously.In this paper, we propose the Sparse Group Lasso Convex Clustering (SGLCC) method by adopting the sparse group lasso penalty [23]. The optimization problem is summarized as(2)minX∈ℝn×pfX+γ2∑j=1p1−αujX⋅j2+αX⋅j1,where the tuning parameter γ2 controls the number of informative features and u=u1,u2,…,upT∈ℝp is the weight vector. The form of the penalty term is much like the elastic-net penalty [24], and X⋅j2 and X⋅j1 denote the Euclidean norm and ℓ1 norm of vector X⋅j, respectively. The parameter 0≤α≤1 controls the balance between the group lasso [25, 26] and the lasso [27]. Obviously, convex clustering (1) is a special case of (2) if we set γ2=0. Unlike (1), model (2) can perform variable selection. When α=0, model (2) reduces to sparse convex clustering (SCC) [22]. SCC uses the group lasso penalty to promote group sparsity, that is, all components of the solution X^⋅j are either zero or nonzero. However, it cannot induce ubiquitous within-group sparsity in genetic data. For example, a biological pathway may be implicated in the progression of a particular type of cancer, but not all genes in the pathway need to be active [28]. The sparse group lasso penalty is designed to achieve such within-group sparsity. Therefore, we believe that the proposed clustering model (2) can significantly improve the clustering performance of model (1). It is worth noting that solving the optimization problem (2) is more challenging than convex clustering (1) and SCC. If the optimization algorithms utilized in [10, 22] are adopted to solve (2) directly, the computation efficiency may be poor because the subproblem might have no closed-form solution and the step length is too small. Hence, we developed the steps in this paper along with the theory, algorithm, and applications of SGLCC.The main contributions of the paper are summarized as follows:(1)
For the proposed SGLCC, we obtain a finite sample error bound and feature screening consistency under mild conditions. Our results are not only applicable to the more practical SGLCC model but also subsume the known results of (1) and SCC as special cases.(2)
We adopt the semiproximal alternating direction method of multipliers to solve (2), and its convergence analysis is established without any additional conditions. The subproblem per iteration of this algorithm has a closed-form solution and can be implemented efficiently.(3)
The experimental results on both synthetic and real datasets illustrate that SGLCC provides superior clustering performance and feature selection abilities to other clustering methods.The remainder of this paper is organized as follows. In Section2, we summarize some preliminaries and notations. We study the finite sample error bound and feature screening consistency of our proposed method in Section 3. In Section 4, we introduce an efficient optimization algorithm for solving (2). After that, in Section 5, we conduct numerical experiments to verify the effectiveness of the proposed method in synthetic datasets and four microarray datasets. Conclusions are given in Section 6. The proof of the main results can be found in the Supplementary Materials (available here).


## 2. Preliminaries

In this section, we introduce some preliminaries that are used later in the paper.For convenience, we start with some necessary notations. For a vectorx∈ℝn and q≥1, xq denotes the ℓq norm of x and 1n∈ℝn denotes a vector with all components equal to one. Let x=x1,x2,…,xnT. For a matrix X∈ℝn×m, Xi⋅X⋅j is the ith row (the jth column) of X, XF denotes the Frobenius norm of X, X2 denotes the spectral norm of X, and the vectorization of X is defined by vecX=X11,X21,…,Xn1,X12,…,XnmT. The symbols ο and ⊗ denote the Hadamard product and Kronecker product, respectively. The linear operator ℋ:ℝn×p⟶ℝΘ×p is given by(3)ℋX=Xℓ1⋅−Xℓ2⋅ℓ1,ℓ2∈Θ=HX.Theℓth row of H is denoted as Hℓ⋅=eℓ1−eℓ2T,ℓ=ℓ1,ℓ2∈Θ, and eℓ1 is a vector with the ℓ1th component equal to one and the rest equal to zero. The adjoint operator of ℋ is given by ℋ∗:ℝΘ×p⟶ℝn×p, that is, ℋ∗Y=HTY.Next, some definitions and results from convex analysis [29] are provided. Proximal mapping is important for designing optimization algorithms and has been well studied. For a proper closed convex function g:ℝn⟶−∞;+∞, the proximal mapping Proxtgx for g at any x∈ℝn with t>0 is defined by(4)Proxtgx=argminytgy+12x−y22.In particular, ifg is an indicator function of the set C (gx=0 if x∈C; otherwise gx=∞); then, the proximal mapping reduces to the projection operator onto C. The Fenchel conjugate of g is defined by(5)g∗y=supxx,y−gx.A key property that connects the proximal mapping and Fenchel conjugate is the so-called Moreau decomposition:(6)x=Proxtgx+tProxt−1g∗t−1x.For example, ifgx=xq with q≥1, from Moreau decomposition, we obtain(7)Proxtgx=x−tProxt−1g∗t−1x=x−Πtℬx,where ℬ=xxq∗≤1 and 1/q+1/q∗=1.


## 3. Statistical Properties

In this section, we study the finite sample error bound and feature screening consistency of sparse group lasso convex clustering (2). We start with the following conditions that facilitate the technical proofs.A1: assume that vecA=x∗+ε, where x∗=vecX∗∈ℝnp is the mean vector and ε∈ℝnp is an error vector of independent sub-Gaussian random variables with mean zero and variance δ2A2: only ss≪p features are informative; the informative feature set and its complementary set are then denoted as ℐ and ℐc, respectivelyThe conditionA1 implies that the data matrix A is composed of sub-Gaussian random variables and is commonly used in the literature, see, e.g., [15, 22]. In high-dimensional scenarios, only a few features are active in general. Thus, condition A2 has been widely used in high-dimensional data analysis, see, e.g., [22, 30]. For simplicity, as described in the literature [12, 13, 15, 22], we consider the case with ωℓ=1.


### 3.1. Bounds for Prediction Error

The following theorems provide the finite sample bounds for prediction error of model (2), and our theoretical results subsume the existing results for CC and SCC as special cases.Theorem 1.
Suppose that condition A1 is satisfied. LetX^ be the solution of model (2) with q=1. If γ1≥2δlogpΘ/n, then there exists positive constants b1 and b2 such that(8)Pc1npX^−X∗F2≤4γ1npvecHX∗1+2γ2αnpvecX∗1+c2+c3≥1−c4.Here,(9)c1=1+γ2α−1,c2=2δ21n+lognpn2p+2np,c3=γ21−αnpu22,c4=exp−minb1lognp,b2plognp+γ21−αu1+γ2α2δ2np+2pΘ.From Theorem1, we can observe that the average prediction error is bounded, and c2 tends to zero as n,p⟶∞. We know that vecX∗1=On and vecHX∗1=On2 by combining condition A2 with Θ=nn−1/2. Model (2) with q=1 has prediction consistency only if 0≤γ2<1/1−α, γ2=O1, γ21−α/npu12=o1, and nlogpΘ/p2=o1. To be specific, we find that c1>0,c1=O1 and 2γ2α/npvecX∗1 tend to zero as n,p⟶∞ by 0≤γ2<1/1−α and γ2=O1. If γ21−α/npu12=o1, then(10)0≤c3≤γ21−αnpu12=o1,and c4 decays to zero as n,p⟶∞. Additionally, by including γ1≥2δlogpΘ/n, the condition nlogpΘ/p2=o1 implies that 4γ1/npvecHX∗1=o1.The next theorem presents the finite sample bound for the prediction error of model (2) with q=2.Theorem 2.
Suppose that condition A1 is satisfied. LetX^ be the solution of model (2) with q=2. If γ1≥2δplogpΘ/n, then there exists positive constants b1 and b2 such that(11)Pc1npX^−X∗F2≤4γ1np∑ℓ∈ΘHℓ⋅X∗2+2γ2αnpvecX∗1+c2+c3≥1−c4.We only discuss the first term on the right-hand side of (11) in Theorem 2 because all other terms are the same as (8) in Theorem 1. We find that Hℓ⋅X∗2=O1 by condition A2. Note that Θ=nn−1/2, and thus ∑ℓ∈ΘHℓ⋅X∗2=On2. Suppose that nlogpΘ/p=o1, we obtain that 4γ1/np∑ℓ∈ΘHℓ⋅X∗2=o1 by combining γ1≥2δplogpΘ/n.Remark 1.
(1)
In high-dimensional scenarios, the number of featuresp is considerably larger than the number of observations n. If p=Onr,r>1, then conditions nlogpΘ/p2=o1 and nlogpΘ/p=o1 hold automatically, and this condition can be found in [15].(2)
This intuitive illustration implies that the conditionγ21−α/npu12=o1 requires that the weight u not be too large.(3)
Theorems1 and 2 subsume the known results for CC and SCC by taking γ2=0 and α=0, respectively. In addition, we obtain encouraging results provided in Corollary 1 when α=1.
The following corollary follows directly from Theorems1 and 2 by letting α=1.Corollary 1.
Suppose that conditions A1-A2 are satisfied. LetX^ be the solution of model (2) with α=1. We obtain the following statements.(1)
Whenq=1, if γ1≥2δlogpΘ/n, γ2=O1, and nlogpΘ/p2=o1, then, as n,p⟶∞,(12)P1npX^−X∗F2≤4γ1npvecHX∗1+2γ2npvecX∗1+c2⟶1.(2)
Whenq=2, if γ1≥2δplogpΘ/n, γ2=O1, and logpΘ/np=o1, then, as n,p⟶∞,(13)P1npX^−X∗F2≤4γ1np∑ℓ∈ΘHℓ⋅X∗2+2γ2npvecX∗1+c2⟶1.As we can see from Corollary1, SGLCC with α=1 has prediction consistency and only need some conditions that are easy to satisfy.


### 3.2. Feature Screening Consistency

Feature screening consistency is a desirable property in high-dimensional scenarios where many noninformative features need to be eliminated. This section demonstrates the asymptotic feature screening consistency of SGLCC.Theorem 3.
LetX^ be the solution of model (2) with q=1. Suppose that conditions A1-A2 are hold. If γ1≥2δlogpΘ/n, 4γ1/npvecHX∗1=o1, 0≤γ2<1/1−α,γ2=O1, and γ21−α/npu12=o1, then PX^⋅j2=0⟶1, for any j∈ℐc.Theorem 4.
LetX^ be the solution of model (2) with q=2. Suppose that conditions A1-A2 are hold. If γ1≥2δplogpΘ/n, 4γ1/np∑ℓ∈ΘHℓ⋅X∗2=o1, 0≤γ2<1/1−α,γ2=O1, and γ21−α/npu12=o1, then PX^⋅j2=0⟶1, for any j∈ℐc.Theorems3 and 4 give the asymptotic consistency of feature selection of SGLCC with q=1 and q=2, respectively. In this sense, SGLCC possesses the ability to correctly eliminate the noninformative features that distort clustering performance. Hence, we believe that the clustering performance of SGLCC is promising.


## 3.1. Bounds for Prediction Error

The following theorems provide the finite sample bounds for prediction error of model (2), and our theoretical results subsume the existing results for CC and SCC as special cases.Theorem 1.
Suppose that condition A1 is satisfied. LetX^ be the solution of model (2) with q=1. If γ1≥2δlogpΘ/n, then there exists positive constants b1 and b2 such that(8)Pc1npX^−X∗F2≤4γ1npvecHX∗1+2γ2αnpvecX∗1+c2+c3≥1−c4.Here,(9)c1=1+γ2α−1,c2=2δ21n+lognpn2p+2np,c3=γ21−αnpu22,c4=exp−minb1lognp,b2plognp+γ21−αu1+γ2α2δ2np+2pΘ.From Theorem1, we can observe that the average prediction error is bounded, and c2 tends to zero as n,p⟶∞. We know that vecX∗1=On and vecHX∗1=On2 by combining condition A2 with Θ=nn−1/2. Model (2) with q=1 has prediction consistency only if 0≤γ2<1/1−α, γ2=O1, γ21−α/npu12=o1, and nlogpΘ/p2=o1. To be specific, we find that c1>0,c1=O1 and 2γ2α/npvecX∗1 tend to zero as n,p⟶∞ by 0≤γ2<1/1−α and γ2=O1. If γ21−α/npu12=o1, then(10)0≤c3≤γ21−αnpu12=o1,and c4 decays to zero as n,p⟶∞. Additionally, by including γ1≥2δlogpΘ/n, the condition nlogpΘ/p2=o1 implies that 4γ1/npvecHX∗1=o1.The next theorem presents the finite sample bound for the prediction error of model (2) with q=2.Theorem 2.
Suppose that condition A1 is satisfied. LetX^ be the solution of model (2) with q=2. If γ1≥2δplogpΘ/n, then there exists positive constants b1 and b2 such that(11)Pc1npX^−X∗F2≤4γ1np∑ℓ∈ΘHℓ⋅X∗2+2γ2αnpvecX∗1+c2+c3≥1−c4.We only discuss the first term on the right-hand side of (11) in Theorem 2 because all other terms are the same as (8) in Theorem 1. We find that Hℓ⋅X∗2=O1 by condition A2. Note that Θ=nn−1/2, and thus ∑ℓ∈ΘHℓ⋅X∗2=On2. Suppose that nlogpΘ/p=o1, we obtain that 4γ1/np∑ℓ∈ΘHℓ⋅X∗2=o1 by combining γ1≥2δplogpΘ/n.Remark 1.
(1)
In high-dimensional scenarios, the number of featuresp is considerably larger than the number of observations n. If p=Onr,r>1, then conditions nlogpΘ/p2=o1 and nlogpΘ/p=o1 hold automatically, and this condition can be found in [15].(2)
This intuitive illustration implies that the conditionγ21−α/npu12=o1 requires that the weight u not be too large.(3)
Theorems1 and 2 subsume the known results for CC and SCC by taking γ2=0 and α=0, respectively. In addition, we obtain encouraging results provided in Corollary 1 when α=1.
The following corollary follows directly from Theorems1 and 2 by letting α=1.Corollary 1.
Suppose that conditions A1-A2 are satisfied. LetX^ be the solution of model (2) with α=1. We obtain the following statements.(1)
Whenq=1, if γ1≥2δlogpΘ/n, γ2=O1, and nlogpΘ/p2=o1, then, as n,p⟶∞,(12)P1npX^−X∗F2≤4γ1npvecHX∗1+2γ2npvecX∗1+c2⟶1.(2)
Whenq=2, if γ1≥2δplogpΘ/n, γ2=O1, and logpΘ/np=o1, then, as n,p⟶∞,(13)P1npX^−X∗F2≤4γ1np∑ℓ∈ΘHℓ⋅X∗2+2γ2npvecX∗1+c2⟶1.As we can see from Corollary1, SGLCC with α=1 has prediction consistency and only need some conditions that are easy to satisfy.


## 3.2. Feature Screening Consistency

Feature screening consistency is a desirable property in high-dimensional scenarios where many noninformative features need to be eliminated. This section demonstrates the asymptotic feature screening consistency of SGLCC.Theorem 3.
LetX^ be the solution of model (2) with q=1. Suppose that conditions A1-A2 are hold. If γ1≥2δlogpΘ/n, 4γ1/npvecHX∗1=o1, 0≤γ2<1/1−α,γ2=O1, and γ21−α/npu12=o1, then PX^⋅j2=0⟶1, for any j∈ℐc.Theorem 4.
LetX^ be the solution of model (2) with q=2. Suppose that conditions A1-A2 are hold. If γ1≥2δplogpΘ/n, 4γ1/np∑ℓ∈ΘHℓ⋅X∗2=o1, 0≤γ2<1/1−α,γ2=O1, and γ21−α/npu12=o1, then PX^⋅j2=0⟶1, for any j∈ℐc.Theorems3 and 4 give the asymptotic consistency of feature selection of SGLCC with q=1 and q=2, respectively. In this sense, SGLCC possesses the ability to correctly eliminate the noninformative features that distort clustering performance. Hence, we believe that the clustering performance of SGLCC is promising.


## 4. Algorithmic Design

In this section, we adopt the semiproximal alternating direction method of multipliers to solve (2), and its global convergence can be established. For simplicity, we focus on designing an algorithm to solve (2) with q=2. The same algorithmic design and implementation can be applied to cases q=1 and q=∞ without difficulty.


### 4.1. Optimality Conditions

By ignoring the terms withωℓ=0, model (2) can be rewritten as(14)minX∈ℝn×p12X−AF2+γ1∑ℓ∈ΕωℓXℓ1⋅−Xℓ2⋅2+γ2∑j=1p1−αujX⋅j2+αX⋅j1,where ℰ=ℓ=ℓ1,ℓ2ωℓ>0. For the convenience of mathematical expression, we define the linear operator J:ℝn×p⟶ℝℰ×p and its adjoint operator J∗:ℝℰ|×p⟶ℝn×p, respectively, by(15)JX=Xℓ1⋅−Xℓ2⋅ℓ1,ℓ2∈Ε=JXandJ∗Y=JTY,where Jℓ⋅=eℓ1−eℓ2T. Furthermore, model (13) can be reformulated as the following linearly constrained convex optimization problem with a separable structure:(16)minX,Y,Z,WFX,Y,Z,W≔12X−AF2+Q1Y+Q2Z+Q3W,s.t.JX−Y=0,X−Z=0,X−W=0,where Q1Y=γ1∑ℓ∈ℰωℓYℓ⋅2, Q2Z=γ21−α∑j=1pujZ⋅j2, and Q3W=γ2α∑j=1pX⋅j1. The Lagrangian function associated with (16) is defined by(17)LX,Y,Z,W;λ,μ,ξ=12X−AF2+Q1Y+Q2Z+Q3W+λ,JX−Y+μ,X−Z+ξ,X−W,where λ∈ℝℰ×p, μ∈ℝn×p, and ξ∈ℝn×p are Lagrangian multipliers. For a given parameter σ>0, the augmented Lagrangian function for (16) is defined by(18)LσX,Y,Z,W;λ,μ,ξ=LX,Y,Z,W;λ,μ,ξ+σ2JX−YF2+σ2X−ZF2+σ2X−WF2.The Karush–Kuhn–Tucker (KKT) conditions for (16) are given by(19)X−A+JTλ+μ+ξ=0,Y−ProxQ1Y+λ=0,Z−ProxQ2Z+μ=0,W−ProxQ3W+ξ=0,JX−Y=0,X−Z=0,X−W=0.Letℱ:=X,Y,Z,WJX−Y=0,X−Z=0,X−W=0 be a feasible set of (16). Obviously, ℱ≠∅ and FX,Y,Z,W≥0. Then, we know from Corollaries 28.2.2 and 28.3.1 in [29] that X^,Y^,Z^,W^ is an optimal solution for (16) if and only if there exists a Lagrangian multiplier λ^,μ^,ξ^ such that X^,Y^,Z^,W^,λ^,μ^,ξ^ satisfies the KKT conditions (19). Hence, we could design a stopping criterion based on the KKT conditions in terms of guaranteeing the optimality of a point generated by the algorithm proposed in Section 4.2.


### 4.2. Algorithm and Convergence Analysis

Although the optimization problem (16) is a convex optimization problem with linear constraints, it is difficult to minimize all four variables at the same time. Hence, we divide all four variables into two groups and update the two groups of variables alternately. According to the separable structure, we employ the semiproximal alternating direction method of multipliers (sPADMM) to solve (16) and view X as one block and V=Y,Z,W as another.Next, we discuss the iteration schemes of sPADMM for solving the optimization problem (16). Given the kth iteration Xk,Vk,λk,μk,ξk and two self-adjoint positive semidefinite operators P1 and P2, first update variables X and V as follows:(20)Xk+1=argminXLσX,Vk;λk,μk,ξk+σ2X−XkP12,Vk+1=argminVLσXk+1,V;λk,μk,ξk+σ2V−VkP22.Then, update Lagrangian multipliers:(21)λk+1=λk+τσJXk+1−Yk+1,μk+1=μk+τσXk+1−Zk+1,ξk+1=ξk+τσXk+1−Wk+1.Here,τ is step length and X−XkP12:=∑j=1pX⋅j−X⋅jkTP1X⋅j−X⋅jk. It is well known that the sPADMM reduces to the classical ADMM introduced by [31, 32] when P1=0 and P2=0. A comprehensive review of sPADMM can be found in Appendix B of [33].Remark 2.
The choice of the two self-adjoint positive semidefinite linear operatorsP1 and P2 depends greatly on the problem. The general principle is that both P1 and P2 should be as small as possible, while the corresponding subproblems are still relatively easy to solve.
Each block of (20) involves solving a convex optimization problem. Next, we show the closed-form solutions of each block. First, we deal with the first block in (20):(22)Xk+1=argminXLσX,Vk;λk,μk,ξk+σ2X−XkP12=argminX12X−AF2+σ2JX−Yk+λkσF2+σ2X−Zk+μkσF2+σ2X−Wk+ξkσF2+σ2X−XkP12.
TheX-subproblem is a smooth, strongly convex optimization that has a unique global minimizer Xk+1. Finding the minimizer Xk+1 is equivalent to finding the solution of the following linear system:(23)MX=A+σJTYk+Zk+Wk+P1Xk−JTλk−μk−ξk,where M=σJTJ+P1+1+2σIn. To solve this linear system efficiently, we design an appropriate P1. In our experiment, let P1=nIn−1n1nT−JTJ, which is a positive semidefinite matrix. Then, M=1+2σ+nσIn−σ1n1nT, and(24)M−1=11+2σ+nσIn+σ1+2σ1n1nT,by the Sherman–Morrison–Woodbury formula [34]. Hence, the closed-form solution for this linear system can be obtained as(25)Xk+1=M−1A+σJTYk+Zk+Wk+P1Xk−JTλk−μk−ξk.
For the second block in (20), we let the linear operator P2=0 since the linear operator of V in (16) is an identity operator that ensures the well-defined V-subproblem:(26)Vk+1=argminVLσXk+1,V;λk,μk,ξk=argminY,Z,WQ1Y+Q2Z+Q3W+σ2JXk+1+λkσ−YF2+σ2Xk+1+μkσ−ZF2+σ2Xk+1+ξkσ−WF2,which can be divided into three parts:(27)Yk+1=argminYQ1Y+σ2JXk+1+λkσ−YF2,Zk+1=argminZQ2Z+σ2Xk+1+μkσ−ZF2,Wk+1=argminWQ3W+σ2Xk+1+ξkσ−WF2.
Together with the definition of proximal mapping and Moreau decomposition, we show the closed-form solutions of the above subproblems. For anyℓ∈ℰ, the Y-subproblem has a closed-form solution that is given by(28)Yℓ⋅k+1=max1−γ1ωℓσJℓ⋅Xk+1+λℓ⋅k/σ2,0Jℓ⋅Xk+1+λℓ⋅kσ.
Similarly, for anyj=1,2,…,p, the solution of the Z-subproblem is given by(29)Z⋅jk+1=max1−γ2uj1−ασX⋅jk+1+μ⋅jk/σ2,0X⋅jk+1+μ⋅jkσ.
The last subproblem is a soft threshold operator for each column, that is, for anyj=1,2,…,p,(30)W⋅jk+1=signX⋅jk+1+ξ⋅jkσ°maxX⋅jk+1+ξ⋅jkσ−γ2ασ,0.
Based on the above analysis, the algorithm framework of sPADMM for solving (16) is summarized in Algorithm 1.
The following theorem provides the convergence analysis for Algorithm1. The proof is inspired by Theorem B.1 in [33], but no additional conditions are required in our proof.Algorithm 1: sPADMM for solving (16).
(1)
Initialization Let τ=1.618 and σ>0. Given a start point X0,V0,λ0,μ0,ξ0, for k=0,1,2,⋯,(2)
repeat(3)
Step 1. UpdateXk+1 according to (25)(4)
Step 2. UpdateVk+1 according to (28), (29) and (30)(5)
Step 3. Update the Lagrangian multipliers by (21)(6)
until Stopping criterion is satisfied.Theorem 5.
If sequenceXk,Vk,λk,μk,ξk is generated from Algorithm 1, then sequence Xk,Vk converges to an optimal solution for (16) and any accumulation point of Xk,Vk,λk,μk,ξk such that KKT conditions (19) hold.Proof.
We first justify that the conditions required in Theorem B.1 in [33] automatically hold for (16).(a)
The existence of a solution for (16): note that the objective function FX,Y,Z,W tends to infinity as X,Y,Z,WF⟶∞. Therefore, the level set(31)ℒ0≔X,Y,Z,WFX,Y,Z,W<F0.is bounded by Theorem 4.11 in [35]. Hence, we further know that the minimizer of (16) can be attained by Theorem 4.10 in [35].(b)
The constraints of (16) can be rewritten as A∗X−ℬ∗V=0, where(32)A∗X=JIIXandℬ∗V=V.
Hence, we find thatAA∗=JTJ+2I>0, so σP1+AA∗ is positive definite. Moreover, σP2+ℬℬ∗ is positive definite because ℬ∗ is the identity operator.
Based on these ingredients, the remainder of the proofs can be obtained easily (see Appendix B of [33] for details).


## 4.1. Optimality Conditions

By ignoring the terms withωℓ=0, model (2) can be rewritten as(14)minX∈ℝn×p12X−AF2+γ1∑ℓ∈ΕωℓXℓ1⋅−Xℓ2⋅2+γ2∑j=1p1−αujX⋅j2+αX⋅j1,where ℰ=ℓ=ℓ1,ℓ2ωℓ>0. For the convenience of mathematical expression, we define the linear operator J:ℝn×p⟶ℝℰ×p and its adjoint operator J∗:ℝℰ|×p⟶ℝn×p, respectively, by(15)JX=Xℓ1⋅−Xℓ2⋅ℓ1,ℓ2∈Ε=JXandJ∗Y=JTY,where Jℓ⋅=eℓ1−eℓ2T. Furthermore, model (13) can be reformulated as the following linearly constrained convex optimization problem with a separable structure:(16)minX,Y,Z,WFX,Y,Z,W≔12X−AF2+Q1Y+Q2Z+Q3W,s.t.JX−Y=0,X−Z=0,X−W=0,where Q1Y=γ1∑ℓ∈ℰωℓYℓ⋅2, Q2Z=γ21−α∑j=1pujZ⋅j2, and Q3W=γ2α∑j=1pX⋅j1. The Lagrangian function associated with (16) is defined by(17)LX,Y,Z,W;λ,μ,ξ=12X−AF2+Q1Y+Q2Z+Q3W+λ,JX−Y+μ,X−Z+ξ,X−W,where λ∈ℝℰ×p, μ∈ℝn×p, and ξ∈ℝn×p are Lagrangian multipliers. For a given parameter σ>0, the augmented Lagrangian function for (16) is defined by(18)LσX,Y,Z,W;λ,μ,ξ=LX,Y,Z,W;λ,μ,ξ+σ2JX−YF2+σ2X−ZF2+σ2X−WF2.The Karush–Kuhn–Tucker (KKT) conditions for (16) are given by(19)X−A+JTλ+μ+ξ=0,Y−ProxQ1Y+λ=0,Z−ProxQ2Z+μ=0,W−ProxQ3W+ξ=0,JX−Y=0,X−Z=0,X−W=0.Letℱ:=X,Y,Z,WJX−Y=0,X−Z=0,X−W=0 be a feasible set of (16). Obviously, ℱ≠∅ and FX,Y,Z,W≥0. Then, we know from Corollaries 28.2.2 and 28.3.1 in [29] that X^,Y^,Z^,W^ is an optimal solution for (16) if and only if there exists a Lagrangian multiplier λ^,μ^,ξ^ such that X^,Y^,Z^,W^,λ^,μ^,ξ^ satisfies the KKT conditions (19). Hence, we could design a stopping criterion based on the KKT conditions in terms of guaranteeing the optimality of a point generated by the algorithm proposed in Section 4.2.


## 4.2. Algorithm and Convergence Analysis

Although the optimization problem (16) is a convex optimization problem with linear constraints, it is difficult to minimize all four variables at the same time. Hence, we divide all four variables into two groups and update the two groups of variables alternately. According to the separable structure, we employ the semiproximal alternating direction method of multipliers (sPADMM) to solve (16) and view X as one block and V=Y,Z,W as another.Next, we discuss the iteration schemes of sPADMM for solving the optimization problem (16). Given the kth iteration Xk,Vk,λk,μk,ξk and two self-adjoint positive semidefinite operators P1 and P2, first update variables X and V as follows:(20)Xk+1=argminXLσX,Vk;λk,μk,ξk+σ2X−XkP12,Vk+1=argminVLσXk+1,V;λk,μk,ξk+σ2V−VkP22.Then, update Lagrangian multipliers:(21)λk+1=λk+τσJXk+1−Yk+1,μk+1=μk+τσXk+1−Zk+1,ξk+1=ξk+τσXk+1−Wk+1.Here,τ is step length and X−XkP12:=∑j=1pX⋅j−X⋅jkTP1X⋅j−X⋅jk. It is well known that the sPADMM reduces to the classical ADMM introduced by [31, 32] when P1=0 and P2=0. A comprehensive review of sPADMM can be found in Appendix B of [33].Remark 2.
The choice of the two self-adjoint positive semidefinite linear operatorsP1 and P2 depends greatly on the problem. The general principle is that both P1 and P2 should be as small as possible, while the corresponding subproblems are still relatively easy to solve.
Each block of (20) involves solving a convex optimization problem. Next, we show the closed-form solutions of each block. First, we deal with the first block in (20):(22)Xk+1=argminXLσX,Vk;λk,μk,ξk+σ2X−XkP12=argminX12X−AF2+σ2JX−Yk+λkσF2+σ2X−Zk+μkσF2+σ2X−Wk+ξkσF2+σ2X−XkP12.
TheX-subproblem is a smooth, strongly convex optimization that has a unique global minimizer Xk+1. Finding the minimizer Xk+1 is equivalent to finding the solution of the following linear system:(23)MX=A+σJTYk+Zk+Wk+P1Xk−JTλk−μk−ξk,where M=σJTJ+P1+1+2σIn. To solve this linear system efficiently, we design an appropriate P1. In our experiment, let P1=nIn−1n1nT−JTJ, which is a positive semidefinite matrix. Then, M=1+2σ+nσIn−σ1n1nT, and(24)M−1=11+2σ+nσIn+σ1+2σ1n1nT,by the Sherman–Morrison–Woodbury formula [34]. Hence, the closed-form solution for this linear system can be obtained as(25)Xk+1=M−1A+σJTYk+Zk+Wk+P1Xk−JTλk−μk−ξk.
For the second block in (20), we let the linear operator P2=0 since the linear operator of V in (16) is an identity operator that ensures the well-defined V-subproblem:(26)Vk+1=argminVLσXk+1,V;λk,μk,ξk=argminY,Z,WQ1Y+Q2Z+Q3W+σ2JXk+1+λkσ−YF2+σ2Xk+1+μkσ−ZF2+σ2Xk+1+ξkσ−WF2,which can be divided into three parts:(27)Yk+1=argminYQ1Y+σ2JXk+1+λkσ−YF2,Zk+1=argminZQ2Z+σ2Xk+1+μkσ−ZF2,Wk+1=argminWQ3W+σ2Xk+1+ξkσ−WF2.
Together with the definition of proximal mapping and Moreau decomposition, we show the closed-form solutions of the above subproblems. For anyℓ∈ℰ, the Y-subproblem has a closed-form solution that is given by(28)Yℓ⋅k+1=max1−γ1ωℓσJℓ⋅Xk+1+λℓ⋅k/σ2,0Jℓ⋅Xk+1+λℓ⋅kσ.
Similarly, for anyj=1,2,…,p, the solution of the Z-subproblem is given by(29)Z⋅jk+1=max1−γ2uj1−ασX⋅jk+1+μ⋅jk/σ2,0X⋅jk+1+μ⋅jkσ.
The last subproblem is a soft threshold operator for each column, that is, for anyj=1,2,…,p,(30)W⋅jk+1=signX⋅jk+1+ξ⋅jkσ°maxX⋅jk+1+ξ⋅jkσ−γ2ασ,0.
Based on the above analysis, the algorithm framework of sPADMM for solving (16) is summarized in Algorithm 1.
The following theorem provides the convergence analysis for Algorithm1. The proof is inspired by Theorem B.1 in [33], but no additional conditions are required in our proof.Algorithm 1: sPADMM for solving (16).
(1)
Initialization Let τ=1.618 and σ>0. Given a start point X0,V0,λ0,μ0,ξ0, for k=0,1,2,⋯,(2)
repeat(3)
Step 1. UpdateXk+1 according to (25)(4)
Step 2. UpdateVk+1 according to (28), (29) and (30)(5)
Step 3. Update the Lagrangian multipliers by (21)(6)
until Stopping criterion is satisfied.Theorem 5.
If sequenceXk,Vk,λk,μk,ξk is generated from Algorithm 1, then sequence Xk,Vk converges to an optimal solution for (16) and any accumulation point of Xk,Vk,λk,μk,ξk such that KKT conditions (19) hold.Proof.
We first justify that the conditions required in Theorem B.1 in [33] automatically hold for (16).(a)
The existence of a solution for (16): note that the objective function FX,Y,Z,W tends to infinity as X,Y,Z,WF⟶∞. Therefore, the level set(31)ℒ0≔X,Y,Z,WFX,Y,Z,W<F0.is bounded by Theorem 4.11 in [35]. Hence, we further know that the minimizer of (16) can be attained by Theorem 4.10 in [35].(b)
The constraints of (16) can be rewritten as A∗X−ℬ∗V=0, where(32)A∗X=JIIXandℬ∗V=V.
Hence, we find thatAA∗=JTJ+2I>0, so σP1+AA∗ is positive definite. Moreover, σP2+ℬℬ∗ is positive definite because ℬ∗ is the identity operator.
Based on these ingredients, the remainder of the proofs can be obtained easily (see Appendix B of [33] for details).


## 5. Numerical Results

The aim of this section is to illustrate the practical performance of sparse group lasso convex clustering (SGLCC), which performs well in theory. We conduct experiments both on synthetic and real datasets and compare the proposed method to K-means clustering, CC [10] and SCC [22]. The clustering results of K-means clustering can be obtained by built-in functions (kmeans) in MATLAB. All numerical experiments were performed in MATLAB R2017b on a personal computer with a 3.30 GHz processor and 4 GB of RAM.In our experiments, we stop sPADMM based on the relative KKT residual of (16), which is given by(33)ɛ=maxɛ1,ɛ2,ɛ3,ɛ4,ɛ5≤Tol,where(34)ɛ1=JX−YF+X−ZF+X−WF1+YF+ZF+WF,ɛ2=X−A+JTλ+μ+ξF1+AF,ɛ3=Y−ProxQ1Y+λF1+YF,ɛ4=Z−ProxQ2Z+μF1+ZF,ɛ5=W−ProxQ3W+ξF1+WF,and Tol is a given tolerance. If the requirement ɛ≤Tol is not satisfied after a maximum number of iterations (maxiter), we also terminate the algorithm. In our experiments, we set Tol=10−3 and maxiter=10000. Moreover, we use the same method to select the weights ωℓ and the adaptive weights uj in all experiments. For the weights ωℓ, we select the most popular and practical method proposed by [10]:(35)ωℓ=exp−ϕAℓ1⋅−Aℓ2⋅22,ℓ2is amongℓ1′sm−nearest neighbors,0,otherwise,where ϕ>0 and m<n. We set ϕ=0.5 and m=5 in all experiments. Inspired by [36], we let the adaptive weights uj=1/X^⋅j02, where X^⋅j0 is the estimator of X⋅j in (1). As can be seen from models (1) and (2), the tuning parameters γ1 and γ2,α control the number of clusters and the number of informative features, respectively. In our experiments, the principle for selecting γ1 and γ2,α is to make the results match the true number of clusters and the true number of informative features as much as possible, and the Rand index as large as possible. This principle is also used in the literature [21] to determine the tuning parameters.


### 5.1. Synthetic Data

In this section, we evaluate the performance of SGLCC on synthetic datasets by comparing the results with those obtained with K-means clustering, CC and SCC. We consider three synthetic datasets. Each synthetic dataset consists ofn=200 observations with either K=2 or 4 clusters. The number of features is p and ranges from 2000 to 5000, where the percentage of informative features is 2%. The cluster label Lii=1,2,…,n is uniformly sampled from 1,2,…,K. Without loss of generality, we assume that the first s2%p features are informative and the remaining features are noninformative. The first s informative features are generated from an s-dimensional normal distribution with mean x¯KLi and covariance Σs, whose i,j entry is Σsi,j=ρi−j,1≤i,j≤s, and the noninformative features are generated from a p−s-dimensional standard normal distribution. Here,(36)x¯KLi=1sILi=1−1sILi=2,K=2,1s/2,1s/2ILi=1+1s/2,−1s/2ILi=2+−1s/2,1s/2ILi=3+−1s/2,−1s/2ILi=4,K=4,where I⋅ is a function whose value is one if the event happens and zero otherwise. Let us illustrate this definition with an example. Specifically, if K=2, the first s informative features of the observations with cluster labels 1 and 2 are generated from s-dimensional normal distributions with mean x¯21=1s and x¯22=−1s, respectively. In this part, we consider the following three synthetic datasets. Case I: K=2,Σsi,j=0.5i−j; Case II: K=4,Σsi,j=0.5i−j; Case III: K=4,Σsi,j=0.8i−j.We repeat each experiment 50 times and evaluate the performance through four common criteria for cluster analysis. The Rand index [37] and Fowlkes–Mallows index [38] are two measures of the similarity between the estimated clustering result and the underlying true clustering assignment. Their values range from 0 to 1, and a larger value indicates better performance. For simplicity, we shall abbreviate the Rand index and Fowlkes–Mallows index as RI and FMI, respectively. Moreover, the false negative rate (FNR) and false positive rate (FPR) are reported to evaluate the performance of feature screening. All experimental results for the synthetic datasets are summarized in Tables 1–3.Table 1
Results for Case I. Empirical mean and standard deviation (SD) of the RI, FMI, FNR, and FPR based on 50 repetitions.
Case I:K=2,Σsi,j=0.5i−jpMethodsRIFMIFNRFPRMeanSDMeanSDMeanSDMeanSD2000K-means0.9700.1200.9800.0790.0000.0001.0000.000CC0.8160.0770.7890.0960.0000.0001.0000.000SCC1.0000.0001.0000.0000.0000.0000.2230.010SGLCC1.0000.0001.0000.0000.0000.0000.0000.0003000K-means0.9800.0990.9860.0690.0000.0001.0000.000CC0.8040.0760.7730.0960.0000.0001.0000.000SCC1.0000.0001.0000.0000.0000.0000.2210.008SGLCC1.0000.0001.0000.0000.0000.0000.0000.0004000K-means1.0000.0001.0000.0000.0000.0001.0000.000CC0.7840.0830.7470.1040.0000.0001.0000.000SCC0.9990.0040.9990.0050.0000.0000.2230.007SGLCC1.0000.0001.0000.0000.0000.0000.0000.0005000K-means0.9900.0710.9940.0420.0000.0001.0000.000CC0.7620.0710.7170.0930.0000.0001.0000.000SCC0.9980.0060.9980.0060.0000.0000.2240.008SGLCC1.0000.0001.0000.0000.0000.0010.0000.000Table 2
Results for Case II. Empirical mean and standard deviation (SD) of the RI, FMI, FNR, and FPR based on 50 repetitions.
Case II:K=4,Σsi,j=0.5i−jpMethodsRIFMIFNRFPRMeanSDMeanSDMeanSDMeanSD2000K-means0.9070.0730.8480.1170.0000.0001.0000.000CC0.9700.0450.9420.0850.0000.0001.0000.000SCC0.9980.0020.9950.0040.0000.0000.5500.127SGLCC0.9950.0210.9910.0340.0150.0450.0000.0003000K-means0.8870.0870.8150.1380.0000.0001.0000.000CC0.9280.0290.8410.0680.0000.0001.0000.000SCC0.9970.0150.9930.0320.0000.0000.6810.071SGLCC1.0000.0001.0000.0000.0000.0030.0000.0004000K-means0.8920.0760.8230.1220.0000.0001.0000.000CC0.9250.0300.8340.0690.0000.0001.0000.000SCC0.9980.0060.9960.0120.0000.0000.6850.007SGLCC1.0000.0001.0000.0000.0000.0010.0000.0005000K-means0.8850.0770.8100.1230.0000.0001.0000.000CC0.9110.0280.8000.0680.0000.0001.0000.000SCC0.9860.0330.9680.0770.0000.0000.7110.085SGLCC1.0000.0001.0000.0000.0000.0000.0000.000Table 3
Results for Case III. Empirical mean and standard deviation (SD) of the RI, FMI, FNR, and FPR based on 50 repetitions.
Case III:K=4,Σsi,j=0.8i−jpMethodsRIFMIFNRFPRMeanSDMeanSDMeanSDMeanSD2000K-means0.9030.0810.8370.1340.0000.0001.0000.000CC0.9660.0630.9500.0850.0000.0001.0000.000SCC0.9880.0270.9770.0450.0150.0750.5060.115SGLCC0.9930.0060.9870.0120.0290.1150.1280.0703000K-means0.9140.0850.8600.1360.0000.0001.0000.000CC0.9470.0710.9240.0960.0000.0001.0000.000SCC0.9890.0150.9780.0310.0070.0350.7750.125SGLCC0.9930.0190.9860.0300.0290.1140.0890.0414000K-means0.8970.0800.8340.1220.0000.0001.0000.000CC0.9150.0970.7930.2530.0000.0001.0000.000SCC0.9330.0800.8250.2280.0250.1010.8230.231SGLCC0.9830.0600.9790.0700.0300.1190.0770.0285000K-means0.8920.0850.8250.1360.0000.0001.0000.000CC0.9000.0480.7640.1150.0000.0001.0000.000SCC0.9380.0580.8560.1370.0000.0000.7800.214SGLCC0.9990.0030.9990.0060.0100.0680.0030.001From Table1, several interesting observations can be made. (a) When the feature dimensions are high (p=2000,3000,4000,5000), CC does not perform well, and its clustering results worsen as p increases. (b) Because CC and K-means do not have the ability to eliminate noninformative features, their FNR and FPR are zero and one, respectively. (c) The clustering results of SGLCC and SCC are significantly better than those of CC when the feature dimensions are high because the penalty term is incorporated in convex clustering (1). (d) SGLCC has more stable performance than K-means and CC because it has the smallest standard deviation (SD). (e) SGLCC selects informative features with greater accuracy than the SCC, that is, with a lower FPR and almost the same FNR. The phenomenon mentioned above can also be observed from the results for Case II, see Table 2 for details. It is worth noting that SCC and SGLCC have similar performance for Case II. This reflects the ability of SCC to conduct feature selection, especially for data with group sparsity and dimensions that are not very large. However, compared with that of SGLCC, the performance of SCC becomes unstable as the number of dimensions increases. The reason behind this difference is that SGLCC allows a more nuanced selection of informative features.The experimental results for Case III are reported in Figure1 and Table 3. As we can see from Figure 1, SGLCC is significantly better than SCC when there are highly correlated variables in the dataset, and the improvement becomes increasingly evident as p increases. The main reason behind this may be the form of the sparse group lasso penalty, which, much like the elastic net, is good at handling the collinearity problem.Figure 1
Results for Case III:K=4,Σsi,j=0.8i−j.


### 5.2. Real Data

In this section, five public datasets are used to validate the performance of our proposed SGLCC method. These datasets have been broadly applied to numerous studies in the literature, which is the primary reason why we selected them. All real datasets are summarized in Table4. Brain A contains 42 observations with 5597 genes of five clusters that contain 10, 10, 10, 4, and 8 observations. SRBCT contains 63 observations of four clusters, and each observation has 2308 genes. The four clusters contain 23, 8, 12, and 20 observations. Prostate contains 102 observations with 5597 genes of two clusters that contain 50 and 52. Colon contains 62 observations of two clusters that contain 22 and 40 observations, each of which has 2000 genes. ARB contains 590 observations of eight clusters, and each observation has 8266 variables. The eight clusters contain 46, 81, 162, 189, 31, 12, 50, and 19 observations. We used the preprocessing versions of the first four datasets based on [39]. In our experiments, each feature from all of the datasets is centered. The Rand index of K-means, hierarchical clustering (Hclust), spectral clustering (Sclust), CC, SCC, and SGLCC for all five real datasets are reported in Table 5, and the best results are highlighted in bold. Here, the clustering results of Hclust and Sclust are obtained through two built-in functions in MATLAB, namely, clusterdata and spectralcluster.Table 4
Details of the real datasets.
DatasetsFeaturesObservationsClustersPoints of per clusterBrain a [40]559742510, 10, 10, 4, 8SRBCT [41]230863423, 8, 12, 20Prostate [42]6033102250, 52Colon [43]200062222, 40ARB8266590846, 81, 162, 189, 31, 12, 50, 19Table 5
Results for the real datasets.
DatasetsMethodsRand indexBrain AK-means0.798Hclust0.317Sclust0.506CC0.577SCC0.793SGLCC0.890ProstateK-means0.505Hclust0.496Sclust0.497CC0.511SCC0.518SGLCC0.530ARBK-means0.539Hclust0.233Sclust0.635CC0.485SCC0.562SGLCC0.675SRBCTK-means0.593Hclust0.308Sclust0.525CC0.595SCC0.634SGLCC0.770ColonK-means0.556Hclust0.527Sclust0.511CC0.492SCC0.596SGLCC0.640As we can see from Table5, the performance of CC is unsatisfactory and worse than that of K-means for Brain A, Colon, and ARB. A similar phenomenon can also be observed in [15, 22]. SGLCC is significantly superior to K-means, Hclust, Sclust, and CC in term of clustering performance, which implies that feature screening is an indispensable part of high-dimensional data analysis. Moreover, SGLCC is better than SCC because a more reasonable penalty term is used. In summary, our proposed SGLCC method has the largest Rand index among the six methods.


## 5.1. Synthetic Data

In this section, we evaluate the performance of SGLCC on synthetic datasets by comparing the results with those obtained with K-means clustering, CC and SCC. We consider three synthetic datasets. Each synthetic dataset consists ofn=200 observations with either K=2 or 4 clusters. The number of features is p and ranges from 2000 to 5000, where the percentage of informative features is 2%. The cluster label Lii=1,2,…,n is uniformly sampled from 1,2,…,K. Without loss of generality, we assume that the first s2%p features are informative and the remaining features are noninformative. The first s informative features are generated from an s-dimensional normal distribution with mean x¯KLi and covariance Σs, whose i,j entry is Σsi,j=ρi−j,1≤i,j≤s, and the noninformative features are generated from a p−s-dimensional standard normal distribution. Here,(36)x¯KLi=1sILi=1−1sILi=2,K=2,1s/2,1s/2ILi=1+1s/2,−1s/2ILi=2+−1s/2,1s/2ILi=3+−1s/2,−1s/2ILi=4,K=4,where I⋅ is a function whose value is one if the event happens and zero otherwise. Let us illustrate this definition with an example. Specifically, if K=2, the first s informative features of the observations with cluster labels 1 and 2 are generated from s-dimensional normal distributions with mean x¯21=1s and x¯22=−1s, respectively. In this part, we consider the following three synthetic datasets. Case I: K=2,Σsi,j=0.5i−j; Case II: K=4,Σsi,j=0.5i−j; Case III: K=4,Σsi,j=0.8i−j.We repeat each experiment 50 times and evaluate the performance through four common criteria for cluster analysis. The Rand index [37] and Fowlkes–Mallows index [38] are two measures of the similarity between the estimated clustering result and the underlying true clustering assignment. Their values range from 0 to 1, and a larger value indicates better performance. For simplicity, we shall abbreviate the Rand index and Fowlkes–Mallows index as RI and FMI, respectively. Moreover, the false negative rate (FNR) and false positive rate (FPR) are reported to evaluate the performance of feature screening. All experimental results for the synthetic datasets are summarized in Tables 1–3.Table 1
Results for Case I. Empirical mean and standard deviation (SD) of the RI, FMI, FNR, and FPR based on 50 repetitions.
Case I:K=2,Σsi,j=0.5i−jpMethodsRIFMIFNRFPRMeanSDMeanSDMeanSDMeanSD2000K-means0.9700.1200.9800.0790.0000.0001.0000.000CC0.8160.0770.7890.0960.0000.0001.0000.000SCC1.0000.0001.0000.0000.0000.0000.2230.010SGLCC1.0000.0001.0000.0000.0000.0000.0000.0003000K-means0.9800.0990.9860.0690.0000.0001.0000.000CC0.8040.0760.7730.0960.0000.0001.0000.000SCC1.0000.0001.0000.0000.0000.0000.2210.008SGLCC1.0000.0001.0000.0000.0000.0000.0000.0004000K-means1.0000.0001.0000.0000.0000.0001.0000.000CC0.7840.0830.7470.1040.0000.0001.0000.000SCC0.9990.0040.9990.0050.0000.0000.2230.007SGLCC1.0000.0001.0000.0000.0000.0000.0000.0005000K-means0.9900.0710.9940.0420.0000.0001.0000.000CC0.7620.0710.7170.0930.0000.0001.0000.000SCC0.9980.0060.9980.0060.0000.0000.2240.008SGLCC1.0000.0001.0000.0000.0000.0010.0000.000Table 2
Results for Case II. Empirical mean and standard deviation (SD) of the RI, FMI, FNR, and FPR based on 50 repetitions.
Case II:K=4,Σsi,j=0.5i−jpMethodsRIFMIFNRFPRMeanSDMeanSDMeanSDMeanSD2000K-means0.9070.0730.8480.1170.0000.0001.0000.000CC0.9700.0450.9420.0850.0000.0001.0000.000SCC0.9980.0020.9950.0040.0000.0000.5500.127SGLCC0.9950.0210.9910.0340.0150.0450.0000.0003000K-means0.8870.0870.8150.1380.0000.0001.0000.000CC0.9280.0290.8410.0680.0000.0001.0000.000SCC0.9970.0150.9930.0320.0000.0000.6810.071SGLCC1.0000.0001.0000.0000.0000.0030.0000.0004000K-means0.8920.0760.8230.1220.0000.0001.0000.000CC0.9250.0300.8340.0690.0000.0001.0000.000SCC0.9980.0060.9960.0120.0000.0000.6850.007SGLCC1.0000.0001.0000.0000.0000.0010.0000.0005000K-means0.8850.0770.8100.1230.0000.0001.0000.000CC0.9110.0280.8000.0680.0000.0001.0000.000SCC0.9860.0330.9680.0770.0000.0000.7110.085SGLCC1.0000.0001.0000.0000.0000.0000.0000.000Table 3
Results for Case III. Empirical mean and standard deviation (SD) of the RI, FMI, FNR, and FPR based on 50 repetitions.
Case III:K=4,Σsi,j=0.8i−jpMethodsRIFMIFNRFPRMeanSDMeanSDMeanSDMeanSD2000K-means0.9030.0810.8370.1340.0000.0001.0000.000CC0.9660.0630.9500.0850.0000.0001.0000.000SCC0.9880.0270.9770.0450.0150.0750.5060.115SGLCC0.9930.0060.9870.0120.0290.1150.1280.0703000K-means0.9140.0850.8600.1360.0000.0001.0000.000CC0.9470.0710.9240.0960.0000.0001.0000.000SCC0.9890.0150.9780.0310.0070.0350.7750.125SGLCC0.9930.0190.9860.0300.0290.1140.0890.0414000K-means0.8970.0800.8340.1220.0000.0001.0000.000CC0.9150.0970.7930.2530.0000.0001.0000.000SCC0.9330.0800.8250.2280.0250.1010.8230.231SGLCC0.9830.0600.9790.0700.0300.1190.0770.0285000K-means0.8920.0850.8250.1360.0000.0001.0000.000CC0.9000.0480.7640.1150.0000.0001.0000.000SCC0.9380.0580.8560.1370.0000.0000.7800.214SGLCC0.9990.0030.9990.0060.0100.0680.0030.001From Table1, several interesting observations can be made. (a) When the feature dimensions are high (p=2000,3000,4000,5000), CC does not perform well, and its clustering results worsen as p increases. (b) Because CC and K-means do not have the ability to eliminate noninformative features, their FNR and FPR are zero and one, respectively. (c) The clustering results of SGLCC and SCC are significantly better than those of CC when the feature dimensions are high because the penalty term is incorporated in convex clustering (1). (d) SGLCC has more stable performance than K-means and CC because it has the smallest standard deviation (SD). (e) SGLCC selects informative features with greater accuracy than the SCC, that is, with a lower FPR and almost the same FNR. The phenomenon mentioned above can also be observed from the results for Case II, see Table 2 for details. It is worth noting that SCC and SGLCC have similar performance for Case II. This reflects the ability of SCC to conduct feature selection, especially for data with group sparsity and dimensions that are not very large. However, compared with that of SGLCC, the performance of SCC becomes unstable as the number of dimensions increases. The reason behind this difference is that SGLCC allows a more nuanced selection of informative features.The experimental results for Case III are reported in Figure1 and Table 3. As we can see from Figure 1, SGLCC is significantly better than SCC when there are highly correlated variables in the dataset, and the improvement becomes increasingly evident as p increases. The main reason behind this may be the form of the sparse group lasso penalty, which, much like the elastic net, is good at handling the collinearity problem.Figure 1
Results for Case III:K=4,Σsi,j=0.8i−j.


## 5.2. Real Data

In this section, five public datasets are used to validate the performance of our proposed SGLCC method. These datasets have been broadly applied to numerous studies in the literature, which is the primary reason why we selected them. All real datasets are summarized in Table4. Brain A contains 42 observations with 5597 genes of five clusters that contain 10, 10, 10, 4, and 8 observations. SRBCT contains 63 observations of four clusters, and each observation has 2308 genes. The four clusters contain 23, 8, 12, and 20 observations. Prostate contains 102 observations with 5597 genes of two clusters that contain 50 and 52. Colon contains 62 observations of two clusters that contain 22 and 40 observations, each of which has 2000 genes. ARB contains 590 observations of eight clusters, and each observation has 8266 variables. The eight clusters contain 46, 81, 162, 189, 31, 12, 50, and 19 observations. We used the preprocessing versions of the first four datasets based on [39]. In our experiments, each feature from all of the datasets is centered. The Rand index of K-means, hierarchical clustering (Hclust), spectral clustering (Sclust), CC, SCC, and SGLCC for all five real datasets are reported in Table 5, and the best results are highlighted in bold. Here, the clustering results of Hclust and Sclust are obtained through two built-in functions in MATLAB, namely, clusterdata and spectralcluster.Table 4
Details of the real datasets.
DatasetsFeaturesObservationsClustersPoints of per clusterBrain a [40]559742510, 10, 10, 4, 8SRBCT [41]230863423, 8, 12, 20Prostate [42]6033102250, 52Colon [43]200062222, 40ARB8266590846, 81, 162, 189, 31, 12, 50, 19Table 5
Results for the real datasets.
DatasetsMethodsRand indexBrain AK-means0.798Hclust0.317Sclust0.506CC0.577SCC0.793SGLCC0.890ProstateK-means0.505Hclust0.496Sclust0.497CC0.511SCC0.518SGLCC0.530ARBK-means0.539Hclust0.233Sclust0.635CC0.485SCC0.562SGLCC0.675SRBCTK-means0.593Hclust0.308Sclust0.525CC0.595SCC0.634SGLCC0.770ColonK-means0.556Hclust0.527Sclust0.511CC0.492SCC0.596SGLCC0.640As we can see from Table5, the performance of CC is unsatisfactory and worse than that of K-means for Brain A, Colon, and ARB. A similar phenomenon can also be observed in [15, 22]. SGLCC is significantly superior to K-means, Hclust, Sclust, and CC in term of clustering performance, which implies that feature screening is an indispensable part of high-dimensional data analysis. Moreover, SGLCC is better than SCC because a more reasonable penalty term is used. In summary, our proposed SGLCC method has the largest Rand index among the six methods.


## 6. Conclusion

In this paper, we propose the Sparse Group Lasso Convex Clustering method for high-dimensional datasets, which can determine the number of clusters by learning from the data. SGLCC not only divides the observation set but also eliminates noninformative features. For the proposed clustering method, we provide the theoretical finite sample bounds of the prediction error and feature screening consistency. Moreover, an efficient sPADMM is designed to implement our proposed estimator, and each subproblem yields closed-form solutions. Convergence analysis of the sPADMM is established without any additional conditions. The experimental results illustrate that SGLCC provides superior clustering performance and feature selection abilities than other compared clustering methods. In particular, our method is also very effective for real data applications.



---

*Source: 9216351-2020-09-21.xml*