Keywords: orthogonal greedy algorithms, sparse reconstruction, non-negativity, nonnegative least squares, active-set algorithms, exact recovery condition, mutual coherence, NPhardness, 0 minimization

Algorithmes gloutons orthogonaux sous contrainte de positivité

Reconstruction de signaux parcimonieux positifs

De nombreux domaines applicatifs conduisent à résoudre des problèmes inverses où le signal ou l'image à reconstruire sont à la fois parcimonieux et positifs. Citons les exemples de la télédétection [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF], de l'analyse de signaux audio [START_REF] Virtanen | Active-set Newton algorithm for overcomplete non-negative representations of audio[END_REF], de la chimiométrie [Bro and De Jond, 1997], de l'astrophysique [Högbom, 1974] ou de la reconstruction tomographique [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF]. La reconstruction d'un signal parcimonieux et positif est naturellement formulée comme un problème de minimisation 0 sous contrainte de positivité, du type:

min x y -Hx 2 s.c. x ≥ 0, x 0 ≤ K, (1)
où . et . 0 représentent respectivement la norme euclidienne et la "norme" 0 , correspondant au nombre d'éléments non-nuls d'un vecteur. Ce problème est connu pour être NP-difficile [Natarajan, 1995]. On distingue plusieurs catégories d'approches allant des méthodes de relaxation convexe (via les algorithmes proximaux notamment) aux méthodes d'optimisation non-convexe d'un critère continu et non-différentiable approchant [START_REF] Gasso | Recovering sparse signals with a certain family of nonconvex penalties and DC programming[END_REF] ou reformulant [START_REF] Soubies | A continuous exact l0 penalty (cel0) for least squares regularized problem[END_REF] le critère (1), et aux méthodes abordant directement l'optimisation du critère 0 non-modifié via une stratégie d'optimisation exacte [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF] ou sous-optimale. Les algorithmes gloutons font partie de cette dernière catégorie. Leur principe est de partir d'une solution nulle et de construire progressivement une solution K-parcimonieuse en sélectionnant un à un des atomes du dictionnaire H, avec une mise à jour des amplitudes associées et du résidu d'approximation r = y -Hx au fur et à mesure que les nouveaux atomes sont sélectionnés.

Les algorithmes gloutons orthogonaux comme Orthogonal Matching Pursuit (OMP) et Orthogonal Least Squares (OLS) mettent à jour les amplitudes par un calcul de projection orthogonale. En notant S le support courant, les amplitudes associées x(S) sont celles pour lesquelles le résidu d'approximation r S = y -H S x(S) est minimum en norme 2 , où H S est la sousmatrice de H réduite à S. Leur calcul nécessite la résolution d'un problème de moindres carrés non-contraints (ULS pour Unconstrained Least Squares). Les implémentations efficaces des algorithmes gloutons reposent sur l'idée que la solution de ce problème de moindres carrés peut être mise à jour rapidement lorsque le support est augmenté d'un élément (S ← S ∪ { }), modulo la mise à jour d'une factorisation matricielle de type Gram-Schmidt ou de Cholesky [START_REF] Sturm | Comparison of orthogonal matching pursuit implementations[END_REF].

Les algorithmes gloutons généralisés à la reconstruction de signaux positifs possèdent une structure proche des algorithmes gloutons standard, avec à chaque itération, la sélection d'un nouvel atome et la mise à jour des amplitudes en résolvant un problème de moindres carrés sous contrainte de positivité (NNLS pour Non-Negative Least Squares) [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF]. Cependant, leur implémentation pose une difficulté majeure liée au fait que les problèmes NNLS ne possèdent pas de solution explicite. Dans la littérature, ces algorithmes sont réputés lents et des solutions récursives approchées ont été privilégiées [Yaghoobi et al., 2015]. Dans [Nguyen et al., 2019a], nous montrons que des implémentations récursives exactes sont possibles sans nécessiter de surcoût calculatoire substantiel. Cette contribution est présentée dans les chapitres 2 et 3.

L'analyse théorique des algorithmes de reconstruction parcimonieuse vise à caractériser leur capacité à reconstruire le support d'une représentation K-parcimonieuse. Les 8 fin d'algorithmes gloutons en K-itérations cherchent à prouver que chaque itération sélectionne un atome dans le support du vecteur parcimonieux inconnu [Tropp, 2004]. Elles ne s'étendent pas directement aux algorithmes gloutons non-négatifs, car les propriétés géométriques sur lesquelles elles s'appuient (en particulier, l'orthogonalité entre le résidu et l'espace d'approximation) ne sont plus valables. Il est donc nécessaire de repenser ces analyses. Dans [START_REF] Nguyen | Sign preservation analysis of orthogonal greedy algorithms[END_REF], nous proposons une analyse basée sur l'hypothèse de cohérence mutuelle faible µ < 1/(2K -1) pour analyser de façon unifiée différents algorithmes proposés dans la littérature [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF], Yaghoobi et al., 2015]. Les mécanismes sur lesquels cette analyse repose sont synthétisés dans le chapitre 4.

Algorithmes gloutons orthogonaux positifs

La structure générale des algorithmes gloutons orthogonaux positifs (NNOG pour Non-Negative Orthogonal Greedy algorithms) est résumée dans l'algorithme 1. Une itération de l'algorithme consiste en trois étapes:

1. Sélection d'un nouvel atome / ∈ S pour enrichir le support S de la représentation parcimonieuse.

2. Mise à jour des amplitudes liées aux atomes sélectionnés (i ∈ S) via la résolution du problème NNLS associé au support de x:

x+ S = arg min x≥0 y -Hx 2 s.c. supp(x) ⊂ S. (2)
3. Compression du support.

La principale différence structurelle avec les algorithmes gloutons orthogonaux classiques est la troisième étape de compression du support. On remarque en effet que les coefficients de la représentation parcimonieuse obtenus après l'étape NNLS s'annulent lorsque les contraintes de positivité du problème (2) sont activées. L'étape de compression que nous avons proposée dans [START_REF] Nguyen | An optimized version of non-negative omp[END_REF], Nguyen et al., 2019a] vise à supprimer les indices correspondants aux coefficients annulés, pour faire coïncider S avec le support de la solution courante x. La structure de l'algorithme est donc sensiblement différente des algorithmes gloutons classiques, qui possèdent une propriété d'emboîtement entre supports. L'algorithme 1 est conçu comme un algorithme de descente. Ses deux conditions d'arrêt garantissent que le vecteur obtenu est exactement K-parcimonieux (|S| = K) ou qu'aucune Contents xiii décroissance du critère y-Hx 2 n'est possible en ajoutant un nouvel atome au support courant. La condition h t i r S > 0 garantit en effet que l'erreur quadratique peut décroître en sélectionnant l'atome i, ce qui nous a conduit à introduire le concept de "atome descendant" dans le chapitre 2.

La structure des algorithmes NNOG dans le tableau 1 est commune à différents algorithmes comme Non-Negative Orthogonal Matching Pursuit (NNOMP) [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF], défini comme une généralisation d'OMP, Non-Negative Orthogonal Least Squares (NNOLS) et Suboptimal Non-Negative Orthogonal Least Squares (SNNOLS), qui sont deux généralisations d'OLS [Yaghoobi and Davies, 2015]. Dans leur principe, ces algorithmes diffèrent uniquement par leur règle de sélection du nouvel atome:

nnomp ∈ arg max i / ∈S h t i r S , (3)
snnols ∈ arg max i / ∈S b t i r S , (4)
nnols ∈ arg min

i / ∈S r + S∪{i} 2 , (5)
où r + S∪{i} := y -H x+ S∪{i} désigne le "résidu positif" lié à la solution NNLS pour le support S ∪ {i}, et b i représente le projeté orthogonal normalisé de h i sur le sous-espace orthogonal aux atomes h j , j ∈ S.

La règle de sélection de NNOMP est une extension directe de celle d'OMP, qui maximise |h t i r S |. De même, NNOLS suit le mécanisme d'OLS, où l'atome sélectionné est celui qui engendre une décroissance maximale du résidu d'approximation. Finalement, SNNOLS est inspiré d'une vision géométrique d'OLS en terme d'angle maximal entre le projeté orthogonal b i de l'atome candidat et le résidu courant. Si les deux interprétations d'OLS en terme d'optimisation et de géométrie projective sont équivalentes, cette équivalence n'est plus valable dans le cas des algorithmes gloutons positifs, ce qui conduit à deux algorithmes distincts SNNOLS et NNOLS (voir chapitre 2).

NNOMP est clairement le moins coûteux des trois algorithmes, car sa règle de sélection ne repose que sur le calcul d'un produit scalaire par atome. Néanmoins, ses limites de performance pour des problèmes mal posés nous ont poussé à nous intéresser aux algorithmes inspirés d'OLS, dont la complexité calculatoire est accrue. La différence de complexité entre NNOLS et SNNOLS s'avère délicate à analyser à partir des expressions (4) et (5). Nous avons montré que (4) équivaut à (voir chapitre 2):

snnols ∈ arg min i / ∈S min u,v≥0 y -H S u -h i v 2 . ( 6
)
Il apparaît alors que le test de sélection d'atome de SNNOLS est largement moins coûteux que celui de NNOLS, car (6) est un problème de moindres carrés non-contraint par rapport à u, et sous contrainte de positivité par rapport à la variable scalaire v. Par opposition, le calcul de r + S∪{i} dans le cas de NNOLS nécessite de résoudre un problème NNLS avec des contraintes de positivité sur u et v.

Implémentation rapide

Dans notre article de conférence [START_REF] Nguyen | An optimized version of non-negative omp[END_REF], nous avons présenté une implémentation rapide de NNOMP basée sur la résolution récursive des sous-problèmes NNLS par l'algorithme des contraintes actives [START_REF] Lawson | Solving least squares problems[END_REF] avec un démarrage à chaud correspondant à une initialisation par l'itéré courant de NNOMP. L'algorithme des contraintes actives résout le problème NNLS de façon exacte en un nombre fini d'itérations. Comme il possède une structure gloutonne, les implémentations qui résultent de son couplage avec les algorithmes NNOG sont xiv Contents pleinement récursives. Nous avons constaté empiriquement que l'algorithme des contraintes actives avec démarrage à chaud converge en un très faible nombre d'itérations.

Dans le chapitre 2, nous avons proposé une série d'améliorations supplémentaires visant à résoudre les problèmes de sélection d'atomes (3)-(5) de manière rapide. Il s'agit (i) de prédétecter que certains atomes candidats ne sont pas descendants et donc les écarter; (ii) pour les atomes restants et dans le cas de NNOLS (qui requiert la résolution répétée de problèmes NNLS), un précalcul des solutions du problème ULS pour le support candidat S ∪ {i} permet de détecter des atomes non-optimaux au sens de (5), ce qui nous évite des appels à NNLS.

À titre d'exemple, les tests effectués dans le chapitre 3 pour un problème de déconvolution impulsionnelle de taille x ∈ R 1140 avec un filtre passe-bas montrent que pour un support S de taille 80, il y a 30 % d'atomes candidats (i / ∈ S) non testés, 68 % d'atomes candidats pour lesquels seule la solution ULS associée à S ∪ {i} est calculée, et seulement 2% des atomes candidats pour lequel un appel à NNLS est requis.

Analyse d'algorithmes gloutons orthogonaux positifs

L'analyse théorique des algorithmes NNOG en termes de reconstruction exacte du support d'une représentation parcimonieuse y = Ax (avec x ≥ 0) est un problème ouvert. En effet, à notre connaissance, il n'existe pas d'analyse mathématique des algorithmes NNOMP, NNOLS et SNNOLS. Néanmoins, Bruckstein et al. [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] ont conjecturé que NNOMP présente des garanties de reconstruction exacte en K itérations si la cohérence mutuelle du dictionnaire µ = max i =j | h i , h j | vérifie:

µ < 1 2K -1 . ( 7
)
Dans le chapitre 4, nous avons prouvé cette conjecture non seulement pour NNOMP mais aussi pour les algorithmes NNOLS et SNNOLS. Notre technique de preuve s'appuie sur une analyse de préservation de signe par les algorithmes OMP et OLS. Nous avons d'abord établi le résultat suivant.

Theorem 1. [Nguyen et al., 2019d, Cor. III.1] Supposons que µ < 1 2K-1 . Soit y = Hx une représentation K-parcimonieuse avec x ≥ 0. Alors OMP et OLS reconstruisent correctement le support de x en K itérations. De plus, pour chaque itération k = 1, . . . , K, les poids des k atomes sélectionnés sont positifs.

La première partie de ce résultat est bien connue. En effet, (7) est une condition nécessaire [START_REF] Cai | Stable recovery of sparse signals using an oracle inequality[END_REF] et suffisante [Tropp, 2004, Herzet et al., 2013] de reconstruction exacte quelles que soient les amplitudes non nulles des poids dans x. La deuxième partie du théorème indique qu'en régime de reconstruction exacte, l'itéré courant d'OMP/OLS a des coordonnées positives. Comme cet itéré correspond à la solution du problème des moindres carrés (ULS) associé au support courant S, la solution NNLS correspondant au même support coïncide avec la solution ULS. La très grande proximité entre les règles de sélection d'atomes d'OMP et de NNOMP (3) nous permet de conclure que les itérés de NNOMP et OMP coïncident. De façon parallèle, ceux de NNOLS, SNNOLS et OLS coïncident, d'où le résultat suivant.

Theorem 2. [Nguyen et al., 2019d, Cor. III.2] Supposons que µ < 1 2K-1 . Soit y = Hx une représentation K-parcimonieuse avec x ≥ 0. Alors les itérés de NNOMP coïncident avec ceux de OMP. De plus, les itérés de NNOLS et SNNOLS coïncident avec ceux de OLS. Par conséquent, NNOMP, SNNOLS et NNOLS reconstruisent tous le support de x en K itérations.

Contents xv

Notons enfin que NNLS étant des problèmes d'optimisation contraints, certaines contraintes de positivité sont susceptibles d'être activées à chaque itération des algorithmes NNOG, entraînant l'annulation des poids correspondants dans la représentation parcimonieuse. Dans ce cas, l'algorithme NNOG doit effectuer un nombre d'itérations supérieur à K pour pouvoir reconstruire une représentation réellement K-parcimonieuse. Le théorème 2 indique que le phénomène d'annulation des poids ne se produit pas en régime de reconstruction exacte, où la reconstruction exacte se produit en K itérations.

Conclusions

Dans ce document, nous présentons un panorama unifié des algorithmes gloutons orthogonaux non-négatifs (NNOG), conçus comme des algorithmes de descente pour le problème de minimisation 0 (1), ainsi qu'une analyse de leur propriété de reconstruction exacte d'un vecteur K-parcimonieux positif sous l'hypothèse de cohérence mutuelle inférieure à 1 2K-1 . Du point de vue pratique, nous avons développé un logiciel Matlab en accès libre qui intègre les implémentations récursives proposées dans la section 2.4 [START_REF] Nguyen | Non-negative greedy algorithms: Matlab implementation[END_REF]. Nous avons effectué des simulations numériques extensives (voire chapitre 3) visant à comparer les algorithmes NNOMP, SNNOLS et NNOLS pour des problèmes difficiles de type déconvolution impulsionnelle avec un filtre passe-bas pour lesquels la cohérence mutuelle est proche de 1. Il ressort d'abord de ces tests que bien que les algorithmes non-négatifs sont plus coûteux que les versions standards OMP et OLS, le surcoût n'est pas rédhibitoire grâce aux implémentations rapides proposées (le ratio en temps de calcul entre les versions standards et non-négatives n'excède jamais un facteur 5 pour des problèmes de taille 1200 avec K = 80). NNOMP s'avère plus rapide mais nettement moins précis que SNNOLS et NNOLS pour reconstruire le vecteur parcimonieux. La précision est évaluée en termes d'erreur quadratique sur x et de nombre de vrais positifs dans le support trouvé. Ces deux algorithmes fournissent des reconstructions parcimonieuses de qualité comparable, SNNOLS étant sensiblement plus rapide. Plus précisément, on peut noter que le nombre d'itérations de NNOLS pour atteindre un support de taille K est sensiblement plus faible que celui de SNNOLS (en moyenne 120 itérations contre 150 sont nécessaires pour obtenir un support de taille K = 80 pour l'exemple donné dans la section 3.3) mais le coût par itération de NNOLS reste plus important. En revanche, les défauts de SNNOLS semblent être plus prononcés pour des problèmes de type "super-résolution" en utilisant des grilles fines (impliquant des dictionnaires plus gros et plus corrélés), où le nombre d'itérations de SNNOLS augmente significativement, induisant un coût calculatoire qui devient globalement supérieur à celui de NNOLS. This work is supported by the BECOSE project (2016-2019)1 aiming to develop efficient sparse algorithms for ill-conditioned inverse problems. The BECOSE project contains three main tasks. The first task is to develop new sparsity-based algorithms dedicated to ill-posed inverse problems by exploiting additional structural information such as non-negativity or by working on continuous dictionaries. The second task is to enrich the theoretical analysis of sparse algorithms, especially the stepwise orthogonal greedy search algorithms, by taking into account prior knowledge on the signs, coefficient values, and partial support information. The third task is to assess the proposed algorithms in the context of tomographic Particle Image Velocimetry (PIV) application. This thesis addresses partly the first and second tasks of the BECOSE project.

Orthogonal greedy algorithms for non-negative sparse reconstruction

In this introduction, we start by presenting the non-negative sparse reconstruction problem and some motivating applications. We next review the state-of-the-art sparse algorithms with a focus on greedy algorithms and recall the state-of-the-art exact recovery analyses of greedy algorithms. Then we summarize our contributions related to non-negative orthogonal greedy algorithms and exact recovery analysis. We end up with the list of resulting publications.

Non-negative sparse reconstruction

Sparse approximation appears in a wide range of applications, especially in signal and image processing and compressive sensing [Elad, 2010]. Given a signal y ∈ R m and a redundant dictionary H ∈ R m×n , one is interested in finding the smallest set of dictionary columns (also called atoms) that describes y well enough. In other words, one seeks for a vector x ∈ R n that gives the best approximation y ≈ Hx and has the fewest non-zero coefficients (i.e., the sparsest solution). A simple illustration is given in Figure 1.1. Sparse approximation often leads to solving an 0 minimization problem that takes one of the following forms:

min x x 0 s.t. y -Hx 2 ≤ ( 0 C ) min x y -Hx 2 2 s.t. x 0 ≤ K ( 0 C K ) min x y -Hx 2 2 + λ x 0 ( 0 P )
in which , K and λ are positive quantities related to the noise standard deviation, the sparsity level and the regularization strength, respectively. Let us recall that the 0 -"norm" counts the number of non-zero coefficients: x 0 = Card{i : x i = 0}. According to its definition, 0 is not a proper norm since it is not homogeneous. However, it is widely used as a pseudo-norm in statistics, scientific computing and information theory. Hereafter, we call Hx a K-sparse representation if x 0 = K. Without loss of generality, we assume that the dictionary H is normalized (so every atom h i has 2 -norm equal to 1). The 2 -norm • 2 will be also denoted • . In several application fields such as geoscience and remote sensing [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF],Esser et al., 2013], audio [START_REF] Virtanen | Active-set Newton algorithm for overcomplete non-negative representations of audio[END_REF], chemometrics [Bro and De Jond, 1997], bioinformatics [START_REF] Slawski | Isotope pattern deconvolution for peptide mass spectrometry by non-negative least squares/least absolute deviation template matching[END_REF], astrophysics [Högbom, 1974] and computed tomography [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF], the signal or image of interest is sparse, but also non-negative. In such contexts, a common practice is to regularize the inverse problem to favor both sparsity and non-negativity of the sought signal, see, e.g., [Hoyer, 2004, Petra and Schnörr, 2014, Slawski et al., 2012, Rapin et al., 2013]. This regularization leads to non-negative 0 minimization problems as follows.

min x≥0 x 0 s.t. y -Hx 2 ≤ ( 0 C +) min x≥0 y -Hx 2 2 s.t. x 0 ≤ K ( 0 C K +) min x≥0 y -Hx 2 2 + λ x 0 ( 0 P +)
To address these problems, a usual approach is to incorporate the non-negativity into existing sparse solvers. The obtained algorithms are called non-negative sparse solvers. However, the extension of sparse solvers to the non-negative setting is not always obvious. Besides, the nonnegative sparse solvers may be significantly more expensive than the primary sparse solvers. We will discuss this later in Section 1.3. Now, let us take a quick tour on some applications of non-negative sparse reconstruction.

Some applications

In this section, we briefly discuss two among several applications in which both sparsity and nonnegativity are required for the sought solution. The presented applications contain tomographic PIV reconstruction and spectral unmixing. More details on other applications can be found in the references cited in the previous section.

Tomographic PIV reconstruction

Tomographic Particle Image Velocimetry (Tomo-PIV) [START_REF] Elsinga | Tomographic particle image velocimetry[END_REF]] is a method to reconstruct 3-D maps of the velocity field of the flows in fluids. The fluid is seeded with point-like tracer particles whose motion is used to calculate the speed and direction of the flow. A laser with an optical arrangement is used to limit the physical region illuminated. Several (at least four) cameras are placed at different angles to record simultaneous views (radiograph, i.e., projection images) of the illuminated particles. The volume is then reconstructed to yield a discretized 3-D intensity field. The principle of Tomo-PIV is summarized in Figure 1.2. Tomo-PIV contains two steps: (i) reconstruction of 3-D volume from 2-D radiographs and (ii) estimation of the vector field from reconstructed volumes at consecutive times. The reconstruction procedure is an under-determined inverse problem since the number of unknowns (the number of voxels of the 3-D volume) is usually much higher than the number of data (the number of radiograph pixels). The light intensity corresponding to a voxel is nonzero if a particle is present at that voxel and zero otherwise. Since the number of tracer particles is small compared with the number of voxels and the light intensity is non-negative, the sought solution corresponding to the volume intensity is sparse and non-negative [START_REF] Barbu | A new approach for volume reconstruction in tomoPIV with the alternating direction method of multipliers[END_REF]. Therefore, Tomo-PIV reconstruction leads to a non-negative sparse approximation problem.

One example of Tomo-PIV reconstruction is given in . The simulated discretized 3-D volume of dimension 79 × 75 × 12 with 81 tracer particles is reconstructed from 4 quantized radiographs of size 32 × 322 . Here the data vector y is built from the 4 images whose pixels are unfolded and rearranged in lexicographical order. The ground truth vector x is built similarly from the unfolded discretized 3-D volume and the amplitudes represent the light intensities at each voxel of the discretized volume. The dictionary H is built from the projection operator [START_REF] Barbu | A new approach for volume reconstruction in tomoPIV with the alternating direction method of multipliers[END_REF], each column is the projection of each voxel. It should be noted that the choice of the 3-D grid is somewhat arbitrary although there are standard conventions in Tomo-PIV. The choice of a fine grid leads to a high resolution but also a more difficult problem Figure 1.2: Principle of Tomo-PIV (source: [START_REF] Elsinga | Tomographic particle image velocimetry[END_REF]) since the number of voxels dramatically increases. It is hence important to know until which resolution one can handle (see Chapter 3 for further discussion). A recent trend is to work with a continuous dictionary (off-the-grid approach) [Ait Tilat et al., 2019,Elvira et al., 2019]. In this thesis, we will not explore the continuous approach. Instead, we will elaborate on the case of a discrete dictionary corresponding to a fine grid.

Figure 1.4 presents our reconstruction from the simulated data in Figure 1.3 using the SBR algorithm [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF]. The rate of correct detection is 68 %. It should be noted that correct detection is the detection of a particle in the same voxel as the ground truth particle. The rate 68 % might seem low but, as can be seen in Figure 1.4, the reconstruction is very accurate in the qualitative viewpoint: each detected particle is located very close to a ground truth particle. Hyperspectral images are often presented as a cube with two spatial dimensions and one spectral dimension. Each hyperspectral pixel is a spectrum collecting light reflectance at different spectral bands. Very often, each spectrum of a hyperspectral image can be seen as a linear mixture of pure spectra related to several materials [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. This happens when different materials are combined in the scene location or when the spatial resolution of the hyperspectral sensor is low enough for many materials to jointly occupy a single pixel. The linear mixture model assumes that the mixed pixel spectrum can be presented as a linear combination of spectra of constituent materials [START_REF] Keshava | Spectral unmixing[END_REF]. Therefore, decomposing the spectrum helps to detect the materials constituting the measured pixel.

Spectral unmixing

Spectral unmixing is the procedure that decomposes a mixed spectrum into different spectra (endmembers) and the corresponding fractions (abundances). This decomposition is often done using a library (dictionary) of spectra [Mortada, 2018]. Since the number of constituting spectra is often very small compared with the number of spectra in the library and the corresponding fractions are non-negative, spectral unmixing leads to a non-negative sparse approximation problem [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF], Wang et al., 2018]. Using our notation, the data vector y is the mixed spectrum, the dictionary H is the library of spectra in which each column presents a pure spectrum. One looks for the vector x containing the fractions of each pure spectrum in the mixed spectrum y. In Figure 1.5, the mixed spectrum (bottom left) is decomposed as a linear combination of 4 spectra (right) extracted from the large size dictionary (not shown) which may contain many hundreds or many thousands of spectra.

State-of-the-art algorithms

The aforementioned 0 minimization problems can be seen as discrete (combinatorial problem) since the main difficulty is to find the support of the sparse representation S = {i : x i = 0}.

The 0 minimization problems are known to be NP-hard for general dictionaries [Natarajan, 1995, Davis et al., 1997]. This means that when the problem dimension is large, it is hopeless to look for the optimal solution since the exhaustive search might be extremely long. Therefore, one rather relies on suboptimal algorithms which run in polynomial time and return a solution (hopefully) close to the optimal. Such algorithms include greedy algorithms, convex and nonconvex solvers. In this section we make an overview of the three categories with a focus on greedy algorithms and their non-negative extensions. We therefore choose to start the section with a review of greedy algorithms, then go towards non-negative greedy algorithms, then quickly present convex solvers and end up with a short discussion about non-convex solvers.

Greedy algorithms

Greedy algorithms are suboptimal algorithms acting on the 0 -norm of solution. Originally, greedy algorithms start from the zero vector (sparsest solution) which corresponds to an empty support, and gradually change the cardinality of the support one-by-one through iterations, hence the name "greedy". According to this formulation, greedy algorithms can be interpreted as descent algorithms (the residual norm y -Hx is decreasing at each iteration) especially suited for reconstructing highly sparse solutions. In the literature, the name "greedy" also refers to another class of algorithms in which a support of fixed cardinality is produced at each iteration. Unlike the previous strategies, these algorithms are not always descent algorithms and they require some prior knowledge about the cardinality of the solution support. We next discuss some typical greedy algorithms of both strategies. These algorithms and the corresponding problem they address are shown in Table 1.1. Note that in the rest of the thesis we use the term "greedy" to refer to the one-by-one updating strategy only.

Problem

Algorithms ( 0 C ) MP, OMP, OLS ( 0 C K ) MP, OMP, OLS, CoSaMP, SP, IHT, HTP ( 0 P ) Bayesian MP, Bayesian OMP, SBR, CSBR Greedy algorithms increasing the support cardinality one-by-one at each iteration contain Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF], 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] 3 , in increasing order of complexity. The greedy mechanism consists of two steps: (i) new atom selection (h ) yielding a new support S ← S ∪ { }, and (ii) update of the sparse coefficients x S . MP and OMP share the way of choosing the new atom by maximizing the inner product of the current residual and dictionary atoms ∈ arg max

i |r t h i |. (1.1)
However, they differ in the way the coefficients are updated. While MP only updates the coefficient related to the new selected atom x , OMP updates all the coefficients x i , i ∈ S by 1. Introduction performing an orthogonal projection (hence the name "orthogonal"). As a result, the residual r = y -H S x(S) is orthogonal to the subspace spanned by selected atoms span (H S ) (see Fig. Least Squares (ULS) subproblem

x(S) ← arg min z y -H S z 2 . (1.2)
MP iteration is very fast but the number of MP iterations to reach a given residual value might be significantly larger than that of OMP. In addition, OMP is empirically more effective than MP [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. OLS updates the coefficients the same way as OMP but differs in atom selection. While OMP selects the atom having the highest inner product (in magnitude) with the current residual, OLS considers the inner product with the projection of atoms onto the orthogonal complement of selected subspace after normalization (the normalized projected atoms). This is illustrated in Algorithms 2-3 in which S denotes the support of the current solution and gS i denotes the normalized projected atom, i.e.,

hS i = P ⊥ S h i , gS i = hS i hS i . (1.3)
Algorithm 2: OMP [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] Input: y, H, K Output:

x solving ( 0 C K ) 1 Initialization: x ← 0; S ← ∅; r ← y ; 2 while |S| < K do 3 ∈ arg max i / ∈S |r t h i | ; 4 S ← S ∪ { } ; 5 x(S) ← H † S y ; 6 r ← y -H S x(S) ; 7 end
Algorithm 3: OLS [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF]]

Input: y, H, K Output: x solving ( 0 C K ) 1 Initialization: x ← 0; S ← ∅; r ← y ; 2 while |S| < K do 3 ∈ arg max i / ∈S |r t gS i | ; 4 S ← S ∪ { } ; 5 x(S) ← H † S y; 6 r ← y -H S x(S); 7 end
Note that from an optimization viewpoint, the selection rule of OMP (line 3 of Algorithm 2) is equivalent to

∈ arg min i / ∈S min z y -H S x(S) -zh i 2 . (1.4) Indeed, min z y -H S x(S) -zh i 2 = min z r -zh i 2 = r 2 -|r t h i | 2 as h i = 1. Besides,
OLS selection rule (line 3 of Algorithm 3) is equivalent to [START_REF] Blumensath | On the difference between Orthogonal Matching Pursuit and Orthogonal Least Squares[END_REF]]

∈ arg min i / ∈S min z y -H S∪{i} z 2 . (1.5)
From this point of view, it is obvious that the OLS selection rule is more costly than that of OMP, although recursive4 update schemes are available, see, e.g., [Chen et al., 1989, Rebollo-Neira and[START_REF] Lowe ; Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF]. In this thesis we are particularly interested in orthogonal greedy algorithms such as OMP and OLS since they are fast in comparison with non-greedy solvers and empirically effective. Some greedy algorithms change the support cardinality one-by-one at each iteration but in both forward and backward directions [Miller, 2002]. This means these algorithms either add or remove one atom at each iteration. Examples include Bayesian MP (BMP), Bayesian OMP (BOMP) [START_REF] Herzet | Bayesian pursuit algorithms[END_REF], Single Best Replacement (SBR) [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] and CSBR [START_REF] Soussen | Homotopy based algorithms for 0 -regularized least-squares[END_REF]. BMP, BOMP, SBR are forward-backward extensions of MP, OMP, OLS, respectively. They are designed as descent algorithms dedicated to the penalized 0 minimization problem ( 0 P ) with a particular value of the parameter λ. CSBR is the "continuous" version of SBR which solves ( 0 P ) for a continuum of values of λ. Like MP, BMP is fast but the three algorithms BOMP, SBR and CSBR are empirically more effective [Herzet andDrémeau, 2010,Soussen et al., 2011]. Like OMP and OLS, the forward-backward extensions BOMP, SBR and CSBR can be efficiently implemented using recursive updating schemes [Sturm andChristensen, 2012, Soussen et al., 2011].

Fixed-cardinality updating strategy

Instead of changing the support size one-by-one, some algorithms produce a support of fixedcardinality at every iteration. Such algorithms are specifically designed to address ( 0 C K ). They include 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], Iterative Hard Thresholding (IHT) [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF], and Hard Thresholding Pursuit (HTP) [Foucart, 2011], to name but a few. CoSaMP and SP are very similar, as well as IHT and HTP. We hence choose to present them closely for ease of comparison (see .

At each CoSaMP iteration, at most 2K atoms (reselection might occur) are selected to enter the support, then the coefficients are updated by performing an orthogonal projection followed by a thresholding. Similarly, SP selects at most K atoms to enter the support and updates the coefficients by performing an orthogonal projection followed by a thresholding and an additional orthogonal projection. Note that CoSaMP performs one orthogonal projection related to a support of size around 3K at each iteration while SP performs two orthogonal projections related to support sizes around 2K and K. In the literature, CoSaMP and SP are known to have very competitive performances [Needell andTropp, 2009, Dai andMilenkovic, 2009].

IHT's principle is to repeat a gradient descent move followed by a thresholding [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF]. Therefore, IHT is very fast at each iteration but it might need a large number of iterations to converge. To improve the convergence speed, more effective step search strategies were proposed later [START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF]. Hereafter, we use the name IHT to refer to the later scheme (presented in Algorithm 6). It should be noted that a version of IHT dedicated to ( 0 P ) was also proposed but it was reported to be less effective than Algorithm 6 [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF].

HTP [Foucart, 2011] is a variant of IHT in which an orthogonal projection is performed to update the coefficients at each iteration. Note that the orthogonal prjection at line 12 of

Introduction

Algorithm 7 is computed approximately using a certain number of gradient descent iterations (otherwise α = 0 at the next iteration). Each iteration of HTP is more costly than that of IHT. However, HTP is empirically more effective [Foucart, 2011].

Since CoSaMP, SP and HTP may change many of the support elements at each iteration, they are difficult to lend themselves to recursive implementations. Besides, their coefficient update is done by orthogonal projection which might be expensive without recursive implementation. Therefore, their coefficient update step is often implemented approximately (see, e.g., [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]). Finally, it is worth emphasizing that, among greedy algorithms performing orthogonal projection at every stage, OMP and OLS can address both ( 0 C ) and ( 0 C K ) with a slight modification in their stopping rules (replacing |S| < K by r > ). On the contrary, CoSaMP, SP, and HTP are restricted to ( 0 C K ) since they require the sparsity level K as an input argument.

Algorithm 4: CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF] Input: y, H, K Output:

x solving ( 0 C K ) 1 Initialization: x ← 0; S ← ∅;

Non-negative greedy algorithms

Some greedy algorithms such as MP and IHT, due to the simplicity of their iterations, can be turned to non-negative setting quite simply. However, it is not the case for orthogonal greedy algorithms such as OMP and OLS. Taking non-negativity into account within OMP and OLS can raise some non-trivial issues since the structure of the algorithm is deeply modified and, more importantly, proposing efficient implementation is not obvious. Let us now explain the issues and how they are addressed in the literature.

Non-negative extensions of OMP The first non-negative extension of OMP was introduced by Bruckstein et al. under the name Non-negative OMP (NNOMP) [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF]. This extension naturally considers only positive inner products in the selection step (see Algorithm 8). In addition, the coefficient update is done by solving a Non-Negative Least Squares (NNLS) subproblem x(S) ← arg min

z≥0 y -H S z 2 (1.6)
instead of the ULS subproblem (1.2). Since the NNLS solution might differ from the ULS solution when some non-negativity constraints are activated, the orthogonality between the residual vector and the selected subspace might not hold anymore. Besides, unlike the ULS problem, NNLS does not have a closed-form solution, so one needs to run an iterative subroutine to compute the NNLS solution at each iteration. As a result, NNOMP and other non-negative versions of OMP and OLS are expected to be much more expensive than OMP and OLS [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF], Yaghoobi et al., 2015, Yaghoobi and Davies, 2015]. The same argument applies to the nonnegative extensions of CoSaMP, SP and HTP. However, the non-negative extensions of OMP and OLS can still have fully recursive implementation as explained in Chapter 2 while it seems not to be the case for the non-negative extensions of CoSaMP, SP and HTP. The fast implementation of NNOMP was not addressed in [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] but later in [Yaghoobi et al., 2015, Ramamurthy et al., 2014, Kim and Haldar, 2016, Wang et al., 2018]. Yaghoobi et al. [Yaghoobi et al., 2015] proposed a scheme named Fast NNOMP (FNNOMP) which combines the selection and the coefficient update steps in one so that the residual vector is orthogonal to the selected subspace at every iteration. Thus, FNNOMP does not require to solve any NNLS subproblem but it might return a different output from that of NNOMP. In other words, FNNOMP can be seen as a variant of NNOMP without NNLS solving. Another variant of NNOMP was proposed by [START_REF] Wang | Cone-based joint sparse modelling for hyperspectral image classification[END_REF]. This variant, named NN-OMP, selects the new atom the same way as OMP but updates the coefficients by solving a NNLS subproblem. Thus NN-OMP and NNOMP have roughly the same costs but likely different outputs.

Non-negative extensions of OLS [Yaghoobi and Davies, 2015] firstly introduced a nonnegative extension of OLS named Non-Negative OLS (NNOLS). This version (presented in Algorithm 10) replaces all the orthogonal projections (1.5) by NNLS subproblems

∈ arg min i / ∈S min z≥0 y -H S∪{i} z 2 .
(1.7)

As a result, NNOLS requires to solve as many NNLS subproblems as the number of non-selected atoms at each iteration, hence it was considered to be too expensive [Yaghoobi and Davies, 2015].

Yaghoobi and Davies proposed two variants of NNOLS namely Suboptimal NNOLS (SNNOLS) and Fast NNOLS (FNNOLS) [Yaghoobi and Davies, 2015]. SNNOLS (presented in Algorithm 9) is based on OLS selection rule (see Algorithm 3) but it only considers positive inner products.

Besides, SNNOLS requires to solve only one NNLS subproblem related to the augmented support at each iteration. FNNOLS was inspired from SNNOLS the same way as FNNOMP was inspired from NNOMP. The two approximate schemes FNNOMP and FNNOLS are not presented here because we have examples showing that FNNOMP and FNNOLS do not always respect the positivity constraint [START_REF] Nguyen | An optimized version of non-negative omp[END_REF]. NNOMP and SNNOLS (Algorithms 8-9) can be seen as intuitive extensions to the nonnegative setting of OMP and OLS (Algorithms 2-3), respectively. Besides, NNOLS is a direct non-negative extension of OLS in optimization viewpoint. Indeed, NNOLS is designed as a descent algorithm aiming to decrease at most the residual using a single change in the support. It is noteworthy that, while the two interpretations of OLS ((1.5) and line 3-Algorithm 3) are equivalent [START_REF] Blumensath | On the difference between Orthogonal Matching Pursuit and Orthogonal Least Squares[END_REF], NNOLS and SNNOLS significantly differ since (1.5) and line 3-Algorithm 3 are no more equivalent in the non-negative setting. Indeed, the SNNOLS selection rule (line 3-Algorithm 9) can be interpreted as

∈ arg min i / ∈S min u,v≥0
y -H S uvh i 2 while (1.7) can be rewritten as

∈ arg min i / ∈S min u≥0,v≥0
y -H S uvh i 2 (see Chapter 2 for more details). From this point of view, we can see that SNNOLS selection rule is less costly than that of NNOLS. Finally, the three algorithms NNOLS, SNNOLS and NNOMP can address either ( 0 C +) or ( 0 C K +) by a slight modification in the stopping rule (replacing |S| < K by r > ).

Algorithm 8: NNOMP [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] Input: y, H, K Output: x solving ( 0 C K +) [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF]. These extensions (named NNCoSaMP, NNSP and NNHTP, respectively) consider only the positive inner products in atom selection and replace the orthogonal projection by NNLS solving, just as the way NNOMP [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] was designed. NNCoSaMP is displayed in Algorithm 11 for an illustration. It is noteworthy that, while CoSaMP, SP and HTP produce an exactly K-sparse solution at each iteration, their non-negative extension do not. This is because the NNLS solution differs from the ULS solution (and is sparser) when some non-negative constraints are activated. More importantly, there is no obvious way to enforce the K-sparseness of the output of NNCoSaMP, NNSP and NNHTP. This is one drawback when using these algorithms in practice. In addition, NNCoSaMP, NNSP and NNHTP are restricted to address ( 0 C K +) as they requires the prior knowledge on the sparsity level K of the solution. Finally, these algorithms are hard to lend themselves to recursive implementation since several support elements are changed at each iteration.

Algorithm 11: Non-Negative CoSaMP (NNCoSaMP) [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF] Input: y, H, K Output: x solving ( 0 C K +) 1 Initialization: x ← 0; S ← ∅; r ← y ; 2 while not stop do

Convex solvers

Since 0 minimization is NP-hard and several local minimizers may exist, an alternative approach is to replace 0 by 1 to obtain a convex optimization problem which is easier to solve (and there exist fast solvers) and has a single minimizer. However, the 1 solution might be very far from the global minimizer of the 0 problem (except for easy problems where 0 -1 equivalence holds).

1 minimization problem which is often referred to as LASSO [Tibshirani, 1996, Friedman et al., 2007] or Basic Pursuit [Chen et al., 2001, van den Berg and[START_REF] Van Den Berg | Probing the pareto frontier for basis pursuit solutions[END_REF] can take one of the following forms

min x x 1 s.t. y -Hx 2 ≤ ( 1 C ) min x y -Hx 2 2 s.t. x 1 ≤ K ( 1 C K ) min x y -Hx 2 2 + λ x 1 ( 1 P )
Note that the three 1 minimization problems are equivalent in the sense that we can find a value of λ so that ( 1 P ) solution identifies with that of ( 1 C ) or ( 1 C K ). Moreover, it was shown that when λ tends to zero, the solution of ( 1 P ) converges to that of ( 1 C ) for = 0 [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF]. 1 minimization problems can be solved by many methods such as interior point, gradient projection, homotopy, Lagrange multiplier, etc. Popular 1 solvers include ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], FISTA [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF], LARS [START_REF] Efron | Least angle regression[END_REF] and Homotopy [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF]. Originally, LARS [START_REF] Efron | Least angle regression[END_REF] was designed as a forward algorithm. LARS with the LASSO modification [START_REF] Efron | Least angle regression[END_REF] is forward-backward and it is identical to Homotopy. While ADMM an FISTA only solve ( 0 P ) approximately for a particular value of λ, LARS and Homotopy solve it exactly under mild assumptions and produce the entire solution path for all values of λ (for undercomplete case in [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF] and overcomplete case in [START_REF] Efron | Least angle regression[END_REF]). Interestingly, both Homotopy and LARS use greedy strategies and they are shown to be structurally close to OMP [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF]. It is noteworthy that since 1 solvers tend to produce biased coefficents [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF], usually a post-processing step is done to debias amplitudes S ← supp (x), x(S) ← H † S y.

(1.8)

Addressing a well-defined convex problem, 1 solvers can be extended to the non-negative setting without many difficulties. ADMM [START_REF] Barbu | A new approach for volume reconstruction in tomoPIV with the alternating direction method of multipliers[END_REF] and Non-negative LARS (NLARS) [START_REF] Efron | Least angle regression[END_REF], Morup et al., 2008] are typical examples among others. Note that NLARS is the non-negative extension of modified LARS (i.e., non-negative extension of Homotopy). Since we are interested in greedy structure and non-negativity, NLARS is worth to be commented. NLARS returns the solution path of non-negative LASSO

min x≥0 y -Hx 2 2 + λ x 1 ( 1 P +)
for all values of λ. It is notable that

( 1 P +) is equivalent to min x≥0 y -Hx 2 2 + λ1 t x
which is a standard quadratic program for a certain value of λ hence it can also be solved by active-set method [START_REF] Nocedal | Numerical optimization[END_REF]. The main idea of NLARS is to maintain 1 optimal condition when λ is changing. Each iteration consists of either an addition or a removal of one atom from the support and the update of sparse approximation Hx along the direction equiangular to the selected atoms. The main structural difference with LARS is that NLARS considers only the positive inner products r t h i in the atom selection. In terms of complexity, NLARS requires the update of the inverse Gram matrix (H t S H S ) -1 at each iteration but it does not perform orthogonal projection H † S y. Finally, NLARS can also be turned to obtain a K-sparse solution (see Matlab code accompanying to the paper [START_REF] Morup | Approximate l0 constrained non-negative matrix and tensor factorization[END_REF]).

Non-convex solvers

Since 1 solvers tend to produce less sparse solutions and biased estimates for large coefficients [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF], another approach is to replace 0 by a non-convex separable function which is singular at the origin to promote sparsity and has vanished derivative for large values to be unbiased. Such non-convex functions include the smoothly clipped absolute deviation (SCAD) [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF], the p pseudonorm ( x p p with 0 < p < 1) [START_REF] Frank | A statistical view of some chemometrics regression tools[END_REF]] and logarithm i log(|x i | + η) [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF] among others. The SCAD function is piecewise smooth. It is defined in 1-D by clipping three 1-D functions with separable regularization. Intuitively, it looks like 1 norm around the origin and flat at large values. The p pseudonorm (with 0 < p < 1) and logarithm are quasi-smooth approximations of the 0 norm. To minimize the non-convex objective function, a usual technique is to locally approximate it by a convex (linear or quadratic) function then use one or many iterations of proximal gradient methods or Newton's method to solve the resulting convex optimization problem [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF], Huang et al., 2008, Saab et al., 2008]. As a result, non-convex solvers often suffer from a rather high computational complexity compared with 1 solvers and greedy algorithms. Another technique is to solve the non-convex minimization problem through a sequence of weighted 1 minimization problems [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF]. This is a special case of DC programming [Thoai, 1988, Le Thi et al., 2013]. In the literature, non-negative extension of non-convex solvers is quite rare.

State-of-the-art exact recovery analyses of greedy algorithms

Exact recovery with greedy algorithms

OMP is guaranteed to recover the K-sparse representation of a data signal in the noiseless case if the mutual coherence µ = max i =j

|h t i h j |
of the dictionary H satisfies [Tropp, 2004, Herzet et al., 2013] µ < 1 2K -1 .

(1.9)

This condition was proved to be sharp [START_REF] Cai | Stable recovery of sparse signals using an oracle inequality[END_REF] and extensions to noisy cases were considered [Ben-Haim et al., 2010,Cai andWang, 2011]. Besides, relaxations of (1.9) given partial information on the decomposition of the data signal were proposed [START_REF] Herzet | Exact recovery conditions for sparse representations with partial support information[END_REF], Herzet et al., 2016]. Other exact recovery analyses of OMP also consider ERC condition [Tropp, 2004]

max i / ∈S * H † S * h i 1 < 1 (1.10)
(where S * presents the true support of the K-sparse representation) or the Restricted Isometry Property (RIP) condition of order K [START_REF] Davenport | Analysis of orthogonal matching pursuit using the restricted isometry property[END_REF],Satpathi et al., 2013,Li et al., 2015[START_REF]A sharp condition for exact support recovery with orthogonal matching pursuit[END_REF]]

(1 -δ) x 2 ≤ Hx 2 ≤ (1 + δ) x 2 (1.11)
for all K-sparse vectors x (where δ presents the RIP constant). While the mutual coherence based condition (1.9) is intuitive and its computation is feasible, it is stronger than (i.e., it implies) the ERC condition (1.10) [Tropp, 2004] and might never hold for very large size dictionaries [START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF]. On the contrary, evaluating ERC condition (1.10) requires the knowledge of the true support S * which is often unavailable and evaluating RIP condition (1.11) is computationally intractable in general [Tillmann, 2015].

Exact support recovery properties of OLS have been studied more recently. In [START_REF] Herzet | Exact recovery conditions for sparse representations with partial support information[END_REF], the mutual coherence based condition (1.9) was extended to OLS in the noiseless case, while [START_REF] Soussen | Joint kstep analysis of orthogonal matching pursuit and orthogonal least squares[END_REF] considered ERC condition and [Wen et al., 2017a,Wang andLi, 2017] used RIP assumptions.

Exact recovery analyses of CoSaMP, SP, IHT and HTP are all based on RIP condition [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], Dai and Milenkovic, 2009, Blumensath and Davies, 2008, Foucart, 2011].

Exact recovery with non-negative greedy algorithms

It was stated in [Bruckstein et al., 2008, Theorem 3] that the mutual coherence based condition (1.9) is a sufficient condition for exact recovery of any K-sparse representation using NNOMP [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] in the noiseless setting. However, the proof is not provided. Besides, Kim et al. generalized (1.9) for the case of multiple measurement vector (i.e., data y is a matrix instead of a vector) in the noiseless setting but the authors acknowledged that their condition turns out to be very restrictive [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF].

In the literature, the exact recovery analysis of NNOMP in noisy cases has not been addressed yet. The theoretical guarantees of other non-negative greedy algorithms are still missing too.

Contributions

Targeting 0 minimization, we are interested in greedy algorithms since their structure is simple and fast implementation is available (for standard greedy algorithms). On the contrary, the 1 global minimizer may be far from the 0 global minimizer and non-convex solvers are with rather high computational complexity. Orthogonal greedy algorithms such as OMP and OLS often have good performance in sparse reconstruction. Moreover, they can be easily adapted to address the constraint related to the sparsity level x 0 ≤ K and/or noise level y -Hx ≤ . It is not the case for the algorithms such as CoSaMP, SP and HTP which are also difficult to lend themselves to recursive implementation. However, the extension of orthogonal greedy algorithms to the non-negative setting is challenging. The issue arises since the orthogonal projection which is equivalent to ULS solving (1.2) is replaced by NNLS solving (1.6). Since the NNLS problem does not have a closed-form solution, solving NNLS is more computationally demanding than ULS solving. In addition, when some non-negativity constraints are activated, the NNLS solution differs from the ULS solution so the orthogonality between the residual vector and the selected subspace does not hold anymore. Addressing these algorithmic issues is our first main contribution.

Besides, the theoretical guarantee of OMP and OLS has always been a fruitful topic. OMP and OLS are guaranteed to have exact support recovery under several types of conditions such as mutual coherence based condition, Tropp's condition ERC, RIP condition, etc. These results strengthen the interest of OMP and OLS and make them more popular. However, equivalent results for the non-negative extensions of OMP and OLS are rare and limited. Hence our theoretical study of non-negative orthogonal greedy algorithms is a second important contribution.

More precisely, our contributions can be summarized in four main points as follows, where the first contribution is split in two: analyzing structural properties ( §1.5.1) and proposing fast implementation ( §1.5.2) of non-negative greedy algorithms. The second contribution appears in §1.5.3. The last contribution is about NP-hardness of the non-negative 0 minimization and can be found in §1.5.4.

Structural properties of non-negative greedy algorithms

We revisit the non-negative versions of OMP and OLS in a unified way. We introduce a new class of Non-Negative Orthogonal Greedy (NNOG) algorithms which covers many existing nonnegative extensions of OMP and OLS such as NNOMP, NNOLS and SNNOLS, up to some modifications. Moreover, we show that NNOG algorithms share many desirable properties, which are not obvious for existing non-negative extensions of OMP and OLS, such as: the norm of residual is decreasing at each iteration and the residual vector is orthogonal to the selected atoms. This is possible because a so-called support compression step is included in NNOG structure to make the support set identifies to the support of the current NNLS solution at each iteration. This contribution is presented in Chapter 2.

Fast implementation

As the NNLS problem does not have a closed-form solution, solving it is computationally demanding. Existing implementations of non-negative greedy algorithms either solve NNLS subproblems independently (hence in an exact but time-consuming way) or avoid NNLS solving (hence yielding fast but approximate schemes). Unlike any of previous contributions, we propose a fast and exact implementation of NNOG algorithms. This implementation exploits the active-set algorithm [START_REF] Lawson | Solving least squares problems[END_REF] for NNLS solving, and makes use of warm start initialization. Warm start is a popular technique but it is especially efficient when it is combined with active-set NNLS in NNOG implementation. Moreover, we elaborate on recursive implementations of NNOG algorithms and design further types of accelerations. According to our comparisons, the computing time of NNOG algorithms using proposed implementation becomes comparable to that of orthogonal greedy algorithms. This contribution is presented in Chapter 2 with numerical results in Chapter 3.

Theoretical guarantees

In the literature, theoretical analyses of NNOMP are rare and somewhat discordant while the analysis of other NNOG algorithms is missing. In this context, we establish a unified exact support recovery analysis for non-negative extensions of OMP and OLS using mutual coherence based condition. This is a consequence of a more fundamental result about OMP and OLS which is named sign preservation property. More precisely, we show that under mutual coherence based condition (1.9), (i) OMP and OLS outputs identify with the true support of the K-sparse representation in K iterations, and (ii) the weights of selected atoms have the correct sign at any of the K iterations. As a result, NNOG algorithms recover exactly the true support of the K-sparse representation with non-negative weights in K iterations under the condition (1.9). It should be noted that while the very few existing analyses of NNOMP are made in the noise-free case, our analysis is made in the noisy cases and it also applies to other algorithms such as NNOLS and SNNOLS. This contribution is presented in Chapter 4.

NP-hardness of non-negative 0 minimization

The final contribution of the thesis is about the NP-hardness of non-negative 0 minimization problems. The reasoning is the same as that of Natarajan [Natarajan, 1995] and our analysis is a direct extension. However, we discovered that some existing analyses have errors and clarify these issues. Besides, our analysis can apply to more general cases. This contribution is presented in Chapter 5.

List of publications

Preprints

Thanh T. Nguyen, Charles Soussen, Jerome Idier, El-Hadi Djermoune. Sign Preservation Analysis of Orthogonal Greedy Algorithms. January 2019. hal-01971697.

Conference papers

Presentations

• Ph.D student day at L2S, Gif-sur-Yvette, France, June 28, 2018.

• Workshop on Sparsity and Applications, Bordeaux, France, May 3, 2018.

• Ph.D. student day at CRAN, Nancy, France, January 12, 2018.

• Seminar at L2S, Gif-sur-Yvette, France, December 8, 2017.

• Seminar at LS2N, Nantes, France, October 12, 2017.

Chapter 2

Unified view of non-negative orthogonal greedy algorithms In this chapter, we revisit the non-negative versions of OMP and OLS in a unified way. In particular, we introduce a new class of non-negative orthogonal greedy (NNOG) algorithms which covers many existing non-negative versions of OMP and OLS, up to some modifications. We exhibit the structural properties of NNOG algorithms and introduce a fast implementation. The proposed fast implementation of NNOG algorithm will be assessed afterward by a set of numerical results in Chapter 3.

Introduction

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]. OMP and OLS differ in the way the new atom is selected. In both cases, the atom inducing the largest decrease of the norm of the residual 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., [Miller, 2002]. 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].

In many applications, 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. Some classical sparse algorithms can be straightforwardly adapted to deal with non-negativity 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]Pesquet, 2011, Boyd et al., 2011], 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], Le Thi et al., 2013]. On the contrary, the nonnegative extension of greedy algorithms is a challenging issue since the unconstrained leastsquares subproblems are replaced by non-negative 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 [Björck, 1996] such as active-set [START_REF] Lawson | Solving least squares problems[END_REF], interior-point [Wright, 1992], and gradient-projection [START_REF] Benvenuto | Nonnegative least-squares image deblurring: improved gradient projection approaches[END_REF]. The latter two families typically require to tune some stopping criteria empirically, resulting in approximate resolution of the NNLS problem. Here, we are focusing on active-set methods for NNLS solving 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 [Slawski andHein, 2011, Slawski andHein, 2013] and Sparse NNLS or Reverse Sparse NNLS [START_REF] Peharz | Sparse nonnegative matrix factorization with l0 constraints[END_REF].

Several existing contributions deal with orthogonal greedy algorithms in the non-negative case. Non-Negative OMP (NNOMP) was first proposed by [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 [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 (FNNOMP), which avoids solving NNLS subproblems but rather recursively approximates the sought solution using QR matrix factorization [Yaghoobi et al., 2015]. Although FNNOMP is much faster than canonical NNOMP, it is an approximate version likely to deliver a different output. In [Yaghoobi and Davies, 2015], 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., nk NNLS problems have to be solved at iteration k, with n being the size of the dictionary. Yaghoobi et al. [Yaghoobi and Davies, 2015] further proposed two accelerated versions of NNOLS named Suboptimal NNOLS (SNNOLS) and Fast NNOLS (FNNOLS). SNNOLS [Yaghoobi and Davies, 2015] selects the atom that is positively correlated with the current residual whose projection forms a maximum angle with the residual. Besides, 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 FN-NOLS 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 non-negativity 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], Yaghoobi et al., 2015,Yaghoobi and Davies, 2015], 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.

In this chapter, 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 non-negative orthogonal greedy algorithms exploiting the active-set algorithm [START_REF] Lawson | Solving least squares problems[END_REF] for NNLS solving, and based on warm start initialization. Moreover, we elaborate on recursive implementations and we design further types of accelerations.

The chapter is organized as follows. Section 2.2 introduces the family of so-called Nonnegative 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 2.3, we propose a fast implementation based on recursivity and on the use of warm starts for solving the NNLS subproblems. Section 2.4 is devoted to NNOG acceleration. Finally, discussion and perspectives will be found in Section 2.5.

Non-negative greedy algorithms

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 +)
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 (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 hS i = 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 leastsquares (NNLS) solution corresponding to S, any vector x in R n or R n + , respectively, that minimizes y -Hx 2 with the constraint that 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 chapter, 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 [Lawson and Hanson, 1974, Chap. 3], which are recalled next for completeness.

Lemma 2.1. Consider the NNLS problem related to support S:

min x≥0 y -Hx 2 s.t. supp(x) ⊂ S.
(2.2)

x+ S is a solution to (2.
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 (2.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.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.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.
(2.4)

For quadratic programming problems involving positive semidefinite matrices, the KKT conditions are necessary and sufficient conditions of optimality [Nocedal and Wright, 2006, Chap. 16], so x+ S is a solution to (2.2) if and only if (2.5)

     λ = -2H t (y -H x+ S ) λ(C) = 0 λ(S\C) ≥ 0 that is, when (2.3) is satisfied.
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.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 (2.3).

Positive supports will play an important role in our specification of fast non-negative orthogonal greedy algorithms.

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 12. 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.

Algorithm 12: General structure of a non-negative orthogonal greedy algorithm to solve ( 0 +). input : y, H, K output:

x 1 x ← 0 ; S ← ∅ ; r S ← y ; 2 while |S| < K and max i∈ S h t i r S > 0 do 3 Select an index ∈ S by a selection rule S(y, H, S) ; 4 S ← S ∪ { } ; 5 x ← x+ S ; 6 S ← supp(x) ; 7 r S = y -Hx ; 8 end
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 at each iteration 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 ). According to the above definition of NNOG algorithms, distinct algorithms can only differ by the selection rule used to select an atom at each iteration. The design of a selection rule corresponds to the definition of a function S(y, H, S), taking values in S. It is clear that some indices ∈ S correspond to inappropriate choices, in the sense that their selection would produce x+ S∪{ } ( ) = 0, and hence a useless iteration, and possibly an early stopping of the algorithm. In contrast, in the unconstrained case, any selection ∈ S yields a decrease of y -Hx 2 unless h ∈ span(H S ). The capacity of some selection rules to avoid inappropriate selections is examined in the next two subsections.

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 for NNOG algorithms compared to their unconstrained counterparts. However, NNOG algorithms lend themselves to recursive implementations akin to usual orthogonal greedy schemes, as detailed in Section 2.3.

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.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 2.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 2.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 2.1. The descending atom condition i ∈ D S is equivalent to

0 < h t i r + S .
(2.7)

When S is a positive support, it is also equivalent to each condition

0 < h t i r S , (2.8) 0 < gt i r + S ,
(2.9)

0 < gt i r S .
(2.10)

When H S∪{i} is full column rank, it is also equivalent to x+ S∪{i} (i) > 0.
(2.11)

Proof. See Appendix 2.A.1.
Let us remark that from the first item of Proposition 2.1, the selection rule is invoked at Line 3 of Algorithm 12 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 three lemmas have interesting consequences for the practical specification of valid NNOG algorithms.

Lemma 2.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 a positive support, 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 corresponding to S ∪ {i}, and also an NNLS solution corresponding to S ∪ {i} since xS ≥ 0. This implies that i ∈ D S .

Lemma 2.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 2.3 and 2.5, starting with the empty support.

Let us stress that Lemma 2.6 refers to NNOG algorithms strictly conforming to the scheme of Algorithm 12 (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 2.5 applicable. To our best knowledge, such a compression step has not been proposed in any previous contribution about nonnegative greedy schemes targeting 0 minimization2 .

According to Lemma 2.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 12. This simplification is adopted in the rest of the chapter. 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.

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] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF], Ramamurthy et al., 2014, Yaghoobi et al., 2015, Kim and Haldar, 2016]:

S 1 (S) ∈ arg max i / ∈S h t i r S
(2.12)

• Suboptimal NNOLS (SNNOLS, [Yaghoobi and Davies, 2015]) rule:

S 2 (S) ∈ arg max i / ∈S gt i r S (2.13)
• NNOLS rule [Yaghoobi and Davies, 2015]:

S 3 (S) ∈ arg min i / ∈S r + S∪{i} 2 (2.14)
(2.14) is a descent selection rule by definition. The fact that (2.12) and (2.13) are descent selection rules is deduced from recursive application of Proposition 2.1 and Lemma 2.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 (2.13), so we can expect the computing cost of (2.13) to be larger than that of (2.12). Rule (2.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 [Yaghoobi et al., 2015, Yaghoobi andDavies, 2015] rely on a descent selection rule. They will therefore not be further analyzed.

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

Proposition 2.2. Rules (2.12)-(2.14) are equivalent to

S j (S) ∈ arg min i∈D S µ j (S, i),
(2.15)

where µ j are specific to each rule:

µ 1 (S, i) = min v y -H S x S -h i v 2 , (2.16) µ 2 (S, i) = min u,v y -H S u -h i v 2 ,
(2.17)

µ 3 (S, i) = min u≥0,v≥0 y -H S u -h i v 2 .
(2.18)

Proof. Let us first emphasize that because (2.12)-(2.14) are descent selection rules, i / ∈ S in (2.12)-(2.14) can be replaced by i ∈ D S . Clearly, (2.15) with (2.18) simply duplicate (2.14). Using identity (2.1) with r = r S and h = h i , (2.12) can be rewritten as the argmin over i ∈ D S of min

v≥0 y -H S x S -h i v 2 . (2.19)
Likewise, applying (2.1) with r = r S and h = gi and using the fact that gi is the normalized version of hi , (2.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 .
(2.20)

(2.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 )) Figure 2.1: Graphical illustration of (2.16)-(2.18) in a simple case where S contains a single atom, whose weight u varies along the horizontal axis, while the weight v of the candidate atom varies along the vertical axis. The contour plot represents yh S uh i v 2 .

On the usefulness of support compression

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

• 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 (2.12) and (2.14), with r S replaced by the current residual r + S ), Proposition 2.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 (2.13) is a descent selection rule when S is not positive, since the equivalence (2.9) in Proposition 2.1 does not hold anymore.

• Finally, for any NNOG algorithm, Lemma 2.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 2.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.

Active-set NNLS algorithms and recursivity

Let us consider the NNLS problem related to support S in (2.2). 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, active-set NNLS (AS-NNLS) algorithms are well-known schemes that solve NNLS problems in a finite number of iterations [Lawson andHanson, 1974, Chen et al., 2010]. 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 | [END_REF]Pernkopf, 2012, Foucart andKoslicki, 2014].

Fast active-set algorithms

Among many numerical methods for solving (2.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 structure is an essential element to obtain practically fast implementations [START_REF] Lawson | Solving least squares problems[END_REF], Björck, 1996, Nocedal and Wright, 2006]. 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 [Zhang, 2011]), also called stepwise algorithms [Miller, 2002, 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.

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 [Nocedal and Wright, 2006, Chap. 16], Nocedal and Wright proposed an AS-NNLS algorithm with any feasible vector as a possible starting point. Algorithm 13 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 13 identifies with Nocedal and Wright's scheme.

Algorithm 13: Active-set algorithm to solve the NNLS problem related to T , starting from a positive support S. Format: AS_NNLS(y, H, T, S, xS ) input : y, H, target set T , initial support S ⊂ T , xS output: At any iteration of a valid NNOG algorithm, the current solution x+ S can be used as the initial point to compute x+ S∪{ } using Algorithm 13 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 14 is a global view of the active-set implementation of NNOG obtained by integrating the calls to the AS-NNLS solver (Algorithm 13) in the NNOG framework (Algorithm 12). 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 13, Line 7) and then alternates between atom selections and series of deselections. This mechanism is illustrated by a simple example in the next subsection.

V := supp(x + T ), x+ T := xV 1 x ← xS ; V ← S; r V ← y -Hx ; 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 15); 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 \ { -} ;
A reduced number of iterations is obtained because we use Algorithm 13 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 2.B, Algorithm 15 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 13 (Lines 4-5 and 10-11).

2.3.3

Step-by-step illustrative example of NNOG Fig. 2.2 displays 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 ← S k-1 ∪ { }. Therefore, S 1 ⊂ S 2 ⊂ S 3
are of cardinalities 1, 2 and 3, respectively. At iteration 4, S 3 ∪ { } is not positive so AS-NNLS performs two support changes, namely the selection of and a deselection. The next NNOG iterate reads S 4 ← S 3 ∪ { }\{ 1 }. Iteration 5 is more tricky (and unlikely). 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 13, with T ← S 4 ∪ { }) is not met since h 1 is a descending atom. Therefore, 1 is reselected: 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.

Connection between AS-NNLS and NNOMP

In Algorithm 14, 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 l0 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 13). 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 l0 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 2.3. Let us consider Sparse NNLS, i.e., Algorithm 13 initialized with the empty support together with the early stopping rule |S| = K. In any case where no iteration produces two or more successive removals in Lines 6-12, the output of Sparse NNLS identifies with that of NNOMP (i.e., Algorithm 12 with rule S 1 in (2.12)).

Proof. See Appendix 2.A.2.

Acceleration of NNOG algorithms

The NNOMP, SNNOLS and NNOLS selection rules (2.12)-(2.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 2.1) to carry out the optimization tasks (2.13) and (2.14) over i ∈ D S only. This operation is referred to as type-I pruning of the dictionary3 . 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.

Atom selection

The atom selection step of both NNOMP and SNNOLS can be efficiently implemented using vectorized computations. 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 (2.16), so they share the same possibilities of parallel computations. Likewise, OLS and SNNOLS being both based on the minimization of (2.17), they share the same possibilities of parallel implementation. In coherence with our choice of recursive implementation of ULS solutions (see Appendix 2.B), we have adopted a MIL based solution to solve (2.17) in a vectorized way (see Matlab code in supplementary material).

In contrast, the NNOLS selection rule (2.14) 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 each 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 13 yields the same unconstrained minimizer xS∪{i} as the one involved in the SNNOLS rule (2.17). 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.

Coefficient update

Once an atom i is selected, NNOG algorithms need to update the coefficients by solving an NNLS subproblem (Algorithm 12, 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 the AS-NNLS algorithm (Algorithm 13) from the initial set S to the target set S ∪ { }. For SNNOLS, a call to Algorithm 13 is needed only when S ∪ { } is not positive.

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.

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 [Yaghoobi et al., 2015, Yaghoobi andDavies, 2015] 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 l0 constraints[END_REF] for an interesting case of sparse NMF).

In this chapter, our first contribution 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 non-negative 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.

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 CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], 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.

Appendix 2.A Proof of technical results

2.A.1 Proof of Proposition 2.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 (2.1) for r = r + S , we deduce that i ∈ D S since

r + S∪{i} 2 ≤ min v≥0 r + S -vh i 2 < r + S 2 .
Conversely

f (z S∪{i} ) < f (z S ). Since f is convex, one has (z S∪{i} -z S ) t ∇f (z S ) ≤ f (z S∪{i} ) -f (z S ) < 0 (2.21)
where the gradient of f is defined by

∇f (z S ) = 2H t S∪{i} (H S∪{i} z S -y) = -2H t S∪{i} r + S .
Denoting by C := supp(x

(h t i r + S ) x+ S∪{i} (i) > 0 (2.24)
and thus h t i r + S > 0. Let us now assume that S is a positive support. According to Lemma 2.2, we have x+ S = xS and r + S = r S , so (2.7) and (2.8) are identical, as well as (2.9) and (2.10). To show that (2.7)-(2.8) are equivalent to (2.9)-(2.10), 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, (2.8) rereads 0 < ht i r S , which implies that hi = 0 and gi = hi / hi = 0. Hence, (2.8) rereads 0 < gt i r S , which identifies with (2.10). Finally, let us show that i ∈ D S is equivalent to condition (2.11). 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+ 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 .

2.A.2 Proof of Proposition 2.3

Any iteration of AS-NNLS starts with the addition of a new atom to the current support (Line 4 of Algorithm 13). 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 first 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 2.3.2. According to Lemma 2.1, the KKT conditions read:

H t V (y -H xV ) = 0, (2.25a) h t -(y -H xV ) ≤ 0.
(2.25b)

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

Appendix 2.B Recursive implementation of ULS

Algorithm 15 recalls the recursive ULS computation using MIL [Björck, 1996] 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 15. Notation -j refers to all indices except j, and θ j stands for the j-th column of Θ.

For reminding, MIL for block matrix can be read as [Golub and Van Loan, 1996]

A B C D -1 = A -1 + A -1 B(D -CA -1 B) -1 CA -1 -A -1 B(D -CA -1 B) -1 -(D -CA -1 B) -1 CA -1 (D -CA -1 B) -1 (2.26)
When a new atom is added to the support, S = S ∪ { } and

H S = [H S , h ]. Therefore, H t S H S = H t S H S H t S h h t H S h i 2 .
By applying (2.26), one can express the inversion (H t S H S ) -1 by means of (H t S H S ) -1 as formulated at line 7 of Algorithm 15. When an atom is removed from the support, S = S \ { }. Let us firstly consider the case where h is the last column of H S . Then H S = [H S , h ]. By switching S and S in the calculation of the case of atom addition, one can express the inversion (H t S H S ) -1 by means of (H t S H S ) -1 as formulated at line 13 of Algorithm 15. The case where h is not the last column of H S can be turned to the previous case by using a permutation matrix

P j = [e 1 , • • • , e j-1 , e |S| , e j , • • • , e |S|-1 ]
2.B. Recursive implementation of ULS Algorithm 15: Recursive ULS [Björck, 1996]. Format: ULS(y, H, V, fw, , x, Θ, e 2 )

1 if fw then 2 φ ← H t V h ; 3 δ ← (1 -φ t Θφ) -1 ; 4 β ← φ t x(V ) -h t y ; 5 e 2 ← e 2 -δβ 2 ; 6 x(V ∪ { }) ← x(V ∪ { }) + δβ Θφ -1 ; 7 Θ ← Θ 0 0 0 + δ Θφ -1 Θφ -1 t ; 8 V ← V ∪ { } ; 9 else 10 j ← index of in V ; 11 e 2 ← e 2 + (x( )) 2 /θ j (j) ; 12 x(V ) ← x(V ) -x( ) θ j /θ j (j) ; 13 Θ ← Θ(-j, -j) -θ j (-j)θ j (-j) t /θ j (j) ; 14 V ← V \ { } ; 15 end
where j is the index of in S and e k is the vector of size |S| × 1 whose k-th entry is equal to and other entries are equal to 0. Note that P j is an orthogonal matrix and its inversion is also a permutation matrix

P -1 j = P t j = [e 1 , • • • , e j-1 , e j+1 , • • • , e |S| , e j ].
One can easily check that H S P j moves the last column of H S to the j-th column (and keeps the same order for the other columns) and H S P t j moves the j-th column of H S to the last column. The calls to Algorithm 15 for updating ULS solutions at Lines 4-5 and 10-11 of Algorithm take the respective forms:

ULS(y, H, V, 1, + , xV , Θ, r V 2 ), ULS(y, H, V, 0, -, xV , Θ, r V 2 ).
Chapter 3

Numerical validation of NNOG algorithms

Introduction

In the previous chapter, we presented our unified framework of NNOG algorithms and proposed a (recursive) fast implementation. The purpose of this chapter is to assess the proposed fast implementation of NNOG algorithms by a set of numerical results. Firstly, we introduce the problems considered in our tests which include: (i) sparse deconvolution with Gaussian kernel, (ii) sparse deconvolution for super resolution and (iii) decomposition of real-world near-infrared (NIR) spectra into elementary Gaussian features. Secondly, we validate the accelerations of NNOG algorithms proposed in Section 2.4. Thirdly, we compare the performance of several NNOG algorithms. Finally, we compare NNOG algorithms with competing algorithms which include: (i) unconstrained greedy algorithms, (ii) fast approximate non-negative greedy algorithms, (iii) non-negative extensions of CoSaMP, SP, HTP and (iv) non-negative extension of LARS. The results presented in this chapter are achieved on a macOS X system with 16 GB RAM and Intel Core i7 processor at 2.7 GHz. Before going further, let us emphasize that in applications (microscopy, tomography, etc...), the data are discrete. Therefore, two grids will be considered: (i) the sampling grid corresponding to the data domain (e.g., the pixels of the radiograph in Tomo-PIV), and (ii) the grid corresponding to the sparse solution (e.g., the voxels of the discretized 3-D volume in Tomo-PIV). Note that in sparse deconvolution both data and solution domain coincide (e.g., the time domain). However, the practitioner may set the solution grid according to the desired resolution level.

Problems

Sparse deconvolution with Gaussian kernel

We consider a convolution problem with a normalized Gaussian kernel h of standard deviation σ. h is approximated by a finite impulse response of length 6σ by thresholding the smallest values. The discrete convolution h * x * can be rewritten as a matrix-vector product Hx * where H is a Toeplitz matrix containing all possible delayed version of the Gaussian kernel h

= [h 1 , • • • , h p ], i.e., H =             h 1 0 . . . 0 . . . h 1 . . . . . . h p . . . . . . 0 0 h p . . . h 1 . . . . . . . . . . . . 0 . . . 0 h p             . (3.1)
H is a slightly under-complete dictionary (of size m × n with n = mp + 1) since we set the boundary condition so that the supports of all atoms are included in the observation window. Simulated data are generated according to

y = Hx * + (3.2)
where x * and stand for the ground truth coefficients 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 set to either a positive constant or randomly generated from a Gamma distribution. The signal-to-noise ratio is defined by SNR = 10 log 10 (P Hx * /P ) (3.3)

where

P Hx * = Hx * 2 /m (3.4)
is the average power of the noise-free data and P is the noise variance.

A simple illustration of the data decomposition is given in Fig. 3.1. The data y is sampled with step equal to 1 and it is decomposed into three Gaussians h 1 , h 2 , h 3 whose centers are located on the sampling grid. The Gaussians h 1 , h 2 , h 3 are of the same width, thus each Gaussian is a delayed version of the first Gaussian. Fig. 3.2 shows an example of simulated data y corresponding to an SNR = 30 dB. The number of samples is set to m = 300. The Gaussian kernel is of width σ = 5. The ground truth x * contains 10 spikes located at integer values. It should be noted that in this simulation the grid corresponding to the delays (locations of spikes) is the same as the one of the data signal y with step size equal to 1. In the generation of ground truth x * , the delays of Gaussians are set to integer values. Therefore, it is possible to get exact recovery results (y = Hx) in the noise-free case ( = 0) if the delays are correctly found by the algorithm.

All comparisons using this problem 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 * ;

(3.5)

• Residual norm: Euclidean norm of the data residual ( y -Hx ).

Note that here the accuracy criterion (3.5) makes sense because the true spikes are located at integer values. This criterion will be no longer applicable in the simulation of §3.2.2 since the true spikes in that simulation are located at real values.

Sparse deconvolution for super resolution

While it is more convenient to impose that the true spike are located at integer values or on a predefined grid, this situation is not feasible in real-world applications. For instance, in Tomo-PIV application (see Chapter 1- §1.2.1), the spike locations correspond to the 3-D positions of the tracer particles in the volume. Thus, when the volume is discretized by a predefined grid, it is unrealistic to expect the spike locations to be located at the center of voxels. A recent trend is to work with a continuous dictionary (off-the-grid approach) (see, e.g., [Ait Tilat et al., 2019,Elvira et al., 2019]). In the continuous approach, the dictionary contains an infinite number of atoms.

Here we do not exploit the continuous approach. Instead, we elaborate on the case of a discrete dictionary corresponding to a fine grid.

In this simulation, we consider a more realistic setting where the spikes are located at real values. The continuous signal y(τ ) is obtained by a convolution of a Gaussian function

h(τ ) = exp -τ 2 2σ 2 (3.6)
and a sparse spike train

x(τ ) = K j=1 x * j δ(τ -c * j ) (3.7)
where x * j are positive weights, δ(τ ) is the Dirac delta function and c * j ∈ R, j = 1, . . . , K stands for the location of the K spikes. Hence, the noisy continuous signal can be expressed as

y(τ ) = (h * x) (τ ) + noise (3.8) = K j=1
x * j h(τc * j ) + noise.

(3.9)

We assume, without loss of generality, that the sampling step corresponding to the acquisition of the data y is equal to 1, so that the samples y i correspond to the time τ = i. Therefore, for i = 1, • • • , m with m being the number of samples,

y i = K j=1 x * j h(i -c * j ) + noise. (3.10)
Thus, the simulated data are generated according to

y = K j=1 x * j h c * j + (3.11)
where h c * j stands for the discrete Gaussian of width σ and centered at c * j , and stands for the white Gaussian noise. The true centers c * j and true weights x * j are randomly generated with a uniform distribution. Note that we restrict the true centers to

3σ ≤ c * j ≤ m -3σ (3.12)
so that the supports of all Gaussians are included in the observation window.

For sparse reconstruction, it is usual to choose a resolution grid of step size ∆ = 1/α where α is a positive integer. We approximate the data y in (3.11) by a linear combination of normalized discrete Gaussians h c j whose centers c j are located on the resolution grid. It should be noted that in the previous simulation ( §3.2.1), as the true centers c * j are located at integer values, one can expect the exact recovery (i.e., c j = c * j ) in the noise-free case (i.e., = 0). On the contrary, in the current simulation, as the true center c * j are located at real values, one can only expect that the estimated centers will be in the same cell as the true centers (i.e., |c jc * j | < ∆/2). The related dictionary H contains all normalized discrete Gaussians h c i whose centers c i are located on the resolution grid of step size ∆ = 1/α. Thus, H has roughly αm atoms. However, since we set the boundary condition so that the supports of all atoms are included in the observation window, the number of atoms n is restricted to

n = α m -6σ -1 (3.13)
where • denotes the integer part. Note that H can be read as a block Toeplitz matrix up to reordering of its atoms. A simple illustration is shown in Fig. 3.3 where m = 9, σ = 1, n = 5, and step size ∆ = 1/2. H contains 5 normalized atoms:

H = [h 1 , h 2 , h 3 , h 4 , h 5 ] (3.14)
with support of size p = 7. Each atom is a delayed version of the previous one by a step ∆. By reordering the atoms, H can be read as a block Toeplitz matrix H = [H (1) , H (2) ] where

H (1) = [h 1 , h 3 , h 5 ] and H (2) = [h 2 , h 4 ].
Since the dictionary H can be read as a block Toeplitz matrix, it can lead to fast computations [Carcreff, 2014]. It should be noted that in this simulation all atoms are of the same width. On the contrary, in the decomposition of NIR spectra (see §3.2.3), the Gaussian dictionary is built as a block Toeplitz matrix but each block is corresponding to the convolution matrix related to the Gaussian kernel of a different width.

Decomposition of NIR spectra into elementary Gaussian features

In the two previous problems, we have moved the simulated data from a basic framework ( §3.2.1) to a more realistic one ( §3.2.2). Here we deal with real-world near-infrared (NIR) data processing. The data were obtained during Trispirabois project (2013-2016) gathering two public laboratories (CRAN and LCPME), three industrial companies (Egger, Pellenc ST and CrittBois) and funded by the french "FUI EcoIndustries 2012" program. The main objective of the project was to increase the rate of recycled wood in the production of new particleboards. In this respect, more than 300 wood samples with different composition (raw wood, plywood, particle boards, MDF-HDF) and different finition (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 near-infrared (NIR) range 3600-10000 cm -1 (which correspond 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], Wagner et al., 2015]. 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 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.

The principle of wood classification by the decomposition of NIR spectra is summarized in Fig. 3.4. The different wood pieces are placed on a black background (conveyor) and scanned using a spectrometer. The resulted data can be represented as a cube with two spatial dimensions and one spectral dimension. One pixel of the data cube is a NIR spectrum composed of different wavelengths. One decompose the NIR spectrum (into elementary Gaussian features) to identify a subset of informative wavelengths in order to detect recyclable wood samples.

For convenience, let us start by introducing the sparse decomposition of one spectrum into elementary Gaussian features. Fig. 3.5 presents a simple illustration of such a decomposition. The simulated data (in black) is decomposed into three Gaussians of different widths σ (σ = 10 for the red Gaussian, 20 for the green and 30 for the blue). The noisy data is generated such that SNR = 30 dB. Fig. 3.6 presents another illustration of the decomposition with a focus on how the dictionary looks like. The dictionary is a block Toeplitz matrix

H = [H 1 , • • • , H j , • • • , H q ] (3.15)
where each block H j is the Toeplitz dictionary corresponding to the convolution matrix related to the Gaussian impulse response of width σ j . In Fig. 3.6, H contains 3 blocks: red, green, blue (hence q = 3) where the red block corresponds to the smallest width σ and the blue block corresponds to the largest width. Hereafter, we consider the decomposition of 50 NIR spectra, seen as data vectors y of length 1647. The data pre-processing includes baseline removal, offset correction ensuring zero lower bound, and unit energy normalization. One typical spectrum is shown in Fig. 3.7. To decompose the spectra, we build a dictionary H with Gaussian-shaped columns obtained by discretizing the parameters of a Gaussian function (centers and widths). This dictionary is formed by appending the columns of the convolution dictionaries (corresponding to a fixed width σ) used in Subsection 3.2.1 for 60 equally spaced values of σ ∈ [10, 600] cm -1 . The generated dictionary is composed of 2998 atoms. Note that 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 .

Validation of NNOG accelerations

In §2.4 we proposed accelerations for NNOG algorithms. Let us remind that NNOG algorithms can be made fully recursive by using active-set NNLS initialized with warm start. Moreover, several pruning strategies (type-I and type-II) can be combined to reduce the number of tests at the atom selection stage. Finally, the latter step can benefit from vectorized computations.

In this section we compare the accelerated implementation of NNOG algorithms with their standard implementation in terms of computing time in a sparse deconvolution problem with Gaussian kernel (see §3.2.1 for description). Our comparison aims to illustrate the effect of the proposed accelerations in reducing the computing time of NNOG algorithms.

Sparse deconvolution with Gaussian kernel

We generate a convolution dictionary of size 1200 × 1140 (with σ = 10). For each value of 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: ratio, in terms of CPU time, between canonical (non-recursive) and accelerated recursive implementations of NNOG algorithms.

• Non-negativity loss: ratio, in terms of CPU time, 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 if 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 Table 3.1. Accelerated implementations yield a gain in time by a factor greater than 2. For NNOMP, the gain increases with K. Since NNOLS needs to solve many NNLS subproblems per iteration, the gain is much larger. The time gain is intermediate for SNNOLS. 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.5. 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. Note that the number of times the support compression occurs does not identify to the gap between the number of iterations and K, since support compression may remove several atoms per iteration.

Computation burden of SNNOLS and NNOLS

Hereafter, the computation burden of SNNOLS and NNOLS is assessed more thoroughly so as to evaluate the accelerations proposed in Section 2.4. Four indicators are computed at each iteration:

1. ρ ↓ = |D 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.

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

The underlying idea is that 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 candidate that are truly considered in the selection rule (2.15) 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 nonnegative ULS solution. ρ + sel is defined similar to ρ + cand . However, ρ + sel does not apply to the candidate atoms, but to the selected atoms, in order to bring information on the computational burden of the coefficient update stage of SNNOLS.

Using the dictionary of size 1200 × 1140 and the settings K = 80 and SNR = 30 dB, the computational burden of SNNOLS and NNOLS is assessed in Fig. 3.8. 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. 3.8(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 lines and columns of Table 3.1). 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. 3.8(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, . . . , 79}. 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. 3.8(a). One can observe from the curve (1ρ ↓ ) that descending atoms are numerous at early iterations, and then their rate gradually decreases for both SNNOLS and NNOLS. Therefore, type-I pruning is more effective at late iterations, where it discards about half of the atoms.

The effect of type-II pruning in NNOLS is measured by ρ + cand and ρ II indicators in Fig. 3.8(a). Large values of ρ + cand and ρ II are obtained, which indicates that the computing cost of NNOLS is dramatically reduced. Recall that a large value of ρ II implies that type-II pruning eliminates most atoms in D S that do not result in a positive support. The decrease of ρ + cand implies that the rate of descending atoms not resulting in a positive support gets larger as the number of iteration increases. However, type-II pruning remains truly effective. 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, that is, 22% of descending atoms. As a result, only 3% of descending atoms require to proceed 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.

Besides, the decrease of the indicators ρ + sel , ρ + cand and ρ II highlights that both SNNOLS and NNOLS call NNLS more often at the late iterations. A reason for this behavior is that the correlation between atoms becomes larger when the set of considered atoms enlarges.

Conclusion

The accelerated implementation of NNOG algorithms are significantly faster than the standard implementation and becomes comparable with that of unconstrained greedy algorithms. Therefore it is reasonable to use NNOG algorithms in any non-negative sparse problem instead of the traditional OMP and OLS.

Comparison of three NNOG algorithms

In this section we compare the computing time and performance of the three NNOG algorithms (NNOMP, NNOLS and SNNOLS, see Chapter 2 for their descriptions) in two simulated problems: (i) sparse deconvolution with Gaussian kernel, and (ii) sparse deconvolution for super resolution. The readers are referred to §3.2.1- §3.2.2 for the description of these problems.

Sparse deconvolution with Gaussian kernel

The curves shown in Fig. 3.9 are obtained for a convolution dictionary of size 2400 × 2388, with σ = 2 and for various settings of K ∈ [2, 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 constant; it is set to P = 10 -2 . In this simulation condition, it is noticeable in Fig. 3.9(b) that the coefficient inaccuracy xx * 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 see differences in terms of estimation quality.

In this test, NNOMP is faster but less accurate than NNOLS and SNNOLS. SNNOLS has roughly the same performance as NNOLS but with lower cost which is close to that of NNOMP. Illustration on one trial is shown in Fig. 3.10.

Super resolution problem

The noiseless data of size m = 300 are obtained by a linear combination of K = 40 discrete Gaussian features of width σ = 2 and centered at continuous locations c * j . The true centers c * j and weights x * j are randomly generated with uniform distribution. Then the data signal y is obtained by adding a white Gaussian noise corresponding to a SNR=30 dB. The resolution grid is chosen with step size ∆ = 2 -N where N varies from 1 to 7. For each value of N (hence each value of ∆), the corresponding dictionary H is built from all possible delayed versions of the Gaussian of width σ. The centers of atoms are located on the grid of chosen step ∆.

NNOMP, NNOLS and SNNOLS are run until a support of cardinality K is reached. We compare the three algorithms in terms of CPU time, the number of iterations and the reconstruction accuracy which is characterized by three factors:

• Support recovery: ratio of true positives to K

• Relative residual norm: y -Hx / y

• Normalized support distance:

  i min j |c * i -c j | + j min i |c * i -c j |   /(2K) (3.16)
where c j are estimated locations of the true centers c * j .

Note that (3.16) is a symmetric measure of the minimum distance (in average) between the location of one true spike and the corresponding estimated location. Its value is always non-zero since the true spike locations c * i are set to real values, hence, they are almost surely located off the resolution grid containing the estimated locations c i . One example of simulated data and the average performance of the three algorithms over 100 trials are shown in Figure 3.11. NNOMP yields the weakest performance while NNOLS and SNNOLS perform competitively. We observed that, for larger value of N , NNOLS yields a lower support error in comparison with SNNOLS. This means that NNOLS has a better capacity in locating true spikes even when the dictionary is highly correlated (i.e.,, when the grid step size is small, the consecutive atoms are very similar hence it is more difficult to distinguish). Besides, when N is larger than 5, SNNOLS becomes more costly than NNOLS as it requires many more iterations.

Conclusion

Among the three NNOG algorithms, NNOMP is the fastest but yields the weakest performance. NNOLS has the best capacity in locating the true spikes, specially when dealing with highly correlated dictionaries. SNNOLS has the best trade of between computing time and reconstruction accuracy.

Comparison with unconstrained greedy algorithms

In this section we perform comparisons between NNOG algorithms and (unconstrained) OMP and OLS to demonstrate the usefulness of working with NNOG algorithms. The comparisons are done in a simulated sparse deconvolution problem with Gaussian kernel (described in §3.2.1) and in the decomposition of real-world NIR spectra into elementary Gaussian features (described in §3.2.3). 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 nonnegativity constraints are activated. The resulting schemes are denoted by OMP+ and OLS+.

Sparse deconvolution with Gaussian kernel

We use the same setting as in §3.4.1. Figure 3.12 includes a comparison of NNOMP and NNOLS with OMP, OLS, OMP+ and OLS+. The performance of OMP, OLS, OMP+ and OLS+ in terms of support recovery and coefficient accuracy is weaker than those of NNOG algorithms. In Figure 3.12(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 Figure 3.12(cd), 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.

Decomposition of NIR spectra

For each NIR spectrum, NNOMP, NNOLS, OMP, OLS are run until a support of cardinality K = 20 is reached. The simple OMP+ and OLS+ algorithms are considered as well. For each competing algorithm, the CPU time and the relative residual norm y -Hx / y , averaged over 50 spectra, are displayed in Tab. 3.2. First, the time and approximation errors of OMP (resp., OLS) are both lower than those of NNOMP (resp., 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- These results are further illustrated for a specific spectrum in Fig. 3.13, where approximations and sparse recoveries are displayed for all competing algorithms. OMP and OLS yield 7 and 9 negative peaks, respectively. Besides, the OMP+ and OLS+ approximations around 4500-7000 cm -1 and 5000 cm -1 , respectively, are rather poor. NNOLS outperforms NNOMP, in particular, around 4200 and 7000 cm -1 .

Conclusion

OMP and OLS have weak capacities of reconstructing non-negative ground truth. The simple post-processing by removing negative coefficients is not sufficient to convert OMP and OLS to non-negative sparse solvers. NNOMP and NNOLS significant improve the performance of OMP and OLS in the reconstruction of non-negative ground truth. Moreover, the computing time of NNOMP and NNOLS becomes comparable to that of OMP and OLS. Therefore, it is reasonable to use NNOMP and NNOLS in any non-negative sparse problem instead of OMP and OLS.

Comparison with fast approximate non-negative greedy algorithms

In this section we compare NNOG algorithms with FNNOMP [Yaghoobi et al., 2015] and FN-NOLS [Yaghoobi and Davies, 2015]. Let us remind that FNNOMP was introduced in [Yaghoobi et al., 2015] as a fast (but approximative) implementation of NNOMP. Similarly, FNNOLS was introduced in [Yaghoobi and Davies, 2015] as a fast implementation of NNOLS. Our comparisons are done in a simulated sparse deconvolution problem with Gaussian kernel (described in §3.2.1) and in the decomposition of real-world NIR spectra into elementary Gaussian features (described in §3.2.3).

Sparse deconvolution with Gaussian kernel

We use the same setting as in §3.4.1. Figure 3.14 presents a comparison with the FNNOMP and FNNOLS implementations1 of Yaghoobi et al. [Yaghoobi et al., 2015, Yaghoobi andDavies, 2015]. 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. In contrast, our fast implementation of NNOMP is exact and the non-negativity constraint is always satisfied.

Decomposition of NIR spectra

For each NIR spectrum, NNOG algorithms, FNNOMP and FNNOLS are run until a support of cardinality K = 20 is reached. The average computing time and normalized approximation error over 50 spectra are shown in Tab. 3.3. As can be seen, NNOMP and FNNOMP have roughly same cost and performance. FNNOLS is more costly than NNOLS and SNNOLS with a slightly weaker performance. Illustration on a spectrum is shown in Fig. 3.15.

Conclusion

FNNOMP has roughly the same cost and performance as our exact implementation of NNOMP. However, FNNOMP might return some negative coefficients. Besides, FNNOLS is much slower than NNOLS and SNNOLS with slightly weaker performance. Therefore, we conclude that it is not necessary to work on approximative schemes to derive efficient implementations.

Comparison with non-negative extensions of CoSaMP, SP and HTP

In this section we compare NNOG algorithms with non-negative extensions of CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], SP [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] and HTP [Foucart, 2011] which are named NNCoSaMP, NNSP and NNHTP [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF], respectively. Our comparisons are done in a simulated sparse deconvolution problem with Gaussian kernel (described in §3.2.1) and in the decomposition of real-world NIR spectra into elementary Gaussian features (described in §3.2.3).

Let us remind that the principle of NNCoSaMP and NNSP [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF] 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 NNCoSaMP, 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 NNSP is similar, but two NNLS problems (of size roughly 2K and K) are solved before and after thresholding. The NNSP output thus contains at most K nonzero elements depending on the activation of nonnegativity constraints. NNHTP [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF] maintains an estimated support of size K by keeping the K coefficients having the largest magnitudes from a gradient descend move then solving the NNLS problem related to the maintained support. Therefore, NNHTP output also contains at most K nonzero elements.

Sparse deconvolution with Gaussian kernel

We use the same setting as in §3.4.1 and the result is shown in Figure 3.16. As can be seen, the NNOG algorithms yield competitive accuracy performance compared to NNCoSaMP and NNSP while NNHTP yields the weakest performance. 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 NNCoSaMP, NNSP and NNHTP. 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, NNCoSaMP 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. Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, P n = 10 -2 , and a dictionary of size 2400 × 2388.

Decomposition of NIR spectra

For each NIR spectrum, NNOG algorithms, NNCoSaMP, NNSP and NNHTP are run until a support of cardinality K = 20 is reached. The average computing time and normalized approximation error over 50 spectra are shown in Tab. 3.4. NNCoSaMP and NNSP are slower than NNOG algorithms and yield higher errors. NNHTP is at the same cost as NNOLS but with a weaker performance. This result is further illustrated in Fig. 3.17. In this example, SNNOLS and NNOLS outperform NNCoSaMP and NNSP around 9000 cm -1 while NNHTP performance is rather poor. .17: Data approximation and sparse recovery of an NIR spectrum using non-negative greedy algorithms. The Gaussian dictionary contains 2998 atoms and the sparsity level is set to K = 20.

Conclusion

NNOG algorithms yields competitive performances to those of NNCoSaMP, NNSP and NNHTP. However, NNOG algorithms can easily lend themselves to recursive implementation while it is not the case for NNCoSaMP, NNSP and NNHTP.

Comparison with non-negative extension of LARS

In this section we compare NNOG algorithms with non-negative extensions of LARS [START_REF] Efron | Least angle regression[END_REF] named NLARS [START_REF] Morup | Approximate l0 constrained non-negative matrix and tensor factorization[END_REF]. Note that LARS is an 1 solver with a greedy structure. Each iteration of NLARS consists of either a forward move or a backward move and the update of amplitudes along a computed direction. Our comparisons are done in a simulated sparse deconvolution problem with Gaussian kernel (described in §3.2.1) and in the decomposition of real-world NIR spectra into elementary Gaussian features (described in §3.2.3).

Sparse deconvolution with Gaussian kernel

We use the same setting as in §3.4.1 and the result is shown in Figure 3.18. The NNOG algorithms and NLARS yield competitive performances in terms of support recovery and coefficient accuracy. However, NNOG algorithms have a better ability to find a low value of the residual for a given cardinality K while NLARS is faster.

Decomposition of NIR spectra

For each NIR spectrum, NNOG algorithms and NLARS are run until a support of cardinality K = 20 is reached. The average computing time and normalized approximation error over 50 spectra are shown in Tab. 3.5. Similar to the previous comparison ( §3.8.1), NLARS yields a larger approximation error than those of NNOG algorithms. However, in this comparison, NLARS is slower than NNOG algorithms. We noticed that the number of iterations of NLARS is much higher than those of NNOG algorithms (50 iterations in average for NLARS in comparison with 21 for NNOG algorithms). This might come from the fact that the dictionary in this problem is over-complete (with roughly twice columns as rows, hence highly correlated) while the one in §3.8.1 is slightly under-complete. Since, at each iteration of NLARS, the residual is not orthogonal to the selected subspace, NLARS might need several iterations to select and deselect the same atoms before finding a good candidate to append to the current support. Therefore, NLARS can be more costly than NNOG algorithms even though each NLARS iteration might be cheaper than that of NNOG algorithm. This comparison is further illustrated in Fig. 3.19. As can be seen, the approximation resulted by NLARS is much poorer than those of NNOG algorithms.

Conclusion

NNOG algorithms and NLARS yield competitive performances in terms of support recovery and coefficient accuracy. However, NNOG algorithms yield a lower residual value for a given cardinality K. Besides, NLARS tends to be slower when dealing with highly correlated dictionaries.

Conclusion

According to our numerical comparisons, we can conclude that the computing time of NNOMP then becomes comparable to that of OMP, so that the additional cost should not prevent users from skipping from OMP to NNOMP in any 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, we conclude that it is not necessary 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 [Yaghoobi and Davies, 2015]. 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 [Yaghoobi and Davies, 2015]. Our conclusion is that SNNOLS represents a good trade-off between NNOMP and NNOLS: it is structurally simpler than NNOLS and hence significantly faster, with very similar performance in terms of estimation error. However, NNOLS is more robust in dealing with high conditioning.

Our comparisons with NNSP and NNCoSaMP 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 NNSP and NNCoSaMP. An attractive feature of NNOG algorithms is their versatility. Indeed, the implementations presented in this chapter can be used to address both forms of (1) problems with constraints x 0 ≤ K (as addressed in this chapter) 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 12 by r S 2 2 ≤ ε 2 . NNSP, NNCoSaMP, and NNHTP 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].

Finally, our comparison shows that NNOG algorithms and NLARS are very competitive in performance. However, NNOG algorithms have a better ability to find a low value of the residual for a given cardinality K.

Introduction

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]] are well known greedy, iterative algorithms for sparse decomposition. Their common principle is to sequentially select atoms from a given dictionary, that produce a maximal decrease of the residual error in the least square sense. Finding the K atoms that jointly minimize the residual error at a given sparsity level K is an NP-hard subset selection problem [Natarajan, 1995]. Greedy algorithms such as OMP and OLS form a family of approximate schemes, with a relatively low computing cost compared to exact methods [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF], while convex relaxation yields another important branch of approximate methods [Tibshirani, 1996, Donoho andTsaig, 2008].

Orthogonal greedy algorithms are able to achieve exact support recovery under certain conditions. The Mutual Incoherence Property (MIP) condition states that OMP is guaranteed to recover the K-sparse representation of a data signal in the noiseless case if the mutual coherence of the dictionary is less than 1 2K-1 [Tropp, 2004, Herzet et al., 2013]. In the last years, much attention was paid to variations over such a condition. In particular, extensions to noisy cases were considered, see, e.g., [Ben-Haim et al., 2010, Cai andWang, 2011], and relaxations given partial information on the decomposition of the data signal were proposed [START_REF] Herzet | Exact recovery conditions for sparse representations with partial support information[END_REF], Herzet et al., 2016]. Other studies about the theoretical assessment of OMP also consider the Restricted Isometry Property (RIP) condition [START_REF] Davenport | Analysis of orthogonal matching pursuit using the restricted isometry property[END_REF], Satpathi et al., 2013, Li et al., 2015[START_REF]A sharp condition for exact support recovery with orthogonal matching pursuit[END_REF]. Exact support recovery properties of OLS have been studied more recently. The MIP condition was extended to OLS in the noiseless case [START_REF] Herzet | Exact recovery conditions for sparse representations with partial support information[END_REF] as well as noisy case [START_REF] Herzet | Relaxed recovery conditions for OMP/OLS by exploiting both coherence and decay[END_REF], whereas the exact recovery analyses in [Wen et al., 2017a, Wang andLi, 2017] are based on RIP assumptions. In [START_REF] Soussen | Joint kstep analysis of orthogonal matching pursuit and orthogonal least squares[END_REF], Herzet et al., 2013, Herzet et al., 2016], OMP and OLS were treated as two instances of a unique orthogonal greedy scheme, and the generic acronym Oxx was used to refer to both. We will adopt the same convention in this chapter.

In the traditional K-step analysis, exact support recovery holds when each iteration of Oxx necessarily selects an atom in the true support S * , so that S * is recovered after K = |S * | iterations. Our main contribution in this chapter is to unveil an exact sign recovery property regarding the weights of the selected atoms in the best current approximation under the MIP condition and when the noise is bounded in 2 norm: at any iteration k between 1 and K, the k coefficients found by Oxx have the same sign as the k corresponding ones in the true decomposition of the data signal.

As an immediate consequence of this sign preservation property, we can theoretically assess exact support recovery properties for the non-negative version of OMP and other non-negative greedy schemes. Non-negative OMP (NNOMP) was first introduced by Bruckstein et al. [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF]. It relies on the repeated maximization of the positive inner product between the residual error and the dictionary atoms, followed by the resolution of a non negative least square problem. For such a sign constrained version of OMP, existing exact support recovery analyses are rare, and somewhat discordant. On the one hand, Bruckstein et al. conjectured that a MIP type property holds for NNOMP [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF], and that the proof is similar to the one given in [Tropp, 2004, Donoho et al., 2006] for OMP. Specifically, [Bruckstein et al., 2008, Theorem 3] states that µ < 1 2K-1 is a sufficient condition for exact recovery of any K-sparse representation using NNOMP. On the other hand, Kim et al. elaborate a unified MIP analysis of NNOMP and its generalized version in the multiple measurement vector setting [Kim and Haldar, 2016, Theorem 1]. In the specific case of NNOMP, that is, for single measurement vectors, the authors acknowledged that their MIP condition turns out to be very restrictive. Indeed, µ is required to be lower than 1 K-1 -1 2 , which can occur only when K = 1 or 2. In addition, both analyses in [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] and [START_REF] Kim | Greedy algorithms for nonnegativityconstrained simultaneous sparse recovery[END_REF] are made for the noiseless case. To the best of our knowledge, the exact support recovery of non-negative OLS has not been addressed yet and neither does the extension of related analyses to noisy cases.

We have been unable to prove [Bruckstein et al., 2008, Theorem 3] as prescribed by the authors (as a direct extension of the derivations in [Tropp, 2004,Donoho et al., 2006]). The major obstacle is that the NNOMP selection rule performs comparisons between signed inner products, whereas a small mutual coherence condition yields a bound on the unsigned magnitude of inner products (see Section 4.2.3 for further details). Fortunately, the sign preservation property of OMP implies that OMP and NNOMP coincide in the exact support recovery regime, so [Bruckstein et al., 2008, Theorem 3] becomes a trivial byproduct of our sign preservation analysis. Under the same conditions, OLS and NNOLS [Yaghoobi and Davies, 2015] coincide, which allows us to extend [Bruckstein et al., 2008, Theorem 3] to NNOLS. Similarly, the same applies to SNNOLS [Yaghoobi and Davies, 2015]. Furthermore, our analysis does not restrict to the noiseless case but can apply to noisy cases in which the noise is 2 -bounded.

The chapter is organized as follows. Section 4.2 introduces useful notations and briefly recalls known results about orthogonal greedy algorithms and their non-negative versions. Section 4.3 contains our main results on sign preservation of Oxx, as well as their consequences on signconstrained greedy algorithms. The key result, stated as Theorem 4.1, is proved in the same section, most technical steps being postponed in Appendix. Section 4.4 contains a limited set of simulations to illustrate the average behavior of the studied algorithms outside the exact support recovery regime. Finally, conclusion can be found in Section 4.5.

Some notations and background

Notations

In this chapter, for convenience, we denote by x S the subvector indexed by S (and x i being the i-th entry of x). Besides, we denote by P ⊥ S = I m -H S H † S the orthogonal projection onto span(H S ) ⊥ , where I m stands for the identity matrix of size m × m. We place the convention that P ⊥ S = I m when S = ∅. When H S is full column rank, one has a further explicit formulation

H † S = (H t S H S ) -1 H t S . (4.1)
Let us recall that

hS i = P ⊥ S h i (4.2)
and

gS i = hS i / hS i , (4.3)
Later, we also use the generic notation

cS i = hS i (OMP case), gS i (OLS case) (4.4)
to refer to either hS i or gS i depending on the context. We will denote by HS (resp., GS and CS ) the matrix of size m × n formed by gathering all projected atoms hS i (resp., gS i and cS i ). For convenience, the residual vector and the support found by orthogonal greedy algorithms at iteration k will be denoted r k and S k , respectively, with |S k | = k. By extension, we will denote r 0 = y and S 0 = ∅. Whenever unambiguous, we will use the simpler notations r and S.

In our analysis, we make use of the mutual coherence of the dictionary, defined by

µ(H) = max i =j |h t i h j |. (4.5)
This quantity tells us how much the dictionary atoms look alike.

Exact recovery analysis of OMP and OLS

As seen in Chapter 1, OMP and OLS address the following minimization problem:

min x 0 ≤K y -Hx 2 . (4.6)
Both algorithms start with a zero vector solution corresponding to the empty support. Then, at each iteration, a new atom is selected and added to the support. This process is repeated

Algorithm 16: Oxx in K steps. Input: y, H Output: x with x 0 = K 1 Initialization: x ← 0; S ← ∅; 2 repeat 3 r ← y -H S x S ;
4 Choose an atom according to either (4.7) or (4.8);

5 S ← S ∪ { }; 6 x S ← H † S y; 7 until |S| = K;
until the support of cardinality K is reached (see Algorithm 16). OMP and OLS share the same coefficient update and only differ by their selection rule:

OMP: ∈ arg max i / ∈S |r t h i |, (4.7)
OLS:

∈ arg max i / ∈S |r t gS i | (4.8)
where r = P ⊥ S y denotes the current residual vector. Note that (4.7) can also be written as

∈ arg max i / ∈S |r t hS i | (4.9)
since r is orthogonal to span(H S ). These inner product expressions come from the geometrical interpretation of OMP and OLS [START_REF] Blumensath | On the difference between Orthogonal Matching Pursuit and Orthogonal Least Squares[END_REF]. Besides, from an optimization viewpoint, the selection rule of OMP is based on the minimization of rz i h i 2 w.r.t. scalar z i , whereas OLS relies on the minimization of y -H S∪{i} z 2 w.r.t. vector z of length |S| + 1 [START_REF] Blumensath | On the difference between Orthogonal Matching Pursuit and Orthogonal Least Squares[END_REF]. From this point of view, it is obvious that the OLS selection rule is more costly than that of OMP, although recursive update schemes are available, see, e.g., [Chen et al., 1989, Rebollo-Neira and[START_REF] Lowe ; Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF].

Let us state the MIP condition that holds for both OMP and OLS in the noiseless case. Note that under the assumption of Lemma 4.1, the K-sparse representation is unique [Donoho and Elad, 2003, Theorem 7].

Lemma 4.1. [Tropp, 2004, Theorem 3.1-3.5], [Herzet et al., 2013, Theorem 2]. Assume that µ(H) < 1 2K-1 . Let y = Hx * be a K-sparse representation. Then Oxx recover the support of x * in K iterations.

The MIP condition in Lemma 4.1 was extended to noisy cases in which the noise is bounded in 2 norm in [Cai andWang, 2011, Herzet et al., 2016]. We combine these results in a unique statement in Lemma 4.2. One can find back the result in [START_REF] Cai | Orthogonal matching pursuit for sparse signal recovery with noise[END_REF], Theorem 1] when replacing 2 in (4.11) by any upper bound b ≥ 2 . When the non-zero coefficients of x * are sorted in a descending order, the result in [Herzet et al., 2016, Theorem 4] is obtained when replacing (4.11) by and Wang, 2011, Theorem 1], [Herzet et al., 2016, Theorem 4]. Assume that µ(H) < 1 2K-1 . Let y = Hx * + where x * is K-sparse and

2 < 1 -(2K -i)µ(H) 2 |x * i |, ∀ 1 ≤ i ≤ K. (4.10) Lemma 4.2. [Cai
2 < 1 -(2K -1)µ(H) 2 min{|x * i |, x * i = 0}. (4.11)
Then Oxx recover the support of x * in K iterations.

Obviously, Lemma 4.1 is a special case of Lemma 4.2 when = 0. Also, let us emphasize that the mutual coherence in the previous simulations in Chapter 3 does not satisfy the MIP condition hence the Lemmas 4.1-4.2 do not apply. For instance, in the simulation of §3.4.1, µ(H) = 0.94 hence the MIP condition µ(H) < 1 2K-1 only holds for K = 1 and the inequality (4.11) only holds for the noise level P < 10 -7 min{|x * i | 2 , x * i = 0}. Here and throughout this chapter, we consider that in the special case where the Oxx selection rules (4.7)-(4.8) yield several solutions (i.e., two atoms i 1 and i 2 for which |r t cS i 1 | = |r t cS i 2 |) including an atom that does not belong to the support of x * , Oxx makes the wrong decision and hence K-step exact recovery does not occur.

Exact recovery analysis in the non-negative setting

Extension of K-step exact recovery analysis of OMP

Bruckstein et al. [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF] claimed that the K-step exact recovery analysis of NNOMP can be carried out as a straightforward extension of the classical K-step exact recovery analyses of OMP [Tropp, 2004,Donoho et al., 2006]. Here, we argue that this extension does not actually appear to be trivial.

Tropp's reasoning in [Tropp, 2004] consists in lower bounding the absolute value of the inner product |r t h i | between the residual and the correct dictionary atoms based on induced matrix norm identities. Unfortunately, similar lower bounds cannot be directly extended when absolute values are dropped. Donoho et al.'s recursive proof [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF] exploits that for noiseless sparse inputs, the data residual lays in the subspace spanned by the true atoms. So the result at the first iteration (guaranteed selection of a true atom) can be exploited again in the next iterations by replacing the original input signal by the residual. To generalize this reasoning to the non-negative case, one would need to prove that for non-negative sparse decompositions y = Hx * (i.e., inputs laying in the positive span of the true atoms), the residual lays in the positive span of the dictionary atoms as well. This conjecture turns out to be false; if true atoms are selected until iteration k, the residual r = H(x *x(k) ) lays in the subspace spanned by the true atoms, where x(k) denotes the NNOMP iterate. Provided that the latter atoms are linearly independent, it is clear that r lays in their positive span if and only if x(k) ≤ x * . Simple numerical tests show that even for toy problems, x(k) ≤ x * may not hold when x * ≥ 0 and µ(H) < 1 2K-1 : see Fig. . 4.1 for a simple illustration, also § 4.4.2 and Fig. 4.3(a). Therefore, sign-preservation does not apply to the residual vector, hence Donoho et al.'s reasoning cannot be directly extended to the non-negative case.

1 analysis with non-negativity constraints

The theoretical analysis of non-negative versions of Basis Pursuit BP : min turns out to be far simpler than the analysis of orthogonal greedy algorithms. Indeed, it is wellknown that contrary to greedy algorithms, the exact recovery analysis of Basis Pursuit heavily depends on the sign pattern. Fuchs [Fuchs, 2004] proved that when y = Hx * , BP and BPDN (for sufficiently small λ) both have a unique solution under the MIP assumption, and that this solution identifies with x * as long as H S * is full rank and

∀j / ∈ S * , | σ * , H † S * h j | < 1 (4.14) y h 1 h 2 x (1) 1 x (2) 1
x

(2) 2 r where S * := supp(x * ) and σ * := sign(x * S * ) ∈ {-1, 1} K denotes the sign pattern. It is noticeable that the latter condition2 only depends on the sign pattern, but not on the magnitudes of coefficients x * i . Denoting by 1 the "all-ones" vector of size K, condition (4.14) with σ * ← 1 thus guarantees exact recovery of x * for any x * ≥ 0 supported by S * . It follows that x = x * is not only the unique solution of BP/BPDN but also of their non-negative counterpart

min x≥0 y -Hx 2 + λ x 1 , (4.15)
often referred to as the non-negative Garrote in statistics [Miller, 2002]. Indeed, the minimum value of y -Hx 2 + λ x 1 over x ≥ 0 can be reached only for x = x * . Summarizing, the sign preservation property of BP is guaranteed as long as Fuchs condition (4.14) is met.

Extension of K-step exact recovery analysis of 1 homotopy

Homotopy is a popular greedy algorithm dedicated to 1 minimization. It was originally proposed in [START_REF] Efron | Least angle regression[END_REF] in the undercomplete setting, and named modified Least Angle Regression (LARS). It was later renamed homotopy by Donoho and Tsaig, and further analyzed in the overcomplete setting under the MIP assumption in the same paper [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF]. Homotopy aims to solve the BPDN problem (4.13) for a continuum of λ's. The principle of homotopy is to reconstruct the regularization path (defined as the set of solutions for all λ) by solving the 1 problem for gradually decreasing λ's starting from λ = +∞, with the corresponding solution x = 0. Homotopy has a stepwise mechanism with an atom selection or deselection at each iteration. It turns out [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF], Theorem 1] that when µ(H) < 1 2K-1 , K-step recovery of any K-sparse vector x * from noise-free observations y = Hx * is guaranteed. Not only the support but also the magnitudes of x * are found after K iterations (no deselection occurs).

Although homotopy has a greedy structure, a major difference with Oxx algorithms is that homotopy solves an 1 -penalized least square problem at each iteration. Therefore, the related exact recovery analysis significantly differs from that of Oxx. The exact recovery analysis in [Donoho and Tsaig, 2008, Theorem 1] is based on two ingredients: the correct selection of atoms indexed by supp(x * ) and the so-called sign agreement property. The latter ensures that elements from the solution support are never removed, so the deselection event never occurs. Sign agreement is defined as follows: at every iteration k, the homotopy iterate x satisfies sign(x S ) = sign(H t S r), where S = supp(x). Donoho and Tsaig showed that when µ(H) < 1 2K-1 , the magnitudes of the active atoms always increase while λ is decreasing. Since the correct magnitudes are exactly found at iteration K (x = x * ), it follows that their sign is unchanged throughout the iterations. So, a sign preservation result in the 1 case is truly obtained as a byproduct of [Donoho and Tsaig, 2008, Theorem 1]. In the Oxx setting, the sign preservation result will be obtained here using far less straightforward reasoning. Moreover, the sign agreement property "sign(x S ) = sign(H t S r)" does not make sense anymore since H t S r = 0, due to the orthogonality between the residual vector and the selected atoms.

Sign preservation and exact recovery

This section contains our main results concerning sign preservation of Oxx when the exact recovery condition in terms of mutual coherence is met. The cornerstone of our study is Theorem 4.1, while the other results of Subsection 4.3.1 are rather direct consequences. In Subsection 4.3.2, we have decomposed the proof of Theorem 4.1 into distinct steps, most technical elements being postponed in Appendix.

Main results

While the dictionary H satisfies µ(H) < 1 2K-1 , let us define the constant

M x * := 1 -(2K -1)µ(H) 2 min{|x * i |, x * i = 0} (4.16)
for any K-sparse vector x * . The sign preservation property for arbitrary noisy K-sparse representations is stated now.

Theorem 4.1. Assume that µ(H) < 1 2K-1 . Let y = Hx * + be a noisy K-sparse representation with 2 < M x * . Then Oxx recovers the support of x * in K iterations, and at each iteration, the weights of selected atoms are of the same sign as the corresponding magnitudes of x * .

An obvious corollary can be stated in the special case of noisy K-sparse representations with non-negative weights.

Corollary 4.1. Assume that µ(H) < 1 2K-1 . Let y = Hx * + be a noisy K-sparse representation with x * ≥ 0 and 2 < M x * . Then, Oxx recovers the support of x * in K iterations, and at each iteration, the weights of selected atoms are positive.

In the non-negative setting, Corollary 4.1 has interesting implications concerning non-negative versions of Oxx. Let us start with the following lemma. Lemma 4.3. Let y ∈ R m and assume that any combination of K dictionary columns is linearly independent.

• Assume that at every iteration k = 1, . . . , K of OMP, the selection rule yields a unique optimal index and that the weights of selected atoms are all positive. Then NNOMP provides the same iterates as OMP (i.e., with the same support and the same weights at every iteration).

• The same applies if OMP and NNOMP are replaced by OLS and SNNOLS, respectively.

• It also applies if OMP and NNOMP are replaced by OLS and NNOLS, respectively.

This result is intuitive since at each iteration of Oxx, the weights of selected atoms form a vector defined as a least squares solution. Clearly, if an unconstrained least squares solution is positive, then it is also the solution of the corresponding NNLS problem.

Proof. See Appendix 4.B.2. From Corollary 4.1 and Lemma 4.3, we can deduce the following result.

Corollary 4.2. Assume that µ(H) < 1 2K-1 . Let y = Hx * + be a noisy K-sparse representation with x * ≥ 0 and 2 < M x * . Then NNOMP identifies with OMP whereas both NNOLS and SNNOLS identify with OLS. Thus, NNOMP, NNOLS and SNNOLS all recover the support of x * in K iterations.

By Proposition 2.3 in Chapter 2 one can deduce that Sparse NNLS [START_REF] Peharz | Sparse nonnegative matrix factorization with l0 constraints[END_REF] identifies with NNOMP and OMP under the assumptions of Corollary 4.2. Therefore, we have the following result regarding the exact recovery of Sparse NNLS.

Corollary 4.3. Assume that µ(H) < 1 2K-1 . Let y = Hx * + be a noisy K-sparse representation with x * ≥ 0 and 2 < M x * . Then Sparse NNLS identifies with OMP and recovers the support of x * in K iterations.

It is known that the MIP condition µ(H) < 1 2K-1 is not only sufficient but also necessary for uniform (i.e., irrespective of the magnitudes of nonzero coefficients in the sparse representation and of the choice of the dictionary) K-step exact support recovery [START_REF] Cai | Stable recovery of sparse signals using an oracle inequality[END_REF],Herzet et al., 2013] by Oxx. [START_REF] Cai | Stable recovery of sparse signals using an oracle inequality[END_REF] indeed exhibited an equiangular dictionary whose mutual coherence equals µ(H) = 1 2K-1 and a vector y having two K-sparse representations Hx = Hz with distinct supports. K-step exact support recovery does not make sense anymore in this situation, since either the support of x or z cannot be reconstructed in K steps. The same analysis can be made concerning non-negative extensions of Oxx.

Corollary 4.4. The MIP condition µ(H) < 1 2K-1 is necessary for K-step exact recovery of any non-negative K-sparse vector by non-negative orthogonal greedy algorithms, since there exists a dictionary H with µ(H) = 1 2K-1 and a vector y having two non-negative K-sparse representations with distinct supports.

Proof. Consider the dictionary H ∈ R m×n with µ(H) = 1 2K-1 and the vector y ∈ R m proposed in [Cai et al., 2010, Section III], the latter vector having two K-sparse representations Hx = Hz with distinct supports. Since both supports are distinct, one can define the subrogate dictionary H ∈ R m×n as:

h i =    sign(x i ) h i if i ∈ supp(x), sign(z i ) h i if i ∈ supp(z), h i otherwise.

Proof of Theorem 4.1

Let us remark that it is sufficient to consider the case of non-negative weights. Indeed, an arbitrary K-sparse representation can be reduced to the case with non-negative weights, by an obvious transformation where each negative weight is replaced by its opposite value, the corresponding atom being also replaced by the opposite one. Moreover, it is straightforward to check that the list of atoms selected by Oxx is invariant through such a transformation. We therefore simply need to prove the result for non-negative weights, which corresponds to the setting of Corollary 4.1. Before going further, let us denote by

C + K, := {Hx * + , x * ≥ 0, x * 0 = K and 2 < M x * } (4.18)
the set of noisy K-sparse representations with non-negative weights and 2 -bounded noise. Now, consider y = Hx * + ∈ C + K, . From Lemma 4.2, Oxx recovers the support S * of x * in K iterations. At any iteration k ≤ K, the support S k of the current solution x(k) is therefore a subset of S * . Recall that x(k) S k is the unconstrained least squares solution related to support S k , see Algorithm 16, line 6. Let r k = P ⊥ S k y denote the related residual vector, with r 0 = y corresponding to S 0 = ∅.

We proceed in two distinct steps to prove that x(k) S k > 0 for all k ∈ {1, . . . , K}. First, we prove that the weight of each newly selected atom

x(k) S k \S k-1 is positive for any k ≤ K. Then, we show that the updated coefficients x(k) S k-1 remain positive. Let us remark that Theorem 4.1 states a trivial fact at iteration K, since x(K) = x * according to Lemma 4.1.

Let us first characterize the last kj coefficients of x(k) S k , j < k being an arbitrary index in the following lemma. Note that this lemma applies to any real vector y. Lemma 4.4. Let y ∈ R m , and let j and k be two iteration indices such that 0 ≤ j < k. Assume that H S k is full column rank. Then, the k-th iterate of Oxx satisfies The previous section shows that in some specifically favorable situations, greedy algorithms such as OMP have not only the capacity to recover the support of the true solution, but also to recover the correct weight signs. In such conditions, according to Lemma 4.3, implementing nonnegative (or more generally, sign-aware) versions of such greedy algorithms is useless. However, one cannot generalize such a conclusion to more realistic scenarios. On the contrary, one can empirically observe that sign-aware greedy algorithms tend to reach superior performance, which is in agreement with the fact that they exploit more information than usual greedy algorithms. To illustrate this fact, let us consider a dictionary H with 22 regularly spaced, discretized Gaussianshaped atoms, with a constant standard deviation σ = 0.5, and a mutual coherence µ(H) = 0.37.

x(k) S k \S j = HS j S k \S j † r j . ( 4
We randomly choose K = 10 atoms in H, whose location in the dictionary are distributed according to the uniform distribution. The atoms are all equally weighted with a unit weight, and we generate data y = Hx * + with some additive Gaussian noise . Note that µ(H) > 1 2K-1 , which implies that exact support recovery is not mathematically guaranteed even at low noise. Within a certain range of signal-to-noise ratio (SNR, defined by 10 log 10 ( Hx * 2 / 2 )), we have generated average performance for OMP, NNOMP, OLS, SNNOLS and NNOLS in terms of rate of exact support recovery, the stopping condition being that the cardinality of the estimated support should be 10. This experiment has been repeated 1000 times to obtain the average results shown in Fig. 4.2. Several empirical conclusions can be drawn. Some of them are already acknowledged facts. For instance, greedy algorithms keep some exact recovery capacities far beyond the "zero defect" area. In contrast, in the low SNR regime, the exact recovery capacity almost surely vanishes. The intermediate zone is the most interesting. Specifically, one can notice that there is a significant difference of performance between the usual greedy algorithms and their non-negative extensions. We have also performed a sign-preservation test for OMP and OLS that simply consists in checking whether at all iterations, all estimated weights are positive. Fig. 4.2 shows that such a sign-preservation property is rather robust, and that empirically, sign-preservation is always reached with Oxx whenever exact support recovery is found.

Non-monotony of the magnitude variations

As argued in § 4.2.3.3, 1 homotopy is a stepwise greedy algorithm for which the sign preservation guarantee holds whenever µ(H) < 1 2K-1 . In [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF], Donoho and Tsaig proved a stronger result under the mutual coherence condition: the magnitudes x(k) i of the selected atoms keep increasing while k is increasing. In contrast, such monotony property does not hold for Oxx algorithms, since the magnitudes may either increase or decrease during the iterations. This claim can be proven analytically at iteration 2, by using the fact that

x(1) S 1 > 0 and

x(2) S 2 > 0 according to Corollary 4.1. Since x(1) and x(2) have closed-form expressions, a simple calculation (skipped for brevity reasons) shows that the first magnitude is decreasing, i.e.,

x(1) S 1x(2) S 1 > 0 if and only if the inner product between the atoms selected in the first two iterations is positive. So, the magnitude of xS 1 may either increase or decrease depending on the sign of the inner product. We further compare 1 -homotopy with the Oxx algorithms for a toy problem of dimension (m, n) = (5, 5), with an equiangular dictionary H such that µ(H) = 0.9 2K-1 . The columns of H satisfy h i , h j = ±µ(H), where the sign of the inner product is randomly chosen. The ground truth vector x * is K-sparse with K = 4, with nonzero magnitudes drawn from the uniform distribution U([0.6, 1]). Since µ(H) < 1 2K-1 , K-step exact recovery holds for all considered algorithms. In Fig. 4.3, the variation of each entry

x(k) i with respect to k is represented with a specific color. As expected (since exact support recovery holds), the black magnitude, which corresponds to the wrong atom i / ∈ S * , is equal to 0 throughout the iterations.

1 -homotopy yields magnitudes that are indeed increasing with k, which is consistent with the theoretical result in [START_REF] Donoho | Fast solution of 1 -norm minimization problems when the solution may be sparse[END_REF]. On the contrary, the OMP and OLS iterates (which are identical here; the same indices are selected at each iteration) are non-monotonous.

Conclusion

We have established that OMP and OLS satisfy the sign preservation property when dealing with sufficiently incoherent dictionaries: in the exact recovery regime µ(H) < 1 2K-1 , the Oxx estimates are guaranteed to keep the same sign as the ground-truth sparse vector, at all iterations. This interesting property enables us to establish the first K-step recovery analysis of three nonnegative greedy algorithms proposed previously, namely NNOMP, NNOLS and SNNOLS, under the MIP condition µ(H) < 1 2K-1 . This exact recovery condition turns out to be identical for both Oxx algorithms and their non-negative extensions. Moreover, it is not only a sufficient but also a (worst-case) necessary condition of exact recovery. Therefore, one cannot distinguish the performance of Oxx and their nonnegative counterparts from this theoretical viewpoint. There is still room for improvement though in cases where the dictionary atoms are known to be non-negative valued, since then the MIP condition is not guaranteed to be necessary anymore. This setting covers many application fields ranging from sparse deconvolution [Bendory, 2017] to tomographic image reconstruction [START_REF] Petra | Average case recovery analysis of tomographic compressive sensing[END_REF].

Empirically, we observed that the sign preservation property holds for noiseless scenarios when µ(H) is substantially greater than 1 2K-1 . However, for truly coherent dictionaries and for noisy observations, Oxx algorithms do yield iterates with negative entries, thus sign-aware versions are worth being considered. Furthermore, the latter versions yield improvement of empirical exact recovery performance for random dictionaries when µ(H) > 1 2K-1 (see Chapter 3). A perspective of this work is to elaborate on exact recovery conditions when µ(H) ≥ 1 2K-1 . We argued that the MIP condition µ(H) < 1 2K-1 is both necessary and sufficient for uniform recovery of K-sparse signals, that is, irrespective of the magnitudes of the nonzero coefficients in the K-sparse representation. It is well-known that such worst case analysis is pessimistic since practically, algorithms may succeed far beyond the exact recovery regime µ(H) < 1 2K-1 . It was shown in [START_REF] Herzet | Relaxed recovery conditions for OMP/OLS by exploiting both coherence and decay[END_REF] that exact recovery guarantees can be obtained for Oxx when µ(H) ∈ [ 1 2K-1 , 1 K ) provided that the decay of the magnitudes of the nonzero coefficients is fast enough (the minimum rate of decay required to ensure exact support recovery depends on µ(H)). It would then be interesting to generalize our sign preservation analysis to the case µ(H) ∈ [ 1 2K-1 , 1 K ). However, this extension does not appear to be obvious and is left for future work. With the same idea of elaborating on sign preservation under weaker K-step support recovery conditions, Tropp's exact recovery condition (ERC) is worth being considered. Indeed, it is both necessary and sufficient for uniform K-step exact recovery of all representations having a given support S * with both OMP and OLS [Tropp, 2004,Soussen et al., 2013]. In all our empirical tests, we found that sign preservation holds whenever the ERC is met. However, proving this conjecture would necessitate more involved theoretical analysis than those we derived in this document.

Appendix 4.A Useful lemmas

Let us recall some useful lemmas. Lemma 4.7 provides an upper bound on the 1 -norm of the Gram matrix columns by means of mutual coherence. Lemma 4.8 provides lower and upper bounds on the inner product of projected atoms. Lemma 4.9 is related to the full rankness of the matrix of projected atoms. Lemma 4.10 is a simple algebraic manipulation related to the pseudo-inverse. Lemma 4.7. [Foucart and Rauhut, 2013, Theorem 5.3], [Tropp, 2004, Proposition 2.1, Theorem 3.5]. If B is a column-normalized matrix with k columns and µ(B) < 1 k-1 then B is full column rank and

(B t B) -1 1,1 ≤ 1 1 -(k -1)µ(B) (4.20)
wherein • 1,1 equals the maximum absolute column sum of its argument.

Lemma 4.8. [Herzet et al., 2016, Lemmas 2-3]

. If µ(H) ≤ 1 k+1 with k = |S| then ∀i / ∈ S, hS i 2 ≥ β k , (4.21) ∀p = q, |(c S p ) t hS q | ≤ µ k η k , (4.22)
where

β k = (µ(H) + 1)(1 -kµ(H)) 1 -(k -1)µ(H) (4.23) µ k = µ(H) 1 -kµ(H) (4.24) η k = β k for OMP, √ β k for OLS. ( 4

.25)

Lemma 4.9. [Soussen et al., 2013, Lemma 8] If S ∩ S = ∅ and H S∪S is full column rank, then matrices HS S and GS S are full column rank.

Lemma 4.10. Consider a full column rank matrix H = h 1 , H 2 ∈ R p×q where h 1 ∈ R p and H 2 ∈ R p×(q-1) denotes the submatrix formed of the last q -1 column vectors. Then,

∀r ∈ R p , H † r 1 = r, P ⊥ S 2 h 1 P ⊥ S 2 h 1 2 . (4.26)
where the index set S 2 = {2, . . . , q} corresponds to the columns of H 2 .

Proof. Let r ∈ R p . Since H is full column rank, r can be uniquely decomposed as Let us now start with the proof of Lemma 4.4, since this Lemma will be used later in the proofs of Lemmas 4.3,4.5,and 4.6. Since x(k) S k is the unconstrained least squares solution related to subset S k , we have for j < k,

r = p H + p H ⊥ = (H † r) 1 h 1 + r H 2 + p H ⊥ ( 4
x(k) S k = arg min z∈R k y -H S k z 2 = arg min v,w y -H S j v -H S k \S j w 2 . (4.29)
In addition,

min v,w y -H S j v -H S k \S j w 2 = min w min v (y -H S k \S j w) -H S j v 2 = min w P ⊥ S j (y -H S k \S j w) 2 = min w r j - HS j S k \S j w 2 . ( 4

4.B.2 Proof of Lemma 4.3

We prove by induction that the supports found by OMP and NNOMP (respectively, by OLS and NNOLS/SNNOLS) coincide at each iteration.

4.B.2.1 NNOMP vs OMP

Let us first prove the claim at the first iteration. Let us denote by the index of the atom selected by OMP. The first OMP iterate is the one-sparse vector supported by { } such that x(1) = h t y. By assumption, x(1) > 0. Since is the unique solution to (4.7) with r 0 = y, we have:

{ } = arg max i {h t i y}, (4.32)
which implies that is also selected by NNOMP. Since h t y > 0, x(1) ≥ 0 corresponds to both the unconstrained and non-negative solutions related to subset { }, that is, to both NNOMP and OMP very first iterates.

Let us assume that OMP and NNOMP deliver the same support S k-1 and solution after k -1 iterations. The residual vectors r k-1 corresponding to OMP and NNOMP therefore coincide. Let S k = S k-1 ∪ { } denote the support found by OMP at iteration k with the index of the atom selected at iteration k, and let x(k) denote the OMP iterate. Lemma 4.4 with j ← k -1 implies that: span(H S k-1 ). By assumption, x(k) ≥ 0, so (4.34) implies that r t k-1 h ≥ 0. Since is the unique solution to (4.7), we have:

x(k) = hS k-1 † r k-1 = r t k-1 hS k-1 hS k-1 2 (4.33) = r t k-1 h hS k-1 2 . ( 4
{ } = arg max i / ∈S k-1 {r t k-1 h i }. (4.35)
So, is also selected by NNOMP at iteration k, leading to the same subset S k as for OMP.

Because the OMP iterate is non-negative, it is also the NNLS solution corresponding to S k . Hence, NNOMP yields the same iterate as OMP.

4.B.2.2 SNNOLS vs OLS

We can make a similar argument as in the previous case by replacing the dictionary atoms h i by their normalized projections gS i . At the first iteration, g∅ i = h∅ i = h i , so the very first iterates of SNNOLS and OLS respectively identify with those of NNOMP and OMP. They coincide according to § 4.B.2.1. At iteration k, the proof of § 4.B.2.1 can be repeated, where x(k) now denotes the OLS iterate. The right-hand side of (4.33) rereads r

t k-1 gS k-1 / hS k-1 ≥ 0, so r t k-1 gS k-1 ≥ 0.
Finally, since is the only maximizer of the OLS selection rule (4.8), we have

{ } = arg max i / ∈S k-1 r t k-1 gS k-1 i . (4.36)
So, is also selected by SNNOLS. Similar to the NNOMP vs OMP case, we conclude that NNOLS and SNNOLS yield the same iterate at iteration k.

4.B.2.3 NNOLS vs OLS

The very first iterates of OLS and NNOLS identify to those of OMP and NNOMP, respectively. We have proved above that they coincide. Assume that OLS and NNOLS deliver the same support S k-1 after k -1 iterations, and the same iterate

x(k) S k-1 = H † S k-1 y > 0.
Let denote the atom selected by OLS at iteration k, with S k = S k-1 ∪ { }. Since OLS selects the atom yielding the least squared error and the optimal index is unique, we have

∀i / ∈ S k , min z y -H S k z 2 < min z y -H S k-1 ∪{i} z 2 . (4.37) Clearly, min z y -H S k-1 ∪{i} z 2 ≤ min z≥0 y -H S k-1 ∪{i} z 2 . ( 4
y -H S k z 2 < min z≥0 y -H S k-1 ∪{i} z 2 ,
which implies that is also selected by NNOLS. By assumption, H S k is full column rank, so the unconstrained and non-negative least-squares solutions related to S k are unique. (4.39) implies that they coincide. Hence, OLS and NNOLS yield the same iterate at iteration k.

4.B.3 Proof of Lemma 4.5

In this proof, the abridged notations ỹ, hi , gi , ci , ˜ correspond to projected vectors onto span(H S k-1 ) ⊥ . In a somewhat counter-intuitive manner, our proof of Lemma 4.5 is not recursive: to prove that x(k) S k \S k-1 > 0, we exploit that good atoms have been selected at previous iterations, but the current weight signs are not taken into consideration. We denote by the atom selected at iteration k, so that S k \S k-1 = { }.

Let us first consider the case K = 1. Without loss of generality, we assume that y = xh 1 + where x > 0. From Lemma 4.2, h 1 is selected at the first iteration. for i = 1 and i = p. Now h t 1 y + h t p y > 0 implies (4.41) so the case K = 1 is proved. Hereafter we consider the case where K ≥ 2. Also, it should be noticed that µ(H) < 1 2K-1 implies that H S * is full column rank according to Lemma 4.7. Applying Lemma 4.4 for j ← k -1, we have

y = x + xh t 1 h p + t h 1 + t h p (4.42) ≥ x -µ(H)x -| t h 1 | -| t h p | (4.43) ≥ x(1 -µ(H)) -2 (
x(k) S k \S k-1 = ht r k-1 . (4.47)
Therefore, x(k) S k \S k-1 has the same sign as ct r k-1 since c = h or h / h ( h = 0 according to Lemma 4.9). The remaining part of the proof consists in showing that ct r k-1 > 0.

Since h is selected at the k-th iteration, we have

∀j / ∈ S k , |c t r k-1 | ≥ |c t j r k-1 |. (4.48)
To prove that ct r k-1 > 0, we are going to exhibit an index p / ∈ S k such that Let us introduce the decomposition

ct r k-1 > -|c t p r k-1 |. ( 4
r k-1 = P ⊥ S k-1 y = i / ∈S k-1 x * i hi + ˜ . (4.50)
For any j ∈ S * \S k-1 , we deduce

ct j r k-1 = x * j ht j cj + i / ∈S k-1 ∪{j} x * i ht i cj + ˜ t cj . ( 4
ht j cj ≥ η k-1 , (4.53)
| ht i cj | ≤ µ k-1 η k-1 (4.54)
with µ k-1 and η k-1 defined in (4.24)- (4.25). By Cauchy-Schwarz inequality,

|˜ t cj | ≤ ˜ cj ≤ ˜ (4.55) as cj ≤ 1. From (4.51), we get for j ∈ S * \S k-1 , ct j r k-1 ≥ η k-1 x * j -µ k-1 i / ∈S k-1 ∪{j} x * i -˜ . (4.56)
Notice that there are Kk nonzero terms in the sum in (4.56) since x * is K-sparse. The latter sum can then be upper bounded as:

i / ∈S k-1 ∪{j} x * i ≤ x * + (K -k -1)x * p if j = , (K -k)x * p if j = (4.57)
with p defined as:

p ∈ arg max i∈S * \S k x * i . (4.58)
Apply (4.56) twice with j ← p and j ← . We get

ct p r k-1 ≥ η k-1 (1 -(K -k -1)µ k-1 )x * p -µ k-1 x * -˜ , ct r k-1 ≥ η k-1 x * -(K -k)µ k-1 x * p -˜ , from which we deduce ct p r k-1 + ct r k-1 ≥ η k-1 (τ 1 x * + τ 2 x * p ) -2 ˜ , (4.59)
with

τ 1 := 1 -µ k-1 , (4.60) τ 2 := 1 -(2K -2k -1)µ k-1 . (4.61) According to (4.24), τ 1 > 0 if µ(H) < 1 k and τ 2 > 0 if (2K -k -2)µ(H) < 1. Since µ(H) < 1 2K-1 , it holds true that τ 1 > 0 and τ 2 > 0. Let us define q ∈ arg min i∈S * x * i . (4.62) Then η k-1 (τ 1 x * + τ 2 x * p ) -2 ˜ ≥ η k-1 (τ 1 + τ 2 )x * q -2 ˜ . (4.63) According to (4.21), β k-1 ≤ 1 which implies β k-1 ≤ β k-1 ≤ 1 hence η k-1 ≥ β k-1 by definition. Besides, ˜ ≤ < 1 -(2K -1)µ(H) 2 x * q . (4.64) Therefore, η k-1 (τ 1 + τ 2 )x * q -2 ˜ > β k-1 (τ 1 + τ 2 )x * q -(1 -(2K -1)µ(H)) x * q . ( 4
β k-1 (τ 1 + τ 2 ) = 2β k-1 (1 -(K -k)µ k-1 ) (4.66) = 2(1 + µ(H))(1 -(k -1)µ(H)) 1 -(k -2)µ(H) 1 - (K -k)µ(H) 1 -(k -1)µ(H) (4.67) = 2(1 + µ(H))(1 -(k -1)µ(H)) 1 -(k -2)µ(H) 1 -(K -1)µ(H) 1 -(k -1)µ(H) (4.68) = 2(1 + µ(H))(1 -(K -1)µ(H)) 1 -(k -2)µ(H) . (4.69) As k ≥ 1, we have 1 -(k -2)µ(H) ≤ 1 + µ(H) so 2(1 + µ(H))(1 -(K -1)µ(H)) 1 -(k -2)µ(H) ≥ 2(1 -(K -1)µ(H)). (4.70) As a result, β k-1 (τ 1 + τ 2 ) ≥ 2(1 -(K -1)µ(H)). (4.71)
From (4.65) and ( 4.71) we have

η k-1 (τ 1 + τ 2 )x * q -2 ˜ > 2(1 -(K -1)µ(H))x * q -(1 -(2K -1)µ(H))x * q (4.72) = (1 + µ(H))x * q (4.73) > 0. (4.74)
Thus, from (4.63) we get η k-1 (τ 1 x * + τ 2 x * p ) -2 ˜ > 0. Hence, from (4.59) we get ct p r k-1 + ct r k-1 > 0, which implies (4.49). We thus conclude that ct r k-1 and

x(k) S k \S k-1 are positive quantities.

4.B.4 Proof of Lemma 4.6

In this proof, the abridged notations hi , gi , ci correspond to projected vectors onto span(H S j-1 ) ⊥ (the latter space being different from that of Subsection 4.B.3). Similarly, we will use the abridged notations H, G, C to refer to HS j-1 , GS j-1 and CS j-1 , respectively. For any subset S, the Gram matrices of HS and GS will be respectively denoted by Gr OMP S := Ht S HS and Gr OLS S := Gt S GS . For statements that are common to both OMP and OLS, we will use the simpler generic notation Gr S := Ct S CS . The proof of Lemma 4.6 is not recursive. However, in order to prove that x(k) S k-1 > 0, we use the fact that x(j) S j \S j-1 > 0 for j < k, which holds according to Lemma 4.5. Let ∈ S k-1 and denote by j < k the iteration at which the atom indexed by has been selected by Oxx, so that S j \S j-1 = { }. According to Lemma 4.4 and since is the first entry in the ordered set S k \S j-1 , we have

x(k) = H † S k \S j-1 r j-1 1 . ( 4
\S j-1 = CS k \S j-1 ∆,
where ∆ > 0 is a square diagonal matrix whose diagonal elements are either equal to 1 (OMP,

C ← H) or to hi > 0, i ∈ S k \S j-1 (OLS, C ← G).
The positivity property is related to the full rankness of CS k \S j-1 , which deduces from that of H S k according to Lemma 4.9. From standard properties of pseudo-inverses, we have H †

S k \S j-1 = ∆ -1 C † S k \S j-1 . (4.75) yields x(k) ∝ + C † S k \S j-1 r j-1 1 = c , CS k \S j † r j-1 1 (4.76) = r j-1 , P ⊥ S k \S j c P ⊥ S k \S j c 2 (4.77)
where ∝ + indicates proportionality up to a positive factor, P ⊥ S k \S j denotes the orthogonal projection onto the orthogonal complement of span( CS k \S j ), and (4.77) deduces from Lemma 4.10.

(4.77) implies that:

x(k) > 0 ⇐⇒ P ⊥ S k \S j r j-1 , c > 0 (4.78) ⇐⇒ CS k \S j C † S k \S j r j-1 , c < r j-1 , c (4.79) ⇐⇒ C † S k \S j r j-1 , Ct S k \S j c < r j-1 , c . ( 4
S k \S j c ∞ C † S k \S j r j-1 1 ≤ Ct S k \S j c ∞ Gr -1 S k \S j 1,1 Ct S k \S j r j-1 1 ≤ µ( C) Gr -1 S k \S j 1,1 Ct S k \S j r j-1 1 ≤ µ( C) Gr -1 S k \S j 1,1 (k -j) | r j-1 , c |. (4.81)
To obtain the last inequality, we exploit that h has been selected at the j-th iteration of Oxx, therefore:

∀i ∈ S k \S j , | r j-1 , ci | ≤ | r j-1 , c |. (4.82)
From Lemma 4.5, we have x(j) > 0, and as already remarked in the proof of Lemma 4.5, ct r j-1 is of same sign as x(j) , see (4.47). The upper bound in (4.81) thus rewrites

(k -j) µ( C) Gr -1 S k \S j 1,1 r j-1 , c (4.83)
and we deduce from (4.80) that

(k -j)µ( C) Gr -1 S k \S j 1,1 < 1 =⇒ x(k) > 0 . ( 4

.84)

Let us now provide some upper bounds on µ( C) and Gr -1 S k \S j 1,1 in order to show that the LHS of (4.84) holds true.

4.B.4.1 Upper bound on

µ( C) Since µ(H) < 1 2K-1 < 1 j , Lemma 4.8 yields: ∀i / ∈ S j-1 , hi 2 ≥ β j-1 , (4.85) ∀p = q, | ht p hq | ≤ µ j-1 β j-1 , (4.86)
with µ j-1 and β j-1 defined in (4.23)-( 4.24), from which we can easily deduce that

µ( C) = µ( H) ≤ µ j-1 β j-1 (OMP case), µ( G) ≤ µ j-1 (OLS case). ( 4
µ( G) ≤ µ j-1 ≤ 1 2K -j ≤ 1 k -j . ( 4
.88) Lemma 4.7 then applies to matrix GS k \S j ∈ R m×(k-j) , so the latter is full column rank, and

Gr OLS S k \S j -1 1,1 ≤ 1 1 -(k -j -1)µ( GS k \S j ) ≤ 1 1 -(k -j -1)µ( G) . (4.89) It follows from (4.87) that µ( G) Gr OLS S k \S j -1 1,1 ≤ µ j-1 Gr OLS S k \S j -1 1,1 ≤ µ j-1 1 -(k -j -1)µ( G) , ≤ 1 2K -k + 1 (4.90) < 1 k -j (4.91)
where (4.90) follows from the second upper bound in (4.88). ( 4.91) thus implies that the LHS of (4.84) is true in the OLS case.

4.B.4.3 Upper bound on

Gr -1 S k \S j 1,1 in the OMP case Contrary to the OLS case, the diagonal elements of Gr OMP S k \S j ( hi 2 , i ∈ S k \S j ) are not equal to 1, so Lemma 4.7 does not apply. Let ∆ be the square diagonal matrix with the elements hi , i ∈ S k \S j on its diagonal. Clearly, HS k \S j = GS k \S j ∆, hence

Gr OMP S k \S j -1 = ∆ -1 Gr OLS S k \S j -1 ∆ -1 , (4.92)
and thus

Gr OMP S k \S j -1 1,1 ≤ Gr OLS S k \S j -1 1,1 ∆ -1 2 1,1 . (4.93) Moreover, ∆ -1 2 1,1 = 1 min i∈S k \S j hi 2 ≤ 1 β j-1 (4.94)
by (4.85). We have thus

µ( H) Gr OMP S k \S j -1 1,1 ≤ µ( H) β j-1 Gr OLS S k \S j -1 1,1 ≤ µ j-1 Gr OLS S k \S j -1 1,1
according to (4.87). From (4.91), we conclude that the LHS of (4.84) is also true in the OMP case.

Chapter 5

NP-hardness of 0 minimization problems In chapters 2-4, we addressed the algorithmic issues and theoretical analysis of NNOG algorithms. In particular, we exhibited structural properties and proposed a fast implementation of NNOG algorithms. Moreover, we established an exact support recovery analysis of NNOG algorithms under the standard mutual coherence based condition. Note that NNOG algorithms are designed to address the constrained non-negative 0 minimization problems. In this chapter, we discuss about NP-hardness of non-negative 0 minimization problems.

Introduction

Sparse approximation appears in a wide range of applications, especially in signal processing, image processing and compressed sensing [Elad, 2010]. Given a signal data y ∈ R m and a dictionary H of size m × n, the aim is to find a signal x ∈ R n that gives the best approximation y ≈ Hx and has the fewest non-zero coefficients (i.e.,, sparsest solution). This task leads to solving one of the following constrained or penalized 0 minimization problems:

min x x 0 s.t. y -Hx 2 ≤ ( 0 C ) min x y -Hx 2 2 s.t. x 0 ≤ K ( 0 C K ) min x y -Hx 2 2 + λ x 0 ( 0 P )
in which , K and λ are positive quantities related to the noise standard deviation, the sparsity level and regularization strength, respectively. Letters C and P respectively indicate that the problem is constrained or penalized. Depending on application, the appropriate statement will be addressed. It is noteworthy that n and K often depend on m when one considers the size of problem. ( 0 C ) and ( 0 C K ) are well known to be NP-hard [Natarajan, 1995, Davis et al., 1997]. The NP-hardness of ( 0 P ) was claimed to be a particular case of more general complexity analyses in [Chen et al., 2014, Huo andChen, 2010]. However, we point out hereafter that these complexity analyses do not rigorously apply to ( 0 P ) as claimed. Then we provide a new proof for the NP-hardness of ( 0 P ) adapted from Natarajan's construction [Natarajan, 1995].

In several applications the signal or image of interest is both sparse and non-negative. Therefore, many researchers are interested in non-negative 0 minimization problems as follows:

min x≥0 x 0 s.t. y -Hx 2 ≤ ( 0 C +) min x≥0 y -Hx 2 2 s.t. x 0 ≤ K ( 0 C K +) min x≥0 y -Hx 2 2 + λ x 0 ( 0 P +)
Several papers address non-negative 0 minimization problems in the literature (see, e.g.,, [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF],Peharz and Pernkopf, 2012,Yaghoobi et al., 2015,Wang et al., 2018]). However, to the best of our knowledge, the complexity of these problems has not been addressed yet, the question of their NP-hardness being still open. Here we show that these problems are NP-hard and the proof can be derived from the NP-hardness of 0 minimization problems.

Background on NP-hardness and constrained 0 minimization problems

Let us recall that an NP-complete problem is a problem in NP to which any other problem in NP can be reduced in polynomial time. Thus NP-complete problems are identified as the hardest problems in NP. An NP-complete problem is strongly NP-complete if it remains NP-complete when all of its numerical parameters are bounded by a polynomial in the length of the input. NP-hard problems are at least as hard as NP-complete problems. However, NP-hard problems do not need to be in NP and do not need to be decision problems. Formally, a problem is NP-hard (respectively, strongly NP-hard) if a NP-complete (respectively, strongly NP-complete) problem can be reduced in polynomial time to it. The reader is referred to [Garey andJohnson, 1979, Leeuwen, 1990] for more information on this topic.

In the literature, problem ( 0 C ), called SAS in [Natarajan, 1995], is well known to be NPhard [Natarajan, 1995, Theorem 1]. The NP-hardness of ( 0 C ) is a valuable extension of an earlier result: the problem of minimum weight solution to linear equations (equivalent to ( 0 C ) with = 0) is NP-hard [Garey and Johnson, 1979, p. 246]. Davis et al. proved that ( 0 C K ), called M -optimal approximation in [START_REF] Davis | Adaptive greedy approximations[END_REF], is NP-hard for any K < m [Davis et al., 1997, Theorem 2.1]. Both analyses of Natarajan and Davis were made by a polynomial time reduction from the "exact cover by 3-sets" problem2 which is known to be NP-complete [Garey and Johnson, 1979, p. 221].

Existing analyses on penalized 0 minimization

In [Chen et al., 2014, Huo andChen, 2010], the NP-hardness of ( 0 P ) is deduced as a particular case of more general complexity analyses. However, it turns out that the latter do not apply to ( 0 P ), as explained hereafter. Chen et al. [START_REF] Chen | Complexity of unconstrained l2-lp minimization[END_REF] address the unconstrained qp minimization problem, defined by: min x y -Hx q q + λ x p p ( qp )

where λ > 0, q ≥ 1 and 0 ≤ p < 1. The authors showed that problem ( qp ) is NP-hard with any λ > 0, q ≥ 1 and 0 ≤ p < 1 [Chen et al., 2014, Theorem 3]. Obviously, ( 0 P ) is the case where q = 2 and p = 0. The proof was done by i) introducing an invertible transformation which scales any instance of problem ( qp ) to the problem ( qp ) with λ = 1/2, and ii) establishing a polynomial time reduction from the partition problem which is known to be NP-complete [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF] to the problem ( qp ) with λ = 1/2. In other words, they showed that problem ( qp ) with λ = 1/2 is NP-hard and, because there exists an invertible transformation from any problem ( qp ) to the one with λ = 1/2, every problem ( qp ) is NP-hard. Similarly, they showed that ( qp ) is strongly NP-hard [Chen et al., 2014, Theorem 5] by a reduction from the 3-partition problem which is known to be strongly NP-hard [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. The invertible transform used in [START_REF] Chen | Complexity of unconstrained l2-lp minimization[END_REF] is defined by:

x = (2λ) 1/p x, H = (2λ) -1/p H.

(5.1)

Unfortunately, (5.1) is not well-defined when p = 0. Therefore, [Chen et al., 2014, Theorems 3 and 5] do not apply to ( 0 P ) when λ = 1/2. Using a different approach, Huo and Chen's paper [START_REF] Huo | Complexity of penalized likelihood estimation[END_REF] addresses the penalized least-squares problem defined by: min for every 0 ≤ τ < τ 0 .

C3. For the aforementioned τ 0 , if τ 1 , τ 2 < τ 0 then φ(τ 1 ) + φ(τ 2 ) ≥ φ(τ 1 + τ 2 ).

C4. For every 0 ≤ τ < τ 0 , φ(τ ) + φ(τ 0τ ) > φ(τ 0 ).

(5.2)

The proof of [Huo and Chen, 2010, Theorem 3.1] is by a reduction from the NP-complete problem X3C to the decision version of (PLS); this leads to the NP-completeness of the decision version of (PLS) and so the NP-hardness of (PLS) [START_REF] Huo | Complexity of penalized likelihood estimation[END_REF], Appendix 1]. The authors claimed that the 0 penalty function satisfies conditions C1-C4 for τ 0 = d = 1. Therefore, the (PLS) problem with the 0 penalty function is NP-hard [Huo and Chen, 2010, Corollary 3.2]. Unfortunately, it turns out that the 0 penalty does not fulfill condition C4 as claimed. Indeed, for τ = 0 the strict inequality (5.2) becomes φ(0) > 0. Besides, in the proof [START_REF] Huo | Complexity of penalized likelihood estimation[END_REF], Appendix 1], the inputs of the decision problem are not guaranteed to have rational values. This might also violate the polynomiality of the reduction. Therefore, [Huo and Chen, 2010, Theorem 3.1] does not apply to ( 0 P ).

In [START_REF] Huo | Complexity of penalized likelihood estimation[END_REF], the authors also mention an alternate proof of NP-hardness of ( 0 P ) from Huo and Ni's earlier paper [START_REF] Huo | When do stepwise algorithms meet subset selection criteria?[END_REF] as a special case of their results. In this proof [Huo and Ni, 2007, Appendix A.1], the relation between ( 0 P ) and ( 0 C ) is established using the principle of Lagrange multiplier. More precisely, the authors introduce an instance of ( 0 C ) in which is defined from the minimizer of ( 0 P ) and argue that solving ( 0 P ) is equivalent to solving the mentioned instance of ( 0 C ), which is known to be NP-hard [Natarajan, 1995]. There are a number of issues in the NP-hardness proof in [START_REF] Huo | When do stepwise algorithms meet subset selection criteria?[END_REF]. For instance, the proposed transformation between ( 0 P ) and ( 0 C ) is not a polynomial time reduction. Besides, it is well known that ( 0 P ) and ( 0 C ) are not equivalent [Nikolova, 2013].

New analysis on penalized 0 minimization

To prove that a problem T is NP-hard, one must establish a polynomial time reduction (briefly called reduction hereafter) from some known NP-hard or NP-complete problem to T [Leeuwen, 1990]. Roughly speaking, the reduction from a problem T1 to another problem T2 implies that T1 is not harder than T2. Therefore, if there exists a reduction from T1 to T2 and if T1 is NP-hard, T2 must be NP-hard too. The NP-hardness proofs in [Natarajan, 1995] and [START_REF] Davis | Adaptive greedy approximations[END_REF]] use this principle. As an adaptation of Natarajan's construction, we prove the NP-hardness of ( 0 P ) using the same principle as follows.

Theorem 5.1. Problem ( 0 P ) is NP-hard.

The proof is by a reduction from the known NP-complete problem X3C to ( 0 P ). The proof contains three steps: (1) Construct an instance of ( 0 P ) from a given instance of X3C; (2) Construct a solution of ( 0 P ) from a solution of X3C; (3) Construct a solution of X3C from a solution of ( 0 P ).

We only need to consider the case where 0 < λ < 3. For the case where λ ≥ 3, we can always scale3 the problem to the case where λ < 3 by respectively replacing y, H and λ by 1 √ t y, 1 √ t H and λ t where t > 1 is a positive integer satisfying 3 ≤ λ < 3t.

Construction of an instance of ( 0 P ) from a given instance of X3C

Given an instance of X3C: S = {s 1 , s 2 , ..., s m } is a set of m elements. C is a collection of n triplets c j , 1 ≤ j ≤ n. Without loss of generality we can assume that m is a multiple of 3 since otherwise there is trivially no exact cover so no solution of X3C. We now construct an instance of ( 0 P ). Let y = [1, 1, ..., 1] T ∈ R m . Let H = (h ij ) 1≤i≤m,1≤j≤n where h ij = 1 if s i ∈ c j and h ij = 0 otherwise. Let λ ∈ Q, 0 < λ < 3. Let F (x) := y -Hx 2 2 + λ x 0 .

(5.3)

Construction of a solution of ( 0 P ) from a solution of X3C

Assume that there is a subcollection of disjoint triplets Ĉ which exactly covers S. Let x * = [x * 1 , x * 2 , ..., x * n ] T where x * j = 1 if c j ∈ Ĉ and x * j = 0 otherwise. We will prove that x * is a solution of ( 0 P ).

Since Ĉ exactly covers S, | Ĉ| = m/3 and y = Hx * . Thus, x * 0 = m/3 and

F (x * ) = 0 + λ m 3 = λ m 3 .
Suppose that there exists x such that F (x) < F (x * ) = λ m 3 .

(5.4)

Let us show that this leads to a contradiction.

Since F (x) ≥ λ x 0 , from (5.4) we have x 0 < m/3. Therefore, we can rewrite x 0 = m/3q for some q ∈ N, 1 ≤ q < m/3. Note that H x has m entries. Since the number of non-zero entries of H x identifies with the number of elements s i recovered by the subcollection corresponding to x, this number cannot exceed 3 x 0 = m -3q. As a result, the number of zero entries of H x must be between 3q and m. Since y is the all-one vector, y -H x has at least 3q entries valued 1, which implies y -H x 2 2 ≥ 3q.

(5.5) Hence, F (x) ≥ 3q + λ m 3 q = λ m 3 + (3λ)q > λ m 3 , (5.6) which contradicts (5.4). Therefore, x * is a solution of ( 0 P ).

Construction of a solution of X3C from a solution of ( 0 P )

Assume that x * is a solution of ( 0 P ). We will consider four cases as follows.

Case x * 0 > m/3 We deduce that X3C has no solution. Indeed, assume that Ĉ is an exact cover for S. Define x = [x 1 , x 2 , ..., x n ] T where x j = 1 if c j ∈ Ĉ and x i = 0 otherwise. Then we have

F (x) = λ m 3 < λ x * 0 ≤ F (x * )
which contradicts the fact that x * is a solution of ( 0 P ).

Case x * 0 < m/3 We deduce that X3C has no solution. Indeed, assume that Ĉ is an exact cover for S. Let x = [x 1 , x 2 , ..., x n ] T where x j = 1 if c j ∈ Ĉ and x i = 0 otherwise. Then we have F (x) = λ m 3 . Since x * 0 < m/3, we can write x * 0 = m/3q for some q ∈ N and 1 ≤ q < m/3. Similar to (5.6), we have F (x * ) > λ m 3 . Since F (x) = λ m 3 , we obtain F (x * ) > F (x) which contradicts the fact that x * is a solution of ( 0 P ).

Case where x * 0 = m/3 and y = Hx * We deduce that X3C has no solution. Indeed, assume that Ĉ is an exact cover for S. Define x = [x 1 , x 2 , ..., x n ] T where x j = 1 if c j ∈ Ĉ and x i = 0 otherwise. Then we have

F (x) = λ m 3 < y -Hx * 2 2 + λ x * 0 = F (x * )
which contradicts the fact that x * is a solution of ( 0 P ).

Case where x * 0 = m/3 and y = Hx * Let Ĉ be the collection of triplets c j such that the j th entry of x * is non-zero. Obviously, Ĉ is an exact cover for S so a solution of X3C. Thus Theorem 5.1 is proved.

It is notable that the proof above is also valid when F (x) := y -Hx p p + λ x 0 for any p ≥ 1. Indeed, one only needs to check whether (5.5) still holds when the 2 norm is replaced by the p norm with p ≥ 1. This is the case since y -H x has at least 3q entries equal to 1. Therefore, we have the following generalization of Theorem 5.1.

Theorem 5.2. Problem min x y -Hx p p + λ x 0 is NP-hard for p ≥ 1.

Hardness of non-negative 0 minimization problems

The NP-hardness of non-negative 0 minimization problems is a consequence of NP-hard proofs of ( 0 C ) [Natarajan, 1995], ( 0 C K ) [START_REF] Davis | Adaptive greedy approximations[END_REF]] and ( 0 P ) (Theorem 5.1). Indeed, all these proofs consist in a reduction from X3C and the solution that established equivalence is binary. Therefore, the additional non-negativity constraints do not change the validity of these proofs. In other words, one can repeat the same proofs as in [Natarajan, 1995,Davis et al., 1997] and that of Theorem 5.1 for the corresponding non-negative 0 minimization problems ( 0 C +), ( 0 C K +) and ( 0 P +). Another way to prove the NP-hardness of non-negative 0 minimization problems is by a reduction from the corresponding 0 minimization problems which are known to be NP-hard. In this reduction, the instance of non-negative problems is defined by

ỹ = y, H = [H, -H], x = x + x -
where x + = max{x, 0}, x -= max{-x, 0}. Naturally, by this construction, one gets x ≥ 0, x 0 = x 0 and H x = Hx. The proofs (skipped for brevity) contain three steps similar to that of Theorem 5.1.

Therefore, we can state the following theorem without proof.

Theorem 5.3. ( 0 C +), ( 0 C K +) and ( 0 P +) are NP-hard.

In the same spirit and using the same argument as at the end of Section 5.4 one can directly extend Theorem 5.2 to the non-negative setting.

Theorem 5.4. Problem min x≥0 y -Hx p p + λ x 0 is NP-hard for p ≥ 1.

Conclusion

NP-hardness of penalized 0 minimization problems cannot be deduced from previous complexity analyses, as stated in [Chen et al., 2014,Huo andChen, 2010]. Here, we introduced a new proof of NP-hardness of penalized 0 minimization problems by an adaptation of Natarajan's construction [Natarajan, 1995] . Besides, we showed that the 0 minimization problems with non-negative constraints are also NP-hard. Therefore it is of interest to propose effective heuristic algorithms such as NNOG algorithms (see Chapter 2). This work can be extended in several directions. For instance, researchers interested in what makes NP-hard problems even harder might be interested in the strong NP-hardness of the aforementioned optimization problems. As it is widely believed that X3C is strongly NPcomplete, one might easily deduce the strong NP-hardness of ( 0 C ), ( 0 C K ) and other problems which are reduced from X3C. However, to the best of our knowledge, X3C is only proved to be NP-complete [Garey and Johnson, 1979, pp. 53, 221] and the strong NP-completeness has not been rigorously shown yet. Therefore, we believe that the question of strong NP-hardness of (non-negative) 0 minimization problems is not trivial and needs more work in future.

Besides, as (non-negative) 0 minimization problems are NP-hard, it would be interesting to know if the associated decision problems are in NP (so being NP-complete). Let us consider the decision problem associated with ( 0 C ): given y ∈ Q m , H ∈ Q m×n , a positive rational number and a positive integer K, does there exist x ∈ R n such that y -Hx 2 ≤ and x 0 ≤ K? This decision problem should be in NP since if one can guess a rational solution x, it can be verified in polynomial time if y -Hx 2 ≤ and x 0 ≤ K. Similarly, we conjecture that the decision version of other optimization problems mentioned in the chapter are in NP as well.

Another perspective is the approximability of aforementioned NP-hard problems. The hardness of approximating ( 0 C ) was discussed in [START_REF] Amaldi | [END_REF]Kann, 1998, Tillmann, 2015]. It was shown that approximating ( 0 C ) to within a factor of (1α) ln(n), 0 < α < 1 is NP-hard [Tillmann, 2015]. Examining whether similar results can be obtained on other NP-hard problems presented in the chapter would require more involved theoretical analysis, which is left for future work.

Chapter 6 The thesis addresses non-negative sparse approximation by greedy approach. The nonnegative extension of orthogonal greedy algorithms is investigated to target the 0 minimization problem under non-negativity constraints. The thesis contributions cover both algorithmic aspect and theoretical guarantee of non-negative orthogonal greedy (NNOG) algorithms.

Conclusion and perspectives

Summary of the contributions

In Chapter 2 we presented a unified framework of NNOG algorithms. We introduced the concepts of descending atom and descent selection rule and integrated the support compression step into the structure of standard NNOG algorithms. The resulted NNOG algorithms have simple structure but more advantageous structural properties in comparison with existing approximate schemes. More importantly, the proposed NNOG algorithms can naturally lend themselves to a recursive (so fast) implementation that existing approximate schemes can hardly do without some loss in the algorithmic behavior. The proposed recursive implementation uses active-set method with warm start initialization to solve the non-negative least-square (NNLS) subproblems. We also elaborated on this implementation and suggested further accelerations using several pruning strategies and vectorized computation.

The proposed implementation is validated through a rich set of simulations in Chapter 3. The proposed NNOG implementations are compared with the corresponding unconstrained versions and with existing approximate schemes, in terms of computing time and reconstruction accuracy. According to our comparison, the proposed implementation significantly reduces the cost of NNOG algorithms and appears to be more advantageous than existing approximate schemes. We also compared NNOG algorithms with other sparse solvers in their non-negative setting such as NNCoSaMP, NNSP, NNHTP and NLARS. It turns out that NNOG algorithms are quite competitive with other non-negative sparse solver in terms of reconstruction accuracy. However, NNOG algorithms have the advantage of a fast implementation that the competitors hardly apply.

While Chapters 2-3 are dedicated to algorithmic development of NNOG algorithms, Chapter 4 discusses the theoretical guarantee issue. We presented for the first time a unified K-step exact support recovery analysis of NNOG algorithms when the mutual coherence of the dictionary is lower than 1/(2K -1). Note that the mutual coherence value considered in exact recovery analysis is different from that in the simulations in Chapter 3, so the conclusion of Chapter 4 differs from the empirical observation in Chapter 3. More precisely, exact recovery analysis is derived from easy problem with low mutual coherence so that all NNOG algorithms perform the same. On the contrary, the mutual coherence value considered in simulations is close to 1. This yields a highly ill-posed problem for which the performances of NNOG algorithms differ. It is usually difficult to carry out exact recovery analyses when the mutual coherence is high but, even when the mutual coherence is low, the exact recovery analysis in the non-negative setting was not done before. Therefore, our analysis is strongly novel and it can be investigated in the context of higher mutual coherence as a perspective.

Finally we showed in Chapter 5 that non-negative 0 minimization problem is actually NPhard. Therefore, it is of interest to develop efficient heuristic algorithms. The study of NNOG algorithms and their fast implementation in Chapters 2-4 is hence meaningful.

Perspectives

Algorithmic aspect

The proposed framework of NNOG algorithms 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, forward-backward 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. Besides, the combination between NNOG algorithms and 1 solvers might generate a hybrid algorithm with higher performance than both greedy algorithms and 1 solvers as suggested in [START_REF] Wen | A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization, and continuation[END_REF].

Theoretical guarantee

Our contribution related to K-step exact support recovery analysis of NNOG algorithms when the mutual coherence of the dictionary is lower than 1/(2K -1) can have different extensions. For instance, the mutual coherence bound 1/(2K -1) might be weakened by taking into account prior information on the sought signal such as the decay of the magnitudes of the nonzero coefficients [START_REF] Herzet | Relaxed recovery conditions for OMP/OLS by exploiting both coherence and decay[END_REF]. With the same idea of deriving exact recovery analysis under weaker conditions, Tropp's exact recovery condition (ERC) [Tropp, 2004] is worth being considered, especially when this condition has a strong connection with mutual coherence value. Finally it is also interesting to check if the usual analysis of OMP and OLS based on RIP condition [START_REF] Davenport | Analysis of orthogonal matching pursuit using the restricted isometry property[END_REF] can be extended to the non-negative setting in a similar way.

Figure 1

 1 Figure 1.1: A simple illustration of sparse approximation. The red color indicates the non-zero coefficients and the corresponding dictionary atoms.

Figure 1 Figure 1 . 4 :

 114 Figure 1.3: Four quantized images in a Tomo-PIV simulation provided by ONERA.

Figure 1 . 5 :

 15 Figure 1.5: Spectral unmixing (image credit: Boeing)

  Figure 1.6: Orthogonal projection of y onto span (H S ) and the resulting residual r.

3T 4 S

 4 ← set of indices of the 2K largest positive entries of H t r ; ← S ∪ T ; 5 x(S) ← arg min z≥0 y -H S z 2 ; 6 S ← set of indices of the K largest magnitude of x ; 7 Set x i = 0 for all i / ∈ S ; 8 r ← y -H S x(S) ; 9 end

  Nguyen, Jerome Idier, Charles Soussen, El-Hadi Djermoune. Non-Negative Orthogonal Greedy Algorithms. IEEE Transactions on Signal Processing, vol. 67(21), pp. 5643-5658. 2019.

Definition 2. 1 .

 1 Let us call a positive support related to the full-size NNLS problem min x≥0 y -Hx 2 ,

11

  Update xV (call Algorithm 15); 12 end 13 x ← xV ; 14 r V ← y -Hx ; 15 end Algorithm 14: NNOG with active-set implementation. input : y, H, K output: x := x+ S (with S a positive support) 1 x ← 0 ; S ← ∅ ; r S ← y ; 2 while |S| < K and max i∈ S h t i r S > 0 do 3 Select an index ∈ S by a selection rule S(y, H, S) ; 4 Call [C, x] = AS_NNLS(y, H, S ∪ { }, S, x); 5 S ← C ; 6 r S = y -Hx ; 7 end

Figure 2

 2 Figure 2.2: 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.

  , 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

Figure 3

 3 Figure 3.1: A simple illustration of sparse deconvolution with Gaussian kernel. (top) Discrete data are represented with bullets. Sampling step is set to 1.The three spikes (one in red, one in blue, one in green) corresponding to the centers of the constitute Gaussians h 1 , h 2 , h 3 are located on the sampling grid. (bottom) The data y can be represented as a linear combination of three Gaussians h 1 , h 2 , h 3 with corresponding weights x 1 , x 2 , x 3 .

Figure 3

 3 Figure 3.2: An example of simulated data and the corresponding ground truth in a sparse deconvolution problem with Gaussian kernel.

Figure 3

 3 Figure 3.3: A simple illustration of the dictionary in super resolution problem. H = [h 1 , h 2 , h 3 , h 4 ]. Discrete atoms are represented by bullets. Sampling grid has step size 1 while the grid for sparse reconstruction has step ∆ = 0.5. Each atom is a delayed version of the previous atom by a step ∆.

Figure 3 . 4 :Figure 3

 343 Figure 3.4: Principle of wood classification by NIR spectra decomposition.

Figure 3

 3 Figure3.6: A simple illustration of sparse decomposition of one spectrum. The Gaussian dictionary consists of three Toeplitz blocks (red, green and blue) corresponding to the convolution matrices related to the Gaussian impulse response of given widths (smallest width in red and largest width in blue). The spectrum y can be decomposed into three spectra of different widths (one red, one green and one blue) with the corresponding weights in vector x. The three spectral components are highlighted in bold and the non-zero weights are in color (red, green and blue). Figure was kindly shared by Dr. Hassan Mortada.

Figure

  Figure3.6: A simple illustration of sparse decomposition of one spectrum. The Gaussian dictionary consists of three Toeplitz blocks (red, green and blue) corresponding to the convolution matrices related to the Gaussian impulse response of given widths (smallest width in red and largest width in blue). The spectrum y can be decomposed into three spectra of different widths (one red, one green and one blue) with the corresponding weights in vector x. The three spectral components are highlighted in bold and the non-zero weights are in color (red, green and blue). Figure was kindly shared by Dr. Hassan Mortada.

  Figure 3.8: Complexity analysis of SNNOLS and NNOLS for a simulated 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 Lt to L). ρ ↓ : ratio of descending atoms; ρ + cand : ratio of atoms resulting in a positive support in NNOLS; ρ II : ratio of atoms discarded by type-II pruning in NNOLS; ρ + sel : rate of positive supports found by SNNOLS. (b) Histogram of the average number of iterations L required to reach a support of cardinality K = 80.

Figure 3 . 9 :

 39 Figure 3.9: Comparison of NNOG algorithms. Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, P n = 10 -2 , and a dictionary of size 2400 × 2388.

  Figure3.10: One trial of data y and ground truth x * when K = 10 and the reconstruction using NNOMP. NNOMP detects correctly 9 spikes (in red) and yields one false detection (in blue).

Figure 3

 3 Figure 3.11: Comparison between NNOMP, NNOLS and SNNOLS in a super resolution problem where resolution grid step is set to ∆ = 2 -N for N varying between 1 and 7.

Figure 3

 3 Figure3.12: Comparison of NNOG algorithms with unconstrained greedy algorithms. Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, P n = 10 -2 , and a dictionary of size 2400 × 2388.

Figure 3

 3 Figure3.13: Data approximation and sparse recovery of an NIR spectrum using NNOG algorithms and unconstrained greedy algorithms. The Gaussian dictionary contains 2998 atoms and the sparsity level is set to K = 20.

Figure 3

 3 Figure3.14: Comparison of NNOG algorithms with with approximative fast non-negative greedy algorithms. Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, P n = 10 -2 , and a dictionary of size 2400 × 2388.

Figure 3

 3 Figure3.15: Data approximation and sparse recovery of an NIR spectrum using NNOG algorithms and approximative fast versions. The Gaussian dictionary contains 2998 atoms and the sparsity level is set to K = 20.

Figure 3

 3 Figure 3.16: Comparison of NNOG algorithms with non-negative versions of CoSaMP, SP, HTP.Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, P n = 10 -2 , and a dictionary of size 2400 × 2388.

  Figure3.17: Data approximation and sparse recovery of an NIR spectrum using non-negative greedy algorithms. The Gaussian dictionary contains 2998 atoms and the sparsity level is set to K = 20.

Figure 3

 3 Figure 3.18: Comparison of NNOG algorithms with NLARS. Performances are averaged over 400 trials in a simulated sparse deconvolution problem with Gaussian kernel with σ = 2, P n = 10 -2 , and a dictionary of size 2400 × 2388.

Figure 3

 3 Figure3.19: Data approximation and sparse recovery of an NIR spectrum using NLARS and NNOG algorithms. The Gaussian dictionary contains 2998 atoms and the sparsity level is set to K = 20.

Figure 4 . 1 :

 41 Figure 4.1: The positive span assumption does not propagate to the residual.

  |x| and |z| ∈ R n denote the non-negative vectors whose entries are equal to |x i | and |z i |, respectively. Obviously, H |x| = H |z|, and µ(H ) = µ(H) = 1 2K-1 .

4. 4

 4 Numerical study 4.4.1 Comparison of Oxx and their sign-aware versions

Figure 4 . 2 :

 42 Figure 4.2: Rate of exact support recovery (S K = S * ) and non-negativity at all iterations (∀k, x(k) ≥ 0) w.r.t. the signal-to-noise ratio for a simulated data experiment involving Gaussianshaped atoms. SNNOLS (not shown here) yields the exact same curve as NNOLS.

Figure 4 . 3 :

 43 Figure 4.3: Behavior of Oxx and 1 -homotopy in the case of a toy problem (m = n = 5) corresponding to a 4-sparse noiseless representation. (a) The OMP and OLS iterates are identical and yield magnitudes x(k) i with non-monotonous variations throughout the iterations k. (b) On the contrary, the magnitudes of the 1 -homotopy iterates are always increasing. The ground truth magnitudes x * i are represented with bullets for k = 4, and are exactly recovered after K = 4 iterations (x (k) = x * ).

  .34) hS k-1 = 0 in (4.33) follows from Lemma 4.9 and the full rankness assumption of H S k .(4.34) holds because hS k-1h lays in span(H S k-1 ), and the OMP residual r k-1 is orthogonal to

  .49) (4.49) and (4.48) with j ← p indeed imply that ct r k-1 > 0.

  i |), (PLS)where φ is a penalty function mapping non-negative values to non-negative values. The authors showed that (PLS) is NP-hard if the penalty function φ satisfies the following four conditions[Huo and Chen, 2010, Theorem 3.1]: C1. φ(0) = 0 and ∀ 0 ≤ τ 1 < τ 2 , φ(τ 1 ) ≤ φ(τ 2 ).C2. There exists τ 0 > 0 and a constant d > 0 such that φ(τ ) ≥ φ(τ 0 )d(τ 0τ ) 2

Contents 6 . 1

 61 Summary of the contributions . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2.1 Algorithmic aspect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2.2 Theoretical guarantee . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.2.2 Sparse deconvolution for super resolution . . . . . . . . . . . . . . . . . . 3.2.3 Decomposition of NIR spectra into elementary Gaussian features . . . . . 3.3 Validation of NNOG accelerations . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.3.2 Computation burden of SNNOLS and NNOLS . . . . . . . . . . . . . . . 3.3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Comparison of three NNOG algorithms . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.4.2 Super resolution problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Comparison with unconstrained greedy algorithms . . . . . . . . . . . . . . . . . 3.5.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.5.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Comparison with fast approximate non-negative greedy algorithms . . . . . . . . Comparison with non-negative extension of LARS . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents 3.1 Notations and abbreviations 3 Numerical validation of NNOG algorithms
	Notations and abbreviations Notations	ix
	y		Data vector
	Résumé de thèse H Dictionary (observation matrix)	xi
	x		Solution vector
			Error vector
	S		Support of solution
	|S| S	Cardinality of S Complement of S
	. t	• 0 • p •	0 -"norm" counting the number of nonzero entries p -norm Euclidean norm Transpose operator
	. †		Moore-Penrose pseudo-inverse
	. ⊥	Orthogonal complement
	span(H) Column space of H
	H S	Subdictionary indexed by S
	x(S)	Subvector indexed by S
	x S	Subvector indexed by S (in Chapter 4 only)
	xS	ULS solution related to S
	x+

3.6.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.6.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Comparison with non-negative extensions of CoSaMP, SP and HTP . . . . . . . . 3.7.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.7.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 3.8.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 3.8.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 S NNLS solution related to S h i Column of dictionary H (atom) hS i Orthogonal projection of h i onto span(H S ) ⊥ (projected atom) gS

  Structure générale des algorithmes gloutons orthogonaux pour minimiser (1). h i désigne la i-ème colonne du dictionnaire H, et x+ S représente la solution NNLS associée au support S.

	xii	Contents
	Algorithm 1:
		analyses classiques

entrées: y, H, K sorties : x

1 x ← 0 ; S ← ∅ ; r S ← y ; 2 tant que |S| < K et max i / ∈S h t i r S >

0 faire 3 Sélectionner un atome / ∈ S ; 4 S ← S ∪ { } ; 5 x ← x+ S ; 6 S ← supp(x) ; 7 r S = y -Hx ;

  Non-negative sparse reconstruction . . . . . . . . . . . . . . . . . . . . . 2 1.2 Some applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Tomographic PIV reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Spectral unmixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 1
	Introduction
	Contents
	1.1

6 1.3 State-of-the-art algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.1 Greedy algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.2 Non-negative greedy algorithms . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3.3 Convex solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.4 Non-convex solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 State-of-the-art exact recovery analyses of greedy algorithms . . . . . . 15 1.4.1 Exact recovery with greedy algorithms . . . . . . . . . . . . . . . . . . . . . 15 1.4.2 Exact recovery with non-negative greedy algorithms . . . . . . . . . . . . . 15 1.5 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.5.1 Structural properties of non-negative greedy algorithms . . . . . . . . . . . 16 1.5.2 Fast implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.5.3 Theoretical guarantees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.5.4 NP-hardness of non-negative 0 minimization . . . . . . . . . . . . . . . . . 17 1.6 List of publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

Table 1

 1

.1: Typical greedy algorithms and the corresponding problems addressed 1.3.1.1 One-by-one updating strategy

  Thanh T. Nguyen, Charles Soussen, Jerome Idier, El-Hadi Djermoune. Algorithmes gloutons orthognaux pour la reconstruction de signaux parcimonieux positifs. GRETSI, Lille, France. August 2019. hal-02149677. -Thanh T. Nguyen, Charles Soussen, Jerome Idier, El-Hadi Djermoune. An optimized version of non-negative OMP. GRETSI, Juan-les-pins, France. September 2017. hal-01585253. Matlab implementation of Non-Negative Orthogonal Greedy Algorithms (CodeOcean), supplementary material to the paper: Thanh T. Nguyen, Jerome Idier, Charles Soussen, El-Hadi Djermoune. Non-Negative Orthogonal Greedy Algorithms. IEEE Transactions on Signal Processing, vol. 67(21), pp. 5643-5658. 2019.

	Open source software

• International conference with proceedings -Thanh T. Nguyen, Charles Soussen, Jerome Idier, El-Hadi Djermoune. NP-hardness of 0 minimization problems: revision and extension to the non-negative setting. International Conference on Sampling Theory and Applications, Bordeaux, France. July 2019. hal-02112180. • International workshop -Thanh T. Nguyen, Charles Soussen, Jerome Idier, El-Hadi Djermoune. Exact recovery analysis of non-negative Orthogonal Matching Pursuit. Signal Processing with Adaptive Sparse Structured Representations (SPARS) workshop, Toulouse, France, July 2019.

• National conference with proceedings 1. Introduction -

1

  Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Non-negative greedy algorithms . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.1 Basic definitions and notations . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2.2 Non-negative orthogonal greedy algorithms . . . . . . . . . . . . . . . . . . 23 2.2.3 Descending atoms and descent selection rules . . . . . . . . . . . . . . . . . 24 2.2.4 Examples of descent selection rules . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.5 On the usefulness of support compression . . . . . . . . . . . . . . . . . . . 27 2.3 Active-set NNLS algorithms and recursivity . . . . . . . . . . . . . . . . 28 2.3.1 Fast active-set algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.2 Warm start and full recursivity . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.3 Step-by-step illustrative example of NNOG . . . . . . . . . . . . . . . . . . 30 2.3.4 Connection between AS-NNLS and NNOMP . . . . . . . . . . . . . . . . . 31 2.4 Acceleration of NNOG algorithms . . . . . . . . . . . . . . . . . . . . . . 31 2.4.1 Atom selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.2 Coefficient update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4.3 Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.A Proof of technical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.A.1 Proof of Proposition 2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.A.2 Proof of Proposition 2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.B Recursive implementation of ULS . . . . . . . . . . . . . . . . . . . . . . 34

  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+ 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. Lemma 2.3. For any S, let C = supp(x + S ) (where neither x+ S nor C are necessarily unique if H S is not full column rank). Then we have

	x+ S = xC = x+ C .	(2.6)
	Proof. Using Lemma 2.1, we have H t C r +

S

is orthogonal to span(H S ). S = 0. Since r + S = y -H x+ S , we get H t C y = H t C H x+ S . Thus, x+ S is a ULS solution (denoted by xC ) associated to support C. We have also xC = x+ C since xC ≥ 0.

  ⊥ . Finally, the positivity constraint on v in (2.19)-(2.20) turns out to be inactive. Indeed, according to Proposition 2.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 (2.1) that the minimum error norm in (2.19)-(2.20) is then strictly lower than r S . Hence, the optimal variable v in (2.19)-(2.20) 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 . Figure 2.2.4 illustrates the three rules in a simple case where all results differ. 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 (2.16)-(2.17) as those of OMP and OLS, respectively.

	is nonzero.
	Finally, (2.15) with (2.19) and (2.15) with (2.20) rewrite as (2.15) with (2.16) and (2.15) with
	(2.17), respectively.

  Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 38 3.2.2 Sparse deconvolution for super resolution . . . . . . . . . . . . . . . . . . . 40 3.2.3 Decomposition of NIR spectra into elementary Gaussian features . . . . . . 41 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 54 3.6.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.7 Comparison with non-negative extensions of CoSaMP, SP and HTP . 54 3.7.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 56 3.7.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.7.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.8 Comparison with non-negative extension of LARS . . . . . . . . . . . . 58

3.3 Validation of NNOG accelerations . . . . . . . . . . . . . . . . . . . . . . 44 3.3.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 45 3.3.2 Computation burden of SNNOLS and NNOLS . . . . . . . . . . . . . . . . 46 3.3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Comparison of three NNOG algorithms . . . . . . . . . . . . . . . . . . . 48 3.4.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 48 3.4.2 Super resolution problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5 Comparison with unconstrained greedy algorithms . . . . . . . . . . . . 50 3.5.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 50 3.5.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.6 Comparison with fast approximate non-negative greedy algorithms . . 54 3.6.1 3.8.1 Sparse deconvolution with Gaussian kernel . . . . . . . . . . . . . . . . . . 59 3.8.2 Decomposition of NIR spectra . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.8.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Table 3 .

 3 1: Acceleration 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 NNOMP SNNOLS NNOLS NNOMP SNNOLS NNOLS Non-negativity loss Iterations
	20	2.7	10.5	792	1.2	1.0	1.5	20	22	21
	40	4.5	11.5	1110	1.2	1.2	1.9	41	51	43
	60	6.4	11.4	952	1.1	1.2	2.1	65	97	73
	80	8.5	11.5	646	1.1	1.3	4.1	95	152	121

Table 3 .

 3 2: CPU time and normalized approximation error of NNOG algorithms and unconstrained 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 **OLS OLS+ NNOLS
	Time (ms)	23	24	28	25	26	40
	y -Hx / y × 10 2	6.6	16.4	9.6	5.7	16.6	6.4

Table 3

 3

	.3: CPU time and normalized approximation error of different implementations of NNOG
	algorithms for sparse decomposition of NIR spectra. Average over 50 spectra.
	Algorithm	NNOMP FNNOMP SNNOLS NNOLS FNNOLS
	Time (ms)	28	25	30	40	2720
	y -Hx / y × 10 2	9.6	9.6	6.4	6.4	7.0

Table 3 .

 3

	Algorithm	NNOMP SNNOLS NNOLS NNCoSaMP NNSP NNHTP
	Time (ms)	28	30	40	45	57	40
	y -Hx / y × 10 2	9.6	6.4	6.4	10.1	8.2	9.7

4: CPU time and normalized approximation error of several non-negative greedy algorithms for sparse decomposition of NIR spectra. Average over 50 spectra.

Table 3 .

 3 5: CPU time and normalized approximation error of NLARS and NNOG algorithms for sparse decomposition of NIR spectra. Average over 50 spectra.

		Algorithm		NNOMP SNNOLS NNOLS NLARS
		Time (ms)			28	30	40	77
			y -Hx / y × 10 2		9.6	6.4	6.4	21.1
	0.06
	0.05
	0.04
	0.03
	0.02
	0.01
	0
	4000	5000	6000	7000	8000	9000	10000

  .27) where p H and p H ⊥ are the orthogonal projections of r onto span(H) and span(H) ⊥ , respectively, and r H 2 ∈ span(H 2 ). Note though that the decomposition p H = (H † r) 1 h 1 + r H 2 is not orthogonal. Rewriting h 1 as P S 2 h 1 + P ⊥ S 2 h 1 yields the orthogonal decomposition:r = H † r 1 P ⊥ S 2 h 1 + (r H 2 + P S 2 h 1 ) + p H ⊥ . (4.28) (4.26) is obtained directly from (4.28) by writing the inner product r, P ⊥ S 2 h 1 because the latter decomposition is orthogonal. The denominator in (4.26) is nonzero because of the full rankness of H.

  .38) Moreover, since the OLS iterate at iteration k is non-negative, it is also an NNLS solution related to S k . Therefore, min

	(4.39)
	From (4.37)-(4.39), we get
	∀i / ∈ S k , min z≥0

z y -H S k z 2 = min z≥0 y -H S k z 2 .

  .75)In order to exploit the Oxx selection rule which is based on the atoms ci defined in (4.4), let us rewrite (4.75) with respect to matrix C. Obviously, one can rewrite HS k

1

  Contents 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.2 Background on NP-hardness and constrained 0 minimization problems 84 5.3 Existing analyses on penalized 0 minimization . . . . . . . . . . . . . . 84 5.4 New analysis on penalized 0 minimization . . . . . . . . . . . . . . . . . 86 5.5 Hardness of non-negative 0 minimization problems . . . . . . . . . . . . 88 5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Exact support recovery with NNOG algorithms 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Funded by the French National Research Agency (No. ANR-15-CE23-0021). http://becose.univ-lorraine.fr

This work is done within the scope of the collaboration with ONERA in the BECOSE project.

In the literature, OLS is also known as Order Recursive Matching Pursuit (ORMP)[START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF], Optimized Orthogonal Matching Pursuit (OOMP)[START_REF] Lowe ; Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF], and Pure Orthogonal Matching Pursuit[Foucart, 2013].

in the sense that the current iterate is used to accelerate the computation of the next iterate.

r ← y -H S x(S) ;

This chapter is an adaptation of the paper[Nguyen et al.,

2019a]

However, this kind of operation was introduced in[START_REF] Leichner | A strictly improving linear programming Phase I algorithm[END_REF], Morup et al., 2008].

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

kindly provided by Dr. Mehrdad Yaghoobi.

This chapter is an adaptation of the paper[Nguyen et al.,

2019d]

referred to as Fuchs corollary condition in[Plumbley, 2007] as opposed to the sharp (but more involved) Fuchs condition in[Plumbley, 2007, Theorem 2]. The latter condition solely depends on the sign pattern as well.

This chapter is an adaptation of the paper[Nguyen et al.,

2019c]

The latter problem, denoted by X3C in[Natarajan, 1995, Garey andJohnson, 1979], is stated as follows: Given a set S and a collection C of

3-element subsets of S (called triplets), is there a subcollection of disjoint triplets that exactly covers S?

This transformation was kindly suggested by Andreas M. Tillmann at the SampTA conference in July 2019.

Acknowledgments

This thesis was impossible without the guidance of my supervisors Charles Soussen, Jérôme Idier and El-Hadi Djermoune. I am deeply grateful for their kindness and support during the three years of my PhD. I would like to thank the members of my jury: Charles Dossal, Loïc Denis, Christian Heinrich and Laura Rebollo-Neira for their time reviewing my work. It is a great honor for me to present my work to them. Special thanks go to Mehrdad Yaghoobi for providing his Matlab implementation of FNNOLS and Hassan Mortada for sharing the figure of spectra decomposition. I would like to acknowledge to my laboratory CRAN and the Inverse Problem Group (GPI) at L2S for providing the best working condition for my thesis preparation. Finally, I would like to thank my family for their unconditional love and my friends for their endless encouragement.

Chapter 4

Exact support recovery with NNOG algorithms 1

Contents

Résumé

De nombreux domaines applicatifs conduisent à résoudre des problèmes inverses où le signal ou l'image à reconstruire est à la fois parcimonieux et positif. Si la structure de certains algorithmes de reconstruction parcimonieuse s'adapte directement pour traiter les contraintes de positivité, il n'en va pas de même des algorithmes gloutons orthogonaux comme OMP et OLS. Leur extension positive pose des problèmes d'implémentation car les sous-problèmes de moindres carrés positifs à résoudre ne possèdent pas de solution explicite. Dans la littérature, les algorithmes gloutons positifs (NNOG, pour non-negative orthogonal greedy algorithms) sont souvent considérés comme lents, et les implémentations récemment proposées exploitent des schémas récursifs approchés pour compenser cette lenteur.

Dans ce manuscrit, les algorithmes NNOG sont vus comme des heuristiques pour résoudre le problème de minimisation 0 sous contrainte de positivité. La première contribution est de montrer que ce problème est NP-difficile. Deuxièmement, nous dressons un panorama unifié des algorithmes NNOG et proposons une implémentation exacte et rapide basée sur la méthode des contraintes actives avec démarrage à chaud pour résoudre les sous-problèmes de moindres carrés positifs. Cette implémentation réduit considérablement le coût des algorithmes NNOG et s'avère avantageuse par rapport aux schémas approximatifs existants. La troisième contribution consiste en une analyse de reconstruction exacte en K étapes du support d'une représentation K-parcimonieuse par les algorithmes NNOG lorsque la cohérence mutuelle du dictionnaire est inférieure à 1/(2K -1). C'est la première analyse de ce type.

Mots clefs : algorithmes gloutons orthogonaux, reconstruction parcimonieux, contrainte de positivité, moindres carrés positif, l'algorithme des contraintes actives, reconstruction de support exact, cohérence mutuelle, NP-difficile, problème minimisé 0 .

Abstract

Non-negative sparse approximation arises in many applications fields such as biomedical engineering, fluid mechanics, astrophysics, and remote sensing. Some classical sparse algorithms can be straightforwardly adapted to deal with non-negativity constraints. On the contrary, the non-negative extension of orthogonal greedy algorithms is a challenging issue since the unconstrained least square subproblems are replaced by non-negative least squares subproblems which do not have closed-form solutions. In the literature, non-negative orthogonal greedy (NNOG) algorithms are often considered to be slow. Moreover, some recent works exploit approximate schemes to derive efficient recursive implementations.

In this thesis, NNOG algorithms are introduced as heuristic solvers dedicated to 0 minimization under non-negativity constraints. It is first shown that the latter 0 minimization problem is NP-hard. The second contribution is to propose a unified framework on NNOG algorithms together with an exact and fast implementation, where the non-negative least-square subproblems are solved using the active-set algorithm with warm start initialization. The proposed implementation significantly reduces the cost of NNOG algorithms and appears to be more advantageous than existing approximate schemes. The third contribution consists of a unified K-step exact support recovery analysis of NNOG algorithms when the mutual coherence of the dictionary is lower than 1/(2K -1). This is the first analysis of this kind.