Modeling signals as sparse linear combinations of atoms selected from a dictionary has become a popular paradigm in many fields, including signal processing, statistics, and machine learning. This line of research, also known assparse coding, has witnessed the development of several well-founded theoretical frameworks (Tibshirani, 1996; Chen et al., 1998; Mallat, 1999; Tropp, 2004, 2006; Wainwright, 2009; Bickel et al., 2009) and the emergence of many efficient algorithmic tools (Efron et al., 2004; Nesterov, 2007; Beck and Teboulle, 2009; Wright et al., 2009; Needell and Tropp, 2009; Yuan et al., 2010).
In many applied settings, the structure of the problem at hand, such as, e.g., the spatial arrangement of the pixels in an image, or the presence of variables corresponding to several levels of a given factor, induces relationships between dictionary elements. It is appealing to use this a priori knowledge about the problem directly to constrain the possible sparsity patterns. For instance, when the dictionary elements are partitioned into predefined groups corresponding to different types of features, one can enforce a similar block structure in the sparsity pattern—that is, allow only that either all elements of a group are part of the signal decomposition or that all are dismissed simultaneously (see Yuan and Lin, 2006; Stojnic et al., 2009).
This example can be viewed as a particular instance of structured sparsity, which has been lately the focus of a large amount of research (Baraniuk et al., 2010; Zhao et al., 2009; Huang et al., 2009; Jacob et al., 2009; Jenatton et al., 2009; Micchelli et al., 2010). In this paper, we concentrate on a specific form of structured sparsity, which we call hierarchical sparse coding: the dictionary elements are assumed to be embedded in a directed tree , and the sparsity patterns are constrained to form a connected and rooted subtree of (Donoho, 1997; Baraniuk, 1999; Baraniuk et al., 2002, 2010; Zhao et al., 2009; Huang et al., 2009). This setting extends more generally to a forest of directed trees.111A tree is defined as a connected graph that contains no cycle (see Ahuja et al., 1993).
In fact, such a hierarchical structure arises in many applications. Wavelet decompositions lend themselves well to this tree organization because of their multiscale structure, and benefit from it for image compression and denoising (Shapiro, 1993; Crouse et al., 1998; Baraniuk, 1999; Baraniuk et al., 2002, 2010; He and Carin, 2009; Zhao et al., 2009; Huang et al., 2009). In the same vein, edge filters of natural image patches can be represented in an arborescent fashion (Zoran and Weiss, 2009). Imposing these sparsity patterns has further proven useful in the context of hierarchical variable selection, e.g., when applied to kernel methods (Bach, 2008), to log-linear models for the selection of potential orders (Schmidt and Murphy, 2010), and to bioinformatics, to exploit the tree structure of gene networks for multi-task regression (Kim and Xing, 2010)2009, and references therein) have also emerged as a natural structure in several applications, notably to model text documents. In particular, in the context of topic models (Blei et al., 2003), a hierarchical model of latent variables based on Bayesian non-parametric methods has been proposed by Blei et al. (2010) to model hierarchies of topics.
To perform hierarchical sparse coding, our work builds upon the approach of Zhao et al. (2009) who first introduced a sparsity-inducing norm leading to this type of tree-structured sparsity pattern. We tackle the resulting nonsmooth convex optimization problem with proximal methods (e.g., Nesterov, 2007; Beck and Teboulle, 2009; Wright et al., 2009; Combettes and Pesquet, 2010) and we show in this paper that its key step, the computation of the proximal operator, can be solved exactly with a complexity linear, or close to linear, in the number of dictionary elements—that is, with the same complexity as for classical -sparse decomposition problems (Tibshirani, 1996; Chen et al., 1998). Concretely, given an -dimensional signal along with a dictionary composed of atoms, the optimization problem at the core of our paper can be written as
In this formulation, the sparsity-inducing norm encodes a hierarchical structure among the atoms of , where this structure is assumed to be known beforehand. The precise meaning of hierarchical structure and the definition of will be made more formal in the next sections. A particular instance of this problem—known as the proximal problem—is central to our analysis and concentrates on the case where the dictionary is orthogonal.
In addition to a speed benchmark that evaluates the performance of our proposed approach in comparison with other convex optimization techniques, two types of applications and experiments are considered. First, we consider settings where the dictionary is fixed and given a priori, corresponding for instance to a basis of wavelets for the denoising of natural images. Second, we show how one can take advantage of this hierarchical sparse coding in the context of dictionary learning (Olshausen and Field, 1997; Aharon et al., 2006; Mairal et al., 2010a), where the dictionary is learned to adapt to the predefined tree structure. This extension of dictionary learning is notably shown to share interesting connections with hierarchical probabilistic topic models.
To summarize, the contributions of this paper are threefold:
We show that the proximal operator for a tree-structured sparse regularization can be computed exactly in a finite number of operations using a dual approach. Our approach is equivalent to computing a particular sequence of elementary proximal operators, and has a complexity linear, or close to linear, in the number of variables. Accelerated gradient methods (e.g., Nesterov, 2007; Beck and Teboulle, 2009; Combettes and Pesquet, 2010) can then be applied to solve large-scale tree-structured sparse decomposition problems at the same computational cost as traditional ones using the -norm.
We propose to use this regularization scheme to learn dictionaries embedded in a tree, which, to the best of our knowledge, has not been done before in the context of structured sparsity.
Our method establishes a bridge between hierarchical dictionary learning and hierarchical topic models (Blei et al., 2010), which builds upon the interpretation of topic models as multinomial PCA (Buntine, 2002), and can learn similar hierarchies of topics. This point is discussed in Sections 5.5 and 6.
Note that this paper extends a shorter version published in the proceedings of the international conference of machine learning (Jenatton et al., 2010).
Vectors are denoted by bold lower case letters and matrices by upper case ones. We define for the -norm of a vector in as , where denotes the -th coordinate of , and . We also define the -pseudo-norm as the number of nonzero elements in a vector:222Note that it would be more proper to write instead of to be consistent with the traditional notation . However, for the sake of simplicity, we will keep this notation unchanged in the rest of the paper. . We consider the Frobenius norm of a matrix in : , where denotes the entry of at row and column . Finally, for a scalar , we denote .
The rest of this paper is organized as follows: Section 2 presents related work and the problem we consider. Section 3 is devoted to the algorithm we propose, and Section 4 introduces the dictionary learning framework and shows how it can be used with tree-structured norms. Section 5 presents several experiments demonstrating the effectiveness of our approach and Section 6 concludes the paper.
2 Problem Statement and Related Work
Let us consider an input signal of dimension , typically an image described by its pixels, which we represent by a vector in . In traditional sparse coding, we seek to approximate this signal by a sparse linear combination of atoms, or dictionary elements, represented here by the columns of a matrix in . This can equivalently be expressed as for some sparse vector in , i.e, such that the number of nonzero coefficients is small compared to . The vector is referred to as the code, or decomposition, of the signal .
In the rest of the paper, we focus on specific sets of nonzero coefficients—or simply, nonzero patterns—for the decomposition vector . In particular, we assume that we are given a tree333Our analysis straightforwardly extends to the case of a forest of trees; for simplicity, we consider a single tree . whose nodes are indexed by in . We want the nonzero patterns of to form a connected and rooted subtree of ; in other words, if denotes the set of indices corresponding to the ancestors444We consider that the set of ancestors of a node also contains the node itself. of the node in (see Figure 1), the vector obeys the following rule
Informally, we want to exploit the structure of in the following sense: the decomposition of any signal can involve a dictionary element only if the ancestors of in the tree are themselves part of the decomposition.
We now review previous work that has considered the sparse approximation problem with tree-structured constraints (1). Similarly to traditional sparse coding, there are basically two lines of research, that either (A) deal with nonconvex and combinatorial formulations that are in general computationally intractable and addressed with greedy algorithms, or (B) concentrate on convex relaxations solved with convex programming methods.
2.1 Nonconvex Approaches
For a given sparsity level (number of nonzero coefficients), the following nonconvex problem
has been tackled by Baraniuk (1999); Baraniuk et al. (2002) in the context of wavelet approximations with a greedy procedure. A penalized version of problem (2) (that adds to the objective function in place of the constraint ) has been considered by Donoho (1997), while studying the more general problem of best approximation from dyadic partitions (see Section 6 in Donoho, 1997). Interestingly, the algorithm we introduce in Section 3 shares conceptual links with the dynamic-programming approach of Donoho (1997), which was also used by Baraniuk et al. (2010), in the sense that the same order of traversal of the tree is used in both procedures. We investigate more thoroughly the relations between our algorithm and this approach in Appendix A.
Problem (2) has been further studied for structured compressive sensing (Baraniuk et al., 2010), with a greedy algorithm that builds upon Needell and Tropp (2009). Finally, Huang et al. (2009) have proposed a formulation related to (2), with a nonconvex penalty based on an information-theoretic criterion.
2.2 Convex Approach
2.2.1 Hierarchical Sparsity-Inducing Norms
Condition (1) can be equivalently expressed by its contrapositive, thus leading to an intuitive way of penalizing the vector to obtain tree-structured nonzero patterns. More precisely, defining analogously to for in , condition (1) amounts to saying that if a dictionary element is not used in the decomposition, its descendants in the tree should not be used either. Formally, this can be formulated as:
From now on, we denote by the set defined by and refer to each member of as a group (Figure 2). To obtain a decomposition with the desired property (3), one can naturally penalize the number of groups in that are “involved” in the decomposition of , i.e., that record at least one nonzero coefficient of :
While this intuitive penalization is nonconvex (and not even continuous), a convex proxy has been introduced by Zhao et al. (2009). It was further considered by Bach (2008); Kim and Xing (2010); Schmidt and Murphy (2010) in several different contexts. For any vector , let us define
where is the vector of size whose coordinates are equal to those of for indices in the set , and to otherwise555Note the difference with the notation , which is often used in the literature on structured sparsity, where is a vector of size .. The notation stands in practice either for the - or -norm, and denotes some positive weights666For a complete definition of for any -norm, a discussion of the choice of , and a strategy for choosing the weights (see Zhao et al., 2009; Kim and Xing, 2010).. As analyzed by Zhao et al. (2009) and Jenatton et al. (2009), when penalizing by , some of the vectors are set to zero for some .777It has been further shown by Bach (2010) that the convex envelope of the nonconvex function of Eq. (4) is in fact with being the -norm. Therefore, the components of corresponding to some complete subtrees of are set to zero, which exactly matches condition (3), as illustrated in Figure 2.
Note that although we presented for simplicity this
hierarchical norm in the context of a single tree with a single element at
each node, it can easily be extended to the case of forests of trees, and/or
trees containing arbitrary numbers of dictionary elements at each node (with nodes possibly containing no dictionary element).
More broadly, this formulation can be extended with the notion
of tree-structured groups, which we now present:
[Tree-structured set of groups.]
A set of groups is said to be tree-structured in , if and if for all , For such a set of groups, there exists a (non-unique) total order relation such that:
Given such a tree-structured set of groups and its associated norm , we are interested throughout the paper in the following hierarchical sparse coding problem,
where is the tree-structured norm we have previously introduced, the non-negative scalar is a regularization parameter controlling the sparsity of the solutions of (5), and
a smooth convex loss function (see Section3 for more details about the smoothness assumptions on ). In the rest of the paper, we will mostly use the square loss with a dictionary in , but the formulation of Eq. (5) extends beyond this context. In particular one can choose to be the logistic loss, which is commonly used for classification problems (e.g., see Hastie et al., 2009).
Before turning to optimization methods for the hierarchical sparse coding problem, we consider a particular instance. The sparse group Lasso was recently considered by Sprechmann et al. (2010) and Friedman et al. (2010) as an extension of the group Lasso of Yuan and Lin (2006). To induce sparsity both groupwise and within groups, Sprechmann et al. (2010) and Friedman et al. (2010) add an term to the regularization of the group Lasso, which given a partition of in disjoint groups yields a regularized problem of the form
Since is a partition, the set of groups in and the singletons form together a tree-structured set of groups according to definition 2 and the algorithm we will develop is therefore applicable to this problem.
2.2.2 Optimization for Hierarchical Sparsity-Inducing Norms
While generic approaches like interior-point methods (Boyd and Vandenberghe, 2004) and subgradient descent schemes (Bertsekas, 1999) might be used to deal with the nonsmooth norm , several dedicated procedures have been proposed.
In Zhao et al. (2009), a boosting-like technique is used, with a path-following strategy in the specific case where is the -norm. Based on the variational equality
Kim and Xing (2010) follow a reweighted least-square scheme that is well adapted to the square loss function. To the best of our knowledge, a formulation of this type is however not available when is the -norm. In addition it requires an appropriate smoothing to become provably convergent. The same approach is considered by Bach (2008), but built upon an active-set strategy. Other proposed methods consist of a projected gradient descent with approximate projections onto the ball (Schmidt and Murphy, 2010), and an augmented-Lagrangian based technique (Sprechmann et al., 2010) for solving a particular case with two-level hierarchies.
While the previously listed first-order approaches are (1) loss-function dependent, and/or (2) not guaranteed to achieve optimal convergence rates, and/or (3) not able to yield sparse solutions without a somewhat arbitrary post-processing step, we propose to resort to proximal methods888Note that the authors of Chen et al. (2010) have considered proximal methods for general group structure when is the -norm; due to a smoothing of the regularization term, the convergence rate they obtained is suboptimal. that do not suffer from any of these drawbacks.
We begin with a brief introduction to proximal methods, necessary to present our contributions. From now on, we assume that is convex and continuously differentiable with Lipschitz-continuous gradient. It is worth mentioning that there exist various proximal schemes in the literature that differ in their settings (e.g., batch versus stochastic) and/or the assumptions made on . For instance, the material we develop in this paper could also be applied to online/stochastic frameworks (Duchi and Singer, 2009; Hu et al., 2009; Xiao, 2010) and to possibly nonsmooth functions (e.g., Duchi and Singer, 2009; Xiao, 2010; Combettes and Pesquet, 2010, and references therein). Finally, most of the technical proofs of this section are presented in Appendix B for readability.
3.1 Proximal Operator for the Norm
Proximal methods have drawn increasing attention in the signal processing (e.g., Becker et al., 2009; Wright et al., 2009; Combettes and Pesquet, 2010, and numerous references therein) and the machine learning communities (e.g., Bach et al., 2011, and references therein), especially because of their convergence rates (optimal for the class of first-order techniques) and their ability to deal with large nonsmooth convex problems (e.g., Nesterov, 2007; Beck and Teboulle, 2009). In a nutshell, these methods can be seen as a natural extension of gradient-based techniques when the objective function to minimize has a nonsmooth part. Proximal methods are iterative procedures. The simplest version of this class of methods linearizes at each iteration the function
around the current estimate, and this estimate is updated as the (unique by strong convexity) solution of the proximal problem, defined as follows:
The quadratic term keeps the update in a neighborhood where is close to its linear approximation, and is a parameter which is an upper bound on the Lipschitz constant of . This problem can be equivalently rewritten as:
Solving efficiently and exactly this problem is crucial to
enjoy the fast convergence rates of proximal methods.
In addition, when the nonsmooth term is not present, the previous
proximal problem exactly leads to the standard gradient update rule.
More generally, we define the proximal operator:
The proximal operator associated with our regularization term , which we denote by , is the function that maps a vector to the unique solution of
This operator was initially introduced by Moreau (1962) to generalize the projection operator onto a convex set. What makes proximal methods appealing for solving sparse decomposition problems is that this operator can be often computed in closed-form. For instance,
When is the -norm—that is, , the proximal operator is the well-known elementwise soft-thresholding operator,
When is a group-Lasso penalty with -norms—that is, , with being a partition of , the proximal problem is separable in every group, and the solution is a generalization of the soft-thresholding operator to groups of variables:
where denotes the orthogonal projection onto the ball of the -norm of radius .
When is a group-Lasso penalty with -norms—that is, , the solution is also a group-thresholding operator:
where denotes the orthogonal projection onto the -ball of radius , which can be solved in operations (Brucker, 1984; Maculan and Galdino de Paula, 1989). Note that when , we have a group-thresholding effect, with .
More generally, a classical result (see, e.g., Combettes and Pesquet, 2010; Wright et al., 2009) says that the proximal operator for a norm can be computed as the residual of the projection of a vector onto a ball of the dual-norm denoted by , and defined for any vector in by .999It is easy to show that the dual norm of the -norm is the -norm itself. The dual norm of the is the -norm. This is a classical duality result for proximal operators leading to the different closed forms we have just presented. We have indeed that and , where Id stands for the identity operator. Obtaining closed forms is, however, not possible anymore as soon as some groups in overlap, which is always the case in our hierarchical setting with tree-structured groups.
3.2 A Dual Formulation of the Proximal Problem
We now show that Eq. (7) can be solved using a dual approach, as
described in the following lemma. The result relies on conic duality
(Boyd and Vandenberghe, 2004), and does not make any assumption on the choice of the norm :
[Dual of the proximal problem]
Let and let us consider the problem
where and denotes the -th coordinate of the vector in . Then, problems (7) and (8) are dual to each other and strong duality holds. In addition, the pair of primal-dual variables is optimal if and only if is a feasible point of the optimization problem (8), and
where we denote by the orthogonal projection onto the ball . Note that we focus here on specific tree-structured groups, but the previous lemma is valid regardless of the nature of . The rationale of introducing such a dual formulation is to consider an equivalent problem to (7) that removes the issue of overlapping groups at the cost of a larger number of variables. In Eq. (7), one is indeed looking for a vector of size , whereas one is considering a matrix in in Eq. (8) with nonzero entries, but with separable (convex) constraints for each of its columns.
This specific structure makes it possible to use block coordinate ascent (Bertsekas, 1999). Such a procedure is presented in Algorithm 1. It optimizes sequentially Eq. (8) with respect to the variable , while keeping fixed the other variables , for . It is easy to see from Eq. (8) that such an update of a column , for a group in , amounts to computing the orthogonal projection of the vector onto the ball of radius of the dual norm .
3.3 Convergence in One Pass
In general, Algorithm 1 is not guaranteed to solve exactly Eq. (7) in a finite number of iterations. However, when is the - or -norm, and provided that the groups in are appropriately ordered, we now prove that only one pass of Algorithm 1, i.e., only one iteration over all groups, is sufficient to obtain the exact solution of Eq. (7). This result constitutes the main technical contribution of the paper and is the key for the efficiency of our procedure.
Before stating this result, we need to introduce a lemma showing that,
given two nested groups such that ,
if is updated before in Algorithm 1, then
the optimality condition for is not perturbed by the update of .
[Projections with nested groups]
Let denote either the - or -norm, and and be two nested groups—that is, . Let be a vector in , and let us consider the successive projections
with . Let us introduce . The following relationships hold
The previous lemma establishes the convergence in one pass of Algorithm 1 in the case where only contains two nested groups , provided that is computed before . Let us illustrate this fact more concretely. After initializing and to zero, Algorithm 1 first updates with the formula , and then performs the following update: (where we have used that since ). We are now in position to apply Lemma 3.3 which states that the primal/dual variables satisfy the optimality conditions (9), as described in Lemma 3.2. In only one pass over the groups , we have in fact reached a solution of the dual formulation presented in Eq. (8), and in particular, the solution of the proximal problem (7).
In the following proposition, this lemma is extended to general tree-structured
sets of groups :
[Convergence in one pass]
Suppose that the groups in are ordered according to the total order relation of Definition 2, and that the norm is either the - or -norm. Then, after initializing to , a single pass of Algorithm 1 over with the order yields the solution of the proximal problem (7). The proof largely relies on Lemma 3.3 and proceeds by induction. By definition of Algorithm 1, the feasibility of is always guaranteed. We consider the following induction hypothesis
Since the dual variables are initially equal to zero, the summation over is equivalent to a summation over . We initialize the induction with the first group in , that, by definition of , does not contain any other group. The first step of Algorithm 1 easily shows that the induction hypothesis is satisfied for this first group.
We now assume that is true and consider the next group , , in order to prove that is also satisfied. We have for each group ,
Since for , we have
and following the update rule for the group ,
At this point, we can apply Lemma 3.3 for each group , which proves that the induction hypothesis is true. Let us introduce . We have shown that for all in , As a result, the pair satisfies the optimality conditions (9) of problem (8). Therefore, after one complete pass over , the primal/dual pair is optimal, and in particular, is the solution of problem (7). Using conic duality, we have derived a dual formulation of the proximal operator, leading to Algorithm 1 which is generic and works for any norm , as long as one is able to perform projections onto balls of the dual norm . We have further shown that when is the - or the -norm, a single pass provides the exact solution when the groups are correctly ordered. We show however in Appendix C, that, perhaps surprisingly, the conclusions of Proposition 3.3 do not hold for general -norms, if . Next, we give another interpretation of this result.
3.4 Interpretation in Terms of Composition of Proximal Operators
In Algorithm 1, since all the vectors are initialized to , when the group is considered, we have by induction . Thus, to maintain at each iteration of the inner loop one can instead update after updating according to . Moreover, since is no longer needed in the algorithm, and since only the entries of indexed by are updated, we can combine the two updates into , leading to a simplified Algorithm 2 equivalent to Algorithm 1.
Actually, in light of the classical relationship between proximal operator and projection (as discussed in Section 3.1), it is easy to show that each update is equivalent to . To simplify the notations, we define the proximal operator for a group in as for every vector in .
3.5 Efficient Implementation and Complexity
Since Algorithm 2 involves
projections on the dual balls (respectively the - and the -balls
for the - and -norms) of vectors in , in a first approximation, its complexity is at most , because each
of these projections can be computed in
operations (Brucker, 1984; Maculan and Galdino de Paula, 1989).
But in fact, the algorithm performs one projection for each group involving variables, and
the total complexity is therefore .
By noticing that if and are two groups with the same depth in the tree,
then , it is easy to show that the number of variables involved in all the projections
is less than or equal to , where is the depth of the tree:
[Complexity of Algorithm 2]
Algorithm 2 gives the solution of the primal problem Eq. (7) in operations, where is the depth of the tree. Lemma 3.5 should not suggest that the complexity is linear in , since could depend of as well, and in the worst case the hierarchy is a chain, yielding . However, in a balanced tree, . In practice, the structures we have considered experimentally are relatively flat, with a depth not exceeding , and the complexity is therefore almost linear.
Moreover, in the case of the -norm, it is actually possible to propose an algorithm with complexity . Indeed, in that case each of the proximal operators is a scaling operation: . The composition of these operators in Algorithm 1 thus corresponds to performing sequences of scaling operations. The idea behind Algorithm 3 is that the corresponding scaling factors depend only on the norms of the successive residuals of the projections and that these norms can be computed recursively in one pass through all nodes in operations; finally, computing and applying all scalings to each entry takes then again operations.
To formulate the algorithm, two new notations are used: for a group in , we denote by the indices of the variables that are at the root of the subtree corresponding to ,101010As a reminder, is not a singleton when several dictionary elements are considered per node. and by the set of groups that are the children of in the tree. For example, in the tree presented in Figure 2, , , , and . Note that all the groups of are necessarily included in .
The next lemma is proved in Appendix B.
[Correctness and complexity of Algorithm 3]
When is chosen to be the -norm, Algorithm 3 gives the solution of the primal problem Eq. (7) in operations. So far the dictionary was fixed to be for example a wavelet basis. In the next section, we apply the tools we developed for solving efficiently problem (5) to learn a dictionary adapted to our hierarchical sparse coding formulation.
4 Application to Dictionary Learning
We start by briefly describing dictionary learning.
4.1 The Dictionary Learning Framework
Let us consider a set in of signals of dimension . Dictionary learning is a matrix factorization problem which aims at representing these signals as linear combinations of the dictionary elements, that are the columns of a matrix in . More precisely, the dictionary is learned along with a matrix of decomposition coefficients in , so that for every signal .
While learning simultaneously and , one may want to encode specific prior knowledge about the problem at hand, such as, for example, the positivity of the decomposition (Lee and Seung, 1999), or the sparsity of (Olshausen and Field, 1997; Aharon et al., 2006; Lee et al., 2007; Mairal et al., 2010a). This leads to penalizing or constraining and results in the following formulation:
where and denote two convex sets and is a regularization term, usually a norm or a squared norm, whose effect is controlled by the regularization parameter . Note that is assumed to be bounded to avoid any degenerate solutions of Problem (10). For instance, the standard sparse coding formulation takes to be the -norm, to be the set of matrices in whose columns have unit -norm, with (Olshausen and Field, 1997; Lee et al., 2007; Mairal et al., 2010a).
However, this classical setting treats each dictionary element independently from the others, and does not exploit possible relationships between them. To embed the dictionary in a tree structure, we therefore replace the -norm by our hierarchical norm and set in Eq. (10).
A question of interest is whether hierarchical priors are more appropriate in supervised settings or in the matrix-factorization context in which we use it. It is not so common in the supervised setting to have strong prior information that allows us to organize the features in a hierarchy. On the contrary, in the case of dictionary learning, since the atoms are learned, one can argue that the dictionary elements learned will have to match well the hierarchical prior that is imposed by the regularization. In other words, combining structured regularization with dictionary learning has precisely the advantage that the dictionary elements will self-organize to match the prior.
4.2 Learning the Dictionary
Optimization for dictionary learning has already been intensively studied. We choose in this paper a typical alternating scheme, which optimizes in turn and while keeping the other variable fixed (Aharon et al., 2006; Lee et al., 2007; Mairal et al., 2010a).111111Note that although we use this classical scheme for simplicity, it would also be possible to use the stochastic approach proposed by Mairal et al. (2010a). Of course, the convex optimization tools we develop in this paper do not change the intrinsic non-convex nature of the dictionary learning problem. However, they solve the underlying convex subproblems efficiently, which is crucial to yield good results in practice. In the next section, we report good performance on some applied problems, and we show empirically that our algorithm is stable and does not seem to get trapped in bad local minima. The main difficulty of our problem lies in the optimization of the vectors , in , for the dictionary kept fixed. Because of , the corresponding convex subproblem is nonsmooth and has to be solved for each of the signals considered. The optimization of the dictionary (for fixed), which we discuss first, is in general easier.
Updating the dictionary .
We follow the matrix-inversion free procedure of Mairal et al. (2010a) to update the dictionary. This method consists in iterating block-coordinate descent over the columns of . Specifically, we assume that the domain set has the form
or with . The choice for these particular domain sets is motivated by the experiments of Section 5. For natural image patches, the dictionary elements are usually constrained to be in the unit -norm ball (i.e., ), while for topic modeling, the dictionary elements are distributions of words and therefore belong to the simplex (i.e., ). The update of each dictionary element amounts to performing a Euclidean projection, which can be computed efficiently (Mairal et al., 2010a). Concerning the stopping criterion, we follow the strategy from the same authors and go over the columns of only a few times, typically times in our experiments. Although we have not explored locality constraints on the dictionary elements, these have been shown to be particularly relevant to some applications such as patch-based image classification (Yu et al., 2009). Combining tree structure and locality constraints is an interesting future research.
Updating the vectors .
The procedure for updating the columns of is based on the results derived in Section 3.3. Furthermore, positivity constraints can be added on the domain of , by noticing that for our norm and any vector in , adding these constraints when computing the proximal operator is equivalent to solving This equivalence is proved in Appendix B.6. We will indeed use positive decompositions to model text corpora in Section 5. Note that by constraining the decompositions to be nonnegative, some entries may be set to zero in addition to those already zeroed out by the norm . As a result, the sparsity patterns obtained in this way might not satisfy the tree-structured condition (1) anymore.
We next turn to the experimental validation of our hierarchical sparse coding.
5.1 Implementation Details
In Section 3.3, we have shown that the proximal operator
associated to can be computed exactly and efficiently. The problem is therefore amenable to
fast proximal algorithms that are well suited to nonsmooth convex
Specifically, we tried the accelerated scheme from both
Nesterov (2007) and Beck and Teboulle (2009), and finally opted for the latter
since, for a comparable level of precision, fewer calls of the proximal
operator are required. The basic proximal scheme presented in Section 3.1
is formalized by Beck and Teboulle (2009) as an algorithm called ISTA; the same authors propose moreover
an accelerated variant, FISTA, which is a similar procedure, except that the
operator is not directly applied on the current estimate, but on an auxiliary
sequence of points that are linear combinations of past estimates. This latter algorithm has
an optimal convergence rate in the class of first-order techniques, and also
allows for warm restarts, which is crucial in the alternating scheme of dictionary
learning.121212 Unless otherwise specified,
the initial stepsize in ISTA/FISTA is chosen as the maximum eigenvalue of the sampling covariance matrix divided by 100, while
the growth factor in the line search is set to
Unless otherwise specified, the initial stepsize in ISTA/FISTA is chosen as the maximum eigenvalue of the sampling covariance matrix divided by 100, while the growth factor in the line search is set to.
Finally, we monitor the convergence of the algorithm by checking the relative decrease in the cost function.131313We are currently investigating algorithms for computing duality gaps based on network flow optimization tools (Mairal et al., 2010b). Unless otherwise specified, all the algorithms used in the following experiments are implemented in C/C++, with a Matlab interface. Our implementation is freely available at http://www.di.ens.fr/willow/SPAMS/.
5.2 Speed Benchmark
To begin with, we conduct speed comparisons between our approach and other
convex programming methods, in the setting where is chosen to be a
linear combination of -norms. The algorithms that take part in the
following benchmark are:
• Proximal methods, with ISTA and the accelerated FISTA methods (Beck and Teboulle, 2009).
• A reweighted-least-square scheme (Re-), as described by Jenatton et al. (2009); Kim and Xing (2010). This approach is adapted to the square loss, since closed-form updates can be used.141414The computation of the updates related to the variational formulation (6) also benefits from the hierarchical structure of , and can be performed in operations.
• Subgradient descent, whose step size is taken to be equal either to or (respectively referred to as SG and ), where is the iteration number, and are the best151515“The best step size” is understood as being the step size leading to the smallest cost function after 500 iterations. parameters selected on the logarithmic grid .
• A commercial software (Mosek, available at http://www.mosek.com/) for second-order cone programming (SOCP).
Moreover, the experiments we carry out cover various settings, with notably different sparsity regimes, i.e., low, medium and high, respectively corresponding to about and of the total number of dictionary elements. Eventually, all reported results are obtained on a single core of a 3.07Ghz CPU with 8GB of memory.
5.2.1 Hierarchical dictionary of natural image patches
In this first benchmark, we consider a least-squares regression problem regularized by that arises in the context of denoising of natural image patches, as further exposed in Section 5.4. In particular, based on a hierarchical dictionary, we seek to reconstruct noisy -patches. The dictionary we use is represented on Figure 17. Although the problem involves a small number of variables, i.e., dictionary elements, it has to be solved repeatedly for tens of thousands of patches, at moderate precision. It is therefore crucial to be able to solve this problem quickly and efficiently.
We can draw several conclusions from the results of the simulations reported in Figure 6. First, we observe that in most cases, the accelerated proximal scheme performs better than the other approaches. In addition, unlike FISTA, ISTA seems to suffer in non-sparse scenarios. In the least sparse setting, the reweighted- scheme is the only method that competes with FISTA. It is however not able to yield truly sparse solutions, and would therefore need a subsequent (somewhat arbitrary) thresholding operation. As expected, the generic techniques such as SG and SOCP do not compete with dedicated algorithms.
5.2.2 Multi-class classification of cancer diagnosis
The second benchmark explores a different supervised learning setting, whereis no longer the square loss function. The goal is to demonstrate that our optimization tools apply in various scenarios, beyond traditional sparse approximation problems. To this end, we consider a gene expression dataset161616The dataset we use is 14_Tumors, which is freely available at http://www.gems-system.org/. in the context of cancer diagnosis. More precisely, we focus on a multi-class classification problem where the number
of samples to be classified is small compared to the numberof gene expressions that characterize these samples. Each atom thus corresponds to a gene expression across the samples, whose class labels are recorded in the vector in .
The dataset contains samples, variables and classes. In addition, the data exhibit highly-correlated dictionary elements. Inspired by Kim and Xing (2010), we build the tree-structured set of groups
using Ward’s hierarchical clustering(Johnson, 1967) on the gene expressions. The norm built in this way aims at capturing the hierarchical structure of gene expression networks (Kim and Xing, 2010).
Instead of the square loss function, we consider the multinomial logistic loss function that is better suited to deal with multi-class classification problems (see, e.g., Hastie et al., 2009). As a direct consequence, algorithms whose applicability crucially depends on the choice of the loss function are removed from the benchmark. This is the case with reweighted- schemes that do not have closed-form updates anymore. Importantly, the choice of the multinomial logistic loss function leads to an optimization problem over a matrix with dimensions times the number of classes (i.e., a total of variables). Also, due to scalability issues, generic interior point solvers could not be considered here.
The results in Figure 10 highlight that the accelerated proximal scheme performs overall better that the two other methods. Again, it is important to note that both proximal algorithms yield sparse solutions, which is not the case for SG.
5.3 Denoising with Tree-Structured Wavelets
We demonstrate in this section how a tree-structured sparse regularization can improve classical wavelet representation, and how our method can be used to efficiently solve the corresponding large-scale optimization problems. We consider two wavelet orthonormal bases, Haar and Daubechies3 (see Mallat, 1999), and choose a classical quad-tree structure on the coefficients, which has notably proven to be useful for image compression problems (Baraniuk, 1999). This experiment follows the approach of Zhao et al. (2009) who used the same tree-structured regularization in the case of small one-dimensional signals, and the approach of Baraniuk et al. (2010) and Huang et al. (2009) images where images were reconstructed from compressed sensing measurements with a hierarchical nonconvex penalty.
We compare the performance for image denoising of both nonconvex and convex approaches. Specifically, we consider the following formulation
where is one of the orthonormal wavelet basis mentioned above, is the input noisy image, is the estimate of the denoised image, and is a sparsity-inducing regularization. Note that in this case, . We first consider classical settings where is either the -norm— this leads to the wavelet soft-thresholding method of Donoho and Johnstone (1995)— or the -pseudo-norm, whose solution can be obtained by hard-thresholding (see Mallat, 1999). Then, we consider the convex tree-structured regularization defined as a sum of -norms (-norms), which we denote by (respectively ). Since the basis is here orthonormal, solving the corresponding decomposition problems amounts to computing a single instance of the proximal operator. As a result, when is , we use Algorithm 3 and for , Algorithm 2 is applied. Finally, we consider the nonconvex tree-structured regularization used by Baraniuk et al. (2010) denoted here by , which we have presented in Eq. (4); the implementation details for can be found in Appendix A.
nonconvex method (the standard deviations are computed over the 12 images). CPU times (in second) averaged over all images and noise realizations are reported in brackets next to the names of the methods they correspond to.
Compared to Zhao et al. (2009), the novelty of our approach is essentially to be able to solve efficiently and exactly large-scale instances of this problem. We use classical standard test images,171717These images are used in classical image denoising benchmarks. See Mairal et al. (2009b).
and generate noisy versions of them corrupted by a white Gaussian noise of variance. For each image, we test several values of , with taken in a specific range.181818For the convex formulations, ranges in , while in the nonconvex case ranges in . We then keep the parameter giving the best reconstruction error. The factor
is a classical heuristic for choosing a reasonable regularization parameter(see Mallat, 1999). We provide reconstruction results in terms of PSNR in Table 1.191919Denoting by MSE the mean-squared-error for images whose intensities are between and , the PSNR is defined as and is measured in dB. A gain of dB reduces the MSE by approximately . We report in this table the results when is chosen to be a sum of -norms or -norms with weights all equal to one. Each experiment was run times with different noise realizations. In every setting, we observe that the tree-structured norm significantly outperforms the -norm and the nonconvex approaches. We also present a visual comparison on two images on Figure 15, showing that the tree-structured norm reduces visual artefacts (these artefacts are better seen by zooming on a computer screen). The wavelet transforms in our experiments are computed with the matlabPyrTools software.202020http://www.cns.nyu.edu/~eero/steerpyr/.
This experiment does of course not provide state-of-the-art results for image denoising (see Mairal et al., 2009b, and references therein), but shows that the tree-structured regularization significantly improves the reconstruction quality for wavelets. In this experiment the convex setting and also outperforms the nonconvex one .212121 It is worth mentioning that comparing convex and nonconvex approaches for sparse regularization is a bit difficult. This conclusion holds for the classical formulation we have used, but might not hold in other settings such as Coifman and Donoho (1995). We also note that the speed of our approach makes it scalable to real-time applications. Solving the proximal problem for an image with pixels takes approximately seconds on a single core of a 3.07GHz CPU if is a sum of -norms, and seconds when it is a sum of -norms. By contrast, unstructured approaches have a speed-up factor of about 7-8 with respect to the tree-structured methods.
5.4 Dictionaries of Natural Image Patches
This experiment studies whether a hierarchical structure can help dictionaries for denoising natural image patches, and in which noise regime the potential gain is significant. We aim at reconstructing corrupted patches from a test set, after having learned dictionaries on a training set of non-corrupted patches. Though not typical in machine learning, this setting is reasonable in the context of images, where lots of non-corrupted patches are easily available.222222Note that we study the ability of the model to reconstruct independent patches, and additional work is required to apply our framework to a full image processing task, where patches usually overlap (Elad and Aharon, 2006; Mairal et al., 2009b).
|noise||50 %||60 %||70 %||80 %||90 %|
We extracted patches of size pixels from the Berkeley segmentation database of natural images (Martin et al., 2001), which contains a high variability of scenes. We then split this dataset into a training set , a validation set , and a test set , respectively of size , , and patches. All the patches are centered and normalized to have unit -norm.
For the first experiment, the dictionary is learned on using the formulation of Eq. (10), with for as defined in Eq. (11). The validation and test sets are corrupted by removing a certain percentage of pixels, the task being to reconstruct the missing pixels from the known pixels. We thus introduce for each element of the validation/test set, a vector , equal to for the known pixel values and otherwise. Similarly, we define as the matrix equal to , except for the rows corresponding to missing pixel values, which are set to . By decomposing on , we obtain a sparse code , and the estimate of the reconstructed patch is defined as . Note that this procedure assumes that we know which pixel is missing and which is not for every element .
The parameters of the experiment are the regularization parameter used during the training step, the regularization parameter used during the validation/test step, and the structure of the tree. For every reported result, these parameters were selected by taking the ones offering the best performance on the validation set, before reporting any result from the test set. The values for the regularization parameters were selected on a logarithmic scale , and then further refined on a finer logarithmic scale with multiplicative increments of . For simplicity, we chose arbitrarily to use the -norm in the structured norm , with all the weights equal to one. We tested balanced tree structures of depth and , with different branching factors , where is the depth of the tree and , is the number of children for the nodes at depth . The branching factors tested for the trees of depth where , , and for trees of depth , , and , giving possible structures associated with dictionaries with at most elements. For each tree structure, we evaluated the performance obtained with the tree-structured dictionary along with a non-structured dictionary containing the same number of elements. These experiments were carried out four times, each time with a different initialization, and with a different noise realization.
Quantitative results are reported in Table 2. For all fractions of missing pixels considered, the tree-structured dictionary outperforms the “unstructured one”, and the most significant improvement is obtained in the noisiest setting. Note that having more dictionary elements is worthwhile when using the tree structure. To study the influence of the chosen structure, we report in Figure 16 the results obtained with the tested structures of depth , along with those obtained with unstructured dictionaries containing the same number of elements, when of the pixels are missing. For each dictionary size, the tree-structured dictionary significantly outperforms the unstructured one. An example of a learned tree-structured dictionary is presented on Figure 17. Dictionary elements naturally organize in groups of patches, often with low frequencies near the root of the tree, and high frequencies near the leaves.