Robust Multi-subspace Analysis Using Novel Column L0-norm Constrained Matrix Factorization
Binghui Wang Chuang Lin
cs.CV


























 































































































Robust Multi-subspace Analysis Using Novel Column L_0-norm Constrained Matrix Factorization
    Binghui Wang, Student Member, IEEE,
       and Chuang Lin, Member, IEEE
Binghui Wang is with the Electrical and Computer Engineering Department, Iowa State University, Ames, IA, 50010 USA, e-mail: binghuiw@iastate.edu.
Chuang Lin is with the Shenzhen Institute of Advanced Technology, Chinese Academy of Science, Shenzhen, China, e-mail: chuang.lin@siat.ac.cn.
Binghui Wang and Chuang Lin contributed equally to this work.
The preliminary work was done in Dalian University of Technology, China, and was further extended and completed in Iowa State University, USA.
Preliminary results of the proposed method have been published in <cit.>.
    December 30, 2023
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





We study the underlying structure of data (approximately) generated from a union of independent subspaces. Traditional methods learn only one subspace, failing to discover the multi-subspace structure, while state-of-the-art methods analyze the multi-subspace structure using data themselves as the dictionary, which cannot offer the explicit basis to span each subspace and are sensitive to errors via an indirect representation. Additionally, they also suffer from a high computational complexity, being quadratic or cubic to the sample size. 

To tackle all these problems, we propose a method, called Matrix Factorization with Column L_0-norm constraint (MFC_0), that can simultaneously learn the basis for each subspace, generate direct sparse representation for each data sample, as well as removing errors in the data in an efficient way. 
Furthermore, we develop a first-order alternating direction algorithm, whose computational complexity is linear to the sample size, to stably and effectively solve the nonconvex objective function and nonsmooth l_0-norm constraint of MFC_0. 

Experimental results on both synthetic and real-world datasets demonstrate that besides the superiority over traditional and state-of-the-art methods for subspace clustering, data reconstruction, error correction, MFC_0 also shows its uniqueness for multi-subspace basis learning and direct sparse representation.


















Robust multi-subspace analysis, matrix factorization, multi-subspace basis learning, l_0-norm sparse representation, alternating direction algorithm, proximal algorithm.


















§ INTRODUCTION


The observed data are extremely high dimensional in this “big data" era. Typically, these data reside in a much lower-dimensional latent subspace, instead of being uniformly distributed in the high-dimensional ambient space. Thus, it is of great importance to reveal the underlying structure of the data as it helps to reduce the computational cost and enables a compact representation for learning. Such an idea has been successfully applied to various reseach communities, e.g., dimensionality reduction <cit.>, face recognition <cit.>, metric learning <cit.>, etc..

Subspace methods have been widely used to analyze the data, and linear subspace[A vector space is a subset of some other higher-dimensional vector space.] is the most common choice for its simplicity and computational efficiency. Additionally, linear subspace has shown its effectiveness in modeling real-world problems such as motion segmentation <cit.>, face clustering <cit.>, and handwritten digits recognition <cit.>. Consequently, subspace analysis has been paid much attention in the past decades. For instance, principal component analysis (PCA) aims to learn a subspace while retaining maximal variances of the data. Nonnegative matrix factorization (NMF) <cit.> is designed to learn both the nonnegative basis and nonnegative parts-based representation.

Robust PCA (RPCA) <cit.> assumes that the data are approximately drawn from a low-rank subspace while perturbed by sparse noise. The basic assumption of these methods is the single subspace, which is not the case in many practical applications. A more reasonable way is to consider the data as lying on or near a union of linear subspaces. Unfortunately, the generalization to handle multiple subspaces is quite challenging.

Recently, multi-subspace analysis has attracted increasingly interests in visual data analysis <cit.>. Derived from recent advances in compressive sensing <cit.>, the methods including 

<cit.> 
have incorporated sparse and/or low-rank regularization into their formulations to model the mixture of linear structures for clean data[Data points are strictly sampled from the respective subspace.] as well as dealing with errors in data, e.g., noise <cit.>, missed entries <cit.>, corruptions <cit.>, and outliers <cit.>. Among them, sparse subspace clustering (SSC) <cit.> and low-rank representation (LRR) <cit.> stand out as two most popular methods, which formulate the discovery of multi-subspace structure as finding a sparse or low-rank representation of the data samples using data themselves as the dictionary (basis). It has been shown in literatures that these methods can achieve more accurate subspace structure than single subspace analysis methods. However, both SSC and LRR have the following major drawbacks: First, they use the original data as the basis rather than learning the basis explicitly. It is problematic when data contain errors yet are still used for reconstruction or/and clustering. Second, neither sparse representation via l_1 regularization nor low-rank representation in a global constraint provides a direct description for each data sample. In other words, data samples cannot find the subspace where they lie, and their contained errors are unable to be removed. Thirdly, SSC and LRR suffer from the computational complexity that is quadratic and cubic to the sample size, respectively. 



Our work: In this paper, to tackle all above problems, we aim to simultaneously learn the basis for each subspace and directly generate sparse representation for each data sample, as well as removing errors (e.g., random corruptions[A fraction of random entries of data are grossly corrupted.] and sample-specific outliers[A fraction of data are far away from their respective subspaces.]) in the data in an efficient way. To this end, we propose a novel column L_0-norm constrained matrix factorization (MFC_0) method, and develop a first-order alternating direction algorithm to stably and efficiently solve the nonsmooth and nonconvex objective function of MFC_0. 

Specifically, given a collection of data approximately generated from a union of independent subspaces with an equal subspace dimension, we explore their subspace structure from the matrix factorization perspective in the following way: 

First, we restrict the basis matrix to be orthonormal, which guarantees that the learnt basis indeed span multiple subspaces.
Then, we pursue a direct sparse representation of each data sample on only the basis that form its underlying subspace. 


To achieve this goal, we leverage an l_0-norm to confine the number of nonzeros elements of each column of the representation matrix. 
Furthermore, to handle different types of errors contained in the data, we incorporate different regularizations. 



With above constraints, the objective function of MFC_0 is nonconvex on the coupled basis and representation matrix, and is nonsmooth on the l_0-norm constraint. Thus, it is a very challenging optimization problem that cannot be solved via typical gradient (or subgradient) based methods. To overcome the nonconvex problem, we develop a first-order optimization algorithm motivated by ADMM <cit.> to split multivariable objective function into several univariate subproblems and solve each subproblem one by one. To tackle the subproblem on nonsmooth l_0-norm constraint, we design a novel proximal operator, which is inspired by recent proximal algorithms <cit.>, and solve it efficiently and analyticly. 

After learning basis and representation matrix, MFC_0 can accomplish many tasks, e.g., subspace clustering, data reconstruction, error correction, etc.. We evaluate MFC_0 on both synthetic data and real-world datasets. Experimental results demonstrate that besides the superiority and efficiency over traditional and state-of-the-art methods for subspace clustering, data reconstruction, and error correction, MFC0 also demonstrates its uniqueness for multi-subspace basis learning and direct representation learning. 




In summary, our key contributions are as follows:


  * We propose a novel MFC_0 method from the matrix factorization perspective to perform multi-subspace analysis. To the best of our knowledge, this is the first work that considers simultaneously learning the basis, generating direct sparse representation, and correcting errors of data generated from multiple subspaces. 

  * We develop a first-order alterating direction algorithm to stably and effectively solve the nonconvex objective function of MFC_0. 

  * We design a novel proximal operator to analyticly and efficiently solve the nonsmooth l_0-norm constraint.



  * Experimental results demonstrate that MFC_0 shows high resistance to errors,  

achieves high subspace clustering performance, learns multi-subspace basis, and generates direct sparse representation.

























§ RELATED WORK


Matrix factorization and its variants
Matrix Factorization (MF) has been studied for several decades. Classical methods, such as Principal Component Analysis (PCA), Linear Discriminant Analysis (LDA), Independent Component Analysis (ICA), are exemplars of low-rank matrix factorizations. Afterwards, to solve specific problems, many MF variants, including Concept MF <cit.>, Maximum Margin MF <cit.>, Convex MF <cit.>, and Online MF <cit.>, etc., are proposed. These methods share in common that there is no constraint on the sign of elements of factorized matrices. In contrast, a paradigm of factorization, called NMF, was first proposed in <cit.> and has been widely used for its sparse, parts-based, and additive representation. NMF is suitable for many research fields, e.g., data mining, image processing, pattern recognition, and machine learning. Its real-world applications include document clustering <cit.>, image segmentation <cit.>, image recognition <cit.>, speech and audio processing <cit.>, blind source separation <cit.>, to name a few. One can refer to the comprehensive review of NMF in <cit.>.

Approximate sparsity via l_1-norm 
Sparsity in classic NMF occurs is as a by-product rather than a designed objective. To this end, researchers have paid much attention to impose sparsity on NMF explicitly. These methods either penalize <cit.> or constrain <cit.> l_1-norm of the representation matrix to yield a sparse representation <cit.>. As a convex relaxation of l_0-norm, the exact sparsity measure which renders NP-hard combinational optimization problem, l_1-norm is used to improve computational feasibility and efficiency of sparse representation <cit.>. However, as NMF is NP-hard itself <cit.>, we can say that any algorithm for NMF is suboptimal. Thus, a heuristic l_0-norm constrained NMF might be as appropriate and efficient as l_1-norm sparse NMF. 

L_0-norm based sparsity 
Little work has been done using direct l_0-norm, the exact sparsity measure which causes the NP-hard combinational optimization problem. K-SVD <cit.> aims to find an overcomplete dictionary for sparse representation, and the sparse approximation is based on an l_0-norm inequality constraint and is solved via greedy orthogonal matching pursuit (OMP) <cit.>. 

Peharz and Pernkopf <cit.> further proposed sparse NMF that introduced an l_0-norm to constrain the basis or representation matrix of approximate NMF. It is solved via sparse alternating nonnegative least-squares (sANLS).
Wang et al. <cit.> leveraged hierarchical Bayes to model an adaptive sparseness prior, which has a similar properity with an l_0-norm.
However, these methods are under the hidden assumption that data are generated from single subspace. In addition, they cannot handle various errors. In this work, we focus on analyzing data (approximately) drawn from multiple subspaces. Different from using greedy methods, the l_0-norm constraint in MFC_0 is solved analytically via our designed proximal operator, which is easy to be incorporated into our developed alternating direction algorithm. 



§ NOTATIONS AND PROBLEM DEFINITION




 §.§ Main Notations

In this paper, vectors and matrices are written as bold lowercase and uppercase symbols. The i-th entry of vector 𝐮 is u_i. For matrix 𝐖, its (i,j)-th entry, i-th column, and j-th row are denoted as w_i,j, 𝐰_i, and 𝐰^j respectively. 
The horizontal and vertical concatenation of a set of K matrices 𝐖_1, 𝐖_2, ⋯, 𝐖_K along row and column are denoted by [𝐖_1, 𝐖_2, ⋯, 𝐖_K ] and [𝐖_1; 𝐖_2; ⋯; 𝐖_K ], respectively.
The superscript T stands for the transpose. 

For vector 𝐮∈^m, its l_p-norm is 𝐮_p = (∑_i=1^m |u_i|^p)^1/p, and pseudo l_0-norm is the number of nonzero entries, i.e., 𝐮_0 = #{i:u_i ≠ 0}. For matrix 𝐖∈^m × n, its Frobenius norm is 𝐖_F =
√(∑_j=1^m 𝐰^j_2^2); l_2,1-norm is 𝐖_2,1=∑_i=1^n 𝐰_i _2; l_1-norm is 𝐖_1=∑_j=1^m ∑_i=1^n |w_ij|; and l_∞-norm is 𝐖_∞ = max_i,j|w_i,j|. The inner product between two matrices 𝐖, 𝐔 is ⟨𝐖, 𝐔⟩ = tr(𝐖^T 𝐔), where tr is the trace operator defined on a square matrix.

We denote k-th subspace as 𝒮_k and a collection of K subspaces as {𝒮_k}_k=1^K.



 §.§ Problem Definition


We present two definitions related to our problem at first. 


  * Union and sum of subspaces The union of K subspaces {𝒮_k}_k=1^K is defined as ∪_k=1^K 𝒮_k = {𝐭: 𝐭∈𝒮_k, ∀ k }. The sum of K subspaces {𝒮_k}_k=1^K is defined as ∑_k=1^K 𝒮_k = {𝐭: 𝐭 = ∑_k=1^K 𝐭_k, 𝐭_k ∈𝒮_k }.    


  * Independent subspaces K subspaces {𝒮_k}_k=1^K are independent if and only if 𝐒_k ∩∑_j ≠ k𝒮_j = { 0 }. Or to say, dim(∑_k=1^K 𝒮_k) = ∑_k=1^K dim (𝒮_k). 


With above definitions, our problem is defined as
[Robust Multi-subspace Analysis]
Given a collection of data samples strictly/approximately drawn from a union of independent subspaces with an equal subspace dimension, robust multi-subspace analysis is to learn the basis for each subspace, learn the associated representation for each sample upon the basis, or/and correct the error for each data sample.


Formally, given n data samples 𝐙 = {𝐳_i}_i=1^n ∈^m × n that are generated from a basis matrix 𝐗̃ = {𝐱̃_i}_i=1^d ∈^m × d that span K independent subspaces {𝒮_k}_k=1^K, with a corresponding representation matrix 𝐘̃ = {𝐲̃_i}_i=1^n ∈^d × n, or/and contaminated with errors 𝐄̃∈^m × n 

    𝐙 = 𝐗̃𝐘̃ + 𝐄̃.


Then, our purpose is from the observed data 𝐙 to learn the basis 𝐗̃_k for each subspace 𝒮_k, the representation 𝐘̃_k upon the basis 𝐗̃_k, or/and the error 𝐄̃_k such that


    [ 𝐙_1, ⋯, 𝐙_K ] = [ 𝐗̃_1, ⋯, 𝐗̃_K ] 
    [
         [ 𝐘̃_1;    ⋮; 𝐘̃_K;      ]] 
    + [ 𝐄̃_1, ⋯, 𝐄̃_K ],


where 𝐙= [ 𝐙_1, ⋯, 𝐙_K ] with 𝐙_k ∈^m × n_k merely collecting the samples from the k-th subspace 𝒮_k and ∑_k=1^K n_k = n; 𝐄̃ = [ 𝐄̃_1, ⋯, 𝐄̃_K ] with 𝐄̃_k ∈^m × n_k containing the errors of samples in 𝐙_k; 𝐗̃ = [ 𝐗̃_1, ⋯, 𝐗̃_K ] ∈^m × K d_0 with 𝐗̃_k ∈^m × d_0 the basis to span 𝒮_k; 𝐘̃ = [ 𝐘̃_1; ⋯; 𝐘̃_K ] with 𝐘̃_k ∈^d_0 × n the representation upon the basis 𝐗̃_k; d_0 indicates the subspace dimension, d=K d_0.

In the sequel, we first propose a model to deal with data that are clean, i.e., 𝐄̃=0. Then we generalize our model to handle data that are contaminated with random corruptions or sample-specific outliers, i.e., 𝐄̃≠0.















§ THE PROPOSED METHOD




In this section, we first introduce a theorem that inspires us to formulate the problem with clean data. Then, we generalize the model to handle contaminated data. Finally, we present the algorithm to solve our model.  















 §.§ Problem Formulation

















To start with, we assume that n data samples 𝐙 are strictly generated from K independent subspaces. 
In doing so, we can come to the conclusion as follows:



Suppose K subspaces {𝒮_k}_k=1^K are independent and 𝐗̃_k consists of basis that only span subspace 𝒮_k, then the solution of representation matrix 𝐘̃^⋆ to Eq.(<ref>) dealing with clean data is block-diagonal:


    𝐘̃^⋆ =
    [
         [ Ŷ_1   ⋯   0;   ⋮   ⋱   ⋮;   0   ⋯ Ŷ_K;     ]],


where Ŷ_k ∈^d_0 × n_k and 𝐘̃^k = [0, ⋯, Ŷ_k, 0, ⋯].



See Appendix <ref>.


Theorem <ref> shows the condition under which we can achieve the block-diagonal representation matrix, which is of vital importance to support analyzing the structure of multi-subspace. That is, if we know the basis that clean data lie in, then each block of the representation matrix  characterizes each subspace. 
In practice, however, the basis of each subspace are often unknown, Thus, the generated representation matrix cannot be guaranteed to be block-diagonal. Fortunately, Theorem <ref> offers us the hint to construct a satisfactory model. 
In short, we expect to simultaneously learn the basis matrix and learn the representation matrix inspired by Theorem <ref>. 

To be more specific, we first impose an orthonormal constraint on the basis matrix to ensure that its columns can indeed form the “basis" to span independent subspaces. 

In order to generate the block-diagonal representation, we characterize each data sample as a nonnegative combination[We only deal with data with nonnegative entries.] of only the basis that span its underlying subspace. To this end, we leverage a column l_0-norm constrained on the representation matrix and set its number of nonzero elements equal to the subspace dimension. 



Formally, the objective function satisfying above constraints is defined as



    min_𝐗 ,𝐘𝐙-𝐗𝐘_F^2 
       s.t.  𝐗^𝐓𝐗=𝐈_d,  
           𝐘≥ 0, 
           𝐲_i _0 = d_0, ∀ i,


















where Eq.(<ref>) is the orthonormal constraint on 𝐗; 𝐈_d is a d × d identity matrix. Eq.(<ref>) indicates a nonnegative combination of the orthonormal basis; “≥" is taken component-wise. Eq.(<ref>) is the requirement of the number of nonzero coefficients. For writing simplicity, we omit the “∀ i” symbol in the subsequent related equations.



In practice, the observed data are often contaminated with errors, e.g., noises, random corruptions or sample-specific outliers, rendering the data away from their exact underlying subspaces. To handle different types of errors, we adopt different regularizations. Specifically, we generalize the objective function in Eq.(<ref>) to









    min_𝐗,𝐘,𝐄𝐙-𝐗𝐘-𝐄_F^2 + λ𝐄_Δ
       s.t.  𝐗^𝐓𝐗=𝐈_d,  𝐘≥ 0,  𝐲_i _0 = d_0,








where 𝐄_Δ is some type of matrix norm. For random corruptions, it is 𝐄_1; While for sample-specific outliers, it is 𝐄_2,1. λ>0 is the regularization parameter that controls the effect of two terms.

By implementing above two objective functions, we claim that when the data are clean, then 𝐘 is block-diagonal after reordering its columns and rows, i.e., 𝐘 = 𝐘̃^⋆. When the data are contaminated with random corruptions or sample-specific outliers to some extent, our experimental results show that 𝐘 is still approximate block-diagonal.

It is necessary to notice that SSC and LRR can also be boiled down to the matrix factorization framework. They share some similarities with the proposed method, but there exist two major differences: On one hand, they take advantage of what they refer to as “self-expressiveness" property and use the data themselves to form the dictionary (basis). However, leveraging the data as the basis fail to perform data reconstruction when they contain errors to some degree (See Figure <ref>). On the other hand, they respectively utilize an indirect l_1-norm and nuclear norm to obtain the sparse and low-rank representation of the data, which however, are unable to know the exact basis that are used to represent each data sample. Thus, their ability for subspace clustering is also limited (See Figure <ref>).





 §.§ Problem Optimization


In this part, without loss of generality, we just solve for the general case in Eq.(<ref>). Notice that we need to handle the nonconvex multivariable objective function, as well as the nonsmooth l_0-norm constraint. Fortunately, for function with multivariables, the ADMM method is often the choice <cit.>. Moreover, proximal algorithms has recently been a popular tool for solving nonsmooth and nonconvex problems. It's basic operation is evaluating the proximal operator of a univariable function, which often admits a closed-from solution or can be solved very quickly <cit.>.  Inspired by such superiority, we develop a first-order alternating direction algorithm to solve Eq.(<ref>). The overall procedure of the algorithm is first introducing auxiliary variables and quadratic penalties into the objective function, then iteratively minimizing the augmented Lagrangian function with respect to each primal variable, and finally updating the multipliers. To tackle the subproblem on nonsmooth l_0-norm constraint, we design a novel proximal operator to solve it efficiently and analytically.

We begin with introducing an auxiliary variable 𝐕 to variable 𝐘 which exists in both the object function and the constraint. We reformulate Eq.(<ref>) as

    min_𝐗,𝐘,𝐄,𝐕𝐙-𝐗𝐘-𝐄_F^2 + λ𝐄_Δ, 
           s.t.  𝐗^𝐓𝐗=𝐈_d, 𝐘 = 𝐕, 𝐕≥ 0, 𝐯_i _0 = d_0.


Then, the augmented Lagrangian function of Eq.(<ref>) is


    ℒ_𝒜(𝐗,𝐘,𝐕,𝐄,𝐏) = 𝐙-𝐗𝐘-𝐄_F^2 + λ𝐄_Δ
              + ⟨𝐏, 𝐘-𝐕⟩ + β/2𝐘-𝐕_F^2, 
         s.t.  𝐗^𝐓𝐗=𝐈_d, 𝐕≥ 0,  𝐯_i _0 = d_0,


where 𝐏 is the Lagrangian multiplier, β >0 is the quadratic penalty parameter. Note that adding the penalty term does not change the optimal solution, since any feasible solution satisfying the constraint in Eq.(<ref>) vanishes the penalty term.

Finally, our alternating direction algorithm consists of iteratively minimizing Eq.(<ref>) with respect to one of 𝐗,𝐘,𝐕,𝐄 while fixing the others, and updating the multiplier 𝐏. Specifically, we solve for variables and the multiplier step by step. 

Update X: By discarding terms that are irrelevant to 𝐗, we rewrite the subproblem with respect to 𝐗 as

    min_𝐗𝐙-𝐗𝐘 -𝐄_F^2,   s.t.  𝐗^𝐓𝐗=𝐈_d.

We denote the singular value decomposition (SVD) of (𝐙-𝐄)𝐘^𝐓=𝐋 Σ 𝐑^𝐓, then from Theorem 2 below, we can achieve the closed-form solution of 𝐗. That is,

    𝐗 = 𝐋𝐑^𝐓.




Let 𝐀∈^m × n and 𝐁∈^d × n be any two matrices. Denote the SVD of 𝐀𝐁^𝐓 as 𝐀𝐁^𝐓= 𝐋 Σ 𝐑^𝐓. Then the orthonormal constrained minimization problem

    min_𝐃𝐀 - 𝐃𝐁_F^2     s.t.  𝐃^𝐓𝐃=𝐈_d,

has an analytic solution 𝐃=𝐋𝐑^𝐓.


See Appendix <ref>.


Update Y: The subproblem of 𝐘 is

    min_𝐘𝐙-𝐗𝐘-𝐄_F^2 + ⟨𝐏, 𝐘-𝐕⟩ + β/2𝐘-𝐕_F^2.


Setting the derivative with respect to 𝐘 to zero, rearranging the terms, and using the constraint 𝐗^𝐓𝐗=𝐈_d yields






    𝐘 = (1+β)^-1(𝐗^T (𝐙-𝐄) + β𝐕 - 𝐏).


Update E: The subproblem of 𝐄 becomes

    min_𝐄𝐄-(𝐙-𝐗𝐘)_F^2 + λ𝐄_Δ.

For sample-specific outliers, we set 𝐄_Δ = 𝐄_2,1. Then the solution to Eq.(<ref>) is equivalent to solving the proximal operator for l_2,1-norm. From <cit.>, 𝐄 can be obtained with a closed-form. Concretely, denote 𝐆=𝐙-𝐗𝐘, then

    𝐞_i =
            
                (1-λ/2/𝐠_i _2) 𝐠_i     𝐠_i _2 ≥λ/2, 
    
                0    𝐠_i _2 < λ/2.

Similarly, for random corruptions, we set 𝐄_Δ = 𝐄_1. The solution then equals to finding the proximal operator for l_1-norm, and can be achieved via the Soft-Thresholding Operator <cit.>

    𝐄 =
        𝐆 - λ/2   𝐆≥λ/2, 
    0   -λ/2≥𝐆≤λ/2, 
    𝐆 + λ/2   𝐆≤ -λ/2.


Update V: The subproblem associated with 𝐕 is

    min_𝐕𝐕-(𝐘+β^-1𝐏) _F^2, 
          s.t   𝐯_i ≥ 0, 𝐯_i _0 = d_0.






Here, we emphasize that due to the discrete l_0-norm constraint in Eq.(<ref>), 𝐕 cannot be obtained via gradient or subgradient based methods. Benefiting from aforementioned proximal algorithms <cit.>, we can design a simple yet very efficient algorithm to solve 𝐕 column by column. The detail is as follows:

We denote 𝐔 = (𝐘+ β^-1𝐏) and define the operator P_d_0:^d →^d for 𝐮_i as

    P_d_0(𝐮_i) = _𝐯_i{ || 𝐮_i - 𝐯_i ||_2^2 : ||𝐯_i||_0 = d_0 }.


We also define the nonnegative orthant mapping as

    P_+(𝐮_i) = _𝐯_i{ || 𝐮_i - 𝐯_i ||_2^2 : 𝐯_i ≥ 0 } = max{0, 𝐮_i}.

and introduce an indicator function I_𝒱(𝐯) over the set 

    𝒱 = {𝐯∈^d: 𝐯≥ 0, ||𝐯||_0 = d_0}.







Then, we have the following theorem: 


Given above defined operator P_d_0, nonnegative orthant mapping P_+, and indicator function I_𝒱(𝐯). Then the solution 𝐯_i of Eq.(<ref>) equals to solving the proximal operator[Note that there is a factor 1/2 on the quadratic term in the standard form. Here we remove it to simplify the notation.] of the indicator function 

    prox_I_𝒱(𝐮_i)     = _𝐯_i{ || 𝐮_i - 𝐯_i ||_2^2: 𝐯_i ≥ 0, ||𝐯_i||_0 = d_0}
        = _𝐯_i ∈𝒱𝐮_i - 𝐯_i ^2,

which has an analytical form, i.e., 𝐯_i = P_d_0(P_+(𝐮_i)).



See Appendix <ref>. 
 

Based on Theorem 3, 𝐯_i can be acquired by selecting d_0 largest nonnegative entries of 𝐮_i with corresponding indices
𝐪_j=[ q_1,j, ⋯, q_d_0,j]. To be specific, for i=1,⋯,d, we set

    v_i,j =
    
    u_i,j   i ∈𝐪_j, 
    
    0    i ∉𝐪_j.


Update P: Finally, for the multiplier 𝐏, based on the dual optimal condition <cit.>, its updating rule is

    𝐏 := 𝐏 + μ (𝐘-𝐕),

where μ = min(ρμ, μ_max), with pregiven ρ and μ_max.

Algorithm 1 summarizes the algorithmic procedure of our proposed alternating direction algorithm for solving Eq.(<ref>). 

After obtaining 𝐗, 𝐘, and 𝐄, MFC_0 can accomplish the following tasks: 


  * Subspace Clustering: We apply normalized cut <cit.>[One can also use K-means clustering on 𝐘. We choose normalized cut here so as to have a fair comparison with SSC and LRR.] to 𝐘^T𝐘 to cluster all data samples to their respective subspaces. 


  * Data Reconstruction: We treat the product of 𝐗 and 𝐘 as the reconstruction for the data.

  * Error Correction: We regard 𝐄 as the error for the data. 

  * Representation Learning: We permute the rows and/or columns of 𝐘 to make it be (appximate) block-diagonal. Then the block 𝐘̃_k indicates the representation for subspace 𝒮_k.   

  * Basis Learning: We respectively extract d_0 columns of 𝐗, whose indices correspond to the row numbers of 𝐘̃_k (k=1,⋯, K), to form the basis of subspace 𝒮_k (k=1, ⋯, K).  










§ COMPUTATIONAL ANALYSIS








The computational complexity of MFC_0 in each iteration can be summarized as follows:

Calculate 𝐗 
We first perform matrix subtraction of size m × n and matrix multiplication of size m × n, n × d with complexities O(mn) and O(dmn), respectively. Then, we compute the SVD of an m × d matrix and store both the singular vectors—the complexity is O(m^2d+d^3). Finally, we calculate the matrix product of size m × d and d × d with complexity O(d^2m) to obtain the orthonormal basis 𝐗. The total complexity is O(m^2d+d^3)+O(d^2m+dmn).

Compute 𝐘 
We just perform matrix multiplication of size n × m and m × d and matrix addition of size m × d, with the dominant complexity O(dmn).

Compute 𝐄 
The computation consists of matrix multiplication, matrix substraction, and proximal operator calculation[Both 𝐄_2,1 and 𝐄_1 have the same complexity O(mn).]. The complexity of each part is O(dmn), O(mn), and O(mn), and the dominant complexity is O(dmn).

Compute 𝐕 
The complexity of matrix addition is O(dn). The nonnegative orthant mapping and sorting of d-dim vector respectively need O(d) and O(d ln d) (in the average case) calculation amount. Thus, the dominant complexity is O(nd ln d).

As the dimensionality of multi-subspace d is much smaller than the data dimensionality m and sample size n, i.e., m,n >> d, we conclude that the dominant complexity of MFC_0 is O(dm^2+dmn).





Here, we also discuss the computational complexity of SSC and LRR. For SSC, the computational burden is matrix multiplication, with the complexity O(n^2m). For LRR, the computational burden lies in two parts: matrix multiplication and matrix SVD, the complexity of which are O(n^2m) and O(n^3), respectively. One should note that: (1) the matrix size for matrix multiplication of MFC_0  (m × n and n × d ) is smaller than that of SSC and LRR (m × n and n × n); (2) the matrix size for SVD of MFC_0 (m × d) is smaller that of LRR (m × n). The dominate complexity of SSC, LRR, and MFC_0 are listed in Table <ref>.

Compared MFC_0 with SSC and LRR, we observe that when handling data with high dimensionality, SSC and LRR are relatively faster—the complexity in terms of m is quadratic for MFC_0, and linear for SSC and LRR. Actually, we can leverage random projections, a very 
effective <cit.> feature extraction method, as a preprocessing to reduce the dimensionality, thus largely decreasing the complexity of MFC_0. In contrast, if we deal with large data sample size, MFC_0 is more efficient than SSC and LRR—the complexity of SSC and LRR are quadratic and cubic to the sample size n, while MFC_0 is linear.

In practice, the overall computational complexity is determined by not only the complexity of each iteration but also the number of total iterations. Although there is no theoretical guarantees, experimental results (See Tables <ref>-<ref>) show that MFC_0 can converge with 20-40 iterations, while LRR and SSC need 55-70 and 90-110 iterations, respectively. 








§ EXPERIMENTAL RESULTS


In this section, we carry out several experiments on both synthetic data and real-world datasets to test the performance of our proposed MFC_0 method. 


Compared methods We compare MFC_0 with single subspace learning methods, i.e., PCA, NMF <cit.>, and Sparse NMF (SNMF) <cit.>; and state-of-the-art multi-subspace learning methods, i.e., SSC and LRR. MFC_0 can leverage different regularizations to handle different types of errors. In the experiment, we mainly focus on random corruptions (||𝐄||_1) and sample-specific outliers (||𝐄||_2,1). Accordingly, we denote MFC_0 as MFC_01 and MFC_02, respectively.

For PCA, we preserve 95 percent of total data variance; For NMF and SNMF, we set their basis number as K d_0, the same as MFC_0's. For MFC_0, we randomly generate a matrix with size m × K d_0 to initialize the basis matrix. The hyperparameters of all comparable methods are selected via cross-validation. The stopping criterion for each method is either it reaches the maximal iterations 10^3 or the objective function values between neighboring iterations have difference less than 10^-4. All experiments are conducted in Matlab R2010b with the platform CPU 3.10 GHz and RAM 16.0 GB.

Evaluation metrics 
We use clustering accuracy, time, and iterations as the metrics to evaluate the performance of all compared methods.



  §.§.§ Clustering Accuracy
 
The subspace clustering result is evaluated by comparing the estimated label of each data sample with that provided by the ground truth. In the paper, we use clustering accuracy (ACC) to measure the clustering performance. ACC is defined as

    ACC = 1/n∑_i=1^n δ(s_i,map(r_i)),

where r_i and s_i are the estimated label and the ground truth of the i-th point; δ(x,y)=1 if x=y, and δ(x,y)=0 otherwise. map(x) is the permutation mapping function that maps each label r_i to the equivalent label from the entire data. The best mapping can be efficiently computed by the Kuhn-Munkres algorithm <cit.>.



  §.§.§ Time and Iterations

We use three indices to quantize the computational effciency, i.e., total running time T(s), iterations I, and averaged running time per iteration (T/I). T indicates the overall effciency; I reflects the convergence rate to reach local solutions, and T/I shows the efficiency per iteration. 



 §.§ Results on Synthetic Data



  §.§.§ Data reconstruction visualization on 3-dim space


We visualize data reconstruction in 3-dim ambient space in Figure <ref> and compare MFC_0 with SSC <cit.> and LRR <cit.>. The data samples are drawn from three 1-dim independent subspaces, and disturbed by random corruptions and sample-specific outliers, respectively. 

Figure <ref> shows that SSC and LRR are unable to remove the errors, which validates our aforementioned argument; While MFC_0 fully eliminates them. The reason is that both SSC and LRR use the original contaminated data as the basis, which is unreasonable for data reconstruction. In contrast, MFC_0 learns the basis that span the underlying subspace the clean data lie in–Indeed, the learnt basis are [0, √(2)/2, √(2)/2], [√(2)/2, 0, √(2)/2], and [√(2)/2, √(2)/2, 0], respectively; Errors contained in the data are absorbed in the regularization term.



  §.§.§ Results on high-dimensional data

For high-dimensional data analysis, the synthetic data are generated from K=5 independent subspaces with each containing n_k=100 samples. All subspaces have the same dimension d_0=10 embedded in a D=100 dimensional ambient space. The procedure of generating above subspaces is similar to that of <cit.>: the basis 𝐔_k of each subspace are calculated by 𝐔_k+1 = 𝐓𝐔_k, 1 ≤ k ≤ 4, where 𝐓=orth(rand(D))∈^D × D is a random orthonormal matrix and 𝐔_1 ∈^D × d_0 is a random column orthogonal matrix. The data points from each subspace are sampled by 𝐙_k=𝐔_k 𝐑_k, where elements in 𝐑_k ∈^d_0 × n_k are from standard uniform distribution. 


Representation learning and stability We exhibit MFC_0's ability to learn the representation matrix, i.e., the block-diagonal structure, of data without and with errors, and show its efficiency and stability to obtain the local solution. 







The representation matrix of clean data is shown in Figure <ref>(a) and objective function values versus iteration in 10 runs are plotted in Figure <ref>(b). From Figure <ref>, we can see that MFC_0 can precisely discover the block-diagonal structure when the data samples are clean. Moreover, our proposed first-order optimization algorithm can fast converge (nearly 15 iterations) to a very stable value although using different random initializations. These observations demonstrate that MFC_0 is both effective and efficient to analyze the structure of multiple subspaces.

We further consider the cases that data contain errors. For random corruptions and sample-specific outliers with error ratio[The proportion of data samples are contaminated with errors.] = 0.1, 0.3, and 0.5, the corresponding representation matrices are plotted in Figure <ref> and Figure <ref>, respectively. We observe that even the data are contaminated with different types of errors to a high level, the representation matrix learnt via MFC_0 is still close to be block-diagonal.

Clustering accuracy vs. error ratio We compare MFC_0 with all the other methods in terms of subspace clustering to validate its robustness to resist different types of errors. The performance of all compared methods versus different error ratios, ranging from 0 to 0.8, of random corruptions and sample-specific outliers in the data are depicted in Figure <ref>(a) and <ref>(b). It can be seen that the clustering accuracy of MFC_0 is consistently the highest. SSC and LRR come next. Surprisingly, even when the ratio is 0.6, MFC_0 can still obtain >95% accuracy in both cases. Whereas, for LRR and SSC, the accuracies are approximate 70%, 85% for random corruptions, and 60%, 60% for sample-specific outliers, respectively. The performance of PCA, NMF, and SNMF decrease sharply, which validates that single subspace learning methods fail to work when the error ratio increases.  



  §.§.§ Summary
 


  * MFC_0 can accurately reconstruct the data and discover the basis that span each subspace;   


  * MFC_0 is capable of obtaining the exact block-diagonal representation matrix when data is clean; 


  * MFC_0 shows strong ability to resist errors with a high ratio and is more powerful than state-of-the-art LRR and SSC methods for subspace clustering, not to mention single subspace learning methods; 


  * MFC_0 can fast converge to a very stable solution via our proposed first-order optimization algorithm. 












 §.§ Results on Real-World Datasets


In this part, we pay close attention to the tasks: subspace clustering, as well as associated computational complexity; subspace basis learning and subspace recovery. Experiments are carried out on three real-world face datasets: AR[<http://www2.ece.ohio-state.edu/ aleix/ARdatabase.html>], Yale[<http://vision.ucsd.edu/content/yale-face-database>], and Extended Yale B (EYaleB)[<http://vision.ucsd.edu/ leekc/ExtYaleDatabase/ExtYaleB.html>] and the USPS handwritten digit dataset[<http://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/multiclass.html#usps>]. 



  §.§.§ Dataset Description


AR contains over 3,000 images of 126 individuals (70 men and 56 women). They are taken at two different occasions with different facial expressions, illumination conditions, and occlusions (sun glasses and scarf) under strictly controlled conditions. Raw images are of 768x576 pixels and 24 bits of depth. In the experiment, images are resized to 48x48 pixels. 

Yale contains 165 grayscale images of 15 individuals. There are 11 images per subject, one per different facial expression or configuration: center-light, w/glasses, happy, left-light, w/no glasses, normal, right-light, sad, sleepy, surprised, and wink. We simply use the cropped images with size 32x32 pixels. 

EYaleB consists of 2414 images of 38 subjects under 9 poses and 64 illumination conditions. The subset used in the experiment has 20 individuals and around 64 near frontal images under different illuminations per individual. All images are cropped to 32x32 pixels.

USPS contains 9198 handwritten digit images (from 0 to 9) and the image size is 16x16. In the experiment, we choose 100 images of each digit, thus 1000 images in total. 

Some sample images from the four datasets are shown in Figure <ref>. Statistics of the selected subsets are listed in Table <ref>.







  §.§.§ Subspace Dimension Estimation


For MFC_0, we need to know the subspace dimension of each subject in advance. A simple yet effective way is to calculate the number of nonzero singular values of each subject. Figure <ref> depicts the singular values of several randomly chosen subjects in the four datasets. The subspace dimensions d_0 for AR, Yale, EYaleB, and USPS are approximate 12, 10, 10, and 12 respectively. Note that for real-world datasets, the singular values often gradually close to zero, validating that the images are contaminated with errors. The higher dimension of AR reflects its more complicated face distribution than Yale and EYaleB. 
It is necessary to point out that, we cannot distinguish which type of error contained in the images, therefore we report  






  §.§.§ Subspace Clustering


In this part, we perform subspace clustering on real-world datasets: Given face images or handwritten digits of multiple categories, our task is to group them according to their respective subspace. Our purpose is to verify that it is more advantageous to (1) model face distribution using multiple subspaces than using the single one; (2) learn the basis rather than leverage the data themseleves as the basis. 




Tables <ref> and <ref> show the clustering performance on Yale and EYaleB, respectively. To study the effect of category number, we choose K that is ranging from 2 to 11 on Yale and from 2 to 20 with an interval 2 on EYaleB.

From the results, we have the following observations: 

(1) As a whole, MFC_02 performs the best on both datasets. By learning each subject the respective subspace, together with a direct column l_0-norm constraint and an error correction term, MFC_02 is discriminative to separate face images from the subjects. Essentially, this is the key factor that makes it outperform single subspace analysis methods. Moreover, using a column l_0-norm constraint to directly control the sparsity, MFC_02 is superior to the indirect l_1 penalty and avoids tuning the sparse hyperparameter; 

(2) SSC and LRR can also achieve promising results as they are multi-subspace learning methods with error correction. However, that they use contaminated face images as the basis is problematic and irrational, the accuracy gap between them and MFC_02 validates this argument; 

(3) Errors contained in face images are more likely to be sample-specific outliers, as on both Yale and EYaleB, the result of MFC_02 is better than that of MFC_01. A corroborative evidence is that LRR, also using the 𝐄_2,1 term, is slightly superior to SSC, the 𝐄_1 term; 

(4) Sparse constraint for clustering is important. The better results of SNMF than that of NMF in both datasets verify the importance using sparse constraint;

(5) Modeling face images distribution using multi-subspace Since PCA, NMF, and SNMF are single subspace learning methods, so their performance are worse than SSC, LRR, MFC_01, MFC_02. When the sample size of each subject is small in Yale, the structure of multi-subspace is not so obvious and the accuracy gap is not remarkable. However, when the sample size of each subject in EYaleB is sufficient, all of them obtain poor results on the whole. Even worse, they fail to work when K increases to a certain value, say 10. While for remaining methods, they still work well.

We also implement subspace clustering on USPS with K from 2 to 10. Our aim is to see whether handwritten digits are distributed in “independent subspaces", as no pervious work mentions such assumption. 

Subspaces of handwritten digits are probably not independent. From the results listed in Table <ref>, we notice that all PCA, NMF, and SNMF have promising results, even comparable with SSC and LRR in some Ks. Still, our MFC_01 performs the best, but the superiority is not so obvious as that in EYaleB. One possible reason is that the subspaces of handwritten digits are not independent in essence. In this sense, the distribution of handwritten digits modeled by disjoint subspaces (or submanifolds) may be more reasonable.









  §.§.§ Computational Complexity


Another advantage of MFC_0 against LRR and SSC is its efficient computation in terms of sample size. Therefore, our aim in this part is to validate that when increasing the sample size MFC_0 is more efficient than SSC and LRR. As mentioned, we quantize the efficiency of each method using total running time (T(s)), iterations (I), and averaged running time per iteration (T/I). 

The results on Yale, EYaleB, and USPS are listed in Tables <ref>, <ref>, <ref>, respectively. From Table <ref>, we see that 

(1) Running time per iteration T/I of MFC_01 or MFC_02 is bigger than that of SSC and LRR on dataset with small sample size. The averaged T/I of MFC_01 or MFC_02 is approximately 5 and 3 times of SSC and LRR. It is because when the sample size of each subject n_k=15 is small in Yale, the dimensionality of data samples (m=32 × 32) dominates the computational complexity in each iteration. Remembering that the complexity of MFC_0 is quadratic with respect to m, while of LRR and SSC is linear. 

(2) Iterations I of MFC_0 is much smaller than those of SSC and LRR. Despite of higher T/I, the averaged iteration I of MFC_01 or MFC_02 is 24 or 23, much smaller than those of SSC and LRR, which are 93 and 57, respectively. Such observation demonstrats MFC_0's faster convergence rate than LRR and SSC. Ultimately, the total running time T of these methods are similar.

(3) Increasing sample size makes MFC_0 more advantageous. When the sample size of each subject increases to 64 in EYaleB, the results in Table <ref> manifest that T/I values of SSC and LRR are larger than that of MFC_01 or MFC_02 on average. It is due to the fact that the computational complexity of MFC_0 is linear with respect to n, while quadratic and cubic of LRR and SSC. Similar to the results on Yale, compared with SSC (110) and LRR (65), the averaged values I of MFC_01 (31) and MFC_02 (30) are much smaller. Consequently, the total running time T of MFC_01 and MFC_02 are just 1/4, 1/5.7 and 1/6, 1/8.5 of SSC and LRR. To further highlighting the efficiency of MFC0_0, we add samples to 100 in USPS and related results are displayed in Table <ref>. We observe that when the category number is 10, it takes SSC and LRR almost 94 and 255 seconds to implement the clustering task, while both MFC_01 and MFC_02 consume less than 10 seconds. Again, the less I of MFC_01 or MFC_02 than SSC and LRR suggests they can reach the stable solution more quickly.










































  §.§.§ Face Basis Learning and Reconstruction on AR


Subspace basis learning aims to learn the representative basis that provide a compact representation for each data sample. Traditional methods assume that the face images of all subjects can be represented by a single lower-dimensional subspace. However, previous work <cit.> revealed that, under the Lambertian assumption, face images of a subject with a fixed pose and varying illumination lie close to a linear subspace. Thus for multiple subjects, the reasonable way is to learn each subject its underlying basis. To confirm this argument, we conduct an experiment to visualize the learnt face basis (Note that SSC and LRR are unable to learn the basis, thus we do not show their results). Moreover, we also display the reconstructed face image and the associated error. To see the difference with single subspace learning methods, we choose PCA for comparison.
















We randomly choose three subjects (subject 1, 2, and 55, each 15 samples) in AR to learn the basis, which are then used for face reconstruction. Since some face images in AR are grossly corrupted, for instance, wearing sunglass, we therefore treat them as sample-specific outliers and choose MFC_02. Some examples of learnt face basis, reconstructed face images, and errors of PCA and MFC_02 are plotted in Figure <ref>. From it, we have several observations: 

(1) MFC_0 can learn each subpace the corresponding basis. The learnt basis via MFC_02 can be separated into three parts—basis in each row represent one subject, meaning face images of one subject lie on its own underlying subspace. In contrast, the basis of PCA are mixed together, i.e., each basis contains information of face images across the three subjects;

(2) MFC_0 is effective in reconstructing face images. The reconstructioned face images via MFC_02 are very clean, while they are blurry and untidy via PCA. 

(3) MFC_0 is able to correct discriminant errors, which can further be used to perform recognition tasks. The errors via MFC_02 are almost extracted. For Figure <ref>(c) and (d), the errors are  “teeth"—One reasonable explanation is that face images in the training set do not show their teeth. Similarly, the errors in Figure <ref>(e) and (f) are “sunglass". Interestingly, the errors contain discriminant features. In fact, they can be used for facial emotion recognition or object recognition <cit.>. For instance, we can leverage these errors to detect whether one is laughing or wearing sunglass. Nevertheless, the errors obtained by PCA cannot provide any useful information.



  §.§.§ Summary
 


  * MFC_0 (either MFC_01 or MFC_02) is not only effective for subspace clustering, subspace basis learning, but also implemented efficiently.


  * The advantage of MFC_0 against single subspace learning methods can be magnified when handling data generated from multiple subspaces. 


  * The advantage against SSC and LRR can be boosted by handling data contained gross errors or/and data with large sample size.





§ CONCLUSIONS AND FUTURE WORK


In this paper, we propose a novel method, called Column L_0-norm constrained Matrix Factorization (MFC_0), for robustly analyzing the structure of data generated from multiple categories. By learning the basis with an orthonormal constraint, MFC_0 is able to discover the mixture subspace structure and is robust to different types of errors with the specific regularization. Moreover, MFC_0 directly imposes an l_0-norm constraint on the representation matrix, which achieves a (or approximate) block-diagonal structure when the data are clean (or contain errors). We propose a very efficient first-order alternating direction type algorithm to stably solve the nonconvex and nonsmooth objective function of MFC_0. Experimental results on synthetic data and real-world datasets verify that, besides the superiority and efficiency over traditional and state-of-the-art methods for multi-subspace recovery and clustering, MFC_0 also demonstrates its uniqueness for multi-subspace basis learning and direct representation learning.

Although MFC_0 owns these advantages, it also faces several problems that should be investigated in future. First, we assume that the multiple subspaces are independent, however in practice the data may be generated from disjoint subspaces. In this case, how can we incorporate the constraint or regularization into the model? More challenging, how can we process the data that are generated from nonlinear manifolds, instead of linear subspaces? Second, the computational complexity of MFC_0 is linear with respect to the sample size, so its extension to deal with massive dataset will be investigated in future.
























































































































































IEEEtran




theoremsection



 §.§ Proof of Theorem <ref>

 
Let 𝐘^⋆ be an optimal solution of Eq.(<ref>) under the given basis 𝐗. We can then decompose 𝐘^⋆ into two parts, i.e., the diagonal 𝐘^D and non-diagonal 𝐘^N, where

    𝐘^D =
    [
         [ Ŷ_1^⋆     0     ⋯     0;     0 Ŷ_2^⋆     ⋯     0;     ⋮     ⋱     ⋮     ⋮;     0     0     ⋯ Ŷ_K^⋆;       ]],
    𝐘^N =
    [
         [ 0 ⋆ ⋯ ⋆; ⋆ 0 ⋯ ⋆; ⋮ ⋱ ⋮ ⋆; ⋆ ⋆ ⋯ 0;   ]]

with Ŷ_k^⋆∈^d_0 × n_k and ⟨𝐘^D, 𝐘^N⟩ = 0.

Let 𝐳^⋆_t, 𝐳^D_t, 𝐳^N_t respectively denote the t-th column of 𝐗𝐘^⋆, 𝐗𝐘^D, 𝐗𝐘^N. Suppose that 𝐳^⋆_t ∈𝒮_k. It is easy to check 𝐳^D_t ∈𝒮_k, therefore 𝐳^N_t ∈∑_j ≠ k𝒮_j. However, 𝐳^N_t = 𝐳^⋆_t - 𝐳^D_t ∈𝒮_k. As all subspaces are independent, meaning 𝒮_k ∩∑_j ≠ k𝒮_j = { 0 }, we must require 𝐳^N_t = 0. Thus 𝐗𝐘^N = 0, and 𝐙 = 𝐗𝐘^D, 𝐘^⋆=𝐘^D. Therefore, the solution 𝐘^⋆ to Eq.(<ref>) is block-diagonal.






 §.§ Proof of Theorem <ref>

 
We first unfold the objective function:

    min𝐀-𝐃𝐁_F^2 = tr(𝐀^𝐓𝐀) - 2tr(𝐀𝐁^𝐓𝐃^𝐓) + tr(𝐁^𝐓𝐃^𝐓𝐃𝐁).


Using the constraint 𝐃^𝐓𝐃=𝐈 and the SVD of 𝐀𝐁^𝐓= 𝐋 Σ 𝐑^𝐓 transforms above minimization to the maximization problem

    maxtr(𝐀𝐁^𝐓𝐃^𝐓) = tr(𝐋 Σ 𝐑^𝐓 𝐃^𝐓) = tr(Σ (𝐃𝐑)^𝐓 𝐋).

Note that (𝐃𝐑)^𝐓 𝐋_F^2 = tr((𝐃𝐑)^𝐓 𝐋 𝐋^𝐓 (𝐃𝐑)) = tr(𝐈_d) is a constant. As Σ is a diagonal matrix, so the maximum is achieved when (𝐃𝐑)^𝐓 𝐋 is also a diagonal matrix, with positive diagonal elements. Therefore, we have 𝐃𝐑 = 𝐋, and thus 𝐃=𝐋𝐑^𝐓.




 §.§ Proof of Theorem <ref>


For any 𝐮_i, we first split its elements into the negative and positive parts, with corresponding indices defined as
ℐ_+ = { j | u_i,j≥ 0} and ℐ_- = { j | u_i,j < 0 }. Denote 𝐮_i_+^2 = ∑_j ∈ℐ_+ u_i,j^2 and 𝐮_i_-^2 = ∑_j ∈ℐ_- u_i,j^2. Then, we have two useful properties: 

(1) 𝐮_i _2^2 = 𝐮_i _+^2 + 𝐮_i _-^2; 

(2) 𝐯_i - 𝐮_i _+^2 + 𝐯_i_-^2 = 𝐯_i- P_+(𝐮_i) _2^2.

Property (1) is obvious and (2) can be simply proved: 

    𝐯_i- P_+(𝐮_i) _2^2 = 𝐯_i- P_+(𝐮_i) _+^2 + 𝐯_i- P_+(𝐮_i) _-^2 
          = 𝐯_i- 𝐮_i _+^2 + 𝐯_i- 0 _-^2 = 𝐯_i - 𝐮_i _+^2 + 𝐯_i_-^2.


Based on above properties, we have

    prox_I_𝒱(𝐮_i)    = _𝐯_i{ || 𝐮_i - 𝐯_i ||_2^2: 𝐯_i ≥ 0, ||𝐯_i||_0 = d_0}
        = _𝐯_i ∈𝒱{ || 𝐮_i - 𝐯_i ||_2^2 }
        = _𝐯_i ∈𝒱{ || 𝐮_i - 𝐯_i ||_+^2 + || 𝐮_i - 𝐯_i ||_-^2 }
       ⇔_𝐯_i ∈𝒱{ || 𝐮_i - 𝐯_i ||_+^2 + || 𝐯_i ||_-^2 - 2∑_j ∈ℐ_- v_i,j u_i,j}
        = _𝐯_i ∈𝒱{ || 𝐯_i - P_+(𝐮_i) ||_2^2 - 2∑_j ∈ℐ_- v_i,j u_i,j}

Note that we remove the constant || 𝐮_i||_-^2 in the fourth equation. 

In the last equation, when j ∈ℐ_-, u_i,j <0. As v_i,j≥ 0, so -2∑_j ∈ℐ_- v_i,j u_i,j≥ 0. Therefore, the minimum is achieved if and only if v_i,j = 0 for j ∈ℐ_-, resulting in ||𝐯_i||_-^2 = 0. Finally, we have 

    prox_I_𝒱(𝐮_i)
            = _𝐯_i{ || 𝐯_i - P_+(𝐮_i) ||_2^2: ||𝐯_i||_0 = d_0 }
        = P_d_0(P_+(𝐮_i)).

























































