Block Diagonal Least Squares Regression for Subspace Clustering

Abstract

Least squares regression (LSR) is an effective method that has been widely used for subspace clustering. Under the conditions of independent subspaces and noise-free data, coefficient matrices can satisfy enforced block diagonal (EBD) structures and achieve good clustering results. More importantly, LSR produces closed solutions that are easier to solve. However, solutions with block diagonal properties that have been solved using LSR are sensitive to noise or corruption as they are fragile and easily destroyed. Moreover, when using actual datasets, these structures cannot always guarantee satisfactory clustering results. Considering that block diagonal representation has excellent clustering performance, the idea of block diagonal constraints has been introduced into LSR and a new subspace clustering method, which is named block diagonal least squares regression (BDLSR), has been proposed. By using a block diagonal regularizer, BDLSR can effectively reinforce the fragile block diagonal structures of the obtained matrices and improve the clustering performance. Our experiments using several real datasets illustrated that BDLSR produced a higher clustering performance compared to other algorithms.

1. Introduction

High-dimensional data are increasingly common in everyday life, which has greatly affected the methods that are used for data analysis and processing. Performing clustering analyses on high-dimensional data has become one of the hot spots in current research. It is usually assumed that high-dimensional data are distributed within a joint subspace, which is composed of multiple low-dimensional subspaces. This supposition has become the premise of the study of subspace clustering [1].

Over the last several decades, researchers have proposed many useful approaches for solving clustering problems [2,3,4,5,6,7]. Among these methods, spectral-type subspace clustering has become the most popular for use in subspace algorithms because of its outstanding clustering effects when processing high-dimensional data. The spectral-type method mainly utilizes local or global information from data samples to construct affinity matrices. In general, the representation matrices greatly affect the performance of spectral-type subspace clustering. In recent years, a vast number of different methods have been proposed in the attempt to find better representation matrices. Two typical representative methods in linear spectral-type clustering are sparse subspace clustering (SSC) [2,8] and low-rank representation (LRR) [3,9], which are based on linear representation. They emphasize sparseness and low rank when solving coefficient matrices.

Initially, SSC used the L₀ norm to represent the sparsity of the matrices. When considering that the optimization of the L₀ norm is NP-hard, the L₀ norm was finally replaced by the L₁ norm. SSC pursues the sparse representation by minimizing the L₁ norm of coefficient matrices, which can better classify data. However, the obtained matrices are too sparse and miss the important correlation structures of the data. Thus, SSC leads to unsatisfactory clustering results for highly correlated data. To make up for the shortcomings of SSC, LRR employs the nuclear norm to search for low-rank representations with the goal of better analyzing the relevant data. By describing global structures, LRR improves the clustering performance of highly correlated data. However, the coefficient matrices that are obtained using the LRR model may not be sparse enough. To this end, many improved methods have been proposed.

Luo et al. [10] made use of the prior of sparsity and low rank and then presented a multi-subspace representation (MSR) model. Zhuang et al. [11] improved the MSR model by adding non-negative constraints to the coefficient matrices and then proposed a non-negative low-rank sparse representation (NNLRS) model. Both the MSR and NNLRS models ensure that coefficient matrices have both sparsity and grouping effects by using a combination of sparsity and low rank. To explore local and global structures, Zheng et al. [12] incorporated local constraints into LRR and presented a locally constrained low-rank representation (LRR-LC) model. Chen et al. [13] developed another new model by adding symmetric constraints to the LRR so that highly correlated data in the subspaces had consistent representation. To better reveal the essential structures of subspaces, Lu et al. [14] utilized the F norm to constrain the coefficient matrices and proposed the use of least squares regression (LSR) for subspace clustering. This method demonstrated that enforced block diagonal (EBD) conditions were helpful for block diagonalization, i.e., under EBD conditions, the structures of the obtained coefficient matrices are block diagonal. In the meantime, LSR encouraged the grouping effect by clustering highly correlated data together. It has been theoretically proven that LSR tends to shrink the coefficients of related data and is robust to noise. More importantly, the objective function of LSR produces closed solutions, which are easy to solve and reduce the time complexity of the algorithm. Generally, under the independent subspace assumption, the models that use the L₁ norm, nuclear norm or F norm as regularization terms can satisfy EBD conditions. This can help them to learn affinity matrices with block diagonal properties.

Assuming that the data samples are not contaminated or damaged and that the subspaces are independent of each other, the above methods can obtain coefficient matrices with block diagonal structures. However, the structures that are obtained using these methods are unstable because the required assumptions are often not established. When a fragile block diagonal structure is broken, it causes the clustering performance to decline to a great extent. For this reason, many subspace algorithms have been proposed that extend the idea of block diagonal structures. Zhang et al. [15] introduced an effective structure into LRR and proposed a recognition method to distinguish block diagonal low-rank representation (BDLRR), which learned distinguished data representations while shrinking non-block diagonal parts and highlighting block diagonal representations to improve the clustering effects. Feng et al. [16] devised two new subspace clustering methods, which were based on a graph Laplace constraint. Using the above methods, block diagonal coefficient matrices can be precisely constructed. However, their optimization processes are inefficient because each updating cycle contains an iterative projection. Hence, the hard constraints of the coefficient matrices are difficult to satisfy precisely. To settle this matter, Lu et al. [17] used block diagonal regularization to add constraints into the solutions and then put forward a block diagonal representation (BDR) model. These constraints, which were more flexible than direct hard constraints, could help the BDR model to easily obtain solutions with block diagonal properties. Therefore, the regularization was called a soft constraint. This block diagonal regularization caused BDR to be non-convex, but it could be solved via simple and valid means. This model imposed block diagonal structures on the coefficient matrices, which was more likely to lead to accurate clustering results. Due to the good performance of the k-diagonal block within subspace clustering, numerous corresponding extended algorithms [18,19,20] have been proposed.

Under the conditions of independent subspaces and noise-free data, LSR can obtain coefficient matrices that have block diagonal properties, which usually produce exact clustering results. However, the coefficient matrices that are obtained using the LSR method are sensitive to noise or corruption. In real environments, data samples are often contaminated by noise. In this case, the fragile block diagonal structures can become damaged. When the block diagonal structures are violated, the LSR method does not necessarily lead to satisfactory clustering results.

Inspired by the BDR method, we constructed a novel block diagonal least squares regression (BDLSR) subspace clustering method by combining LSR and BDR. Considering that the LSR method meets EBD conditions and has a grouping effect, we directly pursued block diagonal matrices on the basis of LSR. The block diagonal regularizer was a soft constraint that encouraged the block diagonalization of the coefficient matrices and could improve the accuracy of the clustering. However, the appearance of the block diagonal regularizer resulted in the non-convexity of BDLSR. Fortunately, the alternating minimization method could be used to solve the objective function. Meanwhile, the convergence of the objective function could also be proven. Our experimental results from using several public datasets indicated that BDLSR was more effective.

The rest of the paper is organized as follows. The least squares regression model and the block diagonal regularizer are briefly reviewed in Section 2. The BDLSR model is proposed, optimized and discussed explicitly in Section 3. In Section 4, the experimental effectiveness of BDLSR is demonstrated using three datasets. Finally, the conclusions are summarized.

2. Notations and Preliminaries

2.1. Notations

Many notations were used in this study, which are shown in

Table 1. The notations that were used in this study.

Notation	Meaning
A	Any matrix, represented by capital letters
I	Identity matrix
a	Any vector, represented by boldface lowercase letters
1	The all one vector
Ai,:	The i-th row of A
A:j	The j-th column of A
Ai,j	The entry in the i-th row and j-th column of A
diag(A)	A vector with its i-th entry being the i-th diagonal element of A
Diag(a)	A diagonal matrix with its i-th entry on the diagonal being the i-th element of a
$A\succcurlyeq 0$	Positive semidefinite
$A\preccurlyeq B$	$B-A\succcurlyeq 0$
$A\ge 0$	All entries of A are non-negative
Tr(A)	The trace of a square matrix
$<A,B>$	The inner product between A and B
${\left[A\right]}_{+}$	$\max (A,0)$
$\Vert A\Vert =\sqrt{{\displaystyle \sum _{ij}{a}_{ij}^2}}$	The F norm of A

.

2.2. Least Squares Regression

In real data, most data have strong correlations. However, SSC enforces solutions to be sparse and LRR causes the solutions to have low-rank representation; therefore, these methods still cannot obtain the intrinsic structures of the data very well. To better explore these correlations, Lu et al. [14] obtained coefficient matrices using the F norm constraint and proposed a least squares regression (LSR) subspace algorithm. LSR has a grouping effect, which makes it helpful for the better clustering of highly correlated data.

In a given data matrix Y, the data samples are from multiple independent linear subspaces. To determine the clustering number of the subspaces, the data samples need to be divided into their respective subspaces. Ideally, the data are clean. By solving the following optimization problems, LSR can solve the coefficient matrix Z:

$$\min \Vert Z\Vert ,\begin{array}{cc}& s.t.\begin{array}{cc}Y=YZ,& diag(Z)=0\end{array}\end{array}$$

where, $\Vert Z\Vert$ denotes the F norm of Z, which is defined as $\Vert Z\Vert =\sqrt{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{Z}_{ij}^2}}}$. Note that the constraint $diag(Z)=0$ in the above model can help to avoid trivial solutions. In other words, each data sample is only represented by itself and then the case of Z being the identity matrix can be avoided.

When the data contain noise, a noise term is added to enhance the robustness of the model. The extended function of LSR is shown below:

$$\min {\Vert Y-YZ\Vert }^2+\lambda {\Vert Z\Vert }^2,\begin{array}{cc}& s.t.\begin{array}{cc}diag(Z)=0& \end{array}\end{array}$$

(1)

where $\lambda >0$ is a regularization parameter, which is used to balance the two terms within the above objective function. Note that LSR produces analytical solutions, which are easy to solve. Due to these advantages, LSR has attracted the attention of a growing number of researchers.

2.3. The Block Diagonal Regularizer

Existing studies [2,3,14] have shown that block diagonal structures can better reveal the essential attributes of subspaces. Ideally, under the assumption of independent subspaces and noise-free data, the obtained coefficient matrix Z has a k-block diagonal structure, where k represents the number of blocks or subspaces. Then, Z can be expressed as follows:

$$Z=\left[\begin{array}{cccc}{Z}_1& 0& \dots & 0\\ 0& Z_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & Z_k\end{array}\right]$$

Here, each block represents a subspace, the size of blocks represents the dimensions of the subspaces and the data within the same block belong to the same subspace.

However, due to noise in the data, the assumption of independent subspaces is usually incorrect. Hence, the solutions of the model are always far from the k-block diagonal structure. To obtain the above structure, a block diagonal regularizer was proposed in [17]:

$$||B|{|}_{{}_{\boxed{K}}}={\displaystyle \sum _{i=n-k+1}^n{\lambda }_i(L_B)},$$

(2)

where $B\in {ℝ}^{n\times n}$ satisfies the constraints of $B\ge 0$ and $B=B^T$ and is an affinity matrix, $\left|\right|B|{|}_{{}_{\boxed{K}}}$ denotes the k-block diagonal regularizer (which is defined as the sum of the smallest k eigenvalues of $L_B$) and $L_B$ denotes the corresponding Laplacian matrix (which is defined as $L_B=Diag(B1)-B$, where B1 represents the product of matrix B and vector 1). Here, ${\lambda }_i(L_B)$ denotes the i-th eigenvalue of $L_B$ when arranged in decreasing order. Due to $L_B\succcurlyeq 0$, $L_B$ is positive semidefinite and ${\lambda }_i(L_B)\ge 0$ for all values of i. More precisely, B has a k-block diagonal when (and only when):

$${\lambda }_i(L_B)\{\begin{array}{cc}>0& i=1,\dots ,n-k,\\ =0& i=n-k+1,\dots ,n.\end{array}$$

The term in Equation (2) is a soft constraint, which does not directly require exact solutions with block diagonal structures. However, the number of blocks can be controlled and correct clustering results can still be obtained. Hence, the regularizer is more flexible and produces accurate clustering results more easily than hard constraints [16].

3. Subspace Clustering Using BDLSR

In this section, our BDLSR method, which is a combination of least squares regression and block diagonal representation, is proposed for subspace clustering. Then, the objective function of the BDLSR model is presented and the optimization process is described. Finally, the convergence of the BDLSR model is proven.

3.1. The Proposed Model

In theory, under the independent subspace assumption, LSR can obtain coefficient matrices that have block diagonal attributes. In other words, when the similarity between classes is zero and the similarity within classes is more than zero, perfect data clustering can usually be achieved. However, in real environments, the block diagonal structures of the solutions can be easily broken by data noise or corruptions. Considering the advantages of block diagonal representation, we integrated it into the objective function of LSR to obtain solutions with strengthened block diagonal structures, which required the matrices to be simultaneously non-negative and symmetric. Then, we proposed the BDLSR method for subspace clustering. The objective function of BDLSR for handling noisy data was as follows (where $\alpha >0,\gamma >0$ and Y is a data matrix):

$$\begin{array}{l}\underset{B}{\min }\frac{1}{2}||Y-YB|{|}^2+\frac{\alpha }{2}{\Vert B\Vert }^2+\gamma ||B|{|}_{{}_{\boxed{K}}},\\ s.t.diag(B)=0,B\ge 0,B=B^T.\end{array}$$

(3)

where the coefficient matrix B is non-negative and symmetric and its diagonal elements are all zero. The above conditions had to be satisfied to use a block diagonal regularizer. However, the conditions limited the representation capability of B. To alleviate this issue, we needed to relax the restrictions on B. Thus, we introduced an intermediary term Z into the above formula. The new model was rewritten equivalently as follows (where $\lambda >0$):

$$\begin{array}{l}\underset{Z,B}{\min }\frac{1}{2}||Y-YZ|{|}^2+\frac{\alpha }{2}||Z|{|}^2+\frac{\lambda }{2}||Z-B|{|}^2+\gamma ||B|{|}_{{}_{\boxed{K}}},\\ s.t.diag(B)=0,B\ge 0,B=B^T.\end{array}$$

(4)

where $\alpha ,\lambda$ and $\gamma$ are the balance parameters. When $\lambda$ was sufficiently large, Equation (4) was equivalent to Equation (3). The new added item $\frac{\lambda }{2}||Z-B|{|}^2$ was strongly convex for updating Z and B; thus, it was easy to solve closed-form solutions and perform convergence analyses. Then, we could alternatively optimize the variables Z, V and B.

3.2. The Optimization of BDLSR

Next, we solved the minimization problem of the BDLSR model in Equation (4). Because of the block diagonal regularizer, it was observed that Equation (4) was non-convex. The key to solving this issue was the non-convex term $||B|{|}_{{}_{\boxed{K}}}$. Fortunately, we could utilize the related theorem in [17,21] to reformulate $||B|{|}_{{}_{\boxed{K}}}$. The theorem in [17] was presented as follows.

We let $L\in {ℝ}^{n\times n},L\succcurlyeq 0$. Then:

$${\displaystyle \sum _{i=n-k+1}^n{\lambda }_i}(L)=\underset{V}{\min }<L,V>,s.t.0\preccurlyeq V\preccurlyeq I,Tr(V)=k.$$

where $<L,V>=Tr(L^{T}V)$.

Then, $||B|{|}_{{}_{\boxed{K}}}$ could be reformulated as a convex problem:

$$||B|{|}_{{}_{\boxed{K}}}=\underset{V}{\min }<L_B,V>,s.t.0\preccurlyeq V\preccurlyeq I,Tr(V)=k.$$

Due to $L_B$ being defined as $Diag(B1)-B$, Equation (4) was equivalent to the following problem:

$$\begin{array}{l}\underset{Z,B,V}{\min }\frac{1}{2}||Y-YZ|{|}^2+\frac{\alpha }{2}||Z|{|}^2+\frac{\lambda }{2}||Z-B|{|}^2+\gamma <Diag(B1)-B,V>,\\ s.t.diag(B)=0,B\ge 0,B=B^T,0\preccurlyeq V\preccurlyeq I,Tr(V)=k.\end{array}$$

(5)

The above problem involved three variables, which could be obtained by alternating the minimization. Then, we decomposed the problem into three subproblems and developed the update steps separately.

By fixing the variables B and V, we could update Z using:

$$Z^{m+1}=\arg \underset{Z}{\min }\frac{1}{2}||Y-YZ|{|}^2+\frac{\alpha }{2}||Z|{|}^2+\frac{\lambda }{2}||Z-B|{|}^2,$$

(6)

Note that Equation (6) was a convex program. By setting the partial derivative of Z to 0, we could obtain its closed-form solution. The solution was:

$$Z^{m+1}={(Y^{T}Y+\lambda I+\alpha I)}^{-1}(\lambda B+Y^{T}Y)$$

(7)

Similarly, by fixing the variables B and Z, we could update V using:

$$\begin{array}{l}{V}^{m+1}=\arg \underset{V}{\min }<Diag(B1)-B,V>,\\ s.t.0\preccurlyeq V\preccurlyeq I,Tr(V)=k.\end{array}$$

(8)

According to Equation (8), we could obtain the solution for V as follows:

$$V^{m+1}=U{U}^T,$$

(9)

where $U\in {ℝ}^{n\times k}$ consists of k eigenvectors, which are associated with the k smallest eigenvalues of $Diag(B1)-B$.

Then, we could solve B. By fixing the variables Z and V, we could update B using:

$$\begin{array}{l}{B}^{m+1}=\arg \underset{B}{\min }\frac{\lambda }{2}||Z-B|{|}^2+\gamma <Diag(B1)-B,V>,\\ s.t.diag(B)=0,B\ge 0,B=B^T.\end{array}$$

(10)

By converting Equation (10), the following function was obtained:

$$\begin{array}{l}{B}^{m+1}=\arg \underset{B}{\min }\frac{1}{2}||B-Z+\frac{\gamma }{\lambda }(diag(V)1^T-V)|{|}^2,\\ s.t.diag(B)=0,B\ge 0,B=B^T.\end{array}$$

(11)

According to the proposition in [17], we could obtain the closed-form solution of B:

$$B^{m+1}={{\displaystyle \left[\frac{1}{2}(A-Diag(diag(A)))+\frac{1}{2}{(A-Diag(diag(A)))}^T\right]}}_{+}$$

(12)

where:

$$\begin{gathered} A=Z-\frac{\gamma }{\lambda }(diag(V)1^T-V), \\ {[A]}_{+}=\{\begin{array}{c}{a}_{ij},\begin{array}{cc}if& a_{ij}>0\end{array}\\ 0,\begin{array}{cc}if& a_{ij}\le 0\end{array}\end{array} \end{gathered}$$

The complete optimization procedure for the BDLSR model for subspace clustering is summarized in Algorithm 1.

Algorithm 1 The optimization process for solving Equation (5).

Input: Data matrix Y, $\lambda >0$, $\alpha >0, \gamma >0$.

Initialize:$m=0,Z^m=0,V^m=0,B^m=0$.

while not converge do

1: Update $Z^{m+1}$ by solving (6).

2: Update $V^{m+1}$ by solving (8).

3: Update $B^{m+1}$ by solving (11).

4: m = m + 1.

end while

Output:$Z^{*},V^{*},B^{*}$.

3.3. Discussion

To solve Equation (4), we developed an iterative update algorithm, which is shown in this section. For a given data matrix, we could settle the proposed BDLSR problem by using representation matrices Z and B, which could be obtained using Algorithm 1. Then, we used Z and B to construct the affinity matrices, e.g., $P=(|B|+|B^T|)/2$ or $P=(|Z|+|Z^T|)/2$. Finally, we performed spectral-type subspace clustering on the constructed affinity matrices.

We noted that the BDLSR algorithm contained three parameters, which could be utilized to balance the different terms in Equation (4). We could first run BDR for the suitable $\lambda$ and $\gamma$ values or run LSR for the suitable $\alpha$ and $\lambda$ values. After that, we could search the last remaining parameter to achieve the highest clustering accuracy.

Note that the number of subspaces had to be known to construct the affinity matrices. Considering this requirement is required for all spectral-type subspace algorithms, the number always had to be known. Our proposed BDLSR algorithm used the block diagonal structure and its clustering effectiveness was shown in the following experiments.

3.4. Convergence Analysis

Because of the block diagonal regularizer, the objective function in Equation (4) was non-convex. Fortunately, Equation (4) could be solved by alternating the minimization and its solution could be proven to converge to the global optimum. The alternating minimization method is summarized in Algorithm 1. Now, we explain the convergence of Algorithm 1 theoretically.

To guarantee the convergence of BDLSR, the sequence $\left\{Z^m,V^m,B^m\right\}$ that was generated by Algorithm 1 needed to be proven to have at least one limit point. In addition, any limit point $(Z^{*},V^{*},B^{*})$ of $\left\{Z^m,V^m,B^m\right\}$ was a stationary point of Equation (5). Now, with regard to Algorithm 1, we detail the convergence analysis by referring to [17].

The objective function of Equation (5) was denoted as $G(Z^m,V^m,B^m)$, where m represents the m-th iteration. We let $S_V=\{V|0\preccurlyeq V\preccurlyeq I,Tr(V)=k\}$ and $S_B=\{B|diag(B)=0,B\ge 0,B=B^T\}$. The indicator functions of $S_V$ and $S_B$ were denoted as $g_{S_V}(V)$ and $g_{S_B}(B)$, respectively.

First, according to the update process of $Z^{m+1}$ in Equation (6), we obtained:

$$Z^{m+1}=\arg \underset{Z}{\min }G(Z,V^m,B^m)$$

Note that $G(Z,V^m,B^m)$ is $\lambda$-strongly convex with regard to Z. We obtain

$$G(Z^{m+1},V^m,B^m)\le G(Z^m,V^m,B^m)-\frac{\lambda }{2}{\Vert Z^{m+1}-Z^m\Vert }^2$$

(13)

Then, due to the optimality of $V^{m+1}$ in Equation (8), we obtained:

$$G(Z^{m+1},V^{m+1},B^m)+g_{S_V}(V^{m+1})\le G(Z^{m+1},V^m,B^m)+g_{S_V}(V^m)$$

(14)

Next, according to the update process of $B^{m+1}$ in Equation (11), we obtained:

$$B^{m+1}=\arg \underset{B}{\min }G(Z^{m+1},V^{m+1},B)+g_{S_B}(B).$$

Note that $G(Z^{m+1},V^{m+1},B)+g_{S_B}(B)$ was strongly $\lambda$-convex with regard to B. We obtained:

$$G(Z^{m+1},V^{m+1},B^{m+1})+g_{S_B}(B^{m+1})\le G(Z^{m+1},V^{m+1},B^m)+g_{S_B}(B^m)-\frac{\lambda }{2}{\Vert B^{m+1}-B^m\Vert }^2$$

(15)

By combining and simplifying (13)–(15), we obtained:

$$\begin{array}{cc}\begin{array}{l}G(Z^{m+1},V^{m+1},B^{m+1})+g_{S_V}(V^{m+1})+g_{S_B}(B^{m+1})\\ \le G(Z^m,V^m,B^m)+g_{S_V}(V^m)+g_{S_B}(B^m)-\frac{\lambda }{2}{\Vert Z^{m+1}-Z^m\Vert }^2-\frac{\lambda }{2}{\Vert B^{m+1}-B^m\Vert }^2\end{array}& \end{array}$$

(16)

Hence, $G(Z^m,V^m,B^m)+g_{S_V}(V^m)+g_{S_B}(B^m)$ monotonically decreased and was upper bounded. This implied that $\left\{Z^m\right\}$ and $\left\{B^m\right\}$ were bounded. Additionally, $V^m\in S_V$ implied that ${\Vert V^m\Vert }_2\le 1$ and $\left\{V^m\right\}$ were bounded.

Since $V^m$ and $Diag(B^{m}1)-B^m$ were positive semidefinite, we obtained $\begin{array}{cc}<Diag(B^{m}1)-B^m,V^m>& \ge 0\end{array}$. Thus, $G(Z^m,V^m,B^m)+g_{S_V}(V^m)+g_{S_B}(B^m)\ge 0$. Then, by summing Equation (16) over m = 0, 1, 2, …, we obtained:

$${\displaystyle \sum _{m=0}^{+\infty }\frac{\lambda }{2}}\left({\Vert Z^{m+1}-Z^m\Vert }^2+{\Vert B^{m+1}-B^m\Vert }^2\right)\le G(Z^0,V^0,B^0)$$

This implied:

$$Z^{m+1}-Z^m\to 0$$

(17)

and

$$B^{m+1}-B^m\to 0$$

(18)

Using Equation (18) and the update process of $V^{m+1}$ in Equation (8), we obtained:

$$V^{m+1}-V^m\to 0$$

(19)

Then, since $\left\{Z^m\right\}$, $\left\{V^m\right\}$ and $\left\{B^m\right\}$ were bounded, a point $\{Z^{*},V^{*},B^{*}\}$ and a subsequence $\{Z^{mi},V^{mi},B^{mi}\}$ existed such that $Z^{mi}\to Z^{*}$, $V^{mi}\to V^{*}$ and $B^{mi}\to B^{*}$. Then, according to Equations (17)–(19), we obtained $Z^{mi+1}\to Z^{*}$, $V^{mi+1}\to V^{*}$ and $B^{mi+1}\to B^{*}$. On the other hand, by considering the optimality of $V^{mi+1}$ in Equation (8), $Z^{mi+1}$ in Equation (6) and $B^{mi+1}$ in Equation (10), we obtained:

$$\begin{array}{l}0\in \nabla G_Z(Z^{mi+1},V^{mi},B^{mi}),\\ 0\in \nabla G_V(Z^{mi+1},V^{mi+1},B^{mi})+{\partial }_V{g}_{S_V}(V^{mi+1}),\\ 0\in \nabla G_B(Z^{mi+1},V^{mi+1},B^{mi+1})+{\partial }_B{g}_{S_B}(B^{mi+1}).\end{array}$$

(20)

By letting $m\to +\infty$ in Equation (20), we obtained:

$$\begin{array}{l}0\in \nabla G_Z(Z^{*},V^{*},B^{*}),\\ 0\in \nabla G_V(Z^{*},V^{*},B^{*})+{\partial }_V{g}_{S_V}(V^{*}),\\ 0\in \nabla G_B(Z^{*},V^{*},B^{*})+{\partial }_B{g}_{S_B}(B^{*}).\end{array}$$

(21)

Thus, $(Z^{*},V^{*},B^{*})$ was a stationary point of Equation (5).

4. Experiments

4.1. Experimental Settings

4.1.1. Data

To evaluate the BDLSR model, we performed experiments using different datasets: Hopkins 155 [22], ORL [23] and Extended Yale B [17].

The Hopkins 155 dataset contains 156 different sequences, which were extracted from moving objects. Some samples from the Hopkins 155 dataset are shown in Figure 1.

The ORL dataset, which was created by Cambridge University, contains 400 facial photos of 40 subjects with a dark uniform background under different conditions for time, light and facial expression. We resized each photo to 32 × 32. Some samples from the ORL dataset are shown in Figure 2.

The Extended Yale B dataset is composed of images of 38 individuals and can be used for facial recognition. The dataset contains 2414 images, with 64 photos of each individual under different conditions and angles. We resized each photo to 32 × 32. Some samples from the Extended Yale B dataset are shown in Figure 3.

4.1.2. Baseline

We chose representative subspace clustering algorithms for comparison to verify the effectiveness of BDLSR, including SSC [8], LRR [9], LSR [14] and BDR-B(BDR-Z) [17]. These algorithms are spectral-type algorithms and construct affinity matrices by applying coefficient matrices to achieve the final results.

SSC (sparse subspace clustering) is based on sparse representation, which ensures the sparsity of the coefficient matrices by minimizing the L₁ norm. LRR (low-rank representation) is a typical subspace algorithm that searches for low-rank representations of coefficient matrices using the nuclear norm. LSR (least squares regression) is an easy-to-solve clustering subspace algorithm that uses the F norm to constrain coefficients so that there is a grouping effect between the coefficients. BDR (block diagonal representation) is a novel spectral clustering subspace algorithm that uses a block diagonal regularizer, which provides two methods for solving coefficient matrices: BDR-B and BDR-Z.

In this paper, we chose AC and NMI [24] as the evaluation criteria to quantitatively evaluate BDLSR. AC was used to evaluate the accuracy of the clustering and NMI was used to evaluate the correlations between the clustering results and the real results. The larger the AC and NMI values, the better the clustering results.

We supposed that $k_i$ and $k_i^{*}$ were the ground truth label and the obtained label of a data sample $y_i$, respectively, to define AC is as follows (where n is the number of data samples):

$$\mathrm{AC}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\delta (k_i,k_i^{*}})$$

(22)

where the delta function $\delta (a,b)$ = 1 when (and only when) a = b.

We supposed that the two clustering results were D and D’ to define NMI as follows:

$$\mathrm{NMI}(D,D^{\prime })=\frac{\mathrm{MI}(D,D^{\prime })}{\max (\mathrm{H}(D),\mathrm{H}(D^{\prime }))}$$

(23)

where $H(·)$ and $MI(D,D^{\prime })$ represent the entropy and the mutual information of the clustering results, respectively.

All algorithms were implemented using the MATLAB2016a program on Windows 10.

4.2. Experimental Results

Here, we discuss the clustering effects of BDLSR during our experiments using the Hopkins 155, ORL and Extended Yale B datasets. The comparison algorithms in our experiments were SSC, LRR, LSR and BDR-B (BDR-Z). The code was provided by the original authors to ensure the accuracy of the algorithms. The parameters were adjusted to the optimal values, according to the settings in the original papers.

In our experiments, the specific parameter adjustment process was as follows. We first determined the initial value ranges of the parameters $\lambda$ and $\gamma$, according to the parameter settings in the BDR algorithm. Then, we used the parameter values of $\lambda$ and $\gamma$ that were provided in [17] as the initial values to try for another parameter $\alpha$ and determined the value range of this parameter according to the changes in the clustering results. Finally, the grid search method was used for an exhaustive search of the value ranges of the three determined parameters and the values of the parameters were selected when the BDLSR algorithm produced the best performance. In this paper, $\alpha \in [{10}^{-3},{10}^{-6}]$, $\lambda \in [{10}^0,{10}^{-3}]$ and $\gamma \in [{10}^{-1},{10}^{-4}]$ were selected for Hopkins 155, $\alpha \in [{10}^0,{10}^{-3}]$, $\lambda \in [{10}^0,{10}^{-3}]$ and $\gamma \in [{10}^1,{10}^{-2}]$ were selected for ORL and $\alpha \in [{10}^{-2},{10}^{-5}]$, $\lambda \in [{10}^2,{10}^{-1}]$ and $\gamma \in [{10}^1,{10}^{-2}]$ were selected for Extended Yale B.

The Hopkins 155 dataset contains 156 motion sequences, which are classified into groups of two and three motions. In this experiment, we used the mean and median values of AC to evaluate the different algorithms. We adopted two settings to construct data matrices: 2F-dimensional and 4k-dimensional data points.

Table 2. The accuracy (%) when using the Hopkins 155 dataset with 2F-dimensional data.

Algorithm	2 Motions		3 Motions		All Motions
	Mean	Median	Mean	Median	Mean	Median
SSC	98.48	100.00	95.60	98.37	97.82	100.00
LRR	96.35	99.78	90.60	96.01	95.05	99.47
LSR	96.76	100.00	94.06	97.95	96.15	99.55
BDR-B	99.00	100.00	98.05	99.79	98.78	100.00
BDR-Z	99.05	100.00	99.15	99.79	99.07	100.00
BDLSR-B	99.19	100.00	98.87	99.79	99.12	100.00
BDLSR-Z	98.95	100.00	98.32	99.79	98.81	100.00

and

Table 3. The accuracy (%) when using the Hopkins 155 dataset with 4k-dimensional data.

Algorithm	2 Motions		3 Motions		All motions
	Mean	Median	Mean	Median	Mean	Median
SSC	99.17	100.00	95.60	99.44	97.59	100.00
LRR	95.78	99.71	90.57	96.30	94.60	99.47
LSR	96.65	99.71	93.87	97.95	96.03	99.47
BDR-B	98.74	100.00	98.78	99.79	98.75	100.00
BDR-Z	98.96	100.00	98.78	99.80	98.92	100.00
BDLSR-B	98.76	100.00	97.73	99.79	98.53	100.00
BDLSR-Z	98.98	100.00	98.95	99.81	98.97	100.00

show the AC values of each algorithm using the Hopkins 155 dataset. The partial comparative data in Table 1 and Table 2 originated from [17]. In the experiment, we set $\alpha =4\times {10}^{-4}$, $\lambda =9\times {10}^{-2}$, $\gamma =9\times {10}^{-3}$ for 2F-dimensional data and $\alpha =6.5\times {10}^{-5}$, $\lambda =7\times {10}^{-2}$, $\gamma =7\times {10}^{-3}$ for 4k-dimensional data. Using the selected parameter settings, BDLSR could perform well in most cases. From Table 2 and Table 3, it can be seen that the AC value of BDLSR was higher than those of the other algorithms under the two different settings, except for the three motions in Table 2. These results showed that BDLSR had some advantages over other representative clustering algorithms when using the Hopkins 155 dataset.

The ORL dataset includes 40 subjects, with 10 photos of each. We selected the first k (10, 20, 30, 40) subjects of the dataset for the face clustering experiments. The comparison results from using the ORL dataset are shown in

Table 4. The results from the different algorithms when using the ORL dataset (%).

Number of Subjects	Algorithm	AC	NMI
10	SSC	72.00	77.09
	LRR	63.60	70.06
	LSR	76.00	78.11
	BDR-B	75.00	76.81
	BDR-Z	78.00	79.44
	BDLSR-B	76.00	75.65
	BDLSR-Z	85.00	83.04
20	SSC	66.50	82.55
	LRR	63.50	76.23
	LSR	73.50	84.07
	BDR-B	76.00	85.74
	BDR-Z	74.00	83.35
	BDLSR-B	82.50	86.43
	BDLSR-Z	86.50	89.35
30	SSC	70.00	84.56
	LRR	65.90	78.77
	LSR	76.67	84.41
	BDR-B	79.67	88.95
	BDR-Z	75.33	87.21
	BDLSR-B	86.33	91.36
	BDLSR-Z	81.33	88.36
40	SSC	72.75	84.56
	LRR	66.30	81.07
	LSR	73.25	70.30
	BDR-B	65.75	83.20
	BDR-Z	74.75	87.53
	BDLSR-B	70.25	85.36
	BDLSR-Z	81.00	88.68

. We set different parameters for the different numbers of subjects in the experiments. For instance, we set $\alpha =9\times {10}^{-2}, \lambda =7\times {10}^{-2}, \gamma =1.67$ for 40 subjects. From Table 4, it can be seen that BDLSR improved the AC value from 74.75% to 81.00% and the NMI value from 87.53% to 88.68% when using BDLSR-Z. This showed that BDLSR produced the best clustering results out of all of the algorithms for 10, 20, 30 and 40 subjects.

The Extended Yale B dataset includes 2414 photos of 38 individuals. We selected the first k (3, 5, 8, 10) individuals of the dataset to compare the AC and NMI values of the different algorithms. The comparison results from using the Extended Yale B dataset are shown in

Table 5. The results from the different algorithms when using the Extended Yale B dataset (%).

Number of Subjects	Algorithm	AC	NMI
3	SSC	95.31	87.34
	LRR	92.71	79.07
	LSR	90.10	80.14
	BDR-B	96.35	89.32
	BDR-Z	95.31	87.34
	BDLSR-B	95.31	87.34
	BDLSR-Z	98.96	95.11
5	SSC	94.69	90.39
	LRR	90.94	79.43
	LSR	88.75	78.58
	BDR-B	97.13	93.49
	BDR-Z	95.56	92.27
	BDLSR-B	98.44	95.27
	BDLSR-Z	96.88	93.20
8	SSC	88.48	83.50
	LRR	60.55	61.48
	LSR	78.91	67.22
	BDR-B	95.90	91.85
	BDR-Z	92.19	86.85
	BDLSR-B	96.68	93.21
	BDLSR-Z	92.38	87.46
10	SSC	81.56	76.90
	LRR	64.38	65.72
	LSR	76.25	65.70
	BDR-B	93.59	89.82
	BDR-Z	73.13	71.37
	BDLSR-B	94.06	90.10
	BDLSR-Z	73.59	72.52

. We set $\alpha =1\times {10}^{-4}, \lambda =11, \gamma =1$ for 10 subjects in these experiments. From Table 5, it can be seen that BDLSR improved the AC value from 93.59% to 94.06% and the NMI value from 89.82% to 90.10% when using BDLSR-B. This showed that BDLSR achieved better results than the other algorithms for 3, 5, 8 and 10 subjects.

From the results, it can be seen that the BDR algorithm produced good clustering results. However, the proposed BDLSR algorithm had higher AC and NMI values than the BDR algorithm. This demonstrated that the combination of LSR and BDR could enhance the grouping effect on related data and improve the clustering results.

4.3. Parameter Analysis

The BDLSR algorithm mainly involved three parameters $\alpha , \lambda$ and $\gamma$. Now, we further discuss the effects of the parameter settings on the performance of BDLSR. Considering that BDLSR mainly used block diagonal regularization and least squares regression to improve the clustering effect, we mainly discuss the influence of the parameters $\alpha$ and $\gamma$ on the AC and NMI values. Then, we present the clustering results that were achieved by fixing $\lambda$ and setting different values for $\alpha$ and $\gamma$.

Figure 4 displays the AC values when using the Hopkins 155 dataset with different parameter settings. For 2F dimensional data, $\lambda$ = 0.09 and $\alpha =\{3\times {10}^{-4},4\times {10}^{-4},5\times {10}^{-4},6\times {10}^{-4}\}$. For 4k dimensional data, $\lambda$ = 0.07 and $\alpha =\{6.2\times {10}^{-5},6.5\times {10}^{-5},6.8\times {10}^{-5},7.0\times {10}^{-5}\}$. For both settings, $\gamma =\{7\times {10}^{-3},8\times {10}^{-3},9\times {10}^{-3},1\times {10}^{-2}\}$. By taking the results of all motions in the Hopkins 155 dataset as the example, we could observe the changes in the AC values that were caused by the BDLSR-B method for 2F-dimensional data and those that were caused by the BDLSR-Z method for 4k-dimensional data. From Figure 4, it can be seen that the AC values fluctuated very little with the changes to α and $\gamma$. When $\alpha =4\times {10}^{-4}$, $\gamma =9\times {10}^{-3}$ and $\alpha =6.5\times {10}^{-5}$, $\gamma =7\times {10}^{-3}$, the clustering effect was the best.

Figure 5 displays the AC and NMI values when using the ORL dataset with different parameter settings. We set $\lambda$ = 0.07, $\alpha =\{7\times {10}^{-2},8\times {10}^{-2},9\times {10}^{-2},1\times {10}^{-1}\}$ and $\gamma =\{1.61,1.64,1.67,1.70\}$. To observe the changes in the AC and NMI values that were caused by the BDLSR-Z method, we took 40 individuals from the ORL dataset as the example. From Figure 5, it can be seen that the AC values changed significantly and the NMI values changed slightly within a certain value range. When $\alpha =9\times {10}^{-2}$ and $\gamma =1.67$, BDLSR-Z produced the best clustering effect.

Figure 6 displays the AC and NMI values when using the Extended Yale B dataset with different parameter settings, where λ = 50, $\alpha =\{2\times {10}^{-4},5\times {10}^{-4},7\times {10}^{-4},9\times {10}^{-4}\}$ and $\gamma =\{0.7,0.8,0.9,1.0\}$. To observe the changes in the AC and NMI values that were caused by the BDLSR-B method, we took 5 individuals from the Extended Yale B dataset as the example. From Figure 6, it can be seen that the AC and NMI values were relatively stable. On the whole, changes in the parameters had little impact on the clustering effect.

4.4. Ablation Studies

In this section, we reveal the influences of various parameters on the clustering results of the algorithm using ablation studies. Specifically, the original parameter settings remained unchanged when using the three standard datasets, then $\alpha =0$, $\lambda =0$ and $\gamma =0$ were set to perform the ablation studies. The comparison results are shown in

Table 6. A comparison of the accuracy (%) results when using the Hopkins 155 dataset with 2F-dimensional data in the ablation studies.

Algorithm	2 Motions		3 Motions		All Motions
	Mean	Median	Mean	Median	Mean	Median
BDLSR-B ($\alpha =0$)	98.69	100.00	98.02	99.80	98.53	100.00
BDLSR-Z ($\alpha =0$)	98.90	100.00	99.06	99.79	98.94	100.00
BDLSR-B ($\lambda =0$)	69.03	67.09	56.87	56.30	66.29	64.95
BDLSR-Z ($\lambda =0$)	91.26	99.19	89.13	98.00	90.78	98.95
BDLSR-B ($\gamma =0$)	98.82	100.00	97.91	99.72	98.61	100.00
BDLSR-Z ($\gamma =0$)	98.20	100.00	96.49	99.33	97.81	100.00
BDLSR-B	99.19	100.00	98.87	99.79	99.12	100.00
BDLSR-Z	98.95	100.00	98.32	99.79	98.81	100.00

,

Table 7. A comparison of the accuracy (%) results when using the Hopkins 155 dataset with 4k-dimensional data in the ablation studies.

Algorithm	2 Motions		3 Motions		All Motions
	Mean	Median	Mean	Median	Mean	Median
BDLSR-B ($\alpha =0$)	98.74	100.00	98.78	99.79	98.75	100.00
BDLSR-Z ($\alpha =0$)	98.96	100.00	98.78	99.80	98.92	100.00
BDLSR-B ($\lambda =0$)	69.01	67.09	56.92	56.51	66.28	64.95
BDLSR-Z ($\lambda =0$)	90.24	98.44	84.71	93.16	88.99	98.25
BDLSR-B ($\gamma =0$)	98.43	100.00	96.49	99.72	97.99	100.00
BDLSR-Z ($\gamma =0$)	98.04	100.00	96.20	99.33	97.99	100.00
BDLSR-B	98.76	100.00	97.73	99.79	98.53	100.00
BDLSR-Z	98.98	100.00	98.95	99.81	98.97	100.00

,

Table 8. A comparison of the results from the different algorithms when using the ORL dataset (%).

Number of Subjects	Algorithm	AC	NMI
40	BDLSR-B ($\alpha =0$)	68.00	83.74
	BDLSR-Z ($\alpha =0$)	74.75	86.86
	BDLSR-B ($\lambda =0$)	3.50	20.45
	BDLSR-Z ($\lambda =0$)	75.75	85.99
	BDLSR-B ($\gamma =0$)	74.00	85.36
	BDLSR-Z ($\gamma =0$)	74.75	86.28
	BDLSR-B	70.25	85.36
	BDLSR-Z	81.00	88.68

and

Table 9. A comparison of the results from the different algorithms when using the Extended Yale B dataset in the ablatio studies (%).

Number of Subjects	Algorithm	AC	NMI
10	BDLSR-B ($\alpha =0$)	93.91	89.94
	BDLSR-Z ($\alpha =0$)	73.59	72.74
	BDLSR-B ($\lambda =0$)	10.16	2.77
	BDLSR-Z ($\lambda =0$)	11.72	4.51
	BDLSR-B ($\gamma =0$)	90.16	85.44
	BDLSR-Z ($\gamma =0$)	69.53	65.73
	BDLSR-B	94.06	90.10
	BDLSR-Z	73.59	72.52

.

From the comparative data in Table 6, Table 7, Table 8 and Table 9, it can be seen that after the ablation studies, the clustering effects when using the three standard datasets decreased significantly (except the three motions in the Hopkins 155 dataset with 2F-dimensional data). Our experimental results showed that each component of the objective function of the BDLSR algorithm played an effective role. This also showed that the combination of block diagonal structures and LSR could effectively improve the performance of subspace clustering.

5. Conclusions

In this work, to obtain the advantages of both the LSR and BDR algorithms, block diagonal structures were introduced into an LSR algorithm. Then, the BDLSR method was proposed for subspace clustering. The algorithm made use of a regularization item to enhance the fragile block diagonal structures that were obtained using LSR and could produce better clustering results and demonstrate a strong competitiveness. After that, we settled the objective function using the alternating minimization method. In addition, we also theoretically proved the convergence of this algorithm and verified its effectiveness using three standard datasets. Our experimental results showed that the BDLSR method achieved the best clustering results compared to other algorithms in almost all cases when using the three standard datasets. However, the algorithm also had some shortcomings: it involved more parameters (three parameters), which needed more time to be adjusted; additionally, the clustering results from a certain dataset were not the best out of the compared algorithms. To solve these shortcomings, in our future work, we will explore the concept of parameter reduction and how to set the minimum number of parameters while ensuring the best clustering accuracy so as to obtain a clustering algorithm with an even better performance.