$\ell_0$-Motivated Low-Rank Sparse Subspace Clustering

TL;DR

This paper introduces novel $ ext{ extonehalf}$-norm based regularizations for low-rank sparse subspace clustering, addressing limitations of traditional nuclear and $ ext{ extonehalf}$-norms, and demonstrates improved performance on synthetic and real datasets.

Contribution

It proposes two $ ext{ extonehalf}$-norm based regularizations, GMC-LRSSC and $S_0/ ext{ extonehalf}$-LRSSC, with algorithms and convergence analysis, advancing subspace clustering methods.

Findings

Outperforms state-of-the-art methods on synthetic data.

Achieves better clustering accuracy on real-world datasets.

Provides convergence guarantees for the proposed algorithms.

Abstract

In many applications, high-dimensional data points can be well represented by low-dimensional subspaces. To identify the subspaces, it is important to capture a global and local structure of the data which is achieved by imposing low-rank and sparseness constraints on the data representation matrix. In low-rank sparse subspace clustering (LRSSC), nuclear and $\ell_1$ norms are used to measure rank and sparsity. However, the use of nuclear and $\ell_1$ norms leads to an overpenalized problem and only approximates the original problem. In this paper, we propose two $\ell_0$ quasi-norm based regularizations. First, the paper presents regularization based on multivariate generalization of minimax-concave penalty (GMC-LRSSC), which contains the global minimizers of $\ell_0$ quasi-norm regularized objective. Afterward, we introduce the Schatten-0 ($S_0$) and $\ell_0$ regularized objective and approximate the proximal map of the joint solution using a proximal average method ($S_0/\ell_0$-LRSSC). The resulting nonconvex optimization problems are solved using alternating direction method of multipliers with established convergence conditions of both algorithms. Results obtained on synthetic and four real-world datasets show the effectiveness of GMC-LRSSC and $S_0/\ell_0$-LRSSC when compared to state-of-the-art methods.

Figures

Tables

Table 1. TABLE I: Notations and abbreviations

Notation	Definition
$N$	Number of data points
$n$	Dimension of data points
$L$	Number of subspaces
$𝑿\in {ℝ}^{n\times N}$	Data matrix
$𝑪\in {ℝ}^{N\times N}$	Representation matrix
$𝑾\in {ℝ}^{N\times N}$	Affinity matrix
$𝑿=𝑼𝚺{𝑽}^T$	Singular value decomposition of $𝑿$
$\sigma \left(𝑿\right)$	Vector of singular values of $𝑿$

Table 2. TABLE II: Clustering error (%) on the Extended Yale B dataset

L	SSC	LRR	LRSC	SSC-OMP	TSC	NSN	LRSSC	${\ell }_0$-SSC	$S_{2/3}$-LRR	$S_{1/2}$-LRR	GMC-LRSSC	$S_0$/${\ell }_0$-LRSSC
5	4.54	4.38	14.56	7.46	28.38	5.42	21.19	16.27	8.94	8.92	3.97	3.52
10	8.78	7.80	34.74	15.26	39.73	9.19	26.89	28.45	9.45	9.49	4.00	4.45
20	21.52	16.68	28.23	17.23	45.57	15.02	36.78	39.24	11.98	11.58	6.38	7.14
30	26.73	21.27	31.53	20.53	47.10	18.69	33.60	39.54	12.02	11.61	8.65	8.35

Table 3. TABLE III: Clustering error (%) on MNIST and USPS datasets

Digits	Dataset	SSC	LRR	LRSC	SSC-	TSC	NSN	LRSSC	${\ell }_0$-SSC	$S_{2/3}$-	$S_{1/2}$-	GMC-	$S_0$/${\ell }_0$-
					OMP					LRR	LRR	LRSSC	LRSSC
2,4,8	MNIST	7.43	14.14	10.59	11.06	12.09	13.02	7.01	7.82	14.80	15.03	8.66	8.92
	USPS	6.02	10.37	7.61	21.04	8.32	18.67	7.13	4.88	8.83	9.70	4.92	6.40
3,6,9	MNIST	3.89	3.49	4.61	5.69	3.25	2.99	4.15	3.49	5.93	6.61	2.93	3.25
	USPS	2.05	1.57	4.50	23.43	1.43	2.22	8.13	1.09	3.02	3.68	0.97	1.27
1,4,7	MNIST	47.50	45.09	44.14	42.21	33.72	14.05	47.09	45.40	43.42	42.57	34.50	27.33
	USPS	2.21	4.01	4.33	58.07	7.45	9.45	8.61	3.68	4.27	5.37	2.65	3.82
2,4,6,8,9	MNIST	25.54	33.54	29.22	29.43	28.95	27.04	26.78	28.94	34.90	36.29	27.40	27.20
	USPS	15.69	22.30	18.86	53.43	20.35	26.69	18.67	14.31	19.22	19.95	15.37	13.38
0,1,3,5,7	MNIST	53.60	33.61	35.92	37.51	30.16	22.39	46.86	33.74	32.50	33.08	29.80	27.85
	USPS	30.00	22.83	28.47	74.66	25.58	13.36	35.44	30.17	27.41	27.68	24.76	10.89
0-9	MNIST	47.49	45.13	46.88	46.45	40.00	34.81	45.68	39.51	43.19	43.66	38.01	34.89
	USPS	28.28	33.44	28.58	84.01	29.34	28.43	33.10	27.27	28.75	29.24	20.46	22.45

Table 4. TABLE IV: Clustering error (%) on the ISOLET1 dataset

L	SSC	LRR	LRSC	SSC-OMP	TSC	NSN	LRSSC	${\ell }_0$-SSC	$S_{2/3}$-LRR	$S_{1/2}$-LRR	GMC-LRSSC	$S_0$/${\ell }_0$-LRSSC
5	10.98	7.61	10.25	27.79	11.58	8.23	8.63	9.23	7.45	8.06	7.07	6.87
10	17.11	14.65	15.84	44.44	19.35	16.04	14.72	18.21	14.10	14.34	13.92	13.81
15	25.64	23.14	22.88	54.03	27.07	23.61	23.87	25.73	20.37	20.24	20.29	19.90
20	31.05	30.60	27.93	59.89	31.70	27.94	29.90	30.28	26.02	25.24	25.32	25.07

Table 5. TABLE V: Average time (s) on the Extended Yale B, MNIST and ISOLET1 datasets

Dataset	L	SSC	LRR	LRSC	SSC-	TSC	NSN	LRSSC	${\ell }_0$-SSC	$S_{2/3}$-	$S_{1/2}$-	GMC-	$S_0$/${\ell }_0$-
					OMP					LRR	LRR	LRSSC	LRSSC
Yale B	5	30.88	10.21	0.18	0.56	0.44	1.28	10.65	40.11	9.85	10.02	0.39	0.59
	10	54.30	31.52	0.41	1.29	1.00	2.94	69.97	101.00	33.78	34.32	1.86	2.91
	30	180.67	271.52	2.47	6.29	6.20	14.84	642.33	623.94	295.74	297.36	19.19	29.24
MNIST	3	0.56	1.17	0.05	0.15	0.17	0.52	0.37	2.63	1.35	1.37	0.12	0.11
	5	1.31	4.01	0.10	0.27	0.29	0.85	1.40	5.11	5.03	5.10	0.30	0.26
	10	4.59	13.56	0.26	0.60	0.67	1.83	12.44	14.37	17.09	17.41	1.45	1.18
ISOLET1	5	6.83	4.78	0.09	0.20	0.25	0.60	4.24	4.71	5.99	6.12	0.41	2.36
	10	14.23	17.23	0.26	0.47	0.61	1.31	30.05	12.82	21.78	22.21	1.87	12.95
	20	38.37	30.34	0.28	1.17	1.64	3.05	133.86	45.79	34.05	34.87	6.96	50.93

Equations

$$\Theta (x;\lambda ,a)=\{\begin{array}{ll}0,& \text{if}|x|\le \lambda \\ a(|x|-\lambda )/(a-\lambda )sign(x),& \text{if}\lambda \le |x|\le a\\ x,& \text{if}|x|\ge a.\end{array}$$

$$H(x;\lambda )=\{\begin{array}{ll}x,& \text{if}|x|>\sqrt{2\lambda }\\ \{0,x\},& \text{if}|x|=\sqrt{2\lambda }\\ 0,& \text{if}|x|<\sqrt{2\lambda }.\end{array}$$

$$\|\bm{A}\|_{0}=\#\big{\{}a_{ij}\neq 0,\ i=1..N,\ j=1..M\big{\}},$$

$$\Vert \mathbf{A}{\Vert }_{S_0}=\Vert diag(\mathbf{\Sigma }){\Vert }_0,$$

$$\min_{\bm{C}}\big{\|}\bm{C}\big{\|}_{*}\ \ s.t.\ \ \bm{X}=\bm{XC},$$

$$\hat{\mathbf{C}}=\mathbf{V}{\mathbf{V}}^T.$$

$$\min_{\bm{C}}\frac{\lambda}{2}\big{\|}\bm{X}-\bm{XC}\big{\|}_{F}^{2}+\big{\|}\bm{C}\big{\|}_{*},$$

$$\hat{\mathbf{C}}={\mathbf{V}}_1(\mathbf{I}-\frac{1}{\lambda }{\mathbf{\Sigma }}_1^{-2}){{\mathbf{V}}_1}^T,$$

$$\min_{\bm{C}}\big{\|}\bm{C}\big{\|}_{1}\ \ s.t.\ \ \bm{X}=\bm{XC},\ diag(\bm{C})=\bm{0},$$

$$\min_{\bm{C}}\frac{1}{2}\big{\|}\bm{X}-\bm{XC}\big{\|}_{F}^{2}+\tau\big{\|}\bm{C}\big{\|}_{1}\ \ s.t.\ \ diag(\bm{C})=\bm{0},$$

$$\min_{\bm{C}}\lambda\big{\|}\bm{C}\big{\|}_{*}+\tau\big{\|}\bm{C}\big{\|}_{1}\ \ s.t.\ \ \bm{X}=\bm{XC},\ diag(\bm{C})=\bm{0},$$

$$\begin{split}\min_{\bm{C}}\frac{1}{2}\big{\|}\bm{X}&-\bm{XC}\big{\|}_{F}^{2}+\lambda\big{\|}\bm{C}\big{\|}_{*}+\tau\big{\|}\bm{C}\big{\|}_{1}\\ &\ \ s.t.\ \ diag(\bm{C})=\bm{0}.\end{split}$$

$$\mathbf{W}=|\mathbf{C}|+|\mathbf{C}{|}^T.$$

$$\begin{split}\min_{\bm{C}}\frac{1}{2}\big{\|}\bm{X}&-\bm{XC}\big{\|}_{F}^{2}+\lambda g(\bm{C})+\tau f(\bm{C})\\ &\ \ s.t.\ \ diag(\bm{C})=\bm{0},\end{split}$$

$${\psi }_{\mathbf{B}}(\mathbf{z})=\Vert \mathbf{z}{\Vert }_1-S_{\mathbf{B}}(\mathbf{z}),$$

$$S_{\bm{B}}(\bm{z})=\inf_{\textbf{v}\in\mathbb{R}^{N}}\Big{\{}\|\bm{v}\|_{1}+\frac{1}{2}\|\bm{B}(\bm{z}-\bm{v})\|_{2}^{2}\Big{\}}.$$

$$F(\mathbf{z})=\frac{1}{2}\Vert \mathbf{y}-\mathbf{Az}{\Vert }_2^2+\lambda {\psi }_{\mathbf{B}}(\mathbf{z}),$$

$$\mathbf{B}=\sqrt{\gamma /\lambda }\mathbf{A},0\le \gamma \le 1.$$

$${\mathbf{A}}^T\mathbf{A}=diag({\alpha }_1^2,...,{\alpha }_N^2),$$

$${\mathbf{z}}_n^{\text{opt}}=\Theta ([{\mathbf{A}}^T\mathbf{y}{]}_n/{\alpha }_n^2;\lambda /{\alpha }_n^2,\lambda /(\gamma {\alpha }_n^2)),$$

$$f(z_1,z_2,...,z_N)=f(|z_{\pi (1)}|,|z_{\pi (2)}|,...,|z_{\pi (N)}|),$$

$$\partial [({\psi }_{\mathbf{B}}\circ \sigma )(\mathbf{X})]=\mathbf{U}diag(\partial {\psi }_{\mathbf{B}}[{\sigma }_i(\mathbf{X})]){\mathbf{V}}^T,$$

$${\mathbf{B}}^T\mathbf{B}=diag({{\alpha }_1}^2,...,{{\alpha }_N}^2)\Rightarrow {\psi }_{\mathbf{B}}(\mathbf{z})=\sum _{n=1}^N{\varphi }_{{\alpha }_n}(z_n),$$

$${\varphi }_b(y)=\{\begin{array}{ll}|y|-\frac{1}{2}{b}^2{y}^2,& \text{if}|y|\le 1/b^2,\\ \frac{1}{2b^2},& \text{otherwise}.\end{array}$$

$$\begin{split}\min_{\bm{C}}\frac{1}{2}\big{\|}\bm{X}-&\bm{XC}\big{\|}^{2}_{F}+\lambda\psi_{\bm{B}}({\sigma}(\bm{C}))+\tau\psi_{\bm{B}}(\bm{C})\\ &\ \ s.t.\ \ diag(\bm{C})=\bm{0},\end{split}$$

$$\begin{split}\min_{\bm{J},\bm{C}_{1},\bm{C}_{2}}&\frac{1}{2}\big{\|}\bm{X}-\bm{XJ}\big{\|}^{2}_{F}+\lambda\psi_{\bm{B}}(\sigma(\bm{C}_{1}))+\tau\psi_{\bm{B}}(\bm{C}_{2})\ \ \\ &s.t.\ \ \bm{J}=\bm{C}_{1},\ \ \bm{J}=\bm{C}_{2}-diag(\bm{C}_{2}),\end{split}$$

$$\begin{split}&\mathcal{L}_{\mu_{1},\mu_{2}}\big{(}\bm{J},\bm{C}_{1},\bm{C}_{2},\bm{\Lambda}_{1},\bm{\Lambda}_{2}\big{)}=\frac{1}{2}\big{\|}\bm{X}-\bm{XJ}\big{\|}^{2}_{F}+\lambda\psi_{\bm{B}}(\sigma(\bm{C}_{1}))\\ &+\tau\psi_{\bm{B}}(\bm{C}_{2})+\frac{\mu_{1}}{2}\big{\|}\bm{J}-\bm{C}_{1}\big{\|}^{2}_{F}+\frac{\mu_{2}}{2}\big{\|}\bm{J}-\bm{C}_{2}+diag(\bm{C}_{2})\big{\|}^{2}_{F}\\ &+\big{\langle}\bm{\Lambda}_{1},\bm{J}-\bm{C}_{1}\big{\rangle}+\big{\langle}\bm{\Lambda}_{2},\bm{J}-\bm{C}_{2}+diag(\bm{C}_{2})\big{\rangle},\end{split}$$

$$\begin{split}&\min_{\bm{J}}\mathcal{L}_{\mu_{1}^{k},\mu_{2}^{k}}\big{(}\bm{C}_{1}^{k},\bm{C}_{2}^{k},\bm{J},\bm{\Lambda}_{1}^{k},\bm{\Lambda}_{2}^{k}\big{)}=\min_{\bm{J}}\frac{1}{2}\big{\|}\bm{X}-\bm{XJ}\big{\|}_{F}^{2}\\ &+\frac{\mu_{1}}{2}\big{\|}\bm{J}-\bm{C}_{1}^{k}\big{\|}_{F}^{2}+\frac{\mu_{2}}{2}\big{\|}\bm{J}-\bm{C}_{2}^{k}+diag(\bm{C}_{2}^{k})\big{\|}^{2}_{F}+\big{\langle}\bm{\Lambda}_{1}^{k},\bm{J}-\bm{C}_{1}^{k}\big{\rangle}\quad\\ &+\big{\langle}\bm{\Lambda}_{2}^{k},\bm{J}-\bm{C}_{2}^{k}+diag(\bm{C}_{2}^{k})\big{\rangle}.\end{split}$$

$$\bm{J}^{k+1}=\big{[}\bm{X}^{T}\bm{X}+(\mu_{1}^{k}+\mu_{2}^{k})\bm{I}\big{]}^{-1}\big{[}\bm{X}^{T}\bm{X}+\mu_{1}^{k}{\bm{C}_{1}}^{k}+\mu_{2}^{k}{\bm{C}_{2}}^{k}-{\bm{\Lambda}_{1}}^{k}-{\bm{\Lambda}_{2}}^{k}\big{]}.$$

$$\begin{split}&\min_{\bm{C}_{1}}\mathcal{L}_{\mu_{1}^{k},\mu_{2}^{k}}\big{(}\bm{J}^{k+1},\bm{C}_{1},\bm{C}_{2}^{k},\bm{\Lambda}_{1}^{k},\bm{\Lambda}_{2}^{k}\big{)}\\ &=\min_{\bm{C}_{1}}\lambda\psi_{\bm{B}}(\sigma(\bm{C}_{1}))+\frac{\mu_{1}^{k}}{2}\big{\|}\bm{J}^{k+1}-\bm{C}_{1}\big{\|}_{F}^{2}+\big{\langle}\bm{\Lambda}_{1}^{k},\bm{J}^{k+1}-\bm{C}_{1}\big{\rangle}\\ &=\min_{\bm{C}_{1}}\lambda\psi_{\bm{B}}(\sigma(\bm{C}_{1}))+\frac{\mu_{1}^{k}}{2}\Big{\|}\bm{J}^{k+1}+\frac{\bm{\Lambda}_{1}^{k}}{\mu_{1}^{k}}-\bm{C}_{1}\Big{\|}_{F}^{2}.\end{split}$$

Full text

††publicationid: pubid: 978-1-5386-5541-2/18/$31.00 ©2018 IEEE.

I Introduction

High dimensional data analysis is a widespread problem in many applications of machine learning, computer vision, and bioinformatics [1, 2, 3, 4, 5, 6]. However, in many real-world datasets, the intrinsic dimension of high-dimensional data is much smaller than the dimension of the ambient space and data can be well represented as lying close to a union of low-dimensional subspaces. The problem of segmenting data according to the low-dimensional subspaces they are drawn from is known as subspace clustering [7]. Thanks to their capability to handle arbitrarily shaped clusters and their well-defined mathematical principles, spectral based methods [8, 9] are widely used approaches to subspace clustering. These methods solve the subspace clustering problem by relying on the spectral graph theory and cluster eigenvectors of the graph Laplacian matrix corresponding to the smallest eigenvalues [10].

One of the main challenges in subspace clustering is the construction of the affinity matrix that captures well (di)similarities between data points. Among various approaches proposed in the literature, methods based on sparse [11, 12] and low-rank representations [13, 14, 15] have been among the most successful in many applications [16]. These methods exploit the self-expressiveness property of the data and represent each data point as a linear combination of other data points in the dataset. Low-rank representation (LRR) [13, 14, 17] captures the global structure of the data by imposing a low-rank constraint on the data representation matrix. Low-rank implies that representation matrix is described by a weighted sum of small number of outer products of left and right singular vectors. In order to ensure convexity of the related optimization problem, the rank minimization is relaxed as the nuclear or Schatten-1 norm minimization problem [18, 19, 20]. Different from LRR, Sparse Subspace Clustering (SSC) [21, 11] captures local linear relationships by constraining representation matrix to be sparse. Using the tightest convex relaxation of the $\ell_{0}$ quasi-norm, the SSC model solves sparsity maximization problem as $\ell_{1}$ norm minimization problem [22, 23]. Both LRR and SSC guarantee exact clustering when subspaces are independent, but the independence assumption is overly restrictive for many real-world datasets [24, 25]. Under appropriate conditions [26], SSC also succeeds for disjoint subspaces. However, when the number of dimensions is higher than three, SSC can face connectivity problems resulting in a disconnected graph within a subspace [27]. A natural way to construct adaptive model able to capture the global and the local structure of the data is to constrain representation matrix to be low-rank and sparse. In [28, 29, 16, 30] that is done by combining nuclear and $\ell_{1}$ norms as the measures of rank and sparsity, respectively. The motivation lies in the fact that minimization of these norms results in a convex optimization problem.

Although convex, nuclear and $\ell_{1}$ norms are not exact measures of rank and sparsity. Therefore, optimal solution of the nuclear and $\ell_{1}$ norms regularized objective is only approximate solution of the original problem [31]. Proximity operator associated with the nuclear norm overpenalizes large singular values, leading to biased results in low-rank constrained optimization problems [32, 33]. Similarly, in sparsity regularized problems $\ell_{1}$ norm solution systematically underestimates high amplitude components of sparse vectors [34]. Nonconvex regularizations based on $\ell_{p}$ quasi-norms ($0\leq p<1$) or their approximations have been proposed for various low-rank [32, 33, 35, 36, 37, 38, 39] and sparsity regularized problems [40, 34, 41, 42, 43, 44, 45, 46, 47]. Recently, nonconvex approximations of rank and sparsity have also been introduced in subspace clustering problem [48, 49, 50, 51, 52]. Specifically, $\ell_{0}$-induced sparse subspace clustering is introduced in [48]. The corresponding optimization problem is solved using proximal gradient descent which under assumptions on the sparse eigenvalues converges to a critical point. In [49] authors replaced the nuclear norm regularizer with the nonconvex Ky Fan $p$-$k$ norm [53] and proposed proximal iteratively reweighted optimization algorithm to solve the problem. In [50, 54] rank is approximated using Schatten-$q$ quasi-norm regularization ($0<q<1$). The optimization problem in [50] is solved using generalized matrix soft thresholding algorithm [55]. Schatten-$q$ quasi-norm minimization with tractable $q=2/3$ and $q=1/2$ is proposed in [51]. The Schatten-$q$ ($S_{q}$) quasi-norm for $0<q<1$ is equivalent to $\ell_{q}$ quasi-norm on vector of singular values. Compared to the nuclear norm, it makes a closer approximation of the rank function. In this regard, $S_{0}$ quasi-norm can be seen as an equivalent to the $\ell_{0}$ quasi-norm and stands for the definition of the rank function. Furthermore, [52] combines $S_{q}$ regularizer ($0<q<1$) for low-rank and $\ell_{p}$ quasi-norm regularizer ($0<p<1$) for sparsity constraint. However, recent results in [56] show that in $\ell_{p}$ regularized least squares ($0\leq p<1$) smaller values of $p$ lead to more accurate solutions. If $\ell_{1}$ norm is also considered, authors show that for large measurement noises $\ell_{1}$ outperforms $\ell_{p}$, $p<1$. However, for small measurement noises $\ell_{0}$ quasi-norm regularization outperforms $\ell_{p}$, $0<p\leq 1$.

Motivated by the limitations discussed above, we introduce two $S_{0}$/$\ell_{0}$ quasi-norm based nonconvex regularizations for low-rank sparse subspace clustering (LRSSC). First, we propose regularization based on multivariate generalization of the minimax-concave penalty function (GMC), introduced in [34] for sparsity regularized linear least squares. Here, this approach is extended to the rank approximation. The GMC penalty enables to maintain the convexity of low-rank and sparsity constrained subproblems, while achieving better approximation of rank and sparsity than nuclear and $\ell_{1}$ norms. Importantly, this regularization is closely related to the continuous exact $\ell_{0}$ penalty which contains the global minimizers of $\ell_{0}$ quasi-norm regularized least-squares objective[57, 34]. GMC penalty yields solutions of low-rank and sparsity constrained subproblems based on firm thresholding of singular values and coefficients of representation matrix, respectively. The firm thresholding function $\Theta:\mathbb{R}\rightarrow\mathbb{R}$ is is defined as [58]:

[TABLE]

Next, we propose the direct solution of $S_{0}$ and $\ell_{0}$ quasi-norms regularized objective function. The solution of corresponding low-rank and sparsity constrained subproblems is based on iterative application of hard thresholding operator [59, 60, 61] on the singular values and coefficients of the representation matrix, respectively. The hard thresholding function $H:\mathbb{R}\rightarrow\mathbb{R}$ is defined as [59]:

[TABLE]

Simultaneous rank and sparsity regularization is handled using the proximal average method, introduced in [62] for convex problems and extended recently to nonconvex and nonsmooth functions [63, 64]. Proximal average allows us to approximate the proximal map of joint solution by averaging solutions obtained separately from low-rank and sparsity subproblems, leading to a problem with a low computational cost in each iteration. Furthermore, using proximal average method enables us to establish global convergence guarantee for $S_{0}/\ell_{0}$ regularized LRSSC.

Better approximation of rank and sparsity is a consequence of the properties of firm and hard thresholding operators associated with GMC and $\ell_{0}$ regularizations. As opposed to them, the soft thresholding operator underestimates high amplitude coefficients in $\ell_{1}$ norm based sparsity regularized objective, as well as large singular values in low-rank approximation problem. As an example, Fig. 1 shows soft, firm and hard thresholding operators used in LRSSC, GMC-LRSSC and $S_{0}/\ell_{0}$-LRSSC, respectively.

To solve corresponding optimization problems we derive algorithms based on computationally efficient Alternating Direction Method of Multipliers (ADMM) [65]. Although ADMM has been successfully applied for many nonconvex problems [66, 67, 68], only recent theoretical results establish convergence of ADMM for certain nonconvex functions [69, 70, 71, 72]. For GMC regularization, we show that the sequence generated by the algorithm is bounded and prove that any limit point of the iteration sequence is a stationary point. For $S_{0}/\ell_{0}$ regularization with proximal average approach, based on the property that $\ell_{0}$ and $S_{0}$ quasi-norms belong to a class of semialgebraic functions and satisfy the Kurdyka-Łojasiewicz inequality [73], the global convergence of the algorithm can be guaranteed.

Experimental results on synthetic and four real-world datasets demonstrate that the proposed $\ell_{0}$ based low-rank sparse subspace clustering algorithms converge fast and to a point with lower or similar clustering error than the convex approximations with nuclear and $\ell_{1}$ norms. Compared to the state-of-the-art subspace clustering methods, the proposed algorithms perform better on four benchmark datasets.

I-A Contributions

The contributions of this paper are summarized as follows:

We introduce nonconvex generalized minimax-concave penalty in the low-rank sparse subspace clustering problem, such that the global minimizers of the proposed objective coincide with that of a convex function defined using the continuous exact $\ell_{0}$ penalty [57]. The introduced penalty maintains the convexity of the sparsity and low-rank constrained subproblems. The proximal operator of the related GMC penalty function is the firm thresholding function [34]. 2. 2.

We introduce $S_{0}$ and $\ell_{0}$ pseudo-norm regularizations for LRSSC. Using the proximal average method [63, 62], we average the solutions of proximal maps of low-rank and sparsity subproblems, with the hard thresholding function as a proximity operator of the related penalties [59]. 3. 3.

We derive ADMM based optimization algorithms for LRSSC constrained either with a GMC penalty or with $S_{0}$/$\ell_{0}$ quasi-norms. Iterative firm or hard thresholding of singular values and coefficients of representation matrix is used to obtain the solution of rank and sparsity constrained subproblems. 4. 4.

We prove that the sequence generated by the GMC regularized LRSSC algorithm is bounded and that any limit point of the iteration sequence is a stationary point that satisfies Karush-Kuhn-Tucker (KKT) conditions. 5. 5.

We establish the convergence property of the $S_{0}$/$\ell_{0}$ regularized approach with proximal average and show that the algorithm converges regardless of the initialization. To the best of our knowledge, we are the first to show convergence with $S_{0}$ and $\ell_{0}$ penalties in the low-rank and sparsity constrained optimization problem.

The remainder of this paper is organized as follows. Section II gives a brief overview of the related work. In Section III and IV we introduce GMC and $S_{0}/\ell_{0}$ regularized low-rank sparse subspace clustering methods, respectively. We formulate the problem, present optimization algorithms, and analyze convergence and computational complexity. The experimental results on synthetic and four real-world datasets are presented in Section V. Finally, Section VI concludes this paper.

I-B Main Notation

Scalars are denoted by lower case letters, vectors by bold lower-case letters, matrices are denoted by bold capital and subspaces by calligraphic letters. $\|\cdot\|_{F}$ denotes Frobenius norm defined as the square root of the sum of the squares of matrix elements. $\|\cdot\|_{1}$ denotes $\ell_{1}$ norm defined as the sum of absolute values of matrix elements. $\|\cdot\|_{*}$ denotes nuclear norm defined as the sum of singular values of a matrix. $\ell_{0}$ quasi-norm is denoted by $\|\cdot\|_{0}$ and for matrix $\bm{A}\in\mathbb{R}^{N\times M}$ defined as:

[TABLE]

where $\#$ denotes cardinality function. Schatten-[math] quasi norm is denoted by $\|\cdot\|_{S_{0}}$ and defined as:

[TABLE]

where $\bm{A}=\bm{U\Sigma V}^{T}$ is the singular value decomposition (SVD) of matrix $\bm{A}$. Since $\ell_{0}$ quasi-norm does not satisfy homogeneous property it is not a norm, but with a slight abuse of notation we will refer to it as the $\ell_{0}$ norm in the rest of the paper. Null vector is denoted by $\bm{0}$ and $diag(\cdot)$ is the vector of diagonal elements of a matrix. Table I summarizes some notations used in the paper.

II Background

Consider the data matrix $\bm{X}\in\mathbb{R}^{n\times N}$ the columns of which are data points drawn from a union of $L$ linear subspaces $\bigcup_{i=1}^{L}\mathcal{S}_{i}$ of unknown dimensions \big{\{}d_{i}=\textnormal{dim}(\mathcal{S}_{i})\big{\}}_{i=1}^{L} in $\mathbb{R}^{n}$. Let $\bm{X}_{i}\in\mathbb{R}^{n\times N_{i}}$ be a submatrix of $\bm{X}$ of rank $d_{i}$, $0<d_{i}<N_{i}$ and $\sum_{i=1}^{L}N_{i}=N$. Given data matrix $\bm{X}$, subspace clustering segments data points according to the low-dimensional subspaces. The first step is the construction of the affinity matrix $\bm{W}\in\mathbb{R}^{N\times N}$ whose elements represent the similarity between data points. An ideal affinity matrix is block diagonal (up to a permutation): non-zero distance is assigned to the points in the same subspace and zero distance to the points from different subspaces. Spectral clustering algorithm [8, 9] is then applied to the affinity matrix to obtain memberships of data points to the subspaces.

II-A Related Work

Low-Rank Representation (LRR) [13, 14] aims to find a low-rank representation matrix $\bm{C}\in\mathbb{R}^{N\times N}$ for input data matrix $\bm{X}$ by solving the following convex optimization problem:

[TABLE]

where the nuclear norm is used to approximate the rank of $\bm{C}$. Let $\bm{X}=\bm{U\Sigma V}^{T}$ be the SVD of $\bm{X}$. The closed form solution of problem (3) is given by [14]:

[TABLE]

When data points are contaminated by additive white Gaussian noise (AWGN), the following minimization problem is solved:

[TABLE]

where $\lambda$ is the rank regularization constant. The optimal solution of problem (5) is given by [17, 74]:

[TABLE]

where $\bm{U}=[\bm{U}_{1}\ \bm{U}_{2}]$, $\bm{\bm{\Sigma}}=diag(\bm{\bm{\Sigma}}_{1}\ \bm{\bm{\Sigma}}_{2})$ and $\bm{V}=[\bm{V}_{1}\ \bm{V}_{2}]$. Matrices are partitioned according to the sets $\mathcal{I}_{1}=\{i:\sigma_{i}>\frac{1}{\sqrt{\lambda}}\}$ and $\mathcal{I}_{2}=\{i:\sigma_{i}\leq\frac{1}{\sqrt{\lambda}}\}$, where $\sigma_{i}$ denotes $i$th singular value of $\bm{X}$.

Sparse Subspace Clustering (SSC) [11] represents each data point as a sparse linear combination of other data points and solves the following convex optimization problem:

[TABLE]

where constraint $diag(\bm{C})=\bm{0}$ is used to avoid trivial solution of representing a data point as a linear combination of itself.

For data contaminated by the AWGN, the following minimization problem is solved to approximate sparse representation matrix $\bm{C}$:

[TABLE]

where $\tau$ is the sparsity regularization constant. This problem can be solved efficiently using ADMM optimization procedure [65, 11].

Low-Rank Sparse Subspace Clustering (LRSSC) [29] requires that the representation matrix $\bm{C}$ is simultaneously low-rank and sparse. LRSSC solves the following problem:

[TABLE]

where $\lambda$ and $\tau$ are rank and sparsity regularization constants, respectively. For the AWGN corrupted data the following problem needs to be solved to approximate $\bm{C}$:

[TABLE]

After representation matrix $\bm{C}$ is estimated, the affinity matrix $\bm{W}\in\mathbb{R}^{N\times N}$ is calculated as follows:

[TABLE]

In the next two sections, we introduce two nonconvex regularizers for the low-rank sparse subspace clustering. We formulate low-rank sparse subspace clustering problem in the following general form:

[TABLE]

where $g(\bm{C})$ and $f(\bm{C})$ are functions that, respectively, measure rank and sparsity of the data representation matrix $\bm{C}$. The convex formulation used in (10) implies g(\bm{C})=\big{\|}\bm{C}\big{\|}_{*} and f(\bm{C})=\big{\|}\bm{C}\big{\|}_{1}.

III GMC-LRSSC Algorithm

III-A Problem Formulation

We propose to regularize rank and sparsity using multivariate GMC penalty function, introduced in [34] for sparse regularized least-squares. We start with some definitions and results that will be used throughout the paper.

Definition 1** ([34])**

Let $\bm{z}\in\mathbb{R}^{N}$ and $\bm{B}\in\mathbb{R}^{M\times N}$. The GMC penalty function $\psi_{\bm{B}}:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is defined as:

[TABLE]

where $S_{\bm{B}}:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is the generalized Huber function defined as:

[TABLE]

Lemma 1** ([34])**

Let $\bm{z}\in\mathbb{R}^{N}$, $\bm{y}\in\mathbb{R}^{M}$, $\bm{A}\in\mathbb{R}^{M\times N}$ and $\lambda>0$. Define $F:\mathbb{R}^{N}\rightarrow\mathbb{R}$ as:

[TABLE]

where $\psi_{\bm{B}}:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is the GMC penalty. If $\bm{A}^{T}\bm{A}-\lambda\bm{B}^{T}\bm{B}$ is positive semidefinite matrix, $F$ is a convex function. The convexity condition is satisfied by setting:

[TABLE]

The parameter $\gamma$ controls the nonconvexity of the penalty $\psi_{\bm{B}}$. Larger values of $\gamma$ increase the nonconvexity of the penalty. $\ell_{1}$ norm can be seen as a special case of this penalty by setting $\gamma=0$.

Lemma 2** ([34])**

Let $\bm{z}\in\mathbb{R}^{N}$, $\bm{y}\in\mathbb{R}^{M}$, $\bm{A}\in\mathbb{R}^{M\times N}$ and $\lambda>0$. If $\bm{A}^{T}\bm{A}$ is diagonal with positive entries and $\bm{B}$ is given by (16), then for $0<\gamma\leq 1$ the minimizer of $F$ is given by element-wise firm thresholding. Formally, if

[TABLE]

then

[TABLE]

where $\Theta$ stands for the firm thresholding function [58] defined entry-wise in (1).

Definition 2** ([75, 76])**

Function $f:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is an absolutely symmetric function, if:

[TABLE]

holds for any permutation $\pi$ of $\{1,...,N\}$.

Proposition 1

Let $\bm{B}^{T}\bm{B}$ be a diagonal matrix and $\psi_{\bm{B}}$ be the GMC penalty function defined in (13). The subdifferential of singular value function $\psi_{\bm{B}}\circ\sigma$ of a matrix $\bm{X}$ is given by the following equation:

[TABLE]

where $\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{T}$ is the SVD of $\bm{X}$.

Proof:

It follows from [34] that if $\bm{B}^{T}\bm{B}$ is a diagonal matrix, the GMC penalty $\psi_{\bm{B}}$ is separable, comprising a sum of scalar MC penalties:

[TABLE]

where $\phi_{b}:\mathbb{R}\rightarrow\mathbb{R}$ is the scaled MC penalty [77, 34] defined as:

[TABLE]

Therefore, according to Definition 2, $\psi_{\bm{B}}$ is an absolutely symmetric function. The proof of the proposition then follows from the property of the singular value function $f\circ\sigma$ [76], where $f$ is an absolutely symmetric function. ∎

Proposition 1 allows us to use GMC penalty for rank approximation. We formulate GMC penalty regularized objective for low-rank sparse subspace clustering. Let $\bm{B}\in\mathbb{R}^{N\times N}$, and let ${\sigma}(\bm{C})$ denote vector of singular values of $\bm{C}$. By choosing $g(\bm{C})=\psi_{\bm{B}}({\sigma}(\bm{C}))$ as a rank function, and $f(\bm{C})=\psi_{\bm{B}}(\bm{C})$ as a sparsity function in equation (12), we define the following nonconvex objective function:

[TABLE]

where $\psi_{\bm{B}}$ denotes GMC penalty defined in (13), regularized by matrix $\bm{B}$. In the next section we will show that by solving the objective (23) with ADMM, both sparsity and low-rank subproblems can be reduced to the equation (15) with diagonal $\bm{A}^{T}\bm{A}$. In this case, GMC penalty is closely related to the continuous exact $\ell_{0}$ penalty [34, 57], that approximates the convex hull of the least squares with $\ell_{0}$ regularization. Furthermore, diagonal $\bm{A}^{T}\bm{A}$ reduces both subproblems to element-wise firm thresholding function, defined in (1). In low-rank minimization subproblem, the firm thresholding operator needs to be applied to the vector of singular values.

III-B Optimization Algorithm

To solve optimization problem in (23), we introduce auxiliary variables $\bm{J}$, $\bm{C}_{1}$ and $\bm{C}_{2}$ to split variables and solve subproblems independently. The reformulated objective for GMC penalty in (23) is equivalent to:

[TABLE]

The augmented Lagrangian function of (24) is:

[TABLE]

where $\mu_{1}$, $\mu_{2}>0$ are penalty parameters and $\bm{\Lambda}_{1}$, $\bm{\Lambda}_{2}$ are Lagrange multipliers.

Update rule for $\bm{J}^{k+1}$: Given $\bm{C}_{1}^{k}$, $\bm{C}_{2}^{k}$, $\bm{\Lambda}_{1}^{k}$, $\bm{\Lambda}_{2}^{k}$, $\mu_{1}^{k}$, $\mu_{2}^{k}$, we minimize the Lagrangian function in (25) with respect to $\bm{J}$:

[TABLE]

The optimal solution of (26) is given by the following update:

[TABLE]

Update rule for $\bm{C}_{1}^{k+1}$: Given $\bm{J}^{k+1}$, $\bm{\Lambda}_{1}^{k}$, $\mu_{1}^{k}$, we minimize the Lagrangian function in (25) with respect to $\bm{C}_{1}$:

[TABLE]

It can be seen that (28) corresponds to the least squares problem in (15) with $\bm{A}^{T}\bm{A}=\bm{I}$ and therefore, diagonal. It follows from the condition (16) that in order to maintain convexity of the subproblem, we need to set $\bm{B}=\sqrt{\mu_{1}^{k}\gamma/\lambda}\bm{I}$, $0<\gamma\leq 1$. Using Lemma 2 and Proposition 1, (28) can be solved by element-wise firm thresholding of singular values of matrix \big{(}\bm{J}^{k+1}+\bm{\Lambda}_{1}^{k}/\mu_{1}^{k}\big{)}.

Specifically, let $\bm{U\Sigma V}^{T}$ denote the SVD of matrix \big{(}\bm{J}^{k+1}+{\bm{\Lambda}_{1}}^{k}/\mu_{1}^{k}\big{)}. The closed-form solution of (28) is given by:

[TABLE]

where $\Theta$ is the firm thresholding function defined in (1).

Update rule for $\bm{C}_{2}^{k+1}$: Given $\bm{J}^{k+1}$, $\bm{\Lambda}_{2}^{k}$, $\mu_{2}^{k}$, we minimize the objective (25) with respect to $\bm{C}_{2}$:

[TABLE]

with subtraction of diagonal elements of $\bm{C}_{2}^{k+1}$:

[TABLE]

Similarly as the update for $\bm{C}_{1}$, matrix $\bm{A}^{T}\bm{A}$ is diagonal matrix and we can ensure convexity of the subproblem (30) by setting $\bm{B}=\sqrt{\mu_{2}^{k}\gamma/\tau}\bm{I}$, $0<\gamma\leq 1$. The problem (30) is then solved by firm thresholding elements of matrix \big{(}\bm{J}^{k+1}+{\bm{\Lambda}_{2}^{k}}/\mu_{2}^{k}\big{)} and given by:

[TABLE]

Update rules for Lagrange multipliers $\bm{\Lambda}_{1}^{t+1},\bm{\Lambda}_{2}^{t+1}$: Given $\bm{J}^{k+1}$, $\bm{C}_{1}^{k+1}$, $\bm{C}_{2}^{k+1}$, $\mu_{1}^{k}$, $\mu_{2}^{k}$, Lagrange multipliers are updated with the following equations:

[TABLE]

Penalty parameters $\mu_{1}$, $\mu_{2}$ are are in each step $k$ updated according to:

[TABLE]

where $\rho>1$ is step size for adaptively changing $\mu_{1},\mu_{2}$. Due to numerical reasons $\mu_{1}$, $\mu_{2}$ are bounded with $\mu^{max}$, while in the convergence proof we use formulation ${\mu_{i}}^{k+1}=\rho{\mu_{i}}^{k}$, $i=1,2$.

The main steps of the proposed algorithm are summarized in Algorithm 1.

III-C Convergence Analysis

Although choosing $\gamma\in[0,1]$ guarantees convexity of the low-rank and sparsity subproblems and convergence of related subproblems, the objective in (23) is nonconvex. In this section, we analyze convergence of the proposed method and show that any limit point of iteration sequence satisfies Karush-Kuhn-Tucker (KKT) conditions [78].

Proposition 2

The sequences \big{\{}(\bm{J}^{k},\bm{C}_{1}^{k},\bm{C}_{2}^{k},\bm{\Lambda}_{1}^{k},\bm{\Lambda}_{2}^{k})\big{\}} generated by Algorithm 1 are all bounded.

We now state the main theorem related to convergence property of GMC-LRSSC algorithm.

Theorem 1

Let Y^{k}=\big{\{}\big{(}\bm{J}^{k},\bm{C}_{1}^{k},\bm{C}_{2}^{k},\bm{\Lambda}_{1}^{k},\bm{\Lambda}_{2}^{k}\big{)}\big{\}}_{k=1}^{\infty} be a sequence generated by Algorithm 1. Suppose that \lim_{k\rightarrow\infty}\big{(}Y^{k+1}-Y^{k}\big{)}=\bm{0}. Then, any accumulation point of the sequence \{Y^{k}\big{\}}_{k=1}^{\infty} satisfies the Karush-Kuhn-Tucker (KKT) conditions for problem (24). In particular, whenever \{Y^{k}\big{\}}_{k=1}^{\infty} converges, it converges to a point that satisfies KKT conditions.

The proofs of Proposition 2 and Theorem 1 are given in the Appendix.

III-D Stopping Criteria and Computational Complexity

The steps in Algorithm 1 are repeated until convergence or until the maximum number of iterations is exceeded. We check the convergence by verifying the following inequalities at each iteration $k$: \big{\|}\bm{J}^{k}-\bm{C}_{1}^{k}\|_{\infty}\leq\epsilon, \big{\|}\bm{J}^{k}-\bm{C}_{2}^{k}\|_{\infty}\leq\epsilon, \big{\|}\bm{J}^{k}-\bm{J}^{k-1}\|_{\infty}\leq\epsilon. We found that setting error tolerance to $\epsilon=10^{-4}$ works well in practice. In each step we normalize columns of matrix $\mathbf{J}$. This normalization is frequently applied to stabilize convergence of non-negative matrix factorization algorithms [79].

The computational complexity of Algorithm 1 is $O(nN^{2}+TN^{3})$, where $T$ denotes the number of iterations. In the experiments, we set the maximal $T$ to $100$, but on all datasets the algorithm converged within less than 15 iterations. Note that the computational complexity of spectral clustering step is $O(N^{3})$.

IV $S_{0}/\ell_{0}$-LRSSC Algorithm

IV-A Problem Formulation

In addition to the GMC penalty, we propose to directly use $S_{0}$ and $\ell_{0}$ as constrains for low-rank and sparsity. Specifically, by choosing $g(\bm{C})=\|\bm{C}\|_{S_{0}}$ as a rank function, and $f(\bm{C})=\|\bm{C}\|_{0}$ as a measure of sparsity in formulation (12), we obtain the following nonconvex optimization problem:

[TABLE]

The proximity operator $H:\mathbb{R}\rightarrow\mathbb{R}$ of $\|x\|_{0}$ is defined entry-wise as:

[TABLE]

The closed form solution of (36) at $y\in\mathbb{R}$ is the hard thresholding function defined in (2). The proximity operator of \big{\|}\bm{C}\big{\|}_{S_{0}} is the hard thresholding function applied entry-wise to the vector of singular values [60, 61].

IV-B Optimization Algorithm

To solve minimization problem in (35), we split original problem into two variables $\bm{J}$ and $\bm{C}$. That leads to the following objective function:

[TABLE]

The augmented Lagrangian function of (37) is:

[TABLE]

where $\mu$ is penalty parameter and $\bm{\Lambda}$ is Lagrange multiplier.

Update rule for $\bm{J}^{k+1}$: Given $\bm{C}^{k}$, $\bm{\Lambda}^{k}$, $\mu^{k}$, minimization of the Lagrangian function in (38) yields the following update:

[TABLE]

Update rule for $\bm{C}^{k+1}$: Given $\bm{J}^{k+1}$, $\bm{\Lambda}^{k}$, $\mu^{k}$, the following problem needs to be solved:

[TABLE]

When $\lambda\neq 0$ and $\tau=0$, the proximal map reduces to:

[TABLE]

Let $\bm{U\Sigma V}^{T}$ denote the SVD of matrix \big{(}\bm{J}^{k+1}+{\bm{\Lambda}}^{k}/\mu^{k}\big{)}. The closed-form solution of (41) is given by:

[TABLE]

where $H$ is the hard thresholding function defined in (2) and applied entry-wise to $\bm{\Sigma}$. Similarly, when $\lambda=0$ and $\tau\neq 0$ the proximal map is given by:

[TABLE]

Closed-form solution of (43) is obtained by the hard thresholding operator $H$ applied entry-wise to matrix \big{(}\bm{J}^{k+1}+{\bm{\Lambda}^{k}}/\mu^{k}\big{)}:

[TABLE]

Proximal average, introduced recently in [62] and generalized to nonconvex and nonsmooth setting in [63, 64], allows us to efficiently solve problem in (40) when $\lambda\neq 0$ and $\tau\neq 0$. In particular, given that the proximal maps $P_{f}^{\mu}$ and $P_{g}^{\mu}$ can be easily solved using the hard thresholding operator, we approximate the proximal map $P_{f+g}^{\mu}$ by averaging solutions of proximal maps of low-rank and sparse regularizers:

[TABLE]

where parameters $\tau$ and $\lambda$ are set such that $\tau+\lambda=1$.

Furthermore, since $\ell_{0}$ and $S_{0}$ norms belong to the class of semi-algebraic functions [73], the proximal average function $P_{f+g}$ is also a semi-algebraic function [63].

Update rule for Lagrange multiplier $\bm{\Lambda}^{t+1}$: Given $\bm{J}^{k+1}$, $\bm{C}^{k+1}$, $\mu^{k}$, Lagrange multiplier is updated with the following equation:

[TABLE]

The main steps of the proposed algorithm are summarized in Algorithm 2.

IV-C Convergence Analysis

Theorem 2

Let Y^{k}=\big{\{}\big{(}\bm{J}^{k},\bm{C}^{k},\bm{\Lambda}^{k}\big{)}\big{\}}_{k=1}^{\infty} be a sequence generated by Algorithm 2. Then, for any sufficiently large $\mu$, Algorithm 2 converges globally111That is, regardless of the initialization, it generates a bounded sequence that has at least one limit point which is a stationary point of (38)..

Proof:

The results in [63] guarantee convergence of the proximal average method. To guarantee global convergence of the Algorithm 2, we rewrite the problem (37) using the following more general form:

[TABLE]

where $\bm{A}=\bm{I}$, $\bm{B}=\bm{I}$, $f_{1}(\bm{C})=\lambda\|\bm{C}\|_{S_{0}}+\tau\|\bm{C}\|_{0}$, f_{2}(\bm{J})=\frac{1}{2}\big{\|}\bm{X}-\bm{XJ}\big{\|}^{2}_{F}.

We will now show that the assumptions A1-A5 in [70] which guarantee convergence in nonconvex nonsmooth optimization problem are satisfied. $\|\cdot\|_{0}$ and $\|\cdot\|_{S_{0}}$ are nonnegative lower semi-continuous functions and lower bounded. Therefore, $f_{1}$ as a sum of these functions is also lower semi-continuous and lower bounded. Furthermore, $f_{2}$ is coercive and $\bm{B}=\bm{I}$, so assumptions A1 and A4 hold. $\bm{A}=\bm{I}$ and $\bm{B}=\bm{I}$ imply that assumptions A2 and A3 hold. Next, $f_{2}$ is Lipschitz differentiable function so assumption A5 is also satisfied. Therefore, A1-A5 are satisfied and Algorithm 2 converges for any sufficiently large $\mu$ [70]. Of note, by splitting the original problem in (35) in three variables as done in GMC-LRSSC, we could not guarantee convergence since the assumption A2 would not be satisfied.

Furthermore, $\ell_{0}$ and $S_{0}$ norms belong to the class of semi-algebraic functions and satisfy Kurdyka-Łojasiewicz inequality [80, 73]. Sum of semi-algebraic function is again a semi-algebraic function, so \mathcal{L}_{\mu}\big{(}\bm{J},\bm{C},\bm{\Lambda}\big{)} in (38) is a semi-algebraic function and therefore, satisfies Kurdyka-Łojasiewicz inequality. This allows us to establish stronger convergence property, that is, sequence \big{\{}\big{(}\bm{J}^{k},\bm{C}^{k},\bm{\Lambda}^{k}\big{)}\big{\}} generated by Algorithm 2 converges regardless of the initialization to the unique limit point [70]. ∎

IV-D Stopping Criteria and Computational Complexity

The steps in Algorithm 2 are repeated until convergence or when the maximum number of iterations is exceeded. The convergence is achieved when inequalities \big{\|}\bm{J}^{k}-\bm{C}^{k}\|_{\infty}\leq\epsilon and \big{\|}\bm{J}^{k}-\bm{J}^{k-1}\|_{\infty}\leq\epsilon are satisfied. In all experiments error tolerance $\epsilon$ is set to $10^{-4}$.

As in GMC-LRSSC, the computational complexity of Algorithm 2 is $O(nN^{2}+TN^{3})$, where $T$ denotes the number of iterations. We set the maximal number of iterations $T$ to $100$, but the algorithm typically converged within $20$ iterations.

V Experimental Results

In this section, we compare the clustering performance and efficiency of the proposed algorithms with the state-of-the-art subspace clustering algorithms on synthetic and four real-world datasets. The performance is evaluated in terms of clustering error (CE) defined as:

[TABLE]

where $\Pi_{L}$ is the permutation space of $[L]$.

We compare the performance of our algorithms with the state-of-the-art subspace clustering algorithms, including Sparse Subspace Clustering (SSC) [21], Low-Rank Representation (LRR) [13, 14], closed form Low-Rank Subspace Clustering (LRSC) [17], Sparse Subspace Clustering via Orthogonal Matching Pursuit (SSC-OMP) [81], Thresholding based Subspace Clustering (TSC) [82], Nearest Subspace Neighbor (NSN) [83], Low-Rank Sparse Subspace Clustering (LRSSC) [29], $\ell_{0}$-Sparse Subspace clustering ($\ell_{0}$-SSC) [48, 84] and Schatten-$p$ norm minimization based LRR [51] ($S_{2/3}$-LRR and $S_{1/2}$-LRR).

V-A Experimental Setup

In all experiments, we set the parameters of GMC-LRSSC, and $S_{0}$/$\ell_{0}$-LRSSC as follows: $\tau=1-\lambda$, $\rho=3$, $\mu^{max}=10^{6}$, $\epsilon$ in stopping criteria to $10^{-4}$ and maximum number of iterations to $100$. Parameters $\lambda$ and initial value of $\mu$ are tuned more carefully. For GMC-LRSSC $\lambda$ is parameterized using $\alpha$ as $1/(1+\alpha)$, where $\alpha$ is tested in range $[10^{-3},10^{3}]$ with step $10$. Both $\lambda$ and $\tau$ are scaled by $\mu_{2}^{0}$. For $S_{0}/\ell_{0}$-LRSSC parameter $\lambda$ is optimized in range $[0.1,0.9]$ with step $0.1$. After the best $\lambda$ is found, $\mu_{2}^{0}$ in GMC-LRSSC and $\mu^{0}$ in $S_{0}/\ell_{0}$-LRSSC are tested in the set $\{1,3,5,10,20\}$. Initial value of parameter $\mu_{1}$ in GMC-LRSSC is set to $0.1$ in all experiments. For GMC-LRSSC we test nonconvexity parameter $\gamma\in\{0.1,0.6,1\}$. That resulted in $\gamma=1$ on the Extended Yale B dataset, $\gamma=0.6$ on the MNIST dataset and $\gamma=0.1$ on the USPS and ISOLET1 datasets. On synthetic data we test $\gamma$ from $0.1$ to $1$ with step $0.1$.

For other state-of-the-art algorithms, we use the source codes provided by the authors. If the best parameters are available, we set them as reported in the corresponding papers/source codes. Otherwise, we tuned the parameters and retained those with the best performance. Specifically, for SSC parameter $\alpha\in\{10,20,50,80,100,200,500,800,1000\}$, for LRR $\lambda\in\{0.05,0.1,0.3,0.5,1,2,3,4,5,6,7\}$, for LRSC $\tau\in\{0.1,5,10,20,50,80,100,200,500,800\}$ and $\alpha\in\{0.1\tau,0.5\tau,0.9\tau,1.11\tau,2\tau,10\tau\}$, for $S_{2/3}$-LRR and $S_{1/2}$-LRR $\lambda\in\{0.01,0.05,0.1,0.3,0.5,0.8,1,1.5,2,3,5,10\}$, and $\lambda$ in $\ell_{0}$-SSC is tuned in range $[0.1,1]$ with step $0.1$. In NSN and SSC-OMP number of neighbors is chosen in the set $\{2,3,5,8,10,12,15,18,20\}$. For TSC we set $q=max(3,\lceil n/20\rceil).$ For LRSSC we test $\lambda$ parameter in range $[10^{-3},10^{4}]$ with step $10$ on real-world datasets. In order to have completely same setting on the synthetic data, we tune the LRSSC parameters in the same way as for GMC-LRSSC.

Parameters of all algorithms are tuned on $20$ runs and for $L=\{3,5,10\}$ with different random seed than in the final experiment. In the final experiment we run each algorithm $100$ times.

V-B Synthetic Data

In the synthetic data experiment we compare LRSSC, GMC-LRSSC and $S_{0}/\ell_{0}$-LRSSC for different levels of noise and number of samples. We generate three $5$-dimensional disjoint subspaces embedded in the $100$-dimensional space. Subspace bases $\{U_{i}\}_{i=1}^{3}\in\mathbb{R}^{100\times 5}$ are constructed such that rank$([\bm{U}_{1},\bm{U}_{2},\bm{U}_{3}])=10$. We randomly sample $N_{i}$ data points from each subspace by computing $\{\bm{X}_{i}=\bm{U}_{i}\bm{A}_{i}\}_{i=1}^{3}$, where $\{\bm{A}_{i}\}_{i=1}^{3}\in\mathbb{R}^{5\times N_{i}}$ is generated from $\mathcal{N}(0,1)$ distribution. We sample the same number of data points from each subspace, i.e. $N_{1}=N_{2}=N_{3}$. We then add Gaussian noise with zero mean and vary the noise variance. Fig. 2 shows the average clustering error over $10$ runs for different number of samples per subspace and different noise variance.

For $50$ data points per subspace and small measurement noise, $S_{0}/\ell_{0}$-LRSSC performs better than LRSSC and GMC-LRSSC. On the other hand, for larger measurement noise GMC-LRSSC is the best performing algorithm. When we increase the number of data poins to $100$, GMC-LRSSC remains the best performing algorithm for most levels of noise. However, when further increasing number of data points, $S_{0}/\ell_{0}$ performs better except for very large measurement noise. This is in line with results presented in [56] which show that $\ell_{0}$ quasi-norm regularization of least-squares problems outperforms $\ell_{p}$ regularization $0<p\leq 1$ for small measurement noise. Whereas LRSSC and GMC-LRSSC in most cases do not improve performance when increasing the number of data points, $S_{0}/\ell_{0}$ is often able to exploit additional data.

V-C Face Recognition Dataset

The Extended Yale B dataset [85, 86] consists of face images of 38 individuals (subjects). It contains 64 frontal face images of each individual acquired under different illumination conditions. We use down-sampled $48\times 24$ pixel images and consider each vectorized image as one data point. The face images of each individual in Yale B dataset lie approximately in a 9-dimensional subspace [11].

We perform experiments for different number of clusters, ranging from 5 to 30. In each experiment we sample uniformly $L$ clusters from the total number of subjects and compute average of clustering error over 100 random subsets. The results are reported in Table II with the two best results highligted in boldface. Compared to the state-of-the-art methods, GMC-LRSSC and $S_{0}/\ell_{0}$-LRSSC achieve the lowest clustering error on all four clustering tasks. The difference between our algorithms and the second best performing other is significant for 10, 20 and 30 clusters (FDR<$1\%$; Benjamini-Hochberg corrected). Importantly, by increasing the number of clusters, the difference between our $\ell_{0}$-based formulations and convex low-rank and sparse formulations becomes larger.

To check the effect of parameter $\gamma$ in GMC-LRSSC, we vary parameter $\gamma$ from $0.1$ to $1$, where small $\gamma$ means that GMC penalty $\psi_{B}$ is close to convex, and $\gamma=1$ corresponds to maximally nonconvex value of the penalty. Fig. 3 shows performance as a function of $\gamma$ values for $10$, $20$ and $30$ clusters. On all these tasks, larger values of $\gamma$ achieve lower clustering error than the smaller values.

V-D Handwritten Digit Datasets

For handwriting recognition task, we consider two datasets: MNIST and USPS datasets. Both datasets contain pictures of ten digits (0-9), each digit corresponding to one cluster. The MNIST dataset contains $10000$ centered $28\times 28$ pixel images of handwritten digits. The USPS dataset consists of 92898 handwritten digit images, each of $16\times 16$ dimension. The handwritten digits lie approximately in a 12-dimensional subspace [87].

For both datasets, we use a subset of available images, sampling uniformly 50 images per digit in each run and compute average of clustering error over 100 runs. The performance comparisons for different choice of digits are shown in Table III. On both datasets, $S_{0}/\ell_{0}$-LRSSC and GMC-LRSSC are the only algorithms that consistently achieve high performance across varying combinations of digits.

On the MNIST dataset $S_{0}/\ell_{0}$-LRSSC algorithm is among the best performing algorithms, and the difference is increasing for larger number of clusters. GMC-LRSSC has lower performance than $S_{0}/\ell_{0}$-LRSSC for three combinations, but it is still significantly better than SSC, LRR and LRSSC. SSC and LRSSC have the best performance for digit sets $\{2,4,8\}$ and $\{2,4,6,8,9\}$, but they fail to give satisfactory results on other combinations of digits.

On the USPS dataset GMC-LRSSC is slightly better than $S_{0}/\ell_{0}$-LRSSC, except for the combinations of five digits. Specifically, on the digit set $\{2,4,6,8,9\}$ $S_{0}/\ell_{0}$-LRSSC outperforms all other methods. For the hardest problems with $10$ clusters GMC-LRSSC and $S_{0}/\ell_{0}$-LRSSC again have significantly better performance (FDR<$1\%$) than all other methods, that is $7.2\%$ and $4.8\%$ higher than the second best method, respectively. Fig. 4 illustrates derived affinity matrices on the USPS dataset for 10 clusters.

V-E Speech Recognition Dataset

For the speech recognition task, we evaluate algorithms on the ISOLET dataset [88]. The task is to cluster subjects, where each subject spoke the name of each letter of the alphabet twice. We use dataset ISOLET1 containing 26 subjects with 30 data points from each subject. The features include spectral coefficients, contour features, sonorant, presonorant and postsonorant features. To check whether the subspace clustering assumption holds on the ISOLET1 dataset, we compute the singular values of several subjects. Fig. 5 demonstrates that singular values decay rapidly and confirms that data points are drawn from low-dimensional subspaces.

We sample uniformly $L=\{5,10,15,20\}$ clusters from the total number of subjects over 100 random subsets. The average clustering errors are reported in Table IV. For all tested numbers of clusters $S_{0}/\ell_{0}$-LRSSC is the best performing method. GMC-LRSSC is the second best method for $5$ and $10$ clusters, while for $15$ and $20$ it achieves the same result as Schatten-$2/3$ and Schatten-$1/2$ LRR.

V-F Computational Time and Convergence

We further test the convergence behavior of GMC-LRSSC and $S_{0}/\ell_{0}$-LRSSC. The convergence conditions of GMC-LRSSC are satisfied within less than 15 iterations on all four real-world datasets. Fig. 6 illustrates convergence behavior of GMC-LRSSC on the MNIST and ISOLET1 datasets for 10 clusters. $S_{0}/\ell_{0}$-LRSSC converges within 20 iterations on the Extended Yale B, MNIST and USPS datasets. Fig. 7 shows behavior of $S_{0}/\ell_{0}$-LRSSC on the MNIST and ISOLET1 datasets for 10 clusters. Although the maximal number of iterations is exceeded on the ISOLET1 dataset, Fig. 7 shows that the error decays rapidly and within 20 iterations.

The average computational time over 100 runs of each algorithm is shown in Table V. On all datasets, LRSC is consistently the fastest algorithm. GMC-LRSSC is among the fastest algorithms, whereas $S_{0}/\ell_{0}$-LRSSC is among the fastest algorithms on the Extended Yale B and MNIST datasets. As explained above, on ISOLET1 dataset $S_{0}/\ell_{0}$-LRSSC exceeds maximum number, resulting in higher computational time. All experiments were run in MATLAB 2017a environment on the PC with $256$ GB of RAM and Intel Xeon CPU E5-2650 v4 $2$ processors operating with a clock speed of $2.2$ GHz.

VI Conclusion and Discussion

In this paper we have introduced two nonconvex regularizations in low-rank sparse subspace clustering: (i) multivariate generalization of minimax-concave penalty function; and (ii) $S_{0}/\ell_{0}$ regularization with a proximal operator approximated using the proximal average method. Under the proposed framework, we have presented two algorithms based on the alternating direction method of multipliers. We have performed extensive experiments on synthetic and four real world datasets, including face recognition, handwriting recognition and speech recognition tasks. Our experimental results have shown that both proposed methods converge fast and achieve clustering error lower than nuclear and $\ell_{1}$ norms regularized objective. Moreover, for larger number of clusters proposed methods consistently outperform existing subspace clustering methods. That is explained by their more accurate approximation of rank and sparsity of the data representation matrix than it is the case with nuclear and $\ell_{1}$ norms.

The choice of norm should be decided depending on the dataset. In future work we plan to study how different initializations affect the accuracy of the proposed nonconvex regularization based methods. Instead of directly solving NP hard $\ell_{0}$ quasi-norm minimization problem, we are interested in finding a way to gradually build a solution. A possible strategy could be to start with $\ell_{1}$ norm solution. Since analytic formulas for thresholding function exist for $p=1/2$ and $p=2/3$ [89], we can gradually shrink $p$ and use the current solution to initialize the next step of the algorithm. This approach is called $p$-continuation strategy in [56].

[Proof of Proposition 2 and Theorem 1]

In this section, we first prove the boundedness of variables in Algorithm 1. This result helps us to establish convergence property of Algorithm 1. We then prove Theorem 1 in the paper, where we show that any converging point satisfies Karush-Kuhn-Tucker (KKT) conditions [78].

Proof of Proposition 2:

From the first order optimality conditions of Lagrangian function in (25) we have:

[TABLE]

The optimality condition of problem in (30) implies that:

[TABLE]

where \big{[}\partial\psi_{\bm{B}}\big{(}\bm{C}_{2}^{k+1}\big{)}\big{]}_{ij} denotes gradient of the GMC penalty $\psi_{\bm{B}}$ at \big{[}\bm{C}_{2}^{k+1}\big{]}_{ij}. By the definition of the scaled MC penalty in (22) and using $\bm{B}=\sqrt{\mu_{2}^{k}\gamma/\tau}\bm{I}$ we have:

[TABLE]

where $c_{ij}$ denotes \big{[}\bm{C}_{2}\big{]}_{ij}.

If $|c_{ij}|>\tau/(\mu_{2}^{k}\gamma)$, then from (50) and (51) directly follows \big{[}\bm{\Lambda}_{2}^{k+1}\big{]}_{ij}=0.

Otherwise, we get the following equality:

[TABLE]

Since \big{|}\big{[}\bm{C}_{2}^{k+1}\big{]}_{ij}\big{|}\leq\tau/(\mu_{2}^{k}\gamma), it follows \big{|}\big{[}\bm{\Lambda}_{2}^{k+1}\big{]}_{ij}\big{|}\leq 1. Therefore, sequence \big{\{}\bm{\Lambda}_{2}^{k}\big{\}} is bounded.

The optimality condition of problem in (28) implies that:

[TABLE]

Similarly, following the proof for $\bm{C}_{2}$ and using Proposition 1, it can be shown that the sequence \big{\{}\bm{\Lambda}_{1}^{k}\big{\}} is also bounded. Using the definitions of $\bm{J}^{k+1}$, $\bm{C}_{1}^{k+1}$ and $\bm{C}_{2}^{k+1}$ as minimizers, we have the following inequalities [90]:

[TABLE]

Note that the last term equals:

[TABLE]

where $a_{i}=\frac{\mu_{i}^{k}+\mu_{i}^{k-1}}{2{(\mu_{i}^{k-1})}^{2}}$, $i=1,2$. Since $\mu^{k}=\rho\mu^{k-1}$, $\rho>1$, it follows that $\mu^{k}$ is non-decreasing and $\sum_{k=0}^{\infty}\frac{1}{\mu_{i}^{k}}<\infty$, $i=1,2$. We then have:

[TABLE]

Since \big{\{}\bm{\Lambda}_{1}^{k}\big{\}} and \big{\{}\bm{\Lambda}_{2}^{k}\big{\}} are bounded, then \big{\{}\|\bm{\Lambda}_{1}^{k}-\bm{\Lambda}_{1}^{k-1}\|_{F}^{2}\big{\}} and \big{\{}\|\bm{\Lambda}_{1}^{k}-\bm{\Lambda}_{1}^{k-1}\|_{F}^{2}\big{\}} in (55) are also bounded. Therefore, from (54), (55) and (56) follows that \big{\{}\mathcal{L}_{\mu_{1}^{k},\mu_{2}^{k}}\big{(}\bm{J}^{k+1},\bm{C}_{1}^{k+1},\bm{C}_{2}^{k+1},\bm{\Lambda}_{1}^{k},\bm{\Lambda}_{2}^{k}\big{)}\big{\}} is upper-bounded.

Furthermore, it holds that:

[TABLE]

Since \big{\{}\mathcal{L}_{\mu_{1}^{k-1},\mu_{2}^{k-1}}\big{(}\bm{J}^{k},\bm{C}_{1}^{k},\bm{C}_{2}^{k},\bm{\Lambda}_{1}^{k-1},\bm{\Lambda}_{2}^{k-1}\big{)}\big{\}} is upper-bounded and \big{\{}\bm{\Lambda}_{1}^{k}\big{\}} and \big{\{}\bm{\Lambda}_{2}^{k}\big{\}} are both bounded, \big{\{}\bm{J}^{k}\big{\}} is also bounded.

The boundedness of \big{\{}\bm{C}_{1}^{k}\big{\}} and \big{\{}\bm{C}_{2}^{k}\big{\}} follows from the update rules for Lagrange multipliers in $(\ref{Lambda_rule})$ and boundedness of \big{\{}\bm{J}^{k}\big{\}}, \big{\{}\bm{\Lambda}_{1}^{k}\big{\}}, and \big{\{}\bm{\Lambda}_{2}^{k}\big{\}}.

Therefore, sequence \big{\{}\big{(}\bm{J}^{k},\bm{C}_{1}^{k},\bm{C}_{2}^{k},\bm{\Lambda}_{1}^{k},\bm{\Lambda}_{2}^{k}\big{)}\big{\}} generated by Algorithm 1 is bounded. Bolzano-Weierstrass theorem [91] then guarantees the existence of a convergent subsequence.

Proof of Theorem 1:

Let \big{(}\bm{J}^{*},\bm{C}_{1}^{*},\bm{C}_{2}^{*},\bm{\Lambda}_{1}^{*},\bm{\Lambda}_{2}^{*}\big{)} be a critical point of (24). The KKT conditions are derived as follows:

[TABLE]

From the 1st and 4th KKT conditions, it follows:

[TABLE]

Let $\bm{U}_{1}\bm{\Sigma}_{1}\bm{V}_{1}^{T}$ be the SVD of matrix $\bm{C}_{1}^{*}$. Using Proposition 1, the right hand side of (59) equals:

[TABLE]

where $a_{1}=\mu_{1}^{*}/\lambda$ and $b_{1}=\sqrt{\mu_{1}^{*}\gamma/\lambda}$, $0<\gamma\leq 1$. The scalar function $Q_{a,b}$ is defined as $Q_{a,b}(x)\triangleq x+\frac{1}{a}\partial\phi_{b}(x)$, where $\phi_{b}$ is the scaled MC penalty defined in (22). $Q_{a_{1},b_{1}}(x)$ is applied element-wise to singular values of matrix $\bm{C}_{1}^{*}$.

Let $\bm{U}_{2}\bm{\Sigma}_{2}\bm{V}_{2}^{T}$ be the SVD of matrix \big{(}\bm{J}^{*}+\bm{\Lambda}_{1}^{*}/\mu_{1}^{*}\big{)}. From (59), we get the following relation:

[TABLE]

It is easy to verify that $Q_{a_{1},b_{1}}$ is a monotone function [92, 93] and Q_{a_{1},b_{1}}^{-1}(x)=\Theta\big{(}x;\frac{\lambda}{\mu_{1}^{*}},\frac{\lambda}{\gamma\mu_{1}^{*}}\big{)} for $0<\gamma\leq 1$ , where $\Theta$ is the firm thresholding function defined in (1). Applying $Q_{a_{1},b_{1}}^{-1}$ to both sides of (61) and replacing $\bm{U}_{1}\bm{\Sigma}_{1}\bm{V}_{1}^{T}=\bm{C}_{1}^{*}$, we get:

[TABLE]

Similarly, from the 2nd and 5th KKT conditions, we have the following relations:

[TABLE]

where $a_{2}=\mu_{2}^{*}/\tau$ and $b_{2}=\sqrt{\mu_{2}^{*}\gamma/\tau}$, $0<\gamma\leq 1$.

Again, applying $Q_{a_{2},b_{2}}^{-1}$ to both sides of (63), we get the following equations:

[TABLE]

Therefore, the 4th and 5th KKT conditions can be rewritten as:

[TABLE]

We now show that KKT conditions are satisfied when assumptions of Theorem 1 hold. From (33) we have:

[TABLE]

Since by the assumption \big{(}\bm{\Lambda}_{1}^{k+1}-\bm{\Lambda}_{1}^{k}\big{)}\rightarrow\bm{0} and \big{(}\bm{\Lambda}_{2}^{k+1}-\bm{\Lambda}_{2}^{k}\big{)}\rightarrow\bm{0}, then the first two KKT conditions are satisfied.

From the update rule in (27), we have:

[TABLE]

From the the first two conditions, it follows that when $\bm{J}^{k+1}-\bm{J}^{k}\rightarrow\bm{0}$, the 3rd KKT condition is satisfied.

Next, using the update for $\bm{C}_{1}$ in (29) we obtain the following equation:

[TABLE]

where $\bm{U}_{2}\bm{\Sigma}_{2}{\bm{V}_{2}}^{T}$ is the SVD of matrix \big{(}\bm{J}^{k+1}+{\bm{\Lambda}_{1}}^{k}/\mu_{1}^{k}\big{)}.

From the update rule for $\bm{C}_{2}$ in (32) it follows:

[TABLE]

When $\bm{C}_{1}^{k+1}-\bm{C}_{1}^{k}\rightarrow\bm{0}$ and $\bm{C}_{2}^{k+1}-\bm{C}_{2}^{k}\rightarrow\bm{0}$, follow the equations in (65). Since \big{\{}Y^{k}\big{\}}_{k=1}^{\infty} is bounded and equations (66), (67), (68), (69) go to zero, we conclude that the sequence \big{\{}Y^{k}\big{\}}_{k=1}^{\infty} asymptotically satisfies the KKT conditions in (58).

Acknowledgment

The authors would like to thank Shin Matsushima for his helpful comments on improving the statements of Proposition 2 and Theorem 1. The authors would also like to thank the Editor and reviewers for their comments which helped to improve this paper.