Proximal Methods for Hierarchical Sparse Coding

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 $\mathcal{T}$, and the sparsity patterns are constrained to form a connected and rooted subtree of $\mathcal{T}$ 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.A 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.

In this formulation, the sparsity-inducing norm $\Omega$ encodes a hierarchical structure among the atoms of ${\mathbf{D}}$, where this structure is assumed to be known beforehand. The precise meaning of hierarchical structure and the definition of $\Omega$ 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 ${\mathbf{D}}$ 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 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).

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.

Problem Statement and Related Work

In the rest of the paper, we focus on specific sets of nonzero coefficients—or simply, nonzero patterns—for the decomposition vector ${\boldsymbol{\alpha}}$. In particular, we assume that we are given a treeOur analysis straightforwardly extends to the case of a forest of trees; for simplicity, we consider a single tree ${\mathcal{T}}$. ${\mathcal{T}}$ whose $p$ nodes are indexed by $j$ in $\{1,\dots,p\}$. We want the nonzero patterns of ${\boldsymbol{\alpha}}$ to form a connected and rooted subtree of ${\mathcal{T}}$; in other words, if $\text{ancestors}(j)\subseteq\{1,\dots,p\}$ denotes the set of indices corresponding to the ancestorsWe consider that the set of ancestors of a node also contains the node itself. of the node $j$ in ${\mathcal{T}}$ (see Figure 1), the vector ${\boldsymbol{\alpha}}$ obeys the following rule

Informally, we want to exploit the structure of ${\mathcal{T}}$ in the following sense: the decomposition of any signal ${\mathbf{x}}$ can involve a dictionary element ${\mathbf{d}}^{j}$ only if the ancestors of ${\mathbf{d}}^{j}$ in the tree ${\mathcal{T}}$ 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.

For a given sparsity level $s\geq 0$ (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 $\lambda\|{\boldsymbol{\alpha}}\|_{0}$ to the objective function in place of the constraint $\|{\boldsymbol{\alpha}}\|_{0}\leq s$) 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 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.

Condition (1) can be equivalently expressed by its contrapositive, thus leading to an intuitive way of penalizing the vector ${\boldsymbol{\alpha}}$ to obtain tree-structured nonzero patterns. More precisely, defining $\text{descendants}(j)\subseteq\{1,\dots,p\}$ analogously to $\text{ancestors}(j)$ for $j$ in $\{1,\dots,p\}$, 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 $\mathcal{G}$ the set defined by $\mathcal{G}\stackrel{{\scriptstyle\vartriangle}}{{=}}\{\text{descendants}(j);j\in\{1,\dots,p\}\},$ and refer to each member $g$ of $\mathcal{G}$ as a group (Figure 2). To obtain a decomposition with the desired property (3), one can naturally penalize the number of groups $g$ in $\mathcal{G}$ that are “involved” in the decomposition of ${\mathbf{x}}$, i.e., that record at least one nonzero coefficient of ${\boldsymbol{\alpha}}$:

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:

A set of groups $\mathcal{G}\!\stackrel{{\scriptstyle\vartriangle}}{{=}}\!\{g\}_{g\in\mathcal{G}}$ is said to be tree-structured in $\{1,\dots,p\}$, if $\,\bigcup_{g\in\mathcal{G}}\!g=\{1,\dots,p\}$ and if for all $g,h\in\mathcal{G}$, $(g\cap h\neq\emptyset)\Rightarrow(g\subseteq h\ \text{or}\ h\subseteq g).$ For such a set of groups, there exists a (non-unique) total order relation $\preceq$ such that:

Given such a tree-structured set of groups $\mathcal{G}$ and its associated norm $\Omega$, we are interested throughout the paper in the following hierarchical sparse coding problem,

Since $\mathcal{P}$ is a partition, the set of groups in $\mathcal{P}$ and the singletons form together a tree-structured set of groups according to definition 1 and the algorithm we will develop is therefore applicable to this problem.

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 $\Omega$, several dedicated procedures have been proposed.

Optimization

We begin with a brief introduction to proximal methods, necessary to present our contributions. From now on, we assume that $f$ 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 $f$. 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 $f$ (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.

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 $f$ around the current estimate $\hat{{\boldsymbol{\alpha}}}$, 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 $f$ is close to its linear approximation, and $L\!>\!0$ is a parameter which is an upper bound on the Lipschitz constant of $\nabla f$. 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 $\Omega$ is not present, the previous proximal problem exactly leads to the standard gradient update rule. More generally, we define the proximal operator:

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,

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 $\|.\|$:

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 ${\boldsymbol{\xi}}^{g}$, while keeping fixed the other variables ${\boldsymbol{\xi}}^{h}$, for $h\neq g$. It is easy to see from Eq. (8) that such an update of a column ${\boldsymbol{\xi}}^{g}$, for a group $g$ in $\mathcal{G}$, amounts to computing the orthogonal projection of the vector ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-\sum_{h\neq g}\xi_{{{\scriptscriptstyle\mid}g}}^{h}$ onto the ball of radius $\lambda\omega_{g}$ of the dual norm $\|.\|_{\ast}$.

3 Convergence in One Pass

Before stating this result, we need to introduce a lemma showing that, given two nested groups $g,h$ such that $g\subseteq h\subseteq\{1,\dots,p\}$, if ${\boldsymbol{\xi}}^{g}$ is updated before ${\boldsymbol{\xi}}^{h}$ in Algorithm 1, then the optimality condition for ${\boldsymbol{\xi}}^{g}$ is not perturbed by the update of ${\boldsymbol{\xi}}^{h}$.

with $t_{g},t_{h}>0$. Let us introduce ${\mathbf{v}}={\mathbf{u}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}^{h}$. The following relationships hold

The previous lemma establishes the convergence in one pass of Algorithm 1 in the case where $\mathcal{G}$ only contains two nested groups $g\subseteq h$, provided that ${\boldsymbol{\xi}}^{g}$ is computed before ${\boldsymbol{\xi}}^{h}$. Let us illustrate this fact more concretely. After initializing ${\boldsymbol{\xi}}^{g}$ and ${\boldsymbol{\xi}}^{h}$ to zero, Algorithm 1 first updates ${\boldsymbol{\xi}}^{g}$ with the formula ${\boldsymbol{\xi}}^{g}\leftarrow\Pi_{\|.\|_{\ast}\leq\lambda\omega_{g}}({\mathbf{u}}_{{{\scriptscriptstyle\mid}g}})$, and then performs the following update: ${\boldsymbol{\xi}}^{h}\leftarrow\Pi_{\|.\|_{\ast}\leq\lambda\omega_{h}}({\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}-{\boldsymbol{\xi}}^{g})$ (where we have used that ${\boldsymbol{\xi}}^{g}={\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}h}}^{g}$ since $g\subseteq h$). We are now in position to apply Lemma 2 which states that the primal/dual variables $\{{\mathbf{v}},{\boldsymbol{\xi}}^{g},{\boldsymbol{\xi}}^{h}\}$ satisfy the optimality conditions (9), as described in Lemma 1. In only one pass over the groups $\{g,h\}$, 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 $\mathcal{G}$:

Proof The proof largely relies on Lemma 2 and proceeds by induction. By definition of Algorithm 1, the feasibility of ${\boldsymbol{\xi}}$ is always guaranteed. We consider the following induction hypothesis

Since the dual variables ${\boldsymbol{\xi}}$ are initially equal to zero, the summation over $g^{\prime}\preceq h,\ g^{\prime}\neq g$ is equivalent to a summation over $g^{\prime}\neq g$. We initialize the induction with the first group in $\mathcal{G}$, that, by definition of $\preceq$, does not contain any other group. The first step of Algorithm 1 easily shows that the induction hypothesis $\mathcal{H}$ is satisfied for this first group.

We now assume that $\mathcal{H}(h)$ is true and consider the next group $h^{\prime}$, $h\preceq h^{\prime}$, in order to prove that $\mathcal{H}(h^{\prime})$ is also satisfied. We have for each group $g\subseteq h$,

Since ${\boldsymbol{\xi}}^{g}_{{{\scriptscriptstyle\mid}h}^{\prime}}={\boldsymbol{\xi}}^{g}$ for $g\subseteq h^{\prime}$, we have

and following the update rule for the group $h^{\prime}$,

4 Interpretation in Terms of Composition of Proximal Operators

In Algorithm 1, since all the vectors ${\boldsymbol{\xi}}^{g}$ are initialized to ${\mathbf{0}}$, when the group $g$ is considered, we have by induction ${\mathbf{u}}-\sum_{h\neq g}{\boldsymbol{\xi}}^{h}={\mathbf{u}}-\sum_{h\preceq g}{\boldsymbol{\xi}}^{h}$. Thus, to maintain at each iteration of the inner loop ${\mathbf{v}}={\mathbf{u}}-\sum_{h\neq g}{\boldsymbol{\xi}}^{h}$ one can instead update ${\mathbf{v}}$ after updating ${\boldsymbol{\xi}}^{g}$ according to ${\mathbf{v}}\leftarrow{\mathbf{v}}-\,{\boldsymbol{\xi}}^{g}$. Moreover, since ${\boldsymbol{\xi}}^{g}$ is no longer needed in the algorithm, and since only the entries of ${\mathbf{v}}$ indexed by $g$ are updated, we can combine the two updates into ${\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}\leftarrow{\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}-\Pi_{\|.\|_{\ast}\leq\lambda\omega_{g}}({\mathbf{v}}_{{{\scriptscriptstyle\mid}g}})$, leading to a simplified Algorithm 2 equivalent to Algorithm 1.

Thus, Algorithm 2 in fact performs a sequence of $|\mathcal{G}|$ proximal operators, and we have shown the following corollary of Proposition 1:

Let $g_{1}\preccurlyeq\ldots\preccurlyeq g_{m}$ such that $\mathcal{G}=\{g_{1},\ldots,g_{m}\}$. The proximal operator $\text{Prox}_{\lambda\Omega}$ associated with the norm $\Omega$ can be written as the composition of elementary operators:

5 Efficient Implementation and Complexity

Algorithm 2 gives the solution of the primal problem Eq. (7) in $O(pd)$ operations, where $d$ is the depth of the tree.

Lemma 3 should not suggest that the complexity is linear in $p$, since $d$ could depend of $p$ as well, and in the worst case the hierarchy is a chain, yielding $d=p-1$. However, in a balanced tree, $d=O(\log(p))$. In practice, the structures we have considered experimentally are relatively flat, with a depth not exceeding $d=5$, and the complexity is therefore almost linear.

To formulate the algorithm, two new notations are used: for a group $g$ in $\mathcal{G}$, we denote by $\text{root}(g)$ the indices of the variables that are at the root of the subtree corresponding to $g$,As a reminder, $\text{root}(g)$ is not a singleton when several dictionary elements are considered per node. and by $\text{children}(g)$ the set of groups that are the children of $\text{root}(g)$ in the tree. For example, in the tree presented in Figure 2, $\text{root}(\{3,5,6\})\!=\!\{3\}$, $\text{root}(\{1,2,3,4,5,6\})\!=\!\{1\}$, $\text{children}(\{3,5,6\})\!=\!\{\{5\},\{6\}\}$, and $\text{children}(\{1,2,3,4,5,6\})\!=\!\{\{2,4\},\{3,5,6\}\}$. Note that all the groups of $\text{children}(g)$ are necessarily included in $g$.

So far the dictionary ${\mathbf{D}}$ 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 ${\mathbf{D}}$ adapted to our hierarchical sparse coding formulation.

Application to Dictionary Learning

We start by briefly describing dictionary learning.

While learning simultaneously ${\mathbf{D}}$ and ${\mathbf{A}}$, 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 ${\mathbf{A}}$ Olshausen and Field (1997); Aharon et al. (2006); Lee et al. (2007); Mairal et al. (2010a). This leads to penalizing or constraining $({\mathbf{D}},{\mathbf{A}})$ and results in the following formulation:

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.

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 ${\mathbf{D}}$ and ${\mathbf{A}}=[{\boldsymbol{\alpha}}^{1},\ldots,{\boldsymbol{\alpha}}^{n}]$ while keeping the other variable fixed Aharon et al. (2006); Lee et al. (2007); Mairal et al. (2010a).Note 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 ${\boldsymbol{\alpha}}^{i}$, $i$ in $\{1,\ldots,n\}$, for the dictionary ${\mathbf{D}}$ kept fixed. Because of $\Omega$, the corresponding convex subproblem is nonsmooth and has to be solved for each of the $n$ signals considered. The optimization of the dictionary ${\mathbf{D}}$ (for ${\mathbf{A}}$ fixed), which we discuss first, is in general easier.

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 ${\mathbf{D}}$. Specifically, we assume that the domain set $\mathcal{D}$ has the form

Experiments

We next turn to the experimental validation of our hierarchical sparse coding.

In Section 3.3, we have shown that the proximal operator associated to $\Omega$ 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.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 $1.5$.

Finally, we monitor the convergence of the algorithm by checking the relative decrease in the cost function.We 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/.

2 Speed Benchmark

In this first benchmark, we consider a least-squares regression problem regularized by $\Omega$ 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 $16\!\times\!16$-patches. The dictionary we use is represented on Figure 7. Although the problem involves a small number of variables, i.e., $p=151$ 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.

2.2 Multi-class classification of cancer diagnosis

The dataset contains $m=308$ samples, $p=30\,017$ variables and $26$ classes. In addition, the data exhibit highly-correlated dictionary elements. Inspired by Kim and Xing (2010), we build the tree-structured set of groups $\mathcal{G}$ using Ward’s hierarchical clustering Johnson (1967) on the gene expressions. The norm $\Omega$ built in this way aims at capturing the hierarchical structure of gene expression networks Kim and Xing (2010).

The results in Figure 4 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.

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

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.Note 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).

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 6 the results obtained with the $13$ tested structures of depth $3$, along with those obtained with unstructured dictionaries containing the same number of elements, when $90\%$ 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 7. Dictionary elements naturally organize in groups of patches, often with low frequencies near the root of the tree, and high frequencies near the leaves.

5 Text Documents

This last experimental section shows that our approach can also be applied to model text corpora. The goal of probabilistic topic models is to find a low-dimensional representation of a collection of documents, where the representation should provide a semantic description of the collection. Approaching the problem in a parametric Bayesian framework, latent Dirichlet allocation (LDA) Blei et al. (2003) model documents, represented as vectors of word counts, as a mixture of a predefined number of latent topics that are distributions over a fixed vocabulary. LDA is fundamentally a matrix factorization problem: Buntine (2002) shows that LDA can be interpreted as a Dirichlet-multinomial counterpart of factor analysis. The number of topics is usually small compared to the size of the vocabulary (e.g., 100 against $10\,000$), so that the topic proportions of each document provide a compact representation of the corpus. For instance, these new features can be used to feed a classifier in a subsequent classification task. We similarly use our dictionary learning approach to find low-dimensional representations of text corpora.

Discussion

We have applied hierarchical sparse coding in various settings, with fixed/learned dictionaries, and based on different types of data, namely, natural images and text documents. A line of research to pursue is to develop other optimization tools for structured norms with general overlapping groups. For instance, Mairal et al. (2010b) have used network flow optimization techniques for that purpose, and Bach (2010) submodular function optimization. This framework can also be used in the context of hierarchical kernel learning Bach (2008), where we believe that our method can be more efficient than existing ones.

This work establishes a connection between dictionary learning and probabilistic topic models, which should prove fruitful as the two lines of work have focused on different aspects of the same unsupervised learning problem: Our approach is based on convex optimization tools, and provides experimentally more stable data representations. Moreover, it can be easily extended with the same tools to other types of structures corresponding to other norms Jenatton et al. (2009); Jacob et al. (2009). It should be noted, however, that, unlike some Bayesian methods, dictionary learning by itself does not provide mechanisms for the automatic selection of model hyper-parameters (such as the dictionary size or the topology of the tree). An interesting common line of research to pursue could be the supervised design of dictionaries, which has been proved useful in the two frameworks Mairal et al. (2009a); Bradley and Bagnell (2009); Blei and McAuliffe (2008).

This paper was partially supported by grants from the Agence Nationale de la Recherche (MGA Project) and from the European Research Council (SIERRA Project 239993). The authors would like to thank Jean Ponce for interesting discussions and suggestions for improving this manuscript. They also would like to thank Volkan Cevher for pointing out links between our approach and nonconvex tree-structured regularization and for insightful discussions. Finally, we thank the reviewers for their constructive and helpful comments.

A Links with Tree-Structured Nonconvex Regularization

We present in this section an algorithm introduced by Donoho (1997) in the more general context of approximation from dyadic partitions (see Section 6 in Donoho, 1997). This algorithm solves the following problem

We now briefly show how to derive the dynamic programming approach introduced by Donoho (1997). Given a group $g$ in $\mathcal{G}$, we use the same notations $\text{root}(g)$ and children(g) introduced in Section 3.5. It is relatively easy to show that finding a solution of Eq. (12) amounts to finding the support $S\subseteq\{1,\ldots,p\}$ of its solution and that the problem can be equivalently rewritten

with the abusive notation $\delta^{g}(S)=1$ if $g\cap S\neq\emptyset$ and otherwise. We now introduce the quantity

After a few computations, solving Eq. (13) can be shown to be equivalent to minimizing $\psi_{g_{0}}(S)$ where $g_{0}$ is the root of the tree. It is then easy to prove that for any group $g$ in $\mathcal{G}$, we have

which leads to the following dynamic programming approach presented in Algorithm 4.

This algorithm shares several conceptual links with Algorithm 2 and 3. It traverses the tree in the same order, has a complexity in $O(p)$, and it can be shown that the whole procedure actually performs a sequence of thresholding operations on the variable ${\mathbf{v}}$.

B Proofs

We gather here the proofs of the technical results of the paper.

This primal problem is convex and satisfies Slater’s conditions for generalized conic inequalities (i.e., existence of a feasible point in the interior of the domain), which implies that strong duality holds Boyd and Vandenberghe (2004). We now consider the Lagrangian $\mathcal{L}$ defined as

The dual function is obtained by minimizing out the primal variables. To this end, we take the derivatives of $\mathcal{L}$ with respect to the primal variables ${\mathbf{v}}$ and ${\mathbf{z}}$ and set them to zero, which leads to

After simplifying the Lagrangian and flipping (without loss of generality) the sign of ${\boldsymbol{\xi}}$, we obtain the dual problem in Eq. (8). We derive the optimality conditions from the Karush–Kuhn–Tucker conditions for generalized conic inequalities Boyd and Vandenberghe (2004). We have that $\{{\mathbf{v}},{\mathbf{z}},{\boldsymbol{\tau}},{\boldsymbol{\xi}}\}$ are optimal if and only if

Combining the complementary slackness with the definition of the dual norm, we have

Furthermore, using the fact that $\forall g\in\mathcal{G},\ ({\mathbf{v}}_{{{\scriptscriptstyle\mid}g}},z_{g})\in\mathcal{C}$ and $({\boldsymbol{\xi}}^{g},\tau_{g})=({\boldsymbol{\xi}}^{g},\lambda\omega_{g})\in\mathcal{C^{\ast}}$, we obtain the following chain of inequalities

for which equality must hold. In particular, we have ${\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}^{\top}{\boldsymbol{\xi}}^{g}=\|{\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}\|\|{\boldsymbol{\xi}}^{g}\|_{\ast}$ and $z_{g}\|{\boldsymbol{\xi}}^{g}\|_{\ast}=\lambda z_{g}\omega_{g}$. If ${\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}\neq 0$, then $z_{g}$ cannot be equal to zero, which implies in turn that $\|{\boldsymbol{\xi}}^{g}\|_{\ast}=\lambda\omega_{g}$. Eventually, applying Lemma 5 gives the advertised optimality conditions.

Conversely, starting from the optimality conditions of Lemma 1, and making use again of Lemma 5, we can derive the Karush–Kuhn–Tucker conditions displayed above. More precisely, we define for all $g\in\mathcal{G}$,

The only condition that needs to be discussed is the complementary slackness condition. If ${\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}=0$, then it is easily satisfied. Otherwise, combining the definitions of $\tau_{g},\ z_{g}$ and the fact that

we end up with the desired complementary slackness.

B.2 Optimality condition for the projection on the dual ball

Proof When the vector ${\mathbf{w}}$ is already in the ball of $\|.\|_{\ast}$ with radius $t$, i.e., $\|{\mathbf{w}}\|_{\ast}\leq t$, the situation is simple, since the projection $\Pi_{\|.\|_{\ast}\leq t}({\mathbf{w}})$ obviously gives ${\mathbf{w}}$ itself. On the other hand, a necessary and sufficient optimality condition for having ${\boldsymbol{\kappa}}=\Pi_{\|.\|_{\ast}\leq t}({\mathbf{w}})=\operatornamewithlimits{arg\,min}_{\|{\mathbf{y}}\|_{\ast}\leq t}\|{\mathbf{w}}-{\mathbf{y}}\|_{2}$ is that the residual ${\mathbf{w}}-{\boldsymbol{\kappa}}$ lies in the normal cone of the constraint set Borwein and Lewis (2006), that is, for all ${\mathbf{y}}$ such that $\|{\mathbf{y}}\|_{\ast}\!\leq t$, $({\mathbf{w}}-{\boldsymbol{\kappa}})^{\top}\!({\mathbf{y}}-{\boldsymbol{\kappa}})\!\leq 0$. The displayed result then follows from the definition of the dual norm, namely $\|{\boldsymbol{\kappa}}\|_{\ast}\!=\!\max_{\|{\mathbf{z}}\|\leq 1}{\mathbf{z}}^{\top}{\boldsymbol{\kappa}}$.

B.3 Proof of Lemma 2

Proof First, notice that the conclusion ${\boldsymbol{\xi}}^{h}=\Pi_{\|.\|_{\ast}\leq\lambda\omega_{h}}({\mathbf{v}}_{{{\scriptscriptstyle\mid}h}}+{\boldsymbol{\xi}}^{h})$ simply comes from the definition of ${\boldsymbol{\xi}}^{h}$ and ${\mathbf{v}}$, along with the fact that ${\boldsymbol{\xi}}^{g}={\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}h}}^{g}$ since $g\subseteq h$. We now examine ${\boldsymbol{\xi}}^{g}$.

The proof mostly relies on the optimality conditions characterizing the projection onto a ball of the dual norm $\|\cdot\|_{\ast}$. Precisely, by Lemma 5, we need to show that either

Note that the feasibility of ${\boldsymbol{\xi}}^{g}$, i.e., $\|{\boldsymbol{\xi}}^{g}\|_{\ast}\leq t_{g}$, holds by definition of $\kappa^{g}$.

Let us first assume that $\|{\boldsymbol{\xi}}^{g}\|_{\ast}<t_{g}$. We necessarily have that ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}$ also lies in the interior of the ball of $\|.\|_{\ast}$ with radius $t_{g}$, and it holds that ${\boldsymbol{\xi}}^{g}={\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}$. Since $g\subseteq h$, we have that the vector ${\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}-{\boldsymbol{\xi}}^{g}={\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}-{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}$ has only zero entries on $g$. As a result, ${\boldsymbol{\xi}}^{h}_{g}=0$ (or equivalently, ${\boldsymbol{\xi}}^{h}_{{\scriptscriptstyle\mid}g}=0$) and we obtain

which is the desired conclusion. From now on, we assume that $\|{\boldsymbol{\xi}}^{g}\|_{\ast}=t_{g}$. It then remains to show that

We now distinguish two cases, according to the norm used.

Note that the case $\rho_{h}=0$ leads to ${\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}^{h}=0$, and therefore ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}=0$ since $g\subseteq h$, which directly yields the result. The case $\rho_{g}=0$ implies ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}=0$ and therefore ${\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}=0$, yielding the result as well. Now, we can therefore assume $\rho_{h}>0$ and $\rho_{g}>0$. From the first equality of (14), we have ${\boldsymbol{\xi}}^{g}={\boldsymbol{\xi}}^{g}_{{\scriptscriptstyle\mid}g}$ since $(\rho_{g}+1){\boldsymbol{\xi}}^{g}={\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}$. Further using the fact that $g\subseteq h$ in the second equality of (14), we obtain

This implies that ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}=\rho_{g}{\boldsymbol{\xi}}^{g}-\frac{\rho_{g}}{\rho_{h}+1}{\boldsymbol{\xi}}^{g}$, which eventually leads to

The desired conclusion follows ${\boldsymbol{\xi}}^{g\top}({\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h})=\|{\boldsymbol{\xi}}^{g}\|_{2}\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{2}.$

Looking at the same condition for ${\boldsymbol{\xi}}^{h}$, we have that {\boldsymbol{\xi}}^{h}=\Pi_{\|.\|_{\ast}\leq t_{h}}\big{(}{\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}-{\boldsymbol{\xi}}_{g}\big{)} holds if and only if for all ${\boldsymbol{\xi}}^{h}_{j}\neq 0,j\in h$, we have

From those relationships we notably deduce that for all $j\in g$ such that ${\boldsymbol{\xi}}^{g}_{j}\neq 0$, $\operatorname{sign}({\boldsymbol{\xi}}^{g}_{j})=\operatorname{sign}({\mathbf{u}}_{j})=\operatorname{sign}({\boldsymbol{\xi}}^{h}_{j})=\operatorname{sign}({\mathbf{u}}_{j}-{\boldsymbol{\xi}}^{g}_{j})=\operatorname{sign}({\mathbf{u}}_{j}-{\boldsymbol{\xi}}^{g}_{j}-{\boldsymbol{\xi}}^{h}_{j})$. Let $j\in g$ such that ${\boldsymbol{\xi}}^{g}_{j}\neq 0$. At this point, using the equalities we have just presented,

Since $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}\|_{\infty}\geq\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{\infty}$ (which can be shown using the sign equalities above), and $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}^{h}\|_{\infty}\geq\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{\infty}$ (since $g\subseteq h$), we have

and therefore for all ${\boldsymbol{\xi}}^{g}_{j}\neq 0$, $j\in g$, we have ${\mathbf{u}}_{j}-{\boldsymbol{\xi}}^{g}_{j}-{\boldsymbol{\xi}}^{h}_{j}=\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}^{g}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{\infty}\operatorname{sign}({\boldsymbol{\xi}}^{g}_{j}),$ which yields the result.

B.4 Proof of Lemma 4

Proof Notice first that the procedure computeSqNorm is called exactly once for each group $g$ in $\mathcal{G}$, computing a set of scalars $(\rho_{g})_{g\in\mathcal{G}}$ in an order which is compatible with the convergence in one pass of Algorithm 1—that is, the children of a node are processed prior to the node itself. Following such an order, the update of the group $g$ in the original Algorithm 1 computes the variable ${\boldsymbol{\xi}}^{g}$ which updates implicitly the primal variable as follows

It is now possible to show by induction that for all group $g$ in $\mathcal{G}$, after a call to the procedure computeSqNorm($g$), the auxiliary variable $\eta_{g}$ takes the value $\|{\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}\|_{2}^{2}$ where ${\mathbf{v}}$ has the same value as during the iteration $g$ of Algorithm 1. Therefore, after calling the procedure computeSqNorm($g_{0}$), where $g_{0}$ is the root of the tree, the values $\rho_{g}$ correspond to the successive scaling factors of the variable ${\mathbf{v}}_{{{\scriptscriptstyle\mid}g}}$ obtained during the execution of Algorithm 1. After having computed all the scaling factors $\rho_{g}$, $g\in\mathcal{G}$, the procedure recursiveScaling ensures that each variable $j$ in $\{1,\ldots,p\}$ is scaled by the product of all the $\rho_{h}$, where $h$ is an ancestor of the variable $j$.

The complexity of the algorithm is easy to characterize: Each procedure computeSqNorm and recursiveScaling is called $p$ times, each call for a group $g$ has a constant number of operations plus as many operations as the number of children of $p$. Since each child can be called at most one time, the total number of operation of the algorithm is $O(p)$.

B.5 Sign conservation by projection

Proof Let us consider ${\boldsymbol{\kappa}}=\Pi_{\|.\|_{\ast}\leq t}({\mathbf{w}})$. Using essentially the same argument as in the proof of Lemma 5, we have for all ${\mathbf{y}}$ such that $\|{\mathbf{y}}\|_{q^{\prime}}\!\leq t$, $({\mathbf{w}}-{\boldsymbol{\kappa}})^{\top}\!({\mathbf{y}}-{\boldsymbol{\kappa}})\!\leq 0$. Noticing that ${\mathbf{S}}^{\top}{\mathbf{S}}={\mathbf{I}}$ and $\|{\mathbf{y}}\|_{q^{\prime}}=\|{\mathbf{S}}{\mathbf{y}}\|_{q^{\prime}}$, we further obtain $({\mathbf{S}}{\mathbf{w}}-{\mathbf{S}}{\boldsymbol{\kappa}})^{\top}\!({\mathbf{y}}^{\prime}-{\mathbf{S}}{\boldsymbol{\kappa}})\!\leq 0$ for all ${\mathbf{y}}^{\prime}$ with $\|{\mathbf{y}}^{\prime}\|_{q^{\prime}}\!\leq t$. This implies in turn that ${\mathbf{S}}\Pi_{\|.\|_{\ast}\leq t}({\mathbf{w}})=\Pi_{\|.\|_{\ast}\leq t}({\mathbf{S}}{\mathbf{w}})$, which is equivalent to the advertised conclusion. Based on this lemma, note that we can assume without loss of generality that the vector we want to project (in this case, ${\mathbf{w}}$) has only nonnegative entries. Indeed, it is sufficient to store beforehand the signs of that vector, compute the projection of the vector with nonnegative entries, and assign the stored signs to the result of the projection.

B.6 Non-negativity constraint for the proximal operator

which shows that $\hat{{\mathbf{z}}}^{+}$ is the solution of the right-hand side of (15).

We show in this section that there exists a setting for which the conclusion of Lemma 2 does not hold anymore. We first focus on a necessary condition of Lemma 2:

If the conclusion of Lemma 2 holds—that is, we have ${\boldsymbol{\xi}}^{g}=\Pi_{\|.\|_{\ast}\leq t_{g}}({\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h})$, notice that it is not possible to have the following scenarios, as proved below by contradiction:

If $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{q^{\prime}}<t_{g}$, then we would have ${\boldsymbol{\xi}}^{g}={\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}$, which is impossible since $\|{\boldsymbol{\xi}}^{g}\|_{q^{\prime}}=t_{g}$.

If $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{q^{\prime}}=t_{g}$, then we would have for all $j$ in $g$, $|{\boldsymbol{\xi}}_{j}^{h}|^{q^{\prime}}=\rho_{h}|{\mathbf{u}}_{j}-{\boldsymbol{\xi}}_{j}^{g}-{\boldsymbol{\xi}}_{j}^{h}|^{q}=0$, which implies that ${\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}=0$ and $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}\|_{q^{\prime}}=t_{g}$. This is impossible since we assumed $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}\|_{q^{\prime}}>t_{g}$.

We therefore have $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}\|_{q^{\prime}}>t_{g}$ and using again the second optimality conditions of Lemma 5, there exists $\rho>0$ such that for all $j$ in $g$, $|{\boldsymbol{\xi}}_{j}^{g}|^{q^{\prime}}=\rho|{\mathbf{u}}_{j}-{\boldsymbol{\xi}}_{j}^{g}-{\boldsymbol{\xi}}_{j}^{h}|^{q}$. Combined with the previous relation on ${\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}$, we obtain for all $j$ in $g$, $|{\boldsymbol{\xi}}_{j}^{g}|^{q^{\prime}}=\frac{\rho}{\rho_{h}}|{\boldsymbol{\xi}}_{j}^{h}|^{q^{\prime}}.$ Since we can assume without loss of generality that ${\mathbf{u}}$ only has nonnegative entries (see Lemma 6), the vectors ${\boldsymbol{\xi}}^{g}$ and ${\boldsymbol{\xi}}^{h}$ can also be assumed to have nonnegative entries, hence the desired conclusion. We need another intuitive property of the projection $\Pi_{\|.\|_{\ast}\leq t}$ to derive our counterexample:

Proof Let us first notice that given the assumption on $t$, we have $\|{\boldsymbol{\kappa}}\|_{q^{\prime}}=t$. The Lagrangian $\mathcal{L}$ associated with the convex minimization problem underlying the definition of $\Pi_{\|.\|_{\ast}\leq t}$ can be written as

At optimality, the stationarity condition for ${\boldsymbol{\kappa}}$ leads to

We can assume without loss of generality that ${\mathbf{w}}$ only has nonnegative entries (see Lemma 6). Since the components of ${\boldsymbol{\kappa}}$ and ${\mathbf{w}}$ have the same signs (see Lemma 6), we therefore have $|{\boldsymbol{\kappa}}_{j}|={\boldsymbol{\kappa}}_{j}\geq 0$, for all $j$ in $\{1,\dots,p\}$. Note that $\alpha$ cannot be equal to zero because of $\|{\boldsymbol{\kappa}}\|_{q^{\prime}}=t<\|{\mathbf{w}}\|_{q^{\prime}}$.

Let us consider the continuously differentiable function $\varphi_{w}:\kappa\mapsto\kappa-w+\alpha q^{\prime}\kappa^{q^{\prime}-1}$ defined on $(0,\infty)$. Since $\varphi_{w}(0)=-w<0$, $\lim_{\kappa\to\infty}\varphi_{w}(\kappa)=\infty$ and $\varphi_{w}$ is strictly nondecreasing, there exists a unique $\kappa^{\ast}_{w}>0$ such that $\varphi_{w}(\kappa^{\ast}_{w})=0$. If we now take $w<v$, we have

With $\varphi_{v}$ being strictly nondecreasing, we thus obtain $\kappa^{\ast}_{w}<\kappa^{\ast}_{v}$. The desired conclusion stems from the application of the previous result to the stationarity condition of ${\boldsymbol{\kappa}}$.

Based on the two previous lemmas, we are now in position to present our counterexample:

with $t_{g},t_{h}>0$ satisfying $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}\|_{q^{\prime}}>t_{g}$ and $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}h}}\|_{q^{\prime}}>t_{g}+t_{h}$. Then, the conclusion of Lemma 2 does not hold.

Proof We apply the same rationale as in the proof of Lemma 9. Writing the stationarity conditions for ${\boldsymbol{\xi}}^{g}$ and ${\boldsymbol{\xi}}^{h}$, we have for all $j$ in $g$

with Lagrangian parameters $\alpha,\beta>0$. We now proceed by contradiction and assume that ${\boldsymbol{\xi}}^{g}=\Pi_{\|.\|_{\ast}\leq t_{g}}({\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}-{\boldsymbol{\xi}}_{{{\scriptscriptstyle\mid}g}}^{h}).$ According to Lemma 8, there exists $\rho>0$ such that for all $j$ in $g$, ${\boldsymbol{\xi}}_{j}^{h}=\rho{\boldsymbol{\xi}}_{j}^{g}.$ If we combine the previous relations on ${\boldsymbol{\xi}}^{g}$ and ${\boldsymbol{\xi}}^{h}$, we obtain for all $j$ in $g$,

If $C<0$, then we have a contradiction, since the entries of ${\boldsymbol{\xi}}^{g}$ and ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}$ have the same signs. Similarly, the case $C=0$ leads a contradiction, since we would have ${\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}=0$ and $\|{\mathbf{u}}_{{{\scriptscriptstyle\mid}g}}\|_{q^{\prime}}>t_{g}$. As a consequence, it follows that $C>0$ and for all $j$ in $g$, {\boldsymbol{\xi}}_{j}^{g}=\exp\big{\{}\frac{\log(C)}{2-q^{\prime}}\big{\}}, which means that all the entries of the vector ${\boldsymbol{\xi}}_{g}^{g}$ are identical. Using Lemma 9, since there exists $(i,j)\in g\times g$ such that ${\mathbf{u}}_{i}<{\mathbf{u}}_{j}$, we also have ${\boldsymbol{\xi}}_{i}^{g}<{\boldsymbol{\xi}}_{j}^{g}$, which leads to a contradiction.

References