Davies, Michael

Yaghoobi Vaighan, Mehrdad

An important class of sparse signals is the non-negative sparse signals. While numerous greedy techniques have been introduced for low-complexity sparse approximations, there are few non-negative versions. Among such a large class of greedy techniques, one successful method, which is called the Orthogonal Least Squares (OLS) algorithm, is based on the maximum residual energy reduction at each iteration. However, the basic implementation of the OLS is computationally slow. The OLS algorithm has a fast implementation based on the QR matrix factorisation of the dictionary. The extension of such technique to the non-negative domain is possible. In this paper, we present a fast implementation of the non-negative OLS (NNOLS). The computational complexity of the algorithm is compared with the basic implementation, where the new method is faster with two orders of magnitude. We also show that, if the basic implementation of NNOLS is not computationally feasible for moderate size problems, the proposed method is tractable. We also show that the proposed algorithm is even faster than an approximate implementation of the non-negative OLS algorithm.</p

Edinburgh Research Explorer

     Edinburgh Research Explorer                                      Fast non-negative orthogonal least squaresCitation for published version:Yaghoobi Vaighan, M & Davies, M 2015, Fast non-negative orthogonal least squares. in 2015 23rdEuropean Signal Processing Conference, EUSIPCO 2015., 7362429, Institute of Electrical and ElectronicsEngineers Inc., pp. 479-483, 23rd European Signal Processing Conference, EUSIPCO 2015, Nice, France,31/08/15. DOI: 10.1109/EUSIPCO.2015.7362429Digital Object Identifier (DOI):10.1109/EUSIPCO.2015.7362429Link:Link to publication record in Edinburgh Research ExplorerDocument Version:Peer reviewed versionPublished In:2015 23rd European Signal Processing Conference, EUSIPCO 2015General rightsCopyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)and / or other copyright owners and it is a condition of accessing these publications that users recognise andabide by the legal requirements associated with these rights.Take down policyThe University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorercontent complies with UK legislation. If you believe that the public display of this file breaches copyright pleasecontact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately andinvestigate your claim.Download date: 05. Apr. 2019FAST NON-NEGATIVE ORTHOGONAL LEAST SQUARESMehrdad Yaghoobi and Mike E. DaviesInstitute for Digital Communications, the University of Edinburgh, EH9 3JL, UK{m.yaghoobi-vaighan, d.wu and mike.davies}@ed.ac.ukABSTRACTAn important class of sparse signals is the non-negativesparse signals. While numerous greedy techniques have beenintroduced for low-complexity sparse approximations, thereare few non-negative versions. Among such a large class ofgreedy techniques, one successful method, which is called theOrthogonal Least Squares (OLS) algorithm, is based on themaximum residual energy reduction at each iteration. How-ever, the basic implementation of the OLS is computationallyslow. The OLS algorithm has a fast implementation basedon the QR matrix factorisation of the dictionary. The exten-sion of such technique to the non-negative domain is possi-ble. In this paper, we present a fast implementation of thenon-negative OLS (NNOLS). The computational complexityof the algorithm is compared with the basic implementation,where the new method is faster with two orders of magnitude.We also show that, if the basic implementation of NNOLS isnot computationally feasible for moderate size problems, theproposed method is tractable. We also show that the proposedalgorithm is even faster than an approximate implementationof the non-negative OLS algorithm.Index Terms— Non-negative sparse approximations,Orthogonal Least Squares, Efficient Implementations, Non-negative Orthogonal Least Squares and QR Matrix Factoriza-tion1. INTRODUCTIONThe signal processing with a low-dimensional model, is amodern technique to explore signal structures. We thus needsome low-dimensional signal approximation algorithms forthis purpose. Such methods have to be fast, if the dimen-sion of data is large/huge. There are many general purposesparse approximation techniques with fast implementations.In some applications like hyperspectral data decomposition[1] and DNA microarrays analysis [2], we further know thatnon-zero coefficients of the sparse approximation are posi-tive. In this setting, non-negative sparse approximation tech-niques [3], often provide better solutions than the universaltechniques. Universal greedy techniques can be adapted tosuch a setting by confining the solutions to the set of non-negative coefficients. However, extension of the fast imple-mentations of the greedy methods to the non-negative settingis often not straightforward. An example of such extensionsto a greedy method, with a complicated coefficient update e.g.Orthogonal Matching Pursuit (OMP), has already been pre-sented in [4]. We here present a fast implementation of thenon-negative OLS [5] technique.Let y ∈ Rm be a signal with the sparse representationusing a normalised linear generative model Φ ∈ Rm×n, i.e.dictionary, where y = Φx and x is a sparse vector. At each it-eration, OLS selects an element of dictionary, called an atom,which maximally reduces the residual energy of the signal,after adding the atom to the set of previously selected atoms.It then finds the best possible representation of y using se-lected atoms, i.e. orthogonal projection. In the other words,the algorithm has these steps: i) selecting an atom, i.e.i∗ = argmini minsk+1:=sk∪{i}‖y −Φsk+1Φ†sk+1y‖2, (1)where sk is the index set of selected atoms up to the iterationk, Φsk+1 is the subdictionary selected by the index set sk and† represents the More-Penrose pseudoinverse and ii) using theLeast Squares (LS) representation within the selected atoms.As the optimisation problem (1) is computationally expen-sive, Chen et al. [5] presented a technique, which iterativelyupdates the QR factorisation of the selected subdictionary, i.e.similar to the fast QR factorisation based OMP implementa-tion [6]. It has an extra step in which we orthogonalise the restof the dictionary with respect to the selected subdictionary,and then normalise the new dictionary. We refer the readersto [5], for a more detailed description of the OLS algorithm.If y is “non-negative sparse”, i.e. the representation x hasthe property ‖x‖0 ≤ k, x ∈ Rn+, we intend to incorporatethe extra information and adapt OLS to the new setting. Thisadaptation can be done if we combine the two steps of theOLS and enforce the non-negativity of the coefficients, i.e.i∗ = argmini minsk+1 := sk ∪ {i}α ∈ Rk+1+‖y −Φsk+1α‖2. (2)We call such an implementation Basic NNOLS (BNNOLS).The extension of the fast implementation of OLS for the non-negative representation is not straightforward, if we want toexactly solve NNOLS.Algorithm 1 Suboptimal Non-Negative OLS1: initialisation: s = ∅, k = 0, x = 0, Ω = Φ and r0 = y2: while k < K & max(ΩT rk) > 0 do3: ι∗ ← max(ΩT rk)4: s← s ∪ ι∗5: xs ← argminθ>0‖y −Φsθ‖26: rk+1 ← y −Φsxs7: Updating Ω = [ωi], ωi, ∀i ∈ sc, from (3)8: k ← k + 19: end while10: x|s ← xsWe here present a simple technique to extend the fast OLSimplementation in section 2, without solving an expensive op-timisation problem in the loop. We then show that the newimplementation is significantly faster than a basic implemen-tation of NNOLS in section 3.2. NON-NEGATIVE OLSLet Φk ∈ Rm×k be the k selected atoms of Φ at iteration k, inwhich the columns are ordered based on the iteration number,i.e. the jth column is the jth selected atom. We can nowgenerate Φ¯k ∈ Rm×(n−k) as the rest of the dictionary, i.e.including all non-selected atoms. In each iteration of OLS,we need to calculate a matrix Ω = [ωi]1≤i≤n−k, which is thecolumn-normalised version of the orthogonal part of Φ¯k tothe column span of Φk. In other words,ωi :=φ¯⊥i‖φ¯⊥i ‖2, (3)where φ¯⊥i := (I − ΦkΦ†k)φ¯i. We recall that φ¯i’s with theproperty φ¯⊥i = 0, will be eliminated, as they lie in the spanof Φk and there is no need to consider them in the follow-ing iterations. Let rk = y − Φkxk be the residual error ofthe representation at iteration k, and xk is the coefficient vec-tor at this iteration. As Ω is orthogonal to the span of Φkand normalised, we can find the direction with the maximumresidual energy reduction by simply finding the maximum of|ΩTrk|, where | · | is the element by element absolute valueoperator. The optimal direction corresponding to ωi∗ pro-vides us the index of an atom in this iteration. If we nowconsider the non-negative representation, one suboptimal ap-proach is to: i) select an atom with positive contribution, i.emaxi(ωTi rk) > 0, ii) updating the coefficients using NNLS.A pseudo code for this algorithm, which is called SuboptimalNNOLS (SNNOLS) is presented in Algorithm 1. sc is thecomplement set of s. SNNOLS method does not find the so-lution of iterating (2). The reason is related to the fact that theselection criteria does not consider the fact that representationbased on the selected atoms in Φk have to be positive. In thecase which we enforce an already selected atom to have a zerocoefficient value, we may have some suboptimality. We alsoAlgorithm 2 Fast Non-Negative Orthogonal Least Squares1: initialisation: s = z0 = ∅, k = 0, Ω = Φ and r0 = y2: while k < K & max(ΩT rk) > 0 do3: (ζ, ι)← sort↓(ΩT rk)4: p← 15: pc ← 16: zc ← 07: while ∼ Terminate & p < N do8: zt from (6)9: z ← ψTι(p)rk : ψι(p) = Φ¯⊥k /‖Φ¯⊥k ‖2,Φ¯⊥k = (I−ΨkΨTk )φ¯ι(p)10: Update based on Table 111: end while12: s← s ∪ ι(p)13: Update Ψ, Φ¯, Ω and R−114: zk+1 ← [zk; zk+1]15: rk+1 ← rk − zk+1ψk+116: k ← k + 117: end while18: output: x|s ← R−1zKIf Then0 < z ≤ zt, z > zc zk+1 ← z,Terminate0 < z ≤ zt, z ≤ zc zk+1 ← zc, p← pc,Terminatez > zc ≥ zt p← p+ 1zc ≥ z > zt zk+1 ← zc, p← pc,Terminatez > zt > zc zc ← zt, pc ← p, p← p+ 1z < 0 TerminateTable 1. Decision rules to guarantee positivity of the coeffi-cients.need to consider the computational cost of SNNOLS, as ititeratively needs to solve NNLS. The solution which we pro-pose here is to combine the selection and coefficient updatesteps to avoid such suboptimality and reduce the computa-tional cost by not directly solving NNLS.2.1. Fast Non-Negative OLSIf Ψk ∈ Rm×k is iteratively generated using the selectedatoms of Ω, up to the kth iteration, we can QR factorise Φkby ΨkRk, where Ψk is orthonormal and Rk ∈ Rk×k is anupper-diagonal matrix. The reason for such a structure for Ris the fact that Φj is orthogonal to ψp’s, for ∀p : j < p.This QR factorisation is updating throughout the iterations asfollows: i) adding ωi∗ as the new column of Ψ, i.e.Ψk+1 = [Ψk ωi∗ ] (4)and ii) updating Rk using the following equation,Rk+1 =[Rk ν0 µ], (5)0 5 10 15 20 25 30 35 40 45 5000.20.40.60.81Average Exact RecoveryBNNOLS, M = 64FNNOLS, M = 64BNNOLS, M = 128FNNOLS, M = 128Sparsity0 5 10 15 20 25 30 35 40 45 5010 -210 010 2Computational Time (Sec)BNNOLS, M = 64FNNOLS, M = 64BNNOLS, M = 128FNNOLS, M = 128Fig. 1. Exact recovery and computation time. N = 256, Mis 64 or 128 and sparsity K is changing. The proposed tech-nique FNNOLS is roughly 100 times faster than BNNOLS!where ν = ΨTk φ¯i∗ , µ = ‖φ¯⊥i∗‖2 and φ¯⊥i∗ can be easily calcu-lated using φ¯⊥i∗ = (I −ΨkΨTk )φ¯i∗ , i.e. no need to calculatethe pseudoinverse of Φk, comparing to the definition after (3).This process can be interpreted as a type of Graham-Schmidt(GS) orthogonalisation process, which is a well-known QRfactorisation technique [7]. With the new setting, we canrewrite the OLS algorithm as follows: In the first iteration,we do not need any orthogonalisation and we have φ1 = ψ1and R1 = [ 1 ]. In the kth iteration, k ≥ 1, let the best ap-proximation of y, using Φk, be∑ki=1 xiφi =∑ki=1 ziψi. Inthe k + 1 iteration, we have,k+1∑i=1ziψi =k∑i=1ziψi + zk+1ψk+1=k∑i=1xiφi + zk+1ψk+1.We recall that ψi’s are orthogonal and choosing a new ele-mentary function ψk+1 does not change the optimal value ofthe best representation zi, ∀i : i < k + 1, calculated up tothe previous iteration k. This is obviously not true for xi’s asΦk+1 is not orthogonal.As ψk+1 lives in the span of the non-redundant set{φj}j∈[1,k+1], ψk+1 =∑k+1j=1 γjφj for some unique γj’s.Fig. 2. Computation time for the fixed N = 256 and K =24 and 32.We can then have,k+1∑i=1ziψi =k∑i=1xiφi +k+1∑j=1zk+1γjφj=k∑i=1(xi + zk+1γi)φi + zk+1γk+1φk+1.As zk+1γk+1 is always positive, we only need to assure thatxi + zk+1γi ≥ 0 orzk+1 ≤ zt :={mini,γi<0∣∣∣xiγi ∣∣∣ ∃i, γi < 0+∞ Otherwise .(6)In OLS, we only need to keep zi’s in the memory andcalculate xi’s at the end of algorithm [8]. To assure thatxi’s are all non-negative, zi’s should comply the conditionof (6). We then choose the atom that the corresponding zk+1,or shrunk by the upper-bound of (6), has the largest value.We therefore need to track the record of the best possiblevalue for z, if z = ωTi∗rk does not comply with (6). If(ζ, ι) = sort↓(ΩTrk), where sort↓ is the sorting operator in adescent order, ζ is the sorted coefficients with correspondingindices ι. A pseudocode of the fast non-negative orthogonalleast square is presented in Algorithm 2. In this algorithm, zis the current candidate, starting with z = ζ(p), p = 1, in aninternal loop in the kth iteration, we make the decision basedon the rules of Table 1, to guarantee that the coefficients xi’swill remain non-negative throughout iterations. We also needto calculate R−1 at the last iteration, which would be expen-sive. We thus can iteratively generate R−1 using the followingupdate rule,R−1k+1 =[R−1k −R−1k νµ0 1µ]. (7)Fig. 3. Computation time for the fixed M = 128 and K =64 and 96.where ν and µ were defined after (5). For the fast implemen-tation of this algorithm, we have to also efficiently calculateγ. It is not difficult to show that γ is the last column of R−1k+1,i.e., γ = [−R−1k νµ ;1µ ], which is kept in the memory.3. SIMULATIONSWe investigate the runtime of the algorithms on a single coreof an Intel core 2.6 GHz processor. In the first experiment,we randomly generated Φ with 256 atoms and 64 or 128 rowsusing i.i.d. Gaussian random variables, followed by columnnormalisation. y was generated using a Gaussian-Bernoullimodel, i.e. uniformly random support and Normal distribu-tion for the non-zero coefficients. We repeated the simula-tions 100 times, using proposed FNNOLS and BNNOLS anddifferent sparsity levels. The correct recovery of the supportand computational time are shown in Figure 1. While theexact recoveries are very similar, proposed method is signif-icantly faster. The difference is more obvious for large k’s,where we can get a computation gain of 100.In the second experiment, we used the same method togenerate the dictionary, but we fixed N = 256, to investigatethe computation time as a function of M . Here, K = 24 or32 and M is between 32 and 196. The computational time isplotted in Figure 2. The computational time is slowly increas-ing with M . We can observe that the difference in computa-tions is not significant for different K’s using the fast tech-nique, while it is noticeable using the basic implementation.We then fixed M = 128 and changed N from 128 to512 and similarly repeated each simulation 100 times. Thesparsity was set to be 64 or 96, and the results are shown inFigure 3. The computational cost is slowly increasing w.r.t.the parameter N in FNNOLS, while it is increasing faster forthe BNNOLS. However, the difference between the compu-Fig. 4. Exact recovery and computation time. Comparison be-tween SNNOLS and FNNOLS. N = 1024, M = 512 spar-sity K is changing.tational time of the fast and basic implementations are con-siderable, where we yield a gain of 1000 times acceleration,using the FNNOLS algorithm!We also presented an approximate NNOLS technique,called the SNNOLS, which can be seen as an intuitive exten-sion of the fast OLS implementation to the non-negative do-main. As we can see, the computational time of SNNOLS isalso low, which means that we can compare it with FNNOLSin a larger dimension setting. We then chose M = 512,N = 1024 and changed the sparsity level K from 32 to 80and averaged the simulations, with the random settings whichexplained before, over 500 trials. We have plotted the ex-act support recovery (top panel) and the algorithms run-time(bottom panel) in Figure 4. While the exact recoveries arevery similar, the computational times are slightly different.For the example, the computational time of SNNOLS is 25%higher than FNNOLS for K = 80. The similarity betweenthe exact recovery of two algorithms is an interesting fact,which needs more investigation.4. CONCLUSIONA significantly faster implementation for the non-negativeorthogonal least square algorithm was proposed here. Thenew implementation is based on the recursive generation ofa QR factorisation of the dictionary and it allows us to useFNNOLS in a broader application. We also presented analgorithm which approximately solve the OLS problem in anon-negative setting, which is called SNNOLS. This algo-rithm performs similarly in the exact recovery of the sparsitysupport, when the dictionaries are randomly generated, butit is slightly slower than FNNOLS. The performance of theproposed techniques, i.e. FNNOLS and SNNOLS, have tobe investigated, using other types of dictionaries. The com-putational complexity analysis of the FNNOLS is also aninteresting subject, which we left for the future work.AcknowledgmentThe authors acknowledge the support of EPSRC and theMOD University Defence Research Collaboration in SignalProcessing, through the grant number EP/K014277/1.REFERENCES[1] Y Qian, S Jia, J Zhou, and A Robles-Kelly, “Hyper-spectral unmixing via sparsity-constrained nonnegativematrix factorization,” Geoscience and Remote Sensing,IEEE Transactions on, vol. 49, no. 11, pp. 4282–4297,2011.[2] F. Parvaresh, H. Vikalo, S. Misra, and B. Hassibi, “Re-covering sparse signals using sparse measurement matri-ces in compressed DNA microarrays,” IEEE Journal ofSelected Topics in Signal Processing, vol. 2, no. 3, pp.275–285, 2008.[3] PO Hoyer, “Non-negative sparse coding,” in Neural Net-works for Signal Processing, 2002. Proceedings of the2002 12th IEEE Workshop on. 2002, pp. 557–565, IEEE.[4] M. Yaghoobi, D. Wu, and M.E. Davies, “Fast non-negative orthogonal matching pursuit,” IEEE Signal Pro-cessing Letters, vol. 22, no. 9, 2015.[5] S. Chen, S.A. Billings, and W. Luo, “Orthogonal leastsquares methods and their application to non-linear sys-tem identification,” International Journal of Control, vol.50, no. 5, pp. 1873–1896, 1989.[6] B.K. Natarajan, “Sparse approximate solutions to linearsystems,” SIAM Journal of Comput, vol. 24, no. 2, pp.227–234, 1995.[7] G.H. Golub and C.F. Van Loan, Matrix computations,Johns Hopkins University Press Baltimore, 1996.[8] T Blumensath and ME Davies, “On the difference be-tween orthogonal matching pursuit and orthogonal leastsquares,” Tech. Rep., 2007.

Fast non-negative orthogonal least squares

https://www.pure.ed.ac.uk/ws/files/19708315/YD_FastNNOLS_EUSIPCO_2015.pdf

Fast non-negative orthogonal least squares

Abstract

Similar works

Full text

Available Versions

Edinburgh Research Explorer