Sparse nonnegative tensor decomposition using proximal algorithm and inexact block coordinate descent scheme

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Sparse nonnegative tensor decomposition using proximal algorithm and inexact block coordinate descent scheme © The Author(s) 2021 Published version Wang, Deqing; Chang, Zheng; Cong, Fengyu Wang, D., Chang, Z., & Cong, F. (2021). Sparse nonnegative tensor decomposition using proximal algorithm and inexact block coordinate descent scheme. Neural Computing and Applications, 33(24), 17369-17387. https://doi.org/10.1007/s00521-021-06325-8 2021 ORIGINAL ARTICLE Sparse nonnegative tensor decomposition using proximal algorithm and inexact block coordinate descent scheme Deqing Wang 1,2 •Zheng Chang 2,3 •Fengyu Cong 1,2,4,5 Received: 14 October 2020 / Accepted: 8 July 2021 The Author(s) 2021 Abstract Nonnegative tensor decomposition is a versatile tool for multiway data analysis, by which the extracted components are nonnegative and usually sparse. Nevertheless, the sparsity is only a side effect and cannot be explicitly controlled without additional regularization. In this paper, we investigated the nonnegative CANDECOMP/PARAFAC (NCP) decomposition with the sparse regularization item using l1-norm (sparse NCP). When high sparsity is imposed, the factor matrices will contain more zero components and will not be of full column rank. Thus, the sparse NCP is prone to rank deficiency, and the algorithms of sparse NCP may not converge. In this paper, we proposed a novel model of sparse NCP with the proximal algorithm. The subproblems in the new model are strongly convex in the block coordinate descent (BCD) framework. Therefore, the new sparse NCP provides a full column rank condition and guarantees to converge to a stationary point. In addition, we proposed an inexact BCD scheme for sparse NCP, where each subproblem is updated multiple times to speed up the computation. In order to prove the effectiveness and efficiency of the sparse NCP with the proximal algorithm, we employed two optimization algorithms to solve the model, including inexact alternating nonnegative quadratic programming and inexact hierarchical alternating least squares. We evaluated the proposed sparse NCP methods by experiments on synthetic, real-world, small-scale, and large-scale tensor data. The experimental results demonstrate that our proposed algorithms can efficiently impose sparsity on factor matrices, extract meaningful sparse components, and outperform state-of-the-art methods. Keywords Tensor decomposition Nonnegative CANDECOMP/PARAFAC decomposition Sparse regularization  Proximal algorithm Inexact block coordinate descent 1 Introduction 1.1 Background Nonnegative tensor decomposition is a powerful tool in signal processing and machine learning [10,35]. Nonnegative CANDECOMP/PARAFAC (NCP), as an important decomposition method, has been widely applied to processing multiway data, such as hyperspectral data [39], electroencephalograph (EEG) data [11], fluorescence excitation-emission matrix (EEM) data [13], neural data [46], and many other multiway tensor data [30]. In many cases, the extracted components by NCP are not only nonnegative but also sparse. For example, the spectral components from EEG tensor decomposition are usually very sparse, representing the narrow-band frequencies of some brain activities [11]. For another example, after decomposing EEM tensor, a component in the sample mode denotes the concentrations of a compound in all samples [5], which is sometimes also sparse. The nonnegative constraint in NCP will naturally lead to sparse results. However, this sparsity is only a side effect, which cannot be controlled to a certain level [18]. Without properly controlling the sparsity, the intrinsic components in the data cannot be extracted precisely, especially in low signal-to-noise ratio conditions. Therefore, in order to extract meaningful and accurate sparse components, additional sparse regularization is necessary for NCP tensor decomposition. The design of NCP decomposition with explicit sparse regularization (sparse NCP) will benefit a lot from the methods in nonnegative matrix factorization (NMF) cases. Extended author information available on the last page of the article 123 Neural Computing and Applications https://doi.org/10.1007/s00521-021-06325-8(0123456789().,-volV)(0123456789().,-volV) On the one hand, an early study of NMF [18] proposed the method of projecting components into sparse vectors at some sparsity level. However, this method keeps all components at the same fixed sparsity level, which is not in line with the true sparsity of different components in the data. On the other hand, incorporating sparse regularization items into the optimization model is a popular method. The l1-norm is a conventional and effective regularizer to impose sparsity for signal processing [6]. The reason is that, for most underdetermined linear equations, the optimization problem with l1-norm regularization can yield strong sparsity [12]. More information about the sparse regularization can be found in [2,34,50]. Many works have been devoted to the tensor decomposition with sparse regularization, but only a few can be found for NCP. The works of [1,21] and [29] studied the sparse regularization for tensor decomposition using l1- norm and trace norm, but they only focused on the unconstrained CP model without the nonnegative constraint. The works of [14,28,31] and [47] proposed the methods of imposing sparsity by the l1-norm on nonnegative Tucker decomposition. However, these methods are not suitable for large-scale problems [47], and their effectiveness is unknown to NCP. Kim et al. considered solving sparse NCP using ANLS [23]. Nevertheless, ANLS seriously suffers from rank deficiency caused by high sparsity or zero components in the factor matrices. Recently, Huang et al. have proposed an alternating optimization-based ADMM (AO-ADMM) method, which can handle the l1-norm regularization item in NCP [19]. Nevertheless, there is no experimental evaluation on the sparse NCP in [19]. The work [32] proposed a sparse NCP algorithm, which is targeted at the multiway co-clustering. In practical applications, the sparse NCP may face the following two major challenges. One challenge is that when the tensor data are highly sparse or strongly sparse regularization is imposed on the decomposition, more and more zero components will appear in the factor matrices. Thus, the factor matrices are not of full column rank, which will cause the rank deficiency problem. The rank deficiency will further cause a poor convergence of the tensor decomposition algorithm. It is introduced that the proximal algorithm is an excellent method to improve the convergence of a mathematical optimization method [4]. In an optimization problem by iterations, the proximal algorithm is constructed by adding a proximal regularization item to the original model. This proximal item is the squared Frobenius norm of the difference between the current variable and its value in previous iteration [4]. The proximal algorithm can naturally be incorporated into tensor decomposition [26]. The other challenge is that, for large-scale tensor data, the process of sparse NCP decomposition might be inefficient. It is reported that the inexact block coordinate descent scheme could accelerate the convergence and is very beneficial to the large-scale problem [15,40]. Hence, the inexact scheme can be employed in the sparse NCP problem. 1.2 Contribution Firstly, in this paper, we propose a novel sparse NCP method with the l1-norm and the proximal algorithm. The proposed sparse NCP will overcome the rank deficiency and guarantee the decomposition to converge to a stationary point. The block coordinate descent (BCD) is one of the main techniques for tensor decomposition, especially the constrained one [24]. In BCD framework, each factor matrix is updated as a subproblem alternatively while other factor matrices are fixed. By the proximal algorithm, the proximal regularization item can make the subproblems strongly convex [26] and can provide a full column rank condition for the sparse NCP. Secondly, we develop an inexact BCD scheme for the novel sparse NCP model. The inexact scheme will speed up the computation of the sparse NCP, especially in largescale cases. Specifically, in the inexact BCD scheme, the subproblem of the sparse NCP is iterated multiple times for updating a factor matrix. Thirdly, in order to prove the viability of the sparse NCP model with the proximal algorithm and the inexact scheme, we employ two efficient optimization algorithms to solve the model, including inexact alternating nonnegative quadratic programming and inexact hierarchical alternating least squares. We evaluate the proposed sparse NCP methods on synthetic, real-world, small-scale and largescale tensor data. By properly selecting and tuning the sparse regularization, the effectiveness and efficiency of the sparse NCP methods are demonstrated to impose sparsity on factor matrices. 1.3 Organization The rest of this paper is organized as follows. Section 2 introduces some preliminaries. In Sect. 3, we describe the mathematical model of sparse NCP with the proximal algorithm and inexact BCD scheme. Section 4elucidates the solutions to the sparse NCP model using the optimization methods. Section 5describes the detailed experiments on synthetic and real-world datasets. Some critical observations are discussed in Sect. 6. Finally, we conclude our paper in Sect. 7. Neural Computing and Applications 123 2 Preliminaries In this paper, operator represents the outer product of vectors, represents the Khatri-Rao product, * represents the Hadamard product that is the elementwise matrix product, hi represents the inner product, []represents Kruskal operator and [ ] ? represents nonnegative projection. jj jjFdenotes Frobenius norm, and jj jj1denotes l1- norm. Basics of tensor computation and multi-linear algebra can be found in review papers [25,35]. 2.1 Nonnegative CP decomposition Given an Nth-order nonnegative tensor X2RI1I2IN and a positive number R, nonnegative CANDECOMP/ PARAFAC (NCP) is to solve the following minimization problem: min Að1Þ;...;AðNÞ 1 2jjXsAð1Þ;...;AðNÞtjj2 F s.t. AðnÞ>0 for n¼1;...;N;ð1Þ where AðnÞ2RInRfor n¼1;...;Nare the estimated factor matrices in different modes, Inis the size in mode-n, and Ris the initial number of components. We use FtensorA¼FtensorAð1Þ;...;AðNÞto denote the objective function in (8). The estimated factor matrices in Kruskal operator can be represented by the sum of Rrank-1 tensors in outer product form: sAð1Þ;...;AðNÞt¼X R r¼1 Yr¼X R r¼1 að1Þ raðNÞ r;ð2Þ where aðnÞ rrepresents the rth column of AðnÞ. Let XðnÞ2RInQN ~ n¼1;~ n6¼nI~ nrepresent the mode-nunfolding (matricization) of original tensor X. The mode-nunfolding of the estimated tensor in Kruskal operator sAð1Þ;...;AðNÞt can be written as AðnÞBðnÞT, in which BðnÞ¼AðNÞAðnþ1ÞAðn1ÞAð1Þ2RQN ~n¼1;~n6¼nI~nR . In BCD framework, factor AðnÞis updated alternatively by a subproblem in every iteration, which is equal to the following minimization problem: min AðnÞ FAðnÞ¼1 2jjXðnÞAðnÞBðnÞTjj2 F s.t. AðnÞ>0:ð3Þ The partial gradient (or partial derivative) of FAðnÞwith respect to AðnÞis o oAðnÞFAðnÞ¼AðnÞBðnÞTBðnÞXðnÞBðnÞ:ð4Þ In (4), the item XðnÞBðnÞis called the Matricized Tensor Times Khatri-Rao Product (MTTKRP) [35]. The item BðnÞTBðnÞcan be computed efficiently by BðnÞTBðnÞ¼AðNÞTAðNÞ Aðnþ1ÞTAðnþ1ÞAðn1ÞTAðn1Þ Að1ÞTAð1Þ: ð5Þ 2.2 Sparse regularization with l1-norm In order to impose sparsity to the factor matrices, it is natural to incorporate the sparse regularization items using l1-norm [9,47] into the objective function in (1), which leads to the following basic sparse NCP problem: min Að1Þ;...;AðNÞ 1 2jjXsAð1Þ;...;AðNÞtjj2 FþX N n¼1 bnX R r¼1jjaðnÞ rjj1 s.t. AðnÞ>0 for n¼1;...;N; ð6Þ where bnare positive sparse regularization parameters in parameter vectors b2RN1. The subproblem can be written as the following optimization problem min AðnÞ F0AðnÞ¼1 2XðnÞAðnÞBðnÞT2 FþbnX R r¼1aðnÞ r1 s.t. AðnÞ>0: ð7Þ In the objective function of the subproblem, the sparse regularization is imposed on the factor matrix AðnÞby the l1-norm. 3 The proposed sparse NCP model 3.1 Sparse NCP with proximal algorithm The basic sparse NCP in (6) has a serious drawback. When strongly sparse regularization is imposed in (6), many zero columns will appear in the factor matrices AðnÞ. Thus, both AðnÞand BðnÞcannot guarantee to be of full column rank. Therefore, the basic sparse NCP model in (6) will suffer from rank deficiency and cannot guarantee to converge. In order to overcome the above drawback, we propose the following sparse NCP model with proximal algorithm Neural Computing and Applications 123 (a proximal regularization item using squared Frobenius norm): min Að1Þ;...;AðNÞ1 2XsAð1Þ;...;AðNÞt2 F þX N n¼1 an 2e AðnÞAðnÞ2 FþX N n¼1 bnX R r¼1aðnÞ r1 s.t. AðnÞ>0 for n¼1;...;N; ð8Þ where e AðnÞis the value of the factor AðnÞin previous iteration during updating and anare positive regularization parameters in vectors a2RN1. In BCD framework, the subproblem of model (8) can be written in the following minimization problem: min AðnÞ FPROXAðnÞ¼1 2XðnÞAðnÞBðnÞT2 F þan 2e AðnÞAðnÞ2 FþbnX R r¼1aðnÞ r1 s.t. AðnÞ>0: ð9Þ The objective function FPROXAðnÞcan be further represented by the following form: FPROXAðnÞ¼1 2 XT ðnÞ ffiffiffiffiffi an pe AðnÞT 0 @1 A BðnÞ ffiffiffiffiffi an pIR ! AðnÞT  2 F þbnX R r¼1aðnÞ r1: ð10Þ In (10), it is clear to see that the item BðnÞ ffiffiffiffiffi an pIR  must be of full column rank even though BðnÞis not of full column rank. Thus, the proposed sparse NCP with the proximal algorithm can successfully overcome the rank deficiency problem in the objective function. 3.2 Inexact block coordinate descent scheme The BCD is a main framework to solve tensor decomposition. It is reported that the inexact BCD scheme could accelerate the computation [15,40]. Specifically, the factor matrices AðnÞ;n¼1;...;N, are updated alternatively in outer iterations; meanwhile, in the subproblem (9), the factor AðnÞis also updated several times in inner iterations. The procedures of the inexact scheme are listed in Algorithm 1. 3.3 Convergence analysis The proposed sparse NCP method in (8) can guarantee to converge to a stationary point. Proposition 1 Every limit point of the sequence Að1Þ k;...;AðNÞ k no 1 k¼1generated by the sparse NCP in Algorithm 1 is a stationary point of (6). Proof The objective function FPROXAðnÞin (9) with the proximal regularization item is strictly convex [4]. Moreover, FPROXAðnÞis a proximal upper bound [17,33]of the objective function F0AðnÞin (7). Using the inexact block coordinate descent scheme, the subproblem in Algorithm 1 is updated by a finite number of inner iterations. According to the Theorem 2 in [49], every limit point of the sequence Að1Þ k;...;AðNÞ k no 1 k¼1generated by the sparse NCP in Algorithm 1 is a stationary point of (6). h 4 Optimization methods for solving sparse NCP In order to prove the viability and effectiveness of the novel sparse NCP with the proximal algorithm and inexact scheme, we employ the following two optimization methods to solve the model. 4.1 Alternating nonnegative quadratic programming First, we utilize a method that is based on a general form of the alternating nonnegative least squares (ANLS). The classical ANLS is an important tool for NMF and NCP [24]. Many efficient optimization algorithms were proposed to solve the nonnegative least squares (NNLS) subproblems, such as active-set (AS) [20] and block principal pivoting (BPP) [22]. However, there are two limitations to the application of ANLS to sparse NCP. The first limitation is that ANLS is very prone to rank deficiency. The proximal algorithm can tackle this limitation in our sparse NCP model. The second limitation is that the Neural Computing and Applications 123 subproblem of our proposed sparse NCP model cannot be represented in a least squares form due to the l1-norm regularization, which can be clearly seen in (10). Therefore, some new forms of the objective function in (8) should be considered. Inspired by [27], the subproblem of the proposed sparse NCP in (9) can be represented in the nonnegative quadratic programming (NNQP) form as the following problem: min AðnÞX In i¼11 2AðnÞði;:ÞMAðnÞT ði;:ÞþNði;:ÞAðnÞT ði;:Þ þ1 2XðnÞði;:ÞXðnÞT ði;:Þþan 2e AðnÞði;:Þe AðnÞT ði;:Þ s.t. AðnÞ>0; ð11Þ where ði;:Þrepresents the ith row of a matrix, M¼BðnÞTBðnÞþanIR,N¼bnEXðnÞBðnÞane AðnÞand Eis a matrix of all ones. In fact, NNQP is a general form of NNLS. The above-mentioned optimization methods for NNLS can also be used to solve NNQP problem. In this study, we only use block principal pivoting (BPP) [22] as the NNQP solver, which has been proven to be a very efficient method [22,24]. The solver of BPP contains multiple inner iterations. We limited the inner iterations by several times in the inexact scheme. We name the method of solving tensor decomposition using NNQP as alternating nonnegative quadratic programming (ANQP). Furthermore, we abbreviated the method of solving the sparse NCP with the proximal algorithm using ANQP as PROX-ANQP. Algorithm 2 explicates the PROX-ANQP method. 4.2 Inexact hierarchical alternating least squares Second, we employ an inexact hierarchical alternating least squares (iHALS) method for solving the sparse NCP with the proximal algorithm. The conventional HALS is an efficient method of updating each factor column by column [7,9]. However, the HALS method has two major drawbacks to solving the sparse NCP. First, HALS is also very prone to rank deficiency. Specifically, if a column of the factor matrix AðnÞbecomes a zero vector, the HALS will break down [22]. One practical remedy is to replace the zero elements with a small positive value [9], such as 1016. However, by this modification, the obtained factor matrices are not sparse anymore. Second, HALS suffers from the caveat problem (see Section 5.2 in [22]). Specifically, the unbalanced scales will appear in different columns and factors. For example, one column in the first factor might have a scale of 108 and the corresponding column in the second factor might have a scale of 108. At the same time, another column in the first factor might have a scale of 1016 and the corresponding column in the third factor might have a scale of 1016. One common method of controlling the unbalanced scales is to normalize all columns to unit vectors in the factors [9]. However, by factor normalization, the factor columns will never become zeros vectors. Hence it is impossible to impose sparsity efficiently. The proximal algorithm in our sparse NCP will overcome the above drawbacks. We have mentioned that the proximal will guarantee the full column rank in the model. Moreover, the proximal regularization item in sparse NCP can keep all columns in factors on a balanced scale. Next, we will introduce the solution of the model in (8) using the iHALS method. For the sake of simplification, we use arand brinstead of aðnÞ rand bðnÞ rin this part, which are the rth column of AðnÞand BðnÞ, respectively. We also use AðnÞð:;rÞ¼ar2RIn1to represent the column of a matrix, and AðnÞði;rÞ¼aðnÞ ir to represent an element in a matrix. The objective function in (9) can be further represented as FAðnÞ¼1 2XðnÞX R r¼1 arbT r2 F þan 2X R r¼1jjarearjj2 2þbnX R r¼1jjarjj1; ð12Þ where earis the rth column of e AðnÞ. The minimization problem for (12) can be solved iteratively by columnwise subproblems: Neural Computing and Applications 123 min ar Fr¼1 2ZrarbT r2 Fþan 2jjarearjj2 2þbnjjarjj1 s.t. ar>0; ð13Þ for r¼1;...;R, in which Zr¼XðnÞX R ~ r¼1;~ r6¼r a~ rbT ~ r:ð14Þ The partial derivative of Frwith respect to aris oFr oar¼arbT rZrbrþanaranearþbn1; ¼bT rbrþanarZrbrþanearbn1;ð15Þ where 12RIn1is a vector with all elements equal to 1. When oFr oar¼0, nonnegative column vector arcan be updated as ar Zrbrþanearbn1 bT rbrþan "# þ ;ð16Þ which is a closed form solution of (13) according to the Theorem 2 in [24]. A fast HALS method was utilized to solve the largescale problem [7,24]. We use the same idea to solve the sparse NCP problem. Zrin (14) can also be represented as Zr¼XðnÞX R ~ r¼1 a~rbT ~ rþearbT r:ð17Þ Replacing Zrin (16)by(17), we obtain the new update rule for aras shown in (18). ar XðnÞPR ~r¼1a~ rbT ~ rþearbT rbrþanearbn1 bT rbrþanþ ¼earþXðnÞbrPR ~ r¼1a~rbT ~ rbrbn1 bT rbrþanþ ¼"earþXðnÞBðnÞð:;rÞAðnÞBðnÞTBðnÞð:;rÞbn1 BðnÞTBðnÞðr;rÞþan#þ ð18Þ We implement the above procedures using the inexact scheme. We use PROX-iHALS to denote the inexact hierarchical alternating least squares method for solving the sparse NCP with the proximal algorithm. The PROX- iHALS is illustrated in Algorithm 3. 4.3 Stopping conditions 4.3.1 Stopping condition for outer loop We terminate the outer loop according to the change of relative error during iteration. Relative error is related to data fitting. In the kth outer iteration, the relative error [48] of tensor decomposition is defined by RelErrk¼kXsAð1Þ k;...;AðNÞ ktkF kXkF :ð19Þ Based on the relative error, we terminate the outer loop using the following stopping condition jRelErrk1RelErrkj\e:ð20Þ The threshold of ecan be set by a very small positive value, such as 1e8. In addition, we also set a maximum running time for the outer loop. 4.3.2 Stopping condition for inner loop In the lth inner iteration, we define the relative residual of the nth factor matrix AðnÞas rðnÞ l¼AðnÞ lAðnÞ l1F AðnÞ lF :ð21Þ For the PROX-iHALS, we terminate the inner loop by the stopping condition of rðnÞ l\dðnÞ, where dðnÞis a dynamic positive threshold. If there is only one iteration in the inner loop, we update dðnÞby dðnÞ¼dðnÞ=10. We set the initial value by dðnÞ¼0:01. For the PROX-ANQP, the inner loop is terminated according to the columns in the feasible region of the BPP algorithm [22]. Neural Computing and Applications 123 Since we employ the inexact BCD framework, we also set a maximum number of inner iterations (MAX_- INNER_ITER) to terminate the inner loop. We summarize the stopping conditions for both of the outer and inner loop in Algorithm 4. 4.4 Remarks on convergence The PROX-ANQP and PROX-iHALS methods in the inexact BCD framework have outstanding convergence properties. In Sect. 3.3, we have mentioned that proposed sparse NCP using the proximal algorithm and inexact BCD scheme can guarantee to converge to a stationary point. The subproblem with the proximal algorithm in (9)is strongly convex, which can yield a unique minimum [4]. Furthermore, the optimization methods of ANQP and iHALS can stably decrease the subproblem. According to the Proposition 3.7.1 in [4], both the PROX-ANQP and PROX-iHALS can converge to stationary points. 5 Experiments and results We carried out the experiments on synthetic, real-world, dense, sparse, small-scale, and large-scale tensors. We compared the proposed PROX-ANQP and PROX-iHALS methods with three sparse NCP methods listed below. •AO-ADMM: This is the sparse NCP method using AOADMM algorithm [19], which includes multiple inner iterations. The l1-norm is handled by a proximal operator. •iAPG: We extend the APG method in sparse Tucker decomposition [47] to the sparse NCP problem in (6). In order to make a fair comparison, we implement APG in the inexact scheme using multiple inner iterations, which is abbreviated as iAPG [42]. The l1-norm is handled by a proximal operator. •iMU: This is the sparse NCP method using the classical MU algorithm [9]. We implement MU in the inexact scheme using multiple inner iterations, which is abbreviated as iMU. The above three methods can be directly applied to solve the sparse NCP in (6). The l1-norm can be handled by the proximal operator in AO-ADMM and iAPG. Due to the proximal operator, AO-ADMM and iAPG do not suffer from the rank deficiency. Using the multiplicative updating rule, MU does not suffer from the rank deficiency. In Table 1, we summarized the computational complexity of all the sparse NCP methods. Only the multiplicative operations were counted for mode-nin one outer iteration. The main time cost of these algorithms was spent on the calculation of MTTKRP XðnÞBðnÞ, which consists of two parts: Khatri-Rao product BðnÞand matrix product of XðnÞand BðnÞ. The computational complexity of BðnÞ reaches RQN ~ n¼1;~ n6¼nI~ nand that of XðnÞBðnÞreaches RQN n¼1In. Item BðnÞTBðnÞcan be calculated efficiently by (5), whose complexity is R2PN ~ n¼1;~ n6¼nI~ n. For the inner loop of the subproblem,  Kis assumed to be the average iteration number. In Table 1, we can find that the complexity of these algorithms is highly comparable to each other. It can be inferred that the time of convergence is highly related to the number of iterations. Many experimental parameters and settings will affect the performances of a sparse NCP method. Since our purpose in the experiments is only to test the ability to impose sparsity, we fix the following settings for all methods. •Initialization. For PROX-ANQP, PROX-iHALS, AOADMM and iAPG, all factor matrices were initialized using nonnegative random numbers by MATLAB function max(0,randn(In;R)). Only the iMU was Table 1 Computational Complexity of Multiplicative Operations for Subproblem (9) Method XðnÞBðnÞBðnÞT BðnÞInner loop PROX-ANQP  KðInR2þR3Þ PROX-iHALS  KInR2 AO-ADMM RQN ~n¼1;~n6¼nI~ nR2PN ~n¼1;~n6¼nI~ n KðInR2þR3Þ iAPG þRQN ~ n¼1I~n KInR2 iMU  KInR2  Kis assumed to be the average inner iteration number Neural Computing and Applications 123 initialized by max(0,randn(In;R))?0.1. All initialized factors were scaled by AðnÞ 0¼AðnÞ 0 jjAðnÞ 0jjFffiffiffiffiffiffiffiffiffiffiffiffiffi jjXjjF N p. •The factor updating order was fixed by 1;2;...;N. •The maximum inner iteration MAX_INNER_ITER was fixed by 5 according to the default setting in AOADMM [19]. •For the PROX-ANQP and PROX-iHALS method, the proximal regularization parameter anwas fixed by 1e-4 [45]. The l1-norm regularization parameters of bn;n¼1;...;N; in sparse NCP are the key elements to impose sparsity, which are the most crucial testing parameters in the experiments. We selected a sequence of bnvalues in ascending order for each tensor by manual testing. For synthetic tensors, we stop the increase of bnwhen the true sparse components are recovered, while for real-world tensors, we stop the increase of bnwhen the number of nonzero components is reduced to less than half of the initial number. In order to make it convenient to select and test the parameters, we kept bn;n¼1;...;N;the same in all modes of the tensor. After choosing the bn, we calculated and evaluated the sparsity level [44] of the factor matrices by SparsityAðnÞ¼#AðnÞ i;r\Ts InR;ð22Þ where Tsis a small positive number and # fgdenotes the number of elements that are smaller than the threshold Tsin factor matrix AðnÞ. In the synthetic tensor experiments, we used prior sparse matrices to construct the data. After decomposition, the accuracy of the recovered sparse signals should be evaluated. Let SðnÞ¼½s1;...;sR2RLRdenote the mode-n prior sparse matrix, where Ris the real number of components and Lis the length of a component. Let TðnÞ¼ ½t1;...;t~ R2RL~ Rrepresent the mode-nestimated sparse matrix, where the ~ Ris the estimated number of nonzero components. We evaluate the accuracy of the estimated matrix TðnÞcompared with original sparse signals SðnÞby peak-signal-to-noise ratio (PSNR, see Chapter 3 in [9]) PSNR ¼1 ~ RX ~ R r¼1 10log10 L ^ tr^ sc2 2 ;ð23Þ where ^ tris the rth normalized estimated sparse signal, and ^ scis the normalized reference sparse signal. ^ sccomes from SðnÞ, which has the highest correlation coefficient with ^ tr. All the experiments were conducted on the computer with Intel Core i5-4590 3.30 GHz CPU, 8 GB memory, 64-bit Windows 10 and MATLAB R2016b. The fundamental tensor computation was based on Tensor Toolbox 2.6 [3]. The codes are available on the author’s website http://deqing.me/. 5.1 Synthetic tensor data 5.1.1 Size 1000 ·100 ·100 ·5 with one sparse factor In this experiment, we constructed a synthetic fourth-order tensor by 10 channels of simulated sparse and nonnegative signals, as shown in Fig. 1a. The signals come from the file VSparse_rand_10.mat in NMFLAB [8]. There are 1000 points in each channel, so the sparse signal matrix is 1001 200 300 400 500 600 700 800 900 1000 Point(n) 9 3 8 6 4 7 5 2 1 (b) Estimated Sparse Signals 1001 200 300 400 500 600 700 800 900 1000 Point(n) 1 2 3 4 5 6 7 8 9 10 Magnitude of All Channels (a) Original Sparse Signals 10 Fig. 1 Sparse and nonnegative signals used in synthetic tensor. a shows the original ten channels of signals. bshows the estimated ten channels of signals from the synthetic tensor XSYN1 by sparse NCP based on the PROX-ANQP method with bn¼5. The PSNR is 90.2698 according to (23) Neural Computing and Applications 123 NCP methods are also tested for comparison, including the AO-ADMM, iAPG and iMU. We have the following findings: (1) The PROX-ANQP and PROX-iHALS methods particularly have the fast convergence speed and excellent effects of imposing sparsity in all cases compared with other methods. The outstanding performances of the PROX-ANQP and PROX-iHALS are due to two points: the proximal algorithm that overcomes the rank deficiency; and the inexact scheme that increases the efficiency. (2) The iAPG method contains a proximal operator, which can handle the l1-norm. With the proximal operator, iAPG does not suffer from the rank deficiency. In the experiments, we find that iAPG is very efficient to impose sparsity for small-scale dense tensor decomposition. Nevertheless, iAPG runs slowly for large-scale sparse tensor decomposition. (3) The iMU method converges very slowly for dense tensors, but it becomes fast for large-scale sparse tensor. For sparse tensor decomposition, the extracted factor matrices are extremely sparse already. Most elements in the factors are zeros. According to the multiplicative updating rule, once an element becomes zero, it will never change. This property might be the reason why iMU converges fast on the sparse tensor. (4) The AOADMM also contains a proximal operator that can handle l1-norm and overcome the rank deficiency. However, AOADMM converges slowly compared with PROX-ANQP and PROX-iHALS in most cases, particularly in the largescale sparse tensor case. Moreover, AO-ADMM is inferior to iAPG with bn¼0 in many cases. In a word, our proposed PROX-ANQP and PROX-iHALS methods have the best performances for sparse NCP, and have very good generalization for the different types and scales of datasets. In addition to the solving methods, another critical issue of sparse NCP is the selection of the sparse regularization parameter bn. Firstly, we want to mention that one purpose of this paper is to demonstrate the effectiveness of the algorithms to impose sparsity. Therefore, in order to simplify the selection of parameters, we keep bnthe same for all factor matrices in sparse NCP. With the same bnon all modes, the sparse NCP can still recover high sparse components and low sparse (or even dense) components in different factor matrices. In the future, it would be interesting to investigate how to separately control the sparsity levels of different factor matrices using unbalanced sparse regularization parameters. Secondly, the appropriate value of parameter bndepends on the tensor to be decomposed. In this study, we selected bnseparately for each tensor in the experiments. When the sparse regularization parameter is larger, the extracted factor matrices are sparser, and the relative error of decomposition is also larger. The trade-off Table 6 Comparison of Sparse NCPs on NIPS Publications Tensor XSp2 2R248228621403617 Method bnObj RelErr Time(s) Iter NNC Spars1Spars2Spars3Spars4 PROX- 0 3.11e?07 0.9439 1091.8 18.5 100.00 0.999 0.998 0.920 0.941 ANQP 0.1 3.12e?07 0.9447 1161.9 19.9 99.57 0.999 0.998 0.925 0.941 0.3 3.20e?07 0.9571 990.5 17.0 78.43 0.999 0.998 0.950 0.954 0.5 3.38e?07 0.9841 749.3 12.9 24.83 10.999 0.985 0.985 PROX- 0 3.11e?07 0.9439 991.4 17.2 100.00 0.999 0.997 0.920 0.941 iHALS 0.1 3.12e?07 0.9448 1123.2 19.3 99.93 0.999 0.998 0.925 0.941 0.3 3.20e?07 0.9577 1055.5 18.0 79.60 0.999 0.998 0.948 0.953 0.5 3.37e?07 0.9827 804.5 13.9 28.50 10.999 0.982 0.983 AO- 0 3.11e?07 0.9440 1346.7 23.0 100.00 0.999 0.997 0.919 0.938 ADMM 0.1 3.11e?07 0.9446 1191.1 20.6 99.70 0.999 0.998 0.924 0.939 0.3 3.23e?07 0.9611 1358.5 23.4 74.70 0.998 0.998 0.951 0.950 0.5 3.40e?07 0.9865 912.5 15.5 21.93 0.999 0.999 0.986 0.986 iAPG 0 3.11e?07 0.9436 1263.9 21.9 100.00 0.999 0.997 0.920 0.941 0.1 3.12e?07 0.9449 1325.0 22.8 99.47 0.999 0.998 0.925 0.940 0.3 3.22e?07 0.9600 1587.2 27.1 73.37 0.999 0.998 0.951 0.949 0.5 3.38e?07 0.9837 1131.5 19.5 27.23 10.999 0.982 0.979 iMU 0 3.12e?07 0.9448 1126.3 19.4 100.00 0.999 0.998 0.921 0.941 0.1 3.12e?07 0.9456 1197.4 20.5 99.67 0.999 0.998 0.927 0.941 0.3 3.28e?07 0.9685 1317.3 22.5 60.07 0.999 0.999 0.964 0.965 0.5 3.43e?07 0.9916 920.7 15.8 12.73 110.992 0.992 Sparsn= Sparsity level of the mode-nestimated factor Sparsn1 means that the factor is very close to a zero matrix NNC = Number of nonzero components Neural Computing and Applications 123 between the sparse level and the relative error depends on the meanings of real applications. An example of sparse regularization parameter selection of sparse NCP for ongoing EEG can be found in [44]. It is also possible to select an appropriate parameter for a concrete application using model-order selection methods [37], such as the Bayesian information criteria (BIC). The third critical issue of sparse NCP is the sparse regularization item. In this study, we only investigated the l1-norm item. In the future, it is worth trying to incorporate other types of sparse regularization items [2] to our sparse NCP model besides l1-norm, such as the lq-norm (0\q\1) [36] and trace norm [29]. 7 Conclusion In this paper, we have investigated the nonnegative CANDECOMP/PARAFAC tensor decomposition with l1- norm-based sparse regularization (sparse NCP). We have proposed a novel sparse NCP model using the proximal algorithm, which can guarantee the full column rank condition and the property of convergence to a stationary point. In addition, an inexact block coordinate descent scheme was presented to accelerate the computation of sparse NCP. In the inexact scheme, the subproblems are updated using multiple inner iterations. We have employed two algorithms for solving the proposed sparse NCP model with the proximal algorithm, including the inexact alternating nonnegative quadratic programming (PROXANQP) and the inexact hierarchical alternating least squares (PROX-iHALS). The experimental results on all synthetic, real-world, small-scale and large-scale tensors demonstrated that our sparse NCP methods can impose sparsity and extract meaningful sparse components successfully. Both PROX-ANQP and PROX-iHALS have exhibited the faster computational speed and better performances of imposing sparsity compared with other sparse NCP algorithms. The experimental results proved that the proposed sparse NCP with the proximal algorithm and inexact scheme is effective and efficient. Acknowledgements This work was supported by National Natural Science Foundation of China (Grant No.91748105), National Foundation in China (No. JCKY2019110B009 & 2020-JCJQ-JJ-252), the Objective Function Value Time(s) Time(s) Time(s) Time(s) n= 0 n= 0.1 n= 0.3 n= 0.5 107107 107107 3.1 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 3.1 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 3.1 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 3.1 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 0 200 400 600 800 1000 1200 iMU PROX-ANQP PROX-iHALS iAPG AO-ADMM iMU PROX-ANQP PROX-iHALS iAPG AO-ADMM iMU PROX-ANQP PROX-iHALS iAPG AO-ADMM iMU PROX-ANQP PROX-iHALS iAPG AO-ADMM 0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200 Fig. 4 The objective function value curves of sparse NCPs on fourth-order NIPS publications tensor Neural Computing and Applications 123 Fundamental Research Funds for the Central Universities [DUT20LAB303 & DUT20LAB308] in Dalian University of Technology in China, and the scholarship from China Scholarship Council (No. 201600090043). This study is to memorize Prof. Tapani Ristaniemi for his great help to Fengyu Cong, Zheng Chang and Deqing Wang. Prof. Tapani Ristaniemi supervised this work. Funding Open access funding provided by University of Jyva ¨skyla ¨ (JYU). Declarations Conflict of interest The authors declare that they have no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/. References 1. Allen G (2012) Sparse higher-order principal components analysis. In: Lawrence ND, Girolami M (eds) Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics, PMLR, La Palma, Canary Islands, Proceedings of Machine Learning Research, 22: 27–36 2. Bach F, Jenatton R, Mairal J, Obozinski G (2012) Optimization with sparsity-inducing penalties. Found TrendsMach Learn 4(1):1–106. https://doi.org/10.1561/2200000015 3. Bader BW, Kolda TG et al (2015) Matlab tensor toolbox version 2.6. Available online, http://www.sandia.gov/*tgkolda/TensorToolbox/ 4. Bertsekas DP (2016) Nonlinear Programming, 3rd edn. Athena Scientific, Belmont, Massachusetts 5. Bro R, Kiers HAL (2003) A new efficient method for determining the number of components in PARAFAC models. J Chemom 17(5):274–286. https://doi.org/10.1002/cem.801 6. Bruckstein AM, Donoho DL, Elad M (2009) From sparse solutions of systems of equations to sparse modeling of signals and images. SIAM Rev 51(1):34–81. https://doi.org/10.1137/ 060657704 7. Cichocki A, Phan AH (2009) Fast local algorithms for large scale nonnegative matrix and tensor factorizations. IEICE Trans Fundam Electron Commun Comput Sci E92–A(3):708–721. https:// doi.org/10.1587/transfun.e92.a.708 8. Cichocki A, Zdunek R (2006) NMFLAB—MATLAB toolbox for non-negative matrix factorization. http://www.bsp.brain.riken.jp/ ICALAB/nmflab.html. Accessed 22 Nov 2017 9. Cichocki A, Zdunek R, Phan AH, Si Amari (2009) Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. Wiley, New York 10. Cichocki A, Mandic D, Lathauwer LD, Zhou G, Zhao Q, Caiafa C, Phan HA (2015) Tensor decompositions for signal processing applications: from two-way to multiway component analysis. IEEE Signal Process Mag 32(2):145–163. https://doi.org/10. 1109/msp.2013.2297439 11. Cong F, Lin QH, Kuang LD, Gong XF, Astikainen P, Ristaniemi T (2015) Tensor decomposition of EEG signals: a brief review. J Neurosci Methods 248:59–69. https://doi.org/10.1016/j.jneu meth.2015.03.018 12. Donoho DL (2006) For most large underdetermined systems of linear equations the minimal ‘1-norm solution is also the sparsest solution. Commun Pure Appl Math 59(6):797–829. https://doi. org/10.1002/cpa.20132 13. Elcoroaristizabal S, Bro R, Garcı ´a JA, Alonso L (2015) PARAFAC models of fluorescence data with scattering: a comparative study. Chemom Intell LabD Syst 142:124–130. https://doi.org/10. 1016/j.chemolab.2015.01.017 14. Friedlander MP, Hatz K (2008) Computing non-negative tensor factorizations. Optim Methods Softw 23(4):631–647. https://doi. org/10.1080/10556780801996244 15. Gillis N, Glineur F (2012) Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization. Neural Computation 24(4):1085–1105. https://doi.org/10.1162/ NECO_a_00256 16. Globerson A, Chechik G, Pereira F, Tishby N (2007) Euclidean Embedding of Co-occurrence Data. The Journal of Machine Learning Research 8: 2265–2295, http://www.jmlr.org/papers/ volume8/globerson07a/globerson07a.pdf 17. Hong M, Razaviyayn M, Luo ZQ, Pang JS (2016) A unified algorithmic framework for block-structured optimization involving big data: with applications in machine learning and signal processing. IEEE Signal Process Mag 33(1):57–77. https:// doi.org/10.1109/msp.2015.2481563 18. Hoyer PO (2004) Non-negative matrix factorization with sparseness constraints. J Mach Learn Res 5(Nov):1457–1469 19. Huang K, Sidiropoulos ND, Liavas AP (2016) A flexible and efficient algorithmic framework for constrained matrix and tensor factorization. IEEE Trans Signal Process 64(19):5052–5065. https://doi.org/10.1109/tsp.2016.2576427 20. Kim H, Park H (2008) Nonnegative matrix factorization based on alternating nonnegativity constrained least squares and active set method. SIAM J Matrix Anal Appl 30(2):713–730. https://doi. org/10.1137/07069239x 21. Kim HJ, Ollila E, Koivunen V (2013) Sparse regularization of tensor decompositions. In: 2013 IEEE international conference on acoustics, speech and signal processing. IEEE. https://doi.org/ 10.1109/ICASSP.2013.6638376 22. Kim J, Park H (2011) Fast nonnegative matrix factorization: an active-set-like method and comparisons. SIAM J Sci Comput 33(6):3261–3281. https://doi.org/10.1137/110821172 23. Kim J, Park H (2012) Fast nonnegative tensor factorization with an active-set-like method. In: High-performance scientific computing. Springer, London, pp 311–326. https://doi.org/10.1007/ 978-1-4471-2437-5_16 24. Kim J, He Y, Park H (2014) Algorithms for nonnegative matrix and tensor factorizations: a unified view based on block coordinate descent framework. J Global Optim 58(2):285–319. https:// doi.org/10.1007/s10898-013-0035-4 25. Kolda TG, Bader BW (2009) Tensor decompositions and applications. SIAM Rev 51(3):455–500. https://doi.org/10.1137/ 07070111x 26. Li N, Kindermann S, Navasca C (2013) Some convergence results on the regularized alternating least-squares method for tensor decomposition. Linear Algebra Appl 438(2):796–812. https://doi.org/10.1016/j.laa.2011.12.002 27. Li Y, Ngom A (2013) The non-negative matrix factorization toolbox for biological data mining. Source Code Biol Med 8(1):10. https://doi.org/10.1186/1751-0473-8-10 Neural Computing and Applications 123 28. Liu J, Liu J, Wonka P, Ye J (2012) Sparse non-negative tensor factorization using columnwise coordinate descent. Pattern Recogn 45(1):649–656. https://doi.org/10.1016/j.patcog.2011.05.015 29. Liu Y, Shang F, Jiao L, Cheng J, Cheng H (2015) Trace norm regularized CANDECOMP/PARAFAC decomposition with missing data. IEEE Trans Cybern 45(11):2437–2448. https://doi. org/10.1109/tcyb.2014.2374695 30. Mørup M (2011) Applications of tensor (multiway array) factorizations and decompositions in data mining. Wiley Interdiscip Rev Data Mining Knowl Discov 1(1):24–40. https://doi.org/10. 1002/widm.1 31. Mørup M, Hansen LK, Arnfred SM (2008) Algorithms for sparse nonnegative tucker decompositions. Neural Comput 20(8):2112–2131. https://doi.org/10.1162/neco.2008.11-06-407 32. Papalexakis EE, Sidiropoulos ND, Bro R (2013) From k-means to higher-way co-clustering: Multilinear decomposition with sparse latent factors. IEEE Trans Signal Process 61(2):493–506. https:// doi.org/10.1109/tsp.2012.2225052 33. Razaviyayn M, Hong M, Luo ZQ (2013) A unified convergence analysis of block successive minimization methods for nonsmooth optimization. SIAM J Optim 23(2):1126–1153. https://doi. org/10.1137/120891009 34. Selesnick I (2017) Sparse regularization via convex analysis. IEEE Trans Signal Process 65(17):4481–4494. https://doi.org/10. 1109/tsp.2017.2711501 35. Sidiropoulos ND, Lathauwer LD, Fu X, Huang K, Papalexakis EE, Faloutsos C (2017) Tensor decomposition for signal processing and machine learning. IEEE Trans Signal Process 65(13):3551–3582. https://doi.org/10.1109/tsp.2017.2690524 36. Sigurdsson J, Ulfarsson MO, Sveinsson JR (2014) Hyperspectral unmixing with lqregularization. IEEE Trans Geosci Remote Sens 52(11):6793–6806. https://doi.org/10.1109/tgrs.2014.2303155 37. Stoica P, Selen Y (2004) Model-order selection: a review of information criterion rules. IEEE Signal Process Mag 21(4):36–47. https://doi.org/10.1109/msp.2004.1311138 38. Timmerman ME, Kiers HAL (2000) Three-mode principal components analysis: choosing the numbers of components and sensitivity to local optima. British J Math Stat Psychol 53(1):1–16. https://doi.org/10.1348/000711000159132 39. Veganzones MA, Cohen JE, Farias RC, Chanussot J, Comon P (2016) Nonnegative tensor CP decomposition of hyperspectral data. IEEE Trans Geosci Remote Sens 54(5):2577–2588. https:// doi.org/10.1109/tgrs.2015.2503737 40. Vervliet N, Lathauwer LD (2019) Numerical optimization-based algorithms for data fusion. In: Data handling in science and technology. Elsevier, pp 81–128. https://doi.org/10.1016/B978-0- 444-63984-4.00004-1 41. Viswanath B, Mislove A, Cha M, Gummadi KP (2009) On the evolution of user interaction in facebook. In: Proceedings of the 2nd ACM workshop on Online social networks—WOSN ’09. ACM Press. https://doi.org/10.1145/1592665.1592675 42. Wang D, Cong F (2021) An inexact alternating proximal gradient algorithm for nonnegative CP tensor decomposition. Sci China Technol Sci 64(9):1893–1906. https://doi.org/10.1007/s11431- 020-1840-4 43. Wang D, Cong F, Zhao Q, Toiviainen P, Nandi AK, Huotilainen M, Ristaniemi T, Cichocki A (2016) Exploiting ongoing EEG with multilinear partial least squares during free-listening to music. In: IEEE 26th international workshop on machine learning for signal processing (MLSP). IEEE. https://doi.org/10.1109/ mlsp.2016.7738849 44. Wang D, Wang X, Zhu Y, Toiviainen P, Huotilainen M, Ristaniemi T, Cong F (2018) Increasing stability of EEG components extraction using sparsity regularized tensor decomposition. In: Advances in neural networks – ISNN 2018. Springer, pp 789–799. https://doi.org/10.1007/978-3-319-92537-0_89 45. Wang D, Cong F, Ristaniemi T (2019) Higher-order nonnegative CANDECOMP/PARAFAC tensor decomposition using proximal algorithm. In: 2019 IEEE international conference on acoustics, speech and signal processing (ICASSP). IEEE, pp 3457–3461. https://doi.org/10.1109/ICASSP.2019.8683217 46. Williams AH, Kim TH, Wang F, Vyas S, Ryu SI, Shenoy KV, Schnitzer M, Kolda TG, Ganguli S (2018) Unsupervised discovery of demixed, low-dimensional neural dynamics across multiple timescales through tensor component analysis. Neuron 98(6):1099-1115.e8. https://doi.org/10.1016/j.neuron.2018.05. 015 47. Xu Y (2015) Alternating proximal gradient method for sparse nonnegative tucker decomposition. Math Program Comput 7(1):39–70. https://doi.org/10.1007/s12532-014-0074-y 48. Xu Y, Yin W (2013) A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion. SIAM J Imag Sci 6(3):1758–1789. https://doi.org/10.1137/120887795 49. Yang Y, Pesavento M, Luo ZQ, Ottersten B (2020) Inexact block coordinate descent algorithms for nonsmooth nonconvex optimization. IEEE Trans Signal Process 68:947–961. https://doi.org/ 10.1109/tsp.2019.2959240 50. Zhang H, Wang S, Xu X, Chow TWS, Wu QMJ (2018) Tree2vector: learning a vectorial representation for tree-struc- tured data. IEEE Trans Neural Networks Learn Syst 29(11):5304–5318. https://doi.org/10.1109/tnnls.2018.2797060 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Authors and Affiliations Deqing Wang 1,2 •Zheng Chang 2,3 •Fengyu Cong 1,2,4,5 &Deqing Wang [email protected] &Fengyu Cong [email protected] Zheng Chang [email protected] 1 School of Biomedical Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian 116024, China 2 Faculty of Information Technology, University of Jyva ¨skyla ¨, Jyva ¨skyla ¨40100, Finland Neural Computing and Applications 123 3 School of Computer Science and Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China 4 School of Artificial Intelligence, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian 116024, China 5 Key Laboratory of Integrated Circuit and Biomedical Electronic System, Liaoning Province, Dalian University of Technology, Dalian 116024, China Neural Computing and Applications 123