Group-Sparse Model Selection: Hardness and Relaxations

I Introduction

Information in many natural and man-made signals can be exactly represented or well approximated by a sparse set of nonzero coefficients in an appropriate basis . Compressive sensing (CS) exploits this fact to recover signals from their compressive samples, which are dimensionality reducing, non-adaptive random measurements. According to the CS theory, the number of measurements for stable recovery is proportional to the signal sparsity, rather than to its Fourier bandwidth as dictated by the Shannon/Nyquist theorem . Unsurprisingly, the utility of sparse representations also goes well-beyond CS and permeates a lot of fundamental problems in signal processing, machine learning, and theoretical computer science.

Recent results in CS extend the simple sparsity idea to consider more sophisticated structured sparsity models, which describe the interdependency between the nonzero coefficients . There are several compelling reasons for such extensions: The structured sparsity models allow to significantly reduce the number of required measurements for perfect recovery in the noiseless case and be more stable in the presence of noise. Furthermore, they facilitate the interpretation of the signals in terms of the chosen structures, revealing information that could be used to better understand their properties.

An important class of structured sparsity models is based on groups of variables that should either be selected or discarded together . These structures naturally arise in applications such as neuroimaging , gene expression data , bioinformatics and computer vision . For example, in cancer research, the groups might represent genetic pathways that constitute cellular processes. Identifying which processes lead to the development of a tumor can allow biologists to directly target certain groups of genes instead of others . Incorrect identification of the active/inactive groups can thus have a rather dramatic effect on the speed at which cancer therapies are developed.

In this paper, we consider group-based sparsity models, denoted as $\mathfrak{G}$. These structured sparsity models feature collections of groups of variables that could overlap arbitrarily, that is $\mathfrak{G}=\{\mathcal{G}_{1},\ldots,\mathcal{G}_{M}\}$ where each $\mathcal{G}_{j}$ is a subset of the index set $\{1,\ldots,N\}$, with $N$ being the dimensionality of the signal that we model. Arbitrary overlaps mean that we do not restrict the intersection between any two sets from $\mathfrak{G}$.

where $\operatorname{supp}({\bf z})$ is the support of the vector ${\bf z}$. We call such an approximation as G-group-sparse or in short group-sparse. The projection problem is a fundamental step in Model-based Iterative Hard-Thresholding algorithms for solving inverse problems by imposing group structures .

More importantly, we seek to also identify the G-group-support of the approximation $\hat{{\bf x}}$, that is the $G$ groups that constitute its support. We call this the group-sparse model selection problem. The G-group-support of $\hat{{\bf x}}$ allows us to “interpret” the original signal and discover its properties so that we can, for example, target specific groups of genes instead of others or focus more precise imaging techniques on certain brain regions only . In this work, we study under which circumstances we can correctly and tractably identify the $G$-group-support of the approximation of a given signal. In particular, we show that this problem is equivalent to an NP-hard combinatorial problem known as the weighted maximum coverage problem and we propose a novel polynomial time algorithm for finding its solutions for a certain class of group structures.

If the original signal is affected by noise, i.e., if instead of ${\bf x}$, we measure ${\bf z}:={\bf x}+\boldsymbol{\varepsilon}$, where $\boldsymbol{\varepsilon}$ is some random noise, the $G$-group support of $\hat{\bf z}$ may not exactly correspond to the one of $\hat{\bf x}$. Although this is a paramount statistical issue, here we are solely concerned with the computational problem of finding the $G$-group support of a given signal, irrespective of whether it is affected by noise or not, because any group-based interpretation would necessarily require such computation.

Previous work. Recent works in compressive sensing and machine learning with group sparsity have mainly focused on leveraging group structures for lowering the number of samples required for recovering signals . While these results have established the importance of group structures, many of these works have not fully addressed model selection.

For the special case of non-overlapping groups, dubbed the block-sparsity model, the problem of model selection does not present computational difficulties and features a well-understood theory . The first convex relaxations for group-sparse approximation considered only non-overlapping groups. Its extension to overlapping groups , however, selects supports defined as the complement of a union of groups (see also ), which is the opposite of what applications usually require, where groups of variables need to be selected together, instead of discarded.

For overlapping groups, Eldar et al. consider the union of subspaces framework and cast the model selection problem as a block-sparse model selection one by duplicating the variables that belong to overlaps between the groups. Their uniqueness condition [Prop. 1], however, is infeasible for any group structure with overlaps, because it requires that the subspaces intersect only at the origin, while two subspaces defined by two overlapping groups of variables intersect on a subspace of dimension equal to the number of elements in the overlap.

The recently proposed convex relaxations for group-sparse approximations select group-supports that consist of union of groups. However, the group-support recovery conditions in should be taken with care, because they are defined with respect to a particular subset of group-supports and are not general. As we numerically demonstrate in this paper, the group-supports recovered with these methods might be incorrect. Furthermore, the required consistency conditions in are unverifiable a priori. For instance, they require tuning parameters to be known beforehand to obtain the correct group-support.

Huang et al. use coding complexity schemes over sets to encode sparsity structures. They consider linear regression problems where the coding complexity of the support of the solution is constrained to be below a certain value. Inspired by Orthogonal Matching Pursuit, they then propose a greedy algorithm, named StructOMP, that leverages a block-based approximation to the coding complexity. A particular instance of coding schemes, namely graph sparsity, can be used to encode both group and hierarchical sparsity. Their method only returns an approximation to the original discrete problem, as we illustrate via some numerical experiments.

In this work, we take a completely discrete approach and do not rely on relaxations.

Contributions. This paper is an extended version of a prior submission to the IEEE International Symposium on Information Theory (ISIT), 2013. This version contains all the proofs that were previously omitted due to lack of space, refined explanations of the concepts, and provides the full description of the proposed dynamic programming algorithms.

In stark contrast to the existing literature, we take an explicitly discrete approach to identifying group-supports of signals given a budget constraint on the number of groups. This fresh perspective enables us to show that the group-sparse model selection problem is NP-hard: if we can solve the group model selection problem in general, then we can solve any weighted maximum coverage (WMC) problem instance in polynomial time. However, WMC is known to be NP-Hard . Given this connection, we can only hope to characterize a subset of instances which are tractable or find guaranteed and tractable approximations.

We present group structures that lead to computationally tractable problems via dynamic programming. We do so by leveraging a graph-based representation of the groups and exploiting properties of the induced graph. In particular, we present and describe a novel polynomial-time dynamic program that solves the WMC problem for a group structures whose induced graph is a tree or a forest. This result could indeed be of interest by itself.

We identify tractable discrete relaxations of the group-sparse model selection problem that lead to efficient algorithms. Specifically, we relax the constraint on the number of groups into a penalty term and show that if the remaining group constraints satisfy a property related to the concept of total unimodularity , then the relaxed problem can be efficiently solved using linear program solvers. Furthermore, if the graph induced by the group structure is a tree or a forest, we can solve the relaxed problem in linear time by the sum-product algorithm .

We extend the discrete model to incorporate an overall sparsity constraint and allowing to select individual elements from each group, leading to within-group sparsity. Furthermore, we discuss how this extension can be used to model hierarchical relationships between variables. We present a novel polynomial-time dynamic program that solves the hierarchical model selection problem exactly and discuss a tractable discrete relaxation.

We also interpret the implications of our results in the context of other group-based recovery frameworks. For instance, the convex approaches proposed in also relax the discrete constraint on the cardinality of the group support. However, they first need to decompose the approximation into vector atoms whose support consists only of one group and then penalize the norms of these atoms. It has been observed that these relaxations produce approximations that are group-sparse, but their group-support might include irrelevant groups. We concretely illustrate these cases via Pareto frontier examples on two different group structures.

Paper structure. The paper is organized as follows. In Section 2, we present definitions of group-sparsity and related concepts, while in Section III, we formally define the approximation and model-selection problems and connect them to the WMC problem. We present and analyze discrete relaxations of the WMC in Section IV and consider convex relaxations in Section V. In Section VI, we illustrate via a simple example the differences between the original problem and the relaxations. The generalized model is introduced and analyzed in Section VII, while numerical simulations are presented in Section VIII. We conclude the paper with some remarks in Section IX. The appendices contain the detailed descriptions of the dynamic programs.

II Basic Definitions

We start with the definition of totally unimodularity, a property of matrices that will turn out to be key for obtaining efficient relaxations of integer linear programs.

A totally unimodular matrix (TU matrix) is a matrix for which every square non-singular submatrix has determinant equal to $-1$ or $1$.

We now define the main building block of group sparse model selection, the group structure.

A group structure $\mathfrak{G}=\{\mathcal{G}_{1},\ldots,\mathcal{G}_{M}\}$ is a collection of index sets, named groups, with $\mathcal{G}_{j}\subseteq\mathcal{N}$ and $|\mathcal{G}_{j}|=g_{j}$ for $1\leq j\leq M$ and $\bigcup_{\mathcal{G}\in\mathfrak{G}}\mathcal{G}=\mathcal{N}$.

Another useful representation of a group structure is via an intersection graph $(\mathcal{V},\mathcal{E})$ where the nodes $\mathcal{V}$ are the groups $\mathcal{G}\in\mathfrak{G}$ and the edge set $\mathcal{E}$ contains $e_{ij}$ if $\mathcal{G}_{i}\cap\mathcal{G}_{j}\neq\emptyset$, that is an edge connects two groups that overlap. A sequence of connected nodes $v_{1},v_{2},\ldots,v_{n}$, is a cycle if $v_{1}=v_{n}$.

In order to illustrate these concepts, consider the group structure $\mathfrak{G}^{1}$ defined by the following groups, $\mathcal{G}_{1}=\{1\}$, $\mathcal{G}_{2}=\{2\}$, $\mathcal{G}_{3}=\{1,2,3,4,5\}$, $\mathcal{G}_{4}=\{4,6\}$, $\mathcal{G}_{5}=\{3,5,7\}$ and $\mathcal{G}_{6}=\{6,7,8\}$. $\mathfrak{G}^{1}$ can be represented by the variables-groups bipartite graph of Fig. 1 or by the intersection graph of Fig. 2, which is bipartite and contains cycles.

An important class of group structures is given by groups whose intersection graph is acyclic (i.e., a tree or a forest) and we call them acyclic group structures. A necessary, but not sufficient, condition for a group structure to have an acyclic intersection graph is that each element of the ground set occurs in at most two groups, i.e., the groups are at most pairwise overlapping. Note that a tree or a forest is a bipartite graph, where the two partitions contains the nodes that belong to alternate levels of the tree/forest. For example, consider $\mathcal{G}_{1}=\{1,2,3\}$, $\mathcal{G}_{2}=\{3,4,5\}$, $\mathcal{G}_{3}=\{5,6,7\}$, which can be represented by the intersection graph in Fig. 3(Left). If $\mathcal{G}_{3}$ were to include an element from $\mathcal{G}_{1}$, for example $\{2\}$, we would have the cyclic graph of Fig. 3(Right). Note that $\mathfrak{G}^{1}$ is pairwise overlapping, but not acyclic, since $\mathcal{G}_{3},\mathcal{G}_{4},\mathcal{G}_{5}$ and $\mathcal{G}_{6}$ form a cycle.

We anchor our analysis of the tractability of interpretability via selection of groups on covering arguments. Most of the definition we introduce here can be reformulated as variants of set covers on the support of a signal ${\mathbf{x}}$, however we believe it is more natural in this context to talk about group covers of a signal ${\mathbf{x}}$ directly.

The binary vector $\boldsymbol{\omega}$ indicates which groups are active and the constraint $\mathbf{A}^{\mathfrak{G}}\boldsymbol{\omega}\geq\iota({\bf x})$ makes sure that, for every non-zero component of ${\bf x}$, there is at least one active group that covers it. We also say that $\mathcal{S}({\bf x})$ covers ${\bf x}$. Note that the group cover is often not unique and $\mathcal{S}({\bf x})=\mathfrak{G}$ is a group cover for any signal ${\bf x}$. This observation leads us to consider more restrictive definitions of group covers.

A $G$-group cover $\mathcal{S}^{G}({\bf x})\subseteq\mathfrak{G}$ is a group cover for ${\bf x}$ with at most $G$ elements,

It is not guaranteed that a $G$-group cover always exists for any value of $G$. Finding the smallest $G$-group cover lead to the following definitions.

A signal ${\bf x}$ is $G$-group sparse with respect to a group structure $\mathfrak{G}$ if $\|{\bf x}\|_{\mathfrak{G},0}\leq G$.

In other words, a signal is $G$-group sparse if its support is contained in the union of at most $G$ groups from $\mathfrak{G}$.

III Tractability of interpretations

Although real signals may not be exactly group-sparse, it is possible to give a group-based interpretation by finding a group-sparse approximation and identifying the groups that constitute its support. In this section, we establish the hardness of finding group-based interpretations of signals in general and characterize a class of group structures that lead to tractable interpretations. In particular, we present a polynomial time algorithm that finds the correct $G$-group-support of the $G$-group-sparse approximation of ${\bf x}$, given a positive integer $G$ and the group structure $\mathfrak{G}$.

We first define the $G$-group sparse approximation $\hat{{\bf x}}$ and then show that it can be easily obtained from its $G$-group cover $\mathcal{S}^{G}(\hat{\bf x})$, which is the solution of the model selection problem. We then reformulate the model selection problem as the weighted maximum coverage problem. Finally, we present our main result, the polynomial time dynamic program for acyclic group structures.

If we already know the $G$-group cover of the approximation $\mathcal{S}^{G}(\hat{\bf x})$, we can obtain $\hat{\bf x}$ as $\hat{\bf x}_{\mathcal{I}}={\bf x}_{\mathcal{I}}$ and $\hat{\bf x}_{\mathcal{I}^{c}}=0$, where $\mathcal{I}=\bigcup_{\mathcal{G}\in\mathcal{S}^{G}(\hat{\bf x})}\mathcal{G}$ and $\mathcal{I}^{c}=\mathcal{N}\setminus\mathcal{I}$. Therefore, we can solve Problem 1 by solving the following discrete problem.

To show the connection between the two problems, we first reformulate Problem 1 as

The optimal solution is not changed if we introduce a constant, change sign of the objective and consider maximization instead of minimization

The internal maximization is achieved for $\hat{\bf x}$ as $\hat{\bf x}_{\mathcal{I}}={\bf x}_{\mathcal{I}}$ and $\hat{\bf x}_{\mathcal{I}^{c}}=0$, so that we have, as desired,

The following reformulation of Problem 2 as a binary problem allows us to characterize its tractability.

The proof follows along the same lines as the proof in . Note that in (4), $\boldsymbol{\omega}$ and $\bf y$ are binary variables that specify which groups and which variables are selected, respectively. The constraint $\mathbf{A}^{\mathfrak{G}}\boldsymbol{\omega}\geq{\bf y}$ makes sure that for every selected variable at least one group is selected to cover it, while the constraint $\sum_{j=1}^{M}\omega_{j}\leq G$ restricts choosing at most $G$ groups. ∎

Problem (4) can produce all the instances of the weighted maximum coverage problem (WMC), where the weights for each element are given by $x_{i}^{2}$ ($1\leq i\leq N$) and the index sets are given by the groups $\mathcal{G}_{j}\in\mathfrak{G}$ ($1\leq j\leq M$). Since WMC is in general NP-hard and given Lemma 1, the tractability of (3) directly depends on the hardness of (4), which leads to the following result.

The model selection problem (3) is in general NP-hard.

It is possible to approximate the solution of (4) using the greedy WMC algorithm . At each iteration, the algorithm selects the group that covers new variables with maximum combined weight until $G$ groups have been selected. However, we show next that for certain group structures we can find an exact solution.

Our main result is an algorithm for solving (4) for acyclic group structures. The proof is given in Appendix A.

Given an acyclic group structure $\mathfrak{G}$, there exists a polynomial time dynamic programming algorithm that solves (4).

Sets that are included in one another can be excluded because choosing the larger set would be a strictly dominant strategy, making the smaller set redundant. However, the correctness of the dynamic program is unaffected even if such sets are present, as long as the intersection graph remains acyclic.

It is also possible to consider the case where each group $\mathcal{G}_{i}$ has a cost $C_{i}$ and we are given a maximum group cost budget $C$. The problem then becomes the Budgeted Maximum Coverage . However, this problem is NP-hard, even in the non-overlapping case, because it generalizes the knapsack problem. However, similarly to the pseudo-polynomial time algorithm for knapsack , we can easily devise a pseudo-polynomial time algorithm for the weighted group sparse problem, even for acyclic overlaps. The only condition is that the costs must be integers. The time complexity of the resulting algorithm is then polynomial in $C$, the maximum group cost budget. The algorithm is almost the same as the one given in Appendix A: instead of keeping track of selecting $g$ groups, where $g$ varies from $1$ to $G$; we keep track of selecting groups with total weight equal to $c$, where $c$ varies from $1$ to $C$.

IV Discrete relaxations

Relaxations are useful techniques that allow to obtain approximate, or even sometimes exact, solutions while being computationally less demanding. In our case, we relax the constraint on the number of groups in (4) into a regularization term with parameter $\lambda>0$, which amounts to paying a penalty of $\lambda$ for each selected group. We then obtain the following binary linear program

In general, (6) is NP-hard, however, it is well known that if the constraint matrix $\mathbf{C}$ is Totally Unimodular (TU), then it can be solved in polynomial-time. While the concatenation of two TU matrices is not TU in general, the concatenation of the identity matrix with a TU matrix results in a TU matrix. Thus, due to its structure, $\mathbf{C}$ is TU if and only if $\mathbf{A}^{\mathfrak{G}}$ is TU [28, Proposition 2.1].

The next lemma characterizes which group structures lead to totally unimodular constraints.

Group structures whose intersection graph is bipartite lead to constraint matrices $\mathbf{A}^{\mathfrak{G}}$ that are TU.

We first use a result that establishes that if a matrix is TU, then its transpose is also TU [28, Proposition 2.1]. We then apply [28, Corollary 2.8] to $\mathbf{A}^{\mathfrak{G}}$, swapping the roles of rows and columns. Given a $\{0,1,-1\}$ matrix whose columns can be partitioned into two sets, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, and with no more than two nonzero elements in each row, this corollary provides two sufficient conditions for it being totally unimodular:

If two nonzero entries in a row have the same sign, then the column of one is in $\mathcal{S}_{1}$ and the other is in $\mathcal{S}_{2}$.

If two nonzero entries in a row have opposite signs, then their columns are both in $\mathcal{S}_{1}$ or both in $\mathcal{S}_{2}$.

In our case, the columns of $\mathbf{A}^{\mathfrak{G}}$, which represent groups, can be partitioned in two sets, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ because the intersection graph is bipartite. The two sets represents groups which have no common overlap so that each row of $\mathbf{A}^{\mathfrak{G}}$ contains at most two nonzero entries, one in each set. Furthermore, the entries in $\mathbf{A}^{\mathfrak{G}}$ are only or $1$, so that condition 1) is satisfied and condition 2) does not apply. ∎

Acyclic group structures lead to totally unimodular constraints.

Acyclic group structures have an intersection graph which is a tree or a forest, which is bipartite. ∎

The worst case complexity for solving the linear program (6), via a primal-dual method , is $\mathcal{O}(N^{2}(N+M)^{1.5})$, which is greater than the complexity of the dynamic program of Theorem 1. However, in practice, using an off-the-shelf LP solver may still be faster, because the empirical performance is usually much better than the worst case complexity.

Another way of solving the linear program for acyclic group structures is to reformulate it as an energy maximization problem over a tree, or forest. In particular, let $\psi_{i}=\|{\mathbf{x}}_{\mathcal{G}_{i}}\|_{2}^{2}$ be the energy captured by group $\mathcal{G}_{i}$ and $\psi_{ij}=\|{\mathbf{x}}_{\mathcal{G}_{i}\cap\mathcal{G}_{j}}\|_{2}^{2}$ the energy that is double counted if both $\mathcal{G}_{i}$ and $\mathcal{G}_{j}$ are selected, which then needs to be subtracted from the total energy. Problem (5) can then be formulated as

This problem is equivalent to finding the most probable state of the binary variables $\omega_{i}$, where their probabilities can be factored into node and edge potentials. These potentials can be computed in $\mathcal{O}(N)$ time via a single sweep over the elements, then the most probable state can be exactly estimated by the max-sum algorithm in only $\mathcal{O}(M)$ operations, by sending messages from the leaves to the root and then propagating other message from the root back to the leaves .

The next lemma establishes when the regularized solution coincides with the solution of (4).

If the value of the regularization parameter $\lambda$ is such that the solution $(\boldsymbol{\omega}^{\lambda},{\bf y}^{\lambda})$ of (5) satisfies $\sum_{j}\omega_{j}^{\lambda}=G$, then $(\boldsymbol{\omega}^{\lambda},{\bf y}^{\lambda})$ is also a solution for (4).

This lemma is a direct consequence of Prop. 3 below. ∎

However, as we numerically show in Section VIII, given a value of $G$ it is not always possible to find a value of $\lambda$ such that the solution of (5) is also a solution for (4). Let the set of points $\mathcal{P}=\{G,(f(G))\}_{G=1}^{M}$, where $f(G)=\sum_{i=1}^{N}y^{G}_{i}x_{i}^{2}$, be the Pareto frontier of (4). We then have the following characterization of the solutions of the discrete relaxation.

The discrete relaxation (5) yields only the solutions that lie on the intersection between the Pareto frontier of (4), $\mathcal{P},$ and the boundary of the convex hull of $\mathcal{P}$.

The scalarization of (7) yields the following discrete problem, with $\lambda>0$

whose solutions are the same as for (5). Therefore, the relationship between the solutions of (4) and (5) can be inferred by the relationship between the solutions of (7) and (8). It is known that the solutions of (8) are also Pareto optimal solutions of (7), but only the Pareto optimal solutions of (7) that admit a supporting hyperplane for the feasible objective values of (7) are also solutions of (8) [35, Section 4.7]. In other words, the solutions obtainable via scalarization belong to the intersection of the Pareto optimal solution set and the boundary of its convex hull. ∎

V Convex relaxations

One can in general use (9) to find a group-sparse approximation under the chosen group norm

VI Case study: discrete vs. convex interpretability

The following stylized example illustrates situations that can potentially be encountered in practice. In these cases, the group-support obtained by the convex relaxation will not coincide with the discrete definition of group-cover, while the dynamical programming algorithm of Theorem 1 is able to recover the correct group-cover.

Let $\mathcal{N}=\{1,\ldots,11\}$ and let $\mathfrak{G}=\{\mathcal{G}_{1}=\{1,\ldots,5\},~{}\mathcal{G}_{2}=\{4,\ldots,8\},~{}\mathcal{G}_{3}=\{7,\ldots,11\}\}$ be the acyclic group structure structure with $3$ groups of equal cardinality. Its intersection graph is represented in Fig. 4. Consider the $2$-group sparse signal ${\bf x}=[0~{}0~{}1~{}1~{}1~{}0~{}1~{}1~{}1~{}0~{}0]^{\top}$, with minimal group-cover $\mathcal{M}({\bf x})=\{\mathcal{G}_{1},\mathcal{G}_{3}\}$.

The dynamic program of Theorem 1, with group budget $G=2$, correctly identifies the groups $\mathcal{G}_{1}$ and $\mathcal{G}_{3}$. The TU linear program (5), with $0<\lambda\leq 2$, also yields the correct group-cover. Conversely, the decomposition obtained via (9) with unitary weights is unique, but is not group sparse. In fact, we have $\mathcal{S}({\bf x})=\breve{\mathcal{S}}({\bf x})=\mathfrak{G}$. We can only obtain the correct group-cover if we use the weights $[1~{}d~{}1]$ with $d>\frac{2}{\sqrt{3}}$, that is knowing beforehand that $\mathcal{G}_{2}$ is irrelevant.

This is an example where the correct minimal group-cover exists, but cannot be directly found by the Latent Group Lasso approach. There may also be cases where the minimal group-cover is not unique. We leave to future work, to investigate which of these minimal covers are obtained by the proposed dynamic program and characterize the behavior of relaxations.

VII Generalizations

In this section, we first present a generalization of the discrete approximation problem (4) by introducing an additional overall sparsity constraint. Secondly, we show how this generalization encompasses approximation with hierarchical constraints that can be solved exactly via dynamic programming. Finally, we show that the generalized problem can be relaxed into a linear binary problem and that hierarchical constraints lead to totally unimodular matrices for which there exists efficient polynomial time solvers.

In many applications, for example genome-wide association studies , it is desirable to find approximations that are not only group-sparse, but also sparse in the usual sense (see for an extension of the group lasso). To this end, we generalize our original problem (4) by introducing a sparsity constraint $K$ and allowing to individually select variables within a group. The generalized integer problem then becomes

The problem described above is a generalization of the well-known Weighted Maximum Coverage (WMC) problem. The latter does not have a constraint on the number of indices chosen, so we can simulate it by setting $K=N$. WMC is also well-known to be NP-hard, so that our present problem is also NP-hard, but it turns out that it can be solved in polynomial time for the same group structures that allow to solve (4).

Given an acyclic groups structure $\mathfrak{G}$, there exists a dynamic programming algorithm that solves (11) with complexity $\mathcal{O}(M^{2}GK^{2})$.

The dynamic program is described in Appendix A alongside the proof that it has a polynomial running time. ∎

VII-B Hierarchical constraints

The generalized model allows to deal with hierarchical structures, such as regular trees, frequently encountered in image processing (e.g. denoising using wavelet trees). In such cases, we often require to find $K$-sparse approximations such that the selected variables form a rooted connected subtree of the original tree, see Fig. 5. Given a tree $\mathcal{T}$, the rooted-connected approximation can be cast as the solution of the following discrete problem

where $\mathcal{T}_{K}$ denotes all rooted and connected subtrees of the given tree $\mathcal{T}$ with at most $K$ nodes.

This type of constraint can be represented by a group structure, where for each node in the tree we define a group consisting of that node and all its ancestors. When a group is selected, we require that all its elements are selected as well. We impose an overall sparsity constraint $K$, while discarding the group constraint $G$.

For this particular problem, for which relaxed and greedy approximations have been proposed , in Appendix B, we present a dynamic program that runs in polynomial time.

The time complexity of our dynamic program on a general tree is $\mathcal{O}(NK^{2}D)$, where $D$ is the maximum number of children that a node in the tree can have.

While preparing the final version of this manuscript, independently proposed a similar dynamic program for tree projections on $D$-regular trees with time complexity $\mathcal{O}(NKD)$. Following their approach, we improved the time complexity of our algorithm to $\mathcal{O}(NKD)$ for $D$-regular trees. We also prove that its memory complexity is $\mathcal{O}(N\log_{D}K)$. A computational comparison of the two methods, both implemented in Matlab, is provided in Section VIII, showing that our dynamic program can be up to $60\times$ faster, despite having similar worst-case time complexity.

The time complexity of our dynamic program on $D$-regular trees is $\mathcal{O}(NK^{2}D)$.

The space complexity of our dynamic program on $D$-regular trees is $\mathcal{O}(N\log_{D}K)$.

The description of the algorithm and the proof of its complexity, for both general and $D$-regular trees, can be found in Appendix B.

VII-C Additional discrete relaxations

By relaxing both the group budget and the sparsity budget in (11) into regularization terms, we obtain the following binary linear program

where $\mathbf{w}^{\top}=[x_{1}^{2},\ldots,x_{N}^{2},-\lambda_{G}\mathbf{1}_{M}^{\top},-\lambda_{K}\mathbf{1}_{N}^{\top}]$ and $\mathbf{C}=[\mathbf{I}_{N},~{}-\mathbf{A}^{\mathfrak{G}},~{}\mathbf{0}_{N}]$ and $\lambda_{G},\lambda_{K}>0$ are two regularization parameters that indirectly control the number of active groups and the number of selected elements. (13) can be solved in polynomial time if the constraint matrix $\mathbf{C}$ is totally unimodular. Due to its structure, by Proposition 2.1 in and that concatenating a matrix of zeros to a TU matrix preserves total unimodularity, $\mathbf{C}$ is totally unimodular if and only if $\mathbf{A}^{\mathfrak{G}}$ is totally unimodular. The next results proves that the constraint matrix of hierarchical group structures is totally unimodular.

Hierarchical group structures lead to totally unimodular constraints.

We use the fact that a binary matrix is totally unimodular if there exists a permutation of its columns such that in each row the $1$s appear consecutively, which is a combination of Corollary 2.10 and Proposition 2.1 in . For hierarchical group structures, such permutation is given by a depth-first ordering of the groups. In fact, a variable is included in the group that has it as the leaf and in all the groups that contain its descendants. Given a depth-first ordering of the groups, the groups that contain the descendants of a given node will be consecutive. ∎

The regularized hierarchical approximation problem, in particular

for $\lambda\geq 0$, has already been addressed by Donoho as the “complexity penalized residual sum-of-squares” and linked to the CART algorithm, which can found a solution in $\mathcal{O}(N)$ time. The condensing sort and select algorithm (CSSA) , with complexity $\mathcal{O}(N\log N)$, solves the problem where the indicator variable $\mathbf{y}$ is relaxed to be continuous in $$and constrains$\|\mathbf{y}\|_{1}$to be smaller than a given threshold$\gamma$, yielding rooted connected approximations that might have more than$K$ elements.

VIII Pareto Frontier Examples

The purpose of these numerical simulations is to illustrate the limitations of relaxations and of greedy approaches for correctly estimating the $G$-group cover of an approximation.

We consider the problem of finding a $G$-group sparse approximation of the wavelet coefficients of a given image, in our case a view of the Earth from space, see left inset in Fig. 6. We consider a group structure defined over the 2D wavelet tree. The wavelet coefficients of a 2D image can naturally be organized on three regular quad-trees, corresponding to a multi-scale analysis with wavelets oriented vertically, horizontally and diagonally respectively . We define groups consisting of a node and its four children, therefore each group has $5$ elements, apart from the topmost group that contains the scaled DC term and the first nodes of each of the three quad-trees. These groups overlap only pairwisely and their intersection graph is a tree itself, therefore leading to a totally unimodular constraint matrix. An example is given in the right inset in Fig. 7. For computational reasons, we resize the image to $16\times 16$ pixels and compute its Daubechies-4 wavelet coefficients. At this size, there are $64$ groups, but actually $52$ are sufficient to cover all the variables, since it is possible to discard the penultimate layer of groups while still covering the entire ground set.

Figures 6 and 7 show the Pareto frontier of the approximation error $\|{\bf x}-\hat{\bf x}\|_{2}^{2}$ with respect to the group sparsity $G$ for the proposed dynamic program. We also report the approximation error for the solutions obtained via the totally unimodular linear relaxation (TU-relax) (6) and the latent group lasso formulation (Latent GL) (10) with $p=2$, which we solved with the method proposed in . Fig. 7 shows the performance of StructOMP using the same group structure and of the greedy algorithm for solving the corresponding weighted maximum coverage problem.

We observe that there are points in the Pareto frontier of the dynamic program, for $G=5,10,30,31,50$, that are not achievable by the TU relaxation, since they do not belong to its convex hull. Furthermore, the latent group lasso approach often does not yield the optimal selection of groups, leading to a greater approximation error for the same number of active groups and it needs to select all $64$ groups in order to achieve zero approximation error. It is interesting to notice that the greedy algorithm outperforms StructOMP (see inset of Fig. 7), but still does not achieve the optimal solutions of the dynamic program. Furthermore, StructOMP needs to select all $64$ groups for obtaining zero approximation error, while the greedy algorithm can do with one less, namely $63$.

VIII-B Hierarchical constraints

We now consider the problem of finding a $K$-sparse approximation of a signal imposing hierarchical constraints. We generate a piecewise constant signal of length $N=64$, to which we apply the Haar wavelet transformation, yielding a $25$-sparse vector of coefficients ${\bf x}$ that satisfies hierarchical constraints on a binary tree of depth $5$, see Fig. 8(Left).

We compare the proposed dynamic program (DP) to the regularized totally unimodular linear program approach, two convex relaxations that use group-based norms and the StructOMP greedy approach . The first convex relaxation uses the Latent Group Lasso norm (9) with $p=2$ as a penalty and with groups defined as all parent-child pairs in the tree. We call this approach Parent-Child. This formulation will not enforce all hierarchical constraints to be satisfied, but will only ‘favor’ them. Therefore, we also report the number of hierarchical constraint violations. The second convex relaxation considers a hierarchy of groups where $\mathcal{G}_{j}$ contains node $j$ and all its descendants. Hierarchical constraints are enforced by the group lasso penalty $\Omega_{GL}({\bf x})=\sum_{\mathcal{G}\in\mathfrak{G}}\|{\bf x}_{\mathcal{G}}\|_{p}$, where ${\bf x}_{\mathcal{G}}$ is the restriction of ${\bf x}$ to $\mathcal{G}$, and we assess $p=2$ and $p=\infty$. We call this method Hierarchical Group Lasso. As shown in , solving $\min_{\mathbf{x}}\|\mathbf{y}-{\mathbf{x}}\|_{2}^{2}+\lambda\Omega_{GL}({\mathbf{x}})$, for $p=\infty$, is actually equivalent to solving the totally unimodular relaxation with the same regularization parameter. Once we determine the support of the solution, we assign to the components in the support the values of the corresponding components of the original signal. Finally, for the StructOMPWe used the code provided at http://ranger.uta.edu/~huang/R_StructuredSparsity.htm method, we define a block for each node in the tree. The block contains that node and all its ancestors up to the root. By finely varying the regularization parameters for these methods, we obtain solutions with different levels of sparsity.

In Figures 8(Right), we show the approximation error $\|{\bf x}-\hat{\bf x}\|_{2}^{2}$ as a function of the solution sparsity $K$ for the methods. The values of the DP solutions form the discrete Pareto frontier of the optimization problem controlled by the parameter $K$. Note that there are points in the Pareto frontier that do not lie on its convex hull, hence these solutions are not achievable by the TU linear relaxation. As expected, the Hierarchical Group LassoWe used the code provided at http://spams-devel.gforge.inria.fr/. with $p=\infty$ obtains the same solutions as the TU linear relaxation, while with $p=2$ it also misses the solutions for $K=21$ and $K=23$. The Parent-ChildWe used the algorithm proposed in . approach achieves more levels of sparsity (but still missing the solutions for $K=2,13$ and $15$), although at the price of violating some of the hierarchical constraints, i.e., we count one violation when one node is selected but not its parent. The StructOMP approach yields only few of the solutions on the Pareto frontier, but without violating any constraints. These observations lead us to conclude that, in some cases, relaxations of the original discrete problem or other greedy approaches might not be able to find the correct group-based interpretation of a signal.

In Fig. 9, we report a computational comparison between our dynamic program and the one independently proposed by Cartis and Thompson . We consider the problem of finding the $K=200$ sparse rooted connected tree approximation on a binary tree of a signal of length $2^{L}$, with $L=9,\ldots,18$, whose components are randomly and uniformly drawn from $$. Despite the two algorithms have the same computational complexity,$\mathcal{O}(NKD)$and are both implemented in Matlab, our dynamic program is between$20$and$60$ times faster.

IX Conclusions

Many applications benefit from group sparse representations. Unfortunately, our main result in this paper shows that finding a group-based interpretation of a signal is an integer optimization problem, which is in general NP-hard. To this end, we characterize group structures for which a dynamical programming algorithm can find a solution in polynomial time and also delineate discrete relaxations for special structures (i.e., totally unimodular constraints) that can obtain correct solutions.

Our examples and numerical simulations show the deficiencies of relaxations, both convex and discrete, and of greedy approaches. We observe that relaxations only recover group-covers that lie in the convex hull of the Pareto frontier determined by the solutions of the original integer problem for different values of the group budget $G$ (and sparsity budget $K$ for the generalized model). This, in turn, implies that convex and non-convex relaxations might miss some important groups or include spurious ones in the group-sparse model selection. We summarize our findings in Fig. 10.

There remain several interesting open questions which beg for answers. Firstly, there still lacks an intuitive understanding of under which circumstances the relaxations are able to yield the correct solutions. Secondly, our analysis implicitly assumes an orthogonal basis for the description of signals. In many machine learning and compressive sensing applications however, the structures in signals emerge only after representing them onto an overcomplete basis, e.g. shearlets or sparse coding techniques. Therefore, it would be interesting to explore to which extent our results can be generalized to the overcomplete setting.

Appendix A Dynamical programming for solving (11) for loopless pairwise overlapping groups

Here, we give the proof of Theorem 2. The proof of Theorem 1 follows along similar lines. We start by giving an intuitive understanding of the algorithm, followed by a formal description and proofs of correctness and complexity, both in time and space.

Problem (11) can be equivalently described by the following problem:

In this form, the problem described above is a generalization of the well-known Weighted Maximum Coverage (WMC) problem, which is NP-hard. In fact WMC is just a special case of SGSP with $K=N$. Although this makes it intractable in general, we show that this problem has some interesting structure. This structure allows us to build a dynamic program which can obtain the exact solution in polynomial time, for certain special classes of groups. We believe that this algorithm may be of independent interest outside the information theory community.

A-A An intuitive take on the Dynamic Programming Approach

We first present an informal account of the ideas behind our method. The basic idea we use is dynamic programming, i.e., we build the solution to the global optimization problem from solutions to subproblems. In our case, the subproblems correspond to looking at a subset of the group structure, and solving the optimization problem for this case. What we could do is to start from a single group and keep adding more groups one at a time, updating the optimal solution at each step. One may naïvely hope that such an approach would lead to the global solution. Unfortunately, this basic intuitive approach fails, as we illustrate next through some examples.

Example-1: Consider the case of $N=5$, with the weights being the vector $$. For the sake of illustration, let the group structure be$\mathfrak{G}=\{\mathcal{G}_{1},\mathcal{G}_{2}\}$, where$\mathcal{G}_{1}=\{1,2\},\mathcal{G}_{2}=\{2,3,4\}$. We wish to find the optimal solutions for the cases:

The optimal solutions can be found simply by observation.

The optimal solution for $K=2,G=1$, has weight $15$, and involves selecting group $\mathcal{G}_{2}$, and elements $\{2,4\}$.

The optimal solution for $K=3,G=2$, has weight $20$, and involves selecting both groups and elements $\{1,2,4\}$.

Example-2: Consider the case of $N=5$, with the weights being the vector $$. Let the set of groups be$\mathfrak{G}=\{\mathcal{G}_{1},\mathcal{G}_{2},\mathcal{G}_{3}\}$, where$\mathcal{G}_{1}=\{1,2\}$,$\mathcal{G}_{2}=\{2,3,4\}$,$\mathcal{G}_{3}=\{4,5\}$. We wish to find the optimal solution for the case:$K=4,G=2$. Once again, we can see that the optimal solution involves selecting groups$\mathcal{G}_{1}$and$\mathcal{G}_{3}$, and elements$\{1,2,4,5\}$, for a total value of$26$.

In the examples described above, the set of elements and their weights are the same. However, in the first example, any optimal solution for any meaningful values of the parameters G and K, involves $\mathcal{G}_{2}$. Yet, in the second example, we have a situation where the optimal selection does not involve $\mathcal{G}_{2}$, see Figure 11

A-A2 Boundary-cognizant DP

As we illustrated above, the simple DP approach does not work. When look at a subset of groups, some of which overlap with as yet unexplored groups, decisions regarding the overlapping groups are difficult to make. The reason is that the quality of a group in the view of the algorithm may decrease, if the high-weight elements in the group also happen to be contained in another overlapping group, which is seen in the future. While building partial solutions, we then need to consider both possibilities - an overlapping group is either included or excluded from the putative solution. We now introduce some notation which allows us to describe these ideas more concretely.

Our algorithm is heavily based on the intersection graph of the group structure. Thus, we will frequently refer to the groups as ‘nodes’ in our algorithm. Our approach involves exploring the nodes of the intersection graph one at a time and storing a list of optimal values from the explored nodes. These optimal values constitute the optimal weight of a $g$-group, $k$-element selection from the explored groups, for all $1\leq g\leq G$ and for all $1\leq k\leq K$. Further, we need to store these optimal values for each possible selection of the overlapping groups, so that we do not make decisions concerning such groups at the current step. In terms of the intersection graph, these overlapping groups are simply those nodes which belong to our currently explored set, but are also adjacent to some node which is not in the explored set. We call such nodes boundary nodes. Since our algorithm explores the intersection graph keeping track of all possibilities at the boundary nodes, it may fittingly be called a boundary-cognizant Dynamic Program.

Although this trick of being boundary-cognizant helps us get the correct solution, it can be expensive. Suppose we have $b$ boundary nodes at a certain step of the algorithm. Then the table of optimal values we seek to store has size $GK2^{b}$, which is exponential in $b$. For an arbitrary intersection graph, this factor can indeed be exponential; for example, a complete graph with $M$ nodes will always have $M-1$ boundary nodes at the penultimate step. However, if we restrict the intersection graph to be a tree, then it turns out there is a way to explore the graph such that the number of boundary nodes in a graph with $M$ nodes is only $O(\log M)$. This property allows our algorithm to run in polynomial time on such graphs.

A-B Optimal substructure

We expose the optimal substructure of this problem below by highlighting two key properties: Groups-elements dichotomy property and independence given the boundary property. These provide sufficient evidence that an optimal solution to our problem can be efficiently constructed from optimal solutions to subproblems, indicating the correctness of the dynamic programming approach. Further, we will use a slight generalization of property-2 in the proof of correctness of our algorithm.

Groups-elements dichotomy: Suppose we had access to an oracle who told us the set of $G$ groups that comprise the optimal solution to SGSP. Then we can easily recover the full solution using this information, by picking the $K$ largest-weight elements contained in the union of these $G$ groups.

Interestingly, the converse of the above is not true. If the oracle told us the list of $K$ elements contained in the optimal selection, but not the groups, the problem remains hard. Finding the $G$ groups that comprise the optimal solution is equivalent to finding a $G$-group cover for these $K$ elements, given that such a cover exists. If we could solve this task in polynomial time, the same algorithm would also solve the NP-Hard Set Cover problem in polynomial time.The idea described here is not a formal reduction. It is possible that the additional structure possessed by the optimal solution would allow us to recover the groups in polynomial time. However, there seems to be no clear way to use this additional structure, so the only obvious way to recover the groups is to solve a set-cover problem, which is NP-hard.

In a certain sense, the above shows that the difficult part of finding the optimal solution is selecting the groups. However, this does not imply that the element sparsity constraint is insignificant. It is easy to create problem instances where even a small change in $K$ significantly changes the optimal selection.

Independence given the boundary: Let $\mathfrak{G}$ be the complete set of groups, and let $\mathcal{S}\subset\mathfrak{G}$ be a subset of these groups. Let $\mathcal{B}(\mathcal{S})$ be the boundary nodes of $\mathcal{S}$, that is the nodes in $\mathcal{S}$ that are connected to nodes in its complement, $\mathcal{S}^{c}$. Once again, we assume the existence of an oracle who knows the true solution. Suppose this oracle tells us the following information:

The number of groups in $\mathcal{S}$ which are included in the optimal solution. Call this quantity $G_{1}$.

The number of elements in the optimal solution, which occur in any of the groups in $\mathcal{S}$. Call this quantity $K_{1}$.

The boundary nodes included in the optimal solution.

Then this information allows us to recover the optimal solution, by solving two independent optimization problems on the sets $\mathcal{S}$ and $\mathcal{S}^{c}$ respectively.

For ease of explanation, we refer to the set of boundary nodes included in the optimal selection as the set of ‘active boundary nodes’, $\mathcal{B}_{A}(\mathcal{S})$. Note that $\mathcal{B}_{A}(\mathcal{S})$ is known as it is given to us by the oracle. Further, we call the set of elements included in $\mathcal{B}_{A}(\mathcal{S})$ the set of ‘active boundary elements’, or $\mathcal{E}_{A}$.

Recovery Method: In order to recover the global optimal solution, we need to recover the selection of groups and elements in $\mathcal{S}$ and $\mathcal{S}^{c}$ respectively.

We first describe the procedure for $\mathcal{S}$. Consider all possible ways of choosing $K_{1}$ elements contained in $G_{1}$ groups from $\mathcal{S}$, such that the set of chosen groups in $\mathcal{B}(\mathcal{S})$ exactly matches $\mathcal{B}_{A}(S)$. Among these choices, the choice which has the maximum total weight of chosen elements gives us the selections of groups and elements in $\mathcal{S}$.

Proof: The proof of these two statements is straightforward. First, we formally show how to break the true optimal solution into two disjoint components. After this, we argue that the two components constitute optimal solutions to smaller optimization problems.

Let us denote the set of groups and elements in the global optimal solution by $\mathfrak{G}^{*}$ and ${\mathcal{E}}^{*}$, respectively. We create two new group-element selections, roughly corresponding to $\mathcal{S}$ and $\mathcal{S}^{c}$, which we shall denote by $({\mathfrak{G}}_{1},\mathcal{E}_{1})$ and $({\mathfrak{G}}_{2},\mathcal{E}_{2})$ respectively. These two components are constructed as follows:

The set of selected groups in $\mathcal{S}$ and $\mathcal{S}^{c}$ are already disjoint, so these are directly assigned to ${\mathfrak{G}_{1}}$ and ${\mathfrak{G}_{2}}$ respectively.

For any element in ${\cal E}$ which occurs only in (groups in) $\mathcal{S}$, assign it to ${\cal E}_{1}$.

For any element in ${\cal E}$ which occurs only in $\mathcal{S}^{c}$, assign it to ${\cal E}_{2}$.

For any element which occurs in $\mathcal{S}$ as well as $\mathcal{S}^{c}$ (and hence in $\mathcal{B}(\mathcal{S})$), first try to assign it to ${\cal E}_{1}$. That is, check if this element is contained in $\mathfrak{G}_{1}$, and if so assign the element to ${\cal E}_{1}$. If not, we assign it to ${\cal E}_{2}$.

$\mathfrak{G}_{1}$ and $\mathfrak{G}_{2}$ form a partition of $\mathfrak{G}^{*}$, and similarly $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ form a partition of $\mathcal{E}^{*}$.

$({\mathfrak{G}}_{1},\mathcal{E}_{1})$, $({\mathfrak{G}}_{2},\mathcal{E}_{2})$ represent valid group-element selections over the sets of groups $\mathcal{S}$ and $\mathcal{S}^{c}$ respectively (i.e. $\mathcal{E}_{1}$ is contained in the union of groups in $\mathfrak{G}_{1}$, and similarly $\mathcal{E}_{2}$ is contained in $\mathfrak{G}_{2}$.)

$|\mathfrak{G}_{1}|=G_{1}$, $|\mathcal{E}_{1}|=K_{1}$, $|\mathfrak{G}_{2}|=G_{2}$, $|\mathcal{E}_{2}|=K_{2}$, where $G_{1}$, $K_{1}$, $G_{2}$, $K_{2}$ are defined as above.

We are now ready to prove the correctness of the recovery method.

An identical argument shows that ($\mathfrak{G}_{1},\mathcal{E}_{1}$) represents an optimal $G_{1}$-group, $K_{1}$-element selection over $\mathcal{S}$, among all group-element selections for which the set of chosen nodes from $\mathcal{B}(\mathcal{S})$ equals exactly $\mathcal{B}_{A}(\mathcal{S})$. This proves the correctness of our recovery method. ∎

A-C Overview of our Algorithm

Our algorithm explores the acyclic intersection graph one node at a time, storing the optimal solution among the visited nodes and eventually leading to the optimal solution for the entire graph. It is described by two rules: the Value Update Rule and the Graph Exploration Rule.

Graph Exploration Rule: This rule takes as input a given tree graph, and outputs an order of exploring the graph so as to minimize the number of encountered boundary nodes.

Value Update Rule: The Value Update Rule determines how to update the list of optimal values when we explore a new node.

We first describe the Value Update Rule. While doing so, we assume that the nodes of the graph have been labelled $1,2,\dots,M$ in some suitable manner, and explore them in this order. In order to lay the foundation for describing the update rule, we will first define the table of optimal values maintained by our algorithm, and ensure that the given data is in suitable format.

A-D Table of optimal values

We describe the set of optimal solutions stored by our Table of optimal values. Abstractly, this table can be thought of as a mathematical function with $5$ different parameters. These are described below:

Explored Set : $\mathcal{S}\subseteq\mathfrak{G}$ This is any subset of nodes of the intersection graph. It represents the set of nodes currently visited by our algorithm.

Group Count : $g\in\{1,2,\dots,G\}$. This is the maximum number of groups we are allowed to select.

Element Count : $k\in\{1,2,\dots,K\}$. This is the maximum number of elements we are allowed to select.

Boundary Indicator Vector : $\mathbf{I}_{\mathbf{b}}\in\{0,1\}^{B}$. This is a binary vector of size $B$. Given a boundary set vector, $\mathbf{b}$, for each $i\in\{1,\dots,B\}$, the $i$-th component of $\mathbf{I}_{\mathbf{b}}$ is either or $1$, representing whether the group $\mathbf{b}_{i}$ is selected or excluded in the optimal selection. We also allow $\mathbf{I}_{\mathbf{b}}$ to be an empty vector.

We now define our optimal values function as follows.

$\bf F(\mathcal{S},g,k,\mathbf{b},\mathbf{I}_{\mathbf{b}})$ represents the maximum weight obtainable by selecting at most $k$ elements contained in a union of at most $g$ groups from the set $\mathcal{S}$, with the choice of selections among the set of boundary nodes $\mathbf{b}$ given by $\mathbf{I}_{\mathbf{b}}$. This function is defined for the entire range of its arguments mentioned above.

Although the function is defined for all $\mathcal{S}\subseteq\mathfrak{G}$, in practice we explore the nodes one at a time, in serial order. Thus, we only need to keep track of $M$ different sets of explored nodes, where the $i$-th set, $\mathcal{S}_{i}$, consists of groups $\mathcal{G}_{1},\mathcal{G}_{2},\dots,\mathcal{G}_{i}$, for all $i\in\{1,\dots,M\}$. Furthermore, we only see $M$ different sets of boundary nodes for a given intersection graph, $\mathcal{B}(\mathcal{S}_{i})$. In certain intermediate steps we shall find it convenient to use in place of $\mathcal{B}(\mathcal{S}_{i})$, a different set than the actual set of boundary nodes.

A-E Data Format and Notation

Without loss of generality, we can assume that each group has no more than $K$ elements. Further, we will assume that the indices in each group are specified in decreasing order of weights.

In case the above assumptions are not met a-priori, we can do some preprocessing on the given data. Since we know that each group consists of at most $N$ elements, we can pick the largest $K$ elements and then sort them in $O(N+K\log N)$ time. This can be done by building a max-heap of all $N$ elements and then extracting the topmost element $K$ times. Since we need to do this for each one of the $M$ groups, this leads to a total complexity of $O(MN+MK\log N)$.

While describing the complexity of our main algorithm, we will assume that the groups are already represented in the above canonical form. Hence, we will not consider the above term in our expression for time complexity.

Next, we formally define some notation that we use in our description of the value update rule.

Concatenation Operator: Given two vectors ${\mathbf{x}}$ and $\mathbf{y}$ of lengths $m$ and $n$ respectively, we define the vector ‘$x$ concatenated with $y$’, written as ${\mathbf{x}}.\mathbf{y}$, to be an $m+n$-length vector which consists of entries of ${\mathbf{x}}$ followed by entries of $\mathbf{y}$.

Best-k operator: We define a function $H(\mathcal{S},k)$ to represent the optimal value for choosing $k$ elements from a set $\mathcal{S}$. The set $\mathcal{S}$ could be a single group, a union of groups, or any well-defined collection of elements. As noted earlier, $H(\mathcal{S},k)$ simply equals the sum of the $k$ largest weight elements in $\mathcal{S}$.

A-F Value Update Rule

We shall now describe the Value Update Rule. This rule shows us how to find the optimal solution to SGSP, which is represented by the value: $F(\mathfrak{G},G,K,\emptyset,\emptyset)$.

Base Case. We start with $\mathcal{S}_{0}={\emptyset}$. For this case, all values of $F$ are set to : $F(\emptyset,g,k,\emptyset,\emptyset)=0~{}\forall g,k$.

Update. The update case describes how to recompute the list of optimal values when we explore a new node. We shall apply this rule a total of $M$ times, exploring one new node from the graph each time, and updating our table of values. At the end, we can simply read off the solution from the appropriate entry of the table.

Since we explore the nodes in serial order, at the $i$-th step, our explored set will consist of nodes $1,2,\dots,i$. As mentioned earlier, we denote our explored set after the $i$-th step as $\mathcal{S}_{i}$, and the boundary set vector at this time as $\mathbf{b}_{i}$. We use the notation $\mathcal{G}_{j}$ to refer to the $j$-th group, which is also the $j$-th node of the intersection graph as per our chosen ordering. At the end of the $i$-th step, we assume that we have stored the values of $F$ for the explored set $\mathcal{S}_{i}$ and boundary set vector $\mathbf{b}_{i}$ for each possible value of parameters $g$, $k$, and the indicator variable $\mathbf{I}_{\mathbf{b}_{i}}$, in their respective ranges. Thus, the following values are available to us:

Our objective is to extend these values to the case when we have explored the $i+1$-th node. In other words, defining $\mathcal{S}_{i+1}\triangleq\mathcal{S}_{i}\cup\{\mathcal{G}_{i+1}\}$, we wish to obtain the following set of values:

where $\mathbf{b}_{i+1}$ represents the boundary nodes at time $i+1$ in vector form.

We now describe our method for obtaining these values. When we first consider node $i+1$, we treat it as a new boundary node and compute the optimal values for it being included or excluded from the putative solution. After this, we test for boundary nodes that have fallen into the interior of the explored set. For these redundant boundary nodes, we no longer need to store two separate values for the node being included or excluded, so we condense these into a single value. Our update rule thus consists of $3$ steps:

In this case, we are computing the optimal value for selecting $k$ elements contained in a union of $g$ groups among the first $(i+1)$ groups when the $(i+1)$-th group is not selected, and the groups in $\mathcal{B}_{i}$ are selected as per the indicator variables. Since the $(i+1)$-th group is not chosen, all our groups and elements must be chosen from among the first $i$ groups, with the same restrictions on the choice of boundary nodes. Hence, all optimal values for this case are equal to the corresponding values for $\mathcal{S}_{i}$.

for all $1\leq g\leq G$ and $1\leq k\leq K$ and all $\mathbf{I}_{\mathbf{b}_{i}}\in\{0,1\}^{B_{i}}$.

Case (a): The new node is included and does not overlap with any explored node.

for all $1\leq g\leq G$ and $1\leq k\leq K$ and all $\mathbf{I}_{\mathbf{b}_{i}}\in\{0,1\}^{B_{i}}$.

Case (b): The new node is included but overlaps with some explored nodes.

The update rule is the same as for case (a), but the elements in the region of overlap between the new node and the selected explored nodes must not be considered as being part of the new node. For this step, we need to know exactly which nodes have been chosen while computing an optimal value. This is the reason why we need to store separate values for each boundary node.

for all $1\leq g\leq G$ and $1\leq k\leq K$ and all $I_{B_{i}}\in\{0,1\}^{B_{i}}$, where

That is we “clean” $\mathcal{G}_{i+1}$ of the overlap with the currently selected boundary nodes.

After performing the above steps, the number of stored values will be doubled. We can reduce them: for each boundary node which has fallen into the interior of the explored nodes, we combine the optimal values for it being selected or excluded, into a single value by taking the larger of the two values. Each such operation reduces the number of stored values by half and we perform it after each value update. Unlike the earlier steps, this step may have to be performed multiple times in a single update.

The correctness of our algorithm relies on the correctness of the value update rule. Below, we argue for the correctness of this rule for each of its $3$ steps.

Step 1: The correctness of this step is self-evident.

Step 2, case (a): Since this step is a special case of step 2, case (b), it is sufficient to prove correctness of the latter.

Step 2, case (b): We prove the correctness of this step using the optimal substructure property 2 described in section A-B.

Our task is to find the optimal selection of $g$-groups and $k$-elements from the set $\mathcal{S}_{i+1}\triangleq\mathcal{S}_{i}\cup\mathcal{G}_{i+1}$, when $\mathcal{G}_{i+1}$ is selected, and nodes in $\mathbf{b}_{i}$ are selected according to $\mathbf{I}_{\mathbf{b}_{i}}$. We now consider only the graph consisting of nodes in $\mathcal{S}_{i+1}$. With reference to the substructure property, choose the set $\mathcal{S}$ to be equal to $\mathcal{S}_{i}$. Critically, note that all groups in $\mathcal{B}(\mathcal{S})$ are contained in $\mathbf{b}_{i}$, and thus we store optimal values separately for these.

Although the substructure property 2 was derived on a graph with no additional information, it is equally well-applicable when certain groups (such as $\mathcal{G}_{i+1}$, and groups in $\mathbf{b}_{i}$) are constrained to be selected or excluded in the optimal solution. This property had three preconditions, one of which was the knowledge of boundary nodes in the optimal solution. This is trivially true, since in this particular optimization problem, the selection of groups in $\mathbf{b}_{i}$ is already fixed by $\mathbf{I}_{\mathbf{b}_{i}}$. Then, the property shows us that if we also know the number of groups and elements chosen from the two parts of the graph, we can recover the optimal solution over $\mathcal{S}_{i+1}$ by solving two separate optimization problems over $\mathcal{S}_{i}$ and $\mathcal{G}_{i+1}$ respectively.

Step 3: This is the condensation step. The correctness of this step follows from the interpretation of the objective function - $F(\mathcal{S},g,k,\mathbf{b},\mathbf{I}_{\mathbf{b}})$ represents the optimal values for $g$-group $k$-element selections, when the choices of groups in $\mathbf{b}$ are fixed by $\mathbf{I}_{\mathbf{b}}$. Thus, for groups that are not in $\mathbf{b}$, we need to consider both whether the node is included or excluded. Therefore, in order to remove a node from the set $\mathbf{b}$, we simply take the maximum value of the two cases.

Running Time

The running time of our algorithm is determined by 2 steps - Value Update rule and the Graph Exploration algorithm. As we explain later, the exploration rule can be implemented independently and is computationally much faster, so the time complexity is determined by the value update rule. We analyze the complexity of each step of the update rule below.

Complexity of step 1: All optimal values for this case are simply the optimal values computed before the node is explored. Thus, the update in this case corresponds simply to a table-copying operation. In fact, this copying can be avoided entirely by some clever bookkeeping; all we need to do is remember where the appropriate values are stored in memory. Thus, this step is very inexpensive from a computational point of view.

Complexity of step 3: Since condensation removes an explored node from the boundary set forever, it will have to be performed at most $M$ times in the entire algorithm. Since the set of boundary nodes at each step is fully determined by the intersection graph and the exploration ordering, these can be precomputed without significant time cost. Hence, we assume these are available to us and ignore their complexity. Then the complexity of a single condensation step is determined only by the number of values that need to be condensed, and is given by $\mathcal{O}(GK2^{B_{i}^{\prime}})$, which also equals $\mathcal{O}(GK2^{B_{i}})$.

Overall time complexity: Among the above, the most expensive case is step 2, case (b). The complexity of this step as obtained earlier equals $\mathcal{O}(GK^{2}2^{B_{i}}+KB_{i}2^{B_{i}}\log K)$, for the $i+1$-th value update. We need to perform this step $M$ times, with the parameter $i$ varying from to $M-1$ in the above expression.

Let $B^{*}$ be the maximum number of boundary nodes encountered by the algorithm at any step, i.e., $B^{*}=\max_{i}B_{i}$. Then the running time of our update algorithm is bounded by $\mathcal{O}(M(2^{B^{*}}K^{2}G+2^{B^{*}}{B^{*}}K\log K))$. Our graph exploration rule allows us to explore the graph so that $B^{*}$ is logarithmic in $M$, specifically $B^{*}\leq(\log_{2}M+1)$. Hence $2^{B^{*}}=\mathcal{O}(M)$. Using this in our above expression, we see that the complexity becomes $\mathcal{O}(M^{2}K^{2}G+M^{2}K\log M\log K)$ which shows that our algorithm is polynomial time. If we ignore logarithmic terms, we can write the complexity more compactly as $\mathcal{O}(M^{2}K^{2}G)$.

Space Complexity and Backtracking

We now look at the amount of space (memory) required by our algorithm. To account for this, we also need to describe how we will backtrack, i.e., how we find the optimal selection of groups and elements. Note that the method described above yields the optimal value for selecting $K$ elements from $G$ groups, but does not immediately tell us which groups are selected. We chose a backtracking method which is time-efficient, but involves storing a fair amount of data. Specifically, we store the optimal values obtained at each step of the value update rule prior to condensation, i.e., $F(\mathcal{S}_{i},g,k,\mathbf{b}_{i-1}.\{\mathcal{G}_{i}\},\mathbf{I}_{\mathbf{b}_{i-1}}.\{0\})$ and $F(\mathcal{S}_{i},g,k,\mathbf{b}_{i-1}.\{\mathcal{G}_{i}\},\mathbf{I}_{\mathbf{b}_{i-1}}.\{1\})$ for all $1\leq g\leq G$ , $1\leq k\leq K$, $\mathbf{I}_{\mathbf{b}_{i-1}}\in\{0,1\}^{B_{i-1}}$, $i\in\{1,\dots,M\}$. Thus the number of values we shall need to store is at most $MGK2^{B^{*}}$, which can be simplified to $\mathcal{O}(M^{2}KG)$ using $2^{B^{*}}=\mathcal{O}(M)$ (due to our graph exploration algorithm).

Our algorithm for backtracking is as follows: We start from the $M$-th node and work backwards, determining the number of elements selected from each group. For the $M$-th group, we look at the optimal value for $G$ groups and $K$ elements, for the 2 cases when $\mathcal{G}_{M}$ is selected or unselected. The value which is the larger of these two forms our optimal solution, and thus tells us whether or not $\mathcal{G}_{M}$ is chosen in the optimal selection. If the optimal values stored at the $M-1$-th step involve other boundary nodes besides node $M$, we maximize over all selections of these boundary nodes, since we don’t care about any particular nodes being selected in the optimal solution. We also remember the assignment of the indicator variables which allows us to obtain the largest value of $F$, since it tells us which nodes in $\mathbf{b}_{M-1}$ are included in the optimal solution. If we find that $\mathcal{G}_{M}$ is not chosen in the optimal selection, then we can ignore that group and simply find the optimal $G$-group, $K$-element selection on $M-1$ groups.

It can be verified that the running time of the above algorithm is somewhat smaller than the update rule. Thus, the overall expression for time complexity is unchanged even when we account for backtracking.

A-G Graph Exploration Rule

We determine the order with which the nodes are picked by a value associated to each subtree of the graph, which we call the $D$-value. In the following, we describe how it is computed, how it depends logarithmically on the number of nodes in the graph and how the number of boundary nodes is bounded by the $D$-value.

The D-value of a rooted tree graph is a non-negative integer associated with the graph. We will define the D-value algorithmically later.

Computing $D$-values: The procedure for computing the $D$-values is also recursive. If the tree has only one node, $D=1$. Now, assume the $R$ subtrees at a node $Q$ have values $D_{1}\geq\ldots\geq D_{R}$. Then, $D(Q)=\max(D_{1},D_{2}+1)$. In case there is no second subtree, $D(Q)=D_{1}$. We then have the following bound on the $D$-values.

The $D$-value of a rooted tree graph is logarithmic in the number of nodes, i.e. $D(G)\leq\log_{2}(M)+1$.

Let $D$ be a positive integer and $N(D)$ be the minimum number of nodes that a rooted tree must have in order to have $D$-value of D. We prove by induction that

Base case: $D=2$. A tree with only one node will have a $D$-value of 1. So to have a $D$-value of 2, we require a graph with at least 2 nodes. Hence (15) is satisfied.

Inductive case: $D>2$. Let $\mathcal{T}$ be a smallest (i.e. minimum node) rooted tree graph whose $D$-value is equal to $D$. Spread out $\mathcal{T}$ in the form of root and subtrees. Let the subtrees be $\mathcal{T}_{1},\mathcal{T}_{2},\ldots,\mathcal{T}_{k}$, with corresponding $D$-values $D_{1},D_{2},\ldots,D_{k}$. Without loss of generality, assume that $D_{1}\geq D_{2}\geq\ldots\geq D_{k}$. By definition, $D(\mathcal{T})=\max(D_{1},D_{2}+1)$.

By our assumption, $\mathcal{T}$ is a minimum-node graph with $D$-value equal to $D$, hence we cannot have $D_{1}=D(G)=D$, since that would give us a smaller rooted tree graph ($\mathcal{T}_{1}$) with a $D$-value of $D$. This means that $D_{1}<D$, and since $D=max(D_{1},D_{2}+1)$, hence $D_{2}+1=D$, i.e. $D_{2}=D-1$. Since $D_{1}\geq D_{2}=D-1$ and $D_{1}<D$, then $D_{1}=D-1=D_{2}$. Thus, the graph $\mathcal{T}$ has 2 subtrees ($\mathcal{T}_{1}$ and $\mathcal{T}_{2}$), with $D$-values of $D-1$ each. By definition, any rooted subtree with a $D$-value of $D-1$ must have at least $N(D-1)$ nodes. By our induction hypothesis, $N(D-1)\geq 2^{D-2}$ . Therefore, $\mathcal{T}$ has at least $2\times 2^{D-2}=2^{D-1}$ nodes. But since $\mathcal{T}$ was the smallest rooted tree graph with $D$-value of $D$, this means that $N(D)\geq 2^{D-1}$, as required. ∎

We now link the number of boundary nodes visited by the algorithm to the $D$-value of the intersection graph.

The total number of boundary nodes encountered by the graph exploration algorithm cannot exceed the $D$-value of the graph.

Let $\mathcal{T}$ be the given rooted tree graph, with $M$ nodes. We shall consider the number of boundary nodes when there is a ghost node connected to the root node. The ghost node is a hypothetical node which is not really a part of the graph, but still makes adjacent explored nodes count as boundary nodes. The ghost node captures the fact when we are running the algorithm recursively on a subtree, there will be an additional (potentially unexplored) node connected to the root of the subtree, which may lead to the root being counted as a boundary node. Let $B^{*}(\mathcal{T})$ denote the maximum number of boundary nodes encountered on $\mathcal{T}$ when we pick nodes according to our algorithm, and let $B^{*}_{G}(\mathcal{T})$ represent the same when we also have the ghost node. Clearly, $B^{*}_{G}(\mathcal{T})\geq B^{*}(\mathcal{T})$, hence it is enough to prove the following:

We prove this by strong induction on $M$.

Base Case. Suppose the rooted tree graph $\mathcal{T}$ has only $1$ node. Then the maximum number of boundary nodes encountered is obviously $1$, which is equal to the $D$-value of the graph (by definition). Hence $B^{*}_{G}(\mathcal{T})\leq D(\mathcal{T})$.

Inductive Case. When the graph $\mathcal{T}$ consists of $M$ nodes, $M>1$, consider the graph to be spread out in the form of root and subtrees. Compute the $D$-values for each rooted subtree, where w.l.o.g., $D_{1}\geq D_{2}\geq\ldots D_{k}$. Let $\mathcal{T}_{1},\mathcal{T}_{2},\ldots,\mathcal{T}_{k}$ be the corresponding subtrees. By definition, our algorithm explores nodes in the sequence: $\mathcal{T}_{1},\text{root},\mathcal{T}_{2},\mathcal{T}_{3},\ldots\mathcal{T}_{k}$.

Since each subtree has strictly fewer than $M$ nodes, each subtree satisfies (16) by the induction hypothesis. Also, notice that when exploring the subtree $\mathcal{T}_{1}$ of $\mathcal{T}$, the number of boundary nodes encountered is less than or equal to the number of boundary nodes encountered when exploring $\mathcal{T}_{1}$ as a standalone rooted-tree-graph, with a ghost node connected to its root. By definition, this is exactly equal to $B^{*}_{G}(\mathcal{T}_{1})$, which by our induction hypothesis is bounded by $D_{1}$. Therefore, the number of boundary nodes encountered while exploring $\mathcal{T}_{1}$ in $\mathcal{T}$ cannot exceed $D_{1}$. Once we are finished with $\mathcal{T}_{1}$, we pick the root, so the total number of boundary nodes is $1$. We now proceed to pick $\mathcal{T}_{2}$. By a similar argument, the maximum number of boundary nodes in $\mathcal{T}_{2}$ at any point cannot exceed the number of boundary nodes encountered while exploring $\mathcal{T}_{2}$ as a standalone graph with attached ghost node. In addition, the root of $\mathcal{T}$ can contribute at most 1 additional boundary node (In fact, the ghost node for $\mathcal{T}$ ensures that the root, once picked, will always contribute an additional boundary node). Therefore, the total number of boundary nodes in $\mathcal{T}$ while exploring $\mathcal{T}_{2}$ is at most $D_{2}+1$. Similar arguments hold for all other subtrees — the maximum number of boundary nodes while exploring the $k$-th subtree will be at most $D_{k}+1$, which is upper bounded by $D_{2}+1$.

Therefore, the maximum number of boundary nodes encountered at any step while exploring $\mathcal{T}$ is $B_{G}^{*}(\mathcal{T})\leq\max(D_{1},D_{2}+1)$. By definition, $D(\mathcal{T})=\max(D_{1},D_{2}+1)$. Therefore $B_{G}^{*}(\mathcal{T})\leq D(\mathcal{T})$. ∎

Combining Lemmas 3 and 4, we have the following result.

The maximum number of boundary nodes at any step of the algorithm is logarithmic in the number of nodes, i.e., $B\leq\log_{2}(M)+1$.

The previous lemma establishes the polynomial time complexity of the dynamic program for solving the generalized integer problem (11).

We shall now prove that the exploration rule itself requires minimal computation. This will justify our earlier claim that the running time is determined solely by the value update rule.

The running time of the graph exploration rule is $\mathcal{O}(M)$ for an $M$-node graph.

The exploration rule can be algorithmically run in two loops. In the first, we compute all $D$-values of all the required subtrees in the graph. In the second loop, we find the exploration ordering using these $D$-values. Note that the subtrees encountered by our recursive $D$-value computing algorithm are exactly the same set of subtrees encountered by our exploration rule, which makes it possible to compute all the required $D$-values in a single loop.

For computing $D$-values at a particular node, we use the formula $D=\max(D_{1},D_{2}+1)$, where $D_{1}\geq D_{2}\geq D_{3},D_{4},\dots,D_{k}$. Thus, we need to find the largest and second largest $D$-values among the subtrees. For a node with $d$ children, this takes $\mathcal{O}(d)$ time. Since the values $D_{1},D_{2},\dots,D_{k}$ are obtained recursively, this is the only computation which needs to be performed at the current node. Hence, the total time required is proportional to $\sum_{v\in\mathcal{V}}\max(d(v),1)\leq 2M$, where $d(v)$ represents the number of children that node $v$ has. Hence, this loop runs in $\mathcal{O}(M)$ time.

Obtaining the exploration order is similar. We only need to find the subtree with the largest $D$-value at the current node, so that we can pick the subtrees in the right order. This takes $\mathcal{O}(d)$ time for a node with $d$ children, and hence $\mathcal{O}(M)$ time for all nodes. Since both the above steps are $\mathcal{O}(M)$, the graph exploration rule itself runs in $\mathcal{O}(M)$, i.e., linear time. ∎

The proposed dynamic program solves the Weighted Maximum Coverage problem with an additional constraint on element sparsity for acyclic group structures. Its time complexity is $\mathcal{O}(M^{2}GK^{2})$, where $M$ is the number of groups, $G$ is the group sparsity budget and $K$ is the element sparsity budget.

Appendix B Dynamical programming for solving the hierarchical signal approximation problem (12)

Here we describe the dynamic program for solving the hierarchical signal approximation problem (12) and show that its time complexity is $\mathcal{O}(NK^{2}D)$, for general trees with maximum degree $D$ and $\mathcal{O}(NKD)$ for $D$-regular trees. Furthermore, its space complexity for $D$-regular trees is $\mathcal{O}(N\log_{D}K)$.

Problem (12) can be equivalently rephrased as the following optimization problem.

Rooted-Connected Subtree Problem: Given a rooted tree $\mathcal{T}$ with each node having at most $D$ children, a non-negative real number (weight) assigned to every node and a positive integer $K$, choose a subset of its nodes forming a rooted-connected subtree that maximizes the sum of weights of the chosen elements, such that the number of selected nodes does not exceed $K$. In our case, (12), the weight of a node is the square of the value of the component of the signal associated to that node. The proposed algorithm leverages the optimal substructure of the problem.

B-B Optimal substructure

Suppose that a particular node X belongs to the optimal $K$-node rooted-connected subtree. Consider the subtree $\mathcal{T}_{X,d}$ obtained by choosing X, $d$ of its children ($1\leq d\leq D$) and all descendants of these children. Consider the set of nodes $\mathcal{S}$ consisting of all the nodes of $\mathcal{T}_{X,d}$ which are also present in the optimal $K$-node rooted-connected subtree. Suppose there are $L$ nodes in $\mathcal{S}$. Then the nodes in $\mathcal{S}$ form the optimal $L$-node rooted-connected subtree at X, for the subgraph $\mathcal{T}_{X,d}$. See Fig. 13 for an example.

B-C Dynamic Programming method.

For every node X, we store the weight of the optimal $k$-node rooted-connected subtree at X, using only the nodes in the $d$ rightmost children of X and their descendants, for each $k$ and $d$ such that $1\leq k\leq K$ and $1\leq d\leq D$. We define a function $F(X,k,d)$, to store these optimal values. We start from the leaf nodes and move upwards, for each node assessing all its subtrees from right to left, eventually covering the entire tree. At the end, the optimal value will be given by $F(\text{root},K,D)$, that is the value of the best K-node rooted connected subtree of the root considering all its descendants.

Base Case. For every leaf node X and for all $1\leq k\leq K$ and $1\leq d\leq D$, we set $F(X,k,d)=\text{Weight}(X)$.

Inductive Case. By induction, for every non-leaf node X, all the F-values are known for the descendants of X. Let $X_{1},X_{2},\ldots X_{d}$ be the $d$ children of X in the right-to-left order, where $1\leq d\leq D$. Then, we compute the F-values of X using the following update rules.