Convex and Network Flow Optimization for Structured Sparsity

Introduction

We present proximal splitting methods for solving structured sparse regularized problems.

We demonstrate that our methods are relevant for various applications whose practical success is made possible by our algorithmic tools and efficient implementations. First, we introduce a new CUR matrix factorization technique exploiting structured sparse regularization, built upon the links drawn by Bien et al. 2010 between CUR decomposition (Mahoney and Drineas 2009) and sparse regularization. Then, we illustrate our algorithms with different tasks: video background subtraction, estimation of hierarchical structures for dictionary learning of natural image patches (Jenatton et al. 2010a; Jenatton et al. 2011), wavelet image denoising with a structured sparse prior, and topographic dictionary learning of natural image patches (Hyvärinen et al. 2001; Kavukcuoglu et al. 2009; Garrigues and Olshausen 2010).

Note that this paper extends a shorter version published in Advances in Neural Information Processing Systems (Mairal et al. 2010b), by adding new experiments (CUR matrix factorization, wavelet image denoising and topographic dictionary learning), presenting the proximal splitting methods, providing the full proofs of the optimization results, and adding numerous discussions.

The rest of this paper is organized as follows: Section 2 presents structured sparse models and related work. Section 3 is devoted to proximal gradient algorithms, and Section 4 to proximal splitting methods. Section 5 presents several experiments and applications demonstrating the effectiveness of our approach and Section 6 concludes the paper.

Structured Sparse Models

We are interested in machine learning problems where the solution is not only known beforehand to be sparse—that is, the solution has only a few non-zero coefficients, but also to form non-zero patterns with a specific structure. It is indeed possible to encode additional knowledge in the regularization other than just sparsity. For instance, one may want the non-zero patterns to be structured in the form of non-overlapping groups (Turlach et al. 2005; Yuan and Lin 2006; Stojnic et al. 2009; Obozinski et al. 2010), in a tree (Zhao et al. 2009; Bach 2009; Jenatton et al. 2010a; Jenatton et al. 2011), or in overlapping groups (Jenatton et al. 2009; Jacob et al. 2009; Huang et al. 2009; Baraniuk et al. 2010; Cehver et al. 2008; He and Carin 2009), which is the setting we are interested in here.

As for classical non-structured sparse models, 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.

A first approach introduced by Baraniuk et al. 2010 consists in imposing that the sparsity pattern of a solution (i.e., its set of non-zero coefficients) is in a predefined subset of groups of variables ${\mathcal{G}}\subseteq 2^{\{1,\ldots,p\}}$. Given this a priori knowledge, a greedy algorithm (Needell and Tropp 2009) is used to address the following nonconvex structured sparse decomposition problem

Intuitively, this formulation encourages solutions ${\mathbf{w}}$ whose sparsity patterns have a small coding length, meaning in practice that they can be represented by a union of a small number of groups. Even though they are related, this model is different from the one of Baraniuk et al. 2010.

These two approaches are encoding a priori knowledge on the shape of non-zero patterns that the solution of a regularized problem should have. A different point of view consists of modelling the zero patterns of the solution—that is, define groups of variables that should be encouraged to be set to zero together. After defining a set ${\mathcal{G}}\subseteq 2^{\{1,\ldots,p\}}$ of such groups of variables, the following penalty can naturally be used as a regularization to induce the desired property

2 Convex Approaches with Sparsity-Inducing Norms

In this paper, we are interested in convex regularizations which induce structured sparsity. Generally, we consider the following optimization problem

If ${\mathcal{G}}$ is a partition of $\{1,\ldots,p\}$, i.e. the groups do not overlap, variables are selected in groups rather than individually. When the coefficients of the solution are known to be organized in such a way, explicitly encoding the a priori group structure in the regularization can improve the prediction performance and/or interpretability of the learned models (Turlach et al. 2005; Yuan and Lin 2006; Roth and Fischer 2008; Stojnic et al. 2009; Huang and Zhang 2010; Obozinski et al. 2010). Such a penalty is commonly called group-Lasso penalty.

When the groups overlap, $\Omega$ is still a norm and sets groups of variables to zero together Jenatton et al. 2009. The latter setting has first been considered for hierarchies Zhao et al. 2009; Kim and Xing 2010; Bach 2009; Jenatton et al. 2010a; Jenatton et al. 2011, and then extended to general group structures Jenatton et al. 2009. Solving Eq. (2) in this context is a challenging problem which is the topic of this paper.

Note that other types of structured-sparsity inducing norms have also been introduced, notably the approach of Jacob et al. 2009, which penalizes the following quantity

This penalty, which is also a norm, can be seen as a convex relaxation of the regularization introduced by Huang et al. 2009, and encourages the sparsity pattern of the solution to be a union of a small number of groups. Even though both $\Omega$ and $\Omega^{\prime}$ appear under the terminology of “structured sparsity with overlapping groups”, they have in fact significantly different purposes and algorithmic treatments. For example, Jacob et al. 2009 consider the problem of selecting genes in a gene network which can be represented as the union of a few predefined pathways in the graph (groups of genes), which overlap. In this case, it is natural to use the norm $\Omega^{\prime}$ instead of $\Omega$. On the other hand, we present a matrix factorization task in Section 5.3, where the set of zero-patterns should be a union of groups, naturally leading to the use of $\Omega$. Dealing with $\Omega^{\prime}$ is therefore relevant, but out of the scope of this paper.

3 Convex Optimization Methods Proposed in the Literature

Generic approaches to solve Eq. (2) mostly rely on subgradient descent schemes (Bertsekas 1999, see), and interior-point methods Boyd and Vandenberghe 2004. These generic tools do not scale well to large problems and/or do not naturally handle sparsity (the solutions they return may have small values but no “true” zeros). These two points prompt the need for dedicated methods.

Problem (2) has also been addressed with working-set algorithms (Bach 2009; Jenatton et al. 2009; Schmidt and Murphy 2010). The main idea of these methods is to solve a sequence of increasingly larger subproblems of (2). Each subproblem consists of an instance of Eq. (2) reduced to a specific subset of variables known as the working set. As long as some predefined optimality conditions are not satisfied, the working set is augmented with selected inactive variables (Bach et al. 2011, for more details, see).

Optimization with Proximal Gradient Methods

We address in this section the problem of solving Eq. (2) under the following assumptions:

$f$ is differentiable with Lipschitz-continuous gradient. For machine learning problems, this hypothesis holds when $f$ is for example the square, logistic or multi-class logistic loss (Shawe-Taylor and Cristianini 2004, see).

To the best of our knowledge, no dedicated optimization method has been developed for this setting. Following Jenatton et al. 2010a; Jenatton et al. 2011 who tackled the particular case of hierarchical norms, we propose to use proximal gradient methods, which we now introduce.

Proximal methods have drawn increasing attention in the signal processing (e.g., Wright et al. 2009b; 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 (Nesterov 2007; Beck and Teboulle 2009, e.g.,).

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 a upper bound on the Lipschitz constant of $\nabla f$. This problem can be equivalently rewritten as:

Solving efficiently and exactly this problem allows to attain the fast convergence rates of proximal methods, i.e., reaching a precision of $O(\frac{L}{k^{2}})$ in $k$ iterations. Note, however, that fast convergence rates can also be achieved while solving approximately the proximal problem (see Schmidt et al. 2011, for more details). 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 to solve sparse decomposition problems is that this operator can often be computed in closed form. For instance,

2 Dual of the Proximal Operator

We now show that, for a set $\mathcal{G}$ of general overlapping groups, a convex dual of the proximal problem (4) can be reformulated as a quadratic min-cost flow problem. We then propose an efficient algorithm to solve it exactly, as well as a related algorithm to compute the dual norm of $\Omega$. We start by considering the dual formulation to problem (4) introduced by Jenatton et al. 2010a; Jenatton et al. 2011:

Without loss of generality, Let ${\boldsymbol{\xi}}^{\star}$ denote a solution of Eq. (6). Optimality conditions of Eq. (6) derived in Jenatton et al. 2010a; Jenatton et al. 2011 show that for all $j$ in $\{1,\ldots,p\}$, the signs of the non-zero coefficients ${\boldsymbol{\xi}}_{j}^{\star g}$ for $g$ in ${\mathcal{G}}$ are the same as the signs of the entries ${\mathbf{u}}_{j}$. To solve Eq. (6), one can therefore flip the signs of the negative variables ${\mathbf{u}}_{j}$, then solve the modified dual formulation (with non-negative variables), which gives the magnitude of the entries ${\boldsymbol{\xi}}_{j}^{\star g}$ (the signs of these being known). we assume from now on that the scalars ${\mathbf{u}}_{j}$ are all non-negative, and we constrain the entries of ${\boldsymbol{\xi}}$ to be so. Such a formulation introduces $p|\mathcal{G}|$ dual variables which can be much greater than $p$, the number of primal variables, but it removes the issue of overlapping regularization. We now associate a graph with problem (6), on which the variables ${\boldsymbol{\xi}}_{j}^{g}$, for $g$ in ${\mathcal{G}}$ and $j$ in $g$, can be interpreted as measuring the components of a flow.

3 Graph Model

Let $G$ be a directed graph $G=(V,E,s,t)$, where $V$ is a set of vertices, $E\subseteq V\times V$ a set of arcs, $s$ a source, and $t$ a sink. For all arcs in $E$, we define a non-negative capacity constant, and as done classically in the network flow literature (Ahuja et al. 1993; Bertsekas 1998), we define a flow as a non-negative function on arcs that satisfies capacity constraints on all arcs (the value of the flow on an arc is less than or equal to the arc capacity) and conservation constraints on all vertices (the sum of incoming flows at a vertex is equal to the sum of outgoing flows) except for the source and the sink. For every arc $e$ in $E$, we also define a real-valued cost function, which depends on the value of the flow on $e$. We now introduce the canonical graph $G$ associated with our optimization problem:

Let ${\mathcal{G}}\subseteq\{1,\ldots,p\}$ be a set of groups, and $(\eta_{g})_{g\in{\mathcal{G}}}$ be positive weights. The canonical graph $G=(V,E,s,t)$ is the unique graph defined as follows:

$V=V_{u}\cup V_{gr}$, where $V_{u}$ is a vertex set of size $p$, one vertex being associated to each index $j$ in $\{1,\ldots,p\}$, and $V_{gr}$ is a vertex set of size $|{\mathcal{G}}|$, one vertex per group $g$ in ${\mathcal{G}}$. We thus have $|V|=|{\mathcal{G}}|+p$. For simplicity, we identify groups $g$ in ${\mathcal{G}}$ and indices $j$ in $\{1,\ldots,p\}$ with vertices of the graph, such that one can from now on refer to “vertex $j$” or “vertex $g$”.

For every group $g$ in ${\mathcal{G}}$, $E$ contains an arc $(s,g)$. These arcs have capacity $\lambda\eta_{g}$ and zero cost.

For every group $g$ in ${\mathcal{G}}$, and every index $j$ in $g$, $E$ contains an arc $(g,j)$ with zero cost and infinite capacity. We denote by ${\boldsymbol{\xi}}_{j}^{g}$ the flow on this arc.

For every index $j$ in $\{1,\ldots,p\}$, $E$ contains an arc $(j,t)$ with infinite capacity and a cost $\frac{1}{2}({\mathbf{u}}_{j}-{\boldsymbol{\overline{\xi}}}_{j})^{2}$, where ${\boldsymbol{\overline{\xi}}}_{j}$ is the flow on $(j,t)$.

Examples of canonical graphs are given in Figures 1a- for three simple group structures. The flows ${\boldsymbol{\xi}}_{j}^{g}$ associated with $G$ can now be identified with the variables of problem (6). Since we have assumed the entries of ${\mathbf{u}}$ to be non-negative, we can now reformulate Eq. (6) as

the only arcs with a cost are those leading to the sink, which have the form $(j,t)$, where $j$ is the index of a variable in $\{1,\ldots,p\}$. The sum of these costs is $\sum_{j=1}^{p}\frac{1}{2}({\mathbf{u}}_{j}-{\boldsymbol{\overline{\xi}}}_{j})^{2}$, which is the objective function minimized in Eq. (7);

by flow conservation, we necessarily have ${\boldsymbol{\overline{\xi}}}_{j}=\sum_{g\in{\mathcal{G}}}{{\boldsymbol{\xi}}_{j}^{g}}$ in the canonical graph;

the only arcs with a capacity constraints are those coming out of the source, which have the form $(s,g)$, where $g$ is a group in ${\mathcal{G}}$. By flow conservation, the flow on an arc $(s,g)$ is $\sum_{j\in g}{\boldsymbol{\xi}}_{j}^{g}$ which should be less than $\lambda\eta_{g}$ by capacity constraints;

all other arcs have the form $(g,j)$, where $g$ is in ${\mathcal{G}}$ and $j$ is in $g$. Thus, $\operatorname{Supp}({\boldsymbol{\xi}}^{g})\subseteq g$.

Therefore we have shown that finding a flow minimizing the sum of the costs on such a graph is equivalent to solving problem (6). When some groups are included in others, the canonical graph can be simplified to yield a graph with a smaller number of edges. Specifically, if $h$ and $g$ are groups with $h\subset g$, the edges $(g,j)$ for $j\in h$ carrying a flow ${\boldsymbol{\xi}}^{g}_{j}$ can be removed and replaced by a single edge $(g,h)$ of infinite capacity and zero cost, carrying the flow $\sum_{j\in h}{\boldsymbol{\xi}}^{g}_{j}$. This simplification is illustrated in Figure 1d, with a graph equivalent to the one of Figure 1c. This does not change the optimal value of ${\boldsymbol{\overline{\xi}}}^{\star}$, which is the quantity of interest for computing the optimal primal variable ${\mathbf{w}}^{\star}$. We present in Appendix A a formal definition of equivalent graphs. These simplifications are useful in practice, since they reduce the number of edges in the graph and improve the speed of our algorithms.

4 Computation of the Proximal Operator

The general case of overlapping groups is more difficult. Hochbaum and Hong 1995 have shown that quadratic min-cost flow problems can be reduced to a specific parametric max-flow problem, for which an efficient algorithm exists (Gallo et al. 1989). By definition, a parametric max-flow problem consists in solving, for every value of a parameter, a max-flow problem on a graph whose arc capacities depend on this parameter. While this generic approach could be used to solve Eq. (6), we propose to use Algorithm 1 that also exploits the fact that our graphs have non-zero costs only on edges leading to the sink. As shown in Appendix D, it it has a significantly better performance in practice. This algorithm clearly shares some similarities with existing approaches in network flow optimization such as the simplified version of Gallo et al. 1989 presented by Babenko and Goldberg 2006 that uses a divide and conquer strategy. Moreover, an equivalent algorithm exists for minimizing convex functions over polymatroid sets (Groenevelt 1991). This equivalence, a priori non trivial, is uncovered through a representation of structured sparsity-inducing norms via submodular functions, which was recently proposed by Bach 2010.

The intuition behind our algorithm, computeFlow (see Algorithm 1), is the following: since ${\boldsymbol{\overline{\xi}}}=\sum_{g\in\mathcal{G}}{\boldsymbol{\xi}}^{g}$ is the only value of interest to compute the solution of the proximal operator ${\mathbf{w}}={\mathbf{u}}-{\boldsymbol{\overline{\xi}}}$, the first step looks for a candidate value ${\boldsymbol{\gamma}}$ for ${\boldsymbol{\overline{\xi}}}$ by solving the following relaxed version of problem (7):

The cost function here is the same as in problem (7), but the constraints are weaker: Any feasible point of problem (7) is also feasible for problem (8). This problem can be solved in linear time Brucker 1984. Its solution, which we denote ${\boldsymbol{\gamma}}$ for simplicity, provides the lower bound $\|{\mathbf{u}}-{\boldsymbol{\gamma}}\|_{2}^{2}/2$ for the optimal cost of problem (7).

The second step tries to construct a feasible flow $({\boldsymbol{\xi}},{\boldsymbol{\overline{\xi}}})$, satisfying additional capacity constraints equal to $\gamma_{j}$ on arc $(j,t)$, and whose cost matches this lower bound; this latter problem can be cast as a max-flow problem (Goldberg and Tarjan 1986). If such a flow exists, the algorithm returns ${\boldsymbol{\overline{\xi}}}={\boldsymbol{\gamma}}$, the cost of the flow reaches the lower bound, and is therefore optimal. If such a flow does not exist, we have ${\boldsymbol{\overline{\xi}}}\neq{\boldsymbol{\gamma}}$, the lower bound is not achievable, and we build a minimum $(s,t)$-cut of the graph Ford and Fulkerson 1956 defining two disjoints sets of nodes $V^{+}$ and $V^{-}$; $V^{+}$ is the part of the graph which is reachable from the source (for every node $j$ in ${\mathbf{V}}^{+}$, there exists a non-saturated path from $s$ to $j$), whereas all paths going from $s$ to nodes in $V^{-}$ are saturated. More details about these properties can be found at the beginning of Appendix B. At this point, it is possible to show that the value of the optimal min-cost flow on all arcs between $V^{+}$ and $V^{-}$ is necessary zero. Thus, removing them yields an equivalent optimization problem, which can be decomposed into two independent problems of smaller sizes and solved recursively by the calls to computeFlow$(V^{+},E^{+})$ and computeFlow$(V^{-},E^{-})$. A formal proof of correctness of Algorithm 1 and further details are relegated to Appendix B.

The approach of Hochbaum and Hong 1995; Gallo et al. 1989 which recasts the quadratic min-cost flow problem as a parametric max-flow is guaranteed to have the same worst-case complexity as a single max-flow algorithm. However, we have experimentally observed a significant discrepancy between the worst case and empirical complexities for these flow problems, essentially because the empirical cost of each max-flow is significantly smaller than its theoretical cost. Despite the fact that the worst-case guarantees for our algorithm is weaker than theirs (up to a factor $|V|$), it is more adapted to the structure of our graphs and has proven to be much faster in our experiments (see Appendix D). The best theoretical worst-case complexity of a max-flow is achieved by Goldberg and Tarjan 1986 and is $O\big(|V||E|\log(|V|^{2}/|E|)\big)$. Our algorithm achieves the same worst-case complexity when the cuts are well balanced—that is $|V^{+}|\approx|V^{-}|\approx|V|/2$, but we lose a factor $|V|$ when it is not the case. The practical speed of such algorithms is however significantly different than their theoretical worst-case complexities (Boykov and Kolmogorov 2004, see). Some implementation details are also crucial to the efficiency of the algorithm:

Exploiting maximal connected components: When there exists no arc between two subsets of $V$, the solution can be obtained by solving two smaller optimization problems corresponding to the two disjoint subgraphs. It is indeed possible to process them independently to solve the global min-cost flow problem. To that effect, before calling the function computeFlow($V,E$), we look for maximal connected components $(V_{1},E_{1}),\ldots,(V_{N},E_{N})$ and call sequentially the procedure computeFlow($V_{i},E_{i}$) for $i$ in $\{1,\ldots,N\}$.

Efficient max-flow algorithm: We have implemented the “push-relabel” algorithm of Goldberg and Tarjan 1986 to solve our max-flow problems, using classical heuristics that significantly speed it up in practice; see Goldberg and Tarjan 1986 and Cherkassky and Goldberg 1997. We use the so-called “highest-active vertex selection rule, global and gap heuristics” (Goldberg and Tarjan 1986; Cherkassky and Goldberg 1997), which has a worst-case complexity of $O(|V|^{2}|E|^{1/2})$ for a graph $(V,E,s,t)$. This algorithm leverages the concept of pre-flow that relaxes the definition of flow and allows vertices to have a positive excess.

Using flow warm-restarts: The max-flow steps in our algorithm can be initialized with any valid pre-flow, enabling warm-restarts. This is also a key concept in the parametric max-flow algorithm of Gallo et al. 1989.

Improved projection step: The first line of the procedure computeFlow can be replaced by ${\boldsymbol{\gamma}}\leftarrow\operatornamewithlimits{arg\,min}_{\boldsymbol{\gamma}}\sum_{j\in V_{u}}\frac{1}{2}({\mathbf{u}}_{j}-{\boldsymbol{\gamma}}_{j})^{2}~~\text{s.t.}~~\sum_{j\in V_{u}}{\boldsymbol{\gamma}}_{j}\leq\lambda\sum_{g\in V_{gr}}\eta_{g}~\text{and}~|{\boldsymbol{\gamma}}_{j}|\leq\lambda\sum_{g\ni j}\eta_{g}.$ The idea is to build a relaxation of Eq. (7) which is closer to the original problem than the one of Eq. (8), but that still can be solved in linear time. The structure of the graph will indeed not allow ${\boldsymbol{\overline{\xi}}}_{j}$ to be greater than $\lambda\sum_{g\ni j}\eta_{g}$ after the max-flow step. This modified projection step can still be computed in linear time (Brucker 1984), and leads to better performance.

5 Computation of the Dual Norm

is a key quantity to study sparsity-inducing regularizations in many respects. For instance, dual norms are central in working-set algorithms (Jenatton et al. 2009; Bach et al. 2011), and arise as well when proving theoretical estimation or prediction guarantees (Negahban et al. 2009).

In our context, we use it to monitor the convergence of the proximal method through a duality gap, hence defining a proper optimality criterion for problem (2). 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,). A description of the algorithm we use in the experiments (Beck and Teboulle 2009) along with the integration of the computation of the duality gap is given in Appendix C.

We now denote by $f^{*}$ the Fenchel conjugate of $f$ (Borwein and Lewis 2006), defined by $f^{*}({\boldsymbol{\kappa}})\triangleq\sup_{{\mathbf{z}}}[{\mathbf{z}}^{\top}{\boldsymbol{\kappa}}-f({\mathbf{z}})]$. The duality gap for problem (2) can be derived from standard Fenchel duality arguments (Borwein and Lewis 2006) and it is equal to

Therefore, evaluating the duality gap requires to compute efficiently $\Omega^{*}$ in order to find a feasible dual variable ${\boldsymbol{\kappa}}$ (the gap is otherwise equal to $+\infty$ and becomes non-informative). This is equivalent to solving another network flow problem, based on the following variational formulation:

In the network problem associated with (9), the capacities on the arcs $(s,g)$, $g\in{\mathcal{G}}$, are set to $\tau\eta_{g}$, and the capacities on the arcs $(j,t)$, $j$ in $\{1,\ldots,p\}$, are fixed to ${\boldsymbol{\kappa}}_{j}$. Solving problem (9) amounts to finding the smallest value of $\tau$, such that there exists a flow saturating all the capacities ${\boldsymbol{\kappa}}_{j}$ on the arcs leading to the sink $t$. Equation (9) and Algorithm 2 are proven to be correct in Appendix B.

Optimization with Proximal Splitting Methods

The issue of overlapping groups is removed, but new constraints are added, and as in Section 3, the method introduces additional variables which induce a memory cost of $O(\sum_{g\in{\mathcal{G}}}|g|)$.

Minimize ${\mathcal{L}}$ with respect to ${\mathbf{w}}$, keeping the other variables fixed.

Minimize ${\mathcal{L}}$ with respect to the ${\mathbf{z}}^{g}$’s, keeping the other variables fixed. The solution can be obtained in closed form: for all $g$ in $\mathcal{G}$, ${\mathbf{z}}^{g}\leftarrow\text{prox}_{\frac{\lambda\eta_{g}}{\gamma}\|.\|}[{\mathbf{w}}_{g}-\frac{1}{\gamma}\nu^{g}]$.

Take a gradient ascent step on ${\mathcal{L}}$ with respect to the $\nu^{g}$’s: $\nu^{g}\leftarrow\nu^{g}+\gamma({\mathbf{z}}^{g}-{\mathbf{w}}_{g})$.

Such a procedure is guaranteed to converge to the desired solution for all value of $\gamma>0$ (however, tuning $\gamma$ can greatly influence the convergence speed), but solving efficiently step 1 can be difficult. To cope with this issue, we propose two variations exploiting assumptions (A) and (B).

1.2 Dealing with the Design Matrix

Applications and Experiments

In this section, we present various experiments demonstrating the applicability and the benefits of our methods for solving large-scale sparse and structured regularized problems.

In our experiments, the regularization parameter $\lambda$ is chosen to achieve the same level of sparsity ($20\%$). For SG, ADMM and Lin-ADMM, some parameters are optimized to provide the lowest value of the objective function after $1\,000$ iterations of the respective algorithms. For SG, we take the step size to be equal to $a/(k+b)$, where $k$ is the iteration number, and $(a,b)$ are the pair of parameters selected in $\{10^{-3},\dots,10\}\!\times\!\{10^{2},10^{3},10^{4}\}$. Note that a step size of the form $a/(\sqrt{t}+b)$ is also commonly used in subgradient descent algorithms. In the context of hierarchical norms, both choices have led to similar results (Jenatton et al. 2011). The parameter $\gamma$ for ADMM is selected in $\{10^{-2},\ldots,10^{2}\}$. The parameters $(\gamma,\delta)$ for Lin-ADMM are selected in $\{10^{-2},\ldots,10^{2}\}\times\{10^{-1},\ldots,10^{8}\}$. For interior point methods, since problem (2) can be cast either as a quadratic (QP) or as a conic program (CP), we show in Figure 2 the results for both formulations. On three problems of different sizes, with $(n,p)\in\{(100,10^{3}),(1024,10^{4}),(1024,10^{5})\}$, our algorithms ProxFlow, ADMM and Lin-ADMM compare favorably with the other methods, (see Figure 2), except for ADMM in the large-scale setting which yields an objective function value similar to that of SG after $10^{4}$ seconds. Among ProxFlow, ADMM and Lin-ADMM, ProxFlow is consistently better than Lin-ADMM, which is itself better than ADMM. Note that for the small scale problem, the performance of ProxFlow and Lin-ADMM is similar. In addition, note that QP, CP, SG, ADMM and Lin-ADMM do not obtain sparse solutions, whereas ProxFlow does. To reduce the computational cost of this experiment, the curves reported are the results of one single run. Similar types of experiments with several runs have shown very small variability (Bach et al. 2011).

2 Wavelet Denoising with Structured Sparsity

We consider Daubechies3 wavelets (Mallat 1999, see) for the matrix ${\mathbf{X}}$, use $12$ classical standard test images, These images are used in classical image denoising benchmarks. See Mairal et al. 2009. and generate noisy versions of them corrupted by a white Gaussian noise of variance $\sigma^{2}$. For each image, we test several values of $\lambda=2^{\frac{i}{4}}\sigma\sqrt{\log{p}}$, with $i$ taken in the range $\{-15,-14,\dots,15\}$. We then keep the parameter $\lambda$ giving the best reconstruction error on average on the $12$ images. The factor $\sigma\sqrt{\log{p}}$ is a classical heuristic for choosing a reasonable regularization parameter (Mallat 1999, see). We provide reconstruction results in terms of PSNR in Table 1. Denoting by MSE the mean-squared-error for images whose intensities are between $0$ and $255$, the PSNR is defined as $\textrm{PSNR}=10\log_{10}(255^{2}/\textrm{MSE})$ and is measured in dB. A gain of $1$dB reduces the MSE by approximately $20\%$. Unlike Jenatton et al. 2011, who set all the weights $\eta_{g}$ in $\Omega$ equal to one, we tried exponential weights of the form $\eta_{g}=\rho^{k}$, with $k$ being the depth of the group in the wavelet tree, and $\rho$ is taken in $\{0.25,0.5,1,2,4\}$. As for $\lambda$, the value providing the best reconstruction is kept. The wavelet transforms in our experiments are computed with the matlabPyrTools software. http://www.cns.nyu.edu/$\sim$eero/steerpyr/. Interestingly, we observe in Table 1 that the results obtained with $\Omega_{\text{grid}}$ are significantly better than those obtained with $\Omega_{\text{tree}}$, meaning that encouraging spatial consistency in wavelet coefficients is more effective than using a hierarchical coding.

We also note that our approach is relatively fast, despite the high dimension of the problem. Solving exactly the proximal problem with $\Omega_{\text{grid}}$ for an image with $p=512\times 512=262\,144$ pixels (and therefore approximately the same number of groups) takes approximately $\approx 4-6$ seconds on a single core of a 3.07GHz CPU.

3 CUR-like Matrix Factorization

In Mahoney and Drineas 2009, CUR decompositions are computed by a sampling procedure based on the singular value decomposition of ${\mathbf{X}}$. In a recent work, Bien et al. 2010 have shown that partial CUR decompositions, i.e., the selection of either rows or columns of ${\mathbf{X}}$, can be obtained by solving a convex program with a group-Lasso penalty. We propose to extend this approach to the simultaneous selection of both rows and columns of ${\mathbf{X}}$, with the following convex problem:

Problem (12) has a smooth convex data-fitting term and brings into play a sparsity-inducing norm with overlapping groups of variables (the rows and the columns of ${\mathbf{W}}$). As a result, it is a particular instance of problem (2) that can therefore be handled with the optimization tools introduced in this paper. We now compare the performance of the sampling procedure from Mahoney and Drineas 2009 with our proposed sparsity-based approach. To this end, we consider the four gene-expression datasets 9_Tumors, Brain_Tumors1, Leukemia1 and SRBCT, with respective dimensions $(n,p)\in\{(60,5727),(90,5921),(72,5328),(83,2309)\}$. The datasets are freely available at http://www.gems-system.org/. In the sequel, the matrix ${\mathbf{X}}$ is normalized to have unit Frobenius-norm while each of its columns is centered. To begin with, we run our approach More precisely, since the penalties in problem (12) shrink the coefficients of ${\mathbf{W}}$, we follow a two-step procedure: We first run our approach to determine the sets of nonzero rows and columns, and then compute ${\mathbf{W}}_{\text{I}\,\text{J}}={\mathbf{C}}^{+}{\mathbf{X}}{\mathbf{R}}^{+}$. over a grid of values for $\lambda_{\text{row}}$ and $\lambda_{\text{col}}$ in order to obtain solutions with different sparsity levels, i.e., ranging from $|\text{I}|=p$ and $|\text{J}|=n$ down to $|\text{I}|=|\text{J}|=0$. For each pair of values $[|\text{I}|,|\text{J}|]$, we then apply the sampling procedure from Mahoney and Drineas 2009. Finally, the variance explained by the CUR decompositions is reported in Figure 3 for both methods. Since the sampling approach involves some randomness, we show the average and standard deviation of the results based on five initializations. The conclusions we can draw from the experiments match the ones already reported in Bien et al. 2010 for the partial CUR decomposition. We can indeed see that both schemes perform similarly. However, our approach has the advantage not to be randomized, which can be less disconcerting in the practical perspective of analyzing a single run of the algorithm. It is finally worth being mentioned that the convex approach we develop here is flexible and can be extended in different ways. For instance, we can imagine to add further low-rank/sparsity constraints on ${\mathbf{W}}$ thanks to sparsity-promoting convex regularizations.

4 Background Subtraction

5 Topographic Dictionary Learning

6 Multi-Task Learning of Hierarchical Structures

As mentioned in the previous section, Jenatton et al. 2010a have recently proposed to use a hierarchical structured norm to learn dictionaries of natural image patches. In Jenatton et al. 2010a, the dictionary elements are embedded in a predefined tree $\mathcal{T}$, via a particular instance of the structured norm $\Omega$, which we refer to it as $\Omega_{\text{tree}}$, and call ${\mathcal{G}}$ the underlying set of groups. In this case, using the same notation as in Section 5.5, each signal ${\mathbf{y}}^{i}$ admits a sparse decomposition in the form of a subtree of dictionary elements.

Inspired by ideas from multi-task learning Obozinski et al. 2010, we propose to learn the tree structure $\mathcal{T}$ by pruning irrelevant parts of a larger initial tree $\mathcal{T}_{0}$. We achieve this by using an additional regularization term $\Omega_{\text{joint}}$ across the different decompositions, so that subtrees of $\mathcal{T}_{0}$ will simultaneously be removed for all signals ${\mathbf{y}}^{i}$. With the notation from Section 5.5, the approach of Jenatton et al. 2010a is then extended by the following formulation:

Conclusion

In addition to making it possible to resort to accelerated gradient methods, an efficient computation of the proximal operator offers more generally a certain modularity, in that it can be used as a building-block for other optimization problems. A case in point is dictionary learning where proximal problems come up and have to be solved repeatedly in an inner-loop. Interesting future work includes the computation of other structured norms such as the one introduced in Jacob et al. 2009, or total-variation based penalties, whose proximal operators are also based on minimum cost flow problems (Chambolle and Darbon 2009). Several experiments demonstrate that our algorithm can be applied to a wide class of learning problems, which have not been addressed before with convex sparse methods.

Appendix A Equivalence to Canonical Graphs

Formally, the notion of equivalence between graphs can be summarized by the following lemma:

Let $G=(V,E,s,t)$ be the canonical graph corresponding to a group structure ${\mathcal{G}}$. Let $G^{\prime}=(V,E^{\prime},s,t)$ be a graph sharing the same set of vertices, source and sink as $G$, but with a different arc set $E^{\prime}$. We say that $G^{\prime}$ is equivalent to $G$ if and only if the following conditions hold:

Arcs of $E^{\prime}$ outgoing from the source are the same as in $E$, with the same costs and capacities.

Arcs of $E^{\prime}$ going to the sink are the same as in $E$, with the same costs and capacities.

For every arc $(g,j)$ in $E$, with $(g,j)$ in $V_{gr}\times V_{u}$, there exists a unique path in $E^{\prime}$ from $g$ to $j$ with zero costs and infinite capacities on every arc of the path.

Conversely, if there exists a path in $E^{\prime}$ between a vertex $g$ in $V_{gr}$ and a vertex $j$ in $V_{u}$, then there exists an arc $(g,j)$ in $E$.

Then, the cost of the optimal min-cost flow on $G$ and $G^{\prime}$ are the same. Moreover, the values of the optimal flow on the arcs $(j,t)$, $j$ in $V_{u}$, are the same on $G$ and $G^{\prime}$.

We first notice that on both $G$ and $G^{\prime}$, the cost of a flow on the graph only depends on the flow on the arcs $(j,t)$, $j$ in $V_{u}$, which we have denoted by ${\boldsymbol{\overline{\xi}}}$ in $E$.

We will prove that finding a feasible flow $\pi$ on $G$ with a cost $c(\pi)$ is equivalent to finding a feasible flow $\pi^{\prime}$ on $G^{\prime}$ with the same cost $c(\pi)=c(\pi^{\prime})$. We now use the concept of path flow, which is a flow vector in $G$ carrying the same positive value on every arc of a directed path between two nodes of $G$. It intuitively corresponds to sending a positive amount of flow along a path of the graph.

According to the definition of graph equivalence introduced in the Lemma, it is easy to show that there is a bijection between the arcs in $E$, and the paths in $E^{\prime}$ with positive capacities on every arc. Given now a feasible flow $\pi$ in $G$, we build a feasible flow $\pi^{\prime}$ on $G^{\prime}$ which is a sum of path flows. More precisely, for every arc $a$ in $E$, we consider its equivalent path in $E^{\prime}$, with a path flow carrying the same amount of flow as $a$. Therefore, each arc $a^{\prime}$ in $E^{\prime}$ has a total amount of flow that is equal to the sum of the flows carried by the path flows going over $a^{\prime}$. It is also easy to show that this construction builds a flow on $G^{\prime}$ (capacity and conservation constraints are satisfied) and that this flow $\pi^{\prime}$ has the same cost as $\pi$, that is, $c(\pi)=c(\pi^{\prime})$.

Conversely, given a flow $\pi^{\prime}$ on $G^{\prime}$, we use a classical path flow decomposition (see Bertsekas 1998, Proposition 1.1), saying that there exists a decomposition of $\pi^{\prime}$ as a sum of path flows in $E^{\prime}$. Using the bijection described above, we know that each path in the previous sums corresponds to a unique arc in $E$. We now build a flow $\pi$ in $G$, by associating to each path flow in the decomposition of $\pi^{\prime}$, an arc in $E$ carrying the same amount of flow. The flow of every other arc in $E$ is set to zero. It is also easy to show that this builds a valid flow in $G$ that has the same cost as $\pi^{\prime}$. ∎

Appendix B Convergence Analysis

We show in this section the correctness of Algorithm 1 for computing the proximal operator, and of Algorithm 2 for computing the dual norm $\Omega^{\star}$.

We first prove that our algorithm converges and that it finds the optimal solution of the proximal problem. This requires that we introduce the optimality conditions for problem (6) derived from Jenatton et al. 2010a; Jenatton et al. 2011 since our convergence proof essentially checks that these conditions are satisfied upon termination of the algorithm.

The primal-dual variables $({\mathbf{w}},{\boldsymbol{\xi}})$ are respectively solutions of the primal (4) and dual problems (6) if and only if the dual variable ${\boldsymbol{\xi}}$ is feasible for the problem (6) and

Before proving the convergence and correctness of our algorithm, we also recall classical properties of the min capacity cuts, which we intensively use in the proofs of this paper. The procedure computeFlow of our algorithm finds a minimum $(s,t)$-cut of a graph $G=(V,E,s,t)$, dividing the set $V$ into two disjoint parts $V^{+}$ and $V^{-}$. $V^{+}$ is by construction the sets of nodes in $V$ such that there exists a non-saturating path from $s$ to $V$, while all the paths from $s$ to $V^{-}$ are saturated. Conversely, arcs from $V^{+}$ to $t$ are all saturated, whereas there can be non-saturated arcs from $V^{-}$ to $t$. Moreover, the following properties, which are illustrated on Figure 7, hold

There is no arc going from $V^{+}$ to $V^{-}$. Otherwise the value of the cut would be infinite (arcs inside $V$ have infinite capacity by construction of our graph).

There is no flow going from $V^{-}$ to $V^{+}$ (Bertsekas 1998, see).

The cut goes through all arcs going from $V^{+}$ to $t$, and all arcs going from $s$ to $V^{-}$.

Recall that we assume (cf. Section 3.3) that the scalars ${\mathbf{u}}_{j}$ are all non negative, and that we add non-negativity constraints on ${\boldsymbol{\xi}}$. With the optimality conditions of Lemma 5 in hand, we can show our first convergence result.

Algorithm 1 converges in a finite and polynomial number of operations.

Suppose for instance that $V^{-}\!\!=\emptyset$. In this case, the capacity of the min-cut is equal to $\sum_{j\in V_{u}}{\boldsymbol{\gamma}}_{j}$, and the value of the max-flow is $\sum_{j\in V_{u}}{\boldsymbol{\overline{\xi}}}_{j}$. Using the classical max-flow/min-cut theorem Ford and Fulkerson 1956, we have equality between these two terms. Since, by definition of both ${\boldsymbol{\gamma}}$ and ${\boldsymbol{\overline{\xi}}}$, we have for all $j$ in $V_{u}$, ${\boldsymbol{\overline{\xi}}}_{j}\leq{\boldsymbol{\gamma}}_{j}$, we obtain a contradiction with the existence of $j$ in $V_{u}$ such that ${\boldsymbol{\overline{\xi}}}_{j}\neq{\boldsymbol{\gamma}}_{j}$.

Conversely, suppose now that $V^{+}\!\!=\emptyset$. Then, the value of the max-flow is still $\sum_{j\in V_{u}}{\boldsymbol{\overline{\xi}}}_{j}$, and the value of the min-cut is $\lambda\sum_{g\in V_{gr}}\eta_{g}$. Using again the max-flow/min-cut theorem, we have that $\sum_{j\in V_{u}}{\boldsymbol{\overline{\xi}}}_{j}=\lambda\sum_{g\in V_{gr}}\eta_{g}$. Moreover, by definition of ${\boldsymbol{\gamma}}$, we also have $\sum_{j\in V_{u}}{\boldsymbol{\overline{\xi}}}_{j}\leq\sum_{j\in V_{u}}{\boldsymbol{\gamma}}_{j}\leq\lambda\sum_{g\in V_{gr}}\eta_{g}$, leading to a contradiction with the existence of $j$ in $V_{u}$ satisfying ${\boldsymbol{\overline{\xi}}}_{j}\neq{\boldsymbol{\gamma}}_{j}$. We remind the reader of the fact that such a $j\in V_{u}$ exists since the cut is only computed when the current estimate ${\boldsymbol{\overline{\xi}}}$ is not optimal yet. This proof holds for any graph that is equivalent to the canonical one. ∎

After proving the convergence, we prove that the algorithm is correct with the next proposition.

Algorithm 1 solves the proximal problem of Eq. (4).

For a group structure ${\mathcal{G}}$, we first prove the correctness of our algorithm if the graph used is its associated canonical graph that we denote $G_{0}=(V_{0},E_{0},s,t)$. We proceed by induction on the number of nodes of the graph. The induction hypothesis ${\mathcal{H}}(k)$ is the following:

For all canonical graphs $G=(V=V_{u}\cup V_{gr},E,s,t)$ associated with a group structure ${\mathcal{G}}_{V}$ with weights $(\eta_{g})_{g\in{\mathcal{G}}_{V}}$ such that $|V|\leq k$, computeFlow$(V,E)$ solves the following optimization problem:

Since ${\mathcal{G}}_{V_{0}}={\mathcal{G}}$, it is sufficient to show that ${\mathcal{H}}(|V_{0}|)$ to prove the proposition.

We initialize the induction by ${\mathcal{H}}(2)$, corresponding to the simplest canonical graph, for which $|V_{gr}|=|V_{u}|=1$). Simple algebra shows that ${\mathcal{H}}(2)$ is indeed correct.

The algorithm then computes a max-flow, using the scalars ${\boldsymbol{\gamma}}_{j}$ as capacities, and we now have two possible situations:

If ${\boldsymbol{\overline{\xi}}}_{j}={\boldsymbol{\gamma}}_{j}$ for all $j$ in $V_{u}$, the algorithm stops; we write ${\mathbf{w}}_{j}={\mathbf{u}}_{j}-{\boldsymbol{\overline{\xi}}}_{j}$ for $j$ in $V_{u}$, and using Eq. (17), we obtain

Since all the quantities in the previous sum are positive, this can only hold if for all $g\in V_{gr}$,

Moreover, by definition of the max flow and the optimality conditions, we have

By Lemma 5, we have shown that the problem (16) is solved.

Let us now consider the case where there exists $j$ in $V_{u}$ such that ${\boldsymbol{\overline{\xi}}}_{j}\neq{\boldsymbol{\gamma}}_{j}$. The algorithm splits the vertex set $V$ into two parts $V^{+}$ and $V^{-}$, which we have proven to be non-empty in the proof of Proposition 6. The next step of the algorithm removes all edges between $V^{+}$ and $V^{-}$ (see Figure 7). Processing $(V^{+},E^{+})$ and $(V^{-},E^{-})$ independently, it updates the value of the flow matrix ${\boldsymbol{\xi}}^{g}_{j},\ j\in V_{u},\ g\in V_{gr}$, and the corresponding flow vector ${\boldsymbol{\overline{\xi}}}_{j},\ j\in V_{u}$. As for $V$, we denote by $V^{+}_{u}\triangleq V^{+}\cap V_{u}$, $V^{-}_{u}\triangleq V^{-}\cap V_{u}$ and $V^{+}_{gr}\triangleq V^{+}\cap V_{gr}$, $V^{-}_{gr}\triangleq V^{-}\cap V_{gr}$.

Then, we notice that $(V^{+},E^{+},s,t)$ and $(V^{-},E^{-},s,t)$ are respective canonical graphs for the group structures ${\mathcal{G}}_{V^{+}}\triangleq\{g\cap V_{u}^{+}\mid g\in V_{gr}\}$, and ${\mathcal{G}}_{V^{-}}\triangleq\{g\cap V_{u}^{-}\mid g\in V_{gr}\}$.

Writing ${\mathbf{w}}_{j}={\mathbf{u}}_{j}-{\boldsymbol{\overline{\xi}}}_{j}$ for $j$ in $V_{u}$, and using the induction hypotheses ${\mathcal{H}}(|V^{+}|)$ and ${\mathcal{H}}(|V^{-}|)$, we now have the following optimality conditions deriving from Lemma 5 applied on Eq. (16) respectively for the graphs $(V^{+},E^{+})$ and $(V^{-},E^{-})$:

We will now combine Eq. (19) and Eq. (20) into optimality conditions for Eq. (16). We first notice that $g\cap V_{u}^{+}=g$ since there are no arcs between $V^{+}$ and $V^{-}$ in $E$ (see the properties of the cuts discussed before this proposition). It is therefore possible to replace $g^{\prime}$ by $g$ in Eq. (19). We will show that it is possible to do the same in Eq. (20), so that combining these two equations yield the optimality conditions of Eq. (16).

More precisely, we will show that for all $g\in V_{gr}^{-}$ and $j\in g\cap V_{u}^{+}$, $|{\mathbf{w}}_{j}|\leq\max_{l\in g\cap V_{u}^{-}}|{\mathbf{w}}_{l}|$, in which case $g^{\prime}$ can be replaced by $g$ in Eq. (20). This result is relatively intuitive: $(s,V^{+})$ and $(V^{-},t)$ being an $(s,t)$-cut, all arcs between $s$ and $V^{-}$ are saturated, while there are unsaturated arcs between $s$ and $V^{+}$; one therefore expects the residuals ${\mathbf{u}}_{j}-{\boldsymbol{\overline{\xi}}}_{j}$ to decrease on the $V^{+}$ side, while increasing on the $V^{-}$ side. The proof is nonetheless a bit technical.

Let us show first that for all $g$ in $V_{gr}^{+}$, $\left\|{\mathbf{w}}_{g}\right\|_{\infty}\leq\max_{j\in V_{u}}|{\mathbf{u}}_{j}-{\boldsymbol{\gamma}}_{j}|$. We split the set $V^{+}$ into disjoint parts:

As previously, we denote $V^{+-}\!\triangleq V_{u}^{+-}\!\cup V_{gr}^{+-}$ and $V^{++}\triangleq\!\!V_{u}^{++}\cup V_{gr}^{++}$. We want to show that $V_{gr}^{+-}$ is necessarily empty. We reason by contradiction and assume that $V_{gr}^{+-}\neq\varnothing$.

According to the definition of the different sets above, we observe that no arcs are going from $V^{++}$ to $V^{+-}$, that is, for all $g$ in $V_{gr}^{++}$, $g\cap V_{u}^{+-}=\varnothing$. We observe as well that the flow from $V_{gr}^{+-}$ to $V_{u}^{++}$ is the null flow, because optimality conditions (19) imply that for a group $g$ only nodes $j\in g$ such that ${\mathbf{w}}_{j}=\|{\mathbf{w}}_{g}\|_{\infty}$ receive some flow, which excludes nodes in $V_{u}^{++}$ provided $V_{gr}^{+-}\neq\varnothing$; Combining this fact and the inequality $\sum_{g\in V_{gr}^{+}}\lambda\eta_{g}\geq\sum_{j\in V_{u}^{+}}{\boldsymbol{\gamma}}_{j}$ (which is a direct consequence of the minimum $(s,t)$-cut), we have as well

Let $j\in V_{u}^{+-}$, if ${\boldsymbol{\overline{\xi}}}_{j}\neq 0$ then for some $g\in V_{gr}^{+-}$ such that $j$ receives some flow from $g$, which from the optimality conditions (19) implies ${\mathbf{w}}_{j}=\|{\mathbf{w}}_{g}\|_{\infty}$; by definition of $V_{gr}^{+-}$, $\|{\mathbf{w}}_{g}\|_{\infty}>{\mathbf{u}}_{j}-{\boldsymbol{\gamma}}_{j}$. But since at the optimum, ${\mathbf{w}}_{j}={\mathbf{u}}_{j}-{\boldsymbol{\overline{\xi}}}_{j}$, this implies that ${\boldsymbol{\overline{\xi}}}_{j}<{\boldsymbol{\gamma}}_{j}$, and in turn that $\sum_{j\in V_{u}^{+-}}{\boldsymbol{\overline{\xi}}}_{j}=\lambda\sum_{g\in V_{gr}^{+-}}\eta_{g}$. Finally,

We now have that for all $g$ in $V_{gr}^{+}$, $\left\|{\mathbf{w}}_{g}\right\|_{\infty}\leq\max_{j\in V_{u}}|{\mathbf{u}}_{j}-{\boldsymbol{\gamma}}_{j}|$. The proof showing that for all $g$ in $V_{gr}^{-}$, $\left\|{\mathbf{w}}_{g}\right\|_{\infty}\geq\max_{j\in V_{u}}|{\mathbf{u}}_{j}-{\boldsymbol{\gamma}}_{j}|,$ uses the same kind of decomposition for $V^{-}$, and follows along similar arguments. We will therefore not detail it.

To summarize, we have shown that for all $g\in V_{gr}^{-}$ and $j\in g\cap V_{u}^{+}$, $|{\mathbf{w}}_{j}|\leq\max_{l\in g\cap V_{u}^{-}}|{\mathbf{w}}_{l}|$. Since there is no flow from $V^{-}$ to $V^{+}$, i.e., ${\boldsymbol{\xi}}_{j}^{g}=0$ for $g$ in $V_{gr}^{-}$ and $j$ in $V_{u}^{+}$, we can now replace the definition of $g^{\prime}$ in Eq. (20) by $g^{\prime}\triangleq g\cap V_{u}$, the combination of Eq. (19) and Eq. (20) gives us optimality conditions for Eq. (16).

The proposition being proved for the canonical graph, we extend it now for an equivalent graph in the sense of Lemma 4. First, we observe that the algorithm gives the same values of ${\boldsymbol{\gamma}}$ for two equivalent graphs. Then, it is easy to see that the value ${\boldsymbol{\overline{\xi}}}$ given by the max-flow, and the chosen $(s,t)$-cut is the same, which is enough to conclude that the algorithm performs the same steps for two equivalent graphs. ∎

B.2 Computation of the Dual Norm Ω⋆\Omega^{\star}

As for the proximal operator, the computation of dual norm $\Omega^{*}$ can itself be shown to solve another network flow problem, based on the following variational formulation, which extends a previous result from Jenatton et al. 2009:

By definition of $\Omega^{*}({\boldsymbol{\kappa}})$, we have

with the additional $|{\mathcal{G}}|$ conic constraints $\|{\mathbf{z}}_{g}\|_{\infty}\leq\alpha_{g}$. 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 and flipping the sign of ${\boldsymbol{\xi}}$, the dual problem then reduces to

We now prove that Algorithm 2 is correct.

Algorithm 2 computes the value of the dual norm of Eq. (9) in a finite and polynomial number of operations.

The convergence of the algorithm only requires to show that the cardinality of $V$ in the different calls of the function computeFlow strictly decreases. Similar arguments to those used in the proof of Proposition 6 can show that each part of the cuts $(V^{+},V^{-})$ are both non-empty. The algorithm thus requires a finite number of calls to a max-flow algorithm and converges in a finite and polynomial number of operations.

Let us now prove that the algorithm is correct for a canonical graph. We proceed again by induction on the number of nodes of the graph. More precisely, we consider the induction hypothesis ${\mathcal{H}}^{\prime}(k)$ defined as:

for all canonical graphs $G=(V,E,s,t)$ associated with a group structure ${\mathcal{G}}_{V}$ and such that $|V|\leq k$, dualNormAux$(V=V_{u}\cup V_{gr},E)$ solves the following optimization problem:

We first initialize the induction by ${\mathcal{H}}(2)$ (i.e., with the simplest canonical graph, such that $|V_{gr}|=|V_{u}|=1$). Simple algebra shows that ${\mathcal{H}}(2)$ is indeed correct.

We next consider a canonical graph $G=(V,E,s,t)$ such that $|V|=k$, and suppose that ${\mathcal{H}}^{\prime}(k-1)$ is true. After the max-flow step, we have two possible cases to discuss:

If ${\boldsymbol{\overline{\xi}}}_{j}={\boldsymbol{\gamma}}_{j}$ for all $j$ in $V_{u}$, the algorithm stops. We know that any scalar $\tau$ such that the constraints of Eq. (21) are all satisfied necessarily verifies $\sum_{g\in V_{gr}}\tau\eta_{g}\geq\sum_{j\in V_{u}}{\boldsymbol{\kappa}}_{j}$. We have indeed that $\sum_{g\in V_{gr}}\tau\eta_{g}$ is the value of an $(s,t)$-cut in the graph, and $\sum_{j\in V_{u}}{\boldsymbol{\kappa}}_{j}$ is the value of the max-flow, and the inequality follows from the max-flow/min-cut theorem Ford and Fulkerson 1956. This gives a lower-bound on $\tau$. Since this bound is reached, $\tau$ is necessarily optimal.

We now consider the case where there exists $j$ in $V_{u}$ such that ${\boldsymbol{\overline{\xi}}}_{j}\neq{\boldsymbol{\kappa}}_{j}$, meaning that for the given value of $\tau$, the constraint set of Eq. (21) is not feasible for ${\boldsymbol{\xi}}$, and that the value of $\tau$ should necessarily increase. The algorithm splits the vertex set $V$ into two non-empty parts $V^{+}$ and $V^{-}$ and we remark that there are no arcs going from $V^{+}$ to $V^{-}$, and no flow going from $V^{-}$ to $V^{+}$. Since the arcs going from $s$ to $V^{-}$ are saturated, we have that $\sum_{g\in V_{gr}^{-}}\tau\eta_{g}\leq\sum_{j\in V_{u}^{-}}{\boldsymbol{\kappa}}_{j}$. Let us now consider $\tau^{\star}$ the solution of Eq. (21). Using the induction hypothesis ${\mathcal{H}}^{\prime}(|V^{-}|)$, the algorithm computes a new value $\tau^{\prime}$ that solves Eq. (21) when replacing $V$ by $V^{-}$ and this new value satisfies the following inequality $\sum_{g\in V_{gr}^{-}}\tau^{\prime}\eta_{g}\geq\sum_{j\in V_{u}^{-}}{\boldsymbol{\kappa}}_{j}$. The value of $\tau^{\prime}$ has therefore increased and the updated flow ${\boldsymbol{\xi}}$ now satisfies the constraints of Eq. (21) and therefore $\tau^{\prime}\geq\tau^{\star}$. Since there are no arcs going from $V^{+}$ to $V^{-}$, $\tau^{\star}$ is feasible for Eq. (21) when replacing $V$ by $V^{-}$ and we have that $\tau^{\star}\geq\tau^{\prime}$ and then $\tau^{\prime}=\tau^{\star}$.

To prove that the result holds for any equivalent graph, similar arguments to those used in the proof of Proposition 6 can be exploited, showing that the algorithm computes the same values of $\tau$ and same $(s,t)$-cuts at each step. ∎

Appendix C Algorithm FISTA with duality gap

In this section, we describe in details the algorithm FISTA Beck and Teboulle 2009 when applied to solve problem (2), with a duality gap as the stopping criterion. The algorithm, as implemented in the experiments, is summarized in Algorithm 3.

in place of problem (2). Based on Fenchel duality arguments Borwein and Lewis 2006,

is a duality gap for problem (22), where $f^{*}({\boldsymbol{\kappa}})\triangleq\sup_{{\mathbf{z}}}[{\mathbf{z}}^{\top}{\boldsymbol{\kappa}}-f({\mathbf{z}})]$ is the Fenchel conjugate of $f$ (Borwein and Lewis 2006). Given a primal variable ${\mathbf{w}}$, a good dual candidate ${\boldsymbol{\kappa}}$ can be obtained by looking at the conditions that have to be satisfied by the pair $({\mathbf{w}},{\boldsymbol{\kappa}})$ at optimality Borwein and Lewis 2006. In particular, the dual variable ${\boldsymbol{\kappa}}$ is chosen to be

In our experiment, we choose the line-search parameter $\nu$ to be equal to $1.5$.

Appendix D Speed comparison of Algorithm 1 with parametric max-flow algorithms

As shown by Hochbaum and Hong 1995, min-cost flow problems, and in particular, the dual problem of (4), can be reduced to a specific parametric max-flow problem. We thus compare our approach (ProxFlow) with the efficient parametric max-flow algorithm proposed by Gallo et al. 1989 and a simplified version of the latter proposed by Babenko and Goldberg 2006. We refer to these two algorithms as GGT and SIMP respectively. The benchmark is established on the same datasets as those already used in the experimental section of the paper, namely: (1) three datasets built from overcomplete bases of discrete cosine transforms (DCT), with respectively $10^{4},\ 10^{5}$ and $10^{6}$ variables, and (2) images used for the background subtraction task, composed of 57600 pixels. For GGT and SIMP, we use the paraF software which is a C++ parametric max-flow implementation available at http://www.avglab.com/andrew/soft.html. Experiments were conducted on a single-core 2.33 Ghz. We report in the following table the average execution time in seconds of each algorithm for $5$ runs, as well as the statistics of the corresponding problems:

Although we provide the speed comparison for a single value of $\lambda$ (the one used in the corresponding experiments of the paper), we observed that our approach consistently outperforms GGT and SIMP for values of $\lambda$ corresponding to different regularization regimes.

References