Proximal Methods for Hierarchical Sparse Coding

09/11/2010 ∙ by Rodolphe Jenatton, et al. ∙ Inria 0

Sparse coding consists in representing signals as sparse linear combinations of atoms selected from a dictionary. We consider an extension of this framework where the atoms are further assumed to be embedded in a tree. This is achieved using a recently introduced tree-structured sparse regularization norm, which has proven useful in several applications. This norm leads to regularized problems that are difficult to optimize, and we propose in this paper efficient algorithms for solving them. More precisely, we show that the proximal operator associated with this norm is computable exactly via a dual approach that can be viewed as the composition of elementary proximal operators. Our procedure has a complexity linear, or close to linear, in the number of atoms, and allows the use of accelerated gradient techniques to solve the tree-structured sparse approximation problem at the same computational cost as traditional ones using the L1-norm. Our method is efficient and scales gracefully to millions of variables, which we illustrate in two types of applications: first, we consider fixed hierarchical dictionaries of wavelets to denoise natural images. Then, we apply our optimization tools in the context of dictionary learning, where learned dictionary elements naturally organize in a prespecified arborescent structure, leading to a better performance in reconstruction of natural image patches. When applied to text documents, our method learns hierarchies of topics, thus providing a competitive alternative to probabilistic topic models.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 20

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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 as

sparse 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)

. Hierarchies of latent variables, typically used in neural networks and deep learning architectures

(see Bengio, 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).

1.1 Notation

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 .

Figure 1: Example of a tree when . With the rule we consider for the nonzero patterns, if we have , we must also have for in .

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

(1)

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

(2)

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

We now turn to a convex reformulation of the constraint (1), which is the starting point for the convex optimization tools we develop in Section 3.

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:

(3)

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 :

(4)

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.

Figure 2: Left: example of a tree-structured set of groups (dashed contours in red), corresponding to a tree with nodes represented by black circles. Right: example of a sparsity pattern induced by the tree-structured norm corresponding to : the groups and are set to zero, so that the corresponding nodes (in gray) that form subtrees of are removed. The remaining nonzero variables form a rooted and connected subtree of . This sparsity pattern obeys the following equivalent rules: (i) if a node is selected, the same goes for all its ancestors. (ii) if a node is not selected, then its descendant are not selected.

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,

(5)

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 Section 

3 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

(6)

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.

3 Optimization

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: [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

(7)

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

(8)

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

(9)

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 .

  Inputs: and set of groups .
  Outputs: (primal-dual solutions).
  Initialization: .
  while  ( maximum number of iterations not reached )  do
     for   do
        .
     end for
  end while
  .
Algorithm 1 Block coordinate ascent in the dual

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.

  Inputs: and an ordered tree-structured set of groups .
  Outputs: (primal solution).
  Initialization: .
  for  , following the order ,  do
     .
  end for
Algorithm 2 Practical Computation of the Proximal Operator for - or -norms.

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 .

Thus, Algorithm 2 in fact performs a sequence of proximal operators, and we have shown the following corollary of Proposition 3.3: [Composition of Proximal Operators] 
Let such that . The proximal operator associated with the norm can be written as the composition of elementary operators:

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 .

0:   (input vector), set of groups , (positive weights), and (root of the tree).
1:  Variables: in (scaling factors); in (output, primal variable).
2:  computeSqNorm().
3:  recursiveScaling(,).
4:  Return (primal solution).

Procedure computeSqNorm()

1:  Compute the squared norm of the group: .
2:  Compute the scaling factor of the group: .
3:  Return .

Procedure recursiveScaling(,)

1:  .
2:  .
3:  for  do
4:     recursiveScaling(,).
5:  end for
Algorithm 3 Fast computation of the Proximal operator for -norm case.

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:

(10)

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

(11)

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.

5 Experiments

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

.

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.

(a) scale: small, regul.: low
(b) scale: small, regul.: medium
(c) scale: small, regul.: high
Figure 6: Benchmark for solving a least-squares regression problem regularized by the hierarchical norm . The experiment is small scale, , and shows the performances of six optimization methods (see main text for details) for three levels of regularization. The curves represent the relative value of the objective to the optimal value as a function of the computational time in second on a scale. All reported results are obtained by averaging runs.

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.

(a) scale: large, regul.: low
(b) scale: large, regul.: medium
(c) scale: large, regul.: high
Figure 10: Benchmark for solving a large-scale multi-class classification problem for four optimization methods (see details about the datasets and the methods in the main text). Three levels of regularization are considered. The curves represent the relative value of the objective to the optimal value as a function of the computational time in second on a scale. In the highly regularized setting, tuning the step-size for the subgradient turned out to be difficult, which explains the behavior of SG in the first iterations.

5.2.2 Multi-class classification of cancer diagnosis

The second benchmark explores a different supervised learning setting, where

is 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 number

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

Haar
PSNR 34.48 34.78 35.52 35.89 35.79
29.63 30.24 30.74 31.40 31.23
24.44 25.27 25.30 26.41 26.14
21.53 22.37 20.42 23.41 23.05
19.27 20.09 19.43 20.97 20.58
IPSNR -
-
-
-
-
Daub3
PSNR 34.64 34.95 35.74 36.14 36.00
30.03 30.63 31.10 31.79 31.56
25.04 25.84 25.76 26.90 26.54
22.09 22.90 22.42 23.90 23.41
19.56 20.45 19.67 21.40 20.87
IPSNR -
-
-
-
-
Table 1: Top part of the tables: Average PSNR measured for the denoising of standard images, when the wavelets are Haar or Daubechies3 wavelets (see Mallat, 1999), for two nonconvex approaches ( and ) and three different convex regularizations—that is, the -norm, the tree-structured sum of -norms (), and the tree-structured sum of -norms (). Best results for each level of noise and each wavelet type are in bold. Bottom part of the tables: Average improvement in PSNR with respect to the

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

(a) Lena, ,
(b) Lena, ,
(c) Barb., ,
(d) Barb., ,
Figure 15: Visual comparison between the wavelet shrinkage model with the -norm and the tree-structured model, on cropped versions of the images Lena and Barb.. Haar wavelets are used.

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 %
 flat
 tree
Table 2: Quantitative results of the reconstruction task on natural image patches. First row: percentage of missing pixels. Second and third rows: mean square error multiplied by , respectively for classical sparse coding, and tree-structured sparse coding.
Figure 16: Mean square error multiplied by obtained with structures with error bars, sorted by number of dictionary elements from to . Red plain bars represents the tree-structured dictionaries. White bars correspond to the flat dictionary model containing the same number of dictionary as the tree-structured one. For readability purpose, the -axis of the graph starts at .

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.

Figure 17: Learned dictionary with a tree structure of depth . The root of the tree is in the middle of the figure. The branching factors are , , , . The dictionary is learned on patches of size pixels.