Julien Mairal
email: julien.mairal@inria.fr

Francis Bach
email: francis.bach@inria.fr

Jean Ponce
email: jean.ponce@ens.fr

Sparse Modeling for Image and Vision Processing

Keywords: 5, 3 Iterative reweighted-ℓ 1 methods, 5, 4 Iterative reweighted-ℓ 2 methods, 5, 5 Optimization for dictionary learning, 5, 6 Other optimization techniques

ℓ 0 -penalty . . . . . . . . . 5.2 Sparse reconstruction with the ℓ 1 -norm . .

A Short Introduction to Parsimony

In its most general definition, the principle of sparsity, or parsimony, consists of representing some phenomenon with as few variables as possible. It appears to be central to many research fields and is often considered to be inspired from an early doctrine formulated by the philosopher and theologian William of Ockham in the 14th century, which essentially favors simple theories over more complex ones. Of course, the link between Ockham and the tools presented in this monograph is rather thin, and more modern views seem to appear later in the beginning of the 20th century. Discussing the scientific method, [START_REF] Wrinch | XLII. On certain fundamental principles of scientific inquiry[END_REF] introduce indeed a simplicity principle related to parsimony as follows:

The existence of simple laws is, then, apparently, to be regarded as a quality of nature; and accordingly we may infer that it is justifiable to prefer a simple law to a more complex one that fits our observations slightly better.

Remarkably, [START_REF] Wrinch | XLII. On certain fundamental principles of scientific inquiry[END_REF] further discuss statistical modeling of physical observations and relate the concept of "simplicity" to the number of learning parameters; as a matter of fact, this concept is relatively close to the contemporary view of parsimony.

2

Subsequently, numerous tools have been developed by statisticians to build models of physical phenomena with good predictive power. Models are usually learned from observed data, and their generalization performance is evaluated on test data. Among a collection of plausible models, the simplest one is often preferred, and the number of underlying parameters is used as a criterion to perform model selection [START_REF] Mallows | Choosing variables in a linear regression: A graphical aid[END_REF][START_REF] Mallows | Choosing a subset regression[END_REF][START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF][START_REF] Hocking | A Biometrics invited paper. The analysis and selection of variables in linear regression[END_REF][START_REF] Barron | The minimum description length principle in coding and modeling[END_REF][START_REF] Rissanen | Modeling by shortest data description[END_REF][START_REF] Schwarz | Estimating the dimension of a model[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF].

In signal processing, similar problems as in statistics arise, but a different terminology is used. Observations, or data vectors, are called "signals", and data modeling appears to be a crucial step for performing various operations such as restoration, compression, or for solving inverse problems. Here also, the sparsity principle plays an important role and has been successful [START_REF] Mallat | Matching pursuit in a time-frequency dictionary[END_REF][START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]. Each signal is approximated by a sparse linear combination of prototypes called dictionary elements, resulting in simple and compact models.

However, statistics and signal processing remain two distinct fields with different objectives and methodology; specifically, signals often come from the same data source, e.g., natural images, whereas problems considered in statistics are unrelated to each other in general. Then, a long series of works has been devoted to finding appropriate dictionaries for signal classes of interest, leading to various sorts of wavelets [START_REF] Freeman | The design and use of steerable filters[END_REF][START_REF] Simoncelli | Shiftable multiscale transforms[END_REF][START_REF] Donoho | Wedgelets: Nearly minimax estimation of edges[END_REF][START_REF] Candès | Recovering edges in ill-posed inverse problems: Optimality of curvelet frames[END_REF][START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF][START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Even though statistics and signal processing have devised most of the methodology of sparse modeling, the parsimony principle was also discovered independently in other fields. To some extent, it appears indeed in the work of [START_REF] Markowitz | Portfolio selection[END_REF] about portfolio selection in finance, and also in geophysics [START_REF] Claerbout | Robust modeling with erratic data[END_REF]Muir, 1973, Taylor et al., 1979].

In neuroscience, Olshausen andField [1996, 1997] proposed a significantly different approach to sparse modeling than previously established practices. Whereas classical techniques in signal processing were using fixed off-the-shelf dictionaries, the method of Olshausen andField [1996, 1997] consists of learning it from training data. In a pioneer exploratory experiment, they demonstrated that dictionary learning could easily discover underlying structures in natural image patches; later, their approach found numerous applications in many fields, notably in image and audio processing [START_REF] Lewicki | Efficient coding of natural sounds[END_REF], Elad and Aharon, 2006[START_REF] Mairal | Non-local sparse models for image restoration[END_REF], Yang et al., 2010a] and computer vision [START_REF] Raina | Self-taught learning: transfer learning from unlabeled data[END_REF], Yang et al., 2009[START_REF] Zeiler | Adaptive deconvolutional networks for mid and high level feature learning[END_REF], Mairal et al., 2012[START_REF] Song | Sparselet models for efficient multiclass object detection[END_REF], Castrodad and Sapiro, 2012[START_REF] Elhamifar | See all by looking at a few: Sparse modeling for finding representative objects[END_REF][START_REF] Pokrass | Sparse modeling of intrinsic correspondences[END_REF].

The goal of this monograph is to present basic tools of sparse modeling and their applications to visual recognition and image processing. We aim at offering a self-contained view combining pluri-disciplinary methodology, practical advice, and a large review of the literature. Most of the figures in the paper are produced with the software SPAMS1 , and the corresponding Matlab code will be provided on the first author's webpage.

The monograph is organized as follows: the current introductory section is divided into several parts providing a simple historical view of sparse estimation. In Section 1.1, we start with early concepts of parsimony in statistics and information theory from the 70's and 80's. We present the use of sparse estimation within the wavelet framework in Section 1.2, which was essentially developed in the 90's. Section 1.3 introduces the era of "modern parsimony"-that is, the ℓ 1 -norm and its variants, which have been heavily used during the last two decades. Section 1.4 is devoted to the dictionary learning formulation originally introduced by Olshausen andField [1996, 1997], which is a key component of most applications presented later in this monograph. In Sections 1.5 and 1.6, we conclude our introductory tour with some theoretical aspects, such as the concept of "compressed sensing" and sparse recovery that has attracted a large attention in recent years.

With all these parsimonious tools in hand, we discuss the use of sparse coding and related sparse matrix factorization techniques for discovering the underlying structure of natural image patches in Sec-tion 2 . Even though the task here is subjective and exploratory, it is the first successful instance of dictionary learning; the insight gained from these early experiments forms the basis of concrete applications presented in subsequent sections.

Section 3 covers numerous applications of sparse models of natural image patches in image processing, such as image denoising, superresolution, inpainting, or demosaicking. This section is concluded with other related patch-modeling approaches.

Section 4 presents recent success of sparse models for visual recognition, such as codebook learning of visual descriptors, face recognition, or more low-level tasks such as edge detection and classification of textures and digits. We conclude the section with other computer vision applications such as visual tracking and data visualization.

Section 5 is devoted to optimization algorithms. It presents in a concise way efficient algorithms for solving sparse decomposition and dictionary learning problems.

We see our monograph as a good complement of other books and monographs about sparse estimation, which offer different perspectives, such as [START_REF] Mallat | A wavelet tour of signal processing[END_REF], [START_REF] Elad | Sparse and redundant representations: from theory to applications in signal and image processing[END_REF] in signal and image processing, or Bach et al. [2012a] in optimization and machine learning. We also put the emphasis on the structure of natural image patches learned with dictionary learning, and thus present an alternative view to the book of [START_REF] Hyvärinen | Natural Image Statistics: A Probabilistic Approach to Early Computational Vision[END_REF], which is focused on independent component analysis.

Notation. In this monograph, vectors are denoted by bold lower-case letters and matrices by upper-case ones. For instance, we consider in the rest of this paragraph a vector x in R n and a matrix X in R m×n . The columns of X are represented by indexed vectors x 1 , . . . , x n such that we can write X = [x 1 , . . . , x n ]. The i-th entry of x is denoted by x[i], and the i-th entry of the j-th column of X is represented by X[i, j]. For any subset g of {1, . . . , n}, we denote by x[g] the vector in R |g| that records the entries of x corresponding to indices in g. For q ≥ 1, we define the ℓ q -norm of x as x q ( n i=1 |x[i]| q ) 1/q , and the ℓ ∞ -norm as x ∞ lim q→+∞ x q = max i=1,...,n |x[i]|.

For q < 1, we define the ℓ q -penalty as x q n i=1 |x[i]| q , which, with an abuse of terminology, is often referred to as ℓ q -norm. The ℓ 0penalty simply counts the number of non-zero entries in a vector: x 0 ♯{i s.t. x[i] = 0}. For a matrix X, we define the Frobenius

norm X F = m i=1 n j=1 X[i, j] 2 1/2
. When dealing with a random variable X defined on a probability space, we denote its expectation by E[X], assuming that there is no measurability or integrability issue.

Early concepts of parsimony in statistics

A large number of statistical procedures can be formulated as maximum likelihood estimation. Given a statistical model with parameters θ, it consists of minimizing with respect to θ an objective function representing the negative log-likelihood of observed data. Assuming for instance that we observe independent samples z 1 , . . . , z n of the (unknown) data distribution, we need to solve min θ∈R p L(θ) -n i=1 log P θ (z i ) ,

(1.1)

where P is some probability distribution parameterized by θ.

Simple methods such as ordinary least squares can be written as (1.1). Consider for instance data points z i that are pairs (y i , x i ), with y i is an observation in R and x i is a vector in R p , and assume that there exists a linear relation y i = x ⊤ i θ+ε i , where ε i is an approximation error for observation i. Under a model where the ε i 's are independent and identically normally distributed with zero-mean, Eq. (1.1) is equivalent to a least square problem:

2 min θ∈R p n i=1 1 2 y i -x ⊤ i θ 2 .
To prevent overfitting and to improve the interpretability of the learned model, it was suggested in early work that a solution involving only a few model variables could be more appropriate than an exact solution of (1.1); in other words, a sparse solution involving only-let us sayk variables might be desirable in some situations. Unfortunately, such a strategy yields two difficulties: first, it is not clear a priori how to choose k; second, finding the best subset of k variables is NP-hard in general [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. The first issue was addressed with several criterions for controlling the trade-off between the sparsity of the solution θ and the adequacy of the fit to training data. For the second issue, approximate computational techniques have been proposed.

Mallows's C p , AIC, and BIC. For the ordinary least squares problem, [START_REF] Mallows | Choosing variables in a linear regression: A graphical aid[END_REF][START_REF] Mallows | Choosing a subset regression[END_REF] introduced the C p -statistics, later generalized by [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] with the Akaike information criterion (AIC), and then by [START_REF] Schwarz | Estimating the dimension of a model[END_REF] with the Bayesian information criterion (BIC). Using C p , AIC, or BIC is equivalent to solving the penalized ℓ 0 -maximum likelihood estimation problem min θ∈R p L(θ) + λ θ 0 , (1.2)

where λ depends on the chosen criterion [see [START_REF] Hastie | The elements of statistical learning[END_REF], and θ 0 is the ℓ 0 -penalty. Similar formulations have also been derived by using the minimum description length (MDL) principle for model selection [START_REF] Rissanen | Modeling by shortest data description[END_REF][START_REF] Barron | The minimum description length principle in coding and modeling[END_REF]]. As shown by [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF], the problem (1.2) is NP-hard, and approximate algorithms are necessary unless p is very small, e.g., p < 30.

Forward selection and best subset selection for least squares. To obtain an approximate solution of (1.2), a classical approach is the forward selection technique, which is a greedy algorithm that solves a sequence of maximum likelihood estimation problems computed on a subset of variables. After every iteration, a new variable is added to the subset according to the chosen sparsity criterion in a greedy manner. Some variants allow backward steps-that is, a variable can possibly exit the active subset after an iteration. The algorithm is presented in more details in Section 5.1 and seems to be due to [START_REF] Efroymson | Multiple regression analysis[END_REF], according to [START_REF] Hocking | A Biometrics invited paper. The analysis and selection of variables in linear regression[END_REF]. Other approaches considered in the 70's include also the leaps and bounds technique of [START_REF] Furnival | Regressions by leaps and bounds[END_REF], a branch-and-bound algorithm providing the exact solution of (1.2) with exponential worst-case complexity.

Wavelets in signal processing

In signal processing, similar problems as in statistics have been studied in the context of wavelets. In a nutshell, a wavelet basis represents a set of functions φ 1 , φ 2 , . . . that are essentially dilated and shifted versions of each other. Unlike Fourier basis, wavelets have the interesting properties to be localized both in the space and frequency domains, and to be suitable to multi-resolution analysis of signals [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF].

The concept of parsimony is central to wavelets. When a signal f is "smooth" in a particular sense [see [START_REF] Mallat | A wavelet tour of signal processing[END_REF], it can be well approximated by a linear combination of a few wavelets. Specifically, f is close to an expansion i α i φ i where only a few coefficients α i are non-zero, and the resulting compact representation has effective applications in estimation and compression. The wavelet theory is well developed for continuous signals, e.g., f is chosen in the Hilbert space L 2 (R), but also for discrete signals f in R m , making it suitable to modern digital image processing.

Since the first wavelet was introduced by [START_REF] Haar | Zur theorie der orthogonalen funktionensysteme[END_REF], much research has been devoted to designing a wavelet set that is adapted to particular signals such as natural images. After a long quest for finding good orthogonal basis such as the one proposed by [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF], a series of works has focused on wavelet sets whose elements are not linearly independent. It resulted a large number of variants, such as steerable wavelets [START_REF] Simoncelli | Shiftable multiscale transforms[END_REF], curvelets [START_REF] Candès | Recovering edges in ill-posed inverse problems: Optimality of curvelet frames[END_REF], contourlets [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF], or bandlets [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF]. For the purpose of our monograph, one concept related to sparse estimation is particularly important; it is called wavelet thresholding.

Sparse estimation and wavelet thresholding. Let us consider a discrete signal represented by a vector x in R p and an orthogonal wavelet basis set D = [d 1 , . . . , d p ]-that is, satisfying D ⊤ D = I where I is the identity matrix. Approximating x by a sparse linear combination of wavelet elements can be formulated as finding a sparse vector α in R p , say with k non-zero coefficients, that minimizes min α∈R p 1 2

x -Dα 2 2 s.t. α 0 ≤ k.

(1.

3)

The sparse decomposition problem (1.3) is an instance of the best subset selection formulation presented in Section 1.1 where α represents model parameters, demonstrating that similar topics arise in statistics and signal processing. However, whereas (1.3) is NP-hard for general matrices D [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF], we have assumed D to be orthogonal in the context of wavelets. As such, (1.3) is equivalent to

min α∈R p 1 2 D ⊤ x -α 2 2 s.t. α 0 ≤ k,
and admits a closed form. Let us indeed define the vector β D ⊤ x in R p , corresponding to the exact non-sparse decomposition of x onto D-that is, we have x = Dβ since D is orthogonal. To obtain the best k-sparse approximation, we denote by µ the k-th largest value among the set {|β[1]|, . . . , |β[p]|}, and the solution α ht of (1.3) is obtained by applying to β an operator called "hard-thresholding" and defined as (1.4) where 1 |β [i]|≥µ is the indicator function, which is equal to 1 if |β[i]| ≥ µ and 0 otherwise. In other words, the hard-thresholding operator simply sets to zero coefficients from β whose magnitude is below the threshold µ. The corresponding procedure, called "wavelet thresholding", is simple and effective for image denoising, even though it does not perform as well as recent state-of-the-art techniques presented in Section 3. When an image x is noisy, e.g., corrupted by white Gaussian noise, and µ is well chosen, the estimate Dα ht is a good estimate of the clean original image. The terminology "hard" is defined in contrast to an important variant called the "soft-thresholding operator", which was introduced by [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] in the context of wavelets:

α ht [i] = 1 |β[i]|≥µ β[i] = β[i] if |β[i]| ≥ µ, 0 otherwise,
3 α st [i] sign(β[i]) max(|β[i]| -λ, 0) =      β[i] -λ if β[i] ≥ λ, β[i] + λ if β[i] ≤ -λ, 0
otherwise, (1.5) where λ is a parameter playing the same role as µ in (1.4). Not only does the operator set small coefficients of β to zero, but it also reduces the magnitude of the non-zero ones. Both operators are illustrated and compared to each other in Figure 1.1. Interestingly, whereas α ht is the solution of (1.3) when µ corresponds to the entry of β = D ⊤ x with kth largest magnitude, α st is in fact the solution of the following sparse reconstruction problem with the orthogonal matrix D:

min α∈R p 1 2
x -Dα 2 2 + λ α 1 .

(1.6) This formulation will be the topic of the next section for general nonorthogonal matrices. Similar to statistics where choosing the parameter k of the best subset selection was difficult, automatically selecting the best thresholds µ or λ has been a major research topic [see, e.g. [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF], Chang et al., 2000a,b].

Structured group thresholding. Wavelets coefficients have a particular structure since the basis elements d i are dilated and shifted versions of each other. It is for instance possible to define neighborhood relationships for wavelets whose spatial supports are close to each other, or hierarchical relationships between wavelets with same or similar localization but with different scales. For one-dimensional signals, we present in Figure 1.2 a typical organization of wavelet coefficients on a tree with arrows representing such relations. For two-dimensional images, the structure is slightly more involved and the coefficients are usually organized as a collection of quadtrees [see [START_REF] Mallat | A wavelet tour of signal processing[END_REF], for more details]; we present such a configuration in Figure 1.3.

A natural idea has inspired the recent concept of group sparsity that will be presented in the next section; it consists in exploiting the wavelet structure to improve thresholding estimators. Specifically, it is different purpose. Nodes represent wavelet coefficients and their depth in the tree correspond to the scale parameter of the wavelet. We also illustrate the zero-tree coding scheme [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF] in this figure. Empty nodes correspond to zero coefficient: according to the zero-tree coding scheme, their descendants in the tree are also zero. possible to use neighborhood relations between wavelet basis elements to define groups of coefficients that form a partition G of {1, . . . , p}, and use a group-thresholding operator [START_REF] Hall | On the minimax optimality of block thresholded wavelet estimators[END_REF][START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF] defined for every group g in G as

α gt [g] 1 -λ β[g] 2 β[g] if β[g] 2 ≥ λ, 0 otherwise, (1.7)
where β[g] is the vector of size |g| recording the entries of β whose indices are in g. By using such an estimator, groups of neighbor coefficients are set to zero together when their joint ℓ 2 -norm falls below the threshold λ. Interestingly, even though the next interpretation does not appear in early work about group-thresholding [START_REF] Hall | On the minimax optimality of block thresholded wavelet estimators[END_REF][START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF], it is possible to view α gt with β = D ⊤ x as the solution of the following penalized problem min

α∈R p 1 2 x -Dα 2 2 + λ g∈G α[g] 2 , (1.8)
where the closed-form solution (1.7) holds because D is orthogonal [see Bach et al., 2012a]. Such a formulation will be studied in the next section for general matrices. Finally, other ideas for exploiting both structure and wavelet parsimony have been proposed. One is a coding scheme called "zero-tree" wavelet coding [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF], which uses the tree structure of wavelets to force all descendants of zero coefficients to be zero as well. Equivalently, a coefficient can be non-zero only if its parent in the tree is non-zero, as illustrated in Figure 1.2. This idea has been revisited later in a more general context by [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF]. Other complex models have been used as well for modeling interactions between coefficients: we can mention the application of hidden Markov models (HMM) to wavelets by [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF] and the Gaussian scale mixture model of [START_REF] Portilla | Image denoising using scale mixtures of Gaussians in the wavelet domain[END_REF].

Modern parsimony: the ℓ 1 -norm and other variants

The era of "modern" parsimony corresponds probably to the use of convex optimization techniques for solving feature selection or sparse decomposition problems. Even though the ℓ 1 -norm was introduced for that purpose in geophysics [START_REF] Claerbout | Robust modeling with erratic data[END_REF]Muir, 1973, Taylor et al., 1979], it was popularized in statistics with the Lasso estimator of [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] and independently in signal processing with the basis pursuit formulation of [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Given observations x in R n and a matrix of predictors D in R n×p , the Lasso consists of learning a linear model x ≈ Dα by solving the following quadratic program:

min α∈R p 1 2
x -Dα 2 2 s.t. α 1 ≤ µ.

(1.9)

As detailed in the sequel, the ℓ 1 -norm encourages the solution α to be sparse and the parameter µ is used to control the trade-off between data fitting and the sparsity of α. In practice, reducing the value of µ leads indeed to sparser solution in general, i.e., with more zeroes, even though there is no formal relation between the sparsity of α and its ℓ 1norm for general matrices D.

The basis pursuit denoising formulation of [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] is relatively similar but the ℓ 1 -norm is used as a penalty instead of a constraint. It can be written as min α∈R p 1 2

x -Dα 2 2 + λ α 1 , (1.10) which is essentially equivalent to (1.9) from a convex optimization perspective, and in fact (1.10) is also often called "Lasso" in the literature. Given some data x, matrix D, and parameter µ > 0, we indeed know from Lagrange multiplier theory [see, e.g., [START_REF] Borwein | Convex analysis and nonlinear optimization: theory and examples[END_REF]Lewis, 2006, Boyd and[START_REF] Boyd | Convex Optimization[END_REF]] that for all solution α ⋆ of (1.9), there exists a parameter λ ≥ 0 such that α ⋆ is also a solution of (1.10). We note, however, that there is no direct mapping between λ and µ, and thus the choice of formulation (1.9) or (1.10) should be made according to how easy it is to select the parameters λ or µ. For instance, one may prefer (1.9) when a priori information about the ℓ 1 -norm of the solution is available. In Figure 1.4, we illustrate the effect of changing the value of the regularization parameter λ on the solution of (1.10) for two datasets. When λ = 0, the solution is dense; in general, increasing λ sets more and more variables to zero. However, the relation between λ and the sparsity of the solution is not exactly monotonic. In a few cases, increasing λ yields a denser solution.

Another "equivalent" formulation consists of finding a sparse decomposition under a reconstruction constraint:

min α∈R p α 1 s.t. x -Dα 2 2 ≤ ε. (1.11)
This formulation can be useful when we have a priori knowledge about the noise level and the parameter ε is easy to choose. The link between (1.10) and (1.11) is similar to the link between (1.10) and (1.9). For noiseless problems, [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] have also introduced a formulation simply called "basis pursuit" (without the terminology "denoising"), defined as

min α∈R p α 1 s.t. x = Dα, (1.12)
which is related to (1.10) in the sense that the set of solutions of (1.10) converges to the solutions of (1.12) when λ converges to 0 + , whenever the linear system x = Dα is feasible. These four formulations (1.9-1.12) have gained a large success beyond the statistics and signal processing communities. More generally, the ℓ 1 -norm has been used as a regularization function beyond the least-square context, leading to problems of the form min

α∈R p f (α) + λ α 1 , (1.13)
where f : R p → R is a loss function. In the rest of this section, we will present several variants of the ℓ 1 -norm, but before that, we will try to understand why such a penalty is appropriate for sparse estimation.

Why does the ℓ 1 -norm induce sparsity? Even though we have claimed that there is no rigorous relation between the sparsity of α and its ℓ 1 -norm in general, intuition about the sparsity-inducing effect of the ℓ 1 -norm may be obtained from several viewpoints.

Analytical point of view. In the previous section about wavelets, we have seen that when D is orthogonal, the ℓ 1 -decomposition problem (1.10) admits an analytic closed form solution (1.5) obtained by soft-thresholding. As a result, whenever the magnitude of the inner product d ⊤ i x is smaller than λ for an index i, the corresponding variable α ⋆ [i] is equal to zero. Thus, the number of zeroes of the solution α ⋆ monotonically increases with λ.

For non-orthogonal matrices D, such a monotonic relation does not formally hold anymore; in practice, the sparsity-inducing property of the ℓ 1 -penalty remains effective, as illustrated in Figure 1.4. Some intuition about this fact can be gained by studying optimality conditions for the general ℓ 1 -regularized problem (1.13) where f is a differentiable function. The following lemma details these conditions. Lemma 1.1 (Optimality conditions for ℓ 1 -regularized problems).

A vector α ⋆ in R p is a solution of (1.13) if and only if ∀i = 1, . . . , p -∇f (α ⋆ )[i] = λ sign(α ⋆ [i]) if α ⋆ [i] = 0, |∇f (α ⋆ )[i]| ≤ λ otherwise.
(1.14)

Proof. A proof using the classical concept of subdifferential from convex optimization can be found in [see,e.g., Bach et al., 2012a]. Here, we provide instead an elementary proof using the simpler concept of directional derivative for nonsmooth functions, defined as, when the limit exists,

∇g(α, κ) lim t→0 + g(α + tκ) -g(α) t ,
for a function g : R p → R at a point α in R p and a direction κ in R p . For convex functions g, directional derivatives always exist and a classical optimality condition for α ⋆ to be a minimum of g is to have ∇g(α ⋆ , κ) non-negative for all directions κ [START_REF] Borwein | Convex analysis and nonlinear optimization: theory and examples[END_REF]. Intuitively, this means that one cannot find any direction κ such that an infinitesimal move along κ from α ⋆ decreases the value of the objective. When g is differentiable, the condition is equivalent to the classical optimality condition ∇g(α ⋆ ) = 0.

We can now apply the directional derivative condition to the func-

tion g : α → f (α) + λ α 1 , which is equivalent to ∀κ ∈ R p , ∇f (α ⋆ ) ⊤ κ + λ p i=1 sign(α ⋆ [i])κ[i] if α ⋆ [i] = 0, |κ[i]| otherwise ≥ 0.
(1.15) It is then easy to show that (1.15) holds for all κ if and only the inequality holds for the specific values κ = e i and κ = -e i for all i, where e i is the vector in R p with zeroes everywhere except for the i-th entry that is equal to one. This immediately provides an equivalence between (1.15) and (1.14).

Lemma 1.1 is interesting from a computational point of view (see Section 5.2), but it also tells us that when λ ≥ ∇f (0) ∞ , the conditions (1.14) are satisfied for α ⋆ = 0, the sparsest solution possible.

Physical point of view.

In image processing or computer vision, the word "energy" often denotes the objective function of a minimization problem; it is indeed common in physics to have complex systems that stabilize at a configuration of minimum potential energy. The negative of the energy's gradient represents a force, a terminology we will borrow in this paragraph. Consider for instance a one-dimensional ℓ 1 -regularized estimation problem

min α∈R 1 2 (β -α) 2 + λ|α|, (1.16)
where β is a positive constant. Whenever α is non-zero, the ℓ 1penalty is differentiable with derivative λ sign(α). When interpreting this objective as an energy minimization problem, the ℓ 1 -penalty can be seen as applying a force driving α towards the origin with constant intensity λ. Consider now instead the squared ℓ 2 -penalty, also called regularization of [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF], or ridge regression regularization [START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF]:

min α∈R 1 2 (β -α) 2 + λ 2 α 2 .
(1.17)

The derivative of the quadratic energy (λ/2)α 2 is λα. It can be interpreted as a force that also points to the origin but with linear intensity λ|α|. Therefore, the force corresponding to the ridge regularization can be arbitrarily strong when α is large, but if fades away when α gets close to zero. As a result, the squared ℓ 2 -regularization does not have a sparsity-inducing effect. From an analytical point of view, we have seen that the solution of (1.16) is zero when |β| is smaller than λ. In contrast, the solution of (1.17) admits a closed form α ⋆ = β/(1 + λ).

And thus, regardless of the parameter λ, the solution is never zero. We present a physical example illustrating this phenomenon in Figure 1.5. We use springs whose potential energy is known to be quadratic, and objects with a gravitational potential energy that is approximately linear on the Earth's surface.

Geometrical point of view. The sparsity-inducing effect of the ℓ 1norm can also be interpreted by studying the geometry of the ℓ 1 -ball {α ∈ R p : α 1 ≤ µ}. More precisely, understanding the effect of the Euclidean projection onto this set is important: in simple cases where the design matrix D is orthogonal, the solution of (1.9) can indeed be obtained by the projection min is the same as on the left, whereas E2 is a gravitational potential energy, where g is the gravitational constant on the Earth's surface. Both system can evolve according to their initial positions, and stabilize for the value of α ⋆ that minimizes the energy E1 + E2, assuming that some energy can be dissipated by friction forces.

α∈R p 1 2 β -α 2 2 s.t. α 1 ≤ µ, (1.18) E 1 = k 1 2 (β -α) 2 E 2 = k 2 2 α 2 α β α E 1 = k 1 2 (β -α) 2 E 2 = mg|α|, α ≥ 0 (a) Small regularization E 1 = k 1 2 (β -α) 2 E 2 = k 2 2 α 2 α α ⋆ = 0 E 1 = k 1 2 (β -α) 2 E 2 = mg|α|, α ≥ 0 (b) High regularization
On the left, it is possible to show that α ⋆ = βk1/(k1 + k2) and thus, the solution α ⋆ is never equal to zero, regardless of the strength k2 of the bottom spring. On the right, the solution is obtained by soft-thresholding: α ⋆ = max(β -mg/k1, 0). As shown on Figure 1.5(b), when the mass m is large enough, the blue object touches the ground and α ⋆ = 0. Figure adapted from [START_REF]Sparse coding for machine learning, image processing and computer vision[END_REF].

where β = D ⊤ x. When D is not orthogonal, a classical algorithm for solving (1.9) is the projected gradient method (see Section 5.2), which performs a sequence of projections (1.18) for different values of β. Note that how to solve (1.18) efficiently is well studied; it can be achieved in O(p) operations with a divide-and-conquer strategy [START_REF] Brucker | An O(n) algorithm for quadratic knapsack problems[END_REF][START_REF] Duchi | Efficient projections onto the l1-ball for learning in high dimensions[END_REF].

In Figure 1.6, we illustrate the effect of the ℓ 1 -norm projection and compare it to the case of the ℓ 2 -norm. The corners of the ℓ 1 -ball are on the main axes and correspond to sparse solutions. Two of them are represented by red and green dots, with respective coordinates (µ, 0) and (0, µ). Most strikingly, a large part of the space in the figure, represented by red and green regions, ends up on these corners after projection. In contrast, the set of points that is projected onto the blue dot, is simply the blue line. The blue dot corresponds in fact to a dense solution with coordinates (µ/2, µ/2). Therefore, the figure illustrates that the ℓ 1 -ball in two dimensions encourages solutions to be on its corners. In the case of the ℓ 2 -norm, the ball is isotropic, and treats every direction equally. In Figure 1.7, we represent these two balls in three dimensions, where we can make similar observations.

More formally, we can mathematically characterize our remarks about Figure 1.6. Consider a point y in R p on the surface of the ℓ 1 -ball of radius µ = 1, and define the set N {z ∈ R p : π(z) = y}, where π is the projection operator onto the ℓ 1 -ball. Examples of pairs (y, N ) have been presented in Figure 1.6; for instance, when y is the red or green dot, N is respectively the red or green region. It is particularly informative to study how N varies with y, which is the focus of the next proposition.

Proposition 1.1 (Characterization of the set N ). For a non-zero vector y in R p , the set N defined in the previous paragraph can be written as N = y + K, where K is a polyhedral cone of dimension py 0 + 1. Bertsekas, 1999, Proposition B.11] allows us to rewrite N as

Proof. A classical theorem [see

N = {z ∈ R p : ∀ x 1 ≤ 1, (z -y) ⊤ (x -y) ≤ 0} = y + K,
where y + K denotes the Minkowski sum {y + z : z ∈ K} between the set {y} and the cone K defined as

K {d ∈ R p : ∀ x 1 ≤ 1, d ⊤ (x -y) ≤ 0}.
Note that in the optimization literature, K is often called the "normal cone" to the unit ℓ 1 -ball at the point y [START_REF] Borwein | Convex analysis and nonlinear optimization: theory and examples[END_REF]. Equivalently, we have

K = {d ∈ R p : max x 1 ≤1 d ⊤ x ≤ d ⊤ y} = {d ∈ R p : d ∞ ≤ d ⊤ y}, (1.19)
where we have used the fact that quantity max Bach et al., 2012a]. Note now that according to Hölder's inequality, we also have .19). Therefore, the inequalities are in fact equalities. It is then easy to characterize vectors d such that d ⊤ y = d ∞ y 1 and it is possible to show that K is simply the set of vectors d satisfying d

x 1 ≤1 d ⊤ x, called the dual-norm of the ℓ 1 -norm, is equal to d ∞ [see
d ⊤ y ≤ d ∞ y 1 ≤ d ∞ in Eq. (1
[i] = sign(y[i]) d ∞ for all i such that y[i] = 0.
This would be sufficient to conclude the proposition, but it is also possible to pursue the analysis and exactly characterize K by finding a set of generators.4 Let us define the vector s in {-1, 0, +1} p that carries the sparsity pattern of y, more precisely, with s[i] = sign(y[i]) for all i such that y[i] = 0, and s[i] = 0 otherwise. Let us also define the set of indices {i 1 , . . . , i l } corresponding to the l zero entries of y, and e i in R p the binary vector whose entries are all zero but the i-th one that is equal to 1. Then, after a short calculation, we can geometrically characterize the polyhedral cone K:

K = cone (s, s -e i 1 , s + e i 1 , s -e i 2 , s + e i 2 , . . . , s -e i l , s + e i l ) ,
where the notation "cone" is defined in footnote 4.

It is now easy to see that the set K "grows" with the number l of zero entries in y, and that K lives in a subspace of dimension l + 1 for all non-zero vector y. For example, when l = 0-that is, y is a dense vector (e.g., the blue point in Figure 1.6(a)), K is simply a half-line.

To conclude, the geometrical intuition to gain from this section is that the Euclidean projection onto a convex set encourages solutions on singular points, such as edges or corners for polytopes. Such a principle indeed applies beyond the ℓ 1 -norm. For instance, we illustrate the regularization effect of the ℓ ∞ -norm in Figure 1.8, whose corners coordinates have same magnitude.

Non-convex regularization. Even though it is well established that the ℓ 1 -norm encourages sparse solutions, it remains only a convex proxy of the ℓ 0 -penalty. Both in statistics and signal processing, other sparsity-inducing regularization functions have been proposed, in particular continuous relaxations of ℓ 0 that are non-convex [START_REF] Frank | A statistical view of some chemometrics regression tools[END_REF][START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF][START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF][START_REF] Gasso | Recovering sparse signals with non-convex penalties and DC programming[END_REF]. These functions are using a non-decreasing concave function ϕ : R + → R, and the sparsity-inducing penalty is defined as

ψ(α) p i=1 ϕ (|α[i]|) .
For example, the ℓ q -penalty uses ϕ : x → x q [START_REF] Frank | A statistical view of some chemometrics regression tools[END_REF], or an approximation ϕ : x → (x + ε) q ; the reweighted-ℓ 1 algorithm of [START_REF] Fazel | Matrix rank minimization with applications[END_REF], [START_REF] Fazel | Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices[END_REF], Candès et al. [2008] implicitly uses ϕ : x → log(x + ε). These penalties typically lead to intractable estimation problems, but approximate solutions can be obtained with continuous optimization techniques (see Section 5.3).

The sparsity-inducing effect of the penalties ψ is known to be stronger than ℓ 1 . As shown in Figure 1.9(a), the magnitude of the derivative of ϕ grows when one approaches zero because of its concavity. Thus, in the one-dimensional case, ψ can be interpreted as a force driving α towards the origin with increasing intensity when α gets closer to zero. In terms of geometry, we also display the ℓ q -ball in Figure 1.9(b), with the same red, blue, and green dots as in Figure 1.6. The part of the space that is projected onto the corners of the ℓ q -ball is larger than that for ℓ 1 . Interestingly, the geometrical structure of the red and green regions are also more complex. Their combinatorial The balls are represented in gray. All points from the red regions are projected onto the point of coordinates (0, µ) denoted by a red dot. Similarly, the green and blue regions are projected onto the green and blue dots, respectively. For the ℓ1-norm, a large part of the figure is filled by the red and green regions, whose points are projected to a sparse solution corresponding to a corner of the ball. For the ℓ2-norm, this is not the case: any non-sparse point-say, for instance on the blue line-is projected onto a non-sparse solution.

ℓ 1 -ball α 1 ≤ µ α[2] α[1] (a) Effect of the Euclidean projection onto the ℓ1-ball. α[2] α[1] ℓ 2 -ball α 2 ≤ µ (b)

α[2]

α nature makes the projection problem onto the ℓ q -ball more involved when q < 1.

[1] ℓ ∞ -ball α ∞ ≤ µ
The elastic-net. To cope with some instability issues of the estimators obtained with the ℓ 1 -regularization, [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] have proposed to combine the ℓ 1 -and ℓ 2 -norms with a penalty called elasticnet:

ψ(α) α 1 + γ α 2 2 .
The effect of this penalty is illustrated in Figure 1.10. Compared to Figure 1.6, we observe that the red and green regions are smaller for the elastic-net penalty than for ℓ 1 . The sparsity-inducing effect is thus less aggressive than the one obtained with ℓ 1 .

Total variation. The anisotropic total variation penalty [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] for one dimensional signals is simply the ℓ 1 -norm of finite differences

ψ(α) p-1 i=1 |α[i + 1] -α[i]|,
which encourages piecewise constant signals. It is also known in statistics under the name of "fused Lasso" [START_REF] Tibshirani | Sparsity and smoothness via the fused Lasso[END_REF]. The penalty can easily be extended to two-dimensional signals, and has been widely used for regularizing inverse problems in image processing [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF].

Group sparsity.

In some cases, variables are organized into predefined groups forming a partition G of {1, . . . , p}, and one is looking for a solution α ⋆ such that variables belonging to the same group of G are set to zero together. For example, such groups have appeared in Section 1.2 about wavelets, where G could be defined according to neighborhood relationships of wavelet coefficients. Then, when it is known beforehand that a problem solution only requires a few groups of variables to explain the data, a regularization function automatically selecting the relevant groups has been shown to improve the prediction performance or the interpretability of the solution [START_REF] Turlach | Simultaneous variable selection[END_REF],

α ϕ(|α|) = log(|α| + ε) (a) Illustration of a non-convex sparsity-inducing penalty. α[2] α[1] ℓ q -ball α q ≤ µ with q < 1 (b) ℓq-ball with q < 1.
Figure 1.9: Illustration of the sparsity-inducing effect of a non-convex penalty. In (a), we plot the non-convex penalty α → log(|α| + ε), and in (b), we present a similar figure as 1.6 for the ℓq-penalty, when choosing q < 1.

elastic-net ball

(1

-γ) α 1 + γ α 2 2 ≤ µ α[2] α[1]
Figure 1.10: Similar figure as 1.6 for the elastic-net penalty. [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], Obozinski et al., 2009, Huang and Zhang, 2010].

The group sparsity principle is illustrated in Figure 1.11(b).

An appropriate regularization function to obtain a group-sparsity effect is known as "Group-Lasso" penalty and is defined as

ψ(α) = g∈G α[g] q ,
(1.20)

where . q is either the ℓ 2 or ℓ ∞ -norm. To the best of our knowledge, such a penalty appears in the early work of [START_REF] Grandvalet | Outcomes of the equivalence of adaptive ridge with least absolute shrinkage[END_REF] and [START_REF] Bakin | Adaptive regression and model selection in data mining problems[END_REF] for q = 2, and [START_REF] Turlach | Simultaneous variable selection[END_REF] for q = ∞. It has been popularized later by [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF].

The function ψ in (1.20) is a norm, thus convex, and can be interpreted as the ℓ 1 -norm of the vector [ α[g] q ] g∈G of size |G|. Consequently, the sparsity-inducing effect of the ℓ 1 -norm is applied at the group level. The penalty is highly related to the group-thresholding approach for wavelets, since the group-thresholding estimator (1.7) is linked to ψ through Eq. (1.8).

In Figure 1.13(a), we visualize the unit ball of a Group-Lasso norm obtained when G contains two groups G = {{1, 2}, {3}}. The ball has two singularities: the top and bottom corners, corresponding to solutions where variables 1 and 2 are simultaneously set to zero, and the Non-zero coefficients are represented by gray squares. On the left (a), the vectors are obtained with a simple sparsity-inducing penalty, such as the ℓ1-norm, and the non-zero variables are scattered. In the middle figure (b), a group sparsity-inducing penalty with three groups of variables is used. The sparsity patterns respect the group structure. On the right (c), we present the results obtained with a hierarchical penalty consisting of the Group Lasso plus the ℓ1-norm. The group structure is globally respected, but some variables within a group can be discarded.

middle circle, corresponding to solutions where variable 3 only is set to zero. As expected, the geometry of the ball induces the group-sparsity effect.

Structured sparsity. Group-sparsity is a first step towards the more general idea that a regularization function can encourage sparse solutions with a particular structure. This notion is called structured sparsity and has been introduced under a large number of different point of views [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF][START_REF] Jacob | Group Lasso with overlaps and graph Lasso[END_REF], Jenatton et al., 2011a[START_REF] Baraniuk | Model-based compressive sensing[END_REF][START_REF] Huang | Learning with structured sparsity[END_REF]. To some extent, it follows the concept of group-thresholding introduced in the wavelet literature, which we have presented in Section 1.2. In this paragraph, we briefly review some of these works, but for a more detailed review, we refer the reader to [Bach et al., 2012b]. Some penalties are non-convex. For instance, [START_REF] Huang | Learning with structured sparsity[END_REF] and [START_REF] Baraniuk | Model-based compressive sensing[END_REF] propose two different combinatorial approaches based on a predefined set G of possibly overlapping groups of variables. These penalties encourage solutions whose support is in the union of a few number groups, but they lead to NP-hard optimization problems. Other penalties are convex. In particular, [START_REF] Jacob | Group Lasso with overlaps and graph Lasso[END_REF] introduce a sparsity-inducing norm that is exactly a convex relaxation of the penalty of [START_REF] Huang | Learning with structured sparsity[END_REF], even though these two approaches were independently developed at the same time. As a result, the convex penalty of [START_REF] Jacob | Group Lasso with overlaps and graph Lasso[END_REF] encourages a similar structure as the one of [START_REF] Huang | Learning with structured sparsity[END_REF].

By following a different direction, the Group-Lasso penalty (1.20) has been considered when the groups are allowed to overlap [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF], Jenatton et al., 2011a]. As a consequence, variables belonging to the same groups are encouraged to be set to zero together. It was proposed for hierarchical structures by [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF] with the following rule: whenever two groups g and h are in G, they should be either disjoint, or one should be included in another. Examples of such hierarchical group structures are given in Figures 1.11 and 1.12. The effect of the penalty is to encourage sparsity patterns that are rooted subtrees. Equivalently, a variable can be non-zero only if its parent in the tree is non-zero, which is the main property of the zero-tree coding scheme introduced in the wavelet literature [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF], and already illustrated in Figure 1.2.

Finally, Jenatton et al. [2011a] has extended the hierarchical penalty of [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF] to more general group structures, for example when variable are organized on a two-dimensional grid, encouraging neighbor variables to be simultaneously set to zero. We conclude this brief presentation of structured sparsity with Figure 1.13, where we present the unit balls of some sparsity-inducing norms. Each of them exhibits singularities and encourages particular sparsity patterns.

α 1 α 2 α 3 α 4 α 5 α 6 α 7 α 8 α 9 α 10 α 11 α 12 α 13 α 14 α 15
Figure 1.12: Illustration of the hierarchical sparsity of [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF], which generalizes the zero-tree coding scheme of [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF]. The groups of variables correspond to the red rectangles. The empty nodes represent variable that are set to zero. They are contained in three groups: {4, 8, 9}, {13}, {7, 14, 15}.

Spectral sparsity. Another form of parsimony has been devised in the spectral domain [START_REF] Fazel | A rank minimization heuristic with application to minimum order system approximation[END_REF][START_REF] Srebro | Maximum-margin matrix factorization[END_REF]. For estimation problems where model parameters are matrices, the rank has been used as a natural regularization function. The rank of a matrix is equal to the number of non-zero singular values, and thus, it can be interpreted as the ℓ 0 -penalty of the matrix spectrum. Unfortunately, due to the combinatorial nature of ℓ 0 , the rank penalization typically leads to intractable optimization problems. A natural convex relaxation has been introduced in the control theory literature by [START_REF] Fazel | A rank minimization heuristic with application to minimum order system approximation[END_REF] and consists of computing the ℓ 1norm of the spectrum-that is, simply the sum of the singular values.

The resulting penalty appears under different names, the most common ones being the trace, nuclear, or Schatten norm. It is defined for a matrix A in R p×k with k ≥ p as

A * p i=1 s i (A),
where s i (A) is the i-th singular value of A. Traditional applications of the trace norm in machine learning are matrix completion or collaborative filtering [START_REF] Pontil | Multi-task feature learning[END_REF][START_REF] Abernethy | A new approach to collaborative filtering: operator estimation with spectral regularization[END_REF]. These problems have become popular with the need of scalable recommender systems for video streaming providers. The goal is to infer movie preferences for each customer, based on their partial movie ratings. Typically, the matrix is of size p × k, where p is the number of movies and k is the number of users. Each user gives a score for a few movies, corresponding to some entries of the matrix, and the recommender system tries to infer the missing values. Similar techniques have also recently been used in other fields, such as in genomics to infer missing genetic information [START_REF] Chi | Genotype imputation via matrix completion[END_REF].

Dictionary learning

We have previously presented various formulations where a signal x in R m is approximated by a sparse linear combination of a few columns of a matrix D in R m×p . In the context of signal and image processing, this matrix is often called dictionary and its columns atoms. As seen in Section 1.2, a large amount of work has been devoted in the wavelet literature for designing a good dictionary adapted to natural images.

In neuroscience, Olshausen andField [1996, 1997] have proposed a significantly different approach to sparse modeling consisting of adapting the dictionary to training data. Because the size of natural images is too large for learning a full matrix D, they have chosen to learn the dictionary on natural image patches, e.g., of size m = 16 × 16 pixels, and have demonstrated that their method could automatically discover interpretable structures. We discuss this topic in more details in Section 2.

The motivation of Olshausen andField [1996, 1997] was to show that the structure of natural images is related to classical theories of the mammalian visual cortex. Later, dictionary learning found numerous applications in image restoration, and was shown to significantly outperform off-the-shelf bases for signal reconstruction [see, e.g., [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], Mairal et al., 2008c, 2009[START_REF] Protter | Image sequence denoising via sparse and redundant representations[END_REF], Yang et al., 2010a].

Concretely, given a dataset of n training signals X = [x 1 , . . . , x n ], dictionary learning can be formulated as the following minimization problem min

D∈C,A∈R p×n n i=1 1 2 x i -Dα i 2 2 + λψ(α i ), (1.21)
where A = [α 1 , . . . , α n ] carries the decomposition coefficients of the signals x 1 , . . . , x n , ψ is sparsity-inducing regularization function, and C is typically chosen as the following set:

C {D ∈ R m×p : ∀j d j 2 ≤ 1}.
To be more precise, [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF] proposed several choices for ψ; their experiments were for instance conducted with the ℓ 1norm, or with the smooth function ψ(α)

p j=1 log(ε + α[j] 2
), which has an approximate sparsity-inducing effect. The constraint D ∈ C was also not explicitly modeled in the original dictionary learning formulation; instead, the algorithm of [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF] includes a mechanism to control and rescale the ℓ 2 -norm of the dictionary elements. Indeed, without such a mechanism, the norm of D would arbitrarily go to infinity, leading to small values for the coefficients α i and making the penalty ψ ineffective.

The number of samples n is typically large, whereas the signal dimension m is small. The number of dictionary elements p is often chosen larger than m-in that case, the dictionary is said to be overcomplete-even though a choice p < m often leads to reasonable results in many applications. For instance, a typical setting would be to have m = 10 × 10 pixels for natural image patches, a dictionary of size p = 256, and more than 100 000 training patches.

A large part of this monograph is related to dictionary learning and thus we only briefly discuss this matter in this introduction. Section 2 is indeed devoted to unsupervised learning techniques for natural image patches, including dictionary learning; Sections 3 and 4 present a large number of applications in image processing and computer vision; how to solve (1.21) is explained in Section 5.5 about optimization.

Matrix factorization point of view. An equivalent representation of (1.21) is the following regularized matrix factorization problem min

D∈C,A∈R p×n 1 2 X -DA 2 F + λΨ(A), (1.22)
where Ψ(A) = n i=1 ψ(α i ). Even though reformulating (1.21) as (1.22) is simply a matter of using different notation, seeing dictionary learning as a matrix factorization problem opens up interesting perspectives. In particular, it makes obvious some links with other unsupervised learning approaches such as non-negative matrix factorization [START_REF] Paatero | Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values[END_REF], clustering techniques such as Kmeans, and others [see Mairal et al., 2010a]. These links will be further developed in Section 2.

Risk minimization point of view. Dictionary learning can also be seen from a machine learning point of view. Indeed, dictionary learning can be written as

min D∈C f n (D) 1 n n i=1 L(x i , D) ,
where L : R m × R m×p is a loss function defined as

L(x, D) min α∈R p 1 2 x -Dα 2 2 + λψ(α).
The quantity L(x, D) should be small if D is "good" at representing the signal x in a sparse fashion, and large otherwise. Then, f n (D) is called the empirical cost.

However, as pointed out by [START_REF] Bottou | The tradeoffs of large scale learning[END_REF], one is usually not interested in the exact minimization of the empirical cost f n (D) for a fixed n, which may lead to overfitting on the training data, but instead in the minimization of the expected cost, which measures the quality of the dictionary on new unseen data:

f (D) E x [L(x, D)] = lim n→∞ f n (D) a.s.,
where the expectation is taken relative to the (unknown) probability distribution of the data. 5The expected risk minimization formulation is interesting since it paves the way to stochastic optimization techniques when a large amount of data is available [Mairal et al., 2010a] and to theoretical analysis [START_REF] Maurer | k-dimensional coding schemes in Hilbert spaces[END_REF][START_REF] Vainsencher | The sample complexity of dictionary learning[END_REF][START_REF] Gribonval | Sample complexity of dictionary learning and other matrix factorizations[END_REF], which are developed in Sections 5.5 and 1.6, respectively.

Constrained variants.

Following the original formulation of Olshausen andField [1996, 1997], we have chosen to present dictionary learning where the regularization function is used as a penalty, even though it can also be used as a constraint as in (1.9). Then, natural variants of (1.21) are min

D∈C,A∈R p×n n i=1 1 2 x i -Dα i 2 2 s.t. ψ(α i ) ≤ µ. (1.23) or min D∈C,A∈R p×n n i=1 ψ(α i ) s.t. x i -Dα i 2 2 ≤ ε. (1.24)
Note that (1.23) and (1.24) are not equivalent to (1.21). For instance, problem (1.23) can be reformulated using a Lagrangian function [START_REF] Boyd | Convex Optimization[END_REF] as

min D∈C n i=1 max λ i ≥0 min α i ∈R p 1 2 x i -Dα i 2 2 + λ i (ψ(α i ) -µ) ,
where the optimal λ i 's are not necessarily equal to each other, and their relation with the constraint parameter µ is unknown in advance.

A similar discussion can be conducted for (1.24) and it is thus important in practice to choose one of the formulations (1.21), (1.23), or (1.24); the best one depends on the problem at hand and there is no general rule for preferring one instead of another.

Compressed sensing and sparse recovery

Finally, we conclude our historical tour of parsimony with recent theoretical results obtained in signal processing and statistics. We focus on methods based on the ℓ 1 -norm, i.e., the basis pursuit formulation of (1.9)-more results on structured sparsity-inducing norms are presented by Bach et al. [2012b].

Most analyses rely on particular assumptions regarding the problem. We start this section with a cautionary note from [START_REF] Hocking | A Biometrics invited paper. The analysis and selection of variables in linear regression[END_REF]:

The problem of selecting a subset of independent or predictor variables is usually described in an idealized setting. That is, it is assumed that (a) the analyst has data on a large number of potential variables which include all relevant variables and appropriate functions of them plus, possibly, some other extraneous variables and variable functions and (b) the analyst has available "good" data on which to base the eventual conclusions. In practice, the lack of satisfaction of these assumptions may make a detailed subset selection analysis a meaningless exercise.

In this section, we present such theoretical results where the assumptions are often not met in practice, but also results that either (1) can have an impact on the practice of sparse recovery or (2) do not need strong assumptions.

From support recovery to signal denoising. Given a signal x in R m and a dictionary D in R m×p with ℓ 2 -normalized columns, throughout this section, we assume that x is generated as x = Dα ⋆ + ε with a sparse vector α ⋆ in R p and an additive noise ε in R m . For simplicity, we consider α ⋆ and D as being deterministic while the noise is random, independent and identically distributed, with zero mean and finite variance σ 2 .

The different formulations presented earlier in Section 1.3, for instance basis pursuit, provide estimators α of the "true" vector α ⋆ . Then, the three following goals have been studied in sparse recovery, typically in decreasing order of hardness:

• support recovery and sign consistency: we want the support of α (i.e., the set of non-zero elements) to be the same or to be close to the one of α ⋆ . The problem is often called "model selection" in statistics and "support recovery" in signal processing; it is often refined to the estimation of the full sign pattern-that is, among the non-zero elements, we also want the correct sign to be estimated.

• code estimation: the distance αα ⋆ 2 should be small. In statistical terms, this correspond to the "estimation" of α ⋆ .

• signal denoising: regardless of code estimation, we simply want the distance D α -Dα ⋆ 2 to be small; the goal is not to obtain exactly α ⋆ , but simply to obtain a good denoised version of the signal x = Dα ⋆ + ε.

A good code estimation performance does imply a good denoising performance but the converse is not true in general. In most analyses, support recovery is harder than code estimation. As detailed below, the sufficient conditions for good support recovery lead indeed to good estimation.

High-dimensional phenomenon. Without sparsity assumptions, even the simplest denoising task can only achieve denoising errors of the order 1 n D α -Dα ⋆ 2 2 ≈ σ 2 p m , which is attained for ordinary leastsquares, and is the best possible [START_REF] Tsybakov | Optimal rates of aggregation[END_REF]. Thus, in order to have at least a good denoising performance (prediction performance in statistics), either the noise σ is small, or the signal dimension m (the number of samples) is much larger than the number of atoms p (the number of variables to select from).

When making the assumption that the true code α ⋆ is sparse with at most k non zeros, smaller denoising errors can be obtained. In that case, it is possible indeed to replace the scaling σ 2 p m by σ 2 k log p m . Thus, even when p is much larger than m, as long as log p is much smaller than m, we may have good prediction performance. However, this high-dimensional phenomenon currently6 comes at a price: (1) either an exhaustive search over the subsets of size k needs to be performed [START_REF] Massart | Concentration Inequalities and Model Selection: Ecole d'été de Probabilités de Saint-Flour 23[END_REF][START_REF] Bunea | Aggregation for Gaussian regression[END_REF][START_REF] Raskutti | Minimax rates of estimation for high-dimensional linear regression over-balls[END_REF] or (2) some assumptions have to be made regarding the dictionary D, which we now describe.

Sufficient conditions for high-dimensional fast rates. Most sufficient conditions have the same flavor. A dictionary behaves well if the offdiagonal elements of D ⊤ D are small, in other words, if there is little correlation between atoms. However, the notion of coherence (the maximal possible correlation between two atoms) was the first to emerge [see, e.g., [START_REF] Elad | Sparse and redundant representations: from theory to applications in signal and image processing[END_REF]Bruckstein, 2002, Gribonval and[START_REF] Gribonval | Sparse representations in unions of bases[END_REF]], but it is not sufficient to obtain a high-dimensional phenomenon.

In the noiseless setting, [START_REF] Candes | Decoding by linear programming[END_REF] and [START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF] introduced the restricted isometry property (RIP), which states that all submatrices of size k × k of D ⊤ D should be close to isometries, that is, should have all of their eigenvalues sufficiently close to one. With such an assumption, the Lasso behaves well: it recovers the true support and estimates the code α ⋆ and the signal Dα ⋆ with an error of order σ 2 k log p m . The main advantage of the RIP assumption is that one may exhibit dictionaries for which it is satisfied, usually obtained by normalizing a matrix D obtained from independent Gaussian entries, which may satisfy the condition that (k log p)/m remains small. Thus, the sufficient conditions are not vacuous. However, the RIP assumption has two main drawbacks: first, it cannot be checked on a given dictionary D without checking all O(p k ) submatrices of size k; second, it may be weakened if the goal is support recovery or simply estimation performance (code recovery).

There is therefore a need for sufficient conditions that can be checked in polynomial time while ensuring sparse recovery. However, none currently exists with the same scalings between k, p and m [see, e.g., [START_REF] Juditsky | On verifiable sufficient conditions for sparse signal recovery via ℓ 1 -minimization[END_REF]Nemirovski, 2011, d'Aspremont and[START_REF] Ghaoui | Testing the nullspace property using semidefinite programming[END_REF]. When refining to support recovery, [START_REF] Fuchs | Recovery of exact sparse representations in the presence of bounded noise[END_REF], [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF], [START_REF] Wainwright | Sharp thresholds for noisy and high-dimensional recovery of sparsity using ℓ 1 -constrained quadratic programming[END_REF] provide sufficient and necessary conditions of a similar flavor than requiring that all submatrices of size k are sufficiently close to orthogonal. For the tightest conditions, see, e.g., [START_REF] Bühlmann | Statistics for high-dimensional data: methods, theory and applications[END_REF]. Note that these conditions are also typically sufficient for algorithms that are not based explicitly on convex optimization [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF].

Finally, it is important to note that (a) most of the theoretical results advocate a value for the regularization parameter λ proportional to σ √ m log p, which unfortunately depends on the noise level σ (which is typically unknown in practice), and that (b) for orthogonal dictionaries, all of these assumptions are met; however, this imposes p = m.

Compressed sensing vs. statistics. Our earlier quote from [START_REF] Hocking | A Biometrics invited paper. The analysis and selection of variables in linear regression[END_REF] applies to sparse estimation as used in statistics for least-squares regression, where the dictionary D is simply the input data and x the output data. In most situations, there are some variables, represented by columns of D, that are heavily correlated. Therefore, in most practical situations, the assumptions do not apply. However, it does not mean that the high-dimensional phenomenon does not apply in a weaker sense (see the next paragraph for slow rates); moreover it is important to remark that there are other scenarios, beyond statistical variable selection, where the dictionary D may be chosen.

In particular, in signal processing, the dictionary D may be seen as measurements-that is, we want to encode α ⋆ in R p using m linear measurements Dα ⋆ in R m for m much larger than p. What the result of [START_REF] Candes | Decoding by linear programming[END_REF] alluded to earlier shows is that for random measurements, one can recover a k-sparse α ⋆ from (a potentially noisy version of) Dα ⋆ , with overwhelming probability, as long as (k log p)/m remains small. This is the core idea behind compressive sensing. See more details from [START_REF] Donoho | Compressed sensing[END_REF], Candès and Wakin [2008].

High-dimensional slow rates. While sufficient conditions presented earlier are often not met beyond random dictionaries, for the basis pursuit/Lasso formulation from (1.9), the high-dimensional phenomenon may still be observed, but only for the denoising situation and with a weaker result. Namely, as shown by [START_REF] Greenshtein | Best subset selection, persistence in high-dimensional statistical learning and optimization under ℓ 1 constraint[END_REF] and Bühlmann and Van De Geer [2011, Corollary 6.1], without assumptions regarding correlations, we have

1 n D α-Dα * 2 2 ≈ σ 2 k 2 log p m .
Note that this slower rate does not readily extend to non-convex formulations.

Impact on dictionary learning. The dictionary learning framework which we describe in this monograph relies on sparse estimation, that is, given the dictionary D, the estimation of the code α may be analyzed using the tools we have presented in this section. However, the dictionaries that are learned do not exhibit low correlations between atoms and thus theoretical results do not apply (see dedicated results in the next section). However, they suggest that (a) the codes α may not be unique in general and caution has to be observed when representing a signal x by its code α, (b) methods based on ℓ 1 -penalization are more robust as they still provably perform denoising in presence of strong correlations and (c) incoherence promoting may be used in order to obtain better-behaved dictionaries [see, e.g., [START_REF] Ramirez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF].

Theoretical results about dictionary learning

Dictionary learning, as formulated in Eq. (1.21), may be seen from several perspectives, mainly as an unsupervised learning or a matrix factorization problem. While the supervised learning problem from the previous section (sparse estimation of a single signal given the dictionary) comes with many theoretical analyses, there are still few theoretical results of the same kind for dictionary learning. In this section, we present some of them. For simplicity, we assume that we penalize with the ℓ 1 -norm and consider the minimization of min

D∈C,A∈R p×n n i=1 1 2 x i -Dα i 2 2 + λ α i 1 , (1.25)
where A = [α 1 , . . . , α n ] carries the decomposition coefficients of the signals x 1 , . . . , x n , and C is chosen as the following set:

C {D ∈ R m×p : ∀j d j 2 ≤ 1}.
Non-convex optimization problem. After imposing parsimony through the ℓ 1 -norm, given D the objective function is convex in α, given α the objective and constraints are convex in D. However, the objective function is not jointly convex, which is typical of unsupervised learning formulations. Hence, we consider an optimization problem for which it is not possible in general to guarantee that we are going to obtain the global minimum; the same applies to EM-based approaches [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] or K-means [see, e.g., [START_REF] Bishop | Pattern recognition and machine learning[END_REF].

Symmetries. Worse, the problem in Equation (1.25) exhibits several symmetries and admits multiple global optima, and the descent methods that are described in Section 5 will also have the same invariance property. For example, the columns of D and rows of A can be submitted to p! arbitrary (but consistent) permutations. There are also sign ambiguities: in fact, if (D, A) is solution of (1.25), so is (Ddiag(ε), diag(ε)A), where ε is a vector in {-1, +1} p that carries a sign pattern. Therefore, for every one of the p! possible atom orders, the dictionary learning problem admits 2 p equivalent solutions.

In other words, for a solution (D, A), the pair (DΓ, Γ -1 A) is also solution, where Γ is a generalized permutation formed by the product of a diagonal matrix with +1 and -1's on its diagonal with a permutation matrix (in particular, Γ is thus orthogonal). The fact that there are no other transformations Γ such that (DΓ, Γ -1 A) is also solution of Eq. (1.25) for all solutions (D, A) of this problem follows from a general property of isometries of the ℓ q norm for finite values of q such that q ≥ 1 and q = 2 [START_REF] Li | Isometries of the ℓ p norm[END_REF].

A manifold interpretation of sparse coding with projective geometry. The interpretation of sparse coding as a locally linear representation of a non-linear "manifold" is problematic because certain signals/features are best thought of as "points" in some space rather than vectors. For example, what does it mean to "add" two natural image patches? The simplest point structure that one can think of is affine or projective, and we show below that sparse coding indeed admits a natural interpretation in this setting, at least for normalized signals.

Indeed, Let us restrict our attention from now on to unit-norm signals, as is customary in image processing after the usual centering and normalization steps, which will be studied in Section 2.1. 7 Note that the dictionary elements d j in a solution D = [d 1 , . . . , d p ] of Eq. (1.25) also have unit norm by construction.

Let us now consider the "half sphere"

S m-1 + d ∈ S m-1 : the first non-zero coefficient of d is positive , (1.
26) where S m-1 is the unit sphere of dimension m -1 formed by the unit vectors of R m . 8 A direct consequence of the sign ambiguities of dictionary learning discussed in the previous paragraph is that, for any solution (D, A) of Eq. (1.25), there is an equivalent solution (D ′ , A ′ ) with all columns of D ′ in S m-1 + . Indeed, suppose some column d j is not in S m-1 + , and let d j [i] be its first non-zero coefficient (which is necessarily negative since d j / ∈ S m-1

+

). We can replace d j by -d j and the corresponding row of the matrix A by its opposite to construct an equivalent minimum of the dictionary learning problem in S m-1 + ×R p×n . Likewise, we can restrict the signals x i to lie in S m-1 + since replacing x i by its opposite for a given dictionary simply amounts to replacing the code α i by its opposite. Note that this identifies a patch with its "negative", but remember that the sign of its code elements is not uniquely defined in the first place in conventional dictionary learning settings (it is uniquely defined if we insist that the dictionary elements 7 Note that the fact that the individual signals are centered does not imply that the dictionary elements are.

8 Similarly, one may define the set S m-1 -by replacing "positive" by "non-negative" in (1.26). The two sets S m-1 + and S m-1 -form a partition of S m-1 with equal volume, and indeed, each one geometrically corresponds to a half sphere.

belong to S m-1

+

). This allows us to identify both the dictionary elements and the signals with points in the projective space P m-1 = P (R m ). Any k independent column vectors of D define a (k -1)-dimensional projective subspace of P m-1 (see Figure 1.14). In particular, if the data signals are assumed to be sampled from a "noisy manifold" of dimension k -1 embedded in P m-1 , an approximation of some sample x by a sparse linear combination of k elements of D can be thought of as lying in (or near) the k -1 dimensional "tangent plane" there.

Consistency results. Given the dictionary learning problem from a finite number of signals, there are several interesting theoretical questions to be answered. The first natural question is to understand the properties of the cost function that is minimized when the number of signals tends to infinity, and in particular how it converges to the expectation under the signal generating distribution [Vainsencher et al., 2011, Maurer and[START_REF] Maurer | k-dimensional coding schemes in Hilbert spaces[END_REF]. Then, given the non-convexity of the optimization problems, local consistency results may be obtained, by showing that the cost function which is minimized has a local mini-mum around the pairs (D ⋆ , A ⋆ ) that has generated the data. Given RIP-based assumptions on the dictionary D ⋆ and number of non zero elements in the columns of A ⋆ , and the noise level, [START_REF] Gribonval | Sparse and spurious: dictionary learning with noise and outliers[END_REF] show that the cost function defined in Eq. (1.25) has a local minimum around (D ⋆ , A ⋆ ) with high probability, as long as the number of signals n is greater than a constant times mp 3 . In the noiseless case, earlier results have been also obtained [START_REF] Gribonval | Dictionary identification-sparse matrixfactorization via ℓ 1 -minimization[END_REF], [START_REF] Geng | On the local correctness of ℓ 1 -minimization for dictionary learning[END_REF], and recently it has been shown that under additional assumptions, a good initializer could be found so that the previous type of local consistency results can be applied [START_REF] Agarwal | Learning sparsely used overcomplete dictionaries via alternating minimization[END_REF].

Finally, recent algorithms have emerged in the theoretical science community, which are not explicitly based on optimization [see, e.g., [START_REF] Spielman | Exact recovery of sparsely-used dictionaries[END_REF][START_REF] Recht | Factoring nonnegative matrices with linear programs[END_REF][START_REF] Arora | New algorithms for learning incoherent and overcomplete dictionaries[END_REF]. These come with global convergence guarantees (with additional assumptions regarding the signals), but their empirical performance on concrete signal and image processing problems have not yet been demonstrated.

Dictionary learning was first introduced by Olshausen andField [1996, 1997] as an unsupervised learning technique for discovering and visualizing the underlying structure of natural image patches. The results were found impressive by the scientific community, and dictionary learning gained early success before finding numerous applications in image processing [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], Mairal et al., 2008c, 2009, Yang et al., 2010a]. Without making any a priori assumption about the data except a parsimony principle, the method is able to produce dictionary elements that resemble Gabor wavelets-that is, spatially localized oriented basis functions [START_REF] Gabor | Theory of communication[END_REF][START_REF] Daugman | Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters[END_REF]], as illustrated for example in Figure 2.1.

Surprisingly, the goal of automatically learning local structures in natural images was neither originally achieved in image processing, nor in computer vision. The original motivation of Olshausen andField [1996, 1997] was in fact to establish a relation between the statistical structure of natural images and the properties of neurons from area V1 of the mammalian visual cortex. 1 The emergence of Gabor-like patterns from natural images was considered a significant result by neuroscientists. To better understand this fact, we will briefly say a few words about the importance of Gabor models within experimental studies of the visual cortex.

Since the pioneer work of [START_REF] Hubel | Receptive fields and functional architecture of monkey striate cortex[END_REF], it is known that some visual neurons are responding to particular image features, such as oriented edges. Specifically, displaying oriented bars at a particular location in a visual field's subject may elicit neuronal activity in some cells. Later, [START_REF] Daugman | Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters[END_REF] demonstrated that fitting a linear model to neuronal responses given a visual stimuli may produce filters that can be well approximated by a two-dimensional Gabor function. Models based on such filters have been subsequently widely used, and are still present in state-of-the-art predictive models of the neuronal activity in V1 [START_REF] Kay | Identifying natural images from human brain activity[END_REF][START_REF] Nishimoto | Reconstructing visual experiences from brain activity evoked by natural movies[END_REF]. By showing that Gaborlike visual patterns could be automatically obtained from the statistics of natural images, Olshausen andField [1996, 1997] provided support for the classical Gabor model for V1 cells, even though such an intriguing phenomenon does not constitute a strong evidence that the model matches real biological processes. It is indeed commonly admitted that little is known about the early visual cortex [Olshausen andField, 2005, Carandini et al., 2005] despite intensive studies and impressive results achieved by predictive models.

In this chapter, we show how to use dictionary learning for discovering latent structures in natural images, and we discuss other unsupervised learning techniques. We present them under the unified point of view of matrix factorization, which makes explicit links between different learning methods. Given a training set X = [x 1 , . . . , x n ] of signalshere, natural image patches-represented by vectors x i in R m , we wish to find an approximation Dα i = p j=1 α i [j]d j for each signal x i , where D = [d 1 , . . . , d p ] is a matrix whose columns are called "dictionary elements", and the vectors α i are decomposition coefficients. In other words, we wish to find a factorization X ≈ DA, where the matrix A = [α 1 , . . . , α n ] in R p×n carries the coefficients α i . Many unsupervised learning techniques can be cast as ma- trix factorization problems; this includes dictionary learning, principal component analysis (PCA), clustering or vector quantization [see [START_REF] Nasrabadi | Image coding using vector quantization: A review[END_REF]King, 1988, Gersho and[START_REF] Gersho | Vector quantization and signal compression[END_REF], non-negative matrix factorization (NMF) [START_REF] Paatero | Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values[END_REF]Tapper, 1994, Lee and[START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], archetypal analysis [START_REF] Cutler | Archetypal analysis[END_REF], or independent component analysis (ICA) [START_REF] Hérault | Détection de grandeurs primitives dans un message composite par une architecture de calcul neuromimétique en apprentissage non supervisé[END_REF][START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF], 1997[START_REF] Hyvärinen | Independent component analysis[END_REF]. These methods essentially differ in a priori assumptions that are made on D and A (sparse, low-rank, orthogonal, structured), and in the way the quality of the approximation X ≈ DA is measured. We will see in the rest of this chapter how these criterions influence the results obtained on natural image patches.

Pre-processing

When manipulating data, it is often important to choose an appropriate pre-processing scheme. There are several reasons for changing the way data looks like before using an unsupervised learning algorithm. For example, one may wish to reduce the amount of noise, make the data invariant to some transformation, or remove a confounding (un-wanted) factor. Unfortunately, the literature is not clear about which pre-processing step should be applied given a specific problem at hand. In this section, we describe common practices when dealing with natural image patches, such as centering, contrast normalization, or whitening. We study the effect of these procedures on the image itself, and provide visual interpretation that should help the practitioner make his choice.

Centering. The most basic pre-processing scheme consists of removing the mean intensity from every patch. Such an approach is called "centering the data" in statistics and machine learning. Borrowing some terminology from electronics, it is also called "removing the DC component" in the signal processing literature. Formally, it means applying the following update to every signal x i :

x i ← x i -   1 m m j=1 x i [j]   1 m ,
where 1 m is the vector of size m whose entries are all ones. It can also be interpreted as performing a left multiplication on the matrix X:

X ← I - 1 m 1 m 1 ⊤ m X.
When analyzing the structure of natural patches, one is often interested in retrieving geometrical visual patterns such as edges, which is a concept related to intensity variations and thus invariant to the mean intensity. Therefore, centering can be safely used; it makes the data invariant to the mean intensity, and the learned structures are expected to have zero mean as well. The effect of such a pre-processing scheme can be visualized in Figure 2.2(b) for the image man. We extract all overlapping patches from the image, remove their average pixel value, and recombine the centered patches into a new image, which is displayed in Figure 2.2(b). Because of the overlap, every pixel is contained in several patches; the pixel values in the reconstructed image are thus obtained by averaging their counterparts from the centered patches. The procedure is in fact equivalent to applying a high-pass filter to the original image. As a result, the geometrical structures are still present in the centered image, but some low frequencies have disappeared.

In practice, centering the data has been found useful in concrete applications, such as image denoising with dictionary learning [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. First removing the mean component before learning the patches structure, and then adding the mean back after processing, leads to substantially better results than working with the raw patches directly.

Variance -contrast normalization. After centering, another common practice is to normalize the signals to have unit ℓ 2 -norm, or equivalently unit variance. In computer vision, the terminology of "contrast normalization" is also often used. One motivation for such a pre-processing step is to make different regions of an image more homogeneous and to gain invariance to illumination changes. Whereas this is probably not useful for image reconstruction problems, it is a key component of many visual recognition architectures [START_REF] Pinto | Why is real-world visual object recognition hard?[END_REF][START_REF] Jarrett | What is the best multi-stage architecture for object recognition[END_REF]. However, contrast normalization should be applied with care: the naive update x i ← x i / x i 2 is likely to provide poor results for patches from uniform areas-e.g., from the sky in natural scenes. After centering, such patches do not contain any information about the scene anymore but carry residual noise, which will be greatly amplified by the naive ℓ 2 -normalization. A simple way of fixing that issue consists in providing a different treatment for patches with a small norm-say, smaller than a parameter η-and use the following update instead:

x i ← 1 max( x i 2 , η) x i .
We present in Figure 2.2(c) the effect of contrast normalization where η is chosen to be 0.2 times the mean value of x i 2 in the image, following the same patch recombination scheme as in the previous paragraph. Compared to Figure 2.2(c), areas with small intensity variations have been amplified, making the image more homogeneous. Contrast normalization can also be applied after a whitening step, which we now present.

Whitening and dimensionality reduction. Centering consists of removing from data the first-order statistics for every patch. Some early work dealing with natural images patches goes a step further by removing global second-order statistics, a procedure called whitening [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF], 1997[START_REF] Bell | The "independent components" of natural scenes are edge filters[END_REF][START_REF] Hyvärinen | Natural Image Statistics: A Probabilistic Approach to Early Computational Vision[END_REF].

Let us consider a random variable x representing centered patches uniformly sampled at random from natural images. Even though the centering step is performed at the patch level, the pixel values across centered natural images patches have zero mean µ E[x] = 0. The different pixels of a patch are indeed identically distributed, and thus the vector µ is constant; it is then necessarily equal to the constant zero vector since it has zero mean after centering. Therefore, the covariance matrix of x is the quantity Σ E[xx ⊤ ]. Whitening consists of finding a data transformation such that Σ becomes close to the identity matrix. As a result, the pixel values within a whitened patch are uncorrelated.

In practice, whitening is performed by computing the eigenvalue decomposition of the sample covariance matrix (1/n) n i=1 x i x ⊤ i = US 2 U ⊤ , where U is an orthogonal matrix in R m×m , and S in R m×m is diagonal with non-negative entries. Since the sample covariance matrix is positive semi-definite, the eigenvalues are non-negative and S = diag(s 1 , . . . , s m ) carries in fact the singular values of the matrix (1/ √ n)X. Then, whitening consists of applying the following update to every patch x i :

x i ← US † U ⊤ x i ,
where

S † = diag(s † 1 , . . . , s † m ) with s † j = 1/s j if |s j | > ε
and 0 otherwise, where ε is a small threshold to prevent arbitrarily large entries in S † . It is then easy to show that the sample covariance matrix after whitening is diagonal with k = rank(Σ) entries equal to 1. Choosing ε can be interpreted as controlling the dimensionality reduction effect of the whitening procedure. According to [START_REF] Hyvärinen | Natural Image Statistics: A Probabilistic Approach to Early Computational Vision[END_REF], such a step is important for independent component analysis. The experiments presented in this chapter do not include dimensionality reduction for dictionary learning with grayscale image patches-that is, we choose ε = 0, since the resulting whitened images do not seem to suffer from any visual artifact. The case of color patches is slightly different and is discussed in the next paragraph.

Similarly as for centering and contrast normalization, we visualize the effect of whitening on Figure 2.2(d). Removing the spatial correlation essentially results in sharper edges. As a matter of fact, whitening can be shown to amplify high frequencies and reduce low ones [START_REF] Hyvärinen | Natural Image Statistics: A Probabilistic Approach to Early Computational Vision[END_REF]. The reason is related to the nature of the principal components of natural images that form the columns of the matrix U, which will be discussed in Section 2.2. In practice, whitening is often not used in image restoration applications based on dictionary learning since it modifies the nature of the noise.

Pre-processing color patches. Finally, we conclude this section on pre-processing by discussing the treatment of color image patches. We assume that images are encoded in the RGB space, and that patches are obtained by concatenating information from the three color channels. In other words, l × l image patches are represented by vectors x i of size m = 3 × l × l.

The main question about color processing is probably whether or not the RGB color space should be used [START_REF] Pratt | Spatial transform coding of color images[END_REF][START_REF] Faugeras | Digital color image processing within the framework of a human visual model[END_REF][START_REF] Sharma | Digital color imaging[END_REF]. The choice of the three channels R, G, and B directly originates from our first understanding of the nature of color. After discovering the existence of the color spectrum, Newton [1675] already suggests that the rays of light "impinging" on the retina will "affect the sense with various colours, according to their bigness and mixture; the biggest with the strongest colours, reds and yellow; the least with the weakest, blue and violets; the middle with green".

Building upon Newton's findings, it appears that [START_REF] Young | A course of lectures on natural philosophy and the mechanical art[END_REF] was the first to introduce the concept of trichromatic vision. According to [START_REF] Maxwell | On the theory of compound colours, and the relations of the colours of the spectrum[END_REF]:

"Young appears to have originated the theory, that the three elements of colour are determined as much by the constitution of the sense of sight as by anything external to us. He conceives that three different sensations may be excited by light. . . He conjectures that these primary sensations correspond to red, green, and violet."

Finally, von [START_REF] Helmholtz | LXXXI. on the theory of compound colours[END_REF], [START_REF] Grassmann | LXXXVII. On the theory of compound colours[END_REF], and [START_REF] Maxwell | On the theory of compound colours, and the relations of the colours of the spectrum[END_REF] proposed rigorous rules of color compositions, paving the way to the modern treatment of color in vision processing [see Forsyth and Ponce, 2012, Chapter 3]. More than a century after these early discoveries, the existence in the eye of three types of biological photoreceptorsrespectively corresponding to red, green, and blue wavelengths-is established [START_REF] Nathans | Molecular genetics of human color vision: the genes encoding blue, green, and red pigments[END_REF], and the RGB color space is widely used in electronic devices.

The question of changing the color space can be rephrased as follows: should we apply a linear or non-linear transformation f to each RGB pixel value: (u, v, w) = f (r, g, b) and work in the resulting space? One of the early motivation for studying the representation of color was to reduce the bandwidth required by RGB over communication channels. RGB components are indeed highly correlated; thus, removing the redundant information between the different channels allows more efficient coding for transmitting information [START_REF] Pratt | Spatial transform coding of color images[END_REF]]. Another motivation is that the Euclidean distance on RGB is known to be a poor estimation of the perceptual distance by humans. Therefore, a large variety of color spaces have been designed, such as CIELab, CIEXYZ, YIQ, or YCrBr [see [START_REF] Sharma | Digital color imaging[END_REF], which are partially improving upon RGB regarding the decorrelation of the transformed channels and human perception. However, it remains unclear whether these spaces should be used or not in practice, and the optimal choice depends on the task and algorithm. For instance, for color image denoising, [START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF] and Dabov et al. [2007b] choose the YCrBr space, whereas [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], Mairal et al. [2008c] and [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF] successfully use RGB. Changing the color space modifies indeed the noise structure as a side effect, which may be problematic for some methods.

The other pre-processing steps that we have presented for grayscale images can also be applied to color image patches. The first difference with the monochromatic setting lies in the centering scheme when work-ing with RGB. One motivation for centering grayscale image patches is to learn geometric structures that are invariant to the mean intensity. For color images, the quantity 1 m m j=1 x i [j] unfortunately lacks of interpretation, and removing it from a patch does not seem to achieve any interesting property. Instead, we may be looking for geometrical structures in images that are invariant to color patterns, and thus, one may remove the mean color from the patch. In other words, one should center every R,G,B channel independently. As a result, the centering scheme becomes also invariant to any linear transformation of the color space, which can be a desirable property. The effect of such a preprocessing can be visualized in Figure 2.3(b) for the image kodim07 from the Kodak PhotoCD dataset.2 Note that centered images have negative values and thus images are shifted and rescaled for visualization purposes. Regions with relatively uniform colors appear grayish, and thus mostly contain geometrical content; some other areas contain transitions between a color-say, represented by RGB values (r, g, b)and its opponent one-represented by (-r, -g, -b). The colored halos around the pink flowers and the green leaves are thus due to this color/opponent color transitions.

In the experiments of this paper about dictionary learning, whitening comprises a light dimensionality reduction step, where we threshold to zero the smallest singular values of X (the entries on the diagonal of S). More precisely, we keep the k largest singular values s 1 , . . . , s k such that k i=1 s 2 i ≥ 0.995 m i=1 s 2 i ; in other words, we keep 99.5% of the explained variance of the data. Contrast normalization can then be applied without particular attention to the nature of the patch. The effects of whitening and contrast normalization are displayed on Figures 2.3(c) and 2.3(d), respectively. Similar conclusions can be drawn for color image patches as in the monochromatic case.

Principal component analysis

Principal component analysis (PCA), also known as the Karhunen-Loève or Hotelling transform [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF], is probably the most widely used unsupervised data analysis technique. Even though it is often presented as an iterative process finding orthogonal directions maximizing variance in the data, it can be cast as a low-rank matrix factorization problem:

min U∈R m×k ,V∈R n×k X -UV ⊤ 2 F s.t. U ⊤ U = I k ,
where k is the number of principal components we wish to obtain, and I k is the identity matrix in R k×k . We also assume that the rows of the matrix X have zero mean. As a consequence of the theorem of [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF], the matrix U contains the principal components of X corresponding to the k largest singular values. In Figure 2.4(b), we visualize the principal components of n = 400 000 natural image patches of size 16 × 16 pixels, ordered by largest to smallest singular value (form left to right, then top to bottom). The components resemble Fourier basis-that is, product of sinusoids with different frequencies and phases, similarly to the discrete cosine transform (DCT) dictionary [START_REF] Ahmed | Discrete cosine transform[END_REF] presented in Figure 2.4(a). However, there is a good reason for obtaining sinusoids that is simply related to a property of translation invariance, meaning that the patterns observed in Figure 2.4(b) are unrelated to the underlying structure of natural images. This fact is well known and has been pointed earlier by others in the literature [see, e.g., [START_REF] Bossomaier | Why spatial frequency processing in the visual cortex?[END_REF][START_REF] Field | Relations between the statistics of natural images and the response properties of cortical cells[END_REF][START_REF] Simoncelli | Natural image statistics and neural representation[END_REF][START_REF] Hyvärinen | Natural Image Statistics: A Probabilistic Approach to Early Computational Vision[END_REF]. Second-order statistics of natural images patches are commonly assumed to be invariant by translation, such that the correlation of pixel values at locations z 1 = (k 1 , l 1 ) and z 2 = (k 2 , l 2 ) only depends on the displacement z 1z 2 [START_REF] Simoncelli | Natural image statistics and neural representation[END_REF]. This is particularly true for patches that are extracted from larger images at any arbitrary position, and whose distribution is exactly translation invariant. Such signals are called "stationary" and their principal components (equivalently the eigenvectors of the covariance matrix) are often considered to be well approximated by the Discrete Fourier Transform (DFT), as noted by [START_REF] Pearl | On coding and filtering stationary signals by discrete fourier transforms (corresp.)[END_REF], leading to sinusoidal principal components.

Nevertheless, the relation between principal components and si-nusoids only holds rigorously in an asymptotic regime. For instance, consider the case of an infinite one-dimensional signal with covariance Σ[k, l] = σ(kl) for positions k and l, where σ is an even function.

Then, for all frequency ω and phase ϕ, ωk+ϕ) ,

l Σ(k, l)e i(ωl+ϕ) = l σ(l -k)e i(ωl+ϕ) = l ′ σ(l ′ )e iωl ′ e i(
where i denotes the imaginary unit, and the sums are over all integers.

Since the function σ is even, the infinite sum l ′ σ(l ′ )e iωl ′ is real, and the signals [sin(ωk + ϕ)] k∈Z are all eigenvectors of the covariance operator Σ. Equivalently, they are principal components of the input data. The case of two-dimensional signals can be treated similarly but it involves heavier notation. Drawing conclusions about the eigenvectors of a finite covariance matrix computed on natural image patches is nevertheless subject to discussion, and understanding the quality of the approximation of the principal components by the DFT in the finite regime is non-trivial [START_REF] Pearl | On coding and filtering stationary signals by discrete fourier transforms (corresp.)[END_REF].

From an experimental point of view, the fact that the structure of natural images has nothing to do with the patterns displayed in Figure 2.4(b) is easy to confirm. In Figure 2.5, we present principal components computed on all overlapping patches from the image tiger. Even though this image is sketched by hand-and thus, is not naturalwe recover sine waves as in Figure 2.4(b).

Clustering or vector quantization

Clustering techniques have been used for a long time on natural image patches for compression and communication purposes under the name of "vector quantization" [START_REF] Nasrabadi | Image coding using vector quantization: A review[END_REF]King, 1988, Gersho and[START_REF] Gersho | Vector quantization and signal compression[END_REF]. The goal is to find p clusters in the data, by minimizing the following objective: min

D∈R m×p ∀i, l i ∈{1,...,p} n i=1 x i -d l i 2 2 , (2.1)
where the columns of D = [d 1 , . . . , d p ] are called "centroids" and l i is the index of the cluster associated to the data point x i . The al- gorithm K-means [see [START_REF] Hastie | The elements of statistical learning[END_REF] approximately optimizes the non-convex objective (2.1) by alternatively performing exact minimization with respect to the labels l i with D fixed, and with respect to D with the labels fixed.

To make the link between clustering and matrix factorization, it is also possible to reformulate (2.1) as follows min

D∈R m×p A∈{0,1} p×n 1 n n i=1 1 2 x i -Dα i 2 2 s.t. ∀i, p j=1 α i [j] = 1,
where A = [α 1 , . . . , α n ] carries binary vectors that sum to one, and clustering can be subsequently seen as a matrix factorization problem: min

D∈R m×p A∈{0,1} p×n 1 2n X -DA 2 F s.t. ∀i, p j=1 α i [j] = 1.
Then, the algorithm K-means is performing alternate minimization between A and D, each step decreasing the value of the objective. We visualize clustering results in Figure 2.6 after applying 250 iterations of the K-means algorithm on n = 400 000 natural image patches and choosing p = 256 centroids. Without whitening, K-means produces mostly low-frequency patterns, which is not surprising since most of the energy (ℓ 2 -norm) of images is concentrated in low frequencies. It is indeed known that the spectral power of natural images obtained in the Fourier domain typically decreases according to the power law 1/f 2 for a frequency f [see [START_REF] Simoncelli | Natural image statistics and neural representation[END_REF]. After whitening, high-frequency "Gabor-like" patterns and checkerboard patterns emerge from data. It is thus interesting to see that with an appropriate pre-processing procedure, simple unsupervised learning techniques such as K-means are able to discover Gabor features with different orientations, frequencies, and positions within patches. We note that the localized checkerboard patterns might be an artifact of the whitening procedure without dimensionality reduction, which amplifies such patterns appearing among the last principal components presented in

Dictionary learning

We are now interested in the dictionary learning formulation originally introduced by Olshausen andField [1996, 1997] and its application to natural image patches. We consider the corresponding matrix factorization formulation (1.22), which we recall here min

D∈C,A∈R p×n 1 2 X -DA 2 F + λΨ(A), (2.2)
where Ψ(A) = n i=1 ψ(α i ) and C is the set of matrix whose columns have a Euclidean norm smaller than one.

In Figure 2.7, we present visual results obtained with such a formulation where ψ is the ℓ 1 -norm. As in the previous section about clustering, we use n = 400 000 natural image patches, either gray, or extracted from color images in RGB. We center the patches, according to the procedure presented in Section 2.1, and rescale the matrix X such that its columns have unit norm on average. Then, we learn a dictionary D with p = 256 dictionary elements and set the regularization parameter λ to 0.1. We use the online dictionary learning algorithm of Mairal et al. [2010a] available in the software SPAMS, making 10 passes over the data. The dictionaries obtained for gray and RGB patches and with two different preprocessing steps are displayed in Figure 2.7. Note that other dictionary learning approaches can be used for this task and providing similar results. This is for instance the case of the K-SVD algorithm of [START_REF] Aharon | An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], which we present in Section 5.5.

When processing grayscale image patches, a simple centering step produces dictionary elements with both low-frequency elements, and high-frequency, localized, Gabor-like patterns, as shown in Figure 2.7(a). The results after a whitening step are displayed in Figure 2.7(c). Whitening corresponds to applying a high-pass filter to the original images, and we obtain exclusively Gabor-like features as in the early work of Olshausen andField [1996, 1997]. The case of color image patches presented in Figures 2.7(b) and 2.7(d) is also interesting. After centering, about 230 dictionary elements out of 256 are grayscale, or are almost grayscale, which is consistent with the greyish appearance of images after centering (see Section 2.1). An intriguing phenomenon concerns the remaining colored dictionary elements. Those exhibit lowfrequency patterns with two opponent colors, blue and yellow for 10 of them, green and red/purple for 16 of them. Interestingly, such opponent color patterns are typical characteristics of neurons' receptive fields described in the neuroscience literature [START_REF] Livingstone | Anatomy and physiology of a color system in the primate visual cortex[END_REF]Hubel, 1984, Ts'o and[START_REF] Ts | The organization of chromatic and spatial interactions in the primate striate cortex[END_REF].

Structured dictionary learning

Because the ℓ 1 -norm cannot model interactions between dictionary elements, it is natural to consider extensions of dictionary learning where a particular structure is taken into account. The concept of structured sparsity presented in Section 1.3 is a natural tool for this; it consists of replacing the ℓ 1 -norm by a more complex sparsity-inducing penalty. Jenatton et al. [2010aJenatton et al. [ , 2011b]], when a pre-defined hierarchical structure is assumed to exist among the p dictionary elements. The formulation (2.2) is considered where the function Ψ is the hierarchical Group-Lasso penalty of [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF]. A tree is given, e.g., is defined by the user, and one dictionary element is associated to every node of the tree. The penalty constructs a group structure G of subsets of {1, . . . , p}, each group containing one node and all its descendants in the tree. An example is illustrated in Figure 1.12. The resulting penalty Ψ is then defined as

Hierarchical dictionary learning. A first extension has been investigated by

Ψ(A) = n i=1 g∈G α i [g] q ,
(2.3)

where . q can either represent the ℓ ∞ -or the ℓ 2 -norm. As explained in Section 1.3, an effect of the penalty is that a dictionary element can be used in the decomposition of a patch only if its parent in the tree is also used. We present some visual results for q = ∞ in Figure 2.8, after centering the natural image patches, for some manually tuned regularization parameter λ, and two different tree structures. Dictionary elements naturally organize themselves in the tree, often with low frequencies near the root of the tree, and high frequencies near the leaves. We also observe strong correlations between each parent node and their children in the tree, where children often look like their parent, but with higher frequencies and with minor variations.

Topographic dictionary learning. A two-dimensional grid structure has also been used for learning dictionaries [START_REF] Kavukcuoglu | Learning invariant features through topographic filter maps[END_REF], Mairal et al., 2011], which was directly inspired from the topographic independent component analysis formulation of [START_REF] Hyvärinen | Topographic independent component analysis[END_REF].

The principle is similar to the hierarchical case-that is, we can use a penalty Ψ as in (2.3), but the set of groups G takes into account a different structure. We assume here that the dictionary elements are organized on a grid, such that we can define neighborhood relations between them.

For instance, we can organize the p dictionary elements on a √ p× √ p grid, and consider p overlapping groups that are 3 × 3 or 4 × 4 spatial neighborhoods on the grid (to avoid boundary effects, we assume the grid to be cyclic). We represent some results obtained with such a formulation in Figure 2.9, where a function Ψ is defined in (2.3) with q = 2. The regularization parameter λ is also manually tuned. The dictionary elements automatically organize themselves on the grid, and exhibit some intriguing spatial smoothness. Another formulation achieving a similar effect was also proposed by [START_REF] Garrigues | Group sparse coding with a Laplacian scale mixture prior[END_REF], by mixing sparse coding with a probabilistic model involving latent variables for groups of neighbor variables. Another approach was also proposed by [START_REF] Gregor | Structured sparse coding via lateral inhibition[END_REF], which models inhibition effects between dictionary elements.

Other matrix factorization methods

We have focused our study of learning methods for natural image patches on sparse coding and clustering techniques so far. Other matrix factorization formulations have shown to be also effective for this task.

In particular, independent component analysis (ICA) was applied successfully to image patches around the same time as dictionary learning by [START_REF] Bell | The "independent components" of natural scenes are edge filters[END_REF], who also obtained Gabor-like patterns from whitened data.

Independent component analysis.

A significantly different point of view than dictionary learning consists of factorizing whitened data X as a product DA, where the sparsity principle is replaced by an assumption of statistical independence. Therefore, the ICA approach requires by nature a probabilistic interpretation of the data. We assume that x in R m is a random variable representing a signal-here, a natural image patch-and the columns of X = [x 1 , . . . , x n ] are random realizations of the variable x-that is, a set of natural image patches selected uniformly at random. The principle of ICA assumes that there exists a factorization x = Dα, where D is an orthogonal matrix, and α = D ⊤ x is a random vector whose entries are statistically independent [see [START_REF] Bell | The "independent components" of natural scenes are edge filters[END_REF][START_REF] Hyvärinen | Independent component analysis[END_REF], 2009].

Then, the columns of D are called independent components. Given data, the concept of "statistical independence" alone does not immediately yield a precise cost function to optimize, and thus numerous variants of ICA exist in the literature. Ultimately, most variants of ICA try to approximate the statistical independence assumption, by finding an orthogonal transformation D ⊤ such that the probability density p(α), assuming it exists, factorizes into the product of its marginals p j=1 p(α[j]). A popular cost function to compare probability distributions is the Kullback-Leibler distance [see [START_REF] Cover | Elements of information theory[END_REF], denoted by KL, and a natural goal of ICA is to minimize over

D KL   p(α), p j=1 p(α[j])   R p p(α) log p(α) p j=1 p(α[j]) dα,
which is equal to zero if and only if the α[j]'s are independent. The quantity KL above is also called the mutual information between the variables α[j], and can also be simply written in terms of the entropy function H:

3 KL   p(α), p j=1 p(α[j])   = p j=1 H(α[j]) -H(α).
Interestingly, the entropy H(α) can be shown to be independent of D when D is orthogonal and x is whitened with identity covariance, such that the goal of ICA can be expressed as the minimization of the sum of entropies

p j=1 H(α[j]) = p j=1 H(d ⊤ j x).
(2.4)

Unfortunately, entropy remains an abstract quantity that is not computable, and thus (2.4) is not yet an objective function that is easy to minimize. Turning (2.4) into a concrete formulation and algorithm can be achieved by following different strategies. One of them consists of parameterizing the densities p(d ⊤ j x), assuming they admit a known parametric form, and minimizing (2.4) becomes equivalent to performing maximum likelihood estimation [see [START_REF] Hyvärinen | Independent component analysis[END_REF]]. Another approach consists of using non-parametric estimators of the entropy [START_REF] Pham | Fast algorithms for mutual information based independent component analysis[END_REF]. Finally, a large class of methods consists of using approximations of the entropy cost function (2.4) or encourage the distributions of the α[j]'s to be "non-Gaussian" [START_REF] Cardoso | Dependence, correlation and Gaussianity in independent component analysis[END_REF]. The non-gaussianity principle might not seem intuitive at first sight, but it is in fact rather natural when considering a minimization entropy problem. Among all probability distributions with same variance, the Gaussian ones are known to maximize entropy [START_REF] Cover | Elements of information theory[END_REF]]. As such, minimizing (2.4) implies that we are looking for some non-Gaussian random variables α[j].

Extensions where D is not orthogonal and when the data X is not whitened exist, but we have omitted them in this brief presentation for simplicity [see [START_REF] Hyvärinen | Independent component analysis[END_REF], for a better overview].

A remarkable fact is that ICA applied to natural image patches yields very similar results to the sparse coding model of Olshausen andField [1996, 1997], as shown by [START_REF] Bell | The "independent components" of natural scenes are edge filters[END_REF]. In fact, an equivalence exists between the two formulations under particular asymptotic conditions [see the appendix of [START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by V1?[END_REF].

Non-negative matrix factorization. Another popular unsupervised learning technique is the non-negative matrix factorization (NMF) [START_REF] Paatero | Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values[END_REF], which was shown to be able to automatically discover some interpretable features on datasets of human faces [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. When the data matrix X is non-negative, the method consists of finding a product DA that approximates X, and where each factor is also non-negative. The corresponding formulation can be written as

min D∈R m×p ,A∈R p×n X -DA 2 F s.t. D ≥ 0 and A ≥ 0.
Even though this technique could be applied to natural image patches without pre-processing, it has been empirically observed that NMF does not yield any interpretable feature when applied to such data [Mairal et al., 2010a]. Note that NMF is often used with a different loss function than the square loss, in particular for audio applications [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: with application to music analysis[END_REF].

Archetypal analysis. Around the same time as dictionary learning and non-negative matrix factorization, archetypal analysis was introduced by [START_REF] Cutler | Archetypal analysis[END_REF] for discovering latent factors from high-dimensional data. The main motivation was to propose an unsupervised learning technique that is more interpretable than principal component analysis. The method seeks a factorization X = DA, with two symmetrical geometrical constraints. First, each column d j of D, called "archetype" is forced to be a convex combination of a few data points. In other words, for all j ∈ {1, . . . , p}, there exists a vector β j such that d j = Xβ j , and β j is in the simplex ∆ n defined as

∆ n β ∈ R n s.t. β ≥ 0 and n i=1 β[i] = 1 . (2.5)
Second, each data point x i is encouraged to be close to the convex hull of the archetypes, meaning that x i should be close to a product Dα i , where α i is in the simplex ∆ p , defined as in (2.5) when replacing n by p. The resulting formulation is the following min

α i ∈∆p for 1≤i≤n β j ∈∆n for 1≤j≤p n i=1 x i -Dα i 2 2 s.t. ∀j, d j = Xβ j ,
which, in the matrix factorization form, is equivalent to

min α i ∈∆p for 1≤i≤n β j ∈∆n for 1≤j≤p X -XBA 2 F , where A = [α 1 , . . . , α n ], B = [β 1 , . . . , β p ]
and the matrix of archetypes D is equal to the product XB.

Archetypal analysis is closely related to the dictionary learning formulation of Olshausen andField [1996, 1997] since the simplicial constraint α i ∈ ∆ p can be equivalently rewritten as the ℓ 1 -norm constraint α i 1 = 1 associated to a non-negativity one α i ≥ 0. As a result, the coefficients α i are sparse in practice and archetypal analysis provides a sparse decomposition of the data points x i . The main difference with dictionary learning is the set of constraints d j = Xβ j . Because the vectors β j are constrained to be in the simplex ∆ n , they are encouraged to be sparse and each archetype is a convex combination of a few data points only. Such a relation between latent factors d j and the data X is useful whenever interpreting D is important, e.g., in experimental sciences. For example, clustering techniques provide such associations between data and centroids. It is indeed common in genomics to cluster gene expression data from several individuals, and to interpret each centroid by looking for some common physiological traits among individuals of the same cluster [START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF]].

Archetypal analysis is also related to non-negative matrix factorization [START_REF] Paatero | Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values[END_REF]. When the data X is non-negative, it is easy to see that both the matrix A and the matrix D of archetypes are also non-negative. Unfortunately, archetypal analysis did not encountered as much success as dictionary learning or NMF, despite the fact that it provides an elegant methodology for interpreting its output.

One of the reason for this lack of popularity may be that no efficient software has been available for a long time, which may have limited its application to important scientific problems. Based upon such observations, [START_REF] Chen | Fast and robust archetypal analysis for representation learning[END_REF] have recently revisited archetypal analysis for computer vision, with the goal of bringing back this powerful unsupervised learning technique into favor. They made publicly available an efficient implementation of archetypal analysis in the SPAMS software, and demonstrated that it could perform as well as dictionary learning for some classification tasks, while offering natural mechanisms for visualizing large databases of images (see Section 4.6).

Bayesian models for dictionary learning. Early models of dictionary learning were probabilistic [see, for instance, [START_REF] Lewicki | Efficient coding of natural sounds[END_REF][START_REF] Lewicki | Probabilistic framework for the adaptation and comparison of image codes[END_REF][START_REF] Lewicki | Learning overcomplete representations[END_REF]. Each image patch x is considered to be a random variable with a normal distribution given a dictionary D and a latent variable α. More precisely, we have the conditional law

p(x|D, α) ∝ e -1 2σ 2 x-Dα 2 2 ,
(2.6) and the prior distribution on the coefficients α is Laplacian:

4 p(α) ∝ e -λ α 1 . (2.7)
With this probabilistic model, maximizing the posterior p(α|x, D) given some observation x and a dictionary D, or equivalently minimizing the negative log posterior, yields the Lasso formulation. Indeed, by using Bayes' rule, we have

-log p(α|x, D) = -log p(x|D, α) -log p(α) = 1 2σ 2 x -Dα 2 2 + λ α 1 .
(2.8)

When n natural image patches x 1 , . . . , x n are observed, learning a dictionary D can be achieved by (i) modeling each x i according to (2.6) with a latent variable α i associated to x i , (ii) assuming the α i 's to be statistically independent one from each other, (iii) choosing a prior distribution p(D), e.g., uniform on the set of matrices C defined in (2.2), and (iv) computing the point estimate min

D,α 1 ,...,αn -log p(D, α 1 , . . . , α n |x 1 , . . . , x n ).
By using again Bayes' rule and by using the statistical independence assumption, we recover the classical "matrix factorization" formulation presented earlier:

-log p(D, A|X) = -log p(α 1 , . . . , α n |X, D) -log p(D) = n i=1 -log p(α i |x i , D) -log p(D) = n i=1 1 2σ 2 x i -Dα i 2 2 + λ α i 1 -log p(D),
(2.9)

where X = [x 1 , . . . , x n ] and A = [α 1 , . . . , α n ]. Simply maximizing the posterior distribution is however not satisfactory from the point of view of a Bayesian statistician, who is usually not interested in point estimates, but in the full posterior distribution [see [START_REF] Bayarri | The interplay of bayesian and frequentist analysis[END_REF]. 5 In fact, the formulation (2.9) discards any uncertainty regarding the model parameters D and A and would be called a "frequentist" approach in statistics: assuming that one can repeatedly draw natural image patches at random, one wishes to find a dictionary D that is good on average, which is nothing else than the (non-probabilistic) empirical risk minimization point of view presented in Section 1.4. A first step towards a Bayesian dictionary learning formulation consists of integrating instead of maximizing with respect to the latent variables α i , as typically done in Bayesian sparse linear mod-els [START_REF] Park | The Bayesian Lasso[END_REF]Casella, 2008, Seeger, 2008], and maximizing the posterior p(D|X). As a result, we need to minimize

-log p(D|X) = -log p(X|D) -log p(D) = n i=1 -log p(x i |D) -log p(D) = n i=1 -log α i ∈R p p(x i , α i |D)dα i -log p(D) = n i=1 -log α i ∈R p p(x i |α i , D)p(α i )dα i -log p(D).
(2.10)

The formulation is again a pointwise estimator for the dictionary D but it can be argued to be more "Bayesian" due to the treatment of the latent variables α i . Another possibility is to maximize the data likelihood p(X|D), resulting in the same formulation as (2.10) without the prior termlog p(D); this is the strategy adopted for instance by [START_REF] Lewicki | Probabilistic framework for the adaptation and comparison of image codes[END_REF]. Note that optimizing (2.10) is difficult since the integrals do not admit an analytical closed form, and some approximations have to be made [START_REF] Lewicki | Probabilistic framework for the adaptation and comparison of image codes[END_REF].

Recently, a fully Bayesian dictionary learning formulation has been proposed by [START_REF] Zhou | Nonparametric bayesian dictionary learning for sparse image representations[END_REF][START_REF] Zhou | Nonparametric Bayesian dictionary learning for analysis of noisy and incomplete images[END_REF]. By "fully Bayesian", we mean that all model parameters such as dictionary D, coefficients α, but also hyper-parameters such as σ are modeled with probability distributions, and the full posterior density is estimated by using Gibbs sampling. The model is in fact significantly different than the one described by (2.6) and (2.7). In a nutshell, it involves a Beta-Bernoulli process for selecting the non-zero coefficients in α, Gaussian priors with Gamma hyperpriors for the value of these coefficients and for the dictionary elements. Compared to the traditional matrix factorization formulation (2.2), a Bayesian treatment has both advantages and drawbacks, which often appear in passionate discussions between Bayesians and frequentists: on the one hand, the Bayesian formulation is robust to model misspecification and can learn some hyper-parameters such as the noise level σ; on the other hand, Bayesian inference is significantly more involved and computationally costly than matrix factorization.

Convolutional sparse coding. The dictionary learning formulation of Olshausen andField [1996, 1997] is appropriate for processing natural image patches. A natural extension to full images instead of patches is called "convolutional sparse coding". It consists of linearly decomposing an image by using small dictionary elements placed at all possible locations in the image. Such a principle was described early by [START_REF] Zhu | What are textons?[END_REF] without concrete application. It has been revisited recently, and has been shown useful for various recognition tasks [START_REF] Zeiler | Deconvolutional networks[END_REF], 2011[START_REF] Rigamonti | Are sparse representations really relevant for image classification[END_REF], 2013, Kavukcuoglu et al., 2010a] Specifically, let us consider an image x represented by a vector of size l and let the binary matrix R k in {0, 1} m×l be the linear operator such that R k x is the patch of size √ m× √ m from x centered at the pixel indexed by k. The purpose of R k is to "extract" a patch at a specific location, but it is also easy to show that its adjoint operator R ⊤ k is such that R ⊤ k z, for some vector z in R m , is a vector of size l representing an image with zeroes everywhere except the patch k that contains the m entries of z. In other words, R ⊤ k positions at pixel k a small patch in the larger image. To simplify, we assume that there is no boundary effect when the pixel k is close to the image border, e.g., we use zeropadding to define patches that overlap with the image boundary. Then, it becomes natural to decompose x with the Lasso formulation:

min A∈R p×l 1 2 x - l k=1 R ⊤ k Dα k 2 2 + λ l k=1 α k 1 , (2.11)
where D is a dictionary in R m×p , and, as usual, A = [α 1 , . Given now a collection of n images x 1 , . . . , x n , which are assumed to be of the same size for simplicity, the dictionary D can be learned by minimizing min

{A i ∈R p×l } i=1...n D∈C 1 n n i=1   1 2 x i - l k=1 R ⊤ k Dα i,k 2 2 + λ l k=1 α i,k 1   ,
(2.12) where A i = [α i,1 , . . . , α i,k ] for all i = 1, . . . , n. A natural optimization scheme for addressing (2.12) is to alternate between the minimization of D with A fixed and vice versa, as often done for the classical dictionary learning problem (see Section 5). Even though finding the global optimum of (2.12) is not feasible because the objective function is nonconvex, alternate minimization provides a stationary point [START_REF] Bertsekas | Nonlinear programming[END_REF]. Updating the dictionary D with fixed coefficients A can be achieved by projected gradient descent.

In Figure 2.10, we visualize two dictionaries of size p = 100 with elements of size m = 16 × 16 pixels. We learned them on 30 whitened natural images that are rescaled such that each √ m × √ m patch has unit ℓ 2 -norm on average. We perform 300 steps of alternate minimization. For updating A or D, we always use 5 iterations of the algorithm ISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], yielding a small decrease of the objective function at each step. Interestingly, we observe that the learned features have the following properties: (i) in general, they are well centered; (ii) some of them are more complex than the ones obtained with the classical dictionary learning formulation. In other words, they go beyond simple Gabor features, with curvy patterns, corners, and blobs.

Discussion

In this section, we have reviewed a large number of unsupervised learning techniques that are able to discover underlying structures in natural images. We have focused on matrix factorization and sparse coding, but other techniques that are out of the scope of our study are known to achieve similar results. We briefly mention a few of them.

A successful approach is the Gaussian mixture model (GMM). Even though it is classical, it has been shown only recently to be well adapted to represent natural image patches. GMMs for image patches were first investigated by [START_REF] Yu | Solving inverse problems with piecewise linear estimators: from Gaussian mixture models to structured sparsity[END_REF] and then further developed by [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], yielding very good results for several image restoration tasks. The probabilistic model consists of representing the distribution of a patch x as a convex combination of p Gaussian distributions with means µ k and (often diagonal) covariance matrices Σ k for k = 1, . . . , p. Given a database of patches x 1 , . . . , x n , the EM-algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is typically used to learn means and covariances. When displaying the means, small localized Gabor filters typically appear when the data is whitened [START_REF] Coates | An analysis of single-layer networks in unsupervised feature learning[END_REF].

Finally, small oriented localized filters can also be learned from natural image patches by using undirected graphical models such as restricted Boltzmann machines (RBM) [see [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF][START_REF] Bengio | Learning deep architectures for AI[END_REF][START_REF] Coates | An analysis of single-layer networks in unsupervised feature learning[END_REF], and also by using convolutional neural networks [LeCun et al., 1998a, Zeiler and[START_REF] Zeiler | Visualizing and understanding convolutional neural networks[END_REF], which have recently gained some popularity for solving large-scale visual recognition tasks.

We have previously shown that natural image patches could be modeled by using dictionary learning techniques and its variants, but we have not presented any concrete application yet. We have indeed focused so far on the interpretation of the visual patterns obtained by various methods, but one may wonder whether or not these interpretable structures can be useful for prediction tasks. For quite a long time, the answer to this question was uncertain, until the work of [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] on image denoising that achieved state-of-the-art results compared to other approaches at that time.

The current section is devoted to several applications of dictionary learning in image processing. First, we present the simple and yet effective scheme for image denoising introduced by Elad and Aharon [2006], and then move to more complex tasks. In general, we show that dictionary learning can be useful for predicting missing visual information, leading to state-of-the-art results for image inpainting and demosaicking [Mairal et al., 2008c[START_REF] Mairal | Non-local sparse models for image restoration[END_REF], super-resolution and deblurring [Yang et al., 2010a, 2012a[START_REF] Couzinie-Devy | Dictionary learning for deblurring and digital zoom[END_REF], Dong et al., 2011b[START_REF] Zeyde | On single image scale-up using sparse representations[END_REF], Wang et al., 2012], face compression [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF], or for inverting non-linear local trans-formations [Mairal et al., 2012]. We also present natural extensions of dictionary learning to video processing [START_REF] Protter | Image sequence denoising via sparse and redundant representations[END_REF], and finally, we conclude this section with a presentation of other patch-modeling approaches [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Awate | Unsupervised, information-theoretic, adaptive image filtering for image restoration[END_REF], Dabov et al., 2007a[START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], Yu et al., 2012[START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF].

A small part of the material of this section is borrowed from the PhD thesis of the first author [START_REF]Sparse coding for machine learning, image processing and computer vision[END_REF].

Image denoising

Let us consider the classical problem consisting of restoring a noisy image y in R n that has been corrupted by white Gaussian noise with known standard deviation σ. 1 In many cases, image denoising is formulated with an energy to minimize [e.g., [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]:

min x∈R n 1 2 y -x 2 2 + λψ(x), (3.1)
where ψ : R n → R is a regularization function, and the quadratic "data-fitting" term ensures that the estimate x is close to the noisy observation y. When (3.1) is derived from a probabilistic image model, it is often interpreted from the point of view of maximum a posteriori estimation, where λψ(x) is related to some negative log prior distribution on x. Even though (3.1) looks simple at first sight, finding a good regularization function ψ is difficult, and in fact, it is probably one of the most important research topic in image processing nowadays. In early work, various smoothness assumptions about x have led to different functions ψ. Specifically, natural images were assumed to have small total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], or the smoothness between adjacent pixel values was modeled with Markov random fields (MRF) [START_REF] Zhu | Prior learning and gibbs reaction-diffusion[END_REF].

As described in Section 1.2, sparse image models based on wavelets have also been popular for the denoising task [see [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. This line of work has inspired the method we present in this section but the latter differs in two aspects: (i) the approach of [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] is patch-based like other "modern" image processing methods [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], Dabov et al., 2007a[START_REF] Roth | Fields of experts[END_REF];

(ii) it adapts to the image patches with dictionary learning, and thus it does not use any pre-defined wavelet basis.

Patch denoising given a fixed dictionary. [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] have proposed a simple denoising procedure that treats every patch independently. Let us consider the √ m × √ m patch y i of the noisy image y, centered at the pixel indexed by i, and assume that an appropriate dictionary D in R m×p is given.2 Then, the patch y i is denoised by the following steps:

1. center y i (see Section 2.1),

y c i y i -µ i 1 m with µ i 1 n 1 ⊤ m y i ;
2. find a sparse linear combination of dictionary elements that approximates y c i up to the noise level:

min α i ∈R p α i 0 s.t. y c i -Dα i 2 2 ≤ ε, (3.2)
where ε is proportional to the noise variance σ 2 ;

3. add back the mean component to obtain the clean estimate xi :

xi Dα ⋆ i + µ i 1 m ,
where α ⋆ i is the solution, or approximate solution, obtained when addressing (3.2).

The approach we have just described yields two difficulties. First, it is important to choose a good value for ε. Second, problem (3.2) is unfortunately NP-hard. For patches of size m = 8 × 8 = 64 pixels, [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] recommend the value ε = m(1.15σ) 2 . Another effective heuristic proposed by [START_REF] Mairal | Non-local sparse models for image restoration[END_REF] assumes that the targeted residual y i -Dα i behaves as the Gaussian noise of (supposedly known) variance σ 2 , such that the quantity y i -Dα i 2 2 /σ 2 follows a χ-square distribution with m degrees of freedom. Then, the heuristic sets ε = σ 2 F -1 m (τ ) where F -1 m is the inverse cumulative distribution function of the χ 2 m distribution. Selecting the value τ = 0.9 leads to appropriate values of ε in practice for different patch sizes m. For dealing with (3.2), greedy algorithms such as orthogonal matching pursuit [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] are known to provide approximate solutions that are good enough for the denoising task [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. Such algorithms are presented in details in Section 5.1.

After having presented a simple "patch" denoising procedure given a fixed dictionary D, we study three questions that will explain the different steps needed for denoising a full image:

1. how do we reconstruct the full image from the local patch models? 2. which dictionary should we choose? 3. how does ℓ 1 compares with ℓ 0 for image denoising?

From patches to full image estimation. The denoising scheme proposed by [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] processes independently every patch y i from the noisy image y, before constructing the denoised image by averaging the patch estimates xi . More precisely, since the patches y i overlap, each pixel belongs to m different patches and thus admits m estimates. 3 Then, the denoised image x in R n is obtained as follows:

x = 1 m n i=1 R ⊤ i xi , (3.3)
where R i in R n×n is a binary matrix defined as in Section 2.6-that is, R ⊤ i xi is a vector of size n corresponding to an image with zeroes everywhere except for the patch indexed by i that contains xi . In other words, R ⊤ i is a linear operator that "positions" a patch of size m at the pixel i in a larger image with n > m pixels. 4 As a result, Eq. (3.3) is nothing else than an averaging procedure, which is the simplest way of aggregating estimators of the same quantity. We remark that this averaging strategy is related to the translationinvariant wavelet denoising technique of [START_REF] Coifman | Translation-invariant de-noising[END_REF], where a clean image is reconstructed by averaging estimates obtained by denoising several shifted versions of an input image. Even though aggregating estimators by straight averaging might look suboptimal, we are not aware of any other technique, in the context of local sparse linear models, leading to better results for reconstructing the final image from the estimation of overlapping patches.

Use of generic versus global versus image-adaptive dictionary.

The first baseline considered by [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] uses a fixed overcomplete discrete cosine transform (DCT) dictionary, which was presented earlier in Figure 2.4(a). We call such an approach "generic" since the dictionary is pre-defined. The second strategy, dubbed "global", consists of learning a dictionary on an external database of clean patches, simply following the methodology of Section 2. Finally, we also present the "adaptive" strategy of [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], which learns the dictionary from the pool of noisy patches:

min D∈C,A∈R p×n 1 n n i=1 ψ(α i ) s.t. y i -Dα i 2 2 ≤ ε, (3.4)
where C is defined in Section 2, and ψ is a sparsity-inducing penalty. [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] originally use the ℓ 0 -penalty in place of ψ, and optimize (3.4) with the K-SVD algorithm [START_REF] Aharon | An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. Later, additional experiments have shown that learning the dictionary with ℓ 1 could bring some benefits [START_REF] Mairal | Non-local sparse models for image restoration[END_REF]. In fact, we experimentally demonstrate in the sequel that ℓ 1 consistently provides better results for the denoising task, when used for learning the dictionary only.

4 Note that the original averaging scheme of [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] slightly differs from (3.3) since their estimate x is defined as

x = (1-λ)y+λ(1/m) n i=1 R ⊤ i xi,
where y is the noisy image. We have empirically found that the setting λ = 1 leads in fact to very good results and we have thus omitted the parameter λ.

Indeed, we observe that regardless of the way the dictionary is learned, ℓ 0 should always be used for the final reconstruction step in (3.2).

This conclusion is drawn from the following experiment. We consider 12 classical images presented in Figure 3.1, which were used in other benchmarks about image denoising [e.g., [START_REF] Mairal | Non-local sparse models for image restoration[END_REF]. We corrupt every image with white Gaussian noise of standard deviation σ in {5, 10, 15, 20, 25, 50, 100} for pixel values in the range [0; 255], and we measure the reconstruction quality obtained by the three approaches by using the PSNR criterion. 5 We consider dictionaries of size p = 256 elements, different patch sizes m = l × l, with l in {6, 8, 10, 12, 14, 16}, and we use the parameter τ = 0.9 suggested before for choosing the reconstruction threshold ε. For every noise level, the parameter l is selected such that it maximizes the average PSNR obtained on the last 5 images of the dataset.

We make a comparison between the six scenarios described below. In all cases, the dictionaries are learned with the SPAMS software, by performing 10 passes over the available training data.

1. DCT: we use an overcomplete DCT dictionary; the patches are reconstructed by using the ℓ 0 -regularization (3.2).

2. ℓ 0 -global/ℓ 0 : the dictionary is learned on 400 000 natural image patches extracted from the Kodak PhotoCD images6 . Denoting by x 1 , . . . , x n these patches, we learn D by minimizing min

D∈C,A∈R p×n 1 n n i=1 1 2 x i -Dα i 2 2 s.t. α i 0 ≤ s,
where s = 10 is known to approximate well clean natural image patches [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF].

3. ℓ 0 -adapt/ℓ 0 : the dictionary is learned by minimizing (3.4) with ψ = ℓ 0 ; we initialize the learning procedure with the global dictionary obtained in the scenario ℓ 0 -global/ℓ 0 . 4. ℓ 1 -global/ℓ 0 : the scenario is the same as ℓ 0 -global/ℓ 0 , except that the dictionary is learned with the ℓ 1 -penalty, following exactly the methodology of Section 2.4. The final reconstruction is obtained with ℓ 0 as in (3.2).

5. ℓ 1 -adapt/ℓ 0 : same as ℓ 1 -adapt/ℓ 0 , except that the dictionary is learned by minimizing (3.4) with ψ = ℓ 1 ; we use the dictionary obtained in the scenario ℓ 1 -global for initializing the learning procedure. The final reconstruction is still obtained with the ℓ 0penalty.

6. ℓ 1 -adapt/ℓ 1 : same as ℓ 1 -adapt/ℓ 0 , but the patches y c i in the final reconstruction are denoised by replacing the ℓ 0 -penalty by the ℓ 1 -norm in (3.2).

We report the mean PSNR obtained by the different scenarios in Table 3.1, along with the performance achieved by other state-of-the-art approaches of the literature. The conclusions are the following:

• the simple denoising scheme that we have presented performs slightly worse than more recent methods such as Dabov et al. [2007a], [START_REF] Mairal | Non-local sparse models for image restoration[END_REF], [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF] in terms of PSNR;

• adaptive dictionaries yield better results than global ones, as already observed by [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF];

• learning the dictionary with ℓ 1 consistently yields better results than ℓ 0 , even though the ℓ 0 -penalty is used for the final image reconstruction step in (3.2);

• ℓ 0 yields better results than ℓ 1 in the final image reconstruction step, regardless of the way the dictionary was learned.

Even though the conclusion that one should "first learn the dictionary with ℓ 1 , before using ℓ 0 for denoising the image" might seem counterintuitive, we believe that the reason of the good performance of ℓ 1 for dictionary learning might be a better stability of the sparsity patterns than the ones obtained with ℓ 0 . This may subsequently 3.1 for various levels of noise σ. GSM refers to the Gaussian scale mixture model of [START_REF] Portilla | Image denoising using scale mixtures of Gaussians in the wavelet domain[END_REF]; K-SVD refers to the original method of [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]; BM3D refers to the state-of-the-art denoising approach of Dabov et al. [2007a]; LSSC refers to [START_REF] Mairal | Non-local sparse models for image restoration[END_REF], Plow refers to [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF], EPLL to [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], CSR to [START_REF] Dong | Nonlocally centralized sparse representation for image restoration[END_REF], and BM3D-SAPCA to an extension of BM3D with better results [START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF][START_REF] Katkovnik | From local kernel to nonlocal multiple-model image denoising[END_REF]. For all these approaches, publicly available software is available on the corresponding authors' web pages. yield a better behavior in terms of optimization. Whereas the ℓ 1 -scheme guarantees us to obtain a stationary point of the dictionary learning formulation, the ℓ 0 -counterpart does not. A similar experimented was conducted by Mairal [2010, chapter 1.6], with similar conclusions.

Image inpainting

Dictionary learning is well adapted to the presence of missing datathat is, it is appropriate for inpainting [START_REF] Bertalmio | Image inpainting[END_REF] when the missing pixels form small holes that are smaller than the patch sizes. To deal with unobserved information, the dictionary learning formulation can be modified by introducing a binary mask M i for every patch indexed by i [Mairal et al., 2008c,e]. Formally, we define M i as a diagonal matrix in R m×m whose value on the j-th entry of the diagonal is 1 if the pixel y i [j] is observed and 0 otherwise. Then, the dictionary learning formulation becomes min

D∈C,A∈R p×n 1 n n i=1 1 2 M i (y i -Dα i ) 2 2 + λψ(α i ),
where ψ is a sparsity-inducing penalty, M i y i represents the observed pixels from the i-th patch of the image y and Dα i is the estimate of the full patch i. In practice, the binary mask does not drastically change the optimization procedure, and one can still use classical optimization techniques for dictionary learning, e.g., alternating between the optimization of D and A, as described in Section 5. When the image y is only corrupted by missing pixels and not by other noise source, it is also possible to enforce hard reconstruction constraints, leading to the formulation min

D∈C,A∈R p×n n i=1 ψ(α i ) s.t. M i y i = M i Dα i . (3.5)
Similarly to the denoising problem, once the dictionary is learned, we need to sparsely encode all overlapping patches, e.g., by approximately minimizing (3.5) with a greedy algorithm when ψ is the ℓ 0 -penalty. Then, the full image is obtained by averaging using (3.3).

Before showing any inpainting result, we shall comment on when the formulation described here is supposed to work. Our first remark is that it can only handle holes that are smaller than the patch size. Dealing with larger holes might be possible, but with a different formulation that would be able to synthesize new information in empty patches [see, e.g. [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF][START_REF] Peyré | Sparse modeling of textures[END_REF][START_REF] Roth | Fields of experts[END_REF]. Second, one assumes that the noise pattern is unstructured, such that no noise patterns are learned by the dictionary. The demosaicking task from the next section is a typical example with very structured missing patterns, and where an alternative strategy needs to be used.

Finally, we show inpainting results in Figure 3.2. For both images, the noise structure is relatively random when seen at the patch level, and the dictionary learning approach performs well at recovering the missing pixels. In particular, the brick texture in the image house is hardly visible for the human eye, but it is well recovered by the algorithm.

Image demosaicking

The problem of demosaicking consists of reconstructing a color image from the raw information provided by colored-filtered sensors from digital cameras [START_REF] Kimmel | Demosaicing: image reconstruction from color ccd samples[END_REF]. Most consumer-grade cameras use a grid of sensors; for every pixel, one sensor records the amount of light it receives for a short period of time. A color filter, red, green, or blue, is placed in front of the sensor such that only one color channel is measured for each pixel. Reconstructing the missing channels for the full image is a task called demosaicking. It is usually performed on-the-fly by digital cameras, but can also be achieved with a better quality by professional software on a desktop computer.

A typical pattern of filters, called the "Bayer" pattern, is displayed in Figure 3.3. Every row alternates between red and green, or green and blue filters. As a result, information captured by digital sensors in the green channel is twice more important than in the red or blue channels, even though the final post-processed color image observed by the user does not reflect this fact. An example of a mosaicked image with the Bayer pattern is presented in Figure 3.4(a). Reconstructing the color image without producing any visual artifact is difficult since it requires recovering frequencies in ever channel that are greater than the sampling rate. Specifically, the resolution in the red and blue channels needs to be doubled.

Since demosaicking consists of recovering missing information for every pixel, it can be cast as an inpainting problem, but the inpainting procedure described in the previous section cannot be directly applied because the noise pattern is strongly structured. An effective two-step strategy proposed by Mairal et al. [2008c] to overcome the shortcoming of the original inpainting formulation consists of the following steps:

G R G R G R G B G B G B G R G R G R G B G B G B G R G R G R G B G B G B
1. learn a dictionary D 0 on an external database of color patches; 2. use D 0 to obtain a first estimate of the color image, following the image reconstruction procedure of the previous section;

3. retrain the dictionary on the first image estimate, leading to a new dictionary D 1 adapted to the image content.

4. use D 1 to obtain the final estimate of the color image as in 2.

In addition, the patches should also be centered according to the recommendations of Section 2.1. The results achieved by this method can be were shown by Mairal et al. [2008c] to outperform the state of the art in terms of PSNR on the 24 images from the Kodak PhotoCD benchmark, even though some visual artifacts remain visible in areas that are particularly difficult to reconstruct, as shown in Figure 3.4(b).

The demosaicking approach based on dictionary learning has been improved later by [START_REF] Mairal | Non-local sparse models for image restoration[END_REF] by combining sparse estimation with image self-similarities [START_REF] Buades | Self-similarity driven color demosaicking[END_REF]. The resulting method yields a higher PSNR, but more importantly it seems to significantly reduce visual artifacts (see Figure 3.4(c)), and achieves state-of-the-art results in terms of PSNR [START_REF] Menon | Color image demosaicking: an overview[END_REF].

Image up-scaling

Dictionary learning has also gained important success for reconstructing high-resolution images from low-resolution ones, a problem called image up-scaling or digital zooming. 7 The technique originally proposed by Yang et al. [2010a] consists of learning a pair of dictionaries (D l , D h ) that can respectively reconstruct low-and high-resolution patches with the same sparse decomposition coefficients. Then, the relation between D l and D h is exploited for processing new low-resolution images and turn them into high-resolution ones. Similar to the "example-based" approach of [START_REF] Freeman | Example-based superresolution[END_REF], the method requires a database of pairs of training patches (x l i , x h i ) n i=1 , where x l i in R m l is a low-resolution version of the patch x h i in R m h , and m h > m l . The database is built from training images that are generic enough to represent well the variety of natural images. Then, the pairs of patches are used offline to train the dictionaries D l in R m l ×p and D h in R m h ×p , which are subsequently used for image up-scaling in a way that will be described in the sequel. The methodology is related to the global approach for image denoising from Section 3.1, where the dictionary is not adapted to the image at hand that needs to be processed, but is adapted instead to a generic database of images.

The image upscaling method of Yang et al. [2010a] has received a significant amount of attention because of the quality of its results, and has been improved recently in various ways [START_REF] Couzinie-Devy | Dictionary learning for deblurring and digital zoom[END_REF][START_REF] Zeyde | On single image scale-up using sparse representations[END_REF], Dong et al., 2011b, Wang et al., 2012, Yang et al., 2012a]. In this section, we briefly review some of the ideas that have appeared in this line of work, but refer to the corresponding papers for more details. We start with the original approach of Yang et al. [2010a], before moving to two natural variants.

The original approach of Yang et al. [2010a].

The dictionary learning formulation of Yang et al. [2010a] for image up-scaling consists of 7 The terminology of super-resolution is sometimes used but leads to some confusion. "Super-resolution" traditionally means synthesizing high-resolution images from a sequence of low-resolution ones [see the discussion in Elad, 2010, page 341]. Here, we are dealing with a single input image.

learning D l and D h with the following joint minimization problem: min

D l ∈C l D h ∈R m h ×p A∈R p×n 1 n n i=1 1 2m l x l i -D l α i 2 2 + 1 2m h x h i -D h α i 2 2 + λ α i 1 ,
(3.6) where A = [α 1 , . . . , α n ] is the matrix of sparse coefficients, and C l is the set of matrices in R m l ×p whose columns have less than unit ℓ 2 -norm. The formulation encourages the dictionaries to provide a common sparse representation α i for the pair of patches x l i and x h i . Note that, to simplify, we have omitted pre-processing steps for the patches x l i , which lead to substantially improved results in practice [see Yang et al., 2010a[START_REF] Zeyde | On single image scale-up using sparse representations[END_REF], and we have preferred to focus on the main principles of the method.

Once the dictionaries are trained, they can be used for turning new low-resolution images into high-resolution ones. Specifically, let us consider an image y l of size √ n l × √ n l , which is independent from the training database, and assume that the goal is to synthesize a highresolution image y h of larger size

√ n h × √ n h -say, for instance, √ n h = 2
√ n l when one wishes to double the resolution.

According to (3.6), the main underlying assumption is that highresolution and low-resolution variants of the same patch admit a common sparse decomposition onto D h and D l , respectively. Then, given a patch y l i in R m l centered at pixel i in the low-resolution image y l , the method of Yang et al. [2010a] computes a sparse vector β i in R p such that y l i ≈ D l β i , and the high-resolution estimate of the patch is simply D h β i . Such an approach is however not fully satisfactory regarding the formulation (3.6), where each sparse decomposition is obtained by using both the low-and high-resolution patches. In other words, β i should be ideally computed by minimizing min

β i ∈R p 1 2m l y l i -D l β i 2 2 + 1 2m h y h i -D h β i 2 2 + λ β i 1 , (3.7)
but unfortunately only y l i is available at test time since y h i is the unknown quantity to estimate. Such a discrepancy is addressed by Yang et al. [2010a] with a heuristic, where y h i in (3.7) is replaced by an auxiliary variable representing the current prediction of the highresolution patches given previously processed neighbors. Finally, reconstructing the full image from the local patch estimates can be done as usual by straight averaging, even though other possibilities yielding good results have also been proposed [see [START_REF] Zeyde | On single image scale-up using sparse representations[END_REF].

First variation with a regression formulation.

A strategy for overcoming the discrepancy between training and testing in the previous approach has been independently proposed by several authors [START_REF] Zeyde | On single image scale-up using sparse representations[END_REF][START_REF] Couzinie-Devy | Dictionary learning for deblurring and digital zoom[END_REF], Yang et al., 2012a]. In all these approaches, the sparse coefficients α i and β i are always computed by using low-resolution patches only. This motivates the following two-step training approach:

1. compute D l and A with a classical dictionary learning formulation, for instance: min

D l ∈C l ,A∈R p×n 1 n n i=1 1 2 x l i -D l α i 2 2 + λ α i 1 .

obtain D h by solving a multivariate regression problem:

min

D h ∈R m h ×p 1 n n i=1 1 2 x h i -D h α i 2 2 ,
where the α i 's are fixed after the first step. See also [START_REF] Zeyde | On single image scale-up using sparse representations[END_REF] for other variants.

Once the dictionaries are learned, a new low-resolution image y l is processed by decomposing its patches y l i onto D l , yielding sparse coefficients β i , and the estimate of the corresponding high resolution patches are again the quantities D h β i .

Second variation with bilevel optimization. The previous variant addresses one shortcoming of the original approach of Yang et al. [2010a], but unfortunately loses the ability of jointly learning D l and D h . The variant presented in this paragraph is motivated by that issue. It was independently proposed by Couzinie-Devy et al. [2011] and Yang et al.

[2012a] and allows learning D l and D h at the same time. As we shall see, it raises a challenging bilevel optimization problem.

Let us first introduce a notation describing the solution of the Lasso for any vector x and dictionary D:

α ⋆ (x, D) arg min α∈R p 1 2 x -Dα 2 2 + λ α 1 ,
where we assume the solution of the Lasso problem to be unique.8 Then, the joint dictionary learning formulation consists of minimizing min

D l ∈C l D h ∈R m h ×p 1 n n i=1 1 2 x h i -D h α ⋆ (x l i , D l ) 2 2 . (3.8)
The problem is challenging since it is non-convex, non-differentiable, and the dependency of the objective with respect to D l goes through the argmin of the Lasso formulation. Such a setting is sometimes referred to as "bilevel optimization" problems and is known to be difficult. Eq. (3.8) appears under the name of "coupled dictionary learning" in [Yang et al., 2012a], and is in fact a particular case of the more general "task-driven dictionary learning" formulation of Mairal et al.

[2012], which we will present in more details in Section 4.5 about visual recognition. It is possible to show that in the asymptotic regime where n grows to infinity, the cost function to optimize becomes differentiable and that its gradient can be estimated. Based upon this observation, a stochastic gradient descent algorithm can be used as an effective heuristic for finding an approximate solution [see [START_REF] Couzinie-Devy | Dictionary learning for deblurring and digital zoom[END_REF], Mairal et al., 2012]. We conclude this section with a few visual results, which we present in Figure 3.5.

Inverting nonlinear local transformations

Interestingly, the image up-scaling formulations that we have presented previously do not explicitly take into account the type of transformation between the patches x h i and x l i . In other words, they only assume that the pairs of patches x h i and x l i admit a common sparse representation onto two different dictionaries, but they do not make any a priori assumption on the relation between the patches. For this reason, it is appealing to try similar formulations for other problems than image upscaling.

Several extensions have been developed along this line of work. Couzinie-Devy et al. [2011] and Dong et al. [2011b] have for instance developed image deblurring techniques that are effective when the blur is local. Other less-standard applications also appear in the literature. Dong et al. [2011b] propose indeed to learn a mapping between face photographs and hand-drawn sketches. They use a database of images called CUFS [START_REF] Wang | Face photo-sketch synthesis and recognition[END_REF] of about 600 subjects. For each subject, a grayscale photograph is provided, along with a sketch drawn by an artist. The dataset can be used to create pairs of photographsketch patches and a nonlinear mapping is learned between the two patch spaces, providing convincing visual results.

Finally, we illustrate the ability of the previous formulations for learning nonlinear mappings with the problem of reconstructing grayscale images from binary ones, a task called "inverse halftoning". With the use of binary display and printing technologies in the 70's, converting a grayscale continuous-tone image into a binary one that looks perceptually similar to the original one ("halftoning") was of utmost importance. A classical algorithm was developed for that purpose by [START_REF] Floyd | An adaptive algorithm for spatial grey scale[END_REF]. The goal of "inverse halftoning" is to restore these binary images and estimate the original ones. Since building a large database of pairs of grayscale-halftoned image patches is easy, the methodology of the previous section can be used for learning a mapping between binary and grayscale patches. In Figure 3.6, we display a few visual results obtained with the approach of Mairal et al. [2012]. The dictionaries were trained on an external database build with the Floyd-Steinberg algorithm.

Video processing

Extensions of dictionary learning techniques for dealing with videos have been proposed by [START_REF] Protter | Image sequence denoising via sparse and redundant representations[END_REF]. Given a noisy video sequence, the most naive approach consists of processing each frame independently. However, significantly better results can be obtained by exploiting the temporal consistency that exists among consecutive frames. To do so, a few modifications to the classical denoising formulation are sufficient:

• several frames can be processed at the same time by considering three-dimensional patches that involve one temporal dimension.

For instance, a video patch can be of size m = e × e × T , where T is the number of frames in the patch, e.g., e = 10 pixels and T = 5 frames.

• After processing T frames starting at time t, we obtain a dictionary D t adapted to the set of three-dimensional patches at time t. When moving to the next block of T frames at time t + 1 (considering that the blocks overlap in time), we can learn a dictionary D t+1 adapted to the information available a time t + 1, using D t as an initialization to the learning algorithm.

In Figures 3.7 and 3.8, we present two video processing results from [Mairal et al., 2008e], where the original video extension of [START_REF] Protter | Image sequence denoising via sparse and redundant representations[END_REF] has been adapted to the inpainting and color video denoising tasks. The results of the "video processing" approach are significantly better than those obtained by processing independently each frame.

Face compression

We have shown so far that dictionary learning was appropriate for many restoration tasks, but an intuitive use of sparsity is for compression. Following this insight, [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF] have proposed a technique for encoding photographs of human faces. They achieve impressive compression results, where each image is represented with less than 1 000 bytes while suffering from minor loss only in visual quality. The algorithm requires a database of training face images that are processed as follows:

1. input images are geometrically aligned and downsampled to the same size. As a result, all the main face features are at the same position across images. This step is automatically performed by first detecting 13 face features and warping the images to achieve the desired alignment [see [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF];

2. the images are sliced into 15 × 15 non-overlapping patches;

3. for each patch position, a dictionary of size p = 512 elements is learned by using all patches available at the same position from the training images. The K-SVD algorithm is chosen by [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF] for that purpose. Once encoded, each patch can be decoded with a simple matrix vector multiplication, and the inverse geometrical warping is applied to obtain an approximation of the original image. Visual results are presented in Figure 3.9 for face images represented by 550 bytes only.

Unlike generic approaches such as JPEG and JPEG2000, the approach of [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF] can exploit the fact that each patch represent a specific information, e.g., eye, mouth, which can be learned by using adaptive dictionaries. Since the patches do not overlap, the reconstructed images suffer from blocking artifacts, which can be significantly reduced by using a deblocking post-processing algorithm.

Other patch modeling approaches

Dictionary learning was successful in many image processing tasks because of its ability to model well natural image patches. Other patchbased methods can be used for that purpose, and have also shown impressive results for image restoration. For the sake of completeness, we briefly present a few of them that have gained some success, but we refer the reader to the corresponding papers for more details.

Non-local means and non-parametric approaches. [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF] showed that self-similarities inherent to natural images could be used effectively in texture synthesis tasks. Following their insight, [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] have introduced the non-local means approach to image denoising, where the prominence of self-similarities is used as a prior on natural images. 9 Concretely, let us consider a noisy image written as a column vector y in R n , and denote by y[i] the i-th pixel and, as usual, by y i the patch of size m centered on this pixel for some appropriate size m. This approach exploits the simple but very effective idea that two pixels associated with similar patches y i and y j should have similar values y[i] and y[j]. Using y i as an explanatory variable for y[i] leads to the non-local means formulation, where the denoised pixel x[i] is obtained by a weighted average:

x[i] = n j=1 K h (y i -y j ) n l=1 K h (y i -y l ) y[j], (3.9)
and K h is a Gaussian kernel of bandwidth h. The estimator (3.9) is in fact classical in the non-parametric statistics literature, even though it was never applied to image denoising with natural image patches before; it is usually called Nadaraya-Watson estimator [see [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF][START_REF] Wasserman | All of nonparametric statistics[END_REF]. 10

9 This idea has in fact appeared in the literature in various guises and under different equivalent interpretations, e.g., kernel density estimation, mean-shift iterations [START_REF] Awate | Unsupervised, information-theoretic, adaptive image filtering for image restoration[END_REF], diffusion processes on graphs [START_REF] Szlam | Regularization on graphs with function-adapted diffusion processes[END_REF], and long-range random fields [START_REF] Li | Sparse long-range random field and its application to image denoising[END_REF].

10 Interestingly, the same class of non-parametric estimators have also been used independently by [START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF], where the explanatory variables are not full Later, several extensions have been proposed along this line of work. First attempts have focused on improving the computational speed [START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF], or on automatically adapting the kernel bandwidth h to the pixel of interest [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF]. Extensions of the non-local principle for other tasks than denoising have also been developed, e.g., for image upscaling [START_REF] Glasner | Super-resolution from a single image[END_REF], or demosaicking [START_REF] Buades | Self-similarity driven color demosaicking[END_REF]. Later, other estimators exploiting image self-similarities and based on nonparametric density estimation have been proposed and analyzed by [START_REF] Levin | Patch complexity, finite pixel correlations and optimal denoising[END_REF], [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF].

BM3D. Dabov et al. [2007a] proposed a patch-based procedure called BM3D for image denoising. The method exploits several ideas, including image self-similarities and sparse wavelet estimation, and provides outstanding results at a reasonable computational cost. Several years after it was developed, it remains considered as the state of the art. BM3D proceeds with the following pipeline for denoising an image y:

• block matching: like non-local means, BM3D exploits selfsimilarity; it processes all patches by first matching them with similar ones in the noisy image y, stacking them together into a 3D signal block;

• 3D filtering: each block is denoised by using hard or softthresholding with a 3D orthogonal DCT dictionary, following the wavelet thresholding tradition;

• patch averaging: as in Section 3.1, a denoised image x0 is obtained by averaging the estimates of the overlapping patches; 11

• refinement with Wiener filtering: the previous steps are refined by using the intermediate estimate x0 . First, block matching is applied to x0 instead of y in order to obtain a more reliable image patches, but neighbor pixels. It results in Nadaraya-Watson estimators that exploit local information, as opposed to the non-local one of [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF].

11 More precisely, a weighted average is performed. We omit this detail here for simplicity, and refer to [Dabov et al., 2007a] for more details.

matching of noisy patches from y; weights are obtained for the 3D-filtering by using a Wiener filter. After a new patch averaging step, a final estimate x1 is obtained.

A few additional heuristics are also used to further boost the denoising performance [see Dabov et al., 2007a]. Finally, the scheme has proven to be very efficient and gives substantially better results than regular non-local means. Later, it was improved by [START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF] by adding other components: (i) shape-adaptive patches that are non necessarily rectangular; (ii) 3D filtering with adaptive dictionaries obtained with PCA inside each block instead of using simple DCT.

Non-local sparse models. Since BM3D was successful in combining image self-similarities and wavelet/DCT denoising on blocks of similar patches, [START_REF] Mairal | Non-local sparse models for image restoration[END_REF] have proposed to exploit further this idea and combine the non-local means principle with dictionary learning.

One motivation is that non-local means approach has proven to be effective in general, but it fails in some cases. In the extreme, when a patch does not look like any other one in the image, it is impossible to exploit self-similarities to estimate the corresponding pixel value. To some extent, sparse image models based on dictionary learning can handle such situations when the patch at hand admits a sparse decomposition, but they suffer from another drawback: similar patches sometimes admit very different estimates due to the potential instability of sparsity patterns, which can result in practice in noticeable reconstruction artifacts. The idea of the non-local sparse model is that similar patches should admit similar sparse decompositions. By enforcing this principle, one hopes to obtain more stable decompositions and subsequently a better estimation.

Concretely, let us define for each patch y i the set S i :

S i j = 1, . . . , n s.t. y i -y j 2 2 ≤ ξ ,
where ξ is some threshold. The block S i is essentially constructed by following the block matching step of BM3D. What differs is the way S i is subsequently processed. Whereas BM3D performs some wavelet filtering, the non-local sparse model jointly decomposes the patches onto a previously learned dictionary D in R m×p .

To encourage the decompositions to be similar, the patches from the set S i are forced to share a joint sparsity pattern-that is, they should share a common set of nonzero coefficients. Fortunately, this can be achieved with a group-sparsity penalty. Denoting by A i [α ij ] j∈S i the matrix of coefficient in R p×|S i | corresponding to the group S i , the regularizer encourages the matrix A i to have a small number of nonzero rows, yielding a type of sparsity patterns that is illustrated in Figure 3.10.

Formally, the penalty can be written for a matrix A in R p×k as Ψ q (A)

p j=1 α j q 2 ,
where α j denotes the j-th row of A. When q = 1, we obtain the group-Lasso norm already presented in Section 1.3. When q = 0, the penalty simply counts the number of non-zero rows in A and thus plays the same role as the ℓ 0 -penalty in the context of group sparsity [START_REF] Tropp | Algorithms for simultaneous sparse approximation. part I: Greedy pursuit[END_REF]. Then, decomposing the patch y i with the penalty Ψ q on the set S i amounts to solving min

A i ∈R p×|S i | Ψ q (A i ) s.t. j∈S i y j -Dα ij 2 2 ≤ ε i .
As for the classical denoising approach based on dictionary learning presented in Section 3.1, the following ingredients are used to ultimately reconstruct the denoised image: (i) ε i is chosen according to some effective heuristics; (ii) the non-convex penalty Ψ 0 with greedy algorithms [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] is preferred to Ψ 1 for the final image reconstruction task once the dictionary has been learned; (iii) since the patches overlap, the final step of the algorithm averages the estimates of every pixel. A few additional heuristics are also used to improve the speed or quality of the algorithm [see [START_REF] Mairal | Non-local sparse models for image restoration[END_REF], including a refinement step as in BM3D.

The non-local sparse principle was also successfully applied by [START_REF] Mairal | Non-local sparse models for image restoration[END_REF] to image demosaicking (see Figure 3.4), and more generally to image interpolation [START_REF] Romano | Single image interpolation via adaptive non-local sparsity-based modeling[END_REF]. Later, other variants of non-local sparse models have been proposed, notably the centralized sparse representation technique of Dong et al. [2011a[START_REF] Dong | Nonlocally centralized sparse representation for image restoration[END_REF]. There, input patches are clustered and an orthogonal dictionary is learned for each cluster with PCA. Then, patches are sparsely encoded on the PCA dictionaries, and the decompositions coefficients inside each cluster are encouraged to be close to a per-cluster mean representation. As shown in Sections 3.1 and 3.3, the performance achieved by non-local sparse models is in general competitive with the state of the art for different tasks, e.g., denoising or demosaicking.

Gaussian mixture models. We have already mentioned in Section 2.7 that Gaussian mixture models have been successful for modeling natural image patches [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]Weiss, 2011, Yu et al., 2012]. Not surprisingly, such approaches also lead to impressive results for image denoising and other image reconstructions tasks such as image upscaling [START_REF] Yu | Solving inverse problems with piecewise linear estimators: from Gaussian mixture models to structured sparsity[END_REF]. Given a database of image patches, the model parameters are usually learned by using the EM-algorithm or one of its variants, and denoising can be achieved with classical estimators for probabilistic models, e.g., maximum a posteriori estimation.

Because of its ability to model well natural image patches and automatically discover interpretable visual patterns, one may wonder whether dictionary learning can be useful for visual recognition or not. Patches are indeed often used as lowest-level features in image analysis pipelines. For instance, this is the case of popular approaches based on bags of words [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF][START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], which exploit image descriptors that are nothing else than preprocessed natural image patches [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF][START_REF] Tola | DAISY: An efficient dense descriptor applied to wide-baseline stereo[END_REF]. This is also the case of convolutional neural networks [LeCun et al., 1998a], whose first layer performs local convolutions on raw pixel values.

It is thus tempting to believe that modeling patches is an important step for automatic image understanding and that dictionary learning can be a key component of recognition architectures. The different approaches that we review in this section seem to confirm this fact, even though recognition requires other properties that are not originally provided by sparse coding models, such as invariance or stability of image representations to local perturbations.

In Section 4.1, we show how dictionary learning has been used for image classification as an alternative to clustering techniques in bagof-words models, before reviewing in Section 4.2 numerous variants that have been proposed in the literature. Then, we present different approaches where sparse models are used as a classification tool for face recognition (Section 4.3), patch classification, and edge detection (Section 4.4). Finally, we draw some links between dictionary learning and neural networks in Section 4.5 with backpropagation rules and networks for approximating sparse models, before reviewing a few additional computer vision applications in Section 4.6. We show that dictionary learning has been successfully used in various guises for different recognition tasks.

A coding and pooling approach to image modeling

Following [START_REF] Boureau | Learning mid-level features for recognition[END_REF], we show in this section that interleaved steps of (nonlinear) feature coding and pooling are a very simple but common approach to image modeling in visual recognition tasks such as retrieval and categorization. It is, implicitly or explicitly, used in SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], bags of features [Csurka et al., 2004, Sivic and[START_REF] Sivic | Video Google: a text retrieval approach to object matching in videos[END_REF], HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], the pyramid match kernel [START_REF] Grauman | The pyramid match kernel: Discriminative classification with sets of image features[END_REF], spatial pyramids [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], soft quantization [START_REF] Van Gemert | Visual word ambiguity[END_REF], and most importantly for this presentation, sparse coding approaches to recognition [START_REF] Yang | Linear spatial pyramid matching using sparse coding for image classification[END_REF], Boureau et al., 2010[START_REF] Wang | Locality-constrained linear coding for image classification[END_REF], Yang et al., 2010b].

We also show that the coding/pooling approach can often be intuitively justified in terms of (approximate) feature matching. We first consider global pooling approaches that discard all spatial information for every feature, before moving to local pooling that keeps part of it.

Global pooling

Bags of features and sparse coding. This paragraph largely follows [START_REF] Boureau | Learning mid-level features for recognition[END_REF]. Let us consider an image with its associated local image features, represented by (and identified with) a finite set X = [x 1 , . . . , x n ] of n vectors in R m . For instance, each x i may represent a non-rotationally invariant SIFT descriptor computed at location i in the image, but it may as well represent another feature type, e.g., local binary pattern [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] or a color histogram.

Let us also assume that we are given a coding operator, that is a function α : R m → A mapping features in R m onto the corresponding codes in some space A, and a pooling operator β mapping finite subsets of A onto elements of R p . The simplest version of the coding/pooling approach is to model an image represented by a set of descriptors X in R m×n by the vector γ(X) = β[α(X)], where α(X) denotes the set of points [α(x i )] n i=1 in A n . We dub this approach global pooling because it summarizes the codes of all features in the image by a single vector in R p .

The best known example of global pooling is the bag-of-features model inherited from text processing and popularized in computer vision by [START_REF] Sivic | Video Google: a text retrieval approach to object matching in videos[END_REF] and [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF]: let us consider a quantizer ξ : R m → {1, . . . , p}, constructed (for example) using the algorithm K-means on the elements of X. We can use ξ to encode each feature x as the binary vector α(x) of dimension p with all zero entries except for entry number ξ(x), which is equal to one. Note that the code space A is {0, 1} p in this case.

We can now use mean pooling (also called average pooling) to summarize the feature codes α(x i ) of the whole image by their averagethat is, a single vector of dimension p:

1 γ(X) = β[α(X)], where β[α(X)] = 1 n n i=1 α(x i ).
The vector γ(X) is called the bag of features associated with the image. Two images respectively represented by the local descriptors X = [x 1 , . . . , x n ] in R m×n and Y = [y 1 , . . . , y n ′ ] in R m×n ′ can now be compared by computing the dot product of their bags of features, that is, their similarity is computed as s(X, Y) = γ(X) ⊤ γ(Y). This approach is commonly used in both image retrieval (where s is used for ranking) [START_REF] Sivic | Video Google: a text retrieval approach to object matching in videos[END_REF] and categorization (where s is used as a kernel) [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF]. In both cases, the similarity function can be computed efficiently due to the fact that bags of features are typically very sparse for large values of p.

When the quantizer ξ is obtained using the K-means algorithm, the centroids d j (j = 1, . . . , p) of the corresponding clusters can be seen as the elements ("atoms", or "words") of some visual dictionary, and the codes α(x) tell us what atom a given feature is associated with.

This approach readily generalizes to other types of feature coding. In particular, it is possible to replace the vector-quantized binary codes of the features x in X by their sparse code relative to some (given or learned) dictionary D [START_REF] Yang | Linear spatial pyramid matching using sparse coding for image classification[END_REF], Boureau et al., 2010], that is

α(x) ∈ arg min α ′ ∈R p 1 2 x -Dα ′ 2 2 + λ α ′ 1 ,
then use mean pooling as before.

Feature matching interpretation. Here, we follow [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF], and show that, in the case of bags of features constructed with a quantifier, the similarity function s(X, Y) can be justified intuitively by the fact that it measures the proportion of pairs of features in X and Y that match each other [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF]. Maximizing s(X, Y) can thus be thought of as some (approximate) feature matching process. Indeed, given a quantizer ξ, let us define the functions δ j : R m → {0, 1} (j = 1, . . . , p) and δ : R m × R m → {0, 1} by

δ j (x) =    1 when ξ(x) = j, 0 otherwise, and δ(x, y) =    1 when ξ(y) = ξ(x), 0 otherwise. (4.1)
We can now rewrite our similarity function as

s(X, Y) = p j=1 1 n n i=1 δ j (x i ) 1 n ′ n i ′ =1 δ j (y i ′ ) = 1 n n ′ n i=1 n ′ i ′ =1 δ(x i , y i ′ ),
which indeed measures the fraction of pairs of features in the two images that are deemed to match when they admit the same code.

When the quantizer ξ is obtained by running the K-means algorithm, δ(x, y) is nonzero when the features x and y both lie in the same cell of the Voronoi diagram of the p centroids of the corresponding clusters. Related approaches explicitly phrased in term of feature correspondences have been used in both image matching and retrieval using some nearest-neighbor (or NN) decision rule such as ε-search or K-NN. Concretely, one measures again the fraction of matching features as

s(X, Y) = 1 n n ′ n i=1 n ′ i ′ =1 δ(x i , y i ′ ),
but, this time the function δ is defined by

δ(x, y) =    1 when y -x 2 < ε, 0 otherwise, (4.2)
in the ε-search case, or, in the K-NN case, by

δ(x, y) =       
1 when y is one of the K nearest neighbors of x, 1 when x is one of the K nearest neighbors of y, 0 otherwise.

(4.3)
This approach is the basis of the classical image retrieval and image matching techniques of [START_REF] Schmid | Local grayvalue invariants for image retrieval[END_REF] and [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF]. As in the case of bags of features [START_REF] Sivic | Video Google: a text retrieval approach to object matching in videos[END_REF], the matching step if often complemented by some geometrical verification stage in practice. Similar schemes have also been used in image categorization [START_REF] Wallraven | Recognition with local features: the kernel recipe[END_REF].

It should be noted that the coding/pooling models also admits a feature-matching interpretation when images are encoded using sparse coding instead of a quantifier: we can write in this case the similarity s(X, Y) of two images, respectively represented by the sets of descriptors X and Y, as

s(X, Y) = p j=1 1 n n i=1 α j (x i )   1 n ′ n ′ i ′ =1 α j (y i ′ )   = 1 n n ′ n i=1 n ′ i ′ =1 δ(x i , y i ′ ),
where α(x) = [α 1 (x), . . . , α p (x)] and δ(x, y) = α(x) ⊤ α(y). Note that by analogy, δ(x, y) can be rewritten in a similar form to the previous cases (4.1), (4.3), or (4.2):

δ(x, y) =    α(x) ⊤ α(y) when ∃j, α j (x)α j (y) = 0, 0 otherwise.
The function δ encodes in this case the fact that two feature vectors are deemed to match when they both have at least one nonzero code element for some dictionary atom, the contribution of each matching feature pair being proportional to the dot product of their codes.

Max pooling. The construction of the vector γ(X) associated with an image can be seen as pooling the code vectors [α(x i )] n i=1 associated with the image features and summarizing them by a single vector γ(X) which is their average. As noted before, this process is thus often dubbed mean, or average pooling. A simple variant is sum pooling, where the vectors γ(X) are not normalized and simply sum up the vectors α(x i ).

It has proven useful in several classification tasks to replace mean pooling by max pooling [START_REF] Serre | Object recognition with features inspired by visual cortex[END_REF], where

γ j (X) = max i=1,...,n α j (x i ),
and γ(X) = [γ 1 (X), . . . , γ p (X)]. In the case of bags of features, the vectors α(x i ) are binary, and computing the similarity s(X, Y) amounts to counting the number of bins in the quantization associated with at least one feature in each image. Intuitively, this may be justified in an application such as Video Google [START_REF] Sivic | Video Google: a text retrieval approach to object matching in videos[END_REF], where the query image (typically a box drawn around an object of interest such as a tie or a plate) is normally much smaller than the database images (typically depicting an entire scene) so that features occurring often in the scene are not given too much importance.

Max pooling is also used with sparse coding. This is a priori problematic because of the sign ambiguity of dictionary learning discussed in Section 1.6. Indeed, as noted by [START_REF] Murray | Generalized max pooling[END_REF], max pooling should only be used (or perhaps more accurately, is only intuitively justified) when the individual vector entries encode a strength of association between a descriptor and a codeword, and are thus all positive. This is the case for bags of features but not for sparse coding, with max pooling resulting potentially in 2 p different pooled features. The construction of Section 1.6 can be used to remove this ambiguity, but does not provide an intuitive justification of max pooling in this setting. It may then make sense to use a variant of dictionary learning that enforces non-negativity constraints on the sparse codes, or to duplicate the entries of α(x) into negative and positive components, resulting in a vector of size 2p, before applying max pooling. Specifically, each entry α j (x) will be duplicated into two values max(α j (x), 0) and max(-α j (x), 0).

Finally, other nonlinear pooling technique have also been studied in the literature [see [START_REF] Koniusz | Comparison of mid-level feature coding approaches and pooling strategies in visual concept detection[END_REF]; for simplicity, we restrict our presentation to the average and max-pooling strategies, which are the most popular ones.

Local pooling

Spatial pooling: spatial pyramids and their cousins

The global coding/pooling approach yields "orderless" image models [START_REF] Koenderink | The structure of locally orderless images[END_REF], where all spatial information has been discarded. While this affords a great deal of robustness (including a total invariance to image transformations as long as they leave the local features invariant),2 it seems wasteful to discard all spatial information.

The spatial pyramid of [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] addresses this problem by overlaying a coarse pyramid with L levels over the image (the HOG model of [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] is based on a similar idea using a coarse grid instead of a pyramid). There are 2 2l cells at level l (l = 0, . . . , L -1) of the pyramid, and the features falling in each cell of the pyramid are binned separately. Let us define by S kl the indices of features in {1, . . . , n} falling into cell number k at level l; the spatial pyramid image descriptor γ(X) is obtained by concatenating the unnormalized histograms (sum pooling) γ kl (X) = i∈S kl α(x i ) associated with all cells at all levels of the pyramid to form a vector in R pd , where d = 1+. . .+d L-1 = (2 2L -1)/3, and d l = 2 2l for l = 0, . . . , L-1.

Note that the case L = 1 corresponds to the global pooling model, as used in its bag of features or sparse coding instances, except for the fact that it uses sum pooling instead of mean or max pooling, that is, the sum of the features is used instead of their mean or max value. As explained in the next section, the spatial pyramids associated with two images can be compared using their dot product or the pyramid matching kernel of [START_REF] Grauman | The pyramid match kernel: Discriminative classification with sets of image features[END_REF], and they can be interpreted in terms of feature matching in both cases.

Feature matching interpretation.

Given again two set of descriptors X and Y and a spatial pyramid structure leading to respective sets of indices S kl and S ′ kl for the cell k at level l, we can compute the similarity

s(X, Y ) = γ(X) ⊤ γ(Y) = L-1 l=0 d l k=1 i∈S kl i ′ ∈S ′ kl α(x i ) ⊤ α(y i ′ ).
This similarity function measures the number of matches between the two images, where matches are found at different spatial scales as pairs of features falling in the same cell.

As shown by [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], an alternative is provided by the pyramid match kernel, defined as follows. Let us first define the histogram intersection function as

I l (X, Y) = d l k=1 min[γ kl (X), γ kl (Y)],
where γ kl (Y) = i∈S ′ kl α(y i ), γ kl (X) = i∈S kl α(x i ), and where the min operator is applied componentwise. The j-th coordinate of I l (X, Y) measures the number of matches between X and Y corresponding to the atom d j of the dictionary, measured as the number of points from X and Y that have nonzero codes, and fall in the same cell. These points match in the code space because they have a nonzero element in the same spot, and they also match spatially, because they fall into the same cell.

By construction, the matches found at level l of the pyramid also include the matches found at level l+1. The number of new matches found at level l is thus I l (X, Y) -I l+1 (X, Y) for l = 0, . . . , L -1. This suggests the following pyramid match kernel [START_REF] Grauman | The pyramid match kernel: Discriminative classification with sets of image features[END_REF] to compare two images:

K(X, Y) = I L (X, Y) + L-1 l=0 1 2 L-l [I l (X, Y) -I l+1 (X, Y)] = 1 2 L I 0 (X, Y) + L l=0 1 2 L-l+1 I l (X, Y),
where the weight 1 2 L-l , which is inversely proportional to the cell width at level l is used to penalize matches found in larger cells because they involve increasingly dissimilar features. It is easily shown that K is a positive-definite kernel, which allows using the machinery of kernel methods and reproducing kernel Hilbert spaces [Shawe- [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF].

Like the histogram intersection function, the pyramid match kernel can be interpreted as measuring the number of matches between X and Y in both feature space and the spatial domain. Contrary to the previous cases discussed in this presentation, where all pairs of matching features were counted, a point of X may only match a unique point of Y in this case.

Feature space pooling. The pyramid match kernel was originally proposed by [START_REF] Grauman | The pyramid match kernel: Discriminative classification with sets of image features[END_REF] as a method for matching images (or equivalently counting approximate correspondences) in feature space, without retaining any spatial information. With spatial pyramids, [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] argued that it might be better suited to encoding spatial information in the two-dimensional image, using instead traditional vector quantization based on the K-means algorithm to handle matching in the high-dimensional feature space.

The two approaches can be thought of as performing a local form of pooling, over bins defined in feature space in [START_REF] Grauman | The pyramid match kernel: Discriminative classification with sets of image features[END_REF], or over cells defined in image space in [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF]. In practice, the latter method usually gives better results for image categorization, but it may be interesting to combine the two when the feature codes are obtained using sparse coding. This is precisely what [START_REF] Boureau | Ask the locals: Multi-way local pooling for image recognition[END_REF] have proposed to do [see also [START_REF] Gao | Local features are not lonely-Laplacian sparse coding for image classification[END_REF], 2013[START_REF] Wang | Locality-constrained linear coding for image classification[END_REF], Yang et al., 2012b[START_REF] Zhou | Image classification using supervector coding[END_REF][START_REF] Koniusz | Spatial coordinate coding to reduce histogram representations, dominant angle and colour pyramid match[END_REF], for related work that will be partly detailed in the next section]. In this approach, a dictionary is learned on the training data, and the corresponding sparse codes are then clustered using K-means. Let S klj now denote the index set of features that fall in cell number k of the image at level l of the pyramid, and in the Voronoi cell number j of the sparse code space, the unnormalized histogram γ kl (X) of spatial pyramid is simply replaced by γ klj (X) = i∈S klj α(x i ), resulting in a final descriptor γ(X) of size pdK where K is the number of centroids used by the K-means algorithm. Combined with max pooling, this method indeed gives very good results for image categorization on standard benchmarks such as Caltech 101.

The botany of sparse feature coding

To address the limitations of the original dictionary learning formulation of Olshausen andField [1996, 1997] for feature encoding, many different variants have been proposed in the literature. In this section, we go through a few of them that have gained a significant amount of attention and have been reported to improve upon the original codebook learning approach presented in the previous section.

Local coordinate coding. [START_REF] Yu | Nonlinear learning using local coordinate coding[END_REF] have proposed a significantly different view of the dictionary learning problem than the one presented so far in this monograph, where the learned dictionary elements are interpreted as anchor points on a nonlinear smooth manifold. 3 More precisely, these anchor points define a local coordinate system to represent data points, e.g., natural image patches or local descriptors. To make this interpretation relevant, [START_REF] Yu | Nonlinear learning using local coordinate coding[END_REF] have introduced the concept of locality for the dictionary elements, encouraging them to be close to data points. In practice, the resulting formulation, called "local coordinate coding", consists of optimizing the following cost function:

min D∈R m×p ,A∈R p×n 1 n n i=1   1 2 x i -Dα i 2 2 + λ p j=1 x i -d j 2 2 |α i [j]|   ,
(4.4) where, as usual, X = [x 1 , . . . , x n ] in R m×n is the database of training signals, A = [α 1 , . . . , α n ] is the matrix of sparse coefficients, and D = [d 1 , . . . , d p ] is the dictionary. In contrast to the classical dictionary learning formulation, the regularization function is a weighted ℓ 1 -norm, where the quadratic weights x id j 2 2 encourage the selection of dictionary elements that are close to the signal x i in Euclidean norm. Optimizing (4.4), meaning finding a stationary point since the problem is nonconvex, can be achieved with alternate minimization between D and A, which is a classical strategy for dictionary learning (see Section 5).

Note that the link between (4.4) and the manifold assumption is not obvious at first sight. In fact, [START_REF] Yu | Nonlinear learning using local coordinate coding[END_REF] show that (4.4) can be interpreted as a practical heuristic for a slightly different formulation that has precise theoretical guarantees, namely an approximation bound for modeling the supposedly existing manifold. Later, the locality principle was revisited by [START_REF] Yu | Improved local coordinate coding using local tangents[END_REF] with an improved approximation scheme, and by [START_REF] Wang | Locality-constrained linear coding for image classification[END_REF] with a practical feature encoding scheme.

The latter approach, dubbed "locality-constrained linear coding (LLC)", exploits the locality principle in a simple way to select a predefined number K < p of non-zero coefficients for the codes α i instead of using the ℓ 1 -regularization. Given a fixed dictionary D and a signal x i , the LLC scheme forces the support of α i to correspond

exists p vectors D(x) = [d1(x), . . . , dp(x)] in R m×p such that for all x ′ in R m , min α∈R p x ′ -x -D(x)α 2 ≤ C x ′ -x 2 2 .
to the K-nearest dictionary elements to x i . The value of the non-zero coefficients are then computed in analytical form without having to solve a sparse regression problem. Recent reviews and benchmarks have shown that such a coding scheme is competitive for image classification tasks [START_REF] Sánchez | Image classification with the fisher vector: Theory and practice[END_REF].

Laplacian sparse coding. We have seen in Section 3.8 with the nonlocal sparse models that a successful idea for image restoration is to encourage similar signals to share similar sparse decomposition patterns. For visual recognition, [START_REF] Gao | Local features are not lonely-Laplacian sparse coding for image classification[END_REF][START_REF] Gao | Laplacian sparse coding, hypergraph Laplacian sparse coding, and applications[END_REF] have shown that such a principle can be also useful. Given a training set X = [x 1 , . . . , x n ] of local descriptors, their formulation called "Laplacian sparse coding" consists of minimizing a regularized dictionary learning objective function min

D∈C,A∈R p×n 1 n n i=1 1 2 x i -Dα i 2 2 + λ α i 1 + µ i,j W ij α i -α j 2 2 ,
where W = [W ij ] in R n×n is a similarity matrix; W ij should be typically large when x i and x j are similar according to some appropriate metric, encouraging therefore α i to be close to α j in Euclidean norm. The terminology of "Laplacian" comes from the machine learning literature about semi-supervised learning where model variables sit on a graph [Belkin andNiyogi, 2003, 2004] and where such regularization functions are used. [START_REF] Gao | Local features are not lonely-Laplacian sparse coding for image classification[END_REF][START_REF] Gao | Laplacian sparse coding, hypergraph Laplacian sparse coding, and applications[END_REF] further discuss the choice of a good similarity matrix W when the signals x i are local descriptors, and experimentally show that the additional regularization yields better classification results than the traditional sparse coding approach of [START_REF] Yang | Linear spatial pyramid matching using sparse coding for image classification[END_REF].

Convolutional sparse coding. In Section 2.6, we have presented the convolutional sparse coding model applied to natural images without any concrete application. [START_REF] Zeiler | Adaptive deconvolutional networks for mid and high level feature learning[END_REF] have applied such a formulation for visual recognition with a multilayer scheme inspired from convolutional neural networks and from other hierarchical models [START_REF] Serre | Object recognition with features inspired by visual cortex[END_REF].

For each layer of the network, the formulation leverages the concept of "feature maps"-that is, a set of two-dimensional spatial maps that carry coefficients corresponding to a particular feature type at different locations. A feature map x has typically three dimensions, e.g.,

x is in R √ l× √
l×p for a square feature map of spatial size √ l × √ l and p different feature types, but it is sometimes more convenient to represent x as a vector of size pl. It is also convenient to consider three-dimensional spatial patches in the feature maps, e.g., a patch of size √ e× √ e×p that contains all information across feature maps from a spatial window of size √ e × √ e centered at a particular location.

Assuming that the parameters of the network have been already learned, one layer of the hierarchical model processes input feature maps produced by the previous layer (the first input feature map at the bottom of the hierarchy being the image itself), and produces an output feature map. Each layer performs successively two operations, which are illustrated in Figure 4.1:

1. convolutional sparse coding: the input feature map x is encoded by using the convolutional sparse coding formulation already presented in Section 2.6, which we recall now. Assuming that a dictionary D in R m×p has been previously learned, the feature map x-represented as a vector in R l here-is encoded as follows:

min

A∈R p×l 1 2 x - l k=1 R ⊤ k Dα k 2 2 + λ l k=1 α k 1 ,
where R k is the linear operator that extracts the patch centered at the k-th location from x and R ⊤ k positions a small patch at a location k in a larger feature map, using a similar notation as in (3.3). This operation produces sparse codes A = [α 1 , . . . , α l ] in R p×l , with p different coefficients for each position k. Then, the matrix A can be interpreted as an intermediate feature map of size √ l × √ l with p channels, as shown in Figure 4.1.

three-dimensional max pooling:

The spatial resolution and the number of channels of A is typically reduced by 2 with a max- j=1 are obtained after a three-dimensional max pooling step, reducing the spatial resolution and the number of maps from p to p ′ . pooling operation already presented in Section 4.1 with threedimensional pooling regions.

Separable sparse coding.

A variant of the convolutional sparse coding model has also been investigated by [START_REF] Rigamonti | Learning separable filters[END_REF], where the dictionary elements are spatially separable, or partially separable. More precisely, since a dictionary element d j in R m represents a twodimensional square patch, it is possible to reorganize its entries as a matrix Mat(d j ) in R √ m× √ m . Then, a dictionary element d j is said to be separable when there exist some vectors u j and v j in R √ m such that Mat(d j ) = u j v ⊤ j . In other words, separability simply means that the matrix Mat(d j ) is rank one.

The motivation for looking for separable dictionary elements is to accelerate the computation of the inner products d ⊤ j R k x for all overlapping patches R k x from an input image x in R l . This step indeed dominates the computational cost of reconstruction algorithms for the convolutional sparse coding model, and thus it is critical to make it efficient. The key observation is that computing all products d ⊤ j R k x is equivalent to performing a convolution of the two-dimensional filter Mat(d j ) on the input image x, which naively requires O(ml) operations. 4 In the separable case, we remark that we have

d ⊤ j R k x = trace(Mat(d j ) ⊤ Mat(R k x)) = trace(v j u ⊤ j Mat(R k x)) = u ⊤ j Mat(R k x)v j .
It is then easy to see that the inner products can be obtained by convolving twice the input image x: first with the filter u j along the vertical direction, then with v ⊤ j along the horizontal direction. As a result, the complexity drops to O( √ ml) operations. [START_REF] Rigamonti | Learning separable filters[END_REF] further study these filters for threedimensional data, where the gain in complexity is even more important than in the two-dimensional case. They also show how to learn partially separable filters by using the following regularization function for the dictionary:

ϕ(D) p j=1 Mat(d j ) * ,
where . * is the trace norm presented in Section 1.3 As a result, the penalty ϕ encourages the matrices Mat(d j ) to be low-rank instead of exactly rank one, meaning that they can be expressed by a sum of a few rank-one matrices.

Another variant of "separable dictionary learning" has also been proposed at the same time by [START_REF] Hawe | Separable dictionary learning[END_REF] in a non-convolutional setting, where the dictionary is assumed to be the Kronecker product of two smaller dictionaries:

D = B ⊗ C, where B is in R m b ×p b and C is in R mc×pc .
In other words, the matrix D is in R m b mc×p b pc and is a block matrix with m b × p b blocks, such that the i, j-th block is B[i, j]C. Similar to [START_REF] Rigamonti | Learning separable filters[END_REF], the motivation of [START_REF] Hawe | Separable dictionary learning[END_REF] is to make the matrix-vector multiplications Dα or D ⊤ x more efficient by exploiting the properties of Kronecker products [see [START_REF] Golub | Matrix computations[END_REF]. However, it is worth noticing that the concepts of "separability" used by [START_REF] Rigamonti | Learning separable filters[END_REF] and [START_REF] Hawe | Separable dictionary learning[END_REF] are not equivalent to each other.

Multipath sparse coding and hierarchical matching pursuit. Finally, another model called hierarchical matching pursuit and introduced by [START_REF] Bo | Hierarchical matching pursuit for image classification: Architecture and fast algorithms[END_REF][START_REF] Bo | Multipath sparse coding using hierarchical matching pursuit[END_REF] has recently obtained state-of-the-art results in various recognition tasks. Similar to the convolutional approach of [START_REF] Zeiler | Adaptive deconvolutional networks for mid and high level feature learning[END_REF], the scheme of [START_REF] Bo | Hierarchical matching pursuit for image classification: Architecture and fast algorithms[END_REF] produces for each image a sequence of feature maps organized in a multilayer fashion, but the coding scheme is somewhat simpler than in the convolutional sparse coding model. In a nutshell, the main differences are the following:

• sparse coding of independent patches: instead of using the convolutional model, all patches from an input feature map are coded independently;

• orthogonal matching pursuit: instead of using the ℓ 1regularization, the coding algorithm is a greedy method, which is presented in Section 5.1;

• contrast normalization: after each pooling step, patches are contrast-normalized, following the procedure described in Section 2.1.

Later, the multilayer scheme was improved by [START_REF] Bo | Multipath sparse coding using hierarchical matching pursuit[END_REF] by combining several networks that have a different number of layers and thus that have different invariant properties. A simpler version of this scheme with a single layer has also been shown to be effective for replacing lowlevel features such as histograms of gradients [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] for object detection tasks [START_REF] Ren | Histograms of sparse codes for object detection[END_REF].

Face recognition

One important success of sparse estimation for classification tasks is face recognition. [START_REF] Wright | Robust face recognition via sparse representation[END_REF] have indeed proposed a state-ofthe-art classifier for this task; as we will see, however, even though it was first evaluated on face datasets, the classifier is generic and can be applied to other modalities. [START_REF] Wright | Robust face recognition via sparse representation[END_REF] have proposed to sparsely decompose x onto the matrix D, and then measure which one of the sub-matrix D j is the most "used" in the decomposition. The precise procedure is given in Algorithm 1. Interestingly, this approach called "sparse-representation based classifier (SRC)" is related to other non-parametric machine learning techniques such as nearest neighbor classifiers that look for the most similar training samples to x before choosing the class label with a voting scheme [see, e.g., [START_REF] Hastie | The elements of statistical learning[END_REF], or methods looking for the nearest subspace [START_REF] Naseem | Linear regression for face recognition[END_REF], which project the data x onto each span of the submatrices D j , before selecting the subspace that best reconstructs the input data.

For the problem of face recognition, each matrix D j contains face samples of the same subject, and the method of [START_REF] Wright | Robust face recognition via sparse representation[END_REF] has gained a lot of popularity. However, this approach may suffer from several issues, which have been addressed by various extensions:

• occlusion: face images in realistic environments often contain unwanted sources of variations. Subjects may indeed wear clothes, scarves, glasses, which occlude part of the face, and which are not necessarily present in the training data. To improve the robustness of the formulation, [START_REF] Wright | Robust face recognition via sparse representation[END_REF] have introduced an auxiliary variable e in R p whose purpose is to model occluded parts. Then, Eq. (4.5) is replaced by the following linear program:

(α ⋆ , e ⋆ ) ∈ arg min α∈R p ,e∈R m [ α 1 + e 1 s.t. x = Dα + e] ,
Algorithm 1 Sparse-representation based Classifier (SRC).

Require: training instances

D = [D 1 , . . . , D k ] in R m×p , each D j rep- resents signals from class j; One test signal x in R m . 1: normalization: make each column of X of unit ℓ 2 -norm; 2: sparse decomposition: α ⋆ ∈ arg min α∈R p [ α 1 s.t. x = Dα] ;
(4.5)

3: classification:

 ∈ arg min j∈{1,...,k} x -D j δ j (α ⋆ ) 2 2 , (4.6)
where δ j selects the entries of α ⋆ corresponding to the class j; in other words, Dα ⋆ = k j=1 D j δ j (α ⋆ ); 4: return the estimated class  for the signal x.

where the ℓ 1 -norm encourages the variable e to be sparse. Thus, one hopes the variable e to have non-zero coefficients corresponding to occluded areas that are difficult to reconstruct with the training data D. The selection rule (4.6) needs also to be modified, and becomes

 ∈ arg min j∈{1,...,k} x -D j δ j (α ⋆ ) -e ⋆ 2
2 .

• image misalignment: the main assumption of the face recognition system of [START_REF] Wright | Robust face recognition via sparse representation[END_REF] is that a new test image can be sparsely decomposed as a linear combination of training images. Of course, such an assumption does not hold when the training images are not perfectly aligned-that is, when the face features for a subject are not localized at the same positions across training samples. To deal with that issue, Wagner et al.

[2012] have proposed dedicated solutions, which jointly optimize the variables α and e, along with a deformation variable whose purpose is to automatically "realign images".

Among other variants, it is worth mentioning some work using dictio-nary learning to build the matrices D j [see,e.g., [START_REF] Zhang | Discriminative K-SVD for dictionary learning in face recognition[END_REF], and also the use of random projection features [START_REF] Wright | Robust face recognition via sparse representation[END_REF].

Patch classification and edge detection

The capability of dictionary learning for modeling particular types of signals has been exploited for classification tasks in several ways. Specifically, there are two main successful strategies; the first one learns one dictionary per class of signal and compare reconstruction errors obtained by the dictionaries in the same vein as the SRC classifier of the previous section does. The other direction uses sparse codes produced by a single dictionary as a high-level representation of signals for a subsequent classification procedure. We successively present these two lines of research, along with some successful applications.

Finding the dictionary that best represents a class of signals. Let us consider a training set of signals X = [X 1 , . . . , X k ] in R m×n where the columns of each sub-matrix X j in R m×n j correspond to signals of the j-th class. Even though dictionary learning is an unsupervised learning technique, it can be used to model each class independently by learning a dictionary D j "adapted" to the data X j . For instance, it is possible to define D j as the output of a learning algorithm for optimizing min

D j ∈C,A j ∈R p×n j 1 2 X j -D j A j 2 F + λ n j i=1 α i 1 ,
where

A j = [α 1 , . . . , α n j ]. Then, a new test signal x in R m is classified according to the rule  = arg min j∈{1,...,k} min α∈R p 1 2 x -D j α 2 2 + λ α 1 . (4.7)
This simple rule was proposed by [START_REF] Ramirez | Classification and clustering via dictionary learning with structured incoherence and shared features[END_REF], following earlier work of Mairal et al. [2008a] that focused on the ℓ 0 -penalty and discriminative formulations. Even though the dictionaries are learned in an unsupervised manner, (4.7) performs surprisingly well for data that is well modeled by dictionary learning. For instance, the approach of [START_REF] Ramirez | Classification and clustering via dictionary learning with structured incoherence and shared features[END_REF] essentially builds upon (4.7) while also encouraging the dictionaries D j to be mutually incoherent, and obtains competitive results for digit recognition, namely about 1.2% error rate on the MNIST dataset. In contrast, a nonlinear support vector machine (SVM) with a Gaussian kernel achieves 1.4% and a K-nearest neighbor approach 3.9% [LeCun et al., 1998a]. Of course, better results have been obtained by other techniques on the MNIST dataset, but they are in general not based on generic classifiers. More precisely, state-of-the-art algorithms typically encode invariance to image deformations [Ranzato et al., 2007, Bruna and[START_REF] Bruna | Invariant scattering convolution networks[END_REF].

Since dictionary learning is supposed to be well adapted to natural image patches, the classification rule (4.7) can naturally be applied to model the local appearance of particular object classes. For instance, we show in Figure 4.2 some results obtained with a formulation that is closely related to (4.7) for a weakly supervised learning task [see Mairal et al., 2008a, for more details]. Because the local appearance of bicycles is highly discriminative, the method performs well at detecting patches that overlap with bicycles in the test images. [2008a]. Two dictionaries are learned beforehand, one from images containing a bicycle, and one from images corresponding to background only. For each of the four test images displayed here, we also display a confidence map for the presence of a bicycle. For each pixel, a confidence value is computed for a patch centered at the pixel, by comparing the residual of the reconstruction errors obtained with the two dictionaries; yellow bright pixel values represent high confidence, whereas dark red pixel values represent low confidence.

Another successful application of dictionary learning is edge detection. A first attempt by Mairal et al. [2008d] is based on the principles described in this paragraph, and consists of learning two dictionaries, one on patches centered on edges, and one centered on background. The method was trained on the manually annotated Berkeley segmentation dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] and was shown to perform as well as the dedicated approach called "Pb" [START_REF] Martin | Learning to detect natural image boundaries using local brightness, color, and texture cues[END_REF]. Later, stateof-the-art results have been obtained by [START_REF] Xiaofeng | Discriminatively trained sparse code gradients for contour detection[END_REF] by essentially combining two main ideas: (i) learning a single dictionary and exploiting the sparse codes for classification, a paradigm that will be the focus of the next paragraph, (ii) a dedicated way of performing multiscale feature pooling, giving more stability and some invariance to the sparse representations. We refer to [START_REF] Xiaofeng | Discriminatively trained sparse code gradients for contour detection[END_REF] for all details and we present some results obtained with their software package in Figure 4.3.5 In the next paragraph, we now move to the second paradigm, namely exploiting the sparse codes as a nonlinear transformation of the input data for a subsequent linear classifier. Using sparse codes for classification. Given a signal x in R m , the sparse coding principle produces a sparse vector α ⋆ such that x ≈ Dα ⋆ for some dictionary D in R m×p . The process transforming x into α ⋆ is highly nonlinear, and the entries of the vector α ⋆ can be interpreted as the "contributions" of the corresponding columns of D to the reconstruction of x.

When the dictionary atoms represent discriminative features for a classification task, it becomes appealing to use the vectors α ⋆ as new inputs to a classifier, e.g., linear SVM or logistic regression. The strategy is thus significantly different than the one presented in the previous paragraph. The proof of concept was first demonstrated by [START_REF] Huang | Sparse representation for signal classification[END_REF] with a pre-defined fixed dictionary and a cost function inspired from linear discriminant analysis [see [START_REF] Hastie | The elements of statistical learning[END_REF]. Later, dictionary learning was successfully used by [START_REF] Raina | Self-taught learning: transfer learning from unlabeled data[END_REF] for unsupervised feature learning as an effective way of exploiting unlabeled data for various classification tasks with image and text modalities.

Then, joint cost functions that involve both a discriminative term and a classical dictionary learning formulation have been proposed several times in the literature [START_REF] Rodriguez | Sparse representations for image classification: Learning discriminative and reconstructive non-parametric dictionaries[END_REF], Mairal et al., 2008b[START_REF] Pham | Joint learning and dictionary construction for pattern recognition[END_REF][START_REF] Yang | Fisher discrimination dictionary learning for sparse representation[END_REF]. Let us consider a training set (x i , y i ) n i=1 , where the vectors x i in R m are input signals and the scalars y i are associated labels. The simplest cost functions that combine sparse coding and classification have the following form min

D∈C,W∈R p×k A∈R p×n 1 n n i=1 1 2 x i -Dα i 2 2 +λ α i 1 +γL(y i , W, α i ) + µ 2 W 2 F ,
(4.8) where the loss function L measures how far a prediction based on α i is from the correct label y i . The matrix W represents model parameters for the classifier and the quadratic term (µ/2) W 2 F is a regularizer. Typically, L is a convex function and the prediction model is linear. For example, when the labels y i are in {-1, +1} for a binary classification problem, we may want to learn a linear decision rule parameterized by a vector w (meaning the matrix W has a single column), and use one of the following loss functions

L(y, w, α)          log 1 + e -yw ⊤ α (logistic loss), max 0, 1 -yw ⊤ α (hinge loss), 1 2 y -w ⊤ α 2 (square loss).
(4.9)

For simplicity, we omit the fact that the loss function may involve an intercept-that is a constant term b in R such that the linear model is y ≈ sign(w ⊤ α+b) instead of simply y ≈ sign(w ⊤ α). Depending on the chosen loss function L, the formulation (4.8) can combine dictionary learning with a linear SVM, or a logistic regression classifier. We present the three loss functions (4.9) in Figure 4.4. The logistic and hinge loss are very similar; both of them are asymptotically linear and encourage the signs of the products w ⊤ α i to be the same as y i . In all cases, it is possible to iteratively decrease the value of the cost function (4.8) and obtain a stationary point by alternate minimization between the variables D, A, W. As for the classical dictionary learning formulation, it is in fact impossible to find the global optimum; the problem is indeed convex with respect to each variable when keeping the other variables fixed, but it is not jointly convex. The motivation of the formulation (4.8) is to encourage the dictionary to provide sparse decomposition vectors α that are good for discrimination. It does so by encouraging them to have a small classification error on the training set measured by a convex loss function from machine learning. There are however two possible caveats with such an approach, which we discuss now.

The sparse decomposition patterns may be unstable. The first caveat is the lack of stability of the estimator provided by sparse decomposition algorithms. For instance, the solution of the Lasso can significantly change after a small perturbation of the input data-that is, the vector

α ⋆ (x, D) arg min α∈R p 1 2 x -Dα 2 2 + λ α 1 , (4.10)
assuming the optimization problem has a unique solution, can drastically change for small variations of x, or even small variations of λ in pathological cases [START_REF] Meinshausen | Stability selection[END_REF]Bühlmann, 2010, Mairal andYu, 2012]. For classification tasks, it is unfortunately often desirable to have stable representations, and thus improving this aspect may be important [START_REF] Bruna | Invariant scattering convolution networks[END_REF]. Several strategies have been proposed to deal with such an issue:

• promoting incoherent dictionaries: when D ⊤ D is close to the identity matrix, the dictionary is called "incoherent". When it is the case, sparse reconstruction algorithms may benefit from support recovery guarantees (see Section 1.5). Following this insight, [START_REF] Ramirez | Sparse modeling with universal priors and learned incoherent dictionaries[END_REF] have shown that promoting incoherent dictionaries may lead to better generalization properties. Later, this idea was also found useful in other classification schemes [START_REF] Bo | Multipath sparse coding using hierarchical matching pursuit[END_REF].

• pooling features: when a set of signals that are close to each other are available, e.g., a set patches that significantly overlap with each other, feature pooling on the set of corresponding sparse codes can lead to more stable, but less sparse representations. For instance, [START_REF] Ren | Histograms of sparse codes for object detection[END_REF], [START_REF] Bo | Multipath sparse coding using hierarchical matching pursuit[END_REF] adopts such a strategy in an effective image representation for object detection.

• replacing ℓ 1 by the elastic-net: when replacing the ℓ 1 -norm in (4.10) by the Elastic-net penalty α 1 + (γ/2) α 2 2 , it is possible to show that α ⋆ is Lipschitz continuous [Mairal et al., 2012], with a Lipschitz constant proportional to (1/γ) 2 . Here, the parameter γ controls a trade-off between sparsity and stability. As a side effect, it also leads to a unique solution of the sparse reconstruction problem due to the strong convexity of the penalty.

The label y is not available at test time. Interestingly, the formulation (4.8) is related to the super-resolution approach described in Section 3.4. even though the relation is only clear from hindsight. Indeed, the goal of classification is to map an input signal x to an output label y. In contrast, super-resolution can be seen as a multivariate regression problem, whose goal is to learn a mapping between x and a high-resolution signal-say, a signal y. It is then not surprising that the classification formulation (4.8) is very similar to the regression one (3.6). In both cases, learning the mapping involves a dictionary D for representing the input signals x and a model that maps sparse codes α to the variable to predict. In (4.8), the parameters W play for instance the same role as the parameters D h in (3.6).

In Section 3.4, we have discussed some limitations of the superresolution formulation (3.6), which are in fact also relevant for the classification task. Given some fixed D and W, we remark that the sparse codes α i for the training data may be obtained as follows

α i ∈ arg min α∈R p 1 2 x i -Dα 2 2 +λ α 1 +γL(y i , W, α). (4.11)
Unfortunately, given a new test signal x, the label y is not available and the formulation (4.10) has to be used instead to obtain the corresponding sparse code α. Therefore, there is a discrepancy in the way the signals are encoded during the training and test phases. In Section 3.4, we have addressed that issue by presenting an extension involving a bilevel optimization problem. It is also possible to follow the same strategy for the classification task, but we develop this method in the next section since it draws a strong connection between dictionary learning and neural networks.

Connections with neural networks

There are obvious connections between dictionary learning and neural networks. Some of them have already appeared in this monograph with multilayer schemes involving feature maps, which are typical architectures of convolutional neural networks and related work [LeCun et al., 1998a[START_REF] Serre | Object recognition with features inspired by visual cortex[END_REF]. In this section, we present two other important connections. First, we introduce the concept of "backpropagation" for dictionary learning, which is directly borrowed from neural networks [LeCun et al., 1998b], and which addresses an issue that was left opened at the end of the previous section. Then, we focus on particular networks whose purpose is to approximate efficiently the solution of sparse coding problems [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF].

Backpropagation rules for dictionary learning. In Figure 4.5, we present two paradigms for dictionary learning coupled with a prediction task, e.g., regression or classification. In the first one, the dictionary is obtained without supervision and the prediction model is learned in a subsequent step. In the latter paradigm, all parameters, including the dictionary, are learned jointly for the prediction task. This requires a similar concept as "backpropagation" [LeCun et al., 1998b], which consists of using the chain rule in neural networks.

Formally, let us consider a training set of signals X = [x 1 , . . . , x n ] in R m×n associated to measurements Y = [y 1 , . . . , y n ] in R k×n representing the variables to predict at test time. Such a setting is quite general on purpose: the vectors y i may represent for instance the highresolution patches from Section 3.4 or the half-toned patches from Section 3.5, but they may also represent a label for the classification tasks of Section 4.4. In this last case, k = 1 and the y i 's are simply scalars.

In the formulation of interest, the signals x i are encoded with D: (4.12) where the use of the elastic-net penalty instead of ℓ 1 will be made clear in the sequel. Then, we consider a similar prediction model as in the previous section; following the same methodology, we introduce a loss function L, and the corresponding empirical risk minimization problem min

α ⋆ (x, D) = arg min α∈R p 1 2 x -Dα 2 2 + λ α 1 + µ 2 α 2 2 ,
D∈C,W∈R p×k 1 n n i=1 L (y i , W, α ⋆ (x i , D)) + γ 2 W 2 F , (4.13)
Specifically, any loss from the previous section is considered to be appropriate here. The advantage of (4.12) over (4.11) is that the encoding scheme remains the same between training and testing, and yet D and W can be learned jointly for the final prediction task in (4.13). Moreover (4.13) extends the super-resolution formulation of (3.8) to more general problems such as classification.

Even though (4.13) seems appealing, it is unfortunately particularly difficult to solve. In addition to being nonconvex, the relation between D and the objective function goes through the argmin of a nonsmooth optimization problem. As a consequence, the cost function is also nonsmooth with respect to D and it is not clear how to obtain a descent direction.

In an asymptotic regime where enough training data is available, Mairal et al. [2012] show that the cost function becomes differentiable, and they derive a stochastic gradient descent algorithm from this theoretical result. They call this approach "task-driven dictionary learning". More precisely, consider the expected cost

f (D, W) E (y,x) [L (y, W, α ⋆ (x, D))] ,
where the expectation is taken to the unknown probability distribution of the data pairs (y, x). Under mild assumptions on the data distribution and when µ > 0, it is possible to show that f is differentiable and that its gradient admits a closed form:

   ∇ W f (D, W) = E (y,x) ∇L(y, W, α ⋆ ) ⊤ , ∇ D f (D, W) = E (y,x) -Dβ ⋆ α ⋆⊤ + (x -Dα ⋆ )β ⋆⊤ , (4.14)
where α ⋆ is short for α ⋆ (x, D),

β ⋆ [Γ c ] = 0, and β ⋆ [Γ] D ⊤ Γ D Γ + µI -1 ∇ α[Γ] L(y, W, α ⋆ ),
where Γ is the set of nonzero entries of α ⋆ . Here, the role of parameter µ is to ensure that the matrix to be inverted is well conditioned. Since the gradient has the form of an expectation, it is natural to use a stochastic gradient descent algorithm to optimize the cost function [START_REF] Bottou | The tradeoffs of large scale learning[END_REF]. Even though the nonconvexity makes it impossible to find the global optimum, it was shown by Mairal et al. [2012] that following heuristics from the neural network literature [LeCun et al., 1998b], this strategy yields sufficiently good results for many prediction tasks.

The concept of learning D to be good for prediction is similar to the backpropagation principle in neural networks where all the parameters of the network are learned for a prediction task using the chain rule [LeCun et al., 1998b]. To the best of our knowledge, [START_REF] Bradley | Differential sparse coding[END_REF] were the first to use such a terminology for dictionary learning; they proposed indeed a a different supervised dictionary learning formulation with smoothed approximations of sparsity-inducing penalties. In a heuristic fashion, similar rules as (4.14) have also been derived by [START_REF] Boureau | Learning mid-level features for recognition[END_REF] and Yang et al. [2010b] in the context of classification and by Yang et al. [2012a] for superresolution (see Section 3.4).

Fast approximations of sparse coding. Another link between neural networks and dictionary learning has been established by Kavukcuoglu et al. [2010b] and [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF], who have trained neural networks to approximate sparse codes. Given some fixed dictionary D and a training set X, their idea is to learn a simple nonlinear function g(W, x) that approximates α ⋆ (x, D) for any signal x of interest. By using similar notation as in the previous paragraphs, the general formulation proposed by [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF] is the following

min W∈W 1 n n i=1 α ⋆ (x i , D) -g(W, x i ) 2 2 ,
where W is a set of parameters for the prediction function g optimized with stochastic gradient descent and W is an optimization domain. The basic prediction function g introduced by Kavukcuoglu et al. [2010b] is the composition of a linear transformation and a pointwise nonlinearity, which can be interpreted as a one-layer neural network: An input signal x is sparsely encoded by using a dictionary D, producing sparse codes α. Then, a second model with parameters W makes a prediction y based on α. In (a), the dictionary D is learned unsupervised in a first layer-that is, without exploiting output information y. Then, the parameters W of the prediction model are learned afterwards in the second layer. In (b), we illustrate the concept of supervised dictionary learning, where D and W are jointly learned for the prediction task; this requires "backpropagating" information, represented by the red arrow, from the top layer to the bottom one.

g(W, x)[j] γ j tanh(w ⊤ j x + b j ), ( 4
where the scalars γ j and b j are some weights that can be also optimized. The motivation for such an approach is to reduce the cost of computing sparse codes at test time. Whereas obtaining α ⋆ (x, D) requires solving an optimization problem that can be costly (see Section 5), the complexity of computing g(W, w) is the same as a matrix-vector multiplication, making it possible to develop real-time applications.

However, there are two major shortcomings to the original approach of Kavukcuoglu et al. [2010b], as noted by [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF]. First, the choice of the non-linearity tanh does not yield exactly sparse vectors, but it can be replaced by more appropriate shrinkage functions. Second, the scheme (4.15) encodes each entry of the output vector g(W, x) independently and cannot possibly take into account the fact that there exists correlation among the dictionary elements. To address this issue, [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF] have proposed an iterative approach with a higher complexity than (4.15), but with a better approximation (see their paper for more details).

Other applications

Because the research topic has been very active in the last decade, sparse estimation and dictionary learning have been used in numerous ways in computer vision; unfortunately, writing an exhaustive review would be an endless exercise, and we have instead focused so far on a few approaches related to visual recognition in images. In the remaining paragraphs, we aim to be slightly more exhaustive and briefly mention a few other successful applications.

Action recognition in videos. We start with natural extensions to videos of visual recognition pipelines originally developed for images. [START_REF] Castrodad | Sparse modeling of human actions from motion imagery[END_REF] and [START_REF] Guha | Learning sparse representations for human action recognition[END_REF] have proposed two related approaches for action classification in movies. We present here the main principles from the method of [START_REF] Castrodad | Sparse modeling of human actions from motion imagery[END_REF], which has shown competitive results in several standard evaluation benchmarks.

Let us consider n training video sequences represented by vectors

x 1 , x 2 , . . . , x n in R lT , where T is the number of frames in the sequence and l is the number of pixels in each frame. Each sequence contains a subject performing a particular action and the goal is to predict the type of action in subsequent test videos. We assume that there are k different possible actions.

To address this classification problem, [START_REF] Castrodad | Sparse modeling of human actions from motion imagery[END_REF] have introduced a two-layer dictionary learning scheme, after preprocessing the videos to replace the raw pixels by temporal gradients. Then, the two layers of the pipeline can be described as follows:

• local model of actions: it is common in videos to work with three-dimensional patches with a time dimension to describe the local appearance of actions [START_REF] Laptev | On space-time interest points[END_REF]. For each class j, one dictionary D j is learned on the training data to represent the appearance of such patches within the action class j. This step is particularly useful when some actions are locally discriminative.

• encoding of full sequences: once the dictionaries D j are learned, a train or test video is processed by extracting its overlapping three-dimensional patches, and by encoding them with the concatenated dictionary D = [D 1 , . . . , D k ], say with p dictionary elements. This produces a representation α i for every patch i of the sequence. These patches are pooled into a single vector β in R p for the sequence.

• global model of actions: the first layer provides a highdimensional vector for each video. The second layer is exploiting these nonlinear representations of the sequences as inputs to the classification scheme (4.7). As a result, the layer is learning a global model for each action class, by learning again one dictionary per class. The difference with the first layer is that the dictionaries model the appearance of full sequences instead of local three-dimensional patches.

Note that, to simplify, we have also omitted some useful heuristics such as the removing of three-dimensional patches with low energy, corresponding to motionless regions. We have also not provided all mathematical details such as the chosen dictionary learning formulations and algorithms. We refer to [START_REF] Castrodad | Sparse modeling of human actions from motion imagery[END_REF] for a more exhaustive description of the method.

Visual tracking. We now move to a technique for visual tracking introduced by [START_REF] Mei | Robust visual tracking using ℓ 1 -minimization[END_REF], which relies on the sparsity-inducing effect of the ℓ 1 -norm and on some ideas from the robust face recognition system of [START_REF] Wright | Robust face recognition via sparse representation[END_REF]. Given a stream of images x 1 , x 2 , . . ., the goal of tracking is to find a bounding box around a moving object in the sequence. Typically, the bounding box is annotated for the first image, and the task is to find the boxes in the subsequent frames by exploiting temporally consistency.

The approach of [START_REF] Mei | Robust visual tracking using ℓ 1 -minimization[END_REF] is inspired by template matching techniques [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]Kanade, 1981, Matthews et al., 2004]. It assumes that at time t, a set of p templates D = [d 1 , . . . , d p ] in R m×p is available, where each column d j is a small image filled by the object of interest-in other words, a bounding box. Typically, each d j has been either obtained by a ground truth annotation at the beginning of the sequence, or is one of the estimated boxes in the previous frames. Then, estimating the position of the object in frame t consists of scanning all possible bounding boxes in that frame, and comparing how well they can be reconstructed by the templates D.

Concretely, let us assume that there are l possible boxes, and let us denote by R k the linear operator that extracts the k-th box from the image and subsamples it to make a representation with m parameters. Then, the optimization problem is simply (4.16) where e models the occlusion as in the face recognition formulation of [START_REF] Wright | Robust face recognition via sparse representation[END_REF] presented in Section 4.3. The ℓ 1 -norm encourages e to be sparse, whereas the weighted-ℓ 1 -norm also promotes the selection of a few templates from D. The role of the weights η j is to give more importance to some templates, especially the ones that have been selected recently in the previous frames.

k ∈ arg min k∈{1,...,l}   min α∈R p + ,e∈R m 1 2 R k x t -Dα -e 2 2 + λ e 1 + p j=1 η j |α[j]|   ,
Overall, the tracking scheme performs two successive operations for each frame. It selects the new bounding box with (4.16), and then updates the dictionary D and the weights η j , according to an ad-hoc procedure [see Mei and Ling, 2009, for the details]. According to a recent survey [START_REF] Wu | Online object tracking: A benchmark[END_REF], trackers based on sparse representations, which are essentially improved versions of [START_REF] Mei | Robust visual tracking using ℓ 1 -minimization[END_REF], perform well compared to the state of the art, and are able to deal with occlusion effectively.

Data visualization. Finally, we deviates from visual recognition by presenting applications of sparse estimation to data visualization in computer vision [START_REF] Elhamifar | See all by looking at a few: Sparse modeling for finding representative objects[END_REF][START_REF] Chen | Fast and robust archetypal analysis for representation learning[END_REF]. The first approach proposed by [START_REF] Elhamifar | See all by looking at a few: Sparse modeling for finding representative objects[END_REF] can be viewed as an extension of the sparse subspace clustering method of [START_REF] Elhamifar | Sparse subspace clustering[END_REF], and was successfully applied to the problem of video summarization.

More precisely, consider a video sequence X = [x 1 , . . . , x T ] in R l×T where each x t in R l is the t-th frame. The goal of video summarization is to find a set of s ≪ T frames that best "represents" the video. The selection problem is formulated as follows min (4.17) where A = [α 1 , . . . , α T ] and Ψ 0 (A) counts the number of rows of A that are non-zero. The constraints have two simultaneous effects: first, they encourage every frame x t to be close to a linear combination Xα t of other frames with coefficients α t that sum to one; second, since A has at most s non-zero rows, only s frames from X can be used in the approximations Xα t . The video sequence is consequently "summarized" by these s frames.

A∈R T ×T X -XA 2 F s.t. Ψ 0 (A) ≤ s and ∀i, T j=1 α i [j] = 1,
To visualize large image collections, [START_REF] Chen | Fast and robust archetypal analysis for representation learning[END_REF] has used related ideas by exploiting an older unsupervised learning technique called archetypal analysis, which is presented in Section 2.6. Specifically, let us consider a collection of images, denoted again by X = [x 1 , . . . , x n ], but which do not necessarily form a video sequence. Archetypal analysis looks for a dictionary D in R m×p that "represents" well the data, min α i ∈∆p for 1≤i≤n β j ∈∆n for 1≤j≤p

n i=1 x i -Dα i 2 2 s.t. ∀j, d j = Xβ j ,
where ∆ p and ∆ n are simplex sets-that is, sets of vectors with nonnegative entries that sum to one. In contrast to the original dictionary learning formulation, the dictionary elements, called archetypes, are convex combinations of data points and admit a simple interpretation.

In Figure 4.6, we present a result obtained by [START_REF] Chen | Fast and robust archetypal analysis for representation learning[END_REF] on a database of images downloaded on the Flickr website with the request "Paris". Each image is encoded using the Fisher vector representation [START_REF] Perronnin | Improving the Fisher kernel for large-scale image classification[END_REF] and p = 256 archetypes are learned. Each archetype d j is a convex combination Xβ j of a few input images with a sparse vector β j due to the simplicial constraint. For each archetype, it is the natural to visualize images that are used in the sparse combination, as well as the corresponding non-zero coefficients β j [i]. To some extent, the resulting figure "summarizes" the image collection. , we display some unexpected archetypes corresponding to soccer, graffitis, food, flowers and social gatherings; in (c), we see some archetypes that do not seem to have some semantic meaning, but that capture some scene compositions or texture present in the dataset. Best seen by zooming on a computer screen.

Optimization Algorithms

In the previous parts of the monograph, many formulations for dealing with image processing and computer vision tasks have been introduced, leading in general to cost functions to optimize. However, we have not presented yet in great details how the corresponding algorithms should be implemented in practice. The purpose of this section is to introduce basic optimization tools to address sparse estimation and dictionary learning problems. Because the literature is vast [see Bach et al., 2012a, for a review], we focus on a few techniques that we consider to be easy to use and efficient enough. Specifically, we mostly present algorithms that are parameter-free. This means that in addition to model parameters that are inherent to the formulation, they do not require tuning any input value to converge; they should at most require choosing a number of iterations or a stopping criterion. We also put the emphasis on algorithms that are easy to implement: a few lines of a high-level language, e.g., Matlab or Python, are in general sufficient to obtain an effective prototype, even though cutting-edge implementations may be more involved to develop-that is, they may require a low-level language such as C/C++ and optimized third-party libraries for linear algebra such as BLAS and LAPACK. 1First, we address classical ℓ 0 -regularized problems, before treating the case of the ℓ 1 -penalty. Then, we present iterative reweighted ℓ 1 and ℓ 2 techniques, and conclude with dictionary learning algorithms. Many of these methods have been used for the experiments conducted in this monograph, and are available in the SPAMS toolbox.2

Sparse reconstruction with the ℓ 0 -penalty

We start our review of optimization techniques with least square problems that are regularized with the ℓ 0 -penalty. More precisely, we are interested in finding approximate solutions of

min α∈R p x -Dα 2 2 s.t. α 0 ≤ k, (5.1) or min α∈R p α 0 s.t. x -Dα 2 2 ≤ ε, (5.2)
where x is an input signal in R m and D is a dictionary in R m×p . There are mainly two classes of algorithms for dealing with such NP-hard problems: greedy algorithms and iterative hard-thresholding methods.

We start with algorithms of the first kind, which are designed for dealing with both formulations. They are commonly known as matching pursuit and orthogonal matching pursuit in signal processing, and forward selection techniques in statistics.

Coordinate descent algorithm: matching pursuit (MP).

Introduced by [START_REF] Mallat | Matching pursuit in a time-frequency dictionary[END_REF], matching pursuit is a coordinate descent algorithm that iteratively selects one entry of the current estimate α at a time and updates it. It is closely related to projection pursuit algorithms from the statistics literature [START_REF] Friedman | Projection pursuit regression[END_REF], and is presented in Algorithm 2, where D = [d 1 , . . . , d p ] is a dictionary with unit-norm columns. Each iteration performs in fact a one dimensional minimization of the residual with respect to the -th entry of α Algorithm 2 Matching pursuit algorithm.

Require: Signal x in R m , dictionary D in R m×p with unit-norm columns, and stopping criterion k or ε.

1: Initialize α ← 0; 2: while α 0 < k if k is provided or x -Dα 2 2 > ε if ε is provided do 3:
select the coordinate with maximum partial derivative:

 ∈ arg max j=1...p |d ⊤ j (x -Dα)|; 4: update the coordinate α[] ← α[] + d ⊤  (x -Dα);
(5.3) 5: end while 6: return the sparse decomposition α in R p .

when keeping the other entries fixed. The update (5.3) can indeed be interpreted as

α[] ← arg min α∈R x - l = α[l]d l -αd  2 2 .
(5.4) Therefore, the residual monotonically decreases during the optimization, whereas the model sparsity α 0 increases or remains constant until the stopping criterion is satisfied. One drawback of matching pursuit is that the same entry  can be selected several times during the optimization, possibly leading to a large number of iterations. The next algorithm, called "orthogonal matching pursuit", does not suffer from this issue, and usually provides a better sparse approximation with additional computational cost.

Active-set algorithm: orthogonal matching pursuit (OMP). The algorithm is similar in spirit to matching pursuit, but enforces the residual to be always orthogonal to all previously selected variables, which is equivalent to saying that the algorithm reoptimizes the value of all non-zero coefficients once it selects a new variable. As a result, each iteration of OMP is more costly than those of MP, but the procedure is guaranteed to converge in a finite number of iterations. The algorithm appears in the signal processing literature [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], and is called "forward selection" in statistics (see Section 1.1). Several variants of OMP have been proposed that essentially differ in the way a new variable is selected at each iteration. In this monograph, we focus on an effective variant introduced by Gharavi-Alkhansari and [START_REF] Gharavi-Alkhansari | A fast orthogonal matching pursuit algorithm[END_REF], [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF], which is called "order recursive orthogonal matching pursuit". To simplify, we simply use the terminology "OMP" to denote this variant and omit the term "order recursive" in the rest of the monograph. The procedure is presented in Algorithm 3, where D Γ represents the matrix of size R m×|Γ| whose columns are those of D indexed by Γ, and Γ c denotes the complement set of Γ in {1, . . . , p}. The classical OMP algorithm (not the order recursive variant) chooses  exactly as in the MP algorithm. Through our experience, we have found that the order recursive variant provides a slightly better sparse approximation in practice, even though it may seem more complicated to implement at first sight.

Unlike MP, the set Γ of non-zero variables of α strictly increases by one unit after each iteration, and it is easy to show that the residual is always orthogonal to the matrix D Γ of previously selected dictionary elements. Indeed, we have at each step of the algorithm:

α Γ = D ⊤ Γ D Γ -1 D ⊤ Γ x,
(assuming D Γ to be full-rank), and thus

D ⊤ Γ (x -Dα) = D ⊤ Γ (x -D Γ α Γ ) = 0.
It is interesting to note that efficient implementations of Algorithm 3 exist, despite the fact that the algorithm naively requires to minimize |Γ c | quadratic functions-equivalently solve |Γ c | linear systemsat each iteration. The main ingredients for a fast implementation have been discussed by [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF], and are summarized below:

1. if the Cholesky factorization of (D ⊤ Γ D Γ ) is already computed at the beginning of an iteration, finding the index  to select in Algorithm 3 Orthogonal matching pursuit -order recursive variant.

Require: Signal x in R m , dictionary D in R m×p , and stopping criterion k or ε. 1: Initialize the active set Γ = ∅ and α ← 0;

2: while |Γ| < k if k is provided or x -Dα 2 2 > ε if ε is provided do 3:
select a new coordinate  that leads to the smallest residual:

(, β) ∈ arg min j∈Γ c ,β∈R |Γ|+1 x -D Γ∪{j} β 2 2 ;
(5.5)

4:

update the active set and the solution α:

Γ ← Γ ∪ {}; α Γ ← β and α Γ c ← 0;
5: end while 6: return the sparse decomposition α in R p .

the update (5.5) can be found with very few operations by using simple linear algebra relations [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF].

2. given some index , we have

β = (D ′⊤ Γ D ′ Γ ) -1 D ′⊤ Γ x, where Γ ′ = Γ ∪ {}. Whereas the matrix inversion naively costs O(|Γ ′ | 3 ) operations, it can be obtained in O(|Γ| 2 ) operations if the Cholesky factorization (D ⊤ Γ D Γ ) -1 is already computed. 3. if Q D ⊤ D
is precomputed, a significant gain in terms of speed can be obtained by exploiting the fact that the quadratic function x -Dα 2 2 expands as x 2 2 -2q ⊤ α + α ⊤ Qα. This trivial reformulation leads to an implementation that never explicitly computes the residual r = x -Dα, but instead works with the quantity D ⊤ r, initialized with q, and updated at every iteration. This simple strategy may be useful when large set of signals x 1 , . . . , x n needs to be encoded in parallel with the same dictionary. Then, it is worth computing Q once for all, before processing all signals. It is then possible to obtain an implementation of Algorithm 3 whose complexity for encoding these n signals is

O(mp 2 ) compute D ⊤ D once + nO(k 3 ) Cholesky factorizations + nO(pm + pk 2 ) update of D ⊤ r
.

A similar strategy for the classical OMP algorithm that involves a simpler selection rule of the variable  is discussed by [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF].

Even though the algorithm only provides an approximate solution of the ℓ 0 -regularized least square problem in general, OMP was shown by [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF] to be able to reliably recover a sparse signal from random measurements. Interestingly, the conditions that are required on the dictionary are similar to the ones for sparse recovery in the compressed sensing literature (see Section 1.5).

Gradient-descent technique: iterative hard-thresholding.

A significantly different approach than the previous greedy algorithms consists of using a gradient descent principle, an approach called iterative hard-thresholding [START_REF] Starck | Astronomical image representation by the curvelet transform[END_REF][START_REF] Herrity | Sparse approximation via iterative thresholding[END_REF][START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF], designed to find an approximate solution of either (5.1) or min (5.6) for some penalty parameter λ. The method is presented in Algorithm 4 and is an instance of the more general proximal gradient descent algorithms [see Bach et al., 2012a, for a review]. As we will see, when η is smaller than the inverse of the largest eigenvalue of D ⊤ D, the algorithm monotonically decreases the value of the objective function. To see this, we can interpret the iterative hard-thresholding algorithm as a majorization-minimization scheme [START_REF] Lange | Optimization transfer using surrogate objective functions[END_REF], which consists of minimizing at each iteration a locally tight upper-bound of the objective. This principle is illustrated in Figure 5 hard-thresholding:

α∈R p 1 2 x -Dα 2 2 + λ α 0 ,
α[j] ← α[j] if |α[j]| ≥ τ, 0 otherwise. 7:
end for 8: end for 9: return the sparse decomposition α in R p . The majorizing surrogate can be obtained from the following inequality that holds for all α and β in R p :

f (α) g β (α) β = α old α new f (α) ≤ g β (α)
1 2 x -Dα 2 2 = 1 2 x -Dβ -D(α -β) 2 2 = 1 2 x -Dβ 2 2 -(α-β) ⊤ D ⊤ (x-Dβ) + 1 2 D(α -β) 2 2 ≤ 1 2 x -Dβ 2 2 -(α-β) ⊤ D ⊤ (x-Dβ) + 1 2η α -β 2 2 = C β + 1 2η α -β + ηD ⊤ (x-Dβ) 2 2 g β (α) , ( 5.7)
where C β is a term independent of α. All equalities are simple linear algebra relations and the inequality comes from the fact that z ⊤ D ⊤ Dz ≤ 1 η z ⊤ z for all z in R p , since 1/η is assumed to be larger than the largest eigenvalue of D ⊤ D. Under this assumption, we have the relation 1 2 x -Dα 2 2 ≤ g β (α) for all α in R p , with an equality for α = β. Then, we have two possible cases:

• if the problem of interest is (5.1), k is provided to Algorithm 4, and it is easy to show that the value of α at the end of one iteration is the solution of min α∈R p g β (α) s.t. α 0 ≤ k, where β is the value of α at the beginning of the iteration.

• if one is interested in the penalized problem (5.6), the value of α after each iteration of Algorithm 4 is instead the solution of min α∈R p g β (α) + λ α 0 , and the function α → g β (α) + λ α 0 is a majorizing surrogate of the objective 1 2 x -Dα 2 2 + λ α 0 .

In both cases, Algorithm 4 can be interpreted as a majorizationminimization scheme and the value of the objective function monotonically decreases.

Sparse reconstruction with the ℓ 1 -norm

In the previous section, we have studied three optimization problems (5.1), (5.2) and (5.6). Similarly, we consider three formulations of interest involving the ℓ 1 -norm: (5.10)

min α∈R p 1 2 x -Dα 2 2 + λ α 1 , ( 5
We have also presented three optimization strategies to deal with the ℓ 0penalty, namely: matching pursuit, orthogonal matching pursuit, and iterative hard-thresholding algorithms. Interestingly, it is possible to deal with the ℓ 1 -penalty by following similar techniques, even though the links between ℓ 1 -approaches and ℓ 0 -ones are only clear from hindsight.

Coordinate descent. The counterpart of matching pursuit for ℓ 1regularized problems is probably the coordinate descent method. Even though the general scheme is classical in optimization [see [START_REF] Bertsekas | Nonlinear programming[END_REF], it was first proposed by [START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF] for minimizing the Lasso (5.8). Its convergence to the exact solution of the convex non-smooth problem such as (5.8) was shown by [START_REF] Tseng | A coordinate gradient descent method for nonsmooth separable minimization[END_REF], and its empirical efficiency was demonstrated by [START_REF] Wu | Coordinate descent algorithms for lasso penalized regression[END_REF]. The coordinate descent approach is presented in Algorithm 5, where S λ denotes the soft-thresholding operator already introduced in Section 1.2:

S λ : θ → sign(θ) max(|θ| -λ, 0),
and the selection rule for choosing  can be on the following strategies:

• cycle in {1, . . . , p} (the most typical choice);

• pick up  uniformly at random in {1, . . . , p};

Algorithm 5 Coordinate descent algorithm for the Lasso (5.8).

Require: Signal x in R m , dictionary D in R m×p , penalty parameter λ, number of iterations T , initial solution α 0 in R p . 1: Initialize α ← α 0 ; 2: for t = 1, . . . , T do 3: select one coordinate , e.g., by cycling or at random;

4:

update the coordinate by soft-thresholding:

α[] ← S λ α[] + d ⊤  (x -Dα) d  2 2
;

(5.11)

5: end for 6: return the sparse decomposition α in R p .

• choose the entry  exactly as in the matching pursuit algorithm.

It is relatively easy to show that the update (5.11) is in fact minimizing the objective (5.8) with respect to the -th entry of α when keeping the others fixed. More precisely, (5.11) is equivalent to

α[] ← arg min α∈R x - l = α[l]d l -αd  2 2 + λ|α|,
which is very similar to the matching pursuit update (5.4) but with the ℓ 1 -penalty in the formulation. The benefits of the coordinate descent scheme are three-fold: (i) it is extremely simple to implement; (ii) it has no parameter; (iii) it can be very efficient in practice [see Bach et al., 2012a, Section 8, for a benchmark].

Assuming that one chooses the cycling selection rule, the complexity of a naive update of α[j] is O(mp) for computing the matrix-vector multiplication Dα. Fortunately, better complexities can be obtained by adopting one of the two following alternative strategies:

1. during the algorithm, we choose to update an auxiliary variable representing the residual r such that we always have r = x -Dα. Extensions of coordinate descent to deal with more general smooth loss functions and group-Lasso penalties have also been proposed [START_REF] Tseng | A coordinate gradient descent method for nonsmooth separable minimization[END_REF]. These variants are also easy to implement, and exhibit good performance in practice [see Bach et al., 2012a].

Proximal gradient algorithm -iterative soft-thresholding. The second approach that we consider is the ℓ 1 -counterpart of the iterative hard-thresholding algorithm called "iterative softthresholding" [START_REF] Nowak | Fast wavelet-based image deconvolution using the EM algorithm[END_REF][START_REF] Figueiredo | An EM algorithm for wavelet-based image restoration[END_REF][START_REF] Starck | Astronomical image representation by the curvelet transform[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. Originally designed to address (5.8), it can be easily modified to deal with (5.10). It is similar to the hard-thresholding technique of the previous section and is presented in Algorithm 6, where η satisfies the same property as in the previous section-that is, 1/η is larger than the largest eigenvalue of D ⊤ D. Note that the orthogonal projection (5.12) onto the ℓ 1 -ball can be done efficiently in linear time [START_REF] Brucker | An O(n) algorithm for quadratic knapsack problems[END_REF][START_REF] Duchi | Efficient projections onto the l1-ball for learning in high dimensions[END_REF].

Then, the method can also be interpreted under the majorizationminimization point of view illustrated in Figure 5.1. Indeed, let us consider the two cases:

• if the problem of interest is (5.8), the user provides the parameter λ, and it is easy to show that the value of α at the end of one Algorithm 6 Proximal gradient algorithm for (5.8) or (5.10).

Require: Signal x in R m , dictionary D in R m×p , regularization parameter λ or µ, gradient descent step size η, number of iterations T , initial solution α 0 in R p (with α 0 1 ≤ T if T is provided). µ is provided and one performs the orthogonal projection α ← arg min

β∈R p 1 2 α -β 2 2 s.t. β 1 ≤ µ ;
(5.12) 8:

end if 9: end for 10: return the sparse decomposition α in R p . iteration is the solution of min α∈R p g β (α) + λ α 1 , where g β is defined in (5.7), and the function α → g β (α)+λ α 1 is a majorizing surrogate of the objective (1/2) x-Dα 2 2 +λ α 1 .

• if one is interested instead in the penalized problem (5.10), the value of α after each iteration of Algorithm 4 is the solution of min α∈R p g β (α) s.t. α 1 ≤ µ, where β is the value of α at the beginning of the iteration.

The majorization-minimization principle explains why the algorithm monotonically decreases the value of the objective function, but a more refined analysis can provide convergence rates. In general, we know that the objective function value converges to the optimal one with the convergence rate O(1/t), and that some variants such as FISTA admit a better rate O(1/t 2 ), as shown by [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. These methods are relatively easy to implement, and choosing η is not an issue in practice. When an appropriate value for η is not available beforehand, it is possible to automatically find a good one with a line-search scheme.

Similarly as coordinate descent, extensions to more general problems than ℓ 1 -regularized least squares are possible. In general proximal gradient methods are indeed adapted to solving convex problems of the form min

α∈A f (α) + ψ(α),
where A is a convex set in R p , f is differentiable and its gradient is uniformly Lipschitz continuous, and ψ is a convex nonsmooth function such that the so-called proximal operator of ψ

arg min α∈A 1 2 α -β 2 2 + ψ(α),
can be efficiently computed for all β in R p . As a result, proximal gradient methods have been successfully used with the Group Lasso penalty including the hierarchical and structured variants of Section 1.3 [see Jenatton et al., 2011b, Mairal et al., 2010b].

Homotopy and LARS algorithms Finally, the homotopy method [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF][START_REF] Efron | Least angle regression[END_REF] is related to orthogonal matching pursuit, but for solving the ℓ 1 -regularized problems (5.8), (5.9), and (5.10). Similarly to OMP, the algorithm maintains an active-set of variables Γ, initialized with the empty set, and iteratively updates it by one variable at a time.

The homotopy method follows in fact the regularization path of the Lasso-that is, it finds the set of all solutions of (5.8) for all values of the regularization parameter λ, formally defined as

P {α ⋆ (λ) : λ > 0},
where α ⋆ (λ) is the solution of (5.8), which is assumed to be unique under a small technical assumption [see Osborne et al., 2000, Mairal andYu, 2012]. A remarkable property of the path is that it can be shown to be piecewise linear. In other words, it can be described by a finite number k of linear segments, as illustrated earlier in Figure 1.4:

P = k l=1 {θα ⋆ (λ l ) + (1 -θ)α ⋆ (λ l+1 ) : θ ∈ [0, 1]} .
(5.13)

The homotopy method finds a point on the path, computes the direction of the current linear segment, and follows it until the direction of the path changes. The algorithm either returns the full regularization path-that is, the sequence of α ⋆ (λ l ) from Eq. (5.13), or stops if it reaches a solution for a target λ, a target residual ǫ, or ℓ 1 -norm µ. Thus, it can be used for solving the three variants (5.8), (5.9), or (5.10).

The piecewise linearity of the path was discovered by [START_REF] Markowitz | Portfolio selection[END_REF] in the context of portfolio selection, and was turned into an effective algorithm for solving the Lasso by [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF] and [START_REF] Efron | Least angle regression[END_REF]. It should be also noted that solving (5.8) for all values of λ is in fact an instance of parametric quadratic programming, for which path following algorithm encompassing the homotopy method appear early in the optimization literature [START_REF] Ritter | Ein verfahren zur lösung parameterabhängiger, nichtlinearer maximum-probleme[END_REF]. The homotopy turns out to be extremely efficient for solving very sparse medium-scale problems [see Bach et al., 2012a], even though its worstcase complexity has been shown to be exponential [START_REF] Mairal | Complexity analysis of the Lasso regularization path[END_REF]. Similar to the simplex algorithm for solving linear programs, the empirical complexity is therefore extremely bad in the worst-case scenario, but it is most often very good in practical ones.

Following the presentation of [START_REF] Mairal | Complexity analysis of the Lasso regularization path[END_REF], we describe the method in Algorithm 7. The formula for the direction of the path is given in Eq. (5.14), where the quantity η is defined as η sign(D ⊤ (x -Dα ⋆ )) in {-1, 0, 1} p . This relation can be derived from optimality conditions of the Lasso presented in Section 1.3, but the derivation requires technical developments that are beyond the scope of the monograph. We refer instead to [START_REF] Mairal | Complexity analysis of the Lasso regularization path[END_REF] for more details. Interestingly, implementing efficiently Algorithm 7 raises similar difficulties as for implementing orthogonal matching pursuit. We need indeed to be able to update the inverse matrices (D ⊤ Γ D Γ ) -1 at each iteration when the active-set Γ changes. The Cholesky decomposition or the Woodbury formula [see [START_REF] Golub | Matrix computations[END_REF]] can be used for that purpose.

For example, let us give a few details on how the Woodbury formula can be useful. Denoting by Γ ′ Γ ∪ {}, we have the classical relation

D ⊤ Γ ′ D Γ ′ -1 = (D ⊤ Γ D Γ ) -1 + 1 s zz ⊤ -1 s z -1 s z ⊤
Algorithm 7 Homotopy algorithm for the Lasso.

Require: Signal x in R m , dictionary D in R m×p . 1: λ ← D ⊤ x ∞ ; 2: Γ ← {j ∈ {1, .
. . , p} : |d ⊤ j x| = λ}; 3: while λ > 0 or according to another stopping criterion do 4:

compute the current direction of the regularization path:

α ⋆ Γ (λ) = (D ⊤ Γ D Γ ) -1 (D ⊤ Γ x -λη Γ ) and α ⋆ Γ c (λ) = 0. (5.14) 5:
find the smallest τ > 0 such that one of the following events occurs:

• there exists j in Γ c such that |d ⊤ j (x-Dα ⋆ (λ-τ ))| = λ-τ ; then, add j to Γ;

• there exists j in Γ such that α ⋆ (λτ )[j] hits zero; then, remove j from Γ.

It is also easy to show that the value of τ can be obtained in closed form such that one can "jump" from a kink to another;

6:

update λ ← λτ ; record the pair (λ, α ⋆ (λ)); 7: end while 8: Return: sequence of recorded values (λ, α ⋆ (λ)).

These regularization functions typically have the form For example, for minimizing min

ψ(α) = p j=1 ϕ (|α[j]|) , ( 5
α∈R p 1 2 x -Dα 2 2 + λ p j=1 log(|α[j]| + ε),
the resulting majorization-minimization algorithm solves a sequence of weighted Lasso problems at each iteration (5.16) where α old denotes the current estimate at the beginning of the iteration. The update (5.16) is a reweighted-ℓ 1 algorithm, which is known to provide sparser solutions than the regular Lasso. Such majorization-minimization approaches for non-convex sparse estimation have been introduced in various contexts. For example, [START_REF] Fazel | Matrix rank minimization with applications[END_REF], [START_REF] Fazel | Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices[END_REF] and then later Candès et al. [2008] use such a principle for the regularization function ϕ : x → log(x + ε), whereas [START_REF] Figueiredo | A bound optimization approach to wavelet-based image deconvolution[END_REF], [START_REF] Figueiredo | Majorizationminimization algorithms for wavelet-based image restoration[END_REF] propose majorization-minimization algorithms for dealing with the ℓ q -penalty with q < 1.

α new ← arg min α∈R p   1 2 x -Dα 2 2 + λ p j=1 |α[j]| |α old [j]| + ε   ,

Iterative reweighted-ℓ 2 methods

We now focus on another type of algorithm for dealing with convex regularizers such as the ℓ 1 -norm. Similar to the reweighted-ℓ 1 -algorithms of the previous section, the approaches we present here consist of solving a sequence of sub-problems, each one involving a simple penalty.

The key idea is that many non-smooth convex regularizers can be seen as minima of quadratic functions, leading to reweighted ℓ 2 -algorithms.

Variational formulation of the ℓ 1 -norm. We start with the simplest case of the ℓ 1 -norm, which admits the following variational form [START_REF] Grandvalet | Outcomes of the equivalence of adaptive ridge with least absolute shrinkage[END_REF]Canu, 1999, Daubechies et al., 2010]:

α 1 = inf η∈R p + 1 2 p j=1 α[j] 2 η[j] + η[j] ,
with the optimal η ∈ R Fortunately, this problem is jointly convex in (α, η). The minimization with respect to η given α can be done in closed form (by taking absolute values of the components of α); the minimization with respect to α given η is a weighted least-squares problem, which can be solved by the solution of a linear system as:

α = D ⊤ D + λDiag(η) -1 -1 D ⊤ x.
Note that since the linear system above is ill-conditioned when some entries of η are very small, it is better to equivalently reformulate it as α = Diag(η) 1/2 Diag(η) 1/2 D ⊤ DDiag(η) 1/2 + λI -1 Diag(η) 1/2 D ⊤ x.

This naturally leads to alternating minimization procedures where one alternates between minimizing with respect to η as above in closed form, and minimizing with respect to α, which is now easier because the regularizer has become a quadratic function. With quadratic losses, this leads to a linear system for which many existing algorithms exist. However, because the function (α[j], η[j]) → α[j] 2 /η[j] cannot be extended to be continuous at (0, 0), alternating minimization does not converge and it is customary to add a penalty term p j=1 ε 2η[j] with ε > 0 small, which leads to an update for η[j] of the form η[j] = |α[j]| 2 + ε that avoids instabilities. See more details in Bach et al. [2012a, Section 5].

Extensions to ℓ q -quasi-norms. More generally, for any q ∈ (0, 2), and α q = q j=1 |α[j]| q 1/q and r = q/(2q), we have [Jenatton et al., 2010b]:

α q = inf η∈R p + 1 2 p j=1 α[j] 2 η[j] + 1 2 η r ,
with the optimal η in R p + obtained as η[j] = |α[j]| 2-q α q-1 q . Note that similar variational formulations may be obtained for the squared quasi-norms [START_REF] Micchelli | Learning the kernel function via regularization[END_REF]. This gives an alternative to reweighted ℓ 1 -formulations presented in Section 5.3 for q < 1 (the nonconvex case).

Extensions to group norms. The various formulations above can be extended to structured sparsity-inducing norms, such as the ones presented in Section 1.3. For example, for the group-lasso norm in Eq. (1.20), we have:

g∈G α[g] 2 = inf η∈R p + 1 2 g∈G α[g] 2 2 η[g] + η[g] ,
which also leads to a quadratic function of α and to simple algorithms based on solving linear systems when the loss is quadratic. Note that these quadratic variational formulations extend to all norms and lead to multiple kernel learning formulations when used with positive-definite kernels [see more details in Bach et al., 2012a].

Optimization for dictionary learning

Finally, we review a few dictionary learning algorithms that leverage the optimization algorithms for sparse estimation presented in the previous sections. We start with the stochastic gradient descent method [see [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]Yin, 2003, Bottou and[START_REF] Bottou | The tradeoffs of large scale learning[END_REF]references therein] since it is very close to the original algorithm of [START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by V1?[END_REF]. Then, we move to other classical approaches such as alternate minimization [START_REF] Engan | Method of optimal directions for frame design[END_REF], Lee et al., 2007], block coordinate descent, K-SVD [START_REF] Aharon | An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], and then present efficient methods based on stochastic approximations [Mairal et al., 2010a, Skretting and[START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF]. We consider two dictionary learning formulations. For both of them, the goal is to learn a dictionary D in C, where C is the set of matrices in R m×p whose columns are constrained to have less than unit ℓ 2 -norm. Then, we use a sparsity-inducing penalty ψ either as a penalty with empirical good performance for various tasks, e.g., image denoising. In general, these methods have no other guarantee than providing a stationary point (in the best case), but they have been used successfully in practical contexts, as shown in other parts of this monograph.

Stochastic gradient descent. The first algorithm proposed by Olshausen andField [1996, 1997] addresses the formulation (5.17), where ψ is the ℓ 1 -norm or a smooth approximate sparsity-inducing penalty such as ψ(α) = p j=1 log(α[j] 2 + ε). The algorithm can be described as a heuristic stochastic gradient descent, which alternates between two stages: (i) gradient steps with a fixed step size computed on mini-batches of the training set and (ii) rescaling heuristic that prevents the norm of the columns d j to grow out of bounds-thus, taking implicitly the constraint D ∈ C into account.

A less heuristic but very related procedure is the projected stochastic gradient descent method, which is presented in Algorithm 8. This approach relies on the following equivalent formulation to (5.17 Then, the algorithm performs a Euclidean projection Π C onto the convex set C, which corresponds to the renormalization update d j → (1/ max( d j 2 , 1))d j for all j in {1, . . . , p}. A practical variant also consists of using mini-batches instead of drawing a single data point at each iteration. For simplicity, we omit this variant in our description of Algorithm 8. The stochastic gradient descent approach is effective [see Mairal et al., 2010a, for a benchmark], but it raises a few practical difficulties such as choosing well the step sizes η t in a data-independent manner. Classical strategies use for instance step sizes of the form η t = η/(t + t 0 ) γ , but finding an appropriate set of hyper-parameters η, t 0 , and γ may be challenging in practice.

Algorithm 8 Stochastic gradient descent for (5.17 Alternate minimization. The most classical approach for dictionary learning is the alternate minimization scheme, which consists of optimizing (5.17) or (5.18) by alternating between two minimization steps: one with respect to D with the sparse codes A fixed, and one with respect to A with D fixed. In practice, this approach is not as fast as a well-tuned stochastic gradient descent algorithm, but it is parameterfree and we subjectively find it very reliable according to our experience for image processing and computer vision tasks. Alternate minimization was first proposed by [START_REF] Engan | Method of optimal directions for frame design[END_REF] for dealing with the ℓ 0 -penalty with the method of optimal directions (MOD), and was revisited later by [START_REF] Lee | Efficient sparse coding algorithms[END_REF] with the ℓ 1regularization. We present the MOD approach in Algorithm 9 and its ℓ 1 -counterpart in Algorithm 10. As we shall see, they slightly differ in the dictionary update step, which can be simplified for the ℓ 0 -penalty.

The sparse codes α i for MOD are obtained with any algorithm for dealing approximately with ℓ 0 , e.g., any technique from Section 5.1. Since all the vectors α i can be computed in parallel, it is worth mentioning that this sparse decomposition step can benefit from an efficient implementation of OMP including the heuristics of Section 5.1, notably those consisting of pre-computing D ⊤ D before updating all the α i 's in parallel. For dealing with ℓ 1 , the homotopy method can be used similarly, with the same heuristics that apply to OMP. Moving on now to the dictionary update step, we remark that (5.22) and (5.23) are not the same and are not equivalent to each other, even though, at first sight, the optimization problem with respect to D when A is fixed seems to be independent of the regularization ψ. The MOD update indeed minimizes the least-square function (5.23) by discarding the constraint D in C, yielding a solution XA ⊤ (AA ⊤ ) -1 (assuming AA ⊤ to be invertible), before projecting it onto the constraint set C. In constrained optimization [see [START_REF] Bertsekas | Nonlinear programming[END_REF], it is usually not update α i by solving (5.8) with x = x i if λ is provided, or (5.10) if µ is provided. x i -Dα i 2 2 ;

(5.23)

8: end for 9: return the dictionary D in C.

appropriate to discard a constraint, solve the resulting unconstrained formulation, before projecting back the solution onto the constraint set. In general, such a procedure can be indeed arbitrarily bad regarding the original constrained problem. In the context of ℓ 0 -regularized dictionary learning, however, such an approach does fortunately make sense. The objective function (5.18) when ψ = ℓ 0 is indeed invariant to a rescaling operation consisting of multiplying each column of d j by a scalar γ j and the corresponding j-th row of A by (1/γ j ). As a consequence, the minimum value of the objective function is independent of the scale of the columns of D; thus the constraint D ∈ C can always be safely ignored and the normalization D ← Π C [D] can be applied at any moment without changing the quality of the dictionary. In the case of ℓ 1 , the objective function is not invariant to the previous rescaling operation and the MOD update for the dictionary should not be used anymore. Instead, the objective function (5.23) needs to be minimized with another approach. Since it is convex, [START_REF] Lee | Efficient sparse coding algorithms[END_REF] propose to use a Newton method in the dual of (5.23), which requires solving p×p linear systems at every Newton iteration. Another easy-toimplement choice is the block coordinate descent update of Mairal et al. [2010a], which we present in Algorithm 11. Each step of the block coordinate descent approach is guaranteed to decrease the value of the objective function, while keeping D in the constraint set C, and ultimately solve (5.23) since the problem is convex [see [START_REF] Bertsekas | Nonlinear programming[END_REF], for convergence results of block coordinate descent methods]. In practice, (5.23) does not need to be solved exactly and performing a few passes, say 10, over the columns of D seems to perform well in practice.

The update (5.24) can be justified as follows. Minimizing (5.23) with respect to one column d j when keeping the other columns fixed can be formulated as

d j ← arg min d∈R m , d 2 ≤1    n i=1 1 2 x i - l =j α i [l]d l -α i [j]d 2 2    .
Then, this formulation can be rewritten in a matrix form (5.24)

5:

end for 6: until convergence; 7: return D (updated dictionary).

Block coordinate descent.

We have presented in the previous paragraph the alternate minimization method where the dictionary is updated via block coordinate descent. We have also introduced in Section 5.2 a coordinate descent algorithm for dealing with ℓ 1 -penalized problems. It is then natural to combine the two approaches into a block coordinate descent dictionary learning method for (5.17) with ψ = ℓ 1 . We present it in Algorithm 12, where S λ is the soft-thresholding operator applied element-wise to a vector. The main asset of the blockcoordinate descent approach is its simplicity, and ease of implementation. In practice, it can also be very effective.

K-SVD.

For ℓ 0 -regularized dictionary learning problems, one of the most popular approach is the K-SVD algorithm of [START_REF] Aharon | An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. The algorithm is related to the alternate minimization strategy MOD [START_REF] Engan | Method of optimal directions for frame design[END_REF], but the dictionary update step updates the non-zero coefficients of the matrix A at the same time. The full approach is detailed in Algorithm 13. The sparse encoding approach is the same as in MOD, but the dictionary elements are updated in a block coordinate fashion, one at a time.

Algorithm 12 Block coordinate descent for (5.17 update the dictionary D: perform one pass or more of the inner loop of Algorithm 11. 8: end for 9: return the dictionary D in C.

More precisely, when updating the column d j , the algorithm finds the subset of signals Ω that use d j in their current sparse decomposition. Then, the non-zero coefficients of the j-th row of A, denoted by α j Ω , are updated at the same time as d j in (5.26). Such an update is in fact equivalent to performing a rank-one singular value decomposition, which has inspired the name of the algorithm. A typical technique for solving (5.26) is the power method [see [START_REF] Golub | Matrix computations[END_REF], which alternates between the optimization of d with β fixed, and of β with d fixed.

Online dictionary learning. Finally, we describe the online dictionary learning method of Mairal et al. [2010a], which we have used in the experiments of this monograph. The approach is based on stochastic approximations, where the goal is to minimize the expectation min D∈C E x [L(x, D)],

(5.27)

Algorithm 13 K-SVD for (5.18) with ψ = ℓ 0 .

Require: Signals X = [x 1 , . . . , x n ] in R m×n , initial dictionary D 0 in C, regularization parameter µ, number of iterations T . 1: Initialize D ← D 0 ; 2: for t = 1, . . . , T do 3: compute the sparse codes, e.g., with OMP: end for 12: end for 13: return the dictionary D in C.

where L is defined in (5.20), and the expectation is taken with respect to the unknown probability distribution of the data x. With a discrete probability distribution with a finite sample, (5.27) simply amounts to (5.19).

In a nutshell, the algorithm works by using a stochastic version of the majorization-minimization principle [START_REF]Stochastic majorization-minimization algorithms for large-scale optimization[END_REF]. It sequentially builds a quadratic surrogate of the expected cost (5.27), which is minimized at every iteration. It can be shown that the quality of approximation of the surrogate near the current estimate improves during the algorithm iterates, and that the approach provides asymptotically a stationary point of (5.27). The basic procedure is summarized in Algorithm 14 and some variants that significantly improves its practical efficiency are presented in [Mairal et al., 2010a].

Assuming that the training set is composed of i.i.d. samples of a data distribution, the inner loop draws one element x t at a time, as in the stochastic gradient descent algorithm, and alternates between classical sparse coding steps for computing the decomposition vector α t of x t over the current dictionary D t-1 , and dictionary update steps. The new dictionary D t is computed by minimizing over C the quadratic surrogate (5.28) where the vectors α i for i < t have been computed during the previous steps of the algorithm. The motivation of the approach is based on two ideas:

• the quadratic surrogate aggregates the past information with a few sufficient statistics, namely the matrices B and C. One key component of the convergence analysis shows that the difference between the quadratic surrogate and the objective function at the current estimate converges to zero almost surely.

• with warm restart, the dictionary update is very efficient. In practice, one pass over the dictionary elements seems to perform empirically well.

Variants to improve the efficiency of the online dictionary learning algorithm include: (i) the use of mini-batches instead of drawing a single point at every iteration, (ii) rescaling the past data by updating the matrices B ← γ t B + x t α ⊤ t , with a sequence of parameters γ t that converges to one. The implementation of the SPAMS software contains such variants, which have proven to be important in practice [see Mairal et al., 2010a, for additional discussions].

Algorithm 14 Online dictionary learning for (5.17 with ψ = ℓ 1 .

Require: i.i.d signals x 1 , x 2 , . . . ∈ R m , regularization parameter λ, initial dictionary D 0 ∈ C, number of iterations T . 1: Initialization: B ∈ R m×p ← 0, C ∈ R p×p ← 0, D ← D 0 ; 2: for t = 1, . . . , T do 3: draw x t and sparsely encode it using homotopy: Note that along the same lines, a similar algorithm called "recursive least square dictionary learning" was independently developed by [START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF] focusing on the ℓ 0 -penalty.

α t ∈ arg min
Other extensions. Because the literature on dictionary learning algorithms is vast, we have presented a few ones that address the formulations (5.17) or (5.18). Extending these approaches to other settings, such as other loss functions than the square loss, or other regularization functions, is easy in general by replacing the sparse encoding and dictionary update steps by ad hoc optimization techniques.

Other optimization techniques

For simplicity, we have focused on sparse estimation problems involving least-square cost functions, and have left aside structured sparsity penalties such as the group Lasso. Some of the methods we have described can be easily extended to these settings, in particular the proximal gradient and coordinate descent algorithms are relatively flexible. We leave these extensions to more exhaustive reviews [see Bach et al., 2012a, for instance].

Other classes of algorithms are also popular for solving convex sparse estimation problems. For example, proximal splitting techniques [START_REF] Combettes | Fixed-Point Algorithms for Inverse Problems in Science and Engineering, chapter Proximal Splitting Methods in Signal Processing[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF][START_REF] Raguet | A generalized forward-backward splitting[END_REF], such as the alternating direction method of multipliers [see [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], are often used in signal processing. These methods are very flexible, simple to implement, and can be applied to a large range of problems. However, they require tuning (at least) one parameter to obtain fast convergence in practice, which, to the best of our knowledge, cannot be done automatically. In the context of machine learning, stochastic optimization techniques have also been developed [e.g. [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF][START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF]. These methods are adapted to large-scale problems; even though, they may seem simple at first sight, they require a careful choice of a step size parameter and thus remain a bit difficult to use.

Finally, we have also omitted quasi-Newton methods such as variants of L-BFGS [START_REF] Schmidt | Projected Newton-type methods in machine learning[END_REF], which can exploit the curvature of the objective function to achieve fast convergence.

In this monograph, we have reviewed a large number of applications of dictionary learning in image processing and computer vision and presented basic sparse estimation tools. We started with a historical tour of sparse estimation in signal processing and statistics, before moving to more recent concepts such as sparse recovery and dictionary learning. Subsequently, we have shown that dictionary learning is related to matrix factorization techniques, and that it is particularly effective for modeling natural image patches. As a consequence, it has been used for tackling several image processing problems and is a key component of many state-of-the-art methods in visual recognition. We have concluded the monograph with a presentation of optimization techniques that should make dictionary learning easy to use for researchers that are not experts of the field.

With efficient software available, it is remarkable to see that the interest for dictionary learning and sparse estimation is still acute today, and seems to increase in other scientific communities than image processing or computer vision. In particular, new applications of sparse matrix factorization are found in neuroscience [START_REF] Varoquaux | Multisubject dictionary learning to segment an atlas of brain spontaneous activity[END_REF], bioinformatics for processing gene expression data [see Barillot et al., 2012, Section 5.4], audio processing [START_REF] Baby | Coupled dictionary training for exemplar-based speech enhancement[END_REF], and in astrophysics [START_REF] Vinci | Estimating the distribution of galaxy morphologies on a continuous space[END_REF]. We have also put the emphasis on natural images, and we have omitted the application of dictionary learning to other domains, such as hyperspectral data, where significant success have been obtained [START_REF] Charles | Learning sparse codes for hyperspectral imagery[END_REF], Xing et al., 2012], or to disparity maps for stereo imaging [START_REF] Tosic | Learning sparse representations of depth[END_REF].

Exploring the possible use of dictionary learning in these fields would be of great interest and high importance, but unfortunately this goal falls beyond the scope of this monograph.

  Soft-thresholding operator, α st = sign(β) max(|β| -λ, 0).

  Hard-thresholding operator α ht = 1 |β|≥µ β.

Figure 1 . 1 :

 11 Figure 1.1: Soft-and hard-thresholding operators, which are commonly used for signal estimation with orthogonal wavelet basis.

Figure 1 . 2 :

 12 Figure1.2: Illustration of a wavelet tree with four scales for one-dimensional signals. Nodes represent wavelet coefficients and their depth in the tree correspond to the scale parameter of the wavelet. We also illustrate the zero-tree coding scheme[START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF] in this figure. Empty nodes correspond to zero coefficient: according to the zero-tree coding scheme, their descendants in the tree are also zero.

Figure 1 . 3 :

 13 Figure 1.3: Wavelet coefficients displayed for the image lena using the orthogonal basis of Daubechies [1988]. A few coefficients representing a low-resolution version of the image are displayed on the top-left corner. Wavelets corresponding to this low-resolution image are obtained by filtering the original image with shifted versions of a low-pass filter called "scaling function" or "father wavelet". The rest of the coefficients are organized into three quadtrees (on the right, on the left, and on the diagonal). Each quadtree is obtained by filtering the original image with a wavelet at three different scales and at different positions. The value zero is represented by the grey color; negative values appear in black, and positive values in white. The wavelet decomposition and this figure have been produced with the software package matlabPyrTools developed by Eero Simoncelli and available here: http://www.cns.nyu.edu/~lcv/software.php.

Figure 1 . 4 :

 14 Figure 1.4:Two examples of regularization paths for the Lasso/Basis Pursuit. The curves represent the values of the p = 5 entries of the solutions of (1.10) when varying the parameter λ for two datasets. On the left, the relation between λ and the sparsity of the solution is monotonic; On the right, this is not the case. Note that the paths are piecewise linear, see Section 5.2 for more details.

Figure 1 . 5 :

 15 Figure1.5: A physical system illustrating the sparsity-inducing effect of the ℓ1norm (on the right) in contrast to the Tikhonov-ridge regularization (on the left). Three springs are represented in each figure, two on the left, one on the right. Red points are fixed and cannot move. On the left, two springs are linked to each other by a blue point whose position can vary. On the right, a blue object of mass m is attached to the spring. Right and left configurations define two different dynamical systems with energies E1+E2; on the left, E1 and E2 are elastic potential energies; on the right, E1 is the same as on the left, whereas E2 is a gravitational potential energy, where g is the gravitational constant on the Earth's surface. Both system can evolve according to their initial positions, and stabilize for the value of α ⋆ that minimizes the energy E1 + E2, assuming that some energy can be dissipated by friction forces. On the left, it is possible to show that α ⋆ = βk1/(k1 + k2) and thus, the solution α ⋆ is never equal to zero, regardless of the strength k2 of the bottom spring. On the right, the solution is obtained by soft-thresholding: α ⋆ = max(β -mg/k1, 0). As shown on Figure1.5(b), when the mass m is large enough, the blue object touches the ground and α ⋆ = 0. Figure adapted from[START_REF]Sparse coding for machine learning, image processing and computer vision[END_REF].

  Effect of the Euclidean projection onto the ℓ2-ball.

Figure 1 . 6 :

 16 Figure 1.6: Illustration in two dimensions of the projection operator onto the ℓ1ball in Figure (a) and ℓ2-ball in Figure (b).The balls are represented in gray. All points from the red regions are projected onto the point of coordinates (0, µ) denoted by a red dot. Similarly, the green and blue regions are projected onto the green and blue dots, respectively. For the ℓ1-norm, a large part of the figure is filled by the red and green regions, whose points are projected to a sparse solution corresponding to a corner of the ball. For the ℓ2-norm, this is not the case: any non-sparse point-say, for instance on the blue line-is projected onto a non-sparse solution.

Figure 1 . 7 :

 17 Figure 1.7: Representation in three dimensions of the ℓ1-and ℓ2-balls. Figure borrowed from Bach et al. [2012a], produced by Guillaume Obozinski.

Figure 1 . 8 :

 18 Figure 1.8: Similar illustration as Figure 1.6 for the ℓ∞-norm. The regularization effect encourages solution to be on the corners of the ball, corresponding to points with the same magnitude |α[1]| = |α[2]| = µ.

Figure 1 . 11 :

 111 Figure1.11: Illustration of the sparsity, group sparsity, and hierarchical sparsity principles. Each column represents the sparsity pattern of a vector with 12 variables. Non-zero coefficients are represented by gray squares. On the left (a), the vectors are obtained with a simple sparsity-inducing penalty, such as the ℓ1-norm, and the non-zero variables are scattered. In the middle figure (b), a group sparsity-inducing penalty with three groups of variables is used. The sparsity patterns respect the group structure. On the right (c), we present the results obtained with a hierarchical penalty consisting of the Group Lasso plus the ℓ1-norm. The group structure is globally respected, but some variables within a group can be discarded.

  (a) Group-Lasso penaltyψ(α) = α[1, 2] 2 + |α[3]|. (b) Hierarchical penalty ψ(α) = α 2 + |α[1]| + |α[2]|.(c) Structured sparse penalty. (d) Structured sparse penalty.

Figure 1 . 13 :

 113 Figure 1.13: Visualization in three dimensions of unit balls corresponding to various sparsity-inducing norms. (a): Group Lasso penalty; (b): hierarchical penalty of Zhao et al. [2009]; (c) and (d): examples of structured sparsity-inducing penalties of Jacob et al. [2009]. Figure borrowed from Bach et al. [2012a], produced by Guillaume Obozinski.

Figure 1 . 14 :

 114 Figure 1.14: An illustration of the projective interpretation of sparse coding.

  (a) 2D Gabor filter. (b) With shifted phase. (c) With rotation.

Figure 2 . 1 :

 21 Figure 2.1: Three examples of two-dimensional Gabor filters. Red represents positive values and blue negative ones. The first two filters correspond to the function g(x, y) = e -x 2 /(2σ 2 x )-y 2 /(2σ 2 y ) cos(ωx + ϕ), with ϕ = 0 for Figure (a) and ϕ = π/2 for Figure (b).Figure (c) was obtained by rotating Figure (a)-that is, replacing x and y in the function g by x ′ = cos(θ)x + sin(θ)y and y ′ = cos(θ)x -sin(θ)y. The figure is best seen in color on a computer screen.

  Figure (c) was obtained by rotating Figure (a)-that is, replacing x and y in the function g by x ′ = cos(θ)x + sin(θ)y and y ′ = cos(θ)x -sin(θ)y. The figure is best seen in color on a computer screen.

  After centering and ℓ2-normalization. (d) After whitening.

Figure 2 . 2 :

 22 Figure 2.2: The effect of various pre-processing procedures on the image man. For the figures (b), (c) and (d), the value 0 is represented by a gray pixel. Since they contain negative values, pre-processed images are shifted and rescaled for visualization purposes.

  After centering and ℓ2-normalization. (d) After whitening.

Figure 2 . 3 :

 23 Figure 2.3: The effect of various pre-processing procedures on the image kodim07. Pre-processed images are shifted and rescaled for visualization purposes, since they contain negative values. The figure is best seen in color on a computer screen.

Figure 2 . 4 :

 24 Figure 2.4: On the right, we visualize the principal components of 400 000 randomly sampled natural image patches of size 16 × 16, ordered by decreasing variance, from top to bottom and left to right. On the left, we display a discrete cosine transform (DCT) dictionary. Principal components resemble DCT dictionary elements.

Figure 2 . 5 :

 25 Figure 2.5: Visualization of the principal components of all overlapping patches from the image tiger. Even though the image is not natural, its principal components are similar to the ones of Figure 2.4(b).

  Figure 2.4(b).

  (a) With centering. (b) With whitening.

Figure 2 . 6 :

 26 Figure 2.6: Visualization of p = 256 centroids computed with the algorithm Kmeans on n = 400 000 image patches of size m = 16 × 16 pixels. We compare the results with centering (left), and whitening (right). In both cases, we also apply a contrast normalization step. Centroids are ordered from most to least used (from left to right, then top to bottom).

  (a) With centering -gray. (b) With centering -RGB. (c) With whitening -gray. (d) With whitening -RGB.

Figure 2 . 7 :

 27 Figure 2.7: Dictionaries obtained with the formulation (2.2), with two different pre-processing procedures, and with gray or RGB patches. The dictionary elements are ordered from most to least used (from left to right, then top to bottom). Best seen in color on a computer screen.

Figure 2 . 8 :

 28 Figure 2.8: Dictionaries obtained with the formulation (2.2) and a hierarchical structured-sparsity penalty Ψ for natural image patches of size 16 × 16 pixels. The tree structures are pre-defined, and the dictionary elements naturally organize themselves in the tree during learning. For each tree, the root is represented in the middle of the figure. (a): the tree is of depth 4 and the branching factors at depths 1, 2, 3 are respectively 10, 2, 2; (b): the tree is slightly more complex; its depth is 5 and the branching factors at depths 1, 2, 3, 4 are respectively 10, 2, 2, 2. Figures borrowed fromJenatton et al. [2010aJenatton et al. [ , 2011b]]. Best seen by zooming on a computer screen.

  (a) With 3 × 3 neighborhoods. (b) With 4 × 4 neighborhood.

Figure 2 . 9 :

 29 Figure 2.9: Dictionaries obtained with the formulation (2.2) and a structured sparsity penalty Ψ inducing a grid structure in the dictionary. The dictionaries are computed on whitened natural image patches of size 12 × 12 pixels for two different group structures. Figure borrowed from Mairal et al. [2011]. Best seen on a computer screen.

  (a) With λ = 0.2.

Figure 2 . 10 :

 210 Figure 2.10: Visualization of p = 100 dictionary elements learned on 30 whitened natural images. Dictionary elements are ordered from most to least used (from left to right, then top to bottom).

Figure 3 . 1 :

 31 Figure 3.1: Dataset of 12 standard images used in the image denoising benchmarks.

Figure 3 . 2 :

 32 Figure 3.2: Top: inpainting result obtained by using the source code of Mairal et al.[2008c], where the text is automatically removed on the restored image. Bottom: inpainting result fromMairal et al. [2008e], where 80% of the pixels are randomly missing from the original image. The algorithm is able to reconstruct the brick texture on the right, by adapting the dictionary to the available information displayed on the left. The bottom images are under "copyright c 2008 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved". Best seen by zooming on a computer screen.

Figure 3 . 3 :

 33 Figure 3.3: Common Bayer pattern used for digital camera sensors.

Figure 3 . 4 :

 34 Figure 3.4: (a): Mosaicked image; (b): image demosaicked with the approach of Mairal et al. [2008c] presented in Section 3.3. (c): image demosaicked with the non-local sparse image models of[START_REF] Mairal | Non-local sparse models for image restoration[END_REF]. Note that the part of the image displayed here in known to be particularly difficult to reconstruct. Other reconstructions obtained from the database of 24 images is artifact-free in general. Best seen in color by zooming on a computer screen.

Figure 3 . 5 :

 35 Figure 3.5: Image upscaling results obtained with the dictionary learning technique of Couzinie-Devy et al. [2011]. Images produced by Florent Couzinie-Devy.

Figure 3 . 6 :

 36 Figure 3.6: Inverse halftoning results obtained byMairal et al. [2012]. The last row presents a result on a binary image publicly available on the Internet. No ground truth is available for this binary image from an old computer game. Despite the fact that the algorithm that has generated this image is unknown and probably different than the Floyd-Steinberg algorithm, the result is visually convincing. Best seen by zooming on a computer screen.

Figure 3 . 7 :

 37 Figure 3.7: Color video denoising result from Mairal et al. [2008e]. The third column shows the result when each frame is processed independently from the others. Last column shows the result of the video processing approach. Best seen by zooming on a computer screen. "Copyright c 2008 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved".

Figure 3 . 8 :

 38 Figure 3.8: Video inpainting result from Mairal et al. [2008e]. The third column shows the result when each frame is processed independently from the others. Last column shows the result of the video processing approach. Best seen by zooming on a computer screen. "Copyright c 2008 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved".

Figure 3 . 9 :

 39 Figure 3.9: Results of the face compression method of Bryt and Elad [2008] compared to two classical image compression algorithms. All images are compressed into 550 bytes representations. The raw output of the method, denoted by KSVD, is presented in the penultimate column. The last column KSVD+ is a post-processed result obtained by processing the output of KSVD with a deblocking algorithm that reduces visual artifacts. The numbers in parenthesis are the square root of the mean squared errors (RMSE), obtained by comparing the compressed images with the original one (first column). Images kindly provided to us by Ori Bryt and Michael Elad. Best seen by zooming on a computer screen.

Figure 3 . 10 :

 310 Figure 3.10: Sparsity vs. joint sparsity: Grey squares represents non-zeros values in vectors (left) or matrix (right).

Figure 4 . 1 :

 41 Figure 4.1: Illustration representing one layer of the hierarchical convolutional sparse coding model of Zeiler et al. [2011]. The vector x denotes the input feature maps, which are encoded into sparse codes A = [α k ] l k=1 with the convolutional sparse coding formulation. The output feature maps B = [β j ] l ′j=1 are obtained after a three-dimensional max pooling step, reducing the spatial resolution and the number of maps from p to p ′ .

  More specifically, let us first consider a set of training signals from k different classes, represented by the matrix D = [D 1 , . . . , D k ] in R m×p where the columns of each submatrix D j are signals from the j-th class. In order to classify a new test signal x in R m ,

Figure 4 . 2 :

 42 Figure 4.2: Results of pixelwise classification of test images from Mairal et al.[2008a]. Two dictionaries are learned beforehand, one from images containing a bicycle, and one from images corresponding to background only. For each of the four test images displayed here, we also display a confidence map for the presence of a bicycle. For each pixel, a confidence value is computed for a patch centered at the pixel, by comparing the residual of the reconstruction errors obtained with the two dictionaries; yellow bright pixel values represent high confidence, whereas dark red pixel values represent low confidence.

Figure 4 . 3 :

 43 Figure 4.3: Edge detection results obtained with the software package of Xiaofeng and Bo [2012] on the Berkeley segmentation dataset BSD500 [Martin et al., 2001]. Four pairs of images are presented. The left ones represent input images and the right ones confidence maps for the presence of edges. Dark values represent high confidence. Best seen by zooming on a computer screen.

Figure 4 . 4 :

 44 Figure 4.4: Visualization of three loss functions for a positive label y = 1. The plain red curve represents the square loss L(u) = (1/2)(y -u) 2 , the dotted blue one represents the logistic loss L(u) = log(1-e -yu ), and the green dotted one represents the hinge loss L(u) = max(0, 1 -yu).

Figure 4 . 5 :

 45 Figure 4.5: Two paradigms for dictionary learning are illustrated in the figure.An input signal x is sparsely encoded by using a dictionary D, producing sparse codes α. Then, a second model with parameters W makes a prediction y based on α. In (a), the dictionary D is learned unsupervised in a first layer-that is, without exploiting output information y. Then, the parameters W of the prediction model are learned afterwards in the second layer. In (b), we illustrate the concept of supervised dictionary learning, where D and W are jointly learned for the prediction task; this requires "backpropagating" information, represented by the red arrow, from the top layer to the bottom one.

Figure 4 . 6 :

 46 Figure 4.6:A few archetypes learned from 36 600 pictures corresponding to the request "Paris" downloaded from Flickr and sorted by "relevance". The figure is produced by[START_REF] Chen | Fast and robust archetypal analysis for representation learning[END_REF] with the software package SPAMS available at http://spams-devel.gforge.inria.fr/ The numbers in red represent the values of the non-zero coefficients β j [i] for each image. In (a), we obtain expected landmarks such as the Eiffel tower or Notre Dame; in (b), we display some unexpected archetypes corresponding to soccer, graffitis, food, flowers and social gatherings; in (c), we see some archetypes that do not seem to have some semantic meaning, but that capture some scene compositions or texture present in the dataset. Best seen by zooming on a computer screen.

Algorithm 4

 4 .1. Iterative hard thresholding.Require: Signal x in R m , dictionary D in R m×p , target sparsity k or penalty parameter λ, step size η, number of iterations T , initial solution α 0 in R p (with α 0 0 ≤ k if k is provided). 1: Initialize α ← α 0 ; 2: for t = 1, . . . , T do 3:perform one step of gradient descent:α ← α + ηD ⊤ (x -Dα); 4:choose threshold τ to be the k-th largest magnitude among the entries of α if k is provided, or τ = √ 2λ if λ is provided;

Figure 5 . 1 :

 51 Figure 5.1: Illustration of the basic majorization-minimization principle. The objective function f (in blue) needs to be minimized. At the current iteration, a majorizing surrogate g β (in red) is computed that is locally tight at the previous estimate α old . After minimizing the surrogate, we obtain a new estimate α new such that f (α new ) ≤ f (α old ).

1 :

 1 Initialize α ← α 0 ; 2: for t = 1, . . . , T do 3: perform one step of gradient descent: α ← α + ηD ⊤ (x -Dα); : for all j = 1, . . . , p, α[j] ← S λ (α[j]);

  p + obtained as η[j] = |α[j]|. The Lasso optimization problem then becomes: min α∈R p min η∈R p

  t. ψ(α i ) ≤ µ,(5.18) and X = [x 1 , . . . , x n ] in R m×n is the training set of signals. Even though these problems are non-convex, many algorithms have been developed

  , the algorithm selects one signal î and performs one gradient step D ← Dη t ∇ D L(x î, D) with a step size η t . It can indeed be shown that the function D → L(x î, D) is differentiable according to Danskin's theorem [seeBertsekas, 1999, Proposition B.25] and that its gradient admits a closed form∇ D L(x î, D) = -(x î -Dα î)α ⊤î , where α î is defined in (5.21).

  ).Require: SignalsX = [x 1 , . . . , x n ] in R m×n , initial dictionary D 0 in C, penalty parameter λ, number of iterations T . 1: Initialize D ← D 0 ; 2: for t = 1, . . . , T do 3:select one signal x î at random; gradient descent step:D ← Π C Dη t (Dα îx i )α ⊤ î ;6: end for 7: return the sparse decomposition α in R p .

Algorithm 9

 9 Method of optimal directions for ψ = ℓ 0 .Require: SignalsX = [x 1 , . . . , x n ] in R m×n , initial dictionary D 0 in C, regularization parameter µ, number of iterations T . 1: Initialize D ← D 0 ; 2: for t = 1, . . . , T do 3:compute the sparse codes, e.g., with OMP: C [XA ⊤ (AA ⊤ ) -1 ];(5.22)8: end for 9: return the dictionary D in C.

Algorithm 10

 10 Alternate minimization for ψ = ℓ 1 . Require: Signals X = [x 1 , . . . , x n ] in R m×n , initial dictionary D 0 in C, regularization parameter λ or µ, number of iterations T . 1: Initialize D ← D 0 ; 2: for t = 1, . . . , T do 3: compute the sparse codes, e.g., with the homotopy method: 4: for i = 1, . . . , n do 5:

  + d j α jdα j 2 F , where d j on the right side is the value of the variable d j before the update, and α j in R 1×n is the j-th row of the matrix A. After expanding the Frobenius norm and removing the constant term, we obtaind j ← arg min d∈R m , d 2 ≤1 -d ⊤ (X -DA + d j α j )α j⊤ + d∈R m , d 2 ≤1 -d ⊤ (b j -Dc j + d j C[j, j]) B = [b 1 , . . . , b p ] in R m×p and C = [c 1 , . . . , c p ]in R p×p are defined in Algorithm 11. Finally, we see from the last equation that updating d j amounts to performing one orthogonal projection of the vector (1/C[j, j])(b j -Dc j ) + d j onto the unit Euclidean ball, leading to the update (5.24).Algorithm 11 Dictionary update with block coordinate descent.Require: D 0 ∈ C (input dictionary); X ∈ R m×n (dataset); A ∈ R p×n (sparse codes);1: Initialization: D ← D 0 ; B ← XA ⊤ ; C ← AA ⊤ ; 2: repeat 3: for j = 1, . . . , p do 4:update the j-th column to optimize for (5.23):d j ← 1 C[j, j](b j -Dc j ) + d j ,d j ← 1 max( d j 2 , 1) d j .

  ) with ψ = ℓ 1 .Require: SignalsX = [x 1 , . . . , x n ] in R m×n , initial dictionary D 0 in C, initial sparse coefficients A 0 in R p×n (possibly equal to zero), regularization parameter λ, number of iterations T . 1: Initialize D ← D 0 ; A ← A 0 ; 2: for t = 1, . . . , T do 3:update the sparse codes, one row of A at a time; perform at least one pass (possibly several) of the following iterations: 4: for j = 1, . . . , p do 5: α j ← S λ α j +

  Ω {i ∈ {1, . . . , n} : α i [j] = 0}; 10: update d j and the non-zero coefficients of α j : (d j , α j Ω ) ∈ arg min d∈R m , d 2 ≤1 β∈R |Ω|

  ⊤ DC)trace(D ⊤ B) ; (5.28) 7: end for 8: return learned dictionary D in C.

  . . , α l ]. Solving (2.11) is easy with any standard method adapted to large-scale ℓ 1 -regularized convex problems [seeBach et al., 2012a, for a review]. For instance,[START_REF] Rigamonti | Are sparse representations really relevant for image classification[END_REF] uses a proximal gradient method, whereasKavukcuoglu et al. [2010a] uses a coordinate descent scheme. To conduct the experiment of this paragraph, we have implemented the ISTA and FISTA algorithms of[START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] for solving (2.11), and we have made it available in the SPAMS toolbox. 6

  -global/ℓ 0 37.21 33.37 31.41 30.01 28.95 25.65 22.44 ℓ 1 -global/ℓ 0 37.22 33.50 31.56 30.17 29.13 25.77 22.53 ℓ 0 -adapt/ℓ 0 37.49 33.75 31.70 30.40 29.33 26.04 22.64 ℓ 1 -adapt/ℓ 0 37.60 33.90 31.90 30.51 29.43 26.20 22.72 ℓ 1 -adapt/ℓ 1 36.58 32.85 30.77 29.29 28.14 24.47 21.83

					σ
		5	10	15	20	25	50	100
	DCT	37.30 33.38 31.24 29.75 28.66 25.24 22.00
	ℓ 0 GSM	37.05 33.34 31.31 29.91 28.84 25.66 22.80
	K-SVD	37.42 33.62 31.58 30.18 29.10 25.61 22.10
	BM3D	37.62 34.00 32.05 30.73 29.72 26.38 23.25
	LSSC	37.67 34.06 32.12 30.78 29.74 26.57 23.39
	Plow	37.38 32.98 31.38 30.13 29.30 26.38 23.24
	EPLL	37.36 33.64 31.67 30.32 29.29 26.12 23.03
	CSR	37.61 34.00 32.05 30.72 29.70 26.53 23.45
	BM3D-PCA 37.79 34.20 32.27 30.94 29.92 26.75 23.16

Table 3 . 1 :

 31 Denoising performance in PSNR for various methods on the 12 standard images of Figure

  does not change, or cost O(m) otherwise. Finally, the periteration complexity of this strategy is always O(m); 2. assuming that the matrix D ⊤ D is pre-computed, it is possible to obtain an even better per-iteration complexity by updating an auxiliary variable z instead of r such that we always have z = D

Given r, updating α[] only cost O(m) operations because of the inner-product d ⊤  r. Then, updating r either has zero cost if α[] ⊤ (x -Dα). Given z, updating α[j] only cost O(1) operations. Updating z has either zero cost if α[j] does not change, or cost O(p) operations. As a result, the per-iteration complexity is at most O(p), but more interestingly, each time a variable α[j] is equal to zero before its update and remains equal to zero after the update (a typical scenario for very sparse solutions), the cost per iteration is O(1).

  .15)where the functions ϕ are concave, non-decreasing and differentiable on R + . A natural way of dealing with such penalties when they appear in the objective function is to use an approach called DC-programming, where DC stands for "difference of convex" [seeGasso et al., 2009, for a review]. Even though the resulting optimization problems are nonconvex and impossible to solve exactly in a reasonable amount of time, DC-programming is a simple technique to obtain a stationary point by using the majorization-minimization principle already illustrated in Figure5.1. It consists of exploiting a majorizing surrogate of ψ obtained by linearizing the function ϕ in (5.15), which is upper-bounded by its linear approximations because of concavity, as shown in Figure5.2. Illustration of the DC-programming approach. The function ϕ (in blue) is upper-bounded by its linear approximation (in red), which is tight at the point β. The resulting majorizing surrogate is a weighted ℓ1-norm plus a constant C β .

	g β (α) = |α| |β|+ε + C β
	|β|
	ϕ(α) = log(|α| + ε)
	Figure 5.2:

available here http://spams-devel.gforge.inria.fr/.

Note that the Gaussian noise assumption is not necessary to justify the ordinary least square formulation. It is only sufficient to interpret it as maximum likelihood estimation. In fact, as long as the conditional expectation E[y|x] is linear, the ordinary least square estimator is statistically consistent under mild assumptions.

Note that the soft-thresholding operator appears in fact earlier in the statistics literature [seeEfron and Morris, 1971, Bickel, 1984], but it was used there for a

A collection of vectors z1, z2, . . . , z l are called generators for a cone K when K consists of all positive combinations of the vectors zi. In other words, K = { l i=1 αizi : αi ≥ 0}. In that case, we use the notation K = cone(z1, . . . , z l ).

We use "a.s." to denote almost sure convergence.

Note that recent research suggests that this fast rate of σ 2 k log p m cannot be achieved by polynomial-time algorithms[START_REF] Zhang | Lower bounds on the performance of polynomial-time algorithms for sparse linear regression[END_REF].

Neuroscientists have identified several areas in the human visual cortex. Area V1 is considered to be at the earliest stage of visual processing.

available here: http://r0k.us/graphics/kodak/.

The entropy is defined as H(x) p(x) log p(x)dx [see[START_REF] Cover | Elements of information theory[END_REF]. In general, H(x) is not computable exactly since the probability density p is unknown.

Note that choosing the Laplace distribution for the prior does not mean implicitly assuming that the true distribution of the latent variable α is Laplace. In fact,[START_REF] Gribonval | Compressible distributions for high-dimensional statistics[END_REF] have shown that if this was the case, the choice (2.7) would be inappropriate in the context of maximum a posteriori estimation.

It seems that a common misconception in the computer vision and image processing literature is to call "Bayesian" any probabilistic model that uses Bayes' rule such as (2.8). Such a terminology for maximum a posteriori (MAP) estimation is slightly misleading regarding the hundred-year-old debate among frequentists and Bayesian statisticians. In general, the latter do not assume that there exists a "true" parameter that can be obtained with MAP estimation; instead, they treat all parameters as random variables and model their uncertainty [see[START_REF] Bayarri | The interplay of bayesian and frequentist analysis[END_REF].

http://spams-devel.gforge.inria.fr/.

The formulations presented in this section can be also modified to handle other types of noise, such as Poisson noise [seeGiryes and Elad,

2014].

For simplicity, we always assume in this monograph that the patches are square, but all approaches can be easily extended to deal with other patch shapes.

For simplicity, we choose here to neglect boundary effects for pixels that are close to the image border, which technically belong to fewer patches.

Denoting by MSE the mean squared-error for images whose intensities are between 0 and 255, the PSNR is defined as PSNR = 10 log 10 (255 2 /MSE) and is measured in dB. A gain of 1 dB reduces the MSE by approximately

20%. 6 The dataset is available here: http://r0k.us/graphics/kodak/.

When the solution is not unique, the elastic-net regularization[START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] can be used instead of the ℓ1-norm. It consists of replacing the penalty α 1 by α 1 + (µ/2) α 2 2 . The resulting penalty still enjoy sparsity-inducing properties, but is strongly convex and thus yields a unique solution.

Note that the tf-idf renormalization common in image retrieval[Sivic and Zisserman,

2003] can be achieved by reweighting the components of γ(X): let γ i denote the i-th component of γ(X), i.e., the frequency of the corresponding bin in X, replace γ i by γ i log(N/nγ i ), where N is the total number of features in the database being queried.

One should keep in mind that being invariant to within-image transformations does not imply being invariant to, say, rotations in depth since (at least) some features will become occluded, or will be revealed, as the observed object rotates relative to the camera.

Note that the definition of "manifold" from[START_REF] Yu | Nonlinear learning using local coordinate coding[END_REF] slightly differs from the traditional one. More precisely, they define a "manifold" M as a subset of R m such there exists a constant C and "dimension" p such that for all x in M, there

Note that the fast Fourier transform (FFT) could be used to reduce this theoretical complexity. Unfortunately, the FFT induces a computational overhead that make it inefficient in practice when convolving a large image with a small filter, which is the case here.

The software package is available here: http://homes.cs.washington.edu/~xren.

see http://www.netlib.org/blas/.

http://spams-devel.gforge.inria.fr/.

s , where s d ⊤  d d ⊤  (D ⊤ Γ D Γ ) -1 d  is called the Schur complement, and z (D ⊤ Γ D Γ ) -1 d . Thus, computing (D ⊤ Γ ′ D Γ ′ ) -1 given (D ⊤ Γ D Γ ) -1

Acknowledgements

This monograph was partially supported by the European Research Council (SIERRA and VideoWorld projects), the project Gargantua funded by the program Mastodons-CNRS, the Microsoft Research-Inria joint centre, and a French grant from the Agence Nationale de la Recherche (MACARON project, ANR-14-CE23-0003-01).

The authors would like to thank Ori Bryt and Michael Elad for providing the output of their face compression algorithm (see Figure 3.9). Julien Mairal would also like to thank Zaid Harchaoui, Florent Perronnin and Adrien Gaidon for useful discussions leading to improvements of this monograph.

from the previous iteration can be done by a few matrix-vector multiplications only. In the same vein, obtaining

is also a matter of manipulating a few linear algebra relations.

We also remark that the implementation of the homotopy can benefit from the pre-computation of the matrix D ⊤ D in the exact same way as for orthogonal matching pursuit, such that both implementations have in fact the same complexity (up to a constant factor that is slightly larger for the homotopy method).

Iterative reweighted-ℓ 1 methods

We have mentioned in the introduction of this monograph that nonconvex sparsity-inducing penalties are often used in sparse estimation.