monsieur Stéphane Canu

Paul Honeine

Su Ruan

Je

Linlin Jia

Imad Rida

Ismaila Seck

Je tiens à remercier toutes les personnes qui ont contribué au succès de mes recherches et qui m'ont aidée lors de la rédaction de ma thèse. Ils sont

Résumé

Les méthodes parcimonieuses suscitent un vif intérêt pour l'estimation, de par leur capacité à proposer automatiquement un modèle interprétable ayant un très bon pouvoir prédictif. Ces méthodes se formulent comme des problèmes d'optimisation bi-critères où l'on cherche à obtenir un modèle minimisant simultanément un terme d'attache aux données et une pénalité visant à promouvoir la parcimonie via la sélection de caractéristiques pertinentes. Ce type de problème est aussi référencé sous le nom de codage parcimonieux.

Dans ce cadre, les données disponibles sont représentées par une matrice, considérée comme un dictionnaire, dont chaque colonne représente une caractéristique encore appelée atome. Le choix de ce dictionnaire joue un rôle très important et sa détermination est donc primordiale. Certaines études préconisent l'utilisation d'un dictionnaire prédéfini, par exemple à l'aide de polynômes ou d'ondelettes. Mais récemment, pour différentes applications, il s'est avéré plus efficace d'apprendre un dictionnaire dédié plutôt que d'utiliser un dictionnaire prédéfini. C'est le cas notamment en traitement du signal et des images, ainsi qu'en vision et en reconnaissance des formes.

Cette monographie traite de l'étude d'algorithmes d'apprentissage de dictionnaire, lorsque la parcimonie est favorisée par l'utilisation de la norme 0 , c'est-àdire par le contrôle explicite du nombre de caractéristiques à considérer. Elle est organisée en cinq chapitres. Le premier chapitre donne une brève introduction sur le contexte et les motivations du travail. Les contributions de la thèse et le plan du manuscrit sont également présentés.

Le deuxième chapitre porte sur la théorie de l'apprentissage de dictionnaire parcimonieux et notamment sur les problématiques d'optimisation bi-critères associées. Trois différents choix de pénalités pour le codage parcimonieux sont considérés : une mesure de comptage encore appelée pseudo-norme 0 (ou norme par abus de langage), la norme 1 et les normes p avec p > 1 ou 0 < p < 1. L'intérêt de la v norme 0 est qu'elle permet de compter et donc de contrôler explicitement le nombre de composantes d'un modèle. Ainsi, à performance de prédiction analogue, les modèles résultant de l'usage de la norme 0 sont plus parcimonieux que ceux obtenus en utilisant les autres normes 1 ou p (0 < p < 1 et 1 < p). Cependant, les problèmes d'optimisation associés à l'usage de la norme 0 sont combinatoires, non convexes, non différentiables et NP-difficiles. Pour ces raisons, il sont généralement considérés comme difficiles à résoudre et ne passant pas à l'échelle.

Le troisième chapitre donne un état de l'art sur les méthodes permettant de résoudre le problème de codage parcimonieux pour l'apprentissage de dictionnaire.

Il commence par présenter trois stratégies fréquemment considérées. La première est basée sur l'utilisation d'un algorithme glouton qui a donné les méthodes de poursuite et ses variantes (matching pursuit et orthogonal matching pursuit) pour obtenir une bonne approximation de la solution optimale globale du problème. La deuxième porte sur la relaxation de la norme 0 via l'utilisation de la norme 1 et les algorithmes associés comme le LASSO ou la « poursuite de base » (Basis Pursuit). La troisième utilise des méthodes du gradient, comme les méthodes de seuillage ou la méthode du gradient proximal.

Ce chapitre se poursuit par une présentation du problème d'apprentissage de dictionnaire qui consiste à estimer conjointement un dictionnaire pertinent et les coefficients pondérant ces atomes. Le problème d'optimisation associé est à la fois non convexe et NP-difficile. La manière classique d'aborder ce problème est d'utiliser une procédure de relaxation alternée en deux phases à la Gauss Seidel : 1) la phase de codage parcimonieux (sparse coding), qui consiste à estimer les coefficients de pondération en supposant le dictionnaire connu; 2) la phase de réactualisation de dictionnaire, qui consiste à estimer le dictionnaire en supposant cette fois les coefficients de pondération connus. La dernière partie de ce chapitre présente les principaux algorithmes existants permettant de résoudre le problème d'apprentissage de dictionnaire associé à la norme 0 : MOD (Method of Optimal Direction), K-SVD (K-Singular Value Decomposition) et SOUPDIL (Sum of OUter Products Dictionary Learning).

Le quatrième chapitre introduit une nouvelle méthode d'apprentissage de dictionnaire utilisant la norme 0 . Comme le problème d'optimisation associé à cette norme 0 est NP-difficile, les stratégies développées dans la littérature sont basées sur des approximations ne permettant d'obtenir que des solutions locales. Notre idée est de développer une approche alternative permettant d'obtenir la solution vi globale du problème en un temps raisonnable. Pour ce faire, nous proposons de reformuler le problème du codage parcimonieux comme un programme quadratique mixte en nombres entiers (MIQP) et d'utiliser un logiciel d'optimisation, comme par exemple Gurobi ou Cplex, pour obtenir cet optimum global. La principale difficulté de cette approche étant le temps de calcul, nous proposons deux méthodes permettant de la surmonter : l'ajout de contraintes complémentaires accélérant la convergence de l'algorithme et son initialisation par un minimum local obtenu grâce à une méthode du premier ordre de type gradient proximal. La méthode ainsi obtenue, nommée MIQP accéléré (AcMIQP) permet de diminuer significativement les temps de calculs de la procédure d'optimisation et d'augmenter la taille des problèmes pouvant être ainsi traités. Notre méthode d'apprentissage de dictionnaire AcMIQP a été appliquée sur un problème de débruitage d'images démontrant sa faisabilité et sa pertinence. Les résultats montent que notre méthode est plus performante que les méthodes de références, comme le gradient proximal et l'algorithme K-SVD, et atteint des résultats comparables à ceux de SOUPDIL.

Un autre intérêt de la formulation quadratique mixte en nombres entiers est sa flexibilité. Il est facile d'introduire dans le problème d'optimisation MIQP de nouvelles contraintes pourvues qu'elles soient linéaires. Tirant profit de cet avantage, le cinquième chapitre traite de l'apprentissage de dictionnaire à faible cohérence, c'est-à-dire à faible corrélation entre ses atomes. Si, nous l'avons vu, la norme 0 permet le contrôle du niveau de parcimonie, elle ne dit rien en revanche sur la cohérence de ce dictionnaire. Or, plusieurs études théoriques ont montré la pertinence, pour un dictionnaire, d'avoir une faible cohérence. Les méthodes existantes qui s'attaquent à ce problème reposent sur la relaxation des contraintes, par exemple en ajoutant une étape de décorrélation à chaque itération.

Nous proposons une alternative originale, basée sur les travaux du chapitre précédent, visant à résoudre ce problème d'apprentissage de dictionnaire à faible cohérence en intégrant explicitement des contraintes favorisant l'indépendance des atomes. L'estimation du dictionnaire sous contraintes d'incohérence est abordée en combinant la méthode du lagrangien augmenté (ADMM) et la méthode du gradient proximal alterné étendu (Extended Proximal Alternating Linearized Minimization, EPALM), adaptée à des familles de problèmes non convexes. L'efficacité de la méthode AcMIQP+EPALM est démontrée sur une application de reconstruction d'image. Today, people live in a world surrounded by diverse data. By improving smartphone technologies with cameras of higher and higher resolution or by alleviating user-generated content in the so-called data web such as Facebook and YouTube, large volume data of high dimension is produced every day. It makes a lot of demands on the storage capacity of the device, and moreover, to deal with these data, the high computation complexity is foreseen. However, it is noticed that in natural signals, namely audio, image, video, text and document, the number of information that makes sense is very small comparing to the whole signal, that is to say, the underlying model is sparse. The property of sparseness makes it possible to find a representation quite sparse replacing the original signal. As a result, the work on the sparse representation, rather than the signal, can greatly reduce computing complexity and decrease needs on storage capacity, which makes the research on sparse representation crucial with real benefits.

Sparse representation aims at representing a signal as linear combination, of which most coefficients are zeros, over a set of elementary signals where the set is called dictionary and the elementary signals are the atoms. The problem is often raised to find signals' regularity in signal processing or to extract images' features for tasks in pattern recognition [START_REF] Starck | Sparse image and signal processing: Wavelets and related geometric multiscale analysis[END_REF]. When the given dictionary is a basis, the problem involves solving a matrix equation, which is easy by introducing the matrix inverse. However, dictionaries with orthogonal atoms are not often able to take the most of sparse representations and signal reconstruction. Therefore, in most cases, the dictionary is learned as a redundant matrix, for which the resolution above is no more available. In fact, with a redundant dictionary, a solution to the problem may not be unique. Then, many problems will be raised: What is the condition for the uniqueness of the solution with sparsity constraint on the composition coefficients?

How to handle the problem of simultaneously keeping the signal reconstructed to a great extent and adapting it to certain task? Even today, the studies on these problems are still challenging and open.

Evolution in sparse representation

During the past decades, sparse representation has attracted numerous research attentions. Originally, it is inspired by the decomposition of a studied signal over a well-defined basis in signal processing. In the beginning, the Fourier transform projects the signal in a basis formed by sine and cosine functions of different frequencies [BB86]. Thus, the signal defined in the time space can be analyzed in the frequency space. For the signal with simple composition in the frequency space, its coefficients before most basis functions are zeros, namely the representation of the signal in the frequency space is sparse. However, the Fourier analysis is merely on the frequency space, ignoring the information on time space. The natural signals such as music (audio) and video contain information both in time and frequency space. To deal with signals as music, a time-frequency dictionary by windowed Fourier transform is designed with success to analyze the variation of frequency over time, but it is limited to the signals with the same time and frequency resolution. The wavelets are thus proposed with the purpose of overcoming this limitation [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

The wavelet transform, considered to be time-frequency representation of sig- [START_REF] Mallat | A wavelet tour of signal processing[END_REF], that is, the large-amplitude wavelet coefficient is detected in the surrounding of transient structure and singularities [START_REF] Le | Image compression with geometrical wavelets[END_REF]. In image processing, it means the edge and texture features [START_REF] Unser | Texture classification and segmentation using wavelet frames[END_REF]. Based on its properties, the wavelet transform allows to create sparse representations of signals. However, the wavelet basis is predefined and signal-independent, which makes it necessary to construct an orthonormal wavelet basis adapted to each category of signals and the task at hand. For example, the Gabor wavelet models the impulse responses of the visual cortex, curvelet frames for edge detection, bandlet transform designed to adapt to the geometric image regularity, see chapter 5 and chapter 12 in [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for details. The wavelet transform for sparse representation is widely developed and applied with success in compressed sensing [START_REF] Sing | Image representation using 2D Gabor wavelets[END_REF], image denoising [START_REF] Grace | Adaptive wavelet thresholding for image denoising and compression[END_REF][START_REF] Grace | Spatially adaptive wavelet thresholding with context modeling for image denoising[END_REF], multiresolution image representation [START_REF] Minh | The contourlet transform: An efficient directional multiresolution image representation[END_REF], classification [START_REF] Unser | Texture classification and segmentation using wavelet frames[END_REF][START_REF] Chang | Texture analysis and classification with tree-structured wavelet transform[END_REF], etc. Nevertheless, the sparse representation capacity of small dictionaries such as orthonormal wavelet basis is limited. Moreover, there is no off-the-shelf mathematical model for a predefined dictionary to represent any signal, such as natural images that contain for example textures.

In 1996, by studying the properties of receptive field of simple cells in mammalian primary visual cortex, Bruno Olshausen and David Field proposed a learning algorithm for sparse representation of a natural image, which allows to learn an overcomplete dictionary from the image at hand [START_REF] Bruno | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF]. Since then, the overcomplete dictionary learning has become a hot topic in machine learning, as it allows to outperform predefined dictionaries like orthonormal wavelet bases.

Dictionary learning is defined as an optimization problem with respect to two vectors of optimization, sparse code (i.e., the coefficients of the linear model) and dictionary, which makes the problem harder to solve. Generally, researchers take the strategy of reaching the optimal solution by iteratively processing two alternating steps:

• Sparse Coding. It concerns the optimization problem with respect to sparse code by supposing the dictionary is already known. This sub-problem can be diversely formulated by introducing different functions for sparsity controlling, e.g. the 0 -norm function which is even non-smooth and non-convex.

For solving this problem, numerous algorithms are designed such as matching pursuit [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF] and basic pursuit [START_REF] Chen | Examples of basis pursuit[END_REF].

• Dictionary Updating. This step deals with the optimization problem with respect to the dictionary while keeping the sparse code fixed. Unlike the subproblem of sparse coding, the sub-problem of renewing the dictionary is often convex and differentiable. The algorithms representative for updating dictionary contains, the Method of Optimal Direction (MOD) [START_REF] Kjersti Engan | Method of optimal directions for frame design[END_REF] and the coordinate descent algorithm with each atom updated by exploiting the singular value decomposition (SVD) [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF][START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF].

The algorithms of dictionary learning are then developed by combining the algorithm of sparse coding and that of dictionary updating, like K-SVD [AEB + 06, EA06],

proximal method [START_REF] Jenatton | Proximal methods for sparse hierarchical dictionary learning[END_REF], and online dictionary learning method [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]. Furthermore, for certain tasks, the specified dictionary can be leaned by adding an appropriate term of regularization. For instance, a discriminative dictionary is learned by proposing the regularization of classification error [START_REF] Zhang | Discriminative K-SVD for dictionary learning in face recognition[END_REF][START_REF] Mairal | Task-driven dictionary learning[END_REF], and an incoherent dictionary is learned via the introduction of constraints on the off-diagonal entries of Gram matrix of dictionary [START_REF] Li | Incoherent dictionary learning with log-regularizer based on proximal operators[END_REF].

The sparse representation by learning a dictionary has thus gained great success in the field of signal processing, image processing, pattern recognition, and computer vision [MCW05a, MBP + 14]. Specifically, sparse representation with the learned dictionary is successfully used to accomplish tasks like image denoising, deblurring, inpainting, face recognition, visual tracking, and classification (see [MBP + 14] and therein).

Contributions of the research

This monograph concentrates on the study of dictionary learning for sparse representations. More precisely, it corresponds to the optimization problem of minimizing the reconstruction error with limiting the sparsity of decomposition coefficients.

As is known, originally, the sparsity is explicitly expressed by using the 0 -norm.

This function can exactly control the sparsity level, however, with a shortcoming of non-convexity and non-differentiability, which makes it the obstacle to overcome. The existing methods for solving this problem can be roughly grouped in two major classes: Greedy algorithms attains the solution by iteratively providing suboptimal solutions, such as matching pursuit (MP) and its variants [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF][START_REF] Bo | Hierarchical matching pursuit for image classification: Architecture and fast algorithms[END_REF],

or gradient descent based algorithms such as Iterative Hard Thresholding (IHT) algorithm [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF] or proximal method [START_REF] Bao | L0 norm based dictionary learning by proximal methods with global convergence[END_REF]. The second class corresponds to methods that relax the 0 -norm by replacing it with the 1 -norm, which is still nonsmooth but convex and continuous. This leads to a classical problem, often called the problem of LASSO [START_REF] Meinshausen | Lasso-type recovery of sparse representations for high-dimensional data[END_REF] or Basis Pursuit [START_REF] Chen | Examples of basis pursuit[END_REF] (technically, LASSO or Basis pursuit principle are the methods for solving the 1 -norm based sparse representation problem). However, now in the research community, it is called directly LASSO problem or Basis pursuit problem [MBP + 14]). This relaxation simplifies the problem and makes it possible to transform the original problem into a certain standard optimization problem so that one can exploit optimization techniques such as the simplex method [START_REF] Nocedal | Numerical optimization[END_REF]. However, the 1 -norm cannot always guarantee the required sparsity level [START_REF] Atamturk | Sparse and smooth signal estimation: Convexification of l0 formulations[END_REF] and even the sparse representation is produced but with shrinkage [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

When referring to the problem of dictionary learning, as aforementioned, it corresponds to two vectors of optimization and the strategy is to solve it by iteratively alternating two steps, sparse coding and dictionary updating. This will raise the difficulty. The purpose of our research is thus to find a method to learn the dictionary by solving exactly the 0 -norm based problem. The main contributions of our work can be summarized as follows.

We work on the 0 -norm constrained dictionary learning problem. For reaching the optimal solution, we take the frequently used strategy of iterative processing two alternating steps: sparse coding and dictionary updating. Specifically, in the step of sparse coding, we reformulate the original problem as a problem of mixed integer quadratic programming (MIQP) and solve it by the optimization method for MIQP problem, without using any greedy algorithm, or relaxation. The second step takes advantage of the coordinate descent algorithm (precisely, each atom is updated by SVD). Hence, it is the first time that the dictionary is learned with MIQP for sparse coding.

We propose two techniques for accelerating the algorithm of exact sparse coding. One is to offer an initialization when solving the MIQP problem. This value obtained by using the proximal method is considered to be an approximation of the optimal solution of the original problem. The other is to relax the problem by achieving a convex envelop of the region defined by all the constraints. With the help of these two acceleration methods, the MIQP is not limited to dealing with synthetic data problems but real data. For the demonstration, the accelerated MIQP (AcMIQP) is applied for denoising well-known images (Barbara, Lena etc). Furthermore, the results show the superiority of the dictionary learning algorithm with exact sparse coding method AcMIQP comparing to the methods of proximal and OMP in processing the image with high noise level.

Furthermore, we study incoherent dictionary learning. The coherence is regarded as one important quality of dictionary [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. The uniqueness of the solution of the problem of sparse representation is proved to rely on the coherence of dictionary [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF], and even some (greedy) algorithms can be precisely executed under the condition of incoherence of dictionary [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. Untill now, the methods for learning incoherent dictionary are either by adding a decorrelation step following the dictionary update, such as INK-SVD [START_REF] Mailhé | INK-SVD: Learning incoherent dictionaries for sparse representations[END_REF], or by formulating the optimization problem by introducing a regularization term on the dictionary coherence [START_REF] Ramírez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF]. Unlike these methods, we seek to solve the exact incoherent dictionary learning problem, with explicit constraints on the dictionary and the 0 -norm sparse code. The corresponding dictionary learning algorithm operates two steps, as in the classical problem, sparse coding and dictionary updating at each iteration.

In this problem, the sparse coding is solved exactly with the AcMIQP method. For updating the dictionary, we exploit methods of the augmented Lagrangian and the proximal alternating linearized minimization. Moreover, the convergence of the algorithm is guaranteed.

The incoherent dictionary learning algorithm is used for estimating the relationship of image reconstruction with coherence of dictionary. It is worthy noting that with the proposed incoherent dictionary learning algorithm, the best results are obtained comparing with the methods INK-SVD [START_REF] Mailhé | INK-SVD: Learning incoherent dictionaries for sparse representations[END_REF] and the incoherent dictionary learning algorithm by iterative projection and rotation (IPR) [START_REF] Barchiesi | Learning incoherent dictionaries for sparse approximation using iterative projections and rotations[END_REF]. Moreover, the fact is proved, to our best knowledge, for the first time, that an appropriate higher incoherence favors the image reconstruction.

Outline of the monograph

The monograph focuses on the study of dictionary learning with 0 -norm for sparsity promoting. The rest of the monograph is organized as follow.

Chapter 2 introduces the mathematical framework of the sparse model. We talk about the three main functions 0 -norm, 1 -norm and p -norm (0 < p < 1), which are frequently used to control the sparsity. By analyzing the three functions in ability of sparsity controlling, convexity and influence on the accuracy of sparse code, we

show the interest of studies on the 0 -norm based sparse representation problems.

Chapter 3 presents the state-of-the-art algorithms for the sparse representation.

These algorithms are divided into two parts: algorithms for sparse coding and those for dictionary learning. Specifically, we detail greedy algorithms, relaxation methods, gradient descent, proximal method and other optimization methods.

Chapter 4 develops a dictionary learning method in the same spirit of the K-SVD algorithm. The point of innovation is that an exact optimization technique named mixed integer programming with quadratic objective function (MIQP) is proposed for the sparse coding phase. Furthermore, two optimization techniques are introduced in MIQP to accelerate the algorithm, which makes the algorithm feasible in real image processing, such as image denoising.

Chapter 5 focuses on the study of coherence of dictionary, which is one of the dictionary's important properties with proved theoretical results. Our work studies explicitly the relationship between coherence of dictionary and image reconstruction performance. For this purpose, we study the dictionary learning problem with a predefined coherence level and an exact 0 -norm promoting sparsity. Then, an augmented Lagrangian based algorithm combined with proximal alternating linearized minimization is proposed. The model is finally applied in image reconstruction experiments, whose results confirm the theoretical conclusion.

Chapter 6 gives the conclusion and provides some future work.

Publication

The main contributions presented in this monograph can also be found in our publications:

Journal Papers

Notations

In this section, we give some mathematical notations used in this monograph.

The data is considered to be real-valued. Matrices are denoted by an uppercase letter, vectors by a lowercase letter in bold and scalars by a lowercase letter. For example, X ∈ R m×n denotes a matrix with each column represented by the indexed vector x i with the index i ∈ {1, . . . , n}. The j t h entry of x i is written as x i , j . A subset of integer is represented by SI ⊂ {1, . . . , n}, then a matrix X SI is a submatrix whose columns are formed by the columns of X indexed by all elements in SI. The complement of SI is represented by SI c . We use X T to represent the transpose of X. If m = n and X is non-singular, the inverse of X is denoted by X -1 . More generally, for X with m = n, X has its pseudo-inverse X + , also known as Moore-Penrose inverse. t r (X)

indicates the trace of a square matrix X, which is equal to the sum of all entries on its main diagonal. The Frobenius norm of the matrix X is expressed as

X 2 F = i , j x 2 i , j .
For two matrices A ∈ R m×n and B ∈ R m×n , the Frobenius inner product of the two matrices is written as

〈A, B〉 = t r (A T B).
Given a function of X, denoted by f : R n×m → R, its gradient represented by ∇ X f is a matrix of size n × m where the entry in position (i , j ) is calculated via

(∇ X f ) i , j = ∂ f ∂x i , j .
Similarly, a vector x ∈ R m is composed by m elements, its i t h entry is denoted by x i . The 0 -norm of x, represented by x 0 , indicates the number of nonzero elements in x. The 1 -norm of x is the sum of its entries' absolute value, denoted

x 1 = m i =1 |x i |.
The p -norm with is defined by

x p = m i =1 x p i 1 p .
Specially, when p = 2, it is called Euclidean norm. For a function of vectors, written as f : R m → R, the gradient of the function is denoted by ∇ x f . If the function f is not differentibale at x, we calculate the Fréchet subdifferential of f , represented by ∂ f (x) and defined as follows.

Definition 1. ([Kru03]) Let f : R n → R ∪ {+∞} be a proper lower semicontinuous function. For each x ∈ dom f (where dom f = {x ∈ R n | f (x) < +∞} ), the Fréchet (or regular) subdifferential of f at x is ∂ f (x) = x * ∈ X | lim inf z→x f (z)-f (x)-〈x * ,z-x〉 z-x ≥ 0 , (1.1)
where X denotes the topological dual space of R n .

In linear algebra, the restricted isometry property (RIP) characterizes matrices that are nearly orthonormal. Supposing D a matrix of size n×m and k an integer, the restricted isometry constant (RIC) of the matrix D is defined as the minimal value Given a signal x ∈ R n and a matrix D ∈ R n×m , providing a sparse representation consists in finding the solution α α α ∈ R m of the system of linear equation with a regularization term that promotes sparsity, that is arg min

δ k ∈ (0, 1) that (1 -δ k ) x 2 2 ≤ Dx 2 2 ≤ (1 + δ k ) x 2
α α α Ω(α α α) subject to x = Dα α α, (2.1)
where Ω(α α α) measures the degree of sparsity of the solution α α α where smaller values of R(α α α) indicate more sparse solution.

CHAPTER 2. MATHEMATICAL FRAMEWORK OF SPARSE REPRESENTATION

In practice, the dictionary D = [d 1 , . . . , d i , . . . , d m ] is defined in the set D = {D | d i 2 = 1, ∀i = 1, . . . , m}, which allows to avoid too large values in α α α. Moreover, one must assure that the matrix D can span the entire space R n to be sure of the existence of the solution to problem (2.1). In spite of that, it is still difficult to give the condition defining the uniqueness of its solution.

The choice of the function Ω(α α α) is diverse. Some favor functions explicitly promoting sparsity, however, making the problem hard to deal with. Others operate relaxation by taking advantage of convex function to make the problem easy to solve. For instance, the p -norm functions, especially for p = 0 and p = 1, are frequently studied [START_REF] Xu | Unnatural l0 sparse representation for natural image deblurring[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Additionally, the family of functions such as

R(α α α) = (1 -exp(|α i |)), Ω(α α α) = log(1 + |α i |) and Ω(α α α) = |α i |/(1 + |α i |)
are also used to promote sparsity [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF]. These will be introduced in this chapter.

The 0 based sparse representation

Introduction of 0 -norm function

The 0 -norm function is defined based on the concept of limit, which is:

α α α 0 = lim p→0 α α α p p = lim p→0 |α i | p = card(α α α), (2.2)
where card(α α α) indicates the number of non-zero entries in the vector. However, despite the name of 'norm' is given, the 0 -norm is not a strict 'norm' function , since it does not satisfy all the three conditions: 1) the triangle inequality, 2) absolute homogeneity and 3) being positive definite. Indeed, the 0 -norm function, defined from the domain of vectors R m to the naturals N, obeys the triangle inequality, that is, for any arbitrary two vectors u, v ∈ R m , we have

u + v 0 ≤ u 0 + v 0 .
The third condition of positive definiteness is satisfied as well. For a vector v ∈ R m , if the equation v 0 = 0 holds, then it implies that v = 0. For the homogeneity condition, it is easy to see that it is not satisfied, because the equation av 0 = |a| v 0 is not true for any scalar number a ∈ R and vector v ∈ R m .

Even the 0 -norm function is misleading in a sense, it presents powerful capacity in promoting sparsity. Figure 2.1 illustrates some p -norm functions by setting p = is convex and differentiable which makes the problem easy to solve, but loosing the sparsity. For p = 1, the function is convex but not differentiable at α α α = 0 0 0, making it of interest as discussed in Section 2.2 [Tib96, MBPS09, CD95]. When 0 < p < 1, the function becomes non-convex but continuous, and the smaller the p is, the closer the norm will be to the 0 -norm function. Therefore, an p -norm function with 0 < p < 1 can be used to approximate the 0 -norm function [START_REF] Rinaldi | Concave programming for minimizing the zero-norm over polyhedral sets[END_REF].

The 0 -norm function, however, is non-convex and discontinuous. The optimization problem based on 0 -norm is NP-hard [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. It is difficult to solve it exactly. On the other hand, the 0 -norm explicitly measures the sparsity of the representation. Thus, for all p -norm, the 0 -norm shows the best performance in terms of sparsity control. Its proof is given in [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF] (see page 12 for details). (2.3)

The constraint x = Dα α α gives a strict equality relation. In reality, caused by distortion or noise perturbation from the device or information loss during the transmission,

x is not exactly the pure signal, but with some noise, which can be modeled:

x = Dα α α + e, (2.4)
where e ∈ R m means the unfitness noise. Consequently, the problem (2.3) is reformulated by the quadratic inequality constraint rather than the equality constraint, that is, arg min

α α α α α α 0 subject to x -Dα α α 2 2 < , (2.5)
where is a parameter determining the tolerance of data fitting. Intrinsically, the value of depends on the noise level of the signals. This formulation that assumes knowing the noise level is usually used in image denoising [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], and seldom used in classification [FNZ + 15]. Nevertheless, given a signal, seems hardly to be well estimated.

Another formulation consists in minimizing the reconstruct error in the feasible region defined by the 0 -norm constraint, arg min

α α α 1 2 x -Dα α α 2 2 subject to α α α 0 ≤ k, ( 2
α α α 1 2 x -Dα α α 2 2 + λ α α α 0 , (2.7)
where λ > 0 is a hyperparameter balancing the minimization of the reconstruction error and the maximization of sparsity. Nikolova has proved the relationship between the k-sparsity problem (2.6) and its regularized form (2.7) [START_REF] Nikolova | Relationship between the optimal solutions of least squares regularized with l0-norm and constrained by k-sparsity[END_REF]. Specifically, when the choice of k corresponds to a λ in an interval (bounds depending on the reconstruction error decrease with the representation's sparsity), then the optimal solutions to the two formulations become equivalent. Compared with the k-sparsity problem (2.6), the regularized problem (2.7) and its variants are more frequently applied in image processing [START_REF] Mancera | L0-norm-based sparse representation through alternate projections[END_REF] and machine learning [START_REF] Bao | L0 norm based dictionary learning by proximal methods with global convergence[END_REF][START_REF] Yuan | Visual classification with multitask joint sparse representation[END_REF].

Some approximations of the 0 -norm function

As aforementioned, the non-convexity and discontinuity of the 0 -norm cause the corresponding optimization problems to be NP-hard. One feasible method is to approximate the optimal solution by a convergent series emerging from an iterative process. Proposed resolution methods include the greedy algorithm [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF] and the proximal method [START_REF] Bao | L0 norm based dictionary learning by proximal methods with global convergence[END_REF].

Another idea is inspired by reformulating the optimization problem with a smooth function that is not exactly the 0 -norm, but can exceedingly approach it, that is

α α α 0 ← φ(α i ), (2.8)
where φ : R → R + is a smooth function that can measure if the entry α i is zero or not.

The most important is that the optimal solution of this problem must coincides with the one reached when applying the 0 -norm. We talk about some representative reformulations of 0 -norm in the following.

The simplest proposition is to approximate 0 -norm by p -norm. By this thinking, F. Rinaldi et al. proposed two concave formulations tuned by some hyperparameters [START_REF] Rinaldi | Concave programming for minimizing the zero-norm over polyhedral sets[END_REF]. The first formulation is

φ κ (α i ) = (α i + κ) p , (2.9)
where 0 < p < 1 and 0 < κ are two hyperparameters. A variant of (2.9), is

φ κ (α i ) = -(α i + κ) -p .
(2.10)

These two functions are proved to be equivalent to the 0 -norm function when p and k tend to zero, which is based on the fact that lim p→0,κ→0

(α i + κ) p = lim p→0,κ→0 -(α i + κ) -p = α α α 0 .
Thus the solution of the original problem can be reached by solving the reformulated problem, by applying classical optimization techniques, e.g. Frank-Wolfe algorithm [START_REF] Rinaldi | Concave programming for minimizing the zero-norm over polyhedral sets[END_REF].

In [START_REF] G Hosein Mohimani | Fast sparse representation based on smoothed l0 norm[END_REF], Mohimani et al. introduced a function of zero-mean Gaussian family to estimate the 0 -norm. The expression of the function is

φ σ (α i ) = 1 -exp(-α 2 i /2σ 2 ), (2.11)
where σ denotes a bandwidth parameter. This function is smooth and differentiable, which simplifies the problem; however, the effect of σ becomes significant.

To overcome this issue, the proposed algorithm decreases iteratively the value of σ, allowing to exclude the dependency of the solution on the bandwidth value. Hence, the penalty at each coordinate λ α i 0 can be estimated by the so-called continuous exact 0 penalty

φ( d i 2 , λ, α i ) = λ - d i 2 2 2 |α i | -2λ d i 2 2 δ |α i |≤ 2λ d i 2 (α i ), (2.12)
where, given a set C, δ C represents the indicator function

δ C (α) = 1 if α ∈ C 0 else. ( 2
α α α 0 = m i =1 (1 -z i ) under constraints        α α α = α α α + -α α α - z T (α α α + +α α α -) = 0 α α α + ,α α α -≥ 0 and 0 ≤ z ≤ 1 1 1 m , (2.15)
where α α α + and α α α -are the non-negative and non-positive parts of α α α respectively. The second formulation, called half-complementary, takes the same strategy by only changing the above constraints by

0 ≤ z ≤ 1 1 1 m z •α α α = 0, ( 2

Existence and uniqueness of the solution

Hereinabove, we know what the 0 -norm function is, how the problem is formulated with 0 -norm function in sparse representation and some reformulations of the 0 -norm optimization problem. There is a precondition of all these discussions, namely the existence, and more strictly, the uniqueness, of the solution to the problem, which gives sense to the study.

The existence of a solution to the undetermined problem (2.1) is guaranteed under the condition that D spans the feasible solution space. For the optimization problem with regularization (2.7) or the other two constrained formulations (2.6) and (2.5), it is clear that we can find a vector α α α in the feasible region with an objective value, denoted by F( α α α), that satisfies F( α α α) < ∞. Besides, the objective functions are all bounded below. Hence, there must be a solution to these sparse problems.

Nikolova stated in Theorem 4.4 in [START_REF] Nikolova | Description of the minimizers of least squares regularized with l0-norm: Uniqueness of the global minimizer[END_REF] that not only the optimization problem

(2.7) has a global optimal solution, but also the optimal solution is unique. The above theorem offers the sparsity condition for the uniqueness of the optimal solution. In fact, the Spark of a matrix is quite difficult to obtain, which restricts the usage of the conclusion. Nikolova provided in [START_REF] Nikolova | Description of the minimizers of least squares regularized with l0-norm: Uniqueness of the global minimizer[END_REF] the proof of uniqueness of the minimizer by defining the conditions on D and x. Its conclusion for unique global minimizer is described in Theorem 5.6 in [START_REF] Nikolova | Description of the minimizers of least squares regularized with l0-norm: Uniqueness of the global minimizer[END_REF] (For the completeness of the condition expression, the detail is not listed here).

The 1 based sparse representation

The 1 -norm function is a convex and continuous function, and is the closest convex approximation of the 0 -norm, which explicitly expresses the sparsity.

There is no evidence in providing the relationship between 0 -norm and 1norm in sparsity promoting.

The p based sparse representation

The p -norm function with p < 1, as aforementioned, is closer to the 0 -norm, than 1 -norm. This function shows some similar properties to 0 -norm such as its nonconvexity. Besides, it is non-differentiable at 0. As a consequence, minimizing it is a hard problem to solve. That is why the p -norm when p < 1 is less popular than the case where p = 0 and 1.

The optimization problem with p -norm for penalty can be easily formulated

arg min 1 2 x -Dα α α 2 2 + λ α α α p , (2.20)
where the trade-off parameter λ functions the same as that in the 0 -norm or 1norm problems. This problem has obviously a solution. Furthermore, it is proved that the optimal solution of the p -norm based problem is equivalent to that of the 0 -norm based problem [START_REF] Fung | Equivalence of minimal l0-and lpnorm solutions of linear equalities, inequalities and linear programs for sufficiently small p[END_REF] under the condition that α α α 0 is upper bounded by

Ul = f 1 (D) p f 2 (D) 1 + f 1 (D) p ,
where f 2 is the function indicating the spark of D, and f 1 is the minimum of the

f 2 (D) t h descent ordered NULL(D) [MCW04, MCW05a].
The p -norm is usually used as an approximation to the 0 -norm to make the problem tractable. Nevertheless, it risks of heavy computation cost due to its nonconvexity and non-differentiability. To overcome these difficulties, another approximation is needed to reach the solution of the problem, such as

α α α p p = (|α| 2 + ) p/2 , (2.21)
where ≥ 0 is a smoothing parameter [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. With this approximation replacing the p -norm, the original optimization problem (2.20) can be solved by half-quadratic regularization method or iterative method (see [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF][START_REF] Chen | Local and global optimality of lp minimization for sparse recovery[END_REF] and therein for details).

Besides the p -norm with p < 1, some other norms are used, such as the 2-1 -norm used for face recognition [START_REF] Shi | Face recognition by sparse discriminant analysis via joint l2, 1-norm minimization[END_REF] and feature selection [START_REF] Nie | Efficient and robust feature selection via joint l2, 1-norms minimization[END_REF]. The joint use of several norms has been also considered, such as the combination of the 0 and 1 penalties [START_REF] Liu | Variable selection via a combination of the l0 and l1 penalties[END_REF].

Conclusion

This chapter introduced the mathematical framework of sparse representations. We presented the most-known norm functions for sparsity promoting: 0 -norm, 1norm and p -norm with 0 < p < 1. The characteristics of each norm function were discussed. Of particular interest is the 0 -norm, which counts the number of nonzeros of a vector and thus can strictly control the sparsity. Nevertheless, its undeniable shortcoming is that the function is non-smooth and non-convex, making the problem NP-hard. The 1 -norm is regarded as a relaxation of the 0 -norm for the sparsity-promoting. Unlike the 0 -norm, this function is convex and differentiable at all values except at the zero. The price to pay is the risk of missing the optimal sparsest solution. The p -norm with 0 < p < 1 is closer to the 0 -norm function than the 1 -norm, thus it has stronger sparsity control ability than the latter; however, it is a concave function. For each norm function, we detailed the formulations of the optimization problems of sparse representation and also gave a problem transformation survey.

Next chapter presents algorithms for sparse representation learning. In the previous chapter, we presented the mathematical formulations of sparse representation. It is revealed that the problem can be defined by inducing the sparsity with the 0 -norm, 1 -norm or p -norm (with 0 < p < 1) functions. Owing to its power in controlling explicitly the number of non-zero coefficients, the 0 -norm is theoretically the first choice to build sparse representations. However, the 0norm is non-smooth and non-convex, which makes the corresponding optimization problem NP-hard. To overcome the shortcoming of the 0 -norm, it is usually replace with the 1 -norm, which is the closest convex norm function to 0 -norm.

Moreover, it is proved that, in some conditions (e.g. when the sparse code is quite sparse), the optimal solution of the 1 -norm based optimization problem coincides with that of the 0 -norm based problem.

This chapter presents the state-of-the-art algorithms for sparse representation learning. It considers two major problems: sparse coding and dictionary learning.

Sparse coding is the keystone in sparse representation learning. It consists in estimating the sparse code (i.e., the coefficients of the linear model) while the dictionary is already known. This optimization problem can be formulated by introducing sparsity-promoting constraints, mainly the 0 -norm and the 1 -norm defined in the previous chapter. Since the 1990's, extensive research efforts have been made to address to solve the sparse coding optimization problem. Section 3.2 in this chapter presents the prime algorithms for sparse coding, such as greedy algorithms [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] whose major representative is the Matching Pursuit [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF], relaxation methods with Basis Pursuit [START_REF] João | Distributed basis pursuit[END_REF] and LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], Gradient Descent algorithms [START_REF] Kim | A gradient-based optimization algorithm for lasso[END_REF] and Iterative Thresholding methods [START_REF] Kowalski | Thresholding rules and iterative shrinkage/thresholding algorithm: A convergence study[END_REF].

In Sections 3.3 and 3.4, we will focus on the problems of dictionary learning. In complement to estimating the sparse code as in sparse coding problems, dictionary learning aims at estimating also the optimal dictionary, thus improving the performance of the sparse representation. The dictionary learning aims at finding jointly the optimal solution of the sparse code and of the dictionary. The resulting optimization problem is non-convex and NP-hard, and is more difficult than solving the single sparse coding problem. To overcome this issue, a good resolution strategy is to iteratively alternating two steps: sparse coding and dictionary updating.

The sparse coding problem, as aforementioned, handles the estimation of the decomposition coefficients with a fixed dictionary. While fixing these coefficients, dictionary updating aims at estimating the optimal dictionary. In its general form, the

resulting

Sparse coding and dictionary learning

Given a matrix X = [x 1 , . . . , x i , . . . , x ] ∈ R n× of signals of dimension n, and a (predefined or learned) dictionary dictionary D = [d 1 , . . . , d m ] ∈ R n×m the problem of sparse representation consists in finding a matrix A = [α α α 1 , . . . ,α α α ] ∈ R m× of decomposition coefficients that satisfies X ≈ DA. The columns of the dictionary, i.e., d j for j = 1, . . . , m, are called atoms.

In signal or image data processing, the set of signals is typically larger than its dimension, namely n. The sparse representation can be obtained by solving the following optimization problem:

min α α α i ∈R m 1 i =1 1 2 x i -Dα α α i 2 2 + λΩ(α α α i ) . (3.1)
The first term 1 2 x i -Dα α α i 2 2 is the reconstruction error with . 2 being the Euclidean norm. The second one includes the regularization term Ω(α α α i ) to enforce sparsity.

The regularization parameter λ > 0 controls the trade-off between data fitting and sparsity of α α α. For the sake of clarity of this monograph, the reconstruction error is measured with the square loss; generalization to other loss functions, such as the logistic and hinge losses, is straightforward [ŞE13]. Generally, the regularization function Ω is associated to a norm that promotes sparsity and its formulation depends on the task at hand [BJQS14, AEB + 06]. The natural definition of Ω to promote sparsity is, as mentioned in Chapter 2, the 0 -norm.

With the assumption of independence of the signals, the optimization problem (3.1) can be divided into independent subproblems min

α α α i ∈R m 1 2 x i -Dα α α i 2 2 + λΩ(α α α i ). (3.2)
When the 0 -norm is adopted to promote the sparsity, this problem becomes nonconvex, non-smooth, and thus NP-hard, as analyzed in Chapter 2.

The dictionary plays a crucial role in sparse representation [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF][START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF]. Nevertheless, there is no predefined dictionary that can satisfy all tasks in processing natural signals and images. For all these reasons, learning the dictionary from the data is a main building block in sparse representation. And in general, D is an overcomplete dictionary, that is to say n < m, while the situation n > m is allowed for some discrimination tasks [MPS + 09]. To prevent the 2 -norm of dictionary's atoms from being arbitrarily large which leads to arbitrarily small decomposition coefficient in A, the dictionary D is supposed without loss of generality, to satisfy

D = D ∈ R n×m | d j 2 ≤ 1, ∀ j = 1, . . . , m . (3.3)
Beyond these considerations, the dictionary learning problem can be written in the following general form, min

α α α i ∈R m ,D∈D 1 i =1 1 2 x i -Dα α α i 2 2 + λΩ(α α α i ) . (3.4)
This problem of estimating simultaneously A and D is non-convex and belongs to NP-hard problems. It is often solved via an alternating strategy: 1) fixing D and finding the sparse coefficients A, this becomes the problem of sparse coding as aforementioned; 2) fixing A and searching for the solution D, this is the procedure of dictionary updating. Several popular algorithms for sparse coding and dictionary updating are described in the following, as well as the most known combinations of these algorithms.

Algorithms for sparse coding

The optimization problem of sparse coding (3.2) consists in finding the sparse representation α α α ∈ R m of x over the dictionary D, which can be specifically formulated in different forms.

The k-sparse representation problem with maximal number of non-zero entries fixed is modeled by arg min

α α α 1 2 x -Dα α α 2 2 , α α α ∈ S k , (3.5)
where

S k = {v ∈ R m | v 0 ≤ k} means the set of k-sparse vectors.
By considering the constraint on the reconstruction error, the sparsest representation can be obtained via the problem arg min

α α α Ω(α α α) subject to x -Dα α α 2 2 < , (3.6)
where Ω(α α α) is a function of α α α which measures the sparsity of the vector. These constrained optimization problems can be addressed by considering the minimization of the regularized reconstruction error, arg min

α α α 1 2 x -Dα α α 2 2 + λΩ(α α α), (3.7)
where λ introduced previously balances the effects of data fitting and sparsity.

As aforementioned, these three formulations are used for different tasks and they can achieve the equivalent global optima under certain conditions. In contrast, the strategies for tackling the problems of different formulations are quite different. The 0 -norm based problem is non-convex and non-continuous, which makes it hard to be resolved directly. To overcome these difficulties, two major strategies have been widely investigated. The first one use greedy algorithms, such as

Matching Pursuit [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF], Orthogonal Matching Pursuit [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], subspace pursuit [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], and even proximal method [START_REF] Parikh | Proximal algorithms[END_REF] to let it be tractable. The second strategy is relaxation, since the 1 -norm based problem, also sometimes called LASSO problem, becomes convex which is solved by Basis Pursuit [START_REF] Chen | Examples of basis pursuit[END_REF], by gradient descent algorithm [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF], by homotopy algorithm [START_REF] Garrigues | An homotopy algorithm for the lasso with online observations[END_REF] or by thresholding method [START_REF] Meinshausen | Relaxed lasso[END_REF] and very recently, a novel strategy is emerging, by reformulating the problem as MIQP, which makes it possible to apply some optimization methods like the branch-and-bound method [START_REF] Atamturk | Sparse and smooth signal estimation: Convexification of l0 formulations[END_REF][START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF] or relaxing the integer constraints to be continuous [BBF + 16] and so on. The two major strategies with the main algorithms for sparse coding are described in the following, while the novel strategy is presented in next chapter.

Greedy algorithms

A greedy algorithm solves an optimization problem by searching for the optimal solution of a subproblem at each step. For the sparse coding problem with 0 -norm for sparsity promoting, it selects at each step an atom over which the residual has the maximal projection. Greedy algorithms provide a good sparse approximation [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] rather than the optimal sparse solution. In this part, we present the most known greedy algorithms: matching pursuit and its variants.

Matching pursuit

Matching pursuit (MP) is proposed in [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF] to deal with the problem of signal decomposition. The algorithm is inspired by the fact that, giving a signal x and a dictionary D with each atom of unit norm, denoted by d i with i = 1, . . . , m, then the signal decomposition by orthogonal projection onto some selected atom d t 1 is given by

x = x T d t 1 d t 1 + r 1 , (3.8)
where r 1 is the representation residual orthogonal to d t 1 ,the latter being selected based on the residual. This leads to the result,

x 2 = |x T d t 1 | 2 + r 1 2 . (3.9)
Then it is conducted that by minimizing r 1 2 , namely maximizing |x T d t 1 |. The signal x can be maximally approximated in the space spanned by d t 1 . By successively doing k times orthogonal projections on selected atoms of D, the signal can be approximated by

x = k j =1 x T d t j + r k , (3.10)
where the sequence t j contains the indices of the selected atoms. With lim k→∞ r k 2 = 0, the signal x can be approximated with high precision. It induces also the convergence of the algorithm.

The Matching Pursuit algorithm is exhibited in Algorithm 1. It is worthy to notice that MP algorithm guarantees in each iteration the orthogonality between the i t h selected atom d t i and the residual r i . Nevertheless, the space spanned by the selected atoms [d t 1 , . . . , d t i ] after i approximations is not guaranteed to be orthogonal to r i . This risks of getting the sub-optimal solution in each step, which may cause large error after a finite number of steps [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF].

Orthogonal matching pursuit

The Orthogonal Matching Pursuit (OMP) algorithm [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] is proposed to overcome the drawback of MP algorithm. Intrinsically, OMP inherits the spirit of MP, which selects an atom at each step to decrease the residual. Thus, OMP is still a greedy algorithm. However, unlike MP, OMP recalculates the composition in the space spanned by the selected atoms, which guarantees always the orthogonality of

Algorithm 1 Matching Pursuit Algorithm (MP)

Input: The given signal x, the fixed dictionary D, the stop criterion of sparsity k, the stop criterion of reconstruct error . Output: Sparse approximation α α α function SPARSE CODING Initialization i = 1, the residual r 0 = x, the approximation α α α = 0 and the initial solution support T 0 = while c ar d (T i -1 ) < k + 1 or x -Dα α α 2 2 > do Selecting the atom over which the residual r i -1 has the maximal magnitude of orthogonal projection ĵ = arg max

j =1,...,m |〈r i -1 , d j 〉|; Updating support T i = T i -1 ∪ ĵ
Updating the residual after i t h projection

r i = r i -1 -〈r i -1 , d ĵ 〉d ĵ ;
Updating the coefficient corresponding to the ĵ t h atom

α α α( ĵ ) = 〈r i -1 , d ĵ 〉; i = i + 1; end while end function
representation residual and the spanned space.

The details of OMP algorithm are listed in Algorithm 2. Compared with MP algorithm, OMP algorithm carries out one more operation (3.11) to lead to the best approximation over the selected support [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. Tropp gave the theoretical support that, when the signal is sparse enough, OMP performs well in signal recovery and approximation [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF], even with random measurement matrix such as Gaussian and Bernoulli [START_REF] Joel | Signal recovery from random measurements via orthogonal matching pursuit[END_REF].

It is observed that OMP needs no more than k iterations to achieve the approximate solution. Nevertheless, the computing cost in atom identification and coefficient update (3.11) is unnegligible. A progressive Cholesky process is adopted to handle the problem (3.11) and some implementation skills are used to accelerate the algorithm [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF].

Furthermore, based on the OMP algorithm, some variants are proposed to reduce computation complexity by optimizing the identification step. The Generalized Orthogonal Matching Pursuit algorithm (gOMP) [START_REF] Wang | Generalized orthogonal matching pursuit[END_REF] speeds up the com-

Algorithm 2 Orthogonal Matching Pursuit Algorithm (OMP)

Input: The signal for sparse coding x, the known dictionary D, the stopping criterion of sparsity k or that of reconstruct error . Output: The sparse representation α α α.

1: function SPARSE CODING 2:

Initialization i = 1, the residual r 0 = x, the solution support T 0 = , the sparse approximation α α α = 0 3:

while i < k + 1 or r i -1 2 2 > do 4:
Finding the atom the highest correlated with the current residual, that is

t i = arg max j ∈{1,...,m}\T i -1 |〈r i -1 , d j 〉|;
where j ∈ {1, . . . , m} \ T i -1 means that j ∈ {1, . . . , m} and j ∈ T i -1 .

5:

Adding the new atom index into the support

T i = T i -1 ∪ t i ; 6: Updating coefficient vector α α α T i = (D T T i D T i ) -1 D T T i x; (3.11)
Here we recall the notation of α α α T i and D T i representing the sub-vector and submatrix indexed by the elements in the set T i .

7:

Updating residual

r i = x -D T i α α α T i ; 8: Increment i = i + 1 9:
end while 10: end function putation by selecting more than one atom at each step. Through the method, a k-sparse signal can be approximated by less than k iterations. However, gOMP is not suitable for signals that are not strictly sparse.

The Regularized Orthogonal Matching Pursuit algorithm (ROMP), as is called, is regularized version of OMP [START_REF] Needell | Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit[END_REF]. ROMP takes the same strategy as in gOMP to recover sparse signals. At every iteration, all the atoms with comparable coordinates are selected. This process named regularization is realized by:

|α α α i | ≤ 2|α α α j | ∀i , j ∈ Λ s ,
where Λ s denotes a subset of support Λ at i t h iteration, which is formed by the k atoms with largest correlation to the sparse signal. ROMP gains success in dealing with sparse signal reconstruction under the Restricted Isometry Condition (RIC).

Even for natural signals, which are not strictly sparse, the algorithm shows stability in sparse approximation. However, the level of sparsity k can be missed with high probability [START_REF] Wang | Generalized orthogonal matching pursuit[END_REF].

Similarly

Compressive sampling matching pursuit (CoSaMP) and subspace pursuit (SP)

The Compressive Sampling Matching Pursuit (CoSaMP) [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF] and the Subspace Pursuit (SP) [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] are two other sparse coding algorithms developed on OMP.

Nevertheless, the two algorithms are inspired by the fact that the energy in each component of u = D T Dα α α approximates the energy in the corresponding component of α α α, when D satisfies the RIC, namely the largest k elements in u point to the largest k entries in α α α. Thus, the k support of α α α can be predicted by u, which is called here proxy. Then, by iteratively regulating support to minimize the residual, the ksparsity approximation or recovery problem can be solved.

The pseudocode of the CoSaMP algorithm is given in Algorithm 3. This algorithm, in contrast with gOMP, StOMP and ROMP, can guarantee that the obtained sparsity level is k. To this end, the rigorous condition of RIC on D needs to be established. Thus, signal-to-noise ratio (SNR) reduction occurs in each iteration, which induces the convergence of the algorithm. Moreover, based on analysis on SNR reduction, the number of iterations can be predicted. Even it is undeniable that the

Algorithm 3 Compressive Sampling Matching Pursuit (CoSaMP)

Input: A given signal x, the fixed overcomplete dictionary D, maximal number of iteration N i t er , the sparsity level k. Output: The sparse representation α α α.

1: function SPARSE CODING 2:

Initialization i = 1, the residual r 0 = x, the solution support T 0 = , the sparse approximation α α α = 0 3:

while i < N i t er + 1 do 4:
Calculating 2k largest elements in the proxy u = D T r i -1 , its support is thus obtained T = {t 1 , . . . , t 2k };

where t i = arg max u {1,...,m}\{t 1 ,...,t i -1 } pointing to the i t h largest entry in u.

5:

Merging T and the current support T i to update support

T i = T i -1 ∪ T; 6:
Updating projection value in the new support

α α α T i = (D T T i D T i ) -1 D T T i x; (3.12) 7:
Pruning elements in support to keep only the k atoms retaining the largest correlation with x. The updated support is denoted by T i . The algorithm SP is nearly the same as CoSaMP except for the number of identified atoms [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF]. This algorithm also well performs under RIC with adequate constant requirement. However, in SP, Dai and Milenkovic located only the k largest components in proxy, rather than 2k ones. With this improvement, new theoretical analysis is made in [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] to guarantee the residual reduction in each iteration and the convergence was proved with a limited number of iterations. In this case, the halting criterion can adopt the difference between the residuals before and after each iteration, that is to say, if r i 2 > r i -1 2 , set T i = T i -1 and quit the iteration.

Other matching pursuit

There are other variants of the MP, such as Hierarchical Matching Pursuit (HMP) [START_REF] Bo | Hierarchical matching pursuit for image classification: Architecture and fast algorithms[END_REF], and (fast) Bayesian Matching Pursuit (BMP) [START_REF] Schniter | Fast Bayesian matching pursuit[END_REF][START_REF] Wang | High resolution radar imaging based on compressed sensing and fast Bayesian matching pursuit[END_REF][START_REF] Masood | Sparse reconstruction using distribution agnostic Bayesian matching pursuit[END_REF].

HMP was proposed by Bo et al. for high level image feature extraction [START_REF] Bo | Multipath sparse coding using hierarchical matching pursuit[END_REF][START_REF] Bo | Hierarchical matching pursuit for image classification: Architecture and fast algorithms[END_REF]. In contrast with classical OMP algorithms, HMP is designed with three crucial factors: tree-structured dictionary for sparse coding with algorithm OMP, spatial pyramid max pooling and normalization. HMP contains multi-level and each level operates all the three above-mentioned processes. In this framework, the higher level uses the output feature in the last level as the input data for tree OMP, and its output passes to the next level for further learning. Image feature extracted by HMP is then applied in classification. The performance exhibited in [START_REF] Bo | Hierarchical matching pursuit for image classification: Architecture and fast algorithms[END_REF] proved its advantage even comparing with convolutional neural networks.

BMP regards the sparse representation of a signal as a random variable that satisfies certain probabilistic distribution, rather than a deterministic one. For instance, in [START_REF] Schniter | Fast Bayesian matching pursuit[END_REF], an i.i.d. random variable α α α is drawn from a Q-ary Gaussian mixture parameterized by s = [s 1 , . . . , s m ] T , that is, each entry α i corresponds to a certain Gaussian distribution (s i = q ) with mean µ q and variance σ 2 q ,

α i | s i = q ∼ CN(µ q , σ 2 i ).
With an appropriate Gaussian model (four models were introduced in [SPZ08]:

zeros-means binary prior, nonzero-mean binary prior, zero-mean tenary prior and Q-ary circular prior), zero entries and non-zero entries in α α α can be exactly modeled (s i = 0 ⇒ α i = 0 and s i = 0 ⇒ α i = 0). The probability p(s, x) is estimated by a variable v(s, x) called model selection metric, with p(s, x) = e v(s,x) (see [START_REF] Schniter | Fast Bayesian matching pursuit[END_REF] for details). The determination of the support of x is thus transformed into a problem of finding the Gaussian mixture model with the largest probability. This problem is formulated as a tree search with each node representing a possible mixture s, specifically, the p-level node contains the mixtures with s 0 = p. By iteratively adding a non-zero entry in s via selecting the one with the largest probability, the Gaussian model can be finally reached. Then, the sparse code is computed by

α α α = p(s, x)E(α α α | x, s),
where the expectation can be computed from the Gaussian model assumptions.

p -norm relaxation

By replacing the 0 -norm by another p -norm for sparsity prompting, the optimization problems (3.5), (3.6) and (3.7) are transformed to become continuous, and even convex (for p ≥ 1), and differentiable (for p > 1). The 1 -norm, the nearest norm function to the 0 -norm that is convex, attracted most attention [CDS01, MBPS09, Don06, MKCP17]. The resulting sparse representation problem based on 1 -norm is often solved with the Basis Pursuit [START_REF] Shaobing | Atomic decomposition by basis pursuit[END_REF] and LASSO method [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

Basis pursuit

The Basis Pursuit [START_REF] João | Distributed basis pursuit[END_REF][START_REF] Shaobing | Atomic decomposition by basis pursuit[END_REF] focuses on the linear constrained optimization problem (3.6). As mentioned in [START_REF] Shaobing | Atomic decomposition by basis pursuit[END_REF], the Basis Pursuit is a principle rather than an algorithm. It reformulates the problem by separating the non-negative and the non-positive parts of α α α, namely α α α = α α α + -α α α -with α α α + ≥ 0 and α α α -≥ 0. The optimization problem is thus written by

arg min c T v subject to [D, -D]v = x, v ≥ 0, (3.13)
where c = 1 1 1 2m is the vector of ones of size 2m, and variable v of size 2m is composed by (α α α + ,α α α -). The problem (3.13) is obviously the standard form of linear programming. For denoising, the corresponding optimization problem, formulated with the error constraint (3.6) or regularization term (3.7) [START_REF] Chen | Examples of basis pursuit[END_REF], can be reformulated as a problem of linear programming with quadratic constraints or quadratic programming respectively. In this way, the 1 -norm based problem can be solved via optimization techniques, such as the simplex algorithm and the interior point method.

LASSO

Initialization α α α 0 = (D T D) -1 D T x, defining E = {i | g i T α α α = λ k } where g i is the i t h row of G and setting initial E = {i 0 | g i 0 = sign(α α α 0 )}. 3: do 4:
Finding solution α α α to the problem of minimization, arg min

α α α 1 2 x -Dα α α 2 subject to G E α α α ≤ λ k 1 1 1 |E| (3.14) 5: updating E = E ∪ {i | g i = sign( α α α)} 6:
while α α α 1 > λ k 7: end function is proposed in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. The outline of the algorithm is summarized in Algorithm (4), where the optimization problem (3.14) is a quadratic programming with linear inequality constraints, which can be solved with optimization problem methods such as seeking the solution satisfying the Kuhn-Tucker conditions [START_REF] Bertsekas | Nonlinear programming[END_REF].

In addition to reformulating the 1 -norm regularized least square minimization problem as a classical optimization problem, this problem can also be solved by algorithms such as Least Angle Regression (LARS), homotopy algorithm or some gradient based algorithms [KKK08, GK09](gradient algorithms will be presented in the next subsection). The LARS algorithm updates the estimation iteratively along the least angle direction until all non-zero elements are found [EHJ +

d α d α d α = ∂ α α α 1 = v ∈ R m v i = sign(α i ) if |α i | > 0 v i ∈ [-1, 1] if α i = 0.
Thus, with knowing the support and sign of the entries α α α, the problem parameterized by a certain λ can be easily solved. The Lagrangian method turns the constrained optimization problem into a regularized least-square minimization problem. This optimization problem lies then on the determination of the support and sign of α α α. Homotopy is the algorithm offering the way to find these two important information about α α α. In [START_REF] Dmitry M Malioutov | Homotopy continuation for sparse signal representation[END_REF], the authors regarded the optimal solution to the problem as a function of λ, expressed by α α α(λ). This function shows a piecewise linear continuity, and the support and sign of α α α keep unchangeable in a range [START_REF] Garrigues | An homotopy algorithm for the lasso with online observations[END_REF][START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF]. With this characteristic, the support and sign can be suc- proved the algorithm by rewriting the optimal solution of the problem as a function of an auxiliary variable t and parameter λ, namely α α α(t , λ). Moreover, λ is not directly set but reached by an increasing sequence {λ n }. With these notations, they design a path from α α α at the i t h iteration α α α i = α α α(0, λ i ) to its value at iteration i + 1 with α α α i +1 = α α α(1, λ i +1 ), which is composed of two steps: 1) updating α α α(0, λ i +1 ) from α α α(0, λ i ) with the method mentioned above; 2) fixing λ and varying t from 0 to 1.

Relaxation formulations by p -norm with 0 < p < 1 are also studied in references [FM11, RSS10, SBFA17, ZMZ + 13]. However, the corresponding optimization problem is not convex, which makes it hard to solve.

Gradient descent and iterative thresholding algorithms

Iterative Thresholding methods [START_REF] Kowalski | Thresholding rules and iterative shrinkage/thresholding algorithm: A convergence study[END_REF][START_REF] Blumensath | Iterative hard thresholding and l0 regularisation[END_REF] Initialization i = 1, the coefficient vector α α α 0 = 0.

3:

while x -Dα α α 2 2 > do 4: Gradient Descent to reduce error x -Dα α α 2 2 , α α α m = α α α i -1 + 1 γ D T (x -Dα α α i -1 ). (3.16) 5:
Sorting all entries of |α α α| in descending order, that is,

|α s(1) | ≥ • • • ≥ |α s(m) |. Then, taking the k t h largest one |α s(k) | the threshold. The α α α is thus updated by α α α i = T |α s(k) | (α α α),
where T c (α α α) means that all entries of α α α larger than c in magnitude will be retained while the others are set to 0. 6:

Increment i = i + 1 7:
end while 8: end function

Gradient descent with projection on constraint set

Iterative Hard Thresholding (IHT), also called Gradient Descent with Sparsification (GraDeS) in [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF], creates a sequence of α α α i by iteratively reducing the square error along the gradient direction by a step of length 1/γ followed by retraining the largest k entries of α α α by thresholding. Its framework is described in Algorithm 5.

The IHT algorithm is proved to produce sparse approximation whenever the signal is noisy or not, under the condition that the dictionary D satisfies the RIC with isometric constant δ 2k < 1/3. Furthermore, the algorithm is proved to converge in a limited number of iterations by setting γ = 1 [BD08, BD09, HGT06] or γ = 1 + γ 2k [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF]. Moreover, the update step size γ is allowed to vary according to iteration [START_REF] Foucart | Hard thresholding pursuit: An algorithm for compressive sensing[END_REF]. Foucart [START_REF] Foucart | Hard thresholding pursuit: An algorithm for compressive sensing[END_REF] and Yuan et al [START_REF] Yuan | Gradient hard thresholding pursuit[END_REF] proposed a variant, Hard Thresholding Pursuit (HTP), which combines the IHT method with CoSaMP. At each iteration, following the thresholding step, HTP calculates the maximal projection in the space spanned by the support of α α α i . Compared with IHT, HTP improves the performance on error minimizing and holds on the high-speed computation.

Besides, the IHT is intrinsically equivalent to the proximal method [START_REF] Bao | L0 norm based dictionary learning by proximal methods with global convergence[END_REF][START_REF] Parikh | Proximal algorithms[END_REF] when dealing with a k-sparsity constrained problem (details given in section 4.1.2), although these two methods are developed in two different frameworks.

Iterative thresholding algorithm for regularized sparse problem

For the regularized optimization problem, the solution depends on the value of trade-off parameter λ, which implies that the threshold in the algorithm is a function of λ, namely s(λ). The resolution of this problem is done with the aid of a sur-

α α α,a 1 2 x -Dα α α 2 2 + λ α α α p -Dα α α -Da 2 2 + α α α -a 2 2 , (3.17)
where a is a variable closed to α α α. Developing the problem (3.17), we have arg min

α α α,a j (α 2 j -2α j (a + D T x -D T Da) j + λ|α j | p ) + x 2 + a 2 -Da 2 . ( 3.18)
When p = 0, the objective function is non-differentiable. Considering the problem separately for the cases α j = 0 and α j = 0, then, its solution is of the form

α α α = T λ 0.5 (a + D T (x -Da)),
where T λ 0.5 is the thresholding operator defined as

T λ 0.5 (α) = 0 if |α| ≤ λ 0.5 α if |α| > λ 0.5 .
Hence, replacing a by α α α i , the solution can be reached by a convergent sequence produced by using the IHT algorithm

α α α i +1 = T λ 0.5 a i + D T (x -Da i ) .
It is noticed that, in this algorithm, the descent step size is set to be 1 which could cause unstability. A sufficient condition for stability and convergence is having the eigenvalues of I -D T D strictly between 0 and 1 [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF][START_REF] Blumensath | Iterative hard thresholding and l0 regularisation[END_REF].

When p = 1, the surrogate function becomes differentiable only if α j = 0. Then, except for α j = 0, letting the derivative of the function be equal to 0, the solution is expressed as

α j = TS λ/2 (a j + (D T (x -Da)) j ),
where TS λ/2 is the soft thresholding operator with the definition of

TS λ/2 (α) =        α -λ/2 if α ≥ λ/2 0 if |α| < λ/2 α + λ/2 if α ≤ -λ/2.
By replacing a j by α i j , which is the solution obtained in the i t h iteration, and α j by the updated solution α i +1 j after i + 1 iterations, the Iterative Soft Thresholding algorithm (IST) is obtained with

α i +1 j = TS λ/2 (α i j + (D T (x -Dα i )) j ). (3.19)
The proof of convergence can be found in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

By soft thresholding, the non-zero coefficients are shrunken, namely there is a bias in the coefficient estimation compared to the least-squares solution. Accordingly, relaxed LASSO estimator was proposed in [START_REF] Meinshausen | Relaxed lasso[END_REF], which tunes λ by another parameter varying from 1 to 0. Fujiwara et al. [START_REF] Fujiwara | Fast lasso algorithm via selective coordinate descent[END_REF] improved the efficency of the IST algorithm by a predetermined step in which the necessary updated coefficients are selected and the unnecessary ones are pruned. Moreover, Two-Step Iterative Shrinkage/Thresholding Algorithm (TwIST) updates α α α i +1 depending both on α α α i and α α α i -1 , which shows good performance on denoising and stability [START_REF] José | A new twist: Twostep iterative shrinkage/thresholding algorithms for image restoration[END_REF].

Gradient descent algorithm

Due to the differentiability of the LASSO problem in all α i , except where α i = 0, the gradient descent can be applied to the problem. However, as aforementioned, the step size and its direction are difficult to determine in order to keep the updated coefficient in the acceptance constraint region. In this section, we present the gradient LASSO proposed by Kim et al. in [START_REF] Kim | A gradient-based optimization algorithm for lasso[END_REF], which is stable and convergent.

Similar to the iterative thresholding algorithm, the gradient descent LASSO is composed of two steps: addition step and deletion step. The optimization problem is transformed, by using the variable change ω ω ω = α α α/λ, into arg min

ω ω ω∈S 1 2 x -λDω ω ω 2 2 ,
where

S = {v ∈ R m | v 1 ≤ 1}.
The idea of the gradient LASSO is to find a path of descent ω ω ω to cause reduction of the square error and guarantee ω ω ω is always in the set S, which is realized in [START_REF] Kim | A gradient-based optimization algorithm for lasso[END_REF] with two steps: the first step decreases the square error and the second step corrects ω ω ω in the set S. The outline of the gradient descent LASSO is presented in Algorithm 6. Next, we detail the addition and deletion steps.

The addition step takes advantage of the coordinate gradient descent algorithm (CGD), namely at each iteration, an entry is updated. For instance, Kim et al. select the one with the largest gradient component

ĵ = max j |∇ ω j f (ω ω ω)|,
where

f ω ω ω = 1 2 x -λDω ω ω 2 2 .
Then ω ĵ is updated along the direction in S,

ω i +1 ĵ = ω i ĵ + γ(v -ω i ĵ ),
f (1 -γ)ω ω ω + γv .
The deletion step handles the problem by considering two possible cases: 1) the updated ω ω ω i +1 does not achieve the optimal sparsity level, namely there is non-zero entries in ω ω ω i +1 that ought to be zero in an optimal solution, 2) ω ω ω i +1 has the desirable sparsity level but not optimal. The first case occurs when

〈∇ ω ω ω i +1 I f ,θ θ θ I 〉 < 0
and ω ω ω i +1 1 = 1, where I denotes the set of indices of the active entries I = { j | ω j = 0}, where θ θ θ is a vector that contains all signs of updated coefficient vector θ θ θ = sign(ω ω ω i +1 ). The notation of ω ω ω i +1 I , as a component of ω ω ω i +1 , is given in Chapter 1. Then the direction v was proved to be

v I = -∇ ω ω ω i +1 I f + 〈θ θ θ I ,∇ ω ω ω i +1 I f 〉 ω ω ω i +1 0 θ θ θ I .
Algorithm 6 Gradient LASSO Input: the signal x and the dictionary D, parameter λ, stopping criterion . Output: The sparse representation α α α.

1: function SPARSE CODING 2:

Initialization i = 1, the coefficient vector α α α 0 = 0 3:

while ω ω ω i -ω ω ω i -1 2 2 > do 4:

Addition step

5:

Iteratively updating ω i by gradient descent algorithm

ω m j = (1 -γ)ω i j + vγ, where v = -sign(∇ ω j x -λDω ω ω 2 2 ) and γ = arg min γ∈[0,1] f (ω ω ω i +1 + γv), where v ∈
R m with the j t h entry v j = v.

6:

Deletion step

7:

Finding the index set I which contains the indices of non-zero elements in ω ω ω.

8:

Computing derivative of function f at the sub-vector ω ω ω I , namely Setting v the descent direction with all entries indexed by I the vector v I .

∇ ω ω ω i +1 I f . 9: if 〈∇ ω ω ω i +1 I f ,θ θ θ I 〉 < 0 or ω ω ω i +1 1 = 1 then 10: v I = -∇ ω ω ω i +1 I f + 〈θ θ θ I ,∇ ω ω ω i +1 I f 〉 ω ω ω i +1 0 θ θ θ I ; 11:
15:

Updating ω ω ω i +1 = ω ω ω i +1 + γv, where γ = arg min γ∈J f (ω ω ω i +1 + γv) 16:

Increment i = i + 1 17:
end while 18: end function

The second case is handled in the classical way of gradient descent, that is, v I = -∇ ω ω ω i +1 I f . Then formulate v by assigning all entries indicated by I the corresponding value and the others to zeros. The step size is obtained by minimizing the objective

function but is limited to the interval J = [0, min{-ω i +1 j /v j | ω i +1 j v j < 0, j ∈ I}].
In this way, some coordinates of ω ω ω i +1 would decrease to 0, namely some entries are deleted.

With a finite number of iterations, the gradient based LASSO is ensured to be convergent and with a bounded error [START_REF] Kim | A gradient-based optimization algorithm for lasso[END_REF]. Furthermore, this algorithm has no RIC assumption but only requires the Lipschitz continuity of ∇ ω ω ω f .

Algorithms for dictionary updating

The subproblem of dictionary updating, aims at finding the local optimal solution of the problem (3.4) while fixing the matrix of decomposition coefficients A of for-

mulation min D∈D 1 i =1 1 2 x i -Dα α α i 2 2 + λΩ(α α α i ).
(3.20)

By removing all the constant terms that are associated with the matrix of decomposition coefficients A, this subproblem can be reduced to the following optimization problem:

min D∈D 1 2 X -DA 2 F , (3.21)
where • F denotes the Frobenius norm. This problem is differentiable and convex, which is thus much easier to be solved. In this section, we present two algorithms for dictionary updating: the Method of Optimal Directions (MOD) [START_REF] Kjersti Engan | Method of optimal directions for frame design[END_REF] and the method of Coordinate Descent based on SVD [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF].

The method of optimal directions (MOD)

The MOD algorithm proposed by Engan et al [START_REF] Kjersti Engan | Method of optimal directions for frame design[END_REF] updates the dictionary D by iteratively moving along the optimal direction that minimizes the representation residual.

Supposing that the sparse representation A of the signal matrix X is obtained, its representation residual can be expressed by

r (D) = i r i (D) = i x i -Dα α α i 2 2 = X -DA 2 F . (3.22)
The objective is to update the dictionary D to D + ∆, where the increment matrix ∆ is obtained by minimizing this residual, namely arg min

∆∈R n×m r (D + ∆). (3.23)
Moreover, r (D + ∆) < r (D) would be satisfied to guarantee convergence of the algorithm. By computing the gradient of the above cost function and setting it to zero, we achieve the result

∆AA T = (X -DA)A T .
Thus, D can be updated by D + ∆,

namely D i -1 + (X -D i -1 A)A T (AA T ) -1 , which can be simplified to D i = XA T (AA T ) -1 [SE12].
The MOD algorithm is usually used for dictionary updating in sparse model researches [MLB + 08], Moreover, its variant, Method of Optimal COherence-COnstrained Direction (MOCOD) is developed for incoherent dictionary learning [START_REF] Ramírez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF].

The algorithm of coordinate descent based on SVD

While the MOD algorithm updates the whole dictionary at each iteration, the Co- thus the matrix multiplication DA is represented by the sum of outer products j d j b T j . When updating d i with fixing all the other atoms, for keeping the coefficient matrix in sparsity constraint, only the signals with j t h entry non-zero are considered. Let the set of these indices be denoted by J, thus the set of the signals can be denoted by X J . Supposing E = X Jj =i d j b T j ,J , the problem with respect to

d i is now di = arg min d i 2 =1 E -d i b T i ,J 2 F .
This problem aims at finding a normalized vector di over which E achieves the largest projection. The SVD algorithm is able to find the closest rank-1 matrix (i.e., vector) that approximates E with great extent. Thus, d i and b i ,J will be simultaneously updated in this phase, this is unique in dictionary learning algorithm.

Algorithms for dictionary learning

In this section, we restrict the presentation to the classical dictionary learning, We will not present all the algorithms of dictionary learning for each special task but only the classical one, that is, the signal denoising or signal reconstruction problem. We rewrite the optimization problem of dictionary learning (3.4), arg min

A∈R m×l ,D∈D 1 2 X -DA 2 F + λ i Ω(α α α i ), (3.25)
where the regularization Ω(α α α i ) can be α α α i 0 or α α α i 1 , or any other sparse promoting penalty. Like the sparse coding problem, the above dictionary learning problem has also its constrained formulations, which are the sparsity constrained problem arg min

A∈R m×l ,D∈D 1 2 X -DA 2 F , α α α i ∈ S k ∀i = 1, . . . , l , (3.26)
and the error constrained problem arg min

A∈R m×l ,D∈D i Ω(α α α i ) subject to X -DA 2 F < . (3.27)
The latter optimization problem, jointly on A and D, is NP-hard. It is difficult to solve this problem by directly using optimization techniques. Thus, we takes strategies of tackling the problem by two alternating steps: sparse coding, which is a problem with respect to A by fixing D, and dictionary updating, which fixes A and aims at computing D. The combination of sparse coding algorithms (in Section 3.2) and dictionary updating algorithm (in Section 3.3) allows to define dictionary learning algorithms. This section introduces the most frequently used algorithms for i = 1 to N do 4:

Sparse Coding

5:

Computing sparse representation α α α i j for each signal x j with j varying from 1 to l by applying OMP algorithm. Finding all signals with j t h entry non-zero which is indexed by J, denoted by X J 10:

Calculating E = X J -p = j d p (b i ) T p,J .
11:

Doing SVD decomposition of E = USV T , then

d i j = u 1 and (b i ) T j ,J = s(1, 1)v 1 .
12: end for 13: end for 14: end function for dictionary learning.

K-SVD

The K-SVD algorithm is a representative dictionary learning algorithm that is proposed as a generalization of the k-means algorithm, but regarding a signal as a linear combination of several atoms [AEB + 06].

Since it was proposed, efforts have been conducted on the improvement of implementation [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF] and some theoretical foundations [START_REF] Schnass | On the identifiability of overcomplete dictionaries via the minimisation principle underlying K-SVD[END_REF], which proved the certainty of reaching a local minimum under the condition of coefficient sufficient decay. Hence, the algorithm is widely applied with success in image denoising [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], face recognition [START_REF] Zhang | Discriminative K-SVD for dictionary learning in face recognition[END_REF] and classification [START_REF] Jiang | Label consistent K-SVD: Learning a discriminative dictionary for recognition[END_REF]. K-SVD deals with the optimization problem (3.26) by iteratively alternating the two phases, sparse coding and dictionary updating as follows: it exploits OMP for sparse coding, and dictionary update is operated using the Coordinate Descent Method. While the MOD method updates the whole dictionary in each step, the K-SVD updates each atom of the dictionary successively, through the singular value decomposition (SVD) [AEB + 06], as introduced in Section 3.3.2.

The complete K-SVD algorithm is shown in Algorithm 7. In practice, for speed consideration, the atoms of the dictionary are not updated in order but obeying a

Sum of outer products dictionary learning

The Sum of OUter Products Dictionary Learning (SOUPDIL) algorithm, recently proposed by Ravishankar et al [START_REF] Ravishankar | Efficient sum of outer products dictionary learning (soup-dil) and its application to inverse problems[END_REF], is inspired by the K-SVD algorithm. Specifically, the sparsity regularized dictionary learning problem (3.25) is rewritten as an optimization problem with respect to the matrix B defined in the K-SVD algorithm, rather than the original sparse representation A, that is, arg min Initialization the dictionary D 0 and the matrix B 0 = 0 3:

B∈R l ×m ,D∈D 1 2 X -DB T 2 F + λ B 0 , (3.28
for t = 1 to N do 4:
for i = 1 to m do 5:

Sparse Coding

6:

Calculating E i E i = X - j <i d t j b t j T - j >i d t -1 j b t -1 j T ; 7: Updating b t i by b i = min(|T s (E i T d t -1 i )|, L) sign(T s (E t i T d t -1 i )); 8:
Dictionary updating 9:

Computing E i b t i .

10:

Updating atom d t i by

d i = E i b i E i b i 2 if b i = 0 v else; 11:
end for 12:

end for 13: end function non-coercive objective. This leads to the resulting optimization problem

arg min b∈B E i -d i b T i 2 F + λ b i 0 .
(3.30)

This problem can be handled by using the IHT algorithm with the threshold set to s = λ. The closed-form solution is thus

b i = min(|T s (E T i d i )|, L) sign(T s (E T i d i )),
where indicates the element-wise multiplication.

For updating an atom d i while fixing all the other variables, the problem be-comes arg min

d i 2 =1 E i -d i b T i 2 F . (3.31)
The solution is obtained by

d i = E i b i E i b i 2 if b i = 0 v else,
where v is a unit-norm vector. In [START_REF] Ravishankar | Efficient sum of outer products dictionary learning (soup-dil) and its application to inverse problems[END_REF], the first column of the identity matrix is used.

The SOUPDIL algorithm takes advantage of the manner of updating alternatively d i and b i , which is proved to be convergent in a limited number of iterations.

We summarize SOUPDIL in Algorithm 8.

Proximal method

The proximal method [PB14] is a tool for dealing with non-smooth, large-scale problems. Furthermore, Bolte et al proved that the proximal method can reach convergent result in non-convex problems [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. Due to the advantages of the proximal method, it is used in image processing and machine learning [BJQS14, JMOB10, LDL15]. Indeed, the proximal method is also widely used in dictionary learning [JMOB10, CPR13], and moreover, it is proved to be globally convergent [START_REF] Bao | L0 norm based dictionary learning by proximal methods with global convergence[END_REF].

The dictionary learning problem can be rewritten in a generalized formulation of non-convex and non-smooth problem with respect to the two variables A and D arg min

A,D f (A) + q(A, D) + g (D), (3.32)
where f = A 0 (or f = A 0 + δ A (A) ) and g (D) = δ D (D) are proper lower semicontinuous functions, and q(A, D) = 1 2 X -DA 2 F is a differentiable function with a Lipschitz continuous with determined Lipschitz constants. In these expressions, the indicator functions are defined on the sets

A = {M ∈ R m×l | M 1,∞ ≤ c} and D = {D ∈ R n×m | d i 2 = 1, ∀i = 1, . . . , m}.
The proximal method solves the optimization problem (3.32) by iteratively updating A and D. Specifically, the linearized proximal minimization problem with respect to A, when fixing D, is

 = arg min A f (A) + 〈A -A , ∇ A q(A , D )〉 + t 1 2 A -A 2 F , (3.33)
where ∇ A q(A , D ) denotes the partial derivative with respect to q of A at point (A , D ) and t 1 an appropriate step size. Similarly, the optimization problem with respect to D is

D = arg min D g (D) + 〈D -D , ∇ D q(A , D )〉 + t 2 2 D -D 2 F , (3.34)
where ∇ D q(A , D ) is the partial derivative of q with respect to D at the point (A , D )

and t 2 the corresponds to the appropriately chosen step size. Thus, by replacing Â, A , A , D, D by A t +1 , A t , A t +1 , D t +1 , D t respectively, processing iteratively the problems (3.33) and (3.34) produces a convergent sequence {(A t , D t )} (see proof in [START_REF] Bao | L0 norm based dictionary learning by proximal methods with global convergence[END_REF]).

In fact the problem (3.33) can be solved by the IHT method and the solution of problem (3.34) can be reached by optimization methods such as the Lagrangian method [START_REF] Bertsekas | Nonlinear programming[END_REF]. Specifically, Bao et al used a block coordinate descent algorithm to get the optimal D t +1 . Here, we provide directly the closed-form solutions of the optimization problems without the detail of calculation

   A t +1 = min T 2λ/t 1 (A t -1 t 1 ∇ A q(A t , D t )), c d t +1 j = d m j / d m j 2 , ∀ j = 1, . . . , m,
where d m j is calculated by a gradient descent, that is,

d m j = d t j -1 t 2 ∇ d j q(A t +1 , D t j ),
where D t j is the dictionary in the (t + 1) t h iteration where all its atoms for i ≤ j are updated atoms d t +1 i and the other unchanged

d t i , namely D t j = [d t +1 1 , . . . , d t +1 j -1 , d t j +1 , . . . , d t m ]
. Besides, the parameters determination is quite important in implementation of proximal method. In general, one should ensure that t 1 and t 2 are not larger than the Lipschitz constant of the gradient function with respect to A and d j , respectively [START_REF] Parikh | Proximal algorithms[END_REF]. Hence, for proximal method, it is the direction and descent step, which jointly ensure the reduction in the objective function. This is not the same with Algorithm 9 Proximal Method for 0 -norm based Dictionary Learning Input: the input signals X, parameter λ, the step size parameter t 1 and t 2 , the maximal number of outer iterations N. Output: The sparse representation A and the dictionary D.

1: function DICTIONARY LEARNING 2:

Initialization the dictionary D 0 and the matrix A 0 = 0 3:

for t = 1 to N do 4:
Sparse Coding 5:

Calculating A t A t +1 = min T 2λ/t 1 (A t -1 t 1 ∇ A q(A t , D t )), c ; 6:

Dictionary updating

7:

Updating iteratively atom d t i by

d m j = d t j -1 t 2 ∇ d j q(A t +1 , D t j ); d t +1 j = d m j / d m j 2 ; 8:
end for 9: end function greedy and IHT algorithms, since they require a strict condition on the constant of the Restricted Isometry Condition (RIC).

Conclusion

This chapter presented the methods to solve the optimization problems of sparse coding and dictionary learning for sparse representation. Specifically, the dictionary learning problem, estimating jointly the sparse code and the dictionary, is NPhard and intractable. The frequently exploited strategy is to iteratively process two alternating phases: sparse coding and dictionary updating. This optimization problem can thus be transformed into two tractable subproblems. This chapter presented the state-of-the-art algorithms of sparse coding, of dictionary updating and also of dictionary learning. In Section 3.2, the sparse coding algorithms were discussed within three categories: 1) matching pursuit algorithms for dealing with 0 based problems, 2) relaxation methods by replacing 0 with p , specially the case p = 1 addressed with Basis Pursuit and LASSO, and 3) thresholding algorithms that, in fact, cover both the first problem of 0 penalty or regularization, and the second one, namely LASSO problem and Basis Pursuit. In Section 3.3, the main dictionary updating algorithms were introduced, with the representative algorithm MOD and the Coordinate Descent algorithm based on SVD. In the last section, we presented the main algorithms of dictionary learning that combine the appropriate sparse coding algorithm and dictionary updating algorithm. The well-known algorithms are K-SVD, SOUPDIL, and proximal method.

In the next chapter, we will propose a novel strategy to address the exact 0 -norm optimization problem for dictionary learning. The problem of dictionary learning plays an important role in sparse representation. This problem with quadratic objective function involving two optimisation variables, the sparse code A and the dictionary D, is intractable. To address this problem, one takes usually the strategy of iteratively alternating the two steps of sparse coding and dictionary updating. The subproblem corresponding to the dictionary updating is convex, which can be easily solved by algorithms such as MOD and gradient descent with SVD, as described in Section 3.3. Nevertheless, the subproblem of the sparse coding with the 0 -norm constraint, due to the characteristics of the 0 -norm function, is non-smooth, non-convex, non-differentiable and thus NP-hard [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF]. Generally, researchers take strategies like greedy algorithms [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], relaxation optimization techniques (LASSO [START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF] and concave function p with 0 < p < 1 [START_REF] Rinaldi | Concave programming for minimizing the zero-norm over polyhedral sets[END_REF]) or iteration shrinking/thresholding algorithms [START_REF] Yang | Ways to sparse representation: An overview[END_REF]. However, all these methods risk of missing the global optimal solution [BJQS14, ZXY + 15].

In this chapter, we address the sparse coding with its original 0 -norm formulation. To this end, we reformulate the sparse representation problem as a mixed integer programming (MIP) problem [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF]. Thereby, the optimization problem with 0 constraint is transformed into an optimization problem that can be solved with off-the-shelf optimization software packages, such as CPLEX and Gurobi Optimizer. Moreover, we propose two acceleration techniques that allow to take on real datasets. Finally, we demonstrate the relevance of the proposed MIP based dictionary learning in well-known image denoising tasks. In this section, we will explore a novel approach to address the sparse coding problem with the 0 -norm constraint. By reformulating the 0 constrained problem, it is transformed into a problem of MIP, which makes possible to exactly solve the 0 based problem.

Optimized algorithm for exact

Mixed integer programming (MIP)

Mixed Integer Programming (MIP) refers to optimization problems where some variables are restricted to be integer while others not. problem. In the following, we firstly give an introduction of MIP. Then, we reformulate the problem of sparse coding as MIP. Finally, some advanced optimization methods for dealing with MIP problems will be presented.

Introduction of MIP

The history of MIP dates back to the 1940s [START_REF] Robert | A brief history of linear and mixed-integer programming computation[END_REF]. Its development relies on the algorithmic improvement on Linear Programming (LP) [START_REF] Solow | Linear and nonlinear programming[END_REF][START_REF] Nocedal | Numerical optimization[END_REF], like the simplex method [START_REF] John | A simplex method for function minimization[END_REF], barrier methods [GMS + 86], and the computational capability of hardware. According to [START_REF] Robert | A brief history of linear and mixed-integer programming computation[END_REF], from 1988 to 2004, during these sixteen years, the improvement factor exceeds six orders of magnitude. This significant progress in LP solvers has been injected in solving MIP problems. This is possible thanks to techniques like cutting plans, which bridges the gap between MIP and LP [START_REF] Sen | Decomposition with branch-andcut approaches for two-stage stochastic mixed-integer programming[END_REF].

Hence, MIP is becoming tractable in real-world problems [START_REF] Schouwenaars | Mixed integer programming for multi-vehicle path planning[END_REF]. However, its application is still limited to small size problems [START_REF] Tawarmalani | Global optimization of mixed-integer nonlinear programs: A theoretical and computational study[END_REF][START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF].

The standard formulation of MIP problem with respect to a variable v ∈ R n is arg min Specially, if Q is semi-definite positive, the objective function in (4.2) becomes convex. It is noticed that if Q i = 0, the problem turns into a MIQP. Moreover, if Q is also equal to 0, the problem becomes a MILP.

v v T Qv + c T v (4.2) subject to A i n v ≤ b i n (4.3) lb ≤ v ≤ ub (4.4) v T Q i v + c T i v ≤ b i (4.5) v I ∈ Z,

Reformulation of sparse coding in MIP problem

The first attempt to tackle sparse coding as a MIP problem goes to Jokar and Pfetsch [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF]. They proposed to use the method of branch-and-cut for getting the exact solution, which was regarded as the criterion to evaluate the performance of the heuristic methods mentioned in [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF] such as the BP and the OMP. Even though this method is of high computational complexity, which limits its application to small-scale instances. Nevertheless, it is worthy noting that the exact solution can be used to verify the theoretical conclusion, namely when the optimal solution is sparse enough, it can be found by the heuristics, e.g. greedy algorithms. However, for the cases where the optimal solution has more non-zeros, the heuristics perform usually bad. Quite recently, Bourguignon et al reformulated in [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF] the 0 based sparse coding problem as MIP problem, by replacing the logic relation with a big-M reformulation. Thus, MIP solver was applied in decomposition of small scale synthetic data with added noise [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF]. In our work, MIP solver is used for sparse coding, which is further applied in dictionary learning [START_REF] Liu | Mixed integer programming for sparse coding: Application to image denoising[END_REF]. Furthermore, this model is not limited to small optimization problems but applied in image processing thanks to two proposed acceleration techniques.

In this section, we focus on the phase of sparse coding of problem (4.1) which is the estimation of optimal A with fixing D. We recall that, by assuming independence of signals in X, the sparse coding problem can be regarded as l subproblems with respect to α α α i respectively. For expression simplification, hereafter, x i and α α α i are expressed by x and α α α respectively. Hence the problem is now arg min

α α α∈R m 1 2 x -Dα α α 2 2 , subject to α α α 0 ≤ k. (4.7)
This constrained optimization problem can thus be reformulated, with all the entries of the sparse vector α α α indicated by a binary variable z ∈ {0, 1} m , which can be explained by the logical relation:

α i = 0 if z i = 0 α i = 0 if z i = 1 , (4.8)
where z i and α i indicate the i -th entries of the vectors z and α α α respectively, i = 1, . . . , m.

Since such logical relation cannot be easily integrated into the objective function, we recast the sparsity condition into a linear inequality by introducing a sufficiently large value M > 0 ensuring that α α α ∞ < M for any desirable solution α α α, where

• ∞ means the maximal magnitude of the vector's entries. A too large M value will result in an increased feasible region, which will make the problem less computational efficient. An appropriate value of M improves the performance. A method providing a lower value for M to obtain tight bounds is crucial in algorithm improvement.

Now the indicative function of z is ensured by satisfying the constraints:

-

z i M ≤ α i ≤ z i M, . . .∀i ∈ {1, . . . , m}. (4.9)
Then, the sparsity constraint α α α 0 ≤ T can be depicted by z as:

p i =1 z i ≤ k. (4.10)
As a consequence, the 0 -based sparse coding problem can have a 'big-M' reformu-lation, that is, for a given M large enough:

min α α α∈R m ,z∈{0,1} m 1 2 x -Dα α α 2 2 subject to -zM ≤ α α α ≤ zM 1 T m z ≤ k, (4.11)
where 1 m is the column vector of size m with all elements equal to one. In this formulation, the optimization variables α α α and z are respectively continuous and integer vectors. Hence, the sparse coding problem is reformulated as a mixed-integer programming (MIP) problem.

The above reformulation (4.11) of the sparsity constrained problem (4.7) is valid, as proven by the following proposition.

Proposition 1. The reformulated problem (4.11) is equivalent to its original problem (4.7).

Proof. For proving the equivalence of the two problems, (4.11) and (4.7), we need to prove both the equivalence of the objective functions and of the feasible regions of α α α defined by their constraints.

It is obviously that the two objective functions are identical. Thus, it is sufficient to prove the equivalence of the two feasible regions that are defined respectively by

A 1 = α α α ∈ R m | α α α 0 ≤ k ,
and

A 2 = α α α ∈ R m | z ∈ [0, 1] m , -zM ≤ α α α ≤ zM, 1 T m z ≤ k .
First, we prove that A 1 ⊆ A 2 . Supposing arbitrary v ∈ A 1 , thus, we have v 0 ≤ k. By considering the definition of z, and there exists M sufficiently large and M ≥ v ∞ , then it yields -zM ≤ v ≤ zM and 1 T m z ≤ k. Hence, v ∈ A 2 . Then, we prove the reciprocal, namely A 2 ⊆ A 1 . Let v be an arbitrary vector in A 2 , then it satisfies the condition -zM ≤ v ≤ zM with z ∈ [0, 1] m . It is easy to imply that, if z i = 0, the corresponding entry v i will be equal to 0. By considering the condition 1 T m z ≤ k, we have that there are no more than k non-zero entries in v. That is to say

v 0 ≤ k, namely v ∈ A 1 .
Finally, by combining A 1 ⊆ A 2 and A 2 ⊆ A 1 , we conclude that A 1 = A 2 . Moreover, our optimization problem defined in (4.11) with continuous and integer optimization variables, the objective function is quadratic and all the constraints are linear. Consequently, the sparse coding can be interpreted as a MIQP. In the following, we will write this optimization problem in the standard form expressed in (4.2) ∼ (4.6) in order to use off-the-shelf solvers. Specifically, in the standard MIQP reformulation of sparse coding problem, v is obtained by combining the vectors α α α and z, that is, let v = (α α α T , z T ) T , then the sparse coding problem is of form

min v 1 2 v T Qv + c T v subject to A i n v ≤ b i n v j ∈ {0, 1} ∀ j ∈ I,
A i n =     -I m -MI m I m MI m 0 T m 1 T m     ,
where I m denotes the identity matrix of size m × m, and the (2m + 1) × 1 column vector in the right-hand-size of the inequality is b i n = (0 T 2m , k) T . Finally, the set I in (4.12) indicates the integer components in the MIQP, namely

I = {m + 1, m + 2, . . . , 2m}.
During the implementation, the variable's type is specialized as continuous or binary (namely variable's type is a parameter in MIQP construction) in the input for the solver at hand.

As an alternatve to interpreting the logical relation (

α α α T (1 1 1 m -z) = 0,
where we can find that, when z i = 0, the entry α i is forced to be 0 (thus the equality constraint is active), while z i = 1 disables the constraint. Specifically, in this problem, the indicator function z switches the constraint in a "complementary" way. By adopting this interpretation, the sparse coding problem (4.7) is thus reformulated as min

α α α∈R m ,z∈{0,1} m 1 2 x -Dα α α 2 2 subject to α α α T (1 1 1 m -z) = 0 1 T m z ≤ k. (4.13)
This problem is an MIQP problem but with a nonlinear constraint. It is also called disjunctive programming [START_REF] Bonami | On mathematical programming with indicator constraints[END_REF]. Unlike the big-M reformulation of which the difficulty lies on the determination of an appropriate value for M, this reformulation can be solved by exploiting the lift-and-project cutting technique [START_REF] Fischetti | On the separation of disjunctive cuts[END_REF]. For some special problems, for instance the SVM problem [BBF + 16], the MIQP problem can degenerate into its continuous version. However, this is not the case for the sparse coding problem.

This reformulation risks to highly increase the computational complexity. In Section 4.1.2, we provide two techniques to accelerate the resolution of the above MIQP problem. The first one determines a relevant value of M, thus the big-M reformulation can be more efficiently implemented compared to that of disjunctive programming. The second method considers a convex envelop of the feasible continuous variables, which allows to solve the problem more efficiently. But before, we present recent advanced optimization methods for MIP.

Advanced optimization methods for MIP

This section gives an introduction to advanced optimization techniques for solving the MIP problem. By making clear the principle of the resolution of the MIP problem, it will help us in finding the way to accelerate the MIP solver. Thus, it offers some theoretical support of the inspirations and the motivations for the acceleration method we will propose in the following section.

Because of the optimization over both discrete and continuous variables, the MIP problems cannot be solved by common LP solvers, such as simplex method +∞. This problem can be solved by using a conventional LP solver, such as simplex method [START_REF] Juan | A lifted linear programming branch-and-bound algorithm for mixedinteger conic quadratic programs[END_REF] or interior point method [START_REF] Borchers | Using an interior point method in a branch and bound algorithm for integer programming[END_REF].

2) Stopping criterion testing: If the active node set is empty, the obtained solution v with vI ∈ Z is optimal and it yields the incumbent objective value which updates the upper bound ub. Otherwise, there is no feasible solution to the problem.

3) Subproblem selection: By an adequate searching method [LS99], e.g. depth first search strategy [START_REF] Borchers | Using an interior point method in a branch and bound algorithm for integer programming[END_REF], the subproblem is successively selected and solved.

4) Node deletion: When the objective value of the subproblem has the minimal value lb > ub or there is no feasible solution produced, then this node will be deleted.

5) Upper bound updating or branching:

When a feasible solution is found with its objective value l b, then if lb < ub, update ub = lb [START_REF] Gleixner | Learning and propagating Lagrangian variable bounds for mixed-integer nonlinear programming[END_REF]. However, when the optimal solution with v i , i ∈ I is a fractional number, then branch the problem and continue to solve this subproblem.

The branch-and-bound method makes the MIP problem tractable. However, by only applying this method, the MIP problem is difficult to be solved efficiently.

Hence, another optimization technique, cutting plans [MMWW02] is integrated

with the branch-and-bound method, which is called branch-and-cut [START_REF] Mitchell | Branch-and-cut algorithms for combinatorial optimization problems[END_REF][START_REF] Karuppiah | A Lagrangean based branch-and-cut algorithm for global optimization of nonconvex mixed-integer nonlinear programs with decomposable structures[END_REF].

The cutting plans method provides a way to relax the original problem [START_REF] Mitchell | Branch-and-cut algorithms for combinatorial optimization problems[END_REF]. Suppose that the feasible region is a polyhedron defined by all linear inequality constraints. If such a convex hull of feasible region is found, then the problem will turn into an LP problem, which makes it much easier to be solved. The cutting plans offer the method for iteratively generating this convex hull. When the cutting plans method is used in branch-and-bound, it can be applied for all the problems or just a subproblem in certain nodes [START_REF] Karuppiah | A Lagrangean based branch-and-cut algorithm for global optimization of nonconvex mixed-integer nonlinear programs with decomposable structures[END_REF].

Branch and cut algorithm makes a theoretically improvement in MIP solver.

Moreover, with development in hardware, the computational capability has been greatly improved [START_REF] Bliek1ú | Solving mixed-integer quadratic programming problems with IBM-CPLEX: A progress report[END_REF].

It is worthy noting that in MIP solver, the upper bound ub has an important effect on the computational efficiency. A samllest upper bound means fewer subproblems will be selected and solved. Besides, a good relaxation is crucial in approximating the convex hull of the feasible region. These two methods will be discussed in next section where we provide acceleration of MIQP solver for sparse coding.

For implementation of the MIQP solver, various optimization software packages can be explored. The most known are CPLEX and Gurobi Optimizer. IBM develops CPLEX1 which integrates the latest MIP solvers to solve larger MILP problems [START_REF] Bliek1ú | Solving mixed-integer quadratic programming problems with IBM-CPLEX: A progress report[END_REF]. The recently developed Gurobi Optimizer2 can have an equivalent performance to CPLEX while the latest releases get some improvements [START_REF] Robert | A brief history of linear and mixed-integer programming computation[END_REF]. The developed tools make it possible to apply MIQP into image processing, but by considering its computational complexity, some effort should be done to improve it, as described next.

Accelerated mixed integer programming (AcMIQP)

MIQP problem, as aforementioned, can be solved via bound and cut algorithm.

However, the problem of computational cost has been raised, especially when dealing with real data, such as real images, even though the most advanced devices show powerful ability in computing. In this section, we propose a method of initialization and a relaxation technique for speeding up the MIQP based sparse coding algorithm.

Initialization by proximal method

A good initialization offers a warm start of the tree search, which will help to reduce the size of the problem. Furthermore, with good initial solution, we can get a tighter upper bound by calculating its corresponding objective value and an adequate value M which defines the scale of feasible region (a lower value of M means a smaller search region). Thus, the initialization method acts importantly in the acceleration of the MIQP solver.

In the following, the proximal method will be used for generating an initial solution [PB14, BJQS14], which is theoretically explicated by solving the problem, arg min

α α α 1 2 x -Dα α α 2 2 + δ S k (α α α), (4.14)
where δ S (α α α) denotes the indicator function on the set S of k-sparse vectors, namely

S k = v ∈ R m | v 0 ≤ k .
The proximal method can be regarded intrinsically as finding the minimal upper bound of the objective function by first-order approximation, which is based on the fact that there exists a real value ρ that satisfies

1 2 x -Dα α α 2 2 ≤ 1 2 x -Dα α α i 2 2 + 〈α α α -α α α i , D T (Dα α α i -x)〉 + 1 ρ α α α -α α α i 2 2 . (4.15)
Hence, the solution to the problem (4.14) will be reached by a convergent sequence created by

α α α i +1 = arg min α α α 〈α α α -α α α i , D T (Dα α α i -x)〉 + 1 ρ α α α -α α α i 2 2 + δ S k (α α α), (4.16)
We solve this problem by introducing the proximal operator [START_REF] Parikh | Proximal algorithms[END_REF] defined by

prox t h (u) = arg min x h(x) + t 2 x -u 2 ,
where h defines a proper and lower semi-continuous function, and t > 0 is a step size parameter. See [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] for more detail. In problem (4.16), h(x) is the indicator function δ S k to make sure that the feasible region is in the space S k of k-sparsity.

The proximal operator boils down to the projection onto the sparse space S k :

P S k (u) = arg min x∈S k ( x -u 2 ).
The solution of this problem can be easily obtained by keeping the k largest absolute value components of u and setting the other components to zero:

P S k (u) = u j if j ∈ {(
u = α α α i + ρD T (x -Dα α α i ).
Thus, by applying a proximal algorithm, the sparse representation problem can be solved through the iterative update process:

α α α k+1 ∈ P S k α α α k + ρD T (x -Dα α α i ) .
(4.17)

After a finite number of iterations n i t er , the α α α n i t er will be much closer to the optimal solution of the exact 0 problem.

By considering the definition of 'big-M', the constraints in the problem (4.11) related to M can be well determined by an approximation of the optimal solution. A simple method to determine an appropriate value for M can be:

M = η x n i t er ∞ . (4.18)
The much tighter bound defined by M and an approximate initialization allow to speed up the algorithm. Its performance will be discussed in the following section.

It is noticed that the proximal method for 0 constrained sparse coding problem has exactly the same operation in the IHT algorithm (Algorithm 5) introduced in Section 3.2, that is, the solution is reached by iteratively processing two alternating steps, gradient descent and projection. However, the two algorithms are built on different theoretically foundations. The IHT [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF] assumes that the dictionary D satisfies the RIC with the isometric constant δ 2k < 1/3. And the step size ρ depends on δ 2k . In contrast, the proximal method has no assumption on D, except that the gradient of the function 1 2 x -Dα α α 2 2 with respect to α α α is Lipchitz continuous with constant L. The step size in the proximal method cannot be larger than this constant.

Relaxation

The developments of the MIQP solvers have been following the progress in LP theory. The advanced-start capabilities of simplex algorithms in the branch-andbound [START_REF] Oskar Von | Decomposition of mixed-integer optimal control problems using branch and bound and sparse direct collocation[END_REF] (or now more correctly, branch-and-cut [START_REF] Neumaier | Safe bounds in linear and mixed-integer linear programming[END_REF]) search tree are well exploited by MIQP solvers. No matter which optimization technique is used, the search process remains the main time consumption factor. The searching time heavily relies on the feasible region determined by the constraints. Hence, the effort spent on getting a good formulation of the constraints do help to accelerate the resolution of the optimization problem.

Hoffman and Ralphs have proven in [START_REF] Karla | Integer and combinatorial optimization[END_REF] that, if a feasible solution is obtained by a relaxation, then it must also be the optimal solution to the original problem.

Especially, in the ideal case, if the convex envelope is found, a mixed integer programming will be transformed to the classical linear programming problem. However, it is an NP-hard problem to find the constraints defining the convex envelop [START_REF] Amir | NP-hardness of deciding convexity of quartic polynomials and related problems[END_REF]. The viable strategy is to create a convex envelop of the continuous variables

C = α α α ∈ R m z ∈ {0, 1} m , m j =1 z j ≤ k, |α α α j | ≤ z j k, ,
by adding the constraint about 1 -norm and ∞ -norm of α α α:

     m i =1 |α i | < kM |α i | < M ∀i = 1, . . . , m.
The absolute values can be formulated as linear programs. To do so, we replace each unrestricted variable α i , for i = 1, . . . , m, with the difference of two restricted variables,

α i = α + i -α - i α + i , α - i 0, namely in matrix form α α α = α α α + -α α α - α α α + ,α α α -0.
Then, the absolute value of α i in the above constraints can be represented in the linear program as:

|α i | = α + i + α - i ∀i = 1, . . . , m.
Thus, the constraints for MIQP can be summarized as:

             m i =1 α + i + α - i < kM -zM < α α α + -α α α -< zM 0 ≤ α + i , α - i < Mz 1 T p z ≤ kz . (4.19)
With the new constraints, MIQP can be reformulated as the standard formulation by introducing as updated optimization variable v = (α α α + T ,α α α -T , z T ) T . Accordingly, the model components Q, c, A i n , b i n , l and u are updated respectively as following: The matrix Q becomes the 3m × 3m matrix

Q =     D T D -D T D 0 m,m -D T D D T D 0 m,m 0 m,m 0 m,m 0 m,m .     .
By this representation, we know Q is symmetric and semi-definite positive. It ensures the convexity of the objective function. The vector c becomes the vector of size 3m

c =     -D T x D T x 0 m,1 .   
and the right-hand-side of the inequality constraint becomes

b i n =     kM k 0 T 2m .     .
The two bounds of the new variables v are now defined respectively as u = (M1 T 2m , 1 T ) T and l = 0 3m , and

I = {2m + 1, 2m + 2, . . . , 3m}.
With the new formulation, the problem can be solved more efficiently.

Performance of the AcMIQP

In this section, we will give some experimental evaluations of the performance of the proposed algorithm AcMIQP. The experiments are carried out on synthetic data in order to assess the AcMIQP sparse coding algorithm. The classical MIQP solver, proximal method and OMP are chosen for comparison. By analysis on the accuracy of sparse representation and the computational complexity, the advantage of the proposed AcMIQP becomes remarkable, which provides the support to use AcMIQP in the more complex problem of dictionary learning and further apply it in image processing.

For evaluating the AcMIQP solver in solving the sparse coding problem (4.7), a sparse matrix A ∈ R 128×10000 is created with a column-wise maximum sparsity level of 6, and the dictionary D ∈ R 64×128 , which is column-wise normalized matrix randomly generated from a Gaussian distribution. With the available matrices A and D, the training data X is finally produced by the following equation

X = DA + κE, (4.20)
where E is a randomly generated zero-mean white Gaussian noise matrix and κ a parameter controlling the noise level, set to κ = 0.01 in the experiments.

For statistical purpose, data of size 10000 is divided into 100 units. For each unit X i ∈ R 64×100 , a sparse code matrix A i is estimated. This allows to provide the median, the 5th and 95th percentiles.

The two MIQP sparse coding algorithms, with and without acceleration, are compared to the OMP and proximal methods. The performances of the sparse coding methods are evaluated with three criteria: the difference between A i and the estimated Âi , i.e., A i -Âi F , the reconstruction error, i.e., X i -D Âi 2 F , and the percentage of zero and non-zero elements of the sparse code being found in the right positions.

Table 4.1 presents the computational time and results, in terms of the reconstruction error, the accuracy of the sparse coding estimation and the percentage of number of zero and non-zero elements being recovered in the right position. It shows that the errors obtained by MIQP are far less than that the ones of OMP and proximal methods. Furthermore, the introduction of optimization techniques (including initialization and relaxation) for acceleration has a little effect on the accuracy, while the computational cost is reduced by a factor of 5. In addition, more nonzero elements are found in the right positions. These advantages make AcMIQP of great interest to be used as a sparse coding algorithm and, in conjunction with a dictionary updating rule, as a dictionary learning algorithm. In spite of the over-all strength of MIQP, its Achilles' heel is the excessive computational complexity, making it difficult to use for large-scale problems. However, as aforementioned, the proposed acceleration opens the possibility to apply the MIQP-based dictionary learning algorithm on large-scale problems, such as in image denoising. Note that, in practice, to get an improvement over the proximal method, there is no need to run the optimization until the global minimum. Whatever the computing budget is allocated, the AcMIQP formulation allows to use it to improve the results.

Dictionary learning with AcMIQP for sparse coding

In the previous section, the proposed AcMIQP algorithm was used with success in solving sparse coding problem. The experiments on synthetic data showed the feasibility of the AcMIQP algorithm in processing data with high dimension. In this section, AcMIQP will further be applied in a more sophisticated task: dictionary learning, namely the problem (4.1) given at the beginning of Section 3.4. To this end, the problem is addressed by iteratively processing two alternating steps: sparse coding and dictionary updating.

In the sparse coding phase, the problem (4.7) is reformulated as an MIQP problem. By exploiting the AcMIQP algorithm, the sparse coefficient matrix can be estimated at each iteration. Let A q be the estimated coefficient matrix at the q t h iteration. For updating the dictionary, the problem with respect to D becomes convex, that is, at the q t h iteration, there is,

D q = arg min D∈D 1 2 X -DA q 2 F .
By writing DA q in the form of sum of outer product (3.24), we solve this problem via a coordinate descent algorithm. More precisely, for updating each atom, the method of SVD is applied.(See Section 3.3 for details). Simultaneously, this process updates all the non-zero coefficients in the i t h row of A q . Therefore, the proposed dictionary learning with AcMIQP for sparse coding can be regarded as an exact 0norm resolution of the K-SVD algorithm.

The summary of this method is illustrated in Algorithm 10. The proof of convergence of the algorithm is out of scope due to the difficulty in the analysis of the used MIQP solver. However, it is guaranteed that at each step, with sparse coding and with dictionary updating, the objective value reduces.

Algorithm 10 Dictionary learning algorithm via AcMIQP.

Input: Signals for training X, target sparsity k, step size ρ for updating approximate A by proximal method, coefficient η for optimizing M, number of iteration for dictionary learning N d and for proximal method N p Output: Optimized dictionary D and sparse coefficient matrix A.

1: Initialization of the dictionary D 0 and the coefficients matrix A 0 = 0, 2: for q = 1 to N d do 3:

Initializing X by proximal method:

4:

A q-1,0 = A q-1 5:

for p = 1 to N p do 6:

Updating sparse approximation via

A q-1,p = prox δ S k A q-1,p-1 -ρ (D q-1 ) T D q-1 A q-1,p-1 -(D q-1 ) T X 7:
end for

8: A q-1 = A q-1,N p 9: M = η A 1,∞ , 10:
Optimization of A via MIQP solver:

11:

for i = 1, . . . , l do 12:

Initialization of α α α + i , α α α - i , z i and then v, 13:

α α α + i = max(0,α α α q-1 i ) 14:
α α α - i = max(0, -α α α q-1 i

) 15:

z i = abs(sign(α α α q-1 i )) 16: v = (α α α + i T ,α α α - i T , z T i ) T 17: Solving MIQP problem min v 1 2 v T Qv + c T v subject to A i n v ≤ b i n , l ≤ v ≤ u, v j ∈ {0, 1}, . . .∀ j ∈ I. 18:
Computing the solution of α α α

q i = v[1 : m] -v[m + 1 : 2m] 19:
end for 20:

Update D q with the coordinate descent algorithm (introduced in Algorithm 7) 21: end for

Image denoising based on the exact 0

In the previous section, we introduced the AcMIQP algorithm for solving problem based on exact 0 -norm. Furthermore, a dictionary learning method by applying AcMIQP algorithm for sparse coding was proposed. We showed that the good performance in the signal recovery and that the learned dictionary is of good quality.

In this section, we address image processing with real well-known natural images, such as Barbara and Lena.

Sparse representation is proved to be a good model in image processing, including image denoising [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], image inpainting [START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF], and image deblurring [START_REF] Dong | Image deblurring and super-resolution by adaptive sparse domain selection and adaptive regularization[END_REF]. The performance of sparse model on application of image processing becomes an important criterion for evaluation of dictionary learned by solving specialized sparse representation problem. Thus, hereafter, the proposed AcMIQP for exact 0 based dictionary learning will be used for image denoising on well-known images.

Sparse representation for image denoising

In a natural image, spatial smooth structures appear much more frequently than highly non-smooth and discontinuous structures. This fact acts as a proof of sparseness in image representation [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF]. Hence, image can be formulated as a linear combination of some few atoms of dictionary. However, because of imperfection in imaging device, poor illumination and information loss during signal transmission, the observed image is contaminated by noise [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF]. Thus, the application of image denoising is of great interest, and moreover, being an inverse problem, it provides an evidence for testing image processing algorithms and techniques, with a particular interest in demonstrating the relevance of sparse representations.

Considering an image x contaminated by a white Gaussian noise with zeromean and standard deviation σ, the observed image can be expressed by

y = x + e,
where y indicates the image of observation and e the Gaussian noise. When applying sparse representation model, the image x can be modeled as a linear combination of only a few atoms in a dictionary D, namely x = Dα α α with α α α 0 ≤ k where α α α denotes the sparse coefficient vector and D can be explained as "atomic" images (which is interpret as periodic-table in field of chemistry [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF]). By using this relation, the observed image y can be modeled by

y = Dα α α + e.
With knowledge on the energy bound of noise, namely e 2 ≤ , where is defined by = 1.15nσ 2 in [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF], the problem of image denoising can be formulated as arg min

α α α α α α 0 subject to y -Dα α α 2 2 ≤ . (4.21)
Therefore, this problem consists in finding the sparsest representation of the underlying image x under the condition that the reconstructed one obtained by Dα α α has no more distance than with the noisy observation y. In fact, the representation is quite sparse, namely α α α 0 ≤ k. In this way, the image will be component-wise denoised which is proved to be more efficient [START_REF] Hyvärinen | Image denoising by sparse code shrinkage[END_REF].

Alternatively, one can formulate the image denoising problem by adopting the Lagrangian function associated to the constrained problem (4.22), that is, arg min

α α α λ α α α 0 + 1 2 y -Dα α α 2 2 , (4.22)
where λ is the parameter balancing the data-fitting and the sparsity of α α α. The bigger λ is, the sparser the representation α α α is, which may remove some useful structural information of the image, in addition to noise. As a result, an adequate λ setting is important for achieving good performance in image denoising.

Besides, the corresponding sparsity constrained problem can be also reformulated as arg min

α α α, y-Dα α α 2 2 ≤ 1 2 y -Dα α α 2 2 , subject to α α α 0 ≤ k. (4.23)
where the sparsity level k should be explicitly prefixed by the user.

Problems (4.21), (4.23) and (4.22) are exactly the problems of sparse representation with 0 promoting sparsity. When the dictionary D is given, these problems can be solved via algorithms introduced in Section 3.2. If the optimal value α α α is obtained, then, the underlying image can be constructed directly by

x = D α α α. (4.24)
In the above discussion, the dictionary D is assumed to be known. However, in practice, the determination of D plays an important role in image denoising [ZXY + 15]. At the beginning of the emergence of sparse representations, dictionaries were predefined. Some special designed wavelet transforms [START_REF] Grace | Adaptive wavelet thresholding for image denoising and compression[END_REF] are used for forming a dictionary, for instance, Gabor wavelet [START_REF] Sing | Image representation using 2D Gabor wavelets[END_REF] contains diverse shifts of all entries in an orthogonal wavelet which will produce a shift-invariant representation; or curvelet [START_REF] Starck | The curvelet transform for image denoising[END_REF] is used for extracting discontinuous structures of an image. The application of those dictionaries gained significant improvement in image denoising (see [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF] and therein). Then, with the theoretical development in dictionary learning [MBPS09, AEB + 06, BJQS14], data-driven dictionaries [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] have been extensively used for image denoising.

Learned dictionary processes an image not as a whole but on small patches (8×8 Then by averaging the pixels value, the whole image can be reconstructed.

is
The above method reconstructs a whole image by removing noise in each patch.

Moreover, Elad et al [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] proposed a method for handling directly the entire image, by formulating the problem of image denoising as

arg min A,X,D λ X -Y 2 2 + i ω i α α α 0 + i Dα α α i -R i X 2 2 , (4.26)
where the first term X -Y 2 2 measures the proximity of the observed image Y and its reconstructed one X, and R i is the matrix that projects the whole image to the i t h patch. This problem is solved by firstly applying the K-SVD algorithm to learn an optimal dictionary D. Then, with the known D, the method alternatively solves the problem with respect to A with fixed X, and with respect to X with fixed A. The closed-form solution of X is given by

X = λI + i R T i R i -1 λY + i R T i Dα α α i , ( 4

Large-scale dictionary learning

This part considers the set of high-quality images in order to construct a unique global dictionary that will serve to denoise every image. More than ≈ 1.6 × 10 4 overlapping patches of size n = 8 × 8 from the images are extracted to get a single training dataset denoted Y. The number of atoms is set to p = 256 and the sparsity level is k = 20 (these parameters are determined by preliminary experiments and corroborated by other studies, such as [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]).

The dictionary learning algorithm with AcMIQP for sparse coding and coordinate descent algorithm for dictionary updating (Algorithm 10) is executed for 30 iterations to learn the global dictionary. This number of iterations is more than enough for convergence, as illustrated in Figure 4.2. This figure illustrates the learning curve, namely the evolution of the objective value at each iteration. It shows how the AcMIQP based dictionary learning algorithm converges faster than the other methods, namely K-SVD with OMP, proximal method and SOUPDIL algorithm. Moreover, the limit objective value (here indicating the value obtained after 30 iterations, which is proved to be more than sufficient for converging) of the AcMIQP based algorithm is the smallest comparing to the other methods. Consequently, we can conclude that AcMIQP can yield more exact solution with few iterations.

To measure the quality of the dictionaries, we consider the correlation, measured with the inner product between each pair of atoms of dictionary, thus measuring how much two atoms in the dictionary are similar. This fundamental information allows to define more powerful measures, such as the coherence and Babel function [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF][START_REF] Honeine | Analyzing sparse dictionaries for online learning with kernels[END_REF]. The coherence measure of a given dictionary, defined by the maximum absolute inner product between two distinct atoms, provides strong insights on the capacity of the dictionary to recover sparse signals. For instance, it is shown in [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] that a µ-coherence dictionary can recover a k-sparse signal if µ < 1 2k-1 . It is well known that the OMP algorithm (e.g. K-SVD) often provides dictionaries with high coherence, and most atoms are highly correlated. To over- set to assure the sparse level.

In order to understand the influence of the noise level on the results of the proposed method, we consider additive Gaussian noise of different standard deviations In tests, we compare the following methods: sparse coding using OMP, proximal method and AcMIQP, and dictionary updating with coordinate descent algorithm introduced in the K-SVD algorithm. The performance of the reconstruction accuracy is given in Table 4.3. These results show that the proposed method outperforms OMP and proximal methods in a high noise level. For the large-scale (global) dictionary learning setting, MIQP provides important enhancements, with an average improvement of 1.79 with regard to the proximal method, and 3.73 with regard to the OMP algorithm. It is worthy noting that the enhancement is significant since, the parameters were optimized for OMP, as recommended in [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. According to the experiments, we can conclude that the AcMIQP can reconstruct the signal more accurately with a high tolerance to noise.

Adapted dictionary learning

In this part, the dictionary is trained on the corrupted image under scrutiny, and then used to denoise it; the dictionary is then "adapted" to the image at hand. As in the last section, the signal matrix is created in the same way using overlapping patches. For each corrupted image X, an adapted dictionary is trained on it and then used for denoising the same image.

All three methods, OMP, proximal method and SOUPDIL, are compared with the proposed AcMIQP based dictionary learning method. Moreover, we consider also a variant of K-SVD with OMP, where the signals are pre-centered (subtracting the image mean) prior to learning the dictionary [MBP + 14]; connections between centered and uncentered data are studied in [START_REF] Honeine | An eigenanalysis of data centering in machine learning[END_REF]. In the experiments, SOUPDIL is In this part, the error-constrained optimization problem (4.23) is used for sparse coding. The method of realization is described in the large-scale dictionary learning. In order to ensure the sparsity of the signals, the upper bound k is set to 20 for the proximal and AcMIQP methods, as in the first setting. By fixing the dictionary updating method to coordinate descent algorithm (Algorithm 7, SVD method is used for updating each atom in the dictionary) in all the methods, this allows to have a fairly comparable setting to analyze and compare the performance of the sparse coding methods. The number of atoms is set to p = 256 for OMP, as suggested in [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] where extensive experiments were conducted. The number of atoms for the proximal method is set to p = 65, which is obtained from a set of 14 candidate values {50, 55, 60, 65, . . . , 110, 150, 200, 256, 300} that encloses the most used values in the literature. The same value is used for AcMIQP, which is a less favorable situation for our method. The total number of iterations is still 30 for the two-step sparse coding and dictionary updating. The SOUPDIL method uses the same parameter setting as recommended in [START_REF] Ravishankar | Efficient sum of outer products dictionary learning (soup-dil) and its application to inverse problems[END_REF] after extensive experimental analysis.

With the learned dictionaries, the same reconstruction model (4.27) is used for obtaining the denoised image. Table 4.4 gives the denoising accuracy in terms of PSNR by using the three aforementioned dictionary learning methods. We notice that the influence of data pre-centering is not always positive. It is observable that MIQP can outperform the K-SVD and proximal methods almost in all cases. On average over all five images, the proposed method carries out an improvement of 0.45 with respect to OMP, and 0.08 with respect to the proximal method. These improvements are important since, on one hand, PSNR is a logarithmic-scale measure and, on the other hand, the parameters were optimized for OMP (e.g. , c = 1.15, p = 256) and for the proximal method (p = 65). Even compared with the state-of-the-art dictionary learning algorithm SOUPDIL, MIQP has comparable performance.

Analysis on the computational complexity

In spite of the great performance of the proposed AcMIQP method on all images and compared to all the other methods, it has high computational complexity in implementation. Because we have different sizes of the training data in each setting (global dictionary learning and adapted dictionary training for each image), the training time is not comparable. In the following, we focus on the average time of a single image. While the OMP algorithm and the proximal method require only a couple of minutes for completing the dictionary learning, AcMIQP needs about one hour. See also Table 4.1 for results obtained on synthetic data. However, recent advances in MIQP solvers allow to reduce this gap.

Indeed, while the computational complexity remains the Achilles' heel of such methods, great improvements are being carried out these days on MIQP solvers.

For instance, the new Gurobi Optimizer v8.1 (released in October 2018) is 2.8 times faster overall on MIQP problems, than v8.0 (released in May 2018), which is more than 220% faster than the one used in this thesis (v7.0 released in October 2016).

Moreover, new advances in solvers are exploiting more and more the modern architectures and multi-core processors. Finally, currently available off-the-shelf solvers, such as Gurobi and CPLEX, do not have GPU implementations, which could also provide important computational improvements.

Conclusion

In this chapter, we proposed an exact optimization method AcMIQP for the exact 0 based sparse coding. Thanks to recent advances in linear programming techniques, as well as more powerful hardware, the speed of solving MIQP problems has been greatly improved. Furthermore, by introducing additive constraints and an appropriate proximal initialization, it was proved that it is feasible to use MIQP for sparse coding with the proposed AcMIQP method. In conjunction with a dictionary update, such as coordinate descent method used in K-SVD, we proposed a dictionary learning algorithm with exact 0 based sparsity. Though, the AcMIQP method had much more time complexity in implementation comparing with the approximate methods, the feasibility of the method was proved for large-scale well-known images. Moreover, the image denoising experiments showed the advantage of the proposed AcMIQP method. Furthermore, the high noise-tolerance of our method was demonstrated on both the large-scale and the adapted dictionary learning settings.

Therefore, we have demonstrated that the exact 0 optimization problem in dictionary learning can be solved for image processing, working on real images. While having good performance amelioration, the Achilles' heel of the proposed method is the computational complexity. However, great improvements are being carried out these days on MIQP solvers, with more than 220% speed enhancement in a single year (e.g. Gurobi Optimizer v8.0 versus v7.0).

In the following chapter, we will study the problem of dictionary learning by considering learning incoherent dictionary. For handling the resulting constrained optimization problem, we will exploit the classic strategy of iteratively alternating two steps, sparse coding and dictionary updating. With the similarity of optimization problem formulation in sparse coding, the proposed AcMIQP algorithm will still be used. When updating the dictionary, a new algorithm will be proposed to address the incoherence constraint. The coherence of a dictionary corresponds to the largest correlation between its atoms (e.g. null coherence for dictionaries of orthogonal elements). Beyond being elementary and very simple to compute, the coherence is intimately related to the sparsity level and the relevance of the resulting sparse representation. For example, the uniqueness of solution in the sparse representation problem lies on the condition of incoherence of the dictionary [GJB15, Tro04, DH01]. Moreover, the coherence was used as a criterion in evaluation of the dictionary learned based on the proposed AcMIQP algorithm (in Section 4.1.2). In sparse representation, the coherence measure is crucial and acts as a fundamental measure to characterize a dictionary. Indeed, several theoretical studies have demonstrated the prominence of having incoherent dictionaries, namely dictionaries having a low coherence measure [Hon15b, [START_REF] Honeine | Analyzing sparse dictionaries for online learning with kernels[END_REF][START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. Incoherent dictionary learning, as an extension of the generic dictionary learning, aims at minimizing the reconstruction error by imposing simultaneously the sparsity on the coefficients and the incoherence of the dictionary. For this purpose, several incoherent dictionary learning algorithms have been proposed, within three major strategies: either adding a decorrelation step after dictionary updating at each iteration, such as INK-SVD and related algorithms [START_REF] Mailhé | INK-SVD: Learning incoherent dictionaries for sparse representations[END_REF][START_REF] Barchiesi | Learning incoherent dictionaries for sparse approximation using iterative projections and rotations[END_REF], or introducing an additional regularization on the coherence in the optimization problem [START_REF] Ramírez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF][START_REF] Bao | A convergent incoherent dictionary learning algorithm for sparse coding[END_REF].

In this chapter, we formulate the dictionary learning problem with an explicit constraint on the coherence. Thus, we consider solving an optimization problem with quadratic objective function and quadratic constraints based on the 0norm. As the generic dictionary learning, the solution can be reached by iteratively processing two alternating processes: sparse coding and dictionary updating. For sparse coding, we use the proposed AcMIQP algorithm. For the dictionary updating, we propose a novel resolution method that combines the alternating direction method of multipliers (ADMM) and the method of extended proximal alternating linearized minimization. The relevance of the proposed incoherent dictionary learning method is demonstrated with experiments on real data, compared to wellknown methods. Moreover, we provide a theoretical analysis on the convergence.

Coherence and sparse representation

Sparse representations require finding the sparest code of the given signal over a determined dictionary. The value of the coherence of the dictionary has a significant effect in the resolution of sparse coding and in the estimation of the dictionary, as well as the evaluation of its quality. The coherence is thus obviously a major tool for the dictionary analysis and synthesize in sparse representations.

Definition of coherence

The coherence is defined as the greatest correlation, in absolute value, between two distinct atoms of the dictionary under scrutiny. When dealing with unit-norm atoms, the coherence of a given dictionary D = [d 1 , . . . , d m ] ∈ R n×m , is defined as ,

µ = max i = j |〈d i , d j 〉| = max i = j |d T i d j |. (5.1)
For a dictionary, the coherence is a fundamental quality of assessment which estimates how much two atoms in dictionary are correlated. Specially, when µ = 0, the dictionary become orthogonal. For redundant representations, the dictionary is overcomplete with m > n, then, the coherence will be bounded by

µ ≥ m-n n(m-1) .
(5.2)

The limit of dictionary coherence is met when the dictionary is an optimal Grassmannian frame [WGL17, SWDS15]. Consider for example, a dictionary of 256 atoms living in a 64 dimensional space, its coherence cannot be less than 0.108, i.e., the angle between any two distinct atoms of the dictionary cannot be bigger than 83.8 • .

Letting G = D T D be the Gram matrix of D, then the coherence of D is obtained by searching the maximal absolute value of off-diagonal elements in G, namely µ =

sup i = j |G i j |.
The coherence of the dictionary is also related to the RIC defined in Chapter 1. The relationship is addressed in [START_REF] Honeine | Analyzing sparse dictionaries for online learning with kernels[END_REF] with

δ k = (m -1)µ.
The coherence considers only the extreme correlation of atoms rather the global situation. For avoiding this shortcoming, another definition called cumulative coherence (or babel function) is created [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] µ 1 (g ) = max

|J|=g max i ∈J J |〈d i , d j 〉|, (5.3)
for any integer g , and where J is a subset of {1, . . . , m}. This function is nondecreasing. When the function varies little with g increasing, the dictionary is more likely incoherent. For a given dictionary, its coherence µ and cumulative coherence µ 1 (g ) satisfy the relation µ 1 (g ) ≤ g µ.

Coherence of dictionary in sparse representations

The dictionary coherence has a great effect in sparse representations. First of all, the uniqueness of solution to the problem of sparse representation can be ensured by satisfying a condition on µ, that is, Theorem 5.1.1. (Uniqueness-Coherence [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF]) If a problem of signal estimation

x = Dα α α has a solution α α α with α α α 0 < 1 2 (1 + 1 µ )
, then this solution is necessarily the sparsest possible.

This theorem can be regarded as an inference of Theorem 2.1.1 by considering the relation between spark(D) and µ(D),

spark(D) ≥ 1 + 1 µ(D) .
Via Theorem 5.1.1, the uniqueness of the solution to the problem (4.7) is ensured when the sparsity obeys that k < 1 2 (1 + 1 µ ). Similarly, the uniqueness of the solution can also be guaranteed by the cumulative coherence of D by considering the rela-

tionship spark(D) ≤ min 1≤g ≤m {g | µ 1 (g -1) ≥ 1}.
The importance of the coherence measure to characterize dictionaries has been demonstrated in several works [Hon15b, [START_REF] Honeine | Analyzing sparse dictionaries for online learning with kernels[END_REF][START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. For example, it is proven in [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] that orthogonal matching pursuit and basis pursuit can correctly recover the signal under the condition

max i ∈J D T J d i < 1.
Or if the signal is not strictly sparse, it can be approximated with a k-sparse vector under the condition that k ≤ 1 3µ or µ 1 (k) ≤ 1 3 . Furthermore, Gribonval et al give the cumulative coherence condition to make sure the signal with some noise can be exactly recovered (see Proposition 3 in [START_REF] Gribonval | Sparse and spurious: Dictionary learning with noise and outliers[END_REF]). Although these conditions may not be the same for whichever the sparse coding methods are, the importance of incoherent dictionary learning is undoubted.

More often, the low coherence between atoms or sub-blocks of dictionary makes the dictionary gain the capability of discrimination in classification [WJY + 16, TLZ + 19, LLF14]. In fact, some other discriminative dictionary learning via sparse model can be regarded as implicitly learning a dictionary with low coherence [ZL10, MLB + 08, MPS + 09]. In this sense, incoherent dictionary learning attracts lots of attentions. Great achievement on learning incoherent dictionary is gained and there is still much space for further developments.

Methods of learning incoherent dictionary

The classical dictionary learning problem, with given matrix X = [x 1 , . . . , x i , . . .

α i ∈R p 1 i =1 1 2 x i -Dα α α i 2 2 subject to α α α i 0 ≤ k, i = 1, . . . , , (5.4)
where the dictionary D is restricted in the constraint

C = D ∈ R n×m d T j d j = 1, ∀ j = 1, . . . , p ,
in order to prevent the 2 -norm of dictionary's atoms from being arbitrarily large, which leads to arbitrarily small decomposition coefficients in X. However, by this problem formulation, the resulting dictionary does not guarantee excellent performance, because its atoms can be arbitrarily correlated. Thus, to learn dictionary with low coherence is of great interest.

For learning incoherent dictionary, two main strategies have been considered: 1) Adding a decorrelation step following the phase of dictionary updating at each iteration. This strategy is applied in dictionary learning algorithms such as INK-SVD and the incoherent dictionary learning algorithm via iteratively projection and ro-

p i =1 ( d i 2 -1) 2 [RLS09]
or by a normalization step following the dictionary updating [START_REF] Bao | A convergent incoherent dictionary learning algorithm for sparse coding[END_REF]. In this chapter, we investigate a new strategy that models the problem by proposing the explicit constraints on the coherence of the dictionary. By this method, the problem of incoherent dictionary learning becomes a constrained optimization problem with quadratic objective function and quadratic constraints.

Dictionary learning with additive decorrelation

The conventional dictionary learning problem (5.4) is solved, via iteratively processing two alternating steps, sparse coding and dictionary updating as described in Chapter 3. A first approach to decorrelate was proposed in the implementation1 of the conventional K-SVD (Algorithm 7), where the authors removed the atoms that are highly correlated with each other. However, the resulting algorithm does not ensure an optimal reconstruction performance.

For overcoming the difficulties, the leading method INK-SVD [START_REF] Mailhé | INK-SVD: Learning incoherent dictionaries for sparse representations[END_REF] and incoherent dictionary learning algorithm via iterative projections and rotations [START_REF] Barchiesi | Learning incoherent dictionaries for sparse approximation using iterative projections and rotations[END_REF] are developed from the well-known K-SVD algorithm. The two algorithms take the same strategy that adds a supplementary step following the phase of dictionary updating at each iteration. These two algorithms are described in the following.

The INK-SVD is defined as follows. After obtaining the updated dictionary D t at the t t h iteration, Mailhé et al proposed in [START_REF] Mailhé | INK-SVD: Learning incoherent dictionaries for sparse representations[END_REF] to find the closest dictionary to D t with unit-norm atoms and a coherence below a predefined threshold µ c . The problem is formulated as

D t = arg min D∈D µ D -D t 2 F ,
where

D µ = {D ∈ R n×m | µ(D) ≤ µ c , d i 2 = 1, ∀i = 1, .
. . , m}. To find the optimal solution, an iterative algorithm is proposed by identifying the sub-dictionary (in the same spirit as K-SVD) and decorrelating pairs of atoms with a greedy algorithm.

More precisely, at each iteration, a pair of atoms (d i , d j ) with inner product larger than µ c is selected and rotated symmetrically with respect to their mean vector until the angle between them reaches θ, with cos(θ) = µ c . The dictionary D t can be reached within a finite number of iterations. The algorithm of decorrelation is summarized in Algorithm 11.

The iterative projections and rotation (IPR), proposed by Barchiesi et al in The resolution of this problem consists of two steps:

• Decorrelation of atoms by projection method.

In this step, the Gram matrix G = D T D is considered, rather than processed directly with D. For obtaining the dictionary D with the required low coherence, it is necessary to project G in the space defined by two constraints, named in [START_REF] Barchiesi | Learning incoherent dictionaries for sparse approximation using iterative projections and rotations[END_REF] the structural constraint

G µ c = {G ∈ R m×m | G T = G, diag(G) = 1 m , max i = j g i j ≤ µ c }, (5.5)
and the spectral constraint

G λ = {G ∈ R m×m | G T = G, eig(G) ≥ 0, rank(G) ≤ n}.
(5.6)

For a given G, the projection on G µ c can be easily gained by the operations:

g i j = 1 if i = j sign(g i j ) min(|g i j |, µ c ) if i = j.
For getting the projection on G λ , it computes the eigenvalues λ i , i = 1, . . . , m, by eigendecomposition of G = V∆V T where D is a diagonal matrix formed by the eigvenvalues of G and V formed by the corresponding eigenvectors, and then, retains the largest n ones with ensuring their positivity and setting to zero the remaining ones. However, after projection on G λ , the incoherence of D cannot be guaranteed any more. Hence, the projection on the intersection of the two constraint spaces is not reached by only one projection, but by an iterative process.

When the Gram matrix G belongs to this intersection, namely satisfying G ∈ G µ c ∩ G λ , then D can be recovered by

D = ∆ 1 2 V T .
(5.7)

• Dictionary rotation for minimizing the reconstruction error while keeping the dictionary with the required coherence level.

The dictionary rotation is realized by introducing an orthogonal transform matrix W, by which the coherence quality is still held due to the fact that

Regularized incoherent dictionary learning

Methods from the second strategy seek to learn an incoherent dictionary by minimizing a regularized objective function, where the regularization term constrains the coherence. We describe in the following the most-known methods. The problem is formulated in [START_REF] Ramírez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF] as

arg min D∈R n×m ,A∈R m×l X-DA 2 F +λ 1 j i log(|α i j |+β)+λ 2 D T D-I m 2 F +λ 3 k (d 2 k -1) 2 , (5.8)
where the second term denotes the sparsity-promoting. This formulation is obtained by using a universal models [START_REF] Ramírez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF]. The β is the parameter regulating the universal model. The coefficients λ 1 , λ 2 and λ 3 are trade-off parameters. However, here we focus on the step of dictionary updating. The optimization problem with respect to D is of the form arg min

D∈R n×m X -DA 2 F + λ 2 D T D -I m 2 F + λ 3 k (d 2 k -1) 2 ,
where the second term measures the correlation of two distinct atoms in D, which are the off-diagonal elements of D T D, and the last term makes sure that columns of D are of norm close to 1. This optimization problem can be solved with the method of optimal coherence-constraint directions (MOCOD) [START_REF] Ramírez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF], inspired from the method of optimal direction (MOD), which updates D by

D t +1 = X(A t +1 ) T + 2(λ 2 + λ 3 )D t A t +1 (A t +1 ) T + 2λ 2 (D t ) T D t + 2λ 3 diag((D t ) T D t ) -1
.

The MOCOD method is proved to outperform the MOD method in image reconstruction.

In [START_REF] Bao | A convergent incoherent dictionary learning algorithm for sparse coding[END_REF], the incoherent dictionary learning problem is formulated by introducing only the coherence regularization, namely the Frobenius norm of the difference of Gram matrix and identity, arg min

D∈D,A∈R m×l 1 2 X -DA 2 F + λ 1 A 0 + λ 2 2 D T D -I m 2 F .
The incoherent dictionary is learned via a hybrid alternating proximal method, of which the dictionary is updated atom-by-atom using the operation

d t +1 j = prox δ D λ t j d t -1 λ t j ∇ d j Q(A t , D t j ) , where Q(A, D) = 1 2 X-DA 2 F + λ 2 2 D T D-I m 2 F and D t j = [d t +1 1 , . . . , d t +1 j -1 , d t j , . . . , d t m ],
and then the dictionary is normalized after the dictionary updating at each iteration.

Similarly, Abolghasemi et al [START_REF] Abolghasemi | Fast and incoherent dictionary learning algorithms with application to fMRI[END_REF] tackled the problem with the same coherence regularization as in [START_REF] Bao | A convergent incoherent dictionary learning algorithm for sparse coding[END_REF]. However, they proposed an incoherent dictionary learning algorithm with dictionary updating by a gradient descent method. In addition, the coherence regularization was also measured by the sum of 1 -norm of every two different atoms [START_REF] Li | An efficient algorithm for incoherent analysis dictionary learning based on proximal operator[END_REF][START_REF] Li | Incoherent dictionary learning with log-regularizer based on proximal operators[END_REF]. Moreover, for some tasks such as classification, the Fisher criterion [START_REF] Li | Learning low-rank and discriminative dictionary for image classification[END_REF] can be regarded as a coherence regularization, which makes the sub-dictionary of different class coherent. Incoherent dictionary learning algorithms of the second strategy achieve good performance in data reconstruction [START_REF] Abolghasemi | Fast and incoherent dictionary learning algorithms with application to fMRI[END_REF][START_REF] Li | Incoherent dictionary learning with log-regularizer based on proximal operators[END_REF], classification [BQJ14, LLZ12, LHT + 18] and object recognition [START_REF] Li | Learning low-rank and discriminative dictionary for image classification[END_REF][START_REF] Zhang | Learning structured low-rank representations for image classification[END_REF]. However, they suffer from a major issue: it is not possible to constraint exactly the coherence level to a fixed value, because the relation between it and the regularization trade-off parameter is unknown.

Towards incoherence-constrained dictionary learning

The third strategy considers the simplest way to formulate the problem, by adding the constraints of coherence and unit norm of the dictionary elements into the generic dictionary learning problem. The coherence of the dictionary is constrained with the inequality |d T p d q | ≤ µ c , ∀p, q ∈ {1, 2, . . . , m}, p = q, (5.9)

where µ c is the predefined coherence level. The unit norm of the dictionary's atoms is modeled by the equality d T p d p = 1 ∀p = 1, 2, . . . , m.

(5.10) Thus, the problem of incoherent dictionary learning can be resumed as a constrained optimization problem with quadratic objective function and quadratic constraints min

D∈R n×m ,α α α i ∈R m 1 l l i =1 1 2 x i -Dα α α i 2 2 subject to        |d T p d q | ≤ µ c , ∀p, q ∈ {1, 2, . . . , m}, q = p d T p d p = 1, p = 1, . . . , m x i 0 ≤ k, i = 1, . . . , l .
(5.11)

The problem of estimating simultaneously A and D is non-convex and non-smooth because of the sparsity-prompting 0 -norm and the constraints. To the best of our knowledge, there is no work on incoherent dictionary learning by solving the problem with explicit constraints on dictionary coherence and its unit norm. To solve this constrained optimization problem, we take advantage of recent developments in optimization problem with orthogonality constraints, with the augmented Lagrangian method and the alternating proximal minimization method.

The optimization problem with orthogonality constraints has been recently addressed in physics [UWY + 15], mathematics [START_REF] Lin | The augmented Lagrange multiplier method for exact recovery of corrupted low-rank X BIBLIOGRAPHY matrices[END_REF] and information science [START_REF] Chen | An augmented Lagrangian method for 1-regularized optimization problems with orthogonality constraints[END_REF]. The Lagrangian multiplier method [START_REF] Bertsekas | Nonlinear programming[END_REF] is frequently used to deal with such a problem [LCWM09, UWY + 15]. However, it is not always easy to solve the Lagrangian function by satisfying the first order optimal condition. In [UWY + 15], the Kohn-Sham problem was reformulated by the Lagrangian multiplier method, and the proximal gradient method was then proposed to solve the Lagrange function.

Moreover, it was proven that the algorithm has good convergence property. The orthogonality constrained optimization problem were also solved via the augmented Lagrangian method [START_REF] Chen | An augmented Lagrangian method for 1-regularized optimization problems with orthogonality constraints[END_REF][START_REF] Zhu | Nonconvex and nonsmooth optimization with generalized orthogonality constraints: An approximate augmented Lagrangian method[END_REF]. Compared with the Lagrangian method, the penalty method shows more stability [START_REF] Bertsekas | Nonlinear programming[END_REF]. However, the reformulated problem can be non-convex and non-smooth, which makes the problem hard to tackle. In [START_REF] Chen | An augmented Lagrangian method for 1-regularized optimization problems with orthogonality constraints[END_REF], the alternating proximal method was combined with the augmented Lagrangian method and the existence of the sub-sequence to a KKT point was proven.

The new proposed algorithm was then applied in compressed mode for variational problems in physics, illustrating the effectiveness and efficiency of the algorithm. In [START_REF] Zhu | Nonconvex and nonsmooth optimization with generalized orthogonality constraints: An approximate augmented Lagrangian method[END_REF], an extended proximal alternating linearized minimization method was introduced to solve the Lagrangian function, and its convergence was proven based on the theory of the Kurdyka-Łojasiewicz inequality property [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF].

Exact 0 based incoherent dictionary learning

This section gives a method to solving the exact 0 based incoherent dictionary learning problem (5.11). Similar to all other dictionary learning algorithms presented in Section 3.4, we learn the incoherent dictionary via iteratively processing two alternate steps: sparse coding and dictionary learning. The sparse coding problem of (5.11) is the 0 -norm constrained problem, which is exactly the same as in the generic sparse coding (4.7). Therefore, we use the AcMIQP algorithm proposed in Section 4.1.2 to solve exactly the 0 constrained sparse coding problem. With a fixed sparse code, the problem with respect to the dictionary D becomes a non-convex constrained optimization problem. For solving this problem, the augmented Lagrange method and proximal alternating linearized minimization method are used.

The convergence of the proposed algorithm is analyzed in the following.

Augmented Lagrangian method

The augmented Lagrangian method is used in this part to reformulate the incoherence-constrained dictionary learning problem. The optimization problem with respect to D is: where

min D∈R n×m 1 l l i =1 1 2 x i -Dα α α i 2 2 subject to        |d T p d q | ≤ µ c , ∀p, q ∈ {1,
S G = G ∈ R m×m |g i j | ≤ µ c , i , j = {1, 2, . . . , m}, i = j .
Let δ S G (G) be the indicator function on this set, namely

δ S G (G) = 0, if G ∈ S G
+∞, otherwise.

(5.14)

The constrained optimization problem can be solved by considering the augmented Lagrangian function:

L (c1,c2) (D, G,λ λ λ, H) = 1 2 X -DA 2 F + m p=1 λ p (d T p d p -1) + c 1 2 m p=1 (d T p d p -1) 2 + tr(H(G -D T D)) + c 2 2 G -D T D 2 F + δ S G (G), (5.15)
where λ λ λ = [λ 1 , . . . , λ m ] and H are respectively the vector and matrix associated to the equality constraints on the diagonal of D T D and on G, c 1 and c 2 are the positive penalty parameters (the augmentation). When these parameters grow into infinity, the optimal solution of the original problem (5.13) can be reached.

Therefore, the optimization problem becomes: min D∈R n×m ,G∈R m×m L (c1,c2) (D, G,λ λ λ, H).

(5.16)

It is not the standard augmented Lagrangian method (where the objective function is convex and has only one term, in most case, the constraints are closed convex set).

While our problem is non-convex and non-smooth, it is still reasonable to consider the inexact ADMM framework [START_REF] Bertsekas | Nonlinear programming[END_REF]. The resulting algorithm is illustrated in Algorithm 12.

As presented in Algorithm 12, the inexact augmented Lagrangian method operates in three alternating steps: In Step 1, the primal variables are computed by solving, as explained in next section, the optimization problem,

(D i , G i ) = arg min D,G L (c i 1 ,c i 2 ) (D, G,λ λ λ i , H i ), (5.17

Proximal alternating linearized minimization

The problem (5.17

∇ D h(D, G) -∇ D h( D, G) F ≤ L D D -D F ∇ G h(D, G) -∇ G h(D, Ḡ) F ≤ L G G -Ḡ F , (5.19)
for all (D, D) and (G, Ḡ).

(ii) L (c i 1 ,c i 2 ) (D, G,λ λ λ i , H i ) satisfies the Kurdyka-Łojasiewicz inequality [START_REF] Zhu | Nonconvex and nonsmooth optimization with generalized orthogonality constraints: An approximate augmented Lagrangian method[END_REF]. The proof is given in Section 5.3.5.

The problem (5.17) can now be solved by alternating the optimization problems with respect to D and G, respectively:

                                         D i , j = arg min D∈R n×m f (D) + h(D i , j -1 , G i , j -1 ) +tr (D -D i , j -1 ) T ∇ D h(D i , j -1 , G i , j -1 ) + t1 2 D -D i , j -1 2 F G i , j = arg min G∈R m×m g (G) + h(D i , j , G i , j -1 ) +〈G -G i , j -1 , ∇ G h(D i , j , G i , j -1 )〉 + t 2 2 G -G i , j -1 2 F ,
(5.20)

where 〈M 1 , M 2 〉 = tr(M T 1 M 2 ) is defined as the scalar product in the matrix space M n (R), t1 and t 2 are the coefficients associated respectively to the second order approximation terms. By considering that the function f is differentiable, the problem of estimating D can be expressed as:

D i , j = arg min D∈R n×m f (D i , j -1 ) + h(D i , j -1 , G i , j -1 ) +tr (D -D i , j -1 ) T ∇ D ( f (D i , j -1 ) + h(D i , j -1 , G i , j -1 ) + t 1 2 D -D i , j -1 2 F .
(5.21)

To solve this optimization problem, it is easy to update D by the method of gradient descent. As for the problem of estimating G in (5.20), the proximal method is applicable. By combining both steps, the solution of D and G can be achieved by the following process:

       D i , j = D i , j -1 -1 t 1 ∇ D f (D i , j -1 ) + h(D i , j -1 , G i , j -1 ) G i , j = prox 1 t 2 g G i , j -1 -1 t 2 ∇ G (h(D i , j , G i , j -1 )) , (5.22)
where the notation prox u f denotes the proximal operator, defined in Section 4.1.2, of the scaled function u f (also called the proximal operator of g with parameter u),

namely prox 1 t 2 f (v) = arg min x f (x) + 1 2t 2 x -v 2 .
(5.23)

The partial derivatives of the three parts of (5.18) are:

           ∇ D f (D i , j -1 ) = -(Y -D i , j -1 X)X T + 2D i , j -1 diag(λ i λ i λ i ) + 2c i 1 D i , j -1 diag(v j-1 ) ∇ D h(D i , j -1 , G i , j -1 ) = -D i , j -1 (H i + (H i ) T ) -c i 2 D i , j -1 G i , j -1 + (G i , j -1 ) T -2(D i , j -1 ) T D i , j -1 ∇ G h(D i , j , G i , j -1 ) = H i + c i 2 (G i , j -1 -(D i , j ) T D i , j ), (5.24)
where v j-1 = diag (D i , j -1 ) T D i , j -1 -I m denotes the vector with the entries valued by the diagonal of the matrix (D i , j -1 ) T D i , j -1 -I p . The two expressions diag(λ i λ i λ i ) and diag(v j-1 ) return the matrix with the diagonal filled by the elements in the vectors λ i λ i λ i and v j-1 , respectively. Moreover, G i , j can be computed as:

G i , j (i x , i y ) =        G(i x , i y ) if | G(i x , i y )| ≤ µ c sign( G(i x , i y ))µ c otherwise, (5.25)
where G = G i , j -1 -1 t 2 ∇ G h(D i , j , G i , j -1 ) , sign( G(i x , i y )) is the sign of the G(i x , i y ), i.e., +1 if G(i x , i y ) ≥ 0 and -1 otherwise. In these expressions, i x and i y are respectively the row and column indices of the matrix G, with i x , i y = 1, 2, . . . , m.

The proposed EPALM method for estimating D and G is summarized in Algorithm 13. For the completeness of the algorithm, we provide next a convergence analysis in terms of the subdifferential of the objective function in (5.17), as well as the choice of the parameters.

A necessary but not sufficient condition for x ∈ R n to be a minimizer of f is 0 ∈ ∂ f (x). Back to our optimization problem, the subdifferential of L (c i

1 ,c i 2 ) (D, G,λ λ λ i , H i ) at (D i , G i ), denoted by Θ i = (Θ i D , Θ i G )
and expressed as

Θ i = ∂L (c i 1 ,c i 2 ) (D i , G i ,λ λ λ i , H i ),
can be computed directly and the result can be written in form of:

   Θ i D = ∇ D f (D i ) + ∇ D h(D i , G i ) Θ i G = t 2 (G i -1 -G i ).
(5.26) Thus a solution of our optimization problem can be found when Θ i ∞ → 0. According to the formulation of Θ i in (5.26), D i is exactly the local optimal solution of the subproblem with respect to D, and the sequence G i , j is convergent, since G i -G i -1 F → 0. Besides, it is noticed that to guarantee that every bounded sequence generated by the proposed method converges to a critical point of L (c i 1 ,c i 2 ) (D i , G i ,λ λ λ i , H i ), the parameters c 1 , c 2 and the stepsizes t 1 , t 2 need to be appropriately chosen. The following can be noted:

• The initial positive penalty parameter c 0 1 and c 0 2 should be carefully chosen to avoid ill-conditioning, i.e., they should satisfy the second-order sufficient condition:

∇ 2 DD L(D i , G i ,λ i λ i λ i , H i ) > 0.
Due to the complexity of the derivative of a matrix function with respect to a matrix (the derivative of the function with respect to each element of the matrix being a matrix), we do not give the detail here.

• The convergence of the algorithm requires that the descent stepsize, i.e., 1 t 1 and 1 t 2 , should not be too large, satisfying t 1 > L D and t 2 > L G . noise, the output of the algorithm cannot be exactly D * . For this reason, the objective function for µ c > µ * are always lower than that when µ c = 0.6. However, when µ c < µ * , the situation is totally different, because µ * is out of the feasible region.

Therefore, a solution that satisfies the coherence constraint can be found, but the price to pay is an increase of the objective function, as well as the computational cost to converge. Consequently, by appropriately choosing the coherent parameter, an incoherent dictionary can be produced by this algorithm. Moreover, the smaller the target coherent parameter is, the greater the computational complexity will be.

Incoherent dictionary learning algorithm

We get an incoherent dictionary learning method by combining the dictionary updating method introduced in the previous section with a sparse coding method, such as the proximal method or the AcMIQP introduced in Chapter 4. By considering AcMIQP, we get an incoherent dictionary learning algorithm with an exact 0 optimization. The outline is illustrated in Algorithm 14. In the reminder of this section, we examine the relevance of two combinations: prox+EPALM, which combines the proximal method with the hybrid algorithm of ADMM and EPALM, and AcMIQP+EPALM, which combines AcMIQP with the proposed hybrid algorithm.

Algorithm 14 Exact 0 based incoherent dictionary learning algorithm (5.17)

Input:

The training data X, number of iterations N, all parameters and stropping criterion needed in inexact ADMM algorithm (Algorithm 12) and EPALM algorithm (Algorithms 13). Output: The solution D and A function INCOHERENT DICTIONARY LEARNING Initialisation of dictionary D 0 and sparse representation A = 0, for i = 1 to N do 1. Sparse coding with AcMIQP algorithm. 2. Dictionary updating by ADMM with EPALM (Algorithms 12-13).

end for end function

In the following, the dictionary learning algorithm is evaluated on the segment of image Barbara presented in Figure 4.1. The overlapping patches of size 8 × 8 (namely a signal is a vector of size 64) form the set of signals X. With the signals, a dictionary D is learned by using the proposed method (the hybrid algorithm of the ADMM and EPALM for dictionary updating, and either the proximal method or the AcMIQP for sparse coding) and compared to the other two comparative incoherent dictionary learning algorithms, INK-SVD and IPR. When both D and X are known, the sparse code A can be easily obtained using a sparse coding method, namely proximal method and AcMIQP for our algorithm, OMP algorithm for the other two methods. Then, the reconstructed image is obtained by doing the matrix multiplication X = DA. Consequently, we compare their performance by computing the peak signal-to-noise ratio (PSNR). The MIQP problem is solved by the software Gurobi Optimizer 8.1.0. The parameter settings of Gurobi are fixed as the default values except that the time limit is set to 0.5 seconds and the maximal iteration number to 1000. The initialization of MIQP is given by running the proximal method by setting the maximal iteration number to 200. When only the proximal method is used to sparse code, the iteration number is 1000. The number of atoms is set to p = 256 and the sparsity level k = 20 (the active atoms is less than 8%). The iteration number for learning a dictionary is set to 30, which is sufficient for the algorithms to converge, as shown in Figure 5.2. For the other two comparative methods, the parameter values are chosen as in the original papers [START_REF] Mailhé | INK-SVD: Learning incoherent dictionaries for sparse representations[END_REF][START_REF] Barchiesi | Learning incoherent dictionaries for sparse approximation using iterative projections and rotations[END_REF]. sparse coding and EPALM for dictionary updating has the fastest convergence and the value of limit is the smallest. It is worth pointing out that 30 iterations is sufficient for the algorithms to converge, even though the IPR algorithm shows some convergence unstability.

Figure 5.2 shows the distribution of the absolute inner product between each two atoms in the learned dictionary. Combined with the statistics given in Table 5.2, we notice that independently of the used sparse coding algorithm, the proposed method can achieve a dictionary with almost the target coherence parameter value, which is not the case of IPR. The proximal method combined with EPALM provides the smallest absolute average, which is an important property related to the socalled Babel function whose theory is well established [Hon15b, [START_REF] Honeine | Analyzing sparse dictionaries for online learning with kernels[END_REF][START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF].

However, this algorithm cannot beat the one with MIQP for sparse coding in terms of variance. The INK-SVD algorithm outputs as well a dictionary with almost the target coherence value, but with a higher variance. Nevertheless, INK-SVD updates the dictionary without considering the reconstruction error (see next paragraph).

For the IPR algorithm, the target coherence parameter value cannot be obtained even though it shows the least variance. Considering the distribution of absolute inner products between each two atoms in the learned dictionary, as illustrated in Figure 5.2, it is hard to tell which of the proximal method or the MIQP is better in combination with the proposed algorithm.

To analyze the reconstruction errors, we study seven coherence values µ c = {0.996, 0.966, 0.866, 0.707, 0.500, 0.259, 0.122}, i.e., the angle between any two atoms learning with the proposed method of ADMM and EPALM, the reconstruction performance improves with the coherence of dictionary decreasing, with the best results when µ c = cos(45 • ) with proximal method for sparse coding and µ c = cos(30 • ) with AcMIQP; afterwards, the reconstruction performance begins to decrease. This is different from the results of INK-SVD and IPR algorithms whose performances monotonically decrease with the coherence (i.e., the incoherence of the dictionary is increasing). Hence, our algorithm increases the dictionary incoherence without the risk of loss in reconstruction accuracy. Furthermore, our method proves that an appropriate incoherent dictionary helps to improve the performance. However, one point should be noticed, incoherent dictionary learning algorithm with AcMIQP for sparse coding has the highest computing complexity comparing to the other meth-

Convergence analysis

In this part, we give the proof of convergence of the proposed algorithm for incoherent dictionary learning. As aforementioned, the proposed algorithm aims at tackling the constrained optimization problem by transforming the problem into an unconstrained optimization problem via the augmented Lagrangian method. In each iteration of the augmented Lagrangian method, the minimization problem with respect to the primal variables is solved by the EPALM algorithm. Thus, for proving the convergence of the algorithm, we need to prove the convergence of the augmented Lagrangian method and that of the EPALM.

Convergence of the augmented Lagrangian method

Before proceeding and for completeness, we give here the convergence of the augmented Lagrangian method [START_REF] Bertsekas | Nonlinear programming[END_REF]. Consider the general expression of an equality constrained problem:

min q(x) subject to p(x) = 0, ∀x ∈ X, (5.29)
where X is a closed set, and q and p are continuous functions in X.

Proposition 3 (Proposition 4.2.1 in [START_REF] Bertsekas | Nonlinear programming[END_REF]). Assume q and p are continuous functions, X is a closed set, and the constraint set {x ∈ Then, x * is the solution of the problem (5.30) when λ k here is λ * . Setting the first order derivative ∂(q + λ * p)(x) to zero, the solution x * can be reached.

X | p(x) = 0} is nonempty. For k = 0, 1, • • • , let x k

Convergence of the EPALM algorithm

The analysis of the convergence of the EPALM algorithm is based on the Kurdyka-Łojasiewicz (KL) equality. We begin by introducing its definition (i) φ(0) = 0;

(ii) φ is C 1 on (0, η);

(iii) for all s ∈ (0, η), φ (s) > 0;

(iv) for all x ∈ U ∩ [ f (x * ) < f < f (x * ) + η],
the Kurdyka-Łojasiewicz inequality holds: We now study the convergence property of the algorithm, that is, the convergence of the sequence generated by the proposed algorithm in this chapter. We will prove that, with a and b two fixed positive value, the proposed algorithm generates a sequence x k k∈N that satisfies the following conditions: H1. (Sufficient decrease condition). For each k ∈ N,

φ ( f (x) -f (x * ))d i st (0, ∂ f (x)) ≥ 1. ( 5
f (x k+1 ) + a x k+1 -x k 2 ≤ f (x k ); H2. (Relative error condition). For each k ∈ N, there exists w k+1 ∈ ∂ f (x k+1 ) such that w k+1 ≤ b x k+1 -x k ;

H3. (Continuity condition).

There exist a subsequence x k j j ∈N and x such that x k j → x and f (x k j ) → f ( x), when j → ∞.

(5.32)

Then the following theorem (Theorem 2.9 in [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: Proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF]) will be used to prove the convergence of the proposed algorithm. x k+1x k < +∞.

In the following, we begin with the proof of satisfaction of assumption on functions.

Proposition 4. The objective function

L (c i 1 ,c i 2 ) (D, G,λ λ λ i , H i ) is a KL function.
Proof. The objective function L (c i 1 ,c i 2 ) (D, G,λ λ λ i , H i ) can be written in form of (5.18), namely f (D) + h(D, G) + g (G). According to [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF] and therein, it is easy to prove that f and h are KL functions. Moreover, g is also a KL function because it is the indicator function of a semi-algebraic set. Hence, the sum of the KL functions, i.e.,

L (c i 1 ,c i 2 ) (D, G,λ λ λ i , H i ), is a KL function.
Proposition 5. In problem (5.29), if q(x) is a proper semicontinuous function in a closed set X and p(x) is a proper lower continuous function in X, then the augmented Lagrangian function L c k (x, λ k ) is a proper lower semicontinuous function.

Proof. Firstly, if p(x) is a continuous function in X, then {x | p(x) < ∞} = X. Moreover, q(x) is a proper function in X. L c k (x, λ k ) = q(x) + λ k p(x) + c k 2 p 2 (x) is consequently a proper function in X.

Secondly, it is evident that if p(x) is a continuous function in X, then λ k p(x) and c k 2 p 2 (x) are continuous functions in X. The sum of a semicontinuous function in X, the function q(x) and a continuous function λ k p(x) + c k 2 p 2 (x), is still a semicontinuous function, i.e., L c k (x, λ k ) is a semicontinuous function.

Finally, p and q are both lower-bounded functions in X, that is, ∀x ∈ X, p(x) > -∞ and q(x) > -∞. λ k p(x) + c k 2 p 2 (x) is a convex function because c k > 0, then, λ k p(x) + c k 2 p 2 (x) > -∞, ∀x ∈ X. Hence, L c k (x, λ k ) is the sum of two lower-bounded functions in X.

Therefore, L c k (x, λ k ) is a proper lower semicontinuous function.

In our optimization problem, q(x) is a proper lower semicontinuous function dedicating here to an indicator term, and p(x) represents all the constraints that can be linear functions and quadratic functions, and which are all proper lower continuous functions. By applying the above Proposition 5, we can deduce that the augmented Lagrangian function in our problem is a proper lower semicontinuous function. Now, to prove the convergence of the algorithm, we still need to prove that the generated sequence (D i , j , G i , j ) satisfies the conditions H1, H2 and H3. The se-quence is generated from the process D i , j +1 = arg min D∈R n×m P 1 (D) G i , j +1 = arg min G∈R m×m P 2 (G),

(5.33)

where the functions P 1 and P 2 are defined as:

     P 1 (D) = f (D i , j ) + h(D i , j , G i , j ) + tr (D -D i , j ) T ∇ D ( f (D i , j ) + h(D i , j , G i , j )) + t 1 2 D -D i , j 2 F P 2 (G) = g (G) + h(D i , j +1 , G i , j ) + tr (G -G i , j ) T ∇ G h(D i , j , G i , j ) + t 2 2 G -G i , j 2 F .
(5.34) Proposition 6. The process P 1 produces a sequence {D i , j } that respects the conditions H1, H2 and H3.

Proof. The three functions ∇ f : R n×m → R n×m , ∇ D h : R n×m → R n×m and ∇ G h :

R m×m → R m×m are all Lipchitz continuous functions on their own domain. Then, there exists a Lipchitz constant L 1 = L + L D , where L is the Lipchitz constant for the function ∇ f and L D defined in (5.19), that is f (D i , j +1 ) + h(D i , j +1 , G i , j ) ≤ f (D i , j ) + h(D i , j , G i , j ) + L 1 2 D i , j +1 -D i , j 2 F +tr (D i , j +1 -D ( i , j )) T ∇ D ( f (D i , j ) + h(D i , j , G i , j ) .

(5.35)

The minimization of the optimization problems (5.33) requires that tr (D -D i , j ) T ∇ D ( f (D i , j ) + h(D i , j , G i , j ) + t 1 2 D -D i , j 2 F ≤ 0, (5.36) which assures descent in the objective function. By combining inequality (5.36) and the inequality (5.35), we obtain the following result:

f (D i , j +1 ) + h(D i , j +1 , G i , j ) + t 1 -L 1 2 D i , j +1 -D i , j 2 F ≤ f (D i , j ) + h(D i , j , G i , j ).
(5.37)

The satisfaction of condition H1 can be easily proven by choosing a t 1 greater than the Lipchitz constant L 1 .

We now begin to prove the condition H2. A large-enough b can be found such applying the triangle inequality, the following is deduced:

∇ D f (D i , j +1 ) + ∇ D h(D i , j +1 , G i , j ) ≤ ∇ D f (D i , j ) -h(D i , j , G i , j )
+ ∇ D f (D i , j +1 ) + h(D i , j +1 , G i , j ) -∇ D f (D i , j )h(D i , j , G i , j ) ≤(b + L 1 ) D i , j +1 -D i , j , (5.39) which is the relative error condition H2.

The continuity condition H3 is satisfied because of the continuity of the functions f and h with respect to D.

Proposition 7. The process P 2 produces a sequence {G i , j } having the properties introduced in conditions H1, H2 and H3.

Proof. The minimization of the second subproblem in (5.33) assures that, g (G i , j +1 ) + tr (G i , j +1 -G i , j ) T ∇ G h(D i , j , G i , j ) + t 2 2 G i , j +1 -G i , j 2 F ≤ g (G i , j ). (5.40)

The function G → h(D, G) is a L D -Lipchitz continuous function. Here, for simplification, let L 2 = L G . Thus, the inequality (5.40) becomes g (G i , j +1 ) + -L 2 + t 2 2 G i , j +1 -G i , j 2 F ≤ g (G i , j ).

(5.41)

When t 2 > L 2 , the condition H1 is satisfied.

We prove the satisfaction of condition H2 by using its first order necessary condition:

∂g (G i , j +1 ) + t 2 (G i , j +1 -G i , j + 1 t 2 ∇ G h(D i , j +1 , G i , j )) = 0.

(5.42) Equivalently:

∂g (G i , j +1 ) + ∇ G h(D i , j +1 , G i , j ) = -t 2 (G i , j +1 -G i , j ).

(5.43)

Taking the norm on both sides, the following equality holds:

∂g (G i , j +1 ) + ∇ G h(D i , j +1 , G i , j ) = t 2 G i , j +1 -G i , j .

(5.44)

Then by applying the triangle inequality, the condition H2 can be proven: ∂g (G i , j +1 ) + ∇ G h(D i , j +1 , G i , j +1 ) ≤ ∂g (G i , j +1 ) + ∇ G h(D i , j +1 , G i , j ) + ∇ G h(D i , j +1 , G i , j +1 ) -∇ G h(D i , j +1 , G i , j ) ≤(t 2 + L 2 ) G i , j +1 -G i , j .

(5.45)

H3 is satisfied for the continuous function h and the semicontinuous function

g in S G .
Proposition 8. The iterative process P 1 and P 2 produces a sequence {(D i , j , G i , j )} that satisfies the conditions H1, H2 and H3.

Proof. The Lipchitz continuity of the gradient of G → h(D, G) and the inequality (5.45) infer that there exists an L < 0 that verifies h(D i , j +1 , G i , j +1 )h(D i , j +1 , G i , j ) ≤ L G i , j +1 -G i , j 2 F .

(5.46)

By summing the inequalities (5.37) and (5.41), we get: f (D i , j +1 ) + h(D i , j +1 , G i , j ) + g (G i , j +1 ) + t 1 -L 1 2 D i , j +1 -D i , j 2 F + t 2 -L 2 2 G i , j +1 -G i , j 2 F ≤ f (D i , j ) + h(D i , j , G i , j ) + g (G i , j ). (5.47)

Using the result of (5.46), the inequality becomes:

L c i 1 ,c i 2 (D i , j +1 , G i , j +1 ,λ λ λ i , H i ) + t 1 -L 1 2 D i , j +1 -D i , j 2 F + t 2 -L 2 -2L 2 G i , j +1 -G i , j 2 F ≤ L c i 1 ,c i 2 (D i , j , G i , j ,λ λ λ i , H i ).
(5.48)

Setting a = min t 1 -L 1 2 , t 2 -L 2 -2L 2 , we obtain

L c i 1 ,c i 2 (D i , j +1 , G i , j +1 ,λ λ λ i , H i ) + a (D i , j +1 , G i , j +1 ) -(D i , j , G i , j ) 2 F ≤ L c i 1 ,c i 2 (D i , j , G i , j ,λ λ λ i , H i ).
(5.49) Thus, the sequence (D i , j , G i , j ) j ∈N satisfies the condition H1.

To prove the condition H2, it is necessary to compute the subdifferential of the function (5.18) with respect to the pair of matrix variables (D i , j +1 , G i , j +1 ), which is denoted by ∂L (c i 1 ,c i 2 ) (D i , j +1 , G i , j +1 ,λ λ λ i , H i ). With the results obtained in (5.39) and

(5.45), we use again the triangle inequality, then

∂L (c i 1 ,c i
2 ) (D i , j +1 , G i , j +1 ,λ λ λ i , H i ) = ∇ f (D i , j +1 ) + ∇h(D i , j +1 , G i , j +1 ) + ∂g (G i , j +1 ) ≤ ∇ D f (D i , j +1 ) + ∇ D h(D i , j +1 , G i , j +1 )

+ ∂g (G i , j +1 ) + ∇ G h(D i , j +1 , G i , j +1 ) ≤ ∇ D f (D i , j +1 ) + ∇ D h(D i , j +1 , G i , j )

+ ∂g (G i , j +1 ) + ∇ G h(D i , j +1 , G i , j +1 ) + ∇ D h(D i , j +1 , G i , j +1 ) -∇ D h(D i , j +1 , G i , j ) .

Using the expressions of the partial derivatives in (5.24), then the following inequality holds ∇ D h(D i , j +1 , G i , j +1 ) -∇ D h(D i , j +1 , G i , j ) ≤ ∇h(D i , j +1 , G i , j +1 ) -∇h(D i , j +1 , G i , j )

≤ L (D i , j +1 , G i , j +1 ) -(D i , j , G i , j ) ,
where L is the Lipchitz constant of the function h. Combining the inequalities (5.39) and (5.44), the condition H2 of the global sequence (D i , j , G i , j ) j ∈N is obtained

∂L (c i 1 ,c i
2 ) (D i , j +1 , G i , j +1 ,λ λ λ i , H i ) ≤(L 1 + b) D i , j +1 -D i , j + (t 2 + L 2 ) G i , j +1 -G i , j + L (D i , j +1 , G i , j +1 ) -(D i , j , G i , j ) .

Let t = max(L 1 + b + L, L 2 + t 2 + L), then

∂L (c i 1 ,c i
2 ) (D i , j +1 , G i , j +1 ,λ λ λ i , H i ) ≤ t (D i , j +1 , G i , j +1 ) -(D i , j , G i , j ) .

The condition H3 is straightforward by considering the continuity of the function.

Conclusion

This section investigated the exact incoherent dictionary learning problem, where the constraint of coherence is explicitly added to formulate the 0 -norm based constrained optimisation problem. To solve this constrained optimization problem, we introduced the new dictionary update algorithm EPALM that combines the proximal alternating minimization method and augmented Lagrangian method are in-troduced. This algorithm was used for dictionary learning together with a sparse coding algorithm, such as the proximal method and AcMIQP. We showed firstly the feasibility of the algorithm on synthetic data, examining the performance of the dictionary learning independently of the sparse coding algorithm. And then, the incoherent dictionary algorithm was used for real image reconstruction. We studied the statistics of the resulting dictionary, and the reconstruction performance for a large set of target coherence parameters. It was proven that the combination of EPALM for dictionary updating and MIQP for sparse coding always outperformed the other methods in terms of the reconstruction results. The relevance of having an incoherent dictionary was also demonstrated.

Chapter 6

Conclusion and Future Work

Conclusion

This monograph focused on the optimization problem of sparse representation.

The problem was formulated as the minimization of the mean square error with a constraint on the sparsity of the decomposition coefficients. Generally, the sparsity of the decomposition coefficients can be promoted by the p -norm function for 0 ≤ p ≤ 1, specially, p = 0 or 1. The 0 -norm can explicitly count the number of non-zero coefficients, nevertheless with the characteristics of non-smoothness and non-convexity, which make the problem NP-hard. The nearest convex norm function to 0 -norm is 1 -norm, but is not differentiable at 0. Additionally, it can achieve the optimal solution under some conditions, such as the signal is sufficiently sparse [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. However, the resolution of the 1 -norm based optimization problem risks in missing the optimal sparse solution [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF]. Hence, the resolution of 0 -norm based sparse representation problems is of great interest, which was the main topic in this monograph.

The generic problem of sparse representation aims at estimating the decomposition coefficients with the predefined dictionary. This is the so-called sparse coding problem. Moreover, since the quality of the dictionary greatly affects the performance of the sparse representation, researchers have been increasingly interested in data-driven dictionary learning, instead of a predefined dictionary. Thus the dictionary learning problem emerged and began to play a significant role in sparse representation. The dictionary learning estimates jointly the sparse code and the dictionary. This problem is still NP-hard and much more difficult than the single problem of sparse coding. Nevertheless, this optimization problem can be solved by the strategy of iteratively alternating two steps: sparse coding and dictionary updating. In this monograph, we studied the problem of sparse coding and dictionary learning optimization problems based on the exact 0 -norm.

Chapter 3 presented the state-of-the-art methods for sparse coding and dictionary learning. For tackling the problem of sparse coding, three major classes of methods were proposed. The greedy methods, such as (Orthogonal) Matching Pursuit and its variants, solve the problem by iteratively finding the local optimal solution of the 0 -norm optimization problem. Another class is relaxation methods that replace the 0 -norm by the 1 -norm to make the problem convex and tractable, such as Basis Pursuit and LASSO. The last class brings together gradient descent methods and iterative thresholding methods, such as Iterative Soft Thresholding and Iterative

Hard Thresholding (IHT). However, these algorithms risk in missing the global optimal solution.

In Chapter 4, we reformulated the sparse coding problem by introducing an auxiliary variable to indicate if the corresponding coefficient is zero or not. Then, the subproblem of sparse coding was transformed to a mixed integer quadratic programming (MIQP). This method can exactly solve the problem by exploiting optimization strategies such as the branch-and-cut. The optimizer solver Gurobi was used for the implementation. Moreover, it is worthy noting that this algorithm had no assumption on the dictionary, which is unlike the IHT or greedy algorithms that need to satisfy conditions such as the RIC. Furthermore, two optimization techniques were proposed to accelerate the MIQP for sparse coding. The first one offered an appropriate initialization. The method of proximal method was used here for producing the initial estimation, since it is considered to be a good approximation of the optimal solution, and thus the searching time was greatly reduced. The with the coordinate descent algorithm based on SVD for dictionary updating, the resulting dictionary learning method proved its advantage in image recovery with high-level noise, compared to state-of-the-art methods K-SVD, proximal method, and SOUPDIL.

In Chapter 5, we studied the problem of incoherent dictionary learning, since incoherent dictionaries bring in great improvement in sparse representation. In general, this optimization problem can be tackled with two strategies. The first introduces an additive decorrelation step following the step of dictionary updating at each iteration, such as the algorithm INK-SVD. The other reformulates the classical dictionary learning problem by introducing a term of regularization that measures the dictionary coherence, such as the incoherent dictionary learning by proximal method. However, the additive decorrelation method did not show the improved performance with more incoherent dictionary, and the regularized formulation cannot explicitly measure the relationship of performance with the coherence of dictionary varying. In this monograph, we proposed a third strategy, formulating the dictionary learning problem with explicit constraints on the coherence. We proposed a novel method to solve the resulting optimization problem, using the strategy of iteratively alternating sparse coding and dictionary updating. The sparse coding was solved by applying the AcMIQP method. For dictionary updating, we faced an optimization problem with quadratic inequality constraints. For dealing with this problem, we firstly took advantage of Alternating Direction Method of Multipliers (ADMM). Then, with the method of augmented Lagrangian method, the optimization problem was transformed into a unconstrained optimization problem. However, it is a non-convex problem. To overcome this difficulty, the proximal alternating linearized minimization was exploited. This algorithm was used in the application of image reconstruction. It proved for the first time, to our best knowledge, the accuracy improvement of image recovery by appropriate reduction of coherence of the dictionary.

Future Work

The study of the exact 0 -norm optimization is a promising research topic in sparse representation. In its genesis with the work of Jokar and Pfetsch in [START_REF] Jokar | Exact and approximate sparse solutions of underdetermined linear equations[END_REF], the exact solution of the sparse coding problem was used to verify some propositions and to provide a criterion to evaluate the relevance of existing greedy and relaxation algorithms. More recently, the resolution of the exact 0 -norm sparse coding problem was studied by Bourguignon et al. in [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF], demonstrating that it can tackle a small-scaled problem with synthetic data. In this monograph, we proposed the algorithm AcMIQP which makes it possible to achieve the exact optimal solution and to be applied to real well-known signal and image processing. Besides sparse coding, dictionary learning was proposed with AcMIQP. Therefore, in the future, the exact 0 -norm could be used to deal with more complex problems. With the exact sparse results, the performance of the other algorithms could be evaluated from one more dimension. Hence, the theoretical conclusion could be verified more rigorously.

Moreover, we proposed a flexible model that is easy to extend. In this monograph, we extended the classical dictionary learning problem to the incoherent dictionary learning problem by adding the constraint of the dictionary. Similarly, for the task of classification, the constraint on the classifier may be added. For solving this problem, one needs just to add one more step to update the classifier following dictionary updating. Besides, if there is more demands on the decomposition coefficients or on the dictionary (e.g. the positivity or low-rank of the dictionary), the formulation of the problem is simple, namely, adding the corresponding constraints. Even for solving these novel problems, the strategy is to iteratively alternating the steps of finding the solution of one variable with fixing all the others.

For each subproblem, it may be a classical optimization problem or much more complex, which needs many steps to be transformed into a classical optimization problem. Then, the optimizer solver Cplex or Gurobi could be used to solve the problem, or one could exploit some methods that are less generic and more specific to the sparse representation task at hand.

As already pointed, the shortcoming of the exact 0 -norm method is its high computational complexity. Fortunately, we can benefit from the fact that the optimization theory is well developing and the computational capability of the machine is highly improving. Of particular interest, one may take advantage of the GPU and operate parallel computing. It can thus be predicted that the research on the acceleration of the exact 0 -norm can be an important tendency.

Finalement 1 Introduction"

 1 , nous conclurons nos travaux et donnerons des perspectives de notre recherche. vii viii Chapter Scientists must use the simplest means of arriving at their results and exclude everything not perceived by the senses. Outline of the monograph . . . . . . . . . . . . . . . . . . . . . . . . 7

Figure 1 . 1 -

 11 Figure 1.1 -Example of wavelet functions

  The 0 based sparse representation . . . . . . . . . . . . . . . . . . 12 2.1.1 Introduction of 0 -norm function . . . . . . . . . . . . . . . . 12 2.1.2 The 0 -norm function promoting sparsity . . . . . . . . . . . 14 2.1.3 Some approximations of the 0 -norm function . . . . . . . . 16 2.1.4 Exact reformulation with complementary constraints . . . . 18 2.1.5 Existence and uniqueness of the solution . . . . . . . . . . . . 18 2.2 The 1 based sparse representation . . . . . . . . . . . . . . . . . . 19 2.3 The p based sparse representation . . . . . . . . . . . . . . . . . . 21 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Figure 2 . 1 -

 21 Figure 2.1 -Illustration of α α α p p in one dimension with p = 0, 0.1, 0.5, 0.8, 1, 2

  Despite the difficulty, a large number of researchers paid attention to the 0 based optimization problem [Nik13, AGH18, XZJ13]. The developments have been made both in theory and applications. Tropp analyzed greedy algorithms to solve 0 based sparse representation and summarized the sufficient condition for obtaining an optimal solution [Tro04]. Soubies et al. provided the sufficient and necessary condition for continuous function approximating 0 -norm [SBFA17]. The 0 norm is, thus, utilized in image processing [XZJ13, EA06, MP06], signal processing [Mal99, LPM00], machine learning [BJQS14], etc. 2.1.2 The 0 -norm function promoting sparsity By setting the function Ω(α α α) = α α α 0 , the optimization problem (2.1) becomes: arg min α α α α α α 0 subject to x = Dα α α.

  Figure 2.2 -Illustration in 2D (i.e., α α α ∈ R 2 ) of the minimization of the quadratic reconstruction error (parabola in blue) in the feasible region defined by the p -norm constraint ( p -ball region in red) with sparsity level k = 1.

  r i = x -D T i α α α T i ; number of iterations depends on the properties of the signal. Algorithm 3 limits the iteration number to halt the loop. Besides, other stopping criteria have been presented, such as the reconstruct error [NT09, Apprendix].

LASSO provides a new

  point of view to solve the 1 -norm constrained square error minimization problem (3.5), or its regularized problem (3.7)[START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. The method of LASSO is realized by introducing an auxiliary matrix G whose rows contain the signs of the corresponding α α α. Consequently, the inequality constraint can be represented by the linear inequality Gα α α ≤ λ k 1 1 1. With considering the sign of a real value being either +1 or -1, there are totally 2 p kinds of possible structures. The number of rows of G is thus 2 p . By eliminating calculations on G of too big size, an iterative process Algorithm 4 LASSO Input: the signal x, the dictionary D, the sparsity λ k Output: The sparse representation α α α.1: function SPARSE CODING 2:

  04].The homotopy algorithms[START_REF] Dmitry M Malioutov | Homotopy continuation for sparse signal representation[END_REF][START_REF] Garrigues | An homotopy algorithm for the lasso with online observations[END_REF] are developed based on the fact that the solution to the regularized problem can be computed directly by letting the zero vector belong to the subdifferential of the objective function, which was presented in section 2 of [GG09], namely it exists a d α d α d α that satisfies D T (Dα α αx) + λd α d α d α = 0, (3.15) where d α d α d α denotes the subdifferential of the objective function of the regularized problem, which takes the form

  cessively added or removed by getting the critical value of λ, namely the value that makes a non-zero entries of α α α turn to zero and the subdifferential corresponding to zero entries of α α α reaching the limit d α d α d α ∞ = 1. Thereupon, by decreasing λ from D T x ∞ to zero, we get the optimal solution corresponding to each λ. Choosing the solution with the desirable number of zeros in α results in the final solution. The searching path is given in [GG09, MY12, PH07]. Garrigues and Ghaoui [GG09] im-

  rogate function [BYD07, BD08, DDDM04]. The surrogate function is obtained from reformulating the original problem by introducing an auxiliary variable a ∈ R m as follows arg min

  where v = sign(∇ ω ĵ f (ω ω ω)) and γ ∈ [0, 1]. It is noticed that the direction of descent is v -ω ω ω where v is the vector whose ĵ t h component is v, namely v ∈ S. With assumption of ω ω ω i ∈ S and γ varying from 0 to 1, each update is still in the set S. The descent step γ is determined by minimizing γ = arg min γ∈[0,1]

ordinateF=

  Descent algorithm based on Singular Value Decomposition (SVD) updates successively the atoms of the dictionary. In the following, we introduce the algorithm proposed in [AEB + 06]. Let B = [b 1 , . . . , b m ] ∈ R l ×m with B = A T , namely b j denotes the j t h row of A. The optimization problem with respect to D is then formulated by D = arg min D∈D X -DA 2 arg min D∈D (Xj =i d j b T j )d i b T

  which has been extended to other tasks. Task-specified dictionaries are designed by reformulating the problem with a special regularization term. For instance, discriminative dictionaries are learned by introducing penalty on the classification error (e.g. measured by linear classifier [JLD13, ZL10], by logistic loss [MPS + 09], or by Fisher criterion [HA07]). Low-rank dictionaries are obtained by learned problem with a regularization of the nuclear norm [ZJD13]. Incoherent dictionary learning is realized by adding a regularization term about the difference between D T D and the identity matrix of size m × m, denoted by I m [MBP12, TLZ + 19, AFS15, BQJ14] (this problem will be discussed in detail in Chapter 5 of this monograph).

Algorithm 7

 7 Algorithm of K-SVDInput: Signals set X, sparsity level k, maximal number of iteration N. Output: The sparse representation A and dictionary D.1: Initialization The dictionary D 0 .2: function DICTIONARY LEARNING 3:

m

  = the number of column of D 0 . 8: for j = 1 to m do 9:

  random selection rule. Due to the great achievements of the K-SVD algorithm [Sch14, RZE08, MBP12, ZL10], some variants are developed on the algorithm. For example, the kernel K-SVD provides a nonlinear variant of K-SVD [GE16, VNPNC12, ZSJ + 16, KD16]. The INK-SVD (for INcoherent K-SVD) allows to delete atoms by adding a decorrelation step following the dictionary update [MBP12] (See Chapter 5 for details).

  ) where B 0 corresponds to the number of non-zero elements in matrix B. Then applying the formulation of sum of outer products DB T = i d i b T i , noting that the problem restricts the total number of non-zero elements in all signal representations, in contrast with problem (3.25) that constrains only the sparsity of each signal representation. Thus, this formulation allows a flexible sparsity for different signals.For solving the optimization problem (3.29), a block coordinate descent method is applied. In SOUPDIL, the number of inner iterations is m, which indicates the row size of dictionary D or that of matrix B. Thereby, in each iteration i , the i t h column of B, i.e., b i , and the i t h column of D, i.e., d i , are successively updated. Hereafter, we focus on the algorithm details in each iteration.For sparse coding, consider the problem of estimating a column b i of B, with D and all b j with j = i being fixed.Let E i = Xj =i , d j b T j and B = {v ∈ R l | v ∞ < L}be the set of admissible solutions of b i , where L is defined in[START_REF] Ravishankar | Efficient sum of outer products dictionary learning (soup-dil) and its application to inverse problems[END_REF] to avoid the Algorithm 8 Sum of OUter Product Dictionary Learning (SOUPDIL)Input: the input signals X, parameter λ, the upper bound L, maximal number of outer iterations N, the unit-norm vector v with the first entry one and the rest zeros. Output: The transpose of sparse representation B and dictionary D.1: function DICTIONARY LEARNING 2:

  algorithm for exact 0 penalty problem . . . . . . . . . 54 4.1.1 Mixed integer programming (MIP) . . . . . . . . . . . . . . . . 55 4.1.2 Accelerated mixed integer programming (AcMIQP) . . . . . . 63 4.1.3 Dictionary learning with AcMIQP for sparse coding . . . . . . 69 4.2 Image denoising based on the exact 0 . . . . . . . . . . . . . . . . . 71 4.2.1 Sparse representation for image denoising . . . . . . . . . . . 71 4.2.2 Large-scale dictionary learning . . . . . . . . . . . . . . . . . 75 4.2.3 Adapted dictionary learning . . . . . . . . . . . . . . . . . . . 78 4.2.4 Analysis on the computational complexity . . . . . . . . . . . 81 4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

  0 penalty problem This section focuses on the problem (3.26) with Ω the 0 -norm function. We recall the problem that, given a set of l signals X = [x 1 , . . . , x ] ∈ R n×l of length n and a given sparsity level k, sparse representation aims to find the optimal a linear combination over an overcomplete dictionary D = [d 1 , . . . , d m ] ∈ D where D = {D ∈ R n×m | d i 2 = 1, ∀i = 1, . . . , m} and sparse coefficient matrix A = [α α α 1 , . . . ,α α α ] ∈ R m×l , namely by solving the problem arg min A∈R m×l ,D∈D 1 2 X -DA 2 F , subject to α α α i 0 ≤ k ∀i = 1, . . . , l . (4.1) This optimization problem of estimating simultaneously the two variables A and D is NP-hard. A strategy is to solve A and D alternatingly. When fixing A, D can be effectively obtained by dictionary updating algorithms such as MOD. When fixing D, the optimization problem with respect to A is the sparse coding problem, which has convex objective functions but with a non-convex 0 -norm constraint that makes the optimization problem NP-hard. For the resolution of this problem, methods that rely on a greedy algorithm [MZ93, PRK93, NT09], may produce a local optimal [Tro04, BJQS14], while relaxation formulations, such as with LASSO may sometimes achieve a solution not as sparse as that of the 0 constrained problem.

  where Q ∈ R n×n is a symmetric matrix and the vector c ∈ R n corresponds to the linear part of the objective function. The (4.3) ∼ (4.6) are constraints. Specifically, (4.3) is a linear inequality constraint where A i n is a matrix of n columns with its rows depending on the number of linear inequalities, b i n and b i are the right sides of the linear constraint and quadratic constraint respectively, and (4.4) defines the bounds for variable v where vectors lb ∈ R n and ub ∈ R n are respectively given lower bound and upper bound. The constraints determined by the two bounds are fundamental in the complexity sense; without bound constraints, the problem becomes undecidable[START_REF] Robert | There cannot be any algorithm for integer programming with quadratic constraints[END_REF]. In contrast, by raising the lower bounds and reducing the upper bounds, the computation complexity can be easily decreased[START_REF] Neumaier | Safe bounds in linear and mixed-integer linear programming[END_REF]. For (4.5), it represents the quadratic constraint where Q i ∈ R n×n and c i ∈ R n . If Q i = 0, the quadratic constraint (4.5) degenerates into a linear constraint. The last constraint (4.6) defines the type of each element in v, namely I denotes a proper subset of {1, . . . , n} that identifies which entries of v are restricted to be integers.

  (4.12) where the symmetric matrix Q, defining quadratic part of the objective function, of size 2m × 2m is made up of four sub-matrices, namely Q = D T D 0 m,m 0 m,m 0 m,m , with 0 m,m the zero matrix of size m × m, and c, defining the linear part in the objective function, is a column vector of size 2m that is defined by c = -D T x 0 m,1 . Now, consider the inequality constraint. The matrix of size (2m+1)×2m in the linear inequality constraint is of the form

  4.8) by the big-M reformulation, the indicator function was recently proposed to deal with the logical relation problem, as explained next [BBF + 16, BLTW15]. The logical relation (4.8) can be reformulated by introducing the indicator function as

[

  NM65], dual-simplex method, interior point method[START_REF] Nocedal | Numerical optimization[END_REF] or barrier method [GMS + 86]. An effective method is the branch-and-bound technique [BM91, BM97, BE07]. The branch-and-bound method processes the MIP problem by solving a sequence of relaxed LP problems. These relaxation subproblems are organized in a tree structure. Each node of the tree contains a subproblem. Then the solution searching process [VAN08, BM91] can be summarized as: 1) Initialization: Setting the parent node as an LP problem by relaxing all the integer variables to become continuous. The initial upper bound ub is set to

  often chosen to conduct the experiments). Then by vectorization of all the pixels in the patches, the matrix of signals Y = [y 1 , . . . , y i , . . . , y ] is produced. Then, taking Y as the input signals, the optimal sparse representation A and dictionary D can be learned as mentioned in Section 3.4. Each underlying patch can be recovered by X = DA. (4.25)

Figure 4 . 1 -

 41 Figure 4.1 -Examples in the USC-SIPI Image Database

Figure 4

 4 Figure 4.2 -Convergence of the proposed algorithm and the comparison with K-SVD using OMP, proximal method and SOUPDIL

Figure 4 .

 4 Figure 4.3 provides the histogram of the correlation between each pair of atoms of the learned dictionaries, for each of the four methods under investigation. It is observed that the distribution of the correlations of the obtained dictionaries can be roughly ordered as follows

Figure 4 . 3 -

 43 Figure 4.3 -Histogram of the correlation of the learned dictionaries

here 4 .

 4 For the other three methods, the experiments details and the parameter settings are given in the following. When dealing with noisy data in the training phase, the knowledge about the noise level σ is used for restricting the reconstruction error, as shown in the constraint in the optimization problem (4.23) and the parameter setting = c n σ 2 with c = 1.15. These values, optimized for OMP in[START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], are used here for both proximal method and AcMIQP, thus putting our method in a less favorable situation.

  Coherence and sparse representation . . . . . . . . . . . . . . . . . 86 5.1.1 Definition of coherence . . . . . . . . . . . . . . . . . . . . . . 87 5.1.2 Coherence of dictionary in sparse representations . . . . . . 88 5.2 Methods of learning incoherent dictionary . . . . . . . . . . . . . . 89 5.2.1 Dictionary learning with additive decorrelation . . . . . . . . 90 5.2.2 Regularized incoherent dictionary learning . . . . . . . . . . 92 5.2.3 Towards incoherence-constrained dictionary learning . . . . 94 5.3 Exact 0 based incoherent dictionary learning . . . . . . . . . . . . 95 5.3.1 Augmented Lagrangian method . . . . . . . . . . . . . . . . . 96 5.3.2 Proximal alternating linearized minimization . . . . . . . . . 97 5.3.3 Experiments to assess the dictionary updating algorithm . . 5.3.4 Incoherent dictionary learning algorithm . . . . . . . . . . . . 5.3.5 Convergence analysis . . . . . . . . . . . . . . . . . . . . . . . 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  , x ] ∈ R n× of signals of dimension n, aims at finding a sparse representation of the coefficient matrix A for a decomposition of the form X = DA, where the matrix A = [α α α 1 , . . . ,α α α ] ∈ R m× containing the decomposition coefficients is sparse, and the matrix D = [d 1 , . . . , d m ] ∈ R n×m is the dictionary with each column called atom. The optimization problem is written as min D∈C,α α

  tation (IPR) [MBP12, DM13]. 2) The second strategy learns an incoherent dictionary by introducing the regularization term D T D -I m 2 F , where I m is the identity matrix of size m × m [RLS09, BQJ14]. The normalization is realized by adding the regularization term

[ DM13 ]

 DM13 , optimizes the INK-SVD algorithm by considering simultaneously the minimization of the residual error of sparse approximation when learning the dictionary Algorithm 11 Algorithm of decorrelation in INK-SVD [MBP12] Input: Initial dictionary D and required coherence µ c . Output: Approximate dictionary D with maximal inner product of two atoms equal to µ c . 1: while max i = j |〈d i , d j 〉| > do 4: Selection of a pair of atoms (d i , d j ) = arg max |D T D -I m |; 5: Decorrelation of the two atoms by rotating them until angle between them reaches θ with cos(θ) = µ c . with a fixed target coherence level, that is, arg min D∈D X -DA 2 F subject to µ(D) ≤ µ c .

(

  WD) T (WD) = G. Consequently, the reconstruction error minimizing problem is turned into an optimization problem with respect to W, that is, W = arg min W∈N X -WDA 2 F , where N denotes the set of orthogonal matrices of size m × m. This problem has a closed-form solution W = SU T , where U and S are the left and right unitary matrices in the singular value decomposition (SVD) of DAX T , namely DAX T = UΣS T .

G

  2, . . . , m}, p = q d T p d p = 1, p = 1, . . . , m. (5.12) By introducing a new variable G ∈ R m×m that satisfies the identity G = D T D, the problem can be written in the form 2 min D∈R n×m ,G∈R m×m = D T D, G ∈ S G d T p d p = 1, p = 1, . . . , m, (5.13)

  ) where D i , G i λ λ λ i and H i are the values in the i -th iteration; In Step 2, the Lagrangian multipliers λ λ λ and H are updated; And in Step 3, the penalty parameters c 1 and c 2 are increased. It is proven that the two parameters c 1 and c 2 can stay much smaller to solve such optimization problem[START_REF] Bertsekas | Nonlinear programming[END_REF].

  ) is a non-convex and non-smooth optimization problem. It is unsolvable by satisfying the Karush-Kuhn-Tucker (KKT) conditions. We propose to use an alternating strategy to address this optimization problem. The optimal ma-The problem is well defined with f : R n×m → (-∞, +∞], h : R n×m × R m×m → (-∞, +∞] being a C 1 function (i.e., continuously differentiable), and inf f (D) > -∞, inf h(D, G) > -∞, g : R p×p → [0, ∞] a proper and lower semicontinuous function. In addition, for guaranteeing the convergence of the EPALM method, the following assumptions should be satisfied: (i) The functions D → h(D, G) and G → h(D, G) have their gradients globally Lipschitz continuous with moduli L D and L G , respectively. In other words, the partial gradients of h with respect to D and G verify the property:

Figure 5 .Figure 5

 55 Figure 5.1 presents the convergence property of the algorithms with the coherent parameter set to µ c = 0.6, which corresponds to having angles between any two atoms greater than 53 • . It is observable that the dictionary algorithm with MIQP for

Figure 5

 5 Figure 5.2 -The distribution of the coherence between each two atoms of the proposed algorithms and its comparison to the INK-SVD and IPR algorithms

Figure 5 .

 5 Figure 5.4 presents the best results of image reconstruction obtained respectively by algorithms of INK-SVD, IPR, prox+EPALM and AcMIQP+EPALM by setting the coherence parameter respectively to µ c = 0.996, 0.996, 0.707, 0.866. Patches are recovered by X = DA and then by averaging pixel value, the final reconstructed image is obtained. It is observable that the combination of AcMIQP for sparse coding and the proposed algorithm for dictionary update outperforms the other methods by maintaining most of the details in the image.

  Figure 5.3 -The reconstruction errors in PSNR for each method at each target coherence value

Figure 5

 5 Figure 5.4 -Reconstructed images produced respectively by algorithms (a) INK-SVD, (b) IPR, (c) Prox+EPALM and (d) AcMIQP+EPALM

  [ABS13]: Definition 2. (Kurdyka-Łojasiewicz function) (a) The function f : R n → R ∪ {+∞} is said to have the Kurdyka-Łojasiewicz property at x * ∈ dom ∂ f if there exist η ∈ (0, +∞], a neighborhood U of x * and a continuous concave function φ : [0, η) → R + such that:

  .31) (b) The proper lower semicontinuous functions that satisfy the Kurdyka-Łojasiewicz inequality at each point of dom ∂ f are called KL functions.

Theorem 5.3. 1 .

 1 (Convergence to a critical point) Let f : R n → R ∪ {+∞} be a proper lower semicontinuous function. Consider a sequence x k k∈N that satisfies the conditions H1, H2, H3. If f has the Kurdyka-Łojasiewicz property at the cluster point x specified in H3, then the sequence x k k∈N converges to x, and x → x as k goes to infinity, and x is a critical point of f . Moreover the sequence x k k∈N has a finite length, i.e., +∞ k=0

  that,∇ D ( f (D i , j ) + h(D i , j , G i , j ) ≤ b D k+1 -D k+1 .(5.38)By considering the Lipchitz continuity of the function D → ∇( f (D) + h(D, G)) and

  second technique was the method of relaxation with additive constraints that consider the convex envelop of the continuous variables. With these two acceleration techniques, the proposed Accelerated MIQP (AcMIQP) was investigated for sparse coding and further for handling the dictionary learning problem. AcMIQP allowed to break the limitation of previous work where the MIQP was only applied on smallscale synthetic data. The proposed algorithm can be used to deal with classical well-known data in signal and image processing. In our experiments, AcMIQP was used in image reconstruction and denoising. More precisely, when dealing with different formulations of the sparse representation problem, experiments on AcMIQP showed that the sparsity constrained based model reached the best results compar-ing to the error constrained one. Besides, by combining AcMIQP for sparse coding

integer pro- gramming for sparse coding: Application to image denoising

	national Workshop on Machine Learning for Signal Processing (MLSP), Aal-
	borg, 2018
	• Liu, Yuan, Canu, Stéphane, Honeine, Paul and Ruan, Su, Apprentissage de
	dictionnaire faiblement cohérent par programmation quadratique mixte.
	Actes du 27ème colloque du GRETSI en traitement du signal et des images,
	Lille, 2019.

• Liu, Yuan, Canu, Stéphane, Honeine, Paul and Ruan, Su, Mixed . IEEE Transactions on Computational Imaging, vol. 5, no. 3, pp. 354-365, Sept. 2019. • Liu, Yuan, Canu, Stéphane, Honeine, Paul and Ruan, Su, Incoherent dictionary learning via mixed-integer programming and hybrid augmented lagrangian. Submitted to Digital Signal Processing, Sept. 2019. Conference Papers • Liu, Yuan, Canu, Stéphane, Honeine, Paul and Ruan, Su, Une véritable approche 0 pour l'apprentissage de dictionnaire. Actes du 26ème colloque du GRETSI en traitement du signal et des images, Juan Les Pins, 2017. • Liu, Yuan, Canu, Stéphane, Honeine, Paul and Ruan, Su, K-SVD with a real 0 optimization: Application to image denoising. Proc. of the 28th IEEE Inter-

.4 Exact reformulation of 0 -norm with complementary con- straints

	2.1Very recently, a new strategy [SIBD11, FMP + 18] inspired by the work of [KM82] is
	emerging, where the 0 -norm is reformulated by introducing complementary con-
	straints. Specifically, this strategy consists in indicating if an entry of the vector α α α is
	zero or not by an auxiliary variable z ∈ R m with the definition
	z i = 1 ⇐⇒ α i = 0 ∀i = 1, . . . , m	(2.14)
	where z i means the i t h entry of z. Feng et al. presented two formulations in
	[FMP + 18]. The full-complementary formulation is
		.13)
	This formulation retains the equivalence of the global minimizer with the original
	problem. Moreover, the equivalence between critical minimizers of the two formu-
	lations is proved in Theorems 4.5, 4.7 and 4.8 in [SBFA17].

Theorem 2.1.1. Uniqueness-Spark (Theorem 2.4 [Ela10]) If there

	The uniqueness of the solution depends on a concept named 'Spark' of the dic-
	tionary, which is defined as the smallest number of columns from D that are linearly
	dependent. Theorem 2.4 in [Ela10] (or Theorem 1.2 in [DH01]) states the unique-
	ness of the solution to the equality constrained sparse model:
	is a solution α α α to
	problem (2.1) with α α α 0 < Spark(D)/2, then α α α must be the sparsest one.

Chapter 3 Sparse Representation Algorithms Sommaire 3.1 Sparse coding and dictionary learning . . . . . . . . . . . . . . . . 25 3.2 Algorithms for sparse coding . . . . . . . . . . . . . . . . . . . . . . 26

	3.2.1 Greedy algorithms . . . . . . . . . . . . . . . . . . . . . . . . . 27
	3.2.2 p -norm relaxation . . . . . . . . . . . . . . . . . . . . . . . . . 34
	3.2.3 Gradient descent and iterative thresholding algorithms . . . 36
	3.3

Algorithms for dictionary updating . . . . . . . . . . . . . . . . . . 42

	3.3.1 The method of optimal directions (MOD) . . . . . . . . . . . 42
	3.3.2 The algorithm of coordinate descent based on SVD . . . . . . 43

3.4 Algorithms for dictionary learning . . . . . . . . . . . . . . . . . . 44

	3.4.1 K-SVD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
	3.4.2 Sum of outer products dictionary learning . . . . . . . . . . . 46
	3.4.3 Proximal method . . . . . . . . . . . . . . . . . . . . . . . . . . 48
	3.

5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

  , Donoho et al. proposed in [DDTS06] the algorithm Stagewise Orthogonal Matching Pursuit (StOMP). In this algorithm, numerous atoms are determined at each stage by thresholding. The threshold b is defined as a function of the noise level σ, that is,

b = pσ,

where p is a predefined parameter. With adequate threshold setting, StOMP is proved to reach the exact sparse signal recovery as OMP but with higher speed. The algorithm is also suitable to deal with noisy signals. However, in this case, the results are not the one get by OMP.

OMP and its variants have been frequently used for sparse coding for dictionary learning [AEB + 06, MBP12, Sch14]. Likewise, algorithms applying the same spirit of OMP are used in tasks such as classification [RS08, HA07, ZL10, MLB + 08], face recognition

[START_REF] Zhang | Discriminative K-SVD for dictionary learning in face recognition[END_REF]

and image denoising

[START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]

.

  Iterative Hard Thresholding Algorithm for 0 Regularization Input: the signal x, the dictionary D, the sparsity level k, the stepsize γ, the stopping criterion .

	often take advantage of threshold-ing gradient descent, which is used to solve sparse representation problems due to its fast convergence and theoretical sparse approximation guarantee [Fou11, YLZ18, GK09, KKK08]. It is noticed that Iterative Thresholding Algorithm can be adopted to 1: function SPARSE CODING solve Output: The sparse representation α α α. 2:

0 -norm based optimization problems [BYD07, GK09, YLZ14], and 1 -norm based optimization problems [FISI16, Mei07, WNF09, KXAH15]. Algorithm 5

  Commonly, MIP refers to optimization problems with linear objective function and under linear constraints, which is thus specially called Mixed Integer Linear Programming (MILP). When the optimization problem has a quadratic objective function, it is called Mixed Integer Quadratic Programming (MIQP). When even the constraints are quadratic, it is given the name Mixed Integer Quadratically Constrained Programming (MIQCP)

Table 4 .

 4 1 -Computational time and accuracy results (100-batch median, 5th and 95th percentiles) on synthetic data[START_REF] Liu | Mixed integer programming for sparse coding: Application to image denoising[END_REF].

			Results Method	OMP	Proximal	MIQP	AcMIQP
	Computational time	P 5 median	0.019 0.020	0.023 0.024	2002.5 2002.6	286.98 415.34
			P 95	0.032	0.032	2002.7	543.90
	Reconstruction	P 5	14.00	13.20	2.36	1.58
	error		median	15.90	15.14	2.84	2.74
	X i -D Âi	2 F	P 95	17.45	16.19	2.86	3.51
	Sparse coding	P 5	39.17	34.68	5.41	3.60
	error		median	44.84	42.33	6.45	6.06
	Âi -A i F	P 95	51.13	46.40	6.46	7.99
	Position accuracy	P 5	98.16	98.42	99.73	99.94
	of non-zero	median	98.52	98.62	99.74	99.97
	elements (%)	P 95	98.81	98.98	99.78	99.98

Table 4 .

 4 2 -Accuracy of the denoising in terms of the PSNR in the large-scale (global) dictionary learning, for each of the five images at several noise levels, comparing the sparsityconstrained formulation (4.23) and the error-constrained formulation (4.21)(the higher, the better)

	Image	Sparse coding formulation σ =10 σ =20 σ = 50
	Barbara	error-constrained	24.71	23.75	20.79
		sparsity-constrained	26.77	25.24	20.14
	Cameraman	error-constrained	24.93	23.90	20.16
		sparsity-constrained	27.70	25.75	20.19
	Elaine	error-constrained	26.78	25.64	21.57
		sparsity-constrained	29.87	27.81	21.14
	Lena	error-constrained	26.05	24.98	21.22
		sparsity-constrained	28.83	26.93	20.92
	Man	error-constrained	24.67	23.68	20.80
		sparsity-constrained	27.60	25.97	20.10
	Average	error-constrained	25.43	24.39	20.91
		sparsity-constrained	28.15	21.15	20.49

Table 4 .

 4 3 -Accuracy of the reconstruction in terms of the PSNR (the higher, the better)

	Image	Method	Large-scale	Large-scale
			dictionary learning,	dictionary learning,
			reconstruction with	reconstruction with
			(4.25)	(4.27)
		OMP	19.71	20.03
	Barbara	Prox	20.71	21.03
		AcMIQP	22.73	22.73
		OMP	19.46	19.78
	Cameraman	Prox	21.11	21.43
		AcMIQP	22.30	22.62
		OMP	19.73	20.05
	Elaine	Prox	22.87	23.19
		AcMIQP	24.20	24.52
		OMP	19.79	20.11
	Lena	Prox	22.12	22.44
		AcMIQP	24.20	24.52
		OMP	19.68	20.00
	Man	Prox	21.26	21.58
		AcMIQP	23.59	23.91
		OMP	19.674	19.994
	Average	Prox	21.614	21.934
		AcMIQP	23.404	23.724
	implemented using the original Matlab code provided by its authors and available

Table 4 .

 4 4 -Denoising results in the adapted dictionary learning setting, for each of the five images, as well as the average results (the higher, the better)

	Image	Method	PSNR
		OMP	22.04
		OMP (pre-centering) 21.97
	Barbara	proximal	22.54
		SOUPDIL	22.30
		AcMIQP	22.59
		OMP	22.54
		OMP (pre-centering) 22.63
	Cameramen	proximal	22.49
		SOUPDIL	22.79
		AcMIQP	22.58
		OMP	23.00
		OMP (pre-centering) 22.91
	Elaine	proximal	23.29
		SOUPDIL	23.43
		AcMIQP	23.39
		OMP	22.48
		OMP (pre-centering) 22.51
	Lena	proximal	23.08
		SOUPDIL	23.20
		AcMIQP	23.09
		OMP	21.23
		OMP (pre-centering) 21.32
	Man	proximal	21.70
		SOUPDIL	21.67
		MIQP	21.86

Table 5

 5

	.1 -Accuracy results and computing time on synthetic data
	Coherence parameter µ c	0.5	0.55	0.6	0.7	0.8	0.9	1.0
	Initial objective function
	value: 1 2 X -D 0 A 2 F		9523	9592	9318	9643 9340 9483 9446
	Final objective function
	value: 1 2 X -DA 2 F		216.84 91.80	2.01	1.89	1.89	1.89	1.89
	Accuracy: D -D *	F	1.480	1.060	0.058 0.100 0.100 0.100 0.100
	Outer iteration number:
	N i t er in Algorithm 12	50	22	9	4	4	4	4
	Inner iteration number:
	N i i t er in Algorithm 13	972	877	692	312	319	319	318
	max |d T p d p -1|		0.015	0.001	0.003 0.007 0.007 0.007 0.007
	max |G -D T D|		0.052 0.0057 0.0084 10 -8 10 -8 10 -8 10 -8
	Computing time		155.78 47.31	16.33	3.38	3.42	3.41	3.39

Table 5 .

 5 2 -Statistics on the resulting dictionary

		µ	Average of	Variance of
			{|d T i d j | | i = j } {|d T i d j | | i = j }
	INK-SVD	0.601	0.368	0.0177
	IPR	0.711	0.557	0.0073
	Prox+EPALM	0.608	0.352	0.0176
	AcMIQP+EPALM 0.609	0.382	0.0146
	is bigger than θ

c = {5 • , 15

  • , 30 • , 45• , 60• , 75 • , 83 • }. Table 5.3 and Figure 5.3 present the reconstruction errors for each of these target values. For the incoherent dictionary

Table 5 .

 5 3 -The reconstruction errors in PSNR (in dB) by using the dictionary with different coherence parameter values µ c = cos(θ c ) It is remarked that there is no result of INK-SVD or IPR when θ c = 83 • . That is caused by the non-convergence of the two algorithms when dealing with some high coherence values. By contrast, the hybrid algorithm of ADMM and EPALM presents good performance in stability.

			Largest angle θ c between two atoms
		5 •	15 •	30 •	45 •	60 •	75 •	83 •
	INK-SVD	36.46 36.56 36.26 34,83 34.04 30.15	-
	IPR	36.82 36.57 35.72 31.51 30.60 27.84	-
	Prox+EPALM	27.40 27.42 28.06 29.80 29.31 29.75 22.97
	AcMIQP+EPALM 37.60 37.26 38.89 38.55 36.52 35.31 33.97
	ods.

, ∀x ∈ R m and x 0 ≤ k.If the RIC δ k exists, we say that the matrix D satisfies the k-restricted isometry property.

-norm can reach the same solution obtained by exploiting 0 -norm. Therefore, it is frequently used as a relaxation of the 0 -norm.

https://www.ibm.com/support/knowledgecenter/SSSA5P_12.9.0/ilog.odms.cplex.help/CPLEX/ UsrMan/topics/discr_optim/mip_quadratic/02_introMIQP.html

https://www.gurobi.com/products/gurobi-optimizer

http://sipi.usc.edu/database/database.php?volume=misc

https://gitlab.eecs.umich.edu/fessler/soupdil_dinokat

http://www.cs.technion.ac.il/ ronrubin/software.html

Another formulation can be proposed, by putting the unit-norm constraint in the space definition[START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF][START_REF] Yang | Linearized augmented Lagrangian and alternating direction methods for nuclear norm minimization[END_REF]. However, this formulation suffers from low convergence. For this reason, we consider the optimization problem (5.13) where the explicit unit-norm constraint will be addressed using the augmented Lagrangian method.

Remerciements

Algorithm 12

The inexact ADMM framework for solving (5.16)

Input:

The training data X and coefficient matrix A, the stop criteria ( , N i t er ). Output: The optimal solution D * function DICTIONARYUPDATING Initialization the parameters (λ λ λ 0 , D 0 , c 0 1 , H 0 , G 0 , c 0 2 , ρ 1 , ρ 2 ) for all i = 0 to N i t er -1 do 1. Computing the optimal solution (D i , G i ):

2. Updating the Lagrangian multiplier (λ λ λ i , H i ):

3. Updating the penalty parameters (c i 1 , c i 2 ):

trices of D and G are obtained by alternating the gradient descent method and the proximal method, which can be regarded as the special case of the extended proximal alternating linearized minimization (EPALM) [START_REF] Zhu | Nonconvex and nonsmooth optimization with generalized orthogonality constraints: An approximate augmented Lagrangian method[END_REF].

To investigate the EPALM method, we rewrite the objective function in problem (5.17) in the form of three additive parts:

with the definition of:

Algorithm 13 EPALM algorithm for solving subproblem (5.17)

Input:

The training data X and coefficient matrix A, the parameters (λ i λ i λ i , H i , c i 1 and c i 2 , the step size t 1 and t 2 ), the stop criteria ( i , N i i t er , the subdifferential

). while j < N i i t er and Θ i > i do 1. Updating D i , j by computing:

2. Computing:

3. Projecting the G in the space S G :

4. Calculating the subdifferential Θ i (D i , j , G i , j ). 5. j = j + 1.

end while end function

Proposition 2. To sum up, a sequence (D i , j , G i , j ), j ∈ N is generated by using the proposed method, then the following conditions will be satisfied:

(5.27)

Experiments to assess the dictionary updating algorithm

To assess the performance of the dictionary updating algorithm, it is tested on synthetic data. Specifically, a sparse matrix A ∈ R 20×100 with the maximal column-wise sparsity level 3 is manually created. A dictionary D ∈ R 5×20 is generated from the IPR incoherent dictionary learning algorithm [START_REF] Barchiesi | Learning incoherent dictionaries for sparse approximation using iterative projections and rotations[END_REF] on an arbitrary image, with the coherence parameter set to 0.6; The obtained dictionary has a coherence computed by (5.1) of 0.608. Then, in each test, the set of signals X can be generated with

where the second term in the right-hand-side corresponds to the unfitness noise, where E a white Gaussian zero-mean matrix with a noise level set to ω = 0.1.

To provide an overall evaluation of the proposed algorithm, several coherence parameter values µ c = {0.5, 0.55, 0.6, 0.7, 0.8, 0.9, 1} are used. It is worth noting that values below 0.4 cannot be reached due to geometric constraints, namely the coherence of an overcomplete dictionary of size n × m is bounded by (5.2). The other parameters values are set as follows: For Algorithm 12, we set: ADMM algorithm, the maximal outer iteration number N i t er = 50, the coefficient to update the penalty parameter ρ 1 = ρ 2 = 1.5, the stopping criterion = 0.01; For Algorithm 13, we set:

EPALM algorithm, the maximal inner iteration number N i i t er = 1000, the stopping criterion i = 0 = 0.01. For each coherence value, five independent Monte Carlo simulations are conducted.

We analyze the algorithm through the accuracy D -D * F and the objective function value 1 2 X -DA 2 F , as well as the computing time, where D is the output of the algorithm and D * is the known optimal solution. Moreover, we also study the iteration numbers with different coherence parameter settings. The results are listed in Table 5.1. These results show that, with the decrease of the coherent parameter µ c , more iterations are needed to converge, and thus more time.

On the other hand, as the coherence parameter increases, the stopping criteria max (|d T p d p -1|) ≤ 0.01 and max |G -D T D| ≤ 0.01 can be easily satisfied. It is observable that when µ c = 0.6, which is the closest value to the coherence of the target dictionary D * (i.e., µ * = 0.607), the results have the greatest accuracy of 0.058. For the other values of µ c , the results remain consistent but with a deduced accuracy. This is easy to understand since, for µ c > µ * , the optimal solution D * is in the feasible region, which should also be the output of the algorithm. But, influenced by the

Résumé

Abstract

In this monograph, we study the exact 0 based sparse representation problem. For the classical dictionary learning problem, the solution is obtained by iteratively processing two steps: sparse coding and dictionary updating. However, even the problem associated with sparse coding is non-convex and NP-hard. The method for solving this is to reformulate the problem as mixed integer quadratic programming (MIQP). Then by introducing two optimization techniques, initialization by proximal method and relaxation with augmented contraints, the algorithm is greatly speed up (which is thus called AcMIQP) and applied in image denoising, which shows the good performance. Moreover, the classical problem is extended to learn an incoherent dictionary. For dealing with this problem, AcMIQP or proximal method is used for sparse coding. As for dictionary updating, augmented Lagrangian method (ADMM) and extended proximal alternating linearized minimizing method are combined. This exact 0 based incoherent dictionary learning is applied in image recovery, which illustrates the improved performance with a lower coherence. XXIII