A Fast Algorithm for Separated Sparsity via Perturbed Lagrangians

arXiv:1712.08130v1 [cs.DS] 21 Dec 2017

Aleksander Madry
˛
MIT
madry@mit.edu

Slobodan Mitrović
EPFL
slobodan.mitrovic@epfl.ch

Ludwig Schmidt
MIT
ludwigs@mit.edu

Abstract
Sparsity-based methods are widely used in machine learning, statistics, and signal processing. There
is now a rich class of structured sparsity approaches that expand the modeling power of the sparsity
paradigm and incorporate constraints such as group sparsity, graph sparsity, or hierarchical sparsity. While
these sparsity models offer improved sample complexity and better interpretability, the improvements
come at a computational cost: it is often challenging to optimize over the (non-convex) constraint sets
that capture various sparsity structures. In this paper, we make progress in this direction in the context of
separated sparsity – a fundamental sparsity notion that captures exclusion constraints in linearly ordered
data such as time series. While prior algorithms for computing a projection onto this constraint set
required quadratic time, we provide a perturbed Lagrangian relaxation approach that computes provably
exact projection in only nearly-linear time. Although the sparsity constraint is non-convex, our perturbed
Lagrangian approach is still guaranteed to find a globally optimal solution. In experiments, our new
algorithms offer a 10× speed-up already on moderately-size inputs.

1

1

Introduction

Over the past two decades, sparsity has become a widely used tool in several fields including signal processing,
statistics, and machine learning. In many cases, sparsity is the key concept that enables us to capture important
structure present in real-world data while making the resulting problem computationally tractable and suitable
for mathematical analysis. Among the many applications of sparsity are sparse linear regression, compressed
sensing, sparse PCA, and dictionary learning.
The first wave of sparsity-based techniques focused on the standard notion of sparsity that only constrains
the number of non-zeros. Over time, it became apparent that extending the notion of sparsity to encompass
more complex structures present in real-world data can offer significant benefits. Specifically, utilizing such
additional structure often improves the statistical efficiency in estimation problems and the interpretability of
the final result. There is now a large body of work on structured sparsity that has introduced popular models
such as group sparsity and hierarchical sparsity [1, 2, 3, 4, 5, 6, 7, 8, 9]. These statistical improvements,
however, come at a computational cost: the resulting optimization problems are often much harder to
solve. The key reason is that the combinatorial sparsity structures give rise to non-convex constraints.
Consequently, many of the resulting algorithms have significantly worse running time than their “standard
sparsity” counterparts. This trade-off raises an important question: can we design algorithms for structured
sparsity that match the time complexity of commonly used algorithms for standard sparsity?
In this paper, we address this question in the context of the separated sparsity model, a popular sparsity
model for data with a known minimum distance between large coefficients [10, 11, 12, 13, 14]. In the
one-dimensional case, such as time series data, neuronal spike trains are a natural example. Here, a minimum
refractory period ensures separation between consecutive spikes. In two dimensions, separation constraints
arise in the context of astronomical images or super-resolution applications [12, 15].
We introduce new algorithms for separated sparsity that run in nearly-linear time. This significantly
improves over prior work that required at least quadratic time, which is quickly prohibitive for large data
sets. An important consequence of our fast running time is that it enables methods that utilize separated
sparsity yet are essentially as fast as their counterparts based on standard sparsity only. For instance, when
we instantiate our algorithm in compressive sensing, the running time of our method matches that of common
methods such as IHT or CoSaMP.
Our algorithms stem from a primal-dual linear programming (LP) perspective on the problem. Our
theoretical findings reveal a rich structure behind the separated sparsity model, which we utilize to obtain
efficient methods. Interestingly, our final algorithm has a very simple form that can be interpreted as a
Lagrangian relaxation of the sparsity constraint. In spite of the non-convexity of the constraint, our algorithm
is still guaranteed to find the globally optimal solution.
We also show that these algorithmic and theoretical contributions directly translate into empirical
efficiency. Specifically, we demonstrate that, compared to the state of the art procedures, our methods yield
an order of magnitude speed-up already on moderate-size inputs. We run experiments on synthetic data and
real world neuronal spike train signals.

2

Separated sparsity and applications

In this section, we formally define separated sparsity and the corresponding algorithmic problems. As a
concrete application of our algorithms, we instantiate them in a sparse recovery context that is representative
for many statistical problems such as compressed sensing and sparse linear regression.
First, we briefly introduce our notation. As usual, [d] denotes the set {1, . . . , d}. We say that a vector

1

θ ∈ Rd is k-sparse if θ contains at most k non-zero coefficients. We define the support of θ as the set of
indices corresponding to non-zero coefficients, i.e., supp(θ) = {i ∈ [d] | θi 6= 0}. We let kθk denote the
`2 -norm of a vector θ ∈ Rd .
Separated sparsity. Sparsity models are a natural way to formalize structure beyond “standard” sparsity
[3]. In this work we focus on the the separated sparsity model, defined as follows [10]. For a support Ω ⊆ [d],
let sep(Ω) = mini6=j∈Ω |i − j| be the minimum separation of two indices in the support. We define the
following two sets of supports: set M∆ = {Ω ⊆ [d] | sep(Ω) ≥ ∆}, and ∆-separated sparsity supports
Mk,∆ = {Ω ∈ M∆ | |Ω| = k}. That is, Mk,∆ is the set of support patterns containing k non-zeros with at least
∆ − 1 zero entries between consecutive non-zeros.
In order to employ separated sparsity in statistical problems, we often want to add constraints based
on the support set Mk,∆ to optimization problems such as empirical risk minimization. A standard way
of incorporating constraints into gradient-based algorithms is via a projection operator. In the context of
separated sparsity, this corresponds to the following problem.
Problem 1 For a given input vector x ∈ Rd , our goal is to project x onto the set Mk,∆ , i.e., to find a vector x̂
such that
kx − x0 k .
(1)
x̂ ∈
arg min
x0 ∈Rd : supp(x0 )∈Mk,∆

Problem 1 is the main algorithmic problem we address in this paper.
Sparse recovery. Structured sparsity has been employed in a variety of machine learning tasks. In order to
keep the discussion coherent, we present our results in the context of the well-known sparse linear model:
y = Xθ∗ + e

(2)

where y ∈ Rn are the observations/measurements, X ∈ Rn×d is the design or measurement matrix, and
e ∈ Rn is a noise vector. The goal is to find a good estimate θ̂ of the unknown parameters θ∗ up to the noise
level.
The authors of [3] give an elegant framework for incorporating structured sparsity into the estimation
problem outlined above. They design a general recovery algorithm that relies on a model-specific projection
oracle. In the case of separated sparsity, this oracle is required to solve precisely the Problem 1 stated above.

3

The algorithm and our results

Given an arbitrary vector x ∈ Rd , Problem 1 requires us to find a vector x̂ such that x̂ ∈ Mk,∆ and kx − x̂k
is minimized. We now slightly reformulate the problem. Let c ∈ Rd be a vector such that ci = x2i for all
i. Then it is not hard to see that this problem is equivalent to finding a set of k entries in the vector c such
that each of these entries is separated by at least ∆ and the sum of these entries is maximized. Hence our
main algorithmic problem is to find a set of k entries in an non-negative input vector c so that the entries are
∆-separated and their sum is maximized. More formally, our goal is to find a support Ŝ such that
X
Ŝ ∈ arg max
ci .
(3)
S∈Mk,∆ i∈S

2

In the following, we also consider a relaxed version of (3) called P ROJ L AGR, which is parametrized by a
trade-off parameter λ and a vector c̃:
X
P ROJ L AGR(λ, c̃) := arg max
c̃i + λ (k − |S|) .
S∈M∆

i∈S

Intuitively, P ROJ L AGR represents a Lagrangian relaxation of the sparsity constraint in Equation (3).

3.1

Algorithm

Our main contribution is a new algorithm for Problem 1 that we call Lagrangian Approach to the Separated
Sparsity Problem (LASSP). The pseudo code is given in Algorithm 1.
LASSP is a Las Vegas algorithm: it always returns a correct answer, but the running time of the algorithm
is randomized. Concretely, LASSP repeats a main loop until a stopping criterion is reached. Every iteration
of LASSP first adds a small perturbation to the coefficients c (see Line 3). This perturbation has only a
small effect on the solution but improves the “conditioning” of the corresponding non-convex Lagrangian
relaxation P ROJ L AGR so that it returns a globally optimal solution that almost satisfies the constraint. As
we show in Section 5.2, we can solve this relaxation (Line 4) in nearly-linear time. After the algorithm has
solved the Lagrangian relaxation, it obtains the final support Ŝ in line 5 by solving P ROJ L AGR on a slightly
shifted λ̂ to ensure that the constraint is satisfied with good probability.
Remark: We assume that the bit precision γ required to represent the coefficients c is finite, and we
provide our results as a function of γ. For practical purposes, γ is usually a constant. Since the solution to
Equation (3) is invariant under scaling by a positive integer and γ is finite, without loss of generality we
assume that c ∈ Zd .
Algorithm 1 LASSP:
Input: c ∈ Zd , k ∈ N+
Output: A solution Ŝ to (3)
1: repeat
2:
Let X ∈ Zd be a vector such that Xi is chosen uniformly at random from {0, . . . , d3 − 1}, ∀i
3:
Define vector c̃ := d4 c + X
4:
Let λ̂ := arg minλ∈Z P ROJ
 L AGR(λ, c̃).
 In case of ties, maximize λ̂.
5:

Choose Ŝ ∈ P ROJ L AGR λ̂ −

6:

until Ŝ = k

7:

return Ŝ

3.2

1
d+1 , c̃

Main results

As our main result, in the following theorem we show that LASSP runs in nearly linear time and solves
Problem 1.
Theorem 1 Let c ∈ Rd , and let γ ∈ N+ be the maximal number of bits needed to store any ci . There is an
implementation of LASSP that for every c computes a solution ĉ to Problem 1. With probability 1 − 1/d, the
algorithm runs in time O(d(γ + log d)).

3

Combined with the framework of [3], we get the following.
Corollary 1 Let y, X, θ∗ , and e be as in the sparse linear model in Equation (2). We assume that supp(θ∗ ) ∈
Mk,∆ and that X satisfies the model-RIP for Mk,∆ . There is an algorithm that for every y and e returns an
estimate θ̂ such that
kθ̂ − θ∗ k2 ≤ Ckek2 .
e X + d), where TX is the time of multiplying the matrices X and
Moreover, the algorithm runs in time O(T
T
X by a vector.
The corollary shows that, up to logarithmic factors, the running time is dominated by TX + d. This
matches the time complexity of standard sparse recovery and shows that we can utilize separated sparsity
without a significant increase in time complexity. Many measurement matrices in compressive sensing enable
fast multiplication with X (e.g., a subsampled Fourier matrix), in which case the total running time becomes
e
O(d).
We validate these theoretical findings in Section 7 by showing that LASSP runs significantly faster
than the state of the art algorithm used for sparse recovery with separation constraints, while retaining the
same accuracy of the recovered signal.
Our algorithm LASSP is randomized. However, we also design a deterministic nearly-linear time
algorithm and prove the following theorem. For clarity of exposition, the statement of the deterministic
algorithm and the proof of the theorem are deferred to Appendix J.
Theorem 2 Let c ∈ Rd be the input vector and let γ be as in Theorem 1. Then there is an algorithm that
computes a solution ĉ satisfying Equation (1) and runs in time O(d(γ + log k) log d log ∆).
Remark. The ∆-separated sparsity projection can be reduced to the problem of finding a minimum-weight
path of length k on an edge-weighted directed graph. After this work was done, it was pointed to us by Arturs
Backurs and Christos Tzamos that, by properly designing the graph, the edge-weights of the corresponding
directed graph satisfy the concave Monge property. This further implies that the separated-sparsity problem
can be solved in nearly-linear time [16]. However, the algorithm presented in [16] is more complex than the
algorithm we provide. The goal of our work was to design a simple, and yet very efficient, algorithm by
applying minor but crucial modifications to a widely popular approach, i.e. to Lagrangian relaxation.

3.3

Further results

In Appendix C we present a dynamic programming approach that for a specific, but also natural, family of
instances solves the separated sparsity problem in even linear time.
We also consider a natural extension to the 2D-variant of the separated sparsity projection problem and
show that it is NP-hard in Appendix D. Moreover, in Appendix E, we extend our model to allow for blocks of
separated variables and show that our algorithms for Mk,∆ also applies to the more general variant. Finally,
separated sparsity can be used to model signals in which a longer pattern is repeated multiple times so that
any two patterns are at least ∆ apart. This model is called disjoint pulse streams [12]. Again, the algorithmic
core remains the same and algorithms for Mk,∆ can also be used for this generalization.

3.4

Additional related work

The papers [10, 14] are closely related to our work. The paper [10] proposed the separated sparsity model,
provided a sample complexity upper bound, and gave an LP-based model-projection algorithm. However, they
4

resorted to a black-box approach for solving the LP, that lead to a fairly prohibitive O(d3.5 ) time complexity.
Recently, [14] provided a faster dynamic program for this problem with a time complexity of O(d2 ) and also
showed a sample complexity lower bound. The algorithmic aspect of these papers is the main difference from
our work: we exploit structure in both the primal and dual formulations of the LP and give an algorithm that
provably runs in nearly-linear time.
Beside the papers addressing the core algorithmic question of projecting onto separated sparse vectors,
there is much of work utilizing the sparsity model for applications in neural signal processing [11, 12, 13]
and recovery with coherent dictionaries [17, 18]. In the latter application, the separated sparsity constraint
enforces that the signal representation only consists of incoherent dictionary atoms. We expect that our
algorithmic techniques will also lead to improvements in the context of these applications.
In addition to separated sparsity, a large body of work on structured sparsity has emerged over the past few
years. We refer the reader to the surveys [19, 20, 21, 22] for an overview. The line of work most relevant to
our paper is the model-based compressive sensing framework introduced in [3], which is also the starting point
for [10, 14]. While the framework provides a general recovery scheme based on a model-projection oracle
satisfying Equation (1), it does not provide any guidance on how to design such oracles for a specific sparsity
model. We address precisely this problem for separated sparsity with our nearly-linear time algorithm.
Recently, two papers have proposed a fairly general framework for deriving model-projection oracles
via graph sparsity, i.e., sparsity structures that can be defined through connected components in a graph on
the signal coefficients [5, 9]. This framework generalizes several previously studied sparsity models such
as block sparsity, tree sparsity, and cluster sparsity. Moreover, the paper [9] gives projections that run in
nearly-linear time. Interestingly, these tools do not apply to the separated sparsity model we study in our
paper. Intuitively, graph sparsity captures structures in which the non-zero coefficients are clustered together,
while the separated sparsity model achieves a reduction in sample complexity for the opposite reason: the
non-zero coefficients are far apart. Moreover, the algorithms in [9] are approximate and project into a
sparsity model with a relaxed sparsity constraint. As we explain in Appendix B, the separated sparsity model
requires more careful control over the output sparsity in order to achieve a meaningful sample complexity
improvement over “standard” sparsity. We circumvent this issue by providing an exact projection onto the
separated sparsity model.

4

Proof of correctness and the roadmap

We begin our analysis by proving that LASSP returns a correct result if the algorithm terminates. As we will
see later, establishing termination is the crucial part of the analysis. Nevertheless, the following lemma is a
useful warm-up for understanding how the different pieces of our algorithm fit together.
Lemma 1 When LASSP terminates, it outputs a support Ŝ such that x restricted to Ŝ is a solution to
Problem 1.
Proof. As we have argued above, the problems in Equations (1) and (3) are equivalent. So, we show that
LASSP outputs a solution to the problem in Equation (3).
Let Ŝ be the set returned by LASSP. By the condition of the loop in Line 6, we have |Ŝ| = k. So, as
Ŝ ∈ M∆ (see the definition of P ROJ L AGR), we have Ŝ ∈ Mk,∆ .
Now, towards a contradiction, assume
P that support
PŜ is not a solution to the problem in Equation (3),
?
while support S is. This implies that i∈S ? ci > i∈Ŝ ci . Now since, without loss of generality, we

5

assumed that c ∈ Zd , the last inequality implies
X

P

i∈S ? ci

d4 ci ≥ d4 +

i∈S ?

≥1+
X

P

i∈Ŝ ci ,

and hence

d4 ci .

(4)

i∈Ŝ

Observe that for any support S ∈ Mk,∆ , the term λ (k − |S|) equals zero, and recall that Ŝ, S ? ∈ Mk,∆ .
Furthermore, by the definition of the random vector X in line 2 and from (4)
X
X
X

d4 ci + Xi ≥
d4 ci ≥ d4 +
d4 ci
i∈S ?

i∈S ?

>

X

i∈Ŝ
4



d ci + Xi .

i∈Ŝ

Since S ? ∈ M∆ , this chain
of inequalities contradicts
Line 5 of LASSP which chooses Ŝ as an optimal


4
solution to P ROJ L AGR λ̂ − 1/(d + 1), d c + X . This further implies that Ŝ is a solution to Problem 1. 2

4.1

Roadmap

Lemma 1 shows that LASSP outputs the right answer if it terminates. But does LASSP terminate on every
input? Answering this question is the most intricate part of this paper. We split the proof in two main pieces.
The first part is Section 5, where we provide an alternative view on the separated sparsity problem based on
linear programming duality. The duality view paves the way towards proving our main results. In particular,
we show that the subroutines in LASSP can be implemented quickly.
Lemma 2 Single iteration of LASSP can be implemented to run in time O(d(γ + log d)).
The second part is Section 6, where we further study the duality view on separated sparsity. We show that
after perturbing the input instance in Line 3 of LASSP, the support obtained with a shifted λ̂ in Line 5 has
cardinality k with high probability.
Lemma 3 Algorithm LASSP runs only a single iteration with probability at least 1 − 1/d.
Together with Lemma 1, these results yield Theorem 1.

5

Part I – To duality and further

We now analyze the running time of a single iteration of LASSP. We provide a series of equivalences, as
illustrated in Figure 1, in order to exploit structure in the separated sparsity problem. More precisely, we start
with a linear programming (LP) view on separated sparsity. It has already been shown that this viewpoint
yields a totally unimodular LP [10], which implies that the LP has an integral solution. Hence solving the LP
solves the separated sparsity projection in Problem 1. However, prior work did not utilize this connection to
reason about the power of the Lagrangian relaxation approach to the problem.
We begin our detailed analysis of the LP with the dual program D. By strong duality, the value of D
equals the value of the primal LP. Then we cast D as minimization of LP Dλ over λ. This reduction will play
the central role in our analysis and connect Dλ to Line 4 of LASSP.

6

5.1

The LP Perspective

We start with a linear programming view on problem (3) by considering its LP relaxation denoted by P:
cT u

maximize

d
X

subject to

ui = k

i=1
min{i+∆−1,n}

X

uj ≤ 1

∀i = 1 . . . d

ui ≥ 0

∀i = 1 . . . d

j=i

Given an LP A we use VAL A to denote its optimal objective value. As already noted, P is totally unimodular
and thus there always exists an optimal solution to it that is integral [23]. For completeness, we provide the
proof of total unimodularity in Appendix F.
Remark: This implies that a solution to P can be used to obtain a solution to the separated sparsity model
projection: if u? is an optimal solution to P, we can derive the optimal support of (3) from the non-zero
entries among the LP variables u?1 , . . . , u?d . It is unclear, however, if there is a way to directly solve this LP
fast, e.g. it is not known how to solve P directly in time matching the running time of our algorithm LASSP.
A key step in our approach is understanding the
separated sparsity structure from the dual point of
TUM
view. The dual LP to P, denoted by D, is given as
Problem (3)
P
LP P
max
c
S2Mk;∆
follows
i2S i
minimize

w0 k +

subject to

X

d
X

,

,

wi

i=1

dual LP D

wj ≥ ci

∀i = 1 . . . d

j : j≥1 and
j≤i≤j+∆−1

,

w0 +

strong duality

our algorithm

wi ≥ 0

∀i = 1 . . . d

convexity
+ TUM
minλ2Z Dλ

,

w0 ∈ R
Then, as P is integral, so is D.
line 4 of LASSP

TUM

,

strong duality

min

P

(λ)

λ2Z no k
minλ2Z ProjLagr(λ; c)
Corollary 2 For an integer k and a vector of integers c, there exists ŵ such that ŵ is an optimum of
Figure 1: The values of the problems in the diagram
D and ŵ0 ∈ Z.
are equal. Every equivalence relation carries structural
We also define Dλ as the LP D in which the information that we utilize in our analysis.
variable w0 is set to λ. Now, it is not hard to show
the following lemma, whose proof is deferred to
Lemma 10 in Appendix F.

Lemma 4 Dλ is convex with respect to λ.

7

From the definition of D, the following equality holds VAL D = minλ∈R VAL Dλ . Furthermore, Corollary 2
implies that it is sufficient to consider λ in Z only, i.e.
VAL D = min VAL Dλ .

(5)

λ∈Z

Now, Lemma 4 implies that we can obtain VAL D by applying ternary search over λ on function VAL Dλ .

5.2

Implementing one iteration of LASSP efficiently

We now derive the final connection between D and LASSP which will enable us to obtain λ̂ at line 4 in
nearly-linear time. To that end, consider Pno-k (λ) defined as
maximize

cT u + λ(k − 1T u)
min{i+∆−1,d}

subject to

X

uj ≤ 1

∀i = 1 . . . d

ui ≥ 0

∀i = 1 . . . d

j=i

Observe that compared to P, Pno-k (λ) does not contain the sparsity constraint. Furthermore, the LP Pno-k (λ)
is a relaxed version of P ROJ L AGR (λ, c). Also, as P is, then Pno-k (λ) is totally unimodular. Now this
sequence of conclusions results in the following.
Corollary 3 Problems Pno-k (λ) and P ROJ L AGR (λ, c) are equivalent.
To obtain the final connection, we consider LP given by
LP := min Pno-k (λ).
λ∈R

Next, note that the dual of Pno-k (λ) is Dλ . Hence, the strong duality implies VAL Dλ = VAL Pno-k (λ). That
together with VAL D = minλ∈R VAL Dλ yields that D and LP coincide as functions in λ. Furthermore, since
we can solve D by applying ternary search over integral values of λ and Dλ , we can solve LP by applying
ternary search over integral values of λ and function Pno-k (λ). But since Pno-k (λ) and P ROJ L AGR (λ, c) are
equivalent, we can also obtain λ̂ at line 4 of LASSP by applying ternary search over λ.
Now, it is very easy to see that for an optimal solution w? of D we have w0? ≤ maxi |ci |. It is also not
hard to show that there is an optimal solution such that w0? ≥ −(k − 1) maxi |ci | (see Lemma 18). Therefore,
in order to find optimal λ̂ it suffices to execute O (log maxi |ci | + log k) = O(γ + log d) iterations of ternary
search.
Every iteration of the ternary search invokes P ROJ L AGR, which can be implemented to run in linear time.
Lemma 5 Given λ and ĉ ∈ Rd , there is an algorithm that finds support Ŝ ∈ P ROJ L AGR (λ, ĉ) in time O(d).
0
Proof. Observe
P that for a fixed λ, solving P ROJ L AGR (λ, ĉ) is equivalent to finding support S ∈ M∆ that
maximize i∈S 0 (ĉi − λ). Hence, we can reinterpret P ROJ L AGR (λ, ĉ) as follows: given a vector c̃ := ĉ−λ1,
select a subset of [d] of indices (not necessarily k of them) so that (i) every two indices are at least ∆ apart,
and (ii) the sum of the values of c̃ at the selected indices is maximized.
This task can be solved by standard dynamic programming in the following way. For every i, we define
si to be the maximum value of the described task restricted to the first i indices of c̃. Then, it is easy to

8

see that si+1 = max{si , si+1−∆ + c̃i+1 }. Namely, we can either decide not to select index i + 1, in which
case the best is already contained in si ; or, we can decide to select index i + 1 which has value c̃i+1 and for
the rest we consider si+1−∆ . Therefore, sd can be obtained in time O(d). Now it is easy to reconstruct the
corresponding support in linear time. For completeness, we provide a full algorithm and a detailed proof in
Lemma 12, Appendix G. 2
Putting all together proves Lemma 2.

6

Part II – Active constraints

Note that the chain of equivalences present in Figure 1 shows that minλ P ROJ L AGR(λ, c) outputs the sum of
coordinates of an optimal solution of problem 1. Prior to our work, it was not even known how to obtain
this value in time faster than O(dk), while our result shows we can compute it in nearly-linear time. So, it is
natural to ask whether the same relaxation also outputs a support of size k? The answer is, unfortunately, no.
To see that, consider the example: c = (4, 7, 5, 0, 0, 5, 8, 5), ∆ = 2, and k = 3. For λ > 2 every solution u?λ
to P ROJ L AGR(λ, c) is such that u?λ has less than k non-zeros. On the other hand, for every λ ≤ 2 the exists a
solution u?λ such that u?λ contains more than k non-zeros. Therefore, there is no λ, neither λ − 1/(d + 1), for
which P ROJ L AGR(λ, c) provably outputs a support of cardinality k. This also suggests that the perturbation
we apply in lines 2-3 is essential!
Instead of studying lines 2-5 of LASSP directly, we shift our focus to Dλ . In particular, we exhibit very
close connection between its structure and the sparsity of the primal solution, which we present via the notion
of "active constraints". Then we use these findings in our analysis to show that slight perturbation of the input
instance, while not affecting the value of the solution, makes it possible to obtain a solution to problem 3 by
applying Lagrangian relaxation.

6.1

Solving Dλ

Observe that once we fixed the value of w0 , all remaining constraints in Dλ are “local” since they only
affect a known interval of length ∆. They are also ordered in a natural way. As a result, we can solve Dλ
by making a single pass over these variables. Starting with w1 and all variables set to 0, we consider each
constraint from left to right and increase the variables to satisfy these constraints in a lazy manner. That is, if
in our pass we reach a constraint with index i that is still not satisfied, we increase the value of wi until that
constraint becomes satisfied and then move to the next constraint. Given c and λ, algorithm Dual-Greedy, i.e.
Algorithm 2, formalizes this approach whose analysis appears in Appendix F.

6.2

Tracking the change of Dλ

Next we introduce the key concept that we need for relating the solution of dual to the sparsity of Lagrangian
relaxation of the primal: the notion of active constraints. Let w be a vector obtained by Dual-Greedy(c, λ)
and let w0 be a vector obtained by Dual-Greedy(c, λ − ε), for some fixed λ ∈ Z and small ε ∈ (0, 1). Then,
the set of coordinates that are for ε larger in w0 than in w are called active constraints. Figure 2 provides an
illustration of this concept. Intuitively, the active constraints correspond to those variables of D that increase
when w0 decreases by some small value. Hence, one can interpret active constraints as gradients of Dλ with
respect to the variable λ. This concept appears to be very useful in characterizing the optimal solution of D
in an alternative way. In particular, the following lemma holds.

9

Algorithm 2 Dual-Greedy:
Input: c ∈ Zd , λ ∈ R
Output: an optimal solution w to D such that w0 = λ
1: w ← 0;
w0 ← λ
# initialize the output
2: sum ∆ ← 0
# store the sum of the most recent ∆ variables of w
3: for i := 1 . . . d do
4:
if i − ∆ ≥ 1 then sum ∆ ← sum ∆ − wi−∆
5:
diff ← ci − (w0 + sum ∆ )
# compute "how far" constraint i is from being tight
6:
if diff > 0 then wi ← diff
# if constraint i is not satisfied, make it tight
7:
sum ∆ ← sum ∆ + wi
# update the sum of the ∆ most recent variables of w
8: return w
# solution to D s.t. w0 = λ

c1 = 4 c2 = 6 c3 = 4 c4 = 3 c5 = 3 c6 = 5 c7 = 6 c8 = 2 c9 = 2 c10 = 2

c1 = 4 c2 = 6 c3 = 4 c4 = 3 c5 = 3 c6 = 5 c7 = 6 c8 = 2 c9 = 2 c10 = 2
w0 = 2 − ε

w0 = 2
w5 = 1

w1 = 2

w10 = ε

w5 = 1 + ε

w1 = 2 + ε

w6 = 2

w6 = 2

w2 = 2

w2 = 2

w7 = 1

w7 = 1

Figure 2: A sketch of an instance of D: c = (4, 6, 4, 3, 3, 5, 6, 2, 2, 2) and ∆ = 3. The left and right figure
depicts w obtained by Dual-Greedy(c, 2) and Dual-Greedy(c, 2 − ε) for 0 < ε < 1, respectively. Constraints
1, 5, and 10 are active, i.e., if w0 on the left is decreased by ε then only w1 , w5 , and w10 increase by ε, as
shown on the right. Every wi , for i ≥ 1, covers ∆ c-poles to its right.
Lemma 6 Let c ∈ Zd , λ ∈ Z, and w ← Dual-Greedy(c, λ). Then, if w has exactly k active constraints the
vector w is an optimal solution to D.
A full proof of a statement stronger than Lemma 6 along with its proof appears in Lemma 13, Appendix H,
while in this section we provide a proof sketch. Let w(ε) = Dual-Greedy(c, λ − ε), for some small ε ∈ (0, 1).
By the definition of active constraints and the fact that w has k many, there are exactly k coordinates that
are larger by ε in w(ε) than in w. In addition, w(ε)0 = λ − ε and w0 =P
λ. It is not hard to show
Pd that all the
d
other coordinates of w(ε) and w are the same, which we can express as i=1 w(ε)i = kε + i=1 wi . Now,
P
recall that the objective function of dual D with respect to vector w equals w0 k + di=1 wi . Then we have
w0 k +

d
X

wi = (w0 − ε)k + kε +

i=1

d
X

wi = w(ε)0 k +

i=1

d
X

w(ε)i .

i=1

Hence, the objective values of dual D for vectors w(ε) and w are equal, for all the values ε ∈ (0, 1). As
D is convex in the value of variable w0 and Dual-Greedy(c, λ) provides an optimal solution to D such that
w0 = λ, then w is an optimal solution to D.

6.3

Wrapping up – perturbation and optimal sparsity

Now we use Lemma 6 to prove the following, which essentially justifies line 5 of LASSP.
10

Lemma 7 Let c ∈ Zd , λ ∈ Z, and w̃ ← Dual-Greedy(c, λ). Assume that w̃ has exactly k active constraints.
Then, any optimal support S ? of P ROJ L AGR(λ − ε, c), for 0 < ε < 1/d, has cardinality exactly k.
So, if we produce λ as in Lemma 7, we will solve problem (3). These steps are implemented by lines 4-6 of
LASSP. However, as illustrated in the beginning of the section, λ as in Lemma 7 might not exist. Intuitively,
this situation happens when the number of active constraints jumps from a value smaller than k to a value
larger than k for a very small change of λ. In such a case, we are unable to obtain λ as in Lemma 7. A key
component of our analysis is showing that there is an efficient way of randomly altering c, and obtaining c̃,
so that with high probability c̃ is such that: λ̂ is obtained as at line 4; and, w̃ ← Dual-Greedy(c̃, λ̂) has the
same property as in Lemma 7. Lines 2 and 3 of LASSP implement this random perturbation. Intuitively, the
perturbation achieved by X variables adds noise to our input instance so that the number of active constraints
changes by at most one as λ slides over the integer domain. Following this intuition we obtain a proof of
Lemma 3. However, due to the space limitation, we present its proof in Appendix I.
We conclude the section by giving a proof of Lemma 7.
Proof of Lemma 7. Recall that by our assumption there are exactly k active constraints defined by w̃. Then
from Lemma 6 it follows that w̃ is a minimizer of D, i.e. VAL Dλ = VAL D. By the integrality of D we
have that VAL D ∈ Z. Let λ0 = λ − ε, for some 0 < ε < 1/d. Then, it holds VAL Dλ0 = VAL D as only
the k variables correponding to active constraints increased by ε while w0 = λ0 decreased by ε. From the
strong duality we also have VAL Pno-k (λ0 ) = VAL Dλ0 = VAL D.
Let ũ be an integral optimal solution of Pno-k (λ0 ). Following the definition we have
VAL Pno-k (λ0 ) = (cT − λ1T )ũ + λk + ε(1T ũ − k).
Now we have the following properties: ũ is integral; 0 < ε|1T ũ − k| < 1 whenever 1 ≤ |1T ũ − k| ≤ d;
c ∈ Zd ; v ∈ Z; and VAL Pno-k (λ0 ) ∈ Z. Therefore, we have 1T ũ − k = 0, and hence 1T ũ = k. Since we
showed the equivalence between Pno-k (λ0 ) and P ROJ L AGR(λ0 , c) the lemma follows. 2

7

Experiments

We empirically validate the claims outlined in the previous sections. To that end, we compare LASSP with
the O(dk)-time dynamic program (DP) described in Section C as a baseline. Note that this DP already
has a better time complexity than the best previously published algorithm from [14]. Both algorithms are
implemented in the Julia programming language (version 0.5.0), which typically achieves performance close
to C/C++ for combinatorial algorithms.

7.1

Synthetic data

We perform experiments with synthetic data in order to investigate how the algorithms scale as a function
of the input size. We study two different setups: (i) the running time of the projection algorithms on their
own, and (ii) the overall running time of a sparse recovery algorithm using the projection algorithms as a
subroutine. For the latter, we use the structure-aware variant of the popular CoSaMP algorithm [24, 3].
Figures 3(a)-(b). Given a problem size d, we set the sparsity to k = d/50 and generate a random separated
sparse vector with parameter ∆ = (d − 5(k + 1))/k − 1. The non-zero coefficients are i.i.d. ±1. For the
projection-only benchmark, we add Gaussian noise with σ = 1/10 to all coordinates in order to make the
problem non-trivial. For each problem size, we run 10 independent trials and report their mean.
11

100
50
200

150
100

LASSP
DP
no model

50
0

400

Size d (in 102 )

100

200

300

Running time (s ·102 )

150

Running time (s)

Speed-up of LASSP

(b)

(a)

(c)
102
101
100
10−1

Size d (in 102 )
(e)

(d)

LASSP
DP

0.2

0.2

0.2

0

0

0

−0.2

−0.2

−0.2

104.5 105
Size d
(f)

105.5

Figure 3: Separated sparsity experiments. In the top row we plot running times of our algorithm LASSP
relative to the previous work. Plots (a) and (c) are obtained by projecting signals on the separated sparsity
model. Plot (b) compares the running times of CoSaMP with three different projection operators. The variant
"no-model" makes no structural assumptions and uses hard thresholding as projection operator. In the bottom
row, in (e) we plot the recovery of the signal (d) obtained by the algorithm "no model". Plot (f) shows the
recovery of LASSP. The both procedures use the same number of measurements.
Figure 3(a) shows the speed-up obtained by our nearly-linear time projection relative to the DP baseline.
We observe that LASSP is up to 150× faster. This confirms our expectation that LASSP scales gracefully
with the problem size, while the DP essentially becomes a quadratic-time algorithm.
Figure 3(b) compares the running times of CoSaMP with three different projection operators. The first
variant makes no structural assumptions and uses hard thresholding as projection operator. The other two
variants use a projection for the separated sparsity model, relying on the DP baseline and LASSP, respectively.
The results show that the version of CoSaMP using LASSP instead of the DP is significantly faster. Moreover,
CoSaMP with a simple sparse projection has similar running time to CoSaMP with our structured projection.
Finally, we note that CoSaMP with separated sparsity requires 1.5× fewer measurements to achieve the same
recovery quality as CoSaMP with standard sparsity.
Figure 3(c). We fix the sparsity k = 100 and vary the length d of the signal. The signal is obtained in
the same way as for plots Figure 3(a)-(b). We observe that our algorithm runs 10x faster than the baseline.
Furthermore, the plot shows that the running times of both the baseline DP and LASSP scale linearly with the
signal length d. This behavior is expected for the DP. On the other hand, our theoretical findings predict that
the running time of LASSP scales as d log d. This suggests that the empirical performance of our algorithm
is even (slightly) better than what our proofs state.
In Appendix A we report results of additional experiments.

12

7.2

Neuronal signals

We also test our algorithm on neuron spike train data from [12]. See Figure 3(d) for this input data. First, we
run CoSaMP with a “standard sparsity” projection. The recovered signal is depicted in Figure 3(e). Next, we
run the convolutional sparsity CoSaMP of [12] and use our fast projection algorithm. The recovered signal is
given in Figure 3(f). For the both experiments we use n = 250 measurements.
We also run the convolutional sparsity CoSaMP on neuron spike train data of length 105 , comparing the
running time of LASSP and DP as projection operators. CoSaMP with our algorithm runs 2× faster in this
context. We do not compare the running time relative to CoSaMP with a standard sparsity projection as it
requires 10× more measurements to achieve accurate recovery.

8

Conclusions

We have designed a nearly-linear time algorithm for projecting onto the set of separated sparse vectors. The
core technique in our algorithm is Lagrangian relaxation. One of the key insights here is that even though there
are separated sparsity instances for which the Lagrangian relaxation does not provide an optimal solution, it is
still possible to obtain an optimal solution if the original input instance is only slightly perturbed.Furthermore,
this perturbation does not change the final output, but rather drives the algorithm to choose an optimal solution
of interest even after the hard sparsity constraint is relaxed. Our experiments show that our algorithm is not
only of theoretical significance, but also outperforms the state of the art in practice. Exploring the power of
perturbed Lagrangian relaxations for other non-convex constraint sets is an important direction for future
work. We believe that our framework will enable simple and efficient algorithms for other problems as well.
Acknowledgments. We thank Arturs Backurs for insightful discussions. L. Schmidt thanks Chinmay
Hegde for providing a dataset for some of our experiments. A. Madry
˛
was supported in part by an Alfred.
P. Sloan Research Fellowship, Google Research Award and the NSF grant CCF-1553428. S. Mitrović was
supported by Swiss NSF (grant number P1ELP2_161820). Part of this work was carried out while S. Mitrović
was visiting MIT.

13

References
[1] M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” Journal of
the Royal Statistical Society: Series B (Statistical Methodology), 2006.
[2] Y. Eldar and M. Mishali, “Robust recovery of signals from a structured union of subspaces,” IEEE
Transactions on Information Theory, vol. 55, no. 11, pp. 5302–5316, 2009.
[3] R. G. Baraniuk, V. Cevher, M. F. Duarte, and C. Hegde, “Model-based compressive sensing,” IEEE
Transactions on Information Theory, vol. 56, no. 4, pp. 1982–2001, 2010.
[4] J. Mairal, R. Jenatton, G. Obozinski, and F. Bach, “Convex and network flow optimization for structured
sparsity,” The Journal of Machine Learning Research, vol. 12, pp. 2681–2720, 2011.
[5] J. Huang, T. Zhang, and D. Metaxas, “Learning with structured sparsity,” The Journal of Machine
Learning Research, vol. 12, pp. 3371–3412, 2011.
[6] N. S. Rao, B. Recht, and R. D. Nowak, “Universal measurement bounds for structured sparse signal
recovery.” in AISTATS, ser. JMLR Proceedings, vol. 22, 2012, pp. 942–950.
[7] S. N. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu, “A unified framework for high-dimensional
analysis of M -estimators with decomposable regularizers,” Statistical Science, vol. 27, no. 4, pp. 538–
557, 11 2012.
[8] B. Bah, L. Baldassarre, and V. Cevher, “Model-based sketching and recovery with expanders,” in SODA,
2014, pp. 1529–1543.
[9] C. Hegde, P. Indyk, and L. Schmidt, “A nearly-linear time framework for graph-structured sparsity,” in
ICML. JMLR Workshop and Conference Proceedings, 2015, pp. 928–937.
[10] C. Hegde, M. F. Duarte, and V. Cevher, “Compressive sensing recovery of spike trains using a structured
sparsity model,” in SPARS’09-Signal Processing with Adaptive Sparse Structured Representations,
2009.
[11] E. L. Dyer, M. F. Duarte, D. H. Johnson, and R. G. Baraniuk, “Recovering spikes from noisy neuronal
calcium signals via structured sparse approximation,” in LVA/ICA, 2010, pp. 604–611.
[12] C. Hegde and R. G. Baraniuk, “Sampling and recovery of pulse streams,” IEEE Transactions on Signal
Processing, vol. 59, no. 4, pp. 1505–1517, 2011.
[13] E. L. Dyer, C. Studer, J. T. Robinson, and R. G. Baraniuk, “A robust and efficient method to recover
neural events from noisy and corrupted data,” in 6th International IEEE/EMBS Conference on Neural
Engineering (NER), 2013, pp. 593–596.
[14] S. Foucart, M. F. Minner, and T. Needham, “Sparse disjointed recovery from noninflating measurements,”
Applied and Computational Harmonic Analysis, vol. 39, no. 3, pp. 558 – 567, 2015.
[15] Q. Huang and S. M. Kakade, “Super-resolution off the grid,” in Conference on Neural Information
Processing Systems (NIPS), 2015, pp. 2665–2673.

14

[16] A. Aggarwal, B. Schieber, and T. Tokuyama, “Finding a minimum-weightk-link path in graphs with the
concave monge property and applications,” Discrete & Computational Geometry, vol. 12, no. 3, pp.
263–280, 1994.
[17] M. F. Duarte and R. G. Baraniuk, “Spectral compressive sensing,” Applied and Computational Harmonic
Analysis, vol. 35, no. 1, pp. 111 – 129, 2013.
[18] T. Needham, “Dictionary-sparse and disjointed recovery,” in International Conference on Sampling
Theory and Applications (SampTA), 2015, pp. 278–282.
[19] M. F. Duarte and Y. C. Eldar, “Structured compressed sensing: From theory to applications,” IEEE
Transactions on Signal Processing, vol. 59, no. 9, pp. 4053–4085, 2011.
[20] F. Bach, R. Jenatton, J. Mairal, and G. Obozinski, “Structured sparsity through convex optimization,”
Statistical Science, vol. 27, no. 4, pp. 450–468, 11 2012.
[21] M. J. Wainwright, “Structured regularizers for high-dimensional problems: Statistical and computational
issues,” Annual Review of Statistics and Its Application, vol. 1, no. 1, pp. 233–253, 2014.
[22] C. Hegde, P. Indyk, and L. Schmidt, “Fast algorithms for structured sparsity,” Bulletin of the European
Association for Theoretical Computer Science, vol. 117, p. 197–228, 2015.
[23] G. L. Nemhauser and L. A. Wolsey, Integer and Combinatorial Optimization.
1988.

Wiley-Interscience,

[24] D. Needell and J. A. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate
samples,” Applied and Computational Harmonic Analysis, vol. 26, no. 3, pp. 301–321, 2009.
[25] C. Hegde, P. Indyk, and L. Schmidt, “Approximation algorithms for model-based compressive sensing,”
IEEE Transactions on Information Theory, vol. 61, no. 9, 2015.
[26] ——, “Nearly linear-time model-based compressive sensing,” in ICALP, 2014, vol. 8572.
[27] R. J. Fowler, M. S. Paterson, and S. L. Tanimoto, “Optimal packing and covering in the plane are
np-complete,” Information processing letters, vol. 12, no. 3, pp. 133–137, 1981.
[28] A. Schrijver, Combinatorial optimization: polyhedra and efficiency.
Media, 2002, vol. 24.

15

Springer Science & Business

150
100
50
20 40 60 80 100

(b)
15
10
5
20

Size d (in 103 )

40

60

80

Running time (s ·102 )

200

Speed-up of LASSP

Speed-up of LASSP

(a)

(c)
102

LASSP
DP

101
100
101 102 103
Sparsity k

Size d (in 103 )

Figure 4: The plots compare the running times of our algorithm LASSP and the DP baseline. Plots (a) and
(c) are obtained by projecting signals directly on the separated sparsity model. Plot (b) shows the speed-up of
CoSaMP recovery obtained by using LASSP as a projection operator relative to using the DP for projection.

100

LASSP
DP
no model

50
0

50

100

150

Size d (in 102 )

(c)

200
150
100

LASSP
DP
no model

50
0

100

200

Size d (in 102 )

Running time (s)

150

(b)
Running time (s)

Running time (s)

(a)

200
150
100

LASSP
DP
no model

50
0

200

400

Size d (in 102 )

Figure 5: The plots compare the running times of CoSaMP with three different projection operators. The
variant "no-model" makes no structural assumptions and uses hard thresholding as projection operator. The
other two variants use a projection for the separated sparsity model, relying on the DP baseline and LASSP.

A

Additional experiments

In this section we present experiments additional to those present in Section 7.
In plots Figure 4(a)-(b) we run experiments on real signals. The signals are obtained from signal in
Figure 3(d) by concatenating its copies. Figure 4(a) is obtained by projecting the signal directly to the
separated sparsity model. Then we plot the ratio of the running times of the baseline DP and our algorithm
LASSP. Figure 4(b) represents the same type of speed-up ratio but this time for the CoSaMP recovery of
signal measurements.
For the plots Figure 4(c) and Figure 5(a)-(c) we make the setup similar to the one for Figures 3(a)-(b).
Namely, given a problem size d, parameters α and β we set the sparsity to k = d/α and generate a random
separated sparse vector with parameter ∆ = (d − β(k + 1))/k − 1. The non-zero coefficients are i.i.d. ±1.
We add Gaussian noise with σ to all coordinates in order to make the problem non-trivial. For each problem
size, we run 10 independent trials and report their mean.
In Figure 4(c) we compare the running times of projections of the baseline DP and our algorithm LASSP
on synthetic data for the following parameters: d = 5 · 104 ; varying the sparsity k; β = 5; and σ = 0.5. This
plot confirms our theoretical findings that the running time of LASSP scales more gracefully with the growth
of the sparsity than the running time of DP does.

16

Similar to Figure 3(b), Figure 5(a)-(c) compare the running times of CoSaMP with three different
projection operators. The first variant makes no structural assumptions and uses hard thresholding as
projection operator. The other two variants use a projection for the separated sparsity model, relying on the
DP baseline and LASSP, respectively. We use the following parameters:
• Figure 5(a): α = 10; β = 2; and σ = 0.01
• Figure 5(b): α = 20; β = 8; and σ = 0.95
• Figure 5(c): α = 100; β = 30; and σ = 0.2.
The results show that CoSaMP with a simple sparse projection has similar running time to CoSaMP with our
structured projection. Moreover, CoSaMP using LASSP as the structured projection is usually significantly
faster than CoSaMP using the baseline DP. We note that even in the favorable case for the DP structured
projection, like in Figure 5(c) where the sparsity k is 100× smaller than the dimension d, our algorithm again
shows better performance.

B

An approximate projection counterexample

Multiple recent algorithms for structured sparse projections build on the approximation-tolerant framework of
[25]. In this framework, it suffices to design approximate projections instead of solving the model projection
problem exactly. For some sparsity structures, this approach has led to significantly faster algorithms [26, 9].
Hence it is interesting to see whether the approximation route is also helpful for separated sparsity. In fact, it
is easy to design the following 2-approximation algorithm that runs in nearly-linear time.
Partition the vector c into blocks of length ∆. Then, number those blocks 1 through dd/∆e in the order
they appear in c. In each block, choose an index corresponding to the largest coordinate of c within that block.
Split those indices into two groups G and H based on the parity of the corresponding block. Formally, define
G and H as follows


G := arg
max
cj | 1 ≤ i ≤ dd/∆e and i is odd ,
j=(i−1)∆+1... min{i∆,d}

and similarly

H :=

arg

max
j=(i−1)∆+1... min{i∆,d}


cj | 1 ≤ i ≤ dd/∆e and i is even .

Let g1 , . . . , g|G| be an ordering of the elements of G so that cgi ≥ cgj whenever i ≤ j. Similarly, let
h1 , . . . , h|H| be an ordering of the elements of H such that cgi ≥ cgj whenever i ≤ j.
Next, define G0 := {gi | 1 ≤ i ≤ min{k, |G|}} and H 0 := {hi | 1 ≤ i ≤ min{k, |H|}}. Now, it is not
hard to see that the set attaining larger value among G0 and H 0 has a solution value that is at least OP T /2,
where OP T is the maximum sum attainable with a ∆-separated and k-sparse vector (in the language of [25],
this is an approximate head projection).
While the above algorithm runs in nearly-linear time and achieves a constant-factor approximation,
there is a catch. In particular, the returned support pattern might contain fewer than k indices (as a simple
example, consider the case k = 2, ∆ = 2, and c = (1, 100, 1)). This raises the question of the sample
complexity for this relaxed sparsity model. Let k 0 denote the number of output indices of the 2-approximation
algorithm. As we have described in the preliminaries, the sample complexity of the separated sparsity model
is O(k log (n/k − ∆)). For concreteness, we now consider the case k 0 = k/2 and d/k − ∆ ∈ O(1) (the
17

latter is the important regime where the separated sparsity model achieves a sample complexity of m = O(k)).
Then we have d/k 0 − ∆ ∈ O(1) + ∆. Therefore, if we would like to apply the 2-approximation algorithm in
the case k 0 = k/2, the sample complexity would be
O(k log (2d/k − ∆)) = O(k log ∆) = O(k log d/k)
where we used the assumption that d/k − ∆ ∈ O(1). It is important to note that this sample complexity is
worse than O(k) and falls back to the O(k log d/k) sample complexity of “standard” k-sparse recovery. So
by using an approximate projection we have lost the sample complexity advantage of the separated sparsity
model. Since we already know how to recover k-sparse vectors in nearly-linear time without resorting to
structured sparsity, the approximation algorithm does not provide a novel trade-off.

C

Dynamic programming

Before we introduce our dynamic programming (DP) algorithm for the separated sparsity problem, we briefly
review a variant of the DP given in [14]. We remark that the variant below has a time complexity of O(dk),
which is already an improvement over the O(d2 ) of [14] when k = o(d). However, the authors are not aware
of a work the describes the O(dk) DP approach, and hence refer to it as folklore. We improve the analysis of
this DP further and give a faster variant that runs in time O(k(d − (k − 1)∆)). In the regime where where√the
e d),
slack d/k − ∆ is constant, the running time simplifies to O(d + k 2 ). So for the very sparse case k ∈ O(
our improved DP already achieves a nearly-linear running time.

C.1

The basic dynamic program

The folklore dynamic program fills the following table DP [i][j]. Value DP [i][j] can be interpreted as
follows: Maximal value of choosing j coordinates of c with positions in 1 . . . i, such that any two chosen
coordinates are at distance ∆ at least.
DP [i][0] = 0
DP [i][j] = −∞,

if i < 1

DP [i][j] = max {DP [i − 1][j], ci + DP [i − ∆][j − 1]} , otherwise
The following claim follows easily:
Theorem 3 (Folklore) DP [n][k] is an optimal value to the ∆-separated problem and can be computed in
time O(dk). The elements that constitute the optimal value can be output from DP in time O(d + k).

C.2

An improved dynamic program

The dynamic program outlined above already achieves an improvement over the algorithms in [10, 14].
However, the dynamic program is still wasteful with the state space it considers. Consider case in which there
is only one possible configuration, i.e. d = (k − 1)∆ + 1 for some k, ∆ ≥ 1. Regardless of the cost vector c,
an input with such parameters has only one valid solution. Nevertheless, the natural implementation of the
dynamic program above runs in Θ(dk) time. Note that the definition of DP does not take into account the
mandatory distance required by the separated sparsity model. By mandatory we refer to ∆ distance between

18

any two chosen coordinates of c. This observation gives raise to new dynamic programming definition, that
we denote by DP , which directly implements mandatory distances. Let
s := d − ((k − 1)∆ + 1),
which can be seen as the slack distance that is not mandatory. For instance, if we choose two coordinates of c
with indices i and i + ∆ + 5, and no other coordinate between them, then we say that ∆ distance between
them is mandatory, and the remaining distance of 5 is slack. Now, DP is defined as
DP [i][0] = 0
DP [i][j] = −∞,
if i < 0

DP [i][j] = max DP [i − 1][j], cpos(i,j) + DP [i][j − 1] , otherwise
where pos(i, j) = d − (s − i) − (k − j)∆. Value pos(i, j) gives us the smallest index of c that we can chose
if to the left of that index there are j other chosen indices, and slack distance of i is used. With the definition
of DP in hands one can easily show the following theorem.
Theorem 4 Value DP [s][k] is an optimal value to the ∆-separated problem and can be computed in time
O(ks). The elements that constitute the optimal value can be output from DP in time O(d + k).
Comparing Theorem 3 and Theorem 4, we see that for certain input parameters, e.g. when s ∈ O(d), we did
not make much progress by devising DP . However, as we have already mentioned, the setting we care the
most is when d/k − ∆ is a constant, as with such a choice of parameters the sample complexity is O(k). But
in that case we do achieve a significant improvement over Theorem 3.
Corollary 4 If d/k −∆ ∈ O(1), then there is an algorithm that solves the ∆-separated problem in O(k 2 +d)
time.

Two dimensional ∆-separated problem is NP-Hard

D

In this section we consider a natural extension of the ∆-separated problem in which vector c is two dimensional, i.e. c ∈ Rd×m . Formally, given c, sparsity parameter k, and integer ∆, the goal is to output k pairs of
integers (it , j t ) such that:
• 1 ≤ it ≤ d, and 1 ≤ j t ≤ m, for all t = 1 . . . k;

• for every t 6= s we have min it − is , j t − j s
Pk
•
t=1 cit ,j t is maximized.

≥ ∆; and

We refer by 2D∆-separated(c, k, ∆) to that problem, and show it is NP-hard by reducing it to problem
Box-Pack studied in [27]. Let us start by recalling the definition of Box-Pack.
Let k be a number of identical 3 × 3 squares.1 Our goal is to pack the k squares into a region of plane
defined by set R. Set R consists of pairs of integers, where every pair represents the point at which a square
can be placed, e.g. the upper-left corner of a square. Every square can only be placed so that its sides are
1

We choose sides to be of length 3 as already that setting is sufficient to prove NP-hardness of Box-Pack, see the proof of
Theorem 2 in [27].

19

parallel to the axis. By Box-Pack(k, R) we refer to the problem of answering whether for given k and R
one can place all the k squares in the described way, such that no two squares overlap. In [27] is proved the
following result.
Theorem 5 (Theorem 2 in [27]) Box-Pack is N P -complete.
Now we utilize Theorem 5 to show that 2D∆-separated(c, k, ∆) is NP-hard as well.
Theorem 6 Problem 2D∆-separated(c, k, ∆) is NP-hard.
Proof. We provide a polynomial time reduction of Box-Pack(k, R) to a sequence of 2D∆-separated instances.
Let RE , denoting "extended R", be the set of all integers points that can be occupied by a square, not
necessarily its corners. Observe that RE ≤ 9|R|. We say that two points in RE are adjacent if they share
the same x- or the same y-coordinate. Consider a connected component F of RE .
Assume that Box-Pack(i, F ) can be reduced to an 2D∆-separated instance. Then, let us show that
Box-Pack(k, R) can be reduced to polynomially many 2D∆-separated instances.
Let F be the set of all connected components of RE . Then, for every C ∈ F define valF as follows
valF := arg

max

i∈{0,...,min{|F |,k}}

Box-Pack(i, F ) equals true.

P
Now, clearly, if F ∈F valF ≥ k, then Box-Pack(k, R) equals true. Observe that in order to compute all
valF , we need to compute only polynomially many instances of Box-Pack, where all the input parameters
are bounded by k and R. So, it remains to reduce Box-Pack(i, F ) to an instance of 2D∆-separated.
To that end, consider F ∈ F. Let h and w be the smallest values such that F is enclosed by an axisparallel sides rectangle T of height h and width w. With loss of generality, assume that T and F are translated
so that the bottom-left corner of T at (1, 1). Define vector c ∈ {0, 1}h×w to be an indicator vector of points
of C at which can be placed square, i.e.
(
1 if 3 × 3 square with upper-left corner at (i, j) is in F
F
ci,j =
0 otherwise
Note that hw ≤ |F |2 . Next, define optF as follows
optF :=

max

i∈{1,...,min{|F |,k}}

2D∆-separated(cF , i, 3).

Then, we claim optF equals valF . Now, it is easy to see that the claim is true. Value valF represents the
maximal number of squares that can be packed within F . Each such square corresponds to an entry of cF
which has value 1. So, we have 2D∆-separated(cF , valF , 3) = valF , and hence optF ≥ valF . On the other
hand, by the construction of cF , optF represents a number of squares (not necessarily maximum, though)
that can be packed in F . So, we have valF ≥ optF , and therefore optF = valF as claimed.
This concludes our proof. 2

20

E

Uniform size neuronal spike trains

We propose a generalization of the ∆-separated model. In this model, we assume that neuronal spike trains
correspond to blocks of uniform size. We use b to denote the size of blocks. Formally, we define model
Mk,∆,b to be the set of all vector in Rd such that non-zero entries are grouped in blocks of size b, there are
exactly k blocks, and every two blocks are separated by at least ∆.
In order to apply the existing framework for recovering measured vectors that belong to Mk,∆,b , we
develop a method that solves the following projection problem. Given vector c ∈ Rd , the goal is to output k
indices p1 , . . . , pk such that
• 1 ≤ pi ≤ d − b + 1, for every i = 1 . . . k;
• |pi − pj | ≥ ∆ + b − 1, for all i 6= j; and,
• the sum

k X
b−1
X

cpi +j

i=1 j=0

is maximized.
We refer to that problem by (∆, b)-separated. Now, it is easy to show the following claim.
Theorem 7 Given vector c ∈ Nd+ , sparsity k, and block size b, define vector cb ∈ Nd−b+1
as
+
cbi

:=

i+b−1
X

cj .

j=i

Let S be a set of indices corresponding to an optimal solution for (∆ + b − 1)-separated problem for cost
vector cb and sparsity k. Then, S is an optimal solution to (∆, b)-separated problem for cost vector c and
sparsity k.
Proof. Clearly, S is feasible choice of indices for (∆, b)-separated problem. Towards a contradiction, assume
that there exists another set of indices S 0 that achieves larger value than S.
But then, S 0 is a feasible choice of indices for (∆ + b − 1)-separated problem. Next, notice that S,
and also S 0 , achieve the same value for both (∆, b)- and (∆ + b − 1)-separated problem. However, as S 0
achieves higher value than S, this contradicts our assumption that S is an optimal choice of indices for
(∆ + b − 1)-separated problem. 2
Now we can solve (∆, b)-separated problem in nearly-linear time.
Corollary 5 A solution to (∆, b)-separated problem can be computed in nearly-linear time.
Proof. We use the reduction from Theorem 7 to reduce (∆, b)-separated problem to (∆ + b − 1)-separated
one. The reduction can be applied in linear time. Then, the claim follows by Theorem 12. 2

21

E.1

Sample complexity

Next we analyze the sample complexity of Mk,∆,b model. To that end, we count the number of support
P
in that model. Now, all the supports are captured by (α0 , α1 , . . . , αk ), where αi ≥ 0 and ki=0 αi =
d − kb − (k − 1)∆. The idea is that k blocks split the vector in k + 1 regions. All but the first and last region
correspond to coordinates between two neighboring blocks. Value α0 and αk correspond to coordinates
before the first block and after the last block, respectively. Value αi , for 1 ≤ i < k correspond to the slack
distance between block i and i + 1.
Let Mk,∆,b be the family of support patterns of Mk,∆,b . Then, we have:


d − kb − (k − 1)∆ + k
|Mk,∆,b | =
.
k
Now, from the result that there are RIP matrices with O (k + log |Mk,∆,b |) rows, [3], the following claim
follows.
Theorem 8 The sample complexity of Mk,∆,b is O(k log (d/k − (b + ∆) + (∆ − 1)/k)). In particular, if
(∆ − 1)/k ∈ O(1), then the sample complexity is O(k log (d/k − (b + ∆)).

F

Omitted proofs from Section 5.1 and study of D

Our running-time results are given with respect to γ, where γ is the maximal number of bits needed to
store any ci . As P and D are invariant under shifting and multiplication of c, for the sake of clarity of our
exposition and without loss of generality we assume c is an integral vector with non-negative entries.
Lemma 8 The constraint matrix of problem P is totally unimodular (TUM).
Proof. We first rewrite the constraints of P to be in the form Bu ≤ b, u ≥ 0, where Bu ≤ b is defined as
d
X

ui ≤ k

i=1
d
X

−ui ≤ −k

i=1
min{i+∆−1,d}

X

uj ≤ 1

∀i = 1 . . . d

j=i

Let D be a square submatrix of B. If we show that det D ∈ {−1, 0, 1}, then by the definition of TUM the
lemma follows.
Now, D is either a binary interval matrix, or it has one row where all the non-zero entries are consecutive
0 :=
and equal −1 and the other rows constitute a binary interval matrix. Let D0 be a matrix defined as Di,j
0
0
|Di,j |. Then, we have |det D| = |det D |. Also, we have that D is a binary interval matrix. However, it is
well known that binary interval matrices are TUM, see [23]. Hence, det D ∈ {−1, 0, 1} 2
Corollary 2 For an integer k and a vector of integers c, there exists ŵ such that ŵ is an optimum of D and
ŵ0 ∈ Z.
22

Proof. By Lemma 8, problem P is totally unimodular. By [28] we have that D is totally unimodular as well.
This implies that D has an integral optimal solution. 2
The crucial property of the LP D is the following: for a fixed value of the variable w0 , we can solve D in
linear time. Specifically, let us define Dv to be the LP D in which the variable w0 is set to v. (Note that by
Corollary 2, we can restrict our attention to integer values of v.) Observe that once we fixed the value of w0 ,
all remaining constraints in Dv are “local” since they only affect a known interval of length ∆. They are also
ordered in a natural way. As a result, we can solve Dv by making a single pass over these variables. Starting
with w1 and all variables set to 0, we consider each constraint from left to right and increase the variables to
satisfy these constraints in a lazy manner. That is, if reach in our pass a constraint with index i that is still
not satisfied, we increase the value of wi until that constraint becomes satisfied and then move to the next
constraint and index. Algorithm 2, called Dual-Greedy, present in Section 5 formalizes this approach.
Lemma 9 For any c0 , and k 0 , the algorithm Dual-Greedy(c0 , v) computes an optimal solution Dv (c0 , k 0 ) in
linear time.
0

Proof. Let c0 ∈ Nd+ . We show that the algorithm Dual-Greedy(c0 , v) outputs vector w̃ in O(d0 ) time such that
Dv (c0 , k 0 ) is minimized and w̃0 = v. To show the running time, observe that every iteration in the for loop
takes O(1) time, so the total algorithm runs in O(d0 ) time.
Next, it is easy to see that w̃ is a feasible solution to Dv (c0 , k 0 ). We show that it is a minimal as well.
Towards a contradiction, assume there is a vector w0 such that w00 = w̃0 = v, w0 is a feasible solution to
D(c0 , k 0 ), and kw0 k1 < kw̃k1 . Then, there exists an index j such that wj0 6= w̃j and wi0 = w̃i for all i < j. In
case there are multiple vectors w0 , let w0 be one that maximizes the location of mismatch j. Consider the two
possible cases: wj0 < w̃j , and wj0 > w̃j .
Case wj0 < w̃j . Observe that w̃j is chosen as a function of w̃0 , . . . , w̃j−1 as a minimal value so that w̃ is
feasible. Therefore, if wj0 < w̃j and wi0 = w̃i for all i < j, then w0 could not be feasible.
Case wj0 > w̃j . First, note that j < d0 , as otherwise kw0 k1 > kw̃k1 . Now, we construct w00 as follows. We
00
0
set wi00 = wi0 for all i different than j and j + 1. Set wj00 = w̃j and wj+1
= wj+1
+ w̃j − wj0 . Clearly, w00 is
0
0
00
0
00
also a feasible solution to D(c , k ). Furthermore, kw k1 = kw k1 and w̃ and w agrees on first j coordinates,
contradicting our choice of w0 .
This concludes the proof. 2
Now, to obtain an algorithm that is also able to solve the original dual LP D (instead of only Dv ) it
suffices to provide a procedure for choosing the optimal value of v. A priori, there can be many possible
choices of v and thus an exhaustive search would be prohibitive. Fortunately, Dv is actually convex in v, and
so we can use a ternary search over v to find such an optimal value.
Lemma 10 Dv (c0 , k 0 ) is convex in v.
Proof. Let w̃ and ŵ be w-solution to Dv̂ (c0 , k 0 ) and Dṽ (c0 , k 0 ), respectively. Note that ŵ0 = v̂ and w̃0 = ṽ.
Now if we show that
 Dv̂ (c0 , k 0 ) + Dṽ (c0 , k 0 )
D v̂+ṽ c0 , k 0 ≤
,
2
2
the convexity will follow. We start by showing that ŵ+2 w̃ is a feasible w-vector for the dual, assuming that
both ŵ and w̃ are feasible. We consider the feasibility of each of the constraints.
23

(1) From
X

ŵ0 +

ŵj ≥ c0i

j : j≤i≤j+∆−1 and j≥1

and
X

w̃0 +

w̃j ≥ c0i

j : j≤i≤j+∆−1 and j≥1

we have
X

(ŵ0 + w̃0 ) +

(ŵj + w̃j ) ≥ 2c0i ,

j : j≤i≤j+∆−1 and j≥1

which implies
ŵ0 + w̃0
+
2

X
j : j≤i≤j+∆−1 and j≥1

(2) Also, from ŵj ≥ 0 and w̃j ≥ 0 we have
(3) Trivially,
So, indeed

ŵ0 +w̃0
2

ŵ+w̃
2

ŵj +w̃j
2

ŵj + w̃j
≥ c0i .
2

≥ 0.

∈ R.

is feasible. Hence, we have
d0

D v̂+ṽ

≤

2

=
=

ŵ0 + w̃0 0 X ŵj + w̃j
k +
2
2
i=1




Pd0
P
d0
0
0
ŵ0 k + i=1 ŵj + w̃0 k + i=1 w̃j
2
Dv̂ (c0 , k 0 ) + Dṽ (c0 , k 0 )
.
2

This completes the proof. 2
Putting these pieces together yields the main theorem of this section. As a reminder, we assume that c is
an integral vector with non-negative entries.
Theorem 9 There exists an algorithm Opt-Value-of-D(c0 , k 0 ) that, for any c0 and k 0 , outputs an optimal
solution w? to D(c0 , k 0 ) in time O (d0 (log c0max + log k 0 )).
Proof.
Bounding w0? . Following Lemma 10, we can use ternary search to find w0? . However, to be able to do that,
we have to know the interval [a, b] we are searching over for w0? . By the following lemma we give an upper
bound on b. However, in order to provide a lower bound on a, we need to develop some more machinery. So,
we defer its proof to later sections, and in Lemma 18 we show that a can be lower bounded by −(k − 1)cmax .
Lemma 11 Let w? be a vector that minimizes D. Then, w0? ≤ cmax .

24

Proof. Towards a contradiction, let w? be an optimal vector such that w0? > cmax . Now, as wi ≥ 0 for
all i ≥ 1, we have that the corresponding objective is at least kw0? > kcmax . On the other hand, consider
vector ŵ such that ŵ0 = cmax and ŵi = 0 for all i ≥ 1. Clearly, ŵ is a feasible solution to D. However,
the objective function corresponding to ŵ is kcmax < kw0? , which contradicts our assumption that w? is a
minimizer of D. 2

A nearly-linear time algorithm. Now we provide an algorithm that computes the optimal value of D(c0 , k 0 )
in nearly linear time. It employs ternary search over the interval provided by Lemma 18 and Lemma 11 in
order to find w0 that optimizes D(c0 , k 0 ). At every step of the search, it uses the result from Lemma 9 to find
an optimal solution to Dv (c0 , k 0 ), for v chosen at the current search step.
Algorithm 3 Opt-Value-of-D:
0

Input: c0 ∈ Nd+ , sparsity k 0 , lower bound lb (if not specified, the default value is −(k 0 − 1)c0max ), upper
bound ub (if not specified, the default value is c0max )
Output: a minimizer wbest to D(c0 , k 0 ) constrained to w0best ∈ [lb, ub]; if lb and ub are not specified, wbest is
an optimal solution to D(c0 , k 0 )
1. s ← lb, e ← ub
2. While s ≤ e


 e−s 
1. l ← s + e−s
,
r
←
e
−
3
3
2. (activel , wl ) ← Dual-Greedy(c, l)
3. (activer , wr ) ← Dual-Greedy(c, r)
P0
P0
4. If kw0l + di=1 wil ≤ kw0r + di=1 wir then e ← r − 1, wbest ← wl
5. Else s ← l + 1, wbest ← wr
3. Return wbest
Lemma 9, Lemma 10, Corollary 2 and bounds on a and b, the theorem follows directly. 2

G

A detailed proof of Lemma 5

Algorithm ∆-Separated-no-Sparsity solves Pno-k (λ) (and so P ROJ L AGR).

25

Algorithm 4 ∆-Separated-no-Sparsity:
Input: c ∈ Rd , λ ∈ R
Output: an integral optimal solution of Pno-k (λ)
1: c̃ ← c − 1λ
2: if c̃1 ≥ 0 then
3:
best1 ← c̃1 , pick1 ← 1
4: else
5:
best1 ← 0, pick1 ← −1
6: for i ← 2 . . . ∆ do
7:
if c̃i > besti−1 then
8:
besti ← c̃i , picki ← i
9:
else
10:
besti ← besti−1 , picki ← picki−1
11: for i ← ∆ + 1 . . . d do
12:
if c̃i + besti−∆ > besti−1 then
13:
besti ← c̃i + besti−∆ , picki ← i
14:
else
15:
besti ← besti−1 , picki ← picki−1
16: r ← ∅, i ← d
17: while i ≥ 1 do
18:
if picki ≥ 1 then r ← r ∪ {picki }
19:
i ← picki − ∆
20: return r
Lemma 12 Algorithm ∆-Separated-no-Sparsity(c, λ) solves Pno-k (λ) (and so P ROJ L AGR(λ, c) in O(d)
time.
Proof. Let us rewrite Pno-k (λ) as follows
maximize

(cT − λ1T )u + λk
min{i+∆−1,d}

subject to

X

uj ≤ 1

∀i = 1 . . . d

ui ≥ 0

∀i = 1 . . . d

j=i

Observe that for a given λ and k, the term λk is constant, so our objective becomes (cT − λ1T )u = c̃T u
subject to the given set of constraints. Now we can show that ∆-Separated-no-Sparsity(c, λ) solves this
formulation, which is equivalent to the original one of Pno-k (λ). We proceed by induction. We show that
besti stores the maximum sum of the elements of {c̃1 , . . . , c̃i } so that any two chosen elements c̃i and c̃j are
such that |i − j| ≥ ∆, and picki keeps what is the largest index of element that should be taken to achieve
besti (if there is no such element, then picki equals −1).
Base of induction, i = 1. For i = 1, the lines 2-5 set best1 and pick1 properly.
Inductive step, i > 1. For i = 2 . . . ∆, the loop at line 6 sets besti and picki as required. That is, either
c̃i is the largest among {c̃1 , . . . , c̃i }, or the largest element is among {c̃1 , . . . , c̃i−1 } which is properly set in
26

besti−1 and picki−1 . If there is no element with positive value, then as besti = 0 and picki = −1 as desired.
Similarly, for i > ∆ either we choose c̃i and obtain the remaining of the output over {c̃1 , . . . , c̃i−∆ }, or
we take the best solution over {c̃1 , . . . , c̃i−1 } without including c̃i .
Correctness of the output vector r. Since besti and picki are set as described above, the way r is obtained
trivially satisfies the constraints of Pno-k (λ). Furthermore, as bestd maximizes the objective of Pno-k (λ), r
maximizes the objective of Pno-k (λ) as well.
Running time. Line 1 takes O(d) time. The total number of iterations of loop at line 6 and loop at line 11
is O(d). Every iterations takes O(1) time. The loop at line 17 starts with i being d, and decreases i at every
iteration (under natural assumption ∆ ≥ 1). The fact that i gets decreases at every iteration comes from the
property that picki ≤ i, for every i. Therefore, the while loop takes O(d) iterations as well, while every
iteration taking O(1) time. This completes the analysis. 2

H

Omitted proofs from Section 6.2

In this section we prove some properties of active constraints. We begin by introducing some notation. In
what follows, we will be interested in cost vectors and sparsity parameters other than c and k, respectively.
Hence, whenever this is the case, we will denote the corresponding dual LP by D(c0 , k 0 ).
We next define active constraints algorithmicaly.
Algorithm 5 Active-Constraints:
0

0

Input: c0 ∈ Nd+ , w ∈ Rd +1 a feasible vector of D(c0 , k 0 )
Output: active constraints of D(c0 , k 0 ) for w
1: active ← ∅;
sum ∆ ← 0; last_active ← −∞
0
2: for i := 1 . . . d do
3:
if i − ∆ ≥ 1 then sum ∆ ← sum ∆ − wi−∆
4:
sum ∆ ← sum ∆ + wi
5:
if last_active ≤ i − ∆ and c0i − (w0 + sum ∆ ) = 0 then
6:
active ← active ∪ {i}
7:
last_active ← i
8: return active
Active constraints provide insight into the structure of D(c0 , k 0 ) that we can leverage to optimally
distribute our sparsity budget over the recursive subproblems. The optimality condition of the dual program
can be described in the language of active constraints as follows.
Lemma 13 Let w1 ← Dual-Greedy(c0 , v), active1 ← Active-Constraints(c0 , w1 ), w2 ← Dual-Greedy(c0 , v+
1), and active2 ← Active-Constraints(c0 , w2 ), for some integer v. Then, w1 is an optimal solution of D(c0 , k 0 )
iff |active1 | ≥ k 0 ≥ |active2 |.
Proof. We first prove two properties of actives constraints. First, observe that from the way algorithm
Active-Constraints(c0 , w) outputs active, it corresponds to tight constraints of D(c0 , k 0 ) for a given w. More
precisely, every such tight constraint either is in active, or there is another tight constraint in active which is
27

at most ∆ "far to the left". Furthermore, as c0 is integral, it is easy to see that tight constraints for w0 = v, for
some integer v, and for w0 = v − δ, for δ ∈ [0, 1), are the same. Hence, active constraints for w0 = v and
w0 = v − δ are also the same. Putting these observations together, we get the following claim.
Lemma 14 Let v be an integer and δ ∈ [0, 1). Also, let ŵ ← Dual-Greedy(c0 , v) and ŵ0 ← Dual-Greedy(c0 , v−
δ). Then,
Active-Constraints(c0 , ŵ) = Active-Constraints(c0 , ŵ0 ).
We point out that one can show even stronger statement about tight constraints, not necessarily active
tough. Namely, it holds that if a constraint i is tight for ŵ being Dual-Greedy(c0 , v), then it is tight for any
Dual-Greedy(c0 , v 0 ) such that v 0 ≤ v. It follows from the property that for any value v 0 there is at most one
active constraint j in Active-Constraints(c0 , Dual-Greedy(c0 , v 0 )) such that 0 ≤ i − j ≤ ∆. Therefore, if
w0 = v 0 gets decreased by "a very small" δ, then the variable, e.g. wj , corresponding to active constraint will
decrease by δ as well. Which in turn results i still being a tight constraint. Hence, the following lemma holds,
which we utilize in the sequel.
Lemma 15 Let ŵ ← Dual-Greedy(c0 , v) and w̃ ← Dual-Greedy(c0 , v 0 ), for v 0 ≤ v. Then, if a constraint i
is tight with respect to ŵ, it is tight with respect to w̃ as well.
Using Lemma 14 we can show how Dual-Greedy(c0 , v 0 ) changes for v 0 ∈ (v − 1, v].
Lemma 16 Let v be an integer and δ ∈ [0, 1). By ŵ denote the output of Dual-Greedy(c0 , v), and by ŵδ the
output of Dual-Greedy(c0 , v − δ). Let active be returned by Active-Constraints(c0 , ŵ). Then
X
X
k 0 ŵ0δ +
ŵiδ = k 0 ŵ0 +
ŵi + δ(|active| − k 0 ).
i≥1

i≥1

Proof. From Lemma 14 we have that Active-Constraints(c0 , Dual-Greedy(c0 , v − δ)) = active for all
δ ∈ [0, 1). Now, by the construction of active, it holds
X
X
X
k 0 ŵ0δ +
ŵiδ = k 0 (ŵ0 − δ) +
(ŵi + δ) +
ŵi
i∈active

i≥1
0

= k ŵ0 +

X

i≥1 and i∈active
/
0

ŵi + δ(|active| − k ),

i≥1

as desired. 2
We are now ready to finalize the proof of the lemma. Let us break the equivalence stated in the lemma
into two implications, and show they are true.
2 | be true. By Lemma 16 and the choice of v, we have that D (c0 , k 0 )
( ⇐= ) Let |active1 | ≥ k 0 ≥ |active
z

1
is non-decreasing for z ∈ v, v + 2 and non-increasing for z ∈ [v − 12 , v]. As Dz (c0 , k 0 ) is convex, we have
Dz (c0 , k 0 ) is minimized for z := v.2
2

We use the fact that from the convexity of Dz (c0 , k0 ) we have that Dz (c0 , k0 ) is continuous.

28

( =⇒ ) Let w1 is an optimal solution of D(c0 , k 0 ). Recall that Dz (c0 , k 0 ) is a convex function in z. Then, as
w1 is an optimum of D(c0 , k 0 ), Dx (c0 , k 0 ) is non-increasing in z on interval (−∞, v] and non-decreasing on
[v, ∞). But then from Lemma 16 we conclude that it can only happen if |active1 | ≥ k 0 ≥ |active2 |.
2

I

Omitted proofs from Section 6.3

In this section we finalize the proof of the correctness of our randomized algorithm. Before we delve into
details, we introduce some notation. In what follows, we will be interested in cost vectors and sparsity
parameters other than c and k, respectively. Hence, whenever this is the case, we will write P(c0 , k 0 ) to refer
to the program P for cost vector c0 and sparsity k 0 . Similarly, whenever we consider some different cost
vector c0 and sparsity parameter k 0 , we will denote the corresponding dual LP by D(c0 , k 0 ).
Also, as pointed out in other sections, our running-time results are given with respect to γ, where γ is the
maximal number of bits needed to store any ci . As P and D are invariant under shifting and multiplication of
c, for the sake of clarity of our exposition and without loss of generality we assume c is an integral vector
with non-negative entries.
Lemma 3 Algorithm LASSP runs only a single iteration with probability at least 1 − 1/d.

I.1

A proof of Lemma 3

As pointed out already, without loss of generality in this proof we assume that c is an integral non-negative
vector. We also recall that we showed equivalence between problem (1) and P, and also between Pno-k and
P ROJ L AGR, so in this proof we work with the LP formulations.
The choice of λ is optimal. Consider w̃ as in Lemma 7. First we want to show that if such w̃ exists, then λ
obtained at line 4 of Algorithm 1 is such that it also defines k active constraints. This in turn would imply, by
Lemma 7, that support Ŝ obtained at line 5 has cardinality k, and hence is an optimal solution to P(c̃).
If w̃0 = λ, then we are done. Otherwise, assume that w̃0 6= λ. By Lemma 13, w̃ is an optimal solution
of D(c̃). On the other hand, as we discussed in Section 5, λ is such that VAL Dλ (c̃) = VAL D(c̃), and
w̃0 < λ by the choice of λ. Therefore, by Lemma 13 it holds that Dual-Greedy(c̃, λ) defines at least k active
constraints. Furthermore, as D(c̃) is convex w.r.t. to the variable w0 , then for every λ0 ∈ [w̃0 , λ] we have
VAL Dλ0 (c̃) = VAL D(c̃). In other words, VAL Dλ0 (c̃) remains constant over the given interval. Hence,
from Lemma 16 we conclude that Dual-Greedy(c̃, λ) defines exactly k active constraints.
However, for c̃ given on the input there might not exists any λ such that Dual-Greedy(c̃, λ) defines exactly
k active constraints. Our goal is to show that the randomization we apply assures that for the obtained c̃ it is
always the case that there is some λ so that Dual-Greedy(c̃, λ) has exactly k active constraints.
The evolution of active constraints. Now we want to show that the randomization we apply will result in
an existence of w̃ as described in Lemma 7. We start by studying the evolution of active constraints defined
by the output of Dual-Greedy(c̃, w0 ) as w0 decreases.
First, recall that by Lemma 15 we have that if a constraint becomes tight with respect to some ŵ ←
Dual-Greedy(c̃, λ), it remains tight with respect to ŵ0 ← Dual-Greedy(c̃, λ0 ) for every λ0 ≤ λ. Let T denote
the set of tight constraints with respect to ŵ, and T 0 with respect to ŵ0 . By our discussion T ⊆ T 0 .
29

Let A and A0 be the set of active constraints with respect to Dual-Greedy(c̃, λ) and Dual-Greedy(c̃, λ0 ),
respectively. Clearly A ⊆ T and A0 ⊆ T 0 . We claim that |T ∩ A0 | ≤ |A|. Observe that A is a minimum set
of constraints so that every tight constraint is covered (covered in the natural way). However, A is also a
maximum set of constraints of T that can be chosen so that no two of them overlap, i.e. so that every two
of them are at least ∆ apart. That means if one would choose a subset of T larger than |A| then some two
constraints would overlap. Hence, such a subset can not consist of only active constraints, and therefore
|T ∩ A0 | ≤ |A|.
This implies that if the number of active constraints increases by a > 0 at certain point, then there are at
least some a constraints that became active, but also tight, for the first time. Now we want to study what is the
probability that two or more non-tight constraints become tight with respect to Dual-Greedy(c̃, λ), for any λ.
Let us focus on a single constraint j. Fix randomness of all the Xi for i 6= j, i.e. fix c̃i for all i 6= j. For
Xj = 0, there are at most d − 1 different values of w0 when any of those constraints becomes tight for the
first time. Let W denote the set of these w0 values. So |W | < d. Now, construct set Wj as follows. For each
λ ∈ W:
• If constraint j is already tight with respect to Dual-Greedy(c̃, λ) for Xj = 0, do nothing.
• Define ĉi = c̃i for i 6= j, and ĉj = c̃j + xλ , where xλ is defined as the least value so that constraint
j becomes tight for the first time with respect to Dual-Greedy(ĉ, λ). Observe that xλ ≥ 0 and xλ is
integral. If xλ < d3 , add xλ to Wj .
We have |Wj | ≤ |W | < d. Each of the value of Wj correspond to some value of Xj . Also, observe that for
given λ, a constraints can become tight for the first time for at most one value of Xj . In addition, as long as
constraint j is not tight it does not affect when the other constraints will become tight, as wj = 0 so wj has
no affect on other constraints. This all implies that there are at most d − 1 distinct values of Xj , out of d3 of
them, when constraint j and some other constraint become tight. Therefore,
P r [constraint j and any other constraint become tight for the same value of w0 = λ] <

d
1
= 2.
3
d
d

Now we can apply union bound to conclude
P r [two constraints become tight for the same value of w0 = λ] < d

1
1
= .
d2
d

Therefore, after applying randomness, for every value of w0 = λ at most one constraint becomes tight with
probability at least 1 − 1/d. Following our discussion above, this in turn implies that the number of active
constraints increases by at most 1 after decreasing the value of w0 by 1. Therefore, there is w̃ as in Lemma 7
with probability 1 − 1/d at least.
Required randomness. Every ci we perturb by one out of d3 different values, for which O(log d) random
bits suffices. Therefore, in total we need O(d log d) random bits.
This concludes the proof.

J

Deterministic worst-case nearly-linear time algorithm

In this section we describe our deterministic algorithm. For the sake of clarity, we repeat some of the content
presented in earlier sections.
30

Our algorithm stems from a linear programming view on the separated sparsity recovery. It has already
been shown that this LP is totally unimodular [10], which implies that solving the LP provides an integral
solution and hence solves separated sparsity recovery. However, the previous work resorts to a black-box
approach for solving this LP leading to a prohibitive O(d3.5 ) time complexity. We make a step forward, and
study the dual LP, obtaining a method that takes nearly-linear time to find its optimal solution. By the strong
duality, the value of the dual is the value of the primal LP as well. However, unfortunately, it is not hard to
see that a generic relation between primal and dual LP solutions known as "complementary slackness" does
not lead to recovery of the primal optimal solution itself.
To cope with that shortcoming, we further analyze the properties of dual solution. We exhibit very close
connection between its structure and the sparsity of the primal solution, which we present via the notion of
"active constraints". Intuitively, this structure allows us to characterize the cases in which the complementary
slackness in fact provides the primal from a dual optimal solution. Then we use these findings in our algorithm
to slightly perturb the input instance, while not affecting the value of the solution, so that it is possible to
obtain a primal from a dual solution in the general case, and hence solve the separated sparsity recovery.
In Section 5.1 we defined LP P, which corresponds to the problem (3), as follows:
cT u

maximize

d
X

subject to

ui = k

i=1
min{i+∆−1,n}

X

uj ≤ 1

∀i = 1 . . . d

ui ≥ 0

∀i = 1 . . . d

j=i

We also defined the dual LP to P, denoted by D, as
w0 k +

minimize

d
X

wi

i=1

subject to

w0 +

X

wj ≥ ci

∀i = 1 . . . d

wi ≥ 0

∀i = 1 . . . d

j : j≥1 and
j≤i≤j+∆−1

w0 ∈ R
As Theorem 9, which proofs appears in Section F, this dual LP has a combinatorial structure that enables us
to solve it in nearly-linear time.
Before delving into details, we introduce some notation. In what follows, we will be interested in cost
vectors and sparsity parameters other than c and k, respectively. Hence, whenever this is the case, we will
write P(c0 , k 0 ) to refer to the program P for cost vector c0 and sparsity k 0 . Similarly, whenever we consider
some different cost vector c0 and sparsity parameter k 0 , we will denote the corresponding dual LP by D(c0 , k 0 ).
As already mentioned, our running-time results are given with respect to γ, where γ is the maximal number
of bits needed to store any ci . As P and D are invariant under shifting and multiplication of c, for the sake of
clarity of our exposition and without loss of generality we assume c is an integral vector with non-negative
entries.
31

At a high-level, we develop our algorithm in three main steps. In Section 5.1 we considered the dual of
the LP P and show that its combinatorial structure can be exploited to compute an optimal solution to it in
nearly-linear time. By strong duality, this optimal dual solution gives us then the value of the optimal primal
solution. But, unfortunately, it does not give us the optimal primal solution itself.
To alleviate this issue, we develop a divide-and-conquer procedure for extracting that optimal primal
solution. First, in Section J.1, we demonstrate that by analyzing answers of the dual oracle on perturbed
versions of the original problem we can quickly recover a single non-zero entry of the optimal primal solution.
That entry can be used to partition our instance into two smaller subproblems. Then, in Section J.2, we show
how to make these two subproblems fully independent by devising an optimal split of the sparsity constraint
that they share. This optimal split is extracted from the structure of the dual solutions.
With these components in place, we assemble our final algorithm in Section J.3.

J.1

Recovering a segment of an optimal solution

At this point, we developed a way of computing the optimal value OP T of P. However, this is not sufficient
for our purposes as the separated sparsity problem also requires us to provide the corresponding solution, i.e.,
a binary vector u? corresponding to OP T .
One might hope that this primal solution u? can be inferred from the dual solution ŵ that our algorithm
Opt-Value-of-D(c0 , k 0 ) (see Algorithm 2) provides. It is not hard to see, however, that the generic relationship
between the optimal primal and optimal dual solutions that the so-called “complementary slackness” provides
is not sufficient here.
Therefore, we instead design a problem-specific algorithm for finding the desired primal solution vector
u? . As a first step, we focus on the task of recovering a single segment of u? . Specifically, given some target
segment S of at most ∆ consecutive entries, we want to either find a single index i ∈ S such that u?i = 1, for
some optimal primal solution u? ; or to conclude that u?i = 0 for all i ∈ S.
To recover such a segment, we define in an adaptive manner a family of cost vectors c1 , . . . , ct , for some
t ∈ O(log ∆) and invoke our dual LP solver Opt-Value-of-D(c0 , k 0 ) on these cost vectors. As we show, the
solutions to these perturbed instances allow us to infer u?j for all j ∈ S.
To provide more details, let us fix some ∆-length segment S = [js , je ], i.e., je = min{js + ∆ − 1, n}.
We now want to decide whether there is an optimal solution û to P such that ûi = 1 for some i ∈ [js , je ].
Observe that there might be another optimal solution ũ to P. Furthermore, it might be the case that ũj = 0
for every index j ∈ [js , je ], while there is an index i ∈ [js , je ] such that ûi = 1. So while designing the
algorithm, we distinguish two cases. First, we analyse the case in which there exists an optimal solution û to
P so that ûi = 0 for every i ∈ [js , je ]. Next, we consider the complementary case in which for every optimal
solution û, we have that ûi = 1 for some i ∈ [js , je ]. In the latter case, we recover an index j such that there
is a solution û for which ûj = 1 holds.
The first case is captured by the following claim.
0

Lemma 17 Let S ⊆ [d0 ] be a set of indices, let c0 ∈ Nd+ be a coefficient vector, and let k 0 be the sparsity.
0
Define a vector c00 ∈ Nd+ as follows:
(
1
if i ∈ S
.
c00i =
0
1 + ci otherwise
Then, Opt-Value-of-D(c00 , k 0 ) equals Opt-Value-of-D(c0 , k 0 ) + k 0 iff there exists an optimal solution u? to
P(c0 , k 0 ) such that u?i = 0 whenever i ∈ S.
32

Proof. Let us show the two direction of equivalence separately.
( =⇒ ) Assume that Opt-Value-of-D(c00 , k 0 ) equals Opt-Value-of-D(c0 , k 0 ) + k 0 . Let û be an optimal
solution to P(c0 , k 0 ) for the cost vector given by c00 . Now we want to show that ûi = 0 for every i ∈
a contradiction,
the case, i.e. there exists j ∈ S such that ûj = 1. But
S. Towards
P
Passume it00 is not
0
then, i : ûi =1 (ci + 1) > i : ûi =1 ci as ci + 1 > c00i for i ∈ S, and hence Opt-Value-of-D(c00 , k 0 ) <
Opt-Value-of-D(c0 , k 0 ) + k 0 , contradicting our assumption.
( ⇐= ) First, observe that Opt-Value-of-D(c00 , k 0 ) ≤ Opt-Value-of-D(c0 , k 0 ) + k 0 as c00 ≤ c0 + 1.
On the other hand, if there exists an optimal solution u? to P(c0 , k 0 ) such that u?i = 0 whenever i ∈
S, then u? achieves value Opt-Value-of-D(c0 , k 0 ) + k 0 in P(c0 , k 0 ) for the costs given by c00 . In other
words, restricted to the set of indices outside of S, c0 is a shifted by 1 variant of c00 . Therefore, we also
have Opt-Value-of-D(c00 , k 0 ) ≥ Opt-Value-of-D(c0 , k 0 ) + k 0 . Now this implies Opt-Value-of-D(c00 , k 0 ) =
Opt-Value-of-D(c0 , k 0 ) + k 0 , as desired. 2
Next we provide algorithm ∆-recovery that we use in the proof of the next theorem.
Algorithm 6 ∆-recovery:
0

Input: c0 ∈ Nd+ , sparsity k 0 , index js
Output: index r as described in Theorem 10
1. je ← min{js + ∆ − 1, d0 }
2. s ← js(; e = je ; OP T ← Opt-Value-of-D(c0 , k 0 )

1
if i ∈ [js , je ]
0
1 + ci otherwise
4. If Opt-Value-of-D(c00 , k 0 ) = OP T + k 0 then return -1
5. While s < e
 s+e 
(a) mid ←
( 2
1
if i ∈ [js , je ] \ [s, mid]
(b) c00i ←
0
1 + ci otherwise
(c) If Opt-Value-of-D(c00 , k 0 ) = OP T + k 0 then e ← mid
(d) Else s ← mid + 1
3. c00i ←

6. Return s

0

Theorem 10 There exists an algorithm that given c0 ∈ Nd+ , sparsity k 0 , and an index js ∈ [1, d0 ], outputs an
integer r having the following properties:
• If for every optimal solution û to P(c0 , k 0 ) there is an index i such that ûi = 1 and i ∈ [js , min{js +
∆ − 1, d0 }], then r is set to be an index in [js , min{js + ∆ − 1, d0 }] such that there is an optimal
solution u? for which we have u?r = 1.
• Otherwise, r is set to -1.
Furthermore, if the algorithm runs in time O ((log ∆)d0 (log (c0max + 1) + log k 0 )).
33

Proof. On line 4, ∆-recovery first checks whether we can simply ignore all the entries indexed by
{i, . . . , min{i+∆−1, d0 }}. And if yes, it returns -1. (This step is formalized in the statement of Theorem 10.)
However, entries from those interval can only be disregarded if there exists an optimal solution u? to P(c0 , k 0 )
such that for every j ∈ {i, . . . , min{i + ∆ − 1, d0 }} we have u?j = 0. So, if there is no such u? , i.e. the
entries corresponding to that interval have to be considered, the lines following line 4 of ∆-recovery serve
to pinpoint an entry j ∈ {i, . . . , min{i + ∆ − 1, d0 }} so that there exists an optimal solution û to P(c0 , k 0 )
for which holds ûj = 1. The way it is done is by applying a binary search over the interval that can not be
ignored, tracked via variables s and e of the algorithm.
Correctness. Algorithm ∆-recovery outputs −1 correctly by Lemma 17 for S = {js , . . . , je }.
Next, we show that when the binary search loop starting at line 5 ends we have s = e, i.e. [s, e]
corresponds to a single index. As long as ∆ ≥ 1 and js ≤ d0 , at line 1 and line 2 values s and e are initialized
so that s ≤ e, so initially [s, e] is indeed a non-empty interval. Furthermore, as we have s < e in each
iteration, it holds mid < e. But it also holds mid ≥ s (in fact mid equals s only when s + 1 equals e). So,
updating e to mid at line 5c, or s to mid + 1 at line 5d, [s, e] remains a non-empty interval in any iteration.
Notice that the update rules also guarantee that we either increase s or decrease e at each iteration, and
therefore [s, e] shrinks its size by 1 at least at each iteration. Putting this together, and taking into account the
loop-termination condition at line 5, we conclude that after the loop ends it holds s = e.
We say that interval [s, e] has property Rs,e if there is an optimal solution û to P(c0 , k 0 ) such that ûi = 1
for some i ∈ [s, e]. Now, if we show that after line 4 at every step of the algorithm interval [s, e] has property
Rs,e , then the proof will follow immediately. That is exactly how we proceed. Namely, we show that if [s, e]
has property Rs,e , then after updating s or e at line 5c or line 5d obtaining s0 and e0 , respectively, then interval
[s0 , e0 ] will have property Rs0 ,e0 .
Observe that at the beginning of the very first iteration of the loop property Rs,e , equivalent to Rjs ,je ,
holds as by Lemma 17 every optimal solution û is such that ûi = 1 for some i ∈ [js , je ]. Next, let [s, e] be
updated to [s0 , e0 ]. We want to show that property Rs0 ,e0 holds as well.
First, assume that line 5c gets executed, i.e. s0 = s and e0 = mid. Then, by Lemma 17, we have that
there exists an optimal solution û to P(c0 , k 0 ) such that ûi = 0 for every i ∈ [js , je ] \ [s, mid]. But we also
have that there is i ∈ [js , je ] such that ûi = 1. So, putting it together, we conclude that property Rs0 ,e0 holds.
Next, assume that the if condition at line 5c does not hold. So, line 5d is executed, i.e. s0 = mid + 1 and
0
e = e. Then again by Lemma 17, and as a consequence of both line 4 and line 5c, we have that for every
optimal solution û we have ûi = 1 for some ∈ [js , je ] \ [s, mid]. Also, as je − js + 1 ≥ ∆, we also have
ûj = 0 for every j ∈ [s, mid]. But then, as by our assumption there is an optimal solution ũ and an index i
such that ũi = 1 and i ∈ [s, e], we have that i ∈ [mid + 1, e]. Hence property Rmid+1,e , which is equivalent
to Rs0 ,e0 , holds.
Running time. If ∆-recovery outputs -1 at line 4, it invokes Opt-Value-of-D for the cost vector c0 and c00 .
Note that c0max ≤ c0max + 1. Then, by Theorem 9 that case takes O (d0 (log (c0max + 1) + log k 0 )) time.
If the method does not output -1, then it enters the while loop. The loop applies a standard binary search
over js and je which are set so that it holds je − js + 1 ≤ ∆. So, the loop iterates for O(log ∆) times. Every
iteration invokes Opt-Value-of-D for vector c00 as defined at line 5b, resulting in the total running time of
O((log ∆)d0 (log (c0max + 1) + log k 0 )). 2

34

J.2

Distributing sparsity

The algorithm we presented in the previous section enables us to quickly recover a target length-∆ segment
of some optimal primal solution u? . We now would like to build on this procedure to develop a divide-andconquer approach to recovering the primal solution u? in full.
Our intention is to use the procedure from Section J.1 to split our input problem into two (smaller
and approximately equally sized) subproblems by recovering a “middle” segment of the solution u? and
then to proceed recursively on each of these subproblems. The difficulty here, however, is that these two
resulting subproblems are not really independent, even though they correspond to separate cost vectors. The
subproblems still share the sparsity constraint. That is, to make these subproblems truly independent, we
must also specify the split of our sparsity “budget” k 0 between them.
In this section, we analyze properties of the dual LP and exhibit very close ties between its structure, as
captured by active constraints, and the optimal sparsity distribution for our two subproblems. Specifically, we
establish the following theorem.
0

Theorem 11 Let c0 ∈ Nd+ , let k 00 be a sparsity parameter, and let indices s and e be such that e − s + 1 ≥ ∆
0
or e equals d0 . Define a vector c00 ∈ Nd+ as follows:
(
−∞ if i ∈ [s, e]
00
ci =
c0i
otherwise
Assume that there is a ∆-separated choice of coordinates of c00 such that k 00 coordinates are chosen, and
none of them has an index in [s, e], i.e., the instance is feasible when restricted to coordinates outside of [s, e].
0 with the property that there exists an optimal ∆-separated
Then, there is an algorithm that outputs kL0 and kR
00
00
solution for cost vector c and sparsity k such that it chooses kL0 coordinates of c00 with indices less than s,
0 coordinates of c00 with indices greater than e.
and kR
The algorithm runs in time O (d0 (log c0max + log k 00 )).
Observe that Theorem 11 essentially gives an algorithm that distributes sparsity in an optimal way and
thus enables us to implement the desired divide-and-conquer approach.
Let us illustrate how to use the algorithm ∆-recovery (Algorithm 6) together with Active-Constraints to
recover an optimal solution for the example c = (3, 2, 1, 4, 1, 1, 2, 1), k = 4 and ∆ = 2. Let u? be an optimal
solution to P. As described above, we proceed by splitting c into two subvectors and solve the separated
sparsity problem on each of them independently. The first step is to find the exact splitting point in c. To
achieve this, we invoke ∆-recovery(c, 4, k). As c4 = 4 is part of any optimal solution, the method returns
index 4 and hence we have u?4 = 1. So at this moment we have learned one index of the optimal solution.
In fact, it also implies that u?3 = u?5 = 0. Hence it remains to learn the remaining k − 1 indices that define
u? . However, we have to distribute the remaining sparsity k − 1 over the two subproblems. To that end, we
define a new vector c0 as follows
(
−∞ if i ∈ {3, 4, 5}
c0 =
ci
otherwise.
The definition of c0 enforces that no optimal solution to the separated sparsity problem for sparsity k − 1 will
choose c3 , c4 , or c5 (setting c0i to the dummy value −∞ is for illustration purposes only). Next, we invoke
Active-Constraints(c0 , Opt-Value-of-D(c0 , k − 1)) and denote its output by active0 . Observe that, as long as
Opt-Value-of-D(c, k) has finite value, an optimal solution to Opt-Value-of-D(c0 , k − 1) is finite as well. We

35

have that active0 = {1, 6, 8}. Now we count the number of elements in active0 that are “to the left” and “to
the right” of u?4 , denoting these quantities by k̂ and k̃ respectively, i.e.,
k̂ = |{i ∈ active0 | i < 4}|,

and

k̃ = |{i ∈ active0 | i > 4}| .

In a similar way, we split the vector c0 into ĉ and c̃ as follows:
ĉ = (c1 , c2 ), and c̃ = (c6 , c7 , c8 ) .
Finally, we solve two ∆-separated instances independently, one for ĉ and k̂, and the second one for c̃ and k̃.
Proof of Theorem 11.
Let us start by presenting an algorithm the distributes the sparsity, and continue by analyzing it.
Algorithm 7 Distribute-Sparsity:
0

Input: c0 ∈ Nd+ , sparsity k 00 to be distributed, indices s and e as described in Theorem 11
0 for the right side
Output: sparsity kL0 for the left side, sparsity kR
00
0
1. shif t (
← 1 + (k cmax + 1)
1
if i ∈ [s, e]
0
ci + shif t otherwise
1
w ← Opt-Value-of-D(c00 , k 00 , −(k 00 − 1)c0max + shif t, c0max + shif t)
active1 ← Active-Constraints(c00 , w1 )
active2 ← Active-Constraints(c00 , Dual-Greedy(c00 , w01 + 1))
kL1 ← |{i ∈ active1 | i < s}|
kL2 ← |{i ∈ active2 | i < s}|
2 ← |{i ∈ active2 | i > e}|
kR
2}
kL0 ← kL2 + min{kL1 − kL2 , k 00 − kL2 − kR
0 ← k 00 − k 0
kR
L
0 )
Return (kL0 , kR

2. c00i ←
3.
4.
5.
6.
7.
8.
9.
10.
11.

Optimality for D(c00 , k 00 ). We argue that w1 defined on line 3 is an optimal solution to D(c00 , k 00 ) although
w01 is restricted to belong to interval [−(k 00 − 1)c0max + shif t, c0max + shif t] = [c0max + 2, c0max + shif t].
First, observe that for w01 ≥ c0max + 2 no constraint with index in [s, e] is tight, and hence no such constraint
can be active. Furthermore, as no constraints in [s, e] is tight, we have wi1 = 0 for all i ∈ [e, s]. Now,
consider cL = (c01 , . . . , c0s−1 ) and cR = (c0e+1 , . . . , c0d0 ) as defined in the algorithm. Then, by following
00 0
the proof
Lemma
18,and as e −
 0 of

 s + 1 ≥ ∆, we conclude that Dual-Greedy(c , cmax + 2) has at least
n−e
s−1
0
min k , ∆
+ min k , ∆
active constraints which, by the construction and the assumption that
00
the input instance is feasible, is k at least. So, by the monotonicity of the size of active constraints given by
Lemma 19 and the optimality condition provided via Lemma 13 we have that w1 is an optimal solution to
D(c00 , k 00 ).
Correctness of distributed sparsity. Let active1 denote Active-Constraints(c00 , w1 ) and active2 denote
Active-Constraints(c00 , w2 ), where w2 ← Dual-Greedy(c00 , w01 + 1). Next, split w1 and w2 as follows. Let
1 be a zero-indexed s-dimensional and w 1 be a zero-indexed (d0 − e + 1)-dimensional vector defined as
wL
R
1
1
1
1
wL
← (w01 , w11 , . . . , ws−1
), and wR
← (w01 , we+1
, . . . , wd10 ).

36

1 and w 1 , define
Intuitively, letter ’L’ stands for left to s and letter ’R’ stands for right to e. Similarly to wL
R
2 and w 2 as
wL
R
2
2
2
2
wL
← (w02 , w12 , . . . , ws−1
), and wR
← (w02 , we+1
, . . . , wd20 ).

Also, split active1 and active2 with respect to s and e in the obvious way
active1L ← {i ∈ active1 | i < s}, and active1R ← {i ∈ active1 | i > e},
and
active2L ← {i ∈ active2 | i < s}, and active2R ← {i ∈ active2 | i > e}.
As a reminder, there is no i ∈ [s, e] such that i ∈ active1 . Furthermore, as there is no tight constraint in [s, e]
for w1 , we have wi1 = 0 for every i ∈ [s, e]. Observe that the same holds for active2 and w2 . In addition, we
have e − s + 1 ≥ ∆. From this, we can derive the following list of equalities, which essentially allows us to
1;
split the input problem into two independent subproblems. So, we have: Dual-Greedy(cL , w01 ) equals wL
1
2
1
1
1
Dual-Greedy(cL , w0 +1) equals wL ; Dual-Greedy(cR , w0 ) equals wR ; and, Dual-Greedy(cR , w0 +1) equals
2 . But also, we have: Active-Constraints(c , w 1 ) equals active1 ; Active-Constraints(c , w 2 ) equals
wR
L
L
L
L
L
1 ) equals active1 ; and, Active-Constraints(c , w 2 ) equals active2 .
active2L ; Active-Constraints(cR , wR
R
R
R
R
Now, by Lemma 19 we conclude that |active1L | ≥ |active2L | and |active1R | ≥ |active2R |. This in turn
0 in Distribute-Sparsity are derived correctly. (Note that k 0 + k 0 = k 00 .) But then, by
implies that kL0 and kR
L
R
1 is an optimal solution to D(c , k 0 ) and w 1 is an optimal solution to D(c , k 0 ),
Lemma 13 we have that wL
L L
R R
R
0 ) knowing that the optimal solution, i.e.
so we can independently solve Recover(cL , kL0 ) and Recover(cR , kR
0 ← Opt-Value-of-D(c , k ) might differ from w 1 , and
its value, remains the same. We point out that wL
L L
L
0
1 . In fact, w 0 and w 0 might be such that
similarly wR ← Opt-Value-of-D(cR , kR ) might differ from wR
L
R
0 6= w 0 , although we have that w 1 is equal to w 1 .
wL
0
R0
L0
R0
Running time. Every line, except maybe line 3, take O(d0 ) time. On the other hand, line 3 takes time
O (d0 log {c0max + shif t − (−(k 00 − 1)c0max + shif t)}) which is equal to O (d0 (log c0max + log k 00 )). 2

J.3

Separated sparsity in nearly-linear time

We now have all components in place to state our final algorithm Recover(c0 , k 0 ) that produces an optimal
integral solution for P(c0 , k 0 ) (see Algorithm 8).
We prove the following result for our algorithm.
Theorem 12 The algorithm Recover(c0 , k 0 ) solves the model projection problem for separated sparsity in
time O (n(log (cmax + 1) + log k) log n log ∆).
We start by showing some technical lemmas. Recall that by Lemma 11 we provided an upper bound on
any w0? such that w? is an optimum of D. On the other hand, there is no lower bound on w0? . To see that,
consider a very simple example c = (1), ∆ = k = 1. Nevertheless, we can provide a lower bound in the
following form.
Lemma 18 There exists an optimal solution w? to D(c0 , k 0 ) such that w0? ≥ −(k 0 − 1)c0max .
Proof. To prove the lemma, we utilize the optimality condition described by Lemma 13 and the following
claim.
37

Algorithm 8 Recover:
0

Input: c0 ∈ Nd+ , sparsity k 0
Output: a primal solution d0 -dimensional vector u?
1. If k 0 = 0 return 0 j
k
0
2. r ← ∆-recovery(c0 , d 2+1 , k 0 )
3. If r = −1
j then
k
0
1. s ← d 2+1 ,
2. k 00 ← k 0

e ← min{s + ∆ − 1, d0 }

4. Else
1. s ← max{r − ∆ + 1, 1},
2. k 00 ← k 0 − 1

e ← min{r + ∆ − 1, d0 }

0 ) ← Distribute-Sparsity(c0 , k 00 , s, e)
5. (kL0 , kR
6. cL ← (c01 , . . . , c0s−1 ), cR ← (c0e+1 , . . . , c0d0 )
0 )
7. uL ←
Recover(cL , kL0 ), uR ← Recover(cR , kR
uL
if i < s

i


1
if i ∈ [s, e] and r = i
8. u?i =
R

ui−e if i > e



0
otherwise
?
9. Return u

Lemma 19 Let active1 be the output of Active-Constraints(c0 , Dual-Greedy(c0 , v)) and active2 be the
output of Active-Constraints(c0 , Dual-Greedy(c0 , v − t)), for any integer v and a positive integer t. Then, it
holds
|active1 | ≤ |active2 |.
Proof. Let S and S 0 be tight constraints for Dual-Greedy(c0 , v) and Dual-Greedy(c0 , v − t), respectively.
Then, as we have discussed, S ⊆ S 0 . So, all we have to show is that there are at least as many active
constraints formed from S 0 as there are formed from S. We do that by induction on the size of S under the
assumption that S ⊆ S 0 .
Base of induction: |S| = 0. As the number of active constraints is non-negative, and if |S| = 0 there is no
active constraint, the claim follows.
Inductive step: |S| > d0 and S ⊆ S 0 , for d0 ≥ 0. Let imin and i0min be the smallest index of S and S 0 ,
respectively. Constraint imin is active for Dual-Greedy(c0 , v), and i0min is active for Dual-Greedy(c0 , v − t).
Define T = S \ {imin , . . . , imin + ∆ − 1} and T 0 = S 0 \ {i0min , . . . , i0min + ∆ − 1}. Since i0min ≤ imin it
holds T ⊆ T 0 . Now, as for every active constraint j there is no other active one in the neighborhood of ∆ − 1
around j and except that neighborhood the other constraints are not affected by j, we have that
active1 = {imin } ∪ {active constraints for T }, and
active2 = {imin0 } ∪ {active constraints for T 0 }.

38

Now, as T ⊆ T 0 and |T | < |S|, by inductive hypothesis we have
|{active constraints for T }| ≤ |{active constraints for T 0 }|,
and hence the lemma follows. 2
Let ŵ ← Dual-Greedy(c0 , −(k − 1)cmax ) and active ← Active-Constraints(c0 , ŵ). We show that
|active| ≥ k 0 . Furthermore, we show that l∆ + 1 ∈ active, for all l = 0, . . . , k 0 − 1. Once it is shown, the
claim follows by Lemma 19 and Lemma 13. Precisely, Lemma 19 shows that the number of active constraints
for solutions corresponding to Dual-Greedy(c0 , −(k 0 − 1)c0max − t), for t ≥ 0, is at least k 0 . But then, from
Lemma 13 we have that Dual-Greedy(c0 , −(k 0 − 1)c0max ) or Dual-Greedy(c0 , −(k 0 − 1)c0max + t), for some
t > 0, outputs an optimal solution.
So, it only remain to show |active| ≥ k. We prove that by induction, showing the following property.
Let
i = l∆ + 1 be an index for some integer l such that 0 ≤ l ≤ k 0 − 1. If ŵ0 = −(k 0 − 1)c0max ,
Pi−1
0
j=max{1,i−∆+1} ŵj ≤ lcmax and no constraint in {(l − 1)∆ + 2, . . . , l∆} is active, then constraints
l∆ + 1, (l + 1)∆ + 1, . . . , (k 0 − 1)∆ + 1 are active. The induction is applied in a downward fashion on l, i.e.
the base case is l = k 0 − 1, and our goal is to show it holds for l = 0.
P
0
0
Base of induction: l = k 0 − 1. Let i = (k 0 − 1)∆ + 1. As i−1
j=max{1,i−∆+1} ŵj ≤ (k − 1)cmax and
P
0
c0i ≥ 0, we have i−1
j=max{1,i−∆+1} ŵj ≤ ŵ0 + ci , and hence constraint i is tight. Furthermore, as no
constraint in {i − ∆ + 1, . . . , i − 1} is active, constraint i is an active one.
P
0
Inductive step: 0 ≤ l < k 0 − 1. Let i = l∆ + 1. First, if i−1
j=max{1,i−∆+1} ŵj ≤ lcmax we have that
0
constraint i is tight, and as beforePactive as well, and also we have ŵi ≥ −ŵ0 − lcmax .
0
0
Next, we want to show that i+∆
j=i+1 ŵj ≤ (l + 1)cmax (note that i + ∆ ≤ d ), so that we can use the
0
inductive hypothesis for l + 1. First, we show that for every index i it holds
0

i
X

≤ c0max − ŵ0 .

(6)

j=max{1,i0 −∆+1}

P
Recall that ŵ0 is negative. Towards a contradiction, assume that there exists index q such that qj=max{1,q−∆+1} >
c0max − ŵ0 . In case of tie, let q be the smallest such index, which implies ŵq > 0. But then, it
contradicts
the greedy choice of Dual-Greedy algorithm, as by the greedy choice we have ŵq = 0 or
Pq
0
0
j=max{1,q−∆+1} = cq − ŵ0 ≤ cmax − ŵ0 .
Now we combine (6), for i0 = i + ∆ − 1 and ŵi ≥ −ŵ0 − lc0max to obtain
i+∆
X

ŵj ≤ c0max − ŵ0 − ŵi ≤ (l + 1)c0max ,

j=i+1

as desired. In addition, as constraint i is active, we have that no constraints in {i + 1, . . . , i + ∆ − 1} is
active, and hence we can use the inductive hypothesis.
This concludes the proof. 2

39

J.4

Proof of Theorem 12

The algorithm Recover(c0 , k 0 ) utilizes ∆-recovery and Distribute-Sparsity to split the input problem into two
subproblems and then recurses on them. Let us now analyze the correctness and the running time of this
algorithm.
Base case. If the input sparsity is 0, the algorithm outputs a zero-vector at line 1.
∆-middle entries. Next, the algorithm invokes ∆-recovery over the ∆ middle entries. Depending on r, it
sets s and e to correspond to the interval of c0 that should be removed from consideration in the recursive
calls. In the same time, that intervals serves to break the input problem into two independent subproblems. It
also sets k 00 that represents the sparsity distributed outside interval [s, e]. Correctness of this call is guaranteed
by Theorem 10. Value shif t, and in turn vector c00 , is set so that we have a guarantee that no element from
the interval [s, e] is chosen as part of an optimal solution. To see that, consider vector c(3) ← c00 − shif t · 1.
(3)
(3)
Then, ci = c0i for i ∈
/ [s, e] and ci = −(k 0 c0max + 1). In other words, as long as k 00 ≤ k 0 entries can be
chosen outside of the interval [s, e], no optimal solution should choose any entry within the interval [s, e].
Distributing the sparsity. The correctness of line 5 follows by Theorem 11.
Running time. Line 2 of algorithm Recover runs in time O ((log ∆)d0 (log (c0max + 1) + log k 0 )) which is
a subset of O ((log ∆)d0 (log (cmax + 1) + log k)). Vector c0 is split into cL and cR at line 6 in O(d0 ) time.
Afterwards, the method recurses on the two subproblems, and combines their outputs into u? . Obtaining u?
at line 8 also takes O(d0 ). So, it only to remain to discuss the recursion.
From Theorem 11 we have that line 5 take O (d0 (log c0max + log k 00 )). Furthermore, as k 00 ≤ k and
c0max ≤ cmax , we have O (d0 (log c0max + log k 0 )) = O (d0 (log cmax + log k)).
Every recursive step shrinks the corresponding c0 vector by half at least. So, the recursion has depth
of O(log d). At every level of recursion are considered vectors c0 of total length O(d) – this follows from
the recursive call at line 7 and the fact the cL and cR represent disjoint pieces of c0 . So, this implies the
total running time of O (log d(d (log cmax + log k) + (log ∆)d(log (cmax + 1) + log k))), which in a more
compact way can be written as O (d(log (cmax + 1) + log k) log d log ∆).

40

