A Nearly-Linear Time Framework for Graph-Structured Sparsity

1 Introduction

Over the past decade, sparsity has emerged as an important tool in several fields including signal processing, statistics, and machine learning. In compressive sensing, sparsity reduces the sample complexity of measuring a signal, and statistics utilizes sparsity for high-dimensional inference tasks. In many settings, sparsity is a useful ingredient because it enables us to model structure in high-dimensional data while still remaining a mathematically tractable concept. For instance, natural images are often sparse when represented in a wavelet basis, and objects in a classification task usually belong to only a small number of classes.

Due to the success of sparsity, a natural question is how we can refine the notion of sparsity in order to capture more complex structures. There are many examples where such an approach is applicable: (i) large wavelet coefficients of natural images tend to form connected trees, (ii) active genes can be arranged in functional groups, and (iii) approximate point sources in astronomical data often form clusters. In such cases, exploiting this additional structure can lead to improved compression ratio for images, better multi-label classification, or smaller sample complexity in compressive sensing and statistics. Hence an important question is the following: how can we model such sparsity structures, and how can we make effective use of this additional information in a computationally efficient manner?

There has been a wide range of work addressing these questions, e.g., [Yuan and Lin, 2006; Jacob et al., 2009; He and Carin, 2009; Kim and Xing, 2010; Bi and Kwok, 2011; Huang et al., 2011; Duarte and Eldar, 2011; Bach et al., 2012b; Rao et al., 2012; Negahban et al., 2012; Simon et al., 2013; El Halabi and Cevher, 2015]. Usually, the proposed solutions offer a trade-off between the following conflicting goals:

Generality What range of sparsity structures does the approach apply to?

Statistical efficiency What statistical performance improvements does the use of structure enable?

Computational efficiency How fast are the resulting algorithms?

In this paper, we introduce a framework for sparsity models defined through graphs, and we show that it achieves a compelling trade-off between the goals outlined above. At a high level, our approach applies to data with an underlying graph structure in which the large coefficients form a small number of connected components (optionally with additional constraints on the edges). Our approach offers three main features: (i) Generality: the framework encompasses several previously studied sparsity models, e.g., tree sparsity and cluster sparsity. (ii) Statistical efficiency: our sparsity model leads to reduced sample complexity in sparse recovery and achieves the information-theoretic optimum for a wide range of parameters. (iii) Computational efficiency: we give a nearly-linear time algorithm for our sparsity model, significantly improving on prior work both in theory and in practice. Due to the growing size of data sets encountered in science and engineering, algorithms with (nearly-)linear running time are becoming increasingly important.

We achieve these goals by connecting our sparsity model to the prize collecting Steiner tree (PCST) problem, which has been studied in combinatorial optimization and approximation algorithms. To establish this connection, we introduce a generalized version of the PCST problem and give a nearly-linear time algorithm for our variant. We believe that our sparsity model and the underlying algorithms are useful beyond sparse recovery, and we have already obtained results in this direction. To keep the presentation in this paper coherent, we focus on our results for sparse recovery and briefly mention further applications in Section 4.

We give an overview of our theoretical results in Section 2 and refer the reader to the full version of this paper for proofs and further details [Hegde et al., 2015b]. In Section 3, we complement our theoretical results with an empirical evaluation on both synthetic and real data (a background-subtracted image, a cerebral angiogram, and an image of text).

⇤Authors ordered alphabetically.

Basic notation Let [d] be the set ¶1↪ 2↪ ▷ ▷ ▷ ↪ d♦. We say that a vector β 2 Rd is s-sparse if at most s of its coefficients are nonzero. The support of β contains the indices corresponding to nonzero entries in β, i.e., supp(β) = ¶i 2 [d] ♣ βi 6= 0♦. Given a subset S ✓ [d], we write βS for the restriction of β to indices in S: we have (βS)i = βi for i 2 S and (βS)i = 0 otherwise. The `2-norm of β is kβk = qPi2[d] β2 i .

Sparsity models In some cases, we have more information about a vector than only "standard" s-sparsity. A natural way of encoding such additional structure is via sparsity models [Baraniuk et al., 2010]: let M be a family of supports, i.e., M = ¶S1↪ S2↪ ▷ ▷ ▷ ↪ SL♦ where Si ✓ [d]. Then the corresponding sparsity model M is the set of vectors supported on one of the Si:

2 Our contributions

We state our main contributions in the context of sparse recovery (see Section 4 for further applications). Here, the goal is to estimate an unknown s-sparse vector β 2 Rd from observations of the form

where X 2 Rn⇥d is the design matrix, y 2 Rn are the observations, and e 2 Rn is an observation noise vector. By imposing various assumptions on X and e, sparse recovery encompasses problems such as sparse linear regression and compressive sensing.

2.1 Weighted graph model (WGM)

The core of our framework for structured sparsity is a novel, general sparsity model which we call the weighted graph model. In the WGM, we use an underlying graph G = (V↪ E) defined on the coefficients of the unknown vector β, i.e., V = [d]. Moreover, the graph is weighted and we denote the edge weights with w : E ! N. We identify supports S ✓ [d] with subgraphs in G that are forests (unions of individual trees). Intuitively, the WGM captures sparsity structures with a small number of connected components in G. In order to control the sparsity patterns, the WGM offers three parameters:

• ÷ s, the total sparsity of S.
• ÷ g, the maximum number of connected components formed by the forest F corresponding to S.
• ÷ B, the bound on the total weight w(F) of edges in the forest F corresponding to S.

More formally, let γ(H) be the number of connected components in a graph H. Then we can define the WGM:

Definition 1. The (G↪ s↪ g↪ B)-WGM is the set of supports

Fig. 1 shows how two sparsity structures can be encoded with the WGM. Since our sparsity model applies to arbitrary graphs G, it can describe a wide range of structures. In particular, the model generalizes several previously studied sparsity models, including 1D-clusters, (wavelet) tree hierarchies, the Earth Mover Distance (EMD) model, and the unweighted graph model (see the full paper [Hegde et al., 2015b] for a detailed comparison).

2.2 Recovery of vectors in theWGM

We analyze the statistical efficiency of our framework in the context of sparse recovery. In particular, we prove that the sample complexity of recovering vectors in the WGM is provably smaller than the sample complexity for "standard" s-sparse vectors. To formally state this result, we first introduce a key property of graphs.

Definition 2. Let G = (V↪ E) be a weighted graph with edge weights w : E ! N. Then the weight-degree ⇢(v) of a node v is the largest number of adjacent nodes connected by edges with the same weight, i.e.,

We define the weight-degree of G to be the maximum weightdegree of any v 2 V .

Note that for graphs with uniform edge weights, the weight-degree of G is the same as the maximum node degree. Intuitively, the (weight) degree of a graph is an important property for quantifying the sample complexity of the WGM because the degree determines how restrictive the bound on the number of components g is. In the extreme case of a complete graph, any support can be formed with only a single connected component (see Figure 1). Using Definitions 1 and 2, we now state our sparse recovery result (see the full version of this paper for a more general statement of this result [Hegde et al., 2015b]).

Theorem 3. Let β 2 Rd be in the (G↪ s↪ g↪ B)-WGM. Then

i.i.d. Gaussian observations suffice to estimate β. More precisely, let e 2 Rn be an arbitrary noise vector and let y 2 Rn be defined as in Eq. 2 where X is an i.i.d. Gaussian matrix. Then we can efficiently find an estimate b β such that

where C is a constant independent of all variables above.

Note that in the noiseless case (e = 0), we are guaranteed to recover β exactly. Moreover, our estimate b β is in a slightly enlarged WGM for any amount of noise. Our bound (5) can be instantiated to recover previous sample complexity results, e.g., the n = O(s log ds ) bound for "standard" sparse recovery, which is tight [Do Ba et al., 2010]. For the image grid graph example in Figure 1, Equation (5) becomes n = O(s+g log dg ), which matches the information-theoretic optimum n = O(s) as long as the number of clusters is not too large, i.e., g = O(s◁ log d).

2.3 Efficient projection into theWGM

The algorithmic core of our sparsity framework is a computationally efficient procedure for projecting arbitrary vectors into the WGM. More precisely, the model-projection problem is the following: given a vector b 2 Rd and a WGMM, find the best approximation to b inM, i.e.,

If such a model-projection algorithm is available, one can instantiate the framework of [Baraniuk et al., 2010] in order to get an algorithm for sparse recovery with the respective sparsity model. However, solving Problem (7) exactly is NP-hard for the WGM due to a reduction from the classical Steiner tree problem. To circumvent this hardness result, we use the approximation-tolerant framework of [Hegde et al., 2015a]. Instead of solving (7) exactly, the framework requires two algorithms with the following complementary approximation guarantees.

Tail approximation: Find an S 2 M such that

Head approximation: Find an S 2 M such that

Here, cT > 1 and cH < 1 are arbitrary, fixed constants. Note that a head approximation guarantee does not imply a tail guarantee (and vice versa). In fact, stable recovery is not possible with only one type of approximate projection guarantee [Hegde et al., 2015a]. We provide two algorithms for solving (8) and (9) (one per guarantee) which both run in nearlylinear time.

Our model-projection algorithms are based on a connection to the prize-collecting Steiner tree problem (PCST), which is a generalization of the classical Steiner tree problem. Instead of finding the cheapest way to connect all terminal nodes in a given weighted graph, we can instead omit some terminals from the solution and pay a specific price for each omitted node. The goal is to find a subtree with the optimal trade-off between the cost paid for edges used to connect a subset of the nodes and the price of the remaining, unconnected nodes (see [Goemans and Williamson, 1995] for a formal definition).

We make the following three main algorithmic contributions. Due to the wide applicability of the PCST problem, we believe that these algorithms can be of independent interest (see Section 4).

• ÷ We introduce a variant of the PCST problem in which the goal is to find a set of g trees instead of a single tree. We call this variant the prize-collecting Steiner forest (PCSF) problem and adapt the algorithm of [Goemans andWilliamson, 1995] for this variant.
• ÷ We reduce the projection problems (8) and (9) to a small set of adaptively constructed PCSF instances.
• ÷ We give a nearly-linear time algorithm for the PCSF problem and hence also the model projection problem.

2.4 Comparison to related work

In addition to "point-solutions" for individual sparsity models, there has been a wide range of work on general frameworks for utilizing structure in sparse recovery. The approach most similar to ours is [Baraniuk et al., 2010], which gives a framework underlying several recovery algorithms for individual sparsity models. However, the framework has one important drawback: it does not come with a full recovery algorithm. Instead, the authors only give a recovery scheme that assumes the existence of a model-projection algorithm satisfying (7). Such an algorithm must be constructed from scratch for each model, and the techniques that have been used for various models so far are quite different. Our contribution can be seen as complementing the framework of [Baraniuk et al., 2010] with a nearly-linear time projection algorithm that is applicable to a wide range of sparsity structures. This answers a question raised by the authors of [Huang et al., 2011], who also give a framework for structured sparsity with a universal and complete recovery algorithm. Their framework is applicable to a wide range of sparsity models, but the corresponding algorithm is significantly slower than ours, both in theory as well as in practice, as demonstrated by our experiments in Section 3. Moreover, our recovery algorithm shows more robust performance across different shapes of graph clusters.

Both of the approaches mentioned above use iterative greedy algorithms for sparse recovery. There is also a large body of work on combining M-estimators with convex regularizers that induce structured sparsity, e.g., see the surveys [Bach et al., 2012a] and [Wainwright, 2014]. The work closest to ours is [Jacob et al., 2009], which uses an overlapping group Lasso to enforce graph-structured sparsity (graph Lasso). In contrast to their approach, our algorithm gives more fine-grained control over the number of clusters in the graph. Moreover, our algorithm has better computational complexity, and to the best of our knowledge there are no formal results relating the graph structure to the sample complexity of the graph Lasso. Empirically, our algorithm recovers an unknown vector with graph structure faster and from fewer observations than the graph Lasso (see the supplementary material of the full paper [Hegde et al., 2015b]).

3 Experiments

We focus on the performance of our proposed algorithm (Graph-CoSaMP) for the task of recovering 2D data with clustered sparsity. Multiple methods have been proposed for this problem, and our theoretical analysis shows that our algorithm should improve upon the state of the art. We compare our results to StructOMP [Huang et al., 2011] and the heuristic Lattice Matching Pursuit (LaMP) [Cevher et al., 2009]. The implementations were supplied by the authors and we used the default parameter settings. Moreover, we ran two common recovery algorithms for "standard" s-sparsity: Basis Pursuit [Candès et al., 2006] and CoSaMP [Needell and Tropp, 2009].

We follow a standard evaluation procedure for sparse recovery / compressive sensing: we record n observations y = Xβ of the (vectorized) image β 2 Rd using a subsampled Fourier matrix X. We assume that all algorithms possess prior knowledge of the sparsity s and the number of connected-components g in the true support of the image β. We declare a trial successful if the squared `2-norm of the recovery error is at most 5% of the squared `2-norm of the original vector β. The probability of successful recovery is then estimated by averaging over 50 trials. We perform several experiments with varying oversampling ratios n◁s and three different images. See the supplementary material of the full paper [Hegde et al., 2015b] for a description of the dataset, experiments with noise, and a comparison with the graph Lasso.

Figure 2 demonstrates that Graph-CoSaMP yields consistently competitive phase transitions and exhibits the best sample complexity for images with "long" connected clusters, such as the angiogram image (b) and the text image (c). While StructOMP performs well on "blob"-like images such as the background-subtracted image (a), its performance is poor in our other test cases. For example, it can successfully recover the text image only for oversampling ratios n◁s > 15. Note that the performance of Graph-CoSaMP is very consistent: in all three examples, the phase transition occurs between oversampling ratios 3 and 4. Other methods show significantly more variability.

We also investigate the computational efficiency of GraphCoSaMP. We consider resized versions of the angiogram image and record n = 6s observations for each image size d. Figure 2(d) displays the recovery times (averaged over 50 trials) as a function of d. We observe that the runtime of Graph-CoSaMP scales nearly linearly with d, comparable to the conventional sparse recovery methods. Moreover, GraphCoSaMP is about 20⇥ faster than StructOMP.

4 Further applications

We have introduced a general framework for structured sparsity that encompasses several previously studied sparsity models, but still allows information-theoretically optimal sample complexity and nearly-linear time algorithms.

Our algorithms have found applications beyond sparse recovery. In [Schmidt et al., 2015], we have used our PCST algorithm for a feature extraction task in seismic image processing. Other potential applications include the work of [Rozenshtein et al., 2014], who use a PCST algorithm for event detection in social networks.

Acknowledgements

We thank Jayadev Acharya, Stefanie Jegelka, Youssef Mroueh, Devavrat Shah, and the anonymous reviewers for many helpful comments on earlier versions of this paper. This work was supported by grants from the MITEI-Shell program, the MADALGO center, and the Simons Investigator Award.