Fast Fusion Clustering via Double Random Projection

Abstract

In unsupervised learning, clustering is a common starting point for data processing. The convex or concave fusion clustering method is a novel approach that is more stable and accurate than traditional methods such as k-means and hierarchical clustering. However, the optimization algorithm used with this method can be slowed down significantly by the complexity of the fusion penalty, which increases the computational burden. This paper introduces a random projection ADMM algorithm based on the Bernoulli distribution and develops a double random projection ADMM method for high-dimensional fusion clustering. These new approaches significantly outperform the classical ADMM algorithm due to their ability to significantly increase computational speed by reducing complexity and improving clustering accuracy by using multiple random projections under a new evaluation criterion. We also demonstrate the convergence of our new algorithm and test its performance on both simulated and real data examples.

1. Introduction

Clustering is a pivotal technique in unsupervised learning, applied extensively across various scientific and technological fields that handle large datasets. Clustering also plays a crucial role in data labelling, which sets the stage for the application of artificial intelligence and machine learning models [1,2] on the organized data to perform predictive analytics and classification tasks. Traditional clustering algorithms like k-means, Gaussian mixture models, and hierarchical clustering often face stability challenges due to their concave optimization formulations, which can lead to variability in results due to factors such as initial conditions or data outliers [3,4,5]. Recent advancements in convex or concave fusion methods have shown promise in enhancing stability, achieving more consistent global or local optimality and reliable estimation of cluster centers and counts through sparse-inducing penalties on pairwise centers [6,7,8,9]. For clustering high-dimensional data, the data can be mapped into a high-dimensional feature space (kernel space) for processing [10], or clustering can be achieved by optimizing a smooth and continuous objective function that is based on robust statistics [11]. This paper introduces a comprehensive empirical validation of these methods across simulation studies and real data analysis, detailing their improved stability over traditional methods and the practical implications of these advancements.

In fusion clustering, p-dimensional observations ${\mathit{X}}_i$, $i=1,\dots ,n$ are each parameterized by their own centroid ${\mathit{\mu }}_i$. These centroids are estimated under the assumption that all observations can be grouped into K clusters ${\mathcal{G}}_1,\dots ,{\mathcal{G}}_K$, such that for $i\in {\mathcal{G}}_k$, ${\mathit{\mu }}_i={\mathit{\rho }}_k$, where ${\mathit{\rho }}_k$ represents the cluster center for observations in cluster ${\mathcal{G}}_k$. Fusion clustering aims to concurrently estimate the cluster centroids ${\mathit{\rho }}_k$ and the partitions ${\mathcal{G}}_k$ by minimizing the following objectives

$$\frac{1}{2}\sum _{i=1}^n\parallel {\mathit{X}}_i-{\mathit{\mu }}_i{\parallel }^2+\sum _{i<j}{p}_{\lambda }(\parallel {\mathit{\mu }}_i-{\mathit{\mu }}_j{\parallel }_{\tau }).$$

(1)

The penalty function $p_{\lambda }(\parallel ·\parallel \tau )$ is used to control the complexity of the model, and it is determined by the tuning parameter $\lambda$. The form of the norm used is represented by $\parallel ·\parallel \tau$. This penalty function is typically used in fusion clustering to encourage sparsity in the estimated cluster centroids.

The penalty function $p_{\lambda }(\parallel ·\parallel \tau )$ controls the complexity of the model and is determined by the tuning parameter $\lambda$. The norm used is $\parallel ·\parallel \tau$. The penalty function is typically used in fusion clustering to promote sparsity in cluster centroids.

Convex fusion clustering methods have been widely studied due to their computational simplicity and ability to find global optima. These methods often employ ${\ell }_1$, ${\ell }_2$, or ${\ell }_{\infty }$ penalties as the penalty function $p_{\lambda }{(\parallel ·\parallel }_{\tau })$ [12,13,14,15,16,17]. However, convex fusion can lead to biased estimates of the individual centroids, resulting in solutions with a large number of dense clusters [18,19]. To address this issue, researchers have proposed using concave fusion clustering methods, such as those using minimax concave penalties (MCPs) [20], truncated Lasso penalties (TLPs) [8], and arbitrary concave penalties.

While robust, convex and concave fusion clustering methods are computationally demanding with a $\mathcal{O}\left(n^{2}p\right)$ complexity, which can limit their practicality in scenarios involving large sample sizes n and high-dimensional datasets p. This article proposes a strategy for overcoming this limitation using random projection techniques [21,22,23,24]. The approach involves the construction of a random diagonal matrix whose diagonal elements are sourced from a binary distribution. This matrix is then projected onto the pairwise component of the fusion method. By doing so, the number of pairwise differences between individual centroids, $\parallel {\mathit{\mu }}_i-{\mathit{\mu }}_j\parallel$, is substantially reduced. This reduction not only decreases the computational load but also maintains the integrity of the clustering process, enhancing the algorithm’s scalability without excessively increasing the operational overhead. We provide empirical evidence demonstrating that this method significantly reduces the computational time while preserving the clustering quality, as shown in our simulation section.

In unsupervised learning, rapid clustering processes are crucial for handling large datasets efficiently. Our study introduces a novel approach to fusion clustering to enhance computational speed without compromising accuracy. Our contributions are summarized as follows: (1) We propose using random projection techniques to simplify the fusion aspect of clustering, effectively diminishing the pairwise centroids discrepancies and significantly boosting computational efficiency by minimizing the fusion step’s complexity. (2) We have developed a novel double recursive random projection ADMM method designed for efficient high-dimensional fusion clustering, improving the accuracy of clustering.

In the remainder of this paper, the proposed new ADMM algorithm will be described in Section 2. This section will also include an analysis of the computational complexity and convergence of the algorithm. It will also include a strategy for improving cluster accuracy. The finite-sample properties of the proposed new ADMM algorithm will be evaluated through simulation studies in Section 3, and the method will be demonstrated using a real data example in Section 4. Concluding remarks will be presented in Section 5, and technical proofs will be provided in the Appendix A and Appendix B.

2. Methodology

To improve convex or concave fusion clustering efficiency, we propose an extension of the classical ADMM algorithm based on a random projection called RP-ADMM. A random projection can significantly reduce the time and computational resources needed to analyze high-dimensional data, making it suitable for large datasets and real-time processing. In this section, we will discuss the RP-ADMM algorithm’s computational complexity and convergence.

2.1. Random Projection Based ADMM

Previous ADMM algorithms for convex or concave fusion clustering [6,8] have suffered from a high computational burden due to the need to consider all $n(n-1)/2$ pairwise differences between individual centroids. This is represented by the fusion matrix $\mathcal{E}={\{({\mathit{e}}_i-{\mathit{e}}_j),i<j\}}_{\frac{n(n-1)}{2}\times n}^T$, where ${\mathit{e}}_i$ is the ith unit vector with a 1 in the ith position and 0s elsewhere, and ${\mathit{e}}_i-{\mathit{e}}_j$ can be interpreted as the difference between the ith and jth individual centroids. The computational complexity of this approach is $\mathcal{O}\left(n^2\right)$, which becomes infeasible for large sample sizes n.

• Bernoulli distribution-based random projections ADMM  

It is worth noting that pairwise differences between individual centroids can be deduced from other differences. For example, if we know that ${\mathit{\mu }}_1-{\mathit{\mu }}_2=\mathit{0}$ and ${\mathit{\mu }}_2-{\mathit{\mu }}_3=\mathit{0}$, we can conclude that ${\mathit{\mu }}_1-{\mathit{\mu }}_3=\mathit{0}$. This means that it may be unnecessary to consider the row ${\mathit{e}}_1-{\mathit{e}}_3$ in $\mathcal{E}$. To reduce the computational burden of convex or concave fusion clustering, we propose a random projection approach. This only considers a small subset of the $n(n-1)/2$ pairwise differences between individual centroids. This is achieved by generating indicators ${\pi }_{ij}$ from a Bernoulli distribution with probability $\alpha$. We then form a random matrix $\mathbf{\Pi }$, which is a diagonal matrix with diagonal elements ${({\pi }_{12},\cdots ,{\pi }_{1n},{\pi }_{23},\cdots ,{\pi }_{2n},\cdots ,{\pi }_{(n-1)n})}^T$. If ${\pi }_{ij}=1$, the difference between ${\mathit{\mu }}_i$ and ${\mathit{\mu }}_j$ is taken into account; if ${\pi }_{ij}=0$, it is not considered. The probability $\alpha$ controls the size of the subset of pairwise differences considered. The matrix $\mathbf{\Pi }\mathcal{E}$ can be seen as a projection of $\mathcal{E}$ onto a sparse matrix. This is with about $n(n-1)(1-\alpha )/2$ rows being zero vectors and about $n(n-1)\alpha /2$ ones being nonzero vectors. This projection is based on a Bernoulli distribution. Finally, we form a new fusion matrix $\mathbf{\Omega }$ by deleting the rows of zero vectors in $\mathbf{\Pi }\mathcal{E}$. The new fusion matrix is given by $\mathbf{\Omega }={({\mathbf{\Omega }}_1,\cdots ,{\mathbf{\Omega }}_{\kappa })}^T$, where ${\mathbf{\Omega }}_j,$$j=1,\cdots ,\kappa ,$ denotes jth row vector of $\mathbf{\Omega }$.

 We just consider $\tau =2$ in (1) for simplicity and propose a random projection-based fusion criterion by

$${\ell }_p(\mathit{\mu };\lambda )=\frac{1}{2}\sum _{i=1}^n\parallel {\mathit{X}}_i-{\mathit{\mu }}_i{\parallel }^2+\sum _{i<j}{\pi }_{ij}{p}_{\lambda }(\parallel {\mathit{\mu }}_i-{\mathit{\mu }}_j\parallel ),$$

(2)

where $\mathit{\mu }={({\mathit{\mu }}_1,\cdots ,{\mathit{\mu }}_n)}_{n\times p}^T$. Furthermore, the objective function in (2) is equivalent to

$$\begin{array}{cc}\hfill {\tilde{\ell }}_p(\mathit{\mu },\mathit{\varphi };\lambda )& =\frac{1}{2}{\parallel \mathit{X}-\mathit{\mu }\parallel }_F^2+\sum _{j=1}^{\kappa }{p}_{\lambda }\left(\right||{\mathit{\varphi }}_j\left|\right|),\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \mathbf{\Omega }\mathit{\mu }-\mathit{\varphi }=0,\hfill \end{array}$$

(3)

where $\mathit{X}={({\mathit{X}}_1,\cdots ,{\mathit{X}}_n)}^T$, $\mathit{\varphi }={({\mathit{\varphi }}_1,\cdots ,{\mathit{\varphi }}_{\kappa })}_{\kappa \times p}^T$. Under the constraints in (3), the augmented Lagrangian $Q(\mathit{\mu },\mathit{\varphi },\mathit{\eta };\lambda )$ has the form

$$\begin{array}{c}\hfill {\tilde{\ell }}_p(\mathit{\mu },\mathit{\varphi };\lambda )+\sum _{j=1}^{\kappa }{\mathit{\eta }}_j^T({\mathit{\mu }}^T{\mathbf{\Omega }}_j-{\mathit{\varphi }}_j)+\frac{\phi }{2}{\parallel \mathbf{\Omega }\mathit{\mu }-\mathit{\varphi }\parallel }_F^2,\end{array}$$

(4)

where the dual variables $\mathit{\eta }={({\mathit{\eta }}_1,\cdots ,{\mathit{\eta }}_{\kappa })}_{\kappa \times p}^T$ are Lagrange multipliers, and $\phi$ is a tuning parameter. Under the iterative value ${\mathit{\mu }}^{\left(m\right)}$ and ${\mathit{\eta }}^{\left(m\right)}$ at the mth step, we conduct the Bernoulli distribution-based random projection ADMM (RP-ADMM) iterative algorithm and compute the estimates of $(\mathit{\varphi },\mathit{\eta },\mathit{\mu })$ as follows:

$$\begin{array}{ccc}\hfill {\mathit{\varphi }}^{(m+1)}& =& arg\underset{\mathit{\varphi }}{min}L(\mathit{\varphi },{\mathit{\mu }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)};\lambda ),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\mathit{\eta }}^{(m+1)}& =& {\mathit{\eta }}^{\left(m\right)}+\phi (\mathbf{\Omega }{\mathit{\mu }}^{\left(m\right)}-{\mathit{\varphi }}^{(m+1)}),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill {\mathit{\mu }}^{(m+1)}& =& arg\underset{\mathit{\mu }}{min}Q(\mathit{\mu },{\mathit{\varphi }}^{(m+1)},{\mathit{\eta }}^{(m+1)};\lambda ),\hfill \end{array}$$

(7)

where $L(\mathit{\varphi },{\mathit{\mu }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)};\lambda )$ equals

$$\frac{\phi }{2}{\parallel \mathbf{\Omega }{\mathit{\mu }}^{\left(m\right)}-\mathit{\varphi }+{\phi }^{-1}{\mathit{\eta }}^{\left(m\right)}\parallel }_F^2+\sum _{j=1}^{\kappa }{p}_{\lambda }\left(∥{\mathit{\varphi }}_j∥\right),$$

(8)

and $Q(\mathit{\mu },{\mathit{\varphi }}^{(m+1)},{\mathit{\eta }}^{(m+1)};\lambda )$ equals

$$\begin{array}{cc}\hfill & {\tilde{\ell }}_p(\mathit{\mu },{\mathit{\varphi }}^{(m+1)})+\frac{\phi }{2}{\parallel \Omega \mathit{\mu }-{\mathit{\varphi }}^{(m+1)}\parallel }_F^2\hfill \\ \hfill & +\sum _{j=1}^{\kappa }{\mathit{\eta }}_j^{T(m+1)}({\mathit{\mu }}^T{\mathbf{\Omega }}_j-{\mathit{\varphi }}_j^{(m+1)}).\hfill \end{array}$$

(9)

Ma and Huang (2017) [18] have argued that under (8), the element ${\mathit{\varphi }}_j^{(m+1)}$ of ${\mathit{\varphi }}^{(m+1)}$ is the minimizer of $\frac{\phi }{2}\left|\right|{\mathit{\zeta }}_j^{\left(m\right)}-{\mathit{\varphi }}_j{\left|\right|}^2+p_{\lambda }\left(\right||{\mathit{\varphi }}_j\left|\right|)$, where ${\mathit{\zeta }}_j^{\left(m\right)}={\mathbf{\Omega }}_j^T{\mathit{\mu }}^{\left(m\right)}+{\phi }^{-1}{\mathit{\eta }}_j^{\left(m\right)}$. For different thresholding operator $p_{\lambda }(·)$, the estimate ${\mathit{\varphi }}_j^{(m+1)}$ has different results. Such as,

• For the Lasso penalty [25],

  $$\begin{array}{c}{\mathit{\varphi }}_j^{(m+1)}=S({\mathit{\zeta }}_j^{\left(m\right)},\lambda /\phi );\hfill \\ S(\mathit{w},t)=\left\{\begin{array}{cc}(1-t/|\left|\mathit{w}\right|\left|\right)\mathit{w},\hfill & \mathrm{if}t/\left|\right|\mathit{w}\left|\right|<1;\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$
• For SCAD penalty [26] with $a>1/\phi +1$,

  $$\begin{array}{c}{\mathit{\varphi }}_j^{(m+1)}=\left\{\begin{array}{cc}S({\mathit{\zeta }}_j^{\left(m\right)},\lambda /\phi ),\hfill & \mathrm{if}\left|\right|{\mathit{\zeta }}_j^{\left(m\right)}\left|\right|\le \lambda +\lambda /\phi ;\hfill \\ {\mathit{\zeta }}_j^{\left(m\right)},\hfill & \mathrm{if}\left|\right|{\mathit{\zeta }}_j^{\left(m\right)}\left|\right|>a\lambda ;\hfill \\ \frac{S({\mathit{\zeta }}_j^{\left(m\right)},a\lambda /\left((a-1)\phi \right))}{1-1/\left(\right(a-1\left)\phi \right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$
• For the MCP [27] with $a>1/\phi$,

  $$\begin{array}{c}{\mathit{\varphi }}_j^{(m+1)}=\left\{\begin{array}{cc}\frac{S({\mathit{\zeta }}_j^{\left(m\right)},\lambda /\phi )}{1-1/\left(a\phi \right)},\hfill & \mathrm{if}\left|\right|{\mathit{\zeta }}_j^{\left(m\right)}\left|\right|\le a\lambda ;\hfill \\ {\mathit{\zeta }}_j^{\left(m\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$
• For the TLP [8] with $a>1$,

  $$\begin{array}{c}{\mathit{\varphi }}_j^{(m+1)}=\left\{\begin{array}{cc}S({\mathit{\zeta }}_j^{\left(m\right)},\lambda /\phi ),\hfill & \mathrm{if}\left|\right|{\mathit{\zeta }}_j^{\left(m\right)}\left|\right|\le a\lambda ;\hfill \\ {\mathit{\zeta }}_j^{\left(m\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

Through some algebra, the problem of (9) is equivalent to the minimization of the function $h(\mathit{\mu },{\mathit{\varphi }}^{(m+1)},{\mathit{\eta }}^{(m+1)})$, which has the from

$$\begin{array}{c}\hfill \frac{1}{2}{\parallel \mathit{X}-\mathit{\mu }\parallel }_F^2+\frac{\phi }{2}\left|\right|\mathbf{\Omega }\mathit{\mu }-{\mathit{\varphi }}^{(m+1)}+{\phi }^{-1}{\mathit{\eta }}^{(m+1)}{\left|\right|}_F^2.\end{array}$$

Under the given value of ${\mathit{\varphi }}^{(m+1)}$, ${\mathit{\eta }}^{(m+1)}$, the updated ${\mathit{\mu }}^{(m+1)}$ are

$$\begin{array}{ccc}\hfill {\mathit{\mu }}^{(m+1)}& =& {(\phi {\mathbf{\Omega }}^T\mathbf{\Omega }+{\mathit{I}}_n)}^{-1}\left(\mathit{X}+\phi {\mathbf{\Omega }}^T({\mathit{\varphi }}^{(m+1)}-{\phi }^{-1}{\mathit{\eta }}^{(m+1)})\right)\hfill \end{array}$$

where ${\mathit{I}}_n$ is $n\times n$ identity matrix. ${\mathit{\mu }}^{(m+1)}$ and ${\mathit{\varphi }}^{(m+1)}$ are updated according to the random projection ADMM iterative algorithm (5)–(7) until the input of some convergence criteria, such as both dual and primal residuals being close to zero [28] in our practice. The convergence time of ADMM is highly related to the penalty parameter $\phi$. A poor selection of $\phi$ can result in a slow convergence for the ADMM algorithm [29] and thus RP-ADMM. In this paper, we fix $\phi =1$ throughout for simplicity.

To facilitate the updates of $({\mathit{\varphi }}^{(m+1)},{\mathit{\eta }}^{(m+1)},{\mathit{\mu }}^{(m+1)})$ at the $(m+1)$th step in (5) to (7) of the RP-ADMM iterative algorithm, we need to specify a proper initial value (warm start). Here, we set ${\mathit{\eta }}^{\left(0\right)}=\mathit{0}$, ${\mathit{\varphi }}^{\left(0\right)}=\mathbf{\Omega }{\mathit{\mu }}^{\left(0\right)}$ and obtain the initial estimators ${\mathit{\mu }}^{\left(0\right)}={({\lambda }^{★}{\mathbf{\Omega }}^T\mathbf{\Omega }+{\mathit{I}}_n)}^{-1}\mathit{X}$ as the minimizer of a ridge fusion criterion

$$\begin{array}{c}\hfill \frac{1}{2}{\parallel \mathit{X}-\mathit{\mu }\parallel }_F^2+\frac{{\lambda }^{★}}{2}{\parallel \mathbf{\Omega }\mathit{\mu }\parallel }^2.\end{array}$$

(10)

We summarize the above analysis in Algorithm 1.

Algorithm 1 RP-ADMM for fusion clustering

Input: data ${\mathit{X}}_1,\cdots ,{\mathit{X}}_n$; Initialize ${\mathit{\mu }}^{\left(0\right)}$, ${\mathit{\eta }}^{\left(0\right)}$; tuning parameter, $\lambda$ Output: an estimate of $\mathit{\mu }$  for $m=0,1,2,\cdots$do   compute ${\mathit{\varphi }}^{(m+1)}$ using (5)   compute ${\mathit{\eta }}^{(m+1)}$ using (6)   compute ${\mathit{\mu }}^{(m+1)}$ using (7)   if convergence criterion is met, then    Stop and denote the last iteration by $\widehat{\mathit{\mu }}\left(\lambda \right)$,   else    $m=m+1$   end if   end for

Practically, we would not want to conduct the RP-ADMM updates comprehensively until convergence to save computing time in the first iterations. Another trick is to adopt the initial values of subsequent convex relaxations as optimal values from the previous relaxed convex problem, which significantly reduces the number of RP-ADMM iterations.

2.2. Selection of Optimal Tuning Parameter

For a given $\lambda$, the converging value $\widehat{\mathit{\mu }}\left(\lambda \right)$ of the above RP-ADMM procedure is defined as

$$\widehat{\mathit{\mu }}\left(\lambda \right)={\mathrm{argmin}}_{\mathit{\mu }}{\ell }_p(\mathit{\mu };\lambda ),$$

(11)

where ${\ell }_p(\mathit{\mu };\lambda )$ is defined in (2) and the optimal value of $\lambda$ can be selected via a properly constructed data-driven criterion. In particular, we partition the support of $\lambda$ into a grid of ${\lambda }_{\min }={\lambda }_0<{\lambda }_1<\cdots <{\lambda }_J={\lambda }_{\max }$, and for each ${\lambda }_j$, we compute a solution path of $\widehat{\mathit{\mu }}\left({\lambda }_j\right)$ and obtain $\widehat{K}\left({\lambda }_j\right)$ distinct cluster centroids $\{{\widehat{\mathit{\rho }}}_1\left({\lambda }_j\right),\dots ,{\widehat{\mathit{\rho }}}_{\widehat{K}\left({\lambda }_j\right)}\left({\lambda }_j\right)\}$, The optimal $\widehat{\lambda }$ is selected by minimizing a data-driven BIC, i.e., $\widehat{\lambda }={\mathrm{argmin}}_{{\lambda }_j;j=1,\dots ,J}\mathrm{BIC}\left({\lambda }_j\right)$, where

$$\begin{array}{cc}\hfill \mathrm{BIC}\left(\lambda \right)& =log\left\{\frac{1}{np}{\parallel \mathit{X}-\widehat{\mathit{\mu }}\left(\lambda \right)\parallel }_F^2\right\}\hfill \\ \hfill & +(log\left(np\right)+2log\left(p\right))\widehat{K}\left(\lambda \right)/n.\hfill \end{array}$$

(12)

Subsequently, we obtain the estimator $\widehat{\mathit{\mu }}=\widehat{\mathit{\mu }}\left(\widehat{\lambda }\right)$, and the individuals can be separated into $\widehat{K}=\widehat{K}\left(\widehat{\lambda }\right)$ clusters accordingly, i.e., ${\widehat{\mathcal{G}}}_k=\{i:{\widehat{\mathit{\mu }}}_i={\widehat{\mathit{\rho }}}_k,i=1,\dots n\}$, $k=1,\dots ,\widehat{K}$.

Other methods for tuning parameters in clustering, such as generalized degrees of freedom with generalized cross-validation [8] and stability-based cross validation [25,30] can provide good results but may require extensive computation or the specification of a hyperparameter perturbation size [8]. In contrast, the proposed BIC is easy to compute and performs well in estimating cluster centroids and the true number of clusters (K). Figure 1 shows the change in BIC values against $log\left(\lambda \right)$ and the cluster number of the simulation. Across all cases with different values of n and p, we observe that BIC($\lambda$) decreases as the value of $log\left(\lambda \right)$ increases. With recovering the true cluster number $K=3$, BIC($\widehat{\lambda }$) reaches a minimum at the optimal $\widehat{\lambda }$. Moreover, when $log\left(\lambda \right)$ keeps increasing, the cluster centroids are continuously integrated, and BIC($\lambda$) is enlarged. However, further research is needed to fully prove the consistency of the BIC in combination with the objective function (2).

2.3. Recursive RP-ADMM and Cluster Matrix

In the above cluster analysis, the effect of randomness on the clustering results was not considered. However, empirical analysis has shown that the impact of this randomness on the estimated cluster centers and numbers is minimal (i.e., ${\widehat{\mathit{\rho }}}_k$’s and $\widehat{K}$’s). However, the impact on the final partitioning results (i.e., which observations are grouped into a single cluster) can be significant. In response to this, we propose the Recursive RP-ADMM (RRP-ADMM) procedure, which performs multiple RP-ADMM cluster analyses by generating M random matrices (i.e., ${\mathbf{\Omega }}_m$’s, $m=1,\cdots ,M$) and repeatedly conducting the analysis.

Once the multiple RP-ADMM cluster analyses have been completed, we must summarize the results. We define a $n\times n$ symmetric cluster matrix $\mathcal{C}$ where ${\mathcal{C}}_{ij}=1$ denotes that the ith and jth observations belong to the same cluster; otherwise, ${\mathcal{C}}_{ij}=0$. Another $n\times n$ symmetric matrix $\widehat{\mathcal{D}}$ is introduced, with element ${\widehat{\mathcal{D}}}_{ij}$ representing the relative frequency of the ith and jth observations belonging to the same cluster over the M independent RP-ADMM clustering procedures. The decision of whether the ith and jth observations should be grouped into a single cluster or not can then be treated as a classification problem, with the two possible class labels being 1 (belong to the same cluster) or 0 (do not belong to the same cluster). We can use an indicator function to transform the relative frequency into class labels and generate an estimator for the cluster matrix $\widehat{\mathcal{C}}$, i.e.,

$$\begin{array}{c}\hfill \widehat{\mathcal{C}}=\{{\widehat{\mathcal{C}}}_{ij}:{\widehat{\mathcal{C}}}_{ij}={\mathbb{1}}_{({\widehat{\mathcal{D}}}_{ij}\ge 0.5)}\},\end{array}$$

(13)

where ${\mathbb{1}}_{(·)}$ denotes the indicator function. We summarize the above procedure in Algorithm 2. This transformation can be understood as a voting-based aggregation strategy, similar to the one proposed by [31], which aims to reduce misclassification errors and improve the accuracy of the clustering. To evaluate the accuracy of the clustering results, we define a new measure called the similarity index (SI) between two data clusterings:

$$\begin{array}{c}\hfill \mathrm{SI}=\frac{1}{n^2-n}\parallel \widehat{\mathcal{C}}{-\mathcal{C}\parallel }_1=\frac{1}{n^2-n}\sum _{i=1}^n\sum _{j=1}^n|{\widehat{\mathcal{C}}}_{ij}-{\mathcal{C}}_{ij}|.\end{array}$$

(14)

Like the Rand Index (RI) measure [32], the newly introduced evaluation criterion can be seen as a measure of the percentage of correct decisions made by some algorithm. The SI values also range from 0 to 1, with lower values indicating better algorithm performance.

Algorithm 2 RRP-ADMM for fusion clustering

Input: data ${\mathit{X}}_1,\cdots ,{\mathit{X}}_n$; M; Initialize ${\mathit{\mu }}^{\left(0\right)}$, ${\mathit{\eta }}^{\left(0\right)}$; tuning parameter, $\lambda$ Output: an estimate of $\mathit{\mu }$   for$m=0,1,\cdots$, M do   compute ${\widehat{\mathit{\mu }}}^{\left(m\right)}$ using RP-ADMM  end for  while$1\le i\le n$ do   compute ${\widehat{\mathcal{D}}}_{ij}$ and ${\widehat{\mathcal{C}}}_{ij}$ from (13)  end while

The classical convex or concave fusion clustering procedure in (1) requires $\mathcal{O}\left(n^{2}p\right)$ operations and $\mathcal{O}(n^{2}p+np)$ of storage for a single round of ADMM updates with primal and dual residual calculations, because all pairs of centroids are shrunk together in this method.

The RP-ADMM algorithm significantly improves computational efficiency compared to classical ADMM algorithm. It requires only $\mathcal{O}(\kappa p+np)$ of storage, compared to $\mathcal{O}(n^{2}p+np)$ for the classical ADMM algorithm, because the variables $\mathit{\eta }$ and $\mathit{\varphi }$ have only $\kappa$ columns rather than $n(n-1)/2$. Additionally, the RP-ADMM algorithm requires only $\mathcal{O}\left(\kappa p\right)$ operations for its most computationally demanding step, in comparison to $\mathcal{O}\left(n^{2}p\right)$ for the classical ADMM algorithm. The RP-ADMM algorithm also requires $\mathcal{O}\left(\kappa n\right)$ operations to conduct Cholesky factorization in every iteration, in comparison to $\mathcal{O}\left(n^3\right)$ for the classical ADMM algorithm. This efficient Cholesky factorization is computed only once and reused across repeated RP-ADMM updates.

At the end of this subsection, we will demonstrate the convergence of the RP-ADMM algorithm by showing that the sequence generated by the algorithm contains a subsequence that converges to a stationary point.

Lemma 1.

Let ${\{{\mathit{\mu }}^{\left(m\right)},{\mathit{\varphi }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)}\}}_{k=1}^{\infty }$ be the sequence generated by Algorithm 1, then for some constant $c>0$,

$$\begin{array}{cc}\hfill & Q({\mathit{\mu }}^{(m+1)},{\mathit{\varphi }}^{(m+1)},{\mathit{\eta }}^{(m+1)})-Q({\mathit{\mu }}^{\left(m\right)},{\mathit{\varphi }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)})\hfill \\ & \le -\frac{c}{2}\parallel {\mathit{\mu }}^{(m+1)}-{\mathit{\mu }}^{\left(m\right)}{\parallel }^2+\psi {\parallel {\mathit{\eta }}^{(m+1)}-{\mathit{\eta }}^{\left(m\right)}\parallel }^2\hfill \end{array}$$

(15)

In order to prove that the sequence ${\{{\mathit{\mu }}^{\left(m\right)},{\mathit{\varphi }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)}\}}_{k=1}^{\infty }$ is convergent, we need to assume that ${\mathit{\varphi }}^{\left(m\right)}$ is bounded and $\psi \parallel {\mathit{\eta }}^{(m+1)-{\mathit{\eta }}^{\left(m\right)}}\parallel \to 0$ which are often observed in numerical tests.

Theorem 1.

If ${\left\{{\mathit{\varphi }}^{\left(m\right)}\right\}}_{k=1}^{\infty }$ are bounded and ${\psi }_2\parallel {\mathit{\nu }}^{(m+1)}-{\mathit{\nu }}^{\left(m\right)}{\parallel }_F+{\psi }_1{\parallel {\mathit{\eta }}^{(m+1)}-{\mathit{\eta }}^{\left(m\right)}\parallel }_F\to 0$, then ${\{{\mathit{\mu }}^{\left(m\right)},{\mathit{\varphi }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)}\}}_{k=1}^{\infty }$ is bounded. Moreover, there exist a subsequence ${\{{\mathit{\mu }}^{\left(k_j\right)},{\mathit{\varphi }}^{\left(k_j\right)},{\mathit{\eta }}^{\left(k_j\right)}\}}_{k_j=1}^{\infty }$, such that

$$\begin{array}{cc}\hfill \underset{k_j\to \infty }{lim}(\parallel {\mathit{\mu }}^{(k_j+1)}-{\mathit{\mu }}^{\left(k_j\right)}\parallel & +\parallel {\mathit{\varphi }}^{(k_j+1)}-{\mathit{\varphi }}^{\left(k_j\right)}\parallel \hfill \\ & +\parallel {\mathit{\eta }}^{(k_j+1)}-{\mathit{\eta }}^{\left(k_j\right)}\parallel )=0,\hfill \end{array}$$

and thus, ${\{{\mathit{\mu }}^{\left(m\right)},{\mathit{\varphi }}^{\left(m\right)},{\mathit{\eta }}^{\left(m\right)}\}}_{k=1}^{\infty }$ has a subsequence which converges to the stationary point.

3. Simulation

In this part of the study, simulation experiments were conducted to compare the performance of the extended and classical ADMM clustering algorithms in terms of computational time and clustering accuracy, using the evaluation criterion in (14). The Lasso-based fusion method often leads to the formation of dense clusters with a minor penalty for small differences in $\parallel {\mathit{\varphi }}_j\parallel$, which can result in the formation of many spurious clusters with very small differences among them [6]. In contrast, the concave penalty method tends to produce a clear cluster structure and a well-defined number of clusters [8]. Therefore, in this study, we focus on the MCP-based fusion method [27] which compares the conventional ADMM’s clustering performance and the proposed new ADMM algorithm.

3.1. Low-Dimensional Setting

In this part, we evaluated the clustering performance of the classical ADMM, RP-ADMM, and RRP-ADMM algorithms on low-dimensional synthetic data generated from three overlapping convex clusters with the same spherical shape in some number of dimensions p and sample size n. The synthetic data were generated from three populations ${\mathcal{P}}_k=\mathcal{N}({\rho }_k,\mathbf{\Sigma })$, $k=1,\cdots ,K$ with $K=3$, ${\rho }_1={\mathit{3}}_p$, ${\rho }_2={\mathit{0}}_p$, ${\rho }_3=-{\mathit{3}}_p$ and $\mathbf{\Sigma }={\left({\sigma }_{kj}\right)}_{p\times p}$ with ${\sigma }_{jj}=1$ and ${\sigma }_{kj}=0.1^{|k-j|}$ for $k\ne j$. This setting was chosen deliberately to allow overlap in the sample sets generated from clusters proximal to each other, thereby increasing the complexity of the clustering task. As illustrated in Figure 2c, the clustering performance using a single random projection (RP-ADMM) was suboptimal, indicating challenges with cluster separability under this setup. Conversely, Figure 2b demonstrates that recursive random projection (RRP-ADMM) significantly improved clustering results. The recursive times for the RP-ADMM and RRP-ADMM algorithms were set to $M=10$.

To evaluate the accuracy of the RP-ADMM, relax-and-split approach [33] (RS-ADMM) and RRP-ADMM algorithms in recovering the true cluster matrix, we generated a random sample of $n=60$ observations with 1–20 drawn from ${\mathcal{P}}_1$, 21–40 drawn from ${\mathcal{P}}_2$, and 41–60 drawn from ${\mathcal{P}}_3$, and set the number of dimensions to $p=5$. The probability $\alpha$ of generating a 1 in the random matrix was set to $\alpha =c\frac{log\left(n\right)}{n}$, where c controls the probability size. The level plots in Figure 2 use colour to visualize the values of 1’s and 0’s in the cluster matrix. The results show that both RP-ADMM and RRP-ADMM can accurately recover the true cluster matrix, with RRP-ADMM showing more accurate gradation than the true cluster matrix. Single random projection (RP-ADMM) can cause high variance in clustering outcomes due to the randomness of the sampling process. To mitigate this issue, we have adopted the voting-based pooling technique [31], which reduces variance by averaging results from recursive random projection (RRP-ADMM).

To further evaluate the performance of the algorithms, we calculated the values of the index SI defined in (14) after 100 replicates under different c choices. We depicted the results as boxplots in Figure 3. These results show that RRP-ADMM consistently improves clustering accuracy compared to RP-ADMM, as evidenced by the smaller median and standard error of SI values.

Next, we will compare the performance of classical ADMM and RRP-ADMM in terms of computation time per iteration and the SI after 100 trials. The sample size is varied with $n=60,150,240,360$ points and $\alpha =4\frac{log\left(n\right)}{n}$, while $p=2$ is kept constant. In this study, we have limited the number of points to 360, as the classical ADMM algorithm requires a significant amount of computation time for a single realization with more points. We will also compare the performance of the Similarity Index (SI) and Rand Index (RI) in evaluating the clustering results. Therefore, we should calculate the partitioning structure of all points based on the estimated cluster matrix graph. This process involves first identifying the point $a_1$ with the most neighbors and aggregating the connected points with point $a_1$ as cluster 1, then finding the second point $a_2$ with the most edges to form cluster 2, and repeating this process until there are no more points remaining.

Table 1. The mean values of Similarity index (SI), Rand Index (RI) and run time in seconds against different sample sizes and different methods after 100 replicates.

	ADMM				RRP-ADMM				RS-ADMM
Sample Size	SI	RI	Time		SI	RI	Time		SI	RI	Time
$n=60$	0.081	0.921	7		0.059	0.933	2		0.080	0.925	10
$n=150$	0.058	0.945	88		0.046	0.957	7		0.056	0.947	121
$n=240$	0.049	0.962	352		0.045	0.974	17		0.047	0.966	551
$n=360$	0.042	0.973	1582		0.040	0.986	41		0.042	0.978	1864

Note: ‘SI’ represents the similarity index defined in (14), ‘RI’ denotes Rand Index [32]. TIME is the required time in seconds in a single round of ADMM.

shows the mean values of the SI, RI, and the consumed time in seconds for different sample sizes under different methods after 100 replicates. Based on the data in Table 1, we can observe the following: (i) The proposed RRP-ADMM significantly reduces the time required for convex or concave fusion clustering, especially when the sample size increases. (ii) RRP-ADMM produces smaller SI and larger RI values, possibly due to the voting-based pooling technique improving cluster accuracy. (iii) As the sample size increases, the SI and RI values decrease. The boxplots in Figure 4 and Figure 5 demonstrate the superiority of the RRP-ADMM algorithm over the classical ADMM algorithm in terms of both the SI values and the square root of run time, as seen in the results obtained from 100 replicates with four different sample sizes. These results further reinforce our belief in the effectiveness of the RRP-ADMM algorithm.

3.2. High-Dimensional Setting

In this part, we investigate using the double random projection-based alternating direction method of multiplier (DRP-ADMM and DRRP-ADMM) algorithms for clustering high-dimensional data sets. We employ a recursive Gaussian distribution-based random projection strategy in the first step to mitigate the impact of randomness on cluster results. Since the classical ADMM algorithm is computationally intensive in high-dimensional settings, we focus on evaluating the performance of the DRP-ADMM and DRRP-ADMM algorithms with recursive times $M=9$, using three Gaussian random projections in the outer layer and three binary random projections in the inner layer. The simulated data sets consist of two overlapping convex clusters with the same spherical shape. They are generated using a population ${\mathcal{P}}_k=\mathcal{N}({\rho }_k,\mathbf{\Sigma })$, $k=1,2$ with ${\rho }_1={\mathit{1}}_p$, ${\rho }_2=-{\mathit{1}}_p$. Furthermore, $\mathbf{\Sigma }={\left({\sigma }_{kj}\right)}_{p\times p}$ with ${\sigma }_{jj}=1$ and ${\sigma }_{kj}=0.1^{|k-j|}$ for $k\ne j$. We consider four high-dimensional cases with $p=1000,2000,3000,5000$ and a fixed sample size of $n=100$.

We evaluate the accuracy of the DRP-ADMM and DRRP-ADMM algorithms in recovering the true cluster matrix. To do this, we first generate a Gaussian random matrix $\mathit{R}$ with dimensions $p\times q$ in the first projection. The elements of $\mathit{R}$ correspond to $\mathcal{N}(0,1/\sqrt{q})$. We set $q=\lceil \frac{\kappa }{{\epsilon }^2/2-{\epsilon }^3/3}log\left(n\right)\rceil$ with $\epsilon =1$ and $\kappa =\frac{5}{6}$. See [21,23] for the number of projections. In the second step, we generate a diagonal binary random matrix with probability $\alpha =4\frac{log\left(n\right)}{n}$ of equaling one. Then, we calculate the values of the SI index defined in Equation (14) and plot the results as boxplots in Figure 6 after 100 replicates for different values of p. The results show that the DRRP-ADMM algorithm consistently outperforms the DRP-ADMM algorithm regarding the median and standard error of the SI values for all values of p, indicating that the DRRP-ADMM algorithm improves clustering accuracy.

4. Real Data Analysis

In this study, we use the DrivFace dataset to demonstrate the effectiveness of our proposed clustering procedure. The DrivFace database consists of $n=606$ images of 640,480 pixels each, captured from four drivers (two women and two men) over different days and containing $p=17$ facial features such as glasses and beards. Each driver’s images containing similar facial features can be grouped into one cluster, resulting in a total of $K=4$ clusters as shown in Figure 7a. Firstly, we know the true labels of the dataset; that is, there are four clusters, and we also know which observations belong to the common cluster. Secondly, because the similarity among observations in the pictures is very high across different clusters, it is challenging to separate them. Therefore, we can use this dataset to evaluate our proposed clustering method.

Due to the large sample size of the DrivFace dataset, we do not use the classical ADMM algorithm, which would require $606\times (606-1)\times 17/2$ operations in a single ADMM iteration. Instead, we first scale the samples by each feature and apply the RP-ADMM procedure to estimate individual centers using a grid of $\lambda$ values. We plot the $fusiongrams$ of four selected variables in Figure 8, and the scrutiny of Figure 8a implies that some outlying points (influential points) cause the clusters to be dense. We then remove these 55 points and plot a new $fusiongram$ in Figure 8b. The optimal $\lambda$ value, as determined by the developed BIC criterion in Equation (12), is $1.38$, indicating that the estimated number of clusters is four, the same as the number of drivers. We apply the proposed RRP-ADMM algorithm with a Bernoulli-distribution-based random projection procedure to further improve the cluster accuracy using $\alpha =10log\left(n\right)/n$ and a recursive number $M=20$. Using the estimated optimal tuning parameter of $1.38$, we obtain the estimated cluster matrix in Figure 7b, which closely resembles the true cluster matrix in Figure 7a. The calculated similarity index (SI) value is $0.098$. Moreover, the value of Adjusted Rand Index (ARI) is 0.672.

5. Conclusions

We propose using the recursive random projection-based ADMM (RRP-ADMM) method to improve the speed and accuracy of convex and nonconvex fusion clustering. In simulations and real data examples, the RRP-ADMM method demonstrates superior performance in fast calculation and accurate clustering results. The RRP-ADMM algorithm is scalable and can be applied to deal with heterogeneous issues in any setting that involves fusion techniques.

However, some challenges still need to be addressed in this field. One challenge is efficiently transforming the cluster matrix graph into the target partitioning structure and determining the optimal number of clusters. Another challenge is using prior information about which points are more likely to be integrated into a single cluster to reduce the number of pairwise comparisons. Additionally, a further study is needed to determine the theoretical probability of achieving a probability of one in binary random projection. Another future research direction involves performing clustering simultaneously with feature selection, using techniques such as incorporating feature weights [34] or introducing sparsity [14].