Thi Thanh Nguyen

Member, IEEE Jérôme Idier
email: jerome.idier@ls2n.fr

Member, IEEE Charles Soussen
email: charles.soussen@centralesupelec.fr

El-Hadi Djermoune
email: el-hadi.djermoune@univ-lorraine.fr

Thanh T Nguyen
email: thi-thanh.nguyen@univ-lorraine.fr

Member, IEEE El-Hadi Djermoune

T T Nguyen

E.-H Djermoune

Non-negative orthogonal greedy algorithms

Keywords: Orthogonal greedy algorithms, sparse reconstruction, non-negativity, non-negative least-squares, active-set algorithms

I. INTRODUCTION

In the last decade, sparse approximation has received considerable attention in the signal and image processing community, in connection with compressive sensing. Greedy algorithms for sparse signal reconstruction are very popular iterative schemes. Their principle is to repeatedly (i) enrich the sparse support by selecting a new dictionary atom, and then (ii) update the sparse approximation coefficients. In orthogonal greedy algorithms, the sparse approximation signal is computed as the orthogonal projection of the data vector onto the subspace spanned by the selected atoms. Therefore, the coefficients can be estimated by solving an Unconstrained Least Squares (ULS) problem. Popular orthogonal greedy algorithms include Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] and Orthogonal Least Squares (OLS) [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF], also known as forward selection [START_REF] Miller | Subset selection in regression[END_REF], Order Recursive Matching Pursuit (ORMP) [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF], Optimized Orthogonal Matching Pursuit (OOMP) [START_REF] Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF] and Pure Orthogonal Matching Pursuit [START_REF] Foucart | Stability and robustness of weak orthogonal matching pursuits[END_REF]. OMP and OLS differ in the way the new atom is selected. In both cases, the atom inducing the largest decrease of the residual norm is selected. However, all nonzero atom weights are optimally tuned in OLS whereas only the new atom weight is considered in OMP, which amounts to selecting the atom having the largest inner product with the current residual. The computational complexity of OLS is obviously higher, since the selection rule requires to solve as many ULS problems as the number of candidate atoms. Fortunately, the ULS solutions have a closed-form expression, which can be recursively (fastly) updated when the support is enriched by a new element, see e.g., [START_REF] Miller | Subset selection in regression[END_REF]. Specifically, both OMP and OLS implementations are recursive and make use of matrix factorization, such as Gram-Schmidt orthogonalization, the Cholesky factorization or techniques utilizing the matrix inversion lemma [START_REF] Sturm | Comparison of orthogonal matching pursuit implementations[END_REF].

Many application fields such as geoscience and remote sensing [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF], [START_REF] Esser | A method for finding structured sparse solutions to nonnegative least squares problems with applications[END_REF], audio [START_REF] Virtanen | Active-set Newton algorithm for overcomplete non-negative representations of audio[END_REF], chemometrics [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF], bioinformatics [START_REF] Slawski | Isotope pattern deconvolution for peptide mass spectrometry by non-negative least squares/least absolute deviation template matching[END_REF], astrophysics [START_REF] Högbom | Aperture synthesis with a non-regular distribution of interferometer baselines[END_REF] and computed tomography [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF], give rise to inverse problems where the signal or image of interest is sparse, but also non-negative. In such contexts, a common practice is to regularize the inverse problem in order to favor both sparsity and non-negativity of the sought signal, see, e.g., [START_REF] Slawski | Isotope pattern deconvolution for peptide mass spectrometry by non-negative least squares/least absolute deviation template matching[END_REF], [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF]- [START_REF] Rapin | Sparse and non-negative BSS for noisy data[END_REF]. Some classical sparse algorithms can be straightforwardly adapted to deal with nonnegativity constraints. This is the case of proximal splitting algorithms and the Alternating Direction Method of Multipliers (ADMM) for convex optimization [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], and of the DC algorithm (Difference of Convex functions) for nonconvex optimization [START_REF] Gasso | Recovering sparse signals with a certain family of nonconvex penalties and DC programming[END_REF], [START_REF] Le Thi | Sparse signal recovery by difference of convex functions algorithms[END_REF]. On the contrary, the non-negative extension of greedy algorithms is a challenging issue since the unconstrained least-squares subproblems are replaced by nonnegative least-squares (NNLS) subproblems which do not have closed-form solutions anymore, so a subroutine solving NNLS is needed. There are different methods for NNLS solving [START_REF] Björck | Numerical Methods for Least Squares Problems[END_REF] such as active-set [START_REF] Lawson | Solving least squares problems[END_REF], interior-point [START_REF] Wright | Interior methods for constrained optimization[END_REF], and gradientprojection [START_REF] Benvenuto | Nonnegative leastsquares image deblurring: Improved gradient projection approaches[END_REF]. The latter two families typically require to tune some stopping criteria in an empirical manner, resulting in approximate resolution of the NNLS problem. Here, we are focusing on active-set methods since such methods have a greedy structure and exactly solve NNLS problems after a finite number of iterations. Although our focus will be on extensions of orthogonal greedy schemes to the non-negative case, let us mention that several other non-negative sparse methods have been elaborated on the basis of the active-set NNLS algorithm, e.g., hard-thresholded NNLS [START_REF] Slawski | Sparse recovery by thresholded non-negative least squares[END_REF], [START_REF] Slawski | Non-negative least squares for highdimensional linear models: Consistency and sparse recovery without regularization[END_REF] and Sparse NNLS or Reverse Sparse NNLS [START_REF] Peharz | Sparse nonnegative matrix factorization with 0 constraints[END_REF].

The classical CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], Subspace Pursuit (SP) [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] and Hard Thresholding Pursuit [START_REF] Foucart | Stability and robustness of weak orthogonal matching pursuits[END_REF] algorithms were extended to the non-negative setting in [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]. Similar to orthogonal greedy algorithms, the non-negative versions NN-CoSaMP, NN-SP and NN-HTP are based on straightforward replacements of ULS resolutions by NNLS resolutions. The principle of these algorithms is to update a support of fixed size at consecutive iterations using operations such as support expansion and support shrinkage. Since two consecutive supports are possibly substantially different, one can hardly exploit the knowledge of the current NNLS solution to compute the solution related to the new support in a very efficient manner. On the contrary, the supports found at consecutive iterations of orthogonal greedy algorithms are nested, allowing one to use warm start initialization. The nested property is therefore crucial to design fast implementations.

Several existing contributions deal with orthogonal greedy algorithms in the non-negative case. Non-Negative OMP (NNOMP) was first proposed by Bruckstein et al. [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] as a direct generalization of OMP. At each iteration, the atom having the maximum positive inner product with the current residual is selected. Contrary to OMP, negative inner products are discarded. Then, the sparse approximation coefficients are updated by solving the NNLS problem related to the augmented subset. The canonical (i.e., non-recursive) NNOMP implementation of Bruckstein et al. [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] solves NNLS subproblems independently. Later, Yaghoobi et al. proposed an accelerated version named Fast Non-Negative OMP (FN-NOMP), which avoids solving NNLS subproblems but rather recursively approximates the sought solution using QR matrix factorization [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF]. Although FNNOMP is much faster than canonical NNOMP, it is an approximate version likely to deliver a different output. In [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF], Yaghoobi et al. introduced a canonical version of Non-Negative OLS (NNOLS), defined as a direct non-negative extension of OLS. The principle of NNOLS is to select the atom for which the positive residual (the residual corresponding to the NNLS solution) is minimum. This selection rule appears to be time-consuming since one needs to compute as many positive residuals as the number of candidate atoms, i.e., n -k NNLS problems have to be solved at iteration k, with n the size of the dictionary. Yaghoobi et al. [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF] further proposed two accelerated versions of NNOLS named Suboptimal NNOLS (SNNOLS) and Fast NNOLS (FNNOLS). These proposals can be understood from a dual interpretation of OLS in the unconstrained setting. It is well-known than the OLS selection rule can be formulated in two equivalent ways [START_REF] Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF], [START_REF] Blumensath | On the difference between Orthogonal Matching Pursuit and Orthogonal Least Squares[END_REF]: OLS selects the atom (i) inducing the minimum value of the data residual in 2 -norm; (ii) whose orthogonal projection onto the subspace orthogonal to the active atoms has a maximum angle with the current data residual. When dealing with non-negativity, SNNOLS [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF] selects the atom that is positively correlated with the current residual whose projection forms a maximum angle with the residual. However, the equivalence that holds for OLS is not true anymore since maximizing the angle for positive inner products is not equivalent to minimizing the 2 -norm of the positive residual. Therefore, SNNOLS iterates do not generally identify with NNOLS iterates, hence the suboptimal algorithm denomination. In SNNOLS, the projected atoms have to be recomputed whenever the support is updated. Moreover, an NNLS problem must be solved at each iteration to update the sparse approximation coefficients. FNNOLS is a recursive implementation in the same spirit as FNNOMP, where no NNLS problem needs to be solved anymore. It shall be noticed that FNNOLS and SNNOLS do not necessarily deliver the same iterates, and that both can be viewed as approximate versions of NNOLS.

Generally speaking, the "orthogonal" denomination of NNOMP and NNOLS is somewhat abusive since when the support set S is updated, the related NNLS solution may not identify with the orthogonal projection of the data vector onto the span of atoms indexed by S, but rather with its projection onto their positive span. Both projected vectors differ as soon as the NNLS solution has zero entries, i.e., when some nonnegativity constraints become active. Therefore, in NNOMP, NNOLS and their derived versions [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF]- [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF], the support of the sparse vector at the current iteration is a subset of the current support set S (which is expanded at each iteration) and may not identify with it. In turn, more than K iterations may be necessary to reach a truly K-sparse representation.

Our contributions are twofold. First, non-negative orthogonal greedy algorithms are revisited in a unified way. The algorithms under study share four desirable properties:

1) The norm of the data residual is always decreasing when a new atom enters the solution support.

2) The algorithm does not stop while additional atom selections would make it decrease, unless an explicit stopping condition is reached. 3) A so-called compression step is included to shrink the support set by removing the atoms having zero coefficients, so the support set identifies to the support of the current NNLS solution. 4) The residual vector is orthogonal to the selected atoms.

In other words, the sparse approximation vector identifies with the orthogonal projection of the data vector onto the span of the selected atoms.

These structural properties are exhibited and compared to those of existing non-negative greedy algorithms. The second contribution is a fast and exact implementation of nonnegative orthogonal greedy algorithms exploiting the active-set algorithm [START_REF] Lawson | Solving least squares problems[END_REF] for NNLS solving, based on warm start initialization. Moreover, we elaborate on recursive implementations and design further types of accelerations.

The paper is organized as follows. Section II introduces the family of so-called Non-negative Orthogonal Greedy (NNOG) algorithms. The different members of the family differ by the selection rule to pick a new atom at each iteration. It includes NNOLS, SNNOLS and NNOLS up to a modification of their structure, namely the compression step mentioned above. In Section III, a fast implementation is proposed based on recursivity and on the use of warm starts for solving the NNLS subproblems. Section IV is devoted to NNOG acceleration. Sections V and VI propose numerical validations on a simulated example of sparse deconvolution and on the decomposition of real-world near-infrared (NIR) spectra. Finally, discussion and perspectives will be found in Section VII.

II. NON-NEGATIVE GREEDY ALGORITHMS

A. Basic definitions and notations

Given a data vector y ∈ R m and a dictionary H ∈ R m×n , we are interested in finding a K-sparse non-negative weight vector x ∈ R n + yielding an accurate approximation y ≈ Hx. This can be formulated as the constrained minimization program:

min x≥0 y -Hx 2 2 s.t.
x 0 ≤ K.

( 0 +)
where x ≥ 0 means that each entry of x is nonnegative, x 0 is the 0 -"norm" counting the number of nonzero entries in x, and the quadratic fidelity-to-data term y -Hx 2 2 measures the quality of approximation. The 2 -norm • 2 will be also denoted • . Without loss of generality, we assume that H is column-normalized. Each dictionary column h i , i = 1, . . . , n is called an atom. We have the following useful identity for any two vectors r, h of same length, h being normalized:

min v≥0 r -hv 2 = r 2 -max{h t r, 0} 2 (1)
where . t stands for the transpose operator. We denote by S = supp(x) = {i : x(i) = 0} the support of x (x(i) being the i-th entry of x), S the complement of S, |S| the cardinality of S, H S and x(S) the subdictionary and subvector indexed by S, respectively. Finally, H † and span(H) are the pseudo-inverse and the column space of H, respectively. Let hi,S = h i -H S H † S h i stand for the orthogonal projection of h i onto the orthogonal complement (span(H S )) ⊥ , which will be simply denoted hi whenever unambiguous, and gi = hi / hi denote the normalized projected atom if h i ∈ span(H S ), i.e., hi = 0. If h i ∈ span(H S ), it will be convenient to set gi = 0.

For any support S, let us call an unconstrained least-squares (ULS) and a nonnegative least-squares (NNLS) solution corresponding to S, any vector x in R n or R n + , respectively, that minimizes y -Hx 2 subject to supp(x) ⊂ S. Such vectors will be denoted xS and x+ S , respectively. The following notations will be also useful:

r S = y -H xS , r + S = y -H x+ S . When H S is full column rank, y -H S z 2 is a strictly convex function of z ∈ R |S| , so x+
S and xS are then uniquely defined. Throughout the paper, we will denote by C ⊂ S the so-called compressed support, defined as the support of x+ S . The NNLS optimal solutions can be characterized using the Karush-Kuhn-Tucker (KKT) conditions [START_REF] Lawson | Solving least squares problems[END_REF]Chap. 3], which are recalled next for completeness. (

) 2
x+ S is a solution to (2) if and only if the KKT conditions are satisfied:

H t C (y -H x+ S ) = 0 H t S\C (y -H x+ S ) ≤ 0 (3)
where C := supp(x + S ) ⊂ S.

Proof: From the definition of C, it is clear that x+ S (C) > 0, so the active constraints in (2) are indexed by C. Let λ ∈ R n gather the Lagrange multipliers related to both equality and inequality constraints. The Lagrangian function induced by ( 2) is defined as L(x; λ) = y -Hx 2 -λ t x and the KKT conditions for optimal variables (x + S , λ) read:

   ∇ x L(x + S ; λ) = 2H t (H x+ S -y) -λ = 0, ∀i ∈ S, λ(i)x + S (i) = 0 with x+ S (i) ≥ 0, λ(i) ≥ 0, ∀i / ∈ S, x+ S (i) = 0. (4)
For quadratic programming problems involving positive semidefinite matrices, the KKT conditions are necessary and sufficient conditions of optimality [START_REF] Nocedal | Numerical optimization[END_REF]Chap. 16], so x+ S is a solution to (2) if and only if

   λ = -2H t (y -H x+ S ) λ(C) = 0 λ(S\C) ≥ 0 that is, when (3) is satisfied.
Definition 1 Let us call a positive support related to the fullsize NNLS problem

min x≥0 y -Hx 2 , (5)
any index set S such that H S is full column rank and x+ S (S) > 0. By extension, the empty support S = ∅ will also be considered as a positive support.

Lemma 2 S is a positive support if and only if H S is full rank and xS (S) > 0. Moreover, when S is a positive support, x+ S = xS , r + S = r S and H t S r S = 0.

Proof: When S is a positive support, ULS and NNLS solutions xS and x+ S are uniquely defined and coincide. The orthogonality property H t S r S = 0 follows from (3). Positive supports will play an important role in our specification of fast non-negative orthogonal greedy algorithms.

B. Non-negative orthogonal greedy algorithms

Let us define the class of non-negative orthogonal greedy (NNOG) algorithms, sharing the following general structure. We start from the empty support S = ∅. At each iteration, an atom is moved from S to S. A new NNLS solution x+ S is then computed to optimally adapt the weights to the newly extended support. The algorithm stops when the desired cardinality K is reached or when the norm of the residual cannot decrease anymore. The general structure of NNOG algorithms is given by Algorithm 1. Some aspects will be made clear later, such as the role of the test h t i r S > 0 with respect to the decrease of the norm of the residual.

The NNOG class is a direct adaptation of the family of orthogonal greedy algorithms from the unconstrained case to the nonnegative one. At first glance, the two families only differ by the fact that an NNLS solution is computed rather than a ULS one to update the weights at each iteration. However, some important features differ between the two cases, which require non trivial adaptations.

In both cases, the greedy character corresponds to the fact that a unique atom is added to the current support per iteration. However, a distinct feature of NNOG algorithms is that the support size may be smaller than the current iteration index, because some components of x+ S (S) may vanish due to the activation of the corresponding nonnegativity constraints. In the unconstrained case, some components of xS (S) may also vanish, but such events are fortuitous and do not need any specific consideration.

In the unconstrained case, r S is orthogonal to span(H S ). This geometrical property does not necessarily hold for r + S because NNLS is an inequality constrained problem. Fortunately, provided that the indices of zero components of x+ S are moved to S, it remains true that y -H x+ S is orthogonal to span(H S ). This is a direct consequence of the following lemma, which states that x+ S reads as a ULS solution related to the compressed version of support S. Finally, a practically important aspect is the computing cost of NNOG algorithms. It is computationally more demanding to solve an NNLS problem than the corresponding ULS problem, so one must expect a larger computing cost. However, NNOG algorithms lend themselves to recursive implementations akin to usual orthogonal greedy schemes, as detailed in Section III.

C. Descending atoms and descent selection rules

Greedy algorithms can be interpreted as descent algorithms dedicated to the minimization of the residual norm using supports of growing size. Contrary to the unconstrained case, only the selection of some atoms in S may produce a decrease of the residual at a given iteration of an NNOG algorithm. For the rest of the atoms in S, the residual norm decrease is possible only if the nonnegativity constraint is violated. Such a specificity has both formal and practical consequences on the design of NNOG algorithms, as examined in this subsection.

Definition 2 For a given support S, let us define the set of indices corresponding to descending atoms as follows:

D S = i ∈ {1, . . . , n}, r + S∪{i} < r + S .
Obviously, we have D S ⊂ S. In what follows, we focus on selection rules ensuring the selection of a descending atom at any iteration. The latter rules are referred to as descent selection rules.

Definition 3 A descent selection rule is a function S(y, H, S) that takes its values in D S .
Clearly, NNOG algorithms relying on a descent selection rule are descent algorithms and we have the following property.

Lemma 4 NNOG algorithms relying on a descent selection rule terminate after a finite number of iterations.

Proof: The error norm decreases at each iteration, and there is a finite number of supports with a cardinality not exceeding K, and thus a finite number of solutions to visit.

The following proposition allows one to characterize whether an atom is descending.

Proposition 1 The descending atom condition i ∈ D S is equivalent to 0 < h t i r + S . ( 7
)
When S is a positive support, it is also equivalent to each condition

0 < h t i r S , (8) 0 < gt i r + S , (9) 0 < gt i r S . (10)
When H S∪{i} is full column rank, it is also equivalent to

x+ S∪{i} (i) > 0. (11)
Proof: See Appendix A-A. Let us remark that from the first item of Proposition 1, the selection rule is invoked at Line 3 of Algorithm 1 only if D S = ∅, otherwise the stopping condition of Line 2 is activated. Hence, we do not need to define S(y, H, S) when D S = ∅.

The following lemmas have an interesting impact on the practical specification of NNOG algorithms.

Lemma 5 If S is a positive support and i ∈ D S , then matrix H S∪{i} is full column rank.

Proof: Assume that S is positive, so H S is full column rank, and x+ S = xS . Let us also assume that H S∪{i} is not full column rank. Then h i ∈ span(H S ), so xS is a ULS solution related to S ∪ {i}, and also an NNLS solution related to S ∪ {i} since xS ≥ 0. This implies that i ∈ D S .

Lemma 6 After each iteration of an NNOG algorithm relying on a descent selection rule, it holds that the support of the current solution is positive.

Proof: The proof is immediate by recursive application of Lemmas 3 and 5, starting with the empty support.

Let us stress that Lemma 6 refers to NNOG algorithms strictly conforming to the scheme of Algorithm 1 (with an additional restriction to descent selection rules at Line 3). In particular, the support compression step performed at Line 6 is necessary to make Lemma 5 applicable. To our best knowledge, such a compression step has not been proposed in any previous contribution about nonnegative greedy schemes targeting 0 minimization1 .

According to Lemma 6, the restriction to a descent selection rule implies r + S = r S at any iteration, which justifies that we have dropped the '+' sign in Algorithm 1. This simplification is adopted in the rest of the paper. Moreover, the termination rule max i h t i r S ≤ 0 is used at Line 2 since in this case, there are no descending atoms anymore, so the residual cannot decrease by selection of a new atom.

D. Examples of descent selection rules

In what follows, selection rules are denoted S(S), the dependence on y and H being implicit. Let us introduce three important selection rules by their distinct ways of picking an index in D S when D S = ∅.

• NNOMP rule [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]- [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF], [START_REF] Ramamurthy | Recovering nonnegative and combined sparse representations[END_REF]:

S 1 (S) ∈ arg max i / ∈S h t i r S (12)
• Suboptimal NNOLS (SNNOLS, [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF]) rule:

S 2 (S) ∈ arg max i / ∈S gt i r S (13)
• NNOLS rule [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF]:

S 3 (S) ∈ arg min i / ∈S r + S∪{i} 2 (14)
(14) is a descent selection rule by definition. The fact that ( 12) and ( 13) are descent selection rules is deduced from recursive application of Proposition 1 and Lemma 6. As regards the latter rules, r S is the current residual vector, i.e., a readily available quantity. On the other hand, projected atoms gi enter rule ( 13), so we can expect the computing cost of ( 13) to be larger than that of [START_REF] Slawski | Isotope pattern deconvolution for peptide mass spectrometry by non-negative least squares/least absolute deviation template matching[END_REF]. Rule ( 14) needs the solution of NNLS problems on supports S ∪ {i}, which is even more demanding.

Note that [START_REF] Wang | Cone-based joint sparse modelling for hyperspectral image classification[END_REF] introduced another version of non-negative OMP named NN-OMP in which the selection rule is that of OMP. Clearly, this version does not rely on a descent selection rule. On the other hand, it is unclear whether the FNNOMP and FNNOLS algorithms proposed in [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF], [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF] rely on a descent selection rule. They will therefore not be further analyzed.

The following proposition makes it possible to compare the three rules ( 12)-( 14) by relating them to the minimization of a residual norm.

Proposition 2 Rules (12)-( 14) are equivalent to

S j (S) ∈ arg min i∈D S µ j (S, i), (15)
where µ j are specific to each rule:

µ 1 (S, i) = min v y -H S x S -h i v 2 , (16)
µ 2 (S, i) = min u,v y -H S u -h i v 2 , (17)
µ 3 (S, i) = min u≥0,v≥0 y -H S u -h i v 2 . (18)
Proof: Let us first emphasize that because ( 12)-( 14) are descent selection rules, i / ∈ S in ( 12)-( 14) can be replaced by i ∈ D S . Clearly, [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] with [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] simply duplicate [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF]. Using identity (1) with r = r S and h = h i , ( 12) can be rewritten as the argmin over i ∈ D S of

min v≥0 y -H S x S -h i v 2 . ( 19
)
Likewise, applying (1) with r = r S and h = gi and using the fact that gi is the normalized version of hi , (13) can be rewritten as the argmin over i ∈ D S of

min v≥0 r S -hi v 2 = min u,v≥0 y -H S u -h i v 2 . (20)
(20) follows from explicit minimization with respect to u and from the fact that r S and hi read as the orthogonal projections of y and h i onto (span(H S )) ⊥ . Finally, the positivity constraint on v in ( 19)-( 20) turns out to be inactive. Indeed, according to Proposition 1, we have that h t i r S > 0 and ht i r S > 0 for i ∈ D S . It is easy to see from (1) that the minimum error norm in ( 19)-( 20) is then strictly lower than r S . Hence, the optimal variable v in ( 19)-( 20) is nonzero.

Finally, [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] with ( 19) and ( 15) with (20) rewrite as ( 15) with ( 16) and ( 15) with [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], respectively.

For j ∈ {1, 2, 3}, an alternate way of viewing the descending character of each rule consists in noticing that r S∪{i} 2 ≤ µ j (S, i) < r S 2 for all i ∈ D S . It is interesting to see that, by restricting the selection to the set of descending atoms, the selection rules of NNOMP and SNNOLS rely on the same criteria ( 16)-( 17) as those of OMP and OLS, respectively.

E. On the usefulness of support compression

As already mentioned, rules ( 12)-( 14) were introduced in existing works devoted to problem ( 0 +) without support compression. Let us analyze the impact of removing the support compression step at Line 6 in Algorithm 1.

• The subset S then becomes an extended support, with x = x+ S and supp(x) ⊂ S. As a consequence, when supp(x) = S, the current residual r + S = y -Hx does not identify to r S anymore since x does not read as a ULS solution related to S. • In the case of NNOMP and NNOLS (selection rules [START_REF] Slawski | Isotope pattern deconvolution for peptide mass spectrometry by non-negative least squares/least absolute deviation template matching[END_REF] and ( 14), with r S replaced by the current residual r + S ), Proposition 1 can still be applied, since condition h t i r + S > 0 still characterizes descending atoms without support compression. As for SNNOLS, we cannot guarantee that S(y, H, S) defined in ( 13) is a descent selection rule when S is not positive, since the equivalence (9) in Proposition 1 does not hold anymore.

• Finally, for any NNOG algorithm, Lemma 5 is no more valid when S is not a positive support, so H S might become rank deficient at any subsequent iteration, x+ S not being well defined anymore. We consider the last point as the main justification to apply support compression. An open question would be to determine whether the assumptions of Lemma 5 can be weakened, the positivity of S being replaced by the full rankness of H S . However, we lack concrete elements in this sense. Therefore, at this stage, NNOG with compression present formal guarantees that NNOG without compression do not.

III. ACTIVE-SET NNLS ALGORITHMS AND RECURSIVITY

Let us consider the NNLS problem related to support S in [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF]. When matrix H S is full column rank, it is a special case of a strictly convex quadratic program. The successive resolution of possibly many NNLS problems for nested supports S is a basic ingredient of NNOG algorithms. NNLS problems do not admit closed-form solutions in general. However, activeset NNLS (AS-NNLS) algorithms are well-known schemes that solve NNLS problems in a finite number of iterations [START_REF] Lawson | Solving least squares problems[END_REF], [START_REF] Chen | Lower bound theory of nonzero entries in solutions of 2p minimization[END_REF]. Moreover, the computation of NNLS solutions can be efficiently accelerated by exploiting that warm start solutions are available and using classical recursive updates of matrix factorizations.

This section first contains a short reminder on active-set NNLS algorithms and on their efficient implementation. Then, we show how to preserve computational efficiency when NNLS subproblems are solved within an NNOG scheme. We also analyze the similarity between the structures of AS-NNLS algorithm and NNOMP, already pointed out in [START_REF] Peharz | Sparse nonnegative matrix factorization with 0 constraints[END_REF], [START_REF] Foucart | Sparse recovery by means of nonnegative least squares[END_REF].

A. Fast active-set algorithms Among many numerical methods for solving (2), AS-NNLS algorithms correspond to greedy schemes, since the solution is found by incremental modifications of its support V ⊂ S, the active set being defined as the complementary set S\V (by reference to the active constraints). Such an incremental Algorithm 2: Active-set algorithm to solve the NNLS problem related to T , starting from a positive support S.

input : y, H, target set T , initial support S ⊂ T , xS output:

V := supp(x + T ), x+ T := xV 1 x ← xS ; V ← S ; 2 while max{h t i r V , i ∈ T \V } > 0 do 3 + ← arg max{h t i r V , i ∈ T \V } ; 4 V ← V ∪ { + }; 5 Update xV (call Algorithm 4); 6 while min(x V ) < 0 do 7 -∈ arg min {i∈V :x V (i)<0}
x(i)/(x(i) -xV (i)) ;

8 α ← x( -)/(x( -) -xV ( -)) ; 9 x ← x + α(x V -x) ; 10 V ← V \ { -} ; 11 Update xV (call Algorithm 4); 12 end 13 x ← xV ; 14 r V ← y -Hx ; 15 end
structure is an essential element to obtain practically fast implementations [START_REF] Björck | Numerical Methods for Least Squares Problems[END_REF], [START_REF] Lawson | Solving least squares problems[END_REF], [START_REF] Nocedal | Numerical optimization[END_REF]. Whenever V is modified, the corresponding ULS solution xV is updated. Each iteration requires at least one ULS solution of the selection type, i.e., H † V ∪{ } y. Some iterations also need to compute deselection type ULS solutions H † V \{ } y. For the sake of computational efficiency, it is crucial to compute both types of solutions recursively, given that xV is already available. In this respect, the situation is identical to that of bi-directional greedy algorithms in the unconstrained case (e.g., Bayesian OMP [START_REF] Herzet | Bayesian pursuit algorithms[END_REF], SBR [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] or FoBa [START_REF] Zhang | Sparse recovery with orthogonal matching pursuit under RIP[END_REF]), also called stepwise algorithms [3, Chapter 3], for which selections and deselections are implemented recursively. Fast recursive implementations require specific computations and storage of quantities related to the Gram matrix G V = H t V H V . Efficient selection steps can be obtained using QR or Cholesky matrix factorizations applied to G V , or the Matrix Inversion Lemma (MIL), with roughly the same cost [START_REF] Sturm | Comparison of orthogonal matching pursuit implementations[END_REF]. MIL consists in storing and updating the inverse of G V . It appears to be the cheapest concerning deselections, so our default choice here is based on the MIL.

B. Warm start and full recursivity

Since NNOG algorithms are based on iterated calls to an active-set scheme, we must carefully consider the way we initialize the latter. The starting point of AS-NNLS is usually defined as the zero solution, associated with the empty support. This is the case in the Lawson-Hanson algorithm, which is the reference AS-NNLS scheme [START_REF] Lawson | Solving least squares problems[END_REF]. In [START_REF] Nocedal | Numerical optimization[END_REF]Chap. 16], Nocedal and Wright proposed an AS-NNLS algorithm with any feasible vector as a possible starting point. Algorithm 2 is a generalized version of the Lawson-Hanson algorithm to address the NNLS problem related to an augmented target set T ⊃ S, where the initial point is not restricted to be the zero vector (the rest of the scheme being unaltered). Specifically, the initial point is set as the ULS solution xS ≥ 0, S being a positive support.

For this specific initial point, it can easily be checked that Algorithm 2 identifies with Nocedal and Wright's scheme. At any iteration of an NNOG algorithm, the current solution x+ S can be used as the initial point to compute x+ S∪{ } using Algorithm 2 with T = S ∪ { }. In this way, the initial support is V = S, so T \V = { } and the first iteration of AS-NNLS begins by selecting . In practice, the NNLS algorithm is expected to terminate after a single iteration (if no deselection is performed in the second part of it), or otherwise after very few iterations. Algorithm 3 is a global view of the active-set implementation of NNOG obtained by integrating the calls to the AS-NNLS solver (Algorithm 2) in the NNOG framework (Algorithm 1). Whenever a new atom is selected, AS-NNLS starts by computing xS∪{ } . If xS∪{ } ≥ 0, then x+ S∪{ } = xS∪{ } and AS-NNLS stops after one support change. Otherwise, AS-NNLS deselects at least one atom from S ∪ { } (Algorithm 2, Line 7) and then alternates between atom selections and series of deselections. This mechanism is illustrated by a simple example in the next subsection.

A reduced number of iterations is obtained because we use Algorithm 2 with a warm start. To further improve the overall numerical efficiency, we also need to reduce the computing cost of the ULS solution at Lines 5 and 11. These are selection and deselection-type ULS problems, respectively, that can be solved recursively provided that the inverse of the Gram matrix G V be both an input and an output quantity of the NNLS algorithm. In Appendix B, Algorithm 4 is a pseudo-code to implement ULS updates in the forward (V ← V ∪ { }) and backward scenarios (V ← V \{ }). This implementation enables us to obtain a fully recursive implementation of AS-NNLS as well by updating the Gram matrix inverse at each call to ULS in Algorithm 2 (Lines 4-5 and 10-11). Here again, S 4 ∪ { } is not a positive support. The first deselection does not yield a positive support either, so another deselection is carried out, yielding V ← S 4 ∪ { }\{ 1 , 2 }, V being a positive support. However, the stopping condition of AS-NNLS (Line 2 of Algorithm 2, with T ← S 4 ∪ { }) is not met since h 1 is a descending atom. Therefore, 1 is reselected:

C. Step-by-step illustrative example of NNOG

V ← V ∪ { 1 }.
Since V is now a positive support and there are no descending atoms anymore in T \V , the convergence of AS-NNLS is reached. In the last two NNOG iterations, S k-1 ∪ { } are positive supports, so single selection moves are done within AS-NNLS.

D. Connection between AS-NNLS and NNOMP

In Algorithm 3, the structure of the NNOG main iteration consists of one atom selection followed by a variable number of updates (selections or deselections) of the support S when the ULS solution xS has some negative entries. The first of these updates is necessarily an atom deselection. Interestingly, [START_REF] Peharz | Sparse nonnegative matrix factorization with 0 constraints[END_REF] and [START_REF] Foucart | Sparse recovery by means of nonnegative least squares[END_REF] pointed out that the AS-NNLS algorithm initialized with the zero vector, has a structure similar to NNOMP: each main NNLS iteration consists of an atom selection followed by a series of atom deselections (respectively, Line 3 and Lines 6-12 of Algorithm 2). This connection led Peharz and Pernkopf to propose AS-NNLS with an early stopping rule |S| = K as a stand-alone NNOG algorithm (called Sparse NNLS in [START_REF] Peharz | Sparse nonnegative matrix factorization with 0 constraints[END_REF]). Given the strong similarity of Sparse NNLS and NNOMP, an interesting question is to determine whether their iterates always coincide. It turns out that this is not always true. However, as long as Sparse NNLS performs only simple support changes, both algorithms yield the same iterates, according to the following proposition.

Proposition 3 Let us consider Sparse NNLS, i.e., Algorithm 2 initialized with the empty support together with the early stopping rule |S| = K. When each iteration produces at most one removal in Lines 6-12, the output of Sparse NNLS identifies with that of NNOMP (i.e., Algorithm 1 with rule S 1 in ( 12)).

Proof: See Appendix A-B.

IV. ACCELERATION OF NNOG ALGORITHMS

The NNOMP, SNNOLS and NNOLS selection rules ( 12)-( 14) all read as the optimization of a criterion with respect to the candidate index i / ∈ S. Since all three ensure the selection of an atom in D S , an obvious acceleration of SNNOLS and NNOLS consists of pre-selecting the descending atoms according to h t i r S > 0 (see Proposition 1) to carry out the optimization tasks ( 13) and ( 14) over i ∈ D S only. This operation is referred to as type-I pruning of the dictionary2 . Testing the sign of h t i r S requires O(m) operations. This is much less than the (recursive) computation of the criteria gt i r S and r + S∪{i} 2 , which costs at least O(|S| 2 + km) operations.

A. Atom selection

The atom selection step of both NNOMP and SNNOLS can be efficiently implemented using vectorized computations, which allow one to benefit from the inherent parallelism of SIMD processors. The NNOMP case is the simplest, since H t S r S directly yields the expected set of inner products h t i r S . Indeed, the selection steps of OMP and NNOMP are both based on the minimization of ( 16), so they share the same possibilities of parallel computations. Likewise, OLS and SNNOLS being both based on the minimization of ( 17), they share the same possibilities of parallel implementation. In coherence with our choice of recursive implementation of ULS solutions (see Appendix B), we have adopted a MIL based solution to solve [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] in a vectorized way (see Matlab code in supplementary material).

In contrast, the NNOLS selection rule [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF] does not lend itself to fully vectorized computations, since we have as many NNLS subproblems to solve as candidate atoms, with a variable number of subiterations of AS-NNLS for of them. Fortunately, the structure of the NNOLS rule can be made closer to that of SNNOLS, with the benefit of vectorized computations for the largest part. The key point is that for each candidate atom, the initial step of AS-NNLS corresponding to Lines 1-5 of Algorithm 2 yields the same unconstrained minimizer xS∪{i} as the one involved in the SNNOLS rule [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Hence, these vectors can be obtained using vectorized computations that exactly identify to the main step of SNNOLS atom selection. The extra computations induced by NNOLS reside in the additional AS-NNLS iterations required for each non-positive support S ∪ {i}. According to our empirical tests, only a small minority of atoms needs more than one iteration. Moreover, a lot of them can be pruned without actually performing any additional AS-NNLS iterations. Let us denote by e opt the smallest residual error produced by atoms i for which S ∪ {i} is a positive support. Since r + S∪{i} ≥ r S∪{i} for all i, one can immediately ignore the atoms i for which r S∪{i} ≥ e opt . Moreover, since we sequentially visit the remaining ones, the threshold e opt can be gradually lowered whenever a new atom is found to improve the smallest current residual error. This operation called type-II pruning is specific to NNOLS implementation.

B. Coefficient update

Once an atom i is selected, NNOG algorithms need to update the coefficients by solving an NNLS subproblem (Algorithm 1, Line 5). However, in the case of NNOLS, the update is already being performed in the selection step. In the case of NNOMP, one needs to call AS-NNLS (Algorithm 2) from the initial set S to the target set S ∪ { }. For SNNOLS, a call to Algorithm 2 is needed only when S ∪ { } is not positive.

C. Software

An open source Matlab implementation of the acceleration strategies detailed above is provided on CodeOcean [START_REF] Nguyen | Non-negative greedy algorithms: Matlab implementation[END_REF] (https://doi.org/10.24433/CO.2445681.v1). This software contains a fully recursive, vectorized version of NNOMP, NNOLS and SNNOLS together with test programs.

V. SIMULATED SPARSE DECONVOLUTION A. Data generation

In order to assess the behavior of NNOG algorithms, we consider a convolution problem with a Gaussian kernel h of standard deviation σ. h is approximated by a finite impulse response of length 6σ by thresholding the smallest values. The corresponding normalized dictionary H is a Toeplitz matrix. Simulated data are generated according to y = Hx * + n where x * and n stand for the ground truth and white Gaussian noise, respectively. The support S * of x * , of cardinality K, is randomly generated with a uniform distribution whereas the non-zero coefficients of x * are either set to a positive constant or randomly generated from an i.i.d. folded normal distribution. The signal-to-noise ratio is defined by SNR = 10 log 10 (P Hx * /P n ) where P Hx * = Hx * 2 /m is the average power of the noise-free data and P n is the noise variance. The results presented below are achieved on a macOS X system with 16 GB RAM and Intel Core i7 processor at 2.7 GHz.

Hereafter, we first evaluate the computational cost of NNOG algorithms with the proposed implementations. In subsections V-B and V-C, the computational complexity of NNOMP, SNNOLS and NNOLS is compared with that of OMP and OLS. Subsection V-D contains a comparison of NNOG algorithms with other non-negative sparse algorithms proposed in the literature including FNNOMP, FNNOLS [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF], [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF], NN-CoSaMP, and NN-SP [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]. Algorithms are evaluated in terms of sparse reconstruction accuracy and reconstruction time.

B. Validation of NNOG accelerations

We generate a convolution dictionary of size 1200 × 1140 with σ = 10. For each K ∈ {20, 40, 60, 80}, 200 trials are carried out in which the support S * is generated. The non-zero coefficients of x * are set to 1 and the SNR is set to 30 dB. OMP, OLS, NNOMP, SNNOLS and NNOLS are run until the support cardinality equals K. We consider the recursive implementation of OMP and OLS utilizing the matrix inversion lemma. We chose this implementation to make it consistent with our proposed NNOG implementation. Note that NNOG algorithms may need more than K iterations because of support compression. The following quantities are computed and averaged out over the number of trials:

• Acceleration gain: CPU time ratio between the canonical and accelerated implementations of NNOG algorithms. • Non-negativity loss: CPU time ratio between accelerated implementation of NNOG algorithms and the corresponding unconstrained versions. • Iterations: average number of iterations of NNOG algorithms needed to yield a support of cardinality K. This number is larger than K when support compression occurs. Before going further, let us clarify that the so-called "canonical implementations" refer to the following settings. The AS-NNLS algorithm is called from scratch from the initial zero solution. Moreover, obvious accelerations are taken into account such as type-I pruning, computation of the ULS solutions for candidate supports S ∪ {i} so as to avoid calling AS-NNLS when the augmented support is positive, and recursive computation of the AS-NNLS iterates. On the other hand, advanced accelerations such as type-II pruning and warm start initialization are not included. Support compression is included in such a way that both canonical and accelerated versions of a given NNOG algorithm yield the same iterates.

The scores can be found in Tab. I. Accelerated implementations yield a gain in time that increases with K. Since NNOLS needs to solve many NNLS subproblems per iteration, the gain is much larger. The time gain is intermediate for NNOMP. We further notice that using accelerated implementations, the cost of NNOMP and SNNOLS becomes comparable with that of OMP and OLS, respectively, the non-negativity loss remaining below 1.3. Regarding NNOLS vs OLS, the non-negativity loss remains lower than 5 in these simulations. At last, the fact that the number of iterations is often larger than K reveals that support compression is happening quite often.

C. Computation burden of SNNOLS and NNOLS

The purpose of this subsection is to examine thoroughly the computation burden of SNNOLS and NNOLS so as to evaluate the accelerations proposed in Section IV. Four indicators are computed at each iteration: sel applies to the selected atoms, which brings information on the computational burden of the coefficient update stage of SNNOLS.

1) ρ ↓ = |D S |/
Using the same setting as in the previous subsection (H ∈ R 1200×1140 and SNR = 30 dB) with K = 80, the computational burden indicators of SNNOLS and NNOLS are shown in Fig. 2. It is noticeable that the number of iterations L required to reach a support of cardinality K is larger than K because of support compression. Specifically, the histograms of Fig. 2(b) show that on average, L is larger than K for NNOLS and even larger for SNNOLS (the average values are given in the last columns of Tab. I). This is consistent with the fact that the NNOLS selection rule is more involved but more reliable. Moreover, the standard deviation of L corresponding to the histograms of Fig. 2(b) is 17 and 12 for SNNOLS and NNOLS, respectively, which indicates that the size of the support found after k iterations may significantly vary between trials. In order to get meaningful evaluations, we choose to compute the average values of each indicator over the last t iterations, with t ∈ {0, . . . , K -1}. When t = 0, only the last iteration is taken into account, so the current support is of size K. For larger values of t, the supports found during the last t iterations have varying sizes but the averaging operation remains meaningful, especially for the last iterations, which are the most costly. The curves displaying the average of each indicator over 200 trials and over the last t iterations are shown in Fig. 2(a). One can observe from the curve (1-ρ ↓ ) that the rate of non-descending atoms gradually increases for both SNNOLS and NNOLS. Therefore, type-I pruning is more effective at late iterations, where half of the atoms are discarded.

The efficiency of NNOLS implementation is measured by ρ + cand and ρ II in Fig. 2(a). Large values of ρ + cand and ρ II are obtained, which indicates that the computing cost of NNOLS is dramatically reduced. The score of ρ + cand always remains above 70%. This means that the computation of the NNLS solutions reduces to a single ULS update for more than 70% of candidate atoms. Moreover, many of the remaining atoms can be discarded due to type-II pruning. For instance, in the last iterations, around 25% of descending atoms do not yield a positive support and type-II pruning eliminates 90% of them. As a result, only 3% of descending atoms require to proceed with a complete NNLS solving. Let us stress that the extra cost of NNOLS as compared with SNNOLS essentially comes from the number of atoms related to a complete NNLS solving, so it directly depends on the efficiency of type-II pruning.

Regarding SNNOLS, the ratio ρ + sel remains above 20%, which indicates that the computation of the SNNOLS iterate reduces to a single ULS update for at least 20% of the trials. Generally speaking, the slight decrease of ρ + sel , ρ + cand and ρ II shows that both SNNOLS and NNOLS call AS-NNLS more often at the late iterations. A possible reason for this behavior might be that the average correlation between selected atoms increases with the cardinality of the support.

D. Performance comparison between NNOG and competing algorithms

NNOG algorithms NNOMP, SNNOLS and NNOLS, integrating support compression and the proposed acceleration strategies, are compared with the corresponding unconstrained versions OMP and OLS, where a post-processing is applied to obtain non-negative coefficients. NNOG algorithms are also compared to the fast versions FNNOMP and FNNOLS [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF], [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF], as well as the non-negative extensions of CoSaMP and Subspace Pursuit proposed in [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]. All comparisons are based on reconstruction accuracy and CPU time. Reconstruction accuracy is quantified by three factors:

• Support recovery: ratio of true positives to K; • Coefficient inaccuracy: relative error for the recovered coefficients ( xx * / x * ); • Residual norm: Euclidean norm of the data residual ( y -Hx ). The curves shown in Fig. 3 are obtained for a convolution dictionary of size 2400 × 2388, with σ = 2 and for various settings of K ∈ [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF]80]. For each K, 400 trials are carried out in which the support S * is drawn according to the uniform distribution and the non-zero coefficients of x * are drawn from a Gamma distribution with shape and scale parameters set to 1 and 2, respectively. Note that for increasing K, the density of spikes increases, and hence the difficulty of the problem. Working at a given level of SNR would have even increased the difficulty at large values of K, so we have preferred to keep the noise level P n constant; it is set to P n = 10 -2 . In this simulation condition, it is noticeable in Fig. 3(b,f) that the coefficient inaccuracy x-x becomes somewhat large when the spike density is high. For easier problems (lower value of σ and lower density of spikes), many algorithms would provide good results, so it would be difficult to observe differences in terms of estimation quality.

1) Comparison with unconstrained greedy algorithms:

We first perform a comparison of the NNOG algorithms with (unconstrained) OMP and OLS to demonstrate the usefulness of working with NNOG algorithms. In order to obtain non-negative coefficients from the OMP (respectively, OLS) outputs, a simple possibility is to apply a post-processing step using the solution support yielded after K iterations of OMP (resp., OLS). This post-processing consists of solving the NNLS problem related to the support found by OMP/OLS, with possible annealing of coefficients when non-negativity constraints are activated. The resulting schemes are denoted by OMP+ and OLS+.

Fig. 3(a-d) includes a comparison of NNOMP and NNOLS with OMP, OLS, OMP+ and OLS+. Their support recovery and coefficient accuracy performances are weaker than those of NNOG algorithms. In Fig. 3(a), the scores of OMP+ (OLS+) are identical to those of OMP (OLS). This is because the support recovery measure considers the rate of true positives, which is unchanged since the post-processing step of OMP+ (OLS+) essentially removes false positives. On the contrary, the coefficient inaccuracy ratio is improved due to the latter removal. In Fig. 3(c-d), the time and error scores of OMP (resp., OLS) are both lower than those of NNOMP (resp., SNNOLS and NNOLS). This is not a surprise since unconstrained algorithms are simpler, and the obtained solutions are expected to reach lower values of the residual since they do not satisfy the non-negativity constraint. The simple post-processing in OMP+/OLS+ does not allow to yield a residual norm as small as the one obtained using NNOG. This shows that NNOG algorithms do improve the performance of (unconstrained) OMP and OLS.

2) Comparison with approximate fast non-negative greedy algorithms: Fig. 3(e-h) includes a comparison with the FNNOMP and FNNOLS implementations 3 of Yaghoobi et al. [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF], [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF]. FNNOLS turns out to be much slower and less accurate than our fast implementations of SNNOLS and NNOLS. On the other hand, the statistical performances of NNOMP and FNNOMP are very close both in terms of computational time and accuracy. However, FNNOMP may return some negative coefficients in the reconstructed sparse vector. FNNOMP was introduced as a fast but approximate version of the canonical NNOMP algorithm by Bruckstein et al. In contrast, our fast implementation of NNOMP is exact and the non-negativity constraint is always satisfied.

3) Comparison with non-negative versions of CoSaMP and SP: In Fig. 3(e-h), NNOG implementations are also compared with NN-CoSaMP and NN-SP [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]. NN-HTP [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF] yields the weakest performance, so it is not further considered here.

The principle of NN-CoSaMP and NN-SP is to maintain an estimated support of size K. To do so, three kinds of operations are performed at each iteration including support merging, NNLS estimation of coefficients and support pruning. In NN-CoSaMP, the current support of size K is merged with a support of size 2K, and an NNLS problem is solved with this augmented support. Then, thresholding is performed by keeping the K coefficients having the largest magnitudes. The structure of NN-SP is similar, but two NNLS problems (of size roughly 2K and K) are solved before and after thresholding. The NN-SP output thus contains at most K nonzero elements depending on the activation of nonnegativity constraints.

In Fig. 3(e-h), the NNOG algorithms yield competitive accuracy performance compared to NN-SP and NN-CoSaMP. 3 kindly provided by Dr. Mehrdad Yaghoobi.

In particular, NNOG algorithms have a better ability to find a low value of the residual norm for a given cardinality K. This said, no algorithm outperforms all competitors in terms of coefficient accuracy for all scenarios. Nonetheless, the times of computation of NNOG with our recursive implementations are always lower than those of NN-CoSaMP and NN-SP. Let us stress that the structure of the latter algorithms does not easily lend itself to recursive implementations. Indeed, there is no nested property between the supports found at consecutive iterations. Consequently, solving the NNLS problems in a very efficient manner does not seem obvious anymore. For instance, NN-CoSaMP needs to solve NNLS problems for augmented supports of size roughly 3K. The current K-sparse vector obtained after thresholding may be used as a warm start (although it is usually not an NNLS solution), but the number of inner NNLS iterations is expected to be much larger than for NNOG algorithms.

VI. SPARSE DECOMPOSITION OF NIR SPECTRA

A. Real data and generated dictionary

More than 300 wood samples with different compositions (raw wood, plywood, particle boards, MDF-HDF) and finitions (raw, painted, varnished) were collected on a wood waste park and scanned using a Nicolet 8700 FTIR spectrometer. The resulting reflectance spectra are composed of 1647 wavenumbers covering the NIR range 3600-10000 cm -1 (which corresponds to wavelengths in the interval [1, 2.8] µm). The aim is to design a binary classifier based on the identification of a subset of informative wavelengths for detecting the so-called non recyclable samples, i.e., MDF-HDF [START_REF] Belmerhnia | Simultaneous regularized sparse approximation for wood wastes NIR spectra features selection[END_REF], [START_REF] Wagner | Analysis of selected properties of fibreboard panels manufactured from wood and leather using the near infrared spectroscopy[END_REF]. As the final objective is to design an industrial (fast) sorting system, the number of selected wavelengths has to be as small as possible (typically between 16 and 32). Sparse modeling of the (i) NN-SP Fig. 4. Data approximation and sparse recovery of an NIR spectrum using various algorithms. The Gaussian dictionary contains 2998 atoms and the sparsity level is set to K = 20.

NIR spectra is a common approach to select the informative wavelengths. For instance, Turlach et al. [START_REF] Turlach | Simultaneous variable selection[END_REF] proposed an algorithm for variable selection based on an extension of the LASSO and applied it to near infrared measurements.

Hereafter, we consider the decomposition of 50 NIR spectra, seen as data vectors y of length 1647. Data pre-processing includes baseline removal, offset correction ensuring zero lower bound, and unit energy normalization. To decompose the spectra, a dictionary H is built with Gaussian-shaped columns obtained by discretizing the parameters of a Gaussian function (centers and widths). It is formed by appending the columns of the convolution dictionaries (corresponding to a fixed width σ) used in Section V for 60 equally spaced values of σ ∈ [10, 600] cm -1 . The generated dictionary is composed of 2998 atoms. The centers of Gaussian atoms of same width σ are sampled with a step equal to σ, whereas the sampling step of the input signals y equals 4 cm -1 .

B. Decomposition results

For each spectrum, NNOMP, SNNOLS, NNOLS, Sparse NNLS, OMP, OLS are run until a support of cardinality K = 20 is found. The simple OMP+ and OLS+ algorithms introduced in § V-D are also considered, together with FNNOMP, FNNOLS, NN-CoSaMP, and NN-SP. For each algorithm, the CPU time and the relative residual norm y -Hx / y , averaged over 50 spectra, are displayed in Tab. II.

First, the time and approximation errors of OMP (resp., OLS) are both lower than those of NNOMP (resp., SNNOLS and NNOLS). This is not a surprise since unconstrained algorithms are expected to reach lower criterion values. When one performs a post-processing by removing these negative spikes such as with OMP+/OLS+, the error scores become larger than those of NNOG algorithms, and the sparse approximation accuracy is weaker. This confirms the weak capacity of OMP and OLS to reconstruct correct supports from non-negative sparse representations. NNOMP, Sparse NNLS, and FNNOMP have roughly the same cost and performance. Furthermore, SNNOLS and NNOLS yield the lowest approximation errors.

Similar to the simulations of Section V, the computing times of NN-CoSaMP and NN-SP are greater than those of NNOG implementations. NNOLS is significantly slower than NNOMP and SNNOLS, while SNNOLS has the same performance as NNOLS and a computing time closer to that of NNOMP. Finally, the CPU times of NNOG algorithms vs their related unconstrained versions given in Tab. II are consistent with the non-negativity losses gathered in Tab. I.

These approximation and sparse recovery results are further illustrated for a given spectrum in Fig. 4. Note that the Sparse NNLS result (not shown) identifies with that of NNOMP. OMP and OLS (not shown) yield 7 and 9 negative peaks. Besides, the OMP+ and OLS+ approximations around 4500-7000 cm -1 and 5000 cm -1 , respectively, are rather poor. SNNOLS and NNOLS outputs almost coincide. They outperform that of NNOMP, in particular, around 4200 and 7000 cm -1 and those of NN-CoSaMP and NN-SP around 9000 cm -1 . FNNOLS is more costly and slightly less accurate than SNNOLS and NNOLS.

VII. CONCLUSION

Until now, greedy algorithms dedicated to non-negative sparse signal reconstruction have been considered as slow schemes, requiring the repeated resolution of constrained least square problems. In order to accelerate the computation, approximate schemes have been proposed [START_REF] Yaghoobi | Fast non-negative orthogonal matching pursuit[END_REF], [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF] at the price of some loss of control on the algorithmic behavior, and possibly degraded performance. Another commonly found option has been to replace the 0 "norm" by the 1 -norm and to solve the resulting convex programming problem, with a possible loss in terms of performance (see [START_REF] Peharz | Sparse nonnegative matrix factorization with 0 constraints[END_REF] for an interesting case of sparse NMF).

The first contribution of this paper is to provide a unified framework to define non-negative orthogonal greedy algorithms in a formal way, ensuring well-defined iterations. The second and probably most important one in terms of practical impact, is to show that the additional cost of non-negative greedy algorithms to handle the sign constraint can be strongly reduced using three ingredients. The main one is that nonnegative greedy algorithms can be made fully recursive. Moreover, several pruning strategies can be combined to reduce the number of tests at the atom selection stage. Finally, the latter step can benefit from vectorized computations. According to our practical tests with Gaussian dictionaries, the computing time of NNOMP becomes comparable to that of OMP. Hence, we believe that the additional cost should not prevent users from skipping from OMP to NNOMP in nonnegative sparse problem. The performances of NNOMP and FNNOMP are very close both in terms of computational time and accuracy. Since FNNOMP was introduced as a fast but approximate version of NNOMP, it does not seem mandatory to work on approximate schemes to derive efficient implementations. On the other hand, we have obtained a dramatic acceleration of NNOLS compared to both canonical and fast versions proposed in [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF]. The computing cost of NNOLS remains larger than that of OLS. However, whenever OLS can be used, our exact implementation of NNOLS is a realistic option to handle nonnegativity constraints, given that the potential of performance gain between NNOMP and NNOLS is comparable to the one between OMP and OLS in the unconstrained case. We have also studied SNNOLS, which is a suboptimal version of NNOLS originally introduced in [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF]. Our conclusion is that SNNOLS represents a good trade-off between NNOMP and NNOLS: it is structurally simpler than NNOLS and possibly faster, with very similar performance in terms of estimation error. Likewise, we have compared NNOMP with Sparse NNLS and concluded that Sparse NNLS is structurally simpler than NNOMP with almost the same performance as NNOMP, and a slightly reduced computing cost. Let us point out that the computation time and accuracy found for different algorithms might potentially follow different trends with different kinds of dictionaries. The latter conclusions (e.g., comparisons between NNOLS and SNNOLS) found with Gaussian dictionaries could therefore be different when working with other problems. We conjecture though that the performance gain between SNNOLS/NNOLS with respect to NNOMP would be observed as well for other problems involving ill-posed dictionaries, since this conclusion is consistent with the practitioner knowledge of OLS and OMP.

Our comparisons with NN-SP and NN-CoSaMP show that NNOG algorithms are very competitive in terms of accuracy, although no algorithm outperforms the others in all scenarios. The proposed recursive implementations yield a substantial reduction of computation compared to NN-SP and NN-CoSaMP. An attractive feature of NNOG algorithms is their versatility. Indeed, the implementations presented in this paper can be used to address both forms of ( 0 +) problems with constraints x 0 ≤ K (as addressed in this paper) or y -Hx 2 2 ≤ ε 2 , where ε 2 is related to the noise variance. In the latter case, one simply needs to replace the stopping condition |S| = K in Algorithm 1 by r S 2 2 ≤ ε 2 . NN-SP, NN-CoSaMP, and NN-HTP are apparently less versatile, although the proposal of novel versions dedicated to the constraint y -Hx 2 2 ≤ ε 2 could be addressed as a perspective to the work in [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF].

Our contributions can be extended in several directions. A straightforward generalization can be made to deal with nonnegativity-constrained simultaneous sparse decomposition, which is useful in several applications such as hyperspectral imaging [START_REF] Wang | Cone-based joint sparse modelling for hyperspectral image classification[END_REF], dynamic PET [START_REF] Lin | Sparsity constrained mixture modeling for the estimation of kinetic parameters in dynamic PET[END_REF], and diffusion MRI [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]. On the other hand, other greedy algorithms such as BOMP [START_REF] Herzet | Bayesian pursuit algorithms[END_REF] and SBR [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] could also be extended to the non-negative setting using similar principles and using a recursive implementation. Finally, let us mention that [START_REF] Nguyen | Sign preservation analysis of orthogonal greedy algorithms[END_REF] establishes the first K-step recovery analysis of NNOMP, NNOLS and SNNOLS under the Mutual Incoherence Property condition.

APPENDIX A PROOF OF TECHNICAL RESULTS

A. Proof of Proposition 1

Let us first prove that h t i r + S > 0 implies i ∈ D S . Let h t i r + S > 0. According to (1) for r = r + S , we deduce that i ∈ D S since Since f is convex, one has

(z S∪{i} -z S ) t ∇f (z S ) ≤ f (z S∪{i} ) -f (z S ) < 0 (21)
where the gradient of f is defined by ∇f (z S ) = 2H

and thus h t i r + S > 0. Let us now assume that S is a positive support. According to Lemma 2, we have x+ S = xS and r + S = r S , so [START_REF] Sturm | Comparison of orthogonal matching pursuit implementations[END_REF] and ( 8) are identical, as well as ( 9) and [START_REF] Virtanen | Active-set Newton algorithm for overcomplete non-negative representations of audio[END_REF]. To show that ( 7)-( 8) are equivalent to ( 9)- [START_REF] Virtanen | Active-set Newton algorithm for overcomplete non-negative representations of audio[END_REF], we first notice that r S ∈ (span(H S )) ⊥ , thus ht i r S = h t i r S since hih i ∈ span(H S ). Therefore, (8) rereads 0 < ht i r S , which implies that hi = 0 and gi = hi / hi = 0. Hence, (8) rereads 0 < gt i r S , which identifies with [START_REF] Virtanen | Active-set Newton algorithm for overcomplete non-negative representations of audio[END_REF].

Finally, let us show that i ∈ D S is equivalent to condition [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF] x(V ) ← x(V ) -x( ) θ j /θ j (j) ; 13 Θ ← Θ(-j, -j) -θ j (-j)θ j (-j) t /θ j (j) ; Any iteration of AS-NNLS starts with the addition of a new atom to the current support (Line 4 of Algorithm 2). Then, a variable number of atoms are removed from it one after the other (Lines 6 to 12). Let r denote the number of removals at the current AS-NNLS iteration. Let also V ⊂ {1, . . . , n} and xV respectively stand for the current support and solution obtained at Lines 4 and 5, and let V ⊂ V and xV denote the corresponding quantities after r removals. Let us show that any AS-NNLS iteration for which r = 0 or r = 1 yields a solution of the NNLS problem restricted to the support V .

If xV ≥ 0, then r = 0, so V = V and xV = x+ V ≥ 0. Otherwise, we have min(x V ) < 0 and r > 0. If r = 1, a single index -is removed at Lines 6-12, so that V = V \{ -}, and xV ≥ 0. Let us remark that we have xV ( -) < 0 according to Line 7. Let us then prove that xV = x+ V by showing that KKT conditions are satisfied at xV for the NNLS problem related to support V . Note that the NNLS solution is unique since the supports generated by AS-NNLS are such that H V is full column rank, as pointed out in subsection III-B. According to Lemma 1, the KKT conditions read:

H t V (y -H xV ) = 0, (25a) h t -(y -H xV ) ≤ 0. (25b)
(25a) is obviously satisfied. On the other hand, remark that xV = x+ V since xV ≥ 0 and that according to Proposition 1, xV ( -) < 0 implies that h t -r V ≤ 0, which identifies with (25b). This concludes the proof.

APPENDIX B RECURSIVE IMPLEMENTATION OF ULS

Algorithm 4 recalls the recursive ULS computation using MIL [START_REF] Björck | Numerical Methods for Least Squares Problems[END_REF] for selection (forward move V ← V ∪ { }) and deselection (backward move V ← V \{ }) operations. The ULS solution x := xV is updated. Moreover, Θ := (H t V H V ) -1

refers to the inverse of the Gram matrix related to subset V . The Boolean entry fw is set to true and false for selection and deselection updates, respectively. Finally, e 2 stands for the squared residual error r V 2 . All these factors are updated in Algorithm 4. Notation -j refers to all indices except j, and θ j stands for the j-th column of Θ.

The calls to Algorithm 4 for updating ULS solutions at Lines 4-5 and 10-11 of Algorithm 2 take the respective forms: ULS(y, H, V, 1, + , xV , Θ, r V 2 ), ULS(y, H, V, 0, -, xV , Θ, r V 2 ).

Lemma 1

 1 Consider the NNLS problem related to support S: min x≥0 y -Hx 2 s.t. supp(x) ⊂ S.

Fig. 1 .

 1 Fig. 1. Step-by-step illustration of NNOG after each change of support. Bullets represent supports corresponding to the first seven NNOG iterates whereas other intermediate supports found during the calls to AS-NNLS are represented without bullets.

Fig. 1 4 Call 5 S ← C ; 6 r

 1456 Fig.1displays a schematic step-by-step illustration of NNOG. NNOG iterates S k are represented with bullets. NNOG starts with the empty support. In this example, it turns out that during the first three iterations, NNOG yields a positive support S k-1 ∪ { } (x S k-1 ∪{ } ≥ 0), hence S k ←

Fig. 2 .

 2 Fig. 2. Complexity analysis of SNNOLS and NNOLS for a sparse deconvolution problem with Gaussian kernel (σ = 10) and a SNR of 30 dB. The dictionary size is 1200 × 1140, K = 80 and 200 trials are performed. (a) Evolution of complexity factors during the last t iterations (from L -t to L). (1 -ρ ↓ ): atoms discarded by type-I pruning; ρ + sel : positive supports found by SNNOLS; ρ + cand : candidate atoms yielding a positive support in NNOLS; ρ II : atoms discarded by type-II pruning in NNOLS. (b) Histogram of the average number of iterations L required to reach a support of cardinality K = 80.

Fig. 3 .

 3 Fig. 3. Comparison of NNOG algorithms with unconstrained greedy algorithms and other non-negative sparse solvers. Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, Pn = 10 -2 , and a dictionary of size 2400 × 2388. (a)-(d) NNOG vs unconstrained greedy algorithms; (e)-(h) NNOG vs other sparse non-negative algorithms.

2 .

 2 Conversely, let i ∈ D S . Define f (z) = y -H S∪{i} z 2 where z ∈ R |S|+1 . Let us also define the subvectors z S := x+ S (S ∪ {i}) and z S∪{i} := x+ S∪{i} (S ∪ {i}), x+ S and x+ S∪{i} being two NNLS solutions related to S and S ∪ {i}. Condition i ∈ D S reads f (z S∪{i} ) < f (z S ).

12

 12

14 V

 14 ← V \ { } ; 15 end B. Proof of Proposition 3

TABLE I ACCELERATION

 I GAIN OF NNOG ALGORITHMS FOR A SPARSE DECONVOLUTION PROBLEM WITH GAUSSIAN KERNEL (σ = 10). MEAN OVER 200 TRIALS. THE DICTIONARY SIZE IS 1200 × 1140.

	K	Acceleration gain NNOMP SNNOLS NNOLS	Non-negativity loss NNOMP SNNOLS NNOLS	Iterations NNOMP SNNOLS	NNOLS
	20	2.7	1.2	9	1.2	1.0	1.5	20	22	21
	40	4.5	2.1	59	1.2	1.2	1.9	41	51	43
	60	6.4	3.4	120	1.1	1.2	2.1	65	97	73
	80	8.5	4.0	128	1.1	1.3	4.1	95	152	121

  The computational cost of an NNOG iteration is closely related to the values of these ratios. Indeed, large scores indicate that the cost of testing candidate atoms is dramatically reduced. Specifically, ρ ↓ is the rate of candidates that are truly considered in the selection rule[START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] of NNOG algorithms. Both ratios ρ + cand and ρ II quantify the computational burden of the NNOLS selection step: large values of ρ + cand and ρ II indicate that AS-NNLS has to be run for a few candidate atoms only, since other atoms are either pruned or yield a non-negative ULS solution. ρ + sel is defined similar to ρ + cand . However, ρ +

	4) ρ + sel : rate of selected atoms yielding a positive support
	S ∪ { }.

| S|: rate of descending atoms. The rate of discarded atoms after type-I pruning reads 1 -ρ ↓ . 2) ρ + cand : rate of descending candidate atoms for which S ∪ {i} is a positive support. 3) ρ II : rate of candidate atoms discarded by type-II pruning among all atoms i for which S ∪ {i} is not a positive support.

TABLE II CPU

 II TIME AND NORMALIZED APPROXIMATION ERROR OF GREEDY ALGORITHMS FOR SPARSE DECOMPOSITION OF NIR SPECTRA. AVERAGE OVER 50 SPECTRA. SYMBOL ** INDICATES THAT THE CONSIDERED ALGORITHM DOES NOT ENFORCE THE NON-NEGATIVE CONSTRAINT.

	Algorithm			**OMP OMP+ NNOMP FNNOMP Sparse NNLS **OLS OLS+ SNNOLS NNOLS FNNOLS NN-CoSaMP NN-SP
	Time (ms)			23		24	28	25		25		25	26	30		40	2720		45		57
	y -Hx / y × 10 2	6.6		16.4	9.6	9.6		9.6		5.7	16.6	6.4		6.4	7.0		10.1		8.2
	0.06						0.06							0.06
	0.05						0.05							0.05
	0.04						0.04							0.04
	0.03						0.03							0.03
	0.02						0.02							0.02
	0.01						0.01							0.01
	0						0							0
	4000	5000	6000	7000	8000	9000	10000	4000	5000	6000	7000	8000	9000	10000	4000	5000	6000	7000	8000	9000	10000
			(a) OMP+						(b) NNOMP						(c) FNNOMP
	0.06						0.06							0.06
							0.05							0.05
	0.04						0.04							0.04
	0.03						0.03							0.03
	0.02						0.02							0.02
	0.01						0.01							0.01
	0						0							0
	4000	5000	6000	7000	8000	9000	10000	4000	5000	6000	7000	8000	9000	10000	4000	5000	6000	7000	8000	9000	10000
			(d) OLS+						(e) SNNOLS						(f) NNOLS
	0.06						0.06							0.06
	0.05						0.05							0.05
	0.04						0.04							0.04
	0.03						0.03							0.03
	0.02						0.02							0.02
	0.01						0.01							0.01
	0						0							0
	4000	5000	6000	7000	8000	9000	10000	4000	5000	6000	7000	8000	9000	10000	4000	5000	6000	7000	8000	9000	10000
			(g) FNNOLS						(h) NN-CoSaMP

  t S∪{i} (H S∪{i} z S -y) = -2H t S∪{i} r + S . Denoting by C := supp(x + S ) the compressed support, we have from Lemma 1 that H t S\C r + S ≤ 0 and H t C r + S = 0. Since x+ S is supported by C, the latter equality implies that z t S ∇f (z S ) = 0. (21) yields (z S∪{i} ) t H t S∪{i} r + S > 0, i.e.,

	(x + S∪{i} (S ∪ {i})) t H t S∪{i} r + S > 0.	(22)
	Since H t C r + S = 0, (22) rereads:
	(x + S∪{i} (S\C)) t H t S\C r + S + (h t i r + S ) x+ S∪{i} (i) > 0 (23)
	and since H t S\C r + S ≤ 0 and x+ S∪{i} ≥ 0, (23) implies that
	(h t i r + S ) x+

S∪{i} (i) > 0

  . Consider the function f (z) = y -H S∪{i} z 2 and the notations z S and z S∪{i} defined above. Assuming that H S∪{i} is full column rank, we have that f is strictly convex, so f admits a unique minimizer z S∪{i} . If x+ Algorithm 4: Recursive ULS[START_REF] Björck | Numerical Methods for Least Squares Problems[END_REF]. Format: ULS(y, H, V, fw, , x, Θ, e 2 ) e 2 ← e 2 + (x( )) 2 /θ j (j) ;

	1 if fw then
	2	φ ← H t V h ;
	6	x(V ∪ { }) ← x(V ∪ { }) + δβ	Θφ -1	;
	7	Θ ←	Θ 0 0 0	+ δ	Θφ -1	Θφ -1	t	;
	8	V ← V ∪ { } ;
	9 else
	10	j ← index of in V ;

S∪{i} (i) > 0, then z S∪{i} = z S and

r + S∪{i} 2 = f (z S∪{i} ) < r + S 2 = f (z S ), that is, i ∈ D S . Conversely, x+ S∪{i} (i) = 0 implies x+ S∪{i} = x+ S , hence r + S∪{i} = r + S and i / ∈ D S . 3 δ ← (1 -φ t Θφ) -1 ; 4 β ← φ t x(V ) -h t y ; 5

e 2 ← e 2 -δβ 2 ;

11

However, this kind of operation was introduced in[START_REF] Leichner | A strictly improving linear programming Phase I algorithm[END_REF],[START_REF] M∅rup | Approximate L0 constrained non-negative matrix and tensor factorization[END_REF].

Contrary to screening techniques, see e.g.,[START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and grouplasso[END_REF], note that the atoms indexed by i / ∈ D S are pruned from the dictionary for the current NNOG iteration only, but they are considered again in further iterations.

ACKNOWLEDGMENT

The authors would like to thank Dr. Mehrdad Yaghoobi (University of Edinburgh) for kindly sharing his Matlab implementation of FNNOLS.

This work was supported by the project ANR 15-CE23-0021 BECOSE. It was carried out in part while T. T. Nguyen was visiting L2S during the academic years 2017-2019.