Sparse Nonnegative CANDECOMP/PARAFAC Decomposition in Block Coordinate Descent Framework: A Comparison Study
Deqing Wang Fengyu Cong Tapani Ristaniemi











































Journal of  Class Files, Vol. XX, No. X, Month 2018
Shell et al.: Bare Demo of IEEEtran.cls for IEEE Journals


























Sparse Nonnegative CANDECOMP/PARAFAC Decomposition in Block Coordinate Descent Framework: A Comparison Study
    Deqing Wang, Student Member, IEEE,
        Fengyu Cong, Senior Member, IEEE,
        and Tapani Ristaniemi, Senior Member, IEEE
This work was supported in part by the National Natural Science Foundation of China (Grant No. 81471742), in part by the Fundamental Research Funds for the Central Universities [DUT16JJ(G)03] in Dalian University of Technology in China, and in part by the scholarship from China Scholarship Council (No. 201600090043). (Corresponding author: Fengyu Cong)
D. Wang and F. Cong are with the School of Biomedical Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian 116024, China, and also with the Faculty of Information Technology, University of Jyväskylä, Jyväskylä 40100, Finland (e-mail: deqing.wang@foxmail.com; cong@dlut.edu.cn).
T. Ristaniemi is with the Faculty of Information Technology, University of Jyväskylä, Jyväskylä 40100, Finland (e-mail: tapani.e.ristaniemi@jyu.fi).
Manuscript received Month Day, Year; revised Month Day, Year.
    December 30, 2023
=============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






Nonnegative CANDECOMP/PARAFAC (NCP) decomposition is an important tool to process nonnegative tensor. Sometimes, additional sparse regularization is needed to extract meaningful nonnegative and sparse components. Thus, an optimization method for NCP that can impose sparsity efficiently is required. In this paper, we construct NCP with sparse regularization (sparse NCP) by l_1-norm. Several popular optimization methods in block coordinate descent framework are employed to solve the sparse NCP, all of which are deeply analyzed with mathematical solutions. We compare these methods by experiments on synthetic and real tensor data, both of which contain third-order and fourth-order cases. After comparison, the methods that have fast computation and high effectiveness to impose sparsity will be concluded. In addition, we proposed an accelerated method to compute the objective function and relative error of sparse NCP, which has significantly improved the computation of tensor decomposition especially for higher-order tensor.




Tensor decomposition, CANDECOMP/PARAFAC (CP) decomposition, nonnegative constraint, sparse regularization, block coordinate descent.





















§ INTRODUCTION
















Nonnegative tensor decomposition and nonnegative matrix factorization are powerful tools in signal processing and machine learning <cit.>. Due to the nonsubtractive property and part-based representation <cit.>, they have been widely applied to hyperspectral unmixing <cit.>, cognitive neuroscience <cit.>, chemometrics <cit.>, and many other areas <cit.>. Nonnegative CANDECOMP/PARAFAC (NCP) tensor decomposition is an extension of nonnegative matrix factorization (NMF) from two-way to multi-way. NMF and NCP are data factorization/decomposition methods with nonnegative constraints, which are described as follows.

The NMF problem. Given a nonnegative matrix V∈ℝ^I_W× I_H and a positive integer R<min{ I_W, I_H}, NMF is to solve the following minimization problem:

    min_W,H1/2   V-WH_F^2
    s.t. W   ⩾ 0, H⩾ 0,

where matrices W∈ℝ^I_W× R and H∈ℝ^R × I_H are two estimated nonnegative factors.

The NCP problem. Given a nonnegative Nth-order tensor 𝒳∈ℝ^I_1× I_2×…× I_N and a positive number R, NCP is to solve the following minimization problem:

    min_A^(1),…,A^(N)1/2   𝒳-A^(1),…,A^(N)_F^2
    s.t. A^(n)   ⩾ 0  for  n=1, …,N,

where A^(n)∈ℝ^I_n× R for n=1,…,N are estimated factors on different modes, I_n is the size on mode-n, and R can be seen as the selected rank-1 tensor number (initial components number). The estimated factors in Kruskal operator can be represented by sum of R rank-1 tensors in outer product form:

    A^(1),…,A^(N) = ∑_r=1^R𝒳_r = ∑_r=1^Ra_r^(1)∘…∘a_r^(N),

where a_r^(n) represents the rth column of A^(n).

Both NMF and NCP are non-convex and non-linear optimization problems, therefore finding their golobal minimums is NP-hard <cit.>. Conventionally, problems (<ref>) and (<ref>) can be solved in block coordinate descent (BCD) framework <cit.>, in which each factor is updated alternatively as a subproblem with other factors fixed. These subproblems are usually convex. In NMF case, a lot of optimization methods have been proposed to solve these subproblems. Lee et al. proposed the multiplicative update (MU) method <cit.>, which is the most popular and widely applied method for NMF. Berry et al. introduced the alternating least square (ALS) method in the review paper <cit.>. Cichocki et al. proposed the hierarchical alternating least squares (HALS) method for large-scale NMF problems <cit.>. Xu et al. proposed the alternating proximal gradient (APG) method for NMF with detailed mathematical convergence proofs <cit.>. The same idea as APG was also proposed in <cit.> for NMF independently, which was called NeNMF. Recently, the alternating direction method of multipliers (ADMM) has been employed for NMF problems <cit.>. In addition, the alternating nonnegative least squares (ANLS) method was deeply analyzed in Lin's seminal paper with strong optimization properties <cit.>, which has a significant influence on NMF. When ANLS method is utilized, the subproblems of NMF appears as the nonnegative least squares (NNLS) problems. Many efficient methods have been devoted to solve the NNLS subproblems, such as Lin's project gradient method <cit.>, quasi-Newton method <cit.>, active-set method <cit.>, block principal pivoting method <cit.>, inertial projection neural network <cit.>, and proximal function based method <cit.>. Most of above methods for NMF can be naturally extended to NCP problems <cit.>.

If no other constraint or regularization is imposed, (<ref>) and (<ref>) with only the nonnegative constraints can be seen as bound-constrained optimization problems <cit.>. The nonnegative constraints can guarantee physically meaningful results. Sometimes, incorporating specific regularization can yield more interpretable and accurate results, e.g. sparse regularization, orthogonal regularization, smooth regularization <cit.>, manifold regularization <cit.>, and so on. Due to the intrinsic sparsity in many types of data, such as face image data <cit.>, microarray data <cit.>, text data <cit.>, hyperspectral image data <cit.> and EEG data <cit.>, sparse results are always required in NMF and NCP problems. Nonnegativity constraint will naturally lead to sparsity in the results, but this sparsity in only a side effect, which can't be controlled to a certain level <cit.>. Therefore, explicit sparse regularization is needed. In NMF case, many methods has been proposed to impose sparsity by projection or regularization iterms. Hoyer proposed to project all components in NMF factor into vectors that has desired sparsity degree <cit.>, which can be solved by MU or projected gradient methods <cit.>. However, this method makes all components have the same fixed sparsity degree, which doesn't reveal the true sparsity distribution in data. On the other hand, l_1-norm is a conventional and effective sparse regularizer for signal processing <cit.>. The reason is that for most underdetermined linear equations the minimization problem with l_1-norm regularization can yield strong sparsity <cit.>. Consequently, l_1-norm is a very promising regularization item for NMF and NCP to impose sparsity.

A multitude of works have been devoted to incorporate sparse regularization to NMF <cit.>, but rare works can be found for tensor decomposition. As far as we know, only a few studies had focused on imposing sparsity by l_1-norm regularization to Tucker decomposition <cit.>. In this study, we investigate nonnegative CANDECOMP/PARAFAC (CP) decomposition with sparse regularization, which is abbreviated to “sparse NCP" for convenience. We design the mathematical model of sparse NCP using l_1-norm as the engine to impose sparsity explicitly. In order to prevent rank deficiency and increase the stability, squared Frobenius norm is also added as an auxiliary regularization. The popular optimization methods of multiplicative update (MU), alternating least squares (ALS), hierarchical altering least squares (HALS), alternating proximal gradient (APG) and alternating nonnegative least squares (ANLS) are employed to solve the sparse NCP model, which are in block coordinate descent framework <cit.>. The sparse NCP implemented by all the above optimization methods are tested on synthetic tensor data and real tensor data, both of which contain third-order and fourth-order case. The abilities of these methods to impose sparsity are carefully compared. We want to mention HALS is a special method that includes normalization constraints on factor matrices in addition to nonnegative constraint. With additional normalization constraints, NMF and NCP no longer remain bounded problems <cit.>. Despite the drawback of HALS, we still test its performance on the tensor data for imposing sparsity due to its popularity in the past <cit.>.

Data in tensor, especially higher-order tensor (order ⩾ 4), often consist of a huge amount of points, for which the decomposition might process slowly in some limited conditions. The iterating of tensor decomposition is usually terminated by checking the change of objective function value or data fitting, which is sometimes time consuming. A special strategy based on trace of matrix was used to accelerate the calculation of the objective function at each iteration <cit.>. We extend this strategy to our sparse NCP problem, in which the l_1-norm sparse regularization and Frobenius regularization items can be easily handled.

The rest of this paper is organized as follows. In Section II, we describe the mathematical model of the designed sparse NCP. Section III elucidates the solutions to sparse NCP model using MU, ALS, HALS, APG, and ANLS method. In Section IV, we introduce our strategy to speed up the computation of sparse NCP. Section V describes the detailed experiments on synthetic and real datasets. Finally, we discuss several key issues related to sparse NCP in Section VI and conclude our paper in Section VI.


 



§ SPARSE NONNEGATIVE CANDECOMP/PARAFAC DECOMPOSITION

In this paper, operator ∘ represents outer product of vectors, ⊙ represents the Khatri-Rao product,  ∗ represents the Hadamard product that is the elementwise matrix product, ⟨ ⟩ represents inner product, and  represents Kruskal oprator. _F denotes Frobenius norm, and _1 denotes l_1-norm. Basics of tensor computation and multi-linear algebra can be found in review <cit.>.

We present the designed sparse nonnegative CANDECOMP/PARAFAC decomposition (sparse NCP) in this section.

Given a nonnegative Nth-order tensor 𝒳∈ℝ^I_1× I_2×…× I_N and a positive number R, we design the sparse NCP as the following minimization problem:

    min_A^(1),…,A^(N)   𝒪= 1/2𝒳-A^(1),…,A^(N)_F^2
        +∑_n=1^Nα_n/2A^(n)_F^2+∑_n=1^Nβ_n∑_r=1^Ra_r^(n)_1
    s.t.    A^(n)⩾ 0  for  n=1, …,N,

where A^(n)∈ℝ^I_n× R for n=1,…,N are the estimated factors in different modes, α_n and β_n are positive regularization parameters in parameter vectors α∈ℝ^N× 1 and β∈ℝ^N× 1, I_n is the size in mode-n, a_r^(n) represents the rth column of A^(n), and R is the initial number of components.

Let X_(n)∈ℝ^I_n ×∏_ñ=1,ñ≠ n^N I_ñ represent the mode-n unfolding (matricization) of original tensor 𝒳. And the mode-n unfolding of the estimated tensor in Kruskal operator A^(1),…,A^(N)  can be written as A^(n)(B^(n))^T, in which B^(n)=(A^(N)⊙⋯⊙A^(n+1)⊙A^(n-1)⊙⋯⊙A^(1))∈ℝ^∏_ñ=1,ñ≠ n^N I_ñ× R. In block coordinate descent framework, factor A^(n) is updated alternatively by a subproblem in every iteration, which equals to the following minimization problem:

    min_A^(n)ℱ(A^(n)) =    1/2X_(n)-A^(n)(B^(n))^T_F^2
        +α_n/2A^(n)_F^2+β_n∑_r=1^Ra_r^(n)_1
    s.t.    A^(n)⩾ 0.

As previously mentioned, the sparse NCP problem in (<ref>) is non-convex, therefore finding its global minimum is NP-hard. However, the subproblem (<ref>) with Frobenius norm and l_1-norm is convex <cit.>, which is the key point to solve (<ref>). All the optimization methods introduced in the introduction can be applied to (<ref>). Furthermore, the objective function in (<ref>) can be represented as

    ℱ(A^(n))     = 1/2tr[ X_(n)^T X_(n)] - tr[(A^(n))^T X_nB^(n)]
        + 1/2tr[   A^(n)(B^(n))^T    B^(n)(A^(n))^T    ]
        + α_n/2tr[ (A^(n))^T  A^(n)] + β_n tr[ E^T A^(n)],

where E∈ℝ^I_n × R is a matrix whose elements are all equal to 1. The partial gradient (or partial derivative) of ℱ with respect to A^(n) is always used during computation,

    ∇_A^(n)ℱ(A^(n)) = ∂/∂A^(n)ℱ(A^(n))
        = -X_(n)B^(n) + A^(n)( B^(n))^T B^(n) + α_n A^(n) + β_n E
        =  A^(n)[ ( B^(n))^T B^(n) + α_n I_R ] -X_(n)B^(n) + β_n E,

where X_(n)B^(n) is called the Matricized Tensor Times Khatri-Rao Product (MTTKRP) <cit.>. The item ( B^(n))^T B^(n) can be computed efficiently by

    ( B^(n))^T B^(n) = [ ( A^(N))^T A^(N)] ∗⋯
       ∗[ ( A^(n+1))^T A^(n+1)] ∗[ ( A^(n-1))^T A^(n-1)]
       ∗⋯∗[ ( A^(1))^T A^(1)].




§ OPTIMIZATION METHODS FOR SOLVING SPARSE NCP

In this section, we present the solutions to the sparse NCP problem in (<ref>) by all the optimization methods of MU, ALS, HALS, APG, and ANLS. To the best of our knowledge, there are no existing solution proposed directly for (<ref>). Some optimization methods had been employed to solve NMF and NCP with Frobenius norm or l_1-norm regularization item separately <cit.>, so it is natural to combine the solution for these two items together.


 §.§ Multiplicative Update

Multiplicative update (MU) was first proposed by Lee et al for NMF <cit.>. Cai et al. proposed a straightforword way using lagrange multiplier to solve NMF subproblems <cit.>, where the same update rules can be obtained as Lee's method. We extend Cai's method to tensor case. Represent tensor factor A^(n) elementwisely by A^(n)=[ a_ir^(n)] for i=1,…,I_n and r=1,…,R. Let ψ_ir be the Lagrange multiplier for constraint a_ir^(n)⩾ 0, and Ψ=[ψ_ir], based on (<ref>) the Lagrange ℒ is

    ℒ =    ℱ(A^(n))     +    tr[Ψ(A^(n))^T]
    
    =    1/2tr[ X_(n)^T X_(n)] - tr[(A^(n))^T X_nB^(n)]
        + 1/2tr[   A^(n)(B^(n))^T    B^(n)(A^(n))^T    ]
        + α_n/2tr[ (A^(n))^T  A^(n)] + β_n tr[ E^T A^(n)]
        + tr[Ψ(A^(n))^T].

The partial derivative of ℒ with respect to A^(n) is

    ∂ℒ/∂A^(n) =     -X_(n)B^(n) + A^(n)( B^(n))^T B^(n)
        + α_n A^(n) + β_n E + Ψ.

Using KKT condition ψ_ira_ir^(n)=0, we obtain the following equation for a_ir^(n):

    -[ X_(n)B^(n)]_ir a_ir^(n) + 
       [ A^(n)( B^(n))^T B^(n) + α_n A^(n) + β_n E]_ir a_ir^(n) = 0.

This equation leads to the following multiplicative updating rule:

    a_ir^(n)← a_ir^(n)[ X_(n)B^(n)]_ir/[ A^(n)(  ( B^(n))^T B^(n) + α_n I_R ) + β_n E]_ir,

where I_R∈ℝ^R× R is a identity matrix.
The implementation of MU method is listed in Algorithm <ref>.




 §.§ Alternating Least Squares

Cichochi et al. comprehensively introduced the ALS method to solve NMF and NCP problems with diverse regularization items <cit.>. ALS method is quite easy to implement for the sparse NCP problem (<ref>).

Let the partial derivative of the objective function ℱ(A^(n)) in (<ref>) equal to 0,

    ∂/∂A^(n)ℱ(A^(n))
        = A^(n)[ ( B^(n))^T B^(n) + α_n I_R ] - [ X_(n)B^(n) - β_n E] 
        = 0,

we have the following updating rule for A^(n),

    A^(n)←[ X_(n)B^(n) - β_n E/( B^(n))^T B^(n) + α_n I_R]_+ ,

where [ ]_+ is a half-wave rectifying nonlinear projection to enforce nonnegativity. Algorithm <ref> shows the implementation of ALS method.




 §.§ Hierarchical Alternating Least Squares

Hierarchical alternating least squares (HALS) is a method to update each factor column by column. For the sake of simplification, we use a_r and b_r instead of a_r^(n) and b_r^(n) in this part, which are the rth column of A^(n) and B^(n) respectively. We also use [ A^(n)]_(:,r)=a_r∈ℝ^I_n× 1 to represent the column of a matrix, and [ A^(n)]_(i,r)=a_ir^(n) to represent an element in a matrix. The objective function ℱ in (<ref>) can be represent as

    ℱ=1/2X_(n)-∑_r=1^Ra_rb_r^T_F^2+α_n/2∑_r=1^Ra_r_2^2+β_na_r_1.

The minimization problem in (<ref>) can be solved iteratively by subproblems of columns:

    min_a_rℱ_r=1/2Z_r-a_rb_r^T_F^2+    α_n/2a_r_2^2+β_na_r_1
    s.t. a_r⩾ 0   ,

for r=1,…,R, and

    Z_r=X_(n)-∑_r̃=1,r̃≠ r^Ra_r̃b_r̃^T.

The partial derivative of ℱ_r with respect to a_r is

    ∂ℱ_r/∂a_r    =(a_rb_r^T-Z_r)b_r+α_na_r+β_n1,
        = (b_r^Tb_r+α_n )a_r-(Z_rb_r-β_n1),

where 1∈ℝ^I_n× 1 is a vector with all elements equaling to 1. When ∂ℱ_r/∂a_r=0, nonnagetive column vector a_r can be updated as

    a_r←[ Z_rb_r-β_n1]_+/b_r^Tb_r+α_n,

which is a closed-form solution <cit.>.

A fast HALS method was proposed to solve large-scale NMF problem <cit.>. We use the same idea to solve the sparse NCP problem. Z_r in (<ref>) can also be represented as

    Z_r=X_(n)-∑_r̃=1^Ra_r̃b_r̃^T+a_rb_r^T.

Replacing Z_r in (<ref>) by (<ref>), we obtain the new efficient update rule for a_r as is shown in (<ref>).
HALSUpdate


However, as mentioned above, a_r in different factors is usually poorly scaled, therefore after updating in (<ref>) it has to be normalized as

    a_r^(n)←a_r^(n)/a_r^(n)_2,

for n=1,…,N-1, and r=1,…,R. The procedures of HALS are illustrated in Algorithm <ref>.




 §.§ Alternating Proximal Gradient

The mathematical properties of alternating proximal gradient (APG) method were thoroughly analyzed in Xu's works <cit.>. APG method has exhibited excellent performances on both NMF and NCP problems, which is also efficient to cope with l_1 sparse regularization <cit.>. Supposing updating A^(n) in (<ref>) at the kth iteration, APG is computed as the following.

Calculate block-partial gradient of ℱ(A^(n)) in (<ref>) as ∇_A^(n)ℱ(A^(n))=∂/∂A^(n)ℱ(A^(n)). We take

    L_k-1^(n)=( B_k-1^(n))^T B_k-1^(n)_2

as Lipschitz constant of ∇_A^(n)ℱ with respect to A^(n), where A_2 is spectral norm of matrix. 








We also select

    ω_k-1^(n)=min( ω̂_k-1,δ_ω√(L_k-2^(n)/L_k-1^(n))),

where δ_ω<1 is predefined, and ω̂_k-1=t_k-1-1/t_k with t_0=1, t_k=1/2(1+√(1+4t_k-1^2)).

Let

    A_k-1^(n)=A_k-1^(n)+ω_k-1^(n)(A_k-1^(n)-A_k-2^(n))

denote an extrapolated point where ω_k-1^(n) is the extrapolation weight, and let

    G_k-1^(n)=    A_k-1^(n)[ ( B_k-1^(n))^T B_k-1^(n) + α_nI_R ]
        - X_(n)B^(n)+ β_n E

represent the block-partial gradient of ℱ(A^(n)) at A_k-1^(n). Factor A^(n) at iteration k is updated by

    A_k^(n)=_A^(n)⩾ 0   ⟨G_k-1^(n),A^(n)-A_k-1^(n)⟩
        + L_k-1^(n)/2A^(n)-A_k-1^(n)_F^2.

The closed form of (<ref>) can be written as

    A_k^(n) =max(0,     A_k-1^(n)-G_k-1^(n)/L_k-1^(n))
    
    = max(       0,  A_k-1^(n)    -  A_k-1^(n)[ ( B_k-1^(n))^T B_k-1^(n) + α_nI_R ]/L_k-1^(n)
       +  X_(n)B^(n) - β_n E/L_k-1^(n)).


APG method for sparse NCP can be implemented by procedures in Algorithm <ref>.




 §.§ Alternating Nonnegative Least Squares


Alternating Nonnegative Least Squares (ANLS) is an important and efficient method for NMF problems <cit.>. Kim, H. and Kim, J. utilized ANLS to solve NMF with squared Frobenius-norm reqularization and squared l_1-norm sparse regularization <cit.>. This idea can be naturally extended to the following sparsity regularized NCP problem:

    min_A^(1),…,A^(N)   𝒪_ANLS= 1/2𝒳-A^(1),…,A^(N)_F^2
    
    +∑_n=1^N    α_n/2A^(n)_F^2+∑_n=1^Nβ_n/2∑_r=i^I_n[A^(n)]_(i,:)_1^2
    s.t.    A^(n)⩾ 0  for  n=1, …,N,

where [A^(n)]_(i,:) denotes the ith row of A^(n). Factors A^(n) can be updated alternatively by the following minimization subproblem:

    min_A^(n)ℱ_NNLS(A^(n)) =    1/2X_(n)-A^(n)(B^(n))^T_F^2
    
    +α_n/2A^(n)_F^2     +β_n/2∑_r=i^I_n[A^(n)]_(i,:)_1^2
    s.t.    A^(n)⩾ 0.

By some simple mathematical operations, (<ref>) is equivalent to the following minimization problem:

    min_A^(n)   ℱ_NNLS(A^(n)) =
       1/2[        B^(n);    √(α_n)I_R; √(β_n)1_1× R ](A^(n))^T-
    [  X_(n)^T; 0_R× I_n; 0_1× I_n ]_F^2.
       s.t. A^(n)⩾ 0,

Afterwards, the partial derivative of ℱ_NNLS(A^(n)) to A^(n) is

    ∂/∂A^(n)ℱ_NNLS(A^(n))
        =  A^(n)[ ( B^(n))^T B^(n) + α_n I_R + β_n E] -X_(n)B^(n).

From (<ref>) and (<ref>), it is clear to see that the subproblem in (<ref>) still satisfies the basic nonnegative least squares (NNLS) structure, which can be solve conveniently by those optimization methods for NNLS, such as active-set (AS) <cit.> and block principal pivoting (BPP) <cit.>.

In (<ref>), the squared l_1-norm of the rows in A^(n) is used to impose sparsity in order to satisfy the NNLS structure, which is different from the sparse NCP in (<ref>). Algorithm <ref> explicates the ANLS method for sparsity regularized NCP in (<ref>).





§ STOPPING CONDITION AND ACCELERATING STRATEGY

The optimization procedures for tensor decomposition are implemented by iterations. For NCP, a sequence of {A_k^(1), …, A_k^(N)}_k=1^∞ is produced at each iteration. It is necessary to terminate the iteration until some stopping condition is satisfied. Common stopping conditions includes the following: predefined maximum number of iterations, predefined maximum running time, the change of objective function value, the change of relative error (data fitting) <cit.>.



 §.§ Objective Function and Relative Error


The computations of objective function value and relative error (data fitting) are highly correlated. In the kth iteration, the objective function value of NCP problem (<ref>) is

    𝒪_NCPk = 1/2𝒳-A_k^(1),…,A_k^(N)_F^2,

and the relative error <cit.> is defined by

    RelErr_k=𝒳 - A_k^(1),…,A_k^(N)_F/𝒳_F.

Comparing (<ref>) and (<ref>), the relative error can also be computed from the objective function value directly:

    RelErr_k= √(2𝒪_NCPk)/𝒳_F.

Meanwhile, the data fitting can be computed by

    Fit_k=1-RelErr_k.


Based on the objective function value and the relative error, the stopping condition can be set by

    |RelErr_k-1-RelErr_k | < ϵ

or

    |𝒪_NCPk-1-𝒪_NCPk| < ϵ.

The threshold of ϵ can be set by a very small positive value, such as 1e-8.



 §.§ Accelerated Computation of Objective Function

By mode-n unfolding of tensor, the objective function of NCP in (<ref>) at the kth iteration can be represented equivalently as the following:

    𝒪_NCPk= 1/2X_(n)-A_k^(n)(B_k^(n))^T_F^2.

The works of <cit.> and <cit.> introduced a convenient idea to compute the objective function base on the trace computation of matrix. Inspired by this idea, we further represent the objective function in (<ref>) by

    𝒪_NCPk
       =1/2tr{X_(n) - A_k^(n)(B_k^(n))^T ^T X_(n) - A_k^(n)(B_k^(n))^T }
       =1/2{𝒳_F^2 - 2trA_k^(n)(X_(n)B_k^(n))^T 
       +tr( (A_k^(n))^T A_k^(n)) ( (B_k^(n))^T B_k^(n)) }.

Furthermore, the objective function equals to

    𝒪_NCPk=1/2{𝒳_F^2 - 2∑_j=1^R∑_i=1^I_nN_i,j + ∑_j=1^R ∑_i=1^R M_i,j},

where N=A_k^(n)∗(X_(n)B_k^(n))∈ℝ^I_n × R and M=( (A_k^(n))^T A_k^(n)) ∗( (B_k^(n))^T B_k^(n)) ∈ℝ^R × R. Here, ∗ is the Hadamard product.

Since tensor data usually consist of a large amount of data points, the computation of the objective function at each iteration will be time consuming by (<ref>) or (<ref>). For example, in (<ref>) the computational complexity of A_k^(n)(B_k^(n))^T is O(R ×∏_n=1^N I_n). On the other hand, by observing Algorithm <ref> to <ref>, we find that the items of X_(n)B_k^(n) and (B_k^(n))^T B_k^(n) have been computed in advance in order to update A_k^(n). Therefore, these two items can be employed directly to compute the objective function by (<ref>). In (<ref>), the computational complexity of N and M is only O(I_nR+R^2), which has been reduced significantly.



 §.§ Computation of Objective Function for Sparse NCP

The objective function and relative error for sparse NCP should be carefully considered, due to the extra Frobenius norm and l_1-norm regularization items. Based on the relationship introduced in Subsection IV-A and the accelerating method in Subsection IV-B, the computation of objective function and relative error for sparse NCP is explicated in Algorithm <ref>.



For ANLS method, the step 13 in Algorithm <ref> should be changed into

    𝒪_k=𝒪_k + α_n/2A_k^(n)_F^2 + β_n/2∑_r=i^I_n[A_k^(n)]_(i,:)_1^2

according to (<ref>).




§ EXPERIMENTS AND RESULTS

We carried out the experiments on synthetic tensor data and real tensor data, both of which contain third-order and fourth-order cases. We compared the the abilities of sparse NCP methods to impose sparsity implemented by MU, ALS, HALS, APG, ANLS-AS and ANLS-BPP.

In all sparse NCP experiments, the factor matrices were initialized using nonnegative normally distributed random numbers by MATLAB function . We keep the Frobenius norm regularization in the sparse NCP model to make a fair comparison for all methods, especially the comparison of the objective function, although it is not necessary for some optimization methods, such as ANLS-AS, APG, HALS and MU. We set the Frobenius norm parameters by α_1=α_2=, …, =α_N=1e-6. The l_1-norm regularization parameters of β_n, n=1,…,N, in sparse NCP are the key elements to impose sparsity, which are the most important testing parameters in the experiments. In order to make it convenient to select and test the parameters, we also kept β_n, n=1,…,N, the same in all modes of the tensor. After selecting the β_n, we calculated and evaluated the sparsity level <cit.> of the factor matrices by

    Sparsity_A^(n)=#{A^(n)_i,r<T_s }/I_n × R,

where #{·} denotes the number of elements that are smaller than the threshold T_s in factor matrix A^(n). In the experiments, we selected T_s=1e-3.

All the experiments were conducted on computer with Intel Core i5-4590 3.30GHz CPU, 8GB memory, 64-bit Windows 10 and MATLAB R2016b. The fundamental tensor computation was based on Tensor Toolbox 2.6 <cit.>.



 §.§ Third-Order Synthetic Data


In the first experiment, we constructed a synthetic third-order tensor by 10 channels of simulated sparse and nonnegative signals[The sparse signals come from the file of  included in NMFLAB, which can be downloaded from http://www.bsp.brain.riken.jp/ICALAB/nmflab.htmlhttp://www.bsp.brain.riken.jp/ICALAB/nmflab.html], as shown in Fig. <ref>(a). There are 1000 points in each channel, so the signal matrix is S^(1)=[s_1,…,s_10]∈ℝ^1000 × 10. Two uniformly distributed random matrices A^(2), A^(3)∈ℝ^100 × 10 were employed as mixing matrices, which were generated by  function in MATLAB. Afterwards, we synthesized the third-order tensor by 𝒳_Syn-3rd=S^(1),A^(2),A^(3)∈ℝ^1000 × 100 × 100. Next, nonnegative Gaussian noise was added to the tensor with SNR of 40dB, which was generated by MATLAB code .









For all sparse NCP methods, we used the stop condition of relative error change in (<ref>), in which the threshold is ϵ=1e-8. The maximum running time was set by 180 seconds. We selected a larger value of 20 as the number of components for tensor decomposition[Since 10 channels of signals are mixed in the tensor, naturally, 10 should be selected as the optimal component number. The number of components might also be estimated by some classical methods, such as DIFFIT <cit.>. However, we selected 20 according to the purpose of the experiment.]. The reason is that we intend to recover the 10 channels of true signal just by imposing sparse regularization during decomposition, even though we don't know the exact optimal number of components. We selected values of β_n=0,0.1,0.5,1,2,3 for all the optimization methods to evaluate their abilities to impose sparsity. The selection of sparse regularization parameters depends on the tensor data. After tensor decomposition, the values of objective function value, relative error, running time, iteration number, the sparsity level and nonzero component number of the estimated signal factor matrix were recorded as the performance evaluation criteria. For all optimization methods with each β_n, the sparse NCP was run 30 times, and the average values of all criteria were computed. The results are shown in Table <ref>.

From Table <ref>, it can be found that all ANLS-AS, ANLS-BPP, APG, MU, ALS methods can impose sparsity with proper sparse regularization parameter β_n. With certain sparse regularization, 10 nonzero components are retained in the mode-1 factor matrix, which represent the 10 channels of sparse signals extracted from the mixed tensor. ANLS-AS, ANLS-BPP, APG, and ALS exhibit comparatively low relative errors and less running time, whereas MU converges very slowly. On the other hand, HALS fails to impose sparsity even with a larger β_n and can hardly reduce the number of nonzero components. HALS also has very large objective function value and high relative error with larger β_n.

It can be inferred from Table <ref> that, after properly tuning the sparse regularization parameter β_n, weak components will be removed (set to 0), weak elements in strong components will be prohibited, and the true 10 channels of sparse signals will be recovered. One of the recovered sparse signal matrix by ANLS-BPP is shown in Fig. <ref>(b). Afterwards, the accurary of the recovered signals should be evaluated. Let T^(1)=[t_1,…,t_R̃]∈ℝ^L ×R̃ represent the estimated matrix of sparse signals, in which R̃ is the number of nonzero components and L is the length of component (L=1000 in this experiment). We evaluate the accuracy of the estimated matrix T^(1) compared with original sparse signals S^(1) by Peak Signal-to-Noise-Ratio (PSNR, see Chapter 3 in <cit.>)

    PSNR=1/R̃∑_r=1^R̃10log_10( L/t̂_r-ŝ_c_2^2),

where t̂_r is the rth normalized estimated sparse signal, and ŝ_c is the normalized reference sparse signal. ŝ_c comes from S^(1), which has the highest correlation coefficient with t̂_r. For all the optimization methods with each β_n, all the PSNR values were recorded after 30 times of running of sparse NCP. Subsequently, box plot of PSNR for MU, ALS, HALS, APG, ANLS-AS, ANLS-BPP was drawn in Fig. <ref>.

In Fig. <ref>, it is clear that APG and ALS methods have higher PSNR at larger sparse regularization parameters of β_n=2,3, with which they recover the 10 channels of sparse signals more precisely. For ANLS-AS and ANLS-BPP, their PSNR values decrease when β_n>0.5. The reason might be that ANLS framework uses the squared l_1-norm as the sparse regularization item, which is more sensitive to the sparse regularization parameter β_n (smaller β_n will cause higher sparsity, and larger β_n might spoil the decomposition). Fig. <ref> also shows that HALS performs poorly with sparse sparse regularization.

We also recorded the objective function values of all sparse NCP methods within the first 30 seconds. The results of these methods with β_n=0,0.1,1,2 are shown in Fig. <ref>. In the results, we find that MU converges very slowly in all cases, but it can minimize the objective function to a low level gradually. ALS method is not stable, which sometimes can't ensure that the objective function decreases. However, the stability of ALS will be improved with higher sparse regularization parameters. HALS runs much fast, but its objective function becomes very large when sparse regularization is added. APG shows excellent performances in both running speed and the ability to minimize the objective function in the third-order tensor case. ANLS-AS and ANLS-BPP can minimize the objective function very fast at the beginning, but their objective function values are a bit higher when sparse regularization is added. This is due to the squared l_1-norm item in ANLS framework.







 §.§ Fourth-Order Synthetic Data


In the second experiment, we synthesized a fourth-order tensor by the same 10 channels of simulated sparse and nonnegative signals as that in last section. The tensor was synthesized by 𝒳_Syn-4th=S^(1),A^(2),A^(3),A^(4)∈ℝ^1000 × 100 × 100 × 5, in which S^(1) is the simulated sparse signals, and A^(2),A^(3)∈ℝ^100 × 10,A^(4)∈ℝ^5 × 10 are random matrices in uniform distribution. Nonnegative Gaussian noise was added to this fourth-order tensor with SNR of 40dB.

Processing higher-order (order ⩾ 4) is very challenging due to the huge amount of data and the high complexity of optimization method. In this fourth-order tensor experiment, we selected ϵ=1e-7 as the threshold of stopping condition based on the relative error change in (<ref>), and set 1800 seconds as the maximum running time. Other parameter settings and experimental methods are the same as those in subsection <ref>. The averages of objective function, relative error, running time, iteration number, the sparsity level and nonzero component number of the estimated signal factor matrix after 30 times of running are recorded in Table <ref>. The box plot of PSNR for all methods after 30 times of run is shown in Fig. <ref>. The results of the fourth-order synthetic tensor decomposition in Table <ref> and Fig. <ref> have proven once more the effectiveness of MU, ALS, APG, ANLS-AS, ANLS-BPP to impose sparsity and estimate the true sparse signals.

We recorded the objective function values of all sparse NCP methods for the fourth-order tensor within the first 600 seconds. The results of these methods with β_n=0,0.1,1,2 are shown in Fig. <ref>. It is clear to see that ALS method is still not stable, which can't guarantee the objective function to decrease. MU method still shows slow convergence. Surprisingly, the APG method that presents excellent performance for third-order tensor runs very slowly for the fourth-order tensor, which performs even more poorly than MU at the first tens of seconds. HALS has very large objective function value when sparse regularization is imposed. However, both ANLS-AS and ANLS-BPP still converge very fast for the fourth-order sparse NCP. When sparse regularization is imposed, ANLS-AS and ANLS-BPP have slightly higher objective function values due to the squared l_1-norm item in ANLS framework.





 §.§ Third-Order ERP Data


In the third experiment, we utilized an open preprocessed ERP tensor[Data website: http://www.escience.cn/people/cong/AdvancedSP_ERP.htmlhttp://www.escience.cn/people/cong/AdvancedSP_ERP.html]. Original data are organized in fourth-order form, whose size is channel × frequency × time × subject-group = 9 × 71 × 60 × 42. The 9 channel points denote the 9 electrodes on the scalp, the 71 frequency points show the spectrum within 1-15Hz, the 60 time points illustrate the temporal energy between 0-300ms, and the 42 subject-group points include 21 subjects with reading disability (RD) and 21 subjects with attention deficit (AD) <cit.>. In this experiment, we merged the modes of frequency and time, by which the original tensor was reshaped into a third-order tensor with size channel × frequency-time × subject-group = 9 × 4260 ×  42.

For this third-order tensor, the number of components is set by 40 according to previous study <cit.>. For all methods, the threshold of stopping condition was set by 1e-8 based on the relative error change, and the maximum running time was 240s. The values of β_n=0,10,50,100,200,300 were tested for all methods. We recorded the objective function, relative error, running time, iteration number, the sparsity level and nonzero component number of the frequency-time factor matrix for all methods. The average values after 30 times of run are recorded in Table <ref>.

It can be found from Table <ref> that ANLS-AS, ANLS-BPP and APG are very effective to impose sparsity on the factor matrix and reduce components number by adjusting β_n, which costs comparatively less time. MU is not very effective to impose sparsity, which often reaches the limit of running time with slow convergence. Both HALS and ALS fail to impose proper sparsity by adjusting β_n, which also show high objective function values and relative errors with large β_n.

We also recorded the objective function values of all methods within the first 120 seconds. The results of these methods with β_n=0,10,50,100 are displayed in Fig. <ref>. Fig. <ref> shows that MU still converges very slowly and ALS is not stable. When sparse regularization is added, HALS has very large objective function value. ANLS-AS and ANLS-BPP converge fast at the beginning, but have slightly higher objective function values with β_n > 0 due to the squared l_1-norm regularization item. APG has excellent performances in both speed and convergence for this third-order ERP tensor.









 §.§ Fourth-Order ERP Data


In the fourth experiment, we applied all the sparse NCP algorithms to the original fourth-order ERP tensor introduced in subsection <ref>. For this experiment, the maximum running time was set by 120 seconds. Other settings are the same as those in the third-order case. We recorded the values of objective function, relative error, running time, iteration number, the sparsity level and nonzero component number of the temporal factor matrix. The average values after 30 times of run are recorded in Table <ref>. We also recorded the objective function values of all methods with β_n=0,10,50,100 within the first 30 seconds as shown in Fig. <ref>.

For ANLS-AS, ANLS-BPP, MU, HALS and ALS, the results of the fourth-order ERP tensor in Table <ref> and Fig. <ref> reveal similar properties as those in the third-order case. For the fourth-order case, MU is not very effective to impose sparsity on the factor matrices, whose sparsity are very close to the results without sparse regularization. Hence MU has smaller relative error and lower objective function value as shown in Fig. <ref>. However, APG minimizes the objective function slowly even than MU for the fourth-order ERP tensor, which is quite different from the situation in the third-order case.

Despite the slow decrease of the objective function, APG can guarantee to converge and has exhibited effectiveness to impose sparsity. Two groups of components extracted by APG method without and with sparse regularization are shown in Fig. <ref>. It is clear to see from the temporal components that, with proper sparse regularization parameter β_n (β_n=10 in this case), some weak and redundant elements are suppressed.





§ DISCUSSION

From the results and analyses of the sparse NCP experiments, we have the following findings.

MU is a common method for NMF, although it shows slow convergence. It can also be extended to NCP due to the flexibility to handle regularization. However, MU still convergences very slowly in NCP compared with other methods, which sometimes is not very sensitive to sparse regularization.

ALS method is very simple to implement, but it is not stable for NCP with sparse regularization. By ALS, the nonnegative constraint is just obtained by projecting all negative elements of factor matrices to zeros in the subproblems. Hence the solution by projecting is not an accurate solution to subproblem. Therefore, ALS method can not guarantee the convergence of NCP.

HALS is a very efficient method for NCP. Nevertheless, the required normalization procedures will complicate the optimization problem. In our experiments, we find that HALS does not perform well to impose sparsity. High objective function value and relative error occur with large sparse regularization parameters due to the normalization procedures.

APG has proved to be a convergent and stable method for tensor decomposition due to the proximal gradient, which exhibits high efficiency for third-order tensor. Meanwhile, APG can also handle sparse regularization flexibly, which yields closed-form solution. In spite of the convergence property, APG turns out very slowly for the fourth-order tensor case. In future, further studies are needed to accelerate APG for higher-order (order⩾4) nonnegative tensor decomposition.

ANLS framework converges fast and stably for sparse NCP benefiting from the advantages of many NNLS optimization methods, such as Active Set (AS) and Block Principal Pivoting (BPP). In addition, ANLS exhibits excellent performance to impose sparsity by tuning the sparse regularization parameter β_n. However, in order to be compatible with the ANLS framework, squared l_1-norm has to be employed as the sparse regularization item. Therefore, the method to impose sparsity in ANLS is slightly different from those in other sparse NCP methods. It is very interesting to find the way to solve sparse NCP with non-squared l_1-norm by NNLS, which will further minimize the objective function value as we believe.

At last, we want to mention the convergence properties of above optimization methods. The Frobenius norm regularization in sparse NCP (<ref>) can prevent rank deficiency of the matrix B^(n) in subproblem (<ref>), which will improve the convergence of sparse NCP to a globally optimal solution <cit.>. ANLS-BPP and ALS method require the B^(n) to be of full rank <cit.>, while the methods of ANLS-AS, APG, HALS and MU don't require the full rank condition. Both ANLS framework and APG method have very good convergence properties, which have been explained in <cit.> and <cit.> respectively. The subproblem of HALS method has closed form solution <cit.>, but the normalization procedures might spoil the convergence property of HALS <cit.>. In addition, the convergence of MU and ALS can't be guaranteed. We don't make deep analysis of the convergence properties of these optimization methods in this paper. More information about the convergence properties of the general block coordinate descent (BCD) method can be found in <cit.>.



§ CONCLUSION

In this paper, we investigated CANDECOMP/PARAFAC tensor decomposition with both nonnegative constraint and sparse regularization (sparse NCP). The methods of MU, ALS, HALS, APG and ANLS in block coordinate descent framework were deeply analyzed to solve sparse NCP. We compared all these methods by experiments on synthetic and real tensor data, both of which contain third-order and fourth-order cases. We find that APG and ANLS methods are stably convergent and highly effective to impose sparsity for the tensor decomposition compared with other methods. Meanwhile, ANLS framework is very efficient to process higher-order (order⩾4) tensor data. The proposed accelerated method to compute the objective function and relative error had further improved the efficiency of the sparse NCP.



























































IEEEtran










































