Beyond Linear Subspace Clustering: A Comparative Study of Nonlinear Manifold Clustering Algorithms

Constructing a pairwise affinity matrix that reflects the multiple subspace structure of the data points is the main challenge in the linear SC algorithms based on spectral clustering. The affinity matrix is often constructed by exploiting the self-expressiveness property of the data points. This property, which is known as collaborative representation in the sparse representation literature [141], is based on the fact that each data point can be represented as a linear combination of other data points in the same subspace. Mathematically, for all $i\in\{1,\dots,n\}$,

where $C\in\mathbb{R}^{n\times n}$ is the coefficient matrix, $C(i,i)$ is the $i$-th diagonal entry of $C$, and $X(:,i)$ is the $i$-th column of $X$. The condition $C(i,i)=0\ \text{for all }i\in\{1,\dots.,n\}$ eliminates the trivial solution of expressing each point by itself. However, there are typically infinitely many solutions to (1) and many of them might not be subspace preserving, that is, all data point might not be expressed using a linear combination of data points belonging only to the same subspace. Let us formally define this concept.

To fulfill the goal of linear SC and reveal the association of each data point to the underlying subspace, the coefficient matrix should be (approximately) subspace preserving. Hence, to ensure that the obtained coefficient matrix is subspace preserving, several regularizations on the coefficient matrix are used in the literature, including sparsity [21] and low-rankness [59]. In general, the linear SC problem based on self-expressiveness is formulated as follows:

where $f$ is a regularization function for the coefficient matrix (see below), $E\in\mathbb{R}^{d\times n}$ models the noise, $g$ is a regularization function for modeling the noise (for example the $\ell_{1}$ norm for sparse gross noise, or the Frobenius norm for Gaussian noise), and $\operatorname{diag}(C)$ is the vector containing the diagonal entries of $C$. The conditions under which the obtained coefficient matrix is subspace preserving strongly depends on the regularization function $f$. Table 1 summarizes the most common regularization functions and the theoretical conditions for the corresponding coefficient matrix to be subspace preserving. In this table $||\cdot||_{0}$, $||\cdot||_{1}$, $||\cdot||_{*}$ and $||\cdot||_{F}^{2}$ indicate the $\ell_{0}$ norm (the number of nonzero entries), the component-wise $\ell_{1}$ norm (the sum of the absolute value of the matrix entries), the nuclear norm (the sum of singular values), and the squared Frobenius norm (the sum of the squares of the matrix entries), respectively.

In order to enhance the quality of the coefficient matrix $C$, several post-processing strategies exist. The goal of these approaches is to decrease the number of wrong between-subspaces connections and/or strengthening the within-subspaces connections. Even though there is no generally accepted strategy in the literature, a few common post-processing methods are as follows:

1. 1.

   Normalizing the columns of the coefficient matrix as $C(:,j)\leftarrow C(:,j)/||C(:,j)||_{\infty}$ where $||\cdot||_{\infty}$ is the infinity norm [21]. This can be helpful when some data points have drastically different norms compared to other data points.

2. 2.

   Applying a hard thresholding operator on each column of the coefficient matrix by keeping only the $k$ largest entries in absolute value [84]. This post-processing is based on the property of Intra-subspace Projection Dominance (IPD) which states that the entries corresponding to data points in the same subspace are larger (in absolute value) than the ones corresponding to data points in different subspaces.

3. 3.

   Performing an ad-hoc post-processing strategy using multiple steps. First, a percentage of the top entries in each column of the coefficient matrix are preserved. Next the shape interaction matrix method [15] is applied to remove the noise; it consists in a low-rank approximation of $C$ based on the skinny singular value decomposition. Finally each entry is raised to some power larger than one (the value 4 is the often used) in order to intensify the dominant connections. This strategy depends on several parameters with no explicit theoretical justification. However, it is widely used in nonlinear SC approaches based on neural networks [41, 42].

4. 4.

   Finding good neighbors which correspond to the key connections in the coefficient matrix [125]. This approach is based on three parameters: $k>q>p$. First, only the $k$ largest coefficients are kept for each data point. Then, only a subset of these $k$ connections is preserved based on the notion of good neighbors. A good neighbor is defined as a data point with at least $p$ (usually $p=1$) common data points inside the $k$ largest connections. The $q$ good neighbors of each data point with the largest coefficients form the final selected subset. A disadvantage of this approach is that it depends on parameters which are difficult to tune. (Note that using $q=k$ and $p=0$ reduces to the second post-processing approach described above.)

The methods in this category exploit the local structure of the data to construct the matrix $C$. To do so, a pairwise local similarity matrix $S\in\mathbb{R}^{n\times n}$ is first constructed from the data points. The matrix $S$ is symmetric, and its $(i,j)$-th entry measures how similar or close the data points $X(:,i)$ and $X(:,j)$ are in the ambient space. The entries of $S$ are computed in different ways, including the following:

1. 1.

   Binary k-nearest neighbors (K-NN): $S(i,j)$ is equal to one if $i$ is among the $k$ nearest neighbor of $j$, or vice versa [37], that is,

   $$S(i,j)=\begin{cases}1\qquad\text{if}\ X(:,j)\in\mathcal{N}_{k}(X(:,i))\ \text{ or}\ X(:,i)\in\mathcal{N}_{k}(X(:,j)),\\ 0\qquad\text{otherwise},\end{cases}$$

   where $\mathcal{N}_{k}(x)$ computes the $k$ nearest samples of the data point $x$.

2. 2.

   Continuous similarity functions: the value $S(i,j)$ is computed using various similarity functions, such as Gaussain kernels [129, 67, 139] given by

   $$S(i,j)=e^{\frac{-||X(:,i)-X(:,j)||_{2}^{2}}{\sigma^{2}}},$$

   where $\sigma$ is a parameter. Other data dependent similarity functions such as the geodesic distance for data on Grassmann manifolds can also be used [112].

3. 3.

   Weighted K-NN is a combiantion of the two strategies described above [60]:

   $$S(i,j)=\begin{cases}\text{sim}(X(:,i),X(:,j))\qquad X(:,j)\in\mathcal{N}_{k}(X(:,i))\ \text{ or}\ X(:,i)\in\mathcal{N}_{k}(X(:,j)),\\ 0\qquad\text{otherwise},\end{cases}$$

   where $\text{sim}\big{(}x,y\big{)}$ is a similarity function, such as Gaussian kernels, between two input vectors $x$ and $y$.

For a recent comprehensive overview of similarity and neighborhood construction algorithms, we refer the interested reader to [85].

Once the matrix $S$ is constructed, the following regularizer for $C$ is constructed using this local similarity information:

Let the Laplacian matrix $L\in\mathbb{R}^{n\times n}$ corresponding to $S$ be defined as $L=D-S$ where $D$ is the diagonal matrix with $D_{ii}=\sum_{j}S(i,j)$ for all $i$. The function (3) can be rewritten as:

where $\operatorname{tr}(\cdot)$ denotes the trace of a matrix. This regularizer $\mathcal{L}(C)$ is known as the Graph Laplacian or manifold regularization. It is used within existing global linear SC optimization problems, that is, it is added in the objective function of the optimization problem (2) with a proper penalty parameter. Manifold regularization promotes a grouping effect, which is defined as follows.

Intuitively, the grouping effect encourages locally similar points (that is, pairs of points whose corresponding value in $S$ is large) to have similar coefficient representations.

The grouping effect of Graph Laplacian regularizers generates graphs that better represent the structure of the data and hence are better connected [37]. This might benefit SC approaches based on self-expressiveness as they use spectral clustering for the final step [110]. Note that the Graph Laplacian regularization is based on the assumption that the representation coefficient vector corresponding to each sample, that is, $C(:,i)$ for each $i=1,...,n$, is changing smoothly on each manifold [7]. Hence, nearby data points in the ambient space have similar coefficient vectors. This regularization also implicitly assumes that the data points from different clusters are not likely to be near each other and that sufficiently many data points are sampled from each manifold.

Graph Laplacian regularization can be used exclusively with the reconstruction error, $||X-XC||_{F}^{2}$, as in the Smooth Representation (SMR) algorithm [37]. SMR solves the following optimizatin problem:

Note that SMR does not enforce $\operatorname{diag}(C)=0$ explicitly, because the regularization $\operatorname{tr}(CLC^{\top})$ prevents the identity matrix to be an optimal solution, for $\lambda$ sufficiently large, since the identity matrix $C=I$ satisfies $\|C(:,i)-C(:,j)\|_{2}=1$ for all $i\neq j$.

It is proven that Graph Laplacian regularization of the coefficient matrix leads to subspace preserving coefficients for noiseless linear independent subspaces [37]. Graph Laplacian regularization can be used in combination with other coefficient regularizations such as the $\ell_{1}$ norm [129], or the nuclear norm [60, 131]. In particular, using the $\ell_{1}$ norm along with Graph Laplacian regularization leads to the algorithm referred to as Laplacian Regularized $\ell_{1}$-SSC (LR$\ell_{1}$-SSC) which solves the following optimization problem [129]:

However, as we will see in the numerical experiments in Section 5, promoting the smoothness of the representation over the manifold might not be sufficient to recover the complex nonlinear structures. Moreover, these methods are rather sensitive to the choice of the similarity function used to construct the pairwise local similarity matrix $S$.

Instead of encouraging locally nearby points to have similar coefficient representations, as presented in the previous section, one could encourage the entries of $C$ corresponding to faraway points to have small/zero values. In other words, data points should not use faraway points in their representation, which we refer to as avoiding cannot-links. The name “cannot-links” is borrowed from the graph learning and spectral clustering literature to indicate the pairwise constraints on the links that are encouraged/enforced to be avoided [68, 119].

Avoiding cannot-links attracted less attention compared to the Graph Laplacian regularization methods. Generally, the cannot-links refer to graph links/coefficients corresponding to distant points that are specified based on a dissimilarity criterion. There are two strategies to integrate pairwise cannot-links constraints:

1. 1.

   Explicit elimination of faraway samples: A simple approach to prevent cannot-links is to explicitly add them as proper constraints in the optimization problem; for example, in [155], Zhuang et al. use the following constraint: for all $j$,

   $$C(i,j)=0\;\quad\text{ for all }X(:,i)\notin\mathcal{N}_{k}\big{(}X(:,j)\big{)}.$$

   A similar approach was proposed for linear SSC [31]. Instead of additional explicit constraints on the entries of the coefficient matrix, the self-expressiveness term is modified such that each sample is represented by the locally nearby data points.

2. 2.

   Implicit penalizing of distant pairwise connections: Cannot-links information can be added as a weighted penalty term. For example, in [148], the following regularization term was proposed:

   $$\sum_{i,j}||X(:,i)-X(:,j)||^{2}_{2}\ |C(i,j)|.$$

   This term is similar to (3). However, for any two faraway data points $X(:,i)$ and $X(:,j)$, the corresponding $C(i,j)$ is encouraged to be small while this is not the case for graph Laplacian regularization (for which $S(i,j)$ will be close to zero). In other words, in graph Laplacian regularization, the nearby samples are encouraged to have similar coefficient representation but there is no explicit penalty on the coefficient representations corresponding to faraway data points. In fact, graph Laplacian regularization tend to produce dense coefficient matrices (see the numerical results in Section 5.1), hence do not have the disadvantage of oversegmentation that avoiding cannot-links may have [2].

   A similar idea was presented in [149], together with non-negativity constraint on the coefficient matrix $C$. Another penalty term to penalize “non-local” connectivities was proposed in [12], namely $\sum_{i,j}e^{\operatorname{dist}(X(:,i),X(:,j))^{2}}C_{i,j}^{2}$, where $\operatorname{dist}$ is a normalized distance metric.

Avoiding non-local data points in self-expression might be intuitively similar to encouraging spatially close connections, as with the Graph Laplacian regularization. However, they lack the grouping effect of Graph Laplacian regularizations.

The approaches in this category rely on the estimation of the local tangent space for each data point, by fitting an affine subspace in the neighborhood of each data point. A simple strategy to approximate such tangent spaces is to consider a fixed predefined number of nearest neighbors for each sample, and approximate the tangent space by fitting an affine subspace to the data point and its neighbors. The crucial question is how to define the neighborhood of each data point in order to balance two goals: large enough to capture the local geometry of the manifolds, and small enough to avoid data points from other manifolds/clusters so as to preserve the curvature information of the underlying manifold.

Despite the fact that choosing a fixed value for neighborhood size with no prior knowledge is challenging, Deng et al. [19] proposed a multi-manifold embedding approach based on the angles between approximated tangent spaces using a fixed neighborhood size. In particular, they assumed that if two data points belong to the same manifold, the angle between their corresponding tangent space is small whereas the angle is large if they belong to different manifolds. In contrast to estimating tangent space using fixed neighborhood size, a different strategy is to allow the size of the neighborhood to be arbitrary. This can be achieved via sparse representations in order to estimate the neighboring samples and simultaneously calculate the contribution weight of each neighborhood sample. These weights would be helpful to reflect the clustering information as well. In this section, we focus on this particularly effective strategy.

Patel et al. [79] proposed a kernelized SC approach by using the kernel trick in the (affine) sparse SC formulation. Their algorithm, dubbed Kernelized Sparse Subspace Clustering (KSSC), solves the following optimization problem

where $\Phi(X)=[\phi(X(:,1)),\phi(X(:,2)),\dots,\phi(X(:,n))]$. The kernel trick can be used to avoid computing $\Phi(X)$ explicitly, using

where the kernel $\mathcal{K}_{XX}$ is given in (7). As for SSC, the authors used sparse regularization for the coefficient matrix. Finally, (8) can be reformulated as follows

In order to reduce the computational cost of calculating the coefficient matrix for high-dimensional data, the same authors [78] extended KSSC to combine dimensionality reduction and kernelized SC. To this end, their approach learns the implicit mapping of the data onto a high-dimensional feature space using the kernel method, and then projecting onto a low-dimensional space via a linear projection.

The ease of use of the kernel trick in the self-expressive based linear SC initiated the development of other nonlinear SC algorithms with additional regularizations, motivated by different purposes. These algorithms include for example the kernelized multi-view SC approaches from [138, 121]. Recently Kang et al. [48] proposed additional structure preserving regularization for KSSC in order to minimize the inconsistency between inner products of the transferred data, $\Phi(X)$, and inner products for reconstructed transferred data, $\Phi(X)C$.

However, the traditional challenges of kernel methods are still present for kernel based SC. The two main challenges are:

1. 1.

   Using the Frobenius norm to compute the representation error facilitates the use of the kernel trick; see (9). However, the Frobenius norm implicitly assumes that the noise follows a Gaussian distribution. How can the kernel trick be applied on non-Gaussian noise models including, e.g., gross corruptions and outliers?

2. 2.

   With no prior knowledge, how can we select an appropriate kernel function and its parameters?

In the rest of this section, we review the few approaches that tried to tackle these challenging limitations.

Even though kernel methods are considered principled approaches to provide non-linearity in linear models, they rely on selecting and tuning a kernel. In many (unsupervised) applications, with no prior knowledge, this is a challenging problem. Multiple Kernel Learning (MKL) methods overshadow this issue, by learning a consensus kernel $\mathcal{K}$ from a set of predefined candidate kernels $\{\mathcal{K}^{(i)}_{G}\}_{i=1}^{k}$, where $k$ is the number of given kernels. The sought kernel $\mathcal{K}$ is regularized towards the predefined candidate kernels via a penalty function $h\left(\mathcal{K},\{\mathcal{K}_{G}^{(i)}\}_{i=1}^{k},w\right)$ which is added as a penalty term in the objective function. The vector $w\in\mathbb{R}^{k}_{+}$ contains the weights/importance of the predefined kernels. This regularization function is typically defined in two ways:

1. 1.

   Centroid-based MKL [49] encourages the consensus kernel $\mathcal{K}$ to be close to every predefined kernel by minimizing

   $\displaystyle h\left(\mathcal{K},\{\mathcal{K}_{G}^{(i)}\}_{i=1}^{k},w\right)=\sum_{i=1}^{k}w(i)\left\|\mathcal{K}-K^{(i)}_{G}\right\|_{F}^{2}.$

   (15)

2. 2.

   Linear combination MKL [89] encourages the consensus kernel to be approximately a convex combination of the predefined kernels:

   $\displaystyle h\left(\mathcal{K},\{\mathcal{K}_{G}^{(i)}\}_{i=1}^{k},w\right)=\left\|\mathcal{K}-\sum_{i=1}^{k}w(i)\mathcal{K}_{G}^{(i)}\right\|_{F}^{2}\quad\text{ where }\ w\geq 0,\ \sum_{i=1}^{k}w(i)=1.$

   (16)

The main difference between MKL based approaches is the enforced property for the ideal consensus kernel, using different regularizations. This relates to the criteria used to define the “goodness” of the consensus kernel, and the two most widely used ones promote

• •

  a block-diagonal coefficient matrix, or

• •

  a low-rank consensus kernel.

Let us discuss these two strategies in the next two paragraphs.

Early attempts for employing neural network for nonlinear SC were mainly limited to providing better, that is, more discriminative, features (representations) for the data points. In fact, several SC papers in the literature (including linear and nonlinear methodologies) still rely on manually-designed/handcrafted features, such as Local Binary Patterns (LBP) [77], Histogram of Oriented Gradients (HOG) [17], and Scale-Invariant Feature Transform (SIFT) [61], for satisfactory results on real-world data sets. However, there is no theoretical justification on how any of these handcrafted features might affect the geometric structure of the data lying on multiple (nonlinear) subspaces. From the perspective of unsupervised feature learning, features extracted by neural networks have an important advantage over traditional handcrafted features: they are specifically learned from the data.

The focus of neural network based approaches for SC is to learn proper representations/features from the data points for the specific task of SC. The main shared characteristic of these approaches is that they treat feature learning and clustering as separate tasks, which is suboptimal. They are divided into two main subcategories: (i) transforming the unsupervised problem into a (self)-supervised problem, and (ii) using the connectivity information from the output of a classical linear SC algorithm as prior knowledge. Let us discuss, in the next paragraphs, the approaches within each category in detail.

Separating the two modules of representation/feature learning and (subspace) clustering assignment typically leads to suboptimal results. To overcome this limitation, a group of approaches explicitly include the self-expressiveness term in the network loss function. These approaches are based on the intuition that minimizing self-expressiveness in the embedded space encourages the neural network to learn the appropriate embedding for the task of SC. In general, the optimization problem of these networks has the following form:

where $\Phi_{e}$ and $\Theta$ denote the network embedding and parameters, respectively. The embedding of the network, which is denoted by $Z_{\Theta}$, is encouraged to have a union of linear subspaces structure by explicitly minmizing the self-expressive representation term. As in the rest of the paper, $f(C)$ is the regularization on the coefficient matrix. The extra regularization $h(X,Z_{\Theta},C)$ plays a critical role in removing trivial solutions and specifying the properties of the nonlinear mapping and the embedding space. In particular, looking for self-expressiveness in the transformed space is not sufficient. There exists infinite nonlinear transformations that can lead to small/zero self-expressive error, that is, $||Z_{\Theta}-Z_{\Theta}C||_{F}\approx 0$. Hence, the main difference between these approaches lies on the regularization function $h$ that reduces the possible nonlinear mappings space and forms the geometric structure of the data embedding. There are mainly four choices for $h$ in the literature:

1. 1.

   Instant normalization [82] promotes the norm of the embedded data points in the latent space to be close to 1 using the regularization $\frac{1}{2}\sum_{i=1}^{n}||Z_{\Theta}(:,i)^{\top}Z_{\Theta}(:,i)-1||_{2}^{2}$. This avoids an arbitrary scaling factor in the embedding space.

2. 2.

   Locality preservation [153] captures the local structure of data through Graph Laplacian regularization, similarly as done in Section 3.1.1. This regularization is defined as $\frac{1}{2}\operatorname{tr}(CLC^{\top})$, where $L$ is the Laplacian matrix defined from a similarity matrix over the data points $X$.

3. 3.

   Autoencoder reconstruction loss [42] is based on an autoencoder architecture and uses a decoder network to reconstruct the original data from the self-expressive representation, that is, $Z_{\Theta}C$. More precisely, the regularization is defined as $||X-\tilde{X}_{\Theta}||_{F}^{2}$ where $\tilde{X}_{\Theta}=\Phi_{d}(ZC,\Theta)$, and $\Phi_{d}$ represents the decoding mapping which depends on the parameters $\Theta$. This regularization that minimizes the reconstruction error of the input matrix ensures that the self-expressive based embedding preserves sufficient information to recover the original data.

4. 4.

   Restricted transformations [71] is based on learning stacked linear transformations through multiple layers that are connected via non-negativity constraints inspired by the rectified linear unit (ReLu) that sets negative values to zero after each layer [30, 25]. Specifically, Maggu et al. [71] proposed a three layered transformation for SC, solving the following optimization problem:

   	$\displaystyle\min_{\Theta=\{T_{i}\}_{i=1}^{3},C}\$	$\displaystyle\ \underbrace{||T_{3}T_{2}T_{1}X-Z_{\Theta}||_{F}^{2}}_{\text{deep transformed data: $Z_{\Theta}\approx\Phi_{e}(X,\{T_{i}\}_{i=1}^{3}$)}}+\underbrace{\lambda_{1}||Z_{\Theta}-Z_{\Theta}C||_{F}^{2}+f(C)}_{\text{self-expressiveness in latent space}}$		(21)
   		$\displaystyle+\underbrace{\lambda_{2}\sum_{i=1}^{3}\big{(}||T_{i}||^{2}_{F}-\operatorname{logdet}T_{i}\big{)},}_{\text{avoid trivial transformations}}\quad\text{ such that }\ \underbrace{T_{2}T_{1}X\geq 0\ \text{and}\ T_{1}X\geq 0}_{\text{nonnegative representations}},$

   where $\{T_{i}\}_{i=1}^{3}$ are the linear transformations such that the dimension of each layer output is reduced. To avoid trivial or degenerate solutions ($T_{i}\rightarrow 0\ \text{or}\ T_{i}\rightarrow\infty$), the linear transformation corresponding to each layer is penalized by minimizing the Frobenius norm minus the logarithm of the determinant of the transformation matrices. Note that the non-negativity constraints are inspired by the ReLU activation but they are inherently different.