PMSSC: Parallelizable multi-subset based self-expressive model for subspace clustering

In the past few years, there has been a surge of spectral clustering-based algorithms that segment a set of data points by performing spectral clustering. Previous classical methods, such as $k$-subspaces [15] and median $k$-flats [16], assume that the dimensionalities of the underlying subspaces are given in advance. This latent knowledge is generally hard to access in many real-world applications. In addition, these methods are usually non-convex and thus sensitive to initialization [17, 18]. Aiming to relax the limitations of the $k$-subspace algorithm, the majority of modern subspace clustering methods explored have turned to spectral clustering [19, 20], which segment data using an affinity matrix that captures whether a certain pair of data points lie on the same subspace. While many early methods [21, 22, 12, 23] achieve better segmentation than classical clustering algorithms even without the latent knowledge, these methods produce erroneous segmentation results for data points near the intersection of two subspaces due to the dense sampling of points lying on the subspace [24]. We now introduce previous subspace clustering approaches based on spectral clustering, then describe various techniques of scalable subspace clustering methods for dealing with large-scale datasets, which are closer to our proposed method.

Most subspace clustering approaches based on spectral clustering consist of two phases: (i) computing an affinity matrix based on the nonzero coefficients that appear in the representation of each data point as a combination of other points, and (ii) segmenting data points from the computed affinity matrix by applying spectral clustering. The key to the success of segmentation is the phase of computing the affinity matrix. Therefore, many methods have been proposed to compute the affinity matrix. For example, local subspace affinity [24] and spectral curvature clustering (SCC) [25] find neighborhoods based on the observation that a point and its $k$-nearest neighbors often lie on the same subspace. However, the computational complexity of finding multi-way similarity in these methods grows exponentially with the number of subspace dimensions, motivating the use of a sampling strategy to lower the computational complexity [9]. Recently, the self-expressive model, which employs the self-expressive property in Eq. (1), has become the most popular one. In particular, SSC takes advantage of sparsity [26] by adopting $\ell_{1}$ norm regularization of the coefficient vector to achieve high clustering performance. This idea has motivated many methods, using the $\ell_{2}$ norm in least squares regression [27], the nuclear norm in low rank representation (LRR) [28], the $\ell_{1}$ plus $\ell_{2}$ norm in elastic net subspace clustering (EnSC) [29], and the Frobenius norm in efficient dense subspace clustering [30]. In practice, however, solving the $\ell_{p}$ norm minimization problem for large-scale data may be prohibitive. Also, the memory required becomes larger as the amount of data increases.

When constructing the affinity matrix, several methods based on spectral clustering suffer from high computational complexity. To reduce the computational complexity of this phase, a sparse self-expressive model adopting a greedy algorithm was proposed in [31, 10]. However, these approaches lead to unsatisfactory clustering results if the nonzero elements do not contain sufficient connections within each optimized coefficient vector [32]. Other popular approaches to alleviate the computational and memory loads were inspired by a sampling strategy. In [13], scalable sparse subspace clustering (SSSC) is computationally efficient, using a subset generated by random sampling. However, because the random sampling method relies on a single subset, data points from the same subspace will not be represented by the self-expressive model if they are not appropriately sampled. Exemplar-based subspace clustering [33, 12] is an efficient sampling technique that iteratively selects the least well-represented point as a subset to address the problem. Selective sampling-based scalable sparse subspace clustering (S⁵C) [14], which generates a subset by selective sampling, provides approximation guarantees of the subspace-preserving property. In [34], the subspace-preserving representations are found by solving a consensus problem with multiple subsets to improve the connectivity of the affinity matrix. In [35], a divided-and-conquer framework using multiple subsets obtained by separating the entire dataset is proposed. While this approach can deal with large-scale data, final segmentation results depend on the self-expressive properties of the optimized self-expressive coefficient vectors of each subset. Our method differs significantly from [35] and [34] in that our proposed self-representation model is designed to minimize the difference from the original data points by combining the self-expressive property of multiple subsets. Lastly, in this paper, we limit our discussion to non-deep learning approaches, which are more mathematically straightforward to explain and rely less on parameter tuning.

As a problem definition, our final goal is to find the self-expressive coefficient vector $\boldsymbol{c}_{i}$, which satisfies the subspace-preserving representation in Eq. (1). That is, the self-expressive residual can be obtained by solving the following optimization problem,

where $\|\cdot\|_{0}$ is the $\ell_{0}$ pseudo-norm that returns the number of nonzero entries in the vector. This optimization problem has been shown [36, 37] to recover provably subspace-preserving solutions using the orthogonal matching pursuit (OMP) algorithm [38]. $s$ is a tuning parameter for the OMP algorithm, which controls the sparsity of the solution by selecting up to $s$ entries in the coefficient vector $\boldsymbol{c}_{i}$. Although the OMP algorithm is computationally efficient and is guaranteed to give subspace-preserving solutions under mild conditions, it is unable to produce a subspace-preserving solution with a number of nonzero entries exceeding the dimensionality of the subspace [10]. This leads to poor clustering performance with too sparse affinity between data points, especially when the density of data points lying on the subspace is low.

We propose a novel subspace clustering algorithm with a parallelizable multi-subset based self-expressive model, as illustrated in Fig. 1. Sec. III-B introduces our proposed self-expressive model that extends the model in Eq. (2) to multiple subsets via a sampling technique. Sec. III-C then explains the solution of our self-expressive model by the OMP algorithm. Finally, we summarize the proposed subspace clustering algorithm in Algorithm 4.

To deal with large-scale data, we first generate $T$ index subsets from the whole dataset by weighted random sampling [39] as follows:

where $\mathcal{I}^{(t)}$ is the index set of the $t$-th subset that is sampled with probability proportional to the elements of the weight vector $\boldsymbol{w}^{(t)}\in\mathbb{R}^{N}$, $[N]$ is $N$ indices $\{1,\ldots,N\}$, $0<\delta\leq 1$ is the sampling rate, and $n(\cdot)$ is the cardinality function that is a measure of the number of elements. The $t$-th selected element of $\boldsymbol{w}^{(t)}$ is updated as $\boldsymbol{w}_{i}^{(t+1)}=0.1\boldsymbol{w}_{i}^{(t)}$. Then, in each sampled $t$-th subset, the optimization problem in Eq. (2) can be expressed as follows:

where $X^{(t)}\in\mathbb{R}^{D\times N}$ is the data matrix of the randomly sampled $t$-th subset. $\boldsymbol{c}_{i}^{(t)}\in\mathbb{R}^{N}$ is the self-expressive coefficient vector for each data point $\boldsymbol{x}_{i}^{(t)}$ in the $t$-th subset. Note that to ensure the dimensionality of $\boldsymbol{c}_{i}^{(t)}$ is $N$, the columns of each data matrix $X^{(t)}$ corresponding to the non-sampled indices are replaced by zero-vectors: $\boldsymbol{x}_{i}^{(t)}=\boldsymbol{0},\ \forall i\notin\mathcal{I}^{(t)}$. From each optimized coefficient vector $\boldsymbol{c}_{i}^{\ast(t)}$, each data point $\boldsymbol{x}_{i}^{(t)}$ can be represented by a self-expressive model, given by:

where $\boldsymbol{y}_{i}^{(t)}$ is the data point computed by the self-expressive model from the $t$-th subset. In practice, however, the data point $\boldsymbol{y}_{i}^{(t)}$ in Eq. (5) generally has an error term $\boldsymbol{z}_{i}$, i.e., $\boldsymbol{y}_{i}^{(t)}=\boldsymbol{x}_{i}+\boldsymbol{z}_{i}$, because of the limitations of using $X^{(t)}$ as a dictionary for reconstruction. To minimize $\boldsymbol{z}_{i}$, we first represent $\boldsymbol{x}_{i}$ as a linear combination of $\boldsymbol{y}_{i}^{(t)}$, as follows:

where $\boldsymbol{b}^{(i)}\in\mathbb{R}^{T}$ is the weight coefficient vector to represent $\boldsymbol{x}_{i}$, and $b_{t}^{(i)}\in\mathbb{R}$ is the $t$-th entry of $\boldsymbol{b}^{(i)}$. The coefficient vector $\boldsymbol{b}^{(i)}$ of the linear combination in Eq. (6) can be obtained by solving the following optimization problem,

For simplicity, we introduce a data matrix $Y=[\boldsymbol{y}_{i}^{(1)},\ldots,\boldsymbol{y}_{i}^{(T)}]\in\mathbb{R}^{D\times T}$ with each data point $\boldsymbol{y}_{i}^{(t)}$ from Eq. (5) as columns, and rewrite Eq. (7) as

This is the formulation of the optimization problem for subspace clustering in Eq. (2), and can be further described as:

Unlike in Eq. (1), here $Y$ is the data matrix computed from each subset to represent $\boldsymbol{x}_{i}$. Thus, no constraint is required to avoid the trivial solution of representing a point as a linear combination of itself. To explicitly express Eq. (2), the self-expressive coefficient vector $\boldsymbol{c}_{i}^{\ast}$ corresponding to $X$ is obtained by

It is worth noting that each $\boldsymbol{c}^{\ast{(t)}}$ can be determined independently from each subset, so can be computed in parallel for speed.

In this section, we show that the parameters of the proposed PMS model can be determined by dividing the optimization problem into two small optimization problems as summarized in Algorithm 1. Overall, Eq. (10) can be determined by solving the minimization problems in Eqs. (4) and (8). We introduce both of the minimization procedures below. Specifically, initially, $T$ subset data matrices $\{X^{(t)}\}_{t=1}^{T}$ are generated based on the sampling set $\mathcal{I}^{(t)}$ in Eq. (3).

To efficiently solve for $\boldsymbol{c}_{i}^{\ast(t)}$ in Eq. (4), we introduce Algorithm 2 based on the OMP algorithm. The support set and the residual are initialized to $\mathcal{S}=\emptyset$ and $\boldsymbol{r}=\boldsymbol{x}_{i}^{(t)}$, respectively. $\mathcal{S}$ denotes the index set, which is updated on each iteration by adding one index $j^{\ast}$. $j^{\ast}$ is computed using

Then, using the updated $\mathcal{S}$, the self-expressive coefficient vector $\boldsymbol{c}_{i}^{\ast(t)}$ is found by solving the following problem:

where $\mathrm{supp}(\boldsymbol{c}_{i}^{(t)})$ is the support function that returns the subgroup of the domain containing the elements not mapped to zero. $\boldsymbol{r}$ is updated using:

This process is repeated until the number of iterations $k$ reaches its limit $s$ or $\boldsymbol{r}$ is smaller than the error $\epsilon$.

To find $\boldsymbol{b}^{\ast(i)}$, Eq. (8) can also be solved via the OMP algorithm, as shown in Algorithm 3. The input data matrix $Y\in\mathbb{R}^{D\times T}$ is generated by Eq. (5) from $\boldsymbol{c}_{i}^{\ast(t)}$. Note that the size of $Y$ depends on the number of subsets $T$ and is much smaller than the number of data points. The maximum number of repetitions is $T$, and $\mathcal{S}$ is updated by finding the index $j^{\ast}$ satisying

In addition, the weight coefficient vector $\boldsymbol{b}^{\ast(i)}$ and update of $\boldsymbol{r}$ are determined by solving

For clarity, we summarize the entire framework of our proposed subspace clustering approach in Algorithm 4, calling it the parallelizable multi-subset based sparse subspace clustering (PMSSC) method. Given $X$ and parameters $T$, $\delta$, $s$, and $\epsilon$, the optimal solution $C^{\ast}$ can be found using Algorithm 1. We thus define the affinity matrix as $A=|C^{\ast}|+|C^{\ast T}|$ using the computed $C^{\ast}$; the final clustering results can be obtained by applying spectral clustering to $A$ via normalized cut [19].

We compare our approach to the following eight methods: SCC [25], LRR [28], thresholding-based subspace clustering (TSC) [40], low rank subspace clustering (LRSC) [41], SSSC [13], EnSC [29], SSC-OMP [10], and S⁵C [14]. Tests for all comparative methods used provided source code, and each parameter was carefully tuned to give the best clustering accuracy. For spectral clustering, except for SCC, S⁵C, and SSSC, we applied normalized cut [19] to the affinity matrix $A=|C|+|C^{T}|$. (SCC and S⁵C have their own spectral clustering phase, while SSSC obtains clustering results from the data split into two parts). Unlike SSC-OMP, our method, which involves independent calculation for each subset, can be implemented in parallel with multi-core processing. All algorithms ran on an AMD Ryzen 7 3700x processor with 32 GB RAM. Following [10], as quantitative evaluation metrics, we evaluated each algorithm using clustering accuracy (acc: $a\%$), subspace-preserving representation error (sre: $e\%$), connectivity (conn: $c$), and runtime ($t$ seconds). Clustering accuracy represents the percentage of correctly labeled data points:

where $\pi$ is a permutation of the $L$ cluster groups. $Q^{\mathrm{est}}$ and $Q^{\mathrm{true}}$ are the estimated labeling result and ground-truth, respectively, scoring one in the $(i,j)$th entry if a data point $j$ belongs to the $i$-th cluster and zero otherwise. The subspace-preserving representation error indicates the average fraction of affinities formed from other subspaces in each $\boldsymbol{c}_{j}$,

where $\omega_{ij}\in\{0,1\}$ is the true affinity, and $\|\cdot\|_{1}$ returns the $\ell_{1}$ norm. The connectivity shows the average connection of the affinity matrix with $L$ cluster groups as follows:

where $\lambda_{2}$ is the second smallest eigenvalue of the normalized Laplacian for each of the $L$ subgraph, and $\lambda_{2}^{(i)}$ indicates the algebraic connectivity for each cluster. If $c=0$, then at least one subgraph is not connected [42].