Structured sparsity-inducing norms through submodular functions

Introduction

The concept of parsimony is central in many scientific domains. In the context of statistics, signal processing or machine learning, it takes the form of variable or feature selection problems, and is commonly used in two situations: First, to make the model or the prediction more interpretable or cheaper to use, i.e., even if the underlying problem does not admit sparse solutions, one looks for the best sparse approximation. Second, sparsity can also be used given prior knowledge that the model should be sparse. In these two situations, reducing parsimony to finding models with low cardinality turns out to be limiting, and structured parsimony has emerged as a fruitful practical extension, with applications to image processing, text processing or bioinformatics (see, e.g., and Section 4). For example, in , structured sparsity is used to encode prior knowledge regarding network relationship between genes, while in , it is used as an alternative to structured non-parametric Bayesian process based priors for topic models.

This is done for a particular ensemble of set-functions $F$, namely nondecreasing submodular functions. Submodular functions may be seen as the set-function equivalent of convex functions, and exhibit many interesting properties that we review in Section 2—see for a tutorial on submodular analysis and for other applications to machine learning. This paper makes the following contributions:

$-$ By selecting specific submodular functions in Section 4, we recover and give a new interpretation to known norms, such as those based on rank-statistics or grouped norms with potentially overlapping groups , and we define new norms, in particular ones that can be used as non-factorial priors for supervised learning (Section 4). These are illustrated on simulation experiments in Section 7, where they outperform related greedy approaches .

Review of submodular function theory

Throughout this paper, we consider a nondecreasing submodular function $F$ defined on the power set $2^{V}$ of $V=\{1,\dots,p\}$, i.e., such that:

Moreover, we assume (without loss of generality) that $F(\varnothing)=0$. These set-functions are often referred to as polymatroid set-functions or $\beta$-functions . Also, without loss of generality, we may assume that $F$ is strictly positive on singletons, i.e., for all $k\in V$, $F(\{k\})>0$. Indeed, if $F(\{k\})=0$, then by submodularity and monotonicity, if $A\ni k$, $F(A)=F(A\backslash\{k\})$ and thus we can simply consider $V\backslash\{k\}$ instead of $V$.

Instead of solving a linear program with $p+2^{p}$ contraints, a solution $s$ may be obtained by the following “greedy algorithm”: order the components of $w$ in decreasing order $w_{j_{1}}\geqslant\dots\geqslant w_{j_{p}}$, and then take for all $k\in\{1,\dots,p\}$, $s_{j_{k}}=F(\{j_{1},\dots,j_{k}\})-F(\{j_{1},\dots,j_{k-1}\}).$

Stable sets. A set $A$ is said stable if it cannot be augmented without increasing $F$, i.e., if for all sets $B\supset A$, $B\neq A\Rightarrow F(B)>F(A)$. If $F$ is strictly increasing, then all sets are stable. Stable sets are also sometimes referred to as flat or closed . The set of stable sets is closed by intersection , and will correspond to the set of allowed sparsity patterns (see Section 6.2). For the cardinality function, all sets are stable.

Submodular functions are particularly interesting because they can be minimized in polynomial time. In this paragraph, we consider a non-monotonic submodular function $G$ (otherwise finding the minimum is trivial). Most algorithms for minimizing submodular functions rely on the following strong duality principle :

Definition and properties of structured norms

Assume that the set-function $F$ is submodular, nondecreasing, and strictly positive for all singletons. Define $\Omega:w\mapsto f(|w|)$. Then:

(iii) the dual norm (see, e.g., ) of $\Omega$ is equal to $\Omega^{\ast}(s)=\max_{A\subset V}\frac{\|s_{A}\|_{1}}{F(A)}=\max_{A\in\mathcal{T}}\frac{\|s_{A}\|_{1}}{F(A)}$.

We provide examples of submodular set-functions and norms in Section 4, where we go from set-functions to norms, and vice-versa. From the definition of the Lovász extension in Eq. (1), we see that $\Omega$ is a polyhedral norm (i.e., its unit ball is a polyhedron). The following proposition gives the set of extreme points of the unit ball (see proof in the appendix and examples in Figure 1):

The extreme points of the unit ball of $\ \Omega$ are the vectors $\frac{1}{F(A)}s$, with $s\in\{-1,0,1\}^{p}$, ${\rm Supp}(s)=A$ and $A$ a stable inseparable set.

This proposition shows, that depending on the number and cardinality of the inseparable stable sets, we can go from $2p$ (only singletons) to $3^{p}-1$ extreme points (all possible sign vectors). We show in Figure 1 examples of balls for $p=2$, as well as sets of extreme points. These extreme points will play a role in concentration inequalities derived in Section 6.

Examples of nondecreasing submodular functions

We consider three main types of submodular functions with potential applications to regularization for supervised learning. Some existing norms are shown to be examples of our frameworks (Section 4.1, Section 4.3), while other novel norms are designed from specific submodular functions (Section 4.2). Other examples of submodular functions, in particular in terms of matroids and entropies, may be found in and could also lead to interesting new norms. Note that set covers, which are common examples of submodular functions are subcases of set-functions defined in Section 4.1 (see, e.g., ).

Hierarchical norms. Hierarchical norms defined on directed acyclic graphs correspond to the set-function $F(A)$ which is the cardinality of the union of ancestors of elements in $A$. These have been applied to bioinformatics , computer vision and topic models .

Another interesting new norm may be defined from the groups in the right side of Figure 2. Indeed, it corresponds to the function $F(A)$ equal to $|A|$ plus the number of intervals of $A$. Note that this also favors contiguous patterns but is not limited to selecting a single interval (like the norm obtained from groups in the left side of Figure 2). Note that it is to be contrasted with the total variation (a.k.a. fused Lasso penalty ), which is a relaxation of the number of jumps in a vector $w$ rather than in its support. In 2D or 3D, this extends to the notion of perimeter and area, but we do not pursue such extensions here.

2 Spectral functions of submatrices

In a frequentist setting, the Mallows $C_{L}$ penalty depends on the degrees of freedom, of the form $\mathop{\rm tr}X_{A}^{\top}X_{A}(X_{A}^{\top}X_{A}+\lambda I)^{-1}$. This is a non-factorial prior but unfortunately it does not lead to a submodular function. In a Bayesian context however, it is shown by that penalties of the form $\log\det(X_{A}^{\top}X_{A}+\lambda I)$ (which lead to submodular functions) correspond to marginal likelihoods associated to the set $A$ and have good behavior when used within a non-convex framework. This highlights the need for non-factorial priors which are sub-linear functions of the eigenvalues of $X_{A}^{\top}X_{A}$, which is exactly what nondecreasing submodular function of submatrices are. We do not pursue the extensive evaluation of non-factorial convex priors in this paper but provide in simulations examples with $F(A)=\mathop{\rm tr}(X_{A}^{\top}X_{A})^{1/2}$ (which is equal to the trace norm of $X_{A}$ ).

3 Functions of cardinality

Convex analysis and optimization

In this section we provide algorithmic tools related to optimization problems based on the regularization by our novel sparsity-inducing norms. Note that since these norms are polyhedral norms with unit balls having potentially an exponential number of vertices or faces, regular linear programming toolboxes may not be used.

Subgradient. From $\Omega(w)=\max_{s\in\mathcal{P}}s^{\top}|w|$ and the greedy algorithm The greedy algorithm to find extreme points of the submodular polyhedron should not be confused with the greedy algorithm (e.g., forward selection) that we consider in Section 7. presented in Section 2, one can easily get in polynomial time one subgradient as one of the maximizers $s$. This allows to use subgradient descent, with, as shown in Figure 4, slow convergence compared to proximal methods.

In the proof, it is shown how a solution for one problem may be obtained from a solution to the other problem. Moreover, any algorithm for minimizing submodular functions allows to get directly the support of the unique solution of the proximal problem and that with a sequence of submodular function minimizations, the full solution may also be obtained. Similar links between convex optimization and minimization of submodular functions have been considered (see, e.g., ). However, these are dedicated to symmetric submodular functions (such as the ones obtained from graph cuts) and are thus not directly applicable to our situation of non-increasing submodular functions.

Sparsity-inducing properties

We study the sparsity-inducing properties of solutions of Eq. (6), i.e., we determine in Section 6.2 which patterns are allowed and in Section 6.3 which sufficient conditions lead to correct estimation. Like recent analysis of sparsity-inducing norms , the analysis provided in this section relies heavily on decomposability properties of our norm $\Omega$.

We can now prove the following decomposition properties, which show that under certain circumstances, we can decompose the norm $\Omega$ on subsets $J$ and their complements:

Given $J\subset V$ and $\Omega_{J}$ and $\Omega^{J}$ defined as above, we have:

2 Sparsity patterns

In this section, we do not make any assumptions regarding the correct specification of the linear model. We show that with probability one, only stable support sets may be obtained (see proof in the appendix). For simplicity, we assume invertibility of $X^{\top}X$, which forbids the high-dimensional situation $p\geqslant n$ we consider in Section 6.3, but we could consider assumptions similar to the ones used in .

3 High-dimensional inference

Let $z$ be a normal variable with covariance matrix $Q$. Let $\mathcal{T}$ be the set of stable inseparable sets. Then $\textstyle P(\Omega^{\ast}(z)>t)\leqslant\sum_{A\in\mathcal{T}}2^{|A|}\exp\big(-\frac{t^{2}F(A)^{2}/2}{1^{\top}Q_{AA}1}\big).$

Experiments

Proximal methods vs. subgradient descent. For the submodular function $F(A)=|A|^{1/2}$ (a simple submodular function beyond the cardinality) we compare three optimization algorithms described in Section 5, subgradient descent and two proximal methods, ISTA and its accelerated version FISTA , for $p=n=1000$, $k=100$ and $\lambda=0.1$. Other settings and other set-functions would lead to similar results than the ones presented in Figure 4: FISTA is faster than ISTA, and much faster than subgradient descent.

Conclusions

We have presented a family of sparsity-inducing norms dedicated to incorporating prior knowledge or structural constraints on the support of linear predictors. We have provided a set of common algorithms and theoretical results, as well as simulations on synthetic examples illustrating the good behavior of these norms. Several avenues are worth investigating: first, we could follow current practice in sparse methods, e.g., by considering related adapted concave penalties to enhance sparsity-inducing norms, or by extending some of the concepts for norms of matrices, with potential applications in matrix factorization or multi-task learning (see, e.g., for application of submodular functions to dictionary learning). Second, links between submodularity and sparsity could be studied further, in particular by considering submodular relaxations of other combinatorial functions, or studying links with other polyhedral norms such as the total variation, which are known to be similarly associated with symmetric submodular set-functions such as graph cuts .

Acknowledgements. This paper was partially supported by the Agence Nationale de la Recherche (MGA Project) and the European Research Council (SIERRA Project). The author would like to thank Edouard Grave, Rodolphe Jenatton, Armand Joulin, Julien Mairal and Guillaume Obozinski for discussions related to this work.

Appendix A Properties of the norm

Thus, for all $w$ such that $\|w\|_{\infty}\leqslant 1$,

Note that $F$ non-increasing implies that $f$ is non-increasing with respect to all of its components. (ii) We have $\displaystyle\Omega(w)=f(|w|)=\max_{s\in\mathcal{P}}s^{\top}|w|=\max_{|s|\in\mathcal{P}}s^{\top}w=\max_{\|s_{A}\|_{1}\leqslant F(A),\ A\subset V}s^{\top}w=\max_{\max_{A\subset V}\frac{\|s_{A}\|_{1}}{F(A)}\leqslant 1}s^{\top}w$, which implies the desired result. Note that the maximization may indeed be limited to the stable inseparable sets $A\in\mathcal{T}$.

A.2 Proof of Proposition 2

We have seen in Section 2 that for $A\in\mathcal{T}$ (set of stable inseparable sets), then $\{x(A)=F(A)\}$ is a face of $\mathcal{P}$ (and those sets are the only ones for which this happens). We get to the desired result by considering potential different signs.

Appendix B Convex optimization results

We first prove an additional result related to decomposition of subdifferentials. Note that the exact subdifferential for the non-zero components of $w$ is rather complicated when $w$ has components with equal magnitude. If this is not the case, i.e., $|w_{j_{1}}|>\cdots>|w_{j_{k}}|>0$, where $k=|J|$, then the subdifferential $\partial\Omega_{J}(w_{J})$ is reduced to a point $s$ such that $s_{j_{k}}=F(\{j_{1},\dots,j_{k}\})-F(\{j_{1},\dots,j_{k-1}\})$. For more details on the subdifferential for nonzero components, see .

Following , without loss of generality, we assume that $z$ has nonnegative components. We have by convex duality (which is applicable here because of Slater’s condition):

where the (unique) optimal $w$ is obtained from the optimal $s$ by $w=z-\lambda s$. $s$ is defined constrained to satisfy $\Omega^{\ast}(s)\leqslant 1$, which is equivalent to $|s|\in\mathcal{P}$. Since $z$ has nonnegative components, the minimum restricted to $|s|\in\mathcal{P}$ is the same as the minimum restricted to $s\in\mathcal{P}$, and also the same as the one restricted to the submodular polyhedron without constraints on positivity, i.e., our problem reduces to $\min_{\forall A\subset V,s(A)\subset F(A)}\|s-z/\lambda\|_{2}^{2}$, which is also equivalent to

Up to the constraints $s(V)=F(V)-\lambda^{-1}z(V)$, this is the minimum-norm point problem for the submodular function $G:A\mapsto F(A)-\lambda^{-1}z(A)$. We can then follow two approaches: the first one is to apply directly the minimum-norm point algorithm to the problem in Eq. (5), which we have followed in simulations. The second approach is to consider the regular minimum point algorithm; we can then follow [12, Lemma 7.4]: if $t$ is the minimum-norm solution for the submodular function $G$, then we can obtain $s$ as $\lambda^{-1}z$ plus the negative part of $t$. From $s$ we then get $w$ through $w=z-\lambda s$.

If another algorithm is used for submodular function minimization, then, following [12, Lemma 7.4], we know which components of the (unique) optimal value $t^{\ast}$ are negative and which of them are equal to zero (which corresponds to zero components of $w^{\ast}$). Then, following , if we add a constant vector with components equal to $\alpha$ to $z$, we may obtain level sets of $w^{\ast}$. With several values of $\alpha$, we can then obtain the full solution $w^{\ast}$. However, the minimum norm point algorithm remains the most efficient and allows to obtain directly the solution of the proximal problem.

Appendix C Sparse estimation

(ii) This is immediate from the expression of the Lovász extension in Eq. (1). Indeed, the order within $J$ and the one within $J^{c}$ do not interact. Note that this case includes cases where we some of the components of $|w_{J}|$ are equal to some of $|w_{J^{c}}|$.

(iii) $\Omega^{J}$ corresponds to the submodular function obtained as the contraction of $F$ by $J$. It is thus a norm as soon as $F^{J}$ is positive on all singletons, which is itself equivalent to the stability of $J$. The equivalence of being a norm with stability of the set $J$ is then straightforward.

C.2 Proof of Proposition 5

If $J$ is not a stable set, then, by Proposition 1, this will implies that there exists $j\in J^{c}$ such that $Q_{jJ}[(Q_{JJ}^{-1}-A_{JJ})r_{J}+b_{J}]-r_{j})=0$, i.e.,

The row vector $Q_{jJ}Q_{JJ}^{-1}X_{J}^{\top}-X_{j}^{\top}-Q_{jJ}A_{JJ}X_{J}^{\top}$ cannot be equal to zero, otherwise,

What remains to be shown is the affine representation of $\hat{w}_{J}$ when the support is given; it is essentially equivalent to showing that the path is piecewise affine, which is not surprising for a polyhedral norm . We use the representation $\Omega_{J}(w_{J})=\max_{z\in B}z^{\top}w_{J}$ where $B$ is the finite set of $z$ such that $|z|$ in an extreme point of the submodular polyhedron associated with $\Omega_{J}$.

Necessary optimality conditions for such the problem in Eq. (6) is the existence of $\eta_{z}\geqslant 0$ (for each $z\in B$) such that (1) $\sum_{z\in B}\eta_{z}=1$, (2) $\eta_{z}=0$ if $z$ is not a maximizer of $\max_{z\in B}z^{\top}w_{J}$, and (3) $w_{J}$ is a minimizer of $\frac{1}{2}w_{J}^{\top}Q_{JJ}w_{J}-r_{J}^{\top}w_{J}+\lambda w^{\top}\sum_{z\in A}\eta_{z}z$, i.e., $Q_{JJ}w_{J}+\lambda\sum_{z\in A}\eta_{z}z=r_{J}$. Moreover, by Carathéodory’s theorem , the number $k$ of non-zero $\eta$ may be taken to be less than $|J|+1$.

It is then a simple linear algebra exercise to show that if $k\leqslant|J|+1$, then $w_{J}$ is of the desired form.

C.3 Proof of Proposition 6

C.4 Proof of Proposition 7

Like for the proof of Proposition 6, we have $\Omega(x)\geqslant\Omega_{J}(x_{J})+\Omega^{J}(x_{J^{c}})\geqslant\Omega_{J}(x_{J})+\rho(J)\Omega_{J^{c}}(x_{J^{c}})\geqslant\rho(J)\Omega(x)$. Thus, if we assume $\Omega^{\ast}(q)\leqslant\lambda\rho(J)/2$, then $\Omega_{J}^{\ast}(q_{J})\leqslant\lambda/2$ and $(\Omega^{J})^{\ast}(q_{J^{c}})\leqslant\lambda/2$. Let $\Delta=\hat{w}-w^{\ast}$.

We follow the proof from by using the decomposition property of the norm $\Omega$. We have, by optimality of $\hat{w}$:

Thus $\Omega^{J}(\Delta_{J^{c}})\leqslant 3\Omega_{J}(\Delta_{J})$, which implies $\Delta^{\top}Q\Delta\geqslant\kappa\|\Delta_{J}\|_{2}^{2}$ (we have assumed a restricted eigenvalue condition). Moreover, we have:

This implies that $\frac{\kappa}{c(J)^{2}}\Omega_{J}(\Delta_{J})^{2}\leqslant\kappa\|\Delta_{J}\|_{2}^{2}\leqslant\Delta^{\top}Q\Delta\leqslant\frac{6\lambda}{\rho(J)}\Omega_{J}(\Delta_{J})$, and thus $\Omega_{J}(\Delta_{J})\leqslant\frac{6c(J)^{2}\lambda}{\kappa\rho(J)}$, which leads to the desired result, given the previous inequalities.

C.5 Proof of Proposition 8

We have $\Omega^{\ast}(z)=\max_{\Omega(w)\leqslant 1}w^{\top}z$; the maximum can be taken over the set of extreme points of the unit ball, which leads to the desired result given Proposition 2.

References