Structured Variable Selection with Sparsity-Inducing Norms

Introduction

These real world examples motivate the need for the design of sparsity-inducing regularization schemes, capable of encoding more sophisticated prior knowledge about the expected sparsity patterns.

In this paper, we consider all possible sets of groups and characterize exactly what type of prior knowledge can be encoded by considering sums of norms of overlapping groups of variables. Before describing how to go from groups to nonzero patterns (or equivalently zero patterns), we show that it is possible to “reverse-engineer” a given set of nonzero patterns, i.e., to build the unique minimal set of groups that will generate these patterns. This allows the automatic design of sparsity-inducing norms, adapted to target sparsity patterns. We give in Section 3 some interesting examples of such designs in specific geometric and structured configurations, which covers the type of prior knowledge available in the real world applications described previously.

As will be shown in Section 3, for each set of groups, a notion of hull of a nonzero pattern may be naturally defined. For example, in the particular case of the two-dimensional planar grid considered in this paper, this hull is exactly the axis-aligned bounding box or the regular convex hull. We show that, in our framework, the allowed nonzero patterns are exactly those equal to their hull, and that the hull of the relevant variables is consistently estimated under certain conditions, both in low and high-dimensional settings. Moreover, we present in Section 4 an efficient active set algorithm that scales well to high dimensions. Finally, we illustrate in Section 6 the behavior of our norms with synthetic examples on specific geometric settings, such as lines and two-dimensional grids.

Regularized Risk Minimization

We focus on a general family of sparsity-inducing norms that allow the penalization of subsets of variables grouped together. Let us denote by $\mathcal{G}$ a subset of the power set of $\{1,\dots,p\}$ such that $\bigcup_{G\in\mathcal{G}}\!G=\{1,\dots,p\}$, i.e., a spanning set of subsets of $\{1,\dots,p\}$. Note that $\mathcal{G}$ does not necessarily define a partition of $\{1,\dots,p\}$, and therefore, it is possible for elements of $\mathcal{G}$ to overlap. We consider the norm $\Omega$ defined by

This general formulation has several important sub-cases that we present below, the goal of this paper being to go beyond these, and to consider norms capable to incorporate richer prior knowledge.

Hierarchical norms: when the set $\{1,\dots,p\}$ is embedded into a tree (Zhao et al., 2009) or more generally into a directed acyclic graph (Bach, 2009b), then a set of $p$ groups, each of them composed of descendants of a given variable, is considered.

We study the following regularized problem:

Groups and Sparsity Patterns

We now study the relationship between the norm $\Omega$ defined in Eq. (1) and the nonzero patterns the estimated vector $\hat{w}$ is allowed to have. We first characterize the set of nonzero patterns, then we provide forward and backward procedures to go back and forth from groups to patterns.

We can equivalently use $\mathcal{P}$ or $\mathcal{Z}$ by taking the complement of each element of these sets.

Let $Q$ denote the Gram matrix $\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\top}$. We consider the optimization problem in Eq. (2) with $\mu>0$. If $Q$ is invertible or if $\{1,\dots,p\}$ belongs to $\mathcal{G}$, then the problem in Eq. (2) admits a unique solution.

In other words, when $Y=(y_{1},\dots,y_{n})^{\top}$ is a realization of an absolutely continuous probability distribution, the sparse solutions have a zero pattern in $\mathcal{Z}=\left\{\bigcup_{G\in\mathcal{G}^{\prime}}G;\ \mathcal{G}^{\prime}\subseteq\mathcal{G}\right\}$ almost surely. As a corollary of our two results, if the Gram matrix $Q$ is invertible, the problem in Eq. (2) has a unique solution, whose zero pattern belongs to $\mathcal{Z}$ almost surely. Note that with the assumption made on $Y$, Theorem 2 is not directly applicable to the classification setting. Based on these previous results, we can look at the following usual special cases from Section 2 (we give more examples in Section 3.5):

Hierarchical norms: the set of patterns $\mathcal{P}$ is then all sets $J$ for which all ancestors of elements in $J$ are included in $J$ (Bach, 2009b).

Two natural questions now arise: (1) starting from the groups $\mathcal{G}$, is there an efficient way to generate the set of nonzero patterns $\mathcal{P}$; (2) conversely, and more importantly, given $\mathcal{P}$, how can the groups $\mathcal{G}$—and hence the norm $\Omega(w)$—be designed?

2 General Properties of 𝒢𝒢\mathcal{G}, 𝒵𝒵\mathcal{Z} and 𝒫𝒫\mathcal{P}

We now study the different properties of the set of groups $\mathcal{G}$ and its corresponding sets of patterns $\mathcal{Z}$ and $\mathcal{P}$.

Minimality.

If a group in $\mathcal{G}$ is the union of other groups, it may be removed from $\mathcal{G}$ without changing the sets $\mathcal{Z}$ or $\mathcal{P}$. This is the main argument behind the pruning backward algorithm in Section 3.3. Moreover, this leads to the notion of a minimal set $\mathcal{G}$ of groups, which is such that for all $\mathcal{G}^{\prime}\subseteq\mathcal{Z}$ whose union-closure spans $\mathcal{Z}$, we have $\mathcal{G}\subseteq\mathcal{G}^{\prime}$. The existence and uniqueness of a minimal set is a consequence of classical results in set theory (Doignon and Falmagne, 1998). The elements of this minimal set are usually referred to as the atoms of $\mathcal{Z}$.

Minimal sets of groups are attractive in our setting because they lead to a smaller number of groups and lower computational complexity—for example, for 2 dimensional-grids with rectangular patterns, we have a quadratic possible number of rectangles, i.e., $|\mathcal{Z}|=O(p^{2})$, that can be generated by a minimal set $\mathcal{G}$ whose size is $|\mathcal{G}|=O(\sqrt{p})$.

Hull.

Given a set of groups $\mathcal{G}$, we can define for any subset $I\subseteq\{1,\dots,p\}$ the $\mathcal{G}$-adapted hull, or simply hull, as:

Graphs of patterns.

We consider the directed acyclic graph (DAG) stemming from the Hasse diagram (Cameron, 1994) of the partially ordered set (poset) $(\mathcal{G},\supset)$. By definition, the nodes of this graph are the elements $G$ of $\mathcal{G}$ and there is a directed edge from $G_{1}$ to $G_{2}$ if and only if $G_{1}\supset G_{2}$ and there exists no $G\in\mathcal{G}$ such that $G_{1}\supset G\supset G_{2}$ (Cameron, 1994). We can also build the corresponding DAG for the set of zero patterns $\mathcal{Z}\supset\mathcal{G}$, which is a super-DAG of the DAG of groups (see Figure 3 for examples). Note that we obtain also the isomorphic DAG for the nonzero patterns $\mathcal{P}$, although it corresponds to the poset $(\mathcal{P},\subset)$: this DAG will be used in the active set algorithm presented in Section 4.

Prior works with nested groups (Zhao et al., 2009; Bach, 2009b; Kim and Xing, 2009) have also used a similar DAG structure, with the slight difference that in these works, the corresponding hierarchy of variables is built from the prior knowledge about the problem at hand (e.g., the tree of wavelets in Zhao et al. (2009), the decomposition of kernels in Bach (2009b) or the hierarchy of genes in Kim and Xing (2009)). The DAG we introduce here on the set of groups naturally and always comes up, with no assumption on the variables themselves (for which no DAG is defined in general).

3 From Patterns to Groups

4 From Groups to Patterns

5 Examples

We now present several examples of sets of groups $\mathcal{G}$, especially suited to encode geometric and temporal prior information.

Given $p$ variables organized in a sequence, if we want only contiguous nonzero patterns, the backward algorithm will lead to the set of groups which are intervals $[1,k]_{k\in\{1,\dots,p-1\}}$ and $[k,p]_{k\in\{2,\dots,p\}}$, with both $|\mathcal{Z}|=O(p^{2})$ and $|\mathcal{G}|=O(p)$ (see Figure 4). Imposing the contiguity of the nonzero patterns is for instance relevant for the diagnosis of tumors, based on the profiles of arrayCGH (Rapaport et al., 2008).

Two-dimensional grids.

In Section 6, we notably consider for $\mathcal{P}$ the set of all rectangles in two dimensions, leading by the previous algorithm to the set of axis-aligned half-spaces for $\mathcal{G}$ (see Figure 5), with $|\mathcal{Z}|=O(p^{2})$ and $|\mathcal{G}|=O(\sqrt{p})$. This type of structure is encountered in object or scene recognition, where the selected rectangle would correspond to a certain box inside an image, that concentrates the predictive power for a given class of object/scene (Harzallah et al., 2009).

Extensions.

The sets of groups presented above can be straightforwardly extended to more complicated topologies, such as three-dimensional spaces discretized in cubes or spherical volumes discretized in slices. Similar properties hold for such settings. For instance, if all the axis-aligned half-spaces are considered for $\mathcal{G}$ in a three-dimensional space, then $\mathcal{P}$ is the set of all possible rectangular boxes with $|\mathcal{P}|=O(p^{2})$ and $|\mathcal{G}|=O(p^{1/3})$. Such three-dimensional structures may be interesting to retrieve discriminative and local sets of voxels from fMRI/MEEG responses (Gramfort and Kowalski, 2009; Xiang et al., 2009). Moreover, while the two-dimensional rectangular patterns described previously are adapted to find bounding boxes in static images (Harzallah et al., 2009), scene recognition in videos requires to deal with a third temporal dimension (Dalal et al., 2006). This may be achieved by designing appropriate sets of groups, embedded in the three-dimensional space obtained by tracking the frames over time.

Representation and computation of 𝒢𝒢\mathcal{G}.

The sets of groups described so far can actually be represented in a same form, that lends itself well to the analysis of the next section. When dealing with a discrete sequence of length $l$ (see Figure 4), we have

with $g_{-}^{k}=\{i;\;1\leq i\leq k\}$ and $g_{+}^{k}=\{i;\;l\geq i\geq k\}$. In other words, the set of groups $\mathcal{G}$ can be rewritten as a partitionNote the subtlety: the sets $\mathcal{G}_{\theta}$ are disjoint, that is $\mathcal{G}_{\theta}\cap\mathcal{G}_{\theta^{\prime}}=\varnothing$ for $\theta\neq\theta^{\prime}$, but groups in $\mathcal{G}_{\theta}$ and $\mathcal{G}_{\theta^{\prime}}$ can overlap. in two sets of nested groups, $\mathcal{G}_{\text{left}}$ and $\mathcal{G}_{\text{right}}$.

The same goes for a two-dimensional grid, with dimensions $h\!\times\!l$ (see Figure 5 and Figure 6). In this case, the nested groups we consider are defined based on the following groups of variables

Optimization and Active Set Algorithm

For moderate values of $p$, one may obtain a solution for Eq. (2) using generic toolboxes for second-order cone programming (SOCP) whose time complexity is equal to $O(p^{3.5}+|\mathcal{G}|^{3.5})$ (Boyd and Vandenberghe, 2004), which is not appropriate when $p$ or $|\mathcal{G}|$ are large. This time complexity corresponds to the computation of Eq. (2) for a single value of the regularization parameter $\mu$.

We present in this section an active set algorithm (Algorithm 3) that finds a solution for Eq. (2) by considering increasingly larger active sets and checking global optimality at each step. When the rectangular groups are used, the total complexity of this method is in $O(s\max\{p^{1.75}\!,s^{3.5}\})$, where $s$ is the size of the active set at the end of the optimization. Here, the sparsity prior is exploited for computational advantages. Our active set algorithm needs an underlying black-box SOCP solver; in this paper, we consider both a first order approach (see Appendix H) and a SOCP methodThe C$++$/Matlab code used in the experiments may be downloaded from the authors website. — in our experiments, we use SDPT3 (Toh et al., 1999; Tütüncü et al., 2003). Our active set algorithm extends to general overlapping groups the work of Bach (2009b), by further assuming that it is computationally possible to have a time complexity polynomial in the number of variables $p$.

We primarily focus here on finding an efficient active set algorithm; we defer to future work the design of specific SOCP solvers, e.g., based on proximal techniques (see, e.g., Tseng, 2009, and numerous references therein), adapted to such non-smooth sparsity-inducing penalties.

Let $J\subseteq\{1,\dots,p\}$. The following two problems

are dual to each other and strong duality holds. The pair of primal-dual variables $\{w_{J},\kappa_{J}\}$ is optimal if and only if we have

As a brief reminder, the duality gap of a minimization problem is defined as the difference between the primal and dual objective functions, evaluated for a feasible pair of primal/dual variables (Boyd and Vandenberghe, 2004, see Section 5.5). This gap serves as a certificate of (sub)optimality: if it is equal to zero, then the optimum is reached, and provided that strong duality holds, the converse is true as well (Boyd and Vandenberghe, 2004, see Section 5.5).

The previous proposition enables us to derive the duality gap for the optimization problem Eq. (6), that is reduced to the active set of variables $J$. In practice, this duality gap will always vanish (up to the precision of the underlying SOCP solver), since we will sequentially solve Eq. (6) for increasingly larger active sets $J$. We now study how, starting from the optimality of the problem in Eq. (6), we can control the optimality, or equivalently the duality gap, for the full problem Eq. (5). More precisely, the duality gap of the optimization problem Eq. (6) is

which is a sum of two nonnegative terms, the nonnegativity coming from the Fenchel-Young inequality (Borwein and Lewis, 2006; Boyd and Vandenberghe, 2004, Proposition 3.3.4 and Section 3.3.2 respectively). We can think of this duality gap as the sum of two duality gaps, respectively relative to $L_{J}$ and $\Omega_{J}$. Thus, if we have a primal candidate $w_{J}$ and we choose $\kappa_{J}=-\nabla L_{J}(w_{J})$, the duality gap relative to $L_{J}$ vanishes and the total duality gap then reduces to

In order to check that the reduced solution $w_{J}$ is optimal for the full problem in Eq. (5), we pad $w_{J}$ with zeros on ${J^{c}}$ to define $w$ and compute $\kappa=-\nabla L(w)$, which is such that $\kappa_{J}=-\nabla L_{J}(w_{J})$. For our given candidate pair of primal/dual variables $\{w,\kappa\}$, we then get a duality gap for the full problem in Eq. (5) equal to

Computing this gap requires computing the dual norm which itself is as hard as the original problem, prompting the need for upper and lower bounds on $\Omega^{*}$ (see Propositions 4 and 5 for more details).

2 Active set algorithm

We can interpret the active set algorithm as a walk through the DAG of nonzero patterns allowed by the norm $\Omega$. The parents $\Pi_{\mathcal{P}}(J)$ of $J$ in this DAG are exactly the patterns containing the variables that may enter the active set at the next iteration of Algorithm 3. The groups that are exactly at the boundaries of the active set (referred to as the fringe groups) are $\mathcal{F}_{J}=\{G\in(\mathcal{G}_{J})^{c}\ ;\ \nexists G^{\prime}\in(\mathcal{G}_{J})^{c},\ G\subseteq G^{\prime}\}$, i.e., the groups that are not contained by any other inactive groups.

In simple settings, e.g., when $\mathcal{G}$ is the set of rectangular groups, the correspondence between groups and variables is straightforward since we have $\mathcal{F}_{J}=\bigcup_{K\in\Pi_{\mathcal{P}}(J)}\mathcal{G}_{K}\backslash\mathcal{G}_{J}$ (see Figure 7). However, in general, we just have the inclusion $(\bigcup_{K\in\Pi_{\mathcal{P}}(J)}\mathcal{G}_{K}\backslash\mathcal{G}_{J})\subseteq\mathcal{F}_{J}$ and some elements of $\mathcal{F}_{J}$ might not correspond to any patterns of variables in $\Pi_{\mathcal{P}}(J)$ (see Figure 8).

We now present the optimality conditions (see proofs in Appendix E) that monitor the progress of Algorithm 3:

If $w$ is optimal for the full problem in Eq. (5), then

then $w$ is an approximate solution for Eq. (5) whose duality gap is less than $\varepsilon\geq 0$.

Note that for the Lasso, the conditions $(N)$ and $(S_{0})$ (i.e., the sufficient condition taken with $\varepsilon=0$) are both equivalent (up to the squaring of $\Omega$) to the condition $\|\nabla L(w)_{J^{c}}\|_{\infty}\leq-w^{\top}\nabla L(w)$, which is the usual optimality condition (Wainwright, 2009; Tibshirani, 1996). Moreover, when they are not satisfied, our two conditions provide good heuristics for choosing which $K\in\Pi_{\mathcal{P}}(J)$ should enter the active set.

More precisely, since the necessary condition ($N$) directly deals with the variables (as opposed to groups) that can become active at the next step of Algorithm 3, it suffices to choose the pattern $K\in\Pi_{\mathcal{P}}(J)$ that violates most the condition.

The heuristics for the sufficient condition ($S_{\varepsilon}$) implies to go from groups to variables. We simply consider the group $G\in\mathcal{F}_{J}$ that violates most the sufficient condition and then take all the patterns of variables $K\in\Pi_{\mathcal{P}}(J)$ such that $\ K\cap G\neq\varnothing$ to enter the active set. If $G\cap(\bigcup_{K\in\Pi_{\mathcal{P}}(J)}K)=\varnothing$, we look at all the groups $H\in\mathcal{F}_{J}$ such that $H\cap G\neq\varnothing$ and apply the scheme described before (see Algorithm 4).

A direct consequence of this heuristics is that it is possible for the algorithm to jump over the right active set and to consider instead a (slightly) larger active set as optimal. However if the active set is larger than the optimal set, then (it can be proved that) the sufficient condition $(S_{0})$ is satisfied, and the reduced problem, which we solve exactly, will still output the correct nonzero pattern.

Moreover, it is worthwhile to notice that in Algorithm 3, the active set may sometimes be increased only to make sure that the current solution is optimal (we only check a sufficient condition of optimality).

The procedure described in Algorithm 3 can terminate in two different states. If the procedure stops because of the limit on the number of active variables $s$, the solution might be suboptimal. Note that, in any case, we have at our disposal a upperbound on the duality gap.

Otherwise, the procedure always converges to an optimal solution, either (1) by validating both the necessary and sufficient conditions (see Propositions 4 and 5), ending up with fewer than $p$ active variables and a precision of (at least) $\varepsilon$, or (2) by running until the $p$ variables become active, the precision of the solution being given by the underlying solver.

Algorithmic complexity.

We analyze in detail the time complexity of the active set algorithm when we consider sets of groups $\mathcal{G}$ such as those presented in the examples of Section 3.5. We recall that we denote by $\Theta$ the set of orientations in $\mathcal{G}$ (for more details, see Section 3.5).

For such choices of $\mathcal{G}$, the fringe groups $\mathcal{F}_{J}$ reduces to the largest groups of each orientation and therefore $|\mathcal{F}_{J}|\leq|\Theta|$. We further assume that the groups in $\mathcal{G}_{\theta}$ are sorted by cardinality, so that computing $\mathcal{F}_{J}$ costs $O(|\Theta|)$.

Given an active set $J$, both the necessary and sufficient conditions require to have access to the direct parents $\Pi_{\mathcal{P}}(J)$ of $J$ in the DAG of nonzero patterns. In simple settings, e.g., when $\mathcal{G}$ is the set of rectangular groups, this operation can be performed in $O(1)$ (it just corresponds to scan the (up to) four patterns at the edges of the current rectangular hull).

However, for more general orientations, computing $\Pi_{\mathcal{P}}(J)$ requires to find the smallest nonzero patterns that we can generate from the groups in $\mathcal{F}_{J}$, reduced to the stripe of variables around the current hull. This stripe of variables can be computed as \big{[}\bigcup_{G\in(\mathcal{G}_{J})^{c}\backslash\mathcal{F}_{J}}G\big{]}^{c}\backslash J, so that getting $\Pi_{\mathcal{P}}(J)$ costs $O(s2^{|\Theta|}+p|\mathcal{G}|)$ in total.

Thus, if the number of active variables is upper bounded by $s\!\ll\!p$ (which is true if our target is actually sparse), the time complexity of Algorithm 3 is the sum of:

the computation of the gradient, $O(snp)$ for the square loss.

if the underlying solver called upon by the active set algorithm is a standard SOCP solver, $O(s\max_{J\in\mathcal{P},|J|\leq s}|\mathcal{G}_{J}|^{3.5}+s^{4.5})$ (note that the term $s^{4.5}$ could be improved upon by using warm-restart strategies for the sequence of reduced problems).

$t_{1}$ times the computation of $(N)$, that is $O(t_{1}(s2^{|\Theta|}+p|\mathcal{G}|+sn_{\theta}^{2})+p|\mathcal{G}|)=O(t_{1}p|\mathcal{G}|)$.

During the initialization (i.e., $J=\varnothing$), we have $|\Pi_{\mathcal{P}}(\varnothing)|=O(p)$ (since we can start with any singletons), and $|\mathcal{G}_{K}\backslash\mathcal{G}_{J}|=|\mathcal{G}_{K}|=|\mathcal{G}|$, which leads to a complexity of $O(p|\mathcal{G}|)$ for the sum $\sum_{G\in\mathcal{G}_{K}\backslash\mathcal{G}_{J}}=\sum_{G\in\mathcal{G}_{K}}$. Note however that this sum does not depend on $J$ and can therefore be cached if we need to make several runs with the same set of groups $\mathcal{G}$.

$t_{2}$ times the computation of $(S_{\varepsilon})$, that is $O(t_{2}(s2^{|\Theta|}+p|\mathcal{G}|+|\Theta|^{2}+|\Theta|p+p|\mathcal{G}|))=O(t_{2}p|\mathcal{G}|)$, with $t_{1}+t_{2}\leq s$.

We finally get complexity with a leading term in $O(sp|\mathcal{G}|+s\max_{J\in\mathcal{P},|J|\leq s}|\mathcal{G}_{J}|^{3.5}+s^{4.5})$, which is much better than $O(p^{3.5}+|\mathcal{G}|^{3.5})$, without an active set method. In the example of the two-dimensional grid (see Section 3.5), we have $|\mathcal{G}|=O(\sqrt{p})$ and $O(s\max\{p^{1.75}\!,s^{3.5}\})$ as total complexity. The simulations of Section 6 confirm that the active set strategy is indeed useful when $s$ is much smaller than $p$. Moreover, the two extreme cases where $s\approx p$ or $p\ll 1$ are also shown not to be advantageous for the active set strategy, since either it is cheaper to use the SOCP solver directly on the $p$ variables, or we uselessly pay the additional fixed-cost of the active set machinery (such as computing the optimality conditions). Note that we have derived here the theoretical complexity of the active set algorithm when we use a SOCP method as underlying solver. With the first order method presented in Appendix H, we would instead get a total complexity in $O(sp^{1.5})$.

3 Intersecting Nonzero Patterns

We have seen so far how overlapping groups can encore prior information about a desired set of (non)zero patterns. In practice, controlling these overlaps may be delicate and hinges on the choice of the weights $(d^{\scriptscriptstyle G})_{G\in\mathcal{G}}$ (see the experiments in Section 6). In particular, the weights have to take into account that some variables belonging to several overlapping groups are penalized multiple times.

However, it is possible to keep the benefit of overlapping groups whilst limiting their side effects, by taking up the idea of support intersection (Bach, 2008a; Meinshausen and Bühlmann, 2008). First introduced to stabilize the set of variables recovered by the Lasso, we reuse this technique in a different context, based on the fact that $\mathcal{Z}$ is closed under union.

If we deal with the same sets of groups as those considered in Section 3.5, it is natural to rewrite $\mathcal{G}$ as $\bigcup_{\theta\in\Theta}\!\mathcal{G}_{\theta}$, where $\Theta$ is the set of the orientations of the groups in $\mathcal{G}$ (for more details, see Section 3.5). Let us denote by $\hat{w}$ and \mbox{\hat{w}}^{\theta} the solutions of Eq. (5), where the regularization term $\Omega$ is respectively defined by the groups in $\mathcal{G}$ and by the groupsTo be more precise, in order to regularize every variable, we add the full group $\{1,\dots,p\}$ to $\mathcal{G}_{\theta}$, which does not modify $\mathcal{P}$. in $\mathcal{G}_{\theta}$.

The main point is that, since $\mathcal{P}$ is closed under intersection, the two procedures described below actually lead to the same set of allowed nonzero patterns:

Simply considering the nonzero pattern of $\hat{w}$.

Taking the intersection of the nonzero patterns obtained for each \mbox{\hat{w}}^{\theta}, $\theta$ in $\Theta$.

With the latter procedure, although the learning of several models \mbox{\hat{w}}^{\theta} is required (a number of times equals to the number of orientations considered, e.g., 2 for the sequence, 4 for the rectangular groups and more generally $|\Theta|$ times), each of those learnings involves a smaller number of groups (that is, just the ones belonging to $\mathcal{G}_{\theta}$). In addition, this procedure is a variable selection technique that therefore needs a second step for estimating the loadings (restricted to the selected nonzero pattern). In the experiments, we follow Bach (2008a) and we use an ordinary least squares (OLS). The simulations of Section 6 will show the benefits of this variable selection approach.

Pattern Consistency

In this section, we analyze the model consistency of the solution of Eq. (2) for the square loss. Considering the set of nonzero patterns $\mathcal{P}$ derived in Section 3, we can only hope to estimate the correct hull of the generating sparsity pattern, since Theorem 2 states that other patterns occur with zero probability. We derive necessary and sufficient conditions for model consistency in a low-dimensional setting, and then consider a high-dimensional result.

We begin with the low-dimensional setting where $n$ is tending to infinity with $p$ fixed. In addition, we also assume that the design is fixed and that the Gram matrix $Q=\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\top}$ is invertible with positive-definite (i.e., invertible) limit

In this setting, the noise is the only source of randomness. We denote by $\mathbf{r}_{\mathbf{J}}$ the vector defined as

In the Lasso and group Lasso setting, the vector $\mathbf{r}_{\mathbf{J}}$ is respectively the sign vector \textrm{sign}(\mbox{{\mathbf{w}}}_{\mathbf{J}}) and the vector defined by the blocks (\frac{\mbox{{\mathbf{w}}}_{G}}{\left\|\mbox{{\mathbf{w}}}_{G}\right\|_{2}})_{G\in\mathcal{G}_{\mathbf{J}}}.

We define $\Omega_{\mathbf{J}}^{c}(w_{\mathbf{J}^{c}})=\sum_{G\in(\mathcal{G}_{\mathbf{J}})^{c}}\left\|d^{\scriptscriptstyle G}_{\mathbf{J}^{c}}\circ w_{\mathbf{J}^{c}}\right\|_{2}$ (which is the norm composed of inactive groups) with its dual norm $(\Omega_{\mathbf{J}}^{c})^{\ast}$; note the difference with the norm reduced to $\mathbf{J}^{c}$, defined as $\Omega_{\mathbf{J}^{c}}(w_{\mathbf{J}^{c}})=\sum_{G\in\mathcal{G}}\left\|d^{\scriptscriptstyle G}_{\mathbf{J}^{c}}\circ w_{\mathbf{J}^{c}}\right\|_{2}$.

The following Theorem gives the sufficient and necessary conditions under which the hull of the generating pattern is consistently estimated. Those conditions naturally extend the results of Zhao and Yu (2006) and Bach (2008b) for the Lasso and the group Lasso respectively (see proof in Appendix F).

Assume $\mu\rightarrow 0,\ \mu\sqrt{n}\rightarrow\infty$ in Eq. (2). If the hull is consistently estimated, then $(\Omega_{\mathbf{J}}^{c})^{\ast}[\mathbf{Q}_{\mathbf{J}^{c}\mathbf{J}}\mathbf{Q}_{\mathbf{J}\mathbf{J}}^{-1}\mathbf{r}_{\mathbf{J}}]\leq 1$. Conversely, if $(\Omega_{\mathbf{J}}^{c})^{\ast}[\mathbf{Q}_{\mathbf{J}^{c}\mathbf{J}}\mathbf{Q}_{\mathbf{J}\mathbf{J}}^{-1}\mathbf{r}_{\mathbf{J}}]<1$, then the hull is consistently estimated, i.e.,

The two previous propositions bring into play the dual norm $(\Omega_{\mathbf{J}}^{c})^{\ast}$ that we cannot compute in closed form, but requires to solve an optimization problem as complex as the initial problem in Eq. (5). However, we can prove bounds similar to those obtained in Propositions 4 and 5 for the necessary and sufficient conditions.

2 High-Dimensional Analysis

We prove a high-dimensional variable consistency result (see proof in Appendix G) that extends the corresponding result for the Lasso (Zhao and Yu, 2006; Wainwright, 2009), by assuming that the consistency condition in Theorem 6 is satisfied.

Assume that $Q$ has unit diagonal, $\kappa=\lambda_{\min}(Q_{\mathbf{J}\mathbf{J}})>0$ and $(\Omega_{\mathbf{J}}^{c})^{\ast}[Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1}\mathbf{r}_{\mathbf{J}}]<1-\tau$, with $\tau>0$. If $\tau\mu\sqrt{n}\geq\sigma C_{3}(\mathcal{G},\mathbf{J}),$ and $\mu|\mathbf{J}|^{1/2}\leq C_{4}(\mathcal{G},\mathbf{J}),$ then the probability of incorrect hull selection is upper bounded by:

where $C_{1}(\mathcal{G},\mathbf{J})$, $C_{2}(\mathcal{G},\mathbf{J})$, $C_{3}(\mathcal{G},\mathbf{J})$ and $C_{4}(\mathcal{G},\mathbf{J})$ are constants defined in Appendix G, which essentially depend on the groups, the smallest nonzero coefficient of ${\mathbf{w}}$ and how close the support \{j\in\mathbf{J}:\mbox{{\mathbf{w}}}_{j}\neq 0\} of ${\mathbf{w}}$ is to its hull $\mathbf{J}$, that is the relevance of the prior information encoded by $\mathcal{G}$.

In the Lasso case, we have $C_{1}(\mathcal{G},\mathbf{J})=O(1)$, $C_{2}(\mathcal{G},\mathbf{J})=O(|\mathbf{J}|^{-2})$, $C_{3}(\mathcal{G},\mathbf{J})=O((\log p)^{1/2})$ and $C_{4}(\mathcal{G},\mathbf{J})=O(|\mathbf{J}|^{-1})$, leading to the usual scaling $n\approx\log p$ and $\mu\approx\sigma(\log p/n)^{1/2}$.

We can also give the scaling of these constants in simple settings where groups overlap. For instance, let us consider that the variables are organized in a sequence (see Figure 4). Let us further assume that the weights $(d^{\scriptscriptstyle G})_{G\in\mathcal{G}}$ satisfy the following two properties:

The weights take into account the overlaps, that is,

with $t\mapsto\beta(t)\in(0,1]$ a non-increasing function such that $\beta(0)=1$,

is upper bounded by a constant $\mathcal{K}$ independent of $p$.

Note that we consider such weights in the experiments (see Section 6). Based on these assumptions, some algebra directly leads to

We thus obtain a scaling similar to the Lasso (with, in addition, a control of the allowed nonzero patterns).

With stronger assumptions on the possible positions of $\mathbf{J}$, we may have better scalings, but these are problem-dependent (a careful analysis of the group-dependent constants would still be needed in all cases).

Experiments

In this section, we carry out several experiments to illustrate the behavior of the sparsity-inducing norm $\Omega$. We denote by Structured-lasso, or simply Slasso, the models regularized by the norm $\Omega$. In addition, the procedure (introduced in Section 4.3) that consists in intersecting the nonzero patterns obtained for different models of Slasso will be referred to as Intersected Structured-lasso, or simply ISlasso.

uniform weights, $d^{\scriptscriptstyle G}_{j}=1,$

weights depending on the size of the groups, $d^{\scriptscriptstyle G}_{j}=|G|^{-2},$

weights that take into account overlapping groups, $d^{\scriptscriptstyle G}_{j}=\rho^{\,|\{H\in\mathcal{G}\,;\,H\ni j,\ H\subset G\mbox{ and }H\neq G\}|},$ for some $\rho\in(0,1)$.

For each orientation in $\mathcal{G}$, the third type of weights (W3) aims at reducing the unbalance caused by the overlapping groups. Specifically, given a group $G\in\mathcal{G}$ and a variable $j\in G$, the corresponding weight $d^{\scriptscriptstyle G}_{j}$ is all the more small as the variable $j$ already belongs to other groups with the same orientation.

Unless otherwise specified, we use the third type of weights (W3) with $\rho=0.5$. In the following experiments, the loadings $w_{\mathbf{J}}$, as well as the design matrices, are generated from a standard Gaussian distribution with identity covariance matrix. The positions of $\mathbf{J}$ are also random and are uniformly drawn.

We show in this experiment that the prior information we put through the norm $\Omega$ improves upon the predictive power. We are looking at two situations where we can express a structural prior through $\Omega$, namely (1) the selection of a contiguous pattern on a sequence and (2) the selection of a convex pattern on a grid (see Figure 9).

In what follows, we consider $p=400$ variables with generating patterns ${\mathbf{w}}$ whose hulls have a constant size of $|\mathbf{J}|=24$ variables. In order to evaluate the relevance of the contiguous (or convex) prior, we also vary the number of zero variables that are contained in the hull (see Figure 9). We then compute the prediction error for different sample sizes $n\in\{250,500,1000\}$. The prediction error is understood here as

where $\hat{w}$ denotes the estimate of the OLS, performed on the nonzero pattern found by the model considered (i.e., either Lasso, Slasso or ISlasso)Even though we make comparisons based on prediction errors, the experiments illustrate the results from Section 5 since we first use our method as a variable selection step.. The regularization parameter is chosen by 5-fold cross-validation and the test set consists of 500 samples. For each value of $n$, we display on Figure 11 and Figure 12 the median errors over 50 random settings \{\mathbf{J},\mbox{{\mathbf{w}}},X,\varepsilon\}, for respectively the sequence and the grid.

First and foremost, the simulations highlight how important the weights $(d^{\scriptscriptstyle G})_{G\in\mathcal{G}}$ are. In particular, the uniform (W1) and size-dependent weights (W2) perform poorly since they do not take into account the overlapping groups. The models learned with such weights do not manage to recover the correct nonzero patterns (and even worse, they tend to select every variable—see the right column of Figure 11).

Although groups that moderately overlap do help (e.g., see the Slasso with the weights (W3) on the left column of Figure 11), it remains delicate to handle groups with many overlaps, even with an appropriate choice of $(d^{\scriptscriptstyle G})_{G\in\mathcal{G}}$ (e.g., see the right column of Figure 12 where Slasso considers up to 8 overlaps on the grid). The ISlasso procedure does not suffer from this issue since it reduces the number of overlaps whilst keeping the desirable effects of overlapping groups. Another way to yield a better level of sparsity, even with many overlaps, would be to consider non-convex alternatives to $\Omega$ (see, e.g., Jenatton et al., 2010). Naturally, the benefit of ISlasso is more significant on the grid than on the sequence as the latter deals with fewer overlaps. Moreover, adding the $\pm\pi/4$-groups to the rectangular groups enables to recover a nonzero pattern closer to the generating pattern. This is illustrated on the left column of Figure 12 where the error of ISlasso with only rectangular groups (in black) corresponds to the selection of the smallest rectangular box around the generating pattern.

On the other hand, and more importantly, the experiments show that if the prior about the generating pattern is relevant, then our structured approach performs better that the standard Lasso. Indeed, as displayed on the left columns of Figure 11 and Figure 12, as soon as the hull of the generating pattern does not contain too many zero variables, Slasso/ISlasso outperform Lasso. In fact, the sample complexity of the Lasso depends on the number of nonzero variables in ${\mathbf{w}}$ as opposed to the size of the hull for Slasso/ISlasso. This also explains why the error for Slasso/ISlasso is almost constant with respect to the number of nonzero variables (since the hull has a constant size). Note finally that, even though our structured approach does not always dramatically outperform the standard unstructured approach in terms of prediction, it has the advantage of being more interpretable.

Active set algorithm.

We finally focus on the active set algorithm (see Section 4) and compare its time complexity to the SOCP solver when we are looking for a sparse structured target. More precisely, for a fixed level of sparsity $|\mathbf{J}|=24$ and a fixed number of observations $n=3500$, we analyze the complexity with respect to the number of variables $p$ that varies in $\{100,225,400,900,1600,2500\}$.

We consider the same experimental protocol as above except that we display the median CPU time based onlyNote that it already corresponds to several hundreds of runs for both the SOCP and the active set algorithms since we compute a 5-fold cross-validation for each regularization parameter of the (approximate) regularization path. on 5 random settings \{\mathbf{J},\mbox{{\mathbf{w}}},X,\varepsilon\}.

We assume that we have a rough idea of the level of sparsity of the true vector and we set the stopping criterion $s=4|\mathbf{J}|$ (see Algorithm 3), which is a rather conservative choice. We show on Figure 10 that we considerably lower the computational cost for the same level of performanceWe have not displayed this second figure that is just the superposition of the error curves for the SOCP and the active set algorithms.. As predicted by the complexity analysis of the active set algorithm (see the end of Section 4), considering the set of rectangular groups with or without the $\pm\pi/4$-groups results in the same complexity (up to constant terms). We empirically obtain an average complexity of $\,\approx O(p^{2.13})$ for the SOCP solver and of $\,\approx O(p^{0.45})$ for the active set algorithm.

Not surprisingly, for small values of $p$, the SOCP solver is faster than the active set algorithm, since the latter has to check its optimality by computing necessary and sufficient conditions (see Algorithm 3 and the discussion in the algorithmic complexity paragraph of Section 4).

Conclusion

A natural extension to this work is to consider bootstrapping since this may improve theoretical guarantees and result in better variable selection (Bach, 2008a; Meinshausen and Bühlmann, 2008). In order to deal with broader families of (non)zero patterns, it would be interesting to combine our approach with recent work on structured sparsity: for instance, Baraniuk et al. (2008); Jacob et al. (2009) consider union-closed collections of nonzero patterns, He and Carin (2009) exploit structure through a Bayesian prior while Huang et al. (2009) handle non-convex penalties based on information-theoretic criteria.

More generally, our regularization scheme could also be used for various learning tasks, as soon as prior knowledge on the structure of the sparse representation is available, e.g., for multiple kernel learning (Micchelli and Pontil, 2006), multi-task learning (Argyriou et al., 2008; Obozinski et al., 2009; Kim and Xing, 2009) and sparse matrix factorization problems (Mairal et al., 2010; Jenatton et al., 2010).

Finally, although we have mostly explored in this paper the algorithmic and theoretical issues related to these norms, this type of prior knowledge is of clear interest for the spatially and temporally structured data typical in bioinformatics, computer vision and neuroscience applications (see, e.g., Jenatton et al., 2010).

Acknowledgments

We would like to thank the anonymous reviewers for their constructive comments that improve the clarity and the overall quality of the manuscript. We also thank Julien Mairal and Guillaume Obozinski for insightful discussions. This work was supported in part by a grant from the Agence Nationale de la Recherche (MGA Project) and a grant from the European Research Council (SIERRA Project).

A Proof of Proposition 1

Still by using that a sum of convex functions is constant on a segment if and only if the functions are linear on this segment, the proof can be extended in order to replace the alternative assumption “$\{1,\dots,p\}$ belongs to $\mathcal{G}$” by the weaker but more involved assumption: for any $(j,k)\in\{1,\dots,p\}^{2}$, there exists a group $G\in\mathcal{G}$ which contains both $j$ and $k$.

B Proof of Theorem 2

the Hessian of $L_{J}$, which is positive definite (still from the same argument as in the proof of Theorem 1),

with $\psi=(\psi_{j})_{j\in J}$ a continuously differentiable function satisfying the matricial relation

C Proof of the minimality of the Backward procedure (see Algorithm 1)

There are essentially two points to show:

The first point can be shown by a proof by recurrence on the depth of the DAG. At step $t$, the base $\mathcal{G}^{(t)}$ verifies $\{\bigcup_{G\in\mathcal{G}^{\prime}}G,\ \forall\mathcal{G}^{\prime}\subseteq\mathcal{G}^{(t)}\}=\{G\in\mathcal{Z},|G|\leq t\}$ because an element $G\in\mathcal{Z}$ is either the union of itself or the union of elements strictly smaller. The initialization $t=\min_{G\in\mathcal{Z}}|G|$ is easily verified, the leafs of the DAG being necessarily in $\mathcal{G}$.

As for the second point, we proceed by contradiction. If there exists another base $\mathcal{G}^{*}$ that spans $\mathcal{Z}$ such that $\mathcal{G}^{*}\subset\mathcal{G}$, then

By definition of the set $\mathcal{Z}$, there exists in turn $\mathcal{G}^{\prime}\subseteq\mathcal{G}^{*},\ \mathcal{G}^{\prime}\neq\{e\}$ (otherwise, $e$ would belong to $\mathcal{G}^{*}$), verifying $e=\bigcup_{G\in\mathcal{G}^{\prime}}G$, which is impossible by construction of $\mathcal{G}$ whose members cannot be the union of elements of $\mathcal{Z}$.

D Proof of Proposition 3

The proposition comes from a classic result of Fenchel Duality (Borwein and Lewis, 2006, Theorem 3.3.5 and Exercise 3.3.9) when we consider the convex function

whose Fenchel conjugate $h^{*}_{J}$ is given by $\kappa_{J}\mapsto\frac{1}{2\lambda}\left[\Omega_{J}^{*}(\kappa_{J})\right]^{2}$ (Boyd and Vandenberghe, 2004, example 3.27). Since the set

is not empty, we get the first part of the proposition. Moreover, the primal-dual variables $\{w_{J},\kappa_{J}\}$ is optimal if and only if

where $\partial\Omega_{J}(w_{J})$ denotes the subdifferential of $\Omega_{J}$ at $w_{J}$. The differentiability of $L_{J}$ at $w_{J}$ then gives $\partial L_{J}(w_{J})=\{\nabla L_{J}(w_{J})\}$. It now remains to show that

As a starting point, the Fenchel-Young inequality (Borwein and Lewis, 2006, Proposition 3.3.4) gives the equivalence between Eq. (8) and

Thus, if $\kappa_{J}\in\lambda\Omega_{J}(w_{J})\partial\Omega_{J}(w_{J})$ then $w_{J}^{\top}\kappa_{J}=\lambda\left[\Omega_{J}(w_{J})\right]^{2}.$ Combined with Eq. (10), we obtain $w_{J}^{\top}\kappa_{J}=\frac{1}{\lambda}\left[\Omega_{J}^{*}(\kappa_{J})\right]^{2}.$

Reciprocally, starting from Eq. (9), we notably have

In light of Eq. (11), it suffices to check that $\Omega_{J}^{*}(\kappa_{J})\leq\lambda\Omega_{J}(w_{J})$ in order to have Eq. (8). Combining Eq. (9) with the definition of the dual norm, it comes

which concludes the proof of the equivalence between Eq. (8) and Eq. (9).

E Proofs of Propositions 4 and 5

In order to check that the reduced solution $w_{J}$ is optimal for the full problem in Eq. (5), we complete with zeros on ${J^{c}}$ to define $w$, compute $\kappa=-\nabla L(w)$, which is such that $\kappa_{J}=-\nabla L_{J}(w_{J})$, and get a duality gap for the full problem equal to

By designing upper and lower bounds for $\Omega^{*}(\kappa)$, we get sufficient and necessary conditions.

By combining Lemma 13 and the fact that $\mathcal{G}_{K\backslash J}\cap(\mathcal{G}_{J})^{c}=\mathcal{G}_{K}\backslash\mathcal{G}_{J}$, we have for all $G\in\mathcal{G}_{K}\backslash\mathcal{G}_{J}$, $K\backslash J\subseteq G$ and therefore $u_{G\cap K\backslash J}=u_{K\backslash J}$. Since we cannot compute the dual norm of $u_{K\backslash J}\mapsto\|d^{\scriptscriptstyle G}_{K\backslash J}\circ u_{K\backslash J}\|_{2}$ in closed-form, we instead use the following upperbound

Finally, Proposition 3 gives \lambda\Omega(w^{*})=\big{\{}\!\!-\lambda{w^{*}}^{\top}\nabla L(w^{*})\big{\}}^{\frac{1}{2}}, which leads to the desired result.

E.2 Proof of Proposition 5

The goal of the proof is to upper bound the dual norm $\Omega^{*}(\kappa)$ by taking advantage of the structure of $\mathcal{G}$; we first show how we can upper bound $\Omega^{*}(\kappa)$ by $(\Omega_{J}^{c})^{\ast}[\kappa_{J^{c}}]$. We indeed have:

where in the last line, we use Lemma 15. Thus the duality gap is less than

and a sufficient condition for the duality gap to be smaller than $\varepsilon$ is

Using Proposition 3, we have $-\lambda w^{\top}\nabla L(w)=\left[\Omega_{J}^{*}(\kappa_{J})\right]^{2}$ and we get the right-hand side of Proposition 5. It now remains to upper bound $(\Omega_{J}^{c})^{\ast}[\kappa_{J^{c}}]$. To this end, we call upon Lemma 11 to obtain:

Among all groups $G\in(\mathcal{G}_{J})^{c}$, the ones with the maximum values are the largest ones, i.e., those in the fringe groups $\mathcal{F}_{J}=\{G\in(\mathcal{G}_{J})^{c}\ ;\ \nexists G^{\prime}\in(\mathcal{G}_{J})^{c},G\subseteq G^{\prime}\}$. This argument leads to the result of Proposition 5.

F Proof of Theorem 6

Since $\mu\rightarrow 0$, we also have by taking the directional derivative of $\Omega$ at ${\mathbf{w}}$ in the direction of $\Delta$

By assumption, with probability tending to one, we have $\mathbf{J}=\{j\in\{1,\dots,p\},\hat{w}_{j}\neq 0\}$, hence for any $j\in\mathbf{J}^{c}$ \mu\hat{\Delta}_{j}=(\mbox{\hat{w}}-\mbox{{\mathbf{w}}})_{j}=0. This implies that the nonrandom vector $\Delta^{\ast}$ verifies $\Delta_{\mathbf{J}^{c}}^{\ast}=0$.

Sufficient condition: We turn to the sufficient condition. We first consider the problem reduced to the hull $\mathbf{J}$,

It remains to show that $\hat{w}$ is indeed optimal for the full problem (that admits a unique solution due to the positiveness of $Q$). By construction, the optimality condition (see Lemma 14) relative to the active variables $\mathbf{J}$ is already verified. More precisely, we have

Since we assume $(\Omega_{\mathbf{J}}^{c})^{\ast}[\mathbf{Q}_{\mathbf{J}^{c}\mathbf{J}}\mathbf{Q}_{\mathbf{J}\mathbf{J}}^{-1}\mathbf{r}_{\mathbf{J}}]<1$, we obtain

which proves the optimality condition of Lemma 14 relative to the inactive variables: $\hat{w}$ is therefore optimal for the full problem.

G Proof of Theorem 7

Since our analysis takes place in a finite-dimensional space, all the norms defined on this space are equivalent. Therefore, we introduce the equivalence parameters $a(\mathbf{J}),A(\mathbf{J})>0$ such that

For any matrix $\Gamma$, we also introduce the operator norm $\|\Gamma\|_{m,s}$ defined as

Following Bach (2008b) and Nardi and Rinaldo (2008), we consider the reduced problem on $\mathbf{J}$,

with solution $\hat{w}_{\mathbf{J}}$, which can be extended to $\mathbf{J}^{c}$ with zeros. From optimality conditions (see Lemma 14), we know that

hence from (12) and the definition of $A(\mathbf{J})$,

Thus, if we assume $\mu\leq\frac{\kappa\nu}{3|\mathbf{J}|^{1/2}A(\mathbf{J})}$ and

so that for all $G\in\mathcal{G}_{\mathbf{J}}$, \left\|\mbox{\hat{w}}_{G}\right\|_{\infty}\geq\frac{\nu}{3}, hence the hull is indeed selected.

This also ensures that $\hat{w}_{\mathbf{J}}$ satisfies the equation (see Lemma 14)

We now prove that the $\hat{w}$ padded with zeros on $\mathbf{J}^{c}$ is indeed optimal for the full problem with high probability. According to Lemma 14, since we have already proved (16), it suffices to show that

Defining $q_{\mathbf{J}^{c}|\mathbf{J}}=q_{\mathbf{J}^{c}}-Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1}q_{\mathbf{J}}$, we can write the gradient of $L$ on $\mathbf{J}^{c}$ as

which leads us to control the difference $\hat{r}_{\mathbf{J}}-\mathbf{r}_{\mathbf{J}}$. Using Lemma 12, we get

where w=t_{0}\mbox{\hat{w}}+(1-t_{0})\mbox{{\mathbf{w}}}\, for some $\,t_{0}\in(0,1)$.

Let {\mathbf{\overline{J}}}=\{k\in\mathbf{J}:\mbox{{\mathbf{w}}}_{k}\neq 0\} and let $\varphi$ be defined as

The term $\varphi$ basically measures how close $\mathbf{J}$ and ${\mathbf{\overline{J}}}$ are, i.e., how relevant the prior encoded by $\mathcal{G}$ about the hull $\mathbf{J}$ is. By using (15), we have

Introducing \alpha=\frac{18\varphi^{3/2}\|\mbox{{\mathbf{w}}}\|_{\infty}}{\nu^{2}}\sum_{G\in\mathcal{G}_{\mathbf{J}}}\left\|d^{\scriptscriptstyle G}_{\mathbf{J}}\right\|_{2}, we thus have proved

By writing the Schur complement of $Q$ on the block matrices $Q_{\mathbf{J}^{c}\mathbf{J}^{c}}$ and $Q_{\mathbf{J}\mathbf{J}}$, the positiveness of $Q$ implies that the diagonal terms ${\rm diag}(Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1}Q_{\mathbf{J}\mathbf{J}^{c}})$ are less than one, which results in $\|Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1/2}\|_{\infty,2}\leq 1$. We then have

where the last line comes from Eq. (13) and (17). We get

Thus, if the following inequalities are verified

Combined with earlier constraints, this leads to the first part of the desired proposition.

We now need to make sure that the conditions (14), (23) and (24) hold with high probability. To this end, we upperbound, using Gaussian concentration inequalities, two tail-probabilities. First, $q_{\mathbf{J}^{c}|\mathbf{J}}$ is a centered Gaussian random vector with covariance matrix

where $Q_{\mathbf{J}^{c}\mathbf{J}^{c}|\mathbf{J}}=Q_{\mathbf{J}^{c}\mathbf{J}^{c}}-Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1}Q_{\mathbf{J}\mathbf{J}^{c}}$. In particular, $(\Omega_{\mathbf{J}}^{c})^{\ast}[q_{\mathbf{J}^{c}|\mathbf{J}}]$ has the same distribution as $\psi(W)$, with $\psi:u\mapsto(\Omega_{\mathbf{J}}^{c})^{\ast}(\sigma n^{-1/2}Q_{\mathbf{J}^{c}\mathbf{J}^{c}|\mathbf{J}}^{1/2}u)$ and $W$ a centered Gaussian random variable with unit covariance matrix.

Since for any $u$ we have $u^{\top}Q_{\mathbf{J}^{c}\mathbf{J}^{c}|\mathbf{J}}u\leq u^{\top}Q_{\mathbf{J}^{c}\mathbf{J}^{c}}u\leq\left\|Q^{1/2}\right\|_{2}^{2}\left\|u\right\|_{2}^{2}$, by using Sudakov-Fernique inequality (Adler, 1990, Theorem 2.9), we get:

On the other hand, since $Q$ has unit diagonal and $Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1}Q_{\mathbf{J}\mathbf{J}^{c}}$ has diagonal terms less than one, $Q_{\mathbf{J}^{c}\mathbf{J}^{c}|\mathbf{J}}$ also has diagonal terms less than one, which implies that $\|Q_{\mathbf{J}^{c}\mathbf{J}^{c}|\mathbf{J}}^{1/2}\|_{\infty,2}\leq 1$. Hence $\psi$ is a Lipschitz function with Lipschitz constant upper bounded by $\sigma n^{-1/2}a(\mathbf{J}^{c})^{-1}$. Thus by concentration of Lipschitz functions of multivariate standard random variables (Massart, 2003, Theorem 3.4), we have for $t>0$:

We can apply classical inequalities for standard random variables (Massart, 2003, Theorem 3.4) that directly lead to

where we recall the definitions: $W$ a centered Gaussian random variable with unit covariance matrix, {\mathbf{\overline{J}}}=\{j\in\mathbf{J}:\mbox{{\mathbf{w}}}_{j}\neq 0\}, \nu=\min\{|\mbox{{\mathbf{w}}}_{j}|;\ j\in{\mathbf{\overline{J}}}\},

$\kappa=\lambda_{\min}(Q_{\mathbf{J}\mathbf{J}})>0$ and $\tau>0$ such that $(\Omega_{\mathbf{J}}^{c})^{\ast}[Q_{\mathbf{J}^{c}\mathbf{J}}Q_{\mathbf{J}\mathbf{J}}^{-1}\mathbf{r}]<1-\tau$.

H A first order approach to solve Eq. (2) and Eq. (5)

whose minimum is uniquely attained for $z_{j}=|x_{j}|/\left\|x\right\|_{1}$. Similarly, we have

whose minimun is uniquely obtained for $z_{j}=|x_{j}|$. Thus, we can equivalently rewrite Eq. (2) as

with $\zeta_{j}=(\sum_{G\ni j}(d^{\scriptscriptstyle G}_{j})^{2}(\eta^{\scriptscriptstyle G})^{-1})^{-1}$. In the same vein, Eq. (5) is equivalent to

where $\zeta_{j}$ is defined as above. The reformulations Eq. (28) and Eq. (29) lend themselves well to a simple alternating optimization scheme between $w$ (for instance, $w$ can be computed in closed-form when the square loss is used) and $(\eta^{\scriptscriptstyle G})_{G\in\mathcal{G}}$ (whose optimal value is always a closed-form solution).

This first order approach is computationally appealing since it allows warm-restart, which can dramatically speed up the computation over regularization paths.

I Technical lemmas

In this last section of the appendix, we give several technical lemmas. We consider $I\subseteq\{1,\dots,p\}$ and $\mathcal{G}_{I}=\{G\in\mathcal{G};\ G\cap I\neq\varnothing\}\subseteq\mathcal{G}$, i.e., the set of active groups when the variables $I$ are selected.

We begin with a dual formulation of $\Omega^{*}$ obtained by conic duality (Boyd and Vandenberghe, 2004):

Proof By definiton of $(\Omega_{I})^{*}[u_{I}]$, we have

which is a second-order cone program (SOCP) with $|\mathcal{G}_{I}|$ second-order cone constraints. This primal problem is convex and satisfies Slater’s conditions for generalized conic inequalities, which implies that strong duality holds (Boyd and Vandenberghe, 2004). We now consider the Lagrangian $\mathcal{L}$ defined as

After simplifying the Lagrangian, the dual problem then reduces to

which is equivalent to the displayed result. Since we cannot compute in closed-form the solution of the previous optimization problem, we focus on a different but closely related problem, i.e., when we replace the objective $\max_{G\in\mathcal{G}_{I}}\left\|\xi_{I}^{\scriptscriptstyle G}\right\|_{2}$ by $\max_{G\in\mathcal{G}_{I}}\left\|\xi_{I}^{\scriptscriptstyle G}\right\|_{\infty}$, to obtain a meaningful feasible point:

is minimized for $(\xi_{j}^{\scriptscriptstyle G})^{*}=-\dfrac{u_{j}}{\sum_{H\in j,H\in\mathcal{G}_{I}}\!\!d^{\scriptscriptstyle H}_{j}}.$

Proof We proceed by contradiction. Let us assume there exists $(\xi_{I}^{\scriptscriptstyle G})_{G\in\mathcal{G}_{I}}$ such that

where we denote by $j_{0}$ an argmax of the latter maximization. We notably have for all $G\ni j_{0}$:

By multiplying both sides by $d^{\scriptscriptstyle G}_{j_{0}}$ and by summing over $G\ni j_{0}$, we get

We now give an upperbound on $\Omega^{*}$ based on Lemma 9 and Lemma 10:

Proof We simply plug the minimizer obtained in Lemma 10 into the problem of Lemma 9.

We now derive a lemma to control the difference of the gradient of $\Omega_{J}$ evaluated in two points:

Given an active set $J\subseteq\{1,\dots,p\}$ and a direct parent $K\in\Pi_{\mathcal{P}}(J)$ of $J$ in the DAG of nonzero patterns, we have the following result:

For all $G\in\mathcal{G}_{K}\backslash\mathcal{G}_{J}$, we have

Proof We proceed by contradiction. We assume there exists $G_{0}\in\mathcal{G}_{K}\backslash\mathcal{G}_{J}$ such that $K\backslash J\nsubseteq G_{0}$. Given that $K\in\mathcal{P}$, there exists $\mathcal{G}^{\prime}\subseteq\mathcal{G}$ verifying $K=\bigcap_{G\in\mathcal{G}^{\prime}}G^{c}$. Note that $G_{0}\notin\mathcal{G}^{\prime}$ since by definition $G_{0}\cap K\neq\varnothing$.

which is impossible by definition of $K$.

We give below an important Lemma to characterize the solutions of (2).

In addition, the solution $\hat{w}$ satisfies

Some algebra leads to the following equivalent formulation

The first part of the lemma then comes from the projections on $\hat{J}$ and $\hat{J}^{c}$.

We end up with a lemma regarding the dual norm of the sum of two disjoint norms (see Rockafellar, 1970):

References