\section{Introduction}

Cluster analysis is an unsupervised learning method and aims to assign
observations into a number of clusters such that observations in the
same group are similar to each other. Traditional clustering methods
such as k-means clustering, hierarchical clustering, and Gaussian mixture
models take a greedy approach and suffer from instabilities due to
their non-convex optimization formulations.

To overcome the instability issues of these traditional clustering methods, a new clustering algorithm, \textit{Convex Clustering},
has been recently proposed \citep{Pelckmans.ea:2005,Lindsten.ea:2004,Hocking2011}. Let
$\mathbf{X}\in\mathbb{R}^{n\times p}$ be a data matrix with $n$
observations $X_{i\cdot}$, $i=1,\cdots,n$, and $p$ features. Convex
clustering for these $n$ observations solves the following minimization
problem:
\begin{equation}
\min_{\mathbf{A}\in \mathbb{R}^{n\times p}}\frac{1}{2}\sum_{i=1}^{n}||X_{i\cdot}-A_{i\cdot}||_{2}^{2}+\gamma \sum_{i_{1}<i_{2}}||A_{i_{1}\cdot}-A_{i_{2}\cdot}||_{q},\label{eq:obj-cc}
\end{equation}
where $A_{i\cdot}$ is the $i$-th row of $\mathbf{A}$ and $\|\cdot\|_{q}$ is the $L_{q}$-norm of a vector with $q \in \{1,2,\infty\}$. Note that both k-means clustering and hierarchical
clustering consider $L_{0}$-norm in the second term, which leads to
a non-convex optimization problem \citep{Hocking2011,Tan2015}. Therefore, convex
clustering can be viewed as a convex relaxation of k-means clustering and hierarchical
clustering, and the convex relaxation ensures that it achieves a unique
global minimizer.

Due to the fused-lasso penalty \citep{Tibshirani.ea:2005} in the second term of (\ref{eq:obj-cc}),
the above formulation encourages that some of the rows of the solution $\widehat{{\bf A}}$
are identical. If $\widehat{A}_{i_{1}\cdot}=\widehat{A}_{i_{2}\cdot}$,
then observation $i_{1}$ and observation $i_{2}$ are said to belong
to the same cluster. The tuning parameter $\gamma$ in (\ref{eq:obj-cc})
controls the number of unique rows of $\widehat{{\bf A}}$, that is,
the number of estimated clusters. When $\gamma=0$, $\widehat{{\bf A}}={{\bf X}}$,
and therefore each observation by itself is a cluster. As $\gamma$
increases, some of the rows of $\widehat{{\bf A}}$ become identical,
which demonstrates a fusion process. For sufficiently large $\gamma$,
all the rows of $\widehat{{\bf A}}$ will be identical, implying that
all the observations are estimated to belong to a single cluster.
Compared to traditional non-convex clustering methods, the
solution $\widehat{{\bf A}}$ from convex clustering is unique for
each given $\gamma$ since the objective function in (\ref{eq:obj-cc})
is strictly convex.

In recent years, the computational and statistical properties of convex
clustering have been investigated. In particular, \citet{Zhu2014} provided conditions for convex clustering
to recover the true clusters, \citet{Chi2015} proposed efficient
and scalable implementations for convex clustering, and \citet{Tan2015}
studied several statistical properties of convex clustering. While convex clustering enjoys nice theoretical properties and is
computationally efficient, its performance can be severely deteriorated
when clustering high-dimensional data where the number of features
becomes large and many of them may contain no information about the
clustering structure. Our extensive experimental studies demonstrate
that in high-dimensional scenarios the performance of convex clustering
is unsatisfactory when the uninformative features are included in
the clustering. To overcome such a difficulty, a more appropriate convex
clustering algorithm that can simultaneously perform cluster analysis
and select informative variables is in demand.

In this article, we introduce a new clustering method, \textit{Sparse
Convex Clustering}, to incorporate the sparsity into convex clustering of high dimensional data. The key idea is to formulate convex clustering
in a form of regularization, with an adaptive group-lasso penalty
term on cluster centers to encourage the sparsity. Despite its simplicity,
this regularization operator demands more challenging computational
and statistical analysis than those in original convex clustering.
In particular, computationally, we need to reformulate the sparse convex
clustering into a few sub-optimization problems and then solve each
individual one via a pseudo regression formulation. To prove an unbiased
estimator for the degrees of freedom of the proposed sparse convex clustering
method, we need to carefully quantify the impact of variable selection
due to the group lasso penalty. Moreover, we provide a non-asymptotic analysis for the prediction error of our sparse convex clustering estimator. Under a high-dimensional scenario where the dimension diverges with the sample size, our estimator is further shown to be consistent in variable selection.
Note that our method is not only theoretical sound, but also practically promising. The superior performance
of our procedure is demonstrated in extensive simulated examples and
a real application of hand movement clustering.

We demonstrate the superior performance of the proposed method using a dataset generated from the fourth simulation setting in Section 5. In this dataset, there are 60 subjects from 4 clusters and 500 features, among which the first 20 features are informative. Figure \ref{heatmap} compares the performance of convex clustering (indicated as AMA on the left panel) and sparse convex clustering (indicated as S-AMA on the right panel), by visualizing the regularized feature matrix $\widehat{{\bf A}}$ estimated in each of the two methods. The heap maps show that sparse convex clustering screens out those uninformative features and therefore improves the clustering performance.

\begin{figure}[!htb]
\protect\caption{The heat maps of $p\times n$ matrix $\widehat{{\bf A}}^{\rm T}$, estimated from Convex clustering and Sparse Convex clustering, respectively. Data are generated using Simulation Setting 4 where there are 60 subjects in 4 clusters and 500 features, among which the first 20 features are informative.}
\centering \includegraphics[scale=0.4]{Figure1_colored.pdf}\label{heatmap}
\end{figure}


\subsection{Related Work}

A related paper on convex clustering is its efficient implementations proposed by \citet{Chi2015} and its extension to convex biclustering has been developed by \cite{Chi2016convex}.  Two efficient algorithms ADMM and AMA are introduced while they are mainly designed for clustering low-dimensional data. In order
to address high dimensionality, one key ingredient of our sparse convex
clustering method is a new regularization penalty built upon their
ADMM and AMA algorithms to encourage the sparsity structure of the
clustering centers. As will be shown in experimental studies, such
regularization step is able to significantly improve the clustering
accuracy in high-dimensional clustering problems.

Another line of research focuses on simultaneous clustering and feature
selection. Some approaches are model-based clustering methods, such
as \citet{Raftery.Dean:2006}, \citet{Pan.Shen:2007}, \citet{Wang.Zhu:2008}, \citet{Xie.ea:2008}, and \citet{Guo.ea:2010}. In contrast,
some approaches are model-free, such as \citet{Witten.Tibshirani:2010}, \citet{Sun2012}, and \citet{Wang.ea:2013}. One common building block of
these sparse clustering approaches is the usage of a lasso-type penalty
for feature selection. For example, \citet{Witten.Tibshirani:2010}
developed a unified framework for feature selection in clustering using the
lasso penalty \citep{tibshirani:1996}. \citet{Sun2012} proposed a sparse k-means using the group-lasso penalty \citep{Yuan2006}.
We refer readers to \citet{Alelyani.ea:2013} for a thorough overview. In spite of their good numeric performance, these
sparse clustering procedures still suffer from instabilities due to
the non-convex optimization formulations. To overcome it, our sparse
convex clustering solves a convex optimization problem and ensures
a unique global solution.

\subsection{Paper Organization}

The rest of the manuscript is organized as follows. Section \ref{sec:main}
introduces the sparse convex clustering as well as its two efficient algorithms. Section \ref{sec:thm}
studies its statistical properties and Section \ref{sec:practical}
discusses some practical issues in the proposed implementations.
Section \ref{sec:experiment} evaluates the superior numeric performance
of the proposed methods through extensive simulations and a real data
application. Section \ref{sec:summary} finishes this paper with a brief summary. Technical details are provided in Appendix or online supplementary.


\section{Sparse Convex Clustering}

\label{sec:main}

This section presents the main results. We propose our new method, sparse convex clustering in Section \ref{sec:model}, and then we develop two efficient algorithms to implement the method in Section \ref{sec:algorithm}

\subsection{Model}

\label{sec:model}

To allow an adaptive penalization, we consider a modification of convex
clustering (\ref{eq:obj-cc}),

\begin{equation}
\min_{\mathbf{A}\in \mathbb{R}^{n\times p}}\frac{1}{2}\sum_{i=1}^{n}||X_{i\cdot}-A_{i\cdot}||_{2}^{2}+\gamma\sum_{i_{1}<i_{2}}w_{i_{1},i_{2}}||A_{i_{1}\cdot}-A_{i_{2}\cdot}||_{q},\label{eq:obj-cc-w}
\end{equation}
where the weight $w_{i_{1},i_{2}}\ge0$. \citet{Hocking2011} considered a pairwise
affinity weight $w_{i_{1},i_{2}}=\exp(-\phi\|X_{i_{1}\cdot}-X_{i_{2}\cdot}\|_{2}^{2})$
and \citet{Chi2015} suggested $w_{i_{1},i_{2}}=\iota_{i_{1},i_{2}}^{m}\exp(-\phi\|X_{i_{1}\cdot}-X_{i_{2}\cdot}\|_{2}^{2})$,
where $\iota_{i_{1},i_{2}}^{m}$ is 1 if observation $i_{2}$ is among
$i_{1}$'s $m$ nearest neighbors or vice verse, and 0 otherwise.

To introduce a reformulation of (\ref{eq:obj-cc-w}), we write the
data matrix $\mathbf{X}$ in feature-level as column vector $\mathbf{X}=(\mathbf{x}_{1},\cdots,\mathbf{x}_{p})$,
where $\mathbf{x}_{j}=(X_{1j},\cdots,X_{nj})^{^{\rm T}}$, $j=1,\ldots,p$
and denote $\mathbf{A}$ in feature-level as column vector $\mathbf{A}=(\mathbf{a}_{1},\cdots,\mathbf{a}_{p})$.
Without loss of generality, we assume the feature vectors are centered,
i.e., $\sum_{i=1}^{n}X_{ij}=0$ for each $j=1,\ldots, p$. Simple algebra implies that (\ref{eq:obj-cc-w})
can be reformulated as
\begin{equation} \label{eq:obj-cc-w-new}
\min_{\mathbf{A}\in \mathbb{R}^{n\times p}}\frac{1}{2}\sum_{j=1}^{p}||\mathbf{x}_{j}-\mathbf{a}_{j}||_{2}^{2}+\gamma\sum_{l\in{\cal E}}w_{l}||A_{i_{1}\cdot}-A_{i_{2}\cdot}||_{q},
\end{equation}
where ${\cal E}=\{l=(i_{1},i_{2}):1\leq i_{1}<i_{2}\leq n\}$.

For a given $\gamma$, let $\widehat{\mathbf{A}}=(\widehat{A}_{1\cdot},\cdots,\widehat{A}_{n.})^{^{\rm T}}=(\widehat{\mathbf{a}}_{1},\ldots,\widehat{\mathbf{a}}_{p})$
be the solution to (\ref{eq:obj-cc-w-new}). The clustering structure
is implied by the observation-level estimates, $\widehat{A}_{i\cdot}$,
$i=1,\ldots,n$; that is, if $\widehat{A}_{i_{1}\cdot}=\widehat{A}_{i_{2}\cdot}$,
then observations $i_{1}$ and $i_{2}$ are estimated to belong to
the same cluster. The feature importance is implied by the feature-level
estimates, $\widehat{\mathbf{a}}_{j}$, $j=1,\cdots,p$; that is, if the
components of a feature-level estimate $\widehat{\mathbf{a}}_{j}$ are
identical, then the corresponding feature $j$ is not informative
for clustering. Remind that the feature vectors are centered, then
feature $j$ is not informative if and only if $\|\widehat{\mathbf{a}}_{j}\|_{2}^{2}=\sum_{i=1}^{n}\widehat{A}_{ij}^{2}=0$.

In high-dimensional clustering, it is desired to have a sparse solution
$\widehat{{\bf A}}$ with some of its column vectors being exact ${\bf 0}$'s.
Motivated by the importance of excluding non-informative features, we propose a new sparse convex clustering by incorporating
an adaptive group-lasso penalty \citep{Yuan2006,Wang.Leng:2008} into
the convex clustering objective function (\ref{eq:obj-cc-w-new}).
In particular, sparse convex clustering solves

\begin{equation}\label{eq:obj_constraint}
\min_{\mathbf{A}\in \mathbb{R}^{n\times p}}\frac{1}{2}\sum_{j=1}^{p}||\mathbf{x}_{j}-\mathbf{a}_{j}||_{2}^{2}+\gamma_{1}\sum_{l\in{\cal E}}w_{l}||A_{i_{1}\cdot}-A_{i_{2}\cdot}||_{q}+\gamma_{2}\sum_{j=1}^{p}u_{j}||\mathbf{a}_{j}||_{2},
\end{equation}
where tuning parameter $\gamma_{1}$ controls the cluster size and
tuning parameter $\gamma_{2}$ controls the number of informative
features. In the group-lasso penalty, the weight $u_{j}$ plays
an important role to adaptively penalize the features. Detailed discussions on practical choices of tuning parameters and weights can be found in Section \ref{sec:practical}.


\begin{remark}
In the objective function (\ref{eq:obj_constraint}) of sparse convex clustering, the second group-lasso-type penalty enforces the global sparsity condition; that is, the elements of each column vector ${\bf a}}  \def{\bf A}{{\bf A}_j$ would be all zero or all nonzero. Such penalty is considered for the feature selection purpose. This global sparsity condition can be relaxed in two directions. First, we can replace the second penalty, $\gamma_2\sum_{j=1}^{p}u_{j}||\mathbf{a}_{j}||_{2}$, by a lasso type of penalty, $\gamma_2\sum_{j=1}^{p}||\mathbf{a}_{j}||_{1}$. Second, we can also add another penalty, $\sum_{j=1}^{p}||\mathbf{a}_{j}||_{1}$, to the objective function \ref{eq:obj_constraint}. This results in a new penalty, $\gamma_2\sum_{j=1}^{p}u_{j}||\mathbf{a}_{j}||_{2}+\gamma_3\sum_{j=1}^{p}||\mathbf{a}_{j}||_{1}$, which is the so-called sparse-group-lasso penalty \citep{Friedman2010}.
\end{remark}


\subsection{Algorithms}

\label{sec:algorithm}

This subsection discusses two efficient optimization approaches to
solve the sparse convex clustering by adopting a similar computational
strategy used in \citet{Chi2015}. Our two approaches are based on
the alternating direction method of multipliers (ADMM) algorithm \citep{Boyd2011,Gabay1976,Glowinski1975}
and the alternating minimization algorithm (AMA) \citep{Tseng1991},
and are referred as sparse ADMM (S-ADMM) and sparse AMA (S-AMA), respectively.

To implement the S-ADMM and S-AMA algorithms, we rewrite the
convex clustering problem in formula (\ref{eq:obj_constraint}) as
\begin{eqnarray*}
\min_{{\bf A}\in \mathbb{R}^{n\times p}} && \frac{1}{2}\sum_{j=1}^{p}||\mathbf{x}_{j}-\mathbf{a}_{j}||_{2}^{2}+\gamma_{1}\sum_{l\in{\cal E}}w_{l}||\mathbf{v}_{l}||_{q}+\gamma_{2}\sum_{j=1}^{p}u_{j}||\mathbf{a}_{j}||_{2}, \\
{\rm s.t.} && A_{i_1\cdot}-A_{i_2\cdot}-{\bf v}}  \def{\bf V}{{\bf V}_l={\bf 0}.
\end{eqnarray*}
This is equivalent to minimize the following augmented Lagrangian
function,
\begin{eqnarray*}
\mathcal{L}_{\nu}(\mathbf{A},\mathbf{V},\mathbf{\Lambda}) & = & \frac{1}{2}\sum_{j=1}^{p}\|\mathbf{x}_{j}-\mathbf{a}_{j}\|_{2}^{2}+\gamma_{1}\sum_{l\in{\cal E}}w_{l}\|\mathbf{v}_{l}\|_{q}+\gamma_{2}\sum_{j=1}^{p}u_{i}\|\mathbf{a}_{j}\|_{2}\\
 &  & +\sum_{l\in\mathbf{{\cal E}}}\langle\mathbf{\lambda}_{l},\mathbf{v}_{l}-A_{i_{1}\cdot}+A_{i_{2}\cdot}\rangle+\frac{\nu}{2}\sum_{l\in{\cal E}}\|\mathbf{v}_{l}-A_{i_{1}\cdot}+A_{i_{2}\cdot}\|_{2}^{2},
\end{eqnarray*}
where $\nu$ is a small constant, $\mathbf{V}=(\mathbf{v}_{1},\ldots,\mathbf{v}_{|{\cal E}|})$,
and $\mathbf{\Lambda}=(\mathbf{\lambda}_{1},\ldots,\mathbf{\lambda}_{|{\cal E}|})$.
Compared with the original algorithms proposed in \citet{Chi2015},
it becomes challenging to deal with the feature-level and observation-level
vectors in the new objective function simultaneously.


\subsubsection{S-ADMM}

S-ADMM minimizes the augmented Lagrangian problem by alternatively
solving one block of variables at a time. In particular, S-ADMM solves
\begin{eqnarray}
\mathbf{A}^{m+1} & = & \mathop{\rm argmin}_{\mathbf{A}}{\cal L}_{\nu}(\mathbf{A},\mathbf{V}^{m},\mathbf{\Lambda}^{m}),\nonumber \\
\mathbf{V}^{m+1} & = & \mathop{\rm argmin}_{\mathbf{V}}{\cal L}_{\nu}(\mathbf{A}^{m+1},\mathbf{V},\mathbf{\Lambda}^{m}),\\
{\mbox{\boldmath $\lambda$}}_{l}^{m+1} & = & {\mbox{\boldmath $\lambda$}}_{l}^{m}+\nu({\bf v}}  \def{\bf V}{{\bf V}_{l}^{m+1}-A_{i_{1}\cdot}^{m+1}+A_{i_{2}\cdot}^{m+1}), \ l \in {\cal E}.\nonumber
\end{eqnarray}


Next we discuss the detailed updating implementations for $\mathbf{A},\mathbf{V}$
and $\mathbf{\Lambda}$ in three steps. A summary of the S-ADMM algorithm
is shown in Algorithm \ref{algADMM}.

\textbf{Step 1: update $\mathbf{A}$.} Denote $\widetilde{{\bf v}}  \def{\bf V}{{\bf V}}_{l}={\bf v}}  \def{\bf V}{{\bf V}_{l}+\frac{1}{\nu}{\mbox{\boldmath $\lambda$}}_{l}$.
Updating ${\bf A}$ is equivalent to minimizing
\begin{eqnarray}
f({\bf A})=\frac{1}{2}\sum_{j=1}^{p}\|{\bf x}}  \def{\bf X}{{\bf X}_{j}-{\bf a}}  \def{\bf A}{{\bf A}_{j}\|_{2}^{2}+\frac{\nu}{2}\sum_{l\in{\cal E}}\|\widetilde{{\bf v}}  \def{\bf V}{{\bf V}}_{l}-A_{i_{1}\cdot}+A_{i_{2}\cdot}\|_{2}^{2}+\gamma_{2}\sum_{j=1}^{p}u_{j}\|{\bf a}}  \def{\bf A}{{\bf A}_{j}\|_{2}.\label{eqn:update_A}
\end{eqnarray}


This optimization problem is challenging because the objective function involves both rows and columns
of the matrix ${\bf A}$. To tackle this difficulty, the following key lemma associates $(\ref{eqn:update_A})$
with a group-lasso regression problem which can be efficiently
solved via standard packages.

\begin{lemma} \label{thm:equivalent}
Let ${\bf I}_{n}$ be an $n\times n$ identity
matrix, ${\bf 1}_{n}$ be an $n$-dimensional
vector with each component being 1, and ${\bf e}}  \def{\bf E}{{\bf E}_{i}$ be an $n$-dimensional
vector with each component being 0 but its $i$-th component being
1. Define ${\bf N}^{-1}=(1+n\nu)^{-1/2}[{\bf I}_{n}+n^{-1}(\sqrt{1+n\nu}-1){\bf 1}_{n}{\bf 1}_{n}^{\rm T}]$
and denote ${\bf y}}  \def{\bf Y}{{\bf Y}_{j}={\bf N}^{-1}[{\bf x}}  \def{\bf X}{{\bf X}_{j}+\nu\sum_{l\in{\cal E}}\widetilde{v}_{jl}({\bf e}}  \def{\bf E}{{\bf E}_{i_{1}}-{\bf e}}  \def{\bf E}{{\bf E}_{i_{2}})]$
with $\widetilde{v}_{jl}$ the $j$-th element of $\widetilde{{\bf v}}  \def{\bf V}{{\bf V}}_{l}$. Then,
minimizing $(\ref{eqn:update_A})$ is equivalent to
\begin{eqnarray*}
\min_{{\bf a}}  \def{\bf A}{{\bf A}_{j}}\frac{1}{2}\|{\bf y}}  \def{\bf Y}{{\bf Y}_{j}-{\bf N}{\bf a}}  \def{\bf A}{{\bf A}_{j}\|_{2}^{2}+\gamma_{2}u_{j}\|{\bf a}}  \def{\bf A}{{\bf A}_{j}\|_{2}, \textrm{~for~each~} j=1,\ldots,p.
\end{eqnarray*}
\end{lemma}

The proof of Lemma \ref{thm:equivalent} is discussed in Appendix.
The key ingredient in the proof is a newly established property of
a permutation matrix, i.e., Proposition \ref{prop1}. Based on this
property, we are able to solve the minimization of $f({\bf A})$ by $p$
separate sub-optimization problems. This together with the property
of group-lasso penalty leads to desirable results. Recall that we require the feature vectors are centered, so we center corresponding estimates during each iteration.

\textbf{Step 2: update ${\bf V}$.} For any $\sigma>0$ and norm $\Omega(\cdot)$,
we define a proximal map,
\begin{eqnarray*}
\textrm{prox}_{\sigma\Omega}({\bf u}}  \def{\bf U}{{\bf U})=\mathop{\rm argmin}_{{\bf v}}  \def{\bf V}{{\bf V}}\left[\sigma\Omega({\bf v}}  \def{\bf V}{{\bf V})+\frac{1}{2}\|{\bf u}}  \def{\bf U}{{\bf U}-{\bf v}}  \def{\bf V}{{\bf V}\|_{2}^{2}\right].
\end{eqnarray*}
In S-ADMM, $\Omega(\cdot)$ is a $q$-norm $\|\cdot\|_{q}$ with $q=1,2$,
or $\infty$, and $\sigma=\gamma_{1}w_{l}/\nu$. We refer the readers to Table 1 of \citet{Chi2015} for the explicit formulations
of the proximal map of $q$-norm for $q=1,2$ and $\infty$. Because vectors ${{\bf v}}  \def{\bf V}{{\bf V}}_{l}$ are
separable, they can be solved via proximal maps, that is
\begin{eqnarray*}
{\bf v}}  \def{\bf V}{{\bf V}_{l} & = & \mathop{\rm argmin}_{{\bf v}}  \def{\bf V}{{\bf V}_{l}}\frac{1}{2}\|{\bf v}}  \def{\bf V}{{\bf V}_{l}-(A_{i_{1}\cdot}-A_{i_{2}\cdot}-\nu^{-1}{\mbox{\boldmath $\lambda$}}_{l})\|_{2}^{2}+\frac{\gamma_{1}w_{l}}{\nu}\|{\bf v}}  \def{\bf V}{{\bf V}_{l}\|_{q}\\
 & = & \textrm{prox}_{\sigma_{l}\|\cdot\|_{q}}(A_{i_{1}\cdot}-A_{i_{2}\cdot}-\nu^{-1}{\mbox{\boldmath $\lambda$}}_{l}).
\end{eqnarray*}


\textbf{Step 3: update ${\mbox{\boldmath $\Lambda$}}$.} Finally, ${\mbox{\boldmath $\lambda$}}_{l}$ can
be updated by ${\mbox{\boldmath $\lambda$}}_{l}={\mbox{\boldmath $\lambda$}}_{l}+\nu({\bf v}}  \def{\bf V}{{\bf V}_{l}-A_{i_{1}\cdot}+A_{i_{2}\cdot})$.

\begin{algorithm}[!htb]
\protect\caption{\quad{}S-ADMM \label{algADMM}}

\begin{enumerate}
\item Initialize $\mathbf{V}^{0}$ and $\mathbf{\Lambda}^{0}$. For $m=1,2,\ldots$
\item For $j=1,\ldots,p$, do
\begin{eqnarray*}
\widetilde{{\bf v}}  \def{\bf V}{{\bf V}}_{l}^{m-1} & = & {\bf v}}  \def{\bf V}{{\bf V}_{l}^{m-1}+\frac{1}{\nu}{\mbox{\boldmath $\lambda$}}_{l}^{m-1}, l\in{\cal E} \\
{\bf y}}  \def{\bf Y}{{\bf Y}_{j}^{m-1} & = & {\bf N}^{-1}\left({\bf x}}  \def{\bf X}{{\bf X}_{j}+\nu\sum_{l\in{\cal E}}\widetilde{v}_{lj}^{m-1}({\bf e}}  \def{\bf E}{{\bf E}_{i_{1}}-{\bf e}}  \def{\bf E}{{\bf E}_{i_{2}})\right),\\
{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m} & = & \mathop{\rm argmin}_{{\bf a}}  \def{\bf A}{{\bf A}_{j}}\frac{1}{2}\|{\bf y}}  \def{\bf Y}{{\bf Y}_{j}^{m-1}-{\bf N}{\bf a}}  \def{\bf A}{{\bf A}_{j}\|_{2}^{2}+\gamma_{2}u_{j}\|{\bf a}}  \def{\bf A}{{\bf A}_{j}\|_{2}, \\
{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m} &=& {\bf a}}  \def{\bf A}{{\bf A}_{j}^{m} - \widebar{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m}{\bf 1}_n, \ {\rm where} \ \widebar{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m}={\bf 1}_n^T {\bf a}}  \def{\bf A}{{\bf A}_{j}^{m}/n.
\end{eqnarray*}

\item For $l\in{\cal E}$, do
\begin{eqnarray*}
{\bf v}}  \def{\bf V}{{\bf V}_{l}^{m}=\textrm{prox}_{\sigma_{l}\|\cdot\|_{q}}(A_{i_{1}\cdot}^{m}-A_{i_{2}\cdot}^{m}-\nu^{-1}{\mbox{\boldmath $\lambda$}}_{l}^{m-1}).
\end{eqnarray*}

\item For $l\in{\cal E}$, do
\begin{eqnarray*}
{\mbox{\boldmath $\lambda$}}_{l}^{m}={\mbox{\boldmath $\lambda$}}_{l}^{m-1}+\nu({\bf v}}  \def{\bf V}{{\bf V}_{l}^{m}-A_{i_{1}\cdot}^{m}+A_{i_{2}\cdot}^{m}).
\end{eqnarray*}

\item Repeat Steps 2-4 until convergence. \end{enumerate}
\end{algorithm}



\subsubsection{S-AMA}

To increase the computational efficiency, we introduce another algorithm S-AMA for implementing sparse convex clustering. S-AMA is different
from S-ADMM in the update of ${\bf A}$. In particular, S-AMA solves
${\bf A}$ by treating $\nu=0$, i.e., ${\bf A}^{m+1}=\mathop{\rm argmin}_{{\bf A}}{\cal L}_{0}({\bf A},{\bf V}^{m},{\mbox{\boldmath $\Lambda$}}^{m})$.
When $\nu=0$, we have $\mathbf{N}=\mathbf{I}_{n}$ and ${\bf y}}  \def{\bf Y}{{\bf Y}_{j}={\bf x}}  \def{\bf X}{{\bf X}_{j}$.
According to Lemma \ref{thm:equivalent}, updating ${\bf A}$ requires to solve $p$
group-lasso problems:
\begin{equation}\label{eq:g1}
\min_{{\bf a}}  \def{\bf A}{{\bf A}_{j}}\frac{1}{2}\|\mathbf{x}_{j}-\mathbf{a}_{j}\|_{2}^{2}+\gamma_{2}u_{j}\|\mathbf{a}_{j}\|_{2}, j=1,\ldots,p.
\end{equation}
By Karush-Kuhn-Tucker (KKT) conditions of the group lasso problem \citep{Yuan2006}, the solution
to (\ref{eq:g1}) has a closed form as
\begin{eqnarray*}
\widehat{\bf a}}  \def{\bf A}{{\bf A}_{j}=\left(1-\frac{\gamma_{2}u_{j}}{\|{\bf z}}  \def{\bf Z}{{\bf Z}_{j}\|_{2}}\right)_{+}{\bf z}}  \def{\bf Z}{{\bf Z}_{j},
\end{eqnarray*}
where ${\bf z}}  \def{\bf Z}{{\bf Z}_{j}={\bf x}}  \def{\bf X}{{\bf X}_{j}+\sum_{l\in{\cal E}}\lambda_{jl}({\bf e}}  \def{\bf E}{{\bf E}_{i_{1}}-{\bf e}}  \def{\bf E}{{\bf E}_{i_{2}})$
and $(z)_{+}=\max\{0,z\}$. See the detailed derivations in online Supplementary. Still, we center $\widehat{\bf a}}  \def{\bf A}{{\bf A}_j$ for each $j$.
The above formula significantly reduces the computational cost by solving
$p$ group-lasso problem analytically in each iteration. Note that
the above update of ${\bf A}$ is independent of ${\bf V}$, which indicates
that S-AMA algorithm does not need to compute the update of ${\bf V}$.
Therefore S-AMA is much more efficient than S-ADMM algorithm.

Next, we discuss the update of ${\mbox{\boldmath $\Lambda$}}$. Define ${\cal P}_{B}({\bf z}}  \def{\bf Z}{{\bf Z})$
as a projection onto $B=\{{\bf y}}  \def{\bf Y}{{\bf Y}:\|{\bf y}}  \def{\bf Y}{{\bf Y}\|_{\dag}\leq1\}$ of the norm
$\|\cdot\|_{\dag}$, where $\|\cdot\|_{\dag}$ is the dual norm of
$\|\cdot\|_{q}$, which defines the fusion penalty. We show in
online Supplementary that the update of ${\mbox{\boldmath $\Lambda$}}$ reduces to ${\mbox{\boldmath $\lambda$}}_{l}^{m}= {\cal P}_{C_{l}}[{\mbox{\boldmath $\lambda$}}_{l}^{m-1}-\nu(A_{i_{1}\cdot}^{m}-A_{i_{2}\cdot}^{m})]$
with $C_{l}=\{{\mbox{\boldmath $\lambda$}}_{l}:\|{\mbox{\boldmath $\lambda$}}_{l}\|_{\dag}\leq\gamma_{1}w_{l}\}$.
The S-AMA algorithm is summarized in Algorithm \ref{algAMA}.

\begin{algorithm}[!htb]
\protect\caption{\quad{}S-AMA \label{algAMA}}

\begin{enumerate}
\item Initialize $\mathbf{\Lambda}{}^{0}$. For $m=1,2,\ldots$
\item For $j=1,\ldots,p$, do
\begin{eqnarray*}
{\bf z}}  \def{\bf Z}{{\bf Z}_{j}^{m} & = & {\bf x}}  \def{\bf X}{{\bf X}_{j}+\sum_{l \in{\cal E}}\lambda_{lj}^{m-1}({\bf e}}  \def{\bf E}{{\bf E}_{i_{1}}-{\bf e}}  \def{\bf E}{{\bf E}_{i_{2}}),\\
{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m} & = & \left(1-\frac{\gamma_{2}u_{i}}{\|{\bf z}}  \def{\bf Z}{{\bf Z}_{i}^{m}\|_{2}}\right)_{+}{\bf z}}  \def{\bf Z}{{\bf Z}_{j}^{m},\\
{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m} &=& {\bf a}}  \def{\bf A}{{\bf A}_{j}^{m} - \widebar{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m}{\bf 1}_n, \ {\rm where} \ \widebar{\bf a}}  \def{\bf A}{{\bf A}_{j}^{m}={\bf 1}_n^T {\bf a}}  \def{\bf A}{{\bf A}_{j}^{m}/n.
\end{eqnarray*}

\item For $l\in{\cal E}$, do
\[
{\mbox{\boldmath $\lambda$}}_{l}^{m}={\cal P}_{C_{l}}[{\mbox{\boldmath $\lambda$}}_{l}^{m-1}-\nu(A_{i_{1}\cdot}^{m}-A_{i_{2}\cdot}^{m})],
\]
where $C_{l}=\{{\mbox{\boldmath $\lambda$}}_{l}:\|{\mbox{\boldmath $\lambda$}}_{l}\|_{\dag}\leq\gamma_{1}w_{l}\}$.
\item Repeat Steps 2-3 until convergence. \end{enumerate}
\end{algorithm}



\subsubsection{Algorithmic Convergence}

This subsection discusses the convergence of the proposed S-ADMM and
S-AMA algorithms. \citet{Chi2015} and the references therein provided sufficient conditions for the convergence
of the following general optimization problem,
\begin{eqnarray}
\min_{\xi,\zeta}\ f(\xi)+g(\zeta),\textrm{\ \ s.t.\ \ }A\xi+B\zeta=c.\label{general_obj}
\end{eqnarray}
They verified that the ADMM and AMA algorithms for
convex clustering, as two special cases of \eqref{general_obj}, satisfied
the sufficient conditions under which the convergence was guaranteed.

The convergence of our S-ADMM and S-AMA algorithms follows similar
arguments. Note that the only difference between the objective function
in \eqref{eq:obj_constraint} and its counterpart in \citet{Chi2015}
is a convex penalty term $\gamma_{2}\sum_{j=1}^{p}u_{j}\|\mathbf{a}_{j}\|_{2}$.
Define the summation of the first and third terms of the objective
function in \eqref{eq:obj_constraint} as $f(\cdot)$, and the second
term as $g(\cdot)$. This indicates that the optimization problem
\eqref{eq:obj_constraint} is a special case of \eqref{general_obj}.
Simple algebra implies that $f(\cdot)$ is strongly convex. According
to \citet{Chi2015}, one can show that, under mild regularization conditions, the convergence
of S-ADMM is guaranteed for any $\nu>0$, and the convergence
of S-AMA algorithm is guaranteed provided that positive
constant $\nu$ is not too large.

\subsubsection{Computational Consideration}

Step 2 in both Algorithms \ref{algADMM} and \ref{algAMA} involves $p$ sub-optimization problems. Therefore, S-ADMM and S-AMA merit from the distributed optimization, and they can handle large-scale problems efficiently. To be specific, Step 2 can be distributed to different processors to obtain estimates of ${\bf a}}  \def{\bf A}{{\bf A}_j$'s which are then gathered to update ${\bf A}$. In addition, Steps 3-4 in Algorithm \ref{algADMM} or Step 3 in Algorithm \ref{algAMA} can also be distributed to different processors to obtain fast updates.

It is worth pointing out that the computation of S-AMA is comparable to AMA in \cite{Chi2015}, while S-ADMM is computationally more expensive than ADMM in \cite{Chi2015} and S-AMA. This is because Step 2 in S-ADMM does not have a closed-form formula and it requires solving $p$ group-lasso problems assisted by iterations. Furthermore, Step 3 in S-AMA only requires updates for ${\mbox{\boldmath $\lambda$}}_l$ for those $l$ such that $w_l > 0, l \in {\cal E}$. With a suitable selection of $w_l$ discussed in Section \ref{sec:weights}, the size of working set of $l$ can be dramatically reduced from $n(n-1)/2$ to a much smaller number. Our limited experience in numerical studies also confirms the superiority of S-AMA over S-ADMM in terms of the computational cost.

We have developed an $\textsf{R}$ package ``\textit{scvxclustr}" to facilitate the implementation of proposed methods. Table \ref{tb_time} compares the computational time in seconds of our package with \cite{Chi2015}'s method via ``\textit{cvxclustr}" for the setting with 4 clusters described in Section \ref{sec:simulation} with given tuning parameters. The computer is equipped with a CPU i3-4170 (3.70GHz) and 8G memory.

\begin{table}[!htb]
\centering \protect\caption{Timing comparison under the 4 cluster setting described in Section \ref{sec:simulation} with given tuning parameters under various combinations of $n$ and $p$. The time is in second.}
\label{tb_time} \vspace{6pt}
 \centering  \small
\begin{tabular}{L{1.7cm}R{0.7cm}R{0.7cm}R{0.7cm}p{0.01cm}R{0.7cm}R{0.7cm}R{0.7cm}p{0.01cm}R{0.8cm}R{0.8cm}R{0.8cm}R{0.01cm}R{0.7cm}R{0.7cm}R{0.7cm}}
\hline
$n$  & \multicolumn{3}{c}{60} & \  & \multicolumn{3}{c}{120} & \  & \multicolumn{3}{c}{300} & \  & \multicolumn{3}{c}{400}\tabularnewline
\hline
$p$ & 150  & 300  & 500  & \  & 150  & 300  & 500  & \  & 150  & 300  & 500 & \ & 150  & 300  & 500 \tabularnewline
\hline
\hline
AMA	& 0.6 &	3 &	5 & \ &	3 &	7 &	13 & \ &	11 &	23 &	37 & \ &	15 &	32 &	51 \tabularnewline
ADMM &	6 &	17 &	24 & \ &	19 &	71 &	82 & \ &	157 &	375 &	548 & \ &	329 &	1482	& 1796 \tabularnewline
S-AMA &	0.1 &	0.6 &	2 & \ &	0.8 &	2 &	6 & \ &	3 & 10 &	25 & \ &	6 &	17 &	34 \tabularnewline
S-ADMM	& 16	& 47 &	60 & \ &	60 &	152 &	215 & \ &	593 &	2051 &	2781 & \ &	1252 &	3639 &	7925 \tabularnewline
\hline
\end{tabular}
\end{table}

\section{Theoretical Properties} \label{sec:thm}

In this section, we study statistical properties of the proposed sparse convex clustering. In particular, we provide finite sample bounds for prediction error of the proposed sparse convex clustering estimator and establish its variable selection consistency. In Section S.3 of online Supplementary, we also develop unbiased estimators for the degrees of freedoms of sparse convex clustering.



Assume ${\bf x}}  \def{\bf X}{{\bf X}={\bf a}}  \def{\bf A}{{\bf A}_0 + {\mbox{\boldmath $\varepsilon$}}$, where ${\mbox{\boldmath $\varepsilon$}} \in {\mathbb R}^{np}$ is a vector of
independent sub-Gaussian noise terms with mean zero and variance $\sigma^2$, and ${\bf a}}  \def{\bf A}{{\bf A}_0={\rm vec}({\bf A}_0)=({\bf a}}  \def{\bf A}{{\bf A}_{01}^{\rm T},\ldots, {\bf a}}  \def{\bf A}{{\bf A}_{0p}^{\rm T})^{\rm T}$
is a $np$-dimensional mean vector. Without loss of generality, we assume that only the first $p_0 < p$
features are informative, i.e., $\| {\bf a}}  \def{\bf A}{{\bf A}_{0j} \|_2 \neq 0$ for $j \leq p_0$ and $\| {\bf a}}  \def{\bf A}{{\bf A}_{0j} \|_2 = 0$ for $j > p_0$.
The informative feature set is denoted as ${\cal A}=\{1,\ldots,p_0 \}$ and the noninformative feature set
is ${\cal A}^c = \{p_0+1,\ldots,p \}$. For simplicity, we consider the case with $w_l=1$.

Our sparse convex clustering in (\ref{eq:obj_constraint}) can be reformulated as the following problem:
\begin{eqnarray} \label{eq:obj_bound}
\widehat{{\bf a}}  \def{\bf A}{{\bf A}} = \mathop{\rm argmin}_{{\bf a}}  \def{\bf A}{{\bf A} \in {\mathbb R}^{np}} \frac{1}{2} \| {\bf x}}  \def{\bf X}{{\bf X} - {\bf a}}  \def{\bf A}{{\bf A} \|_2^2 + \gamma_1 \sum_{l \in {\cal E}} \|{\bf C}_l {\bf a}}  \def{\bf A}{{\bf A} \|_q
+ \gamma_2 \sum_{j=1}^{p} u_j \| {\bf a}}  \def{\bf A}{{\bf A}_j \|_2,
\end{eqnarray}
where ${\bf C}_l={\bf I}_p \otimes ({\bf e}}  \def{\bf E}{{\bf E}_{i_1}-{\bf e}}  \def{\bf E}{{\bf E}_{i_2})^{\rm T}$ and hence ${\bf C}_l {\bf a}}  \def{\bf A}{{\bf A} = A_{i_1 \cdot}- A_{i_2 \cdot}$. Define ${\bf C} = ( {\bf C}_1 ^{\rm T}, \ldots, {\bf C}_{|{\cal E}|}^{\rm T} )^{\rm T}$ and denote ${\bf u}}  \def{\bf U}{{\bf U}=(u_1,\ldots,u_p)^{\rm T}$.

The following two theorems provide the finite sample bounds for prediction error of our sparse convex clustering estimator with parameter $q \in \{1,2\}$, respectively.

\begin{theorem} \label{thm:bound1}
Let $\widehat{\bf a}}  \def{\bf A}{{\bf A}$ be the estimate of \eqref{eq:obj_bound} with $q=1$. If $\gamma_1 > 4 \sigma \sqrt{\frac{\hbox{log} (p \cdot \binom{n}{2})}{n}}$, then
\begin{eqnarray*}
\frac{1-\gamma_2}{2np}  \| \widehat{\bf a}}  \def{\bf A}{{\bf A} - {\bf a}}  \def{\bf A}{{\bf A}_0 \|_2^2  \leq
\frac{3\gamma_1}{2np} \| {\bf C} {\bf a}}  \def{\bf A}{{\bf A}_0 \|_1 + \frac{\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_2^2 }{2np} +\sigma^2 \left[
\frac{1}{n} + \sqrt{\frac{\hbox{log} (np)}{n^2 p}} \right]  + \frac{1}{np}
\end{eqnarray*}
holds with probability at least $1-c_3$, where
\begin{eqnarray*}
c_3=\frac{2}{p \cdot \binom{n}{2}} +\exp\{ - \min (c_1 \hbox{log} (np),c_2 \sqrt{p \hbox{log} (np)} )\} + 2 \exp \big( -np/(2 \sigma^2 \gamma_2^2  \| {\bf u}}  \def{\bf U}{{\bf U} \|_1^2) \big),
\end{eqnarray*}
for some positive constants $c_1$ and $c_2$ defined in Lemma \ref{lemma:HW71}.
\end{theorem}



\begin{theorem} \label{thm:bound2}
Let $\widehat{\bf a}}  \def{\bf A}{{\bf A}$ be the estimate of \eqref{eq:obj_bound} with $q=2$. If $\gamma_1 > 4 \sigma \sqrt{\frac{\hbox{log} (p \cdot \binom{n}{2})}{n}}$, then
\begin{eqnarray*}
\frac{1-\gamma_2}{2np}  \| \widehat{\bf a}}  \def{\bf A}{{\bf A} - {\bf a}}  \def{\bf A}{{\bf A}_0 \|_2^2  \leq
\frac{3\gamma_1}{2np} \sum_{l\in{\cal E}} \| {\bf C}_l {\bf a}}  \def{\bf A}{{\bf A}_0 \|_2+ \frac{\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_2^2 }{2np} +\sigma^2 \left[
\frac{1}{n} + \sqrt{\frac{\hbox{log} (np)}{n^2 p}} \right]  + \frac{1}{np}
\end{eqnarray*}
holds with probability at least $1-c_3$, where $c_3$ is defined in Theorem \ref{thm:bound1}.
\end{theorem}



\begin{remark}
Based on the sparsity assumption of features, we know true underlying clusters differ only with respect to the first $p_0$ features. Thus, $\| {\bf C}_l {\bf a}}  \def{\bf A}{{\bf A}_0 \|_1=O(1)$ and $\| {\bf C}_l {\bf a}}  \def{\bf A}{{\bf A}_0 \|_2=O(1)$, $\forall l \in {\cal E}$. Note that $| {\cal E}| = \binom{n}{2} $, and hence $\|{\bf C} {\bf a}}  \def{\bf A}{{\bf A}_0 \|_1 = O(n^2)$ and $\sum_{l \in {\cal E}} \| {\bf C}_l {\bf a}}  \def{\bf A}{{\bf A}_0 \|_2 =O(n^2)$. In order to obtain a general prediction consistency based on Theorems \ref{thm:bound1}-\ref{thm:bound2}, we need $\gamma_2 \to 0, c_3 \to 0$ and the first two terms on right-hand side of the inequalities in Theorems \ref{thm:bound1}-\ref{thm:bound2} vanish. Assume $\gamma_1 \| {\bf C} {\bf a}}  \def{\bf A}{{\bf A}_0 \|_1/ (2np)=o(1)$, and then we have $\sqrt{ n \hbox{log} (p \cdot \binom{n}{2})/p^2}=o(1)$. Additionally, $c_3 \to 0$ and $\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_2^2 / (2np)$ are equivalent to requiring $\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_1^2 /(np) \to 0$ by noting that $\| {\bf u}}  \def{\bf U}{{\bf U} \|_2 \leq \| {\bf u}}  \def{\bf U}{{\bf U} \|_1$.

In particular, we discuss following two choices for the weight $u_j$, which lead to prediction consistency. First, we assume a non-adaptive weight, i.e., $u_j\equiv1$. It is easy to see $\| {\bf u}}  \def{\bf U}{{\bf U} \|_1=p$. Additionally, we require $\gamma_2 \to 0$ and $\gamma_2 p n^{-1} \to 0$. Second, more generally, we assume $ \| {\bf u}}  \def{\bf U}{{\bf U} \|_1$ is bounded above by $p^{\zeta}$, where $\zeta$ is a constant, and require $\gamma_2 \to 0$, $\gamma_2 p^{2\zeta-1} n^{-1} \to 0$. Therefore, in either case $\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_1^2 /(np) \to 0$ and then $\widehat{\bf a}}  \def{\bf A}{{\bf A}$ is prediction consistent with $q=1$ or $q=2$.
\end{remark}




Next, we establish the asymptotic selection consistency of the proposed sparse convex clustering, which is a desirable property in high-dimensional cluster analysis where many features are non-informative.

\begin{theorem} \label{thm:selection_consistency}
If $\gamma_1 > 4 \sigma \sqrt{\hbox{log} (p \cdot \binom{n}{2})/n}$, $\gamma_1 \| {\bf C} {\bf a}}  \def{\bf A}{{\bf A}_0 \|_1/ (2np)=o(1)$, $\gamma_2 \to 0$ and $\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_1^2 /(np) \to 0$ as $n,p \to \infty$, then $P(\|\widehat{\bf a}}  \def{\bf A}{{\bf A}_j \|_2=0)\to 1$ for any $j \in {\cal A}^c$, with the solution $\widehat{\bf a}}  \def{\bf A}{{\bf A}$ to \eqref{eq:obj_bound} with either $q=1$ or $q=2$.
\end{theorem}

\begin{remark}
Condition $\gamma_2 \| {\bf u}}  \def{\bf U}{{\bf U} \|_1^2 /(np) \to 0$ generally implies that the adaptive weights cannot be too large. For example, uniform weights satisfy this condition. Note that Conditions $\gamma_1 > 4 \sigma \sqrt{\hbox{log} (p \cdot \binom{n}{2})/n}$ and $\gamma_1 \| {\bf C} {\bf a}}  \def{\bf A}{{\bf A}_0 \|_1/ (2np)=o(1)$ imply $\sqrt{ n \hbox{log} (p \cdot \binom{n}{2})/p^2}=o(1)$. This is derived to ensure the estimation consistency in Theorems \ref{thm:bound1} and \ref{thm:bound2}, in particular, to ensure the first term on the right-hand side of Theorem \ref{thm:bound1} or \ref{thm:bound2} to converge to zero. This condition requires that $n = o(p^2)$ up to a $\hbox{log}$-term, and hence is satisfied as long as $p$ diverges not too slow, which is typically true in high-dimensional scenarios. Note that in this case $p$ can still be smaller than $n$. Similar phenomenon has also been found in \citet{Tan2015}.
\end{remark}

Theorem \ref{thm:selection_consistency} establishes the asymptotic selection consistency in the sense that
the proposed sparse convex clustering can eliminate the non-informative variables in
the estimated cluster centers with probability tending to one.This variable selection consistency is illustrated in the motivation example shown in Figure \ref{heatmap}.


Proofs of Theorems \ref{thm:bound1}-\ref{thm:selection_consistency} are provided in online Supplementary.


\section{Practical Issues}

\label{sec:practical}

In Section \ref{sec:algorithm}, the S-ADMM and S-AMA algorithms rely
on the choice of weights and the tuning parameters $\gamma_{1}$ and
$\gamma_{2}$. In this section, we discuss how to choose these parameters
in practice.


\subsection{Selection of Weights}
\label{sec:weights}

This subsection introduces practical selections of the weights $w_{i_{1},i_{2}}$,
$(i_{1},i_{2})\in{\cal E}$, in the fused-lasso penalty, and the factors
$u_{j}$, $j=1,\cdots,p$, in the adaptive group-lasso penalty.

Following \citet{Chi2015}, we choose weights by incorporating the
m-nearest-neighbors method with Gaussian kernel. To be specific, the
weight between the sample pair $(i_{1},i_{2})$ is set as $w_{i_{1},i_{2}}=\iota_{i_{1},i_{2}}^{m}\exp(-\phi\|X_{i_{1}\cdot}-X_{i_{2}\cdot}\|_{2}^{2})$,
where $\iota_{i_{1},i_{2}}^{m}$ equals 1 if observation $i_{2}$
is among observation $i_{1}$'s $m$ nearest neighbors or vice versa, and
0 otherwise. This choice of weights works well for a wide range of
$\phi$ when $m$ is small. In our numerical results, $m$ is fixed
at $5$ and $\phi$ is fixed at 0.5.

Next we consider the selection of factor $u_{j}$. As suggested
by \citet{zou:2006}, $u_{j}$ can be chosen as $1/\|\widehat{{\bf a}}  \def{\bf A}{{\bf A}}_{j}^{(0)}\|_{2}$,
where $\widehat{\bf a}}  \def{\bf A}{{\bf A}_{j}^{(0)}$ is the estimate of ${\bf a}}  \def{\bf A}{{\bf A}_{j}$ in \eqref{eq:obj_constraint}
with $\gamma_{2}=0$. Such choice of factors penalizes less on informative
features and penalizes more on uninformative features, and hence leads
to improved clustering accuracy and variable selection performance
than its non-adaptive counterpart.

Finally, in order to ensure that the optimal tuning parameters $\gamma_{1}$
and $\gamma_{2}$ lie in relatively robust intervals regardless of
feature dimension and sample size, weights $w_{i_{1},i_{2}}$ and
factors $u_{j}$ are re-scaled to sum to $1/\sqrt{p}$ and $1/\sqrt{n}$,
respectively. Such re-scaling is only for convenience and does not
affect the final clustering path.


\subsection{Selection of Tuning Parameters}

\label{sec:tuning}

This subsection provides a selection method for tuning parameters
$\gamma_{1}$ and $\gamma_{2}$. Remind that $\gamma_{1}$ controls
the number of estimated clusters and $\gamma_{2}$ controls the number
of selected informative features.

We first illustrate via a toy example the effectiveness of tuning
parameter $\gamma_{2}$ on variable selection accuracy. In this
example, $60$ observations with $p=500$ features are generated from
4 clusters. Among all the features, only $20$ variables differ between
clusters. See detailed simulation setup in Section \ref{sec:simulation}.
By fixing $\widehat\gamma_{1}=2.44$ and varying $\gamma_{2}$ from $e^{-5.0}$
to $e^{7.0}$, we plot the path of false negative rate (FNR) and the
path of false positive rate (FPR) of the final estimator. As shown
in Figure \ref{demo-graph}, when $\gamma_{2}$ is close to zero,
all features are included, and when $\gamma_{2}$ increases to some
ranges of intervals, all and only uninformative features are excluded,
i.e., perfect variable selection performance. This illustrates the
sensitivity of $\gamma_{2}$ to the variable selection performance
of the final estimator. In practice, we aim to estimate a suitable
$\gamma_{2}$ that leads to satisfactory variable selection.

\begin{figure}[!htb]
\protect\caption{Illustration of the effectiveness of $\gamma_{2}$ on variable selection
accuracy. The solid curve is the path of false negative rate (FNR),
and the dashed curve is the path of false positive rate (FPR).}
\centering \includegraphics[scale=0.5]{Figure2.pdf}\label{demo-graph}
\end{figure}

In literature, \citet{Wang:2010} and \citet{Fang2012} proposed stability
selection to estimate the tuning parameters in clustering models.
The idea behind stability selection is that a good tuning parameter
should produce clustering results that are stable with respect to
a small perturbation to the training samples. Stability selection well suits
the model selection in cluster analysis because cluster labels are unavailable and the cross-validation method is not applicable in this case.

In this paper, we propose
to use stability selection in \citet{Fang2012} to tune both parameters
$\gamma_{1}$ and $\gamma_{2}$. To be specific, for any given $\gamma_{1}$
and $\gamma_{2}$, based on two sets of bootstrapped samples, two
clustering results can be produced via \eqref{eq:obj_constraint}, and
then the stability measurement \citep{Fang2012} can be computed to
measure the agreement between the two clustering results. In order
to enhance the robustness of the stability selection method, we repeat
this procedure $50$ times and then compute the averaged stability
value. Finally, the optimal parameter is selected as the one achieving
maximum stability. Our extensive numerical studies show that the selection of
important features is less sensitive to the clustering path, i.e., import features stand out in almost all clustering structures.
Thus, to speed up tuning process, stability path can be computed over
of a coarse grid of $\gamma_{1}$ and a fine grid of $\gamma_{2}$.


\section{Numerical Results}

\label{sec:experiment}

This section demonstrates the superior performance of our sparse convex
clustering in simulated examples in Section \ref{sec:simulation}
and a real application of hand movement clustering in Section \ref{sec:real}.


\subsection{Simulation Studies}

\label{sec:simulation}

In this subsection, simulations studies are conducted to evaluate
the performance of sparse convex clustering methods (S-ADMM and S-AMA). They
are compared to the k-means clustering and two convex clustering algorithms: ADMM and AMA \citep{Chi2015}.

First, we consider four spherical settings. Each simulated dataset consists
of $n=60$ observations with the number of clusters either $K=2$
or 4, and the number of features either $p=150$ or $500$. In each
setting, only the first 20 features are informative and remaining
features are non-informative. The samples $X_{i\cdot}\in{\mathbb R}^{p},i=1,\ldots,n$,
are generated as follows. Denote a $p$-dimensional multivariate normal distribution as $\textrm{MVN}_{p}$. For each $i$, a cluster label $Z_{i}$
is uniformly sampled from $\{1,\ldots,K\}$, and then the first 20
informative features are generated from $\textrm{MVN}_p({\mbox{\boldmath $\mu$}}_{K}(Z_{i}),{\bf I}_{20})$,
where ${\mbox{\boldmath $\mu$}}_{K}(Z_{i})$ is defined as follows:
\begin{itemize}
\item If $K=2$, ${\mbox{\boldmath $\mu$}}_{2}(Z_{i})=\mu{\bf 1}_{20}I(Z_{i}=1)-\mu{\bf 1}_{20}I(Z_{i}=2)$;
\item If $K=4$, ${\mbox{\boldmath $\mu$}}_{4}(Z_{i})=(\mu{\bf 1}_{10}^{\rm T},-\mu{\bf 1}_{10}^{\rm T})^{\rm T} I(Z_{i}=1)+(-\mu{\bf 1}_{10}^{\rm T},-\mu{\bf 1}_{10}^{\rm T})^{\rm T} I(Z_{i}=2)+\ \ \ \ \ \ \ \ \ \ \ \ \ (-\mu{\bf 1}_{10}^{\rm T},\mu{\bf 1}_{10}^{\rm T})^{\rm T} I(Z_{i}=3)+(\mu{\bf 1}_{10}^{\rm T},\mu{\bf 1}_{10}^{\rm T})^{\rm T} I(Z_{i}=4)$,
\end{itemize}
where $\mu$ controls the
distance between cluster centers. Here a large $\mu$ indicates that
clusters are well-separated, whereas a small $\mu$ indicates that
clusters are overlapped. Finally, the rest $p-20$ noise features
are generated from ${\cal N}(0,1)$.

Second, we consider a non-spherical setting with two half moons. Each simulated dataset consists of $n=100$ observations with $K=2$ clusters and $p=40$ features. Only the first two features are informative, and the rest 38 noisy features are generated from ${\cal N} (0, 0.5)$. It is a relatively hard setting because the number of non-informative features are 19 times more than that of informative features. Figure \ref{twohalfmoons} shows one example of two interlocking half moons with the first two features. For a comparison purpose, we also apply the spectral clustering (SPECC, see \cite{ng2002spectral}) via an $\textsf{R}$ package ``\textit{kernlab}"  because SPECC can internally tackle non-spherical clusters.

\begin{figure}[!htb]
\protect\caption{The plot of the first two features for one example of two interlocking half moons.}
\centering \includegraphics[scale=0.45]{Figure3.pdf} \label{twohalfmoons}
\end{figure}

In summary, five simulation settings are considered. Setting 1: $K=2,n=60, p=150$,
and $\mu=0.6$; Setting 2: $K=2,n=60,p=500$, and $\mu=0.7$; Setting 3:
$K=4,n=60,p=150$, and $\mu=0.9$; Setting 4: $K=4,n=60,p=500$, and $\mu=1.2$; Setting 5: two half moons with $K=2, n=100, p=40$.
For each setting, we run 200 repetitions.

The RAND index \citep{Rand1971} is used to measure the agreement
between the estimated clustering result and the underlying true clustering
assignment. The RAND index ranges between 0 and 1, and a higher value
indicates better performance. Note that the true cluster labels are known in
simulation studies, and thus it is feasible to know how well the candidate methods can
perform if they are tuned by maximizing the RAND index.
To ensure fair comparisons, for each
repetition, separate validation samples are generated and used to
select an optimal $k$ in k-means, an optimal $\gamma$ in ADMM and AMA, and optimal $\gamma_{1}$
and $\gamma_{2}$ in S-ADMM and S-AMA. To evaluate the performance
of variable selection, two measurements are reported: the false negative
rate (FNR) and the false positive rate (FPR). All the simulation results
are summarized in Table \ref{tb_simulation_1}. Due to its relatively
expensive computational costs, S-ADMM is not evaluated for Settings
2 and 4 where $p=500$.

In the first four simulation settings, the centers are spherical and
hence k-means always performs well in clustering accuracy, i.e., large
RAND index. The goals of these simulations are to justify that (1)
convex clustering does not perform well when the feature dimension
is high; (2) sparse convex clustering performs very well when the
feature dimension is high; and (3) sparse convex clustering selects
informative features with great clustering accuracy. The goal of Setting 5 is to show under the non-spherical setting, sparse convex clustering can still perform well, even better than the SPECC. All these claims
are justified by the results presented in Table \ref{tb_simulation_1}.


\begin{table}[!htb]
\centering \protect\caption{Empirical mean and standard deviation (SD) of the RAND index, false
positive rate (FPR), and false negative rate (FNR) based on 200 repetitions.
Setting 1: $K=2,n=60,p=150$, and $\mu=0.6$; Setting 2: $K=2,n=60,p=500$,
and $\mu=0.7$; Setting 3: $K=4,n=60,p=150$, and $\mu=0.9$; Setting 4:
$K=4,n=60,p=500$, and $\mu=1.2$; Setting 5: two half moons with $K=2, n=100, p=40$. The best RAND index in each scenario is shown
in bold. }


\label{tb_simulation_1} \vspace{6pt}
 \centering %
\begin{tabular}{p{2cm}p{2cm}p{1.2cm}p{1.2cm}p{0.1cm}p{1.2cm}p{1.2cm}p{0.1cm}p{1.2cm}p{1.2cm}}
\hline
\  & \  & \multicolumn{2}{c}{RAND} & \  & \multicolumn{2}{c}{FNR} & \  & \multicolumn{2}{c}{FPR}\tabularnewline
\hline
 & Algorithm  & mean  & SD  & \  & mean  & SD  & \  & mean  & SD \tabularnewline
\hline
\hline
Setting 1  & k-means  & 0.95  & 0.06  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & ADMM  & 0.53  & 0.39  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & AMA  & 0.66  & 0.40  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & S-ADMM  & 0.82  & 0.24  & \  & 0.04  & 0.05  & \  & 0.25  & 0.16 \tabularnewline
\  & S-AMA  & \textbf{0.96}  & 0.06  & \  & 0.03  & 0.07  & \  & 0.30  & 0.21 \tabularnewline
\hline
Setting 2  & k-means  & 0.95  & 0.11  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & ADMM  & 0.14  & 0.20  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & AMA  & 0.08  & 0.21  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & S-AMA  & \textbf{0.97}  & 0.07  & \  & 0.07  & 0.09  & \  & 0.11  & 0.10 \tabularnewline
\hline
Setting 3  & k-means  & 0.83  & 0.15  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & ADMM  & 0.56  & 0.22  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & AMA  & 0.47  & 0.21  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & S-ADMM  & 0.82  & 0.14  & \  & 0.04  & 0.06  & \  & 0.25  & 0.24 \tabularnewline
\  & S-AMA  & \textbf{0.84}  & 0.13  & \  & 0.02  & 0.04  & \  & 0.11  & 0.18 \tabularnewline
\hline
Setting 4  & k-means  & 0.89  & 0.14  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & ADMM  & 0.31  & 0.23  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & AMA  & 0.31  & 0.20  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & S-AMA  & \textbf{0.94}  & 0.09  & \  & 0.01  & 0.02  & \  & 0.01  & 0.03 \tabularnewline
\hline
Setting 5  & k-means  & 0.51  & 0.07  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & ADMM  & 0.54  & 0.08   & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & AMA  & 0.53  &  0.09  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00 \tabularnewline
\  & S-AMA  & \textbf{0.57} & 0.07  & \  & 0.00  &  0.00  & \  & 0.34 & 0.27  \tabularnewline
\  & SPECC & 0.52 & 0.08  & \  & 0.00  & 0.00  & \  & 1.00  & 0.00   \tabularnewline
\hline
\end{tabular}
\end{table}


First, convex clustering does not perform well when the feature dimension
is high, even much worse than k-means. Similar phenomenon was also
observed in the simulation studies conducted in \citet{Tan2015}.
This is the motivation for developing sparse convex clustering. Second,
sparse convex clustering improves convex clustering significantly.
Sparse convex clustering (S-AMA) performs as well as k-means when
$p=150$, and performs better than k-means when $p=500$. Third, sparse
convex clustering selects informative feature with great accuracy,
that is, with low FNR and FPR. The feature selection performance of
sparse convex clustering is very promising for settings where $p=500$. Compared with S-ADMM, S-AMA is computationally faster and also delivers slightly better accuracy. Therefore, we recommend
S-AMA in practice. Fourth, for Setting 5, due to its non-spherical essential, convex clustering outperforms k-means. Surprisingly, SPECC does not outstand in terms of the RAND index. The existence of many noninformative features undermines the SPECC and it confirms the necessity of selecting informative features. Sparse convex clustering performs the best by selecting informative features.



\subsection{Application}

\label{sec:real}

We evaluate the performance of sparse convex clustering in LIBRAS
movement data from the Machine Learning Repository \citep{Lichman2013}.
The original dataset contains 15 classes with each class referring
to a hand movement type. Each class contains 24 observations, and
each observation has 90 features consisting of the coordinates of
hand movements. We use this dataset without the clustering assignments
to evaluate each clustering algorithms and then compare the results
with the true classes to compute the RAND index. Before cluster analysis,
each feature is centered. In our S-AMA algorithm, we set $m=5$ and
$\phi=1$ for weight $w_{i_{1},i_{2}}$.

Note that some of the original 15 clusters indicate similar hand movements,
such as curved/vertical swing and horizontal/vertical straight-line.
By plotting the first two principal components of the 90 features,
one can see that some clusters are severely overlapped. Therefore,
for evaluation purpose, six clusters, including vertical swing (labeled as 3),
anti-clockwise arc (labeled as 4), clockwise arc (labeled as 5), horizontal straight-line (labeled as 7),
horizontal wavy (labeled as 11), and
vertical wavy (labeled as 12) in the original dataset are selected. The left panel of Figure \ref{truelabel_cvxpath} displays the
plot of the first principal component (PC1) against the second principal
component (PC2) of 90 features for the selected six clusters with
the true cluster labels.

\begin{figure}[!htb]
\protect\caption{Left panel shows the plot of the first principal component against the second principal component of 90 features for the selected six clusters with true cluster
labels; Right panel shows the clustering path of convex clustering (AMA) using all 90 features by plotting the first principal component (PC1) against the second principal component (PC2).}
\centering \includegraphics[scale=0.8]{Figure4_colored.pdf} \label{truelabel_cvxpath}
\end{figure}


We first display the clustering path of convex clustering (AMA) using
all 90 features in the right panel of Figure \ref{truelabel_cvxpath}. Clearly, convex clustering
is only able to distinguish clusters 4 and 5 and treat the rest clusters
as one class. This phenomenon shows the curse of dimensionality in
high-dimensional clustering and motivates the need to conduct feature
selection for improved clustering performance.

We use S-AMA to solve sparse convex clustering. The tuning parameters
are selected according to the stability selection in Section \ref{sec:tuning}.
Table \ref{tb_realdata_1} reports the number of estimated clusters,
the number of selected features, and the RAND index between the estimated
cluster membership and the true cluster membership for k-means,
AMA algorithm, and our S-AMA algorithm. Clearly, both convex clustering
(AMA) and sparse convex clustering (S-AMA) perform better than k-means,
which indicates that the performance of convex clustering or sparse
convex clustering is less sensitive to the assumption of spherical
clustering centers. In addition, by using only $13$ informative features,
our S-AMA is able to improve the clustering accuracy of convex clustering
(AMA) by $45\%$. This indicates the importance of variable selection
in high-dimensional clustering.

\begin{table}[!htb]
\centering \protect\caption{The number of estimated clusters, the number of selected features,
and the RAND index for k-means, AMA, and our S-AMA algorithm.}


\label{tb_realdata_1} \vspace{6pt}
 \centering %
\begin{tabular}{c|ccc}
\hline
Algorithm  & \# of clusters  & \# of features  & RAND index \\
\hline
k-means  & 2  & 90  & 0.06 \\
AMA  & 3  & 90  & 0.31 \\
S-AMA  & 3  & 13  & \textbf{0.45}\\
\hline
\end{tabular}
\end{table}


Next we demonstrate the clustering path of sparse convex clustering
(S-AMA) with the 13 selected features in Figure \ref{scvxpath}.
Figure \ref{scvxpath} displays three big clusters, which is consistent
with the number of estimated clusters shown in Table \ref{tb_realdata_1}.
As tuning parameter $\gamma_{1}$ increases,
the clustering path of S-AMA tends to merge clusters 3, 7 and 12 into
one big cluster, merge cluster 4 and 5 into another big cluster, and
identify cluster 11 as the third cluster. This finding is displayed in the final clustering
path of S-AMA executed at the selected $\gamma_{1}$ and $\gamma_{2}$
as shown in Figure \ref{truebestlabel}. In the plot, the left-panel
graph shows the true cluster labels and the right-panel graph shows
the three estimated clusters using S-AMA.

\begin{figure}[!htb]
\protect\caption{The clustering path of sparse convex clustering (S-AMA) using only
13 selected features by plotting PC1 against PC2. These 13 features
are selected via stability selection.}
\centering \includegraphics[scale=0.5]{Figure5.pdf}\label{scvxpath}
\end{figure}


\begin{figure}[!htb]
\protect\caption{The left-panel graph shows the true cluster labels by plotting PC1
against PC2 using only 13 selected features. The right-panel graph
shows the estimated cluster membership using S-AMA at the selected
tuning parameters.}
\centering \includegraphics[width=16cm,height=8cm]{Figure6_colored.pdf}\label{truebestlabel}
\end{figure}


\section{Summary} \label{sec:summary}

In this paper, an extension of convex clustering, sparse convex clustering, is proposed to simultaneously cluster observations and conduct feature selection. Two algorithms, S-ADMM and S-AMA, are developed to implement the new method. The numerical results show that S-AMA is computationally faster and delivers better performance than S-ADMM. In addition, the numerical results show that the selection of tuning parameters in sparse convex clustering is important and the tuning method based on clustering stability performs well. Moreover, this work can motivate future work. \cite{Chi2016convex} presented a convex formulation of the biclustering problem, which seeks to cluster observations and features at the same time. Similarly, we can extend convex biclustering to sparse convex biclustering, in order to conduct biclustering and feature selection simultaneously.



\bibliographystyle{apa}

