Task-Driven Dictionary Learning

Introduction

The linear decomposition of data using a few elements from a learned dictionary instead of a predefined one—based on wavelets for example—has recently led to state-of-the-art results in numerous low-level signal processing tasks such as image denoising , audio processing , as well as classification tasks . Unlike decompositions based on principal component analysis (PCA) and its variants, these sparse models do not impose that the dictionary elements be orthogonal, allowing more flexibility to adapt the representation to the data.

This classical data-driven approach to dictionary learning is well adapted to reconstruction tasks, such as restoring a noisy signal. These dictionaries, which are good at reconstructing clean signals, but bad at reconstructing noise, have indeed led to state-of-the-art denoising algorithms . Unsupervised dictionary learning has also been used for other purposes than pure signal reconstruction, such as classification , but recent works have shown that better results can be obtained when the dictionary is tuned to the specific task (and not just data) it is intended for. Duarte-Carvajalino and Sapiro have for instance proposed to learn dictionaries for compressed sensing, and in dictionaries are learned for signal classification. In this paper, we will refer to this type of approach as task-driven dictionary learning.

Whereas purely data-driven dictionary learning has been shown to be equivalent to a large-scale matrix factorization problem that can be effectively addressed with several methods , its task-driven counterpart has proven to be much more difficult to optimize. Presenting a general efficient framework for various task-driven dictionary learning problems is the main topic of this paper. Even though it is different from existing machine learning approaches, it shares similarities with many of them.

For instance, Blei et al. have proposed to learn a latent topic model intended for document classification. In a different context, Argyriou et al. introduced a convex formulation for multi-task classification problems where an orthogonal linear transform of input features is jointly learned with a classifier. Learning compact features has also been addressed in the literature of neural networks, with restricted Boltzmann machines (RBM’s) and convolutional neural networks for example (see and references therein). Interestingly, the question of learning the data representation in an unsupervised or supervised way has also been investigated for these approaches. For instance, a supervised topic model is proposed in and tuning latent data representations for minimizing a cost function is often achieved with backpropagation in neural networks .

This paper makes three main contributions:

It introduces a supervised formulation for learning dictionaries adapted to various tasks instead of dictionaries only adapted to data reconstruction.

It shows that the resulting optimization problem is smooth under mild assumptions, and empirically that stochastic gradient descent addresses it efficiently.

It shows that the proposed formulation is well adapted to semi-supervised learning, can exploit unlabeled data when they admit sparse representations, and leads to state-of-the-art results for various machine learning and signal processing problems.

2 Notation

The rest of this paper is organized as follows: Section 2 presents the data-driven dictionary learning framework. Section 3 is devoted to our new task-driven framework, and Section 4 to efficient algorithms to addressing the corresponding optimization problems. Section 5 presents several dictionary learning experiments for signal classification, signal regression, and compressed sensing.

Data-Driven Dictionary Learning

As pointed out by Bottou and Bousquet , one is usually not interested in a perfect minimization of the empirical cost $g_{n}({\mathbf{D}})$, but instead in the minimization with respect to ${\mathbf{D}}$ of the expected cost

where the expectation is taken relative to the (unknown) probability distribution $p({\mathbf{x}})$ of the data, and is supposed to be finite.We use “a.s.” (almost surely) to denote convergence with probability one. In practice, dictionary learning problems often involve a large amount of data. For instance when the vectors ${\mathbf{x}}$ represent image patches, $n$ can be up to several millions in a single image. In this context, online learning techniques have shown to be very efficient for obtaining a stationary point of this optimization problem . In this paper, we propose to minimize an expected cost corresponding to a supervised dictionary learning formulation, which we now present.

Proposed Formulation

We introduce in this section a general framework for learning dictionaries adapted to specific supervised tasks, e.g., classification, as opposed to the unsupervised formulation of the previous section, and present different extensions along with possible applications.

Obtaining a good performance in classification tasks is often related to the problem of finding a good data representation. Sparse decompositions obtained with data-driven learned dictionaries have been used for that purpose in and , showing promising results for audio data and natural images. We present in this section a formulation for learning a dictionary in a supervised way for regression or classification tasks.

where ${\mathbf{W}}$ are model parameters which we want to learn, ${\mathcal{W}}$ is a convex set, $\nu$ is a regularization parameter, and $f$ is a convex function defined as

where ${\mathcal{D}}$ is a set of constraints defined in Eq. (2), and $f$ has the form

The main difficulty of this optimization problem comes from the non-differentiability of ${\boldsymbol{\alpha}}^{\star}$, which is the solution of a nonsmooth optimization problem (4). Bradley and Bagnell have tackled this difficulty by introducing a smooth approximation of the sparse regularization which leads to smooth solutions, allowing the use of implicit differentiation to compute the gradient of the cost function they have introduced. This approximation encourages some coefficients in ${\boldsymbol{\alpha}}^{\star}$ to be small, and does not produce true zeros. It can be used when “true” sparsity is not required. In a different formulation, Mairal et al. have used nonsmooth sparse regularization, but used heuristics to tackle the optimization problem. We show in Section 4 that better optimization tools than these heuristics can be used, while keeping a nonsmooth regularization for computing ${\boldsymbol{\alpha}}^{\star}$.

A difference between supervised and unsupervised dictionary learning is that overcompleteness—that is, the dictionaries have more elements than the signal dimension, has not empirically proven to be necessary. It is indeed often advocated for image processing applications that having $p>m$ provides better reconstruction results , but for discriminative tasks, perfect reconstruction is not always required as long as discriminative features are captured by the sparse coding procedure.

Before presenting extensions and applications of the formulation we have introduced, let us first discuss the assumptions under which our analysis holds.

The data $({\mathbf{y}},{\mathbf{x}})$ admits a probability density $p$ with a compact support $K_{{\mathcal{Y}}}\times K_{{\mathcal{X}}}\subseteq{\mathcal{Y}}\times{\mathcal{X}}$. This is a reasonable assumption in audio, image, and video processing applications, where it is imposed by the data acquisition process, where values returned by sensors are bounded. To simplify the notation we assume from now on that ${\mathcal{X}}$ and ${\mathcal{Y}}$ are compact.Even though images are acquired in practice after a quantization process, it is a common assumption in image processing to consider pixel values in a continuous space.

Assumptions (B) and (C) allow us to use several loss functions such as the square, logistic, or softmax losses.

2 Extensions

We now present two extensions of the previous formulations. The first one includes a linear transform of the input data, and the second one exploits unlabeled data in a semi-supervised setting.

In this section, we add to our basic formulation a linear transform of the input features, represented by a matrix ${\mathbf{Z}}$. Our motivation for this is twofold: It can be appealing to reduce the dimension of the feature space via such a linear transform, and/or it can make the model richer by increasing the numbers of free parameters. The resulting formulation is the following:

where $\nu_{1}$ and $\nu_{2}$ are two regularization parameters, ${\mathcal{Z}}$ is a convex set and

It is worth noticing that the formulations of Eq. (7) and Eq. (9) can also be extended to the case of a cost function depending on several dictionaries involving several sparse coding problems, such as the one used in for signal classification. Such a formulation is not developed here for simplicity reasons, but algorithms to address it can easily be derived from this paper.

2.2 Semi-supervised Learning

As shown in , sparse coding techniques can be effective for learning good features from unlabeled data. The extension of our task-driven formulation to the semi-supervised learning setting is natural and takes the form

3 Applications

For illustration purposes, we present a few applications of our task-driven dictionary learning formulations. Our approach is of course not limited to these examples.

In this setting, ${\mathcal{Y}}$ is a subset of a $q$-dimensional real vector space, and the task is to predict variables ${\mathbf{y}}$ in ${\mathcal{Y}}$ from the observation of vectors ${\mathbf{x}}$ in ${\mathcal{X}}$. A typical application is for instance the restoration of clean signals ${\mathbf{y}}$ from observed corrupted signals ${\mathbf{x}}$. Classical signal restoration techniques often focus on removing additive noise or solving inverse linear problems . When the corruption results from an unknown nonlinear transformation, we formulate the restoration task as a general regression problem. This is the case for example in the experiment presented in Section 5.3.

We define the task-driven dictionary learning formulation for regression as follows:

At test time, when a new signal ${\mathbf{x}}$ is observed, the estimate of the corresponding variable ${\mathbf{y}}$ provided by this model is ${\mathbf{W}}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$ (plus possibly an intercept which we have omitted here for simplicity reasons). Note that we here propose to use the square loss for estimating the difference between ${\mathbf{y}}$ and its estimate ${\mathbf{W}}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$, but any other twice differentiable loss can be used.

3.2 Binary Classification

Once ${\mathbf{D}}$ and ${\mathbf{w}}$ have been learned, a new signal ${\mathbf{x}}$ is classified according to the sign of ${\mathbf{w}}^{\top}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$. For simplicity reasons, we have omitted the intercept in the linear model, but it can easily be included in the formulation. Note that instead of the logistic regression loss, any other twice differentiable loss can be used.

As suggested in , it is possible to extend this approach with a bilinear model by learning a matrix ${\mathbf{W}}$ so that a new vector ${\mathbf{x}}$ is classified according to the sign of ${\mathbf{x}}^{\top}{\mathbf{W}}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$. In this setting, our formulation becomes

This bilinear model requires learning $pm$ parameters as opposed to the $p$ parameters of the linear one. It is therefore richer and can sometimes offer a better classification performance when the linear model is not rich enough to explain the data, but it might be more subject to overfitting.

Note that we have naturally presented the binary classification task using the logistic regression loss, but as we have experimentally observed, the square loss is also an appropriate choice in many situations.

3.3 Multi-class Classification

When ${\mathcal{Y}}$ is a finite set of labels in $\{1,\ldots,q\}$ with $q>2$, extending the previous formulation to the multi-class setting can be done in several ways, which we briefly describe here. The simplest possibility is to use a set of binary classifiers presented in Section 3.3.2 in a “one-vs-all” or “one-vs-one” scheme. Another possibility is to use a multi-class cost function such as the soft-max function, to find linear predictors ${\mathbf{w}}_{k}$, $k$ in $\{1,\ldots,q\}$ such that for a vector ${\mathbf{x}}$ in ${\mathcal{X}}$, the quantities ${\mathbf{w}}_{y}^{\top}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$ are encouraged to be greater than ${\mathbf{w}}_{k}^{\top}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$ for all $k\neq y$. Another possibility is to turn the multi-class classification problem into a regression one and consider that ${\mathcal{Y}}$ is a set of $q$ binary vectors of dimension $q$ such that the $k-$th vector has $1$ on its $k$-th coordinate, and elsewhere. This allows using the regression formulation of Section 3.3.1 to solve the classification problem.

We remark that for classification tasks, scalability issues should be considered when choosing between a one-vs-all scheme (learning independent dictionaries for every class) and using a multi-class loss function (learning a single dictionary shared between all classes). The one-vs-all scheme requires keeping into memory $qpm$ parameters, where $q$ is the number of classes, which is feasible when $q$ is reasonably small. For classifications problems with many classes (for instance $q\geq 1\,000$), using a single (larger) dictionary and a multi-class loss function is more appropriate, and would in addition allow feature sharing between the classes.

3.4 Compressed sensing

In a nutshell, the recovery of ${\mathbf{x}}$ has been proven to be possible when ${\mathbf{x}}$ admits a sparse representation on a dictionary ${\mathbf{D}}$, and the sensing matrix ${\mathbf{Z}}$ is incoherent with ${\mathbf{D}}$, meaning that the rows of ${\mathbf{Z}}$ are sufficiently uncorrelated with the columns of ${\mathbf{D}}$ (see for more details).The assumption of “incoherence” between ${\mathbf{D}}$ and ${\mathbf{Z}}$ can be replaced with a different but related hypothesis called restricted isometry property. Again the reader should refer to for more details. To ensure that this condition is satisfied, ${\mathbf{Z}}$ is often chosen as a random matrix, which is incoherent with any dictionary with high probability.

The choice of a random matrix is appealing for many reasons. In addition to the fact that it provides theoretical guarantees of incoherence, it is well suited to the case where $m$ is large, making it impossible to store a deterministic matrix ${\mathbf{Z}}$ into memory, whereas it is sufficient to store the seed of a random process to generate a random matrix. On the other hand, large signals can often be cut into smaller parts that still admit sparse decompositions, e.g., image patches, which can be treated independently with a deterministic smaller matrix ${\mathbf{Z}}$. When this is the case or when $m$ has a reasonable size, the question of whether to use a deterministic matrix ${\mathbf{Z}}$ or a random one arises, and it has been empirically observed that learned matrices ${\mathbf{Z}}$ can outperform random projections: For example, it is shown in that classical dimensionality reduction techniques such as principal component analysis (PCA) or independent component analysis (ICA) could do better than random projections in noisy settings, and in that jointly learning sensing matrices and dictionaries can do even better in certain cases. A Bayesian framework for learning sensing matrices in compressed sensing applications is also proposed in .

Following the latter authors, we study the case where ${\mathbf{Z}}$ is not random but learned at the same time as the dictionary, and introduce a formulation which falls into out task-driven dictionary learning framework:

where we learn ${\mathbf{D}}$, ${\mathbf{W}}$ and ${\mathbf{Z}}$ so that the variable ${\mathbf{y}}$ should be well reconstructed when encoding the “sensed” signal ${\mathbf{Z}}{\mathbf{x}}$ with a dictionary ${\mathbf{D}}$. In a noiseless setting, ${\mathbf{y}}$ is naturally set to the same value as ${\mathbf{x}}$. In a noisy setting, it can be a corrupted version of ${\mathbf{x}}$.

After having presented our general task-driven dictionary learning formulation, we present next a strategy to address the corresponding nonconvex optimization problem.

Optimization

We first show that the cost function $f$ of our basic formulation (7) is differentiable and compute its gradient. Then, we refine the analysis for the different variations presented in the previous section, and describe an efficient online learning algorithm to address them.

We analyze the differentiability of $f$ as defined in Eq. (7) with respect to its two arguments ${\mathbf{D}}$ and ${\mathbf{W}}$. We consider here the case where ${\mathcal{Y}}$ is a compact subset of a finite dimensional real vector space, but all proofs and formulas are similar when ${\mathcal{Y}}$ is a finite set of labels. The purpose of this section is to show that even though the sparse coefficients ${\boldsymbol{\alpha}}^{\star}$ are obtained by solving a non-differentiable optimization problem, $f$ is differentiable on ${\mathcal{W}}\times{\mathcal{D}}$, and one can compute its gradient.

The main argument in the proof of Propositions 1 and 2 below is that, although the function ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$ is not differentiable, it is uniformly Lipschitz continuous, and differentiable almost everywhere. The only points where ${\boldsymbol{\alpha}}^{\star}$ is not differentiable are points where the set of nonzero coefficients of ${\boldsymbol{\alpha}}^{\star}$ change (we always denote this set by $\Lambda$ in this paper). Considering optimality conditions of the elastic-net formulation of Eq. (1), these points are easy to characterize. The details of the proof have been relegated to the Appendix (Lemma 1 and Proposition 3) for readability purposes. With these results in hand, we then show that $f$ admits a first-order Taylor expansion meaning that it is differentiable, the sets where ${\boldsymbol{\alpha}}^{\star}$ is not differentiable being negligible in the expectation from the definition of $f$ in Eq. (8). We can now state our main result:

Assume $\lambda_{2}>0$, (A), (B) and (C). Then, the function $f$ defined in Eq. (7) is differentiable, and

where $\Lambda$ denotes the indices of the nonzero coefficients of ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$.

The proof of this proposition is given in Appendix. We have shown that the function defined in Eq. (7) is smooth, and computed its gradients. The same can be done for the more general formulation of Eq. (10):

Assume $\lambda_{2}>0$, (A), (B) and (C). Then, the function $f$ defined in Eq. (10) is differentiable. The gradients of $f$ are

where ${\boldsymbol{\alpha}}^{\star}$ is short for ${\boldsymbol{\alpha}}^{\star}({\mathbf{Z}}{\mathbf{x}},{\mathbf{D}})$, and ${\boldsymbol{\beta}}^{\star}$ is defined in Eq. (17).

The proof is similar to the one of Proposition 1 in Appendix, and uses similar arguments.

2 Algorithm

Stochastic gradient descent algorithms are typically designed to minimize functions whose gradients have the form of an expectation as in Eq. (16). They have been shown to converge to stationary points of (possibly nonconvex) optimization problems under a few assumptions that are a bit stricter than the ones satisfied in this paper (see and references therein).As often done in machine learning, we use stochastic gradient descent in a setting where it is not guaranteed to converge in theory, but is has proven to behave well in practice, as shown in our experiments. The convergence proof of Bottou for non-convex problems indeed assumes three times differentiable cost functions. As noted in , these algorithms are generally well suited to unsupervised dictionary learning when their learning rate is well tuned.

The method we propose here is a projected first-order stochastic gradient algorithm (see ), and it is given in Algorithm 1. It sequentially draws i.i.d samples $({\mathbf{y}}_{t},{\mathbf{x}}_{t})$ from the probability distribution $p({\mathbf{y}},{\mathbf{x}})$. Obtaining such i.i.d. samples may be difficult since the density $p({\mathbf{y}},{\mathbf{x}})$ is unknown. At first approximation, the vectors $({\mathbf{y}}_{t},{\mathbf{x}}_{t})$ are obtained in practice by cycling over a randomly permuted training set, which is often done in similar machine learning settings .

At each iteration, the sparse code ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}}_{t},{\mathbf{D}})$ is computed by solving the elastic-net formulation of . We have chosen to use the LARS algorithm, a homotopy method , which was originally developed to solve the Lasso formulation—that is, $\lambda_{2}=0$, but which can be modified to solve the elastic-net problem. Interestingly, it admits an efficient implementation that provides a Cholesky decomposition of the matrix $({\mathbf{D}}_{\Lambda}^{\top}{\mathbf{D}}_{\Lambda}+\lambda_{2}{\mathbf{I}})^{-1}$ (see ) as well as the solution ${\boldsymbol{\alpha}}^{\star}$. In this setting, ${\boldsymbol{\beta}}^{\star}$ can be obtained without having to solve from scratch a new linear system.

The learning rate $\rho_{t}$ is chosen according to a heuristic rule. Several strategies have been presented in the literature (see and references therein). A classical setting uses a learning rate of the form $\rho/t$, where $\rho$ is a constant.A $1/t$-asymptotic learning rate is usually used for proving the convergence of stochastic gradient descent algorithms . However, such a learning rate is known to decrease too quickly in many practical cases, and one sometimes prefers a learning rate of the form $\rho/(t+t_{0})$, which requires tuning two parameters. In this paper, we have chosen a learning rate of the form $\min(\rho,\rho t_{0}/t)$—that is, a constant learning rate $\rho$ during $t_{0}$ iterations, and a $1/t$ annealing strategy when $t>t_{0}$, a strategy used by for instance. Finding good parameters $\rho$ and $t_{0}$ also requires in practice a good heuristic. The one we have used successfully in all our experiments is $t_{0}=T/10$, where $T$ is the total number of iterations. Then, we try several values of $\rho$ during a few hundreds of iterations and keep the one that gives the lowest error on a small validation set.

In practice, one can also improve the convergence speed of our algorithm with a mini-batch strategy—that is, by drawing $\eta>1$ samples at each iteration instead of a single one. This is a classical heuristic in stochastic gradient descent algorithms and, in our case, this is further motivated by the fact that solving $\eta$ elastic-net problems with the same dictionary ${\mathbf{D}}$ can be accelerated by the precomputation of the matrix ${\mathbf{D}}^{\top}{\mathbf{D}}$ when $\eta$ is large enough. Such a strategy is also used in for the classical data-driven dictionary learning approach. In practice, the value $\eta=200$ has given good results in all our experiments (a value found to be good for the unsupervised setting as well).

As many algorithms tackling non-convex optimization problems, our method for learning supervised dictionaries can lead to poor results if is not well initialized. The classical unsupervised approach of dictionary learning presented in Eq. (3) has been found empirically to be better behaved than the supervised one, and easy to initialize . We therefore have chosen to initialize our dictionary ${\mathbf{D}}$ by addressing the unsupervised formulation of Eq. (3) using the SPAMS toolbox .http://www.di.ens.fr/willow/SPAMS/ With this initial dictionary ${\mathbf{D}}$ in hand, we optimize with respect to ${\mathbf{W}}$ the cost function of Eq (5), which is convex. This procedure gives us a pair $({\mathbf{D}},{\mathbf{W}})$ of parameters which are used to initialize Algorithm 1.

3 Extensions

We here present the slight modifications to Algorithm 1 necessary to address the two extensions discussed in Section 3.2.

The last step of Algorithm 1 updates the parameters ${\mathbf{D}}$ and ${\mathbf{W}}$ according to the gradients presented in Eq. (18). Modifying the algorithm to address the formulation of Section 3.2.1 also requires updating the parameters ${\mathbf{Z}}$ according to the gradient from Proposition 2:

where $\Pi_{{\mathcal{Z}}}$ denotes the orthogonal projection on the set ${\mathcal{Z}}$.

The extension to the semi-supervised formulation of Section 3.2.2 assumes that one can draw samples from the marginal distribution $p({\mathbf{x}})$. This is done in practice by cycling over a randomly permuted set of unlabeled vectors. Extending Algorithm 1 to this setting requires the following modifications: At every iteration, we draw one pair $({\mathbf{y}}_{t},{\mathbf{x}}_{t})$ from $p({\mathbf{y}},{\mathbf{x}})$ and one sample ${\mathbf{x}}^{\prime}_{t}$ from $p({\mathbf{x}})$. We proceed exactly as in Algorithm 1, except that we also compute ${\boldsymbol{\alpha}}^{\star\prime}\stackrel{{\scriptstyle\vartriangle}}{{=}}{\boldsymbol{\alpha}}^{\star}({\mathbf{x}}^{\prime}_{t},{\mathbf{D}})$, and replace the update of the dictionary ${\mathbf{D}}$ by

Experimental Validation

Before presenting our experiments, we briefly discuss the question of choosing the parameters in our formulation.

Performing cross-validation on the parameters $\lambda_{1}$, $\lambda_{2}$ (elastic-net parameters), $\nu$ (regularization parameter) and $p$ (size of the dictionary) would of course be cumbersome, and we use a few simple heuristics to either reduce the search space for these parameters or fix arbitrarily some of them. We have proceeded in the following way:

Since we want to exploit sparsity, we often set $\lambda_{2}$ to , even though $\lambda_{2}>0$ is necessary in our analysis for proving the differentiability of our cost function. This has proven to give satisfactory results in most experiments, except for the experiment of Section 5.5, where choosing a small positive value for $\lambda_{2}$ was necessary for our algorithm to converge.

When there is a lot of training data, which is often the case for natural image patches, the regularization with $\nu$ becomes unnecessary and this parameter can arbitrarily set to a small value, e.g., $\nu=10^{-9}$ for normalized input data. When there are not many training points, this parameter is set up by cross-validation.

We have also observed that a larger dictionary usually means a better performance, but a higher computational cost. Setting the size of the dictionary is therefore often a trade-off between results quality and efficiency. In our experiments, we often try the values $p$ in $\{50,100,200,400\}$.

We show in this section several applications of our method to real problems, starting with handwritten digits classification, then moving to the restoration of images damaged by an unknown nonlinear transformation, digital art authentification, and compressed sensing.

2 Handwritten Digits Classification

We consider here a classification task using the MNIST and USPS handwritten datasets. MNIST contains $70\,000$ $28\times 28$ images, $60\,000$ for training, $10\,000$ for testing, whereas USPS has $7\,291$ training images and $2\,007$ test images of size $16\times 16$.

Most effective digit recognition techniques use features with shift invariance properties . Since our formulation is less sophisticated than for instance the convolutional network architecture of and does not enjoy such properties, we have artificially augmented the size of the training set by considering versions of the digits that are shifted by one pixel in every direction. This is of course not an optimal way of introducing shift invariance in our framework, but it is fairly simple.

After choosing the parameters using the validation set, we retrain our model on the full training set. Each experiment is performed with $40\,000$ iterations of our algorithm with a mini-batch of size $200$. We report the performance on the test set achieved for different dictionary sizes, with $p$ in $\{50,100,200,300\}$ for the two datasets, and observe that learning ${\mathbf{D}}$ in a supervised way significantly improves the performance of the classification. Moreover our method achieves state-of-the-art results on MNIST with a $0.54\%$ error rate, which is similar to the $0.60\%$ error rate of .It is also shown in that better results can be achieved by considering deformations of the training set. Our $2.84\%$ error rate on USPS is slightly behind the $2.4\%$ error rate of .

We remark that a digit recognition task was also carried out in , where a similar performance is reported.The error rates in are slightly higher but the dataset used in their paper is not augmented with shifted versions of the digits. Our conclusions about the advantages of supervised versus unsupervised dictionary learning are consistent with , but our approach has two main advantages. First it is much easier to use since it does not requires complicated heuristic procedures to select the parameters, and second it applies to a wider spectrum of applications such as to regression tasks.

Our second experiment follows , where only a few samples are labelled. We use the semi-supervised formulation of Section 3.2.2 which exploits unlabeled data. Unlike the first experiment where the parameters are chosen using a validation set, and following , we make a few arbitrary choices. Indeed, we use $p=300$, $\lambda_{1}=0.075$, and $\nu=10^{-5}$, which were the parameters chosen in the previous experiment. As in the previous experiment, we have observed that these parameters lead to sparse vectors ${\boldsymbol{\alpha}}^{\star}$ with about $15$ non-zero coefficients. The dictionaries associated with each digit class are initialized using the unsupervised formulation of Section 2. To test our algorithm with different values of $\mu$, we use a continuation strategy: Starting with $\mu=1.0$, we sequentially decrease its value by $0.1$ until we have $\mu=0$, learning with $10\,000$ iterations for each new value of $\mu$. We report the error rates in Figure 1, showing that our approach offers a competitive performance similar to . The best error rates of our method for $n=300,1000,5000$ labeled data are respectively $5.81,3.55$ and $1.81\%$, which is similar to who has reported $7.18,3.21$ and $1.52\%$ with the same sets of labeled data.

3 Learning a Nonlinear Image Mapping

We now illustrate our method in a regression context by considering a classical image processing task called “inverse halftoning”. With the development of several binary display technologies in the 70s (including, for example, printers and PC screens), the problem of converting a grayscale continuous-tone image into a binary one that looks perceptually similar to the original one (“halftoning”) was posed to the image processing community. Examples of halftoned images obtained with the classical Floyd-Steinberg algorithm are presented in the second column of Figure 2, with original images in the first column. Restoring these binary images to continuous-tone ones (“inverse halftoning”) has become a classical problem (see and references therein).

Unlike most image processing approaches that explicitly model the halftoning process, we formulate it as a regression problem, without exploiting any prior on the task. We use a database of $36$ images; $24$ are high-quality images from the Kodak PhotoCD datasethttp://r0k.us/graphics/kodak/ and are used for training, and $12$ are classical images often used for evaluating image processing algorithms;The list of these images can be found in , where they are used for the problem of image denoising. the first four (house, peppers, cameraman, lena) are used for validation and the remaining eight for testing.

We apply the Floyd-Steinberg algorithm implemented in the LASIP Matlab toolboxhttp://www.cs.tut.fi/~lasip/ to the grayscale continuous-tone images in order to build our training/validation/testing set. We extract all pairs of patches from the original/halftoned images in the training set, which provides us with a database of approximately $9$ millions of patches. We then use the “signal regression” formulation of Eq. (12) to learn a dictionary ${\mathbf{D}}$ and model parameters ${\mathbf{W}}$, by performing two passes of our algorithm over the $9$ million training pairs.

At this point, we have learned how to restore a small patch from an image, but not yet how to restore a full image. Following other patch-based approaches to image restoration , we extract from a test image all patches including overlaps, and restore each patch independently so that we get different estimates for each pixel (one estimate for each patch the pixel belongs to). These estimates are then averaged to reconstruct the full image, which has proven to give very good results in many image restoration tasks (see, e.g., ). The final image is then post-processed using the denoising algorithm of to remove possible artefacts.

We then measure how well it reconstructs the continuous-tone images from the halftoned ones in the test set. To reduce the number of hyperparameters, we have made a few arbitrary choices: We first use the Lasso formulation for encoding the signals—that is, we set $\lambda_{2}=0$. With millions of training samples, our model is unlikely to overfit and the regularization parameter $\nu$ is set to as well. The remaining free parameters are the size $m$ of the patches, the size $p$ of the dictionary and the regularization parameter $\lambda_{1}$. These parameters are selected by minimizing the mean-squared error reconstruction on the validation set. We have tried patches of size $m=l\times l$, with $l\in\{6,8,10,12,14,16\}$, dictionaries of sizes $p=100$, $250$ and $500$ , and determined $\lambda_{1}$ by first trying values on the logarithmic scale $10^{i}$, $i=-3,2$, then refining this parameter on the scale $0.1,0.2,0.3,\ldots,1.0$. The best parameters found are $m=10\times 10$, $p=500$ and $\lambda_{1}=0.6$. Since the test procedure is slightly different from the training one (the test includes an averaging step to restore a full image whereas the training one does not), we have refined the value of $\lambda_{1}$, trying different values one an additive scale $\{0.4,0.45,\ldots,0.75,0.8\}$, and selected the value $\lambda_{1}=0.55$, which has given the best result on the validation set.

Note that the largest dictionary has been chosen, and better results could potentially be obtained using an even larger dictionary, but this would imply a higher computational cost. Examples of results are presented in Figure 2. Halftoned images are binary but look perceptually similar to the original image. Reconstructed images have very few artefacts and most details are well preserved. We report in Table II a quantitative comparison between our approach and various ones from the literature, including the state-of-the-art algorithm of , which had until now the best results on this dataset. Even though our method does not explicitly model the transformation, it leads to better results in terms of PSNR.Denoting by MSE the mean-squared-error for images whose intensities are between and $255$, the PSNR is defined as $\textrm{PSNR}=10\log_{10}(255^{2}/\textrm{MSE})$ and is measured in dB. A gain of $1$dB reduces the MSE by approximately $20\%$. We also present in Figure 3 the results obtained by applying our algorithm to various binary images found on the web, from which we do not know the ground truth, and which have not necessarily been obtained with the Floyd-Steinberg algorithm. The results are qualitatively rather good.

From this experiment, we conclude that our method is well suited to (at least, some) nonlinear regression problems on natural images, paving the way to new applications of sparse image coding.

4 Digital Art Authentification

Recognizing authentic paintings from imitations using statistical techniques has been the topic of a few recent works . Classical methods compare, for example, the kurtosis of wavelet coefficients between a set of authentic paintings and imitations , or involve more sophisticated features . Recently, Hugues et al. have considered a dataset of $8$ authentic paintings from Pieter Bruegel the Elder, and $5$ imitations.The origin of these paintings is assessed by art historians. They have proposed to learn dictionaries for sparse coding, and use the kurtosis of the sparse coefficients as discriminative features. We use their dataset, which they kindly provided to us.It would have been interesting to use the datasets used in , but they are not publicly available.

The supervised dictionary learning approach we have presented is designed for classifying relatively small signals, and should not be directly applicable to the classification of large images, for which classical computer vision approaches based on bags of words may be better adapted (see for such approaches). However, we show that, for this particular dataset, a simple voting scheme based on the classification of small image patches with our method leads to good results.

The experiment we carry out consists of finding which painting is authentic, and which one is fake, in a pair known to contain one of each.This task is of course considerably easier than classifying each painting as authentic or fake. We do not claim to propose a method that readily applies to forgery detection. We proceed in a leave-one-out fashion, where we remove for testing one authentic and one imitation paintings from the dataset, and learn on the remaining ones. Since the dataset is small and does not have an official training/test set, we measure a cross-validation score, testing all possible pairs. We consider $12\times 12$ color image patches, without any pre-processing, and classify each patch from the test images independently. Then, we use a simple majority vote among the test patches to decide which image is the authentic one in the pair test, and which one is the imitation.Note that this experimental setting is different from , where only authentic paintings are used for training (and not imitations). We therefore do not make quantitive comparison with this work.

For each pair of authentic/imitation paintings, we build a dataset containing $200\,000$ patches from the authentic images and $200\,000$ from the imitations. We use the formulation from Eq. (13) for binary classification, and arbitrarily choose dictionaries containing $p=100$ dictionary elements. Since the training set is large, we set the parameter $\nu$ to , also choose the Lasso formulation for decomposing the patches by setting $\lambda_{2}=0$, and cross-validate on the parameter $\lambda_{1}$, trying values on a grid $\{10^{-4},10^{-3},\ldots,10^{0}\}$, and then refining the result on a grid with a logarithmic scale of $2$. We also compare Eq. (13) with the logistic regression loss and the basic formulation of Eq. (5) where ${\mathbf{D}}$ is learned unsupervised.

For classifying individual patches, the cross-validation score of the supervised formulation is a classification rate of $54.04\pm 2.26\%$, which slightly improves upon the “unsupervised” one that achieves $51.94\pm 1.92\%$. The task of classifying independently small image patches is difficult since there is significant overlap between the two classes. On the other hand, finding the imitation in a pair of (authentic,imitation) paintings with the voting scheme is easier and the “unsupervised formulation” only fails for one pair, whereas the supervised one has always given the right answer in our experiments.

5 Compressed Sensing

In this experiment, we apply our method to the problem of learning dictionaries and projection matrices for compressed sensing. As explained in Section 3.3.4, our formulation and the conclusions of this section hold for relatively small signals, where the sensing matrix can be stored into memory and learned. Thus, we consider here small image patches of natural images of size $m=10\times 10$ pixels. To build our training/validation/test set, we have chosen the Pascal VOC 2006 database of natural images : Images $1$ to $3000$ are used for training; images $3001$ to $4000$ are used for validation, and the remaining $1304$ images are kept for testing. Images are downsampled by a factor $2$ so that the JPEG compression artefacts present in this dataset become visually imperceptible, thereby improving its quality for our experiment.

We consider the case, which we call RANDOM, where the entries of ${\mathbf{Z}}$ are i.i.d. samples of the Gaussian distribution ${\mathcal{N}}(0,1/\sqrt{m})$. Since the purpose of ${\mathbf{Z}}$ is to reduce the dimensionality of the input data, it is also natural to consider the case where ${\mathbf{Z}}$ is obtained by principal component analysis on natural image patches (PCA). Finally, we also learn ${\mathbf{Z}}$ with the supervised learning formulation of Eq. (15), (SL), but consider the case where it is initialized randomly (SL1) or by PCA (SL2).

The matrix ${\mathbf{D}}$ can either be fixed or learned. A typical setting would be to set ${\mathbf{D}}={\mathbf{Z}}{\mathbf{D}}^{\prime}$, where ${\mathbf{D}}^{\prime}$ is discrete-cosine-transform matrix (DCT), often used in signal processing applications . It can also be learned with an unsupervised learning formulation (UL), or a supervised one (SL).

${\mathbf{W}}$ is always learned in a supervised way.

Since our training set is very large (several millions of patches), we arbitrarily set the regularization parameters $\nu_{1}$ and $\nu_{2}$ to . We measure reconstruction errors with dictionaries of various levels of overcompleteness by choosing a size $p$ in $\{100,200,400\}$. The remaining free parameters are the regularization parameters $\lambda_{1}$ and $\lambda_{2}$ for obtaining the coefficients ${\boldsymbol{\alpha}}^{\star}$. We try the values $\lambda_{1}=10^{i}$, with $i$ in $\{-5,\ldots,0\}$. Unlike what we have done in the experiments of Section 5.3, it is absolutely necessary in this setting to use $\lambda_{2}>0$ (according to the theory), since using a zero value for this parameter has led to instabilities and prevented our algorithm from converging. We have tried the values $\lambda_{2}=10^{i}\lambda_{1}$, with $i$ in $\{-2,-1,0\}$. Each learning procedure is performed by our algorithm in one pass on $10$ millions of patches randomly extracted from our training images. The pair of parameters $(\lambda_{1},\lambda_{2})$ that gives the lowest reconstruction error on the validation set is selected, and we report in Table III the result obtained on a test set of $500\,000$ patches randomly extracted from the $1304$ test images. The conclusions of this compressed sensing experiment on natural image patches are the following:

When ${\mathbf{Z}}$ is initialized as a Gaussian random matrix (case RANDOM), learning ${\mathbf{D}}$ and ${\mathbf{Z}}$ significantly improves the reconstruction error (case SL1). A similar observation was made in .

Results obtained with PCA are in general much better than those obtained with random projections, which is consistent with the conclusions of .

However, PCA does better than SL1. When PCA is used for initializing our supervised formulation, results can be slightly improved (case SL2). This illustrates either the limits of the non-convex optimization procedure, or that PCA is particularly well adapted to this problem.

Learned dictionaries (cases UL and SL) outperform classical DCT dictionaries.

Conclusion

We have presented in this paper a general formulation for learning sparse data representations tuned to specific tasks. Unlike classical approaches that learn a dictionary adapted to the reconstruction of the input data, our method learns features in a supervised way. We have shown that this approach is effective for solving classification and regression tasks in a large-scale setting, allowing the use of millions of training samples, and is able of exploiting successfully unlabeled data, when only a few labeled samples are available. Experiments on handwritten digits classification, non-linear inverse image mapping, digital art authentification, and compressed sensing have shown that our method leads to state-of-the-art results for several real problems. Future work will include adapting our method to various image processing problems such as image deblurring and image super-resolution, and other inverse problems. [Proofs and Lemmas] Before giving the proof of Proposition 1, we present two general results on the elastic net formulation of .

The vector ${\boldsymbol{\alpha}}^{\star}$ is a solution of Eq. (4) if and only if for all $j$ in $\{1,\ldots,p\}$,

Denoting by $\Lambda\stackrel{{\scriptstyle\vartriangle}}{{=}}\{j\in\{1,\ldots,p\}~{}~{}\text{s.t.}~{}~{}{\boldsymbol{\alpha}}^{\star}[j]\neq 0\}$ the active set, we also have

where ${\mathbf{s}}_{\Lambda}$ in $\{-1;+1\}^{|\Lambda|}$ carries the signs of ${\boldsymbol{\alpha}}^{\star}_{\Lambda}$.

Equation (20) can be obtained by considering subgradient optimality conditions as done in for the case $\lambda_{2}=0$. These can be written as

The next proposition exploits these optimality conditions to characterize the regularity of ${\boldsymbol{\alpha}}^{\star}$.

The function ${\boldsymbol{\alpha}}^{\star}$ is uniformly Lipschitz on ${\mathcal{X}}\times{\mathcal{D}}$.

Let ${\mathbf{D}}$ be in ${\mathcal{D}}$, $\varepsilon$ be a positive scalar and ${\mathbf{s}}$ be a vector in $\{-1,0,+1\}^{p}$, and define $K_{\mathbf{s}}({\mathbf{D}},\varepsilon)\subseteq{\mathcal{X}}$ as the set of vectors ${\mathbf{x}}$ satisfying for all $j$ in $\{1,\ldots,p\}$,

where ${\boldsymbol{\alpha}}^{\star}$ is shorthand for ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$.

Then, there exists $\kappa>0$ independent of ${\mathbf{s}}$, ${\mathbf{D}}$ and $\varepsilon$ so that for all ${\mathbf{x}}$ in $K_{\mathbf{s}}({\mathbf{D}},\varepsilon)$, the function ${\boldsymbol{\alpha}}^{\star}$ is twice continuously differentiable on $B_{\kappa\varepsilon}({\mathbf{x}})\times B_{\kappa\varepsilon}({\mathbf{D}})$, where $B_{\kappa\varepsilon}({\mathbf{x}})$ and $B_{\kappa\varepsilon}({\mathbf{D}})$ denote the open balls of radius $\kappa\varepsilon$ respectively centered on ${\mathbf{x}}$ and ${\mathbf{D}}$.

The first point is proven in . The proof uses the strong convexity induced by the elastic-net term, when $\lambda_{2}>0$, and the compactness of ${\mathcal{X}}$ from Assumption (A).

For the second point, we study the differentiability of ${\boldsymbol{\alpha}}^{\star}$ on sets that satisfy conditions which are more restrictive than the optimality conditions of Eq. (20). Concretely, let ${\mathbf{D}}$ be in ${\mathcal{D}}$, $\varepsilon>0$ and ${\mathbf{s}}$ be in $\{-1,0,+1\}^{p}$. The set $K_{\mathbf{s}}({\mathbf{D}},\varepsilon)$ characterizes the vectors ${\mathbf{x}}$ so that ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$ has the same signs as ${\mathbf{s}}$ (and same set of zero coefficients), and ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}},{\mathbf{D}})$ satisfies the conditions of Eq. (20), but with two additional constraints: (i) The magnitude of the non-zero coefficients in ${\boldsymbol{\alpha}}^{\star}$ should be greater than $\varepsilon$. (ii) The inequalities in Eq. (20) should be strict with a margin $\varepsilon$. The reason for imposing these assumptions is to restrict ourselves to points ${\mathbf{x}}$ in ${\mathcal{X}}$ that have a stable active set—that is, the set of non-zero coefficients $\Lambda$ of ${\boldsymbol{\alpha}}^{\star}$ should not change for small perturbations of $({\mathbf{x}},{\mathbf{D}})$, when ${\mathbf{x}}$ is in $K_{\mathbf{s}}({\mathbf{D}},\varepsilon)$.

Proving that there exists a constant $\kappa>0$ satisfying the second point is then easy (if a bit technical): Let us assume that $K_{\mathbf{s}}({\mathbf{D}},\varepsilon)$ is not empty (the case when it is empty is trivial). Since ${\boldsymbol{\alpha}}^{\star}$ is uniformly Lipschitz with respect to $({\mathbf{x}},{\mathbf{D}})$, so are the quantities ${\mathbf{d}}_{j}^{\top}({\mathbf{x}}-{\mathbf{D}}{\boldsymbol{\alpha}}^{\star})-\lambda_{2}{\boldsymbol{\alpha}}^{\star}[j]$ and ${\mathbf{s}}[j]{\boldsymbol{\alpha}}^{\star}[j]$, for all $j$ in $\{1,\ldots,p\}$. Thus, there exists $\kappa>0$ independent of ${\mathbf{x}}$ and ${\mathbf{D}}$ such that for all $({\mathbf{x}}^{\prime},{\mathbf{D}}^{\prime})$ satisfying $\|{\mathbf{x}}-{\mathbf{x}}^{\prime}\|_{2}\leq\kappa\varepsilon$ and $\|{\mathbf{D}}-{\mathbf{D}}^{\prime}\|_{F}\leq\kappa\varepsilon$, we have for all $j$ in $\{1,\ldots,p\}$,

where ${\boldsymbol{\alpha}}^{\star\prime}$ is short-hand for ${\boldsymbol{\alpha}}^{\star}({\mathbf{x}}^{\prime},{\mathbf{D}}^{\prime})$, and ${\mathbf{x}}^{\prime}$ is therefore in $K_{{\mathbf{s}}}({\mathbf{D}}^{\prime},\varepsilon/2)$. It is then easy to show that the active set $\Lambda$ of ${\boldsymbol{\alpha}}^{\star}$ and the signs of ${\boldsymbol{\alpha}}^{\star}$ are stable on $B_{\kappa\varepsilon}({\mathbf{x}})\times B_{\kappa\varepsilon}({\mathbf{D}})$, and that ${\boldsymbol{\alpha}}_{\Lambda}^{\star}$ is given by the closed form of Eq. (21). ${\boldsymbol{\alpha}}^{\star}$ is therefore twice differentiable on $B_{\kappa\varepsilon}({\mathbf{x}})\times B_{\kappa\varepsilon}({\mathbf{D}})$. ∎

With this proposition in hand, we can now present the proof of Proposition 1:

Let now choose ${\mathbf{W}}$ in ${\mathcal{W}}$ and ${\mathbf{D}}$ in ${\mathcal{D}}$. We have characterized in Lemma 3 the differentiability of ${\boldsymbol{\alpha}}^{\star}$ on some subsets of ${\mathcal{X}}\times{\mathcal{D}}$. We consider the set

where $g$ has the form given by Eq. (16). This shows that $f$ is differentiable with respect to ${\mathbf{D}}$, and its gradient $\nabla_{\mathbf{D}}f$ is $g$. ∎

Acknowledgments

This paper was supported in part by ANR under grant MGA ANR-07-BLAN-0311 and the European Research Council (SIERRA and VideoWorld Project). Julien Mairal is now supported by the NSF grant SES-0835531 and NSF award CCF-0939370. The authors would like to thank J.M. Hugues and D.J. Graham and D.N. Rockmore for providing us with the Bruegel dataset used in Section 5.4, and Y-Lan Boureau and Marc’Aurelio Ranzato for providing their experimental setting for the digit recognition task.

References